text:
#=
  CyclingTime1 Example from
  "Infer.NET 101 A sample-based introduction to the basics of
  Microsoft Infer.NET programming", page 12ff.
  """
  averageTimePosterior: Gaussian(15.29, 1.559)
  trafficNoisePosterior: Gamma(1.458, 0.3944)[mean=0.5751]
  ...
  Tomorrows predicted time: 15.29 plus or minus 2.66
  ...
  Probability that the trip takes less than 18 min: 0.85
  """

  Summary Statistics
                   parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
                       Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

                  averageTime   15.3182    1.0817     0.0108    0.0357    816.4845    0.9999      209.3550
                 trafficNoise    1.8035    0.5902     0.0059    0.0223    761.6531    1.0011      195.2957
                tomorrowsTime   15.3090    2.1760     0.0218    0.0398   2891.2802    1.0000      741.3539
  probTripLongerThan18Minutes    0.9088    0.2879     0.0029    0.0046   4217.2755    0.9999     1081.3527

  Cf ~/webppl/cycling_time1.wppl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function cycling_time1(travelTimes=[13,17,16])
    averageTime ~ Normal(15,sqrt(100.0))
    trafficNoise ~ Gamma(2.0,1/2.0)

    # travelTimeMonday    ~ Normal(averageTime,trafficNoise)
    # travelTimeTuesday   ~ Normal(averageTime,trafficNoise)
    # travelTimeWednesday ~ Normal(averageTime,trafficNoise)
    # Better:
    travelTimes ~ filldist(Normal(averageTime,trafficNoise),3)
    tomorrowsTime       ~ Normal(averageTime,trafficNoise)

    probTripLongerThan18Minutes ~ Dirac(tomorrowsTime < 18.0)

end

model = cycling_time1()

num_chains = 4

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
# chns = sample(model, PG(5), 10_000)
chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 100_000)

display(chns)
# display(plot(chns))

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
abstract type InitializationStrategyC <: InitializationStrategy end

struct FixedCStrategy <: InitializationStrategyC
    C
end

calculate_C(model, strategy::FixedCStrategy) = strategy.C

"""
Original publication:
Tax, David MJ, and Robert PW Duin. "Support vector data description." Machine learning 54.1 (2004): 45-66.
"""
struct TaxErrorEstimate <: InitializationStrategyC
    target_outlier_percentage::Float64
end

calculate_C(model, strategy::TaxErrorEstimate) = 1 / (size(model.data, 2) * strategy.target_outlier_percentage)

struct BoundedTaxErrorEstimate <: InitializationStrategyC
    target_outlier_percentage::Float64
    lower_bound::Float64
    upper_bound::Float64
end

function calculate_C(model, strategy::BoundedTaxErrorEstimate)
    C = calculate_C(model, TaxErrorEstimate(strategy.target_outlier_percentage))
    return min(max(strategy.lower_bound, C), strategy.upper_bound)
end

struct BinarySearchCStrategy <: InitializationStrategyC
    target_outlier_percentage
    C_init
    max_iter
    eps
    solver
end

BinarySearchCStrategy(target_outlier_percentage, solver) = BinarySearchCStrategy(target_outlier_percentage, 0.5, 5, 0.01, solver)

function calculate_C(model, strategy::BinarySearchCStrategy, kernel)
    iteration = 1
    C_current = strategy.C_init
    C_min = 0
    C_max = 1
    debug(LOGGER, "[BINARY_SEARCH] Searching for parameter C with $(typeof(strategy)) with kernel = $(kernel).")
    debug(LOGGER, "[BINARY_SEARCH] Search parameters: C_min = $(C_min), C_max = $(C_max), max_iter = $(strategy.max_iter).")

    m = deepcopy(model)
    init_strategy = SVDD.FixedParameterInitialization(kernel, C_current)
    SVDD.initialize!(m, init_strategy)

    while true
        SVDD.fit!(m, strategy.solver)
        classification = countmap(SVDD.classify.(SVDD.predict(m, m.data)))
        current_outlier_percentage = haskey(classification, :outlier) ? classification[:outlier] / size(m.data,2) : 0.0

        if abs(current_outlier_percentage - strategy.target_outlier_percentage) < strategy.eps || iteration > strategy.max_iter
            debug(LOGGER, "[BINARY_SEARCH] Search result: C = $(C_current)")
            return C_current
        end

        if current_outlier_percentage > strategy.target_outlier_percentage
            C_min = C_current
            C_current = C_current + (C_max - C_current) / 2
        else
            C_max = C_current
            C_current = C_current - (C_current - C_min) / 2
        end
        SVDD.set_C!(m, C_current)
        iteration +=1
    end
end

ext: jl


text:
using CBindingGen

const LIBCLANG_PATH = ARGS[1]  # LLVM_jll.libclang_path
const LIBCLANG_VERSION = let
	dir = joinpath(dirname(dirname(LIBCLANG_PATH)), "lib", "clang")
	entries = readdir(dir)
	ind = findfirst(x -> !isnothing(match(r"^\d+\.\d+\.\d+$", x)) && isdir(joinpath(dir, x)), entries)
	isnothing(ind) && error("Failed to determine the version of libclang")
	entries[ind]
end

incdir = joinpath(dirname(dirname(LIBCLANG_PATH)), "include")
hdrs = map(hdr -> joinpath("clang-c", hdr), readdir(joinpath(incdir, "clang-c")))

cvts = convert_headers(hdrs, args = ["-I", incdir]) do cursor
	header = CodeLocation(cursor).file
	name   = string(cursor)

	# only wrap the libclang headers
	startswith(header, "$(incdir)/") || return false

	# ignore function that uses time_t since we don't know what time_t is yet
	name == "clang_getFileTime" && return false

	return true
end

open(joinpath(@__DIR__, "libclang-$(LIBCLANG_VERSION).jl"), "w+") do io
	generate(io, :LIBCLANG_PATH => cvts, relto = dirname(dirname(LIBCLANG_PATH)))
end


text:
import Base.Random: rand

export tournsel, tournsel!

@doc """tournsel(p::Population, k::Int64, fs::Vector{Float64})

Conduct tournament selection on the population and return a new population. For
now, we assume `p=1`, so the best individual in the tournament is selected as a
member of the new population.  Each tournament involves `k` members of the
population.

A vector of fitnesses may be passed as a third parameter to avoid re-computing
them since doing so reduces performance a bit.

Note that sampling is done with replacement, so even if `k` is equal to the
population size, there is no guarantee that every individual will participate
in each tournament.

References:

https://en.wikipedia.org/wiki/Tournament_selection
"""
function tournsel(p::Population, k::Int64, fs::Vector{Float64})
  np = Population(p)
  tournsel!(np, k, fs)
  return np
end

tournsel(p::Population, k::Int64) = tournsel(p, k, popfits(p))

function tournsel(mp::MetaPopulation, k::Int64)
  np = MetaPopulation(mp)
  tournsel!(np, k)
  return np
end

@doc """tournsel!(p::Population, k::Int64, fs::Vector{Float64})

Conduct tournament selection in-place.
"""
function tournsel!(p::Population, k::Int64, fs::Vector{Float64})
  n = popsize(p)
  selected = zeros(Int64, n)
  for i = 1:n
    contestants = rand(1:n, k)
    winner = fs[contestants] |> indmax
    selected[i] = contestants[winner]
  end
  p.genotypes[:] = [Genotype(p.genotypes[selected[i]]) for i = 1:n]
end

tournsel!(p::Population, k::Int64) = tournsel!(p, k, popfits(p))

function tournsel!(mp::MetaPopulation, k::Int64)
  for p = mp.populations
    tournsel!(p, k)
  end
end


text:
using Revise
using Random
using POMDPs
using Spot
using POMDPModelChecking
using BeliefUpdaters
using SARSOP
using RockSample
using POMDPSimulators
import Cairo
using POMDPGifs
using ProgressMeter
using Statistics
using POMDPModelTools

# pomdp = RockSamplePOMDP{2}(map_size=(4,4),
#                            rocks_positions=[(2,3), (3,1)])

# pomdp = RockSamplePOMDP{3}(map_size=(5,5),
#                             rocks_positions=[(2,3), (4,4), (4,2)])

pomdp = RockSamplePOMDP{8}(map_size=(7,7),
                           rocks_positions=[(1,2), (2,8), (3,1), (3,5), (4,2), (4,5), (6,6), (7,4)])

@show length(states(pomdp))
@show length(actions(pomdp))
@show length(observations(pomdp))

## Probability of getting at least one good rock

function POMDPModelChecking.labels(pomdp::RockSamplePOMDP, s::RSState, a::Int64)
    if a == RockSample.BASIC_ACTIONS_DICT[:sample] && in(s.pos, pomdp.rocks_positions) # sample
        rock_ind = findfirst(isequal(s.pos), pomdp.rocks_positions) # slow ?
        if s.rocks[rock_ind]
            # return ()
            return (:good_rock,)
        else
            # return (:bad_rock,)
        end
    end
    if isterminal(pomdp, s)
        return (:exit,)
        # return ()
    end
    return ()
end

# prop = ltl"G !bad_rock"
# prop = ltl"F good_rock & F exit"
prop = ltl"F good_rock & G !bad_rock & F exit"
# prop = ltl" (!bad_rock U good_rock) && (!bad_rock U exit)"

