text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.0
#   kernelspec:
#     display_name: Julia 1.4.2
#     language: julia
#     name: julia-1.4
# ---

# # Gradient descent (勾配法)
#
# 勾配法の練習として下記の簡単な設定のもとで勾配法を逐次適用することで目的関数の極値問題を解く際の点の動きを可視化する.
#
# $A$ を 2x2 の実対称行列とし, $f=f(x,y)$ を次で定める二変数関数とする.
#
# $$
# f(x,y) =  \begin{bmatrix}x & y\end{bmatrix} A \begin{bmatrix} x \\ y \end{bmatrix}
# $$
#
# このとき $f$ の勾配ベクトルは $\nabla f$ は
#
# $$
# \nabla f= \begin{bmatrix} \frac{\partial f}{\partial x} \\ \frac{\partial f}{\partial y} \end{bmatrix} = A \begin{bmatrix} x \\ y \end{bmatrix}
# $$
#
# と書ける. 初期位置を $v_0\in\mathbb{R}^2$, 学習率を $\eta$ とすれば $i$ 番目のステップで得られる点は次のようになる.
#
# $$
# v_{i} = v_{i-1} - \eta \nabla f|_{(x,y)=v_{i-1}} \quad (i=1,2,\dots)
# $$

using Colors
using LinearAlgebra
using Plots

# +
A = [
    2.0  0.2
    0.2  1.0
]

f(v) =  dot(v,A,v)
f(x, y) = f([x, y])
∇f(v) = A * v
η = 0.05 # learning rate
# -

# ## 等高線の描画

x = range(-10,stop=10,length=50)
y = range(-10,stop=10,length=50)
p = contour(x, y, f)

# ## 勾配法を試す

# +
x = range(-5, stop=5, length=50)
y = range(-5, stop=5, length=50)
p = contour(x, y, f, levels=range(0, stop=100, step=3))

v = [-5., -5.] # initial point
for i in 1:100
    v .-= η .* ∇f(v)
    scatter!(
        p, [v[1]], [v[2]],
        label=:none, color=:blue, msize=3,
    )
end

plot(p)
# -

# ## ランダムに点を配置しそれを初期値とする勾配法

# +
x = range(-5, stop=5, length=50)
y = range(-5, stop=5, length=50)
p = contour(x, y, f, levels=range(0,stop=100,step=3))
for _ in 1:10
    v = rand(range(-5,stop=5,step=0.01), 2)
    c =　RGB(rand(3)...)
    for _ in 1:30 # num steps
        v .-= η .* ∇f(v)
        scatter!(
            p, [v[1]], [v[2]],
            label=:none, color=c, msize=3,
        )
    end
end

plot(p)


text:
using MonteCarloMeasurements
using Test, LinearAlgebra, Statistics, Random
import MonteCarloMeasurements: ±, ∓
using MonteCarloMeasurements: nakedtypeof, build_container, build_mutable_container, has_particles, particle_paths
using ControlSystems, Test, GenericLinearAlgebra
ControlSystems.TransferFunction(matrix::Array{<:ControlSystems.SisoRational,2}, Ts, ::Int64, ::Int64) = TransferFunction(matrix,Ts)

@testset "deconstruct" begin
    @info "Testing deconstruct"
    unsafe_comparisons()
    N = 50
    P = tf(1 +0.1StaticParticles(N), [1, 1+0.1StaticParticles(N)])
    f = x->c2d(x,0.1)
    w = Workspace(f,P)
    @time Pd = w(P)
    @test !MonteCarloMeasurements.has_mutable_particles(Pd)
    @test MonteCarloMeasurements.has_mutable_particles(MonteCarloMeasurements.build_mutable_container(Pd))
    # See benchmarks below
    # @profiler Workspace(P)

    tt = function (P)
        f = x->c2d(x,0.1)
        w = Workspace(f,P)
        @time Pd = w(P)
    end
    @test_throws MethodError tt(P) # This causes a world-age problem. If this tests suddenly break, it would be nice and we can get rid of the intermediate workspace object.
    tt = function (P)
        f = x->c2d(x,0.1)
        w = Workspace(f,P)
        @time Pd = w(P,true)
    end
    @test tt(P) == Pd == with_workspace(f,P)
    p = 1 ± 0.1
    @test mean_object(p) == pmean(p)
    @test mean_object([p,p]) == pmean.([p,p])
    @test mean_object(P) ≈ tf(tf(1,[1,1])) atol=1e-2

    @test nakedtypeof(P) == TransferFunction
    @test nakedtypeof(typeof(P)) == TransferFunction
    @test typeof(P) == TransferFunction{ControlSystems.Continuous, ControlSystems.SisoRational{StaticParticles{Float64,N}}}
    P2 = build_container(P)
    @test typeof(P2) == TransferFunction{ControlSystems.Continuous, ControlSystems.SisoRational{Float64}}
    @test typeof(build_mutable_container(P)) == TransferFunction{ControlSystems.Continuous, ControlSystems.SisoRational{Particles{Float64,N}}}
    @test has_particles(P)
    @test has_particles(P.matrix)
    @test has_particles(P.matrix[1])
    @test has_particles(P.matrix[1].num)
    @test has_particles(P.matrix[1].num.coeffs)
    @test has_particles(P.matrix[1].num.coeffs[1])

    # P = tf(1 ± 0.1, [1, 1±0.1])
    # @benchmark foreach(i->c2d($(tf(1.,[1., 1])),0.1), 1:N) # 1.7 ms 1.2 Mb

    bP  = bode(P, exp10.(LinRange(-3, log10(10π), 50)))[1] |> vec
    bPd = bode(Pd, exp10.(LinRange(-3, log10(10π), 50)))[1] |> vec

    @test mean(abs2, pmean.(bP) - pmean.(bPd)) < 1e-4

    A = randn(2,2)
    Ap = A .± 0.1
    hfun = A->Matrix(hessenberg(A))
    @test all(ℝⁿ2ℝⁿ_function(hfun, Ap) .≈ Matrix(hessenberg(A)))

    Ap = A .+ 0.1 .* StaticParticles(1)
    @test_nowarn hessenberg(Ap)

    # bodeplot(P, exp10.(LinRange(-3, log10(10π), 50)))
    # bodeplot!(Pd, exp10.(LinRange(-3, log10(10π), 50)), linecolor=:blue)

    let paths = particle_paths(P), P2 = build_container(P), buffersetter = MonteCarloMeasurements.get_buffer_setter(paths)
        Pres = @unsafe build_mutable_container(f(P)) # Auto-created result buffer
        resultsetter = MonteCarloMeasurements.get_result_setter(Pres)
        @test all(1:paths[1][3]) do i
            buffersetter(P,P2,i)
            P.matrix[1].num.coeffs[1].particles[i] == P2.matrix[1].num.coeffs[1] &&
            P.matrix[1].den.coeffs[2].particles[i] == P2.matrix[1].den.coeffs[2]
            P2res = f(P2)
            resultsetter(Pres, P2res, i)
            Pres.matrix[1].num.coeffs[1].particles[i] == P2res.matrix[1].num.coeffs[1] &&
            Pres.matrix[1].den.coeffs[2].particles[i] == P2res.matrix[1].den.coeffs[2]
        end
    end

    @test mean_object(complex(1. ± 0.1, 1.)) isa ComplexF64
    @test mean_object(complex(1. ± 0.1, 1.)) ≈ complex(1,1) atol=1e-3

    Ps = MonteCarloMeasurements.make_scalar(P)
    @test MonteCarloMeasurements.particletypetuple(Ps.matrix[1].num.coeffs[1]) == (Float64,1,Particles)
    @test MonteCarloMeasurements.particletypetuple(MonteCarloMeasurements.restore_scalar(Ps,50).matrix[1].num.coeffs[1]) == (Float64,50,Particles)

    unsafe_comparisons(false)

    N = 50
    P = tf(1 +0.1StaticParticles(N), [1, 1+0.1StaticParticles(N)])
    f = x->c2d(x,0.1)
    res = MonteCarloMeasurements.array_of_structs(f, P)
    @test length(res) == N
    @test res isa Vector{<:TransferFunction}

end

# julia> @benchmark Pd = w(f) # different world age
# BenchmarkTools.Trial:
#   memory estimate:  1.63 MiB
#   allocs estimate:  19178
#   --------------
#   minimum time:     2.101 ms (0.00% GC)
#   median time:      2.199 ms (0.00% GC)
#   mean time:        2.530 ms (10.36% GC)
#   maximum time:     7.969 ms (53.42% GC)
#   --------------
#   samples:          1973
#   evals/sample:     1

# julia> @benchmark Pd = w(f)  # invokelatest
#   BenchmarkTools.Trial:
# memory estimate:  1.64 MiB
# allocs estimate:  19378
# --------------
# minimum time:     2.204 ms (0.00% GC)
# median time:      2.742 ms (0.00% GC)
# mean time:        3.491 ms (13.77% GC)
# maximum time:     17.103 ms (80.96% GC)
# --------------
# samples:          1429
# evals/sample:     1
#
# julia> @benchmark with_workspace($f,$P) # It seems the majority of the time is spent building the workspace object, so invokelatest really isn't that expensive.
# BenchmarkTools.Trial:
#   memory estimate:  7.90 MiB
#   allocs estimate:  148678
#   --------------
#   minimum time:     158.073 ms (0.00% GC)
#   median time:      165.134 ms (0.00% GC)
#   mean time:        165.842 ms (1.75% GC)
#   maximum time:     180.133 ms (4.80% GC)
#   --------------
#   samples:          31
#   evals/sample:     1

ext: jl


text:
"""
`module LinearMapsAA`

Provides `AbstractArray` (actually `AbstractMatrix`)
or "Ann Arbor" version of `LinearMap` objects
"""
module LinearMapsAA

export LinearMapAA, LinearMapAM, LinearMapAO, LinearMapAX

Indexer = AbstractVector{Int}

include("types.jl")
include("multiply.jl")
include("ambiguity.jl")
include("kron.jl")
include("cat.jl")
include("getindex.jl")
include("setindex.jl")
include("block_diag.jl")
include("lm-aa.jl")
include("identity.jl")

end # module


text:
using ArgParse, JLD

include("./src/simulator.jl")

