text:
# x_ij is a keypoint corresponding to the i-th 3D point
# observed from the j-th camera (viewpoint).
# Assume we have a keypoint array X = {x_ij} observed in condition
# n_points = 4, n_viewpoints = 3,
# where x_22 x_31 x_43 are missing.
#
# X            [x_11 x_12 x_13 x_21 x_23 x_32 x_33 x_41 x_42]
# indices          1    2    3    4    5    6    7    8    9
# point index      1    1    1    2    2    3    3    4    4
# viewpoint index  1    2    3    1    3    2    3    1    2

# viewpoints_by_point[1] == [1, 2, 3]  # x_11 x_12 x_13
# viewpoints_by_point[2] == [4, 5]     # x_21 x_23
# viewpoints_by_point[3] == [6, 7]     # x_32 x_33
# viewpoints_by_point[4] == [8, 9]     # x_41 x_42

# points_by_viewpoint[1] == [1, 4, 8]  # x_11 x_21 x_41
# points_by_viewpoint[2] == [2, 6, 9]  # x_12 x_32 x_42
# points_by_viewpoint[3] == [3, 5, 7]  # x_13 x_23 x_33
#
# mask = [
#     1 1 1;  # x_11 x_12 x_13
#     1 0 1;  # x_21      x_23
#     0 1 1;  #      x_32 x_33
#     1 1 0;  # x_41 x_42
# ]

struct Indices
    mask::BitArray
    viewpoints_by_point::Vector{Vector{Int}}
    points_by_viewpoint::Vector{Vector{Int}}
end

n_points(indices::Indices) = length(indices.viewpoints_by_point)
n_viewpoints(indices::Indices) = length(indices.points_by_viewpoint)

function Indices(viewpoint_indices::Array{Int}, point_indices::Array{Int})
    @assert length(viewpoint_indices) == length(point_indices)

    n_points = maximum(point_indices)
    n_viewpoints = maximum(viewpoint_indices)
    mask = BitArray(undef, n_points, n_viewpoints)

    viewpoints_by_point = [Array{Int}([]) for i in 1:n_points]
    points_by_viewpoint = [Array{Int}([]) for j in 1:n_viewpoints]

    unique_points = Set([])
    unique_viewpoints = Set([])

    for (index, (i, j)) in enumerate(zip(point_indices, viewpoint_indices))
        push!(viewpoints_by_point[i], index)
        push!(points_by_viewpoint[j], index)
        mask[i, j] = 1

        push!(unique_points, i)
        push!(unique_viewpoints, j)
    end

    # they cannot be true if some point / viewpoint indices are missing
    # ex. raises AssertionError if n_viewpoints == 4 and
    # unique_viewpoints == [1, 2, 4]  (3 is missing)
    @assert length(unique_viewpoints) == n_viewpoints
    @assert length(unique_points) == n_points

    Indices(mask, viewpoints_by_point, points_by_viewpoint)
end

function points_by_viewpoint(indices::Indices, j::Int)
    """
    'points_by_viewpoint(j)' should return indices of 3D points
    observable from a viewpoint j
    """

    indices.points_by_viewpoint[j]
end

function viewpoints_by_point(indices::Indices, i::Int)
    """
    'viewpoints_by_point(i)' should return indices of viewpoints
    that can observe a point i
    """

    indices.viewpoints_by_point[i]
end

function shared_point_indices(indices::Indices, j::Int, k::Int)
    """
    j, k: viewpoint indices
    This function returns two indices of points commonly observed from both viewpoints.
    These two indices are corresponding to the first and second view respectively
    """
    mask_j = view(indices.mask, :, j)
    mask_k = view(indices.mask, :, k)

    indices_j = Array{Int}([])
    indices_k = Array{Int}([])

    index_j = 0
    index_k = 0
    for (bit_j, bit_k) in zip(mask_j, mask_k)
        if bit_j == 1
            index_j += 1
        end

        if bit_k == 1
            index_k += 1
        end

        if bit_j & bit_k == 1
            push!(indices_j, index_j)
            push!(indices_k, index_k)
        end
    end

    if length(indices_j) == 0  # (== length(indices_k))
        return nothing  # no shared points found between j and k
    end

    (indices.points_by_viewpoint[j][indices_j],
     indices.points_by_viewpoint[k][indices_k])
end


text:
module Xva

import ..Clarus

function fva(;params...)
  return Clarus.api_request("XVA","FVA";params...)
end

function mva(;params...)
  return Clarus.api_request("XVA","FVA";params...)
end

function mvaattribution(;params...)
  return Clarus.api_request("XVA","MVAAttribution";params...)
end

function sensitivity(;params...)
  return Clarus.api_request("XVA","Sensitivity";params...)
end

end #MODULE-END

ext: jl


text:
"""
# script test_parameterVar

- Julia version: 1.4
- Author: niallcullinane
- Date: 2020-08-20

Script to vary parameter of LAMINART and plot layers and activation.
"""

using DrWatson
@quickactivate "MEngProject"
using MEngProject,
    CUDA,
    DifferentialEquations,
    PyPlot,
    NNlib,
    ImageFiltering,
    Images,
    ImageIO,
    MEngProject.LaminartKernels,
    MEngProject.LaminartInitFunc,
    MEngProject.Utils,
    BenchmarkTools,
    Test

using OrdinaryDiffEq,
    ParameterizedFunctions, LSODA, Sundials, DiffEqDevTools, Noise

batch = 1001

files = ["kan_sq_cont_l.png"]

tspan = (0.0f0, 800f0)

batch_ = string(batch, "_", rand(1000:9999))
mkdir(plotsdir(string("paraVar", batch_)))

for file in files
    para_sets = [
        (
            C_1 = 1.5f0,
            C_2 = 0.075f0,
            ϕ = 2.0f0,
            Γ = 0.2f0,
            η_p = 2.1f0,
            η_m = 1.5f0,
            λ = 1.5f0,
            ψ = 0.5f0,
        ), #base
        (
            C_1 = 0.0f0,
            C_2 = 0.0f0,
            ϕ = 2.0f0,
            Γ = 0.2f0,
            η_p = 2.1f0,
            η_m = 1.5f0,
            λ = 1.5f0,
            ψ = 0.5f0,
        ), #C_1, C_2 down full
        (
            C_1 = 1.5f0,
            C_2 = 1.0f0,
            ϕ = 2.0f0,
            Γ = 0.2f0,
            η_p = 2.1f0,
            η_m = 1.5f0,
            λ = 1.5f0,
            ψ = 0.5f0,
        ), #C_2 up
        (
            C_1 = 2.5f0,
            C_2 = 0.075f0,
            ϕ = 2.0f0,
            Γ = 0.2f0,
            η_p = 2.1f0,
            η_m = 1.5f0,
            λ = 1.5f0,
            ψ = 0.5f0,
        ), #C_1 up
        (
            C_1 = 1.5f0,
            C_2 = 0.075f0,
            ϕ = 0.0f0,
            Γ = 0.2f0,
            η_p = 2.1f0,
            η_m = 1.5f0,
            λ = 1.5f0,
            ψ = 0.5f0,
        ), #ϕ down full
        (
            C_1 = 1.5f0,
            C_2 = 0.075f0,
            ϕ = 4.0f0,
            Γ = 0.2f0,
            η_p = 2.1f0,
            η_m = 1.5f0,
            λ = 1.5f0,
            ψ = 0.5f0,
        ), #ϕ up double
        (
            C_1 = 1.5f0,
            C_2 = 0.075f0,
            ϕ = 2.0f0,
            Γ = 0.1f0,
            η_p = 2.1f0,
            η_m = 1.5f0,
            λ = 1.5f0,
            ψ = 0.5f0,
        ), #Γ down half
        (
            C_1 = 1.5f0,
            C_2 = 0.075f0,
            ϕ = 2.0f0,
            Γ = 0.4f0,
            η_p = 2.1f0,
            η_m = 1.5f0,
            λ = 1.5f0,
            ψ = 0.5f0,
        ), #Γ up double
        (
            C_1 = 1.5f0,
            C_2 = 0.075f0,
            ϕ = 2.0f0,
            Γ = 0.2f0,
            η_p = 3.0f0,
            η_m = 0.5f0,
            λ = 1.5f0,
            ψ = 0.5f0,
        ), #η_p up, η_m down
        (
            C_1 = 1.5f0,
            C_2 = 0.075f0,
            ϕ = 2.0f0,
            Γ = 0.2f0,
            η_p = 1.0f0,
            η_m = 2.5f0,
            λ = 1.5f0,
            ψ = 0.5f0,
        ), #η_p down, n_m up
        (
            C_1 = 1.5f0,
            C_2 = 0.075f0,
            ϕ = 2.0f0,
            Γ = 0.2f0,
            η_p = 2.1f0,
            η_m = 1.5f0,
            λ = 2.5f0,
            ψ = 0.5f0,
        ), #λ up
        (
            C_1 = 1.5f0,
            C_2 = 0.075f0,
            ϕ = 2.0f0,
            Γ = 0.2f0,
            η_p = 2.1f0,
            η_m = 1.5f0,
            λ = 0.5f0,
            ψ = 0.5f0,
        ), #λ down
        (
            C_1 = 1.5f0,
            C_2 = 0.075f0,
            ϕ = 2.0f0,
            Γ = 0.2f0,
            η_p = 2.1f0,
            η_m = 1.5f0,
            λ = 1.5f0,
            ψ = 1.0f0,
        ), #ψ up
        (
            C_1 = 1.5f0,
            C_2 = 0.075f0,
            ϕ = 2.0f0,
            Γ = 0.2f0,
            η_p = 2.1f0,
            η_m = 1.5f0,
            λ = 1.5f0,
            ψ = 0.0f0,
        ), #ψ down full
        #         (
        #             C_1 = 0.5f0,
        #             C_2 = 0.037f0,
        #             ϕ = 1.0f0,
        #             Γ = 0.2f0,
        #             η_p = 2.1f0,
        #             η_m = 1.5f0,
        #             λ = 1.5f0,
        #             ψ = 1.5f0,
        #         ), #all fb down half
        #         (
        #             C_1 = 0.0f0,
        #             C_2 = 0.0f0,
        #             ϕ = 0.0f0,
        #             Γ = 0.2f0,
        #             η_p = 2.1f0,
        #             η_m = 1.5f0,
        #             λ = 1.5f0,
        #             ψ = 1.5f0,
        #         ), #all fb down full
        #         (
        #             C_1 = 3.0f0,
        #             C_2 = 0.125f0,
        #             ϕ = 4.0f0,
        #             Γ = 0.2f0,
        #             η_p = 2.1f0,
        #             η_m = 1.5f0,
        #             λ = 1.5f0,
        #             ψ = 1.5f0,
        #         ), #all fb up double
    ]
    test_name = [
        "base",
        "C1C2_0",
        "C2_1",
        "C1_25",
        "phi0",
        "phi4",
        "Gamma01",
        "Gamma04",
        "etaP3_etaM05",
        "etaP1_etaM25",
        "lamda05",
        "lamda15",
        "psi1",
        "psi0",
        #         "fb_half",
        #         "fb0",
        #         "fb_doub",
    ]

    test_name_plt = [
        "Base",
        "\$C_1, C_2 =0\$",
        "\$C_2=1\$",
        "\$C_1= 2.5\$",
        "\$ϕ =0\$",
        "\$ϕ = 4\$",
        "\$Γ= 0.101\$",
        "\$Γ=0.4\$",
        "\$η^+ =3,η^-=0.5\$",
        "\$η^+=1, \$η^-=2.5\$",
        "\$λ=0.5\$",
        "\$λ=1.5\$",
        "\$ψ = 1\$",
        "\$ψ=0\$",
        #         "All feedback halved",
        #         "All feedback \$=0\$",
        #         "All feedback doubled",
    ]

    for para_test ∈ enumerate(para_sets)
        try
            p = LaminartInitFunc.parameterInit_conv_gpu(
                datadir("img", file),
                Parameters.para_var(para_test[2]),
            )

            u0 = cu(reshape(
                zeros(Float32, p.dim_i, p.dim_j * (5 * p.K + 2)),
                p.dim_i,
                p.dim_j,
                5 * p.K + 2,
                1,
            ))

            arr1 = similar(@view u0[:, :, 1:2, :])
            arr2 = similar(@view u0[:, :, 1:1, :])

            f = LaminartFunc.LamFunction(
                arr1, #x
                similar(arr1), #m
                similar(arr1), #s
                arr2, #x_lgn,
                similar(arr1), #C,
                similar(arr1), #H_z,
                similar(arr1), # dy_temp,
                similar(arr1), # dm_temp,
                similar(arr1), # dz_temp,
                similar(arr1), # ds_temp,
                similar(arr2), # dv_temp,
                similar(arr1), # H_z_temp,
                similar(arr2), #  V_temp_1,
                similar(arr2), #  V_temp_2,
                similar(arr1), #  Q_temp,
                similar(arr1), #   P_temp
            )

            prob = ODEProblem(f, u0, tspan, p)
            # 	@benchmark sol = solve(prob)
            sol = solve(prob)

            for t ∈ [25, 50, 100, 200, 400, 800]
                # 				for t ∈ [25,50,100]
                @inbounds begin
                    v0 = @view sol(t)[:, :, :, 1]
                    axMax = findmax(v0)[1]

                    for k ∈ 1:2:10
                        fig, ax = plt.subplots()

                        v1 = @view sol(t)[:, :, k, 1]
                        v2 = @view sol(t)[:, :, k+1, 1]
                        im = ax.imshow(
                            v1,
                            cmap = matplotlib.cm.PRGn,
                            vmax = axMax,
                            vmin = -axMax,
                        )
                        im2 = ax.imshow(
                            v2,
                            cmap = matplotlib.cm.RdBu_r,
                            vmax = axMax,
                            vmin = -axMax,
                            alpha = 0.5,
                        )

                        cbar = fig.colorbar(im2, shrink = 0.9, ax = ax)
                        cbar.ax.set_xlabel("\$k=2\$")
                        cbar = fig.colorbar(im, shrink = 0.9, ax = ax)
                        cbar.set_alpha(0.5)
                        cbar.draw_all()
                        cbar.ax.set_xlabel("\$k=1\$")
                        layer = Utils.layers[k]
                        plt.title(string(
                            "Layer: $layer, \$t=$t\$, ",
                            test_name_plt[para_test[1]],
                        ))
                        plt.axis("off")
                        fig.tight_layout()
                        plt.savefig(plotsdir(
                            string("paraVar", batch_),
                            string(
                                file,
                                "_para_",
                                test_name[para_test[1]],
                                "_t",
                                t,
                                "_",
                                Utils.la[k],
                                ".png",
                            ),
                        ))
                        close("all")
                    end
                end

                k = 11
                fig, ax = plt.subplots()
                v1 = @view sol[:, :, k, 1, t]
                v2 = @view sol[:, :, k+1, 1, t]
                im = ax.imshow(
                    v1,
                    cmap = matplotlib.cm.PRGn,
                    vmax = axMax,
                    vmin = -axMax,
                )
                im2 = ax.imshow(
                    v2,
                    cmap = matplotlib.cm.RdBu_r,
                    vmax = axMax,
                    vmin = -axMax,
                    alpha = 0.5,
                )

                cbar = fig.colorbar(im2, shrink = 0.9, ax = ax)
                cbar.ax.set_xlabel("\$v^-\$")
                cbar = fig.colorbar(im, shrink = 0.9, ax = ax)
                cbar.set_alpha(0.5)
                cbar.draw_all()
                cbar.ax.set_xlabel("\$v^+\$")

                layer = Utils.layers[k]
                plt.title(string(
                    "Layer: $layer, \$t=$t\$, ",
                    test_name_plt[para_test[1]],
                ))
                plt.axis("off")
                fig.tight_layout()

                plt.savefig(plotsdir(
                    string("paraVar", batch_),
                    string(
                        file,
                        "_para_",
                        test_name[para_test[1]],
                        "_t",
                        t,
                        "_",
                        Utils.la[k],
                        ".png",
                    ),
                ))
            end
        catch
        finally
            close("all")
        end

        # time plot
        fig, axs = plt.subplots()
        @inbounds begin
            for k ∈ 1:12
                v3 = @view sol[:, :, k, 1, end]
                v4 = @view sol[findmax(v3)[2][1], findmax(v3)[2][2], k, 1, :]
                layer = Utils.layers_1[k]
                axs.plot(
                    sol.t,
                    v4,
                    Utils.lines[k],
                    label = "$layer",
                    alpha = 0.8,
                )
            end
            axs.set_xlabel("Time")
            axs.set_ylabel("Activation")
            plt.title(test_name_plt[para_test[1]])
            plt.legend()
            fig.tight_layout()
            plt.savefig(plotsdir(
                string("paraVar", batch_),
                string(file, "_para_", test_name[para_test[1]], "_time.png"),
            ))
            close("all")
            u0 = nothing
            p = nothing
            arr1 = nothing
            arr2 = nothing
            f = nothing
            prob = nothing
            sol = nothing
            CUDA.reclaim()
        end
    end
