text:
struct KineticMixing <: AbstractDMModel
    mx::Float64
    mv::Float64
    gvxx::Float64
    eps::Float64
    width_v::Float64
end

function KineticMixing(mx::Float64, mv::Float64, gvxx::Float64, eps::Float64)
    model = KineticMixing(mx, mv, gvxx, eps, 0.0)
    width = width_v(model)
    KineticMixing(mx, mv, gvxx, eps, width)
end

dm_mass(model::KineticMixing) = model.mx
dm_dof(::KineticMixing) = 2.0
dm_spin2(::KineticMixing) = 2

# =================================
# Partial widths of vector mediator
# =================================

"""
    width_v_xx(model::KineticMixing)

Compute the partial width of the vector mediator into dark matter.
"""
function width_v_xx(model::KineticMixing)
    2model.mx >= model.mv && return 0.0
    (
        model.gvxx^2 *
        sqrt(-4 * model.mx^2 + model.mv^2) *
        (2 * model.mx^2 + model.mv^2)
    ) / (12 * model.mv^2 * pi)
end

"""
    width_v_νν(model::KineticMixing)

Compute the partial width of the vector mediator into neutrinos.
"""
function width_v_νν(model::KineticMixing)
    (ALPHA_EM * model.eps^2 * model.mv) / (24 * COS_THETA_WEAK^2)
end

"""
    width_v_ll(ml::Float64, model::KineticMixing)

Compute the partial width of the vector mediator into generic leptons.
"""
function width_v_ll(ml::Float64, model::KineticMixing)
    2ml >= model.mv && return 0.0
    (
        ALPHA_EM *
        model.eps^2 *
        sqrt(-4 * ml^2 + model.mv^2) *
        (7 * ml^2 + 5 * model.mv^2)
    ) / (24 * COS_THETA_WEAK^2 * model.mv^2)
end

"""
    width_v_ee(model::KineticMixing)

Compute the partial width of the vector mediator into electrons.
"""
width_v_ee(model::KineticMixing) = width_v_ll(ELECTRON_MASS, model)

"""
    width_v_μμ(model::KineticMixing)

Compute the partial width of the vector mediator into muons.
"""
width_v_μμ(model::KineticMixing) = width_v_ll(MUON_MASS, model)

"""
    width_v_ττ(model::KineticMixing)

Compute the partial width of the vector mediator into taus.
"""
width_v_ττ(model::KineticMixing) = width_v_ll(TAU_MASS, model)

"""
    width_v_ququ(model::KineticMixing)

Compute the partial width of the vector mediator into generic up-type quarks.
"""
function width_v_ququ(mu::Float64, model::KineticMixing)
    2mu >= model.mv && return 0.0
    (
        ALPHA_EM *
        model.eps^2 *
        sqrt(-4 * mu^2 + model.mv^2) *
        (7 * mu^2 + 17 * model.mv^2)
    ) / (72 * COS_THETA_WEAK^2 * model.mv^2)
end

"""
    width_v_uu(model::KineticMixing)

Compute the partial width of the vector mediator into generic up quarks.
"""
width_v_uu(model::KineticMixing) = width_v_ququ(UP_QUARK_MASS, model)

"""
    width_v_cc(model::KineticMixing)

Compute the partial width of the vector mediator into generic charms quarks.
"""
width_v_cc(model::KineticMixing) = width_v_ququ(CHARM_QUARK_MASS, model)

"""
    width_v_tt(model::KineticMixing)

Compute the partial width of the vector mediator into generic top quarks.
"""
width_v_tt(model::KineticMixing) = width_v_ququ(TOP_QUARK_MASS, model)

"""
    width_v_qdqd(model::KineticMixing)

Compute the partial width of the vector mediator into generic down-type quarks.
"""
function width_v_qdqd(md::Float64, model::KineticMixing)
    2md >= model.mv && return 0.0
    (
        ALPHA_EM *
        model.eps^2 *
        sqrt(-4 * md^2 + model.mv^2) *
        (-17 * md^2 + 5 * model.mv^2)
    ) / (72 * COS_THETA_WEAK^2 * model.mv^2)
end

"""
    width_v_dd(model::KineticMixing)

Compute the partial width of the vector mediator into generic down quarks.
"""
width_v_dd(model::KineticMixing) = width_v_qdqd(DOWN_QUARK_MASS, model)

"""
    width_v_ss(model::KineticMixing)

Compute the partial width of the vector mediator into generic strange quarks.
"""
width_v_ss(model::KineticMixing) = width_v_qdqd(STRANGE_QUARK_MASS, model)

"""
    width_v_bb(model::KineticMixing)

Compute the partial width of the vector mediator into generic bottom quarks.
"""
width_v_bb(model::KineticMixing) = width_v_qdqd(BOTTOM_QUARK_MASS, model)

"""
    width_v_hz(model::KineticMixing)

Compute the partial width of the vector mediator into generic Higgs + Z.
"""
function width_v_hz(model::KineticMixing)
    HIGGS_MASS + Z_BOSON_MASS >= model.mv && return 0.0

    (
        ALPHA_EM^2 *
        model.eps^2 *
        sqrt(
            -HIGGS_MASS^2 +
            (HIGGS_MASS^2 + model.mv^2 - Z_BOSON_MASS^2)^2 / (4 * model.mv^2),
        ) *
        (
            2 * model.mv^2 * Z_BOSON_MASS^2 +
            (-HIGGS_MASS^2 + model.mv^2 + Z_BOSON_MASS^2)^2 / 4.0
        ) *
        pi *
        HIGGS_VEV^2
    ) / (6 * COS_THETA_WEAK^4 * model.mv^4 * Z_BOSON_MASS^2 * SIN_THETA_WEAK^2)
end

"""
    width_v(model::KineticMixing)

Compute the total width of the vector mediator.
"""
function width_v(model::KineticMixing)
    (
        width_v_xx(model) +
        width_v_ee(model) +
        width_v_μμ(model) +
        width_v_ττ(model) +
        width_v_uu(model) +
        width_v_cc(model) +
        width_v_tt(model) +
        width_v_dd(model) +
        width_v_ss(model) +
        width_v_tt(model) +
        width_v_hz(model) +
        3width_v_νν(model)
    )
end

# =========================================================================
# ==== Cross sections =====================================================
# =========================================================================

"""
    cross_section_xx_νν(Q::Real, model::KineticMixing)

Compute cross section for DM annihilating into neutrinos at
center-of-mass energy `Q`.
"""
function cross_section_xx_νν(Q::Real, model::KineticMixing)
    Q <= 2model.mx && return zero(typeof(Q))
    (
        ALPHA_EM *
        model.eps^2 *
        model.gvxx^2 *
        sqrt(Q^2) *
        (2 * model.mx^2 + Q^2)
    ) / (
        24 *
        COS_THETA_WEAK^2 *
        sqrt(-4 * model.mx^2 + Q^2) *
        ((model.mv^2 - Q^2)^2 + model.mv^2 * model.width_v^2)
    )
end

"""
    cross_section_xx_ll(Q::Real, model::KineticMixing)

Compute cross section for DM annihilating into generic leptons at
center-of-mass energy `Q`.
"""
function cross_section_xx_ll(Q::Real, ml::Float64, model::KineticMixing)
    (Q <= 2model.mx || Q <= 2ml) && return zero(typeof(Q))
    (
        ALPHA_EM *
        model.eps^2 *
        model.gvxx^2 *
        (2 * model.mx^2 + Q^2) *
        sqrt(-4 * ml^2 + Q^2) *
        (7 * ml^2 + 5 * Q^2)
    ) / (
        24 *
        COS_THETA_WEAK^2 *
        Q^2 *
        sqrt(-4 * model.mx^2 + Q^2) *
        ((model.mv^2 - Q^2)^2 + model.mv^2 * model.width_v^2)
    )
end

"""
    cross_section_xx_ee(Q::Real, model::KineticMixing)

Compute cross section for DM annihilating into electons at
center-of-mass energy `Q`.
"""
cross_section_xx_ee(Q::Real, model::KineticMixing) =
    cross_section_xx_ll(Q, ELECTRON_MASS, model)

"""
    cross_section_xx_μμ(Q::Real, model::KineticMixing)

Compute cross section for DM annihilating into muon at
center-of-mass energy `Q`.
"""
cross_section_xx_μμ(Q::Real, model::KineticMixing) =
    cross_section_xx_ll(Q, MUON_MASS, model)

"""
    cross_section_xx_ττ(Q::Real, model::KineticMixing)

Compute cross section for DM annihilating into taus at
center-of-mass energy `Q`.
"""
cross_section_xx_ττ(Q::Real, model::KineticMixing) =
    cross_section_xx_ll(Q, TAU_MASS, model)

"""
    cross_section_xx_ququ(Q::Real, model::KineticMixing)

Compute cross section for DM annihilating into generic up-type quarks at
center-of-mass energy `Q`.
"""
function cross_section_xx_ququ(Q::Real, mu::Float64, model::KineticMixing)
    (Q <= 2model.mx || Q <= 2mu) && return zero(typeof(Q))
    (
        ALPHA_EM *
        model.eps^2 *
        model.gvxx^2 *
        (2 * model.mx^2 + Q^2) *
        sqrt(-4 * mu^2 + Q^2) *
        (7 * mu^2 + 17 * Q^2)
    ) / (
        72 *
        COS_THETA_WEAK^2 *
        Q^2 *
        sqrt(-4 * model.mx^2 + Q^2) *
        ((model.mv^2 - Q^2)^2 + model.mv^2 * model.width_v^2)
    )
end

"""
    cross_section_xx_uu(Q::Real, model::KineticMixing)

Compute cross section for DM annihilating into up quarks at
center-of-mass energy `Q`.
"""
cross_section_xx_uu(Q::Real, model::KineticMixing) =
    cross_section_xx_ququ(Q, UP_QUARK_MASS, model)

"""
    cross_section_xx_cc(Q::Real, model::KineticMixing)

Compute cross section for DM annihilating into charm quarks at
center-of-mass energy `Q`.
"""
cross_section_xx_cc(Q::Real, model::KineticMixing) =
    cross_section_xx_ququ(Q, CHARM_QUARK_MASS, model)

"""
    cross_section_xx_tt(Q::Real, model::KineticMixing)

Compute cross section for DM annihilating into top quarks at
center-of-mass energy `Q`.
"""
cross_section_xx_tt(Q::Real, model::KineticMixing) =
    cross_section_xx_ququ(Q, TOP_QUARK_MASS, model)

"""
    cross_section_xx_qdqd(Q::Real, model::KineticMixing)

Compute cross section for DM annihilating into generic down-type quarks at
center-of-mass energy `Q`.
"""
function cross_section_xx_qdqd(Q::Real, md::Float64, model::KineticMixing)
    (Q <= 2model.mx || Q <= 2md) && return zero(typeof(Q))
    (
        ALPHA_EM *
        model.eps^2 *
        model.gvxx^2 *
        sqrt(-4 * md^2 + Q^2) *
        (2 * model.mx^2 + Q^2) *
        (-17 * md^2 + 5 * Q^2)
    ) / (
        72 *
        COS_THETA_WEAK^2 *
        Q^2 *
        sqrt(-4 * model.mx^2 + Q^2) *
        ((model.mv^2 - Q^2)^2 + model.mv^2 * model.width_v^2)
    )
end

"""
    cross_section_xx_dd(Q::Real, model::KineticMixing)

Compute cross section for DM annihilating into down quarks at
center-of-mass energy `Q`.
"""
cross_section_xx_dd(Q::Real, model::KineticMixing) =
    cross_section_xx_qdqd(Q, DOWN_QUARK_MASS, model)

"""
    cross_section_xx_ss(Q::Real, model::KineticMixing)

Compute cross section for DM annihilating into strange quarks at
center-of-mass energy `Q`.
"""
cross_section_xx_ss(Q::Real, model::KineticMixing) =
    cross_section_xx_qdqd(Q, STRANGE_QUARK_MASS, model)

"""
    cross_section_xx_bb(Q::Real, model::KineticMixing)

Compute cross section for DM annihilating into bottom quarks at
center-of-mass energy `Q`.
"""
cross_section_xx_bb(Q::Real, model::KineticMixing) =
    cross_section_xx_qdqd(Q, BOTTOM_QUARK_MASS, model)

"""
    cross_section_xx_hz(Q::Real, model::KineticMixing)

Compute cross section for DM annihilating into Higgs + Z at
center-of-mass energy `Q`.
"""
function cross_section_xx_hz(Q::Real, model::KineticMixing)
    (Q <= 2model.mx || Q <= HIGGS_MASS + Z_BOSON_MASS) && return zero(typeof(Q))
    (
        ALPHA_EM^2 *
        model.eps^2 *
        model.gvxx^2 *
        pi *
        sqrt(
            (
                (HIGGS_MASS - Z_BOSON_MASS - Q) *
                (HIGGS_MASS + Z_BOSON_MASS - Q) *
                (HIGGS_MASS - Z_BOSON_MASS + Q) *
                (HIGGS_MASS + Z_BOSON_MASS + Q)
            ) / (-4 * model.mx^2 + Q^2),
        ) *
        (2 * model.mx^2 + Q^2) *
        (
            HIGGS_MASS^4 + Z_BOSON_MASS^4 + 10 * Z_BOSON_MASS^2 * Q^2 + Q^4 -
            2 * HIGGS_MASS^2 * (Z_BOSON_MASS^2 + Q^2)
        ) *
        HIGGS_VEV^2
    ) / (
        48 *
        COS_THETA_WEAK^4 *
        Z_BOSON_MASS^2 *
        Q^5 *
        SIN_THETA_WEAK^2 *
        ((model.mv^2 - Q^2)^2 + model.mv^2 * model.width_v^2)
    )
end