# overwriding discount, to trigger SARSOP iterations
POMDPs.discount(problem::Union{ProductMDP, ProductPOMDP}) = 0.999

# run(`rm model.pomdpx`)
solver = ModelCheckingSolver(property = prop,
                      solver=SARSOPSolver(precision=1e-3), verbose=true);

policy = solve(solver, pomdp);

## visualize policy

# first simulate the product pomdp
rng = MersenneTwister(2)
up = DiscreteUpdater(policy.problem)
b0 = initialize_belief(up, initialstate(policy.problem))
hr = HistoryRecorder(max_steps=50)
product_hist = simulate(hr, policy.problem, policy, up, b0);

# create a new history with the pomdp state and action, to be replayed and visualized
hist = SimHistory([(s=s.s, a=a) for (s, a) in eachstep(product_hist, "(s,a)")],
                  discount(pomdp), nothing, nothing)

makegif(pomdp, hist, filename="test.gif", spec="(s,a)")

#=

## Monte carlo simulation
function sim(policy)
    n_sim = 300
    successes = zeros(n_sim)
    fails = zeros(n_sim)
    @showprogress for i=1:n_sim
        up = DiscreteUpdater(policy.problem)
        b0 = initialize_belief(up, initialstate_distribution(policy.problem))
        hr = HistoryRecorder(max_steps=50)
        hist = simulate(hr, policy.problem, policy, up, b0)
        println("discounted reward ", discounted_reward(hist))
        println("undiscounted_reward", undiscounted_reward(hist))
        mu_succ, std_mu_succ = running_stats(successes[1:i])
        println("mean: ", mu_succ[end])
        successes[i] = undiscounted_reward(hist) > 0.
        fails[i] = undiscounted_reward(hist) <= 0.
    end
    return successes, fails
end

successes, fails = sim(policy)

function running_stats(vec)
    mu = mean.(vec[1:i] for i=1:length(vec))
    std_mu = std.(mu[1:i] for i=1:length(vec))./collect(1:length(vec))
    return mu, std_mu
end

mu_succ, std_mu_succ = running_stats(successes)
# mu_fail, std_mu_fail = running_stats(fails)

LBOUND = 0.874038; UBOUND=  0.874175

using Plots
pgfplots()

p = plot(mu_succ, yerror=std_mu_succ, label="MC estimate", legend=:bottom, linewidth=3);
plot!(1:length(mu_succ), LBOUND*ones(length(mu_succ)), c="black", linestyle=:dot, label="SARSOP lower bound");
plot!(1:length(mu_succ), UBOUND*ones(length(mu_succ)), c="black", linestyle=:dash, label="SARSOP upper bound")

savefig(p, "succ.png")

dra = DeterministicRabinAutomata(prop)
sink_state = ProductState(first(states(pomdp)), -1)
pmdp = ProductPOMDP(pomdp, dra, sink_state)

mecs = maximal_end_components(pmdp, verbose=true)

states(pmdp)[801]

s = ProductState(pomdp.terminal_state, 1)
labels(pmdp.problem, s.s, 1)

nextstate(pmdp.automata, 1, (:exit,))

using LightGraphs, MetaGraphs
dra = pmdp.automata

nextstate(dra, 1, (:exit,))

lab = (:exit,)

neighbors(dra.transition, 1)

edge_it = filter_edges(dra.transition,
                        (g, e) -> (src(e) == 1) && (lab ∈ props(g, e)[:cond] || (:true_constant,) ∈ props(g, e)[:cond]))

collect(edge_it)

props(dra.transition, 1, 1)

si = stateindex(pmdp, s)

for a in actions(pmdp)
    println(length(support(transition(pmdp, s, a))))
    println(first(support(transition(pmdp, s, a))))
end

for (i, mec) in enumerate(mecs)
    if si in mec
        println(i)
    end
end

last(states(pmdp))

dra.acc_sets

dra = DeterministicRabinAutomata(prop)

translator = LTLTranslator(deterministic=true, buchi=true, state_based_acceptance=true)
aut = translate(translator, prop)

dra = SpotAutomata(aut)
@assert is_deterministic(dra)
states = 1:num_states(dra)
initial_state = get_init_state_number(dra)
APs = atomic_propositions(dra)
edgelist, llabels = get_edges_labels(dra)
sdg = SimpleDiGraph(num_states(dra))
for e in edgelist
    add_edge!(sdg, e)
end
transition = MetaDiGraph(sdg)
conditions = label_to_array.(llabels)

rng = MersenneTwister(2)
up = DiscreteUpdater(policy.problem)
b0 = initialize_belief(up, initialstate_distribution(policy.problem))
hr = HistoryRecorder(max_steps=50)
hist = simulate(hr, policy.problem, policy, up, b0);
prod_state_hist = hist.state_hist
state_hist = [s.s for s in hist.state_hist];
hist = POMDPHistory(state_hist, [getfield(hist, f) for f in fieldnames(POMDPHistory) if f != :state_hist]...);
makegif(pomdp, hist, filename="test.gif", spec="(s,a)")

a = Spot.translate(LTLTranslator(buchi=true, deterministic=true, state_based_acceptance=true), prop)
spot.to_generalized_rabin(a)
display(a)

a = Spot.translate(LTLTranslator(buchi=false, deterministic=true, state_based_acceptance=true), prop)

=#


text:
module UnitfulCooking

export VOL_DATA, volume_to_mass, mass_to_volume

using DataFrames, Query
using CSV
using Unitful

function read_volume_data()
  parse_euro_num(s) = parse(Float64, replace(s, r"\." => ""))
  path = joinpath(dirname(@__FILE__), "..", "data", "volume-weight-conversion-v2.csv")
  vd = CSV.read(path, DataFrame)[3:end,:] |>
    @map({ingredient=_.Column16, grams_per_liter=parse_euro_num(_[Symbol(" l ")])}) |>
    DataFrame
end

const VOL_DATA = read_volume_data()

function grams_per_liter(ingredient)
  matches = VOL_DATA |> @filter(occursin(ingredient, _.ingredient[])) |> DataFrame
  if length(Tables.rows(matches)) == 0
    error("404 Ingredient not found")
  end
  gpl = matches[1,:grams_per_liter]
end

function volume_to_mass(ingredient, vol, output_units=u"g")
  vol_in_liters = convert(Float64, vol / u"l")
  mass = grams_per_liter(ingredient) * vol_in_liters * u"g"
  uconvert(output_units, mass)
end

function mass_to_volume(ingredient, mass, output_units=u"l")
  mass_in_grams = convert(Float64, mass / u"g")
  vol = (mass_in_grams / grams_per_liter(ingredient)) * u"l"
  uconvert(output_units, vol)
end

end # module


text:
# sum will give the integral over the field
function Base.sum(field::Union{Field, Base.Broadcast.Broadcasted{<:FieldStyle}})
    Base.reduce(
        RecursiveApply.radd,
        Base.Broadcast.broadcasted(
            RecursiveApply.rmul,
            weighted_jacobian(field),
            todata(field),
        ),
    )
end


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function testlabels(fname)
    pagewidth, pageheight = 800, 800

    Drawing(pagewidth, pageheight, fname)
    origin() # move 0/0 to center
    background("ivory")
    sethue("black")

    # compasspoint version
    compass = [:E, :SE, :S, :SW, :W, :NW, :N, :NE]

    fontsize(7)
    setline(1)
    for i in 0:pi/64:2pi
        pt = polar(250, i)
        line(between(O, pt, 0.5), pt, :stroke)
        circle(pt, 2, :fill)
        n = mod1(convert(Int, div(i + pi/8, pi/4) + 1), 8)
        @layer begin
            setline(.5)
            label(string(compass[n]),
                compass[n],
                pt,
                offset=30,
                leaderoffsets=[0.2, 0.8],
                leader = isodd(n) ? true : false)
        end
    end

    # angle version
    fontsize(4)
    setline(0.5)
    pg = ngon(O, 80, 32, 0, vertices=true)
    poly(pg, :stroke, close=true)
    for pt in pg
        circle(pt, 3, :fill)
        s = slope(O, pt)
        label(string(round(s, digits=1)), s, pt, leader=true, leaderoffsets=[0.3, 0.9], offset=15)
        label(string(round(s, digits=1)), s, pt, offset=-15)
    end

    @test finish() == true
    println("...finished labels test, saved in $(fname)")
end

fname = "label-test.svg"
testlabels(fname)


text:
using Random

mutable struct EpsilonGreedy <: ExplorationStrategy
    epsilon::Number
end

function get_action_probs(EG::EpsilonGreedy, qs)
    m = maximum(qs)
    probs = zeros(length(qs))

    maxes_ind = findall(x-> x==m, qs)
    for ind in maxes_ind
        probs[ind] += (1 - EG.epsilon) / size(maxes_ind)[1]
    end
    probs .+= EG.epsilon/ length(qs)
    return probs
end

function step!(EG::EpsilonGreedy)
end

(e::EpsilonGreedy)(qs) = get_action_probs(e, qs)
(e::EpsilonGreedy)(qs, action) = get_action_probs(e, qs)[action]