function parallel_plot()
    s = ArgParseSettings(description = "arg parse here for parallel experimentation")
    @add_arg_table! s begin
        "--location"  # data file location
        "--topology", "-t"  # topology name
        "--te", "-a"  # TE algorithm
        "--plottype", "-p"  # type of plot needed
        "--option", "-o"  # if availability, option number, otherwise just placeholder
        "--scale", "-s"  # if throughput, scale number, otherwise just placeholder
        "--plotname", "-n"  # the name of plot
        "--inverselog", "-i"  # inverse log for availability
        "--lotteryticketset", "-l"  # set of lottery tickets for optimization runtime plot
    end
    parsed_args = parse_args(s) # the result is a Dict{String,Any}
    # println("Parsed args:")
    # for (key,val) in parsed_args
    #     println("  $key  =>  $(repr(val))")
    # end
    # println(parsed_args)

    dir = parsed_args["location"]
    topology = parsed_args["topology"]
    AllAlgorithms = split(parsed_args["te"], ",")
    plottype = parsed_args["plottype"]
    option_num = parse(Int64, parsed_args["option"])  # for demand scaling
    scale_num = parse(Float64, parsed_args["scale"])  # for ticket scaling
    plotname = parsed_args["plotname"]
    inverse_log = parse(Bool, parsed_args["inverselog"])
    lotteryticket_set = split(parsed_args["lotteryticketset"], ",")

    ## demand scaling
    scales = []
    availability_AllTraffic = []
    if topology == "B4"
        scales = collect(1.0:0.5:6.0)
        availability_AllTraffic = collect(1:1:30)
        scenarios = collect(1:1:5)
        topology_indexes = [1]
    elseif topology == "IBM"
        scales = collect(1.0:0.5:7.0)
        availability_AllTraffic = collect(1:1:30)
        scenarios = collect(1:1:3)
        topology_indexes = [1]
    end

    ## ticket number scaling
    tickets = []
    tickets_AllTraffic = []
    if topology == "B4"
        tickets = vcat(collect(1:1:9), collect(10:5:80))
        tickets_AllTraffic = collect(1:1:30)
        tickets_scenarios = collect(1:1:5)
        tickets_topology_indexes = [1]
    elseif topology == "IBM"
        tickets = vcat(collect(1:1:9), collect(10:5:90))
        tickets_AllTraffic = collect(1:1:30)
        tickets_scenarios = collect(1:1:3)
        tickets_topology_indexes = [1]
    end

    ## optimization runtime
    if topology == "B4"
        runtime_tickets = [1,10,40,80,90,120]
        runtime_AllTraffic = collect(1:1:30)
        runtime_scenarios = collect(1:1:5)
        runtime_topology_indexes = [1]
    elseif topology == "IBM"
        runtime_tickets = [1,10,40,80,90,120]
        runtime_AllTraffic = collect(1:1:30)
        runtime_scenarios = collect(1:1:3)
        runtime_topology_indexes = [1]
    end

    Scenario_Availability = Dict{String,Array{Float64,4}}()
    Flow_Availability = Dict{String,Array{Float64,4}}()
    Bandwidth_Availability = Dict{String,Array{Float64,4}}()
    conditional_Scenario_Availability = Dict{String,Array{Float64,4}}()
    conditional_Flow_Availability = Dict{String,Array{Float64,4}}()
    conditional_Bandwidth_Availability = Dict{String,Array{Float64,4}}()
    DirectThroughput = Dict{String,Array{Float64,4}}()
    SecureThroughput = Dict{String,Array{Float64,4}}()

    ticket_DirectThroughput = Dict{String,Array{Float64,4}}()
    ticket_SecureThroughput = Dict{String,Array{Float64,4}}()

    Algo_RouterPorts = Dict{String,Array{Float64,4}}()
    Algo_SolveTime = Dict{String,Array{Float64,4}}()

    progress = ProgressMeter.Progress(length(AllAlgorithms)*length(topology_indexes)*length(scenarios)*length(availability_AllTraffic)*length(scales), .1, "Computing $(plottype)-$(plotname) plotting data...", 50)
    for algorithm in AllAlgorithms
        ## jld data aggregation and processing
        if plottype == "availability"
            Scenario_Availability[algorithm] = zeros(length(topology_indexes), length(scenarios), length(availability_AllTraffic), length(scales))
            Flow_Availability[algorithm] = zeros(length(topology_indexes), length(scenarios), length(availability_AllTraffic), length(scales))
            Bandwidth_Availability[algorithm] = zeros(length(topology_indexes), length(scenarios), length(availability_AllTraffic), length(scales))
            conditional_Scenario_Availability[algorithm] = zeros(length(topology_indexes), length(scenarios), length(availability_AllTraffic), length(scales))
            conditional_Flow_Availability[algorithm] = zeros(length(topology_indexes), length(scenarios), length(availability_AllTraffic), length(scales))
            conditional_Bandwidth_Availability[algorithm] = zeros(length(topology_indexes), length(scenarios), length(availability_AllTraffic), length(scales))
            DirectThroughput[algorithm] = zeros(length(topology_indexes), length(scenarios), length(availability_AllTraffic), length(scales))
            SecureThroughput[algorithm] = zeros(length(topology_indexes), length(scenarios), length(availability_AllTraffic), length(scales))

            for tid in 1:length(topology_indexes)
                topology_index = topology_indexes[tid]
                for sid in 1:length(scenarios)
                    scenario_id = scenarios[sid]
                    for t in 1:length(availability_AllTraffic)
                        traffic_id = availability_AllTraffic[t]
                        for s in 1:length(scales)
                            scale_id = scales[s]

                            if plotname == "scenario_availability" || plotname == "all"
                                Scenario_Availability_FileName = "$(dir)/$(topology)/jld/$(algorithm)/Scenario_Availability_$(topology_index)_$(scenario_id)_$(traffic_id)_$(scale_id)_$(option_num).jld"
                                if isfile(Scenario_Availability_FileName) == true
                                    current_Scenario_Availability = load(Scenario_Availability_FileName, "Scenario_Availability")
                                    Scenario_Availability[algorithm][tid, sid, t, s] = current_Scenario_Availability[algorithm][1,1]  # assume parallel for single traffic and single scale
                                else
                                    printstyled("Missing datafile $(Scenario_Availability_FileName)\n", color=:yellow)
                                    Scenario_Availability[algorithm][tid, sid, t, s] = -1
                                end
                            end

                            if plotname == "flow_availability" || plotname == "all"
                                Flow_Availability_FileName = "$(dir)/$(topology)/jld/$(algorithm)/Flow_Availability_$(topology_index)_$(scenario_id)_$(traffic_id)_$(scale_id)_$(option_num).jld"
                                if isfile(Flow_Availability_FileName) == true
                                    current_Flow_Availability = load(Flow_Availability_FileName, "Flow_Availability")
                                    Flow_Availability[algorithm][tid, sid, t, s] = current_Flow_Availability[algorithm][1,1]  # assume parallel for single traffic and single scale
                                else
                                    printstyled("Missing datafile $(Flow_Availability_FileName)\n", color=:yellow)
                                    Flow_Availability[algorithm][tid, sid, t, s] = -1
                                end
                            end

                            if plotname == "bandwidth_availability" || plotname == "all"
                                Bandwidth_Availability_FileName = "$(dir)/$(topology)/jld/$(algorithm)/Bandwidth_Availability_$(topology_index)_$(scenario_id)_$(traffic_id)_$(scale_id)_$(option_num).jld"
                                if isfile(Bandwidth_Availability_FileName) == true
                                    current_Bandwidth_Availability = load(Bandwidth_Availability_FileName, "Bandwidth_Availability")
                                    Bandwidth_Availability[algorithm][tid, sid, t, s] = current_Bandwidth_Availability[algorithm][1,1]  # assume parallel for single traffic and single scale
                                else
                                    printstyled("Missing datafile $(Bandwidth_Availability_FileName)\n", color=:yellow)
                                    Bandwidth_Availability[algorithm][tid, sid, t, s] = -1
                                end
                            end

                            if plotname == "conditional_scenario_availability" || plotname == "all"
                                conditional_Scenario_Availability_FileName = "$(dir)/$(topology)/jld/$(algorithm)/conditional_Scenario_Availability_$(topology_index)_$(scenario_id)_$(traffic_id)_$(scale_id)_$(option_num).jld"
                                if isfile(conditional_Scenario_Availability_FileName) == true
                                    conditional_current_Scenario_Availability = load(conditional_Scenario_Availability_FileName, "conditional_Scenario_Availability")
                                    conditional_Scenario_Availability[algorithm][tid, sid, t, s] = conditional_current_Scenario_Availability[algorithm][1,1]  # assume parallel for single traffic and single scale
                                else
                                    printstyled("Missing datafile $(conditional_Scenario_Availability_FileName)\n", color=:yellow)
                                    conditional_Scenario_Availability[algorithm][tid, sid, t, s] = -1
                                end
                            end

                            if plotname == "conditional_flow_availability" || plotname == "all"
                                conditional_Flow_Availability_FileName = "$(dir)/$(topology)/jld/$(algorithm)/conditional_Flow_Availability_$(topology_index)_$(scenario_id)_$(traffic_id)_$(scale_id)_$(option_num).jld"
                                if isfile(conditional_Flow_Availability_FileName) == true
                                    conditional_current_Flow_Availability = load(conditional_Flow_Availability_FileName, "conditional_Flow_Availability")
                                    conditional_Flow_Availability[algorithm][tid, sid, t, s] = conditional_current_Flow_Availability[algorithm][1,1]  # assume parallel for single traffic and single scale
                                else
                                    printstyled("Missing datafile $(conditional_Flow_Availability_FileName)\n", color=:yellow)
                                    conditional_Flow_Availability[algorithm][tid, sid, t, s] = -1
                                end
                            end

                            if plotname == "conditional_bandwidth_availability" || plotname == "all"
                                conditional_Bandwidth_Availability_FileName = "$(dir)/$(topology)/jld/$(algorithm)/conditional_Bandwidth_Availability_$(topology_index)_$(scenario_id)_$(traffic_id)_$(scale_id)_$(option_num).jld"
                                if isfile(conditional_Bandwidth_Availability_FileName) == true
                                    conditional_current_Bandwidth_Availability = load(conditional_Bandwidth_Availability_FileName, "conditional_Bandwidth_Availability")
                                    conditional_Bandwidth_Availability[algorithm][tid, sid, t, s] = conditional_current_Bandwidth_Availability[algorithm][1,1]  # assume parallel for single traffic and single scale
                                else
                                    printstyled("Missing datafile $(conditional_Bandwidth_Availability_FileName)\n", color=:yellow)
                                    conditional_Bandwidth_Availability[algorithm][tid, sid, t, s] = -1
                                end
                            end

                            if plotname == "direct_throughput" || plotname == "all"
                                DirectThroughput_FileName = "$(dir)/$(topology)/jld/$(algorithm)/Direct_Throughput_$(topology_index)_$(scenario_id)_$(traffic_id)_$(scale_id)_$(option_num).jld"
                                if isfile(DirectThroughput_FileName) == true
                                    current_DirectThroughput = load(DirectThroughput_FileName, "DirectThroughput")
                                    DirectThroughput[algorithm][tid, sid, t, s] = current_DirectThroughput[algorithm][1,1]  # assume parallel for single traffic and single scale
                                else
                                    printstyled("Missing datafile $(DirectThroughput_FileName)\n", color=:yellow)
                                    DirectThroughput[algorithm][tid, sid, t, s] = -1
                                end
                            end

                            if plotname == "secure_throughput" || plotname == "all"
                                SecureThroughput_FileName = "$(dir)/$(topology)/jld/$(algorithm)/Secure_Throughput_$(topology_index)_$(scenario_id)_$(traffic_id)_$(scale_id)_$(option_num).jld"
                                if isfile(SecureThroughput_FileName) == true
                                    current_SecureThroughput = load(SecureThroughput_FileName, "SecureThroughput")
                                    SecureThroughput[algorithm][tid, sid, t, s] = current_SecureThroughput[algorithm][1,1]  # assume parallel for single traffic and single scale
                                else
                                    printstyled("Missing datafile $(SecureThroughput_FileName)\n", color=:yellow)
                                    SecureThroughput[algorithm][tid, sid, t, s] = -1
                                end
                            end

                            ProgressMeter.next!(progress, showvalues = [])
                        end
                    end
                end
            end
        end

        if plottype == "throughput"
            ticket_DirectThroughput[algorithm] = zeros(length(tickets_topology_indexes), length(tickets_scenarios), length(tickets_AllTraffic), length(tickets))
            ticket_SecureThroughput[algorithm] = zeros(length(tickets_topology_indexes), length(tickets_scenarios), length(tickets_AllTraffic), length(tickets))
            if algorithm == "ARROW" || algorithm == "ARROW_BIN"
                for tid in 1:length(tickets_topology_indexes)
                    topology_index = tickets_topology_indexes[tid]
                    for sid in 1:length(tickets_scenarios)
                        scenario_id = tickets_scenarios[sid]
                        for t in 1:length(tickets_AllTraffic)
                            traffic_id = tickets_AllTraffic[t]
                            for k in 1:length(tickets)
                                ticket_id = tickets[k]
                                ticket_DirectThroughput_FileName = "$(dir)/$(topology)/jld/$(algorithm)/Direct_Throughput_$(topology_index)_$(scenario_id)_$(traffic_id)_$(scale_num)_$(ticket_id).jld"
                                if isfile(ticket_DirectThroughput_FileName) == true
                                    ticket_current_DirectThroughput = load(ticket_DirectThroughput_FileName, "DirectThroughput")
                                    ticket_DirectThroughput[algorithm][tid, sid, t, k] = ticket_current_DirectThroughput[algorithm][1,1]  # assume parallel for single traffic and single ticket num
                                else
                                    printstyled("Missing datafile $(ticket_DirectThroughput_FileName)\n", color=:yellow)
                                    ticket_DirectThroughput[algorithm][tid, sid, t, k] = -1
                                end

                                ticket_SecureThroughput_FileName = "$(dir)/$(topology)/jld/$(algorithm)/Secure_Throughput_$(topology_index)_$(scenario_id)_$(traffic_id)_$(scale_num)_$(ticket_id).jld"
                                if isfile(ticket_SecureThroughput_FileName) == true
                                    ticket_current_SecureThroughput = load(ticket_SecureThroughput_FileName, "SecureThroughput")
                                    ticket_SecureThroughput[algorithm][tid, sid, t, k] = ticket_current_SecureThroughput[algorithm][1,1]  # assume parallel for single traffic and single ticket num
                                else
                                    printstyled("Missing datafile $(ticket_SecureThroughput_FileName)\n", color=:yellow)
                                    ticket_SecureThroughput[algorithm][tid, sid, t, k] = -1
                                end
                                ProgressMeter.next!(progress, showvalues = [])
                            end
                        end
                    end
                end
            else  # FFC,TEAVAR,ARROW_NAIVE
                for tid in 1:length(tickets_topology_indexes)
                    topology_index = tickets_topology_indexes[tid]
                    for sid in 1:length(tickets_scenarios)
                        scenario_id = tickets_scenarios[sid]
                        for t in 1:length(tickets_AllTraffic)
                            traffic_id = tickets_AllTraffic[t]
                            for k in 1:length(tickets)
                                ticket_id = 1
                                ticket_DirectThroughput_FileName = "$(dir)/$(topology)/jld/$(algorithm)/Direct_Throughput_$(topology_index)_$(scenario_id)_$(traffic_id)_$(scale_num)_$(ticket_id).jld"
                                if isfile(ticket_DirectThroughput_FileName) == true
                                    ticket_current_DirectThroughput = load(ticket_DirectThroughput_FileName, "DirectThroughput")
                                    ticket_DirectThroughput[algorithm][tid, sid, t, k] = ticket_current_DirectThroughput[algorithm][1,1]  # assume parallel for single traffic and single ticket num
                                else
                                    printstyled("Missing datafile $(ticket_DirectThroughput_FileName)\n", color=:yellow)
                                    ticket_DirectThroughput[algorithm][tid, sid, t, k] = -1
                                end

                                ticket_SecureThroughput_FileName = "$(dir)/$(topology)/jld/$(algorithm)/Secure_Throughput_$(topology_index)_$(scenario_id)_$(traffic_id)_$(scale_num)_$(ticket_id).jld"
                                if isfile(ticket_SecureThroughput_FileName) == true
                                    ticket_current_SecureThroughput = load(ticket_SecureThroughput_FileName, "SecureThroughput")
                                    ticket_SecureThroughput[algorithm][tid, sid, t, k] = ticket_current_SecureThroughput[algorithm][1,1]  # assume parallel for single traffic and single ticket num
                                else
                                    printstyled("Missing datafile $(ticket_SecureThroughput_FileName)\n", color=:yellow)
                                    ticket_SecureThroughput[algorithm][tid, sid, t, k] = -1
                                end
                                ProgressMeter.next!(progress, showvalues = [])
                            end
                        end
                    end
                end
            end
        end

        if plottype == "routerports"
            Algo_RouterPorts[algorithm] = zeros(length(topology_indexes), length(scenarios), length(availability_AllTraffic), length(scales))
            SecureThroughput[algorithm] = zeros(length(topology_indexes), length(scenarios), length(availability_AllTraffic), length(scales))

            for tid in 1:length(topology_indexes)
                topology_index = topology_indexes[tid]
                for sid in 1:length(scenarios)
                    scenario_id = scenarios[sid]
                    for t in 1:length(availability_AllTraffic)
                        traffic_id = availability_AllTraffic[t]
                        for s in 1:length(scales)
                            scale_id = scales[s]
                            Algo_RouterPorts_Filename = "$(dir)/$(topology)/jld/$(algorithm)/Algo_RouterPorts_$(topology_index)_$(scenario_id)_$(traffic_id)_$(scale_id)_$(option_num).jld"
                            if isfile(Algo_RouterPorts_Filename) == true
                                current_Algo_RouterPorts = load(Algo_RouterPorts_Filename, "Algo_RouterPorts")
                                Algo_RouterPorts[algorithm][tid, sid, t, s] = sum(current_Algo_RouterPorts[algorithm][1,1,:])  # assume parallel for single traffic and single scale
                                # println(Algo_RouterPorts[algorithm][tid, sid, t, s])
                            else
                                printstyled("Missing datafile $(Algo_RouterPorts_Filename)\n", color=:yellow)
                                Algo_RouterPorts[algorithm][tid, sid, t, s] = -1
                                # println(Algo_RouterPorts[algorithm][tid, sid, t, s])
                            end

                            SecureThroughput_FileName = "$(dir)/$(topology)/jld/$(algorithm)/Secure_Throughput_$(topology_index)_$(scenario_id)_$(traffic_id)_$(scale_id)_$(option_num).jld"
                            if isfile(SecureThroughput_FileName) == true
                                current_SecureThroughput = load(SecureThroughput_FileName, "SecureThroughput")
                                SecureThroughput[algorithm][tid, sid, t, s] = current_SecureThroughput[algorithm][1,1]  # assume parallel for single traffic and single scale
                            else
                                printstyled("Missing datafile $(SecureThroughput_FileName)\n", color=:yellow)
                                SecureThroughput[algorithm][tid, sid, t, s] = -1
                            end
                            ProgressMeter.next!(progress, showvalues = [])
                        end
                    end
                end
            end
        end

        if plottype == "optimizationtime"
            Algo_SolveTime[algorithm] = zeros(length(topology_indexes), length(scenarios), length(availability_AllTraffic), length(scales))

            for tid in 1:length(runtime_topology_indexes)
                topology_index = runtime_topology_indexes[tid]
                for sid in 1:length(runtime_scenarios)
                    scenario_id = runtime_scenarios[sid]
                    for t in 1:length(runtime_AllTraffic)
                        traffic_id = runtime_AllTraffic[t]
                        for k in 1:length(runtime_tickets)
                            runtime_option_num = runtime_tickets[k]
                            Algo_Runtime_Filename = "$(dir)/$(topology)/jld/$(algorithm)/Algo_Runtime_$(topology_index)_$(scenario_id)_$(traffic_id)_$(scale_num)_$(runtime_option_num).jld"
                            if isfile(Algo_Runtime_Filename) == true
                                current_Algo_SolveTime = load(Algo_Runtime_Filename, "Algo_Runtime")
                                Algo_SolveTime[algorithm][tid, sid, t, k] = current_Algo_SolveTime[algorithm][1,1]  # assume parallel for single traffic and single scale
                            else
                                printstyled("Missing datafile $(Algo_Runtime_Filename)\n", color=:yellow)
                                Algo_SolveTime[algorithm][tid, sid, t, k] = -1
                            end
                        end
                    end
                end
            end
        end
    end

    ## plotting
    if plottype == "availability"
        if plotname == "scenario_availability" || plotname == "all"
            println("plotting $(plotname)")
            xname = "Demand scales"
            yname = "Scenario availability"
            figname = "$(dir)/$(topology)/scenario_availability.png"
            figname2 = "$(dir)/$(topology)/scenario_availability_ribbon.png"
            line_plot(scales, Scenario_Availability, xname, yname, figname, figname2, AllAlgorithms, false, true, inverse_log, false)
            figname_zoom = "$(dir)/$(topology)/scenario_availability_zoom.png"
            figname_zoom2 = "$(dir)/$(topology)/scenario_availability_zoom_ribbon.png"
            line_plot(scales, Scenario_Availability, xname, yname, figname_zoom, figname_zoom2, AllAlgorithms, true, true, inverse_log, false)
            figname_med = "$(dir)/$(topology)/scenario_availability_med.png"
            figname2_med = "$(dir)/$(topology)/scenario_availability_ribbon_med.png"
            line_plot(scales, Scenario_Availability, xname, yname, figname_med, figname2_med, AllAlgorithms, false, false, inverse_log, false)
            figname_zoom_med = "$(dir)/$(topology)/scenario_availability_zoom_med.png"
            figname_zoom2_med = "$(dir)/$(topology)/scenario_availability_zoom_ribbon_med.png"
            line_plot(scales, Scenario_Availability, xname, yname, figname_zoom_med, figname_zoom2_med, AllAlgorithms, true, false, inverse_log, false)
        end

        if plotname == "flow_availability" || plotname == "all"
            println("plotting $(plotname)")
            xname = "Demand scales"
            yname = "Flow availability"
            figname = "$(dir)/$(topology)/flow_availability.png"
            figname2 = "$(dir)/$(topology)/flow_availability_ribbon.png"
            line_plot(scales, Flow_Availability, xname, yname, figname, figname2, AllAlgorithms, false, true, inverse_log, false)
            figname_zoom = "$(dir)/$(topology)/flow_availability_zoom.png"
            figname_zoom2 = "$(dir)/$(topology)/flow_availability_zoom_ribbon.png"
            line_plot(scales, Flow_Availability, xname, yname, figname_zoom, figname_zoom2, AllAlgorithms, true, true, inverse_log, false)
            figname_med = "$(dir)/$(topology)/flow_availability_med.png"
            figname2_med = "$(dir)/$(topology)/flow_availability_ribbon_med.png"
            line_plot(scales, Flow_Availability, xname, yname, figname_med, figname2_med, AllAlgorithms, false, false, inverse_log, false)
            figname_zoom_med = "$(dir)/$(topology)/flow_availability_zoom_med.png"
            figname_zoom2_med = "$(dir)/$(topology)/flow_availability_zoom_ribbon_med.png"
            line_plot(scales, Flow_Availability, xname, yname, figname_zoom_med, figname_zoom2_med, AllAlgorithms, true, false, inverse_log, false)
        end

        if plotname == "bandwidth_availability" || plotname == "all"
            println("plotting $(plotname)")
            xname = "Demand scales"
            yname = "Bandwidth availability"
            figname = "$(dir)/$(topology)/bandwidth_availability.png"
            figname2 = "$(dir)/$(topology)/bandwidth_availability_ribbon.png"
            line_plot(scales, Bandwidth_Availability, xname, yname, figname, figname2, AllAlgorithms, false, true, inverse_log, false)
            figname_zoom = "$(dir)/$(topology)/bandwidth_availability_zoom.png"
            figname_zoom2 = "$(dir)/$(topology)/bandwidth_availability_zoom_ribbon.png"
            line_plot(scales, Bandwidth_Availability, xname, yname, figname_zoom, figname_zoom2, AllAlgorithms, true, true, inverse_log, false)
            figname_med = "$(dir)/$(topology)/bandwidth_availability_med.png"
            figname2_med = "$(dir)/$(topology)/bandwidth_availability_ribbon_med.png"
            line_plot(scales, Bandwidth_Availability, xname, yname, figname_med, figname2_med, AllAlgorithms, false, false, inverse_log, false)
            figname_zoom_med = "$(dir)/$(topology)/bandwidth_availability_zoom_med.png"
            figname_zoom2_med = "$(dir)/$(topology)/bandwidth_availability_zoom_ribbon_med.png"
            line_plot(scales, Bandwidth_Availability, xname, yname, figname_zoom_med, figname_zoom2_med, AllAlgorithms, true, false, inverse_log, false)
        end

        if plotname == "conditional_scenario_availability" || plotname == "all"
            println("plotting $(plotname)")
            xname = "Demand scales"
            yname = "conditional_Scenario availability"
            figname = "$(dir)/$(topology)/conditional_scenario_availability.png"
            figname2 = "$(dir)/$(topology)/conditional_scenario_availability_ribbon.png"
            line_plot(scales, conditional_Scenario_Availability, xname, yname, figname, figname2, AllAlgorithms, false, true, inverse_log, false)
            figname_zoom = "$(dir)/$(topology)/conditional_scenario_availability_zoom.png"
            figname_zoom2 = "$(dir)/$(topology)/conditional_scenario_availability_zoom_ribbon.png"
            line_plot(scales, conditional_Scenario_Availability, xname, yname, figname_zoom, figname_zoom2, AllAlgorithms, true, true, inverse_log, false)
            figname_med = "$(dir)/$(topology)/conditional_scenario_availability_med.png"
            figname2_med = "$(dir)/$(topology)/conditional_scenario_availability_ribbon_med.png"
            line_plot(scales, conditional_Scenario_Availability, xname, yname, figname_med, figname2_med, AllAlgorithms, false, false, inverse_log, false)
            figname_zoom_med = "$(dir)/$(topology)/conditional_scenario_availability_zoom_med.png"
            figname_zoom2_med = "$(dir)/$(topology)/conditional_scenario_availability_zoom_ribbon_med.png"
            line_plot(scales, conditional_Scenario_Availability, xname, yname, figname_zoom_med, figname_zoom2_med, AllAlgorithms, true, false, inverse_log, false)
        end

        if plotname == "conditional_flow_availability" || plotname == "all"
            println("plotting $(plotname)")
            xname = "Demand scales"
            yname = "conditional_Flow availability"
            figname = "$(dir)/$(topology)/conditional_flow_availability.png"
            figname2 = "$(dir)/$(topology)/conditional_flow_availability_ribbon.png"
            line_plot(scales, conditional_Flow_Availability, xname, yname, figname, figname2, AllAlgorithms, false, true, inverse_log, false)
            figname_zoom = "$(dir)/$(topology)/conditional_flow_availability_zoom.png"
            figname_zoom2 = "$(dir)/$(topology)/conditional_flow_availability_zoom_ribbon.png"
            line_plot(scales, conditional_Flow_Availability, xname, yname, figname_zoom, figname_zoom2, AllAlgorithms, true, true, inverse_log, false)
            figname_med = "$(dir)/$(topology)/conditional_flow_availability_med.png"
            figname2_med = "$(dir)/$(topology)/conditional_flow_availability_ribbon_med.png"
            line_plot(scales, conditional_Flow_Availability, xname, yname, figname_med, figname2_med, AllAlgorithms, false, false, inverse_log, false)
            figname_zoom_med = "$(dir)/$(topology)/conditional_flow_availability_zoom_med.png"
            figname_zoom2_med = "$(dir)/$(topology)/conditional_flow_availability_zoom_ribbon_med.png"
            line_plot(scales, conditional_Flow_Availability, xname, yname, figname_zoom_med, figname_zoom2_med, AllAlgorithms, true, false, inverse_log, false)
        end

        if plotname == "conditional_bandwidth_availability" || plotname == "all"
            println("plotting $(plotname)")
            xname = "Demand scales"
            yname = "conditional Bandwidth availability"
            figname = "$(dir)/$(topology)/conditional_bandwidth_availability.png"
            figname2 = "$(dir)/$(topology)/conditional_bandwidth_availability_ribbon.png"
            line_plot(scales, conditional_Bandwidth_Availability, xname, yname, figname, figname2, AllAlgorithms, false, true, inverse_log, false)
            figname_zoom = "$(dir)/$(topology)/conditional_bandwidth_availability_zoom.png"
            figname_zoom2 = "$(dir)/$(topology)/conditional_bandwidth_availability_zoom_ribbon.png"
            line_plot(scales, conditional_Bandwidth_Availability, xname, yname, figname_zoom, figname_zoom2, AllAlgorithms, true, true, inverse_log, false)
            figname_med = "$(dir)/$(topology)/conditional_bandwidth_availability_med.png"
            figname2_med = "$(dir)/$(topology)/conditional_bandwidth_availability_ribbon_med.png"
            line_plot(scales, conditional_Bandwidth_Availability, xname, yname, figname_med, figname2_med, AllAlgorithms, false, false, inverse_log, false)
            figname_zoom_med = "$(dir)/$(topology)/conditional_bandwidth_availability_zoom_med.png"
            figname_zoom2_med = "$(dir)/$(topology)/conditional_bandwidth_availability_zoom_ribbon_med.png"
            line_plot(scales, conditional_Bandwidth_Availability, xname, yname, figname_zoom_med, figname_zoom2_med, AllAlgorithms, true, false, inverse_log, false)
        end

        if plotname == "direct_throughput" || plotname == "all"
            println("plotting $(plotname)")
            xname = "Demand scales"
            yname = "Direct throughput"
            figname = "$(dir)/$(topology)/DirectThroughput.png"
            figname2 = "$(dir)/$(topology)/DirectThroughput_ribbon.png"
            line_plot(scales, DirectThroughput, xname, yname, figname, figname2, AllAlgorithms, false, true, false, true)
            figname_med = "$(dir)/$(topology)/DirectThroughput_med.png"
            figname2_med = "$(dir)/$(topology)/DirectThroughput_ribbon_med.png"
            line_plot(scales, DirectThroughput, xname, yname, figname_med, figname2_med, AllAlgorithms, false, false, false, true)
        end

        if plotname == "secure_throughput" || plotname == "all"
            println("plotting $(plotname)")
            xname = "Demand scales"
            yname = "Guaranteed throughput under Availability"
            figname = "$(dir)/$(topology)/SecureThroughput.png"
            figname2 = "$(dir)/$(topology)/SecureThroughput_ribbon.png"
            line_plot(scales, SecureThroughput, xname, yname, figname, figname2, AllAlgorithms, false, true, false, true)
            figname_med = "$(dir)/$(topology)/SecureThroughput_med.png"
            figname2_med = "$(dir)/$(topology)/SecureThroughput_ribbon_med.png"
            line_plot(scales, SecureThroughput, xname, yname, figname_med, figname2_med, AllAlgorithms, false, false, false, true)
        end

        open("$(dir)/$(topology)/Availability.txt", "w+") do io
            writedlm(io, ("ScenarioAvailability",))
            for alg in AllAlgorithms
                writedlm(io, (alg, Scenario_Availability[alg]))
            end
            writedlm(io, ("FlowAvailability",))
            for alg in AllAlgorithms
                writedlm(io, (alg, Flow_Availability[alg]))
            end
            writedlm(io, ("BandwidthAvailability",))
            for alg in AllAlgorithms
                writedlm(io, (alg, Bandwidth_Availability[alg]))
            end
        end

        if plotname == "direct_throughput" || plotname == "secure_throughput" || plotname == "all"
            open("$(dir)/$(topology)/Throughput.txt", "w+") do io
                writedlm(io, ("Guaranteed Throughput",))
                for alg in AllAlgorithms
                    writedlm(io, (alg, SecureThroughput[alg]))
                end
                writedlm(io, ("Direct Throughput",))
                for alg in AllAlgorithms
                    writedlm(io, (alg, DirectThroughput[alg]))
                end
            end
        end
    end

    if plottype == "throughput"
        xname = "Number of Lottery Tickets"
        yname = "Guaranteed throughput under Availability"
        figname = "$(dir)/$(topology)/SecureThroughput.png"
        figname2 = "$(dir)/$(topology)/SecureThroughput_ribbon.png"
        line_plot(tickets, ticket_SecureThroughput, xname, yname, figname, figname2, AllAlgorithms, false, true, false, true)
        figname_med = "$(dir)/$(topology)/SecureThroughput_med.png"
        figname2_med = "$(dir)/$(topology)/SecureThroughput_ribbon_med.png"
        line_plot(tickets, ticket_SecureThroughput, xname, yname, figname_med, figname2_med, AllAlgorithms, false, false, false, true)

        xname = "Number of Lottery Tickets"
        yname = "Direct throughput"
        figname = "$(dir)/$(topology)/DirectThroughput.png"
        figname2 = "$(dir)/$(topology)/DirectThroughput_ribbon.png"
        line_plot(tickets, ticket_DirectThroughput, xname, yname, figname, figname2, AllAlgorithms, false, true, false, true)
        figname_med = "$(dir)/$(topology)/DirectThroughput_med.png"
        figname2_med = "$(dir)/$(topology)/DirectThroughput_ribbon_med.png"
        line_plot(tickets, ticket_DirectThroughput, xname, yname, figname_med, figname2_med, AllAlgorithms, false, false, false, true)

        open("$(dir)/$(topology)/Throughput.txt", "w+") do io
            writedlm(io, ("Guaranteed Throughput",))
            for alg in AllAlgorithms
                writedlm(io, (alg, ticket_SecureThroughput[alg]))
            end
            writedlm(io, ("Direct Throughput",))
            for alg in AllAlgorithms
                writedlm(io, (alg, ticket_DirectThroughput[alg]))
            end
        end
    end

    if plottype == "routerports"
        colormap = Dict(
            "HYPERTHETICAL" => "black",
            "ARROW" => "red",
            "ARROW_NAIVE" => "blue",
            "ARROW_BIN" => "red",
            "FFC1" => "green",
            "FFC2" => "yellow",
            "TEAVAR" => "purple",
            "ECMP" => "orange",
        )
        all_sum_routerports = zeros(length(AllAlgorithms))
        nbars = length(AllAlgorithms)
        barWidth = 1/(nbars + 1)
        open("$(dir)/$(topology)/routerports.txt", "w+") do io
            writedlm(io, ("Router Ports",))
            for s in 1:length(scales)
                scale_id = scales[s]
                average_routerports = zeros(length(AllAlgorithms))
                for aa in 1:length(AllAlgorithms)
                    algorithm = AllAlgorithms[aa]
                    vectorized = reshape(Algo_RouterPorts[algorithm][:,:,:,s], size(Algo_RouterPorts[algorithm][:,:,:,s],1)*size(Algo_RouterPorts[algorithm][:,:,:,s],2)*size(Algo_RouterPorts[algorithm][:,:,:,s],3))
                    null_position = findall(x->x==-1, vectorized)
                    vectorized = deleteat!(vectorized, null_position)
                    nan_position = findall(x->isnan(x), vectorized)
                    real_vectorized = deleteat!(vectorized, nan_position)
                    average_routerports[aa] = sum(real_vectorized)
                end
                average_routerports /= length(availability_AllTraffic)*length(scenarios)

                PyPlot.clf()
                sum_routerports = zeros(length(AllAlgorithms))
                for aa in 1:length(AllAlgorithms)
                    # calculate secure throughput
                    vectorized = reshape(SecureThroughput[AllAlgorithms[aa]][:,:,:,s], size(SecureThroughput[AllAlgorithms[aa]][:,:,:,s],1)*size(SecureThroughput[AllAlgorithms[aa]][:,:,:,s],2)*size(SecureThroughput[AllAlgorithms[aa]][:,:,:,s],3))
                    null_position = findall(x->x==-1, vectorized)
                    vectorized = deleteat!(vectorized, null_position)
                    nan_position = findall(x->isnan(x), vectorized)
                    real_vectorized = deleteat!(vectorized, nan_position)
                    throughput_avg = round(sum(real_vectorized)/length(real_vectorized), digits=4)
                    if throughput_avg == 0 || isnan(average_routerports[aa]) || isnan(throughput_avg)
                        sum_routerports[aa] = 0
                    else
                        # println("average_routerports[aa] ", average_routerports[aa], " throughput_avg ", throughput_avg)
                        sum_routerports[aa] = average_routerports[aa] / throughput_avg
                    end
                    all_sum_routerports[aa] += sum_routerports[aa]
                    writedlm(io, (scale_id, AllAlgorithms[aa], average_routerports[aa], throughput_avg, sum_routerports[aa]))
                end
                ## Normalized with respect to max router ports
                # max_routerports = maximum(sum_routerports)
                # normalized_sum_routerports = sum_routerports ./ max_routerports

                ## Normalized with respect to hyperthetical TE
                hyperthetical_routerports = sum_routerports[1]  # fix the first one to be hyperthetical
                normalized_sum_routerports = sum_routerports ./ hyperthetical_routerports

                # println("normalized_sum_routerports ", normalized_sum_routerports)
                for aa in 1:length(AllAlgorithms)
                    PyPlot.bar(AllAlgorithms[aa], normalized_sum_routerports[aa], width=barWidth, alpha = 0.8, label=AllAlgorithms[aa], color=colormap[AllAlgorithms[aa]])
                end
                PyPlot.xlabel("TE algorithms")
                PyPlot.ylabel("Router ports (overall subscribed link bw)")
                PyPlot.xticks(rotation=45)
                figname = "$(dir)/$(topology)/router_ports_$(scale_id).png"
                PyPlot.title(topology)
                PyPlot.savefig(figname)
            end
        end

        PyPlot.clf()
        ## Normalized with respect to max router ports
        # max_all_routerports = maximum(all_sum_routerports)
        # normalized_all_sum_routerports = all_sum_routerports ./ max_all_routerports

        ## Normalized with respect to hyperthetical TE
        hyperthetical_all_routerports = all_sum_routerports[1]  # fix the first one to be hyperthetical
        normalized_all_sum_routerports = all_sum_routerports ./ hyperthetical_all_routerports

        println("normalized_all_sum_routerports ", normalized_all_sum_routerports)
        for aa in 1:length(AllAlgorithms)
            PyPlot.bar(AllAlgorithms[aa], normalized_all_sum_routerports[aa], width=barWidth, alpha = 0.8, label=AllAlgorithms[aa], color=colormap[AllAlgorithms[aa]])
        end
        PyPlot.xlabel("TE algorithms")
        PyPlot.ylabel("Router ports (overall subscribed link bw)")
        PyPlot.xticks(rotation=45)
        figname = "$(dir)/$(topology)/router_ports_all.png"
        PyPlot.title(topology)
        PyPlot.savefig(figname)

        open("$(dir)/$(topology)/routerports_agg.txt", "w+") do io
            writedlm(io, ("Normalized Router Ports",))
            for aa in 1:length(AllAlgorithms)
                writedlm(io, (AllAlgorithms[aa], normalized_all_sum_routerports[aa]))
            end
        end
    end

    if plottype == "optimizationtime"
        average_runtime = zeros(length(AllAlgorithms),length(lotteryticket_set))
        for aa in 1:length(AllAlgorithms)
            algorithm = AllAlgorithms[aa]
            for k in 1:length(lotteryticket_set)
                vectorized = reshape(Algo_SolveTime[algorithm][:,:,:,k], size(Algo_SolveTime[algorithm][:,:,:,k],1)*size(Algo_SolveTime[algorithm][:,:,:,k],2)*size(Algo_SolveTime[algorithm][:,:,:,k],3))
                null_position = findall(x->x==-1, vectorized)
                vectorized = deleteat!(vectorized, null_position)
                nan_position = findall(x->isnan(x), vectorized)
                real_vectorized = deleteat!(vectorized, nan_position)
                average_runtime[aa,k] = round(sum(real_vectorized)/length(real_vectorized), digits=4)
            end
        end

        PyPlot.clf()
        nbars = length(AllAlgorithms)
        barWidth = 1/(nbars + 1)
        for aa in 1:length(AllAlgorithms)
            for k in 1:length(lotteryticket_set)
                PyPlot.bar(lotteryticket_set[k], average_runtime[aa,k], width=barWidth, alpha = 0.8, label=lotteryticket_set[k])
            end
        end
        PyPlot.xlabel("Number of Lottery Tickets")
        PyPlot.ylabel("Optimization Solve Time (seconds)")
        PyPlot.xticks(rotation=45)
        figname = "$(dir)/$(topology)/optimization_solve_time.png"
        PyPlot.title(topology)
        PyPlot.savefig(figname)

        open("$(dir)/$(topology)/optimization_solve_time.txt", "w+") do io
            writedlm(io, ("Optimization solve time",))
            for aa in 1:length(AllAlgorithms)
                writedlm(io, (AllAlgorithms[aa], average_runtime[aa]))
            end
        end
    end

    printstyled("\nTE simulation completed! Please find $(topology)'s [$(plottype) ($(plotname))] results in $(dir)\n\n", color=:blue)
