text:
module neuro
using Random

include("activation_functions.jl")
include("sample_container.jl")
include("error_measure.jl")
include("training_options.jl")
include("tools.jl")

# =========================================================================== #
# Neural Network
# =========================================================================== #

struct NeuralNetwork{ErrorFunction}
    n_neurons_per_layers::Vector{Int}
    activation_function::Vector{AbstractActivationFunction}
    error_function::ErrorFunction
    weight_matrices::WeightMatrix{Float64}
    biases::Bias{Float64}

    function NeuralNetwork{ErrorFunction}(
        n_neurons_per_layers::Vector{Int},
        activation_function,
        error_function::ErrorFunction = LeastSquaresErrorFunction()
        ) where     { ErrorFunction <: AbstractErrorFunction}

        weight_matrices = WeightMatrix(n_neurons_per_layers)
        biases          = Bias(n_neurons_per_layers)

        return new(
            n_neurons_per_layers,
            activation_function,
            error_function,
            weight_matrices,
            biases)
    end
end

function NeuralNetwork(
    n_neurons_per_layers::Vector{Int},
    activation_function,
    error_function::ErrorFunction = LeastSquaresErrorFunction()
    ) where {ErrorFunction <: AbstractErrorFunction}

    return NeuralNetwork{ErrorFunction}(
        n_neurons_per_layers,
        activation_function,
        error_function)
end

function copy_in_bias!(nn::NeuralNetwork, b_vec, index::Int)
    for kk = 1:length(nn.biases.array[index])
        nn.biases.array[index][kk] = b_vec[kk]
    end
end

function copy_in_weight_matrix!(nn::NeuralNetwork, weight_mat, index::Int)
    for kk = 1:length(nn.weight_matrices.array[index])
        nn.weight_matrices.array[index][kk] = weight_mat[kk]
    end
end

function n_layers(nn::NeuralNetwork)
    return length(nn.n_neurons_per_layers)
end

function fw_prop(nn::NeuralNetwork, x_input)
    N_layers = n_layers(nn)
    y        = x_input
    for layer = 1:(N_layers-1)
        W = nn.weight_matrices.array[layer]
        # u = W * y
        # x = u + nn.biases.array[layer]
        # y = nn.activation_function[layer].(x)
        y = nn.activation_function[layer].( W * y + nn.biases.array[layer] )
    end

    return y
end

function fw_prop(nn::NeuralNetwork, x_input, x_vec)
    N_layers = n_layers(nn)
    y        = x_input
    for layer = 1:(N_layers-1)
        W = nn.weight_matrices.array[layer]
        u = W * y
        x_vec.array[layer] = u + nn.biases.array[layer]
        y = nn.activation_function[layer].(x_vec.array[layer])
    end

    return y
end

function bw_prop(nn::NeuralNetwork, x_vec, delta_e)
    N_layers = n_layers(nn)
    dy = delta_e
    for layer = N_layers:-1:2
        act_fcn = nn.activation_function[layer-1]

        Ldy = length(dy)
        ds  = zeros(Ldy)
        for ii = 1:Ldy
            ds[ii] = differentiate(act_fcn, x_vec.array[layer-1][ii])
        end
        dx = dy .* ds
        W  = nn.weight_matrices.array[layer-1]
        dy = W' * dx
    end
    return dy
end

function train2!(nn::NeuralNetwork, sample_container::SC, options::TrainingOptions) where {SC <: AbstractSampleContainer}

    # generals
    N_layers  = n_layers(nn)
    N_samples = n_samples(sample_container)

    # create indices
    sample_indices = collect(1:N_samples)

    # pre allocate containers
    x_vec        = IntermediateVectorResult(nn.n_neurons_per_layers)
    y_vec        = IntermediateVectorResult(nn.n_neurons_per_layers)
    delta_x      = IntermediateVectorResult(nn.n_neurons_per_layers)
    delta_y      = IntermediateVectorResult(nn.n_neurons_per_layers)
    delta_bias   = IntermediateVectorResult(nn.n_neurons_per_layers)
    dactifcn     = IntermediateVectorResult(nn.n_neurons_per_layers)
    delta_w_mat  = WeightMatrix(nn.n_neurons_per_layers)
    delta_e      = zeros(nn.n_neurons_per_layers[end])

    iter = 0
    for outer iter = 1:options.max_iter

        # shuffle the indices
        shuffle!(sample_indices)

        for bb = 1:options.batch_size

            # select the sample index
            sample_index = sample_indices[bb]

            # forward propagate and save the inputs
            (x_meas, y_meas) = get_sample(sample_container, sample_index)

            y = x_meas
            layer = 0
            for outer layer = 1:(N_layers-1)
                W = nn.weight_matrices.array[layer]
                x_vec.array[layer] = W * y + nn.biases.array[layer]
                y_vec.array[layer] = nn.activation_function[layer].(x_vec.array[layer])
                y                  = y_vec.array[layer]
            end

            # evaluate the derivative of the error function
            y_out = y_vec.array[layer]
            gradient!(nn.error_function, y_meas, y_out, delta_e)

            # backward propagate the errors and update the deltas
            delta_y.array[layer] = delta_e
            for layer = (N_layers-1):-1:1
                act_fcn = nn.activation_function[layer]

                Ldy = length(x_vec.array[layer])
                for ii = 1:Ldy
                    dactifcn.array[layer][ii] = differentiate(act_fcn, x_vec.array[layer][ii])
                end
                delta_x.array[layer] = delta_y.array[layer] .* dactifcn.array[layer]
                W  = nn.weight_matrices.array[layer]
                if layer >= 2
                    delta_y.array[layer-1] = W' * delta_x.array[layer]
                end
            end

            for layer = 1:N_layers-1
                delta_bias.array[layer]  += delta_x.array[layer]
            end
            delta_w_mat.array[1] += delta_x.array[1] * x_meas'
            for layer = 2:N_layers-1
                delta_w_mat.array[layer] += delta_x.array[layer] * y_vec.array[layer-1]'
            end

        end # batch

        for layer = 1:N_layers-1
            # max norm
            # options.max_norm
            # n = norm()
            nn.biases.array[layer]          -= options.learning_rate .* delta_bias.array[layer]
            nn.weight_matrices.array[layer] -= options.learning_rate .* delta_w_mat.array[layer]
        end

        zeros!(delta_bias)
        zeros!(delta_w_mat)

    end # iter

    return nn
end

# function evaluate_container(nn::NeuralNetwork, cv_cont::MatrixVectorSampleContainer)
function evaluate_container(nn::NeuralNetwork, cv_cont)

    N_layers  = n_layers(nn)
    N_samples = n_samples(cv_cont)
    # N_dim     = n_dim(cv_cont)
    y_out     = zeros(nn.n_neurons_per_layers[end], N_samples)

    for kk = 1:N_samples
        (xm, ym) = get_sample(cv_cont, kk)
        y = xm
        for layer = 1:(N_layers-1)
            W = nn.weight_matrices.array[layer]
            y = nn.activation_function[layer].(W * y .+ nn.biases.array[layer])
        end
        y_out[:, kk] = y
    end

    return y_out
end

# function evaluate_performance(nn::NeuralNetwork, cv_cont::MatrixVectorSampleContainer)
function evaluate_performance(nn::NeuralNetwork, cv_cont)

    N_layers  = n_layers(nn)
    N_samples = n_samples(cv_cont)
    y_out     = zeros(N_samples)

    for kk = 1:N_samples
        (xm, ym) = get_sample(cv_cont, kk)
        y = xm
        for layer = 1:(N_layers-1)
            W = nn.weight_matrices.array[layer]
            y = nn.activation_function[layer].(W * y .+ nn.biases.array[layer])
        end
        # p = nn.error_function(ym, y[1:nn.n_neurons_per_layers[end]])
        p = nn.error_function(ym, y)
        y_out[kk] = p
    end

    return y_out
end

#=
function train!(nn::NeuralNetwork, sample_container::SC, options::TrainingOptions) where {SC <: AbstractSampleContainer}

    # generals
    N_layers  = n_layers(nn)
    N_samples = n_samples(sample_container)

    # create indices
    sample_indices = collect(1:N_samples)

    # pre allocate containers
    x_vec        = IntermediateVectorResult(nn.n_neurons_per_layers)
    y_vec        = IntermediateVectorResult(nn.n_neurons_per_layers)
    delta_x      = IntermediateVectorResult(nn.n_neurons_per_layers)
    delta_y      = IntermediateVectorResult(nn.n_neurons_per_layers)
    delta_bias   = IntermediateVectorResult(nn.n_neurons_per_layers)
    dactifcn     = IntermediateVectorResult(nn.n_neurons_per_layers)
    delta_w_mat  = WeightMatrix(nn.n_neurons_per_layers)
    delta_e      = zeros(nn.n_neurons_per_layers[end])

    iter = 0
    for outer iter = 1:options.max_iter

        # shuffle the indices
        shuffle!(sample_indices)

        for bb = 1:options.batch_size

            # select the sample index
            sample_index = sample_indices[bb]

            # forward propagate and save the inputs
            (x_meas, y_meas) = get_sample(sample_container, sample_index)

            y = x_meas
            for layer = 1:(N_layers-2)
                W                    = nn.weight_matrices.array[layer]
                x_vec.array[layer]   = W * y
                x_vec.array[layer] .+= nn.biases.array[layer]
                y_vec.array[layer]   = nn.activation_function.(x_vec.array[layer])
                y                    = y_vec.array[layer]
            end

            # output layer
            layer                = N_layers-1
            W                    = nn.weight_matrices.array[layer]
            x_vec.array[layer]   = W * y_vec.array[layer-1]
            x_vec.array[layer] .+= nn.biases.array[layer]
            y_vec.array[layer]   = nn.output_activation_function.(x_vec.array[layer])

            # evaluate the derivative of the error function
            y_out = y_vec.array[layer]
            gradient!(nn.error_function, y_meas, y_out, delta_e)

            # backward propagate the errors and update the deltas
            delta_y.array[layer] = delta_e

            # output layer
            for ii = 1:nn.n_neurons_per_layers[end]
                dactifcn.array[layer][ii] = differentiate(nn.output_activation_function, x_vec.array[layer][ii])
            end
            delta_x.array[layer]      = delta_y.array[layer] .* dactifcn.array[end]

            for layer = (N_layers-1):-1:2
                W                      = nn.weight_matrices.array[layer]
                delta_y.array[layer-1] = W' * delta_x.array[layer]
                for ii = 1:nn.n_neurons_per_layers[layer-1]
                    dactifcn.array[layer-1][ii] = differentiate(nn.activation_function, x_vec.array[layer-1][ii])
                end
                delta_x.array[layer-1]    = delta_y.array[layer-1] .* dactifcn.array[layer-1]
                delta_bias.array[layer]  .+= delta_x.array[layer]

                delta_w_mat.array[layer] .+= delta_x.array[layer] * y_vec.array[layer-1]'
            end
            delta_bias.array[1]  .+= delta_x.array[1]
            delta_w_mat.array[1] .+= delta_x.array[1] * x_meas'

        end

        # update the weights
        # b <- b - alpha * delta_x
        # W <- W - alpha * delta_W
        for layer = 1:N_layers-1
            bias        = nn.biases.array[layer]
            bias       .-= 1/options.batch_size .* options.learning_rate .* delta_bias.array[layer]
            weight_mat  = nn.weight_matrices.array[layer]
            weight_mat .-= 1/options.batch_size .* options.learning_rate .* delta_w_mat.array[layer]
        end

        # here could be the performance evaluated

    end

    return (nn, iter)
end
=#

#=
function evaluate(nn::NeuralNetwork, x_input)
    N_layers = n_layers(nn)
    y        = x_input
    for layer = 1:(N_layers-2)
        W = nn.weight_matrices.array[layer]
        y = nn.activation_function.(W * y .+ nn.biases.array[layer])
    end
    # output layer
    layer = N_layers-1
    W     = nn.weight_matrices.array[layer]
    y     = nn.output_activation_function.(W * y .+ nn.biases.array[layer])
    return y
end
=#

end # module

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
using SnoopCompile

@snoopi_bot "BinaryBuilder" begin
    using BinaryBuilder
    # Do an actual build
    products = Product[
        ExecutableProduct("hello_world_c", :hello_world_c),
        ExecutableProduct("hello_world_cxx", :hello_world_cxx),
        ExecutableProduct("hello_world_fortran", :hello_world_fortran),
        ExecutableProduct("hello_world_go", :hello_world_go),
        ExecutableProduct("hello_world_rust", :hello_world_rust),
    ]

    # First, do the build, but only output the meta json, since we definitely want that to be fast
    build_tarballs(
        ["--meta-json=/dev/null"],
        "testsuite",
        v"1.0.0",
        # No sources
        DirectorySource[],
        # Build the test suite, install the binaries into our prefix's `bin`
        raw"""
        # Build testsuite
        make -j${nproc} -sC /usr/share/testsuite install
        # Install fake license just to silence the warning
        install_license /usr/share/licenses/MIT
        """,
        [HostPlatform()],
        products,
        # No dependencies
        Dependency[];
    )

    # Next, actually do a build, since we want that to be fast too.
    build_tarballs(
        ["--verbose"],
        "testsuite",
        v"1.0.0",
        # Add some sources, so that we actually download them
        [
            ArchiveSource("https://github.com/staticfloat/small_bin/raw/master/socrates.tar.gz",
                          "e65d2f13f2085f2c279830e863292312a72930fee5ba3c792b14c33ce5c5cc58"),
            DirectorySource("src"),
        ],
        # Build the test suite, install the binaries into our prefix's `bin`
        raw"""
        # Build testsuite
        make -j${nproc} -sC /usr/share/testsuite install
        # Install fake license just to silence the warning
        install_license /usr/share/licenses/MIT
        """,
        [HostPlatform()],
        products,
        # Add a dependency on Zlib_jll, our favorite scapegoat
        Dependency[
            Dependency("Zlib_jll"),
        ];
        compilers=[:c, :rust, :go],
    )

    rm("build"; recursive=true, force=true)
    rm("products"; recursive=true, force=true)
end

ext: jl


text:
input = readlines("day13/input.txt")
# input = readlines("day13/small.txt")
depart_time = parse(Int, input[1])
buses = convert.(String, split(input[2], ","))

min_mod, min_bus = 10000000000, 0
for bus in buses
    global min_mod,min_bus
    if tryparse(Int, bus) !== nothing
        minimize = parse(Int, bus) - mod(depart_time, parse(Int, bus))
        if minimize < min_mod
            min_mod = minimize
            min_bus = parse(Int, bus)
        end
    end
end

println(min_mod)
println(min_bus)

println(min_mod*min_bus)


text:
"""
...
"""
process_coarse!(problem, level) = nothing

