text:
"""
    make_moons(;    n_samples::Union{Tuple{Int, Int}, Int} = 100,
                    shuffle = true,
                    noise = nothing,
                    random_state = nothing)::DataFrame
Make two interleaving half circles. Sklearn interface to make_moons.
# Arguments
- `n_samples::Union{Tuple{Int, Int}, Int} = 100`: If int, the total number of points generated. If two-element tuple, number of points in each of two moons.
- `shuffle::Bool = true`: Whether to shuffle the samples.
- `noise::Union{Nothing, Float64} = nothing`: Standard deviation of Gaussian noise added to the data.
- `random_state::Union{Int, Nothing} = nothing`: Determines random number generation for dataset shuffling and noise.
Reference: [link](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_moons.html)
"""
function make_moons(;   n_samples::Union{Tuple{Int, Int}, Int} = 100,
                        shuffle::Bool = true,
                        noise::Union{Nothing, Float64} = nothing,
                        random_state::Union{Int, Nothing} = nothing)::DataFrame

    (features, labels) = datasets.make_moons(   n_samples=n_samples,
                                                shuffle = shuffle,
                                                noise = noise,
                                                random_state = random_state)

    return convert(features, labels)
end

"""
    make_blobs(;    n_samples::Union{Int, Array{Int, 1}} = 100,
                    n_features::Int = 2,
                    centers::Union{Int, Union{Nothing, Array{Float64, 2}}} = nothing,
                    cluster_std::Union{Float64, Array{Float64, 1}} = 1.0,
                    center_box = (-10.0, 10.0),
                    shuffle::Bool = true,
                    random_state::Union{Int, Nothing} = nothing)::DataFrame
Generate isotropic Gaussian blobs for clustering. Sklearn interface to make_blobs.
# Arguments
- `n_samples = 100`: If int, it is the total number of points equally divided among clusters. If array-like, each element of the sequence indicates the number of samples per cluster.
- `n_features = 2`: The number of features for each sample.
- `centers::Union{Int, Union{Nothing, Array{Float64, 2}}} = nothing`: The number of centers to generate, or the fixed center locations. If n_samples is an int and centers is None, 3 centers are generated. If n_samples is array-like, centers must be either None or an array of length equal to the length of n_samples.
- `cluster_std::Union{Float64, Array{Float64, 1}} = 1.0`: The standard deviation of the clusters.
- `center_box::Tuple{Float64, Float64} = (-10.0, 10.0)`: The bounding box for each cluster center when centers are generated at random.
- `shuffle::Bool = true`: Shuffle the samples.
- `random_state::Union{Int, Nothing} = nothing`: Determines random number generation for dataset shuffling and noise.
Reference: [link](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_blobs.html)
"""
function make_blobs(;   n_samples::Union{Int, Array{Int, 1}} = 100,
                        n_features::Int = 2,
                        centers::Union{Int, Union{Nothing, Array{Float64, 2}}} = nothing,
                        cluster_std::Union{Float64, Array{Float64, 1}} = 1.0,
                        center_box::Tuple{Float64, Float64} = (-10.0, 10.0),
                        shuffle::Bool = true,
                        random_state::Union{Int, Nothing} = nothing)::DataFrame

    (features, labels) = datasets.make_blobs(   n_samples = n_samples,
                                                n_features = n_features,
                                                centers = centers,
                                                cluster_std = cluster_std,
                                                center_box = center_box,
                                                shuffle = shuffle,
                                                random_state = random_state,
                                                return_centers = false)

    return convert(features, labels)
end

"""
    make_s_curve(;  n_samples::Int = 100,
                    noise = nothing,
                    random_state = nothing)::DataFrame
Generate an S curve dataset. Sklearn interface to make_s_curve.
# Arguments
- `n_samples::Int = 100`: The number of sample points on the S curve.
- `noise::Union{Nothing, Float64} = nothing`: Standard deviation of Gaussian noise added to the data.
- `random_state::Union{Int, Nothing} = nothing`: Determines random number generation for dataset creation. Pass an int for reproducible output across multiple function calls.
Reference: [link](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_s_curve.html)
"""
function make_s_curve(; n_samples::Int = 100,
                        noise::Float64 = 0.0,
                        random_state::Union{Int, Nothing} = nothing)::DataFrame

    (features, labels) = datasets.make_s_curve( n_samples = n_samples,
                                                noise = noise,
                                                random_state = random_state)

    return convert(features, labels)
end

"""
    function make_circles(; n_samples::Int = 100,
                            shuffle::Bool = true,
                            noise::Float64 = 0.0,
                            random_state::Union{Int, Nothing} = nothing,
                            factor::Float64 = 0.8)::DataFrame
Make a large circle containing a smaller circle in 2d. Sklearn interface to make_circles.
# Arguments
- `n_samples::Union{Int, Tuple{Int, Int}} = 100`: If int, it is the total number of points generated. For odd numbers, the inner circle will have one point more than the outer circle. If two-element tuple, number of points in outer circle and inner circle.
- `shuffle::Bool = true`: Whether to shuffle the samples.
- `noise::Union{Nothing, Float64} = nothing`: Standard deviation of Gaussian noise added to the data.
- `random_state::Union{Int, Nothing} = nothing`: Determines random number generation for dataset shuffling and noise. Pass an int for reproducible output across multiple function calls.
- `factor::Float64 = 0.8`: Scale factor between inner and outer circle.
Reference: [link](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_circles.html)

"""
function make_circles(; n_samples::Union{Int, Tuple{Int, Int}} = 100,
                        shuffle::Bool = true,
                        noise::Union{Nothing, Float64} = nothing,
                        random_state::Union{Int, Nothing} = nothing,
                        factor::Float64 = 0.8)::DataFrame

    (features, labels) = datasets.make_circles( n_samples = n_samples,
                            shuffle = shuffle,
                            noise = noise,
                            random_state = random_state,
                            factor = factor)

    return convert(features, labels)
end

"""
    make_regression(;   n_samples::Int = 100,
                        n_features::Int = 100,
                        n_informative::Int = 10,
                        n_targets::Int = 1,
                        bias::Float64 = 0.0,
                        effective_rank::Union{Int, Nothing} = nothing,
                        tail_strength::Float64 = 0.5,
                        noise::Float64 = 0.0,
                        shuffle::Bool = true,
                        coef::Bool = false,
                        random_state::Union{Int, Nothing}= nothing)
Generate a random regression problem. Sklearn interface to make_regression.
# Arguments
- `n_samples::Int = 100`: The number of samples.
- `n_features::Int = 2`: The number of features.
- `n_informative::Int = 10`: The number of informative features, i.e., the number of features used to build the linear model used to generate the output.
- `n_targets::Int = 1`: The number of regression targets, i.e., the dimension of the y output vector associated with a sample. By default, the output is a scalar.
- `bias::Float = 0.0`: The bias term in the underlying linear model.
- `effective_rank::Union{Int, Nothing} = nothing`: If not `nothing`, the approximate number of singular vectors required to explain most of the input data by linear combinations. Using this kind of singular spectrum in the input allows the generator to reproduce the correlations often observed in practice. If `nothing`, the input set is well conditioned, centered and gaussian with unit variance.
- `tail_strength::Float = 0.5`: The relative importance of the fat noisy tail of the singular values profile if effective_rank is not None.
- `noise::Union{Nothing, Float64} = nothing`: Standard deviation of Gaussian noise added to the data.
- `shuffle::Bool = true`: Shuffle the samples and the features.
- `coef::Bool = false`: If `true`, the coefficients of the underlying linear model are returned.
- `random_state::Union{Int, Nothing} = nothing`: Determines random number generation for dataset shuffling and noise.
Reference: [link](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_regression.html)
"""
function make_regression(;  n_samples::Int = 100,
                            n_features::Int = 100,
                            n_informative::Int = 10,
                            n_targets::Int = 1,
                            bias::Float64 = 0.0,
                            effective_rank::Union{Int, Nothing} = nothing,
                            tail_strength::Float64 = 0.5,
                            noise::Float64 = 0.0,
                            shuffle::Bool = true,
                            coef::Bool = false,
                            random_state::Union{Int, Nothing}= nothing)

    (features, labels) = datasets.make_regression(  n_samples = n_samples,
                                                    n_features = n_features,
                                                    n_informative = n_informative,
                                                    n_targets = n_targets,
                                                    bias = bias,
                                                    effective_rank = effective_rank,
                                                    tail_strength = tail_strength,
                                                    noise = noise,
                                                    shuffle = shuffle,
                                                    coef = coef,
                                                    random_state = random_state)

    return convert(features, labels)
end

"""
    function make_classification(;  n_samples::Int = 100,
                                    n_features::Int = 20,
                                    n_informative::Int = 2,
                                    n_redundant::Int = 2,
                                    n_repeated::Int = 0,
                                    n_classes::Int = 2,
                                    n_clusters_per_class::Int = 2,
                                    weights::Union{Nothing, Array{Float64,1}} = nothing,
                                    flip_y::Float64 = 0.01,
                                    class_sep::Float64 = 1.0,
                                    hypercube::Bool = true,
                                    shift::Union{Nothing, Array{Float64,1}} = 0.0,
                                    scale::Union{Nothing, Array{Float64,1}} = 1.0,
                                    shuffle::Bool = true,
                                    random_state::Union{Int, Nothing} = nothing)
Generate a random n-class classification problem. Sklearn interface to make_classification.
#Arguments
- `n_samples::Int = 100`: The number of samples.
- `n_features::Int = 20`: The total number of features. These comprise `n_informative` informative features, `n_redundant` redundant features, `n_repeated` duplicated features and `n_features-n_informative-n_redundant-n_repeated` useless features drawn at random.
- `n_informative::Int = 2`: The number of informative features. Each class is composed of a number of gaussian clusters each located around the vertices of a hypercube in a subspace of dimension `n_informative`. For each cluster, informative features are drawn independently from N(0, 1) and then randomly linearly combined within each cluster in order to add covariance. The clusters are then placed on the vertices of the hypercube.
- `n_redundant::Int = 2`: The number of redundant features. These features are generated as random linear combinations of the informative features.
- `n_repeated::Int = 0`: The number of duplicated features, drawn randomly from the informative and the redundant features.
- `n_classes::Int = 2`: The number of classes (or labels) of the classification problem.
- `n_clusters_per_class::Int = 2`: The number of clusters per class.
- `weights::Union{Nothing, Array{Float64,1}} = nothing`:
- `flip_y::Float64 = 0.01`: The fraction of samples whose class is assigned randomly. Larger values introduce noise in the labels and make the classification task harder. Note that the default setting flip_y > 0 might lead to less than n_classes in y in some cases.
- `class_sep::Float64 = 1.0`: The factor multiplying the hypercube size. Larger values spread out the clusters/classes and make the classification task easier.
- `hypercube::Bool = true`: If True, the clusters are put on the vertices of a hypercube. If False, the clusters are put on the vertices of a random polytope.
- `shift::Union{Nothing, Array{Float64,1}} = 0.0`: Shift features by the specified value. If None, then features are shifted by a random value drawn in [-class_sep, class_sep].
- `scale::Union{Nothing, Array{Float64,1}} = 1.0`: Multiply features by the specified value. If None, then features are scaled by a random value drawn in [1, 100]. Note that scaling happens after shifting.
- `shuffle::Bool = true`: Shuffle the samples and the features.
- `random_state::Union{Int, Nothing} = nothing`: Determines random number generation for dataset creation. Pass an int for reproducible output across multiple function calls. See Glossary.
Reference: [link](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_classification.html)
"""
function make_classification(;  n_samples::Int = 100,
                                n_features::Int = 20,
                                n_informative::Int = 2,
                                n_redundant::Int = 2,
                                n_repeated::Int = 0,
                                n_classes::Int = 2,
                                n_clusters_per_class::Int = 2,
                                weights::Union{Nothing, Array{Float64,1}} = nothing,
                                flip_y::Float64 = 0.01,
                                class_sep::Float64 = 1.0,
                                hypercube::Bool = true,
                                shift::Union{Nothing, Float64, Array{Float64,1}} = 0.0,
                                scale::Union{Nothing, Float64, Array{Float64,1}} = 1.0,
                                shuffle::Bool = true,
                                random_state::Union{Int, Nothing} = nothing)

    (features, labels) = datasets.make_classification(  n_samples = n_samples,
                                                        n_features = n_features,
                                                        n_informative = n_informative,
                                                        n_redundant = n_redundant,
                                                        n_repeated = n_repeated,
                                                        n_classes = n_classes,
                                                        n_clusters_per_class = n_clusters_per_class,
                                                        weights = weights,
                                                        flip_y = flip_y,
                                                        class_sep = class_sep,
                                                        hypercube = hypercube,
                                                        shift = shift,
                                                        scale = scale,
                                                        shuffle = shuffle,
                                                        random_state = random_state)

    return convert(features, labels)
