text:
@testset "Constructors and basic properties" begin
    # numbers
    for T in (Float64, ComplexF64), V in (BinningAnalysis.Variance, BinningAnalysis.FastVariance)
        B = LogBinner(T; accumulator=V)
        B2 = LogBinner(T; accumulator=V)

        @test length(B) == 0
        @test ndims(B) == 0
        @test isempty(B)
        @test eltype(B) == T
        @test capacity(B) == 2^32 - 1
        @test BinningAnalysis.nlevels(B) == 32
        @test B == B2
        @test B ≈ B2
        @test !(B != B2)

        B3 = LogBinner(B, capacity=10_000)
        @test B3 == B
        @test B3 ≈ B
        @test capacity(B3) == 16383

        buffer = rand(rng, T, 1000)
        append!(B, buffer)
        @test length(B) == 1000
        @test !isempty(B)
        @test !(B == B2)
        @test !(B ≈ B2)
        @test B != B2

        @test B3 != B
        @test isempty(B3)
        B3 = LogBinner(B, capacity=1000)
        @test B3 == B
        @test B3 ≈ B
        @test !isempty(B3)

        append!(B2, buffer)
        @test B == B2
        @test B ≈ B2
        @test !(B != B2)

        append!(B, rand(rng, T, 24))
        @test_throws OverflowError LogBinner(B, capacity=1000)

        empty!(B)
        @test length(B) == 0
        @test isempty(B)
    end

    # arrays
    for T in (Float64, ComplexF64), V in (BinningAnalysis.Variance, BinningAnalysis.FastVariance)
        B = LogBinner(zeros(T, 2, 3), accumulator=V)
        B2 = LogBinner(zeros(T, 2, 3), accumulator=V)

        @test length(B) == 0
        @test ndims(B) == 2
        @test isempty(B)
        @test eltype(B) == Array{T, 2}
        @test capacity(B) == 2^32 - 1
        @test BinningAnalysis.nlevels(B) == 32
        @test B == B2
        @test B ≈ B2
        @test !(B != B2)

        B3 = LogBinner(B, capacity=10_000)
        @test B3 == B
        @test B3 ≈ B
        @test capacity(B3) == 16383

        buffer = [rand(rng, T, 2, 3) for _ in 1:1000]
        append!(B, buffer)
        @test length(B) == 1000
        @test !isempty(B)
        @test !(B == B2)
        @test !(B ≈ B2)
        @test B != B2

        @test B3 != B
        @test isempty(B3)
        B3 = LogBinner(B, capacity=1000)
        @test B3 == B
        @test B3 ≈ B
        @test !isempty(B3)

        append!(B2, buffer)
        @test B == B2
        @test B ≈ B2
        @test !(B != B2)

        append!(B, [rand(rng, T, 2, 3) for _ in 1:24])
        @test_throws OverflowError LogBinner(B, capacity=1000)

        empty!(B)
        @test length(B) == 0
        @test isempty(B)
    end

    # Constructor arguments
    B = LogBinner(capacity=12345)
    @test capacity(B) == 16383
    @test_throws ArgumentError LogBinner(capacity=0)
    @test_throws ArgumentError LogBinner(capacity=-1)

    # Test error on overflow
    B = LogBinner(capacity=1)
    push!(B, 1.0)
    @test_throws OverflowError push!(B, 2.0)

    B = LogBinner(zeros(2,2), capacity=1)
    push!(B, rand(rng, 2,2))
    @test_throws OverflowError push!(B, rand(rng, 2,2))

    # time series constructor (#26)
    x = rand(rng, 10)
    B = LogBinner(x)
    @test length(B) == 10
    @test mean(B) ≈ mean(x)

    # Test equality of different sizes
    B2 = LogBinner()
    @test !(B == B2)
    @test !(B ≈ B2)
    @test !(B2 == B)
    @test !(B2 ≈ B)
    @test B != B2
    @test B2 != B
    append!(B2, x)
    @test B == B2
    @test B ≈ B2
    @test B2 == B
    @test B2 ≈ B
    @test !(B != B2)
    @test !(B2 != B)

    # Again for Array inputs
    x = [rand(rng, 2,3) for _ in 1:5]
    B = LogBinner(x)
    @test length(B) == 5
    @test mean(B) ≈ mean(x)

    B2 = LogBinner(zeros(2,3))
    @test !(B == B2)
    @test !(B ≈ B2)
    @test !(B2 == B)
    @test !(B2 ≈ B)
    @test B != B2
    @test B2 != B
    append!(B2, x)
    @test B == B2
    @test B ≈ B2
    @test B2 == B
    @test B2 ≈ B
    @test !(B != B2)
    @test !(B2 != B)
end

@testset "Check variance for complex values" begin
    # NOTE
    # Due to the different (mathematically equivalent) versions of the variance
    # calculated here, the values are only approximately the same. (Float error)
    StableRNGs.seed!(rng, 123)
    xs = rand(rng, ComplexF64, 1_000_000)
    BA = LogBinner(ComplexF64)

    # Test small set (off by one errors are large here)
    for x in xs[1:10]; push!(BA, x) end
    @test var(BA, 1) ≈ var(xs[1:10])
    @test varN(BA, 1) ≈ var(xs[1:10])/10

    # Test full set
    for x in xs[11:end]; push!(BA, x) end
    @test var(BA, 1) ≈ var(xs)
    @test varN(BA, 1) ≈ var(xs)/1_000_000

    # all_* methods
    @test isapprox(all_vars(BA), [0.1665847540639013, 0.08321852969184897, 0.0416386170603128, 0.02080121032765401, 0.010404232627786522, 0.005208087488833146, 0.0025893712817839105, 0.0012979310682925164, 0.0006498645790020864, 0.0003304988207069092, 0.00016135937775941, 8.244734583734034e-5, 4.1904380844737295e-5, 2.1161425558227823e-5, 1.0185807653474429e-5, 5.794287254586299e-6, 2.4120659277772666e-6, 2.132350340200331e-6, 2.6572093711370505e-7])
    @test isapprox(all_varNs(BA), zero(all_varNs(BA)), atol=1e-6)
    @test isapprox(all_taus(BA), [0.0, -0.0004433019126909299, -9.090214413765008e-5, -0.0005254725849703767, -0.0003512686979507129, 0.00022224596475073355, -0.0025962520837867764, -0.001317479480883732, -0.0006266293029093539, 0.007928314422850802, -0.0037744232961079427, 0.0070976636615537325, 0.01547113876712125, 0.020618795006264512, 0.001187681531748308, 0.07971364057757968, -0.01734960354334969, 0.4143137944503138, -0.23414840330888914])
    @test isapprox(all_std_errors(BA), [0.00040814795609423463, 0.00040796698320292777, 0.0004081108528834429, 0.00040793342915386584, 0.00040800456130365054, 0.0004082386552528566, 0.00040708692196405705, 0.00040760987484761305, 0.00040789211844220833, 0.00041137115440979214, 0.0004066045146218508, 0.00041103464473798035, 0.00041441436480546865, 0.0004164784602525399, 0.0004086324183611065, 0.00043948027087255004, 0.00040100423337559803, 0.0005519252459993797, 0.0002976132485367013])

    @test isapprox(tau(BA), 0.001187681531748308)
    @test isapprox(std_error(BA), 0.0004086324183611065)
end

@testset "Check variance for complex vectors" begin
    StableRNGs.seed!(rng, 123)
    xs = [rand(rng, ComplexF64, 3) for _ in 1:1_000_000]
    BA = LogBinner(zeros(ComplexF64, 3))

    # Test small set (off by one errors are large here)
    for x in xs[1:10]; push!(BA, x) end
    @test var(BA, 1) ≈ var(xs[1:10])
    @test varN(BA, 1) ≈ var(xs[1:10])/10

    # Test full set
    for x in xs[11:end]; push!(BA, x) end
    @test var(BA, 1) ≈ var(xs)
    @test varN(BA, 1) ≈ var(xs)/1_000_000

    # all_std_errors for <:AbstractArray
    @test all(isapprox.(all_std_errors(BA), Ref(zeros(3)), atol=1e-2))

    @test all(isapprox.(tau(BA), [0.025099907279786615, 0.07567404578645565, 0.021586917849341525], atol=1e-6))
    @test all(isapprox.(std_error(BA), [0.000418454645620472, 0.0004382193124499972, 0.00041696433722483093], atol=1e-6))
end

@testset "Checking converging data (Complex)" begin
    BA = LogBinner(ComplexF64)

    # block of maximally correlated values:
    N_corr = 16

    # number of blocks
    N_blocks = 131_072 # 2^17

    for _ in 1:N_blocks
        x = rand(rng, ComplexF64)
        for __ in 1:N_corr
            push!(BA, x)
        end
    end

    # BA must diverge until 16 values are binned
    for lvl in 1:4
        @test !has_converged(BA, lvl)
    end

    # Afterwards it must converge (for a while)
    for lvl in 5:8
        @test has_converged(BA, lvl)
    end
    # Later values may fluctuate due to small samples / small threshold in
    # has_converged

    # means should be consistent
    means = BinningAnalysis.all_means(BA)
    for x in means
        @test x ≈ means[1]
    end
end

@testset "Checking converging data (Vector)" begin
    BA = LogBinner(zeros(3))

    # block of maximally correlated values:
    N_corr = 16

    # number of blocks
    N_blocks = 131_072 # 2^17

    for _ in 1:N_blocks
        x = rand(rng, Float64, 3)
        for __ in 1:N_corr
            push!(BA, x)
        end
    end

    # BA must diverge until 16 values are binned
    for lvl in 1:4
        @test !has_converged(BA, lvl)
    end

    # Afterwards it must converge (for a while)
    for lvl in 5:8
        @test has_converged(BA, lvl)
    end
    # Later values may fluctuate due to small samples / small threshold in
    # has_converged

    # means should be consistent
    means = BinningAnalysis.all_means(BA)
    for i in eachindex(means)
        @test means[i] ≈ means[1]
    end
end

@testset "Type promotion" begin
    Bf = LogBinner(zero(1.)) # Float64 LogBinner
    Bc = LogBinner(zero(im)) # Float64 LogBinner

    # Check that this doesn't throw (TODO: is there a better way?)
    @test (append!(Bf, rand(rng, 1:10, 10000)); true)
    @test (append!(Bc, rand(rng, 10000)); true)
end

@testset "Sum-type heuristic" begin
    # numbers
    @test typeof(LogBinner(zero(Int64))) == LogBinner{Float64, 32, BinningAnalysis.Variance{Float64}}
    @test typeof(LogBinner(zero(ComplexF16))) == LogBinner{ComplexF64, 32, BinningAnalysis.Variance{ComplexF64}}

    # arrays
    @test typeof(LogBinner(zeros(Int64, 2,2))) == LogBinner{Matrix{Float64}, 32, BinningAnalysis.Variance{Matrix{Float64}}}
    @test typeof(LogBinner(zeros(ComplexF16, 2,2))) == LogBinner{Matrix{ComplexF64}, 32, BinningAnalysis.Variance{Matrix{ComplexF64}}}
end

@testset "Indexing Bounds" begin
    BA = LogBinner(zero(Float64); capacity=1)
    for func in [:var, :varN, :mean, :tau]
        # check if func(BA, 0) throws BoundsError
        # It should as level 1 is now the initial level
        @test_throws BoundsError @eval $func($BA, 0)
        # Check that level 1 exists
        @test (@eval $func($BA, 1); true)
        # Check that level 2 throws BoundsError
        @test_throws BoundsError @eval $func($BA, 2)
    end
end

@testset "_reliable_level" begin
    BA = LogBinner()
    # Empty Binner
    @test BinningAnalysis._reliable_level(BA) == 1
    @test isnan(std_error(BA, BinningAnalysis._reliable_level(BA)))

    # One Element should still return NaN (due to 1/(n-1))
    push!(BA, rand(rng, ))
    @test BinningAnalysis._reliable_level(BA) == 1
    @test isnan(std_error(BA, BinningAnalysis._reliable_level(BA)))

    # Two elements should return some value
    push!(BA, rand(rng, ))
    @test BinningAnalysis._reliable_level(BA) == 1
    @test !isnan(std_error(BA, BinningAnalysis._reliable_level(BA)))

    # same behavior up to (including) 63 values (31 binned in first binned lvl)
    append!(BA, rand(rng, 61))
    @test BinningAnalysis._reliable_level(BA) == 1
    @test !isnan(std_error(BA, BinningAnalysis._reliable_level(BA)))

    # at 64 or more values, the lvl should be increasing
    push!(BA, rand(rng, ))
    @test BinningAnalysis._reliable_level(BA) == 2
    @test !isnan(std_error(BA, BinningAnalysis._reliable_level(BA)))
end

@testset "Cosmetics (show, print, etc.)" begin
    B = LogBinner();
    # empty binner
    oldstdout = stdout
    (read_pipe, write_pipe) = redirect_stdout()
    println(B) # compact
    show(write_pipe, MIME"text/plain"(), B) # full
    redirect_stdout(oldstdout);
    close(write_pipe);

    # compact
    @test readline(read_pipe) == "LogBinner{Float64,32}()"
    # full
    @test readline(read_pipe) == "LogBinner{Float64,32}"
    @test readline(read_pipe) == "| Count: 0"
    @test length(readlines(read_pipe)) == 0
    close(read_pipe);

    # filled binner
    StableRNGs.seed!(rng, 123)
    append!(B, rand(rng, 1000))
    (read_pipe, write_pipe) = redirect_stdout()
    show(write_pipe, MIME"text/plain"(), B)
    redirect_stdout(oldstdout);
    close(write_pipe);
    @test readline(read_pipe) == "LogBinner{Float64,32}"
    @test readline(read_pipe) == "| Count: 1000"
    @test readline(read_pipe) == "| Mean: 0.49387"
    @test readline(read_pipe) == "| StdError: 0.00952"
    @test length(readlines(read_pipe)) == 0
    close(read_pipe);
end

max_stars_repo_name: crstnbr/BinningAnalysis

max_issues_repo_name: crstnbr/BinningAnalysis

max_forks_repo_name: crstnbr/BinningAnalysis


text:
# This file was generated, do not modify it. # hide
yMode = [mode(ŷ[i]) for i in 1:length(ŷ)]
y = coerce(y[:,1], OrderedFactor)
yMode = coerce(yMode, OrderedFactor)
confusion_matrix(yMode, y)

ext: jl


text:
mutable struct Action
  name::AbstractString
  value::Float64
end

generate_arbitrary_actions(num_actions::Integer) =[Action(string(i), i) for i in 1:num_actions]


text:
formulations = ["objscale10-4_BLOCKS_cholesky_ref0", "objscale10-4_BLOCKS_MD_ref0",
"objscale10-4_BLOCKS_cholesky_real_ref0", "objscale10-4_BLOCKS_cholesky_fusion_Dan_nb_max_10%_ref0", "objscale10-4_BLOCKS_cholesky_fusion50_x6_ref0",
"objscale10-4_BLOCKS_cholesky_fusion50_x5_ref0","objscale10-4_BLOCKS_cholesky_fusion50_x4_ref0", "objscale10-4_BLOCKS_cholesky_fusion50_x3_ref0",
"objscale10-4_BLOCKS_cholesky_fusion50_x2_ref0","objscale10-4_BLOCKS_cholesky_fusion50_x1_ref0"]