end

parallel_plot()

max_stars_repo_name: hipersys-team/arrow

max_issues_repo_name: arrow-bow/arrow

max_forks_repo_name: arrow-bow/arrow


text:
module LabView0mqJlExamples

# using LabView0mqJl
using JSON3
using Colors, TestImages, ImageIO, ImageShow, FileIO
using Statistics
include("./user_functions.jl")
include("./simple_examples.jl")
include("./image_examples.jl")

export show_img, myf1, myf2, myf3, test_bin2nums,numarr_loopback, numarrs_lpbk, test_rgbimg

# from simple_examples
export arr_amplitude, arr_length, arr_mean, arr_median, arr_rms, vec_amplify, vec_attenuate, vec_power

# from image_examples
export get2d_images, send_testimage

end


text:
"""
Renvoie un solution triée par ordre d'atterrissage.

Les contraintes de séparation de deux avions consécutifs sont respectées.

Return:
    Solution
"""
function greedy(instance::Instance)
    n = instance.nb_planes

    targets = [plane.target for plane in instance.planes]

    σ = sortperm(targets)
    sorted_planes = instance.planes[σ]
    targets = targets[σ]

    sep = [get_sep(instance, sorted_planes[i], sorted_planes[i+1])
           for i in 1:n-1]
    x = zeros(Int, n)
    x[1] = targets[1]
    for i=2:n
        # dans l'idéal il faudrait boucler sur TOUS les avions précédant i,
        # car on peut respecter la distance avec i-1 mais pas i-2 par exemple
        x[i] = max(targets[i], x[i-1] + sep[i-1])
    end
    return Solution(instance, sorted_planes, x)
