text:
function test_detector(detector)
    @testset "$detector" begin
        # create data
        is_supervised = isa(detector, SupervisedDetector)

        # specify test parameters
        rng = MersenneTwister(0)
        dim = trainDim
        fraction_train = 0.7
        n_samples = 100
        n_train = Int(n_samples * fraction_train)
        n_test = n_samples - n_train

        # test all different input formats
        y = rand(rng, (-1,0,1), n_samples)
        X_raw = rand(rng, dim, n_samples)
        X_mat = collect(X_raw')
        X_df = table(X_mat)
        X_dfs = source(X_df)
        ys = source(y)
        train, test  = partition(eachindex(y), fraction_train, rng=rng);

        # raw detector with matrix input
        model = is_supervised ? fit(detector, X_raw[:, train], y[train]) :
                                fit(detector, X_raw[:, train])
        train_raw, test_raw = score(detector, model, X_raw[:, test])

        # raw detector with table input
        model_df = is_supervised ? fit(detector, table(X_mat[train, :]), y[train]) :
                                   fit(detector, table(X_mat[train, :]))
        train_raw_df, test_raw_df = score(detector, model_df, table(X_mat[test, :]))

        # MLJ with matrix input
        detector_mat = is_supervised ? machine(detector, X_mat, y) : machine(detector, X_mat)
        fit!(detector_mat, rows=train)
        train_mat, test_mat = predict(detector_mat, rows=test)

        # MLJ with table input
        detector_df = is_supervised ? machine(detector, X_df, y) : machine(detector, X_df)
        fit!(detector_df, rows=train)
        train_df, test_df = predict(detector_df, rows=test)

        @testset "scores have appropriate dimensions" begin
            @test length(train_raw) == length(train_raw_df) == length(train_mat) == length(train_df) == n_train
            @test length(test_raw) == length(test_raw_df) == length(test_mat) == length(test_df) == n_test
        end

        @testset "scores have appropriate values" begin
            @test all(-Inf .< train_raw .< Inf)
            @test all(-Inf .< train_raw_df .< Inf)
            @test all(-Inf .< train_mat .< Inf)
            @test all(-Inf .< train_df .< Inf)

            @test all(-Inf .< test_raw .< Inf)
            @test all(-Inf .< test_raw_df .< Inf)
            @test all(-Inf .< test_mat .< Inf)
            @test all(-Inf .< test_df .< Inf)
        end
    end
end


text:
"""
    NNI(root::T, target::T, lor::Bool)::Int64   where T<:AbstractNode

This function does a nearest neighbour interchange (NNI) move on the tree specified
by `root`. The parameter `target` specifies the node which performs the interchange
move using the left or right child of the target node. If the left child should
be used `lor=true`.

The function returns 1 if the move was successful and 0 else.

* `root` : root node of tree on which to perform the NNI.

* `target` : specific node of tree to interchange.

* `lor` : Bool; "true" uses the left child of `target,` "false," the right child.
"""
function NNI!(root::T, target::T, lor::Bool)::Int64  where T<:AbstractNode
    # NNI move would be illegal
    if target.nchild == 0 || target.root
        return 0
    end # if

    parent::T = get_mother(target)
    sister::T = get_sister(target)

    ychild::T = remove_child!(target, lor)
    xchild::T = remove_child!(parent, sister)

    add_child!(target, sister)
    add_child!(parent, ychild)

    set_binary!(root)

    return 1

end # function

"""
    NNI!(root::T, target::Int64)::Int64  where T<:AbstractNode

This function does a nearest neighbour interchange (NNI) move on the tree specified
by `root`. The target is identified by the number of the target node.

The function returns 1 if the move was successful and 0 else.

* `root` : root node of tree on which to perform the NNI.

* `target` : specific node of tree to interchange.
"""
function NNI!(root::T, target::Int64, lor::Bool)::Int64  where T<:AbstractNode
   tn::T = find_num(root, target)
   NNI!(root, tn, lor)
end #function

"""
    NNI!(root::T, target::Int64)::Int64  where T<:AbstractNode

This function does a nearest neighbour interchange (NNI) move on the tree specified
by `root`. The target is identified by the number of the target node.
The function returns 1 if the move was successfull and 0 else.
"""
function NNI!(root::T, target::Int64)::Int64  where T<:AbstractNode
   tn::T = find_num(root, target)
   lor::Bool = 0.5 > rand()
   NNI!(root, tn, lor)
end #function

"""
    NNI!(root::T)::Int64  where T<:AbstractNode

This function does a nearest neighbour interchange (NNI) move on the tree specified
by `root`. The target is identified by the number of the target node.

The function returns 1 if the move was successful and 0 else.

* `root` : root node of tree on which to perform the NNI.

"""
function NNI!(root::T)::Int64  where T<:AbstractNode
    n = rand(1:size(root)[1])
    tn::T = find_num(root, n)
    lor::Bool = 0.5 > rand()
    NNI!(root, tn, lor)
end #function

"""
    NNI(root::T)::T  where T<:AbstractNode

This function does a nearest neighbour interchange (NNI) move on the tree specified
by `root.`

Returns a mutated copy while leaving the original tree intact.

* `root` : root node of tree on which to perform the NNI.

"""
function NNI(root::T)::T where T<:AbstractNode
    new_root = deepcopy(root)
    NNI!(new_root)
    return new_root
end

ext: jl


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

# helper, builds jacobian (d dx / d x) for FMUs
function _build_jac_dx_x(fmu::FMU2, rdx, rx)
    mat = zeros(length(rdx), length(rx))

    for i in 1:length(rdx)
        for j in 1:length(rx)
            mat[i,j] = fmi2GetDirectionalDerivative(fmu, rdx[i], rx[j])
        end
    end

    mat
end

"""
Performs the equivalent of fmiDoStep for ME-FMUs (note, that fmiDoStep is for CS-FMUs only).
Currently no event handling supported (but often not necessary).
Note, that event-handling during ME simulation in FMI.jl itself is supported.

Optional, FMU-values can be set via keyword arguments `setValueReferences` and `setValues`.
Optional, FMU-values can be retrieved by keyword argument `getValueReferences`.

Function returns an array with state derivatives and optionally the FMU-values for `getValueReferences`.
"""
# function fmi2DoStepME(fmu::FMU2,
#                       x::Array,
#                       t::Real = -1.0;
#                       setValueReferences=[],
#                       setValues=[],
#                       getValueReferences=[])
#     if t < 0.0
#         t = NeuralFMUGetCachedTime(fmu)
#     end
#     _fmi2DoStepME(fmu, x, t, setValueReferences, setValues, getValueReferences)
# end
function fmi2DoStepME(fmu::FMU2,
    x::Array,
    t::Real = -1.0,
    setValueReferences=[],
    setValues=[],
    getValueReferences=[])
    if t < 0.0
        t = NeuralFMUGetCachedTime(fmu)
    end
    _fmi2DoStepME(fmu, x, t, setValueReferences, setValues, getValueReferences)
end

# helper because keyword arguments are not supported by Zygote
function _fmi2DoStepME(fmu::FMU2,
                       x,
                       t::Real,
                       setValueReferences = [],
                       setValues = [],
                       getValueReferences = [])

    @assert fmi2IsModelExchange(fmu) ["fmi2DoStepME(...): As in the name, this function only supports ME-FMUs."]

    setter = (length(setValueReferences) > 0)
    getter = (length(getValueReferences) > 0)

    if setter
        @assert setValues != nothing ["fmi2DoStepME(...): `setValues` is nothing but `setValueReferences` is not!"]
        @assert length(setValueReferences) == length(setValues) ["fmi2DoStepME(...): `setValueReferences` and `setValues` need to be the same length!"]
    end

    fmi2SetTime(fmu, t)

    if setter
        fmi2SetReal(fmu, setValueReferences, setValues)
    end

    fmi2SetContinuousStates(fmu, x)

    fmi2CompletedIntegratorStep(fmu, fmi2True)

    y = []
    if getter
        y = fmi2GetReal(fmu, getValueReferences)
    end

    dx = fmi2GetDerivatives(fmu)

    [dx..., y...] # vcat(dx, y)
end

# The gradient for the function fmi2SetDoStepMEGet.
function _fmi2DoStepME_Gradient(c̄,
                                fmu::FMU2,
                                x,
                                t::Real,
                                setValueReferences = [],
                                setValues = [],
                                getValueReferences = [])

    setter = (length(setValueReferences) > 0)
    getter = (length(getValueReferences) > 0)

    if setter
        @assert setValues != nothing ["_fmi2DoStepME_Gradient(...): `setValues` is nothing but `setValueReferences` is not!"]
        @assert length(setValueReferences) == length(setValues) ["_fmi2DoStepME_Gradient(...): `setValueReferences` and `setValues` need to be the same length!"]
    end

    fmi2SetTime(fmu, t)
    fmi2SetContinuousStates(fmu, x)

    rdx = vcat(fmu.modelDescription.derivativeValueReferences, getValueReferences) #rdx = [fmu.modelDescription.derivativeValueReferences..., getValueReferences...]
    rx = fmu.modelDescription.stateValueReferences
    ru = setValueReferences

    n_dx_x = zeros(Float64, 0)
    if length(x) > 0
        mat_dx_x = _build_jac_dx_x(fmu, rdx, rx)
        n_dx_x = mat_dx_x' * c̄
    end

    n_dx_u = zeros(Float64, 0)
    if setter
        mat_dx_u = _build_jac_dx_x(fmu, rdx, ru)
        n_dx_u = mat_dx_u' * c̄
    end

    svr = 0.0 # zeros(Float64, length(setValueReferences), length(rdx)) * c̄
    gvr = 0.0 #zeros(Float64, length(getValueReferences), length(rdx)) * c̄

    return tuple(0.0, n_dx_x, 0.0, svr, n_dx_u, gvr)
end

# The adjoint connection between ME-function and function gradient.
@adjoint _fmi2DoStepME(fmu, x, t, setValueReferences, setValues, getValueReferences) = _fmi2DoStepME(fmu, x, t, setValueReferences, setValues, getValueReferences), c̄ -> _fmi2DoStepME_Gradient(c̄, fmu, x, t, setValueReferences, setValues, getValueReferences)

"""
Performs a fmiDoStep for CS-FMUs (note, that fmiDoStep is for CS-FMUs only).

Optional, FMU-values can be set via keyword arguments `setValueReferences` and `setValues`.
Optional, FMU-values can be retrieved by keyword argument `getValueReferences`.

Function returns the FMU-values for the optional `getValueReferences`.
"""
function fmi2DoStepCS(fmu::FMU2, dt; setValueReferences::Array = [], setValues::Array = [], getValueReferences::Array = [])
    _fmi2DoStepCS(fmu, dt, setValueReferences, setValues, getValueReferences)
end

function _fmi2DoStepCS(fmu::FMU2, dt, setValueReferences::Array = [], setValues::Array = [], getValueReferences::Array = [])

    @assert fmu.modelDescription.isCoSimulation == fmi2True ["fmi2DoStepCS(...): As in the name, this function only supports CS-FMUs."]
    @assert length(setValueReferences) == length(setValues) ["fmi2DoStepCS(...): `setValueReferences` ($(length(setValueReferences))) and `setValues` ($(length(setValues))) need to be the same length!"]

    if length(setValueReferences) > 0
        fmi2SetReal(fmu, setValueReferences, setValues)
    end

    t = fmu.t
    fmi2DoStep(fmu, t, dt)
    fmu.t += dt

    y = zeros(Float64, 0)

    if length(getValueReferences) > 0
        y = fmi2GetReal(fmu, getValueReferences)
    end

    y
end

# The gradient for the function fmi2InputDoStepCSOutput.
function _fmi2DoStepCS_Gradient(c̄, fmu::FMU2, dt, setValueReferences = [], setValues = [], getValueReferences = [])

    rdx = getValueReferences
    rx = setValueReferences

    n = zeros(Float64, 0)

    if length(getValueReferences) > 0
        mat =  _build_jac_dx_x(fmu, rdx, rx)
        n = mat' * c̄
    end

    nvr = zeros(Float64, length(setValueReferences))
    gvr = zeros(Float64, length(getValueReferences))

    tuple(nothing, 0.0, nvr, n, gvr)
end

# The adjoint connection between CS-function and function gradient.
@adjoint _fmi2DoStepCS(fmu, dt, setValueReferences, setValues, getValueReferences) = _fmi2DoStepCS(fmu, dt, setValueReferences, setValues, getValueReferences), c̄ -> _fmi2DoStepCS_Gradient(c̄, fmu, dt, setValueReferences, setValues, getValueReferences)

"""
Sets all FMU inputs to u, performs a ´´´fmi2DoStep´´´ and returns all FMU outputs.
"""
function fmi2InputDoStepCSOutput(fmu::FMU2, dt, u)
    _fmi2DoStepCS(fmu, dt,
                 fmu.modelDescription.inputValueReferences,
                 u,
                 fmu.modelDescription.outputValueReferences)
end

function fmi2InputDoStepCSOutput_Gradient(c̄, fmu::FMU2, dt, u)
    _fmi2DoStepCS_Gradient(c̄, fmu, dt,
                           fmu.modelDescription.inputValueReferences,
                           u,
                           fmu.modelDescription.outputValueReferences)
end

# The adjoint connection between CS-function and function gradient.
@adjoint fmi2InputDoStepCSOutput(fmu, dt, u) = fmi2InputDoStepCSOutput(fmu, dt, u), c̄ -> fmi2InputDoStepCSOutput_Gradient(c̄, fmu, dt, u)


text:
"""
`sciml_train`

### Unconstrained Optimization

```julia
function sciml_train(loss, _θ, opt = DEFAULT_OPT, adtype = DEFAULT_AD,
                     _data = DEFAULT_DATA, args...;
                     cb = (args...) -> false, maxiters = get_maxiters(data),
                     kwargs...)
```

### Box Constrained Optimization

```julia
function sciml_train(loss, θ, opt = DEFAULT_OPT, adtype = DEFAULT_AD,
                     data = DEFAULT_DATA, args...;
                     lower_bounds, upper_bounds,
                     cb = (args...) -> (false), maxiters = get_maxiters(data),
                     kwargs...)
```

## Optimizer Choices and Arguments

For a full definition of the allowed optimizers and arguments, please see the
[GalacticOptim.jl](https://galacticoptim.sciml.ai/dev/) documentation. As
sciml_train is an interface over GalacticOptim.jl, all of its optimizers and
arguments can be used from here.

## Loss Functions and Callbacks

Loss functions in `sciml_train` treat the first returned value as the return.
For example, if one returns `(1.0, [2.0])`, then the value the optimizer will
see is `1.0`. The other values are passed to the callback function. The callback
function is `cb(p, args...)` where the arguments are the extra returns from the
loss. This allows for reusing instead of recalculating. The callback function
must return a boolean where if `true`, then the optimizer will prematurely end
the optimization. It is called after every successful step, something that is
defined in an optimizer-dependent manner.

## Default AD Choice

The current default AD choice is dependent on the number of parameters.
For <50 parameters both ForwardDiff.jl and Zygote.jl gradients are evaluated
and the fastest is used. If both methods fail, finite difference method
is used as a fallback. For ≥50 parameters Zygote.jl is used.
More refinements to the techniques are planned.

## Default Optimizer Choice

By default, if the loss function is deterministic than an optimizer chain of
ADAM -> BFGS is used, otherwise ADAM is used (and a choice of maxiters is required).
"""
function sciml_train(loss, θ, opt=nothing, adtype=nothing, args...;
                     lower_bounds=nothing, upper_bounds=nothing,
                     maxiters=nothing, kwargs...)
    if adtype === nothing
        if length(θ) < 50
            fdtime = try
                ForwardDiff.gradient(x -> first(loss(x)), θ)
                @elapsed ForwardDiff.gradient(x -> first(loss(x)), θ)
            catch
                Inf
            end
            zytime = try
                Zygote.gradient(x -> first(loss(x)), θ)
                @elapsed Zygote.gradient(x -> first(loss(x)), θ)
            catch
                Inf
            end

            if fdtime == zytime == Inf
                @warn "AD methods failed, using numerical differentiation. To debug, try ForwardDiff.gradient(loss, θ) or Zygote.gradient(loss, θ)"
                adtype = GalacticOptim.AutoFiniteDiff()
            elseif fdtime < zytime
                adtype = GalacticOptim.AutoForwardDiff()
            else
                adtype = GalacticOptim.AutoZygote()
            end

        else
            adtype = GalacticOptim.AutoZygote()
        end
    end

    optf = GalacticOptim.OptimizationFunction((x, p) -> loss(x), adtype)
    optfunc = GalacticOptim.instantiate_function(optf, θ, adtype, nothing)
    optprob = GalacticOptim.OptimizationProblem(optfunc, θ; lb=lower_bounds, ub=upper_bounds, kwargs...)
    if opt !== nothing
        if maxiters !== nothing
            GalacticOptim.solve(optprob, opt, args...; maxiters, kwargs...)
        else
            GalacticOptim.solve(optprob, opt, args...; kwargs...)
        end
    else
        deterministic = first(loss(θ)) == first(loss(θ))
        if (!isempty(args) || !deterministic) && maxiters === nothing
            error("Automatic optimizer determination requires deterministic loss functions (and no data) or maxiters must be specified.")
        end

        if isempty(args) && deterministic
            # If determinsitic then ADAM -> finish with BFGS
            if maxiters === nothing
                res1 = GalacticOptim.solve(optprob, ADAM(0.01), args...; maxiters=300, kwargs...)
            else
                res1 = GalacticOptim.solve(optprob, ADAM(0.01), args...; maxiters, kwargs...)
            end

            optprob2 = GalacticOptim.OptimizationProblem(
                optfunc, res1.u; lb=lower_bounds, ub=upper_bounds, kwargs...)
            res1 = GalacticOptim.solve(
                optprob2, BFGS(initial_stepnorm=0.01), args...; maxiters, kwargs...)
        else
            res1 = GalacticOptim.solve(optprob, ADAM(0.1), args...; maxiters, kwargs...)
        end
    end
end


text:
# Activate necessary environment
import Pkg
Pkg.activate("splitdeliveries")

# Load all necessary packages (repositories)
using DataFrames
using DelimitedFiles
using CSV
using Combinatorics
using Distributions
using JuMP
using GAMS
using Octavian
using LoopVectorization
using Plots
using CUDA

# Load all necessary packages from the standard library module
using LinearAlgebra
using Random
using Statistics
using SparseArrays

## import all basic functions
include("functions/functions_basic.jl")
include("functions/functions_catalan.jl")
include("functions/functions_k-links.jl")
include("functions/functions_chisquare.jl")
include("functions/new_transactions.jl")

## import the heuristic functions
include("heuristics/heuristic_qmkp.jl")
include("heuristics/optimisation_equalcap.jl")
include("heuristics/optimisation_buffercap.jl")
include("heuristics/heuristic_klinks.jl")
include("heuristics/heuristic_greedyseeds.jl")
include("heuristics/heuristic_greedypairs.jl")
include("heuristics/heuristic_bestselling.jl")
include("heuristics/heuristic_chisquare.jl")

## import the main comparison function
include("main_benchmark.jl")

max_stars_repo_name: beyondatlas/splitdeliveries

max_issues_repo_name: beyondatlas/splitdeliveries

max_forks_repo_name: beyondatlas/splitdeliveries


text:
module MultivariateDiscretization

using DataFrames
using CSV
using StatsBase
using Discretizers
using Distributions
using Statistics
using LinearAlgebra
using MultivariateStats
using FreqTables
using Clustering
using InformationMeasures
using Distances

include("DoaneRule.jl")
include("BayesianBlocks.jl")
include("CPD.jl")
include("IPD.jl")
include("UniformWidth.jl")
include("benchmark_measures.jl")

export
    greedy_IPD,
    CPD,
    CPD_clustered,
    BayesianBlocks,
    DoaneRule,
    CJS_empirical,
    UW,
    greedy_IPD_clustered

end


text:
using Random
using YAML

export amplitudes_uniform, amplitudes_all
export generate_simulation_files, run_simulation
export generate_parameter_file, generate_dsl_files

"""
    amplitudes_all(qubits::Int)

Return generator for all the amplitudes of 'qubits' qubits
"""
function amplitudes_all(qubits::Int)
    (bitstring(x-1)[end-qubits+1:end] for x in 1:2^qubits)
end

"""
    amplitudes_uniform(qubits::Int, seed::Union{Int, Nothing}, number_amplitudes::Int)

Return generator to generate uniformly distributed bit strings of 'qubits' bits
(with replacement).
"""
function amplitudes_uniform(qubits::Int, seed::Union{Int, Nothing}, number_amplitudes::Int)
    rng = MersenneTwister(seed)
    bitstring = n -> join(rand(rng, ['0', '1'], n))
    (bitstring(qubits) for _ in 1:number_amplitudes)
end

"""
    generate_simulation_files(circ::QXZoo.Circuit.Circ,
                              output_prefix::String="simulation_input",
                              number_bonds_to_slice::Int=2;
                              decompose::Bool=true,
                              seed::Union{Int64, Nothing}=nothing,
                              output_args::Union{OrderedDict, Nothing}=nothing,
                              kwargs...)

Function to generate files required by qxrun to simulate the given circuit. This includes
a .tl file with a list of the operations, a .jld file with the initial tensors and a .yml
file with the parameters to use during the simulation.

# Keywords
- `number_bonds_to_slice::Int=2`: the number of edges to slice.
- `output_prefix::String="simulation_input"`: the prefix to be used for the simulation files.
- `num_amplitudes::Union{Int64, Nothing}=nothing`: the number of amplitudes to compute.
- `seed::Union{Int64, Nothing}=nothing`: the seed to be used by flow cutter to find a tree
                                         decomposition and for randomly selecting amplitudes to compute..
- `decompose::Bool=true`: set if two qubit gates should be decompoed when the circuit is converted to a tensor network.
- `kwargs`: all other kwargs are passed to `contraction_scheme` when it is called.
"""
function generate_simulation_files(circ::QXZoo.Circuit.Circ,
                                   output_prefix::String="simulation_input",
                                   number_bonds_to_slice::Int=2;
                                   decompose::Bool=true,
                                   seed::Union{Int64, Nothing}=nothing,
                                   output_args::Union{OrderedDict, Nothing}=nothing,
                                   kwargs...)

    @info("Convert circuit to tensor network")
    tnc = convert_to_tnc(circ; decompose=decompose)
    @info("Tensor network created with $(length(tnc)) tensors and $(length(bonds(tnc))) bonds")

    @info("Get contraction plan and edges to slice using QXGraphDecompositions")
    fc_seed = (seed === nothing) ? -1 : seed # flow cutter seed expects -1 in place of nothing.
    bond_groups, plan, metadata = contraction_scheme(tnc.tn, number_bonds_to_slice;
                                                     seed=fc_seed,
                                                     kwargs...)

    compute_tree = build_compute_graph(tnc, plan, bond_groups)

    @info("Prepare DSL and data files")
    generate_dsl_files(compute_tree, output_prefix; force=true, metadata=metadata)

    if output_args === nothing output_args = output_params_dict(qubits(tnc)) end
    generate_parameter_file(output_prefix, output_args)
end

"""
    single_amplitude(tnc::TensorNetworkCircuit, plan::Array{<:Index, 1}, amplitude::Union{String, Nothing}=nothing)

Contract the given tensor network using the given plan to calculate the amplitude of the given bitstring.
Creates a copy of the tensor network, replaces the outputs with those corresponding to the bitstring and then
contracts the network and returns the scalar amplitude.
"""
function single_amplitude(tnc::TensorNetworkCircuit, plan::Array{NTuple{3, Symbol}, 1}, amplitude::Union{String, Nothing}=nothing)
    sim_tnc = copy(tnc)
    add_output!(sim_tnc, amplitude)
    output = contract_tn!(sim_tnc, plan)
    output[1]
end

"""
    run_simulation(circ::QXZoo.Circuit.Circ;
                   num_amplitudes::Union{Int64, Nothing}=nothing,
                   seed::Union{Int64, Nothing}=nothing)

Function to run a simulation for the given circuit. If no number of amplitudes are given
then values for all possible amplitudes are found (only feasible for small systems). Returns
a dictionary of bitstring, amplitude pairs.
"""
function run_simulation(circ::QXZoo.Circuit.Circ;
                        num_amplitudes::Union{Int64, Nothing}=nothing,
                        seed::Union{Int64, Nothing}=nothing)
    @info("Convert circuit to tensor network")
    tnc = convert_to_tnc(circ)
    @info("Tensor network created with $(length(tnc)) tensors and $(length(bonds(tnc))) bonds")

    @info("Convert tensor network to graph")
    g = convert_to_graph(tnc)
    @info("Graph created: $(g)")

    @info("Get contraction plan and edges to slice using QXGraphDecompositions")
    plan = flow_cutter_contraction_plan(tnc; hypergraph=true)

    # to prevent memory issues don't attempt to get all amplitudes for large numbers of qubits
    # unless specifically demanded
    if num_amplitudes === nothing && qubits(tnc) > 30
        num_amplitudes = 1000
    end

    if num_amplitudes === nothing
        amplitudes = amplitudes_all(qubits(tnc))
    else
        amplitudes = amplitudes_uniform(qubits(tnc), seed, num_amplitudes)
    end

    results = OrderedDict{String, ComplexF64}()
    for amplitude in unique(amplitudes)
        results[amplitude] = single_amplitude(tnc, plan, amplitude)
    end
    results
end

"""
    generate_parameter_file(filename_prefix::String,
                            output_parameters)

Generate a yml file with details of how outputs are sampled

output:
    output_method: rejection
    params:
        fix_M: false
        M: 0.001
        num_samples: 10
        seed: ~
"""
function generate_parameter_file(filename_prefix::String,
                                 output_parameters)
    config = Dict()
    config["output"] = output_parameters

    YAML.write_file("$(filename_prefix).yml", config)
end


text:
function interp1_lin(fx::Array{Float64,1},fy::Array{Float64,1}, xx::Array{Float64,1} )
    if(!issorted(fx)); error("fx must be sorted ascending"); end
    yy = Array{Float64}(length(xx))
    for zi in 1:length(xx)
        xi = xx[zi]
        pxa = findlast( fx .<= xi  , true )
        pxb = findfirst( fx .>= xi  , true )
        if(pxa==pxb)
            yy[zi] = fy[pxa]
        else
            yy[zi] = fy[pxa] + (xi- fx[pxa] ) *  (fy[pxb]-fy[pxa]) / (fx[pxb]-fx[pxa])
        end
    end
    return yy
end

function bin2d( x::Array{Float64,2} , nbx::Int64 , nby::Int64 )
    xb = Array{Float64}(nbx,nby)
    limx = [ minimum(x[1,:]) ; maximum(x[1,:]) ]
    limy = [ minimum(x[2,:]) ; maximum(x[2,:]) ]
    bbx  = linspace(limx[1],limx[2],nbx+1)
    bby  = linspace(limy[1],limy[2],nby+1)
    for zi in 1:nbx
        for zj in 1:nby
            idx_ij = ( x[1,:].>= bbx[zi] ) .& ( x[1,:] .< bbx[zi+1] ) .& ( x[2,:] .>= bby[zj] ) .& ( x[2,:] .< bby[zj+1] )
            xb[zi,zj] = sum(idx_ij)
        end
    end
    return ( xb , 0.5.*(bbx[2:end]+bbx[1:end-1]) , 0.5.*(bby[2:end]+bby[1:end-1]) )
end

# sample x[:,i] has function value y[i]. The y corresponding to y coordinates will be binned
function bin2d( x::Array{Float64,2} , y::Array{Float64,1} , nbx::Int64 , nby::Int64 )
    xb = Array{Vector{Float64}}(nbx,nby)
    limx = [ minimum(x[1,:]) ; maximum(x[1,:]) ]
    limy = [ minimum(x[2,:]) ; maximum(x[2,:]) ]
    bbx  = linspace(limx[1],limx[2],nbx+1)
    bby  = linspace(limy[1],limy[2],nby+1)
    for zi in 1:nbx
        for zj in 1:nby
            idx_ij = ( x[1,:].>= bbx[zi] ) .& ( x[1,:] .< bbx[zi+1] ) .& ( x[2,:] .>= bby[zj] ) .& ( x[2,:] .< bby[zj+1] )
            xb[zi,zj] = Vector{Float64}()
            append!(xb[zi,zj], y[idx_ij] )
        end
    end
    return (xb, 0.5.*(bbx[2:end]+bbx[1:end-1]) , 0.5.*(bby[2:end]+bby[1:end-1]) )
end

function bin2d( x::Array{Float64,2} , y::Array{Float64,1} , nbx::Int64 , nby::Int64 , op::Any)
    b2dr_a = bin2d(x::Array{Float64,2} , y::Array{Float64,1} , nbx::Int64 , nby::Int64)
    b2dr = b2dr_a[1]
    opx  = zeros(nby,nbx)
    for zi=1:nbx
        for zj=1:nby
            opx[zi,zj] = op(b2dr[zi,zj])
        end
    end
    return (opx , b2dr_a[2] , b2dr_a[3] )
end


text:
module MaxHelpingHandFlagToggleSmoothCameraOffsetTrigger

using ..Ahorn, Maple

@mapdef Trigger "MaxHelpingHand/FlagToggleSmoothCameraOffsetTrigger" FlagToggleSmoothCameraOffsetTrigger(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight,
    offsetXFrom::Number=0.0, offsetXTo::Number=0.0, offsetYFrom::Number=0.0, offsetYTo::Number=0.0, positionMode::String="NoEffect", onlyOnce::Bool=false, xOnly::Bool=false, yOnly::Bool=false,
    flag::String="flag_toggle_smooth_camera_offset", inverted::Bool=false)

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

function Ahorn.editingOptions(trigger::FlagToggleSmoothCameraOffsetTrigger)
    return Dict{String, Any}(
        "positionMode" => Maple.trigger_position_modes
    )
end

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
function AutomotiveDrivingModels.observe!(lidar::LidarSensor, scene::Scene, env::OccludedEnv, egoid::Int64)
    state_ego = scene[findfirst(egoid, scene)].state
    ego_vel = polar(state_ego.v, state_ego.posG.θ)
    for (i,angle) in enumerate(lidar.angles)
        ray_angle = state_ego.posG.θ + angle
        ray_vec = polar(1.0, ray_angle)
        ray = VecSE2(state_ego.posG.x, state_ego.posG.y, ray_angle)

        range = lidar.max_range
        range_rate = 0.0
        for veh in scene
            if veh.id != egoid
                to_oriented_bounding_box!(lidar.poly, veh)

                range2 = AutomotiveDrivingModels.get_collision_time(ray, lidar.poly, 1.0)
                if !isnan(range2) && range2 < range
                    range = range2
                    relative_speed = polar(veh.state.v, veh.state.posG.θ) - ego_vel
                    range_rate = proj(relative_speed, ray_vec, Float64)
                end
            end
        end
        range_obs = lidar.max_range
        for obs in env.obstacles
            range_obs2 = AutomotiveDrivingModels.get_collision_time(ray, obs, 1.0)
            if !isnan(range_obs2) && range_obs2 < range_obs
                range_obs = range_obs2
                relative_speed = polar(0., 0.) - ego_vel
                range_rate = proj(relative_speed, ray_vec, Float64)
            end
        end

        if range_obs < range
            lidar.ranges[i] = range_obs
            lidar.range_rates[i] = range_rate
        else
            lidar.ranges[i] = range
            lidar.range_rates[i] = range_rate
        end
    end

    lidar
end

mutable struct LidarOverlay <: SceneOverlay
    sensor::LidarSensor
end

function render_ray!(rendermodel::RenderModel, ray::VecSE2;
    color::Colorant=convert(RGB, HSV(rand(1:360), 0.8, 0.8)),
    origin_radius::Float64 = 0.1,
    length::Float64 = 2.0,
    line_width::Float64 = 0.05,
    )

    pos2 = ray + polar(length, ray.θ)
    add_instruction!(rendermodel, render_line_segment, (ray.x, ray.y, pos2.x, pos2.y, color, line_width))
    add_instruction!(rendermodel, render_circle, (ray.x, ray.y, origin_radius, color, color, 0.0))
    rendermodel
end

function render_lidar!(rendermodel::RenderModel, lidar::LidarSensor, posG::VecSE2;
    color::Colorant=colorant"white",
    line_width::Float64 = 0.05,
    )
    for (angle, range, range_rate) in zip(lidar.angles, lidar.ranges, lidar.range_rates)
        ray = VecSE2(posG.x, posG.y, posG.θ + angle)
#         println(ray, range)
        t = 2 / (1 + exp(-range_rate))
        ray_color = t > 1.0 ? lerp(RGBA(1.0,1.0,1.0,0.5), RGBA(0.2,0.2,1.0,0.5), t/2) :
                              lerp(RGBA(1.0,0.2,0.2,0.5), RGBA(1.0,1.0,1.0,0.5), t)
        render_ray!(rendermodel::RenderModel, ray, color=ray_color, length=range, line_width=line_width)
    end
    rendermodel
end

function AutoViz.render!(rendermodel::RenderModel, overlay::LidarOverlay, scene::Scene, env::OccludedEnv)
    ego = scene[findfirst(EGO_ID, scene)]
    observe!(overlay.sensor, scene, env, ego.id)
    return render_lidar!(rendermodel, overlay.sensor, ego.state.posG)
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
using ClimateMachine
using ClimateMachine.Atmos
using ClimateMachine.BalanceLaws
using ClimateMachine.ConfigTypes
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.FVReconstructions: FVLinear
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.TemperatureProfiles
using ClimateMachine.GenericCallbacks
using ClimateMachine.Mesh.Geometry
using ClimateMachine.Mesh.Grids
using ClimateMachine.Mesh.Topologies
using ClimateMachine.MPIStateArrays
using ClimateMachine.ODESolvers
using ClimateMachine.Orientations
using ClimateMachine.Thermodynamics
using ClimateMachine.TurbulenceClosures
using ClimateMachine.VariableTemplates

using CLIMAParameters
using CLIMAParameters.Planet: planet_radius
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

using Test, MPI, Logging, StaticArrays, LinearAlgebra, Printf, Dates

function main()
    ClimateMachine.init()
    ArrayType = ClimateMachine.array_type()

    mpicomm = MPI.COMM_WORLD

    polynomialorder = (4, 0)
    FT = Float64
    NumericalFlux = RoeNumericalFlux
    @info @sprintf """Configuration
                      ArrayType     = %s
                      FT            = %s
                      NumericalFlux = %s
                      """ ArrayType FT NumericalFlux

    numelem_horz = 10
    numelem_vert = 32

    @testset for domain_type in (:box, :sphere)
        err = test_run(
            mpicomm,
            ArrayType,
            polynomialorder,
            numelem_horz,
            numelem_vert,
            NumericalFlux,
            FT,
            domain_type,
        )
        @test err < FT(6e-12)
    end
end

function test_run(
    mpicomm,
    ArrayType,
    polynomialorder,
    numelem_horz,
    numelem_vert,
    NumericalFlux,
    FT,
    domain_type,
)
    domain_height = 10e3
    if domain_type === :box
        domain_width = 20e3
        horz_range =
            range(FT(0), length = numelem_horz + 1, stop = FT(domain_width))
        vert_range = range(0, length = numelem_vert + 1, stop = domain_height)
        brickrange = (horz_range, horz_range, vert_range)
        periodicity = (true, true, false)
        topology =
            StackedBrickTopology(mpicomm, brickrange; periodicity = periodicity)
        meshwarp = (x...) -> identity(x)
    elseif domain_type === :sphere
        _planet_radius::FT = planet_radius(param_set)
        vert_range = grid1d(
            _planet_radius,
            FT(_planet_radius + domain_height),
            nelem = numelem_vert,
        )
        topology = StackedCubedSphereTopology(mpicomm, numelem_horz, vert_range)
        meshwarp = equiangular_cubed_sphere_warp
    end

    grid = DiscontinuousSpectralElementGrid(
        topology,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = polynomialorder,
        meshwarp = meshwarp,
    )

    T0 = FT(300)
    temp_profile = IsothermalProfile(param_set, T0)
    ref_state = HydrostaticState(temp_profile; subtract_off = false)

    physics = AtmosPhysics{FT}(
        param_set;
        ref_state = ref_state,
        turbulence = ConstantDynamicViscosity(FT(0)),
        moisture = DryModel(),
    )

    problem = AtmosProblem(;
        physics = physics,
        init_state_prognostic = initialcondition!,
    )

    if domain_type === :box
        configtype = AtmosLESConfigType
        source = (Gravity(),)
    elseif domain_type === :sphere
        configtype = AtmosGCMConfigType
        source = (Gravity(), Coriolis())
    end

    model =
        AtmosModel{FT}(configtype, physics; problem = problem, source = source)

    dg = DGFVModel(
        model,
        grid,
        HBFVReconstruction(model, FVLinear()),
        NumericalFlux(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient(),
    )

    timeend = FT(100)

    # determine the time step
    cfl = 0.2
    dz = step(vert_range)
    dt = cfl * dz / soundspeed_air(param_set, T0)
    nsteps = ceil(Int, timeend / dt)
    dt = timeend / nsteps

    Q = init_ode_state(dg, FT(0))
    lsrk = LSRK54CarpenterKennedy(dg, Q; dt = dt, t0 = 0)

    eng0 = norm(Q)
    @info @sprintf """Starting
                      domain_type   = %s
                      numelem_horz  = %d
                      numelem_vert  = %d
                      dt            = %.16e
                      norm(Q₀)      = %.16e
                      """ domain_type numelem_horz numelem_vert dt eng0

    # Set up the information callback
    starttime = Ref(now())
    cbinfo = EveryXWallTimeSeconds(60, mpicomm) do (s = false)
        if s
            starttime[] = now()
        else
            energy = norm(Q)
            @views begin
                ρu = extrema(Array(Q.data[:, 2, :]))
                ρv = extrema(Array(Q.data[:, 3, :]))
                ρw = extrema(Array(Q.data[:, 4, :]))
            end
            runtime = Dates.format(
                convert(DateTime, now() - starttime[]),
                dateformat"HH:MM:SS",
            )
            @info @sprintf """Update
                              simtime = %.16e
                              runtime = %s
                              ρu = %.16e, %.16e
                              ρv = %.16e, %.16e
                              ρw = %.16e, %.16e
                              norm(Q) = %.16e
                              """ gettime(lsrk) runtime ρu... ρv... ρw... energy
        end
    end
    callbacks = (cbinfo,)

    solve!(Q, lsrk; timeend = timeend, callbacks = callbacks)

    # final statistics
    Qe = init_ode_state(dg, timeend)
    engf = norm(Q)
    engfe = norm(Qe)
    errf = euclidean_distance(Q, Qe)
    errr = errf / engfe
    @info @sprintf """Finished
    norm(Q)                 = %.16e
    norm(Q) / norm(Q₀)      = %.16e
    norm(Q) - norm(Q₀)      = %.16e
    norm(Q - Qe)            = %.16e
    norm(Q - Qe) / norm(Qe) = %.16e
    """ engf engf / eng0 engf - eng0 errf errr
    errr
end

function initialcondition!(problem, bl, state, aux, coords, t, args...)
    state.ρ = aux.ref_state.ρ
    state.ρu = SVector(0, 0, 0)
    state.energy.ρe = aux.ref_state.ρe
end

main()


text:
using LinearAlgebra, FastSolversForWeightedTV, Flux, TestImages, CUDA, PyPlot
CUDA.allowscalar(false)
using Random; Random.seed!(123)

# Random data
# dim = 2
dim = 3
n = tuple(256*ones(Int64, dim)...)
y_orig = Float32.(TestImages.shepp_logan(n[1:2]...))
dim == 3 && (y_orig = repeat(y_orig; outer=(1,1,n[3])))
y_orig = y_orig |> gpu
y_orig = y_orig/norm(y_orig, Inf)

# Weight for structural norm
# η = 0.1f0*structural_mean(y_orig)
# P = structural_weight(y_orig; η=η)
P = nothing

# Gradient norm
N = (2,1); method = "TV"
# N = (2,2); method = "L2V"
# N = (2,Inf); method = "LInfV"
~isnothing(P) && (method = string(method, "w"))
method = string(method, "_", dim, "D")
h = tuple(ones(Float32, dim)...)
g = gradient_norm(N..., n, h; weight=P) |> gpu
ε_orig = g(y_orig)
dim == 2 && (steplength = 1f0/8f0)
dim == 3 && (steplength = 1f0/12f0)
opt = opt_fista(steplength; niter=1000, tol_x=nothing, nesterov=true)

# Artificial noise
y = y_orig+gpu(0.1f0*randn(Float32, n))

# Proxy
λ = 0.5f0*norm(y)^2/g(y)
xproxy = proxy(y, λ, g, opt)
err_proxy = norm(xproxy-y_orig)/norm(y_orig)
println("Proxy rel err: ", err_proxy)

# Projection test
ε = 0.5f0*ε_orig
xproj = project(y, ε, g, opt)
err_proj = norm(xproj-y_orig)/norm(y_orig)
println("Proj rel err: ", err_proj)

# Plot proxy
vmin = 0f0
vmax = 1f0
close("all")
figure()
subplot(1,2,1)
dim == 2 && imshow(cpu(y); cmap="gray", vmin=vmin, vmax=vmax)
dim == 3 && imshow(cpu(y[:,:,128]); cmap="gray", vmin=vmin, vmax=vmax)
title("Noisy phantom")
axis("off")
subplot(1,2,2)
dim == 2 && imshow(cpu(xproxy); cmap="gray", vmin=vmin, vmax=vmax)
dim == 3 && imshow(cpu(xproxy[:,:,128]); cmap="gray", vmin=vmin, vmax=vmax)
title(string(method, " proxy, ", L"$\mathrm{err}_{\mathrm{rel}}$ = ", string(err_proxy)))
axis("off")
savefig(string("./plots/", method, "_prox.png"), dpi=300, transparent=false, bbox_inches="tight")

# Plot proxy
figure()
subplot(1,2,1)
dim == 2 && imshow(cpu(y); cmap="gray", vmin=vmin, vmax=vmax)
dim == 3 && imshow(cpu(y[:,:,128]); cmap="gray", vmin=vmin, vmax=vmax)
title("Noisy phantom")
axis("off")
subplot(1,2,2)
dim == 2 && imshow(cpu(xproj); cmap="gray", vmin=vmin, vmax=vmax)
dim == 3 && imshow(cpu(xproj[:,:,128]); cmap="gray", vmin=vmin, vmax=vmax)
title(string(method, " proj., ", L"$\mathrm{err}_{\mathrm{rel}}$ = ", string(err_proj)))
axis("off")
savefig(string("./plots/", method, "_proj.png"), dpi=300, transparent=false, bbox_inches="tight")

max_stars_repo_name: grizzuti/FastSolversForWeightedTV

max_issues_repo_name: grizzuti/FastSolversForWeightedTV

max_forks_repo_name: grizzuti/FastSolversForWeightedTV


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

JLLWrappers.@generate_main_file_header("glmnet")
JLLWrappers.@generate_main_file("glmnet", UUID("78c6b45d-5eaf-5d68-bcfb-a5a2cb06c27f"))
end  # module glmnet_jll


text:
# returns a julia centric version of mujoco's model and data fields
# that allows direct access to Ptr array fields as julia vectors
function wrap_array(p, s1::Int, s2::Int, f_type)
   len = s1 * s2
   raw = unsafe_wrap(Array, p, len)

   if f_type == Matrix{eltype(f_type)}
      raw = reshape(raw, s2, s1)
   end
   return raw
end

function mapmodel(pm::Ptr{mjModel})
   c_model= unsafe_load(pm)
   mapmodel(c_model, pm)
end
function mapmodel(c_model::mjModel, pm::Ptr{mjModel})
   margs = Vector{Any}()

   m_fields = fieldnames(jlModel)[2:end] # drop :m pointer
   m_sizes = getmodelsize(c_model)
   jminfo = structinfo(jlModel)
   maxmodelmemptr = convert(UInt64, getfield(c_model, :names))
   for f in m_fields
      adr = convert(UInt64, getfield(c_model, f))
      if adr == 0x0 || adr > maxmodelmemptr # bad pointer
         m_off, m_type = jminfo[f]
         push!(margs, m_type(0) )
      else
         raw = wrap_array(getfield(c_model, f), Int(m_sizes[f][1]), Int(m_sizes[f][2]), fieldtype(jlModel, f))
         push!(margs, raw)
      end
   end
   return jlModel(Base.RefValue(unsafe_load(pm)), margs...)
end

function mapdata(pm::Ptr{mjModel}, pd::Ptr{mjData})
   c_model = unsafe_load(pm)
   mapdata(c_model, pd)
end
function mapdata(c_model::mjModel, pd::Ptr{mjData})
   c_data = unsafe_load(pd)

   dargs = Vector{Any}()

   d_fields = fieldnames(jlData)[2:end] # drop :d pointer
   d_sizes = getdatasize(c_model, c_data)
   for f in d_fields
      raw = wrap_array(getfield(c_data, f), Int(d_sizes[f][1]), Int(d_sizes[f][2]), fieldtype(jlData, f))
      push!(dargs, raw)
   end
   return jlData(Base.RefValue(unsafe_load(pd)), dargs...)
end

function mapmujoco(pm::Ptr{mjModel}, pd::Ptr{mjData})
   c_model = unsafe_load(pm)
   return mapmodel(c_model, pm), mapdata(c_model, pd)
end

# struct manipulation and access
structinfo(T) = Dict(fieldname(T,i)=>(fieldoffset(T,i), fieldtype(T,i)) for i = 1:fieldcount(T))
const minfo = structinfo(mjModel)
const dinfo = structinfo(mjData)

const mjstructs = Dict(mjContact     => structinfo(mjContact),
                       mjWarningStat => structinfo(mjWarningStat),
                       mjTimerStat   => structinfo(mjTimerStat),
                       mjSolverStat  => structinfo(mjSolverStat),

                       mjrContext    => structinfo(mjrContext),

                       mjVFS         => structinfo(mjVFS),
                       mjOption      => structinfo(mjOption),
                       #_global       => structinfo(_global),
                       #quality       => structinfo(quality),
                       #headlight     => structinfo(headlight),
                       #map           => structinfo(map),
                       #scale         => structinfo(scale),
                       #rgba          => structinfo(rgba),
                       mjVisual      => structinfo(mjVisual),
                       mjStatistic   => structinfo(mjStatistic),
                       mjModel       => structinfo(mjModel),

                       mjvPerturb    => structinfo(mjvPerturb),
                       mjvCamera     => structinfo(mjvCamera),
                       mjvGLCamera   => structinfo(mjvGLCamera),
                       mjvGeom       => structinfo(mjvGeom),
                       mjvLight      => structinfo(mjvLight),
                       mjvOption     => structinfo(mjvOption),
                       mjvScene      => structinfo(mjvScene),
                       mjvFigure     => structinfo(mjvFigure))

# access mujoco struct fields through the julia version of model and data
myconvert(T, m::Base.RefValue{mjModel})::T = convert(T, unsafe_convert(Ptr{mjModel}, m))
myconvert(T, d::Base.RefValue{mjData})::T = convert(T, unsafe_convert(Ptr{mjData}, d))
addr_to_val(T::Type{Float64}, p::UInt64)::Float64 = unsafe_load(convert(Ptr{Float64}, p))
addr_to_val(T::Type{Float32}, p::UInt64)::Float32 = unsafe_load(convert(Ptr{Float32}, p))
addr_to_val(T::Type{Int32},   p::UInt64)::Int32   = unsafe_load(convert(Ptr{Int32},   p))
addr_to_val(T::Type{UInt8},   p::UInt64)::UInt8   = unsafe_load(convert(Ptr{UInt8},   p))
addr_to_val(T::Type{Int8},    p::UInt64)::Int8    = unsafe_load(convert(Ptr{Int8},    p))

function get(m::jlModel, field::Symbol)
   f_off, f_type = minfo[field]
   pntr = myconvert(UInt64, m.m) + f_off
   addr_to_val(f_type, pntr)
end
function get(m::jlModel, fstruct::Symbol, field::Symbol)
   s_off, s_type = minfo[fstruct]
   #@assert s_type in (mjOption, mjVisual, mjStatistic)

   f_off, f_type = mjstructs[s_type][field]
   pntr = myconvert(UInt64, m.m) + s_off + f_off
   addr_to_val(f_type, pntr)
end
function get(d::jlData, field::Symbol)
   f_off, f_type = dinfo[field]
   pntr = myconvert(UInt64, d.d) + f_off
   addr_to_val(f_type, pntr)
end
function get(p::Ptr{T}, field::Symbol) where {T}
   f_off, f_type = mjstructs[T][field]
   pntr = Ptr{f_type}(p)
   return unsafe_load(pntr+f_off, 1)
end
function get(p::Ptr{T}, field::Symbol, i::Int) where {T}
   f_off, f_type = mjstructs[T][field]
   ET = eltype(f_type)
   @assert f_type <: SVector
   pntr = Ptr{ET}(p)
   unsafe_load(pntr+f_off, i)
end
function get(p::Ptr{T}, field::Symbol, i::Int, j::Int) where T # does row-col conversion
   f_off, f_type = mjstructs[T][field]
   ET = eltype(eltype(f_type))
   @assert f_type <: SVector && eltype(f_type) <: SVector
   #r, c = size(f_type)
   c = sizeof(eltype(f_type))[1]
   #@assert i <= r && i >= 1
   #@assert j <= c && i >= 1
   #idx = (i-1)*c + (j-1)*sizeof(ET)
   idx = (j-1)*sizeof(ET)
   pntr = Ptr{ET}(p)
   unsafe_load(pntr+f_off + (i-1)*c + (j-1)*sizeof(ET), 1)
end

function update_ptr(p::Ptr, offset::Integer, val::Integer)
   unsafe_store!(convert(Ptr{Cint}, (p + offset)), convert(Cint, val))
end
function update_ptr(p::Ptr, offset::Integer, val::mjtNum)
   unsafe_store!(convert(Ptr{mjtNum}, (p + offset)), val)
end
function update_ptr(p::Ptr, offset::Integer, val::SVector)
   T = typeof(val[1])
   for i=1:length(val)
      unsafe_store!(convert(Ptr{T}, p+offset+(i-1)*sizeof(T)),
                    val[i])
   end
end

# mutate mujoco struct fields through the julia version of model and data
#function set(d::jlData, field::Symbol, val::Union{Integer, mjtNum})
#   f_off, f_type = dinfo[field]
#   update_ptr(d.d, f_off, convert(f_type, val))
#end
#function set(m::jlModel, field::Symbol, val::Union{Integer, mjtNum})
#   f_off, f_type = minfo[field]
#   update_ptr(m.m, f_off, convert(f_type, val))
#end
function set(p::Ptr{T}, field::Symbol, val::Union{Integer, mjtNum, SVector}) where T
   f_off, f_type = mjstructs[T][field]
   update_ptr(p, f_off, convert(f_type, val))
end
function set(p::Ptr{T}, field::Symbol, val, i::Int) where T # write to element in SVector
   f_off, f_type = mjstructs[T][field]
   ET = eltype(f_type)
   @assert f_type <: SVector
   #@assert typeof(val) == ET
   v = convert(ET, val) # use this as a check
   @assert i <= length(f_type) && i >= 1
   unsafe_store!(convert(Ptr{ET}, (p+f_off+(i-1)*sizeof(ET))), v)
end
function set(p::Ptr{T}, field::Symbol, val, i::Int, j::Int) where T # write to element in SVector
   f_off, f_type = mjstructs[T][field]
   @assert f_type <: SVector && eltype(f_type) <: SVector
   ET = eltype(eltype(f_type))
   v = convert(ET, val) # use this as a check
   #r = size(f_type)[1]
   c = sizeof(eltype(f_type))[1]
   #@assert i <= r && i >= 1
   #@assert j <= c && i >= 1
   #idx = (i-1) + (j-1)*r
   idx = p+f_off + (i-1)*c + (j-1)*sizeof(ET)
   unsafe_store!(convert(Ptr{ET}, idx), v)
end

# set struct within model struct
function set(m::jlModel, fstruct::Symbol, field::Symbol, val::Union{Integer, mjtNum, SVector})
   s_off, s_type = minfo[fstruct]
   @assert s_type in (mjOption, mjVisual, mjStatistic)

   f_off, f_type = mjstructs[s_type][field]
   update_ptr(m.m, s_off+f_off, convert(f_type, val))
end
function set(p::Ptr{T}, fstruct::Symbol, field::Symbol, val::Union{Integer, mjtNum}) where T
   s_off, s_type = mjstructs[T][fstruct]
   f_off, f_type = mjstructs[s_type][field]
   update_ptr(p, f_off, convert(f_type, val))
end

#################################### Name Wrappers

function name2idx(m::jlModel, num::Integer, names::Vector{Cint})
   sname = String(copy(m.names))
   idx = names[1] + 1
   split_names = split(sname[idx:end], '\0', limit=(num+1))[1:num]
   d = Dict{Symbol, Integer}(Symbol(split_names[i]) => i for i=1:num)
   return d
end

function name2range(m::jlModel, names::Vector{Cint}, addresses::Vector{Cint})
   return name2range(m, names, addresses, ones(Cint, length(addresses)))
end

function name2range(m::jlModel, num::Integer,
                    names::Vector{Cint}, addresses::Vector{Cint}, dims::Vector{Cint})
   sname = String(copy(m.names))
   idx = names[1] + 1
   split_names = split(sname[idx:end], '\0', limit=(num+1))[1:num]
   d = Dict{Symbol, AbstractRange}(Symbol(split_names[i]) => (addresses[i]+1):(addresses[i]+dims[i]) for i=1:num)
   return d
end


text:
const STATIC_FILE = joinpath(@__DIR__, "fixtures", "static.txt")
const STATIC_DOCUMENTER = [
    PackageSpec(; name="DocStringExtensions", version=v"0.8.1"),
    PackageSpec(; name="Documenter", version=v"0.24.2"),
    PackageSpec(; name="JSON", version=v"0.21.0"),
    PackageSpec(; name="Parsers", version=v"0.3.10"),
]

PT.user_view(::Citation, ::Template, ::AbstractString) = Dict("MONTH" => 8, "YEAR" => 2019)
PT.user_view(::License, ::Template, ::AbstractString) = Dict("YEAR" => 2019)

function pin_documenter(project_dir::AbstractString)
    @suppress PT.with_project(project_dir) do
        foreach(Pkg.add, STATIC_DOCUMENTER)
    end
    toml = joinpath(project_dir, "Project.toml")
    project = TOML.parsefile(toml)
    filter!(p -> p.first == "Documenter", project["deps"])
    open(io -> TOML.print(io, project), toml, "w")
end

function test_all(pkg::AbstractString; kwargs...)
    t = tpl(; kwargs...)
    with_pkg(t, pkg) do pkg
        pkg_dir = joinpath(t.dir, pkg)
        PT.hasplugin(t, Documenter) && pin_documenter(joinpath(pkg_dir, "docs"))
        foreach(readlines(`git -C $pkg_dir ls-files`)) do f
            reference = joinpath(@__DIR__, "fixtures", pkg, f)
            observed = read(joinpath(pkg_dir, f), String)
            @test_reference reference observed
        end
    end
end

@testset "Reference tests" begin
    @testset "Default package" begin
        test_all("Basic"; authors=USER)
    end

    @testset "All plugins" begin
        test_all("AllPlugins"; authors=USER, plugins=[
            AppVeyor(), CirrusCI(), Citation(), Codecov(), CompatHelper(), Coveralls(),
            Develop(), Documenter(), DroneCI(), GitHubActions(), GitLabCI(), TravisCI(),
        ])
    end

    @testset "Documenter (TravisCI)" begin
        test_all("DocumenterTravis"; authors=USER, plugins=[
            Documenter{TravisCI}(), TravisCI(),
        ])
    end

    @testset "Documenter (GitHubActions)" begin
        test_all("DocumenterGitHubActions"; authors=USER, plugins=[
            Documenter{GitHubActions}(), GitHubActions(),
        ])
    end

    @testset "Wacky options" begin
        test_all("WackyOptions"; authors=USER, julia=v"1.2", plugins=[
            AppVeyor(; x86=true, coverage=true, extra_versions=[v"1.1"]),
            CirrusCI(; image="freebsd-123", coverage=false, extra_versions=["1.1"]),
            Citation(; readme=true),
            Codecov(; file=STATIC_FILE),
            CompatHelper(; cron="0 0 */3 * *"),
            Coveralls(; file=STATIC_FILE),
            Documenter{GitLabCI}(
                assets=[STATIC_FILE],
                makedocs_kwargs=Dict(:foo => "bar", :bar => "baz"),
                canonical_url=(_t, _pkg) -> "http://example.com",
            ),
            DroneCI(; amd64=false, arm=true, arm64=true, extra_versions=["1.1"]),
            Git(; ignore=["a", "b", "c"], manifest=true),
            GitHubActions(; x86=true, linux=false, coverage=false),
            GitLabCI(; coverage=false, extra_versions=[v"0.6"]),
            License(; name="ISC"),
            ProjectFile(; version=v"1"),
            Readme(; inline_badges=true),
            TagBot(;
                cron="0 0 */3 * *",
                token=Secret("MYTOKEN"),
                ssh=Secret("SSHKEY"),
                ssh_password=Secret("SSHPASS"),
                changelog="Line 1\nLine 2\n\nLine 4",
                changelog_ignore=["foo", "bar"],
                gpg=Secret("GPGKEY"),
                gpg_password=Secret("GPGPASS"),
                registry="Foo/Bar",
                branches=false,
                dispatch=true,
                dispatch_delay=20,
            ),
            Tests(; project=true),
            TravisCI(;
                coverage=false,
                windows=false,
                x86=true,
                arm64=true,
                extra_versions=["1.1"],
            ),
        ])
    end
end


text:
# @generated function recursofteq(x::T, y::T) where T
#   fieldnames(T)
#   all((softeq(x[]) )
# end

function recursofteq(x::T, y::T) where T
  allₛ([softeq(getfield(x, fn), getfield(y, fn)) for fn in fieldnames(T)])
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
include("orderbook.jl")
include("blotter.jl")
include("trades.jl")
include("portfolio.jl")
include("account.jl")
include("utilities.jl")


text: Dict{String, Any}("mixed" => Any[Any[Dict{String, Any}("value" => "1", "type" => "integer"), Dict{String, Any}("value" => "2", "type" => "integer")], Any[Dict{String, Any}("value" => "a", "type" => "string"), Dict{String, Any}("value" => "b", "type" => "string")], Any[Dict{String, Any}("value" => "1.1", "type" => "float"), Dict{String, Any}("value" => "2.1", "type" => "float")]])


text:
# Methods from Moukalled et al. FVM - OpenFOAM, Matlab

mutable struct Cell
    faceIndices::Vector{Int64}
    pointIndices::Vector{Int64}
end

######################### Mesh/Cell Geometry ###########################
function triangleCentroid(points)
    center = [ 0.0, 0.0, 0.0 ]
    nPts = size(points, 1)
    for pt in 1:nPts
        center = center .+ points[pt]
    end
    center /= nPts

    return center
end

#Alternate name for same function
geometricCenter = triangleCentroid

function triangleArea(points::Array{Array{Float64, 1}})
    side1 = points[2] .- points[1]
    side2 = points[3] .- points[1]
    fAVec = cross(side1, side2) ./ 2
    return fAVec
end

#=
    Calculates face area vector and centroid from the points that make up the face
    Points must be ordered sequentially

    How it works:
        1. Splits face into subtriangles
        2. Area and centroid is computed for each subtriangle
        3. Areas vectors are summed and returned
        4. The centroid returned is obtained from an area-weighted of sum of the subtriangle centroids
=#
function faceAreaCentroid(points::Array{Array{Float64, 1}})
    gC = geometricCenter(points)
    nPts = size(points, 1)

    fAVec = [ 0.0, 0.0, 0.0 ]
    centroid = [ 0.0, 0.0, 0.0 ]

    for i in 1:nPts
        if i < nPts
            subTriPts = [ gC, points[i], points[i+1] ]
        else
            subTriPts = [ gC, points[i], points[1] ]
        end

        triCentroid = triangleCentroid(subTriPts)
        subFAVec = triangleArea(subTriPts)

        fAVec += subFAVec
        centroid += triCentroid .* mag(subFAVec)
    end

    centroid /= mag(fAVec)

    return fAVec, centroid
end

#=
    Calculates cell volume (scalar) and centroid (vector) from the points and face area vectors (fAVecs) that make up the cell
        fAVecs can be computed using the faceAreaCentroids function

    How it works:
        1. Splits cell into polygonal pyramids, each incorporating a single face and the geometric center of the cell
        2. Computes volume and centroid of each sub-pyramid
        3. Resulting volume is sum of sub-pyramid volumes, centroid is the volume-weighted sum of sub-pyramid centroids
=#
function cellVolCentroid(points, fAVecs, faceCentroids)
    gC = geometricCenter(points)
    nFaces = size(fAVecs,1)

    vol = 0.0
    centroid = [ 0.0, 0.0, 0.0 ]

    for f in 1:nFaces
        cellCenterVec = faceCentroids[f] .- gC
        subPyrVol = abs(sum(fAVecs[f] .* cellCenterVec) / 3)
        subPyrCentroid = 0.75.*faceCentroids[f] .+ 0.25.*gC

        vol += subPyrVol
        centroid += subPyrCentroid .* subPyrVol
    end

    centroid /= vol

    return vol, centroid
end

# For each face, calculates the vectors from its owner and neighbour cell centers to its own center
function cellCentroidToFaceVec(faceCentroids::Array{Array{Float64, 1}}, cellCentroids::Array{Array{Float64, 1}})
    #TODO: Understand ordering to make sure we're accepting faces in the order we need to pass them back in??
    nFaces = size(faceCentroids, 1)

    cellToFaceVec = zeros(nFaces, 3)

    for f in 1:nFaces
        cellToFaceVec[f,:] = faceCentroids[f] - cellCentroids
    end

    return cellToFaceVec
end

######################### Utility Functions ###########################
# Returns basic info about the mesh: nCells, nFaces, nBoundaries, nBdryFaces
function unstructuredMeshInfo(mesh::Mesh)
    nCells = size(mesh.cells, 1)
    nFaces = size(mesh.faces, 1)
    nBoundaries = size(mesh.boundaryFaces, 1)

    # Count boundary faces
    nBdryFaces = 0
    for bdry in 1:nBoundaries
        nBdryFaces += size(mesh.boundaryFaces[bdry], 1)
    end

    return nCells, nFaces, nBoundaries, nBdryFaces
end

# Checks whether a string represents a number (using a regular expression)
function isNumber(str)
    re = r"^[+-]?([0-9]+([.][0-9]*)?|[.][0-9]+)$"
    return occursin(re, str)
end

######################### Parse OpenFOAM Meshes ###########################
#=
    In the points, faces, owner, and neighbour files, the beginning of useful information in the file looks like this:

        ...
        object      faces;
        }
        // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

        5135
        (
        4(79 1516 1380 37)
        4(1381 1422 1249 39)
        ...

    The 5135 indicates how many lines of information follow it.
    This function would return 5135, and the line number of the first piece of information in the file ((line number of 5135) + 2)
=#
function OFFile_FindNItems(fileLines)
    nLines = size(fileLines, 1)
    lineCounter = 1
    itemCount = 0
    startLine = 0
    for line in fileLines
        if isNumber(strip(line))
            itemCount = parse(Int64, line)
            startLine = lineCounter+2
            break
        end
        lineCounter += 1
    end
    return startLine, itemCount
end

#=
    Takes the file path to an OpenFOAM mesh 'points' file, returns a 2D array of point coordinates:
            X   Y   Z
        P1  x1  y1  z1
        P2  x2  y2  z2
        ...
=#
function readOFPointsFile(filePath)
    f = open(filePath)
    pointsLines = readlines(f)
    startLine, pCount = OFFile_FindNItems(pointsLines)

    points = zeros(pCount, 3)
    for line in startLine:(startLine+pCount-1)
        pLine = pointsLines[line]

        bracketsRemoved = pLine[2:end-1]
        coords = split(bracketsRemoved)
        for d in 1:3
            points[line+1-startLine, d] = parse(Float64, coords[d])
        end
    end

    return points
end

#=
    Takes the file path to an OpenFOAM mesh 'faces' file, returns a 1D array of 1D arrays.
    Each entry in the upper level array represents a single face
    Each face-array contains the indices of the points that make up that face
    Note that all indices from the file are incremented by 1, to switch from OpenFOAM's 0-based indexing to Julia's 1-based indexing

    Example return value:
        [ [face1_pt1Index, face1_pt2Index, face1_pt3Index], [face2_pt1Index, face2_pt2Index, face2_pt3Index, face2_pt4Index], ... ]
=#
function readOFFacesFile(filePath)
    f = open(filePath)
    facesLines = readlines(f)
    startLine, fCount = OFFile_FindNItems(facesLines)

    faces = []
    for line in startLine:(startLine+fCount-1)
        fLine = facesLines[line]
        facePts = []

        bracketIndex = findfirst("(", fLine)[1]
        bracketsRemoved = fLine[bracketIndex+1:end-1]
        ptNumbers = split(bracketsRemoved)
        nPts = length(ptNumbers)
        for p in 1:nPts
            # Add one to convert to 1-based indexing
            push!(facePts, parse(Int64, ptNumbers[p])+1)
        end
        push!(faces, facePts)
    end

    return faces
end

#=
    Takes the file path to an OpenFOAM 'owner' file, returns a 1D array of owner cell indices
    Example return value:
        [ face1_ownerCellIndex, face2_ownerCellIndex, ...  ]
=#
function readOFOwnerFile(filePath)
    f = open(filePath)
    ownerLines = readlines(f)
    startLine, oCount = OFFile_FindNItems(ownerLines)

    faceOwnerCells = zeros(Int64, oCount)
    for line in startLine:(startLine+oCount-1)
        # Add one to convert to 1-based indexing
        faceOwnerCells[line-startLine+1] = parse(Int64, ownerLines[line])+1
    end

    return faceOwnerCells
end

#=
    Takes the file path to an OpenFOAM 'neighbour' file, returns a 1D array of neighbour cell indices
    Note that all indices from the file are incremented by 1, to switch from OpenFOAM's 0-based indexing to Julia's 1-based indexing

    Example return value:
        [ face1_neighbourCellIndex, face2_neighbourCellIndex, ...  ]
=#
function readOFNeighbourFile(filePath)
    f = open(filePath)
    neighbourLines = readlines(f)
    startLine, nCount = OFFile_FindNItems(neighbourLines)

    faceNeighbourCells = zeros(Int64, nCount)
    for line in startLine:(startLine+nCount-1)
        # Add one to convert to 1-based indexing
        faceNeighbourCells[line-startLine+1] = parse(Int64, neighbourLines[line])+1
    end

    return faceNeighbourCells
end

#=
    Returns index of first line containing str
    Otherwise returns -1
=#
function findInLines(str, lines, startLine)
    nLines = size(lines, 1)
    for i in startLine:nLines
        if occursin(str, lines[i])
            return i
        end
    end
    return -1
end

#=
    Takes the file path to an OpenFOAM 'boundary' file, returns 3 arrays:
    boundaryNames:          Array of boundary names (strings)
    boundaryNumFaces:       Array containing the number of faces in each boundary
    boundaryStartFaces:     Array containing the index of the first face in the boundary (boundaries always occupy a continuous sequence of face indices, and are usually numbered last)
=#
function readOFBoundaryFile(filePath)
    f = open(filePath)
    bLines = readlines(f)
    startLine, bCount = OFFile_FindNItems(bLines)
    nLines = size(bLines, 1)

    boundaryNames = Array{String, 1}(undef, bCount)
    boundaryStartFaces = Array{Int64, 1}(undef, bCount)
    boundaryNumFaces = Array{Int64, 1}(undef, bCount)
    for i in 1:bCount
        bNameLine = findInLines("{", bLines, startLine)-1
        boundaryNames[i] = strip(bLines[bNameLine])

        bNFacesLine = findInLines("nFaces", bLines, startLine)
        boundaryNumFaces[i] = parse(Int64, split(bLines[bNFacesLine])[2][1:end-1])

        bStartFaceLine = findInLines("startFace", bLines, startLine)
        boundaryStartFaces[i] = parse(Int64, split(bLines[bStartFaceLine])[2][1:end-1])+1

        startLine = findInLines("}", bLines, startLine)+1
    end

    return boundaryNames, boundaryNumFaces, boundaryStartFaces
end

#=
    Input: Path to an OpenFOAM mesh FOLDER.
    Output: Calls readOFPoints/Faces/Owner/NeighbourFile and returns all of their results (reads all OpenFOAM's mesh file data into arrays)
=#
function readOpenFOAMMesh(polyMeshPath)
    println("Reading mesh: $polyMeshPath")
    pointsFilePath = "$polyMeshPath/points"
    points = readOFPointsFile(pointsFilePath)
    facesFilePath = "$polyMeshPath/faces"
    faces = readOFFacesFile(facesFilePath)
    ownerFilePath = "$polyMeshPath/owner"
    owner = readOFOwnerFile(ownerFilePath)
    neighbourFilePath = "$polyMeshPath/neighbour"
    neighbour = readOFNeighbourFile(neighbourFilePath)
    boundaryFilePath = "$polyMeshPath/boundary"
    boundaryNames, boundaryNumFaces, boundaryStartFaces = readOFBoundaryFile(boundaryFilePath)

    return points, faces, owner, neighbour, boundaryNames, boundaryNumFaces, boundaryStartFaces
end

#=
    Function used to find all the points in each cell.
    Not required for CFD, which is face-based, but required for cell-based .vtk file output.

    Input: Path to an OpenFOAM mesh FOLDER.
    Returns:
        points:         Array of points obtained from readOFPoints()
        cellPtIndices:  Array of arrays, where each subarrary represents a cell, and each entry in a cell's array is the index of one of its points
=#
function OpenFOAMMesh_findCellPts(polyMeshPath)
    pointLocations, pointIndicesByFace, faceOwnerCellIndices, faceNeighborCellIndices, boundaryNames, boundaryNumFaces, boundaryStartFaces = readOpenFOAMMesh(polyMeshPath)
    nCells = maximum(faceOwnerCellIndices)
    nFaces = size(pointIndicesByFace, 1)
    nBoundaries = size(boundaryNames, 1)

    # Output
    cells = Array{Cell, 1}(undef, nCells)
    for i in eachindex(cells)
        cells[i] = Cell(Vector{Int64}(undef, 0), Vector{Int64}(undef, 0))
    end

    ### Populate the cells array ###
    function addCellFaceIndices(adjacentCellList)
        nFaces = length(adjacentCellList)
        for f in 1:nFaces
            ownerCellIndex = adjacentCellList[f]
            push!(cells[ownerCellIndex].faceIndices, f)
        end
    end

    # Start by a just adding face indices, want to order the points appropriately before adding them
    addCellFaceIndices(faceOwnerCellIndices)
    addCellFaceIndices(faceNeighborCellIndices)

    function addAllNewPoints!(cell::Cell, faceIndex)
        for pointIndex in pointIndicesByFace[faceIndex]
            if !any(x->x==pointIndex, cell.pointIndices)
                push!(cell.pointIndices, pointIndex)
            end
        end
    end

    function disjoint(f1Points, f2Points)
        allPoints = vcat(f1Points, f2Points)
        pointSet = Set(allPoints)
        return length(allPoints) == length(pointSet)
    end

    function intersection(points1, points2)
        return intersect(Set(points1), Set(points2))
    end

    function index(element, array)
        for i in eachindex(array)
            if element == array[i]
                return i
            end
        end
        return -1
    end

    function addEdges!(cellPoints, facePoints, endFacePoints, endFacePointsSet, unusedFaces, pointOffset=4)
        oppositeFaceIndex = -1
        for i in 1:length(unusedFaces)
            fi = unusedFaces[i]
            fiPoints = pointIndicesByFace[fi]
            commonPoints = intersection(fiPoints, facePoints)

            if length(commonPoints) == 2
                # These two points form an edge connecting faces one and two
                points = collect(commonPoints)

                if points[1] in endFacePointsSet
                    p1 = points[1]
                    p2 = points[2]
                else
                    p1 = points[2]
                    p2 = points[1]
                end

                position = index(p1, endFacePoints) # Find position over the point on the face one side
                if position == -1
                    throw(ErrorException("Point $p1 not found in $endFacePoints"))
                end
                cellPoints[ position + pointOffset ] = p2 # The point on the face two side goes in the matching spot
            else
                oppositeFaceIndex = i # This is the face opposite to the one containing facePoints
            end
        end
        return oppositeFaceIndex
    end

    function populatePointIndices_Tet!(cell::Cell)
        # For a tetrahedron, the order of points is unimportant
        # And all points from an arbitrary face
        addAllNewPoints!(cell, cell.faceIndices[1])

        # Choose another face, it will contain the fourth point we need to complete the tetrahedron
        addAllNewPoints!(cell, cell.faceIndices[2])
    end

    function populatePointIndices_Pyramid!(cell::Cell)
        # Pyramids need to have the points that make up their square base ordered 1-4, with the tip of the pyramid coming last
        # Find the quad face
        for faceIndex in cell.faceIndices
            if length(pointIndicesByFace[faceIndex]) == 4
                addAllNewPoints!(cell, faceIndex)

                # There is only one square face, and all other faces contain the pyramid tip
                # Find any other face and add its points to make the pyramid tip the final point
                if faceIndex != cell.faceIndices[1]
                    addAllNewPoints!(cell, cell.faceIndices[1])
                else
                    addAllNewPoints!(cell, cell.faceIndices[2])
                end

                # All done
                break
            end
        end
    end

    function populatePointIndices_Wedge!(cell::Cell)
        # Wedges need the triangular faces numbered 1-3 and 4-6 respectively, where 1 is aligned with 4, 2 with 5, and 3 with 6
        # Using similar strategy it to that used for hex cells: Identify two end faces, then use the faces that connect them to ensure their points are aligned

        # First identify the two triangular faces, keep track of one and get rid of the other
        unusedFaces = deepcopy(cell.faceIndices)
        t1 = -1 # Triangular face one, we will base the ordering off this one

        for j in 1:2
            for i in eachindex(unusedFaces)
                faceIndex = unusedFaces[i]

                if length(pointIndicesByFace[faceIndex]) == 3
                    # Found a triangular face
                    if t1 == -1
                        t1 = faceIndex
                    end
                    deleteat!(unusedFaces, i)
                    break
                end
            end
        end

        if length(unusedFaces) != 3
            throw(ErrorException("Expected three quadrilateral faces remaining after removal of two triangular end faces, got $(length(unusedFaces))"))
        end

        t1Points = pointIndicesByFace[t1]
        t1PointsSet = Set(t1Points)
        cell.pointIndices = vcat(t1Points, [0, 0, 0]) # Still have to determine proper ordering of points in the other triangular end face, 0's will be replaced below

        # Now all the remaining faces (3) will be quadrilaterals connecting the two triangular end faces
        q1 = pop!(unusedFaces) # Select one at random
        q1Points = pointIndicesByFace[q1]
        addEdges!(cell.pointIndices, q1Points, t1Points, t1PointsSet, unusedFaces, 3) # Use it to populate 2/3 of the unknown points

        q2 = pop!(unusedFaces) # Select another at random
        q2Points = pointIndicesByFace[q2]
        addEdges!(cell.pointIndices, q2Points, t1Points, t1PointsSet, unusedFaces, 3) # Use it to populate the last unknown point
    end

    function populatePointIndices_Hex!(cell::Cell)
        # Need to order point such that the points that make up a single quadrilateral face are ordered 1-4
            # The face opposite to it must have its points ordered 5-8, with point 1 connected to/aligned with point 5, 2 with 6, etc...
            # We will accomplish this by first picking a starting face at random, this face's points will be numbered 1-4
            # Then we will identify its opposite face and remove it
            # Then all the remaining faces form the connections between the first face we picked and its opposite Face
            # Pairs of these faces will contain 0 or 2 points in common
                # When zero, they are opposite each other
                # When two, they are adjacent and the two points they have in common form an edge connecting the two end faces
                    # This indicates that these two points should be aligned, and is used to properly order the remaining points

        # Start by choosing an arbitrary face
        unusedFaces = deepcopy(cell.faceIndices)
        f1 = pop!(unusedFaces)
        f1Points = pointIndicesByFace[f1]
        f1PointsSet = Set(f1Points)

        # Find the opposite face (no points in common) and remove it
        for i in 1:length(unusedFaces)
            fi = unusedFaces[i]
            fiPoints = pointIndicesByFace[fi]

            if disjoint(f1Points, fiPoints)
                deleteat!(unusedFaces, i)
                break
            end
        end

        ### Now find the correct orientation
            # Point f1_i needs to be aligned (spatially) with Point f2_i
            # We can check for this by using one of the other faces as a guide, since it forms part of the connection between f1 and f2
        cell.pointIndices = vcat(f1Points, [0, 0, 0, 0])

        # Pick an arbitrary other face
        f3 = pop!(unusedFaces)
        f3Points = pointIndicesByFace[f3]
        lastFace = addEdges!(cell.pointIndices, f3Points, f1Points, f1PointsSet, unusedFaces)

        # Add points from the edges of the face opposite f3
        lastFacePoints = pointIndicesByFace[unusedFaces[lastFace]]
        deleteat!(unusedFaces, lastFace)
        noResult = addEdges!(cell.pointIndices, lastFacePoints, f1Points, f1PointsSet, unusedFaces)

        if noResult != -1
            throw(ErrorException("Failure to appropriately order hexahedral cell points for writing to .vtk"))
        end
    end

    # Now use the face indices to gather and order each cell's points appropriately
        # Correct point ordering is important for .vtk output: (See figure 2) https://vtk.org/wp-content/uploads/2015/04/file-formats.pdf
    for cell in cells
        nFaces = length(cell.faceIndices)
        if nFaces == 4
            populatePointIndices_Tet!(cell)
        elseif nFaces == 5
            quadFaceCount = 0

            for f in cell.faceIndices
                if length(pointIndicesByFace[f]) == 4
                    quadFaceCount += 1
                end
            end

            if quadFaceCount == 1
                populatePointIndices_Pyramid!(cell)
            elseif quadFaceCount == 3
                populatePointIndices_Wedge!(cell)
            else
                throw(ErrorException("Unrecognized cell type, cell: $cell. Expecting only hex, tet, wedge, or pyramid cells for vtk output."))
            end
        elseif nFaces == 6
            populatePointIndices_Hex!(cell)
        end
    end

    return pointLocations, cells
end

#=
    Reads an OpenFOAM mesh and returns a Mesh object suitable for calculations in FvCFD.jl.

    Mesh defined in dataStructures.jl, documented in dataStructureDefintions.md
=#
function OpenFOAMMesh(polyMeshPath)
    points, OFfaces, owner, neighbour, boundaryNames, boundaryNumFaces, boundaryStartFaces = readOpenFOAMMesh(polyMeshPath)
    nCells = maximum(owner)
    nFaces = size(OFfaces, 1)
    nBoundaries = size(boundaryNames, 1)

    cells = Array{Array{Int64, 1}, 1}(undef, nCells)                # Indices of faces that make up a cell
    cVols = zeros(nCells)                                           # Cell Volumes
    cCenters = Array{Array{Float64, 1}, 1}(undef, nCells)           # Cell Centroids
    cellSizes = Matrix{Float64}(undef, nCells, 3)                   # Cell sizes (x,y,z) direction

    faces = Array{Array{Int64, 1}, 1}(undef, nFaces)                # Indices of cells adjacent to each face always owner first, then neighbour
    fAVecs = Array{Array{Float64, 1}, 1}(undef, nFaces)             # Face Area Vectors
    fCenters = Array{Array{Float64, 1}, 1}(undef, nFaces)           # Face Centroids
    boundaryFaces = Array{Array{Int64, 1}, 1}(undef, nBoundaries)   # Indices of faces that make up each boundary

    # Calculate face area vectors and face centroids from list of points in each face
    for f in 1:nFaces
        fPts = [ points[pt,:] for pt in OFfaces[f] ]
        fAVecs[f], fCenters[f] = faceAreaCentroid(fPts)
    end
    #fAVecs and fCenters now complete

    # Build 'cells' and 'faces' arrays from info from 'owner' file
    nOwners = size(owner, 1)
    for f in 1:nOwners
        ownerCell = owner[f]
        if isassigned(cells, ownerCell)
            push!(cells[ownerCell], f)
        else
            cells[ownerCell] = [f,]
        end

        faces[f] = [ownerCell, -1]
    end

    # Finish building 'cells' and 'faces' arrays, adding info from 'neighbour' file
    nNeighbours = size(neighbour, 1)
    for f in 1:nNeighbours
        neighbourCell = neighbour[f]
        if isassigned(cells, neighbourCell)
            push!(cells[neighbourCell], f)
        else
            cells[neighbourCell] = [f,]
        end

        faces[f][2] = neighbourCell
    end
    # fAVecs, fCenters, faces, cells now complete

    # For each cell, use the face area vectors and face centroids for all of its faces, to calculate it's volume and centroid location
    for c in 1:nCells
        pts = []
        for f in cells[c]
            for pt in OFfaces[f]
                if !(points[pt,:] in pts)
                    push!(pts, points[pt,:])
                end
            end
        end
        fCs = [ fCenters[f] for f in cells[c] ]
        cell_fAVecs = [ fAVecs[f] for f in cells[c] ]

        cVols[c], cCenters[c] = cellVolCentroid(pts, cell_fAVecs, fCs)
    end
    # fAVecs, fCenters, faces, cells, cVols, cCenters now complete

    # Create boundaryFaces array
    for b in 1:nBoundaries
        startF = boundaryStartFaces[b]
        endF = startF + boundaryNumFaces[b] - 1
        boundaryFaces[b] = Array(startF:endF)
    end

    # Compute cell sizes in the x, y, z directions
    # In the past, cell sizes were used for a shitty CFL calculation, may be able to get rid of them now
    maxCoords = [ -1000000.0, -1000000.0, -1000000.0 ]
    minCoords = [ 1000000.0, 1000000.0, 1000000.0 ]
    for c in 1:nCells
        fill!(maxCoords, -10000000)
        fill!(minCoords, 10000000)

        pts = Vector{Vector{Float64}}()

        # Add points
        for f in cells[c]
            for pt in OFfaces[f]
                if !any(x->x==pt, pts)
                    push!(pts, points[pt, :])
                end
            end
        end

        for pt in pts
            for d in 1:3
                maxCoords[d] = max(maxCoords[d], pt[d])
                minCoords[d] = min(minCoords[d], pt[d])
            end
        end

        for d in 1:3
            cellSizes[c,d] = maxCoords[d] - minCoords[d]
        end
    end

    return Mesh(cells, cVols, cCenters, cellSizes, faces, fAVecs, fCenters, boundaryFaces)
end

max_stars_repo_name: henrystoldt/fvCFD


text:
using MPIReco

@testset "multi-patch in-memory reconstruction" begin
  bSF = MultiMPIFile([joinpath(datadir, "calibrations", "12.mdf")])
  dirs = ["1.mdf", "2.mdf", "3.mdf", "4.mdf"]
  b = MultiMPIFile(joinpath.(datadir, "measurements", "20211226_203916_MultiPatch", dirs))
  names = (:color, :x, :y, :z, :time)
  values1 = (1:1,
	    -27.5u"mm":1.25u"mm":27.5u"mm",
	    -27.5u"mm":1.25u"mm":27.5u"mm",
	    0.0u"mm":1.0u"mm":0.0u"mm",
	    0.0u"ms":0.6528u"ms":0.0u"ms")
  values2 = (values1[1], -19.375u"mm":1.25u"mm":19.375u"mm",
	     -19.375u"mm":1.25u"mm":19.375u"mm", values1[4:5]...)

  # basic multi-patch reconstruction
  c1 = reconstruction(bSF, b;
			    SNRThresh=5, frames=1, minFreq=80e3,
			    recChannels=1:2,iterations=1,
			    spectralLeakageCorrection=false)
  @test axisnames(c1) == names
  @test axisvalues(c1) == values1
  exportImage(joinpath(imgdir, "MultiPatch1.png"), arraydata(c1[1,:,:,1,1]))

  # TODO test description
  c2 = reconstruction(bSF, b;
			    SNRThresh=5, frames=1, minFreq=80e3,
			    recChannels=1:2,iterations=1, roundPatches=true,
			    spectralLeakageCorrection=false)
  @test axisnames(c2) == names
  @test axisvalues(c2) == values1
  exportImage(joinpath(imgdir, "MultiPatch2.png"), arraydata(c2[1,:,:,1,1]))

  # multi-patch reconstruction using multiple system matrices
  dirs = ["8.mdf", "9.mdf", "10.mdf", "11.mdf"]
  bSFs = MultiMPIFile(joinpath.(datadir, "calibrations", dirs))
  c3 = reconstruction(bSFs, b;
			    SNRThresh=5, frames=1, minFreq=80e3,
			    recChannels=1:2,iterations=1,
			    spectralLeakageCorrection=false)
  @test axisnames(c3) == names
  @test axisvalues(c3) == values2
  exportImage(joinpath(imgdir, "MultiPatch3.png"), arraydata(c3[1,:,:,1,1]))

  # flexible multi-patch reconstruction
  dirs = ["8.mdf", "9.mdf", "10.mdf", "11.mdf"]
  bSFs = MultiMPIFile(joinpath.(datadir, "calibrations", dirs))
  mapping = [1,2,3,4]
  freq = filterFrequencies(bSFs, SNRThresh=5, minFreq=80e3)
  S = [getSF(SF,freq,nothing,"kaczmarz", bgcorrection=false)[1] for SF in bSFs]
  SFGridCenter = zeros(3,4)
  FFPos = zeros(3,4)
  FFPos[:,1] = [-0.008, 0.008, 0.0]
  FFPos[:,2] = [-0.008, -0.008, 0.0]
  FFPos[:,3] = [0.008, 0.008, 0.0]
  FFPos[:,4] = [0.008, -0.008, 0.0]
  c4 = reconstruction(bSFs, b;
			    SNRThresh=5, frames=1, minFreq=80e3,
			    recChannels=1:2,iterations=1,
			    spectralLeakageCorrection=false, mapping=mapping,
			    systemMatrices = S, SFGridCenter=SFGridCenter,
			    FFPos=FFPos, FFPosSF=FFPos)
  @test axisnames(c4) == names
  @test axisvalues(c4) == values2
  exportImage(joinpath(imgdir, "MultiPatch4.png"), arraydata(c4[1,:,:,1,1]))
  # TODO the last test shows odd results
end


text:
using Documenter, CancerImagingArchive

makedocs(;
    modules=[CancerImagingArchive],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    pages=[
        "Introduction" => "index.md",
        "Guide" => Any[
            "Formats" => "guide/formats.md",
            "Queries" => "guide/queries.md",
            "Data Download" => "guide/downloading.md"
        ],
        "Functions" => "functions.md"
    ],
    repo="https://github.com/notZaki/CancerImagingArchive.jl/blob/{commit}{path}#L{line}",
    sitename="CancerImagingArchive.jl"
)

deploydocs(;
    repo="github.com/notZaki/CancerImagingArchive.jl",
)


text:
#=
max_percent_diff.jl
Copyright 2000-9-16, Jeff Fessler, University of Michigan
=#

export max_percent_diff

"""
    d = max_percent_diff(s1, s2, [options])

Compute the "maximum percent difference" between two signals: `s1, s2`.

Default is to normalize by `maximum(abs.(s1))`.

options
- `maxboth::Bool` use max of both arguments to normalize; default `false`
- `normalize::Bool` normalize each before comparing; default `false`

"""
function max_percent_diff(s1, s2 ;
        maxboth::Bool = false, normalize::Bool = false)

    # check to make sure that we have comparable signals
    size(s1) != size(s2) &&
        throw(DomainError("Dimension mismatch $(size(s1)) $(size(s2))"))

    any(isnan.(s1)) && @warn("NaN values in s1 input to max_percent_diff")
    any(isnan.(s2)) && @warn("NaN values in s2 input to max_percent_diff")

    if normalize
        normFactor = sum(abs.(s2)) / sum(abs.(s1))
        isinteractive() && (@info "Normalization factor (matrix 2 / matrix 1): $normFactor \n")
        s1 /= sum(abs.(s1))
        s2 /= sum(abs.(s2))
    end

    if maxboth
        denom = max(maximum(abs.(s1)), maximum(abs.(s2)))
        d = denom == 0 ? 0 : maximum(abs.(s1 - s2)) / denom
    else
        denom = maximum(abs.(s1))
        denom = denom == 0 ? denom : maximum(abs.(s2))
        d = denom == 0 ? 0 : maximum(abs.(s1 - s2)) / denom
    end

    return 100 * d
end


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

ServerDiary.writeDiary()

max_stars_repo_name: Agapanthus/server-diary

max_issues_repo_name: Agapanthus/server-diary

max_forks_repo_name: Agapanthus/server-diary


text:
#= -------------------------------------------------------------------------
# @file testing_transformations.jl
#
# @date 08/07/19 08:34:43
# @author Martin Noblia
# @email mnoblia@disroot.org
#
# @brief
#
# @detail
#
#  Licence:
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or (at
# your option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License

---------------------------------------------------------------------------=#
# testing
using Robotics, LinearAlgebra
# NOTE(elsuizo:2019-08-07): no usar == cuando se comparan Floats
@testset "transformations" begin
   @testset "rotx" begin
      R1 = rotx(deg2rad(30))
      R2 = rotx(deg2rad(30 + 360))
      @test R1 ≈ R2
      @test det(R1) ≈ det(R2) ≈ 1.0
   end
   @testset "roty" begin
      R3 = roty(deg2rad(30))
      R4 = roty(deg2rad(30 + 360))
      @test R3 ≈ R4
      @test det(R3) ≈ det(R4) ≈ 1.0
   end
   @testset "rotz" begin
      R5 = rotz(deg2rad(30))
      R6 = rotz(deg2rad(30 + 360))
      @test R5 ≈ R6
      @test det(R5) ≈ det(R6) ≈ 1.0
   end
end

max_stars_repo_name: elsuizo/Robot

max_issues_repo_name: elsuizo/Robot

max_forks_repo_name: elsuizo/Robot


text:
let
import CuArrays
hasmethod(args...; kwargs...) = Base.hasmethod(args...; kwargs...)
if VERSION < v"1.2"
	# Backport keyword checking support from https://github.com/JuliaLang/julia/pull/30712
   function hasmethod(@nospecialize(f), @nospecialize(t), kwnames::Tuple{Vararg{Symbol}}; world=typemax(UInt))
      Base.hasmethod(f, t, world=world) || return false
      isempty(kwnames) && return true
      m = which(f, t)
      Base.max_world(m) <= world || return false
      kws = Base.kwarg_decl(m, Core.kwftype(typeof(f)))
      for kw in kws
          endswith(String(kw), "...") && return true
      end
      issubset(kwnames, kws)
   end
end

# Backport https://github.com/JuliaGPU/CuArrays.jl/pull/181 if needed
if !hasmethod(CuArrays.CUDNN.conv!, Tuple{CuArrays.CuArray{Float32}, CuArrays.CuArray{Float32}, CuArrays.CuArray{Float32}}, (:flipkernel,))
@eval CuArrays.CUDNN begin
 function conv!(y::CuArray{T}, x::CuArray{T}, w::CuArray{T};
               pad=0, stride=1, flipkernel=0, alpha=1, dilation=1,
               workspace::Union{CuVector, Nothing}=nothing, algo=0) where T<:CUDNNFloat
  if CUDNN_VERSION < 6000
    all(x -> x == 1, dilation) || error("Only dilation = 1 is supported in cuDNN version < 6")
  end
  if workspace === nothing
    workspace_size =
      cudnnGetConvolutionForwardWorkspaceSize(y, x, w, padding=pad, stride=stride, dilation=dilation,
                                              algo=algo, mode=flipkernel)
    workspace = workspace_size != 0 ? conv_workspace(workspace_size) : workspace
  else
    workspace_size = length(workspace[])
  end
  cudnnConvolutionForward(y, x, w, padding=pad, stride=stride, dilation=dilation, mode=flipkernel,
			  alpha=alpha, algo=algo, workspace=workspace, workspace_size=workspace_size)
end

function ∇conv_filter!(dw::CuArray{T}, dy::CuArray{T}, x::CuArray{T}, w::CuArray{T};
                       pad=0, stride=1, flipkernel=0, alpha=1, dilation=1,
                       workspace::Union{CuVector, Nothing}=nothing, algo=0) where T<:CUDNNFloat
  if CUDNN_VERSION < 6000
    all(x -> x == 1, dilation) || error("Only dilation = 1 is supported in cuDNN version < 6")
  end
  if workspace === nothing
    workspace_size =
      cudnnGetConvolutionBackwardFilterWorkspaceSize(dw, x, w, dy, padding=pad, stride=stride,
					             dilation=dilation, algo=algo, mode=flipkernel)
    workspace = workspace_size != 0 ? conv_workspace(workspace_size) : workspace
  else
    workspace_size = length(workspace[])
  end
  cudnnConvolutionBackwardFilter(dw, x, w, dy, padding=pad, stride=stride, dilation=dilation,
				 mode=flipkernel, alpha=alpha, algo=algo, workspace=workspace,
                                 workspace_size=workspace_size)
end

function ∇conv_data!(dx::CuArray{T}, dy::CuArray{T}, x::CuArray{T}, w::CuArray{T};
                     pad=0, stride=1, flipkernel=0, alpha=1, dilation=1,
                     workspace::Union{CuVector, Nothing}=nothing, algo=0) where T<:CUDNNFloat
  if CUDNN_VERSION < 6000
    all(x -> x == 1, dilation) || error("Only dilation = 1 is supported in cuDNN version < 6")
  end
  if workspace === nothing
    workspace_size =
      cudnnGetConvolutionBackwardDataWorkspaceSize(dx, x, w, dy, padding=pad, stride=stride,
                                                   dilation=dilation, algo=algo, mode=flipkernel)
    workspace = workspace_size != 0 ? conv_workspace(workspace_size) : workspace
  else
    workspace_size = length(workspace[])
  end
  cudnnConvolutionBackwardData(dx, x, w, dy, padding=pad, stride=stride, dilation=dilation,
			       mode=flipkernel, alpha=alpha, algo=algo, workspace=workspace,
                               workspace_size=workspace_size)
end
end
end
end

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
export regularize

"derivative of regularizar"
function regularize(W::Array{Float32})
  reg=W # L1:sign(W)
end


text:
using HardTestProblems
using Test

@testset "PMM" begin
    for i in 1:6
        for (uldim, lldim) in zip([2,5,10], [3,5,10])
            # objective function and problem configuration
            F, f, conf = PMM_get_problem(i, uldim = uldim, lldim = lldim)
            cond1 = length(conf[:xmin]) == uldim == length(conf[:xmax])
            cond2 = length(conf[:ymin]) == lldim == length(conf[:ymax])
            cond3 = length(conf[:xbest]) == uldim && length(conf[:ybest]) == lldim
            cond4 = conf[:leader_optimum] ≈ 0 && conf[:follower_optimum] ≈ 0
            !cond4 && display(conf)
            @test cond1 && cond2 && cond3 && cond4
        end
    end
end

@testset "SMD" begin
    for i in 1:12
        for (uldim, lldim) in zip([2,5,10], [3,5,10])
            # objective function and problem configuration
            F, f, conf = SMD_get_problem(i, uldim = uldim, lldim = lldim)
            cond1 = length(conf[:xmin]) == uldim == length(conf[:xmax])
            cond2 = length(conf[:ymin]) == lldim == length(conf[:ymax])
            cond3 = length(conf[:xbest]) == uldim && length(conf[:ybest]) == lldim
            @test cond1 && cond2 && cond3
        end
    end
end

@testset "CEC2020" begin
    for fnum in 1:10
        for D in [5,10,15,20,30,50,100]
            f, conf = get_cec2020_problem(fnum)
            F =  f(rand(77, D))
            fmin = conf[:minimum]
            @test sum(F .< fmin) == 0
        end
    end
end

@testset "RW-MOP-2021" begin
    for i in 1:50
        # objective function and problem configuration
        f, conf = get_RW_MOP_problem(i)

        _, conf_tmp = get_RW_MOP_problem(conf[:function])
        # check get_RW_MOP_problem
        @test conf[:function] == conf_tmp[:function]

        # check dimension
        @test length(conf[:xmin]) == length(conf[:xmax]) == conf[:n]

        # check objectives
        @test length(conf[:nadir]) == conf[:fn]

        fx, gx, hx = f(conf[:xmin])

        # constraint violation sum
        vio = sum(max.(0.0, gx .<= 0)) + sum(max.(0.0, abs.(hx) .- 1e-2 .<= 0))

        # infeasible or feasible and not worst than nadir point.
        @test vio > 0 || (vio == 0 && sum(fx .> conf[:nadir]) == 0)
    end

end


text:
function get_points_spherical_random(n::Int)
    r = ones(n)
    θ = acos.(1-2rand(n))
    ϕ = 2pi*rand(n)

    return [r θ ϕ]
end

function get_points_spherical_evenly(n::Int)
    r = ones(n)
    θ = acos.(1-2linspace(0.0, 1.0, n))
    ϕ = 2pi*linspace(0.0, 1.0, n)

    return [r θ ϕ]
end

function get_points_cartesian_random(n::Int)
    points = get_points_spherical_random(n)

    return [sin.(points[:,2]).*cos.(points[:,3]) sin.(points[:,2]).*sin.(points[:,3]) cos.(points[:,2])]
end

function random_sphere_pack(r::AbstractFloat, R::AbstractFloat, n::Int)

    step = r/10
    x = [ collect(-step : -step : -R); collect(0 : step : R)]
    y = [ collect(-step : -step : -R); collect(0 : step : R)]
    z = [ collect(-step : -step : -R); collect(0 : step : R)]
    len = length(x)
    places = zeros(len, len, len)

    for i = 1:len
        for j = 1:len
            for k = 1:len
                if sqrt(x[i]^2 + y[j]^2 + z[k]^2) < R
                    places[i,j,k] = 1
                end
            end
        end
    end

    points = zeros(3, n)
    count = 0

    while count < n
        vacant_points = find(places)
        if isempty(vacant_points)
            error("there are no more vacant points in the space of the cluster at `count`=$count")
        end

        random_vacant_point = vacant_points[rand(1:length(vacant_points))]
        (I,J,K) = ind2sub(places, random_vacant_point)

        if places[I,J,K] != 1
            display("Something went wrong! places(p(1),p(2),p(3))!=1")
        end

        for i = 1:len
            for j = 1:len
                for k = 1:len
                    if sqrt((x[i]-x[I])^2 + (y[j]-y[J])^2 + (z[k]-z[K])^2) < 2r
                        places[i,j,k] = 0
                    end
                end
            end
        end

        count += 1;
        points[:, count] = [x[I] y[J] z[K]]
    end

    return points
end