text:
using DelayEmbeddings
using DelimitedFiles
using StatsBase
using Random

include("../../../src/pecuzal_method.jl")
include("../../../src/data_analysis_functions.jl")

## We analyze the reconstruction from standard time delay embedding, the
# MDOP embedding algorithm and the Garcia & Almeida method and compare it to
# our propsed method: The pecuzal method.
# For evaluation of the different reconstructions we consider Uzal's L-statistic,
# the fraction of recurrence rates from the JRP of the reconstruction and the
# reference and only the reference. As a third criterion we look at the
# recurrence time entropy, since it is related to the Kolmogorov-entropy as a
# dynamical invariant. We also look at two other RQA-measures and, moreover, the
# generalized mutual false nearest neighbors.

# Here we look at real data of a chaotic electrochemical oscillator (Ni) from
# the Kiss group with two different resistances R₁ = 1,5kΩ, R₂ = 2,5kΩ leading to
# "more" chaoticity.

## Load the data
Random.seed!(1234)
data1 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/s3ch4.txt")
data2 = readdlm("./scripts/Experimental data/Standard electrochemical chaos/s1ch3.txt")
# length of the whole time series
M1 = length(data1)
M2 = length(data2)
M = minimum(hcat(M1, M2))

# downsampling
data1 = data1[1:2:M]
data2 = data2[1:2:M]
M = length(data1)

data1 = regularize(data1 .+ .00000001*randn(length(data1)))
data2 = regularize(data2 .+ .00000001*randn(length(data2)))

# the possible lags
taus = 0:200

# set recurrence threshold for further analysis (and min line length)
ε = 0.08
lmin = 2

# set length of time series sample
N = 500

# number of time series samples
samplesize = 1

# draw samples
idxs1 = sample(vec(1:M-N-1), samplesize, replace = false)
idxs2 = sample(vec(1:M-N-1), samplesize, replace = false)

# compute reference (full time series length)
println("Start compute reference from full time series")

# w1 = estimate_delay(data1, "mi_min")
# w2 = estimate_delay(data2, "mi_min")
#
# #Standard TDE
# println("Compute reference TDE")
# Y_tde1, τ_tde1, _ = optimal_traditional_de(data1; τs = taus)
# Y_tde2, τ_tde2, _ = optimal_traditional_de(data2; τs = taus)
#
# dim_tde_ref1 = size(Y_tde1,2)
# dim_tde_ref2 = size(Y_tde2,2)
# τ_tdes1 = [(i-1)*τ_tde1 for i = 1:(size(Y_tde1,2))]
# τ_tdes2 = [(i-1)*τ_tde2 for i = 1:(size(Y_tde2,2))]
# Tw_tde_ref1 = sum(τ_tdes1)
# Tw_tde_ref2 = sum(τ_tdes2)
#
# #MDOP embedding
# println("Compute reference MDOP")
# tw1 = mdop_maximum_delay(data1)
# tw2 = mdop_maximum_delay(data2)
# lm1 = DelayEmbeddings.findlocalminima(tw1[2])
# if lm1[1] ≤ 3
#     lmm1 = try lm1[2]
#     catch
#         lm1[1]
#     end
# else
#     lmm1 = lm1[1]
# end
# lm2 = DelayEmbeddings.findlocalminima(tw2[2])
# if lm2[1] ≤ 3
#     lmm2 = try lm2[2]
#     catch
#         lm2[1]
#     end
# else
#     lmm2 = lm2[1]
# end
# tau_max1 = lmm1
# tau_max2 = lmm2
# taus_mdop1 = 0:tau_max1
# taus_mdop2 = 0:tau_max2
#
# Y_mdop1, τ_vals_mdop1, ts_vals_mdop1, FNNs_mdop1 , βs_mdop1 = mdop_embedding(data1;
#                                                     τs = taus_mdop1 , w = w1)
# dim_mdop_ref1 = size(Y_mdop1,2)
# Tw_mdop_ref1 = sum(τ_vals_mdop1)
# Y_mdop2, τ_vals_mdop2, ts_vals_mdop2, FNNs_mdop2 , βs_mdop2 = mdop_embedding(data2;
#                                                     τs = taus_mdop2 , w = w2)
# dim_mdop_ref2 = size(Y_mdop2,2)
# Tw_mdop_ref2 = sum(τ_vals_mdop2)
#
# #Garcia & Almeida
# println("Compute reference GA")
# Y_GA1, τ_vals_GA1, ts_vals_GA1, FNNs_GA1 , ns_GA1 = garcia_almeida_embedding(data1;
#                                                     τs = taus , w = w1, T = w1)
# dim_GA_ref1 = size(Y_GA1,2)
# Tw_GA_ref1 = sum(τ_vals_GA1)
#
# Y_GA2, τ_vals_GA2, ts_vals_GA2, FNNs_GA2 , ns_GA2 = garcia_almeida_embedding(data2;
#                                                     τs = taus , w = w2, T = w2)
# dim_GA_ref2 = size(Y_GA2,2)
# Tw_GA_ref2 = sum(τ_vals_GA2)
#
# #Pecuzal
# println("Compute reference PECUZAL")
# Y_pec1, τ_vals_pec1, ts_vals_pec1, Ls_pec1 , εs_pec1 = pecuzal_embedding_update(data1;
#                                                             τs = taus , w = w1)
# dim_pec_ref1 = size(Y_pec1,2)
# Tw_pec_ref1 = sum(τ_vals_pec1)
#
# Y_pec2, τ_vals_pec2, ts_vals_pec2, Ls_pec2 , εs_pec2 = pecuzal_embedding_update(data2;
#                                                             τs = taus , w = w2)
# dim_pec_ref2 = size(Y_pec2,2)
# Tw_pec_ref2 = sum(τ_vals_pec2)
#
# ## Make all reconstructions equally long
# N1 = length(Y_tde1)
# N2 = length(Y_GA1)
# N3 = length(Y_mdop1)
# N4 = length(Y_pec1)
# NN1 = minimum(hcat(N1, N2, N3, N4))
#
# N1 = length(Y_tde2)
# N2 = length(Y_GA2)
# N3 = length(Y_mdop2)
# N4 = length(Y_pec2)
# NN2 = minimum(hcat(N1, N2, N3, N4))
#
# ## RQA
# R11 = RecurrenceMatrix(Y_tde1[1:NN1,:], ε; fixedrate = true)
# R21 = RecurrenceMatrix(Y_GA1[1:NN1,:], ε; fixedrate = true)
# R31 = RecurrenceMatrix(Y_mdop1[1:NN1,:], ε; fixedrate = true)
# R41 = RecurrenceMatrix(Y_pec1[1:NN1,:], ε; fixedrate = true)
#
# R12 = RecurrenceMatrix(Y_tde2[1:NN2,:], ε; fixedrate = true)
# R22 = RecurrenceMatrix(Y_GA2[1:NN2,:], ε; fixedrate = true)
# R32 = RecurrenceMatrix(Y_mdop2[1:NN2,:], ε; fixedrate = true)
# R42 = RecurrenceMatrix(Y_pec2[1:NN2,:], ε; fixedrate = true)
#
# RQA11 = rqa(R11; theiler = w1, lmin = lmin)
# RQA21 = rqa(R21; theiler = w1, lmin = lmin)
# RQA31 = rqa(R31; theiler = w1, lmin = lmin)
# RQA41 = rqa(R41; theiler = w1, lmin = lmin)
#
# RQA12 = rqa(R12; theiler = w2, lmin = lmin)
# RQA22 = rqa(R22; theiler = w2, lmin = lmin)
# RQA32 = rqa(R32; theiler = w2, lmin = lmin)
# RQA42 = rqa(R42; theiler = w2, lmin = lmin)
#
# RQA_tde_ref1 = zeros(4)
# RQA_tde_ref2 = zeros(4)
# RQA_GA_ref1 = zeros(4)
# RQA_GA_ref2 = zeros(4)
# RQA_mdop_ref1 = zeros(4)
# RQA_mdop_ref2 = zeros(4)
# RQA_pec_ref1 = zeros(4)
# RQA_pec_ref2 = zeros(4)
#
# RQA_tde_ref1[1] = RQA11.ENTR
# RQA_tde_ref1[2] = RQA11.LAM
# RQA_tde_ref1[3] = RQA11.RTE
# RQA_tde_ref1[4] = RQA11.TRANS
# RQA_tde_ref2[1] = RQA12.ENTR
# RQA_tde_ref2[2] = RQA12.LAM
# RQA_tde_ref2[3] = RQA12.RTE
# RQA_tde_ref2[4] = RQA12.TRANS
#
# RQA_GA_ref1[1] = RQA21.ENTR
# RQA_GA_ref1[2] = RQA21.LAM
# RQA_GA_ref1[3] = RQA21.RTE
# RQA_GA_ref1[4] = RQA21.TRANS
# RQA_GA_ref2[1] = RQA22.ENTR
# RQA_GA_ref2[2] = RQA22.LAM
# RQA_GA_ref2[3] = RQA22.RTE
# RQA_GA_ref2[4] = RQA22.TRANS
#
# RQA_mdop_ref1[1] = RQA31.ENTR
# RQA_mdop_ref1[2] = RQA31.LAM
# RQA_mdop_ref1[3] = RQA31.RTE
# RQA_mdop_ref1[4] = RQA31.TRANS
# RQA_mdop_ref2[1] = RQA32.ENTR
# RQA_mdop_ref2[2] = RQA32.LAM
# RQA_mdop_ref2[3] = RQA32.RTE
# RQA_mdop_ref2[4] = RQA32.TRANS
#
# RQA_pec_ref1[1] = RQA41.ENTR
# RQA_pec_ref1[2] = RQA41.LAM
# RQA_pec_ref1[3] = RQA41.RTE
# RQA_pec_ref1[4] = RQA41.TRANS
# RQA_pec_ref2[1] = RQA42.ENTR
# RQA_pec_ref2[2] = RQA42.LAM
# RQA_pec_ref2[3] = RQA42.RTE
# RQA_pec_ref2[4] = RQA42.TRANS

println("Start sampling the time series")
# sampling the time series

#preallocation
RQA_tde1 = zeros(4,samplesize)
RQA_tde2 = zeros(4,samplesize)
RQA_GA1 = zeros(4,samplesize)
RQA_GA2 = zeros(4,samplesize)
RQA_mdop1 = zeros(4,samplesize)
RQA_mdop2 = zeros(4,samplesize)
RQA_pec1 = zeros(4,samplesize)
RQA_pec2 = zeros(4,samplesize)

dim_tde1 = zeros(Int, samplesize)
dim_tde2 = zeros(Int, samplesize)
Tw_tde1 = zeros(samplesize)
Tw_tde2 = zeros(samplesize)

dim_GA1 = zeros(Int, samplesize)
dim_GA2 = zeros(Int, samplesize)
Tw_GA1 = zeros(samplesize)
Tw_GA2 = zeros(samplesize)

dim_mdop1 = zeros(Int, samplesize)
dim_mdop2 = zeros(Int, samplesize)
Tw_mdop1 = zeros(samplesize)
Tw_mdop2 = zeros(samplesize)

dim_pec1 = zeros(Int, samplesize)
dim_pec2 = zeros(Int, samplesize)
Tw_pec1 = zeros(samplesize)
Tw_pec2 = zeros(samplesize)

#for i = 1:samplesize
i = 1
display("run: $i")
# draw the actual sample
idx1 = idxs1[i]
idx2 = idxs2[i]
x1 = data1[idx1:idx1+N]
x2 = data2[idx2:idx2+N]

w1 = estimate_delay(x1, "mi_min")
w2 = estimate_delay(x2, "mi_min")

## Perform embeddings
#Standard TDE
Y_tde1, τ_tde1, _ = optimal_traditional_de(x1; τs = taus)
Y_tde2, τ_tde2, _ = optimal_traditional_de(x2; τs = taus)

dim_tde1[i] = size(Y_tde1,2)
dim_tde2[i] = size(Y_tde2,2)
τ_tdes1 = [(i-1)*τ_tde1 for i = 1:(size(Y_tde1,2))]
τ_tdes2 = [(i-1)*τ_tde2 for i = 1:(size(Y_tde2,2))]
Tw_tde1[i] = sum(τ_tdes1)
Tw_tde2[i] = sum(τ_tdes2)

#MDOP embedding
tw1 = mdop_maximum_delay(x1)
tw2 = mdop_maximum_delay(x2)
lm1 = DelayEmbeddings.findlocalminima(tw1[2])
if lm1[1] ≤ 3
    lmm1 = try lm1[2]
    catch
        lm1[1]
    end
else
    lmm1 = lm1[1]
end
lm2 = DelayEmbeddings.findlocalminima(tw2[2])
if lm2[1] ≤ 3
    lmm2 = try lm2[2]
    catch
        lm2[1]
    end
else
    lmm2 = lm2[1]
end
tau_max1 = lmm1
tau_max2 = lmm2
taus_mdop1 = 0:tau_max1
taus_mdop2 = 0:tau_max2

Y_mdop1, τ_vals_mdop1, ts_vals_mdop1, FNNs_mdop1 , βs_mdop1 = mdop_embedding(x1;
                                                    τs = taus_mdop1 , w = w1)
dim_mdop1[i] = size(Y_mdop1,2)
Tw_mdop1[i] = sum(τ_vals_mdop1)

Y_mdop2, τ_vals_mdop2, ts_vals_mdop2, FNNs_mdop2 , βs_mdop2 = mdop_embedding(x2;
                                                    τs = taus_mdop2 , w = w2)
dim_mdop2[i] = size(Y_mdop2,2)
Tw_mdop2[i] = sum(τ_vals_mdop2)

#Garcia & Almeida
Y_GA1, τ_vals_GA1, ts_vals_GA1, FNNs_GA1 , ns_GA1 = garcia_almeida_embedding(x1;
                                                    τs = taus , w = w1, T = w1)
dim_GA1[i] = size(Y_GA1,2)
Tw_GA1[i] = sum(τ_vals_GA1)

Y_GA2, τ_vals_GA2, ts_vals_GA2, FNNs_GA2 , ns_GA2 = garcia_almeida_embedding(x2;
                                                    τs = taus , w = w2, T = w2)
dim_GA2[i] = size(Y_GA2,2)
Tw_GA2[i] = sum(τ_vals_GA2)

#Pecuzal
Y_pec1, τ_vals_pec1, ts_vals_pec1, Ls_pec1 , εs_pec1 = pecuzal_embedding_update(x1;
                                                            τs = taus , w = w1)
dim_pec1[i] = size(Y_pec1,2)
Tw_pec1[i] = sum(τ_vals_pec1)

Y_pec2, τ_vals_pec2, ts_vals_pec2, Ls_pec2 , εs_pec2 = pecuzal_embedding_update(x2;
                                                            τs = taus , w = w2)
dim_pec2[i] = size(Y_pec2,2)
Tw_pec2[i] = sum(τ_vals_pec2)

## Make all reconstructions equally long
N1 = length(Y_tde1)
N2 = length(Y_GA1)
N3 = length(Y_mdop1)
N4 = length(Y_pec1)
NN1 = minimum(hcat(N1, N2, N3, N4))

