text:
include("./controller.jl")

## calculate value at risk for failure scenarios
function VAR(losses, probabilities, beta)
    usorted = []
    psorted = []
    umodified = losses
    pmodified = probabilities
    while length(umodified) > 0
        s = argmin(umodified)
        push!(usorted, umodified[s])
        push!(psorted, pmodified[s])
        umodified = umodified[1:end .!= s]
        pmodified = pmodified[1:end .!= s]
    end

    total = 0
    loss = 0
    for s in 1:size(usorted, 1)
        total += psorted[s]
        loss = usorted[s]
        if total >= beta
            break
        end
    end
    return loss
end

function TrafficReAssignment(links_utilization, edges, capacity, demand, flows, T, Tf, k, TunnelBw, FlowBw, scenarios, restored, algorithm, verbose; progress=false)
    nedges = length(edges)
    nflows = length(flows)
    ntunnels = length(T)
    nscenarios = length(scenarios)
    RouterPorts = deepcopy(links_utilization)
    arrow_algorithms = ["ARROW", "ARROW_NAIVE", "ARROW_BIN"]

    scenario_lost = []
    scenario_restored = []
    if verbose
        if in(algorithm, arrow_algorithms)
            for q in 1:nscenarios
                scenario_capacity = capacity .* scenarios[q]
                push!(scenario_lost, sum(capacity)-sum(scenario_capacity))
                push!(scenario_restored, sum(restored[q])-sum(scenario_capacity))
            end
            println("per scenario lost capacity ", scenario_lost)
            println("per scenario sum restored ", scenario_restored)
            println("$(algorithm) Scenario restoration ratio", scenario_restored ./ scenario_lost)
        else
            for q in 1:nscenarios
                scenario_capacity = capacity .* scenarios[q]
                push!(scenario_lost, sum(capacity)-sum(scenario_capacity))
                push!(scenario_restored, sum(restored[q]))
            end
            if sum(scenario_restored) > 0
                println("$(algorithm) selected per scenario sum restored ", scenario_restored)
                println("$(algorithm) Scenario restoration ratio", scenario_restored ./ scenario_lost)
                println("$(algorithm) sum scenario restored ", sum(scenario_restored))
                printstyled("Error! $(algorithm) has restorable capacity!\n", color=:red)
            end
        end
    end
    if verbose println("reassign demand: ", sum(demand), " ", demand) end
    if verbose println("reassign capacity: ", sum(capacity), " ", capacity) end
    if verbose println("reassign restored: ", sum(restored), " ", restored) end

    #CREATE TUNNEL SCENARIO MATRIX
    X  = ones(nscenarios,ntunnels)
    for s in 1:nscenarios
        for t in 1:ntunnels
            if size(T[t],1) == 0
                X[s,t] = 0
            else
                for e in 1:nedges
                    if scenarios[s][e] == 0 && restored[s][e] == 0  # both failed and non-restorable
                        back_edge = findfirst(x -> x == (edges[e][2],edges[e][1],edges[e][3]), edges)
                        if in(e, T[t]) || in(back_edge, T[t])
                            X[s,t] = 0
                        end
                    end
                end
            end
        end
    end

    #CREATE TUNNEL EDGE MATRIX
    L = zeros(ntunnels, nedges)
    for t in 1:ntunnels
        for e in 1:nedges
            if in(e, T[t])
                L[t,e] = 1
            end
        end
    end

    routed = zeros(nscenarios, nflows, k)
    allowed = zeros(nscenarios)
    scenarioAffectedFlow = zeros(nscenarios)
    scenarioCongestionLoss = zeros(nscenarios)

    # SCENARIO LOSS PER FLOW
    for s in 1:nscenarios
        if verbose printstyled("== Process scenario $(s)\n", color=:yellow) end
        scenario_network_capacity = 0
        scenario_tunnel_bw = deepcopy(TunnelBw)  ## every failure scenario should start with the initial TE splitratio results
        scenario_flow_bw = deepcopy(FlowBw)
        affected_flows = zeros(nflows)

        ## healthy state should respect tunnel reservations (FFC, TEAVAR, ARROW), overflow traffic should be put on other 0 allocated tunnels
        ## in this case, the tunnel allocation will be modified because we are accepting more traffic than formulation output
        for f in 1:nflows
            if occursin("TEAVAR", algorithm)  || occursin("FFC", algorithm) || occursin("ARROW", algorithm) || occursin("ARROW_NAIVE", algorithm)  # reservation based TE
                if scenario_flow_bw[f] < demand[f]  # assign overflow traffic to 0 unallocated tunnels
                    overflow = demand[f] - scenario_flow_bw[f]
                    zero_unallocated_tunnel = zeros(size(Tf[f],1))
                    for t in 1:size(Tf[f],1)
                        if scenario_tunnel_bw[f,t] == 0
                            zero_unallocated_tunnel[t] = 1
                        end
                    end
                    if sum(zero_unallocated_tunnel) > 0  # there exist 0 unallocated tunnel
                        if verbose println("Flow $(f)'s overflow $(overflow) load balancing to $(sum(zero_unallocated_tunnel)) unallocated tunnels") end
                        for t in 1:size(Tf[f],1)
                            if zero_unallocated_tunnel[t] == 1
                                scenario_tunnel_bw[f,t] = round(overflow/sum(zero_unallocated_tunnel), digits=16)  # actual allocation without reservation
                            end
                        end
                    end
                end
            end
        end

        if verbose println("scenario_tunnel_bw ", scenario_tunnel_bw) end

        ## calculate demand loss due to local proportional routing under failure scenarios
        for f in 1:nflows
            ## for reservation based TE, use reservation first, if all reservation fails, use other 0 allocated tunnels with local proportional routing for traffic
            totalup = 0
            for t in 1:size(Tf[f],1)
                totalup += round(scenario_tunnel_bw[f,t] * X[s,Tf[f][t]], digits=16)  # split ratio sum of residual (or restorable) tunnels
            end
            if totalup == 0
                scenario_tunnel_bw[f,:] = scenario_tunnel_bw[f,:] .+ 1
                if verbose println("ECMP for local proportional routing because of no/failed reserved tunnels of flow $(f)") end
                for t in 1:size(Tf[f],1)
                    totalup += round(scenario_tunnel_bw[f,t] * X[s,Tf[f][t]], digits=16)  # sum split ratios of survived tunnels
                end
            end
            for t in 1:size(Tf[f],1)
                if totalup != 0
                    routed[s,f,t] = round(scenario_tunnel_bw[f,t] * X[s,Tf[f][t]] / totalup, digits=16) * round(demand[f], digits=16)  # local proportional routing
                else
                    if verbose println("flow $(f)'s original $(t) tunnels all fail (original bw $(scenario_tunnel_bw[f,t]), $(TunnelBw[f,t]))") end
                end
            end
            if round((sum(routed[s,f,:]) - demand[f])/demand[f],  digits=6) < 0
                affected_flows[f] = 1
            end
        end
        allowed[s] = round(Float64(sum(routed[s,:,:])),  digits=16)  # expected throughput via local proportional routing (congestion included)
        if verbose println("$(algorithm) allowed: $(allowed), demand: $(demand)") end

        ## calculate demand loss due to link congestion after local proportional routing
        congestion_loss = 0
        for e in 1:nedges
            edge_utilization = 0
            for f in 1:nflows  # edge utilization after local proportional routing
                edge_utilization += round(Float64(sum(routed[s,f,t] * L[Tf[f][t],e] * X[s,Tf[f][t]] for t in 1:size(Tf[f],1))),  digits=16)
            end
            # println("Edge: ", edges[e])
            # println(max(0, edge_utilization - capacity[e]))
            if restored[s][e] > 0 && scenarios[s][e] == 0
                per_link_current_capacity = restored[s][e]
            elseif restored[s][e] > 0 && scenarios[s][e] > 0
                per_link_current_capacity = capacity[e]
            elseif restored[s][e] == 0 && scenarios[s][e] > 0
                per_link_current_capacity = capacity[e]
            elseif restored[s][e] == 0 && scenarios[s][e] == 0
                per_link_current_capacity = 0
            end
            scenario_network_capacity += per_link_current_capacity

            # RouterPorts[e] = max(RouterPorts[e], max(0, min(edge_utilization, per_link_current_capacity)-restored[s][e]*scenarios[s][e]))
            RouterPorts[e] = max(RouterPorts[e], min(edge_utilization, per_link_current_capacity))
            congestion_loss += max(0, round(Float64(edge_utilization - per_link_current_capacity), digits=6))

            # if there is a link congested, then all tunnels (flows) on that link will be affected (uniformly)
            if round(Float64(edge_utilization - per_link_current_capacity), digits=6) > 0
                for f in 1:nflows
                    flow_affect = 0
                    for t in 1:size(Tf[f],1)
                        if L[Tf[f][t],e] > 0 && scenario_tunnel_bw[f,t] > 0  # all non-zero tunnel after local proportional routing running on the failed link
                            flow_affect = 1  # t is the global index of all tunnels
                        end
                    end
                    affected_flows[f] = max(flow_affect, affected_flows[f])
                end
                if verbose printstyled("$(algorithm) network $(s) scenario's $(e) link: restored-$(restored[s][e])/capacity-$(capacity[e])/fail-$(scenarios[s][e]) edge_utilization $(edge_utilization) per_link_current_capacity $(per_link_current_capacity) affected flow $(sum(affected_flows)) \n", color=:red) end
            end
        end
        if verbose println("scenario_network_capacity ", scenario_network_capacity) end
        scenarioAffectedFlow[s] = sum(affected_flows)  # number of affected flows
        scenarioCongestionLoss[s] += congestion_loss   # amount of affected traffic
    end

    if scenarioCongestionLoss[1] > 0 && verbose
        printstyled("healthy $(algorithm) network $(scenarioAffectedFlow[1]) flows has congestion loss $(scenarioCongestionLoss[1]): demand too much!\n", color=:yellow)
    end
    if verbose println("allowed ", allowed) end

    umax = []
    for l in 1:nscenarios
        if sum(demand) == 0
            push!(umax, 1)  # if demand overscaling, no bw accessible, loss is equivalent 100%
            scenarioAffectedFlow[l] = nflows  # all flow loss
        else
            push!(umax, 1-round((allowed[l] - scenarioCongestionLoss[l])/sum(demand), digits=6))
        end
    end

    umax_zero_indices = findall(x -> x == 0, umax)
    flows_zero_indices = findall(x -> x == 0, scenarioAffectedFlow)

    if umax_zero_indices != flows_zero_indices && verbose
        printstyled("Losses don't match on zero indices\n", color=:red)
        printstyled("$(algorithm) allowed  $(allowed)\n", color=:red)
        printstyled("$(algorithm) scenarioCongestionLoss $(scenarioCongestionLoss)\n", color=:red)
        printstyled("$(algorithm) Loss ratio $(umax)\n", color=:red)
        printstyled("$(algorithm) scenarioAffectedFlow $(scenarioAffectedFlow)\n", color=:red)
    end

    return umax, scenarioAffectedFlow, RouterPorts
end

## percentage of scenarios that sustain 100% throughput
function ScenarioAvailability(losses, probabilities, sla, conditional)
    usorted = []
    psorted = []
    if conditional
        umodified = deepcopy(losses[2:end])
        pmodified = deepcopy(probabilities[2:end])
        pmodified = pmodified / (1 - probabilities[1])
    else
        umodified = deepcopy(losses)
        pmodified = deepcopy(probabilities)
    end

    while length(umodified) > 0
        s = argmin(umodified)
        push!(usorted, umodified[s])
        push!(psorted, pmodified[s])
        umodified = umodified[1:end .!= s]
        pmodified = pmodified[1:end .!= s]
    end

    total = 0
    loss = 0
    for s in 1:length(usorted)
        loss = usorted[s]
        if loss > sla
            break
        end
        total += psorted[s]  # all-or-nothing scenario level availability
    end

    return total
end

## percentage of flows weighted by scenarios that sustain 100% throughput
function FlowAvailability(losses, affectflows, flows, probabilities, sla, conditional)
    usorted = []
    psorted = []
    if conditional
        umodified = deepcopy(losses[2:end])
        pmodified = deepcopy(probabilities[2:end])
        pmodified = pmodified / (1 - probabilities[1])
    else
        umodified = deepcopy(losses)
        pmodified = deepcopy(probabilities)
    end

    while length(umodified) > 0
        s = argmin(umodified)
        push!(usorted, umodified[s])
        push!(psorted, pmodified[s])
        umodified = umodified[1:end .!= s]
        pmodified = pmodified[1:end .!= s]
    end
    total = 0
    loss = 0
    for s in 1:length(usorted)
        loss = usorted[s]
        if loss > sla  # there is loss
            total += (psorted[s]) * (1 - (affectflows[s] / length(flows)))  # flow level availability
        else
            total += (psorted[s])
        end
    end
    return total
end

# percentage of un-affected traffic weighted by sccenario that sustain 100% throughput
function BandwidthAvailability(losses, probabilities, sla, conditional)
    usorted = []
    psorted = []
    if conditional
        umodified = deepcopy(losses[2:end])
        pmodified = deepcopy(probabilities[2:end])
        pmodified = pmodified / (1 - probabilities[1])
    else
        umodified = deepcopy(losses)
        pmodified = deepcopy(probabilities)
    end

    while length(umodified) > 0
        s = argmin(umodified)
        push!(usorted, umodified[s])
        push!(psorted, pmodified[s])
        umodified = umodified[1:end .!= s]
        pmodified = pmodified[1:end .!= s]
    end

    total = 0
    loss = 0
    for s in 1:length(usorted)
        loss = usorted[s]
        if loss > sla  # there is loss
            total += (psorted[s]) * (1 - loss)  # bandwidth level availability
        else
            total += (psorted[s])
        end
    end

    return total