end


text:
@testset "W447" begin
    msgs = lintstr("a = Dict{:Symbol, Any}")
    @test messageset(msgs) == Set([:W447])
    @test occursin("type parameter for Dict", msgs[1].message)

    msgs = lintstr("a = Dict{:Symbol, Any}()")
    @test messageset(msgs) == Set([:W447])
    @test occursin("type parameter for Dict", msgs[1].message)

    @test isempty(lintstr("a = Set{Tuple{Int, Int}}()"))

    msgs = lintstr("""
    b = :Symbol
    a = Dict{b, Any}()
    """)
    @test messageset(msgs) == Set([:W447])
    @test occursin("type parameter for Dict", msgs[1].message)

    msgs = lintstr("a = Array{2, Int64}()")
    @test messageset(msgs) == Set([:W447])
    @test occursin("type parameter for Array", msgs[1].message)
end

@testset "Curly" begin
    s = """
    a = Set{(Int, Int)}()
    """
    msgs = lintstr(s)
    @test msgs[1].code == :W441
    @test occursin("probably illegal use inside curly", msgs[1].message)

    s = """
    a = Array{Int64, 5, 5}()
    """
    msgs = lintstr(s)
    @test msgs[1].code == :W446
    @test occursin("too many type parameters", msgs[1].message)

    s = """
    a = Ptr{Cvoid}
    """
    msgs = lintstr(s)
    @test isempty(msgs)

    s = """
    @traitfn ft1{X,Y; Cmp{X,Y}}(x::X, y::Y) = x > y ? 5 : 6
    """
    msgs = lintstr(s)
    @test isempty(msgs)
end


text:
@testset "intermodes" begin
    original_image = testimage("lena")
    for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64})
        img = T.(original_image)
        img₀₁ = binarize(Intermodes(), img)

        # Check original image is unchanged.
        @test img == T.(testimage("lena"))

        # Check that the image only has ones or zeros.
        non_zeros = findall(x -> x != 0.0 && x != 1.0, img₀₁)
        @test length(non_zeros) == 0

        # Check type of binarized image.
        @test typeof(img₀₁) == Array{Gray{Bool},2}

        # Check that ones and zeros have been assigned to the correct side of the threshold.
        maxval, maxpos = findmax(Gray.(img))
        @test img₀₁[maxpos] == 1
        minval, minpos = findmin(Gray.(img))
        @test img₀₁[minpos] == 0
    end
end


text:
using ADI
using BenchmarkTools
using CSV
using DataFrames
using HCIDatasets: BetaPictoris, HR8799
using InvertedIndices
using PyCall

results = []

# julia benchmarks

julia_algs = (Median(), PCA(20), NMF(20))

for dataset in (BetaPictoris, HR8799)
    @info "Benchmarking - $dataset ADI Reduction"
    cube = dataset[:cube]
    angles = dataset[:pa]
    for (alg, name) in zip(julia_algs, ("median", "pca_20", "nmf_20"))
        (alg isa NMF && dataset === HR8799) && continue
        @info "Julia - $alg"

        time_elapsed = @belapsed $alg($cube, $angles)

        @info "Julia - $alg" time=time_elapsed
        push!(results, (framework="ADI.jl", alg=name, N=length(cube), time=time_elapsed))
    end
end

# python benchmarks

vip = pyimport("vip_hci")
py_algs = (
    (args...; kwargs...) -> vip.medsub.median_sub(args..., verbose=false, kwargs...),
    (args...; kwargs...) -> vip.pca.pca(args...; ncomp=20, verbose=false, kwargs...),
    (args...; kwargs...) -> vip.nmf.nmf(args...; ncomp=20, verbose=false, kwargs...),
)

for dataset in (BetaPictoris, HR8799)
    @info "Benchmarking - $dataset ADI Reduction"
    cube = dataset[:cube]
    angles = dataset[:pa]
    for (alg, name) in zip(py_algs, ("median", "pca_20", "nmf_20"))
        (name == "nmf_20" && dataset === HR8799) && continue
        @info "Python - $alg"

        time_elapsed = @belapsed $alg($cube, $angles)

        @info "Python - $alg" time=time_elapsed
        push!(results, (framework="VIP", alg=name, N=length(cube), time=time_elapsed))
    end
end

# load results file and update
path = joinpath(@__DIR__, "adi_benchmarks.csv")
df = isfile(path) ? CSV.File(path) |> DataFrame : DataFrame()
out = append!(DataFrame(results), df)
unique!(out, Not(:time)) |> CSV.write(path)


text:
module Mandelbrot

using Base.Threads
using Plots
using ProgressMeter
using ArrayFire

include("fractal_struct.jl")
include("computation_routines.jl")
include("visual_utils.jl")

export FractalData,
       computeMandelbrot!,
       move_center!,
       move_up!,
       move_down!,
       move_left!,
       move_right!,
       display_fractal,
       preview_fractal

end # end module


text:
"""
Class for variational Gaussian Processes models (non-sparse)

```julia
 VGP(X::AbstractArray{T1,N1},y::AbstractArray{T2,N2},
     kernel::Union{Kernel,AbstractVector{<:Kernel}},
     likelihood::LikelihoodType,inference::InferenceType;
     verbose::Integer=0,Autotuning::Bool=true,
     atfrequency::Integer=1,IndependentPriors::Bool=true,
     ArrayType::UnionAll=Vector)
```

Argument list :

**Mandatory arguments**

 - `X` : input features, should be a matrix N×D where N is the number of observation and D the number of dimension
 - `y` : input labels, can be either a vector of labels for multiclass and single output or a matrix for multi-outputs (note that only one likelihood can be applied)
 - `kernel` : covariance function, can be either a single kernel or a collection of kernels for multiclass and multi-outputs models
 - `likelihood` : likelihood of the model, currently implemented : Gaussian, Bernoulli (with logistic link), Multiclass (softmax or logistic-softmax) see [`Likelihood Types`](@ref likelihood_user)
 - `inference` : inference for the model, can be analytic, numerical or by sampling, check the model documentation to know what is available for your likelihood see the [`Compatibility Table`](@ref compat_table)

**Keyword arguments**

 - `verbose` : How much does the model print (0:nothing, 1:very basic, 2:medium, 3:everything)
 - `Autotuning` : Flag for optimizing hyperparameters
 - `atfrequency` : Choose how many variational parameters iterations are between hyperparameters optimization
 - `IndependentPriors` : Flag for setting independent or shared parameters among latent GPs
 - `ArrayType` : Option for using different type of array for storage (allow for GPU usage)
"""
mutable struct VGP{L<:Likelihood,I<:Inference,T<:Real,V<:AbstractVector{T}} <: AbstractGP{L,I,T,V}
    X::Matrix{T} #Feature vectors
    y::LatentArray #Output (-1,1 for classification, real for regression, matrix for multiclass)
    nSample::Int64 # Number of data points
    nDim::Int64 # Number of covariates per data point
    nFeature::Int64 # Number of features of the GP (equal to number of points)
    nLatent::Int64 # Number pf latent GPs
    IndependentPriors::Bool # Use of separate priors for each latent GP
    nPrior::Int64 # Equal to 1 or nLatent given IndependentPriors
    μ::LatentArray{V}
    Σ::LatentArray{Symmetric{T,Matrix{T}}}
    η₁::LatentArray{V}
    η₂::LatentArray{Symmetric{T,Matrix{T}}}
    μ₀::LatentArray{MeanPrior{T}}
    Knn::LatentArray{Symmetric{T,Matrix{T}}}
    invKnn::LatentArray{Symmetric{T,Matrix{T}}}
    kernel::LatentArray{Kernel{T}}
    likelihood::Likelihood{T}
    inference::Inference{T}
    verbose::Int64 #Level of printing information
    Autotuning::Bool
    atfrequency::Int64
    Trained::Bool
end

function VGP(X::AbstractArray{T1,N1},y::AbstractArray{T2,N2},kernel::Union{Kernel,AbstractVector{<:Kernel}},
            likelihood::LikelihoodType,inference::InferenceType;
            verbose::Integer=0,Autotuning::Bool=true,atfrequency::Integer=1,mean::Union{<:Real,AbstractVector{<:Real},MeanPrior}=ZeroMean(),
            IndependentPriors::Bool=true,ArrayType::UnionAll=Vector) where {T1<:Real,T2,N1,N2,LikelihoodType<:Likelihood,InferenceType<:Inference}

            X,y,nLatent,likelihood = check_data!(X,y,likelihood)
            @assert check_implementation(:VGP,likelihood,inference) "The $likelihood is not compatible or implemented with the $inference"

            nPrior = IndependentPriors ? nLatent : 1
            nFeature = nSample = size(X,1); nDim = size(X,2);
            kernel = ArrayType([deepcopy(kernel) for _ in 1:nPrior])

            μ = LatentArray([zeros(T1,nFeature) for _ in 1:nLatent]); η₁ = deepcopy(μ)
            Σ = LatentArray([Symmetric(Matrix(Diagonal(one(T1)*I,nFeature))) for _ in 1:nLatent]);
            η₂ = -0.5*inv.(Σ);
            μ₀ = []
            if typeof(mean) <: Real
                μ₀ = [ConstantMean(mean) for _ in 1:nPrior]
            elseif typeof(mean) <: AbstractVector{<:Real}
                μ₀ = [EmpiricalMean(mean) for _ in 1:nPrior]
            else
                μ₀ = [mean for _ in 1:nPrior]
            end
            Knn = LatentArray([deepcopy(Σ[1]) for _ in 1:nPrior]);
            invKnn = copy(Knn)

            likelihood = init_likelihood(likelihood,inference,nLatent,nSample)
            inference = init_inference(inference,nLatent,nSample,nSample,nSample)

            VGP{LikelihoodType,InferenceType,T1,ArrayType{T1}}(X,y,
                    nFeature, nDim, nFeature, nLatent,
                    IndependentPriors,nPrior,μ,Σ,η₁,η₂,
                    μ₀,Knn,invKnn,kernel,likelihood,inference,
                    verbose,Autotuning,atfrequency,false)
end

function Base.show(io::IO,model::VGP{<:Likelihood,<:Inference,T}) where T
    print(io,"Variational Gaussian Process with a $(model.likelihood) infered by $(model.inference) ")
end


text:
## Basic tests for simple utilities within BB
using BinaryBuilder, Test, Pkg, UUIDs
using BinaryBuilder: preferred_runner, resolve_jlls, CompilerShard, preferred_libgfortran_version, preferred_cxxstring_abi, gcc_version, available_gcc_builds, getversion, generate_compiler_wrappers!, getpkg, build_project_dict
using BinaryBuilder.BinaryBuilderBase
using BinaryBuilder.Wizard