end

"""
    function make_friedman1(;   n_samples::Int = 100,
                                n_features::Int = 10,
                                noise::Float64 = 0.0,
                                random_state::Union{Int, Nothing} = nothing)::DataFrame
Generate the “Friedman #1” regression problem. Sklearn interface to make_regression.
#Arguments
- `n_samples::Int = 100`: The number of samples.
- `n_features::Int = 10`: The number of features. Should be at least 5.
- `noise::Union{Nothing, Float64} = nothing`: The standard deviation of the gaussian noise applied to the output.
- `random_state::Union{Int, Nothing} = nothing`: Determines random number generation for dataset noise. Pass an int for reproducible output across multiple function calls.
Reference: [link](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_friedman1.html)
"""
function make_friedman1(;   n_samples::Int = 100,
                            n_features::Int = 10,
                            noise::Float64 = 0.0,
                            random_state::Union{Int, Nothing} = nothing)::DataFrame

    (features, labels) = datasets.make_friedman1(   n_samples = n_samples,
                                                    n_features = n_features,
                                                    noise = noise,
                                                    random_state = random_state)

    return convert(features, labels)
end

"""
    function make_friedman2(;   n_samples::Int = 100,
                                noise::Float64 = 0.0,
                                random_state::Union{Int, Nothing} = nothing)::DataFrame
Generate the “Friedman #2” regression problem. Sklearn interface to make_friedman2.
#Arguments
- `n_samples::Int = 100`: The number of samples.
- `n_features::Int = 10`: The number of features. Should be at least 5.
- `noise::Union{Nothing, Float64} = nothing`: The standard deviation of the gaussian noise applied to the output.
- `random_state::Union{Int, Nothing} = nothing`: Determines random number generation for dataset noise. Pass an int for reproducible output across multiple function calls.
Reference: [link](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_friedman2.html)
"""
function make_friedman2(;   n_samples::Int = 100,
                            noise::Float64 = 0.0,
                            random_state::Union{Int, Nothing} = nothing)::DataFrame

    (features, labels) = datasets.make_friedman2(   n_samples = n_samples,
                                                    noise = noise,
                                                    random_state = random_state)

    return convert(features, labels)
end

"""
    function make_friedman3(;   n_samples::Int = 100,
                                noise::Float64 = 0.0,
                                random_state::Union{Int, Nothing} = nothing)::DataFrame
Generate the “Friedman #3” regression problem. Sklearn interface to make_friedman3.
#Arguments
- `n_samples::Int = 100`: The number of samples.
- `noise::Union{Nothing, Float64} = nothing`: The standard deviation of the gaussian noise applied to the output.
- `random_state::Union{Int, Nothing} = nothing`: Determines random number generation for dataset noise. Pass an int for reproducible output across multiple function calls.
Reference: [link](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_friedman3.html)
"""
function make_friedman3(;   n_samples::Int = 100,
                            noise::Float64 = 0.0,
                            random_state::Union{Int, Nothing} = nothing)::DataFrame

    (features, labels) = datasets.make_friedman3(   n_samples = n_samples,
                                                    noise = noise,
                                                    random_state = random_state)

    return convert(features, labels)
end

"""
    function make_low_rank_matrix(; n_samples::Int =100,
                                    n_features::Int =100,
                                    effective_rank::Int =10,
                                    tail_strength::Float64 =0.5,
                                    random_state::Union{Int, Nothing} = nothing)
Generate a mostly low rank matrix with bell-shaped singular values
#Arguments
- `n_samples::Int = 100`: The number of samples.
- `n_features::Int = 20`: The total number of features. These comprise `n_informative` informative features, `n_redundant` redundant features, `n_repeated` duplicated features and `n_features-n_informative-n_redundant-n_repeated` useless features drawn at random.
- `effective_rank::Int = 10`: The approximate number of singular vectors required to explain most of the data by linear combinations.
- `tail_strength::Float64 = 0.5`: The relative importance of the fat noisy tail of the singular values profile.
- `random_state::Union{Int, Nothing} = nothing`: Determines random number generation for dataset creation. Pass an int for reproducible output across multiple function calls. See Glossary.
Reference: [link](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_low_rank_matrix.html)
"""
function make_low_rank_matrix(; n_samples::Int = 100,
                                n_features::Int = 100,
                                effective_rank::Int = 10,
                                tail_strength::Float64 = 0.5,
                                random_state::Union{Int, Nothing} = nothing)

    features = datasets.make_low_rank_matrix(   n_samples = n_samples,
                                                n_features = n_features,
                                                effective_rank = effective_rank,
                                                tail_strength = tail_strength,
                                                random_state = random_state)
   return features
end

"""
function make_swiss_roll(;  n_samples::Int = 100,
                            noise::Float64 = 0.0,
                            random_state::Union{Int,Nothing} = nothing)
Generate a swiss roll dataset.
#Arguments
- `n_samples::Int = 100`: The number of samples.
- `noise::Float64 = 0.0 : Standard deviation of Gaussian noise added to the data.
- `random_state::Union{Int, Nothing} = nothing`: Determines random number generation for dataset creation. Pass an int for reproducible output across multiple function calls. See Glossary.
Reference: [link](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_swiss_roll.htmll)
"""
function make_swiss_roll(;  n_samples::Int = 100,
                            noise::Float64 = 0.0,
                            random_state::Union{Int,Nothing} = nothing)

   (features, labels) = datasets.make_swiss_roll(  n_samples = n_samples,
                                                   noise = noise,
                                                   random_state = random_state)

   return convert(features, labels)
end

"""
    function make_hastie_10_2(;  n_samples::Int = 12000,
                                    random_state::Union{Int,Nothing} = nothing)
Generates data for binary classification used in Hastie et al. 2009, Example 10.2.
#Arguments
- `n_samples::Int = 100`: The number of samples..
- `random_state::Union{Int, Nothing} = nothing`: Determines random number generation for dataset creation. Pass an int for reproducible output across multiple function calls. See Glossary.
Reference: [link](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_hastie_10_2.html)
"""
function make_hastie_10_2(; n_samples::Int = 12000,
                            random_state::Union{Int,Nothing} = nothing)

   (features, labels) = datasets.make_hastie_10_2(  n_samples = n_samples,
                                                    random_state = random_state)

   return convert(features, labels)
end

"""
    function make_gaussian_quantiles(;  mean::Array{<:Union{Number, Nothing}, 1} = [nothing],
                                        cov::Float64 = 1,
                                        n_samples::Int = 100,
                                        n_features::Int = 2,
                                        n_classes::Int = 3,
                                        shuffle::Bool = true,
                                        random_state::Union{Int, Nothing} = nothing)

Generate isotropic Gaussian and label samples by quantile.
#Arguments
- `mean::Union{Array{<:Number, 1}, Nothing} = nothing`: The mean of the multi-dimensional normal distribution. If None then use the origin (0, 0, …).
- `cov::Float64 = 1`: The covariance matrix will be this value times the unit matrix.
- `n_samples::Int = 100`: The total number of points equally divided among classes.
- `n_features::Int = 2`: The number of features for each sample.
- `n_classes::Int = 3`: The number of classes.
- `shuffle::Bool = true`: Shuffle the samples.
- `random_state::Union{Int, Nothing} = nothing`: Determines random number generation for dataset creation. Pass an int for reproducible output across multiple function calls. See Glossary.
Reference: [link](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_gaussian_quantiles.html)
"""
function make_gaussian_quantiles(;  mean::Union{Array{<:Number, 1}, Nothing} = nothing,
                                    cov::Float64 = 1.0,
                                    n_samples::Int = 100,
                                    n_features::Int = 2,
                                    n_classes::Int = 3,
                                    shuffle::Bool = true,
                                    random_state::Union{Int, Nothing} = nothing)

    typeof(mean) != Nothing && length(mean) != n_features && throw(DimensionMismatch("length of mean must be equal to n_features."))

    (features, labels) = datasets.make_gaussian_quantiles(  mean = mean,
                                                            cov = cov,
                                                            n_samples = n_samples,
                                                            n_features = n_features,
                                                            n_classes = n_classes,
                                                            shuffle = shuffle,
                                                            random_state = random_state)

    return convert(features, labels)
end


text:
using Documenter, SamplingHS3

makedocs(modules=[SamplingHS3],
        doctest=true)

deploydocs(deps   = Deps.pip("mkdocs", "python-markdown-math"),
    repo = "github.com/lithom/SamplingHS3.git",
    julia  = "0.6.2",
    osname = "Windows")

ext: jl


text:
@testset "Known gradients" begin
    x = ones(Float64, (2, 2, 1, 1))
    grid = Array{Float64}(undef, 2, 2, 2, 1)
    grid[:, 1, 1, 1] .= (-1, -1)
    grid[:, 2, 1, 1] .= (1, -1)
    grid[:, 1, 2, 1] .= (-1, 1)
    grid[:, 2, 2, 1] .= (1, 1)

    ∇grid_true = Array{Float64}(undef, size(grid))
    ∇grid_true[:, :, 1, 1] = [[0.0, 0.0] [-0.5, 0.0]]
    ∇grid_true[:, :, 2, 1] = [[0.0, -0.5] [-0.5, -0.5]]

    padding_mode = :zeros
    sampled = grid_sample(x, grid; padding_mode=padding_mode)
    @test x == sampled
    @test eltype(sampled) == Float64
    external_grad = ones(size(sampled))
    ∇input, ∇grid = ∇grid_sample(external_grad, x, grid; padding_mode=padding_mode)
    @test ∇input == x
    @test ∇grid == ∇grid_true
    @test eltype(∇input) == Float64
    @test eltype(∇grid) == Float64

    # ∇grid from FiniteDifferences is incorrent in case when 0-padding.
    # gradtest(grid_sample, x, grid; fkwargs=(padding_mode=padding_mode,))

    padding_mode = :border
    fill!(∇grid_true, 0.0)
    sampled = grid_sample(x, grid; padding_mode=padding_mode)
    @test x == sampled
    @test eltype(sampled) == Float64
    external_grad = ones(size(sampled))
    ∇input, ∇grid = ∇grid_sample(external_grad, x, grid; padding_mode=padding_mode)
    @test ∇input == x
    @test ∇grid == ∇grid_true
    @test eltype(∇input) == Float64
    @test eltype(∇grid) == Float64

    gradtest(grid_sample, x, grid; fkwargs=(padding_mode=padding_mode,))
end

@testset "Test out-of-bounds for different paddings" begin
    x = ones(Float64, (2, 2, 1, 1))
    grid = Array{Float64}(undef, 2, 3, 2, 1)
    grid[:, 1, 1, 1] .= (-3, -1)
    grid[:, 2, 1, 1] .= (0, -1)
    grid[:, 3, 1, 1] .= (3, -1)
    grid[:, 1, 2, 1] .= (-1, 3)
    grid[:, 2, 2, 1] .= (0, 1)
    grid[:, 3, 2, 1] .= (1, 3)

    # With 0-padding, out-of-bound values are will contribute nothing to
    # the output values, because they are too far from any bound.
    y = grid_sample(x, grid; padding_mode=:zeros)
    y_true = reshape(Float64[[0, 1, 0] [0, 1, 0]], size(y))
    @test y_true == y

    # With border-padding, out-of-bound values simly become border values
    # and the result should be all ones.
    y = grid_sample(x, grid; padding_mode=:border)
    y_true = ones(Float64, size(y))
    @test y_true == y
end


text:
#=
0. gpioWaveClear
1. gpioWaveAddGeneric ( gpioWaveAddNew, gpioWaveAddSerial )
2. gpioWaveCreate (gpioWaveCreatePad, )
3. gpioWaveTxSend (gpioWaveChain, )

gpioWaveTxAt
gpioWaveDelete
gpioWaveTxStop
=#

"""
    mutable struct DisplayBCD <: DisplayNoBuffer
        digits_pins::AbstractVector{Int}      # pin numbers to control digits: [LSD ... MSD]
        sectors_pins::Tuple{Int,Int,Int,Int,Int,Int,Int,Int}  # pin numbers for binary code: (d0,d1,d2,d3,-,-,-,DP)
        buffer::AbstractVector{UInt8}         # internal storage of digits values
        usDelay::Int                          # duration of period when digit is on
        inverted_digits::Bool                 # if digit pins states must be inverted, i.e. 1 means LOW pin state
        inverted_sectors::Bool                # if input pins and dp pin states must be inverted, i.e. 1 means LOW pin state
        intensity::Int                        # 16 max
        limit::Int
        test_mode::Bool
    end
"""
mutable struct DisplayBCD <: DisplayNoBuffer
    digits_pins::AbstractVector{Int}
    sectors_pins::Tuple{Int,Int,Int,Int,Int,Int,Int,Int}
    buffer::AbstractVector{UInt8}
    usDelay::Real
    inverted_digits::Bool
    inverted_sectors::Bool
    intensity::Int
    limit::Int
    test_mode::Bool