end

function greedy!(sol::Solution; presort=:ARGS)
    n = sol.inst.nb_planes
    initial_sort!(sol, presort=presort, solve=false)

    sep = [get_sep(sol.inst, first_plane, second_plane)
           for (first_plane, second_plane) in zip(sol.planes[1:end-1], sol.planes[2:end])]

    sol.x[1] = sol.planes[1].target
    for i=2:n
        sol.x[i] = max(sol.planes[i].target, sol.x[i-1] + sep[i-1])
    end
    update_costs!(sol)
end


text:
using QuadGK
using FastClosures
"""
	WharmonicR(t, tf, coeffs)
Harmonic correction Hamiltonian in the rotating frame.
"""
function WharmonicR(t, tf, coeffs)
	ε1, ε2, Δ = coeffs[1], coeffs[2], coeffs[3]
	cos2t, sin2t = cos(2*t), sin(2*t)
	g = (ε1*cos2t + ε2*sin2t)*(1 - cos(2*pi*t/tf))
	W = zeros(Float64, 3)
	W[1] = g*cos2t
	W[2] = g*sin2t
	W[3] = Δ/2
	return W
end
"""
	WharmonicI(t, tf, coeffs)
Harmonic correction Hamiltonian in the interaction picture.
"""
function WharmonicI(t, tf, coeffs)
	ε1, ε2, Δ = coeffs[1], coeffs[2], coeffs[3]
	cos2t, sin2t = cos(2*t), sin(2*t)
	g = (ε1*cos2t + ε2*sin2t)*(1 - cos(2*pi*t/tf))
	Λt = Λ(t,tf)
	cosh2Λ, sinh2Λ = cosh(2*Λt), sinh(2*Λt)
	W = zeros(Float64, 3)
	W[1] = g*cos2t*cosh2Λ + (Δ/2 + g)*sinh2Λ
	W[2] = g*sin2t
	W[3] = g*cos2t*sinh2Λ + (Δ/2 + g)*cosh2Λ
	return W