@testset "File Collection" begin
    temp_prefix() do prefix
        # Create a file and a link, ensure that only the one file is returned by collect_files()
        f = joinpath(prefix, "foo")
        f_link = joinpath(prefix, "foo_link")
        touch(f)
        symlink(f, f_link)
        d = joinpath(prefix, "bar")
        d_link = joinpath(prefix, "bar_link")
        mkpath(d)
        symlink(d, d_link)

        files = collect_files(prefix)
        @test length(files) == 3
        @test realpath(f) in files
        @test realpath(f_link) in files
        @test !(realpath(d) in files)
        @test !(realpath(d_link) in files)
        @test d_link in files

        collapsed_files = collapse_symlinks(files)
        @test length(collapsed_files) == 1
        @test realpath(f) in collapsed_files

        files = collect_files(prefix, exclude_dirs = false)
        @test length(files) == 4
        @test realpath(f) in files
        @test realpath(f_link) in files
        @test realpath(d) in files
        @test realpath(d_link) in files

        files = collect_files(prefix, islink)
        @test length(files) == 2
        @test !(realpath(f) in files)
        @test f_link in files
        @test !(realpath(d) in files)
        @test d_link in files

        files = collect_files(prefix, islink, exclude_dirs = false)
        @test length(files) == 2
        @test !(realpath(f) in files)
        @test f_link in files
        @test !(realpath(d) in files)
        @test d_link in files
    end
end

@testset "environment and history saving" begin
    mktempdir() do temp_path
        # This is a litmus test, to catch any errors before we do a `@test_throws`
        @test_logs (:error, r"^Unable to find valid license file") match_mode=:any autobuild(
            temp_path,
            "this_will_pass",
            v"1.0.0",
            # No sources to speak of
            FileSource[],
            # Just exit with code 0
            """
            exit 0
            """,
            # Build for this platform
            [platform],
            # No products
            Product[],
            # No depenedencies
            Dependency[],
        )

        @test_throws ErrorException autobuild(
            temp_path,
            "this_will_fail",
            v"1.0.0",
            FileSource[],
            # Simple script that just sets an environment variable
            """
            MARKER=1
            exit 1
            """,
            [platform],
            Product[],
            Dependency[],
        )

        # build_path is the nonce'd build directory
        build_path = joinpath(temp_path, "build", triplet(platform))
        build_path = joinpath(build_path, first(readdir(build_path)))

        # Ensure that we get a metadir, and that our history and .env files are in there!
        metadir = joinpath(build_path, "metadir")
        @test isdir(metadir)

        hist_file = joinpath(metadir, ".bash_history")
        env_file = joinpath(metadir, ".env")
        @test isfile(hist_file)
        @test isfile(env_file)

        # Test that exit 1 is in .bash_history
        @test occursin("\nexit 1\n", read(open(hist_file), String))

        # Test that MARKER=1 is in .env:
        @test occursin("\nMARKER=1\n", read(open(env_file), String))

        # Delete the build path
        rm(build_path, recursive = true)
    end
end

@testset "Wizard Utilities" begin
    # Make sure canonicalization does what we expect
    zmq_url = "https://github.com/zeromq/zeromq3-x/releases/download/v3.2.5/zeromq-3.2.5.tar.gz"
    @test Wizard.canonicalize_source_url(zmq_url) == zmq_url
    this_url = "https://github.com/JuliaPackaging/BinaryBuilder.jl/blob/1fee900486baedfce66ddb24872133ef36b9d899/test/wizard.jl"
    this_url_ans = "https://raw.githubusercontent.com/JuliaPackaging/BinaryBuilder.jl/1fee900486baedfce66ddb24872133ef36b9d899/test/wizard.jl"
    @test Wizard.canonicalize_file_url(this_url) == this_url_ans

    # Make sure normalization does what we expect
    @test Wizard.normalize_name("foo/libfoo.tar.gz") == "libfoo"
    @test Wizard.normalize_name("foo/libfoo-2.dll") == "libfoo"
    @test Wizard.normalize_name("libfoo") == "libfoo"

    # with_gitcreds
    local creds_outer = nothing
    Wizard.with_gitcreds("user", "password") do creds
        @test creds isa LibGit2.UserPasswordCredential
        @test hasproperty(creds, :user)
        @test hasproperty(creds, :pass)
        creds_outer = creds # assign to parent scope, so that we can check on it later
        @test creds.user == "user"
        @test String(read(creds.pass)) == "password"
        @test !Base.isshredded(creds.pass)
    end
    @test creds_outer isa LibGit2.UserPasswordCredential
    @test creds_outer.user == ""
    @test Base.isshredded(creds_outer.pass)
    @test eof(creds_outer.pass)
    # in case it throws:
    creds_outer = nothing
    @test_throws ErrorException Wizard.with_gitcreds("user", "password") do creds
        creds_outer = creds
        error("...")
    end
    @test creds_outer isa LibGit2.UserPasswordCredential
    @test creds_outer.user == ""
    @test Base.isshredded(creds_outer.pass)
    @test eof(creds_outer.pass)
end

@testset "State serialization" begin
    state = Wizard.WizardState()
    state.step = :step34
    state.platforms = [Platform("x86_64", "linux")]
    state.source_urls = ["http://127.0.0.1:14444/a/source.tar.gz"]
    state.source_files = [BinaryBuilder.SetupSource{ArchiveSource}("/tmp/source.tar.gz", bytes2hex(sha256("a")), "")]
    state.name = "libfoo"
    state.version = v"1.0.0"
    state.dependencies = [Dependency(PackageSpec(;name="Zlib_jll")),
                          Dependency(PackageSpec(;name="CompilerSupportLibraries_jll"))]
    state.history = "exit 1"

    io = Dict()
    Wizard.serialize(io, state)
    new_state = Wizard.unserialize(io)

    for field in fieldnames(Wizard.WizardState)
        @test getfield(state, field) == getfield(new_state, field)
    end
end

# Test that updating Yggdrasil works
@testset "Yggdrasil" begin
    Core.eval(Wizard, :(yggdrasil_updated = false))
    @test_logs (:info, r"Yggdrasil") Wizard.get_yggdrasil()
end

@testset "Registration utils" begin
    name = "CGAL"
    version = v"1"
    dependencies = [Dependency("boost_jll"), Dependency("GMP_jll"),
                    Dependency("MPFR_jll"), Dependency("Zlib_jll")]
    dict = build_project_dict(name, version, dependencies)
    @test dict["name"] == "$(name)_jll"
    @test dict["version"] == "1.0.0"
    @test dict["uuid"] == "8fcd9439-76b0-55f4-a525-bad0597c05d8"
    @test dict["compat"] == Dict{String,Any}("julia" => "1.0", "JLLWrappers" => "1.2.0")
    @test all(in.(
        (
            "Pkg"       => "44cfe95a-1eb2-52ea-b672-e2afdf69b78f",
            "Libdl"     => "8f399da3-3557-5675-b5ff-fb832c97cbdb",
            "GMP_jll"   => "781609d7-10c4-51f6-84f2-b8444358ff6d",
            "MPFR_jll"  => "3a97d323-0669-5f0c-9066-3539efd106a3",
            "Zlib_jll"  => "83775a58-1f1d-513f-b197-d71354ab007a",
            "boost_jll" => "28df3c45-c428-5900-9ff8-a3135698ca75",
        ), Ref(dict["deps"])))
    project = Pkg.Types.Project(dict)
    @test project.name == "$(name)_jll"
    @test project.uuid == UUID("8fcd9439-76b0-55f4-a525-bad0597c05d8")
    # Make sure that a `BuildDependency` can't make it to the list of
    # dependencies of the new JLL package
    @test_throws MethodError build_project_dict(name, version, [BuildDependency("Foo_jll")])

    version = v"1.6.8"
    next_version = BinaryBuilder.get_next_wrapper_version("Xorg_libX11", version)
    @test next_version.major == version.major
    @test next_version.minor == version.minor
    @test next_version.patch == version.patch

    # Ensure passing a Julia dependency bound works
    dict = build_project_dict(name, version, dependencies, "1.4")
    @test dict["compat"] == Dict{String,Any}("julia" => "1.4", "JLLWrappers" => "1.2.0")

    dict = build_project_dict(name, version, dependencies, "~1.4")
    @test dict["compat"] == Dict{String,Any}("julia" => "~1.4", "JLLWrappers" => "1.2.0")

    @test_throws ErrorException build_project_dict(name, version, dependencies, "nonsense")

    # Ensure passing compat bounds works
    dependencies = [
        Dependency(PackageSpec(name="libLLVM_jll"), compat="=9.0.0"),
    ]
    dict = build_project_dict("Clang", v"9.0.1+2", dependencies)
    @test dict["compat"]["julia"] == "1.0"
    @test dict["compat"]["libLLVM_jll"] == "=9.0.0"

    dependencies = [
        Dependency(PackageSpec(name="libLLVM_jll"), compat="8.3 - 10"),
    ]
    dict = build_project_dict("Clang", v"9.0.1+2", dependencies)
    @test dict["compat"]["julia"] == "1.0"
    @test dict["compat"]["libLLVM_jll"] == "8.3 - 10"

    dependencies = [
        Dependency(PackageSpec(name="libLLVM_jll"), compat="8.3"),
    ]
    dict = build_project_dict("Clang", v"9.0.1+2", dependencies)
    @test dict["compat"]["libLLVM_jll"] == "8.3"
end


text:
"""

General Framework for the data augmented Gaussian Processes

"""
module AugmentedGaussianProcesses

const AGP = AugmentedGaussianProcesses
export AGP
export AbstractGPModel, GP, VGP, SVGP, VStP, MCGP, MOVGP, MOSVGP, MOARGP, OnlineSVGP # All models
export AbstractLikelihood, RegressionLikelihood, BernoulliLikelihood, MultiClassLikelihood # All categories of likelihoods
export GaussianLikelihood, StudentTLikelihood, LaplaceLikelihood, HeteroscedasticLikelihood # Regression Likelihoods
export LogisticLikelihood, BayesianSVM # Classification Likelihoods
export SoftMaxLikelihood, LogisticSoftMaxLikelihood # Multiclass Classification Likelihoods
export PoissonLikelihood, NegBinomialLikelihood # Event Likelihoods
export AbstractInference, Analytic, AnalyticVI, AnalyticSVI # Inference objects
export GibbsSampling, HMCSampling # Sampling inference
export NumericalVI,
    NumericalSVI, MCIntegrationVI, MCIntegrationSVI, QuadratureVI, QuadratureSVI # Numerical inference
export PriorMean, ZeroMean, ConstantMean, EmpiricalMean, AffineMean # Prior means
#Useful functions
export train!, sample
export predict_f, predict_y, proba_y
export fstar
export ELBO
export covariance, diag_covariance, prior_mean
export @augmodel

#General modules
using Reexport
@reexport using KernelFunctions
@reexport using Optimisers
@reexport using InducingPoints

using AbstractMCMC: AbstractMCMC, step, sample
# using AdvancedHMC
using ChainRulesCore: ChainRulesCore, NoTangent
using Distributions:
    Distributions,
    Distribution,
    dim,
    cov,
    mean,
    var,
    pdf,
    logpdf,
    loglikelihood,
    Normal,
    Poisson,
    NegativeBinomial,
    InverseGamma,
    Laplace,
    MvNormal,
    Gamma
using FastGaussQuadrature: gausshermite
using ForwardDiff
using GPLikelihoods: GPLikelihoods, AbstractLink, AbstractLikelihood
using GPLikelihoods:
    BernoulliLikelihood, PoissonLikelihood, HeteroscedasticGaussianLikelihood
using GPLikelihoods: LogisticLink, SoftMaxLink
using KernelFunctions: ColVecs, RowVecs
using LinearAlgebra
using ProgressMeter
using Random
using StatsBase
using SimpleTraits
using StatsFuns
using SpecialFunctions
using Zygote

#Include custom module for additional distributions
include("ComplementaryDistributions/ComplementaryDistributions.jl")
using .ComplementaryDistributions

# Main classes
abstract type AbstractInference{T<:Real} end
abstract type VariationalInference{T} <: AbstractInference{T} end
abstract type SamplingInference{T} <: AbstractInference{T} end
abstract type AbstractLatent{T<:Real,Tpr,Tpo} end

include("mean/priormean.jl")
include("data/datacontainer.jl")
include("functions/utils.jl")

# Models
include("models/AbstractGP.jl")
include("gpblocks/latentgp.jl")
include("models/GP.jl")
include("models/VGP.jl")
include("models/MCGP.jl")
include("models/SVGP.jl")
include("models/VStP.jl")
include("models/MOSVGP.jl")
include("models/MOVGP.jl")
include("models/OnlineSVGP.jl")
include("models/single_and_multi_output_utils.jl")

include("inference/inference.jl")
include("likelihood/likelihood.jl")
include("likelihood/generic_likelihood.jl")

include("functions/KLdivergences.jl")
include("functions/ELBO.jl")
include("data/utils.jl")
include("functions/plotting.jl")

# Training and prediction functions
include("training/states.jl")
include("training/training.jl")
include("training/sampling.jl")
include("training/onlinetraining.jl")
include("hyperparameter/autotuning.jl")
include("training/predictions.jl")
include("ar_predict.jl")

end #End Module


text:
abstract type IsingModel end

mutable struct IsingModel2D <: IsingModel
  h::Float64
  J::Float64
  T::Float64
  n::Int64
  m::Int64
  data::AbstractArray
  rng::AbstractRNG
  function IsingModel2D(h::Float64,J::Float64,T::Float64,n::Int64,m::Int64;seed::UInt64=rand(UInt64))
    data = zeros(Int64,n,m)
    rng = MersenneTwister(seed)
    new(h,J,T,n,m,data,rng)
  end
end
IsingModel2D(h::Float64,J::Float64,T::Float64,n::Int64;seed::UInt64=rand(UInt64)) = IsingModel2D(h,J,T,n,n;seed=seed)
IsingModel2D(n::Int64,m::Int64;seed::UInt64=rand(UInt64)) = IsingModel2D(1.0,1.0,1.0,n,m;seed=seed)
IsingModel2D(n::Int64;seed::UInt64=rand(UInt64)) = IsingModel2D(1.0,1.0,1.0,n,n;seed=seed)

