text:
using OpenStreetMapX
using OpenStreetMapXPlot
using LightGraphs
using MetaGraphs
using Luxor
using Colors
import Plots

#=
visualization.jl

this script implements visualization alternatives for LPBasic problems and solutions

=#

#both of these work ok, but are quite slow and the osm information isn't really useful
export PlotInstance, PlotSolution

"""
    function PlotInstance(instance::ProblemInstance, osmpath::String)
Return a plot of the given instance using the given Open Street Map information.

Substantially slower than the PlotInstance method not using OSM information, but the result is better looking.
"""
function PlotInstance(instance::ProblemInstance, osmpath::String)
    map_data = OpenStreetMapX.get_map_data(osmpath)

    Plots.gr() #using the GR back-end
    p = OpenStreetMapXPlot.plotmap(map_data,width=600,height=400);

    # "convert" from our instance's network to a format frienldy to the OpenStreetMapX library

    for line in instance.lines

        #a route is simply an array of MapData nodes
        route = Array{Int64, 1}() #where T = Int64
        for stop_tuple in line.walk
            from, to = stop_tuple
            #latitude and longitude for to and from stops
            from_coord = (get_prop(instance.network, vertices(instance.network)[from], :latitude), get_prop(instance.network, vertices(instance.network)[from], :longitude))
            to_coord = (get_prop(instance.network, vertices(instance.network)[to], :latitude), get_prop(instance.network, vertices(instance.network)[to], :longitude))

            #the DataMap nodes
            from_node = OpenStreetMapXPlot.point_to_nodes(from_coord, map_data)
            to_node = OpenStreetMapXPlot.point_to_nodes(to_coord, map_data)

            #append to route
            push!(route, from_node)

        end

        #add route to plot
        addroute!(p,map_data,route,route_color="red",start_name="",end_name="");

    end

    return p
end

"""
    PlotSolution(instance::ProblemInstance, solution::Array{Bool, 1}, osmpath::String)
Identical to PlotInstance(instance::ProblemInstance, osmpath::String), but paint unused lines differently

Substantially slower than the PlotSolution method not using OSM information, but the result is better looking.
"""
function PlotSolution(instance::ProblemInstance, solution::Array{Bool, 1}, osmpath::String)
    map_data = OpenStreetMapX.get_map_data(osmpath)

    Plots.gr() #using the GR back-end
    p = OpenStreetMapXPlot.plotmap(map_data,width=600,height=400);

    # "convert" from our instance's network to a format frienldy to the OpenStreetMapX library

    line_i = 1
    for line in instance.lines

        #a route is simply an array of MapData nodes
        route = Array{Int64, 1}() #where T = Int64
        for stop_tuple in line.walk
            from, to = stop_tuple
            #latitude and longitude for to and from stops
            from_coord = (get_prop(instance.network, vertices(instance.network)[from], :latitude), get_prop(instance.network, vertices(instance.network)[from], :longitude))
            to_coord = (get_prop(instance.network, vertices(instance.network)[to], :latitude), get_prop(instance.network, vertices(instance.network)[to], :longitude))

            #the DataMap nodes
            from_node = OpenStreetMapXPlot.point_to_nodes(from_coord, map_data)
            to_node = OpenStreetMapXPlot.point_to_nodes(to_coord, map_data)

            #append to route
            push!(route, from_node)

        end

        #add route to plot
        if solution[line_i]
            addroute!(p,map_data,route,route_color="red",start_name="",end_name="")
        else
            addroute!(p,map_data,route,route_color="blue",start_name="",end_name="")
        end

        line_i +=1

    end

    return p

end

#visualization alternative using Luxor
"""
    PlotInstance(instance::ProblemInstance; out_path::String = "./", should_preview = true)
Return a simple Luxor drawing of the given instance

Substantially faster than the PlotInstance method using Open Street Map information, but draws only simple lines on a blank background.
"""
function PlotInstance(instance::ProblemInstance; out_path::String = "./", should_preview = true)
    PlotSolution(instance, [true for i in 1:length(instance.lines)], out_path = out_path, should_preview = should_preview)
end

"""
    PlotSolution(instance::ProblemInstance, solution::Array{Bool, 1};  <keyword arguments>)
Return a simple Luxor drawing of the given instance, painting lines used in solution differently

# Keywork arguments:
- [out_path]: optional path determining where to save the image. If this isn't set, a time-stamped name will be used
- should_preview=true: should the image be previewed in the default image viewer?
"""
function PlotSolution(instance::ProblemInstance, solution::Array{Bool, 1}; out_path::Union{String, Nothing} = nothing, should_preview = true)
    #min and max of latitudes and longitudes, used for scaling
    min_lat, max_lat = Inf, -Inf
    min_lon, max_lon = Inf, -Inf
    for i in 1:nv(instance.network)
        lat, lon = get_prop(instance.network, vertices(instance.network)[i], :latitude), get_prop(instance.network, vertices(instance.network)[i], :longitude)
        min_lat = min(min_lat, lat)
        max_lat = max(max_lat, lat)
        min_lon = min(min_lon, lon)
        max_lon = max(max_lon, lon)
    end

    #not sure if this is correct
    if max_lat - min_lat > max_lon - min_lon
        #latitude has bigger span than longitude
        slat = 600 / (max_lat - min_lat)
        slon = slat
    else
        #longitude has bigger span than latitude
        slon = 600 / (max_lon - min_lon)
        slat = slon
    end

    Drawing(600, 600, :png, if out_path != nothing out_path else "./$(Dates.minute(now())).png" end)
    origin()
     #default size is 600x600

        cx, cy = get_prop(instance.network, :buffer_center_lat), get_prop(instance.network, :buffer_center_lon)

        #random, different Colors
        # if seed != nothing
        #     dcolors = Colors.distinguishable_colors(length(instance.lines), seed)
        # else
        #     dcolors = Colors.distinguishable_colors(length(instance.lines))
        # end

        sf = 1
        println("slat: $(slat)  e slon:$(slon)")
        for i in 1:(length(instance.lines))
            line = instance.lines[i]
            Luxor.setcolor(if solution[i] "red" else "black" end) #red if included in solution, blue otherwise
            Luxor.setdash("solid")
            #Luxor.setcolor(if solution[i] dcolors[i] else "black" end)
            #Luxor.setdash(if solution[i] "solid" else "dashed" end)

            #Luxor.setcolor(dcolors[i])
            #Luxor.setopacity(if solution[i] 1.0 else 0.5 end)

            #if solution[i] action = :fillstroke else action = :dotted end
            #action = if (solution[i]) :fillstroke else :dotted end
            #Luxor.setcolor("black")
            for j in 1:(length(line.walk) - 1) #stop at -1 because we'll draw a line between every two stops
                from, to = line.walk[j]
                p1_x = slat * sf * (get_prop(instance.network, vertices(instance.network)[from], :latitude) - cx)
                p1_y = slon * -sf * (get_prop(instance.network, vertices(instance.network)[from], :longitude) -cy) #negated bc y axis points downwards

                p2_x = slat * sf * (get_prop(instance.network, vertices(instance.network)[to], :latitude) -cx)
                p2_y = slon * -sf * (get_prop(instance.network, vertices(instance.network)[to], :longitude) -cy)

                #might raise ERROR: cant draw arrow between two identical points. Maybe bc of imprecision?
                try
                    #Luxor.arrow(Luxor.Point(p1_x, p1_y), Luxor.Point(p2_x, p2_y))
                    Luxor.line(Luxor.Point(p1_x, p1_y), Luxor.Point(p2_x, p2_y), :stroke)
                catch e
                    @warn "ERROR: cant draw arrow between two identical points????"
                end
            end
        end

    finish()
    !should_preview || preview()
end

#to-do: plotting based on edge frequency, and not lines


text:
using Test
using SparseArrays
using ForwardDiff

@testset "GloVe" begin

    include("../src/GloVe.jl")
    test_corpus_path = "test_corpus.txt"
    test_cooc_path = "test_cooc.bin"
    embedding_path = "../embeddings/vectors-C0-V20-W8-D25-R0.05-E15-S1.bin"
    full_cooc_path = "../embeddings/cooc-C0-V20-W8.bin"

    @testset "Model loading" begin
        global M = GloVe.load_model(embedding_path)
        @test M.D == 25 # vector dim
        @test M.d == 8 # window
        # Test against some random examples from `head vocab.txt`
        @test M.ivocab[3] == "and"
        @test M.vocab["is"] == GloVe.WORD_INFO(6, 213234)
    end

    @testset "Vector saving" begin
        vec_path = "test_vecs.tmp"
        (V, D) = (100, 10)
        (W, b_w, U, b_u) = (rand(V, D), rand(V), rand(V, D), rand(V))
        GloVe.save_bin_vectors(vec_path, W, b_w, U, b_u)
        (W_1, b_w_1, U_1, b_u_1) = GloVe.load_bin_vectors(vec_path, V)
        @test isapprox(W, W_1)
        @test isapprox(b_w, b_w_1)
        @test isapprox(U, U_1)
        @test isapprox(b_u, b_u_1)
        rm(vec_path)
    end

    @testset "Cooc Matrix parsing/creation" begin
        X = GloVe.load_cooc(test_cooc_path, M.V)
        X_1 = GloVe.docs2cooc(eachline(test_corpus_path), M.vocab, M.d)
        @test isapprox(X, X_1)
        X_2 = spzeros(M.V, M.V)
        for line in eachline(test_corpus_path)
            X_2 += GloVe.doc2cooc(line, M.vocab, M.d)
        end
        @test isapprox(X, X_2)
        test_words = ["april", "the", "and", "having"]
        indices = [M.vocab[word].index for word in test_words]
        X_3 = GloVe.load_cooc(test_cooc_path, M.V, indices)
        @test X_3[indices, :] == X[indices, :]
        @test X_3[:, indices] == X[:, indices]
    end

    @testset "Cooc Matrix saving/loading" begin
        V = 100
        test_path = "test_cooc_2.tmp"
        X = sprand(V, V, 0.1)
        GloVe.save_coocs(test_path, X)
        X_1 = GloVe.load_cooc(test_path, V)
        @test isapprox(X, X_1)
        rm(test_path)
    end

    @testset "Cooc perturbation" begin
        V = 100
        test_path = "test_cooc_2.tmp"
        pert_path = "test_pert.tmp"

        X = 100 * sprand(V, V, 0.5)
        GloVe.save_coocs(test_path, X)
        δX = sprand(V, V, 0.05)

        GloVe.perturb_coocs(test_path, pert_path, -δX)
        X1 = GloVe.load_cooc(pert_path, V)
        Y = X - δX
        Y[Y .< 0] .= 0
        dropzeros(Y)

        @test isapprox(X1, Y)

        rm(pert_path)
        rm(test_path)
    end

    @testset "Gradients and Hessians" begin
        ∇ = ForwardDiff.gradient
        ∇² = ForwardDiff.hessian

        # Make some very toy examples
        (V, D) = (25, 4)
        (W, b_w, U, b_u) = (rand(V, D), rand(V), rand(V, D), rand(V))
        X = sprand(V, V, 0.5); X += X'

        ∇J = ∇(w->GloVe.J(w, b_w, U, b_u, X), W)
        ∇²J = ∇²(w->GloVe.J(w, b_w, U, b_u, X), W)

        # Convert ForwardDiff's indexing
        reindex(i) = V  * collect(0:D-1) .+ i

        gradients_match = true
        hessians_match = true

        for i in 1:V
            gradients_match &= isapprox(∇J[i,:], GloVe.∇Li(W, b_w, U, b_u, X, i))
            hessians_match &= isapprox(∇²J[reindex(i), reindex(i)], GloVe.∇²Li(U, X, i))
        end

        @test gradients_match
        @test hessians_match

    end

    @testset "Influence Function approximation" begin
        document = "the quick brown fox jumped over the lazy dog"
        X = GloVe.load_cooc(full_cooc_path, M.V)
        deltas = GloVe.compute_IF_deltas(document, M, X)
        # Check only words in document are affected
        @test Set(keys(deltas)) == Set([M.vocab[word].index for word in split(document)])
        all_non_zeros = true
        for (idx, delta) in pairs(deltas)
            all_non_zeros &= (length(delta) == M.D) && (delta != zeros(M.D))
        end
        @test all_non_zeros

        # Now test more efficient version
        target_words = split("jumped dog excluded")
        target_indices = [M.vocab[word].index for word in target_words]
        H = GloVe.inv_hessians_for(target_indices, M, X)
        G = GloVe.gradients_for(target_indices, M, X)

        select_deltas = GloVe.compute_IF_deltas(document, M, X, target_indices, H, G)
        @test length(select_deltas) == 2  # for "jumped" and "dog"

        same_values = true
        for (idx, delta) in pairs(select_deltas)
            same_values &= delta == deltas[idx]
        end
        @test same_values
    end

end

ext: jl


