text:
## ------------  ParamInfo

export ParamInfo
# export set_value!

# function ParamInfo(pName :: Symbol, startIdx :: Integer,
#     valueV :: T1, lbV :: T2, ubV :: T2) where {T1 <: AbstractArray, T2 <: Real}

#     sz = size(valueV);
#     return ParamInfo{T1}(pName, startIdx, valueV, fill(lbV, sz), fill(ubV, sz));
# end

function Base.isapprox(p1 :: ParamInfo{F1}, p2 :: ParamInfo{F1};
    atol :: Real = 1e-8) where F1
    return isequal(p1.pName, p2.pName)  &&
        # isapprox(p1.valueV, p2.valueV; atol)  &&
        isapprox(p1.lbV, p2.lbV; atol)  &&
        isapprox(p1.ubV, p2.ubV; atol);
end

Base.isapprox(p1 :: ParamInfo{F1}, p2 :: ParamInfo{F2}) where {F1, F2} = false;

lb(pInfo :: ParamInfo{F1}) where F1 = pInfo.lbV;
ub(pInfo :: ParamInfo{F1}) where F1 = pInfo.ubV;
n_values(pInfo :: ParamInfo{F1}) where F1 = length(pInfo.lbV);
indices(pInfo :: ParamInfo{F1}) where F1 <: Real = pInfo.startIdx;
indices(pInfo :: ParamInfo{F1}) where F1 =
    pInfo.startIdx .+ (0 : (n_values(pInfo) - 1));

random_values(pInfo :: ParamInfo{F1}, rng :: AbstractRNG) where F1 =
    pInfo.lbV .+ (pInfo.ubV .- pInfo.lbV) .* rand(rng, size(pInfo.lbV));
random_values(pInfo :: ParamInfo{F1}, rng :: AbstractRNG) where F1 <: Real =
    pInfo.lbV .+ (pInfo.ubV .- pInfo.lbV) .* rand(rng);

# Reshape an input to match dimensions of ParamInfo.
reshape_vector(pInfo :: ParamInfo{F1}, v :: F1) where F1 <: Real = v;
reshape_vector(pInfo :: ParamInfo{F1}, v) where F1 <: AbstractArray{<:Real} =
    reshape(v, size(pInfo.lbV));

# function set_value!(pInfo :: ParamInfo{F1}, v :: Real) where F1 <: Real
#     pInfo.valueV = v;
# end

# function set_value!(pInfo :: ParamInfo{F1}, v :: AbstractArray{F2}) where {F1 <: Real, F2 <: Real}
#     pInfo.valueV = only(v);
# end

# function set_value!(pInfo :: ParamInfo{<: AbstractArray{F1}}, v :: AbstractArray{F2}) where {F1 <: Real, F2 <: Real}
#     pInfo.valueV = reshape(v, size(pInfo.lbV));
# end

# function value_from_vector!(pInfo :: ParamInfo{F1},
#     valueV :: AbstractVector{F1}) where F1 <: Real
#     pInfo.valueV = valueV[pInfo.startIdx];
# end

# function value_from_vector!(pInfo :: ParamInfo{<:AbstractArray{F1}},
#     valueV :: AbstractVector{F1}) where F1 <: Real
#     pInfo.valueV = reshape(valueV[indices(pInfo)], size(pInfo.lbV));
# end

function make_test_param_info(pName, sz; startIdx = 10)
    if isempty(sz)
        lbV = 0.3;
        ubV = 5.8;
        # valueV = 3.5;
    else
        lbV = fill(0.3, sz);
        ubV = fill(5.8, sz);
        # valueV = 0.3 .* lbV .+ 0.7 .* ubV;
    end
    ParamInfo(pName, startIdx, lbV, ubV);
end

## ------------  ValueVector

# ValueVector
# export ValueVector, set_values!, values, lb, ub

ValueVector{F1}() where F1 = ValueVector{F1}(OrderedDict{Symbol, ParamInfo}());

function make_test_value_vector(n)
    vv = ValueVector{ValueType}();
    startIdx = 1;
    for j = 1 : n
        pName = Symbol("p$j");
        if j == 1
            sz = ();
        elseif j == 2
            sz = (j, );
        else
            sz = (j, 2);
        end
        vv.d[pName] = make_test_param_info(pName, sz; startIdx);
        nElem = prod(sz);
        startIdx += nElem;
    end
    return vv
end

function validate_param_match(pvec :: ParamVector, vv :: ValueVector{F1}) where F1
    isValid = true;
    nCal, _ = n_calibrated_params(pvec);
    if length(vv) != nCal
        isValid = false;
        @warn """
            No of calibrated params differs for $pvec
            pvec: $(nCal)
            vv:   $(length(vv))
            """;
    end

    if isValid  &&  (length(vv) > 0)
        pList = calibrated_params(pvec);
        for p in pList
            pName = name(p);
            if haskey(vv.d, pName)
                pInfo = vv.d[pName];
                if size(calibrated_value(p)) != size(lb(pInfo))
                    isValid = false;
                    @warn "Size mismatch for $pName";
                end
            else
                isValid = false;
                @warn "Param $pName missing from ValueVector";
            end
        end
    end
    return isValid
end

function validate_vv(vv :: ValueVector{F1}) where F1
    isValid = true;
    isValid = isValid  &&  validate_indices(vv);
    return isValid
end

function validate_indices(vv :: ValueVector{F1}) where F1
    isValid = true;
    nValues = n_values(vv);
    if nValues > 0
        takenV = zeros(Int, nValues);
        for (_, pInfo) in vv.d
            # This can be a Range of scalar.
            idxV = indices(pInfo) .- start_index(vv) .+ 1;
            if idxV[1] < 1
                @warn "Indices below start index";
                isValid = false;
            end
            if last(idxV) > nValues
                @warn "Indices past last index";
                isValid = false;
            end
            # Must make this a range for broadcasting to work
            isValid  &&  (takenV[first(idxV) : last(idxV)] .+= 1);
        end
        isValid = all(takenV .== 1);
    end
    return isValid
end

Base.isempty(vv :: ValueVector{F1}) where F1 =
    n_values(vv) == 0;
Base.length(vv :: ValueVector{F1}) where F1 = length(vv.d);

"""
	$(SIGNATURES)

Total number of scalar values across all parameters.
"""
function n_values(vv :: ValueVector{F1}) where F1
    nValues = 0;
    for (_, pInfo) in vv.d
        nValues += n_values(pInfo);
    end
    return nValues
end

function start_index(vv :: ValueVector{F1}) where F1
    (_, pInfo) = first(vv.d);
    return pInfo.startIdx;
end

Base.last(vv :: ValueVector{F1}) where F1 = vv.d[last(collect(keys(vv.d)))];

function last_index(vv :: ValueVector{F1}) where F1
    pInfo = last(vv);
    return last(indices(pInfo))
end

# values(vv :: ValueVector{F1}) where F1 = vv.valueV;
# lb(vv :: ValueVector{F1}) where F1 = vv.lbV;
# ub(vv :: ValueVector{F1}) where F1 = vv.ubV;
# pnames(vv :: ValueVector{F1}) where F1 = vv.pNameV;
# Base.length(vv :: ValueVector{F1}) where F1 = Base.length(vv.valueV);

function Base.isapprox(vv1 :: ValueVector{F1}, vv2 :: ValueVector{F1};
    atol :: Real = 1e-8) where F1

    if length(vv1.d) != length(vv2.d)
        return false;
    end
    for (pName, pInfo1) in vv1.d
        if haskey(vv2.d, pName)
            pInfo2 = vv2.d[pName]
            if !isapprox(pInfo1, pInfo2; atol)
                return false
            end
        else
            return false
        end
    end
    return true
end

# """
# 	$(SIGNATURES)

# Set values from a Vector{Real}. These are usually the transformed values for the entire `ModelObject`.
# """
# function set_values!(vv :: ValueVector{F1}, valueV :: AbstractVector{F1}) where F1
#     @assert length(valueV) >= last_index(vv)  "Vector too short";
#     for (_, pInfo) in vv.d
#         set_value!(pInfo, valueV[indices(pInfo)]);
#     end
#     return nothing
# end

"""
	$(SIGNATURES)

Get *transformed* values for all parameters in the `ValueVector`.
"""
function get_values(pvec :: ParamVector, vv :: ValueVector{F1}) where F1
    # @assert isequal(get_object_id(pvec), get_object_id(vv))  """
    #     ObjectId mismatch:
    #     $(get_object_id(pvec))  vs  $(get_object_id(vv))
    #     """;
    valueV = Vector{F1}();
    for (pName, pInfo) in vv.d
        p = retrieve(pvec, pName);
        @assert size(lb(pInfo)) == size(calibrated_value(p));
        append!(valueV, transform_param(pvec.pTransform, p));
    end

    return valueV
end

function random_guess(vv :: ValueVector{F1}, rng :: AbstractRNG) where F1
    guessV = Vector{F1}();
    for (_, pInfo) in vv.d
        rv = random_values(pInfo, rng);
        if rv isa Real
            push!(guessV, rv);
        else
            append!(guessV, vec(rv));
        end
    end
    return guessV
end

# ----------------


text:
module MLJ

## METHOD IMPORT

# from the Standard Library:
import Distributed: @distributed, nworkers, pmap
import Pkg
import Pkg.TOML