end

"""
    function DisplayBCD(
        digits_pins::AbstractVector{Int},
        sectors_pins::Tuple{Int,Int,Int,Int,Int}, # A, B, C, D, DP
        scan_rate::Real = 800,                    # Hz
        inverted_digits::Bool = false,
        inverted_sectors::Bool = false
    )

Creates device representing numerical display with several digits under control of the BCD chip.
The number of display digits equal to `digits_pins` count.

## Arguments

- `digits_pins` : Vector of GPIO pin numbers connected to anode. The HIGH state means the digit is on. LOW means off.
    The first pin in array should manage the less significant digit.

- `sectors_pins` : Tuple of length 5 consisting of GPIO numbers representing the 4-bit code of a digit.
    The sequence of pins is the following: (d0, d1, d2, d3, DP).

- `scan_rate` : refresh rate of digits in Hz.
    The digits in display are controlled by impulses of `digits_pins`.
    This argument sets the time period for display of one digit.
    If `scan_rate=1000` the width will be recalculated as `1/1000 = 1e-3` second or `1e3` microsecond.
    The default value is 800 Hz.

- `inverted_digits` : This option inverts `digit_pins` active states.

- `inverted_sectors` : This option inverts `sectors_pins` active states.
"""
function DisplayBCD(
    digits_pins::AbstractVector{Int},
    sectors_pins::Tuple{Int,Int,Int,Int,Int},
    scan_rate::Real = 800, # Hz
    inverted_digits::Bool = false,
    inverted_sectors::Bool = false
)
    @assert 1 <= length(digits_pins) <= 8 "The package supports up to 8 digits, got $(length(digits_pins))"

    if PiGPIOC.gpioInitialise() < 0
        throw("pigpio initialisation failed.")
    else
        @info "pigpio initialised okay."
    end

    # set frequency
    # TODO: display real frequency
    min_rate = ceil(Int, 10^6 / PiGPIOC.gpioWaveGetMaxMicros())
    max_rate = 10^6 # based on minimal usDelay = 1
    @assert min_rate <= scan_rate <= max_rate "Scan rate must be between $min_rate and $max_rate, got $scan_rate"
    usDelay = 10^6 / scan_rate # conversion to us

    # init pins
    PiGPIOC.gpioSetMode(digits_pins, PiGPIOC.PI_OUTPUT)
    PiGPIOC.gpioWrite(digits_pins, 0)
    PiGPIOC.gpioSetMode(sectors_pins, PiGPIOC.PI_OUTPUT)
    PiGPIOC.gpioWrite(sectors_pins, 0)

    buffer = fill(0b10001111, length(digits_pins)) # set empty display

    DisplayBCD(
        digits_pins,
        (sectors_pins[1], sectors_pins[2], sectors_pins[3], sectors_pins[4], -1, -1, -1, sectors_pins[5]),
        buffer,
        usDelay,
        inverted_digits,
        inverted_sectors,
        5,                   # default intensity
        length(digits_pins), # default limit,
        false
    )
end

decode(::DisplayBCD) = 0b00000000

# this function is not avalable for BCD displays
function decode_mode(d::DisplayBCD, decode::UInt8 = 0b1111_1111)
    @warn "`decode_mode` method doesn't work for DisplayBCD"
end

# generate test wave with all sectors are active and maximal intensity
function generate_test_wave(d::DisplayBCD)
    pulse = PiGPIOC.gpioPulse_t[]

    # set active state
    activePeriod = ceil(Int, d.usDelay)
    for i in 1:size(d)
        gpioOn = 0x0
        gpioOff = 0x0

        # digits
        for j in 1:size(d)
            if d.digits_pins[j] >= 0
                if xor(i == j, d.inverted_digits)
                    gpioOn |= 1 << d.digits_pins[j]
                else
                    gpioOff |= 1 << d.digits_pins[j]
                end
            end
        end

        # sectors
        value = 0b1000_1000 # representation of 8. = 0x1000
        for j in 1:8
            if d.sectors_pins[j] >= 0
                if xor(value % 2 == 1, d.inverted_sectors)
                    gpioOn |= 1 << d.sectors_pins[j]
                else
                    gpioOff |= 1 << d.sectors_pins[j]
                end
            end
            value >>= 1
        end

        push!(pulse, PiGPIOC.gpioPulse_t(gpioOn, gpioOff, activePeriod)) # on, off, usDelay
    end

    PiGPIOC.gpioWaveAddGeneric(size(d), pulse)
    wave_id = PiGPIOC.gpioWaveCreate()
    if wave_id < 0
        # return Upon success a wave id greater than or equal to 0 is returned,
        # otherwise PI_EMPTY_WAVEFORM, PI_TOO_MANY_CBS, PI_TOO_MANY_OOL,
        # or PI_NO_WAVEFORM_ID
        throw("Error in 'PiGPIOC.gpioWaveCreate()' with code: $(wave_id)")
    end

    return wave_id
end


text:
__precompile__()
module Librsync

const depsfile = joinpath(dirname(dirname(@__FILE__)), "deps", "deps.jl")
if isfile(depsfile)
  include(depsfile)
else
  error("Librsync not properly installed. Please run Pkg.build(\"Librsync\")")
end

include("librsync_types.jl")
include("librsync_core.jl")

# const librsync_version = unsafe_string(unsafe_load(cglobal((:rs_librsync_version, librsync), Ptr{Cchar})))
const MD4_SUM_LENGTH = unsafe_load(cglobal((:RS_MD4_SUM_LENGTH, Librsync.librsync), Cint))
const BLAKE2_SUM_LENGTH = unsafe_load(cglobal((:RS_BLAKE2_SUM_LENGTH, Librsync.librsync), Cint))
const inbuflen = unsafe_load(cglobal((:rs_inbuflen, Librsync.librsync), Cint))
const outbuflen = unsafe_load(cglobal((:rs_outbuflen, Librsync.librsync), Cint))

end # module


text:
import Base: copy, copyto!, ==
import RandomNumbers: AbstractRNG, seed_type

"""
```julia
WrappedRNG{R, T1, T2} <: AbstractRNG{T2}
WrappedRNG(base_rng, T2)
WrappedRNG(R, T2, args...)
```

Wrap a RNG which originally provides output in T1 into a RNG that provides output in T2.

# Examples
```jldoctest
julia> r = Xorshifts.Xorshift128Star(123);

julia> RandomNumbers.output_type(r)
UInt64

julia> r1 = WrappedRNG(r, UInt32);

julia> RandomNumbers.output_type(r1)
UInt32

julia> r2 = WrappedRNG(Xorshifts.Xorshift128Star, UInt32, 123);

julia> RandomNumbers.output_type(r2)
UInt32

julia> @Test.test rand(r1, UInt32, 3) == rand(r2, UInt32, 3)
Test Passed
```
"""
mutable struct WrappedRNG{R<:AbstractRNG, T1<:BitTypes, T2<:BitTypes} <: AbstractRNG{T2}
    base_rng::R
    x::T1
    p::Int
    WrappedRNG{R, T1, T2}() where {R <: AbstractRNG, T1 <: BitTypes, T2 <: BitTypes} =
        (@assert T1 ≠ T2; new())
end

function WrappedRNG(base_rng::AbstractRNG{T1}, ::Type{T2}) where {T1 <: BitTypes, T2 <: BitTypes}
    wr = WrappedRNG{typeof(base_rng), T1, T2}()
    wr.base_rng = copy(base_rng)
    if sizeof(T1) > sizeof(T2)
        wr.x = rand(wr.base_rng, T1)
    end
    wr.p = 0
    wr
end

function WrappedRNG(::Type{R}, ::Type{T2}, args...) where {R <: AbstractRNG, T2 <: BitTypes}
    base_rng = R(args...)
    WrappedRNG(base_rng, T2)
end

WrappedRNG(base_rng::WrappedRNG{R, T1, T2}, ::Type{T3}) where
    {R <: AbstractRNG, T1 <: BitTypes, T2 <: BitTypes, T3 <: BitTypes} = WrappedRNG(base_rng.base_rng, T3)

seed_type(::Type{WrappedRNG{R, T1, T2}}) where {R, T1, T2} = seed_type(R)

function copyto!(dest::R, src::R) where R <: WrappedRNG
    copyto!(dest.base_rng, src.base_rng)
    dest.x = src.x
    dest.p = src.p
    dest
end

function copy(src::R) where R <: WrappedRNG
    wr = R()
    wr.base_rng = copy(src.base_rng)
    wr.x = src.x
    wr.p = src.p
    wr
end

==(r1::R, r2::R) where R <: WrappedRNG = r1.base_rng == r2.base_rng && r1.x == r2.x && r1.p == r2.p

function seed!(wr::WrappedRNG, seed...)
    seed!(wr.base_rng, seed...)
        if sizeof(T1) > sizeof(T2)
            wr.x = rand(wr.base_rng, T1)
        end
    wr.p = 0
    wr
end

@inline function rand(rng::WrappedRNG{R, T1, T2}, ::Type{T2}) where
        {R <: AbstractRNG, T1 <: BitTypes, T2 <: BitTypes}
    s1 = sizeof(T1)
    s2 = sizeof(T2)
    if s2 >= s1
        t = rand(rng.base_rng, T1) % T2
        for i in 2:(s2 ÷ s1)
            t |= (rand(rng.base_rng, T1) % T2) << ((s1 << 3) * (i - 1))
        end
    else
        t = rng.x % T2
        rng.p += 1
        if rng.p == s1 ÷ s2
            rng.p = 0
            rng.x = rand(rng.base_rng, T1)
        else
            rng.x >>= s2 << 3
        end
    end
    return t
end


text:
using Test

include("./test.jl")


text:
module OperationsTest

import Random: randstring
import LibGit2
using Test
using UUIDs

using Pkg3
using Pkg3.Types

import Random: randstring
import LibGit2

function temp_pkg_dir(fn::Function)
    local env_dir
    local old_load_path
    local old_depot_path
    try
        old_load_path = copy(LOAD_PATH)
        old_depot_path = copy(DEPOT_PATH)
        empty!(LOAD_PATH)
        empty!(DEPOT_PATH)
        mktempdir() do env_dir
            mktempdir() do depot_dir
                pushfirst!(LOAD_PATH, env_dir)
                pushfirst!(DEPOT_PATH, depot_dir)
                # Add the standard library paths back
                vers = "v$(VERSION.major).$(VERSION.minor)"
                push!(LOAD_PATH, abspath(Sys.BINDIR, "..", "local", "share", "julia", "site", vers))
                push!(LOAD_PATH, abspath(Sys.BINDIR, "..", "share", "julia", "site", vers))
                fn(env_dir)
            end
        end
    finally
        append!(LOAD_PATH, old_load_path)
        append!(DEPOT_PATH, old_depot_path)
    end
end

const TEST_PKG = (name = "Example", uuid = UUID("7876af07-990d-54b4-ab0e-23690620f79a"))
isinstalled(pkg) = Base.locate_package(Base.PkgId(pkg.uuid, pkg.name)) !== nothing

