text:
# the packages we need
include("CostumeFunctions.jl")

using Gtk, Graphics, Logging, Printf

# the main window
win = GtkWindow("SO(3)")

render_scale = 50

function init_canvas(h,w)
    # create the drawing canvas
    c = @GtkCanvas(h,w)
    # create the initial drawing inside a try/catch loop
    @guarded draw(c) do widget
        #  draw the background with the canvas drawing context
        # the code for this comes later
        draw_the_canvas(c)
    end

    show(c)
    return c
end

# Conversion from meter to pixel and pixel to meter
PIXEL_TO_METER = 0.01
METER_TO_PIXEL = 100

# make the canvas
the_canvas = init_canvas(800,800)

# Our origin
origin2D = [50 600]

# --------- part 2 -------------
# define all the widgets

# bool that determines if we are at camera viewpoint or global viewpoint
viewpoint_global = false

# bool that determines if we show lines between projection and camera
show_lines = true

# Array of original circle points
circle_points = 0

# Array of circle points view from the camera
camera_points = 0

# Array of circle points on the focal plane
spatial_points = 0

V1 = [1 0 0]
V2 = [0 1 0]
V3 = [0 0 1]

function find_by_name(list, name)
    for item in list
        if get_gtk_property(item, :name, String) == name
            return item
        end
    end
    @warn name, "not found in list"
    @warn "available names are"
    for item in list
        @warn get_gtk_property(item, :name, String)
    end
end

# Make a vector of dimensions x an homogenic vector of dimensions x+1
function homog(v, f=1.0)
    if (typeof(f) != Float64) f = parse(Float64, f) end
    push!(v, f)
    return v
end

# Make a vector of dimenions x a dehmog vector of dimensions x-1
function dehomog(v)
    vector = zeros(Float64, length(v)-1)
    for i = 1:length(v)-1
        vector[i] = v[i]/v[length(v)]
    end
    return vector
end

function button_callback(widget)
    name = get_gtk_property(widget, :name, String)

    if name[1] == 'P'
        global viewpoint_global = !viewpoint_global
    else
        global show_lines = !show_lines
    end

    # actually draw the changes
    draw_the_canvas(the_canvas)
    reveal(the_canvas)

end

function button_box(label_string)

    button = GtkButton(label_string)
    set_gtk_property!(button, :name, label_string)
    set_gtk_property!(button, "halign", 3)
    set_gtk_property!(button, "valign", 3)

    signal_connect(button_callback, button, "clicked")

    return button

end

function bold_label(label_string)
    label = GtkLabel("")
    GAccessor.markup(label, """<b>""" * label_string * """</b>""")
    return label
end

function init_window(win, canvas)

    # make a vertically stacked box for the data entry widgets
    control_box = GtkBox(:v)
    push!(control_box, button_box("Press to change view"))
    push!(control_box, button_box("View Lines"))

    # make another box for the drawing canvas
    canvas_box = GtkBox(:v)
    push!(canvas_box, canvas)

    # make a containing box that will stack the widgets and the canvas side by side
    global_box = GtkBox(:v)
    push!(global_box, canvas_box)
    push!(global_box, GtkLabel("   ")) # a very basic separator
    push!(global_box, control_box)
    push!(global_box, GtkLabel("   ")) # a very basic separator

    # put it all inside the window
    push!(win, global_box)

end

# Draw an axis on a point on screen
function draw_axis(origin2D, R1, R2, R3, ctx)
    # Set line
    set_line_width(ctx, 5)
    set_source_rgb(ctx, 0, 0, 1)
    # Move to Origin
    move_to(ctx, origin2D[1], origin2D[2])
    line_to(ctx, origin2D[1] + R1[1] * render_scale, origin2D[2] - R1[2] * render_scale)
    stroke(ctx)
    # Circle to determine direction
    circle(ctx, origin2D[1] + R1[1] * render_scale, origin2D[2] - R1[2] * render_scale, 5)
    fill(ctx)

    # Repeat on the other two vectors
    set_line_width(ctx, 5)
    set_source_rgb(ctx, 1, 0, 0)
    move_to(ctx , origin2D[1], origin2D[2])
    line_to(ctx, origin2D[1] + R2[1] * render_scale, origin2D[2] - R2[2] * render_scale)
    stroke(ctx)

    circle(ctx, origin2D[1] + R2[1] * render_scale, origin2D[2] - R2[2] * render_scale, 5)
    fill(ctx)

    set_line_width(ctx, 5)
    set_source_rgb(ctx, 0, 1, 0)
    move_to(ctx, origin2D[1], origin2D[2])
    line_to(ctx, origin2D[1] + R3[1] * render_scale, origin2D[2] - R3[2] * render_scale)
    stroke(ctx)

    circle(ctx, origin2D[1] + R3[1] * render_scale, origin2D[2] - R3[2] * render_scale, 5)
    fill(ctx)

end

function resize_matrix(A, rows, cols)
    resized = similar(A, rows, cols)
    resized[1:size(A, 1), 1:size(A, 2)] = A
    return resized
end

# Read circle.txt file
function read_circle()
    # Get path:
    path = joinpath(@__DIR__, "")

    input_path = string(path,"circle.txt")

    # Points array
    points = Float64[]

    open(input_path,"r") do file
        # line_number
        line = 0
        # read till end of file
        while ! eof(file)
            # read a new / next line for every iteration
            read = readline(file)
            line += 1

            # context = current txt line context
            context = split(read,"\t")

            # Get current point into points array
            for i = 1:3
                push!(points,parse(Float64, context[i]))
            end

        end
    end

    # Get number of points based on dimension (3)
    size_points = Int(length(points) / 3)

    # Create a matrix of dimensions (size_points, 3)
    global circle_points = zeros(Float64, (size_points, 3))

    counter = 1

    # Put all the points into the circle_points matrix
    for i = 1:size_points
        for j = 1:3
            circle_points[i,j] = points[counter]
            counter = counter + 1
        end
    end
end

function calculate_circle_coordinates()

    # Calculate camera viewpoint
    camera_pos = [1; 6; 1]
    ay = 90
    az = -20
    focal_distance = 1.0

    # Get rotation
    Ry = [cosd(ay) 0 sind(ay);
           0      1      0;
         -sind(ay) 0 cosd(ay)]

    Rz = [cosd(az) -sind(az) 0;
          sind(az) cosd(az)  0;
             0      0      1]

    rotation_matrix = Rz * Ry

    #Create affine camera base
    camera_affine_base = [rotation_matrix camera_pos]
    camera_affine_base = [camera_affine_base;[0 0 0 1]]

    camera_affine_base_inverted = inv(camera_affine_base)

    # Create projection matrix
    projection_matrix = [focal_distance 0.0 0.0 0.0; 0.0 focal_distance 0.0 0.0; 0.0 0.0 1.0 0.0]

    # Create matrix of the same dimensions as circle_points to store the points on the plane
    global spatial_points = copy(circle_points)

    # Create a matrix of the size of circle_points but 2 dimensional to store the points seen from camera
    global camera_points = zeros(Float64, (size(spatial_points)[1], 2))

    # Calculate camera points and spatial_points
    for i = 1:size(circle_points)[1]
        # Camera Points-----------
        point = [copy(circle_points[i,1]); copy(circle_points[i,2]); copy(circle_points[i,3])]
        point = dehomog(projection_matrix * camera_affine_base_inverted * homog(point))

        for k = 1:2
            camera_points[i,k] = point[k]
        end

        # Spatial points ------------------------------------
        spatial_point = dehomog(camera_affine_base*homog(homog(point, focal_distance)))
        for j = 1:3
            spatial_points[i,j] = spatial_point[j]
        end
    end

end

# Draw the scene in 3D
function draw_canvas_global(ctx)

    # Draw a global axis---------------
    R1 = [10 0 0]
    R2 = [0 10 0]
    R3 = [0 0 10]

    R1 = to_2d(R1)
    R2 = to_2d(R2)
    R3 = to_2d(R3)

    draw_axis(origin2D, R1, R2, R3, ctx)

    # Draw circle on camera plane-----------------------------
    for i = 1:size(spatial_points)[1]
        point = [spatial_points[i,1]; spatial_points[i,2]; spatial_points[i,3]]
        point = to_2d(point)

        set_source_rgb(ctx, 1, 0, 0)
        circle(ctx, origin2D[1] + (point[1] * METER_TO_PIXEL), origin2D[2] + (point[2] * METER_TO_PIXEL), 1)
        fill(ctx)
    end

    # Draw original circle --------------------------------------
    for i = 1:size(circle_points)[1]

        original_point = [circle_points[i,1]; circle_points[i,2]; circle_points[i,3]]
        original_point = to_2d(original_point)
        set_source_rgb(ctx, 0, 0, 1)
        circle(ctx, origin2D[1] + (original_point[1] * METER_TO_PIXEL), origin2D[2] + (original_point[2] * METER_TO_PIXEL), 1)
        fill(ctx)

        # Draw lines from circle points to camera position
        if (show_lines == true)
            camera_pos = to_2d([1;6;1])
            set_line_width(ctx, 1)
            # Move to Origin
            move_to(ctx, origin2D[1] + camera_pos[1] * METER_TO_PIXEL, origin2D[2] + camera_pos[2] * METER_TO_PIXEL)
            line_to(ctx, origin2D[1] + original_point[1] * METER_TO_PIXEL, origin2D[2] + original_point[2] * METER_TO_PIXEL)
            stroke(ctx)
        end
    end

    # Draw camera axis------------------------------------------------------
    camera_pos = to_2d([1;6;1])

    # Get camera pos in 2D
    camera_pos_pixels = [origin2D[1] + (camera_pos[1] * METER_TO_PIXEL); origin2D[2] + (camera_pos[2] * METER_TO_PIXEL)]

    # Get camera axis vectors
    V1 = [1;0;0]
    V2 = [0;1;0]
    V3 = [0;0;1]

    # Rotate vectors on the Y axis
    V1 = rotate_phi(V1, 90, [0;1;0])
    V2 = rotate_phi(V2, 90, [0;1;0])
    V3 = rotate_phi(V3, 90, [0;1;0])

    # Rotate vectors on the resultant Z axis
    V1 = rotate_phi(V1, -20, V3)
    V2 = rotate_phi(V2, -20, V3)
    V3 = rotate_phi(V3, -20, V3)

    V1 = to_2d(V1)
    V2 = to_2d(V2)
    V3 = to_2d(V3)

    # Draw the axis
    draw_axis(camera_pos_pixels, V1, V2, V3, ctx)

end

# Draw the scene in 2D (camera point of view)
function draw_canvas_camera(ctx)
    # Draw circle seen from camera
    for i = 1:size(camera_points)[1]
        set_source_rgb(ctx, 1, 0, 0)
        circle(ctx, 700 + camera_points[i,1] * render_scale * 20, camera_points[i,2] * render_scale * 20  , 1)
        fill(ctx)
    end
end

function draw_the_canvas(canvas)

    h   = height(canvas)
    w   =  width(canvas)
    ctx =  getgc(canvas)
    # clear the canvas
    rectangle(ctx, 0, 0, w, h)
    set_source_rgb(ctx, 1, 1, 1)
    fill(ctx)

    # Check if we are seeing the 3D or 2D point of view
    if (viewpoint_global == true)
        draw_canvas_global(ctx)
    else
        draw_canvas_camera(ctx)
    end
end
println("----------------------------------START--------------------------------")

# prepare and  the initial widgets
init_window(win, the_canvas)
showall(win)
read_circle()
calculate_circle_coordinates()


text:
import Base: show, display, hash, ==, isequal, eltype, keys, values
import Base: iterate, length # Dict stuff

# ========== Mapping/Predicate =========

# TODO: limit types for arguments of f
struct Predicate
  f::Function
end

struct Mapping
  f::Function
end

