text:
using GraphRecipes, Impero, Test
using Plots
theme(:default)
default(size=(500, 500))
plot_flag = false

@wrapper α=1 β=1 ψ=1 ϕ=1

@testset "Testing to_expr function" begin
    impero_expr = α + β + sin(ϕ)^2
    julia_expr = to_expr(impero_expr)
    @test impero_expr == eval(julia_expr)
end

if plot_flag
    impero_expr = α + β + sin(ϕ)
    julia_expr = to_expr(impero_expr)
    plot(julia_expr)
end


text:
function getTableauVSPARK(s, σ, o, tsym, g, h, lq, lp, ω, d=Nothing)
    # α_q_1 and α_p_2 need to be conjugate symplectic
    # α_q_1 and α_p_3 need to be conjugate symplectic
    # α_q_2 and α_p_2 need to be conjugate symplectic
    # α_q_2 and α_p_3 need to be conjugate symplectic
    # -> α_q_1 and α_q_2 need to be identical
    # -> α_p_2 and α_p_3 need to be identical

    α_q_1 = lq.a
    α_q_2 = lq.a
    α_p_2 = lp.a
    α_p_3 = lp.a

    b_q_1 = lq.b
    b_q_2 = lq.b
    b_p_2 = lp.b
    b_p_3 = lp.b

    # α_p_1 is free, but determines a_q_1 and a_q_2
    α_p_1 = h.a

    # a_p_1, a_p_2 and a_p_3 are free (they are not used in VSPARKpSecondary integrators anymore)
    a_p_1 = g.a
    b_p_1 = g.b

    a_p_2 = g.a
    a_p_3 = g.a

    β_q_1 = lq.b
    β_q_2 = lq.b

    β_p_1 = lp.b
    β_p_2 = lp.b
    β_p_3 = lp.b

    # a_q_1_ij = b_q_1_j / b_p_1_i * (b_p_1_i - α_p_1_ji)
    a_q_1 = zeros(s,σ)
    for i in 1:s
        for j in 1:σ
            a_q_1[i,j] = b_q_1[j] / b_p_1[i] * (b_p_1[i] - α_p_1[j,i])
        end
    end

    # a_q_2_ij = b_q_2_j / b_p_1_i * (b_p_1_i - α_p_1_ji)
    a_q_2 = zeros(s,σ)
    for i in 1:s
        for j in 1:σ
            a_q_2[i,j] = b_q_2[j] / b_p_1[i] * (b_p_1[i] - α_p_1[j,i])
        end
    end

    a_q = (a_q_1, a_q_2)
    b_q = (b_q_1, b_q_2)
    c_q = g.c

    a_p = (a_p_1, a_p_2, a_p_3)
    b_p = (b_p_1, b_p_2, b_p_3)
    c_p = g.c

    α_q = (α_q_1, α_q_2)
    β_q = (β_q_1, β_q_2)
    γ_q = lq.c

    α_p = (α_p_1, α_p_2, α_p_3)
    β_p = (β_p_1, β_p_2, β_p_3)
    γ_p = lp.c

    coeff_q = CoefficientsSPARK(tsym, o, s, σ, a_q, b_q, c_q)
    coeff_p = CoefficientsSPARK(tsym, o, s, σ, a_p, b_p, c_p)
    coeff_q̃ = CoefficientsSPARK(tsym, o, σ, s, α_q, β_q, γ_q)
    coeff_p̃ = CoefficientsSPARK(tsym, o, σ, s, α_p, β_p, γ_p)

    TableauVSPARKsecondary(tsym, o, s, σ, coeff_q, coeff_p, coeff_q̃, coeff_p̃, ω, d)
end

function TableauVSPARKLobattoIII(s, lq, lp; name = Symbol("VSPARKLobattoIII"))
    o = 2s-2
    getTableauVSPARK(s, s, o, name, lq, lp, lq, lp, get_lobatto_ω_matrix(s), get_lobatto_nullvector(s))
end

function TableauVSPARKLobattoIIIAB(s)
    lq = TableauLobattoIIIA(s)
    lp = TableauLobattoIIIB(s)
    TableauVSPARKLobattoIII(s, lq, lp; name = Symbol("VSPARKLobattoIIIAIIIB"))
end

function TableauVSPARKLobattoIIIBA(s)
    lq = TableauLobattoIIIB(s)
    lp = TableauLobattoIIIA(s)
    TableauVSPARKLobattoIII(s, lq, lp; name = Symbol("VSPARKLobattoIIIBIIIA"))
end

function TableauVSPARKLobattoIIICC̄(s)
    lq = TableauLobattoIIIC(s)
    lp = TableauLobattoIIIC̄(s)
    TableauVSPARKLobattoIII(s, lq, lp; name = Symbol("VSPARKLobattoIIICIIIC̄"))
end

function TableauVSPARKLobattoIIIC̄C(s)
    lq = TableauLobattoIIIC̄(s)
    lp = TableauLobattoIIIC(s)
    TableauVSPARKLobattoIII(s, lq, lp; name = Symbol("VSPARKLobattoIIIC̄IIIC"))
end

function TableauVSPARKLobattoIIID(s)
    l = TableauLobattoIIID(s)
    TableauVSPARKLobattoIII(s, l, l; name = Symbol("VSPARKLobattoIIID"))
end

function TableauVSPARKLobattoIIIE(s)
    l = TableauLobattoIIIE(s)
    TableauVSPARKLobattoIII(s, l, l; name = Symbol("VSPARKLobattoIIIE"))
end

function TableauVSPARKGLRKLobattoIII(s, σ, lq, lp; name = Symbol("VSPARKGLRKLobattoIII"))
    o = 2s
    g = TableauGauss(s)
    getTableauVSPARK(s, σ, o, name, g, get_lobatto_glrk_coefficients(s, σ), lq, lp, get_GLRK_ω_matrix(σ), get_lobatto_nullvector(σ))
end

function TableauVSPARKGLRKLobattoIIIAB(s, σ=s+1)
    TableauVSPARKGLRKLobattoIII(s, σ, TableauLobattoIIIA(σ), TableauLobattoIIIB(σ); name = Symbol("VSPARKGLRKLobattoIIIAIIIB"))
end

function TableauVSPARKGLRKLobattoIIIBA(s, σ=s+1)
    TableauVSPARKGLRKLobattoIII(s, σ, TableauLobattoIIIB(σ), TableauLobattoIIIA(σ); name = Symbol("VSPARKGLRKLobattoIIIBIIIA"))
end

function TableauVSPARKGLRKLobattoIIICC̄(s, σ=s+1)
    TableauVSPARKGLRKLobattoIII(s, σ, TableauLobattoIIIC(σ), TableauLobattoIIIC̄(σ); name = Symbol("VSPARKGLRKLobattoIIICIIIC̄"))
end

function TableauVSPARKGLRKLobattoIIIC̄C(s, σ=s+1)
    TableauVSPARKGLRKLobattoIII(s, σ, TableauLobattoIIIC̄(σ), TableauLobattoIIIC(σ); name = Symbol("VSPARKGLRKLobattoIIIC̄IIIC"))
end

function TableauVSPARKGLRKLobattoIIID(s, σ=s+1)
    TableauVSPARKGLRKLobattoIII(s, σ, TableauLobattoIIID(σ), TableauLobattoIIID(σ); name = Symbol("VSPARKGLRKLobattoIIID"))
end

function TableauVSPARKGLRKLobattoIIIE(s, σ=s+1)
    TableauVSPARKGLRKLobattoIII(s, σ, TableauLobattoIIIE(σ), TableauLobattoIIIE(σ); name = Symbol("VSPARKGLRKLobattoIIIE"))
end

ext: jl


text:
import StatsBase:rmsd

@sk_import linear_model: LinearRegression

"""
    get_molecules(data::Dict)

Returns molecules from data file.
This is usability function that is meant to simplify things.
"""
function get_molecules(data::Dict)
    m1=MoleculeIdenticalInformation{AtomOnlySymbol}(data["cluster1"].atoms)
    m2=MoleculeIdenticalInformation{AtomOnlySymbol}(data["cluster2"].atoms)
    return m1, m2
end

"""
    FitData

Structure to help potential parameters fitting.

# Fields
- `variables` : variables
- `E`         : energy
- `w`         : weights

# Creation
    FitData(mpp, points, energy)
    FitData(mpp, data::Dict...)

- `data::Dict`   :  Dict that is returned by `load_data_file`

"""
mutable struct FitData
    "Variables"
    variables::Vector{Any}
    "Energy"
    E::Vector{Float64}
    "Weights"
    w::Vector{Float64}
    function FitData(mpp::MoleculePairPotential, points, energy)
        new( potential_variables(mpp,points) , vec(energy), ones(length(energy)))
    end
end

function FitData(mpp::MoleculePairPotential, data::Dict...)
    function _samemolecule(c1, c2)
        if length(c1) != length(c2)
            return false
        elseif all(c1.atoms .== c2.atoms)
            return true
        else
            return false
        end
    end
    for x in data
        @assert haskey(x, "Energy")
        @assert haskey(x, "Points")
    end
    @assert all( x-> _samemolecule(x["cluster1"], data[1]["cluster1"]), data)
    @assert all( x-> _samemolecule(x["cluster2"], data[1]["cluster2"]), data)

    ptmp = [x["Points"] for x in data]
    points = hcat(ptmp...)
    etmp = [x["Energy"] for x in data]
    energy = hcat(etmp...)

    return FitData(mpp, points, energy)
end

"""
give_as_potential(T, data)

Returns MoleculePairPotential{T} from given data that must have `"c1_molecule"`
and `"c2_molecule"` fields
"""
function give_as_potential(T, data)
    return MoleculePairPotential{T}(data["c1_molecule"],data["c2_molecule"])
end

"""
setweight_e_more!(data::FitData, w, e; unit="cm-1")

Sets weigth when energy is more than given one.

# Arguments
- `data::FitData`  : data where weigth is adjusted
- `w`              : new weigth
- `e`              : energy
- `unit="cm-1"`    : energy unit
"""
function setweight_e_more!(data::FitData, w, e; unit="cm-1")
    ec = energy_from(e, unit)
    data.w[data.E .> ec] .= w
end

"""
setweight_e_more!(data::FitData, w, e; unit="cm-1")

Sets weigth when energy is more than given one.

# Arguments
- `data::FitData`  : data where weigth is adjusted
- `w`              : new weigth
- `e`              : energy
- `unit="cm-1"`    : energy unit
"""
function setweight_e_less!(data::FitData, w, e, unit="cm-1")
    ec = energy_from(e, unit)
    data.w[data.E .< ec] .= w
end

"""
fit_potential!(model, mpp::MoleculePairPotential, fdata::FitData)

Fits potential using given model

# Arguments
- `model`                       : ScikitLearn model
- `mpp::MoleculePairPotential`  : potential
- `fdata::FitData`              : data used in fitting
"""
function fit_potential!(model, mpp::MoleculePairPotential, fdata::FitData)
    r = hcat(fdata.variables...)
    fit!(model, r, fdata.E, fdata.w)
    tl =  [ size(v)[2]  for v in fdata.variables ]
    ir=[]
    i = 1
    for x in tl
        push!(ir, i:i+x-1)
        i += x
    end
    for i in eachindex(ir)
        get_potential!(mpp.topology[i].potential, model.coef_[ir[i]]...)
    end
    mpp
end

"""
predict_potential(model, mpp::MoleculePairPotential, points)

Uses `model` to predict potential on given points
"""
function predict_potential(model, mpp::MoleculePairPotential, points)
    l1 = length(mpp.mol1)
    l2 = length(mpp.mol2)
    c1 = map(x->x[1:l1],  points)
    c2 = map(x->x[l1+1:l1+l2],  points)
    return predict_potential(moldel, mpp, c1, c2)
end

"""
predict_potential(model, mpp::MoleculePairPotential, cluster1, cluster2)

Uses `model` to predict potential on given cluster points
"""
function predict_potential(model, mpp::MoleculePairPotential, cluster1, cluster2)
    r = hcat(potential_variables(mpp,cluster1,cluster2)...)
    return predict(model, r)
end

"""
rmsd(points, energy, mpp::MoleculePairPotential; emax=0, unit="cm^-1")

Calculates root mean square error for potential `mpp`.

# Atributes
- `points`                      : points where potential is tested
- `energy`                      : referece energy for given points
- `mpp::MoleculePairPotential`  : potential
- `emax=0`                      : cut points where energy is larger than this
- `unit="cm^-1"`                : unit for energy
"""
function rmsd(points, energy, mpp::MoleculePairPotential; emax=0, unit="cm^-1")
    @assert size(points) == size(energy) "points and energy need to have same size"
    e = energy_from(emax,unit)
    i = energy .< e
    ec = mpp.(points)
    return energy_to(rmsd(ec[i], energy[i]),unit)
end


text:
using WebAssembly, Test
using WebAssembly.Instructions
using WebAssembly: parsewast, interpretwasm, interpret_module, WType, Func, Module, FuncType, Func, Table, Mem, Global, Elem, Data, Import, Export, i32, i64, f32, f64

@testset "WebAssembly" begin

  b = Block([Nop(), Nop()]) |> WebAssembly.nops
  @test isempty(b.body)

  @testset "Interpreter" begin
    include("interpret.jl")
  end

  @testset "IR" begin
    include("ir.jl")
  end

  @testset "Import-export" begin
    m = WebAssembly.Module(FuncType[], Func[], Table[], Mem[], Global[], Elem[], Data[], Ref(0),
               [Import(:env, :mathfun, :func, [i32, f64], f64)],
               [Export(:fun, :fun, :func)])
  end

end


text:
import Pkg
Pkg.activate(".")
Pkg.instationate()
include("test-notebook.jl")

max_stars_repo_name: Jblew/towards-julia-notebook

max_issues_repo_name: Jblew/towards-julia-notebook

max_forks_repo_name: Jblew/towards-julia-notebook


text:
# using Plots
using RecipesBase
# IDEA: use MAKIE for faster compile?
# IDEA: 2d and 3d recipies
# default range
plotting_range(G::Gaussian) = range(-1, 1, length = 256)
function plotting_range(G::Gaussian{<:ConditionalMean})
	range(minimum(G.μ.x) - std(G.μ.x), maximum(G.μ.x) + std(G.μ.x), length = 256)
end

@recipe function f(G::Gaussian, x::AbstractVector = plotting_range(G))
	seriestype  :=  :path
	μ = mean(G, x)
	σ = std(G).(x)
	ribbon := 2σ
	(x, μ)
end

# IDEA: add a few samples to the plot (currently requires a dense factorization)
# if n > 0
# 	s = sample(f, n, tol = 1e-6)
# 	plot!(xs, s, seriesalpha = .7)
# end


text:
### A Pluto.jl notebook ###
# v0.14.6

using Markdown
using InteractiveUtils

# ╔═╡ a4028042-2072-11eb-0aa6-77e29a756648
begin
	using TextAnalysis
	using Arpack
	using Clustering
	using Languages
	using Plots
end

# ╔═╡ 2abb6d30-2079-11eb-2d32-47713a476605
md"
# Clustering Textbook

We will try to cluster 144 sections of textbook _Matematika 1_ - a first semester calculus with some linear algebra and analytic geometry using spectral clustering on a `DocumentTermMatrix()`.

Package `TextAnalysis.jl` has some needed functionality.
"