temp_pkg_dir() do project_path

    @testset "simple add and remove with preview" begin
        Pkg3.init(project_path)
        Pkg3.add(TEST_PKG.name; preview = true)
        @test !isinstalled(TEST_PKG)
        Pkg3.add(TEST_PKG.name)
        @test isinstalled(TEST_PKG)
        @eval import $(Symbol(TEST_PKG.name))
        Pkg3.rm(TEST_PKG.name; preview = true)
        @test isinstalled(TEST_PKG)
        Pkg3.rm(TEST_PKG.name)
        @test !isinstalled(TEST_PKG)
    end

    @testset "package with wrong UUID" begin
        @test_throws ResolverError Pkg3.add(PackageSpec(TEST_PKG.name, UUID(UInt128(1))))
    end

    @testset "adding and upgrading different versions" begin
        # VersionNumber
        Pkg3.add(PackageSpec(TEST_PKG.name, v"0.3"))
        @test Pkg3.installed()[TEST_PKG.name] == v"0.3"
        Pkg3.add(PackageSpec(TEST_PKG.name, v"0.3.1"))
        @test Pkg3.installed()[TEST_PKG.name] == v"0.3.1"
        Pkg3.rm(TEST_PKG.name)

        # VersionRange
        Pkg3.add(PackageSpec(TEST_PKG.name, VersionSpec(VersionRange("0.3.0-0.3.2"))))
        @test Pkg3.installed()[TEST_PKG.name] == v"0.3.2"
        Pkg3.up(; level = UpgradeLevel(:patch))
        @test Pkg3.installed()[TEST_PKG.name] == v"0.3.3"
        Pkg3.up(; level = UpgradeLevel(:minor))
        @test Pkg3.installed()[TEST_PKG.name].minor != 3
        Pkg3.rm(TEST_PKG.name)
    end

    @testset "testing" begin
        # TODO: Check that preview = true doesn't actually execute the test
        # TODO: Test-only dependencies
        Pkg3.add(TEST_PKG.name)
        Pkg3.test(TEST_PKG.name; coverage=true)
        pkgdir = Base.locate_package(Base.PkgId(TEST_PKG.uuid, TEST_PKG.name))
        # No coverage files being generated?
        @test_broken TEST_PKG.name * ".cov" in readdir(pkgdir)
        Pkg3.rm(TEST_PKG.name)
    end

    @testset "pinning / freeing" begin
        Pkg3.add(TEST_PKG.name)
        old_v = Pkg3.installed()[TEST_PKG.name]
        Pkg3.pin(PackageSpec(TEST_PKG.name, v"0.2"))
        @test Pkg3.installed()[TEST_PKG.name].minor == 2
        Pkg3.up(TEST_PKG.name)
        @test Pkg3.installed()[TEST_PKG.name].minor == 2
        Pkg3.free(TEST_PKG.name)
        Pkg3.up()
        @test Pkg3.installed()[TEST_PKG.name] == old_v
        Pkg3.rm(TEST_PKG.name)
    end

    @testset "checking out / freeing" begin
        Pkg3.add(TEST_PKG.name)
        old_v = Pkg3.installed()[TEST_PKG.name]
        Pkg3.rm(TEST_PKG.name)
        mktempdir() do devdir
            Pkg3.checkout(TEST_PKG.name; path = devdir)
            @test isinstalled(TEST_PKG)
            @test Pkg3.installed()[TEST_PKG.name] > old_v
            test_pkg_main_file = joinpath(devdir, TEST_PKG.name, "src", TEST_PKG.name * ".jl")
            @test isfile(test_pkg_main_file)
            # Pkg3 #152
            write(test_pkg_main_file,
                """
                module Example
                    export hello, domath
                    const example2path = joinpath(@__DIR__, "..", "deps", "deps.jl")
                    if !isfile(example2path)
                        error("Example is not installed correctly")
                    end
                    hello(who::String) = "Hello, \$who"
                    domath(x::Number) = x + 5
                end
                """)
            mkpath(joinpath(devdir, TEST_PKG.name, "deps"))
            write(joinpath(devdir, TEST_PKG.name, "deps", "build.jl"),
                """
                touch("deps.jl")
                """
            )
            Pkg3.build(TEST_PKG.name)
            @test isfile(joinpath(devdir, TEST_PKG.name, "deps", "deps.jl"))
            Pkg3.test(TEST_PKG.name)
            Pkg3.free(TEST_PKG.name)
            @test Pkg3.installed()[TEST_PKG.name] == old_v
        end
        mktempdir() do tmp
            withenv("JULIA_PKG_DEVDIR" => tmp) do
                @test_throws CommandError Pkg3.checkout(TEST_PKG.name, "nonexisting_branch",)
            end
        end
    end

    @testset "package name in resolver errors" begin
        try
            Pkg3.add([PackageSpec(TEST_PKG.name, VersionSpec(v"55"))])
        catch e
            @test contains(sprint(showerror, e), TEST_PKG.name)
        end
    end

    @testset "check logging" begin
        usage = Pkg3.TOML.parse(String(read(joinpath(Pkg3.logdir(), "manifest_usage.toml"))))
        @test any(x -> startswith(x, joinpath(project_path, "Manifest.toml")), keys(usage))
    end

    @testset "adding nonexisting packages" begin
        nonexisting_pkg = randstring(14)
        @test_throws CommandError Pkg3.add(nonexisting_pkg)
        @test_throws CommandError Pkg3.up(nonexisting_pkg)
    end

    Pkg3.rm(TEST_PKG.name)
end

end # module

max_stars_repo_name: blegat/julia

max_issues_repo_name: blegat/julia

max_forks_repo_name: blegat/julia


text:
# Benchmark the counting of alanine residues in a PDB file

using MIToS.PDB

pdb_filepath = "data/1AKE.pdb"
struc = read(pdb_filepath, PDBFile)

counter() = count(res -> res.id.name == "ALA", struc)

# Run to JIT compile
counter()

elapsed = @elapsed counter()

println(elapsed)

max_stars_repo_name: wojdyr/pdb-benchmarks

max_issues_repo_name: wojdyr/pdb-benchmarks

max_forks_repo_name: wojdyr/pdb-benchmarks


text:
# grep "#define" cpxconst.h | grep "CPX_PARAM_|CPXPARAM" | awk '{ printf("\"%s\" => %s,\n",$2,$3) }'
const paramName2Indx = Dict(
"CPX_PARAM_ADVIND" => 1001,
"CPX_PARAM_AGGFILL" => 1002,
"CPX_PARAM_AGGIND" => 1003,
"CPX_PARAM_CLOCKTYPE" => 1006,
"CPX_PARAM_CRAIND" => 1007,
"CPX_PARAM_DEPIND" => 1008,
"CPX_PARAM_DPRIIND" => 1009,
"CPX_PARAM_PRICELIM" => 1010,
"CPX_PARAM_EPMRK" => 1013,
"CPX_PARAM_EPOPT" => 1014,
"CPX_PARAM_EPPER" => 1015,
"CPX_PARAM_EPRHS" => 1016,
"CPX_PARAM_SIMDISPLAY" => 1019,
"CPX_PARAM_ITLIM" => 1020,
"CPX_PARAM_ROWREADLIM" => 1021,
"CPX_PARAM_NETFIND" => 1022,
"CPX_PARAM_COLREADLIM" => 1023,
"CPX_PARAM_NZREADLIM" => 1024,
"CPX_PARAM_OBJLLIM" => 1025,
"CPX_PARAM_OBJULIM" => 1026,
"CPX_PARAM_PERIND" => 1027,
"CPX_PARAM_PERLIM" => 1028,
"CPX_PARAM_PPRIIND" => 1029,
"CPX_PARAM_PREIND" => 1030,
"CPX_PARAM_REINV" => 1031,
"CPX_PARAM_SCAIND" => 1034,
"CPX_PARAM_SCRIND" => 1035,
"CPX_PARAM_SINGLIM" => 1037,
"CPX_PARAM_TILIM" => 1039,
"CPX_PARAM_PREDUAL" => 1044,
"CPX_PARAM_PREPASS" => 1052,
"CPX_PARAM_DATACHECK" => 1056,
"CPX_PARAM_REDUCE" => 1057,
"CPX_PARAM_PRELINEAR" => 1058,
"CPX_PARAM_LPMETHOD" => 1062,
"CPX_PARAM_QPMETHOD" => 1063,
"CPX_PARAM_WORKDIR" => 1064,
"CPX_PARAM_WORKMEM" => 1065,
"CPX_PARAM_THREADS" => 1067,
"CPX_PARAM_CONFLICTALG" => 1073,
"CPX_PARAM_CONFLICTDISPLAY" => 1074,
"CPX_PARAM_SIFTDISPLAY" => 1076,
"CPX_PARAM_SIFTALG" => 1077,
"CPX_PARAM_SIFTITLIM" => 1078,
"CPX_PARAM_MPSLONGNUM" => 1081,
"CPX_PARAM_MEMORYEMPHASIS" => 1082,
"CPX_PARAM_NUMERICALEMPHASIS" => 1083,
"CPX_PARAM_FEASOPTMODE" => 1084,
"CPX_PARAM_PARALLELMODE" => 1109,
"CPX_PARAM_TUNINGMEASURE" => 1110,
"CPX_PARAM_TUNINGREPEAT" => 1111,
"CPX_PARAM_TUNINGTILIM" => 1112,
"CPX_PARAM_TUNINGDISPLAY" => 1113,
"CPX_PARAM_WRITELEVEL" => 1114,
"CPX_PARAM_RANDOMSEED" => 1124,
"CPX_PARAM_DETTILIM" => 1127,
"CPX_PARAM_FILEENCODING" => 1129,
"CPX_PARAM_APIENCODING" => 1130,
"CPX_PARAM_OPTIMALITYTARGET" => 1131,
"CPX_PARAM_CLONELOG" => 1132,
"CPX_PARAM_TUNINGDETTILIM" => 1139,
"CPX_PARAM_CPUMASK" => 1144,
"CPX_PARAM_SOLUTIONTYPE" => 1147,
"CPX_PARAM_SIFTSIM" => 1158,
"CPX_PARAM_WORKERALG" => 1500,
"CPX_PARAM_BENDERSSTRATEGY" => 1501,
"CPX_PARAM_BENDERSFEASCUTTOL" => 1509,
"CPX_PARAM_BENDERSOPTCUTTOL" => 1510,
"CPX_PARAM_BRDIR" => 2001,
"CPX_PARAM_BTTOL" => 2002,
"CPX_PARAM_CLIQUES" => 2003,
"CPX_PARAM_COEREDIND" => 2004,
"CPX_PARAM_COVERS" => 2005,
"CPX_PARAM_CUTLO" => 2006,
"CPX_PARAM_CUTUP" => 2007,
"CPX_PARAM_EPAGAP" => 2008,
"CPX_PARAM_EPGAP" => 2009,
"CPX_PARAM_EPINT" => 2010,
"CPX_PARAM_MIPDISPLAY" => 2012,
"CPX_PARAM_MIPINTERVAL" => 2013,
"CPX_PARAM_INTSOLLIM" => 2015,
"CPX_PARAM_NODEFILEIND" => 2016,
"CPX_PARAM_NODELIM" => 2017,
"CPX_PARAM_NODESEL" => 2018,
"CPX_PARAM_OBJDIF" => 2019,
"CPX_PARAM_MIPORDIND" => 2020,
"CPX_PARAM_RELOBJDIF" => 2022,
"CPX_PARAM_STARTALG" => 2025,
"CPX_PARAM_SUBALG" => 2026,
"CPX_PARAM_TRELIM" => 2027,
"CPX_PARAM_VARSEL" => 2028,
"CPX_PARAM_BNDSTRENIND" => 2029,
"CPX_PARAM_HEURFREQ" => 2031,
"CPX_PARAM_MIPORDTYPE" => 2032,
"CPX_PARAM_CUTSFACTOR" => 2033,
"CPX_PARAM_RELAXPREIND" => 2034,
"CPX_PARAM_PRESLVND" => 2037,
"CPX_PARAM_BBINTERVAL" => 2039,
"CPX_PARAM_FLOWCOVERS" => 2040,
"CPX_PARAM_IMPLBD" => 2041,
"CPX_PARAM_PROBE" => 2042,
"CPX_PARAM_GUBCOVERS" => 2044,
"CPX_PARAM_STRONGCANDLIM" => 2045,
"CPX_PARAM_STRONGITLIM" => 2046,
"CPX_PARAM_FRACCAND" => 2048,
"CPX_PARAM_FRACCUTS" => 2049,
"CPX_PARAM_FRACPASS" => 2050,
"CPX_PARAM_FLOWPATHS" => 2051,
"CPX_PARAM_MIRCUTS" => 2052,
"CPX_PARAM_DISJCUTS" => 2053,
"CPX_PARAM_AGGCUTLIM" => 2054,
"CPX_PARAM_MIPCBREDLP" => 2055,
"CPX_PARAM_CUTPASS" => 2056,
"CPX_PARAM_MIPEMPHASIS" => 2058,
"CPX_PARAM_SYMMETRY" => 2059,
"CPX_PARAM_DIVETYPE" => 2060,
"CPX_PARAM_RINSHEUR" => 2061,
"CPX_PARAM_SUBMIPNODELIM" => 2062,
"CPX_PARAM_LBHEUR" => 2063,
"CPX_PARAM_REPEATPRESOLVE" => 2064,
"CPX_PARAM_PROBETIME" => 2065,
"CPX_PARAM_POLISHTIME" => 2066,
"CPX_PARAM_REPAIRTRIES" => 2067,
"CPX_PARAM_EPLIN" => 2068,
"CPX_PARAM_EPRELAX" => 2073,
"CPX_PARAM_FPHEUR" => 2098,
"CPX_PARAM_EACHCUTLIM" => 2102,
"CPX_PARAM_SOLNPOOLCAPACITY" => 2103,
"CPX_PARAM_SOLNPOOLREPLACE" => 2104,
"CPX_PARAM_SOLNPOOLGAP" => 2105,
"CPX_PARAM_SOLNPOOLAGAP" => 2106,
"CPX_PARAM_SOLNPOOLINTENSITY" => 2107,
"CPX_PARAM_POPULATELIM" => 2108,
"CPX_PARAM_MIPSEARCH" => 2109,
"CPX_PARAM_MIQCPSTRAT" => 2110,
"CPX_PARAM_ZEROHALFCUTS" => 2111,
"CPX_PARAM_POLISHAFTEREPAGAP" => 2126,
"CPX_PARAM_POLISHAFTEREPGAP" => 2127,
"CPX_PARAM_POLISHAFTERNODE" => 2128,
"CPX_PARAM_POLISHAFTERINTSOL" => 2129,
"CPX_PARAM_POLISHAFTERTIME" => 2130,
"CPX_PARAM_MCFCUTS" => 2134,
"CPX_PARAM_MIPKAPPASTATS" => 2137,
"CPX_PARAM_AUXROOTTHREADS" => 2139,
"CPX_PARAM_INTSOLFILEPREFIX" => 2143,
"CPX_PARAM_PROBEDETTIME" => 2150,
"CPX_PARAM_POLISHAFTERDETTIME" => 2151,
"CPX_PARAM_LANDPCUTS" => 2152,
"CPX_PARAM_RAMPUPDURATION" => 2163,
"CPX_PARAM_RAMPUPDETTILIM" => 2164,
"CPX_PARAM_RAMPUPTILIM" => 2165,
"CPX_PARAM_LOCALIMPLBD" => 2181,
"CPX_PARAM_BQPCUTS" => 2195,
"CPX_PARAM_RLTCUTS" => 2196,
"CPX_PARAM_BAREPCOMP" => 3002,
"CPX_PARAM_BARGROWTH" => 3003,
"CPX_PARAM_BAROBJRNG" => 3004,
"CPX_PARAM_BARALG" => 3007,
"CPX_PARAM_BARCOLNZ" => 3009,
"CPX_PARAM_BARDISPLAY" => 3010,
"CPX_PARAM_BARITLIM" => 3012,
"CPX_PARAM_BARMAXCOR" => 3013,
"CPX_PARAM_BARORDER" => 3014,
"CPX_PARAM_BARSTARTALG" => 3017,
"CPX_PARAM_BARCROSSALG" => 3018,
"CPX_PARAM_BARQCPEPCOMP" => 3020,
"CPX_PARAM_QPNZREADLIM" => 4001,
"CPX_PARAM_CALCQCPDUALS" => 4003,
"CPX_PARAM_QPMAKEPSDIND" => 4010,
"CPX_PARAM_QTOLININD" => 4012,
"CPX_PARAM_NETITLIM" => 5001,
"CPX_PARAM_NETEPOPT" => 5002,
"CPX_PARAM_NETEPRHS" => 5003,
"CPX_PARAM_NETPPRIIND" => 5004,
"CPX_PARAM_NETDISPLAY" => 5005
)