end
"""
	correction_matrix(tf)
Correction matrix for the harmonic correction.
"""
function correction_matrix(tf)
	M = zeros(Float64, 3, 3)
	h(t) = 1 - cos(2*pi*t/tf)

	func11 = @closure t-> h(t)*((1 + cos(4*t))*cosh(2*Λ(t,tf))
								+ 2*cos(2*t)*sinh(2*Λ(t,tf)))/2
	M[1, 1], err = quadgk(func11, 0.0, tf, rtol=1e-10, atol=1e-14)
	func12 = @closure t-> h(t)*(sin(4*t)*cosh(2*Λ(t,tf))
								+ 2*sin(2*t)*sinh(2*Λ(t,tf)))/2
	M[1, 2], err = quadgk(func12, 0.0, tf, rtol=1e-10, atol=1e-14)
	func13 = @closure t-> sinh(2*Λ(t,tf))/2
	M[1, 3], err = quadgk(func13, 0.0, tf, rtol=1e-10, atol=1e-14)
	func21 = t-> h(t)*sin(4*t)/2
	M[2, 1], err = quadgk(func21, 0.0, tf, rtol=1e-10, atol=1e-14)
	func22 = t-> h(t)*(1 - cos(4*t))/2
	M[2, 2], err = quadgk(func22, 0.0, tf, rtol=1e-10, atol=1e-14)
	func31 = @closure t-> h(t)*((1 + cos(4*t))*sinh(2.0*Λ(t,tf))
								+ 2.0*cos(2*t)*cosh(2.0*Λ(t,tf)))/2
	M[3, 1], err = quadgk(func31, 0.0, tf, rtol=1e-10, atol=1e-14)
	func32 = @closure t-> h(t)*(sin(4*t)*sinh(2.0*Λ(t,tf))
								+ 2.0*sin(2*t)*cosh(2.0*Λ(t,tf)))/2
	M[3, 2], err = quadgk(func32, 0.0, tf, rtol=1e-10, atol=1e-14)
	func33 = @closure t-> cosh(2*Λ(t,tf))/2
	M[3, 3], err = quadgk(func33, 0.0, tf, rtol=1e-10, atol=1e-14)

	return M