# ╔═╡ 30498550-2071-11eb-0f6b-33372a37909e
# Uncoment this block to download files. You may need to crete directory
# "../files/Mat1/node$i.html" manually.
#=
for i=1:144
	download("http://www.mathematics.digital/matematika1/predavanja/node$i.html","../files/Mat1/node$i.html")
end
=#

# ╔═╡ 621ed340-2072-11eb-1742-3b4d0b042fb9
# Number of documents / nodes
n=144

# ╔═╡ b59d17a0-2100-11eb-27de-c3196ecc1268
md"
__N.B.__ Added `ii,iii,iv,v,nbsp,times,home` to `Languages/.../data/stopwords/Croatian.txt`. Correct the filename if it is `Croation.txt`.
"

# ╔═╡ 08057a10-20dd-11eb-300e-f52ff0eaf35f
begin
	# Test on one file
	f=open("./files/Mat1/node10.html","r")
	fr=readlines(f)
	close(f)
	frjoin=join(fr,"  ")
	sdf=StringDocument(frjoin)
	language!(sdf,Languages.Croatian())
	remove_corrupt_utf8!(sdf)
	prepare!(sdf, strip_html_tags | strip_case | strip_numbers | strip_whitespace |  strip_punctuation | strip_stopwords)
end

# ╔═╡ 68174819-9f40-4627-ab86-5bda4f6cc652
frjoin

# ╔═╡ 6c615a00-20de-11eb-0841-c30b5cf7d68a
language(sdf)

# ╔═╡ 1c2c2b0e-20dd-11eb-39b6-a7a56b61391c
TextAnalysis.text(sdf)

# ╔═╡ 022dae30-2076-11eb-225e-233a63507457
a=Array{StringDocument{String}}(undef,n)

# ╔═╡ 2f61f232-2076-11eb-0606-351e586f37f1
# Process all files
for i=1:n
	f=open("./files/Mat1/node$i.html","r")
	fr=readlines(f)
	close(f)
	frjoin=join(fr,"  ")
	sdf=StringDocument(frjoin)
	language!(sdf,Languages.Croatian())
	remove_corrupt_utf8!(sdf)
	prepare!(sdf, strip_html_tags | strip_case | strip_numbers | strip_whitespace |  strip_punctuation | strip_stopwords)
	a[i]=sdf
end

# ╔═╡ 0b67ef4e-2077-11eb-2716-5fed149d4e28
c=Corpus(a)

# ╔═╡ 168f01c0-2077-11eb-20ee-e381bbb1aa61
update_lexicon!(c)

# ╔═╡ 24e5a940-2077-11eb-3071-21fb31744527
update_inverse_index!(c)

# ╔═╡ 307d09b2-2077-11eb-2f26-2d00a7bae004
c

# ╔═╡ 3d92a5b0-2077-11eb-39ba-5df49551535e
c.documents

# ╔═╡ 4296d3ae-2077-11eb-0ecb-13e1aff1268c
c.lexicon

# ╔═╡ 4996d1b0-2077-11eb-2ac0-b7296b769be5
c.h

# ╔═╡ 4c8ce030-2077-11eb-0176-9de66dab00b3
c.inverse_index

# ╔═╡ 53a7b932-2077-11eb-2bf4-1feca7dfa435
c.total_terms

# ╔═╡ 58f17bb0-2077-11eb-0939-3dff13742e14
M = DocumentTermMatrix(c)

# ╔═╡ 6ed98bc2-2077-11eb-15fa-3f753c8ed3c7
D=dtm(M)

# ╔═╡ 78d31ab0-2077-11eb-0688-7d70444e8687
T = tf_idf(D)

# ╔═╡ 588e73f4-4544-4683-92da-29fad5da9b45
T[10,:]

# ╔═╡ 88b8c6a0-2077-11eb-3f48-2f453e82c0cb
md"
We are ready for spectral clustering. We try to find 6 clusters.
"

# ╔═╡ 8cb43730-2077-11eb-08ac-451458efac8d
S,rest=svds(T,nsv=6)

# ╔═╡ a70064ae-2077-11eb-2b51-3d76c0e7e933
size(S.U)

# ╔═╡ b34b84be-2077-11eb-3a72-91c2fc5739dd
outU=kmeans(Matrix(transpose(S.U)),6)

# ╔═╡ 3fb27600-2107-11eb-2a03-ed723e052828
scatter(outU.assignments,xlabel="Documents",ylabel="Clusters",legend=false)

# ╔═╡ Cell order:
# ╟─2abb6d30-2079-11eb-2d32-47713a476605
# ╠═a4028042-2072-11eb-0aa6-77e29a756648
# ╠═30498550-2071-11eb-0f6b-33372a37909e
# ╠═621ed340-2072-11eb-1742-3b4d0b042fb9
# ╟─b59d17a0-2100-11eb-27de-c3196ecc1268
# ╠═08057a10-20dd-11eb-300e-f52ff0eaf35f
# ╠═68174819-9f40-4627-ab86-5bda4f6cc652
# ╠═6c615a00-20de-11eb-0841-c30b5cf7d68a
# ╠═1c2c2b0e-20dd-11eb-39b6-a7a56b61391c
# ╠═022dae30-2076-11eb-225e-233a63507457
# ╠═2f61f232-2076-11eb-0606-351e586f37f1
# ╠═0b67ef4e-2077-11eb-2716-5fed149d4e28
# ╠═168f01c0-2077-11eb-20ee-e381bbb1aa61
# ╠═24e5a940-2077-11eb-3071-21fb31744527
# ╠═307d09b2-2077-11eb-2f26-2d00a7bae004
# ╠═3d92a5b0-2077-11eb-39ba-5df49551535e
# ╠═4296d3ae-2077-11eb-0ecb-13e1aff1268c
# ╠═4996d1b0-2077-11eb-2ac0-b7296b769be5
# ╠═4c8ce030-2077-11eb-0176-9de66dab00b3
# ╠═53a7b932-2077-11eb-2bf4-1feca7dfa435
# ╠═58f17bb0-2077-11eb-0939-3dff13742e14
# ╠═588e73f4-4544-4683-92da-29fad5da9b45
# ╠═6ed98bc2-2077-11eb-15fa-3f753c8ed3c7
# ╠═78d31ab0-2077-11eb-0688-7d70444e8687
# ╟─88b8c6a0-2077-11eb-3f48-2f453e82c0cb
# ╠═8cb43730-2077-11eb-08ac-451458efac8d
# ╠═a70064ae-2077-11eb-2b51-3d76c0e7e933
# ╠═b34b84be-2077-11eb-3a72-91c2fc5739dd
# ╠═3fb27600-2107-11eb-2a03-ed723e052828

max_stars_repo_path: src/Module C - Applications/L12a Matematika1.jl

max_stars_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_issues_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_forks_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course


text:
# @variables(model, begin
#     z[i = 1 : ncontacts, j = 1 : nregions, n = 1 : nsteps], Bin
#     r̂[k = 1 : 2, i = 1 : ncontacts, j = 1 : nregions, n = 1 : nsteps]
#     f̂[k = 1 : 2, i = 1 : ncontacts, j = 1 : nregions, n = 1 : nsteps]
# end)
# for j = 1 : nregions
#     region = regions[j]
#     for i = 1 : ncontacts, n = 1 : nsteps
#         polyhedron_constraints(model, region.position, r̂[:, i, j, n])
#         polyhedron_constraints(model, region.force, f̂[:, i, j, n])
#     end
#     lr, ur = axis_aligned_bounding_box(region.position)
#     lf, uf = axis_aligned_bounding_box(region.force)
#     @constraints(model, begin
#         [i = 1 : ncontacts, k = 1 : 2, n = 1 : nsteps], z[i, j, n] * lr[k] <= r̂[k, i, j, n]
#         [i = 1 : ncontacts, k = 1 : 2, n = 1 : nsteps], r̂[k, i, j, n] <= z[i, j, n] * ur[k]
#         [i = 1 : ncontacts, k = 1 : 2, n = 1 : nsteps], z[i, j, n] * lf[k] <= f̂[k, i, j, n]
#         [i = 1 : ncontacts, k = 1 : 2, n = 1 : nsteps], f̂[k, i, j, n] <= z[i, j, n] * uf[k]
#     end)
# end
# @constraints(model, begin
#     [i = 1 : ncontacts, k = 1 : 2, n = 1 : nsteps], sum(r̂[k, i, j, n] for j = 1 : nregions) == rc[k, i, n]
#     [i = 1 : ncontacts, k = 1 : 2, n = 1 : nsteps], sum(f̂[k, i, j, n] for j = 1 : nregions) == f[k, i, n]
#     [i = 1 : ncontacts, n = 1 : nsteps], sum(z[i, j, n] for j = 1 : nregions) == 1
# end)


text:
function test_camera(camera::Raylib.RayCamera2D, ans::Vector)
    return all([
        camera.offset == ans[1],
        camera.target == ans[2],
        camera.rotation == ans[3],
        camera.zoom == ans[4]
    ])
end

function test_camera(camera::Raylib.RayCamera3D, ans::Vector)
    return all([
        camera.position == ans[1],
        camera.target == ans[2],
        camera.up == ans[3],
        camera.fovy == ans[4],
        camera.projection == ans[5]
    ])
end

@testset "set camera properties 2d" begin
    camera = Raylib.RayCamera2D(
        Raylib.rayvector(1, 2), # camera offset
        Raylib.rayvector(3, 4), # camera target
        0, # camera rotation
        1, # camera zoom
    )
    @test test_camera(camera, [
        Raylib.rayvector(1, 2),
        Raylib.rayvector(3, 4),
        0,
        1
    ])

    camera.offset += Raylib.rayvector(5, 6)
    camera.target += Raylib.rayvector(7, 8)
    camera.rotation = 45
    camera.zoom = 3
    @test test_camera(camera, [
        Raylib.rayvector(1, 2) + Raylib.rayvector(5, 6),
        Raylib.rayvector(3, 4) + Raylib.rayvector(7, 8),
        45,
        3
    ])
end

@testset "set camera properties 3d" begin
    camera = Raylib.RayCamera3D(
        Raylib.rayvector(0, 10, 10), # camera position
        Raylib.rayvector(0, 0, 0), # camera looking at point
        Raylib.rayvector(0, 1, 0), # camera up vector (rotation towards target)
        45, # camera field-of-view Y
        Raylib.CAMERA_PERSPECTIVE, # camera mode type
    )
    @test test_camera(camera, [
        Raylib.rayvector(0, 10, 10),
        Raylib.rayvector(0, 0, 0),
        Raylib.rayvector(0, 1, 0),
        45,
        Raylib.CAMERA_PERSPECTIVE
    ])

    camera.position += Raylib.rayvector(1, 2, 3)
    camera.target += Raylib.rayvector(3, 2, 1)
    camera.up += Raylib.rayvector(4, 5, 6)
    camera.fovy = 60
    camera.projection = Raylib.CAMERA_ORTHOGRAPHIC
    @test test_camera(camera, [
        Raylib.rayvector(0, 10, 10) + Raylib.rayvector(1, 2, 3),
        Raylib.rayvector(0, 0, 0) + Raylib.rayvector(3, 2, 1),
        Raylib.rayvector(0, 1, 0) + Raylib.rayvector(4, 5, 6),
        60,
        Raylib.CAMERA_ORTHOGRAPHIC
    ])

    Raylib.SetCameraMode(camera, Int(Raylib.CAMERA_FIRST_PERSON))
    Raylib.UpdateCamera!(camera)
    @test !test_camera(camera, [
        Raylib.rayvector(0, 10, 10) + Raylib.rayvector(1, 2, 3),
        Raylib.rayvector(0, 0, 0) + Raylib.rayvector(3, 2, 1),
        Raylib.rayvector(0, 1, 0) + Raylib.rayvector(4, 5, 6),
        60,
        Raylib.CAMERA_ORTHOGRAPHIC
    ])

    Raylib.SetCameraMode(camera, Int(Raylib.CAMERA_CUSTOM))
    Raylib.UpdateCamera!(camera)
    camera.position = Raylib.rayvector(1, 2, 3)
    camera.target = Raylib.rayvector(3, 2, 1)
    camera.up = Raylib.rayvector(4, 5, 6)
    camera.fovy = 60
    camera.projection = Raylib.CAMERA_ORTHOGRAPHIC
    @test test_camera(camera, [
        Raylib.rayvector(1, 2, 3),
        Raylib.rayvector(3, 2, 1),
        Raylib.rayvector(4, 5, 6),
        60,
        Raylib.CAMERA_ORTHOGRAPHIC
    ])

    Raylib.SetCameraMode(camera, Int(Raylib.CAMERA_FIRST_PERSON))
    camera = Raylib.UpdateCamera(camera)
    @test !test_camera(camera, [
        Raylib.rayvector(1, 2, 3),
        Raylib.rayvector(3, 2, 1),
        Raylib.rayvector(4, 5, 6),
        60,
        Raylib.CAMERA_FIRST_PERSON
    ])
end


text:
import ACE: OneParticleBasis, AbstractState
import ACE.OrthPolys: TransformedPolys

@doc raw"""
`struct Scal1pBasis <: OneParticleBasis`

One-particle basis of the form $P_n(x_i)$ for a general scalar, invariant
input `x`. This type basically just translates the `TransformedPolys` into a valid
one-particle basis.
"""
mutable struct Scal1pBasis{T, TT, TJ} <: OneParticleBasis{T}
   P::TransformedPolys{T, TT, TJ}
   varsym::Symbol
   idxsym::Symbol
end

scal1pbasis(varsym::Symbol, idxsym::Symbol, args...; kwargs...) =
            Scal1pBasis( ACE.OrthPolys.transformed_jacobi(args...; kwargs...),
                         varsym, idxsym )

# ---------------------- Implementation of Scal1pBasis

Base.length(basis::Scal1pBasis) = length(basis.P)

get_spec(basis::Scal1pBasis) =
      [  NamedTuple{(basis.idxsym,)}(n) for n = 1:length(basis) ]

==(P1::Scal1pBasis, P2::Scal1pBasis) =  ACE._allfieldsequal(P1, P2)

write_dict(basis::Scal1pBasis{T}) where {T} = Dict(
      "__id__" => "ACE_Scal1pBasis",
          "P" => write_dict(basis.P) ,
          "varsym" => string(varsym),
          "idxsym" => string(idxsym) )

read_dict(::Val{:ACE_Scal1pBasis}, D::Dict) =
      Scal1pBasis(read_dict(D["P"]), Symbol(D["varsym"]), Symbol(D["idxsym"]))

fltype(basis::Scal1pBasis{T}) where T = T

gradtype(::Scal1pBasis{T}) where {T} = SVector{3, T}

symbols(basis::Scal1pBasis) = [basis.idxsym]

indexrange(basis::Scal1pBasis) = NamedTuple{(basis.idxsym,)}((1:length(basis),))

_getidx(b, basis::Scal1pBasis) = b[basis.idxsym]

isadmissible(b, basis::Scal1pBasis) = (1 <= _getidx(b, basis) <= length(basis))

degree(b, basis::Scal1pBasis) = _getidx(b, basis) - 1

get_index(basis::Scal1pBasis, b) = _getidx(b, basis)

rand_radial(basis::Scal1pBasis) = rand_radial(basis.P)

# ---------------------------  Evaluation code
#