text:
function TwoBodiesBaryFrame(arg0::CelestialBody, arg1::CelestialBody)
    return TwoBodiesBaryFrame((CelestialBody, CelestialBody), arg0, arg1)
end


text:
# roseplot

#JL Plots
using Random, Plots
Random.seed!(2018)

n = 24
R = rand(n+1)
θ = 0:2pi/n:2pi
Plots.plot(θ, R, proj=:polar, line=:steppre, lims=(0,1), c=:steelblue, linewidth=2)

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
import Base.reset
export mock, calls, reset, call

global _calls = Dict{Function, Vector}()

function mock(;return_value::Any=nothing, side_effect::Union(Function, Nothing)=nothing)
  function impl(args::Any...; kwargs...)
    _calls[impl] = vcat(calls(impl), [(args, kwargs)])
    if side_effect==nothing
      return return_value
    else
      return side_effect(args..., kwargs...)
    end
  end
  return impl
end

function calls(fn::Function)
  return get(_calls, fn, [])
end

function reset(fn::Function)
  _calls[fn] = []
  return
end

function call(args::Any...; kwargs...)
  return args, kwargs
end


text:
#===============================================================================
    Definitions of abstract types reflecting conceptual inheritence structure
===============================================================================#

"""
Supertype of all MCMC updates
"""
abstract type MCMCUpdate end

"""
    MCMCParamUpdate <: MCMCUpdate

Supertype of all updates that make changes to the main MCMC parameter called
`state` in the `GlobalWorkspace`.
"""
abstract type MCMCParamUpdate <: MCMCUpdate end

"""
    MCMCGradientBasedUpdate <: MCMCParamUpdate

Supertype of all updates that require local gradient information to perform
update of `state`.
"""
abstract type MCMCGradientBasedUpdate <: MCMCParamUpdate end

"""
    MCMCConjugateUpdate <: MCMCUpdate

Supertype of all conjugate updates, for which sampling can be done directly,
without resorting to a Metropolis-Hastings algorithm
"""
abstract type MCMCConjugateUpdate <: MCMCUpdate end

"""
    MCMCConjugateParamUpdate <: MCMCParamUpdate

Supertype of all conjugate updates that update parameters and for which sampling
can be done directly, without resorting to a Metropolis-Hastings algorithm
"""
abstract type MCMCConjugateParamUpdate <: MCMCParamUpdate end

"""
    MCMCImputation <: MCMCUpdate

Supertype of all updates that do not make any changes to the main parameter
`state` in `GlobalWorkspace`, but instead, perform sampling on any auxiliary
variables that are not of direct interest to the MCMC chain.
"""
abstract type MCMCImputation <: MCMCUpdate end

"""
Supertype of all decorators to update schemes. In this context, we use the word
`decorator` to refer to any additional information that needs to be conveyed to
the MCMC sampler that is not directly extractable from the update-objects
themselves (for instance, a change in the delimitation of blocks that is made
in-between updates).
"""
abstract type MCMCUpdateDecorator end

"""
    isdecorator(u)

Returns true if `u` is a subtype of decorators.
"""
isdecorator(u) = (typeof(u) <: MCMCUpdateDecorator)

"""
Supertype of all workspaces—i.e. of structs that gather in one place various
objects that the MCMC sampler operates on.
"""
abstract type Workspace end

"""
    GlobalWorkspace{T} <: Workspace

Supertype of all global workspaces. Each MCMC sampler must have a unique global
workspace, which contains `state`, `state_history`, `state_proposal_history`,
`acceptance_history` and `data`. `state` is the paramater vector that the MCMC
sampling is done for (other names being self-explanatory).
"""
abstract type GlobalWorkspace{T} <: Workspace end

"""
    LocalWorkspace{T} <: Workspace

Supertype of all local workspaces. Local workspace should contain any additional
gathering of objects that are needed by specific updates, but are not are not
already in a global workspace. Each MCMC update has its own `LocalWorkspace`.
"""
abstract type LocalWorkspace{T} <: Workspace end

"""
Supertype of all transition kernels that perform updates on the main `state`
of the global workspace. Conceptually, these are used by subtypes of
`MCMCParamUpdate` to perform actual sampling and instances of types inheriting
from `TransitionKernel` are usually member objects of the instances inheriting
from `MCMCParamUpdate`.
"""
abstract type TransitionKernel end

"""
Supertype for all adaptation schemes.
"""
abstract type Adaptation end

"""
Supertype of all backends for the MCMC sampler.
"""
abstract type MCMCBackend end

"""
    GenericMCMCBackend <: MCMCBackend

A flag that no specific backend is passed.
"""
struct GenericMCMCBackend <: MCMCBackend end

abstract type ChainStats end

struct Previous end
struct Proposal end
struct PreMCMCStep end
struct PostMCMCStep end

#NOTE also appears in DiffusionDefinition.jl, but not used there.
"""
    remove_curly(::Type{K}) where K

Utility function that removes all type-specifiers listed in the curly brackets.

# Examples
```julia-repl
julia> remove_curly(Array{Float64,1})
Array
```
"""
@generated function remove_curly(::Type{K}) where K
    name_without_curly = Meta.parse(string(K)).args[1]
    :( $name_without_curly )
end