end
"""
	coefficients_W(Omega, M)
Calculate the coefficients of the harmonic correction.
"""
function coefficients_W(Omega, M)
	ih_vec = Omega
	coeffs = inv(M)*ih_vec
	return coeffs
end
"""
	get_coeffs(M, coeffs, order, tf)
Calculate the coefficients of the Harmonic correction.
"""
function get_coeffs(M, coeffs, order, tf)
	A = @closure (t-> -V_I(t, tf) .- WharmonicI(t, tf, coeffs))

    Ndims = 3
    Omegaf = solveOmega(tf, order, A, comm!, Ndims)

    vec = sum(Omegaf, dims=2)
    coeffs = coefficients_W(vec, M)
    return coeffs
end
"""
	get_coeffs(M, order, tf)
Calculate the coefficients of the Harmonic correction.
"""
function get_coeffs(M, order, tf)
    coeffs = zeros(3)
    for k = 1: order
		coeffs2 = get_coeffs(M, coeffs, order, tf)
        coeffs = coeffs .+ coeffs2
    end
    return coeffs
end

max_stars_repo_name: thalesfr/MagnusCorrection

max_issues_repo_name: thalesfr/MagnusCorrection

max_forks_repo_name: thalesfr/MagnusCorrection


text:
"""
```
immutable Dimension{D}
    power::Rational{Int}
    Dimension(p) = new(p)
    Dimension(t,p) = new(p)
end
```

Description of a dimension. The name of the dimension `D` is a symbol, e.g.
`:Length`, `:Time`, `:Mass`, etc.

`Dimension{D}` objects are collected in a tuple, which is used for the type
parameter `N` of a [`Dimensions{N}`](@ref) object.

The two-argument constructor ignores the first argument and is used only in the
function [`*(::Unitlike,::Unitlike...)`](@ref).
"""
immutable Dimension{D}
    power::Rational{Int}
    Dimension(p) = new(p)
    Dimension(t,p) = new(p)
end

"""
```
immutable Unit{U}
    tens::Int
    power::Rational{Int}
end
```

Description of a physical unit, including powers-of-ten prefixes and powers of
the unit. The name of the unit `U` is a symbol, e.g. `:Meter`, `:Second`,
`:Gram`, etc. `Unit{U}` objects are collected in a tuple, which is used for the
type parameter `N` of a [`Units{N,D}`](@ref) object.
"""
immutable Unit{U}
    tens::Int
    power::Rational{Int}
end

"""
```
abstract Unitlike
```

Abstract container type for units or dimensions, which need similar
manipulations for collecting powers and sorting. This abstract type is probably
not strictly necessary but facilitates code reuse (see
[`*(::Unitlike,::Unitlike...)`](@ref)).
"""
abstract Unitlike

"""
```
abstract DimensionedUnits{D} <: Unitlike
```

Abstract type used in dispatching on the dimension of a unit.
"""
abstract DimensionedUnits{D} <: Unitlike

"""
```
immutable Units{N,D} <: DimensionedUnits{D}
```

Instances of this object represent units, possibly combinations thereof.
Example: the unit `m` is actually a singleton of type
`Units{(Unit{:Meter}(0,1),), typeof(u"𝐋")}`.
After dividing by `s`, a singleton of type
`Units{(Unit{:Meter}(0,1),Unit{:Second}(0,-1)), typeof(u"𝐋")}` is returned.
"""
immutable Units{N,D} <: DimensionedUnits{D} end

"""
```
immutable Dimensions{N} <: Unitlike
```

Instances of this object represent dimensions, possibly combinations thereof.
"""
immutable Dimensions{N} <: Unitlike end

"""
```
abstract DimensionedQuantity{D} <: Number
```

Super-type of [`Unitful.Quantity`](@ref) types. Used in dispatch on quantities
of a particular dimension, without having to specify the units. The type
parameter `D <: ` [`Unitful.Dimensions`](@ref).
"""
abstract DimensionedQuantity{D} <: Number

"""
```
immutable Quantity{T,D,U} <: DimensionedQuantity{D}
```

A quantity, which has dimensions and units specified in the type signature.
The dimensions and units are allowed to be the empty set, in which case a
dimensionless, unitless number results.

The type parameter `T` represents the numeric backing type. The type parameters
`D <: ` [`Unitful.Dimensions`](@ref) and `U <: ` [`Unitful.Units`](@ref).
Of course, the dimensions follow from the units, but the type parameters are
kept separate to permit convenient dispatch on dimensions.
"""
immutable Quantity{T,D,U} <: DimensionedQuantity{D}
    val::T
    Quantity(v::Number) = new(v)
    Quantity(v::Quantity) = convert(Quantity{T,D,U}, v)
end

"""
```
typealias DimensionlessQuantity{T,U} Quantity{T, Dimensions{()},U}
```

Useful for dispatching on [`Unitful.Quantity`](@ref) types that may have units
but no dimensions. (Units with differing power-of-ten prefixes are not canceled
out.)

Example:
```jldoctest
julia> isa(1.0u"mV/V", DimensionlessQuantity)
true
```
"""
typealias DimensionlessQuantity{T,U} Quantity{T, Dimensions{()}, U}

"""
```
@generated function Quantity(x::Number, y::Units)
```

Outer constructor for `Quantity`s. This is a generated function to avoid
determining the dimensions of a given set of units each time a new quantity is
made.
"""
@generated function Quantity(x::Number, y::Units)
    if y == typeof(NoUnits)
        :(x)
    else
        u = y()
        d = dimension(u)
        :(Quantity{typeof(x), typeof($d), typeof($u)}(x))
    end
end


text:
export Euler

"""
    Euler()

This uses the explicit Euler method for prediction, also known as the
tangent predictor.
"""
struct Euler <: AbstractPredictor end
struct EulerCache{T} <: AbstractPredictorCache
    A::Matrix{T}
    b::Vector{T}
end

cache(::Euler, H, x, ẋ, t) = EulerCache(jacobian(H, x, t), dt(H, x, t))

function predict!(
    xnext,
    cache::EulerCache,
    H::HomotopyWithCache,
    x,
    t,
    Δt,
    ẋ,
    ::JacobianMonitor,
)
    @inbounds for i = 1:length(x)
        xnext[i] = x[i] + Δt * ẋ[i]
    end
    nothing
end

order(::EulerCache) = 2


text:
"""
    targetPos(target,tank)

Is a function that returns the coordinates of the target in a coordinate system which has the origin on the tank at the base of the the canon
"""
function targetPos(target::AbstractTarget, tank::Tank)
    rdoelprimprim = [target.ρ; 0; 0]
    Rhull = angle_to_dcm(-deg2rad(tank.hull.Φ), -deg2rad(tank.hull.ψ), deg2rad(tank.hull.χ), :XZY)
    RLos = angle_to_dcm(0, -deg2rad(tank.sight.θ), deg2rad(tank.sight.ξ), :XZY)
    rdoelprim = Rhull*RLos*rdoelprimprim
    rdoel = rdoelprim + [0; tank.rWY; 0]
    return rdoel
    end
"""
    muzzlePos(tank)

Returns the coordinates of the muzzle exit in a coordinate system which has the origin on the tank at the base of the the canon
"""
function muzzlePos(tank::Tank)
    rprimprim = [tank.canon.lw; 0; 0]
    Rhull = angle_to_dcm(-deg2rad(tank.hull.Φ), -deg2rad(tank.hull.ψ), deg2rad(tank.hull.χ), :XZY)
    Rloop = angle_to_dcm(0, -deg2rad(tank.canon.θ), 0, :XZY)
    Rtoren = angle_to_dcm(0, 0, deg2rad(tank.turret.ξ), :XZY)
    rprim = Rhull*Rtoren*Rloop*rprimprim
    r = rprim +[0; tank.rWY; 0]
    return r
    end
"""
    muzzleVel(tank)

Returns the muzzle velocity vector in a coordinate system which has the origin on the tank at the base of the the canon
"""

    function muzzleVel(tank::Tank)
    uprimprim = [tank.canon.u₀; 0; 0]
    Rhull = angle_to_dcm(-deg2rad(tank.hull.Φ), -deg2rad(tank.hull.ψ), deg2rad(tank.hull.χ), :XZY)
    Rloop = angle_to_dcm(0, -deg2rad(tank.canon.θ), 0, :XZY)
    Rtoren = angle_to_dcm(0, 0, deg2rad(tank.turret.ξ), :XZY)
    uprim = Rhull*Rtoren*Rloop*uprimprim
    u =uprim
    return u
    end

"""
    wind(tank,w)

Returns the wind (w) vector in a coordinate system which has the origin on the tank at the base of the the canon
"""

    function wind(tank::Tank, w::Wind)
        w = [-w.cross*sind(tank.turret.ξ+tank.hull.χ)-w.range*cosd(tank.turret.ξ+tank.hull.χ);0.0;w.cross*cosd(tank.turret.ξ+tank.hull.χ)-w.range*sind(tank.turret.ξ+tank.hull.χ)]
    end


text:
#Here we have two different methods of generating arrays of RTN.
using LinearAlgebra

function generateRTN(time::StepRangeLen, mu, noise_number::Int=1)
    # this one is the faster: we generate the time interval after which
    # we switch the RTN value
    timedim = time.len

    temp = Array{Float64}(undef, timedim, noise_number)
    rMu = mu * Float64(time.step) # Rescaled mu

    if rMu > jumpProb
        @warn "The given timestep is too large compared to the switching rate. Some jumps may be missed."
    end

    for n=1:noise_number    # we generate different values of RTN for each indipendent link of the graph

        i::Integer=1    # index of the noise array
        t::Integer=1    # starting time in units of timedim

        s = 1. * sign(1. - rand() * 2.)  # initial RTN value (either 1 or -1 with probability=0.5)
        # Notice that sign(0.) = 0 but here we can't get it since rand() is in [0, 1)

        # we push forward the time of a random value given by t=-floor(log(p0)/rescaledMu)
        # where p0 is a random number in [0,1], so that t is a random variable distributed
        # according to p(t)=Exp(-rescaledMu t) (no switches)

        t += Int(fld(-log(1. - rand()), rMu))

        while t < timedim

            view(temp, i:t, n) .= s      # we fill the noise array until t

            i=t+1
            s=-s                         # then we switch value
            t += 1 + Int(fld(-log(1. - rand()),rMu))   # again we push forward the time
        end

        view(temp, i:timedim, n) .= s        # last filling

    end

    return temp
end

max_stars_repo_name: matteoacrossi/QuantumSpatialSearch