alloc_B(basis::Scal1pBasis, args...) = alloc_B(basis.P)

alloc_dB(basis::Scal1pBasis, args...) = zeros(fltype(basis.P), length(basis))

alloc_temp(basis::Scal1pBasis, args...) = alloc_temp(basis.P)

alloc_temp_d(basis::Scal1pBasis, args...) = alloc_temp_d(basis.P)

_getval(X::AbstractState, basis::Scal1pBasis) = getproperty(X, basis.valsym)
_getval(x::Number, basis::Scal1pBasis) = x

evaluate!(B, tmp, basis::Scal1pBasis, X::Union{AbstractState, Number}) =
      evaluate!(B, tmp, basis.P, _getval(X, basis))

function evaluate_d!(dB, tmpd, basis::Scal1pBasis, X::Union{AbstractState, Number})
   evaluate_d!(dB, tmpd, basis.P, _getval(X, basis))
   return dB
end

function evaluate_ed!(B, dB, tmpd, basis::Scal1pBasis, X::Union{AbstractState, Number})
   x = _getval(X, basis)
   evaluate!(B, tmpd, basis.P, x)
   evaluate_d!(dB, tmpd, basis.P, x)
   return nothing
end


text:
using NeXLUncertainties
using NeXLCore
using Colors
using ImageCore: colorview
using Statistics

"""
The members in common between KRatio and KRatios
"""
abstract type KRatioBase
    # element::Element
    # xrays::Vector{CharXRay} # Which CharXRays were measured?
    # unkProps::Dict{Symbol,Any} # Beam energy, take-off angle, coating, ???
    # stdProps::Dict{Symbol,Any} # Beam energy, take-off angle, coating, ???
    # standard::Material
end

"""
    Base.findfirst(krs::AbstractVector{<:KRatioBase}, cxr::CharXRay)

Find the first KRatio or KRatios in which the .xrays field contains the cxr::CharXRay.
"""
Base.findfirst(krs::AbstractVector{<:KRatioBase}, cxr::CharXRay) =
    krs[findfirst(kr -> cxr in kr.xrays, krs)]
element(kr::KRatioBase) = kr.element
xrays(krs::KRatioBase) = krs.xrays
standard(krs::KRatioBase) = krs.standard
"""
    elms(krs::Vector{KRatio})::Set{Element}

Return a set containing the elements present in krs.
"""
function elms(krs::Vector{<:KRatioBase})::Set{Element}
    res = Set{Element}()
    for kr in krs
        push!(res, kr.element)
    end
    return res
end

"""
    KRatio

The k-ratio is the ratio of two similar intensity measurements - one on
a material of unknown composition and one on a standard with known
composition. Each measurement has properties like :BeamEnergy (req),
:TakeOffAngle (req), :Coating (opt) that characterize the measurement.
A minimal set of properties includes:

Properties: (These `Symbol`s are intentionally the same used in NeXLSpectrum)

    :BeamEnergy incident beam energy in eV
    :TakeOffAngle in radians
    :Coating A NeXLCore.Film object or Film[] detailing a conductive coating

Some algorithms may require additional properties.

k-ratios are associated with characteristic X-rays (`CharXRay`) from a single element.
WDS k-ratios are typically associated with a single `CharXRay` while EDS
measurements may be associated with many `CharXRay` that are similar in energy.
k-ratios are always relative to another material.  Usually the composition of this
`Material` is well-known.  However, when a k-ratio is restandardized, it is possible
for the intermediate material to be less well-known.
"""
struct KRatio <: KRatioBase
    element::Element
    xrays::Vector{CharXRay} # Which CharXRays were measured?
    unkProps::Dict{Symbol,Any} # Beam energy, take-off angle, coating, ???
    stdProps::Dict{Symbol,Any} # Beam energy, take-off angle, coating, ???
    standard::Material
    kratio::UncertainValue

    function KRatio(
        xrays::AbstractVector{CharXRay},
        unkProps::Dict{Symbol,<:Any},
        stdProps::Dict{Symbol,<:Any},
        standard::Material,
        kratio::AbstractFloat,
    )
        if length(xrays) < 1
            error("A k-ratio must specify at least one characteristic X-ray.")
        end
        elm = element(xrays[1])
        if !all(element(l) == elm for l in xrays)
            error(
                "The characteristic X-rays in a k-ratio must all be from the same element.",
            )
        end
        if value(standard[elm]) <= 1.0e-6
            error("The standard $standard does not contain the element $(elm).")
        end
        if haskey(unkProps, :TakeOffAngle) &&
           haskey(stdProps, :TakeOffAngle) &&
           (
               !isapprox(
                   unkProps[:TakeOffAngle],
                   stdProps[:TakeOffAngle],
                   atol = deg2rad(0.1),
               )
           )
            @warn "The unknown and standard take-off angles do not match for $elm in $standard and $xrays."
        end
        return new(
            elm,
            sort(xrays, rev = true),
            copy(unkProps),
            copy(stdProps),
            standard,
            convert(UncertainValue, kratio),
        )
    end
end

NeXLUncertainties.value(kr::KRatio) = value(kr.kratio)
NeXLUncertainties.σ(kr::KRatio) = σ(kr.kratio)
nonnegk(kr::KRatio) = value(kr.kratio) < 0.0 ? uv(0.0, σ(kr.kratio)) : kr.kratio
Statistics.mean(krs::AbstractVector{KRatio})::UncertainValue =
    mean((kr.kratio for kr in krs)...)

function Base.getindex(krs::AbstractVector{KRatio}, cxr::CharXRay)
    i = findfirst(kr->cxr in kr.xrays, krs)
    return isnothing(i) ? zero(UncertainValue) : krs[i].kratio
end

"""
    strip(krs::AbstractVector{KRatio}, els::Element...)::Vector{KRatio}

Creates a new Vector{KRatio} containing all the KRatio objects in `krs` except those associated with the specified elements.
"""
Base.strip(krs::AbstractVector{KRatio}, els::Element...) =
    collect(filter(k -> !(element(k) in els), krs))

Base.show(io::IO, kr::KRatio) =
    print(io, "k[$(name(kr.xrays)), $(name(kr.standard))] = $(round(kr.kratio))")

function NeXLUncertainties.asa(::Type{DataFrame}, krs::AbstractVector{KRatio})::DataFrame
    return DataFrame(
        Symbol("X-rays") => [ repr(kr.xrays) for kr in krs ],
        Symbol("Standard") => [ name(kr.standard) for kr in krs ],
        Symbol("C[std]") => [ value(kr.standard[kr.element]) for kr in krs ],
        Symbol("ΔC[std]") => [ σ(kr.standard[kr.element]) for kr in krs ],
        Symbol("E₀[unk]") => [ get(kr.unkProps, :BeamEnergy, missing) for kr in krs ],
        Symbol("θ[unk]") => [ get(kr.unkProps, :TakeOffAngle, missing) for kr in krs ],
        Symbol("E₀[std]") => [ get(kr.stdProps, :BeamEnergy, missing) for kr in krs ],
        Symbol("θ[std]") => [ get(kr.stdProps, :TakeOffAngle, missing) for kr in krs ],
        Symbol("k") => [ value(kr.kratio) for kr in krs ],
        Symbol("Δk") => [ σ(kr.kratio) for kr in krs ],
    )
end

"""
`KRatios` represents the hyper-spectral equivalent of the KRatio type.  Each pixel in the `KRatios` object
must be characterized by the same unknown and standard properties, the same X-ray lines and the other
properties.
"""
struct KRatios{T} <: KRatioBase where {T <: AbstractFloat}
    element::Element
    xrays::Vector{CharXRay} # Which CharXRays were measured?
    unkProps::Dict{Symbol,Any} # Beam energy, take-off angle, coating, ???
    stdProps::Dict{Symbol,Any} # Beam energy, take-off angle, coating, ???
    standard::Material
    kratios::Array{T}

    function KRatios(
        xrays::Vector{CharXRay},
        unkProps::Dict{Symbol,<:Any},
        stdProps::Dict{Symbol,<:Any},
        standard::Material,
        kratios::Array{T},
    ) where { T <: AbstractFloat }
        if length(xrays) < 1
            error("A k-ratio must specify at least one characteristic X-ray.")
        end
        elm = element(xrays[1])
        if !all(element(l) == elm for l in xrays)
            error(
                "The characteristic X-rays in a k-ratio must all be from the same element.",
            )
        end
        if value(standard[elm]) <= 1.0e-6
            error("The standard $standard does not contain the element $(elm).")
        end
        if haskey(unkProps, :TakeOffAngle) &&
           haskey(stdProps, :TakeOffAngle) &&
           (
               !isapprox(
                   unkProps[:TakeOffAngle],
                   stdProps[:TakeOffAngle],
                   atol = deg2rad(0.1),
               )
           )
            @warn "The unknown and standard take-off angles do not match for $elm in $standard and $xrays."
        end
        return new{T}(elm, xrays, copy(unkProps), copy(stdProps), standard, kratios)
    end
end

Base.show(io::IO, kr::KRatios) = print(
    io,
    "k[$(name(kr.standard)), $(name(kr.xrays))] = $(eltype(kr.kratios))[ $(size(kr.kratios)) ]",
)

Base.getindex(krs::KRatios, idx::Int...) = KRatio(
    krs.xrays,
    krs.unkProps,
    krs.stdProps,
    krs.standard,
    getindex(krs.kratios, idx...),
)

Base.getindex(krs::KRatios, ci::CartesianIndex) = KRatio(
    krs.xrays,
    krs.unkProps,
    krs.stdProps,
    krs.standard,
    getindex(krs.kratios, ci),
)

Base.size(krs::KRatios) = size(krs.kratios)
Base.size(krs::KRatios, idx::Int) = size(krs.kratios, idx)
Base.CartesianIndices(krs::KRatios) = CartesianIndices(krs.kratios)

"""
    LinearAlgebra.normalize(krs::AbstractVector{<:KRatios}; norm::Float32=1.0f)::Vector{KRatios}
    LinearAlgebra.normalize(krs::AbstractVector{KRatio}; norm::Float32=1.0f)::Vector{KRatio}

Computes the pixel-by-pixel normalized k-ratio for each point in the KRatios data array. `norm` specifies normalization
constants other than 1.0.
"""
function LinearAlgebra.normalize(
    krs::AbstractVector{<:KRatios};
    norm::Float64 = 1.0,
    minsum::Float64 = 0.0
)
    sz = size(krs[1].kratios)
    @assert all((sz == size(kr.kratios) for kr in krs[2:end]))
    s = map(CartesianIndices(krs[1].kratios)) do ci
        ss = sum( max(0.0, kr.kratios[ci]) for kr in krs)
        ss<=minsum ? 1.0e100 : ss
    end
    return map(krs) do kr
        KRatios(kr.xrays, copy(kr.unkProps), copy(kr.stdProps), kr.standard, norm .* (kr.kratios ./ s))
    end
end
function LinearAlgebra.normalize(
    krs::AbstractVector{KRatio};
    norm::Float64 = 1.0,
    minsum::Float64 = 0.0
)
    s = sum(kr->max(0.0, kr.kratio), krs)
    s = s<=minsum ? 1.0e100 : s
    map(krs) do kr
        KRatio(kr.xrays, copy(kr.unkProps), copy(kr.stdProps), kr.standard, (norm/s) * kr.kratio)
    end
end

"""
    brightest(krs::Union{KRatios, KRatio})

Returns a new KRatios (referencing same basic data as krs) but with a single CharXRay in the `lines` field.
"""
brightest(krs::KRatios) = KRatios([brightest(krs.xrays)], krs.unkProps, krs.stdProps, krs.standard, krs.kratios )
brightest(krs::KRatio) = KRatio([brightest(krs.xrays)], krs.unkProps, krs.stdProps, krs.standard, krs.kratio )

"""
    colorize(krs::AbstractVector{<:KRatios}, red::Element, green::Element, blue::Element, normalize=:All[|:Each])
    colorize(krs::AbstractVector{<:KRatios}, elms::AbstractVector{Element}, normalize=:All)

Create RGB colorized images from up to three `Element`s.  The elements
are normalized relative to all `KRatios` in `krs`. The resulting images are scaled by the factor
`scale` to allow visualization of trace elements.
"""
function colorize(krs::AbstractVector{<:KRatios}, red::Element, green::Element, blue::Element, normalize=:All)
    colorize(krs, [red, green, blue], normalize)
end
function colorize(krs::AbstractVector{<:KRatios}, elms::AbstractVector{Element}, normalize=:All)
    idx  = collect( findfirst(kr->isequal(kr.element, elm), krs) for elm in elms )
    if normalize==:All
        # Normalize relative to sum of KRatios at each point
        s = map(CartesianIndices(krs[1])) do ci
            ss = sum( max(0.0, kr.kratios[ci]) for kr in krs)
            ss<=0.0 ? 1.0 : ss
        end
        clip(skr) = min(skr, 1.0)
        norm(kr) = clip.(kr.kratios ./ s)
        colorview(RGB,
            norm(krs[idx[1]]),
            length(idx)>1 ? norm(krs[idx[2]]) : zeroarray,
            length(idx)>2 ? norm(krs[idx[3]]) : zeroarray)
    else
        # Normalize relative to max of each KRatios independently
        each(kr) = kr.kratios / maximum(kr.kratios)
        colorview(RGB,
            each(krs[idx[1]]),
            length(idx)>1 ? each(krs[idx[2]]) : zeroarray,
            length(idx)>2 ? each(krs[idx[3]]) : zeroarray)
    end
end

function Base.getindex(krs::AbstractVector{<:KRatios}, elm::Element)
    colorize(krs, [ elm, elm, elm ])
end
function Base.getindex(krs::AbstractVector{<:KRatios}, red::Element, green::Element)
    colorize(krs, [ red, green ])
end
function Base.getindex(krs::AbstractVector{<:KRatios}, red::Element, green::Element, blue::Element)
    colorize(krs, [ red, green, blue ])
end


text:
using Test
using Sabermetrics

# --- Mike Trout 2019 ---
const PA = 600
const AB = 470
const HITS = 137
const DOUBLE = 27
const TRIPLE = 2
const HR = 45
const SB = 11
const CS = 2
const BB = 110
const SO = 120
const HBP = 16
const DP = 5
const SH = 0
const SF = 4
const IBB = 14

@test basehit(HITS, DOUBLE, TRIPLE, HR) == 63
@test avg(HITS, AB) == 0.291
@test obp(HITS, BB, HBP, AB, SF) == 0.438
@test slg(HITS, DOUBLE, TRIPLE, HR, AB) == 0.645
@test ops(HITS, DOUBLE, TRIPLE, HR, BB, HBP, AB, SF) == 1.083
@test iso(HITS, DOUBLE, TRIPLE, HR, AB) == 0.354
@test rc(HITS, DOUBLE, TRIPLE, HR, BB, AB) == 129.036
@test lazyrc(HITS, DOUBLE, TRIPLE, HR, BB, HBP, AB, SF) == 132.714
@test rc27(HITS, DOUBLE, TRIPLE, HR, BB, AB, CS, SH, SF, DP) == 10.128
@test babip(HITS, HR, AB, SO, SF) == 0.298
@test woba(HITS, DOUBLE, TRIPLE, HR, AB, BB, HBP, SF) == 0.455
@test wraa(HITS, DOUBLE, TRIPLE, HR, AB, PA, BB, HBP, SF, 0.33) == 62.5