"""
    cross_section_xx_vv(Q::Real, model::KineticMixing)

Compute cross section for DM annihilating into vector mediator at
center-of-mass energy `Q`.
"""
function cross_section_xx_vv(Q::Real, model::KineticMixing)
    (Q <= 2model.mx || Q <= 2model.mv) && return zero(typeof(Q))
    (
        model.gvxx^4 * (
            (
                -2 *
                sqrt(-4 * model.mx^2 + Q^2) *
                sqrt(-4 * model.mv^2 + Q^2) *
                (4 * model.mx^4 + 2 * model.mv^4 + model.mx^2 * Q^2)
            ) / (model.mv^4 + model.mx^2 * (-4 * model.mv^2 + Q^2)) +
            (2 * model.mx^2 - 2 * model.mv^2 + Q^2) * log(
                (
                    -2 * model.mv^2 +
                    Q^2 +
                    sqrt(-4 * model.mx^2 + Q^2) * sqrt(-4 * model.mv^2 + Q^2)
                )^2 /
                (
                    2 * model.mv^2 - Q^2 +
                    sqrt(-4 * model.mx^2 + Q^2) * sqrt(-4 * model.mv^2 + Q^2)
                )^2,
            ) +
            (
                2 *
                (
                    -4 * model.mx^4 +
                    2 * model.mv^2 * Q^2 +
                    model.mx^2 * (-2 * model.mv^2 + Q^2)
                ) *
                log(
                    (
                        -2 * model.mv^2 +
                        Q^2 +
                        sqrt(-4 * model.mx^2 + Q^2) *
                        sqrt(-4 * model.mv^2 + Q^2)
                    )^2 /
                    (
                        2 * model.mv^2 - Q^2 +
                        sqrt(-4 * model.mx^2 + Q^2) *
                        sqrt(-4 * model.mv^2 + Q^2)
                    )^2,
                )
            ) / (-2 * model.mv^2 + Q^2)
        )
    ) / (16 * pi * Q^2 * (-4 * model.mx^2 + Q^2))
end

"""
    dm_annihilation_cross_section(Q::Real, model::KineticMixing)

Compute total annihilation cross section for DM center-of-mass energy `Q`.
"""
function dm_annihilation_cross_section(Q::Real, model::KineticMixing)
    (
        3cross_section_xx_νν(Q, model) +
        cross_section_xx_ee(Q, model) +
        cross_section_xx_μμ(Q, model) +
        cross_section_xx_ττ(Q, model) +
        cross_section_xx_uu(Q, model) +
        cross_section_xx_cc(Q, model) +
        cross_section_xx_tt(Q, model) +
        cross_section_xx_dd(Q, model) +
        cross_section_xx_ss(Q, model) +
        cross_section_xx_bb(Q, model) +
        cross_section_xx_hz(Q, model) +
        cross_section_xx_vv(Q, model)
    )
end

"""
    dm_thermal_cross_section(Q::Float64, model::KineticMixing)

Compute thermally-averaged annihilation cross section for DM at a scaled
temperatue `x=mx/T`
"""
function dm_thermal_cross_section(x::Float64, model::KineticMixing)
    den = 2 * besselkx(2, x)
    pf = x / den^2

    integrand(z::Real) = (
        dm_annihilation_cross_section(model.mx * z, model) *
        z^2 *
        (z^2 - 4) *
        besselkx(1, x * z) *
        exp(-x * (z - 2))
    )

    return pf * quadgk(integrand, 2, Inf)[1]
end

"""
    dm_thermal_cross_section(Q::Real, model::KineticMixing)

Compute thermally-averaged annihilation cross section for DM at a scaled
temperatue `x=mx/T`
"""
function dm_thermal_cross_section(x::Real, model::KineticMixing)
    den = 2 * besselk(2, x)
    pf = x / den^2

    integrand(z::Real) = (
        dm_annihilation_cross_section(model.mx * z, model) *
        z^2 *
        (z^2 - 4) *
        besselk(1, x * z)
    )
    return pf * quadgk(integrand, 2, Inf)[1]
end

max_stars_repo_name: LoganAMorrison/AsymptoticBoltzmann

max_issues_repo_name: LoganAMorrison/AsymptoticBoltzmann

max_forks_repo_name: LoganAMorrison/AsymptoticBoltzmann


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).
using BAT
using Test
using Plots
using StatsBase

@testset "bin centers & local modes" begin

    data1 = [1, 2, 3, 4, 5, 5, 6, 7, 8, 9, 9, 10]
    hist_1d = fit(Histogram, data1, nbins = 10, closed = :left)
    bathist_1d = BATHistogram(hist_1d)

    @testset "1D BATHistogram" begin
        @test BAT.get_bin_centers(bathist_1d) == [collect(1.5:1:10.5)]
        @test BAT.find_marginalmodes(bathist_1d) == [[5.5], [9.5]]
    end

    data2 = [10, 20, 30, 40, 50, 50, 60, 70, 80, 90, 90, 100]
    hist_2d = fit(Histogram, (data1, data2), nbins = 10, closed = :left)
    bathist_2d = BATHistogram(hist_2d)

    @testset "2D BATHistogram" begin
        @test BAT.get_bin_centers(bathist_2d) == [collect(1.5:1:10.5), collect(15:10:105)]
        @test BAT.get_bin_centers(bathist_2d.h) == [collect(1.5:1:10.5), collect(15:10:105)]
        @test BAT.find_marginalmodes(bathist_2d) == [[5.5, 55.0], [9.5, 95.0]]
    end
end

ext: jl


text:
module InverseFieldsTests

using Gridap.Fields
using Gridap.Arrays
using Gridap.TensorValues

using Test

b0 = Point(0,0)
m0 = TensorValue(1,0,0,1)
id = AffineMap(m0,b0)

b1 = Point(1,1)
m1 = TensorValue(2,0,0,2)
h1 = AffineMap(m1,b1)

b2 = Point(3,3)
m2 = TensorValue(4,0,0,4)
h2 = AffineMap(m2,b2)

h = h2 ∘ h1

# (4x+3) ∘ (2x+1) = 8x+7
b3 = Point(7,7)
m3 = TensorValue(8,0,0,8)
h3 = AffineMap(m3,b3)

h1inv = inverse_map(h1)
h2inv = inverse_map(h2)
hinv = inverse_map(h)
hinvinv = inverse_map(hinv)

function fields_are_equal(f::Field, g::Field)
    local xs = Point{2,Float64}[Point(0,0), Point(1,0), Point(0,1)]
    # all(f(x) ≈ g(x) for x in xs)
    all(f(xs) ≈ g(xs))
end

@test fields_are_equal(id ∘ id, id)

@test fields_are_equal(h, h3)

@test fields_are_equal(id ∘ h, h)
@test fields_are_equal(h ∘ id, h)

@test fields_are_equal(hinv ∘ h, id)
@test fields_are_equal(h ∘ hinv, id)

@test fields_are_equal(hinvinv, h)

@test fields_are_equal(hinv, h1inv ∘ h2inv)

end # module


text:
# Inspired by
# - https://gitlab.com/tkpapp/GitlabJuliaDemo.jl by Tamas K. Papp
# - https://github.com/tpapp/LocalCoverage.jl by Tamas K. Papp

using Coverage

const report_dir = "coverage_report"
const lcov_info_file = "lcov.info"

# Change path to root directory
cd(joinpath(@__DIR__, "..", "..")) do
  # Process coverage files
  processed = process_folder("src")

  # Uncomment the following line once Codecov support is enabled
  # Codecov.submit_local(processed)

  # Calculate coverage
  covered_lines, total_lines = get_summary(processed)
  percentage = covered_lines / total_lines * 100

  # Print coverage in a format that can be easily parsed
  println("($(percentage)%) covered")

  # Try to generate a coverage report
  isdir(report_dir) || mkdir(report_dir)
  tracefile = joinpath(report_dir, lcov_info_file)
  Coverage.LCOV.writefile(tracefile, processed)
  branch = strip(read(`git rev-parse --abbrev-ref HEAD`, String))
  commit = strip(read(`git rev-parse --short HEAD`, String))
  title = "commit $(commit) on branch $(branch)"
  run(`genhtml -t $(title) -o $(report_dir) $(tracefile)`)

  # Clean up .cov files
  clean_folder("src")
end


text:
import ..UncertainDatasets:
    AbstractUncertainValueDataset

import ..UncertainStatistics: mean, std, median, quantile

@recipe function f(udata::AbstractUncertainValueDataset, quants::Vector{Float64} = [0.33, 0.67])
    n_points = length(udata)

    for i = 1:n_points
        med = median(udata[i], 10000)
        lower = quantile(udata[i], minimum(quants), 10000)
        upper = quantile(udata[i], maximum(quants), 10000)

        @series begin
            seriescolor --> :black
            label --> ""
            yerror --> ([med - lower], [upper - med])
            [i], [med]

        end
    end

end


text:
"""
    parse_onda("name-of-calibration-file.txt")

Translates calibration text files from Onda
into the format required by this package.
"""
function parse_onda(name)
    file = readlines(name)
    info = Dict()
    count = 0

    function found(s, i)
        return occursin(s, i) && !occursin("NONE", i)
    end

    function keep(s, i)
        if found(s, i)
            row = split(i)
            if length(row) == 2
                info[s] = row[2]
            else
                info[s] = row
            end
        end
    end

    for i in file
        count += 1
        keep("MFG", i)
        keep("MODEL", i)
        keep("SN", i)
        if occursin("HEADER_END", i)
            break
        end
    end

    info["DATA"] = file[count + 1:end]
    device = "$(info["MFG"])_$(info["MODEL"])_$(info["SN"])"
    device = replace(device, "-"=>"_")
    open("$device.jl", "w") do io
       println(io, "@calibration :$device [")
       for i in info["DATA"]
           println(io, i)
       end
       println(io, "]")
    end
end


text:
module GeneticAlgorithms

# -------

importall Base

export  Entity,
        GAmodel,

        runga,
        freeze,
        defrost,
        generation_num,
        population

# -------

abstract Entity

isless(lhs::Entity, rhs::Entity) = lhs.fitness < rhs.fitness

fitness!(ent::Entity, fitness_score) = ent.fitness = fitness_score

# -------

type EntityData
    entity
    generation::Int

    EntityData(entity, generation::Int) = new(entity, generation)
    EntityData(entity, model) = new(entity, model.gen_num)
end

# -------

type GAmodel
    initial_pop_size::Int
    gen_num::Int

    population::Array
    pop_data::Array{EntityData}
    freezer::Array{EntityData}

    rng::AbstractRNG

    ga

    GAmodel() = new(0, 1, Any[], EntityData[], EntityData[], MersenneTwister(time_ns()), nothing)
end

global _g_model

# -------

function freeze(model::GAmodel, entity::EntityData)
    push!(model.freezer, entity)
    println("Freezing: ", entity)
end

function freeze(model::GAmodel, entity)
    entitydata = EntityData(entity, model.gen_num)
    freeze(model, entitydata)
end

freeze(entity) = freeze(_g_model, entity)

function defrost(model::GAmodel, generation::Int)
    filter(model.freezer) do entitydata
        entitydata.generation == generation
    end
end

defrost(generation::Int) = defrost(_g_model, generation)

generation_num(model::GAmodel = _g_model) = model.gen_num

population(model::GAmodel = _g_model) = model.population

function runga(mdl::Module; initial_pop_size = 128)
    model = GAmodel()
    model.ga = mdl
    model.initial_pop_size = initial_pop_size

    runga(model)
end

function runga(model::GAmodel)
    reset_model(model)
    create_initial_population(model)

    while true
        evaluate_population(model)

        grouper = @task model.ga.group_entities(model.population)
        groupings = Any[]
        while !istaskdone(grouper)
            group = consume(grouper)
            group != nothing && push!(groupings, group)
        end

        if length(groupings) < 1
            break
        end

        crossover_population(model, groupings)
        mutate_population(model)
    end

    model
end

# -------

function reset_model(model::GAmodel)
    global _g_model = model

    model.gen_num = 1
    empty!(model.population)
    empty!(model.pop_data)
    empty!(model.freezer)
end

function create_initial_population(model::GAmodel)
    for i = 1:model.initial_pop_size
        entity = model.ga.create_entity(i)

        push!(model.population, entity)
        push!(model.pop_data, EntityData(entity, model.gen_num))
    end
end

function evaluate_population(model::GAmodel)
    scores = pmap(model.ga.fitness, model.population)
    for i in 1:length(scores)
        fitness!(model.population[i], scores[i])
    end

    sort!(model.population; rev = true)
end

function crossover_population(model::GAmodel, groupings)
    old_pop = model.population

    model.population = Any[]
    sizehint!(model.population, length(old_pop))

    model.pop_data = EntityData[]
    sizehint!(model.pop_data, length(old_pop))

    model.gen_num += 1

    for group in groupings
        parents = Any[ old_pop[i] for i in group ]
        entity = model.ga.crossover(parents)

        push!(model.population, entity)
        push!(model.pop_data, EntityData(entity, model.gen_num))
    end
end

function mutate_population(model::GAmodel)
    for entity in model.population
        model.ga.mutate(entity)
    end
end

end


text:
using openbuildnet

# Create the node
node = OBNNode("extnode", "testext", length(ARGS)>0?ASCIIString(ARGS[1]):"tcp://localhost:1883")

# Create the ports
udouble = create_input(node, "input_scalardouble", Float64)
ydouble = create_output(node, "output_scalardouble", Float64)

uvector = create_input(node, "input_vectordouble", Vector{Float64})
yvector = create_output(node, "output_vectordouble", Vector{Float64})

umatrix = create_input(node, "input_matrixdouble", Matrix{Float64})
ymatrix = create_output(node, "output_matrixdouble", Matrix{Float64})

# Create Callbacks
on_init(node) do
  println("INIT")
end

on_term(node) do
  println("TERMINATE")
end

const MAINBLOCK = 0
d = 1.0
v = [2.1, 0.9]
m = [1.0 2.0; 3.0 4.0]

on_block_output(node, MAINBLOCK) do
  println("Current time is ", sim_time(node), "s.")
  #println("Current time is ", wallclock_time(node))
  println("Scalar double: ", get(udouble))
  println("Vector double: ", get(uvector))
  println("Matrix double: ", get(umatrix))

  global d
  set(ydouble, d)
  d += 1.0

  global v
  set(yvector, v)
  v += 0.5

  global m
  set(ymatrix, m)
  m += eye(2)

  schedule(node, update_mask(MAINBLOCK), sim_time(node)+0.5, :s, 1.0)
