text:
using PowerModelsInterface
using PowerSystems
import InfrastructureSystems
using PowerSystemCaseBuilder
import PowerSystemCaseBuilder: PSITestSystems
using PowerModels
using Dates
import JuMP
import Ipopt
import Cbc
import SCS
import Memento

# Test Packages
using Test
using Logging

# Code Quality Tests
import Aqua
Aqua.test_unbound_args(PowerModelsInterface)
Aqua.test_undefined_exports(PowerModelsInterface)
Aqua.test_ambiguities(PowerModelsInterface)

ipopt_solver = JuMP.optimizer_with_attributes(
    Ipopt.Optimizer,
    "tol" => 1e-4,
    "print_level" => 0,
    "max_iter" => 100000,
)
ipopt_win_solver = JuMP.optimizer_with_attributes(
    Ipopt.Optimizer,
    "tol" => 1e-4,
    "print_level" => 0,
    "max_iter" => 1000000,
)

cbc_solver = JuMP.optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
# juniper_solver = JuMP.optimizer_with_attributes(Juniper.Optimizer, "nl_solver"=>JuMP.optimizer_with_attributes(Ipopt.Optimizer, "tol"=>1e-4, "print_level"=>0), "log_levels"=>[])
scs_solver = JuMP.optimizer_with_attributes(
    SCS.Optimizer,
    "max_iters" => 1000000,
    "eps" => 1e-4,
    "verbose" => 0,
)

const MODEL_SOLVER_MAP = Dict(
    ACPPowerModel => Sys.iswindows() ? ipopt_win_solver : ipopt_solver,
    SOCWRConicPowerModel => scs_solver,
)

const PM = PowerModels
const PSY = PowerSystems
const PMI = PowerModelsInterface
const PSB = PowerSystemCaseBuilder
const IS = InfrastructureSystems

PM_BASE_DIR = abspath(joinpath(dirname(Base.find_package("PowerModels")), ".."))
PM_DATA_DIR = joinpath(PM_BASE_DIR, "test", "data")

LOG_FILE = "power-models-interface.log"
LOG_LEVELS = Dict(
    "Debug" => Logging.Debug,
    "Info" => Logging.Info,
    "Warn" => Logging.Warn,
    "Error" => Logging.Error,
)

memento_logger = Memento.config!("error")

"""
Copied @includetests from https://github.com/ssfrr/TestSetExtensions.jl.
Ideally, we could import and use TestSetExtensions.  Its functionality was broken by changes
in Julia v0.7.  Refer to https://github.com/ssfrr/TestSetExtensions.jl/pull/7.
"""

"""
Includes the given test files, given as a list without their ".jl" extensions.
If none are given it will scan the directory of the calling file and include all
the julia files.
"""
macro includetests(testarg...)
    if length(testarg) == 0
        tests = []
    elseif length(testarg) == 1
        tests = testarg[1]
    else
        error("@includetests takes zero or one argument")
    end

    quote
        tests = $tests
        rootfile = @__FILE__
        if length(tests) == 0
            tests = readdir(dirname(rootfile))
            tests = filter(
                f ->
                    startswith(f, "test_") && endswith(f, ".jl") && f != basename(rootfile),
                tests,
            )
        else
            tests = map(f -> string(f, ".jl"), tests)
        end
        println()
        for test in tests
            print(splitext(test)[1], ": ")
            include(test)
            println()
        end
    end
end

function get_logging_level_from_env(env_name::String, default)
    level = get(ENV, env_name, default)
    return IS.get_logging_level(level)
end

function run_tests()
    logging_config_filename = get(ENV, "SIIP_LOGGING_CONFIG", nothing)
    if logging_config_filename !== nothing
        config = IS.LoggingConfiguration(logging_config_filename)
    else
        config = IS.LoggingConfiguration(
            filename = LOG_FILE,
            file_level = Logging.Info,
            console_level = Logging.Error,
        )
    end
    console_logger = ConsoleLogger(config.console_stream, config.console_level)

    IS.open_file_logger(config.filename, config.file_level) do file_logger
        levels = (Logging.Info, Logging.Warn, Logging.Error)
        multi_logger =
            IS.MultiLogger([console_logger, file_logger], IS.LogEventTracker(levels))
        global_logger(multi_logger)

        if !isempty(config.group_levels)
            IS.set_group_levels!(multi_logger, config.group_levels)
        end

        # Testing Topological components of the schema
        @time @testset "Begin PowerModelsInterface tests" begin
            @includetests ARGS
        end

        @test length(IS.get_log_events(multi_logger.tracker, Logging.Error)) == 1
        @info IS.report_log_summary(multi_logger)
    end
end

logger = global_logger()

try
    run_tests()
finally
    # Guarantee that the global logger is reset.
    global_logger(logger)
    nothing
end

max_stars_repo_name: NREL-SIIP/PowerModelsInterface

max_issues_repo_name: NREL-SIIP/PowerModelsInterface

max_forks_repo_name: NREL-SIIP/PowerModelsInterface


text:
#=
day7:
- Julia version: 1.7.0
- Author: Paul.Mealor
- Date: 2021-12-07
=#

using AoC, Test

parsePositions(line) = parse.(Int, split(line, ","))
@test parsePositions("5,34,3") == [5,34,3]

totalDistance(candidate, positions) = sum(p -> abs(candidate - p), positions)
@test totalDistance(4,[4,4]) == 0
@test totalDistance(1,[16,1,2,0,4,2,7,1,2,14]) == 41
@test totalDistance(2,[16,1,2,0,4,2,7,1,2,14]) == 37
@test totalDistance(3,[16,1,2,0,4,2,7,1,2,14]) == 39
@test totalDistance(10,[16,1,2,0,4,2,7,1,2,14]) == 71

totalDistance2(candidate, positions) = sum(positions) do p
    distance = abs(candidate - p)
    return distance * (distance + 1) ÷ 2
end
@test totalDistance2(4,[4,4]) == 0
@test totalDistance2(5,[16,1,2,0,4,2,7,1,2,14]) == 168
@test totalDistance2(2,[16,1,2,0,4,2,7,1,2,14]) == 206

function bestPosition(positions)
    bestDistance = typemax(eltype(positions))
    result = -1
    for candidate = min(positions...):max(positions...)
        thisDistance = totalDistance(candidate, positions)
        if thisDistance < bestDistance
            bestDistance = thisDistance
            result = candidate
        end
    end
    return result
end

function bestPosition2(positions)
    bestDistance = typemax(eltype(positions))
    result = -1
    for candidate = min(positions...):max(positions...)
        thisDistance = totalDistance2(candidate, positions)
        if thisDistance < bestDistance
            bestDistance = thisDistance
            result = candidate
        end
    end
    return result
end

@test bestPosition([16,1,2,0,4,2,7,1,2,14]) == 2
@test bestPosition2([16,1,2,0,4,2,7,1,2,14]) == 5

function part1(line)
    positions = parsePositions(line)
    position = bestPosition(positions)
    return totalDistance(position, positions)
end

@test part1(firstExampleLine(7,1)) == 37

function part2(line)
    positions = parsePositions(line)
    position = bestPosition2(positions)
    return totalDistance2(position, positions)
end

@test part2(firstExampleLine(7,1)) == 168

firstLine(7) |> l -> @time part1(l) |> show
firstLine(7) |> l -> @time part1(l) |> show
firstLine(7) |> l -> @time part2(l) |> show
firstLine(7) |> l -> @time part2(l) |> show

ext: jl


text:
################################################################################
# Define wrapped functions
module AddIntegersCxx
using CxxInterface

const libAddIntegersCxx = joinpath(pwd(), "libAddIntegersCxx")

eval(cxxsetup())
eval(cxxnewfile("AddIntegersCxx.cxx", ""))

eval(cxxfunction(FnName(:add_int, "add_int", libAddIntegersCxx), FnResult(Cint, "int", Int, expr -> :(convert(Int, $expr))),
                 [FnArg(:x, Cint, "x", "int", Integer, identity), FnArg(:y, Cint, "y", "int", Integer, identity)], "return x + y;"))
end

################################################################################
# Extract C++ code
using CxxInterface
using .AddIntegersCxx
AddIntegersCxx.cxx_write_code!()

################################################################################
# Compile C++ code
# (This fails if there is no C++ compiler available)
using Libdl: dlext
run(`c++ -fPIC -c AddIntegersCxx.cxx`)
run(`c++ -shared -o libAddIntegersCxx.$dlext AddIntegersCxx.o`)

# Please, DO NOT call a C++ compiler manually in your own Julia
# packages. This works only in very controlled environments such as on
# CI infrastructure. If you do, your package will be fragile, and will
# create lots of headaches for your users in the wild. Instead, use
# [BinaryBuilder](https://binarybuilder.org) and store your build
# recipes on [Yggdrasil](https://github.com/JuliaPackaging/Yggdrasil).

################################################################################
# Call the wrapped function
# (This fails if the C++ compiler is not compatible with the currently
# running Julia executable)

# Only test cases where the Github CI environment supports this (:i686 is not supported)
if Sys.ARCH ≡ :x86_64
    @test AddIntegersCxx.add_int(2, 3) == 5
end


text:
# This file is part of the IntervalArithmetic.jl package; MIT licensed

using IntervalArithmetic
using Base.Test

setprecision(Interval, Float64)
# # setrounding(Interval, :narrow)

