text:
@testset "custom deterministic generative function with custom update and gradient" begin

    struct MyDeterministicGFState
        sum::Float64
        prev::AbstractArray
        my_param::Float64
    end

    # sum(arr) * my_param
    # arr is the argument
    # my_param is a trainable parameter

    mutable struct MyDeterministicGF <: CustomDetermGF{Float64,MyDeterministicGFState}
        my_param::Float64
        my_grad::Float64
    end

    MyDeterministicGF() = MyDeterministicGF(1., 0.)

    function Gen.apply_with_state(gen_fn::MyDeterministicGF, args::Tuple)
        arr = args[1]
        retval = sum(arr) * gen_fn.my_param
        state = MyDeterministicGFState(retval, arr, gen_fn.my_param)
        retval, state
    end

    function Gen.update_with_state(gen_fn::MyDeterministicGF, state, args, argdiffs)
        arr = args[1]
        retval = sum(arr) * state.my_param
        state = MyDeterministicGFState(retval, arr, state.my_param)
        state, retval, UnknownChange()
    end

    function Gen.update_with_state(::MyDeterministicGF, state, args, argdiffs::Tuple{NoChange})
        state, state.sum * state.my_param, NoChange()
    end

    function Gen.update_with_state(::MyDeterministicGF, state, args, argdiffs::Tuple{VectorDiff})
        arr = args[1]
        retval = state.sum * state.my_param
        for i in keys(argdiffs[1].updated)
            retval += (arr[i] - state.prev[i])
        end
        prev_length = length(state.prev)
        new_length = length(arr)
        for i=prev_length+1:new_length
            retval += arr[i]
        end
        for i=new_length+1:prev_length
            retval -= arr[i]
        end
        state = MyDeterministicGFState(retval, arr, state.my_param)
        state, retval, UnknownChange()
    end

    Gen.has_argument_grads(::MyDeterministicGF) = (true,)

    Gen.accepts_output_grad(::MyDeterministicGF) = true

    function Gen.gradient_with_state(::MyDeterministicGF, state, args, retgrad)
        arr_gradient = fill(retgrad * state.my_param, length(state.prev))
        (arr_gradient,)
    end

    function Gen.accumulate_param_gradients_determ!(
            gen_fn::MyDeterministicGF, state, args, retgrad, scaler)
        gen_fn.my_grad += (retgrad * state.sum) * scaler
        arr_gradient = fill(retgrad * state.my_param, length(state.prev))
        (arr_gradient,)
    end

    # simulate
    gen_fn = MyDeterministicGF()
    trace = simulate(gen_fn, ([1, 2, 3],))
    @test get_retval(trace) == 1 + 2 + 3
    @test get_score(trace) == 0.
    @test project(trace, EmptySelection()) == 0.
    @test isempty(get_choices(trace))
    @test get_args(trace) == ([1, 2, 3],)
    @test get_gen_fn(trace) == gen_fn

    # generate
    trace, w = generate(MyDeterministicGF(), ([1, 2, 3],))
    @test w == 0.
    @test get_retval(trace) == 1 + 2 + 3

    # update (UnknownChange)
    trace = simulate(MyDeterministicGF(), ([1, 2, 3],))
    new_trace, w, retdiff = update(trace, ([1, 2, 4],), (UnknownChange(),), EmptyChoiceMap())
    @test w == 0.
    @test get_retval(new_trace) == 1 + 2 + 4
    @test get_args(new_trace) == ([1, 2, 4],)
    @test retdiff == UnknownChange()

    # update (VectorDiff)
    trace = simulate(MyDeterministicGF(), ([1, 2, 3],))
    diff = VectorDiff(4, 3, Dict(3 => UnknownChange()))
    new_trace, w, retdiff = update(trace, ([1, 2, 4, 5],), (diff,), EmptyChoiceMap())
    @test w == 0.
    @test get_retval(new_trace) == 1 + 2 + 4 + 5
    @test get_args(new_trace) == ([1, 2, 4, 5],)
    @test retdiff == UnknownChange()

    # regenerate
    trace = simulate(MyDeterministicGF(), ([1, 2, 3],))
    new_trace, w, retdiff = regenerate(trace, ([1, 2, 4],), (UnknownChange(),), EmptySelection())
    @test w == 0.
    @test get_retval(new_trace) == 1 + 2 + 4
    @test get_args(new_trace) == ([1, 2, 4],)
    @test retdiff == UnknownChange()

    # choice gradients
    trace = simulate(MyDeterministicGF(), ([1, 2, 3],))
    arg_grads, _, _ = choice_gradients(trace, EmptySelection(), 2.)
    @test arg_grads[1] == [2., 2., 2.]

    # accumulate parameter gradients
    gen_fn = MyDeterministicGF()
    trace = simulate(gen_fn, ([1, 2, 3],))
    arg_grads = accumulate_param_gradients!(trace, 2., 3.)
    @test arg_grads[1] == [2., 2., 2.] # not scaled by 3.
    @test gen_fn.my_grad == 2. * (1 + 2 + 3) * 3. # scaled by 3.
end

@testset "custom gradient GF" begin

    struct CustomPlus <: CustomGradientGF{Float64} end
    Gen.has_argument_grads(::CustomPlus) = (true, true)
    Gen.apply(::CustomPlus, args::Tuple{Float64,Float64}) = args[1] + args[2]
    Gen.gradient(::CustomPlus, args, retval, retgrad) = (retgrad, retgrad)
    custom_plus = CustomPlus()
    trace = simulate(custom_plus, (1., 2.))

    # choice gradients
    arg_grads, _, _ = choice_gradients(trace, EmptySelection(), 2.)
    @test arg_grads == (2., 2.)

    # accumulate parameter gradients
    arg_grads = accumulate_param_gradients!(trace, 2., 3.)
    @test arg_grads == (2., 2.) # note arg grads are not scaled
end

@testset "custom update GF" begin

    mutable struct MyState
        prev_arr::Vector{Float64}
        sum::Float64
    end

    struct MySum <: CustomUpdateGF{Float64,MyState} end
    mysum = MySum()

    function Gen.apply_with_state(::MySum, args)
        arr = args[1]
        s = sum(arr)
        state = MyState(arr, s)
        (s, state)
    end

    function Gen.update_with_state(::MySum, state, args, argdiffs::Tuple{VectorDiff})
        arr = args[1]
        prev_sum = state.sum
        retval = prev_sum
        for i in keys(argdiffs[1].updated)
            retval += (arr[i] - state.prev_arr[i])
        end
        prev_length = length(state.prev_arr)
        new_length = length(arr)
        for i=prev_length+1:new_length
            retval += arr[i]
        end
        for i=new_length+1:prev_length
            retval -= arr[i]
        end
        state = MyState(arr, retval)
        (state, retval, UnknownChange())
    end

    Gen.num_args(::MySum) = 1

    # simulate
    trace = simulate(mysum, ([1, 2, 3],))
    @test get_retval(trace) == 1 + 2 + 3
    @test get_score(trace) == 0.

    # update (UnknownChange)
    new_trace, w, retdiff, discard = update(trace, ([1, 2, 4],), (UnknownChange(),), EmptyChoiceMap())
    @test w == 0.
    @test get_retval(new_trace) == 1 + 2 + 4
    @test get_args(new_trace) == ([1, 2, 4],)
    @test retdiff == UnknownChange()
    @test isempty(discard)

    # update (VectorDiff)
    diff = VectorDiff(4, 3, Dict(3 => UnknownChange()))
    new_trace, w, retdiff, discard = update(trace, ([1, 2, 4, 5],), (diff,), EmptyChoiceMap())
    @test w == 0.
    @test get_retval(new_trace) == 1 + 2 + 4 + 5
    @test get_args(new_trace) == ([1, 2, 4, 5],)
    @test retdiff == UnknownChange()
    @test isempty(discard)
end

max_stars_repo_name: georgematheos/Gen

max_issues_repo_name: georgematheos/Gen

max_forks_repo_name: georgematheos/Gen


text:
"""
    disparity(M, ft; refGrp=nothing, func=/)

Computes disparity for fairness tensor `ft` with respect to an array of metrics `M`.

For any class A and a reference Group B, `disparity = func(metric(A), metric(B))`.
By default `func` is `/` .

A dataframe is returned with disparity values for all combinations of metrics and classes.
It contains a column named labels for the classes and has a column for disparity of each metric in M.
The column names are metric names appended with `_disparity`.

## Keywords

* `refGrp=nothing` : The reference group
* `func=/` : The function used to evaluate disparity. This function should take 2 arguments.
The second argument shall correspond to reference group.

Please note that division by 0 will result in NaN
"""
function disparity(measures::Vector, ft::FairTensor{C}; refGrp=nothing, func=/) where C
    refGrp!==nothing || throw(ArgumentError("Value of reference group needs to be provided"))
    refGrpIdx = _ftIdx(ft, refGrp)
    df = DataFrame(labels=ft.labels)
    for measure in measures
        colName = MMI.name(measure) * "_disparity"
        colDisparity = Symbol(colName)
        # TODO : _ftIdx is repeatedly called internally. Instead the value can be stored and reused.
        baseVal = measure(ft, grp=refGrp)
        arr = zeros(Float64, C)
        for i in 1:C
            arr[i] = func(measure(ft, grp=ft.labels[i]), baseVal)
        end
        df[:, colDisparity] = arr
    end
    return df
end

"""
    Disparity

Disparity uses the `disparity` function and has information about the protected attribute, reference Group and custom function
It will help in automatic evaluation using MLJ.evaluate
"""
struct Disparity <: MLJBase.Measure
    measure::Measure
    grp::Symbol
    refGrp
    func
end

"""
    Disparity(measure, grp=:class, refGrp=nothing, func=/)

Instantiates the struct Disparity.
"""
function Disparity(measure::Measure; grp::Symbol=:class, refGrp=nothing, func=/)
    Disparity(measure, grp, refGrp, func)
end

"""
    Disparities(measures, grp=:class, refGrp=nothing, func=/)

Creates instances of Disparity struct for each measure in measures.
"""
function Disparities(measures::Vector{<:Measure}; grp::Symbol, refGrp=nothing, func=/)
    wrappedMetrics = []
    for measure in measures
        push!(wrappedMetrics, Disparity(measure; grp=grp, refGrp=refGrp, func=func))
    end
    return wrappedMetrics
end

function (D::Disparity)(ŷ, X, y)
    if typeof(ŷ) <: UnivariateFiniteArray
        ŷ = StatsBase.mode.(ŷ)
    end
    if typeof(y) <: UnivariateFiniteArray
        y = StatsBase.mode.(y)
    end
    grps = X[:, D.grp]
    n_grps = length(levels(grps))

    ft = fair_tensor(categorical(ŷ), categorical(y), categorical(grps))

    valuesDict = Dict()
    df = disparity([D.measure], ft,refGrp=D.refGrp, func=D.func)
    vals = df[!, names(df)[2]]
    for i in 1:n_grps
        valuesDict[df.labels[i]]=vals[i]
    end

    valuesDict["overall"] = D.func(D.measure(ft), D.measure(ft, grp=D.refGrp))
    return valuesDict
end

MLJBase.name(::Type{<:Disparity}) = "Disparity"
# MLJBase.target_scitype(::Type{<:Disparity}) = AbstractArray{Multiclass{2},1}
MLJBase.supports_weights(::Type{<:Disparity}) = false # for now
# MLJBase.prediction_type(::Type{<:Disparity}) = :deterministic # Not specifying it to have check_measures false
MLJBase.orientation(::Type{<:Disparity}) = :other # other options are :score, :loss
MLJBase.reports_each_observation(::Type{<:Disparity}) = false
MLJBase.aggregation(::Type{<:Disparity}) = Mean()
MLJBase.is_feature_dependent(::Type{<:Disparity}) = true

# To display the name of actual measure while printing instead of the name Disparity
Base.show(stream::IO, D::Disparity) = print(stream, string(D.measure)*"_disparity")
#------------------------------------------------------------------
# Helper for parity function. Default criteria to calculate parity
_parity_func(ϵ) = x -> isnan(x) ? NaN : (1-ϵ) <= x <= 1/(1-ϵ)

"""
    parity(df, ϵ=nothing; func=nothing)

Takes the dataframe `df` returned from `disparity` function and adds columns for parity values corresponding to each disparity column.
It then calculates the parity values for measures that were passed for disparity calculation.

Parity is a boolean value indicating whether a fairness constraint is satisfied by disparity values.

The default fairness criteria for a disparity value `x` and fairness threshold `ϵ` for a group is:

    (1-ϵ) <= x <= 1/(1-ϵ)

Here ϵ is the fairness threshold which is required if default fairness constraint is used.

## Keywords

* `func=nothing` : single argument function specifying custom Fairness constraint for disparity instead of the default criteria
"""
function parity(df::DataFrame, ϵ=nothing; func=nothing)
    if func==nothing
        ϵ!=nothing || throw(ArgumentError("For using default fairness criteria, ϵ i.e. fairness threshold should be passed."))
        func = _parity_func(ϵ)
    end
    cols = names(df)
    M = string.(cols)
    for i in 1:length(cols)
        if !endswith(M[i], "_disparity")
            continue
        end
        m = replace(M[i], "disparity" => "parity")
        m = Symbol(m)
        df[:, m] = func.(df[:, cols[i]])
    end
    return df
end

ext: jl


text:
using DataArrays
dv = @data([7, 3, NA, 5, 42])
    #> 5-element DataArray{Int64,1}: 7  3   NA  5 42