end

## compute link utilization in units of bw
function computeLinksUtilization(edges, capacity, demand, flows, T, Tf, k, splittingratios_original; progress=false)
    nedges = length(edges)
    nflows = length(flows)
    ntunnels = length(T)

    #CREATE TUNNEL EDGE MATRIX
    L = zeros(ntunnels, nedges)
    for t in 1:ntunnels
        for e in 1:nedges
            if in(e, T[t])
                L[t,e] = 1
            end
        end
    end

    splittingratios = deepcopy(splittingratios_original)

    # compute how reassignment is done under no failure scenario
    routed = zeros(nflows, ntunnels)
    for f in 1:nflows
        totalup = 0
        for t in 1:size(Tf[f],1)
            totalup += Float64(splittingratios[f,t])
        end
        ## if the only tunnel fails, other original 0 allocated tunnels should be used
        # if totalup == 0
        #     splittingratios[f,:] = splittingratios[f,:] .+ 1
        #     for t in 1:size(Tf[f],1)
        #         totalup += Float64(splittingratios[f,t])
        #     end
        # end

        for t in 1:size(Tf[f],1)
            if totalup != 0
                routed[f,t] = Float64(splittingratios[f,t] / totalup * demand[f])
            else
                routed[f,t] = 0  # no traffic routed for this flow
            end
        end
    end

    links_utilization = [0.0 for e in 1:nedges]
    for e in 1:nedges
        e_utilization = sum(routed[f,t] * L[Tf[f][t],e] for f in 1:nflows for t in 1:size(Tf[f],1))
        links_utilization[e] = e_utilization
    end

    return links_utilization
end

max_stars_repo_name: hipersys-team/arrow

max_issues_repo_name: arrow-bow/arrow

max_forks_repo_name: arrow-bow/arrow


text:
"""

    fdt(λ::Real, ν::Real, dt::Real)

Compute the variance that satisfies the fluctuation-dissipation theorem.

# External links

* [Fluctuation-dissipation theorem on Wikipedia](https://en.wikipedia.org/wiki/Fluctuation-dissipation_theorem)

# Examples
```jldoctest
julia> fdt(1.0, 1.0, 0.1)
1.378404875209022

```

"""
function fdt(λ::Real, ν::Real, dt::Real)
    return sqrt(ν * (2*λ*dt - λ^2*dt^2) / dt)
end

"""

    CAM1D(d::SGS, n::Integer; dt::Real=1/24, λ::Real=1.0)

Create a timeseries using the one-dimensional CAM noise model of Sardeshmukh and
Sura (2009) for a particular SGS distribution.

See also:

>Sardeshmukh, P. D., and P. Sura, 2009: Reconciling Non-Gaussian Climate Statistics with Linear Dynamics. *Journal of Climate*, **22**, 1193–1207, [https://doi.org/10.1175/2008JCLI2358.1](https://doi.org/10.1175/2008JCLI2358.1).
>
>Sardeshmukh, P. D., G. P. Compo, and C. Penland, 2015: Need for Caution in Interpreting Extreme Weather Statistics. *Journal of Climate*, **28**, 9166–9187, [https://doi.org/10.1175/JCLI-D-15-0020.1](https://doi.org/10.1175/JCLI-D-15-0020.1).

# Examples
```jldoctest
julia> d = fit(SGS, 1, 1, 5)
SGS{Float64}(E=0.6236095644623235, b=1.1709106481844573, g=0.48997894350611143)

julia> x = CAM1D(d, 1000, seed=42)
retcode: Success
Interpolation: 1st order linear
t: 1001-element Array{Float64,1}:
  0.0
  0.041666666666666664
  0.08333333333333333
  0.125
  0.16666666666666666
  0.20833333333333331
  0.24999999999999997
  0.29166666666666663
  0.3333333333333333
  0.375
  ⋮
 41.333333333333165
 41.37499999999983
 41.416666666666494
 41.45833333333316
 41.49999999999982
 41.54166666666649
 41.58333333333315
 41.624999999999815
 41.625
u: 1001-element Array{Array{Float64,1},1}:
 [0.0]
 [-0.05538507236046931]
 [0.3759078513765659]
 [0.31482891888314446]
 [0.2586995978771358]
 [0.0063144005340343146]
 [0.0819076515551697]
 [0.07880240481554046]
 [-0.09903540075387843]
 [-0.0867782526753085]
 ⋮
 [-0.048583800782025996]
 [0.30770620643597935]
 [0.6914165604620269]
 [1.2048828846162343]
 [0.6995798768325538]
 [0.5417553247799438]
 [0.593919650292867]
 [0.6785936302181321]
 [0.6785940957553152]

```

"""
function CAM1D(d::SGS, n::Integer; dt::Real=1/24, λ::Real=1.0, seed::Real=-1)
    # TODO: What are the best values for the default arguments?

    tspan = (0.0, n*dt-dt)      # time domain
    u0 = [0.0]                  # initial value
    p = [d.E, d.b, d.g, λ]      # parameter array

    prob = SDEProblem(CAM1D_dt, CAM1D_dW, u0, tspan, p, noise_rate_prototype=zeros(1,2))  #,reltol=1e-8,abstol=1e-8)

    # alg = RKMil(interpretation=:Stratonovich)
    alg = EulerHeun()

    if seed > -1
        sol = solve(prob, alg, dt=dt, seed=seed, abstol=1e-10, reltol=1e-10)
    else
        sol = solve(prob, alg, dt=dt, abstol=1e-10, reltol=1e-10) #, seed=Random.seed!())
    end

    return sol

end

function CAM1D(u0::AbstractArray, tspan::Tuple, p::AbstractArray, dt::Real; seed::Real=-1)

    prob = SDEProblem(CAM1D_dt, CAM1D_dW, u0, tspan, p, noise_rate_prototype=zeros(1,2))  #,reltol=1e-8,abstol=1e-8)

    if seed > -1
        sol = solve(prob, EulerHeun(), dt=dt, seed=seed)
    else
        sol = solve(prob, EulerHeun(), dt=dt) #, seed=Random.seed!())
    end

    return sol

end

"""

    Hasselmann1D(n::Integer; dt::Real=1/24, λ::Real=1.0, seed::Real=-1)

Create a timeseries using the one-dimensional noise model of Hasselmann (1976).

See also:

>Hasselmann, K., 1976: Stochastic climate models Part I. Theory. *Tellus*, **28**, 473–485, [https://doi.org/10.1111/j.2153-3490.1976.tb00696.x](https://doi.org/10.1111/j.2153-3490.1976.tb00696.x).

# Examples
```jldoctest
julia> Hasselmann1D(1000, seed=42)
retcode: Success
Interpolation: 1st order linear
t: 1001-element Array{Float64,1}:
  0.0
  0.041666666666666664
  0.08333333333333333
  0.125
  0.16666666666666666
  0.20833333333333331
  0.24999999999999997
  0.29166666666666663
  0.3333333333333333
  0.375
  ⋮
 41.333333333333165
 41.37499999999983
 41.416666666666494
 41.45833333333316
 41.49999999999982
 41.54166666666649
 41.58333333333315
 41.624999999999815
 41.625
u: 1001-element Array{Array{Float64,1},1}:
 [0.0]
 [-0.13005535542166108]
 [0.04828669603525988]
 [0.658691166050971]
 [0.4766431404349038]
 [0.2527905702593073]
 [0.5063085613562422]
 [0.4424101651546515]
 [0.42457583766050366]
 [0.17019337828909467]
 ⋮
 [0.04255470162528974]
 [-0.009396010196968833]
 [-0.4700188956712661]
 [-0.6640224731670598]
 [-0.6412697125734018]
 [-0.8226707386234939]
 [-0.9240097641730844]
 [-0.3550868399241892]
 [-0.35508653697280324]

```

"""
function Hasselmann1D(n::Integer; dt::Real=1/24, λ::Real=1.0, seed::Real=-1)
    # TODO: What are the best values for the default arguments?

    ν = 1.0                     # ASSUMED: variance of the white noise is 1.0
    σ = fdt(λ, ν, dt)           # scale the variance for the Hasselman case.
    p = [λ, σ]                  # DiffEq.jl parameter array
    tspan = (0.0, n*dt-dt)      # time domain
    u0 = [0.0]                  # initial value

    prob = SDEProblem(Hasselmann1D_dt, Hasselmann1D_dW, u0, tspan, p)

    if seed > -1
        sol = solve(prob, EM(), dt=dt, seed=seed)
    else
        sol = solve(prob, EM(), dt=dt)
    end

    return sol

end

function Hasselmann1D(u0::AbstractArray, tspan::Tuple, p::AbstractArray, dt::Real; seed::Real=-1)

    prob = SDEProblem(Hasselmann1D_dt, Hasselmann1D_dW, u0, tspan, p)

    if seed > -1
        sol = solve(prob, EM(), dt=dt, seed=seed)
    else
        sol = solve(prob, EM(), dt=dt)
    end

    return sol

end

function CAM1D_dt(du::AbstractArray, u::AbstractArray, p::AbstractArray, t::Real)

    # Unpack the parameters.
    E = p[1]
    b = p[2]
    g = p[3]
    λ = p[4]

    # Unpack the variables.
    x = u[1]

    du[1] = -((1.0 + 0.5 * E^2) * x + 0.5 * E * g) * λ

    return du

end

function CAM1D_dW(du::AbstractArray, u::AbstractArray, p::AbstractArray, t::Real)

    # Unpack the parameters.
    E = p[1]
    b = p[2]
    g = p[3]
    λ = p[4]

    # Unpack the variables.
    x = u[1]

    # Calculate two temporally uncorrelated standard normal variables.
    # η1 = randn()
    # η2 = randn()

    # du[1] = (b * η1 + (E * x + g) * η2) * sqrt(λ)
    du[1,1] = b * sqrt(λ)
    du[1,2] = (E * x + g) * sqrt(λ)

    return du

end

"""

    CAM1D_dx(d::SGS, x::Real, λ::Real, η1::Float64, η2::Float64, dt::Real)

Compute the change in the stochastically perturbed damped linear Markov process
``x(t)`` in a small time interval, ``dt``, for the CAM noise model.

"""
function CAM1D_dx(d::SGS, x::Real, λ::Real, η1::Float64, η2::Float64, dt::Real)
    (E, b, g) = params(d)

    # Construct the change, dx, of the stochastically perturbed damped linear
    # Markov process x(t) in a small time interval dt.
    term1 = -((1.0 + 0.5 * E^2) * x + 0.5 * E * g) * λ * dt
    term2 = ((b * η1) + ((E * x + g) * η2)) * sqrt(λ * dt)

    return term1 + term2

end

function Hasselmann1D_dt(du::AbstractArray, u::AbstractArray, p::AbstractArray, t::Real)

    # Unpack the parameters, variables.
    λ = p[1]
    σ = p[2]
    x = u[1]

    du[1] = -(λ * x)

    return du

end

function Hasselmann1D_dW(du::AbstractArray, u::AbstractArray, p::AbstractArray, t::Real)

    # Unpack the parameters, variables.
    λ = p[1]
    σ = p[2]
    x = u[1]

    du[1] = σ

    return du

end

"""

    Hasselmann1D_dx(x::Real, λ::Real, σ::Real, η::Float64, dt::Real)

Compute the change in the stochastically perturbed damped linear Markov process
``x(t)`` in a small time interval, ``dt``, for the Hasselmann model.

"""
function Hasselmann1D_dx(x::Real, λ::Real, σ::Real, η::Float64, dt::Real)

    return (-λ * x * dt) + (σ * η * sqrt(dt))

end

function Euler(d::Normal, u0::AbstractArray, tspan::Tuple, p::AbstractArray, dt::Real)

    # Unpack the parameters, variables.
    λ = p[1]
    σ = p[2]
    n = Int(round((tspan[2] - tspan[1]) / dt))     # number of iterations

    η = randn(n)
    u = Array{Float64}(undef, n)
    t = Array{Float64}(undef, n)

    u[1] = u0[1]
    t[1] = tspan[1]

    for i=2:n
        u[i] = u[i-1] + Hasselmann1D_dx(u[i-1], λ, σ, η[i], dt)
        t[i] = t[i-1] + dt
    end

    return (u, t)

end

function Euler(d::SGS, u0::AbstractArray, tspan::Tuple, p::AbstractArray, dt::Real)

    # Unpack the parameters.
    E = p[1]
    b = p[2]
    g = p[3]
    λ = p[4]
    n = Int(round((tspan[2] - tspan[1]) / dt))       # number of iterations

    # Calculate two temporally uncorrelated standard normal variables.
    η1 = randn(n)
    η2 = randn(n)
    u = Array{Float64}(undef, n)
    t = Array{Float64}(undef, n)

    u[1] = u0[1]
    t[1] = tspan[1]

    for i=2:n
        u[i] = u[i-1] + CAM1D_dx(d, u[i-1], λ, η1[i], η2[i], dt)
        t[i] = t[i-1] + dt
    end

    return (u, t)

end

function RK4(d::Normal, u0::AbstractArray, tspan::Tuple, p::AbstractArray, dt::Real)

    # Unpack the parameters, variables.
    λ = p[1]
    σ = p[2]
    n = Int(round((tspan[2] - tspan[1]) / dt))       # number of iterations

    η = randn(n)
    u = Array{Float64}(undef, n)
    t = Array{Float64}(undef, n)

    u[1] = u0[1]
    t[1] = tspan[1]

    for i=2:n
        k11 = Hasselmann1D_dx(u[i-1], λ, σ, η[i], dt)
        k21 = Hasselmann1D_dx(u[i-1] + (k11 / 2.0), λ, σ, η[i], dt)
        k31 = Hasselmann1D_dx(u[i-1] + (k21 / 2.0), λ, σ, η[i], dt)
        k41 = Hasselmann1D_dx(u[i-1] + k31, λ, 1.0, η[i], dt)

        u[i] = u[i-1] + (k11 + 2.0*k21 + 2.0*k31 + k41) / 6.0
        t[i] = t[i-1] + dt
    end

    return (u, t)