@testset "Numeric tests" begin

    a = @interval(0.1, 1.1)
    b = @interval(0.9, 2.0)
    c = @interval(0.25, 4.0)

    ## Basic arithmetic
    @test @interval(0.1) == Interval(9.9999999999999992e-02, 1.0000000000000001e-01)
    @test +a == a
    @test a+b == Interval(9.9999999999999989e-01, 3.1000000000000001e+00)
    @test -a == Interval(-1.1000000000000001e+00, -9.9999999999999992e-02)
    @test a-b == Interval(-1.9000000000000001e+00, 2.0000000000000018e-01)
    @test Interval(1//4,1//2) + Interval(2//3) == Interval(11//12, 7//6)
    @test Interval(1//4,1//2) - Interval(2//3) == Interval(-5//12, -1//6)

    @test 10a == @interval(10a)
    #@test 10Interval(1//10) == one(@interval(1//10))
    @test Interval(-30.0,-15.0) / Interval(-5.0,-3.0) == Interval(3.0, 10.0)
    @test @interval(-30,-15) / @interval(-5,-3) == Interval(3.0, 10.0)
    @test b/a == Interval(8.18181818181818e-01, 2.0000000000000004e+01)
    @test a/c == Interval(2.4999999999999998e-02, 4.4000000000000004e+00)
    @test c/4.0 == Interval(6.25e-02, 1e+00)
    @test c/zero(c) == emptyinterval(c)
    @test Interval(0.0, 1.0)/Interval(0.0,1.0) == Interval(0.0, Inf)
    @test Interval(-1.0, 1.0)/Interval(0.0,1.0) == entireinterval(c)
    @test Interval(-1.0, 1.0)/Interval(-1.0,1.0) == entireinterval(c)
    a = @interval(1.e-20)
    @test a == Interval(1.0e-20, 1.0000000000000001e-20)
    @test diam(a) == eps(1.e-20)
end

@testset "Power tests" begin
    @test @interval(0,3) ^ 2 == Interval(0, 9)
    @test @interval(2,3) ^ 2 == Interval(4, 9)
    @test @interval(-3,0) ^ 2 == Interval(0, 9)
    @test @interval(-3,-2) ^ 2 == Interval(4, 9)
    @test @interval(-3,2) ^ 2 == Interval(0, 9)
    @test @interval(0,3) ^ 3 == Interval(0, 27)
    @test @interval(2,3) ^ 3 == Interval(8, 27)
    @test @interval(-3,0) ^ 3 == @interval(-27., 0.)
    @test @interval(-3,-2) ^ 3 == @interval(-27, -8)
    @test @interval(-3,2) ^ 3 == @interval(-27., 8.)
    @test Interval(0,3) ^ -2 == Interval(1/9, Inf)
    @test Interval(-3,0) ^ -2 == Interval(1/9, Inf)
    @test Interval(-3,2) ^ -2 == Interval(1/9, Inf)
    @test Interval(2,3) ^ -2 == Interval(1/9, 1/4)
    @test Interval(1,2) ^ -3 == Interval(1/8, 1.0)
    @test Interval(0,3) ^ -3 == @interval(1/27, Inf)
    @test Interval(-1,2) ^ -3 == entireinterval()
    @test_throws ArgumentError interval(-1, -2) ^ -3  # wrong way round
    @test Interval(-3,2) ^ (3//1) == Interval(-27, 8)
    @test Interval(0.0) ^ 1.1 == Interval(0, 0)
    @test Interval(0.0) ^ 0.0 == emptyinterval()
    @test Interval(0.0) ^ (1//10) == Interval(0, 0)
    @test Interval(0.0) ^ (-1//10) == emptyinterval()
    @test ∅ ^ 0 == ∅
    @test Interval(2.5)^3 == Interval(15.625, 15.625)
    #@test Interval(5//2)^3.0 == Interval(125//8)

    x = @interval(-3,2)
    @test x^3 == @interval(-27, 8)

    @test @interval(-3,4) ^ 0.5 == @interval(0, 2)
    @test @interval(-3,4) ^ 0.5 == @interval(-3,4)^(1//2)
    @test @interval(-3,2) ^ @interval(2) == Interval(0.0, 4.0)
    @test @interval(-3,4) ^ @interval(0.5) == Interval(0, 2)
    @test @biginterval(-3,4) ^ 0.5 == @biginterval(0, 2)

    @test dist(@interval(1,27)^@interval(1/3), Interval(1., 3.)) < 2*eps(Interval(1,3))
    @test dist(@interval(1,27)^(1/3), Interval(1., 3.)) < 2*eps(Interval(1,3))
    @test Interval(1., 3.) ⊆ @interval(1,27)^(1//3)
    @test @interval(0.1,0.7)^(1//3) == Interval(0.46415888336127786, 0.8879040017426008)
    @test dist(@interval(0.1,0.7)^(1/3),
        Interval(0.46415888336127786, 0.8879040017426008)) < 2*eps(@interval(0.1,0.7)^(1/3))

    setprecision(Interval, 256)
    x = @biginterval(27)
    y = x^(1//3)
    @test (0 < diam(y) < 1e-76)
    y = x^(1/3)
    @test (0 < diam(y) < 1e-76)

end

setprecision(Interval, Float64)

@testset "Exp and log tests" begin
    @test exp(@biginterval(1//2)) ⊆ exp(@interval(1//2))
    @test exp(@interval(1//2)) == Interval(1.648721270700128, 1.6487212707001282)
    @test exp(@biginterval(0.1)) ⊆ exp(@interval(0.1))
    @test exp(@interval(0.1)) == Interval(1.1051709180756475e+00, 1.1051709180756477e+00)
    @test diam(exp(@interval(0.1))) == eps(exp(0.1))

    @test log(@biginterval(1//2)) ⊆ log(@interval(1//2))
    @test log(@interval(1//2)) == Interval(-6.931471805599454e-01, -6.9314718055994529e-01)
    @test log(@biginterval(0.1)) ⊆ log(@interval(0.1))
    @test log(@interval(0.1)) == Interval(-2.3025850929940459e+00, -2.3025850929940455e+00)
    @test diam(log(@interval(0.1))) == eps(log(0.1))

    @test exp2(@biginterval(1//2)) ⊆ exp2(@interval(1//2))
    @test exp2(Interval(1024.0)) == Interval(1.7976931348623157e308, Inf)
    @test exp10(@biginterval(1//2)) ⊆ exp10(@interval(1//2))
    @test exp10(Interval(308.5)) == Interval(1.7976931348623157e308, Inf)

    @test log2(@biginterval(1//2)) ⊆ log2(@interval(1//2))
    @test log2(@interval(0.25, 0.5)) == Interval(-2.0, -1.0)
    @test log10(@biginterval(1//10)) ⊆ log10(@interval(1//10))
    @test log10(@interval(0.01, 0.1)) == @interval(log10(0.01), log10(0.1))
end

@testset "Comparison tests" begin
    d = @interval(0.1, 2)

    @test d < 3
    @test d <= 2
    @test (d < 2) == false
    @test -1 < d
    @test !(d < 0.15)

    # abs
    @test abs(@interval(0.1, 0.2)) == Interval(9.9999999999999992e-02, 2.0000000000000001e-01)
    @test abs(@interval(-1, 2)) == Interval(0, 2)

    # real
    @test real(@interval(-1, 1)) == Interval(-1, 1)
end

@testset "Rational tests" begin

    f = 1 // 3
    g = 1 // 3

    @test @interval(f*g) == Interval(1.1111111111111109e-01, 1.1111111111111115e-01)
    @test big(1.)/9 ∈ @interval(f*g)
    @test @interval(1)/9 ⊆ @interval(f*g)
    @test @interval(1)/9 ≠ @interval(f*g)

    h = 1/3
    i = 1/3

    @test @interval(h*i) == Interval(1.1111111111111109e-01, 1.1111111111111115e-01)
    @test big(1.)/9 ∈ @interval(1/9)

    @test @interval(1/9) == @interval(1//9)
end

@testset "Floor etc. tests" begin
    a = @interval(0.1)
    b = Interval(0.1, 0.1)
    @test dist(a,b) <= eps(a)

    @test floor(@interval(0.1, 1.1)) == Interval(0, 1)
    @test round(@interval(0.1, 1.1), RoundDown) == Interval(0, 1)
    @test ceil(@interval(0.1, 1.1)) == Interval(1, 2)
    @test round(@interval(0.1, 1.1), RoundUp) == Interval(1, 2)
    @test sign(@interval(0.1, 1.1)) == Interval(1.0)
    @test trunc(@interval(0.1, 1.1)) == Interval(0.0, 1.0)
    @test round(@interval(0.1, 1.1), RoundToZero) == Interval(0.0, 1.0)
    @test round(@interval(0.1, 1.1)) == Interval(0.0, 1.0)
    @test round(@interval(0.1, 1.5)) == Interval(0.0, 2.0)
    @test round(@interval(-1.5, 0.1)) == Interval(-2.0, 0.0)
    @test round(@interval(-2.5, 0.1)) == Interval(-2.0, 0.0)
    @test round(@interval(0.1, 1.1), RoundTiesToEven) == Interval(0.0, 1.0)
    @test round(@interval(0.1, 1.5), RoundTiesToEven) == Interval(0.0, 2.0)
    @test round(@interval(-1.5, 0.1), RoundTiesToEven) == Interval(-2.0, 0.0)
    @test round(@interval(-2.5, 0.1), RoundTiesToEven) == Interval(-2.0, 0.0)
    @test round(@interval(0.1, 1.1), RoundTiesToAway) == Interval(0.0, 1.0)
    @test round(@interval(0.1, 1.5), RoundTiesToAway) == Interval(0.0, 2.0)
    @test round(@interval(-1.5, 0.1), RoundTiesToAway) == Interval(-2.0, 0.0)
    @test round(@interval(-2.5, 0.1), RoundTiesToAway) == Interval(-3.0, 0.0)

    # # :wide tests
    # # setrounding(Interval, :wide)
    # setprecision(Interval, Float64)
    #
    # a = @interval(-3.0, 2.0)
    # @test a == Interval(-3.0, 2.0)
    # @test a^3 == Interval(-27, 8)
    # @test Interval(-3,2)^3 == Interval(-27, 8)
    #
    # @test Interval(-27.0, 8.0)^(1//3) == Interval(0, 2.0000000000000004)
    #
    # # setrounding(Interval, :narrow)
end

@testset "Fast power" begin

    @testset "Fast integer powers" begin
        x = 1..2
        @test pow(x, 2) == pow(-x, 2) == Interval(1, 4)
        @test pow(-x, 3) == Interval(-8.0, -1.0)

        @test pow(-1..2, 2) == 0..4
        @test pow(-1..2, 3) == -1..8
        @test pow(-1..2, 4) == 0..16

        @test pow(@biginterval(-1, 2), 2) == 0..4
        @test pow(@biginterval(-1, 2), 3) == -1..8
        @test pow(@biginterval(1, 2), 2) == 1..4

        x = @interval(pi)
        @test x^100 ⊆ pow(x, 100)
        @test x^50 ⊆ pow(x, 50)
        @test isinterior(x^50, pow(x, 50))

        x = Interval(2)
        @test pow(x, 2000) == Interval(realmax(), Inf)
    end

    @testset "Fast real powers" begin
        x = 1..2
        @test pow(x, 0.5) == Interval(1.0, 1.4142135623730951)
        @test pow(x, 0.5) == x^0.5

        @test pow(x, 17.3) != x^17.3

        y = 2..3
        @test pow(y, -0.5) == Interval(0.5773502691896257, 0.7071067811865476)

        y = -2..3
        @test pow(y, 2.1) == Interval(0.0, 10.045108566305146)
        @test pow(y, 2.1) != y^2.1
        @test y^2.1 ⊆ pow(y, 2.1)
    end

    @testset "Fast interval powers" begin
        x = 1..2
        @test x^Interval(-1.5, 2.5) == Interval(0.35355339059327373, 5.656854249492381)

        y = -2..3
        @test pow(y, 2.1) == Interval(0.0, 10.045108566305146)
        @test pow(y, Interval(-2, 3)) == Interval(0, Inf)

        @test pow(y, @interval(2.1)) == Interval(0.0, 10.045108566305146)
    end

    # comment out for now:
    #=
    @testset "Float32 intervals" begin

        a = Interval{Float32}(1e38)
        b = Interval{Float32}(1e2)
        @test a * b == Interval{Float32}(realmax(Float32), Inf)
    end
    =#
end


text:
###############################################################################
#
# Create canvas
#
###############################################################################
"""
    create_canvas(
                id::Union{Int,String};
                ncol::Int = 1,
                nrow::Int = 1,
                axids::Array{Int,1} = Int[],
                figsize::Tuple{Number,Number} = (0.5+ncol*3, 0.5+nrow*3),
                dpi::Number = 100
    )

Create a canvas, given
- `id` ID of the figure
- `ncol` Number of columns in the figure
- `nrow` Number of rows in the figure
- `axids` Given indicies of the subplots in the figure
- `figsize` Given canvas size
- `dpi` Given pixels per inch
"""
function create_canvas(
            id::Union{Int,String};
            ncol::Int = 1,
            nrow::Int = 1,
            axids::Array{Int,1} = Int[],
            figsize::Tuple{Number,Number} = (0.5+ncol*3, 0.5+nrow*3),
            dpi::Number = 100
)
    # create a clean figure
    fig = figure(id, figsize=figsize, dpi=dpi);
    clf();

    # if axids is empty
    if axids==[]
        axids = collect(1:(nrow*ncol));
    end

    # create axes
    axs = [];
    for indx in axids
        ax = subplot(nrow, ncol, indx);
        push!(axs, ax);
    end

    return fig,axs
end

###############################################################################
#
# Save canvas
#
###############################################################################
"""
    save_canvas!(fig::Figure, path::String, saving::Bool = false)

Save a canvas, given
- `fig` Canvas to save
- `path` Path to save the figure
- `saving` Optional. If true, save the figure
"""
function save_canvas!(fig::Figure, path::String, saving::Bool = false)
    fig.set_tight_layout(true);
    if saving
        fig.savefig(path, bbox_inches="tight");
    end

    return nothing
end


text:
function location2idx(::LatticeGrid, f, dir, slice_location::Number)
    m, idx = findmin(map(x->abs(x-slice_location), getproperty(getdomain(f), dir)))
    @debug "Closest index to slice location $slice_location is $idx at $m"

    return idx
end

function location2idx(::ParticleGrid, f, dir, slice_location::Number; ϵ)
    grid = getdomain(f)
    grid_axis = getproperty(grid, dir)
    idxs = filter(i->grid_axis[i] ∈ slice_location ± ϵ, axes(grid_axis, 1))
    @debug "Found $(length(idxs)) indices close to slice location $slice_location"

    return idxs
end

function Base.selectdim(f::T, dir::Symbol, slice_location; kwargs...) where T <: AbstractPICDataStructure
    @debug "Generic slice fallback for AbstractPICDataStructure"
    idx = location2idx(domain_discretization(T), f, dir, slice_location; kwargs...)
    selectdim(scalarness(T), domain_discretization(T), f, dir, idx)
end

function Base.selectdim(::ScalarQuantity, ::LatticeGrid, f, dir, idx)
    @debug "Scalar selectdim along $dir at index $idx"
    dim = dir_to_idx(dir)
    # TODO: Figure out if we can get rid of dir_to_idx
    data = selectdim(unwrapdata(f), dim, idx)
    grid = dropdims(getdomain(f), dims=dir)

    newstruct(f, data, grid)
end

function Base.selectdim(::ScalarQuantity, ::ParticleGrid, f, dir, idx)
    @debug "Scalar selectdim along $dir at index $idx"
    dim = dir_to_idx(dir)
    # TODO: Figure out if we can get rid of dir_to_idx
    data = selectdim(unwrapdata(f), dim, idx)
    grid_lower_dim = dropdims(getdomain(f), dims=dir)

    grid = selectdim(grid_lower_dim, dir, idx)

    newstruct(f, data, grid)
end

function Base.selectdim(::VectorQuantity, ::LatticeGrid, f, dir, idx::Int)
    dim = dir_to_idx(dir)
    full_data = unwrapdata(f)
    grid = getdomain(f)

    sliced_data = selectdim(full_data, dim, idx:idx)
    f_sliced = newstruct(f, sliced_data, grid)
    f_reduced = dropdims(f_sliced, dims=dir)

    dropgriddims(f_reduced, dims=dir)
end

function Base.selectdim(::VectorQuantity, ::ParticleGrid, v, dir, idx)
    @debug "Vector selectdim along $dir at index $idx"
    full_data = unwrapdata(v)
    dim = dir_to_idx(dir)

    grid_lower_dim = dropdims(getdomain(v), dims=dir)
    grid = selectdim(grid_lower_dim, dir, idx)
    @debug grid
    # TODO: Figure out if we can get rid of dir_to_idx
    sliced_data = selectdim(full_data, dim, idx)
    @debug sliced_data
    data = removedims(sliced_data; dims=dir)

    newstruct(v, data, grid)
end

function removedims(data; dims)
    dims = dims isa Symbol ? (dims,) : dims
    selected_dims = (setdiff(propertynames(data),dims)...,)
    @debug "Keeping grid dims: $selected_dims"

    StructArray(component.((data,), selected_dims), names=selected_dims)
end

Base.dropdims(f::T; dims) where T <: AbstractPICDataStructure = dropdims(scalarness(T), f; dims)

function Base.dropdims(::VectorQuantity, f; dims)
    data_rm = removedims(unwrapdata(f); dims)
    dim = dir_to_idx.(dims)
    data = dropdims(data_rm; dims=dim)
    grid = getdomain(f)

    newstruct(f, data, grid)
end

function dropgriddims(f; dims)
    grid = getdomain(f)
    data = unwrapdata(f)
    sliced_grid = dropdims(grid; dims)

    newstruct(f, data, sliced_grid)
end


text:
### struct for defining the Muller-Brown potential
mutable struct CustomTrainableModel <: AbstractTrainablePotential
  dim::Int64
  nparams::Int64
  theta::Vector{Float64}
  pfunc::Union{Function,nothing}
  ffunc::Union{Function,nothing}
  f::Vector{Float64}
  function CustomTrainableModel(dim::Int64,pfunc::Union{Function,nothing},ffunc::Union{Function,nothing})
    new(dim,pfunc,ffunc2,zeros(Float64,dim))
  end
end
CustomTrainableModel(dim::Int64;pfunc::Union{Function,nothing}=nothing,ffunc::Union{Function,nothing}=nothing) = CustomTrainableModel(dim,pfunc,ffunc)

function getdimensionality(model::CustomTrainableModel)
  return model.dim
end

function jacshape(model::CustomTrainableModel)
  return (model.nparams,model.dim)
end

function gradshape(model::CustomTrainableModel)
  return model.nparams
end

"""
Potential functions
"""
function potential(system::AbstractSystem, model::CustomTrainableModel)
  if model.pfunc != nothing
    return model.pfunc(system)
  end
end

"""
Force functions with return value
"""
function force(system::AbstractSystem, model::CustomTrainableModel)
  if model.ffunc != nothing
    return model.ffunc(system)
  else
    return zeros(Float64,model.dim)
  end
end

"""
Force functions with in-place changing
"""
function force!(system::AbstractSystem, model::CustomTrainableModel)
  if model.ffunc != nothing
    model.ffunc(system,model.f)
  else
    modelf .= 0.0
  end
end

"""
Jacobian wrt parameters for neural network model
"""
function jacobian!(jac::Array{Float64},sys::AbstractSystem,model::CustomTrainableModel)
  # TODO
  grads::Grads = jacobian(() -> nnm.ffunc(sys), nnm.pars)
  # store gradients into jacobian array
  bcount::Int64 = 1
  for i in 1:length(nnm.pars)
    # TODO replace with reduce function
    ptup::Tuple = nnm.pinfo[i]
    ecount::Int64 = 1
    for j in 1:length(ptup)
      ecount *= ptup[j]
    end
    ecount += (bcount - 1)
    @. @views jac[bcount:ecount,:] .= grads[nnm.pars[i]]'
    bcount = ecount + 1
  end
end

"""
Update parameters for the model
"""
function updateparams!(lr::Float64,dtheta::Vector{Float64},model::CustomTrainableModel)
    # update coefficients for forces
    model.theta .+= (lr .* dtheta)
end

"""
Set parameters for the model
"""
function setparams!(theta::Vector{Float64},model::CustomTrainableModel)
    # update coefficients for forces
    model.theta .= theta
end


text:
module UtilityFunctionsLH

using ArgCheck, DocStringExtensions
export CRRA, utility, marginal_utility, cons_age1, cons_growth, cons_path,
    euler_deviation, lifetime_utility, mu_wealth,
    pv_consumption
export discounted_sum

# -------------  Types

abstract type UtilityFunction end
abstract type UtilityOneArg <: UtilityFunction end

"""
	$(SIGNATURES)

Utility = c^(1-sigma) / (1-sigma) - 1.0.

Log utility is special case when `sigma == 1.0`.
"""
mutable struct CRRA <: UtilityOneArg
    sigma  ::  Float64
    dbg  ::  Bool
 end

include("crra.jl")

## ---------  Helpers

"""
	$(SIGNATURES)

Discounted sum: `xV[1] + ... + pBeta ^ (T-1) * xV[T]`
"""
function discounted_sum(pBeta :: T1, xV) where T1 <: AbstractFloat
    T = length(xV);
    s = zero(eltype(xV));
    for t = T : -1 : 1
        s = pBeta * s + xV[t];
    end
    return s
end

# Vector with constant growth factor
function const_growth_vector(c1 :: T1,
    growthFactor :: T1, T :: Integer) where T1 <: AbstractFloat

    cV = Vector{T1}(undef, T);
    ct = c1;
    for j in eachindex(cV)
        cV[j] = ct;
        ct *= growthFactor;
    end
    return cV
end

# Sum of a constant growth vector: 1 + growthFactor + ... + growthFactor ^ (T-1)
function const_growth_sum(growthFactor :: T1, T :: Integer) where T1
    if isapprox(growthFactor, one(T1))
        gSum = T1(T);
    else
        gSum = (one(T1) - growthFactor ^ T) / (one(T1) - growthFactor);
    end
    return gSum
end

end # module


text:
module MLJBase

# ===================================================================
# IMPORTS

import Base: ==, precision, getindex, setindex!
import Base.+, Base.*

# Scitype
import MLJScientificTypes.ScientificTypes: TRAIT_FUNCTION_GIVEN_NAME
using MLJScientificTypes
using MLJModelInterface

# Interface
import MLJModelInterface: fit, update, update_data, transform,
    inverse_transform, fitted_params, predict,
    predict_mode, predict_mean, predict_median, predict_joint,
    evaluate, clean!, is_same_except,
    save, restore, is_same_except, istransparent,
    params

# Macros
using Parameters

# Containers & data manipulation
using Tables
import PrettyTables
using DelimitedFiles
using OrderedCollections
using CategoricalArrays
import InvertedIndices: Not
import JLSO
import Dates

# Distributed computing
using Distributed
using ComputationalResources
using ComputationalResources: CPUProcesses
using ProgressMeter
import .Threads

# Operations & extensions
import LossFunctions
import LossFunctions: DistanceLoss, MarginLoss, SupervisedLoss
import StatsBase
import StatsBase: fit!, mode, countmap
import Missings: levels
import Distributions
import Distributions: pdf, logpdf, sampler
const Dist = Distributions

# from Standard Library:
using Statistics, LinearAlgebra, Random, InteractiveUtils

# external loss functions:
import LossFunctions: DWDMarginLoss, ExpLoss, L1HingeLoss, L2HingeLoss,
    L2MarginLoss, LogitMarginLoss, ModifiedHuberLoss, PerceptronLoss,
    SigmoidLoss, SmoothedL1HingeLoss, ZeroOneLoss, HuberLoss, L1EpsilonInsLoss,
    L2EpsilonInsLoss, LPDistLoss, LogitDistLoss, PeriodicLoss, QuantileLoss

# ===================================================================
## EXPORTS

# -------------------------------------------------------------------
# re-exports from MLJModelInterface, MLJScientificTypes
# NOTE: MLJBase does **not** re-export UnivariateFinite to avoid
# ambiguities between the raw constructor (MLJBase.UnivariateFinite)
# and the general method (MLJModelInterface.UnivariateFinite)

# MLJ model hierarchy
export MLJType, Model, Supervised, Unsupervised,
    Probabilistic, JointProbabilistic, Deterministic, Interval, Static,
    ProbabilisticComposite, JointProbabilisticComposite, DeterministicComposite,
    IntervalComposite, UnsupervisedComposite, StaticComposite,
    ProbabilisticSurrogate, JointProbabilisticSurrogate, DeterministicSurrogate,
    IntervalSurrogate, UnsupervisedSurrogate, StaticSurrogate,
    Surrogate, Composite

export UnivariateFinite

# MLJType equality
export is_same_except

# model constructor + metadata
export @mlj_model, metadata_pkg, metadata_model

# model api
export fit, update, update_data, transform, inverse_transform,
    fitted_params, predict,
    predict_mode, predict_mean, predict_median, predict_joint,
    evaluate, clean!

# model/measure matching:
export Checker, matching

# model traits
export input_scitype, output_scitype, target_scitype,
    is_pure_julia, package_name, package_license,
    load_path, package_uuid, package_url,
    is_wrapper, supports_weights, supports_online,
    docstring, name, is_supervised,
    prediction_type, implemented_methods, hyperparameters,
    hyperparameter_types, hyperparameter_ranges

# data operations
export matrix, int, classes, decoder, table,
    nrows, selectrows, selectcols, select

# re-exports from MLJScientificTypes
export Unknown, Known, Finite, Infinite,
    OrderedFactor, Multiclass, Count, Continuous, Textual,
    Binary, ColorImage, GrayImage, Image, Table
export scitype, scitype_union, elscitype, nonmissing, trait
export coerce, coerce!, autotype, schema, info

# -------------------------------------------------------------------
# exports from this module, MLJBase

export DeterministicComposite,
    ProbabilisticComposite,
    UnsupervisedComposite, @load

# computational_resources.jl:
export default_resource

# one_dimensional_ranges.jl:
export ParamRange, NumericRange, NominalRange, iterator, scale

# parameter_inspection.jl:
export params # note this is *not* an extension of StatsBase.params

# data.jl:
export partition, unpack, complement, restrict, corestrict

# utilities.jl:
export flat_values, recursive_setproperty!,
       recursive_getproperty, pretty, unwind

# show.jl
export HANDLE_GIVEN_ID, @more, @constant, @bind, color_on, color_off

# univariate_finite/
export average, UnivariateFiniteArray, UnivariateFiniteVector

# info_dict.jl:
export info_dict

# datasets.jl:
export load_boston, load_ames, load_iris,
       load_reduced_ames, load_crabs, load_smarket,
       @load_boston, @load_ames, @load_iris,
       @load_reduced_ames, @load_crabs, @load_smarket

# sources.jl:
export source, Source, CallableReturning

# machines.jl:
export machine, Machine, fit!, report, fit_only!

# datasets_synthetics.jl
export make_blobs, make_moons, make_circles, make_regression

# composition:
export machines, sources, @from_network, @pipeline,
    glb, @tuple, node, @node, sources, origins, return!,
    nrows_at_source, machine,
    rebind!, nodes, freeze!, thaw!, models, Node, AbstractNode,
    DeterministicSurrogate, ProbabilisticSurrogate, UnsupervisedSurrogate,
    DeterministicComposite, ProbabilisticComposite, UnsupervisedComposite

# aliases to the above,  kept for backwards compatibility:
export  DeterministicNetwork, ProbabilisticNetwork, UnsupervisedNetwork

# resampling.jl:
export ResamplingStrategy, Holdout, CV, StratifiedCV,
       evaluate!, Resampler, PerformanceEvaluation

# openml.jl:
export OpenML

# -------------------------------------------------------------------
# exports from MLJBase specific to Measure (these may go in their
# specific MLJMeasureInterface package in some future)

# measures/registry.jl:
export measures, metadata_measure

# measure/measures.jl:
export orientation, reports_each_observation,
    is_feature_dependent, aggregation,
    aggregate, default_measure, value,
    supports_class_weights, prediction_type, human_name

# measures/continuous.jl:
export mav, mae, mape, rms, rmsl, rmslp1, rmsp, l1, l2, log_cosh,
    MAV, MAE, MeanAbsoluteError, mean_absolute_error, mean_absolute_value,
    LPLoss, RootMeanSquaredProportionalError, RMSP,
    RMS, rmse, RootMeanSquaredError, root_mean_squared_error,
    RootMeanSquaredLogError, RMSL, root_mean_squared_log_error, rmsl, rmsle,
    RootMeanSquaredLogProportionalError, rmsl1, RMSLP,
    MAPE, MeanAbsoluteProportionalError, log_cosh_loss, LogCosh, LogCoshLoss

# measures/confusion_matrix.jl:
export confusion_matrix, confmat, ConfusionMatrix

# measures/finite.jl:
export cross_entropy, BrierScore, brier_score,
    BrierLoss, brier_loss,
    LogLoss, log_loss,
    misclassification_rate, mcr, accuracy,
    balanced_accuracy, bacc, bac, BalancedAccuracy,
    matthews_correlation, mcc, MCC, AUC, AreaUnderCurve,
    MisclassificationRate, Accuracy, MCR, BACC, BAC,
    MatthewsCorrelation

# measures/finite.jl -- Multiclass{2} (order independent):
export auc, area_under_curve, roc_curve, roc

# measures/finite.jl -- OrderedFactor{2} (order dependent):
export TruePositive, TrueNegative, FalsePositive, FalseNegative,
    TruePositiveRate, TrueNegativeRate, FalsePositiveRate,
    FalseNegativeRate, FalseDiscoveryRate, Precision, NPV, FScore,
    NegativePredictiveValue,
    # standard synonyms
    TPR, TNR, FPR, FNR, FDR, PPV,
    Recall, Specificity, BACC,
    # instances and their synonyms
    truepositive, truenegative, falsepositive, falsenegative,
    true_positive, true_negative, false_positive, false_negative,
    truepositive_rate, truenegative_rate, falsepositive_rate,
    true_positive_rate, true_negative_rate, false_positive_rate,
    falsenegative_rate, negativepredictive_value,
    false_negative_rate, negative_predictive_value,
    positivepredictive_value, positive_predictive_value,
    tpr, tnr, fpr, fnr,
    falsediscovery_rate, false_discovery_rate, fdr, npv, ppv,
    recall, sensitivity, hit_rate, miss_rate,
    specificity, selectivity, f1score, fallout

# measures/finite.jl -- Finite{N} - multiclass generalizations of
# above OrderedFactor{2} measures (but order independent):
export MulticlassTruePositive, MulticlassTrueNegative, MulticlassFalsePositive,
      MulticlassFalseNegative, MulticlassTruePositiveRate,
      MulticlassTrueNegativeRate, MulticlassFalsePositiveRate,
      MulticlassFalseNegativeRate, MulticlassFalseDiscoveryRate,
      MulticlassPrecision, MulticlassNegativePredictiveValue, MulticlassFScore,
      # standard synonyms
      MTPR, MTNR, MFPR, MFNR, MFDR, MPPV,
      MulticlassRecall, MulticlassSpecificity,
      # instances and their synonyms
      multiclass_truepositive, multiclass_truenegative,
      multiclass_falsepositive,
      multiclass_falsenegative, multiclass_true_positive,
      multiclass_true_negative, multiclass_false_positive,
      multiclass_false_negative, multiclass_truepositive_rate,
      multiclass_truenegative_rate, multiclass_falsepositive_rate,
      multiclass_true_positive_rate, multiclass_true_negative_rate,
      multiclass_false_positive_rate, multiclass_falsenegative_rate,
      multiclass_negativepredictive_value, multiclass_false_negative_rate,
      multiclass_negative_predictive_value, multiclass_positivepredictive_value,
      multiclass_positive_predictive_value, multiclass_tpr, multiclass_tnr,
      multiclass_fpr, multiclass_fnr, multiclass_falsediscovery_rate,
      multiclass_false_discovery_rate, multiclass_fdr, multiclass_npv,
      multiclass_ppv, multiclass_recall, multiclass_sensitivity,
      multiclass_hit_rate, multiclass_miss_rate, multiclass_specificity,
      multiclass_selectivity, macro_f1score, micro_f1score,
      multiclass_f1score, multiclass_fallout, multiclass_precision,
      # averaging modes
      no_avg, macro_avg, micro_avg

# measures/loss_functions_interface.jl
export dwd_margin_loss, exp_loss, l1_hinge_loss, l2_hinge_loss, l2_margin_loss,
    logit_margin_loss, modified_huber_loss, perceptron_loss, sigmoid_loss,
    smoothed_l1_hinge_loss, zero_one_loss, huber_loss, l1_epsilon_ins_loss,
    l2_epsilon_ins_loss, lp_dist_loss, logit_dist_loss, periodic_loss,
    quantile_loss

# ------------------------------------------------------------------------
# re-export from LossFunctions:
export DWDMarginLoss, ExpLoss, L1HingeLoss, L2HingeLoss, L2MarginLoss,
    LogitMarginLoss, ModifiedHuberLoss, PerceptronLoss, SigmoidLoss,
    SmoothedL1HingeLoss, ZeroOneLoss, HuberLoss, L1EpsilonInsLoss,
    L2EpsilonInsLoss, LPDistLoss, LogitDistLoss, PeriodicLoss,
    QuantileLoss

# -------------------------------------------------------------------
# re-export from Random, StatsBase, Statistics, Distributions,
# OrderedCollections, CategoricalArrays, InvertedIndices:
export pdf, sampler, mode, median, mean, shuffle!, categorical, shuffle,
       levels, levels!, std, Not, support, logpdf, LittleDict

# ===================================================================
## CONSTANTS

# the directory containing this file: (.../src/)
const MODULE_DIR = dirname(@__FILE__)

# horizontal space for field names in `MLJType` object display:
const COLUMN_WIDTH = 24
# how deep to display fields of `MLJType` objects:
const DEFAULT_SHOW_DEPTH = 0
const DEFAULT_AS_CONSTRUCTED_SHOW_DEPTH = 2
const INDENT = 4

const Arr = AbstractArray
const Vec = AbstractVector

const MMI = MLJModelInterface
const FI  = MLJModelInterface.FullInterface
# ===================================================================
# Computational Resource
# default_resource allows to switch the mode of parallelization

default_resource()    = DEFAULT_RESOURCE[]
default_resource(res) = (DEFAULT_RESOURCE[] = res)

# stub for @load (extended by MLJModels)
macro load end

# ===================================================================
# Includes

include("init.jl")
include("utilities.jl")
include("show.jl")
include("info_dict.jl")
include("interface/data_utils.jl")
include("interface/model_api.jl")

include("univariate_finite/types.jl")
include("univariate_finite/methods.jl")
include("univariate_finite/arrays.jl")

include("sources.jl")
include("machines.jl")

include("composition/abstract_types.jl")
include("composition/learning_networks/nodes.jl")
include("composition/learning_networks/inspection.jl")
include("composition/learning_networks/machines.jl")
@static if VERSION ≥ v"1.3.0-"
    include("composition/learning_networks/arrows.jl")
end

include("composition/models/methods.jl")
include("composition/models/from_network.jl")
include("composition/models/inspection.jl")
include("composition/models/pipelines.jl")
include("composition/models/_wrapped_function.jl")

include("operations.jl")
include("resampling.jl")

include("hyperparam/one_dimensional_ranges.jl")
include("hyperparam/one_dimensional_range_methods.jl")

include("data/data.jl")
include("data/datasets.jl")
include("data/datasets_synthetic.jl")

include("matching.jl")
include("measures/measures.jl")
include("measures/measure_search.jl")

include("openml.jl")

end # module


text:
# # Daisyworld
# ![](daisyworld.gif)
#
# Study this example to learn about
# - Simple agent properties with complex model interactions
# - Rolling your own plots
# - Collecting data with the low-level data collection API
# - Simultaneously plotting and collecting data
# - Analyzing the behavior of a model
#
# ## Overview of Daisyworld
#
# This model explores the [Gaia hypothesis](https://en.wikipedia.org/wiki/Gaia_hypothesis),
# which considers the Earth as a single, self-regulating system including both living and
# non-living parts.
#
# Daisyworld is filled with black and white daisies.
# Their albedo's differ, with black daisies absorbing light and heat,
# warming the area around them; white daisies doing the opposite.
# Daisies can only reproduce within a certain temperature range, meaning too much
# (or too little) heat coming from the sun and/or surrounds will ultimately halt daisy
# propagation.
#
# When the climate is too cold it is necessary for the black daisies to propagate in order
# to raise the temperature, and vice versa -- when the climate is too warm, it is
# necessary for more white daisies to be produced in order to cool the temperature.
# The interplay of the living and non living aspects of this world manages to find an
# equilibrium over a wide range of parameter settings, although with enough external
# forcing, the daisies will not be able to self regulate the temperature of the planet
# and eventually go extinct.

# ## Defining the agent type

# The agent here is not so complex. We see it has three values (other than the required
# `id` and `pos` for an agent that lives on a [`GridSpace`](@ref). Each daisy has an `age`,
# confined later by a maximum age set by the user, a `breed` (either `:black` or `:white`)
# and an associated `albedo` value, again set by the user.
using Agents, AgentsPlots, Plots
using Statistics: mean
using Random # hide
gr() # hide

mutable struct Daisy <: AbstractAgent
    id::Int
    pos::Tuple{Int,Int}
    breed::Symbol
    age::Int
    albedo::Float64 # 0-1 fraction
end

# ## World heating

# The surface temperature of the world is heated by its sun, but daisies growing upon it
# absorb or reflect the starlight -- altering the local temperature.

function suface_temperature!(pos::Tuple{Int,Int}, model::ABM{Daisy})
    ids = ids_in_position(pos, model)
    absorbed_luminosity = if isempty(ids)
        ## Set luminosity via surface albedo
        (1 - model.surface_albedo) * model.solar_luminosity
    else
        ## Set luminosity via daisy albedo
        (1 - model[ids[1]].albedo) * model.solar_luminosity
    end
    ## We expect local heating to be 80C for an absorbed luminosity of 1,
    ## approximately 30 for 0.5 and approximately -273 for 0.01.
    local_heating = if absorbed_luminosity > 0
        72 * log(absorbed_luminosity) + 80
    else
        80
    end
    ## Surface temperature is the average of the current temperature and local heating.
    model.temperature[pos] = (model.temperature[pos] + local_heating) / 2
end

function diffuse_temperature!(pos::Int, model::ABM{Daisy}; ratio = 0.5)
    neighbors = nearby_positions(pos, model)
    model.temperature[pos] =
        (1 - ratio) * model.temperature[pos] +
        ## Each neighbor is giving up 1/8 of the diffused
        ## amount to each of *its* neighbors
        sum(model.temperature[neighbors]) * 0.125 * ratio
end
nothing # hide

# ## Initialising Daisyworld

# Here, we construct a function to initialize a Daisyworld. We need to know how many
# daisies of each type to seed the planet with and what their albedo's are. The albedo
# of the planet, as well as how intense the world's star tends to be. Alternatively
# we can provide a `scenario` flag, which alters the stars luminosity in different
# ways.

function daisyworld(;
    griddims = (30, 30),
    max_age = 25,
    init_white = 20, # % cover of the world surface of white breed
    init_black = 20, # % cover of the world surface of black breed
    albedo_white = 0.75,
    albedo_black = 0.25,
    albedo_surface = 0.4,
    solar_luminosity = 0.8,
    scenario = :default,
)
    @assert scenario ∈ [
        :default, # User provided solar_luminosity
        :ramp, # Increase & decrease luminosity over an 850 year period
        :high, # White daisies will prefer this climate
        :low, # Black daisies will prefer this climate
        :ours, # The Sun's equivalent, achieving an equilibrium of daisies
    ]

    space = GridSpace(griddims)
    luminosity = if scenario == :ramp
        0.8
    elseif scenario == :high
        1.4
    elseif scenario == :low
        0.6
    elseif scenario == :ours
        1.0
    else
        solar_luminosity
    end

    properties = Dict(
        :max_age => max_age,
        :temperature => zeros(prod(griddims)),
        :surface_albedo => albedo_surface,
        :solar_luminosity => luminosity,
        :scenario => scenario,
        :tick => 0,
    )
    model = ABM(Daisy, space; properties = properties)
    for _ in 1:(init_white * nv(space) / 100)
        add_agent_single!(model, :white, rand(0:max_age), albedo_white)
    end
    for _ in 1:(init_black * nv(space) / 100)
        add_agent_single!(model, :black, rand(0:max_age), albedo_black)
    end
    for p in positions(model)
        suface_temperature!(p, model)
    end
    return model
end
nothing # hide

# ## Model dynamics

# The final piece of the puzzle is the life-cycle of each daisy. This method defines an
# optimal temperature for growth. If the temperature gets too hot or too cold, daisies
# will not wish to propagate and may even die out. So long as the temperature is favorable,
# daisies compete for land and attempt to spawn a new plant of their `breed` in locations
# close to them.

function propagate!(pos::Tuple{Int,Int}, model::ABM{Daisy})
    agents = collect(agents_in_position(pos, model))
    if !isempty(agents)
        agent = agents[1]
        temperature = model.temperature[pos]
        ## Set optimum growth rate to 22.5C, with bounds of [5, 40]C
        seed_threshold = (0.1457 * temperature - 0.0032 * temperature^2) - 0.6443
        if rand() < seed_threshold
            ## Collect all adjacent position that are empty
            empty_neighbors = Vector{Int}(undef, 0)
            neighbors = nearby_positions(pos, model)
            for n in neighbors
                if isempty(ids_in_position(n, model))
                    push!(empty_neighbors, n)
                end
            end
            if !isempty(empty_neighbors)
                ## Seed a new daisy in one of those position
                seeding_place = rand(empty_neighbors)
                add_agent!(seeding_place, model, agent.breed, 0, agent.albedo)
            end
        end
    end
end
nothing # hide

# Now, we need to write the model and agent step functions for Agents.jl to advance
# Daisyworld's dynamics. Since we have constructed a number of helper functions,
# these methods are quite straightforward.

function solar_activity!(model::ABM{Daisy})
    if model.scenario == :ramp
        if model.tick > 200 && model.tick <= 400
            model.solar_luminosity += 0.005
        end
        if model.tick > 500 && model.tick <= 750
            model.solar_luminosity -= 0.0025
        end
    end
end

function model_step!(model::ABM{Daisy})
    for p in positions(model)
        suface_temperature!(p, model)
        diffuse_temperature!(p, model)
        propagate!(p, model)
    end
    model.tick += 1
    solar_activity!(model)
end

function agent_step!(agent::Daisy, model::ABM{Daisy})
    agent.age += 1
    agent.age >= model.max_age && kill_agent!(agent, model)
end
nothing # hide

# ## Look at the pretty flowers!
# Lets run the model for a bit and see what our world looks like when the solar
# activity is similar to that of our own:

cd(@__DIR__) #src
Random.seed!(165) # hide
model = daisyworld(scenario = :ours)
step!(model, agent_step!, model_step!, 100)
daisycolor(a) = a.breed
plotabm(model; ac = daisycolor, as = 5)

# We can see that this world achieves quasi-equilibrium, where one `breed` does not
# totally dominate the other.

sum(map(a -> [a.breed == :white, a.breed == :black], allagents(model)))

# ---

# Now we'll take a look at some of the complex dynamics this world can manifest.
# Some of these methods are, for the moment, not implemented in
# [AgentsPlots](https://github.com/JuliaDynamics/AgentsPlots.jl), although this does
# give us an opportunity to test out some of the new data collection features in
# Agents.jl v3.0. *Think you have a nice recipe for a plot that would help others?*
# [Send us a pull request](https://github.com/JuliaDynamics/AgentsPlots.jl/pulls)
# or [open an issue](https://github.com/JuliaDynamics/AgentsPlots.jl/issues).

# First, our fluctuating solar luminosity scenario.

model = daisyworld(scenario = :ramp)

# Then, let us initialize some dataframes for our model and agents. We are interested
# in the global surface temperature, the current solar luminosity and populations of
# each daisy breed. Notice that we made sure that `sum` has been given a default value
# since this model is using `kill_agent!` (see [`run!`](@ref) for more details).

global_temperature(model) = mean(model.temperature)
mdata = [global_temperature, :solar_luminosity]
model_df = init_model_dataframe(model, mdata)

white(agent) = agent.breed == :white
black(agent) = agent.breed == :black
total(v) = length(v) == 0 ? 0.0 : sum(v)
adata = [(white, total), (black, total)]
agent_df = init_agent_dataframe(model, adata)
nothing # hide

# Now we can evolve our model and observe what happens

anim = @animate for t in 1:900
    step!(model, agent_step!, model_step!, 1)
    collect_model_data!(model_df, model, mdata, t)
    collect_agent_data!(agent_df, model, adata, t)
    heatmap(
        1:model.space.dimensions[1],
        1:model.space.dimensions[2],
        transpose(reshape(model.temperature, model.space.dimensions));
        clims = (-50, 110),
        colorbar_title = "Temperature",
    )
    scatter!(
        [a.pos for a in allagents(model)];
        marker = (:circle, 5),
        markercolor = [a.breed for a in allagents(model)],
        label = :none,
        showaxis = false,
    )
end
gif(anim, "daisyworld.gif", fps = 10)

# Very interesting! But why is this all happening? Luckily we have collected some useful
# data, so now if we plot our different properties over the same time period, we can see
# how each of the values effect Daisyworld as a whole.

p1 = plot(model_df[!, :solar_luminosity], legend = false, ylabel = "Solar Luminosity")
p2 = plot(model_df[!, :global_temperature], legend = false, ylabel = "Global Temperature")
p3 = plot(
    [agent_df[!, aggname(white, total)], agent_df[!, aggname(black, total)]],
    legend = false,
    ylabel = "Population",
)
plot(p1, p2, p3, layout = (3, 1), size = (500, 800))

# We observe an initial period of low solar luminosity which favors a large population of
# black daisies. The population however is kept in check by competition from white daisies
# and a semi-stable global temperature regime is reached, fluctuating between ~32 and 41
# degrees.
#
# An increase in solar luminosity forces a population inversion, then a struggle for
# survival for the black daisies -- which ultimately leads to their extinction. At
# extremely high solar output the white daisies dominate the landscape, leading to a
# uniform surface temperature.
#
# Finally, as the sun fades back to normal levels, both the temperature and white daisy
# population struggle to find equilibrium. The counterbalancing force of the black daisies
# being absent, Daisyworld is plunged into a chaotic regime -- indicating the strong role
# biodiversity has to play in stabilizing climate.


text:
using PyPlot, DelimitedFiles, Random, Statistics, Test, Revise, transD_GP

Random.seed!(23)

zfixed = [-1e5]
ρfixed = [1e12]
nmax = 100

zstart = 0.
extendfrac, dz = 1.06, 2.
zall, znall, zboundaries = transD_GP.setupz(zstart, extendfrac, dz=dz, n=40, showplot=true)
z, ρ, nfixed = transD_GP.makezρ(zboundaries; zfixed=zfixed, ρfixed=ρfixed)
@info z

## geometry parameters for tempest
zTx = -120
zRx = -80
x_rx = -115.0
y_rx = 0.
rx_roll = 3.
rx_pitch = 6.
rx_yaw = 10.
tx_roll = 5.
tx_pitch = 4.
tx_yaw = -7.

#current ramp and gate times
ramp =  [  # -0.0400066666667    0.5
			-0.0400000000000    0.0
			-0.0399933333333   -0.5
			-0.0200066666667   -0.5
			-0.0200000000000    0.0
			-0.0199933333333    0.5
			-0.0000066666667    0.5
			 0.0000000000000    0.0
			 0.0000066666667   -0.5]

times = sqrt.([
0.0000066667*0.0000200000
0.0000333333*0.0000466667
0.0000600000*0.0000733333
0.0000866667*0.0001266667
0.0001400000*0.0002066667
0.0002200000*0.0003400000
0.0003533333*0.0005533333
0.0005666667*0.0008733333
0.0008866667*0.0013533333
0.0013666667*0.0021000000
0.0021133333*0.0032733333
0.0032866667*0.0051133333
0.0051266667*0.0079933333
0.0080066667*0.0123933333
0.0124066667*0.0199933333
])

## fill in detail in ohm-m
ρ[1:end] .= ρ[1] # no contrasts
## create total field operator (required for nuisance inversion)
tempest = transD_GP.TEMPEST1DInversion.Bfield(
    zTx = zTx, zRx = zRx, x_rx = x_rx, y_rx = y_rx,
    rx_roll = rx_roll, rx_pitch = rx_pitch, rx_yaw = rx_yaw,
    tx_roll = tx_roll, tx_pitch = tx_pitch, tx_yaw = tx_yaw,
	ramp = ramp, times = times,
	z=z,
	ρ=ρ,
	addprimary = true #this ensures that the geometry update actually changes everything that needs to be
)
# plot before testing
transD_GP.TEMPEST1DInversion.plotmodelfield!(tempest,z,ρ)

μ = transD_GP.TEMPEST1DInversion.μ₀
Bx, By, Bz = tempest.Hx[1]*μ*1e15, tempest.Hy[1]*μ*1e15, tempest.Hz[1]*μ*1e15
@testset "Primary field YPR order checks" begin
	@test isapprox(Bx, 30.2751, atol=1e-4)
	@test isapprox(By, -3.7899, atol=1e-4)
	@test isapprox(Bz, -13.8092, atol=1e-4)
end


text:
import YAML
# import Pkg; Pkg.add("OrderedCollections")

# Pkg.add("YAML")

using OrderedCollections
using CSV
using DataFrames
using Dates
using Query
using Statistics

# data = YAML.load_file("../mintsDefinitionsV2.yaml")

yamlFile = YAML.load_file("../mintsDefinitionsV2.yaml"; dicttype=OrderedDict{String,Any})      # loading yaml file and arrange in order
dirctory = yamlFile["dataFolder"]                            # taking the directory of data stored
#sensorNames = yamlFile["liveStack"]                         # all the sensor names
sensorNames = ["AS7262","BME280","MGS001",]

# only 11th node data is saved
nodeIDs_array = yamlFile["nodeIDs"]                         # node's IDs array
node11 = nodeIDs_array[11]                                  # Taking 11th node
nodeId = node11["nodeID"]                                   # taking 11th nodeID

path = dirctory*"/"*nodeId

############# anaylzing and averaging data of AS7262 sensor
function ss(dd, mm, yy, sensr)
    CSV_fileName_arry = readdir(path*"/"*yy*"/"*mm*"/"*dd*"/")       # taking all the file names from the folder
    num_files = length(CSV_fileName_arry)                            # number of files

    CSV_file = "MINTS_"*nodeId*"_"*sensr*"_"*yy*"_"*mm*"_"*dd*".csv"
    if (CSV_file in CSV_fileName_arry)
        df = CSV.read(path*"/"*yy*"/"*mm*"/"*dd*"/"*CSV_file, DataFrame)
        df.dateTime =  SubString.(string.(df.dateTime), 1, 19)
        df.dateTime = DateTime.(df.dateTime,"yyyy-mm-dd HH:MM:SS")

        df.dateTime = map((x) -> round(x, Dates.Second(30)), df.dateTime)
        gdf = groupby(df, :dateTime)                                                                                # making groups by same dateTime
        cgdf = combine(gdf, valuecols(gdf) .=> mean)

        colFullName = names(cgdf)
        for x in 2:length(colFullName)
            num_let = length(colFullName[x])                                        # taking number of characters in one column name
            colName = SubString.(string.(colFullName[x]), 1, num_let-5)             # removing last five characters of each column name
            rename!(cgdf,colFullName[x] => sensr*"_"*colName)                       # adding sensor name to the column name
        end

        return cgdf
    end

end

n = 0

for sensr in sensorNames
    k = 0

    years = readdir(path)

    for y in years
        months = readdir(path*"/"*y)

        for m in months
            days = readdir(path*"/"*y*"/"*m)

            for d in days

                filess = readdir(path*"/"*y*"/"*m*"/"*d)
                CSV_file = "MINTS_"*nodeId*"_"*sensr*"_"*y*"_"*m*"_"*d*".csv"

                if (CSV_file in filess)
                    println("Date: "*y*"_"*m*"_"*d)
                    averaged_df = ss(d, m, y, sensr)
                    df2 = averaged_df
                    k = k + 1

                    if (k > 1)
                        global new_df = outerjoin(new_df, df2, on = intersect(names(new_df), names(df2)))
                    else
                        new_df = averaged_df
                    end
                else
                    println("No data of "*sensr)
                    continue
                end
                #=

                    averaged_df = ss(d, m, y, sensr)
                    df2 = averaged_df

                    k = k + 1

                    if (k > 1)
                        global new_df = outerjoin(new_df, df2, on = intersect(names(new_df), names(df2)))
                    else
                        new_df = averaged_df
                    end

                    =#

            end

        end

    end

    global n = n + 1

    if (n > 1)
        global final_df = outerjoin(final_df, new_df, on = :dateTime, makeunique=true)
    else
        final_df = new_df
    end

end

final_df2 = sort!(final_df)
CSV.write("/home/prabu/dfcomfff.csv", final_df2)

#=

############ read each CSV file and take average
for i in 1:num_files-1
    df = CSV.read("/home/prabu/Research/mintsData/001e06305a6c/2019/07/20/"*CSV_fileName[i], DataFrame)         # Reading ith CSV data file
    df.dateTime =  SubString.(string.(df.dateTime), 1, 19)                                                      # Removing milisecond part of dateTime
    df.dateTime = DateTime.(df.dateTime,"yyyy-mm-dd HH:MM:SS")                                                  # converting to dateTime

    df.dateTime = map((x) -> round(x, Dates.Second(30)), df.dateTime)                                           # rounding dateTime column for 30s
    gdf = groupby(df, :dateTime)                                                                                # making groups by same dateTime
    cgdf = combine(gdf, valuecols(gdf) .=> mean)                                                                # taking average of same group and combine
    push!(dfArray,cgdf)                                                                                # put the DataFrame in the array

    sName =  SubString.(string.(CSV_fileName[i]), 20, 25)                                                       # taking sensor name from CSV file
    push!(sensorName,sName)                                                                                # put sensor name in the array
end

=#

#=
arry = []
for (root, dirs, files) in walkdir(path)
    arr = joinpath.(root, files) # files is a Vector{String}, can be empty
    push!(arry, arr)
end
println(arry[4])
=#

#=
CSV_fileName1 = readdir(path1)       # taking all the file names from the folder
num_files1 = length(CSV_fileName1)                             # number of files
df1 = CSV.read(path1*CSV_fileName1[1], DataFrame)

CSV_fileName2 = readdir(path2)       # taking all the file names from the folder
num_files2 = length(CSV_fileName2)                             # number of files
df2 = CSV.read(path2*CSV_fileName2[1], DataFrame)
=#

#new_df = outerjoin(df1, df2, on = intersect(names(df1), names(df2)))

#CSV.write("/home/prabu/dfcom3.csv", new_df)
#println("Done")


text:
using DataStreamsTest

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

using DataFrames
using CSV

tuple.(CSV.read("test.csv").columns...)


text:
"""
    Interpolate(; limit=nothing)

Performs linear interpolation between the nearest values in an vector.
The current implementation is univariate, so each variable in a table or matrix will
be handled independently.

!!! Missing values at the head or tail of the array cannot be interpolated if there
are no existing values on both sides. As a result, this method does not guarantee
that all missing values will be imputed.

# Keyword Arguments
* `limit::Union{UInt, Nothing}`: Optionally limit the gap sizes that can be interpolated.

# Example
```jldoctest
julia> using Impute: Interpolate, impute

julia> M = [1.0 2.0 missing missing missing 6.0; 1.1 missing missing 4.4 5.5 6.6]
2×6 Matrix{Union{Missing, Float64}}:
 1.0  2.0       missing   missing   missing  6.0
 1.1   missing  missing  4.4       5.5       6.6

julia> impute(M, Interpolate(); dims=:rows)
2×6 Matrix{Union{Missing, Float64}}:
 1.0  2.0  3.0  4.0  5.0  6.0
 1.1  2.2  3.3  4.4  5.5  6.6

julia> impute(M, Interpolate(; limit=2); dims=:rows)
2×6 Matrix{Union{Missing, Float64}}:
 1.0  2.0   missing   missing   missing  6.0
 1.1  2.2  3.3       4.4       5.5       6.6
```
"""
struct Interpolate <: Imputor
    limit::Union{UInt, Nothing}
end

Interpolate(; limit=nothing) = Interpolate(limit)

function _impute!(data::AbstractVector{<:Union{T, Missing}}, imp::Interpolate) where T
    @assert !all(ismissing, data)
    i = findfirst(!ismissing, data) + 1

    while i < lastindex(data)
        if ismissing(data[i])
            prev_idx = i - 1
            next_idx = findnext(!ismissing, data, i + 1)

            if next_idx !== nothing
                gap_sz = (next_idx - prev_idx) - 1

                if imp.limit === nothing || gap_sz <= imp.limit
                    diff = data[next_idx] - data[prev_idx]
                    incr = diff / T(gap_sz + 1)
                    val = data[prev_idx] + incr

                    # Iteratively fill in the values
                    for j in i:(next_idx - 1)
                        data[j] = val
                        val += incr
                    end
                end

                i = next_idx
            else
                break
            end
        end
        i += 1
    end

    return data
end


text:
mutable struct Moon
    x::Int
    y::Int
    z::Int
    i::Int
    j::Int
    k::Int
end

function gravity(moons::Array{Moon,1})
    attract(moons[1], moons[2])
    attract(moons[1], moons[3])
    attract(moons[1], moons[4])
    attract(moons[2], moons[3])
    attract(moons[2], moons[4])
    attract(moons[3], moons[4])
end

function attract(left::Moon, right::Moon)
    if left.x < right.x
        left.i += 1
        right.i += -1
    elseif left.x > right.x
        left.i += -1
        right.i += 1
    else
    end
    if left.y < right.y
        left.j += 1
        right.j += -1
    elseif left.y > right.y
        left.j += -1
        right.j += 1
    else
    end
    if left.z < right.z
        left.k += 1
        right.k += -1
    elseif left.z > right.z
        left.k += -1
        right.k += 1
    else
    end
end

function velocity(moons::Array{Moon,1})
    for moon in moons
        moon.x += moon.i
        moon.y += moon.j
        moon.z += moon.k
    end
end

function potential(moon::Moon)
    return abs(moon.x) + abs(moon.y) + abs(moon.z)
end

function kinetic(moon::Moon)
    return abs(moon.i) + abs(moon.j) + abs(moon.k)
end

function energy(moon::Moon)
    p = potential(moon)
    k = kinetic(moon)
    return p * k
end

function simulate(moons)
    for i in 1:1000
        gravity(moons)
        velocity(moons)
    end
    return sum([energy(moon) for moon in moons])
end

function cycle(moons)
    xs = [(copy(m.x), copy(m.i)) for m in moons]
    x = nothing
    ys = [(copy(m.y), copy(m.j)) for m in moons]
    y = nothing
    zs = [(copy(m.z), copy(m.k)) for m in moons]
    z = nothing
    i = 0
    while true
        i += 1
        gravity(moons)
        velocity(moons)
        if x === nothing
            xt = [(copy(m.x), copy(m.i)) for m in moons]
            if xs == xt
                x = i
            end
        end
        if y === nothing
            yt = [(copy(m.y), copy(m.j)) for m in moons]
            if ys == yt
                y = i
            end
        end
        if z === nothing
            zt = [(copy(m.z), copy(m.k)) for m in moons]
            if zs == zt
                z = i
            end
        end
        x !== nothing && y !== nothing && z !== nothing ? break : continue
    end
    return lcm(lcm(x, y), z)
end

moons = [
    Moon(-3, 10, -1, 0, 0, 0),
    Moon(-12, -10, -5, 0, 0, 0),
    Moon(-9, 0, 10, 0, 0, 0),
    Moon(7, -5, -3, 0, 0, 0)
]

p1 = simulate(moons)
@assert p1 == 10944

moons = [
    Moon(-3, 10, -1, 0, 0, 0),
    Moon(-12, -10, -5, 0, 0, 0),
    Moon(-9, 0, 10, 0, 0, 0),
    Moon(7, -5, -3, 0, 0, 0)
]

p2 = cycle(moons)
@assert p2 == 484244804958744

print("-----------------------------------------------------------------------\n")
print("the n-body problem -- part one\n    total energy : $p1\n")
print("the n-body problem -- part two\n    t-stop       : $p2\n")
print("-----------------------------------------------------------------------\n")


text:
include(joinpath(@__DIR__,"Turtle_struct.jl"))

# based on: https://github.com/cormullion/Lindenmayer.jl
mutable struct LSystem
    rules::Dict{String, String}
    state::Array{Int64, 1}
    initial_state::Array{Int64, 1}
    function LSystem(rules, state_as_string)
        newlsystem = new(rules, string_to_array(state_as_string), string_to_array(state_as_string))
        return newlsystem
    end
end

function string_to_array(str::String)
    temp = Array{Int64, 1}()
    for c in str
        push!(temp, Int(c))
    end
    return temp
end

function array_to_string(arr::Array)
    temp = ""
    for c in arr
        temp = string(temp, Char(abs(c)))
    end
    return temp
end

"""
    evaluate(ls::LSystem, iterations=1)
Apply the rules in the LSystem to the initial state repeatedly. The ls.state array holds
the result.
This must be inefficient, creating a new copy of the state each time......? :(
"""
function evaluate(ls::LSystem, iterations=1; debug=false)
    for i in 1:iterations
        debug && println("iteration $i")
        the_state = Array{Int64, 1}()
        for j in 1:length(ls.state) # each character in state
            s = string(Char(ls.state[j]))
            if haskey(ls.rules, s)
                #  replace it using the rule
                value = ls.rules[s]
                varr = string_to_array(value)
                if ! isempty(value)
                    push!(the_state, varr...)
                end
            else # keep it in
                push!(the_state, ls.state[j])
            end
        end
        ls.state = the_state
        debug == true ? println(array_to_string(ls.state)) : print("")
    end
end

"""
    render(ls::LSystem)
Once the LSystem has been evaluated, the LSystem.state can be drawn.
"""
function render(ls::LSystem, t::Turtles, stepdistance, rotangle, do_till_escaped; debug=false)
    counter = 1
    escaped = false
    for a in ls.state
        command = string(Char(a))
        if command =="F"
            escaped = draw_line(t, 0, stepdistance)
        elseif command =="G"
            escaped = draw_line(t, 0, stepdistance)
        elseif command =="B"
            cw(t, π)
            escaped = draw_line(t, 0, stepdistance)
            cw(t, π)
        elseif command =="V"
            cw(t, π)
            escaped = draw_line(t, 0, stepdistance)
            cw(t, π)
        elseif command =="f"
            escaped = draw_line(t, 0, stepdistance/2)
        elseif command =="b"
            cw(t, π)
            escaped = draw_line(t, 0, stepdistance/2)
        elseif command =="U"
            penup(t)
        elseif command =="D"
            pendown(t)
        elseif command =="+"
            cw(t, rotangle)
        elseif command =="-"
            ccw(t, rotangle)
        elseif command =="r"
            rotangle = deg2rad.([10, 15, 30, 45, 60])[rand(1:end)]
        # elseif command =="T"
        #     randomhue()
        # elseif command =="t"
        #     HueShift(t, 5) # shift hue round the Hue scale (0-360)
        # elseif command =="c"
        #     Randomize_saturation(t) # shift saturation
        # elseif command =="O"
        #     Pen_opacity_random(t)
        elseif command =="l"
            stepdistance = stepdistance + 1 # larger
        elseif command =="s"
            stepdistance = stepdistance - 1 # smaller
        elseif command =="5"
            t.pen.size = 5
        elseif command =="4"
            t.pen.size = 4
        elseif command =="3"
            t.pen.size = 3
        elseif command =="2"
            t.pen.size = 2
        elseif command =="1"
            t.pen.size = 1
        elseif command =="n"
            t.pen.size = 0.5
        # elseif command =="o"
            # Circle(t, stepdistance/4)
        # elseif command =="q"
            # Rectangle(t, stepdistance/4, stepdistance/4)
        elseif command =="["
            push!(t) # push //TODO - need to implimant
        elseif command =="]"
            pop!(t)   # pop //TODO - need to implimant
        end
        counter += 1
        if do_till_escaped && escaped
            break
        end
    end
    return (counter,escaped)
end


text:
# used in shortest path calculations

struct DefaultDistance<:AbstractMatrix{Int}
    nv::Int
    DefaultDistance(nv::Int=typemax(Int)) = new(nv)
end

DefaultDistance(nv::Integer) = DefaultDistance(Int(nv))

show(io::IO, x::DefaultDistance) = print(io, "$(x.nv) × $(x.nv) default distance matrix (value = 1)")
show(io::IO, z::MIME"text/plain", x::DefaultDistance) = show(io, x)

getindex(::DefaultDistance, s::Integer, d::Integer) = 1
getindex(::DefaultDistance, s::UnitRange, d::UnitRange) = DefaultDistance(length(s))
size(d::DefaultDistance) = (d.nv, d.nv)
transpose(d::DefaultDistance) = d
ctranspose(d::DefaultDistance) = d

"""
    eccentricity(g[, v][, distmx])

Return the eccentricity[ies] of a vertex / vertex list `v` or the
entire graph. An optional matrix of edge distances may be supplied; if missing,
edge distances default to `1`.

The eccentricity of a vertex is the maximum shortest-path distance between it
and all other vertices in the graph.

The output is either a single float (when a single vertex is provided) or a
vector of floats corresponding to the vertex vector. If no vertex vector is
provided, the vector returned corresponds to each vertex in the graph.

### Performance
Because this function must calculate shortest paths for all vertices supplied
in the argument list, it may take a long time.

### Implementation Notes
The eccentricity vector returned by `eccentricity()` may be used as input
for the rest of the distance measures below. If an eccentricity vector is
provided, it will be used. Otherwise, an eccentricity vector will be calculated
for each call to the function. It may therefore be more efficient to calculate,
store, and pass the eccentricities if multiple distance measures are desired.
"""
function eccentricity(
    g::AbstractGraph,
    v::Integer,
    distmx::AbstractMatrix{T} = weights(g)
) where T
    e = maximum(dijkstra_shortest_paths(g, v, distmx).dists)
    e == typemax(T) && error("Infinite path length detected")

    return e
end

eccentricity(
    g::AbstractGraph,
    vs::AbstractVector = vertices(g),
    distmx::AbstractMatrix = weights(g)
) = [eccentricity(g, v, distmx) for v in vs]

eccentricity(g::AbstractGraph, distmx::AbstractMatrix) =
    eccentricity(g, vertices(g), distmx)

"""
    diameter(g, distmx=weights(g))
    diameter(eccentricities)

Given a graph and optional distance matrix, or a vector of precomputed
eccentricities, return the maximum eccentricity of the graph.
"""
diameter(eccentricities::Vector) = maximum(eccentricities)
diameter(g::AbstractGraph, distmx::AbstractMatrix = weights(g)) =
    maximum(eccentricity(g, distmx))

"""
    periphery(g, distmx=weights(g))
    periphery(eccentricities)

Given a graph and optional distance matrix, or a vector of precomputed
eccentricities, return the set of all vertices whose eccentricity is
equal to the graph's diameter (that is, the set of vertices with the
largest eccentricity).
"""
function periphery(eccentricities::Vector)
    diam = maximum(eccentricities)
    return filter(x -> eccentricities[x] == diam, 1:length(eccentricities))
end

periphery(g::AbstractGraph, distmx::AbstractMatrix = weights(g)) =
    periphery(eccentricity(g, distmx))

"""
    radius(g, distmx=weights(g))
    radius(eccentricities)

Given a graph and optional distance matrix, or a vector of precomputed
eccentricities, return the minimum eccentricity of the graph.
"""
radius(eccentricities::Vector) = minimum(eccentricities)
radius(g::AbstractGraph, distmx::AbstractMatrix = weights(g)) =
    minimum(eccentricity(g, distmx))

"""
    center(g, distmx=weights(g))
    center(eccentricities)

Given a graph and optional distance matrix, or a vector of precomputed
eccentricities, return the set of all vertices whose eccentricity is equal
to the graph's radius (that is, the set of vertices with the smallest eccentricity).
"""
function center(eccentricities::Vector)
    rad = radius(eccentricities)
    return filter(x -> eccentricities[x] == rad, 1:length(eccentricities))
end

center(g::AbstractGraph, distmx::AbstractMatrix = weights(g)) =
    center(eccentricity(g, distmx))


text:
#    Problem : GROUP A
#    *********
#    Source: problem 10 in
#    Ph.L. Toint,
#    "Test problems for partially separable optimization and results
#    for the routine PSPMIN",
#    Report 83/4, Department of Mathematics, FUNDP (Namur, B), 1983.
#
#    See also Buckley#116.  Note that MGH#21 is the separable version.
#    SIF input: Ph. Toint, Dec 1989.
#
#    problem 29 in
#    L. Luksan, C. Matonoha and J. Vlcek
#    Modified CUTE problems for sparse unconstrained optimization,
#    Technical Report 1081,
#    Institute of Computer Science,
#    Academy of Science of the Czech Republic
#
#    http://www.cs.cas.cz/matonoha/download/V1081.pdf
#
#    classification SUR2-AN-V-0
#
# # Daniel Henderson, 09/2021

f = (x) -> begin # this is correct
    return 100.0*sum((x[i] - x[i - 1]^2)^2 for i = 1:2:(lastindex(x)-1)) + (1.0 - x[1])^2
end

g! = (g, x) -> begin
    @warn "EXTROSEN g! not implemented"
end

fg! = (g, x) -> begin
    @warn "EXTROSEN g! not implemented"
end

init = (n::Int=1000) -> begin
	n < 2 && @warn("EXTROSNB: number of variables must be ≥ 2")
	n = max(n, 2)

    return n, -1.0*ones(n)
end

#TestSet["EXTROSNB"] = UncProgram("EXTROSNB",  f, g!, fg!, init)


text:
#---------------------------------------------------------------------#
#This function plots the solution
#Written by F.X. Giraldo on April 19 2019
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#
function plot_solution(q,coord,space_method,time)

    #Plot Solution
    text=string(space_method,", time= ","$time")
    plot_handle=plot(coord,q,xlabel="x",ylabel="q(x,t)",legend=true,lw=3,label=[space_method],title=text)
    display(plot_handle)

end


text:
if VERSION < v"1.7.0-DEV.849"
else
    @error "" VERSION >= v"1.7.0-DEV.849"
    error()
end

include("generator_utils.jl")

using REPL
write_dict_to_file("repl_emoji_symbols_below_170.jl",
                   "repl_emoji_symbols_below_170",
                   REPL.REPLCompletions.emoji_symbols)


text:
function to_root(io::IO, root, val)
    if root == 1
        write(io, string(val))
        return
    end
    if root == 3
        print(io, '∛')
    elseif root == 4
        print(io, '∜')
    else
        root != 2 && to_superscript(io, root)
        print(io, '√')
    end
    to_overline(io, string(val))
end
to_root(io::IO, val) = to_root(io, 2, val)

function to_root(root, val)
    sprint() do io
        to_root(io, root, val)
    end
end
to_root(val) = to_root(2, val)