"""
    Taken from Matt Schlegel Action RNN repo:
    https://github.com/mkschleg/ActionRNNs.jl/blob/master/src/ActingPolicy.jl
    Commit Hash: 5fa7dd8ace97014d7240d1e857c113d6f0c00e4e

    ϵGreedyDecay{AS}(ϵ_range, decay_period, warmup_steps, action_set::AS)
    ϵGreedyDecay(ϵ_range, end_step, num_actions)
This is an acting policy which decays exploration linearly over time. This api will possibly change overtime once I figure out a better way to specify decaying epsilon.
# Arguments
`ϵ_range::Tuple{Float64, Float64}`: (max epsilon, min epsilon)
`decay_period::Int`: period epsilon decays
`warmup_steps::Int`: number of steps before decay starts
"""
Base.@kwdef mutable struct ϵGreedyDecay{AS} <: ExplorationStrategy
    ϵ_range::Tuple{Float64, Float64}
    decay_period::Int
    warmup_steps::Int
    cur_step::Int = 0
    action_set::AS
    ϵGreedyDecay(ϵ_range, decay_period, warmup_steps, action_set::AS) where {AS} =
        new{AS}(ϵ_range, decay_period, warmup_steps, 0, action_set)
end

ϵGreedyDecay(ϵ_range, end_step, num_actions) = ϵGreedyDecay(ϵ_range, end_step, 1:num_actions)

action_set(ap::ϵGreedyDecay) = ap.action_set

function _get_eps_for_step(ap::ϵGreedyDecay, step=ap.cur_step)
    ϵ_min = ap.ϵ_range[2]
    ϵ_max = ap.ϵ_range[1]

    steps_left = ap.decay_period + ap.warmup_steps - ap.cur_step
    bonus = (ϵ_max - ϵ_min) * steps_left / ap.decay_period
    bonus = clamp(bonus, 0.0, ϵ_max - ϵ_min)
    ϵ_min + bonus
end

function get_action_probs(ap::ϵGreedyDecay, qs, step=ap.cur_step)
    ϵ = _get_eps_for_step(ap, step)

    m = maximum(qs)
    probs = zeros(length(qs))

    maxes_ind = findall(x-> x==m, qs)
    for ind in maxes_ind
        probs[ind] += (1 - ϵ) / size(maxes_ind)[1]
    end
    probs .+= ϵ/ length(qs)
    return probs
end

function step!(ap::ϵGreedyDecay)
    ap.cur_step += 1
end
(e::ϵGreedyDecay)(qs) = get_action_probs(e, qs)
(e::ϵGreedyDecay)(qs, action) = get_action_probs(e, qs)[action]

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
"""
Plane wave single-particle basis states.
"""
module PlaneWaves

export Spin, Down, Up, PlaneWave, dimension, flip, sphere, sphere_with_same_spin, ξ, fractional_spin_polarization, find_fourth_orb_for_kink

using ..CPIMC
using StaticArrays

"""
Type for storing spin projection of particles with spin 1/2.
"""
@enum Spin Down Up

"""
    flip(s::Spin)

return opposite spin projection
"""
function flip(s::Spin)
    if s == Down
        Up
    else
        Down
    end
end

"""
    struct PlaneWave{D} <: Orbital

Representation of plane wave single-particle states with dimensionality `D`.

**Fields**
- `vec  :: StaticVector{D,Int}` -- momentum vector
- `spin :: Spin`                -- spin

"""
struct PlaneWave{D} <: Orbital
    vec :: SVector{D,Int}
    spin :: Spin
end

"""
    PlaneWave(v::Tuple,s=Up)

Outer constructor with default value `spin = Up`.
"""
PlaneWave(v::Tuple,s=Up) = PlaneWave(SVector(v),s)

"""
    flip(o::PlaneWave)
return an PlaneWave with the same momentum vector but opposite spin projection as the given o::PlaneWave """
flip(o::PlaneWave) = PlaneWave(tuple(o.vec...),flip(o.spin))

"""
    dimension(o::PlaneWave{D}) where {D}
return the dimension of the momentum vector of an Orbital"""
dimension(o::PlaneWave{D}) where {D} = D

"""
    dimension(os::Set{<:PlaneWave{D}}) where {D}
return the dimension of the momentum vectors of a Set{<:Orbital{D}}
"""
dimension(os::Set{<:PlaneWave{D}}) where {D} = D

"""
    sphere_with_same_spin(o::PlaneWave{1}; dk=2)
    sphere_with_same_spin(o::PlaneWave{2}; dk=2)
    sphere_with_same_spin(o::PlaneWave{3}; dk=2)

Return the set of all PlaneWaves in a sphere of radius `dk` around `o` in momentum space that also have the same spin as `o`.
"""
function sphere_with_same_spin(o::PlaneWave{1}; dk::Number=2)
    os = Set{PlaneWave{1}}()

    for x in -Int(ceil(dk)):Int(ceil(dk))
        if x*x <= dk*dk
            push!(os, PlaneWave(o.vec+SVector(x),o.spin))
        end
    end
    os
end

function sphere_with_same_spin(o::PlaneWave{2}; dk::Number=2)
    os = Set{PlaneWave{2}}()

    for x in -Int(ceil(dk)):Int(ceil(dk))
        for y in -Int(ceil(dk)):Int(ceil(dk))
            if x*x + y*y <= dk*dk
                push!(os, PlaneWave(o.vec+SVector(x,y),o.spin))
            end
        end
    end
    os
end

function sphere_with_same_spin(o::PlaneWave{3}; dk::Number=2)
    os = Set{PlaneWave{3}}()

    for x in -Int(ceil(dk)):Int(ceil(dk))
        for y in -Int(ceil(dk)):Int(ceil(dk))
            for z in -Int(ceil(dk)):Int(ceil(dk))
                if x*x + y*y + z*z <= dk*dk
                    push!(os, PlaneWave(o.vec+SVector(x,y,z),o.spin))
                end
            end
        end
    end
    os
end

"""
    sphere(o::PlaneWave{1}; dk=2)
    sphere(o::PlaneWave{2}; dk=2)
    sphere(o::PlaneWave{3}; dk=2)

Return the set of all PlaneWaves in a sphere of radius `dk` around `o` in momentum space.
"""
sphere(o::PlaneWave; dk=2) = union(sphere_with_same_spin(o,dk=dk),sphere_with_same_spin(PlaneWave(o.vec,flip(o.spin)),dk=dk))

"""
    fractional_spin_polarization(occ::Set{PlaneWave)
calculate the fractional spin polarization from a set of Orbitals which each must have a field `spin` of type `@enum Spin Down Up`,
the fractional spin polarization is defined as the ratio of the absolute difference of the number of particles which Spin Down and Spin Up and the particle number,
thus here no convention is made as to which spin component should be occupied more often, as opposed to some literature where N↑ > N↓ is used
"""
function fractional_spin_polarization(occ::Set{<:PlaneWave})
    N_up = length(filter(x -> x.spin == Up, occ))
    N_down = length(filter(x -> x.spin == Down, occ))
    N = length(occ)
    return abs(N_up - N_down) / N
end

"""
    ξ(occ::Set{PlaneWave})
calculate the fractional spin polarization from a set of Orbitals which each must have a field `spin` of type `@enum Spin Down Up`,
the fractional spin polarization is defined as the ratio of the absolute difference of the number of particles which Spin Down and Spin Up and the particle number,
thus here no convention is made as to which spin component should be occupied more often, as opposed to some literature where N↑ > N↓ is used
"""
ξ(occ::Set{<:PlaneWave}) = fractional_spin_polarization(occ)

"""
    fractional_spin_polarization(c::Configuration{<:PlaneWave})
calculate the fractional spin polarization from a set of Orbitals which each must have a field `spin` of type `@enum Spin Down Up`,
the fractional spin polarization is defined as the ratio of the absolute difference of the number of particles which Spin Down and Spin Up and the particle number,
thus here no convention is made as to which spin component should be occupied more often, as opposed to some literature where N↑ > N↓ is used
"""
function fractional_spin_polarization(c::Configuration{<:PlaneWave})# TODO: added {<:PlaneWave}
    fractional_spin_polarization(c.occupations)
end

"""
    ξ(c::Configuration{<:PlaneWave})
calculate the fractional spin polarization from a set of Orbitals which each must have a field `spin` of type `@enum Spin Down Up`,
the fractional spin polarization is defined as the ratio of the absolute difference of the number of particles which Spin Down and Spin Up and the particle number,
thus here no convention is made as to which spin component should be occupied more often, as opposed to some literature where N↑ > N↓ is used
"""
ξ(c::Configuration{<:PlaneWave}) = fractional_spin_polarization(c)

"""
    find_fourth_orb_for_kink(same_kind_ladder_operator, other_kind_ladder_operator1, other_kind_ladder_operator2)

Find the fourth orb to build a Kink accoring to spin/momentum conservation, where the Operator thats acts on the
fourth orb is of the same kind (creator/ annihilator) as the one that acts on same_kind_ladder_operator,
While on other_kind_ladder_operator1 and other_kind_ladder_operator2 the other kind of operators will act.
"""
function find_fourth_orb_for_kink(same_kind_ladder_operator, other_kind_ladder_operator1, other_kind_ladder_operator2)
    @assert(in(same_kind_ladder_operator.spin, [other_kind_ladder_operator1.spin, other_kind_ladder_operator2.spin]))
    if other_kind_ladder_operator1.spin == other_kind_ladder_operator2.spin
        fourth_orb = PlaneWave(other_kind_ladder_operator1.vec + other_kind_ladder_operator2.vec - same_kind_ladder_operator.vec, same_kind_ladder_operator.spin)

    else
        fourth_orb = PlaneWave(other_kind_ladder_operator1.vec + other_kind_ladder_operator2.vec - same_kind_ladder_operator.vec, flip(same_kind_ladder_operator.spin))
    end
    return(fourth_orb)