text:
@testset "Lattices" begin

  #
  # Constructors
  #

  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  f = x - 1;
  K, a = number_field(f, "a", cached = false)
  D = matrix(K, 3, 3, [4, 0, 0, 0, 10, 0, 0, 0, 20]);
  gens = Vector{nf_elem}[map(K, [0, 1, 0]), map(K, [0, 1, 0]), map(K, [-5//4, 3//2, 3//4]), map(K, [-5//4, 3//2, 3//4]), map(K, [-5//4, -1//2, -1//4]), map(K, [-5//4, -1//2, -1//4])]
  L = quadratic_lattice(K, gens, gram = D)

  L1 = @inferred quadratic_lattice(base_field(L), pseudo_matrix(L))
  @test pseudo_matrix(L1) == pseudo_matrix(L)
  @test ambient_space(L1) != ambient_space(L)

  L2 = @inferred quadratic_lattice(base_field(L), pseudo_matrix(L), gram = D)
  @test ambient_space(L2) === ambient_space(L)
  @test L == L2

  L3 = @inferred quadratic_lattice(base_field(L), matrix(pseudo_matrix(L)))
  @test matrix(pseudo_matrix(L3)) == matrix(pseudo_matrix(L))
  @test pseudo_matrix(L3) != pseudo_matrix(L)
  @test ambient_space(L3) != ambient_space(L)

  L4 = @inferred quadratic_lattice(base_field(L), generators(L))
  @test ambient_space(L4) != ambient_space(L)

  L5 = @inferred quadratic_lattice(base_field(L), generators(L), gram = D)
  @test ambient_space(L5) === ambient_space(L)

  #
  # Zero lattice
  #

  LL = @inferred lattice(ambient_space(L), [])
  @test ambient_space(LL) === ambient_space(L)
  @test rank(LL) == 0

  @test_throws AssertionError quadratic_lattice(base_field(L), [])

  LL = @inferred quadratic_lattice(base_field(L), [], gram = D)
  @test ambient_space(LL) === ambient_space(L)
  @test rank(LL) == 0

  D = matrix(K, 0, 0, [])
  gens = Vector{nf_elem}[]
  L = @inferred quadratic_lattice(K, gens, gram = D)
  @test isdefinite(L)
  @test rank(L) == 0

end


text:
#f(p,x)=KL(p,softmax(x))

FKLsoftmax(p::Array{Float64,2},x::Array{Float64,2})=begin DN=prod(size(x));Z=sum(exp(x),1);aux=(Z,DN); return ((sum(p.*(log(p)-x))+sum(log(Z)))/DN,aux); end

function FKLsoftmax_inplace(value,auxvalue,p::Array{Float64,2},x::Array{Float64,2})
    DN=prod(size(x));Z=sum(exp(x),1);aux=(Z,DN)
    value=((sum(p.*(log(p)-x))+sum(log(Z)))/DN,aux)
end

function DKLsoftmax(derivativeIDX,f_c,faux_c,grad_c,grad_n,x::Array{Float64,2})
    if derivativeIDX==1
        grad_c.*(1+log(p)-x)
    elseif derivativeIDX==2
        grad_c.*(exp(x)./aux[1]-p)./aux[2]
    end
end

if PROC=="GPU"
end

Derivative[FKLsoftmax]=DKLsoftmax
Inplace[FKLsoftmax]=FKLsoftmax_inplace

KLsoftmax(p::ADnode,x::ADnode)=ADnode(FKLsoftmax,[p x])

export KLsoftmax

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
"""
docstring_ext

Helper function to generate the docstring for a package.
"""
function docstring_ext(T; descr::String="")
    package_name = MLJBase.package_name(T)
    package_url  = MLJBase.package_url(T)
    model_name   = MLJBase.name(T)
    # the message to return
    message      = "$descr"
    message     *= "\n→ based on [$package_name]($package_url)."
    message     *= "\n→ do `@load $model_name pkg=\"$package_name\"` to use the model."
    message     *= "\n→ do `?$model_name` for documentation."
end

"""
metadata_pkg

Helper function to write the metadata for a package.
"""
function metadata_pkg(T; name::String="unknown", uuid::String="unknown",
                      url::String="unknown", julia::Union{Missing,Bool}=missing,
                      license::String="unknown", is_wrapper::Bool=false)
    ex = quote
        MLJBase.package_name(::Type{<:$T}) = $name
        MLJBase.package_uuid(::Type{<:$T}) = $uuid
        MLJBase.package_url(::Type{<:$T}) = $url
        MLJBase.is_pure_julia(::Type{<:$T}) = $julia
        MLJBase.package_license(::Type{<:$T}) = $license
        MLJBase.is_wrapper(::Type{<:$T}) = $is_wrapper
    end
    parentmodule(T).eval(ex)
end

"""
metadata_model

Helper function to write the metadata for a single model of a package (complements
[`metadata_ext`](@ref)).
"""
function metadata_model(T; input=Unknown, target=Unknown,
                        output=Unknown, weights::Bool=false,
                        descr::String="", path::String="")
    if isempty(path)
        path = "MLJModels.$(MLJBase.package_name(T))_.$(MLJBase.name(T))"
    end
    ex = quote
        MLJBase.input_scitype(::Type{<:$T}) = $input
        MLJBase.output_scitype(::Type{<:$T}) = $output
        MLJBase.target_scitype(::Type{<:$T}) = $target
        MLJBase.supports_weights(::Type{<:$T}) = $weights
        MLJBase.docstring(::Type{<:$T}) = MLJBase.docstring_ext($T; descr=$descr)
        MLJBase.load_path(::Type{<:$T}) = $path
    end
    parentmodule(T).eval(ex)
end

"""
metadata_measure

Helper function to write the metadata for a single measure.
"""
function metadata_measure(T; name::String="",
                          target_scitype=Unknown,
                          prediction_type::Symbol=:unknown,
                          orientation::Symbol=:unknown,
                          reports_each_observation::Bool=true,
                          aggregation=Mean(),
                          is_feature_dependent::Bool=false,
                          supports_weights::Bool=false,
                          docstring::String="")
    pred_str = "$prediction_type"
    orientation_str = "$orientation"
    ex = quote
        isempty($name) || (MLJBase.name(::Type{<:$T}) = $name)
        isempty($docstring) || (MLJBase.docstring(::Type{<:$T}) = $docstring)
        MLJBase.target_scitype(::Type{<:$T}) = $target_scitype
        MLJBase.prediction_type(::Type{<:$T}) = Symbol($pred_str)
        MLJBase.orientation(::Type{<:$T}) = Symbol($orientation_str)
        MLJBase.reports_each_observation(::Type{<:$T}) = $reports_each_observation
        MLJBase.aggregation(::Type{<:$T}) = $aggregation
        MLJBase.is_feature_dependent(::Type{<:$T}) = $is_feature_dependent
        MLJBase.supports_weights(::Type{<:$T}) = $supports_weights
    end
    parentmodule(T).eval(ex)
end


text:
# RELOG: Reverse Logistics Optimization
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using DataStructures
using JSON
using JSONSchema
using Printf
using Statistics

"""
    _compress(instance::Instance)

Create a single-period instance from a multi-period one. Specifically,
replaces every time-dependent attribute, such as initial_amounts,
by a list with a single element, which is either a sum, an average,
or something else that makes sense to that specific attribute.
"""
function _compress(instance::Instance)::Instance
    T = instance.time
    compressed = deepcopy(instance)
    compressed.time = 1
    compressed.building_period = [1]

    # Compress products
    for p in compressed.products
        p.transportation_cost = [mean(p.transportation_cost)]
        p.transportation_energy = [mean(p.transportation_energy)]
        for (emission_name, emission_value) in p.transportation_emissions
            p.transportation_emissions[emission_name] = [mean(emission_value)]
        end
    end

    # Compress collection centers
    for c in compressed.collection_centers
        c.amount = [maximum(c.amount) * T]
    end

    # Compress plants
    for plant in compressed.plants
        plant.energy = [mean(plant.energy)]
        for (emission_name, emission_value) in plant.emissions
            plant.emissions[emission_name] = [mean(emission_value)]
        end
        for s in plant.sizes
            s.capacity *= T
            s.variable_operating_cost = [mean(s.variable_operating_cost)]
            s.opening_cost = [s.opening_cost[1]]
            s.fixed_operating_cost = [sum(s.fixed_operating_cost)]
        end
        for (prod_name, disp_limit) in plant.disposal_limit
            plant.disposal_limit[prod_name] = [sum(disp_limit)]
        end
        for (prod_name, disp_cost) in plant.disposal_cost
            plant.disposal_cost[prod_name] = [mean(disp_cost)]
        end
    end

    return compressed
end

max_stars_repo_name: ANL-CEEESA/RELOG

max_issues_repo_name: ANL-CEEESA/RELOG

max_forks_repo_name: ANL-CEEESA/RELOG


text:
export PARSDMM

"""
Computes the projection onto an intersection of sets
min_x 1/2||x-m||2^2 s.t. x in C_1(A_1*x), x in C_2(A_2*x), ... , x in C_p(A_p*x)

input:

#Arguments
- m                   - vector to be projected
- TD_OP               - vector of linear operators [A_1;A_2;...;A_p]
- set_Prop            - set properties, structure where each property is vector which has a lenght of the number of sets
                            (see SetIntersectionProjection.jl and setup_constraints.jl)
- comp_grid           - structure with grid info. comp_grid.d = (dx,dy,dz), comp_grid.n = (nx,ny,nz)
- P_sub               - vector of projection functions onto C: P_sub = [P_C_1(.); P_C_2(.); ... ; P_C_p(.)]
- options             - structure with options, see SetIntersectionProjection.jl
- x                   - output vector
- l                   - vector of vector of Lagrangian multipliers l = [l_1;l_2,...;l_p] (optional)
- y                   - vector of vector of auxiliary vectors y = [y_1;y_2,...;y_p] (optional)

output:
 - x           -   result
 - log_PARSDMM -   constains log information about various quantities per iteration
"""
function PARSDMM(m         ::Vector{TF},
                 AtA       ::Union{Vector{Array{TF, 2}}, Vector{SparseMatrixCSC{TF, TI}}, Vector{joAbstractLinearOperator{TF, TF}}, Vector{Union{Array{TF, 2}, SparseMatrixCSC{TF, TI}, joAbstractLinearOperator{TF, TF}}}},
                 TD_OP     ::Union{Vector{Union{SparseMatrixCSC{TF,TI},joAbstractLinearOperator{TF,TF}}},DistributedArrays.DArray{Union{JOLI.joAbstractLinearOperator{TF,TF}, SparseMatrixCSC{TF,TI}},1,Array{Union{joAbstractLinearOperator{TF,TF}, SparseMatrixCSC{TF,TI}},1}} },
                 set_Prop,
                 P_sub     ::Union{Vector{Any},DistributedArrays.DArray{Any,1,Array{Any,1}}},
                 comp_grid,
                 options,
                 x=zeros(TF,length(m)) ::Vector{TF},
                 l=[],
                 y=[]
                 ) where {TF<:Real,TI<:Integer}

# Create a TimerOutput, this is the main type that keeps track of everything.
to = TimerOutput()

@timeit to "initialization" begin

# Parse default options
convert_options!(options,TF)
@unpack  x_min_solver,maxit,evol_rel_tol,feas_tol,obj_tol,rho_ini,rho_update_frequency,gamma_ini,
adjust_rho,adjust_gamma,adjust_feasibility_rho,Blas_active,
feasibility_only,FL,parallel,zero_ini_guess = options

# Input checks
# PARSDMM can work with imaginary numbers, but we want to keep it real
if isreal(m)==false || isreal(x)==false || isreal(l)==false || isreal(y)==false
    error("input for PARSDMM is not real")
end

# initialize
pp = length(TD_OP);
feasibility_only == false && (pp=pp-1)

(ind_ref,N,TD_OP,AtA,p,rho_update_frequency,adjust_gamma,adjust_rho,adjust_feasibility_rho,gamma_ini,rho,gamma,y,y_0,y_old,l,l_0,l_old,
l_hat_0,x_0,x_old,r_dual,rhs,s,s_0,Q,prox,log_PARSDMM,l_hat,x_hat,r_pri,d_l_hat,d_H_hat,d_l,
d_G_hat,P_sub,Q_offsets,stop,feasibility_initial,set_feas,Ax_out)=PARSDMM_initialize(x,l,y,AtA,TD_OP,set_Prop,P_sub,comp_grid,maxit,rho_ini,gamma_ini,
x_min_solver,rho_update_frequency,adjust_gamma,adjust_rho,adjust_feasibility_rho,m,parallel,options,zero_ini_guess,feasibility_only)

if stop==true #stop if feasibility of input is detected by PARSDMM_initialize
  copy!(x,m)
  if options.Minkowski == true
    if length(x)==length(m)
      x = [x ; zeros(TF,length(x)) ]
    end
  end
  log_PARSDMM.obj             = log_PARSDMM.obj[[1]]
  log_PARSDMM.evol_x          = log_PARSDMM.evol_x[[1]]
  log_PARSDMM.r_pri_total     = log_PARSDMM.r_pri_total[[1]]
  log_PARSDMM.r_dual_total    = log_PARSDMM.r_dual_total[[1]]
  log_PARSDMM.r_pri           = log_PARSDMM.r_pri[[1],:]
  log_PARSDMM.r_dual          = log_PARSDMM.r_dual[[1],:]
  log_PARSDMM.cg_it           = log_PARSDMM.cg_it[[1]]
  log_PARSDMM.cg_relres       = log_PARSDMM.cg_relres[[1]]
  log_PARSDMM.set_feasibility = log_PARSDMM.set_feasibility[[1],:]
  log_PARSDMM.gamma           = log_PARSDMM.gamma[[1],:]
  log_PARSDMM.rho             = log_PARSDMM.rho[[1],:]
  return x, log_PARSDMM, l, y
end

#make sure proper vectors are supplied as initial guess for Minkowski sets.
if options.Minkowski == true
  if length(x)==length(m)
    x = [x ; zeros(TF,length(m)) ]
  end
end

counter = 2

x_solve_tol_ref = TF(1.0) #scalar required to determine tolerance for x-minimization, initial value doesn't matter

end #end initialization timer

for i=1:maxit #main loop

  #form right hand side for x-minimization
  @timeit to "form rhs for linear system" begin
  rhs = rhs_compose(rhs,l,y,rho,TD_OP,p,Blas_active,parallel)
  end #end timer for rhs

  # x-minimization
  @timeit to "argmin x" begin
  copy!(x_old,x);
  (x,iter,relres,x_solve_tol_ref) = argmin_x(Q,rhs,x,x_solve_tol_ref,i,log_PARSDMM,Q_offsets,Ax_out,comp_grid)
  log_PARSDMM.cg_it[i]            = iter
  log_PARSDMM.cg_relres[i]        = relres
  end #end timer for argmin x

  # y-minimization & l-update
  @timeit to "argmin y and l update" begin
  if parallel==true

    #using spawnat
      [ @sync @spawnat pid update_y_l_parallel(x,i,Blas_active,
        y[:L],y_old[:L],l[:L],l_old[:L],rho[:L],gamma[:L],prox[:L],TD_OP[:L],P_sub[:L],
        x_hat[:L],r_pri[:L],s[:L],set_feas[:L],feasibility_only) for pid in y.pids]

    #logging distributed quantities
    for ii=1:p
      log_PARSDMM.r_pri[i,ii] = norm(r_pri[ii])
      #log_PARSDMM.r_dual[i,ii] = norm( rho[ii]*(TD_OP[ii]'*(y[ii]-y_old[ii])) )
    end
    if mod(i,10)==0 #log every 10 it, or whatever number is suitable
      for ii=1:pp
        log_PARSDMM.set_feasibility[counter,ii] = set_feas[ii]
      end
      counter+=1
    end
  else #not distributed
    (y,l,r_pri,s,log_PARSDMM,counter,y_old,l_old) = update_y_l(x,p,i,Blas_active,y,y_old,l,l_old,rho,gamma,prox,TD_OP,log_PARSDMM,P_sub,counter,x_hat,r_pri,s,feasibility_only);
    log_PARSDMM.r_dual_total[i] = sum(log_PARSDMM.r_dual[i,:]);
  end

  #some more logging
  log_PARSDMM.r_pri_total[i] = sum(log_PARSDMM.r_pri[i,:]);
  if options.Minkowski == false
    log_PARSDMM.obj[i] = TF(0.5) .* norm(x .- m)^2
  else
    log_PARSDMM.obj[i] = TF(0.5) .* norm(TD_OP[end]*x .- m)^2
  end

  log_PARSDMM.evol_x[i]  = norm(x_old .- x) ./ norm(x);
  log_PARSDMM.rho[i,:]   = rho;
  log_PARSDMM.gamma[i,:] = gamma;

  end #end timer for argmin y and l-update

  # Stopping conditions
  @timeit to "stopping conditions check" begin
  (stop,adjust_rho,adjust_gamma,adjust_feasibility_rho,ind_ref) = stop_PARSDMM(log_PARSDMM,i,evol_rel_tol,feas_tol,obj_tol,adjust_rho,adjust_gamma,adjust_feasibility_rho,ind_ref,counter);
  if stop==true
    (TD_OP,AtA,log_PARSDMM) = output_check_PARSDMM(x,TD_OP,AtA,log_PARSDMM,i,counter)
    log_PARSDMM.timing = to
    return x, log_PARSDMM, l, y
  end

  end #end timer for stopping conditions

  # adjust penalty parameters rho and relaxation parameters gamma
  @timeit to "adjust rho and gamma" begin
  if i==1 #only at the first iteration
    if parallel==true
      [@spawnat pid l_hat[:L][1] = l_old[:L][1] .+ rho[:L][1] .* ( -s[:L][1] .+ y_old[:L][1] ) for pid in l_hat.pids]
      [@spawnat pid copy!(l_hat_0[:L][1],l_hat[:L][1] ) for pid in l_hat.pids]
      [@spawnat pid copy!(y_0[:L][1] , y[:L][1] ) for pid in y.pids]
      [@spawnat pid copy!(s_0[:L][1] , s[:L][1] ) for pid in s.pids]
      [@spawnat pid copy!(l_0[:L][1] , l[:L][1] ) for pid in l.pids]
    else
      for ii = 1:p
        l_hat[ii]   = l_old[ii] .+ rho[ii] .* ( -s[ii] .+ y_old[ii] );
        copy!(l_hat_0[ii],l_hat[ii])
        copy!(y_0[ii],y[ii])
        copy!(s_0[ii],s[ii])
        copy!(l_0[ii],l[ii])
      end
    end
  end

  if (adjust_rho == true || adjust_gamma == true) && mod(i,rho_update_frequency)==0#TF(0)
    if parallel==true
      #using spawnat
      [@sync @spawnat pid adapt_rho_gamma_parallel(gamma[:L],rho[:L],adjust_gamma,adjust_rho,y[:L],
      y_old[:L],s[:L],s_0[:L],l[:L],l_hat_0[:L],l_0[:L],l_old[:L],y_0[:L],l_hat[:L],d_l_hat[:L],d_H_hat[:L],d_l[:L],d_G_hat[:L]) for pid in y.pids]

    else
      (rho,gamma,l_hat,d_l_hat,d_H_hat,d_l,d_G_hat) = adapt_rho_gamma(i,gamma,rho,adjust_gamma,adjust_rho,y,y_old,s,s_0,l,l_hat_0,l_0,l_old,y_0,p,l_hat,d_l_hat,d_H_hat,d_l,d_G_hat);
    end

    if i>1 #after the first iteration
      if parallel==true
        [@spawnat pid copy!(l_hat_0[:L][1],l_hat[:L][1] ) for pid in l_hat.pids]
        [@spawnat pid copy!(y_0[:L][1] , y[:L][1] ) for pid in y.pids]
        [@spawnat pid copy!(s_0[:L][1] , s[:L][1] ) for pid in s.pids]
        [@spawnat pid copy!(l_0[:L][1] , l[:L][1] ) for pid in l.pids]
      else
        for ii=1:p
          copy!(l_hat_0[ii],l_hat[ii])
          copy!(y_0[ii],y[ii])
          copy!(s_0[ii],s[ii])
          copy!(l_0[ii],l[ii])
        end
      end #end if parallel
    end
  end #end adjust rho and gamma

  #adjust rho to set-feasibility estimates
  if parallel==true
    rho = convert(Vector{TF},rho); #gather rho
  end
  if adjust_feasibility_rho == true && mod(i,10)==TF(0) #&& norm(rho-log_PARSDMM.rho[i,:])<(10*eps(TF))#only update rho if it is not already updated
    ## adjust rho feasibility
    #if primal residual for a set is much larger than for the other sets
    #and the feasibility error is also much larger, increase rho to lower
    #primal residual and (hopefully) feasibility error
    (max_set_feas,max_set_feas_ind) = findmax(log_PARSDMM.set_feasibility[counter-1,:])
    sort_feas = sort(log_PARSDMM.set_feasibility[counter-1,:]);
    if i>10
      rho[max_set_feas_ind] = TF(2.0) .* rho[max_set_feas_ind]
    end
  end #end adjust_feasibility_rho

  #enforce max and min values for rho, to prevent the condition number of Q -> inf
  rho = max.(min.(rho,TF(1e4)),TF(1e-2)) #hardcoded bounds
  end #end timer for rho and gamma update

  @timeit to "Q-update" begin
  ind_updated = findall(rho .!= log_PARSDMM.rho[i,:]) # locate changed rho index
  ind_updated = convert(Array{TI,1},ind_updated)

  #re-assemble total transform domain operator as a matrix
  if isempty(findall((in)(ind_updated),p))==false
    if parallel==true && feasibility_only==false
      prox    = convert(Vector{Any},prox); #gather rho
      prox[p] = input -> prox_l2s!(input,rho[p],m)
      prox    = distribute(prox)
    elseif feasibility_only==false
      prox[p] = input -> prox_l2s!(input,rho[p],m)
    end
  end
  Q = Q_update!(Q,AtA,set_Prop,rho,ind_updated,log_PARSDMM,i,Q_offsets)
  if parallel==true
    rho = distribute(rho) #distribute again (gather -> process -> distribute is a ugly hack, fix this later)
  end
end #end Q-update timer

  if i==maxit
    println("PARSDMM reached maxit")
    (TD_OP,AtA,log_PARSDMM) = output_check_PARSDMM(x,TD_OP,AtA,log_PARSDMM,i,counter)
  end

 end #end main loop

log_PARSDMM.timing = to
return x , log_PARSDMM, l , y
end #end funtion

# output checks
function output_check_PARSDMM(x,TD_OP,AtA,log_PARSDMM,i,counter)
  if isreal(x)==false
      @warn "Result of PARSDMM is not real"
  end
  log_PARSDMM.obj             = log_PARSDMM.obj[1:i]
  log_PARSDMM.evol_x          = log_PARSDMM.evol_x[1:i]
  log_PARSDMM.r_pri_total     = log_PARSDMM.r_pri_total[1:i]
  log_PARSDMM.r_dual_total    = log_PARSDMM.r_dual_total[1:i]
  log_PARSDMM.r_pri           = log_PARSDMM.r_pri[1:i,:]
  log_PARSDMM.r_dual          = log_PARSDMM.r_dual[1:i,:]
  log_PARSDMM.cg_it           = log_PARSDMM.cg_it[1:i]
  log_PARSDMM.cg_relres       = log_PARSDMM.cg_relres[1:i]
  log_PARSDMM.set_feasibility = log_PARSDMM.set_feasibility[1:counter,:]
  log_PARSDMM.gamma           = log_PARSDMM.gamma[1:i,:]
  log_PARSDMM.rho             = log_PARSDMM.rho[1:i,:]

  return TD_OP , AtA , log_PARSDMM
end #end output_check_PARSDMM

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
module CBindingGen
	import Clang_jll

	const LIBCLANG_PATH = Clang_jll.libclang_path
	const LIBCLANG_VERSION = let
		dir = joinpath(dirname(dirname(LIBCLANG_PATH)), "lib", "clang")
		entries = readdir(dir)
		ind = findfirst(x -> !isnothing(match(r"^\d+\.\d+\.\d+$", x)) && isdir(joinpath(dir, x)), entries)
		isnothing(ind) && error("Failed to determine the version of libclang")
		entries[ind]
	end

	baremodule LibClang
		using CBinding: 𝐣𝐥
		using ..CBindingGen: LIBCLANG_PATH, LIBCLANG_VERSION

		const size_t = 𝐣𝐥.Csize_t

		𝐣𝐥.Base.include((𝐣𝐥.@__MODULE__), 𝐣𝐥.joinpath(𝐣𝐥.dirname(𝐣𝐥.@__DIR__), "deps", "libclang-$(LIBCLANG_VERSION).jl"))
	end

	using CBinding
	import Markdown

	export LibClang
	export Converted, CodeLocation, CodeRange
	export convert_header, convert_headers, generate

	function _string(func::Function, args...)
		cxstr = nothing
		try
			cxstr = func(args...)
			ptr = LibClang.clang_getCString(cxstr)
			return ptr == C_NULL ? nothing : unsafe_string(ptr)
		finally
			isnothing(cxstr) || LibClang.clang_disposeString(cxstr)
		end
	end

	Base.string(cursor::LibClang.CXCursor) = _string(LibClang.clang_getCursorSpelling, cursor)
	Base.string(tu::LibClang.CXTranslationUnit, token::LibClang.CXToken) = _string(LibClang.clang_getTokenSpelling, tu, token)

	function children(cursor::LibClang.CXCursor)
		Cvisitor = @ccallback function visitor(c::LibClang.CXCursor, p::LibClang.CXCursor, cs::LibClang.CXClientData)::LibClang.CXChildVisitResult
			push!(unsafe_pointer_to_objref(cs)::Vector{LibClang.CXCursor}, c)
			return LibClang.CXChildVisit_Continue
		end

		cursors = LibClang.CXCursor[]
		Bool(LibClang.clang_visitChildren(cursor, Cvisitor, pointer_from_objref(cursors))) && error("Failed to visit cursor children")
		return cursors
	end

	struct CodeLocation
		file::String
		line::Int
		col::Int

		function CodeLocation(loc::LibClang.CXSourceLocation)
			cxfile = Ref{LibClang.CXFile}(C_NULL)
			line = Ref{Cuint}(0)
			col = Ref{Cuint}(0)
			offset = Ref{Cuint}(0)
			LibClang.clang_getSpellingLocation(loc, cxfile, line, col, offset)

			str = cxfile[] == C_NULL ? nothing : _string(LibClang.clang_getFileName, cxfile[])
			file = isnothing(str) ? "<unknown>" : str

			return new(file, line[], col[])
		end
	end
	CodeLocation(cursor::LibClang.CXCursor) = CodeLocation(LibClang.clang_getCursorLocation(cursor))

	Base.isless(a::CodeLocation, b::CodeLocation) = a.file == b.file && (a.line < b.line || (a.line == b.line && a.col < b.col))

	function Base.string(cl::CodeLocation; relto::String = ".")
		file = !startswith(relpath(cl.file, relto), "../") ? "./$(relpath(cl.file, relto))" : !startswith(relpath(cl.file, homedir()), "../") ? "~/$(relpath(cl.file, homedir()))" : cl.file
		return "$(file):$(cl.line):$(cl.col)"
	end

	struct CodeRange <: AbstractRange{CodeLocation}
		start::CodeLocation
		stop::CodeLocation

		function CodeRange(range::LibClang.CXSourceRange)
			start = CodeLocation(LibClang.clang_getRangeStart(range))
			stop  = CodeLocation(LibClang.clang_getRangeEnd(range))
			return new(start, stop)
		end
	end
	function CodeRange(cursor::LibClang.CXCursor)
		cursor.kind == LibClang.CXCursor_NoDeclFound && error("Cannot create a CodeRange for a $(cursor.kind) cursor")
		return CodeRange(LibClang.clang_getCursorExtent(cursor))
	end

	function Base.coalesce(cursors::Vector{LibClang.CXCursor})
		# TODO: switch to OrderedDict and use in last part of this function so its not O(N^2)
		coalesced = Vector{Pair{CodeRange, LibClang.CXCursor}}[]
		startInds = Dict{CodeLocation, Int}()
		for cursor in cursors
			range = CodeRange(cursor)

			if !haskey(startInds, range.start)
				push!(coalesced, [range => cursor])
				startInds[range.start] = length(coalesced)
			else
				# NOTE: insertion sorting the coalesced group so we know where the largest range in the group is
				coal = coalesced[startInds[range.start]]
				ind = findfirst(c -> range in c.first, coal)
				ind = isnothing(ind) ? length(coal) : ind
				insert!(coal, ind+1, range => cursor)
			end
		end

		# once the cursors are coalesced, we need to remove groups that are nested within other groups
		coalesced = filter(inners -> !any(outers -> last(inners) !== last(outers) && last(inners).first in last(outers).first, coalesced), coalesced)
		return map(cs -> map(c -> c.second, cs), coalesced)
	end

	Base.in(inner::LibClang.CXCursor, outer::LibClang.CXCursor) = in(CodeRange(inner), CodeRange(outer))
	function Base.in(inner::CodeRange, outer::CodeRange)
		inner.start.file == outer.start.file || return false
		inner.stop.file == outer.stop.file || return false

		!(inner.start < outer.start) && !(outer.stop < inner.stop)
	end

	struct Comment
		md::Union{Markdown.MD, Nothing}
		locs::Vector{CodeLocation}
	end

	function Base.string(comment::Comment; expr::String = "", relto::String = ".")
		contents = []
		isempty(expr) || push!(contents, Markdown.Code(expr))
		isnothing(comment.md) || append!(contents, comment.md.content)
		push!(contents, Markdown.Header(length(comment.locs) > 1 ? "References" : "Reference", 2))
		for loc in comment.locs
			push!(contents, Markdown.Paragraph(Markdown.Link("$(basename(loc.file)):$(loc.line)", string(loc, relto = relto))))
		end
		return string(Markdown.MD(contents))
	end

	struct Converted
		expr::String
		comments::Dict{String, Comment}
	end

	# NOTE: func returns:  false = don't convert, true = convert as `string(cursor)`, "jlname" = convert as "jlname"
	convert_header(func::Function, header::String; kwargs...) = convert_headers(func, [header]; kwargs...)
	function convert_headers(func::Function, headers::Vector{String}; args::Vector{String} = String[])
		args = vcat(args, [
			"-isystem", joinpath(dirname(dirname(LIBCLANG_PATH)), "lib", "clang", LIBCLANG_VERSION, "include"),
		])

		return mktempdir() do dir
			hdr = joinpath(dir, "parse_headers.h")
			open(hdr, "w+") do file
				for header in headers
					println(file, "#include <$(header)>")
				end
			end

			ind = nothing
			try
				ind = LibClang.clang_createIndex(0, 1)
				ind == C_NULL && error("Failed to create index")

				tu = Ref(LibClang.CXTranslationUnit(C_NULL))
				try
					argsPtrs = map(i -> pointer(args[i]), eachindex(args))
					err = LibClang.clang_parseTranslationUnit2(ind, hdr, argsPtrs, length(args), C_NULL, 0, LibClang.CXTranslationUnit_DetailedPreprocessingRecord, tu)
					err == LibClang.CXError_Success || error("Failed to parse translation unit $(err)")

					for i in 1:LibClang.clang_getNumDiagnostics(tu[])
						diag = nothing
						try
							diag = LibClang.clang_getDiagnostic(tu[], i-1)
							sev = LibClang.clang_getDiagnosticSeverity(diag)
							sev in (LibClang.CXDiagnostic_Error, LibClang.CXDiagnostic_Fatal) && error("Errors encountered parsing headers, unable to proceed with conversion")
						finally
							LibClang.clang_disposeDiagnostic(diag)
						end
					end

					return convert_parsed(func, tu[])
				finally
					tu[] == C_NULL || LibClang.clang_disposeTranslationUnit(tu[])
				end
			finally
				isnothing(ind) || LibClang.clang_disposeIndex(ind)
			end
		end
	end

	function convert_parsed(func::Function, tu::LibClang.CXTranslationUnit)
		root = LibClang.clang_getTranslationUnitCursor(tu)

		# remove compiler-internal cursors and pre-processor cursors
		cursors = filter(c -> CodeLocation(c).file != "<unknown>" && !(c.kind in (
			LibClang.CXCursor_InclusionDirective,
			LibClang.CXCursor_MacroDefinition,
			LibClang.CXCursor_MacroInstantiation,
		)), children(root))

		cvts = Converted[]
		for coal in coalesce(cursors)
			try
				coal[1].kind in (
					LibClang.CXCursor_TypedefDecl,
					LibClang.CXCursor_EnumDecl,
					LibClang.CXCursor_StructDecl,
					LibClang.CXCursor_UnionDecl,
				) && func(coal[1]) === false && continue

				cvt = convert_fields(coal, 0, func)
				isnothing(cvt) || push!(cvts, cvt)
			catch
				@warn "The following exception occurred when converting near $(string(CodeLocation(coal[1])))"
				rethrow()
			end
		end

		exports = mapreduce(cvt -> collect(keys(cvt.comments)), vcat, cvts, init = String[])
		cursors = filter(c -> CodeLocation(c).file != "<unknown>" && c.kind == LibClang.CXCursor_MacroDefinition, children(root))

		macros = Converted[]
		names = Dict{String, Int}()
		for cursor in cursors
			try
				func(cursor) || continue

				cvt = convert_macro(tu, cursor, exports)
				if !isnothing(cvt)
					(name, m) = cvt
					if haskey(names, name)
						macros[names[name]] = m
					else
						push!(macros, m)
						names[name] = length(macros)
					end
				end
			catch
				@warn "The following exception occurred when converting macro near $(string(CodeLocation(cursor)))"
				rethrow()
			end
		end
		append!(cvts, macros)

		return cvts
	end

	include("names.jl")
	include("macros.jl")
	include("types.jl")
	include("fields.jl")
	include("functions.jl")
	include("comments.jl")
	include("generate.jl")
end


text:
using EDFPlus
using Test

edfh = loadfile("EDFPlusTestFile.edf")

@test edfh.gender == "Female"
@test edfh.annotationchannel == 30
@test edfh.filetype == FileStatus(1)

fs = samplerate(edfh,1)
f1 = EDFPlus.recordslice(edfh, 21, 22)[1,:]
f2 = highpassfilter(reshape(f1, length(f1)), fs)
f3 = lowpassfilter(f2, fs)
f4 = notchfilter(f3, fs)
@test round(f4[end-3], digits=2) == 6803.28
@test EDFPlus.recordindexat(edfh, edfh.file_duration - 0.05) == edfh.datarecords
@test EDFPlus.epochmarkers(edfh, 10.525)[7] == (632, 10)

eegpages = epoch_iterator(edfh, 12, channels=[7,8,9])
annots = annotation_epoch_iterator(edfh, 12)

for (pgnum, page) in enumerate(eegpages)
    if pgnum == 1
        @test round(page[3][end], digits=3) == -24.121
    elseif pgnum == 30
        @test round(page[1][100], digits=3) == -87.012
    end
end

for (i, ann) in enumerate(annots)
    if i == 9
        @test ann[1][1].onset == 95.9
    end
end

closefile!(edfh)
@test edfh.filetype == EDFPlus.CLOSED

bdfh = loadfile("samplefrombiosemicom.bdf")

@test "$(bdfh.filetype)" == "BDF" && bdfh.filetype == EDFPlus.BDF
@test trim(bdfh.signalparam[4].label) == "A4"
@test bdfh.startdate_year == 2001
@test bdfh.signalparam[end].physdimension == "Boolean"
@test EDFPlus.recordslice(bdfh, 4, 14)[1,end-3] == 1835263

@test EDFPlus.version() == 0.59
sig = physicalchanneldata(bdfh, 1)
@test round(sig[100], digits=3) == 523818.031
@test EDFPlus.recordindexat(bdfh, bdfh.file_duration - 0.3) == bdfh.datarecords

true


text:
function xl(exportall::Bool = false)
    update!(CACHE["config"])

    files = xlsxfilenames(CACHE["config"])
    if isempty(files)
        print_section("nothing to export."; color=:yellow)
    else
        print_section(
            "exporting xlsx files... ⚒\n" * "-"^(displaysize(stdout)[2] - 4);
            color = :cyan,
        )
        for f in files
            try
                export_xlsxtable(f)
            catch e
                printstyled("$f export failed\n"; color = :red)
            end
        end
        print_section("$(length(files)) xlsx files are exported ☺", "DONE"; color = :cyan)
    end
end
function xl(fname)
    update!(CACHE["config"])

    print_section(
        "exporting xlsx file... ⚒\n" * "-"^(displaysize(stdout)[2] - 4);
        color = :cyan,
    )
    export_xlsxtable(fname)
    print_section("export complete ☺", "DONE"; color = :cyan)

    nothing
end

"""
    export_xlsxtable(file::AbstractString)

export a given excel worsheets to a output file specified in 'config.json'
"""
function export_xlsxtable(fname)
    println("『", fname, "』")
    tb = loadtable(fname)

    for s in sheetnames(tb)
        fname = tb.out[s]
        write_worksheet(fname, tb.data[s])
        localizedata = tb.localizedata[s]
        if !ismissing(localizedata)
            write_localize(fname, localizedata)
        end
    end
    nothing
end
function write_worksheet(fname, jws::JSONWorksheet)
    dir = GAMEENV["OUT"]

    io = joinpath(dir, fname)
    ext = splitext(fname)[2]
    if ext == ".json"
        newdata = JSON.json(jws, 2)
    elseif ext == ".csv"
        newdata = delimit(jws, ',')
    elseif ext == ".tsv"
        newdata = delimit(jws, '\t')
    else
        throw(ArgumentError("\"$ext\" file type is not supported, use \".json\" or \".csv\""))
    end

    # 편집된 시트만 저장
    modified = true
    if isfile(io)
        modified = !issamedata(read(io, String), newdata)
    end
    if modified
        write(io, newdata)
        print(" SAVE => ")
        printstyled(normpath(io), "\n"; color = :blue)
    else
        print("  ⁿ/ₐ => ")
        print(normpath(io), "\n")
    end
end

function delimit(jws::JSONWorksheet, delim)
    # you cannot use column name from xlsx for the type notation
    # colnames = map(el -> '/' * join(el.tokens, '/'), keys(jws))
    s = join(keys(jws[1]), delim) * '\n'
    for i in 1:length(jws)
        s *= join(map(el -> delimit(el, delim), values(jws[i])), delim)
        if i < length(jws)
            s *= '\n'
        end
    end
    return s
end

delimit(x, delim) = string(x)
delimit(x::Missing, delim) = ""
delimit(x::Nothing, delim) = ""
delimit(x::AbstractString, delim) = x
function delimit(x::AbstractArray, delim)
    "[" * join(x, ';') * "]"
end
function delimit(x::AbstractDict, delim)
    s = "{"
    for (i, (k, v)) in enumerate(x)
        s *= "\"$k\":" * delimit(v, delim)
        if i < length(x)
            s*=";"
        end
    end
    s *= "}"
    return s
end

function write_localize(fname, localizedata)
    config = CACHE["config"]["localization"]
    modified = true

    # TODO: proper warnning, when localization setting isn't there
    baselanguage = get(config, "baseLanguage", "kr")

    # only .json is allowed for localization data
    io, ext = splitext(fname)
    io = joinpath(GAMEENV["LOCALIZE"], "$(io)_$(baselanguage).json")
    newdata = JSON.json(localizedata, 2)

    if isfile(io)
        modified = !issamedata(read(io, String), newdata)
    end
    if modified
        write(io, newdata)
        print("  ⨽Localize => ")
        printstyled(normpath(io), "\n"; color = :blue)
    else
        print("  ⨽Localize => ")
        print(normpath(io), "\n")
    end
end


text:
using DiskCollections
using Test

@testset "DiskCollections.jl" begin
    # Write your tests here.
end


text:
function fqBackPredGSB(x::Array{Float64}, T::Int64, gibbsIter::Int64, burnIn::Int64, gEps1::Float64, gEps2::Float64, thLow::Float64,
  thUp::Float64, zpLow::Float64, zpUp::Float64, zLow::Float64, zUp::Float64, papr::Float64, pbpr::Float64, thin::Int64, samplerSeed::Int64)

	srand(samplerSeed)

  filename = "/Results"

  savelocation = string(pwd(), filename, "/seed$samplerSeed")
  mkpath(savelocation)

	# initialization
	n = length(x)
  nn = Int((gibbsIter - burnIn) / thin)

  if T .> 0 # back-pred
    x = append!(zeros(T), x)
    for i in 1:T
      x[i] = rand(Uniform(-1.,1.))
    end
    ss = n + T
    xbpred = zeros(nn, T)
  else # simple GSBR
    ss = n
  end

  theta = zeros(6) # polynomial coefficients vector
  sampledTheta = zeros(nn, 6) # matrix to store sample thetas

  x0s = zeros(nn, 2)
  sx01, sx02 = 0.5, 0.5

  N = collect(1:1:ss)
  d = collect(1:1:ss)
  noise = zeros(nn)
  clusters = zeros(nn)
  ps = zeros(nn)
  p = 0.5
	toler = 1e-06

	display("Starting MCMC ...")

	##########################################################################
	# MCMC
	##########################################################################

	for iter in 1:gibbsIter

      # Compute weights up to ss
      ##########################################################################

      M = maximum(N)
      w = zeros(M)

      for i in 1:M
        w[i] = p * (1 - p)^(i - 1)
      end

      ## Sample precisions
      ##########################################################################

      Lambdas = zeros(M)

      for j in 1:M
        counter = 0.0
        temp = 0.0
        if d[1] .== j
          counter = counter + 1.0
          temp = (x[1] - henon(theta, sx01, sx02)) ^ 2
        end
        if d[2] .== j
          counter = counter + 1.0
          temp += (x[2] - henon(theta, x[1], sx01)) ^ 2
        end
        for i in 3:ss
          if d[i] .== j
            counter = counter + 1.0
            temp = temp + (x[i] - henon(theta, x[i-1], x[i-2])) ^ 2
          end
        end
        shapel = gEps1 + 0.5 * counter
        ratel = gEps2 + 0.5 * temp
        # println("shape: $shapel")
        # println("rate: $ratel")
        Lambdas[j] = rand(Gamma(shapel, 1. / ratel))
      end

      ## Sample the indicator variables d(i) for i=1,...,ss
      ##########################################################################

      nc1 = 0.0
      prob1 = 0.0
      for j in 1:1:N[1]
        nc1 = nc1 + Lambdas[j]^0.5 * exp(-0.5 * Lambdas[j] * (x[1] - henon(theta, sx01, sx02))^2)
      end
      rd = rand()
      for j in 1:1:N[1]
        prob1 = prob1 + Lambdas[j]^0.5 * exp(-0.5 * Lambdas[j] * (x[1] - henon(theta, sx01, sx02))^2) / nc1
        if rd .< prob1
          d[1] = j
          break
        end
      end

      nc2 = 0.0
      prob2 = 0.0
      for j in 1:1:N[2]
        nc2 = nc2 + Lambdas[j]^0.5 * exp(-0.5 * Lambdas[j] * (x[2] - henon(theta, x[1], sx01))^2)
      end
      rd = rand()
      for j in 1:1:N[2]
        prob2 = prob2 + Lambdas[j]^0.5 * exp(-0.5 * Lambdas[j] * (x[2] - henon(theta, x[1], sx01))^2) / nc2
        if rd .< prob2
          d[2] = j
          break
        end
      end

      for i in 3:1:ss
        nc = 0.0
        prob = 0.0
        for j in 1:1:N[i]
          nc = nc + Lambdas[j]^0.5 * exp(-0.5 * Lambdas[j] * (x[i] - henon(theta, x[i-1], x[i-2]))^2)
        end
        rd = rand()
        for j in 1:1:N[i]
          prob = prob + Lambdas[j]^0.5 * exp(-0.5 * Lambdas[j] * (x[i] - henon(theta, x[i-1], x[i-2]))^2) / nc
          if rd .< prob
            d[i] = j
            break
          end
        end
      end

      # Sample N_i auxilliary variables
      ##########################################################################

      N = etgeornd(p, d)

      ## Number of clusters
      ##########################################################################

      ucl = length(unique(d))
      # println("clusters: $ucl")
      # @printf("uCl: %s \n", ucl)

      ## Sample geometricx probability
      ##########################################################################

      p = rand(Beta(papr + 2 * ss, pbpr + sum(N) - ss))

      ## Sample θ₀ = theta[1]
      ##########################################################################

      muth = Lambdas[d[1]] * (x[1] - (theta[5] * sx01^2 + theta[4] * sx01 * sx02 + theta[6] * sx02^2 + theta[2] * sx01 + theta[3] * sx02)) +
             Lambdas[d[2]] * (x[2] - (theta[5] * x[1]^2 + theta[4] * x[1] * sx01 + theta[6] * sx01^2 + theta[2] * x[1] + theta[3] * sx01))

      tauth = Lambdas[d[1]] + Lambdas[d[2]]

      for j = 3:ss
        muth += Lambdas[d[j]] * (x[j] - (theta[5] * x[j-1]^2 + theta[4] * x[j-1] * x[j-2] + theta[6] * x[j-2]^2 + theta[2] * x[j-1] + theta[3] * x[j-2]))
        tauth += Lambdas[d[j]]
      end
      muth = muth / tauth

      temp = -2. / tauth * log.(rand()) + (theta[1] - muth) ^ 2

      theta[1] = rand(Uniform(max(thLow, muth - temp ^ 0.5), min(thUp, muth + temp ^ 0.5)))

      #
      # ## Sample θ₁ = theta[2]
      # ##########################################################################
      #
      muth =  Lambdas[d[1]] * (x[1] * sx01 - (theta[5] * sx01^3 + theta[4] * sx01^2 * sx02 + theta[6] * sx01 * sx02^2 +
                              theta[3] * sx01 * sx02 + theta[1] * sx01)) +
              Lambdas[d[2]] * (x[2] * x[1] - (theta[5] * x[1]^3 + theta[4] * x[1]^2 * sx01 + theta[6] * x[1] * sx01^2 +
                              theta[3] * x[1] * sx01 + theta[1] * x[1]))

      tauth = Lambdas[d[1]] * sx01 ^ 2 + Lambdas[d[2]] * x[1] ^ 2

      for j = 3:ss
        muth += Lambdas[d[j]] * (x[j] * x[j-1] - (theta[5] * x[j-1]^3 + theta[4] * x[j-1]^2 * x[j-2] + theta[6] * x[j-1] * x[j-2]^2 +
                                theta[3] * x[j-1] * x[j-2] + theta[1] * x[j-1]))
        tauth += Lambdas[d[j]] * x[j - 1] ^ 2
      end
      muth = muth / tauth

      temp = -2. / tauth * log.(rand()) + (theta[2] - muth) ^ 2

      theta[2] = rand(Uniform(max(thLow, muth - temp ^ 0.5), min(thUp, muth + temp ^ 0.5)))

      #
      # ## Sample θ₂ = theta[3]
      # ##########################################################################
      #
      muth =  Lambdas[d[1]] * (x[1] * sx02 - (theta[5] * sx01^2 * sx02 + theta[4] * sx01 * sx02^2 + theta[6] * sx02^3 +
                              theta[2] * sx01 * sx02 + theta[1] * sx02)) +
              Lambdas[d[2]] * (x[2] * sx01 - (theta[5] * x[1]^2 * sx01 + theta[4] * x[1] * sx01^2 + theta[6] * sx01^3 +
                              theta[2] * x[1] * sx01 + theta[1] * sx01))

      tauth = Lambdas[d[1]] * sx02 ^ 2 + Lambdas[d[2]] * sx01 ^ 2

      for j = 3:ss
        muth += Lambdas[d[j]] * (x[j] * x[j-2] - (theta[5] * x[j-1]^2 * x[j-2] + theta[4] * x[j-1] * x[j-2]^2 + theta[6] * x[j-2]^3 +
                                theta[2] * x[j-1] * x[j-2] + theta[1] * x[j-2]))
        tauth += Lambdas[d[j]] * x[j - 2] ^ 2
      end

      muth = muth / tauth

      temp = -2. / tauth * log.(rand()) + (theta[3] - muth) ^ 2

      theta[3] = rand(Uniform(max(thLow, muth - temp ^ 0.5), min(thUp, muth + temp ^ 0.5)))

      #
      # ## Sample θ₃ = theta[4]
      # ##########################################################################
      #
      muth = Lambdas[d[1]] * (x[1] * sx01 * sx02 - (theta[5] * sx01^3 * sx02 + theta[6] * sx01 * sx02^3 + theta[2] * sx01^2 * sx02 +
                              theta[3] * sx01 * sx02^2 + theta[1] * sx01 * sx02)) +
             Lambdas[d[2]] * (x[2] * x[1] * sx01 - (theta[5] * x[1]^3 * sx01 + theta[6] * x[1] * sx01^3 + theta[2] * x[1]^2 * sx01 +
                              theta[3] * x[1] * sx01^2 + theta[1] * x[1] * sx01))

      tauth = Lambdas[d[1]] * sx01 ^ 2 * sx02 ^ 2 + Lambdas[d[2]] * x[1] ^ 2 *sx01 ^ 2

      for j = 3:ss
        muth += Lambdas[d[j]] * (x[j] * x[j-1] * x[j-2] - (theta[5] * x[j-1]^3 * x[j-2] + theta[6] * x[j-1] * x[j-2]^3 + theta[2] * x[j-1]^2 * x[j-2] +
                                theta[3] * x[j-1] * x[j-2]^2 + theta[1] * x[j-1] * x[j-2]))
        tauth += Lambdas[d[j]] * x[j - 1] ^ 2 * x[j - 2] ^ 2
      end

      muth = muth / tauth

      temp = -2. / tauth * log.(rand()) + (theta[4] - muth) ^ 2

      theta[4] = rand(Uniform(max(thLow, muth - temp ^ 0.5), min(thUp, muth + temp ^ 0.5)))

      #
      # ## Sample θ₄ = theta[5]
      # ##########################################################################
      #
      muth = Lambdas[d[1]] * (x[1] * sx01^2 - (theta[4] * sx01^3 * sx02 + theta[6] * sx01^2 * sx02^2 + theta[2] * sx01^3 +
                              theta[3] * sx01^2 * sx02 + theta[1] * sx01^2)) +
             Lambdas[d[2]] * (x[2] * x[1]^2 - (theta[4] * x[1]^3 * sx01 + theta[6] * x[1]^2 * sx01^2 + theta[2] * x[1]^3 +
                              theta[3] * x[1]^2 * sx01 + theta[1] * x[1]^2))

      tauth = Lambdas[d[1]] * sx01 ^ 4 + Lambdas[d[2]] * x[1] ^ 4

      for j = 3:ss
        muth += Lambdas[d[j]] * (x[j] * x[j-1]^2 - (theta[4] * x[j-1]^3 * x[j-2] + theta[6] * x[j-1]^2 * x[j-2]^2 + theta[2] * x[j-1]^3 +
                                theta[3] * x[j-1]^2 * x[j-2] + theta[1] * x[j-1]^2))
        tauth += Lambdas[d[j]] * x[j - 1] ^ 4
      end

      muth = muth / tauth

      temp = -2. / tauth * log.(rand()) + (theta[5] - muth) ^ 2

      theta[5] = rand(Uniform(max(thLow, muth - temp ^ 0.5), min(thUp, muth + temp ^ 0.5)))

      #
      # ## Sample θ₅ = theta[6]
      # ##########################################################################
      #
      muth = Lambdas[d[1]] * (x[1] * sx02^2 - (theta[5] * sx01^2 * sx02^2 + theta[4] * sx01 * sx02^3 + theta[2] * sx01 * sx02^2 +
                              theta[3] * sx02^3 + theta[1] * sx02^2)) +
             Lambdas[d[2]] * (x[2] * sx01^2 - (theta[5] * x[1]^2 * sx01^2 + theta[4] * x[1] * sx01^3 + theta[2] * x[1] * sx01^2 +
                              theta[3] * sx01^3 + theta[1] * sx01^2))

      tauth = Lambdas[d[1]] * sx02 ^ 4 + Lambdas[d[2]] * sx01 ^ 4

      for j = 3:ss
        muth += Lambdas[d[j]] * (x[j] * x[j-2]^2 - (theta[5] * x[j-1]^2 * x[j-2]^2 + theta[4] * x[j-1] * x[j-2]^3 + theta[2] * x[j-1] * x[j-2]^2 +
                                theta[3] * x[j-2]^3 + theta[1] * x[j-2]^2))
        tauth += Lambdas[d[j]] * x[j - 2] ^ 4
      end

      muth = muth / tauth

      temp = -2. / tauth * log.(rand()) + (theta[6] - muth) ^ 2

      theta[6] = rand(Uniform(max(thLow, muth - temp ^ 0.5), min(thUp, muth + temp ^ 0.5)))

      ## Sample past unobserved observations
      ##########################################################################

      if T .> 0

        if T .> 1
          for i in 1:1:(T-2)

            a4 = Lambdas[d[T-i+2]] * theta[5] ^ 2 + Lambdas[d[T-i+3]] * theta[6] ^ 2

            a3 = 2. * Lambdas[d[T-i+2]] * (theta[4] * theta[5] * x[T-i] + theta[2] * theta[5]) + 2. * Lambdas[d[T-i+3]] * (theta[4] * theta[6] * x[T-i+2] + theta[3] * theta[6])

            a2 = Lambdas[d[T-i+2]] * (theta[4]^2 * x[T-i]^2 + 2. * theta[5] * theta[6] * x[T-i]^2 + 2. * theta[2] * theta[4] * x[T-i] + 2. * theta[3] * theta[5] * x[T-i] -
                                2. * theta[5] * x[T-i+2] + 2. * theta[1] * theta[5] + theta[2]^2) +
               Lambdas[d[T-i+3]] * (theta[4]^2 * x[T-i+2]^2 + 2. * theta[5] * theta[6] * x[T-i+2]^2 + 2 * theta[2] * theta[6] * x[T-i+2] + 2 * theta[3] * theta[4] * x[T-i+2] -
                                2. * theta[6] * x[T-i+3] + 2. * theta[1] * theta[6] + theta[3]^2) +
                Lambdas[d[T-i+1]]

            a1 = 2. * (Lambdas[d[T-i+2]] * (theta[4] * theta[6] * x[T-i]^3 + (theta[2] * theta[6] + theta[3] * theta[4]) * x[T-i]^2 + (-theta[4] * x[T-i+2] + theta[1] * theta[4] + theta[2] * theta[3]) * x[T-i] -
                                      theta[2] * x[T-i+2] + theta[1] * theta[2]) +
                    Lambdas[d[T-i+3]] * (theta[4] * theta[5] * x[T-i+2]^3 + (theta[2] * theta[4] + theta[3] * theta[5]) * x[T-i+2]^2 -
                                    theta[4] * x[T-i+2] * x[T-i+3] + (theta[1] * theta[4] + theta[2] * theta[3]) * x[T-i+2] - theta[3] * x[T-i+3] + theta[1] * theta[3]) -
                    Lambdas[d[T-i+1]] * (theta[4] * x[T-i-1] * x[T-i] + theta[5] * x[T-i]^2 + theta[6] * x[T-i-1]^2 + theta[2] * x[T-i] + theta[3] * x[T-i-1] + theta[1])  )

            aux = -2. * log.(rand()) + a1 * x[T-i+1] + a2 * x[T-i+1] ^ 2 + a3 * x[T-i+1] ^ 3 + a4 * x[T-i+1] ^ 4

            poly = [-aux / a4; a1 / a4; a2 / a4; a3 / a4; 1.0]
            allRoots = Polynomials.roots(Poly(poly))
            sel = abs.(imag(allRoots)) .< toler # treat as real the roots with imagine part < ϵ
            realRoots = sort(real(allRoots[sel]))
            intervals = rangeIntersection(realRoots, [zpLow; zpUp])

            x[T-i+1] = unifmixrnd(intervals) ##### go backwards ##### !!!!!!!!!!!!!!!!!!!!!!

          end

          # 2nd last unobserved past observation
          a4 = Lambdas[d[3]] * theta[5] ^ 2 + Lambdas[d[4]] * theta[6] ^ 2

          a3 = 2. * Lambdas[d[3]] * (theta[4] * theta[5] * x[1] + theta[2] * theta[5]) + 2. * Lambdas[d[4]] * (theta[4] * theta[6] * x[3] + theta[3] * theta[6])

          a2 = Lambdas[d[3]] * (theta[4]^2 * x[1]^2 + 2. * theta[5] * theta[6] * x[1]^2 + 2. * theta[2] * theta[4] * x[1] + 2. * theta[3] * theta[5] * x[1] -
                                2. * theta[5] * x[3] + 2. * theta[1] * theta[5] + theta[2]^2) +
               Lambdas[d[4]] * (theta[4]^2 * x[3]^2 + 2. * theta[5] * theta[6] * x[3]^2 + 2 * theta[2] * theta[6] * x[3] + 2 * theta[3] * theta[4] * x[3] -
                                2. * theta[6] * x[4] + 2. * theta[1] * theta[6] + theta[3]^2) +
               Lambdas[d[2]]

          a1 = 2. * (Lambdas[d[3]] * (theta[4] * theta[6] * x[1]^3 + (theta[2] * theta[6] + theta[3] * theta[4]) * x[1]^2 + (-theta[4] * x[3] + theta[1] * theta[4] + theta[2] * theta[3]) * x[1] -
                                      theta[2] * x[3] + theta[1] * theta[2]) +
                    Lambdas[d[4]] * (theta[4] * theta[5] * x[3]^3 + (theta[2] * theta[4] + theta[3] * theta[5]) * x[3]^2 -
                                    theta[4] * x[3] * x[4] + (theta[1] * theta[4] + theta[2] * theta[3]) * x[3] - theta[3] * x[4] + theta[1] * theta[3]) -
                    Lambdas[d[2]] * (theta[4] * x[1] * sx01 + theta[5] * x[1]^2 + theta[6] * sx01^2 + theta[2] * x[1] + theta[3] * sx01 + theta[1])  )

          aux = -2. * log.(rand()) + a1 * x[2] + a2 * x[2] ^ 2 + a3 * x[2] ^ 3 + a4 * x[2] ^ 4

          poly = [-aux / a4; a1 / a4; a2 / a4; a3 / a4; 1.0]
          allRoots = Polynomials.roots(Poly(poly))
          sel = abs.(imag(allRoots)) .< toler # treat as real the roots with imagine part < ϵ
          realRoots = sort(real(allRoots[sel]))
          intervals = rangeIntersection(realRoots, [zpLow; zpUp])
          x[2] = unifmixrnd(intervals)

        end

        # Last unobserved past observation
        a4 = Lambdas[d[2]] * theta[5] ^ 2 + Lambdas[d[3]] * theta[6] ^ 2

        a3 = 2. * Lambdas[d[2]] * (theta[4] * theta[5] * sx01 + theta[2] * theta[5]) + 2. * Lambdas[d[3]] * (theta[4] * theta[6] * x[2] + theta[3] * theta[6])

        a2 = Lambdas[d[2]] * (theta[4]^2 * sx01^2 + 2. * theta[5] * theta[6] * sx01^2 + 2. * theta[2] * theta[4] * sx01 + 2. * theta[3] * theta[5] * sx01 -
                              2. * theta[5] * x[2] + 2. * theta[1] * theta[5] + theta[2]^2) +
             Lambdas[d[3]] * (theta[4]^2 * x[2]^2 + 2. * theta[5] * theta[6] * x[2]^2 + 2 * theta[2] * theta[6] * x[2] + 2 * theta[3] * theta[4] * x[2] -
                              2. * theta[6] * x[3] + 2. * theta[1] * theta[6] + theta[3]^2) +
             Lambdas[d[1]]

        a1 = 2. * (Lambdas[d[2]] * (theta[4] * theta[6] * sx01^3 + (theta[2] * theta[6] + theta[3] * theta[4]) * sx01^2 + (-theta[4] * x[2] + theta[1] * theta[4] + theta[2] * theta[3]) * sx01 -
                                    theta[2] * x[2] + theta[1] * theta[2]) +
                  Lambdas[d[3]] * (theta[4] * theta[5] * x[2]^3 + (theta[2] * theta[4] + theta[3] * theta[5]) * x[2]^2 -
                                  theta[4] * x[2] * x[3] + (theta[1] * theta[4] + theta[2] * theta[3]) * x[2] - theta[3] * x[3] + theta[1] * theta[3]) -
                  Lambdas[d[1]] * (theta[4] * sx01 * sx02 + theta[5] * sx01^2 + theta[6] * sx02^2 + theta[2] * sx01 + theta[3] * sx02 + theta[1])  )

        aux = -2. * log.(rand()) + a1 * x[1] + a2 * x[1] ^ 2 + a3 * x[1] ^ 3 + a4 * x[1] ^ 4

        poly = [-aux / a4; a1 / a4; a2 / a4; a3 / a4; 1.0]
        allRoots = Polynomials.roots(Poly(poly))
        sel = abs.(imag(allRoots)) .< toler # treat as real the roots with imagine part < ϵ
        realRoots = sort(real(allRoots[sel]))
        intervals = rangeIntersection(realRoots, [zpLow; zpUp])
        x[1] = unifmixrnd(intervals)

      end

      ## Sample initial conditions x₀,y₀ = x0,x00
      ##########################################################################

      # x₀

      a4 = Lambdas[d[1]] * theta[5] ^ 2 + Lambdas[d[2]] * theta[6] ^ 2

      a3 = 2. * Lambdas[d[1]] * (theta[4] * theta[5] * sx02 + theta[2] * theta[5]) + 2. * Lambdas[d[2]] * (theta[4] * theta[6] * x[1] + theta[3] * theta[6])

      a2 = Lambdas[d[1]] * (theta[4]^2 * sx02^2 + 2. * theta[5] * theta[6] * sx02^2 + 2. * theta[2] * theta[4] * sx02 + 2. * theta[3] * theta[5] * sx02 -
                            2. * theta[5] * x[1] + 2. * theta[1] * theta[5] + theta[2]^2) +
           Lambdas[d[2]] * (theta[4]^2 * x[1]^2 + 2. * theta[5] * theta[6] * x[1]^2 + 2 * theta[2] * theta[6] * x[1] + 2 * theta[3] * theta[4] * x[1] -
                            2. * theta[6] * x[2] + 2. * theta[1] * theta[6] + theta[3]^2)

      a1 = 2. * (Lambdas[d[1]] * (theta[4] * theta[6] * sx02^3 + (theta[2] * theta[6] + theta[3] * theta[4]) * sx02^2 + (-theta[4] * x[1] + theta[1] * theta[4] + theta[2] * theta[3]) * sx02 -
                                  theta[2] * x[1] + theta[1] * theta[2]) +
                Lambdas[d[2]] * (theta[4] * theta[5] * x[1]^3 + (theta[2] * theta[4] + theta[3] * theta[5]) * x[1]^2 -
                                theta[4] * x[1] * x[2] + (theta[1] * theta[4] + theta[2] * theta[3]) * x[1] - theta[3] * x[2] + theta[1] * theta[3]))

      aux = -2. * log.(rand()) + a1 * sx01 + a2 * sx01 ^ 2 + a3 * sx01 ^ 3 + a4 * sx01 ^ 4

      poly = [-aux / a4; a1 / a4; a2 / a4; a3 / a4; 1.0]
      #allRoots = roots(poly, polish=true, epsilon=1e-20)
      allRoots = Polynomials.roots(Poly(poly))

      #println("allRoots = $allRoots")
      # sel = imag(allRoots) .== 0.0
      # realRoots = sort(real(allRoots[sel]))

      sel = abs.(imag(allRoots)) .< toler # treat as real the roots with imagine part < ϵ
      realRoots = sort(real(allRoots[sel]))
      intervals = rangeIntersection(realRoots, [zLow; zUp])
      sx01 = unifmixrnd(intervals)
      #
      # x00 = y₀

      a4  = Lambdas[d[1]] * theta[6] ^ 2

      a3 = 2. * theta[6] * Lambdas[d[1]] * (theta[4] * sx01 + theta[3])

      a2 = Lambdas[d[1]] * ( -2. * theta[6] * (x[1] - theta[1] - theta[2] * sx01 - theta[5] * sx01^2) + (theta[4] * sx01 + theta[3])^2 )

      a1 = -2. * Lambdas[d[1]] * (x[1] - theta[1] - theta[5] * sx01^2 - theta[2] * sx01) * (theta[4] * sx01 + theta[3])

      aux = -2. * log.(rand()) + a1 * sx02 + a2 * sx02 ^ 2 + a3 * sx02 ^ 3 + a4 * sx02 ^ 4

      poly = [-aux / a4; a1 / a4; a2 / a4; a3 / a4; 1.0]
      allRoots = Polynomials.roots(Poly(poly))
      sel = abs.(imag(allRoots)) .< toler # treat as real the roots with imagine part < ϵ
      # sel = imag(allRoots) .== 0.0 # treat as real the roots with imagine part < ϵ
      realRoots = sort(real(allRoots[sel]))
      intervals = rangeIntersection(realRoots, [zLow; zUp])
      sx02 = unifmixrnd(intervals)

      ## After Burn-In period
      ###############################wwwwwwwww###########################################

      if (iter .> burnIn) & ((iter-burnIn) % thin .== 0)

        ii = Int((iter - burnIn)/thin)

        ## Sample noise predictive
        ##########################################################################

        cW = cumsum(w)
        flag = rand()
        if cW[end] .< flag
          pred = rand(Normal(0, sqrt(1 / rand(Gamma(gEps1, 1 / gEps2))))) # draw from the prior
        else
          for j in 1:1:length(cW)
            if flag .< cW[j]
               pred = rand(Normal(0, sqrt(1 / Lambdas[j])))
               break
            end
          end
        end

        # Store values
        sampledTheta[ii, :] = theta
        x0s[ii, :] = [sx01 sx02]
        noise[ii] = pred
        clusters[ii] = ucl
        ps[ii] = p

        if T .> 0
          xbpred[ii,:] = x[1:1:T]
        end

      end

      if iter % 50000 .== 0
         println("MCMC Iterations: $iter")
      end

	end

	display("... MCMC finished !")

	## Write values in .txt files - specific path
	##########################################################################

  if T .> 0
    writedlm(string(savelocation,"/xbpred.txt"), xbpred)
  end
	writedlm(string(savelocation, "/thetas.txt"), sampledTheta)
	writedlm(string(savelocation,"/x0s.txt"), x0s)
  writedlm(string(savelocation,"/noise.txt"), noise, '\n')
	writedlm(string(savelocation,"/ucl.txt"), clusters, '\n')
  writedlm(string(savelocation,"/ps.txt"), ps, '\n')

 # #, xbpred
 return x0s

end

max_stars_repo_name: KosKal/SM-GSBR

max_issues_repo_name: KosKal/SM-GSBR

max_forks_repo_name: KosKal/SM-GSBR


text:
using DiffPrecTest
using Statistics, StatsBase, LinearAlgebra
using SparseArrays
using ProximalBase, CoordinateDescent, CovSel
using Random, Distributions
using JLD

@show gethostname()

rep   = parse(Int,ARGS[1])
ip    = parse(Int,ARGS[2])
ialpha = parse(Int,ARGS[3])
dir = ARGS[4]

alpha = 0.:0.1:1.
pArr = [50, 100, 150]
n = 300
NUM_REP = 1000

p = pArr[ip]
Random.seed!(7689)

# generate X
Ωx = Matrix{Float64}(I, p, p)
for l=1:p-1
  Ωx[l  , l+1] = 0.3
  Ωx[l+1, l  ] = 0.3
end
for l=1:p-2
  Ωx[l  , l+2] = 0.2
  Ωx[l+2, l  ] = 0.2
end
Ωy = Matrix{Float64}(I, p, p) * 1.1
for l=1:p-1
  Ωy[l  , l+1] = 0.3
  Ωy[l+1, l  ] = 0.3
end
for l=1:p-2
  Ωy[l  , l+2] = -0.1
  Ωy[l+2, l  ] = -0.1
end
d = Vector{Float64}(undef, p)
rand!(Uniform(0.5, 2.5), d)
d .= sqrt.(d)
D = Diagonal(d)
Ωx = D * Ωx * D
Ωy = D * Ωy * D

Ωy = (1. - alpha[ialpha]) * Ωx + alpha[ialpha] * Ωy
Δ = Ωx - Ωy

Σx = inv(Symmetric(Ωx))
Σy = inv(Symmetric(Ωy))
dist_X = MvNormal(convert(Matrix, Σx))
dist_Y = MvNormal(convert(Matrix, Σy))

rp = div(p * (p+1), 2)
res = Matrix{DiffPrecResultNormal}(undef, rp, NUM_REP)

# generate data
Random.seed!(12346 + rep*100 + ialpha)
X = rand(dist_X, n)'
Y = rand(dist_Y, n)'

Sx = Symmetric( X'X / n )
Sy = Symmetric( Y'Y / n )

S = BitArray( abs(Δ[r, c]) > 1e-4 for r = 1:p, c = 1:p )
eS_init = Array{BitArray}(undef, div((p + 1)*p, 2))
for j=1:div((p + 1)*p, 2)
    eS_init[j] = S
end

@time boot_res, eS = bootstrap(X, Y; estimSupport=eS_init)

@save "$(dir)/res_$(rep).jld" boot_res


text:
################################################################################
# Box (2D)
################################################################################
struct Box2D end

function box2D(;
    representation=:minimal,
    timestep=0.05,
    gravity=[0.0; 0.0; -9.81],
    friction_coefficient=0.8,
    info=nothing,
    seed=1,
    contact=true,
    vis=Visualizer(),
    name=:robot,
    infeasible_control=false,
    opts_step=SolverOptions(rtol=3.0e-4, btol=3.0e-4, undercut=1.5),
    opts_grad=SolverOptions(rtol=3.0e-4, btol=3.0e-4, undercut=1.5),
    T=Float64)

    mechanism = get_mechanism(:box2D,
        timestep=timestep,
        gravity=gravity,
        friction_coefficient=friction_coefficient)

    initialize!(mechanism, :box2D)

    if representation == :minimal
        nx = minimal_dimension(mechanism)
    elseif representation == :maximal
        nx = maximal_dimension(mechanism)
    end
    nu_inf = input_dimension(mechanism)
    nu = infeasible_control ? nu_inf : nu_inf - 2 # remove last 2 controls
    no = nx

    aspace = BoxSpace(nu,
        low=(-1.0e8 * ones(nu)),
        high=(1.0e8 * ones(nu)))
    ospace = BoxSpace(no,
        low=(-Inf * ones(no)),
        high=(Inf * ones(no)))

    rng = MersenneTwister(seed)
    z = get_maximal_state(mechanism)
    x = representation == :minimal ? maximal_to_minimal(mechanism, z) : z
    fx = zeros(nx, nx)
    fu = zeros(nx, nu)

    u_prev = zeros(nu)
    control_mask = infeasible_control ? I(nu) : [I(nu) zeros(nu, 2)]
    control_scaling = Diagonal(ones(nu))

    build_robot(mechanism, vis=vis, name=name)

    TYPES = [Quadruped, T, typeof(mechanism), typeof(aspace), typeof(ospace), typeof(info)]
    Environment{TYPES...}(mechanism, representation, aspace, ospace,
        x, fx, fu,
        u_prev, control_mask' * control_scaling,
        nx, nu, no,
        info,
        [rng], vis,
        opts_step, opts_grad)
end


text:
# ======================
# Multivariate functions
# ======================
DAISY["himmilbeau"] = Dict("dom" => a×b,
    "numvars" => 2,
    "vars" => "x1 x2",
    "f" => ((x1, x2) -> (x1*x1 + x2 - 11)*(x1 * x1 + x2 - 11) + (x1 + x2*x2 - 7)*(x1 + x2*x2 - 7)),
    "ref" => Interval(85.46830677734748, 221.7338939301446) # MOSEK deg 6
    )

DAISY["kepler1"] = Dict("dom" => a×b×c×d×e×f,
    "numvars" => 6,
    "vars" => "x1 x2 x3 x4 x5 x6",
    "f" => ((x1, x2, x3, x4, x5, x6) -> x2 * x5 + x3 * x6 - x2 * x3 - x5 * x6 + x1 * (-x1 + x2 + x3 - x4 + x5 + x6)),
    "ref" => Interval(-5.255935494810441, 7.321362422825775)  # MOSEK deg 6
    )

DAISY["kepler2"] = Dict("dom" => a×b×c×d,
    "numvars" => 4,
    "vars" => "x1 x2 x3 x4",
    "f" => ((x1, x2, x3, x4) -> x1 * x4 * (-x1 + x2 + x3 - x4) +
                                x2 * (x1 - x2 + x3 + x4) + x3 * (x1 + x2 - x3 + x4) -x2 * x3 * x4 - x1 * x3 - x1 * x2 - x4),
    "ref" => Interval(-195.36974909125482, 78.3669520644375)  # MOSEK deg 6
    )

DAISY["kepler3"] = Dict("dom" => a×b×c×d×e×f,
    "numvars" => 6,
    "vars" => "x1 x2 x3 x4 x5 x6",
    "f" => ((x1, x2, x3, x4, x5, x6) -> x1 * x4 * (-x1 + x2 + x3 - x4 + x5 + x6) +
                                        x2 * x5 * (x1 - x2 + x3 + x4 - x5 + x6) +x3* x6 * (x1 + x2 - x3 + x4 + x5 - x6) -
                                        x2 * x3 * x4 -x1* x3* x5 - x1 * x2 * x6 - x4 * x5 * x6),
    "ref" => Interval(-309.8484155131222, 17.982082401462407) # MOSEK deg 6
    )

DAISY["Rigidbody1"] = Dict("dom" => a×b×c,
    "numvars" => 3,
    "vars" => "x1 x2 x3",
    "f" => ((x1, x2, x3) -> -x1*x2 - 2*x2*x3 - x1 - x3),
    "ref" => Interval(-20.786552979420335, -0.540012836551535) # MOSEK deg 6
    )

DAISY["Rigidbody2"] = Dict("dom" => a×b×c,
    "numvars" => 3,
    "vars" => "x1 x2 x3",
    "f" => ((x1, x2, x3) -> 2*(x1*x2*x3) + (3*x3*x3) - x2*(x1*x2*x3) + (3*x3*x3) - x2),
    "ref" => Interval(68.81138021006673, 359.98566570476504)  # MOSEK deg 6
    )


text:
# ----------- definitions of possible symmetries a 1-p basis may possess

# TODO: attach these properties to the one-particle basis sets

abstract type AbstractSymmetry end

struct EuclideanO3Equivariance  <: AbstractSymmetry end

struct SphericalO3Equivariance  <: AbstractSymmetry end

struct O3Invariance <: AbstractSymmetry end

groupaction(X::AbstractState) = groupaction(typeof(X))

# ----------- This file implements the abstract one-particle basis interface

alloc_B(basis::OneParticleBasis) = zeros( fltype(basis), length(basis) )

alloc_dB(basis::OneParticleBasis, N::Integer = 1) =
      zeros(gradtype(basis), (length(basis), N))

function evaluate!(A, tmp, basis::OneParticleBasis,
                   cfg::AbstractConfiguration)
   fill!(A, 0)
   for X in cfg
      add_into_A!(A, tmp, basis, X)
   end
   return A
end

function evaluate!(A, tmp, basis::OneParticleBasis, X::AbstractState)
   fill!(A, 0)
   add_into_A!(A, tmp, basis, X)
   return A
end

function evaluate_ed!(A, dA, tmpd, basis::OneParticleBasis,
                      cfg::AbstractConfiguration)
   fill!(A, 0)
   fill!(dA, zero(eltype(dA)))  # TODO: this should not be necessary!
   for (j, X) in enumerate(cfg)
      dAview = @view dA[:, j]
      add_into_A_dA!(A, dAview, tmpd, basis, X)
   end
   return dA
end

function evaluate_ed!(A, dA, tmpd, basis::OneParticleBasis, X::AbstractState)
   fill!(A, 0)
   add_into_A_dA!(A, dA, tmpd, basis, X)
   return dA
end

# # TODO: fix the dispatch structure on alloc_dB
# alloc_dB(basis::OneParticleBasis, ::JVec) = alloc_dB(basis, 1)
# alloc_dB(basis::OneParticleBasis, Rs::AbstractVector{<: JVec}, args...) =
#       alloc_dB(basis, length(Rs))
# alloc_dB(basis::OneParticleBasis) = alloc_dB(basis, 1)
#
# function alloc_dB(basis::OneParticleBasis, maxN::Integer)
#    NZ = numz(basis)
#    maxlen = maximum( sum( length(basis, iz, iz0) for iz = 1:NZ )
#                      for iz0 = 1:NZ )
#    T = fltype(basis)
#    return zeros(JVec{T}, (maxlen, maxN))
# end

# function set_Aindices!(basis::OneParticleBasis)
#    NZ = numz(basis)
#    for iz0 = 1:NZ
#       idx = 0
#       for iz = 1:NZ
#          len = length(basis, iz, iz0)
#          basis.Aindices[iz, iz0] = (idx+1):(idx+len)
#          idx += len
#       end
#    end
#    return basis
# end
#
# get_basis_spec(basis::OneParticleBasis, iz0::Integer) =
#       get_basis_spec(basis, i2z(iz0))
#
# get_basis_spec(basis::OneParticleBasis, s::Symbol) =
#       get_basis_spec(basis, atomic_number(s))

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

"""
`struct One1p` : the 1-p basis where everything just maps to 1
"""
struct One1pBasis <: OneParticleBasis{Bool}
end

Base.length(::One1pBasis) = 1

function evaluate!(B, tmp, basis::One1pBasis, args...)
   B[1] = 1
   return B
end

function evaluate_d!(dB, tmp, basis::One1pBasis, args...)
   dB[1] = zero(eltype(dB))
   return dB
end

fltype(::One1pBasis) = Bool

symbols(::One1pBasis) = Symbol[]

indexrange(::One1pBasis) = Dict()

isadmissible(b, ::One1pBasis) = true

degree(b, ::One1pBasis) = 0

get_spec(::One1pBasis) = [(1,)]


text:
immutable Poisson <: DiscreteUnivariateDistribution
    lambda::Float64
    function Poisson(l::Real)
    	l > zero(l) || error("lambda must be positive")
        new(float64(l))
    end
    Poisson() = new(1.0)
end

@_jl_dist_1p Poisson pois

function entropy(d::Poisson)
    λ = d.lambda
    if λ < 50.0
        s = 0.0
        for k in 1:100
            s += λ^k * lgamma(k + 1.0) / gamma(k + 1.0)
        end
        return λ * (1.0 - log(λ)) + exp(-λ) * s
    else
        return 0.5 * log(2 * pi * e * λ) -
               (1 / (12 * λ)) -
               (1 / (24 * λ * λ)) -
               (19 / (360 * λ * λ * λ))
    end
end

insupport(::Poisson, x::Real) = isinteger(x) && zero(x) <= x
insupport(::Type{Poisson}, x::Real) = isinteger(x) && zero(x) <= x

kurtosis(d::Poisson) = 1.0 / d.lambda

mean(d::Poisson) = d.lambda

median(d::Poisson) = quantile(d, 0.5)

function mgf(d::Poisson, t::Real)
    l = d.lambda
    return exp(l * (exp(t) - 1.0))
end

function cf(d::Poisson, t::Real)
    l = d.lambda
    return exp(l * (exp(im * t) - 1.0))
end

mode(d::Poisson) = ifloor(d.lambda)
modes(d::Poisson) = [mode(d)]

skewness(d::Poisson) = 1.0 / sqrt(d.lambda)

var(d::Poisson) = d.lambda

function fit_mle(::Type{Poisson}, x::Array)
    for i in 1:length(x)
        if !insupport(Poisson(), x[i])
            error("Poisson observations must be non-negative integers")
        end
    end
    Poisson(mean(x))
end


text:
# Autogenerated wrapper script for ODEInterface_jll for powerpc64le-linux-gnu-libgfortran5
export libbvp_m_proxy, libbvpsol, libbvpsol_i32, libcolnew, libcolnew_i32, libddeabm, libddeabm_i32, libddebdf, libddebdf_i32, libdop853, libdop853_i32, libdopri5, libdopri5_i32, libodex, libodex_i32, libradau, libradau5, libradau5_i32, libradau_i32, librodas, librodas_i32, libseulex, libseulex_i32

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("ODEInterface")
JLLWrappers.@declare_library_product(libbvp_m_proxy, "bvp_m_proxy.so")
JLLWrappers.@declare_library_product(libbvpsol, "bvpsol.so")
JLLWrappers.@declare_library_product(libbvpsol_i32, "bvpsol_i32.so")
JLLWrappers.@declare_library_product(libcolnew, "colnew.so")
JLLWrappers.@declare_library_product(libcolnew_i32, "colnew_i32.so")
JLLWrappers.@declare_library_product(libddeabm, "ddeabm.so")
JLLWrappers.@declare_library_product(libddeabm_i32, "ddeabm_i32.so")
JLLWrappers.@declare_library_product(libddebdf, "ddebdf.so")
JLLWrappers.@declare_library_product(libddebdf_i32, "ddebdf_i32.so")
JLLWrappers.@declare_library_product(libdop853, "dop853.so")
JLLWrappers.@declare_library_product(libdop853_i32, "dop853_i32.so")
JLLWrappers.@declare_library_product(libdopri5, "dopri5.so")
JLLWrappers.@declare_library_product(libdopri5_i32, "dopri5_i32.so")
JLLWrappers.@declare_library_product(libodex, "odex.so")
JLLWrappers.@declare_library_product(libodex_i32, "odex_i32.so")
JLLWrappers.@declare_library_product(libradau, "radau.so")
JLLWrappers.@declare_library_product(libradau5, "radau5.so")
JLLWrappers.@declare_library_product(libradau5_i32, "radau5_i32.so")
JLLWrappers.@declare_library_product(libradau_i32, "radau_i32.so")
JLLWrappers.@declare_library_product(librodas, "rodas.so")
JLLWrappers.@declare_library_product(librodas_i32, "rodas_i32.so")
JLLWrappers.@declare_library_product(libseulex, "seulex.so")
JLLWrappers.@declare_library_product(libseulex_i32, "seulex_i32.so")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libbvp_m_proxy,
        "lib/bvp_m_proxy.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libbvpsol,
        "lib/bvpsol.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libbvpsol_i32,
        "lib/bvpsol_i32.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcolnew,
        "lib/colnew.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcolnew_i32,
        "lib/colnew_i32.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libddeabm,
        "lib/ddeabm.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libddeabm_i32,
        "lib/ddeabm_i32.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libddebdf,
        "lib/ddebdf.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libddebdf_i32,
        "lib/ddebdf_i32.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libdop853,
        "lib/dop853.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libdop853_i32,
        "lib/dop853_i32.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libdopri5,
        "lib/dopri5.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libdopri5_i32,
        "lib/dopri5_i32.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libodex,
        "lib/odex.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libodex_i32,
        "lib/odex_i32.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libradau,
        "lib/radau.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libradau5,
        "lib/radau5.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libradau5_i32,
        "lib/radau5_i32.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libradau_i32,
        "lib/radau_i32.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        librodas,
        "lib/rodas.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        librodas_i32,
        "lib/rodas_i32.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libseulex,
        "lib/seulex.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libseulex_i32,
        "lib/seulex_i32.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module MaxHelpingHandFlagToggleCameraOffsetTrigger

using ..Ahorn, Maple

@mapdef Trigger "MaxHelpingHand/FlagToggleCameraOffsetTrigger" FlagToggleCameraOffsetTrigger(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight,
    cameraX::Number=0.0, cameraY::Number=0.0, flag::String="flag_toggle_camera_offset", inverted::Bool=false)

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Camera Offset (Flag-Toggled) (max480's Helping Hand)" => Ahorn.EntityPlacement(
        FlagToggleCameraOffsetTrigger,
        "rectangle"
    )
)

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
type ProdKernel <: Kernel
    kerns::Vector{Kernel}
    function ProdKernel(args...)
        kerns = Array(Kernel, length(args))
        for (i,k) in enumerate(args)
            isa(k, Kernel) || throw(ArgumentError("All arguments of ProdKernel must be Kernel objects"))
            kerns[i] = k
        end
        return new(kerns)
    end
end

function show(io::IO, prodkern::ProdKernel, depth::Int = 0)
    pad = repeat(" ", 2 * depth)
    println(io, "$(pad)Type: $(typeof(prodkern))")
    for k in prodkern.kerns
        show(io, k, depth+1)
    end
end

function cov(prodkern::ProdKernel, x::Vector{Float64}, y::Vector{Float64})
    p = 1.0
    for k in prodkern.kerns
        p = p.*cov(k, x, y)
    end
    return p
end

function cov(prodkern::ProdKernel, X::Matrix{Float64})
    d, nobsv = size(X)
    p = ones(nobsv, nobsv)
    for k in prodkern.kerns
        p[:,:] = p .* cov(k, X)
    end
    return p
end

function get_params(prodkern::ProdKernel)
    p = Array(Float64, 0)
    for k in prodkern.kerns
        append!(p, get_params(k))
    end
    p
end

get_param_names(prodkern::ProdKernel) = composite_param_names(prodkern.kerns, :pk)

function num_params(prodkern::ProdKernel)
    n = 0
    for k in prodkern.kerns
        n += num_params(k)
    end
    n
end

function set_params!(prodkern::ProdKernel, hyp::Vector{Float64})
    i, n = 1, num_params(prodkern)
    length(hyp) == num_params(prodkern) || throw(ArgumentError("ProdKernel object requires $(n) hyperparameters"))
    for k in prodkern.kerns
        np = num_params(k)
        set_params!(k, hyp[i:(i+np-1)])
        i += np
    end
end

# This function is extremely inefficient
function grad_kern(prodkern::ProdKernel, x::Vector{Float64}, y::Vector{Float64})
     dk = Array(Float64, 0)
      for k in prodkern.kerns
          p = 1.0
          for j in prodkern.kerns[find(k.!=prodkern.kerns)]
              p = p.*cov(j, x, y)
          end
        append!(dk,grad_kern(k, x, y).*p)
      end
    dk
end

function grad_stack!(stack::AbstractArray, X::Matrix{Float64}, prodkern::ProdKernel)
    d, nobsv = size(X)
    num_kerns = length(prodkern.kerns)

    cross_kerns = Array(Float64, nobsv, nobsv, num_kerns)
    for (i,kern) in enumerate(prodkern.kerns)
        cross_kerns[:,:,i] = cov(kern, X)
    end

    s = 1
    for (i,kern) in enumerate(prodkern.kerns)
        np = num_params(kern)
        grad_stack!(view(stack,:,:,s:(s+np-1)), X, kern)
        for j in 1:num_kerns
            if j != i
                broadcast!(.*, view(stack,:,:,s:(s+np-1)), view(stack,:,:,s:(s+np-1)), view(cross_kerns, :,:,j))
            end
        end
        s += np
    end
    stack
end

# Multiplication operators
function *(k1::ProdKernel, k2::Kernel)
    kerns = [k1.kerns; k2]
    ProdKernel(kerns...)
end
function *(k1::ProdKernel, k2::ProdKernel)
    kerns = [k1.kerns; k2.kerns]
    ProdKernel(kerns...)
end
*(k1::Kernel, k2::Kernel) = ProdKernel(k1,k2)
*(k1::Kernel, k2::ProdKernel) = *(k2,k1)


text:
# https://www.codewars.com/kata/5583d268479559400d000064

function binarytostring(binary)
  if length(binary) < 8
    ""
  else
    ascii = parse(Int, binary[1:8], base=2)
    Char(ascii) * binarytostring(binary[9:end])
  end
end

max_stars_repo_name: hiljusti/codewars-solutions

max_issues_repo_name: hiljusti/codewars-solutions

max_forks_repo_name: hiljusti/codewars-solutions


text:
"""
	project(cmd0::String="", arg1=nothing, kwargs...)

Project data onto lines or great circles, generate tracks, or translate coordinates.

Full option list at [`project`](http://docs.generic-mapping-tools.org/latest/project.html)

Parameters
----------

- **C** | **origin** | **start_pt** :: [Type => list/tuple]    ``Arg = (x,y)``

    Sets the origin of the projection, in Definition 1 or 2.
    (http://docs.generic-mapping-tools.org/latest/project.html#c)

- **A** | **azim** :: [Type => Number]    ``Arg = azimuth``

    Defines the azimuth of the projection (Definition 1).
    (http://docs.generic-mapping-tools.org/latest/project.html#a)
- **E** | **end_pt** :: [Type => list/tuple]    ``Arg = (bx,by)``

    bx,by defines the end point of the projection path (Definition 2).
    (http://docs.generic-mapping-tools.org/latest/project.html#e)
- **F** | **out_flags** :: [Type => Str]    ``Arg = xyzpqrs``

    Specify your desired output using any combination of xyzpqrs, in any order [Default is xyzpqrs].
    (http://docs.generic-mapping-tools.org/latest/project.html#f)
- **G** | **step** | **small_circle**: [Type => Number or list/tuple--    ``Arg = dist[/colat][+h]``

    Generate mode. No input is read. Create (r, s, p) output points every dist units of p. See Q option.
    (http://docs.generic-mapping-tools.org/latest/project.html#g)
- **L** | **length_control** :: [Type => Number or list/tuple]    ``Arg = [w|l_min/l_max]``

    Length controls. Project only those points whose p coordinate is within l_min < p < l_max.
    (http://docs.generic-mapping-tools.org/latest/project.html#l)
- **N** | **flat_earth** :: [Type => Bool or []]

    Flat Earth. Make a Cartesian coordinate transformation in the plane. [Default uses spherical trigonometry.]
    (http://docs.generic-mapping-tools.org/latest/project.html#n)
- **Q** | **km** :: [Type => Bool or []]

    Map type units.
    (http://docs.generic-mapping-tools.org/latest/project.html#q)
- **S** | **sort** :: [Type => Bool or []]

    Sort the output into increasing p order. Useful when projecting random data into a sequential profile.
    (http://docs.generic-mapping-tools.org/latest/project.html#s)
- **T** | **pole** :: [Type => list/tuple]    ``Arg = (px,py)``

    px,py sets the position of the rotation pole of the projection. (Definition 3).
    (http://docs.generic-mapping-tools.org/latest/project.html#t)
- **W** | **width_control** :: [Type => list/tuple]    ``Arg = (w_min,w_max)``

    Width controls. Project only those points whose q coordinate is within w_min < q < w_max.
    (http://docs.generic-mapping-tools.org/latest/project.html#w)

- $(GMT.opt_write)
- $(GMT.opt_append)
- $(GMT.opt_V)
- $(GMT.opt_b)
- $(GMT.opt_d)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_s)
- $(GMT.opt_swap_xy)
"""
function project(cmd0::String="", arg1=nothing; kwargs...)

	length(kwargs) == 0 && return monolitic("project", cmd0, arg1)

	d = KW(kwargs)
	cmd = parse_common_opts(d, "", [:V_params :b :d :e :f :g :h :i :s :yx])
	cmd = parse_these_opts(cmd, d, [[:A :azim], [:C :origin :start_pt], [:E :end_pt], [:F :out_flags], [:G :step :small_circle :no_input], [:L :length_control], [:N :flat_earth], [:Q :km], [:S :sort], [:T :pole],[:W :width_control]])

	if (!occursin("-G", cmd))
		cmd, got_fname, arg1 = find_data(d, cmd0, cmd, arg1)
	end
	common_grd(d, "project " * cmd, arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
project(arg1, cmd0::String=""; kw...) = project(cmd0, arg1; kw...)


text:
export conjugates_init, is_constant, is_squarefree, conjugates, angle, cos,
       sin, abs, abs2, sqrt

function is_constant(f::PolyElem)
  return f.length<2
end

function is_squarefree(f::PolyElem)
  return is_constant(gcd(f, derivative(f)))
end

function conjugates_init(f::Union{fmpz_poly, fmpq_poly})
  if typeof(f) == fmpq_poly
    f = f*den(f)
    g = Array(fmpz, length(f))
    for i = 1:f.length
      g[i] = ZZ(num(coeff(f, i-1)))
    end
    g = PolynomialRing(ZZ, string(var(parent(f))))[1](g)
    f = g
  end
  is_constant(gcd(f, derivative(f))) || error("poly should be square-free")
  c = roots_ctx()
  c.f = f
  r = _roots(f, 100)

  r_d = Array(BigComplex, 0)
  c_d = Array(BigComplex, 0)
  for i = 1:length(r)
    rr = BigComplex(r[i])
    if Base.abs2(imag(rr)) < 1e-20
      push!(r_d, rr)
      continue
    end
    if imag(rr) > 0
      push!(c_d, rr)
      continue
    end
  end
  assert(length(r_d) + 2*length(c_d) == length(r))
  c.r1 = length(r_d)
  c.r2 = length(c_d)
  sort!(r_d, lt = function(a,b) return real(a) < real(b); end)
  sort!(c_d, lt = function(a,b) return angle(a) < angle(b); end)
  c.r_d = vcat(r_d, c_d)

  c.r = Array(BigComplex, 0)

  old = precision(BigFloat)
  setprecision(53)
  for i = 1:length(c.r_d)
    push!(c.r, c.r_d[i])
  end
  setprecision(old)

  return c
end

function evaluate(f::fmpq_poly, r::BigComplex)
  #Horner - not elegant, but workable
  l = f.length-1
  s = BigComplex(BigFloat(coeff(f, l)))
  for i =l-1:-1:0
    s = s*r+BigComplex(BigFloat(coeff(f, i)))
  end
  return s
end

function evaluate{T <: RingElem}(f::fmpq_poly, r::T)
  R = parent(r)
  if iszero(f)
    return zero(R)
  end
  l = length(f) - 1
  s = R(coeff(f, l))
  for i in l-1:-1:0
    s = s*r + R(coeff(f, i))
  end
  return s
end

function evaluate(f::fmpz_poly, r::BigComplex)
  #Horner - not elegant, but workable
  l = f.length-1
  s = BigComplex(coeff(f, l))
  for i =l-1:-1:0
    s = s*r+BigComplex(coeff(f, i))
  end
  return s
end

function hensel_lift(f::fmpz_poly, r::BigComplex)
  return r - evaluate(f, r)/evaluate(derivative(f), r)
end

function conjugates(c::roots_ctx, p::Int)
  prec = precision(c.r[1])
  while prec < 2*p
    prec *= 2
    for i = 1:length(c.r)
      function _h()
        return hensel_lift(c.f, c.r[i])
      end
      c.r[i] = setprecision(_h, prec)
    end
  end
  return set_precision(c.r, p)
end

function set_precision(a::BigFloat, p::Int)
  return setprecision(function() return a*1.0; end, p)
end

function set_precision(a::BigComplex, p::Int)
  return setprecision(function() return BigComplex(a.re*1.0, a.im*1.0); end, p)
end

function set_precision(a::Array{BigComplex, 1}, p::Int)
  b = Array(BigComplex, 0);
  for i = 1:length(a)
    push!(b, set_precision(a[i], p))
  end
  return b
end

function minkowski(c::roots_ctx, a::nf_elem, p::Int)
  r = conjugates(c, p)
  m = Array(BigFloat, 0);
  a = (parent(a).pol.parent)(a)
  s2 = sqrt(BigFloat(2.0))
  for i = 1:c.r1
    push!(m, real(evaluate(a, r[i])))
  end
  for i = 1:c.r2
    z = evaluate(a, r[i+c.r1])
    push!(m, s2*real(z))
    push!(m, s2*imag(z))
  end
  return m
end

function length(c::roots_ctx, a::nf_elem, p::Int = 50)
  m = minkowski(c, a, p)
  return sum([x*x for x in m])
end

function minkowski_mat(c::roots_ctx, p::Int)
  if isdefined(c, :minkowski_mat) && c.minkowski_mat_p == p
    return c.minkowski_mat
  end
  old = precision(BigFloat)
  setprecision(p)
  r = conjugates(c, p)
  d = Array(typeof(r[1]), length(r))
  one = BigComplex(BigFloat(1.0), BigFloat(0.0))
  s2 = Base.sqrt(BigFloat(2.0))
  for i = 1:length(r)
    d[i] = one
  end

  n = degree(c.f)
  m = Array(BigFloat, n, n)

  for i = 1:n
    for j = 1:c.r1
      m[i,j] = real(d[j])
      if i < n
        d[j] = d[j] * r[j]
      end
    end
    for j=1:c.r2
      m[i, 2*j+c.r1-1] = s2*real(d[j+c.r1])
      m[i, 2*j+c.r1] = s2*imag(d[j+c.r1])
      if i < n
        d[j+c.r1] = d[j+c.r1]*r[j+c.r1]
      end
    end
  end
  setprecision(old)
  c.minkowski_mat = m
  c.minkowski_mat_p = p
  return m
end

function minkowski_mat(c::roots_ctx, K::AnticNumberField, p::Int = 50)
  return minkowski_mat(c, p)
end

function mult!(c::Array{BigFloat, 2}, a::fmpz_mat, b::Array{BigFloat, 2})
  s = Base.size(b)
  t = Base.size(c)
  cols(a) == s[1] || error("dimensions do not match")
  rows(a) == t[1] || error("dimensions do not match")
  t[2] == s[2]    || error("dimensions do not match")

  R = RealRing()
  tmp_mpz = R.z1
  tmp_mpz_r = R.t1
  tmp_mpfr = R.t2

  ##CF: careful: one SHOULD use mpfr_mul_z, but this converts to
  ##    mpfr every time. I think mpfr_mul_z should be re-written
  ##    to directly operate. Hoever, this is not going to happen soon.
  for i = 1:rows(a)
    for j = 1:t[2]
      if isdefined(c, i, j)
        ccall((:mpfr_set_zero, :libmpfr), Void, (Ptr{BigFloat}, Int), &c[i,j], 0)
      else
        c[i,j] = BigFloat(0)
      end

    end
    for j = 1:cols(a)
      z = ccall((:fmpz_mat_entry, :libflint), Ptr{fmpz},
        (Ptr{fmpz_mat}, Int, Int), &a, i - 1, j - 1)
      ccall((:flint_mpz_init_set_readonly, :libflint),
        Void, (Ptr{BigInt}, Ptr{fmpz}), &tmp_mpz, z)
      ccall((:mpfr_set_z, :libmpfr), Void, (Ptr{BigFloat}, Ptr{BigInt}, Int32),
        &tmp_mpz_r, &tmp_mpz, __get_rounding_mode())
      for k = 1:s[2]
        ccall((:mpfr_mul, :libmpfr), Int, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &tmp_mpfr, &b[j,k], &tmp_mpz_r, Base.MPFR.ROUNDING_MODE[end])
        ccall((:mpfr_add, :libmpfr), Int, (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigFloat}, Int32), &c[i,k], &tmp_mpfr, &c[i,k], Base.MPFR.ROUNDING_MODE[end])
      end
#      ccall((:flint_mpz_clear_readonly, :libflint), Void, (Ptr{BigInt}), &tmp_mpz)
    end
  end
  return c
end

function *(a::fmpz_mat, b::Array{BigFloat, 2})
  s = Base.size(b)
  cols(a) == s[1] || error("dimensions do not match")

  c = Array(BigFloat, rows(a), s[2])
  return mult!(c, a, b)
end

for (s,f) in ((:trunc, Base.trunc), (:round, Base.round), (:ceil, Base.ceil), (:floor, Base.floor))
  @eval begin
    function ($s)(a::Array{BigFloat, 2})
      s = Base.size(a)
      m = MatrixSpace(ZZ, s[1], s[2])()
      for i = 1:s[1]
        for j = 1:s[2]
          m[i,j] = ZZ(BigInt(($f)(a[i,j])))
        end
      end
      return m
    end
  end
end


text:
function destroy(args::Dict)
  file_type = args["<generator>"]
  file_name = args["<name>"]

  destroy(file_type, file_name)
end

function destroy(file_type, file_name)

  if file_type == "init"
    rm("$(pwd())/config/initializers/$file_name.jl")
    return
  end

  if file_type == "task"
    rm("$(pwd())/lib/tasks/$file_name.jl")
    return
  end

  parent_folders = [
    "src",
    "test"
  ]

  for parent_folder in parent_folders
    suffix = ( parent_folder == "test" ) ? "_test" : ""
    file_path = "$(pwd())/$parent_folder/$(pluralize(file_type))/$file_name$suffix.jl"

    if !isfile(file_path) ; continue ; end
    rm(file_path)

    folder_paths = []

    nested_path, parsed_file_name = contains(file_name, "/") ?
      rsplit(file_name, "/"; limit=2) : ("", file_name)

    parsed_file_name = String(parsed_file_name)

    current_nested_path = ""
    for nested_ancestor in split(nested_path, "/")
      current_nested_path = "$current_nested_path/$nested_ancestor"
      folder_path = "$(pwd())/$parent_folder/$(pluralize(file_type))$current_nested_path"

      push!(folder_paths, folder_path)
    end

    reverse!(folder_paths)

    for nested_path in folder_paths
      if length(readdir(nested_path)) != 0
        break
      end
      rm(nested_path)
    end
  end

end


text:
using HTTP
include("../src/newick_subset_parser.jl")

function run_server()

    query_dict = Dict()

    HTTP.listen("0.0.0.0", 8089, reauseaddr=true) do request::HTTP.Request

        println("******************")
        @show request
        println("******************")

        uri = parse(HTTP.URI, request.target)
        query_dict = HTTP.queryparams(uri)

        headers = Dict{AbstractString,AbstractString}(
            "Server"            => "Julia/$VERSION",
            "Content-Type"      => "text/html; charset=utf-8",
            "Content-Language"  => "en",
            "Date"              => Dates.format(now(Dates.UTC), Dates.RFC1123Format),
            "Access-Control-Allow-Origin" => "*" )

        return HTTP.Response(200, HTTP.Headers(collect(headers)), body = main(nwk, query_dict["id"], query_dict["level"], dname_dict))

    end
end

run_server()

max_stars_repo_name: bcbi/PhytoKB-API

max_issues_repo_name: bcbi/PhytoKB-API

max_forks_repo_name: bcbi/PhytoKB-API