"""
    Predicate(f)(x,y)
Functor definition for Predicate. This allows to use the type Predicate as a function that is guaranteed to
have two arguments and return a Bool.
# Examples
```
Predicate((x,y) -> x>y)(3,4)
julia> true
p = Predicate((x,y) -> x>y);
p(4,3)
julia> false
```
"""

((p::Predicate)(x::T, y::T)::Bool) where T = p.f(x,y)

"""
    Mapping(f)(x,y)
Functor definition for Mapping. This allows to use the type Mapping as a function that is guaranteed to
have two arguments and return a Bool.
# Examples
```
Mapping(x -> [x,2*x])(3)
julia> [3,6]
Mapping(x -> 2:x)(30)
julia> 2:30
```
"""
((m::Mapping)(x::T)::AbstractArray{T,1}) where T = m.f(x)

# ========== ReduceMap =========
"""
    ReduceMap

"""
struct ReduceMap{K, V} <: AbstractDict{K,V}
    map::OrderedDict{K,V}
    function ReduceMap(d::OrderedDict{K1,V1}) where {K1,V1}
        ReduceMap_TConstr(K1,V1)
        new{K1,V1}(d, OrderedDict{K1,UInt32}())
    end
end

function ReduceMap_TConstr(::Type{K}, ::Type{V}) where {K, V}
    if !((K === V) || ((V === UInt32) || (V === Int64) || (V == Tuple{UInt32, UInt32})))
        throw(TypeError(:ReduceMap,"Key and Value must have the same type, OR Value must be UInt32/Int64", Union{K,UInt32}, V))
     end
     return nothing
end

ReduceMap(args...) = ReduceMap(OrderedDict(args...))
ReduceMap{K,V}() where {K,V} = ReduceMap(OrderedDict{K,V}())

# ---------- Dict Wrapper  ---------
@forward ReduceMap.map (Base.show, Base.length, Base.iterate, Base.getindex, Base.setindex!)

# ---------- Auxiliary Functions  ---------

function toIndexMap(rm::ReduceMap{T,T}, orderedSet) where T
    res = ReduceMap{T,UInt32}()
    for (K,V) in rm
        res[K] = UInt32(findfirst(x->x==V, orderedSet))
    end
    return res
end

function toDirectMap(rm::ReduceMap{T,UInt32}, orderedSet) where T
    res = ReduceMap{T,T}()
    for (K,V) in rm
        res[K] = orderedSet[rm[K]]
    end
    return res
end

# ========== ExpandMap =========
"""
    ExpandMap(rm::ReduceMap [,sorted])

Constructs mapping from a set of representatives of classes to full set using
a ReduceMap `rm`.
Each lit of mappings to the full can be ordered by setting `sorted`.
"""
struct ExpandMap{K,V <: Union{K, Int64, UInt32}}
    map::Dict{V,Array{K,1}}
    function ExpandMap(rm::ReduceMap{K1,V1}, orderedSet=nothing; sorted=false) where {K1,V1}
        if (V1 == UInt32 || V1 == Int64) && orderedSet != nothing
            new{K1,K1}(invertDict(toDirectMap(rm, orderedSet).map, sorted=sorted))
        else
            new{K1,V1}(invertDict(rm.map, sorted=sorted))
        end
    end
end

# ---------- Dict Wrapper  ---------
@forward ExpandMap.map (Base.keys, Base.values, Base.show, Base.length, Base.iterate, Base.getindex, Base.setindex!)

# --------------------  Auxiliary functions --------------------

Base.hash(a::ExpandMap, h::UInt) = hash(a.map, hash(:ExpandMap, h))
Base.:(==)(a::ExpandMap{T}, b::Dict{T,Array{T,1}}) where T = (a.map == b)
Base.:(==)(a::ExpandMap, b::ExpandMap) = (a.map == b.map)

# ====================  Auxiliary functions ====================

function display(el::ExpandMap)
    show(stdout, el)
end

function show(io::IO,el::ExpandMap)
    show(io, "ExpandMap{$(eltype(el))}["*string(length(el.map))*"]:\n"*string(el.map))
end

minimal_set(m::ExpandMap; sorted=true) = sorted ? sort(collect(keys(m))) : collect(keys(m))
minimal_set(m::ReduceMap; sorted=true) = sorted ? sort(collect(unique(values(m)))) : collect(unique(values(m)))

"""
    labelsMap(m::ReduceMap)
Transforms values to set of continuous numbers.

# Examples
label(ReduceMap(1=>1, 2=>3, 3=>5))
julia> ReduceMap{Int64,Int64} with 3 entries:
  2 => 3
  3 => 5
  1 => 1
"""
function labelsMap(m::ReduceMap)::ReduceMap
    vv = unique(values(m))
    lookup = OrderedDict{keytype(m), UInt32}(zip(vv, 1:length(vv)))
    res = OrderedDict{keytype(m), UInt32}()
    for (k,v) in m
        res[k] = lookup[v]
    end
    return ReduceMap(res)
end

ext: jl


text:
println("test mel-generalized cepstrum analysis")

function test_mgcep_base(f::Function, order, args...; kargs...)
    srand(98765)
    dummy_input = rand(256)
    @assert applicable(f, dummy_input, order, args...)
    c = f(dummy_input, order, args...)
    @test length(c) == order + 1
    @test all(isfinite.(c))
end

function test_mcep_and_mgcep_consistency(order)
    srand(98765)
    dummy_input = rand(256)
    # mgcep when gamma = 0, the result of mgcep is corresponds to mcep
    mc = mcep(dummy_input, order, 0.41)
    mgc = mgcep(dummy_input, order, 0.41, 0.0)
    @test isapprox(mc, mgc, atol=1.0e-3)
end

function test_mcep_exceptions()
    dummy_input = ones(256)
    @test_throws ArgumentError mcep(dummy_input, itype=-1)
    @test_throws ArgumentError mcep(dummy_input, itype=5)
    @test_throws ArgumentError mcep(dummy_input, eps=-1.0)
    @test_throws ArgumentError mcep(dummy_input, etype=-1)
    @test_throws ArgumentError mcep(dummy_input, etype=-3)
    @test_throws ArgumentError mcep(dummy_input, etype=1, eps=-1.0)
    @test_throws ArgumentError mcep(dummy_input, etype=2, eps=-1.0)
    @test_throws ArgumentError mcep(dummy_input, min_det=-1.0)

    # should have error in computing log periodogram
    @test_throws ErrorException mcep(ones(256), 40, 0.41)
end

function test_gcep_exceptions()
    dummy_input = ones(256)

    # invalid γ
    @test_throws ArgumentError gcep(dummy_input, 40, 0.1)
    @test_throws ArgumentError gcep(dummy_input, 40, -2.1)

    # invalid optinal paramters
    @test_throws ArgumentError gcep(dummy_input, itype=-1)
    @test_throws ArgumentError gcep(dummy_input, itype=5)
    @test_throws ArgumentError gcep(dummy_input, eps=-1.0)
    @test_throws ArgumentError gcep(dummy_input, etype=-1)
    @test_throws ArgumentError gcep(dummy_input, etype=-3)
    @test_throws ArgumentError gcep(dummy_input, etype=1, eps=-1.0)
    @test_throws ArgumentError gcep(dummy_input, etype=2, eps=-1.0)
    @test_throws ArgumentError gcep(dummy_input, min_det=-1.0)

    # should have error in theq
    @test_throws ErrorException gcep(ones(256), 40, 0.0)
end

function test_mgcep_exceptions()
    dummy_input = ones(256)

    # invalid γ
    @test_throws ArgumentError mgcep(dummy_input, 40, 0.41, 0.1)
    @test_throws ArgumentError mgcep(dummy_input, 40, 0.41, -2.0)

    # invalid optinal paramters
    @test_throws ArgumentError mgcep(dummy_input, itype=-1)
    @test_throws ArgumentError mgcep(dummy_input, itype=5)
    @test_throws ArgumentError mgcep(dummy_input, eps=-1.0)
    @test_throws ArgumentError mgcep(dummy_input, etype=-1)
    @test_throws ArgumentError mgcep(dummy_input, etype=-3)
    @test_throws ArgumentError mgcep(dummy_input, etype=1, eps=-1.0)
    @test_throws ArgumentError mgcep(dummy_input, etype=2, eps=-1.0)
    @test_throws ArgumentError mgcep(dummy_input, min_det=-1.0)
    @test_throws ArgumentError mgcep(dummy_input, otype=-1)
    @test_throws ArgumentError mgcep(dummy_input, otype=6)

    # should have error in theq
    @test_throws ErrorException mgcep(ones(256))
end

function test_uels_exceptions()
    dummy_input = ones(256)

    # invalid optinal paramters
    @test_throws ArgumentError uels(dummy_input, itype=-1)
    @test_throws ArgumentError uels(dummy_input, itype=5)
    @test_throws ArgumentError uels(dummy_input, eps=-1.0)
    @test_throws ArgumentError uels(dummy_input, etype=-1)
    @test_throws ArgumentError uels(dummy_input, etype=-3)
    @test_throws ArgumentError uels(dummy_input, etype=1, eps=-1.0)
    @test_throws ArgumentError uels(dummy_input, etype=2, eps=-1.0)

    # should have error in computing log periodogram
    @test_throws ErrorException uels(ones(256), 40)
end

function test_lpc_exceptions()
    @test_throws ArgumentError lpc(ones(256), 40, min_det=-1.0)
    @test_throws ErrorException lpc(zeros(256), 40)
end

println("-- test_mcep")
for order in [20, 22, 24]
    for α in [0.35, 0.41, 0.5]
        println(" where order = $order, α = $α")
        test_mgcep_base(mcep, order, α)
    end
end

test_mcep_exceptions()

println("-- test_gcep")
for order in [20, 22, 24]
    for γ in [-1.0, -0.5, 0.0]
        println(" where order = $order, γ = $γ")
        test_mgcep_base(gcep, order, γ)
    end
end

test_gcep_exceptions()

println("-- test_mgcep")
for order in [15, 20, 25, 30]
    for α in [0.35, 0.41, 0.5]
        for γ in [-1.0, -0.5, 0.0]
            for otype in 1:5
                println(" where order = $order, α = $α, γ = $γ, otype=$otype")
                test_mgcep_base(mgcep, order, α, γ, otype=otype)
            end
        end
    end
end

test_mgcep_exceptions()

println("-- test_fftcep")
for order in [20, 22, 24]
    println(" where order = $order")
    test_mgcep_base(fftcep, order)
end

println("-- test_uels")
for order in [20, 22, 24]
    println(" where order = $order")
    test_mgcep_base(uels, order)
end

test_uels_exceptions()

println("-- test_lpc")
for order in [20, 22, 24]
    println(" where order = $order")
    test_mgcep_base(lpc, order)
end

test_lpc_exceptions()

println("-- test mcep and mgcep consistency")
for order in [20, 22, 24]
    println(" where order = $order")
    test_mcep_and_mgcep_consistency(order)
end


text:
# to be used by Run_UniformGeneration.py

using Pkg
Pkg.activate(".")
# you can comment the command below after running this script for the first time
Pkg.instantiate()

include("CountingMatchings_rosso.jl")
include("UniformGeneration_julia.jl")

# pre compilation with easy example
G,E = Read_Graph("auxiliar_instance")
time_cpu = @elapsed Uniform_Number(G,E,2,0)

# filename given by python
G,E = Read_Graph("Graph_file")
f = open("size")
lines =  readlines(f)
close(f)
lines=split(lines[1])
# gen_NUM is the K, i.e, the number of optimal matchings to be generated
gen_NUM =parse(Int,lines[1])
# gen_NUM given by python
time_cpu = @elapsed Uniform_Number(G,E,gen_NUM,1)
time_cpu = @sprintf "%.2f" time_cpu

# save computational results
fres = open("results.txt","a")
write(fres,"time_Uniform_Gen= ",time_cpu, " \n")
close(fres)