# array(dv) #> ERROR: NAException
# convert(Array, dv) #> ERROR: NAException("Cannot convert DataArray with NA's to desired type")
mean(dv) #> NA
sum(dv) #> NA
dv[5] = NA  #> dv now becomes [7, 3, NA, 5, NA]

# getting rid of NA's:
sum(dropna(dv)) #> 15
repl = -1
sum(array(dv, repl)) #> 13

max_stars_repo_path: Module 1/Chapter08/dataarrays.jl

max_issues_repo_path: Module 1/Chapter08/dataarrays.jl

max_forks_repo_path: Module 1/Chapter08/dataarrays.jl


text:
using Documenter, NarmaxLasso, Plots
println("Generating plots")

plotsDir = (pwd()[end-3:end] == "docs") ? "src/plots" : "docs/src/plots"
if ~isdir(plotsDir)
    mkdir(plotsDir)
end
Plots.pyplot()

## ARMAX example
u = generate_random_input(2000, 5);
armax_mdl(y, u, v) = 0.5*y[1] - 0.5*u[1] + 0.5*v[1];
yterms0 = [1]; uterms0 = [1]; vterms0 = [1];
y = simulate(u, armax_mdl, yterms0, uterms0, vterms0;
             σv=0.3);
ny = 1; nu = 1; nv = 1; order = 1;
mdl = generate_all(NarmaxRegressors, Monomial, ny, nu, nv, order);
result = narmax_lasso(y, u, mdl);
Plots.plot(result)
Plots.savefig(plotsDir*"/ARMAXPaths.png")

## NARMAX example
u = generate_random_input(2000, 5);
narmax_mdl(y, u, v) = 0.1*y[1]^2 - 0.5*u[1] + 0.5*v[1];
yterms0 = [1, 2]; uterms0 = [1]; vterms0 = [1];
y = simulate(u, narmax_mdl, yterms0, uterms0, vterms0;
             σv=0.3);
ny = 2; nu = 1; nv = 1; order = 1:2;
mdl = generate_all(NarmaxRegressors, Monomial, ny, nu, nv, order);
result = narmax_lasso(y, u, mdl);
Plots.plot(result, xscale=:log10)
Plots.savefig(plotsDir*"/NARMAXPaths.png")


text: import Base: &, |, -, !


text:
@testset "GridWorld tests" begin

    next_states = [
        1 5 2 15
        2 6 3 1
        3 7 3 2
        15 8 5 4
        1 9 6 4
        2 10 7 5
        3 11 7 6
        4 12 9 8
        5 13 10 8
        6 14 11 9
        7 15 11 10
        8 12 13 12
        9 13 14 12
        10 14 15 13
        15 15 15 15
    ]

    rewards = -1 * ones(size(next_states))
    rewards[15, :] .= 0

    fmdp = MDP.DeterministicFiniteMDP((s, a) -> next_states[s, a], (s, a) -> rewards[s, a], 4, 15, 15)

    optimal_V = [-1.0, -2.0, -3.0, -1.0, -2.0, -3.0, -2.0, -2.0, -3.0, -2.0, -1.0, -3.0, -2.0, -1.0, 0.0]

    uniform_random_policy = 0.25 * ones(size(next_states)...)
    V_uniform_random_policy = [-14.0, -20.0, -22.0, -14.0, -18.0, -20.0, -20.0, -20.0, -20.0, -18.0, -14.0, -22.0, -20.0, -14.0, 0.0]

    random_deterministic_policy = rand(1:4, size(next_states, 1))

    @testset "Dynamic Programming (DP) tests" begin

        dp_V_tol = 1e-4
        dp_maxiter_eval = 1000
        dp_maxiter_optim = 1000

        @testset "DP: policy evaluation (textbook)" begin
            𝐏 = copy(uniform_random_policy)

            V = MDP.allocate_V(fmdp)
            Q = MDP.allocate_Q(fmdp)

            V[MDP.terminal_state(fmdp)] = 0
            iters_no, converged = MDP.dp_evaluate_policy_textbook!(V, Q, fmdp, 𝐏, 1.0; tol = dp_V_tol, maxiter = dp_maxiter_eval)

            ΔV = abs.((V - V_uniform_random_policy))
            err = max(ΔV...)
            @info "DP: policy evaluation (textbook): max abs err = $err (iters no = $iters_no)"
            @test converged && isapprox(err, 0.0; atol = 1e-1)
        end

        @testset "DP: policy evaluation via V update" begin
            𝐏 = copy(uniform_random_policy)

            V = MDP.allocate_V(fmdp)
            V[MDP.terminal_state(fmdp)] = 0

            converged = false
            iters_no = dp_maxiter_eval
            for i = 1:dp_maxiter_eval
                Δ = MDP.dp_update_V!(V, fmdp, 𝐏, 1.0)
                if Δ < dp_V_tol
                    converged = true
                    iters_no = i
                    break
                end
            end

            ΔV = abs.((V - V_uniform_random_policy))
            err = max(ΔV...)
            @info "DP: policy evaluation (V update): max abs err = $err (iters no = $iters_no)"
            @test converged && isapprox(err, 0.0; atol = 1e-1)
        end

        @testset "DP: policy optimization (textbook GPI)" begin
            𝐩 = copy(random_deterministic_policy)

            V = MDP.allocate_V(fmdp)
            Q = MDP.allocate_Q(fmdp)

            V[MDP.terminal_state(fmdp)] = 0
            converged = false
            iters_no = dp_maxiter_optim
            for i = 1:dp_maxiter_optim
                MDP.dp_evaluate_policy_textbook!(V, Q, fmdp, 𝐩, 1.0; tol = dp_V_tol, maxiter = 5)
                modified = MDP.𝐩_from_Q!(𝐩, Q)
                if !modified
                    iters_no = i
                    converged = true
                    break
                end
            end # for: iterations

            ΔV = abs.((V - optimal_V))
            err = max(ΔV...)
            @info "DP: policy optimization (textbook GPI): max abs err = $err (iters no = $iters_no)"
            @test converged && isapprox(err, 0.0; atol = 1e-1)
        end

        @testset "DP: policy optimization (value iteration)" begin
            V = MDP.allocate_V(fmdp)
            V[MDP.terminal_state(fmdp)] = 0

            converged = false
            iters_no = dp_maxiter_optim
            for i = 1:dp_maxiter_optim
                Δ = MDP.dp_update_V!(V, fmdp, 1.0)
                if Δ < dp_V_tol
                    converged = true
                    iters_no = i
                    break
                end
            end

            ΔV = abs.((V - optimal_V))
            err = max(ΔV...)
            @info "DP: policy optimization (value iteration): max abs err = $err (iters no = $iters_no)"
            @test converged && isapprox(err, 0.0; atol = 1e-1)
        end

    end

    @testset "Monte Karlo (MK) tests" begin

        @testset "MK: policy evaluation" begin
            𝐏 = copy(uniform_random_policy)
            simulator = MDP.create_simulator_from_policy(fmdp, 𝐏, 1000)

            V = MDP.allocate_V(fmdp)
            Q = MDP.allocate_Q(fmdp)

            MDP.mk_evaluate_policy!(Q, simulator, 1.0; maxiter = 10000)

            MDP.V_from_Q!(V, Q, 𝐏)
            ΔV = abs.((V - V_uniform_random_policy))
            err = max(ΔV...)
            @info "MK: policy evaluation: max abs err = $err"
            @test isapprox(err, 0.0; atol = 1)
        end

        @testset "MK: policy evaluation (incremental)" begin
            𝐏 = copy(uniform_random_policy)
            simulator = MDP.create_simulator_from_policy(fmdp, 𝐏, 1000)

            V = MDP.allocate_V(fmdp)
            Q = MDP.allocate_Q(fmdp)

            ΔQ = -Inf
            for _ = 1:10000
                ΔQ = MDP.mk_update_Q!(Q, 0.99, simulator, 1.0)
            end

            MDP.V_from_Q!(V, Q, 𝐏)
            ΔV = abs.((V - V_uniform_random_policy))
            err = max(ΔV...)
            @info "MK: policy evaluation (incremental): max abs err = $err (final ΔQ = $ΔQ)"
            @test isapprox(err, 0.0; atol = 5)
        end

        @testset "MK: ε-greedy simulation" begin
            Q = MDP.allocate_Q(fmdp)
            MDP.Q_from_V!(Q, optimal_V, fmdp, 1.0)

            𝐩 = copy(random_deterministic_policy)
            MDP.𝐩_from_Q!(𝐩, Q)

            simulator = MDP.create_simulator_from_policy(fmdp, 𝐩, 0.05, 100)
            episode = simulator(5, 2)
        end

        @testset "MK: policy optimization" begin
            𝐩 = copy(random_deterministic_policy)
            simulator = MDP.create_simulator_from_policy(fmdp, 𝐩, 0.05, 100)

            Q = MDP.allocate_Q(fmdp)
            MDP.𝐩_from_Q!(𝐩, Q)

            for i = 1:10
                MDP.mk_evaluate_policy!(Q, simulator, 1.0; maxiter = 10000)
                modified = MDP.𝐩_from_Q!(𝐩, Q)
            end # for: iterations

            V = MDP.allocate_V(fmdp)
            MDP.V_from_Q!(V, Q, 𝐩)
            ΔV = abs.(V - optimal_V)
            err = max(ΔV...)
            @info "MK: policy optimization: max abs err = $err"
            @test isapprox(err, 0.0; atol = 1)
        end

        @testset "MK: policy optimization (incremental)" begin
            Q = MDP.allocate_Q(fmdp)
            simulator = MDP.create_simulator_from_Q(fmdp, Q, 0.05, 100)

            ΔQ = -Inf
            for i = 1:10000
                ΔQ = MDP.mk_update_Q!(Q, 0.95, simulator, 1.0)
            end # for: iterations

            V = MDP.allocate_V(fmdp)
            MDP.V_from_Q!(V, Q)
            ΔV = abs.(V - optimal_V)
            err = max(ΔV...)
            @info "MK: policy optimization (incremental): max abs err = $err (final ΔQ = $ΔQ)"
            @test isapprox(err, 0.0; atol = 1)
        end

        @testset "TD(0): policy evaluation" begin
            𝐏 = copy(uniform_random_policy)

            V = MDP.allocate_V(fmdp)
            V[MDP.terminal_state(fmdp)] = 0

            δV = -Inf
            for i = 1:5000
                δV = -Inf
                s = rand(1:MDP.states_no(fmdp)-1)
                if s == MDP.terminal_state(fmdp)
                    continue
                end
                while true
                    a = MDP.draw_action_from_policy(s, 𝐏)
                    sn, r = MDP.act_once(fmdp, s, a)
                    δ = MDP.td_0_update_V!(V, 0.01, s, r, sn, 1.0)
                    δV = max(δV, abs(δ))
                    if sn == MDP.terminal_state(fmdp)
                        break
                    end
                    s = sn
                end
                V[MDP.terminal_state(fmdp)] = 0
            end

            ΔV = abs.(V - V_uniform_random_policy)
            err = max(ΔV...)
            @info "TD(0): policy evaluation: max abs err = $err (final ΔQ = $δV)"
            @test isapprox(err, 0.0; atol = 3)
        end

        @testset "TD SARSA: policy optimization" begin
            𝐏 = copy(uniform_random_policy)

            Q = MDP.allocate_Q(fmdp)

            δQ = -Inf
            for i = 1:5000
                δQ = -Inf
                s = rand(1:MDP.states_no(fmdp)-1)
                if s == MDP.terminal_state(fmdp)
                    continue
                end
                a = rand(1:MDP.actions_no(fmdp))
                while true
                    sn, r = MDP.act_once(fmdp, s, a)
                    an = MDP.draw_action_from_Q(s, Q, 0.05)
                    δ = MDP.td_sarsa_update_Q!(Q, 0.1, s, a, r, sn, an, 1.0)
                    δQ = max(δQ, abs(δ))
                    if sn == MDP.terminal_state(fmdp)
                        break
                    end
                    s = sn
                    a = an
                end
            end

            V = MDP.allocate_V(fmdp)
            MDP.V_from_Q!(V, Q)
            ΔV = abs.(V - optimal_V)
            err = max(ΔV...)
            @info "TD SARSA: policy optimization: max abs err = $err (final ΔQ = $δQ)"
            @test isapprox(err, 0.0; atol = 3)
        end

        @testset "TD Q LEARNING: policy optimization" begin
            𝐏 = copy(uniform_random_policy)

            Q = MDP.allocate_Q(fmdp)

            δQ = -Inf
            for i = 1:5000
                δQ = -Inf
                s = rand(1:MDP.states_no(fmdp)-1)
                if s == MDP.terminal_state(fmdp)
                    continue
                end
                a = rand(1:MDP.actions_no(fmdp))
                while true
                    sn, r = MDP.act_once(fmdp, s, a)
                    δ = MDP.td_ql_update_Q!(Q, 0.1, s, a, r, sn, 1.0)
                    δQ = max(δQ, abs(δ))
                    if sn == MDP.terminal_state(fmdp)
                        break
                    end
                    s = sn
                    a = MDP.draw_action_from_Q(s, Q, 0.05)
                end
            end

            V = MDP.allocate_V(fmdp)
            MDP.V_from_Q!(V, Q)
            ΔV = abs.(V - optimal_V)
            err = max(ΔV...)
            @info "TD Q LEARNING: policy optimization: max abs err = $err (final ΔQ = $δQ)"
            @test isapprox(err, 0.0; atol = 3)
        end

    end