end

# Test the optional parameters given to callbacks
type TheState
  x::Float64
end
mystate = TheState(0.0)

on_block_state(node, MAINBLOCK, mystate) do s
  s.x += get(udouble)
  println("My new state: ", s.x)
end

println("Ready to start simulation...")
status = run(node, 20)
println("Final status: ", status)


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# McCormick.jl
# A McCormick operator library in Julia
# See https://github.com/PSORLab/McCormick.jl
#############################################################################
# src/forward_operators/concave_increasing.jl
# Contains definitions of log, log2, log10, log1p, acosh, sqrt.
#############################################################################

for opMC in (:log, :log2, :log10, :log1p)
   opMC_kernel = Symbol(String(opMC)*"_kernel")
   dop = diffrule(:Base, opMC, :midcc)
   MCexp = quote
              xLc = z.lo
              xUc = z.hi
              xL = x.Intv.lo
              xU = x.Intv.hi
              midcc, cc_id = mid3(x.cc, x.cv, xU)
              midcv, cv_id = mid3(x.cc, x.cv, xL)
              dcv = (xUc > xLc) ? (xUc - xLc)/(xU - xL) : 0.0
              convex = dcv*(midcv - xL) + xLc
              concave = (NaNMath.$opMC)(midcc)
              concave_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*$dop
              convex_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
              convex, concave, convex_grad, concave_grad = cut(xLc, xUc, convex, concave, convex_grad, concave_grad)
              return MC{N, T}(convex, concave, z, convex_grad, concave_grad, x.cnst)
        end
    dop = diffrule(:Base, opMC, :(x.cv))
    dMCexp = quote
               xLc = z.lo
               xUc = z.hi
               xL = x.Intv.lo
               xU = x.Intv.hi
               midcc = mid3v(x.cv, x.cc, xU)
               midcv = mid3v(x.cv, x.cc, xL)
               deltaX = (xU - xL)
               slope = (xUc - xLc)/deltaX
               dcv = (xUc > xLc) ? slope : 0.0
               if deltaX == 0.0
                  convex = xLc
               else
                  convex = (xLc*(xU - midcv) + xUc*(midcv - xL))/deltaX
               end
               concave = (NaNMath.$opMC)(midcc)
               convex_grad = slope*x.cv_grad
               concave_grad = ($dop)*x.cc_grad
               return MC{N, Diff}(convex, concave, z, convex_grad, concave_grad, x.cnst)
         end

     @eval @inline ($opMC_kernel)(x::MC{N, T}, z::Interval{Float64}) where {N,T<:Union{NS,MV}} = $MCexp
     @eval @inline ($opMC_kernel)(x::MC{N, Diff}, z::Interval{Float64}) where {N} = $dMCexp
     @eval @inline ($opMC)(x::MC) = ($opMC_kernel)(x, ($opMC)(x.Intv))
end

@inline function acosh_kernel(x::MC{N, T}, z::Interval{Float64}) where {N, T<:Union{NS, MV}}
     isempty(x) && (return empty(x))
     (x.Intv.lo < 1.0 || x.Intv.hi < 1.0) && (return nan(MC{N,T}))
     xLc = z.lo
     xUc = z.hi
     xL = x.Intv.lo
     xU = x.Intv.hi
     midcc, cc_id = mid3(x.cc, x.cv, xU)
     midcv, cv_id = mid3(x.cc, x.cv, xL)
     dcv = (xUc > xLc) ? (xUc - xLc)/(xU - xL) : 0.0
     convex = dcv*(midcv - xL) + xLc
     concave = NaNMath.acosh(midcc)
     concave_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*inv(NaNMath.sqrt(x.cv - 1.0)*NaNMath.sqrt(x.cv + 1.0))
     convex_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
     convex, concave, convex_grad, concave_grad = cut(xLc, xUc, convex, concave, convex_grad, concave_grad)
     return MC{N, T}(convex, concave, z, convex_grad, concave_grad, x.cnst)
end
@inline function acosh_kernel(x::MC{N, Diff}, z::Interval{Float64}) where N
     isempty(x) && (return empty(x))
     (x.Intv.lo < 1.0 || x.Intv.hi < 1.0) && (return nan(MC{N,Diff}))
     xLc = z.lo
     xUc = z.hi
     xL = x.Intv.lo
     xU = x.Intv.hi
     midcc = mid3v(x.cv, x.cc, xU)
     midcv = mid3v(x.cv, x.cc, xL)
     deltaX = (xU - xL)
     slope = (xUc - xLc)/deltaX
     dcv = (xUc > xLc) ? slope : 0.0
     if deltaX == 0.0
        convex = xLc
     else
        convex = (xLc*(xU - midcv) + xUc*(midcv - xL))/deltaX
     end
     concave = NaNMath.acosh(midcc)
     convex_grad = slope*x.cv_grad
     concave_grad = inv(NaNMath.sqrt(x.cv - 1.0)*NaNMath.sqrt(x.cv + 1.0))*x.cc_grad
     return MC{N, Diff}(convex, concave, z, convex_grad, concave_grad, x.cnst)
end
@inline acosh(x::MC) = acosh_kernel(x, acosh(x.Intv))

@inline function sqrt_kernel(x::MC{N, T}, z::Interval{Float64}) where {N, T<:Union{NS, MV}}
     isempty(x) && (return empty(x))
     (x.Intv.lo < 0.0 || x.Intv.hi < 0.0) && (return nan(MC{N,T}))
     xLc = z.lo
     xUc = z.hi
     xL = x.Intv.lo
     xU = x.Intv.hi
     midcc, cc_id = mid3(x.cc, x.cv, xU)
     midcv, cv_id = mid3(x.cc, x.cv, xL)
     dcv = (xUc > xLc) ? (xUc - xLc)/(xU - xL) : 0.0
     convex = dcv*(midcv - xL) + xLc
     concave = NaNMath.sqrt(midcc)
     concave_grad = mid_grad(x.cc_grad, x.cv_grad, cc_id)*(0.5/NaNMath.sqrt(x.cv))
     convex_grad = mid_grad(x.cc_grad, x.cv_grad, cv_id)*dcv
     convex, concave, convex_grad, concave_grad = cut(xLc, xUc, convex, concave, convex_grad, concave_grad)
     return MC{N, T}(convex, concave, z, convex_grad, concave_grad, x.cnst)
end
@inline function sqrt_kernel(x::MC{N, Diff}, z::Interval{Float64}) where N
     isempty(x) && (return empty(x))
     (x.Intv.lo < 0.0 || x.Intv.hi < 0.0) && (return nan(MC{N,Diff}))
     xLc = z.lo
     xUc = z.hi
     xL = x.Intv.lo
     xU = x.Intv.hi
     midcc = mid3v(x.cv, x.cc, xU)
     midcv = mid3v(x.cv, x.cc, xL)
     deltaX = (xU - xL)
     slope = (xUc - xLc)/deltaX
     dcv = (xUc > xLc) ? slope : 0.0
     if deltaX == 0.0
        convex = xLc
     else
        convex = (xLc*(xU - midcv) + xUc*(midcv - xL))/deltaX
     end
     concave = NaNMath.sqrt(midcc)
     convex_grad = slope*x.cv_grad
     concave_grad = (0.5/NaNMath.sqrt(x.cv))*x.cc_grad
     return MC{N, Diff}(convex, concave, z, convex_grad, concave_grad, x.cnst)
end
@inline sqrt(x::MC) = sqrt_kernel(x, sqrt(x.Intv))

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# variations on the feedback policy

# Stephan's suggestion
type SolveSampledFeedback <: Policy
    problem::OPCSP
    solver::OPSolver
    rng::AbstractRNG
end
updater(p::SolveSampledFeedback) = OPCSPUpdater(p.problem)
type SampledFeedbackSolver <: Solver
    solver::OPSolver
    rng::AbstractRNG
end
POMDPs.solve(s::SampledFeedbackSolver, p::OPCSP) = SolveSampledFeedback(p, s.solver, s.rng)

function sampled_problem(op::OPCSP, b::OPCSPDistribution, rng::AbstractRNG)
    samp = rand!(rng, create_state(op), b)
    openset = cat(1, b.i, collect(b.open))
    samp_op = SimpleOP(op.r[openset]+samp.d[openset],
                       [],
                       b.remaining,
                       1,
                       findfirst(openset, op.stop),
                       op.distances[openset, openset])
    return SubsetOP(openset, samp_op)
end

function action(p::SolveSampledFeedback, b::Union{OPCSPDistribution,OPCSPBelief}, act::OPCSPAction=OPCSPAction(0))
    sp = sampled_problem(p.problem, b, p.rng)
    act.next = first_move(p.solver, sp)
    return act
end

# Try to have a bonus for the influence
type SolveInfluenceBonusFB <: Policy
    problem::OPCSP
    factor::Float64
end
updater(p::SolveInfluenceBonusFB) = OPCSPUpdater(p.problem)
type InfluenceBonusFBSolver <: Solver
    factor::Float64
end
POMDPs.solve(s::InfluenceBonusFBSolver, p::OPCSP) = SolveInfluenceBonusFB(p, s.factor)

function action(p::SolveInfluenceBonusFB, b::OPCSPDistribution, act::OPCSPAction=OPCSPAction(0))
    tail = tail_opcsp(p.problem, b)
    path = build_path(tail.op, solve_influence_bonus(tail.op, p.factor*b.remaining/p.problem.distance_limit))
    act.next = tail.included_nodes[path[2]]
    return act
end

function tail_opcsp(op::OPCSP, b::OPCSPDistribution)
    openset = cat(1, b.i, collect(b.open))
    tail = OPCSP(op.r[openset] + b.dist.mean[openset],
                 [],
                 b.dist.covariance[openset,openset],
                 b.remaining,
                 1,
                 findfirst(openset, op.stop),
                 op.distances[openset,openset])
    return SubsetOP(openset, tail)
end

function solve_influence_bonus(op::OPCSP, influence_factor::Float64)
    m = Model(solver=Gurobi.GurobiSolver(OutputFlag=false))
    N = length(op)

    # calculate influence
    influence = Array(Float64, length(op))
    for i in 1:length(op)
        if op.covariance[i,i] == 0.0
            influence[i] = 0.0
        else
            influence[i] = influence_factor*sum(abs(op.covariance[:,i]))/op.covariance[i,i]
        end
    end

    without_start = [1:op.start-1; op.start+1:N]
    without_stop = [1:op.stop-1; op.stop+1:N]
    without_both = intersect(without_start, without_stop)

    @defVar(m, x[1:N,1:N], Bin)
    @defVar(m, 2 <= u[without_start] <= N, Int)

    # @setObjective(m, Max, sum{ sum{op.r[i]*x[i,j], j=1:N}, i=1:N })
    @setObjective(m, Max, sum{ sum{op.r[i]*x[i,j], j=1:N}, i=1:N } + sum{ influence[j]*x[op.start,j], j=1:N })

    @addConstraint(m, sum{x[op.start,j], j=without_start} == 1)
    @addConstraint(m, sum{x[i,op.stop], i=without_stop} == 1)

    @addConstraint(m, connectivity[k=without_both], sum{x[i,k], i=1:N} == sum{x[k,j], j=1:N})

    @addConstraint(m, once[k=1:N], sum{x[k,j], j=1:N} <= 1)
    @addConstraint(m, sum{ sum{op.distances[i,j]*x[i,j], j=1:N}, i=1:N } <= op.distance_limit)
    @addConstraint(m, nosubtour[i=without_start,j=without_start], u[i]-u[j]+1 <= (N-1)*(1-x[i,j]))

    if op.start != op.stop
        @addConstraint(m, sum{x[op.stop,i],i=1:N} == 0)
    end

    status = JuMP.solve(m)

    if status != :Optimal
        warn("Not solved to optimality:\n$op")
    end

    soln = OPSolution(getValue(x), getValue(u))

    if distance(op,build_path(op,soln)) > op.distance_limit
        warn("Path Length: $(distance(op,build_path(op,soln))), Limit: $(op.distance_limit)")
    end
    return soln

end


text:
#Write a program that prints the next 20 leap years.

using Dates

function calc_leap()
     yr = year(now())
     i = 0
     while i < 20
        if (yr) % 4 == 0
           println(yr)
           i+=1
        end
        yr+=1
     end

end

calc_leap()

max_stars_repo_name: atott/Simple-Programming-Problems-Julia

max_issues_repo_name: atott/Simple-Programming-Problems-Julia

max_forks_repo_name: atott/Simple-Programming-Problems-Julia


text:
include("Conserved.jl")
include("RPCond.jl")

struct HLLC
end

function calcWaveSpeeds(left::RPCond,
	 				    right::RPCond)::Array{Float64,1}
	dl = left.density;
	pl = left.pressure;
	vl = left.velocity;
	cl = left.sound_speed
    dr = right.density;
	pr = right.pressure;
	vr = right.velocity;
	cr = right.sound_speed
	sl = min(vl - cl, vr - cr);
	sr = max(vl + cl, vr + cr);
	ss = (pr - pl + dl*vl*(sl - vl) - dr*vr*(sr - vr)) /
			(dl*(sl - vl) - dr*(sr - vr));
	return [sl, ss, sr];
end

function calcTotalEnergyDensity(state::RPCond)::Float64
	d = state.density
	p = state.pressure
	v = state.velocity
	e = state.energy
	return d*(e+v^2/2)
end

function calcStarredState(state::RPCond,
	 					  sk::Float64,
						  ss::Float64)::Conserved
	dk = state.density
	pk = state.pressure
	vk = state.velocity
	ds = dk*(sk - vk) / (sk - ss)
	ek = calcTotalEnergyDensity(state)
	mass = ds
	momentum = ds*ss
	energy = ek*ds / dk +
			ds*(ss - vk)*(ss + pk / dk / (sk - vk));
	return Conserved(mass, momentum, energy)
end

function calcMassDensity(state::RPCond)::Float64
	return state.density
end

function calcMomentumDensity(state::RPCond)::Float64
	return state.density*state.velocity
end

function primitive2conserved(state::RPCond)::Conserved
	return Conserved(calcMassDensity(state),
					 calcMomentumDensity(state),
					 calcTotalEnergyDensity(state))