max_issues_repo_name: matteoacrossi/QuantumSpatialSearch

max_forks_repo_name: matteoacrossi/QuantumSpatialSearch


text:
println("\nTesting extraction from Terms\n")

# testing extracting stuff from terms

m,n1 = 5,3
x1 = Variable(n1)
A = randn(m,n1)
# single term, single variable
cf = ls(A*x1)
xAll = StructuredOptimization.extract_variables(cf)
@test xAll[1] == x1
L = StructuredOptimization.extract_operators(xAll,cf)
@test typeof(L) <: MatrixOp
La = StructuredOptimization.extract_affines(xAll,cf)
@test typeof(La) <: MatrixOp
f = StructuredOptimization.extract_functions(cf)
@test typeof(f) <: SqrNormL2

# multiple terms, single variable
b1 = randn(n1)
cf = ls(A*x1) + 2.5*norm(x1+b1,1)
xAll = StructuredOptimization.extract_variables(cf)
@test xAll[1] == x1
V = StructuredOptimization.extract_operators(xAll,cf)
@test typeof(V) <: VCAT
@test typeof(V[1]) <: MatrixOp
@test typeof(V[2]) <: Eye
V2 = StructuredOptimization.extract_affines(xAll,cf)
@test typeof(V2) <: VCAT
@test typeof(V2[1]) <: MatrixOp
@test typeof(V2[2]) <: AffineAdd{T} where {T <: Eye}
f = StructuredOptimization.extract_functions(cf)
@test typeof(f) <: SeparableSum
@test typeof(f.fs[1]) <: SqrNormL2
@test typeof(f.fs[2]) <: Postcompose
x = randn(n1)
@test norm(f.fs[2](x) - 2.5*norm(x+b1,1)) < 1e-12

# single term, multiple variables
x2 = Variable(m)
cf = ls(A*x1+x2+20)
xAll = StructuredOptimization.extract_variables(cf)
xAll = (x2,x1) # change the order on pourpose
H = StructuredOptimization.extract_operators(xAll,cf)
@test typeof(H) <: HCAT
@test typeof(H[1]) <: Eye
@test typeof(H[2]) <: MatrixOp
H2 = StructuredOptimization.extract_affines(xAll,cf)
@test typeof(H2[1]) <: AffineAdd{T} where {T <: Eye}
@test typeof(H2[2]) <: AffineAdd{T} where {T <: MatrixOp}
f = StructuredOptimization.extract_functions(cf)
@test typeof(f) <: PrecomposeDiagonal

### multiple terms, multiple variables
n1,n2,n3,n4,n5 = 3,3,4,4,7
A = randn(n5,n1)
x1,x2,x3,x4,x5 = Variable(randn(n1)),Variable(randn(n2)),Variable(randn(n3)),Variable(randn(n4)),Variable(randn(n5))

cf = ls(x1+x2)
xAll = StructuredOptimization.extract_variables(cf)
@test xAll == (x1,x2)

cf = ls(x1+x2)+ls(x1)
xAll = StructuredOptimization.extract_variables(cf)
@test xAll == (x1,x2)

cf = ls(x1+x2)+ls(x3+x4)+ls(x5)+ls(x5+A*x2)+ls(x1)+ls(x5)
xAll = StructuredOptimization.extract_variables(cf)
@test xAll == (x1,x2,x3,x4,x5)

V = StructuredOptimization.extract_operators(xAll,cf)

@test typeof(V[1][1]) <: Eye
@test typeof(V[1][2]) <: Eye
@test typeof(V[1][3]) <: Zeros
@test typeof(V[1][4]) <: Zeros
@test typeof(V[1][5]) <: Zeros

@test typeof(V[2][1]) <: Zeros
@test typeof(V[2][2]) <: Zeros
@test typeof(V[2][3]) <: Eye
@test typeof(V[2][4]) <: Eye
@test typeof(V[2][5]) <: Zeros

@test typeof(V[3][1]) <: Zeros
@test typeof(V[3][2]) <: Zeros
@test typeof(V[3][3]) <: Zeros
@test typeof(V[3][4]) <: Zeros
@test typeof(V[3][5]) <: Eye

@test typeof(V[4][1]) <: Zeros
@test typeof(V[4][2]) <: MatrixOp
@test typeof(V[4][3]) <: Zeros
@test typeof(V[4][4]) <: Zeros
@test typeof(V[4][5]) <: Eye

@test typeof(V[5][1]) <: Eye
@test typeof(V[5][2]) <: Zeros
@test typeof(V[5][3]) <: Zeros
@test typeof(V[5][4]) <: Zeros
@test typeof(V[5][5]) <: Zeros

@test typeof(V[6][1]) <: Zeros
@test typeof(V[6][2]) <: Zeros
@test typeof(V[6][3]) <: Zeros
@test typeof(V[6][4]) <: Zeros
@test typeof(V[6][5]) <: Eye

println("\nTesting splitting Terms\n")

x = Variable(5)
y = Variable(5)
cf = ls(x)+10*norm(x,2)+ls(x+y)

f, g = StructuredOptimization.split_smooth(cf)
@test f[1] == cf[1]
@test f[2] == cf[3]
@test g[1] == cf[2]

cf = ls(x)
f, g = StructuredOptimization.split_smooth((cf,))
@test f == (cf,)
@test g == ()

cf = norm(x,1)+norm(y,2)+norm(randn(5,5)*x+y,Inf)
xAll = StructuredOptimization.extract_variables(cf)
AAc, nonAAc = StructuredOptimization.split_AAc_diagonal(cf)
@test AAc[1] == cf[1]
@test AAc[2] == cf[2]
@test nonAAc[1] == cf[3]

cf = ls(sigmoid(x)) + ls(x)
fq, fs = StructuredOptimization.split_quadratic(cf)
@test fs[1] == cf[1]
@test fq[1] == cf[2]

println("\nTesting extracting Proximable functions\n")
# testing is_proximable
@test StructuredOptimization.is_proximable(AAc) == true
@test StructuredOptimization.is_proximable(nonAAc) == false

cf = norm(x[1:2],1)+norm(x[3:5])
xAll = StructuredOptimization.extract_variables(cf)

@test all(StructuredOptimization.is_AAc_diagonal.(cf)) == true
@test StructuredOptimization.is_proximable(cf) == true

cf = norm(x[1:2],1)+norm(x[3:5])+norm(x,Inf)
xAll = StructuredOptimization.extract_variables(cf)

@test all(StructuredOptimization.is_AAc_diagonal.(cf)) == true
@test StructuredOptimization.is_proximable(cf) == false

# testing extract_proximable
# single variable, single term
x = Variable(randn(5))
b = randn(5)
cf = 10*norm(x-b,1)
xAll = StructuredOptimization.extract_variables(cf)
@test StructuredOptimization.is_proximable(cf) == true

f = StructuredOptimization.extract_proximable(xAll,cf)
@test norm(f(~x) - 10*norm(~x-b,1)) < 1e-12

# single variable, single term, diagonal term
x = Variable(randn(5))
b = randn(5)
d = randn(5)
cf = 10*norm(d.*x-b,1)
xAll = StructuredOptimization.extract_variables(cf)
@test StructuredOptimization.is_proximable(cf) == true

f = StructuredOptimization.extract_proximable(xAll,cf)
@test norm(f(~x) - 10*norm(d.*~x-b,1)) < 1e-12

# single variable, single term, tight frame term
x = Variable(randn(5))
b = randn(5)
d = randn(5)
cf = 10*norm(dct(x)-b,1)
xAll = StructuredOptimization.extract_variables(cf)
@test StructuredOptimization.is_proximable(cf) == true

f = StructuredOptimization.extract_proximable(xAll,cf)
@test norm(f(~x) - 10*norm(dct(~x)-b,1)) < 1e-12

# single variable, single term, tight frame term, fft
# TODO this not working (probably fix needed in ProxOp)
#x = Variable(randn(5))
#b = randn(5)
#d = randn(5)
#cf = 10*norm(fft(x)-b,1)
#xAll = StructuredOptimization.extract_variables(cf)
#@test StructuredOptimization.is_proximable(cf) == true
#
#f = StructuredOptimization.extract_proximable(xAll,cf)
#@test norm(f(~x) - 10*norm(fft(~x)-b,1)) < 1e-12

# single variable, multiple terms with GetIndex
x = Variable(randn(5))
b = randn(2)
cf = 10*norm(x[1:2]-b,1)+norm(x[3:5],2)
xAll = StructuredOptimization.extract_variables(cf)
@test StructuredOptimization.is_proximable(cf) == true
f = StructuredOptimization.extract_proximable(xAll,cf)
@test norm(f(~x) - sum([10*norm((~x)[1:2]-b,1);norm((~x)[3:5],2)])) < 1e-12

# single variable, multiple terms with GetIndex composed with dct
x = Variable(randn(5))
b = randn(2)
cf = 10*norm(x[1:2]-b,1)+norm(dct(x[3:5]),2)
xAll = StructuredOptimization.extract_variables(cf)
@test StructuredOptimization.is_proximable(cf) == true
f = StructuredOptimization.extract_proximable(xAll,cf)
@test norm(f(~x) - sum([10*norm((~x)[1:2]-b,1);norm(dct((~x)[3:5]),2)])) < 1e-12

# multiple variables, multiple terms
x1 = Variable(randn(5))
b1 = randn(5)
x2 = Variable(randn(3))
b2 = randn(3)

cf = 10*norm(x2-b2,1)+norm(x1+b1,2)
xAll = (x1,x2)
@test StructuredOptimization.is_proximable(cf) == true
f = StructuredOptimization.extract_proximable(xAll,cf)
@test norm(f.fs[1](~x1)-norm(~x1+b1,2) ) < 1e-12
@test norm(f.fs[2](~x2)-10*norm(~x2-b2,1) ) < 1e-12

x1 = Variable(randn(5))
b1 = randn(5)
x2 = Variable(randn(5))
b2 = randn(5)

# TODO fix this?
#cf = 10*norm(x2+x1+b2,1)
#xAll = (x1,x2)
#@test StructuredOptimization.is_proximable(cf) == true
#f = StructuredOptimization.extract_proximable(xAll,cf)
# TODO fix this! in ProxOp?
# @test norm(f((~x1,~x2))-10*norm(~x2+~x1+b2,1) ) < 1e-12

# multiple variables, missing terms
x1 = Variable(randn(5))
b1 = randn(5)
x2 = Variable(randn(3))
b2 = randn(3)

cf = 10*norm(x2-b2,1)
xAll = (x1,x2)
@test StructuredOptimization.is_proximable(cf) == true
f = StructuredOptimization.extract_proximable(xAll,cf)
@test f.fs[1](~x1) == 0.
@test norm(f.fs[2](~x2)-10*norm(~x2-b2,1) ) < 1e-12

# multiple variables, multiple terms, with GetIndex
x1 = Variable(randn(5))
b1 = randn(5)
x2 = Variable(randn(3))
b2 = randn(3)

cf = norm(x1[3:5]+b1[3:5],1)+10*norm(x2-b2,1)+norm(x1[1:2]+b1[1:2],2)
xAll = (x1,x2)
@test StructuredOptimization.is_proximable(cf) == true
f = StructuredOptimization.extract_proximable(xAll,cf)
@test norm(f.fs[1](~x1)-norm((~x1)[1:2]+b1[1:2],2)-norm((~x1)[3:5]+b1[3:5],1) ) < 1e-12
@test norm(f.fs[2](~x2)-10*norm(~x2-b2,1) ) < 1e-12