text:
    #
    # Two-photon absorption by monochromatic and equally-polarized photons, usually from the same beam.
    #
    """
    `struct  MultiPhotonDeExcitation.Channel_2pAbsorptionMonochromatic`
        ... defines a type for a two-photon absorption channel for the absorption of monochromatic light with well-defined
            multipolarities.

        + K              ::AngularJ64             ... Rank K of the channel.
        + omega          ::Float64                ... omega.
        + multipole1     ::EmMultipole            ... Multipole M1.
        + multipole2     ::EmMultipole            ... Multipole M2.
        + gauge          ::EmGauge                ... Gauge for dealing with the (coupled) radiation field.
        + Jsym           ::LevelSymmetry          ... Symmetry of the Green function channel used in the summation.
        + amplitude      ::Complex{Float64}       ... reduced two-photon absorption amplitude U^(K, 2gamma, absorption) (..)
                                                      associated with the given channel.
    """
    struct  Channel_2pAbsorptionMonochromatic
        K                ::AngularJ64
        omega            ::Float64
        multipole1       ::EmMultipole
        multipole2       ::EmMultipole
        gauge            ::EmGauge
        Jsym             ::LevelSymmetry
        amplitude        ::Complex{Float64}
    end

    """
    `struct  MultiPhotonDeExcitation.Line_2pAbsorptionMonochromatic`
        ... defines a type for a two-photon absorption line by monochromatic light that may include the definition of channels.

        + initialLevel     ::Level          ... initial-(state) level
        + finalLevel       ::Level          ... final-(state) level
        + omega            ::Float64        ... Energy of the incoming photons.
        + csLinear         ::EmProperty     ... Total cross section for linearly-polarized incident light.
        + csRightCircular  ::EmProperty     ... Total cross section for right-circularly polarized incident light.
        + csUnpolarized    ::EmProperty     ... Total cross section for unpolarized incident light.
        + channels         ::Array{MultiPhotonDeExcitation.Channel_2pAbsorptionMonochromatic,1}
                                            ... List of MultiPhotonDeExcitation.Channel_2pAbsorptionMonochromatic's of this line.
    """
    struct  Line_2pAbsorptionMonochromatic
        initialLevel       ::Level
        finalLevel         ::Level
        omega              ::Float64
        csLinear           ::EmProperty
        csRightCircular    ::EmProperty
        csUnpolarized      ::EmProperty
        channels           ::Array{MultiPhotonDeExcitation.Channel_2pAbsorptionMonochromatic,1}
    end

    # `Base.show(io::IO, line::MultiPhotonDeExcitation.Line_2pAbsorptionMonochromatic)`
    #   ... prepares a proper printout of the variable line::MultiPhotonDeExcitation.Line_2pAbsorptionMonochromatic.
    function Base.show(io::IO, line::MultiPhotonDeExcitation.Line_2pAbsorptionMonochromatic)
        println(io, "initialLevel:      $(line.initialLevel)  ")
        println(io, "finalLevel:        $(line.finalLevel)  ")
        println(io, "omega:             $(line.omega)  ")
        println(io, "csLinear:          $(line.csLinear)  ")
        println(io, "csRightCircular:   $(line.csRightCircular)  ")
        println(io, "csUnpolarized:     $(line.csUnpolarized)  ")
        println(io, "channels:          $(line.channels)  ")
    end

    """
    `MultiPhotonDeExcitation.computeAmplitudesProperties_2pAbsorptionMonochromatic(
                             line::MultiPhotonDeExcitation.Line_2pAbsorptionMonochromatic, grid::Radial.Grid,
                             settings::MultiPhotonDeExcitation.Settings)`
        ... to compute all amplitudes and properties of the given line; a line::MultiPhotonDeExcitation.Line_2pAbsorptionMonochromatic
            is returned for which the amplitudes and properties are now evaluated.
    """
    function  computeAmplitudesProperties_2pAbsorptionMonochromatic(line::MultiPhotonDeExcitation.Line_2pAbsorptionMonochromatic,
                                                                    grid::Radial.Grid, settings::MultiPhotonDeExcitation.Settings)
        newChannels = MultiPhotonDeExcitation.Channel_2pAbsorptionMonochromatic[]
        for channel in line.channels
            amplitude = MultiPhotonDeExcitation.computeReducedAmplitudeAbsorption(channel.K, line.finalLevel, channel.multipole2,
                                            channel.Jsym, channel.omega, channel.multipole1, line.initialLevel, channel.gauge, grid, settings.green)
            push!( newChannels, MultiPhotonDeExcitation.Channel_2pAbsorptionMonochromatic(channel.K, channel.omega, channel.multipole1, channel.multipole2,
                                                                                channel.gauge, channel.Jsym, amplitude) )
        end
        # Calculate the requested cross sections, etc.
        csLinear        = EmProperty(0., 0.)
        csRightCircular = EmProperty(0., 0.)
        csUnpolarized   = EmProperty(0., 0.)
        for property in settings.process.properties
            if      typeof(property) == TotalCsLinear
                totalCs_Cou = MultiPhotonDeExcitation.computeTotalCsLinear(line, EmGauge("Coulomb"), settings)
                totalCs_Bab = MultiPhotonDeExcitation.computeTotalCsLinear(line, EmGauge("Babushkin"), settings)
                csLinear    = EmProperty( totalCs_Cou,  totalCs_Bab)
            elseif  typeof(property) == TotalCsRightCircular
                totalCs_Cou = MultiPhotonDeExcitation.computeTotalCsRightCircular(line, EmGauge("Coulomb"), settings)
                totalCs_Bab = MultiPhotonDeExcitation.computeTotalCsRightCircular(line, EmGauge("Babushkin"), settings)
                csRightCircular = EmProperty( totalCs_Cou,  totalCs_Bab)
            elseif  typeof(property) == TotalCsUnpolarized
                totalCs_Cou = MultiPhotonDeExcitation.computeTotalCsUnpolarized(line, EmGauge("Coulomb"), settings)
                totalCs_Bab = MultiPhotonDeExcitation.computeTotalCsUnpolarized(line, EmGauge("Babushkin"), settings)
                csUnpolarized = EmProperty( totalCs_Cou,  totalCs_Bab)
            end
        end
        line = MultiPhotonDeExcitation.Line_2pAbsorptionMonochromatic( line.initialLevel, line.finalLevel, line.omega,
                                                                       csLinear, csRightCircular, csUnpolarized, newChannels)
        return( line )
    end

    """
    `MultiPhotonDeExcitation.computeLines(process::TwoPhotonAbsorptionMonochromatic, finalMultiplet::Multiplet,
                                          initialMultiplet::Multiplet, grid::Radial.Grid, settings::MultiPhotonDeExcitation.Settings;
                                          output=true)`
        ... to compute the multiphoton transition amplitudes and all properties as requested by the given settings. A list of
            lines::Array{MultiPhotonDeExcitation.Line_2pAbsorptionMonochromatic,1} is returned.
    """
    function  computeLines(process::TwoPhotonAbsorptionMonochromatic, finalMultiplet::Multiplet, initialMultiplet::Multiplet,
                           grid::Radial.Grid, settings::MultiPhotonDeExcitation.Settings; output=true)
        println("")
        printstyled("MultiPhotonDeExcitation.computeLines(::TwoPhotonAbsorptionMonochromatic): The computation of amplitudes starts now ... \n", color=:light_green)
        printstyled("---------------------------------------------------------------------------------------------------------------------- \n", color=:light_green)
        println("")
        #
        lines = MultiPhotonDeExcitation.determineLines_2pAbsorptionMonochromatic(finalMultiplet, initialMultiplet, settings)
        # Display all selected lines before the computations start
        if  settings.printBefore    MultiPhotonDeExcitation.displayLines_2pAbsorptionMonochromatic(lines)    end
        # Calculate all amplitudes and requested properties
        newLines = MultiPhotonDeExcitation.Line_2pAbsorptionMonochromatic[]
        for  line in lines
            newLine = MultiPhotonDeExcitation.computeAmplitudesProperties_2pAbsorptionMonochromatic(line, grid, settings)
            push!( newLines, newLine)
        end
        # Print all results to screen
        MultiPhotonDeExcitation.displayCrossSections_2pAbsorptionMonochromatic(stdout, settings.process.properties, lines)
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary    MultiPhotonDeExcitation.displayCrossSections_2pAbsorptionMonochromatic(iostream, settings.process.properties, lines)   end
        #
        if    output    return( lines )
        else            return( nothing )
        end
    end

    """
    `MultiPhotonDeExcitation.computeReducedAmplitudeAbsorption(K::AngularJ64, finalLevel::Level, multipole2::EmMultipole, Jsym::LevelSymmetry,
                                                                                omega::Float64, multipole1::EmMultipole, initialLevel::Level,
                                                                 gauge::EmGauge, grid::Radial.Grid, greenChannels::Array{AtomicState.GreenChannel,1})`
        ... to compute the reduced amplitude U^{K, 2gamma emission} (K, Jf, multipole2, Jsym, omega, multipole1, Ji) by means of the
            given Green function channels. An amplitude::Complex{Float64} is returned.
    """
    function computeReducedAmplitudeAbsorption(K::AngularJ64, finalLevel::Level, multipole2::EmMultipole, Jsym::LevelSymmetry,
                                                                omega::Float64, multipole1::EmMultipole, initialLevel::Level,
                                                                gauge::EmGauge, grid::Radial.Grid, greenChannels::Array{AtomicState.GreenChannel,1})
        U = Complex(0.);    found = false
        for channel in greenChannels
            if  Jsym == channel.symmetry
                found = true
                for (i, nuLevel) in enumerate(channel.gMultiplet.levels)
                    U = U + PhotoEmission.amplitude("absorption", multipole2, gauge, omega, finalLevel, nuLevel, grid, display=false) *
                            PhotoEmission.amplitude("absorption", multipole1, gauge, omega, nuLevel, initialLevel, grid, display=false) /
                            (initialLevel.energy + omega - nuLevel.energy)
                end
            end
        end

        if    found
              U = U * AngularMomentum.Wigner_6j(initialLevel.J, finalLevel.J, K, AngularJ64(multipole2.L), AngularJ64(multipole1.L), Jsym.J)
              ##x println(" ")
        else  println("No Green function cannel found for amplitude U^{K, 2gamma absorption} (K, Jf, multipole2, Jsym = $Jsym, omega, multipole1, Ji) ")
        end

        return( U )
    end

    """
    `MultiPhotonDeExcitation.computeTotalCsLinear(line::MultiPhotonDeExcitation.Line_2pAbsorptionMonochromatic,
                                                  gauge::EmGauge, settings::MultiPhotonDeExcitation.Settings)`
        ... to compute the total cross sections for linearly-polarized incident light. A tcs::Float64 is returned.
    """
    function computeTotalCsLinear(line::MultiPhotonDeExcitation.Line_2pAbsorptionMonochromatic,
                                  gauge::EmGauge, settings::MultiPhotonDeExcitation.Settings)
        tcs = 0.;   Klist = oplus(line.finalLevel.J, line.initialLevel.J);      omega = line.omega
        symi = LevelSymmetry(line.initialLevel.J, line.initialLevel.parity);    symf = LevelSymmetry(line.finalLevel.J, line.finalLevel.parity)

        for  K in Klist
            qList = AngularMomentum.m_values(K)
            for  q in qList
                amp = ComplexF64(0.)
                for  lambda1  in [-1, 1]
                    for  lambda2  in [-1, 1]
                        for  mp1 in settings.multipoles
                            for  mp2 in settings.multipoles
                                if   mp1.electric   p1 = 1    else    p1 = 0    end
                                if   mp2.electric   p2 = 1    else    p2 = 0    end
                                symmetries  = AngularMomentum.allowedTotalSymmetries(symf, mp2, mp1, symi)
                                for Jsym in symmetries
                                    wa = (1.0im)^(mp1.L - p1 + mp2.L - p2) * (-lambda1)^p1 * (-lambda2)^p2
                                    wb = sqrt( (2*mp1.L + 1)*(2*mp2.L + 1) ) * (2*Basics.twice(K) + 1)
                                    wc = AngularMomentum.Wigner_3j(mp1.L, mp2.L, K, lambda1, lambda2, q)
                                    wd = MultiPhotonDeExcitation.getReducedAmplitudeAbsorption(K, line.finalLevel, mp2, Jsym, omega, mp1,
                                                                                line.initialLevel, gauge, line.channels)
                                    ##x println("computeTotalCsLinear: L1, L2, K, lambda1, lambda2, q = $(mp1.L), $(mp2.L), $K, $lambda1, $lambda2, $q")
                                    ##x println("computeTotalCsLinear: wa, wb, wc, wd = $wa  $wb  $wc  $wd")

                                    amp = amp + (1.0im)^(mp1.L - p1 + mp2.L - p2) * (-lambda1)^p1 * (-lambda2)^p2           *
                                                sqrt( (2*mp1.L + 1)*(2*mp2.L + 1) ) * (2*Basics.twice(K) + 1)       *
                                                AngularMomentum.Wigner_3j(mp1.L, mp2.L, K, lambda1, lambda2, q)             *
                                                MultiPhotonDeExcitation.getReducedAmplitudeAbsorption(K, line.finalLevel, mp2, Jsym, omega, mp1,
                                                                                line.initialLevel, gauge, line.channels)
                                end
                            end
                        end
                    end
                end
                tcs = tcs + abs( amp )^2
            end
        end

        ## println("computeTotalCsLinear: tcs = $tcs")
        tcs = tcs * 8*pi^5 * Defaults.getDefaults("alpha")^2 / (Basics.twice(line.initialLevel.J) + 1) / omega^2

        return( tcs )
    end

    """
    `MultiPhotonDeExcitation.computeTotalCsRightCircular(line::MultiPhotonDeExcitation.Line_2pAbsorptionMonochromatic,
                                                         gauge::EmGauge, settings::MultiPhotonDeExcitation.Settings)`
        ... to compute the total cross sections for right-cicularly polarized incident light. A tcs::Float64 is returned.
    """
    function computeTotalCsRightCircular(line::MultiPhotonDeExcitation.Line_2pAbsorptionMonochromatic,
                                         gauge::EmGauge, settings::MultiPhotonDeExcitation.Settings)
        tcs = 0.;   Klist = oplus(line.finalLevel.J, line.initialLevel.J);      omega = line.omega
        symi = LevelSymmetry(line.initialLevel.J, line.initialLevel.parity);    symf = LevelSymmetry(line.finalLevel.J, line.finalLevel.parity)

        # Need to be filled

        tcs = tcs * 8*pi^5 * Defaults.getDefaults("alpha")^2 / (Basics.twice(line.initialLevel.J) + 1) / omega^2

        return( tcs )
    end

    """
    `MultiPhotonDeExcitation.computeTotalCsUnpolarized(line::MultiPhotonDeExcitation.Line_2pAbsorptionMonochromatic,
                                                       gauge::EmGauge, settings::MultiPhotonDeExcitation.Settings)`
        ... to compute the total cross sections for linearly-polarized incident light. A tcs::Float64 is returned.
    """
    function computeTotalCsUnpolarized(line::MultiPhotonDeExcitation.Line_2pAbsorptionMonochromatic,
                                       gauge::EmGauge, settings::MultiPhotonDeExcitation.Settings)
        tcs = 0.;   Klist = oplus(line.finalLevel.J, line.initialLevel.J);      omega = line.omega
        symi = LevelSymmetry(line.initialLevel.J, line.initialLevel.parity);    symf = LevelSymmetry(line.finalLevel.J, line.finalLevel.parity)

        for  K in Klist
            qList = AngularMomentum.m_values(K)
            for  q in qList
                for  lambda1  in [-1, 1]
                    for  lambda2  in [-1, 1]
                        amp = ComplexF64(0.)
                        for  mp1 in settings.multipoles
                            for  mp2 in settings.multipoles
                                if   mp1.electric   p1 = 1    else    p1 = 0    end
                                if   mp2.electric   p2 = 1    else    p2 = 0    end
                                symmetries  = AngularMomentum.allowedTotalSymmetries(symf, mp2, mp1, symi)
                                for Jsym in symmetries
                                    amp = amp + (1.0im)^(mp1.L - p1 + mp2.L - p2) * (-lambda1)^p1 * (-lambda2)^p2 *
                                                sqrt( (2*mp1.L + 1)*(2*mp2.L + 1) ) * (2*Basics.twice(K) + 1)       *
                                                AngularMomentum.Wigner_3j(mp1.L, mp2.L, K, lambda1, lambda2, q)             *
                                                MultiPhotonDeExcitation.getReducedAmplitudeAbsorption(K, line.finalLevel, mp2, Jsym, omega, mp1,
                                                                                line.initialLevel, gauge, line.channels)
                                end
                            end
                        end
                        tcs = tcs + abs( amp )^2
                    end
                end
            end
        end

        tcs = tcs * 8*pi^5 * Defaults.getDefaults("alpha")^2 / (Basics.twice(line.initialLevel.J) + 1) / omega^2

        return( tcs )
    end

    """
    `MultiPhotonDeExcitation.getReducedAmplitudeEmission(K::AngularJ64, finalLevel::Level, multipole2::EmMultipole, Jsym::LevelSymmetry, omega::Float64,
                                                                                multipole1::EmMultipole, initialLevel::Level,
                                                         gauge::EmGauge, channels::Array{MultiPhotonDeExcitation.Channel_2pAbsorptionMonochromatic,1})`
        ... to get/return the reduced amplitude U^{K, 2gamma emission} (K, Jf, multipole2, Jsym, omega, multipole1, Ji) from the calculated list
            of channels. An amplitude::Complex{Float64} is returned.
    """
    function getReducedAmplitudeAbsorption(K::AngularJ64, finalLevel::Level, multipole2::EmMultipole, Jsym::LevelSymmetry, omega::Float64,
                                                                             multipole1::EmMultipole, initialLevel::Level,
                                           gauge::EmGauge, channels::Array{MultiPhotonDeExcitation.Channel_2pAbsorptionMonochromatic,1})
        U = Complex(0.);    found = false
        for channel in channels
            if  K == channel.K  &&  omega == channel.omega  &&   Jsym == channel.Jsym       &&  multipole1 == channel.multipole1  &&
                multipole2 == channel.multipole2            &&   (gauge == channel.gauge  ||  EmGauge("Magnetic")  == channel.gauge)
                U = channel.amplitude;    found = true
            end
        end

        if    found
              println("U^{K, 2gamma absorption} (..) = $U   ** amplitude found. ")
        else  println("U^{$K, 2gamma absorption} (Jf, mp2=$multipole2, Jsym=$Jsym, omega=$omega, mp1=$multipole1, Ji)  ** NO amplitude found.")
        end

        return( U )
    end

    """
    `MultiPhotonDeExcitation.determineChannels_2pAbsorptionMonochromatic(omega::Float64, finalLevel::Level, initialLevel::Level,
                                                                         settings::MultiPhotonDeExcitation.Settings)`
        ... to determine a list of MultiPhotonDeExcitation.Channel_2pAbsorptionMonochromatic for a transitions from the initial to
            final level and by taking into account the particular settings of for this computation;
            an Array{MultiPhotonDeExcitation.Channel_2pAbsorptionMonochromatic,1} is returned.
    """
    function determineChannels_2pAbsorptionMonochromatic(omega::Float64, finalLevel::Level, initialLevel::Level, settings::MultiPhotonDeExcitation.Settings)
        channels   = MultiPhotonDeExcitation.Channel_2pAbsorptionMonochromatic[];
        symi       = LevelSymmetry(initialLevel.J, initialLevel.parity);    symf = LevelSymmetry(finalLevel.J, finalLevel.parity)
        for  mp1 in settings.multipoles
            for  mp2 in settings.multipoles
                symmetries  = AngularMomentum.allowedTotalSymmetries(symf, mp2, mp1, symi)
                Klist       = oplus(symf.J, symi.J)
                for  symn in symmetries
                    for  gauge in settings.gauges
                        # Include further restrictions if appropriate
                        if     string(mp1)[1] == 'E' || string(mp2)[1] == 'E'  &&   gauge == Basics.UseCoulomb
                            for K in Klist  push!(channels, MultiPhotonDeExcitation.Channel_2pAbsorptionMonochromatic(K, omega, mp1, mp2, Basics.Coulomb, symn, 0.) )     end
                        elseif string(mp1)[1] == 'E' || string(mp2)[1] == 'E'  &&   gauge == Basics.UseBabushkin
                            for K in Klist  push!(channels, MultiPhotonDeExcitation.Channel_2pAbsorptionMonochromatic(K, omega, mp1, mp2, Basics.Babushkin, symn, 0.) )   end
                        elseif string(mp1)[1] == 'M' && string(mp2)[1] == 'M'
                            for K in Klist  push!(channels, MultiPhotonDeExcitation.Channel_2pAbsorptionMonochromatic(K, omega, mp1, mp2, Basics.Magnetic, symn, 0.) )    end
                        end
                    end
                end
            end
        end

        return( channels )
    end

    """
    `MultiPhotonDeExcitation.determineLines_2pAbsorptionMonochromatic(finalMultiplet::Multiplet, initialMultiplet::Multiplet,
                                                                      settings::MultiPhotonDeExcitation.Settings)`
        ... to determine a list of MultiPhotonDeExcitation.Line_2pAbsorptionMonochromatic's for transitions between the levels from the given
            initial- and final-state multiplets and by taking into account the particular selections and settings for this computation;
            an Array{MultiPhotonDeExcitation.Line_2pAbsorptionMonochromatic,1} is returned. Apart from the level specification, all physical
            properties are set to zero during this initialization process.
    """
    function  determineLines_2pAbsorptionMonochromatic(finalMultiplet::Multiplet, initialMultiplet::Multiplet,
                                                       settings::MultiPhotonDeExcitation.Settings)
        lines = MultiPhotonDeExcitation.Line_2pAbsorptionMonochromatic[]
        for  iLevel  in  initialMultiplet.levels
            for  fLevel  in  finalMultiplet.levels
                if  Basics.selectLevelPair(iLevel, fLevel, settings.lineSelection)
                    omega    = (fLevel.energy - iLevel.energy) / 2.
                    channels = MultiPhotonDeExcitation.determineChannels_2pAbsorptionMonochromatic(omega, fLevel, iLevel, settings)
                    push!( lines, MultiPhotonDeExcitation.Line_2pAbsorptionMonochromatic(iLevel, fLevel, omega,
                                                                     EmProperty(0., 0.), EmProperty(0., 0.), EmProperty(0., 0.), channels) )
                end
            end
        end
        return( lines )
    end

    """
    `MultiPhotonDeExcitation.displayLines_2pAbsorptionMonochromatic(lines::Array{MultiPhotonDeExcitation.Line_2pAbsorptionMonochromatic,1})`
        ... to display a list of lines and channels that have been selected due to the prior settings. A neat table of all selected
            transitions and energies is printed but nothing is returned otherwise.
    """
    function  displayLines_2pAbsorptionMonochromatic(lines::Array{MultiPhotonDeExcitation.Line_2pAbsorptionMonochromatic,1})
        nx = 175
        println(" ")
        println("  Selected two-photon absorption lines by monochromatic and equally-polarized photons:")
        println(" ")
        println("  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(18, "i-level-f"; na=0);                         sb = sb * TableStrings.hBlank(18)
        sa = sa * TableStrings.center(18, "i--J^P--f"; na=3);                         sb = sb * TableStrings.hBlank(21)
        sa = sa * TableStrings.center(12, "Energy"; na=2);
        sb = sb * TableStrings.center(12, TableStrings.inUnits("energy"); na=2)
        sa = sa * TableStrings.center(12, "omega"; na=4);
        sb = sb * TableStrings.center(12, TableStrings.inUnits("energy"); na=4)
        sa = sa * TableStrings.flushleft(90, "List of multipoles & intermediate level symmetries"; na=4)
        sb = sb * TableStrings.flushleft(90, "(K-rank, multipole_1, Jsym, multipole_2, gauge), ..."; na=4)
        println(sa);    println(sb);    println("  ", TableStrings.hLine(nx))
        #
        for  line in lines
            sa  = "";    isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                         fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
            sa = sa * TableStrings.center(18, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=2)
            sa = sa * TableStrings.center(18, TableStrings.symmetries_if(isym, fsym); na=4)
            energy = line.finalLevel.energy - line.initialLevel.energy
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", energy))      * "    "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", line.omega))  * "     "
            mpGaugeList = Tuple{AngularJ64, Basics.EmMultipole, LevelSymmetry, Basics.EmMultipole, Basics.EmGauge}[]
            for  channel in  line.channels
                push!( mpGaugeList, (channel.K, channel.multipole1, channel.Jsym, channel.multipole2, channel.gauge) )
            end
            ##x println("displayLines_2pAbsorptionMonochromatic: mpGaugeList = ", mpGaugeList)
            wa = TableStrings.twoPhotonGaugeTupels(105, mpGaugeList)
            if  length(wa) > 0    sb = sa * wa[1];    println( sb )    end
            for  i = 2:length(wa)
                sb = TableStrings.hBlank( length(sa) );    sb = sb * wa[i];    println( sb )
            end
        end
        println("  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

    """
    `MultiPhotonDeExcitation.displayCrossSections_2pAbsorptionMonochromatic(stream::IO,
                                                  properties::Array{MultiPhotonDeExcitation.AbstractMultiPhotonProperty,1},
                                                  lines::Array{MultiPhotonDeExcitation.Line_2pAbsorptionMonochromatic,1})`
        ... to display all results, energies, rates, etc. of the selected lines. A neat table is printed but nothing is
            returned otherwise.
    """
    function  displayCrossSections_2pAbsorptionMonochromatic(stream::IO, properties::Array{AbstractMultiPhotonProperty,1},
                                                             lines::Array{Line_2pAbsorptionMonochromatic,1})
        nx = 75
        println(stream, " ")
        println(stream, "  Two-photon absorption by monochromatic and equally-polarized photons (usually from the same beam):")
        println(stream, " ")
        println(stream, "  Cross sections are given for:")
        noCs = 0  # Number of cross sections to be printed
        for property in properties
            if      typeof(property) == TotalCsLinear
                noCs = noCs + 1;   println(stream, "    + total cross sections for linearly-polarized incident light ($noCs)")
            elseif  typeof(property) == TotalCsRightCircular
                noCs = noCs + 1;   println(stream, "    + total cross sections for right-circularly polarized incident light ($noCs)")
            elseif  typeof(property) == TotalCsUnpolarized
                noCs = noCs + 1;   println(stream, "    + total cross sections for unpolarized incident light ($noCs)")
            end
        end
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx + 34noCs))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(18, "i-level-f"; na=2);                         sb = sb * TableStrings.hBlank(20)
        sa = sa * TableStrings.center(18, "i--J^P--f"; na=4);                         sb = sb * TableStrings.hBlank(22)
        sa = sa * TableStrings.center(10, "Energy"; na=4);
        sb = sb * TableStrings.center(10, TableStrings.inUnits("energy"); na=4)
        sa = sa * TableStrings.center(10, "omega";  na=4);
        sb = sb * TableStrings.center(10, TableStrings.inUnits("energy"); na=7)
        for no = 1:noCs
            sa = sa * TableStrings.center(28, "Cou -- cross section ($no) -- Bab"; na=4);
            sb = sb * TableStrings.center(28, TableStrings.inUnits("cross section") * "          " * TableStrings.inUnits("cross section"); na=8)
        end

        println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx + 34noCs))
        #
        for  line in lines
            sa  = "";    isym = LevelSymmetry( line.initialLevel.J, line.initialLevel.parity)
                         fsym = LevelSymmetry( line.finalLevel.J,   line.finalLevel.parity)
            sa = sa * TableStrings.center(18, TableStrings.levels_if(line.initialLevel.index, line.finalLevel.index); na=4)
            sa = sa * TableStrings.center(18, TableStrings.symmetries_if(isym, fsym); na=4)
            energy = line.finalLevel.energy - line.initialLevel.energy
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", energy))         * "    "
            sa = sa * @sprintf("%.4e", Defaults.convertUnits("energy: from atomic", line.omega))     * "        "
            for property in properties
                if      typeof(property) == TotalCsLinear
                    sa = sa * @sprintf("%.4e", Defaults.convertUnits("cross section: from atomic",   line.csLinear.Coulomb))          * "      "
                    sa = sa * @sprintf("%.4e", Defaults.convertUnits("cross section: from atomic",   line.csLinear.Babushkin))        * "          "
                elseif  typeof(property) == TotalCsRightCircular
                    sa = sa * @sprintf("%.4e", Defaults.convertUnits("cross section: from atomic",   line.csRightCircular.Coulomb))   * "      "
                    sa = sa * @sprintf("%.4e", Defaults.convertUnits("cross section: from atomic",   line.csRightCircular.Babushkin)) * "          "
                elseif  typeof(property) == TotalCsUnpolarized
                    sa = sa * @sprintf("%.4e", Defaults.convertUnits("cross section: from atomic",   line.csUnpolarized.Coulomb))     * "      "
                    sa = sa * @sprintf("%.4e", Defaults.convertUnits("cross section: from atomic",   line.csUnpolarized.Babushkin))   * "          "
                end
            end
            println(stream, sa )
        end
        println(stream, "  ", TableStrings.hLine(nx + 34noCs))
        #
        return( nothing )
    end