"""
...
"""
process_fine!(problem, level) = nothing

"""
...
"""
process_uncoarse!(problem, level) = nothing

"""
...
"""
coarsen!(problem, level) = nothing

"""
...
"""
uncoarsen!(problem, level) = nothing

"""
...
"""
doinitial(problem, level) = false

"""
...
"""
initial!(problem, level) = nothing

"""
...
"""
process_initial!(problem, level) = nothing

"""
...
"""
precycle!(problem, level) = nothing

"""
...
"""
postcycle!(problem, level) = nothing

"""
...
"""
preprocessing!(problem, level) = nothing

"""
...
"""
postprocessing!(problem, level) = nothing


text:
"""
    first_order_model(X, I; kwargs...)

Returns a `BernoulliCodeDistribution` with `p_i` equal to the expected value of
bit `i` in `X[I,:]`.

"""
function first_order_model(X, I=1:size(X,1); kwargs...)
    mu = sum(X[I,:],2) / size(X,2)
    return BernoulliCodeDistribution(mu[:]; indices=I, autocomment="first_order_model", kwargs...)
end

"""
    second_order_model(X, I; kwargs...)

Returns an `IsingDistribution` which is fit to the pairwise correlations in `X[I,:]`.

Some useful keywords particular to this function:
 * `J0` is the initial value used for optimization

keyword argument `algorithm` sets the algorithm:
 * `algorithm = LBFGS()` default is to use the LBFGS algorithm in the `Optim` package.
 * `algorithm = :LD_LBFGS` uses LBFGS algorithm in `NLopt`
 * `algorithm = :naive` uses the function `gradient_optimizer` in `optimizers.jl`
 * `algorithm = :LD_MMA` uses the MMA algorithm as described in the `NLopt` package. Don't use this one it's hella slow.
"""
function second_order_model(X, I=1:size(X,1);
    verbose=0, algorithm=LBFGS(), kwargs...)
    # algorithm=:LD_LBFGS

    if algorithm == :naive
        return Naive_second_order_model(X, I; verbose=verbose, kwargs...)
    # elseif typeof(algorithm) == Symbol
    #     return _NLopt_second_order_model(X, I; verbose=verbose, algorithm=algorithm, kwargs...)
    else
        return Optim_second_order_model(X, I; verbose=verbose, algorithm=algorithm, kwargs...)
    end