end

end


text:
function triangulate(_dataset::AbstractVTKStructuredData, decompose_cell_data::Bool=false)
    dataset = VTKStructuredData(_dataset)
    _dim = dim(dataset)

    if decompose_cell_data
        if _dim == 2
            return triangulate_2d_with_cell_data(dataset)
        elseif _dim == 3
            return triangulate_3d_with_cell_data(dataset)
        else
            throw("Unsupported dimension.")
        end
    else
        if _dim == 2
            return triangulate_2d_no_cell_data(dataset)
        elseif _dim == 3
            return triangulate_3d_no_cell_data(dataset)
        else
            throw("Unsupported dimension.")
        end
    end
end

function triangulate_2d_with_cell_data(dataset::VTKStructuredData)
    cextents = cell_extents(dataset)
    _num_of_points = num_of_points(dataset)
    _num_of_cells = num_of_cells(dataset)
    _dim = dim(dataset)

    point_coords = reshape(dataset.point_coords, (_dim, _num_of_points))
    point_data = empty(dataset.point_data)
    for m in keys(dataset.point_data)
        _var_dim = var_dim(dataset, m, "Point")
        if _var_dim == 1
            point_data[m] = reshape(dataset.point_data[m], (_num_of_points,))
        else
            point_data[m] = reshape(dataset.point_data[m], (_var_dim, _num_of_points))
        end
    end

    _cell_data = empty(dataset.cell_data)
    faces = Vector{Int}[]

    for m in keys(dataset.cell_data)
        _var_dim = var_dim(dataset, m, "Cell")
        if _var_dim == 1
            _cell_data[m] = Float64[]
        else
            _cell_data[m] = Float64[]
        end
    end

    _cell_types = fill(5, 2*_num_of_cells)
    for cind in Iterators.product(1:cextents[1], 1:cextents[2])
        quad_cc = cell_connectivity(dataset, cind)
        push!(faces, [quad_cc[1], quad_cc[2], quad_cc[3]])
        push!(faces, [quad_cc[1], quad_cc[3], quad_cc[4]])
        for m in keys(dataset.cell_data)
            _var_dim = var_dim(dataset, m, "Cell")
            if _var_dim == 1
                append!(_cell_data[m], dataset.cell_data[m][cind...])
                append!(_cell_data[m], dataset.cell_data[m][cind...])
            else
                append!(_cell_data[m], dataset.cell_data[m][:, cind...])
                append!(_cell_data[m], dataset.cell_data[m][:, cind...])
            end
        end
    end

    for m in keys(dataset.cell_data)
        _var_dim = var_dim(dataset, m, "Cell")
        if _var_dim == 1
            _cell_data[m] = reshape(_cell_data[m], (2*_num_of_cells,))
        else
            _cell_data[m] = reshape(_cell_data[m], (_var_dim, 2*_num_of_cells))
        end
    end

    return VTKPolyData(point_coords, _cell_types, faces, point_data, _cell_data)
end

function triangulate_3d_with_cell_data(dataset::VTKStructuredData)
    cextents = cell_extents(dataset)
    _num_of_points = num_of_points(dataset)
    _dim = dim(dataset)

    point_coords = reshape(dataset.point_coords, (_dim, _num_of_points))
    point_data = empty(dataset.point_data)
    for m in keys(dataset.point_data)
        _var_dim = var_dim(dataset, m, "Point")
        if _var_dim == 1
            point_data[m] = reshape(dataset.point_data[m], (_num_of_points,))
        else
            point_data[m] = reshape(dataset.point_data[m], (_var_dim, _num_of_points))
        end
    end

    _cell_data = empty(dataset.cell_data)
    faces = Vector{Int}[]
    for m in keys(dataset.cell_data)
        _cell_data[m] = Float64[]
    end

    for cind in Iterators.product(1:cextents[1], 1:cextents[2], 1:cextents[3])
        hexa_cc = cell_connectivity(dataset, cind)
        push!(faces, [hexa_cc[1], hexa_cc[2], hexa_cc[3]])
        push!(faces, [hexa_cc[1], hexa_cc[3], hexa_cc[4]])
        for m in keys(dataset.cell_data)
            _var_dim = var_dim(dataset, m, "Cell")
            if _var_dim == 1
                if cind[3] != 1
                    _cd = (dataset.cell_data[m][cind...] +
                        dataset.cell_data[m][cind[1:2]...,cind[3]-1]) / 2
                else
                    _cd = dataset.cell_data[m][cind...]
                end
                push!(_cell_data[m], _cd)
                push!(_cell_data[m], _cd)
            else
                if cind[3] != 1
                    _cd = (dataset.cell_data[m][:, cind...] +
                        dataset.cell_data[m][:, cind[1:2]...,cind[3]-1]) / 2
                else
                    _cd = dataset.cell_data[m][:, cind...]
                end
                append!(_cell_data[m], _cd)
                append!(_cell_data[m], _cd)
            end
        end

        push!(faces, [hexa_cc[1], hexa_cc[2], hexa_cc[6]])
        push!(faces, [hexa_cc[1], hexa_cc[6], hexa_cc[5]])
        for m in keys(dataset.cell_data)
            _var_dim = var_dim(dataset, m, "Cell")
            if _var_dim == 1
                if cind[2] != 1
                    _cd = (dataset.cell_data[m][cind...] +
                        dataset.cell_data[m][cind[1],cind[2]-1,cind[3]]) / 2
                else
                    _cd = dataset.cell_data[m][cind...]
                end
                push!(_cell_data[m], _cd)
                push!(_cell_data[m], _cd)
            else
                if cind[3] != 1
                    _cd = (dataset.cell_data[m][:, cind...] +
                        dataset.cell_data[m][:, cind[1],cind[2]-1,cind[3]]) / 2
                else
                    _cd = dataset.cell_data[m][:, cind...]
                end
                append!(_cell_data[m], _cd)
                append!(_cell_data[m], _cd)
            end
        end

        push!(faces, [hexa_cc[4], hexa_cc[1], hexa_cc[5]])
        push!(faces, [hexa_cc[4], hexa_cc[5], hexa_cc[8]])
        for m in keys(dataset.cell_data)
            _var_dim = var_dim(dataset, m, "Cell")
            if _var_dim == 1
                if cind[1] != 1
                    _cd = (dataset.cell_data[m][cind...] +
                        dataset.cell_data[m][cind[1]-1,cind[2:3]...]) / 2
                else
                    _cd = dataset.cell_data[m][cind...]
                end
                push!(_cell_data[m], _cd)
                push!(_cell_data[m], _cd)
            else
                if cind[3] != 1
                    _cd = (dataset.cell_data[m][:, cind...] +
                        dataset.cell_data[m][:, cind[1]-1,cind[2:3]...]) / 2
                else
                    _cd = dataset.cell_data[m][:, cind...]
                end
                append!(_cell_data[m], _cd)
                append!(_cell_data[m], _cd)
            end
        end

        if cind[1] == cextents[1]
            push!(faces, [hexa_cc[2], hexa_cc[3], hexa_cc[7]])
            push!(faces, [hexa_cc[2], hexa_cc[7], hexa_cc[6]])
            for m in keys(dataset.cell_data)
                _var_dim = var_dim(dataset, m, "Cell")
                if _var_dim == 1
                    _cd = dataset.cell_data[m][cind...]
                    push!(_cell_data[m], _cd)
                    push!(_cell_data[m], _cd)
                else
                    _cd = dataset.cell_data[m][:, cind...]
                    append!(_cell_data[m], _cd)
                    append!(_cell_data[m], _cd)
                end
            end
        end
        if cind[2] == cextents[2]
            push!(faces, [hexa_cc[4], hexa_cc[3], hexa_cc[7]])
            push!(faces, [hexa_cc[4], hexa_cc[7], hexa_cc[8]])
            for m in keys(dataset.cell_data)
                _var_dim = var_dim(dataset, m, "Cell")
                if _var_dim == 1
                    _cd = dataset.cell_data[m][cind...]
                    push!(_cell_data[m], _cd)
                    push!(_cell_data[m], _cd)
                else
                    _cd = dataset.cell_data[m][:, cind...]
                    append!(_cell_data[m], _cd)
                    append!(_cell_data[m], _cd)
                end
            end
        end
        if cind[3] == cextents[3]
            push!(faces, [hexa_cc[5], hexa_cc[6], hexa_cc[7]])
            push!(faces, [hexa_cc[5], hexa_cc[7], hexa_cc[8]])
            for m in keys(dataset.cell_data)
                _var_dim = var_dim(dataset, m, "Cell")
                if _var_dim == 1
                    _cd = dataset.cell_data[m][cind...]
                    push!(_cell_data[m], _cd)
                    push!(_cell_data[m], _cd)
                else
                    _cd = dataset.cell_data[m][:, cind...]
                    append!(_cell_data[m], _cd)
                    append!(_cell_data[m], _cd)
                end
            end
        end
    end
    _cell_types = fill(5, length(faces))

    for m in keys(dataset.cell_data)
        _var_dim = var_dim(dataset, m, "Cell")
        if _var_dim == 1
            _cell_data[m] = reshape(_cell_data[m], (length(faces),))
        else
            _cell_data[m] = reshape(_cell_data[m], (_var_dim, length(faces)))
        end
    end

    return VTKPolyData(point_coords, _cell_types, faces, point_data, _cell_data)