end

function calcMassFlux(state::RPCond)::Float64
	return state.density*state.velocity
end

function calcMomentumFlux(state::RPCond)::Float64
	return state.density*state.velocity^2+state.pressure
end

function calcEnergyFlux(state::RPCond)::Float64
	energy_density = calcTotalEnergyDensity(state)
	enthalpy = energy_density+state.pressure
	return enthalpy*state.velocity
end

function primitive2flux(state::RPCond)::Conserved
	return Conserved(calcMassFlux(state),
					 calcMomentumFlux(state),
					 calcEnergyFlux(state))
end

function calcRestFrameFluxes(left::RPCond,
							 right::RPCond)::Conserved

	ws = calcWaveSpeeds(left, right)

	ul = primitive2conserved(left)
	ur = primitive2conserved(right)

	fl = primitive2flux(left)
	fr = primitive2flux(right)

	usl = calcStarredState(left, ws[1], ws[2])
	usr = calcStarredState(right, ws[3], ws[2])

	if ws[1]>0
		return fl
	elseif ws[2]>=0
		return fl + ws[1]*(usl-ul)
	elseif ws[3]>=0
		return fr + ws[3]*(usr-ur)
	else
		return fr
	end
end

function boostFlux(flux::Conserved,
				   velocity::Float64)::Conserved
	return Conserved(flux.mass,
					 flux.momentum+
					 velocity*flux.mass,
					 flux.energy+
					 0.5*flux.mass*velocity^2+
					 flux.momentum*velocity)
end

function boostCell(p::RPCond,
				   velocity::Float64)::RPCond
	return RPCond(p.density,
				  p.pressure,
				  p.velocity-velocity,
				  p.sound_speed,
				  p.energy)
end

function calcFlux(rs::HLLC,
				  left::RPCond,
				  right::RPCond,
				  velocity::Float64,
				  eos)::Conserved

	return boostFlux(calcRestFrameFluxes(boostCell(left,velocity),
										 boostCell(right,velocity)),
					 velocity)
end

max_stars_repo_name: bolverk/Sirion

max_issues_repo_name: bolverk/Sirion

max_forks_repo_name: bolverk/Sirion


text:
# installation management
precompile(__init_toolkit__, ())
precompile(libcuda, ())

# array
precompile(CuArray, (Vector{Int},))

# compilation
precompile(CUDACompilerTarget, (CuDevice,))
precompile(cufunction_compile, (CompilerJob,))
precompile(cufunction_link, (CompilerJob,NamedTuple{(:image, :entry, :external_gvars), Tuple{Vector{UInt8}, String, Vector{String}}}))
precompile(cufunction_cache, (CuContext,))
precompile(create_exceptions!, (CuModule,))
precompile(run_and_collect, (Cmd,))

# launch
precompile(cudaconvert, (Function,))
precompile(Core.kwfunc(cudacall), (NamedTuple{(:threads, :blocks), Tuple{Int64, Int64}},typeof(cudacall),CuFunction,Type{Tuple{}}))
precompile(Core.kwfunc(launch), (NamedTuple{(:threads, :blocks), Tuple{Int64, Int64}},typeof(launch),CuFunction))


text:
using RoME
using Test
# using RoMEPlotting

## list of tested factors
# - PriorPose2
# - PriorPoint2
# - Pose2Pose2
# - Pose2Point2BearingRange
# - Pose2Point2
# - DynPose2VelocityPrior
# - VelPose2VelPose2

@testset "Test PriorPose2 and Pose2Pose2" begin

fg = LightDFG( solverParams=SolverParams(algorithms=[:default, :parametric]))

# Add the first pose :x0
x0 = addVariable!(fg, :x0, Pose2)
prior = addFactor!(fg, [:x0], PriorPose2( MvNormal([10; 10; pi+1e-5], [0.1;0.1;0.05])))

for i in 0:3
    psym = Symbol("x$i")
    nsym = Symbol("x$(i+1)")
    addVariable!(fg, nsym, Pose2)
    pp = Pose2Pose2(MvNormal([10.0;0;pi/2], [0.1;0.1;0.1]))
    addFactor!(fg, [psym;nsym], pp)
end
initAll!(fg)

IIF.initParametricFrom!(fg)

vardict, result, varIds, Σ = IIF.solveGraphParametric(fg, useCalcFactor=true) #autodiff=:finite)

#TODO test +-pi used pi+1e-5
@test isapprox(vardict[:x0].val, [10, 10, -pi], atol = 1e-3)
@test isapprox(vardict[:x1].val, [0, 10, -pi/2], atol = 1e-3)
@test isapprox(vardict[:x2].val, [0, 0, 0], atol = 1e-3)
@test isapprox(vardict[:x3].val, [10, 0, pi/2], atol = 1e-3)
@test isapprox(vardict[:x4].val, [10, 10, -pi], atol = 1e-3)

# IIF.updateParametricSolution(fg, vardict)
# pl = plotSLAM2D(fg; lbls=true, solveKey=:parametric, point_size=4pt, drawPoints=false, drawContour=false)
end

##
# using Random
# Random.seed!(42);
# fg = LightDFG( solverParams=SolverParams(algorithms=[:default, :parametric]))

# fg.solverParams.graphinit = false
# pr_noise = [0.01, 0.01, 0.001]
# od_noise = [0.2;0.2;0.2]

# σ_bearing = 0.008
# σ_range = 0.01

# addVariable!(fg, :x0, Pose2)
# prpo = MvNormal([0,0,-pi], pr_noise)
# addFactor!(fg, [:x0], PriorPose2(MvNormal(rand(prpo), pr_noise)))

# addVariable!(fg, :l1, Point2, tags=[:LANDMARK])
# addVariable!(fg, :l2, Point2, tags=[:LANDMARK])
# # p2br = Pose2Point2BearingRange(Normal(pi/4,0.01),Normal(sqrt(2),0.1))
# p2br = Pose2Point2BearingRange(Normal(pi/4 + rand(Normal(0,σ_bearing)), σ_bearing),
#                                Normal(sqrt(2) + rand(Normal(0,σ_range)), σ_range))
# addFactor!(fg, [:x0; :l1], p2br)

# addVariable!(fg, :x1, Pose2)
# pp = MvNormal([1.0,0,0], od_noise)
# addFactor!(fg, [:x0,:x1], Pose2Pose2(MvNormal(rand(pp), od_noise)))

# p2br = Pose2Point2BearingRange(Normal(pi/2 + rand(Normal(0,σ_bearing)), σ_bearing),
#                                Normal(1 + rand(Normal(0,σ_range)), σ_range))
# addFactor!(fg, [:x1; :l1], p2br)

# p2br = Pose2Point2BearingRange(Normal(-pi/2 + rand(Normal(0,σ_bearing)), σ_bearing),
#                                Normal(1 + rand(Normal(0,σ_range)), σ_range))
# addFactor!(fg, [:x1; :l2], p2br)

# addVariable!(fg, :x2, Pose2)
# pp = MvNormal([1.0,0,0], od_noise)
# addFactor!(fg, [:x1,:x2], Pose2Pose2(MvNormal(rand(pp), od_noise)))

# p2br = Pose2Point2BearingRange(Normal(3pi/4 + rand(Normal(0,σ_bearing)), σ_bearing),
#                                Normal(sqrt(2) + rand(Normal(0,σ_range)), σ_range))
# addFactor!(fg, [:x2; :l1], p2br)

# p2br = Pose2Point2BearingRange(Normal(-3pi/4 + rand(Normal(0,σ_bearing)), σ_bearing),
#                                Normal(sqrt(2) + rand(Normal(0,σ_range)), σ_range))
# addFactor!(fg, [:x2; :l2], p2br)

# initAll!(fg)

# IIF.initParametricFrom!(fg)

# vardict, result, varIds, Σ = IIF.solveGraphParametric(fg, useCalcFactor=true) #autodiff=:finite)
# # vardict, result, varIds, Σ = IIF.solveGraphParametric(fg, useCalcFactor=true, autodiff=:finite)

# IIF.updateParametricSolution(fg, vardict)

# pl = plotSLAM2D(fg; lbls=true, solveKey=:parametric, point_size=4pt, drawPoints=false, drawContour=false)

##

@testset "Test Parametric PriorPose2 and Pose2Point2" begin

fg = LightDFG( solverParams=SolverParams(algorithms=[:default, :parametric]))

addVariable!(fg, :x1, Pose2)
addVariable!(fg, :l1, Point2)

addFactor!(fg, [:x1], PriorPose2(MvNormal([0.,0, 0], [0.01, 0.01, 0.01])))

addFactor!(fg, [:x1; :l1], Pose2Point2(MvNormal([1.0, 1], [0.1,0.1])))

initAll!(fg)

IIF.initParametricFrom!(fg)

vardict, result, varIds, Σ = IIF.solveGraphParametric(fg, useCalcFactor=true) #autodiff=:finite)

@test isapprox(vardict[:x1].val, [0, 0, 0], atol = 1e-3)
@test isapprox(vardict[:l1].val, [1,  1], atol = 1e-3)

# IIF.updateParametricSolution(fg, vardict)
# pl = plotSLAM2D(fg; lbls=true, solveKey=:parametric, point_size=4pt, drawPoints=false, drawContour=false)
# push!(pl, Coord.cartesian(fixed=true),style(background_color=RGB(1,1,1)))
end

##
@testset "Test Parametric PriorPoint2 and Pose2Point2BearingRange" begin
fg = LightDFG( solverParams=SolverParams(algorithms=[:default, :parametric]))

addVariable!(fg, :x1, Pose2)
addVariable!(fg, :l1, Point2)
addVariable!(fg, :l2, Point2)

addFactor!(fg, [:l1], PriorPoint2(MvNormal([1., 1], [0.01, 0.01])))
addFactor!(fg, [:l2], PriorPoint2(MvNormal([1.,-1], [0.01, 0.01])))

addFactor!(fg, [:x1; :l1], Pose2Point2BearingRange(Normal(pi/4, 0.01), Normal(sqrt(2), 0.1)))
addFactor!(fg, [:x1; :l2], Pose2Point2BearingRange(Normal(3pi/4, 0.01), Normal(sqrt(2), 0.1)))

initAll!(fg)

IIF.initParametricFrom!(fg)

vardict, result, varIds, Σ = IIF.solveGraphParametric(fg, useCalcFactor=true) #autodiff=:finite)

@test isapprox(vardict[:x1].val, [2, 0, pi/2], atol = 1e-3)
@test isapprox(vardict[:l1].val, [1,  1], atol = 1e-3)
@test isapprox(vardict[:l2].val, [1, -1], atol = 1e-3)

# IIF.updateParametricSolution(fg, vardict)
# pl = plotSLAM2D(fg; lbls=true, solveKey=:parametric, point_size=4pt, drawPoints=false, drawContour=false)
end

@testset "Test Parametric DynPose2VelocityPrior and VelPose2VelPose2" begin

fg = LightDFG( solverParams=SolverParams(algorithms=[:default, :parametric]))

# add first pose locations
addVariable!(fg, :x0, DynPose2; nanosecondtime=0)

# Prior factor as boundary condition
pp0 = DynPose2VelocityPrior(MvNormal(zeros(3), [0.01; 0.01; 0.001]), MvNormal([10.0;0], [0.1; 0.1]))
addFactor!(fg, [:x0;], pp0)

addVariable!(fg, :x1, DynPose2;  nanosecondtime=1000_000_000)

# conditional likelihood between Dynamic Point2
dp2dp2 = VelPose2VelPose2(MvNormal([10.0;0;0], [0.01;0.01;0.001]), MvNormal([0.0;0], [0.1; 0.1]))
addFactor!(fg, [:x0;:x1], dp2dp2)

initAll!(fg)

vardict, result, varIds, Σ = IIF.solveGraphParametric!(fg)

@test isapprox(vardict[:x0].val, [0, 0, 0, 10, 0], atol = 1e-3)
@test isapprox(vardict[:x1].val, [10, 0, 0, 10, 0], atol = 1e-3)

end

@testset "Test Parametric PriorPoint2 and Point2Point2Range" begin

fg = LightDFG( solverParams=SolverParams(algorithms=[:default, :parametric]))

addVariable!(fg, :x1, Point2)
addVariable!(fg, :l1, Point2)
addVariable!(fg, :l2, Point2)
addVariable!(fg, :l3, Point2)

addFactor!(fg, [:l1], PriorPoint2(MvNormal([0., 0], [0.01, 0.01])))
addFactor!(fg, [:l2], PriorPoint2(MvNormal([1., 0], [0.01, 0.01])))
addFactor!(fg, [:l3], PriorPoint2(MvNormal([0., 1], [0.01, 0.01])))

addFactor!(fg, [:x1; :l1], Point2Point2Range(Normal(sqrt(2), 0.1)))
addFactor!(fg, [:x1; :l2], Point2Point2Range(Normal(1.0, 0.1)))
addFactor!(fg, [:x1; :l3], Point2Point2Range(Normal(1.0, 0.1)))

# initAll!(fg)
#FIXME needs initializaiton from non-parametric to converge
solveTree!(fg)
IIF.initParametricFrom!(fg)

vardict, result, varIds, Σ = IIF.solveGraphParametric(fg)

@test isapprox(vardict[:x1].val, [1, 1], atol = 1e-3)

end


text:
using Distributed

@everywhere include("parallel-binomial.jl")
using PyPlot

n = 100
p = 0.5
reps = 1000
rounds = 1000

println("Starting computation")

@time c = parallel_binomial_counts(n, p, reps, rounds)

a = minimum(keys(c))
b = maximum(keys(c))
s = sum(c)

x = collect(a:b)
y = [c[t] / s for t in x]

bar(x, y)
grid()
title("Binomial distribution B($n, $p)")
savefig("binomial.pdf")


text:
module Ising2D

export NetworkDynamic, IsingModel

using Statistics, ProgressMeter

# function FixBounds(dim::Integer, indices::Vector{Vector{Int}})
#     for i ∈ 1:4
#         for j ∈ 1:2
#             if indices[i][j] == dim + 1
#                 indices[i][j] = 1
#             elseif indices[i][j] == 0
#                 indices[i][j] = dim
#             end
#         end
#     end
#     return indices
# end