end
function Naive_second_order_model(X, I=1:size(X,1); verbose=0, kwargs...)
    dkwargs = Dict(kwargs)
    Xselected = X[I,:]
    N_neurons,N_samples = size(Xselected)

    Jseed = rand(N_neurons,N_neurons); Jseed = (Jseed + Jseed') / (2 * N_neurons)
    Jseed = pop!(dkwargs, :J0, Jseed)
    mu = Xselected * Xselected' / N_samples
    L_X(J,g) = loglikelihood(Xselected, J, g; mu_X = mu)
    K_X(J,g) = MPF_objective(Xselected, J, g)
    fun = (N_neurons > ISING_METHOD_THRESHOLD || pop!(dkwargs, :force_MPF, false)) ? "MPF" : "loglikelihood"
    objective = (fun == "loglikelihood") ? :max : :min

    if verbose > 0
        println("second_order_model[gradient_optimizer]: setting $objective objective $fun")
    end
    # verbosity = verbose + pop!(dkwargs, :more_verbose, false)
    if fun == "loglikelihood"
        if verbose > 0
            println("\tApproximate cost: 2^$N_neurons + $N_samples per evalution")
        end
        (F_opt, J_opt, stop) = gradient_optimizer(L_X, Jseed[:]; objective=objective, verbose=verbose, dkwargs...)
    else
        if verbose > 0
            println("\tApproximate cost: $N_neurons * $N_samples per evaluation")
        end
        (F_opt, J_opt, stop) = gradient_optimizer(K_X, Jseed[:]; objective=objective, verbose=verbose, dkwargs...)
    end
    J_opt = reshape(J_opt, N_neurons, N_neurons)
    return IsingDistribution(J_opt; indices=I, autocomment="second_order_model[gradient_optimizer|$fun]", opt_val=F_opt, opt_ret=stop, dkwargs...)
    # return (F_opt, J_opt, stop, Jseed, mu, (fun == "loglikelihood" ? L_X : K_X))
end
function Optim_second_order_model(X, I=1:size(X,1); verbose=0, res=[], kwargs...)
    dkwargs = Dict(kwargs)
    _X = sortcols(X[I,:]) # sorting columns for branch prediction and memory access
    N_neurons,N_samples = size(_X)

    # J_prev_K = similar(J0)
    # J0_K = similar(J0)
    # J0_K[:] = J0[:]
    # J0_K[1:(N_neurons+1):end] = 0.0
    # buf_K = [_X, 2 * _X - 1, similar(_X, Float64), diag(J0), J0_K]
    # K_X(J) = K_MPF(J, J_prev_K, buf_K)
    # dK_X!(G, J) = dK_MPF!(G, J, J_prev_K, buf_K)

    alg = pop!(dkwargs, :algorithm, LBFGS())

    # Configure various options for Optim.jl
    # show_trace = pop!(dkwargs, :show_trace, verbose >= 2)
    options = Optim.Options(
        show_trace = pop!(dkwargs, :show_trace, verbose >= 2),
        x_tol = pop!(dkwargs, :x_tol, 0.0),
        f_tol = pop!(dkwargs, :f_tol, 0.0),
        allow_f_increases = pop!(dkwargs, :allow_f_increases, false),
        iterations = pop!(dkwargs, :iterations, 500),
        show_every = pop!(dkwargs, :show_every, 10),
        store_trace = pop!(dkwargs, :store_trace, false)
        )

    if verbose > 0
        println("second_order_model[Optim/$(summary(alg))]: preparing to run optimization")
        # println("  objective: $fun")
        print("  algorithm: $(summary(alg))")
        if typeof(alg) <: LBFGS
            println(" m = $(alg.m)")
        else
            println()
        end
        println("  N_neurons: $(size(_X,1))")
        println("  N_samples: $(size(_X,2))")
        #TODO potentially other information should be printed
    end

    # println("Using Minimum Probability Flow to estimate parameters...")
    # tic()
    # res_K = optimize(K_X, dK_X!, J0, alg, options)
    # t_K = toq()
    # J_K = Optim.minimizer(res_K)
    # println("Done! MPF $(Optim.converged(res_K) ? "did" : "did not") converge in $t_K s")
    # show(res_K)
    # println()

    J0 = -rand(N_neurons,N_neurons)
    J0 = (J0 + J0')/2
    J0 = pop!(dkwargs, :J0, J0)
    W0 = lt2v(J0)

    W_prev_L = similar(W0)
    mu_X = _X * _X' / N_samples
    buf_L = [_X, mu_X, IsingDistribution(J0), similar(W0)]
    L_X(W) = negloglikelihood(W, W_prev_L, buf_L)
    dL_X!(G, W) = dnegloglikelihood!(G, W, W_prev_L, buf_L)
    HL_X!(H, W) = Hnegloglikelihood!(H, W, W_prev_L, buf_L)

    if verbose > 0
        println("Maximizing likelihood...")
    end
    tic()
    res_L = optimize(L_X, dL_X!, HL_X!, W0, alg, options)
    t_L = toq()
    push!(res, res_L)
    println("Done! MLE $(Optim.converged(res_L) ? "did" : "did NOT") converge in $t_L s")
    show(res_L)
    println()
    # J_opt = reshape(Optim.minimizer(res_L), N_neurons, N_neurons)
    J_opt = v2lt(Optim.minimizer(res_L), N_neurons)
    P2 = IsingDistribution(J_opt;
        indices=I,
        autocomment="second_order_model[Optim/$(summary(alg))|MLE]",
        # MPF_converged=Optim.converged(res_K),
        MLE_converged=Optim.converged(res_L),
        # MPF_iterations=Optim.iterations(res_K),
        MLE_iterations=Optim.iterations(res_L),
        MLE_iteration_limit_reached=Optim.iteration_limit_reached(res_L),
        # J0_K=Optim.initial_state(res_K),
        # J0_L=Optim.initial_state(res_L),
        dkwargs...)
    # hide_metadata!(P2, :J0_L)
    # hide_metadata!(P2, :J0_K)

    return P2
end

"""
    data_model(X, I; kwargs...)

Returns a `DataDistribution` based on the specified data and indices.
"""
function data_model(X::Union{Matrix{Bool},BitMatrix}, I=1:size(X,1); kwargs...)
    return DataDistribution(X[I,:]; autocomment="data_model", indices=I, kwargs...)
end

max_stars_repo_name: sekunder/BinaryVectorProbability

max_issues_repo_name: sekunder/BinaryVectorProbability

max_forks_repo_name: sekunder/BinaryVectorProbability


text:
#========================================================================================#
#	Laboratory 500
#
# Welcome to course 500: An Introduction to Agent-Based Systems!
#
# Author: Niall Palfreyman, 24/04/2022
#========================================================================================#
[
	Activity(
		"""
		???
		""",
		"???",
		x -> true
	),
]


text:
# AUTO GENERATED FILE - DO NOT EDIT

export notificationhandler

"""
    notificationhandler(;kwargs...)

A NotificationHandler component.
Mantine notifications system. For more information, see: https://mantine.dev/others/notifications/
Keyword arguments:
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks
- `task` (optional): Task for notification handler along with notification props. task has the following type: lists containing elements 'command', 'id', 'props'.
Those elements have the following types:
  - `command` (a value equal to: "hide", "show", "update"; required)
  - `id` (String; required)
  - `props` (optional): . props has the following type: lists containing elements 'color', 'style', 'title', 'loading', 'message', 'autoClose', 'disallowClose'.
Those elements have the following types:
  - `color` (String; optional)
  - `style` (Dict; optional)
  - `title` (String; optional)
  - `loading` (Bool; optional)
  - `message` (String; optional)
  - `autoClose` (Real | a value equal to: false; optional)
  - `disallowClose` (Bool; optional)
"""
function notificationhandler(; kwargs...)
        available_props = Symbol[:id, :task]
        wild_props = Symbol[]
        return Component("notificationhandler", "NotificationHandler", "dash_mantine_components", available_props, wild_props; kwargs...)
end


text:
function numb!(r)
    function change(side, num)
        return HorizonSide((Int(side) + num) % 4)
    end

    function go_around(side)
        num = 0
        while isborder(r, side)
            if !isborder(r, change(side, 1))
                move!(r, change(side, 1))
                num += 1
            else
                brea
            end
        end
        if isborder(r, side)
            for i in 1:num
                if !isborder(r, change(side, 3))
                    move!(r, change(side, 3))
                end
            end
            return false
        else
            move!(r, side)
            for i in 1:num
                if !isborder(r, change(side, 3))
                    move!(r, change(side, 3))
                end
            end
            return true
        end
    end

    function move_and_num(side, num)
        flag = false
        ok = true
        count = 0
        while true
            if (isborder(r, change(side, num)) && !flag)
                flag = true
                count += 1
            elseif !isborder(r, change(side, num))
                flag = false
            end
            if (isborder(r, side))
                ok = go_around(side)
            end
            if ok
                if !isborder(r, side)
                    move!(r, side)
                end
            else
                break
            end
        end
        return count
    end

    function num(side, move_side, ok)
        count = 0
        n = 1
        if ok
            n = 1
        else
            n = 3
        end
        while !isborder(r, side)
            count += move_and_num(move_side, n)
            move_side = change(move_side, 2)
            n = (n + 2) % 4
            move!(r, side)
        end
        return count
    end

    function solve()
        count = 0
        count += num(Nord, Ost, true)
        while !isborder(r, Sud)
            move!(r, Sud)
        end
        count += num(Ost, Nord, false)
        println(count)
    end

    solve()
end

max_stars_repo_name: dmitrii-a-piskunov/mirea-progs

max_issues_repo_name: dmitrii-a-piskunov/mirea-progs

max_forks_repo_name: dmitrii-a-piskunov/mirea-progs


text:
mutable struct TAMrates{T}
    binding::MVector{4, T} # fwd/rev binding rate for Ig1, then Ig2
    xRev::MVector{6, T} # xRev 1, 2, 3, 4, 5, 6
    expression::T # AXL expression rate.
    xFwd6::T
end

mutable struct hetRates{T}
    xRev::MVector{10, T} # xRev 7, 8, 9, 10, 11, 12, 13, 14, 15, 16
    xFwd15::T
    xFwd16::T
end

TAMsType{T} = LabelledArrays.@SLVector TAMrates{T} (:Axl, :MerTK, :Tyro3)
hetRType{T} = LabelledArrays.@SLVector hetRates{T} (:AM, :MT, :AT)

mutable struct Rates{T}
    TAMs::TAMsType{T}
    internalize::T # Non-pY species internalization rate.
    pYinternalize::T # pY species internalization rate.
    fElse::T # Recycling fraction for non-D2 species.
    kRec::T # Recycling rate.
    kDeg::T # Degradation rate.
    xFwd::T
    gasCur::T
    hetR::hetRType{T}
end

mutable struct comprates{T}
    TAMs::TAMsType{T}
    diff::T # Receptor diffusivity
    gasPart::T # Partitioning of ligand into PtdSer region
    gasCur::T
    xFwd::T
    hetR::hetRType{T}
end

mutable struct Lsrates{T}
    GBinding::MVector{4, T} # Gas6 fwd/rev rate for Ig1, then Ig2
    PBinding::MVector{4, T} # ProS fwd/rev rate for Ig1, then Ig2
    xRev::MVector{16, T} # Different receptor-receptor unbinding rates

    xFwd27::T
    xFwd29::T

    kRec::T # Recycling rate.
    kDeg::T # Degradation rate.
    fElse::T # Recycling fraction for non-D2 species.
    internalize::T # Non-pY species internalization rate.
    pYinternalize::T # pY species internalization rate.
    expression::T # Receptor expression rate.
    curL::MVector{2, T}
    xFwd::T
end


text:
using Distributed
using Dates
# addprocs(6) # Call instead with: julia -O3 -p18
using RoME
using RoMEPlotting
using Gadfly
@everywhere using RoME, RoMEPlotting, Gadfly

# Make sure plots look a bit nicer.
latex_fonts = Theme(major_label_font="CMU Serif", major_label_font_size=16pt,
                    minor_label_font="CMU Serif", minor_label_font_size=14pt,
                    key_title_font="CMU Serif", key_title_font_size=12pt,
                    key_label_font="CMU Serif", key_label_font_size=10pt)
Gadfly.push_theme(latex_fonts)

# Parse the arguments.
qfl_length = parse(Int, ARGS[1])

# Let's load the Manhattan scenario using the g2o file.
file = (normpath(Base.find_package("RoME"), "../..", "examples", "manhattan_incremental.g2o"))
global instructions = importG2o(file)

# Let's also get the file for batch factor graph at step 500.
fg_name = "fg-after-solve.tar.gz"
global fg_file = (normpath(Base.find_package("RoME"), "../..", "examples", fg_name))

function go_fixedlag_frombatch(qfl_length_arg::Integer)
    # Choose where to save the step's data.
    qfl_length = qfl_length_arg # Fixed lag window size.
    data_logpath = ENV["HOME"]*"/Documents/wafr/manhattan-frombatch-b$(qfl_length)-$(now())"

    # Instead of creating new graph, load the batch one.
    fg = LightDFG{SolverParams}(solverParams=SolverParams(logpath=data_logpath))
    loadDFG(fg_file, Main, fg)
    tree = emptyBayesTree()

    # Set up the fixed lag smoothing.
    getSolverParams(fg).logpath = data_logpath
    getSolverParams(fg).isfixedlag = true
    getSolverParams(fg).qfl = qfl_length
    getSolverParams(fg).limitfixeddown = true
    getSolverParams(fg).dbg = true
    # CHECK I think we need to manually set all variables to not solvable.
    for vsym in ls(fg)
        getVariable(fg, vsym).solvable = 0
    end
    # But release the last 10 poses.
    activevars = []
    for i in 350:360
        push!(activevars, Symbol("x", i))
        getVariable(fg, Symbol("x", i)).solvable = 1
    end
    # CHECK As well as all of the factors.
    for fsym in lsf(fg)
        getFactor(fg, fsym).solvable = 0
    end
    # And activate the factors corresponding to the active variables.
    sym_act_vars = Array{Symbol, 1}(activevars)
    for fsym in getFactorsAmongVariablesOnly(fg, sym_act_vars)
        getFactor(fg, fsym).solvable = 1
    end

    # Add the next---or initial offset of---measurements to the graph.
    padded_step = lpad(501, 4, "0")
    parseG2oInstruction!(fg, instructions[501])

    # Solve the graph, and save a copy of the tree.
    saveDFG(fg, "$(getLogPath(fg))/fg-before-solve$(padded_step)")
    getSolverParams(fg).maxincidence = 1000
    tree = solveTree!(fg)
    saveDFG(fg, "$(getLogPath(fg))/fg-after-solve$(padded_step)")
    saveTree(tree, "$(getLogPath(fg))/tree$(padded_step).jld2")
    drawTree(tree, show=false, filepath="$(getLogPath(fg))/bt$(padded_step).pdf")

    # Analyze clique counts.
    fid = open("$(getLogPath(fg))/clique-counts.txt", "w")
    nCliqs, nMarg, nReused, nBoth = calcCliquesRecycled(tree)
    println(fid, "$(padded_step), $(nCliqs), $(nMarg), $(nReused), $(nBoth)")

    # Just store some quick plots.
    pl1 = drawPoses(fg, spscale=0.6)
    Gadfly.draw(PDF("$(getLogPath(fg))/poses$(padded_step).pdf", 20cm, 10cm), pl1)

    # Solver stride.
    solveStride = 0
    # Run the loop for the remaining time steps.
    for i in 502:5453
        # Add the next measurement to the graph.
        parseG2oInstruction!(fg, instructions[i])
        padded_step = lpad(i, 4, "0")

        # Store each graph.
        saveDFG(fg, "$(getLogPath(fg))/fg-before-solve$(padded_step)")

        # Just store some quick plots, on another process
        remotecall((fgl, padded_stepl) -> begin
          @info "drawPoses, $(padded_stepl), for fg num variables=$(length(ls(fgl)))."
          pl1 = plotSLAM2DPoses(fgl, dyadScale=0.6, lbls=false)
          pl1 |> PDF("$(getLogPath(fgl))/poses$(padded_stepl).pdf", 20cm, 10cm)
        end, rand(Categorical(nprocs()-1))+1, fg, padded_step)

        # Only solve every 10th instruction.
        solveStride += 1
        if solveStride % 10 != 0
            @info "solveStride=$solveStride"
            continue
        end
        @info "Going for solve"

        # Solve the graph, and save a copy of the tree.
        getSolverParams(fg).maxincidence = 1000
        tree = solveTree!(fg, tree)
        saveDFG(fg, "$(getLogPath(fg))/fg-after-solve$(padded_step)")
        saveTree(tree, "$(getLogPath(fg))/tree$(padded_step).jld2")
        drawTree(tree, show=false, filepath="$(getLogPath(fg))/bt$(padded_step).pdf")

        # Analyze clique number.
        nCliqs, nMarg, nReused, nBoth = calcCliquesRecycled(tree)
        println(fid, "$(padded_step), $(nCliqs), $(nMarg), $(nReused), $(nBoth)")
        flush(fid)

        # Run the garbage collector.
        GC.gc()
    end
    # Final plot.
    padded_step = lpad(final_timestep+1, 4, "0")
    pl1 = drawPoses(fg, spscale=0.6)
    pl1 |> PDF("$(getLogPath(fg))/poses$(padded_step).pdf", 20cm, 10cm)

    close(fid)
end

# Run within a function to avoid undefined variable errors.
go_fixedlag_frombatch(qfl_length)


text:
@testset "Lagrangian Dynamics" begin
	# Code gen
	dir = joinpath(module_dir(), "src", "dynamics", "rigidbody")
	# include(joinpath(dir, "model.jl"))

	# Analytical Model
	model = deepcopy(ContactImplicitMPC.rigidbody)

	path_base = joinpath(dir, "test_dynamics/base.jld2")
	path_dyn = joinpath(dir, "test_dynamics/dynamics.jld2")

	expr_base = generate_base_expressions(model,
		M_analytical = true,
		mapping = ContactImplicitMPC.G_func,
		nv = model.dim.q - 1)

	save_expressions(expr_base, path_base, overwrite=true)
	instantiate_base!(model, path_base)

	expr_dyn = generate_dynamics_expressions(model, nv = model.dim.q - 1)
	save_expressions(expr_dyn, path_dyn, overwrite=true)
	instantiate_dynamics!(model, path_dyn)

	# Lagrangian Model
	model_lag = deepcopy(ContactImplicitMPC.rigidbody)
	path_base_lag = joinpath(dir, "test_dynamics/base_lag.jld2")
	path_dyn_lag = joinpath(dir, "test_dynamics/dynamics_lag.jld2")

	expr_base_lag = generate_base_expressions(model_lag,
		M_analytical = false,
		mapping = ContactImplicitMPC.G_func,
		nv = model.dim.q - 1)

	save_expressions(expr_base, path_base_lag, overwrite=true)
	instantiate_base!(model_lag, path_base_lag)

	expr_dyn_lag = generate_dynamics_expressions(model_lag, nv = model.dim.q - 1)
	save_expressions(expr_dyn_lag, path_dyn_lag, overwrite=true)
	instantiate_dynamics!(model_lag, path_dyn_lag)

	# Data
	env = environment_3D_flat()
	nq = model.dim.q
	nv = model.dim.q - 1
	nu = model.dim.u
	nw = model.dim.w
	nc = model.dim.c
	h = 0.01

	# initial conditions
	r0 = [0.0; 0.0; 1.0]
	v0 = [7.5; 5.0; 0.0]

	quat0 = [1.0; 0.0; 0.0; 0.0]
	ω0 = [0.0; 0.0; 0.0]

	q0 = SVector{model.dim.q}([r0; quat0])
	q1 = SVector{model.dim.q}([r0 + v0 * h; 0.5 * h * L_multiply(quat0) * [sqrt((2.0 / h)^2.0 - ω0' * ω0); ω0]])
	q2 = SVector{model.dim.q}([r0 + v0 * 2h; 0.5 * 2h * L_multiply(quat0) * [sqrt((2.0 / (2h))^2.0 - ω0' * ω0); ω0]])
	q̇0 = zeros(SVector{nv})
	u1 = zeros(SVector{nu})
	w1 = zeros(SVector{nw})
	λ1 = zeros(SVector{nc * dim(env)})

	# Test equality
	@test norm(ContactImplicitMPC.M_fast(model, q0) -
		       ContactImplicitMPC.M_fast(model_lag, q0)) < 1e-10
	@test norm(ContactImplicitMPC.C_fast(model, q0, q̇0) -
			   ContactImplicitMPC.C_fast(model_lag, q0, q̇0)) < 1e-10
	# @test norm(ContactImplicitMPC.d_fast(model, [h], q0, q1, u1, w1, λ1, q2) -
	# 		   ContactImplicitMPC.d_fast(model_lag, [h], q0, q1, u1, w1, λ1, q2)) < 1e-10 # Needs to be fixed
end


text:
include("../src/TableSchema.jl")
using TableSchema

# Either import the functions, or use TableSchema.read(<Table>) in the code
import TableSchema: read, is_valid, validate

import DelimitedFiles: readdlm

t = Table()

s = Schema("../data/schema_invalid_empty.json")
if is_valid(s) == false; println("An invalid Schema was found"); end

s = Schema("../data/schema_valid_missing.json")
if is_valid(s); println("A valid Schema is ready"); end

t.schema = s
source = readdlm("../data/data_types.csv", ',')
tr = read(t, data=source) # 5x5 Array{Any,2}
println( "The length is ", length(tr[:,1]) ) # 5
println( "Sum of column 2 is ", sum([ row for row in tr[2] ]) ) # 51.0

if validate(t); println("The table is valid according to the Schema"); end

t2 = Table("../data/data_constraints.csv", s)
if validate(t2) == false; println("This other table is not valid"); end
for err in t2.errors
    println(string(err.field.name, " has (expected) error: ", err.message))
end

max_stars_repo_name: loleg/tableschema-jl

max_issues_repo_name: loleg/tableschema-jl

max_forks_repo_name: loleg/tableschema-jl


text:
using Pkg;
Pkg.activate("");
using CUTEst

nmax = 1000000 # maximum is 1000000
_pnames = CUTEst.select(contype = "unc", max_var = nmax)

#Remove all the problems ending by NE as Ipopt cannot handle them.
pnamesNE = _pnames[findall(x -> occursin(r"NE\b", x), _pnames)]
pnames = setdiff(_pnames, pnamesNE)

open("list_problems_$nmax.dat", "w") do io
    for name in pnames
        write(io, name * "\n")
    end
end


text:
@testset "$TEST $G" begin

f, fio = mktemp()

g = G(10,0)
writegraph(f, g, :gexf)
ga = readgraph(f, :gexf, G)
@test g == ga

g = G(10,20)
writegraph(f, g, :gexf)
ga = readgraph(f, :gexf, G)
@test g == ga

g = DG(10,20)
writegraph(f, g, :gexf)
ga = readgraph(f, :gexf, G)
@test g == ga

if G <: ANetwork ##################

g = readnetwork(joinpath(testdir, "testdata", "lesmis.gexf"), G)
writenetwork(f, g, :gexf)
h = readnetwork(f, :gexf, G)
test_networks_eq(g, h)

for gname in [:lesmis, :karate]
    g = readnetwork(gname, G)
    rem_gprop!(g, "readme") # TODO should strip characters " " in the readme string
    rem_gprop!(g, "description")
    rem_vprop!(g, "pos") # TODO vector attributes not supported
    writenetwork(f, g, :gexf)
    h = readnetwork(f, :gexf, G)
    test_networks_eq(g, h)
end

end #if network ####################

isfile(f) && rm(f)
 end #testset


text:
using MLDatasets
using Flux
using Flux: train!, @epochs, mse, throttle, params
using ConditionalDists
using ConditionalDists: SplitLayer
using GenerativeModels

train_x, _ = MNIST.traindata()
train_x = reshape(train_x, :, size(train_x, 3)) |> gpu
data = Flux.Data.DataLoader(train_x, batchsize=200, shuffle=true)

encoder = Chain(
  Dense(size(flat_x, 1), 512, relu),
  Dense(512, 256, relu),
  SplitLayer(256, [2, 2], [identity, softplus])
) |> ConditionalMvNormal
decoder = Chain(
  Dense(2, 256, relu),
  Dense(256, 512, relu),
  SplitLayer(512, [size(flat_x, 1), 1], σ)
) |> ConditionalMvNormal
model = VAE(zlength, encoder, decoder) |> gpu

loss(x) = -elbo(model,x)
ps = params(model)
opt = ADAM()

@epochs 50 begin
  @info "Epoch $e" loss(flat_x)
  Flux.train!(loss, ps, data, opt)
end


text:
"""
```
eqcond(m::Schorf)
```

Expresses the equilibrium conditions in canonical form using Γ0, Γ1, C, Ψ, and Π matrices.
Using the mappings of states/equations to integers defined in schorf.jl, coefficients are
specified in their proper positions.

### Outputs

* `Γ0` (`n_states` x `n_states`) holds coefficients of current time states.
* `Γ1` (`n_states` x `n_states`) holds coefficients of lagged states.
* `C`  (`n_states` x `1`) is a vector of constants
* `Ψ`  (`n_states` x `n_shocks_exogenous`) holds coefficients of iid shocks.
* `Π`  (`n_states` x `n_states_expectational`) holds coefficients of expectational states.
"""
function eqcond(m::AnSchorfheide)
    endo = m.endogenous_states
    exo  = m.exogenous_shocks
    ex   = m.expected_shocks
    eq   = m.equilibrium_conditions

    Γ0 = zeros(n_states(m), n_states(m))
    Γ1 = zeros(n_states(m), n_states(m))
    C  = zeros(n_states(m))
    Ψ  = zeros(n_states(m), n_shocks_exogenous(m))
    Π  = zeros(n_states(m), n_shocks_expectational(m))

    ### ENDOGENOUS STATES ###

    ### 1. Consumption Euler Equation

    Γ0[eq[:eq_euler], endo[:y_t]]  = 1.
    Γ0[eq[:eq_euler], endo[:R_t]] = 1/m[:τ]
    Γ0[eq[:eq_euler], endo[:g_t]] = -(1-m[:ρ_g])
    Γ0[eq[:eq_euler], endo[:z_t]] = -m[:ρ_z]/m[:τ]
    Γ0[eq[:eq_euler], endo[:Ey_t]] = -1
    Γ0[eq[:eq_euler], endo[:Eπ_t]] = -1/m[:τ]

    ### 2. NK Phillips Curve

    Γ0[eq[:eq_phillips], endo[:y_t]] = -m[:κ]
    Γ0[eq[:eq_phillips], endo[:π_t]] = 1
    Γ0[eq[:eq_phillips], endo[:g_t]] = m[:κ]
    Γ0[eq[:eq_phillips], endo[:Eπ_t]] = -1/(1+m[:rA]/400)

    ### 3. Monetary Policy Rule

    Γ0[eq[:eq_mp], endo[:y_t]] = -(1-m[:ρ_R])*m[:ψ_2]
    Γ0[eq[:eq_mp], endo[:π_t]] = -(1-m[:ρ_R])*m[:ψ_1]
    Γ0[eq[:eq_mp], endo[:R_t]] = 1
    Γ0[eq[:eq_mp], endo[:g_t]] = (1-m[:ρ_R])*m[:ψ_2]
    Γ1[eq[:eq_mp], endo[:R_t]] = m[:ρ_R]
    Ψ[eq[:eq_mp], exo[:rm_sh]] = 1

    ### 4. Output lag

    Γ0[eq[:eq_y_t1], endo[:y_t1]] = 1
    Γ1[eq[:eq_y_t1], endo[:y_t]] = 1

    ### 5. Government spending

    Γ0[eq[:eq_g], endo[:g_t]] = 1
    Γ1[eq[:eq_g], endo[:g_t]] = m[:ρ_g]
    Ψ[eq[:eq_g], exo[:g_sh]] = 1

    ### 6. Technology

    Γ0[eq[:eq_z], endo[:z_t]] = 1
    Γ1[eq[:eq_z], endo[:z_t]] = m[:ρ_z]
    Ψ[eq[:eq_z], exo[:z_sh]] = 1

    ### 7. Expected output

    Γ0[eq[:eq_Ey], endo[:y_t]] = 1
    Γ1[eq[:eq_Ey], endo[:Ey_t]] = 1
    Π[eq[:eq_Ey], ex[:Ey_sh]] = 1

    ### 8. Expected inflation

    Γ0[eq[:eq_Eπ], endo[:π_t]] = 1
    Γ1[eq[:eq_Eπ], endo[:Eπ_t]] = 1
    Π[eq[:eq_Eπ], ex[:Eπ_sh]] = 1

    return Γ0, Γ1, C, Ψ, Π
end


text:
"""
    upgma(dm::Array{Float64,2}, Array{String,1})

This function returns a phylogenetic tree by using UPGMA based on a
given distance matrix and an array of leaf names.

Returns a node of the resulting tree, from which it can be traversed.

* `dm` : Matrix from which to create the tree.

* `leaf_names` : array of strings containing names of leaf nodes.
"""
function upgma(dm::Array{Float64,2}, leaf_names::Array{String,1})
    n = size(dm)[1]
    if n != size(leaf_names, 1)
        throw("Distance Matrix and leaf Names size do not match")
    end # end if
    # build array of leaves from leaf names
    leaves = Array{Tuple{Node,Float64,Int64},1}(undef, size(dm, 1))
    for (ind, leaf) in enumerate(leaf_names)
        new_leaf = Node(leaf)
        leaves[ind] = (new_leaf, 0.0, 1)
    end # end for
    upgma_int(dm, leaves)
end

"""
    upgma(dm::Array{Float64,2})

This function returns a phylogenetic tree by using UPGMA based on a
given distance matrix. Creates an array of nodes to be used as leaves.

Returns a node of the resulting tree, from which it can be traversed.

* `dm` : Matrix from which to create the tree.
"""
function upgma(dm::Array{Float64,2})
    n = size(dm)[1]
    leaves = Array{Tuple{Node,Float64,Int64},1}(undef, n)
    # build array of dummy leaves
    for i in 1:n
        new_leaf = Node("leaf_$i")
        leaves[n] = (new_leaf, 0.0, 1)
    end # end for
    upgma_int(dm, leaves)
end

"""
    upgma_int(dm::Array{Float64,2},leaves::Vector{Node})

--- INTERNAL ---
Internal function that is called by both UPGMA methods. Contains the
actual UPGMA algorithm, i.e. builds a phylogenetic tree from the
given distance matrix and array of leaves. Returns a node of that tree, from
which it can be traversed.
"""
function upgma_int(
    dm::Array{Float64,2},
    clusters::Array{Tuple{Node,Float64,Int64},1},
)

    n = size(dm)[1]
    # count for node names
    count = 1
    while true
        # set 0 entries to infinity, to enable use of findmin method
        dm[CartesianIndex.(1:n, 1:n)] .= Inf
        index = findmin(dm)[2]
        first_cluster = clusters[index[2]]
        second_cluster = clusters[index[1]]
        # create new Nodecluster
        new_cluster = Node("Cluster_$count")
        count += 1
        # calculate length of new branches
        total_path_length = dm[index] / 2
        first_cluster[1].inc_length = total_path_length - first_cluster[2]
        second_cluster[1].inc_length = total_path_length - second_cluster[2]
        # add children to the new Nodecluster
        add_child!(new_cluster, first_cluster[1])
        add_child!(new_cluster, second_cluster[1])
        # return (root) node when algorithm is finished
        n -= 1
        if n == 1
            set_binary!(new_node)
            number_nodes!(new_node)
            return new_cluster
        end
        # initalize cluster weights
        first_cluster_weight = clusters[index[2]][3]
        second_cluster_weight = clusters[index[1]][3]
        # initalize next distance matrix
        next_dm = zeros(Float64, n, n)
        # fill first row and column of next distance matrix
        j = 1
        for i in 2:n
            while j == index[1] || j == index[2]
                j += 1
            end # end while
            next_dm[i, 1] =
                next_dm[1, i] =
                    (
                        dm[index[2], j] * first_cluster_weight +
                        dm[index[1], j] * second_cluster_weight
                    ) / (first_cluster_weight + second_cluster_weight)
            j += 1
        end # end for
        # copy values of the last distance matrix to fill
        # the remaining rows and columns of the new one
        array = [1:1:n+1;]
        deleteat!(array, [index[2], index[1]])
        next_dm[2:end, 2:end] .= dm[array, array]
        dm = next_dm
        # calculate weight of the next cluster and add it to the leaves list
        new_cluster_weight = first_cluster_weight + second_cluster_weight
        # update array with leaves
        deleteat!(clusters, [index[2], index[1]])
        insert!(
            clusters,
            1,
            (new_cluster, total_path_length, new_cluster_weight),
        )
    end # end while
end # end function upgma

"""
    neighbor_joining(dm::Array{Float64,2}, Array{String,1})

This function returns a phylogenetic tree by using neighbor-joining based on a
given distance matrix and an array of leaf names.

Returns a node of the resulting tree, from which it can be traversed.

* `dm` : Matrix used to create Tree.

* `leaf_names` : Array containing names of leaf nodes.
"""
function neighbor_joining(dm::Array{Float64,2}, leaf_names::Array{String,1})
    n = size(dm)[1]
    if n != size(leaf_names, 1)
        throw("Distance Matrix and leaf names array size do not match")
    end # end if
    # build array of leaves from leaf names
    leaves = Array{Node,1}(undef, size(dm, 1))
    for (ind, leaf) in enumerate(leaf_names)
        new_leaf = Node(leaf)
        leaves[ind] = new_leaf
    end # end for
    neighbor_joining_int(dm, leaves)
end

"""
    neighbor_joining(dm::Array{Float64,2})

This function returns a phylogenetic tree by using neighbor-joining based on a
given distance matrix. Creates an array of nodes to be used as leaves.

Returns a node of the resulting tree, from which it can be traversed.

* `dm` : Matrix from which to create tree.
"""
function neighbor_joining(dm::Array{Float64,2})
    n = size(dm)[1]
    leaves = Array{Node,1}(undef, n)
    # build array of dummy leaves
    for i in 1:n
        new_leaf = Node("leaf_$i")
        leaves[i] = new_leaf
    end # end for
    neighbor_joining_int(dm, leaves)
end

"""
    neighbor_joining_int(dm::Array{Float64,2},leaves::Vector{Node})

--- INTERNAL ---
Internal function that is called by both neighbor_joining methods. Contains the
actual neighbor-joining algorithm, i.e. builds a phylogenetic tree from the
given distance matrix and array of leaves. Returns a node of that tree, from
which it can be traversed.
"""
function neighbor_joining_int(dm::Array{Float64,2}, leaves::Vector{Node})

    n = size(dm)[1]
    # count for node names
    count = 1
    while n > 1
        # build Q1 matrix
        q1 = zeros(Float64, n, n)
        for i in 1:n
            for j in i+1:n
                q1[i, j] =
                    q1[j, i] =
                        (n - 2) * dm[i, j] - sum(dm[i, :]) - sum(dm[j, :])
            end # end for
        end # end for
        # locate minimum of Q1
        index = findmin(q1)[2]
        first_node = leaves[index[2]]
        second_node = leaves[index[1]]
        new_node = Node("Node_$count")
        count += 1
        first_node.inc_length =
            0.5 * dm[index] +
            (1 / (2 * (n - 2))) * (sum(dm[index[2], :]) - sum(dm[index[1], :]))
        second_node.inc_length = dm[index] - first_node.inc_length
        add_child!(new_node, first_node)
        add_child!(new_node, second_node)
        # update array with leaves
        deleteat!(leaves, [index[2], index[1]])
        insert!(leaves, 1, new_node)
        # initalize next distance matrix
        n -= 1
        next_dm = zeros(Float64, n, n)
        # fill first row and column of next distance matrix
        j = 1
        for i in 2:n
            while j == index[1] || j == index[2]
                j += 1
            end # end while
            next_dm[i, 1] =
                next_dm[1, i] =
                    0.5 *
                    (dm[index[1], j] + dm[index[2], j] - dm[index[1], index[2]])
            j += 1
        end # end for
        # copy values of last distance matrix to finish filling the next one
        array = [1:1:n+1;]
        deleteat!(array, [index[2], index[1]])
        next_dm[2:end, 2:end] .= dm[array, array]
        # add final node to tree, if only 2 nodes are left in the list of nodes
        dm = next_dm
        if n == 2
            final_leaf = pop!(leaves)
            final_leaf.inc_length = dm[1, 2]
            # add third child to new node created previously
            add_child!(new_node, final_leaf)
            set_binary!(new_node)
            number_nodes!(new_node)
            return new_node
        end # end if
    end # end while
end # end function neighbor_joining


text:
import Base.isempty

export AbstractPolytope,
       vertices_list,
       singleton_list,
       isempty,
       minkowski_sum

"""
    AbstractPolytope{N<:Real} <: AbstractPolyhedron{N}

Abstract type for compact convex polytopic sets.

### Notes

Every concrete `AbstractPolytope` must define the following functions:
- `vertices_list(::AbstractPolytope{N})::Vector{Vector{N}}` -- return a list of
    all vertices

```jldoctest; setup = :(using LazySets: subtypes)
julia> subtypes(AbstractPolytope)
4-element Array{Any,1}:
 AbstractCentrallySymmetricPolytope
 AbstractPolygon
 HPolytope
 VPolytope
```

A polytope is a bounded polyhedron (see [`AbstractPolyhedron`](@ref)).
Polytopes are compact convex sets with either of the following equivalent
properties:
1. They are the intersection of a finite number of closed half-spaces.
2. They are the convex hull of finitely many vertices.
"""
abstract type AbstractPolytope{N<:Real} <: AbstractPolyhedron{N} end

# =============================================
# Common AbstractPolytope functions
# =============================================

"""
    isbounded(P::AbstractPolytope)::Bool

Determine whether a polytopic set is bounded.

### Input

- `P` -- polytopic set

### Output

`true` (since a polytope must be bounded).
"""
function isbounded(::AbstractPolytope)::Bool
    return true
end

"""
    singleton_list(P::AbstractPolytope{N})::Vector{Singleton{N}} where {N<:Real}

Return the vertices of a polytopic set as a list of singletons.

### Input

- `P` -- polytopic set

### Output

List containing a singleton for each vertex.
"""
function singleton_list(P::AbstractPolytope{N}
                       )::Vector{Singleton{N}} where {N<:Real}
    return [Singleton(vi) for vi in vertices_list(P)]
end

"""
    isempty(P::AbstractPolytope)::Bool

Determine whether a polytope is empty.

### Input

- `P` -- abstract polytope

### Output

`true` if the given polytope contains no vertices, and `false` otherwise.

### Algorithm

This algorithm checks whether the `vertices_list` of the given polytope is empty
or not.
"""
function isempty(P::AbstractPolytope)::Bool
    return isempty(vertices_list(P))
end

function default_polyhedra_backend(P, N)
    require(:Polyhedra; fun_name="default_polyhedra_backend")
    error("no default backend for numeric type $N")
end

# given a polytope P, apply the linear map P to each vertex of P
# it is assumed that the interface function `vertices_list(P)` is available
@inline function _linear_map_vrep(M::AbstractMatrix{N}, P::AbstractPolytope{N}) where {N<:Real}
    vertices = broadcast(v -> M * v, vertices_list(P))
    m = size(M, 1) # output dimension
    if m == 1
        # TODO: substitute with Interval(convex_hull(vertices)...): convex hull 1D
        return Interval(minimum(vertices)[1], maximum(vertices)[1])
    elseif m == 2
        return VPolygon(vertices)
    else
        return VPolytope(vertices)
    end
end

@inline function _linear_map_hrep(M::AbstractMatrix{N}, P::AbstractPolytope{N},
                          use_inv::Bool) where {N<:Real}
    constraints = _linear_map_hrep_helper(M, P, use_inv)
    m = size(M, 1) # output dimension
    if m == 1
        return convert(Interval, HPolygon(constraints))
    elseif m == 2
        return HPolygon(constraints)
    else
        return HPolytope(constraints)
    end
end

"""
    minkowski_sum(P1::AbstractPolytope{N}, P2::AbstractPolytope{N};
                  [apply_convex_hull]=true,
                  [backend]=default_polyhedra_backend(P1, N),
                  [solver]=default_lp_solver(N)) where {N<:Real}

Compute the Minkowski sum between two polytopes using their vertex representation.

### Input

- `P1`                -- polytope
- `P2`                -- another polytope
- `apply_convex_hull` -- (optional, default: `true`) if `true`, post-process the
                         pairwise sums using a convex hull algorithm
- `backend`           -- (optional, default: `default_polyhedra_backend(P1, N)`)
                         the backend for polyhedral computations used to
                         post-process with a convex hull
- `solver`            -- (optional, default: `default_lp_solver(N)`) the linear programming
                         solver used in the backend

### Output

A new polytope in vertex representation whose vertices are the convex hull of
the sum of all possible sums of vertices of `P1` and `P2`.
"""
function minkowski_sum(P1::AbstractPolytope{N}, P2::AbstractPolytope{N};
                       apply_convex_hull::Bool=true,
                       backend=nothing,
                       solver=nothing) where {N<:Real}

    @assert dim(P1) == dim(P2) "cannot compute the Minkowski sum between a polyotope " *
        "of dimension $(dim(P1)) and a polytope of dimension $((dim(P2)))"

    vlist1 = _vertices_list(P1, backend)
    vlist2 = _vertices_list(P2, backend)
    n, m = length(vlist1), length(vlist2)
    Vout = Vector{Vector{N}}()
    sizehint!(Vout, n * m)
    for vi in vlist1
        for vj in vlist2
            push!(Vout, vi + vj)
        end
    end
    if apply_convex_hull
        if backend == nothing
            backend = default_polyhedra_backend(P1, N)
            solver = default_lp_solver(N)
        end
        convex_hull!(Vout, backend=backend, solver=solver)
    end
    return VPolytope(Vout)
end

# the "backend" argument is ignored, used for dispatch
function _vertices_list(P::AbstractPolytope, backend)
    return vertices_list(P)
end

# =============================================
# Functions that require Polyhedra
# =============================================

function load_polyhedra_abstractpolytope() # function to be loaded by Requires
return quote
import .Polyhedra: polyhedron
export polyhedron
using .Polyhedra: HRep, VRep,
                  removehredundancy!, removevredundancy!,
                  hreps, vreps,
                  intersect,
                  convexhull,
                  hcartesianproduct, vcartesianproduct,
                  points,
                  default_library
import JuMP, GLPK

function default_polyhedra_backend(P, N::Type{<:AbstractFloat})
    return default_library(LazySets.dim(P), Float64)
end

function default_polyhedra_backend(P, N::Type{<:Rational})
    return default_library(LazySets.dim(P), Rational{Int})
end

function default_lp_solver(N::Type{<:AbstractFloat})
    return JuMP.with_optimizer(GLPK.Optimizer)
end

function default_lp_solver(N::Type{<:Rational})
    return JuMP.with_optimizer(GLPK.Optimizer, method=:Exact)
end

end # quote
end # function load_polyhedra_hpolytope()


text:
function set_wally(nodes, bdcon, wally, cellcenter, cellxmax, cellymax, cellzmax, icell)
	# serch wall
	swith_wall = zeros(6)   # x-, x+, y-, y+, z-, z+を検索する
							# 1:壁
	for i in 1:6
		if Int(bdcon[i][1]) == 2 || Int(bdcon[i][1]) == 7
			swith_wall[i] = 1
		end
	end

	# cal number of ponits
	nop = 0
	if swith_wall[1] == 1
		nop = nop + (cellymax-1)*(cellzmax-1)
	end
	if swith_wall[2] == 1
		nop = nop + (cellymax-1)*(cellzmax-1)
	end
	if swith_wall[3] == 1
		nop = nop + (cellzmax-1)*(cellxmax-1)
	end
	if swith_wall[4] == 1
		nop = nop + (cellzmax-1)*(cellxmax-1)
	end
	if swith_wall[5] == 1
		nop = nop + (cellxmax-1)*(cellzmax-1)
	end
	if swith_wall[6] == 1
		nop = nop + (cellxmax-1)*(cellzmax-1)
	end

	wallpoint = zeros(nop, 3) # (格子数, 座標)
	ite = 1

	# x-
	if swith_wall[1] == 1
		for k in 1+icell:cellzmax-icell
			for j in 1+icell:cellymax-icell
				wallpoint[ite,1] = nodes[icell+1,j,k,1]
				wallpoint[ite,2] = nodes[icell+1,j,k,2]
				wallpoint[ite,3] = nodes[icell+1,j,k,3]
				ite = ite + 1
			end
		end
	end
	# x+
	if swith_wall[2] == 1
		for k in 1+icell:cellzmax-icell
			for j in 1+icell:cellymax-icell
				wallpoint[ite,1] = nodes[cellxmax-icell+1,j,k,1]
				wallpoint[ite,2] = nodes[cellxmax-icell+1,j,k,2]
				wallpoint[ite,3] = nodes[cellxmax-icell+1,j,k,3]
				ite = ite + 1
			end
		end
	end
	# y-
	if swith_wall[3] == 1
		for k in 1+icell:cellzmax-icell
			for i in 1+icell:cellxmax-icell
				wallpoint[ite,1] = nodes[i,icell+1,k,1]
				wallpoint[ite,2] = nodes[i,icell+1,k,2]
				wallpoint[ite,3] = nodes[i,icell+1,k,3]
				ite = ite + 1
			end
		end
	end
	# y+
	if swith_wall[4] == 1
		for k in 1+icell:cellzmax-icell
			for i in 1+icell:cellxmax-icell
				wallpoint[ite,1] = nodes[i,cellymax-icell+1,k,1]
				wallpoint[ite,2] = nodes[i,cellymax-icell+1,k,2]
				wallpoint[ite,3] = nodes[i,cellymax-icell+1,k,3]
				ite = ite + 1
			end
		end
	end
	# z-
	if swith_wall[5] == 1
		for j in 1+icell:cellymax-icell
			for i in 1+icell:cellxmax-icell
				wallpoint[ite,1] = nodes[i,j,icell+1,1]
				wallpoint[ite,2] = nodes[i,j,icell+1,2]
				wallpoint[ite,3] = nodes[i,j,icell+1,3]
				ite = ite + 1
			end
		end
	end
	# z+
	if swith_wall[6] == 1
		for j in 1+icell:cellymax-icell
			for i in 1+icell:cellxmax-icell
				wallpoint[ite,1] = nodes[i,j,cellzmax-icell+1,1]
				wallpoint[ite,2] = nodes[i,j,cellzmax-icell+1,2]
				wallpoint[ite,3] = nodes[i,j,cellzmax-icell+1,3]
				ite = ite + 1
			end
		end
	end

	# 距離格納用
	distance_list = zeros(nop, 3)
	distance = 0

	# 最大値を2つ抽出
	for k in 1+icell:cellzmax-icell
		for j in 1+icell:cellymax-icell
			for i in 1+icell:cellxmax-icell
				for n in 1:nop
					distance = ((cellcenter[i,j,k,1]-wallpoint[n,1])^2 +
								(cellcenter[i,j,k,2]-wallpoint[n,2])^2 +
								(cellcenter[i,j,k,3]-wallpoint[n,3])^2
								)^0.5
					distance_list[n,1] = n
					distance_list[n,2] = distance
				end

				qdistance_list = copy(distance_list)
				qdistance_list = quicksort(qdistance_list,1,nop,2,2)

				# wall_point
				point1 = Int(qdistance_list[1,1])
				point2 = Int(qdistance_list[2,1])

				dis_bottom = ((wallpoint[point1, 1] - wallpoint[point2, 1])^2 +
							  (wallpoint[point1, 2] - wallpoint[point2, 2])^2 +
							  (wallpoint[point1, 3] - wallpoint[point2, 3])^2)^0.5

				# O型格子の関係で同じ点が二つあるため
				if dis_bottom == 0.0
					point1 = Int(qdistance_list[1,1])
					point2 = Int(qdistance_list[3,1])

					dis_bottom = ((wallpoint[point1, 1] - wallpoint[point2, 1])^2 +
				   			      (wallpoint[point1, 2] - wallpoint[point2, 2])^2 +
							      (wallpoint[point1, 3] - wallpoint[point2, 3])^2)^0.5
				end

				dis1 = distance_list[point1, 2]
				dis2 = distance_list[point2, 2]

				# ヘロンの公式から距離を算出
				# 正確な平面への投影はしてないので、注意
				s = (dis_bottom + dis1 + dis2)*0.5
				surface = (s*(s-dis_bottom)*(s-dis1)*(s-dis2))^0.5
				wally[i,j,k] = surface/dis_bottom * 2
			end
		end
	end

	# 仮想セルへの代入
	for k in 1+icell:cellzmax-icell
		for j in 1:icell
			for i in 1+icell:cellxmax-icell
				wally[i,j,k] = wally[i,icell+1,k]
				wally[i,cellymax-j+1,k] = wally[i,cellymax-icell,k]
			end
		end
	end
	for k in 1+icell:cellzmax-icell
		for j in 1+icell:cellymax-icell
			for i in 1:icell
				wally[i,j,k] = wally[icell+1,j,k]
				wally[cellxmax-i+1,j,k] = wally[cellxmax-icell,j,k]
			end
		end
	end
	for k in 1:icell
		for j in 1+icell:cellymax-icell
			for i in 1+icell:cellxmax-icell
				wally[i,j,k] = wally[i,j,icell+1]
				wally[i,j,cellzmax-k+1] = wally[i,j,cellzmax-icell]
			end
		end
	end

	return wally, swith_wall
end

function cal_yplus(yplus, Qbase, wally, swith_wall, mu, cellxmax, cellymax, cellzmax, vecAx, vecAy, vecAz, volume, icell)
	tvA = zeros(3)

	# 壁に並行な方向の速度を求める
	# 速度ベクトル内壁に垂直な成分は、ux = vec{u}*vecAx
	# u_purpose = (u^2+v^2+w^2 - ux^2)^2
	for k in 1+icell:cellzmax-icell
		for j in 1+icell:cellymax-icell
			for i in 1+icell:cellzmax-icell
				if swith_wall[1] == 1 || swith_wall[2] == 1
					tvA[1] = 0.5*(vecAx[i-1,j,k,1] + vecAx[i,j,k,1])
					tvA[2] = 0.5*(vecAx[i-1,j,k,2] + vecAx[i,j,k,2])
					tvA[3] = 0.5*(vecAx[i-1,j,k,3] + vecAx[i,j,k,3])
				elseif swith_wall[3] == 1 || swith_wall[4] == 1
					tvA[1] = 0.5*(vecAy[i,j-1,k,1] + vecAy[i,j,k,1])
					tvA[2] = 0.5*(vecAy[i,j-1,k,2] + vecAy[i,j,k,2])
					tvA[3] = 0.5*(vecAy[i,j-1,k,3] + vecAy[i,j,k,3])
				elseif swith_wall[5] == 1 || swith_wall[6] == 1
					tvA[1] = 0.5*(vecAz[i,j,k-1,1] + vecAz[i,j,k,1])
					tvA[2] = 0.5*(vecAz[i,j,k-1,2] + vecAz[i,j,k,2])
					tvA[3] = 0.5*(vecAz[i,j,k-1,3] + vecAz[i,j,k,3])
				end

				sqA = (tvA[1]^2 + tvA[2]^2 + tvA[3]^2)^0.5
				tvA[1] = tvA[1] / sqA
				tvA[2] = tvA[2] / sqA
				tvA[3] = tvA[3] / sqA

				nu = mu[i,j,k] / Qbase[i,j,k,1]
				u = tvA[1]*Qbase[i,j,k,2] + tvA[2]*Qbase[i,j,k,3] + tvA[3]*Qbase[i,j,k,4]
				ubar = (Qbase[i,j,k,2]^2 + Qbase[i,j,k,3]^2 + Qbase[i,j,k,4]^2 - u^2)^0.5
				yplus[i,j,k], up = cal_wall_val_spalding(ubar, wally[i,j,k], nu)
			end
		end
	end

	return yplus
end

function cal_wall_val_spalding(u,y,nu)
	# spalding 定数
	k = 0.4
	B = 5.5

	# ニュートン法
	up = 200.0
	for i in 1:100
		old = up

		up = up - spalding(up,u,y,nu,k,B)/spalding_dash(up,u,y,nu,k,B)

		if abs(up-old)/up < 10^(-4)
			break
		end
	end

	yplus = u*y/nu/up
	return yplus, up
end

function spalding(up,u,y,nu,k,B)
	t = k*up
	F = - 1/up* u*y/nu + up + exp(-k*B)*(exp(t)-1-t-t^2/2-t^3/6)
	return F
end

function spalding_dash(up,u,y,nu,k,B)
	t = k*up
	Fdash = 1/up^2* u*y/nu + 1 + exp(-k*B)*(k*exp(t)-k-k^2*up-k^3*up^2/2)
	return Fdash
end

# test
#=
uu = 6.64037328
yy = 0.06299976860
nunu = 1.406107133e-5
println(cal_wall_val_spalding(uu,yy,nunu))
throw(UndefVarError(:x))
=#

function wallf_Van_Driest(yplus)
	f = 1 - exp(-yplus/26)
	return f
end


text:
export pixel_shuffle

"""
    pixel_shuffle(x, r)

Pixel shuffling operation. `r` is the upscale factor for shuffling.
The operation converts an input of size [W,H,r²C,N] to size [rW,rH,C,N]
Used extensively in super-resolution networks to upsample
towards high resolution features.

Reference : https://arxiv.org/pdf/1609.05158.pdf
"""
function pixel_shuffle(x::AbstractArray, r::Integer)
    @assert ndims(x) > 2
    d = ndims(x) - 2
    sizein = size(x)[1:d]
    cin, n = size(x, d+1), size(x, d+2)
    @assert cin % r^d == 0
    cout = cin ÷ r^d
    # x = reshape(x, sizein..., fill(r, d)..., cout, n) # bug https://github.com/FluxML/Zygote.jl/issues/866
    x = reshape(x, sizein..., ntuple(i->r, d)..., cout, n)
    perm = [d+1:2d 1:d]' |> vec  # = [d+1, 1, d+2, 2, ..., 2d, d]
    x = permutedims(x, (perm..., 2d+1, 2d+2))
    return reshape(x, ((r .* sizein)..., cout, n))
end


text:
# This file was generated, do not modify it. # hide
x = randn(1_000) # 1_000 points iid from a N(0, 1)
μ = mean(x)
σ = std(x)
@show (μ, σ)


text:
@testset "Continous Power Analysis" begin
    @testset "Bernoulli model with probability to beat all with same distribution" begin
        Random.seed!(1234)
        truedists = [Bernoulli(0.2), Bernoulli(0.2)]

        modelA = ConjugateBernoulli(1, 1)
        modelB = ConjugateBernoulli(1, 1)
        stoppingrule = ProbabilityBeatAllThresh(0.99)
        experiment = ExperimentAB([modelA, modelB], stoppingrule)

        simulation = Simulation(
            experiment=experiment,
            datagendists=truedists,
            maxsteps=30,
            onestepsizes=[10000, 10000],
            minsteps=5
        )

        numsteps, winners, _ = runsequential(simulation, numsamples=1000, numsims=50)
        @test sum(winners .== "control") == 3
        @test sum(winners .== "variant 1") ==  2
        @test sum(winners .== "nothing") == 45
    end

    @testset "Bernoulli model with expected loss with same distribution" begin
        Random.seed!(1234)
        truedists = [Bernoulli(0.2), Bernoulli(0.2)]
        modelA = ConjugateBernoulli(1, 1)
        modelB = ConjugateBernoulli(1, 1)
        stoppingrule = ExpectedLossThresh(1e-3)
        experiment = ExperimentAB([modelA, modelB], stoppingrule)

        simulation = Simulation(
            experiment=experiment,
            datagendists=truedists,
            maxsteps=30,
            onestepsizes=[1000, 1000],
            minsteps=5
        )

        numsteps, winners, _ = runsequential(simulation, numsamples=5000, numsims=50)
        @test sum(winners .== "control") == 27
        @test sum(winners .== "variant 1") ==  22
        @test sum(winners .== "nothing") == 1
    end

    @testset "ChainedModel with expected loss with same distribution" begin
        Random.seed!(1234)
        truedists = [
            [Bernoulli(0.1), LogNormal(6.0, 0.8)],
            [Bernoulli(0.1), LogNormal(6.0, 0.8)]
        ]

        modelA = ChainedModel(
            [ConjugateBernoulli(1, 1), ConjugateLogNormal(0.0, 1.0, 0.001, 0.001)],
        )
        modelB = ChainedModel(
            [ConjugateBernoulli(1, 1), ConjugateLogNormal(0.0, 1.0, 0.001, 0.001)],
        )

        stoppingrule = ExpectedLossThresh(0.1)
        experiment = ExperimentAB([modelA, modelB], stoppingrule)

        simulation = Simulation(
            experiment=experiment,
            datagendists=truedists,
            maxsteps=30,
            onestepsizes=[[2000, 400], [1000, 200]]
        )

        numsteps, winners, _ = runsequential(simulation, numsamples=1000, numsims=20)
        @test sum(winners .== "control") ==  5
        @test sum(winners .== "variant 1") == 6
        @test sum(winners .== "nothing") ==  9
    end

end


text:
# functions related to displaying types

function show{M <: HMesh}(io::IO, m::M)
    println(io, "HomogenousMesh(")
    for (key,val) in attributes(m)
        print(io, "    ", key, ": ", length(val), "x", eltype(val), ", ")
    end
    println(io, ")")
end


text:
using ImageFiltering

@testset "SSIM" begin
    @info "test: SSIM"

    iqi_a = SSIM()
    iqi_b = SSIM(KernelFactors.gaussian(1.5, 11))
    iqi_c = SSIM(KernelFactors.gaussian(1.5, 11), (1.0, 1.0, 1.0))
    @test (iqi_a.kernel == iqi_b.kernel == iqi_c.kernel) &&
          (iqi_a.W == iqi_b.W == iqi_c.W)

    iqi = SSIM()
    sz_img_3 = (3, 3, 3)

    # numerical test
    img1 = testimage("cameraman")
    img2 = testimage("lena_gray_512")
    @test ssim(img1, img2) ≈ 0.3595 atol=1e-4 # MATLAB built-in ssim result
    iqi_δ = SSIM(KernelFactors.gaussian(1.5, 11), (1.0+1e-5, 1.0, 1.0))
    @test assess(iqi_δ, img1, img2) ≈ assess(SSIM(), img1, img2) atol = 1e-4

    # Gray image
    type_list = generate_test_types([Bool, Float32, N0f8], [Gray])
    A = [1.0 1.0 1.0; 1.0 1.0 1.0; 0.0 0.0 0.0]
    B = [1.0 1.0 1.0; 0.0 0.0 0.0; 1.0 1.0 1.0]
    for T in type_list
        test_ndarray(iqi, sz_img_3, T)

        a = A .|> T
        b = B .|> T

        @test ssim(a, b) == assess(SSIM(), a, b) == SSIM()(a, b)
        @test ssim(a, a) ≈ 1.0

        # FIXME: the result of Bool type is not strictly equal to others
        eltype(T) <: Bool && continue
        test_numeric(iqi, a, b, T)
        test_numeric(iqi, channelview(a), channelview(b), T)
    end
    test_cross_type(iqi, A, B, type_list)

    # RGB image
    img1 = testimage("mandril_color")
    img2 = testimage("lena_color_512")
    @test ssim(img1, img2) ≈ 0.0664 atol=1e-4 # MATLAB built-in ssim result

    type_list = generate_test_types([Float32, N0f8], [RGB])
    A = [RGB(0.0, 0.0, 0.0) RGB(0.0, 1.0, 0.0) RGB(0.0, 1.0, 1.0)
        RGB(0.0, 0.0, 1.0) RGB(1.0, 0.0, 0.0) RGB(1.0, 1.0, 0.0)
        RGB(1.0, 1.0, 1.0) RGB(1.0, 0.0, 1.0) RGB(0.0, 0.0, 0.0)]
    B = [RGB(0.0, 0.0, 0.0) RGB(0.0, 0.0, 1.0) RGB(1.0, 1.0, 1.0)
        RGB(0.0, 1.0, 0.0) RGB(1.0, 0.0, 0.0) RGB(1.0, 0.0, 1.0)
        RGB(0.0, 1.0, 1.0) RGB(1.0, 1.0, 0.0) RGB(0.0, 0.0, 0.0)]
    for T in type_list
        test_ndarray(iqi, sz_img_3, T)

        a = A .|> T
        b = B .|> T

        @test ssim(a, b) == assess(SSIM(), a, b) == SSIM()(a, b)
        @test ssim(a, b) == ssim(channelview(a), channelview(b))
        @test ssim(a, a) ≈ 1.0

        # RGB is treated as 3d gray image
        test_numeric(iqi, a, b, T)
        test_numeric(iqi, channelview(a), channelview(b), T; filename="references/SSIM_2d_RGB")
    end
    type_list = generate_test_types([Float32, N0f8], [RGB, BGR])
    test_cross_type(iqi, A, B, type_list)

    # general Color3 images
    type_list = generate_test_types([Float32], [Lab, HSV])
    A = [RGB(0.0, 0.0, 0.0) RGB(0.0, 1.0, 0.0) RGB(0.0, 1.0, 1.0)
        RGB(0.0, 0.0, 1.0) RGB(1.0, 0.0, 0.0) RGB(1.0, 1.0, 0.0)
        RGB(1.0, 1.0, 1.0) RGB(1.0, 0.0, 1.0) RGB(0.0, 0.0, 0.0)]
    B = [RGB(0.0, 0.0, 0.0) RGB(0.0, 0.0, 1.0) RGB(1.0, 1.0, 1.0)
        RGB(0.0, 1.0, 0.0) RGB(1.0, 0.0, 0.0) RGB(1.0, 0.0, 1.0)
        RGB(0.0, 1.0, 1.0) RGB(1.0, 1.0, 0.0) RGB(0.0, 0.0, 0.0)]
    for T in type_list
        a = A .|> T
        b = B .|> T

        @test_nowarn ssim(A, b), ssim(a, B)

        @test ssim(a, b) == assess(SSIM(), a, b) == SSIM()(a, b)
        @test ssim(A, A) ≈ 1.0

        # conversion to RGB first differs from no conversion
        @test ssim(a, b) ≠ ssim(channelview(a), channelview(b))
    end
    @test ssim(A, B) ≈ ssim(Lab.(A), B) atol=1e-4
end


text:
type Light
	position::Vec3f
	color::Vec4f
	brightness::Vec4f

	Light(position=Vec3f(),color=Vec4f(),brightness=ones(Vec4f)) = new(position,color,brightness)
end

getLightData(light::Light) = [
	light.position.x,light.position.y,light.position.z,
	light.color.x,light.color.y,light.color.z,light.color.a,
	light.brightness.x,light.brightness.y,light.brightness.z,light.brightness.a
	]

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
using NDTensors
using LinearAlgebra
using Test

@testset "BlockSparseTensor basic functionality" begin

  # Indices
  indsA = ([2, 3], [4, 5])

  # Locations of non-zero blocks
  locs = [(1, 2), (2, 1)]

  A = BlockSparseTensor(locs, indsA...)
  randn!(A)

  @test blockdims(A, (1, 2)) == (2, 5)
  @test blockdims(A, (2, 1)) == (3, 4)
  @test nnzblocks(A) == 2
  @test nnz(A) == 2 * 5 + 3 * 4
  @test inds(A) == ([2, 3], [4, 5])
  @test isblocknz(A, (2, 1))
  @test isblocknz(A, (1, 2))
  @test !isblocknz(A, (1, 1))
  @test !isblocknz(A, (2, 2))

  # Test different ways of getting nnz
  @test nnz(blockoffsets(A), inds(A)) == nnz(A)

  B = 2 * A
  @test B[1, 1] == 2 * A[1, 1]
  @test nnz(A) == 2 * 5 + 3 * 4
  @test nnz(B) == 2 * 5 + 3 * 4
  @test nnzblocks(A) == 2
  @test nnzblocks(B) == 2

  B = A / 2
  @test B[1, 1] == A[1, 1] / 2
  @test nnz(A) == 2 * 5 + 3 * 4
  @test nnz(B) == 2 * 5 + 3 * 4
  @test nnzblocks(A) == 2
  @test nnzblocks(B) == 2

  A[1, 5] = 15
  A[2, 5] = 25

  @test A[1, 1] == 0
  @test A[1, 5] == 15
  @test A[2, 5] == 25

  D = dense(A)

  @test D == A

  for I in eachindex(A)
    @test D[I] == A[I]
  end

  A12 = blockview(A, (1, 2))

  @test dims(A12) == (2, 5)

  for I in eachindex(A12)
    @test A12[I] == A[I + CartesianIndex(0, 4)]
  end

  B = BlockSparseTensor(undef, locs, indsA)
  randn!(B)

  C = A + B

  for I in eachindex(C)
    @test C[I] == A[I] + B[I]
  end

  Ap = permutedims(A, (2, 1))

  @test blockdims(Ap, (1, 2)) == (4, 3)
  @test blockdims(Ap, (2, 1)) == (5, 2)
  @test nnz(A) == nnz(Ap)
  @test nnzblocks(A) == nnzblocks(Ap)

  for I in eachindex(C)
    @test A[I] == Ap[NDTensors.permute(I, (2, 1))]
  end

  A = BlockSparseTensor(ComplexF64, locs, indsA)
  randn!(A)
  @test conj(data(store(A))) == data(store(conj(A)))
  @test typeof(conj(A)) <: BlockSparseTensor

  @testset "Random constructor" begin
    T = randomBlockSparseTensor([(1, 1), (2, 2)], ([2, 2], [2, 2]))
    @test nnzblocks(T) == 2
    @test nnz(T) == 8
    @test eltype(T) == Float64
    @test norm(T) ≉ 0

    Tc = randomBlockSparseTensor(ComplexF64, [(1, 1), (2, 2)], ([2, 2], [2, 2]))
    @test nnzblocks(Tc) == 2
    @test nnz(Tc) == 8
    @test eltype(Tc) == ComplexF64
    @test norm(Tc) ≉ 0
  end

  @testset "Complex Valued Operations" begin
    T = randomBlockSparseTensor(ComplexF64, [(1, 1), (2, 2)], ([2, 2], [2, 2]))
    rT = real(T)
    @test eltype(rT) == Float64
    @test nnzblocks(rT) == nnzblocks(T)
    iT = imag(T)
    @test eltype(iT) == Float64
    @test nnzblocks(iT) == nnzblocks(T)
    @test norm(rT)^2 + norm(iT)^2 ≈ norm(T)^2

    cT = conj(T)
    @test eltype(cT) == ComplexF64
    @test nnzblocks(cT) == nnzblocks(T)
  end

  @testset "BlockSparseTensor setindex! add block" begin
    T = BlockSparseTensor([2, 3], [4, 5])

    for I in eachindex(C)
      @test T[I] == 0.0
    end
    @test nnz(T) == 0
    @test nnzblocks(T) == 0
    @test !isblocknz(T, (1, 1))
    @test !isblocknz(T, (2, 1))
    @test !isblocknz(T, (1, 2))
    @test !isblocknz(T, (2, 2))

    T[1, 1] = 1.0

    @test T[1, 1] == 1.0
    @test nnz(T) == 8
    @test nnzblocks(T) == 1
    @test isblocknz(T, (1, 1))
    @test !isblocknz(T, (2, 1))
    @test !isblocknz(T, (1, 2))
    @test !isblocknz(T, (2, 2))

    T[4, 8] = 2.0

    @test T[4, 8] == 2.0
    @test nnz(T) == 8 + 15
    @test nnzblocks(T) == 2
    @test isblocknz(T, (1, 1))
    @test !isblocknz(T, (2, 1))
    @test !isblocknz(T, (1, 2))
    @test isblocknz(T, (2, 2))

    T[1, 6] = 3.0

    @test T[1, 6] == 3.0
    @test nnz(T) == 8 + 15 + 10
    @test nnzblocks(T) == 3
    @test isblocknz(T, (1, 1))
    @test !isblocknz(T, (2, 1))
    @test isblocknz(T, (1, 2))
    @test isblocknz(T, (2, 2))

    T[4, 2] = 4.0

    @test T[4, 2] == 4.0
    @test nnz(T) == 8 + 15 + 10 + 12
    @test nnzblocks(T) == 4
    @test isblocknz(T, (1, 1))
    @test isblocknz(T, (2, 1))
    @test isblocknz(T, (1, 2))
    @test isblocknz(T, (2, 2))
  end

  @testset "Add with different blocks" begin
    # Indices
    inds = ([2, 3], [4, 5])

    # Locations of non-zero blocks
    locsA = [(1, 1), (1, 2), (2, 2)]
    A = BlockSparseTensor(locsA, inds...)
    randn!(A)

    locsB = [(1, 2), (2, 1)]
    B = BlockSparseTensor(locsB, inds...)
    randn!(B)

    R = A + B

    @test nnz(R) == dim(R)
    for I in eachindex(R)
      @test R[I] == A[I] + B[I]
    end
  end

  @testset "permutedims!! with different blocks" begin
    # Indices
    indsA = ([2, 3], [4, 5])

    # Locations of non-zero blocks
    locsA = [(1, 2), (2, 1)]
    A = BlockSparseTensor(locsA, indsA...)
    randn!(A)

    perm = (2, 1)

    locsB = [(2, 1)]
    indsB = NDTensors.permute(indsA, perm)
    B = BlockSparseTensor(locsB, indsB...)
    randn!(B)

    R = permutedims!!(B, A, perm)

    @test nnz(R) == nnz(A)
    for I in eachindex(A)
      @test R[NDTensors.permute(I, perm)] == A[I]
    end
  end

  @testset "Contract" begin
    indsA = ([2, 3], [4, 5])
    locsA = [(1, 1), (2, 2), (2, 1), (1, 2)]
    A = BlockSparseTensor(locsA, indsA...)
    randn!(A)

    indsB = ([4, 5], [3, 2])
    locsB = [(1, 2), (2, 1), (1, 1)]
    B = BlockSparseTensor(locsB, indsB...)
    randn!(B)

    R = contract(A, (1, -1), B, (-1, 2))

    DA = dense(A)
    DB = dense(B)
    DR = contract(DA, (1, -1), DB, (-1, 2))

    for I in eachindex(R)
      @test R[I] ≈ DR[I]
    end
  end

  @testset "reshape" begin
    indsA = ([2, 3], [4, 5])
    locsA = [(2, 1), (1, 2)]
    A = BlockSparseTensor(locsA, indsA...)
    randn!(A)

    indsB = ([8, 12, 10, 15],)
    B = reshape(A, indsB)

    @test nnzblocks(A) == nnzblocks(B)
    @test nnz(A) == nnz(B)
    for (bA, bB) in zip(eachnzblock(A), eachnzblock(B))
      blockA = blockview(A, bA)
      blockB = blockview(B, bB)
      @test reshape(blockA, size(blockB)) == blockB
    end
  end

  @testset "permute_combine" begin
    indsA = ([2, 3], [4, 5], [6, 7, 8])
    locsA = [(2, 1, 1), (1, 2, 1), (2, 2, 3)]
    A = BlockSparseTensor(locsA, indsA...)
    randn!(A)

    B = NDTensors.permute_combine(A, 3, (2, 1))

    @test nnzblocks(A) == nnzblocks(B)
    @test nnz(A) == nnz(B)

    Ap = permutedims(A, (3, 2, 1))

    for (bAp, bB) in zip(eachnzblock(Ap), eachnzblock(B))
      blockAp = blockview(Ap, bAp)
      blockB = blockview(B, bB)
      @test reshape(blockAp, size(blockB)) == blockB
    end
  end

  @testset "svd" begin
    @testset "svd example 1" begin
      A = BlockSparseTensor([(2, 1), (1, 2)], [2, 2], [2, 2])
      randn!(A)
      U, S, V = svd(A)
      @test isapprox(norm(array(U) * array(S) * array(V)' - array(A)), 0; atol=1e-14)
    end

    @testset "svd example 2" begin
      A = BlockSparseTensor([(1, 2), (2, 3)], [2, 2], [3, 2, 3])
      randn!(A)
      U, S, V = svd(A)
      @test isapprox(norm(array(U) * array(S) * array(V)' - array(A)), 0.0; atol=1e-14)
    end

    @testset "svd example 3" begin
      A = BlockSparseTensor([(2, 1), (3, 2)], [3, 2, 3], [2, 2])
      randn!(A)
      U, S, V = svd(A)
      @test isapprox(norm(array(U) * array(S) * array(V)' - array(A)), 0.0; atol=1e-14)
    end

    @testset "svd example 4" begin
      A = BlockSparseTensor([(2, 1), (3, 2)], [2, 3, 4], [5, 6])
      randn!(A)
      U, S, V = svd(A)
      @test isapprox(norm(array(U) * array(S) * array(V)' - array(A)), 0.0; atol=1e-13)
    end

    @testset "svd example 5" begin
      A = BlockSparseTensor([(1, 2), (2, 3)], [5, 6], [2, 3, 4])
      randn!(A)
      U, S, V = svd(A)
      @test isapprox(norm(array(U) * array(S) * array(V)' - array(A)), 0.0; atol=1e-13)
    end
  end

  @testset "exp" begin
    A = BlockSparseTensor([(1, 1), (2, 2)], [2, 4], [2, 4])
    randn!(A)
    expT = exp(A)
    @test isapprox(norm(array(expT) - exp(array(A))), 0.0; atol=1e-13)

    # Hermitian case
    A = BlockSparseTensor(ComplexF64, [(1, 1), (2, 2)], ([2, 2], [2, 2]))
    randn!(A)
    Ah = BlockSparseTensor(ComplexF64, undef, [(1, 1), (2, 2)], ([2, 2], [2, 2]))
    for bA in eachnzblock(A)
      b = blockview(A, bA)
      blockview(Ah, bA) .= b + b'
    end
    expTh = exp(Hermitian(Ah))
    @test array(expTh) ≈ exp(Hermitian(array(Ah))) rtol = 1e-13

    A = BlockSparseTensor([(2, 1), (1, 2)], [2, 2], [2, 2])
    @test_throws ErrorException exp(A)
  end
end

nothing


text:
# Unique Work Queue

struct WorkQueue{T}
  items::Vector{T}
end

WorkQueue{T}() where T = WorkQueue(T[])

WorkQueue() = WorkQueue{Any}()

WorkQueue(xs) = WorkQueue(collect(xs))

function Base.push!(q::WorkQueue, x)
  i = findfirst(==(x), q.items)
  i === nothing || (deleteat!(q.items, i))
  push!(q.items, x)
  return q
end

Base.pop!(q::WorkQueue) = pop!(q.items)
Base.isempty(q::WorkQueue) = isempty(q.items)
Base.length(q::WorkQueue) = length(q.items)
Base.empty!(q::WorkQueue) = empty!(q.items)


text:
__precompile__()

module Jute

using Pkg
using Random
using BenchmarkTools

include("utils.jl")
export benchmark_result

include("options.jl")

include("fixture_factory.jl")

include("fixtures.jl")

include("testcases.jl")

include("macros.jl")
export @testcase
export @testgroup
export @global_fixture
export @local_fixture
export @produce

include("inspect.jl")

include("run_testcase.jl")
export @test
export @test_throws
export @test_broken
export @test_skip
export @test_result
export @test_fail
export @inferred
export @test_warn
export @test_nowarn
export @critical

include("reporting.jl")

include("runtests.jl")
export runtests

include("builtin_fixtures.jl")
export run_options
export temporary_dir
export fixed_rng

end


text:
using NLsolve

include("functions_withG.jl")

hh = Household();
r = 0.019;
A = 1.72;
T = 0.093;
ϕ = 342;
B = 1.0;

@unpack θ,δ,Amat,Ymat,na,nϵ,τ,γ,η,bc = hh
w = A_guess*(1-θ)*((r+δ)/(A*θ))^(θ/(θ-1));
cnext = @. r*Amat+w*Ymat+T

function getcBinding(a,c,ϵ)
    @. (1-τ)*w*ϵ*getl(c,ϵ,γ,ϕ,η,τ,w)+T+(1+(1-τ)*r)*a-bc-c;
end

cbinding = similar(Amat);

for i = 1:na
    for j = 1:nϵ
        cbinding[i,j] = NewtonRoot(c -> getcBinding(Amat[i,j],c,Ymat[i,j]),1.0);
    end
end

cpol,apol,lpol = egm(hh,T_guess,ϕ_guess,w=w,cnext = cnext,cbinding=cbinding,r=r,rnext=r)

cpol,apol,lpol = iterate_egm(hh,A_guess,T_guess,ϕ_guess;r=r_guess)

plot(hh.Amat[1:10,1], hh.Amat[1:10,1], label = "", color = :black, linestyle = :dash)
plot!(hh.Amat[1:10,1], apol[1:10,1], label = "ϵ = $(round(hh.Ymat[1,1], digits = 3))")
plot!(hh.Amat[1:10,1], apol[1:10,5], label = "ϵ = $(round(hh.Ymat[1,5], digits = 3))")

plot(hh.Amat[:,1], lpol[:,1], label = "ϵ = $(round(hh.Ymat[1,1], digits = 3))")
plot!(hh.Amat[:,1], lpol[:,2], label = "ϵ = $(round(hh.Ymat[1,2], digits = 3))")

r,K,N,w = market_clearing(hh,A_guess,r_guess,T_guess,B_guess,ϕ_guess)
plot_market_clearing(hh,A_guess,T_guess,B_guess,ϕ_guess)

function f!(F,x)
    r = x[1]
    A = x[2]
    T = x[3]
    ϕ = x[4]
    Y = 1.0
    θ = 0.3
    G = 0.2
    B = 1.0
    τ = 0.4
    r_eq,K,N,w = market_clearing(hh,A,r,T,B,ϕ)
    F[1] = r_eq-r
    F[2] = Y-A*(K^θ)*(N^(1-θ))
    F[3] = G+T+r-τ*(w*N+r*(K+B))
    F[4] = N-0.28
end

SS = nlsolve(f!,[0.012,5.0,0.05,30.0])


text:
@testset "autodiff" begin
    Random.seed!(0)

    # Should throw, as :wah is not a proper autodiff choice
    @test_throws ErrorException OnceDifferentiable(x->x, rand(10); autodiff=:wah)
    @test_throws ErrorException OnceDifferentiable(x->x, rand(10), 0.0; autodiff=:wah)
    @test_throws ErrorException TwiceDifferentiable(x->x, rand(10); autodiff=:wah)
    @test_throws ErrorException TwiceDifferentiable(x->x, rand(10), 0.0; autodiff=:wah)
    #@test_throws ErrorException TwiceDifferentiable(x->x, rand(10), 0.0, rand(10); autodiff=:wah)
    @test_throws ErrorException TwiceDifferentiable(x->x, x->x, rand(10); autodiff=:wah)
    @test_throws ErrorException TwiceDifferentiable(x->x, x->x, rand(10), 0.0; autodiff=:wah)
    #@test_throws ErrorException TwiceDifferentiable(x->x, x->x, rand(10), 0.0, rand(10); autodiff=:wah)

    for T in (OnceDifferentiable, TwiceDifferentiable)
        odad1 = T(x->5.0, rand(1); autodiff = :finite)
        odad2 = T(x->5.0, rand(1); autodiff = :forward)
        gradient!(odad1, rand(1))
        gradient!(odad2, rand(1))
        #    odad3 = T(x->5., rand(1); autodiff = :reverse)
        @test gradient(odad1) == [0.0]
        @test gradient(odad2) == [0.0]
        #    @test odad3.g == [0.0]
    end

    for a in (1.0, 5.0)
        xa = rand(1)
        odad1 = OnceDifferentiable(x->a*x[1], xa; autodiff = :finite)
        odad2 = OnceDifferentiable(x->a*x[1], xa; autodiff = :forward)
    #    odad3 = OnceDifferentiable(x->a*x[1], xa; autodiff = :reverse)
        gradient!(odad1, xa)
        gradient!(odad2, xa)
        @test gradient(odad1) ≈ [a]
        @test gradient(odad2) == [a]
    #    @test odad3.g == [a]
    end
    for a in (1.0, 5.0)
        xa = rand(1)
        odad1 = OnceDifferentiable(x->a*x[1]^2, xa; autodiff = :finite)
        odad2 = OnceDifferentiable(x->a*x[1]^2, xa; autodiff = :forward)
    #    odad3 = OnceDifferentiable(x->a*x[1]^2, xa; autodiff = :reverse)
        gradient!(odad1, xa)
        gradient!(odad2, xa)
        @test gradient(odad1) ≈ 2.0*a*xa
        @test gradient(odad2) == 2.0*a*xa
    #    @test odad3.g == 2.0*a*xa
    end
    nx = 2
    x = rand(nx)
    f(x) = sum(x.^3)
    fx = f(x)
    g(G, x) = copyto!(G, 3 .* x.^2)
    gx = g(NLSolversBase.alloc_DF(x, 0.0), x)
    h(H, x) = copyto!(H, Diagonal(6 .* x))
    hx = h(fill(0.0, nx, nx), x)
    for dtype in (OnceDifferentiable, TwiceDifferentiable)
        for autodiff in (:finite, :forward)
            # :forward should be exact, but :finite will not be
            differentiable = dtype(f, copy(x); autodiff = autodiff)
            value!(differentiable, copy(x))
            @test isapprox(value(differentiable), fx)
            clear!(differentiable)
            value_gradient!(differentiable, x)
            @test isapprox(value(differentiable), fx)
            @test isapprox(gradient(differentiable), gx)
            clear!(differentiable)
            gradient!(differentiable, x)
            @test isapprox(gradient(differentiable), gx)
            clear!(differentiable)
            if dtype == TwiceDifferentiable
                hessian!(differentiable, x)
                if autodiff == :finite
                    # we have to increase the tolerance here, as the hessian is
                    # not very accurate
                    @test isapprox(hessian(differentiable), hx; atol = 1e-6)
                else
                    @test hessian(differentiable) == hx
                end
            end
        end
    end
    for autodiff in (:finite, :forward)
        td = TwiceDifferentiable(x->sum(x), (G, x)->copyto!(G, fill!(copy(x),1)), copy(x); autodiff = autodiff)
        value(td)
        value!(td, x)
        value_gradient!(td, x)
        gradient!(td, x)
        hessian!(td, x)
    end
    for autodiff in (:finite, :forward)
        for nd = (NonDifferentiable(x->sum(x), copy(x)), NonDifferentiable(x->sum(x), copy(x), 0.0))
            td = TwiceDifferentiable(nd; autodiff = autodiff)
            value(td)
            value!(td, x)
            value_gradient!(td, x)
            gradient!(td, x)
            hessian!(td, x)
        end
        for od = (OnceDifferentiable(x->sum(x), (G, x)->copyto!(G, fill!(copy(x),1)), copy(x)), OnceDifferentiable(x->sum(x), (G, x)->copyto!(G, fill!(copy(x),1)), copy(x), 0.0))
            td = TwiceDifferentiable(od; autodiff = autodiff)
            value(td)
            value!(td, x)
            value_gradient!(td, x)
            gradient!(td, x)
            hessian!(td, x)
        end
    end
    @testset "autodiff ℝᴺ → ℝᴺ" begin
        function f!(F, x)
            F[1] = 1.0 - x[1]
            F[2] = 10.0*(x[2] - x[1]^2)
            F
        end
        function j!(J, x)
            J[1,1] = -1.0
            J[1,2] = 0.0
            J[2,1] = -20.0*x[1]
            J[2,2] = 10.0
            J
        end
        # Some random x
        x = rand(2)
        # A type and shape-correct F vector
        F = similar(x)
        # A type and shape-correct J matrix
        J = NLSolversBase.alloc_DF(x, F)
        # Default, should be :central
        od_fd = OnceDifferentiable(f!, x, F)
        value_jacobian!(od_fd, x)
        @test value(od_fd) == f!(F, x)
        # Can't test equality here
        @test jacobian(od_fd) ≈ j!(J, x)
        # Specifically :central
        od_fd_2 = OnceDifferentiable(f!, x, F, :central)
        @test value!(od_fd_2, x) == f!(F, x)
        # Can't test equality here
        @test jacobian!(od_fd_2, x) ≈ j!(J, x)
        # Test that they're identical -> they used the same scheme
        @test jacobian(od_fd) == jacobian(od_fd_2)

        od_ad = OnceDifferentiable(f!, x, F, :forward)
        @test value!(od_ad, x) == f!(F, x)
        # Can't test equality here
        @test jacobian!(od_ad, x) ≈ j!(J, x)

        od_ad_2 = OnceDifferentiable(f!, x, F, :forward)
        @test value!(od_ad_2, x) == f!(F, x)
        # Can't test equality here
        @test jacobian!(od_ad_2, x) ≈ j!(J, x)
        # Test that they're identical -> they used the same scheme
        @test jacobian(od_ad) == jacobian(od_ad_2)
        @testset "error handling" begin
            x = rand(2)
            F = similar(x)
            # Wrong symbol
            @test_throws ErrorException OnceDifferentiable(f!, x, F, :foo)
            # Wrong bool
            @test_throws ErrorException OnceDifferentiable(f!, x, F, false)
        end
    end
end
@testset "value/grad" begin
    a = 3.0
    x_seed = rand(1)
    odad1 = OnceDifferentiable(x->a*x[1]^2, x_seed)
    value_gradient!(odad1, x_seed)
    @test gradient(odad1) ≈ 2 .* a .* (x_seed)
    @testset "call counters" begin
        @test f_calls(odad1) == 1
        @test g_calls(odad1) == 1
        @test h_calls(odad1) == 0
        value_gradient!(odad1, x_seed .+ 1.0)
        @test f_calls(odad1) == 2
        @test g_calls(odad1) == 2
        @test h_calls(odad1) == 0
    end
    @test gradient(odad1) ≈ 2 .* a .* (x_seed .+ 1.0)
end

@testset "residual function" begin
    function f!(F, x)
        F[1] = 2x[1]+x[2]
        F[2] = x[1]+x[2]^2
        F[3] = x[1]^2+x[2]^2
        F
    end

    function j!(J, x)
        J[1, 1] = 2
        J[1, 2] = 1
        J[2, 1] = 1
        J[2, 2] = 2 * x[2]
        J[3, 1] = 2 * x[1]
        J[3, 2] = 2 * x[2]
        J
    end

    F = zeros(3)
    J = zeros(3, 2)
    x_init = [1., 1.]
    od = OnceDifferentiable(f!, x_init,  F)
    value_jacobian!(od, x_init)
    @test length(value(od)) == length(F)
    @test value(od) ≈ f!(F, x_init)
    @test jacobian(od) ≈ j!(J, x_init)
end


text:
using TaylorSeries, TaylorIntegration, Distributions
using Plots, LaTeXStrings
pyplot()

println("Using TaylorSeries, TaylorIntegration, Distributions, DataFrames, Plots, LaTeXStrings packages... \n")

## Custom libraries ###
include("JTFunctions.jl")
using .JTFunctions
include("ThreeBody.jl")
using .ThreeBody
import ThreeBody: μ

include("params_CR3BP.jl")

max_stars_repo_name: blas-ko/tesis

max_issues_repo_name: blas-ko/tesis

max_forks_repo_name: blas-ko/tesis


text:
# Return the vector in the box and the number of periods moved.
function period_check(vec::T, box::T) where T <: MosiVector
    t = round.(vec ./ box)
    vec -= t .* box
    vec, t
end

function update_periods!(vecs::AbstractVector{T}, periods::AbstractVector{T}, box::T) where T <: MosiVector
    N = length(vecs)
    @inbounds @simd for i in 1:N
        vec, t = period_check(vecs[i], box)
        vecs[i] = vec
        periods[i] += t
    end
    vecs, periods
end

default_distance(a::T, b::T) where T <: MosiVector = a - b
pbc_distance(pbox::T) where T <: MosiVector = (a::T, b::T) -> period_check(a - b, pbox)[1]

original_vector(vec::T, period::T, box::T) where T <: MosiVector = vec + period .* box

function original_vectors(vecs::AbstractVector{T}, periods::AbstractVector{T}, box::T) where T <: MosiVector
    N = length(vecs)
    if length(periods) ≡ 0
        return copy(vecs)
    end
    @inbounds [
        original_vector(vecs[i], periods[i], box)
        for i = 1:N
    ]
end

function pbc_box(N::Integer, ρ::Real, is3d = true)
    L = (N / ρ) ^ (1 / (is3d ? 3 : 2))
    is3d ? Vector3(L, L, L) : Vector2(L, L)
end


text:
using GnuplotRecipes
using Test
using Gnuplot, DataFrames, Measurements, PNGFiles, ColorTypes, FixedPointNumbers

function compare_plot(name)
    fn = "img/$name.png"
    fn_expected = "img/$name.expect.png"
    fn_diff = "img/$name.diff.png"
    isdir("img") || mkdir("img")
    save(term="pngcairo size 640,360 fontscale 0.8", output=fn)
    cur = PNGFiles.load(fn);
    if isfile(fn_expected)
        exp = PNGFiles.load(fn_expected);
    else
        exp = cur;
        cp(fn, fn_expected)
    end
    if exp != cur
        diff = RGB{N0f8}(1,1,1) .- (RGB{Float16}.(cur) .- RGB{Float16}.(exp) .|> abs .|> RGB{N0f8})
        PNGFiles.save(fn_diff, diff);
    end

    return exp == cur
end

macro gpok(args...)
    gp_call = :(@gp $(esc.(args)...))
    if ! Gnuplot.options.gpviewer
        gp_call = :(display($gp_call))
    end
    return quote
        $gp_call
        compare_plot($("gp_" * join(string.(args))))
    end
end

orig_term = Gnuplot.options.term
if !isinteractive()
    Gnuplot.options.term = "unknown"
end

@testset "GnuplotRecipes.jl" begin
    # Plain numbers
    dfnum = DataFrame(names=["a", "b", "c"], temp=10:12, speed=4:-1:2)
    @test @gpok bars(dfnum)
    # Measurement values
    df = DataFrame(names=["very long label", "b", "c"], temp=10:12, speed=collect(5:-1.5:2) .± 1)
    @test @gpok bars(df)
    @test @gpok bars(df, errorbars="lw 5 lc rgb '#ff0000'")
    @test @gpok bars(df, fill_style="pattern 1")
    @test @gpok bars(df, gap=2)
    @test @gpok bars(df, box_width=1)
    @test @gpok bars(df, label_rot=0)
    @test @gpok bars(df, label_rot=90) # visually check that labels are top-aligned
    @test @gpok bars(df, y2cols=[:speed])
    @test @gpok bars(df, y2cols=[:speed, :temp])
    @test @gpok bars(df, hist_style = "rowstacked")
    @test_logs (:warn, r"columnstacked") @gp bars(df, hist_style = "columnstacked")
    @test_throws AssertionError @gp bars(df, y2cols=[:bad])
    # Symbols as lables
    @test @gpok bars(DataFrame(l=[:a, :b], v=[2, 1]))
    # xlabelcol
    @test @gpok bars(DataFrame(v=[2, 1, 0.5], l=[:a, :b, :c]), xlabelcol=2)
    @test @gpok bars(DataFrame(v=[2, 1, 0.5], l=[:a, :b, :c]), xlabelcol=:l)
end

Gnuplot.options.term = orig_term


text:
#=

#Bagging of multiple boosted trees
function bagged_boosted_tree(mappings::Array{Array{String,1},1},candMatWOMaxValues::Array{Array{Float64,1},1},sett::ModelSettings,actualNumerator::Array{Float64,1},denominator::Array{Float64,1},weight::Array{Float64,1}, numfeatures::Array{PooledArray,1},charfeatures::Array{PooledArray,1},actualNumeratorVAL::Array{Float64,1},denominatorVAL::Array{Float64,1},weightVAL::Array{Float64,1}, numfeaturesVAL::Array{PooledArray,1},charfeaturesVAL::Array{PooledArray,1})
    error("this will not work as it has not been updated for DecisionTrees.")
	@assert sett.boolRandomizeOnlySplitAtTopNode=true #this must be true because of the selection of features in the tree construction ->    inds=randomFeatureSelection(settings.number_of_num_features,settings.number_of_char_features,settings.subsampling_features_prop)
	obs,obsVAL,meanobservedvalue,empty_rows_after_iteration_stats,showTimeUsedByEachIteration,chosen_apply_tree_fn,startAtMean,adaptiveLearningRate,moderationvector,iterations,nameOflistofpredictorsSheet,nameOfpredictorsSheet,nameOfModelStatisticsSheet,nameOfScoresSheet,nameOfOverviewSheet,nameOfSettingsSheet,nameOfValidationSheet,xlData,showProgressBar_time=initSettingsWhichAreTheSameForBoostingAndBagging(actualNumerator,actualNumeratorVAL,denominator,sett)
	cBB_niterBoosting=sett.cBB_niterBoosting
	cBB_niterBagging=sett.cBB_niterBagging

	iterationsPerCore=calcIterationsPerCore(cBB_niterBagging,Distributed.nprocs())
	sampleSizeForEachTreeCanBeNEG=convert(Int,round(sett.subsampling_prop*length(denominator)))

	vecBoostedTrees = Distributed.@distributed (vcat) for ii=1:Distributed.nprocs()
	#	ii=1;vecTreesWErrs=
		create_Subsamples_and_BoostedTrees(iterationsPerCore[ii],sampleSizeForEachTreeCanBeNEG,meanobservedvalue,candMatWOMaxValues,mappings,sett,actualNumerator,denominator,weight,numfeatures,charfeatures)
    end

end

function function create_Subsamples_and_BoostedTrees(itr::Int,sampleSizeCanBeNEGATIVE::Int,meanobservedvalue::Float64,candMatWOMaxValues::Array{Array{Float64,1},1},mappings::Array{Array{String,1},1},originalSett::ModelSettings,numerator::Array{Float64,1},denominator::Array{Float64,1},weight::Array{Float64,1},numfeatures::Array{PooledArray,1},charfeatures::Array{PooledArray,1})
#NOTE: we want to avoid that the individual trees of each boosted tree select variables
#Sett sampling settings for individual tree
	sett=deepcopy(originalSett)
	fixednuminds,fixedcharinds=randomFeatureSelection(length(numfeatures),length(charfeatures),originalSett.subsampling_features_prop)
	sett.fixedcharinds=deepcopy(fixedcharinds)
	sett.fixednuminds=deepcopy(fixednuminds)
	sett.subsampling_features_prop=1.0
#init sampledata
	abssampleSize,sampleVector,num,denom,w,numf,charf,ooBagnum,ooBagdenom,ooBagw,ooBagnumf,ooBagcharf,ooBagsize=initBoostrapSample(sampleSizeCanBeNEGATIVE,numerator,denominator,weight,charfeatures,numfeatures)
	leafNrooBag=Array{Int}(undef,ooBagsize)
	ooBagEstimates=Array{Float64}(undef,ooBagsize)

	for i=1:itr
	#sample data
    error("this will not work as it has not been updated for DecisionTrees.")
		unusedSamplePart=sampleData!(sampleSizeCanBeNEGATIVE,numerator,denominator,weight,numfeatures,charfeatures,sampleVector,num,denom,w,numf,charf,ooBagnum,ooBagdenom,ooBagw,ooBagnumf,ooBagcharf)
	#build boosted tree
		xlData,vectorOfLeafNumbersTrn,vectorOfLeafArrays,vectorOfRulePathsToLeavesArrays,rawObservedRatioPerScore,est_matrixFromScores,est_matrixFromScoresVAL,stats,scoresval,est_matrix_val,estimateUnsmoothedVal,estimateSmoothedVal,estimateFromRelativitiesVal,estimateUnsmoothedTrn,estimateSmoothedTrn,estimateFromRelativitiesTrn,resultEnsemble=boosted_tree(mappings,candMatWOMaxValues,sett,num,denom,w,numf,charf,ooBagnum,ooBagdenom,ooBagw,ooBagnumf,ooBagcharf)
	end
    error(".... in the works....")
	return true
end

#@show "this is the end of the file bagged_boosted_tree.jl" =#


text:
"""
Bridges `CP.GlobalCardinality{CP.VARIABLE_COUNTED_VALUES, CP.OPEN_COUNTED_VALUES, T}` to `CP.Count`.
"""
struct GlobalCardinalityVariableOpen2CountBridge{T} <: MOIBC.AbstractBridge
    cons_count::Vector{MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, CP.Count{MOI.EqualTo{T}}}}
end

function MOIBC.bridge_constraint(
    ::Type{GlobalCardinalityVariableOpen2CountBridge{T}},
    model,
    f::MOI.VectorOfVariables,
    s::CP.GlobalCardinality{CP.VARIABLE_COUNTED_VALUES, CP.OPEN_COUNTED_VALUES, T},
) where {T}
    return MOIBC.bridge_constraint(
        GlobalCardinalityVariableOpen2CountBridge{T},
        model,
        MOI.VectorAffineFunction{T}(f),
        s,
    )
end

function MOIBC.bridge_constraint(
    ::Type{GlobalCardinalityVariableOpen2CountBridge{T}},
    model,
    f::MOI.VectorAffineFunction{T},
    s::CP.GlobalCardinality{CP.VARIABLE_COUNTED_VALUES, CP.OPEN_COUNTED_VALUES, T},
) where {T}
    f_scalars = MOIU.scalarize(f)
    f_array = f_scalars[1:s.dimension]
    f_counts = f_scalars[(s.dimension + 1) : (s.dimension + s.n_values)]
    f_values = f_scalars[(s.dimension + s.n_values + 1) : (s.dimension + 2 * s.n_values)]

    cons_count = MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, CP.Count{MOI.EqualTo{T}}}[
        MOI.add_constraint(
            model,
            MOIU.vectorize(
                MOI.ScalarAffineFunction{T}[
                    f_counts[i],
                    (f_array .- f_values[i])...,
                ]
            ),
            CP.Count(s.dimension, MOI.EqualTo(zero(T)))
        )
        for i in 1:s.n_values
    ]

    return GlobalCardinalityVariableOpen2CountBridge(cons_count)
end

function MOI.supports_constraint(
    ::Type{GlobalCardinalityVariableOpen2CountBridge{T}},
    ::Union{Type{MOI.VectorOfVariables}, Type{MOI.VectorAffineFunction{T}}},
    ::Type{CP.GlobalCardinality{CP.VARIABLE_COUNTED_VALUES, CP.OPEN_COUNTED_VALUES, T}},
) where {T}
    return true
end

function MOIB.added_constrained_variable_types(::Type{GlobalCardinalityVariableOpen2CountBridge{T}}) where {T}
    return Tuple{Type}[]
end

function MOIB.added_constraint_types(::Type{GlobalCardinalityVariableOpen2CountBridge{T}}) where {T}
    return [
        (MOI.VectorAffineFunction{T}, CP.Count{MOI.EqualTo{T}}),
    ]
end

function MOI.get(::GlobalCardinalityVariableOpen2CountBridge{T}, ::MOI.NumberOfVariables) where {T}
    return 0
end

function MOI.get(
    b::GlobalCardinalityVariableOpen2CountBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{T}, CP.Count{MOI.EqualTo{T}},
    },
) where {T}
    return length(b.cons_count)
end

function MOI.get(
    b::GlobalCardinalityVariableOpen2CountBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VectorAffineFunction{T}, CP.Count{MOI.EqualTo{T}},
    },
) where {T}
    return copy(b.cons_count)
end