N1 = length(Y_tde2)
N2 = length(Y_GA2)
N3 = length(Y_mdop2)
N4 = length(Y_pec2)
NN2 = minimum(hcat(N1, N2, N3, N4))

## RQA
R11 = RecurrenceMatrix(Y_tde1[1:NN1,:], ε; fixedrate = true)
R21 = RecurrenceMatrix(Y_GA1[1:NN1,:], ε; fixedrate = true)
R31 = RecurrenceMatrix(Y_mdop1[1:NN1,:], ε; fixedrate = true)
R41 = RecurrenceMatrix(Y_pec1[1:NN1,:], ε; fixedrate = true)

R12 = RecurrenceMatrix(Y_tde2[1:NN2,:], ε; fixedrate = true)
R22 = RecurrenceMatrix(Y_GA2[1:NN2,:], ε; fixedrate = true)
R32 = RecurrenceMatrix(Y_mdop2[1:NN2,:], ε; fixedrate = true)
R42 = RecurrenceMatrix(Y_pec2[1:NN2,:], ε; fixedrate = true)

RQA11 = rqa(R11; theiler = w1, lmin = lmin)
RQA21 = rqa(R21; theiler = w1, lmin = lmin)
RQA31 = rqa(R31; theiler = w1, lmin = lmin)
RQA41 = rqa(R41; theiler = w1, lmin = lmin)

RQA12 = rqa(R12; theiler = w2, lmin = lmin)
RQA22 = rqa(R22; theiler = w2, lmin = lmin)
RQA32 = rqa(R32; theiler = w2, lmin = lmin)
RQA42 = rqa(R42; theiler = w2, lmin = lmin)

RQA_tde1[1,i] = RQA11.ENTR
RQA_tde1[2,i] = RQA11.LAM
RQA_tde1[3,i] = RQA11.RTE
RQA_tde1[4,i] = RQA11.TRANS
RQA_tde2[1,i] = RQA12.ENTR
RQA_tde2[2,i] = RQA12.LAM
RQA_tde2[3,i] = RQA12.RTE
RQA_tde2[4,i] = RQA12.TRANS

RQA_GA1[1,i] = RQA21.ENTR
RQA_GA1[2,i] = RQA21.LAM
RQA_GA1[3,i] = RQA21.RTE
RQA_GA1[4,i] = RQA21.TRANS
RQA_GA2[1,i] = RQA22.ENTR
RQA_GA2[2,i] = RQA22.LAM
RQA_GA2[3,i] = RQA22.RTE
RQA_GA2[4,i] = RQA22.TRANS

RQA_mdop1[1,i] = RQA31.ENTR
RQA_mdop1[2,i] = RQA31.LAM
RQA_mdop1[3,i] = RQA31.RTE
RQA_mdop1[4,i] = RQA31.TRANS
RQA_mdop2[1,i] = RQA32.ENTR
RQA_mdop2[2,i] = RQA32.LAM
RQA_mdop2[3,i] = RQA32.RTE
RQA_mdop2[4,i] = RQA32.TRANS

RQA_pec1[1,i] = RQA41.ENTR
RQA_pec1[2,i] = RQA41.LAM
RQA_pec1[3,i] = RQA41.RTE
RQA_pec1[4,i] = RQA41.TRANS
RQA_pec2[1,i] = RQA42.ENTR
RQA_pec2[2,i] = RQA42.LAM
RQA_pec2[3,i] = RQA42.RTE
RQA_pec2[4,i] = RQA42.TRANS
#end

# save results
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_tde_ref1.csv", RQA_tde_ref1)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_tde_ref2.csv", RQA_tde_ref2)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_GA_ref1.csv", RQA_GA_ref1)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_GA_ref2.csv", RQA_GA_ref2)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_mdop_ref1.csv", RQA_mdop_ref1)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_mdop_ref2.csv", RQA_mdop_ref2)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_pec_ref1.csv", RQA_pec_ref1)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_pec_ref2.csv", RQA_pec_ref2)
#
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_tde_ref1.csv", dim_tde_ref1)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_tde_ref2.csv", dim_tde_ref2)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_tde_ref1.csv", Tw_tde_ref1)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_tde_ref2.csv", Tw_tde_ref2)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_GA_ref1.csv", dim_GA_ref1)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_GA_ref2.csv", dim_GA_ref2)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_GA_ref1.csv", Tw_GA_ref1)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_GA_ref2.csv", Tw_GA_ref2)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_mdop_ref1.csv", dim_mdop_ref1)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_mdop_ref2.csv", dim_mdop_ref2)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_mdop_ref1.csv", Tw_mdop_ref1)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_mdop_ref2.csv", Tw_mdop_ref2)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_pec_ref1.csv", dim_pec_ref1)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_pec_ref2.csv", dim_pec_ref2)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_pec_ref1.csv", Tw_pec_ref1)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_pec_ref2.csv", Tw_pec_ref2)
#
#
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_tde1.csv", RQA_tde1)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_tde2.csv", RQA_tde2)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_GA1.csv", RQA_GA1)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_GA2.csv", RQA_GA2)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_mdop1.csv", RQA_mdop1)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_mdop2.csv", RQA_mdop2)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_pec1.csv", RQA_pec1)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/RQA_pec2.csv", RQA_pec2)
#
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_tde1.csv", dim_tde1)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_tde2.csv", dim_tde2)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_tde1.csv", Tw_tde1)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_tde2.csv", Tw_tde2)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_GA1.csv", dim_GA1)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_GA2.csv", dim_GA2)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_GA1.csv", Tw_GA1)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_GA2.csv", Tw_GA2)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_mdop1.csv", dim_mdop1)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_mdop2.csv", dim_mdop2)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_mdop1.csv", Tw_mdop1)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_mdop2.csv", Tw_mdop2)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_pec1.csv", dim_pec1)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/dim_pec2.csv", dim_pec2)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_pec1.csv", Tw_pec1)
# writedlm("./scripts/Experimental data/Standard electrochemical chaos/RQA results/Tw_pec2.csv", Tw_pec2)

max_stars_repo_path: scripts/Experimental data/Standard electrochemical chaos/comm_analyze_standard_ecoscillator_2nd_try.jl


text:
# include("sei3r.jl")
# include("evolution.jl")

using DiseaseOutbreak
# TODO: make numerical tests to replace plot tests
using Test

# default parameters for RCAN model
d0 = getParams(5.0,6.0,0.5,0.1,0.1,0.15,0.06,6.0,8.0,0.4,SEI3R())
# social distancing reducing transmission probs by 1/2
d1 = getParams(5.0,6.0,0.5/2,0.1/2,0.1/2,0.15,0.06,6.0,8.0,0.4,SEI3R())

N = 1000.0
state0 = initialize(1.0/N,d0)
(states0,deltaStates0) = evolve(N,state0,d0,300)
plotEvolution(N,states0,d0)
# savefig("dist300.pdf")
state0 = initialize(1.0/N,d0)
(states0,deltaStates0) = evolve(N,state0,d0,60)
(states1,deltaStates1) = evolve(N,states0[:,60],d1,241)
plotEvolution(N,[states0 states1[:,2:end]],d0)
# savefig("dist60.pdf")
state0 = initialize(1.0/N,d0)
(states0,deltaStates0) = evolve(N,state0,d0,30)
(states1,deltaStates1) = evolve(N,states0[:,30],d1,271)
plotEvolution(N,[states0 states1[:,2:end]],d0)
# savefig("dist30.pdf")
state1 = initialize(1.0/N,d1)
(states1,deltaStates1) = evolve(N,state1,d1,300)
plotEvolution(N,states1,d1)
# savefig("dist0.pdf")

ext: jl


text:
module TestKernels
using GaussianProcesses, Calculus
using Test, LinearAlgebra, Statistics, Random
using ForwardDiff
using GaussianProcesses: EmptyData, update_target_and_dtarget!,
      cov_ij, dKij_dθp, dKij_dθ!,
      get_params, set_params!, num_params, StationaryARD, WeightedEuclidean
import Calculus: gradient

Random.seed!(1)
const d, n, n2 = 3, 6, 3
function testkernel(kern::Kernel)
    X = randn(d, n)
    X2 = randn(d, n2)
    y = randn(n)

    # Random columns
    i, j = rand(1:n), rand(1:n2)
    Xi = view(X, :, i)
    Xj = view(X, :, j) # works since n2 < n
    X2j = view(X2, :, j)

    # Preallocate some matrices
    cK = zeros(n, n)
    cK2 = zeros(n, n2)
    @test length(GaussianProcesses.get_param_names(kern)) ==
        length(get_params(kern)) ==
        num_params(kern)

    @testset "get and set params" begin
        kcopy = deepcopy(kern)
        params_1 = randn(num_params(kcopy))
        set_params!(kcopy, params_1)
        params_2 = get_params(kcopy)
        @test params_1 ≈ params_2
    end
    @testset "Variance" begin
        spec = cov(kern, X)
        @test spec ≈ invoke(cov, Tuple{Kernel, Matrix{Float64}}, kern, X)
        @test spec[i, j] ≈ cov(kern, Xi, Xj)

        key = GaussianProcesses.kernel_data_key(kern, X, X)
        @test typeof(key) == String
        # check we've overwritten the default if necessary
        kdata = GaussianProcesses.KernelData(kern, X, X)
        if typeof(kdata) != EmptyData
            @test key != "EmptyData"
        end
    end

    data = GaussianProcesses.KernelData(kern, X, X)
    data12 = GaussianProcesses.KernelData(kern, X, X2)

    @testset "Covariance" begin
        spec = cov(kern, X, X2)
        @test spec[i,j] ≈ cov(kern, Xi, X2j)
        spec = cov(kern, X, X2, data12)
        @test spec[i,j] ≈ cov(kern, Xi, X2j)
    end

    @testset "Gradient" begin
        nparams = num_params(kern)
        init_params = Vector(get_params(kern))
        dK = zeros(nparams)
        i, j = 3, 5
        dKij_dθ!(dK, kern, X, X, data, i, j, d, nparams)
        dK1 = copy(dK)
        dKij_dθ!(dK, kern, X, X, EmptyData(), i, j, d, nparams)
        dK2 = copy(dK)
        @test dK1 ≈ dK2
        for p in 1:nparams
            @test dK[p] ≈ dKij_dθp(kern, X, X, data,        i, j, p, d)
            @test dK[p] ≈ dKij_dθp(kern, X, X, EmptyData(), i, j, p, d)
            try
                dkp = dKij_dθp(kern, X, X, i, j, p, d)
                @test dkp == dK[p]
            catch
                # that's OK too
                continue
            end
        end
        if nparams > 0
            numer_grad = Calculus.gradient(init_params) do params
                set_params!(kern, params)
                t = cov_ij(kern, X, X, data, i, j, d)
                set_params!(kern, init_params)
                t
            end
            theor_grad = dK
            @test numer_grad ≈ theor_grad rtol=1e-3 atol=1e-3
        end
    end
    @testset "Gradient stack X1 ≠ X2" begin
        nparams = num_params(kern)
        init_params = Vector(get_params(kern))
        stack1 = Array{Float64}(undef, n, n2, nparams)
        stack2 = Array{Float64}(undef, n, n2, nparams)

        GaussianProcesses.grad_stack!(stack1, kern, X, X2, data12)

        theor_grad = vec(sum(stack1; dims=[1,2]))
        if nparams > 0
            numer_grad = Calculus.gradient(init_params) do params
                set_params!(kern, params)
                t = sum(cov(kern, X, X2, data12))
                set_params!(kern, init_params)
                t
            end
            @test theor_grad ≈ numer_grad rtol=1e-2 atol=1e-2
        end

        GaussianProcesses.grad_stack!(stack2, kern, X, X2, EmptyData())
        # invoke(GaussianProcesses.grad_stack!,
               # Tuple{AbstractArray, Kernel, Matrix{Float64}, Matrix{Float64},
                     # EmptyData},
               # stack2, kern, X, X2, EmptyData())
        @test stack1 ≈ stack2 rtol=1e-3 atol=1e-3
    end

    @testset "Gradient stack" begin
        nparams = num_params(kern)
        init_params = Vector(get_params(kern))
        stack1 = Array{Float64}(undef, n, n, nparams)
        stack2 = Array{Float64}(undef, n, n, nparams)

        GaussianProcesses.grad_stack!(stack1, kern, X, X, data)
        invoke(GaussianProcesses.grad_stack!,
               Tuple{AbstractArray, Kernel, Matrix{Float64}, Matrix{Float64},
                     EmptyData},
               stack2, kern, X, X, EmptyData())
        @test stack1 ≈ stack2

        theor_grad = vec(sum(stack1; dims=[1,2]))
        if nparams > 0
            numer_grad = Calculus.gradient(init_params) do params
                set_params!(kern, params)
                t = sum(cov(kern, X))
                set_params!(kern, init_params)
                t
            end
            @test theor_grad ≈ numer_grad rtol=1e-1 atol=1e-2
        end
    end

    @testset "dtarget" begin
        nparams = num_params(kern)
        gp = GPE(X, y, MeanConst(0.0), kern, -3.0)
        init_params = Vector(get_params(gp))
        update_target_and_dtarget!(gp)
        theor_grad = copy(gp.dtarget)
        if nparams > 0
            numer_grad = Calculus.gradient(init_params) do params
                set_params!(gp, params)
                update_target!(gp)
                t = gp.target
                set_params!(gp, init_params)
                t
            end
            @test theor_grad ≈ numer_grad rtol=1e-3 atol=1e-3
        end
    end

    @testset "EmptyData" begin
        gp = GPE(X, y, MeanConst(0.0), kern, -3.0)
        gp_empty = GPE(X, y, MeanConst(0.0), kern, -3.0, EmptyData())
        update_target_and_dtarget!(gp_empty)
        update_target_and_dtarget!(gp)
        @test gp.dmll ≈ gp_empty.dmll rtol=1e-6 atol=1e-6
    end

    @testset "predict gradient" begin
        gp = GPE(X, y, MeanConst(0.0), kern, -3.0)
        f = x -> sum(predict_y(gp, reshape(x, :, 1)))[1]
        z = rand(d)
        autodiff_grad = ForwardDiff.gradient(f, z)
        numer_grad = Calculus.gradient(f, z)
        @test autodiff_grad ≈ numer_grad rtol = 1e-3 atol = 1e-3
    end
end