function FindEnergy(Network::Matrix{Int8}, dim::Integer)
    E = 0.0
    for i ∈ 1:dim
        for j ∈ 1:dim
            Neighbors = [[1 + i % dim, j], [i, 1 + j % dim], [(-2 + i + dim) % dim + 1, j], [i, (-2 + j + dim) % dim + 1]]
            E += -sum(Network[i, j] * Network[CartesianIndex.(Tuple.(Neighbors))]) / 2
        end
    end
    return E / dim ^2
end

function FindMagnetization(Network::Matrix{Int8})
    return mean(Network)
end

function InitialCondition(dim::Integer)
    return rand([Int8(-1), Int8(1)], dim, dim)
end

function NetworkDynamic(β::Real; dim::Integer = 150, MCLSize::Integer = 200, NLSize::Integer = 5000, ProgBar::Bool = true)
    if ProgBar == true
        Prog = Progress(MCLSize * NLSize)
    end
    Networks = Array{Matrix{Int8},1}(undef, MCLSize)
    Network = InitialCondition(dim)
    for MCL ∈ 1:MCLSize
        for NL ∈ 1:NLSize
            if ProgBar == true
                next!(Prog)
            end
            i, j = rand(1:dim, 2)
            Neighbors = [
                [1 + i % dim, j],
                [i, 1 + j % dim],
                [(-2 + i + dim) % dim + 1, j],
                [i, (-2 + j + dim) % dim + 1]
            ]
            ΔE = 2 * Network[i, j] * sum(Network[CartesianIndex.(Tuple.(Neighbors))])
            if ΔE <= 0
                Network[i, j] *= -1
            elseif rand() < exp(-ΔE * β)
                Network[i, j] *= -1
            end
        end
        Networks[MCL] = copy(Network)
    end
    return Networks
end

function MonteCarloStep(StepSize::Integer, Network::Matrix{Int8}, dim::Integer, β::Real)
    ΔEₙₑₜ = 0
    ΔMₙₑₜ = 0
    for NL ∈ 1:StepSize
        i, j = rand(1:dim, 2)
        Neighbors = [
            [1 + i % dim, j],
            [i, 1 + j % dim],
            [(-2 + i + dim) % dim + 1, j],
            [i, (-2 + j + dim) % dim + 1]
        ]
        ΔE = 2 * Network[i, j] * sum(Network[CartesianIndex.(Tuple.(Neighbors))])
        if ΔE <= 0 || rand() < exp(-ΔE * β)
            Network[i, j] *= -1
            ΔEₙₑₜ += ΔE / dim^2
            ΔMₙₑₜ += 2 * Network[i, j] / dim^2
        end
    end
    return ΔEₙₑₜ, ΔMₙₑₜ
end

function IsingModel(β::Real; dim::Integer = 15, MCLSize::Integer = 1000, NLSize::Integer = 1000, SkipNum::Integer = 200, ProgBar::Bool = true)
    if ProgBar == true
        Prog = Progress(MCLSize)
    end
    Network = InitialCondition(dim)
    E = FindEnergy(Network, dim)
    M = FindMagnetization(Network)
    EₙₑₜList = Array{Float64,1}(undef, MCLSize)
    MₙₑₜList = Array{Float64,1}(undef, MCLSize)
    for MCL ∈ 1:MCLSize
        if ProgBar == true
            next!(Prog)
        end
        ΔEₙₑₜ, ΔMₙₑₜ = MonteCarloStep(NLSize, Network, dim, β)
        E += ΔEₙₑₜ
        M += ΔMₙₑₜ
        EₙₑₜList[MCL] = E
        MₙₑₜList[MCL] = abs(M)
    end
    Ē, M̄ = mean(EₙₑₜList[SkipNum:end]), mean(MₙₑₜList[SkipNum:end])
    Cᵥ, χ = β^2 * var(EₙₑₜList[SkipNum:end]) * dim^2, β * var(MₙₑₜList[SkipNum:end])
    return Ē, M̄, Cᵥ, χ
end

end


text:
function stemleaf(; y::AbstractArray{Int} = Int[])

    #Stem and Leaf calculations
    sort!(y)
    leaf = Int[x % 10 for x in y]
    stem = Int[fld(x, 10) for x in y]

    position = Int[]
    i = stem[1]
    counter = 1

    for elem in stem
        if elem == i
            push!(position, counter)
            counter += 1
        else
            counter = 1
            push!(position, counter)
            i = elem
            counter += 1
        end
    end
    #End Stem and Leaf calculations

    v = scatterplot(x = position, y = stem, group = leaf)
    ylim!(v, min = minimum(stem), max = maximum(stem), reverse = true)
    legend!(v, show = false)
    colorscheme!(v, palette = "black")

    #Switch from circles to text
    v.marks[1]._type = "text"
    v.marks[1].properties.enter.shape = VegaValueRef(value = "")
    v.marks[1].properties.enter.size = VegaValueRef(value = "")
    v.marks[1].properties.enter.text = VegaValueRef(field = "group")
    v.marks[1].properties.enter.fillOpacity = VegaValueRef(value = 1)

    #Text offsets
    v.marks[1].properties.enter.baseline = VegaValueRef(value = "middle")
    v.marks[1].properties.enter.y = VegaValueRef(scale = "y", field = "y", offset = 0.5)
    v.marks[1].properties.enter.x = VegaValueRef(scale = "x", field = "x", offset = -10)

    #Hide x ticks, axis, labels
    xlab!(v, title = " ", properties = Dict{Any, Any}("ticks" => Dict{Any, Any}("stroke" => VegaValueRef(value = "transparent")),
                                         "axis" => Dict{Any, Any}("stroke" => VegaValueRef(value = "transparent")),
                                         "labels" => Dict{Any, Any}("fill" => VegaValueRef(value = "transparent"))
                                        )
         )

    #Hide y ticks
    ylab!(v, title = " ", ticks = length(unique(stem)), tickSizeEnd = 0, properties = Dict{Any, Any}("ticks" => Dict{Any, Any}("stroke" => VegaValueRef(value = "transparent"))
                                        )
         )

    return v

end

max_stars_repo_name: jpthor/vega


text:
function extract_media(src)
    r = ZipFile.Reader(src)
    media_files = filter(f -> occursin(r"^word/media", f.name), r.files)
    Dict(
        replace(f.name, "word/" => "") => read(f)
        for f
        in r.files
        if occursin(r"^word/media", f.name)
    )
end


text:
using .StatsBase

@inline @inbounds function skewedGaussianQuadratic(x::AbstractArray, p::Function, t::Real)
    (μ, η, κ) = p(t)
    x = μ*x[1] + μ*x[1]*exp(-κ*x[1])*(κ*x[1]-2.0)
    return SVector{1}(x)
end

@inline @inbounds function skewedGaussianQuadratic_σ(σ::AbstractArray, p::Function, t::Real)
    (μ, η, κ) = p(t)
    σ = η
end

"""Skewed Pitchfork Bifurcation"""
function skewedGaussianQuadratic(P::Process)
    seed(P.solver_rng)
    prob = SDEProblem(P.process, skewedGaussianQuadratic_σ, P.X0, (P.transient_t0, P.tmax), tuplef2ftuple(P.parameter_profile, P.parameter_profile_parameters))
    sol = dsolve(prob, P.alg; dt = P.dt, saveat=P.savedt, P.solver_opts...)
end

skewedGaussianQuadraticSim = Process(
    process = skewedGaussianQuadratic,
    X0 = [0.0],
    parameter_profile = (constantParameter, constantParameter, rampInterval),
    parameter_profile_parameters = ((-1.0,), (2.0,), (0.01, 0.05, 0.0, 1000.0)),
    transient_t0 = -100.0,
    t0 = 0.0,
    dt = 0.001,
    savedt = 0.001,
    tmax = 1000.0,
    alg = EM())
export skewedGaussianQuadraticSim

"""A system combining multiple classes of features. Combine Chen's system with gaussian noise and an exponential measurement function"""
function seeIfItSticks(P::Process)
    seed(P.solver_rng)
    prob = ODEProblem(chensSystem, P.X0, (P.transient_t0, P.tmax), tuplef2ftuple(P.parameter_profile[1:3], P.parameter_profile_parameters[1:3]), jac=chensSystem_J)
    sol = dsolve(prob, P.alg; dt = P.dt, saveat=P.savedt, P.solver_opts...)
    η = parameterseries(P, p=[4], transient=true)
    α = parameterseries(P, p=[5], transient=true)
    sol = sol .+ η.*randn(size(sol))
    sol = (sol.*exp.(α.*(sol.-mean(sol, dims=2))))'
end

seeIfItSticksSim = Process(
    process = seeIfItSticks,
    X0 = [-3.0, 2.0, 20.0],
    parameter_profile = (ramp, ramp, ramp, constant, constant),
    parameter_profile_parameters = ((42.0, 42.0, 0.0, 500.0), (11.0, 11.0, 0.0, 500.0), (28.0, 28.0, 0.0, 500.0), 2.0, 0.01), # (𝑎, 𝑏, 𝑐, η, α)
    transient_t0 = -100.0,
    t0 = 0.0,
    dt = 0.001,
    savedt = 0.025,
    tmax = 500.0,
    alg = RK4(),
    solver_opts = Dict(:adaptive => false))
export seeIfItSticksSim


text:
"""
Implementation of Approximate Bayesian Computation (ABC) Algorithms

Includes standard ABC algorithm and Population Monte Carlo (PMC)variation.
Functions to make model definition easier are included as well as a variety
of methods for computing non parametric distances between data.
"""
module ApproximateBayesianComputation

import StatsBase, Distributions, JLD
import Base.copy

export ParticleDimension,
    Univariate,
    Multivariate,
    DistanceMeasure,
    SingleMeasure,
    MultiMeasure,
    AbcPmcStep,
    UnivariateAbcPmc,
    MultivariateAbcPmc,
    SingleMeasureAbcPmc,
    MultiMeasureAbcPmc,
    USAbcPmc,
    UMAbcPmc,
    MSAbcPmc,
    MMAbcPmc,
    AbcPmc
export kernel_weights, pmc_start, pmc_step, abc_pmc, abc_pmc_warmstart
export abc_standard, ABCResult, abc_standard1D, abc_standardMultiD
export boxcar_accept, gaussian_accept
export maxdist_bandwidth, gaussianprob_bandwidth
export weightedstd, weightedstd2
export make_normal_kernel, make_jointnormal_kernel, make_truncatednormal_kernel,
    make_joint_kernel
export make_sample_prior, make_density_prior, make_model_prior
export sample_particle_distance, ppmc_start, ppmc_step, pabc_pmc
export pabc_standard, find_particle, pabc_pmc_warmstart
export totalsamples_thresholds
export perm_distance, sumrank_distance, maxrank_distance, normrank_distance, minarea_distance
export quantile_acceptbw, independentquantile_acceptbw, jointquantile_acceptbw
export duration_to_string

include("AbcPmcStep.jl")
include("abc_pmc.jl")
include("abc_standard.jl")
include("accept_reject.jl")
include("initial_acceptbw.jl")
include("kernel_bandwidth.jl")
include("make_kernel_model.jl")
include("make_prior_model.jl")
include("parallel_abc_pmc.jl")
include("parallel_abc_standard.jl")
include("post_processing.jl")
include("rank_distances.jl")
include("shrink_acceptbw.jl")
include("utils.jl")

end # module


text:
include("drawBlocks.jl")

# draws the diamond and the block of true code together
function drawCondition(env::cairo_env,condition::String,x::Float64,y::Float64,true_code::String)

    # draw the diamond
    dia_size = drawTextDiamond(env,condition,x,y)

    # draw the arrow going down from the diamond
    y+=dia_size

    # later return where the no arrow is
    false_x = x
    false_y =y +line_dist_y

    # draw a line from the bottom of the diamond down
    drawArrowLineCairo(env,x,y,0.0,line_dist_y, arrow_width, arrow_height)

    # move coordinates to the right point of the diamond
    x += dia_size/2
    y -= dia_size/2

    # make an arrow from the diamond to a block of true code
    drawArrowLineCairo(env,x,y, line_dist_x,0.0, arrow_width,arrow_height)

    # set the x,y to the end of the arrow
    x += line_dist_x

    # draw the box at the end of the arrow
    rect_w,rect_h = drawTextBox(env,true_code,x,y,false)

    # return the x,y of the end of the true code block
    x = x + rect_w/2 - line_dist_x
    y = y + rect_h

    # we need to return all these since the program needs to be able to
    # figure out how to backtrack to the "No" arrow
    return x,y,false_x,false_y
end

function connectWArrowLines(env,origin_x,origin_y,end_x,end_y)

    # get the distance to draw
    x_dist = end_x - origin_x
    y_dist = end_y - origin_y

    # draw the first relative line before the joining line
    drawRelLineCairo(env,origin_x,origin_y,y_dist,0.0)

    drawArrowLineCairo(env,origin_x,origin_y+y_dist,0.0,x_dist)

end