end

function RK4(d::SGS, u0::AbstractArray, tspan::Tuple, p::AbstractArray, dt::Real)

    # Unpack the parameters.
    E = p[1]
    b = p[2]
    g = p[3]
    λ = p[4]
    n = Int(round((tspan[2] - tspan[1]) / dt))       # number of iterations

    # Calculate two temporally uncorrelated standard normal variables.
    η1 = randn(n)
    η2 = randn(n)
    u = Array{Float64}(undef, n)
    t = Array{Float64}(undef, n)

    u[1] = u0[1]
    t[1] = tspan[1]

    for i=2:n
        # Do not mulitply kij by dt?
        k11 = CAM1D_dx(d, u[i-1], λ, η1[i], η2[i], dt)
        k21 = CAM1D_dx(d, u[i-1] + (k11 / 2.0), λ, η1[i], η2[i], dt)
        k31 = CAM1D_dx(d, u[i-1] + (k21 / 2.0), λ, η1[i], η2[i], dt)
        k41 = CAM1D_dx(d, u[i-1] + k31, λ, η1[i], η2[i], dt)

        u[i] = u[i-1] + (k11 + 2.0*k21 + 2.0*k31 + k41) / 6.0
        t[i] = t[i-1] + dt
    end

    return (u, t)

end

ext: jl


text:
# COVID 19 Global Deaths
# Building a situation graph regarding global deaths due to covid19
# Achieving this graph a time series csv file cab be downloaded from
# John Hopkins' github repository.

# Import all needed packages using the 'using' keyword followed by the package
# name i.e. using Plots.

# IMPORTTANT. To run these packages or any line of code click on the package name
# or line of code and press Shift + Enter once. This will excecute the code.

using Plots
using CSV
using Dates
using Interact
using DataFrames
using PlotThemes
using ORCA

# The COVID-19 Data Repository contains the up to date
# time_series_covid19_deaths_global.csv file which should do the job
# a. Head to the github weblink: https://github.com/CSSEGISandData/COVID-19
# b. Look for the directory named csse_covid_19_data and click once to open it
# c. Then look for the csse_covid_19_time_series directory and click once to open it
# d. The time_series_covid19_deaths_global.csv is the file you will need to work with.
# e. On the right-hand side click on the option Raw
# d. Finallym in the web browser's search bar highlight the webaddress and copy it

# Respecting Julia's naming convention; All variable names, parameters, and
# function names will be in lowercase. Camel casing names should only be used
# if the variable name is be hard to read if its in lowercase.

# Creat a variable named url or any name relevant to your project and
# assign this new varialbe name with the webaddress copied in the previous steps
url = "https://raw.githubusercontent.com/CSSEGISandData/COVID-19/81edec24de608560c4a67d497248a7f658a1fb0b/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv"

# Download this csv_covid_19_data file using the download()
# The download function takes two parameter;
# A variable name and a name for the downloaded file
# Note that this csv file will be downloaded in the working directory
# This is where jl extensio file is located
download(url, "csv_covid_19_global_death.csv")

# Lets find out if the file has been downloaded using the function readdir()
# You should see the newly named csv file -> csv_covid_19_global_death.csv
readdir()

# We can assign this downloaded csv file with a variable name which can
# make it more convenient to work with. We can assign it the name data.
# data = CSV.read("csv_covid_19_global_death.csv")
data  = DataFrame(CSV.File("csv_covid_19_global_death.csv"))
data

# After viewing the data, we can rename certain columns
# This may provide a good way to manipulate data with certain columns
globaldata = rename!(data, 1=> "province", 2=>"country")

# After renaming the desired colum(s) we can view the csv file again
globaldata
# We may also need the original data when things goes wrong.
# It is always good practice to make a copy of the original data before
# renaming colum(s). In this instance there is no need to make a copy
# because the original data is readily available on the internet

# The purpose of this project is to demonstrate global death due to the coronavirus
# So lets start extracting data from columns we might need to present our graph with
# We will extract all countries and territories from the column we renamed to country
# and put it in a new variable named countries
countries = globaldata[:,2] # colon means read all rows
                    # and the 2 means the second column
                    # Take note unlike in other languages where
                    # indices start at 0, in Julia indexing start at 1

# Now lets collect these countries and put them gain in the varialbe allcountries
# using the collect();
allcountries = collect(globaldata[:, 2])

# We need the date of death (sadly) and we will use the names() assign it to
# the variable name dates
# Read the Julia documentation about the names()
dates = names(globaldata)

# If we want to see that death rate of a specific conuntry
# we can make use of the function startswith(), lets try South Africa
startswith("South", "S")
# The return type startswith() is a boolean.
# We can take it a step futher by looking up the column country
# filtering out all possible entries that might return a result
# starting with South and assign this to with a new variable name
sa = [startswith(country, "South") for country in allcountries]

globaldata[sa, :] #We will look for all possible fields that starts with South

# Have you noticed in the previous steps that some sort of table is returned
# showing a table? This is achieved with the DataFrames packages, pretty neat!

allcountries .== "South Africa"
# period before == is "broadcasting": apply operation to each element of a vector

# We will now look for the entry South Africa in allcountries and assign it to
# a sadata
sadata = [startswith(country, "South Africa") for country in allcountries]

# findfirst() is able to find South Africa
sadatarow = findfirst(allcountries .== "South Africa")

# Once found we can assign this to southafricadata which is the start to plotting
# our graph
southafricadata = globaldata[sadatarow, :]

# Starting from the 5th column; this is where available data begins
saplot = convert(Vector, southafricadata[5:end])

# After we have the countries and days of death we need to collect
# this data by calling the String() with the confirmed days
# We need to apply the String() to each element
# we will call the days:
confirmeddate = String.(names(globaldata))[5:end]
                                #I noticed the days start from
                                #the 5th colum till the last
#Let see what date we can retrieve  on the first row
confirmeddate[1]

# Notice that the date are formatted as mm/dd/yy
formatdate = Dates.DateFormat("m/d/y")

# We can parse these dates: meaning making the dates which are strings
# into numbers.
# The first argument is the Date() from the Dates package, then we need
# to tell pars() which object we want the Date() to apply to and the format
parse(Date, confirmeddate[1], formatdate) +Year(2000)

# You can assign these dates into a new varialbe name which will form the x-axis
# .+Year() is a way to fill in the year for centuries
deceaseddays = parse.(Date, confirmeddate, formatdate) .+Year(2000)

# Next, we will uses theh GR framwork
# gr is a universal framework for cross-platform visualization applications.
# It represents a 2D x,y graph
gr()

#We will use the keyword plot()and as argument our final plotting variable name
plot(saplot)
plotlyjs()
# We will provide attrributes to our plotting here.
# Read more about these atrributes
defaultplot = plot(deceaseddays, legend=:bottomleft, label="time-line",
     saplot,
     xticks=deceaseddays[1:14:end],
     xrotation=45,
     leg=:topleft,
     grid=false,
     #size=(1200,800)
     )

    scatter!(deceaseddays,
             saplot,
            linewidths=1,
            label="death",
            marker=:o,
            markersize=3,
            markeralpha = 7.0,
            markercolor="brown",
            hover=saplot)
theme(:juno)
xlabel!("Confirmed death every 14 days")
ylabel!("Total deaths")
title!("COVID-19 death")


text:
using VisualRegressionTests
# using ExamplePlots

if isinteractive()
    @eval Main import Gtk
end

# import DataFrames, RDatasets

# don't let pyplot use a gui... it'll crash
# note: Agg will set gui -> :none in PyPlot
# ENV["MPLBACKEND"] = "Agg"
# try
#   @eval import PyPlot
#   @info("Matplotlib version: $(PyPlot.matplotlib[:__version__])")
# end

using Plots
# using StatsPlots
import PlotReferenceImages
using Random
using Test

default(size=(500,300))

# TODO: use julia's Condition type and the wait() and notify() functions to initialize a Window, then wait() on a condition that
#       is referenced in a button press callback (the button clicked callback will call notify() on that condition)

import Plots._current_plots_version

function image_comparison_tests(pkg::Symbol, idx::Int; debug = false, popup = isinteractive(), sigma = [1,1], tol = 1e-2)
    Plots._debugMode.on = debug
    example = Plots._examples[idx]
    @info("Testing plot: $pkg:$idx:$(example.header)")
    backend(pkg)
    backend()

    # ensure consistent results
    Random.seed!(1234)

    # reference image directory setup
    # refdir = joinpath(Pkg.dir("ExamplePlots"), "test", "refimg", string(pkg))
    refdir = joinpath(dirname(pathof(PlotReferenceImages)), "..", "Plots", string(pkg))
    fn = "ref$idx.png"

    # firgure out version info
    vns = filter(x->x[1] != '.', readdir(refdir))
    versions = sort(VersionNumber.(vns), rev = true)
    versions = filter(v -> v <= _current_plots_version, versions)
    # @show refdir fn versions

    newdir = joinpath(refdir, string(_current_plots_version))
    newfn = joinpath(newdir, fn)

    # figure out which reference file we should compare to, by finding the highest versioned file
    reffn = nothing
    for v in versions
        tmpfn = joinpath(refdir, string(v), fn)
        if isfile(tmpfn)
            reffn = tmpfn
            break
        end
    end

    # now we have the fn (if any)... do the comparison
    # @show reffn
    if reffn == nothing
        reffn = newfn
    end
    # @show reffn
    # return

    # test function
    func = (fn, idx) -> begin
        map(eval, example.exprs)
        png(fn)
    end

    # try
    #     run(`mkdir -p $newdir`)
    # catch err
    #     display(err)
    # end
    # # reffn = joinpath(refdir, "ref$idx.png")

    # the test
    vtest = VisualTest(func, reffn, idx)
    test_images(vtest, popup=popup, sigma=sigma, tol=tol, newfn = newfn)
end

function image_comparison_facts(pkg::Symbol;
                                skip = [],      # skip these examples (int index)
                                only = nothing, # limit to these examples (int index)
                                debug = false,  # print debug information?
                                sigma = [1,1],  # number of pixels to "blur"
                                tol = 1e-2)     # acceptable error (percent)
  for i in 1:length(Plots._examples)
    i in skip && continue
    if only == nothing || i in only
      @test image_comparison_tests(pkg, i, debug=debug, sigma=sigma, tol=tol) |> success == true
    end
  end
end


text:
#####
##### apply flux boundary conditions to tracer tendency `Gc`
#####

# a positive west flux is associated with an *increase* in `Gc` near the west boundary.
# same for south and bottom.
@kernel function apply_west_bcs_kernel!(Gc, grid, west_bc, iter, ΔT)
    j, k = @index(Global, NTuple)
    jj = j + grid.Hy
    kk = k + grid.Hz
    @inbounds Gc[1+grid.Hx, jj, kk] += getbc(west_bc, j, k, iter) * ΔT * Ax(1+grid.Hx, jj, kk, grid) / volume(1+grid.Hx, jj, kk, grid)
end

@kernel function apply_south_bcs_kernel!(Gc, grid, south_bc, iter, ΔT)
    i, k = @index(Global, NTuple)
    ii = i + grid.Hx
    kk = k + grid.Hz
    @inbounds Gc[ii, 1+grid.Hy, kk] += getbc(south_bc, i, k, iter) * ΔT * Ay(ii, 1+grid.Hy, kk, grid) / volume(ii, 1+grid.Hy, kk, grid)
end

@kernel function apply_bottom_bcs_kernel!(Gc, grid, bottom_bc, iter, ΔT)
    i, j = @index(Global, NTuple)
    ii = i + grid.Hx
    jj = j + grid.Hy
    @inbounds Gc[ii, jj, grid.Nz+grid.Hz] += getbc(bottom_bc, i, j, iter) * ΔT * Az(ii, jj, grid.Nz+grid.Hz, grid) / volume(ii, jj, grid.Nz+grid.Hz, grid)
end

# a positive east flux is associated with an *decrease* in `Gc` near the east boundary.
# same for north and top.
@kernel function apply_east_bcs_kernel!(Gc, grid, east_bc, iter, ΔT)
    j, k = @index(Global, NTuple)
    jj = j + grid.Hy
    kk = k + grid.Hz
    @inbounds Gc[grid.Nx+grid.Hx, jj, kk] -= getbc(east_bc, j, k, iter) * ΔT * Ax(grid.Nx+grid.Hx, jj, kk, grid) / volume(grid.Nx+grid.Hx, jj, kk, grid)
end

@kernel function apply_north_bcs_kernel!(Gc, grid, north_bc, iter, ΔT)
    i, k = @index(Global, NTuple)
    ii = i + grid.Hx
    kk = k + grid.Hz
    @inbounds Gc[ii, grid.Ny+grid.Hy, kk] -= getbc(north_bc, i, k, iter) * ΔT * Ay(ii, grid.Ny+grid.Hy, kk, grid) / volume(ii, grid.Ny+grid.Hy, kk, grid)
end

@kernel function apply_top_bcs_kernel!(Gc, grid, top_bc, iter, ΔT)
    i, j = @index(Global, NTuple)
    ii = i + grid.Hx
    jj = j + grid.Hy
    @inbounds Gc[ii, jj, 1+grid.Hz] -= getbc(top_bc, i, j, iter) * ΔT * Az(ii, jj, 1+grid.Hz, grid) / volume(ii, jj, 1+grid.Hz, grid)