max_stars_repo_name: mxmmargarida/KEG

max_issues_repo_name: mxmmargarida/KEG

max_forks_repo_name: mxmmargarida/KEG


text:
module BayesianModels

using ExpFamilyDistributions

export Model
export BayesianParameter
export AbstractConjugateParameter
export ConjugateParameter
export HierarchicalConjugateParameter
export elbo
export getconjugateparams

"""
    BayesianParameter

Supertype for all Bayesian parameters, i.e. parameter having a prior
and a posterior.
"""
abstract type BayesianParameter end

"""
    AbstractConjugateParameter <: BayesianParameter

Bayesian parameter having the same type of prior and posterior.
"""
abstract type AbstractConjugateParameter <: BayesianParameter end

Base.show(io::IO, p::AbstractConjugateParameter) = print(io, "$(typeof(p))")

"""
    ConjugateParameter <: AbstractConjugateParameter

Standard implementation of a conjugate parameter.
"""
mutable struct ConjugateParameter{T <: ExpFamilyDistribution} <: AbstractConjugateParameter
    prior::T
    posterior::T
    stats::Function

    function ConjugateParameter(prior::ExpFamilyDistribution,
                                post::ExpFamilyDistribution)
        new{typeof(post)}(prior, post)
    end
end

"""
    getconjugateparams(model)

Return the list of the conjugate parameters of an object.
"""
function getconjugateparams(model)
    params = []
    for name in fieldnames(typeof(model))
        field = getfield(model, name)
        if typeof(field) <: ConjugateParameter
            push!(params, field)
        end
    end
    return params
end

"""
    HierarchicalConjugateParameter <: AbstractConjugateParameter
Σ
Conjugate parameter whose prior depends on another parameter.
"""
mutable struct HierarchicalConjugateParameter{T <: ExpFamilyDistribution} <: AbstractConjugateParameter
    getprior::Function
    posterior::T
    stats::Function

    function HierarchicalConjugateParameter(getprior::Function,
                                            post::ExpFamilyDistribution)
        new{typeof(post)}(getprior, post)
    end
end

function Base.getproperty(p::HierarchicalConjugateParameter, sym::Symbol)
    if sym ≡ :prior
        return p.getprior()
    else
        return getfield(p, sym)
    end
end

"""
    Model

Supertype for JuBeer models.
"""
abstract type Model end

# Imitate pytorch printing style of model.
function Base.show(io::IO, model::Model; pad = 0, padshift = 2, prefix = "")
    originalpad = pad
    modelname = "$(prefix)$(typeof(model))("
    println(io, lpad(modelname, length(modelname) + pad))
    pad += padshift
    for name in propertynames(model)
        pval = getproperty(model, name)
        if isa(pval, Model)
            pstr = "($(name)): "
            print(io, lpad(pstr, length(pstr) + pad))
            Base.show(io, pval, pad=pad)
        elseif isa(pval, BayesianParameter)
            pstr = "($(name)): $(pval)"
            println(io, lpad(pstr, length(pstr) + pad))
        elseif isa(pval, AbstractVector{<:Model})
            pstr = "($(name)): ["
            println(io, lpad(pstr, length(pstr) + pad))
            for (i, m) in enumerate(pval)
                Base.show(io, m, pad = pad + padshift, prefix = "($i) ")
            end
            pstr = "]"
            println(io, lpad(pstr, length(pstr) + pad))
        end
    end
    println(io, lpad(")", 1 + originalpad))
    return nothing
end

"""
    elbo(model)

Compute the evidence lower-bound of the model.
"""
function elbo(model::Model, X...)
    KL = sum([kldiv(param.posterior, param.prior)
              for param in getconjugateparams(model)])
    return sum(model(X...)) - KL
end

include("LogisticRegression.jl")
include("ARNormals.jl")

# Helper functions for plotting parameters of the models.
include("Plotting.jl")

end # BayesianModels

max_stars_repo_name: fnlandini/BayesianModels

max_issues_repo_name: fnlandini/BayesianModels

max_forks_repo_name: fnlandini/BayesianModels


text:
function plmdca_sym(Z::Array{T,2},W::Vector{Float64};
                decimation::Bool=false,
                fracmax::Float64 = 0.3,
                fracdec::Float64 = 0.1,
                blockdecimate::Bool=true, # decimate on a per-block base (all J[:,:,i,j] = 0)
                remove_dups::Bool = true,
                min_separation::Int = 1,
                max_gap_fraction::Real = 0.9,
                theta = :auto,
                lambdaJ::Real=0.01,
                lambdaH::Real=0.01,
                epsconv::Real=1.0e-5,
                maxit::Int=1000,
                verbose::Bool=true,
                method::Symbol=:LD_LBFGS) where T<: Integer

    all(x->x>0,W) || throw(DomainError("vector W should normalized and with all positive elements"))
    isapprox(sum(W),1) || throw(DomainError("sum(W) ≠ 1. Consider normalizing the vector W"))

    N,M = size(Z)
    M = length(W)
    q = Int(maximum(Z))

    plmalg = PlmAlg(method,verbose, epsconv ,maxit)
    plmvar = PlmVar(N,M,q,q*q,lambdaJ,lambdaH,Z,W)

    Jmat,pseudolike = if decimation == false
        MinimizePLSym(plmalg,plmvar)
    else
        if blockdecimate
            decvar = DecVar{1}(fracdec, fracmax, blockdecimate, ones(Bool, binomial(N,2)))
        else
            decvar = DecVar{1}(fracdec, fracmax, blockdecimate, ones(Bool, binomial(N,2)*q*q+N*q))
        end
        DecimateSym!(plmvar, plmalg, decvar)
    end
    score, Jtens, htens = ComputeScoreSym(Jmat, plmvar, min_separation)

    return PlmOut(pseudolike, Jtens, htens, score)
end

function plmdca_sym(filename::String;
                theta::Union{Symbol,Real}=:auto,
                max_gap_fraction::Real=0.9,
                remove_dups::Bool=true,
                kwds...)
    W,Z,N,M,q = ReadFasta(filename,max_gap_fraction, theta, remove_dups)
    plmdca_sym(Z,W; kwds...)
end

function MinimizePLSym(alg::PlmAlg, var::PlmVar)

    N  = var.N
    q  = var.q
    q2 = var.q2
    Z = var.Z
    Nc2 = binomial(N,2)
    LL  = Nc2 * q2  + N * q

    x0 = zeros(Float64, LL)
    batchidx = myrange(Z)
    opt = Opt(alg.method, length(x0))
    ftol_abs!(opt, alg.epsconv)
    ftol_rel!(opt, alg.epsconv)
    xtol_rel!(opt, alg.epsconv)
    xtol_abs!(opt, alg.epsconv)
    maxeval!(opt, alg.maxit)
    min_objective!(opt, (x,g)->like_grad!(g,x,var,batchidx,alg.verbose))
    elapstime = @elapsed  (minf, minx, ret) = optimize(opt, x0)
    alg.verbose && @printf("pl = %.4f\t time = %.4f\t exit status = ", minf, elapstime)
    alg.verbose && println(ret)

    return minx, minf
end

function myrange(q::DenseMatrix)
    nchunks = nprocs()
    nchunks > 2 || (return [1:size(q,2),])
    splits = [round(Int, s) for s in range(0, stop=size(q,2), length=nchunks)]
    [splits[idx]+1:splits[idx+1] for idx in 1:nchunks-1]
end

reduce_res(x::Tuple) = x
(reduce_res(x::Vector{T}) where T<:Tuple) = broadcast(+,x...)

function L2norm_sym(vec::AbstractVector, var::PlmVar)

    q = var.q
    N = var.N
    lambdaJ = var.lambdaJ
    lambdaH = var.lambdaH

    LL = length(vec)

    mysum1 = 0.0
    @inbounds @simd for i in 1:(LL-N*q)
        mysum1 += vec[i] * vec[i]
    end
    mysum1 *= lambdaJ

    mysum2 = 0.0
    @inbounds @simd for i in (LL - N*q + 1):LL
        mysum2 += vec[i] * vec[i]
    end
    mysum2 *= 2lambdaH
# 	mysum2 *= lambdaH
    return mysum1+mysum2
end

function add_l2grad!(grad::Vector,vecJ::AbstractVector,plmvar::PlmVar)
    N = plmvar.N
    LL = length(grad)
    q = plmvar.q
    lambdaJ = plmvar.lambdaJ
    lambdaH = plmvar.lambdaH
    @inbounds begin
        @simd for i in 1:LL-N*q
            grad[i] += 1.0 * vecJ[i] * lambdaJ
        end
        @simd for i in (LL-N*q + 1):LL
            grad[i] += 2.0 * vecJ[i] * lambdaH
#			grad[i] += 2.0 * vecJ[i] * lambdaH
        end
    end
    nothing
end
function like_grad!(g::Vector,vecJ::AbstractVector,plmvar::PlmVar,vec_chunk,verbose)

    g === nothing && (g = zero(vecJ))
    sJ = SharedArray(vecJ)
    res = @distributed vcat for chunk in vec_chunk
        pl_chunk, gr_chunk = plm_site_grad(sJ, plmvar, chunk)
    end
    pl,gr = reduce_res(res)
    pl += L2norm_sym(vecJ, plmvar)
    verbose && println("pl = $pl")
    add_l2grad!(gr,vecJ, plmvar)
    g .= gr
    return pl
end

function plm_site_grad(vecJ::AbstractVector, plmvar::PlmVar,chunk)

    LL = length(vecJ)
    q2 = plmvar.q2
    q = plmvar.q
    N = plmvar.N
    M = plmvar.M
    Z = plmvar.Z
    W = plmvar.W
    grad = zero(vecJ)
    pseudolike = 0.0
    for a in chunk
        pseudolike += ComputePatternPLSym!(grad, vecJ, Z[:,a],W[a], N, q, q2)
    end

    return pseudolike,grad
end

function ComputePatternPLSym!(grad::Array{Float64,1}, vecJ::AbstractVector, Z::AbstractArray{Int,1}, Wa::Float64, N::Int, q::Int, q2::Int)
    vecene = zeros(Float64,q)
    expvecenesunorm = zeros(Float64,q)
    pseudolike = 0.0
    offset = mygetindex(N-1, N, q, q, N, q, q2)

    @inbounds for site=1:N    # site < i
        fillvecenesym!(vecene, vecJ, Z, site, q,N)

  	    norm = sumexp(vecene)
        expvecenesunorm .= exp.(vecene .- log(norm))
        pseudolike -= Wa * ( vecene[Z[site]] - log(norm) )
		@simd for i = 1:(site-1)
             for s = 1:q
                grad[ mygetindex(i, site, Z[i], s, N, q, q2) ] += 0.5 * Wa * expvecenesunorm[s]
            end
            grad[ mygetindex(i, site , Z[i], Z[site],  N,q,q2)] -= 0.5 * Wa
        end
		@simd for i = (site+1):N
             for s = 1:q
                grad[ mygetindex(site, i , s,  Z[i], N,q,q2) ] += 0.5 * Wa * expvecenesunorm[s]
            end
            grad[ mygetindex(site, i , Z[site], Z[i], N,q,q2)] -= 0.5* Wa
        end
        @simd for s = 1:q
            grad[ offset + s ] += Wa *  expvecenesunorm[s]
        end
		grad[ offset + Z[site] ] -= Wa
 		offset += q
    end
    return pseudolike
end