@testset "kernel shortcuts" begin
    x1 = [0.1, 0.2]
    x2 = [1.1, 1.0]

    for pairs in [
                  (SEIso(1.0, 1.2), SE(1.0, 1.2)),
                  (SEArd([1.0, 1.5], 1.3), SE([1.0, 1.5], 1.3)),
                  (RQIso(1.0, 1.2, 0.5), RQ(1.0, 1.2, 0.5)),
                  (RQArd([1.0, 1.5], 1.3, 0.5), RQ([1.0, 1.5], 1.3, 0.5)),
                  (Matern(1/2, 1.0, 1.2), Mat12Iso(1.0, 1.2)),
                  (Matern(3/2, 1.0, 1.2), Mat32Iso(1.0, 1.2)),
                  (Matern(5/2, 1.0, 1.2), Mat52Iso(1.0, 1.2)),
                  (Matern(1/2, [1.0, 1.5], 1.2), Mat12Ard([1.0, 1.5], 1.2)),
                  (Matern(3/2, [1.0, 1.5], 1.2), Mat32Ard([1.0, 1.5], 1.2)),
                  (Matern(5/2, [1.0, 1.5], 1.2), Mat52Ard([1.0, 1.5], 1.2)),
                  (Lin(1.0), LinIso(1.0)),
                  (Lin([1.0, 1.5]), LinArd([1.0, 1.5])),
                 ]
        @test cov(pairs[1], x1, x2) == cov(pairs[2], x1, x2)
    end

end

@testset "Kernels" begin
    ll = rand(d)
    kernels = [# Isotropic kernels
               SEIso(1.0, 1.0), Mat12Iso(1.0,1.0), Mat32Iso(1.0,1.0), Mat52Iso(1.0,1.0),
               RQIso(1.0, 1.0, 1.0), Periodic(1.0, 1.0, 2π),
               # Non-isotropic
               Lin(1.0), Poly(0.0, 0.0, 2), Noise(1.0),
               # Constant kernel
               Const(1.0),
               # ARD kernels
               SEArd(ll, 1.0), Mat12Ard(ll, 1.0), Mat32Ard(ll, 1.0), Mat52Ard(ll, 1.0),
               RQArd(ll, 0.0, 2.0), LinArd(ll),
               # Composite kernels
               SEIso(1.0, 1.0) + Mat12Iso(1.0, 1.0),
               (SEIso(1.0, 1.0) + Mat12Iso(1.0, 1.0)) + Lin(1.0),
               SEIso(1.0, 1.0) * Mat12Iso(1.0, 1.0),
               (SEIso(1.0, 1.0) * Mat12Iso(1.0, 1.0)) * Lin(1.0),
               # Fixed Kernel
               fix(RQIso(1.0, 1.0, 1.0), :lσ),
               fix(RQIso(1.0, 1.0, 1.0)),
               # Sum and Product and Fix
               SEIso(1.0, 1.0) * Mat12Iso(1.0, 1.0) +
               Lin(1.0) * fix(RQIso(1.0, 1.0, 1.0), :lσ)]
    @testset for kern in kernels
        println("\tTesting ", nameof(typeof(kern)), "...")
        testkernel(kern)
    end
    @testset "Masked" for kernel in kernels
        println("\tTesting masked", nameof(typeof(kernel)), "...")
        if isa(kernel, LinArd) || isa(kernel, GaussianProcesses.StationaryARD)
            par = get_params(kernel)
            k_masked = typeof(kernel).name.wrapper([par[1]], par[d+1:end]...)
            kern = Masked(k_masked, [1])
        else
            kern = Masked(kernel, [1])
        end
        println("\tTesting masked ", nameof(typeof(kern)), "...")
        testkernel(kern)
    end
    @testset "autodiff" for kernel in kernels[1:end-1]
        println("\tTesting autodiff ", nameof(typeof(kernel)), "...")
        testkernel(autodiff(kernel))
    end

end # testset
end # module


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM, Test, LinearAlgebra

x1 = [2.0, 3.0, 4.0]
x2 = [6.0, 3.0, 2.0]
x3 = [2.0, 5.0, 1.0]
x4 = [4.0, 3.0, 6.0]
x5 = 0.5*(x1+x2)
x6 = 0.5*(x2+x3)
x7 = 0.5*(x3+x1)
x8 = 0.5*(x1+x4)
x9 = 0.5*(x2+x4)
x10 = 0.5*(x3+x4)
X = Dict(
    1 => x1, 2 => x2, 3 => x3, 4 => x4, 5 => x5,
    6 => x6, 7 => x7, 8 => x8, 9 => x9, 10 => x10)
u = Dict(i => zeros(3) for i in 1:10)

