text:
"""
`barplot(text, heights; nargs...)` → `Plot`

Description
============

Draws a horizontal barplot. It uses the first parameter (`text`) to denote
the names for the bars, and the second parameter (`heights`) as their values.
This means that the two vectors have to have the same length.
Alternatively, one can specify a barplot using a dictionary.
In that case, the keys will be used as the names and the values,
which have to be numeric, will be used as the heights of the bars.

Usage
======

    barplot(text, heights; border = :solid, title = "", margin = 3, padding = 1, color = :blue, width = 40, labels = true, symb = "▪")

    barplot(dictionary; nargs...)

Arguments
==========

- **`text`** : The labels/captions of the bars

- **`heights`** : The values/heights of the bars

- **`dictionary`** : A dictonary in which the keys will be used as `text`
and the values will be utilized as `heights`.

- **`border`** : The style of the bounding box of the plot.
Supports `:solid`, `:bold`, `:dashed`, `:dotted`, `:ascii`, and `:none`.

- **`title`** : Text to display on the top of the plot.

- **`margin`** : Number of empty characters to the left of the whole plot.

- **`padding`** : Space of the left and right of the plot between the labels and the canvas.

- **`color`** : Color of the drawing. Can be any of `:blue`, `:red`, `:yellow`

- **`width`** : Number of characters per row that should be used for plotting.

- **`labels`** : Can be used to hide the labels by setting `labels=false`.

- **`symb`** : Specifies the character that should be used to render the bars

Returns
========

A plot object of type `Plot{BarplotGraphics}`

Author(s)
==========

- Christof Stocker (Github: https://github.com/Evizero)

Examples
=========

    julia> barplot(["Paris", "New York", "Moskau", "Madrid"],
                   [2.244, 8.406, 11.92, 3.165],
                   title = "Population")

                            Population
             ┌────────────────────────────────────────┐
       Paris │▪▪▪▪▪▪ 2.244                            │
    New York │▪▪▪▪▪▪▪▪▪▪▪▪▪▪▪▪▪▪▪▪▪▪▪ 8.406           │
      Moskau │▪▪▪▪▪▪▪▪▪▪▪▪▪▪▪▪▪▪▪▪▪▪▪▪▪▪▪▪▪▪▪▪▪ 11.92 │
      Madrid │▪▪▪▪▪▪▪▪▪ 3.165                         │
             └────────────────────────────────────────┘

See also
=========

`Plot`, `histogram`, `BarplotGraphics`
"""
function barplot{T<:AbstractString,N<:Real}(
        text::Vector{T}, heights::Vector{N};
        border = :solid,
        title::AbstractString = "",
        margin::Int = 3,
        padding::Int = 1,
        color::Symbol = :blue,
        width::Int = 40,
        labels::Bool = true,
        symb = "▪")
    margin >= 0 || throw(ArgumentError("Margin must be greater than or equal to 0"))
    length(text) == length(heights) || throw(DimensionMismatch("The given vectors must be of the same length"))
    minimum(heights) >= 0 || throw(ArgumentError("All values have to be positive. Negative bars are not supported."))
    width = max(width, 5)

    area = BarplotGraphics(heights, width, color = color, symb = symb)
    new_plot = Plot(area, title = title, margin = margin,
                   padding = padding, border = border, labels = labels)
    for i in 1:length(text)
        annotate!(new_plot, :l, i, text[i])
    end
    new_plot
end

"""
`barplot!(plot, text, heights; nargs)` → `Plot`

Mutating variant of `barplot`, in which the first parameter (`plot`) specifies
the existing plot to draw on.

See `barplot` for more information.
"""
function barplot!{C<:BarplotGraphics,T<:AbstractString,N<:Real}(
        plot::Plot{C},
        text::Vector{T},
        heights::Vector{N};
        args...)
    length(text) == length(heights) || throw(DimensionMismatch("The given vectors must be of the same length"))
    !isempty(text)|| throw(ArgumentError("Can't append empty array to barplot"))
    curidx = nrows(plot.graphics)
    addrow!(plot.graphics, heights)
    for i = 1:length(heights)
        annotate!(plot, :l, curidx + i, text[i])
    end
    plot
end

function barplot!{C<:BarplotGraphics,T<:AbstractString,N<:Real}(
        plot::Plot{C},
        text::T,
        heights::N;
        args...)
    text == "" && throw(ArgumentError("Can't append empty array to barplot"))
    curidx = nrows(plot.graphics)
    addrow!(plot.graphics, heights)
    annotate!(plot, :l, curidx + 1, text)
    plot
end

function barplot{T,N<:Real}(dict::Dict{T,N}; args...)
    barplot(collect(keys(dict)), collect(values(dict)); args...)
end

function barplot{T,N<:Real}(labels::Vector{T}, heights::Vector{N}; args...)
    labels_str = map(string, labels)
    barplot(labels_str, heights; args...)
end


text:
module ATSMod

using Pkg.Artifacts

# All this module will do is reference its `arty` Artifact.
arty = artifact"arty"

function do_test()
    return isfile(joinpath(arty, "bin", "socrates"))
end

end # module

ext: jl


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule MDAL_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("MDAL")
JLLWrappers.@generate_main_file("MDAL", UUID("ba355c7d-02a9-511b-bb8a-ec390ce7a57a"))
end  # module MDAL_jll


text:
using Clustering
using GaussianMixtures
using LinearAlgebra

mutable struct MOE{X,Y,L,U,S,K,M,V,W} <: AbstractSurrogate
    x::X
    y::Y
    lb::L
    ub::U
    local_surr::S
    k::K
    means::M
    varcov::V
    weights::W
end

#Radial structure:
function RadialBasisStructure(;radial_function,scale_factor,sparse)
    return (name = "RadialBasis", radial_function = radial_function, scale_factor = scale_factor, sparse = sparse)
end

#Kriging structure:
function KrigingStructure(;p,theta)
    return (name = "Kriging", p = p, theta = theta)
end

#Linear structure
function LinearStructure()
    return (name = "LinearSurrogate")
end

#InverseDistance structure
function InverseDistanceStructure(;p)
    return (name = "InverseDistanceSurrogate", p = p)
end

#Lobachesky structure
function LobacheskyStructure(;alpha,n,sparse)
    return (name = "LobacheskySurrogate", alpha = alpha, n = n, sparse = sparse)
end

#Neural structure
function NeuralStructure(;model,loss,opt,n_echos)
    return (name ="NeuralSurrogate", model = model ,loss = loss,opt = opt,n_echos = n_echos)
end

#Random forest structure
function RandomForestStructure(;num_round)
    return (name = "RandomForestSurrogate", num_round = num_round)
end

#Second order poly structure
function SecondOrderPolynomialStructure()
    return (name = "SecondOrderPolynomialSurrogate")
end

#Wendland structure
function WendlandStructure(; eps, maxiters, tol)
    return (name = "Wendland", eps = eps, maxiters = maxiters, tol = tol)
end

#Polychaos structure
function PolyChaosStructure(; op)
    return (name = "PolynomialChaosSurrogate", op = op)
end

function MOE(x,y,lb::Number,ub::Number; k::Int = 2, local_kind = [RadialBasisStructure(radial_function = linearRadial, scale_factor=1.0,sparse = false),RadialBasisStructure(radial_function = cubicRadial, scale_factor=1.0, sparse = false)])
    if k != length(local_kind)
        throw("Number of mixtures = $k is not equal to length of local surrogates")
    end
    n = length(x)
    # find weight, mean and variance for each mixture
    # For GaussianMixtures I need nxd matrix
    X_G = reshape(x,(n,1))
    moe_gmm = GaussianMixtures.GMM(k,X_G)
    weights = GaussianMixtures.weights(moe_gmm)
    means = GaussianMixtures.means(moe_gmm)
    variances = moe_gmm.Σ

    #cluster the points
    #For clustering I need dxn matrix
    X_C = reshape(x,(1,n))
    KNN = kmeans(X_C, k)
    x_c = [ Array{eltype(x)}(undef,0) for i = 1:k]
    y_c = [ Array{eltype(y)}(undef,0) for i = 1:k]
    a = assignments(KNN)
    @inbounds for i = 1:n
        pos = a[i]
        append!(x_c[pos],x[i])
        append!(y_c[pos],y[i])
    end

    local_surr = Dict()
    for i = 1:k
        if local_kind[i][1] == "RadialBasis"
            #fit and append to local_surr
            my_local_i = RadialBasis(x_c[i],y_c[i],lb,ub,rad = local_kind[i].radial_function, scale_factor = local_kind[i].scale_factor, sparse = local_kind[i].sparse)
            local_surr[i] = my_local_i

        elseif local_kind[i][1] == "Kriging"
            my_local_i = Kriging(x_c[i], y_c[i],lb,ub, p = local_kind[i].p, theta = local_kind[i].theta)
            local_surr[i] = my_local_i

        elseif local_kind[i] == "LinearSurrogate"
            my_local_i = LinearSurrogate(x_c[i], y_c[i],lb,ub)
            local_surr[i] = my_local_i

        elseif local_kind[i][1] == "InverseDistanceSurrogate"
            my_local_i = InverseDistanceSurrogate(x_c[i], y_c[i],lb,ub, local_kind[i].p)
            local_surr[i] = my_local_i

        elseif local_kind[i][1] == "LobacheskySurrogate"
            my_local_i = LobacheskySurrogate(x_c[i], y_c[i],lb,ub,alpha = local_kind[i].alpha , n = local_kind[i].n, sparse = local_kind[i].sparse)
            local_surr[i] = my_local_i

        elseif local_kind[i][1] == "NeuralSurrogate"
            my_local_i = NeuralSurrogate(x_c[i], y_c[i],lb,ub, model = local_kind[i].model , loss = local_kind[i].loss ,opt = local_kind[i].opt ,n_echos = local_kind[i].n_echos)
            local_surr[i] = my_local_i

        elseif local_kind[i][1] == "RandomForestSurrogate"
            my_local_i = RandomForestSurrogate(x_c[i], y_c[i],lb,ub, num_round = local_kind[i].num_round)
            local_surr[i] = my_local_i

        elseif local_kind[i] == "SecondOrderPolynomialSurrogate"
            my_local_i = SecondOrderPolynomialSurrogate(x_c[i], y_c[i],lb,ub)
            local_surr[i] = my_local_i

        elseif local_kind[i][1] == "Wendland"
            my_local_i = Wendand(x_c[i], y_c[i],lb,ub, eps = local_kind[i].eps, maxiters = local_kind[i].maxiters, tol = local_kind[i].tol)
            local_surr[i] = my_local_i
        elseif local_kind[i][1] == "PolynomialChaosSurrogate"
            my_local_i = PolynomialChaosSurrogate(x,y,lb,ub, op = local_kind[i].op)
            local_surr[i] = my_local_i

        else
            throw("A surrogate with name provided does not exist or is not currently supported with MOE.")
        end
    end
    return MOE(x,y,lb,ub,local_surr,k,means,variances,weights)
end