function fillvecenesym!(vecene::AbstractArray, vecJ::AbstractVector, Z::AbstractArray{Int64,1}, site::Int, q::Int ,N::Int)

    q2 = q*q
    LL=length(vecJ)

    @inbounds begin
        @simd for l = 1:q
            offset::Int = 0
            scra::Float64 = 0.0

            for i=1:site-1
                scra += vecJ[ mygetindex(i, site, Z[i], l,  N, q, q2)]
            end
    	    for i = site+1:N
                scra += vecJ[ mygetindex(site, i, l, Z[i], N, q, q2)]
            end # End sum_i \neq site J
           	#scra *= 0.5
            #offset = mygetindex(N-1, N, q, q, N, q, q2)  + ( site - 1) * q  # last J element + (site-1)*q

            scra += vecJ[LL + ( site - 1 - N ) * q + l] # sum H
            vecene[l] = scra
        end
    end
end

@inline function mygetindex( i::Int, j::Int, coli::Int, colj::Int, N::Int, q::Int, q2::Int)
    offset_i = ( (i-1) * N  - ( (i * ( i -1 ) ) >> 1 ) )  # (i-1) N q2 + i (i-1) q2 / 2
    offset_j = (j - i - 1 )
    return ( offset_i + offset_j ) * q2 + coli + q * (colj - 1)
end

function ComputeScoreSym(Jvec::Array{Float64,1}, var::PlmVar, min_separation::Int)

    LL = length(Jvec)
    N=var.N
    q=var.q
    Nc2 = binomial(N,2)

    Jtens=reshape(Jvec[1:LL-N*q],q,q,Nc2)
    htens=fill(0.0,q,N)
    J = zeros(q,q,Nc2)

    htens=reshape(Jvec[LL-N*q + 1:end],q,N)

    for l=1:Nc2
        J[:,:,l] = Jtens[:,:,l] -
					repeat(mean(Jtens[:,:,l],dims=1),q,1) -
					repeat(mean(Jtens[:,:,l],dims=2),1,q) .+
					mean(Jtens[:,:,l])
    end

    FN = zeros(Float64, N,N)
    l = 1
    for i=1:N-1
        for j=i+1:N
            FN[i,j] = norm(J[:,:,l],2)
            FN[j,i] = FN[i,j]
            l+=1
        end
    end
    FN = correct_APC(FN)
    score = compute_ranking(FN,min_separation)
	GC.gc() # something wrong with SharedArray
    return score, inflate_matrix(Jtens,N),htens
end

max_stars_repo_name: giovact/PlmDCA

max_issues_repo_name: giovact/PlmDCA

max_forks_repo_name: giovact/PlmDCA


text:
using CSV, ModelSelection

data = CSV.read("visitors.csv")

data = ModelSelection.gsr(
    "australia china japan uk",
    data,
    intercept=true,
    time=:date,
    fe_sqr=[:uk, :china],
    fe_log=[:japan],
    fe_inv=:uk,
    preliminaryselection=:lasso,
    outsample=10,
    criteria=[:aic, :aicc],
    ttest=true,
    modelavg=true,
    residualtest=true,
    orderresults=true,
    kfoldcrossvalidation=true,
    numfolds=10,
    exportcsv="visitors_output.csv",
    exportlatex="latex.zip"
)


text:
using MetidaStats
using MetidaBase
using StatsBase
using Test
using DataFrames, CSV

path     = dirname(@__FILE__)
io       = IOBuffer();

@testset "  Descripitive statistics                                  " begin
    ds  = CSV.File(path*"/csv/ds.csv") |> DataFrame
    di  = MetidaStats.dataimport(ds, vars = [:var1, :var2], sort = [:col, :row])
    des = MetidaStats.descriptives(di; skipmissing = true, skipnonpositive = true, stats = MetidaStats.STATLIST)

    di  = MetidaStats.dataimport(ds, vars = [:var1, :var2])

    mt  = MetidaBase.metida_table(des; stats = [:mean, :geom])

    mt  = MetidaBase.metida_table(des; stats = [:mean, :geom], id = [:Variable,:row])

    des2 = MetidaStats.descriptives(ds, [:var1, :var2], [:col, :row]; skipmissing = true, skipnonpositive = true, stats = MetidaStats.STATLIST)

    @test des[:, :mean] == des2[:, :mean]

    show(io, des)

    sort!(des2, [:col, :row, :Variable])

    @test des2[:, :n]  == [24
    22
    44
    44
    20
    20
    83
    83]

    @test des2[:, :mean]  ≈ [51.84339580144830000
    -3.2427526832130800
    58.06263961606340000
    1.8800364385656400
    51.84107641126520000
    0.4353627340913350
    47.25778648204360000
    -3.2515999935069800]

    @test des2[:, :sd]  ≈ [30.67890817616400000
    26.7334282254075000
    26.95184403646440000
    28.9591749375126000
    25.29979111399270000
    27.5437061096263000
    27.16599075359280000
    28.8185876098308000]

    @test des2[:, :se]  ≈ [6.26230590810839000
    5.6995860482981700
    4.06314336714229000
    4.3657598866360500
    5.65720527474328000
    6.1589599213400700
    2.98185487195489000
    3.1632509429411400]

    @test des2[1, :geom]  ≈ 39.7551197516893
    @test des2[3, :geom]  ≈ 48.4385401710072
    @test des2[5, :harmmean] ≈ 38.7909224529887
    @test des2[7, :harmmean] ≈ 15.394283582287837

    di  = MetidaStats.dataimport(ds, vars = [:var1, :var2], sort = [:col, :row])
    sort!(di, [:col, :row, :Variable])

    des2[1, :skew] ≈ skewness(di[1].obs)

    des2[1, :kurt] ≈ kurtosis(di[1].obs)

end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    @assert precompile(HyperRectangle{2,Float32}, (Int, Int, Int, Int))
    @assert precompile(==, (HyperRectangle{2,Float32}, HyperRectangle{2,Float32}))
    @assert precompile(normal_mesh, (Tesselation{3,Float32,Cylinder{3,Float32},1},))
    @assert precompile(normal_mesh, (Tesselation{3,Float32,HyperSphere{3,Float32},1},))
    @assert precompile(normal_mesh, (HyperSphere{3,Float32},))
    @assert precompile(Tuple{typeof(*),SMatrix{4, 4, Float32, 16},HyperRectangle{3, Float32}})   # time: 0.11091917
    @assert precompile(Tuple{typeof(coordinates),HyperRectangle{2, Float32},Tuple{Int64, Int64}})   # time: 0.08693867
    @assert precompile(union, (HyperRectangle{3, Float32}, HyperRectangle{3, Float32}))
    @assert precompile(Tuple{typeof(decompose),Type{Point{2, Float32}},HyperRectangle{2, Float32}})   # time: 0.026609203
    @assert precompile(Tuple{Type{HyperRectangle{3, Float32}},HyperRectangle{2, Float32}})   # time: 0.023717888
    @assert precompile(Tuple{typeof(+),HyperRectangle{3, Float32},Point{3, Float32}})   # time: 0.006633118
    @assert precompile(Tuple{Type{Rect2D{T} where T},Float32,Float32,Float32,Float32})   # time: 0.001636267
    @assert precompile(Tuple{typeof(*),HyperRectangle{2, Float32},Float32})   # time: 0.001057589

    if Base.VERSION >= v"1.6.0-DEV.1083"
        @assert precompile(triangle_mesh, (Polygon{2, Float32, Point2f0, LineString{2, Float32, Point2f0,
                                        Base.ReinterpretArray{Line{2, Float32}, 1, Tuple{Point2f0, Point2f0}, TupleView{Tuple{Point2f0, Point2f0}, 2, 1, Vector{Point2f0}}, false}},
                                        Vector{LineString{2, Float32, Point2f0, Base.ReinterpretArray{Line{2, Float32}, 1, Tuple{Point2f0, Point2f0}, TupleView{Tuple{Point2f0, Point2f0}, 2, 1, Vector{Point2f0}}, false}}}},))
    else
        @assert precompile(triangle_mesh, (Polygon{2, Float32, Point2f0, LineString{2, Float32, Point2f0,
                                        Base.ReinterpretArray{Line{2, Float32}, 1, Tuple{Point2f0, Point2f0}, TupleView{Tuple{Point2f0, Point2f0}, 2, 1, Vector{Point2f0}}}},
                                        Vector{LineString{2, Float32, Point2f0, Base.ReinterpretArray{Line{2, Float32}, 1, Tuple{Point2f0, Point2f0}, TupleView{Tuple{Point2f0, Point2f0}, 2, 1, Vector{Point2f0}}}}}},))
    end

    @assert precompile(split_intersections, (Vector{Point2f0},))
end


text:
module RubiSymbolics
using SpecialFunctions, Symbolics
using Symbolics: value, get_variables
using SymbolicUtils: Symbolic, Sym, @rule, @theory, Prewalk, Postwalk, PassThrough, Fixpoint, Chain
using Metatheory: @matchable, DynamicRule
using Metatheory.Syntax: rmlines
import TermInterface
import Elliptic, HypergeometricFunctions

@matchable struct Antiderivative
    expr::Symbolic
    x::Sym
    Antiderivative(expr::Num, x) = new(value(expr), x)
end

identifier(s::String) = Regex("(?<![a-zA-Z_])$s(?![a-zA-Z0-9_])")
identifier(sym::Symbol) = identifier(string(sym))

mathematica_version = 13    # just a dummy value, no real meaning yet

syntax_repl_dict = Dict(
    identifier("E") => "ℯ",
    identifier("Pi") => "pi",
    identifier("I") => "im",
    identifier("True") => "true",
    identifier("False") => "false",
    identifier("Null") => "nothing",
    "(*"    =>  "#=",
    "*)"    =>  "=#",
    "{"     =>  "[",
    "}"     =>  "]",
    "["     =>  "(",
    "]"     =>  ")",
    r"(?<=([a-zA-Z0-9_\)\]])|(_\.)) (?=[a-zA-Z0-9_\(\[])"=>  "*",   # space as *
    "x -("              =>  "x - (",                # space required inside [ ]
    r"(\d)(?=\.[^\d])"  => s"\1.0",                 # float fix
    "(a_. + b_.*x_)!"   => "factorial(~a' + ~b'*~x)",# fix in 8.6 Gamma Functions
    "(a + b*x)!"        => "factorial(a + b*x)",    # fix in 8.6 Gamma Functions
    # "\$VersionNumber"   => string(mathematica_version),  # remove $ character
    # r"(?<base>\(\h*-\h*\d+\)\h*\^\h*(?<exp>()|())" => s"\g<base>"
    "/;"    =>  "<--",   # temp. replacement for conditionals
    "/."    =>  "←",     # temp. replacement for replaceall
    "/:"    =>  "↶",     # temp. replacement for TagSetDelayed
    "->"    =>  "=>",
    "=!="   =>  "!==",
    "f_'''" =>  "Derivative(3)(f_)",    # Derivatives
    "f_''"  =>  "Derivative(2)(f_)",
    "f_'"   =>  "Derivative(1)(f_)",
    "g_'''" =>  "Derivative(3)(g_)",
    "g_''"  =>  "Derivative(2)(g_)",
    "g_'"   =>  "Derivative(1)(g_)",
    "_."    =>  "_'",    # "." --> "'"
    r"(?<expr>(#\d*))\[\[(?<idxs>\h*(\d+)|(\w+)\h*(,\h*\d+\h*)*)\]\]" => s"Part(Symbol(\"\g<expr>\"), \g<idxs>)",
    r"(?<expr>\w[a-zA-Z0-9]*)\[\[(?<idxs>\h*(\d+)|(\w+)\h*(,\h*\d+\h*)*)\]\]" => s"Part(\g<expr>, \g<idxs>)",
    r"(?<expr>(\w[a-zA-Z0-9]*)\[[^\[]*\])\[\[(?<idxs>\h*(\d+)|(\w+)\h*(,\h*\d+\h*)*)\]\]" => s"Part(\g<expr>, \g<idxs>)",
    r"#(\d)"=>  s"Symbol(\"#\1\")",
    "#"     =>  "Symbol(\"#1\")",
    "\$" => "DOLLAR",
)