element = Element(Tet10, (1, 2, 3, 4, 5, 6, 7, 8, 9, 10))
update!(element, "youngs modulus", 480.0)
update!(element, "poissons ratio", 1/3)
update!(element, "geometry", X)
update!(element, "displacement", u)
problem = Problem(Elasticity, "tet10", 3)
add_element!(problem, element)
time = 0.0
assemble!(problem, time)
eigs = real(eigvals(Matrix(problem.assembly.K)))
eigs_expected = [8809.45, 4936.01, 2880.56, 2491.66, 2004.85,
                 1632.49, 1264.32, 1212.42, 817.905,
                 745.755, 651.034, 517.441, 255.1, 210.955,
                 195.832, 104.008, 72.7562, 64.4376, 53.8515,
                 23.8417, 16.6354, 9.54682, 6.93361, 2.22099,
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
@test isapprox(eigs, eigs_expected; atol=1.0e-2)


text:
###################################
# Destination location selectors
###################################

"""
Destination location selector based on DP (Demographic Profile)

Selects destination DA_work for an agent by randomly choosing the company he works in

**Arguments**
* `agent_profile` : agent demographic profile
* `business_data` : array with dictionaries  with business location, industry and estimated number of employees
* `industry` : dictionary matching industry demographic data from business_data with the ones selected

**Assumptions based on agent demographic profile
- agents work in the business in accordance with their work_industry
- exact number of employess in businesses is estimated in each iteration based on "Number of employees" intervals
- differences between maximum number of employees and actual number of employees are used as a probability weights for each workplace location
"""
function destination_location!(agent_profile::DataFrames.DataFrame,
                            business_data::Array{Dict{Symbol,Union{String, Int,UnitRange{Int}}},1},
                            industry::Dict{String,Array{String,1}})
    indices  = findall((in)(Set(industry[agent_profile.work_industry[1]])),[ds[:ICLS_DESC] for ds in business_data])
    weights = StatsBase.fweights([maximum(business_data[i][:IEMP_DESC]) - business_data[i][:no_of_workers] for i in indices])
    index = StatsBase.sample(indices,weights)
    DA_work = business_data[index][:DA_ID]
    business_data[index][:no_of_workers] += 1
    agent_profile.DA_work[1] = DA_work
    return nothing
end

"""
Destination location selector based on JM (Jurney Matrix)

Selects destination DA_work for an agent randomly weighted by Pij Journey Matrix

**Arguments**
* `agent_profile` : agent demographic profile
* `flow_dictionary ` :  Dictionary mapping Pij Journey Matrix columns to DAs
* `flow_matrix ` :  Pij Journey Matrix with *FlowVolume* from *DA_I* to *DA_J*
"""
function destination_location!(agent_profile::DataFrames.DataFrame,
                                flow_dictionary::Dict{Int,Int},
                                flow_matrix::SparseMatrixCSC{Int,Int})
    row = flow_dictionary[agent_profile.DA_home[1]]
    column = StatsBase.sample(StatsBase.fweights(flow_matrix[row,:]))
    agent_profile.DA_work[1] = collect(keys(flow_dictionary))[something(findfirst(isequal(column), collect(values(flow_dictionary))),rand(1:length(flow_dictionary)))]
    return nothing
end

###################################
## Selection based on Both Journey Matrix and DP (Demographic Profile)

"""
Destination location selector based on both, JM (Jurney Matrix) and DP (Demographic Profile)

Selects destination DA_work for an agent choosing a proper industry and then by randomly weighted DAs with Pij Journey Matrix

**Arguments**
* `agent_profile` : agent demographic profile
* `sim_data ` :  `SimData` object
* `industry` : dictionary matching industry demographic data from business_data with the ones selected

**Assumptions based on agent demographic profile
- agents work in the business in accordance with their work_industry
- list of possible businesses locations is returned based on agent's work_industry
- then probability of selecting one of possible businesses is weighted by flows from DA_home to DA corresponding to each business
"""
function destination_location!(agent_profile::DataFrames.DataFrame,
                            sim_data::OpenStreetMapXSim.SimData;
                            industry::Dict{String,Array{String,1}} = OpenStreetMapXSim.industry)
    row = sim_data.DAs_flow_dictionary[agent_profile.DA_home[1]]
    indices  = findall((in)(Set(industry[agent_profile.work_industry[1]])), [ds[:ICLS_DESC] for ds in sim_data.business_data])
    columns = [sim_data.DAs_flow_dictionary[sim_data.business_data[index][:DA_ID]] for index in indices if haskey(sim_data.DAs_flow_dictionary,sim_data.business_data[index][:DA_ID])]
    weights = StatsBase.fweights([sim_data.DAs_flow_matrix[row,column] for column in columns])
    index = StatsBase.sample(indices,weights)
    DA_work = sim_data.business_data[index][:DA_ID]
    sim_data.business_data[index][:no_of_workers] += 1
    agent_profile.DA_work[1] = DA_work
    return nothing
end


text:
function xml_to_dict(node::XMLElement; path="")
    dict = Dict{Symbol, Any}()
    for a=attributes(node)
        bt = base_type_parser("$(name(node)).$(name(a))")
        val = value(a)
        if !(bt isa Nothing)
            val = bt(val)
        end
        dict[Symbol(name(a))] = val
    end
    for n=child_elements(node)
        nxt = xml_to_dict(n)
        if name(n) in keys(dict)
            cur = dict[Symbol(name(n))]
            if cur isa Array
                push!(cur, nxt)
            else
                dict[Symbol(name(n))] = [cur, nxt]
            end
        else
            dict[Symbol(name(n))] = nxt
        end
    end
    return dict
end

xsd_types = Dict("xsd:byte"=>(x -> parse(Int8, x)),
                 "xsd:decimal"=>(x -> parse(Float64, x)),
                 "xsd:int"=>(x -> parse(Int32, x)),
                 "xsd:integer"=>(x -> parse(Int64, x)),
                 "xsd:long"=>(x -> parse(Int64, x)),
                 "xsd:short"=>(x -> parse(Int16, x)),
                 "xsd:unsignedLong"=>(x -> parse(UInt64, x)),
                 "xsd:unsignedInt"=>(x -> parse(UInt32, x)),
                 "xsd:unsignedShort"=>(x -> parse(UInt16, x)),
                 "xsd:unsignedByte"=>(x -> parse(UInt8, x)),
                 "xsd:dateTime"=>(x->DateTime(x, Date.DateFormat("yyyy-mm-ddTHH:MM:SS.sss"))),
                 "xsd:float"=>(x -> parse(Float64, x)),
                 "xsd:double"=>(x -> parse(Float64, x)),
                 "xsd:boolean"=>(x -> parse(Bool, x))
                )

function parse_schema(doc::Nothing)
end

function parse_schema(doc::XMLDocument)
    children = child_elements(root(doc))
    return Dict(attribute(e, "name")=>parse_schema(e) for e=children)
end

function parse_schema(element::XMLElement)
    if name(element) == "simpleType"
        r = find_element(element, "restriction")
        return attribute(r, "base")
    elseif name(element) == "complexType"
        return Dict(attribute(e, "name")=>parse_schema(e) for e=element["attribute"])
    elseif name(element) == "attribute"
        if has_attribute(element, "type")
            return attribute(element, "type")
        else
            return parse_schema(find_element(element, "simpleType"))
        end
    elseif name(element) == "element"
        try
            return parse_schema([element["simpleType"];element["complexType"]][1])
        catch
        end
    end
end

function base_type_parser(attribute)
    parts = split(attribute, ".")
    element = schema
    try
        for p=parts
            element = element[p]
        end
        while element in keys(schema)
            element = schema[element]
        end
        return xsd_types[element]
    catch KeyError
    end
end

schema_path = joinpath(@__DIR__, "..", "deps", "ome.xsd")
schema = parse_schema(LightXML.parse_file(schema_path))


text:
preorder(node::T) where {T} = set_preorder!( T[], node )
postorder(node::T) where {T} = set_postorder!( T[], node )

function set_preorder!(order, node)
   push!(order, node)
   for child in node.children
      set_preorder!(order, child)
   end
   return order
end

function set_postorder!(order, node)
   for child in node.children
      set_postorder!(order, child)
   end
   push!(order, node)
   return order
end


text:
using MRIgeneralizedBloch
using Test
using LinearAlgebra
using StaticArrays
include("OCT_finite_difference_gradients.jl")

##
Npulse = 100
α = abs.(π/2 * sin.(π/2 * (1:Npulse) / Npulse))
TRF = 300e-6 .+ 200e-6 * cos.(π * (1:Npulse) / Npulse)
α[1] = π
TRF[1] = 500e-6
ω1 = α ./ TRF

R2slT = precompute_R2sl(ω1_max = 1.1 * maximum(ω1))
TR = 3.5e-3

B1 = 1
ω0 = 0
m0s = 0.25
R1f = 0.3
R1s = 2
R2f = 1 / 65e-3
T2s = 10e-6
Rx = 20

grad_list = [grad_m0s(), grad_R1f(), grad_R2f(), grad_Rx(), grad_R1s(), grad_T2s(), grad_ω0(), grad_B1()]
w = transpose([1/m0s;1/R1f;1/R2f;0;0;0;0;0;0].^2)

## ########################################################################
# Test dCRBdm
###########################################################################
# m = calculatesignal_linearapprox(α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT; grad_list=grad_list)
# CRB_fd, d_fd = dCRBdm_fd(m,w)

# m = calculatesignal_linearapprox(α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT; grad_list=grad_list, output=:realmagnetization)
# CRB, d = MRIgeneralizedBloch.dCRBdm(m, w)

# _dCRBdm    = [d(t,r,g)[i]    for t=1:size(m,1), r=1:size(m,2), g=1:size(m,3), i ∈ 1:11]
# _dCRBdm_fd = [d_fd(t,r,g)[i] for t=1:size(m,1), r=1:size(m,2), g=1:size(m,3), i ∈ 1:11]
# # _dCRBdm_fd[isnan.(_dCRBdm_fd)]
# @test _dCRBdm ≈ _dCRBdm_fd rtol = 1e-2

## ########################################################################
# Test OCT gradients
###########################################################################
# OCT_gradient: analytical
(F0, grad_ω1, grad_TRF) = MRIgeneralizedBloch.OCT_gradient(ω1, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT, grad_list, w)

## Finite Difference gradient: ω1
_grad_ω1_fd = grad_ω1_fd(w, ω1, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT, grad_list)
@test grad_ω1 ≈ _grad_ω1_fd rtol = 1e-1

## Finite Difference gradient: TRF
_grad_TRF_fd = grad_TRF_fd(w, ω1, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT, grad_list)
@test grad_TRF ≈ _grad_TRF_fd rtol = 1e-1

## OCT_TV_gradient: analytical
λ = 1
(F0, grad_ω1, grad_TRF) = MRIgeneralizedBloch.OCT_TV_gradient(ω1, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT, grad_list, w, λ)

## FD gradient: ω1
_,_grad_ω1_fd = grad_TV_ω1_fd(w, ω1, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT, grad_list, λ)
@test grad_ω1 ≈ _grad_ω1_fd rtol = 1e-1


text:
### File declaring all the various reaction networks the package is tested on ###
#A ;arge number of reaction networks are added to test as many different cases as possible.
#More test networks can be added favourably.

#Decalares 3 vectors to contain the various test sets (and a fourth with fixed concentrations coupled with the last test set).
reaction_networks_standard = Vector{DiffEqBase.AbstractReactionNetwork}(undef,10)
reaction_networks_hill = Vector{DiffEqBase.AbstractReactionNetwork}(undef,10)
reaction_networks_fixed_conc = Vector{DiffEqBase.AbstractReactionNetwork}(undef,10)
fixed_conc = Vector{Vector{Expr}}(undef,10)

#Standard reaction networks.
reaction_networks_standard[1] = @reaction_network begin
    (p1,p2,p3), ∅ → (X1,X2,X3)
    (k1,k2), X2 ⟷ X1 + 2X3
    (k3,k4), X1 ⟷ X3
    (d1,d2,d3), (X1,X2,X3) → ∅
end p1 p2 p3 k1 k2 k3 k4 d1 d2 d3

reaction_networks_standard[2] = @reaction_network begin
    mmR(X2,v1,K1), ∅ → X1
    mm(X1,v2,K2), ∅ → X2
    d, X1+X2 → ∅
end v1 K1 v2 K2 d

reaction_networks_standard[3] = @reaction_network begin
    mm(X2,v1,K1), ∅ → X1
    mm(X3,v2,K2), ∅ → X2
    (k1,k2), X1 ⟷ X4
    (k3,k4), X4 + X2 ⟷ X3 +X1
    d, (X1,X2,X3,X4) → ∅
end v1 K1 v2 K2 k1 k2 k3 k4 d

reaction_networks_standard[4] = @reaction_network begin
    mmR(X4,v1,K1), ∅ → X1
    mmR(X1,v2,K2), ∅ → X2
    mmR(X2,v3,K3), ∅ → X3
    mmR(X3,v4,K4), ∅ → X4
    (d1,d2,d3,d4), (X1,X2,X3,X4) → ∅
end v1 K1 v2 K2 v3 K3 v4 K4 d1 d2 d3 d4

reaction_networks_standard[5] = @reaction_network begin
    p, ∅ → X1
    (k1,k2), X1 ⟷ X2
    (k3,k4), X2 ⟷ X3
    (k5,k6), X3 ⟷ X4
    d, X4 → ∅
end p k1 k2 k3 k4 k5 k6 d

reaction_networks_standard[6] = @reaction_network begin
    (p1,p2), ∅ → (X1,X2)
    (k1,k2), 2X1 ⟷ X3
    (k3,k4), X2 + X3 ⟷ 4X4
    (k5,k6), X4 + X1 ⟷ 2X3
    d, (X1,X2,X3,X4) → ∅
end p1 p2 k1 k2 k3 k4 k5 k6 d

reaction_networks_standard[7] = @reaction_network begin
    (p1,p2,p3,p4,p5), ∅ → (X1,X2,X3,X4,X5)
    (k1,k2), X1 + 3X2 ⟷ 2X3 + 2X4
    (k3,k4), X2 + X3 ⟷ X5
    (mm(X1,v1,K1),k5), X5 ⟷ X1
    (mm(X2,v2,K2),k6), X1 ⟷ X2
    (d1,d2,d3,d4,d5), (X1,X2,X3,X4,X5) → ∅
end p1 p2 p3 p4 p5 k1 k2 k3 k4 k5 k6 v1 K1 v2 K2 d1 d2 d3 d4 d5

reaction_networks_standard[8] = @reaction_network begin
    (p1,p2), ∅ → (X1,X2)
    k1, X1 + X2 → X3
    (k2, k3), X3 ⟷ X4
    d, X4 → ∅
end p1 p2 k1 k2 k3 d

reaction_networks_standard[9] = @reaction_network begin
    (p1,p2,p3), ∅ ⟶ (X1,X2,X3)
    (d1,d2,d3), (X1,X2,X3) ⟶ ∅
    (k1,k2), X1 + X2 ⟷ X3
    (k3,k4), X3 ⟷ X4
    (k5,k6), X4 ⟷ X1 + X2
end p1 p2 p3 d1 d2 d3 k1 k2 k3 k4 k5 k6

reaction_networks_standard[10] = @reaction_network begin
    p, ∅ ⟶ X1
    (k1,k2), X1 ⟷ X3+ 2X2
    (k3,k4), X2 ⟷ X4
    (k5,k6), 3X3 + X4 ⟷ X5
    d, X5  ⟶ ∅
end p k1 k2 k3 k4 k5 k6 d

#Hill type reaction networks (these contaisn a parameter in the exponent and should have integer parameters).
reaction_networks_hill[1] = @reaction_network begin
    hillR(X2,v1,K1,n1),   ∅ → X1
    hillR(X1,v2,K2,n2),   ∅ → X2
    (d1,d2),          (X1,X2) → ∅
end v1 v2 K1 K2 n1 n2 d1 d2

reaction_networks_hill[2] = @reaction_network begin
    hillR(X3,v1,K1,n1),   ∅ → X1
    hillR(X1,v2,K2,n2),   ∅ → X2
    hillR(X2,v3,K3,n3),   ∅ → X3
    (d1,d2,d3),          (X1,X2,X3) → ∅
end v1 v2 v3 K1 K2 K3 n1 n2 n3 d1 d2 d3

reaction_networks_hill[3] = @reaction_network begin
    hillR(X2,v1,K1,n1), ∅ → X1
    hill(X1,v2,K2,n2), ∅ → X2
    d, X1+X2 → ∅
end v1 K1 n1 v2 K2 n2 d

reaction_networks_hill[4] = @reaction_network begin
    hillR(X2,v1,K1,n1)*hillR(X3,v1,K1,n1), ∅ → X1
    hillR(X1,v2,K2,n2)*hillR(X3,v2,K2,n2), ∅ → X2
    hillR(X1,v3,K3,n3)*hillR(X2,v3,K3,n3), ∅ → X3
    (d1,d2,d3), (X1,X2,X3)  ⟶ ∅
end v1 K1 n1 v2 K2 n2 v3 K3 n3 d1 d2 d3

reaction_networks_hill[5] = @reaction_network begin
    hillR(X2,v1,K1,n1)*hill(X4,v1,K1,n1), ∅ → X1
    hill(X5,v2,K2,n2), ∅ → X2
    hill(X3,v3,K3,n3), ∅ → X3
    hillR(X1,v4,K4,n4), ∅ → X4
    hill(X2,v5,K5,n5), ∅ → X5
    (k1,k2), X2 ⟷ X1 + 2X4
    (k3,k4), X4 ⟷ X3
    (k5,k6), 3X5 + X1 ⟷ X2
    (d1,d2,d3,d4,d5), (X1,X2,X3,X4,X5)  ⟶ ∅
end v1 K1 n1 v2 K2 n2 v3 K3 n3 v4 K4 n4 v5 K5 n5 k1 k2 k3 k4 k5 k6 d1 d2 d3 d4 d5

reaction_networks_hill[6] = @reaction_network begin
    v/10+hill(X1,v,K,n), ∅ → X1
    d, X1 → ∅
end v K n d

reaction_networks_hill[7] = @reaction_network begin
    v/10 + hill(X1,v,K,n), ∅ → X1+X2
    (k1,k2), X1 + X2 ↔ X3
    k3, X3 → X1
    d, (X1,X2,X3) → ∅
end v K n k1 k2 k3 d

reaction_networks_hill[8] = @reaction_network begin
    hill(X2,v1,K1,n1), ∅ → X1
    hillR(X1,v2,K2,n2)*hill(X3,v3,K3,n3), ∅ → X2
    hill(X2,v4,K4,n4), ∅ → X3
    (d1,d2,d3), (X1,X2,X3) → ∅
end v1 K1 n1 v2 K2 n2 v3 K3 n3 v4 K4 n4 d1 d2 d3

reaction_networks_hill[9] = @reaction_network begin
    hill(X1,v1,K1,n1)*hillR(X1,v2,K2,n2), ∅ → X1
    d, X1 → ∅
end v1 K1 n1 v2 K2 n2 d

reaction_networks_hill[10] = @reaction_network begin
    hill(X2,v1,K1,n1), ∅ → X1
    hillR(X4,v2,K2,n2), ∅ → X2
    (k1,k2), 2X1 + X2 ⟷ X3
    (k3,k4), 2X2 + X3 ⟷ X4
    (k5,k6), X1 + X2 + X3 + X4 ⟷ X5 + X6
    (d1,d2), (X5,X6) → ∅
end v1 K1 n1 v2 K2 n2 k1 k2 k3 k4 k5 k6 d1 d2

#Reaction networks with fixed concentrations.
#The fixed concentration(s) is declared underneth the reaction network.
reaction_networks_fixed_conc[1] = @reaction_network begin
    (k1,k2), X1 ↔ X2
    (k3,k4), X2 ↔ X3
    (k5,k6), X3 ↔ X1
end k1 k2 k3 k4 k5 k6 FC
fixed_conc[1] = [:(X1+X2+X3=FC)]

reaction_networks_fixed_conc[2] = @reaction_network begin
    (k1,k2), X1 ↔ 2X1
    (k3,k4), X1 + X2 ↔ X3
    (k5,k6), X3 ↔ X2
end k1 k2 k3 k4 k5 k6 FC
fixed_conc[2] = [:(X2+X3=FC)]

reaction_networks_fixed_conc[3] = @reaction_network begin
    (k1,k2*X5), X1 ↔ X2
    (k3*X5,k4), X3 ↔ X4
    (p+k5*X2*X3,d), ∅ ↔ X5
end k1 k2 k3 k4 p k5 d FC1 FC2
fixed_conc[3] = [:(X1+X2=FC1),:(X3+X4=FC2)]

reaction_networks_fixed_conc[4] = @reaction_network begin
    (k1,k2), X1 + X2 ↔ X3
    (mm(X3,v,K),d), ∅ ↔ X4
end k1 k2 v K d FC1 FC2
fixed_conc[4] = [:(X1+X3=FC1),:(X2+X3=FC2)]

reaction_networks_fixed_conc[5] = @reaction_network begin
    (k1,k2), X1 ↔ 2X2
    (k3,k4), 2X2 ↔ 3X3
    (k5,k6), 3X3 ↔ 4X4
end k1 k2 k3 k4 k5 k6 FC
fixed_conc[5] = [:(24X1+12X2+8X3+6X4=FC)]

reaction_networks_fixed_conc[6] = @reaction_network begin
    mmR(X1,v1,K1), X1 → X2
    mmR(X2,v2,K2), X2 → X3
    mmR(X3,v3,K3), X3 → X1
end v1 K1 v2 K2 v3 K3 FC
fixed_conc[6] = [:(X1+X2+X3=FC)]

reaction_networks_fixed_conc[7] = @reaction_network begin
    (k1,k2), X1 + X2 ↔ X3
    (mm(X3,v,K),d), ∅ ↔ X2
    (k3,k4), X2 ↔ X4
end k1 k2 k3 k4 v K d FC
fixed_conc[7] = [:(X1+X3=FC)]

reaction_networks_fixed_conc[8] = @reaction_network begin
    (k1,k2), X1 + X2 ↔ X3
    (mm(X3,v1,K1),mm(X4,v2,K2)), X3 ↔ X4
end k1 k2 v1 K1 v2 K2 FC1 FC2
fixed_conc[8] = [:(X1+X2+X4=FC1),:(X2+X3+X4=FC2)]

reaction_networks_fixed_conc[9] = @reaction_network begin
    (k1,k2), X1 + X2 ↔ X3
    (k3,k4), X3 + X4 ↔ X5
    (k5,k6), X5 + X6 ↔ X7
end k1 k2 k3 k4 k5 k6 FC1 FC2 FC3 FC4
fixed_conc[9] = [:(X1+X3+X5+X7=FC1),:(X2+X3+X5+X7=FC2),:(X4+X5+X7=FC3),:(X6+X7=FC4)]

reaction_networks_fixed_conc[10] = @reaction_network rnType begin
    kDeg,       (w,w2,w2v,v,w2v2,vP,σB,w2σB) ⟶ ∅
    kDeg,       vPp ⟶ phos
    (kBw,kDw),  2w ⟷ w2
    (kB1,kD1),  w2 + v ⟷ w2v
    (kB2,kD2),  w2v + v ⟷ w2v2
    kK1,        w2v ⟶ w2 + vP
    kK2,        w2v2 ⟶ w2v + vP
    (kB3,kD3),  w2 + σB ⟷ w2σB
    (kB4,kD4),  w2σB + v ⟷ w2v + σB
    (kB5,kD5),  vP + phos ⟷ vPp
    kP,         vPp ⟶ v + phos
    v0*((1+F*σB)/(K+σB)),     ∅ ⟶ σB
    λW*v0*((1+F*σB)/(K+σB)),  ∅ ⟶ w
    λV*v0*((1+F*σB)/(K+σB)),  ∅ ⟶ v
end kBw kDw kD kB1 kB2 kB3 kB4 kB5 kD1 kD2 kD3 kD4 kD5 kK1 kK2 kP kDeg v0 F K λW λV pTot;
fixed_conc[10] = [:(vPp+phos=pTot)]

max_stars_repo_name: TorkelE/BRNEquilibrate

max_issues_repo_name: TorkelE/BRNEquilibrate

max_forks_repo_name: TorkelE/BRNEquilibrate


text:
### Measurements.jl ---  Uncertainty propagation library
#
# Copyright (C) 2016 Mosè Giordano.
#
# Maintainer: Mosè Giordano <mose AT gnu DOT org>
# Keywords: uncertainty, error propagation, physics
#
# This file is a part of Measurements.jl.
#
# License is MIT "Expat".
#
### Commentary:
#
# This file is the entry point of the package: it defines the new `Measurement'
# type, some functions to handle it within Julia and the new functions provided
# by the package and exposed to users.
#
### Code:

__precompile__()

module Measurements

# This is used to calculate numerical derivatives in "@uncertain" macro.
using Calculus

# Function to handle new type
import Base: show

# Functions provided by this package and exposed to users
export Measurement, measurement, ±

# Define the "Derivatives" type, used inside "Measurement" type.  This should be
# a lightweight and immutable dictionary.
include("derivatives-type.jl")

##### New Type: Measurement
# Definition.  The Measurement type is composed by the following fields:
#   * val: the nominal value of the measurement
#   * err: the uncertainty, assumed to be standard deviation
#   * tag: a (hopefully) unique identifier, it is used to identify a specific
#     measurement in the list of derivatives.  This is usually created with
#     `rand'.
#   * der: the list of derivates.  It is a lightweight dictionary, whose keys
#     are the tuples (nominal value, uncertainty, tag) of all independent
#     variables from which the object has been derived, the corresponding value
#     is the partial derivative of the object with respect to that independent
#     variable.  This dictionary is useful to trace the contribution of each
#     measurement and propagate the uncertainty in the case of functions with
#     more than one argument (in order to deal with correlation between
#     arguments).
immutable Measurement{T<:AbstractFloat} <: AbstractFloat
    val::T
    err::T
    tag::Float64
    der::Derivatives{Tuple{T, T, Float64}, T}
end

# The constructor that users are going to use.  As a Julia convention, the
# lowercase version of a type constructor does something more than the
# constructor itself.
function measurement(val::Real, err::Real=zero(float(val)))
    val, err, der = promote(float(val), float(err), one(float(val)))
    tag = rand()
    return Measurement(val, err, tag, Derivatives((val, err, tag)=>der))
end
@vectorize_2arg Real measurement
const ± = measurement

# Type representation
function show(io::IO, measure::Measurement)
    print(io, measure.val, " ± ", measure.err)
end

include("conversions.jl")
include("comparisons-tests.jl")
include("utils.jl")
include("math.jl")

end # module


text:
# strangle_stu.py
#
# This is a MOS model: Mudd Open Source only.
# This is a deep strangle model based upon the Taboga concept (student version).
# Much of the documentation can be found in the core Taboga model or finutil.py
# Designed in Python 3.6 on for Economics 136
# Original version developed on May 30, 2017 by Prof Evans
# This is version 2.1stu, dated April 11, 2017
#

# This model draws upon utility model finutil_stu. That model must either be in
# the same directory as this model or the path command below must be used to
# point to its location.
# import sys
# sys.path.append('c:/Users/Prof Gary Evans/Dropbox/PyGo/PyFi')
import finutil_stu
fu = finutil_stu
# This version calculates IDV for the put and call, but also requires that you
# supply historical volatility (252 and 30 day). You can include an earnings
# date if date is relevant.
#
# Load data
aname = "Prof Evans"
stosym = "MSFT"
tnow = Dates.today()
# Set earnings to True if this is an earnings-based strangle.
earnings = true
# If you don't need to calculate the days to expiry and earnings (optional)
# then calcdays is False (default).
calcdays = false
# BEGIN of optional time data. Load below ONLY if calcdays = True
# Expiration date:
exyear = 2018
exmonth = 2
exday = 16
# Earnings date (load only if relevant):
eayear = 2018
eamonth = 1
eaday = 31
# END of optional time data load
#
# This model does not use the alpha.
# stock = 15.15
# callstrike = 16
# callbid = 0.90
# callask = 0.93
# putstrike = 15
# putbid = 1.150
# putask = 1.18
# sigma252 = 0.043534
# sigma30 = 0.025
###########################
stock = 166.86
callstrike = 170
callbid = 3.20
callask = 3.25
putstrike = 165
putbid = 4.50
putask = 4.60
sigma252 = 0.011358
sigma30 = 0.016387
if calcdays
    days = fu.daysto(exyear,exmonth,exday)
    if earnings
        daystoearn = fu.daysto(eayear,eamonth,eaday)
    end
else
    # Here you have to supply the days if you don't calculate them.
    days = 14.0
    if earnings
        daystoearn = 11.0
    end
end
#
moname = fu.monthname(exmonth)
# Set effective options prices (PEG)
peg = 0.6
cspread = callask - callbid
callpr = callbid + (peg*cspread)
pspread = putask - putbid
putpr = putbid + (peg*pspread)
cost = callpr + putpr
# Note: callbe and putbe should be calced based upon the prices they would
# have to be the day after earnings, especially for more distant options.
callbe = stock + cost
putbe = stock - cost
callbegr = log(callbe/stock)*100
putbegr = log(putbe/stock)*100
# Consider adding log percent breakeven to current value
#
# Calculate call and put ideal prices based upon historical volatilities.
hv252 = sigma252*fu.durvol(days)
cop252 = fu.otranche(stock,callstrike,hv252,true)
pop252 = fu.otranche(stock,putstrike,hv252,false)
hv30 = sigma30*fu.durvol(days)
cop30 = fu.otranche(stock,callstrike,hv30,true)
pop30 = fu.otranche(stock,putstrike,hv30,false)
#
# Calculate call and put IDV based upon actual prices
# See the documentation for finutil and tabogacidv to see what these are doing.
# Start with the call:
#
seedsigma = 1e-6
cutoff = 1e-4
tempop = 0.00
price2 = 0.00
derive = 0.00
cdurseed = seedsigma*fu.durvol(days)
cdaysigma = 0.06
cdursigma = cdaysigma*fu.durvol(days)
#
# The convergence shown here (Newton's method) was designed by Alec Griffith '17
#
while abs(tempop - callpr) > cutoff
    tempop = fu.otranche(stock,callstrike,cdursigma,true)
    price2 = fu.otranche(stock,callstrike,cdursigma+cdurseed,true)
    deriv = (price2-tempop)/seedsigma
    cdaysigma -= (tempop-callpr)/deriv
    cdursigma = cdaysigma*fu.durvol(days)
end
# End of Loop!
#
# Repeat for the put:
#
seedsigma = 1e-6
cutoff = 1e-4
tempop = 0.00
price2 = 0.00
derive = 0.00
pdurseed = seedsigma*fu.durvol(days)
pdaysigma = 0.01
pdursigma = pdaysigma*fu.durvol(days)
call = false
#
#
while abs(tempop - putpr) > cutoff
    tempop = fu.otranche(stock,putstrike,pdursigma,false)
    price2 = fu.otranche(stock,putstrike,pdursigma+pdurseed,false)
    deriv = (price2-tempop)/seedsigma
    pdaysigma -= (tempop-putpr)/deriv
    pdursigma = pdaysigma*fu.durvol(days)
end
# End of Loop!
callsigmaratio = cdaysigma/sigma252
putsigmaratio = pdaysigma/sigma252
#
dateString = Dates.format(now(), "e u d y")
println()
println("Strangle Model")
println("Analyst's name: $(aname)")
println( "Date: $(dateString)")
println( "Stock symbol: $(stosym)")
println("Option chain: $(moname), $(exday)")
println("Days to expiry: $(days)")
println( "Stock price: $(stock)")
println( "Call strike price: $(callstrike)")
println( "Call Bid: $(callbid)")
println( "Call Ask: $(callask)")
println( "Call price: $(callpr)")
println( "Put strike price: $(putstrike)")
println( "Put Bid: $(putbid)")
println( "Put Ask: $(putask)")
println( "Put price: $(putpr)")
println( "Position cost: $(cost)")
if earnings
    println("Days to earnings: $(daystoearn)")
println( "Call breakeven price: $(callbe)")
println( "Call breakeven ln percent: $(callbegr)")
println( "Put breakeven price: $(putbe)")
println( "Put breakeven ln percent: $(putbegr)")
println( "Call value based upon 252-day HV ($(sigma252)) : $(cop252)")
println( "Call value based upon 30-day HV ($(sigma30)) : $(cop30)")
println( "Put value based upon 252-day HV ($(sigma252)) : $(pop252)")
println( "Put value based upon 30-day HV ($(sigma30)) : $(pop30)")
println( "Call implied duration volatility: $(cdursigma)")
println( "Put implied duration volatility: $(pdursigma)")
println("CRITICAL")
println( "Call implied daily volatility: $(cdaysigma)")
println( "Put implied daily volatility: $(pdaysigma)")
println( "Call sigma ratio: $(callsigmaratio)")
println( "Put sigma ratio: $(putsigmaratio)")
#
# Calculate one day time decay using our calculated IDV.
#
if days >= 2.0
    days = days - 1.0
    cdursigma = cdaysigma*fu.durvol(days)
    cprice1d = fu.otranche(stock,callstrike,cdursigma,true)
    pdursigma = pdaysigma*fu.durvol(days)
    pprice1d = fu.otranche(stock,putstrike,pdursigma,false)
    timedecay = (callpr - cprice1d) + (putpr - pprice1d)
    println( "One day time decay: $(timedecay)")
else
    println( "No time decay with one day remaining.")
end
end
#


text:
using GeometryBasics
using StaticArrays
using LinearAlgebra: normalize, ⋅

####
#### Type definitions
####

# type defs from Julius Krumbiegels PR
# https://github.com/JuliaPlots/Makie.jl/pull/979
abstract type AbstractPath{PT<:AbstractPoint} end

# simple lines get there own type. This helps with type stability
# for graphs without curvy edges and improves performance
struct Line{PT} <: AbstractPath{PT}
    p0::PT
    p::PT
end

abstract type PathCommand{PT<:AbstractPoint} end

struct BezierPath{PT} <: AbstractPath{PT}
    commands::Vector{PathCommand{PT}}
end

struct MoveTo{PT} <: PathCommand{PT}
    p::PT
end

struct LineTo{PT} <: PathCommand{PT}
    p::PT
end

struct CurveTo{PT} <: PathCommand{PT}
    c1::PT
    c2::PT
    p::PT
end

ptype(::Union{AbstractPath{PT}, Type{<:AbstractPath{PT}}}) where {PT} = PT

####
#### Helper functions to work with bezier pathes
####

"""
    interpolate(p::AbstractPath, t)

Parametrize path `p` from `t ∈ [0, 1]`. Return postion at `t`.

TODO: Points are not necessarily evenly spaced!
"""
function interpolate(p::BezierPath{PT}, t) where PT
    @assert p.commands[begin] isa MoveTo
    N = length(p.commands) - 1

    tn = N*t
    seg = min(floor(Int, tn), N-1)
    tseg = tn - seg

    p0 = p.commands[seg+1].p
    return _interpolate(p.commands[seg+2], p0, tseg)
end

interpolate(l::Line{PT}, t) where PT = l.p0 + t*(l.p - l.p0) |> PT

"""
    tangent(p::AbstractPath, t)

Parametrize path `p` from `t ∈ [0, 1]`. Return tangent at `t`.
"""
function tangent(p::BezierPath, t)
    @assert p.commands[begin] isa MoveTo
    N = length(p.commands) - 1

    tn = N*t
    seg = min(floor(Int, tn), N-1)
    tseg = tn - seg

    p0 = p.commands[seg+1].p
    return _tangent(p.commands[seg+2], p0, tseg)
end
tangent(l::Line, _) = normalize(l.p - l.p0)

"""
    discretize(path::AbstractPath)

Return vector of points which represent the given `path`.
"""
function discretize(path::BezierPath{T}) where {T}
    v = Vector{T}()
    for c in path.commands
        _discretize!(v, c)
    end
    return v
end

discretize(l::Line) = [l.p0, l.p]

"""
    _interpolate(c::PathCommand, p0, t)

Returns positions along the path `c` starting from `p0` in range `t ∈ [0, 1]`.
"""
_interpolate(c::LineTo{PT}, p0, t) where {PT} = p0 + t*(c.p - p0) |> PT
function _interpolate(c::CurveTo{PT}, p0, t) where {PT}
    p1, p2, p3 = c.c1, c.c2, c.p
    (1 - t)^3 * p0 + 3(t - 2t^2 + t^3) * p1 + 3(t^2 -t^3) * p2 + t^3 * p3 |> PT
end

"""
    _tangent(c::PathCommand, p0, t)

Returns tangent vector along the path `c` starting from `p0` in range `t ∈ [0, 1]`.
"""
_tangent(c::LineTo, p0, _) = normalize(c.p - p0)
function _tangent(c::CurveTo{PT}, p0, t) where PT
    p1, p2, p3 = c.c1, c.c2, c.p
    normalize(-3(1 - t)^2 * p0 + 3(1 - 4t + 3t^2) * p1 + 3(2t -3t^2) * p2 + 3t^2 * p3) |> PT
end

"""
    _discretize!(v::Vector{AbstractPint}, c::PathCommand)

Append interpolated points of path `c` to pos vector `v`
"""
_discretize!(v::Vector{<:AbstractPoint}, c::Union{MoveTo, LineTo}) = push!(v, c.p)
function _discretize!(v::Vector{<:AbstractPoint}, c::CurveTo)
    N0 = length(v)
    p0 = v[end]
    N = 60 # TODO: magic number of points for discrtization
    resize!(v, N0 + N)
    dt = 1.0/N
    for (i, t) in enumerate(dt:dt:1.0)
        v[N0 + i] = _interpolate(c, p0, t)
    end
end

"""
    waypoints(p::BezierPath)

Returns all the characteristic points of the path. For debug reasons.
"""
function waypoints(p::BezierPath{PT}) where {PT}
    v = PT[]
    for c in p.commands
        if c isa CurveTo
            push!(v, c.c1)
            push!(v, c.c2)
        end
        push!(v, c.p)
    end
    return v
end

"""
    isline(p::AbstractPath)

True if the AbstractPath just represents a straight line.
"""
isline(p::Line) = true
isline(p::BezierPath) = length(p.commands)==2 && p.commands[1] isa MoveTo && p.commands[2] isa LineTo

####
#### Special constructors to create abstract pathes
####

"""
    Path(P::Vararg{PT, N}; tangents, tfactor=.5) where {PT<:AbstractPoint, N}

Create a bezier path by natural cubic spline interpolation of the points `P`.
If there are only two points and no tangents return a straight line.

The `tangents` kw allows you pass two vectors as tangents for the first and the
last point. The `tfactor` affects the curvature on the start and end given some
tangents.
"""
function Path(P::Vararg{PT, N}; tangents=nothing, tfactor=.5) where {PT<:AbstractPoint, N}
    @assert N>2

    # cubic_spline will work for each dimension separatly
    pxyz = _cubic_spline(map(p -> p[1], P)) # get first dimension
    for i in 2:length(PT) # append all other dims
        pxyz = hcat(pxyz, _cubic_spline(map(p -> p[i], P)))
    end

    # create waypoints from waypoints in separat dementions
    WP = SVector{length(P)-1, PT}(PT(p) for p in eachrow(pxyz))

    commands = Vector{PathCommand{PT}}(undef, N)
    commands[1] = MoveTo(P[1])

    # first command, recalculate WP if tangent is given
    first_wp = WP[1]
    if tangents !== nothing
        p1, p2, t = P[1], P[2], normalize(tangents[1])
        dir = p2 - p1
        d = tfactor * norm(dir ⋅ t)
        first_wp = PT(p1+d*t)
    end
    commands[2] = CurveTo(first_wp,
                          2*P[2] - WP[2],
                          P[2])
    # middle commands
    for i in 3:(N-1)
        commands[i] = CurveTo(WP[i-1],
                              2*P[i] - WP[i],
                              P[i])
    end
    # last command, recalculate last WP if tangent is given
    last_wp = (P[N] + WP[N-1])/2
    if tangents !== nothing
        p1, p2, t = P[N-1], P[N], normalize(tangents[2])
        dir = p2 - p1
        d = tfactor * norm(dir ⋅ t)
        last_wp = PT(p2-d*t)
    end
    commands[N] = CurveTo(WP[N-1],
                          last_wp,
                          P[N])

    BezierPath(commands)
end

# same function as above but for just 2 points specificially
function Path(P::Vararg{PT, 2}; tangents=nothing, tfactor=.5) where {PT<:AbstractPoint}
    if tfactor isa NTuple{2, <:Number}
        tf1, tf2 = tfactor
    else
        tf1 = tf2 = tfactor
    end

    p1, p2 = P
    if tangents === nothing
        return Line(p1, p2)
    else
        t1 = normalize(Pointf(tangents[1]))
        t2 = normalize(Pointf(tangents[2]))
        len = norm(p2 - p1)
        return BezierPath([MoveTo(p1),
                           CurveTo(PT(p1+len*tf1*t1),
                                   PT(p2-len*tf2*t2),
                                   p2)])
    end
end

"""
    Path(radius::Real, p::Vararg{PT, N}) where {PT<:AbstractPoint,N}

Draw straight lines through the points `p`. Within `radius` of each
point the line will be smoothly connected.
"""
function Path(radius::Real, p::Vararg{PT, N}) where {PT<:AbstractPoint,N}
    if iszero(radius)
        commands = PathCommand{PT}[MoveTo(PT(p[1]))]
        for pos in p[2:end]
            push!(commands, LineTo(PT(pos)))
        end
        return BezierPath(commands)
    else
        points = collect(reverse(p))
        pos = pop!(points)
        commands = PathCommand{PT}[MoveTo(PT(pos))]

        while length(points) > 1
            mid = pop!(points)
            dst = points[end] # due to reverse...
            dir1 = normalize(mid - pos)
            dir2 = normalize(dst - mid)
            r1 = mid - radius * dir1
            r2 = mid + radius * dir2
            c1 = mid - .25 * radius * dir1
            c2 = mid + .25 * radius * dir2
            push!(commands, LineTo(PT(r1)))
            push!(commands, CurveTo(PT(c1),
                                    PT(c2),
                                    PT(r2)))
            pos = r2
        end

        push!(commands, LineTo(points[end]))
        return BezierPath(commands)
    end
end

"""
    _cubic_spline(p)

Given a number of points in one dimension calculate waypoints between them.

    _cubic_spline([x1, x2, x3])

Will return the x coordinates of the waypoints `wp1` and `wp2`.
Those are the first waypoints between in the cubic bezier sense.
"""
function _cubic_spline(p)
    N = length(p) - 1

    M = SMatrix{N,N}(if i==j # diagonal
                         if i==1; 2; elseif i==N; 7; else 4 end
                     elseif i==j+1 # lower
                         if i==N; 2; else 1 end
                     elseif i==j-1 # upper
                         1
                     else
                         0
                     end for i in 1:N, j in 1:N)

    b = SVector{N}(if i == 1
                       p[i] + 2p[i+1]
                   elseif i == N
                       8p[i] + p[i+1]
                   else
                       4p[i] + 2p[i+1]
                   end for i in 1:N)
    return M \ b
end


text:
#Stochastic collocation towards a polynomial
using Polynomials
using Roots

import AQFED.Math: normcdf, normpdf, norminv, FitResult
import AQFED.Black: blackScholesFormula, blackScholesVega
using AQFED.Bachelier
#using MINPACK #slower
using LeastSquaresOptim
#using SCS #slower
export Polynomial, IsotonicCollocation, solveStrike, priceEuropean, density, makeIsotonicCollocation, weightedPrices

struct IsotonicCollocation{T,U}
    p1::AbstractPolynomial{T}
    p2::AbstractPolynomial{T}
    forward::U
    # coeffs::AbstractArray{T} #Polynomials.Polynomial{T,:x}(Val(false), coeffs)
    #  function IsotonicCollocation{T, U}(p1::AbstractPolynomial{T},
    #      p2::AbstractPolynomial{T},
    #      forward::U) where {T, U}
    #      new{T,U}(p1, p2, forward, zeros(T,degree(p1)+degree(p2)+3))
    #  end
end

Base.broadcastable(p::IsotonicCollocation) = Ref(p)

function Polynomial(iso::IsotonicCollocation; minSlope = 1e-5)
    p = integrate(iso.p1^2 + iso.p2^2 + minSlope * iso.forward) #add a minimum slope such that density is not too spiky and collocation not too flat
    theoForward = hermiteIntegral(p)
    p[0] = iso.forward - theoForward
    #fixing the forward is important to avoid the case where C(K=0) < 0.
    return p
end

function solveStrike(p::AbstractPolynomial{T}, pd::AbstractPolynomial{T},pd2::AbstractPolynomial{T}, strike::U)::T where {T,U}
    function objHalley(x::W)::Tuple{W,W,W} where {W}
        u = p(x)-strike
        du =  pd(x)
        d2u = pd2(x)
        return (u, u / du, d2u == 0 ? Inf : du / d2u)
    end
    guess = 0.0
    if abs(p[1]) > eps(strike)
        guess = (strike - p[0]) / p[1]
    end
    try
        return find_zero(
            objHalley,
            guess,
            Roots.SuperHalley(), #seems to be (much) more robust around the spike in the density.
            atol = 100 * eps(strike),
            maxevals = 16,
            verbose = false,
        )
    catch err
        function obj(x::W)::W where {W}
            p(x) - strike
        end
        return find_zero(
            obj,
            (guess-10,guess+10),
            Bisection(),
            atol = 100 * eps(strike)
        )
    end
end
function solveStrike(p::AbstractPolynomial{T}, strike::U; useHalley = true)::T where {T,U}
    if useHalley
        pd = derivative(p)
        pd2 = derivative(p, 2)
        solveStrike(p, pd, pd2, strike)
    else
        pk = p - strike
        r = roots(pk)
        for ri in r
            if abs(imag(ri)) < eps(strike)
                return real(ri)
            end
        end
        throw(DomainError(strike, "no real roots at the given strike"))
    end
end

function priceEuropean(p::AbstractPolynomial{T}, pd::AbstractPolynomial{T}, pd2::AbstractPolynomial{T}, isCall::Bool, strike::U, forward::U, discountDf::U)::T where {T,U}
    ck = solveStrike(p, pd, pd2, strike)
    valuef = hermiteIntegralBounded(p, ck)
    valuek = normcdf(-ck)
    callPrice = valuef - strike * valuek
    putPrice = -(forward - strike) + callPrice
    if isCall
        return callPrice * discountDf
    else
        return putPrice * discountDf
    end
end

function priceEuropean(p::AbstractPolynomial{T}, isCall::Bool, strike::U, forward::U, discountDf::U)::T where {T,U}
    ck = solveStrike(p, strike)
    valuef = hermiteIntegralBounded(p, ck)
    valuek = normcdf(-ck)
    callPrice = valuef - strike * valuek
    putPrice = -(forward - strike) + callPrice
    if isCall
        return callPrice * discountDf
    else
        return putPrice * discountDf
    end
end

function density(p::AbstractPolynomial, strike::Number)::Number
    ck = solveStrike(p, strike)
    return normpdf(ck) / derivative(p)(ck)
end

function rawMoment(p::AbstractPolynomial, moment::Int)::Number
    if moment == 0
        return 1.0
    end
    q = p
    for i = 2:moment
        q *= p
    end
    return hermiteIntegral(q)
end

#return mean, standard dev, skew, kurtosis of the collocation
function stats(p::AbstractPolynomial)
    μ = hermiteIntegral(p)
    μ2 = hermiteIntegral((p - μ)^2)
    μ3 = hermiteIntegral((p - μ)^3)
    μ4 = hermiteIntegral((p - μ)^4)

    skew = μ3 / μ2^1.5
    kurtosis = μ4 / μ2^2
    return μ, sqrt(μ2), skew, kurtosis
end
function hermiteIntegral(p::AbstractPolynomial{T})::T where {T}
    m0 = 1.0
    nx0 = 0.0
    sum = 0.0
    sum += p[0] * m0
    m1 = nx0
    sum += p[1] * m1
    for i = 2:degree(p)
        m2 = m0 * (i - 1)
        sum += p[i] * m2
        m0 = m1
        m1 = m2
    end
    return sum
end

function hermiteIntegralBounded(p::AbstractPolynomial{T}, ck::U)::T where {T,U}
    x0 = ck
    m0 = normcdf(-x0)
    nx0 = normpdf(x0)
    sum = 0.0
    sum += p[0] * m0
    m1 = nx0
    sum += p[1] * m1
    for i = 2:degree(p)
        m2 = m0 * (i - 1) + nx0 * x0^(i - 1)
        sum += p[i] * m2
        m0 = m1
        m1 = m2
    end
    return sum
end

function weightedPrices(
    isCall::Bool,
    strikes::Vector{T},
    vols::Vector{T},
    weights::Vector{T},
    forward::T,
    discountDf::T,
    tte::T;
    vegaFloor = T(1e-3), #vega floored at 1e-3*forward
)::Tuple{Vector{T},Vector{T}} where {T}
    prices = @. blackScholesFormula(isCall, strikes, forward, vols^2 * tte, 1.0, discountDf)
    w = @. weights / max(vegaFloor * forward, blackScholesVega(strikes, forward, vols^2 * tte, 1.0, discountDf, tte))
    return prices, w
end

function makeIsotonicCollocationGuess(
    strikes::Vector{T},
    callPrices::Vector{T},
    weights::Vector{T},
    τ::T,
    forward::T,
    discountDf::T;
    deg = 3, #1 = Bachelier (trivial, smoother if least squares iterations small). otherwise 3 is good.
) where {T}
    if deg >= 3
        strikesf, pricesf, weightsf = filterConvexPrices(strikes, callPrices ./ discountDf, weights, forward)
        strikesf, pif, xif = makeXFromUndiscountedPrices(strikesf, pricesf)
        cubic = Polynomials.fit(xif, strikesf, 3; weights = weightsf)
        if !isCubicMonotone(cubic)
            cubic = LeastSquaresCubicMurrayForward(xif, strikesf, weightsf, forward)
        end
        if deg == 3
            return IsotonicCollocation(cubic, forward)
        else
            return fitMonotonic(xif, strikesf, weightsf, forward, cubic, deg = deg)
        end
    else
        return fitBachelier(strikes, callPrices, weights, τ, forward, discountDf)
    end
end

function makeIsotonicCollocation(
    strikes::Vector{T},
    callPrices::Vector{T},
    weights::Vector{T},
    τ::T,
    forward::T,
    discountDf::T;
    deg = 3, #degree of collocation. 5 is usually best from a stability perspective.
    degGuess = 3, #1 = Bachelier (trivial, smoother if least squares iterations small). otherwise 3 is good.
    minSlope = 1e-6,
    penalty = 0.0,
)::Tuple{IsotonicCollocation,FitResult} where {T} #return collocation and error measure
    isoc = makeIsotonicCollocationGuess(strikes, callPrices, weights, τ, forward, discountDf, deg = degGuess)
    #optimize towards actual prices
    isoc, m =
        fit(isoc, strikes, callPrices, weights, forward, discountDf, deg = deg, minSlope = minSlope, penalty = penalty)
    return isoc, m
end

function fitBachelier(strikes, prices, weights, τ, forward, discountDf)
    m = length(strikes)
    i = findfirst(x -> x > forward, strikes)
    if i == nothing
        i = m
    elseif i == 1
        i = 2
    end
    price = (prices[i] * (forward - strikes[i-1]) + prices[i-1] * (strikes[i] - forward)) / (strikes[i] - strikes[i-1])
    strike = forward
    bvol = Bachelier.bachelierImpliedVolatility(price, true, strike, τ, forward, discountDf)
    isoc = IsotonicCollocation(Polynomials.Polynomial([sqrt(bvol * sqrt(τ))]), Polynomials.Polynomial([0.0]), forward)
    return isoc
end

function fit(
    isoc::IsotonicCollocation,
    strikes::AbstractArray{U},
    prices::AbstractArray{U},
    weights::AbstractArray{U},
    forward::U,
    discountDf::U;
    deg = 3,
    minSlope = 1e-4,
    penalty = 0.0,
) where {U}
    q = trunc(Int, (deg + 1) / 2)
    iter = 0
    function obj!(fvec::Z, c::AbstractArray{W})::Z where {Z,W}
        p1 = Polynomials.Polynomial(c[1:q])
        p2 = Polynomials.Polynomial(c[q+1:2*q])
        isoc = IsotonicCollocation(p1, p2, forward)
        p = Polynomial(isoc, minSlope = minSlope)
        pd = derivative(p,1)
        pd2 = derivative(p,2)
        iter += 1
        n = length(strikes)
        @. fvec[1:n] = weights * (priceEuropean(p, pd, pd2, true, strikes, forward, discountDf) - prices)
        if penalty > 0
            ip = hermiteIntegral(derivative(p, 2)^2)
            pvalue = penalty * ip
            fvec[n+1] = pvalue
        end
        fvec
    end
    c0 = zeros(Float64, 2 * q)
    c1 = Polynomials.coeffs(isoc.p1)
    for i = 1:min(q, length(c1))
        c0[i] = c1[i]
    end
    c2 = Polynomials.coeffs(isoc.p2)
    for i = 1:min(q - 1, length(c2))
        c0[q+i] = c2[i]
    end
    for i = 1:length(c0)
        if c0[i] == 0
            c0[i] = eps(forward)  #zero would break automatic differentiation
        end
    end
    outlen = length(strikes)
    if penalty > 0
        outlen += 1
    end
    #TODO optimize the number of allocations. This would mean to abandon Polynomials package and (re)use Vectors instead
    fit = optimize!(
        LeastSquaresProblem(x = c0, f! = obj!, autodiff = :forward,
        # g! = jac!, #useful to debug issue with ForwardDiff NaNs
        output_length = outlen),
        LevenbergMarquardt();
        iterations = deg*300,
    )
    fvec = zeros(Float64, outlen)
    obj!(fvec,fit.minimizer)
    #fit = optimize(obj, c0, LevenbergMarquardt(); show_trace = false, autodiff = :forward, iterations = deg * 300) #autodiff breaks without Halley. Would need custom
    # function obj!(fvec, x)
    #     fvec[:] = obj(x)
    #     fvec
    # end
    # fit = fsolve(obj!, c0, length(strikes); show_trace=true, method=:lm, tol=1e-10) #fit.x
    #println(iter, " fit ", fit)
    c0 = fit.minimizer
    measure = fit.ssr
    return IsotonicCollocation(Polynomials.Polynomial(c0[1:q]), Polynomials.Polynomial(c0[q+1:2*q]), forward),  FitResult(fit.ssr, iter, fit.minimizer, fvec, fit)
end

function IsotonicCollocation(cubic::AbstractPolynomial, forward::Number)
    if degree(cubic) > 3
        throw(DomainError(degree(cubic), "expected a cubic"))
    end
    e = (cubic[1] - cubic[2]^2 / (3 * cubic[3]))
    e = max(e, forward * 1e-6)
    p2 = Polynomials.Polynomial([sqrt(e)])
    p1 = Polynomials.Polynomial([cubic[2] / sqrt(3 * cubic[3]), sqrt(3 * cubic[3])])
    return IsotonicCollocation(p1, p2, forward)
end

function fitMonotonic(xif, strikesf, w1, forward, cubic; deg = 3)
    xmin = -3.0
    if deg > 3 && xif[1] > xmin && strikesf[1] > 0
        #stabilize fit
        xif = [xmin; xif]
        strikesf = [0.0; strikesf]
        w1 = [w1[1]; w1]
    end
    q = trunc(Int, (deg + 1) / 2)
    iter = 0
    function obj(c)
        p1 = Polynomials.Polynomial(c[1:q])
        p2 = Polynomials.Polynomial(c[q+1:2*q])
        isoc = IsotonicCollocation(p1, p2, forward)
        p = Polynomial(isoc)
        iter += 1
        return @. w1 * (p(xif) - strikesf)
    end

    isocubic = IsotonicCollocation(cubic, forward)
    c1 = Polynomials.coeffs(isocubic.p1)
    c2 = Polynomials.coeffs(isocubic.p2)
    c0 = zeros(Float64, 2 * q)
    c0[1:min(q, length(c1))] = c1
    c0[q+1:min(2 * q, q + length(c2))] = c2
    xamax = max(abs(xif[1]), abs(xif[end]))
    for i = 1:length(c0)
        if c0[i] == 0
            c0[i] = eps(forward) / (xamax)^(i - 1)  #zero would break automatic differentiation
        end
    end
    fit = optimize(obj, c0, LevenbergMarquardt(); show_trace = false, autodiff = :forward)
    c0 = fit.minimizer
    isoc = IsotonicCollocation(Polynomials.Polynomial(c0[1:q]), Polynomials.Polynomial(c0[q+1:2*q]), forward)
    #println(iter, " fitMonotonic ", Polynomial(isoc), " ", fit)
    return isoc
end

function isCubicMonotone(cubic::AbstractPolynomial)::Bool
    c = Polynomials.coeffs(cubic)
    if c[4] == 0 && c[3] == 0
        return c[2] > 0
    else
        return (c[3]^2 - 3 * c[2] * c[4] < 0) && (c[4] > 0)
    end
end

#LeastSquaresMurray retuns cubic coefficients a_0 + a_1 x + a_3 x^3 (a_2=0) of the least squares cubic going through x with weights w
function LeastSquaresCubicMurrayForward(x::Vector{T}, y::Vector{T}, w::Vector{T}, forward::T) where {T}
    rhs = [sum(@. w * x * y), sum(@. w * x^3 * y)]
    m = [sum(@. w * x^2) sum(@. w * x^4); sum(@. w * x^4) sum(@. w * x^6)]
    result = m \ rhs #solve linear system
    return Polynomials.Polynomial([forward, abs(result[1]), zero(T), abs(result[2])])
end


text:
module FlowWorkspace

    using GigaSOM,DataFrames
    using EzXML,MetaGraphs

    using PolygonOps, StaticArrays
    using Glob: glob
    using HTTP: escapeuri

    include("Groups.jl")
    include("Gating.jl")
    include("Utils.jl")

    export load
end


text:
module forprop

relu(Z) = max(0, Z)

function softmax(Z)
    exped = exp.(Z)
    return exped ./ sum(exped, dims = 1)
end

sigmoid(Z) = 1 / (1 + exp(-Z))

end

max_stars_repo_name: prenc/MiraculousNN

max_issues_repo_name: prenc/MiraculousNN

max_forks_repo_name: prenc/MiraculousNN


text:
using Literate

# directory where the markdown files are put
md_dir = joinpath(@__DIR__, "md")
Literate.markdown("car_travels.jl", md_dir, execute=true, documenter=false, credit=false)


text:
using Test: Serialization
using Test
using Graphs
using SimpleWeightedGraphs
using Serialization
using ..ChemistryFeaturization: FeaturizedAtoms

@testset "FeaturizedAtoms" begin

    local fnames = ["X", "Block", "Atomic mass"]

    # TODO - Actually test the `featurized_atoms.encoded_features` value. Probably define a custom type for this?
    @testset "Construct and encode" begin

        efds = ElementFeatureDescriptor.(fnames)
        fzn = GraphNodeFeaturization(efds)

        triangle_C = AtomGraph(Float32.([0 1 1; 1 0 1; 1 1 0]), ["C", "C", "C"])

        featurized_atoms = FeaturizedAtoms(triangle_C, fzn)

        # @test sef.encoded_features == ???
    end

    @testset "Decode" begin
        fzn = GraphNodeFeaturization(fnames, nbins = 2)
        F2 = AtomGraph(Float32.([0 1; 1 0]), ["F", "F"])
        featurized_atoms = featurize(F2, fzn)

        decoded_matrix = decode(fzn, featurized_atoms.encoded_features)

        decoded_ag = decode(featurized_atoms)
        @test all(
            map(d -> d[1]["Block"] == d[2]["Block"] == "p", [decoded_matrix, decoded_ag]),
        )

        # check if `encoded_features` are generated correctly using encode itself
        fzn2 = GraphNodeFeaturization(fnames, nbins = [2, 4, 2])
        @test all(encode(fzn2, F2) .== featurized_atoms.encoded_features)

    end

end


text:
###################################################
## RoutingPaths type
###################################################

"""
    Stores path information in network
    - for now, just compute the shortest paths
    - if no_sp = true, returns an incomplete object (in case want to run shortest paths later)
"""
mutable struct RoutingPaths
    # Required attributes
    "A reference to the network"
    network::Network

    "road timing list: usually sparse matrix"
    times::AbstractArray{Float64,2}

    # optional attributes
    "if computed, contains timings from anywhere to anywhere"
    pathTimes::AbstractArray{Float64,2}
    "if computed, contains information to recontruct the paths"
    pathPrevious::AbstractArray{Int64,2}

    "Constructor: reference to network and timing information"
    function RoutingPaths(n::Network, times::AbstractArray{Float64,2}; no_sp=false)
        # a few tests
        if size(times) != (length(n.nodes),length(n.nodes))
            error("The given timings do not fit the network")
        end
        for (i,j) in keys(n.roads)
            if times[i,j] < 0.
                error("$i => $j : negative timing $(times[i,j])")
            end
        end
        obj = new()
        obj.network = n
        obj.times = times

        if !no_sp
            if nworkers() > 1
                parallelShortestPaths!(obj)
            else
                shortestPaths!(obj)
            end
        end
        return obj
    end
end

"""
    Without timings, uses distances
"""
RoutingPaths(n::Network) = RoutingPaths(n,roadDistances(n))

function Base.show(io::IO, r::RoutingPaths)
    g = r.network.graph
    println("Network routing times")
    if isdefined(r,:pathTimes)
        println(io, "with path and time for all $(nv(g)*(nv(g)-1)) possible trips")
    else
        println(io, "for all $(ne(g)) roads (no path computed yet)")
    end
end

"""
    Returns a RoutingTimes object where the "timings" are the distances
"""
function roadDistances(n::Network)
    dists = spzeros(length(n.nodes), length(n.nodes))
    for ((i,j),r) in n.roads
        dists[i,j] = r.distance
    end
    return dists
end

"""
    Returns walking road distances (undirected)
"""
function walkingDistances(n::Network)
    dists = roadDistances(n)
    return max(dists,dists')
end

"""
    Gives access to the timing matrix of the graph after path computation
"""
function getPathTimes(r::RoutingPaths)
    if !isdefined(r,:pathTimes)
        error("The paths have not been computed yet")
    end
    return r.pathTimes
end

"""
    `traveltime(r,i,j,k)` : travel time from i to j coming from node k
    if k == i: start in i
"""
traveltime(r::RoutingPaths, i::Int, j::Int, k::Int = i) = r.pathTimes[i,j]

"""
    `getPath(RoutingPaths, orig, dest, prev)`
    Returns path between origin and destination (list of node ids) coming from prev
    if prev = orig then shortest origin
"""
function getPath(r::RoutingPaths, orig::Int, dest::Int, prev::Int = orig)
    path = Int[dest]
    lastNode = dest
    while lastNode != orig
        lastNode = r.pathPrevious[orig,lastNode]
        push!(path, lastNode)
    end
    return path[end:-1:1]
end

"""
    `getPathEdges(RoutingPaths, orig, dest)`
    Returns path between origin and destination as a list of edge IDs (i,j)
"""
function getPathEdges(r::RoutingPaths, orig::Int, dest::Int)
    # empty path case
    if orig == dest
        return Edge[]
    end
    edgeEnd = dest
    edgeStart = r.pathPrevious[orig, edgeEnd]
    path = [Edge(edgeStart, edgeEnd)]
    while edgeStart != orig
        edgeEnd = edgeStart
        edgeStart = r.pathPrevious[orig, edgeEnd]
        push!(path, Edge(edgeStart, edgeEnd))
    end
    # don't forget to reverse path
    return path[end:-1:1]
end

"""
    `getPathWithTime(RoutingPaths, orig, dest, prev)`
    Returns path and times between origin and destination (list of node ids) coming from prev
    if prev = orig then shortest origin
"""
function getPathWithTimes(r::RoutingPaths, orig::Int, dest::Int, prev::Int = orig; startTime::Float64=0.)
    path = Int[dest]
    lastNode = dest
    while lastNode != orig
        lastNode = r.pathPrevious[orig,lastNode]
        push!(path, lastNode)
    end
    path = path[end:-1:1]

    times = Array(Tuple{Float64,Float64},length(path)-1)
    t = startTime
    for i in 1:length(path) - 1
        rdTime = r.times[path[i],path[i+1]]
        times[i] = (t, t+rdTime)
        t += rdTime
    end
    return path, times
end

"""
    Returns the given path time (just sum the link times on path)
"""
function pathTime(times::AbstractArray{Float64,2}, path::Vector{Int})
    time = 0.
    for i in 1:length(path) - 1
        time+=times[path[i],path[i+1]]
    end
    return time
end
pathTime(r::RoutingPaths, path::Vector{Int}) = pathTime(r.times,path)

"""
    Returns the given path time from path edge list (just sum the link times on path)
"""
function pathEdgesTime(times::AbstractArray{Float64,2}, pathEdges::Vector{Edge})
    time = 0.
    for edge in pathEdges
        time += times[src(edge),dst(edge)]
    end
    return time
end
pathEdgesTime(r::RoutingPaths, pathEdges::Vector{Edge}) = pathEdgesTime(r.times,pathEdges)

"""
    `longestPathTime`, longuest path time in network
"""
longestPathTime(r::RoutingPaths) = maximum(r.pathTimes)

max_stars_repo_name: andrew-zhm/RoutingNetworks

max_issues_repo_name: andrew-zhm/RoutingNetworks

max_forks_repo_name: andrew-zhm/RoutingNetworks


text:
module mutex

using PDDL
using Combinatorics

include("types.jl")
include("graph.jl")

export get_all_actions,
       create_graph,
       Action

end


text:
using CSV, DataFrames, StatsPlots, StanSample
using Distributions, Statistics, Random

ProjDir = @__DIR__

# Simulate the data

Random.seed!(123)

N = 15
obspairs = ((30.0, 53), (35.0, 45), (40.0, 28), (45.0, 26), (50.0, 25))

ppu = []; quantity = []
for i in 1:N in
	obs = rand(obspairs, 1)
	append!(ppu, [rand(Normal(obs[1][1], 1.0))])
	append!(quantity, [rand(Normal(obs[1][2], 1.0))])
end

df = DataFrame()
df[!, :ppu] = ppu
df[!, :quantity] = quantity
df[!, :ppu_l] = log.(ppu)
df[!, :quantity_l] = log.(quantity)

CSV.write(joinpath(ProjDir, "data", "data.csv"), df)

opt_01 = "
data {
	int<lower=0> N;
	vector[N] q; // Outcome
	vector[N] p; // Predictor
}
parameters {
	real a;
	real b;
	real<lower=0> sigma;
}
model {
	vector[N] mu;            		// mu is a vector
	a ~ cauchy(0.0, 5.0);			// intercept positive
	b ~ cauchy(0.0, 5.0);			// demand drops as price increases
	sigma ~ uniform(0, 15);
	mu = a + b * (p - mean(p));
	q ~ normal(mu, sigma);
}
";

mean_q_l = mean(df[:, :quantity_l])
mean_p_l = mean(df[:, :ppu_l])
data = Dict(
  :N => size(df, 1),
  :q => df[:, :quantity_l], # .- mean_q_l,
  :p => df[:, :ppu_l]		# .- mean_p_l
)

p1 = scatter(df[:, :ppu], df[:, :quantity],
	xlabel="ppu", ylabel="quantity sold", leg=false)
p2 = scatter(data[:p], data[:q],
	xlabel="log(ppu)", ylabel="log(quantity sold)", leg=false)
plot(p1, p2, layout=(2, 1))
savefig(joinpath(ProjDir, "plots", "data.png"))

sm = SampleModel("opt_01", opt_01)

rc = stan_sample(sm, data=data)

if success(rc)
	chns = read_samples(sm; output_format=:mcmcchains)
	show(chns)
	plot(chns)
	savefig(joinpath(ProjDir, "plots", "chains.png"))

	dfs = read_samples(sm; output_format=:dataframe)

	dfsa = DataFrame(chns)
	#density(dfsa[:, Symbol("mu_pred.1")], label="mu_pred.1")
	#for i in 2:3:N
	#	density!(dfsa[:, Symbol("mu_pred.$i")], label="mu_pred.$i")
	#end
	#savefig(joinpath(ProjDir, "plots", "density.png"))

	p = 3.4:0.01:4.0
	ps = collect(p)
	p1 = scatter(data[:p], data[:q],
		#xlims=(3.4, 4.0), ylims=(3.0, 4.0),
		xlabel="ppu_l", ylabel="quantity_l sold", leg=false)
	a_mean = mean(dfs[:, :a])
	b_mean = mean(dfs[:, :b])
	mu_q = a_mean .+ b_mean .* (p .- mean(p))
	plot!(p1, p, mu_q)
	savefig(joinpath(ProjDir, "plots", "regression.png"))

	price_points_l = data[:p]
	plot(xlabel="Price", ylabel="Quantity", leg=false)
	selected_rows = rand(1:size(dfs, 1), 75)
	for row in eachrow(dfs[selected_rows,:])
		qs = exp.(row[:a] .+ row[:b] .* (ps .- mean(ps)))
		plot!(exp.(p), qs, color=:gray)
	end
	qqs = exp.(a_mean .+ b_mean .* (price_points_l .- mean(ps)))
	scatter!(exp.(price_points_l), qqs)
	savefig(joinpath(ProjDir, "plots", "sampling.png"))

	scatter(dfs[:, :a], dfs[:, :b],
		#xlims = (-0.1, 1.0), ylims= (-10.0, 10.0),
		xlabel="a", ylabel="b", leg=false)
	savefig(joinpath(ProjDir, "plots", "correlations.png"))
end


text:
@testset ExtendedTestSet "squared_euclidean" begin
    @testset ExtendedTestSet "squared_euclidean" begin
        x = CUDA.CuArray(
            boolean_indicator(
                [
                    0 1 1 1 0
                    1 1 1 1 1
                    1 0 0 0 1
                    1 0 0 0 1
                    1 0 0 0 1
                    1 1 1 1 1
                    0 1 1 1 0
                ]
            ),
        )
        dt = CuArray{Float32}(undef, size(x))
        v = CUDA.ones(Int64, size(x))
        z = CUDA.zeros(Float32, size(x) .+ 1)
        tfm = DistanceTransforms.SquaredEuclidean(x, dt, v, z)
        answer = CUDA.CuArray(
            [
                1.0 0.0 0.0 0.0 1.0
                0.0 0.0 0.0 0.0 0.0
                0.0 1.0 1.0 1.0 0.0
                0.0 4.0 4.0 4.0 0.0
                0.0 1.0 1.0 1.0 0.0
                0.0 0.0 0.0 0.0 0.0
                1.0 0.0 0.0 0.0 1.0
            ],
        )
        @test transform!(x, tfm) == answer
    end
end


text:
# 定義の仕方がpython (:) と違って =>　を つかう

fruits = Dict("apple"=> 1, "banana"=> 2, "orange"=>3)
println(fruits)

# アクセスはキー
println(fruits["apple"])

fruits["mango"] = 4
println(fruits)

# 削除
pop!(fruits,"banana")
println(fruits)

# 削除part2
delete!(fruits,"apple")
println(fruits)

# 順序という概念がないから数値インデックスでアクセスはできない
# println(fruits[1])

# 辞書型の操作
println(haskey(fruits, "apple"))

# get(dict, key, default)
# ないときはdefault
println(get(fruits,"apple",42))
println(get(fruits,"mango",42))

# キーたちを取得したい時
println(keys(fruits))

# valueたちを取得したいとき
println(values(fruits))

# 複数の辞書もマージできるよ
a = Dict("apple"=>1,"banana"=>2)
b = Dict("banana"=>1,"peach"=>2)

# 後から来た方は上書きされる

println(merge(a,b))
println(merge(b,a))


text:
using DrWatson, Test

# Test commit function
com = gitdescribe(@__DIR__)
@test com isa String

com = gitdescribe(dirname(@__DIR__))
@test com !== nothing
@test typeof(com) <: String
@test com[1] == 'v' # test that it has a version tag
# tag!
d1 = Dict(:x => 3, :y => 4)
d2 = Dict("x" => 3, "y" => 4)
for d in (d1, d2)
    d = tag!(d, gitpath=@__DIR__)

    @test haskey(d, keytype(d)(:gitcommit))
    @test d[keytype(d)(:gitcommit)] |> typeof <: String
end

# @tag!
for d in (d1, d2)
    d = @tag!(d, gitpath=@__DIR__)
    @test haskey(d, keytype(d)(:gitcommit))
    @test d[keytype(d)(:gitcommit)] |> typeof <: String
    @test split(d[keytype(d)(:script)], '#')[1] == basename(@__FILE__)
end

# Tag kw-functions

d = Dict(:x => 3, :y => 4)
d_new = tag!(d)

@test d == d_new

d_new = tag!(d,gitpath=@__DIR__,source="foo")
@test endswith(d_new[:script],"foo")

d_new = @tag!(d,gitpath=@__DIR__)
@test split(d_new[keytype(d_new)(:script)], '#')[1] == basename(@__FILE__)

ex = @macroexpand @tag!(d,gitpath="path")

@test ex.args[1].name == :tag!
@test ex.args[2] == :d
@test ex.args[3].head == :kw
@test ex.args[3].args[1] == :gitpath
@test ex.args[3].args[2] == "path"
@test ex.args[4].head == :kw

# Test force kw

d = Dict(:x => 3, :y => 4, :gitcommit => "")
@test tag!(d,gitpath=@__DIR__)[:gitcommit] == ""
@test tag!(d,gitpath=@__DIR__,force=true)[:gitcommit] == com
d = Dict(:x => 3, :y => 4, :gitcommit => "")
@test (@tag!(d, gitpath=@__DIR__,force=true))[:gitcommit] == com

# Test dictionary expansion
c = Dict(:a => [1, 2], :b => 4);
c1 = [ Dict(:a=>1,:b=>4)
 Dict(:a=>2,:b=>4)]

v1 = dict_list(c)
for el in c1
    @test el ∈ v1
end
@test keytype(eltype(v1)) == Symbol

c[:c] = "test"; c[:d] = ["lala", "lulu"];
c2 = [ Dict(:a=>1,:b=>4,:d=>"lala",:c=>"test")
 Dict(:a=>2,:b=>4,:d=>"lala",:c=>"test")
 Dict(:a=>1,:b=>4,:d=>"lulu",:c=>"test")
 Dict(:a=>2,:b=>4,:d=>"lulu",:c=>"test")]

v2 = dict_list(c)
for el in c2
    @test el ∈ v2
end
@test keytype(eltype(v2)) == Symbol

c[:e] = [[1, 2], [3, 5]];
c3 = [
Dict(:a=>1,:b=>4,:d=>"lala",:e=>[1, 2],:c=>"test")
Dict(:a=>2,:b=>4,:d=>"lala",:e=>[1, 2],:c=>"test")
Dict(:a=>1,:b=>4,:d=>"lulu",:e=>[1, 2],:c=>"test")
Dict(:a=>2,:b=>4,:d=>"lulu",:e=>[1, 2],:c=>"test")
Dict(:a=>1,:b=>4,:d=>"lala",:e=>[3, 5],:c=>"test")
Dict(:a=>2,:b=>4,:d=>"lala",:e=>[3, 5],:c=>"test")
Dict(:a=>1,:b=>4,:d=>"lulu",:e=>[3, 5],:c=>"test")
Dict(:a=>2,:b=>4,:d=>"lulu",:e=>[3, 5],:c=>"test")
]

v3 = dict_list(c)
@test dict_list_count(c) == length(c3) == length(v3)
for el in c3
    @test el ∈ v3
end
@test keytype(eltype(v3)) == Symbol

v4 = dict_list(Dict(:a => 1, :b => 2.0)) # both non-iterable
@test keytype(eltype(v4)) == Symbol

v5 = dict_list(Dict(:a => [1], :b => 2.0)) # one non-iterable
@test keytype(eltype(v5)) == Symbol

v6 = dict_list(Dict(:a => [1], :b => [2.0])) # both iterable
@test keytype(eltype(v6)) == Symbol

struct SolverA end
struct SolverB end

# Simple example with no chained dependencies

p = Dict(:α => 1,
         :solver => [SolverA,SolverB],
         :c => @onlyif(:solver == SolverA , [100,200]),
         :d => @onlyif(:solver == SolverB, 1)
        )

@test Set([ Dict(:α => 1, :solver => SolverA, :c => 100),
           Dict(:α => 1, :solver => SolverA, :c => 200),
           Dict(:α => 1, :solver => SolverB, :d => 1),
          ]) == Set(dict_list(p))

# Advanced example with chained dependency. SolverA => :c => :d
p = Dict(:α => 1,
         :solver => [SolverA,SolverB],
         :c => @onlyif(:solver == SolverA , [100,200]),
         :d => @onlyif(:c == 100, 1)
        )

@test Set([
           Dict(:α => 1, :solver => SolverA, :c => 100, :d => 1),
           Dict(:α => 1, :solver => SolverA, :c => 200),
           Dict(:α => 1, :solver => SolverB),
          ]) == Set(dict_list(p))

# Advanced condition definition

test_param = @onlyif(begin
                d = Dict( :f => (conds...)->all(conds) )
                cond1 = :b == :c
                cond2 = :d == :something
                cond3 = "d" == :d
                cond4 = :α^2 == 1
                return d[:f](cond1, cond2, cond3, cond4)
            end, nothing)

dummy_dict = Dict(:α=>1, :b=>1, :c=>1, :d=>:something, "d"=>:something)

@test test_param.condition(dummy_dict, Dict(:α=>1, :b=>1, :c=>1, :d=>:something, "d"=>:something))
@test !test_param.condition(dummy_dict, Dict(:α=>2, :b=>1, :c=>1, :d=>:something, "d"=>:something))
@test !test_param.condition(dummy_dict, Dict(:α=>1, :b=>2, :c=>1, :d=>:something, "d"=>:something))
@test !test_param.condition(dummy_dict, Dict(:α=>1, :b=>1, :c=>1, :d=>:foo, "d"=>:something))
@test !test_param.condition(dummy_dict, Dict(:α=>1, :b=>1, :c=>1, :d=>:something, "d"=>:foo))
@test !test_param.condition(dummy_dict, Dict(:b=>1, :c=>1, :d=>:something, "d"=>:something))

module TestMod
    struct Foo end
end

dummy_dict = Dict(:solver => TestMod.Foo)
test_param = @onlyif(:solver == TestMod.Foo,[100,200])

@test test_param[1].condition(dummy_dict,Dict(:solver=>TestMod.Foo))
@test !test_param[1].condition(dummy_dict,Dict(:solver=>:Foo))

# partially restricted and mixed keytypes parameters

p = Dict(
         :a => :a1,
         "b" => [:b1,:b2],
         :c => [:c1,@onlyif("b" == :b2, :c2)],
        )

@test Set(dict_list(p)) == Set([
                          Dict(:a=>:a1, "b"=>:b1, :c=>:c1),
                          Dict(:a=>:a1, "b"=>:b2, :c=>:c1),
                          Dict(:a=>:a1, "b"=>:b2, :c=>:c2),
                         ])

# Every value restriced, but always at least 1 value available

p = Dict(
         :a => [1,2],
         :b => @onlyif(:a==1,[1,2]),
         :c => @onlyif(:a==3,[1,2]),
        )

@test Set(dict_list(p)) == Set([
                                Dict(:a => 2)
                                Dict(:a => 1,:b => 1)
                                Dict(:a => 1,:b => 2)
                               ])

p = Dict(
         :a => [1,2],
         :b => @onlyif(:a==1,[1,2]),
         :c => [@onlyif(:b==1,1), @onlyif(:b==2,2)],
        )

@test Set(dict_list(p)) == Set([
                                Dict(:a => 1,:b => 1,:c => 1)
                                Dict(:a => 1,:b => 2,:c => 2)
                                Dict(:a => 2)
                               ])

p = Dict(
         :a => [1,2],
         :b => @onlyif(:a==1,[1,2]),
         :c => [@onlyif(:b==1,1), @onlyif(:b==1,2)],
        )

@test Set(dict_list(p)) == Set([
                                Dict(:a => 1,:b => 1,:c => 1)
                                Dict(:a => 1,:b => 1,:c => 2)
                                Dict(:a => 2)
                                Dict(:a => 1,:b => 2)
                               ])

@test Set(dict_list(Dict(
   :a => [1,2,3],
   :b => [@onlyif(:a==1, 10), @onlyif(:a==2, [20]), @onlyif(:a==3, [30,30])]))) == Set(
           [Dict(:a => 1,:b => 10),
            Dict{Symbol,Any}(:a => 2,:b => [20]),
            Dict{Symbol,Any}(:a => 3,:b => [30, 30])])

# Testing nested @onlyif calls
@test Set(dict_list(Dict(
                   :a=>[1,2],
                   :b => [3,4],
                   :c => @onlyif( :a == 2, [5, @onlyif(:b == 4, 6)])
                  ))) == Set([Dict(:a => 1,:b => 3),
                              Dict(:a => 2,:b => 3,:c => 5),
                              Dict(:a => 1,:b => 4),
                              Dict(:a => 2,:b => 4,:c => 5),
                              Dict(:a => 2,:b => 4,:c => 6)])

# Test dict_list retaining original types

dlist = dict_list(Dict(
       :n => [1],
       :h => [10, 15, 20],
       :a => 0.01,:d=>0.001))

@test typeof(dlist[1][:a]) == Float64
@test typeof(dlist[1][:n]) == Int
@test typeof(dlist[1][:h]) == Int
@test typeof(dlist[1][:d]) == Float64

### tmpsave ###
tmpdir = joinpath(@__DIR__, "tmp")
ret = tmpsave(v3, tmpdir)
for r in ret
    @test isfile(joinpath(tmpdir, r))
    a = load(joinpath(tmpdir, r))
    @test a ∈ v3
end
rm(tmpdir, force = true, recursive = true)
@test !isdir(tmpdir)


text:
module Lale

using PyCall

# load abstract super-types and utils
using AMLPipelineBase
using AMLPipelineBase.AbsTypes
export fit, fit!, transform, transform!, fit_transform!

using AMLPipelineBase
using AMLPipelineBase: AbsTypes, Utils, BaselineModels, Pipelines
using AMLPipelineBase: BaseFilters, FeatureSelectors, DecisionTreeLearners
using AMLPipelineBase: EnsembleMethods, CrossValidators
using AMLPipelineBase: NARemovers

export Machine, Learner, Transformer, Workflow, Computer
export holdout, kfold, score, infer_eltype, nested_dict_to_tuples,
       nested_dict_set!, nested_dict_merge, create_transformer,
       mergedict, getiris, getprofb,
       skipmean,skipmedian,skipstd,
       aggregatorclskipmissing,
       find_catnum_columns,
       train_test_split

export Baseline, Identity
export Imputer,OneHotEncoder,Wrapper
export PrunedTree,RandomForest,Adaboost
export VoteEnsemble, StackEnsemble, BestLearner
export FeatureSelector, CatFeatureSelector, NumFeatureSelector, CatNumDiscriminator
export crossvalidate
export NARemover
export @pipeline, @pipelinex
export Pipeline, ComboPipeline

import AMLPipelineBase.AbsTypes: fit!, transform!

# ------
module LaleAbsTypes
   using ..AbsTypes
   using DataFrames: DataFrame
   using PyCall
   using PyCall: PyIterator, inspect

   export LaleOperator, fit, transform, predict, lalepropertynames
   export pretty_print

   abstract type LaleOperator <: Learner end

   #fit(o::Machine, x::DataFrame, y::Vector=Vector()) = fit!(o,x,y)
   #transform(o::Machine, x::DataFrame) = transform!(o,x)
   predict(o::Machine, x::DataFrame) = transform!(o,x)
   pretty_print(o::LaleOperator,args...) = nothing

   lalepropertynames(o::PyObject) = ispynull(o) ? Symbol[] : map(x->Symbol(first(x)), PyIterator{PyObject}(pycall(inspect."getmembers", PyObject, o)))
end

using .LaleAbsTypes

include("lalelibop.jl")
using .LaleLibOps
export LalePipeOptimizer, lalepipeoptimizers
export LalePipe, visualize, pretty_print
export >>, +, |, |>, &

include("laleop.jl")
using .LaleOps
export LaleOp, skops, autogenops, lalelibops, lalesummary

export laleoperator
export fit!, transform!, fit, transform, predict

function laleoperator(name::String,type::String="sklearn"; args...)
   try
      obj = LaleOp(name,type;args...)
      return obj
   catch ArgumentError
      sk = keys(LaleOps.sk_dict)
      ag = keys(LaleOps.ag_dict)
      ll = keys(LaleOps.ll_dict)
      mb = keys(LaleOps.mb_dict)
      println("Please choose among these pipeline elements:")
      println()
      println("sklearn: ", [sk...])
      println()
      println("autogen: ", [ag...])
      println()
      println("lale: ",    [ll...])
      println()
      println("imblearn: ",[mb...])
      println()
      throw(ArgumentError("$name does not exist"))
   end
end

end