function randomdata!(model::IsingModel)
  model.data = rand(Int64[-1,1],(model.n,model.m))
end

function allup!(model::IsingModel)
  model.data = ones(Int64,(model.n,model.m))
end

function alldown!(model::IsingModel)
  model.data = -ones(Int64,(model.n,model.m))
end

function flipbit(ind::Int64,model::IsingModel)
  return model.data[ind] == 1 ? -1 : 1
end

function flipbit!(ind::Int64,model::IsingModel)
  model.data[ind] == 1 ? model.data[ind] = -1 : model.data[ind] = 1
end

function localediff(ind::Int64,model::IsingModel)
  fbit::Int64 = flipbit(ind,model) - model.data[ind]
  # nearest-neighbor interactions
  x::Int64 = ind%model.m
  x = x != 0 ? x : model.m
  y::Int64 = ceil(Int64,ind/model.n)
  nn::Int64 = 0
  # check x-direction w/ pbc
  nn += (x+1) <= model.n ? model.data[x+1,y] : model.data[1,y]
  nn += (x-1) >= 1 ? model.data[x-1,y] : model.data[model.n,y]
  # check y-direction w/ pbc
  nn += (y+1) <= model.m ? model.data[x,y+1] : model.data[x,1]
  nn += (y-1) >= 1 ? model.data[x,y-1] : model.data[x,model.m]
  return -(model.h*fbit + model.J*fbit*nn)
end


text:
function count_unambiguous_digits(text_list, digit_list, pattern_dict)

    digit_counts = Dict((i, 0) for i in digit_list)
#    targets = Dict((length(pattern_dict[i]),i) for i in digit_list)
    # text list is a list of space separated strings
    for line in text_list
        lengths = [length(i) for i in split(line)]
        for digit in keys(digit_counts)
            digit_counts[digit] += count(x->x==length(pattern_dict[digit]), lengths)
        end
    end
    return digit_counts
end

patterns = Dict("0"=>"abcefg",
                "1"=>"cf",
                "2"=>"acdeg",
                "3"=>"acdfg",
                "4"=>"bcdf",
                "5"=>"abdfg",
                "6"=>"abdefg",
                "7"=>"acf",
                "8"=>"abcdefg",
                "9"=>"abcdfg")

pattern_length = Dict((i, length(patterns[i])) for i in keys(patterns))

input_lines = readlines("input.txt")

signal_patterns = map(line->split(line, "|")[1], input_lines)
outputs = map(line->split(line, "|")[2], input_lines)

count_dict = count_unambiguous_digits(outputs, ["1","4","7","8"], patterns)
println("Count dict ", count_dict)
answer = sum(values(count_dict))
println("total is ",answer)


text:
module ExitBlock

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Exit Block" => Ahorn.EntityPlacement(
        Maple.ExitBlock,
        "rectangle",
        Dict{String, Any}(),
        Ahorn.tileEntityFinalizer
    ),
)

Ahorn.editingOptions(entity::Maple.ExitBlock) = Dict{String, Any}(
    "tileType" => Ahorn.tiletypeEditingOptions()
)

Ahorn.minimumSize(entity::Maple.ExitBlock) = 8, 8
Ahorn.resizable(entity::Maple.ExitBlock) = true, true

Ahorn.selection(entity::Maple.ExitBlock) = Ahorn.getEntityRectangle(entity)

Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.ExitBlock, room::Maple.Room) = Ahorn.drawTileEntity(ctx, room, entity, alpha=0.7)

end


text:
struct Kpca
    X::Array{Float64}
    Kt::Array{Float64}
    T::Array{Float64}
    P::Array{Float64}
    sv::Vector{Float64}
    eig::Vector{Float64}
    D::Array{Float64}
    DKt::Array{Float64}
    vtot::Array{Float64}
    weights::Vector{Float64}
    kern
    dots
end

"""
    kpca(X, Y, weights = ones(size(X, 1)); nlv , kern = "krbf", kwargs...)
Kernel PCA  (Scholkopf et al. 1997, Scholkopf & Smola 2002, Tipping 2001).

* `X` : X-data.
* `weights` : vector (n,).
* `nlv` : Nb. principal components (PCs), or collection of nb. PCs, to consider.
* `kern` : Type of kernel used to compute the Gram matrices.
    Possible values are "krbf" of "kpol" (see respective functions `krbf` and `kpol`).
* `kwargs` : Named arguments to pass in the kernel function.

The method is implemented by SVD factorization of the weighted Gram matrix
D^(1/2) * Phi(`X`) * Phi(`X`)' * D^(1/2), where D is a diagonal matrix of weights for
the observations (rows of X).

The kernel Gram matrices are internally centered.

## References

Scholkopf, B., Smola, A., MÃ¼ller, K.-R., 1997. Kernel principal component analysis,
in: Gerstner, W., Germond, A., Hasler, M., Nicoud, J.-D. (Eds.), Artificial Neural Networks,
ICANN 97, Lecture Notes in Computer Science. Springer, Berlin, Heidelberg,
pp. 583-588. https://doi.org/10.1007/BFb0020217

Scholkopf, B., Smola, A.J., 2002. Learning with kernels: support vector machines, regularization,
optimization, and beyond, Adaptive computation and machine learning. MIT Press, Cambridge, Mass.

Tipping, M.E., 2001. Sparse kernel principal component analysis. Advances in neural information
processing systems, MIT Press. http://papers.nips.cc/paper/1791-sparse-kernel-principal-component-analysis.pdf

"""
function kpca(X, weights = ones(size(X, 1)); nlv, kern = "krbf", kwargs...)
    X = ensure_mat(X)
    n = size(X, 1)
    nlv = min(nlv, n)
    weights = mweights(weights)
    fkern = eval(Meta.parse(kern))
    K = fkern(X, X; kwargs...)     # In the future: fkern!(K, X, X; kwargs...)
    D = Diagonal(weights)
    Kt = K'
    DKt = D * Kt
    vtot = sum(DKt, dims = 1)
    Kc = K .- vtot' .- vtot .+ sum(D * DKt')
    sqrtD = sqrt.(D)
    Kd = sqrtD * Kc * sqrtD
    res = LinearAlgebra.svd(Kd)
    U = res.V[:, 1:nlv]
    eig = res.S
    eig[eig .< 0] .= 0
    sv = sqrt.(eig)
    P = sqrtD * scale(U, sv[1:nlv])     # In the future: scale!
    T = Kc * P       # T = Kc * P = D^(-1/2) * U * Delta
    Kpca(X, Kt, T, P, sv, eig, D, DKt, vtot, weights, kern, kwargs)
end

"""
    summary(object::Kpca, X)
Summarize the maximal (i.e. with maximal nb. PCs) fitted model.
* `object` : The fitted model.
* `X` : The X-data that was used to fit the model.
"""
function Base.summary(object::Kpca)
    nlv = size(object.T, 2)
    TT = object.D * (object.T).^2
    tt = vec(sum(TT, dims = 1))
    sstot = sum(object.eig)
    pvar = tt / sstot
    cumpvar = cumsum(pvar)
    explvar = DataFrame(pc = 1:nlv, var = tt, pvar = pvar, cumpvar = cumpvar)
    (explvar = explvar,)
end

"""
    transform(object::Kpca, X; nlv = nothing)
Compute PCs (scores matrix "T") from a fitted model and X-data.
* `object` : The maximal fitted model.
* `X` : X-data for which PCs are computed.
* `nlv` : Nb. PCs to consider. If nothing, it is the maximum nb. PCs.
"""
function transform(object::Kpca, X; nlv = nothing)
    a = size(object.T, 2)
    isnothing(nlv) ? nlv = a : nlv = min(nlv, a)
    fkern = eval(Meta.parse(object.kern))
    K = fkern(X, object.X; object.dots...)
    DKt = object.D * K'
    vtot = sum(DKt, dims = 1)
    Kc = K .- vtot' .- object.vtot .+ sum(object.D * object.DKt')
    T = Kc * @view(object.P[:, 1:nlv])
    T
end


text:
"""
    Gibbs(n_iters, alg_1, alg_2)

Compositional MCMC interface.

Example:
```julia
alg = Gibbs(1000, HMC(1, 0.2, 3, :v1), PG(20, 1, :v2))
```
"""
mutable struct Gibbs{A} <: InferenceAlgorithm
    n_iters   ::  Int     # number of Gibbs iterations
    algs      ::  A   # component sampling algorithms
    thin      ::  Bool    # if thinning to output only after a whole Gibbs sweep
    gid       ::  Int
end
Gibbs(n_iters::Int, algs...; thin=true) = Gibbs(n_iters, algs, thin, 0)
Gibbs(alg::Gibbs, new_gid) = Gibbs(alg.n_iters, alg.algs, alg.thin, new_gid)

const GibbsComponent = Union{Hamiltonian,MH,PG}

function Sampler(alg::Gibbs)
    n_samplers = length(alg.algs)
    samplers = Array{Sampler}(undef, n_samplers)

    space = Set{Symbol}()

    for i in 1:n_samplers
        sub_alg = alg.algs[i]
        if isa(sub_alg, GibbsComponent)
            samplers[i] = Sampler(typeof(sub_alg)(sub_alg, i))
        else
            @error("[Gibbs] unsupport base sampling algorithm $alg")
        end
        space = union(space, sub_alg.space)
    end

    # Sanity check for space
    @assert issubset(Turing._compiler_[:pvars], space) "[Gibbs] symbols specified to samplers ($space) doesn't cover the model parameters ($(Turing._compiler_[:pvars]))"

    if Turing._compiler_[:pvars] != space
        @warn("[Gibbs] extra parameters specified by samplers don't exist in model: $(setdiff(space, Turing._compiler_[:pvars]))")
    end

    info = Dict{Symbol, Any}()
    info[:samplers] = samplers

    Sampler(alg, info)
end

function sample(
                model::Function,
                alg::Gibbs;
                save_state=false,         # flag for state saving
                resume_from=nothing,      # chain to continue
                reuse_spl_n=0             # flag for spl re-using
                )

    # Init the (master) Gibbs sampler
    spl = reuse_spl_n > 0 ? resume_from.info[:spl] : Sampler(alg)

    @assert typeof(spl.alg) == typeof(alg) "[Turing] alg type mismatch; please use resume() to re-use spl"

    # Initialize samples
    sub_sample_n = []
    for sub_alg in alg.algs
        if isa(sub_alg, GibbsComponent)
            push!(sub_sample_n, sub_alg.n_iters)
        else
            @error("[Gibbs] unsupport base sampling algorithm $alg")
        end
    end

    # Compute the number of samples to store
    n = reuse_spl_n > 0 ? reuse_spl_n : alg.n_iters
    sample_n = n * (alg.thin ? 1 : sum(sub_sample_n))

    # Init samples
    time_total = zero(Float64)
    samples = Array{Sample}(undef, sample_n)
    weight = 1 / sample_n
    for i = 1:sample_n
        samples[i] = Sample(weight, Dict{Symbol, Any}())
    end

    # Init parameters
    varInfo = if resume_from == nothing
        vi_ = VarInfo()
        Base.invokelatest(model, vi_, HamiltonianRobustInit())
        vi_
    else
        resume_from.info[:vi]
    end

    n = spl.alg.n_iters; i_thin = 1

    # Gibbs steps
    PROGRESS[] && (spl.info[:progress] = ProgressMeter.Progress(n, 1, "[Gibbs] Sampling...", 0))
    for i = 1:n
        @debug "Gibbs stepping..."

        time_elapsed = zero(Float64)
        lp = nothing; epsilon = nothing; lf_num = nothing

        for local_spl in spl.info[:samplers]
            last_spl = local_spl
      # PROGRESS[] && haskey(spl.info, :progress) && (local_spl.info[:progress] = spl.info[:progress])

            @debug "$(typeof(local_spl)) stepping..."

            if isa(local_spl.alg, GibbsComponent)
                if isa(local_spl.alg, Hamiltonian)  # clean cache
                    local_spl.info[:grad_cache] = Dict{UInt64,Vector}()
                    local_spl.info[:reverse_diff_cache] = Dict()
                end

                for _ = 1:local_spl.alg.n_iters
                    @debug "recording old θ..."
                    time_elapsed_thin = @elapsed varInfo = step(model, local_spl, varInfo, i==1)

                    if ~spl.alg.thin
                        samples[i_thin].value = Sample(varInfo).value
                        samples[i_thin].value[:elapsed] = time_elapsed_thin
                        if ~isa(local_spl.alg, Hamiltonian)
                            # If statement below is true if there is a HMC component which provides lp and epsilon
                            if lp != nothing samples[i_thin].value[:lp] = lp end
                            if epsilon != nothing samples[i_thin].value[:epsilon] = epsilon end
                            if lf_num != nothing samples[i_thin].value[:lf_num] = lf_num end
                        end
                        i_thin += 1
                    end
                    time_elapsed += time_elapsed_thin
                end

                if isa(local_spl.alg, Hamiltonian)
                    lp = getlogp(varInfo)
                    epsilon = local_spl.info[:wum][:ϵ][end]
                    lf_num = local_spl.info[:lf_num]
                end
            else
                @error("[Gibbs] unsupport base sampler $local_spl")
            end
        end

        time_total += time_elapsed

        if spl.alg.thin
            samples[i].value = Sample(varInfo).value
            samples[i].value[:elapsed] = time_elapsed
            # If statement below is true if there is a HMC component which provides lp and epsilon
            if lp != nothing samples[i].value[:lp] = lp end
            if epsilon != nothing samples[i].value[:epsilon] = epsilon end
            if lf_num != nothing samples[i].value[:lf_num] = lf_num end
        end

        if PROGRESS[]
            if haskey(spl.info, :progress)
                ProgressMeter.update!(spl.info[:progress], spl.info[:progress].counter + 1)
            end
        end
    end

    @info("[Gibbs] Finished with")
    @info("  Running time    = $time_total;")

    if resume_from != nothing   # concat samples
        pushfirst!(samples, resume_from.value2...)
    end
    c = Chain(0, samples)       # wrap the result by Chain

    if save_state               # save state
        save!(c, spl, model, varInfo)
    end

    return c