end


text:
@testset "DataLoader" begin
    X2 = reshape([1:10;], (2, 5))
    Y2 = [1:5;]

    d = DataLoader(X2, batchsize=2)
    @test_broken @inferred(first(d)) isa Array
    batches = collect(d)
    @test_broken  eltype(d) == typeof(X2)
    @test eltype(batches) == typeof(X2)
    @test length(batches) == 3
    @test batches[1] == X2[:,1:2]
    @test batches[2] == X2[:,3:4]
    @test batches[3] == X2[:,5:5]

    d = DataLoader(X2, batchsize=2, partial=false)
    # @inferred first(d)
    batches = collect(d)
    @test_broken eltype(d) == typeof(X2)
    @test length(batches) == 2
    @test batches[1] == X2[:,1:2]
    @test batches[2] == X2[:,3:4]

    d = DataLoader((X2,), batchsize=2, partial=false)
    # @inferred first(d)
    batches = collect(d)
    @test_broken eltype(d) == Tuple{typeof(X2)}
    @test eltype(batches) == Tuple{typeof(X2)}
    @test length(batches) == 2
    @test batches[1] == (X2[:,1:2],)
    @test batches[2] == (X2[:,3:4],)

    d = DataLoader((X2, Y2), batchsize=2)
    # @inferred first(d)
    batches = collect(d)
    @test_broken eltype(d) == Tuple{typeof(X2), typeof(Y2)}
    @test eltype(batches) == Tuple{typeof(X2), typeof(Y2)}
    @test length(batches) == 3
    @test length(batches[1]) == 2
    @test length(batches[2]) == 2
    @test length(batches[3]) == 2
    @test batches[1][1] == X2[:,1:2]
    @test batches[1][2] == Y2[1:2]
    @test batches[2][1] == X2[:,3:4]
    @test batches[2][2] == Y2[3:4]
    @test batches[3][1] == X2[:,5:5]
    @test batches[3][2] == Y2[5:5]

    # test with NamedTuple
    d = DataLoader((x=X2, y=Y2), batchsize=2)
    # @inferred first(d)
    batches = collect(d)
    @test_broken eltype(d) == NamedTuple{(:x, :y), Tuple{typeof(X2), typeof(Y2)}}
    @test eltype(batches) == NamedTuple{(:x, :y), Tuple{typeof(X2), typeof(Y2)}}
    @test length(batches) == 3
    @test length(batches[1]) == 2
    @test length(batches[2]) == 2
    @test length(batches[3]) == 2
    @test batches[1][1] == batches[1].x == X2[:,1:2]
    @test batches[1][2] == batches[1].y == Y2[1:2]
    @test batches[2][1] == batches[2].x == X2[:,3:4]
    @test batches[2][2] == batches[2].y == Y2[3:4]
    @test batches[3][1] == batches[3].x == X2[:,5:5]
    @test batches[3][2] == batches[3].y == Y2[5:5]

    # test iteration
    X3 = zeros(2, 10)
    d  = DataLoader(X3)
    for x in d
        @test size(x) == (2,1)
    end

    # test iteration
    X3 = ones(2, 10)
    Y3 = fill(5, 10)
    d  = DataLoader((X3, Y3))
    for (x, y) in d
        @test size(x) == (2,1)
        @test y == [5]
    end

    @testset "shuffle & rng" begin
        X4 = rand(2, 1000)
        d1 = DataLoader(X4, batchsize=2; shuffle=true)
        d2 = DataLoader(X4, batchsize=2; shuffle=true)
        @test first(d1) != first(d2)
        Random.seed!(17)
        d1 = DataLoader(X4, batchsize=2; shuffle=true)
        x1 = first(d1)
        Random.seed!(17)
        d2 = DataLoader(X4, batchsize=2; shuffle=true)
        @test x1 == first(d2)
        d1 = DataLoader(X4, batchsize=2; shuffle=true, rng=MersenneTwister(1))
        d2 = DataLoader(X4, batchsize=2; shuffle=true, rng=MersenneTwister(1))
        @test first(d1) == first(d2)
    end

    # numobs/getobs compatibility
    d = DataLoader(CustomType(), batchsize=2)
    @test first(d) == [1, 2]
    @test length(collect(d)) == 8

    @testset "Dict" begin
        data = Dict("x" => rand(2,4), "y" => rand(4))
        dloader = DataLoader(data, batchsize=2)
        @test_broken eltype(dloader) == Dict{String, Array{Float64}}
        c = collect(dloader)
        @test eltype(c) == Dict{String, Array{Float64}}
        @test c[1] == Dict("x" => data["x"][:,1:2], "y" => data["y"][1:2])
        @test c[2] == Dict("x" => data["x"][:,3:4], "y" => data["y"][3:4])

        data = Dict("x" => rand(2,4), "y" => rand(2,4))
        dloader = DataLoader(data, batchsize=2)
        @test_broken eltype(dloader) == Dict{String, Matrix{Float64}}
        @test eltype(collect(dloader)) == Dict{String, Matrix{Float64}}
    end

    @testset "range" begin
        data = 1:10

        dloader = DataLoader(data, batchsize=2)
        c = collect(dloader)
        @test eltype(c) == UnitRange{Int64}
        @test c[1] == 1:2

        dloader = DataLoader(data, batchsize=2, shuffle=true)
        c = collect(dloader)
        @test eltype(c) == Vector{Int}
    end
end


text:
"""
Computes ̂ρ , the distance between the satellite and the assumed user position

$SIGNATURES
´ξ´: Combination of estimated user position and time correction
´rSat´: Single satellite Positions
"""
function ρ_hat(ξ, rSats)
    rₙ = ξ[1:3]
    tc = ξ[4]

    ρ_hat = map(i -> norm(rSats[i] - rₙ) + tc, 1:length(rSats))
end

"""
Computes e, the direction vector from satellite to user

$SIGNATURES
´ξ´: Combination of estimated user position and time correction
´rSat´: Single satellite Positions
"""
function e(ξ, rSat)
    rₙ = ξ[1:3]
    e = (rₙ - rSat) / norm(rₙ - rSat)
end

"""
Computes Geometry Matrix H

$SIGNATURES
´ξ´: Combination of estimated user position and time correction
´rSats´: Matrix of satellite Positions
"""
function H(ξ, rSats)
    H = mapreduce(rSat -> [transpose(e(ξ, rSat)) 1], vcat, rSats)
end