end

##### apply flux boundary conditions in x direction (west and east)
function apply_west_bcs!(Gc, grid, west_bc, iter, ΔT, arch, dep)
    kernel! = apply_west_bcs_kernel!(device(arch), (16,16), (grid.Ny, grid.Nz))
    event = kernel!(Gc, grid, west_bc, iter, ΔT, dependencies=dep)
    wait(device(arch), event)
    return nothing
end
function apply_east_bcs!(Gc, grid, east_bc, iter, ΔT, arch, dep)
    kernel! = apply_east_bcs_kernel!(device(arch), (16,16), (grid.Ny, grid.Nz))
    event = kernel!(Gc, grid, east_bc, iter, ΔT, dependencies=dep)
    wait(device(arch), event)
    return nothing
end

##### apply flux boundary conditions in y direction (south and north)
function apply_south_bcs!(Gc, grid, south_bc, iter, ΔT, arch, dep)
    kernel! = apply_south_bcs_kernel!(device(arch), (16,16), (grid.Nx, grid.Nz))
    event = kernel!(Gc, grid, south_bc, iter, ΔT, dependencies=dep)
    wait(device(arch), event)
    return nothing
end
function apply_north_bcs!(Gc, grid, north_bc, iter, ΔT, arch, dep)
    kernel! = apply_north_bcs_kernel!(device(arch), (16,16), (grid.Nx, grid.Nz))
    event = kernel!(Gc, grid, north_bc, iter, ΔT, dependencies=dep)
    wait(device(arch), event)
    return nothing
end

##### apply flux boundary conditions in z direction (bottom and top)
function apply_bottom_bcs!(Gc, grid, bottom_bc, iter, ΔT, arch, dep)
    kernel! = apply_bottom_bcs_kernel!(device(arch), (16,16), (grid.Nx, grid.Ny))
    event = kernel!(Gc, grid, bottom_bc, iter, ΔT, dependencies=dep)
    wait(device(arch), event)
    return nothing
end
function apply_top_bcs!(Gc, grid, top_bc, iter, ΔT, arch, dep)
    kernel! = apply_top_bcs_kernel!(device(arch), (16,16), (grid.Nx, grid.Ny))
    event = kernel!(Gc, grid, top_bc, iter, ΔT, dependencies=dep)
    wait(device(arch), event)
    return nothing
end

##### apply flux boundary conditions in all three directions, six faces.
function apply_bcs!(Gcs, nuts, grid, iter, ΔT, arch)
    barrier = Event(device(arch))
    events = []
    for name in nut_names
        west_event   =   apply_west_bcs!(Gcs[name].data, grid, nuts[name].bc.west,   iter, ΔT, arch, barrier)
        east_event   =   apply_east_bcs!(Gcs[name].data, grid, nuts[name].bc.east,   iter, ΔT, arch, barrier)
        south_event  =  apply_south_bcs!(Gcs[name].data, grid, nuts[name].bc.south,  iter, ΔT, arch, barrier)
        north_event  =  apply_north_bcs!(Gcs[name].data, grid, nuts[name].bc.north,  iter, ΔT, arch, barrier)
        bottom_event = apply_bottom_bcs!(Gcs[name].data, grid, nuts[name].bc.bottom, iter, ΔT, arch, barrier)
        top_event    =    apply_top_bcs!(Gcs[name].data, grid, nuts[name].bc.top,    iter, ΔT, arch, barrier)
        push!(events, west_event, east_event, south_event, north_event, bottom_event, top_event)
    end

    # filter out ::Nothing
    events = filter(e -> typeof(e) <: Event, events)

    wait(device(arch), MultiEvent(Tuple(events)))
    return nothing
end

# Avoid some computation / memory accesses when no flux boundary conditions will be applied
  apply_west_bcs!(Gc, grid, ::Nothing, iter, ΔT, arch, dep) = nothing
  apply_east_bcs!(Gc, grid, ::Nothing, iter, ΔT, arch, dep) = nothing
 apply_north_bcs!(Gc, grid, ::Nothing, iter, ΔT, arch, dep) = nothing
 apply_south_bcs!(Gc, grid, ::Nothing, iter, ΔT, arch, dep) = nothing
apply_bottom_bcs!(Gc, grid, ::Nothing, iter, ΔT, arch, dep) = nothing
   apply_top_bcs!(Gc, grid, ::Nothing, iter, ΔT, arch, dep) = nothing


text:
################################################################################
#
#  The main class group part
#
################################################################################

function show(io::IO, c::ClassGrpCtx)
  println(io, "Ctx for class group of ", order(c.FB.ideals[1]))
  println(io, "Factorbase with ", length(c.FB.ideals), " ideals of norm up to ", norm(c.FB.ideals[1]))
  println(io, "Relations module: ", c.M)
end

function order(c::ClassGrpCtx)
  return ring(c.FB)
end

function nf(c::ClassGrpCtx)
  return nf(order(c))
end

function class_group_init(FB::NfFactorBase, T::DataType = SMat{fmpz}; add_rels::Bool = true, use_aut::Bool = false)
  O = order(FB.ideals[1])
  n = degree(O)
  clg = ClassGrpCtx{T}()

  clg.FB = FB

  clg.bad_rel = 0
  clg.rel_cnt = 0
  clg.last = 0

  clg.M = ModuleCtx_fmpz(length(FB.ideals))
  clg.R_gen = Vector{nf_elem}()
  clg.R_rel = Vector{nf_elem}()

  clg.c = conjugates_init(nf(O).pol)
  add_rels && for I in clg.FB.ideals
    a = I.gen_one
    class_group_add_relation(clg, nf(O)(a), fmpq(abs(a)^n), fmpz(1), orbit = false)
    b = nf(O)(I.gen_two)
    bn = norm_div(b, fmpz(1), 600)
    if nbits(numerator(bn)) < 550
      class_group_add_relation(clg, b, abs(bn), fmpz(1), orbit = false)
    end
  end

  l = zero_matrix(FlintZZ, n, 1+clg.c.r2)
  for i = 1:n
    l[i,1] = 1
  end
  for i = 1:clg.c.r2
    l[clg.c.r1+2*i, i+1] = 1
    l[clg.c.r1+2*i-1, i+1] = -1
  end
  # what I want is a lll-reduced basis for the kernel
  # it probably should be a sep. function
  # however, there is nullspace - which is strange...
  l, t = hnf_with_transform(l)
  if 1 + clg.c.r2 + 1 > nrows(l)
    t = zero_matrix(FlintZZ, 0, 0)
  else
    t = view(t, (1+clg.c.r2+1):nrows(l), 1:nrows(l))
  end
  l = lll(t)
  clg.val_base = l

  if use_aut
    au = automorphisms(nf(O), copy = false)
    clg.aut_grp = class_group_add_auto(clg, au)
    clg.normCtx = NormCtx(O, div(nbits(discriminant(O)), 2) + 20,
                                                     length(au) == degree(O))
  else
    clg.normCtx = NormCtx(O, div(nbits(discriminant(O)), 2) + 20, false)
  end

  return clg
end

function class_group_init(O::NfOrd, B::Int; min_size::Int = 20, add_rels::Bool = true,
                          use_aut::Bool = false,
                          complete::Bool = true, degree_limit::Int = 0, T::DataType = SMat{fmpz})
  @vprint :ClassGroup 2 "Computing factor base ...\n"

  @assert B > 0

  FB = NfFactorBase()
  while true
    FB = NfFactorBase(O, B, complete = complete, degree_limit = degree_limit)
    if length(FB.ideals) > min_size
      break
    end
    B *= 2
    @vprint :ClassGroup 2 "Increasing bound to $B ...\n"
  end
  @vprint :ClassGroup 2 " done\n"
  return class_group_init(FB, T, add_rels = add_rels, use_aut = use_aut)
end

function _get_autos_from_ctx(ctx::ClassGrpCtx)
  return ctx.aut_grp::Vector{Tuple{NfToNfMor, Perm{Int}}}
end

################################################################################
#
#  Conversion to Magma
#
################################################################################
function to_magma(f::IOStream, clg::ClassGrpCtx)
  print(f, "K<a> := NumberField(", nf(order(clg.FB.ideals[1])).pol, ");\n");
  print(f, "M := MaximalOrder(K);\n");
  print(f, "fb := [ ")
  for i=1:clg.FB.size
    to_magma(f, clg.FB.ideals[i], "M")
    if i < clg.FB.size
      print(f, ",\n")
    end
  end
  print(f, "];\n")

  R = vcat(clg.R_gen, clg.R_rel)
  print(f, "R := [ ")
  for i = 1:length(R)
    print(f, R[i])
    if i < length(R)
      print(f, ",\n")
    end
  end
  print(f, "];\n")

  to_magma(f, clg.M, name = "MM")
end

function to_magma(s::String, c::ClassGrpCtx)
  f = open(s, "w")
  to_magma(f, c)
  close(f)
end

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

function to_hecke(f::IOStream, clg::ClassGrpCtx; field_name = "K")
  O = order(clg.FB.ideals[1])
  L = nf(O)
  varr = string(L.S)
  print(f, "$field_name, $vvar = number_field(", nf(order(clg.FB.ideals[1])).pol, ", \"$varr\");\n");
  O = order(clg.FB.ideals[1])
  to_hecke(f, basis(O))
  print(f, "O = Order($field_name, map($field_name, R))\n")
  print(f, "O.ismaximal = 1\n")

  print(f, "c = Hecke.class_group_init(O, ", norm(clg.FB.ideals[1]), ")\n")

  #_print_nf_elem_array_serialize(a)

  R = vcat(clg.R_gen, clg.R_rel)
  for i = 1:length(R)
    print(f, "Hecke.class_group_add_relation(c, $field_name(", R[i], "))\n")
  end

end

function to_hecke(s::String, c::ClassGrpCtx; T...)
  f = open(s, "w")
  to_hecke(f, c; T...)
  close(f)
end


text:
module _Inplace

using Brainstorm._Algorithm: binarysearch
using Brainstorm._DataStructure: PreAllocatedStack
using Brainstorm._Math._GMP: set!, add!, mul!, div!, pow!

include("binarysplitting.jl")

preallocatedstack(size) = PreAllocatedStack(0, [(big(0), big(0)) for i in 1:size])

# Compute int(e * 10^digits)
# This is done using Taylor series with binary splitting
euler(digits::Integer) = begin
    # How many terms to compute
    #
    # http://www.johndcook.com/blog/2011/06/10/stirling-approximation/
    #
    #             k! > 10^digits
    # =>      log k! > digits
    # => k log k - k > digits
    #
    # with log n! = sum_k=1^n log k
    #            >= int_x=1^n log x dx
    #             = n log n - n + 1
    #             > n log n - n
    #
    n = binarysearch(termfunc, digits + 1, one(digits), digits + 100)[2]

    depth = 1
    while (1 << depth) <= n; depth += 1; end
    s = preallocatedstack(depth + 1)

    # calculate P(0, N) and Q(0, N)
    binarysplitting!(s, zero(n), n)
    p, q = s[1]
    o, _ = s[2]

    set!(o, 10)
    pow!(o, digits) # one

    mul!(p, o)
    div!(p, q)
    add!(p, o)
end

termfunc(k) = k * (ndigits(k) - 2)

compute!(s, a, b) = begin
    p, q = s[1]

    set!(p, 1)
    set!(q, b)
end

combine!(s) = begin
    p₁, q₁ = s[1]
    p₂, q₂ = s[2]

    # p = p₁ * q₂ + p₂
    mul!(p₁, q₂)
    add!(p₁, p₂)

    # q = q₁ * q₂
    mul!(q₁, q₂)
end

end


text:
@doc raw"""
    rsi(ta::TimeArray, n = 14; wilder = false)

Relative Strength Index

# Formula

```math
\text{RSI} = \frac{\text{EMA}(\text{Up}, n)}{\text{EMA}(\text{Up}, n) + \text{EMA}(\text{Dn}, n)}
```
"""
function rsi(ta::TimeArray, n::Int = 14; wilder::Bool = false)
  ret = diff(values(ta), dims=1)
  ups = zeros(size(values(ta), 1) - 1, size(values(ta), 2))
  dns = zeros(size(values(ta), 1) - 1, size(values(ta), 2))

  #skip NaNs (by explicitly comparing on both sides of zero)
  @inbounds for i in 1:size(values(ta), 1) - 1
    for j in 1:size(values(ta), 2)
      if ret[i, j] >= 0
        ups[i, j] += ret[i, j]
      elseif ret[i, j] < 0
        dns[i, j] += ret[i, j]
      end
    end
  end

  upsema = ema(ups, n, wilder=wilder)
  dnsema = abs.(ema(dns, n, wilder=wilder))
  rs     = upsema ./ dnsema

  res = @. 100 - (100 / (1 + rs))

  cname = [Symbol(c, "_rsi_", n) for c ∈ colnames(ta)]

  TimeArray(timestamp(ta)[n+1:end], res, cname, meta(ta))
end

@doc raw"""
    cci(ohlc::TimeArray, ma = 20, c = 0.015)

Commodity Channel Index

# Formula

```math
CCI = \frac{P^{typical} - \text{SMA}(P^{typical})}{c \times \sigma(P^{typical})}
```

# References

- [Wikipedia](https://en.wikipedia.org/wiki/Commodity_channel_index)

"""
function cci(ohlc::TimeArray, ma::Int = 20, c::AbstractFloat = 0.015)
  pt = typical(ohlc)
  cci = safediv.((pt .- sma(pt, ma)), moving(mean_abs_dev, pt, ma)) ./ c
  rename(cci, :cci)