function_repl_dict = Dict(
    identifier("Expand") => "expand",
    identifier("Int") => "Antiderivative",
    identifier("Function") => "MFunction",

    identifier("Exp") => "exp",
    identifier("Sqrt") => "sqrt",
    identifier("Log") => "log",

    identifier("Sin") => "sin",
    identifier("Sinh") => "sinh",
    identifier("ArcSin") => "asin",
    identifier("ArcSinh") => "asinh",
    identifier("Cos") => "cos",
    identifier("Cosh") => "cosh",
    identifier("ArcCos") => "acos",
    identifier("ArcCosh") => "acosh",
    identifier("Tan") => "tan",
    identifier("Tanh") => "tanh",
    identifier("ArcTan") => "atan",
    identifier("ArcTanh") => "atanh",
    identifier("Cot") => "cot",
    identifier("Coth") => "coth",
    identifier("ArcCot") => "acot",
    identifier("ArcCoth") => "acoth",
    identifier("Sec") => "sec",
    identifier("Sech") => "sech",
    identifier("ArcSec") => "asec",
    identifier("ArcSech") => "asech",
    identifier("Csc") => "csc",
    identifier("Csch") => "csch",
    identifier("ArcCsc") => "acsc",
    identifier("ArcCsch") => "acsch",

    identifier("EllipticK") => "Elliptic.K",
    identifier("EllipticF") => "Elliptic.F",
    identifier("EllipticE") => "Elliptic.E",
    identifier("EllipticPi") => "Elliptic.Pi",
    identifier("Hypergeometric2F1") => "HypergeometricFunctions._₂F₁",
    identifier("Erf") => "erf",
    identifier("Erfi") => "erfi",
    # # identifier("PolyLog") => "Polylogarithms.polylog",
)
@register_symbolic Base.exp(x::Complex{Real})::Complex false
@register_symbolic Base.sqrt(x::Complex{Real})::Complex false
@register_symbolic Base.log(x::Complex{Real})::Complex false
@register_symbolic Base.sin(x::Complex)::Complex false
@register_symbolic Base.sinh(x::Complex{Real})::Complex false
@register_symbolic Base.asin(x::Complex{Real})::Complex false
@register_symbolic Base.asinh(x::Complex{Real})::Complex false
@register_symbolic Base.cos(x::Complex)::Complex false
@register_symbolic Base.cosh(x::Complex{Real})::Complex false
@register_symbolic Base.acos(x::Complex{Real})::Complex false
@register_symbolic Base.acosh(x::Complex{Real})::Complex false
@register_symbolic Base.tan(x::Complex{Real})::Complex false
@register_symbolic Base.tanh(x::Complex)::Complex false
@register_symbolic Base.atan(x::Complex{Real})::Complex false
@register_symbolic Base.atanh(x::Complex)::Complex false
@register_symbolic Base.cot(x::Complex)::Complex false
@register_symbolic Base.coth(x::Complex)::Complex false
@register_symbolic Base.acot(x::Complex)::Complex false
@register_symbolic Base.acoth(x::Complex)::Complex false
@register_symbolic Base.sec(x::Complex)::Complex false
@register_symbolic Base.sech(x::Complex)::Complex false
@register_symbolic Base.asec(x::Complex)::Complex false
@register_symbolic Base.asech(x::Complex)::Complex false
@register_symbolic Base.csc(x::Complex)::Complex false
@register_symbolic Base.csch(x::Complex)::Complex false
@register_symbolic Base.acsc(x::Complex)::Complex false
@register_symbolic Base.acsch(x::Complex)::Complex false
@register_symbolic Elliptic.K(m)
@register_symbolic Elliptic.K(m::Complex)::Complex false
@register_symbolic Elliptic.F(phi, m)
@register_symbolic Elliptic.F(phi::Complex, m::Complex)::Complex false
@register_symbolic Elliptic.E(phi, m)
# @register_symbolic Elliptic.E(phi, m)::Complex
@register_symbolic Elliptic.E(m) false
@register_symbolic Elliptic.E(m::Complex)::Complex false
@register_symbolic Elliptic.Pi(n, phi, m)
# @register_symbolic Elliptic.Pi(n, phi, m)::Complex
@register_symbolic HypergeometricFunctions._₂F₁(a, b, c, z)
@register_symbolic HypergeometricFunctions._₂F₁(a, b, c, z::Complex)::Complex false
@register_symbolic HypergeometricFunctions._₂F₁(a, b::Complex, c::Complex, z::Complex)::Complex false

integrate(exp, x) = x   # TODO:fix dummy
If(cond, true_stmt, false_stmt=nothing) = cond ? true_stmt : false_stmt
LogIntegral(z) = nothing
SinIntegral(z) = nothing
SinhIntegral(z) = nothing
CosIntegral(z) = nothing
CoshIntegral(z) = nothing
ExpIntegralE(n, z) = nothing
ExpIntegralEi(z) = nothing
FresnelS(z) = nothing
FresnelC(z) = nothing
PolyLog(n, z) = nothing
ProductLog(z) = nothing
Zeta(z) = nothing
Gamma(z) = nothing
Gamma(a, z) = nothing
Gamma(a, z0, z1) = nothing
AppellF1(a, b1, b2, c, x, y) = nothing

Unintegrable(ex) = nothing
Unintegrable(x, y) = nothing
CannotIntegrate(ex) = nothing
CannotIntegrate(x, y) = nothing

function D end
function Integrate end
function Sum end
function Product end
function Dif end

@register_symbolic LogIntegral(z)
@register_symbolic SinIntegral(z)
@register_symbolic SinhIntegral(z)
@register_symbolic CosIntegral(z)
@register_symbolic CoshIntegral(z)
@register_symbolic ExpIntegralE(n, z)
@register_symbolic ExpIntegralEi(z)
@register_symbolic FresnelS(z)
@register_symbolic FresnelC(z)
@register_symbolic PolyLog(n, z)
@register_symbolic PolyLog(n, z::Complex)::Complex false
@register_symbolic ProductLog(z)
@register_symbolic Zeta(z)
@register_symbolic Gamma(z)
@register_symbolic Gamma(a, z) false
@register_symbolic Gamma(a, z::Complex)::Complex false
@register_symbolic Gamma(a, z0, z1) false
@register_symbolic AppellF1(a, b1, b2, c, x, y)
@register_symbolic AppellF1(a, b1, b2, c, x::Complex, y::Complex)::Complex false
# @register_symbolic AppellF1(a, b1, b2, c, x, y::Complex)::Complex false
# @register_symbolic AppellF1(a, b1, b2, c, x::Complex, y)::Complex false
@register_symbolic Unintegrable(ex)
@register_symbolic Unintegrable(x, y) false
@register_symbolic CannotIntegrate(ex)
@register_symbolic CannotIntegrate(x, y) false

DOLLARTimeLimit = nothing

NumberQ(::Any) = false
NumberQ(::Number) = true

Not(b::Bool) = !b

function FreeQ(expr, x::Sym)
    vars = get_variables(expr)
    return !any(isequal.(x, vars))
end
FreeQ(expr, x::Num) = FreeQ(expr, value(x))
FreeQ(expr_lst::AbstractArray, x) = all(FreeQ.(expr_lst, x))

Refine(expr) = simplify(expr)   # no assumptionsystem -> should be OK, speed ?
Quiet(expr) = expr  # nothing to mute here ?!
PossibleZeroQ(x) = false    # Would check numerically for zero -> use simplify straight away

abstract type AbstractAttribute end
struct HoldFirstType <: AbstractAttribute end
HoldFirst = HoldFirstType()
struct HoldAllType <: AbstractAttribute end
HoldAll = HoldAllType()
SetAttributes(obj, a::AbstractAttribute) = nothing
function Hold end
function HoldForm end
function Defer end
function Pattern end
Clear(x) = nothing
ClearAll(args...) = nothing
MemberQ(lst, form) = nothing

const util_rules = DynamicRule[]
expand_utils = Prewalk(Chain(util_rules))

cond_repl = Prewalk(PassThrough(@rule (~a <-- ~b) --> (expand_utils(~a) where expand_utils(~b))))
macro apply_utils(expr::Expr)
    expr = cond_repl(expr)  # here because printing with 'where' fails
    return :($(esc(expr)))
end

rm_block(expr) = expr
function rm_block(expr::Expr)
    if expr.head == :block && length(expr.args) == 2 && expr.args[1] isa LineNumberNode
        return expr.args[2]
    else
        return Expr(expr.head, map(rm_block, expr.args)...)
    end
end

macro rm_blocks(expr::Expr)
    return esc(rm_block(expr))
end

r_if1 = @rule @rm_blocks If(~cond, ~then_var = ~then_val) --> (~cond ? ~then_var = ~then_val : nothing)
r_if2 = @rule If(~cond, ~then) --> (~cond ? ~then : nothing)
r_if3 = @rule If(~cond, ~then, ~else_stmt) --> (~cond ? ~then : ~else_stmt)
if_simplifier = Prewalk(Chain([r_if1, r_if2, r_if3]))

r_anon_sym = @theory begin
    Symbol("#") --> :__args[1]
    Symbol("#1") --> :__args[1]
    Symbol("#2") --> :__args[2]
    Symbol("#3") --> :__args[3]
    Symbol("#4") --> :__args[4]
    Symbol("#5") --> :__args[5]
    Symbol("#6") --> :__args[6]
    Symbol("#7") --> :__args[7]
    Symbol("#8") --> :__args[8]
    Symbol("#9") --> :__args[9]
    Symbol("#10") --> :__args[10]
end
convert_ht_syms = Prewalk(Chain(r_anon_sym))
r_lambda = @rule MFunction(~body) => :((__args...) -> $(convert_ht_syms(~body)))
lambda_simplifier = Prewalk(PassThrough(r_lambda))
r_multi_stmt = @rule Block(~lst, ~expr1; ~expr2) --> Block(~lst, begin ~expr1; ~expr2 end)
multi_stmt_simplifier = Prewalk(PassThrough(r_multi_stmt))

r_with = @rule With([~assignments...], ~expr) => begin
    rlst = []
    for ass in ~assignments
        @assert ass.head == :(=)
        r = eval(:(@rule $(Meta.quot(:($(ass.args[1])))) --> $(Meta.quot(:($(ass.args[2]))))))
        push!(rlst, r)
    end
    return Prewalk(Chain(rlst))(~expr)
end
with_simplifier = Prewalk(PassThrough(r_with))

r_module = @rule (Module([~assignments...], ~expr)) => Expr(:let, Expr(:block, ~assignments...), ~expr)
module_simplifier = Prewalk(PassThrough(r_module))

macro util_collect(expr::Expr)
    e = rmlines(expr)
    e = macroexpand(__module__, e)
    e = if_simplifier(e)
    e = lambda_simplifier(e)
    e = multi_stmt_simplifier(e)
    e = module_simplifier(e)
    e = with_simplifier(e)
    if e.head == :block
        rule_defs = filter(x-> x.head == :(:=) , e.args)
        declarations = Expr[]
        for stmt in rule_defs
            @assert stmt.args[1].head == :call
            name = stmt.args[1].args[1]
            args = stmt.args[1].args
            num_args = length(args) - 1
            decl = quote
                $name(args::Real...) = SymbolicUtils.Term{Real}($name, args)
                @register_symbolic $name($(Tuple(:x for n in 1:num_args)...)) false
                push!(util_rules, @rule $(stmt.args[1]) => $(stmt.args[2]))
            end
            pushfirst!(declarations, decl)
        end
        prepend!(e.args, declarations)
        filter!(x -> x.head != :(:=), e.args)
        return esc(e)
    else
        error("theory is not in form begin a => b; ... end")
    end
end

export  integrate, Antiderivative
export  If, LogIntegral, SinIntegral, SinhIntegral, CosIntegral,
        CoshIntegral, ExpIntegralE, ExpIntegralEi, FresnelS, FresnelC, PolyLog,
        ProductLog, Zeta, Gamma, AppellF1, Unintegrable, CannotIntegrate