end


text:
include("death_rate_zero.jl")
include("birth_rate_zero.jl")
include("equal_rates.jl")
include("unequal_rates.jl")

"""
    gradient_hessian(η, i, j, t)

Compute the gradient and Hessian of the log-transition probability of a simple
birth and death process evaluated at the point ``η = (λ, μ)^{\\prime}``.
Variable `i` is the initial population size, `j` is the final population size,
and `t` is the elapsed time.
"""
function gradient_hessian(
  η::Vector{F},
  i::I,
  j::I,
  t::F
)::Tuple{Vector{F}, Symmetric{F, Matrix{F}}} where {
  F <: AbstractFloat,
  I <: Integer
}
  ϵ = floatmin(F)

  if t < ϵ
    fill(zero(F), 2), Symmetric(fill(zero(F), (2, 2)))
  elseif (η[1] < ϵ) && (η[2] < ϵ)
    if i == j
      fill(- i * t, 2), Symmetric([zero(F) (i * t)^2; (i * t)^2 zero(F)])
    else
      fill(F(NaN), 2), Symmetric(fill(F(NaN), (2, 2)))
    end
  elseif η[1] < ϵ
    gradient_hessian_λ_zero(η[2], i, j, t)
  elseif η[2] < ϵ
    gradient_hessian_μ_zero(η[1], i, j, t)
  elseif η[2] ≈ η[1]
    # TODO: Improve numerical accuracy of formulas
    if (F === BigFloat) || (η[1] * t > 0.001)
      gradient_hessian_equal_rates(η[1], F(i), F(j), t)
    else
      setprecision(BigFloat, 256) do
        l = BigFloat(η[1])
        a = BigFloat(i)
        b = BigFloat(j)
        s = BigFloat(t)

        g, H = gradient_hessian_equal_rates(l, a, b, s)

        F[g[1], g[2]], Symmetric(F[H[1, 1] H[2, 1]; H[1, 2] H[2, 2]])
      end
    end
  else
    if (F === BigFloat) || (log(η[1] / η[2]) / ((η[1] - η[2]) * t) < 1_000)
      gradient_hessian_unequal_rates(η, F(i), F(j), t)
    else
      setprecision(BigFloat, 256) do
        l = BigFloat(η[1])
        m = BigFloat(η[2])
        a = BigFloat(i)
        b = BigFloat(j)
        s = BigFloat(t)

        g, H = gradient_hessian_unequal_rates([l, m], a, b, s)

        F[g[1], g[2]], Symmetric(F[H[1, 1] H[2, 1]; H[1, 2] H[2, 2]])
      end
    end
  end
end

function gradient_hessian(
  η::Vector{F},
  x::ObservationDiscreteTimeEven
)::Tuple{Vector{F}, Symmetric{F, Matrix{F}}} where {
  F <: AbstractFloat
}
  ∇ = fill(zero(Float64), 2)
  H = fill(zero(Float64), (2, 2))

  for n = 1:size(x.state, 2), s = 2:size(x.state, 1)
    a, b = gradient_hessian(η, x.state[s - 1, n], x.state[s, n], x.u)
    ∇ .+= a
    H .+= b
  end

  ∇, Symmetric(H)
end


text:
function eval_balmer_pseudocont(Temp, Ne, fwhm)
    function readfile(file, Ne)
        lower = Vector{Float64}()
        upper = Vector{Float64}()
        norm  = Vector{Float64}()
        density = NaN
        current_upper = 0
        for l in readlines(file)
            l = strip(l)
            while occursin(" =", l); l = replace(l, " =" => "="); end
            while occursin("= ", l); l = replace(l, "= " => "="); end
            s = split(l)
            if s[1] == "DENS"
                density = NaN
            elseif !isfinite(density)
                density = parse(Float64, s[1])
            elseif density == Ne
                if (length(s[1]) >= 4)  &&  (s[1][1:4] == "E_NU")
                    current_upper = parse(Int, split(s[1],'=')[2])
                elseif occursin(r"\D+", s[1])
                    current_upper = 0
                elseif (occursin(r"\d+", s[1]))  &&  (current_upper > 0)
                    @assert mod(length(s), 2) == 0
                    for ii in 1:2:length(s)
                        push!(upper, current_upper)
                        push!(lower, parse(Int    , s[ii]))
                        push!(norm , parse(Float64, s[ii+1]))
                    end
                end
            end
        end
        norm ./= maximum(norm)
        return (lower, upper, norm)
    end

    # List of available temperatures in SH95 [K/100]
    Tavail = [5, 10, 30, 50, 75, 100, 125, 150, 200, 300]
    iTemp = argmin(abs.(Temp./100. .- Tavail))
    if Temp != Tavail[iTemp]*100
        @info "Balmer pseudo-continuum: T (requested)=$Temp, T (available) = $(Tavail[iTemp]*100)"
    end

    # List of available electron densities in SH95 [log (Ne / cm^-3)]
    logNeAvail = [2., 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14]
    iNe = argmin(abs.(log10(Ne) .- logNeAvail))
    if Ne != 10^logNeAvail[iNe]
        @info "Balmer pseudo-continuum: Ne (requested)=$Ne, Ne (available) = $(10^logNeAvail[iNe])"
    end

    # Read appropriate file
    path = qsfit_data() * "/SH1995/"
    file = path * @sprintf("r1b%04d.d", Tavail[iTemp])
    electronDensity = 10^logNeAvail[iNe]

    (lower, upper, norm) = readfile(file, electronDensity)

    ryd = 911.2671808036101 # Rydberg in A
    edge = 4 * ryd          # Balmer (2^2=4) edge in Angstrom

    wave = ryd ./ (1 ./ lower.^2 - 1 ./ upper.^2)

    ii = findall((lower .== 2)  .&  (upper .> 6))
    wave = wave[ii]
    norm = norm[ii]

    # Prepare output variables
    λ = 10 .^range(log10(3400.), stop=log10(4300), length=500) # Angstrom
    cont = fill(0., length(λ))       # Lum. density (arbitrary units)

    # Broadening of high order Balmer lines ==> pseudo-continuum
    σ = (fwhm / 3.e5) .* wave ./ 2.355
    for i in 1:length(wave)
        cont += norm[i] * gauss.(λ, wave[i], σ[i])
    end
    cont ./= Spline1D(λ, cont, k=1, bc="error")(edge)
    return (λ, cont)
end

function eval_balmer_continuum(Temp, Tau, fwhm)
    ryd = 911.2671808036101 # Rydberg in A
    edge = 4 * ryd          # Balmer (2^2=4) edge in Angstrom

    λ = 10 .^range(log10(912.), stop=log10(5000), length=600) .* 1e-8 # cm

    # Planck function
    b = planck(λ, Temp)

    # Take into account optical depth
    λ .*= 1e8
    cont = b .* (1. .- exp.(-(Tau .* (λ ./ edge).^3.)))
    cont[findall(λ .>= edge)] .= 0.
    cont ./= maximum(cont)

    # Broadening
    σ = fwhm / 3.e5 / 2.355 * edge
    kernel = gauss(λ, mean(λ), σ)
    kernel = kernel[findall(kernel .> maximum(kernel) / 1e3)]
    (mod(length(kernel), 2) == 0)  &&  (push!(kernel, 0.))
    cont = convol(cont, kernel)

    # Normalize Balmer continuum at 3000A
    cont ./=     Spline1D(λ, cont, k=1, bc="extrapolate")(3000.)
    contAtEdge = Spline1D(λ, cont, k=1, bc="extrapolate")(3000.)
    return (λ, cont, contAtEdge)
end

# ____________________________________________________________________
# balmercont
#
mutable struct balmercont <: AbstractComponent
    norm::Parameter
    ratio::Parameter
    c1::Vector{Float64}
    c2::Vector{Float64}

    function balmercont(norm::Number, ratio::Number)
        out = new(Parameter(norm), Parameter(ratio), Vector{Float64}(), Vector{Float64}())
        out.norm.low = 0
        out.ratio.low = 0
        return out
    end
end

function prepare!(comp::balmercont, domain::Domain{1})
    T = 15000.
    Ne = 1e9
    Tau = 1.
    fwhm = 5000.
    (λ1, c1, contAtEdge) = eval_balmer_continuum(T, Tau, fwhm)
    (λ2, c2) = eval_balmer_pseudocont(T, Ne, fwhm)
    c2 .*= contAtEdge
    comp.c1 = Spline1D(λ1, c1, k=1, bc="extrapolate")(domain[:])
    comp.c2 = Spline1D(λ2, c2, k=1, bc="extrapolate")(domain[:])
    comp.c1[findall(comp.c1 .< 0.)] .= 0.
    comp.c2[findall(comp.c2 .< 0.)] .= 0.
    return fill(NaN, length(domain))
end

function evaluate!(buffer, comp::balmercont, domain::Domain{1},
                   norm, ratio)
    buffer .= norm .* (comp.c1 .+ ratio .* comp.c2)
end

#=
x = Domain(1000.:4500)
model = Model(x, :c => QSFit.balmercont(0.1, 0.5))
@gp x[:] model() "w l"

@gp "set grid" :-
@gp :- xlab="Wavelength [A]" ylab="Lum. density [arb.units]" :-
@gp :- 3645.07.*[1,1] [0, 1.1] "w l dt 4 lc rgb 'black'" :-
model[:c].norm.val = 0.5
model[:c].ratio.val = 1
evaluate!(model)
@gp :- x[:] model() "w l lw 2 lc rgb 'red'"
model[:c].ratio.val = 0.5
evaluate!(model)
@gp :- x[:] model() "w l lw 2 lc rgb 'blue'"
=#


text:
"""
this function solves Dantzig-Wolfe ralaxed main model.
"""
function solve_dw_main_model(cgd::CGdata, optimizer::MOI.AbstractOptimizer)
    λ_counter       = cgd.λ_counter
    μ_counter       = cgd.μ_counter
    (λ_counter > 0) && (extr_ptn_sol = cgd.extr_ptn_sol[λ_counter].solution)
    (μ_counter > 0) && (extr_dir_sol = cgd.extr_dir_sol[μ_counter].solution)
    disc2var_idx    = cgd.pricing_sub_disc2var_idx
    cons1           = cgd.dw_main_constr_ref[1]
    cons2           = cgd.dw_main_constr_ref[2]
    cons3           = cgd.dw_main_constr_ref[3]
    nitr            = cgd.num_oa_iter
    nbin            = length(cgd.jump_bin_var)
    constr_varidx   = cgd.dw_main_constr_varidx
    constr_coef     = cgd.dw_main_constr_coef
    λ               = cgd.dw_main_x[1]
    μ               = cgd.dw_main_x[2]
    (λ_counter > 0) && print("\n", extr_ptn_sol, "\n")
    (μ_counter > 0) && print("\n", extr_dir_sol, "\n")
    #update coefficient of constr set "z[1] >= sum(u_p[1,i]*L[1]*λ[i] for i in Κ) + sum(u_d[1,j]*L[1]*μ[j] for j in Τ)"
    if cgd.psm_feasible == true
        for i in 1:nbin
            #val = cgd.obj_terms_lower_bound[i]*extr_ptn_sol[disc2var_idx[i]]
            val = -1*cgd.obj_terms_lower_bound[i]*extr_ptn_sol[i]
            JuMP.set_normalized_coefficient(cons1[i], λ[λ_counter], val)
        end
    elseif cgd.psm_extr_dir_feasible == true
        for i in 1:nbin
            #val = cgd.obj_terms_lower_bound[i]*extr_dir_sol[disc2var_idx[i]]
            val = -1*cgd.obj_terms_lower_bound[i]*extr_dir_sol[i]
            JuMP.set_normalized_coefficient(cons1[i], μ[μ_counter], val)
        end
    end
    #update coefficient of constr set "z[1]+∇f1(̄x)*̄x+U[1] >= sum((u_p[1,i]*U[1]+∇f1(̄x)*v_p[i])*λ[i] for i in Κ)+sum((u_d[1,j]*U[1]+∇f1(̄x)*v_d[j])*μ[j] for j in Τ)"
    number = 0
    if cgd.psm_feasible == true
        for i in 1:nitr
            for j in 1:nbin
                number += 1
                val = dot(constr_coef[i][j],extr_ptn_sol[constr_varidx[i][j]])
                JuMP.set_normalized_coefficient(cons2[number], λ[λ_counter], val)
            end
        end
    elseif cgd.psm_extr_dir_feasible == true
        for i in 1:nitr
            for j in 1:nbin
                number += 1
                val = dot(constr_coef[i][j],extr_dir_sol[constr_varidx[i][j]])
                JuMP.set_normalized_coefficient(cons2[idx], μ[μ_counter], val)
            end
        end
    end
    #update for the third set of constraints
    if cgd.psm_feasible == true
        JuMP.set_normalized_coefficient(cons3[1], λ[λ_counter], 1.0)
    end
    #solve the dw_main_model
    JuMP.set_optimizer(cgd.dw_main_model, optimizer.options.lp_solver)
    JuMP.optimize!(cgd.dw_main_model)
    if (JuMP.termination_status(cgd.dw_main_model) ==  MOI.OPTIMAL) || (JuMP.termination_status(cgd.dw_main_model) ==  MOI.LOCALLY_SOLVED)
        obj_val              = JuMP.objective_value(cgd.dw_main_model)
        dw_main_x            = JuMP.all_variables(cgd.dw_main_model)
        cgd.dw_main_solution = JuMP.value.(dw_main_x)
        dw_sol_obj           = dwmSolutionObj(cgd.dw_main_solution, obj_val)
        push!(cgd.dw_main_sol, dw_sol_obj)
        cgd.λ_val = JuMP.value.(cgd.dw_main_x[1])
        cgd.μ_val = JuMP.value.(cgd.dw_main_x[2])
    else
        @warn "The dw_main_model is infeasible.\n"
    end
    #get dual solutions
    dual_sol_exists = JuMP.has_duals(cgd.dw_main_model)
    dual_solution = []
    if dual_sol_exists == true
        w1 = JuMP.dual.(cons1)
        w2 = JuMP.dual.(cons2)
        w3 = JuMP.dual.(cons3)
        w = [w1, w2, w3]
        for d in w
            for d_val in d
                push!(dual_solution, d_val)
            end
        end
    end
    cgd.dw_main_dual_solution = dual_solution
    #return dual solution
    return nothing