end

function triangulate_2d_no_cell_data(dataset::VTKStructuredData)
    cextents = cell_extents(dataset)
    _num_of_points = num_of_points(dataset)
    _num_of_cells = num_of_cells(dataset)
    _dim = dim(dataset)

    point_coords = reshape(dataset.point_coords, (_dim, _num_of_points))
    point_data = empty(dataset.point_data)
    for m in keys(dataset.point_data)
        _var_dim = var_dim(dataset, m, "Point")
        if _var_dim == 1
            point_data[m] = reshape(dataset.point_data[m], (_num_of_points,))
        else
            point_data[m] = reshape(dataset.point_data[m], (_var_dim, _num_of_points))
        end
    end

    _cell_data = empty(dataset.cell_data)
    faces = Vector{Int}[]

    _cell_types = fill(5, 2*_num_of_cells)
    for cind in Iterators.product(1:cextents[1], 1:cextents[2])
        quad_cc = cell_connectivity(dataset, cind)
        push!(faces, [quad_cc[1], quad_cc[2], quad_cc[3]])
        push!(faces, [quad_cc[1], quad_cc[3], quad_cc[4]])
    end

    return VTKPolyData(point_coords, _cell_types, faces, point_data, _cell_data)
end

function triangulate_3d_no_cell_data(dataset::VTKStructuredData)
    cextents = cell_extents(dataset)
    _num_of_points = num_of_points(dataset)
    _dim = dim(dataset)

    point_coords = reshape(dataset.point_coords, (_dim, _num_of_points))
    point_data = empty(dataset.point_data)
    for m in keys(dataset.point_data)
        _var_dim = var_dim(dataset, m, "Point")
        if _var_dim == 1
            point_data[m] = reshape(dataset.point_data[m], (_num_of_points,))
        else
            point_data[m] = reshape(dataset.point_data[m], (_var_dim, _num_of_points))
        end
    end

    _cell_data = empty(dataset.cell_data)
    faces = Vector{Int}[]

    for cind in Iterators.product(1:cextents[1], 1:cextents[2], 1:cextents[3])
        hexa_cc = cell_connectivity(dataset, cind)
        push!(faces, [hexa_cc[1], hexa_cc[2], hexa_cc[3]])
        push!(faces, [hexa_cc[1], hexa_cc[3], hexa_cc[4]])

        push!(faces, [hexa_cc[1], hexa_cc[2], hexa_cc[6]])
        push!(faces, [hexa_cc[1], hexa_cc[6], hexa_cc[5]])

        push!(faces, [hexa_cc[4], hexa_cc[1], hexa_cc[5]])
        push!(faces, [hexa_cc[4], hexa_cc[5], hexa_cc[8]])

        if cind[1] == cextents[1]
            push!(faces, [hexa_cc[2], hexa_cc[3], hexa_cc[7]])
            push!(faces, [hexa_cc[2], hexa_cc[7], hexa_cc[6]])
        end
        if cind[2] == cextents[2]
            push!(faces, [hexa_cc[4], hexa_cc[3], hexa_cc[7]])
            push!(faces, [hexa_cc[4], hexa_cc[7], hexa_cc[8]])
        end
        if cind[3] == cextents[3]
            push!(faces, [hexa_cc[5], hexa_cc[6], hexa_cc[7]])
            push!(faces, [hexa_cc[5], hexa_cc[7], hexa_cc[8]])
        end
    end
    _cell_types = fill(5, length(faces))

    return VTKPolyData(point_coords, _cell_types, faces, point_data, _cell_data)
end


text:
mutable struct MCIntegrationVI{T<:Real,N,Tx,Ty} <: NumericalVI{T}
    nMC::Int64 #Number of samples for MC Integrations
    clipping::T
    ϵ::T #Convergence criteria
    nIter::Integer #Number of steps performed
    stoch::Bool # Flag for stochastic optimization
    nSamples::Int #Number of samples of the data
    nMinibatch::Int #Size of mini-batches
    ρ::T # Scaling coeff. for stoch. opt.
    NaturalGradient::Bool
    HyperParametersUpdated::Bool # Flag for updating kernel matrix
    vi_opt::NTuple{N,NVIOptimizer}
    MBIndices::Vector{Int} #Indices of the minibatch
    xview::Tx
    yview::Ty
    function MCIntegrationVI{T}(
        ϵ::T,
        nMC::Int,
        optimiser,
        Stochastic::Bool,
        clipping::Real,
        nMinibatch::Int,
        natural::Bool,
    ) where {T<:Real}
        return new{T,1,Vector{T},Vector{T}}(
            nMC,
            clipping,
            ϵ,
            0,
            Stochastic,
            0,
            nMinibatch,
            one(T),
            natural,
            true,
            (NVIOptimizer{T}(0, 0, optimiser),),
        )
    end
    function MCIntegrationVI{T}(
        ϵ::Real,
        Stochastic::Bool,
        nMC::Int,
        clipping::Real,
        nFeatures::Vector{<:Int},
        nSamples::Int,
        nMinibatch::Int,
        nLatent::Int,
        optimiser,
        natural::Bool,
        xview::Tx,
        yview::Ty,
    ) where {T,Tx,Ty}
        vi_opts = ntuple(i -> NVIOptimizer{T}(nFeatures[i], nMinibatch, optimiser), nLatent)
        return new{T,nLatent,Tx,Ty}(
            nMC,
            clipping,
            ϵ,
            0,
            Stochastic,
            nSamples,
            nMinibatch,
            T(nSamples / nMinibatch),
            natural,
            true,
            vi_opts,
            1:nMinibatch,
            xview,
            yview,
        )
    end
end

"""
    MCIntegrationVI(;ϵ::T=1e-5, nMC::Integer=1000, clipping::Real=Inf, natural::Bool=true, optimiser=Momentum(0.001))

Variational Inference solver by approximating gradients via MC Integration.
It means the expectation `E[log p(y|f)]` as well as its gradients is computed
by sampling from q(f).

## Keyword arguments
- `ϵ::Real` : convergence criteria, which can be user defined
- `nMC::Int` : Number of samples per data point for the integral evaluation
- `clipping::Real` : Limit the gradients values to avoid overshooting
- `natural::Bool` : Use natural gradients
- `optimiser` : Optimiser used for the variational updates. Should be an Optimiser object from the [Flux.jl](https://github.com/FluxML/Flux.jl) library, see list here [Optimisers](https://fluxml.ai/Flux.jl/stable/training/optimisers/) and on [this list](https://github.com/theogf/AugmentedGaussianProcesses.jl/tree/master/src/inference/optimisers.jl). Default is `Momentum(0.01)`
"""
MCIntegrationVI

function MCIntegrationVI(;
    ϵ::T=1e-5,
    nMC::Integer=1000,
    optimiser=Momentum(0.01),
    clipping::Real=Inf,
    natural::Bool=true,
) where {T<:Real}
    return MCIntegrationVI{T}(ϵ, nMC, optimiser, false, clipping, 1, natural)
end

"""
    MCIntegrationSVI(nMinibatch::Int; ϵ::Real=1e-5, nMC::Integer=1000, clipping=Inf, natural=true, optimiser=Momentum(0.0001))

Stochastic Variational Inference solver by approximating gradients via Monte Carlo integration when using minibatches
See [`MCIntegrationVI`](@ref) for more explanations.

## Argument

-`nMinibatch::Integer` : Number of samples per mini-batches

## Keyword arguments

- `ϵ::T` : convergence criteria, which can be user defined
- `nMC::Int` : Number of samples per data point for the integral evaluation
- `clipping::Real` : Limit the gradients values to avoid overshooting
- `natural::Bool` : Use natural gradients
- `optimiser` : Optimiser used for the variational updates. Should be an Optimiser object from the [Flux.jl](https://github.com/FluxML/Flux.jl) library, see list here [Optimisers](https://fluxml.ai/Flux.jl/stable/training/optimisers/) and on [this list](https://github.com/theogf/AugmentedGaussianProcesses.jl/tree/master/src/inference/optimisers.jl). Default is `Momentum()` (ρ=(τ+iter)^-κ)
"""
MCIntegrationSVI

function MCIntegrationSVI(
    nMinibatch::Integer;
    ϵ::T=1e-5,
    nMC::Integer=200,
    optimiser=Momentum(0.001),
    clipping::Real=0.0,
    natural::Bool=true,
) where {T<:Real}
    return MCIntegrationVI{T}(ϵ, nMC, optimiser, true, clipping, nMinibatch, natural)
end

function tuple_inference(
    i::MCIntegrationVI{T},
    nLatent::Int,
    nFeatures::Vector{<:Int},
    nSamples::Int,
    nMinibatch::Int,
    xview::Tx,
    yview::Ty,
) where {T,Tx,Ty}
    return MCIntegrationVI{T}(
        conv_crit(i),
        isStochastic(i),
        i.nMC,
        i.clipping,
        nFeatures,
        nSamples,
        nMinibatch,
        nLatent,
        i.vi_opt[1].optimiser,
        i.NaturalGradient,
        xview,
        yview,
    )
end