end

@doc raw"""
    chaikinoscillator(ohlcv::TimeArray, fast = 3, slow = 10; h = :High, l = :Low, c = :Close)

Chaikin Oscillator

Developed by Marc Chaikin.

# Formula

```math
\text{Chaikin OSC} = \text{EMA}(\text{ADL}, \text{fast}) - \text{EMA}(\text{ADL}, \text{slow})
```

where the [`adl`](@ref) is the Accumulation/Distribution Line.

# References

- [StockCharts]
  (http://stockcharts.com/school/doku.php?id=chart_school:technical_indicators:chaikin_oscillator)
"""
function chaikinoscillator(ohlcv::TimeArray, fast::Integer = 3, slow::Integer = 10;
                           h = :High, l = :Low, c = :Close, v = :Volume)
  _adl = adl(ohlcv, h=h, l=l, c=c, v=v)
  rename(ema(_adl, fast) .- ema(_adl, slow), [:chaikinoscillator])
end

@doc raw"""
    macd(ta, fast=12, slow=26, signal=9; wilder=false)

Moving Average Convergence / Divergence

# Formula

```math
\begin{align*}
  \text{MACD Bar} & = \text{DIF} - \text{DEM} \\
  \text{DIF} & = \text{EMA}(P^\text{Close}, \text{fast}) - \text{EMA}(P^\text{Close}, \text{slow}) \\
  \text{DEM} & = \text{EMA}(\text{DIF}, 9) \tag{signal}
\end{align*}
```

# Return

`TimeArray` with 3 columns `[:macd, :dif, :signal]`.

If the input is a multi-column `TimeArray`, the new column names will be
`[:A_macd, :B_macd, :A_dif, :B_dif, :A_signal, :B_signal]`.

"""
function macd(ta::TimeArray{T,N}, fast::Integer = 12, slow::Integer = 26, signal::Integer = 9;
              wilder::Bool = false) where {T, N}
  dif = ema(ta, fast, wilder=wilder) .- ema(ta, slow, wilder=wilder)
  sig = ema(dif, signal, wilder=wilder)
  osc = dif .- sig

  cols = [:macd, :dif, :signal]
  new_cols = (N > 1) ? gen_colnames(colnames(ta), cols) : cols

  merge(merge(osc, dif), sig, colnames=new_cols)
end

@doc raw"""
    roc(ta::TimeArray, n)

Rate of Change

# Formula

```math
\text{ROC} = \frac{P_{t} - P_{t-n}}{P_{t-n}}
```

# References

- [Wikipedia](https://en.wikipedia.org/wiki/Momentum_(technical_analysis))

"""
function roc(ta::TimeArray, n::Integer)
  prev = lag(ta, n)
  rename((ta .- prev) ./ prev, [Symbol("$(c)_roc_$(n)") for c ∈ colnames(ta)])
end

@doc raw"""
    aroon(ohlc::TimeArray, n = 25; h = :High, l = :Low)

Aroon Oscillator

# Formula

```math
\begin{align*}
  \text{up}   & = \frac{\mathop{argmax}(\text{High}_{t-n} \dots \text{High}_t)}{n} \times 100 \\
  \text{down} & = \frac{\mathop{argmin}(\text{Low}_{t-n} \dots \text{Low}_t)}{n} \times 100 \\
  \text{osc}  & = \text{up} - \text{down}
\end{align*}
```

# References

- [StockCharts]
  (http://stockcharts.com/school/doku.php?id=chart_school:technical_indicators:aroon_oscillator)
"""
function aroon(ohlc::TimeArray, n::Integer = 25; h = :High, l = :Low)
  up = rename(moving(nanargmax, ohlc[h], n) ./ n .* 100, :up)
  dn = rename(moving(nanargmin, ohlc[l], n) ./ n .* 100, :down)
  osc = rename(up .- dn, :osc)

  merge(merge(up, dn), osc)
end

@doc raw"""
    adx(ohlc::TimeArray, n = 14; h = :High, l = :Low, c = :Close)

Average Directional Movement Index

Developed by J. Welles Wilder.
This implementation follows StockCharts.

# References

- [wikipedia]
  (https://en.wikipedia.org/wiki/Average_directional_movement_index)

- [StockCharts]
  (http://stockcharts.com/school/doku.php?id=chart_school:technical_indicators:average_directional_index_adx)
"""
function adx(ohlc::TimeArray, n::Integer = 14; h = :High, l = :Low, c = :Close)
  dm = relu(merge(ohlc[h] .- lag(ohlc[h]), lag(ohlc[l]) .- ohlc[l]))

  _dm_values = values(dm)
  _dm_values[findmin(_dm_values, dims=2)[2]] .= 0

  dm = TimeArray(timestamp(dm), _dm_values, colnames(dm))

  dm = wilder_smooth(dm, n)

  tr = wilder_smooth(truerange(ohlc, h=h, l=l, c=c), n)

  di = (dm ./ tr) .* 100
  di = rename(di, [:di_plus, :di_minus])

  dx = @. abs((di[:di_plus] - di[:di_minus]) / (di[:di_plus] + di[:di_minus])) * 100
  adx = wilder_smooth(dx, n, dx=true)

  rename(merge(adx, merge(dx, di)), ([:adx, :dx, :di_plus, :di_minus]))
end

@doc raw"""
    stochasticoscillator(ohlc, n = 14, fast_d = 3, slow_d = 3; h = :High, l = :Low, c = :Close)

Stochastic Oscillator

A.k.a *%K%D*, or *KD*

# Arguments

- `n`: period of fast(raw) `%K`

- `fast_d`: MA period of fast `%D`

- `slow_d`: MA period of slow `%D`

# Formula

```math
\begin{align*}
  \text{fast %K} & = \frac{\text{Close}_t - \max(\text{High}_{t-n}, \dots, \text{High}_t)}
        {\max(\text{High}_{t-n}, \dots, \text{High}_t) - \min(\text{Low}_{t-n}, \dots, \text{Low}_t)}
        \times 100 \\
  \text{fast %D} & = \text{SMA}(\text{fast %K}) \\
  \text{slow %D} & = \text{SMA}(\text{fast %D})
\end{align*}
```

# References

- [Wikipedia]
  (https://en.wikipedia.org/wiki/Stochastic_oscillator)

- [FMLabs]
  (http://www.fmlabs.com/reference/default.htm?url=StochasticOscillator.htm)
"""
function stochasticoscillator(ohlc::TimeArray, n::Integer = 14, fast_d::Integer = 3,
                              slow_d::Integer = 3; h = :High, l = :Low, c = :Close)
  high = moving(nanmax, ohlc[h], n)
  low = moving(nanmin, ohlc[l], n)
  fast_k = rename((ohlc[c] .- low) ./ (high .- low) .* 100, :fast_k)
  fast_d = rename(sma(fast_k, fast_d), :fast_d)
  slow_d = rename(sma(fast_d, slow_d), :slow_d)
  stch = merge(merge(fast_k, fast_d), slow_d)
end

"""
    vortex(ohlc::TimeArray, n = 14; h = :High, l = :Low, c = :Close)

Vortex Indicator

It consists of two oscillators that capture positive and negative trend
movement. A bullish signal triggers when the positive trend indicator
crosses above the negative trend indicator or a key level.

# References

- [StockCharts]
  (http://stockcharts.com/school/doku.php?id=chart_school:technical_indicators:vortex_indicator)
"""
function vortex(ohlc::TimeArray, n::Integer = 14; h = :High, l = :Low, c = :Close)
  _h = ohlc[h]
  _l = ohlc[l]
  _c = ohlc[c]

  _lagc = lagfill(_c, 1, nanmean(values(_c)))

  tr = max.(_h, _lagc) .- min.(_l, _lagc)
  trn = moving(nansum, tr, n)

  _lagl = lagfill(_l, 1, nanmean(values(_l)))
  _lagh = lagfill(_h, 1, nanmean(values(_h)))

  vp = moving(nansum, abs.(_h .- _lagl), n) ./ trn
  vm = moving(nansum, abs.(_l .- _lagh), n) ./ trn

  vi = rename(merge(vp, vm), [:v_plus, :v_minus])
end

"""
    trix(ta::TimeArray, n = 14)

TRIX

Shows the percent rate of change of a triple exponentially smoothed moving
average.

# References

- [StockCharts]
  (http://stockcharts.com/school/doku.php?id=chart_school:technical_indicators:trix)
"""
function trix(ta::TimeArray, n::Integer = 14)
  ema1 = ema(ta, n)
  ema2 = ema(ema1, n)
  ema3 = ema(ema2, n)

  _lagema3 = lagfill(ema3, 1, nanmean(values(ema3))[1])

  trix = (ema3 .- _lagema3) ./ _lagema3
  trix = rename(100 .* trix, [:trix])
end

"""
    massindex(ohlc::TimeArray, n = 14, n2 = 25; h = :High, l = :Low)

Mass Index

It uses the high-low range to identify trend reversals based on range
expansions. It identifies range bulges that can foreshadow a reversal of
the current trend.

# References
- [StockCharts]
  (http://stockcharts.com/school/doku.php?id=chart_school:technical_indicators:mass_index)
"""
function massindex(ohlc::TimeArray, n::Integer = 14, n2::Integer = 25; h = :High, l = :Low)
  _h = ohlc[h]
  _l = ohlc[l]

  amplitude = _h .- _l

  ema1 = ema(amplitude, n)
  ema2 = ema(ema1, n)
  mass = ema1 ./ ema2
  mass = rename(moving(nansum, mass, n2), [:mi])
end

"""
    dpo(ta::TimeArray, n = 14)

Detrended Price Oscillator (DPO)

Is an indicator designed to remove trend from price and make it easier to
identify cycles.

# References
- [StockCharts]
  (http://stockcharts.com/school/doku.php?id=chart_school:technical_indicators:detrended_price_osci)
"""
function dpo(ta::TimeArray, n::Integer = 14)
  lagPeriod = Int(floor(0.5*n) + 1)
  _lagc = lagfill(ta, lagPeriod, nanmean(values(ta)))

  dpo = rename(_lagc .- moving(nanmean, ta, n), [:dpo])
end

"""
    kst(ta::TimeArray, r1 = 10, r2 = 15, r3 = 20, r4 = 30,
        n1 = 10, n2 = 10, n3 = 10, n4 = 15, nsig = 9)

KST Oscillator (KST)

It is useful to identify major stock market cycle junctures because its
formula is weighed to be more greatly influenced by the longer and more
dominant time spans, in order to better reflect the primary swings of stock
market cycle.

# References
- [Wikipedia](https://en.wikipedia.org/wiki/KST_oscillator)
"""
function kst(ta::TimeArray, r1::Integer = 10, r2::Integer = 15, r3::Integer = 20, r4::Integer = 30,
             n1::Integer = 10, n2::Integer = 10, n3::Integer = 10, n4::Integer = 15, nsig::Integer = 9)
  _meanc = nanmean(values(ta))
  _lagcr1 = lagfill(ta, r1, _meanc)
  _lagcr2 = lagfill(ta, r2, _meanc)
  _lagcr3 = lagfill(ta, r3, _meanc)
  _lagcr4 = lagfill(ta, r4, _meanc)

  rocma1 = (ta .- _lagcr1) ./ moving(nanmean, _lagcr1, n1, padding = true)
  rocma2 = (ta .- _lagcr2) ./ moving(nanmean, _lagcr2, n2, padding = true)
  rocma3 = (ta .- _lagcr3) ./ moving(nanmean, _lagcr3, n3, padding = true)
  rocma4 = (ta .- _lagcr4) ./ moving(nanmean, _lagcr4, n4, padding = true)

  kst = 100 .* (rocma1 .+ (2 .* rocma2) .+ (3 .* rocma3) .+ (4 .* rocma4))
  kst_sig = moving(nanmean, kst, nsig)

  rename(merge(kst, kst_sig), [:kst, :signal])
end

"""
Ichimoku Kinkō Hyō (Ichimoku)

    It identifies the trend and look for potential signals within that trend.

    http://stockcharts.com/school/doku.php?id=chart_school:technical_indicators:ichimoku_cloud
"""
function ichimoku(ohlc::TimeArray, n1::Integer = 9, n2::Integer = 26, n3::Integer = 52;
                  visual::Bool = false, h = :High, l = :Low)
  _h = ohlc[h]
  _l = ohlc[l]

  conv = 0.5 .* (moving(nanmax, _h, n1) + moving(nanmin, _l, n1))
  base = 0.5 .* (moving(nanmax, _h, n2) + moving(nanmin, _l, n2))

  spana = 0.5 .* (conv + base)
  spanb = 0.5 .* (moving(nanmax, _h, n3) + moving(nanmin, _l, n3))

  if visual
      spana = lagfill(spana, n2, nanmean(values(spana)))
      spanb = lagfill(spana, n2, nanmean(values(spanb)))
  end

  ichimoku = rename(merge(spana, spanb), [:ichimoku_a, :ichimoku_b])
end

