text:
module Fitting

using GenericSVD
using Test
using Unitful: EnergyUnits, PressureUnits, ustrip, @u_str
using UnitfulAtomic
using YAML

using EquationsOfStateOfSolids:
    Parameters,
    Murnaghan1st,
    BirchMurnaghan2nd,
    BirchMurnaghan3rd,
    BirchMurnaghan4th,
    PoirierTarantola3rd,
    Vinet,
    EnergyEquation
using EquationsOfStateOfSolids.Fitting: nonlinfit, linfit

import Unitful

Unitful.promote_unit(::S, ::T) where {S<:EnergyUnits,T<:EnergyUnits} = u"eV"
Unitful.promote_unit(::S, ::T) where {S<:PressureUnits,T<:PressureUnits} = u"GPa"

# Do not export! Only for internal use!
function _isapprox(a::T, b::T; kwargs...) where {T<:Parameters}
    x, y = [], []
    for f in fieldnames(T)
        pa, pb = map(ustrip, promote(getfield(a, f), getfield(b, f)))
        push!(x, pa)
        push!(y, pb)
    end
    return isapprox(x, y; kwargs...)
end

# Data from https://github.com/materialsproject/pymatgen/blob/19c4d98/pymatgen/analysis/tests/test_eos.py#L17-L73
@testset "Test data from Pymatgen" begin
    data = open("data/si.yml", "r") do io
        YAML.load(io)
    end
    @testset "With `Float64`" begin
        volumes, energies = data["volume"], data["energy"]
        @test _isapprox(
            nonlinfit(EnergyEquation(BirchMurnaghan3rd(40, 0.5, 4)), volumes, energies),
            BirchMurnaghan3rd(
                40.98926572528106,
                0.5369258245417454,
                4.178644235500821,
                -10.842803908240892,
            ),
        )
        @test _isapprox(
            linfit(EnergyEquation(BirchMurnaghan3rd(40, 0.5, 4)), volumes, energies),
            BirchMurnaghan3rd(
                40.98926572528106,
                0.5369258245417454,
                4.178644235500821,
                -10.842803908240892,
            ),
        )
        @test _isapprox(
            nonlinfit(EnergyEquation(Murnaghan1st(41, 0.5, 4)), volumes, energies),
            Murnaghan1st(
                41.13757930387086,
                0.5144967693786603,
                3.9123862262572264,
                -10.836794514626673,
            ),
        )
        @test _isapprox(
            nonlinfit(EnergyEquation(PoirierTarantola3rd(41, 0.5, 4)), volumes, energies),
            PoirierTarantola3rd(
                40.86770643373908,
                0.5667729960804602,
                4.331688936974368,
                -10.851486685041658,
            ),
        )
        @test _isapprox(
            linfit(EnergyEquation(PoirierTarantola3rd(41, 0.5, 4)), volumes, energies),
            PoirierTarantola3rd(
                40.86770643373908,
                0.5667729960804602,
                4.331688936974368,
                -10.851486685041658,
            ),
        )
        @test _isapprox(
            nonlinfit(EnergyEquation(Vinet(41, 0.5, 4)), volumes, energies),
            Vinet(
                40.916875663779784,
                0.5493839425156859,
                4.3051929654936885,
                -10.846160810560756,
            ),
        )
    end

    @testset "`linfit` for `BigFloat` parameters" begin
        volumes, energies = data["volume"], data["energy"]
        @test _isapprox(
            linfit(
                EnergyEquation(BirchMurnaghan3rd(40, 0.5, 4)),
                big.(volumes),
                big.(energies),
            ),
            BirchMurnaghan3rd{BigFloat}(
                40.98926572528106,
                0.5369258245417454,
                4.178644235500821,
                -10.842803908240892,
            );
            atol = 1e-8,
        )
        @test _isapprox(
            linfit(EnergyEquation(BirchMurnaghan3rd{BigInt}(40, 1, 4)), volumes, energies),
            BirchMurnaghan3rd{BigFloat}(
                40.98926572528106,
                0.5369258245417454,
                4.178644235500821,
                -10.842803908240892,
            );
            atol = 1e-8,
        )
        @test _isapprox(
            linfit(
                EnergyEquation(PoirierTarantola3rd(41, 0.5, 4)),
                big.(volumes),
                big.(energies),
            ),
            PoirierTarantola3rd{BigFloat}(
                40.86770643373908,
                0.5667729960804602,
                4.331688936974368,
                -10.851486685041658,
            );
            atol = 1e-8,
        )
        @test _isapprox(
            linfit(
                EnergyEquation(PoirierTarantola3rd{BigFloat}(41, 1 // 2, 4)),
                volumes,
                energies,
            ),
            PoirierTarantola3rd{BigFloat}(
                40.86770643373908,
                0.5667729960804602,
                4.331688936974368,
                -10.851486685041658,
            );
            atol = 1e-8,
        )
    end
end

# Data from https://github.com/materialsproject/pymatgen/blob/19c4d98/pymatgen/analysis/tests/test_eos.py#L92-L167
@testset "Test Mg dataset" begin
    data = open("data/mp153.yml", "r") do io
        YAML.load(io)
    end
    @testset "without unit" begin
        volumes, energies, known_energies_vinet =
            data["volume"], data["energy"], data["known_energy_vinet"]
        fitted_eos = nonlinfit(EnergyEquation(Vinet(23, 0.5, 4)), volumes, energies)
        @test _isapprox(
            fitted_eos,
            Vinet(22.957645593, 0.22570911414, 4.06054339, -1.59442926062),
        )
        @test isapprox(
            map(EnergyEquation(fitted_eos), volumes),
            known_energies_vinet;
            atol = 1e-6,
        )
    end

    @testset "with units" begin
        volumes, energies, known_energies_vinet = data["volume"] * u"angstrom^3",
        data["energy"] * u"eV",
        data["known_energy_vinet"] * u"eV"
        fitted_eos = nonlinfit(
            EnergyEquation(Vinet(23u"angstrom^3", 36.16u"GPa", 4)),
            volumes,
            energies,
        )
        @test _isapprox(
            fitted_eos,
            Vinet(
                22.957645593u"angstrom^3",
                36.16258657649159u"GPa",  # https://github.com/materialsproject/pymatgen/blob/19c4d98/pymatgen/analysis/tests/test_eos.py#L181
                4.06054339,
                -1.59442926062u"eV",
            ),
        )
        @test isapprox(
            map(EnergyEquation(fitted_eos), volumes),
            known_energies_vinet;
            atol = 1e-6u"eV",
        )
    end
end

# Data from https://github.com/materialsproject/pymatgen/blob/19c4d98/pymatgen/analysis/tests/test_eos.py#L185-L260
@testset "Test Si dataset" begin
    data = open("data/mp149.yml", "r") do io
        YAML.load(io)
    end
    @testset "without unit" begin
        volumes, energies, known_energies_vinet =
            data["volume"], data["energy"], data["known_energy_vinet"]
        fitted_eos = nonlinfit(EnergyEquation(Vinet(20, 0.5, 4)), volumes, energies)
        @test _isapprox(
            fitted_eos,
            Vinet(20.446696754, 0.55166385214, 4.32437391, -5.42496338987),
        )
        @test isapprox(
            map(EnergyEquation(fitted_eos), volumes),
            known_energies_vinet;
            atol = 1e-5,
        )
    end

    @testset "with units" begin
        volumes, energies, known_energies_vinet = data["volume"] * u"angstrom^3",
        data["energy"] * u"eV",
        data["known_energy_vinet"] * u"eV"
        fitted_eos = nonlinfit(
            EnergyEquation(Vinet(20u"angstrom^3", 88.39u"GPa", 4)),
            volumes,
            energies,
        )
        @test _isapprox(
            fitted_eos,
            Vinet(
                20.446696754u"angstrom^3",
                88.38629264585195u"GPa",  # https://github.com/materialsproject/pymatgen/blob/19c4d98/pymatgen/analysis/tests/test_eos.py#L274
                4.32437391,
                -5.42496338987u"eV",
            ),
        )
        @test isapprox(
            map(EnergyEquation(fitted_eos), volumes),
            known_energies_vinet;
            atol = 1e-6u"eV",
        )
    end
end

# Data from https://github.com/materialsproject/pymatgen/blob/19c4d98/pymatgen/analysis/tests/test_eos.py#L278-L353
@testset "Test Ti dataset" begin
    data = open("data/mp72.yml", "r") do io
        YAML.load(io)
    end
    @testset "without unit" begin
        volumes, energies, known_energies_vinet =
            data["volume"], data["energy"], data["known_energy_vinet"]
        fitted_eos = nonlinfit(EnergyEquation(Vinet(17, 0.5, 4)), volumes, energies)
        @test _isapprox(
            fitted_eos,
            Vinet(17.1322302613, 0.70297662247, 3.638807756, -7.89741495912),
        )
        @test isapprox(
            map(EnergyEquation(fitted_eos), volumes),
            known_energies_vinet;
            atol = 1e-5,
        )
    end

    @testset "with units" begin
        volumes, energies, known_energies_vinet = data["volume"] * u"angstrom^3",
        data["energy"] * u"eV",
        data["known_energy_vinet"] * u"eV"
        fitted_eos = nonlinfit(
            EnergyEquation(Vinet(17u"angstrom^3", 112.63u"GPa", 4)),
            volumes,
            energies,
        )
        @test _isapprox(
            fitted_eos,
            Vinet(
                17.1322302613u"angstrom^3",
                112.62927094503254u"GPa",  # https://github.com/materialsproject/pymatgen/blob/19c4d98/pymatgen/analysis/tests/test_eos.py#L367
                3.638807756,
                -7.89741495912u"eV",
            ),
        )
        @test isapprox(
            map(EnergyEquation(fitted_eos), volumes),
            known_energies_vinet;
            atol = 1e-7u"eV",
        )
    end
end

@testset "Test `w2k-lda-na.dat` from `Gibbs2`" begin
    data = open("data/w2k-lda-na.yml", "r") do io
        YAML.load(io)
    end
    @testset "without unit" begin
        volumes, energies = data["volume"], data["energy"]
        # See https://github.com/aoterodelaroza/asturfit/blob/4de9b41/test/test03.out#L117-L122
        @test _isapprox(
            nonlinfit(EnergyEquation(Murnaghan1st(224, 6e-4, 4)), volumes, energies),
            Murnaghan1st(224.501825, 6.047952315268776e-4, 3.723835, -323.417686);
            atol = 1e-5,
        )
        # No reference data, I run on my computer.
        @test _isapprox(
            nonlinfit(EnergyEquation(BirchMurnaghan2nd(224, 6e-4)), volumes, energies),
            BirchMurnaghan2nd(223.7192539523166, 6.268341030294978e-4, -323.4177121144877);
            atol = 1e-8,
        )
        # No reference data, I run on my computer.
        @test _isapprox(
            linfit(EnergyEquation(BirchMurnaghan2nd(224, 6e-4)), volumes, energies),
            BirchMurnaghan2nd(223.7192539523166, 6.268341030294978e-4, -323.4177121144877),
        )
        # See https://github.com/aoterodelaroza/asturfit/blob/4de9b41/test/test03.out#L15-L20
        @test _isapprox(
            nonlinfit(EnergyEquation(BirchMurnaghan3rd(224, 6e-4, 4)), volumes, energies),
            BirchMurnaghan3rd(224.444565, 6.250619105057268e-4, 3.740369, -323.417714),
        )
        @test _isapprox(
            linfit(EnergyEquation(BirchMurnaghan3rd(224, 6e-4, 4)), volumes, energies),
            BirchMurnaghan3rd(224.444565, 6.250619105057268e-4, 3.740369, -323.417714),
        )
        # See https://github.com/aoterodelaroza/asturfit/blob/4de9b41/test/test03.out#L30-L36
        @test _isapprox(
            nonlinfit(
                EnergyEquation(BirchMurnaghan4th(224, 6e-4, 4, -5460)),  # bohr^3, Ry/bohr^3, 1, bohr^3/Ry, Ry
                volumes,
                energies,
            ),
            BirchMurnaghan4th(
                224.457562,  # bohr^3
                6.229381129795094e-4,  # Ry/bohr^3
                3.730992,
                -5322.7030547560435,  # bohr^3/Ry
                -323.417712,  # Ry
            );
            rtol = 1e-5,
        )
        # See https://github.com/aoterodelaroza/asturfit/blob/4de9b41/test/test03.out#L30-L36
        # @test _isapprox(
        #     linfit(
        #         EnergyEquation(BirchMurnaghan4th(224, 6e-4, 4, -5460)),  # bohr^3, Ry/bohr^3, 1, bohr^3/Ry, Ry
        #         volumes,
        #         energies,
        #     ),
        #     BirchMurnaghan4th(
        #         224.457562,  # bohr^3
        #         6.229381129795094e-4,  # Ry/bohr^3
        #         3.730992,
        #         -5322.7030547560435,  # bohr^3/Ry
        #         -323.417712,  # Ry
        #     );
        #     rtol = 1e-5,
        # )
        # See https://github.com/aoterodelaroza/asturfit/blob/4de9b41/test/test03.out#L98-L105
        # @test _isapprox(
        #     nonlinfit(
        # EnergyEquation(        #         BirchMurnaghan5th(224.445371, 0.0006, 4, -5500, 3.884535907971559e7, -323)),
        #         volumes,
        #         energies,
        #     ),
        #     BirchMurnaghan5th(
        #         224.451813,
        #         0.0006228893043314733,
        #         3.736723,
        #         -5292.414119096362,
        #         6.3542116050611705e7,
        #         -323.417712,
        #     );
        #     atol = 1,
        # )  # FIXME: result is wrong
        # # No reference data, I run on my computer.
        @test _isapprox(
            nonlinfit(EnergyEquation(Vinet(224, 6e-4, 4)), volumes, energies),
            Vinet(
                224.45278665796354,
                6.313500637481759e-4,
                3.7312381477678853,
                -323.4177229576912,
            ),
        )
        # See https://github.com/aoterodelaroza/asturfit/blob/4de9b41/test/test03.out#L66-L71
        @test _isapprox(
            nonlinfit(EnergyEquation(PoirierTarantola3rd(224, 6e-4, 4)), volumes, energies),
            PoirierTarantola3rd(224.509208, 6.3589226415983795e-4, 3.690448, -323.41773);
            atol = 1e-5,
        )
        # See https://github.com/aoterodelaroza/asturfit/blob/4de9b41/test/test03.out#L66-L71
        @test _isapprox(
            linfit(EnergyEquation(PoirierTarantola3rd(224, 6e-4, 4)), volumes, energies),
            PoirierTarantola3rd(224.509208, 6.3589226415983795e-4, 3.690448, -323.41773);
            atol = 1e-5,
        )
        # # FIXME: This cannot go through
        # @test _isapprox(
        #     nonlinfit(
        # EnergyEquation(        #         PoirierTarantola4th(220, 0.0006, 3.7, -5500, -323)),
        #         volumes,
        #         energies,
        #     ),
        #     PoirierTarantola4th(
        #         224.430182,
        #         0.0006232241765069493,
        #         3.758360,
        #         -5493.859729817176,
        #         -323.417712,
        #     ),
        # )
        # See https://github.com/aoterodelaroza/asturfit/blob/4de9b41/test/test03.out#L98-L105
        # @test _isapprox(
        #     nonlinfit(
        #         EnergyEquation(PoirierTarantola5th(
        #             224.445371,
        #             0.0006,
        #             3.8,
        #             -5500,
        #             6e7,
        #             -323,
        #         )),
        #         volumes,
        #         energies,
        #     ),
        #     PoirierTarantola5th(
        #         224.451250,
        #         0.000622877204137392,
        #         3.737484,
        #         -5283.999708607125,
        #         6.296000262990379e7,
        #         -323.417712,
        #     );
        #     rtol = 0.05,
        # )
    end

    @testset "with units" begin
        volumes, energies = data["volume"] * u"bohr^3", data["energy"] * u"Ry"
        # See https://github.com/aoterodelaroza/asturfit/blob/4de9b41/test/test03.out#L117-L122
        @test _isapprox(
            nonlinfit(
                EnergyEquation(Murnaghan1st(224u"bohr^3", 9u"GPa", 4)),
                volumes,
                energies,
            ),
            Murnaghan1st(224.501825u"bohr^3", 8.896845u"GPa", 3.723835, -323.417686u"Ry");
            atol = 1e-5,
        )
        # See https://github.com/aoterodelaroza/asturfit/blob/4de9b41/test/test03.out#L15-L20
        @test _isapprox(
            nonlinfit(
                EnergyEquation(BirchMurnaghan3rd(224u"bohr^3", 9u"GPa", 4)),
                volumes,
                energies,
            ),
            BirchMurnaghan3rd(
                224.444565u"bohr^3",
                9.194978u"GPa",
                3.740369,
                -323.417714u"Ry",
            ),
        )
        # See https://github.com/aoterodelaroza/asturfit/blob/4de9b41/test/test03.out#L15-L20
        @test _isapprox(
            linfit(
                EnergyEquation(BirchMurnaghan3rd(224u"bohr^3", 9u"GPa", 4)),
                volumes,
                energies,
            ),
            BirchMurnaghan3rd(
                224.444565u"bohr^3",
                6.250619009766436e-4u"Ry/bohr^3",
                3.740369,
                -323.417714u"Ry",
            ),
        )
        # See https://github.com/aoterodelaroza/asturfit/blob/4de9b41/test/test03.out#L30-L36
        @test _isapprox(
            nonlinfit(
                EnergyEquation(BirchMurnaghan4th(224u"bohr^3", 9u"GPa", 4, -0.3u"1/GPa")),
                volumes,
                energies,
            ),
            BirchMurnaghan4th(
                224.457562u"bohr^3",
                9.163736u"GPa",
                3.730992,
                -0.361830u"1/GPa",
                -323.417712u"Ry",  # Ry
            );
            rtol = 1e-7,
        )
        # See https://github.com/aoterodelaroza/asturfit/blob/4de9b41/test/test03.out#L66-L71
        @test _isapprox(
            nonlinfit(
                EnergyEquation(PoirierTarantola3rd(224u"bohr^3", 9u"GPa", 4)),
                volumes,
                energies,
            ),
            PoirierTarantola3rd(
                224.509208u"bohr^3",
                9.354298u"GPa",
                3.690448,
                -323.41773u"Ry",
            ),
        )
        # See https://github.com/aoterodelaroza/asturfit/blob/4de9b41/test/test03.out#L66-L71
        @test _isapprox(
            linfit(
                EnergyEquation(PoirierTarantola3rd(224u"bohr^3", 9u"GPa", 4)),
                volumes,
                energies,
            ),
            PoirierTarantola3rd(
                224.509208u"bohr^3",
                6.3589226415983795e-4u"Ry/bohr^3",
                3.690448,
                -323.41773u"Ry",
            );
            atol = 1e-5,
        )
    end
end

@testset "Test `w2k-lda-k.dat` from `Gibbs2`" begin
    data = open("data/w2k-lda-k.yml", "r") do io
        YAML.load(io)
    end
    @testset "without unit" begin
        volumes, energies = data["volume"], data["energy"]
        @test _isapprox(
            nonlinfit(EnergyEquation(Murnaghan1st(430, 3e-4, 4)), volumes, energies),
            Murnaghan1st(
                435.05782299050884,
                2.8297159355249786e-4,
                3.5705032675000785,
                -1201.2082739321822,
            ),
        )
        @test _isapprox(
            nonlinfit(EnergyEquation(BirchMurnaghan2nd(430, 3e-4)), volumes, energies),
            BirchMurnaghan2nd(430.10027687726716, 3.02451215462375e-4, -1201.2083221436026),
        )
        @test _isapprox(
            nonlinfit(EnergyEquation(BirchMurnaghan3rd(430, 3e-4, 4)), volumes, energies),
            BirchMurnaghan3rd(
                432.67139080209046,
                3.0508544859901674e-4,
                3.7894868450211923,
                -1201.2083959943404,
            ),
        )
        @test _isapprox(
            nonlinfit(
                EnergyEquation(BirchMurnaghan4th(432, 3e-4, 3.8, -11773)),
                volumes,
                energies,
            ),
            BirchMurnaghan4th(
                432.8012195854224,
                3.041889904166284e-4,
                3.774020919355492,
                -11773.192574765615,
                -1201.2083912308235,
            );
            rtol = 1e-4,
        )
        @test _isapprox(
            nonlinfit(EnergyEquation(Vinet(432, 3e-4, 3.8)), volumes, energies),
            Vinet(
                432.04609865398015,
                3.137631070690569e-4,
                3.837666939407128,
                -1201.2084453225773,
            ),
        )
    end

    @testset "with units" begin
        volumes, energies = data["volume"] * u"bohr^3", data["energy"] * u"Ry"
        @test _isapprox(
            nonlinfit(
                EnergyEquation(BirchMurnaghan3rd(430u"bohr^3", 3e-4u"Ry/bohr^3", 4)),
                volumes,
                energies,
            ),
            BirchMurnaghan3rd(
                432.6713907942206u"bohr^3",
                3.0508544859901674e-4u"Ry/bohr^3",
                3.789486849598267,
                -1201.208395994332u"Ry",
            ),
        )
        @test _isapprox(
            linfit(
                EnergyEquation(BirchMurnaghan3rd(430u"bohr^3", 3e-4u"Ry/bohr^3", 4)),
                volumes,
                energies,
            ),
            BirchMurnaghan3rd(
                432.6713907942206u"bohr^3",
                3.0508544859901674e-4u"Ry/bohr^3",
                3.789486849598267,
                -1201.208395994332u"Ry",
            ),
        )
    end
end

# Data from http://tutorials.crystalsolutions.eu/tutorial.html?td=eos&tf=eos_tut
@testset "Test MgO data" begin
    volumes = [17.789658, 18.382125, 18.987603, 19.336585, 19.606232, 20.238155, 20.883512]
    energies = [
        -275.3023029137,
        -275.3039525749,
        -275.3047829172,
        -275.3049167772,
        -275.3048617406,
        -275.3042539881,
        -275.3030141001,
    ]
    @testset "without unit" begin
        @test _isapprox(
            nonlinfit(EnergyEquation(Murnaghan1st(19, 100, 4)), volumes, energies),
            Murnaghan1st(19.3620, 0.035786498967406696, 3.70, -275.30491639);
            atol = 1e-2,
        )
        @test _isapprox(
            nonlinfit(EnergyEquation(BirchMurnaghan3rd(19, 100, 4)), volumes, energies),
            BirchMurnaghan3rd(19.3618, 0.035859897760315104, 3.72, -275.30491678);
            atol = 1e-2,
        )
        @test _isapprox(
            nonlinfit(EnergyEquation(PoirierTarantola3rd(19, 100, 4)), volumes, energies),
            PoirierTarantola3rd(19.3617, 0.0358988908690477, 3.73, -275.30491699);
            atol = 1e-2,
        )
        @test _isapprox(
            nonlinfit(EnergyEquation(Vinet(19, 100, 4)), volumes, energies),
            Vinet(19.3617, 0.035880541170820596, 3.73, -275.30491690);
            atol = 1e-4,
        )
    end

    @testset "with units" begin
        volumes *= u"angstrom^3"
        energies *= u"Eh_au"
        @test _isapprox(
            nonlinfit(
                EnergyEquation(Murnaghan1st(19u"angstrom^3", 100u"GPa", 4)),
                volumes,
                energies,
            ),
            Murnaghan1st(19.3620u"angstrom^3", 156.02u"GPa", 3.70, -275.30491639u"Eh_au");
            rtol = 1e-6,
        )
        @test _isapprox(
            nonlinfit(
                EnergyEquation(BirchMurnaghan3rd(19u"angstrom^3", 100u"GPa", 4)),
                volumes,
                energies,
            ),
            BirchMurnaghan3rd(
                19.3618u"angstrom^3",
                156.34u"GPa",
                3.72,
                -275.30491678u"Eh_au",
            );
            rtol = 1e-6,
        )
        @test _isapprox(
            nonlinfit(
                EnergyEquation(PoirierTarantola3rd(19u"angstrom^3", 100u"GPa", 4)),
                volumes,
                energies,
            ),
            PoirierTarantola3rd(
                19.3617u"angstrom^3",
                156.51u"GPa",
                3.73,
                -275.30491699u"Eh_au",
            );
            rtol = 1e-6,
        )
        @test _isapprox(
            nonlinfit(
                EnergyEquation(Vinet(19u"angstrom^3", 100u"GPa", 4)),
                volumes,
                energies,
            ),
            Vinet(19.3617u"angstrom^3", 156.43u"GPa", 3.73, -275.30491690u"Eh_au");
            rtol = 1e-6,
        )
        @test _isapprox(
            linfit(
                EnergyEquation(BirchMurnaghan3rd(19u"angstrom^3", 100u"GPa", 4)),
                volumes,
                energies,
            ),
            BirchMurnaghan3rd(
                19.3618u"angstrom^3",
                0.035859897760315104u"Eh_au / angstrom^3",
                3.72,
                -275.30491678u"Eh_au",
            );
            rtol = 1e-4,
        )
        @test _isapprox(
            linfit(
                EnergyEquation(PoirierTarantola3rd(19u"angstrom^3", 150u"GPa", 3.7)),
                volumes,
                energies,
            ),
            PoirierTarantola3rd(
                19.3617u"angstrom^3",
                0.0358988908690477u"Eh_au / angstrom^3",
                3.73,
                -275.30491699u"Eh_au",
            );
            rtol = 1e-5,
        )
    end
end

end


text:
##############################################################################
##
## Object constructed by the user
##
##############################################################################

# Object constructed by the user
type InteractiveFixedEffectModel
    id::Symbol
    time::Symbol
    rank::Int64
end

abstract AbstractFactorModel{T}
abstract AbstractFactorSolution{T}
##############################################################################
##
## Factor Model
##
##############################################################################

type FactorModel{Rank, W, Rid, Rtime} <: AbstractFactorModel{Rank}
    y::Vector{Float64}
    sqrtw::W
    idrefs::Vector{Rid}
    timerefs::Vector{Rtime}
end

function FactorModel{W, Rid, Rtime}(y::Vector{Float64}, sqrtw::W, idrefs::Vector{Rid}, timerefs::Vector{Rtime}, rank::Int)
    FactorModel{rank, W, Rid, Rtime}(y, sqrtw, idrefs, timerefs)
end

rank{Rank}(::FactorModel{Rank}) = Rank

type FactorSolution{Rank, Tid, Ttime} <: AbstractFactorSolution{Rank}
    idpool::Tid
    timepool::Ttime
end

function FactorSolution{Tid, Ttime}(idpool::Tid, timepool::Ttime)
    r = size(idpool, 2)
    @assert r == size(timepool, 2)
    FactorSolution{r, Tid, Ttime}(idpool, timepool)
end

function view(f::AbstractFactorSolution, I::Union{AbstractArray,Colon,Int64}...)
    FactorSolution(view(f.idpool, I...), view(f.timepool, I...))
end

## subtract_factor! and subtract_b!
function subtract_factor!(fm::AbstractFactorModel, fs::AbstractFactorSolution)
    for r in 1:rank(fm)
        subtract_factor!(fm, view(fs, :, r))
    end
end

function subtract_factor!(fm::AbstractFactorModel, fs::FactorSolution{1})
    @inbounds @simd for i in 1:length(fm.y)
        fm.y[i] -= fm.sqrtw[i] * fs.idpool[fm.idrefs[i]] * fs.timepool[fm.timerefs[i]]
    end
end

## rescale a factor model
function reverse{R}(m::Matrix{R})
    out = similar(m)
    for j in 1:size(m, 2)
        invj = size(m, 2) + 1 - j
        @inbounds @simd for i in 1:size(m, 1)
            out[i, j] = m[i, invj]
        end
    end
    return out
end
function rescale!(fs::FactorSolution{1})
    out = norm(fs.timepool)
    scale!(fs.idpool, out)
    scale!(fs.timepool, 1/out)
end
# normalize factors and loadings so that F'F = Id, Lambda'Lambda diagonal
function rescale!(newfs::AbstractFactorSolution, fs::AbstractFactorSolution)
    U = eigfact!(Symmetric(At_mul_B(fs.timepool, fs.timepool)))
    sqrtDx = diagm(sqrt(abs(U[:values])))
    A_mul_B!(newfs.idpool,  fs.idpool,  U[:vectors] * sqrtDx)
    V = eigfact!(At_mul_B(newfs.idpool, newfs.idpool))
    A_mul_B!(newfs.idpool, fs.idpool, reverse(U[:vectors] * sqrtDx * V[:vectors]))
    A_mul_B!(newfs.timepool, fs.timepool, reverse(U[:vectors] * (sqrtDx \ V[:vectors])))
    return newfs
end

rescale(fs::FactorSolution) = rescale!(similar(fs), fs)

## Create dataframe from pooledfactors
function getfactors(fp::AbstractFactorModel, fs::AbstractFactorSolution)
    # partial out Y and X with respect to i.id x factors and i.time x loadings
    newfes = FixedEffect[]
    for r in 1:rank(fp)
        idinteraction = build_interaction(fp.timerefs, view(fs.timepool, :, r))
        idfe = FixedEffect(fp.idrefs, size(fs.idpool, 1), fp.sqrtw, idinteraction, :id, :time, :(idxtime))
        push!(newfes, idfe)
        timeinteraction = build_interaction(fp.idrefs, view(fs.idpool, :, r))
        timefe = FixedEffect(fp.timerefs, size(fs.timepool, 1), fp.sqrtw, timeinteraction, :time, :id, :(timexid))
        push!(newfes, timefe)
    end
    # obtain the residuals and cross
    return newfes
end

function build_interaction(refs::Vector, pool::AbstractVector)
    interaction = Array(Float64, length(refs))
    @inbounds @simd for i in 1:length(refs)
        interaction[i] = pool[refs[i]]
    end
    return interaction
end

function DataFrame(fp::AbstractFactorModel, fs::AbstractFactorSolution, esample::BitVector)
    df = DataFrame()
    anyNA = all(esample)
    for r in 1:rank(fp)
        # loadings
        df[convert(Symbol, "loadings$r")] = build_column(fp.idrefs, fs.idpool[:, r], esample)
        df[convert(Symbol, "factors$r")] = build_column(fp.timerefs, fs.timepool[:, r], esample)
    end
    return df
end
function build_column(refs::Vector, pool::Vector, esample::BitVector)
    T = eltype(refs)
    newrefs = fill(zero(T), length(esample))
    newrefs[esample] = refs
    return PooledDataArray(RefArray(newrefs), pool)
end

##############################################################################
##
## Interactive Fixed Effect Models
##
##############################################################################

type InteractiveFixedEffectsModel{Rank, W, Rid, Rtime} <: AbstractFactorModel{Rank}
    y::Vector{Float64}
    sqrtw::W
    X::Matrix{Float64}
    idrefs::Vector{Rid}
    timerefs::Vector{Rtime}
end

function InteractiveFixedEffectsModel{W, Rid, Rtime}(y::Vector{Float64}, sqrtw::W, X::Matrix{Float64}, idrefs::Vector{Rid}, timerefs::Vector{Rtime}, rank::Int)
    InteractiveFixedEffectsModel{rank, W, Rid, Rtime}(y, sqrtw, X, idrefs, timerefs)
end

rank{Rank}(::InteractiveFixedEffectsModel{Rank}) = Rank

convert{Rank, W, Rid, Rtime}(::Type{FactorModel}, f::InteractiveFixedEffectsModel{Rank, W, Rid, Rtime}) = FactorModel{Rank, W, Rid, Rtime}(f.y, f.sqrtw, f.idrefs, f.timerefs)

type InteractiveFixedEffectsSolution{Rank, Tb, Tid, Ttime} <: AbstractFactorSolution{Rank}
    b::Tb
    idpool::Tid
    timepool::Ttime
end
function InteractiveFixedEffectsSolution{Tb, Tid, Ttime}(b::Tb, idpool::Tid, timepool::Ttime)
    r = size(idpool, 2)
    r == size(timepool, 2) || throw("factors and loadings don't have same dimension")
    InteractiveFixedEffectsSolution{r, Tb, Tid, Ttime}(b, idpool, timepool)
end
convert(::Type{FactorSolution}, f::InteractiveFixedEffectsSolution) = FactorSolution(f.idpool, f.timepool)

type InteractiveFixedEffectsSolutionT{Rank, Tb, Tid, Ttime} <: AbstractFactorSolution{Rank}
    b::Tb
    idpool::Tid
    timepool::Ttime
end
function InteractiveFixedEffectsSolutionT{Tb, Tid, Ttime}(b::Tb, idpool::Tid, timepool::Ttime)
    r = size(idpool, 1)
    r == size(timepool, 1) || throw("factors and loadings don't have same dimension")
    InteractiveFixedEffectsSolutionT{r, Tb, Tid, Ttime}(b, idpool, timepool)
end

function rescale(fs::InteractiveFixedEffectsSolution)
    fss = FactorSolution(fs.idpool, fs.timepool)
    newfss = similar(fss)
    rescale!(newfss, fss)
    InteractiveFixedEffectsSolution(fs.b, newfss.idpool, newfss.timepool)
end

type HalfInteractiveFixedEffectsModel{Rank, W, Rid, Rtime} <: AbstractFactorModel
    y::Vector{Float64}
    sqrtw::W
    X::Matrix{Float64}
    idrefs::Vector{Rid}
    timerefs::Vector{Rtime}
    timepool::Matrix{Float64}
    size::Tuple{Int, Int}
end
function HalfInteractiveFixedEffectsModel{W, Rid, Rtime}(y::Vector{Float64}, sqrtw::W, X::Matrix{Float64}, idrefs::Vector{Rid}, timerefs::Vector{Rtime}, timepool::Matrix{Float64}, size, rank::Int)
    HalfInteractiveFixedEffectsModel{rank, W, Rid, Rtime}(y, sqrtw, X, idrefs, timerefs, timepool, size)
end

type HalfInteractiveFixedEffectsSolution{Tb, Tid} <: AbstractFactorSolution
    b::Tb
    idpool::Tid
end

##############################################################################
##
## Results
##
##############################################################################'

type FactorResult
    esample::BitVector
    augmentdf::DataFrame

    ess::Float64
    iterations::Int64
    converged::Bool
end

# result
type InteractiveFixedEffectsResult <: AbstractRegressionResult
    coef::Vector{Float64}   # Vector of coefficients
    vcov::Matrix{Float64}   # Covariance matrix

    esample::BitVector      # Is the row of the original dataframe part of the estimation sample?
    augmentdf::DataFrame

    coefnames::Vector       # Name of coefficients
    yname::Symbol           # Name of dependent variable
    formula::Formula        # Original formula

    nobs::Int64             # Number of observations
    df_residual::Int64      # degree of freedoms

    r2::Float64             # R squared
    r2_a::Float64           # R squared adjusted
    r2_within::Float64      # R within

    ess::Float64
    iterations::Int         # Number of iterations
    converged::Bool         # Has the demeaning algorithm converged?

end

predict(::InteractiveFixedEffectsResult, ::AbstractDataFrame) = error("predict is not defined for linear factor models. Use the option save = true")
residuals(::InteractiveFixedEffectsResult, ::AbstractDataFrame) = error("residuals is not defined for linear factor models. Use the option save = true")
title(::InteractiveFixedEffectsResult) = "Linear Factor Model"
top(x::InteractiveFixedEffectsResult) = [
            "Number of obs" sprint(showcompact, nobs(x));
            "Degree of freedom" sprint(showcompact, nobs(x) - df_residual(x));
            "R2"  @sprintf("%.3f", x.r2);
            "R2 within"  @sprintf("%.3f", x.r2_within);
            "Iterations" sprint(showcompact, x.iterations);
            "Converged" sprint(showcompact, x.converged)
            ]

ext: jl


text:
function voigt(x, σ, γ)
    faddeeva(z) = erfcx(-im * z)
    z = (x + im * γ) / σ / sqrt(2)
    return real(faddeeva(z)) / (σ * sqrt(2pi))
end

function voigt_fwhm(σ, γ)
    # J.J.Olivero and R.L. Longbothum in Empirical fits to the Voigt line width: A brief review, JQSRT 17, P233, 1977
    # http://snst-hu.lzu.edu.cn/zhangyi/ndata/Voigt_profile.html
    fg = σ * 2.355
    fl = 2 * γ
    0.5346 * fl + sqrt(0.2166 * fl^2 + fg^2)
    # fl * (0.5346 + sqrt(0.2166 + (fg/fl)^2))
end

#=
@gp "unset grid" xr=[-10,10]
x = -100:0.01:100
σ = 1.53;  γ = 0. ;  vp = voigt.(x, σ, γ);  @gp :- x vp "w l"; println(int_tabulated(x, vp)); print(voigt_fwhm(σ, γ) - QSFit.estimate_fwhm(x, vp))
σ = 1.30;  γ = 0.5;  vp = voigt.(x, σ, γ);  @gp :- x vp "w l"; println(int_tabulated(x, vp)); print(voigt_fwhm(σ, γ) - QSFit.estimate_fwhm(x, vp))
σ = 0.01;  γ = 1.8;  vp = voigt.(x, σ, γ);  @gp :- x vp "w l"; println(int_tabulated(x, vp)); print(voigt_fwhm(σ, γ) - QSFit.estimate_fwhm(x, vp))
σ = 1.  ;  γ = 1. ;  vp = voigt.(x, σ, γ);  @gp :- x vp "w l"; println(int_tabulated(x, vp)); print(voigt_fwhm(σ, γ) - QSFit.estimate_fwhm(x, vp))
σ = 1.  ;  γ = 5. ;  vp = voigt.(x, σ, γ);  @gp :- x vp "w l"; println(int_tabulated(x, vp)); print(voigt_fwhm(σ, γ) - QSFit.estimate_fwhm(x, vp))
=#

# ____________________________________________________________________
# SpecLineVoigt
#

mutable struct SpecLineVoigt <: AbstractComponent
    norm::Parameter
    center::Parameter
    fwhm::Parameter
    gamma::Parameter
    voff::Parameter
    norm_integrated::Bool

    function SpecLineVoigt(center::Number)
        out = new(Parameter(1),
                  Parameter(center),
                  Parameter(3000),
                  Parameter(3000),
                  Parameter(0),
                  true)

        @assert center > 0
        out.norm.low = 0
        out.center.low = 0
        out.fwhm.low = 0
        out.gamma.low = 0
        out.voff.low = 0
        out.center.fixed = true
        return out
    end
end

function prepare!(comp::SpecLineVoigt, domain::Domain{1})
    return fill(NaN, length(domain))
end

function evaluate!(buffer, comp::SpecLineVoigt, x::Domain{1},
                   norm, center, fwhm, gamma, voff)

    x0 = center - (voff / 3.e5) * center
    σ = fwhm  / 3.e5 * center / 2.355
    γ = gamma / 3.e5 * center / 2.
    buffer .= norm .* voigt.(x .- x0, σ, γ)
    comp.norm_integrated  ||  (buffer .*= (σ * sqrt(2pi)))
end

#=
    x = Domain(500:1:1500.)
    comp = QSFit.SpecLineVoigt(1000.)
    comp.fwhm.val  = 3e4
    comp.gamma.val = 3e4
    ceval = GFit.CompEval(comp, x)
    evaluate!(ceval)
    @gp x[:] ceval.buffer ./ maximum(ceval.buffer) "w l"
=#


text:
module Abaco

using Dates

const DEFAULT_TYPE = ""

include("logger.jl")
include("types.jl")
include("context.jl")
include("extractor.jl")
include("formula.jl")
include("node.jl")
include("poll.jl")
include("time.jl")

export abaco_init, formula, add_formulas, node,
       ingest, ingest!,
       last_value, last_point, get_values, get_collected, sum_collected,
       setup_settings, oncomplete, nowts,
       Context, Node, dependents, etype, span

"""
    abaco_init(onresult; handle=nothing, interval::Int=900, ages::Int=4, emitone=true)::Context

Initialize the abaco context:

* `onresult`: function callback that gets called each time a formula value is computed.

* `handle`: user defined object.
            If handle is defined it is the first argument of `onresult`, default to `nothing`.

* `interval`: the span interval in seconds, default to 900 seconds (15 minutes).
              If interval is equal to -1 there is just one infinite time span.

* `ages`: the number of active rops managed by the abaco. Default to 4.
          If `interval` is equal to -1 `ages`  is not applicable because it loses meaning.

* `emitone`: if `true` emits for each time span at most 1 formula value,
             otherwise emits a new result at every new inputs. Defaut to `true`

Example 1: defining `onresult` callback that uses of an handle object.

```
# the handle object is a socket
sock = connect(3001)

function onresult(handle, ts, ne, formula_name, value, inputs)
    # build a pkt message from ts, ne, ...
    pkt = ...
    write(sock, pkt)
end
```

Example 2: defining `onresult` callback that doesn't use an handle object.

```julia
abaco = abaco_init(onresult, handle=sock)

function onresult(ts, ne, formula_name, value, inputs)
    @info "[\$ts][\$ne] function \$fname=\$value"
end

abaco = abaco_init(onresult)
```
"""
function abaco_init(onresult;
                    handle=nothing,
                    interval::Int=-1,
                    ages::Int=1,
                    emitone::Bool=false)::Context
    DEBUG = get(ENV, "ABACO_DEBUG", "0")
    logging(debug = DEBUG=="0" ? [] : [@__MODULE__, Main])
    @debug "resetting abaco ..."
    cfg = SnapsSetting(handle, emitone, onresult)
    Context(interval, ages, Dict(DEFAULT_TYPE=>cfg), onresult)
end

abaco_init() = abaco_init() do ts, ne, name, value, inputs
    println("$name(ts:$ts, ne:$ne) = $value")
end

function oncomplete(onresult, abaco::Context)
    abaco.cfg[DEFAULT_TYPE].oncomplete = onresult
end

function abaco_multi()
end

function loginit()
    DEBUG = get(ENV, "ABACO_DEBUG", "0")
    Abaco.logging(debug = DEBUG=="0" ? [] : [Abaco, Main])
end

end


text:
using Test

import MIDI
import SampledSignals

using Revise
using CureMIDI

@testset "basic synth" begin

    # pitch, velocity, position, duration,channel
    C = MIDI.Note(60, 96, 0, 3000)
    notes = MIDI.Notes()
    push!(notes, C)

    track = MIDI.MIDITrack()
    MIDI.addnotes!(track, notes)

    tpq = Int16(960)
    bpm = 60
    sample_rate = 44100
    sampled_buf = synth(track, tpq, bpm, sample_rate, "violin_sample.sf2")

    # import PortAudio
    # stream = PortAudio.PortAudioStream(0, 2)
    # PortAudio.write(stream, sampled_buf)
end

@testset "pitch range" begin
    test_note = [MIDI.Note(0, 127, 0, 100), MIDI.Note(127, 127, 0, 100)]
    for note in test_note
        notes = MIDI.Notes()
        push!(notes, note)

        track = MIDI.MIDITrack()
        MIDI.addnotes!(track, notes)
        tpq = Int16(960)
        bpm = 60
        sample_rate = 44100
        sampled_buf = synth(track, tpq, bpm, sample_rate, "violin_sample.sf2")
        #import PortAudio
        #stream = PortAudio.PortAudioStream(0, 2)
        #PortAudio.write(stream, sampled_buf)
    end
end


text:
@testset "nested mock call" begin
    readfile(filename) = (@mock isfile(filename)) ? read((@mock open(filename)), String) : ""

    # Testing with both generic and anonymous functions
    patches = Patch[
        @patch isfile(f::AbstractString) = f == "foo"
        @patch open(f::AbstractString) = IOBuffer("bar")
    ]

    @test readfile("foo") == ""

    apply(patches) do
        @test readfile("foo") == "bar"
    end

    @test readfile("foo") == ""
end


text:
# ImagePFM.jl
# FileIO.jl integration for PFM image files.
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

using Test

using ImagePFM
import ImagePFM: little_endian, _read_line, _read_type, _parse_endianness, _parse_int, _parse_img_size, _TypeStream, _read_matrix
import ColorTypes: RGB
import FileIO: @format_str

# testset variables
endian_f        = little_endian ? ltoh : ntoh
test_matrix     = RGB{Float32}[RGB(1.0f1, 2.0f1, 3.0f1) RGB(1.0f2, 2.0f2, 3.0f2)
                               RGB(4.0f1, 5.0f1, 6.0f1) RGB(4.0f2, 5.0f2, 6.0f2)
                               RGB(7.0f1, 8.0f1, 9.0f1) RGB(7.0f2, 8.0f2, 9.0f2)]
test_matrix_err = RGB{Float64}[RGB(1.0e1, 2.0e1, 3.0e1) RGB(1.0e2, 2.0e2, 3.0e2)
                               RGB(4.0e1, 5.0e1, 6.0e1) RGB(4.0e2, 5.0e2, 6.0e2)
                               RGB(7.0e1, 8.0e1, 9.0e1) RGB(7.0e2, 8.0e2, 9.0e2)]
expected_output_file = little_endian ? "reference_le.pfm" : "reference_be.pfm"

# test color write to IO
@testset "write" begin
    image = test_matrix
    io = IOBuffer()
    ImagePFM.write(io, format"PFM", image)
    @test take!(io) == open(read, expected_output_file)
end

# test _parse_endianness
@testset "_parse_endianness" begin
    @test _parse_endianness("1.0") == ntoh
    @test _parse_endianness("-1.0") == ltoh
    @test_throws InvalidPfmFileFormat _parse_endianness("abba")
    @test_throws InvalidPfmFileFormat _parse_endianness("2.0")
end

# test _parse_int
@testset "_parse_int" begin
    @test _parse_int("12") === UInt(12)
    @test_throws InvalidPfmFileFormat _parse_int("abba")
    @test_throws InvalidPfmFileFormat _parse_int("-1")
    @test_throws InvalidPfmFileFormat _parse_int("1.0")
end

# test _parse_img_size
@testset "_parse_img_size" begin
    @test _parse_img_size("1920 1080") == UInt[1920, 1080]

    # test exceptions
    @test_throws InvalidPfmFileFormat _parse_img_size("abba 1920")
    @test_throws InvalidPfmFileFormat _parse_img_size("1920 -1080")
    @test_throws InvalidPfmFileFormat _parse_img_size("1920 1080 256")
    @test_throws InvalidPfmFileFormat _parse_img_size("1920")
end

# test _read_line
@testset "_read_line" begin
    io = IOBuffer(b"hello\nworld")
    @test _read_line(io) == "hello\n"
    @test _read_line(io) == "world"
    @test _read_line(io) === nothing

    # test exceptions
    io = IOBuffer(b"è")
    @test_throws InvalidPfmFileFormat _read_line(io)
end

# test _read_matrix
@testset "_read_matrix" begin
    io = IOBuffer()
    write(io, test_matrix)
    seekstart(io)
    @test all(_read_matrix(io, eltype(test_matrix), size(test_matrix)...) .== test_matrix)

    # test exceptions
    io = IOBuffer()
    write(io, test_matrix[begin:end-1])
    seekstart(io)
    @test_throws EOFError _read_matrix(io, eltype(test_matrix), size(test_matrix)...)
end

# test read(io, ::FE"PFM")
@testset "read" begin
    img = ImagePFM.read(IOBuffer(open(read, expected_output_file)), format"PFM")
    @test size(img) == size(test_matrix)
    @test all(img == test_matrix)

    # test exceptions
    @test_throws InvalidPfmFileFormat ImagePFM.read(IOBuffer(b"PF\n3 2\n-1.0\nstop"), format"PFM")
end

# test write/read compatibility
@testset "write/read compatibility" begin
    img = test_matrix
    io = IOBuffer()
    ImagePFM.write(io, format"PFM", img)
    seekstart(io)
    @test all(ImagePFM.read(io, format"PFM") == img)
end

# test throw RGB eltype error
@testset "throw RGB eltype error" begin
    img = test_matrix_err
    io = IOBuffer()
    @test_throws InvalidRGBEltype ImagePFM.write(io, format"PFM", img)
end


text:
@testset "Test ForecastCache" begin
    sys = IS.SystemData()
    name = "Component1"
    component = IS.TestComponent(name, 5)
    IS.add_component!(sys, component)
    file = joinpath(FORECASTS_DIR, "DateTimeAsColumnDeterministic.csv")
    forecast = IS.Deterministic("test", file, component, Dates.Hour(1))
    initial_timestamp = IS.get_initial_timestamp(forecast)
    initial_times = collect(IS.get_initial_times(forecast))
    interval = IS.get_interval(forecast)
    IS.add_time_series!(sys, component, forecast)

    cache = IS.ForecastCache(IS.Deterministic, component, "test")
    @test cache.in_memory_count == 168
    @test IS.get_next_time(cache) == initial_timestamp
    @test length(cache) == cache.common.num_iterations == 168

    # Iterate over all initial times with default cache size.
    cache = IS.ForecastCache(IS.Deterministic, component, "test")
    for (i, ta) in enumerate(cache)
        it = initial_times[i]
        @test TimeSeries.timestamp(ta) ==
              IS.get_time_series_timestamps(component, forecast, it)
        @test TimeSeries.values(ta) == IS.get_time_series_values(component, forecast, it)
    end

    IS.reset!(cache)
    for it in initial_times
        ta = IS.get_next_time_series_array!(cache)
        @test first(TimeSeries.timestamp(ta)) == it
        @test TimeSeries.timestamp(ta) ==
              IS.get_time_series_timestamps(component, forecast, it)
        @test TimeSeries.values(ta) == IS.get_time_series_values(component, forecast, it)
    end
    @test IS.get_next_time(cache) === nothing

    # Iterate over all initial times with custom cache size.
    cache = IS.ForecastCache(IS.Deterministic, component, "test"; cache_size_bytes=1024)
    @test length(cache) == cache.common.num_iterations == 168
    for (i, ta) in enumerate(cache)
        it = initial_times[i]
        @test TimeSeries.timestamp(ta) ==
              IS.get_time_series_timestamps(component, forecast, it)
        @test TimeSeries.values(ta) == IS.get_time_series_values(component, forecast, it)
    end

    IS.reset!(cache)
    for it in initial_times
        ta = IS.get_next_time_series_array!(cache)
        @test TimeSeries.timestamp(ta) ==
              IS.get_time_series_timestamps(component, forecast, it)
        @test TimeSeries.values(ta) == IS.get_time_series_values(component, forecast, it)
    end

    # Start at an offset.
    cache = IS.ForecastCache(
        IS.Deterministic,
        component,
        "test";
        start_time=Dates.DateTime("2020-01-02T00:00:00"),
    )
    for (i, ta) in enumerate(cache)
        it = initial_times[i + 24]
        @test TimeSeries.timestamp(ta) ==
              IS.get_time_series_timestamps(component, forecast, it)
        @test TimeSeries.values(ta) == IS.get_time_series_values(component, forecast, it)
    end

    # Test caching internals.
    cache = IS.ForecastCache(IS.Deterministic, component, "test"; cache_size_bytes=1024)
    @test cache.in_memory_count == 5
    @test IS.get_next_time(cache) == initial_timestamp
    for it in initial_times[1:(cache.in_memory_count)]
        ta = IS.get_next_time_series_array!(cache)
        @test IS._get_last_cached_time(cache) == initial_times[5]
        @test TimeSeries.timestamp(ta) ==
              IS.get_time_series_timestamps(component, forecast, it)
        @test TimeSeries.values(ta) == IS.get_time_series_values(component, forecast, it)
    end

    # The next access should trigger a read.
    ta = IS.get_next_time_series_array!(cache)
    @test IS._get_last_cached_time(cache) == initial_times[10]
    @test TimeSeries.timestamp(ta) ==
          IS.get_time_series_timestamps(component, forecast, initial_times[6])
    @test TimeSeries.values(ta) ==
          IS.get_time_series_values(component, forecast, initial_times[6])
    @test IS.get_next_time(cache) == initial_times[7]
end

@testset "Test StaticTimeSeriesCache" begin
    sys = IS.SystemData()
    name = "Component1"
    component = IS.TestComponent(name, 5)
    IS.add_component!(sys, component)

    initial_timestamp = Dates.DateTime("2020-09-01")
    resolution = Dates.Hour(1)

    len = 96
    data = TimeSeries.TimeArray(
        range(initial_timestamp; length=len, step=resolution),
        rand(len),
    )
    ts = IS.SingleTimeSeries(data=data, name="test")
    IS.add_time_series!(sys, component, ts)

    cache = IS.StaticTimeSeriesCache(IS.SingleTimeSeries, component, "test")
    @test cache.in_memory_rows == len
    @test IS.get_next_time(cache) == initial_timestamp
    @test length(cache) == cache.common.num_iterations == len

    # Iterate over all initial times with default cache size.
    cache = IS.StaticTimeSeriesCache(IS.SingleTimeSeries, component, "test")
    for (i, ta) in enumerate(cache)
        it = initial_timestamp + (i - 1) * resolution
        @test TimeSeries.timestamp(ta) ==
              IS.get_time_series_timestamps(component, ts, it, len=1)
        @test TimeSeries.values(ta) == IS.get_time_series_values(component, ts, it, len=1)
    end

    ta = IS.get_next_time_series_array!(cache)
    @test first(TimeSeries.timestamp(ta)) == initial_timestamp
    @test TimeSeries.timestamp(ta) ==
          IS.get_time_series_timestamps(component, ts, initial_timestamp, len=1)
    @test TimeSeries.values(ta) ==
          IS.get_time_series_values(component, ts, initial_timestamp, len=1)

    # Iterate over all initial times with custom cache size.
    cache_size_bytes = 96
    cache = IS.StaticTimeSeriesCache(
        IS.SingleTimeSeries,
        component,
        "test";
        cache_size_bytes=cache_size_bytes,
    )
    @test cache.in_memory_rows == cache_size_bytes / 8
    @test length(cache) == cache.common.num_iterations == len
    ta = IS.get_next_time_series_array!(cache)
    @test first(TimeSeries.timestamp(ta)) == initial_timestamp
    @test length(ta) == 1
    @test TimeSeries.values(ta) == [TimeSeries.values(data)[1]]

    IS.reset!(cache)
    for (i, ta) in enumerate(cache)
        it = initial_timestamp + (i - 1) * resolution
        @test TimeSeries.timestamp(ta) ==
              IS.get_time_series_timestamps(component, ts, it; len=1)
        @test TimeSeries.values(ta) == IS.get_time_series_values(component, ts, it; len=1)
    end

    IS.reset!(cache)
    for i in 1:3
        ta = IS.get_next_time_series_array!(cache)
        it = initial_timestamp + (i - 1) * resolution
        @test TimeSeries.timestamp(ta) ==
              IS.get_time_series_timestamps(component, ts, it; len=1)
        @test TimeSeries.values(ta) == IS.get_time_series_values(component, ts, it; len=1)
    end

    cache_size_bytes = 96
    start_time = Dates.DateTime("2020-09-02T08:00:00")
    cache = IS.StaticTimeSeriesCache(
        IS.SingleTimeSeries,
        component,
        "test";
        start_time=start_time,
        cache_size_bytes=cache_size_bytes,
    )
    @test cache.in_memory_rows == cache_size_bytes / 8
    @test cache.common.num_iterations ==
          len - (start_time - initial_timestamp) / Dates.Millisecond(resolution)
    @test IS._get_length_remaining(cache) == cache.common.num_iterations
    for i in 1:2
        ta = IS.get_next_time_series_array!(cache)
        it = start_time + (i - 1) * resolution
        @test TimeSeries.timestamp(ta) ==
              IS.get_time_series_timestamps(component, ts, it; len=1)
        @test TimeSeries.values(ta) == IS.get_time_series_values(component, ts, it; len=1)
    end
end

@testset "Test DeterministicSingleTimeSeries with ForecastCache" begin
    sys = IS.SystemData()
    component = IS.TestComponent("Component1", 5)
    IS.add_component!(sys, component)

    resolution = Dates.Minute(5)
    dates = create_dates("2020-01-01T00:00:00", resolution, "2020-01-31T23:00:00")
    data = collect(1:length(dates))
    ta = TimeSeries.TimeArray(dates, data, [IS.get_name(component)])
    name = "val"
    ts = IS.SingleTimeSeries(name, ta)
    IS.add_time_series!(sys, component, ts)
    horizon = 24
    interval = Dates.Hour(1)
    IS.transform_single_time_series!(
        sys,
        IS.DeterministicSingleTimeSeries,
        horizon,
        interval,
    )

    forecast = IS.get_time_series(IS.AbstractDeterministic, component, name)
    initial_times = collect(IS.get_initial_times(forecast))
    cache =
        IS.ForecastCache(IS.AbstractDeterministic, component, name; cache_size_bytes=1024)

    for (i, ta) in enumerate(cache)
        @test TimeSeries.timestamp(ta) ==
              IS.get_time_series_timestamps(component, forecast, initial_times[i])
        @test TimeSeries.values(ta) ==
              IS.get_time_series_values(component, forecast, initial_times[i])
    end
end

@testset "Test Probabilistic with ForecastCache" begin
    initial_time = Dates.DateTime("2020-09-01")
    interval = Dates.Hour(1)
    resolution = Dates.Hour(1)
    name = "test"
    horizon = 24
    data = SortedDict{Dates.DateTime, Matrix{Float64}}()
    for (i, it) in enumerate(range(initial_time, step=interval, length=100))
        data[it] = ones(horizon, 99) * i
    end
    sys = IS.SystemData()
    component_name = "Component1"
    component = IS.TestComponent(component_name, 5)
    IS.add_component!(sys, component)
    forecast = IS.Probabilistic(name, data, ones(99), resolution)
    IS.add_time_series!(sys, component, forecast)

    # Iterate over all initial times with custom cache size.
    sz = 1024 * 1024
    cache = IS.ForecastCache(IS.Probabilistic, component, "test"; cache_size_bytes=sz)
    initial_times = collect(keys(data))
    @test cache.in_memory_count == trunc(Int, sz / (99 * 8 * 24))
    for (i, ta) in enumerate(cache)
        it = initial_times[i]
        @test TimeSeries.timestamp(ta) ==
              IS.get_time_series_timestamps(component, forecast, it)
        @test TimeSeries.values(ta) == IS.get_time_series_values(component, forecast, it)
    end
end


text:
"""
Sensor model with gaussian noise in position and velocity measurement
"""
struct SimpleSensor
    pos_noise::Float64
    vel_noise::Float64
end

@with_kw mutable struct POMDPSensor
    pomdp::SingleOCPOMDP = SingleOCPOMDP()
    sensor::SimpleSensor = SimpleSensor(0., 0.)
end

const OFF_KEY = -1

"""
    measure(scene::Scene, sensor::POMDPSensor, env::CrosswalkEnv)
return a measurement of a given scene from the point of view of the ego car
"""
function measure(scene::Scene, sensor::POMDPSensor, env::CrosswalkEnv)
    pos_noise = sensor.sensor.pos_noise
    vel_noise = sensor.sensor.vel_noise
    pomdp = sensor.pomdp
    o = Dict{Int64, SingleOCObs}()
    ego = scene[findfirst(EGO_ID, scene)].state
    for i=1:scene.n
        veh = scene[i]
        if veh.id == 1
            continue
        end
        car = veh.state
        if is_observable_fixed(ego, car, env)
            y_pos = min(pomdp.y_goal, car.posG.y + pos_noise*randn())
            obs_state = VehicleState(VecSE2(car.posG.x + pos_noise*randn(),
                                                      y_pos,
                                                      car.posG.θ),
                                                      env.roadway,
                                                      car.v + vel_noise*randn())
            collision = collision_checker(ego, veh.state, pomdp.ego_type, veh.def)
            o[veh.id] = SingleOCState(collision, ego, veh.state)
        end
    end
    o[OFF_KEY] = SingleOCState(false, ego, get_off_the_grid(pomdp))
    return o
end

"""
    measure(scene::Scene, sensor::SimpleSensor, env::CrosswalkEnv)
"""
function measure(scene::Scene, sensor::SimpleSensor, env::CrosswalkEnv)
    pos_noise = sensor.pos_noise
    vel_noise = sensor.vel_noise
    o = Vehicle[]
    ego = scene[findfirst(EGO_ID, scene)]
    push!(o, ego)
    for veh in scene
        if veh.id == ego.id
            continue
        end
        car = veh.state
        if is_observable_fixed(ego.state, car, env)
            obs_state = VehicleState(VecSE2(car.posG.x + pos_noise*randn(),
                                            car.posG.y + pos_noise*randn(),
                                            car.posG.θ),
                                     env.roadway,
                                     car.v + vel_noise*randn())
            push!(o, Vehicle(obs_state, VehicleDef(), veh.id))
        end
    end
    return o
end

"""
    measure(ego::VehicleState, cars::EntityFrame{VehicleState, VehicleDef, Int64}, model::SimpleSensor, env::CrosswalkEnv)
return a measurement of a given scene from the point of view of the ego car
"""
function measure(ego::Vehicle, cars::Vector{Entity{VehicleState, VehicleDef, Int64}}, model::SimpleSensor, env::CrosswalkEnv)
    observed = Vehicle[] # sizehint!(Vector{Vehicle(0)}, 5)
    for veh in cars
        car = veh.state
        if is_observable_fixed(ego.state, car, env)
            obs_state = VehicleState(VecSE2(car.posG.x + model.pos_noise*randn(),
                                                      car.posG.y + model.pos_noise*randn(),
                                                      car.posG.θ),
                                                      env.roadway,
                                                      car.v + model.vel_noise*randn())
            push!(observed, Vehicle(obs_state, veh.def, veh.id))
        end
    end
    insert!(observed, 1, ego)
    return observed
end


text:
#    Modified Bessel functions of the second kind of order zero and one
#                       besselk0, besselk1
#
#    Scaled modified Bessel functions of the second kind of order zero and one
#                       besselk0x, besselk1x
#
#    (Scaled) Modified Bessel functions of the second kind of order nu
#                       besselk, besselkx
#
#
#    Calculation of besselk0 is done in two branches using polynomial approximations [1]
#
#    Branch 1: x < 1.0
#              besselk0(x) + log(x)*besseli0(x) = P7(x^2)
#                            besseli0(x) = [x/2]^2 * P6([x/2]^2) + 1
#    Branch 2: x >= 1.0
#              sqrt(x) * exp(x) * besselk0(x) = P22(1/x) / Q2(1/x)
#    where P7, P6, P22, and Q2 are 7, 6, 22, and 2 degree polynomials respectively.
#
#
#
#    Calculation of besselk1 is done in two branches using polynomial approximations [2]
#
#    Branch 1: x < 1.0
#              besselk1(x) - log(x)*besseli1(x) - 1/x = x*P8(x^2)
#                            besseli1(x) = [x/2]^2 * (1 + 0.5 * (*x/2)^2 + (x/2)^4 * P5([x/2]^2))
#    Branch 2: x >= 1.0
#              sqrt(x) * exp(x) * besselk1(x) = P22(1/x) / Q2(1/x)
#    where P8, P5, P22, and Q2 are 8, 5, 22, and 2 degree polynomials respectively.
#
#
#    The polynomial coefficients are taken from boost math library [3].
#    Evaluation of the coefficients using Remez.jl is prohibitive due to the increase
#    precision required in ArbNumerics.jl.
#
#    Horner's scheme is then used to evaluate all polynomials.
#
#    Calculation of besselk and besselkx can be done with recursion starting with
#    besselk0 and besselk1 and using upward recursion for any value of nu (order).
#
#                    K_{nu+1} = (2*nu/x)*K_{nu} + K_{nu-1}
#
#    When nu is large, a large amount of recurrence is necesary.
#    We consider uniform asymptotic expansion for large orders to more efficiently
#    compute besselk(nu, x) when nu is larger than 100 (Let's double check this cutoff)
#    The boundary should be carefully determined for accuracy and machine roundoff.
#    We use 10.41.4 from the Digital Library of Math Functions [5].
#    This is also 9.7.8 in Abramowitz and Stegun [6].
#    K_{nu}(nu*z) = sqrt(pi / 2nu) *exp(-nu*n)/(1+z^2)^1/4 * sum((-1^k)U_k(p) /nu^k)) for k=0 -> infty
#    The U polynomials are the most tricky. They are listed up to order 4 in Table 9.39
#    of [6]. For Float32, >=4 U polynomials are usually necessary. For Float64 values,
#    >= 8 orders are needed. However, this largely depends on the cutoff of order you need.
#    For moderatelly sized orders (nu=50), might need 11-12 orders to reach good enough accuracy
#    in double precision.
#
#    However, calculation of these higher order U polynomials are tedious. These have been hand
#    calculated and somewhat crosschecked with symbolic math. There could be errors. They are listed
#    in src/U_polynomials.jl as a reference as higher orders are impossible to find while being needed for any meaningfully accurate calculation.
#    For large orders these formulas will converge much faster than using upward recurrence.
#
#
# [1] "Rational Approximations for the Modified Bessel Function of the Second Kind
#     - K0(x) for Computations with Double Precision" by Pavel Holoborodko
# [2] "Rational Approximations for the Modified Bessel Function of the Second Kind
#     - K1(x) for Computations with Double Precision" by Pavel Holoborodko
# [3] https://github.com/boostorg/math/tree/develop/include/boost/math/special_functions/detail
# [4] "Computation of Bessel Functions of Complex Argument and Large Order" by Donald E. Amos
#      Sandia National Laboratories
# [5] https://dlmf.nist.gov/10.41
# [6] Abramowitz, Milton, and Irene A. Stegun, eds. Handbook of mathematical functions with formulas, graphs, and mathematical tables.
#     Vol. 55. US Government printing office, 1964.
#

"""
    besselk0(x::T) where T <: Union{Float32, Float64}

Modified Bessel function of the second kind of order zero, ``K_0(x)``.
"""
function besselk0(x::T) where T <: Union{Float32, Float64}
    x <= zero(T) && return throw(DomainError(x, "`x` must be nonnegative."))
    if x <= one(T)
        a = x * x / 4
        s = muladd(evalpoly(a, P1_k0(T)), inv(evalpoly(a, Q1_k0(T))), T(Y_k0))
        a = muladd(s, a, 1)
        return muladd(-a, log(x), evalpoly(x * x, P2_k0(T)))
    else
        s = exp(-x / 2)
        a = muladd(evalpoly(inv(x), P3_k0(T)), inv(evalpoly(inv(x), Q3_k0(T))), one(T)) * s / sqrt(x)
        return a * s
    end
end
"""
    besselk0x(x::T) where T <: Union{Float32, Float64}

Scaled modified Bessel function of the second kind of order zero, ``K_0(x)*e^{x}``.
"""
function besselk0x(x::T) where T <: Union{Float32, Float64}
    x <= zero(T) && return throw(DomainError(x, "`x` must be nonnegative."))
    if x <= one(T)
        a = x * x / 4
        s = muladd(evalpoly(a, P1_k0(T)), inv(evalpoly(a, Q1_k0(T))), T(Y_k0))
        a = muladd(s, a, 1)
        return muladd(-a, log(x), evalpoly(x * x, P2_k0(T))) * exp(x)
    else
        return muladd(evalpoly(inv(x), P3_k0(T)), inv(evalpoly(inv(x), Q3_k0(T))), one(T)) / sqrt(x)
    end
end
"""
    besselk1(x::T) where T <: Union{Float32, Float64}

Modified Bessel function of the second kind of order one, ``K_1(x)``.
"""
function besselk1(x::T) where T <: Union{Float32, Float64}
    x <= zero(T) && return throw(DomainError(x, "`x` must be nonnegative."))
    if x <= one(T)
        z = x * x
        a = z / 4
        pq = muladd(evalpoly(a, P1_k1(T)), inv(evalpoly(a, Q1_k1(T))), Y_k1(T))
        pq = muladd(pq * a, a, (a / 2 + 1))
        a = pq * x / 2
        pq = muladd(evalpoly(z, P2_k1(T)) / evalpoly(z, Q2_k1(T)), x, inv(x))
        return muladd(a, log(x), pq)
    else
        s = exp(-x / 2)
        a = muladd(evalpoly(inv(x), P3_k1(T)), inv(evalpoly(inv(x), Q3_k1(T))), Y2_k1(T)) * s / sqrt(x)
        return a * s
    end
end
"""
    besselk1x(x::T) where T <: Union{Float32, Float64}

Scaled modified Bessel function of the second kind of order one, ``K_1(x)*e^{x}``.
"""
function besselk1x(x::T) where T <: Union{Float32, Float64}
    x <= zero(T) && return throw(DomainError(x, "`x` must be nonnegative."))
    if x <= one(T)
        z = x * x
        a = z / 4
        pq = muladd(evalpoly(a, P1_k1(T)), inv(evalpoly(a, Q1_k1(T))), Y_k1(T))
        pq = muladd(pq * a, a, (a / 2 + 1))
        a = pq * x / 2
        pq = muladd(evalpoly(z, P2_k1(T)) / evalpoly(z, Q2_k1(T)), x, inv(x))
        return muladd(a, log(x), pq) * exp(x)
    else
        return muladd(evalpoly(inv(x), P3_k1(T)), inv(evalpoly(inv(x), Q3_k1(T))), Y2_k1(T)) / sqrt(x)
    end
end
#=
If besselk0(x) or besselk1(0) is equal to zero
this will underflow and always return zero even if besselk(x, nu)
is larger than the smallest representing floating point value.
In other words, for large values of x and small to moderate values of nu,
this routine will underflow to zero.
For small to medium values of x and large values of nu this will overflow and return Inf.
=#
"""
    besselk(x::T) where T <: Union{Float32, Float64}

Modified Bessel function of the second kind of order nu, ``K_{nu}(x)``.
"""
function besselk(nu, x::T) where T <: Union{Float32, Float64, BigFloat}
    T == Float32 ? branch = 20 : branch = 50
    if nu < branch
        return up_recurrence(x, besselk0(x), besselk1(x), nu)[1]
    else
        return besselk_large_orders(nu, x)
    end
end

"""
    besselk(x::T) where T <: Union{Float32, Float64}

Scaled modified Bessel function of the second kind of order nu, ``K_{nu}(x)*e^{x}``.
"""
function besselkx(nu::Int, x::T) where T <: Union{Float32, Float64}
    T == Float32 ? branch = 20 : branch = 50
    if nu < branch
        return up_recurrence(x, besselk0x(x), besselk1x(x), nu)[1]
    else
        return besselk_large_orders_scaled(nu, x)
    end
end
function besselk_large_orders(v, x::T) where T <: Union{Float32, Float64, BigFloat}
    S = promote_type(T, Float64)
    x = S(x)
    z = x / v
    zs = hypot(1, z)
    n = zs + log(z) - log1p(zs)
    coef = SQPIO2(S) * sqrt(inv(v)) * exp(-v*n) / sqrt(zs)
    p = inv(zs)
    p2  = v^2/fma(max(v,x), max(v,x), min(v,x)^2)

    return T(coef*Uk_poly_Kn(p, v, p2, T))
end
function besselk_large_orders_scaled(v, x::T) where T <: Union{Float32, Float64, BigFloat}
    S = promote_type(T, Float64)
    x = S(x)
    z = x / v
    zs = hypot(1, z)
    n = zs + log(z) - log1p(zs)
    coef = SQPIO2(S) * sqrt(inv(v)) * exp(-v*n + x) / sqrt(zs)
    p = inv(zs)
    p2  = v^2/fma(max(v,x), max(v,x), min(v,x)^2)

    return T(coef*Uk_poly_Kn(p, v, p2, T))
end


text:
"""
    Transform

Abstract type defining a transformation of the input.
"""
abstract type Transform end

Base.map(t::Transform, x::AbstractVector) = _map(t, x)
_map(t::Transform, x::AbstractVector) = t.(x)

"""
    IdentityTransform()

Transformation that returns exactly the input.
"""
struct IdentityTransform <: Transform end

(t::IdentityTransform)(x) = x
_map(::IdentityTransform, x::AbstractVector) = x

### TODO Maybe defining adjoints could help but so far it's not working

# @adjoint function ScaleTransform(s::T) where {T<:Real}
#     @check_args(ScaleTransform, s, s > zero(T), "s > 0")
#     ScaleTransform{T}(s),Δ->ScaleTransform{T}(Δ)
# end
#
# @adjoint function ScaleTransform(s::A) where {A<:AbstractVector{<:Real}}
#     @check_args(ScaleTransform, s, all(s.>zero(eltype(A))), "s > 0")
#     ScaleTransform{A}(s),Δ->begin; @show Δ,size(Δ); ScaleTransform{A}(Δ); end
# end

# @adjoint transform(t::ScaleTransform{<:AbstractVector{<:Real}},x::AbstractVector{<:Real}) = transform(t,x),Δ->(ScaleTransform(nothing),t.s.*Δ)
#
#     @adjoint transform(t::ARDTransform{<:Real},X::AbstractMatrix{<:Real},obsdim::Int) = transform(t,X,obsdim),Δ->begin
#     @show Δ,size(Δ);
#     return (obsdim == 1 ? ARD()Δ'.*X : ScaleTransform()Δ.*X,transform(t,Δ,obsdim),nothing)
#     end
#
# @adjoint transform(t::ScaleTransform{T},x::AbstractVecOrMat,obsdim::Int) where {T<:Real} = transform(t,x), Δ->(ScaleTransform(one(T)),t.s.*Δ,nothing)


text:
using Fst

lexicon = read_wfst(readall("fsts/paper_wfsts/lexicon.txt"))
create_pdf(lexicon, "fsts/paper_wfsts/lexicon.pdf")

lm = read_wfst(readall("fsts/paper_wfsts/lm.txt"))
create_pdf(lm, "fsts/paper_wfsts/lm.pdf")

lexicon_lm = compose_epsilon(lexicon, lm)
create_pdf(lexicon_lm, "fsts/paper_wfsts/lexicon_lm.pdf")

tm = read_wfst(readall("fsts/paper_wfsts/tm.txt"))
create_pdf(tm, "fsts/paper_wfsts/tm.pdf")

lexicon_tm = compose_epsilon(lexicon, tm)
lexicon_tm_trim = trim(lexicon_tm)
create_pdf(lexicon_tm_trim, "fsts/paper_wfsts/lexicon_tm_trim.pdf")

lexicon_tm_noeps = compose(lexicon, tm)
lexicon_tm_noeps_trim = trim(lexicon_tm_noeps)
create_pdf(lexicon_tm_noeps_trim, "fsts/paper_wfsts/lexicon_tm_noeps_trim.pdf")


text:
using ApproximateBayesianComputing
const ABC = ApproximateBayesianComputing
using Distributions

if VERSION >= v"0.7"
  using Random
  Random.seed!(1234)
  using Statistics
  using Distributed
  using Test
  #import Statistics: mean, median, maximum, minimum, quantile, std, var, cov, cor
else
  using Base.Test
  srand(1234)
end

include("test1.jl")

@time @test test1()


text:
function FFTOp(in,adj;normalize=true)
	norm = normalize ? sqrt(length(in[:])) : 1.
	if (adj)
		out = fft(in)/norm
	else
		out = bfft(in)/norm
	end

	return out
end


text:
"""
Module: DD (Decision Diagram)
"""

module DD

include("_bdd.jl")
# include("_bdd2.jl")

end


text:
using Base.Test, JuLIP

if hasase
   println("Testing `jbulk == bulk` ...")
   for sym in [:Fe, :W, :Al, :Cu, :Si, :C, :Mg, :Be]
      @test isapprox(bulk(sym), Atoms(ASE.bulk(string(sym))), tol = 1e-12)
      @test isapprox(bulk(sym, cubic=true), Atoms(ASE.bulk(string(sym), cubic=true)), tol = 1e-12)
   end
end


text:
# Based on atomics_base.jl which was copied from base/atomics.jl in julia

for typ in atomictypes
    lt = llvmtypes[typ]
    ilt = llvmtypes[inttype(typ)]
    rt = "$lt, $lt*"
    irt = "$ilt, $ilt*"
    @eval Base.getindex(ref::AtomicRef{$typ}) = let __x = ref.x; GC.@preserve __x begin
        llvmcall($"""
                 %ptr = inttoptr i$WORD_SIZE %0 to $lt*
                 %rv = load atomic $rt %ptr acquire, align $(gc_alignment(typ))
                 ret $lt %rv
                 """, $typ, Tuple{Ptr{$typ}}, ref.pointer)
    end end
    @eval Base.setindex!(ref::AtomicRef{$typ}, v::$typ) = let __x = ref.x; GC.@preserve __x begin
        llvmcall($"""
                 %ptr = inttoptr i$WORD_SIZE %0 to $lt*
                 store atomic $lt %1, $lt* %ptr release, align $(gc_alignment(typ))
                 ret void
                 """, Cvoid, Tuple{Ptr{$typ}, $typ}, ref.pointer, v)
    end end
    # Note: atomic_cas! succeeded (i.e. it stored "new") if and only if the result is "cmp"
    if typ <: Integer
        @eval Threads.atomic_cas!(ref::AtomicRef{$typ}, cmp::$typ, new::$typ) = let __x = ref.x; GC.@preserve __x begin
            llvmcall($"""
                     %ptr = inttoptr i$WORD_SIZE %0 to $lt*
                     %rs = cmpxchg $lt* %ptr, $lt %1, $lt %2 acq_rel acquire
                     %rv = extractvalue { $lt, i1 } %rs, 0
                     ret $lt %rv
                     """, $typ, Tuple{Ptr{$typ},$typ,$typ},
                     ref.pointer, cmp, new)
        end end
    else
        @eval Threads.atomic_cas!(ref::AtomicRef{$typ}, cmp::$typ, new::$typ) = let __x = ref.x; GC.@preserve __x begin
            llvmcall($"""
                     %iptr = inttoptr i$WORD_SIZE %0 to $ilt*
                     %icmp = bitcast $lt %1 to $ilt
                     %inew = bitcast $lt %2 to $ilt
                     %irs = cmpxchg $ilt* %iptr, $ilt %icmp, $ilt %inew acq_rel acquire
                     %irv = extractvalue { $ilt, i1 } %irs, 0
                     %rv = bitcast $ilt %irv to $lt
                     ret $lt %rv
                     """, $typ, Tuple{Ptr{$typ},$typ,$typ},
                     ref.pointer, cmp, new)
        end end
    end

    arithmetic_ops = [:add, :sub]
    for rmwop in [arithmetic_ops..., :xchg, :and, :nand, :or, :xor, :max, :min]
        rmw = string(rmwop)
        fn = Symbol("atomic_", rmw, "!")
        if (rmw == "max" || rmw == "min") && typ <: Unsigned
            # LLVM distinguishes signedness in the operation, not the integer type.
            rmw = "u" * rmw
        end
        if rmwop in arithmetic_ops && !(typ <: ArithmeticTypes) continue end
        if typ <: Integer
            @eval Threads.$fn(ref::AtomicRef{$typ}, v::$typ) = let __x = ref.x; GC.@preserve __x begin
                llvmcall($"""
                         %ptr = inttoptr i$WORD_SIZE %0 to $lt*
                         %rv = atomicrmw $rmw $lt* %ptr, $lt %1 acq_rel
                         ret $lt %rv
                         """, $typ, Tuple{Ptr{$typ}, $typ}, ref.pointer, v)
            end end
        else
            rmwop === :xchg || continue
            @eval Threads.$fn(ref::AtomicRef{$typ}, v::$typ) = let __x = ref.x; GC.@preserve __x begin
                llvmcall($"""
                         %iptr = inttoptr i$WORD_SIZE %0 to $ilt*
                         %ival = bitcast $lt %1 to $ilt
                         %irv = atomicrmw $rmw $ilt* %iptr, $ilt %ival acq_rel
                         %rv = bitcast $ilt %irv to $lt
                         ret $lt %rv
                         """, $typ, Tuple{Ptr{$typ}, $typ}, ref.pointer, v)
            end end
        end
    end
end


text:
using DataFrames
using CSV
using DataFramesMeta
using Query
using DataFramesMeta
using GlobalDecls
using Utils
using ONSCodes
@enum AggLevel national regional la

function loadone( name :: AbstractString ) :: DataFrame
    df = CSV.File(
        name,
        delim=',',
        missingstrings=["x","","-"],
        types=make_type_block(6:1000)) |> DataFrame
    lcnames = Symbol.(lowercase.(string.(names(df))))
    names!(df,lcnames)
    df
end

#
# DFE Raw data on care movements 2015-2018; see
# only 2017 and 2018 have Staying Put data
#
function getunderlyingdata()::NamedTuple
    (
    c17_18_2017=loadone(
            DATADIR*"underlying_data/2017/SFR50_CareLeavers17182017.csv" ),
    c19_21_2017 = loadone(
            DATADIR*"underlying_data/2017/SFR50_CareLeavers19to212017.csv" ),
    c17_18_2018 = loadone(
            DATADIR*"underlying_data/2018/CareLeavers17182018_amended.csv" ),
    c19_21_2018 = loadone(
            DATADIR*"underlying_data/2018/CareLeavers_Acc_StayPut19to212018.csv" ),
    c17_18_2019 = loadone(
            DATADIR*"underlying_data/2019/Care_Leavers_17182019.csv" ),
    c19_21_2019 = loadone(
            DATADIR*"underlying_data/2019/CareLeavers_Acc_StayPut19to212019.csv" )
    )
end

UDATA = getunderlyingdata()

"""
extract the number leaving care and the number staying put
for the given target (LA, country, etc.)
age (18,19,20) and year (2017,2018)
"""
function doquery(
    year   :: Integer,
    target :: AbstractString,
    age    :: Integer )
    as = age == 18 ? "17_18" : "19_21"
    # e.g. c17_18_2017
    field = Symbol( "c$(as)_$(year)" )
    f1 = Symbol( "cl_stayput_$age" )
    f2 = Symbol( "cl_stayput_ffc_$age" )
    println( "f1 $f1 f2 $f2 field $field" )
    q = @from i in UDATA[field] begin
        @where i.new_geog_code == target
        @select i[f1], i[f2]
        @collect
    end
    # @assert length( q ) == 1 no - 3 londons ffs we use London not Inner/Outer London
    println( "q1 $q[1]" )
    # return 1st row
    q[1]
end

"""
 return ONS code of either the LA, its region, or UK, depending on `AggLevel`
"""
function gettargetla(
    lacode :: AbstractString,
    agglev :: AggLevel ) :: AbstractString
    target = lacode
    if agglev == regional
        target = region_code_from_code( lacode )
    elseif agglev == national
        target = "E92000001" # england
    end
    target
end

"""
return the exit rates for the given LA - a named tuple a18,a19,a20
  will be the regional one if the calculation is not possible for an LA
  rates are proportions for that age - so 0.5,0.5,0.5 means 50% stay 1 year, 25% 2 12.5% 3 and so on
"""
function get_exit_rates(
    lacode    :: AbstractString,
    agglev    :: AggLevel,
    poolyears :: Bool  ) :: Vector
    println( "get_exit_rates for lacode $lacode agglev $agglev" );
    avprop = zeros(3)
    target = gettargetla( lacode, agglev )
    ages = [18,19,20]
    years = poolyears ? [2017,2018,2019] : [GlobalDecls.SIMULATION_YEAR]
    for year in years
        i = 0
        numleaving=zeros(3)
        numstaying=zeros(3)
        prop=zeros(3)
        for age in ages
            i+=1
            q = doquery( year, target, age )
            if zeroormissing( q )
                @assert agglev == la
                target = gettargetla( lacode, regional )
                println( "zeros detected! lacode $lacode new target $target $year age $age " );
                q = doquery( year, target, age )
            end
            @assert (! zeroormissing( q ))
            numleaving[i] += get(q[1],0)
            numstaying[i] += get(q[2],0)
        end
        prop = numstaying./numleaving
        prop[2] /= prop[1]
        prop[3] /= prop[2]
        println( "prop $prop " );
        println( "numstaying $numstaying numleaving $numleaving" )
        avprop .+= prop
    end
    avprop /= length( years )
    # t = (a18=avprop[1],a19=avprop[2],a20=avprop[3])
    return avprop
end

# df meta example:
# @linq las |> where( :LAD19NM.=="City of London" ) |> select( :LAD19CD )

@enum Region london se rest_of_england

function pay_class_from_region_code( regcode :: AbstractString ) :: Region
    # E12000001  #	A 	North East England
    # E12000002  # 	B 	North West England
    # E12000003  #	D 	Yorkshire and the Humber
    # E12000004  #	E 	East Midlands
    # E12000005  # 	F 	West Midlands
    # E12000006  # 	G 	East of England
    if regcode == "E12000007"
        r = london # 	H 	London
    elseif regcode == "E12000008"
        r = se
        # J 	South East England
    else
        r = rest_of_england
    end
    # E12000009 	K 	South West England
    r
end

function is_aggregate( ccode )
    if ccode === missing
        return true
    end
    if ccode == Nothing
        return true
    end
    if length(ccode)==0
        return true
    end
    if ccode[1:3] in ["E12", "E92" ]
        return true
    end
    return false

#     return name in [
#         "East of England",
#         "England",
#         "North East England",
#         "North West England",
#         "Yorkshire and The Humber",
#         "East Midlands",
#         "West Midlands",
#         "London",
#         "Inner London",
#         "Outer London",
#         "South East England",
#         "South West England",
#         "London Tri-borough" ]
end


text:
function test_mathfuns(::Type{T}) where T <: CasadiSymbolicObject
    @testset "$( string("Math functions for ", T, "                            ") )" begin
        x = randn()

        @test all(to_julia.( sincos( T(x) ) ) .≈ sincos(x))
        @test to_julia( sinc( T(x) ) )    ≈ sinc(x/π)
        @test to_julia( abs( T(x) ) )     ≈ abs(x)
    end
end


text:
# ---
# title: 188. Best Time to Buy and Sell Stock IV
# id: problem188
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/best-time-to-buy-and-sell-stock-iv/description/>
# hidden: true
# ---
#
# You are given an integer array `prices` where `prices[i]` is the price of a
# given stock on the `ith` day.
#
# Design an algorithm to find the maximum profit. You may complete at most `k`
# transactions.
#
# **Notice** that you may not engage in multiple transactions simultaneously
# (i.e., you must sell the stock before you buy again).
#
#
#
# **Example 1:**
#
#
#
#     Input: k = 2, prices = [2,4,1]
#     Output: 2
#     Explanation: Buy on day 1 (price = 2) and sell on day 2 (price = 4), profit = 4-2 = 2.
#
#
# **Example 2:**
#
#
#
#     Input: k = 2, prices = [3,2,6,5,0,3]
#     Output: 7
#     Explanation: Buy on day 2 (price = 2) and sell on day 3 (price = 6), profit = 6-2 = 4. Then buy on day 5 (price = 0) and sell on day 6 (price = 3), profit = 3-0 = 3.
#
#
#
#
# **Constraints:**
#
#   * `0 <= k <= 109`
#   * `0 <= prices.length <= 1000`
#   * `0 <= prices[i] <= 1000`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
"""
    evaluate_error(eq::Equation, scenario::Scenario, grid::Grid, t::Float64)

Prints out the L1/L2/L∞ integral norms of the error for
equation `eq`, scenario (with analytical solution) `scenario`, grid
`grid` at time `t`.
Error is defined as difference between analytical solution and numerical
approximation.
"""
function evaluate_error(eq::Equation, scenario::Scenario, grid::Grid, t::Float64)
    @assert is_analytical_solution(eq, scenario)
    l1_error = similar(grid.dofs, get_ndofs(eq))
    l2_error = similar(l1_error)
    linf_error = similar(l1_error)
    l1_error .= 0
    l2_error .= 0
    linf_error .= 0

    # TODO "Project solution to basis" = comp. average
    # Compare cell average with exact average

    #dg_order = size(grid.basis,1)
    dg_order = 1
    quad_order = 3
    quadpoints, quadweights = get_quadpoints(quad_order)
    #quad_order = max(dg_order + 2, 6)
    #quad_basis = Basis(quad_order, grid.basis.dimensions)
    number_of_quadpoints = length(quadpoints)

    for cell in grid.cells
        # Constant cell averages -> only one value
        for var=1:get_ndofs(eq)
            @views value = grid.dofs[var, cell.dataidx]
            for i=1:number_of_quadpoints
                for j=1:number_of_quadpoints
                    x = quadpoints[i]
                    y = quadpoints[j]
                    quadweight = quadweights[i] * quadweights[j]
                    weight = volume(cell) * quadweight
                    pxg, pyg = globalposition(cell, (x,y))
                    analytical = get_initial_values(eq, scenario, (pxg,pyg), t=t)
                    error = abs(value - analytical[var])
                    l1_error[var] += weight * error
                    l2_error[var] += weight * error^2
                    linf_error[var] = max(linf_error[var], error)
                end
            end
        end
    end

    l2_error .= sqrt.(l2_error)
    println("Errors for each variable.")
    println("Var\tL1\tL2\tL∞")
    for i=1:get_ndofs(eq)
        var_name = string(get_variable_name(eq,i))
        @printf "%s\t%15.6e\t%15.6e\t%15.6e\n" var_name l1_error[i] l2_error[i] linf_error[i]
    end

end


text:
@testset "file1 gets run" begin
    global file1_run
    @test !file1_run
    file1_run = true
    @test file1_run
end


text: Dict{String, Any}("albums" => Dict{String, Any}("songs" => Any[Dict{String, Any}("name" => Dict{String, Any}("value" => "Glory Days", "type" => "string"))]))


text:
##  Copyright 2021 Myrl Marmarelis
##
##  Licensed under the Apache License, Version 2.0 (the "License");
##  you may not use this file except in compliance with the License.
##  You may obtain a copy of the License at
##
##    http://www.apache.org/licenses/LICENSE-2.0
##
##  Unless required by applicable law or agreed to in writing, software
##  distributed under the License is distributed on an "AS IS" BASIS,
##  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
##  See the License for the specific language governing permissions and
##  limitations under the License.

using StatsBase
using StaticArrays
using PDMats
using Distributions
using LinearAlgebra

# `N` is the number of candidates, summed up over all issues (should there be more than one issue).
# `M` is the number of mixtures. `V` is the number of voters.
# `R` is the size of the ranking (i.e. number of choices.)

## possibly for the future, note that one could use `... where C <: SMatrix{N,N,T}`
## in order to avoid specifying the `L` explicitly, which can't be computed within the type
## arguments .. (no "dependent types"). can't just omit the type because that would make it
## dynamically resolved for every instance..

## if `matrix` is an SMatrix, no new allocations (should) occur
#Base.convert(::Type{PDMat{T,SMatrix{N,N,T,L}}}, matrix::AbstractMatrix{T}) where {N,T,L} =
#  PDMat{T,SMatrix{N,N,T,L}}(N, matrix, cholesky(matrix)) # assert L == N*N?
# these things graduated to the level of a massive pain...

const VoterCohort{ N, T, L } =
  MvNormal{ T, PDMat{ T, SMatrix{ N, N, T, L } }, SVector{ N, T } }
# ^ allow structures and functions that operate thereupon to fully specialize
# careful introspection revealed that this enabled massive performance gains
# by stabilizing the types in critical loops

# nice-to-have: enable cohorts to assign values to more than one set of candidates, i.e. question.
#   instead of sampling for each separately, it would be valuable to force the cohorts to correspond
#   across the two or more questions. alternatively, find post-hoc couplings between independent samples
#   my past attempt was to sample realizations for each issue separately, conditioning on other issues' realizations.
#   how this precisely squares with the cross-issue covariance structure is tenuous.
#   I have a hunch that it extends, and abides fully with my Gibbs scheme.
struct VoterMixture{ N, M, T, L } # all locally and statically allocated
  cohorts :: SVector{ M, VoterCohort{ N, T, L } }
  shares :: SVector{ M, T } # must sum to unit. the composition hereof
end

# ranked choices come in vectors that need not be as long as the total number of candidates.
# in that case, the constraints regarding the unincluded (~not quite excluded~) candidates are
# that their specific utility is less than that of the included ones.
# ordered list of candidate indicies, descending in utility
const RankedChoice{ R } = SVector{ R, Int } #where { R, T <: Real } rather than RankedChoice{ R, T }
# better semantics than the implicity of a matrix?
# strongly typed and explicitly-declared type parameters

struct IssueVote{ R }
  choices :: Vector{ RankedChoice{ R } }
  n_candidates :: Int # all issue `n_candidates` sum to the cohorts' big `N`
end

# encode in such a way that the entire nested type-structure
# (e.g. each array's heterogeneous `R`) becomes concrete at runtime
const MultiIssueVote = Tuple{ Vararg{ IssueVote } }

const Utility{ N, T } = SVector{ N, T }

const Membership = Int # if only I could easily parametrize the integer to live within a certain range, and define operations and promotion rules with it

struct VoterRealization{ N, M, T }
  membership :: Membership # Membership{M} ?
  utility :: Utility{ N, T }
end

Base.convert(::Type{VoterRealization{N,M,T}},
    params::Tuple{Int, AbstractVector{T}}) where {N,M,T} = # to allow easy array instantiations
  VoterRealization{N,M,T}(params[1], Utility{N,T}(params[2]))

Base.zero(::Type{VoterRealization{N,M,T}}) where {N,M,T} =
  VoterRealization{N,M,T}(0, zero(Utility{N,T}))

function validate_choices(rankings::AbstractVector{RankedChoice{R}}, N::Int)::Bool where R
  all( all( (ranking .>= 1) .& (ranking .<= N) ) for ranking in rankings )
end

# the shipped inv(::PDMat) is broken for SMatrices (that are not ::StridedVecOrMat)
Distributions.invcov(cohort::VoterCohort) = cohort.Σ.chol |> inv

## predictive posterior of these voter realizations
# for issue A and B:
# KL(joint(A,B) || single(A) x single(B)) is a sort of mutual information?
# consider this quantity more deeply.
# different from just taking one joint and splitting it because we *may*
# seek to study coupled vs. uncoupled


text:
module Metaprogramming

import Base.+, Base.-

export True, False, TrueOrFalse
export VoidTuple, MInt, MetaPair, Vals

# Make bools a bit easier to deal with
typealias True Val{true}
typealias False Val{false}
typealias TrueOrFalse Union{Type{Val{true}},Type{Val{false}}}

typealias VoidTuple{N} NTuple{N,Void}

# Stolen from @mbauman but doesn't work in-general
#check_Tuple(T) = (T===Tuple || T===NTuple) && throw(ArgumentError("parameters of $T are undefined"))
function concatenate{T<:Tuple, S<:Tuple}(::Type{T}, ::Type{S})
    Base.@_pure_meta
    Tuple{T.parameters..., S.parameters...}
end

# Some instantiators - people can define a "default" instantiator for a type which won't require too much computation for Julia
instantiate{T}(::Type{T}) = T()
instantiate{T<:Number}(::Type{T}) = zero(T)
instantiate{T<:AbstractString}(::Type{T}) = T(fieldtype(T,1)[])
instantiate(::Type{Symbol}) = symbol("")
instantiate{T}(::Type{Type{T}}) = T # ?

#"""
#By default, Julia provides limited constructors for `Tuple`-types. This function
#provides
#"""
#function instantiate{T<:Tuple}(::Type{T})

"""
A meta-integer `MInt{N}` is a special Julia type that naturally represents
unsigned integers `N` and can be added, multiplied, etc. Internally it is
manipulated using a tuple of `N` `nothing` values in order to make all
operations type-safe and have zero run-time penalty.

NOTE: Current Julia limitations mean this approach only works until N=8
"""
immutable MInt{N}; end

MInt_from_VoidTuple{N}(::VoidTuple{N}) = MInt{N}
MInt_from_VoidTuple(x::Void...) = MInt_from_VoidTuple(x)
VoidTuple_from_MInt{N}(::Type{MInt{N}}) = _voidtuple(Val{N},())

 _voidtuple{N}(::Type{Val{N}},x::VoidTuple{N}) = x
 _voidtuple{N,M}(::Type{Val{N}},x::VoidTuple{M}) = _voidtuple(Val{N},(nothing,x...))

Base.convert{N}(::Type{Val},::Type{MInt{N}}) = Val{N}
Base.convert{N}(::Type{Val{N}},::Type{MInt{N}}) = Val{N}
Base.convert{N}(::Type{MInt},::Type{Val{N}}) = MInt{N}
Base.convert{N}(::Type{MInt{N}},::Type{Val{N}}) = MInt{N}

# Base.promote_type(::Type{Val},::Type{MInt}) = ?

# Generated version.
#@generated +{N,M}(x::Type{MInt{N}},y::Type{MInt{M}}) = :(MInt{$(N+M)})
#@generated -{N,M}(x::Type{MInt{N}},y::Type{MInt{M}}) = :(MInt{$(N-M)})

# Pure version for Julia 0.5 (not working yet as of 3/3/2016)
function +{N,M}(x::Type{MInt{N}},y::Type{MInt{M}})
    Base.@_pure_meta
    MInt{N+M}
end
function -{N,M}(x::Type{MInt{N}},y::Type{MInt{M}})
    Base.@_pure_meta
    MInt{N-M}
end

# Hacky method making use of NTuple{N,Void}, Brittle and
#+{N,M}(x::Type{MInt{N}},y::Type{MInt{M}}) = MInt_from_VoidTuple(VoidTuple_from_MInt(x)...,VoidTuple_from_MInt(y)...)
#-{N,M}(x::Type{MInt{N}},y::Type{MInt{M}}) = MInt_from_VoidTuple(_subtract(Val{M},(),VoidTuple_from_MInt(x)...))

_subtract{N}(::Type{Val{N}},::VoidTuple{N},nothing,out::Void...) = (nothing, out...)
_subtract{N,M}(::Type{Val{N}},in::VoidTuple{M},nothing,out::Void...) = _subtract(Val{N},(nothing,in...),out...)
_subtract{N}(::Type{Val{N}},in::VoidTuple{N}) = ()
_subtract{N,M}(::Type{Val{N}},in::VoidTuple{M}) = error("Subtraction error: MInt must be positive")

"""
The `MetaPair` type is used for recursive-tree representation of one-or-more
type parameters, particularly used internally in `Vals`
"""
immutable MetaPair{T,M}; end
MetaPair() = Void
MetaPair{T}(::Type{T}) = MetaPair{T,Void}
MetaPair{T}(::Type{T},x::Type...) = MetaPair{T,MetaPair(x...)}

Base.show{T}(io::IO,::Type{MetaPair{T,Void}}) = show(io::IO,T)
function Base.show{T1,T2}(io::IO,::Type{MetaPair{T1,T2}})
    show(io::IO,T1)
    print(io::IO,",")
    show(io::IO,T2)
end

# The above was a little sneaky, so perhaps the user should be able to see the "real" thing when necessary
Base.showall{T}(io::IO,::Type{MetaPair{T,Void}}) = print(io,"MetaPair{$T,Void}")
function Base.showall{T1,T2}(io::IO,::Type{MetaPair{T1,T2}})
    print(io,"MetaPair{")
    showall(io,T1)
    print(io,",")
    showall(io,T2)
    print(io,"}")
end

_length{T,N}(::Type{MetaPair{T,Void}}, ::Type{MInt{N}}) = MInt{N} + MInt{1}
_length{T1,T2,N}(::Type{MetaPair{T1,T2}}, ::Type{MInt{N}}) = _length(T2, MInt{N} + MInt{1})

_getindex{T1,T2,N,M}(::Type{MetaPair{T1,T2}}, ::Type{Val{N}}, counter::Type{MInt{M}}) = _getindex(T2,Val{N},counter+MInt{1},X)
_getindex{T,N}(::Type{MetaPair{T,Void}}, ::Type{Val{N}}, counter::Type{MInt{N}}) = T
_getindex{T,N,M}(::Type{MetaPair{T,Void}}, ::Type{Val{N}}, counter::Type{MInt{M}}) = throw(BoundsError(X,N))

_getindex{T1,T2,N,X}(::Type{MetaPair{T1,T2}}, ::Type{Val{N}}, ::Type{MInt{N}},::Type{X}) = T1
_getindex{T1,T2,N,M,X}(::Type{MetaPair{T1,T2}}, ::Type{Val{N}}, counter::Type{MInt{M}},::Type{X}) = _getindex(T2,Val{N},counter+MInt{1},X)
_getindex{T,N,X}(::Type{MetaPair{T,Void}}, ::Type{Val{N}}, counter::Type{MInt{N}},::Type{X}) = T
_getindex{T,N,M,X}(::Type{MetaPair{T,Void}}, ::Type{Val{N}}, counter::Type{MInt{M}},::Type{X}) = throw(BoundsError(X,N))

_push!{T,X}(::Type{MetaPair{T,Void}},::Type{X}) = MetaPair{T,MetaPair{X,Void}}
_push!{T1,T2,X}(::Type{MetaPair{T1,T2}},::Type{X}) = MetaPair{T1,_push!(T2,X)}

_pop!{T}(::Type{MetaPair{T,Void}}) = (T,Void)
_pop!{T1,T2}(::Type{MetaPair{T1,T2}}) = _pop!(T2,MetaPair{T1,Void})
_pop!{T,Prev<:MetaPair}(::Type{MetaPair{T,Void}},::Type{Prev}) = (T,Prev)
_pop!{T1,T2,Prev<:MetaPair}(::Type{MetaPair{T1,T2}},::Type{Prev}) = _pop!(T2,_push!(Prev,T1))

_splat{T}(::Type{MetaPair{T,Void}}) = (T,)
_splat{T,M<:MetaPair}(::Type{MetaPair{T,M}}) = (T,_splat(M)...)

_splat_instances{T}(::Type{MetaPair{T,Void}}) = (instantiate(T),)
_splat_instances{T,M<:MetaPair}(::Type{MetaPair{T,M}}) = (instantiate(T),_splat_instances(M)...)

"""
`Vals` represents a meta-collection, with some similarities to `Tuple`, but
designed to represent a collection of `Val` or or other types. Convenience
functions for building, indexing, etc are provided.
"""
immutable Vals{T<:Union{MetaPair,Void}}; end
Vals() = Vals{Void}
Vals{T}(::Type{T}) = Vals{MetaPair{T,Void}}
Vals{T}(::Type{T},x::Type...) = Vals{MetaPair{T,MetaPair(x...)}}

function Base.showall{T}(io::IO,::Type{Vals{T}})
    print("Vals{")
    showall(T)
    print("}")
end

Base.getindex{N}(::Type{Vals{Void}},::Union{Type{Val{N}},Type{MInt{N}}}) = throw(BoundsError(Vals{Void},N))
Base.getindex{T,N}(x::Type{Vals{T}},::Union{Type{Val{N}},Type{MInt{N}}}) = _getindex(T,Val{N},MInt{1},Vals{T})

Base.endof(::Type{Vals{Void}}) = Val{0}
Base.endof{T}(::Type{Vals{T}}) = _length(T, MInt{0})

Base.push!{T}(::Type{Vals{Void}},::Type{T}) = Vals{MetaPair{T,Void}}
Base.push!{M,T}(::Type{Vals{M}},::Type{T}) = Vals{_push!(M,T)}

Base.pop!(::Type{Vals{Void}}) = throw(ArgumentError("Vals{Void} must be non-empty"))
function Base.pop!{M}(::Type{Vals{M}})
    v,m = _pop!(M)
    v,Vals{m}
end

Base.shift!(::Type{Vals{Void}}) = throw(ArgumentError("Vals{Void} must be non-empty"))
Base.shift!{T,M}(::Type{Vals{MetaPair{T,M}}}) = (T,Vals{M})

Base.unshift!{T,M}(::Type{Vals{M}},::Type{T}) = Vals{MetaPair{T,M}}

# These seem to have strange return type:: Tuple{DataType,DataType,...} and not Tuple{Type{Int64},...} or whatever...
# Also generate a LOT of code... (maybe that's splatting overhead that disappears in Julia 0.5?)
tuple_from_Vals(::Type{Vals{Void}}) = ()
tuple_from_Vals{M<:MetaPair}(::Type{Vals{M}}) = (_splat(M)...)

# OLD VERSION: It seemed to not be type-stable (possibly a Julia deficiency?)
#Tuple_from_Vals(::Type{Vals{Void}}) = Tuple{}
#Tuple_from_Vals{T}(::Type{Vals{MetaPair{T,Void}}}) = Tuple{T}
#Tuple_from_Vals{M<:MetaPair}(::Type{Vals{M}}) = Tuple{_splat(M)...}

# NEW VERSION: Is type-stable but only generates a no-op for isbits types excluding symbols (so symbols, arrays, strings, Ref{}, etc cause overhead)
# TODO: See if we can make this work through some C-calls or other hacking
Tuple_from_Vals(::Type{Vals{Void}}) = Tuple{}
Tuple_from_Vals{M<:MetaPair}(::Type{Vals{M}}) = typeof(_splat_instances(M))

end # module


text:
using InstantiateFromURL
# optionally add arguments to force installation: instantiate = true, precompile = true
github_project("QuantEcon/quantecon-notebooks-julia", version = "0.8.0")

using LinearAlgebra, Statistics
using Plots, Distributions, Random, Statistics
gr(fmt = :png, size = (900, 500))

function ksl(distribution, n = 100)
    title = nameof(typeof(distribution))
    observations = rand(distribution, n)
    sample_means = cumsum(observations) ./ (1:n)
    μ = mean(distribution)
    plot(repeat((1:n)', 2),
         [zeros(1, n); observations'], label = "", color = :grey, alpha = 0.5)
    plot!(1:n, observations, color = :grey, markershape = :circle,
          alpha = 0.5, label = "", linewidth = 0)
    hline!([μ], color = :black, linewidth = 1.5, linestyle = :dash, grid = false, label = ["Mean"])
    plot!(1:n, sample_means, linewidth = 3, alpha = 0.6, color = :green,
          label = "Sample mean")
    return plot!(title = title)
end

# Example distributions
distributions = [TDist(10), Beta(2,2), Gamma(5,2), Poisson(4), LogNormal(0.5), Exponential(1)]

Random.seed!(0); #reproducible results
ksl(Normal())
plot(ksl.(sample(distributions, 3, replace = false))..., layout = (3, 1), legend = false)

# Cauchy distribution may not converge as E|X| < ∞ not satisfied (heavy tailed)

function plotmeans(distribution::Distribution{Univariate,Continuous};n=1000)
    sample_mean  = cumsum(rand(distribution,n)) ./ (1:n)
    plot(1:n, sample_mean, label = "sample mean")
    hline!([0], color = :black, linestyle = :dash)
end

Random.seed!(0);
plotmeans(Cauchy(),n=10000)

using StatsPlots

function simulation1(distribution::Distribution{Univariate,Continuous}; n=250, k = 10_000)
    σ = std(distribution)
    y = rand(distribution,n,k)
    y .-= mean(distribution)
    y = mean(y,dims=1)
    y = √n *vec(y)
    density(y,label = "Empirical Distribution")
    return plot!(Normal(0,σ), linestyle = :dash, label = "Normal(0.0, $(σ^2))")
end

Random.seed!(0);
simulation1(Exponential(0.5), n=100)

function simulation2(distribution = Beta(2, 2), n = 5, k = 10_000)
    y = rand(distribution, k, n)
    for col in 1:n
        y[:,col] += rand([-0.5, 0.6, -1.1], k)
    end
    y = (y .- mean(distribution)) ./ std(distribution)
    y = cumsum(y, dims = 2) ./ sqrt.(1:5)' # return grid of data
end

ys = simulation2()
plots = [] # would preallocate in optimized code
for i in 1:size(ys, 2)
    p = density(ys[:, i], linealpha = i, title = "n = $i")
    push!(plots, p)
end

plot(plots..., legend = false)

# Exercise 1

function exercise1(a,b; g = sin, g′ = cos, n=250, k=10_000)
    distribution = Uniform(a,b)
    μ, σ = mean(distribution), std(distribution)
    gμ = sin.(mean(distribution))
    X = rand(distribution,n,k)
    X = mean(X, dims = 1)
    X = vec(X)
    y = √n .* vec(g.(X) .- g.(μ))
    density(y,label = "Empirical Distribution")
    return plot!(Normal(0,g′(μ).*σ), linestyle = :dash, label = "Asymptotic")
end

exercise1(0.0,π/2,k=15000)

# Exercise 2
function exercise2(;n = 250, k = 50_000, dw = Uniform(-1, 1), du = Uniform(-2, 2))
    vw = var(dw)
    vu = var(du)
    Σ = [vw vw; vw vu+vw]
    Q = inv(sqrt(Σ))

    function generate_data(dw, du, n)
        dw = rand(dw, n)
        X = [dw dw + rand(du, n)]
        return sqrt(n) * mean(X, dims = 1)
    end

    X = mapreduce(x -> generate_data(dw, du, n), vcat, 1:k)
    X = sqrt(n) .* mean(X, dims=1)
    X = Q*X'
    X = sum(abs2, X, dims = 1)
    X = vec(X)
    density(X, label = "", xlim = (0, 10))
    return plot!(Chisq(2), color = :black, linestyle = :dash,
                 label = "Chi-squared with 2 degrees of freedom", grid = false)
end
exercise2()

max_stars_repo_name: shanemcmiken/quantecon-notebooks-julia

max_issues_repo_name: shanemcmiken/quantecon-notebooks-julia

max_forks_repo_name: shanemcmiken/quantecon-notebooks-julia