# draws the control flow in the array
#input is an array of control blocks
#
# control flow
#  if the control block is an if, draw from the true_code block
#
#  An else if block draws from the 'No' arrow
#
# an else block also draws from the 'No' arrow
# an else will rejoin the true_code arrow after the else code block
function drawFlow(input,filename::String, line_width::Float64)

    current_x = 350.0
    current_y = 10.0
    dia_size = 0.0
    rect_w = 0.0
    rect_h = 0.0

    # holds coordinates for the true and false paths
    true_x = 350.0
    true_y = 10.0
    false_x =350.0
    false_y = 10.0

    rect_x = 0.0
    rect_y = 0.0

    # index and indicator of program scope
    scope_stack = Vector{ Vector{Tuple{Int64,Int64}} }(undef, 0)
    scope_index::Int64 = 0
    block_index::Int64 = 0

    num_blocks = length(input)
    size_mult = 400.0

    # TODO make the size of the svg proportional to the number of elements
    # in input
    env = initCairo(filename,size_mult*num_blocks,size_mult*num_blocks)

    setLineWidthCairo(env ,line_width)

    for i in 1:length(input)

        if typeof(input[i]) == if_block

            # draw the block and store where the true_code block ends
            current_x,current_y,false_x,false_y = drawCondition(env,input[i].condition,current_x,current_y,input[i].true_code)

            # increase the scope or the loop if the else needs to be carried on
            if input[i].has_else == true

                push!(scope_stack,(current_x,current_y))
                scope_index+=1
                block_index+=1
            else # connect the previous lines

            end

        elseif typeof(input[i]) == if_else_block

            # if the previous element is not an if, then an error should be
            # thrown
            # I am assuming that the previous element should be an if
            if i > 1 && typeof(input[i-1]) == if_block
                # backtrack to the 'No' section, and draw another condition from there
                current_x,current_y,false_x,false_y = drawCondition(env,input[i].condition,false_x,false_y,input[i].true_code)

            # TODO output error and stop
            else

            end

            # if the else continues, do not connect the line from the true code
            if input[i].has_else == true
                push!(scope_stack[scope_index],(current_x,current_y))
                block_index +=1

            else # connect the previous lines

            end

        elseif typeof(input[i]) == else_block

            # if the previous element is not an if, then an error should be
            # thrown
            # I am assuming that the previous element should be an if
            if i > 1 && (typeof(input[i-1]) == if_block || typeof(input[i-1]) == if_else_block)

                # just draw the true_code at for the else statement
                rect_w,rect_h = drawTextBox(env,input[i].true_code,false_x,false_y,true)

                false_y+=rect_h

            # TODO output error and stop
            else

            end

            # draw the arrow if it is not at the end
            if i != num_blocks
                drawArrowLineCairo(env,current_x,current_y,0.0,line_dist_y,arrow_width,arrow_height)
            end

            current_y+=line_dist_y

            # connect the blocks with arrow lines

        # this is for regular blocks of code, with no conditions
        elseif typeof(input[i]) == String

            # just draw the true_code at for the else statement
            rect_w,rect_h = drawTextBox(env,input[i],false_x,false_y,true)

            false_y+=rect_h

            current_y = false_y

            # draw the arrow if it is not at the end
            if i != num_blocks
                drawArrowLineCairo(env,current_x,current_y,0.0,line_dist_y,arrow_width,arrow_height)
            end

            current_y+=line_dist_y

        else
            # type in the array was not provisioned for
            throw(ArgumentError("Element ",i," type not recognized."))
        end

    end

    destroyCairo(env)
end

max_stars_repo_name: aklapatch/CodeVisualized

max_issues_repo_name: aklapatch/CodeVisualized

max_forks_repo_name: aklapatch/CodeVisualized


text:
"""
A basic cost function, where the computed cost is the size
(number of children) of the current expression.
"""
function astsize(n::ENode, g::EGraph, an::Type{<:AbstractAnalysis})
    cost = 1 + arity(n)
    for id ∈ n.args
        eclass = geteclass(g, id)
        !hasdata(eclass, an) && (cost += Inf; break)
        cost += last(getdata(eclass, an))
    end
    return cost
end

"""
A basic cost function, where the computed cost is the size
(number of children) of the current expression, times -1.
Strives to get the largest expression
"""
astsize_inv(n::ENode, g::EGraph, an::Type{<:AbstractAnalysis}) = -1 * astsize(n, g, an)

"""
An [`AbstractAnalysis`](@ref) that computes the cost of expression nodes
and chooses the node with the smallest cost for each E-Class.
This abstract type is parametrised by a function F.
This is useful for the analysis storage in [`EClass`](@ref)
"""
abstract type ExtractionAnalysis{F} <: AbstractAnalysis end

make(a::Type{ExtractionAnalysis{F}}, g::EGraph, n::ENode) where F = (n, F(n, g, a))

join(a::Type{<:ExtractionAnalysis}, from, to) = last(from) <= last(to) ? from : to

islazy(a::Type{<:ExtractionAnalysis}) = true

function rec_extract(g::EGraph, an::Type{<:ExtractionAnalysis}, id::EClassId)
    eclass = geteclass(g, id)
    anval = getdata(eclass, an, missing)
    if anval === missing
        analyze!(g, an, id)
        anval = getdata(eclass, an)
    end
    (cn, ck) = anval
    (!isterm(termtype(cn)) || ck == Inf) && return cn.head

    extractnode(g, cn, an; eclass=eclass)
end

function extractnode(g::EGraph, n::ENode, an::Type{<:ExtractionAnalysis}; eclass=nothing)
    children = map(n.args) do a
        rec_extract(g, an, a)
    end

    meta = nothing
    if !isnothing(eclass)
        meta = getdata(eclass, MetadataAnalysis, nothing)
    end
    T = termtype(n)
    similarterm(T, n.head, children; metadata = meta)
end

# TODO CUSTOMTYPES document how to for custom types
# TODO maybe extractor can just be the array of extracted children?
function extractnode(g::EGraph, n::ENode{Expr}, extractor::Function)::Expr
    return Expr(n.head, map(extractor, n.args)...)
end

function extractnode(g::EGraph, n::ENode{T}, extractor::Function) where T
    if arity(n) > 0
        error("ENode extraction is not defined for non-literal type $T")
    end
    return n.head
end

"""
Given an [`ExtractionAnalysis`](@ref), extract the expression
with the smallest computed cost from an [`EGraph`](@ref)
"""
function extract!(g::EGraph, a::Type{ExtractionAnalysis{F}} where F; root=-1)
    # @show root g.root
    if root == -1
        root = g.root
    end
    # @show root g.root
    analyze!(g, a, root)
    !(a ∈ g.analyses) && error("Extraction analysis is not associated to EGraph")
    rec_extract(g, a, root)
end

"""
Given a cost function, extract the expression
with the smallest computed cost from an [`EGraph`](@ref)
"""
function extract!(g::EGraph, costfun::Function; root=-1)
    extran = ExtractionAnalysis{costfun}
    extract!(g, extran; root=root)
end

macro extract(expr, theory, costfun)
    quote
        let g = EGraph($expr)
            saturate!(g, $theory)
            ex = extract!(g, $costfun)
            (g, ex)
        end
    end |> esc
end


text:
"""
    VectorSpace

Abstract type for all vector spaces.
"""
abstract type VectorSpace end

Base.:(==)(::VectorSpace, ::VectorSpace) = false
Base.issubset(::VectorSpace, ::VectorSpace) = false
Base.intersect(s₁::VectorSpace, s₂::VectorSpace) = throw(MethodError(intersect, (s₁, s₂)))
Base.union(s₁::VectorSpace, s₂::VectorSpace) = throw(MethodError(union, (s₁, s₂)))

# Parameter space

"""
    ParameterSpace <: VectorSpace

Space of a parameter corresponding to a commutative field.
"""
struct ParameterSpace <: VectorSpace end

Base.:(==)(::ParameterSpace, ::ParameterSpace) = true
Base.issubset(::ParameterSpace, ::ParameterSpace) = true
Base.intersect(::ParameterSpace, ::ParameterSpace) = ParameterSpace()
Base.union(::ParameterSpace, ::ParameterSpace) = ParameterSpace()

dimension(::ParameterSpace) = 1
_firstindex(::ParameterSpace) = 1
_lastindex(::ParameterSpace) = 1
indices(::ParameterSpace) = Base.OneTo(1)

_findposition(i, ::ParameterSpace) = i

# Sequence spaces

"""
    SequenceSpace <: VectorSpace

Abstract type for all sequence spaces.
"""
abstract type SequenceSpace <: VectorSpace end

"""
    BaseSpace <: SequenceSpace

Abstract type for all sequence spaces that are not a [`TensorSpace`](@ref) but can be interlaced to form one.
"""
abstract type BaseSpace <: SequenceSpace end

"""
    TensorSpace{T<:NTuple{N,BaseSpace} where {N}} <: SequenceSpace

[`SequenceSpace`](@ref) resulting from the tensor product of some [`BaseSpace`](@ref).

Fields:
- `spaces :: T`
"""
struct TensorSpace{T<:NTuple{N,BaseSpace} where {N}} <: SequenceSpace
    spaces :: T
end

spaces(s::TensorSpace) = s.spaces

⊗(s₁::BaseSpace, s₂::BaseSpace) = TensorSpace((s₁, s₂))
⊗(s₁::TensorSpace, s₂::TensorSpace) = TensorSpace((s₁.spaces..., s₂.spaces...))
⊗(s₁::TensorSpace, s₂::BaseSpace) = TensorSpace((s₁.spaces..., s₂))
⊗(s₁::BaseSpace, s₂::TensorSpace) = TensorSpace((s₁, s₂.spaces...))

Base.@propagate_inbounds Base.getindex(s::TensorSpace, i::Int) = getindex(s.spaces, i)
Base.@propagate_inbounds Base.getindex(s::TensorSpace, u::AbstractRange{Int}) = TensorSpace(getindex(s.spaces, u))
Base.@propagate_inbounds Base.getindex(s::TensorSpace, u::AbstractVector{Int}) = TensorSpace(getindex(s.spaces, u))
Base.@propagate_inbounds Base.getindex(s::TensorSpace, c::Colon) = TensorSpace(getindex(s.spaces, c))

Base.front(s::TensorSpace) = TensorSpace(Base.front(s.spaces))
Base.tail(s::TensorSpace) = TensorSpace(Base.tail(s.spaces))

#

function Base.:(==)(s₁::TensorSpace{<:NTuple{N,BaseSpace}}, s₂::TensorSpace{<:NTuple{N,BaseSpace}}) where {N}
    s₁[1] == s₂[1] && return Base.tail(s₁) == Base.tail(s₂)
    return false
end
Base.:(==)(s₁::TensorSpace{<:Tuple{BaseSpace}}, s₂::TensorSpace{<:Tuple{BaseSpace}}) =
    s₁[1] == s₂[1]
function Base.issubset(s₁::TensorSpace{<:NTuple{N,BaseSpace}}, s₂::TensorSpace{<:NTuple{N,BaseSpace}}) where {N}
    issubset(s₁[1], s₂[1]) && return issubset(Base.tail(s₁), Base.tail(s₂))
    return false
end
Base.issubset(s₁::TensorSpace{<:Tuple{BaseSpace}}, s₂::TensorSpace{<:Tuple{BaseSpace}}) =
    issubset(s₁[1], s₂[1])
Base.intersect(s₁::TensorSpace{<:NTuple{N,BaseSpace}}, s₂::TensorSpace{<:NTuple{N,BaseSpace}}) where {N} =
    TensorSpace(map(intersect, s₁.spaces, s₂.spaces))
Base.union(s₁::TensorSpace{<:NTuple{N,BaseSpace}}, s₂::TensorSpace{<:NTuple{N,BaseSpace}}) where {N} =
    TensorSpace(map(union, s₁.spaces, s₂.spaces))

dimension(s::TensorSpace) = mapreduce(dimension, *, s.spaces)
dimension(s::TensorSpace, i::Int) = dimension(s.spaces[i])
dimensions(s::TensorSpace) = map(dimension, s.spaces)
_firstindex(s::TensorSpace) = map(_firstindex, s.spaces)
_lastindex(s::TensorSpace) = map(_lastindex, s.spaces)
indices(s::TensorSpace) = Base.Iterators.ProductIterator(map(indices, s.spaces))

_findindex_constant(::TensorSpace{<:NTuple{N,BaseSpace}}) where {N} = ntuple(i -> 0, Val(N))

@generated function _findposition(α::NTuple{N,Int}, s::TensorSpace{<:NTuple{N,BaseSpace}}) where {N}
    # follows column major convention
    idx = :(_findposition(α[1], s.spaces[1]))
    n = 1
    for i ∈ 2:N
        n = :(dimension(s.spaces[$i-1]) * $n)
        idx = :($n * (_findposition(α[$i], s.spaces[$i]) - 1) + $idx)
    end
    return idx
end
_findposition(u::NTuple{N,Any}, s::TensorSpace{<:NTuple{N,BaseSpace}}) where {N} =
    _findposition(Base.Iterators.ProductIterator(map(_colon2indices, u, s.spaces)), s)
_colon2indices(u, s) = u
_colon2indices(c::Colon, s) = indices(s)
function _findposition(u::Base.Iterators.ProductIterator, s::TensorSpace)
    v = Vector{Int}(undef, length(u))
    @inbounds for (i, uᵢ) ∈ enumerate(u)
        v[i] = _findposition(uᵢ, s)
    end
    return v
end
_findposition(u::AbstractVector{NTuple{N,Int}}, s::TensorSpace{<:NTuple{N,BaseSpace}}) where {N} =
    map(α -> _findposition(α, s), u)
_findposition(::NTuple{N,Colon}, ::TensorSpace{<:NTuple{N,BaseSpace}}) where {N} = Colon()
_findposition(c::Colon, ::TensorSpace) = c

# order, frequency

order(s::TensorSpace) = map(order, s.spaces)
order(s::TensorSpace, i::Int) = order(s.spaces[i])

frequency(s::TensorSpace) = map(frequency, s.spaces)
frequency(s::TensorSpace, i::Int) = frequency(s.spaces[i])

# promotion

Base.convert(::Type{T}, s::T) where {T<:TensorSpace} = s
Base.convert(::Type{TensorSpace{T}}, s::TensorSpace) where {T} =
    TensorSpace{T}(convert(T, s.spaces))

Base.promote_rule(::Type{T}, ::Type{T}) where {T<:TensorSpace} = T
Base.promote_rule(::Type{TensorSpace{T}}, ::Type{TensorSpace{S}}) where {T,S} =
    TensorSpace{promote_type(T, S)}

#

"""
    Taylor <: BaseSpace

Taylor sequence space whose elements are Taylor sequences of a prescribed order.

Fields:
- `order :: Int`
"""
struct Taylor <: BaseSpace
    order :: Int
    function Taylor(order::Int)
        order < 0 && return throw(DomainError(order, "Taylor is only defined for positive orders"))
        return new(order)
    end
end

order(s::Taylor) = s.order