text:
#=
signomial minimization problem adapted from
"Signomial and Polynomial Optimization via Relative Entropy and Partial Dualization"
(2019) by Murray, Chandrasekaran, & Wierman

a signomial f = Sig(A, c) of variables
    x in R^n
with parameters
    c in R^m
    A in R^{m, n}
takes values
    f(x) = sum_{i in [m]} c_i exp(A_{i, :} * x)

for signomials f, (g_p)_{p in 1:q}, the signomial minimization problem is
    inf_{x in R^n} f(x) : g_p(x) >= 0, p in 1:q

the convex SAGE relaxation problem is
    max_{γ in R, μ in R_+^q} γ : c - γ - sum_{p in 1:q} μ_p * g_p in C_SAGE(A)

a vector d in R^m belongs to SAGE cone C_SAGE(A) if
    exists C in R^{m, m} :
    d = sum_{k in [m]} C_{k, :}
    C_{k, :} in C_AGE(A, k)
which is equivalent to the feasibility problem over C in R^{m, m}, V in R^{m, m-1}
    d = sum_{k in [m]} C_{k, :}
    for k in [m]:
        (A_{\k, :} - [1,...,1] * A_{k, :}')' * V_{k, :} == [0,...,0]
        [C_{k, k} + sum(V_{k, :}), C_{k, \k}, V_{k, :}] in RelEntr(1 + 2(m - 1))
=#

include(joinpath(@__DIR__, "data.jl"))

struct SignomialMinJuMP{T <: Real} <: ExampleInstanceJuMP{T}
    fc::Vector
    fA::AbstractMatrix
    gc::Vector
    gA::Vector
    obj_ub::Real
end

SignomialMinJuMP{Float64}(sig_name::Symbol) =
    SignomialMinJuMP{Float64}(signomialmin_data[sig_name]...)

SignomialMinJuMP{Float64}(m::Int, n::Int) =
    SignomialMinJuMP{Float64}(signomialmin_random(m, n)...)

function build(inst::SignomialMinJuMP{T}) where {T <: Float64}
    (fc, fA, gc, gA) = (inst.fc, inst.fA, inst.gc, inst.gA)
    (fm, n) = size(fA)
    @assert length(fc) == fm
    q = length(gc)
    @assert length(gA) == q
    @assert all(size(gA_p, 2) == n for gA_p in gA)
    @assert all(size(gA_p, 1) == length(gc_p) for (gc_p, gA_p) in zip(gc, gA))

    # find unique terms
    unique_terms = Dict{Vector, Tuple{Vector{Int}, Vector{Tuple{Int, Int}}}}()
    for k in 1:size(fA, 1)
        row_k = fA[k, :]
        if row_k in keys(unique_terms)
            push!(unique_terms[row_k][1], k)
        else
            unique_terms[row_k] = (Int[k], Tuple{Int, Int}[])
        end
    end
    for (p, gA_p) in enumerate(gA)
        for k in 1:size(gA_p, 1)
            p_row_k = gA_p[k, :]
            if p_row_k in keys(unique_terms)
                push!(unique_terms[p_row_k][2], (p, k))
            else
                unique_terms[p_row_k] = (Int[], Tuple{Int, Int}[(p, k)])
            end
        end
    end
    A = vcat((row_k' for row_k in keys(unique_terms))...)
    m = size(A, 1)

    model = JuMP.Model()
    JuMP.@variable(model, γ)
    JuMP.@objective(model, Max, γ)
    JuMP.@variable(model, μ[1:q] >= 0)

    d = zeros(JuMP.AffExpr, m)
    const_found = false
    for k in 1:m
        row_k = A[k, :]
        (ks, pks) = unique_terms[row_k]
        for k2 in ks
            d[k] += fc[k2]
        end
        for (p, k2) in pks
            d[k] -= μ[p] * gc[p][k2]
        end
        if iszero(norm(row_k))
            # row is a constant term
            d[k] -= γ
            @assert !const_found
            const_found = true
        end
    end
    if !const_found
        A = vcat(A, zeros(A, 1, n))
        d = vcat(d, -γ)
    end

    # setup SAGE constraints
    notk = [[l for l in 1:m if l != k] for k in 1:m]
    JuMP.@variable(model, C[1:m, 1:m])
    JuMP.@variable(model, V[1:m, 1:(m - 1)])
    JuMP.@constraint(model, [k in 1:m], d[k] == sum(C[:, k]))
    JuMP.@constraint(model, [k in 1:m, i in 1:n],
        dot(A[notk[k], i] .- A[k, i], V[k, :]) == 0)
    JuMP.@constraint(model, [k in 1:m], vcat(C[k, k] +
        sum(V[k, :]), C[k, notk[k]], V[k, :]) in MOI.RelativeEntropyCone(2m - 1))

    return model
end

function test_extra(inst::SignomialMinJuMP{T}, model::JuMP.Model) where T
    stat = JuMP.termination_status(model)
    @test stat == MOI.OPTIMAL
    if (stat == MOI.OPTIMAL) && !isnan(inst.obj_ub)
        # check objective value is correct
        tol = eps(T)^0.2
        @test JuMP.objective_value(model) <= inst.obj_ub + tol
    end
end


text:
export Vortex

"""
$(TYPEDEF)

Defines a point vortex with x-position `x`, y-position `y`, and strength `Γ`.

# Fields

$(TYPEDFIELDS)
"""
mutable struct Vortex
    """x: x-coordinate of the vortex position."""
    x::Float64
    """y: y-coordinate of the vortex position."""
    y::Float64
    """Γ: Strength of the vortex. Positive if counter-clockwise."""
    Γ::Float64
end

"""
$(TYPEDSIGNATURES)

Returns the strengths of all point vortices in the `StructArray` `vl` as `ScalarData`.
"""
function getstrengths(vl::StructArray{Vortex})
    return ScalarData(vl.Γ)
end

"""
$(TYPEDSIGNATURES)

Returns the positions of all point vortices in the `StructArray` `vl` as `VectorData`.
"""
function getpositions(vl::StructArray{Vortex})
    return VectorData(vl.x,vl.y)
end

"""
$(TYPEDSIGNATURES)

Sets the `x` and `y` fields of the point vortices in the `StructArray` `vl` to the entries of `xnew` and `ynew`, respectively.
"""
function setpositions!(vl::StructArray{Vortex}, xnew, ynew)
    vl.x .= xnew
    vl.y .= ynew
end

"""
$(TYPEDSIGNATURES)

Sets the `Γ` field of the point vortices in the `StructArray` `vl` to the entries of `Γnew`.
"""
function setstrengths!(vl::StructArray{Vortex}, Γnew, idx=nothing)
    if isnothing(idx)
        vl.Γ .= Γnew
    else
        setindex!(vl.Γ, Γnew, idx)
    end
end


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

JLLWrappers.@generate_main_file_header("libevdev")
JLLWrappers.@generate_main_file("libevdev", UUID("2db6ffa8-e38f-5e21-84af-90c45d0032cc"))
end  # module libevdev_jll


text:
dir = @__DIR__

function genMarkdownFiles(scriptName::String)
    open(joinpath(dir, "$(splitext(scriptName)[1]).md"), "w") do outf
        redirect_stdout(outf) do
            toFile = """
            @def title = "$(scriptName)"
            @def hascode = true
            @def tags = ["syntax", "code"]
            ## $(scriptName)
            \\fig{/_assets/scripts/output/$(scriptName).png}
            \\input{julia}{/_assets/scripts/output/$(scriptName).jl}
            ~~~
            <span style="color:#e53e00;"> <strong> Dependencies $(scriptName) </strong> </span>
            ~~~
            \\prettyshow{/_assets/scripts/$(scriptName)}
            """
            print(toFile)
        end
    end
end
genMarkdownFiles("testmdgen")

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
using JuMP, Xpress
using Random
using Compat.Test

# Based on original model implemented in Mosel, courtesy of Truls Flatberg
function wordhunt(Words::Array{String,1}, D = [:E,:S,:W, :N, :SE, :NE], Gridsize=7,woptimizer=with_optimizer(Xpress.Optimizer,MIPRELSTOP=0.1),printres=true)

    Maxlength = maximum(length.(Words))
    M = 1:Gridsize
    N = uppercase.(Words)
    L = unique(join(N,""))

    model = Model(woptimizer)

    @variable(model, x[M,M,L],Bin)  # Letter in location M,M
    @variable(model, y[N,M,M,D])    # Word in location M,M in direction D

    # Allowed start positions for each direction in D
    for n in N, i in M, j in M, d in D
        if (d == :E) && (j > Gridsize - length(n) + 1)
            @constraint(model, y[n,i,j,d] <=0)
        elseif (d == :S) && (i > Gridsize - length(n) + 1 )
            @constraint(model, y[n,i,j,d] <=0)
        elseif (d == :W) && (j < length(n))
            @constraint(model, y[n,i,j,d] <=0)
        elseif (d == :N) && (i < length(n))
            @constraint(model, y[n,i,j,d] <=0)
        elseif (d == :SE) && ((i > Gridsize - length(n) + 1) || (j > Gridsize - length(n) + 1))
            @constraint(model, y[n,i,j,d] <=0)
	    elseif (d == :NE) && ((j > Gridsize - length(n) + 1) || (i < length(n) ))
            @constraint(model, y[n,i,j,d] <=0)
        else
            set_binary(y[n,i,j,d])
        end
    end

    # Each position has maximum one letter
    for i in M, j in M
	    @constraint(model,sum( x[i,j,l] for l in L) <= 1)
    end

    # Each word is allowed one position and direction (if inserted)
    for n in N
	    @constraint(model,sum( y[n,i,j,d] for i in M, j in M, d in D) <= 1)
    end

    # Placement of words
    for n in N, i in M, j in M, d in D
	    ii = i
	    jj = j
        for l in 1:length(n)
            if  ii in M && jj in M
                @constraint(model, x[ii,jj,n[l]] >= y[n,i,j,d] )
        	    if d == :E
		    	    jj = jj + 1
                elseif d == :S
	    	 	    ii = ii + 1
                elseif d == :SE
			        ii = ii + 1
			        jj = jj + 1
                elseif d == :NE
           	        ii = ii -1
           	        jj = jj +1
                elseif d == :W
           	        jj = jj - 1
                elseif d == :N
           	        ii = ii - 1
                end
            end
        end
    end

    # Objective function
    most_words = sum( 2 * Maxlength * y[n,i,j,d] for n in N, i in M, j in M, d in D)
    least_letters = sum( x[i,j,l] for i in M, j in M, l in L)

    P  = [:W,:N,:SE]
    if length(intersect(P,D)) > 0
        pref = sum( y[n,i,j,d] for n in N, i in M, j in M, d in intersect(P,D))
    else
        pref = 0
    end

    @objective(model, Max, most_words - least_letters + 0.1 * pref )

    optimize!(model)

    function printSol(x,L,fillrand=true)
        letters = "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
        for i in 1:size(x,1)
            for j in 1:size(x,2)
                gotsol = false
                for l in L
                    if JuMP.value(x[i,j,l])>0.5
                        print(l)
                        gotsol = true
                    end
                end
                if !gotsol
                    if fillrand
                        print(letters[rand(1:end)])
                    else
                        print(' ')
                    end
                end
                print(' ')
            end
            println()
        end
    end

    if printres
        # Show status and objective value
        print(JuMP.termination_status(model),"\t")
        println(JuMP.objective_value(model))
        println()

        # Print solution to console
        printSol(x,L,false)
        println()
        printSol(x,L)
    end

    return JuMP.termination_status(model)
end

const Words = ["Ada","Claire","Hugo","Idris","Julia","Karel","Mary","Max","Maya","Miranda"]

# Demo
# @time wordhunt(Words,[:SE,:NE]);

# Tests:
@testset "Solver parameters" begin
    # Easy problem, should solve sub second
    @test wordhunt(Words,[:E],4,with_optimizer(Xpress.Optimizer,MAXTIME=15),false) == MOI.OPTIMAL
    # Check time limit
    @test wordhunt(Words,[:E,:S,:SE,:NE],7,with_optimizer(Xpress.Optimizer,MAXTIME=1),false) == MOI.TIME_LIMIT
    # Check miprelstop
    @test wordhunt(Words,[:E,:S,:SE,:NE],7,with_optimizer(Xpress.Optimizer,MIPRELSTOP=.9),false) == MOI.OPTIMAL
    # Test combinations of miprelstop and maxtime
    @test wordhunt(Words,[:E,:W,:S,:SE,:NE],7,with_optimizer(Xpress.Optimizer,MIPRELSTOP=.01,MAXTIME=1),false) == MOI.TIME_LIMIT
    @test wordhunt(Words,[:E,:W,:S,:SE,:NE],7,with_optimizer(Xpress.Optimizer,MIPRELSTOP=.9,MAXTIME=15),false) == MOI.OPTIMAL
end