end

"""
This function solves the dw pricing sub model.
"""
function solve_pricing_sub_model(cgd::CGdata, optimizer::MOI.AbstractOptimizer)
    dwmd_sol     = cgd.dw_main_dual_solution
    model        = cgd.pricing_sub_model
    x            = JuMP.all_variables(model)
    obj_terms_lb = cgd.obj_terms_lower_bound
    varid        = cgd.dw_main_constr_varidx
    coef         = cgd.dw_main_constr_coef
    bin_vars     = cgd.jump_bin_var
    #num_constr   = cgd.dw_main_num_constr
    nitr         = cgd.num_oa_iter
    nbin = length(bin_vars)
    num_constr   = nbin+nbin*nitr+1
    #set obj func
    (cgd.pricing_sub_obj_sense == :MAX_SENSE) && JuMP.@objective(model, Max, sum(bin_vars[i]*obj_terms_lb[i]*dwmd_sol[i] for i in 1:nbin)-sum(sum(dot(coef[i][j],x[varid[i][j]])*dwmd_sol[nbin+(i-1)*nbin+j] for j in 1:nbin) for i in 1:nitr)-dwmd_sol[num_constr])
    (cgd.pricing_sub_obj_sense == :MIN_SENSE) && JuMP.@objective(model, Min, sum(bin_vars[i]*obj_terms_lb[i]*dwmd_sol[i] for i in 1:nbin)-sum(sum(dot(coef[i][j],x[varid[i][j]])*dwmd_sol[nbin+(i-1)*nbin+j] for j in 1:nbin) for i in 1:nitr)-dwmd_sol[num_constr])
    #set lp solver
    JuMP.set_optimizer(model, optimizer.options.lp_solver)
    JuMP.optimize!(model)
    if (JuMP.termination_status(model) ==  MOI.OPTIMAL) || (JuMP.termination_status(model) ==  MOI.LOCALLY_SOLVED)
        cgd.pricing_sub_solution   = JuMP.value.(x);
        objective_value            = JuMP.objective_value(model);
        psm_sol_obj                = psmSolutionObj(cgd.pricing_sub_solution, objective_value)
        cgd.psm_feasible           = true
        cgd.psm_extr_dir_feasible  = false
        cgd.λ_counter             += 1
        push!(cgd.extr_ptn_sol, psm_sol_obj)
    else
        @info "\n The pricing sub problem is infeasible. \n"
    end
end

"""
This function solves the pricing sub problem for extreme directions
"""
function solve_pricing_sub_extm_dir_model(cgd::CGdata, optimizer::MOI.AbstractOptimizer)
    dwmd_sol     = cgd.dw_main_dual_solution
    obj_terms_lb = cgd.obj_terms_lower_bound
    varid        = cgd.dw_main_constr_varidx
    coef         = cgd.dw_main_constr_coef
    num_constr   = cgd.dw_main_num_constr
    siz          = length(cgd.pricing_sub_disc2var_idx)
    model        = cgd.pricing_sub_extm_dir_model
    x            = cgd.pricing_sub_extm_dir_x[1]
    nitr         = cgd.num_oa_iter
    #define obj function
    (cgd.pricing_sub_obj_sense == :MAX_SENSE) && JuMP.@objective(model, Max, sum(x[cgd.pricing_sub_disc2var_idx[i]]*obj_terms_lb[i]*dwmd_sol[i] for i in 1:siz)-sum(sum(dot(coef[i][j],x[varid[i][j]])*dwmd_sol[siz+(i-1)*nbin+j] for j in 1:nbin) for i in 1:nitr))
    (cgd.pricing_sub_obj_sense == :MIN_SENSE) && JuMP.@objective(model, Min, sum(x[cgd.pricing_sub_disc2var_idx[i]]*obj_terms_lb[i]*dwmd_sol[i] for i in 1:siz)-sum(sum(dot(coef[i][j],x[varid[i][j]])*dwmd_sol[siz+(i-1)*nbin+j] for j in 1:nbin) for i in 1:nitr))
    #set lp solver
    JuMP.set_optimizer(model, optimizer.options.lp_solver)
    JuMP.optimize!(model)
    if (JuMP.termination_status(model) ==  MOI.OPTIMAL) || (JuMP.termination_status(model) ==  MOI.LOCALLY_SOLVED)
        cgd.pricing_sub_extr_dir_solution   = JuMP.value.(x)
        objective_value                     = JuMP.objective_value(model)
        psm_extm_dir_sol_obj                = psmExtrmDirSolutionObj(pricing_sub_extr_dir_solution, objective_value)
        cgd.psm_feasible                    = false
        cgd.psm_extr_dir_feasible           = true
        cgd.μ_counter                      += 1
        push!(cgd.extr_dir_sol, psm_extm_dir_sol_obj)
    else
        @info "\n The pricing sub problem for extreme direction is infeasible. \n"
    end
end

"""
The column generation algorithm is stated here.
"""
function cg_algorithm(model::JuMP.Model, oad::OAdata, op::HybridProblem, optimizer::MOI.AbstractOptimizer)
    #initiate cg data storage
    cgd = CGdata()
    #populate cgd with oad data
    init_cg_data(cgd, oad)
    #decompose the model into main and sub problems
    decompose_lp_model(model, optimizer, op, oad, cgd)
    duration = 100
    start = time()
    while((time()-start) < duration)
        cgd.num_cg_iter += 1
        #solve dw main problem
        if cgd.num_cg_iter <= 1
            #solve pricing sub problem for a feasible solution
            ps_model = cgd.pricing_sub_model
            x        = JuMP.all_variables(ps_model)
            #set obj func
            (cgd.pricing_sub_obj_sense == :MAX_SENSE) && JuMP.@objective(ps_model, Max, 0)
            (cgd.pricing_sub_obj_sense == :MIN_SENSE) && JuMP.@objective(ps_model, Min, 0)
            #set lp solver
            JuMP.set_optimizer(ps_model, optimizer.options.lp_solver)
            print(ps_model)
            JuMP.optimize!(ps_model)
            if (JuMP.termination_status(ps_model) ==  MOI.OPTIMAL) || (JuMP.termination_status(ps_model) ==  MOI.LOCALLY_SOLVED)
                cgd.pricing_sub_solution   = JuMP.value.(x)
                objective_value            = JuMP.objective_value(ps_model)
                psm_sol_obj                = psmSolutionObj(cgd.pricing_sub_solution, objective_value)
                cgd.psm_feasible           = true
                cgd.psm_extr_dir_feasible  = false
                cgd.λ_counter             += 1
                push!(cgd.extr_ptn_sol, psm_sol_obj)
            else
                @info "\n The pricing sub problem is infeasible. \n"
                #solve pricing sub problem for feasible extreme directions
                ps_dir_model   = cgd.pricing_sub_extm_dir_model
                x              = cgd.pricing_sub_extm_dir_x[1]
                #define obj function
                (cgd.pricing_sub_obj_sense == :MAX_SENSE) && JuMP.@objective(ps_dir_model, Max, 0)
                (cgd.pricing_sub_obj_sense == :MIN_SENSE) && JuMP.@objective(ps_dir_model, Min, 0)
                #set lp solver
                JuMP.set_optimizer(ps_dir_model, optimizer.options.lp_solver)
                JuMP.optimize!(ps_dir_model)
                if (JuMP.termination_status(ps_dir_model) ==  MOI.OPTIMAL) || (JuMP.termination_status(ps_dir_model) ==  MOI.LOCALLY_SOLVED)
                    cgd.pricing_sub_extr_dir_solution   = JuMP.value.(x)
                    objective_value                     = JuMP.objective_value(ps_dir_model)
                    psm_extm_dir_sol_obj                = psmExtrmDirSolutionObj(pricing_sub_extr_dir_solution, objective_value)
                    cgd.psm_feasible                    = false
                    cgd.psm_extr_dir_feasible           = true
                    cgd.μ_counter                      += 1
                    push!(cgd.extr_dir_sol, psm_extm_dir_sol_obj)
                else
                    @info "\n The pricing sub problem for extreme direction is infeasible. \n"
                end
            end
            #solve dw_main_model
            solve_dw_main_model(cgd, optimizer)
        elseif (cgd.num_cg_iter >1) && (cgd.num_cg_iter <= 5)
            solve_dw_main_model(cgd, optimizer)
        #elseif TO DO: for n>50 need to update basis of the dw_model instead of adding columns
        elseif cgd.num_cg_iter > 99
            @info "maximum number of column generation iteration reached. \n"
            break
        end
        #solve dw pricing sub problem
        solve_pricing_sub_model(cgd, optimizer)
        #check cg optimality
        if cgd.psm_feasible && (cgd.num_cg_iter > 3)
            if (cgd.extr_ptn_sol[cgd.λ_counter].obj_value >= 0)
                @info "The optimality condition has been satisfied. \n"
                break
            end
        end
        #solve pricing sub for extreme direction in case psm infeasible
        !cgd.psm_feasible && solve_pricing_sub_extm_dir_model(cgd, optimizer)
    end
    ExtremePoints = cgd.extr_ptn_sol
    ExtremeDirections = cgd.extr_dir_sol
    dw_solution = cgd.dw_main_sol[cgd.num_cg_iter-1].solution
    (cgd.λ_counter > 0) && (lp_ext_sol = sum(ExtremePoints[i].solution*cgd.λ_val[i] for i in 1:cgd.λ_counter))
    (cgd.μ_counter > 0) && (lp_ext_dir_sol = sum(ExtremeDirections[i].solution*cgd.μ_val[i] for i in 1:cgd.μ_counter))
    (cgd.λ_counter > 0) && (cgd.μ_counter > 0) && (lp_solution = lp_ext_sol+lp_ext_dir_sol)
    (cgd.λ_counter > 0) && !(cgd.μ_counter > 0) && (lp_solution = lp_ext_sol)
    cgd.best_solution = lp_solution
    #return solution of the input model
    print("The optimal lp solution using column generation method is :\n", lp_solution, "\n Thank you for having patience.\n")
    return nothing
end

max_stars_repo_name: MdMahbubarRahman/Duran

max_issues_repo_name: MdMahbubarRahman/Duran

max_forks_repo_name: MdMahbubarRahman/Duran


text:
module TestDefaults

using Mimi
using Test

@defcomp A begin
    p1 = Parameter(default = 1)
    p2 = Parameter()
end

m = Model()
set_dimension!(m, :time, 1:10)
add_comp!(m, A)
set_param!(m, :p2, 2)

# So far only :p2 is in the model definition's dictionary
@test length(m.md.external_params) == 1

run(m)

# During build, :p1's value is set to it's default
@test m[:A, :p1] == 1

# But the original model definition does not have :p1 in it's external parameters
@test length(m.md.external_params) == 1
@test length(m.mi.md.external_params) == 2      # But the model instance's md is where the default value was set
@test ! (:p1 in keys(m.md.external_params))
@test :p1 in keys(m.mi.md.external_params)

# This errors because p1 isn't in the model definition's external params
@test_throws ErrorException update_param!(m, :p1, 10)

# Need to use set_param! instead
set_param!(m, :p1, 10)

# Now there is a :p1 in the model definition's dictionary
@test :p1 in keys(m.md.external_params)

run(m)
@test m[:A, :p1] == 10
update_param!(m, :p1, 11)   # Now we can use update_param!

end


text:
using Test

tests = [
  "utils",
  "sparseIterate",
  "proximal_functions",
  "differentiable_functions"
]

for t in tests
	f = "$t.jl"
	println("* running $f ...")
	t = @elapsed include(f)
    println("done (took $t seconds).")
end


text:
import ..Common.test_constructor

function run_test_constructor()
    test_constructor(
        Dict(
            :index => 1,
            :atom_name => "A",
            :residue_name => "RES1",
            :residue_id => 1,
            :segment_name => "SEG1",
            :chain_name => "C",
            :system_name => "S",
            :coordinate => [0., 1., 0.]
        ),

        Atom(;
            index=1,
            residue_name="RES1",
            residue_id=1,
            segment_name="SEG1",
            chain_name="C",
            system_name="S",
            atom_name="A",
            coordinate=[0., 1., 0.]
        )
    )
end