"""
Money Flow Index (MFI)

Uses both price and volume to measure buying and selling pressure. It is
positive when the typical price rises (buying pressure) and negative when
the typical price declines (selling pressure). A ratio of positive and
negative money flow is then plugged into an RSI formula to create an
oscillator that moves between zero and one hundred.

http://stockcharts.com/school/doku.php?id=chart_school:technical_indicators:money_flow_index_mfi
"""
function moneyflowindex(ohlcv::TimeArray, n::Integer = 14;
                        c = :Close, h = :High, l = :Low, v = :Volume)
  _v = ohlcv[v]
  typ = typical(ohlcv, c = c, h = h, l = l)

  rmf = typ .* _v
  values_rmf = values(rmf)

  lagtyp = lagfill(typ, 1, NaN)

  mf_pos_flag = values(typ .> lagtyp)
  mf_neg_flag = values(typ .< lagtyp)

  # Positive rate of money flow
  values_rmf_p = copy(values_rmf)
  values_rmf_p[mf_pos_flag] .= 0
  rmf_p = moving(nansum, TimeArray(timestamp(ohlcv), values_rmf_p, [:rmf_p]), n, padding = true)

  # Negative rate of money flow
  values_rmf_n = copy(values_rmf)
  values_rmf_n[mf_neg_flag] .= 0
  rmf_n = moving(nansum, TimeArray(timestamp(ohlcv), values_rmf_n, [:rmf_n]), n, padding = true)

  mfr = rmf_p ./ rmf_n

  mfi_values = 100 .- 100 ./ (1 .+ values(mfr))

  mfi = TimeArray(timestamp(ohlcv), mfi_values, [:mfi])
end

"""
TRUE STRENGTH INDEX

    Shows both trend direction and overbought/oversold conditions.

    https://en.wikipedia.org/wiki/True_strength_index
"""
function tsi(ohlc::TimeArray, r::Integer = 25, s::Integer = 13; c = :Close)
  _c = ohlc[c]
  _lagc = lagfill(_c, 1, nanmean(values(_c)))

  m = _c .- _lagc

  m1 = ema(ema(m, r), s)
  m2 = ema(ema(abs.(m), r), s)

  tsi = rename(100 .* (m1 ./ m2), [:tsi])
end

@doc raw"""
    ultimateoscillator(ohlc::TimeArray, s = 7, m = 14, len = 28,
                       ws = 4.0, wm = 2.0, wl = 1.0;
                       c = :Close, h = :High, l = :Low)

Ultimate Oscillator

Larry Williams' (1976) signal, a momentum oscillator designed to capture
momentum across three different timeframes.

```math
\begin{align*}
  \text{BP}        & = P^\text{Close} - \min(P^\text{Low}, P_{t-1}^\text{Close}) \\
  \text{TR}        & = \max(P^\text{High}, P_{t-1}^\text{Close}) -
                       \min(P^\text{Low},  P_{t-1}^\text{Close}) \\
  \text{Average7}  & = \frac{\text{7-period BP Sum}}{\text{7-period TR Sum}} \\
  \text{Average14} & = \frac{\text{14-period BP Sum}}{\text{14-period TR Sum}} \\
  \text{Average28} & = \frac{\text{28-period BP Sum}}{\text{28-period TR Sum}} \\
  \text{UO}        & = 100 \times
                       \frac{(4 \times \text{Average7}) +
                             (2 \times \text{Average14}) +
                             \text{Average28}}
                            {4 + 2 + 1}
\end{align*}
```

# References

- [StockCharts](http://stockcharts.com/school/doku.php?id=chart_school:technical_indicators:ultimate_oscillator)
"""
function ultimateoscillator(ohlc::TimeArray, s::Integer = 7, m::Integer = 14, len::Integer = 28,
                            ws::AbstractFloat = 4.0, wm::AbstractFloat = 2.0, wl::AbstractFloat = 1.0;
                            c = :Close, h = :High, l = :Low)
  _c = ohlc[c]
  _h = ohlc[h]
  _l = ohlc[l]

  _lagc = lagfill(_c, 1, nanmean(values(_c)))

  min_l_or_pc = min.(_l, _lagc)
  max_h_or_pc = max.(_h, _lagc)

  bp = _c .- min_l_or_pc
  tr = max_h_or_pc .- min_l_or_pc

  avg_s = moving(nansum, bp, s, padding = true) ./ moving(nansum, tr, s, padding = true)
  avg_m = moving(nansum, bp, m, padding = true) ./ moving(nansum, tr, m, padding = true)
  avg_l = moving(nansum, bp, len, padding = true) ./ moving(nansum, tr, len, padding = true)

  uo = rename(100.0 .* ((ws .* avg_s) .+ (wm .* avg_m) .+ (wl .* avg_l)) ./ (ws + wm + wl), [:uo])
end

@doc raw"""
    williamsr(ohlc::TimeArray, n = 14; c = :Close, h = :High, l = :Low)

Williams %R

Developed by Larry Williams, Williams %R is a momentum indicator that is
the inverse of the Fast Stochastic Oscillator. Also referred to as %R,
Williams %R reflects the level of the close relative to the highest high
for the look-back period. In contrast, the Stochastic Oscillator reflects
the level of the close relative to the lowest low. %R corrects for the
inversion by multiplying the raw value by -100. As a result, the Fast
Stochastic Oscillator and Williams %R produce the exact same lines, only
the scaling is different. Williams %R oscillates from 0 to -100 [1].

Readings from 0 to -20 are considered overbought. Readings from -80 to -100
are considered oversold.

Unsurprisingly, signals derived from the Stochastic Oscillator are also
applicable to Williams %R.

```math
\text{%R} = \frac{\text{Highest High} - P^\text{Close}}
                 {\text{Highest High} - \text{Lowest Low}}
                 \times -100
```

Lowest Low = lowest low for the look-back period.

Highest High = highest high for the look-back period.

%R is multiplied by -100 correct the inversion and move the decimal.

The Williams %R oscillates from 0 to -100. When the indicator produces
readings from 0 to -20, this indicates overbought market conditions. When
readings are -80 to -100, it indicates oversold market conditions [2].

# References

1. [StockCharts](http://stockcharts.com/school/doku.php?id=chart_school:technical_indicators:williams_r)
2. [Investopedia](https://www.investopedia.com/terms/w/williamsr.asp)

"""
function williamsr(ohlc::TimeArray, n::Integer = 14; c = :Close, h = :High, l = :Low)
  _c = ohlc[c]
  _h = ohlc[h]
  _l = ohlc[l]

  hh = moving(nanmax, _h, n) # highest high over lookback period lbp
  ll = moving(nanmin, _l, n) # lowest low over lookback period lbp

  wr = rename(-100 .* (hh .- _c) ./ (hh .- ll), [:wr])
end

@doc raw"""
    awesomeoscillator(ohlc::TimeArray, s = 5, n = 34; h = :High, l = :Low)

Awesome Oscillator

The Awesome Oscillator is an indicator used to measure market momentum. AO
calculates the difference of a 34 Period and 5 Period Simple Moving
Averages. The Simple Moving Averages that are used are not calculated
using closing price but rather each bar's midpoints. AO is generally used
to affirm trends or to anticipate possible reversals [1].

Awesome Oscillator is a 34-period simple moving average, plotted through
the central points of the bars (H+L)/2, and subtracted from the 5-period
simple moving average, graphed across the central points of the bars
``\frac{H + L}{2}`` [2].

```math
\begin{align*}
  \text{Median Price} & = \frac{P^\text{High} + P^\text{Low}}{2} \\
  \text{AO}           & = \text{SMA}(\text{Median Price}, 5) -
                          \text{SMA}(\text{Median Price}, 34)
\end{align*}
```

# References

1. [TradingView](https://www.tradingview.com/wiki/Awesome_Oscillator_(AO))
2. https://www.ifcm.co.uk/ntx-indicators/awesome-oscillator
"""
function awesomeoscillator(ohlc::TimeArray, s::Integer = 5, n::Integer = 34;
                           h = :High, l = :Low)
  _h = ohlc[h]
  _l = ohlc[l]

  mp = 0.5 .* (_h .+ _l)

  ao = rename(moving(nanmean, mp, s) .- moving(nanmean, mp, n), [:ao])
end


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

"""
    StochasticProgram

An instance of a stochastic optimization problem.
"""
struct StochasticProgram{N, S <: NTuple{N, Stage}, ST <: AbstractStochasticStructure{N}}
    stages::S
    decisions::Decisions{N}
    structure::ST
    generator::Dict{Symbol, Function}
    problemcache::Dict{Symbol, JuMP.Model}
    solutioncache::Dict{Symbol, SolutionCache}
    optimizer::StochasticProgramOptimizer

    function StochasticProgram(stages::NTuple{N, Stage},
                               scenario_types::ScenarioTypes{M},
                               instantiation::StochasticInstantiation,
                               optimizer_constructor) where {N, M}
        N >= 2 || error("Stochastic program needs at least two stages.")
        M == N - 1 || error("Inconsistent number of stages $N and number of scenario types $M")
        S = typeof(stages)
        decisions = Decisions(Val(N))
        optimizer = StochasticProgramOptimizer(optimizer_constructor)
        structure = StochasticStructure(decisions, scenario_types, default_structure(instantiation, optimizer.optimizer))
        ST = typeof(structure)
        return new{N, S, ST}(stages,
                             decisions,
                             structure,
                             Dict{Symbol, Function}(),
                             Dict{Symbol, JuMP.Model}(),
                             Dict{Symbol, SolutionCache}(),
                             optimizer)
    end

    function StochasticProgram(stages::NTuple{N, Stage},
                               scenarios::NTuple{M, Vector{<:AbstractScenario}},
                               instantiation::StochasticInstantiation,
                               optimizer_constructor) where {N, M}
        N >= 2 || error("Stochastic program needs at least two stages.")
        M == N - 1 || error("Inconsistent number of stages $N and number of scenario types $M")
        S = typeof(stages)
        decisions = Decisions(Val(N))
        optimizer = StochasticProgramOptimizer(optimizer_constructor)
        structure = StochasticStructure(decisions, scenarios, default_structure(instantiation, optimizer.optimizer))
        ST = typeof(structure)
        return new{N, S, ST}(stages,
                             decisions,
                             structure,
                             Dict{Symbol, Function}(),
                             Dict{Symbol, JuMP.Model}(),
                             Dict{Symbol, SolutionCache}(),
                             optimizer)
    end
end
TwoStageStochasticProgram{S <: Tuple{Stage, Stage}, ST <: AbstractStochasticStructure{2}} = StochasticProgram{2, S, ST}

# Constructors #
# ========================== #
# Two-stage
# ========================== #
"""
    StochasticProgram(first_stage_params::Any,
                      second_stage_params::Any,
                      instantiation::StochasticInstantiation,
                      optimizer_constructor=nothing) where T <: AbstractFloat

Create a new two-stage stochastic program with stage data given by `first_stage_params` and `second_stage_params`. After construction, scenarios of type `Scenario` can be added through `add_scenario!`. Optionally, a capable `optimizer_constructor` can be supplied to later optimize the stochastic program. If multiple Julia processes are available, the resulting stochastic program will automatically be memory-distributed on these processes. This can be avoided by setting `procs = [1]`.
"""
function StochasticProgram(first_stage_params::Any,
                           second_stage_params::Any,
                           instantiation::StochasticInstantiation,
                           optimizer_constructor = nothing)
    stages = (Stage(first_stage_params), Stage(second_stage_params))
    return StochasticProgram(stages, (Scenario,), instantiation, optimizer_constructor)
end
"""
    StochasticProgram(first_stage_params::Any,
                      second_stage_params::Any,
                      ::Type{Scenario},
                      instantiation::StochasticInstantiation,
                      optimizer_constructor=nothing) where Scenario <: AbstractScenario

Create a new two-stage stochastic program with stage data given by `first_stage_params` and `second_stage_params`. After construction, scenarios of type `S` can be added through `add_scenario!`. Optionally, a capable `optimizer_constructor` can be supplied to later optimize the stochastic program. If multiple Julia processes are available, the resulting stochastic program will automatically be memory-distributed on these processes. This can be avoided by setting `procs = [1]`.
"""
function StochasticProgram(first_stage_params::Any,
                           second_stage_params::Any,
                           ::Type{Scenario},
                           instantiation::StochasticInstantiation,
                           optimizer_constructor = nothing) where Scenario <: AbstractScenario
    stages = (Stage(first_stage_params), Stage(second_stage_params))
    return StochasticProgram(stages, (Scenario,), instantiation, optimizer_constructor)
end
"""
    StochasticProgram(::Type{Scenario},
                      instantiation::StochasticInstantiation,
                      optimizer_constructor=nothing) where Scenario <: AbstractScenario

Create a new two-stage stochastic program with scenarios of type `Scenario` and no stage data. Optionally, a capable `optimizer_constructor` can be supplied to later optimize the stochastic program.
"""
function StochasticProgram(::Type{Scenario},
                           instantiation::StochasticInstantiation,
                           optimizer_constructor = nothing) where Scenario <: AbstractScenario
    stages = (Stage(nothing), Stage(nothing))
    return StochasticProgram(stages, (Scenario,), instantiation, optimizer_constructor; procs = procs)
end
"""
    StochasticProgram(first_stage_params::Any,
                      second_stage_params::Any,
                      scenarios::Vector{<:AbstractScenario},
                      instantiation::StochasticInstantiation,
                      optimizer_constructor = nothing)

Create a new two-stage stochastic program with a given collection of `scenarios`. Optionally, a capable `optimizer_constructor` can be supplied to later optimize the stochastic program. If multiple Julia processes are available, the resulting stochastic program will automatically be memory-distributed on these processes. This can be avoided by setting `procs = [1]`.
"""
function StochasticProgram(first_stage_params::Any,
                           second_stage_params::Any,
                           scenarios::Vector{<:AbstractScenario},
                           instantiation::StochasticInstantiation,
                           optimizer_constructor = nothing)
    stages = (Stage(first_stage_params), Stage(second_stage_params))
    return StochasticProgram(stages, (scenarios,), instantiation, optimizer_constructor)