function read_mosek_log_more_info(log_file)
    lines = readlines(log_file)
    i = 1
    while split(lines[i], ':')[1]!= "nb_coupling_constraints primal " && split(lines[i], ':')[1]!="nb_coupling_constraints "
        i+=1
    end
    nlc = parse(Int, split(lines[i], ':')[2])
    while split(lines[i], ':')[1] != "  Constraints            "
        i+=1
    end
    n_total_constraints = parse(Int,split(lines[i], ':')[2])
    m = n_total_constraints - nlc
    while split(lines[i], ':')[1] != "Optimizer terminated. Time" && i < length(lines)
        i+=1
    end
    time = split(lines[i], ':')[2][2:end]
    nb_iter = split(lines[i-1], ' ')[1]
    objective = lines[i+2][18:end]
    return time, nb_iter, objective, m, nlc
end

nb_formulations = length(formulations)

DATA_PATH = joinpath("..","data","decompositions")

recap = open("data_linear_regression\\recap_linear_regression_objscale10-4.csv", "w")
recap_pareto = open("data_linear_regression\\recap_error_pareto_objscale10-4.csv", "w")
write(recap, "Instance ; alpha ; beta ; constant ; error ; with LOG10 ; Instance ; alpha ; beta ; constant ; error \n")
write(recap_pareto, "alpha * log10( sum |Ci|^3) + (1-alpha) * log10((m+nlc)^3) \n
Instance ; 0.1 ; 0.2 ; 0.3 ; 0.4 ; 0.5 ; 0.6 ; 0.7 ; 0.8 ; 0.9 \n")

for instance in instances
    println(instance)
    i = 0
    t = zeros(nb_formulations,1)
    q1 = zeros(nb_formulations,1)
    q2 = zeros(nb_formulations,1)
    e = ones(nb_formulations,1)
    errors_pareto = zeros(9,1)
    f = open("data\\parameters_linear_regression\\data_objscale10-4_$instance.csv", "w")
    write(f, "FORMULATION ; m ; nlc ; sum |Ci|^3 ; time per iteration  ;  ; log10( sum |Ci|^3) ;  log10((m+nlc)^3) ; ; alpha 0.1 ; alpha 0.2 ; alpha 0.3\n")
    for formulation in formulations
        # println(formulation)
        i += 1
        log_file = joinpath(pwd(), "..", "data", "logs", repo, "$(instance)_$(formulation).log")
        time = nb_iter = objective =  m =  nlc = "#"
        time_per_iteration = "#"
        remove_row = false
        try
            time, nb_iter, objective, m, nlc = read_mosek_log_more_info(log_file)
            time_per_iteration = parse(Float64,time)/parse(Int64,nb_iter)
            t[i,1] = time_per_iteration
            q2[i,1] = (m+nlc)^3
        catch
            println("WARNING : pb read log")
            remove_row = true
            time = nb_iter = objective =  m =  nlc = 0
            time_per_iteration = 0
            # t = t[1:end-1, 1]
            # q2 = q2[1:end-1, 1]
            # q1 = q1[1:end-1, 1]
            # e = e[1:end-1, 1]
        end
        #sum |C_i|
        #open blocks file
        if formulation[1:8] == "objscale"
            blocks_formulation = formulation[14:end]
        else
            blocks_formulation = formulation
        end
        if formulation[end-3:end-1] == "ref"
            blocks_formulation = blocks_formulation[1:end-5]
        end

        if blocks_formulation[1:end-3]=="BLOCKS_cholesky_fusion50"
            blocks_formulation = "blocks_cholesky_fusion_size_max_50_$(blocks_formulation[end-1:end])"
        end
        BLOCKS = readdlm(joinpath(DATA_PATH, blocks_formulation ,instance*"_sdp_blocks.txt"))
        block_var = create_list_of_variables_per_block(BLOCKS)
        sum_size_cliques_cube = 0
        for (block, list_var) in block_var
            size_block = length(list_var)
            sum_size_cliques_cube += size_block^3
        end
        if remove_row == false
            q1[i,1] = sum_size_cliques_cube
        else
            t=t[setdiff(1:end, i), :]
            q1=q1[setdiff(1:end, i), :]
            q2=q2[setdiff(1:end, i), :]
            e=e[setdiff(1:end, i), :]
            i-=1
        end

        write(f, "$formulation ; $m ; $nlc ; $sum_size_cliques_cube ; $time_per_iteration  ; ;  $(log10(sum_size_cliques_cube));  $(log10((m+nlc)^3)); ;")
         for i in 1:9
             alpha = 0.1*i
             criterion = alpha * log10(sum_size_cliques_cube) + (1 - alpha)*log10((m+nlc)^3)
             write(f, " $(criterion/time_per_iteration) ;")
             errors_pareto[i] += 0.5*(time_per_iteration-criterion)^2
        end
        write(f, "\n")
    end
    A = [q1'*q1 q1'*q2 q1'*e ; q2'*q1 q2'*q2 q2'*e ; e'*q1 e'*q2 e'*e ]
    b = [ q1'*t ; q2'*t ; e'*t]
    x = A\b
    error =  0.5*(t-(x[1]*q1+x[2]*q2+x[3]*e))'*(t-(x[1]*q1+x[2]*q2+x[3]*e))
    println("alpha : $(x[1]) ; beta : $(x[2]) ; constant : $(x[3]) \n")
    println("error 0,5 ||t-(alpha q1 + beta q2 + c)||^2 : ", error)
    write(f, "\n alpha ; beta ; constant \n $(x[1]) ; $(x[2]) ; $(x[3]) \n")
    write(f, "error 0,5 ||t-(alpha q1 + beta q2 + c)||^2 :  ; $(error[1]) \n")
    write(recap, "$instance ; $(x[1]) ; $(x[2]) ; $(x[3]) ; $(error[1]) ; ")
    write(f, "\n WITH LOG10 \n")
    q1 = log10.(q1)
    q2 = log10.(q2)
    A = [q1'*q1 q1'*q2 q1'*e ; q2'*q1 q2'*q2 q2'*e ; e'*q1 e'*q2 e'*e ]
    b = [ q1'*t ; q2'*t ; e'*t]
    x = A\b
    error =  0.5*(t-(x[1]*q1+x[2]*q2+x[3]*e))'*(t-(x[1]*q1+x[2]*q2+x[3]*e))
    println("alpha : $(x[1]) ; beta : $(x[2]) ; constant : $(x[3]) \n")
    println("error 0,5 ||t-(alpha log10(q1) + beta log10(q2) + c)||^2 : ", error)
    write(f, "\n alpha ; beta ; constant \n $(x[1]) ; $(x[2]) ; $(x[3]) \n")
    write(f, "error 0,5 ||t-(alpha log10(q1) + beta log10(q2) + c)||^2 :  ; $(error[1]) \n")
    close(f)
    write(recap, " WITHLOG10 ; $instance ; $(x[1]) ; $(x[2]) ; $(x[3]) ; $(error[1]) \n ")
    write(recap_pareto, "$instance ; ")
    for i in 1:9
        write(recap_pareto, "$(errors_pareto[i]) ; ")
    end
    write(recap_pareto, "\n")
end

close(recap)
close(recap_pareto)

max_stars_repo_name: JulieSliwak/CliqueMerging

max_issues_repo_name: JulieSliwak/CliqueMerging

max_forks_repo_name: JulieSliwak/CliqueMerging


text:
using LinearAlgebra
using StaticArrays

include("./io.jl")

"""
    repulsion(npoints, dims, fileopts...; niter = 1e7, convergence_tol = 1e-10)

Generates a uniform grid of `npoints` points in `dims` dimensions, distributed uniformly on the surface of the unit sphere
using the REPULSION algorithm

The function iterates, up to a maximum of `niter` times, until the algorithm converges. The convergence criterion
`convergence_tol` is interpreted as the target for the maximum displacement between points in two successive iterations. The points
inherit their floating point precision from `convergence_tol`.
The points are returned as a vector of vectors of length `npoints`, with each element containing the angles `[α, β, γ]` and the
weight corresponding to each point (currently uniform weights are assigned)

The function optionally writes the resulting grid of points to a file  by calling `writedata`, which uses the contents of the optional
`fileopts` parameters to determine the file format and filename. The first fileopt parameter is a String containing the file format,
and the second is a String containing the filename. If the filename is not provided, a default is chosen using the current system time.
The file formats currently supported are: "CSV", "TSV", "MAT", "JLD", "JLD2". The [MAT](https://github.com/JuliaIO/MAT.jl) and
[JLD](https://github.com/JuliaIO/JLD.jl)/[JLD2](https://github.com/JuliaIO/JLD2.jl) file types require packages not in the standard library
that must be installed and loaded by the user.

This function is adapted from its counterpart in Spinach (https://spindynamics.org) for MATLAB. The Spinach
function, repulsion.m, is Copyright Ilya Kuprov and Frederic Mentink-Vigier
"""
function repulsion(npoints, dims, fileopts...; niter = 1e7, convergence_tol::T = 1e-10) where {T <: AbstractFloat}
    R_raw = rand(T, dims, npoints) .- one(T)/2 # generating the first set of random points
    R = [SVector{dims, T}(ntuple(i -> R_raw[n+i], Val(dims))) for n in 0:dims:(dims*npoints-1)]
    R_n = copy(R) # pre-allocate
    F = zeros(SVector{dims, T},npoints) # pre-allocate "forces"

    for i in 1:niter # starting the conjugate gradient optimization
        F .= zero(T) .* F # re-initializing F at each step
        for k in eachindex(R)
            F .-= dist_vector.(R,[R[k]]) # compute the distance between points, in place (speed-up)
        end
        R_n = normalize.(R_n - dims*F/npoints) # normalize
        max_diff = maximum(norm.(R-R_n)) # get the maximum distance
        R = R_n
	    if mod(i, 1000) == 0
	        @debug "Maximum difference " max_diff "At iteration " i # returns the convergence value every 1000 steps
	    end
        if max_diff <= convergence_tol
            @info "Converged at iteration " i # stops and return number of steps when converged
            break
        end
    end
    angles = @. vcat(toangles(R), one(T)/npoints) # computing the angles and weights

    if !isempty(fileopts)
        writedata(angles, fileopts...) # Optionally, write angles to file
    end

    return angles
end

"""
    dist_vector(a::SVector{N,T}, b::SVector{N,T}) where N where T

Computes the normalised distance between two vectors `a` and `b`, scaled by the cosine of the angle between them
"""
function dist_vector(a::SVector{N,T}, b::SVector{N,T}) where N where T
    dd = normalize(a - b)
    if isnan(first(dd))
        zeros(typeof(a))
    else
        dd*dot(a, b)
    end
end

"""
    toangles(p)

Converts a point `p` on the unit sphere to a vector of angles `[α, β, γ]`

The interpretation of `p` and the conversion depend on the number of dimensions:
- Two dimensions: `p` is treated as 2D cartesian co-ordinates, and the polar angle is returned as `β`
- Three dimensions: `p` is treated as 3D cartesian co-ordinates, and the spherical angles `θ` and `φ` are returned as `β` and `γ`
- Four dimensions: `p` is treated as a quaternion, and the ZYZ Euler angles are returned
All angles are shifted to lie in the ranges ``α ∈ [0,2π)``, ``β ∈ [0,π]``, ``γ ∈ [0,2π)``. Missing angles are returned as 0s
"""
function toangles(p::SVector{N,T}) where N where T
    angles = zeros(T, 3)
    if N == 2
	    angles[2] = atan(p[2], p[1])
    elseif N == 3
	    angles[3] = atan(p[2], p[1])
	    angles[2] = pi/(2*one(T)) + atan(p[3], hypot(p[1], p[2]))
    elseif N == 4
	    angles = quat_to_euler(p)
    end
    angles
end

"""
    quat_to_euler(q)

Converts a normalised quaternion to its corresponding Euler angles `[α, β, γ]` about the
ZYZ rotation axes.

An active rotation convention is assumed - invert the input quaternion
for passive rotations. This code is adapted from EasySpin's quat2euler implementation at:
https://github.com/StollLab/EasySpin/blob/main/easyspin/quat2euler.m
EasySpin is Copyright Stefan Stoll and other contributors
"""
function quat_to_euler(q::AbstractVector{T}) where T
    if length(q) != 4
        e = DomainError(q, "A Quaternion must have four elements")
        throw(e)
    end

    sy = 2*sqrt((q[3]*q[4] + q[2]*q[1])^2 + (q[2]*q[4] - q[3]*q[1])^2)
    beta = real(atan(sy, one(T) - 2*q[2]^2 - 2*q[3]^2))
    if sy < 1e-10 # Singularity check: 1e-10 ≈ 0 for normalised quantities
        alpha = zero(T)
        gamma = atan(-2*(q[4]*q[1] - q[2]*q[3]), one(T) - 2*(q[2]^2 - q[4]^2))
    else
        alpha = atan(2*(q[3]*q[4] - q[1]*q[2]), 2*(q[2]*q[4] + q[1]*q[3]))
        gamma = atan(2*(q[3]*q[4] + q[1]*q[2]), -2*(q[2]*q[4] - q[1]*q[3]))
        if alpha < 0
            alpha = alpha + pi + pi # Maintaining floating point precision
        end
        if gamma < 0
            gamma = gamma + pi + pi
        end
    end
    [alpha; beta; gamma]
end


text:
using MAT
fid = matopen("../data/nstest_in.mat")
b = read(fid, "b")
N = read(fid, "N")
alt = read(fid, "alt")
r = read(fid, "r")
eta = read(fid, "eta")

fidX = matopen("../data/X_in.mat")
X_in = read(fidX,"X")

#include("../src/AltMinFrob.jl")
using DistributedArrays
@everywhere include("../src/dclrMN_sync.jl")

@time ownerLMN,ownerRMN = dclrMN(b+N, eta, convert(Int,r), convert(Int,alt))

# Gather L and R results
L = fetch(@spawnat ownerL Main.L)
R = fetch(@spawnat ownerR Main.R)

@everywhere include("../src/dclr.jl")
@time ownerL,ownerR = dclr(b, convert(Int,r), convert(Int,alt))
LMN = fetch(@spawnat ownerLMN Main.L)
RMN = fetch(@spawnat ownerRMN Main.R)


text:
const ReComp = Union{Real,Complex}

struct Dual{T<:ReComp} <: Number
    value::T
    epsilon::T
end
Dual(x::S, y::T) where {S<:ReComp,T<:ReComp} = Dual(promote(x,y)...)
Dual(x::ReComp) = Dual(x, zero(x))
Dual{T}(x::ReComp) where T<:ReComp = Dual{T}(T(x), zero(T))

const ɛ = Dual(false, true)
const imɛ = Dual(Complex(false, false), Complex(false, true))

const Dual128 = Dual{Float64}
const Dual64  = Dual{Float32}
const Dual32  = Dual{Float16}
const DualComplex256 = Dual{ComplexF64}
const DualComplex128 = Dual{ComplexF32}
const DualComplex64  = Dual{ComplexF16}

Base.convert(::Type{Dual{T}}, z::Dual{T}) where {T<:ReComp} = z
Base.convert(::Type{Dual{T}}, z::Dual) where {T<:ReComp} = Dual{T}(convert(T, value(z)), convert(T, epsilon(z)))
Base.convert(::Type{Dual{T}}, x::Number) where {T<:ReComp} = Dual{T}(convert(T, x), convert(T, 0))
Base.convert(::Type{T}, z::Dual) where {T<:ReComp} = (epsilon(z)==0 ? convert(T, value(z)) : throw(InexactError()))

Base.promote_rule(::Type{Dual{T}}, ::Type{Dual{S}}) where {T<:ReComp,S<:ReComp} = Dual{promote_type(T, S)}
Base.promote_rule(::Type{Dual{T}}, ::Type{S}) where {T<:ReComp,S<:ReComp} = Dual{promote_type(T, S)}
Base.promote_rule(::Type{Dual{T}}, ::Type{T}) where {T<:ReComp} = Dual{T}

Base.widen(::Type{Dual{T}}) where {T} = Dual{widen(T)}

value(z::Dual) = z.value
epsilon(z::Dual) = z.epsilon

value(x::Number) = x
epsilon(x::Number) = zero(typeof(x))

dual(x::ReComp, y::ReComp) = Dual(x, y)
dual(x::ReComp) = Dual(x)
dual(z::Dual) = z

const realpart = value
const dualpart = epsilon

Base.isnan(z::Dual) = isnan(value(z))
Base.isinf(z::Dual) = isinf(value(z))
Base.isfinite(z::Dual) = isfinite(value(z))
isdual(x::Dual) = true
isdual(x::Number) = false
Base.eps(z::Dual) = eps(value(z))
Base.eps(::Type{Dual{T}}) where {T} = eps(T)

function dual_show(io::IO, z::Dual{T}, compact::Bool) where T<:Real
    x, y = value(z), epsilon(z)
    if isnan(x) || isfinite(y)
        compact ? show(IOContext(io, :compact=>true), x) : show(io, x)
        if signbit(y)==1 && !isnan(y)
            y = -y
            print(io, compact ? "-" : " - ")
        else
            print(io, compact ? "+" : " + ")
        end
        compact ? show(IOContext(io, :compact=>true), y) : show(io, y)
        printtimes(io, y)
        print(io, "ɛ")
    else
        print(io, "Dual{",T,"}(", x, ",", y, ")")
    end
end

function dual_show(io::IO, z::Dual{T}, compact::Bool) where T<:Complex
    x, y = value(z), epsilon(z)
    xr, xi = reim(x)
    yr, yi = reim(y)
    if isnan(x) || isfinite(y)
        compact ? show(IOContext(io, :compact=>true), x) : show(io, x)
        if signbit(yr)==1 && !isnan(y)
            yr = -yr
            print(io, " - ")
        else
            print(io, " + ")
        end
        if compact
            if signbit(yi)==1 && !isnan(y)
                yi = -yi
                show(IOContext(io, :compact=>true), yr)
                printtimes(io, yr)
                print(io, "ɛ-")
                show(IOContext(io, :compact=>true), yi)
            else
                show(IOContext(io, :compact=>true), yr)
                print(io, "ɛ+")
                show(IOContext(io, :compact=>true), yi)
            end
        else
            if signbit(yi)==1 && !isnan(y)
                yi = -yi
                show(io, yr)
                printtimes(io, yr)
                print(io, "ɛ - ")
                show(io, yi)
            else
                show(io, yr)
                print(io, "ɛ + ")
                show(io, yi)
            end
        end
        printtimes(io, yi)
        print(io, "imɛ")
    else
        print(io, "Dual{",T,"}(", x, ",", y, ")")
    end
end

function dual_show(io::IO, z::Dual{T}, compact::Bool) where T<:Bool
    x, y = value(z), epsilon(z)
    if !value(z) && epsilon(z)
        print(io, "ɛ")
    else
        print(io, "Dual{",T,"}(", x, ",", y, ")")
    end
end

function dual_show(io::IO, z::Dual{Complex{T}}, compact::Bool) where T<:Bool
    x, y = value(z), epsilon(z)
    xr, xi = reim(x)
    yr, yi = reim(y)
    if !xr
        if xi*!yr*!yi
            print(io, "im")
        elseif !xi*yr*!yi
            print(io, "ɛ")
        elseif !xi*!yr*yi
            print(io, "imɛ")
        end
    else
        print(io, "Dual{",T,"}(", x, ",", y, ")")
    end
end

function printtimes(io::IO, x::Real)
    if !(isa(x,Integer) || isa(x,Rational) ||
         isa(x,AbstractFloat) && isfinite(x))
        print(io, "*")
    end
end

Base.show(io::IO, z::Dual) = dual_show(io, z, get(IOContext(io), :compact, false))

function Base.read(s::IO, ::Type{Dual{T}}) where T<:ReComp
    x = read(s, T)
    y = read(s, T)
    Dual{T}(x, y)
end
function Base.write(s::IO, z::Dual)
    write(s, value(z))
    write(s, epsilon(z))
end

## Generic functions of dual numbers ##

Base.convert(::Type{Dual}, z::Dual) = z
Base.convert(::Type{Dual}, x::Number) = Dual(x)

Base.:(==)(z::Dual, w::Dual) = value(z) == value(w)
Base.:(==)(z::Dual, x::Number) = value(z) == x
Base.:(==)(x::Number, z::Dual) = value(z) == x

Base.isequal(z::Dual, w::Dual) = isequal(value(z),value(w)) && isequal(epsilon(z), epsilon(w))
Base.isequal(z::Dual, x::Number) = isequal(value(z), x) && isequal(epsilon(z), zero(x))
Base.isequal(x::Number, z::Dual) = isequal(z, x)

Base.isless(z::Dual{<:Real},w::Dual{<:Real}) = value(z) < value(w)
Base.isless(z::Real,w::Dual{<:Real}) = z < value(w)
Base.isless(z::Dual{<:Real},w::Real) = value(z) < w

Base.hash(z::Dual) = (x = hash(value(z)); epsilon(z)==0 ? x : bitmix(x,hash(epsilon(z))))

Base.float(z::Union{Dual{T}, Dual{Complex{T}}}) where {T<:AbstractFloat} = z
Base.complex(z::Dual{<:Complex}) = z

Base.floor(z::Dual) = floor(value(z))
Base.ceil(z::Dual)  = ceil(value(z))
Base.trunc(z::Dual) = trunc(value(z))
Base.round(z::Dual) = round(value(z))
Base.floor(::Type{T}, z::Dual) where {T<:Real} = floor(T, value(z))
Base.ceil( ::Type{T}, z::Dual) where {T<:Real} = ceil( T, value(z))
Base.trunc(::Type{T}, z::Dual) where {T<:Real} = trunc(T, value(z))
Base.round(::Type{T}, z::Dual) where {T<:Real} = round(T, value(z))

for op in (:real, :imag, :conj, :float, :complex)
    @eval Base.$op(z::Dual) = Dual($op(value(z)), $op(epsilon(z)))
end

Base.abs(z::Dual) = sqrt(abs2(z))
Base.abs2(z::Dual) = real(conj(z)*z)

Base.real(z::Dual{<:Real}) = z
Base.abs(z::Dual{<:Real}) = z ≥ 0 ? z : -z

Base.angle(z::Dual{<:Real}) = z ≥ 0 ? zero(z) : one(z)*π
function Base.angle(z::Dual{Complex{T}}) where T<:Real
    if z == 0
        if imag(epsilon(z)) == 0
            Dual(zero(T), zero(T))
        else
            Dual(zero(T), convert(T, Inf))
        end
    else
        real(log(sign(z)) / im)
    end
end

Base.flipsign(x::Dual,y::Dual) = y == 0 ? flipsign(x, epsilon(y)) : flipsign(x, value(y))
Base.flipsign(x, y::Dual) = y == 0 ? flipsign(x, epsilon(y)) : flipsign(x, value(y))
Base.flipsign(x::Dual, y) = dual(flipsign(value(x), y), flipsign(epsilon(x), y))

# algebraic definitions
conjdual(z::Dual) = Dual(value(z),-epsilon(z))
absdual(z::Dual) = abs(value(z))
abs2dual(z::Dual) = abs2(value(z))

# algebra

Base.:+(z::Dual, w::Dual) = Dual(value(z)+value(w), epsilon(z)+epsilon(w))
Base.:+(z::Number, w::Dual) = Dual(z+value(w), epsilon(w))
Base.:+(z::Dual, w::Number) = Dual(value(z)+w, epsilon(z))

Base.:-(z::Dual) = Dual(-value(z), -epsilon(z))
Base.:-(z::Dual, w::Dual) = Dual(value(z)-value(w), epsilon(z)-epsilon(w))
Base.:-(z::Number, w::Dual) = Dual(z-value(w), -epsilon(w))
Base.:-(z::Dual, w::Number) = Dual(value(z)-w, epsilon(z))

# avoid ambiguous definition with Bool*Number
Base.:*(x::Bool, z::Dual) = ifelse(x, z, ifelse(signbit(real(value(z)))==0, zero(z), -zero(z)))
Base.:*(x::Dual, z::Bool) = z*x

Base.:*(z::Dual, w::Dual) = Dual(value(z)*value(w), epsilon(z)*value(w)+value(z)*epsilon(w))
Base.:*(x::Number, z::Dual) = Dual(x*value(z), x*epsilon(z))
Base.:*(z::Dual, x::Number) = Dual(x*value(z), x*epsilon(z))

Base.:/(z::Dual, w::Dual) = Dual(value(z)/value(w), (epsilon(z)*value(w)-value(z)*epsilon(w))/(value(w)*value(w)))
Base.:/(z::Number, w::Dual) = Dual(z/value(w), -z*epsilon(w)/value(w)^2)
Base.:/(z::Dual, x::Number) = Dual(value(z)/x, epsilon(z)/x)

for f in [:(Base.:^), :(NaNMath.pow)]
    @eval function ($f)(z::Dual{T1}, w::Dual{T2}) where {T1, T2}
        T = promote_type(T1, T2) # for type stability in ? : statements
        val = $f(value(z), value(w))

        ezvw = epsilon(z) * value(w) # for using in ? : statement
        du1 = iszero(ezvw) ? zero(T) : ezvw * $f(value(z), value(w) - 1)
        ew = epsilon(w) # for using in ? : statement
        # the float is for type stability because log promotes to floats
        du2 = iszero(ew) ? zero(float(T)) : ew * val * log(value(z))
        du = du1 + du2

        Dual(val, du)
    end
end

Base.mod(z::Dual, n::Number) = Dual(mod(value(z), n), epsilon(z))

# introduce a boolean !iszero(n) for hard zero behaviour to combat NaNs
function pow(z::Dual, n::AbstractFloat)
    return Dual(value(z)^n, !iszero(n) * (epsilon(z) * n * value(z)^(n - 1)))
end
function pow(z::Dual{T}, n::Integer) where T
    iszero(n) && return Dual(one(T), zero(T)) # avoid DomainError Int^(negative Int)
    isone(z) && return Dual(one(T), epsilon(z) * n)
    return Dual(value(z)^n, epsilon(z) * n * value(z)^(n - 1))
end
# these first two definitions are needed to fix ambiguity warnings
for T1 ∈ (:Integer, :Rational, :Number)
    @eval Base.:^(z::Dual{T}, n::$T1) where T = pow(z, n)
end

NaNMath.pow(z::Dual{T}, n::Number) where T = Dual(NaNMath.pow(value(z),n), epsilon(z)*n*NaNMath.pow(value(z),n-1))
NaNMath.pow(z::Number, w::Dual{T}) where T = Dual(NaNMath.pow(z,value(w)), epsilon(w)*NaNMath.pow(z,value(w))*log(z))

Base.inv(z::Dual) = dual(inv(value(z)),-epsilon(z)/value(z)^2)

# force use of NaNMath functions in derivative calculations
function to_nanmath(x::Expr)
    if x.head == :call
        funsym = Expr(:.,:NaNMath,Base.Meta.quot(x.args[1]))
        return Expr(:call,funsym,[to_nanmath(z) for z in x.args[2:end]]...)
    else
        return Expr(:call,[to_nanmath(z) for z in x.args]...)
    end
end
to_nanmath(x) = x

for (funsym, exp) in Calculus.symbolic_derivatives_1arg()
    funsym == :exp && continue
    funsym == :abs2 && continue
    funsym == :inv && continue
    if isdefined(SpecialFunctions, funsym)
        @eval function SpecialFunctions.$(funsym)(z::Dual)
            x = value(z)
            xp = epsilon(z)
            Dual($(funsym)(x),xp*$exp)
        end
    elseif isdefined(Base, funsym)
        @eval function Base.$(funsym)(z::Dual)
            x = value(z)
            xp = epsilon(z)
            Dual($(funsym)(x),xp*$exp)
        end
    end
    # extend corresponding NaNMath methods
    if funsym in (:sin, :cos, :tan, :asin, :acos, :acosh, :atanh, :log, :log2, :log10,
          :lgamma, :log1p)
        funsym = Expr(:.,:NaNMath,Base.Meta.quot(funsym))
        @eval function $(funsym)(z::Dual)
            x = value(z)
            xp = epsilon(z)
            Dual($(funsym)(x),xp*$(to_nanmath(exp)))
        end
    end
end

# only need to compute exp/cis once
Base.exp(z::Dual) = (expval = exp(value(z)); Dual(expval, epsilon(z)*expval))
Base.cis(z::Dual) = (cisval = cis(value(z)); Dual(cisval, im*epsilon(z)*cisval))

Base.exp10(x::Dual) = (y = exp10(value(x)); Dual(y, y * log(10) * epsilon(x)))

## TODO: should be generated in Calculus
Base.sinpi(z::Dual) = Dual(sinpi(value(z)),epsilon(z)*cospi(value(z))*π)
Base.cospi(z::Dual) = Dual(cospi(value(z)),-epsilon(z)*sinpi(value(z))*π)

Base.checkindex(::Type{Bool}, inds::AbstractUnitRange, i::Dual) = checkindex(Bool, inds, value(i))


text:
# Transforming a graph By adding and removing nodes and edges.

using UUIDs

export transform!, applyRule!, transformingGraph!

function transform!(graph, addarcs, removearcs)
    (e -> remove_edge!(graph, e)).(removearcs)
    (e -> add_edge!(graph, e)).(addarcs)
    graph
end

"""
    applyRule!(graph, rule)
Apply rule to each node of graph and modify the graph's edges
as determined by rule.  Rule is a function that takes the graph
and a node as arguments and returns the set of edges to add and
the set of edges to remove.
"""
function applyRule!(graph, rule)
    nodes(graph) .|>
	(node -> rule(graph, node)) .|>
	(add_remove -> transform!(graph, add_remove...))
    nothing
end

struct NonlocalTransfer <: Exception
    uid

    function NonlocalTransfer()
        new(UUIDs.uuid1())
    end
end

function Base.showerror(io::IO, e::NonlocalTransfer)
    print(io, "NonlocalTransfer not caught!")
end

function transformingGraph!(f)
    addset = Set()
    removeset = Set()
    exittag = NonlocalTransfer()
    add(arc) = push!(addset, arc)
    remove(arc) = push!(removeset, arc)
    function check(condition)
	if !condition
	    throw(exittag)
	end
    end
    try
	f(check, add, remove)
    catch e
	if e != exittag
	    rethrow(e)
	end
    end
    return addset, removeset
end


text:
using UnicodePlots

out(plt) = println(IOContext(stdout, :color=>true), plt)

function print_mem_usage()
	out = read(`top -bn1 -p $(getpid())`, String)
	res = split(split(out,  "\n")[end-1])
	println("RES: ", res[6], "   %MEM: ", res[10], "   %CPU: ", res[9])
	parse(Float64, replace(res[10], "," => "."))
end

module MyFunctions
function f1(x::Matrix{Float64})
    out = copy(x)
    out[x .< 0.5] .= 0.0
    out[x .>= 0.5] .= 1.0
    out
end
function f2(x::Matrix{Float64})
	reshape([x[i] < 0.5 ? 0.0 : 1.0 for i in eachindex(x)], size(x))
end
end

struct Foo
	functions::Vector{Function}
end

function Foo(function_module::Module, fname::String)
	functions = Array{Function}(undef, 1)
	functions[1] = getfield(function_module, Symbol(fname))
	Foo(functions)
end

mutable struct Bar
	f::Function
end

function Bar(foo::Foo)
	f = foo.functions[1]
	Bar(f)
end

function test(bar::Bar)
	mem_usage = Vector{Float64}()
	for i in 1:100
		x = rand(1000, 1000)
		bar.f(x)
		mem = print_mem_usage()
		push!(mem_usage, mem)
		out(lineplot(mem_usage, title = "%MEM"))
	end
end

function_module = MyFunctions
fname = "f1"
foo = Foo(function_module, fname)
bar = Bar(foo)

test(bar)


text:
using AbstractGPs
using TemporalGPs
using Distributions
using StatsFuns: logistic

# In this example we are showing how to work with non-Gaussian likelihoods,
# here a binary classification problem with a logistic link,
# by using variable augmentations as in [Galy-Fajou et al, 2020](https://arxiv.org/abs/2002.114510).
# We base our example on the "exact_time_inference" example

# Load up the separable kernel from TemporalGPs.
using TemporalGPs: RegularSpacing

# Build a GP as per usual, and wrap it inside a TemporalGPs.jl object.
f_raw = GP(Matern52Kernel());
f = to_sde(f_raw, SArrayStorage(Float64));

# Specify a collection of inputs. Must be increasing.
T = 1_000;
x = RegularSpacing(0.0, 1e-1, T);

# Generate some synthetic data from the GP and get a random binary output
σ²_noise = 0.01;
f_true = rand(f(x, σ²_noise));
y = rand.(Bernoulli.(logistic.(f_true)));
y_sign = sign.(y .- 0.5)

# We are using the augmentation trick and look for the optimal
# expected value of the variable ω

η₁_like(ω, g, β) = g .+ β .* ω # The augmented likelihood first natural parameter
Λ_like(ω, γ) = 2 * γ * ω # The augmented likelihood precision (diagonal)

# These are the values corresponding to our problem
γ = 0.5
g = 0.5 * y_sign
β = 0

function compute_optimal_expectation(f, x, g, β, γ; n_iter=5)
    T = length(x)
    ω̄ = rand(T) # preallocation
    c = zeros(T) # preallocation
    for i in 1:n_iter
        Λ = Λ_like(ω̄, γ)
        p_f = marginals(posterior(f(x, inv.(Λ)), inv.(Λ) .* η₁_like(ω̄, g, β))(x))
        @. c = sqrt(var(p_f) + abs2(mean(p_f))) / 2
        @. ω̄ = 0.5 * tanh(c) / c
    end
    return ω̄
end

ω̄ = compute_optimal_expectation(f, x, g, β, γ)
Λ = Λ_like(ω̄, γ)
f_post = posterior(f(x, inv.(Λ)), inv.(Λ) .* η₁_like(ω̄, g, β))

# Specify some locations at which to make predictions.
T_pr = 1_200_000;
x_pr = RegularSpacing(0.0, 1e-4, T_pr);

# Compute the exact posterior marginals at `x_pr`.
f_post_marginals = marginals(f_post(x_pr));
m_post_marginals = mean.(f_post_marginals);
σ_post_marginals = std.(f_post_marginals);

# Generate a few posterior samples. Not fantastically-well optimised at present.
f_post_samples = [rand(f_post(x_pr)) for _ in 1:5];

# Visualise the posterior. The if block is just to prevent it running in CI.
if get(ENV, "TESTING", "FALSE") == "FALSE"
    using Plots
    plt = plot()
    plot!(plt, f_post(x_pr); ribbon_scale=3.0, label="");
    plot!(plt, x_pr, f_post_samples; color=:red, alpha=0.3, label="");
    plot!(plt, x, f_true; label="", lw=2.0, color=:blue); # Plot the true latent GP on top
    scatter!(plt, x, y; label="", markersize=1.0, alpha=1.0); # Plot the data
    savefig(plt, "posterior.png");
end


text:
"""
    get_optimized_1tree(cost)

Obtain an optimized 1tree given the cost matrix. This calls [`get_1tree`](@ref) and then adds extra
benefits for each point which has only one edge to make it more likely in the next run of [`get_1tree`](@ref)
that the point will have two edges. A negative benefit is used for points which have more than two edges.

Return the tree as a list of edges and the cost of the tree.
"""
function get_optimized_1tree(original_cost; runs=10)
    cost = deepcopy(original_cost)
    N = size(cost, 1)
    tree, lb = get_1tree(cost)
    degrees = zeros(Int, N)
    point_benefit = zeros(Float64, N)
    extra_cost = 0.0
    non_inf_cost = filter(c->!isinf(c), cost)
    max_cost = maximum(non_inf_cost)

    cost_factor = max_cost/N
    for _ in 1:runs
        degrees .= 0
        point_benefit .= 0.0
        for edge in tree
            degrees[edge.src] += 1
            degrees[edge.dst] += 1
        end
        for i in 1:N
            point_benefit[i] = cost_factor*(2-degrees[i])
            cost[i, :] .-= point_benefit[i]
            cost[:, i] .-= point_benefit[i]
        end
        extra_cost = 2*sum(point_benefit)
        tree, lb = get_1tree(cost)
        cost_factor *= 0.9
    end

    return tree, lb+extra_cost
end

"""
    get_1tree(cost)

Get a 1tree and cost of it based on a given cost matrix.
Create an MST first for the first 1 to N-1 points and then add two edges which connects the Nth point to its two nearest neighbors.

Can be optimized by adding extra costs with [`get_optimized_1tree`](@ref) which obtains a better lower bound.

Return the tree as a list of edges and the cost which can be used as the lower bound for the TSP problem.
"""
function get_1tree(cost)
    N = size(cost, 1)

    tree = prim_mst(complete_graph(N-1), cost)

    mst_cost = get_edges_cost(cost, tree)

    # extra cost for the two edges from the last point to the two nearest neighbors
    nearest_neighbors = partialsortperm(cost[N,:], 1:3)
    n_actual_neighbors = 0
    extra_costs = 0.0
    for neighbor in nearest_neighbors
        n_actual_neighbors == 2 && break
        if neighbor != N
            n_actual_neighbors += 1
            push!(tree, Edge(N, neighbor))
            extra_costs += cost[N, neighbor]
        end
    end

    lb = mst_cost + extra_costs

    return tree, lb
end


text:
"Return an array of the unique angles for three vectors `a`, `b`, and `c`."
function angles_abc(a, b, c)
    return ang(a, b), ang(a, c), ang(b, c)
end # function

"""
Score distortion of cell angles from 90 degrees as a sum of squares.

Higher scores indicate larger deviations from 90 degrees.
"""
function score_angles_cell(cell)
    ab, ac, bc = angles_abc(cell[1,:], cell[2,:], cell[3,:])
    return (ab-π/2)^2 + (ac-π/2)^2 + (bc-π/2)^2
end

"""
Score the lengths of the unique vectors `a`, `b`, and `c` pairwise as sum of square differences.
"""
function score_lengths_abc(a, b, c)
    ab = norm(a) - norm(b)
    ac = norm(a) - norm(c)
    bc = norm(b) - norm(c)
    return (ab^2 + ac^2 + bc^2)
end # function

"""
Score the lengths of a cell's vectors pairwise as sum of square differences.
"""
score_lengths_cell(cell) = score_lengths_abc(cell[1,:], cell[2,:], cell[3,:])

const M1TO1 = [-1, 0, 1]
"""
Score based on the maximum radius of a sphere needed to contain the whole cell.

The smallest sphere (lowest score) gives the most compact cell. The radius of
the sphere needed is given by the maximum distance of any cell corner to the
center of the cell.

As given by equation (3) by Militzer.
"""
function score_sphere(cell)
    score = -Inf
    for i=M1TO1, j=M1TO1, k=M1TO1
        candidate = 0.5*norm(i*cell[1,:]+j*cell[2,:]+k*cell[3,:])
        score = max(score, candidate)
    end # for
    return score
end #function

const PREV_NAN = MVector{3}(Array{Float64}(undef, 3))
"""
Redefine the lattice vectors of a cell to point to the next nearest image
See equations (5)-(7) by Militzer.
"""
function closestimagevectors(cell)
    a, b, c = cell[1,:], cell[2,:], cell[3,:]

    prev_b = PREV_NAN
    prev_c = PREV_NAN
    while b != prev_b && c != prev_c
        a, b, c = sort([a, b, c], by=norm)
        prev_b[:] = b
        prev_c[:] = c

        b = b - a*round(dot(b,a)/dot(a,a), RoundNearestTiesUp)
        c = c - a*round(dot(c,a)/dot(a,a), RoundNearestTiesUp)
        c = c - b*round(dot(c,b)/dot(b,b), RoundNearestTiesUp)
    end # while

    return a, b, c
end # function

"""
Score based on the minimum distance to the next image.

The cell with the largest minimum image distance (maximizing peroidic distances)
is optimal. Solves equation (4) given by Militzer.
"""
function score_image_distance(cell)
    a, b, c = closestimagevectors(cell)
    score = Inf
    for i=M1TO1, j=M1TO1, k=M1TO1
        if (i^2 + j^2 + k^2) > 0
            candidate = norm(i*a + j*b + k*c)
            score = min(score, candidate)
        end # if
    end # for
    return score
end # function

"""
The score is the negative of the minimum distance to the next image, so that the largest distance is the minimum.
"""
score_image_distance_min(cell) = -score_image_distance(cell)

"""
Compare the scores for two differenct cells by the scoring function.

Returns true if an optimal cell has been assigned, false if no optimal cell is found (scores are tied).

Modifies the score arrays and the optimal cell and matrix in place.
"""
function compare_cell_scores!(scorefunc, optcell, candcell, optscore, candscore, optmatrix, candmatrix, idx)
    optsc = isnan(optscore[idx]) ? scorefunc(optcell) : optscore[idx]
    optscore[idx] = optsc
    candsc = scorefunc(candcell)
    candscore[idx] = candsc
    if optsc != candsc
        # we can choose the cell with the minimum sphere size
        if argmin((optsc, candsc)) == 1
            # optimum cell stays the same
            return true
        else
            # candidate cell is better
            optcell[:,:] = candcell
            optmatrix[:,:] = candmatrix
            optscore[:] = candscore[:]
            return true
        end # if
    else
        return false
    end # if
end # function

const SCOREFUNCS = [score_sphere, score_image_distance_min, score_angles_cell, score_lengths_cell]
const NSCOREFUNCS = length(SCOREFUNCS)
const CANDSCORES = Array{Float64}(undef, NSCOREFUNCS)
"""
Loop through scoring criteria functions to find the optimal cell in a pairwise comparison.

Generally speaking, scoring should be as lazy as possible and only progress
through the scoring criteria as necessary.

This method modifies the optimal cell and score with the optimal values found.

"""
function cellchoose!(optcell, optscore, candcell, best_matrix, cur_matrix)
    candscore = CANDSCORES
    for i=1:NSCOREFUNCS
        candscore[i] = NaN
    end
    for i=1:NSCOREFUNCS
        optimum_found = compare_cell_scores!(SCOREFUNCS[i], optcell, candcell, optscore, candscore, best_matrix, cur_matrix, i)
        if optimum_found
            return
        end # if
    end #for
    # these cells are equivalent, just return (keeping the original best cell)
    return
end # function


text:
using BioSimulator
include("winData/multiomics_networks_simulation/julia_functions.jl")

function res2df(output :: BioSimulator.SimData)
    t, data = get_data(output)
    id2ind  = output.id2ind

    d1, d2, d3 = size(data)

    ids = collect(keys(id2ind))

    df = DataFrame()
    df[:time] = repeat(collect(t), outer=[d3])
    for i = ids
        df[i] = vec(convert(Array, view(data, id2ind[i], :, :)))
    end
    df[:trial] = repeat(collect(1:d3), inner=[d2])
    return df
end

model = BioSimulator.Network("coucou")

model <= BioSimulator.Species("X", 25)
model <= BioSimulator.Reaction("transcription", 0.02, "0 --> X")
model <= BioSimulator.Reaction("decay", 0.00069, "X --> 0")

result = BioSimulator.simulate(model, algorithm = SSA, time = 1000.0, epochs = 1000)
resdf = res2df(result)

plot(resdf[:, :time], resdf[:, :X])

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

stochmodel = Dict("reactionsnames"=>Any["transcription1GCN1", "RNAdecayR1GCN1", "transcription1GCN2", "RNAdecayR1GCN2", "transcription1GCN3", "RNAdecayR1GCN3", "transcription1GCN4", "RNAdecayR1GCN4", "transcription2GCN1", "RNAdecayR2GCN1", "transcription2GCN2", "RNAdecayR2GCN2", "transcription2GCN3", "RNAdecayR2GCN3", "transcription2GCN4", "RNAdecayR2GCN4", "transcription3GCN1", "RNAdecayR3GCN1", "transcription3GCN2", "RNAdecayR3GCN2", "transcription3GCN3", "RNAdecayR3GCN3", "transcription3GCN4", "RNAdecayR3GCN4", "transcription4GCN1", "RNAdecayR4GCN1", "transcription4GCN2", "RNAdecayR4GCN2", "transcription4GCN3", "RNAdecayR4GCN3", "transcription4GCN4", "RNAdecayR4GCN4", "transcription5GCN1", "RNAdecayR5GCN1", "transcription5GCN2", "RNAdecayR5GCN2", "transcription5GCN3", "RNAdecayR5GCN3", "transcription5GCN4", "RNAdecayR5GCN4", "transcription6GCN1", "RNAdecayR6GCN1", "transcription6GCN2", "RNAdecayR6GCN2", "transcription6GCN3", "RNAdecayR6GCN3", "transcription6GCN4", "RNAdecayR6GCN4", "transcription7GCN1", "RNAdecayR7GCN1", "transcription7GCN2", "RNAdecayR7GCN2", "transcription7GCN3", "RNAdecayR7GCN3", "transcription7GCN4", "RNAdecayR7GCN4", "transcription8GCN1", "RNAdecayR8GCN1", "transcription8GCN2", "RNAdecayR8GCN2", "transcription8GCN3", "RNAdecayR8GCN3", "transcription8GCN4", "RNAdecayR8GCN4", "transcription9GCN1", "RNAdecayR9GCN1", "transcription9GCN2", "RNAdecayR9GCN2", "transcription9GCN3", "RNAdecayR9GCN3", "transcription9GCN4", "RNAdecayR9GCN4", "transcription10GCN1", "RNAdecayR10GCN1", "transcription10GCN2", "RNAdecayR10GCN2", "transcription10GCN3", "RNAdecayR10GCN3", "transcription10GCN4", "RNAdecayR10GCN4", "translation1GCN1", "proteindecayP1GCN1", "translation1GCN2", "proteindecayP1GCN2", "translation1GCN3", "proteindecayP1GCN3", "translation1GCN4", "proteindecayP1GCN4", "translation2GCN1", "proteindecayP2GCN1", "translation2GCN2", "proteindecayP2GCN2", "translation2GCN3", "proteindecayP2GCN3", "translation2GCN4", "proteindecayP2GCN4", "translation3GCN1", "proteindecayP3GCN1", "translation3GCN2", "proteindecayP3GCN2", "translation3GCN3", "proteindecayP3GCN3", "translation3GCN4", "proteindecayP3GCN4", "translation4GCN1", "proteindecayP4GCN1", "translation4GCN2", "proteindecayP4GCN2", "translation4GCN3", "proteindecayP4GCN3", "translation4GCN4", "proteindecayP4GCN4", "translation5GCN1", "proteindecayP5GCN1", "translation5GCN2", "proteindecayP5GCN2", "translation5GCN3", "proteindecayP5GCN3", "translation5GCN4", "proteindecayP5GCN4", "translation7GCN1", "proteindecayP7GCN1", "translation7GCN2", "proteindecayP7GCN2", "translation7GCN3", "proteindecayP7GCN3", "translation7GCN4", "proteindecayP7GCN4"],
    "propensities"=>Any[:(0.020032 * ((QTLeffects["GCN1"])["qtlTCrate"])[1]), :(0.000694444 * ((QTLeffects["GCN1"])["qtlRDrate"])[1]), :(0.020032 * ((QTLeffects["GCN2"])["qtlTCrate"])[1]), :(0.000694444 * ((QTLeffects["GCN2"])["qtlRDrate"])[1]), :(0.020032 * ((QTLeffects["GCN3"])["qtlTCrate"])[1]), :(0.000694444 * ((QTLeffects["GCN3"])["qtlRDrate"])[1]), :(0.020032 * ((QTLeffects["GCN4"])["qtlTCrate"])[1]), :(0.000694444 * ((QTLeffects["GCN4"])["qtlRDrate"])[1]), :(0.0435213 * ((QTLeffects["GCN1"])["qtlTCrate"])[2]), :(0.000534759 * ((QTLeffects["GCN1"])["qtlRDrate"])[2]), :(0.0435213 * ((QTLeffects["GCN2"])["qtlTCrate"])[2]), :(0.000534759 * ((QTLeffects["GCN2"])["qtlRDrate"])[2]), :(0.0435213 * ((QTLeffects["GCN3"])["qtlTCrate"])[2]), :(0.000534759 * ((QTLeffects["GCN3"])["qtlRDrate"])[2]), :(0.0435213 * ((QTLeffects["GCN4"])["qtlTCrate"])[2]), :(0.000534759 * ((QTLeffects["GCN4"])["qtlRDrate"])[2]), :(0.0766692 * ((QTLeffects["GCN1"])["qtlTCrate"])[3]), :(0.00310559 * ((QTLeffects["GCN1"])["qtlRDrate"])[3]), :(0.0766692 * ((QTLeffects["GCN2"])["qtlTCrate"])[3]), :(0.00310559 * ((QTLeffects["GCN2"])["qtlRDrate"])[3]), :(0.0766692 * ((QTLeffects["GCN3"])["qtlTCrate"])[3]), :(0.00310559 * ((QTLeffects["GCN3"])["qtlRDrate"])[3]), :(0.0766692 * ((QTLeffects["GCN4"])["qtlTCrate"])[3]), :(0.00310559 * ((QTLeffects["GCN4"])["qtlRDrate"])[3]), :(0.0897525 * ((QTLeffects["GCN1"])["qtlTCrate"])[4]), :(0.000286451 * ((QTLeffects["GCN1"])["qtlRDrate"])[4]), :(0.0897525 * ((QTLeffects["GCN2"])["qtlTCrate"])[4]), :(0.000286451 * ((QTLeffects["GCN2"])["qtlRDrate"])[4]), :(0.0897525 * ((QTLeffects["GCN3"])["qtlTCrate"])[4]), :(0.000286451 * ((QTLeffects["GCN3"])["qtlRDrate"])[4]), :(0.0897525 * ((QTLeffects["GCN4"])["qtlTCrate"])[4]), :(0.000286451 * ((QTLeffects["GCN4"])["qtlRDrate"])[4]), :(0.071785 * ((QTLeffects["GCN1"])["qtlTCrate"])[5]), :(0.00247525 * ((QTLeffects["GCN1"])["qtlRDrate"])[5]), :(0.071785 * ((QTLeffects["GCN2"])["qtlTCrate"])[5]), :(0.00247525 * ((QTLeffects["GCN2"])["qtlRDrate"])[5]), :(0.071785 * ((QTLeffects["GCN3"])["qtlTCrate"])[5]), :(0.00247525 * ((QTLeffects["GCN3"])["qtlRDrate"])[5]), :(0.071785 * ((QTLeffects["GCN4"])["qtlTCrate"])[5]), :(0.00247525 * ((QTLeffects["GCN4"])["qtlRDrate"])[5]), :(0.0514678 * ((QTLeffects["GCN1"])["qtlTCrate"])[6]), :(0.000312891 * ((QTLeffects["GCN1"])["qtlRDrate"])[6]), :(0.0514678 * ((QTLeffects["GCN2"])["qtlTCrate"])[6]), :(0.000312891 * ((QTLeffects["GCN2"])["qtlRDrate"])[6]), :(0.0514678 * ((QTLeffects["GCN3"])["qtlTCrate"])[6]), :(0.000312891 * ((QTLeffects["GCN3"])["qtlRDrate"])[6]), :(0.0514678 * ((QTLeffects["GCN4"])["qtlTCrate"])[6]), :(0.000312891 * ((QTLeffects["GCN4"])["qtlRDrate"])[6]), :(0.0220153 * ((QTLeffects["GCN1"])["qtlTCrate"])[7]), :(0.000307692 * ((QTLeffects["GCN1"])["qtlRDrate"])[7]), :(0.0220153 * ((QTLeffects["GCN2"])["qtlTCrate"])[7]), :(0.000307692 * ((QTLeffects["GCN2"])["qtlRDrate"])[7]), :(0.0220153 * ((QTLeffects["GCN3"])["qtlTCrate"])[7]), :(0.000307692 * ((QTLeffects["GCN3"])["qtlRDrate"])[7]), :(0.0220153 * ((QTLeffects["GCN4"])["qtlTCrate"])[7]), :(0.000307692 * ((QTLeffects["GCN4"])["qtlRDrate"])[7]), :(0.0721572 * ((QTLeffects["GCN1"])["qtlTCrate"])[8]), :(0.000376506 * ((QTLeffects["GCN1"])["qtlRDrate"])[8]), :(0.0721572 * ((QTLeffects["GCN2"])["qtlTCrate"])[8]), :(0.000376506 * ((QTLeffects["GCN2"])["qtlRDrate"])[8]), :(0.0721572 * ((QTLeffects["GCN3"])["qtlTCrate"])[8]), :(0.000376506 * ((QTLeffects["GCN3"])["qtlRDrate"])[8]), :(0.0721572 * ((QTLeffects["GCN4"])["qtlTCrate"])[8]), :(0.000376506 * ((QTLeffects["GCN4"])["qtlRDrate"])[8]), :(0.0658958 * ((QTLeffects["GCN1"])["qtlTCrate"])[9]), :(0.000824402 * ((QTLeffects["GCN1"])["qtlRDrate"])[9]), :(0.0658958 * ((QTLeffects["GCN2"])["qtlTCrate"])[9]), :(0.000824402 * ((QTLeffects["GCN2"])["qtlRDrate"])[9]), :(0.0658958 * ((QTLeffects["GCN3"])["qtlTCrate"])[9]), :(0.000824402 * ((QTLeffects["GCN3"])["qtlRDrate"])[9]), :(0.0658958 * ((QTLeffects["GCN4"])["qtlTCrate"])[9]), :(0.000824402 * ((QTLeffects["GCN4"])["qtlRDrate"])[9]), :(0.0930883 * ((QTLeffects["GCN1"])["qtlTCrate"])[10]), :(0.000316957 * ((QTLeffects["GCN1"])["qtlRDrate"])[10]), :(0.0930883 * ((QTLeffects["GCN2"])["qtlTCrate"])[10]), :(0.000316957 * ((QTLeffects["GCN2"])["qtlRDrate"])[10]), :(0.0930883 * ((QTLeffects["GCN3"])["qtlTCrate"])[10]), :(0.000316957 * ((QTLeffects["GCN3"])["qtlRDrate"])[10]), :(0.0930883 * ((QTLeffects["GCN4"])["qtlTCrate"])[10]), :(0.000316957 * ((QTLeffects["GCN4"])["qtlRDrate"])[10]), :(3.34576 * ((QTLeffects["GCN1"])["qtlTLrate"])[1]), :(7.62369e-5 * ((QTLeffects["GCN1"])["qtlPDrate"])[1]), :(3.34576 * ((QTLeffects["GCN2"])["qtlTLrate"])[1]), :(7.62369e-5 * ((QTLeffects["GCN2"])["qtlPDrate"])[1]), :(3.34576 * ((QTLeffects["GCN3"])["qtlTLrate"])[1]), :(7.62369e-5 * ((QTLeffects["GCN3"])["qtlPDrate"])[1]), :(3.34576 * ((QTLeffects["GCN4"])["qtlTLrate"])[1]), :(7.62369e-5 * ((QTLeffects["GCN4"])["qtlPDrate"])[1]), :(3.00529 * ((QTLeffects["GCN1"])["qtlTLrate"])[2]), :(0.00012364 * ((QTLeffects["GCN1"])["qtlPDrate"])[2]), :(3.00529 * ((QTLeffects["GCN2"])["qtlTLrate"])[2]), :(0.00012364 * ((QTLeffects["GCN2"])["qtlPDrate"])[2]), :(3.00529 * ((QTLeffects["GCN3"])["qtlTLrate"])[2]), :(0.00012364 * ((QTLeffects["GCN3"])["qtlPDrate"])[2]), :(3.00529 * ((QTLeffects["GCN4"])["qtlTLrate"])[2]), :(0.00012364 * ((QTLeffects["GCN4"])["qtlPDrate"])[2]), :(1.15585 * ((QTLeffects["GCN1"])["qtlTLrate"])[3]), :(0.00010269 * ((QTLeffects["GCN1"])["qtlPDrate"])[3]), :(1.15585 * ((QTLeffects["GCN2"])["qtlTLrate"])[3]), :(0.00010269 * ((QTLeffects["GCN2"])["qtlPDrate"])[3]), :(1.15585 * ((QTLeffects["GCN3"])["qtlTLrate"])[3]), :(0.00010269 * ((QTLeffects["GCN3"])["qtlPDrate"])[3]), :(1.15585 * ((QTLeffects["GCN4"])["qtlTLrate"])[3]), :(0.00010269 * ((QTLeffects["GCN4"])["qtlPDrate"])[3]), :(3.5602 * ((QTLeffects["GCN1"])["qtlTLrate"])[4]), :(0.000100725 * ((QTLeffects["GCN1"])["qtlPDrate"])[4]), :(3.5602 * ((QTLeffects["GCN2"])["qtlTLrate"])[4]), :(0.000100725 * ((QTLeffects["GCN2"])["qtlPDrate"])[4]), :(3.5602 * ((QTLeffects["GCN3"])["qtlTLrate"])[4]), :(0.000100725 * ((QTLeffects["GCN3"])["qtlPDrate"])[4]), :(3.5602 * ((QTLeffects["GCN4"])["qtlTLrate"])[4]), :(0.000100725 * ((QTLeffects["GCN4"])["qtlPDrate"])[4]), :(0.567232 * ((QTLeffects["GCN1"])["qtlTLrate"])[5]), :(9.17936e-5 * ((QTLeffects["GCN1"])["qtlPDrate"])[5]), :(0.567232 * ((QTLeffects["GCN2"])["qtlTLrate"])[5]), :(9.17936e-5 * ((QTLeffects["GCN2"])["qtlPDrate"])[5]), :(0.567232 * ((QTLeffects["GCN3"])["qtlTLrate"])[5]), :(9.17936e-5 * ((QTLeffects["GCN3"])["qtlPDrate"])[5]), :(0.567232 * ((QTLeffects["GCN4"])["qtlTLrate"])[5]), :(9.17936e-5 * ((QTLeffects["GCN4"])["qtlPDrate"])[5]), :(4.44623 * ((QTLeffects["GCN1"])["qtlTLrate"])[7]), :(0.00010015 * ((QTLeffects["GCN1"])["qtlPDrate"])[7]), :(4.44623 * ((QTLeffects["GCN2"])["qtlTLrate"])[7]), :(0.00010015 * ((QTLeffects["GCN2"])["qtlPDrate"])[7]), :(4.44623 * ((QTLeffects["GCN3"])["qtlTLrate"])[7]), :(0.00010015 * ((QTLeffects["GCN3"])["qtlPDrate"])[7]), :(4.44623 * ((QTLeffects["GCN4"])["qtlTLrate"])[7]), :(0.00010015 * ((QTLeffects["GCN4"])["qtlPDrate"])[7])],
    "reactions"=>Any["0 --> R1GCN1", "R1GCN1 --> 0", "0 --> R1GCN2", "R1GCN2 --> 0", "0 --> R1GCN3", "R1GCN3 --> 0", "0 --> R1GCN4", "R1GCN4 --> 0", "0 --> R2GCN1", "R2GCN1 --> 0", "0 --> R2GCN2", "R2GCN2 --> 0", "0 --> R2GCN3", "R2GCN3 --> 0", "0 --> R2GCN4", "R2GCN4 --> 0", "0 --> R3GCN1", "R3GCN1 --> 0", "0 --> R3GCN2", "R3GCN2 --> 0", "0 --> R3GCN3", "R3GCN3 --> 0", "0 --> R3GCN4", "R3GCN4 --> 0", "0 --> R4GCN1", "R4GCN1 --> 0", "0 --> R4GCN2", "R4GCN2 --> 0", "0 --> R4GCN3", "R4GCN3 --> 0", "0 --> R4GCN4", "R4GCN4 --> 0", "0 --> R5GCN1", "R5GCN1 --> 0", "0 --> R5GCN2", "R5GCN2 --> 0", "0 --> R5GCN3", "R5GCN3 --> 0", "0 --> R5GCN4", "R5GCN4 --> 0", "0 --> R6GCN1", "R6GCN1 --> 0", "0 --> R6GCN2", "R6GCN2 --> 0", "0 --> R6GCN3", "R6GCN3 --> 0", "0 --> R6GCN4", "R6GCN4 --> 0", "0 --> R7GCN1", "R7GCN1 --> 0", "0 --> R7GCN2", "R7GCN2 --> 0", "0 --> R7GCN3", "R7GCN3 --> 0", "0 --> R7GCN4", "R7GCN4 --> 0", "0 --> R8GCN1", "R8GCN1 --> 0", "0 --> R8GCN2", "R8GCN2 --> 0", "0 --> R8GCN3", "R8GCN3 --> 0", "0 --> R8GCN4", "R8GCN4 --> 0", "0 --> R9GCN1", "R9GCN1 --> 0", "0 --> R9GCN2", "R9GCN2 --> 0", "0 --> R9GCN3", "R9GCN3 --> 0", "0 --> R9GCN4", "R9GCN4 --> 0", "0 --> R10GCN1", "R10GCN1 --> 0", "0 --> R10GCN2", "R10GCN2 --> 0", "0 --> R10GCN3", "R10GCN3 --> 0", "0 --> R10GCN4", "R10GCN4 --> 0", "R1GCN1 --> R1GCN1 + P1GCN1", "P1GCN1 --> 0", "R1GCN2 --> R1GCN2 + P1GCN2", "P1GCN2 --> 0", "R1GCN3 --> R1GCN3 + P1GCN3", "P1GCN3 --> 0", "R1GCN4 --> R1GCN4 + P1GCN4", "P1GCN4 --> 0", "R2GCN1 --> R2GCN1 + P2GCN1", "P2GCN1 --> 0", "R2GCN2 --> R2GCN2 + P2GCN2", "P2GCN2 --> 0", "R2GCN3 --> R2GCN3 + P2GCN3", "P2GCN3 --> 0", "R2GCN4 --> R2GCN4 + P2GCN4", "P2GCN4 --> 0", "R3GCN1 --> R3GCN1 + P3GCN1", "P3GCN1 --> 0", "R3GCN2 --> R3GCN2 + P3GCN2", "P3GCN2 --> 0", "R3GCN3 --> R3GCN3 + P3GCN3", "P3GCN3 --> 0", "R3GCN4 --> R3GCN4 + P3GCN4", "P3GCN4 --> 0", "R4GCN1 --> R4GCN1 + P4GCN1", "P4GCN1 --> 0", "R4GCN2 --> R4GCN2 + P4GCN2", "P4GCN2 --> 0", "R4GCN3 --> R4GCN3 + P4GCN3", "P4GCN3 --> 0", "R4GCN4 --> R4GCN4 + P4GCN4", "P4GCN4 --> 0", "R5GCN1 --> R5GCN1 + P5GCN1", "P5GCN1 --> 0", "R5GCN2 --> R5GCN2 + P5GCN2", "P5GCN2 --> 0", "R5GCN3 --> R5GCN3 + P5GCN3", "P5GCN3 --> 0", "R5GCN4 --> R5GCN4 + P5GCN4", "P5GCN4 --> 0", "R7GCN1 --> R7GCN1 + P7GCN1", "P7GCN1 --> 0", "R7GCN2 --> R7GCN2 + P7GCN2", "P7GCN2 --> 0", "R7GCN3 --> R7GCN3 + P7GCN3", "P7GCN3 --> 0", "R7GCN4 --> R7GCN4 + P7GCN4", "P7GCN4 --> 0"],
    "initialconditions"=>Any[:(((28.8461 * ((QTLeffects["GCN1"])["qtlTCrate"])[1]) / ((QTLeffects["GCN1"])["qtlRDrate"])[1]) * ((Initvar["GCN1"])["R"])[1]), :(((28.8461 * ((QTLeffects["GCN2"])["qtlTCrate"])[1]) / ((QTLeffects["GCN2"])["qtlRDrate"])[1]) * ((Initvar["GCN2"])["R"])[1]), :(((28.8461 * ((QTLeffects["GCN3"])["qtlTCrate"])[1]) / ((QTLeffects["GCN3"])["qtlRDrate"])[1]) * ((Initvar["GCN3"])["R"])[1]), :(((28.8461 * ((QTLeffects["GCN4"])["qtlTCrate"])[1]) / ((QTLeffects["GCN4"])["qtlRDrate"])[1]) * ((Initvar["GCN4"])["R"])[1]), :(((81.3848 * ((QTLeffects["GCN1"])["qtlTCrate"])[2]) / ((QTLeffects["GCN1"])["qtlRDrate"])[2]) * ((Initvar["GCN1"])["R"])[2]), :(((81.3848 * ((QTLeffects["GCN2"])["qtlTCrate"])[2]) / ((QTLeffects["GCN2"])["qtlRDrate"])[2]) * ((Initvar["GCN2"])["R"])[2]), :(((81.3848 * ((QTLeffects["GCN3"])["qtlTCrate"])[2]) / ((QTLeffects["GCN3"])["qtlRDrate"])[2]) * ((Initvar["GCN3"])["R"])[2]), :(((81.3848 * ((QTLeffects["GCN4"])["qtlTCrate"])[2]) / ((QTLeffects["GCN4"])["qtlRDrate"])[2]) * ((Initvar["GCN4"])["R"])[2]), :(((24.6875 * ((QTLeffects["GCN1"])["qtlTCrate"])[3]) / ((QTLeffects["GCN1"])["qtlRDrate"])[3]) * ((Initvar["GCN1"])["R"])[3]), :(((24.6875 * ((QTLeffects["GCN2"])["qtlTCrate"])[3]) / ((QTLeffects["GCN2"])["qtlRDrate"])[3]) * ((Initvar["GCN2"])["R"])[3]), :(((24.6875 * ((QTLeffects["GCN3"])["qtlTCrate"])[3]) / ((QTLeffects["GCN3"])["qtlRDrate"])[3]) * ((Initvar["GCN3"])["R"])[3]), :(((24.6875 * ((QTLeffects["GCN4"])["qtlTCrate"])[3]) / ((QTLeffects["GCN4"])["qtlRDrate"])[3]) * ((Initvar["GCN4"])["R"])[3]), :(((313.326 * ((QTLeffects["GCN1"])["qtlTCrate"])[4]) / ((QTLeffects["GCN1"])["qtlRDrate"])[4]) * ((Initvar["GCN1"])["R"])[4]), :(((313.326 * ((QTLeffects["GCN2"])["qtlTCrate"])[4]) / ((QTLeffects["GCN2"])["qtlRDrate"])[4]) * ((Initvar["GCN2"])["R"])[4]), :(((313.326 * ((QTLeffects["GCN3"])["qtlTCrate"])[4]) / ((QTLeffects["GCN3"])["qtlRDrate"])[4]) * ((Initvar["GCN3"])["R"])[4]), :(((313.326 * ((QTLeffects["GCN4"])["qtlTCrate"])[4]) / ((QTLeffects["GCN4"])["qtlRDrate"])[4]) * ((Initvar["GCN4"])["R"])[4]), :(((29.0011 * ((QTLeffects["GCN1"])["qtlTCrate"])[5]) / ((QTLeffects["GCN1"])["qtlRDrate"])[5]) * ((Initvar["GCN1"])["R"])[5]), :(((29.0011 * ((QTLeffects["GCN2"])["qtlTCrate"])[5]) / ((QTLeffects["GCN2"])["qtlRDrate"])[5]) * ((Initvar["GCN2"])["R"])[5]), :(((29.0011 * ((QTLeffects["GCN3"])["qtlTCrate"])[5]) / ((QTLeffects["GCN3"])["qtlRDrate"])[5]) * ((Initvar["GCN3"])["R"])[5]), :(((29.0011 * ((QTLeffects["GCN4"])["qtlTCrate"])[5]) / ((QTLeffects["GCN4"])["qtlRDrate"])[5]) * ((Initvar["GCN4"])["R"])[5]), :(((164.491 * ((QTLeffects["GCN1"])["qtlTCrate"])[6]) / ((QTLeffects["GCN1"])["qtlRDrate"])[6]) * ((Initvar["GCN1"])["R"])[6]), :(((164.491 * ((QTLeffects["GCN2"])["qtlTCrate"])[6]) / ((QTLeffects["GCN2"])["qtlRDrate"])[6]) * ((Initvar["GCN2"])["R"])[6]), :(((164.491 * ((QTLeffects["GCN3"])["qtlTCrate"])[6]) / ((QTLeffects["GCN3"])["qtlRDrate"])[6]) * ((Initvar["GCN3"])["R"])[6]), :(((164.491 * ((QTLeffects["GCN4"])["qtlTCrate"])[6]) / ((QTLeffects["GCN4"])["qtlRDrate"])[6]) * ((Initvar["GCN4"])["R"])[6]), :(((71.5497 * ((QTLeffects["GCN1"])["qtlTCrate"])[7]) / ((QTLeffects["GCN1"])["qtlRDrate"])[7]) * ((Initvar["GCN1"])["R"])[7]), :(((71.5497 * ((QTLeffects["GCN2"])["qtlTCrate"])[7]) / ((QTLeffects["GCN2"])["qtlRDrate"])[7]) * ((Initvar["GCN2"])["R"])[7]), :(((71.5497 * ((QTLeffects["GCN3"])["qtlTCrate"])[7]) / ((QTLeffects["GCN3"])["qtlRDrate"])[7]) * ((Initvar["GCN3"])["R"])[7]), :(((71.5497 * ((QTLeffects["GCN4"])["qtlTCrate"])[7]) / ((QTLeffects["GCN4"])["qtlRDrate"])[7]) * ((Initvar["GCN4"])["R"])[7]), :(((191.649 * ((QTLeffects["GCN1"])["qtlTCrate"])[8]) / ((QTLeffects["GCN1"])["qtlRDrate"])[8]) * ((Initvar["GCN1"])["R"])[8]), :(((191.649 * ((QTLeffects["GCN2"])["qtlTCrate"])[8]) / ((QTLeffects["GCN2"])["qtlRDrate"])[8]) * ((Initvar["GCN2"])["R"])[8]), :(((191.649 * ((QTLeffects["GCN3"])["qtlTCrate"])[8]) / ((QTLeffects["GCN3"])["qtlRDrate"])[8]) * ((Initvar["GCN3"])["R"])[8]), :(((191.649 * ((QTLeffects["GCN4"])["qtlTCrate"])[8]) / ((QTLeffects["GCN4"])["qtlRDrate"])[8]) * ((Initvar["GCN4"])["R"])[8]), :(((79.9315 * ((QTLeffects["GCN1"])["qtlTCrate"])[9]) / ((QTLeffects["GCN1"])["qtlRDrate"])[9]) * ((Initvar["GCN1"])["R"])[9]), :(((79.9315 * ((QTLeffects["GCN2"])["qtlTCrate"])[9]) / ((QTLeffects["GCN2"])["qtlRDrate"])[9]) * ((Initvar["GCN2"])["R"])[9]), :(((79.9315 * ((QTLeffects["GCN3"])["qtlTCrate"])[9]) / ((QTLeffects["GCN3"])["qtlRDrate"])[9]) * ((Initvar["GCN3"])["R"])[9]), :(((79.9315 * ((QTLeffects["GCN4"])["qtlTCrate"])[9]) / ((QTLeffects["GCN4"])["qtlRDrate"])[9]) * ((Initvar["GCN4"])["R"])[9]), :(((293.694 * ((QTLeffects["GCN1"])["qtlTCrate"])[10]) / ((QTLeffects["GCN1"])["qtlRDrate"])[10]) * ((Initvar["GCN1"])["R"])[10]), :(((293.694 * ((QTLeffects["GCN2"])["qtlTCrate"])[10]) / ((QTLeffects["GCN2"])["qtlRDrate"])[10]) * ((Initvar["GCN2"])["R"])[10]), :(((293.694 * ((QTLeffects["GCN3"])["qtlTCrate"])[10]) / ((QTLeffects["GCN3"])["qtlRDrate"])[10]) * ((Initvar["GCN3"])["R"])[10]), :(((293.694 * ((QTLeffects["GCN4"])["qtlTCrate"])[10]) / ((QTLeffects["GCN4"])["qtlRDrate"])[10]) * ((Initvar["GCN4"])["R"])[10]), :(((1.26595e6 * ((QTLeffects["GCN1"])["qtlTCrate"])[1] * ((QTLeffects["GCN1"])["qtlTLrate"])[1]) / (((QTLeffects["GCN1"])["qtlRDrate"])[1] * ((QTLeffects["GCN1"])["qtlPDrate"])[1])) * ((Initvar["GCN1"])["P"])[1]), :(((1.26595e6 * ((QTLeffects["GCN2"])["qtlTCrate"])[1] * ((QTLeffects["GCN2"])["qtlTLrate"])[1]) / (((QTLeffects["GCN2"])["qtlRDrate"])[1] * ((QTLeffects["GCN2"])["qtlPDrate"])[1])) * ((Initvar["GCN2"])["P"])[1]), :(((1.26595e6 * ((QTLeffects["GCN3"])["qtlTCrate"])[1] * ((QTLeffects["GCN3"])["qtlTLrate"])[1]) / (((QTLeffects["GCN3"])["qtlRDrate"])[1] * ((QTLeffects["GCN3"])["qtlPDrate"])[1])) * ((Initvar["GCN3"])["P"])[1]), :(((1.26595e6 * ((QTLeffects["GCN4"])["qtlTCrate"])[1] * ((QTLeffects["GCN4"])["qtlTLrate"])[1]) / (((QTLeffects["GCN4"])["qtlRDrate"])[1] * ((QTLeffects["GCN4"])["qtlPDrate"])[1])) * ((Initvar["GCN4"])["P"])[1]), :(((1.97821e6 * ((QTLeffects["GCN1"])["qtlTCrate"])[2] * ((QTLeffects["GCN1"])["qtlTLrate"])[2]) / (((QTLeffects["GCN1"])["qtlRDrate"])[2] * ((QTLeffects["GCN1"])["qtlPDrate"])[2])) * ((Initvar["GCN1"])["P"])[2]), :(((1.97821e6 * ((QTLeffects["GCN2"])["qtlTCrate"])[2] * ((QTLeffects["GCN2"])["qtlTLrate"])[2]) / (((QTLeffects["GCN2"])["qtlRDrate"])[2] * ((QTLeffects["GCN2"])["qtlPDrate"])[2])) * ((Initvar["GCN2"])["P"])[2]), :(((1.97821e6 * ((QTLeffects["GCN3"])["qtlTCrate"])[2] * ((QTLeffects["GCN3"])["qtlTLrate"])[2]) / (((QTLeffects["GCN3"])["qtlRDrate"])[2] * ((QTLeffects["GCN3"])["qtlPDrate"])[2])) * ((Initvar["GCN3"])["P"])[2]), :(((1.97821e6 * ((QTLeffects["GCN4"])["qtlTCrate"])[2] * ((QTLeffects["GCN4"])["qtlTLrate"])[2]) / (((QTLeffects["GCN4"])["qtlRDrate"])[2] * ((QTLeffects["GCN4"])["qtlPDrate"])[2])) * ((Initvar["GCN4"])["P"])[2]), :(((2.77873e5 * ((QTLeffects["GCN1"])["qtlTCrate"])[3] * ((QTLeffects["GCN1"])["qtlTLrate"])[3]) / (((QTLeffects["GCN1"])["qtlRDrate"])[3] * ((QTLeffects["GCN1"])["qtlPDrate"])[3])) * ((Initvar["GCN1"])["P"])[3]), :(((2.77873e5 * ((QTLeffects["GCN2"])["qtlTCrate"])[3] * ((QTLeffects["GCN2"])["qtlTLrate"])[3]) / (((QTLeffects["GCN2"])["qtlRDrate"])[3] * ((QTLeffects["GCN2"])["qtlPDrate"])[3])) * ((Initvar["GCN2"])["P"])[3]), :(((2.77873e5 * ((QTLeffects["GCN3"])["qtlTCrate"])[3] * ((QTLeffects["GCN3"])["qtlTLrate"])[3]) / (((QTLeffects["GCN3"])["qtlRDrate"])[3] * ((QTLeffects["GCN3"])["qtlPDrate"])[3])) * ((Initvar["GCN3"])["P"])[3]), :(((2.77873e5 * ((QTLeffects["GCN4"])["qtlTCrate"])[3] * ((QTLeffects["GCN4"])["qtlTLrate"])[3]) / (((QTLeffects["GCN4"])["qtlRDrate"])[3] * ((QTLeffects["GCN4"])["qtlPDrate"])[3])) * ((Initvar["GCN4"])["P"])[3]), :(((1.10747e7 * ((QTLeffects["GCN1"])["qtlTCrate"])[4] * ((QTLeffects["GCN1"])["qtlTLrate"])[4]) / (((QTLeffects["GCN1"])["qtlRDrate"])[4] * ((QTLeffects["GCN1"])["qtlPDrate"])[4])) * ((Initvar["GCN1"])["P"])[4]), :(((1.10747e7 * ((QTLeffects["GCN2"])["qtlTCrate"])[4] * ((QTLeffects["GCN2"])["qtlTLrate"])[4]) / (((QTLeffects["GCN2"])["qtlRDrate"])[4] * ((QTLeffects["GCN2"])["qtlPDrate"])[4])) * ((Initvar["GCN2"])["P"])[4]), :(((1.10747e7 * ((QTLeffects["GCN3"])["qtlTCrate"])[4] * ((QTLeffects["GCN3"])["qtlTLrate"])[4]) / (((QTLeffects["GCN3"])["qtlRDrate"])[4] * ((QTLeffects["GCN3"])["qtlPDrate"])[4])) * ((Initvar["GCN3"])["P"])[4]), :(((1.10747e7 * ((QTLeffects["GCN4"])["qtlTCrate"])[4] * ((QTLeffects["GCN4"])["qtlTLrate"])[4]) / (((QTLeffects["GCN4"])["qtlRDrate"])[4] * ((QTLeffects["GCN4"])["qtlPDrate"])[4])) * ((Initvar["GCN4"])["P"])[4]), :(((1.7921e5 * ((QTLeffects["GCN1"])["qtlTCrate"])[5] * ((QTLeffects["GCN1"])["qtlTLrate"])[5]) / (((QTLeffects["GCN1"])["qtlRDrate"])[5] * ((QTLeffects["GCN1"])["qtlPDrate"])[5])) * ((Initvar["GCN1"])["P"])[5]), :(((1.7921e5 * ((QTLeffects["GCN2"])["qtlTCrate"])[5] * ((QTLeffects["GCN2"])["qtlTLrate"])[5]) / (((QTLeffects["GCN2"])["qtlRDrate"])[5] * ((QTLeffects["GCN2"])["qtlPDrate"])[5])) * ((Initvar["GCN2"])["P"])[5]), :(((1.7921e5 * ((QTLeffects["GCN3"])["qtlTCrate"])[5] * ((QTLeffects["GCN3"])["qtlTLrate"])[5]) / (((QTLeffects["GCN3"])["qtlRDrate"])[5] * ((QTLeffects["GCN3"])["qtlPDrate"])[5])) * ((Initvar["GCN3"])["P"])[5]), :(((1.7921e5 * ((QTLeffects["GCN4"])["qtlTCrate"])[5] * ((QTLeffects["GCN4"])["qtlTLrate"])[5]) / (((QTLeffects["GCN4"])["qtlRDrate"])[5] * ((QTLeffects["GCN4"])["qtlPDrate"])[5])) * ((Initvar["GCN4"])["P"])[5]), :(((3.17649e6 * ((QTLeffects["GCN1"])["qtlTCrate"])[7] * ((QTLeffects["GCN1"])["qtlTLrate"])[7]) / (((QTLeffects["GCN1"])["qtlRDrate"])[7] * ((QTLeffects["GCN1"])["qtlPDrate"])[7])) * ((Initvar["GCN1"])["P"])[7]), :(((3.17649e6 * ((QTLeffects["GCN2"])["qtlTCrate"])[7] * ((QTLeffects["GCN2"])["qtlTLrate"])[7]) / (((QTLeffects["GCN2"])["qtlRDrate"])[7] * ((QTLeffects["GCN2"])["qtlPDrate"])[7])) * ((Initvar["GCN2"])["P"])[7]), :(((3.17649e6 * ((QTLeffects["GCN3"])["qtlTCrate"])[7] * ((QTLeffects["GCN3"])["qtlTLrate"])[7]) / (((QTLeffects["GCN3"])["qtlRDrate"])[7] * ((QTLeffects["GCN3"])["qtlPDrate"])[7])) * ((Initvar["GCN3"])["P"])[7]), :(((3.17649e6 * ((QTLeffects["GCN4"])["qtlTCrate"])[7] * ((QTLeffects["GCN4"])["qtlTLrate"])[7]) / (((QTLeffects["GCN4"])["qtlRDrate"])[7] * ((QTLeffects["GCN4"])["qtlPDrate"])[7])) * ((Initvar["GCN4"])["P"])[7])],"species"=>Any["R1GCN1", "R1GCN2", "R1GCN3", "R1GCN4", "R2GCN1", "R2GCN2", "R2GCN3", "R2GCN4", "R3GCN1", "R3GCN2", "R3GCN3", "R3GCN4", "R4GCN1", "R4GCN2", "R4GCN3", "R4GCN4", "R5GCN1", "R5GCN2", "R5GCN3", "R5GCN4", "R6GCN1", "R6GCN2", "R6GCN3", "R6GCN4", "R7GCN1", "R7GCN2", "R7GCN3", "R7GCN4", "R8GCN1", "R8GCN2", "R8GCN3", "R8GCN4", "R9GCN1", "R9GCN2", "R9GCN3", "R9GCN4", "R10GCN1", "R10GCN2", "R10GCN3", "R10GCN4", "P1GCN1", "P1GCN2", "P1GCN3", "P1GCN4", "P2GCN1", "P2GCN2", "P2GCN3", "P2GCN4", "P3GCN1", "P3GCN2", "P3GCN3", "P3GCN4", "P4GCN1", "P4GCN2", "P4GCN3", "P4GCN4", "P5GCN1", "P5GCN2", "P5GCN3", "P5GCN4", "P7GCN1", "P7GCN2", "P7GCN3", "P7GCN4"])

QTLeffects = Dict("GCN3"=>Dict("qtlPDregbind"=>[1.0, 0.903064, 1.0, 1.01958, 0.908109, 0.0, 0.985851, 0.0, 0.0, 0.0],"qtlTLrate"=>[1.0, 0.98387, 1.12259, 0.947592, 1.02447, 0.0, 0.897265, 0.0, 0.0, 0.0],"qtlTCrate"=>[1.0, 1.17192, 1.0, 0.973951, 1.07116, 1.11191, 1.0003, 1.04185, 0.96825, 1.0],"qtlRDbindreg"=>[1.0, 0.969047, 0.883935, 1.16487, 1.18382, 0.926073, 1.13709, 1.16358, 1.20015, 1.0],"qtlTCregbind"=>[1.0, 0.939463, 1.0, 0.996446, 1.11144, 1.0, 0.92738, 1.0, 1.11756, 1.0],"qtlactivity"=>[1.0, 1.07163, 1.0, 0.948701, 1.00975, 1.0, 1.00126, 0.875307, 1.03423, 1.0],"qtlRDrate"=>[1.0, 1.0, 1.0, 1.08387, 1.02613, 1.0, 1.20753, 1.06082, 1.0, 1.0],"qtlTLregbind"=>[1.0, 0.969293, 1.0, 0.9764, 0.932668, 0.0, 0.918416, 0.0, 0.0, 0.0],"qtlPDrate"=>[1.0, 1.14071, 1.0, 0.93931, 1.02331, 0.0, 1.05153, 0.0, 0.0, 0.0]),"GCN2"=>Dict("qtlPDregbind"=>[1.0, 0.903064, 1.0, 0.904002, 1.0, 0.0, 0.947351, 0.0, 0.0, 0.0],"qtlTLrate"=>[0.888289, 0.98387, 1.0, 0.633331, 0.927585, 0.0, 1.03055, 0.0, 0.0, 0.0],"qtlTCrate"=>[1.0, 1.17192, 1.0, 1.12524, 1.0, 1.0, 0.993596, 0.99752, 1.0, 1.0],"qtlRDbindreg"=>[1.0, 0.969047, 1.0, 1.04387, 1.0, 0.888253, 1.0, 0.843412, 1.0, 1.0],"qtlTCregbind"=>[1.0, 0.939463, 1.0, 1.15848, 1.0, 0.919143, 1.0, 0.813784, 0.975534, 1.0],"qtlactivity"=>[1.0, 1.07163, 1.0, 0.864578, 1.0, 1.13318, 0.918793, 1.08566, 1.0, 1.0],"qtlRDrate"=>[1.0, 1.0, 1.0, 1.01624, 1.0, 1.0, 1.0, 0.926551, 0.941337, 1.0],"qtlTLregbind"=>[1.0, 0.969293, 1.0, 1.05993, 1.0, 0.0, 0.992744, 0.0, 0.0, 0.0],"qtlPDrate"=>[1.0, 1.14071, 1.0, 1.00784, 1.0, 0.0, 0.990836, 0.0, 0.0, 0.0]),"GCN4"=>Dict("qtlPDregbind"=>[1.0, 0.903064, 1.0, 1.0, 1.0, 0.0, 0.985851, 0.0, 0.0, 0.0],"qtlTLrate"=>[0.888289, 0.98387, 1.04376, 1.0, 1.0, 0.0, 0.897265, 0.0, 0.0, 0.0],"qtlTCrate"=>[1.0, 1.17192, 1.0, 1.0, 1.0, 1.0, 1.0003, 1.0, 0.934712, 1.09404],"qtlRDbindreg"=>[1.0, 0.969047, 1.08212, 1.0, 1.0, 1.0, 1.13709, 1.0, 0.953415, 1.08835],"qtlTCregbind"=>[1.0, 0.939463, 1.0, 1.0, 1.0, 1.0, 0.92738, 1.0, 1.0, 0.874129],"qtlactivity"=>[1.0, 1.07163, 0.985766, 1.0, 1.0, 1.0, 1.00126, 1.0, 0.87951, 0.774213],"qtlRDrate"=>[1.0, 1.0, 0.923331, 1.0, 1.0, 1.0, 1.20753, 1.0, 1.0, 1.0],"qtlTLregbind"=>[1.0, 0.969293, 0.947754, 0.790886, 1.0, 0.0, 0.918416, 0.0, 0.0, 0.0],"qtlPDrate"=>[1.0, 1.14071, 1.0, 1.0, 1.0, 0.0, 1.05153, 0.0, 0.0, 0.0]),"GCN1"=>Dict("qtlPDregbind"=>[1.15029, 1.16139, 1.0, 1.0, 0.965733, 0.0, 1.0, 0.0, 0.0, 0.0],"qtlTLrate"=>[1.0, 0.961051, 1.12259, 1.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0],"qtlTCrate"=>[1.0, 1.24792, 1.0, 0.894408, 1.0, 1.0, 1.0, 0.99752, 0.96825, 1.0],"qtlRDbindreg"=>[1.0, 0.83549, 0.883935, 1.07553, 1.0, 0.870481, 1.0, 0.843412, 1.20015, 1.0],"qtlTCregbind"=>[1.0, 1.08633, 1.0, 1.06498, 1.03092, 1.0, 1.0, 0.813784, 1.11756, 1.0],"qtlactivity"=>[1.0, 0.791183, 1.0, 1.0, 0.948967, 1.0, 1.0, 1.08566, 1.03423, 1.0],"qtlRDrate"=>[1.0, 0.818898, 1.0, 1.16617, 1.09283, 1.0, 1.0, 0.926551, 1.0, 1.0],"qtlTLregbind"=>[1.0, 0.9524, 1.0, 0.970783, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0],"qtlPDrate"=>[1.11645, 0.940561, 1.0, 0.949767, 0.837376, 0.0, 1.0, 0.0, 0.0, 0.0]))
Initvar = Dict("GCN3"=>Dict("P"=>[0.922647, 0.827737, 0.962663, 1.07635, 0.848211, 1.01387, 1.15418, 0.997701, 0.989415, 1.06221],"R"=>[1.01076, 0.903517, 0.931234, 0.948392, 1.08051, 1.07374, 0.843925, 0.956626, 1.1761, 1.09528]),"GCN2"=>Dict("P"=>[0.853842, 0.93087, 1.12092, 0.964489, 1.08834, 1.05661, 0.833048, 1.06602, 1.16125, 1.10232],"R"=>[0.97554, 0.932472, 1.0511, 0.898645, 1.05261, 1.05966, 0.7716, 0.90272, 0.807083, 1.06358]),"GCN4"=>Dict("P"=>[0.990594, 1.09941, 1.0421, 0.930396, 0.81033, 0.90876, 0.86762, 1.13135, 0.93453, 1.04141],"R"=>[1.00492, 1.1564, 1.14431, 0.93374, 0.970863, 1.13996, 1.1174, 0.978497, 1.10544, 1.00269]),"GCN1"=>Dict("P"=>[1.07192, 1.15348, 1.06973, 0.926324, 1.04767, 0.955735, 1.08041, 0.796943, 0.988385, 0.919596],"R"=>[0.855099, 0.818183, 0.857508, 0.910632, 1.00853, 0.99178, 0.934635, 0.879292, 0.843702, 0.908478]))

model = Network("coucou")

for i in 1:length(stochmodel["species"])
    i0 = round(Int, eval(stochmodel["initialconditions"][i]))
    #println(stochmodel["species"][i]* "\t"*string(i0))
    model <= BioSimulator.Species(stochmodel["species"][i], i0)
end

for i in eachindex(stochmodel["reactions"])
    #println(join([stochmodel["reactionsnames"][i], eval(stochmodel["propensities"][i]), stochmodel["reactions"][i]], "\t", "\n"))
    model <= BioSimulator.Reaction(stochmodel["reactionsnames"][i], eval(stochmodel["propensities"][i]), stochmodel["reactions"][i])
end

result = simulate(model, algorithm=SSA, time = 1000.0, epochs = 1000, trials = 1)

test = res2df(result)
plot(test[:, :time], test[:, :R1GCN1])

genes = collect(1:10)
gcnList = ["GCN1", "GCN2", "GCN3", "GCN4"]

for g in genes
##    toplotRNA = [test[t, Symbol("R"*string(g).*gcn)] for t in 1:size(test)[1], gcn in gcnList]
##    toplotProt = [test[t, Symbol("P"*string(g).*gcn)] for t in 1:size(test)[1], gcn in gcnList]
##    gui(plot(test[:, :time], hcat(toplotRNA, toplotProt), show = true, reuse = false))

    toplotRNA = [test[t, Symbol("R"*string(g).*gcn)] for t in 1:size(test)[1], gcn in gcnList]
    if in(g, [1,2,3,4,5,7])
        toplotProt = [test[t, Symbol("P"*string(g).*gcn)] for t in 1:size(test)[1], gcn in gcnList]
        gui(Plots.plot(test[:, :time], hcat(toplotRNA, toplotProt)[:, [1,5,2,6,3,7,4,8]], layout = @layout([a b]), show = true, reuse = false))
    else
        gui(Plots.plot(test[:, :time], toplotRNA, show = true, reuse = false))
    end
end

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

QTLeffects = Dict("GCN3"=>Dict("qtlPDregbind"=>[0.0, 1.0, 0.996807, 0.868407, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0],"qtlTLrate"=>[0.0, 1.06523, 1.04637, 0.996554, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0],"qtlTCrate"=>[1.00504, 1.0, 1.0, 1.0519, 1.0, 0.939048, 1.0, 1.0, 1.0, 1.0],"qtlRDbindreg"=>[1.01377, 1.0, 1.15503, 1.0, 1.16826, 1.01473, 1.0, 1.0, 1.0, 1.0],"qtlTCregbind"=>[0.980469, 1.0, 0.957135, 1.14658, 1.0, 0.987476, 1.0, 1.0, 1.0, 1.0],"qtlactivity"=>[0.914794, 1.0, 0.949882, 0.923203, 1.0, 1.01069, 1.0, 1.0, 1.0, 1.0],"qtlRDrate"=>[0.921899, 1.0, 0.988992, 0.899535, 1.0, 0.912394, 1.0, 1.0, 1.0, 1.0],"qtlTLregbind"=>[0.0, 1.0, 0.933088, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0],"qtlPDrate"=>[0.0, 1.0, 1.00773, 1.0762, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0]),"GCN2"=>Dict("qtlPDregbind"=>[0.0, 1.0, 0.996807, 0.868407, 0.0, 0.0, 1.0, 0.97277, 1.0, 1.0],"qtlTLrate"=>[0.0, 1.06523, 1.04637, 0.996554, 0.0, 0.0, 1.0, 1.12512, 1.0, 1.0],"qtlTCrate"=>[1.00504, 1.0, 1.0, 1.0519, 1.09237, 1.08641, 1.0, 0.901656, 1.04507, 1.0],"qtlRDbindreg"=>[1.01377, 1.0, 1.15503, 1.0, 0.959294, 0.818953, 1.0, 1.0482, 1.0, 1.0],"qtlTCregbind"=>[0.980469, 1.0, 0.957135, 1.14658, 1.22001, 1.02222, 1.0, 1.15073, 1.0, 1.0],"qtlactivity"=>[0.914794, 1.0, 0.949882, 0.923203, 1.07234, 1.01743, 1.0, 0.974383, 1.0, 1.0],"qtlRDrate"=>[0.921899, 1.0, 0.988992, 0.899535, 1.02111, 1.04598, 1.0, 1.06141, 1.01365, 1.0],"qtlTLregbind"=>[0.0, 1.0, 0.933088, 1.0, 0.0, 0.0, 1.0, 1.07071, 1.0, 1.0],"qtlPDrate"=>[0.0, 1.0, 1.00773, 1.0762, 0.0, 0.0, 1.0, 1.03252, 1.0, 1.0]),"GCN4"=>Dict("qtlPDregbind"=>[0.0, 1.0, 1.13171, 0.868407, 0.0, 0.0, 1.05037, 1.0, 1.0, 1.0],"qtlTLrate"=>[0.0, 1.0, 0.969644, 0.996554, 0.0, 0.0, 1.07574, 1.0, 0.956107, 1.0],"qtlTCrate"=>[1.0, 1.0, 1.0, 1.0519, 1.09237, 0.939048, 0.922822, 1.0, 1.0, 1.0],"qtlRDbindreg"=>[1.0, 1.0, 1.0, 1.0, 0.959294, 1.01473, 0.927262, 1.0, 0.853407, 1.0],"qtlTCregbind"=>[1.0, 1.0, 1.0, 1.14658, 1.22001, 0.987476, 0.931289, 1.0, 1.0, 1.0],"qtlactivity"=>[1.0, 1.0, 1.09505, 0.923203, 1.07234, 1.01069, 0.902183, 1.0, 1.03853, 1.0],"qtlRDrate"=>[1.0, 1.0, 1.0, 0.899535, 1.02111, 0.912394, 0.836054, 1.0, 1.0, 1.0],"qtlTLregbind"=>[0.0, 1.0, 1.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0],"qtlPDrate"=>[0.0, 1.0, 1.03975, 1.0762, 0.0, 0.0, 0.92047, 1.0, 1.0, 1.0]),"GCN1"=>Dict("qtlPDregbind"=>[0.0, 1.0, 1.13171, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0],"qtlTLrate"=>[0.0, 1.05908, 0.969644, 1.06483, 0.0, 0.0, 1.0, 0.796604, 1.0, 1.0],"qtlTCrate"=>[1.08969, 1.0, 1.0, 1.0, 1.14363, 1.04624, 1.0, 1.01885, 1.04507, 1.0],"qtlRDbindreg"=>[1.08142, 0.991458, 1.0, 1.0, 0.959343, 1.0, 1.0, 0.980113, 1.0, 1.0],"qtlTCregbind"=>[1.0, 1.0, 1.0, 0.778718, 1.0, 1.01828, 1.0, 1.15897, 1.0, 1.0],"qtlactivity"=>[0.94031, 1.0, 1.09505, 0.946715, 1.0, 1.0, 1.0, 0.945421, 1.0, 1.0],"qtlRDrate"=>[1.0, 1.09692, 1.0, 1.0, 1.0, 1.04338, 1.0, 0.896263, 1.01365, 1.0],"qtlTLregbind"=>[0.0, 1.0, 1.0, 1.0, 0.0, 0.0, 1.05649, 1.03906, 1.0, 1.0],"qtlPDrate"=>[0.0, 1.0, 1.03975, 1.0, 0.0, 0.0, 1.0, 1.01075, 1.0, 1.0]))

InitVar = Dict("GCN3"=>Dict("P"=>[0.80262, 0.96886, 0.928253, 0.922065, 0.879601, 0.990862, 1.05296, 0.888656, 1.17134, 1.09654],"R"=>[0.904391, 1.08354, 1.1342, 1.06773, 0.962775, 1.08633, 1.09988, 1.04768, 1.09627, 0.980728]),"GCN2"=>Dict("P"=>[1.03636, 0.852727, 0.907479, 1.13289, 0.689987, 1.1046, 1.072, 1.10135, 1.00242, 0.848639],"R"=>[0.889878, 1.03319, 1.03879, 0.960242, 0.953933, 0.95715, 0.939481, 0.883476, 0.851746, 1.15431]),"GCN4"=>Dict("P"=>[1.07095, 1.09463, 1.1466, 0.954728, 0.901106, 1.17743, 1.14382, 0.993785, 1.14216, 0.974802],"R"=>[1.0076, 0.996241, 1.21269, 0.931376, 1.0446, 0.95608, 0.964789, 1.12328, 1.05084, 1.08357]),"GCN1"=>Dict("P"=>[1.09452, 1.07919, 0.866409, 1.09506, 1.03763, 1.09511, 0.907524, 0.962947, 0.867508, 1.17224],"R"=>[1.16395, 0.96986, 0.988619, 0.992444, 0.940217, 0.985759, 0.777827, 0.929596, 0.832051, 0.883847]))

nod = Dict{String,Any}(Pair{String,Any}("TargetReaction", String["TL", "TL", "TC", "TC", "TL", "RD", "TC", "TL", "TL", "TL"]),Pair{String,Any}("TCrate", [0.0876407, 0.0590708, 0.0874417, 0.017104, 0.0708739, 0.0741194, 0.0804994, 0.0102489, 0.0217458, 0.0622714]),Pair{String,Any}("RDrate", [0.000417537, 0.000398565, 0.000361141, 0.00048852, 0.000350877, 0.000358551, 0.00458716, 0.00129199, 0.000702741, 0.000365364]),Pair{String,Any}("id", [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]),Pair{String,Any}("ActiveForm", String["R1", "Pm2", "Pm3", "P4", "R5", "R6", "Pm7", "Pm8", "P9", "Pm10"]),Pair{String,Any}("coding", String["NC", "PC", "PC", "PC", "NC", "NC", "PC", "PC", "PC", "PC"]),Pair{String,Any}("TLrate", [0.0, 4.70877, 2.47364, 3.30288, 0.0, 0.0, 4.7373, 2.23753, 1.76261, 3.7403]),Pair{String,Any}("PTMform", String["0", "1", "1", "0", "0", "0", "1", "1", "0", "1"]),Pair{String,Any}("PDrate", [0.0, 8.8739e-5, 0.000151953, 7.41565e-5, 0.0, 0.0, 8.75197e-5, 9.34842e-5, 8.85897e-5, 0.000114916]))

simtime = 100
nepochs = -1
ntrials = 1
simalgorithm = "SSA"
modelname = "coucou"

genes = collect(1:10)
gcnList = ["GCN1", "GCN2", "GCN3", "GCN4"]

for g in genes
##    toplotRNA = [abundancedf[t, Symbol("R"*string(g).*gcn)] for t in 1:size(abundancedf)[1], gcn in gcnList]
##    toplotProt = [abundancedf[t, Symbol("P"*string(g).*gcn)] for t in 1:size(abundancedf)[1], gcn in gcnList]
##    gui(plot(abundancedf[:, :time], hcat(toplotRNA, toplotProt), show = true, reuse = false))

    toplotRNA = [abundancedf[t, Symbol("R"*string(g).*gcn)] for t in 1:size(abundancedf)[1], gcn in gcnList]
    if nod["coding"][g] == "PC"
        toplotProt = [abundancedf[t, Symbol("P"*string(g).*gcn)] for t in 1:size(abundancedf)[1], gcn in gcnList]
        gui(Plots.plot(abundancedf[:, :time], hcat(toplotRNA, toplotProt)[:, [1,5,2,6,3,7,4,8]], layout = @layout([a b]), show = true, reuse = false, title = "Gene "*string(g)))

    else
        gui(Plots.plot(abundancedf[:, :time], toplotRNA, show = true, reuse = false, title = "Gene "*string(g)))
    end

end

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

using JLD
using BioSimulator
include("winData/multiomics_networks_simulation/julia_functions.jl")

mechant = load("/home/oangelin/Documents/errmodel.jld")
QTLeffects = mechant["QTLeffects"]
InitVar = mechant["InitVar"]
stochmodel = mechant["stochmodel"]

gentil = load("/home/oangelin/Documents/goodmodel.jld")
QTLeffects = gentil["QTLeffects"]
InitVar = gentil["InitVar"]
stochmodel = gentil["stochmodel"]
nod = load("/home/oangelin/Documents/goodmodelnod.jld")["nod"]

model = BioSimulator.Network("coucouBad")

## Add the species in the model, with their initial abundance
for i in 1:length(stochmodel["species"])
    i0 = replace(stochmodel["initialconditions"][i], "QTLeffects", "$QTLeffects")
    i0 = replace(i0, "InitVar", "$InitVar")
    i0 = eval(parse(i0))
    #println(stochmodel["species"][i]* "\t"*string(i0))
    model <= BioSimulator.Species(stochmodel["species"][i], round(Int, i0))

    if !isa(i0, Number)
        println(stochmodel["initialconditions"][i])
        error("Pb i0")
    end
    if typeof(stochmodel["species"][i]) != String
        println(stochmodel["species"][i])
        error("Pb species")
    end
end

## Add the reactions in the model, with their name and rate
for i in eachindex(stochmodel["reactions"])
    prop = replace(stochmodel["propensities"][i], "QTLeffects", "$QTLeffects")
    #println(prop)
    prop = eval(parse(prop))
    reacname = replace(stochmodel["reactionsnames"][i], "de-", "de")
    #model <= BioSimulator.Reaction(stochmodel["reactionsnames"][i], prop, stochmodel["reactions"][i])
    model <= BioSimulator.Reaction(reacname, prop, stochmodel["reactions"][i])

    if !isa(prop, Number)
        println(stochmodel["propensities"][i])
        error("Pb prop")
    end
    if typeof(stochmodel["reactionsnames"][i]) != String
        println(stochmodel["reactionsnames"][i])
        error("Pb reactionsnames")
    end
    if typeof(stochmodel["reactions"][i]) != String
        println(stochmodel["reactions"][i])
        error("Pb reactions")
    end

end

#result = simulate(model, algorithm = SSA, time = 10.0, epochs = 1, trials = 1)

time = 10.0
epochs = 2
trials = 1
algorithm = SSA
D = Array
T = SSA

# extract model information
c = n_species(model)
d = n_reactions(model)

species = species_list(model)
reactions = reaction_list(model)

# create simulation data structures
x0, rxn, id, id2ind = BioSimulator.make_datastructs(species, reactions, c, d)
xt     = copy(x0)
alg    = algorithm(time)

output = BioSimulator.SimData(
    id2ind,
    linspace(0.0, time, epochs + 1),
    #SharedArray{eltype(x0)}(c, epochs + 1, trials),
    D{eltype(x0)}(c, epochs + 1, trials),
    alg.stats
)

BioSimulator.init!(alg, xt, rxn)

#simulate_wrapper!(output, xt, x0, alg, rxn)

data = output.data

i = 1
len = div(size(data, 3), Threads.nthreads())
domain = ((i-1)*len+1):i*len

#simulate_chunk!(output, xt, x0, alg, rxn, domain)
# here Xt = xt, X0 = x0, algorithm = alg, reactions = rxn, trial_set = domain

Xt = xt
X0 = x0

a = BioSimulator.propensities(rxn)

# for trial in trial_set
trial = 1

copy!(Xt, X0)

BioSimulator.update_all_propensities!(a, rxn, Xt)
BioSimulator.reset!(alg, a)

#simulate!(output, Xt, algorithm, reactions, trial)

epoch = 1

#props = zeros(Float64, length(a), 10000)
count = 0
tic(); while !done(alg)
#tic(); for count in 1:size(props, 2)
    t     = BioSimulator.get_time(alg)
    epoch = BioSimulator.update!(output, Xt, t, epoch, trial)

    BioSimulator.step!(alg, Xt, rxn)

#    props[:, count] = BioSimulator.propensities(rxn)

    count = count+1
    if count == 100000
        count = 0
        println(t)
        println(findmax(BioSimulator.propensities(rxn).cache))
    end

end; toc()

#BioSimulator.step!(alg, Xt, rxn)


text:
module ArgusWrapper

using CxxWrap
@wrapmodule(joinpath(@__DIR__, "..", "build", "lib", "libargus_julia_wrapper"), :define_argus_wrapper)

function __init__()
    @initcxx
end

end # module