text:
function qnwbeta(n::Int, a::T, b::S) where {T <: Real, S <: Real}
    a -= 1
    b -= 1

    maxit = 25

    x = zeros(n)
    w = zeros(n)

    for i=1:n
        if i == 1
            an = a / n
            bn = b / n
            r1 = (1 + a) * (2.78 / (4 + n * n) + 0.768an / n)
            r2 = 1 + 1.48 * an + 0.96bn + 0.452an*an + 0.83an*bn
            z = 1 - r1 / r2

        elseif i == 2
            r1 = (4.1 + a) / ((1 + a) * (1 + 0.156a))
            r2 = 1 + 0.06 * (n - 8) * (1 + 0.12a) / n
            r3 = 1 + 0.012b * (1 + 0.25 * abs(a)) / n
            z = z - (1 - z) * r1 * r2 * r3

        elseif i == 3
            r1 = (1.67 + 0.28a) / (1 + 0.37a)
            r2 = 1 + 0.22 * (n - 8) / n
            r3 = 1 + 8 * b / ((6.28 + b) * n * n)
            z = z - (x[1] - z) * r1 * r2 * r3

        elseif i == n - 1
            r1 = (1 + 0.235b) / (0.766 + 0.119b)
            r2 = 1 / (1 + 0.639 * (n - 4) / (1 + 0.71 * (n - 4)))
            r3 = 1 / (1 + 20a / ((7.5+ a ) * n * n))
            z = z + (z - x[n-3]) * r1 * r2 * r3

        elseif i == n
            r1 = (1 + 0.37b) / (1.67 + 0.28b)
            r2 = 1 / (1 + 0.22 * (n - 8) / n)
            r3 = 1 / (1 + 8 * a / ((6.28+ a ) * n * n))
            z = z + (z - x[n-2]) * r1 * r2 * r3

        else
            z = 3 * x[i-1] - 3 * x[i-2] + x[i-3]
        end

        ab = a + b

        for its = 1:maxit
            temp = 2 + ab
            p1 = (a - b + temp * z) / 2
            p2 = 1
            for j=2:n
              p3 = p2
              p2 = p1
              temp = 2 * j + ab
              aa = 2 * j * (j + ab) * (temp - 2)
              bb = (temp - 1) * (a * a - b * b + temp * (temp - 2) * z)
              c = 2 * (j - 1 + a) * (j - 1 + b) * temp
              p1 = (bb * p2 - c * p3) / aa
            end
            pp = (n * (a - b - temp * z) * p1 +
                  2 * (n + a) * (n + b) * p2) / (temp * (1 - z * z))
            z1 = z
            z = z1 - p1 ./ pp
            if abs(z - z1) < 3e-14 break end
        end

        if its >= maxit
            error("Failure to converge in qnwbeta1")
        end

        x[i] = z
        w[i] = temp / (pp * p2)
    end

    x = (1 - x) ./ 2
    w = w * exp(gammaln(a + n) +
                gammaln(b + n) -
                gammaln(n + 1) -
                gammaln(n + ab + 1))
    w = w / (2 * exp(gammaln(a + 1) +
                     gammaln(b + 1) -
                     gammaln(ab + 2)))

    return x, w
end


text:
function start!(mcmc::MCMC{T},
                n_chains::Int64;
                attempts::Int64=1000) where T <: EpidemicModel
  pmeter = Progress(n_chains * attempts, "Initialization progress")
  pchannel = RemoteChannel(()->Channel{Bool}(10), 1)
  mc_Futures = Future[]
  for i = 1:n_chains
    @logmsg LogLevel(-3000) "Creating Markov chain $i initialization Future"
    push!(mc_Futures, @spawn initialize(MarkovChain, mcmc, pchannel, attempts=attempts))
  end
  @logmsg LogLevel(-2000) "$n_chains Markov chain initialization Futures successfully created"
  while !all(isready.(mc_Futures))
    take!(pchannel) && next!(pmeter)
  end
  close(pchannel)
  finish!(pmeter)
  markov_chains = [fetch(i) for i in mc_Futures]
  append!(mcmc.markov_chains, markov_chains[nothing.!==markov_chains])
  @logmsg LogLevel(-3000) "$(sum(nothing.!==markov_chains)) Markov chains successfully initialized"
  return mcmc
end

function start!(mcmc::MCMC{T};
                attempts::Int64=1000) where T <: EpidemicModel
  # return start!(mcmc, 1, attempts)
  mc = initialize(MarkovChain, mcmc, attempts=attempts)
  if mc !== nothing
    push!(mcmc.markov_chains, mc)
  else
    @error "Failed to initialize a Markov chain"
  end
  return mcmc
end


text:
using Test;

@

max_stars_repo_name: daviddoret/julia-cybersecurity

max_issues_repo_name: daviddoret/julia-cybersecurity

max_forks_repo_name: daviddoret/julia-cybersecurity


text:
# check the Bernoulli numbers and polynomials
#    this is mostly just a check that code is correct
#    most errors are arising because of errors in the functions being used
using Polylogarithms
using BenchmarkTools
using DataFrames, CSV
using PyPlot
include("utilities.jl")

##############################################
# Bernoulli numbers
B = Symbol("B_n")
data1 = CSV.read("../data/bernoulli_test_data_1.csv", DataFrame)
m = size(data1,1)
B1 = zeros(Rational{Int64}, m)
B1_ref = zeros(Rational{Int64}, m)
error1 = zeros(Float64, m)
rel_error1 = zeros(Float64, m)
for i=1:m
    B1[i] = bernoulli( data1[i,:n] )
    B1_ref[i] = parse(Rational{Int64},data1[i,B])
    error1[i] = Float64( B1[i] - B1_ref[i]  )
    rel_error1[i] = error1[i]/Float64( B1_ref[i] )
end
println("max abs. error in Bernoulli numbers = $(maximum( abs.(error1) ))")

##############################################
# Bernoulli polynomials: identities
m = 20
B0 = zeros(Float64, m)
B0_ref = zeros(Rational{Int64}, m)
error0 = zeros(Float64, m)
for i=1:m
    B0_ref[i] = bernoulli( i )
    B0[i] = bernoulli( i, 0.0 )
    error0[i] = Float64( B0[i] - B0_ref[i]  )
end
figure(0)
clf()
semilogy( 1:m, abs.(error0) )
xlabel("n, (x=0.0)")
xticks(1:m)
ylabel("absolute error")
# ylim( 1.0e-17, 1.0e-15)
title("Bernoulli polynomials")
savefig("plots/bernoulli_bench_0.pdf")

##############################################
# Bernoulli polynomials: data v x
B = Symbol("B_n(x)")
data2 = CSV.read("../data/bernoulli_test_data_2.csv", DataFrame; delim=",", type=String)
data2[!,:n] = parse.(Float64, data2[!,:n] )
data2[!,:x] = parse.(Float64, data2[!,:x] )
data2[!,B] = parse.(Float64, data2[!,B] )
m = size(data2,1)
B2 = zeros(Float64, m)
B2_ref = zeros(Float64, m)
error2 = zeros(Float64, m)
rel_error2 = zeros(Float64, m)
for i=1:m
    B2[i] = bernoulli( Int(data2[i,:n]),  data2[i,:x])
    B2_ref[i] = data2[i,B]
    error2[i] = Float64( B2[i] - B2_ref[i]  )
    rel_error2[i] = error2[i] / B2_ref[i]
end
figure(2)
clf()
semilogy( data2[!,:x], abs.(rel_error2) )
xlabel("x, (n=12)")
ylabel("absolute rel. error")
# ylim( 1.0e-17, 1.0e-15)
title("Bernoulli polynomials")
savefig("plots/bernoulli_bench_2.pdf")

##############################################u
# Bernoulli polynomials: data v n
B = Symbol("B_n(x)")
data3 = CSV.read("../data/bernoulli_test_data_3.csv", DataFrame)
m = size(data3,1)
B3 = zeros(Float64, m)
B3_ref = zeros(Float64, m)
error3 = zeros(Float64, m)
rel_error3 = zeros(Float64, m)
for i=1:m
    B3[i] = bernoulli( Int(data3[i,:n]),  data3[i,:x])
    B3_ref[i] = data3[i,B]
    error3[i] = Float64( B3[i] - B3_ref[i]  )
    rel_error3[i] = error3[i] / B3_ref[i]
end
figure(3)
clf()
semilogy( data3[!,:n], abs.(rel_error3) )
xlabel("n, (x=1.5)")
ylabel("absolute rel. error")
# ylim( 1.0e-17, 1.0e-15)
title("Bernoulli polynomials")
savefig("plots/bernoulli_bench_3.pdf")


text:
###############################################################################
#
# Update soil from water content or pressure
#
###############################################################################
"""
    soil_moisture_swc!(node::SPACSimple{FT}, swc::FT) where {FT<:AbstractFloat}

Update soil moisture and soil matrix potential, given
- `node` [`SPACSimple`](@ref) type struct
- `swc` Given soil water content
"""
function soil_moisture_swc!(
            node::SPACSimple{FT},
            swc::FT
) where {FT<:AbstractFloat}
    if swc < node.mswc
        node.swc = swc;
    else
        node.swc = node.mswc;
    end

    node.p_soil = soil_p_25_swc(node.hs.root.sh, node.swc) * node.hs.root.f_st;

    return nothing
end

"""
    soil_moisture_p!(node::SPACSimple{FT}, p::FT) where {FT<:AbstractFloat}

Update soil moisture and soil matrix potential, given
- `node` [`SPACSimple`](@ref) type struct
- `p` Given soil maxtrix potential
"""
function soil_moisture_p!(
            node::SPACSimple{FT},
            p::FT
) where {FT<:AbstractFloat}
    node.p_soil = p;
    p_25        = p / node.hs.root.f_st;
    node.swc    = soil_swc(node.hs.root.sh, p_25);

    return nothing
end

"""
    soil_moisture_p25!(
                node::SPACSimple{FT},
                p_25::FT
    ) where {FT<:AbstractFloat}

Update soil moisture and soil matrix potential, given
- `node` [`SPACSimple`](@ref) type struct
- `p_25` Given soil maxtrix potential at 25 Celcius
"""
function soil_moisture_p25!(
            node::SPACSimple{FT},
            p_25::FT
) where {FT<:AbstractFloat}
    node.p_soil = p_25 * node.hs.root.f_st;
    node.swc    = soil_swc(node.hs.root.sh, p_25);

    return nothing
end

"""
    soil_moisture!(
                node::SPACSimple{FT},
                flow::FT,
                Δt::FT = FT(1)
    ) where {FT<:AbstractFloat}

Update soil moisture and soil matrix potential, given
- `node` [`SPACSimple`](@ref) type struct
- `flow` Mean outlet flow rate in `[Kg h⁻¹]`
- `Δt` Time period in `[h]`
"""
function soil_moisture!(
            node::SPACSimple{FT},
            flow::FT,
            Δt::FT = FT(1)
) where {FT<:AbstractFloat}
    # 1. positive flow means out, flow in Kg h⁻¹, dt in h
    m_all  = node.swc * node.gaba * node.h_soil * ρ_H₂O(FT);
    m_out  = flow * Δt;
    m_all -= m_out;
    swc    = m_all / (node.gaba * node.h_soil * ρ_H₂O(FT));
    swc    = max(swc, FT(1e-6));

    soil_moisture_swc!(node, swc);

    return nothing
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
# AUTO GENERATED FILE - DO NOT EDIT

export ''_fefferycaptcha