end
"""
    StochasticProgram(scenarios::Vector{<:AbstractScenario},
                      instantiation::StochasticInstantiation,
                      optimizer_constructor = nothing)

Create a new two-stage stochastic program with a given collection of `scenarios` and no stage data. Optionally, a capable `optimizer_constructor` can be supplied to later optimize the stochastic program. If multiple Julia processes are available, the resulting stochastic program will automatically be memory-distributed on these processes. This can be avoided by setting `procs = [1]`.
"""
function StochasticProgram(scenarios::Vector{<:AbstractScenario},
                           instantiation::StochasticInstantiation,
                           optimizer_constructor = nothing)
    stages = (Stage(nothing), Stage(nothing))
    return StochasticProgram(stages, (scenarios,), instantiation, optimizer_constructor)
end

function Base.copy(src::StochasticProgram{N}; instantiation = UnspecifiedInstantiation(), optimizer = nothing) where N
    stages = ntuple(Val(N)) do i
        Stage(stage_parameters(src, i))
    end
    scenario_types = ntuple(Val(N-1)) do i
        scenario_type(src, i+1)
    end
    dest = StochasticProgram(stages, scenario_types, instantiation, optimizer)
    merge!(dest.generator, src.generator)
    return dest
end

# Printing #
# ========================== #
function Base.show(io::IO, stochasticprogram::StochasticProgram{N}) where N
    plural(n) = (n == 1 ? "" : "s")
    stage(sp, s) = begin
        stage_key = Symbol(:stage_,s)
        ndecisions = num_decisions(sp, s)
        nscenarios = num_scenarios(sp, s)
        if s == 1
            if N == 2
                return " * $(ndecisions) decision variable$(plural(ndecisions))"
            else
                return " * Stage $s:\n   * $(ndecisions) decision variable$(plural(ndecisions))"
            end
        elseif s == 2 && N == 2
            stype = typename(scenario_type(stochasticprogram))
            return " * $(ndecisions) recourse variables\n * $(nscenarios) scenario$(plural(nscenarios)) of type $stype"
        else
            stype = typename(scenario_type(stochasticprogram, s))
            if distributed(stochasticprogram, s)
                if s == N
                    return " * Distributed stage $s:\n   * $(ndecisions) recourse variables\n * $(nscenarios) scenario$(plural(nscenarios)) of type $stype"
                else
                    return " * Distributed stage $s:\n   * $(ndecisions) decision variable$(plural(ndecisions))\n   * $(nscenarios) scenario$(plural(nscenarios)) of type $stype"
                end
            else
                if s == N
                    return " * Stage $s:\n   * $(ndecisions) recourse variables\n   * $(nscenarios) scenario$(plural(nscenarios)) of type $stype"
                else
                    return " * Stage $s:\n   * $(ndecisions) decision variable$(plural(ndecisions))\n   * $(nscenarios) scenario$(plural(nscenarios)) of type $stype"
                end
            end
        end
    end
    if deferred(stochasticprogram)
        n = num_scenarios(stochasticprogram)
        if n == 0
            return print(io, "Deferred stochastic program")
        else
            return print(io, "Deferred stochastic program with $n scenario$(plural(n))")
        end
    end
    if N == 2 && distributed(stochasticprogram)
        println(io, "Distributed stochastic program with:")
    else
        println(io, "Stochastic program with:")
    end
    print(io, stage(stochasticprogram, 1))
    for s = 2:N
        println(io,"")
        print(io, stage(stochasticprogram, s))
    end
    println(io,"")
    print(io, "Structure: ")
    print(io, structure_name(stochasticprogram))
    println(io,"")
    print(io, "Solver name: ")
    print(io, optimizer_name(stochasticprogram))
end
function Base.print(io::IO, stochasticprogram::StochasticProgram)
    if !deferred(stochasticprogram)
        # Delegate printing according to stochastic structure
        print(io, structure(stochasticprogram))
        print(io, "Solver name: ")
        print(io, optimizer_name(stochasticprogram))
    else
        # Just give summary if the stochastic program has not been initialized yet
        show(io, stochasticprogram)
    end
end
# ========================== #

# MOI #
# ========================== #
function MOI.get(stochasticprogram::StochasticProgram, attr::Union{MOI.TerminationStatus, MOI.PrimalStatus, MOI.DualStatus})
    # Check if there is a cached solution
    cache = solutioncache(stochasticprogram)
    if haskey(cache, :solution)
        # Returned cached solution if possible
        try
            return MOI.get(cache[:solution], attr)
        catch
        end
    end
    if haskey(cache, :node_solution_1)
        # Value was possibly only cached in first-stage solution
        try
            return MOI.get(cache[:node_solution_1], attr)
        catch
        end
    end
    return MOI.get(optimizer(stochasticprogram), attr)