export erf, erfi
export FreeQ

end # module


text:
"""
    estimation_input_test(estim::VMASettings, benchmark_data::Tuple)

Check whether the settings in `estim` are equal to those in `benchmark_data`.

    estimation_input_test(estim_a::VMASettings, estim_b::VMASettings)

Check whether the settings in `estim_a` are equal to those in `estim_b`.
"""
function estimation_input_test(estim::VMASettings, benchmark_data::Tuple)

    benchmark_Y, benchmark_n, benchmark_T, benchmark_r, benchmark_nr, benchmark_m, benchmark_λ,
        benchmark_α, benchmark_β, benchmark_Γ, benchmark_ε, benchmark_tol, benchmark_max_iter,
            benchmark_prerun, benchmark_verb = benchmark_data;

    return ~(false in [estim.Y === benchmark_Y;
                      estim.n == benchmark_n;
                      estim.T == benchmark_T;
                      estim.r == benchmark_r;
                      estim.nr == benchmark_nr;
                      estim.m == benchmark_m;
                      estim.λ == benchmark_λ;
                      estim.α == benchmark_α;
                      estim.β == benchmark_β;
                      estim.Γ == benchmark_Γ;
                      estim.ε == benchmark_ε;
                      estim.tol == benchmark_tol;
                      estim.max_iter == benchmark_max_iter;
                      estim.prerun == benchmark_prerun;
                      estim.verb == benchmark_verb]);
end

function estimation_input_test(estim_a::VMASettings, estim_b::VMASettings)

    return ~(false in [estim_a.Y === estim_b.Y;
                      estim_a.n == estim_b.n;
                      estim_a.T == estim_b.T;
                      estim_a.r == estim_b.r;
                      estim_a.nr == estim_b.nr;
                      estim_a.m == estim_b.m;
                      estim_a.λ == estim_b.λ;
                      estim_a.α == estim_b.α;
                      estim_a.β == estim_b.β;
                      estim_a.Γ == estim_b.Γ;
                      estim_a.ε == estim_b.ε;
                      estim_a.tol == estim_b.tol;
                      estim_a.max_iter == estim_b.max_iter;
                      estim_a.prerun == estim_b.prerun;
                      estim_a.verb == estim_b.verb]);
end

"""
    vma_estimation_test(Y::JArray, r::Int64, λ::Number, α::Number, β::Number, benchmark_data::Tuple)

Run basic tests on VMASettings.
"""
function vma_estimation_test(Y::JArray, r::Int64, λ::Number, α::Number, β::Number, benchmark_data::Tuple)

    #=
    Tests on VMASettings
    =#

    # First group of settings

    estim1 = VMASettings(Y, r, λ, α, β);
    estim2 = VMASettings(Y, r, λ, α, β, ε=1e-4);
    estim3 = VMASettings(Y, r, λ, α, β, ε=1e-4, tol=1e-4);
    estim4 = VMASettings(Y, r, λ, α, β, ε=1e-4, tol=1e-4, max_iter=1000);
    estim5 = VMASettings(Y, r, λ, α, β, ε=1e-4, tol=1e-4, max_iter=1000, prerun=2);
    estim6 = VMASettings(Y, r, λ, α, β, ε=1e-4, tol=1e-4, max_iter=1000, prerun=2, verb=true);

    @test estimation_input_test(estim1, benchmark_data);
    @test estimation_input_test(estim1, estim2);
    @test estimation_input_test(estim1, estim3);
    @test estimation_input_test(estim1, estim4);
    @test estimation_input_test(estim1, estim5);
    @test estimation_input_test(estim1, estim6);

    # Second group of settings

    estim1b = VMASettings(Y, r, λ, α, β, ε=1e-8);
    estim2b = VMASettings(Y, r, λ, α, β, ε=1e-8, tol=1e-4);
    estim3b = VMASettings(Y, r, λ, α, β, ε=1e-8, tol=1e-4, max_iter=1000);
    estim4b = VMASettings(Y, r, λ, α, β, ε=1e-8, tol=1e-4, max_iter=1000, prerun=2);
    estim5b = VMASettings(Y, r, λ, α, β, ε=1e-8, tol=1e-4, max_iter=1000, prerun=2, verb=true);

    @test estimation_input_test(estim1b, estim2b);
    @test estimation_input_test(estim1b, estim3b);
    @test estimation_input_test(estim1b, estim4b);
    @test estimation_input_test(estim1b, estim5b);

    # Third group of settings

    estim1c = VMASettings(Y, r, λ, α, β, ε=1e-8, tol=1e-16);
    estim2c = VMASettings(Y, r, λ, α, β, ε=1e-8, tol=1e-16, max_iter=1000);
    estim3c = VMASettings(Y, r, λ, α, β, ε=1e-8, tol=1e-16, max_iter=1000, prerun=2);
    estim4c = VMASettings(Y, r, λ, α, β, ε=1e-8, tol=1e-16, max_iter=1000, prerun=2, verb=true);

    @test estimation_input_test(estim1c, estim2c);
    @test estimation_input_test(estim1c, estim3c);
    @test estimation_input_test(estim1c, estim4c);

    # Fourth group of settings

    estim1d = VMASettings(Y, r, λ, α, β, ε=1e-8, tol=1e-16, max_iter=10000);
    estim2d = VMASettings(Y, r, λ, α, β, ε=1e-8, tol=1e-16, max_iter=10000, prerun=2);
    estim3d = VMASettings(Y, r, λ, α, β, ε=1e-8, tol=1e-16, max_iter=10000, prerun=2, verb=true);

    @test estimation_input_test(estim1d, estim2d);
    @test estimation_input_test(estim1d, estim3d);

    # Fifth group of settings

    estim1e = VMASettings(Y, r, λ, α, β, ε=1e-8, tol=1e-16, max_iter=10000, prerun=3);
    estim2e = VMASettings(Y, r, λ, α, β, ε=1e-8, tol=1e-16, max_iter=10000, prerun=3, verb=true);

    @test estimation_input_test(estim1e, estim2e);

    # Sixth group of settings

    estim1f = VMASettings(Y, r, λ, α, β, ε=1e-8, tol=1e-16, max_iter=10000, prerun=3, verb=false);

    @test estim1f.verb == false;

    #=
    Set default settings
    =#

    estim = estim1;

    #=
    Tests on initialisation
    =#

    sspace, B_star, C_star, Q_view, coordinates_measurement_states, coordinates_transition_current, coordinates_transition_lagged, coordinates_transition_PPs, coordinates_transition_P0, coordinates_free_params_B, coordinates_free_params_C = MessyTimeSeriesOptim.initialise(estim);

    @test sspace.B[:, 1:estim.n] == Matrix(1.0I, estim.n, estim.n);
    @test sspace.C[1:estim.n, :] == zeros(estim.n, estim.n+estim.m);

    #=
    Tests on ksmoother_ecm!
    =#

    # Run Kalman filter
    status = SizedKalmanStatus(sspace.Y.T);
    kfilter_full_sample!(sspace, status);

    # Compute Kalman stats from the ksmoother function
    X_sm, P_sm, X0_sm, P0_sm = ksmoother(sspace, status);

    # Memory pre-allocation
    smoother_arrays_benchmark = MessyTimeSeriesOptim.SmootherArrays(estim, sspace, coordinates_measurement_states, coordinates_transition_current, coordinates_transition_lagged, coordinates_transition_PPs);

    # Compute F_benchmark, G_benchmark and H_benchmark
    for t=length(X_sm)-1:-1:1
        copyto!(smoother_arrays_benchmark.Xs_leading, X_sm[t+1]);
        copyto!(smoother_arrays_benchmark.Ps_leading, P_sm[t+1]);
        Xs = X_sm[t];
        Ps = P_sm[t];
        MessyTimeSeriesOptim.update_ecm_stats_transition!(smoother_arrays_benchmark.F, smoother_arrays_benchmark.G, estim, smoother_arrays_benchmark, Xs, Ps, coordinates_transition_current, coordinates_transition_lagged, coordinates_transition_PPs);
    end

    copyto!(smoother_arrays_benchmark.Xs_leading, X_sm[1]);
    copyto!(smoother_arrays_benchmark.Ps_leading, P_sm[1]);
    MessyTimeSeriesOptim.update_ecm_stats_transition!(smoother_arrays_benchmark.F, smoother_arrays_benchmark.G, estim, smoother_arrays_benchmark, X0_sm, P0_sm, coordinates_transition_current, coordinates_transition_lagged, coordinates_transition_PPs);

    # Compute Kalman stats from the ksmoother_ecm! function
    smoother_arrays = MessyTimeSeriesOptim.SmootherArrays(estim, sspace, coordinates_measurement_states, coordinates_transition_current, coordinates_transition_lagged, coordinates_transition_PPs);
    MessyTimeSeriesOptim.ksmoother_ecm!(estim, sspace, status, smoother_arrays, coordinates_measurement_states, coordinates_transition_current, coordinates_transition_lagged, coordinates_transition_PPs, coordinates_transition_P0);

    # Tests
    @test smoother_arrays.F == smoother_arrays_benchmark.F;
    @test smoother_arrays.G == smoother_arrays_benchmark.G;
    @test smoother_arrays.H == smoother_arrays_benchmark.H;
    @test sspace.X0 == X0_sm;
    @test sspace.P0 == P0_sm;
end