# --- Gerritt Cole 2019 ---
const ER = 59
const IP = 212.1
const K = 326
const PBB = 48
const PIBB = 0
const PHBP = 3
const PHR = 29
const PH = 142

@test era(ER, IP) == 2.504
@test whip(PH, PBB, IP) == 0.896
@test fip(PHR, PBB, PIBB, PHBP, K, IP) == 2.425
@test k9(K, IP) == 13.833
@test bb9(PBB, IP) == 2.037
@test h9(PH, IP) == 6.025
@test hr9(PHR, IP) == 1.231


text:
include("common.jl")
@testset "Backward integration" begin
    h(p, t) = 1.0
    u₀ = 1.0
    tspan = (2.0, 0.0)

    f(u, h, p,t) = h(p, t + 1)
    # solution on [0,2]
    function f_analytic(u₀, p, t)
        if t < 1
            return (t^2 - 1)/2
        else
            return t - 1
        end
    end
    dde_f = DDEFunction(f, analytic = f_analytic)

    @testset "Without lags" begin
        sol = solve(DDEProblem(dde_f, u₀, h, tspan), MethodOfSteps(RK4()))
        @test sol.errors[:l∞] < 1e-5
    end

    @testset "Constant lags" begin
        # incorrect lags
        prob1 = DDEProblem(dde_f, u₀, h, tspan; constant_lags = [1.0])
        @test_throws ErrorException solve(prob1, MethodOfSteps(Tsit5()))

        prob2 = DDEProblem(dde_f, u₀, h, tspan; constant_lags = [-1.0])
        dde_int = init(prob2, MethodOfSteps(Tsit5()))
        @test dde_int.opts.d_discontinuities.valtree == [Discontinuity(1.0, 2)]
        sol = solve!(dde_int)
        @test sol.errors[:l∞] < 1e-12
        @test dde_int.tracked_discontinuities == [Discontinuity(2.0, 1),
                                                  Discontinuity(1.0, 2)]
    end

    @testset "Dependent lags" begin
        prob = DDEProblem(dde_f, u₀, h, tspan; dependent_lags = ((u, p, t) -> -1.0,))
        dde_int = init(prob, MethodOfSteps(Tsit5()))
        @test isempty(dde_int.opts.d_discontinuities)
        sol = solve!(dde_int)
        @test sol.errors[:l∞] < 1e-12
        @test dde_int.tracked_discontinuities == [Discontinuity(2.0, 1),
                                                  Discontinuity(1.0, 2)]
    end

    @testset "dt and dtmax" begin
        prob = DDEProblem(dde_f, u₀, h, tspan)

        dde_int = init(prob, MethodOfSteps(RK4()); dt = 0.1, dtmax = 0.5)
        @test dde_int.dt == -0.1 && dde_int.opts.dtmax == -0.5
        sol1 = solve!(dde_int)

        sol2 = solve(prob, MethodOfSteps(RK4()); dt = -0.1, dtmax = -0.5)
        @test sol1.t == sol2.t && sol1.u == sol2.u
    end
end


text:
using libm

N = 1_000_000
srand(1)
X = exp(randn(N))
function maxerr{T}(X::AbstractArray{T})
    e = 0.0
    for x in X
        a = libm.log1p_tang(x)
        b = log1p(big(x))
        d = abs(T(b-a))
        r = d / eps(T(b))
        e = max(e,r)
    end
    e
end

@printf "Float64 max rel. error %f ulps\n" maxerr(X)
@printf "Float32 max rel. error %f ulps\n" maxerr(float32(X))


text:
#using DifferentialEquations
make_plots = false
if make_plots
using Plots
end
using Statistics
using FFTW
FFTW.set_provider!("mkl")
#import Pkg; Pkg.add("FileIO")
using FileIO
using Base.Threads
using Distributions
#Pkg; Pkg.add("DistributedArrays")
println("nthreads = $(nthreads())")

include("aux_functions.jl")

undamped = false
damped = false
#undamped = true
damped = true

if undamped
    run_name = "landau_undamped_"
    const L = 39.738
    const J = 3522
elseif damped
    run_name = "landau_damped_"
    const L = 7.455
    const J = 930
end
order = 5
#const N = 80_000
#const N = 800_000
const N = 8_000_000
exp_Theta = 3
exp_t = 1
θ = 10.0^(-exp_Theta)
t = 0.0
t_f = 1.5*10.0^(exp_t)
M = 1501
M_g = 1500 + 1 #number of outputs, starting from the initial data
dt = t_f / (M-1)
t_i = 0.0
#M = convert(Int64,t_f/dt)
#M=1
const κ = 2π/L # for Fourier Transform
dx = L/J
x = [dx*(i-1) for i in 1:J] ;
p = (L, N, J, κ, dx, order)

println("t_f = $(t_f), M = $M, dt = $(dt), exp_Theta = $(exp_Theta)")

animation = false
phase_space_show = false

run_name = run_name * "t$(convert(Int,t_f))_L$(L)_N8_6_J$(J)_M$(M)_o$(order)_T$(exp_Theta)"
println(run_name)

#test_parameters(M, M_g, dt, 0.0, t_f)

E = zeros(J)
ϕ = zeros(J)
n = zeros(J) #charge density
S = zeros(J) #carge current
du = zeros(2*N+J); # contains r, v and E

if undamped
    #par_dis = load("Initial_Distributions/par_dis_8_5_abig_undamping.jld2", "par_dis");
    #par_dis = load("Initial_Distributions/par_dis_8_5_undamping.jld2", "par_dis");
elseif damped
    #par_dis = load("Initial_Distributions/par_dis_8_5_abig_damping.jld2", "par_dis");
    #par_dis = load("Initial_Distributions/par_dis_8_5_damping.jld2", "par_dis");
    #par_dis = load("Initial_Distributions/par_dis_8_4_damping.jld2", "par_dis");
    par_dis = load("Initial_Distributions/par_dis_8_6_damping.jld2", "par_dis");
end

#length(par_dis)
@assert length(par_dis) ÷ 2 == N

@time get_density!(par_dis, n, p)
n0 = N/L
#get_ϕ!(ϕ, -n/n0 .+ 1., κ) # chenge the sign here to make it consistent with charge conservation and the time derivative of E
@time get_ϕ!(ϕ, n/n0 .- 1., κ)
@time get_E_from_ϕ!(ϕ,E,dx)
u = [par_dis;E];
length(u)

@time get_current!(u, S, p)

println(maximum(u[1:N]))
println(minimum(u[1:N]))

@time Coordinate_test(u[1:N],L)

println("n_total = $(sum(n .- n0))")
println("v_total = $(sum(u[N+1:2N]))")
println("E_total = $(sum(E))")

println("S_total = $(sum(S))")

if make_plots
    plot(layout=(2,2))
    plot!(subplot=1,x,n/n0, title = "density", legend = :false)
    plot!(subplot=2,x,ϕ, title = "potential", legend = :false)
    plot!(subplot=3,x,E, title = "Electric Field", legend = :false)
    plot!(subplot=4,x,S, title = "Current", legend = :false)

    plot(layout=(2,2))
    histogram!(subplot=1,u[1:N], title = "density", legend = :false, bins = 300)
    histogram!(subplot=2,u[N+1:2*N], title = "velocity", legend = :false)

    histogram!(subplot=3,S, title = "current", legend = :false)

    if phase_space_show
        Plots.scatter(u[1:N],u[N+1:2N]
        , thickness_scaling = 0.3
        , markersize = 0.3
        , title = "phase-space", legend =:false)
    end
end

if nthreads() > 1
TS = zeros(J, nthreads())
p_RHSC = (N, J, L, dx, order, n, S, du, get_density!, get_current_threads!, Interpolate_2, TS) ;
else
p_RHSC = (N, J, L, dx, order, n, S, du, get_density!, get_current!, Interpolate_2) ;
end

Energy_K = zeros(M_g)
Energy_E = zeros(M_g)
E_T = zeros(M_g)
v_T = zeros(M_g)
D_T = zeros(M_g)
S_T = zeros(M_g)
T = zeros(M_g)
if animation
par = zeros(M_g,2N)
end
Energy_K[1], Energy_E[1]  = get_energy(u,(L, N, J))
E_T[1] = sum(u[2N+1:end])
v_T[1] = sum(u[N+1:2N])
get_density!(u, n, p)
get_current!(u, S, p)
D_T[1] = sum(n)/n0/J - 1
S_T[1] = sum(S)/n0/J
T[1] = var(u[N+1:2N])
if animation
par[1,:] = u[1:2N]
end

t = 0.0
j = 1

for k in 2:M
  RK4_Step!(RHSC,u,t,dt,p_RHSC)
  global u = [make_periodic!(u[1:N],L); u[N+1:end]]
  #filter_constant!(u[2N+1:end])
  global t = t + dt
  if (k-1) % (M÷(M_g-1)) == 0
    j = (k-1)÷(M÷(M_g-1))+1
    #scatter(plt, u[1:N], u[N+1:2*N])
    Energy_K[j], Energy_E[j] = get_energy(u,(L, N, J))
    E_T[j] = sum(u[2N+1:end])
    v_T[j] = sum(u[N+1:2N])
    get_density!(u, n, p)
    get_current!(u, S, p)
    D_T[j] = sum(n)/n0/J - 1
    S_T[j] = sum(S)/n0/J
    T[j] = var(u[N+1:2N])
    println("j = $j , t = $t, k = $k, nthreads = $(nthreads())")
    if animation
    par[j,:] = u[1:2N]
    end
  end
end

n_F = zeros(J)
S_F = zeros(J)
get_density!(u, n_F, p)
get_current!(u, S_F, p)
run = Dict("run_name" => run_name, "p" => p, "Energy_E" => Energy_E, "Energy_K" => Energy_K, "E_f" => u[2N+1:end], "n_F" => n_F, "S_F" => S_F, "E_T"=> E_T, "v_T" => v_T, "S_T" => S_T, "D_T" => D_T, "T" => T)
save(run_name * "th$(nthreads())_results.jld2", run)

if make_plots
    plot(abs.(Energy_K[2:end] .- Energy_K[1]), title = "Energy conservation", label = "Kinetic Energy"
    #, legend = :outertopright
    , legend = :bottomright)
    plot!(abs.(Energy_E[2:end] .- Energy_E[1]), label = "|Electric Energy|")
    plot!(abs.(Energy_K[2:end]  + Energy_E[2:end] .- (Energy_K[1]+Energy_E[1])) ./ (Energy_K[1]+Energy_E[1]) , yscale=:log10
    #, xscale=:log10
    , label = "Total Energy / Initial Energy -1 ")
    #png(run_name * "energy_conservation")

    abs.(Energy_K[end]  + Energy_E[end] .- (Energy_K[1]+Energy_E[1])) ./ (Energy_K[1]+Energy_E[1])

    #ω = 1
    #tv = [dt*(i-1) for i in 1:M]
    plot(T, label= "θ", title = "Temperature = var", legend = false)
    #plot!(sin.(ω*tv))
    #png(run_name * "temperature")
end

Th = rfft(T)
plot(real.(Th))
argmin(real.(Th))*t_f/length(T) * sqrt(θ)
#length(T)

T[end]

if make_plots
    plot(layout=(2,2))
    plot!(subplot=1,E_T,title="Total Electric Field")
    plot!(subplot=2,v_T./N, title="Total velocity")
    plot!(subplot=3,D_T,title = "Total density")
    plot!(subplot=4,S_T,title = "Total Current")
    #png(run_name * "totals")
end

E_F = zeros(J)
ϕ_F = zeros(J)
n_F = zeros(J)
get_density!(u, n_F, p)
n0 = N/L
get_ϕ!(ϕ_F, n_F/n0 .+ 1, κ)
#get_ϕ!(ϕ_F, n_F .- n0, κ)
get_E_from_ϕ!(ϕ_F,E_F,dx)

if make_plots
    plot(x,u[2N+1:end], label = "dynamical", title="Electric Field with order = $(order)")
    plot!(x,E_F,label="from constraint", ls=:dash)
    #t_f = 40
    #png(run_name * "E")
end

println("averaged total E field = $(sum(E_F))")
println("E_total = $(sum(u[2N+1:end])/J)")
println("Total velocity = $(sum(u[N+1:2N])/N)")
println("Total Charge = $(sum(n_F .- n0))")
println("Final Energy = $(get_energy(u,(L, N, J)))")
get_current!(u, S, p)
println("Total_current = $(sum(S)/J)")

if make_plots
plot(x,u[2N+1:end]-E_F, label = "difference", title="Electric Field")
#png(run_name * "electric_diff")
if make_plots

if phase_space_show
Plots.scatter(u[1:N],u[N+1:2N]
#, thickness_scaling = 0.3
, markersize = 0.3
, title = "phase-space", legend =:false)
#png(run_name * "ps")
end

if animation
    anim = @animate for i = 1:M_g
        Plots.scatter(par[i,1:N], par[i,N+1:2N]
        , markersize = 0.3
        , title = "phase-space"
        , legend=false
        , ylim = (-0.3,0.3)
        , xlim = (0,5)
        )
    end

    gif(anim, run_name * "ps_fps5.gif", fps = 5)
end

#load("Results/" * run_name * "results.jld2", "S_F")

#(p, Energy_E, Energy_K, E_f, n_F, S_F, E_T, v_T, S_T, D_T, T) = load("Results/" * run_name * "results.jld2", "p", "Energy_E", "Energy_K", "E_f", "n_F",  "S_F", "E_T", "v_T", "S_T", "D_T", "T");

#using LsqFit


text:
incircle(x) = transpose(x) * x < 1

computeπ(tol) = foldl(right,
       Map(incircle)
    |> Mean()
    |> Drop(10)
    |> Map(x -> x*4)
    |> TakeWhile(x -> x.err > tol)
    , Quasi(2)
    )

computeπ(0.001)


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE

using Test, Tensors

let E = 200.0e3,
    nu = 0.3,
    yield_strength = 100.0,
    parameters = PerfectPlasticParameterState(youngs_modulus=E,
                                              poissons_ratio=nu,
                                              yield_stress=yield_strength),  # yield in shear = R0 / sqrt(3)
    mat = PerfectPlastic(parameters=parameters),
    times = [0.0],
    loads = [0.0],
    dt = 1.0,
    G = 0.5*E/(1+nu),
    # vonMises = sqrt(3 J_2) = sqrt(3/2 tr(s^2)) = sqrt(3) |tau| = sqrt(3)*G*|gamma|
    # gamma = 2 e12
    # set  vonMises = Y
    gamma_yield = yield_strength/(sqrt(3)*G)

    # Go to elastic border
    push!(times, times[end] + dt)
    push!(loads, loads[end] + gamma_yield*dt)
    # Proceed to plastic flow
    push!(times, times[end] + dt)
    push!(loads, loads[end] + gamma_yield*dt)
    # Reverse direction
    push!(times, times[end] + dt)
    push!(loads, loads[end] - gamma_yield*dt)
    # Continue and pass yield criterion
    push!(times, times[end] + dt)
    push!(loads, loads[end] - 2*gamma_yield*dt)

    stresses = [copy(tovoigt(mat.variables.stress))]
    for i=2:length(times)
        dtime = times[i] - times[i-1]
        dstrain31 = loads[i] - loads[i-1]
        dstrain_voigt = [0.0, 0.0, 0.0, 0.0, 0.0, dstrain31]
        dstrain_tensor = fromvoigt(Symm2{Float64}, dstrain_voigt; offdiagscale=2.0)
        mat.ddrivers = PerfectPlasticDriverState(time=dtime, strain=dstrain_tensor)
        integrate_material!(mat)
        update_material!(mat)
        push!(stresses, copy(tovoigt(mat.variables.stress)))
    end

    for i in 1:length(times)
        @test isapprox(stresses[i][1:5], zeros(5); atol=1e-6)
    end

    let y = yield_strength/sqrt(3.0)
        s31 = [s[6] for s in stresses]
        s31_expected = [0.0, y, y, 0.0, -y]
        @test isapprox(s31, s31_expected; rtol=1.0e-2)
    end