Base.:(==)(s₁::Taylor, s₂::Taylor) = s₁.order == s₂.order
Base.issubset(s₁::Taylor, s₂::Taylor) = s₁.order ≤ s₂.order
Base.intersect(s₁::Taylor, s₂::Taylor) = Taylor(min(s₁.order, s₂.order))
Base.union(s₁::Taylor, s₂::Taylor) = Taylor(max(s₁.order, s₂.order))

dimension(s::Taylor) = s.order + 1
_firstindex(::Taylor) = 0
_lastindex(s::Taylor) = s.order
indices(s::Taylor) = 0:s.order

_findindex_constant(::Taylor) = 0

_findposition(i::Int, ::Taylor) = i + 1
_findposition(u::AbstractRange{Int}, ::Taylor) = u .+ 1
_findposition(u::AbstractVector{Int}, s::Taylor) = map(i -> _findposition(i, s), u)
_findposition(c::Colon, ::Taylor) = c

#

"""
    Fourier{T} <: BaseSpace

Fourier sequence space whose elements are Fourier sequences of a prescribed order and frequency.

Fields:
- `order :: Int`
- `frequency :: T`
"""
struct Fourier{T<:Real} <: BaseSpace
    order :: Int
    frequency :: T
    function Fourier{T}(order::Int, frequency::T) where {T<:Real}
        order < 0 && return throw(DomainError(order, "Fourier is only defined for positive orders"))
        return new{T}(order, frequency)
    end
end

Fourier(order::Int, frequency::T) where {T<:Real} = Fourier{T}(order, frequency)

order(s::Fourier) = s.order

frequency(s::Fourier) = s.frequency

Base.:(==)(s₁::Fourier, s₂::Fourier) = (s₁.frequency == s₂.frequency) & (s₁.order == s₂.order)
Base.issubset(s₁::Fourier, s₂::Fourier) = (s₁.frequency == s₂.frequency) & (s₁.order ≤ s₂.order)
function Base.intersect(s₁::Fourier{T}, s₂::Fourier{S}) where {T<:Real,S<:Real}
    s₁.frequency == s₂.frequency || return throw(DomainError)
    R = promote_type(T, S)
    return Fourier(min(s₁.order, s₂.order), convert(R, s₁.frequency))
end
function Base.union(s₁::Fourier{T}, s₂::Fourier{S}) where {T<:Real,S<:Real}
    s₁.frequency == s₂.frequency || return throw(DomainError)
    R = promote_type(T, S)
    return Fourier(max(s₁.order, s₂.order), convert(R, s₁.frequency))
end

dimension(s::Fourier) = 2s.order + 1
_firstindex(s::Fourier) = -s.order
_lastindex(s::Fourier) = s.order
indices(s::Fourier) = -s.order:s.order

_findindex_constant(::Fourier) = 0

_findposition(i::Int, s::Fourier) = i + s.order + 1
_findposition(u::AbstractRange{Int}, s::Fourier) = u .+ (s.order + 1)
_findposition(u::AbstractVector{Int}, s::Fourier) = map(i -> _findposition(i, s), u)
_findposition(c::Colon, ::Fourier) = c

# promotion

Base.convert(::Type{T}, s::T) where {T<:Fourier} = s
Base.convert(::Type{Fourier{T}}, s::Fourier) where {T<:Real} =
    Fourier{T}(s.order, convert(T, s.frequency))

Base.promote_rule(::Type{T}, ::Type{T}) where {T<:Fourier} = T
Base.promote_rule(::Type{Fourier{T}}, ::Type{Fourier{S}}) where {T<:Real,S<:Real} =
    Fourier{promote_type(T, S)}

#

"""
    Chebyshev <: BaseSpace

Chebyshev sequence space whose elements are Chebyshev sequences of a prescribed order.

Fields:
- `order :: Int`
"""
struct Chebyshev <: BaseSpace
    order :: Int
    function Chebyshev(order::Int)
        order < 0 && return throw(DomainError(order, "Chebyshev is only defined for positive orders"))
        return new(order)
    end
end

order(s::Chebyshev) = s.order

Base.:(==)(s₁::Chebyshev, s₂::Chebyshev) = s₁.order == s₂.order
Base.issubset(s₁::Chebyshev, s₂::Chebyshev) = s₁.order ≤ s₂.order
Base.intersect(s₁::Chebyshev, s₂::Chebyshev) = Chebyshev(min(s₁.order, s₂.order))
Base.union(s₁::Chebyshev, s₂::Chebyshev) = Chebyshev(max(s₁.order, s₂.order))

dimension(s::Chebyshev) = s.order + 1
_firstindex(::Chebyshev) = 0
_lastindex(s::Chebyshev) = s.order
indices(s::Chebyshev) = 0:s.order

_findindex_constant(::Chebyshev) = 0

_findposition(i::Int, ::Chebyshev) = i + 1
_findposition(u::AbstractRange{Int}, ::Chebyshev) = u .+ 1
_findposition(u::AbstractVector{Int}, s::Chebyshev) = map(i -> _findposition(i, s), u)
_findposition(c::Colon, ::Chebyshev) = c

# Cartesian spaces

"""
    CartesianSpace <: VectorSpace

Abstract type for all cartesian spaces.
"""
abstract type CartesianSpace <: VectorSpace end

_firstindex(::CartesianSpace) = 1
_lastindex(s::CartesianSpace) = dimension(s)
indices(s::CartesianSpace) = Base.OneTo(dimension(s))

_findposition(i, ::CartesianSpace) = i

_component_findposition(u::AbstractRange{Int}, s::CartesianSpace) =
    mapreduce(i -> _component_findposition(i, s), union, u)
_component_findposition(u::AbstractVector{Int}, s::CartesianSpace) =
    mapreduce(i -> _component_findposition(i, s), union, u)
_component_findposition(c::Colon, s::CartesianSpace) = c

"""
    CartesianPower{T<:VectorSpace} <: CartesianSpace

Cartesian space resulting from the cartesian products of a [`VectorSpace`](@ref).

Fields:
- `space :: T`
- `n :: Int`
"""
struct CartesianPower{T<:VectorSpace} <: CartesianSpace
    space :: T
    n :: Int
    function CartesianPower{T}(space::T, n::Int) where {T<:VectorSpace}
        n < 1 && return throw(DomainError(n, "CartesianPower is only defined for strictly positive integers"))
        return new{T}(space, n)
    end
end

CartesianPower(space::T, n::Int) where {T<:VectorSpace} =
    CartesianPower{T}(space, n)

space(s::CartesianPower) = s.space

spaces(s::CartesianPower) = fill(s.space, s.n)

nb_cartesian_product(s::CartesianPower) = s.n

Base.:^(s::VectorSpace, n::Int) = CartesianPower(s, n)

Base.@propagate_inbounds function Base.getindex(s::CartesianPower, i::Int)
    @boundscheck((1 ≤ i) & (i ≤ s.n) || throw(BoundsError(s, i)))
    return s.space
end
Base.@propagate_inbounds function Base.getindex(s::CartesianPower, u::AbstractRange{Int})
    @boundscheck((1 ≤ first(u)) & (last(u) ≤ s.n) || throw(BoundsError(s, u)))
    return CartesianPower(s.space, length(u))
end
Base.@propagate_inbounds function Base.getindex(s::CartesianPower, u::AbstractVector{Int})
    @boundscheck(all(i -> (1 ≤ i) & (i ≤ s.n), u) || throw(BoundsError(s, u)))
    return CartesianPower(s.space, length(u))
end
Base.@propagate_inbounds Base.getindex(s::CartesianPower, ::Colon) = s

#

Base.:(==)(s₁::CartesianPower, s₂::CartesianPower) =
    (s₁.n == s₂.n) & (s₁.space == s₂.space)
Base.issubset(s₁::CartesianPower, s₂::CartesianPower) =
    (s₁.n == s₂.n) & issubset(s₁.space, s₂.space)
function Base.intersect(s₁::CartesianPower, s₂::CartesianPower)
    s₁.n == s₂.n || return throw(DimensionMismatch("dimensions must match: s₁ has dimension $(s₁.n), s₂ has dimension $(s₂.n)"))
    return CartesianPower(intersect(s₁.space, s₂.space), s₁.n)
end
function Base.union(s₁::CartesianPower, s₂::CartesianPower)
    s₁.n == s₂.n || return throw(DimensionMismatch("dimensions must match: s₁ has dimension $(s₁.n), s₂ has dimension $(s₂.n)"))
    return CartesianPower(union(s₁.space, s₂.space), s₁.n)
end

dimension(s::CartesianPower) = dimension(s.space)*s.n
function dimension(s::CartesianPower, i::Int)
    (1 ≤ i) & (i ≤ s.n) || return throw(BoundsError(s, i))
    return dimension(s.space)
end
dimensions(s::CartesianPower) = fill(dimension(s.space), s.n)

# order, frequency

order(s::CartesianPower) = fill(order(s.space), s.n)
function order(s::CartesianPower, i::Int)
    (1 ≤ i) & (i ≤ s.n) || return throw(BoundsError(s, i))
    return order(s.space)
end

frequency(s::CartesianPower) = fill(frequency(s.space), s.n)
function frequency(s::CartesianPower, i::Int)
    (1 ≤ i) & (i ≤ s.n) || return throw(BoundsError(s, i))
    return frequency(s.space)
end

#

function _component_findposition(i::Int, s::CartesianPower)
    dim = dimension(s.space)
    x = (i-1)*dim
    return 1+x:dim+x
end
function _component_findposition(u::UnitRange{Int}, s::CartesianPower)
    dim = dimension(s.space)
    x = (first(u)-1)*dim
    return 1+x:dim*length(u)+x
end

# promotion

Base.convert(::Type{T}, s::T) where {T<:CartesianPower} = s
Base.convert(::Type{CartesianPower{T}}, s::CartesianPower) where {T} =
    CartesianPower{T}(convert(T, s.space), s.n)

Base.promote_rule(::Type{T}, ::Type{T}) where {T<:CartesianPower} = T
Base.promote_rule(::Type{CartesianPower{T}}, ::Type{CartesianPower{S}}) where {T,S} =
    CartesianPower{promote_type(T, S)}

"""
    CartesianProduct{T<:NTuple{N,VectorSpace} where {N}} <: CartesianSpace

Cartesian space resulting from `N` cartesian products of some [`VectorSpace`](@ref).

Fields:
- `spaces :: T`
"""
struct CartesianProduct{T<:NTuple{N,VectorSpace} where {N}} <: CartesianSpace
    spaces :: T
end

spaces(s::CartesianProduct) = s.spaces

nb_cartesian_product(s::CartesianProduct{<:NTuple{N,VectorSpace}}) where {N} = N

LinearAlgebra.:×(s₁::VectorSpace, s₂::VectorSpace) = CartesianProduct((s₁, s₂))
LinearAlgebra.:×(s₁::CartesianProduct, s₂::CartesianProduct) = CartesianProduct((s₁.spaces..., s₂.spaces...))
LinearAlgebra.:×(s₁::CartesianProduct, s₂::VectorSpace) = CartesianProduct((s₁.spaces..., s₂))
LinearAlgebra.:×(s₁::VectorSpace, s₂::CartesianProduct) = CartesianProduct((s₁, s₂.spaces...))

Base.@propagate_inbounds Base.getindex(s::CartesianProduct, i::Int) = getindex(s.spaces, i)
Base.@propagate_inbounds Base.getindex(s::CartesianProduct, u::AbstractRange{Int}) = CartesianProduct(getindex(s.spaces, u))
Base.@propagate_inbounds Base.getindex(s::CartesianProduct, u::AbstractVector{Int}) = CartesianProduct(getindex(s.spaces, u))
Base.@propagate_inbounds Base.getindex(s::CartesianProduct, c::Colon) = CartesianProduct(getindex(s.spaces, c))

Base.front(s::CartesianProduct) = CartesianProduct(Base.front(s.spaces))
Base.tail(s::CartesianProduct) = CartesianProduct(Base.tail(s.spaces))

#

function Base.:(==)(s₁::CartesianProduct{<:NTuple{N,VectorSpace}}, s₂::CartesianProduct{<:NTuple{N,VectorSpace}}) where {N}
    s₁[1] == s₂[1] && return Base.tail(s₁) == Base.tail(s₂)
    return false
end
Base.:(==)(s₁::CartesianProduct{<:Tuple{VectorSpace}}, s₂::CartesianProduct{<:Tuple{VectorSpace}}) =
    s₁[1] == s₂[1]
function Base.issubset(s₁::CartesianProduct{<:NTuple{N,VectorSpace}}, s₂::CartesianProduct{<:NTuple{N,VectorSpace}}) where {N}
    issubset(s₁[1], s₂[1]) && return issubset(Base.tail(s₁), Base.tail(s₂))
    return false
end
Base.issubset(s₁::CartesianProduct{<:Tuple{VectorSpace}}, s₂::CartesianProduct{<:Tuple{VectorSpace}}) =
    issubset(s₁[1], s₂[1])
Base.intersect(s₁::CartesianProduct{<:NTuple{N,VectorSpace}}, s₂::CartesianProduct{<:NTuple{N,VectorSpace}}) where {N} =
    CartesianProduct(map(intersect, s₁.spaces, s₂.spaces))
Base.union(s₁::CartesianProduct{<:NTuple{N,VectorSpace}}, s₂::CartesianProduct{<:NTuple{N,VectorSpace}}) where {N} =
    CartesianProduct(map(union, s₁.spaces, s₂.spaces))

dimension(s::CartesianProduct) = mapreduce(dimension, +, s.spaces)
dimension(s::CartesianProduct, i::Int) = dimension(s.spaces[i])
dimensions(s::CartesianProduct) = map(dimension, s.spaces)

# order, frequency

order(s::CartesianProduct) = map(order, s.spaces)
order(s::CartesianProduct, i::Int) = order(s.spaces[i])

frequency(s::CartesianProduct) = map(frequency, s.spaces)
frequency(s::CartesianProduct, i::Int) = frequency(s.spaces[i])

#

function _component_findposition(i::Int, s::CartesianProduct)
    dims = dimensions(s)
    dim = dims[i]
    x = mapreduce(j -> dims[j], +, 1:i-1; init=0)
    return 1+x:dim+x