"""
    naive_vma_simulation(true_coefficient::FloatMatrix, n::Int64, T::Int64, r::Int64, λ::Number, α::Number, β::Number; no_sim::Int64=500, verb::Bool=false)

Run a VMA(r) naive simulation.
"""
function naive_vma_simulation(true_coefficient::FloatMatrix, n::Int64, T::Int64, r::Int64, λ::Number, α::Number, β::Number; no_sim::Int64=500, verb::Bool=false)

    # Run simulation
    rng = StableRNG(1);
    B_sim = zeros(n, n*r, no_sim);
    params_accuracy = zeros(no_sim);
    oos_fc_error = zeros(no_sim);

    for sim=1:no_sim

        # Simulate data with mean zero and unit standard deviation
        v = randn(rng, n, T+r);
        Y = zeros(n, T+r);
        for t=1:T+r
            if t <= r
                Y[:,t] .= 0.0;
            else
                Y[:,t] = true_coefficient*vcat([v[:,t-lag][:] for lag=1:r]...) + v[:,t];
            end
        end
        Y = standardise(Y[:, r+1:end]);

        # Estimate
        t0 = floor(T/2) |> Int64;
        estim = VMASettings(Y[:,1:t0], r, λ, α, β, verb=false);
        sspace = ecm(estim, output_sspace_data=Y);

        # Run Kalman filter
        status = kfilter_full_sample(sspace);

        # Forecast
        forecast = mapreduce(Xt -> sspace.B*Xt, hcat, status.history_X_prior);

        # Store estimated parameters
        B_sim[:,:,sim] = sspace.B[1:n, n+1:n+n*r];

        # Store params_accuracy
        difference_B_truth = @views B_sim[:,:,sim] - true_coefficient;
        params_accuracy[sim] = sqrt(tr(difference_B_truth'*difference_B_truth));

        # Store oos_fc_error
        oos_fc_error[sim] = @views mean(mean(abs.(Y[:,t0+1:end] - forecast[:,t0+1:end]), dims=2)./mean(abs.(Y[:,t0+1:end]), dims=2));
    end

    if verb == true
        println("=============================================================================================================================================")
        println("True parameters:")
        println(true_coefficient);
        println("=============================================================================================================================================")
        println("Estimated parameters with λ=$(λ), α=$(α), and β=$(β):")
        println(round.(median(B_sim, dims=3)[:,:], digits=5));
        println("=============================================================================================================================================")
        println("");
    end

    return B_sim, params_accuracy, oos_fc_error;
end

"""
    vma_simulation_tests(model_id::String, λ::Number, α::Number, β::Number)

Default tests based on VMA simulations.
"""
function vma_simulation_tests(model_id::String, λ::Number, α::Number, β::Number)

    # Run simulations
    _, bivariate_vma1_full, _ = naive_vma_simulation([0.3 -0.2; -0.2 0.5], 2, 200, 1, λ, α, β, no_sim=10);
    _, bivariate_vma2_full, _ = naive_vma_simulation([0.5 0.4 0.2 0.1; 0.3 -0.2 0.1 0.1], 2, 200, 2, λ, α, β, no_sim=10);
    _, bivariate_vma2_sparse, _ = naive_vma_simulation([0.3 0.0 0.0 0.0; 0.0 -0.1 0.0 0.0], 2, 200, 2, λ, α, β, no_sim=10);
    _, trivariate_vma2_groups, _ = naive_vma_simulation([0.5 0.4 0.0 0.2 0.1 0.0; 0.3 -0.2 0.0 0.1 0.1 0.0; 0.0 0.0 0.9 0.0 0.0 0.0], 3, 200, 2, λ, α, β, no_sim=10);

    # Load benchmarks
    benchmark_bivariate_vma1_full = read_test_input("./input/vma/$(model_id)_bivariate_vma1_full");
    benchmark_bivariate_vma2_full = read_test_input("./input/vma/$(model_id)_bivariate_vma2_full");
    benchmark_bivariate_vma2_sparse = read_test_input("./input/vma/$(model_id)_bivariate_vma2_sparse");
    benchmark_trivariate_vma2_groups = read_test_input("./input/vma/$(model_id)_trivariate_vma2_groups");

    # Perform tests
    @test norm(bivariate_vma1_full-benchmark_bivariate_vma1_full, Inf) <= 1e-8;
    @test norm(bivariate_vma2_full-benchmark_bivariate_vma2_full, Inf) <= 1e-8;
    @test norm(bivariate_vma2_sparse-benchmark_bivariate_vma2_sparse, Inf) <= 1e-8;
    @test norm(trivariate_vma2_groups-benchmark_trivariate_vma2_groups, Inf) <= 1e-8;
end

@testset "VMA basic functionalities" begin

    # Initialise data and state-space parameters
    Y = [0.72 missing 1.86 missing missing 2.52 2.98 3.81 missing 4.36;
         0.95 0.70 missing missing missing missing 2.84 3.88 3.84 4.63];

    r = 2;
    nr = 4;
    m = 6;
    λ = 0.0;
    α = 0.5;
    β = 1.0;
    ε = 1e-4;
    tol = 1e-4;
    max_iter = 1000;
    prerun = 2;
    verb = true;

    Γ_1 = build_Γ(2, r, λ, β);
    Γ_2 = build_Γ(2, r, λ+1, β);
    Γ_3 = build_Γ(2, r, λ+1, β+1);

    # Benchmark data
    benchmark_data_1 = (Y, size(Y,1), size(Y,2), r, nr, m, λ, α, β, Γ_1, ε, tol, max_iter, prerun, verb);
    benchmark_data_2 = (Y, size(Y,1), size(Y,2), r, nr, m, λ+1, α, β, Γ_2, ε, tol, max_iter, prerun, verb);
    benchmark_data_3 = (Y, size(Y,1), size(Y,2), r, nr, m, λ+1, α, β+1, Γ_3, ε, tol, max_iter, prerun, verb);

    # Run tests
    vma_estimation_test(Y, r, λ, α, β, benchmark_data_1);
    vma_estimation_test(Y, r, λ+1, α, β, benchmark_data_2);
    vma_estimation_test(Y, r, λ+1, α, β+1, benchmark_data_3);
end

@testset "VMA simulations: MLE" begin
    vma_simulation_tests("mle", 0.0, 0.0, 1.0);
end

@testset "VMA simulations: penalised MLE" begin

    models_triplets = Dict("ridge" => (10.0, 0.0, 1.0), "lasso" => (10.0, 1.0, 1.0), "en" => (10.0, 0.5, 1.0));

    for (key, value) in models_triplets
        vma_simulation_tests(key, value...);
    end
end


text:
#!/usr/bin/env julia

# Uses the EnsembleSampler to explore the posterior.
# In contrast to codes like mach_three.jl, this architecture means that within each likelihood call, the global lnprob for all channels is evaluated in *serial*, while the walkers themselves are parallelized.

using ArgParse
using Distributed
@everywhere using Printf
@everywhere using Test

s = ArgParseSettings()
@add_arg_table s begin
    "--p", "-p"
    help = "number of processes to add"
    arg_type = Int
    default = 0
    "--run_index", "-r"
    help = "Output run index"
    arg_type = Int
    "--config"
    help = "a YAML configuration file"
    default = "config.yaml"
    "--cpus"
    help = "Which CPUS to add"
    arg_type = Array{Int,1}
    eval_arg = true
    "--optim"
    help = "Optimize instead of sample?"
    action = :store_true
    "--test"
    help = "Is this a test run of venus.jl? Allow using many fewer walkers for eval purposes."
    action = :store_true
    "--MPI"
    help = "Run the script using MPI"
    action = :store_true
end

parsed_args = parse_args(ARGS, s)

cpus = parsed_args["cpus"]

if cpus != nothing
    using ClusterManagers
    np = length(cpus)
    addprocs(LocalAffinityManager(;np = np, affinities = cpus))
end

# Since we've made this a #!/usr/bin/env julia script, we can no longer specify the extra
# processors via the julia executable, so we need to ghost this behavior back in.
if parsed_args["p"] > 0
    addprocs(parsed_args["p"])
end

# Use MPI as the transport mechanism--good for a cluster environment
if parsed_args["MPI"]
    import MPI
    MPI.Init()
    rank = MPI.Comm_rank(MPI.COMM_WORLD)
    size_workers = MPI.Comm_size(MPI.COMM_WORLD)

    manager = MPI.start_main_loop(MPI.MPI_TRANSPORT_ALL)
end

import YAML
config = YAML.load(open(parsed_args["config"]))

outfmt(run_index::Int) = config["out_base"] * @sprintf("run%02d/", run_index)

# This code is necessary for multiple simultaneous runs on a high performance cluster
# so that different runs do not write into the same output directory
if parsed_args["run_index"] == nothing
    global run_index = 0
    global outdir = outfmt(run_index)
    while ispath(outdir)
        println(outdir, " exists")
        global run_index += 1
        global outdir = outfmt(run_index)
    end
    # are we starting in a fresh directory?
    fresh = true
else
    global run_index = parsed_args["run_index"]
    global outdir = outfmt(run_index)
    if ispath(outdir)
        # we are not starting in a fresh directory, so try loading pos0.npy from this directory.
        fresh = false
    else
        fresh = true
    end
end

if fresh
    # make the output directory
    println("Creating ", outdir)
    mkdir(outdir)
    # Copy the config.yaml file here for future reference
    cp(parsed_args["config"], outdir * parsed_args["config"])
else
    println("$outdir exists, using current.")
end

@everywhere using AbstractFFTs
@everywhere using DiskJockey.constants
@everywhere using DiskJockey.visibilities
@everywhere using DiskJockey.image
@everywhere using DiskJockey.gridding
@everywhere using DiskJockey.model

# Determine if we will be including the User-defined prior
if isfile("prior.jl")
    println("Including user-defined prior from prior.jl")
    @everywhere include(pwd() * "/prior.jl")
    # Make a copy to the outdirectory for future reference
    cp("prior.jl", outdir * "prior.jl")
end

# load data and figure out how many channels
dvarr = DataVis(config["data_file"])
nchan = length(dvarr)

lam0 = lam0s[config["species"] * config["transition"]]
# calculate the velocities corresponding to dvarr
lams = Float64[dv.lam for dv in dvarr]
vels = c_kms * (lams .- lam0) / lam0

# create a mask that excludes channels we don't want
if haskey(config, "exclude")
    exclude = config["exclude"]
    vel_mask = generate_vel_mask(exclude, vels)
else
    vel_mask = trues(nchan)
end

println("Using the following channels (velocities in km/s): ", vels[vel_mask])

# This program is meant to be started with the -p option.
nchild = nworkers()
println("Workers allocated ", nchild)

# make the values of run_index and config available on all processes
for process in procs()
    @spawnat process global run_id = run_index
    @spawnat process global cfg = config
    @spawnat process global vmask = vel_mask
end

# Convert the parameter values from the config.yaml file from Dict{String, Float64} to
# Dict{Symbol, Float64} so that they may be splatted as extra args into the convert_vector method
@everywhere xargs = Dict{Symbol,Float64}(Symbol(index) => value for (index, value) in pairs(cfg["parameters"]))

# Read the fixed and free parameters from the config.yaml file
@everywhere fix_params = cfg["fix_params"]

# Set the correct model, and the parameters that will be fixed
@everywhere function convert_p(p::Vector{Float64})
    return convert_vector(p, cfg["model"], fix_params; xargs...)
end

# Now, redo this to only load the dvarr for the keys that we need, and conjugate
# @everywhere dvarr = DataVis(cfg["data_file"], kl)
@everywhere dvarr = DataVis(cfg["data_file"], vmask)
@everywhere visibilities.conj!(dvarr)
@everywhere nchan = length(dvarr)

@everywhere const global species = cfg["species"]

@everywhere const global npix = cfg["npix"] # number of pixels, can alternatively specify x and y separately

# Keep track of the current home directory
@everywhere const global homedir = pwd() * "/"

# Create the model grid
@everywhere grd = cfg["grid"]
@everywhere grid = Grid(grd)

# calculate the interpolation closures, since we are keeping the angular size of the image
# fixed thoughout the entire simulation.
# calculate dl and dm (assuming they are equal).
@everywhere dl = cfg["size_arcsec"] * arcsec / cfg["npix"] # [radians/pixel]

# Determine dRA and dDEC from the image and distance
# These will be used to correct for the half-pixel offset
@everywhere half_pix = dl / (arcsec * 2.) # [arcsec/half-pixel]

@everywhere uu = fftshift(fftfreq(npix, dl)) * 1e-3 # [kλ]
@everywhere vv = fftshift(fftfreq(npix, dl)) * 1e-3 # [kλ]

# For each channel, calculate the interpolation closures
@everywhere int_arr = Array{Function}(undef, nchan)
@everywhere for (i, dset) in enumerate(dvarr)
    int_arr[i] = plan_interpolate(dset, uu, vv)
end

# This function is fed to the EnsembleSampler
# That means, using the currently available global processes, like the data visibilities,
# and it must create it's own temporary directory to write the necessary files for
# RADMC to run.
@everywhere function fprob(p::Vector{Float64})
    # Each walker needs to create it's own temporary directory
    # where all RADMC-3D files will reside and be driven from
    # It only needs to last for the duration of this function, so we use a tempdir
    keydir = mktempdir() * "/"

    lnp = try
        # Convert the vector to a specific type of pars, e.g. ParametersStandard, ParametersTruncated, etc, which will be used for multiple dispatch from here on out.
        pars = convert_p(p) # will raise an error if can't convert

        lnpr = lnprior(pars, grid)
        (sizeau_desired, sizeau_command) = size_au(cfg["size_arcsec"], pars.dpc, grid) # [AU]

        # Copy all relevant configuration scripts to the keydir so that RADMC-3D can run.
        # these are mainly setup files that will be static throughout the run
        # they were written by DJInitialize.jl and write_grid()
        for fname in ["radmc3d.inp", "wavelength_micron.inp", "lines.inp", "molecule_" * molnames[species] * ".inp"]
            run(`cp $(homedir)$fname $keydir`)
        end

        # change the subprocess to reside in this directory for the remainder of the run
        # where it will drive its own independent RADMC3D process for a subset of channels
        cd(keydir)

        # write the amr_grid.inp file
        write_grid(keydir, grid)

        # Compute disk structure files using model.jl, write to disk in current directory
        # Based upon typeof(pars), the subroutines within this function will dispatch to the correct model.
        write_model(pars, keydir, grid, species)

        # Doppler shift the dataset wavelengths to rest-frame wavelength
        beta = pars.vel / c_kms # relativistic Doppler formula
        lams = Array{Float64}(undef, nchan)
        for i = 1:nchan
            lams[i] =  dvarr[i].lam * sqrt((1. - beta) / (1. + beta)) # [microns]
        end

        write_lambda(lams, keydir) # write into current directory

        # Run RADMC-3D, redirect output to /dev/null
        run(pipeline(`radmc3d image incl $(pars.incl) posang $(pars.PA) npix $npix loadlambda sizeau $sizeau_command`, devnull))

        # Read the RADMC-3D images from disk
        img = try
            imread() # we should already be in the sub-directory, no path required
        # If the synthesized image is screwed up, just say there is zero probability.
        catch y
            if isa(y, SystemError)
                throw(ImageException("Failed to read synthesized image for parameters $pars"))
            else
                println("Unexpected error for $pars, exiting")
                throw(y)
            end
        end

        phys_size = img.pixsize_x * npix / AU

        # Convert raw images to the appropriate distance
        skim = imToSky(img, pars.dpc)

        # Apply the gridding correction function before doing the FFT
        # No shift needed, since we will shift the resampled visibilities
        corrfun!(skim)

        lnprobs = Array{Float64}(undef, nchan)
        # Do the Fourier domain stuff per channel
        for i = 1:nchan
            dv = dvarr[i]
            # FFT the appropriate image channel
            vis_fft = transform(skim, i)

            # Interpolate the `vis_fft` to the same locations as the DataSet
            mvis = int_arr[i](dv, vis_fft)

            # Apply the phase shift here, which includes a correction for the half-pixel offset
            # inherent to any image synthesized with RADMC3D.
            phase_shift!(mvis, pars.mu_RA + half_pix, pars.mu_DEC - half_pix)

            # Calculate the likelihood between these two using our chi^2 function.
            lnprobs[i] = lnprob(dv, mvis)
        end

        # Sum them all together
        # this quantity is assigned to the value of lnp at the start of the try/catch block
        sum(lnprobs) + lnpr

    catch y
        # If an error occured through some of the known routines in this package handle it and return -Inf
        if isa(y, DiskJockeyException)
            # print the error message and move on, returning -Inf
            println(y)
        else
            # if the error is not coming from one of the channels we
            # designed, actually terminate the program after printing the
            # parameter values at which the error occured
            println("Unforeseen error with parameter vector $p")
            println("Unforeseen error with parameter type $pars")
            throw(y)
        end

        # -Inf is returned as the value of lnp
        -Inf
    finally
        # regardless of whether we ignored the error or allowed it
        # to terminate the program, change back to the home directory and
        # then remove the temporary directory
        cd(homedir)
        run(`rm -rf $keydir`)
    end

    return lnp

end

# Decide if we're going to optimize the parameters or sample the posteriors.
if parsed_args["optim"]
    sparams = registered_params[cfg["model"]]

    # fit_params are the ones in sample_params that are not in fix_params
    fit_params = filter(x->!in(x, fix_params), sparams)

    # note the Float64 is how we prepend a type to the comprehension
    # p0 = Float64[cfg["parameters"][parname] for parname in fit_params]

    parrangesdict = cfg["parameter_ranges"]
    MaxFuncEvals = cfg["MaxFuncEvals"]

    # get the ranges
    sranges = [tuple(convert(Array{Float64,1}, cfg["parameter_ranges"][par])...) for par in fit_params]

    # get the number of parameters, starting pos, and bounding ranges
    # bounding ranges can be a dictionary in YAML
    nll = x->-fprob(x)
    @everywhere using BlackBoxOptim
    opt = bbsetup(nll; Method = :xnes, SearchRange = sranges, MaxFuncEvals = MaxFuncEvals, Workers = workers())
    res = bboptimize(opt)

    println(res)

else
    # instead, sample using MCMC
    using NPZ

    if fresh
        pos0 = npzread(config["pos0"])
    else
        pos0 = npzread(outdir * config["pos0"])
    end

    ndim, nwalkers = size(pos0)

    # Use the EnsembleSampler to do the optimization
    using DiskJockey.EnsembleSampler

    sampler = Sampler(nwalkers, ndim, fprob, parsed_args["test"])

    run_schedule(sampler, pos0, config["samples"], config["loops"], outdir)

    write_samples(sampler, outdir)

    if parsed_args["MPI"]
        MPI.stop_main_loop(manager)
    end
end


text:
module TestCore
using PhaseSpaceIO
using PhaseSpaceIO.Testing: arbitrary
using Test
using Base.Meta

@testset "show" begin
    p_iaea = IAEAParticle(typ=photon,E=1.23,x=1,y=2,z=3,
                 u=0,v=0,w=1)
    p_egs = EGSParticle(latch=Latch(charge=1), E=1.23, x=1, y=2,
                        u=0,v=0,w=1)
    for p in [p_egs, p_iaea,]
        io = IOBuffer()
        s = sprint(show, p)
        for prop in propertynames(p)
            sprop = string(prop)
            @test occursin(sprop, s)
        end

        ex = Meta.parse(s)
        @test eval(ex) === p
    end
end

end#module


text:
struct RegularGrid{T <: AbstractFloat} <: ApproximationGrid
    a::T
    b::T
    n::Int
end

function RegularGrid(a::A, b::B, n) where {A, B}
    T = promote_type(float(A), float(B))
    return RegularGrid{T}(a, b, n)
end

function nodes(g::RegularGrid{T}) where {T}
    n = g.n
    h = (g.b - g.a) / (n - 1)
    return g.a .+ h * (0:(n - 1))
end

function weights(g::RegularGrid)
    @assert g.n ≥ 4
    m = g.n % 3

    w = fill(oftype(g.a, 9 // 8), g.n)
    w[1] = 3 // 8
    for i = 4:3:(g.n - 3)
        w[i] = 3 // 4
    end
    if m == 1
        w[end    ] =  3 // 8
    elseif m == 0
        w[end - 2] = 17 // 24
        w[end - 1] =  4 // 3
        w[end    ] =  1 // 3
    elseif m == 2
        w[end - 3] =  7 // 6
        w[end - 2] = 11 // 12
        w[end - 1] =  7 // 6
        w[end    ] =  3 // 8
    end

    return w
end

jacobian(g::RegularGrid) = (g.b - g.a) / (g.n - 1)

function interpolate(g::RegularGrid, v::AbstractVector, x′)
    @assert (g.a ≤ x′ ≤ g.b)
    r = x′ .- range(g.a, stop = g.b, length = g.n)

    i = findfirst(x -> x ≤ 0, r)
    if r[i] == 0
        return v[i] / one(r[i]) # in case there are type instabilities
    end
    k = 3 * div(i - 1, 3)
    k = (k + 4 > g.n) ? (g.n - 4) : k

    # Barycentric interpolation weights
    w₁ =  1 / r[k + 1]
    w₂ = -3 / r[k + 2]
    w₃ =  3 / r[k + 3]
    w₄ = -1 / r[k + 4]
    # Interpolation function reference values
    v₁, v₂, v₃, v₄ = v[k + 1], v[k + 2], v[k + 3], v[k + 4]

    return (w₁ * v₁ + w₂ * v₂ + w₃ * v₃ + w₄ * v₄) / (w₁ + w₂ + w₃ + w₄)
end


text:
export
    Table

"""
mutable struct Table{NumInputs, I <: NTuple{NumInputs, Any}, J, K, O, S <: SFunc{J,O}} <: Conditional{I, J, K, O, S}

`Table` defines a `Conditional` whose function is given by a multi-dimensional table of type
`Array{Q, NumInputs}`, where `Q` is the output type of the internal *sfunc* and `NumInputs`
is the count of incoming parent values.

See also:  [`Conditional`](@ref), [`DiscreteCPT`](@ref), [`CLG`](@ref)
"""
mutable struct Table{NumInputs, I <: NTuple{NumInputs, Any}, J, K, O, S <: SFunc{J,O}} <: Conditional{I, J, K, O, S}
    icombos :: Vector{I}
    iranges :: NTuple{NumInputs, Array{Any, 1}}
    isizes :: NTuple{NumInputs, Int}
    imults :: NTuple{NumInputs, Int}
    inversemaps :: NTuple{NumInputs, Dict{Any, Int}}
    sfs :: Array{S, NumInputs}
    """
        function Table(J, O, NumInputs::Int, paramdict, sfmaker:Function) where {J, O, S <: SFunc{J,O}}

    `Table` constructor.

    # Arguments
    - `J` the type of inputs to the sfuncs in the table
    - 'O' the type of outputs from the sfuncs in the table
    - `NumInputs::Int`  the count of incoming parent values
    - `paramdict` see [`DiscreteCPT`](@ref) and [`CLG`](@ref) for examples
    - `sfmaker` a function from Q to S
    """
    function Table(J, O, NumInputs::Int, paramdict::Dict{I,Q}, sfmaker::Function) where {I, Q}
        K = extend_tuple_type(I,J)
        icombos = keys(paramdict)
        iranges :: Array{Array{Any,1}} = [unique(collect([combo[k] for combo in icombos])) for k in 1:NumInputs]
        isizes = tuple([length(irange) for irange in iranges]...)
        m = 1
        ims = zeros(Int, NumInputs)
        for k in NumInputs:-1:1
            ims[k] = m
            m *= isizes[k]
        end
        imults = tuple(ims...)
        # TODO: Fix ordering of Dict
        inversemaps = tuple([Dict([x => i for (i,x) in enumerate(irange)]) for irange in iranges]...)
        sortedcombos = [tuple(p...) for p in cartesian_product(iranges)]
        S = SFunc{J,O}
        sfs = Array{S, NumInputs}(undef, isizes)
        for k in 1:length(sortedcombos)
            is = sortedcombos[k]
            q = paramdict[is]
            sfs[k] = sfmaker(q)
        end
        new{NumInputs, I, J, K, O, S}(sortedcombos, tuple(iranges...), isizes, imults, inversemaps, sfs)
    end
end

#=
function dict2tableparams(sf::Table{NumInputs,I,J,K,O,S}, p::Dict{I,Q}) where {NumInputs,I,J,K,O,Q,S}
    k1 = collect(keys(p))[1]
    @assert NumInputs == length(k1)

    icombos = keys(p)
    iranges :: Array{Array{Any,1}} = [unique(collect([combo[k] for combo in icombos])) for k in 1:NumInputs]
    sortedcombos = [tuple(p...) for p in cartesian_product(iranges)]
    isizes = tuple([length(irange) for irange in iranges]...)
    array = Array{Q, NumInputs}(undef, isizes)
    for k in 1:length(sortedcombos)
        is = sortedcombos[k]
        array[k] = get(p, is, sf.default) # The input dict does not have to have all combos
    end
    array
end
=#

function gensf(t::Table{N,I,J,K,O,S}, parvals::NTuple{N,Any}) where {N,I,J,K,O,Q,S}
    inds = tuple([t.inversemaps[k][parvals[k]] for k in 1:length(parvals)]...)
    i = 1
    for k in 1:N
        i += (inds[k] - 1) * t.imults[k]
    end
    return t.sfs[i]# Change this to just index on the tuple of indices rather than do the calculation ourselves
end
#=
function do_maximize_stats(t::Table{N,I,J,K,O,Q,S}, sfmaximizers) where {N,I,J,K,O,Q,S}
    result = Array{Q, N}(undef, t.isizes)
    # Since this is a table, the new parameters have to have an entry for each of the original parent values
    for k in 1:length(t.params)
        is = t.icombos[k]
        result[k] = get(sfmaximizers, is, t.default)
    end
    return result
end
=#