end


text:
#####
##### Codegen
#####

struct BrutusCompilerTarget <: AbstractCompilerTarget end
GPUCompiler.llvm_triple(::BrutusCompilerTarget) = Sys.MACHINE
GPUCompiler.llvm_machine(::BrutusCompilerTarget) = tm[]

module Runtime
    # the runtime library
    signal_exception() = return
    malloc(sz) = Base.Libc.malloc(sz)
    report_oom(sz) = return
    report_exception(ex) = return
    report_exception_name(ex) = return
    report_exception_frame(idx, func, file, line) = return
end

@enum DumpOption::UInt8 begin
    DumpIRCode        = 0
    DumpTranslated    = 1
    DumpCanonicalized = 2
    DumpLoweredToStd  = 4
    DumpLoweredToLLVM = 8
    DumpTranslateToLLVM = 16
end

struct BrutusCompilerParams <: AbstractCompilerParams
    emit_fptr::Bool
    dump_options::Vector{DumpOption}
end

GPUCompiler.ci_cache(job::CompilerJob{BrutusCompilerTarget}) = GLOBAL_CI_CACHE
GPUCompiler.runtime_module(job::CompilerJob{BrutusCompilerTarget}) = Runtime
GPUCompiler.isintrinsic(::CompilerJob{BrutusCompilerTarget}, fn::String) = true
GPUCompiler.can_throw(::CompilerJob{BrutusCompilerTarget}) = true
GPUCompiler.runtime_slug(job::CompilerJob{BrutusCompilerTarget}) = "brutus"

function find_invokes(IR)
    callees = Core.MethodInstance[]
    for stmt in IR.stmts
        if stmt isa Expr
            if stmt.head == :invoke
                mi = stmt.args[1]
                push!(callees, mi)
            end
        end
    end
    return callees
end

# Emit MLIR IR to stdout
function emit(job::CompilerJob)
    ft = job.source.f
    tt = job.source.tt
    emit_fptr = job.params.emit_fptr
    dump_options = job.params.dump_options
    name = (ft <: Function) ? nameof(ft.instance) : nameof(ft)

    # get first method instance matching signature
    entry_mi = get_methodinstance(Tuple{ft, tt.parameters...})
    IR, rt = code_ircode(entry_mi)

    if DumpIRCode in dump_options
        println("return type: ", rt)
        println("IRCode:\n")
        println(IR)
    end

    worklist = [IR]
    methods = Dict{Core.MethodInstance, Tuple{Core.Compiler.IRCode, Any}}(
        entry_mi => (IR, rt)
    )

    while !isempty(worklist)
        code = pop!(worklist)
        callees = find_invokes(code)
        for callee in callees
            if !haskey(methods, callee)
                _code, _rt = code_ircode(callee)

                methods[callee] = (_code, _rt)
                push!(worklist, _code)
            end
        end
    end

    # generate LLVM bitcode and load it
    dump_flags = reduce(|, map(UInt8, dump_options), init=0)
    fptr = ccall((:brutus_codegen, "libbrutus"),
                 Ptr{Nothing},
                 (Any, Any, Cuchar, Cuchar),
                 methods, entry_mi, emit_fptr, dump_flags)
    return (fptr, rt)
end

function emit(@nospecialize(ft), @nospecialize(tt);
              emit_fptr::Bool=true,
              dump_options::Vector{DumpOption}=DumpOption[])
    fspec = GPUCompiler.FunctionSpec(ft, Tuple{tt...}, false, nothing)
    target = BrutusCompilerTarget()
    params = BrutusCompilerParams(emit_fptr, dump_options)
    job = CompilerJob(target, fspec, params)
    return emit(job)
end

max_stars_repo_name: stephenneuendorffer/brutus

max_issues_repo_name: stephenneuendorffer/brutus

max_forks_repo_name: stephenneuendorffer/brutus


text:
include("./BinaryIO.jl")

module ProgramTunnel

using Formatting
using ..BinaryIO

export recvText, sendText, recvBinary!, sendBinary!,
       mkTunnel, defaultTunnelSet, getTunnelFilename!, reverseRole!

mutable struct Tunnel
    fns      :: AbstractArray{AbstractString}
    next_ptr :: Integer

    function Tunnel(name::AbstractString)
        fns = Array{AbstractString}(undef, 2)
        for i in 1:length(fns)
            fns[i] = format("_{}_{:d}.fifo", name, i)
        end
        return new(fns, 1)
    end
end

mutable struct TunnelSet

    tnls :: Dict{Symbol, Tunnel}

    function TunnelSet(;
        recv_txt :: AbstractString,
        send_txt :: AbstractString,
        recv_bin :: AbstractString,
        send_bin :: AbstractString,
        path :: AbstractString = ".",
    )

        tnls = Dict{Symbol, Tunnel}(
            :recv_txt => Tunnel(recv_txt),
            :send_txt => Tunnel(send_txt),
            :recv_bin => Tunnel(recv_bin),
            :send_bin => Tunnel(send_bin),
        )

        for (k, tnl) in tnls
            for i in 1:length(tnl.fns)
                tnl.fns[i] = normpath(joinpath(path, tnl.fns[i]))
            end

        end

        return new(tnls)
    end

end

function getTunnelFilename!(tnl::Tunnel)
    next_ptr = tnl.next_ptr
    tnl.next_ptr = mod(next_ptr, length(tnl.fns)) + 1
    return tnl.fns[next_ptr]
end

function getTunnelFilename!(ts::TunnelSet, key::Symbol)
    return getTunnelFilename!(ts.tnls[key])
end

function defaultTunnelSet(;path::AbstractString=".")
    return TunnelSet(
        send_txt = "X2Y_txt",
        recv_txt = "Y2X_txt",
        send_bin = "X2Y_bin",
        recv_bin = "Y2X_bin",
        path = path,
    )
end

function mkTunnel(TS::TunnelSet)

    for (k, tnl) in TS.tnls
        for fn in tnl.fns
            if !isfifo(fn)
                println(fn, " is not a fifo or does not exist. Remove it and create a new one.")
                rm(fn, force=true)
                run(`mkfifo $fn`)
            end
        end
    end
end

function reverseRole!(TS::TunnelSet)
    TS.tnls[:send_txt], TS.tnls[:recv_txt] = TS.tnls[:recv_txt], TS.tnls[:send_txt]
    TS.tnls[:send_bin], TS.tnls[:recv_bin] = TS.tnls[:recv_bin], TS.tnls[:send_bin]
end

function recvText(TS::TunnelSet)
    local result

    open(getTunnelFilename!(TS, :recv_txt), "r") do io
        result = strip(read(io, String))
    end

    return result
end

function sendText(TS::TunnelSet, msg::AbstractString)

    open(getTunnelFilename!(TS, :send_txt), "w") do io
        write(io, msg)
    end
end

function recvBinary!(
    TS        :: TunnelSet,
    #msg       :: AbstractString,
    arr       :: AbstractArray{Float64},
    buffer    :: AbstractArray{UInt8};
    endianess :: Symbol=:little_endian,
)
    #recv_msg = recvText(TS)
    #println("[", recv_msg, "]")
    #if recv_msg != msg
    #    throw(ErrorException("Expect message ", msg, ", but get ", recv_msg))
    #end

    BinaryIO.readBinary!(
        getTunnelFilename!(TS, :recv_bin),
        arr,
        buffer;
        endianess=endianess
    )
end

function sendBinary!(
    TS        :: TunnelSet,
    #msg       :: AbstractString,
    arr       :: AbstractArray{Float64},
    buffer    :: AbstractArray{UInt8};
    endianess :: Symbol=:little_endian,
)

    #sendText(TS, msg)
    BinaryIO.writeBinary!(
        getTunnelFilename!(TS, :send_bin),
        arr,
        buffer;
        endianess=endianess
    )
end

function hello(TS::TunnelSet)
    sendText(TS, "<<TEST>>")
    recv_msg = recvText(TS)
    if recv_msg != "<<TEST>>"
        throw(ErrorException("Weird message: " * recv_msg))
    end
end

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
const DATASETPARAMS = [:is_sparse, :max_bin, :data_random_seed, :categorical_feature]

const BOOSTERPARAMS = [
    :application, :boosting, :learning_rate, :num_leaves, :max_depth, :tree_learner,
    :num_threads, :histogram_pool_size, :min_data_in_leaf, :min_sum_hessian_in_leaf,
    :lambda_l1, :lambda_l2, :min_gain_to_split, :feature_fraction, :feature_fraction_seed,
    :bagging_fraction, :bagging_freq, :bagging_seed, :early_stopping_round, :sigmoid,
    :is_unbalance, :drop_rate, :max_drop, :skip_drop, :xgboost_dart_mode, :uniform_drop,
    :drop_seed, :top_rate, :other_rate, :metric, :is_training_metric, :ndcg_at, :num_machines,
    :local_listen_port, :time_out, :machine_list_file, :num_class, :device_type,
]

const INDEXPARAMS = [:categorical_feature]

const MAXIMIZE_METRICS = ["auc", "ndcg"]

# LOGGING Funcs
function log_fatal(verbosity::Integer, msg...)
    warn(msg...)
end

function log_warning(verbosity::Integer, msg...)
    verbosity >= 0 && warn(msg...)
end

function log_info(verbosity::Integer, msg...)
    verbosity >= 1 && print(msg...)
end

function log_debug(verbosity::Integer, msg...)
    verbosity >= 2 && print(msg...)
end

"""
    savemodel(estimator, filename; [num_iteration = -1])

Save the fitted model in `estimator` as `filename`.

# Arguments
* `estimator::LGBMEstimator`: the estimator to use in the prediction.
* `filename::String`: the name of the file to save the model in.
* `num_iteration::Integer`: keyword argument that sets the number of iterations of the model that
    should be saved. `< 0` for all iterations.
"""
function savemodel(estimator::LGBMEstimator, filename::String; num_iteration::Integer = -1,start_iteration::Integer=0)
    @assert(estimator.booster.handle != C_NULL, "Estimator does not contain a fitted model.")
    LGBM_BoosterSaveModel(estimator.booster, start_iteration , num_iteration , filename)
    return nothing
end

"""
    loadmodel(estimator, filename)

Load the fitted model `filename` into `estimator`. Note that this only loads the fitted model—not
the parameters or data of the estimator whose model was saved as `filename`.

# Arguments
* `estimator::LGBMEstimator`: the estimator to use in the prediction.
* `filename::String`: the name of the file that contains the model.
"""
function loadmodel(estimator::LGBMEstimator, filename::String)
    estimator.booster = LGBM_BoosterCreateFromModelfile(filename)
    return nothing
end

function shrinkresults!(results, last_retained_iter::Integer)
    for test_key in keys(results)
        test = results[test_key]
        for metric_key in keys(test)
            test[metric_key] = test[metric_key][1:last_retained_iter]
        end
    end
    return nothing
end

function tryload!(estimator::LGBMEstimator)

    if estimator.booster.handle == C_NULL
        # first check for a serialised model
        if length(estimator.model) == 0
            throw(ErrorException("Estimator does not contain a fitted model."))
        end
        # load it
        estimator.booster = LGBM_BoosterLoadModelFromString(estimator.model)
    end

    return nothing
end

function get_iter_number(estimator::LGBMEstimator)

    if estimator.booster.handle == C_NULL
        # We cannot call LGBM_BoosterGetCurrentIteration without an initialised booster
        throw(ErrorException("Estimator does not contain any form of booster"))
    end

    return LGBM_BoosterGetCurrentIteration(estimator.booster)

end

function feature_importance_wrapper(estimator::LGBMEstimator, importance_type::Integer, num_iteration::Integer)
    # reason why main func isn't marked as a mutator because it isnt an "important" mutation
    tryload!(estimator)
    return LGBM_BoosterFeatureImportance(estimator.booster, num_iteration, importance_type)

end

"""
    gain_importance(estimator, num_iteration)
    gain_importance(estimator)

    Returns the importance of a fitted booster in terms of information gain across
    all boostings, or up to `num_iteration` boostings
"""
gain_importance(estimator::LGBMEstimator, num_iteration::Integer) = feature_importance_wrapper(estimator, 1, num_iteration)
gain_importance(estimator::LGBMEstimator) = feature_importance_wrapper(estimator, 1, 0)

"""
    split_importance(estimator, num_iteration)
    split_importance(estimator)

    Returns the importance of a fitted booster in terms of number of times feature was
    used in a split across all boostings, or up to `num_iteration` boostings
"""
split_importance(estimator::LGBMEstimator, num_iteration::Integer) = feature_importance_wrapper(estimator, 0, num_iteration)
split_importance(estimator::LGBMEstimator) = feature_importance_wrapper(estimator, 0, 0)


text:
@testset "event.jl" begin

    @event test_event() begin
        return true
    end

    @test test_event().name == "test_event()"
    @test test_event().proc() === true
    @test typeof(test_event()) === Event

    @testset "Event Errors" begin
        @test_throws ArgumentError @m_throw @event foo begin true end
    end

    @testset "Macro Side Effects" begin
        @event foo() 4

        @test foo().proc() === 4
    end
end


text:
using Test
using Scruff
using Scruff.Operators
using Scruff.MultiInterface
using Scruff.SFuncs
using Scruff.Utils
import Distributions

function test_support(sf::SFunc{I,O}, parranges, target, quality;
        size = 100, curr = O[]) where {I,O}
    s = support(sf, parranges, size, curr)
    @test Set(s) == Set(target)
    @test support_quality(sf, parranges) == quality
end

function test_sample(sf::SFunc{I,O}, parvals, range, probs;
        num_samples = 1000, tolerance = 0.05) where {I,O}
    d = Dict{O, Int}()
    for i in 1:num_samples
        x = sample(sf, parvals)
        d[x] = get(d,x,0) + 1
    end
    for (x,p) in zip(range, probs)
        @test isapprox(d[x] / num_samples, p; atol = tolerance)
    end
end

function test_pi(pi, range, probs)
    for (x,p) in zip(range, probs)
        @test isapprox(cpdf(pi, (), x), p)
    end
end