"""
    ''_fefferycaptcha(;kwargs...)

A FefferyCaptcha component.

Keyword arguments:
- `id` (optional)
- `bgColor` (optional)
- `captcha` (optional)
- `charNum` (optional)
- `className` (optional)
- `fontSize` (optional)
- `height` (optional)
- `loading_state` (optional)
- `setProps` (optional): Dash-assigned callback that should be called to report property changes
to Dash, to make them available for callbacks.
- `style` (optional)
- `width` (optional)
"""
function ''_fefferycaptcha(; kwargs...)
        available_props = Symbol[:id, :bgColor, :captcha, :charNum, :className, :fontSize, :height, :loading_state, :style, :width]
        wild_props = Symbol[]
        return Component("''_fefferycaptcha", "FefferyCaptcha", "feffery_utils_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: CNFeffery/feffery-utils-components

max_issues_repo_name: CNFeffery/feffery-utils-components

max_forks_repo_name: CNFeffery/feffery-utils-components


text:
export nanfill!

"""
  nanfill!(S::SeisData)
  nanfill!(C::SeisChannel)

Replace NaNs in `:x` with mean of non-NaN values.
"""
function nanfill!(S::GphysData)
  for i = 1:S.n
    if !isempty(S.x[i])
      nanfill!(S.x[i])
      note!(S, i, "nanfill!")
    end
  end
  return nothing
end
nanfill!(C::GphysChannel) = (nanfill!(C.x); note!(C, "nanfill!"))


text:
import sbp
err = sbp.MMS
using Test
using LinearAlgebra
using DataStructures

@testset "extract_vars" begin
        n = [4, 8, 12]
        data = Array(1:sum(n))
        vars = OrderedDict("u" => n[1], "v" => n[2], "w" => n[3])
        ans_u = Array(1:n[1])
        ans_v = Array(n[1]+1:n[1] + n[2])
        ans_w = Array(n[1]+n[2]+1:sum(n))
        expected = Dict("u" => ans_u,
                         "v" => ans_v,
                         "w" => ans_w)
        extracted_vars = err.extract_vars(data, vars)
        @test extracted_vars == expected
end

@testset "l2_errors" begin
        u = ones(4)
        v = ones(4)
        H = Matrix(I, 4, 4)
        solution = Dict("u" => u, "v" => v)
        norms = Dict("u" => H, "v" => H)
        mms = Dict("u" => u, "v" => v)
        errors = err.l2_errors(solution, mms, norms)
        expected = Dict("u" => 0.0, "v" => 0.0)
        @test errors == expected
end

@testset "functional_errors" begin
        u = ones(4)
        v = ones(4)
        H = Matrix(I, 4, 4)
        solution = Dict("u" => u, "v" => v)
        norms = Dict("u" => H, "v" => H)
        mms = Dict("u" => 4.0, "v" => 4.0)
        errors = err.functional_errors(solution, mms, norms)
        expected = Dict("u" => 0.0, "v" => 0.0)
        @test errors == expected
end

@testset "log2_convergence_rates" begin
        u = [16, 8, 4, 2]
        q = err.log2_convergence_rates(u)
        ans = [Inf, 1, 1, 1]
        @test q == ans
end

max_stars_repo_name: ooreilly/sbpjl

max_issues_repo_name: ooreilly/sbpjl

max_forks_repo_name: ooreilly/sbpjl


text:
using GeometryTypes, GLVisualize, GLAbstraction, Reactive, GLWindow, GLFW
using ModernGL, Colors
msh = loadasset("cat.obj")
const w = GLPlot.viewing_screen
v = vertices(msh)
f = faces(msh)
colors = RGBA{Float32}[RGBA{Float32}(rand(), rand(), rand(), 1.) for i=1:length(v)]
colored_mesh = GLNormalVertexcolorMesh(
    vertices=v, faces=f,
    color=colors
)
model = translationmatrix(Vec3f0(2, 0, 0))
cat_robj = glplot(colored_mesh, model=model).children[]
indices = decompose(Face{2, GLuint, -1}, colored_mesh)
line_robj = glplot(
    cat_robj[:vertices], :linesegment,
    indices=reinterpret(GLuint, indices),
    thickness=2.5f0, color=RGBA{Float32}(1,1,1,1),
    preferred_camera=:perspective,
    model=model
).children[]
const gpu_position = line_robj[:vertex]

point_robj = glplot(
    (Circle(Point2f0(0), 0.01f0), gpu_position),
    billboard=true,
    preferred_camera=:perspective,
    model=model
).children[]

const m2id = mouse2id(w)
# interaction
@materialize mouse_buttons_pressed, mouseposition = w.inputs
isoverpoint = const_lift(is_same_id, m2id, point_robj)

# single left mousekey pressed (while no other mouse key is pressed)
key_pressed = const_lift(GLAbstraction.singlepressed, mouse_buttons_pressed, GLFW.MOUSE_BUTTON_LEFT)
# righ
right_pressed = const_lift(GLAbstraction.singlepressed, mouse_buttons_pressed, GLFW.MOUSE_BUTTON_RIGHT)
# dragg while key_pressed. Drag only starts if isoverpoint is true
mousedragg = GLAbstraction.dragged(mouseposition, key_pressed, isoverpoint)
camera = w.cameras[:perspective]

# use mousedrag and mouseid + index to actually change the gpu array with the positions
function apply_drag(v0, dragg)
    if dragg == Vec2f0(0) # if drag just started. Not the best way, maybe dragged should return a tuple of (draggvalue, started)
        id, index = value(m2id)
        if id==point_robj.id && length(gpu_position) >= index
            prj_view = value(camera.projectionview)
            p0 = Point4f0(prj_view * Vec4f0(gpu_position[index], 1))
        else
            p0 = v0[3]
        end
    else
        id, index, p0 = v0
        if id==point_robj.id && length(gpu_position) >= index
            prj_view_inv = inv(value(camera.projectionview))
            area = value(camera.window_size)
            cam_res = Vec2f0(widths(area))
            dragg_clip_space = (Vec2f0(dragg)./cam_res) * p0[4] * 2
            pos_clip_space = p0 + Point4f0(dragg_clip_space, 0, 0)
            p_world_space = Point3f0(prj_view_inv * Vec4f0(pos_clip_space))
            gpu_position[index] = p_world_space
        end

    end
    return id, index, p0
end

x = preserve(foldp(apply_drag, (value(m2id)..., Point4f0(0)), mousedragg))


text:
downloaddir = joinpath(dirname(@__FILE__),"downloads")
bindir = joinpath(dirname(@__FILE__),"usr","lib")

const weber_sound_version = 3

# remove any old build files
for d in [downloaddir,bindir]
  rm(d,recursive=true,force=true)
  mkpath(d)
end

################################################################################
# install SDL2 and plugins

SDL2 = "UNKNOWN"
SDL2_ttf = "UNKNOWN"

@static if is_windows()
  # WinRPM lacks SDL2_ttf and SDL2_mixer binaries, so I'm just directly
  # downloading them from the SDL website.
  function setupbin(library,uri)
    libdir = joinpath(downloaddir,library)
    zipfile = joinpath(downloaddir,library*".zip")
    try
      download(uri,zipfile)
      run(`$(joinpath(JULIA_HOME, "7z.exe")) x $zipfile -y -o$libdir`)
      for lib in filter(s -> endswith(s,".dll"),readdir(libdir))
        cp(joinpath(libdir,lib),joinpath(bindir,lib))
      end
      replace(joinpath(bindir,library*".dll"),"\\","\\\\")
    finally
      rm(libdir,recursive=true,force=true)
      rm(zipfile,force=true)
    end
  end
  try
    SDL2 = setupbin("SDL2","https://www.libsdl.org/release/SDL2-2.0.5-win32-x64.zip")
    SDL2_ttf = setupbin("SDL2_ttf","https://www.libsdl.org/projects/SDL_ttf/release/SDL2_ttf-2.0.14-win32-x64.zip")
  finally
    rm(downloaddir,recursive=true,force=true)
  end
elseif is_apple()
  using Homebrew

  Homebrew.add("sdl2")
  Homebrew.add("sdl2_ttf")

  prefix = joinpath(Homebrew.prefix(),"lib")
  SDL2 = joinpath(prefix,"libSDL2-2.0.0.dylib")
  SDL2_ttf = joinpath(prefix,"libSDL2_ttf-2.0.0.dylib")
elseif is_linux()
    error("Weber does not support Linux. You can try manually installing ",
          "SDL2, SDL2_ttf and portaudio and then creating an appropriate ",
          "deps.jl file in $(dirname(@__FILE__)). Be warned however that ",
          "I have encountered strange llvm runtime errors with the linux ",
          "implementation, possibly related to SDL2 problems.")
  # function getdir(lib)
  #   dir = readlines(`dpkg -L $lib-0`)
  #   libs = filter(x -> ismatch(Regex(lowercase(lib)*".so"),lowercase(x)),dir)
  #   chomp(last(libs))
  # end
  # SDL2 = getdir("libsdl2-2.0")
  # SDL2_mixer = getdir("libsdl2-mixer-2.0")
  # SDL2_ttf = getdir("libsdl2-ttf-2.0")
else
  error("Unsupported operating system.")
end

@assert SDL2 != "UNKNOWN"
@assert SDL2_ttf != "UNKNOWN"

deps = joinpath(dirname(@__FILE__),"deps.jl")
open(deps,"w") do s
  for (var,val) in [:weber_SDL2 => SDL2,
                    :weber_SDL2_ttf => SDL2_ttf]
    println(s,"const $var = \"$val\"")
  end
end


text:
## --- Simple linear interpolations

    function _linterp1(x, y, xq::Number, ::Symbol)
        knot_index = searchsortedfirst(x, xq, Base.Order.ForwardOrdering()) - 1
        𝔦₋ = min(max(knot_index, firstindex(x)), lastindex(x) - 1)
        𝔦₊ = 𝔦₋ + 1
        x₋, x₊ = x[𝔦₋], x[𝔦₊]
        y₋, y₊ = y[𝔦₋], y[𝔦₊]
        f = (xq - x₋) / (x₊ - x₋)
        return f*y₊ + (1-f)*y₋
    end

    function _linterp1(x, y, xq::AbstractArray, ::Symbol)
        i₁, iₙ = firstindex(x), lastindex(x) - 1
        knot_index = searchsortedfirst_vec(x, xq) .- 1
        Tₓ = promote_type(eltype(x), eltype(xq))
        T = promote_type(eltype(y), Base.promote_op(/, Tₓ, Tₓ))
        yq = similar(xq, T, size(xq))
        @inbounds for i=1:length(knot_index)
            knot_index[i] = min(max(knot_index[i], i₁), iₙ)
        end
        @inbounds for i=1:length(knot_index)
            𝔦₋ = knot_index[i]
            𝔦₊ = 𝔦₋ + 1
            x₋, x₊ = x[𝔦₋], x[𝔦₊]
            y₋, y₊ = y[𝔦₋], y[𝔦₊]
            f = (xq[i] - x₋)/(x₊ - x₋)
            yq[i] = f*y₊ + (1-f)*y₋
        end
        return yq
    end

    function _linterp1(x, y::AbstractArray{<:AbstractFloat}, xq::AbstractArray, ::Symbol)
        i₁, iₙ = firstindex(x), lastindex(x) - 1
        knot_index = searchsortedfirst_vec(x, xq) .- 1
        Tₓ = promote_type(eltype(x), eltype(xq))
        T = promote_type(eltype(y), Base.promote_op(/, Tₓ, Tₓ))
        yq = similar(xq, T, size(xq))
        @turbo for i=1:length(knot_index)
            knot_index[i] = min(max(knot_index[i], i₁), iₙ)
        end
        @turbo for i=1:length(knot_index)
            𝔦₋ = knot_index[i]
            𝔦₊ = 𝔦₋ + 1
            x₋, x₊ = x[𝔦₋], x[𝔦₊]
            y₋, y₊ = y[𝔦₋], y[𝔦₊]
            f = (xq[i] - x₋)/(x₊ - x₋)
            yq[i] = f*y₊ + (1-f)*y₋
        end
        return yq
    end

    function _linterp1(x, y, xq::Number, extrapolate::Number)
        i₁, iₙ = firstindex(x), lastindex(x) - 1
        knot_index = searchsortedfirst(x, xq, Base.Order.ForwardOrdering()) - 1
        Tₓ = promote_type(eltype(x), eltype(xq))
        T = promote_type(eltype(y), Base.promote_op(/, Tₓ, Tₓ))
        if i₁ <= knot_index <= iₙ
            𝔦₋ = knot_index
            𝔦₊ = 𝔦₋ + 1
            x₋, x₊ = x[𝔦₋], x[𝔦₊]
            y₋, y₊ = y[𝔦₋], y[𝔦₊]
            f = (xq - x₋) / (x₊ - x₋)
            return f*y₊ + (1-f)*y₋
        else
            return T(extrapolate)
        end
    end

    function _linterp1(x, y, xq::AbstractArray, extrapolate::Number)
        i₁, iₙ = firstindex(x), lastindex(x) - 1
        knot_index = searchsortedfirst_vec(x, xq) .- 1
        Tₓ = promote_type(eltype(x), eltype(xq))
        T = promote_type(eltype(y), Base.promote_op(/, Tₓ, Tₓ))
        yq = similar(xq, T, size(xq))
        @inbounds for i=1:length(knot_index)
            𝔦 = knot_index[i]
            if i₁ <= 𝔦 <= iₙ
                𝔦₋ = 𝔦
                𝔦₊ = 𝔦₋ + 1
                x₋, x₊ = x[𝔦₋], x[𝔦₊]
                y₋, y₊ = y[𝔦₋], y[𝔦₊]
                f = (xq[i] - x₋)/(x₊ - x₋)
                yq[i] = f*y₊ + (1-f)*y₋
            else
                yq[i] = extrapolate
            end
        end
        return yq
    end

    # Vectorization-friendly searchsortedfirst implementation from Interpolations.jl
    # https://github.com/JuliaMath/Interpolations.jl
    Base.@propagate_inbounds function searchsortedfirst_exp_left(v, xx, lo, hi)
        for i in 0:4
            ind = lo + i
            ind > hi && return ind
            xx <= v[ind] && return ind
        end
        n = 3
        tn2 = 2^n
        tn2m1 = 2^(n-1)
        ind = lo + tn2
        while ind <= hi
            xx <= v[ind] && return searchsortedfirst(v, xx, lo + tn2 - tn2m1, ind, Base.Order.Forward)
            tn2 *= 2
            tn2m1 *= 2
            ind = lo + tn2
        end
        return searchsortedfirst(v, xx, lo + tn2 - tn2m1, hi, Base.Order.Forward)
    end

    function searchsortedfirst_vec(v::AbstractVector, x::AbstractVector)
        issorted(x) || return searchsortedfirst.(Ref(v), x)
        out = zeros(Int, length(x))
        lo = 1
        hi = length(v)
        @inbounds for i in 1:length(x)
            xx = x[i]
            y = searchsortedfirst_exp_left(v, xx, lo, hi)
            out[i] = y
            lo = min(y, hi)
        end
        return out
    end

## --- Linear interpolation, top-level functions

    """
    ```julia
    yq = linterp1(x::AbstractArray, y::AbstractArray, xq; extrapolate=:Linear)
    ```
    Simple linear interpolation in one dimension. Given a vector of knots `x`
    and values `y`, find the corresponding `y` values at position(s) `xq`.

    Knots `x` must be sorted in increasing order.

    If the optional keyword argument `extrapolate` is set to `:Linear` (default),
    `xq` values outside the range of `x` will be extrapolated using a linear
    extrapolation of the closest two `x`-`y` pairs. Otherwise, if `extrapolate`
    is set to a `Number` (e.g., `0`, or `NaN`), that number will be used
    instead.
    """
    function linterp1(x::AbstractArray, y::AbstractArray, xq; extrapolate=:Linear)
        issorted(x) || error("knot-vector `x` must be sorted in increasing order")
        return _linterp1(x, y, xq, extrapolate)
    end
    export linterp1

    """
    ```julia
    yq = linterp1s(x::AbstractArray, y::AbstractArray, xq; extrapolate=:Linear)
    ```
    As as `linterp1` (simple linear interpolation in one dimension), but will sort
    the knots `x` and values `y` pairwise if `x` if not already sorted in
    increasing order.
    """
    function linterp1s(x::AbstractArray, y::AbstractArray, xq; extrapolate=:Linear)
        sI = sortperm(x) # indices to construct sorted array
        return _linterp1(x[sI], y[sI], xq, extrapolate)
    end
    export linterp1s

    # Linearly interpolate vector y at index i, returning outboundsval if outside of bounds
    function linterp_at_index(y::AbstractArray, i::Number, extrapolate=float(eltype(y))(NaN))
        if firstindex(y) <= i < lastindex(y)
            𝔦₋ = floor(Int, i)
            𝔦₊ = 𝔦₋ + 1
            f = i - 𝔦₋
            return f*y[𝔦₊] + (1-f)*y[𝔦₋]
        else
            return extrapolate
        end
    end
    export linterp_at_index

## --- Resize and interpolate arrays of colors

    # Linearly interpolate array of colors at positions xq
    function linterp1(x::AbstractArray, image::AbstractArray{<:Color}, xq)
        # Interpolate red, green, and blue vectors separately
        r_interp = linterp1(x, image .|> c -> c.r, xq)
        g_interp = linterp1(x, image .|> c -> c.g, xq)
        b_interp = linterp1(x, image .|> c -> c.b, xq)
        # Convert back to a color
        return RGB.(r_interp,g_interp,b_interp)
    end

    function resize_colormap(cmap::AbstractArray{<:Color}, n::Integer)
        cNum = length(cmap)
        if n<2
            cmap[1:1]
        else
            linterp1(1:cNum,cmap,collect(range(1,cNum,length=n)))
        end
    end
    export resize_colormap

max_stars_repo_name: brenhinkeller/StatGeochemBase

max_issues_repo_name: brenhinkeller/StatGeochemBase

max_forks_repo_name: brenhinkeller/StatGeochemBase