end
function _component_findposition(u::UnitRange{Int}, s::CartesianProduct)
    dims = dimensions(s)
    dim = mapreduce(j -> dims[j], +, u)
    x = mapreduce(j -> dims[j], +, 1:first(u)-1; init=0)
    return 1+x:dim+x
end

# promotion

Base.convert(::Type{T}, s::T) where {T<:CartesianProduct} = s
Base.convert(::Type{CartesianProduct{T}}, s::CartesianProduct) where {T} =
    CartesianProduct{T}(convert(T, s.spaces))

Base.promote_rule(::Type{T}, ::Type{T}) where {T<:CartesianProduct} = T
Base.promote_rule(::Type{CartesianProduct{T}}, ::Type{CartesianProduct{S}}) where {T,S} =
    CartesianProduct{promote_type(T, S)}

#

_deep_nb_cartesian_product(::VectorSpace) = 1
_deep_nb_cartesian_product(s::CartesianPower) = s.n * _deep_nb_cartesian_product(s.space)
_deep_nb_cartesian_product(s::CartesianProduct) = sum(_deep_nb_cartesian_product, s.spaces)

#

_iscompatible(::VectorSpace, ::VectorSpace) = false
_iscompatible(::ParameterSpace, ::ParameterSpace) = true
_iscompatible(s₁::TensorSpace{<:NTuple{N,BaseSpace}}, s₂::TensorSpace{<:NTuple{N,BaseSpace}}) where {N} =
    @inbounds _iscompatible(s₁[1], s₂[1]) & _iscompatible(Base.tail(s₁), Base.tail(s₂))
_iscompatible(s₁::TensorSpace{<:Tuple{BaseSpace}}, s₂::TensorSpace{<:Tuple{BaseSpace}}) =
    @inbounds _iscompatible(s₁[1], s₂[1])
_iscompatible(::Taylor, ::Taylor) = true
_iscompatible(s₁::Fourier, s₂::Fourier) = frequency(s₁) == frequency(s₂)
_iscompatible(::Chebyshev, ::Chebyshev) = true
_iscompatible(s₁::CartesianPower, s₂::CartesianPower) =
    (nb_cartesian_product(s₁) == nb_cartesian_product(s₂)) & _iscompatible(space(s₁), space(s₂))
_iscompatible(s₁::CartesianProduct{<:NTuple{N,VectorSpace}}, s₂::CartesianProduct{<:NTuple{N,VectorSpace}}) where {N} =
    @inbounds _iscompatible(s₁[1], s₂[1]) & _iscompatible(Base.tail(s₁), Base.tail(s₂))
_iscompatible(s₁::CartesianProduct{<:Tuple{VectorSpace}}, s₂::CartesianProduct{<:Tuple{VectorSpace}}) =
    @inbounds _iscompatible(s₁[1], s₂[1])
_iscompatible(s₁::CartesianPower, s₂::CartesianProduct) =
    (nb_cartesian_product(s₁) == nb_cartesian_product(s₂)) & all(s₂ᵢ -> _iscompatible(space(s₁), s₂ᵢ), spaces(s₂))
_iscompatible(s₁::CartesianProduct, s₂::CartesianPower) =
    (nb_cartesian_product(s₁) == nb_cartesian_product(s₂)) & all(s₁ᵢ -> _iscompatible(s₁ᵢ, space(s₂)), spaces(s₁))

# show

Base.show(io::IO, ::MIME"text/plain", s::VectorSpace) = print(io, string_space(s))

string_space(::ParameterSpace) = "𝕂"

string_space(s::TensorSpace) = string_space(s[1]) * " ⨂ " * string_space(Base.tail(s))
string_space(s::TensorSpace{<:NTuple{2,BaseSpace}}) = string_space(s[1]) * " ⨂ " * string_space(s[2])
string_space(s::TensorSpace{<:Tuple{BaseSpace}}) = "TensorSpace(" * string_space(s[1]) * ")"
string_space(::TensorSpace{Tuple{}}) = "TensorSpace(())"

string_space(s::Taylor) = "Taylor(" * string(order(s)) * ")"
string_space(s::Fourier) = string(typeof(s)) * "(" * string(order(s)) * ", " * string(frequency(s)) * ")"
string_space(s::Chebyshev) = "Chebyshev(" * string(order(s)) * ")"

string_space(s::CartesianPower) = string_space(space(s)) * _supscript(nb_cartesian_product(s))
string_space(s::CartesianPower{<:TensorSpace}) = "(" * string_space(space(s)) * ")" * _supscript(nb_cartesian_product(s))
string_space(s::CartesianPower{<:CartesianSpace}) = "(" * string_space(space(s)) * ")" * _supscript(nb_cartesian_product(s))

string_space(s::CartesianProduct) = cartesian_string_space(s[1]) * " × " * string_space(Base.tail(s))
string_space(s::CartesianProduct{<:NTuple{2,VectorSpace}}) = cartesian_string_space(s[1]) * " × " * cartesian_string_space(s[2])
string_space(s::CartesianProduct{<:Tuple{VectorSpace}}) = "CartesianProduct(" * string_space(s[1]) * ")"
string_space(::CartesianProduct{Tuple{}}) = "CartesianProduct(())"
cartesian_string_space(s::VectorSpace) = string_space(s)
cartesian_string_space(s::TensorSpace) = "(" * string_space(s) * ")"
cartesian_string_space(s::CartesianProduct) = "(" * string_space(s) * ")"

function _supscript_digit(i::Int)
    if i == 0
        return "⁰"
    elseif i == 1
        return "¹"
    elseif i == 2
        return "²"
    elseif i == 3
        return "³"
    elseif i == 4
        return "⁴"
    elseif i == 5
        return "⁵"
    elseif i == 6
        return "⁶"
    elseif i == 7
        return "⁷"
    elseif i == 8
        return "⁸"
    else
        return "⁹"
    end
end
function _supscript(n::Int)
    x = ""
    for i ∈ reverse!(digits(n))
        x *= _supscript_digit(i)
    end
    return x
end


text:
# Estimate for the number of weights that don't underflow.
# This is only an estimate, not an upper bound, so code has to be able to cope
# with the estimate being off.
# The estimate here is similar to the Laguerre one, but more general with the
# additional argument m.
estimate_freud_reduced_n(n, alpha, m) = round(typeof(n), min(exp(exp(1/m)*1.05)*n^(1-1/2/m), n))

"""
Compute the n-point Gaussian quadrature rule for the Freud-type weight function
`x^α exp(-q_m x^m)`.
The computations are based on asymptotic expansions for the corresponding orthogonal
polynomials.
"""
function asy_gaussfreud(n::Integer, alpha = 0.0, m = 1, qm = 1.0; reduced = false)
    if alpha <= -1
        error("The parameter α <= -1 corresponds to a nonintegrable weight function")
    end
    if n < 0
        error("gaussfreud($n,$alpha,$m,$qm) not defined: n must be positive.")
    end
    if alpha^2/n > 1
        warn("A large alpha may lead to inaccurate results because the weight is low and R(z) is not close to identity.")
    end

    ELT = promote_type(typeof(float(alpha)), typeof(qm))

    T = ceil(Int, 34/log(n) ) # Heuristic for number of terms, should be scaled by the logarithm of eps(ELT) over the machine precision.
    UQ0 = getUQ(alpha, qm, m, T)
    n_alloc = reduced ? estimate_freud_reduced_n(n, alpha, m) : n

    n_pre = max(ceil(Int, 3.6*n^0.188), 7)
    # Heuristics to switch between Bessel, extrapolation and Airy initial guesses.
    igatt = ceil(Int, n_alloc + 1.31*n^0.4 - n)

    A = zeros(ELT, m+1)
    for k in 0:m
        A[k+1] = prod((2*(1:k) .- 1)/2 ./ (1:k))
    end
    softEdge = (n*2/m/qm/A[m+1] )^(1/m)
    # Use finite differences for derivative of polynomial when not x^alpha*exp(-x) and use other initial approximations
    useFinDiff = (m != 1) || (qm != 1.0)
    bes = besselroots(alpha, n_pre).^2 # [Tricomi 1947 pg. 296]
    w = zeros(ELT, n_alloc)

    # Find initial values for x
    if useFinDiff
        x = [bes*(2*m-1)^2/16/m^2/n^2*softEdge ; zeros(ELT, n_alloc-n_pre) ]
    else
        ak = [-13.69148903521072; -12.828776752865757; -11.93601556323626;    -11.00852430373326; -10.04017434155809; -9.02265085340981; -7.944133587120853;    -6.786708090071759; -5.520559828095551; -4.08794944413097; -2.338107410459767]
        t = 3*pi/2*( (igatt:-1:12).-0.25) # [DLMF (9.9.6)]
        ak = [-t.^(2/3).*(1 .+ 5/48 ./ t.^2 - 5/36 ./t.^4 .+ 77125/82944 ./t.^6 .- 10856875/6967296 ./t.^8); ak[max(1,12-igatt):11] ]
        nu = 4*n+2*alpha+2 # [Gatteshi 2002 (4.9)]
        air = (nu .+ ak*(4*nu)^(1/3) .+ ak.^2*(nu/16)^(-1/3)/5 .+ (11/35 .- alpha^2 .- 12/175*ak.^3)/nu + (16/1575*ak .+ 92/7875*ak.^4)*2^(2/3)*nu^(-5/3) .- (15152/3031875*ak.^5 .+ 1088/121275*ak.^2)*2^(1/3)*nu^(-7/3))
        x = [ bes/(4*n + 2*alpha+2) .* (1 .+ (bes .+ 2*(alpha^2 - 1) )/(4*n + 2*alpha+2)^2/3 ) ; zeros(ELT, n_alloc - n_pre -max(igatt,0) ) ; air]
    end

    if !useFinDiff
        UQ1 = getUQ(alpha+1, qm, m, T)
        factor0 = 1-4im*4^alpha*sum((UQ0[1,2,1:(T-1),1,2] + UQ0[1,2,1:(T-1),1,1])./n.^reshape(1:(T-1), (1,1,T-1)) )
        factor1 = 1-4im*4^(alpha+1)*sum((UQ1[1,2,1:(T-1),1,2] + UQ1[1,2,1:(T-1),1,1])./n.^reshape(1:(T-1), (1,1,T-1)) )
        factorx = real(sqrt(factor0/factor1 )/2/(1 - 1/n)^(1+alpha/2))
        factorw = real( -(1 - 1/(n + 1) )^(n + 1+ alpha/2)*(1 - 1/n)^(1 + alpha/2)*exp(1 + 2*log(2) )*4^(1+alpha)*pi*n^alpha*sqrt(factor0*factor1)*(1 + 1/n)^(alpha/2) )
    end

    noUnderflow = true
    for k = 1:n
        if useFinDiff && (k > n_pre)
            # Use linear extrapolation for the initial guesses for robustness in generalised weights.
            x[k] = 2*x[k-1] -x[k-2]
        elseif k > n_pre
            # Use sextic extrapolation for the initial guesses.
            x[k] = 7*x[k-1] -21*x[k-2] +35*x[k-3] -35*x[k-4] +21*x[k-5] -7*x[k-6] +x[k-7]
        end
        step = x[k]
        l = 0 # Newton-Raphson iteration number
        max_iter = 20
        ov = floatmax(ELT) # Previous/old value
        ox = x[k] # Old x
        # Accuracy of the expansions up to machine precision would lower this bound.
        while ( abs(step) > eps(ELT)*40*x[k] ) && ( l < max_iter)
            l = l + 1
            pe = polyAsyRHgen(n, x[k], alpha, T, qm, m, UQ0)
            if abs(pe) >= abs(ov)*(1-35*eps(ELT))
                # The function values do not decrease enough any more due to roundoff errors.
                x[k] = ox # Set to the previous value and quit.
                break
            end
            if useFinDiff
                hh = max(sqrt(eps(ELT))*x[k], sqrt(eps(ELT)) )
                step = pe*hh/(polyAsyRHgen(n, x[k]+hh, alpha, T, qm, m, UQ0) - pe)
            else
                # poly' = (p*exp(-Q/2) )' = exp(-Q/2)*(p' -p/2) with orthonormal p.
                step = pe/(polyAsyRHgen(n-1, x[k], alpha+1, T, qm, m, UQ1)*factorx - pe/2)
            end
            ox = x[k]
            x[k] = x[k] -step
            ov = pe
        end
        if (x[k] < 0) || (l == max_iter) || ( (k != 1) && (x[k-1] >= x[k]) ) || isnan(x[k])
            # Print some debugging information and throw an error.
            print(x[k], " = x[k], k = ", k, ", l = ", l, ", x[k-1] = ", x[k-1], ", x[k-2] = ", x[k-1], ", step = ", step, ", ox=", ox, ", ov = ", ov, ".\n")
            error("Newton method may not have converged.")
        elseif x[k] > softEdge
            warn("Node is outside the support of the measure: inaccuracy is expected.")
        end
        if noUnderflow
            if useFinDiff
                hh = max(sqrt(eps(ELT))*x[k], sqrt(eps(ELT)) )
                w[k] = hh/(polyAsyRHgen(n,x[k]+hh,alpha,T,qm,m,UQ0) -polyAsyRHgen(n,x[k],alpha,T,qm,m,UQ0))/polyAsyRHgen(n-1,x[k],alpha,T,qm,m,UQ0)/exp(qm*x[k]^m) # This leaves out a constant factor, given by a ratio of leading order coefficients and normalising constants
            else noUnderflow
                w[k] = factorw/polyAsyRHgen(n-1, x[k], alpha+1, T, qm, m, UQ1)/polyAsyRHgen(n+1, x[k], alpha, T,qm,m,UQ0)/exp( x[k] )
            end
        end
        if noUnderflow && (w[k] < underflow_threshold(ELT))
            # Frome here on after weights will no longer be computed
            noUnderflow = false
        end
        if reduced
            if (k > 1) && !noUnderflow
                x = x[1:k-1]
                w = w[1:k-1]
                return x, w
            end
            if k == n_alloc
                # We have to allocate a bigger array
                n_alloc *= 2
                x1 = x
                w1 = w
                x = zeros(T, n_alloc)
                w = zeros(T, n_alloc)
                x[1:k] = x1
                w[1:k] = w1
            end
        end
    end
    x, w
end