@testset "SFuncs" begin

    @testset "Constant" begin
        c = Constant(2)
        test_support(c, (), [2], :CompleteSupport)
        test_sample(c, (), [2], 1.0)
        @test logcpdf(c, (), 2) == 0.0
        @test logcpdf(c, (), 1) == -Inf
        (lfs, ufs) = make_factors(c, [2], (), 1, ())
        @test length(lfs) == 1
        lf = lfs[1]
        @test lf.keys == (1,)
        @test lf.dims == (1,)
        @test lf.entries == [1.0]
        @test ufs == lfs
        #=
        s1 = initial_stats(c)
        @test s1 === nothing
        s2 = expected_stats(c, [2], (), (), SoftScore([2], [1.0]))
        @test s2 === nothing
        s3 = accumulate_stats(c, s1, s2)
        @test s3 === nothing
        ps = maximize_stats(c, s3)
        @test ps === nothing
        =#
        pi = compute_pi(c, [1,2], (), ())
        test_pi(pi, [1,2], [0.0, 1.0])
    end

    @testset "Cat" begin
        c = Cat([:a,:b,:c], [0.2, 0.3, 0.5])
        c2 = Cat([:a => 0.2, :b => 0.3, :c => 0.5])
        c3 = Cat([1,1,2], [0.1, 0.3, 0.6]) # must handle duplicates in range correctly
        @test c2.range == [:a,:b,:c]
        @test c2.params == [0.2, 0.3, 0.5]
        test_support(c, (), [:a,:b,:c], :CompleteSupport)
        test_support(c3, (), [1,2], :CompleteSupport)
        test_sample(c, (), [:a,:b,:c], [0.2, 0.3, 0.5])
        @test isapprox(logcpdf(c, (), :a), log(0.2))
        @test isapprox(logcpdf(c, (), :b), log(0.3))
        @test isapprox(logcpdf(c, (), :c), log(0.5))
        @test logcpdf(c,(),:d) == -Inf
        @test isapprox(logcpdf(c3, (), 1), log(0.1 + 0.3))
        @test isapprox(logcpdf(c3, (), 2), log(0.6))
        (lfs, ufs) = make_factors(c, [:a, :b, :c], (), 1, ())
        @test length(lfs) == 1
        lf = lfs[1]
        @test lf.keys == (1,)
        @test lf.dims == (3,)
        @test lf.entries == [0.2, 0.3, 0.5]
        @test ufs == lfs
        #=
        @test initial_stats(c) == [0, 0, 0]
        chlam1 = SoftScore(c.range, [0.3, 0.1, 0.2])
        chlam2 = SoftScore(c.range, [0.2, 0.3, 0.1])
        s1 = expected_stats(c, [:a,:b,:c], (), (), chlam1)
        @test isapprox(s1, [0.06, 0.03, 0.1])
        s2 = expected_stats(c, [:a,:b,:c], (), (), chlam2)
        @test isapprox(s2, [0.04, 0.09, 0.05])
        s3 = accumulate_stats(c, s1, s2)
        @test isapprox(s3[1], 0.06 + 0.04)
        @test isapprox(s3[2], 0.03 + 0.09)
        @test isapprox(s3[3], 0.1 + 0.05)
        ps2 = maximize_stats(c, s3)
        z = sum(s3)
        @test isapprox(ps2[1], s3[1] / z)
        @test isapprox(ps2[2], s3[2] / z)
        @test isapprox(ps2[3], s3[3] / z)
        =#
        ps3 = compute_pi(c, [:a, :b, :c], (), ())
        test_pi(ps3, [:a, :b, :c], [0.2, 0.3, 0.5])

        cf = Cat(["abc", "defg"], [0.1, 0.9])
        e = 3 * 0.1 + 4 * 0.9
        @test isapprox(f_expectation(cf, (), length), e)
    end

    @testset "Cat with duplicates" begin
        c = Cat([:a,:b,:a], [0.2, 0.3, 0.5])
        test_support(c, (), [:a,:b], :CompleteSupport)
        test_sample(c, (), [:a,:b], [0.7, 0.3])
    end

    @testset "Flip" begin
        f = Flip(0.7)
        test_support(f, (), [false, true], :CompleteSupport)
        test_sample(f, (), [false, true], [0.3, 0.7])
        (lfs, ufs) = make_factors(f, [false, true], (), 7, ())
        @test length(lfs) == 1
        lf = lfs[1]
        @test lf.keys == (7,)
        @test lf.dims == (2,)
        @test isapprox(lf.entries, [0.3, 0.7])
        @test ufs == lfs
        #=
        @test initial_stats(f) == [0.0, 0.0]
        chlam1 = SoftScore(f.range, [0.2, 0.3])
        chlam2 = SoftScore(f.range, [0.4, 0.5])
        s1 = expected_stats(f, [false, true], (), (), chlam1)
        @test isapprox(s1[1], 0.06)
        @test isapprox(s1[2], 0.21)
        s2 = expected_stats(f, [false, true], (), (), chlam2)
        @test isapprox(s2[1], 0.12)
        @test isapprox(s2[2], 0.35)
        s3 = accumulate_stats(f, s1, s2)
        @test isapprox(s3[1], 0.18)
        @test isapprox(s3[2], 0.56)
        ps = maximize_stats(f, s3)
        @test isapprox(ps[2], 0.56 / 0.74)
        set_params!(f, ps)
        =#
        @test isapprox(compute_pi(f, [false, true], (), ()).params, [0.3, 0.7])
    end

    @testset "Normal" begin
        n = SFuncs.Normal(-1.0,1.0)
        dist = Distributions.Normal(-1.0, 1.0)
        empty = Vector{Float64}()

        test_support(n, (), [-1.0], :IncrementalSupport; size = 1)
        test_support(n, (), [-2.0, -1.0, 0.0], :IncrementalSupport; size = 3)
        test_support(n, (), [1.0, 2.0, 3.0, 4.0], :IncrementalSupport;
            size = 3, curr = [1.0, 2.0, 3.0, 4.0])
        test_support(n, (), [-2.0, -1.0, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0],
            :IncrementalSupport; size = 5, curr = [1.0, 2.0, 3.0, 4.0])
        range = [-2.0, -1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5,
                4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5]
        test_support(n, (), range,
            :IncrementalSupport; size = 15, curr = [1.0, 2.0, 3.0, 4.0])

        c = 0
        tot = 1000
        for i = 1:tot
            if sample(n, ()) < 0.5
                c += 1
            end
        end
        @test isapprox(c / tot, Distributions.cdf(dist, 0.5); atol = 0.05)

        @test isapprox(logcpdf(n, (), 0.5), Distributions.logpdf(dist, 0.5))

        (lfs, ufs) = make_factors(n, range, (), 2, ())
        @test length(lfs) == 1
        @test length(ufs) == 1
        lf = lfs[1]
        uf = ufs[1]
        @test lf.keys == (2,)
        @test uf.keys == lf.keys
        @test lf.dims == (length(range),)
        @test uf.dims == lf.dims
        ls = lf.entries
        us = uf.entries
        @test length(ls) == length(range)
        @test length(us) == length(ls)
        for i = 1:length(range)
            @test ls[i] >= 0
            @test us[i] <= 1
            @test ls[i] <= us[i]
            otherl = 0
            otheru = 0
            for j = 1:length(range)
                if j != i
                    otherl += ls[j]
                    otheru += ls[j]
                end
            end
            @test us[i] <= 1 - otherl
            @test otheru <= 1 - ls[i]
        end

        pi = compute_pi(n, range, (), ())
        probs = [Distributions.pdf(dist, x) for x in range]
        test_pi(pi, range, probs)
    end

    @testset "Det" begin
        f(i :: Float64, j :: Float64) = Int(floor(i+j))
        d = Det(Tuple{Float64, Float64}, Int, f)
        parranges = ([1.1, 2.2], [3.3, 4.4, 5.5])
        pis = (Cat([1.1, 2.2], [0.4, 0.6]), Cat([3.3, 4.4, 5.5], [0.2,0.3,0.5]))

        test_support(d, parranges, [4,5,6,7], :CompleteSupport)
        c1 = Constant(1)
        c2 = Constant(4)
        test_sample(d, (1.1, 4.4), [5], 1.0)
        @test logcpdf(d, (1.1, 4.4), 5) == 0.0
        @test logcpdf(d, (1.1, 4.4), 4) == -Inf
        a = [1.0,0.0,0.0,0.0, 0.0,1.0,0.0,0.0, 0.0,0.0,1.0,0.0, 0.0,1.0,0.0,0.0, 0.0,0.0,1.0,0.0, 0.0,0.0,0.0,1.0]
        (lfs, ufs) = make_factors(d, [4,5,6,7], parranges, 1, (2, 3))
        @test length(lfs) == 1
        lf = lfs[1]
        @test lf.keys == (2,3,1)
        @test lf.dims == (2,3,4)
        @test lf.entries == a
        @test ufs == lfs
        #=
        s1 = initial_stats(d)
        @test s1 === nothing
        s2 = expected_stats(d, [4,5,6,7], parranges, (), SoftScore(Vector{Int}(), Vector{Float64}()))
        @test s2 === nothing
        s3 = accumulate_stats(d, s1, s2)
        @test s3 === nothing
        ps = maximize_stats(d, s3)
        @test ps === nothing
        set_params!(d, ps)
        =#
        p4 = 0.4 * 0.2
        p5 = 0.4 * 0.3 + 0.6 * 0.2
        p6 = 0.4 * 0.5 + 0.6 * 0.3
        p7 = 0.6 * 0.5

        pi1 = compute_pi(d, [4,5,6,7], parranges, pis)
        test_pi(pi1, [4,5,6,7], [p4, p5, p6, p7])

        chlam1 = SoftScore([4,5,6,7], [0.1, 0.2, 0.3, 0.4])
        lam11 = send_lambda(d, chlam1, [4,5,6,7], parranges, pis, 1)
        lam12 = send_lambda(d, chlam1, [4,5,6,7], parranges, pis, 2)
        l11 = 0.2 * 0.1 + 0.3 * 0.2 + 0.5 * 0.3 # pi(v2) * chlam(f(1.1, v2))
        l12 = 0.2 * 0.2 + 0.3 * 0.3 + 0.5 * 0.4 # pi(v2) * chlam(f(2.2, v2))
        l23 = 0.4 * 0.1 + 0.6 * 0.2 # pi(v1) * chlam(f(v1, 3.3))
        l24 = 0.4 * 0.2 + 0.6 * 0.3 # pi(v1) * chlam(f(v1, 4.4))
        l25 = 0.4 * 0.3 + 0.6 * 0.4 # pi(v1) * chlam(f(v1, 5.5))
        @test isapprox(get_score(lam11, 1.1), l11)
        @test isapprox(get_score(lam11, 2.2), l12)
        @test isapprox(get_score(lam12, 3.3), l23)
        @test isapprox(get_score(lam12, 4.4), l24)
        @test isapprox(get_score(lam12, 5.5), l25)

        # test incremental support
        @test issubset([4], support(d, parranges, 3, [4])) == true
        @test issubset([4,5], support(d, parranges, 3, [4,5])) == true
        @test issubset([4,5,6], support(d, parranges, 3, [4,5,6])) == true
        @test issubset([4,5,6,7], support(d, parranges, 3, [4,5,6,7])) == true
        # test size in support
        @test length(support(d, parranges, 3, [4])) == 3
        @test length(support(d, parranges, 3, [4,5])) == 3
        @test length(support(d, parranges, 3, [4,5])) == 3
        @test length(support(d, parranges, 3, [4,5,5,5])) == 3
        @test length(support(d, parranges, 3, [4,5,6,7])) == 4
        @test length(support(d, parranges, 50, collect(1:60))) == 4
        @test length(support(d, parranges, 50, collect(1:30))) == 4
    end

    @testset "DiscreteCPT" begin
        d = Dict((:x,1) => [0.1,0.2,0.7], (:x,2) => [0.2,0.3,0.5], (:x,3) => [0.3,0.4,0.3],
                (:y,1) => [0.4,0.5,0.1], (:y,2) => [0.5,0.1,0.4], (:y,3) => [0.6,0.2,0.2])
        # A bug was uncaught because this range was originally in alphabetical order!
        range = ['c', 'b', 'a'] # note reverse order
        c = DiscreteCPT(range, d)
        pis = ([0.4,0.6], [0.2,0.3,0.5])
        picat1 = Cat([:x, :y], pis[1])
        picat2 = Cat([1,2,3], pis[2])
        picats = (picat1, picat2)
        # Cat range is in arbitrary order so we need to get it directly from the Cat
        parranges = (picat1.range, picat2.range)
        ks = collect(keys(d))
        ks1 = unique(first.(ks))
        ks2 = unique(last.(ks))
        test_support(c, parranges, range, :CompleteSupport)
        test_sample(c, (:x,2), ['a', 'b', 'c'], [0.5, 0.3, 0.2]) # note reverse order

        l3 = logcpdf(c, (:x,2), 'b')
        l4 = logcpdf(c, (:x,2), 'd')
        @test isapprox(l3, log(0.3))
        @test isapprox(l4, -Inf)

        k1 = nextkey()
        k2 = nextkey()
        k3 = nextkey()
        (lfs, ufs) = make_factors(c, range, parranges, k1, (k2, k3))
        @test length(ufs) == length(lfs)
        for i = 1:length(ufs)
            @test ufs[i].entries == lfs[i].entries
        end
        @test length(lfs) == 8 # six for the cases and two for the switches
        switchfact1 = lfs[7]
        switchkeys1 = switchfact1.keys
        @test length(switchkeys1) == 2
        @test switchkeys1[1] == k2
        switchkey = switchkeys1[2]
        @test switchfact1.dims == (2,6)
        @test switchfact1.entries ==
            [
                1.0, 1.0, 1.0, 0.0, 0.0, 0.0,
                0.0, 0.0, 0.0, 1.0, 1.0, 1.0,
            ]
        switchfact2 = lfs[8]
        @test switchfact2.keys == (k3, switchkey)
        @test switchfact2.dims == (3,6)
        @test switchfact2.entries ==
            [
                1.0, 0.0, 0.0, 1.0, 0.0, 0.0,
                0.0, 1.0, 0.0, 0.0, 1.0, 0.0,
                0.0, 0.0, 1.0, 0.0, 0.0, 1.0
            ]
        for i = 1:2
            for j = 1:3
                switchval = (i-1)*3 + j
                infact = lfs[switchval]
                @test infact.keys == (k1, switchkey)
                @test infact.dims == (3,6)
                es = infact.entries
                x = parranges[1][i]
                y = parranges[2][j]
                ps = d[(x,y)]
                for k = 1:3
                    for l = 1:6
                        n = (k-1)*6 + l
                        if l == switchval
                            @test es[n] == ps[k]
                        else
                            @test es[n] == 1.0
                        end
                    end
                end
            end
        end

        chlam1 = SoftScore(range, [0.9,0.8,0.7])
        chlam2 = SoftScore(range, [0.6,0.5,0.4])
        #=
        @test isempty(initial_stats(c))

        s1 = expected_stats(c, range, parranges, picats, chlam1)

        ks1 = collect(keys(s1))
        @test length(ks1) == 6
        @test (:x,1) in ks1
        @test (:x,2) in ks1
        @test (:x,3) in ks1
        @test (:y,1) in ks1
        @test (:y,2) in ks1
        @test (:y,3) in ks1
        pix2 = cpdf(picats[1], (), :x) * cpdf(picats[2], (), 2)
        nm1 = [cpdf(c, (:x,2), range[i]) * get_score(chlam1, range[i]) for i in 1:3]
        @test isapprox(s1[(:x,2)], pix2 .* nm1)

        s2 = expected_stats(c, range, parranges, picats, chlam2)
        nm2 = [cpdf(c, (:x,2), range[i]) * get_score(chlam2, range[i]) for i in 1:3]
        @test isapprox(s2[(:x,2)], pix2 .* nm2)
        s3 = accumulate_stats(c, s1, s2)
        ks3 = collect(keys(s3))
        @test length(ks1) == 6
        @test (:x,1) in ks1
        @test (:x,2) in ks1
        @test (:x,3) in ks1
        @test (:y,1) in ks1
        @test (:y,2) in ks1
        @test (:y,3) in ks1
        @test isapprox(s3[(:x,2)], s1[(:x,2)] .+ s2[(:x,2)])
        ps = maximize_stats(c, s3)
        for i in 1:2
            for j in 1:3
                x = parranges[1][i]
                y = parranges[2][j]
                m = c.inversemaps[1][x]
                n = c.inversemaps[2][y]
                @test ps[(m-1)*3+n] == normalize(s3[(x,y)])
            end
        end