function grad_expectations!(m::AbstractGP{T,L,<:MCIntegrationVI{T,N}}) where {T,L,N}
    raw_samples = randn(T, inference(m).nMC, nLatent(m))
    samples = similar(raw_samples)
    μ = mean_f(m)
    σ² = var_f(m)
    nSamples = length(MBIndices(m))
    for j in 1:nSamples # Loop over every data point
        samples .= raw_samples .* sqrt.([σ²[k][j] for k in 1:nLatent(m)])' .+ [μ[k][j] for k in 1:nLatent(m)]'
        grad_samples(m, samples, j) # Compute the gradient for data point j
    end
end

function expec_loglikelihood(
    l::AbstractLikelihood, i::MCIntegrationVI{T,N}, y, μ, σ²
) where {T,N} # μ and σ² are tuples of vectors
    raw_samples = randn(T, i.nMC, N) # dimension nMC x nLatent
    # samples = similar(raw_samples)
    nSamples = length(MBIndices(i))
    tot = 0.0
    for j in 1:nSamples # Loop over every data point
        samples = raw_samples .* sqrt.([σ²[k][j] for k in 1:N])' .+ [μ[k][j] for k in 1:N]'
        # samples is of dimension nMC x nLatent again
        y_j = getindex.(y, j) # Obtain the label for data point j
        for f in eachrow(samples)
            # return sum(eachrow(samples)) do f # We now compute the loglikelihood over every sample
            tot += loglikelihood(l, y_j, f)
        end
    end
    return tot / i.nMC
end

# function expec_loglikelihood(
#     l::AbstractLikelihood, i::MCIntegrationVI{T,N}, y, μ, σ²
# ) where {T,N}
#     raw_samples = randn(T, i.nMC, N) # dimension nMC x nLatent
#     # samples = similar(raw_samples)
#     nSamples = length(MBIndices(i))
#     return sum(expec_loglikelihood_datapoint.(l, Ref(raw_samples), Ref(y), Ref(μ), Ref(σ²), 1:nSamples, N)) / i.nMC
#     # return sum(1:nSamples) do j # Loop over every data point
#     #     samples = raw_samples .* sqrt.([σ²[k][j] for k in 1:N])' .+ [μ[k][j] for k in 1:N]'
#     #     # samples is of dimension nMC x nLatent again
#     #     y_j = getindex.(y, j) # Obtain the label for data point j
#     #     return sum(eachrow(samples)) do f # We now compute the loglikelihood over every sample
#     #         return loglikelihood(l, y_j, f)
#     #     end
#     # end / i.nMC
# end

# function expec_loglikelihood_datapoint(l::AbstractLikelihood, raw_samples, y, μ, σ², j, N)
#     samples = raw_samples .* sqrt.([σ²[k][j] for k in 1:N])' .+ [μ[k][j] for k in 1:N]'
#     # samples is of dimension nMC x nLatent again
#     y_j = getindex.(y, j) # Obtain the label for data point j
#     return sum((f->loglikelihood(l, y_j, f)).(eachrow(samples))) # We now compute the loglikelihood over every sample
#     # return
# # end
# end


text:
using WAVEFORM, LinearAlgebra, Test

# Test gmres
n = 100
A = randn(n,n)
A = A'*A
x = randn(100)
b = A*x

# Solve
x0 = similar(x) * 0.0
x0, rec = FGMRES(A, b, x0; tol=1e-9)

@test norm(x - x0) / norm(x0) < 1e-4

# Solve w/ preconditioner
Ai = A^-1   # perfect preconditioner
precond(x) = Ai*x

x0 = similar(x) * 0.0
x0, rec = FGMRES(A, b, x0; tol=1e-9, precond=precond)

@test norm(x - x0) / norm(x0) < 1e-4

# Test linear solver wrapper: gmres
opts = LinSolveOpts(; tol=1e-9, maxit=1000, solver=:fgmres, precond=:identity)
solver = solvesystem(A, opts)
x0 = randn(n)
x0 = solver(b, x0, true)

@test norm(x - x0) / norm(x0) < 1e-4

# Test linear solver wrapper: lufact
opts = LinSolveOpts(; tol=1e-9, maxit=1000, solver=:lufact, precond=:identity)
solver = solvesystem(A, opts)
x0 = randn(n)
x0 = solver(b, x0, true)

@test norm(x - x0) / norm(x0) < 1e-4


text:
@testset "Unit" begin
    @info "Running unit tests..."
    @test true
end


text:
# linear phase FIR filter design that optimizes maximum error
# in the frequency domain

using DSP

"""
    remez_wrap(numtaps::Integer,
          bands::Vector,
          desired::Vector;
          weight::Vector=[],
          Hz::Real=1.0,
          filter_type::RemezFilterType=filter_type_bandpass,
          maxiter::Integer=25,
          grid_density::Integer=16)

Calculate the minimax optimal filter using the Remez exchange algorithm [^McClellan1973a] [^McClellan1973b].

Calculate the filter-coefficients for the finite impulse response
(FIR) filter whose transfer function minimizes the maximum error
between the desired gain and the realized gain in the specified
frequency bands using the Remez exchange algorithm.

`remez_wrap` is equivalent to `remez` in the DSP package, except it uses
the C code ported from scipy.

# Arguments
- `numtaps::Integer`: The desired number of taps in the filter.
    The number of taps is the number of terms in the filter, or the filter
    order plus one.
- `bands::Vector`: A monotonic sequence containing the band edges in Hz.
    All elements must be non-negative and less than half the sampling
    frequency as given by `Hz`.
- `desired::Vector`:A sequence half the size of bands containing the desired
    gain in each of the specified bands.
- `weight::Vector`: (optional)
    A relative weighting to give to each band region. The length of
    `weight` has to be half the length of `bands`.
- `Hz::Real`: The sampling frequency in Hz. Default is 1.
- `filter_type::RemezFilterType`: Default is filter_type_bandpass.
    The type of filter:
      filter_type_bandpass : flat response in bands. This is the default.
      filter_type_differentiator : frequency proportional response in bands.
        Assymetric as in filter_type_hilbert case, but with a linear sloping
        desired response.
      filter_type_hilbert : filter with odd symmetry, that is, type III
                  (for even order) or type IV (for odd order)
                  linear phase filters.
- `maxiter::Integer`: (optional)
    Maximum number of iterations of the algorithm. Default is 25.
- `grid_density:Integer`: (optional)
    Grid density. The dense grid used in `remez` is of size
    ``(numtaps + 1) * grid_density``. Default is 16.

# Returns
- `h::Array{Float64,1}`: A rank-1 array containing the coefficients of the optimal
    (in a minimax sense) filter.

[^McClellan1973a]:
J. H. McClellan and T. W. Parks, A unified approach to the
design of optimum FIR linear phase digital filters,
IEEE Trans. Circuit Theory, vol. CT-20, pp. 697-701, 1973.

[^McClellan1973b]:
J. H. McClellan, T. W. Parks and L. R. Rabiner, A Computer
Program for Designing Optimum FIR Linear Phase Digital
Filters, IEEE Trans. Audio Electroacoust., vol. AU-21,
pp. 506-525, 1973.

# Examples
Construct a length 35 filter with a passband at 0.15-0.4 Hz
(desired response of 1), and stop bands at 0-0.1 Hz and 0.45-0.5 Hz
(desired response of 0). Note: the behavior in the frequency ranges between
those bands - the transition bands - is unspecified.

```julia-repl
julia> bpass = remez_wrap(35, [0 0.1 0.15 0.4 0.45 0.5], [0 1 0])
```

You can trade-off maximum error achieved for transition bandwidth.
The wider the transition bands, the lower the maximum error in the
bands specified. Here is a bandpass filter with the same passband, but
wider transition bands.

```julia-repl
julia> bpass2 = remez_wrap(35, [0 0.08 0.15 0.4 0.47 0.5], [0 1 0])
```

Here we compute the frequency responses and plot them in dB.

```julia-repl
using PyPlot
b = DSP.Filters.PolynomialRatio(bpass, [1.0])
b2 = DSP.Filters.PolynomialRatio(bpass2, [1.0])
f = linspace(0, 0.5, 1000)
plot(f, 20*log10.(abs.(freqz(b,f,1.0))))
plot(f, 20*log10.(abs.(freqz(b2,f,1.0))))
grid()
```
"""
function remez_wrap(numtaps::Integer, bands::Vector, desired::Vector;
               weight::Vector=[],
               Hz::Real=1.0,
               filter_type::RemezFilterType=filter_type_bandpass,
               maxiter::Integer=25,
               grid_density::Integer=16)

    # Convert type
    if (filter_type == filter_type_bandpass)
        tnum = 1
    elseif (filter_type == filter_type_differentiator)
        tnum = 2
    elseif (filter_type == filter_type_hilbert)
        tnum = 3
    end
    if (length(weight)==0)
        weight = ones(desired)
    end

    h = zeros(numtaps)
    #return sigtools._remez(numtaps, bands, desired, weight, tnum, Hz,
                           #maxiter, grid_density)
    #int pre_remez(double *h2, int numtaps, int numbands, double *bands,
    #              double *response, double *weight, int type, int maxiter,
    #              int grid_density);
    bands = copy(bands)/Hz
    ccall((:pre_remez, "remez"), Integer, (Ptr{Cdouble}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cint, Cint),
        h, numtaps, length(bands)/2, bands, desired, weight, tnum, maxiter, grid_density)
    return h

end


text:
using ArgParse
global args