# from the MLJ universe:
using MLJBase
import MLJBase.save
using MLJEnsembles
using MLJTuning
using MLJModels
using OpenML
using MLJSerialization
import MLJSerialization.save # but not re-exported
using MLJIteration
import MLJIteration.IterationControl

using Tables, CategoricalArrays
import Distributions
import Distributions: pdf, mode
import Statistics, StatsBase, LinearAlgebra, Random
import Random: AbstractRNG, MersenneTwister
using ProgressMeter
using ComputationalResources
using ComputationalResources: CPUProcesses

# to be extended:
import MLJBase: fit, update, clean!, fit!, predict, fitted_params,
                show_as_constructed, ==
import MLJModels: models
import ScientificTypes

## CONSTANTS

const srcdir = dirname(@__FILE__)
const TRAITS_NOT_EXPORTED = [
    :supports_online, # no models support this
    :name,            # likely conflict with other uses
    :abstract_type,   # for advanced development only
    :hyperparameter_ranges # not implemented and will probably be deprecated
                           # in favour of hyperparameter priors
]

## INCLUDE FILES

include("version.jl")       # defines MLJ_VERSION constant
include("scitypes.jl")      # extensions to ScientificTypesBase.scitype

## METHOD EXPORT

export MLJ_VERSION

## METHOD RE-EXPORT

# traits for models and measures:
using MLJBase.StatisticalTraits
for trait in setdiff(StatisticalTraits.TRAITS, TRAITS_NOT_EXPORTED)
    :(export $trait) |> eval
end

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

# re-exports from (MLJ)ScientificTypesBase via MLJBase
export Scientific, Found, 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

# re-export from MLJBase:
export nrows, color_off, color_on,
    selectrows, selectcols, restrict, corestrict, complement,
    training_losses,
    predict, predict_mean, predict_median, predict_mode, predict_joint,
    transform, inverse_transform, evaluate, fitted_params, params,
    @constant, @more, HANDLE_GIVEN_ID, UnivariateFinite,
    classes, table, report, rebind!,
    partition, unpack,
    @load_boston, @load_ames, @load_iris, @load_reduced_ames, @load_crabs,
    load_boston, load_ames, load_iris, load_reduced_ames, load_crabs,
    Machine, machine, AbstractNode, @node,
    source, node, fit!, freeze!, thaw!, Node, sources, origins,
    machines, sources, anonymize!, @from_network, fitresults,
    @pipeline, Stack, Pipeline, TransformedTargetModel,
    ResamplingStrategy, Holdout, CV, TimeSeriesCV,
    StratifiedCV, evaluate!, Resampler, iterator,
    default_resource, pretty,
    make_blobs, make_moons, make_circles, make_regression,
    fit_only!, return!, int, decoder

# MLJBase/composition/abstract_types.jl:
for T in vcat(MLJBase.MLJModelInterface.ABSTRACT_MODEL_SUBTYPES,
    MLJBase.COMPOSITE_TYPES,
    MLJBase.SURROGATE_TYPES)
    @eval(export $T)
end
export Surrogate, Composite
export DeterministicNetwork, ProbabilisticNetwork, UnsupervisedNetwork

# MLJBase/measures:
# measure names:
for m in MLJBase.MEASURE_TYPES_ALIASES_AND_INSTANCES
    :(export $m) |> eval
end
export measures,
    aggregate, default_measure, value, skipinvalid,
    roc_curve, roc,
    no_avg, macro_avg, micro_avg

# re-export from MLJEnsembles:
export EnsembleModel

# re-export from MLJTuning:
export Grid, RandomSearch, Explicit, TunedModel, LatinHypercube,
    learning_curve!, learning_curve

# re-export from MLJModels:
export models, localmodels, @load, @iload, load, info,
    ConstantRegressor, ConstantClassifier,     # builtins/Constant.jl
    FeatureSelector, UnivariateStandardizer,   # builtins/Transformers.jl
    Standardizer, UnivariateBoxCoxTransformer,
    OneHotEncoder, ContinuousEncoder, UnivariateDiscretizer,
    FillImputer, matching, BinaryThresholdPredictor

# re-export from MLJIteration:
export MLJIteration
export IteratedModel
for control in MLJIteration.CONTROLS
    eval(:(export $control))
end
export IterationControl

# re-export from MLJOpenML
const MLJOpenML = OpenML

export OpenML, MLJOpenML

# re-export from MLJSerialization:
export Save # control, not method

# re-export from ComputaionalResources:
export CPU1, CPUProcesses, CPUThreads

end # module

ext: jl