=#
        p = compute_pi(c, range, parranges, picats)
        q1 = pis[1][1] .* (pis[2][1] .* d[(:x, 1)] .+ pis[2][2] .* d[(:x, 2)] .+ pis[2][3] .* d[(:x, 3)])
        q2 = pis[1][2] .* (pis[2][1] .* d[(:y, 1)] .+ pis[2][2] .* d[(:y, 2)] .+ pis[2][3] .* d[(:y, 3)])
        q = q1 .+ q2
        test_pi(p, range, q)

        # FIXME cannot test send_lambda
        l1 = send_lambda(c, chlam1, range, parranges, picats, 1)
        l2 = send_lambda(c, chlam1, range, parranges, picats, 2)
        b1x = 0.0
        b1y = 0.0
        chl1 = [get_score(chlam1, i) for i in range]
        for j = 1:3
            p = cpdf(picats[2], (), parranges[2][j])
            qx = [cpdf(c, (:x, parranges[2][j]), r) for r in range]
            qy = [cpdf(c, (:y, parranges[2][j]), r) for r in range]
            b1x += p * sum(qx .* chl1)
            b1y += p * sum(qy .* chl1)
        end
        @test isapprox(get_score(l1, :x), b1x)
        @test isapprox(get_score(l1, :y), b1y)
        b21 = 0.0
        b22 = 0.0
        b23 = 0.0
        for i = 1:2
            p = cpdf(picats[1], (), parranges[1][i])
            q1 = [cpdf(c, (parranges[1][i], 1), r) for r in range]
            q2 = [cpdf(c, (parranges[1][i], 2), r) for r in range]
            q3 = [cpdf(c, (parranges[1][i], 3), r) for r in range]
            b21 += p * sum(q1 .* chl1)
            b22 += p * sum(q2 .* chl1)
            b23 += p * sum(q3 .* chl1)
        end
        @test isapprox(get_score(l2, 1), b21)
        @test isapprox(get_score(l2, 2), b22)
        @test isapprox(get_score(l2, 3), b23)
    end

    @testset "LinearGaussian" begin
        lg = LinearGaussian((-1.0, 1.0, 2.0), 3.0, 1.0)
        pars = ([0.0, 1.0], [2.0], [3.0, 4.0, 5.0])
        v1 = support(lg, pars, 10, Vector{Float64}())
        v2 = support(lg, pars, 100, v1)
        @test support_quality(lg, pars) == :IncrementalSupport
        @test length(v1) >= 10
        @test length(v2) >= 100
        @test all(v -> v in v2, v1)
    end

    @testset "CLG" begin
        d = Dict((:x,1) => ((-1.0, 1.0, 2.0), 3.0, 1.0), (:x,2) => ((-2.0, 4.0, 2.0), 3.0, 1.0),
                (:x,3) => ((-3.0, 2.0, 2.0), 3.0, 1.0), (:y,1) => ((-4.0, 5.0, 2.0), 3.0, 1.0),
                (:y,2) => ((-5.0, 3.0, 2.0), 3.0, 1.0), (:y,3) => ((-6.0, 6.0, 2.0), 3.0, 1.0))
        clg = CLG(d)
        pars = ([:x, :y], [1, 2, 3], [0.0, 1.0], [2.0], [3.0, 4.0, 5.0])
        v1 = support(clg, pars, 10, Vector{Float64}())
        v2 = support(clg, pars, 100, v1)
        v3 = support(clg, pars, 1000, v2)
        @test support_quality(clg, pars) == :IncrementalSupport
        @test length(v1) >= 10
        @test length(v2) >= 100
        @test length(v3) >= 1000
        @test all(v -> v in v2, v1)
        @test all(v -> v in v3, v2)

        # CLG with 1 discrete and 0 continuos parents
        d2 = Dict((:x,) => ((), 0.0, 0.1), (:y,) => ((), 0.5, 0.2), (:z,) => ((), 1.5, 0.2))
        clg2 = CLG(d2)
        pars = ([:x, :y, :z],)
        v2 = support(clg2, pars, 100, Float64[])
    end

    @testset "Mixture" begin
        s1 = Flip(0.9)
        s2 = Cat([true, false], [0.2, 0.8]) # order of values reversed
        mx1 = Mixture([s1, s2], [0.4, 0.6])

        d1 = DiscreteCPT([1,2], Dict(tuple(false) => [0.1, 0.9], tuple(true) => [0.2, 0.8]))
        d2 = DiscreteCPT([1,2], Dict(tuple(false) => [0.3, 0.7], tuple(true) => [0.4, 0.6]))
        mx2 = Mixture([d1, d2], [0.6, 0.4])

        d3 = DiscreteCPT([:a, :b], Dict((1,1) => [0.1, 0.9], (1,2) => [0.2, 0.8], (2,1) => [0.3, 0.7], (2,2) => [0.4, 0.6]))
        d4 = DiscreteCPT([:a, :b], Dict((1,1) => [0.9, 0.1], (1,2) => [0.8, 0.2], (2,1) => [0.7, 0.3], (2,2) => [0.6, 0.4]))
        mx3 = Mixture([d3, d4], [0.6, 0.4])

        vr = support(mx1, tuple(), 100, Vector{Bool}())
        @test support_quality(mx1, tuple()) == :CompleteSupport
        @test length(vr) == 2
        @test false in vr
        @test true in vr

        total = 1000
        n = 0
        for i in 1:total
            if sample(mx1, ())
                n += 1
            end
        end
        pt = 0.3 * 0.2 + 0.7 * 0.6
        @test isapprox(n / total, pt, atol = 0.05)

        p1givenf = 0.6 * 0.1 + 0.4 * 0.3
        p1givent = 0.6 * 0.2 + 0.4 * 0.4
        @test isapprox(logcpdf(mx1, tuple(), true), log(pt))
        @test isapprox(logcpdf(mx2, tuple(false), 1), log(p1givenf))
        @test isapprox(logcpdf(mx2, tuple(true), 1), log(p1givent))

        pi = compute_pi(mx2, [1,2], ([false, true],), (Cat([false, true], [0.7, 0.3]),))
        p1 = cpdf(pi, (), 1)
        p2 = cpdf(pi, (), 2)
        @test isapprox(p1, 0.7 * p1givenf + 0.3 * p1givent)
        @test isapprox(p2, 1 - p1)
        p1given11 = 0.6 * 0.1 + 0.4 * 0.9
        p1given12 = 0.6 * 0.2 + 0.4 * 0.8
        p1given21 = 0.6 * 0.3 + 0.4 * 0.7
        p1given22 = 0.6 * 0.4 + 0.4 * 0.6
        q1 = 0.3 * p1given11 + 0.7 * p1given12 # takes into account pi message from parent 2
        q2 = 0.3 * p1given21 + 0.7 * p1given22 # takes into account pi message from parent 2
        lam1 = send_lambda(mx3,
            SoftScore([:a, :b], [0.4, 0.3]), [:a, :b], ([1,2], [1,2]),
            (Cat([1,2], [0.1, 0.9]), Cat([1,2], [0.3, 0.7])), 1)
        @test isapprox(get_score(lam1, 1), 0.4 * q1 + 0.3 * (1 - q1))
        @test isapprox(get_score(lam1, 2), 0.4 * q2 + 0.3 * (1 - q2))
    end

    @testset "Separable" begin
        cpt1 = Dict((:x,) => [0.1, 0.9], (:y,) => [0.2, 0.8], (:z,) => [0.3, 0.7])
        cpt2 = Dict((1,) => [0.4, 0.6], (2,) => [0.5, 0.5])
        cpt3 = Dict(('a',) => [0.6, 0.4], ('b',) => [0.7, 0.3])
        s = Separable([false, true], [0.2, 0.3, 0.5], [cpt1, cpt2, cpt3])

        parranges = ([:x, :y, :z], [1, 2], ['a', 'b'])
        myrange = [false, true]
        @test support(s, parranges, 100, Bool[]) == myrange
        @test support_quality(s, parranges) == :CompleteSupport

        parvals = (:y, 1, 'b')
        c1 = Cat([:x, :y, :z], [0.2, 0.3, 0.5])
        c2 = Cat([1, 2], [0.1, 0.9])
        c3 = Cat(['a', 'b'], [0.2, 0.8])
        n = 0
        tot = 1000
        for i = 1:tot
            if sample(s, (:y, 1, 'b'))
                n += 1
            end
        end
        p = 0.2 * 0.8 + 0.3 * 0.6 + 0.5 * 0.3
        @test isapprox(n / tot, p; atol = 0.05)
        @test isapprox(logcpdf(s, parvals, true), log(p))

        k1 = nextkey()
        k2 = nextkey()
        k3 = nextkey()
        k4 = nextkey()
        (lfs,ufs) = make_factors(s, myrange, parranges, k1, (k2, k3, k4))
        @test length(ufs) == length(lfs)
        for i = 1:length(ufs)
            @test ufs[i].entries == lfs[i].entries
        end
        @test length(lfs) == 11 # each component has (num parent values + 1) factors, and one for the mixture
        mixfact = lfs[11]
        @test length(mixfact.keys) == 1
        mixkey = mixfact.keys[1]
        @test mixfact.dims == (3,)
        @test mixfact.entries == [0.2, 0.3, 0.5]
        comp1facts = lfs[1:4]
        comp2facts = lfs[5:7]
        comp3facts = lfs[8:10]
        comp1sw = last(comp1facts)
        comp2sw = last(comp2facts)
        comp3sw = last(comp3facts)
        @test length(comp1sw.keys) == 3
        @test length(comp2sw.keys) == 3
        @test length(comp3sw.keys) == 3
        @test comp1sw.keys[1] == k2
        @test comp2sw.keys[1] == k3
        @test comp3sw.keys[1] == k4
        sw1key = comp1sw.keys[2]
        sw2key = comp2sw.keys[2]
        sw3key = comp3sw.keys[2]
        @test comp1sw.keys[3] == mixkey
        @test comp2sw.keys[3] == mixkey
        @test comp3sw.keys[3] == mixkey
        @test comp1sw.dims == (3,3,3)
        @test comp2sw.dims == (2,2,3)
        @test comp3sw.dims == (2,2,3)
        @test comp1sw.entries ==
        [
            1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 0.0, 1.0, 1.0,
            0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0,
            0.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0
        ]
        @test comp2sw.entries ==
        [
            1.0, 1.0, 1.0, 1.0, 0.0, 1.0,
            1.0, 0.0, 1.0, 1.0, 1.0, 1.0
        ]
        @test comp3sw.entries ==
        [
            1.0, 1.0, 1.0, 1.0, 1.0, 0.0,
            1.0, 1.0, 0.0, 1.0, 1.0, 1.0
        ]
    end

    @testset "Switch" begin
        s1 = LinearSwitch(2, Symbol)
        s2 = If{Symbol}()

        prs1 = ([2,1], [:a, :b], [:c, :b]) # different order
        prs2 = ([false, true], [:a, :b], [:c, :b])
        range = [:c, :b, :a] # different order

        v1 = support(s1, prs1, 100, Symbol[])
        v2 = support(s2, prs2, 100, Symbol[])
        @test length(v1) == 3
        @test :a in v1
        @test :b in v1
        @test :c in v1
        @test length(v2) == 3
        @test :a in v2
        @test :b in v2
        @test :c in v2

        ns1 = Dict(:a => 0, :b => 0, :c => 0)
        ns2 = Dict(:a => 0, :b => 0, :c => 0)
        tot = 1000
        for j in 1:tot
            ns1[sample(s1, (2,:a,:b))] += 1
            ns2[sample(s2, (false,:a,:b))] += 1
        end
        @test ns1[:a] == 0
        @test ns1[:b] == tot
        @test ns1[:c] == 0
        @test ns2[:a] == 0
        @test ns2[:b] == tot
        @test ns2[:c] == 0

        @test isapprox(logcpdf(s1, (2,:a,:b), :a), -Inf)
        @test isapprox(logcpdf(s1, (2,:a,:b), :b), 0.0)

        @test support_quality(s1, prs1) == :CompleteSupport
        @test support_quality(s2, prs2) == :CompleteSupport
        @test support_quality(s1, ([2], [:a, :b], [:b, :c])) == :BestEffortSupport
        @test support_quality(s1, ([false], [:a, :b], [:b, :c])) == :BestEffortSupport

        incoming_pis1 = (Cat([2,1], [0.4, 0.6]), Cat([:a, :b], [0.1, 0.9]), Cat([:c, :b], [0.8, 0.2]))
        incoming_pis2 = (Cat([false, true], [0.4, 0.6]), Cat([:a, :b], [0.1, 0.9]), Cat([:c, :b], [0.8, 0.2]))
        pi1 = compute_pi(s1, range, prs1, incoming_pis1)
        pi2 = compute_pi(s2, range, prs2, incoming_pis2)
        p2ac = 0.4 * 0.1 * 0.8
        p2ab = 0.4 * 0.1 * 0.2
        p2bc = 0.4 * 0.9 * 0.8
        p2bb = 0.4 * 0.9 * 0.2
        p1ac = 0.6 * 0.1 * 0.8
        p1ab = 0.6 * 0.1 * 0.2
        p1bc = 0.6 * 0.9 * 0.8
        p1bb = 0.6 * 0.9 * 0.2
        pc = p2ac + p2bc
        pb = p2ab + p2bb + p1bc + p1bb
        pa = p1ac + p1ab
        test_pi(pi1, range, [pc, pb, pa])
        test_pi(pi2, range, [pc, pb, pa])

        chlam = SoftScore(range, [0.1, 0.2, 0.7])
        lam1 = send_lambda(s1, chlam, range, prs1, incoming_pis1, 1)
        lam2 = send_lambda(s1, chlam, range, prs1, incoming_pis1, 2)
        lam3 = send_lambda(s1, chlam, range, prs1, incoming_pis1, 3)
        @test Set(keys(lam1.logscores)) == Set(prs1[1])
        @test Set(keys(lam2.logscores)) == Set(prs1[2])
        @test Set(keys(lam3.logscores)) == Set(prs1[3])
        pac = 0.1 * 0.8 # ignoring first parent for lambda message to first parent
        pab = 0.1 * 0.2
        pbc = 0.9 * 0.8
        pbb = 0.9 * 0.2
        l2ac = pac * 0.1
        l2ab = pab * 0.2
        l2bc = pbc * 0.1
        l2bb = pbb * 0.2
        l1ac = pac * 0.7
        l1ab = pab * 0.7
        l1bc = pbc * 0.2
        l1bb = pbb * 0.2
        l2 = [l2ac + l2bc, l2ab + l2bb]
        l1 = [l1ac + l1ab, l1bc + l1bb]
        @test isapprox(get_score(lam1, 2), sum(l2))
        @test isapprox(get_score(lam1, 1), sum(l1))
        # Since this is a LinearSwitch, parent 2 is the choice for input 1, which is second in the first parent's range
        con2 = 0.4 * sum(l2)
        @test isapprox(get_score(lam2, :a), 0.6 * 0.7 + con2)
        @test isapprox(get_score(lam2, :b), 0.6 * 0.2 + con2)
        # And parent 3 is the choice for input 2, which is first in the first parent's range
        con3 = 0.6 * sum(l1)
        @test isapprox(get_score(lam3, :c), 0.4 * 0.1 + con3)
        @test isapprox(get_score(lam3, :b), 0.4 * 0.2 + con3)
    end

    @testset "Generate" begin
        frng1 = [Flip(0.1), Flip(0.2)]
        f1 = Flip(0.1)
        f2 = Flip(0.2)
        g = Generate{Bool}()
        vs = support(g, ([f1, f2],), 0, Bool[])
        @test support_quality(g, ([f1, f2],)) == :CompleteSupport
        @test length(vs) == 2
        @test Set(vs) == Set([false, true])
        @test isapprox(cpdf(g, (f1,), true), 0.1)
        total = 1000
        n = 0
        for i = 1:total
            if sample(g, (f1,))
                n += 1
            end
        end
        @test isapprox(n / total, 0.1, atol = 0.05)
    end

    @testset "Apply" begin
        frng1 = [Flip(0.1), Flip(0.2)]
        l1 = LinearGaussian((1.0,), 0.0, 1.0)
        l2 = LinearGaussian((1.0,), 1.0, 1.0)
        @test typeof(l2) <: SFunc{<:Tuple{Vararg{Float64}}, Float64}
        frng2 = [l1, l2]
        jrng2 = [(1.0,), (2.0,), (3.0,)]
        a1 = Apply{Tuple{}, Bool}()
        a2 = Apply{Tuple{Float64}, Float64}()
        vs = support(a1, (frng1, Vector{Tuple}[]), 0, Bool[])
        @test support_quality(a1, (frng1, Tuple[])) == :CompleteSupport
        @test support_quality(a2, (frng2, jrng2)) == :IncrementalSupport
        @test length(vs) == 2
        @test Set(vs) == Set([false, true])
        @test isapprox(logcpdf(a2, (l1, (1.0,)), 1.0), log(1 / sqrt(2 * pi)))
        total = 1000
        n = 0
        for i = 1:total
            if sample(a2, (l1, (1.0,))) < 1.0
                n += 1
            end
        end
        @test isapprox(n / total, 0.5, atol = 0.05)
    end

    @testset "Chain" begin
        @testset "With simple I and no J" begin
            sf = Chain(Int, Int, i -> Constant(i+1))
            @test sample(sf, (1,)) == 2
        end

        @testset "With tuple I and no J" begin
            sf = Chain(Tuple{Int}, Int, i -> Constant(i[1]+1))
            @test sample(sf, (1,)) == 2
        end

        @testset "With simple I and tuple J" begin
            sf = Chain(Int, Tuple{Int}, Int,
                i -> Det(Tuple{Int}, Int, j -> j[1] + i))
            @test sample(sf, (1,2)) == 3
        end
    end

    @testset "Invertible" begin
        i = Invertible{Int,Int}(i -> i + 1, o -> o - 1)
        @test support(i, ([1,2],), 100, Int[]) == [2,3]
        @test support_quality(i, ([1,2],)) == :CompleteSupport
        @test sample(i, (1,)) == 2
        @test cpdf(i, (1,), 2) == 1.0
        @test cpdf(i, (1,), 3) == 0.0
        ps = [1.0,0.0,0.0,0.0,1.0,0.0]
        (bps1, bps2) = Scruff.Operators.bounded_probs(i, [2,3,4], ([1,2],))
        @test bps1 == ps
        @test bps2 == bps1
        (facts1, facts2) = make_factors(i, [2,3,4], ([1,2],), 5, (7,))
        @test facts1 == facts2
        @test length(facts1) == 1
        fact1 = facts1[1]
        @test fact1.keys == (7,5)
        @test fact1.dims == (2,3)
        @test fact1.entries == ps
        parpis = (Cat([1,2], [0.1,0.9]),)
        chlam = SoftScore([2,3,4], [0.2,0.3,0.5])
        pi = compute_pi(i, [2,3,4], ([1,2],), parpis)
        @test cpdf(pi, (), 2) == 0.1
        @test cpdf(pi, (), 3) == 0.9
        @test cpdf(pi, (), 4) == 0.0
        lam = send_lambda(i, chlam, [2,3,4], ([1,2],), parpis, 1)
        @test get_score(lam, 1) == 0.2
        @test get_score(lam, 2) == 0.3
        @test get_score(lam, 3) == 0.0 # even though it maps to 4, it's not in the parent range
        #=
        @test initial_stats(i) == nothing
        @test accumulate_stats(i, nothing, nothing) == nothing
        @test expected_stats(i, [2,3,4], ([1,2],), parpis, chlam) == nothing
        @test maximize_stats(i, nothing) == nothing
        =#
    end

    @testset "Serial" begin
        sf1 = DiscreteCPT([:a, :b],
            Dict((1,1) => [0.1, 0.9], (1,2) => [0.2, 0.8],
                 (2,1) => [0.3, 0.7], (2,2) => [0.4, 0.6]))
        sf2 = DiscreteCPT([false, true], Dict((:a,) => [0.6, 0.4], (:b,) => [0.8, 0.2]))
        sf3 = Invertible{Bool, Int}(b -> b ? 5 : 6, i -> i == 5)
        ser = Serial(Tuple{Int,Int}, Int, (sf1,sf2,sf3))

        total = 1000
        n = 0
        for i in 1:total
            if sample(ser, (1,2)) == 6
                n += 1
            end
        end
        @test isapprox(n / total, 0.2 * 0.6 + 0.8 * 0.8; atol = 0.05)

        prs = ([1,2], [2,1])
        sup = support(ser, prs, 10, Int[])
        @test Set(sup) == Set([5,6])
        @test support_quality(ser, prs) == :CompleteSupport

        @test isapprox(cpdf(ser, (1,2), 6), 0.2 * 0.6 + 0.8 * 0.8)

        bps = [
                0.2 * 0.6 + 0.8 * 0.8, 0.2 * 0.4 + 0.8 * 0.2,
                0.1 * 0.6 + 0.9 * 0.8, 0.1 * 0.4 + 0.9 * 0.2,
                0.4 * 0.6 + 0.6 * 0.8, 0.4 * 0.4 + 0.6 * 0.2,
                0.3 * 0.6 + 0.7 * 0.8, 0.3 * 0.4 + 0.7 * 0.2
              ]
        @test isapprox(bounded_probs(ser, [6,5], prs)[1], bps)

        facts = make_factors(ser, [6,5], prs, 3, (1,2))[1]
        fs1 = make_factors(sf1, [:a,:b], prs, 4, (1,2))[1]
        fs2 = make_factors(sf2, [false, true], ([:a,:b],), 5, (4,))[1]
        fs3 = make_factors(sf3, [6,5], ([false, true],), 3, (5,))[1]
        @test length(facts) == length(fs1) + length(fs2) + length(fs3)
        fs = copy(fs1)
        append!(fs, fs2)
        append!(fs, fs3)
        for (computed, actual) in zip(facts, fs)
            @test computed.entries == actual.entries
        end

        parpis = (Cat([1,2], [0.9, 0.1]), Cat([2,1], [0.8, 0.2]))
        pi = compute_pi(ser, [6,5], prs, parpis)
        pi6 = 0.9 * 0.8 * bps[1] + 0.9 * 0.2 * bps[3] +
                0.1 * 0.8 * bps[5] + 0.1 * 0.2 * bps[7]
        pi5 = 0.9 * 0.8 * bps[2] + 0.9 * 0.2 * bps[4] +
                0.1 * 0.8 * bps[6] + 0.1 * 0.2 * bps[8]
        @test isapprox(cpdf(pi, (), 6), pi6)
        @test isapprox(cpdf(pi, (), 5), pi5)

        chlam = SoftScore([6,5], [0.7, 0.3])
        lam1 = send_lambda(ser, chlam, [6,5], prs, parpis, 1)
        lam2 = send_lambda(ser, chlam, [6,5], prs, parpis, 2)
        lf = 0.7
        lt = 0.3
        la = 0.6 * lf + 0.4 * lt
        lb = 0.8 * lf + 0.2 * lt
        l11 = 0.8 * (0.2 * la + 0.8 * lb) + 0.2 * (0.1 * la + 0.9 * lb)
        l12 = 0.8 * (0.4 * la + 0.6 * lb) + 0.2 * (0.3 * la + 0.7 * lb)
        @test isapprox(get_score(lam1, 1), l11)
        @test isapprox(get_score(lam1, 2), l12)
        l22 = 0.9 * (0.2 * la + 0.8 * lb) + 0.1 * (0.4 * la + 0.6 * lb)
        l21 = 0.9 * (0.1 * la + 0.9 * lb) + 0.1 * (0.3 * la + 0.7 * lb)
        @test isapprox(get_score(lam2, 1), l21)
        @test isapprox(get_score(lam2, 2), l22)

        #=
        is1 = initial_stats(sf1)
        is2 = initial_stats(sf2)
        is3 = initial_stats(sf3)
        istats = initial_stats(ser)
        @test istats == (is1,is2,is3)

        es1 = expected_stats(sf1, [:a,:b], prs, parpis, SoftScore([:a,:b], [la,lb]))
        pi1 = compute_pi(sf1, [:a,:b], prs, parpis)
        es2 = expected_stats(sf2, [false, true], ([:a,:b],), (pi1,),
                SoftScore([false, true], [lf,lt]))
        pi2 = compute_pi(sf2, [false, true], ([:a,:b],), (pi1,))
        es3 = expected_stats(sf3, [6,5], ([false, true],), (pi2,), chlam)
        estats = expected_stats(ser, [6,5], prs, parpis, chlam)
        @test length(estats) == 3
        @test keys(estats[1]) == keys(es1)
        for k in keys(estats[1])
            @test isapprox(estats[1][k], es1[k])
        end
        @test keys(estats[2]) == keys(es2)
        for k in keys(estats[2])
            @test isapprox(estats[2][k], es2[k])
        end
        @test estats[3] == es3

        as1 = accumulate_stats(sf1, is1, es1)
        as2 = accumulate_stats(sf2, is2, es2)
        as3 = accumulate_stats(sf3, is3, es3)
        astats = accumulate_stats(ser, istats, estats)
        @test length(astats) == 3
        @test keys(astats[1]) == keys(as1)
        for k in keys(astats[1])
            @test isapprox(astats[1][k], as1[k])
        end
        @test keys(astats[2]) == keys(as2)
        for k in keys(astats[2])
            @test isapprox(astats[2][k], as2[k])
        end
        @test astats[3] == as3

        mp1 = maximize_stats(sf1, as1)
        mp2 = maximize_stats(sf2, as2)
        mp3 = maximize_stats(sf3, as3)
        mparams = maximize_stats(ser, astats)
        @test length(mparams) == 3
        @test isapprox(mparams[1], mp1)
        @test isapprox(mparams[2], mp2)
        @test mparams[3] == mp3
        =#
    end