"""
Parses the arguments given by the user
"""
function parse_commandline()

    settings = ArgParseSettings()

    # define arguments (type, description, default value)
    @add_arg_table! settings begin
        "--external_solver", "-x"
            help = "External solver: \"cbc\", \"cplex\" or \"glpk\""
            arg_type = String
            default = "cplex"

        "--input", "-i"
            help = "Instance file path"
            arg_type = Union{Nothing,String}
            default = "data/n20_3.txt"

        "--method", "-m"
            help = "Solution method: \"sequential\", \"single-flow\", \"std-constr-gen\" or \"adv-constr-gen\""
            arg_type = String
            default = "single-flow"

        "--verbosity", "-v"
            help = "Verbosity level"
            arg_type = Int
            default = 1
    end

    parsed_args = parse_args(settings)

    # check validity of arguments
    if !(parsed_args["external_solver"] in ["cbc", "cplex", "glpk"])
        error("Unexpected value for --external_solver: $(parsed_args["external_solver"])")
    elseif !(parsed_args["input"] === nothing) && !isfile(parsed_args["input"])
        error("Invalid path for --input: $(parsed_args["infile"])")
    elseif !(parsed_args["method"] in ["sequential", "single-flow", "std-constr-gen", "adv-constr-gen"])
        error("Unexpected value for --method: $(parsed_args["method"])")
    end

    return parsed_args

end

"""
Outputs the active parameters
"""
function display_active_args()
    ln0("\n========= Active parameters =========\n")
    global args
    for arg in args
        if arg[1] in ["method", "input"]
            ln0("\t$(arg[1])\t\t= $(arg[2])")
        else
            ln0("\t$(arg[1])\t= $(arg[2])")
        end
    end
end

max_stars_repo_name: antcrepin/flight-route-optimizer

max_issues_repo_name: antcrepin/flight-route-optimizer

max_forks_repo_name: antcrepin/flight-route-optimizer


text:
using Compat # until https://github.com/JuliaCI/BenchmarkTools.jl/issues/83 is fixed
using BenchmarkTools
using ByteParsers
using TextParse
import TextParse: tryparsenext

sep = UInt8(';')
lopt = TextParse.LocalOpts(';', '\"', '\\', false, false)

string_and_bytes(s) = s, convert(Vector{UInt8}, s)

# parsing integers
s1, b1 = string_and_bytes("123456789;")
p1 = PositiveInteger()
t1 = TextParse.Numeric(Int)

@benchmark parsenext($p1, $b1, 1, $sep)                    # 15 ns
@benchmark tryparsenext($t1, $s1, 1, $(length(s1)), $lopt) # 21 ns

# parsing dates
s2, b2 = string_and_bytes("19800101;")
p2 = DateYYYYMMDD()
t2 = TextParse.DateTimeToken(Date, DateFormat("yyyymmdd"))

@benchmark parsenext($p2, $b2, 1, $sep)                    # 29 ns
@benchmark tryparsenext($t2, $s2, 1, $(length(s2)), $lopt) # 65 ns

# parsing strings
s3, b3 = string_and_bytes("foobar;")
p3 = ViewBytes()
t3 = TextParse.StringToken(TextParse.StrRange)

@benchmark parsenext($p3, $b3, 1, $sep)                    # 19 ns
@benchmark tryparsenext($t3, $s3, 1, $(length(s3)), $lopt) # 25 ns

# skipping string (use previous b3)

p4 = Skip()
@benchmark parsenext($p4, $b3, 1, $sep)                    # 6 ns

######################################################################
# a mess of various benchmarks below, FIXME: organize
######################################################################

str = b"1980;1990;19800101;"
sep = UInt8(';')

p1 = PositiveInteger(Int)
p2 = PositiveFixedInteger(4)
pd = DateYYYYMMDD()
r1 = Line(p1, p1)
r2 = Line(p1, p1, pd)

parsenext(p1, str, 1, sep)
parsenext(p2, str, 1, sep)
parsenext(r1, str, 1, sep)
parsenext(r2, str, 1, sep)
@code_typed parsenext(r1, str, 1, sep)

@benchmark parsenext($p1, $str, 1, $sep) # 12 ns
@benchmark parsenext($r1, $str, 1, $sep) # 19 ns
@benchmark parsenext($p2, $str, 1, $sep) # 11 ns
@benchmark parsenext($r2, $str, 1, $sep) # 44 ns
@benchmark parsenext($pd, $str, 11, $sep) # 30 ns

str2 = b"irrelevant;"
ps = Skip()
pv = ViewBytes()
@benchmark parsenext($ps, $str2, 1, $sep) # 8 ns
@benchmark parsenext($pv, $str2, 1, $sep) # 20 ns

## skipping
skiptest = Line(PositiveInteger(), Skip(), Skip(), Skip(), ViewBytes(), ViewBytes())
skipstr = b"1234567;skipped;skipped;skipped;relevant;super-relevant;"
parsenext(skiptest, skipstr, 1, sep)
@benchmark parsenext($skiptest, $skipstr, 1, $sep) # 80ns

## how does Date parsing compare?
df = DateFormat("yyyymmdd")
str2 = "19800101"
@benchmark tryparse(Date, $str2, $df) # 69 ns

parsenext(PositiveInteger(), b"111;", 1, sep)
parsenext(PositiveInteger(), b"111", 1, sep)


text:
abstract type AbstractRiskRegion end

@deprecate RiskRegion EllipticalRiskRegion

mutable struct EllipticalRiskRegion <: AbstractRiskRegion
    μ::Vector{Float64}
    K::Cone
    inv_P::Matrix{Float64}
    α::Float64
    function EllipticalRiskRegion(μ::Vector{Float64}, Σ::Matrix{Float64},
                        K::FiniteCone, α::Float64)
        checkRiskRegionArgs(μ, Σ, K, α)
        P = LinearAlgebra.cholesky(Σ).U
        inv_P = inv(P)
        new_K = FiniteCone(P*(-K.A))
        new(μ, new_K, inv_P, α)
    end

    function EllipticalRiskRegion(μ::Vector{Float64}, Σ::Matrix{Float64},
                        K::PolyhedralCone, α::Float64)
        checkRiskRegionArgs(μ, Σ, K, α)
        P = LinearAlgebra.cholesky(Σ).U
        inv_P = inv(P)
        new_K = PolyhedralCone{Float64}((-K.A) * inv_P) # quickfix
        new(μ, new_K, inv_P, α)
    end
end

# function EllipticalRiskRegion(μ::Distributions.ZeroVector, Σ::Matrix{Float64}, K::Cone, α::Float64)
#     EllipticalRiskRegion(zeros(μ.len), Σ, K, α)
# end
EllipticalRiskRegion(dist::AbstractMvNormal, K::Cone, β::Float64) = EllipticalRiskRegion(dist.μ, dist.Σ.mat, K, quantile(Normal(), β))
EllipticalRiskRegion(dist::Distributions.AbstractMvTDist, K::Cone, β::Float64) = EllipticalRiskRegion(dist.μ, dist.Σ.mat, K, quantile(TDist(dist.df), β))
# EllipticalRiskRegion(dist::MvSkewTDist, K::Cone, β::Float64) = EllipticalRiskRegion(dist.ξ, dist.Ω.mat, K, quantile(TDist(dist.df), β))

function transformed_size(Ω::EllipticalRiskRegion, x::Vector{Float64})
    x_trans = x - Ω.μ
    y = Ω.inv_P'x_trans
    norm(project(Ω.K, y))
end

function ∈(x::Vector{Float64}, Ω::EllipticalRiskRegion)
    transformed_size(Ω::EllipticalRiskRegion, x::Vector{Float64}) >= Ω.α
end

length(Ω::EllipticalRiskRegion) = length(Ω.μ)

mutable struct MonotonicEllipticalRiskRegion <: AbstractRiskRegion
    # Note that cost function for portfolio selection is negative
    # monotonic, so roles of risk and non-risk frontiers are reversed
    Ω::EllipticalRiskRegion
    nonrisk_frontier::LinkedList{VecF64}
    risk_frontier::LinkedList{VecF64}
    max_frontier_points::Int
    function MonotonicEllipticalRiskRegion(args...)
        Ω=EllipticalRiskRegion(args...)
        new(Ω, nil(VecF64), nil(VecF64), length(Ω)*10)
    end
end

function ∈(x::VecF64, Ω::MonotonicEllipticalRiskRegion)
    # Note that cost function for portfolio selection is negative
    # monotonic, so roles of risk and non-risk frontiers are reversed
    if above_risk_frontier(x, Ω.risk_frontier) return false
    elseif below_nonrisk_frontier(x, Ω.nonrisk_frontier) return true
    else
        if x ∈ Ω.Ω
            if length(Ω.nonrisk_frontier) < Ω.max_frontier_points
                add_to_nonrisk_frontier!(Ω, x)
            end
            return true
        else
            if length(Ω.risk_frontier) < Ω.max_frontier_points
                add_to_risk_frontier!(Ω, x)
            end
            return false
        end
    end
end