text:
# Autogenerated wrapper script for booster_jll for armv6l-linux-gnueabihf
export booster

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("booster")
JLLWrappers.@declare_executable_product(booster)
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_executable_product(
        booster,
        "bin/booster",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using ShowLint
using Test

SL = ShowLint

@testset "ShowLint" begin
    SL.write_configs()
    include("patterns.jl")
    include("repositories.jl")

    SL.clone_repositories(; production=false)
    SL.create_repo_pages(; production=false)
end

max_stars_repo_name: rikhuijzer/display-lint

max_issues_repo_name: rikhuijzer/display-lint

max_forks_repo_name: rikhuijzer/ShowLint


text:
using SymPy: N, Sym, diff, expand, solve, subs

"""
Axisymmetric Solov'ev equilibra with X-point in (R/R₀,Z/R₀,phi) coordinates.
Based on Cerfon & Freidberg, Physics of Plasmas 17, 032502, 2010.

Parameters:
 * `R₀`: position of magnetic axis
 * `B₀`: B-field at magnetic axis
 * `ϵ`:  inverse aspect ratio
 * `κ`:  elongation
 * `δ`:  triangularity
 * `a`:  free constant, determined to match a given beta value
 * `xₛₑₚ`: x position of the X point
 * `yₛₑₚ`: y position of the X point
"""
struct SolovevXpoint{T <: Number} <: AbstractSolovevEquilibrium
    name::String
    R₀::T
    B₀::T
    ϵ::T
    κ::T
    δ::T
    a::T
    xₛₑₚ::T
    yₛₑₚ::T
    c::Vector{T}

    function SolovevXpoint{T}(R₀::T, B₀::T, ϵ::T, κ::T, δ::T, a::T, xₛₑₚ::T, yₛₑₚ::T, c::Vector{T}) where T <: Number
        new("SolovevXpointEquilibrium", R₀, B₀, ϵ, κ, δ, a, xₛₑₚ, yₛₑₚ, c)
    end
end

function SolovevXpoint(R₀::T, B₀::T, ϵ::T, κ::T, δ::T, a::T, xₛₑₚ::T, yₛₑₚ::T) where T <: Number

    n = 12

    x = [Sym("x" * string(i)) for i in 1:3]
    c = [Sym("c" * string(i)) for i in 1:n]

    ψ = (ψ₀(x,a) + c[1]  * ψ₁(x)
                 + c[2]  * ψ₂(x)
                 + c[3]  * ψ₃(x)
                 + c[4]  * ψ₄(x)
                 + c[5]  * ψ₅(x)
                 + c[6]  * ψ₆(x)
                 + c[7]  * ψ₇(x)
                 + c[8]  * ψ₈(x)
                 + c[9]  * ψ₉(x)
                 + c[10] * ψ₁₀(x)
                 + c[11] * ψ₁₁(x)
                 + c[12] * ψ₁₂(x) )

    eqs = [
        expand(subs(subs(ψ, x[1], 1+ϵ), x[2], 0)),
        expand(subs(subs(ψ, x[1], 1-ϵ), x[2], 0)),
        expand(subs(subs(ψ, x[1], 1-δ*ϵ), x[2], κ*ϵ)),
        expand(subs(subs(ψ, x[1], xₛₑₚ), x[2], yₛₑₚ)),
        expand(subs(subs(diff(ψ, x[2]), x[1], 1+ϵ), x[2], 0)),
        expand(subs(subs(diff(ψ, x[2]), x[1], 1-ϵ), x[2], 0)),
        expand(subs(subs(diff(ψ, x[1]), x[1], 1-δ*ϵ), x[2], κ*ϵ)),
        expand(subs(subs(diff(ψ, x[1]), x[1], xₛₑₚ), x[2], yₛₑₚ)),
        expand(subs(subs(diff(ψ, x[2]), x[1], xₛₑₚ), x[2], yₛₑₚ)),
        expand(subs(subs(diff(ψ, x[2], 2), x[1], 1+ϵ), x[2], 0) - (1 + asin(δ))^2 / (ϵ * κ^2) * subs(subs(diff(ψ, x[1]), x[1], 1+ϵ), x[2], 0)),
        expand(subs(subs(diff(ψ, x[2], 2), x[1], 1-ϵ), x[2], 0) + (1 - asin(δ))^2 / (ϵ * κ^2) * subs(subs(diff(ψ, x[1]), x[1], 1-ϵ), x[2], 0)),
        expand(subs(subs(diff(ψ, x[1], 2), x[1], 1-δ*ϵ), x[2], κ*ϵ) - κ / (ϵ * (1 - δ^2)) * subs(subs(diff(ψ, x[2]), x[1], 1-δ*ϵ), x[2], κ*ϵ))
    ]

    csym = solve(eqs, c)
    cnum = [N(csym[c[i]]) for i in 1:n]

    SolovevXpoint{T}(R₀, B₀, ϵ, κ, δ, a, xₛₑₚ, yₛₑₚ, cnum)
end

SolovevXpointITER() = SolovevXpoint(6.2, 5.3, 0.32, 1.7, 0.33, -0.155, 0.88, -0.60)
SolovevXpointNSTX() = SolovevXpoint(0.85, 0.3, 0.78, 2.0, 0.35, -0.05, 0.70, -1.71)

function Base.show(io::IO, equ::SolovevXpoint)
    print(io, "Solovev Xpoint Equilibrium with\n")
    print(io, "  R₀ = ", equ.R₀, "\n")
    print(io, "  B₀ = ", equ.B₀, "\n")
    print(io, "  ϵ  = ", equ.ϵ, "\n")
    print(io, "  κ  = ", equ.κ, "\n")
    print(io, "  δ  = ", equ.δ, "\n")
    print(io, "  a  = ", equ.a, "\n")
    print(io, "  xₛₑₚ  = ", equ.xₛₑₚ, "\n")
    print(io, "  yₛₑₚ  = ", equ.yₛₑₚ)
end

function A₃(x::AbstractArray{T,1}, equ::SolovevXpoint) where {T <: Number}
    (ψ₀(x, equ.a) + equ.c[1]  * ψ₁(x)
                  + equ.c[2]  * ψ₂(x)
                  + equ.c[3]  * ψ₃(x)
                  + equ.c[4]  * ψ₄(x)
                  + equ.c[5]  * ψ₅(x)
                  + equ.c[6]  * ψ₆(x)
                  + equ.c[7]  * ψ₇(x)
                  + equ.c[8]  * ψ₈(x)
                  + equ.c[9]  * ψ₉(x)
                  + equ.c[10] * ψ₁₀(x)
                  + equ.c[11] * ψ₁₁(x)
                  + equ.c[12] * ψ₁₂(x) )
end

macro solovev_xpoint_equilibrium(R₀, B₀, ϵ, κ, δ, a, xₛₑₚ, yₛₑₚ)
    generate_equilibrium_code(SolovevXpoint(R₀, B₀, ϵ, κ, δ, a, xₛₑₚ, yₛₑₚ); output=false)
end


text:
# Parsed file utility
"""
A cache of file content, parsed with an arbitrary parser function.

This is a modified and generalized version of `Base.CachedTOMLDict`.

Getting the value of the cache with `f[]` will automatically update the parsed
value whenever the file changes.
"""
mutable struct CachedParsedFile{T}
    path::String
    inode::UInt64
    mtime::Float64
    size::Int64
    hash::Vector{UInt8}
    d::T      # The cached value.
    parser    # Parser function used to update the cached value
end

function CachedParsedFile{T}(parser::Function, path::String) where T
    s = stat(path)
    if ispath(s)
        content = read(path)
        hash = sha1(content)
    else
        hash = UInt8[]
        content = nothing
    end
    d = parser(content)
    return CachedParsedFile{T}(
        path,
        s.inode,
        s.mtime,
        s.size,
        hash,
        d,
        parser,
    )
end

function Base.getindex(f::CachedParsedFile)
    s = stat(f.path)
    time_since_cached = time() - f.mtime
    rough_mtime_granularity = 0.1 # seconds
    # In case the file is being updated faster than the mtime granularity,
    # and have the same size after the update we might miss that it changed. Therefore
    # always check the hash in case we recently created the cache.
    if #==# time_since_cached < rough_mtime_granularity ||
            s.inode != f.inode || s.mtime != f.mtime || f.size != s.size
        if ispath(s)
            content = read(f.path)
            new_hash = sha1(content)
        else
            content = nothing
            new_hash = UInt8[]
        end
        if new_hash != f.hash
            f.inode = s.inode
            f.mtime = s.mtime
            f.size = s.size
            f.hash = new_hash
            @debug "Cache of file $(repr(f.path)) invalid, reparsing..."
            return f.d = f.parser(content)
        end
    end
    return f.d
end

function Base.show(io::IO, m::MIME"text/plain", f::CachedParsedFile)
    println(io, "Cache of file $(repr(f.path)) with value")
    show(io, m, f[])
end

# Parse Data.toml into DataProject which updates when the file does.
function parse_and_cache_project(sys_path::AbstractString)
    sys_data_dir = dirname(sys_path)
    CachedParsedFile{DataProject}(sys_path) do content
        if isnothing(content)
            DataProject()
        else
            _load_project(String(content), sys_data_dir)
        end
    end
end

#------------------------------------------------------------------------------
abstract type AbstractTomlFileDataProject <: AbstractDataProject end

function Base.get(proj::AbstractTomlFileDataProject, name::AbstractString, default)
    get(_get_cached(proj), name, default)
end

function Base.keys(proj::AbstractTomlFileDataProject)
    keys(_get_cached(proj))
end

function Base.iterate(proj::AbstractTomlFileDataProject, state=nothing)
    # This is a little complex because we want iterate to work even if the
    # active project changes concurrently, which means wrapping up the initial
    # result of _get_cached with the iterator state.
    if isnothing(state)
        cached_values = values(_get_cached(proj))
        if isnothing(cached_values)
            return nothing
        end
        wrapped_itr = iterate(cached_values)
    else
        (cached_values, wrapped_state) = state
        wrapped_itr = iterate(cached_values, wrapped_state)
    end
    if isnothing(wrapped_itr)
        return nothing
    else
        (data, wrapped_state) = wrapped_itr
        (data, (cached_values, wrapped_state))
    end
end

Base.pairs(proj::AbstractTomlFileDataProject) = pairs(_get_cached(proj))

data_drivers(proj::AbstractTomlFileDataProject) = data_drivers(_get_cached(proj))

#-------------------------------------------------------------------------------
"""
Data project which automatically updates based on a TOML file on the local
filesystem.
"""
mutable struct TomlFileDataProject <: AbstractTomlFileDataProject
    path::String
    cache::CachedParsedFile{DataProject}
end

function TomlFileDataProject(path::String)
    cache = parse_and_cache_project(path)
    TomlFileDataProject(path, cache)
end

function _get_cached(proj::TomlFileDataProject)
    proj.cache[]
end

project_name(proj::TomlFileDataProject) = proj.path

#------------------------------------------------------------------------------
"""
Data project, based on the location of the current explicitly selected Julia
Project.toml, as reported by `Base.active_project(false)`.

Several factors make the implementation a bit complicated:
  * The active project may change at any time without warning
  * The active project may be `nothing` when no explicit project is selected
  * There might be no Data.toml for the active project
  * The user can change Data.toml interactively and we'd like that to be
    reflected within the program.
"""
mutable struct ActiveDataProject <: AbstractTomlFileDataProject
    active_project_path::Union{Nothing,String} # Detects when Base.active_project changes
    cache::Union{DataProject,CachedParsedFile{DataProject}}
end

function ActiveDataProject()
    proj = ActiveDataProject(nothing, DataProject())
    _get_cached(proj)
    proj
end

function _active_project_data_toml(project_path=Base.active_project(false))
    isnothing(project_path) ?
        nothing :
        joinpath(dirname(project_path), "Data.toml")
end

function _get_cached(proj::ActiveDataProject)
    active_project = Base.active_project(false)
    if proj.active_project_path != active_project
        # The unusual case: active project has changed.
        if isnothing(active_project)
            proj.cache = DataProject()
        else
            data_toml = _active_project_data_toml(active_project)
            # Need to re-cache
            proj.cache = parse_and_cache_project(data_toml)
        end
        proj.active_project_path = active_project
    end
    proj.cache isa DataProject ? proj.cache : proj.cache[]
end

project_name(::ActiveDataProject) = _active_project_data_toml()

#-------------------------------------------------------------------------------

function _load_project(content::AbstractString, sys_data_dir)
    toml_str = _fill_template(sys_data_dir, content)
    config = TOML.parse(toml_str)
    load_project(config)
end


text:
function print_iter_head(;with_gnorm::Bool=true)
    if with_gnorm
        @printf("%5s   %12s   %12s   %12s   %12s\n",
            "Iter", "f.value", "f.change", "g.norm", "step")
        println("==================================================================")
    else
        @printf("%5s   %12s   %12s   %12s\n",
            "Iter", "f.value", "f.change", "step")
        println("======================================================")
    end
end

function print_iter(t::Int, v::Real)
    @printf("%5d   %12.4e\n", t, v)
end

function print_iter(t::Int, v::Real, v_pre::Real, g::StridedArray, α::Real)
    @printf("%5d   %12.4e   %12.4e   %12.4e   %12.4e\n", t, v, v - v_pre, vecnorm(g), α)
end

function print_iter(t::Int, v::Real, v_pre::Real, α::Real)
    @printf("%5d   %12.4e   %12.4e   %12.4e\n", t, v, v - v_pre, α)
end

function print_final(t::Int, v::Real, converged::Bool)
    if converged
        @printf("Converged with %d iterations @ f.value = %.4e\n", t, v)
    else
        @printf("Terminated without convergence with %d iterations @ f.value = %.4e\n", t, v)
    end
end


text:
using StatsBase
export get_aggregate_adj_mat

function StatsBase.counts(data::Vector{T}, ref::Vector{T}) where {T}
    mapper = Dict(val=>i for (i,val) in enumerate(ref))
    out = zeros(Int, length(ref))
    for x in data
        out[mapper[x]] += 1
    end
    return out
end

function get_aggregate_adj_mat(data::Vector{Vector{T}}, ref::Vector{T}) where {T}
    mapper = Dict(val=>i for (i,val) in enumerate(ref))
    A = zeros(Int, length(ref), length(ref))
    for path in data
        for i in Iterators.rest(eachindex(path),1)
            A[mapper[path[i-1]], mapper[path[i]]] += 1
        end
    end
    return A
end


text:
function schedule(
    batch::Batch,
    scheduled_jobs::Vector{Job},
    machines::Vector{Machine},
    order_time::DateTime,
    beta;
    threads = 1,
    max_solutions = -1,
    big_m = 1e8,
)

    ################
    # Model
    ################

    scheduler = Model(Cbc.Optimizer)
    set_optimizer_attribute(scheduler, "threads", threads)
    set_optimizer_attribute(scheduler, "preProcess", "off")
    if max_solutions > 0
        set_optimizer_attribute(scheduler, "maxSolutions", max_solutions)
    end

    ################
    # Constants
    ################

    # machines
    M = length(machines)
    T = [m.T for m in machines]

    # operations
    N = get_N(batch)
    P = get_P(batch)
    C = get_C(batch)
    tau = get_tau(batch)

    # TCMB
    constraints = get_constraints(batch)

    # convert order_time to epoch sec
    t0 = Int64(Dates.datetime2epochms(order_time) / 1000)

    # scheduled operations
    scheduled_operations = vcat([job.operations for job in scheduled_jobs]...)
    N_prime = length(scheduled_operations)
    tau_prime = [op.tau for op in scheduled_operations]
    S_prime = [op.S - t0 for op in scheduled_operations]
    E_prime = [op.E for op in scheduled_operations]

    ################
    # Variables
    ################

    # start time
    @variable(scheduler, S[1:N] >= 0)

    # binary notation of which machine to process an operation
    F = Dict()
    for a = 1:N, m = 1:M
        if C[a] == T[m]
            F[a, m] = @variable(scheduler, binary = true)
        end
    end

    # shared-machine precedence relationship
    Q = Dict()
    for a = 1:N, b = 1:N, m = 1:M
        if a != b && C[a] == C[b] == T[m]
            Q[a, b, m] = @variable(scheduler, binary = true)
        end
    end

    # shared-machine precedence between a scheduled and an yet-to-be-scheduled operation
    R_precede = Dict()
    R_follow = Dict()
    for a = 1:N, b_prime = 1:N_prime
        if T[E_prime[b_prime]] == C[a]
            R_precede[a, b_prime] = @variable(scheduler, binary = true)
            R_follow[a, b_prime] = @variable(scheduler, binary = true)
        end
    end

    # variable for detecting lastly processed operation's end time
    @variable(scheduler, Omega, Int)

    ################
    # Constraints
    ################

    # 1-1: the end time of the entire schedule must larger than the end time of any operation
    for a = 1:N
        @constraint(scheduler, Omega >= S[a] + tau[a])
    end

    # 2-1: an operation is processed by and only by one machine
    for a = 1:N
        compatible = filter(m -> C[a] == T[m], 1:M)
        @constraint(scheduler, sum(F[a, m] for m in compatible) == 1)
    end

    # 2-2: a shared-machine precedence relationship exists between each pair of operations sharing the common machine
    for m = 1:M, a = 1:N, b = a:N
        if a != b && C[a] == C[b] == T[m]
            # AND circuit
            @constraint(scheduler, Q[a, b, m] + Q[b, a, m] <= F[a, m])
            @constraint(scheduler, Q[a, b, m] + Q[b, a, m] <= F[b, m])
            @constraint(scheduler, Q[a, b, m] + Q[b, a, m] >= F[a, m] + F[b, m] - 1)
        end
    end

    # 2-3: the shared-machine precedence relationship between a certain pair of operations can exist for at most one machine
    for a = 1:N, b = 1:N
        if C[a] == C[b]
            compatible = filter(m -> a != b && C[a] == C[b] == T[m], 1:M)
            @constraint(scheduler, sum(Q[a, b, m] for m in compatible) <= 1)
        end
    end

    # 2-4: a shared-machine precedence relationship must not violate the operation dependency
    for a = 1:N, b = 1:N, m = 1:M
        if P[a, b] == 1 && C[a] == C[b] == T[m]
            @constraint(scheduler, Q[b, a, m] == 0)
        end
    end

    # 2-5: a machine can process at most one operation at a time
    for a = 1:N, b = 1:N, m = 1:M
        if a != b && C[a] == C[b] == T[m]
            @constraint(scheduler, S[a] + tau[a] + beta <= S[b] + big_m * (1 - Q[a, b, m]))
        end
    end

    # 2-6: the dependency among operations must hold
    for a = 1:N, b = 1:N
        if P[a, b] == 1
            @constraint(scheduler, S[a] + tau[a] <= S[b])
        end
    end

    # 3-1: the absolute value of the difference between two operation boundaries is less than or equal to the maximum
    #      tolerable difference. There are four cases depending on the combination of boundaries
    for con in constraints
        a = con.operation_id_1
        b = con.operation_id_2
        if con.boundary_1 == :start && con.boundary_2 == :start
            @constraint(scheduler, S[a] - S[b] <= con.alpha)
            @constraint(scheduler, S[b] - S[a] <= con.alpha)
        elseif con.boundary_1 == :start && con.boundary_2 == :end
            @constraint(scheduler, S[a] - (S[b] + tau[b]) <= con.alpha)
            @constraint(scheduler, (S[b] + tau[b]) - S[a] <= con.alpha)
        elseif con.boundary_1 == :end && con.boundary_2 == :start
            @constraint(scheduler, (S[a] + tau[a]) - S[b] <= con.alpha)
            @constraint(scheduler, S[b] - (S[a] + tau[a]) <= con.alpha)
        elseif con.boundary_1 == :end && con.boundary_2 == :end
            @constraint(scheduler, (S[a] + tau[a]) - (S[b] + tau[b]) <= con.alpha)
            @constraint(scheduler, (S[b] + tau[b]) - (S[a] + tau[a]) <= con.alpha)
        end
    end

    # 4-1: an shared-machine precedence or following relationship must exist for any pair of operations sharing a
    #      common machine to process them, one of which is in a previously scheduled job and the other is in a new job
    for a = 1:N, b_prime = 1:N_prime
        if T[E_prime[b_prime]] == C[a]
            @constraint(scheduler, R_precede[a, b_prime] + R_follow[a, b_prime] == F[a, E_prime[b_prime]])
        end
    end

    # 4-2: a machine can process at most one operation at a time (case: an operation in a new job precedes another
    #      operation in a previously scheduled job)
    for a = 1:N, b_prime = 1:N_prime
        if T[E_prime[b_prime]] == C[a]
            @constraint(
                scheduler,
                S[a] + tau[a] + beta <= S_prime[b_prime] + big_m * (1 - R_precede[a, b_prime])
            )
        end
    end

    # 4-3: a machine can process at most one operation at a time (case: an operation in a new job follows another
    #      operation in a previously scheduled job)
    for a = 1:N, b_prime = 1:N_prime
        if T[E_prime[b_prime]] == C[a]
            @constraint(
                scheduler,
                S_prime[b_prime] + tau_prime[b_prime] + beta <= S[a] + big_m * (1 - R_follow[a, b_prime])
            )
        end
    end

    ################
    # Objective
    ################

    @objective(scheduler, Min, Omega)

    ################
    # Optimize
    ################

    JuMP.optimize!(scheduler)

    ################
    # Out
    ################

    # Time to start each operation
    S = [Int64(round(JuMP.value(S[a]))) + t0 for a = 1:N] # rounding might cause bug in the future

    # Machines selected for each operation
    E = zeros(Int64, N)
    for a = 1:N, m = 1:M
        if C[a] == T[m] && JuMP.value(F[a, m]) == 1
            E[a] = m
        end
    end

    set_schedule(batch, S, E)

    append!(scheduled_jobs, batch.jobs)

    return JuMP.value(Omega)
end


text:
include("column_names.jl")

function demographic_inertia_simulation(initial_values,
                                        runs,
                                        duration)

    # build the data array
    initial_values["duration"] = duration # TODO bad programming style. Need to fix this.
    res = build_simulation_results_array(runs, duration)
    populate_initial_values!(res, initial_values)
    run_model_with_drift!(res, initial_values, runs)
    @bp
    res
end

function build_simulation_results_array(runs,
                                        duration)
    num_rows = runs * duration
    int_df = DataFrame(repeat([Union{Int64, Missing}], size(INTEGER_COLUMNS)[1]),
                INTEGER_COLUMNS,
                runs * duration)

    float_df = DataFrame(repeat([Union{Float64, Missing}], size(FLOAT_COLUMNS)[1]),
                FLOAT_COLUMNS,
                runs * duration)
    #df[!, :scen] .= Union{String, Missing}
    text_df = DataFrame(repeat([Union{String, Missing}], size(TEXT_COLUMNS)[1]),
                TEXT_COLUMNS,
                runs * duration)
    text_df[!, :date_time] .= Dates.format(Dates.now(), "mm/dd/yyyy HH:MM")
    hcat(int_df, float_df, text_df)

end

function populate_initial_values!(ra, iv)
    ra[1, :f1] = iv["number_of_females_1"]
    ra[1, :f2] = iv["number_of_females_2"]
    ra[1, :f3] = iv["number_of_females_3"]
    ra[1, :m1] = iv["number_of_males_1"]
    ra[1, :m2] = iv["number_of_males_2"]
    ra[1, :m3] = iv["number_of_males_3"]
    ra[1, :lev1] = iv["number_of_females_1"] + iv["number_of_males_1"]
    ra[1, :lev2] = iv["number_of_females_2"] + iv["number_of_males_2"]
    ra[1, :lev3] = iv["number_of_females_3"] + iv["number_of_males_3"]
    ra[1, :f] = iv["number_of_females_1"] + iv["number_of_females_2"] + iv["number_of_females_3"]
    ra[1, :m] = iv["number_of_males_1"] + iv["number_of_males_2"] + iv["number_of_males_3"]
    ra[1, :fpct1] = round(ra[1, :f1]/ra[1, :lev1], digits=3)
    ra[1, :fpct2] = round(ra[1, :f2]/ra[1, :lev2], digits=3)
    ra[1, :fpct3] = round(ra[1, :f3]/ra[1, :lev3], digits=3)
    ra[1, :mpct1] = round(1 - ra[1, :fpct1], digits=3)
    ra[1, :mpct2] = round(1 - ra[1, :fpct2], digits=3)
    ra[1, :mpct3] = round(1 - ra[1, :fpct3], digits=3)
    ra[1, :fpct] = round(sum(ra[1,:f])/sum(ra[1, ALL_FACULTY]), digits=3)
    ra[1, :deptn] = sum(ra[1, ALL_FACULTY])
    ra[1, :attr] = 0
    ra[1, :run] = 0
    ra[1, :yr] = 0
    ra[1, :year] = 0
    ra[1, :hire] = 0
    ra[1, :r_attr] = 0.
    ra[1, :r_hire] = 0.
    ra[1, :r_hire_adj] = 0.
    ra[1, :f1_1d] = 0
    ra[1, :f2_1d] = 0
    ra[1, :f3_1d] = 0
    ra[1, :m1_1d] = 0
    ra[1, :m2_1d] = 0
    ra[1, :m3_1d] = 0
    ra[1, :deptn_1d] = 0
    ra[:, :r_fhire1] .= iv["hiring_rate_women_1"]
    ra[:, :r_fhire2] .= iv["hiring_rate_women_2"]
    ra[:, :r_fhire3] .= iv["hiring_rate_women_3"]
    ra[:, :r_mhire1] .= iv["hiring_rate_men_1"]
    ra[:, :r_mhire2] .= iv["hiring_rate_men_2"]
    ra[:, :r_mhire3] .= iv["hiring_rate_men_3"]
    ra[:, :r_fattr1] .= iv["attrition_rate_women_1"]
    ra[:, :r_fattr2] .= iv["attrition_rate_women_2"]
    ra[:, :r_fattr3] .= iv["attrition_rate_women_3"]
    ra[:, :r_mattr1] .= iv["attrition_rate_men_1"]
    ra[:, :r_mattr2] .= iv["attrition_rate_men_2"]
    ra[:, :r_mattr3] .= iv["attrition_rate_men_3"]
    ra[:, :r_fprom1] .= iv["female_promotion_probability_1"]
    ra[:, :r_fprom2] .= iv["female_promotion_probability_2"]
    ra[:, :r_mprom1] .= iv["male_promotion_probability_1"]
    ra[:, :r_mprom2] .= iv["male_promotion_probability_2"]
    ra[:, :prom3] .= 0
    ra[:, :r_fprom3] .= 0.
    ra[:, :r_mprom3] .= 0.
    ra[:, :ss_deptn_ub] .= iv["upperbound"]
    ra[:, :ss_deptn_lb] .= iv["lowerbound"]
    ra[:, :ss_deptn_range] .= ra[:, :ss_deptn_ub] - ra[:, :ss_deptn_lb]
    ra[:, :ss_duration] .= iv["duration"]
    ra[:, :model] .= iv["model_name"]
    ra[:, :language] .= "julia"
    ra[:, :scen] .= "scenario name"
    ra[:, :notes] .= "notes"
end

function run_model_with_drift!(df, initial_values, runs)

    department_size_forecasts = calculate_yearly_dept_size_targets(df[1, :ss_duration], initial_values["growth_rate"])

    dept_attr_rate = mean(df[1, [:r_fattr1,
                                :r_fattr2,
                                :r_fattr3,
                                :r_mattr1,
                                :r_mattr2,
                                :r_mattr3 ]])
    base_hiring_rate = dept_attr_rate

    df[1, :g_yr_rate] = department_size_forecasts[1] # TODO move this within the loop to set initial growth value
    df[:, :runn] .= runs

    @showprogress for j in 1:runs

        for i in 2:df[1,:ss_duration]
            k = i + (j-1)*df[1, :ss_duration]
            if i == 2
                df[k-1, :f1] = df[1, :f1]
                df[k-1, :f2] = df[1, :f2]
                df[k-1, :f3] = df[1, :f3]
                df[k-1, :m1] = df[1, :m1]
                df[k-1, :m2] = df[1, :m2]
                df[k-1, :m3] = df[1, :m3]
                df[k-1, :lev1] = df[1, :lev1]
                df[k-1, :lev2] = df[1, :lev2]
                df[k-1, :lev3] = df[1, :lev3]
                df[k-1, :f] = df[1, :f]
                df[k-1, :m] = df[1, :m]
                df[k-1, :fpct1] = df[1, :fpct1]
                df[k-1, :fpct2] = df[1, :fpct2]
                df[k-1, :fpct3] = df[1, :fpct3]
                df[k-1, :mpct1] = df[1, :mpct1]
                df[k-1, :mpct2] = df[1, :mpct2]
                df[k-1, :mpct3] = df[1, :mpct3]
                df[k-1, :fpct] = df[1, :fpct]
                df[k-1, :deptn] = df[1, :deptn]
                df[k-1, :attr] = df[1, :attr]
                df[k-1, :run] = j-1
                df[k-1, :yr] = i-2 # sets yr to 0.
                df[k-1, :year] = i-2
                df[k-1, :r_attr] = df[1, :attr]
                df[k-1, :r_hire] = df[1, :hire]
                df[k-1, :r_hire_adj] = df[1, :r_hire_adj]
                df[k-1, :f1_1d] = df[1, :f1_1d]
                df[k-1, :f2_1d] = df[1, :f2_1d]
                df[k-1, :f3_1d] = df[1, :f3_1d]
                df[k-1, :m1_1d] = df[1, :m1_1d]
                df[k-1, :m2_1d] = df[1, :m2_1d]
                df[k-1, :m3_1d] = df[1, :m3_1d]
                df[k-1, :deptn_1d] = df[1, :deptn_1d]
            end

            df[k, :yr] = i - 1
            df[k, :year] = i - 1
            df[k, :run] = j - 1

            # process attritions
            df[k, :fattr1] = rand(Binomial(df[k-1, :f1], df[k-1, :r_fattr1]), 1)[1]
            df[k, :fattr2] = rand(Binomial(df[k-1, :f2], df[k-1, :r_fattr2]), 1)[1]
            df[k, :fattr3] = rand(Binomial(df[k-1, :f3], df[k-1, :r_fattr3]), 1)[1]
            df[k, :mattr1] = rand(Binomial(df[k-1, :m1], df[k-1, :r_mattr1]), 1)[1]
            df[k, :mattr2] = rand(Binomial(df[k-1, :m2], df[k-1, :r_mattr2]), 1)[1]
            df[k, :mattr3] = rand(Binomial(df[k-1, :m3], df[k-1, :r_mattr3]), 1)[1]

            # update model numbers after attrition
            df[k, :f1] = df[k-1, :f1 ] - df[k, :fattr1]
            df[k, :f2] = df[k-1, :f2 ] - df[k, :fattr2]
            df[k, :f3] = df[k-1, :f3 ] - df[k, :fattr3]
            df[k, :m1] = df[k-1, :m1 ] - df[k, :mattr1]
            df[k, :m2] = df[k-1, :m2 ] - df[k, :mattr2]
            df[k, :m3] = df[k-1, :m3 ] - df[k, :mattr3]

            # promotion process

            df[k, :fprom2] = rand(Binomial(df[k, :f2], df[k-1, :r_fprom2]), 1)[1]
            df[k, :fprom1] = rand(Binomial(df[k, :f1], df[k-1, :r_fprom1]), 1)[1]
            df[k, :mprom2] = rand(Binomial(df[k, :m2], df[k-1, :r_mprom2]), 1)[1]
            df[k, :mprom1] = rand(Binomial(df[k, :m1], df[k-1, :r_mprom1]), 1)[1]

            df[k, :prom1] = sum(df[k, [:fprom1, :mprom1]])
            df[k, :prom2] = sum(df[k, [:fprom2, :mprom2]])
            df[k, :prom3] = 0
            df[k, :fprom] = sum(df[k, [:fprom1, :fprom2]])
            df[k, :mprom] = sum(df[k, [:mprom1, :mprom2]])
            df[k, :prom] = sum(df[k, [:fprom, :mprom]])

            # update model numbers after promotion
            # add promotions to senior level
            df[k, :f2] += df[k, :fprom1]
            df[k, :f3] += df[k, :fprom2]
            df[k, :m2] += df[k, :mprom1]
            df[k, :m3] += df[k, :mprom2]

            # remove promoted folks from previous levels
            df[k, :f1] -= df[k, :fprom1]
            df[k, :f2] -= df[k, :fprom2]
            df[k, :m1] -= df[k, :mprom1]
            df[k, :m2] -= df[k, :mprom2]

            # hiring new faculty
            base_hr_with_growth = base_hiring_rate + department_size_forecasts[i]
            lower_bound_adjustment = base_hr_with_growth/(1 + exp(df[k-1, :deptn] - df[k, :ss_deptn_lb]))
            upper_bound_adjustment = base_hr_with_growth/(1 + exp(df[k, :ss_deptn_ub] - df[k-1, :deptn]))
            dept_adjusted_hiring_rate = max(base_hr_with_growth + lower_bound_adjustment - upper_bound_adjustment, 0)
            df[k, :r_attr] = dept_attr_rate
            df[k, :r_hire] = base_hiring_rate
            df[k, :r_hire_adj] = dept_adjusted_hiring_rate

            total_vacancies = rand(Binomial(df[k-1, :deptn], dept_adjusted_hiring_rate), 1)[1]
            hiring_probabilities = [df[k, :r_fhire1],
                                    df[k, :r_fhire2],
                                    df[k, :r_fhire3],
                                    df[k, :r_mhire1],
                                    df[k, :r_mhire2],
                                    df[k, :r_mhire3]]

            hires = rand(Multinomial(total_vacancies,(hiring_probabilities/sum(hiring_probabilities))), 1)

            df[k, :f1] += hires[1]
            df[k, :f2] += hires[2]
            df[k, :f3] += hires[3]
            df[k, :m1] += hires[4]
            df[k, :m2] += hires[5]
            df[k, :m3] += hires[6]

            df[k, :fhire1] = hires[1]
            df[k, :fhire2] = hires[2]
            df[k, :fhire3] = hires[3]
            df[k, :mhire1] = hires[4]
            df[k, :mhire2] = hires[5]
            df[k, :mhire3] = hires[6]

            df[k, :hire1] = sum(df[k, [:fhire1, :mhire1]])
            df[k, :hire2] = sum(df[k, [:fhire2, :mhire2]])
            df[k, :hire3] = sum(df[k, [:fhire3, :mhire3]])
            df[k, :mhire] = sum(df[k, [:mhire1, :mhire2, :mhire3]])
            df[k, :fhire] = sum(df[k, [:fhire1, :fhire2, :fhire3]])
            df[k, :hire] = sum(hires)

            # save system state
            df[k, :attr1] = sum(df[k, [:fattr1, :mattr1]])
            df[k, :attr2] = sum(df[k, [:fattr2, :mattr2]])
            df[k, :attr3] = sum(df[k, [:fattr3, :mattr3]])
            df[k, :fattr] = sum(df[k, [:fattr1, :fattr2, :fattr3]])
            df[k, :mattr] = sum(df[k, [:mattr1, :mattr2, :mattr3]])
            df[k, :attr] = sum(df[k, [:fattr, :mattr]])

            df[k, :fpct] = round(sum(df[k, FEMALE_COLUMNS])/sum(df[k, ALL_FACULTY]), digits=3)
            df[k, :fpct1] = round(df[k, :f1]/sum(df[k, [:f1, :m1]]), digits=3)
            df[k, :fpct2] = round(df[k, :f2]/sum(df[k, [:f2, :m2]]), digits=3)
            df[k, :fpct3] = round(df[k, :f3]/sum(df[k, [:f3, :m3]]), digits=3)
            df[k, :mpct1] = round(df[k, :m1]/sum(df[k, [:f1, :m1]]), digits=3)
            df[k, :mpct2] = round(df[k, :m2]/sum(df[k, [:f2, :m2]]), digits=3)
            df[k, :mpct3] = round(df[k, :m3]/sum(df[k, [:f3, :m3]]), digits=3)

            df[k, :deptn] = sum(df[k, ALL_FACULTY])
            df[k, :lev1] = sum(df[k, [:f1, :m1]])
            df[k, :lev2] = sum(df[k, [:f2, :m2]])
            df[k, :lev3] = sum(df[k, [:f3, :m3]])
            df[k, :f] = sum(df[k, FEMALE_COLUMNS])
            df[k, :m] = sum(df[k, MALE_COLUMNS])

            df[k, :g_yr_rate] = department_size_forecasts[i]
            df[k, :f1_1d] = df[k, :f1] - df[k-1, :f1]
            df[k, :f2_1d] = df[k, :f2] - df[k-1, :f2]
            df[k, :f3_1d] = df[k, :f3] - df[k-1, :f3]
            df[k, :m1_1d] = df[k, :m1] - df[k-1, :m1]
            df[k, :m2_1d] = df[k, :m2] - df[k-1, :m2]
            df[k, :m3_1d] = df[k, :m3] - df[k-1, :m3]
            df[k, :deptn_1d] = df[k, :deptn] - df[k-1, :deptn]
        end
    end
    CSV.write("model_simulation_detail_$(Dates.format(Dates.now(), "mm-dd-yyyy HH:MM")).csv", df)
    df
end

function _calculate_forecast_interval(duration, forecasts_list)
    convert(UInt16, ceil(duration/length(forecasts_list)))
end

function calculate_yearly_dept_size_targets(duration, forecasts_list)
    repeat_interval = _calculate_forecast_interval(duration, forecasts_list)
    repeat(forecasts_list, inner=repeat_interval)
end

function calculate_simulation_summary(sim; criterion = 0.25)
    # subset the dataframe with the columns that are to be summarized

    reduced_df = select(sim, Not(unique(vcat(UNSUMMARIZED_COLUMNS, DECORATOR_COLUMNS, :run))))
    unsummarized_df = select(sim, unique(vcat(UNSUMMARIZED_COLUMNS, :yr)))
    text_df = select(sim, vcat(TEXT_COLUMNS, :yr))
    # aggregate over the dataframe

    summary = aggregate(reduced_df, :yr,
                        [mean, std, _025, _975])
    summary = replace_nans.(summary)
    summary = round.(summary, digits=3)
    p_fpct = by(sim, :yr, :fpct =>
            x -> p_fpct = 1 - sum(x .<= criterion)/length(x))
    p_fpct = rename(select(p_fpct, :fpct_function), [:p_fpct])
    summary = hcat(summary, p_fpct)

    unsummarized = aggregate(unsummarized_df, :yr, mean)
    unsummarized = select(unsummarized, Not(:yr))
    unsummarized = update_unsummarized_column_names!(unsummarized)
    unsummarized = round.(unsummarized, digits=3)

    text_summary = aggregate(text_df, :yr, return_first_element)
    text_summary = select(text_summary, Not(:yr))
    text_summary = update_text_column_names!(text_summary)

    summary = hcat(summary, unsummarized)
    summary = update_column_names!(summary)
    summary = hcat(summary, text_summary)
    # create non_summarized colums, but of the same length as the summarized dataframe

    # concatenate the dataframes
    summary = concatenate_actual_data_to_simulation!(summary)

    CSV.write("model_summary_$(Dates.format(Dates.now(), "mm-dd-yyyy HH:MM")).csv", summary)

    simulation = replace_nans.(sim)
    return simulation, summary
end

# function calculate_probability_of_target(criterion::Float64, df::DataFrame)
#     df2 = @byrow! df begin
#         @newcol p_fpct::Array{Float64}
#         :p_fpct = 1 - sum()
# end
# end

function calculate_empirical_probability_of_value(criterion, data_vector)
    emp_prob = 1 - sum(data_vector .<= criterion)/length(data_vector)
    emp_prob
end

function _025(x)
    if any(ismissing, x)
        return NaN
    else
        return quantile(collect(skipmissing(x)), 0.025)
    end
end

function _975(x)
    if any(ismissing, x)
        return NaN
    else
        return quantile(collect(skipmissing(x)), 0.975)
    end
end

function return_first_element(x::AbstractArray)
    if ismissing(x[1])
        return missing
    else
        return x[1]
    end
end

function update_column_names!(df)
    vec_names = string.(names(df))
    vec_names = replace.(vec_names, "__" => "_")
    vec_names = replace.(vec_names, "mean" => "avg")
    sym_names = Symbol.(vec_names)
    rename!(df, sym_names)
end

function update_unsummarized_column_names!(df)
    vec_names = string.(names(df))
    vec_names = replace.(vec_names, "_mean" => "")
    sym_names = Symbol.(vec_names)
    rename!(df, sym_names)
end

function update_text_column_names!(df)
    vec_names = string.(names(df))
    vec_names = replace.(vec_names, "_return_first_element" => "")
    sym_names = Symbol.(vec_names)
    rename!(df, sym_names)
end

function concatenate_actual_data_to_simulation!(df_simulation::DataFrame)
    df_actual = CSV.read(joinpath(@__DIR__, "mgmt_data.csv"))

    if size(df_simulation)[1] > size(df_actual)[1]
        difference = size(df_simulation)[1] - size(df_actual)[1]
        missing_df = missings(difference, size(df_actual)[2]) |> DataFrame
        rename!(missing_df, names(df_actual))
        padded_df = vcat(df_actual, missing_df)
        df_simulation = hcat(padded_df, df_simulation)
    elseif size(df_simulation)[1] == size(df_actual)[1]
        df_simulation = hcat(df_actual, df_simulation)
    else
        difference = size(df_actual)[1] - size(df_simulation)[1]
        df_simulation = hcat(df_actual[1:size(df_simulation)[1], :], df_simulation)
    end
    return df_simulation
end

function find_intersection_of_names(a::AbstractArray, b::AbstractArray)
    findall(in(b),a)
end

function replace_nans(a::Union{Real, Missing})
    if ismissing(a)
        return a
    elseif isnan(a)
        return missing
    else
        return a
    end
end

function replace_nans(a::Union{String, Missing})
    if ismissing(a)
        return a
    else
        return a
    end
end

function replace_nans(a::Missing)
    return missing
end

function replace_nans(a::Union{Real, Missing})
    if ismissing(a)
        return a
    elseif isnan(a)
        return missing
    else
        return a
    end
end


text:
### A Pluto.jl notebook ###
# v0.14.7

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ c7df3890-c5f0-11eb-0ef8-d5f1224128ed
begin
	using CSV
	using DataFrames
	using Plots
	using GLM
	using PlutoUI
	using Statistics
end

# ╔═╡ 45f5bb32-1451-4fa6-9d94-1f1e6078e478
md"# GDP PerCapita Income Analysis and Prediction using Linear Regression"

# ╔═╡ 5ae1d16f-5e5b-4bee-8fa4-e16291e04069
md"## Packages used:"

# ╔═╡ 57589c61-1a9e-4463-8fb4-d716d748265d
md"## Read the data using CSV and DataFrames: You can find the Dataset form [here](https://github.com/chidhvilash/GDP-per-capita-Linear-Regression/tree/main/Dataset)"

# ╔═╡ fd2070bf-e385-4c8b-aec6-5d70218407a2
data = CSV.read("D:\\Workspace\\Machine Learning\\GDP Analysis\\All country\\Data\\World_GDP_2020.csv",DataFrame)

# ╔═╡ 1dc9ec31-b8bf-45ee-af9b-bc34b5d95605
md"**NOTE:**Missing values are replaced with median value of the GDP"

# ╔═╡ 048b3f49-fc6a-4627-8456-e03d4db386be
md"**Lets analyse the data of India,USA and Canada**"

# ╔═╡ 59bf180e-e4dc-4088-b671-0f255e04b8e7
begin
	plot(data.year,data.IND,title="Indian vs USA vs Canada",yaxis="GDP",xaxis="years")
	plot!(data.year,data.CAN,colour="red")
	plot!(data.year,data.USA,colour="green")
end

# ╔═╡ 3e0bc932-ea9a-4102-b117-1aa844b3271b
md"**To make this notebook Interactive I have added more additional stuff to it like selecting countries, Input value to predict etc.... You can even do directly without creating a alternate arry as below.**"

# ╔═╡ c0548b0e-5a96-4ab5-b784-2be5a247377e
coun = ("ABW	AFG	AGO	ALB	AND	ARB	ARE	ARG	ARM	ASM	ATG	AUS	AUT	AZE	BDI	BEL	BEN	BFA	BGD	BGR	BHR	BHS	BIH	BLR	BLZ	BMU	BOL	BRA	BRB	BRN	BTN	BWA	CAF	CAN	CEB	CHE	CHI	CHL	CHN	CIV	CMR	COD	COG	COL	COM	CPV	CRI	CSS	CUB	CUW	CYM	CYP	CZE	DEU	DJI	DMA	DNK	DOM	DZA	EAP	EAR	EAS	ECA	ECS	ECU	EGY	EMU	ERI	ESP	EST	ETH	EUU	FCS	FIN	FJI	FRA	FRO	FSM	GAB	GBR	GEO	GHA	GIB	GIN	GMB	GNB	GNQ	GRC	GRD	GRL	GTM	GUM	GUY	HIC	HKG	HND	HPC	HRV	HTI	HUN	IBD	IBT	IDA	IDB	IDN	IDX	IMN	IND	INX	IRL	IRN	IRQ	ISL	ISR	ITA	JAM	JOR	JPN	KAZ	KEN	KGZ	KHM	KIR	KNA	KOR	KWT	LAC	LAO	LBN	LBR	LBY	LCA	LCN	LDC	LIC	LIE	LKA	LMC	LMY	LSO	LTE	LTU	LUX	LVA	MAC	MAF	MAR	MCO	MDA	MDG	MDV	MEA	MEX	MHL	MIC	MKD	MLI	MLT	MMR	MNA	MNE	MNG	MNP	MOZ	MRT	MUS	MWI	MYS	NAC	NAM	NCL	NER	NGA	NIC	NLD	NOR	NPL	NRU	NZL	OED	OMN	OSS	PAK	PAN	PER	PHL	PLW	PNG	POL	PRE	PRI	PRK	PRT	PRY	PSE	PSS	PST	PYF	QAT	ROU	RUS	RWA	SAS	SAU	SDN	SEN	SGP	SLB	SLE	SLV	SMR	SOM	SRB	SSA	SSD	SSF	SST	STP	SUR	SVK	SVN	SWE	SWZ	SXM	SYC	SYR	TCA	TCD	TEA	TEC	TGO	THA	TJK	TKM	TLA	TLS	TMN	TON	TSA	TSS	TTO	TUN	TUR	TUV	TZA	UGA	UKR	UMC	URY	USA	UZB	VCT	VEN	VGB	VIR	VNM	VUT	WLD	WSM	XKX	YEM	ZAF	ZMB	ZWE")

# ╔═╡ b61a5400-6667-4500-ae24-a30118c8930f
coun1 = split(coun,"\t")

# ╔═╡ 5c5d200f-5177-4850-b097-35228bb8ec3d
md"**Please select the country name for the analysis and Prediction**
$@bind country Select(coun1)
"

# ╔═╡ 4b9f6d4d-bcc1-4512-ab04-919dcbe70e15
md"**Scatter plot of selected countries GDP**"

# ╔═╡ ef6b4555-187e-44aa-9736-97842f7908ff
scatter(data[:,country],yaxis="GDP",xaxis="years",title="Scatter Plot of $country")

# ╔═╡ e22edf5a-3c8a-42b3-8eb5-9a6843e5e5f7
md"**Below table provides us the year and GDP of the selected country**"

# ╔═╡ e5c3a93f-4221-4746-9c72-aed3382c9829
data2 = (y=float(data.year),z = data[:,country])

# ╔═╡ 912c2ce5-9e9c-437c-81a9-68402a687d42
md"**Now let us replace missing values with median of the entire GDP value of that country**"

# ╔═╡ e9d842ee-2020-498b-9ccb-51ce7d417807
replace!(data2.z, missing => median(skipmissing(data2.z)))

# ╔═╡ 13fee9d5-5c4f-4aa3-a99a-15387d3c8418
md"**Now Let us find the Coefficents: Intercept and Slope**"

# ╔═╡ 79b94b7c-2d2e-474f-82c8-e421fcb7b984
reg = fit(LinearModel,@formula(z~y),data2)

# ╔═╡ e34d411c-1569-4bc1-9489-7c27cc71125b
md"**Now let us predict the Line**"

# ╔═╡ 61f49974-7a0d-4083-b714-4e0f942f22d2
linear_regression = predict(reg)

# ╔═╡ a80795ec-1d90-4aa4-a557-e92955302f56
md"**As you can see below the we have the plot of predicted line and the scatter plot of the given dataset**"

# ╔═╡ 2c13d909-78c1-4dbe-8de5-5c6ca3b95791
begin
	scatter(data2.y,data2.z,label="Scatter Plot")
	plot!(data2.y,linear_regression,label="Predicted Line")
end

# ╔═╡ 9c093e8a-142e-4f18-a7ba-2cc2fdf6c8e0
md"**Select year for future prediction**: $@bind yr NumberField(2020:2100)"

# ╔═╡ e86667e0-50a7-4ffa-a9f5-e0304148f179
md"**Predicted GDP of $country in year $yr:** "

# ╔═╡ 99a505e6-6e5b-4c13-a59c-64ae1bdc9c24
pGDP = coef(reg)[2] * yr + coef(reg)[1]

# ╔═╡ 5d1c7d64-fa5b-482a-b3fb-700ba2190c11

# ╔═╡ Cell order:
# ╟─45f5bb32-1451-4fa6-9d94-1f1e6078e478
# ╟─5ae1d16f-5e5b-4bee-8fa4-e16291e04069
# ╠═c7df3890-c5f0-11eb-0ef8-d5f1224128ed
# ╟─57589c61-1a9e-4463-8fb4-d716d748265d
# ╠═fd2070bf-e385-4c8b-aec6-5d70218407a2
# ╟─1dc9ec31-b8bf-45ee-af9b-bc34b5d95605
# ╟─048b3f49-fc6a-4627-8456-e03d4db386be
# ╠═59bf180e-e4dc-4088-b671-0f255e04b8e7
# ╟─3e0bc932-ea9a-4102-b117-1aa844b3271b
# ╟─c0548b0e-5a96-4ab5-b784-2be5a247377e
# ╟─b61a5400-6667-4500-ae24-a30118c8930f
# ╟─5c5d200f-5177-4850-b097-35228bb8ec3d
# ╟─4b9f6d4d-bcc1-4512-ab04-919dcbe70e15
# ╠═ef6b4555-187e-44aa-9736-97842f7908ff
# ╟─e22edf5a-3c8a-42b3-8eb5-9a6843e5e5f7
# ╠═e5c3a93f-4221-4746-9c72-aed3382c9829
# ╟─912c2ce5-9e9c-437c-81a9-68402a687d42
# ╠═e9d842ee-2020-498b-9ccb-51ce7d417807
# ╠═13fee9d5-5c4f-4aa3-a99a-15387d3c8418
# ╠═79b94b7c-2d2e-474f-82c8-e421fcb7b984
# ╟─e34d411c-1569-4bc1-9489-7c27cc71125b
# ╠═61f49974-7a0d-4083-b714-4e0f942f22d2
# ╟─a80795ec-1d90-4aa4-a557-e92955302f56
# ╠═2c13d909-78c1-4dbe-8de5-5c6ca3b95791
# ╟─9c093e8a-142e-4f18-a7ba-2cc2fdf6c8e0
# ╟─e86667e0-50a7-4ffa-a9f5-e0304148f179
# ╠═99a505e6-6e5b-4c13-a59c-64ae1bdc9c24
# ╟─5d1c7d64-fa5b-482a-b3fb-700ba2190c11

max_stars_repo_name: chidhvilash/GDP-PerCapita-Analysis

max_issues_repo_name: chidhvilash/GDP-PerCapita-Analysis

max_forks_repo_name: chidhvilash/GDP-PerCapita-Analysis


text:
#%%
using Glimpse
const Gl = Glimpse
#[%%
Entity(dio, Gl.GuiText("test"))
Entity(dio, Gl.GuiText("test1"))
#%%


text:
module TestAPI

using Test
using TestFunctionRunner

using ..Samples: SampleSimplePass

function test_run_modules()
    TestFunctionRunner.run([SampleSimplePass, SampleSimplePass])
end

end  # module


text:
using RobustAndOptimalControl, ControlSystems

G = tf(200, [10, 1])*tf(1, [0.05, 1])^2     |> ss
Gd = tf(100, [10, 1])                       |> ss
W1 = tf([1, 2], [1, 1e-6])                  |> ss
Gs = G*W1
Ks, γ, info = glover_mcfarlane(Gs, 1.1)
@test info.γmin ≈ 2.34 atol=0.005

if isinteractive()
    bodeplot([G, Gs, Gs*Ks]) |> display

    plot( step(Gd*feedback(1, G*W1), 3))
    plot!(step(Gd*feedback(1, G*W1*Ks), 3)) |> display

    nyquistplot([G*W1, G*W1*Ks], ylims=(-2,1), xlims=(-2, 1), Ms_circles=1.5) |> display
end

## Discrete case
disc = G->c2d(G, 0.01)
G = tf(200, [10, 1])*tf(1, [0.05, 1])^2     |> ss |> disc
Gd = tf(100, [10, 1])                       |> ss |> disc
W1 = tf([1, 2], [1, 1e-6])                  |> ss |> disc
K, γ, info = glover_mcfarlane(G, 1.1; W1)
@test info.γmin ≈ 2.4086 atol=0.005

Gcl0 = let
    tempA = [0.6542 -0.2067 -1.6744 0.2619 -0.3391 -1.0397 -6.283 -0.6716; 0.1309 0.9823 -0.1434 0.0224 -0.029 -0.089 -0.538 -0.0575; 0.0014 0.0199 0.999 0.0002 -0.0002 -0.0006 -0.0037 -0.0004; 0.0 0.0 -0.1271 1.0 -0.0259 -0.0794 -0.4799 -0.0513; 0.0 0.0 -1.8079 0.0 0.315 -1.2464 -6.1494 -0.4097; 0.0 0.0 0.3374 0.0 0.1019 0.8933 -1.0187 -0.0351; 0.0 0.0 0.1454 0.0 0.0012 0.0193 0.8499 -0.0002; 0.0 0.0 0.0181 0.0 -0.0259 -0.0794 -0.6251 0.9487]
    tempB = [-0.1065 0.1309; -0.0091 0.0112; -0.0001 0.0001; -0.0081 0.0; -0.1157 0.0; 0.0216 0.0; 0.0093 0.0; 0.0012 0.0]
    tempC = [0.0 0.0 15.625 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 -12.7105 2.0 -2.5903 -7.941 -47.9883 -5.1294]
    tempD = [1.0 0.0; -0.8135 0.0]
    ss(tempA, tempB, tempC, tempD, 0.01)
end

@test hinfnorm2(Gcl0 - info.Gcl)[1] < 0.02

if isinteractive()
    bodeplot([G, G*K]) |> display
    plot( step(Gd*feedback(1, info.Gs), 3))
    plot!(step(Gd*feedback(1, G*K), 3)) |> display
    nyquistplot([info.Gs, G*K], ylims=(-2,1), xlims=(-2, 1), Ms_circles=1.5) |> display
end

Ko = observer_controller(info)
@test nugap(G*K, info.Gs*Ko)[1] < 1e-6 # The loop transfer is used since the observer controller is for the scaled plant whereas K is for the unscaled plant.
# bodeplot([G*K, info.Gs*Ko])

## Test the observer predictor
Gs = info.Gs
Ph = observer_predictor(info)
u = randn(Gs.nu, 100)
res = lsim(Gs, u)
uy = [u; res.y]
reso = lsim(Ph, uy)

# plot(res)
# plot!(reso) |> display
# plot(res.x')
# plot!(reso.x') |> display

@test res.y ≈ reso.y
@test res.x ≈ reso.x

## Strictly proper case
Ksp, γsp, infosp = glover_mcfarlane(G, 1.1; W1, strictly_proper=true)
@test infosp.γmin > 2.4086
@test iszero(Ksp.D)

Ko = observer_controller(infosp)
@test nugap(G*Ksp, infosp.Gs*Ko)[1] < 1e-6
# bodeplot([G*Ksp, infosp.Gs*Ko])

# Can be manually verified to be "slightly worse" than the controller without strictly proper restriction.
# if isinteractive()
#     plot( step(Gd*feedback(1, info.Gs), 3))
#     plot!(step(Gd*feedback(1, G*K), 3))
#     plot!(step(Gd*feedback(1, G*Ksp), 3)) |> display
#     nyquistplot([info.Gs, G*K, G*Ksp], ylims=(-2,1), xlims=(-2, 1), Ms_circles=1.5) |> display
# end

Gs = info.Gs
Ph = observer_predictor(infosp)
u = randn(Gs.nu, 100)
res = lsim(Gs, u)
uy = [u; res.y]
reso = lsim(Ph, uy)
@test res.y ≈ reso.y
@test res.x ≈ reso.x

##
W1h = hanus(W1)
@test W1h.nu == 2W1.nu

##

P = tf([1, 5], [1, 2, 10])
W1 = tf(1,[1, 0])

K,γ,info = glover_mcfarlane(ss(P); W1)

Gcl0 = let
    tempA = [-2.0 -2.5 2.0 0.0 0.0 0.0; 4.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 -0.7189 0.0453 -2.9175; 0.0493 0.0616 0.0 -2.0493 -2.5616 2.0; 0.399 0.4988 0.0 3.601 -0.4988 0.0; 0.5 0.625 0.0 -1.2189 -0.5797 -2.9175]
    tempB = [0.0 2.0; 0.0 0.0; 0.0 0.0; 0.0985 0.0; 0.7981 0.0; 1.0 0.0]
    tempC = [0.5 0.625 0.0 0.0 0.0 0.0; 0.0 0.0 1.0 0.0 0.0 0.0]
    tempD = [1.0 0.0; 0.0 0.0]
    ss(tempA, tempB, tempC, tempD)
end

K0 = let
    tempA = [0.0 -0.7189 0.0453 -2.9175; 0.0 -2.0493 -2.5616 2.0; 0.0 3.601 -0.4988 0.0; 0.0 -1.2189 -0.5797 -2.9175]
    tempB = [0.0; 0.0985; 0.7981; 1.0;;]
    tempC = [1.0 0.0 0.0 0.0]
    tempD = [0.0;;]
    ss(tempA, tempB, tempC, tempD)
end

Gcl = extended_gangoffour(P, K)
@test hinfnorm2(Gcl-Gcl0)[1] < 1e-4
@test hinfnorm2(info.Gcl-Gcl0)[1] < 1e-4
@test nugap(K, -K0)[1] < 1e-4
@test info.margin ≈ 0.6325 atol=1e-3

@test ncfmargin(P, K)[1] ≈ 0.4472 atol=1e-3

p = ss(tf(4, [1, -0.001]))
cL = 1
cH = 10
@test ncfmargin(p,cL)[1] ≈ 0.7069 atol=1e-3

P = ss(tf(1, [1,1]))
K = ss(1)
S, PS, CS, T = RobustAndOptimalControl.gangoffour2(P,K)
# gangoffourplot(P, K)
# bodeplot!(extended_gangoffour(P, K), plotphase=false)

@test all(isstable.((S, PS, CS, T)))
gof = extended_gangoffour(P, K)
isstable(gof)
@test nugap(S, gof[1,1])[1] < 1e-6
@test nugap(PS, gof[1,2])[1] < 1e-6
@test nugap(CS, -gof[2,1])[1] < 1e-6 # NOTE: slightly disturbing to have - here
@test nugap(T, -gof[2,2])[1] < 1e-6

## 2 DOF GMcF

P = tf([1, 5], [1, 2, 10])
W1 = tf(1,[1, 0]) |> ss

Tref = tf(1, [1, 1]) |> ss

K1dof, γ1, info1 = glover_mcfarlane(ss(P), 1.1; W1)
K2dof, γ2, info2 = RobustAndOptimalControl.glover_mcfarlane_2dof(ss(P), Tref, 1.1, 1.1; W1)

G1 = feedback(P*K1dof)
G2 = info2.Gcl

plot([step(G1, 15), step(G2, 15), step(Tref, 15)], lab=["1-DOF" "2-DOF" "Tref"])
w = exp10.(LinRange(-3, 2, 200))
bodeplot(info2.K1, w, lab="Feedforward filter")

@test dcgain(G2)[] ≈ 1 rtol=1e-4