function MOE(x,y,lb,ub; k::Int = 2,
            local_kind = [RadialBasisStructure(radial_function = linearRadial, scale_factor=1.0, sparse = false),RadialBasisStructure(radial_function = cubicRadial, scale_factor=1.0, sparse = false)])

    n = length(x)
    d = length(lb)
    #GMM parameters:
    X_G = collect(reshape(collect(Base.Iterators.flatten(x)), (d,n))')
    my_gmm = GMM(k,X_G,kind = :full)
    weights = my_gmm.w
    means = my_gmm.μ
    varcov = my_gmm.Σ

    #cluster the points
    X_C = collect(reshape(collect(Base.Iterators.flatten(x)), (d,n)))
    KNN = kmeans(X_C, k)
    x_c = [ Array{eltype(x)}(undef,0) for i = 1:k]
    y_c = [ Array{eltype(y)}(undef,0) for i = 1:k]
    a = assignments(KNN)
    @inbounds for i = 1:n
        pos = a[i]
        x_c[pos] = vcat(x_c[pos],x[i])
        append!(y_c[pos],y[i])
    end

    local_surr = Dict()
    for i = 1:k
        if local_kind[i][1] == "RadialBasis"
            #fit and append to local_surr
            my_local_i = RadialBasis(x_c[i],y_c[i],lb,ub,rad = local_kind[i].radial_function, scale_factor = local_kind[i].scale_factor, sparse = local_kind[i].sparse)
            local_surr[i] = my_local_i

        elseif local_kind[i][1] == "Kriging"
            my_local_i = Kriging(x_c[i], y_c[i],lb,ub, p = local_kind[i].p, theta = local_kind[i].theta)
            local_surr[i] = my_local_i

        elseif local_kind[i] == "LinearSurrogate"
            my_local_i = LinearSurrogate(x_c[i], y_c[i],lb,ub)
            local_surr[i] = my_local_i

        elseif local_kind[i][1] == "InverseDistanceSurrogate"
            my_local_i = InverseDistanceSurrogate(x_c[i], y_c[i],lb,ub, local_kind[i].p)
            local_surr[i] = my_local_i

        elseif local_kind[i][1] == "LobacheskySurrogate"
            my_local_i = LobacheskySurrogate(x_c[i], y_c[i],lb,ub,alpha = local_kind[i].alpha , n = local_kind[i].n, sparse = local_kind[i].sparse)
            local_surr[i] = my_local_i

        elseif local_kind[i][1] == "NeuralSurrogate"
            my_local_i = NeuralSurrogate(x_c[i], y_c[i],lb,ub, model = local_kind[i].model , loss = local_kind[i].loss ,opt = local_kind[i].opt ,n_echos = local_kind[i].n_echos)
            local_surr[i] = my_local_i

        elseif local_kind[i][1] == "RandomForestSurrogate"
            my_local_i = RandomForestSurrogate(x_c[i], y_c[i],lb,ub, num_round = local_kind[i].num_round)
            local_surr[i] = my_local_i

        elseif local_kind[i] == "SecondOrderPolynomialSurrogate"
            my_local_i = SecondOrderPolynomialSurrogate(x_c[i], y_c[i],lb,ub)
            local_surr[i] = my_local_i

        elseif local_kind[i][1] == "Wendland"
            my_local_i = Wendand(x_c[i], y_c[i],lb,ub, eps = local_kind[i].eps, maxiters = local_kind[i].maxiters, tol = local_kind[i].tol)
            local_surr[i] = my_local_i

        elseif local_kind[i][1] == "PolynomialChaosSurrogate"
            my_local_i = PolynomialChaosSurrogate(x,y,lb,ub, op = local_kind[i].op)
            local_surr[i] = my_local_i
        else
            throw("A surrogate with name "* local_kind[i][1] *" does not exist or is not currently supported with MOE.")
        end
    end
    return MOE(x,y,lb,ub,local_surr,k,means,varcov,weights)
end

function _prob_x_in_i(x::Number,i,mu,varcov,alpha,k)
    num = (1/sqrt(varcov[i]))*alpha[i]*exp(-0.5(x-mu[i])*(1/varcov[i])*(x-mu[i]))
    den = sum([(1/sqrt(varcov[j]))*alpha[j]*exp(-0.5(x-mu[j])*(1/varcov[j])*(x-mu[j]))  for j = 1:k])
    return num/den
end

function _prob_x_in_i(x,i,mu,varcov,alpha,k)
    num = (1/sqrt(det(varcov[i])))*alpha[i]*exp(-0.5*(x .- mu[i,:])'*(inv(varcov[i]))*(x .- mu[i,:]))
    den = sum([(1/sqrt(det(varcov[j])))*alpha[j]*exp(-0.5*(x .- mu[j,:])'*(inv(varcov[j]))*(x .- mu[j,:])) for j = 1:k])
    return num/den
end

function (moe::MOE)(val)
    return prod([moe.local_surr[i](val)*_prob_x_in_i(val,i,moe.means,moe.varcov,moe.weights,moe.k) for i = 1:moe.k])
end

function add_point!(my_moe::MOE,x_new,y_new)
    if length(my_moe.x[1]) == 1
        #1D
        my_moe.x = vcat(my_moe.x,x_new)
        my_moe.y = vcat(my_moe.y,y_new)
        n = length(my_moe.x)

        #New mixture parameters
        X_G = reshape(my_moe.x,(n,1))
        moe_gmm = GaussianMixtures.GMM(my_moe.k,X_G)
        my_moe.weights = GaussianMixtures.weights(moe_gmm)
        my_moe.means = GaussianMixtures.means(moe_gmm)
        my_moe.varcov = moe_gmm.Σ

        #Find cluster of new point(s):
        n_added = length(x_new)
        X_C = reshape(my_moe.x,(1,n))
        KNN = kmeans(X_C, my_moe.k)
        a = assignments(KNN)
        #Recalculate only relevant surrogates
        for i = 1:n_added
            pos = a[n-n_added+i]
            add_point!(my_moe.local_surr[i],my_moe.x[n-n_added+i],my_moe.y[n-n_added+i])
        end
    else
        #ND
        my_moe.x = vcat(my_moe.x,x_new)
        my_moe.y = vcat(my_moe.y,y_new)
        n = length(my_moe.x)
        d = length(my_moe.lb)
        #New mixture parameters
        X_G = collect(reshape(collect(Base.Iterators.flatten(my_moe.x)), (d,n))')
        my_gmm = GMM(my_moe.k,X_G,kind = :full)
        my_moe.weights = my_gmm.w
        my_moe.means = my_gmm.μ
        my_moe.varcov = my_gmm.Σ

        #cluster the points
        X_C = collect(reshape(collect(Base.Iterators.flatten(my_moe.x)), (d,n)))
        KNN = kmeans(X_C, my_moe.k)
        a = assignments(KNN)
        n_added = length(x_new)
        for i = 1:n_added
            pos = a[n-n_added+i]
            add_point!(my_moe.local_surr[i],my_moe.x[n-n_added+i],my_moe.y[n-n_added+i])
        end
    end
    nothing
end


text:
using FluxionArrays
using FluxionFields
using Base.Test

eval_get(args...; kwds...) = from_backend(
    evaluate_to_array(default_array_backend(), args...; kwds...))

@testset "Fields" begin

@testset "stage1" begin
    # That's what happens at the start of the integration when the initial condition
    # has to be evaluated into an array.

    x = BasicDimension(RegularGrid(0, 1, 20))
    y = BasicDimension(RegularGrid(0, 1, 30))

    x_arr = vertices(grid(x))
    y_arr = vertices(grid(y))

    f = x + y + 1

    # Then there would be some call like `integrate(eqn, f, ...)`
    # and inside it we will have to obtain the initial array representing `f`.

    f_arr = from_backend(evaluate(f, [x, y]).data)
    @test f_arr == x_arr .+ y_arr' .+ 1

    f_arr = eval_get(f, [x, y])
    @test f_arr == x_arr .+ y_arr' .+ 1

end

@testset "stage2" begin

    # That's what happens inside the integration.
    # An array function is prepared and then evaluated on each step.

    x = BasicDimension(RegularGrid(0, 1, 20))
    y = BasicDimension(RegularGrid(0, 1, 30))
    t = UnknownDimension()
    f = UnknownField([t, x, y])

    f_initial = x + y + 1
    right_side = f + t + x + 2

    # Variant 1: we treat known and unknown dimensions differently,
    # so f(t, x, y) is understood as an unknown array of two (known) dimensions, x and y.
    #
    # right_side_f = as_array_function(right_side, [f, t], [x, y])
    #
    # We need to ensure that the return value of this call can be supplied to the next call.
    # Currently it seems that the order of known dimensions in `f` is fixed
    # (since it is a `Field`, not `FieldExpr`), so all the information is available.
    #
    # What if we want to make a function that takes an array of `t` and a corresponding known
    # field `f(t, x, y)`? How do we distinguish between these cases? Do we actually need this case?

    # Variant 2: known and unknown dimensions are the same. E.g. an unknown field
    # for finding a stationary state of a heat equation has no unknown dimensions.
    # So every solver will do its own postprocessing of the equation,
    # replacing the field if necessary.
    f_step = remove_unknown_dimensions(f)
    @test known_dimensions(f_step) == [x, y]
    @test unknown_dimensions(f_step) == []
    right_side_step = replace_expr(right_side, f => f_step)
    # Do we automatically assume that an unknown dimension is a scalar?
    right_side_f = compile_wrapped_array_function(
        default_array_backend(),
        wrapped_array_function([right_side_step => [x, y]], [f_step, t]))

    # Then we can call the function as
    f_arr = eval_get(f_initial, [x, y])
    t = 3
    new_f_arr = from_backend(right_side_f(f_arr, t))

    x_arr = vertices(grid(x))
    y_arr = vertices(grid(y))

    @test new_f_arr == x_arr .+ y_arr' .+ 1 .+ t .+ x_arr .+ 2

end

@testset "stage3" begin

    # Finally, this is what's happening inside a sampler
    # We have an array that we need to wrap in a field
    # and present to the user for him to do operations on

    x = BasicDimension(RegularGrid(0, 1, 20))
    y = BasicDimension(RegularGrid(0, 1, 30))
    t = BasicDimension(RegularGrid(4, 5, 10))

    x_arr = vertices(grid(x))
    y_arr = vertices(grid(y))
    t_arr = vertices(grid(t))

    f_arr = eval_get(x + y + 1, [x, y])

    f = known_field(to_backend(f_arr), [x, y]) # this is passed to the user

    results = Array{KnownField, 1}()

    # defined by the user
    sampler(f_val, t_val) = f_val + x + t_val

    # constructed by the integrator
    f_unknown = UnknownField([x, y])
    t_unknown = UnknownDimension()
    sampler_expr = sampler(f_unknown, t_unknown)

    # The function will be a part of Fluxion, it's too high-level for this package
    # observable_dims = ordered_dimensions(sampler_expr, [x, y])
    observable_dims = [x, y]

    sampler_func = field_function([sampler_expr => observable_dims], [f_unknown, t_unknown])
    sampler_func = compile_field_function(default_array_backend(), sampler_func)

    # The integration loop
    for t_val in vertices(grid(t))
        observable = sampler_func(f, t_val)
        push!(results, observable)
    end

    samples = join_fields(results, t)

    @test known_dimensions(samples) == [x, y, t]

    samples_data = eval_get(samples, [x, y, t])
    reference_data = (
        reshape(x_arr, (length(x_arr), 1, 1)) .* 2
        .+ reshape(y_arr, (1, length(y_arr), 1))
        .+ 1
        .+ reshape(t_arr, (1, 1, length(t_arr))))

    @test isapprox(samples_data, reference_data)

end

#=
@testset "regular randoms" begin

    # Typical example: adding Wigner noise
    x = BasicDimension(RegularGrid(-5, 5, 64))
    classic_state = 1 / cosh(10 - x)

    k = fourier_space(x)
    classic_state_k = fourier(classic_state, x => k)
    wigner_state_k = classic_state_k + FluxionFields.random_normal([k], 0, 1/2) # TODO: should be complex
    wigner_state = fourier(wigner_state_k, k => x, inverse=true)

    # in the integrate()
    # we have the number of trajectories and a seed that we need to use to instantiate the array

    # Trajectories are passed to sample(), because there may be an expression which has
    # different number of trajectories in subtrees, e.g.
    # mean(sample(expr1, trajectories=100)) + mean(sample(expr2, trajectories=200))
    # But they also can be passed to as_array()

    # Seed is passed to as_array, because we will have to create a random state internally
    # to sample the randoms, and creating random states inside the generated function
    # will be technically difficult.

    wigner_state_sampled = sample(wigner_state, 100)
    arr1 = as_array(wigner_state_sampled, [x, STOCHASTIC_DIMENSION], seed=123)
    @test size(arr1) == (64, 100)

    arr2 = as_array(wigner_state_sampled, [x, STOCHASTIC_DIMENSION], seed=123) # same array
    @test isapprox(arr1, arr2)

    arr3 = as_array(wigner_state_sampled, [x, STOCHASTIC_DIMENSION]) # random seed used
    @test !isapprox(arr1, arr3)

    arr4 = as_array(wigner_state_sampled, [STOCHASTIC_DIMENSION, x]) # random seed used
    @test size(arr4) == (100, 64)

end

@testset "distinguishability of random variables" begin
    x = BasicDimension(RegularGrid(-5, 5, 32))
    rv1 = FluxionFields.random_normal([x], 0, 1/2)
    rv2 = FluxionFields.random_normal([x], 0, 1/2)

    arr1 = as_array(sample(rv1 - rv1, 10), [x, STOCHASTIC_DIMENSION])
    @test isapprox(arr1, zeros(arr1))

    arr2 = as_array(sample(rv1 - rv2, 10), [x, STOCHASTIC_DIMENSION])
    @test !isapprox(arr2, zeros(arr2))

    # Same RV only evaluate to the same random arrays only inside a single sample() call
    arr1 = as_array(sample(rv1, 10) - sample(rv1, 10), [x, STOCHASTIC_DIMENSION])
    @test !isapprox(arr1, zeros(arr1))

end

@testset "returned shape" begin
    x = BasicDimension(RegularGrid(-5, 5, 64))
    y = BasicDimension(RegularGrid(-5, 5, 32))

    @test size(as_array(x, [x, y])) == (64, 32)
    @test size(as_array(1, [x, y])) == (64, 32)

    dims = [x, y, STOCHASTIC_DIMENSION]
    @test size(as_array(sample(1, 10), dims)) == (64, 32, 10)
    @test size(as_array(sample(x, 10), dims)) == (64, 32, 10)
    @test size(as_array(sample(FluxionFields.random_normal([x], 0, 1), 10), dims)) == (64, 32, 10)
    @test size(as_array(sample(y + FluxionFields.random_normal([x], 0, 1), 10), dims)) == (64, 32, 10)
end
=#

#=
@testset "Wiener process" begin
    x = BasicDimension("x", RegularGrid(0, 1, 20))
    y = BasicDimension("y", RegularGrid(0, 1, 30))
    t = UnknownDimension("t")

    z = wiener_normal(x, t)

    func = as_array_function(x + y + z)
    func(x_arr, y_arr, dt, seed=123)

    # Need to handle:
    # - varying dt
    # - single/double step or some kind of noise splitting (check Num.Rec.)
    # - persistent state (for performance reasons)

end
=#

end


text:
export Random, CSV, Dates, DataFrames, Statistics
#export read, write
export DateTime, Date,Dates, Day,Year, Month, Week, Minute, Hour, year, day, month, hour, minute, week, dayofweek, dayofquarter
export DataFrame, DataFrameRow, DataFrameRows, DataFrameColumns, nrow, ncol, eachcol
export seed!
export mean, median, middle, quantile, std, var, cor, cov

module PkgDeps

import CSV: read, write
using Dates: DateTime, Date,Dates, Day,Year, Month, Week, Minute, Hour, year, day, month, hour, minute, week, dayofweek, dayofquarter
using DataFrames: DataFrame, DataFrameRow, DataFrameRows, DataFrameColumns, nrow, ncol, eachcol
using Random: seed!
using Statistics: mean, median, middle, quantile, std, var, cor, cov
import Random, CSV, Dates, DataFrames, Statistics

export Random, CSV, Dates, DataFrames, Statistics
export read, write
export DateTime, Date,Dates, Day,Year, Month, Week, Minute, Hour, year, day, month, hour, minute, week, dayofweek, dayofquarter
export DataFrame, DataFrameRow, DataFrameRows, DataFrameColumns, nrow, ncol, eachcol
export seed!
export mean, median, middle, quantile, std, var, cor, cov

end


text:
using Compose
import Cairo, Fontconfig

imgs = [SVG("arrow.svg", 5cm, 5cm),
        PDF("arrow.pdf", 5cm, 5cm)]

X =  [0.047 0.87 0.95 0.93;
0.22 0.01 0.21 0.7;
0.86 0.85 0.95 0.21]

point_array = [[(x1,y1), (x2,y2)] for (x1,y1,x2,y2) in zip(X[:,1],X[:,2],X[:,3],X[:,4])]
img = compose(context(),
    (context(), line(point_array), stroke(["red","green","deepskyblue"]), arrow())
)

draw.(imgs, [img])


text:
"""
    volume(S::AbstractSet)

Return the volume of the given set.

### Input

- `S` -- set

### Output

A number.

### Notes

The `volume` refers to the volume of the set with respect to some measure.
The actual notion of `volume` adopted is left to the package implementing
this function.
"""
function volume(::AbstractSet) end

"""
    surface_area(S::AbstractSet)

Return the surface area of the given set.

### Input

- `S` -- set

### Output

A number.

### Notes

The `surface_area` refers to the surface area of the set with respect to some
measure. The actual notion of `surface` adopted is left to the package
implementing this function.
"""
function surface_area(::AbstractSet) end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module TestAddmul

using Base: rtoldefault
using Test
using LinearAlgebra
using LinearAlgebra: AbstractTriangular
using SparseArrays
using Random

_rand(::Type{T}) where {T <: AbstractFloat} = T(randn())
_rand(::Type{T}) where {F, T <: Complex{F}} = T(_rand(F), _rand(F))
_rand(::Type{T}) where {T <: Integer} =
    T(rand(max(typemin(T), -10):min(typemax(T), 10)))
_rand(::Type{BigInt}) = BigInt(_rand(Int))

function _rand(A::Type{<:Array}, shape)
    T = eltype(A)
    data = T[_rand(T) for _ in 1:prod(shape)]
    return copy(reshape(data, shape))
end

constructor_of(::Type{T}) where T = getfield(parentmodule(T), nameof(T))

function _rand(A::Type{<: AbstractArray}, shape)
    data = _rand(Array{eltype(A)}, shape)
    T = constructor_of(A)
    if A <: Union{Bidiagonal, Hermitian, Symmetric}
        return T(data, rand([:U, :L]))
        # Maybe test with both :U and :L?
    end
    return T(data)
end

_rand(A::Type{<: SymTridiagonal{T}}, shape) where {T} =
    SymTridiagonal(_rand(Symmetric{T}, shape))

const FloatOrC = Union{AbstractFloat, Complex{<: AbstractFloat}}
const IntegerOrC = Union{Integer, Complex{<: Integer}}
const LTri = Union{LowerTriangular, UnitLowerTriangular, Diagonal}
const UTri = Union{UpperTriangular, UnitUpperTriangular, Diagonal}

needsquare(::Type{<:Matrix}) = false
needsquare(::Type) = true

testdata = []

sizecandidates = 1:4
floattypes = [
    Float64, Float32, ComplexF64, ComplexF32,  # BlasFloat
    BigFloat,
]
inttypes = [
    Int,
    BigInt,
]
# `Bool` can be added to `inttypes` but it's hard to handle
# `InexactError` bug that is mentioned in:
# https://github.com/JuliaLang/julia/issues/30094#issuecomment-440175887
alleltypes = [floattypes; inttypes]
celtypes = [Float64, ComplexF64, BigFloat, Int]

mattypes = [
    Matrix,
    Bidiagonal,
    Diagonal,
    Hermitian,
    LowerTriangular,
    SymTridiagonal,
    Symmetric,
    Tridiagonal,
    UnitLowerTriangular,
    UnitUpperTriangular,
    UpperTriangular,
]

isnanfillable(::AbstractArray) = false
isnanfillable(::Array{<:AbstractFloat}) = true
isnanfillable(A::AbstractArray{<:AbstractFloat}) = parent(A) isa Array

"""
Sample `n` elements from `S` on average but make sure at least one
element is sampled.
"""
function sample(S, n::Real)
    length(S) <= n && return S
    xs = randsubseq(S, n / length(S))
    return length(xs) > 0 ? xs : rand(S, 1)  # sample at least one
end

function inputeltypes(celt, alleltypes = alleltypes)
    # Skip if destination type is "too small"
    celt <: Bool && return []
    filter(alleltypes) do aelt
        celt <: Real && aelt <: Complex && return false
        !(celt <: BigFloat) && aelt <: BigFloat && return false
        !(celt <: BigInt) && aelt <: BigInt && return false
        celt <: IntegerOrC && aelt <: FloatOrC && return false
        if celt <: IntegerOrC && !(celt <: BigInt)
            typemin(celt) > typemin(aelt) && return false
            typemax(celt) < typemax(aelt) && return false
        end
        return true
    end
end
# Note: using `randsubseq` instead of `rand` to avoid repetition.

function inputmattypes(cmat, mattypes = mattypes)
    # Skip if destination type is "too small"
    cmat <: Union{Bidiagonal, Tridiagonal, SymTridiagonal,
                  UnitLowerTriangular, UnitUpperTriangular,
                  Hermitian, Symmetric} && return []
    filter(mattypes) do amat
        cmat <: Diagonal && (amat <: Diagonal || return false)
        cmat <: LowerTriangular && (amat <: LTri || return false)
        cmat <: UpperTriangular && (amat <: UTri || return false)
        return true
    end
end

n_samples = 1.5
# n_samples = Inf  # to try all combinations
for cmat in mattypes,
    amat in sample(inputmattypes(cmat), n_samples),
    bmat in sample(inputmattypes(cmat), n_samples),
    celt in celtypes,
    aelt in sample(inputeltypes(celt), n_samples),
    belt in sample(inputeltypes(celt), n_samples)

    push!(testdata, (cmat{celt}, amat{aelt}, bmat{belt}))
end

@testset "mul!(::$TC, ::$TA, ::$TB, α, β)" for (TC, TA, TB) in testdata
    if needsquare(TA)
        na1 = na2 = rand(sizecandidates)
    else
        na1, na2 = rand(sizecandidates, 2)
    end
    if needsquare(TB)
        nb2 = na2
    elseif needsquare(TC)
        nb2 = na1
    else
        nb2 = rand(sizecandidates)
    end
    asize = (na1, na2)
    bsize = (na2, nb2)
    csize = (na1, nb2)

    @testset for α in Any[true, eltype(TC)(1), _rand(eltype(TC))],
                 β in Any[false, eltype(TC)(0), _rand(eltype(TC))]

        C = _rand(TC, csize)
        A = _rand(TA, asize)
        B = _rand(TB, bsize)

        # This is similar to how `isapprox` choose `rtol` (when
        # `atol=0`) but consider all number types involved:
        rtol = max(rtoldefault.(real.(eltype.((C, A, B))))...,
                   rtoldefault.(real.(typeof.((α, β))))...)

        Cc = copy(C)
        Ac = Matrix(A)
        Bc = Matrix(B)
        returned_mat = mul!(C, A, B, α, β)
        @test returned_mat === C
        @test collect(returned_mat) ≈ α * Ac * Bc + β * Cc  rtol=rtol

        y = C[:, 1]
        x = B[:, 1]
        yc = Vector(y)
        xc = Vector(x)
        returned_vec = mul!(y, A, x, α, β)
        @test returned_vec === y
        @test collect(returned_vec) ≈ α * Ac * xc + β * yc  rtol=rtol

        if TC <: Matrix
            @testset "adjoint and transpose" begin
                @testset for fa in [identity, adjoint, transpose],
                             fb in [identity, adjoint, transpose]
                    fa === fb === identity && continue

                    Af = fa === identity ? A : fa(_rand(TA, reverse(asize)))
                    Bf = fb === identity ? B : fb(_rand(TB, reverse(bsize)))

                    Ac = collect(Af)
                    Bc = collect(Bf)
                    Cc = collect(C)

                    returned_mat = mul!(C, Af, Bf, α, β)
                    @test returned_mat === C
                    @test collect(returned_mat) ≈ α * Ac * Bc + β * Cc  rtol=rtol
                end
            end
        end

        if isnanfillable(C)
            @testset "β = 0 ignores C .= NaN" begin
                parent(C) .= NaN
                Ac = Matrix(A)
                Bc = Matrix(B)
                returned_mat = mul!(C, A, B, α, zero(eltype(C)))
                @test returned_mat === C
                @test collect(returned_mat) ≈ α * Ac * Bc  rtol=rtol
            end
        end

        if isnanfillable(A)
            @testset "α = 0 ignores A .= NaN" begin
                parent(A) .= NaN
                Cc = copy(C)
                returned_mat = mul!(C, A, B, zero(eltype(A)), β)
                @test returned_mat === C
                @test collect(returned_mat) ≈ β * Cc  rtol=rtol
            end
        end
    end
end

end  # module

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
module CodeGen

using ExprTools
using MatchCore
using Comonicon.Types

export ASTCtx, ZSHCompletionCtx
export codegen, rm_lineinfo, prettify, pushmaybe!

include("utils.jl")
include("ast.jl")
include("completion.jl")

end


text:
module LinearAlgebraStudy
using Symbolics
using SymbolicUtils
using Plots
using LinearAlgebra

export drawVector, testVectorsIndependance, transformVectorBasis2D, transformVectorBasis4D

"""
    drawVector(f, V::Vector, origin::Vector = [0,0], xlim=(0,5), ylim=(0,5), title="Vector Plot", label = "Vector 1"))

    Draws a Vector V in the current plot at the origin.

    # Parameters
    f: Plot function
        plot or plot!
    V : Vector
        The vector to be drawn.
    origin : Origin Coordinate
        The origin of the coordinate system.

    # Returns
    plot :
        The plot of the vector.

    # Example
   ```julia
    using LinearAlgebraStudy
    V = Vector([1,2])
    origin = [0,0]
    drawVector(V, origin, (0,5), (0,5), "Vector Plot", "Vector 1")
    ```
"""
function drawVector(f, V::Vector, origin::Vector = [0,0], xlim=(0,5), ylim=(0,5), title="Vector Plot", label = "Vector 1")

    return f([origin[1], V[1] + origin[1]], [origin[2], V[2] + origin[2]],  xlims = xlim, ylims = ylim, line=:arrow, title = title, label = label)

end

"""
testVectorsIndependance(kwargs...)
    This function will take a variable number of vectors and test if they are independant.
    A Matrix "M" will be created with the vectors as columns.
    The rank of the Matrix will then be calculated.
    If the rank is equal to the number of vectors, then the vectors are independant.

    # Parameters
    kwargs : Vectors
        The vectors to be tested.

    # Returns
    result : Bool
        True if the vectors are independant, false otherwise.

"""
function testVectorsIndependance(kwargs...)
    M = hcat(kwargs...)
    rk = rank(M)
    if rk == size(M)[2]
        return true
    else
        return false
    end
end

"""
    transformVectorBasis2D(P::Vector, BX::Vector, BY::Vector)
        Function to transform a given Vector P point to a new Basis. Original Basis is assumed to be
        the standard cartesian basis.
        The transformed Vector is the equivalent of the original Vector P in the new Basis.

        Vectors can be described in terms of scalar multiplication of the basis vectors:
        Example:
            P=[2,1] is the original Vector with standard cartersian basis (X=[1,0],Y=[0,1]).
            This can be written as linear combination of the basis vectors:
                P=2*X+1*Y

            If we want to transform this Vector to a new basis (X'=[1,1],Y'=[0,2]), we can write the
            new vector on the new basis as a linear combination and then solve for the new scalar coefficients:

             ```

            P'=c₁*X'+ c₂*Y'  ⟹ c₁ * [1,1] + c₂ * [0,2] = [2,1]

            let A = [1 0   # Matrix A with new basis vectors
                     1 2]

            let b = [2 1] # Vector b with original coordinates of point P.

            let c = Vector of new coordinates on new Basis corresponding to original coordinates point P.

                c = A\b

            ```

    # Parameters
    b::Vector
        Vector to be transformed.
    BX::Vector
        X-component of the new basis.
    BY::Vector
        Y-component of the new basis.

    # Returns
    p::Vector
        Transformed vector.

    # Examples

"""
function transformVectorBasis2D(b::Vector, BX::Vector, BY::Vector)
    A = [BX BY] # Matrix A with new basis vectors as columns

    c = A\b # Solve for new coordinates on new Basis corresponding to original coordinates point P.

    return c
end

"""
    transformVectorBasis2D(P::Vector, BX::Vector, BY::Vector)
        Function to transform a given Vector P point to a new Basis. Original Basis is assumed to be
        the standard cartesian basis.
        The transformed Vector is the equivalent of the original Vector P in the new Basis.

        Vectors can be described in terms of scalar multiplication of the basis vectors:
        Example:
            P=[2,1] is the original Vector with standard cartersian basis (X=[1,0],Y=[0,1]).
            This can be written as linear combination of the basis vectors:
                P=2*X+1*Y

            If we want to transform this Vector to a new basis (X'=[1,1],Y'=[0,2]), we can write the
            new vector on the new basis as a linear combination and then solve for the new scalar coefficients:

             ```

            P'=c₁*X'+ c₂*Y'  ⟹ c₁ * [1,1] + c₂ * [0,2] = [2,1]

            let A = [1 0   # Matrix A with new basis vectors
                     1 2]

            let b = [2 1] # Vector b with original coordinates of point P.

            let c = Vector of new coordinates on new Basis corresponding to original coordinates point P.

                c = A\b

            ```

    # Parameters
    b::Vector
        Vector to be transformed.
    BX::Vector
        X-component of the new basis.
    BY::Vector
        Y-component of the new basis.

    # Returns
    p::Vector
        Transformed vector.

    # Examples

"""
function transformVectorBasis4D(b::Vector, BX::Vector, BY::Vector, BZ::Vector, BW::Vector)
    A = [BX BY BZ BW] # Matrix A with new basis vectors as columns

    c = A\b # Solve for new coordinates on new Basis corresponding to original coordinates point P.

    return c
end

end

max_stars_repo_name: gjunqueira-sys/LinearAlgebraStudy


text:
#=
Classify Clouds based on attributes and surrounding points
Format should be according to:

    classify!cloud::Cloud; kwargs)
        cloud'
    end

=#

## BASICS
" classification of point in Cloud"
function classify!(cloud::Cloud, i::Integer, class::UInt8)
    cloud[:classification][i] = class
    nothing
end

" classification of of Cloud based on pointfilter"
function classify!(cloud::Cloud, pointfilter, class::UInt8)
    for i in 1:length(cloud)
        (pointfilter(cloud, i) || continue)
        cloud[:classification][i] = class
    end
    nothing
end

"set classification of a point to zero"
function reset_class!(cloud::Cloud, i::Integer)
    classify!(cloud, i, UInt8(0))
    nothing
end

"set classification of all points in Cloud to zero"
function reset_class!(cloud::Cloud;
    pointfilter = nothing)

    for i in 1:length(cloud)
      (pointfilter != nothing) && (pointfilter(cloud, i) || continue)
      classify!(cloud, i, UInt8(0))
    end
    nothing
end

## Raster based
"function to classify minima in pointcloud based per gridcell"
function classify_min!(cloud::Cloud, cellsize::Real, class::Int;
   pointfilter = nothing, # predicate function to filter individual points
   min_dens = 0) # minimum point density in cell to consider for educed cloud

    r = define_raster(cloud, cellsize; pointfilter = pointfilter)
    classify_min!(cloud, r, class; pointfilter = pointfilter, min_dens = min_dens)
end

"function to classify minima in pointcloud based per gridcell"
function classify_min!(cloud::Cloud, r::Raster, class::Int;
    pointfilter = nothing, # predicate function to filter individual points
    min_dens = 0) # minimum point density in cell to consider for educed cloud

    subset = reduce_index(cloud, r; reduceri = reducer_minz_index, pointfilter = pointfilter, min_dens = min_dens)

    # return raster with statistics, density is saved to last layer
    for i in subset
        classify!(cloud, i, UInt8(class))
    end

    nothing
end

"function to classify maxima in pointcloud based per gridcell"
function classify_max!(cloud::Cloud, cellsize::Real, class::Int;
    pointfilter = nothing, # predicate function to filter individual points
    min_dens = 0) # minimum point density in cell to consider for educed cloud

    r = define_raster(cloud, cellsize; pointfilter = pointfilter)
    classify_max!(cloud, r, class; pointfilter = pointfilter, min_dens = min_dens)
end

"function to classify maxima in pointcloud based per gridcell"
function classify_max!(cloud::Cloud, r::Raster, class::Int;
    pointfilter = nothing, # predicate function to filter individual points
    min_dens = 0) # minimum point density in cell to consider for educed cloud

    subset = reduce_index(cloud, r; reduceri = reducer_maxz_index, pointfilter = pointfilter, min_dens = min_dens)

    # return raster with statistics, density is saved to last layer
    for i in subset
        classify!(cloud, i, UInt8(class))
    end

    nothing
end

"classifies points with a cell if z <= to the z-value in surface + tolerance"
function classify_below_surface!(
    cloud::Cloud, r::Raster, surf::AbstractArray{T, 2}, class::Integer;
    tolerance = 0.0,
    pointfilter = nothing,
    nodata = -9999.0) where T <: Real

    ((size(r)[1] == size(surf)[1]) && (size(r)[2] == size(surf)[2])) || error("raster and surface size mismatch")
    for icell in 1:length(r)
        z_max = surf[icell] + tolerance
        z_max == nodata && continue
        for i in r[icell] # indices Cloud
            positions(cloud)[i][3] <= z_max || continue
            (pointfilter != nothing) && (pointfilter(cloud, i) || continue)
            classify!(cloud, i, UInt8(class))
        end
    end

    nothing
end

function classify_below_surface!(cloud::Cloud, surf::Real, class::Integer;
    tolerance = 0.0,
    pointfilter = nothing)

   z_max = surf + tolerance
   for i in 1:length(position(cloud)) # indices Cloud
      positions(cloud)[i][3] <= z_max || continue
      (pointfilter != nothing) && (pointfilter(cloud, i) || continue)
      classify!(cloud, i, UInt8(class))
   end

   nothing
end

"classifies points with a cell if z > to the z-value in surface + buffer"
function classify_above_surface!(cloud::Cloud, r::Raster,
    surf::AbstractArray{T, 2}, class::Integer;
    buffer = 0.0,
    pointfilter = nothing,
    nodata = -9999.0) where T <: Real

    ((size(r)[1] == size(surf)[1]) && (size(r)[2] == size(surf)[2])) || error("raster and surface size mismatch")
    for icell in 1:length(r)
        z_min = surf[icell] + buffer
        z_min == nodata && continue
        for i in r[icell] # indices Cloud
          positions(cloud)[i][3] > z_min || continue
          (pointfilter != nothing) && (pointfilter(cloud, i) || continue)
          classify!(cloud, i, UInt8(class))
        end
    end

    nothing
end

"classify points outside of raster mask where true indicates data not to classify"
function classify_mask!(cloud::Cloud, r::Raster, mask::Union{BitArray, Array{Bool}}, class::Integer)
    size(r) == size(mask) || error("raster and mask size mismatch")
    for icell in 1:length(r)
        mask[icell] && continue  # skip points inside mask
        for i in r[icell] # pointcloud indices
          classify!(cloud, i, UInt8(class))
        end
    end
    nothing
end

function classify_above_surface!(cloud::Cloud, surf::Real, class::Integer;
    buffer = 0.0,
    pointfilter = nothing)

    z_min = surf + buffer
    for i in 1:length(position(cloud)) # indices Cloud
      positions(cloud)[i][3] > z_min || continue
      (pointfilter != nothing) && (pointfilter(cloud, i) || continue)
      classify!(cloud, i, UInt8(class))
    end

    nothing
end

## Slope based
"""Adapted from Vosselman (2000)
slope based classification of pointdata from point cloud"""
function vosselman(cloud::Cloud, i::Integer, neighbors::Vector{U};
    max_slope=0.2,
    tolerance=0.0) where U <: Integer

    if length(neighbors) > 1
        # get positions point i
        xi, yi, zi = positions(cloud)[i]

        # get z values
        z_vec = getz(cloud, neighbors) # z only
        idx_order = sortperm(z_vec) # loop through points from min to max elevation

        # initiate
        isclass = false
        for j in neighbors[idx_order]
            cloud[:classification][j] == 7 && continue # skip comparison with low outlier
            # 1) get horizontal and vertical distance
            xj, yj, zj = positions(cloud)[j]
            di = hypot(xi - xj, yi - yj)
            dz = zj - zi + tolerance
            # 2) calculate max z difference based on distance
            thresh = -1 * di * max_slope
            # 3) compare with actual z diff; if z diff too large -> no ground point
            dz < thresh && break
            # 4) classify if all points below point i have been checked without break
            if dz > 0
                isclass = true
                break
            end
        end
      else
          isclass = true # no neighbors
      end

    isclass
end

## applied
"classify outliers"
function classify_outliers!(cloud::Cloud;
    cellsize=100.0, # cellsize size for regural grid
    dz=1.0, # min distance for outlier classification
    max_outliers=15, # max number of outlier per cell
    max_height=Inf) # maximum heigh above lowest non outlier in cell

    low_noise = UInt8(7) # from ASPRS Standard LIDAR Point Classes
    high_noise = UInt8(18) # from ASPRS Standard LIDAR Point Classes (LAS 1.4)
    do_high_noise = !isinf(max_height)
    # all points below -5 m are always set to low noise
    const cutoff_low = -5.0
    too_low(cloud::Cloud, i::Integer) = getz(cloud, i) < cutoff_low
    classify!(cloud, too_low, low_noise)
    # all points above 100 m are always set to high noise
    const cutoff_high = 100.0
    too_high(cloud::Cloud, i::Integer) = getz(cloud, i) > cutoff_high
    classify!(cloud, too_high, high_noise)
    # points outside the cutoff should not be included in the outlier analysis
    not_cut_off(cloud::Cloud, i::Integer) = cutoff_low <= getz(cloud, i) <= cutoff_high
    r = define_raster(cloud, cellsize, pointfilter=not_cut_off)
    # create dz Vector with length of pointcloud to store as new attribute later
    dz_outliers = zeros(length(positions(cloud)))
    # loop through cells
    @showprogress 5 "Classifying outliers.." for icell in 1:length(r)
        # find points in range
        idxs = r[icell] # get points in rastercell
        np = length(idxs)
        np < (max_outliers*10) && continue # max 10% outliers. make sure <max_outliers> is in balance with <cellsize>

        # get z values of points in cell
        z_vec = getz(cloud, idxs)
        if np >= 1e5 # RadixSort only faster with many points
            idx_order = sortperm(z_vec, alg=RadixSort)
        else
            idx_order = sortperm(z_vec)
        end

        # get difference in z for lowest max_outliers+1 points
        dz_vec = z_vec[idx_order[2:max_outliers+1]] - z_vec[idx_order[1:max_outliers]]
        if sum(dz_vec .>= dz) == 0 && !do_high_noise
            continue # no low outliers
        end
        ilast = findlast(dz_vec .>= dz) # find the index of the last jump
        zmin = z_vec[idx_order[ilast+1]] # lowest non outlier

        if do_high_noise
            # classify points higher than max_height above lowest non outlier as high noise
            zdif = z_vec - zmin
            toohi = zdif .> max_height
            itoohi = idxs[toohi]
            dz_outliers[itoohi] = zdif[toohi]
            for i in itoohi
                classify!(cloud, i, high_noise)
            end
        end

        for i in idxs[idx_order[1:ilast]]
            dz_outliers[i] = zmin - getz(cloud, i) # difference in elevation with lowes non outlier
            classify!(cloud, i, low_noise)
        end
    end
    # update Cloud with dz of outliers
    cloud[:dz] = dz_outliers
    nothing
end

"Apply Vosselman (2000) to classify ground points"
function classify_ground!(cloud::Cloud, spatialindex::NearestNeighbors.NNTree;
    radius=3.0, # radius to check slope condition [m]
    max_slope=0.3, # max allowd slope between ground points [m/m]
    tolerance=0.0, # vertical tolerance [m]
    pointfilter=nothing, # only classify points with class <class_only>
    class=2)

    # check if the spatialindex is 2D Euclidean
    @assert (spatialindex.metric == Euclidean()) && (length(spatialindex.data[1]) == 2)
    n = length(positions(cloud))
    # loop through points
    @showprogress 5 "Classifying low surface points (Vosselman).." for i in 1:n
        # get xyz & class of point i
        (pointfilter != nothing) && (pointfilter(cloud, i) || continue)

        xi, yi, zi = positions(cloud)[i]
        # find number of points in range
        neighbors = inrange(spatialindex, [xi, yi], radius)

        # classify
        isground = vosselman(cloud, i, neighbors; max_slope=max_slope, tolerance=tolerance)

        # classes according to asprs las format
        if isground
            classify!(cloud, i, UInt8(class)) # ground
        end
    end

    nothing
end

"""Apply Vosselman (2000) to classify water points
only apply to points below min_z_tolerance from minima in grid with cellsize=radius"""
function classify_water!(cloud::Cloud, spatialindex::NearestNeighbors.NNTree;
    radius=50.0, # radius to check slope condition [m]
    max_slope=0.01, # max allowd slope between ground points [m/m]
    tolerance=0.05, # vertical tolerance vosselman [m] (due to vairance in water level measurements)
    min_z_tolerance = 0.3, # only apply to points below min_z_tolerance from minima in grid with cellsize=radius
    pointfilter = nothing, # only classify points with class <class_only>
    class = 9)

    r_large = define_raster(cloud, radius; pointfilter = pointfilter)
    z_max = rasterize(cloud, r_large, reducer = reducer_minz) + min_z_tolerance

    # check if the spatialindex is 2D Euclidean
    @assert (spatialindex.metric == Euclidean()) && (length(spatialindex.data[1]) == 2)

    # loop through points
    @showprogress 5 "Classifying water points (Vosselman).." for icell in 1:length(r_large)
        max_z = z_max[icell]
        idx0 = r_large[icell]
        for i in idx0
            # get xyz & class of point i
            (pointfilter != nothing) && (pointfilter(cloud, i) || continue)

            xi, yi, zi = positions(cloud)[i]
            zi <= max_z || continue

            # find number of points in range
            neighbors = inrange(spatialindex, [xi, yi], radius)

            # classify
            iswater = vosselman(cloud, Int(i), neighbors; max_slope=max_slope, tolerance=tolerance)

            # classes according to asprs las format
            if iswater
                classify!(cloud, i, UInt8(class))  # ground
            end
        end
    end

    nothing
end


text:
#!/usr/bin/env julia

include("../src/MsgPackRpcClient.jl")
include("dataset.jl")
include("const.jl")
using MsgPackRpcClient
using Base.Test

function test_data(data, number = 1)
  for d in data
    println("$(number): $d ")
    try
      result = MsgPackRpcClient.call(session, d["method"], d["arg"]; sync = true)
      #println("d[\"arg\"] = ", d["arg"], ", result = $result", ", first(result]) = ", first(result))
      if haskey(d, "expect")
        @test first(result) == d["expect"]
      else
        @test first(result) == d["arg"]
      end
      println("  => OK")
    catch
      println("Error on d[\"arg\"] = ", d["arg"])
    end
    number += 1
  end
end

session = MsgPackRpcClientSession.create()

test_data(get_dataset())

session.destroy()

max_stars_repo_name: keithseahus/msgpack-rpc-julia

max_issues_repo_name: keithseahus/msgpack-rpc-julia

max_forks_repo_name: keithseahus/msgpack-rpc-julia


text:
abstract type AbstractOptimizer end

@params mutable struct MMAOptimizer{T} <: AbstractOptimizer
    model::Model{T}
    mma_alg
    suboptimizer
    obj
    constr
    workspace
    convstate::ConvergenceState
    options::MMA.Options
end
GPUUtils.whichdevice(o::MMAOptimizer) = o.model
multol!(o::MMAOptimizer, m::Real) = o.options.tol *= m

function Functions.maxedfevals(o::MMAOptimizer)
    maxedfevals(o.obj) || maxedfevals(o.constr)
end
@inline function Functions.maxedfevals(c::Tuple{Vararg{Constraint}})
    maxedfevals(c[1]) || maxedfevals(Base.tail(c))
end
function Functions.maxedfevals(c::Vector{<:Constraint})
    all(c -> maxedfevals(c), c)
end

function MMAOptimizer(args...; kwargs...)
    return MMAOptimizer{CPU}(args...; kwargs...)
end
function MMAOptimizer{T}(args...; kwargs...) where T
    return MMAOptimizer(T(), args...; kwargs...)
end
function MMAOptimizer(  device::Tdev,
                        obj::Objective{<:AbstractFunction{T}},
                        constr,
                        opt = MMA.MMA87(),
                        subopt = Optim.ConjugateGradient();
                        options = MMA.Options(),
                        convcriteria = KKTCriteria()
                    ) where {T, Tdev}

    solver = getsolver(obj)
    nvars = length(solver.vars)
    xmin = solver.xmin

    model = Model{Tdev}(nvars, obj)

    box!(model, zero(T), one(T))
    ineq_constraint!(model, constr)

    if Tdev <: CPU && whichdevice(obj) isa GPU
        x0 = Array(solver.vars)
    else
        x0 = solver.vars
    end

    workspace = MMA.Workspace(model, x0, opt, subopt; options = options, convcriteria = convcriteria)

    return MMAOptimizer(model, opt, subopt, obj, constr, workspace, ConvergenceState(T), options)
end

Utilities.getpenalty(o::MMAOptimizer) = getpenalty(o.obj)
function Utilities.setpenalty!(o::MMAOptimizer, p)
    setpenalty!(o.obj, p)
    setpenalty!.(o.constr, p)
end

function (o::MMAOptimizer)(x0::AbstractVector)
    @unpack workspace, options = o
    @unpack model = workspace
    @unpack x, g, ∇g, ∇f_x = workspace.primal_data

    reset!(workspace, x0)
    setoptions!(workspace, options)
    mma_results = MMA.optimize!(workspace)
    o.convstate, mma_results.convstate
    return mma_results
end
function setoptions!(workspace, options)
    @unpack options = workspace
    @unpack store_trace, show_trace, extended_trace, dual_caps = options
    @unpack maxiter, tol, subopt_options = options
    @pack! options = dual_caps, subopt_options, tol
    @pack! options = show_trace, extended_trace, store_trace

    return workspace
end

function reset!(w::Workspace, args...)
    @unpack primal_data, model = w
    # primal_data.r0 = 0
    # Assess multiple types of convergence
    outer_iter, iter, f_calls, g_calls = 0, 0, 0, 0
    @pack! w = f_calls, g_calls
    MMA.update_values!(w, args...)
    w.convstate = MMA.assess_convergence(w)
    @pack! w = outer_iter, iter
    return w
end

# For adaptive SIMP
function (o::MMAOptimizer)(workspace::MMA.Workspace)
    mma_results = MMA.optimize!(workspace)
    return mma_results
end

@params struct MMAOptionsGen
    maxiter
    outer_maxiter
    tol
    s_init
    s_incr
    s_decr
    dual_caps
    store_trace
    show_trace
    extended_trace
    subopt_options
end
function (g::MMAOptionsGen)(i)
    MMA.Options(
        g.maxiter(i),
        g.outer_maxiter(i),
        g.tol(i),
        g.s_init(i),
        g.s_incr(i),
        g.s_decr(i),
        g.dual_caps(i),
        g.store_trace(i),
        g.show_trace(i),
        g.extended_trace(i),
        g.subopt_options(i)
    )
end

function (g::MMAOptionsGen)(options, i)
    MMA.Options(
        optionalcall(g, :maxiter, options, i),
        optionalcall(g, :outer_maxiter, options, i),
        optionalcall(g, :tol, options, i),
        optionalcall(g, :s_init, options, i),
        optionalcall(g, :s_incr, options, i),
        optionalcall(g, :s_decr, options, i),
        optionalcall(g, :dual_caps, options, i),
        optionalcall(g, :store_trace, options, i),
        optionalcall(g, :show_trace, options, i),
        optionalcall(g, :extended_trace, options, i),
        optionalcall(g, :subopt_options, options, i)
    )
end
function optionalcall(g, s, options, i)
    if getproperty(g, s) isa Nothing
        return getproperty(options, s)
    else
        return getproperty(g, s)(i)
    end
end


text:
#Example demonstrating the use of overlap to solve a long horizon control problem
using Plasmo, Ipopt
using KaHyPar
using SchwarzOpt

T = 100              #number of time points
d = sin.(1:T)        #a disturbance vector
imbalance = 0.1      #partition imbalance
distance = 5         #expand distance
n_parts = 6          #number of partitions

#Create the model (an optigraph)
graph = OptiGraph()
@optinode(graph,state[1:T])
@optinode(graph,control[1:T-1])

for (i,node) in enumerate(state)
    @variable(node,x)
    @constraint(node, x >= 0)
    @objective(node,Min,0.001*x^2) #- 2*x*d[i])
end
for node in control
    @variable(node,u)
    @constraint(node, u >= -1000)
    @objective(node,Min,u^2)
end
n1 = state[1]
@constraint(n1,n1[:x] == 0)

for i = 1:T-1
    @linkconstraint(graph, state[i][:x] + control[i][:u] + d[i] == state[i+1][:x],attach = state[i+1])
end

#Partition the optigraph using recrusive bisection over a hypergraph
hypergraph,hyper_map = hyper_graph(graph) #create hypergraph object based on graph

partition_vector = KaHyPar.partition(hypergraph,n_parts,
configuration = (@__DIR__)*"/cut_rKaHyPar_sea20.ini",
imbalance = imbalance)

partition = Partition(hypergraph,partition_vector,hyper_map)
apply_partition!(graph,partition)

#Inspect the graph structure. It should be RECURSIVE_GRAPH
println(Plasmo.graph_structure(graph))

#calculate subproblems using expansion distance
subgraphs = getsubgraphs(graph)
expanded_subgraphs = Plasmo.expand.(graph,subgraphs,distance)
sub_optimizer = optimizer_with_attributes(Ipopt.Optimizer,"print_level" => 0)

#optimize using schwarz overlapping decomposition
SchwarzOpt.optimize!(graph;
subgraphs = expanded_subgraphs,
sub_optimizer = sub_optimizer,
max_iterations = 50)


text:
"""
  BagChain(layers...)
	BagChain multiple layers / functions together, so that they are called in sequence
	on a given input supported by bags.

"""
struct BagChain{T<:Tuple}
  layers::T
  BagChain(xs...) = new{typeof(xs)}(xs)
end

Flux.@treelike BagChain
# Flux.@functor BagChain

Flux.@forward BagChain.layers Base.getindex, Base.first, Base.last, Base.lastindex
Flux.@forward BagChain.layers Base.iterate

applychain(::Tuple{}, x, bags) = x
applychain(fs::Tuple, x, bags) = applychain(Base.tail(fs), first(fs)(x, bags), bags)

(c::BagChain)(x, bags) = applychain(c.layers, x, bags)
(c::BagChain)(x::ArrayNode, bags) = ArrayNode(applychain(c.layers, x.data, bags))
(c::BagChain)(x::BagNode) = applychain(c.layers, x.data, x.bags)

Base.getindex(c::BagChain, i::AbstractArray) = BagChain(c.layers[i]...)

Base.show(io::IO, m::BagChain) = modelprint(io, m)

function modelprint(io::IO, m::BagChain; pad=[])
    c = COLORS[(length(pad)%length(COLORS))+1]
    paddedprint(io, "BagChain (\n", color=c)
    for j in 1:length(m.layers) - 1
        paddedprint(io, "  ├── ", color=c, pad=pad)
	    modelprint(io, m.layers[j], pad=[pad; (c, "  │   ")])
    end
    paddedprint(io, "  └── ", color=c, pad=pad)
    modelprint(io, m.layers[end], pad=[pad; (c, "      ")])
    paddedprint(io, ")\n", color=c, pad=pad)
end

activations(c::BagChain, x, bags) = accumulate((x, m) -> m(x, bags), c.layers, init = x)


text:
"""
Convenient structure for postprocessing. It is less performant than having `SArray`s but is easier to manipulate.
"""
struct SpacecraftTrajectory
    position::Array{Float64, 2}
    attitude::Array{Float64, 2}
    velocity::Array{Float64, 2}
    angular_velocity::Array{Float64, 2}
    thrust::Array{Float64, 2}
    angular_controls::Array{Float64, 2}
    time
    SpacecraftTrajectory(states, controls, time) = new(
        states[1:3, :],
        states[4:7, :],
        states[7:9, :],
        states[10:12, :],
        controls[1:1, :],
        controls[2:4, :],
        time)
end

function Base.write(model::Spacecraft, traj::SpacecraftTrajectory, filename)
    dict = Dict(
        "$name" => getproperty(traj, name) for name ∈ fieldnames(SpacecraftTrajectory)
    )
    open(filename, "w") do io
        write(io, json(dict))
    end
end

"""
Print some information about the solved problem.
"""
function postprocess(solver, x, u)
    println("Cost: ", cost(solver))
    println("Constraint violation: ", max_violation(solver))
    println("Iterations: ", iterations(solver))
    println("Final state: ", x[:, end])
    println("Maximum control: ", maximum(u))
end


text:
using Random
using Plots
#=
Julia Func B:
 This is a sequence of numbers I found on the OEIS and expanded to large porportions.
 I was shocked to find the minimum would always be -199.

 Perhaps I'll come back and redo this one later. i got a little frustrated, so this looks pretty bad
 with how the variables are laid out.
=#

numMax = 200

nums = ( function ()
  global lpMax
  maxiNum::UInt64 = (2^(63)-1)
  arrPos = 1
  numArr = zeros(Int64, numMax)
  itrcnt = 3
  chck = true
  i = 1
  num = 0
  str = 0
  while i <= numMax

    if (chck == true)
     str = string(-(i))
     num = parse(Int64, str[1:2])
     chck = false
   else
     num += parse(Int8, str[itrcnt])
     itrcnt += 1
   end

   if (itrcnt > length(str))
   num ^= 2
   num += -(i)
     if ((string(num))[1] != '-' && num != 0 ) numArr[arrPos] = -(i)
     println(arrPos, " ", -(i))
     arrPos += 1
     end
     chck = true
     itrcnt = 3
     i += 1
   end
   @assert (num >= maxiNum)  ? false : true (string( "Input too high; Num: ", num))

  end
return numArr
end)()

function Plot(maxNums, nums)
  println("Would you like to see the plot?")
  if (readline() == "y")
    println("Enter Character in Console to Exit")
    gui(plot(1:maxNums,nums))
  else
    exit()
  end
  readline()
end

Plot(numMax, nums)
exit()


text:
# This file was generated, do not modify it. # hide
nullspace(W') ≈ nullspace(Q')

max_stars_repo_name: cscherrer/InformativePrior

max_issues_repo_name: cscherrer/InformativePrior

max_forks_repo_name: cscherrer/InformativePrior


text:
# Only for simple rules
"""
    generate_caos_rules(tree_file_path::String, output_directory::String)

Takes a Nexus file and generates all the CAOS rules for the tree.

# Arguments
- `tree_file_path::String`: path to the Nexus file.
- `output_directory::String`: path to the output directory.
"""
function generate_caos_rules(tree_file_path::String, output_directory::String; protein::Bool=false)

    # Create directory
    try
        mkdir("$output_directory")
    catch
    end

    # Parse the tree and prepare data for CA's
    nodes, taxa_labels, character_labels, _ = parse_tree(tree_file_path)

    # Initialize variables for constructing tree
    sPu, sPr = get_sPu_and_sPr(nodes, 1, taxa_labels, character_labels; protein=protein)
    cPu = Array{Dict{String,Any}}(undef, 0)
    cPr = Array{Dict{String,Any}}(undef, 0)
    #cPu,cPr = get_cPu_and_cPr(nodes,1,taxa_labels,character_labels,sPu,sPr)

    # Get CA's for tree
    tree = Node([])
    add_nodes!(tree,sPu,sPr,cPu,cPr,taxa_labels,character_labels,nodes,1,complex=false, protein=protein)

    # Save tree to json
    tree_data = JSON.json(tree)
    open("$output_directory/caos_rules.json", "w") do f
        write(f, tree_data)
    end

    # Change blanks to N's for the database
    character_labels_no_gaps = remove_blanks(character_labels)

    # Create fasta file from dictionary of character labels
    writefasta("$output_directory/char_labels.fasta", character_labels_no_gaps)

    # Create the blast database
    if protein
        run(`makeblastdb -in $output_directory/char_labels.fasta -dbtype prot`)
    else
    run(`makeblastdb -in $output_directory/char_labels.fasta -dbtype nucl`)
    end

    # Save character and taxa labels
    character_data = JSON.json(character_labels)
    open("$output_directory/character_labels.json", "w") do f
        write(f, character_data)
    end
    taxa_data = JSON.json(taxa_labels)
    open("$output_directory/taxa_labels.json", "w") do f
        write(f, taxa_data)
    end

    return tree, character_labels, taxa_labels

end

"""
    classify_new_sequence(tree::Node, character_labels::Dict{String,String}, taxa_labels::Dict{String,String},
    sequence_file_path::String, output_directory::String ; all_CA_weights::Dict{Int64,
    Dict{String,Int64}}=Dict(1=>Dict("sPu"=>1,"sPr"=>1,"cPu"=>1,"cPr"=>1)), occurrence_weighting::Bool=false,
    tiebreaker::Vector{Dict{String,Int64}}=[Dict{String,Int64}()], combo_classification::Bool=false)

Takes a tree (Node) and a sequence, and classifies the new sequence using the CAOS tree.

# Arguments
- `tree::Node`: the tree represented as a Node.
- `character_labels::Dict{String,String}`: a mapping of the character labels to the corresponding sequences.
- `taxa_labels::Dict{String,String}`: a mapping of the taxa labels to the character labels.
- `sequence_file_path::String`: a file path to the sequence to classify.
- `output_directory::String`: path to the output directory.
- `all_CA_weights::Dict{Int64,Dict{String,Int64}}=Dict(1=>Dict("sPu"=>1,"sPr"=>1,"cPu"=>1,"cPr"=>1))`: CA weights to be used.
- `occurrence_weighting::Bool=false`: whether to use occurence weighting in classification.
- `tiebreaker::Vector{Dict{String,Int64}}=[Dict{String,Int64}()]`: tiebreaker to be used in classification.
- `combo_classification::Bool=false`: whether to use a combo of Blast and CAOS for classification.
"""
function classify_new_sequence(tree::Node, character_labels::Dict{String,String}, taxa_labels::Dict{String,String},
    sequence_file_path::String, output_directory::String ; all_CA_weights::Dict{Int64,Dict{String,Int64}}=Dict(1=>Dict("sPu"=>1,"sPr"=>1,"cPu"=>1,"cPr"=>1)),
    occurrence_weighting::Bool=false, tiebreaker::Vector{Dict{String,Int64}}=[Dict{String,Int64}()], combo_classification::Bool=false, protein::Bool=false)

    # Create directory
    try
        mkdir("$output_directory")
    catch
    end

    character_labels_no_gaps = remove_blanks(character_labels, change_to_N=false)

    # Get the new sequence after imputing blanks
    new_seq = add_blanks("$sequence_file_path", "$output_directory/char_labels.fasta", character_labels, character_labels_no_gaps, protein=protein)
    classification = classify_sequence(new_seq, tree, all_CA_weights[1], all_CA_weights, occurrence_weighting, 1,
                                       tiebreaker, combo_classification=combo_classification, protein=protein)

    return classification

end

# Function to read a tree from file and convert it to a node object
"""
    load_tree(directory::String)

Loads a CAOS tree from file.

# Arguments
- `directory::String`: path to directory where tree exists.
"""
function load_tree(directory::String)

    # Load tree from json
    open("$directory/caos_rules.json", "r") do f
        global loaded_tree
        loaded_tree=JSON.parse(f)  # parse and transform data
    end

    open("$directory/character_labels.json", "r") do f
        global character_labels
        character_labels=convert(Dict{String,String}, JSON.parse(f))  # parse and transform data
    end

    open("$directory/taxa_labels.json", "r") do f
        global taxa_labels
        taxa_labels=convert(Dict{String,String}, JSON.parse(f))  # parse and transform data
    end

    tree_load = Node([])

    # Get tree back into proper format
    tree = convert_to_struct(loaded_tree, tree_load)

    return tree, character_labels, taxa_labels

end


text:
module FastAI

using Reexport

@reexport using DLPipelines
@reexport using FluxTraining
@reexport using DataLoaders
@reexport using Flux

using Animations
using AbstractPlotting
using Colors
using DataAugmentation
using DataAugmentation: getbounds, makebounds
using DLPipelines: methoddataset, methodmodel, methodlossfn, methoddataloaders
using LearnBase: getobs, nobs
using FilePathsBase
using FixedPointNumbers
using Flux
using FluxTraining: Learner, handle
using FluxTraining.Events
using MLDataPattern
using Parameters
using StaticArrays
using WordTokenizers: spaces, TokenBuffer, character, isdone, flush!

include("tasks.jl")
include("plotting.jl")
include("learner.jl")

# method implementations and helpers
include("./steps/utils.jl")
include("./steps/step.jl")
include("./steps/spatial.jl")
include("./steps/imagepreprocessing.jl")
include("./methods/imageclassification.jl")
include("./methods/imagesegmentation.jl")

# submodules
include("datasets/Datasets.jl")
using .Datasets

include("models/Models.jl")
using .Models

# training
include("training/utils.jl")
include("training/onecycle.jl")
include("training/finetune.jl")
include("training/lrfind.jl")

export methodlossfn

export
    # submodules
    Datasets,
    Models,
    datasetpath,
    loadtaskdata,
    mapobs,
    groupobs,
    filterobs,
    shuffleobs,

    # method API
    methodmodel,
    methoddataset,
    methoddataloaders,
    methodlossfn,
    getobs,
    nobs,

    # pipeline steps
    ProjectiveTransforms, ImagePreprocessing,

    # tasks
    ImageClassificationTask,
    ImageSegmentationTask,

    # methods
    ImageClassification,
    ImageSegmentation,

    # training
    methodlearner,
    Learner,
    fit!,
    fitonecycle!,
    finetune!,

    gpu

end  # module


text:
using Yao
using YaoExtensions: variational_circuit, Sequence, faithful_grad, numdiff
using QuAlgorithmZoo: Adam, update!
import Yao: tracedist

"""
Quantum GAN.

Reference:
    Benedetti, M., Grant, E., Wossnig, L., & Severini, S. (2018). Adversarial quantum circuit learning for pure state approximation, 1–14.
"""
struct QuGAN{N}
    target::ArrayReg
    generator::AbstractBlock{N}
    discriminator::AbstractBlock
    reg0::ArrayReg
    witness_op::AbstractBlock
    circuit::AbstractBlock

    function QuGAN(target::ArrayReg, gen::AbstractBlock, dis::AbstractBlock)
        N = nqubits(target)
        c = Sequence([gen, addbits!(1), dis])
        witness_op = put(N+1, (N+1)=>ConstGate.P0)
        new{N}(target, gen, dis, zero_state(N), witness_op, c)
    end
end

# INTERFACES
circuit(qg::QuGAN) = qg.circuit
loss(qg::QuGAN) = p0t(qg) - p0g(qg)

function gradient(qg::QuGAN)
    grad_gen = faithful_grad(qg.witness_op, qg.reg0 => qg.circuit)
    grad_tar = faithful_grad(qg.witness_op, qg.target => qg.circuit[2:end])
    ngen = nparameters(qg.generator)
    [-grad_gen[1:ngen]; grad_tar - grad_gen[ngen+1:end]]
end

"""probability to get evidense qubit 0 on generation set."""
p0g(qg::QuGAN) = expect(qg.witness_op, qg.reg0 => qg.circuit) |> real
"""probability to get evidense qubit 0 on target set."""
p0t(qg::QuGAN) = expect(qg.witness_op, qg.target => qg.circuit[2:end]) |> real
"""generated wave function"""
outputψ(qg::QuGAN) = copy(qg.reg0) |> qg.generator

"""tracedistance between target and generated wave function"""
tracedist(qg::QuGAN) = tracedist(qg.target, outputψ(qg))[]

using Test, Random
Random.seed!(2)

nbit = 3
depth_gen = 4
depth_dis = 4

# define a QuGAN
target = rand_state(nbit)
generator = dispatch!(variational_circuit(nbit, depth_gen), :random)
discriminator = dispatch!(variational_circuit(nbit+1, depth_dis), :random)
qg = QuGAN(target, generator, discriminator)

# check the gradient
grad = gradient(qg)
numgrad = numdiff(c->loss(qg), qg.circuit)
@test isapprox(grad, numgrad, atol=1e-4)

# learning rates for the generator and discriminator
g_lr = 0.2
d_lr = 0.5
for i=1:300
    ng = nparameters(qg.generator)
    grad = gradient(qg)
    dispatch!(-, qg.generator, grad[1:ng]*g_lr)
    dispatch!(-, qg.discriminator, -grad[ng+1:end]*d_lr)
    println("Step $i, trace distance = $(tracedist(qg))")
end

@test qg |> loss < 0.1


text:
module Causality
using Docile
@docstrings

# using commands
using Graph, DataFrames, CSVReaders

# abstract declarations
# export commands

# include calls
include("pc_algo.jl")
include("utils.jl")
end


text:
#
# Scalar function lifted to a query combinator.
#

function LiftQuery(fn, argtypes::Tuple{Vararg{Type}}, restype::Type, qs::Vector{Query})
    @assert method_exists(fn, argtypes) "$fn($(join(argtypes, ", "))) is not defined"
    @assert length(argtypes) == length(qs)
    qs′ = Query[]
    for (argtype, q) in zip(argtypes, qs)
        if argtype <: Nullable
            T = eltype(argtype)
            if fits(output(q), Output(T) |> setoptional())
                q = RecordQuery(q)
                q = q >> Query(EncodeNullableSig(), Input(domain(output(q))), Output(Nullable{T}))
            end
        elseif argtype <: Vector
            T = eltype(argtype)
            if fits(output(q), Output(T) |> setoptional() |> setplural())
                q = RecordQuery(q)
                q = q >> Query(EncodeVectorSig(), Input(domain(output(q))), Output(Vector{T}))
            end
        end
        @assert fits(domain(output(q)), Domain(argtype)) "$q does not fit $argtype"
        push!(qs′, q)
    end
    argq = RecordQuery(qs′)
    sig = LiftSig(fn, argtypes, restype)
    ity = Input(domain(output(argq)))
    oty = Output(restype, obound([mode(output(q)) for q in qs′]))
    return argq >> Query(sig, ity, oty)
end

LiftQuery{T}(fn, argtypes::Tuple{Vararg{Type}}, restype::Type{Nullable{T}}, qs::Vector{Query}) =
    let q = invoke(LiftQuery, (Any, Tuple{Vararg{Type}}, Type, Vector{Query}), fn, argtypes, restype, qs)
        q >> Query(DecodeNullableSig(), Input(Nullable{T}), Output(T) |> setoptional())
    end

LiftQuery{T}(fn, argtypes::Tuple{Vararg{Type}}, restype::Type{Vector{T}}, qs::Vector{Query}) =
    let q = invoke(LiftQuery, (Any, Tuple{Vararg{Type}}, Type, Vector{Query}), fn, argtypes, restype, qs)
        q >> Query(DecodeVectorSig(), Input(Vector{T}), Output(T) |> setoptional() |> setplural())
    end

LiftQuery(fn, argtypes::Tuple{Vararg{Type}}, restype::Type, qs::Query...) =
    LiftQuery(fn, argtypes, restype, collect(qs))

LiftQuery(fn, argtypes::Tuple{Vararg{Type}}, qs::Vector{Query}) =
    let restype = Union{Base.return_types(fn, argtypes)...}
        LiftQuery(fn, argtypes, restype, qs)
    end

LiftQuery(fn, argtypes::Tuple{Vararg{Type}}, qs::Query...) =
    LiftQuery(fn, argtypes, collect(qs))

LiftQuery(fn, qs::Vector{Query}) =
    let argtypes = ((datatype(domain(output(q))) for q in qs)...)
        LiftQuery(fn, argtypes, qs)
    end

LiftQuery(fn, qs::Query...) =
    LiftQuery(fn, collect(qs))

immutable LiftSig <: AbstractPrimitive
    fn::Function
    argtypes::Tuple{Vararg{Type}}
    restype::Type
end

describe(io::IO, sig::LiftSig) = print(io, sig.fn, "′")

ev(sig::LiftSig, dv::DataVector) =
    all(isplain, columns(dv)) ?
        plain_lift_impl(sig.fn, sig.restype, length(dv), map(values, columns(dv))...) :
        lift_impl(sig.fn, sig.restype, length(dv), columns(dv)...)

@generated function plain_lift_impl{T}(
        fn::Function,
        otype::Type{T},
        len::Int,
        args::AbstractVector...)
    ar = length(args)
    argvals_vars = ((Symbol("argvals", i) for i = 1:ar)...)
    init = :()
    for i = 1:ar
        init = quote
            $init
            $(argvals_vars[i]) = args[$i]
        end
    end
    return quote
        $init
        ($(argvals_vars...),) = args
        vals = Vector{T}(len)
        @inbounds for k = 1:len
            vals[k] = fn($((:($argvals_var[k]) for argvals_var in argvals_vars)...))
        end
        return PlainColumn(vals)
    end
end

@generated function lift_impl{T}(
        fn::Function,
        otype::Type{T},
        len::Int,
        args::Column...)
    optional = any(isoptional, args)
    plural = any(isplural, args)
    ar = length(args)
    argoffs_vars = ((Symbol("argoffs", i) for i = 1:ar)...)
    argvals_vars = ((Symbol("argvars", i) for i = 1:ar)...)
    idx_vars = ((Symbol("idx", i) for i = 1:ar)...)
    val_vars = ((Symbol("val", i) for i = 1:ar)...)
    init = :()
    for i = 1:ar
        init = quote
            $init
            $(argoffs_vars[i]) = offsets(args[$i])
            $(argvals_vars[i]) = values(args[$i])
        end
    end
    loop = quote
        vals[n] = fn($((val_vars[i] for i = 1:ar)...))
        n += 1
    end
    for i = ar:-1:1
        loop = quote
            for $(idx_vars[i]) = $(argoffs_vars[i])[k]:$(argoffs_vars[i])[k+1]-1
                $(val_vars[i]) = $(argvals_vars[i])[$(idx_vars[i])]
                $loop
            end
        end
    end
    return quote
        $init
        size = 0
        for k = 1:len
            size += *(1, $((:($(argoffs_vars[i])[k+1]-$(argoffs_vars[i])[k]) for i = 1:ar)...))
        end
        offs = Vector{Int}(len+1)
        vals = Vector{T}(size)
        offs[1] = 1
        if size == 0
            for k = 1:len
                offs[k+1] = 1
            end
            return Column{$optional,$plural}(offs, vals)
        end
        n = 1
        for k = 1:len
            $loop
            offs[k+1] = n
        end
        return Column{$optional,$plural}(offs, vals)
    end
end

# Nullable to optional decoder.

immutable DecodeNullableSig <: AbstractPrimitive
end

function ev{T}(::DecodeNullableSig, ivals::AbstractVector{Nullable{T}})
    len = length(ivals)
    sz = 0
    for val in ivals
        if !isnull(val)
            sz += 1
        end
    end
    offs = Vector{Int}(len+1)
    vals = Vector{T}(sz)
    offs[1] = 1
    n = 1
    for k = 1:len
        val = ivals[k]
        if !isnull(val)
            vals[n] = get(val)
            n += 1
        end
        offs[k+1] = n
    end
    return OptionalColumn(offs, vals)
end

# Vector to plural decoder.

immutable DecodeVectorSig <: AbstractPrimitive
end

function ev{T}(::DecodeVectorSig, ivals::AbstractVector{Vector{T}})
    len = length(ivals)
    sz = 0
    for val in ivals
        sz += length(val)
    end
    offs = Vector{Int}(len+1)
    vals = Vector{T}(sz)
    offs[1] = 1
    n = 1
    for k = 1:len
        for val in ivals[k]
            vals[n] = val
            n += 1
        end
        offs[k+1] = n
    end
    return PluralColumn(offs, vals)
end

function ev{V<:AbstractVector}(::DecodeVectorSig, ivals::AbstractVector{V})
    len = length(ivals)
    sz = 0
    for val in ivals
        sz += length(val)
    end
    offs = Vector{Int}(len+1)
    vals = Vector{eltype(V)}(sz)
    offs[1] = 1
    n = 1
    for k = 1:len
        for val in ivals[k]
            vals[n] = val
            n += 1
        end
        offs[k+1] = n
    end
    return PluralColumn(offs, vals)
end

# Optional to Nullable encoder.

immutable EncodeNullableSig <: AbstractPrimitive
end

function ev(::EncodeNullableSig, dv::DataVector)
    return encode_nullable_impl(column(dv, 1))
end

function encode_nullable_impl(col::Column)
    cr = cursor(col)
    T = eltype(col.vals)
    vals = Vector{Nullable{T}}(length(col))
    while !done(col, cr)
        next!(col, cr)
        vals[cr.pos] = length(cr) > 0 ? Nullable{T}(cr[1]) : Nullable{T}()
    end
    return PlainColumn(vals)
end

# Plural to Vector encoder.

immutable EncodeVectorSig <: AbstractPrimitive
end

function ev(::EncodeVectorSig, dv::DataVector)
    return encode_vector_impl(column(dv, 1))
end

function encode_vector_impl(col::Column)
    cr = cursor(col)
    T = eltype(col.vals)
    vals = Vector{Vector{T}}(length(col))
    while !done(col, cr)
        next!(col, cr)
        vals[cr.pos] = cr[1:end]
    end
    return PlainColumn(vals)
end

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
using Test, Pkg

@testset "Land" begin
    all_tests = isempty(ARGS) || "all" in ARGS ? true : false
    for submodule in ["Model"]
        if all_tests ||
           "$submodule" in ARGS ||
           "Land/$submodule" in ARGS ||
           "Land" in ARGS
            include_test(submodule)
        end
    end
end


text:
module FGenerators

# Use README as the docstring of the module:
@doc let path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    replace(read(path, String), r"^```julia"m => "```jldoctest README")
end FGenerators

export @fgenerator, @yield, @yieldfrom, Foldable

import ContextualMacros
using AbstractYieldMacros: @yield, @yieldfrom
using Base.Meta: isexpr
using FLoopsBase: with_extra_state_variables
using MacroTools: @capture, combinedef, splitdef
using Transducers: AdHocFoldable, Foldable, Transducers, foldl_nocomplete

"""
    @yield item

Yield an item from a generator.  This is usable only inside special
contexts such as within [`@fgenerator`](@ref) macro.
"""
:(@yield)

"""
    @yieldfrom foldable

Yield items from a `foldable`.  This is usable only inside special
contexts such as within [`@fgenerator`](@ref) macro.

## Examples

```jldoctest @yieldfrom
julia> using FGenerators

julia> @fgenerator function flatten2(xs, ys)
           @yieldfrom xs
           @yieldfrom ys
       end;

julia> collect(flatten2(1:2, 11:12))
4-element Array{Int64,1}:
  1
  2
 11
 12
```
"""
:(@yieldfrom)

const RF = gensym(:__rf__)
const ACC = gensym(:__acc__)

function _on_yield(ctx)
    if length(ctx.args) != 1
        throw(ArgumentError("`@yield` requires exactly one argument. Got:\n$(ctx.args)"))
    end
    x, = ctx.args
    quote
        $ACC = $Transducers.@next($RF, $ACC, $x)
    end |> esc
end

function _on_yieldfrom(ctx)
    if length(ctx.args) != 1
        throw(ArgumentError("`@yieldfrom` requires exactly one argument. Got:\n$(ctx.args)"))
    end
    foldable, = ctx.args
    quote
        $ACC = $Transducers.@return_if_reduced $foldl_nocomplete($RF, $ACC, $foldable)
    end |> esc
end

function _hack_annon!(ex)
    if isexpr(ex, :function, 2) && isexpr(ex.args[1], :tuple)
        ex.args[1] = Expr(:call, :_, ex.args[1].args...)
        # Update MacroTools.jl after anonymous function is supported:
        # https://github.com/MikeInnes/MacroTools.jl/pull/140
    elseif isexpr(ex, :->, 2) && ex.args[1] isa Symbol
        ex = Expr(:function, Expr(:call, :_, ex.args[1]), ex.args[2])
    elseif isexpr(ex, :->, 2) && isexpr(ex.args[1], :tuple)
        ex = Expr(:function, Expr(:call, :_, ex.args[1].args...), ex.args[2])
    end
    return ex
end

"""
    @fgenerator function f(...) ... end
    @fgenerator f(...) = ...
    @fgenerator(generator::GeneratorType) do; ...; end

Define a function `f` that returns a "generator" usable with
Transducers.jl-compatible APIs (akd _foldable_ interface).  In
`@fgenerator(generator::GeneratorType) do ... end` form, define
Transducers.jl interface for `GeneratorType`.

Use [`@yield`](@ref) in the function body for producing an item.  Use
`return` to finish producing items.

See also [`FGenerators`](@ref).

# Extended help
## Examples

Defining a function that returns a generator:

```jldoctest @fgenerator
julia> using FGenerators

julia> @fgenerator function generate123()
           @yield 1
           @yield 2
           @yield 3
       end;

julia> collect(generate123())
3-element Array{Int64,1}:
 1
 2
 3
```

Defining foldable interface for a pre-existing type:

```jldoctest @fgenerator
julia> struct Counting end;

julia> @fgenerator(itr::Counting) do
           i = 1
           while true
               @yield i
               i += 1
           end
       end;

julia> using Transducers  # for `Take`

julia> Counting() |> Take(3) |> collect
3-element Array{Int64,1}:
 1
 2
 3
```

Note that function such as `collect` and `sum` still dispatches to
`iterate`-based methods (above `Counting` example worked because
`Counting` was wrapped by `Take`).  To automatically dispatch to
`foldl`-based methods, use `Foldable` as the supertype:

```jldoctest @fgenerator
julia> struct Count <: Foldable
           start::Int
           stop::Int
       end;

julia> @fgenerator(itr::Count) do
           i = itr.start
           i > itr.stop && return
           while true
               @yield i
               i == itr.stop && break
               i += 1
           end
       end;

julia> collect(Count(0, 2))
3-element Array{Int64,1}:
 0
 1
 2
```
"""
macro fgenerator(ex)
    ex = _hack_annon!(ex)
    def = splitdef(ex)

    if def[:name] === :_  # see `_hack_annon!`
        def[:name] = gensym(:anonymous)
    end

    allargs = map([def[:args]; def[:kwargs]]) do x
        if isexpr(x, :kw)
            x = x.args[1]
        end
        if @capture(x, args_...)
            args
        elseif @capture(x, a_::T_)
            a::Symbol
        else
            x::Symbol
        end
    end

    body = def[:body]
    funcname = gensym(string(def[:name], "#foldl"))
    folddef = define_foldl(__module__, funcname, NamedTuple, allargs, body)
    if VERSION < v"1.4"
        namedtuple(; kw...) = kw.data
        nt = :($namedtuple(; $((Expr(:kw, a, a) for a in allargs)...)))
    else
        nt = :((; $((Expr(:kw, a, a) for a in allargs)...)))
    end
    def[:body] = :($AdHocFoldable($folddef, $nt))
    esc(combinedef(def))
end

macro fgenerator(fun, coll)
    fun = _hack_annon!(fun)
    def = splitdef(fun)
    if !@capture(coll, collvar_::typename_)
        msg =
            "`@fgenerator(...) do; ...; end` (or two-argument form of `@fgenerator`)" *
            " requires `generator::GeneratorType` as in"
        throw(ArgumentError("""
        $msg

            @fgenerator(generator::GeneratorType) do
                ...
            end

        Instead of the argument of the form `generator::GeneratorType`, got:
        $coll
        """))
    end

    if length(def[:args]) != 0
        throw(ArgumentError(
            "`@fgenerator(...) do; ...; end` (or two-argument form of" *
            " `@fgenerator`) cannot be used with arguments after `do`. Got:\n" *
            join(def[:args], ", "),
        ))
    end
    if !isempty(def[:kwargs])
        # Since `do` block can't generate keyword arguments, assume
        # that it is from two-argument form.
        throw(ArgumentError(
            "Two-argument form of `@fgenerator` requires the first argument to" *
            " be a function with a single argument. Got a function with" *
            " keyword argument(s):\n" *
            string(:(; $(def[:kwargs]...))),
        ))
    end

    return esc(define_foldl(
        __module__,
        :($Transducers.__foldl__),
        typename,
        collvar,
        def[:body],
    ))
end

function is_function(ex)
    if isexpr(ex, :(=)) || isexpr(ex, :->)
        isexpr(ex.args[1], :call) && return true
        isexpr(ex.args[1], :where) && return true
    elseif isexpr(ex, :function)
        return true
    end
    return false
end

function define_foldl(__module__::Module, funcname, structname, allargs, body)
    completion = :(return $Transducers.complete($RF, $ACC))
    function rewrite(body)
        body isa Expr || return body
        is_function(body) && return body
        isexpr(body, :meta) && return body
        if isexpr(body, :return)
            returning = get(body.args, 1, nothing)
            if returning in (:nothing, nothing)
                return completion
            end
            error("Returning non-nothing from a generator: $(body.args[1])")
        end
        return Expr(body.head, map(rewrite, body.args)...)
    end
    body = rewrite(body)
    if allargs isa Symbol
        xs = allargs
        unpack = []
    else
        @gensym xs
        unpack = [:(local $a = $xs.$a) for a in allargs]
    end
    ex = quote
        function $funcname($RF::RFType, $ACC, $xs::$structname) where {RFType}
            $(unpack...)
            $body
            $completion
        end
    end
    with_extra_state_variables([ACC]) do
        return ContextualMacros.expandwith(
            __module__,
            ex;
            yield = _on_yield,
            yieldfrom = _on_yieldfrom,
        )
    end
end

"""
    FGenerators.@with(f(...)) do ... end

Equivalent to `f(...) do ... end` but workarounds the boxing problem
when `@yield` is used inside the body of the `do` block.

!!! warning

    This macro is valid to use only for "context manager" type of `do`
    block use-case where the `do` block is executed *exactly once* and
    the result of the body is returned from `f`.  For example,
    `@with(open(filename)) do` is valid but `@with(map(xs)) do` is
    not.  This is because `map` can call the `do` block arbitrary
    number of times.
"""
macro with(doblock, call)
    @gensym acc1 acc2 ans1 ans2
    true_doblock = Expr(:block, __source__, doblock.args[end])
    doblock.args[end] = quote
        local $ACC = $acc1
        local $ans1 = $true_doblock
        ($ans1, $ACC)
    end
    calldo = Expr(:do, call, doblock)
    quote
        local $ans2, $acc2
        ($ans2, $acc2) = let $acc1 = $ACC
            $calldo
        end
        $ACC = $Transducers.@return_if_reduced $acc2
        $ans2
    end |> esc
end

end


text:
module CuckooFilters
using SHA: sha256

include("bucket.jl")
include("filter.jl")

fingerprint(x::AbstractVector{UInt8}) = sha256(x)[1]
fingerprint(x::Integer) = UInt8(x & 0xff)
fingerprint_type(::Type{<:AbstractVector{UInt8}}) = UInt8
fingerprint_type(::Type{<:Integer}) = UInt8

# Exports
export CuckooFilter
export insert_unique!

end # module


text:
struct FiniteCustomTrigger <: AbstractFiniteTrigger
    f::Function
    n::Int
end

struct InfiniteCustomTrigger <: AbstractInfiniteTrigger
    f::Function
end

"""
    CustomTrigger(f::Function[, n=number_of_times])

A trigger which should trigger a job according a function `f`.

It's generally a better idea (cleaner implementation) to write your own
trigger from  `AbstractTrigger`, `AbstractFiniteTrigger` or `AbstractInfiniteTrigger`
but passing a function to a `CustomTrigger` can be quite handy

Job can be triggered:

- once (`n=1`)
- a finite number of times (`n=number_of_times`)
- indefinitely (without setting `n`)

# Example

    f = (dt_previous_fire, dt_now) -> dt_now + Dates.Minute(5)
    trigger = CustomTrigger(f)

should run a job every 5 minutes
"""
function CustomTrigger(f; n=-1)
    if n < 0
        InfiniteCustomTrigger(f)
    else
        FiniteCustomTrigger(f, n)
    end
end

"""
    Trigger(f::Function[, n=number_of_times])

Return an `CustomTrigger` which should trigger a job according a function `f`.
"""
Trigger(f::Function; kwargs...) = CustomTrigger(f; kwargs...)

function get_next_dt_fire(trigger::Union{FiniteCustomTrigger,InfiniteCustomTrigger}, dt_previous_fire, dt_now)
    trigger.f(dt_previous_fire, dt_now)
end


text:
using UnifiedPseudopotentialFormat
using Test

@testset "UnifiedPseudopotentialFormat.jl" begin
    include("name.jl")
    include("PSlibrary.jl")
end


text:
include("..\\libs\\FunctionsRobot.jl")

function search_out!(r::Robot)
    side = West
    counter = 1

    while isborder(r, Nord)
        moves!(r, side, counter)
        side = reversSide(side)
        counter += 1
    end
end


text:
# ---- Quadrature Methods ---------

"""
    trapz(x, y)

Integrate y w.r.t. x using the trapezoidal method.
"""
function trapz(x, y)

    integral = 0.0
    for i = 1:length(x)-1
        integral += (x[i+1]-x[i])*0.5*(y[i] + y[i+1])
    end
    return integral
end


text:
# This file was generated, do not modify it. # hide
x = [1, 3, 2, 5]
@show x
@show length(x)


text:
# This file is a part of RunStatistics.jl, licensed under the MIT License (MIT).

using RunStatistics
using Test

@testset "t_obs" begin

    @test RunStatistics.t_obs([-1, 1, 3, -2], 0, 1) == (10.0, [2,3])
    @test RunStatistics.t_obs([-1, 1, 3, -2], zeros(4), ones(4)) == (10.0, [2,3])

    @test RunStatistics.t_obs([-1, 1, 3, -2, 1, 3], 0, 1) == (10.0, [[2,3],[5, 6]])
    @test RunStatistics.t_obs([-1, 1, 3, -2, 1, 3], zeros(6), ones(6)) == (10.0, [[2,3],[5, 6]])
end


text:
# -------------------------------------------------------------------------------------------
# X. Gandibleux - Metaheuristiques : TSP TS (version etudiant) - Septembre 2019
# Base de depart « etudiant » pour exercice "TS applique au TSP"

@static if VERSION < v"1.0-"
    error("NOT COMPLIANT WITH JULIA < v1.0.0")
end

using Random
using Printf
using PyPlot # a commenter si le tracage des resultats n'est pas appele

# -------------------------------------------------------------------------------------------
# construit un distancier symetrique aux valeurs de distances uniformément aléatoires dans 1..vSup
function tsp_elaboreDistancierAleatoire(
           n::Int64,    # in : nombre de villes
           vSup::Int64  # in : plus grande distance entre deux villes
         )
    distancier = rand(1:vSup,n,n)
    for i=1:n
        for j=1:n
    		    distancier[i,j] = distancier[j,i]
        end
        distancier[i,i] = 0
    end
    return distancier # out : matrice carree symetrique de distances entieres
end

# -------------------------------------------------------------------------------------------
# calcule la longueur d’un circuit hamiltonien correspondant a une solution TSP admissible
function tsp_calculeLongueurCircuit(
           d::Matrix{Int64},  # in : distancier symetrique de ville a ville
           x::Array{Int64},   # in : solution admissible représentée par codage de permutation des villes
         )
    n::Int64 = length(x)
    longueur::Int64 = d[ x[n] , x[1] ]
    for i = 1:n-1
        longueur += d[ x[i] , x[i+1] ]
    end
    return longueur # out : longueur du circuit correspondant a la solution admissible
end

# -------------------------------------------------------------------------------------------
# construit une solution aléatoire TSP admissible représentée dans un codage de permutation
function tsp_elaboreSolutionAleatoire(
           n::Int64  # in : nombre de villes
         )
    return shuffle(collect(1:n))  # out : une permutation aleatoire des villes
end

# ------------------------------------------------------------
# Plot a graph with the results obtained by the SA (z for all solutions, best, greedy)

function traceResultat(titre, zAllLst, zBestLst, zBestAllLst)
    figure("Examen d'un run",figsize=(6,6))
    title(titre)
    xlabel("Itérations")
    ylabel("valeurs de z(x)")
    grid()
    plot(zAllLst,label="all solutions")
    plot(zBestAllLst, linewidth=2.0, color="green", label="best solutions")
    legend(loc=1, fontsize ="small")
end

# -------------------------------------------------------------------------------------------
function tsp_solveTS(
           d::Matrix{Int64},      # in : distancier symetrique de ville a ville
           x::Array{Int64},       # in : solution TSP admissible représentée pas un codage de permutation des villes
           z::Int64,              # in : valeur de la solution TSP admissible
           iterTabou::Int64,      # in : nombre d'iterations a realiser par la recherche
           dureeTabou::Int64      # in : valeur de la longueur tabou
         )

    # traces d'activite pour le trace de graphiques
    zAllLst     = [] # all values of f(x) collected for each iteration
    zBestLst    = [] # each best values of f(x) collected
    zBestAllLst = [] # best values of f(x) collected for each iteration
    push!(zAllLst, z)
    push!(zBestLst, z)
    push!(zBestAllLst, z) ; zMin = z

    @printf(" parametres TS : iterTabou= %d, dureeTabou= %d\n\n", iterTabou,  dureeTabou)

    n      = length( x )
    xCur   = copy( x ) # solution courante
    zCur   = z
    xBest  = copy( x ) # meilleure solution
    zBest  = z

    # Initiallisation memoire tabou avec  "tous les mouvements sont autorises"
    listeTabou  = zeros(Int64,n,n)

    # Boucle principale de la TS
    for nr_iteration=1:iterTabou

        trouveAmeliorant = false
        deltaBest = typemax(Int64)
        iDbest = 0  ;  iFbest = 0
        jDbest = 0  ;  jFbest = 0

        # recherche meilleur mouvement autorise en appliquant un 2opt sur un voisinage complet
        for i=0:n-1
            for j=0:n-4

                iD = i%n + 1       ;  iF = (i+1)%n + 1
                jD = (i+j+2)%n + 1 ;  jF = (i+j+3)%n + 1

                d1 = d[ xCur[iD] , xCur[iF] ] + d[ xCur[jD] , xCur[jF] ]
                d2 = d[ xCur[iD] , xCur[jD] ] + d[ xCur[iF] , xCur[jF] ]
                delta = d2 - d1

                voisinTrouve = ( delta  <  deltaBest )
                mouvementAutorise = (listeTabou[ xCur[iD] , xCur[jD] ] <= nr_iteration) && (listeTabou[ xCur[iF] , xCur[jF] ] <= nr_iteration)
                critereAspiration = ( zCur + delta < zBest )

                if ( voisinTrouve  && ( mouvementAutorise || critereAspiration ))
                    trouveAmeliorant = true
                    deltaBest = delta
                    iDbest = iD  ;  iFbest = iF
                    jDbest = jD  ;  jFbest = jF
                end
            end
        end

        if trouveAmeliorant

            # maj memoire tabou
            listeTabou[ xCur[iDbest] , xCur[iFbest] ] = dureeTabou + nr_iteration
            listeTabou[ xCur[jDbest] , xCur[jFbest] ] = dureeTabou + nr_iteration
            listeTabou[ xCur[iFbest] , xCur[iDbest] ] = dureeTabou + nr_iteration
            listeTabou[ xCur[jFbest] , xCur[jDbest] ] = dureeTabou + nr_iteration

            # maj incrementale de la valeur de la solution
            zCur = zCur + deltaBest
            push!(zAllLst,zCur)

            # maj solution : 2opt et inverse les segments du chemin entre les 2 segments traites
            if iDbest < jDbest
                xCur[iDbest+1:jDbest]= reverse(xCur[iDbest+1:jDbest])
            else
                xtwice=vcat(xCur, xCur)
                xtwice[iDbest+1:jDbest+n]= reverse(xtwice[iDbest+1:jDbest+n])
                xCur[1:jDbest] = xtwice[n+1:n+jDbest]
                xCur[iDbest:n] = xtwice[iDbest:n]
            end

            # met a jour la meilleure solution si il y a eu amelioration
            if (zCur < zBest)
                zBest = zCur
                xBest = copy(xCur)
                @printf("  numIter = %3d    z(xBest) = %4d\n", nr_iteration, zBest)
                push!(zBestLst,zBest)
                zMin = zBest
            end

        else
            @assert false "Tous les mouvements sont interdits : liste tabou trop longue"
        end

        push!(zBestAllLst,zMin)
    end # indice sur nr_iteration

    return xBest, zBest, zAllLst, zBestLst, zBestAllLst

end

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

function tsp_TS(
           n::Int64            = 20,  # in : nombre de villes
           vSup::Int64         = 25, # in : plus grande distance entre deux villes
           iterTabou::Int64    = 60, # in : nombre d'iterations a realiser par la recherche
           dureeTabou::Int64   = 7   # in : valeur de la longueur tabou
         )

    d = tsp_elaboreDistancierAleatoire(n,vSup) # distancier aleatoire

    xInit = tsp_elaboreSolutionAleatoire(n)
    zInit = tsp_calculeLongueurCircuit(d,xInit)

    println("Init    : z(xInit) = ", zInit)
    println("  ")

    xBest, zBest, zAllLst, zBestLst, zBestAllLst = tsp_solveTS(d, xInit, zInit, iterTabou, dureeTabou)
    println("\nTS      : z(xTS)   = ", zBest)
    println("  ")

    traceResultat("TS : TSP | n=" * string(20) , zAllLst, zBestLst, zBestAllLst)

end

max_stars_repo_name: xgandibleux/metaTSPstu

max_issues_repo_name: xgandibleux/metaTSPstu

max_forks_repo_name: xgandibleux/metaTSPstu