text:
module _rMSEn
export rMSEn
using Statistics: std, mean, median, var
using DSP.Filters: filtfilt, Butterworth, Lowpass, digitalfilter
using Plots

    """
        MSx, CI = rMSEn(Sig, Mobj)

    Returns a vector of refined multiscale entropy values (`MSx`) and the complexity
    index (`CI`) of the data sequence (`Sig`) using the parameters specified by
    the multiscale object (`Mobj`) and the following default parameters:
    Scales = 3, Butterworth LPF Order = 6, Butterworth LPF cutoff frequency
    at scale (T): Fc = 0.5/T.
    If the entropy function specified by `Mobj` is SampEn or ApEn, rMSEn
    updates the threshold radius of the data sequence (Xt) at each scale
    to 0.2std(Xt) if no `r` value is provided by Mobj, or r.std(Xt) if `r`
    is specified.

        MSx, CI = rMSEn(Sig::AbstractArray{T,1} where T<:Real, Mobj::NamedTuple; Scales::Int=3,
                            F_Order::Int=6, F_Num::Float64=0.5, RadNew::Int=0, Plotx::Bool=false)

    Returns a vector of refined multiscale entropy values (`MSx`) and the complexity
    index (`CI`) of the data sequence (`Sig`) using the parameters specified by
    the multiscale object (`Mobj`) and the following 'keyword' arguments:

    # Arguments:
    `Scales`   - Number of temporal scales, an integer > 1 (default = 3) \n
    `F_Order`  - Butterworth low-pass filter order, a positive integer
                 (default: 6) \n
    `F_Num`    - Numerator of Butterworth low-pass filter cutoff frequency,
                 a scalar value in range [0 < `F_Num` < 1]. The cutoff frequency
                 at each scale (T) becomes: Fc = `F_Num/T.  (default: 0.5) \n
    `RadNew`   - Radius rescaling method, an integer in the range [1 4].
                 When the entropy specified by `Mobj` is `SampEn` or `ApEn`,
                 `RadNew` allows the radius threshold to be updated at each
                 time scale (Xt). If a radius value is specified by `Mobj` (`r`),
                 this becomes the rescaling coefficient, otherwise it is set
                 to 0.2 (default). The value of `RadNew` specifies one of the
                 following methods:\n
                 [1]    Standard Deviation          - r*std(Xt)\n
                 [2]    Variance                    - r*var(Xt) \n
                 [3]    Mean Absolute Deviation     - r*mean_ad(Xt) \n
                 [4]    Median Absolute Deviation   - r*med_ad(Xt)\n
    `Plotx`    - When `Plotx` == true, returns a plot of the entropy value at each
                 time scale (i.e. the multiscale entropy curve)  [default: false] \n

    # See also `MSobject`, `MSEn`, `cMSEn`, `hMSEn`, `SampEn`, `ApEn`, `XMSEn`

    # References:
        [1] Madalena Costa, Ary Goldberger, and C-K. Peng,
            "Multiscale entropy analysis of complex physiologic time series."
            Physical review letters
            89.6 (2002): 068102.

        [2] Vadim V. Nikulin, and Tom Brismar,
            "Comment on “Multiscale entropy analysis of complex physiologic
            time series”."
            Physical review letters
            92.8 (2004): 089803.

        [3] Madalena Costa, Ary L. Goldberger, and C-K. Peng.
            "Costa, Goldberger, and Peng reply."
            Physical Review Letters
            92.8 (2004): 089804.

        [4] José Fernando Valencia, et al.,
            "Refined multiscale entropy: Application to 24-h holter
            recordings of heart period variability in healthy and aortic
            stenosis subjects."
            IEEE Transactions on Biomedical Engineering
            56.9 (2009): 2202-2213.

        [5] Puneeta Marwaha and Ramesh Kumar Sunkaria,
            "Optimal selection of threshold value ‘r’for refined multiscale
            entropy."
            Cardiovascular engineering and technology
            6.4 (2015): 557-576.

    """
    function rMSEn(Sig::AbstractArray{T,1} where T<:Real, Mobj::NamedTuple; Scales::Int=3,
        F_Order::Int=6, F_Num::Float64=0.5, RadNew::Int=0, Plotx::Bool=false)

    (size(Sig,1)>10) ? nothing : error("Sig:   must be a numeric vector" )
    (length(Mobj) >= 1) ? nothing :  error("Mobj:    must be a multiscale entropy object created
         with the function EntropyHub.MSobject")
    (Scales>1) ? nothing : error("Scales:     must be an integer > 1")
    (F_Order>1) ? nothing :  error("F_Order:     must be an integer > 1")
    (0 < F_Num < 1) ? nothing : error("F_Num:     must be a scalar in range 0 < F_Num < 1")
    (RadNew==0 || (RadNew in 1:4 && String(Symbol(Mobj.Func)) in ("SampEn","ApEn"))) ? nothing :
    error("RadNew:  must be 0, or an integer in range [1 4] with entropy function `SampEn` or `ApEn`")

    MSx = zeros(Scales)
    Args = NamedTuple{keys(Mobj)[2:end]}(Mobj)
    (RadNew==0 && String(Symbol(Mobj.Func)) in ("SampEn","ApEn")) ? RadNew=1 : nothing
    if RadNew > 0
        if RadNew == 1
            Rnew = x -> std(x, corrected=false)
        elseif RadNew == 2
            Rnew = x -> var(x, corrected=false)
        elseif RadNew == 3
            Rnew = x -> mean(abs.(x .- mean(x)))
        elseif RadNew == 4
            Rnew = x -> median(abs.(x .- median(x)))
        end

        if haskey(Mobj,:r)
            Cx = Mobj.r
        else
            Cy = ("Standard Deviation","Variance","Mean Abs Deviation", "Median Abs Deviation")
            @warn("No radius value provided in Mobj.
                Default set to 0.2*$(Cy[RadNew]) of each new time-series.")
            Cx = .2
        end
    end

    for T = 1:Scales
        print(". ")
        Temp = refined(Sig,T,F_Order,F_Num)
        RadNew > 0 ? Args = (Args..., r=Cx*Rnew(Temp[:])) : nothing
        Tempx = Mobj.Func(Temp[:]; Args...)
        #MSx[T] = Tempx[1][end]
        typeof(Tempx)<:Tuple ? MSx[T] = Tempx[1][end] : MSx[T] = Tempx[end]
    end
    CI = sum(MSx)
    print("\n")
    if any(isnan.(MSx))
        println("Some entropy values may be undefined.")
    end

    if Plotx
        p1 = plot(1:Scales, MSx, c=RGB(8/255, 63/255, 77/255), lw=3)
        scatter!(1:Scales, MSx, markersize=6, c=RGB(1, 0, 1),
        xlabel = "Scale Factor", ylabel = "Entropy Value",
        guidefont = font(12, "arial", RGB(7/255, 54/255, 66/255)),
        tickfontsize = 10, tickfontfamily="arial", legend=false,
        title = "Refined Multiscale $(Mobj.Func)",
        plot_titlefontsize=16, plot_titlefontcolor=RGB(7/255, 54/255, 66/255))
        display(p1)
    end

    return MSx, CI
    end

    function refined(Z,sx,P1,P2)
        Yt = filtfilt(digitalfilter(Lowpass(P2/sx), Butterworth(P1)), Z)
        return Yt[1:sx:end]
    end

end

"""

Copyright 2021 Matthew W. Flood, EntropyHub

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

For Terms of Use see https://github.com/MattWillFlood/EntropyHub

"""


text:
module TestAMLPipeline
using Test
include("test_featureselector.jl")
include("test_util.jl")
include("test_baselinemodel.jl")
include("test_basefilter.jl")
include("test_decisiontree.jl")
include("test_ensemble.jl")
include("test_crossvalidator.jl")
include("test_naremover.jl")
include("test_pipeline.jl")
end


text:
#=
  From https://probmods.org/chapters/conditioning.html
  (WebPPL):
  """
  This classic Bayesian inference task is a special case of conditioning. Kahneman and Tversky,
  and Gigerenzer and colleagues, have studied how people make simple judgments like the following:

  The probability of breast cancer is 1% for a woman at 40 who participates in a routine screening.
  If a woman has breast cancer, the probability is 80% that she will have a positive mammography.
  If a woman does not have breast cancer, the probability is 9.6% that she will also have a
  positive mammography. A woman in this age group had a positive mammography in a routine
  screening. What is the probability that she actually has breast cancer?

  What is your intuition? Many people without training in statistical inference judge the
  probability to be rather high, typically between 0.7 and 0.9. The correct answer is much lower,
  less than 0.1, as we can see by running this WebPPL inference:
  """

  p(breastCancer): 0.0737

  According to the WebPPL model:
    breastCancer:
    false : 0.922360248447205
    true : 0.07763975155279507

=#

using Turing
include("jl_utils.jl")

@model function medical_diagnosis2()
    breastCancer ~ flip(0.01)
    positiveMammogram ~ breastCancer ? flip(0.8) : flip(0.096)

    true ~ Dirac(positiveMammogram)

end

model = medical_diagnosis2()

chns = sample(model, MH(), 10_000)
# chns = sample(model, PG(15), 10_000)
# chns = sample(model, IS(), 10_000)
# chns = sample(model, SMC(), 10_000)

display(chns)


text:
# ATTENTION: This file was generated, do not edit by hand.
module GAP_pkg_ferret

using Pkg.Artifacts
using GAP

include("jll.jl")

function __init__()
    # ensure GAP will load the correct instance of this package
    path = artifact"ferret"
    @debug "GAP package 'ferret' path = " * path
    GAP.Globals.SetPackagePath(GapObj("ferret"), GapObj(path))
end

function load()
    res = GAP.Globals.LoadPackage(GapObj("ferret"))
    if res != true
        error("failed to load GAP package ferret")
    end
end
end # module