end
function MOI.get(stochasticprogram::StochasticProgram, attr::MOI.AbstractModelAttribute)
    if MOI.is_set_by_optimize(attr)
        # Check if there is a cached solution
        cache = solutioncache(stochasticprogram)
        if haskey(cache, :solution)
            # Returned cached solution if possible
            try
                return MOI.get(cache[:solution], attr)
            catch
            end
        end
        if haskey(cache, :node_solution_1)
            # Value was possibly only cached in first-stage solution
            try
                return MOI.get(cache[:node_solution_1], attr)
            catch
            end
        end
        check_provided_optimizer(stochasticprogram.optimizer)
        if MOI.get(stochasticprogram, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
            throw(OptimizeNotCalled())
        end
        return MOI.get(optimizer(stochasticprogram), attr)
    else
        if is_structure_independent(attr)
            # Get attribute from first stage of proxy if structure independent
            return MOI.get(proxy(stochasticprogram, 1), attr)
        else
            # Handle in structure otherwise
            return MOI.get(structure(stochasticprogram), attr)
        end
    end
end
function MOI.get(stochasticprogram::StochasticProgram, attr::ScenarioDependentModelAttribute)
    if MOI.is_set_by_optimize(attr)
        # Check if there is a cached solution
        cache = solutioncache(stochasticprogram)
        key = Symbol(:node_solution_, attr.stage, :_, attr.scenario_index)
        if haskey(cache, key)
            try
                return MOI.get(cache[key], attr.attr)
            catch
            end
        end
        check_provided_optimizer(stochasticprogram.optimizer)
        # Return statuses without checks
        if typeof(attr.attr) <: Union{MOI.TerminationStatus, MOI.PrimalStatus, MOI.DualStatus}
            try
                # Try to get scenario-dependent value directly
                return MOI.get(optimizer(stochasticprogram), attr)
            catch
                # Fallback to resolving scenario-dependence in structure if
                # not supported natively by optimizer
                return MOI.get(structure(stochasticprogram), attr)
            end
        end
        if MOI.get(stochasticprogram, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
            throw(OptimizeNotCalled())
        end
        try
            # Try to get scenario-dependent value directly
            return MOI.get(optimizer(stochasticprogram), attr)
        catch
            # Fallback to resolving scenario-dependence in structure if
            # not supported natively by optimizer
            MOI.get(structure(stochasticprogram), attr)
        end
    else
        if is_structure_independent(attr)
            # Get attribute from first stage of proxy if structure independent
            return MOI.get(proxy(stochasticprogram, 1), attr)
        else
            # Handle in structure otherwise
            return MOI.get(structure(stochasticprogram), attr)
        end
    end
end
function MOI.get(stochasticprogram::StochasticProgram, attr::MOI.AbstractOptimizerAttribute)
    MOI.get(optimizer(stochasticprogram), attr)
end

function MOI.set(sp::StochasticProgram, attr::MOI.AbstractOptimizerAttribute, value)
    MOI.set(optimizer(sp), attr, value)
    return nothing
end
function MOI.set(sp::StochasticProgram, attr::MOI.AbstractModelAttribute, value)
    if is_structure_independent(attr)
        MOI.set(proxy(sp, attr.stage), attr, value)
    end
    MOI.set(structure(sp), attr, value)
    return nothing
end
function MOI.set(sp::StochasticProgram, attr::ScenarioDependentModelAttribute, value)
    if is_structure_independent(attr)
        MOI.set(proxy(sp, attr.stage), attr, value)
    end
    MOI.set(structure(sp), attr, value)
    return nothing
end

function MOI.set(sp::StochasticProgram, attr::MOI.Silent, flag)
    # Ensure that Silent is always passed
    MOI.set(structure(sp), attr, flag)
    # Pass to optimizer anyway
    MOI.set(optimizer(sp), attr, flag)
    return nothing
end

function JuMP.check_belongs_to_model(con_ref::ConstraintRef{<:StochasticProgram}, stochasticprogram::StochasticProgram)
    if owner_model(con_ref) !== model
        throw(ConstraintNotOwned(con_ref))
    end
end

Base.broadcastable(sp::StochasticProgram) = Ref(sp)


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

using Markdown
using InteractiveUtils

# ╔═╡ 755e081d-7176-4325-b481-d5d67db6176c
begin
    import Pkg
    Pkg.activate(mktempdir())
    Pkg.add([
        Pkg.PackageSpec(name="Plots", version="1"),
        Pkg.PackageSpec(name="PlutoUI", version="0.7"),
	    Pkg.PackageSpec(name="DifferentialEquations", version="6"),
	    Pkg.PackageSpec(name="Parameters", version="0.12"),
        Pkg.PackageSpec(name="LabelledArrays", version="1"),
        Pkg.PackageSpec(name="Setfield", version="0.7"),
    ])
    using Plots, PlutoUI, LinearAlgebra, DifferentialEquations, Parameters, LabelledArrays, Setfield

	Plots.gr(fmt=:png, lw=2)

	PlutoUI.TableOfContents()
end

# ╔═╡ 4a209b69-b1b4-4a63-8786-c2b2a1a9b6ee
md"""
# Figure 4.1 to 4.5

## Figure 4.1, 4.2, and 4.3

Steady states and phase plots in an assymetric network
"""

# ╔═╡ 9b0868de-687e-44f0-a3b8-2321cebd3a52
params = (K1=20.0, K2=5.0, K3=5.0, K4=5.0, K5=2.0, N=4)

# ╔═╡ 69f0e704-2a02-430d-97fd-4a7e5b1fc031
u0s = (LVector(A=0.0, B=0.0),
       LVector(A=0.5, B=0.6),
       LVector(A=0.17, B=1.1),
       LVector(A=0.25, B=1.9),
       LVector(A=1.85, B=1.70))

# ╔═╡ 298447c2-7e64-422c-9d25-9fb4dee35441
md"""
## Figure 4.4, 4.5

Vector fields
"""

# ╔═╡ 0b3412c0-bde3-11eb-0869-17548acad535
md"""
## Appendix
Running environment and some auxillary functions.

**You can use [this helper tool](https://fonsp.com/article-test-3/pkghelper.html) to generate these commands!**
"""

# ╔═╡ 3d3f9ef7-4195-413a-a4d4-563f327ccda4
begin
	# Convenience functions
    hill(x, k) = x / (x + k)
    hill(x, k, n) = hill(x^n, k^n)
    exprel(x) = ifelse(x≈zero(x), one(x), x / expm1(x))
end

# ╔═╡ f8ff273a-e1c7-4f00-9f9d-0ecb99f90b1e
function model!(du, u, p, t)
    @unpack K1, K2, K3, K4, K5, N = p
    @unpack A, B = u

    v1 = K1 * hill(1, B, N)
    v5 = K5 * A

    du.A = v1 - v5 - K3 * A
    du.B = K2 + v5 - K4 * B
end

# ╔═╡ 632c9de0-dc67-4a7f-882d-e04403d24d0e
sols = map(u0 -> solve(ODEProblem(model!, u0, 1.5, params)), u0s);

# ╔═╡ 8d6fb7a3-b5b0-4d49-9b07-8dcb1e2105cc
plot(sols[1], xlabel="Time", ylabel="Concentration", title="Fig. 4.2A (Time series)")

# ╔═╡ afed1e2a-40a6-4ddf-8988-6a66661e5063
plot(sols[1], vars=(1, 2), xlabel="[A]", ylabel="[B]", aspect_ratio=:equal,
     title="Fig. 4.2B (Phase plot)", ylims=(0.0, 2.0), xlims=(0.0, 2.0),
     legend=nothing)

# ╔═╡ ff7eed37-4569-45c1-af5c-aebda8fb2003
begin
	p3 = plot()

	for sol in sols
		plot!(p3, sol, linealpha=0.5, legend = nothing)
	end

	plot!(p3, xlabel="Time", ylabel="Concentration", title="Fig. 4.3A (Time series)")
end

# ╔═╡ d9295b32-da22-45e6-910b-6bd1126bb92b
begin
	p4 = plot()

	for sol in sols
		plot!(p4, sol, vars=(1, 2), linealpha=0.7, legend = nothing)
	end

	plot!(p4, aspect_ratio=:equal, title="Fig. 4.3B (Phase plot)", xlabel="[A]", ylabel="[B]", ylims=(0.0, 2.0), xlims=(0.0, 2.0), size=(600, 600))
end

# ╔═╡ a243b6e3-dd6b-47d0-b807-db88622bb3a4
# vector field
function df(x, y, p = params, scale=20)
	u = LVector(A=x, B=y)
	du = similar(u)
	model!(du, u, p, 0.0)
	return du ./ (norm(du)^0.5 * scale)
end

# ╔═╡ 0dd10ca8-71eb-4f68-89d7-a4717e06afc1
begin
	xx = [x for y in 0.0:0.1:2.0, x in 0.0:0.1:2.0]
	yy = [y for y in 0.0:0.1:2.0, x in 0.0:0.1:2.0]
	p1 = quiver(xx, yy, quiver=df, line=(:lightblue))

	for sol in sols
		plot!(p1, sol, vars=(1, 2), legend = nothing, line=(:green))
	end

	plot!(p1, aspect_ratio=:equal, title="Fig. 4.4A (Phase Plot with vector field)",
      xlabel="[A]", ylabel="[B]", xlim=(0.0, 2.0), ylim=(0.0, 2.0), size=(600, 600))
	p1
end

# ╔═╡ 799cb2ff-a7a2-4947-99c2-d23c06f0b920
# Nullclines
begin
	nullcline_a(b, p) = p.K1 / (p.K5 + p.K4)  * hill(1, b, p.N)
	nullcline_a(b) = nullcline_a(b, params)
	nullcline_b(b, p) = (p.K4*b - p.K2) / p.K5
	nullcline_b(b) = nullcline_b(b, params)
end

# ╔═╡ 160b9f34-e089-40b1-8976-b1a168da2044
begin
	# Figure 4.5A
	p45a = plot(aspect_ratio=:equal, title="Fig. 4.5A, Phase plot with nullclines")

	# Phase plots
	for sol in sols
		plot!(p45a, sol, vars=(1, 2), linealpha=0.7, lab=nothing)
	end

	# Parametric plotting for nullcline
	plot!(p45a, nullcline_a, identity, 0.0, 2.0, label="A nullcline", line=(:black, :dot))
	plot!(p45a, nullcline_b, identity, 0.0, 2.0, label="B nullcline", line=(:black, :dash))
	plot!(p45a, xlim=(0.0, 2.0), ylim=(0.0, 2.0), legend=:bottomright, size=(600, 600), xlabel="[A]", ylabel="[B]")

	p45a
end

# ╔═╡ 30ecb42b-ed42-455a-a933-2f00a4495907
begin
	p45b = quiver(xx, yy, quiver=df, line=(:lightblue), title="Fig. 4.5B, Vector field with nullclines", xlabel="[A]", ylabel="[B]")
	plot!(p45b, nullcline_a, identity, 0.0, 2.0, label="A nullcline", line=(:black, :dot))
	plot!(p45b, nullcline_b, identity, 0.0, 2.0, label="B nullcline", line=(:black, :dash))
	plot!(p45b, aspect_ratio=1.0, xlim=(0.0, 2.0), ylim=(0.0, 2.0), legend=:bottomright, size=(600, 600))

	p45b
end

# ╔═╡ Cell order:
# ╠═4a209b69-b1b4-4a63-8786-c2b2a1a9b6ee
# ╠═f8ff273a-e1c7-4f00-9f9d-0ecb99f90b1e
# ╠═9b0868de-687e-44f0-a3b8-2321cebd3a52
# ╠═69f0e704-2a02-430d-97fd-4a7e5b1fc031
# ╠═632c9de0-dc67-4a7f-882d-e04403d24d0e
# ╠═8d6fb7a3-b5b0-4d49-9b07-8dcb1e2105cc
# ╠═afed1e2a-40a6-4ddf-8988-6a66661e5063
# ╠═ff7eed37-4569-45c1-af5c-aebda8fb2003
# ╠═d9295b32-da22-45e6-910b-6bd1126bb92b
# ╠═298447c2-7e64-422c-9d25-9fb4dee35441
# ╠═799cb2ff-a7a2-4947-99c2-d23c06f0b920
# ╠═a243b6e3-dd6b-47d0-b807-db88622bb3a4
# ╠═0dd10ca8-71eb-4f68-89d7-a4717e06afc1
# ╠═160b9f34-e089-40b1-8976-b1a168da2044
# ╠═30ecb42b-ed42-455a-a933-2f00a4495907
# ╠═0b3412c0-bde3-11eb-0869-17548acad535
# ╠═3d3f9ef7-4195-413a-a4d4-563f327ccda4
# ╠═755e081d-7176-4325-b481-d5d67db6176c

max_stars_repo_name: stevengogogo/pluto-notebooks

max_issues_repo_name: stevengogogo/pluto-notebooks

max_forks_repo_name: stevengogogo/pluto-notebooks


text:
# Holger Teichgraeber, Elias Kuepper, 2019

 ######################
 # TimeSeriesClustering
 # Analyzing clustering techniques as input for energy systems optimization
 #
 #####################
module TimeSeriesClustering
  using Reexport
  using LinearAlgebra
  using CSV
  using Clustering
  using DataFrames
  using Distances
  using StatsBase
  @reexport using FileIO
  using JuMP

   #TODO how to make PyPlot, PyCall, and TimeWarp optional? -> only import when needed

   export InputData,
          FullInputData,
          ClustData,
          ClustDataMerged,
          AbstractClustResult,
          ClustResultAll,
          ClustResult,
          SimpleExtremeValueDescr,
          load_timeseries_data,
          combine_timeseries_weather_data,
          extreme_val_output,
          simple_extr_val_sel,
          representation_modification,
          get_sup_kw_args,
          run_clust,
          run_battery_opt,
          run_gas_opt,
          data_type,
          get_EUR_to_USD, #TODO Check which of the following should really be exported
          z_normalize,
          undo_z_normalize,
          sakoe_chiba_band,
          kmedoids_exact,
          sort_centers,
          calc_SSE,
          find_medoids,
          resize_medoids

  include(joinpath("utils","datastructs.jl"))
  include(joinpath("utils","utils.jl"))
  include(joinpath("utils","load_data.jl"))
  include(joinpath("optim_problems","run_opt.jl"))
  include(joinpath("clustering","run_clust.jl"))
  include(joinpath("clustering","exact_kmedoids.jl"))
  include(joinpath("clustering","extreme_vals.jl"))
  include(joinpath("clustering","attribute_weighting.jl"))
  include(joinpath("clustering","intraperiod_segmentation.jl"))
end # module TimeSeriesClustering


text:
#####################################################################
# special.jl
#
# includes non-standard code I needed for papers and whatnot
######################################################################

function norm2(x::LocTuple, y::LocTuple)
	d1 = x[1] - y[1]
	d2 = x[2] - y[2]
	return sqrt(d1*d1 + d2*d2)
end

function sim(j::LocTuple, actions::Vector{Pose}, observations::Vector{Float64}, b::Vector{Matrix{Float64}}, m::SearchDomain, x::Vehicle; show_mean=false, show_cov=false)
	# always start the vehicle in the center
	# NO, don't do that
	#x.x = m.length / 2.0
	#x.y = m.length / 2.0
	theta!(m, j[1], j[2])

	# warn user if jammer is not where it should be
	#if abs(dx) > x.x || abs(dy) > x.y
	#	println("WARNING: jammer outside search domain.")
	#end

	# loop through all observations...
	for (oi,o) in enumerate(observations)
		#update!(f, x, o)
		#plot(m, b[oi], x, show_mean=true, show_cov=true)
		plot(m, b[oi], x, show_mean=show_mean, show_cov=show_cov, obs=o)
		savefig("temp_$(oi).png", format="png")
		hold(false)
		act!(m, x, actions[oi])
	end
end

export sim2
function sim2(j::LocTuple, states::Vector{Pose}, observations::Vector{Float64}, b::Vector{Matrix{Float64}}, m::SearchDomain, x::Vehicle; show_mean=false, show_cov=false, show_path::Bool=true)
	# always start the vehicle in the center
	# NO, don't do that
	#x.x = m.length / 2.0
	#x.y = m.length / 2.0
	theta!(m, j[1], j[2])

	path_x = Float64[]
	path_y = Float64[]

	# loop through all observations...
	for (oi,o) in enumerate(observations)
		if show_path
			#plot path
			plot(path_x, path_y, "k")
		end
		hold(true)
		plot(m, b[oi], x, show_mean=show_mean, show_cov=show_cov, obs=oi)
		#savefig("temp_$(oi).pdf", format="pdf", dpi=300)
		savefig("temp_$(oi).png", format="png")
		hold(false)
		push!(path_x, x.x)
		push!(path_y, x.y)
		x.x = states[oi][1]
		x.y = states[oi][2]
		x.heading = states[oi][3]
		#act!(m, x, actions[oi])
	end
end

# This was for AIAA GNC 2018 paper.
# This probably shouldn't be here.
export eval2
function eval2(j::LocTuple, states::Vector{Pose}, observations::Vector{Float64}, b::Vector{Matrix{Float64}}, m::SearchDomain, x::Vehicle; show_mean=false, show_cov=false, show_path::Bool=true)
	# always start the vehicle in the center
	# NO, don't do that
	#x.x = m.length / 2.0
	#x.y = m.length / 2.0
	theta!(m, j[1], j[2])

    front_cone = Float64[]
    rear_cone = Float64[]
    sides = Float64[]

	# loop through all observations...
	for (oi,o) in enumerate(observations)
        # get true bearing between x
        xt = (x.x,x.y,x.heading)
        rel_bearing = fit_180(xt[3] - true_bearing(xt, m.theta))
        if rel_bearing < 0.0
            rel_bearing = -1.0 * rel_bearing
        end

        cla()
		plot(m, b[oi], x, show_mean=show_mean, show_cov=show_cov, obs=oi)
        #title("o = $o")
        #hold(false)
        if rel_bearing < 60.0
            # expect observation 1
            title("expect 1, got $o")
            push!(front_cone, o)
        elseif rel_bearing < 120.0
            title("expect either, got $o")
            push!(sides, o)
        else
            title("expect 0, got $o")
            push!(rear_cone,o)
        end
        pause(0.15)

        # update vehicle position
		x.x = states[oi][1]
		x.y = states[oi][2]
		x.heading = states[oi][3]
	end
    return front_cone, sides, rear_cone
end

function simulate(m::SearchDomain, uav::SimUnit;
                  video::Bool=true,
                  pause_time=0.3
                 )

    # reset the filter, vehicle, and policy
    # TODO: I think I assume the SimUnit comes in clean and ready to go
    #reset!(uav.f)
    #reset!(m, uav.x)
    #reset!(uav.p)

    # What was the cost to getting this first observation?
    temp_cost = get_cost(uav, m)

    # before doing anything else, we observe
    #  and update filter once
    o = observe(m, uav.x)
    update!(uav, o)

    # This was our first step; steps count number of observations
    step_count = 1

    # plot if need be
    if video
        figure("Simulation")
        plot(m, uav.f, uav.x)
        title("i = $(step_count)")
    end

    while !is_complete(uav.f, uav.tc, step_count)
        # act
        a = action(m, uav, o)
        act!(m, uav.x, a)

        # get cost and update step count
        temp_cost += get_cost(uav, m, a)
        step_count += 1

        # observe and update
        o = observe(m, uav.x)
        update!(uav, o)

        # plot if need be
        if video
            pause(pause_time)
            figure("Simulation")
            cla()
            plot(m, uav.f, uav.x)
            max_b = maximum(uav.f.b)
            title("i = $(step_count), max = $(round(max_b,3))")
        end
    end

    return temp_cost
end


text:
module DoubleArrayTrie

export DATrie

"""
* child = parent.base + key
* child.check == parent
"""
immutable Node
  base::Int
  check::Int
end

isused(n::Node) = n.check > 0

type DATrie{T}
  nodes::Vector{Node}
  values::Vector{T}
  count::Int
  base0::Int
end

Base.length(trie::DATrie) = trie.count

function Base.get(trie::DATrie, keys::Vector{Int}, default)
  nodeid = 1
  for k in keys
    nodeid = get(trie, nodeid, k, 0)
    #println(nodeid)
    nodeid == 0 && return default
  end
  nodeid = get(trie, nodeid, 0, 0)
  @assert nodeid != 0
  trie.values[trie.nodes[nodeid].base]
end

function Base.get(trie::DATrie, nodeid::Int, key::Int, default)
  nextid = trie.nodes[nodeid].base + key
  (nextid <= 0 || nextid > length(trie.nodes)) && return default
  trie.nodes[nextid].check == nodeid ? nextid : default
end

function Base.resize!(trie::DATrie, len::Int)
  n = length(trie.nodes)
  resize!(trie.nodes, len)
  for i = n+1:len
    trie.nodes[i] = Node(0, 0)
  end
  trie
end

"""
    DATrie{T}(keys, values)

Construct DATrie.
"""
function DATrie{T}(keys::Vector{Vector{Int}}, values::Vector{T})
  @assert length(keys) == length(values)
  # initialize
  trie = DATrie([Node(0,1)], values, length(keys), 1)
  resize!(trie, length(keys)*2)
  keys = map(keys) do key
    any(k -> k <= 0, key) && throw("Key must be positive.")
    push!(copy(key), 0) # append '0' to the end
  end

  items = [(1,1:length(trie),1)]
  while length(items) > 0
    nodeid, range, depth = pop!(items)

    subranges = Range[]
    k = first(range)
    for i in range
      keys[i][depth] == keys[k][depth] && continue
      push!(subranges, k:i-1)
      k = i
    end
    push!(subranges, k:last(range))

    child_keys = map(r -> keys[first(r)][depth], subranges)
    base = append!(trie, nodeid, first(range), child_keys)

    for i = length(subranges):-1:1
      r = subranges[i]
      key = keys[first(r)][depth]
      key == 0 && continue # leaf
      push!(items, (base+key, r, depth+1))
    end
  end
  trie
end

function append!(trie::DATrie, nodeid::Int, keyid::Int, keys::Vector{Int})
  nonempty = 0
  base = trie.base0
  nodes = trie.nodes
  while true
    while base + keys[end] > length(nodes)
      resize!(trie, length(nodes)*2)
    end
    all(k -> !isused(nodes[base+k]), keys) && break
    isused(nodes[base+keys[1]]) && (nonempty += 1)
    base += 1
  end
  nodes[nodeid] = Node(base, nodes[nodeid].check)
  for k in keys
    nodes[base+k] = Node(keyid, nodeid)
  end
  alpha = nonempty / (base - trie.base0 + 1)
  alpha > 0.9 && (trie.base0 = base)
  base
end

end