end


text:
# %% Import packages
using Plots, Flux

# %% The training data X and the output y
X = [0f0 0 1 1; 0 1 0 1]
y = [0f0 1 1 0]

# Our model has 2 input, 2 neurons in hidden layer and 1 output. The activation function is a sigmoid
xor_model = Chain(
    Dense(2, 2, sigmoid),
    Dense(2, 1, sigmoid)
)

# Check our model
params(xor_model)

# Check our inputs and output
params(xor_model)[1]
params(xor_model)[3]

# Implements a loss functions using mse
loss_fn(a, b) = Flux.mse(xor_model(a), b)

# We will use ADAM optimizer
opt = ADAM(0.001)

# Number of Epochs to run
N = 5000
loss = zeros(N)

# Train the model
for i in 1:N
    Flux.train!(loss_fn, params(xor_model), [(X, y)], opt)
    loss[i] = loss_fn(X, y)
    if i % 100 == 0
        println("loss $loss[i]")
    end
end

# %% Print loss and predict our data
loss_fn(X, y)
xor_model(X)

# Plot the loss data
loss_plt = plot(1:N, loss, xlabel="Epochs", legend=:none, ylabel="Loss (mse)");
display("image/png", loss_plt)

max_stars_repo_name: jcbritobr/flux-xor-ml

max_issues_repo_name: jcbritobr/flux-xor-ml

max_forks_repo_name: jcbritobr/flux-xor-ml


text:
import Compat.LinearAlgebra: eigvals

eigvals(A::QuMatrix, args...; kwargs...) =
    eigvals(A.data, args...; kwargs...)

function eigen(A::QuMatrix, args...; kwargs...)
    D,V = eigen(A.data, args...; kwargs...)
    return D, [Ket(V[:,i],dims(A)) for i = 1:size(V,2)]
end

function eigs(A::Operator, args...; kwargs...)
    D,V = eigs(A.data, args...; which=:SR, kwargs...)
    return D, [Ket(V[:,i],dims(A)) for i = 1:size(V,2)]
end


text:
function formatinfo(s::Spec, x::Complex)
    s′ = Spec(s, width = nothing, sign = SIGN_NONE)
    size1, meta1 = formatinfo(s′, real(x))
    size2, meta2 = formatinfo(s′, imag(x))
    size = size1 + size2 + 5  # " + " and "im"
    return size, (meta1, meta2)
end

function formatfield(data::Vector{UInt8}, p::Int, s::Spec, x::Complex, (meta1, meta2))
    start = p
    p, _ = sign(data, p, real(x), s.sign)
    s′ = Spec(s, width = nothing, sign = SIGN_NONE)
    p = formatfield(data, p, s′, real(x), meta1)
    if signbit(imag(x))
        p = @copy data p " - "
    else
        p = @copy data p " + "
    end
    p = formatfield(data, p, s′, imag(x), meta2)
    p = @copy data p "im"
    if isspecified(s.width)
        p = aligncontent(data, p, start, p - start, s.fill, default(s.align, ALIGN_RIGHT), s.width)
    end
    return p
end


text:
module TexTables

# package code goes here

# Nice string formattting
using Formatting, DataStructures, DataFrames

# Required for StatsModels Integration
using StatsBase
using Distributions
using GLM: LinearModel
using Parameters
using Compat

export FormattedNumber, FNum, FNumSE, @fmt, TableCol, star!
export TableCol, Table, TexTable, get_vals
export IndexedTable, append_table, join_table, promote_rule
export to_tex, to_ascii, write_tex, regtable
export setcoef!, setmeta!, setstats!, RegCol

# Import from base to extend
import Base.getindex, Base.setindex!, Base.push!
import Base: isless, ==
import Base: getindex, size, hcat, vcat, hvcat, convert, promote_rule
import Base: show, size, print

export summarize, summarize_by, tabulate

# Need this abstract type
abstract type TexTable end

include("FormattedNumbers.jl")
include("TableCol.jl")
include("CompositeTable.jl")
include("Printing.jl")
include("StatsModelsInterface.jl")
include("QuickTools.jl")
end # module