"""
Calculates the dilution of precision for a given geometry matrix H

$SIGNATURES
´H_GEO´: Geometry matrix
"""
function calc_DOP(H_GEO)
    D = inv(H_GEO' * H_GEO)
    TDOP = sqrt(D[4,4]) # temporal dop
    VDOP = sqrt(D[3,3]) # vertical dop
    HDOP = sqrt(D[1,1] + D[2,2]) # horizontal dop
    PDOP = sqrt(D[1,1] + D[2,2] + D[3, 3]) # position dop
    GDOP = sqrt(sum(diag(D))) # geometrical dop

    return DOP(GDOP,PDOP,VDOP,HDOP,TDOP)
end

"""
Computes user position

$SIGNATURES
´rSats´: Array of Satellite positions. needs 3 values per satellite (xyz), size must be (3, N)")
´ρ´: Array of pseudo ranges

Calculates the user position by least squares method. The algorithm is based on the common reception method.
"""
function user_position(prev_pos ,rSats, ρ, accuracy::Float64 = 0.02)

    length(rSats) >= 4 || throw(SmallData("Not enough usable SV Data in user_position estimation."))

    # First Guesses of Position (Center of WGS84, time error = 0)
    #r₀ = [0.0, 0.0, 0.0]
    r₀ = prev_pos
    t = 0.0
    ξ = [r₀; t]

    Δξ = 100
    while norm(Δξ) > accuracy
        Δρ = ρ - ρ_hat(ξ, rSats)
        Δξ = H(ξ, rSats) \ Δρ #(transpose(H(ξ)) * H(ξ)) \ (transpose(H(ξ)) * Δρ) #H(ξ) \ Δρ
        ξ = ξ + Δξ # ξₙ₊₁ = ξₙ + Δξ
    end

    pos = ECEF(ξ[1:3])
    dt_receiver = ξ[4]
    DOP_val = calc_DOP(H(ξ, rSats))
    PVT = PVTSolution(pos, dt_receiver, DOP_val, [], [], [])
    return PVT
end


text:
println( "\nTest speed" )

const BENCH_REP = 200_000

function native_int()
    for i in 1:BENCH_REP
        @sprintf( "%10d", i )
    end
end
@static if VERSION >= v"1.6"
function format_int()
    fmt = Printf.Format( "%10d" )
    for i in 1:BENCH_REP
        Printf.format( fmt, i )
    end
end
end
function runtime_int()
    for i in 1:BENCH_REP
        cfmt( "%10d", i )
    end
end
function runtime_int_spec()
    fs = Format.FmtSpec("%10d")
    for i in 1:BENCH_REP
        cfmt( fs, i )
    end
end
function runtime_int_bypass()
    f = generate_formatter( "%10d" )
    for i in 1:BENCH_REP
        f( i )
    end
end

println( "integer @sprintf speed")
@time native_int()
@static if VERSION >= v"1.6"
    println( "integer format speed")
    @time format_int()
end
println( "integer cfmt speed")
@time runtime_int()
println( "integer cfmt spec speed")
@time runtime_int_spec()
println( "integer cfmt speed, bypass repeated lookup")
@time runtime_int_bypass()

set_seed!(10)
const testflts = [ _erfinv( rand() ) for i in 1:200_000 ]

function native_float()
    for v in testflts
        @sprintf( "%10.4f", v)
    end
end
@static if VERSION >= v"1.6"
function format_float()
    fmt = Printf.Format( "%10.4f" )
    for v in testflts
        Printf.format( fmt, v )
    end
end
end
function runtime_float()
    for v in testflts
        cfmt( "%10.4f", v)
    end
end
function runtime_float_spec()
    fs = Format.FmtSpec("%10.4f")
    for v in testflts
        cfmt( fs, v )
    end
end
function runtime_float_bypass()
    f = generate_formatter( "%10.4f" )
    for v in testflts
        f( v )
    end
end

println()
println( "float64 @sprintf speed")
@time native_float()
@static if VERSION >= v"1.6"
println( "float64 format speed")
@time format_float()
end
println( "float64 cfmt speed")
@time runtime_float()
println( "float64 cfmt spec speed")
@time runtime_float_spec()
println( "float64 cfmt speed, bypass repeated lookup")
@time runtime_float_bypass()


text:
@testset "Clusters" begin

    for int in 1:100
        @test getweight(NoClustering(), int) == 1.0
    end
end

@testset "Hobohm I" begin

    # DAWAEE
    # DAWAEF  83.3
    # DAWAED  83.3
    # DAYCMD  33.3
    # DAYCMT  33.3  83.3
    # DAYCMT  33.3  83.3

    fasta = read(joinpath(pwd(), "data", "Gaoetal2011.fasta"), FASTA)
    clusters = hobohmI(fasta, 62)

    @test nclusters(clusters) == 2
    @test nelements(clusters) == 6
    @test getweight(clusters, 1) == 1/3
    @test getweight(clusters, 6) == 1/3

    @testset "Clusters getters" begin

        @test getweight(clusters) == clusters.weights
        @test assignments(clusters) == clusters.clusters
        @test counts(clusters) == clusters.clustersize
    end

    @testset "Convert to Clusters" begin

        @test convert(Clusters, clusters) == clusters

        distance = convert(Matrix{Float64}, 100.0 .- percentidentity(fasta))
        cr = Clustering.dbscan(distance, 38.0, 2)
        @test convert(Clusters, cr) == clusters
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 41b66620-64d0-11eb-24a9-4f77cf2d6280
using Pkg

# ╔═╡ 7b921bd0-64a0-11eb-1683-998ca0325d4d
md"""# Introduction"""

# ╔═╡ 0d5727e0-64a1-11eb-2f03-ebae9b5bb6e6
md"""## Who are we?

* Lecturer: LCL IMM Ben Lauwens / D30.20 / [ben.lauwens@mil.be](mailto:ben.lauwens@mil.be)

* Assistant: CDT IMM François Dossogne / D30.25 / [francois.dossogne@mil.be](mailto:francois.dossogne@mil.be)"""

# ╔═╡ 4bfff530-64a1-11eb-2baa-270ad32c69f9
md"""## Programming Algorithms and Project

* Theory: 18 Hr 	→ 54 Hr (preparation of exercises)
* Exercises: 18 Hr 	→ 36 Hr (preparation of project)
* Project: 14 Hr 	→ 42 Hr (implementation of project)"""

# ╔═╡ ba22b0c2-64a1-11eb-2849-3d2d7f279e81
md"""## Documentation

* All notebooks can be found on github: [https://github.com/BenLauwens/ES123](https://github.com/BenLauwens/ES123)

* Complete book is available online: [https://benlauwens.github.io/ThinkJulia.jl/latest/book.html](https://benlauwens.github.io/ThinkJulia.jl/latest/book.html)."""

# ╔═╡ 1650e820-64a3-11eb-2e13-7b4fcebcdbc1
md"""## Schedule

### Theory

* 08/02: Ch 1 and Ch 2
* 11/02: Ch 3
* 18/02: Ch 4 and Ch 5
* 25/02: Ch 6 and Ch 7
* 04/03: Ch 8 and Ch 9
* 11/03: Ch 10
* 18/03: Ch 11 and Ch 12
* 25/03: Ch 13 and Ch 14
* 01/04: Ch 15, Ch 16 and Ch 17

### Exercises

* 22/02: Ch 1 and Ch 2
* 01/03: Ch 3
* 08/03: Ch 4 and Ch 5
* 15/03: Ch 6 and Ch 7
* 22/03: Ch 8 and Ch 9
* 29/03: Ch 10
* 12/04: Ch 11 and Ch 12
* 19/04: Ch 13 and Ch 14
* 22/04: Ch 15, Ch 16 and Ch 17

### Project

* 19/04: Project list available
* 22/04: Choice of project (groups of 2)
* 29/04: First contact: Understanding of the problem

We are always available during contact hours"""

# ╔═╡ b227f7b0-64a4-11eb-384b-8ba98563d161
md"""## Evaluation

Tests:
* 18/03 - 2Hr: Exercises
* 06-10/05: Intermediate evaluation of project (group of 2)

Exam: final evaluation of project (individual)"""

# ╔═╡ 325b6010-64a6-11eb-06bb-dd521b29756f
md"""## Julia

* Install Julia on CDN laptop
* Install Notepad++ on CDN laptop with the [Julia plugin](https://github.com/JuliaEditorSupport/julia-NotepadPlusPlus)
* Start Julia REPL and install Pluto:

```julia
julia> using Pkg

julia> pkg"add Pluto"
```

* Start Pluto:

```julia
julia> using Pluto

julia> Pluto.run()
```

* Open [https://localhost:1234]()

* Extra packages you will need:
"""

# ╔═╡ 4b982750-64d0-11eb-04ef-bdf336df89ad
pkg"add PlutoUI"

# ╔═╡ 56a63380-64d0-11eb-33e1-09c879da551c
pkg"add https://github.com/BenLauwens/NativeSVG.jl.git"

# ╔═╡ 045aeb3e-8a91-11eb-05c2-7bfe69d7f97f

# ╔═╡ Cell order:
# ╟─7b921bd0-64a0-11eb-1683-998ca0325d4d
# ╟─0d5727e0-64a1-11eb-2f03-ebae9b5bb6e6
# ╟─4bfff530-64a1-11eb-2baa-270ad32c69f9
# ╟─ba22b0c2-64a1-11eb-2849-3d2d7f279e81
# ╟─1650e820-64a3-11eb-2e13-7b4fcebcdbc1
# ╟─b227f7b0-64a4-11eb-384b-8ba98563d161
# ╟─325b6010-64a6-11eb-06bb-dd521b29756f
# ╠═41b66620-64d0-11eb-24a9-4f77cf2d6280
# ╠═4b982750-64d0-11eb-04ef-bdf336df89ad
# ╠═56a63380-64d0-11eb-33e1-09c879da551c
# ╠═045aeb3e-8a91-11eb-05c2-7bfe69d7f97f


text:
using Revise
using AlgebraicPetri
using OrdinaryDiffEq
using DiffEqFlux, Flux
using Catalyst
using ModelExploration

SIR =  LabelledPetriNet([:S, :I, :R],
                        :inf => ((:S, :I)=>(:I, :I)),
                        :rec => (:I=>:R));

SIR2 =  LabelledPetriNet([:S, :I, :R],
                        :inf => (:S=>:I),
                        :rec => (:I=>:R));

initial_conc=[1e7,1e5,1e-10]
true_rates = [1e-8,1e-4]
data = generate_data(SIR, true_rates,initial_conc)

final_loss, fitted_params = train(SIR, data)

predicted_traj = generate_data(SIR, fitted_params[1:2],fitted_params[3:5])
plot(data; xaxis="days")
plot!(predicted_traj; xaxis="days")

f = eval_petri_fn(SIR, true_rates, initial_conc)
v1= f(SIR)
v2 = f(SIR2)


text:
"
    getdata_arlett_noflow([Int]) -> DataFrame

Disable flow bias/wall-catalase, run the Arlett simulation, and revert controls.

For each walk that ends in a sensor collision, return
    1. the trial #
    2. inner sensor cross talk, average of left and right
    3. outer sensor cross talk, average of left and right
in a DataFrame.

Simulation controls:
    set_NUMBER_OF_WALKS()
    set_MAX_STEPS_PER_WALK()
"
function getdata_arlett_noflow(seed::Int=randseed())
    Random.seed!(seed)
    float_arr = rand(NUMBER_OF_WALKS)
    stub_arr = rand(NUMBER_OF_WALKS) # for consistency with run_sim!()

    init_FLOW_BIAS = FLOW_BIAS
    init_CATALASE = CATALASE_ON_WALLS
    set_FLOW_BIAS(false)
    set_CATALASE_ON_WALLS(false)

    arlett_print(seed)
    n_array = []
    innercrosstalk_array = []
    outercrosstalk_array = []
    center, inner, outer = (0, 0, 0)
    innercrosstalk = 0.0; outercrosstalk = 0.0
    # Record data between trials instead of only the end using run_sim!()
    @showprogress for n in 1:NUMBER_OF_WALKS
        peroxidexy = spawnrandompoint(float_arr[n])
        steps_sofar = 0
        while peroxidexy != undef && steps_sofar < MAX_STEPS_PER_WALK
            peroxidexy, collision = one_step!(peroxidexy)
            steps_sofar += 1
            if collision == "left outer sensor" || collision == "right outer sensor"
                outer += 1
            elseif collision == "left inner sensor" || collision == "right inner sensor"
                inner += 1
            elseif collision == "center sensor"
                center += 1
            end
        end
        if center > 0
            innercrosstalk = 50.0 * inner / center
            outercrosstalk = 50.0 * outer / center
        else
            if inner > 0
                innercrosstalk = 100.0
            elseif outer > 0
                outercrosstalk = 100.0
            end
        end
        append!(n_array, n)
        append!(innercrosstalk_array, innercrosstalk)
        append!(outercrosstalk_array, outercrosstalk)
    end
    df = DataFrame(
        nth_trial = n_array,
        inner_crosstalk = innercrosstalk_array,
        outer_crosstalk = outercrosstalk_array
    )
    current_df(df)
    current_seed(seed)
    current_path("out/noflowdata/")
    set_FLOW_BIAS(init_FLOW_BIAS)
    set_CATALASE_ON_WALLS(init_CATALASE)
    return df
end

max_stars_repo_name: s-y-wu/HMCResearchRandomWalks

max_issues_repo_name: s-y-wu/HMCResearchRandomWalks

max_forks_repo_name: s-y-wu/HMCResearchRandomWalks


text:
module Checkpoint

export write_checkpoint, rm_checkpoint, read_checkpoint

using JLD2
using MPI
using Printf

using ..ODESolvers

"""
    write_checkpoint(solver_config, checkpoint_dir, name, mpicomm, num)

Read in the state and auxiliary arrays as well as the simulation time
stored in the checkpoint file for `name` and `num`.
"""
function write_checkpoint(
    solver_config,
    checkpoint_dir::String,
    name::String,
    mpicomm::MPI.Comm,
    num::Int,
)
    nm = replace(name, " " => "_")
    cname = @sprintf(
        "%s_checkpoint_mpirank%04d_num%04d.jld2",
        nm,
        MPI.Comm_rank(mpicomm),
        num,
    )
    cfull = joinpath(checkpoint_dir, cname)
    @info @sprintf(
        """
Checkpoint
    saving to %s""",
        cfull
    )

    dg = solver_config.dg
    Q = solver_config.Q
    if Array ∈ typeof(Q).parameters
        h_Q = Q.realdata
        h_aux = dg.state_auxiliary.realdata
    else
        h_Q = Array(Q.realdata)
        h_aux = Array(dg.state_auxiliary.realdata)
    end
    t = ODESolvers.gettime(solver_config.solver)
    @save cfull h_Q h_aux t

    return nothing
end

"""
    rm_checkpoint(checkpoint_dir, name, mpicomm, num)

Remove the checkpoint file identified by `solver_config.name` and `num`.
"""
function rm_checkpoint(
    checkpoint_dir::String,
    name::String,
    mpicomm::MPI.Comm,
    num::Int,
)
    nm = replace(name, " " => "_")
    cname = @sprintf(
        "%s_checkpoint_mpirank%04d_num%04d.jld2",
        nm,
        MPI.Comm_rank(mpicomm),
        num,
    )
    rm(joinpath(checkpoint_dir, cname), force = true)

    return nothing
end

"""
    read_checkpoint(checkpoint_dir, name, mpicomm, num)

Read in the state and auxiliary arrays as well as the simulation time
stored in the checkpoint file for `name` and `num`.
"""
function read_checkpoint(
    checkpoint_dir::String,
    name::String,
    array_type,
    mpicomm::MPI.Comm,
    num::Int,
)
    nm = replace(name, " " => "_")
    cname = @sprintf(
        "%s_checkpoint_mpirank%04d_num%04d.jld2",
        nm,
        MPI.Comm_rank(mpicomm),
        num,
    )
    cfull = joinpath(checkpoint_dir, cname)
    if !isfile(cfull)
        error("Cannot restore from checkpoint in %s, file not found")
    end

    @load cfull h_Q h_aux t

    return (array_type(h_Q), array_type(h_aux), t)
end

end # module


text:
using PowerSystems
using Plots

sys = System(0, enable_compression=true)
add_component!(sys, Bus(nothing))
for i in 1:3
    plant = RenewableDispatch(nothing)
    set_name!(plant, "plant$i")
    add_component!(sys, plant)
end

# required_metadata

metadata_file = "data/siip_timeseries_metadata.json"
add_time_series!(sys, metadata_file)
to_json(sys, "basic_system/sys.json", force=true)

ta = get_time_series_array(
    SingleTimeSeries,
    get_component(RenewableDispatch, sys, "plant1"),
    "max_active_power",
)
pyplot()
savefig(plot(ta), "basic_system/plant1_max_active_power.png")

max_stars_repo_name: NREL-SIIP/reV-PowerSystems

max_issues_repo_name: NREL-SIIP/reV-PowerSystems

max_forks_repo_name: NREL-SIIP/reV-PowerSystems


text:
import Photometry.Background: ShepardIDWInterpolator

@testset "IDWInterpolator" begin
    @test IDWInterpolator((2, 3), k = 4)(ones(3, 2)) == ones(6, 6)
    @test IDWInterpolator(2)(ones(3, 4)) == IDWInterpolator((2, 2))(ones(3, 4))
    @test_throws ErrorException IDWInterpolator(2)(ones(3, 2))

    # interface
    @test IDWInterpolator((2, 3)) == IDWInterpolator(2, 3)
    @test IDWInterpolator((2, 3), k = 4) == IDWInterpolator(2, 3, k = 4)
end

@testset "ShepardIDWInterpolator" begin
    # 1D
    x = rand(100)
    y = sin.(x)
    it = ShepardIDWInterpolator(x', y)
    @test it(0.4) ≈ sin(0.4) atol = 1e-2
    @test it.(0.2:0.1:0.5) ≈ sin.(0.2:0.1:0.5) atol = 1e-2

    # 2D
    x = rand(2, 10000)
    y = sin.(x[1, :] .+ x[2, :])
    it = ShepardIDWInterpolator(x, y)
    @test it(0.5, 0.6) ≈ sin(0.5 + 0.6) atol = 1e-2

    # interface

    #= Warning! These are not accurate for use as a standard interpolator,
       but are what we need for our use with images =#
    @test size(it) == (10000,)
    @test axes(it) == axes(y)
end


text:
@testset "Lane" begin

    # Test circle
    x = 1.0
    y = 2.0
    r = 10.0
    c = CircularWall(x,y,r)
    @test typeof(c) <: CircularWall

    c1 = CircularWall()
    @test c1.x == 0.0
    @test c1.y == 0.0
    @test c1.r == 0.0

    # Test Lane
    id = 1
    walls = [Algames.Wall([0.1, 0.2], [0.3, 0.4], [0.5, 0.5]) for i=1:10]
    circles = [CircularWall(0.1, 0.2, 0.3) for i=1:10]
    start = StartingArea()
    f = x -> x
    lane = Lane(id, :lane_1, walls, circles, start, f)
    @test typeof(lane) <:Lane

end


text:
const DO_NOT_PUSH_TO_THESE_DESTINATIONS = String[
    "MirrorUpdater",
    "MirrorUpdater.jl",
    "UnofficialJuliaMirrorSnapshots/MirrorUpdater",
    "UnofficialJuliaMirrorSnapshots/MirrorUpdater.jl",
    "UnofficialJuliaMirrorBot/MirrorUpdater",
    "UnofficialJuliaMirrorBot/MirrorUpdater.jl",
    #
    "RepoSnapshots",
    "RepoSnapshots.jl",
    "UnofficialJuliaMirrorSnapshots/RepoSnapshots",
    "UnofficialJuliaMirrorSnapshots/RepoSnapshots.jl",
    "UnofficialJuliaMirrorBot/RepoSnapshots",
    "UnofficialJuliaMirrorBot/RepoSnapshots.jl",
    ]


text:
export securitydice, normaldice, riskydice

"""
    securitydice(list)

Computes the transition probability matrix for the security dice
with squares types defined in `list`.
"""
function securitydice(list::Vector{Int64})::Array{Float64,2}
    proba = zeros(15, 15)
    for i = 1:14
        proba[i,i] = 0.5
        if i == 3
            proba[i,i+1] = 0.25
            proba[i,i+8] = 0.25
        elseif i == 10
            proba[i,i+5] = 0.5
        else
            proba[i,i+1] = 0.5
        end
    end
    proba[15,15] = 1
    return proba
end

"""
    normaldice(list, circular)

Computes the transition probability matrix for the normal dice
with squares types defined in `list` and circularity of
board defined by `circular`.
"""
function normaldice(list::Vector{Int64}, circular::Bool)::Array{Float64,2}
    proba = zeros(15, 15)
    # Compute first the transition proba matrix as if there were no traps
    p = 1/3
    for i = 1:14
        proba[i,i] = p
        if i == 3
            proba[i,i+1] = p/2
            proba[i,i+2] = p/2
            proba[i,i+8] = p/2
            proba[i,i+9] = p/2
        elseif i == 9
            proba[i,i+1] = p
            proba[i,i+6] = p
        elseif i == 10
            if circular
                proba[i,i+5] = p
                proba[i,(i+6) % 15] = p
            else
                proba[i,i+5] = 2*p
            end
        elseif i == 14
            if circular
                proba[i,i+1] = p
                proba[i,(i+2) % 15] = p
            else
                proba[i,i+1] = 2*p
            end
        else
            proba[i,i+1] = p
            proba[i,i+2] = p
        end
    end
    proba[15,15] = 1

    copy_proba = copy(proba)
    # Remodify each line to take traps into account
    for i = 1:14
        for j = 1:15
            if list[j] == 1 # Trap of type 1 : go back to square 1
                if j != 1
                    proba[i,1] = proba[i,1] + proba[i,j]/2
                    proba[i,j] = proba[i,j]/2
                end
            elseif list[j] == 2 # Trap of type 2 : go back three squares
                proba[i,squaresub(j, 3)] = proba[i,squaresub(j, 3)] + proba[i,j]/2
                proba[i,j] = proba[i,j]/2
            end
        end
    end
    return proba
end

"""
    riskydice(list, circular)

Computes the transition probability matrix for the risky dice
with squares types defined in `list` and circularity of
board defined by `circular`.
"""
function riskydice(list::Vector{Int64}, circular::Bool)::Array{Float64,2}
    proba = zeros(15, 15)
    # Compute first the transition proba matrix as if there were no traps
    for i = 1:14
        proba[i,i] = 0.25
        if i == 3
            proba[i,i+1] = 0.125
            proba[i,i+2] = 0.125
            proba[i,i+3] = 0.125
            proba[i,i+8] = 0.125
            proba[i,i+9] = 0.125
            proba[i,i+10] = 0.125
        elseif i == 8
            proba[i,i+1] = 0.25
            proba[i,i+2] = 0.25
            proba[i,i+7] = 0.25
        elseif i == 9
            proba[i,i+1] = 0.25
            if circular
                proba[i,i+6] = 0.25
                proba[i,i-8] = 0.25
            else
                proba[i,i+6] = 0.5
            end
        elseif i == 10
            if circular
                proba[i,i+5] = 0.25
                proba[i,i-9] = 0.25
                proba[i,i-8] = 0.25
            else
                proba[i,i+5] = 0.75
            end
        elseif i == 13
            proba[i,i+1] = 0.25
            if circular
                proba[i,i+2] = 0.25
                proba[i,i-12] = 0.25
            else
                proba[i,i+2] = 0.5
            end
        elseif i == 14
            if circular
                proba[i,i+1] = 0.25
                proba[i,i-13] = 0.25
                proba[i,i-12] = 0.25
            else
                proba[i,i+1] = 0.75
            end
        else
            proba[i,i+1] = 0.25
            proba[i,i+2] = 0.25
            proba[i,i+3] = 0.25
        end
    end
    proba[15,15] = 1

    copy_proba = copy(proba)
    # Remodify each line to take traps into account
    for i = 1:14
        for j = 1:15
            if list[j] == 1 # Trap of type 1 : go back to square 1
                if j != 1
                    proba[i,1] = proba[i,1] + proba[i,j]
                    proba[i,j] = 0
                end
            elseif list[j] == 2 # Trap of type 2 : go back three squares
                proba[i, squaresub(j, 3)] = proba[i, squaresub(j, 3)] + proba[i,j]
                proba[i,j] = 0
            end
        end
    end
    return proba
end

max_stars_repo_name: qlete/markov-decision

max_issues_repo_name: qlete/markov-decision

max_forks_repo_name: qlete/markov-decision


text:
using BayesNN
using Base.Test

include("constructors.jl")
include("fit.jl")


text:
function invariants_Q6_check(x)
pv=zeros(6,1);

v=zeros(6,1);
 v[1] = 1
 pv[1] = x[1]^2*x[2] + x[1]*x[2]^2 + x[1]^2*x[4] +
x[2]^2*x[4] + x[1]*x[4]^2 + x[2]*x[4]^2 +
x[1]^2*x[3] + x[2]^2*x[3] + x[1]*x[3]^2 +
x[2]*x[3]^2 + x[1]^2*x[5] + x[4]^2*x[5] +
x[3]^2*x[5] + x[1]*x[5]^2 + x[4]*x[5]^2 +
x[3]*x[5]^2 + x[2]^2*x[6] + x[4]^2*x[6] +
x[3]^2*x[6] + x[5]^2*x[6] + x[2]*x[6]^2 +
x[4]*x[6]^2 + x[3]*x[6]^2 + x[5]*x[6]^2
 v[2:2] = pv[1:1]
 pv[2] = x[1]^3*x[2] + x[1]*x[2]^3 + x[1]^3*x[4] +
x[2]^3*x[4] + x[1]*x[4]^3 + x[2]*x[4]^3 +
x[1]^3*x[3] + x[2]^3*x[3] + x[1]*x[3]^3 +
x[2]*x[3]^3 + x[1]^3*x[5] + x[4]^3*x[5] +
x[3]^3*x[5] + x[1]*x[5]^3 + x[4]*x[5]^3 +
x[3]*x[5]^3 + x[2]^3*x[6] + x[4]^3*x[6] +
x[3]^3*x[6] + x[5]^3*x[6] + x[2]*x[6]^3 +
x[4]*x[6]^3 + x[3]*x[6]^3 + x[5]*x[6]^3
 v[3:3] = pv[2:2]
 pv[3] = x[1]^5 + x[2]^5 + x[4]^5 + x[3]^5 + x[5]^5 +
x[6]^5
 v[4:4] = pv[3:3]
 v[5] = pv[1]*pv[1]
 v[6] = pv[2]*pv[3]
 Primary_invariants= [ x[1] + x[2] + x[3] + x[4] + x[5] + x[6], x[1]*x[6] + x[2]*x[5] + x[3]*x[4], x[1]^2 + x[2]^2 + x[3]^2 + x[4]^2 + x[5]^2 + x[6]^2, x[1]*x[2]*x[3] + x[1]*x[4]*x[5] + x[2]*x[4]*x[6] + x[3]*x[5]*x[6], x[1]^3 + x[2]^3 + x[3]^3 + x[4]^3 + x[5]^3 + x[6]^3, x[1]^4 + x[2]^4 + x[3]^4 + x[4]^4 + x[5]^4 + x[6]^4 ]

return Primary_invariants, v

end
x = rand(6)
display(invariants_Q6_check(x))

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
module StrideArraysCore

using VectorizationBase, ArrayInterface, ThreadingUtilities
using ArrayInterface: StaticInt, Zero, One, StaticBool, True, False,
    OptionallyStaticUnitRange, size, strides, offsets, indices,
    static_length, static_first, static_last, axes,
    dense_dims, stride_rank
using VectorizationBase: align, gep, AbstractStridedPointer, AbstractSIMDVector, vnoaliasstore!, staticm1,
    static_sizeof, lazymul, vmul_fast, StridedPointer, gesp, zero_offsets, pause, zstridedpointer,
    val_dense_dims, val_stride_rank, preserve_buffer

export PtrArray, StrideArray, StaticInt

include("closeopen.jl")
include("ptr_array.jl")
include("stridearray.jl")
include("thread_compatible.jl")
include("views.jl")
include("adjoints.jl")

function __init__()
    # @require LoopVectorization="bdcacae8-1622-11e9-2a5c-532679323890" @eval using StrideArrays
end

end


text:
function nicoud(L,z;maxiter=10,tol=0.,relax=1.,n_eig_val=3,v0=[],output=true)
    if output
        println("Launching Nicoud...")
        println("Iter    Res:     dz:     z:")
        println("----------------------------------")
    end
    z0=complex(Inf)
    n=0
    flag=itsol_converged

    M=L(1,oplist=["M"],in_or_ex=true) #get mass matrix
    K=L(1,oplist=["K"],in_or_ex=true)
    C=L(1,oplist=["C"],in_or_ex=true)
    d=size(M)[1] #TODO: implement a dimension method for LinearOperatorFamily
    I=SparseArrays.sparse(SparseArrays.I,d,d)
    O=SparseArrays.spzeros(d,d)
    Y=[I O; O M]
    if v0==[]
        v0=ones(ComplexF64,d)
    end
    v0=[v0; z*v0]
    try
        while abs(z-z0)>tol && n< maxiter
            if output; println(n,"\t\t",abs(z-z0),"\t", z );flush(stdout); end
            z0=z
            Q=L(z,oplist=["Q"],in_or_ex=true)
            X=[O -I; K+Q C]
            z,v0=Arpack.eigs(-X,Y,sigma=z0,v0=v0,nev=n_eig_val)
            indexing=sortperm(z.-z0,by=abs)
            z,v0=z[indexing[1]],v0[:,indexing[1]]
            z=z0+relax*(z-z0)
            #TODO consider relaxation in v0
            n+=1
        end
    catch excp
        if output
            println("Error occured:")
            println(excp)
            println(typeof(excp))
            println("...aborted Nicoud!")
        end

        flag=itsol_unknown
        if typeof(excp) <: Arpack.ARPACKException
           flag=itsol_arpack_exception
           if excp==Arpack.ARPACKException(-9999)
               flag=itsol_arpack_9999
           end
        elseif excp== LinearAlgebra.SingularException(0) #This means that the solution is so good that L(z) cannot be LU factorized...TODO: implement other strategy into perturb
          flag=itsol_singular_exception
          L.params[L.eigval]=z
        end
    end
    if flag==itsol_converged
        L.params[L.eigval]=z
        if output;  println(n,"\t\t",abs(z-z0),"\t", z ); end

        if n>=maxiter
            flag=itsol_maxiter
            if output; println("Warning: Maximum number of iterations has been reached!");end

        elseif abs(z-z0)<=tol
            flag=itsol_converged
            if output; println("Solution has converged!"); end
        elseif isnan(z)
            flag=itsol_isnan
            if output; println("Warning: computer arithmetics problem. Eigenvalue is NaN"); end
        else
            if output; println("Warning: This should not be possible....\n If you can read this contact GAM!");end
            flag=itsol_impossible
        end
        if output
            println("...finished Nicoud!")
            println("#####################")
            println(" Nicoud results ")
            println("#####################")
            println("Number of steps: ",n)
            println("Last step parameter variation:",abs(z0-z))
            println("Eigenvalue:",z)
            println("Eigenvalue/(2*pi):",z/2/pi)
        end
    end
    return Solution(L.params,v0[1:d],[],L.eigval), n, flag
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 6f9a94bd-dd28-4110-a7ca-0ed84e9c7c3f
begin
	import Pkg
    # activate the shared project environment
    Pkg.activate(Base.current_project())
	using Interpolations
	using Plots
	using Dierckx
	using GasChromatographySimulator
	using GasChromatographyTools
	using Plots
	using PlutoUI
	TableOfContents()
end

# ╔═╡ d3c2c2f4-1100-438d-88d6-d7502bcb9310
using Statistics

# ╔═╡ 0b51a792-d9b8-4c29-928f-57aaf41f1c20
plotly()

# ╔═╡ 93ba6afc-4e9a-11ec-08d9-81c0a9bc502e
md"""
# Test of the RT lock functions
The test is necessary, because of failing (infinit loop in case of opt\_ipt="linear", NaN abort in case of opt\_itp="spline").
"""

# ╔═╡ 51ec0223-199a-4a25-8768-cd0b7e9f864d
begin
	function stretched_program(n::Float64, par::GasChromatographySimulator.Parameters)
		# stretch the temperature program in 'par' by a factor 'n'
		if isa(par, Array)==true
			error("Select an element of the array of GC-system parameters.")
		else
			new_tsteps = n.*par.prog.time_steps
			new_T_itp = GasChromatographySimulator.temperature_interpolation(new_tsteps, par.prog.temp_steps, par.prog.gf, par.col.L)
			new_pin_itp = GasChromatographySimulator.pressure_interpolation(new_tsteps, par.prog.pin_steps)
			new_pout_itp = GasChromatographySimulator.pressure_interpolation(new_tsteps, par.prog.pout_steps)
			new_prog = GasChromatographySimulator.Program(new_tsteps, par.prog.temp_steps, par.prog.pin_steps, par.prog.pout_steps, par.prog.gf, par.prog.a_gf, new_T_itp, new_pin_itp, new_pout_itp)
			new_par = GasChromatographySimulator.Parameters(par.col, new_prog, par.sub, par.opt)
			return new_par
		end
	end

	function initial_n(n::Float64, tR_lock::Float64, ii::Int, par::GasChromatographySimulator.Parameters)
	    par_n = stretched_program(n, par)
	    sol_n = GasChromatographySimulator.solve_system_multithreads(par_n)
	    tR_n = sol_n[ii].u[end][1]
	    if n>1
	        while tR_n-tR_lock<0 && n<130.0
	            n = n*2.0
	            par_n = stretched_program(n, par)
	            sol_n = GasChromatographySimulator.solve_system_multithreads(par_n)
	            tR_n = sol_n[ii].u[end][1]
	        end
	    elseif n<1
	        while tR_n-tR_lock>0 && n>0.01
	            n = n*0.5
	            par_n = stretched_program(n, par)
	            sol_n = GasChromatographySimulator.solve_system_multithreads(par_n)
	            tR_n = sol_n[ii].u[end][1]
	        end
	    end
	    if n>130.0
	        error("The choosen retention time for locking is to big.")
	    elseif n<0.01
	        error("The choosen retention time for locking is to small.")
	    else
	        return n
	    end
	end

	function RT_locking(par::GasChromatographySimulator.Parameters, tR_lock::Float64, tR_tol::Float64, solute_RT::String; opt_itp="linear")
		# estimate the factor 'n' for the temperature program to achieve the retention time 'tR_lock' for 'solute_RT' with the GC-system defined by 'par'
		if isa(par, Array)==true
			error("Select an element of the array of GC-systems.")
		else
			# find 'solute_RT' in the substances of 'par'
			name = Array{String}(undef, length(par.sub))
			for i=1:length(par.sub)
				name[i] = par.sub[i].name
			end
			ii = findfirst(name.==solute_RT)
			# calculate the retention time for the original (un-stretched) program
			sol₀ = GasChromatographySimulator.solving_odesystem_r(par.col, par.prog, par.sub[ii], par.opt)
			tR₀ = sol₀.u[end][1]
			# start value for the factor 'n'
			if tR₀-tR_lock<0
				n₁ = initial_n(2.0, tR_lock, ii, par)
			else
				n₁ = initial_n(0.5, tR_lock, ii, par)
			end
			# using a recursive function to estimate 'n'
			n = recur_RT_locking(n₁, [1.0], [tR₀], par, tR_lock, tR_tol, ii; opt_itp=opt_itp)
		end
		return n
	end

	function recur_RT_locking(n::Float64, n_vec::Array{Float64,1}, tR_vec::Array{Float64,1}, par::GasChromatographySimulator.Parameters, tR_lock::Float64, tR_tol::Float64, ii::Int64; opt_itp="linear")
		# recursive function to find the factor 'n' for the temperature program to achieve the retention time 'tR_lock' for solute index 'ii' with the GC-system defined by 'par'
		# calculate the retention time with the input guess 'n'
		par₁ = stretched_program(n, par)
		if par.opt.odesys==true
			sol₁ = GasChromatographySimulator.solving_odesystem_r(par₁.col, par₁.prog, par₁.sub[ii], par₁.opt)
			tR₁ = sol₁.u[end][1]
		else
			sol₁ = GasChromatographySimulator.solving_migration(par₁.col, par₁.prog, par₁.sub[ii], par₁.opt)
			tR₁ = sol₁.u[end]
		end
		if abs(tR₁-tR_lock)<tR_tol
			# if retention time is less than 'tR_tol' from 'tR_lock' we found the factor 'n'
			return n
		else
			# estimate a new factor 'new_n' by linear interpolation of the factors 'n_vec' + 'n' over the corresponding retention times
			new_n_vec = sort([n_vec; n])
			new_tR_vec = sort([tR_vec; tR₁])
			if opt_itp=="spline"
				# Dierckx.jl
				if length(new_tR_vec)<4
					k = length(new_tR_vec)-1
				else
					k = 3
				end
				itp = Spline1D(new_tR_vec, new_n_vec, k=k)
			else # opt_itp=="linear"
				# Interpolations.jl
				itp = LinearInterpolation(sort([tR_vec; tR₁]), sort([n_vec; n]))
			end
			new_n = itp(tR_lock)
			#println("new_n=$(new_n), tR₁=$(tR₁)")
			# use the new factor 'new_n' and call the recursive function again
			return recur_RT_locking(new_n, new_n_vec, new_tR_vec, par, tR_lock, tR_tol, ii; opt_itp=opt_itp)
		end
	end

	md"""
	## Copies of the RT lock functions
	In the package GasChromatographyTools.jl the functions where adapted to solve this problem.
	"""
end

# ╔═╡ 5e642f09-9a8f-4cca-b61f-b27c8433a2e5
begin
	opt = GasChromatographySimulator.Options(OwrenZen5(), 1e-6, 1e-4, "inlet", true)
	L = 4.0
	d = 0.1e-3
	df = 0.1e-6
	sp = "Rxi17SilMS" # ["Rxi17SilMS" -> ERR, "SLB5ms", "SPB50", "Wax", "DB5ms", "Rxi5MS", "genericLB", "genericJL"]
	gas = "He"
	col = GasChromatographySimulator.constructor_System(L, d, df, sp, gas)

	db_path = "/Users/janleppert/Documents/GitHub/Publication_GCsim/data/Databases/"
	db_file = "Database_append.csv"
	first_alkane = "C8"
	last_alkane = "C15"
	sub = GasChromatographySimulator.load_solute_database(db_path, db_file, sp, gas, [first_alkane, last_alkane], zeros(2), zeros(2))

	Tst = 273.15
	Tref = 150.0 + Tst
	pn = 101300.0
	pin = 300000.0 + pn
	pout = 0.0
	dimless_rate = 0.4
	Theat = 1000.0
	Tshift = 40.0
	tMref = GasChromatographySimulator.holdup_time(Tref, pin, pout, L, d, gas)
	rate = dimless_rate*30/tMref
	theat = Theat/rate
	Tstart = sub[1].Tchar-Tst-Tshift

	ΔT = 135.0
	α = 9.0
	prog0 = GasChromatographySimulator.constructor_Program([0.0, theat],[Tstart, Tstart+Theat], [pin, pin],[pout, pout],[ΔT, ΔT], [0.0, 0.0], [L, L], [α, α], opt.Tcontrol, L)

	par0 = GasChromatographySimulator.Parameters(col, prog0, sub, opt)

	tR_lock = 12*tMref
	tR_tol = 1e-3
	md"""
	## Settings
	The following settings produce the problem:
	"""
end

# ╔═╡ 3de48bb6-8eb7-4a33-9a98-d5fe3a19f6c6
md"""
## The ERROR
"""

# ╔═╡ 19d1abb6-32b6-414f-bc95-55635cbaa73a
n = RT_locking(par0, tR_lock, 1e-3, last_alkane; opt_itp="spline")

# ╔═╡ 2cad8185-4fc5-4009-98df-07a2be2133c6
md"""
## Reproducing the error by stepwise calculations
"""

# ╔═╡ 30207849-cafe-4ec6-af8d-ee7b2e2e6de0
begin
	# initial simulation
	sol₀ = GasChromatographySimulator.solving_odesystem_r(par0.col, par0.prog, par0.sub[2], par0.opt)
	tR₀ = sol₀.u[end][1]
end

# ╔═╡ 84e7e869-0fbf-4590-a8de-28855856661f
tR₀-tR_lock

# ╔═╡ c5ace7ce-cfa3-4c15-bcc1-8b66ad1c16e9
begin
	# first stretch of the program
	if tR₀-tR_lock<0
	    n₁ = initial_n(2.0, tR_lock, 2, par0)
	else
	    n₁ = initial_n(0.5, tR_lock, 2, par0)
	end
	# recur function
	par₁ = stretched_program(n₁, par0)
	sol₁ = GasChromatographySimulator.solving_odesystem_r(par₁.col, par₁.prog, par₁.sub[2], par₁.opt)
	tR₁ = sol₁.u[end][1]
	n_vec₁ = sort([1.0; n₁])
	tR_vec₁ = sort([tR₀; tR₁])
	# first interpolation
	itp₁ = Spline1D(tR_vec₁, n_vec₁, k=1)
end

# ╔═╡ c4d53222-03b1-4ce4-a49a-690945347432
tR₁-tR_lock

# ╔═╡ 6acfbf9b-f8ac-4483-8c93-17920d0d9f0e
begin
	# second stretch
	n₂ = itp₁(tR_lock)
	par₂ = stretched_program(n₂, par0)
	sol₂ = GasChromatographySimulator.solving_odesystem_r(par₂.col, par₂.prog, par₂.sub[2], par₂.opt)
	tR₂ = sol₂.u[end][1]
	n_vec₂ = sort([n_vec₁; n₂])
	tR_vec₂ = sort([tR_vec₁; tR₂])
	itp₂ = Spline1D(tR_vec₂, n_vec₂, k=2)

	p1 = plot([tR_lock, tR_lock], [0.4, 0.6], label="tR_lock")
	scatter!(p1, [tR₂, tR₂], [n₂, n₂], label="2")
	plot!(p1, 40.0:0.001:60.0, itp₂.(40.0:0.001:60.0), xlims=(50.0,55.0), ylims=(0.40, 0.60), label="itp₂")
	p1
end

# ╔═╡ 3b31f580-efe2-4a7a-89dc-228b38f2a71e
tR₂-tR_lock

# ╔═╡ b6c5126a-3d41-4809-ab9b-d554262e0668
begin
	n₃ = itp₂(tR_lock)
	par₃ = stretched_program(n₃, par0)
	sol₃ = GasChromatographySimulator.solving_odesystem_r(par₃.col, par₃.prog, par₃.sub[2], par₃.opt)
	tR₃ = sol₃.u[end][1]
	n_vec₃ = sort([n_vec₂; n₃])
	tR_vec₃ = sort([tR_vec₂; tR₃])
	k = length(tR_vec₃)-1
	itp₃ = Spline1D(tR_vec₃, n_vec₃, k=k)

	scatter!(p1, [tR₃, tR₃], [n₃, n₃], label="3")
	plot!(p1, 51.0:0.001:52.0, itp₃.(51.0:0.001:52.0), label="itp₃")
end

# ╔═╡ 6c5443bc-3dab-4317-99c9-1bb32b184fbc
tR₃-tR_lock

# ╔═╡ 3b1bb187-66f9-40bd-a0a5-a3c4e1a23819
begin
	n₄ = itp₃(tR_lock)
	par₄ = stretched_program(n₄, par0)
	sol₄ = GasChromatographySimulator.solving_odesystem_r(par₄.col, par₄.prog, par₄.sub[2], par₄.opt)
	tR₄ = sol₄.u[end][1]
	n_vec₄ = sort([n_vec₃; n₄])
	tR_vec₄ = sort([tR_vec₃; tR₄])
	itp₄ = Spline1D(tR_vec₄, n_vec₄, k=3)

	scatter!(p1, [tR₄, tR₄], [n₄, n₄], label="4")
	plot!(p1, 51.0:0.001:52.0, itp₄.(51.0:0.001:52.0), label="itp₄")
end

# ╔═╡ b81b2905-6756-445b-b1be-c54298df8b3f
tR₄-tR_lock

# ╔═╡ eb2d9f6a-8384-427a-a036-b4f0019d8251
md"""
The proposed settings seam to lead to jumping estimations around the searched retention time.
"""

# ╔═╡ 734bae50-12fa-4717-8446-addb963b8673
begin
	n₅ = itp₄(tR_lock)
	par₅ = GasChromatographyTools.stretched_program(n₅, par0)
	sol₅ = GasChromatographySimulator.solving_odesystem_r(par₅.col, par₅.prog, par₅.sub[2], par₅.opt)
	tR₅ = sol₅.u[end][1]
	n_vec₅ = sort([n_vec₄; n₅])
	tR_vec₅ = sort([tR_vec₄; tR₅])
	itp₅ = Spline1D(tR_vec₅, n_vec₅, k=3)

	scatter!(p1, [tR₅, tR₅], [n₅, n₅], label="5")
	plot!(p1, 51.0:0.001:52.0, itp₅.(51.0:0.001:52.0), xlims=(51.47,51.485), ylims=(0.4995, 0.49975), label="itp₅")
end

# ╔═╡ 8de716b9-ca56-42f5-aebf-39ad467f4613
tR₅-tR_lock

# ╔═╡ dadbefbf-a107-4c89-a4ef-0eb756517c1e
begin
	n₆ = itp₅(tR_lock)
	par₆ = GasChromatographyTools.stretched_program(n₆, par0)
	sol₆ = GasChromatographySimulator.solving_odesystem_r(par₆.col, par₆.prog, par₆.sub[2], par₆.opt)
	tR₆ = sol₆.u[end][1]
	n_vec₆ = sort([n_vec₅; n₆])
	tR_vec₆ = sort([tR_vec₅; tR₆])
	tR_vec₆.-tR_lock
	itp₆ = Spline1D(tR_vec₆, n_vec₆, k=3)

	scatter!(p1, [tR₆, tR₆], [n₆, n₆], label="6")
	plot!(p1, 51.0:0.001:52.0, itp₆.(51.0:0.001:52.0), xlims=(51.47,51.485), ylims=(0.4995, 0.49975), label="itp₆")
end

# ╔═╡ 313d1431-ec39-4d1e-91fd-e025efc1f5c3
tR₆-tR_lock

# ╔═╡ fab59b1e-ba32-49d3-b0f1-564db037400c
begin
	n₇ = itp₆(tR_lock)
	par₇ = GasChromatographyTools.stretched_program(n₇, par0)
	sol₇ = GasChromatographySimulator.solving_odesystem_r(par₇.col, par₇.prog, par₇.sub[2], par₇.opt)
	tR₇ = sol₇.u[end][1]
	n_vec₇ = sort([n_vec₆; n₇])
	tR_vec₇ = sort([tR_vec₆; tR₇])
	tR_vec₇.-tR_lock
	itp₇ = Spline1D(tR_vec₇, n_vec₇, k=3)

	scatter!(p1, [tR₇, tR₇], [n₇, n₇], label="7")
	plot!(p1, 51.0:0.001:52.0, itp₇.(51.0:0.001:52.0), xlims=(51.47,51.485), ylims=(0.4995, 0.49975), label="itp₇")
end

# ╔═╡ 25725568-ae9a-4ab0-8f98-87fec12c867a
begin
	n₈ = itp₇(tR_lock)
	par₈ = GasChromatographyTools.stretched_program(n₈, par0)
	sol₈ = GasChromatographySimulator.solving_odesystem_r(par₈.col, par₈.prog, par₈.sub[2], par₈.opt)
	tR₈ = sol₈.u[end][1]
	n_vec₈ = sort([n_vec₇; n₈])
	tR_vec₈ = sort([tR_vec₇; tR₈])
	tR_vec₈.-tR_lock
	itp₈ = Spline1D(tR_vec₈, n_vec₈, k=3)

	scatter!(p1, [tR₈, tR₈], [n₈, n₈], label="8")
	plot!(p1, 51.0:0.001:52.0, itp₈.(51.0:0.001:52.0), xlims=(51.47,51.485), ylims=(0.4995, 0.49975), label="itp₈")
end

# ╔═╡ 404c1d12-5217-485e-b3a6-6f024d29b544
md"""
The estimates continue to jump around tR_lock. The last estimate is further away than the one before.

Also, there seems to be some form of discontiuity of the simulation result around tR_lock, which produces the problem in the first place.
"""

# ╔═╡ 57634f00-45d5-4cf9-94f5-76319d4e5436
begin
	n₉ = itp₈(tR_lock)
	par₉ = GasChromatographyTools.stretched_program(n₉, par0)
	sol₉ = GasChromatographySimulator.solving_odesystem_r(par₉.col, par₉.prog, par₉.sub[2], par₉.opt)
	tR₉ = sol₉.u[end][1]
	n_vec₉ = sort([n_vec₈; n₉])
	tR_vec₉ = sort([tR_vec₈; tR₉])
	tR_vec₉.-tR_lock
	itp₉ = Spline1D(tR_vec₉, n_vec₉, k=3)

	scatter!(p1, [tR₉, tR₉], [n₉, n₉], label="9")
	plot!(p1, 51.0:0.001:52.0, itp₉.(51.0:0.001:52.0), xlims=(51.47,51.485), ylims=(0.4995, 0.49975), label="itp₉")
end

# ╔═╡ eba919a1-4c2c-4eaf-ba12-2c02a8d0d6fb
n_vec₉

# ╔═╡ ae803a26-c436-497a-b245-a0afec92e46f
md"""
Make simulations for a range of stretch factors ``n``.
"""

# ╔═╡ 1983258e-e84b-4f39-9ce8-0e20e78a0893
begin
	nn = [0.49954, 0.49956, 0.49958, 0.49960, 0.49962, 0.49964, 0.49966]
	ttR = Array{Float64}(undef, length(nn))
	for i=1:length(nn)
	    pars = stretched_program(nn[i], par0)
	    sols = GasChromatographySimulator.solving_odesystem_r(pars.col, pars.prog, pars.sub[2], pars.opt)
	    ttR[i] = sols.u[end][1]
	end
	p2 = plot([tR_lock, tR_lock], [0.4, 0.6], label="tR_lock")
	plot!(p2, ttR, nn, line=(2,:solid), markers=:square, xlims=(51.47,51.49), ylims=(0.4995, 0.49975))
end

# ╔═╡ 89cbd5c0-0df6-43d3-bd5c-d3d55d669f33
begin
	nnn = 0.499625.+collect(0.0000001:0.0000001:0.0000019)
	tttR = Array{Float64}(undef, length(nnn))
	for i=1:length(nnn)
	    pars = stretched_program(nnn[i], par0)
	    sols = GasChromatographySimulator.solving_odesystem_r(pars.col, pars.prog, pars.sub[2], pars.opt)
	    tttR[i] = sols.u[end][1]
	end
	plot!(p2, tttR, nnn, markers=:circle, ylims=(0.499625, 0.499627), xlims=(51.482, 51.485))
end

# ╔═╡ f541b210-3ef8-4003-b358-065e5c5949ad
par8 = stretched_program(nnn[8], par0)

# ╔═╡ 524b256b-0c1e-48bc-ac0b-2ee7fc1eab2b
par9 = stretched_program(nnn[9], par0)

# ╔═╡ 6d15c2a4-6d1e-4b89-b9e8-ecff004c4730
sol8 = GasChromatographySimulator.solving_odesystem_r(par8.col, par8.prog, par8.sub[2], par8.opt)

# ╔═╡ 665f302d-204b-47ac-90df-c5979350707c
sol9 = GasChromatographySimulator.solving_odesystem_r(par9.col, par9.prog, par9.sub[2], par9.opt)

# ╔═╡ 56095d71-6169-44b3-89d1-7ea7f1b6ddfb
sol8.destats

# ╔═╡ a51cd1dc-50bb-4444-a0f7-c0e4229b1257
sol9.destats

# ╔═╡ 225434d5-c753-4476-8f68-f5761a454852
sol8.retcode

# ╔═╡ 9a6f436e-e50e-4697-a1d1-3d2d43fc62fc
sol9.retcode

# ╔═╡ d165d269-58a9-4ad1-b5b2-6704ade60936
tR8 = sol8.u[end][1]

# ╔═╡ 61171239-d6f3-4e1b-9334-349633ee61a6
tR9 = sol9.u[end][1]

# ╔═╡ 3a759d13-d67a-4e9b-a166-dbd4eeb7b61b
tR8-tR9

# ╔═╡ 6040a096-38bd-4789-98d3-a02166e73f49
(tR8-tR9)/tR8

# ╔═╡ 08a0972a-55ff-42c9-838b-1a649afe9a46
md"""
## Decreased relative tolerance

The problem is not originated in the RT_lock algorithm but in the simulation itself. For certain settings a small change (like small strectch in the program) can result in bigger changes in retention time.
"""

# ╔═╡ 2d7da55e-2711-44a4-b8b9-bda6321a4c48
begin
	# decrease the relative tolerance
	opt_1 = GasChromatographySimulator.Options(OwrenZen5(), 1e-6, 1e-4, "inlet", true)
	par_1 = GasChromatographySimulator.Parameters(col, prog0, sub, opt_1)
	# repeat the simulation from above
	nnnn_1 = sort!(rand(0.499625:0.000000001:0.499627, 100))
	ttttR_1 = Array{Float64}(undef, length(nnnn_1))
	for i=1:length(nnnn_1)
	    pars = stretched_program(nnnn_1[i], par_1)
	    sols = GasChromatographySimulator.solving_odesystem_r(pars.col, pars.prog, pars.sub[2], pars.opt)
	    ttttR_1[i] = sols.u[end][1]
	end
	plot!(p2, ttttR_1, nnnn_1, ylims=(0.499625, 0.499627), xlims=(51.482, 51.485), markers=:v, label="reltol=1e-4")
end

# ╔═╡ 544b9b37-e047-49c1-b0d4-c58b4fde741c
Statistics.mean(ttttR_1)

# ╔═╡ b4a880cb-6138-432a-be88-a556908ab207
Statistics.std(ttttR_1)

# ╔═╡ ccd65c3a-2e6a-4540-9d82-e2bc23ef7230
begin
	# decrease the relative tolerance
	opt_2 = GasChromatographySimulator.Options(OwrenZen5(), 1e-7, 1e-4, "inlet", true)
	par_2 = GasChromatographySimulator.Parameters(col, prog0, sub, opt_2)
	# repeat the simulation from above
	nnnn_2 = sort!(rand(0.499625:0.000000001:0.499627, 100))
	ttttR_2 = Array{Float64}(undef, length(nnnn_2))
	for i=1:length(nnnn_2)
	    pars = stretched_program(nnnn_2[i], par_2)
	    sols = GasChromatographySimulator.solving_odesystem_r(pars.col, pars.prog, pars.sub[2], pars.opt)
	    ttttR_2[i] = sols.u[end][1]
	end
	plot!(p2, ttttR_2, nnnn_2, ylims=(0.499625, 0.499627), xlims=(51.482, 51.485), markers=:v, label="abstol=1e-7, reltol=1e-4")
end

# ╔═╡ 70d47e99-fcee-42bf-b502-b311bd347904
Statistics.mean(ttttR_2)

# ╔═╡ ef410f01-1bdb-4e1f-9e09-e1d309d1c546
Statistics.std(ttttR_2)

# ╔═╡ f012e638-ecdb-41f6-8cfa-16a43b503970
begin
	# decrease the relative tolerance
	opt_3 = GasChromatographySimulator.Options(OwrenZen5(), 1e-7, 1e-5, "inlet", true)
	par_3 = GasChromatographySimulator.Parameters(col, prog0, sub, opt_3)
	# repeat the simulation from above
	nnnn_3 = sort!(rand(0.499625:0.000000001:0.499627, 100))
	ttttR_3 = Array{Float64}(undef, length(nnnn_3))
	for i=1:length(nnnn_3)
	    pars = stretched_program(nnnn_3[i], par_3)
	    sols = GasChromatographySimulator.solving_odesystem_r(pars.col, pars.prog, pars.sub[2], pars.opt)
	    ttttR_3[i] = sols.u[end][1]
	end
	plot!(p2, ttttR_3, nnnn_3, ylims=(0.499625, 0.499627), xlims=(51.482, 51.485), markers=:v, label="abstol=1e-7, reltol=1e-5")
end

# ╔═╡ 012230e2-33f1-4658-8f5c-067b3df23a28
Statistics.mean(ttttR_3)

# ╔═╡ 834bfb61-b3cf-4ff9-b041-78b397b254b0
Statistics.std(ttttR_3)

# ╔═╡ 0aeb3ece-93a4-4638-bd90-1fdf13231e2c
md"""
**This is a good test inviroment for variations of abstol and reltol.**

It seems, that, if the result is 'noisy' (which is not always the case), than we have two branches of results for retention times for variing program stretch factors n.

-> new notebook in GasChromatographySimulator.jl
"""

# ╔═╡ 89fa3139-1d10-4dd3-90b1-d9c0f65745c6
begin
	n_3 = RT_locking(par_3, tR_lock, 1e-3, last_alkane; opt_itp="spline")
	par_3_s = stretched_program(n_3, par_3)
	sol_3_s = GasChromatographySimulator.solving_odesystem_r(par_3_s.col, par_3_s.prog, par_3_s.sub[2], par_3_s.opt)
	tR_3_s = sol_3_s.u[end][1]
	#scatter!(p2, [tR_1_s, tR_1_s], [n_1, n_1], ylims=(0.9699739, n_1*1.000001))
	tR_3_s-tR_lock
end

# ╔═╡ e9accaac-587f-498d-ab0c-a8064f573678
begin
	# simulation arround n_1
	n_range = (n_3-0.0001):0.000001:(n_3+0.0001)
	tR_range = Array{Float64}(undef, length(n_range))
	for i=1:length(n_range)
	    pars = stretched_program(n_range[i], par_3)
	    sols = GasChromatographySimulator.solving_odesystem_r(pars.col, pars.prog, pars.sub[2], pars.opt)
	    tR_range[i] = sols.u[end][1]
	end
	#plot!(p2, tR_range, n_range, ylims=(n_range[1], n_range[end]), xlims=(tR_range[1], tR_range[end]))
end

# ╔═╡ a0e85c4b-7332-42a6-8fa0-77a9bb0c89a8
md"""
# Conclusion
The restriction has to be more strict, than only tR\_tol<reltol. Properly tR\_tol<1000*abstol (tR\_tol=1e-3 with abstol=1e-7).
"""

# ╔═╡ f38f4866-0f1b-4a41-a370-42dcb0c9e907
function RT_locking_mod(par::GasChromatographySimulator.Parameters, tR_lock::Float64, tR_tol::Float64, solute_RT::String; opt_itp="linear")
		# estimate the factor 'n' for the temperature program to achieve the retention time 'tR_lock' for 'solute_RT' with the GC-system defined by 'par'
		if isa(par, Array)==true
			error("Select an element of the array of GC-systems.")
			elseif tR_tol<=par.opt.reltol || tR_tol<=1000*par.opt.abstol
		error("The tolerance for retention time locking tR_tol is to small. Use tR_tol > par.opt.reltol ($(par.opt.reltol)), resp. tR_tol > 1000*par.opt.abstol ($(1000*par.opt.abstol)).")
		else
			# find 'solute_RT' in the substances of 'par'
			name = Array{String}(undef, length(par.sub))
			for i=1:length(par.sub)
				name[i] = par.sub[i].name
			end
			ii = findfirst(name.==solute_RT)
			# calculate the retention time for the original (un-stretched) program
			sol₀ = GasChromatographySimulator.solving_odesystem_r(par.col, par.prog, par.sub[ii], par.opt)
			tR₀ = sol₀.u[end][1]
			# start value for the factor 'n'
			if tR₀-tR_lock<0
				n₁ = initial_n(2.0, tR_lock, ii, par)
			else
				n₁ = initial_n(0.5, tR_lock, ii, par)
			end
			# using a recursive function to estimate 'n'
			n = recur_RT_locking(n₁, [1.0], [tR₀], par, tR_lock, tR_tol, ii; opt_itp=opt_itp)
		end
		return n
	end

# ╔═╡ dfd1adb5-4681-446f-9c2e-25f89e56f6fc
RT_locking_mod(par_1, tR_lock, tR_tol, last_alkane; opt_itp="spline")

# ╔═╡ 420cff3b-f302-4751-8c6c-c11a6a45e2a5
RT_locking_mod(par_2, tR_lock, tR_tol, last_alkane; opt_itp="spline")

# ╔═╡ e37c9975-6c51-46a8-8695-c633566a841b
RT_locking_mod(par_3, tR_lock, tR_tol, last_alkane; opt_itp="spline")

# ╔═╡ abade442-bb18-406f-b9f7-6c6f6f04c515
RT_locking_mod(par_1, tR_lock, tR_tol*10, last_alkane; opt_itp="spline")

# ╔═╡ cfe9d6ec-5b22-46b5-b30f-ceed51d49644
function recur_RT_locking_rel(n::Float64, n_vec::Array{Float64,1}, tR_vec::Array{Float64,1}, par::GasChromatographySimulator.Parameters, tR_lock::Float64, tR_tol::Float64, ii::Int64; opt_itp="linear")
		# recursive function to find the factor 'n' for the temperature program to achieve the retention time 'tR_lock' for solute index 'ii' with the GC-system defined by 'par'
		# calculate the retention time with the input guess 'n'
		par₁ = stretched_program(n, par)
		if par.opt.odesys==true
			sol₁ = GasChromatographySimulator.solving_odesystem_r(par₁.col, par₁.prog, par₁.sub[ii], par₁.opt)
			tR₁ = sol₁.u[end][1]
		else
			sol₁ = GasChromatographySimulator.solving_migration(par₁.col, par₁.prog, par₁.sub[ii], par₁.opt)
			tR₁ = sol₁.u[end]
		end
		if abs(tR₁-tR_lock)/tR_lock<tR_tol
			# if retention time is less than 'tR_tol' from 'tR_lock' we found the factor 'n'
			return n
		else
			# estimate a new factor 'new_n' by linear interpolation of the factors 'n_vec' + 'n' over the corresponding retention times
			new_n_vec = sort([n_vec; n])
			new_tR_vec = sort([tR_vec; tR₁])
			if opt_itp=="spline"
				# Dierckx.jl
				if length(new_tR_vec)<4
					k = length(new_tR_vec)-1
				else
					k = 3
				end
				itp = Spline1D(new_tR_vec, new_n_vec, k=k)
			else # opt_itp=="linear"
				# Interpolations.jl
				itp = LinearInterpolation(sort([tR_vec; tR₁]), sort([n_vec; n]))
			end
			new_n = itp(tR_lock)
			#println("new_n=$(new_n), tR₁=$(tR₁)")
			# use the new factor 'new_n' and call the recursive function again
			return recur_RT_locking_rel(new_n, new_n_vec, new_tR_vec, par, tR_lock, tR_tol, ii; opt_itp=opt_itp)
		end
	end

# ╔═╡ 56a2001a-b09b-4a52-a889-c6cb599be242
function RT_locking_rel(par::GasChromatographySimulator.Parameters, tR_lock::Float64, tR_tol::Float64, solute_RT::String; opt_itp="linear")
		# estimate the factor 'n' for the temperature program to achieve the retention time 'tR_lock' for 'solute_RT' with the GC-system defined by 'par'
		if isa(par, Array)==true
			error("Select an element of the array of GC-systems.")
		elseif tR_tol<par.opt.reltol
			error("The relative tolerance for retention time locking tR_tol is to small. Use tR_tol > par.opt.reltol ($(par.opt.reltol)).")
		else
			# find 'solute_RT' in the substances of 'par'
			name = Array{String}(undef, length(par.sub))
			for i=1:length(par.sub)
				name[i] = par.sub[i].name
			end
			ii = findfirst(name.==solute_RT)
			# calculate the retention time for the original (un-stretched) program
			sol₀ = GasChromatographySimulator.solving_odesystem_r(par.col, par.prog, par.sub[ii], par.opt)
			tR₀ = sol₀.u[end][1]
			# start value for the factor 'n'
			if tR₀-tR_lock<0
				n₁ = initial_n(2.0, tR_lock, ii, par)
			else
				n₁ = initial_n(0.5, tR_lock, ii, par)
			end
			# using a recursive function to estimate 'n'
			n = recur_RT_locking_rel(n₁, [1.0], [tR₀], par, tR_lock, tR_tol, ii; opt_itp=opt_itp)
		end
		return n
	end

# ╔═╡ bb863e5c-5830-41a0-929a-001b15df7606
n_rel = RT_locking_rel(par0, tR_lock, 1e-4, last_alkane; opt_itp="spline")

# ╔═╡ 49c332ae-ae5f-4625-a693-b698ea678569
par0_rel = stretched_program(n_rel, par0)

# ╔═╡ 9356daef-2045-4391-afe8-2f67bb646dea
sol_rel = GasChromatographySimulator.solving_odesystem_r(par0_rel.col, par0_rel.prog, par0_rel.sub[2], par0_rel.opt)

# ╔═╡ 1422fb36-1933-4074-aa2c-556b9814fc42
	tR_rel = sol_rel.u[end][1]

# ╔═╡ da60b623-4c87-4798-899c-ec6e0b991cf4
tR_rel-tR_lock

# ╔═╡ 6dc84c22-042e-42ac-aab6-d048f879e228
(tR_rel-tR_lock)/tR_lock

# ╔═╡ 42a1fccd-e782-42fb-91de-18c8b33d507d
md"""
# End
"""

# ╔═╡ Cell order:
# ╠═6f9a94bd-dd28-4110-a7ca-0ed84e9c7c3f
# ╠═0b51a792-d9b8-4c29-928f-57aaf41f1c20
# ╠═93ba6afc-4e9a-11ec-08d9-81c0a9bc502e
# ╠═51ec0223-199a-4a25-8768-cd0b7e9f864d
# ╠═5e642f09-9a8f-4cca-b61f-b27c8433a2e5
# ╠═3de48bb6-8eb7-4a33-9a98-d5fe3a19f6c6
# ╠═19d1abb6-32b6-414f-bc95-55635cbaa73a
# ╠═2cad8185-4fc5-4009-98df-07a2be2133c6
# ╠═30207849-cafe-4ec6-af8d-ee7b2e2e6de0
# ╠═84e7e869-0fbf-4590-a8de-28855856661f
# ╠═c5ace7ce-cfa3-4c15-bcc1-8b66ad1c16e9
# ╠═c4d53222-03b1-4ce4-a49a-690945347432
# ╠═6acfbf9b-f8ac-4483-8c93-17920d0d9f0e
# ╠═3b31f580-efe2-4a7a-89dc-228b38f2a71e
# ╠═b6c5126a-3d41-4809-ab9b-d554262e0668
# ╠═6c5443bc-3dab-4317-99c9-1bb32b184fbc
# ╠═3b1bb187-66f9-40bd-a0a5-a3c4e1a23819
# ╠═b81b2905-6756-445b-b1be-c54298df8b3f
# ╟─eb2d9f6a-8384-427a-a036-b4f0019d8251
# ╠═734bae50-12fa-4717-8446-addb963b8673
# ╠═8de716b9-ca56-42f5-aebf-39ad467f4613
# ╠═dadbefbf-a107-4c89-a4ef-0eb756517c1e
# ╠═313d1431-ec39-4d1e-91fd-e025efc1f5c3
# ╠═fab59b1e-ba32-49d3-b0f1-564db037400c
# ╠═25725568-ae9a-4ab0-8f98-87fec12c867a
# ╠═404c1d12-5217-485e-b3a6-6f024d29b544
# ╠═57634f00-45d5-4cf9-94f5-76319d4e5436
# ╠═eba919a1-4c2c-4eaf-ba12-2c02a8d0d6fb
# ╠═ae803a26-c436-497a-b245-a0afec92e46f
# ╠═1983258e-e84b-4f39-9ce8-0e20e78a0893
# ╠═89cbd5c0-0df6-43d3-bd5c-d3d55d669f33
# ╠═f541b210-3ef8-4003-b358-065e5c5949ad
# ╠═524b256b-0c1e-48bc-ac0b-2ee7fc1eab2b
# ╠═6d15c2a4-6d1e-4b89-b9e8-ecff004c4730
# ╠═665f302d-204b-47ac-90df-c5979350707c
# ╠═56095d71-6169-44b3-89d1-7ea7f1b6ddfb
# ╠═a51cd1dc-50bb-4444-a0f7-c0e4229b1257
# ╠═225434d5-c753-4476-8f68-f5761a454852
# ╠═9a6f436e-e50e-4697-a1d1-3d2d43fc62fc
# ╠═d165d269-58a9-4ad1-b5b2-6704ade60936
# ╠═61171239-d6f3-4e1b-9334-349633ee61a6
# ╠═3a759d13-d67a-4e9b-a166-dbd4eeb7b61b
# ╠═6040a096-38bd-4789-98d3-a02166e73f49
# ╟─08a0972a-55ff-42c9-838b-1a649afe9a46
# ╠═2d7da55e-2711-44a4-b8b9-bda6321a4c48
# ╠═d3c2c2f4-1100-438d-88d6-d7502bcb9310
# ╠═544b9b37-e047-49c1-b0d4-c58b4fde741c
# ╠═b4a880cb-6138-432a-be88-a556908ab207
# ╠═ccd65c3a-2e6a-4540-9d82-e2bc23ef7230
# ╠═70d47e99-fcee-42bf-b502-b311bd347904
# ╠═ef410f01-1bdb-4e1f-9e09-e1d309d1c546
# ╠═f012e638-ecdb-41f6-8cfa-16a43b503970
# ╠═012230e2-33f1-4658-8f5c-067b3df23a28
# ╠═834bfb61-b3cf-4ff9-b041-78b397b254b0
# ╠═0aeb3ece-93a4-4638-bd90-1fdf13231e2c
# ╠═89fa3139-1d10-4dd3-90b1-d9c0f65745c6
# ╠═e9accaac-587f-498d-ab0c-a8064f573678
# ╠═a0e85c4b-7332-42a6-8fa0-77a9bb0c89a8
# ╠═f38f4866-0f1b-4a41-a370-42dcb0c9e907
# ╠═dfd1adb5-4681-446f-9c2e-25f89e56f6fc
# ╠═420cff3b-f302-4751-8c6c-c11a6a45e2a5
# ╠═e37c9975-6c51-46a8-8695-c633566a841b
# ╠═abade442-bb18-406f-b9f7-6c6f6f04c515
# ╠═56a2001a-b09b-4a52-a889-c6cb599be242
# ╠═cfe9d6ec-5b22-46b5-b30f-ceed51d49644
# ╠═bb863e5c-5830-41a0-929a-001b15df7606
# ╠═49c332ae-ae5f-4625-a693-b698ea678569
# ╠═9356daef-2045-4391-afe8-2f67bb646dea
# ╠═1422fb36-1933-4074-aa2c-556b9814fc42
# ╠═da60b623-4c87-4798-899c-ec6e0b991cf4
# ╠═6dc84c22-042e-42ac-aab6-d048f879e228
# ╠═42a1fccd-e782-42fb-91de-18c8b33d507d


text:
push!(LOAD_PATH,"../src/")

using Documenter
using Iris

makedocs(
    sitename = "Iris.jl",
    doctest = false,
    # modules = [Iris, Iris.Common, Iris.Spectral, Iris.Common.Domains],
    pages = [
            "index.md",
            "Manual" => [
                "Building a simulation" => "Common.md",
                "Spectral Analysis" => "Spectral.md",
                "Band dispersion (under construction)" => "Floquet.md",
                "Scattering" => "Scattering.md",
                "S-Matrix" => "SMatrices.md",
                "Lasing" => "Lasing.md",
                "Saturating CPA" => "SaturableCPA.md",
                "FDTD" => "TimeDomain.md",
                ],
            "Examples.md",
            "Library" => [
                "Library/Iris.md",
                "Common" => [
                    "Library/Common/Boundaries.md",
                    "Library/Common/BoundaryConditions.md",
                    "Library/Common/BoundaryLayers.md",
                    "Library/Common/Common.md",
                    "Library/Common/Curlcurls.md",
                    "Library/Common/DielectricFunctions.md",
                    "Library/Common/Dispersions.md",
                    "Library/Dispersions" => [
                        "Library/Common/Dispersions/Kerrs.md",
                        "Library/Common/Dispersions/TwoLevelSystems.md",
                        ],
                    "Library/Common/Domains.md",
                    "Library/Common/ElectricFields.md",
                    "Library/Common/Lattices.md",
                    "Library/Common/LU_Factorizations.md",
                    "Library/Common/Plotting.md",
                    "Library/Common/Points.md",
                    "Library/Common/PumpFunctions.md",
                    "Library/Common/SelfEnergies.md",
                    "Library/Common/Shapes.md",
                    "Library/Common/Simulations.md",
                    ],
                "Library/Spectral.md",
                "Library/Floquet.md",
                "Library/Scattering.md",
                "Library/SMatrices.md",
                "Library/Lasing.md",
                "Library/SaturableCPA.md",
                "Library/TimeDomain.md"
                ]
            ]
    )