text:
"""
      Satellite

      Orientation represented with
      Modified Rodrigues Parameters.
"""

struct Satellite{I, T} <: Model{I, T}
      n::Int
      m::Int
      d::Int

      J      # inertia matrix
end

function f(model::Satellite, z, u, w)
      # states
      r = view(z, 1:3)
      ω = view(z, 4:6)

      # controls
      τ = view(u, 1:3)

      SVector{6}([0.25 * ((1.0 - r' * r) * ω - 2.0 * cross(ω, r) + 2.0 * (ω' * r) * r);
                  model.J \ (τ - cross(ω, model.J * ω))])
end

function kinematics(model::Satellite, q)
	p = @SVector [0.1, 0.0, 0.0]
	k = MRP(view(q, 1:3)...) * p
	return k
end

# model
n, m, d = 6, 3, 0
J = Diagonal(@SVector[1.0, 2.0, 3.0])

model = Satellite{Midpoint, FixedTime}(n, m, d, J)

# visuals
function visualize!(vis, p::Satellite, q; Δt = 0.1)
	setvisible!(vis["/Background"], true)
	setprop!(vis["/Background"], "top_color", RGBA(0.0, 0.0, 0.0, 1.0))
	setprop!(vis["/Background"], "bottom_color", RGBA(0.0, 0.0, 0.0, 1.0))
	setvisible!(vis["/Axes"], false)
	setvisible!(vis["/Grid"], false)

    setobject!(vis[:satellite],
    	Rect(Vec(-0.25, -0.25, -0.25),Vec(0.5, 0.5, 0.5)),
    	MeshPhongMaterial(color = RGBA(1.0, 1.0, 1.0, 1.0)))

    arrow_x = ArrowVisualizer(vis[:satellite][:arrow_x])
    mat = MeshPhongMaterial(color=RGBA(1.0, 0.0, 0.0, 1.0))
    setobject!(arrow_x, mat)
    settransform!(arrow_x,
    	Point(0.0, 0.0, 0.0),
    	Vec(0.75, 0.0, 0.0),
    	shaft_radius=0.05,
    	max_head_radius=0.1)

    arrow_y = ArrowVisualizer(vis[:satellite][:arrow_y])
    mat = MeshPhongMaterial(color=RGBA(0.0, 1.0, 0.0, 1.0))
    setobject!(arrow_y, mat)
    settransform!(arrow_y,
    	Point(0.0, 0.0, 0.0),
    	Vec(0.0, 0.75, 0.0),
    	shaft_radius=0.05,
    	max_head_radius=0.1)

    arrow_z = ArrowVisualizer(vis[:satellite][:arrow_z])
    mat = MeshPhongMaterial(color=RGBA(0.0, 0.0, 1.0, 1.0))
    setobject!(arrow_z, mat)
    settransform!(arrow_z,
    	Point(0.0, 0.0, 0.0),
    	Vec(0.0, 0.0, 0.75),
    	shaft_radius=0.05,
    	max_head_radius=0.1)

    anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

     for t = 1:length(q)
    	 MeshCat.atframe(anim, t) do
    		 settransform!(vis["satellite"],
    			   compose(Translation(0.0 * [-0.25; -0.25; -0.25]...),
    					 LinearMap(MRP(q[t][1:3]...))))
    	 end
     end

     MeshCat.setanimation!(vis, anim)
end

max_issues_repo_name: thowell/DirectMotionPlanning


text:
"""
    plotCounts(fcsRaw, md, group_by = "condition")

Barplot showing the number of cells per sample, used as a guide to identify samples where not enough cells were assayed

# Arguments:
- `fcsRaw`: raw FCS data
- `md`: Metadata table
"""

function plotCounts(fcsRaw, md, group_by = "condition")

    df_barplot = DataFrame(filename = String[], size = Int[], condition = String[])

    for (k,v) in fcsRaw
        sid = md.sample_id[k .== md.file_name]
        # println(sid[1])
        condition = md.condition[k .== md.file_name]
        push!(df_barplot, (string(sid[1]), size(v)[1], condition[1]) )
    end
    sort!(df_barplot)
    bar(df_barplot.filename, df_barplot.size, title="Numer of Cells", group=df_barplot.condition,xrotation=60)
end

"""
    plotPCA(daFrame)

Plotting the PCA of all median marker expression

# Arguments:
- `daFrame`: daFrame containing the fcs data, metadata and panel
"""

function plotPCA(daf, md)
    dfall_median = aggregate(daf.fcstable, :sample_id, Statistics.median)

    T = convert(Matrix, dfall_median)
    samples_ids = T[:,1]
    T_reshaped = permutedims(convert(Matrix{Float64}, T[:, 2:10]), [2, 1])

    my_pca = StatsBase.fit(MultivariateStats.PCA, T_reshaped)

    yte = MultivariateStats.transform(my_pca,T_reshaped)

    df_pca = DataFrame(yte')
    df_pca[:sample_id] = samples_ids

    # get the condition per sample id and add in DF
    v1= df_pca.sample_id; v2=md.sample_id
    idxs = indexin(v1, v2)
    df_pca[:condition] = md.condition[idxs]

    StatsPlots.@df df_pca scatter(:x1, :x2, group=:condition)
end


text:
facts("Interface.") do

    context("Module.") do

        @fact isempty(Docile.Interface.documented()) --> false
        @fact typeof(Docile.Interface.documented())  --> Set{Module}

        @fact Docile.Interface.isdocumented(Docile)           --> true
        @fact Docile.Interface.isdocumented(Docile.Interface) --> true

        @fact typeof(Docile.Interface.metadata(Docile)) --> Docile.Interface.Metadata

    end

    context("Metadata.") do

        meta = Docile.Interface.metadata(Docile)

        @fact typeof(Docile.Interface.modulename(meta)) --> Module

        @fact isempty(Docile.Interface.manual(meta)) --> true
        @fact typeof(Docile.Interface.manual(meta))  --> Vector{String}

        @fact isempty(Docile.Interface.entries(meta))                                    --> false
        @fact isempty(Docile.Interface.entries(Docile.Interface.metadata(Docile.Cache))) --> false

        @fact typeof(Docile.Interface.entries(meta)) --> ObjectIdDict

        dir = Pkg.dir("Docile", "src", "Docile.jl")
        @fact Docile.Interface.root(meta)                                    --> dir

        dir = Pkg.dir("Docile", "src", "Cache", "Cache.jl")
        @fact Docile.Interface.root(Docile.Interface.metadata(Docile.Cache)) --> dir

        @fact isempty(Docile.Interface.metadata(meta)) --> false
        @fact typeof(Docile.Interface.metadata(meta))  --> Dict{Symbol, Any}

        legacy = Docile.Interface.metadata(Docile.Legacy)
        s = symbol("@comment")
        λ = getfield(Docile.Legacy, s)
        @fact Docile.Interface.macroname(Docile.Interface.entries(legacy)[λ]) --> s

        @fact Docile.Interface.name(isempty)                 --> :isempty
        @fact Docile.Interface.name(:isempty)                --> :isempty
        @fact Docile.Interface.name(first(methods(isempty))) --> :isempty
        @fact Docile.Interface.name(Docile)                  --> :Docile
        @fact Docile.Interface.name(Docile.Legacy.Metadata)  --> :Metadata

        @fact Docile.Interface.files(meta)                --> Set{String}()
        @fact Docile.Interface.isloaded(meta)             --> true
        @fact Docile.Interface.isexported(Docile, Docile) --> true

        cache = Docile.Interface.metadata(Docile.Cache)
        comment = first(filter(k -> isa(k, Docile.Legacy.Comment), keys(Docile.Interface.entries(cache))))
        @fact Docile.Interface.isexported(Docile.Cache, comment) --> false

    end

    context("Entry.") do

        meta = Docile.Interface.metadata(Docile.Cache)
        doc =  Docile.Interface.entries(meta)[first(methods(Docile.Cache.clear!))]

        @fact Docile.Interface.category(doc) --> :method

        @fact Docile.Interface.modulename(doc) --> Docile.Cache

        @fact isempty(Docile.Interface.metadata(doc)) --> false
        @fact typeof(Docile.Interface.metadata(doc))  --> Dict{Symbol, Any}

        @fact typeof(Docile.Interface.docs(doc)) --> Docile.Interface.Docs{:md}

    end

    context("Docs.") do

        meta = Docile.Interface.metadata(Docile.Cache)
        ent =  Docile.Interface.entries(meta)[first(methods(Docile.Cache.clear!))]
        doc =  Docile.Interface.docs(ent)

        @fact isempty(Docile.Interface.data(doc)) --> false
        @fact isa(Docile.Interface.data(doc), String) --> true

        @fact Docile.Interface.format(doc) --> :md

        @fact Docile.Interface.parsedocs(Docile.Interface.Docs{:txt}("")) --> ""

        @fact_throws Docile.Interface.parsed(doc)

        @fact_throws Docile.Interface.parsedocs(Docile.Interface.Docs{:md}(""))

    end

    context("Deprecated.") do

        @fact typeof(Docile.Interface.documentation(Docile)) --> Docile.Legacy.Metadata

    end

end


text:
module Constants

    export me_g, mp_g, me, qe_C, qe, g, eps0

    const me_g = 9.11e-28#9.109383e-28 #g
    const mp_g = 1.*1836.*9.11e-28#1836.152674 * me_g #proton mass in grams
    const me = 1.e-3 * me_g #kg
    const qe_C = 1.6e-19#1.602177e-19
    const qe = 4.8032e-10#4.803204e-10 #statcoulomb
    const g = 5./3. #gamma
    const eps0 = 8.854188e-12

end

max_stars_repo_name: cbellei/JUPL

max_issues_repo_name: cbellei/JUPL

max_forks_repo_name: cbellei/JUPL


text:
# Julia1_2.5.jl Strings (Unicode and UTF-8)
#
# ASCII is preserved in the first 128
#  characters of UTF-8 so only 1 byte necessary
# A UTF-8 encoding can take 1 to 4 bytes in a UTF-8 String
#  to get the size in bytes:
sizeof("∑")  # ∑ is upper case sigma

s1="abcπ∑d" # UTF-8 String type
println(sizeof(s1))   # Number of bytes

hg=80
s2=string("xy",3,Char(0x3B1),false,s1,hg) # function converts UTF-8
println(s2)
println(sizeof(s2))   # Number of bytes
println(length(s1))   # number of characters
# translate a UTF-8 string to a vector of UInt8,
#  UInt16, UInt32 or vice-versa
s="abcde"
t=transcode(UInt8,s)

a1=ascii(s)  # Will convert to an ascii string or error
#    if non-ascii (used to send strings to a C routine)
a1=ascii(s1) # will report invalid ascii character

println(s1)
a2=transcode(UInt32,s1)  # can make 32-bit Char Vector.
println(a2)
println(sizeof(a2))   # Number of bytes
println(length(a2))   # number of 32-bit character

# Char() can be used to check out unicode characters to see
#    if they are what you want
#  I always use the hexadecimal representations of the Unicode
Char(0x3B3)
'γ'  # to see the hex Unicode for a character
Char(0x3C6)
Char(0x3C9)  # Greek alphabet 0x3B1 through 0x3C9.
Char(0x263A)


text:
function divedivedive()
    global depth
    depth += 1
    divedivedive()
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode