text:
type ConfusionMatrix
  classes::Vector
  matrix::Matrix{Int}
  accuracy::FloatingPoint
  kappa::FloatingPoint
end

function show(io::IO, cm::ConfusionMatrix)
  print(io, "Classes: ")
  show(io, cm.classes)
  print(io, "\nMatrix: ")
  display(cm.matrix)
  print(io, "\nAccuracy: ")
  show(io, cm.accuracy)
  print(io, "\nκ: ")
  show(io, cm.kappa)
end

## measures ##

function confusion_matrix(actual::Vector, predicted::Vector)
    @assert length(actual) == length(predicted)
    N = length(actual)
    _actual = zeros(Int,N)
    _predicted = zeros(Int,N)
    classes = sort(unique([actual; predicted]))
    N = length(classes)
    for i in 1:N
        _actual[actual .== classes[i]] = i
        _predicted[predicted .== classes[i]] = i
    end
    CM = zeros(Int,N,N)
    for i in zip(_actual, _predicted)
        CM[i[1],i[2]] += 1
    end
    accuracy = trace(CM) / sum(CM)
    prob_chance = (sum(CM,1) * sum(CM,2))[1] / sum(CM)^2
    kappa = (accuracy - prob_chance) / (1.0 - prob_chance)
    return ConfusionMatrix(classes, CM, accuracy, kappa)
end

function mean_squared_error(actual, predicted)
    @assert length(actual) == length(predicted)
    return mean((actual - predicted).^2)
end

function R2(actual, predicted)
    @assert length(actual) == length(predicted)
    ss_residual = sum((actual - predicted).^2)
    ss_total = sum((actual .- mean(actual)).^2)
    return 1.0 - ss_residual/ss_total
end

## helper functions ##

function _kfoldCV{T<:FloatingPoint, U<:Real}(regressor::Symbol, labels::Vector{T}, features::Matrix{U}, nfolds::Integer,
                                             n_components::Integer, nferns::Integer)
    N = length(labels)
    ntest = _int(floor(N / nfolds))
    inds = randperm(N)
    R2s = zeros(nfolds)
    for i in 1:nfolds
        test_inds = falses(N)
        test_inds[(i - 1) * ntest + 1 : i * ntest] = true
        train_inds = !test_inds
        test_features = features[inds[test_inds],:]
        test_labels = labels[inds[test_inds]]
        train_features = features[inds[train_inds],:]
        train_labels = labels[inds[train_inds]]
        if regressor == :fern
            model = build_fern(train_labels, train_features, n_components)
            predictions, w = apply_fern(model, test_features)
        elseif regressor == :ensemble
            model = build_ensemble(train_labels, train_features, n_components, nferns)
            predictions = apply_ensemble(model, test_features)
        end

        err = mean_squared_error(test_labels, predictions)
        corr = cor(test_labels, predictions)
        r2 = R2(test_labels, predictions)
        R2s[i] = r2

        println("\nFold ", i)
        println("Mean Squared Error:     ", err)
        println("Correlation Coeff:      ", corr)
        println("Coeff of Determination: ", r2)
    end
    println("\nMean Coeff of Determination: ", mean(R2s))
    return R2s
end

kfoldCV_fern{T<:FloatingPoint, U<:Real}(labels::Vector{T}, features::Matrix{U}, nfolds::Integer, n_components::Integer)      = _kfoldCV(:fern, labels, features, nfolds, n_components, 1)
kfoldCV_fernEnsemble{T<:FloatingPoint, U<:Real}(labels::Vector{T}, features::Matrix{U}, nfolds::Integer, n_components::Integer, n_ferns::Integer)  = _kfoldCV(:ensemble, labels, features, nfolds, n_components, n_ferns)


text:
"""
    struct Grid1D{Tx<:AbstractTopology} <: AbstractGrid

Returns a one-dimensional staggered `grid` with topology `Tx`.

$(TYPEDFIELDS)
"""
struct Grid1D{Tx<:AbstractTopology} <: AbstractGrid
    "Number of points in x-direction"
    nx::Int
    "Number of halo points in x-direction"
    hx::Int
    "Grid spacing in x-direction"
    dx::Float64
    "Domain extent in x-direction"
    Lx::Float64
    "Cell faces in x-direction"
    xF::AbstractArray
    "Cell centres in x-direction"
    xC::AbstractArray
end

"""
    Grid1D(Tx, nx, x_start, x_end; hx=1)

Construct a one-dimensional staggered `grid` on domain `x ∈ [x_start, x_end]` with topology
`Tx`, with `nx` interior grid points, and `hx` halo points.

Example
=======

```jldoctest
julia> using MixedLayerThermoclineDynamics

julia> grid = Grid1D(Periodic(), 10, 0, 2.0)
1-Dimensional Grid
  ├───────── topology: Periodic
  ├─ domain extent Lx: 2.0
  ├──── resolution nx: 10
  ├── grid spacing dx: 0.2
  └─── halo points nx: 1
```
"""
function Grid1D(Tx, nx, x_start, x_end; hx=1)

    if hx == 0; throw(error("Number of halo points in x cannot be zero")); end

    Lx = x_end - x_start

    dx = Lx/nx

    xF = construct_faces(Tx, nx, hx, dx, Lx, x_start)
    xC = construct_centres(Tx, nx, hx, dx, Lx, x_start)

    return Grid1D{typeof(Tx)}(nx, hx, dx, Lx, xF, xC)
end

"""
    struct Grid2D{Tx<:AbstractTopology, Ty<:AbstractTopology} <: AbstractGrid

Returns a two-dimensional staggered `grid` with topologies `{Tx, Ty}`.

$(TYPEDFIELDS)
"""
struct Grid2D{Tx<:AbstractTopology, Ty<:AbstractTopology} <: AbstractGrid
    "Number of points in x-direction"
    nx::Int
    "Number of points in y-direction"
    ny::Int
    "Number of halo points in x-direction"
    hx::Int
    "Number of halo points in y-direction"
    hy::Int
    "Grid spacing in x-direction"
    dx::Float64
    "Grid spacing in y-direction"
    dy::Float64
    "Domain extent in x-direction"
    Lx::Float64
    "Domain extent in y-direction"
    Ly::Float64
    "Cell faces in x-direction"
    xF::AbstractArray
    "Cell centres in x-direction"
    xC::AbstractArray
    "Cell faces in y-direction"
    yF::AbstractArray
    "Cell centres in y-direction"
    yC::AbstractArray
end

"""
    Grid2D(Tx, Ty, nx, ny, x_start, x_end, y_start, y_end; hx=1, hy=1)

Construct a two-dimensional staggered `grid` on domain `(x, y) ∈ [x_start, x_end] x [y_start, y_end]`
with topologies `{Tx, Ty}`, with `{nx, ny}` interior grid points, and `{hx, hy}` halo points.

Example
=======

```jldoctest
julia> using MixedLayerThermoclineDynamics

julia> grid = Grid2D(Periodic(), Periodic(), 10, 15, 0, 2.0, 0, 3.0)
2-Dimensional Grid
  ├──── topology in x: Periodic
  ├─ domain extent Lx: 2.0
  ├──── resolution nx: 10
  ├── grid spacing dx: 0.2
  ├─── halo points nx: 1
  ├──── topology in y: Periodic
  ├─ domain extent Ly: 3.0
  ├──── resolution ny: 15
  ├── grid spacing dy: 0.2
  └─── halo points ny: 1
```
"""
function Grid2D(Tx, Ty, nx, ny, x_start, x_end, y_start, y_end; hx=1, hy=1)

    if hx == 0; throw(error("Number of halo points in x cannot be zero")); end
    if hy == 0; throw(error("Number of halo points in y cannot be zero")); end

    Lx = x_end - x_start
    Ly = y_end - y_start

    dx, dy = Lx/nx, Ly/ny

    xF = construct_faces(Tx, nx, hx, dx, Lx, x_start)
    xC = construct_centres(Tx, nx, hx, dx, Lx, x_start)
    yF = construct_faces(Ty, ny, hy, dy, Ly, y_start)
    yC = construct_centres(Ty, ny, hy, dy, Ly, y_start)

    return Grid2D{typeof(Tx), typeof(Ty)}(nx, ny, hx, hy, dx, dy, Lx, Ly, xF, xC, yF, yC)
end

##
# helper functions

"""
    construct_faces(T::AbstractTopology, n, h, d, L)

Returns the face locations for a dimension with topology `T`, number of grid points `n` and `h`
halo points, grid spacing `d`, and extent `L`.
"""
function construct_faces(T::AbstractTopology, n, h, d, L, start)
    start_face  = start
      end_face  = isa(T, Periodic) ? start + L - d : start + L
    total_faces = isa(T, Periodic) ? n             : n + 1

    F = range(start_face - h*d, stop = end_face + h*d, length = total_faces + 2h)

    return OffsetArray(F, -h)
end

"""
    construct_centres(T::AbstractTopology, n, h, d, L, start_center)

Returns the cell centers locations for a dimension with topology `T`, number of grid points
`n` and `h` halo points, grid spacing `d`, and extent `L`.
"""
function construct_centres(T::AbstractTopology, n, h, d, L, start)
    start_center = start + d/2
      end_center = start + L - d/2

    C = range(start_center - h*d, stop = end_center + h*d, length = n + 2h)

    return OffsetArray(C, -h)
end

show(io::IO, grid::Grid1D{Tx}) where Tx =
     print(io, "1-Dimensional Grid\n",
               "  ├───────── topology: ", Tx, '\n',
               "  ├─ domain extent Lx: ", grid.Lx, '\n',
               "  ├──── resolution nx: ", grid.nx, '\n',
               "  ├── grid spacing dx: ", grid.dx, '\n',
               "  └─── halo points nx: ", grid.hx)

show(io::IO, grid::Grid2D{Tx, Ty}) where {Tx, Ty} =
     print(io, "2-Dimensional Grid\n",
               "  ├──── topology in x: ", Tx, '\n',
               "  ├─ domain extent Lx: ", grid.Lx, '\n',
               "  ├──── resolution nx: ", grid.nx, '\n',
               "  ├── grid spacing dx: ", grid.dx, '\n',
               "  ├─── halo points nx: ", grid.hx, '\n',
               "  ├──── topology in y: ", Ty, '\n',
               "  ├─ domain extent Ly: ", grid.Ly, '\n',
               "  ├──── resolution ny: ", grid.ny, '\n',
               "  ├── grid spacing dy: ", grid.dy, '\n',
               "  └─── halo points ny: ", grid.hy)

ext: jl


text:
sys.it.personal_allowance=24
sys.it.savings_thresholds = [7000,9000,222000]
sys.lmt.child_tax_credit.disability = 8282.0

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
module PWscf

using Test

using Setfield
using StructArrays: StructArray

using QuantumESPRESSOBase
using QuantumESPRESSOBase.Inputs.PWscf

@testset "Constructing `AtomicSpecies`" begin
    # Data from https://github.com/QEF/q-e/blob/7be27df/PW/examples/gatefield/run_example#L128.
    x = AtomicSpecies("S", 32.066, "S.pz-n-rrkjus_psl.0.1.UPF")
    @test_throws AssertionError @set x.atom = "sulfur"
    @test_throws InexactError @set x.mass = 1im
    @test x == AtomicSpecies('S', 32.066, "S.pz-n-rrkjus_psl.0.1.UPF")
    @test AtomicSpecies(
        AtomicPosition('S', [0.500000000, 0.288675130, 1.974192764]),
        32.066,
        "S.pz-n-rrkjus_psl.0.1.UPF",
    ) == x
end

@testset "Test constructing `AtomicSpeciesCard` from `StructArray`s" begin
    a = ["Al", "As"]
    m = [24590.7655930491, 68285.4024548272]
    pp = ["Al.pbe-n-kjpaw_psl.1.0.0.UPF", "As.pbe-n-kjpaw_psl.1.0.0.UPF"]
    card = AtomicSpeciesCard(StructArray{AtomicSpecies}((a, m, pp)))
    @test card.data == [
        AtomicSpecies("Al", 24590.7655930491, "Al.pbe-n-kjpaw_psl.1.0.0.UPF"),
        AtomicSpecies("As", 68285.4024548272, "As.pbe-n-kjpaw_psl.1.0.0.UPF"),
    ]
    push!(card.data, AtomicSpecies("Si", 25591.1924913552, "Si.pbe-n-kjpaw_psl.1.0.0.UPF"))
    @test card.data == [
        AtomicSpecies("Al", 24590.7655930491, "Al.pbe-n-kjpaw_psl.1.0.0.UPF"),
        AtomicSpecies("As", 68285.4024548272, "As.pbe-n-kjpaw_psl.1.0.0.UPF"),
        AtomicSpecies("Si", 25591.1924913552, "Si.pbe-n-kjpaw_psl.1.0.0.UPF"),
    ]
end

@testset "Constructing `AtomicPosition`" begin
    # Data from https://github.com/QEF/q-e/blob/7be27df/PW/examples/gatefield/run_example#L129-L132.
    x = AtomicPosition("S", [0.500000000, 0.288675130, 1.974192764])
    @test_throws AssertionError @set x.atom = "sulfur"
    @test_throws InexactError @set x.pos = [1im, 2im, 3im]
    @test_throws ErrorException x.posi = [0.1, 0.2, 0.3]  # Given a wrong field name
    @test x.if_pos == [1, 1, 1]
    @test x == AtomicPosition('S', [0.500000000, 0.288675130, 1.974192764])
    @test AtomicPosition(
        AtomicSpecies('S', 32.066, "S.pz-n-rrkjus_psl.0.1.UPF"),
        [0.500000000, 0.288675130, 1.974192764],
    ) == x
end

@testset "Test constructing `AtomicPositionsCard` from `StructArray`s" begin
    # Data from https://github.com/QEF/q-e/blob/7be27df/PW/examples/gatefield/run_example#L129-L132.
    a = ["S", "Mo", "S"]
    pos = [
        [0.500000000, 0.288675130, 1.974192764],
        [0.000000000, 0.577350270, 2.462038339],
        [0.000000000, -0.577350270, 2.950837559],
    ]
    card = AtomicPositionsCard(
        StructArray{AtomicPosition}((a, pos, [[1, 1, 1], [1, 1, 1], [1, 1, 1]])),
        "alat",
    )
    @test card.data == [
        AtomicPosition("S", [0.500000000, 0.288675130, 1.974192764]),
        AtomicPosition("Mo", [0.000000000, 0.577350270, 2.462038339]),
        AtomicPosition("S", [0.000000000, -0.577350270, 2.950837559]),
    ]
end

@testset "Constructing `CellParametersCard`" begin
    #Data from https://gitlab.com/QEF/q-e/blob/master/NEB/examples/neb1.in
    option = "bohr"
    data = [
        12 0 0
        0 5 0
        0 0 5
    ]
    card = CellParametersCard(data, option)
    @test_throws AssertionError @set card.option = "ala"
    @test_throws AssertionError @set card.option = "crystal" # Allowed options are alat, angstrom, bohr
    @test_throws DimensionMismatch @set card.data = [ # Matrix size should be (3, 3)
        1 2
        3 4
    ]
    @test_throws DimensionMismatch @set card.data = [
        1 2 3 4
        5 6 7 8
        4 3 2 1
        8 7 6 5
    ]
    @test CellParametersCard(data).option == "alat" # default option is alat
end

@testset "Constructing `AtomicForce`" begin
    x = AtomicForce("H", [1, 2, 3])
    @test_throws DimensionMismatch @set x.force = [1, 2]
    @test_throws DimensionMismatch @set x.force = [1, 2, 3, 4]
end

@testset "Test constructing `AtomicForce` from `StructArray`" begin
    a = ["H", "O", "H"]
    f = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
    card = AtomicForcesCard(StructArray{AtomicForce}((a, f)))
    @test card.data == [
        AtomicForce("H", [1, 2, 3]),
        AtomicForce("O", [4, 5, 6]),
        AtomicForce("H", [7, 8, 9]),
    ]
end

@testset "Test constructing a `PWInput`: silicon" begin
    # This example is from https://github.com/QEF/q-e/blob/master/PW/examples/example01/run_example.
    for diago in ("david", "cg", "ppcg")
        control = ControlNamelist(
            tstress = true,
            tprnfor = true,
            outdir = raw"$TMP_DIR/",
            prefix = "silicon",
            pseudo_dir = raw"$PSEUDO_DIR/",
        )
        system =
            SystemNamelist(ibrav = 2, celldm = [10.2], nat = 2, ntyp = 1, ecutwfc = 18.0)
        electrons = ElectronsNamelist(conv_thr = 1.0e-8, diagonalization = "$diago")
        atomic_species = AtomicSpeciesCard([AtomicSpecies("Si", 28.086, "Si.pz-vbc.UPF")])
        atomic_positions = AtomicPositionsCard([
            AtomicPosition("Si", [0.0, 0.0, 0.0]),
            AtomicPosition("Si", [0.25, 0.25, 0.25]),
        ])
        k_points = SpecialPointsCard(
            [
                0.125 0.125 0.125 1.0
                0.125 0.125 0.375 3.0
                0.125 0.125 0.625 3.0
                0.125 0.125 0.875 3.0
                0.125 0.375 0.375 3.0
                0.125 0.375 0.625 6.0
                0.125 0.375 0.875 6.0
                0.125 0.625 0.625 3.0
                0.375 0.375 0.375 1.0
                0.375 0.375 0.625 3.0
            ],
        )
        object = PWInput(
            control = control,
            system = system,
            electrons = electrons,
            atomic_species = atomic_species,
            atomic_positions = atomic_positions,
            k_points = k_points,
            cell_parameters = nothing,
        )
    end
end

@testset "Test constructing a `PWInput`: silicon bands" begin
    # This example is from https://github.com/QEF/q-e/blob/master/PW/examples/example01/run_example.
    for diago in ("david", "cg", "ppcg")
        control = ControlNamelist(
            calculation = "bands",
            pseudo_dir = raw"$PSEUDO_DIR/",
            outdir = raw"$TMP_DIR/",
            prefix = "silicon",
        )
        system = SystemNamelist(
            ibrav = 2,
            celldm = [10.2],
            nat = 2,
            ntyp = 1,
            ecutwfc = 18.0,
            nbnd = 8,
        )
        electrons = ElectronsNamelist(diagonalization = "$diago")
        atomic_species = AtomicSpeciesCard([AtomicSpecies("Si", 28.086, "Si.pz-vbc.UPF")])
        atomic_positions = AtomicPositionsCard([
            AtomicPosition("Si", [0.0, 0.0, 0.0]),
            AtomicPosition("Si", [0.25, 0.25, 0.25]),
        ])
        k_points = SpecialPointsCard(
            [
                0.0 0.0 0.0 1.0
                0.0 0.0 0.1 1.0
                0.0 0.0 0.2 1.0
                0.0 0.0 0.3 1.0
                0.0 0.0 0.4 1.0
                0.0 0.0 0.5 1.0
                0.0 0.0 0.6 1.0
                0.0 0.0 0.7 1.0
                0.0 0.0 0.8 1.0
                0.0 0.0 0.9 1.0
                0.0 0.0 1.0 1.0
                0.0 0.0 0.0 1.0
                0.0 0.1 0.1 1.0
                0.0 0.2 0.2 1.0
                0.0 0.3 0.3 1.0
                0.0 0.4 0.4 1.0
                0.0 0.5 0.5 1.0
                0.0 0.6 0.6 1.0
                0.0 0.7 0.7 1.0
                0.0 0.8 0.8 1.0
                0.0 0.9 0.9 1.0
                0.0 1.0 1.0 1.0
                0.0 0.0 0.0 1.0
                0.1 0.1 0.1 1.0
                0.2 0.2 0.2 1.0
                0.3 0.3 0.3 1.0
                0.4 0.4 0.4 1.0
                0.5 0.5 0.5 1.0
            ],
        )
        object = PWInput(
            control = control,
            system = system,
            electrons = electrons,
            atomic_species = atomic_species,
            atomic_positions = atomic_positions,
            k_points = k_points,
            cell_parameters = nothing,
        )
    end
end

@testset "Test constructing a `PWInput`: aluminium" begin
    # This example is from https://github.com/QEF/q-e/blob/master/PW/examples/example01/run_example.
    for diago in ("david", "cg", "ppcg")
        control = ControlNamelist(
            calculation = "scf",
            restart_mode = "from_scratch",
            pseudo_dir = raw"$PSEUDO_DIR/",
            outdir = raw"$TMP_DIR/",
            prefix = "al",
            tprnfor = true,
            tstress = true,
        )
        system = SystemNamelist(
            ibrav = 2,
            celldm = [7.50],
            nat = 1,
            ntyp = 1,
            ecutwfc = 15.0,
            occupations = "smearing",
            smearing = "marzari-vanderbilt",
            degauss = 0.05,
        )
        electrons = ElectronsNamelist(diagonalization = "$diago", mixing_beta = 0.7)
        atomic_species = AtomicSpeciesCard([AtomicSpecies("Al", 26.98, "Al.pz-vbc.UPF")])
        atomic_positions = AtomicPositionsCard([AtomicPosition("Al", [0.0, 0.0, 0.0])])
        k_points = SpecialPointsCard(
            [
                0.0625000 0.0625000 0.0625000 1.00
                0.0625000 0.0625000 0.1875000 3.00
                0.0625000 0.0625000 0.3125000 3.00
                0.0625000 0.0625000 0.4375000 3.00
                0.0625000 0.0625000 0.5625000 3.00
                0.0625000 0.0625000 0.6875000 3.00
                0.0625000 0.0625000 0.8125000 3.00
                0.0625000 0.0625000 0.9375000 3.00
                0.0625000 0.1875000 0.1875000 3.00
                0.0625000 0.1875000 0.3125000 6.00
                0.0625000 0.1875000 0.4375000 6.00
                0.0625000 0.1875000 0.5625000 6.00
                0.0625000 0.1875000 0.6875000 6.00
                0.0625000 0.1875000 0.8125000 6.00
                0.0625000 0.1875000 0.9375000 6.00
                0.0625000 0.3125000 0.3125000 3.00
                0.0625000 0.3125000 0.4375000 6.00
                0.0625000 0.3125000 0.5625000 6.00
                0.0625000 0.3125000 0.6875000 6.00
                0.0625000 0.3125000 0.8125000 6.00
                0.0625000 0.3125000 0.9375000 6.00
                0.0625000 0.4375000 0.4375000 3.00
                0.0625000 0.4375000 0.5625000 6.00
                0.0625000 0.4375000 0.6875000 6.00
                0.0625000 0.4375000 0.8125000 6.00
                0.0625000 0.4375000 0.9375000 6.00
                0.0625000 0.5625000 0.5625000 3.00
                0.0625000 0.5625000 0.6875000 6.00
                0.0625000 0.5625000 0.8125000 6.00
                0.0625000 0.6875000 0.6875000 3.00
                0.0625000 0.6875000 0.8125000 6.00
                0.0625000 0.8125000 0.8125000 3.00
                0.1875000 0.1875000 0.1875000 1.00
                0.1875000 0.1875000 0.3125000 3.00
                0.1875000 0.1875000 0.4375000 3.00
                0.1875000 0.1875000 0.5625000 3.00
                0.1875000 0.1875000 0.6875000 3.00
                0.1875000 0.1875000 0.8125000 3.00
                0.1875000 0.3125000 0.3125000 3.00
                0.1875000 0.3125000 0.4375000 6.00
                0.1875000 0.3125000 0.5625000 6.00
                0.1875000 0.3125000 0.6875000 6.00
                0.1875000 0.3125000 0.8125000 6.00
                0.1875000 0.4375000 0.4375000 3.00
                0.1875000 0.4375000 0.5625000 6.00
                0.1875000 0.4375000 0.6875000 6.00
                0.1875000 0.4375000 0.8125000 6.00
                0.1875000 0.5625000 0.5625000 3.00
                0.1875000 0.5625000 0.6875000 6.00
                0.1875000 0.6875000 0.6875000 3.00
                0.3125000 0.3125000 0.3125000 1.00
                0.3125000 0.3125000 0.4375000 3.00
                0.3125000 0.3125000 0.5625000 3.00
                0.3125000 0.3125000 0.6875000 3.00
                0.3125000 0.4375000 0.4375000 3.00
                0.3125000 0.4375000 0.5625000 6.00
                0.3125000 0.4375000 0.6875000 6.00
                0.3125000 0.5625000 0.5625000 3.00
                0.4375000 0.4375000 0.4375000 1.00
                0.4375000 0.4375000 0.5625000 3.00
            ],
        )
        object = PWInput(
            control = control,
            system = system,
            electrons = electrons,
            atomic_species = atomic_species,
            atomic_positions = atomic_positions,
            k_points = k_points,
            cell_parameters = nothing,
        )
    end
end

end


text:
immutable BlockBandedMatrix{T,RI,CI} <: AbstractMatrix{T}
    data::Vector{T}  # the entries
    rows::RI
    columns::CI

    l::Int # lower bandwidth ≥0
    u::Int # upper bandwidth ≥0

    columnlookup::Vector{Int}  # translates column first data entry
    firstrow::Vector{Int}     # gives first row index in

    function BlockBandedMatrix(dat::Vector{T},rws::RI,cls::CI,l::Int,u::Int)
        n = 0
        for J=1:length(columns),K=max(1,J-u):min(J+l,length(rows))
            n += rows[K]*columns[J]
        end

        if n != length(dat)
            error("Data must be consistant size with number of non-zeros")
        end

        m = 1
        curcol = 1

        columnlookup = Array(Int,0)
        firstrow = Array(Int,0)
        for J=1:length(columns)
            sm = sum(rows[max(1,J-u):min(J+l,length(rows))])
            nrws = sum(rows[max(1,J-u):J-1])
            for j=1:columns[J]
                push!(columnlookup,m)
                push!(firstrow,curcol-j-nrws+1)

                m += sm
                curcol += 1
            end
        end

        new(dat,rws,cls,l,u,columnlookup,firstrow)
    end
end

BlockBandedMatrix(dat::Vector,rws,cls,l,u) =
    BlockBandedMatrix{eltype(dat),typeof(rws),typeof(cls)}(dat,rws,cls,l,u)

Base.size(B::BlockBandedMatrix) = sum(B.rows),sum(B.columns)
Base.linearindexing{BBM<:BlockBandedMatrix}(::Type{BBM}) = Base.LinearSlow()
function Base.getindex(B::BlockBandedMatrix,k::Integer,j::Integer)
    if k > size(B,1) || j > size(B,2)
        throw(BoundsError())
    end

    fr = B.firstrow[j]
    nrws = j != size(B,2) ?
            B.columnlookup[j+1] - B.columnlookup[j] :
            length(B.data) - B.columnlookup[j] + 1

    if k < fr || k > fr + nrws - 1
        zero(eltype(B))
    else
        B.data[B.columnlookup[j] + k - fr]
    end
end

function Base.setindex!(B::BlockBandedMatrix,v,k::Integer,j::Integer)
    if k > size(B,1) || j > size(B,2)
        throw(BoundsError())
    end

    fr = B.firstrow[j]
    nrws = j != size(B,2) ?
            B.columnlookup[j+1] - B.columnlookup[j] :
            length(B.data) - B.columnlookup[j] + 1

    if k < fr || k > fr + nrws - 1
        if v!=0
            error("Cannot change entries outside block bands")
        end
    else
        B.data[B.columnlookup[j] + k - fr] = v
    end
end

for (op,bop) in ((:(Base.rand),:bbrand),(:(Base.zeros),:bbzeros),(:(Base.ones),:bbones))
    @eval begin
        function $bop{T}(::Type{T},rows,columns,l,u)
            n = 0
            for J=1:length(columns),K=max(1,J-u):min((J+l),length(rows))
                n += rows[K]*columns[J]
            end
            BlockBandedMatrix($op(T,n),rows,columns,l,u)
        end
        $bop(rows,columns,l,u) = $bop(Float64,rows,columns,l,u)
    end
end

rows=1:50
columns=1:50
l=u=1

B=bbrand(rows,columns,l,u)

sz=sizeof(eltype(B))

v=pointer(B.data)
α=1.0
Z=bbzeros(rows,columns,l,u+l)
    w=pointer(Z.data)
    for k=1:4
        N=B.columnlookup[k+1]-B.columnlookup[k]
        BLAS.axpy!(N,α,v+(B.columnlookup[k]-1)*sz,1,w+(B.columnlookup[k]-1)*sz,1)
    end

B
unsafe_load(v+1*sz)
Z
using SO
B

B[10,3]=5

B

v=B.data[1:3]
norm(v)

B.columnlookup[2]

B.columnlookup

B.firstrow

size(B)

k=j=1
    B.data[B.columnlookup[j] + k-1]

B.data

columnlookup = Array(Int,0)
    m=1
    push!(columnlookup,m)  # the first column begins at index 1
    for J=1:length(columns)
        sm = sum(rows[max(1,J-u):min(J+l,length(rows))])
        for j=1:columns[J]
            m += sm
            push!(columnlookup,m)
        end
    end

columns[2]

columnlookup

n = 0
    for J=1:length(columns),K=max(1,J-u):(J+l)
        n += rows[K]*columns[J]
    end

    B=BlockBandedMatrix(rand(n),1,1,rows,columns)

B.data

size(B)

T=eltype(data)
@time pointer_to_array(pointer(data)+sizeof(T)*2,(2,2))


text:
"""
    ρ_of(atoms::Shape{T})
Computes Volumetric Density
"""
ρ_of(atoms::Atom{T}) where T = atoms.N / Volume_of(atoms)

"""
    Volume_of(atoms::Shape{Cube})
"""
Volume_of(atoms::Atom{Cube}) = π * atoms.sizes^2 # sizes == length of Cube's Sides

"""
    Volume_of(atoms::Shape{Cube})
"""
Volume_of(atoms::Atom{Circle}) = atoms.sizes^3 # sizes == length of Cube's Sides

"""
    b₀_of(atoms::Shape{Cube})
"""
b₀_of(atoms::Atom{Cube}) = (ρ_of(atoms) .^ 2 * atoms.N / ((4π / 3) * (3 / (16π))^3))^(1 / 3)

"""
    estimate_E₀(laser)
"""
estimate_E₀(laser) = √(laser.s * (1 + 4(laser.Δ / Γ)^2) / 2)


text:
using Pandas
using Test

df = DataFrame(Dict(:name=>["a", "b"], :age=>[27, 30]))
age = values(df.age)
age[2] = 31
@test loc(df)[1, "age"] == 31

query(df, :(age!=27))  # Issue #26

df = read_csv(joinpath(dirname(@__FILE__), "test.csv"))
typeof(df)
@test isa(df, Pandas.DataFrame)

include("test_tabletraits.jl")

@test !isempty(df)
empty!(df)
@test isempty(df)

x = Series([3,5], index=[:a, :b])

@test x.a == 3
@test x[:a] == 3
@test loc(x)[:a] == 3
@test x.b == 5
@test iloc(x)[1] == 3
@test iloc(x)[2] == 5
@test length(x) == 2
@test values(x+1) == [4, 6]
@test sum(x) == 8
@test eltype(x) == Int64
@test all(iloc(x)[1:2]==x)

# Rolling
roll = rolling(Series([1,2,3,4,5]), 3)
@test isequal(values(mean(roll)), [NaN, NaN, 2.0, 3.0, 4.0])

# HDF
mktempdir() do dir
    store = HDFStore("$(dir)/store.h5")
    x = Series(["a", "b"])
    store["x"] = x
    x_fromstore = store["x"]
    @test values(x_fromstore) == ["a", "b"]
end


text:
module CanonicalizeRequest

import ..Layer, ..request
using ..Messages
using ..Strings.tocameldash!

"""
    request(CanonicalizeLayer, method, ::URI, headers, body) -> HTTP.Response

Rewrite request and response headers in Canonical-Camel-Dash-Format.
"""

abstract type CanonicalizeLayer{Next <: Layer} <: Layer end
export CanonicalizeLayer

function request(::Type{CanonicalizeLayer{Next}},
                 method::String, url, headers, body; kw...) where Next

    headers = canonicalizeheaders(headers)

    res = request(Next, method, url, headers, body; kw...)

    res.headers = canonicalizeheaders(res.headers)

    return res
end

canonicalizeheaders(h::T) where T = T([tocameldash!(k) => v for (k,v) in h])

end # module CanonicalizeRequest


text:
"A plain value type that is serializable/deserializable by the Julia `TOML` package."
const TomlTypePlain = Union{Integer, AbstractFloat, AbstractString, Dates.TimeType}

"A simple collection that is serializable/deserializable by the Julia `TOML` package."
const TomlType = Union{TomlTypePlain,
                       AbstractVector{<:TomlTypePlain},
                       AbstractDict{<:AbstractString, <:TomlTypePlain}
                      }

#NOTE: originally this was an alias for TOML.Internals.Printer.TOMLValue,
#    but that causes any kind of dictionary/vector to get caught,
#    even if TOML can't actually serialize it.
# Plus we probably don't want to depend on TOML internals like that.


text:
module Routino

using CCallHelp

export route, distance, Router, find_waypoint, clear_waypoints!, add_waypoint!, quickest_route, shortest_route, get_profile_names, close_router, free_xml, walk_to_distance

const BINARY = "/usr/bin/routino-router"
const DATADIR = "/trip/osm"
const LIB = "/usr/lib/libroutino.so.0"
const ROUTINO_API_VER = Cint(8)
const PROFILES = "/usr/share/routino/profiles.xml"
const TRANSLATIONS = "/usr/share/routino/translations.xml"

const DB = Ptr{Cvoid}
const Profile = Ptr{Cvoid}
const Waypoint = Ptr{Cvoid}
const Translation = Ptr{Cvoid}
const ProgressFunc= Ptr{Cvoid}

const PUchars = Ptr{Cuchar}
const ListString = Ptr{PUchars}

const LoLa = typeof((lo=0.0, la=0.0))
"""
    Router(profile="motorcar"; language="en", path=DATADIR, prefix="GB", profilexml=PROFILES, translationsxml=TRANSLATIONS)
Load the Routino Library and intialize it with its database and options.
This is the best way to use the library but was written after the shell out version
# Arguments
- `profile` - name of the profile to use, must be listed in the profilexml
- `language` - language to use, defaults to "en"
- `path` - directory containing the pem files generated by planetsplit - see the Routino docs
- `prefix` - Country prefix for the pem files defaults to "GB"
- `profilexml` - filename from which to load the profiles, defaults to "/usr/share/routino/profiles.xml"
- `translationsxml` filename from which to load the translations, defaults to "/usr/share/routino/translations.xml"
"""
struct Router
    db::Ptr{Cvoid}
    profile::Ptr{Cvoid}
    translation::Ptr{Cvoid}
    waypoints::Vector{Ptr{Cvoid}}
    function Router(profile="motorcar"; language="en", path=DATADIR, prefix="GB", profilexml=PROFILES, translationsxml=TRANSLATIONS)
        load_xml_profiles(profilexml)
        load_xml_translations(translationsxml)
        db = open_database(path, prefix)
        pf = get_profile(profile)
        validate_profile(db, pf)
        tr = get_translation(language)
        new(db, pf, tr, Ptr{Cvoid}[])
    end
end

function close_router(r)
    @ccall LIB.Routino_UnloadDatabase(r.db::DB)::Cvoid
end

function free_xml()
     @ccall LIB.Routino_FreeXMLProfiles()::Cvoid
     @ccall LIB.Routino_FreeXMLTranslations()::Cvoid
end

"""
    clear_waypoints!(r::Router)
Clear the waypoints in the active router, to reuse it for multiple route calculations
"""
function clear_waypoints!(r::Router)
    foreach(free_ptr, r.waypoints)
    empty!(r.waypoints)
end

#==
 routino-router  --quickest --output-text --output-stdout --dir=/trip/osm --prefix=GB --lon1=-1.7587022270945216 --lat1=53.67050050773988 --lon2=-1.6091688317085435 --lat2=53.74083763716689 --profiles=/usr/share/routino/profiles.xml --translations=/usr/share/routino/translations.xml
 ==#
function shell(wayp1::LoLa, wayp2::LoLa, binary=BINARY; datadir=DATADIR, prefix="GB", profiles=PROFILES, translations=TRANSLATIONS)
    try
        read(`$binary --quiet --quickest --output-text --output-stdout --dir=$datadir --prefix=$prefix --lon1=$(wayp1.lo) --lat1=$(wayp1.la) --lon2=$(wayp2.lo) --lat2=$(wayp2.la) --profiles=$profiles --translations=$translations`)
    catch
        ""
    end
end
"""
    route(wayp1::LoLa, wayp2::LoLa, binary; datadir=DATADIR, prefix="GB")
Get the textual output from the shell invocation of the binary
"""
function route(wayp1::LoLa, wayp2::LoLa, binary; datadir=DATADIR, prefix="GB")
    split(String(shell(wayp1, wayp2, binary; datadir, prefix)), "\n")
end

"""
    distance(wayp1::LoLa, wayp2::LoLa, router=Router(), default=nothing)
    distance(wayp1::LoLa, wayp2::LoLa, binary, datadir=DATADIR, prefix="GB", default=nothing; options=1024) # 1025 for shortest
Calculate the route and return the `(km=0, mins=0)` or `default` using either the `Router` or the binary via the shell
"""
function distance(wayp1::LoLa, wayp2::LoLa, router::Router=Router(), default=nothing; options=1024) # quickest
    clear_waypoints!(router)
    if add_waypoint!(wayp1, router) && add_waypoint!(wayp2, router)
        route = calculate_route(router, options)
        km_mins = walk_to_distance(route)
        delete_route(route)
        return km_mins === nothing ? default : km_mins
    end
    default
end

function distance(wayp1::LoLa, wayp2::LoLa, binary::AbstractString, datadir=DATADIR, prefix="GB", default=nothing)
    tline = filter(t->occursin("Waypt#2", t), route(wayp1, wayp2, binary; datadir, prefix))
    if length(tline) == 1
        parts = split(tline[1], "\t")
        km = parse(Float64, split(strip(parts[5]), " ")[1])
        mins = parse(Float64, split(strip(parts[6]), " ")[1])
        (;km, mins)
    else
        default
    end
end
#==
wayp1 =
route(wayp1, wayp2)

wayp1, wayp2 =  (lo=-1.7587022270945216, la=53.67050050773988),  (lo=-1.6091688317085435, la=53.74083763716689) # # "HD2 1JT", "LS27 0AL"

wayp1, wayp2 = (lo=-4.737703404818423, la=58.50695215591557), (lo=-2.954101, la=53.479544)
rs = [Router() for _ in 1:5]
Threads.@threads for i in 1:length(rs)
    println(distance(wayp1, wayp2, rs[i]))
end

==#

struct COutput
    next::Ptr{COutput}
    lon::Cfloat
    lat::Cfloat
    dist::Cfloat
    time::Cfloat
    speed::Cfloat
    type::Cint
    turn::Cint
    bearing::Cint
    name::PUchars
    desc1::PUchars
    desc2::PUchars
    desc3::PUchars
end

check_api_version(ver=ROUTINO_API_VER) = (@ccall LIB.Routino_Check_API_Version(ver::Cint)::Cint) == 0
open_database(path=DATADIR, prefix="GB") = @ccall LIB.Routino_LoadDatabase(path::Cstring, prefix::Cstring)::DB
load_xml_profiles(filename=PROFILES) = (@ccall LIB.Routino_ParseXMLProfiles(filename::Cstring)::Cint) == 0
load_xml_translations(filename=TRANSLATIONS) = (@ccall LIB.Routino_ParseXMLTranslations(filename::Cstring)::Cint) == 0
"""
    get_profile_names()
Get the profile names currently loaded (via the xml)
"""
get_profile_names() = vec_string_from_ptr(@ccall LIB.Routino_GetProfileNames()::ListString)
get_profile(name::AbstractString) = @ccall LIB.Routino_GetProfile(name::Cstring)::Profile
validate_profile(db, profile) = (@ccall LIB.Routino_ValidateProfile(db::DB, profile::Profile)::Cint) == 0
get_translation(name="en") = @ccall LIB.Routino_GetTranslation(name::Cstring)::Translation
find_waypoint(lo, la, db, profile) = @ccall LIB.Routino_FindWaypoint(db::DB, profile::Profile, la::Cdouble, lo::Cdouble)::Waypoint
"""
    find_waypoint(wp, router=Router())
    find_waypoint(lo, la, router=Router())
Test if a waypoint is in the database.
# Arguments
- `wp` waypoint in (lo=0, la=0) format
- `la`, `lo`  Longitude, Latitude
"""
find_waypoint(wp::LoLa, router::Router=Router()) = find_waypoint(wp.lo, wp.la, router)
find_waypoint(lo::Float64, la::Float64, router::Router=Router()) = find_waypoint(lo, la, router.db, router.profile)

calculate_route(waypoints::Vector{Ptr{Cvoid}}, db, profile, translation, options) = @ccall LIB.Routino_CalculateRoute(db::DB, profile::Profile, translation::Translation, waypoints::Waypoint, length(waypoints)::Cint, options::Cint, C_NULL::Ptr{Cvoid})::Ptr{COutput}

function calculate_route(router, options)
    if length(router.waypoints) > 1
        return calculate_route(router.waypoints, router.db, router.profile, router.translation, options)
    end
end

"""
    quickest_route(router)::Ptr{COutput}
    shortest_route(router)::Ptr{COutput}
Run the calculation if there are sufficient waypoints returning either the quickest or shortest route.
"""
quickest_route(router) = calculate_route(router, 1024)
shortest_route(router) = calculate_route(router, 1025)

"""
    add_waypoint!(wp::LoLa, router::Router)::Bool
    add_waypoint!(lo::Float64, la::Float64, router::Router)::Bool
Add a waypoint to the router and return if it was found
# Arguments
- `wp` waypoint in `(lo=0, la=0)` format
- `lo`, `la` Longitude, Latitude
"""
add_waypoint!(wp::LoLa, router::Router) = add_waypoint!(wp.lo, wp.la, router)
function add_waypoint!(lo::Float64, la::Float64, router::Router)
    wp = find_waypoint(lo, la, router)
    if wp != C_NULL
        push!(router.waypoints, wp)
        return true
    end
    false
end

delete_route(ptr::Ptr{COutput}) = @ccall LIB.Routino_DeleteRoute(ptr::Ptr{COutput})::Cvoid

function walk_to_distance(ptr::Ptr{COutput})
    if ptr == C_NULL
        return
    end
    out = unsafe_load(ptr, 1)
    if out.next != C_NULL
        return walk_to_distance(out.next)
    end
    (;km=round(Int, out.dist), mins=round(Int, out.time))
end

#==
    ct9 = (lo=1.3543992881446998, la=51.37956887853703)
    maid = (lo=0.469877, la=51.294277) # 254 # (km = 76, mins = 1137)
    broad = (lo=1.395803, la=51.362011) # 56 # (km = 4, mins = 61)
    rf = Router("foot")
    add_waypoint!(ct9, rf)
=##
###
end


text:
abstract type HTMLNode end

struct HTMLText <: HTMLNode
    parent::HTMLNode
    text::AbstractString
end

# convenience method for defining without parent
HTMLText(text::AbstractString) = HTMLText(NullNode(), text)

struct NullNode <: HTMLNode end

mutable struct HTMLElement{T} <: HTMLNode
    children::Vector{HTMLNode}
    parent::HTMLNode
    attributes::Dict{AbstractString,AbstractString}
end

# convenience method for defining an empty element
HTMLElement(T::Symbol) = HTMLElement{T}(HTMLNode[],NullNode(),Dict{AbstractString,AbstractString}())

mutable struct HTMLDocument
    doctype::AbstractString
    root::HTMLElement
end

struct InvalidHTMLException <: Exception
    msg::AbstractString
end


text:
#
# Functions to work with exact multiplication and solutions for R, for its Cholesky
# (lower) decomposition L, and for L'.
#

using LinearAlgebra
using SparseArrays
using ToeplitzMatrices

"""
    convolve_same(b, kernel)

Convolve the vector `b` with `kernel`, yielding a vector of the same size
as `b`, effectively padding `b` with `length(kernel)-1` initial zeros. This is
equivalent to multiplying `b` by the banded lower triangular Toeplitz matrix
whose first column begins with `kernel`.
"""
function convolve_same(b::AbstractVector, kernel::AbstractVector)
    x = kernel[1] * Array(b)
    for i=2:length(kernel)
        x[i:end] += kernel[i] * b[1:end+1-i]
    end
    x
end

"""
    deconvolve_same(b, kernel)

Reverse the effect of `convolve_same(b, kernel)`. That is, solve `Kx=b` for
`x` where `K` is the banded lower triangular Toeplitz matrix whose first row
begins with `kernel`.
"""
function deconvolve_same(b::AbstractVector, kernel::AbstractVector)
    Nk, N = length(kernel), length(b)
    lenrek = reverse(kernel[2:end])
    x = Array(b) / kernel[1]
    for i=2:min(Nk,N)
        x[i] = (x[i]-dot(kernel[i:-1:2], x[1:i-1])) / kernel[1]
    end
    for i=Nk+1:N
        x[i] = (x[i]-dot(lenrek, x[i-Nk+1:i-1])) / kernel[1]
    end
    x
end

"""
    ARMASolver(m, N)

An object containing the information needed to perform certain noise-related
operations on a vector within an ARMA model `m::ARMAModel`, including the
following (where `R` is the ARMA model noise covariance matrix and `L` is its
lower Cholesky factor):

1. L\\v  `whiten(solver, v)`
1. L*v  `unwhiten(solver, v)`
1. R\\v  `solve_covariance(solver, v)`
1. R*v  `mult_covariance(solver, v)`
1. inv(R[1:n,1:n]) `inverse_covariance(solver, n)`

Note that the object is constructed with a maximum vector length `N`, so that
these operations require `length(v) <= N` and `n <= N`."""
mutable struct ARMASolver
    p         ::Int
    q         ::Int
    phicoef   ::Vector{Float64}
    covarIV   ::Vector{Float64}
    # RR = Phi*R*Phi', a banded matrix = RRt + RRu
    RRu       ::Matrix{Float64}   # The upper left (p)x(p+q) corner of RR-RRt
    RRt       ::Vector{Float64}   # The banded Toeplitz part of RR
    LL        ::SparseMatrixCSC{Float64}   # LL = lower Cholesky factor of RR, also banded

    function ARMASolver(p, q, phicoef, covarIV, RRu, RRt, LL)
        @assert p>=0
        @assert q>=0
        @assert p+1 == length(phicoef)
        @assert length(covarIV) >= 1+q
        @assert length(covarIV) >= p
        @assert phicoef[1] == 1.0

        new(p, q, phicoef, covarIV, RRu, RRt, LL)
    end
end

function ARMASolver(m::ARMAModel, N::Integer)
    covar = model_covariance(m, 1+m.p+2m.q)
    R_corner  = SymmetricToeplitz(covar)
    Nc = length(covar)
    Nbands = max(m.p, m.q+1)
    x = fill(0.0, Nc)
    y = fill(0.0, Nc)
    x[1:m.p+1] = m.phicoef
    y[1] = x[1]
    Phi = Toeplitz(x, y)
    # Hermitian enforces exact symetry even after round-off error.
    RR_corner = Hermitian(Phi*R_corner*Phi')
    RR_toeplitz = fill(0.0, 2m.q+1)
    RR_toeplitz[1:m.q+1] = RR_corner[end, end-m.q:end]
    RR_toeplitz[m.q+1:end] = RR_corner[end, end:-1:end-m.q]
    RR_rectcorner = RR_corner[1:max(m.p,Nbands-1), :]

    # Now compute LL such that LL*LL' == RR, without ever representing RR as a full
    # matrix in memory. Awesome.
    #
    # Find the nontrivial sparse elements of this lower-triangular banded matrix
    # LL, then construct with zero values, then fill them.
    Ndata = N*Nbands - div(Nbands*(Nbands-1),2)
    rows = Array{Int}(undef, Ndata)
    cols = Array{Int}(undef, Ndata)
    vals = fill(0.0, Ndata)
    i = 1
    for r=1:min(N)
        for c=1+max(0,r-Nbands):r
            rows[i] = r
            cols[i] = c
            i += 1
        end
    end
    @assert i-1 == Ndata
    LL = sparse(rows, cols, vals, N, N)

    # Fill LL with its nontrivial values.
    LL_corner = cholesky(RR_corner).L
    for r=1:min(Nc,N)
        for c=1+max(0,r-Nbands):r
            LL[r,c] = LL_corner[r,c]
        end
    end
    for r = Nc+1:N
        mincol = r-m.q
        LL[r,mincol] = RR_toeplitz[1] / LL[mincol,mincol]
        for c=mincol+1:r-1
            LL[r,c] = (RR_toeplitz[c-r+Nbands] - dot(vec(LL[r,mincol:c-1]), vec(LL[c,mincol:c-1]))) / LL[c,c]
        end
        Rt = RR_toeplitz[m.q+1]
        S = sum( LL[ r,mincol:r-1] .^ 2)
        LL[r,r] = sqrt(Rt - S)
    end
    ARMASolver(m.p, m.q, m.phicoef, covar, RR_rectcorner, RR_toeplitz, LL)
end

"""
    whiten(solver::ARMASolver, v::Vector)

Use `solver` for an ARMA model to whiten the vector `v`. Here, "whiten" means
return `w=L\\v`, where `L` is the lower Cholesky factor of the covariance
matrix. If `v` has expectation zero, then  the expected value of `w*w'` over
noise realizations is the identity matrix.

In place of `v`, a matrix can be used. Its columns will be whitened.
"""
function whiten(solver::ARMASolver, v::AbstractVector)
    Phiv = convolve_same(v, solver.phicoef)
    nv = length(v)
    if nv  < size(solver.LL)[1]
        return solver.LL[1:nv, 1:nv] \ Phiv
    end
    solver.LL \ Phiv
end

function whiten(solver::ARMASolver, M::AbstractMatrix)
    ws(v::AbstractVector) = whiten(solver, v)
    mapslices(ws, M, dims=1)
end

"""
    unwhiten(solver::ARMASolver, w::AbstractVector)

Use `solver` for an ARMA model to unwhiten the vector `w`. Here, "unwhiten"
means return `v=L*w`, where `L` is the lower Cholesky factor of the covariance
matrix. If the expected value of `w` is zeros and of `w*w'` is the identity
matrix, then the expected value of `v*v'` is the data covariance matrix of the
ARMA model.

In place of `v`, a matrix can be used. Its columns whill be unwhitened.
"""
function unwhiten(solver::ARMASolver, w::AbstractVector)
    nw = length(w)
    if nw  < size(solver.LL)[1]
        x = solver.LL[1:nw, 1:nw] * w
    else
        x = solver.LL * w
    end
    deconvolve_same(x, solver.phicoef)
end

function unwhiten(solver::ARMASolver, M::AbstractMatrix)
    uws(v::AbstractVector) = unwhiten(solver, v)
    mapslices(uws, M, dims=1)
end

"""
    mult_covariance(solver::ARMASolver, v::AbstractVector)

Use `solver` for an ARMA model to multiply the vector `v` by the covariance matrix.
"""
function mult_covariance(solver::ARMASolver, v::AbstractVector)
    # Reversing before and after deconvolve ensures that we are solving Phi', not Phi.
    v1 = reverse(deconvolve_same(reverse(v), solver.phicoef))

    # Multiply by RR, the "moving-average-only R matrix".
    Nv = length(v1)
    Nu1,Nu2 = size(solver.RRu)
    q = solver.q
    if Nv > Nu2
        x = similar(v1)
        x[1:Nu1] = solver.RRu * v1[1:Nu2]
        for r = Nu1+1 : Nv-q
            x[r] = dot(solver.RRt, v1[r-q:r+q])
        end
        for r = Nv-q+1 : Nv
            x[r] = dot(solver.RRt[1:end-r+Nv-q], v1[r-q:end])
        end

    elseif Nv > Nu1
        x = similar(v1)
        x[1:Nu1] = solver.RRu[:, 1:Nv] * v1
        for r = Nu1+1 : Nv-q
            x[r] = dot(solver.RRt, v1[r-q:r+q])
        end
        for r = max(Nu1,Nv-q)+1 : Nv
            x[r] = dot(solver.RRt[1:end-r+Nv-q], v1[r-q:end])
        end

    else
        x = solver.RRu[1:Nv,1:Nv] * v1
    end
    deconvolve_same(x, solver.phicoef)
end

function mult_covariance(solver::ARMASolver, M::AbstractMatrix)
    mc(v::AbstractVector) = mult_covariance(solver, v)
    mapslices(mc, M, dims=1)
end

"""
    solve_covariance(solver::ARMASolver, v::AbstractVector)

Use `solver` for an ARMA model to solve `Rx=v` where `R` is the noise covariance
matrix.
"""
function solve_covariance(solver::ARMASolver, v::AbstractVector)
    nv = length(v)
    v1 = whiten(solver, v)
    if nv  < size(solver.LL)[1]
        v2 = solver.LL[1:nv,1:nv]' \ v1
    else
        v2 = solver.LL' \ v1
    end
    reverse(convolve_same(reverse(v2), solver.phicoef))
end

function solve_covariance(solver::ARMASolver, M::AbstractMatrix)
    sc(v::AbstractVector) = solve_covariance(solver, v)
    mapslices(sc, M, dims=1)
end

"""
    inverse_covariance(solver::ARMASolver, N)

Use `solver` for an ARMA model to compute the inverse noise covariance matrix,
to size `[N,N]`.
"""
function inverse_covariance(solver::ARMASolver, N::Integer)
    M = Matrix(1.0I, N, N)
    hcat([solve_covariance(solver, M[:,i]) for i=1:N]...)
end


text:
"""
Generate an identifier from a Vulkan identifier, in lower snake case and without pointer prefixes (such as in `pNext`).
"""
function wrap_identifier(identifier)
    identifier ∈ keys(convention_exceptions) && return convention_exceptions[identifier]
    var_str = @match s = string(identifier) begin
        GuardBy(startswith(r"p+[A-Z]")) => remove_prefix(convert(SnakeCaseLower, CamelCaseLower(s))).value
        _ => nc_convert(SnakeCaseLower, s)
    end
    Symbol(var_str)
end

function wrap_identifier(spec::Union{SpecFuncParam,SpecStructMember})

    # handle the case where two identifiers end up being the same
    # issue caused by VkAccelerationStructureBuildGeometryInfoKHR
    # which has both pGeometries and ppGeometries
    siblings = children(parent_spec(spec))
    id = wrap_identifier(spec.name)
    ids = wrap_identifier.(siblings.name)
    if count(==(id), ids) == 2
        lastidx = findlast(==(id), ids)
        if siblings[lastidx] == spec
            Symbol(id, "_2")
        else
            id
        end
    else
        id
    end
end

function wrap_identifier(spec::SpecHandle)
    # try to get an id from an existing function parameter
    cfs = create_funcs(spec)
    id = nothing
    if !isempty(cfs)
        for cf in cfs
            id = _wrap_identifier(spec, first(cfs).func)
            !isnothing(id) && break
        end
    end
    if isnothing(id)
        for f ∈ spec_funcs
            id = _wrap_identifier(spec, f)
            !isnothing(id) && return id
        end
        wrap_identifier(remove_vk_prefix(spec.name))
    end
end

function _wrap_identifier(spec::SpecHandle, func::SpecFunc)
    for param ∈ children(func)
        if spec.name == param.type
            return wrap_identifier(param)
        end
    end
    nothing
end


text:
using PyCall: PyObject

using PyFortran90Namelists: f90nml, @pyinterface

export Tokenizer, update_chars, lex

mutable struct Tokenizer
    o::PyObject
end
Tokenizer() = Tokenizer(f90nml.tokenizer.Tokenizer())

update_chars(tk::Tokenizer) = tk.update_chars()

lex(tk::Tokenizer, line) = tk.parse(line)

@pyinterface Tokenizer


text:
"""
	$(SIGNATURES)

Give tuition as function of [quality, hs gpa, parental income, t].
This is the common portion of tuition. Does not include type specific shocks.
"""
abstract type AbstractTuitionFunction <: ModelObject end
abstract type AbstractTuitionSwitches end

n_colleges(s :: AbstractTuitionSwitches) = s.nColleges;

function validate_tf(tf :: AbstractTuitionFunction)
    return true  # stub +++
end

make_test_tuition_function(switches :: AbstractTuitionSwitches) =
    init_tuition_fct(ObjectId(:test),  switches);

isfree(tf :: AbstractTuitionFunction, qual :: Integer) = false;
# not_free(tf :: AbstractTuitionFunction, qual :: Integer) = !isfree(tf, qual);

free_factor(tf :: AbstractTuitionFunction, qual :: Integer) =
    isfree(tf, qual)  ?  0.0  :  1.0;
free_factor(tf :: AbstractTuitionFunction, qualV :: AbstractVector{I1}) where I1 =
    [free_factor(tf, qual)  for qual ∈ qualV];

# function zero_out_free_colleges(tf :: AbstractTuitionFunction, tuitionV, qual :: Integer)
#     isfree(tf, qual)  &&  (tuitionV[qual] = 0.0);
# end

# function zero_out_free_colleges!(tf :: AbstractTuitionFunction, tuitionV,
#     qualV :: AbstractVector{I1}) where I1
#     for qual in qualV
#         zero_out_free_colleges!(tf, tuitionV, qual);
#     end
# end

include("tuition_by_qual.jl");
include("tuition_linear.jl");

# ---------


text:
using DifferentiableStateSpaceModels, LinearAlgebra, Test, Zygote
using CSV, DataFrames, Turing, DistributionsAD
using DifferentiableStateSpaceModels.Examples
using FiniteDiff: finite_difference_gradient

@testset "Sequence Simulation, 1st order" begin
    m = @include_example_module(Examples.rbc_observables)
    p_f = (ρ=0.2, δ=0.02, σ=0.01, Ω_1=0.1)
    p_d = (α=0.5, β=0.95)

    c = SolverCache(m, Val(1), p_d)
    sol = generate_perturbation(m, p_d, p_f; cache = c)

    T = 9
    eps_value = [[0.22], [0.01], [0.14], [0.03], [0.15], [0.21], [0.22], [0.05], [0.18]]
    x0 = zeros(m.n_x)
    simul = solve(
        DifferentiableStateSpaceModels.dssm_evolution,
        DifferentiableStateSpaceModels.dssm_volatility,
        x0,
        (0, T),
        sol;
        h = DifferentiableStateSpaceModels.dssm_observation,
        noise = eps_value,
    )
    @test simul.z[2:end] ≈ [[-0.0014843113235688628, 0.0], [-0.001672969226342977, -0.013660616212663025], [-0.0025907902434120613, -0.016424018091334355], [-0.002808352075911724, -0.02507880927301923], [-0.003749820417827369, -0.027731843802628737], [-0.00514124775142971, -0.036595970922849726], [-0.006607765317710925, -0.05006822400565087], [-0.006885970365182856, -0.06457803532896965], [-0.00789049635407183, -0.06822941930528102]]
    @test simul.z ≈ solve(sol, x0, (0, T), DifferentiableStateSpaceModels.LTI(); noise = eps_value).z

    # inference
    @inferred solve(
        DifferentiableStateSpaceModels.dssm_evolution,
        DifferentiableStateSpaceModels.dssm_volatility,
        x0,
        (0, T),
        sol;
        h = DifferentiableStateSpaceModels.dssm_observation,
        noise = eps_value,
    )
    # @inferred generate_perturbation(m, p_d, p_f)
end

# @testset "Gradients, generate_perturbation + simulation, 1st order" begin
#     m = @include_example_module(Examples.rbc_observables)
#     p_f = (ρ=0.2, δ=0.02, σ=0.01, Ω_1=0.1)
#     p_d = (α=0.5, β=0.95)
#     p_d_input = [0.5, 0.95]

#     T = 9
#     ϵ_mat = [0.22, 0.01, 0.14, 0.03, 0.15, 0.21, 0.22, 0.05, 0.18]
#     x0 = zeros(m.n_x)

#     function sum_test_joint_first(p_d_input, ϵ_mat, x0, T, p_f, m; kwargs...)
#         p_d = (α=p_d_input[1], β=p_d_input[2])
#         sol = generate_perturbation(m, p_d, p_f; kwargs...)
#         ϵ = map(i -> ϵ_mat[i:i], 1:T)
#         simul = solve(
#             DifferentiableStateSpaceModels.dssm_evolution,
#             DifferentiableStateSpaceModels.dssm_volatility,
#             x0,
#             (0, T),
#             sol;
#             h = DifferentiableStateSpaceModels.dssm_observation,
#             noise = ϵ,
#         )
#         return sum(sum(simul.z))
#     end
#     settings = PerturbationSolverSettings()
#     # The function I wrote here is not type stable
#     # @inferred sum_test_joint_first(p_d_input, ϵ_mat, x0, T, p_f, m; settings)
#     res_zygote = gradient(
#         (p_d_input, ϵ_mat) -> sum_test_joint_first(p_d_input, ϵ_mat, x0, T, p_f, m; settings),
#         p_d_input,
#         ϵ_mat
#     )
#     # p
#     res_finite = finite_difference_gradient(
#        p_d_input -> sum_test_joint_first(p_d_input, ϵ_mat, x0, T, p_f, m; settings),
#        p_d_input,
#     )
#     @test res_zygote[1] ≈ res_finite
#     # ϵ
#     res_finite = finite_difference_gradient(
#         ϵ_mat -> sum_test_joint_first(p_d_input, ϵ_mat, x0, T, p_f, m; settings),
#         ϵ_mat,
#     )
#     @test res_zygote[2] ≈ res_finite

#     # inference
#     # @inferred sum_test_joint_first(p_d_input, ϵ_mat, x0, T, p_f, m; settings)
# end

function kalman_test(p_d_input, p_f, m, z, tspan)
    p_d = (α=p_d_input[1], β=p_d_input[2])
    sol = generate_perturbation(m, p_d, p_f)
    return solve(sol, sol.x_ergodic, tspan; observables = z).logpdf
end

@testset "Kalman filter and its gradient" begin
    m = @include_example_module(Examples.rbc_observables)

    p_f = (ρ=0.2, δ=0.02, σ=0.01, Ω_1=0.1)
    p_d = (α=0.5, β=0.95)
    p_d_input = [0.5, 0.95]
    z = [
        [-0.6949847708598687, -0.8456988740809867],
        [-0.7804117657996692, 0.07781473603479207],
        [-1.1363021614363802, -2.41253450179418],
        [-0.2140813001516194, -0.10914617826240575],
        [-1.0365874981404577, 0.9869373465251516],
        [-0.7321498641416826, 0.012293325072265942],
        [-0.054809260599132194, -1.8233591236618099],
        [0.5407452466493482, -0.9773559802938866],
        [1.3968232347532277, -2.139194998843768],
        [1.3427856576886874, -0.3100476471887863],
    ]

    tspan = (0, length(z))
    sol = generate_perturbation(m, p_d, p_f)
    res = kalman_test(p_d, p_f, m, z, tspan)
    @test res ≈ -805.7558351251781
    res = gradient(p_d_input -> kalman_test(p_d_input, p_f, m, z, tspan), p_d_input)
    res_finite = finite_difference_gradient(p_d_input -> kalman_test(p_d_input, p_f, m, z, tspan), p_d_input)
    @test res_finite ≈ res[1]

    # inference
    @inferred solve(sol, sol.x_ergodic, tspan; observables = z)
    # @inferred kalman_test(p_d_input, p_f, m, z, tspan)
end

function likelihood_test_joint_first(p_d_input, p_f, ϵ, x0, m, tspan, z)
    p_d = (α=p_d_input[1], β=p_d_input[2])
    sol = generate_perturbation(m, p_d, p_f)
    return solve(
        DifferentiableStateSpaceModels.dssm_evolution,
        DifferentiableStateSpaceModels.dssm_volatility,
        x0,
        tspan,
        sol;
        observables = z,
        h = DifferentiableStateSpaceModels.dssm_observation,
        sol.D,
        noise = ϵ
    ).logpdf
end

function likelihood_test_joint_first_sol(p_d_input, p_f, ϵ, x0, m, tspan, z)
    p_d = (α=p_d_input[1], β=p_d_input[2])
    return solve(generate_perturbation(m, p_d, p_f), x0, tspan; observables = z, noise = ϵ).logpdf
end

@testset "Gradients, generate_perturbation + likelihood, 1st order" begin
    m = @include_example_module(Examples.rbc_observables)
    p_f = (ρ=0.2, δ=0.02, σ=0.01, Ω_1=0.1)
    p_d = (α=0.5, β=0.95)
    p_d_input = [0.5, 0.95]
    ϵ = [[0.22], [0.01], [0.14], [0.03], [0.15], [0.21], [0.22], [0.05], [0.18]]
    z = [
        [-0.6949847708598687, -0.8456988740809867],
        [-0.7804117657996692, 0.07781473603479207],
        [-1.1363021614363802, -2.41253450179418],
        [-0.2140813001516194, -0.10914617826240575],
        [-1.0365874981404577, 0.9869373465251516],
        [-0.7321498641416826, 0.012293325072265942],
        [-0.054809260599132194, -1.8233591236618099],
        [0.5407452466493482, -0.9773559802938866],
        [1.3968232347532277, -2.139194998843768],
    ]
    x0 = zeros(m.n_x)
    tspan = (0, length(z))

    res = gradient((p_d_input, ϵ) -> likelihood_test_joint_first(p_d_input, p_f, ϵ, x0, m, tspan, z), p_d_input, ϵ)
    @test res[1] ≈ [303.7133186356109, 553.6149537473261]
    @test res[2] ≈ [
        [40.62454806083384],
        [39.38899479341156],
        [25.095297618483304],
        [26.06697625612332],
        [33.10959536324157],
        [31.484308705831474],
        [19.172319198105615],
        [11.464791870737214],
        [-0.9477420442978448],
    ]

    #lifting the solution object directly.

    res2 = gradient(
        (p_d_input, ϵ) -> likelihood_test_joint_first_sol(p_d_input, p_f, ϵ, x0, m, tspan, z),
        p_d_input,
        ϵ,
    )
    @test res[1] ≈ res2[1]
    @test res[2] ≈ res2[2]

    # inference
    # sol = generate_perturbation(m, p_d, p_f)
    # @inferred solve(
    #     DifferentiableStateSpaceModels.dssm_evolution,
    #     DifferentiableStateSpaceModels.dssm_volatility,
    #     x0,
    #     tspan,
    #     sol;
    #     observables = z,
    #     h = DifferentiableStateSpaceModels.dssm_observation,
    #     sol.D,
    #     noise = ϵ
    # ).logpdf
    # @inferred likelihood_test_joint_first(p_d_input, p_f, ϵ, x0, m, tspan, z)
end

function minimal_likelihood_test_kalman_first(A, B, C, D, u0, noise, observables)
    return solve(A, B, C, D, u0, (0, length(observables)), DifferentiableStateSpaceModels.LTILikelihood(); noise = nothing, observables).logpdf
end

@testset "FVGQ20 Kalman likelhood derivative in 1st order" begin
    path = joinpath(pkgdir(DifferentiableStateSpaceModels), "test", "data")
    file_prefix = "FVGQ20"
    A = Matrix(DataFrame(CSV.File(joinpath(path, "$(file_prefix)_A.csv"), header = false)))
    B = Matrix(DataFrame(CSV.File(joinpath(path, "$(file_prefix)_B.csv"), header = false)))
    C = Matrix(DataFrame(CSV.File(joinpath(path, "$(file_prefix)_C.csv"), header = false)))
    D_raw =
        Matrix(DataFrame(CSV.File(joinpath(path, "$(file_prefix)_D.csv"), header = false)))
    D = Turing.TuringDiagMvNormal(zero(vec(D_raw)), vec(D_raw))
    observables_raw = Matrix(
        DataFrame(
            CSV.File(joinpath(path, "$(file_prefix)_observables.csv"), header = false),
        ),
    )
    noise_raw = Matrix(
        DataFrame(CSV.File(joinpath(path, "$(file_prefix)_noise.csv"), header = false)),
    )
    observables = [observables_raw[i, :] for i = 1:size(observables_raw, 1)]
    noise = [noise_raw[i, :] for i = 1:size(noise_raw, 1)]
    u0_raw = Matrix(
        DataFrame(CSV.File(joinpath(path, "$(file_prefix)_ergodic.csv"), header = false)),
    )
    u0 = DistributionsAD.TuringDenseMvNormal(
        zeros(size(u0_raw, 1)),
        cholesky(Symmetric(u0_raw)),
    )

    minimal_likelihood_test_kalman_first(A, B, C, D, u0, noise, observables)

    res = gradient(minimal_likelihood_test_kalman_first, A, B, C, D, u0, noise, observables)

    # Some tests
    @test finite_difference_gradient(
        A -> minimal_likelihood_test_kalman_first(A, B, C, D, u0, noise, observables),
        A,
    ) ≈ res[1] rtol = 1e-3
    @test finite_difference_gradient(
        B -> minimal_likelihood_test_kalman_first(A, B, C, D, u0, noise, observables),
        B,
    ) ≈ res[2] rtol = 1e-3
    @test finite_difference_gradient(
        C -> minimal_likelihood_test_kalman_first(A, B, C, D, u0, noise, observables),
        C,
    ) ≈ res[3] rtol = 1e-3

    # missing test for the D = MvNormal.  No finite_difference_gradient support yet.

    # @test finite_difference_gradient(u0 -> minimal_likelihood_test_kalman_first(A, B, C, D, u0, noise, observables), u0) ≈ res[5] rtol=1E-7

    observables_grad = finite_difference_gradient(
        observables_mat -> minimal_likelihood_test_kalman_first(
            A,
            B,
            C,
            D,
            u0,
            noise,
            [observables_mat[i, :] for i = 1:size(observables_mat, 1)],
        ),
        observables_raw,
    )
    @test [observables_grad[i, :] for i = 1:size(observables_raw, 1)] ≈ res[7] rtol = 1e-5

    # inference
    @inferred solve(
        A,
        B,
        C,
        D,
        u0,
        (0, length(observables)),
        DifferentiableStateSpaceModels.LTILikelihood();
        noise = nothing,
        observables,
    )
end

function minimal_likelihood_test_joint_first(A, B, C, D, u0, noise, observables)
    return solve(A, B, C, D, u0, (0, length(observables)), DifferentiableStateSpaceModels.LTILikelihood(); noise, observables).logpdf
end

@testset "FVGQ20 joint likelhood derivative in 1st order" begin
    path = joinpath(pkgdir(DifferentiableStateSpaceModels), "test", "data")
    file_prefix = "FVGQ20"
    A = Matrix(DataFrame(CSV.File(joinpath(path, "$(file_prefix)_A.csv"), header = false)))
    B = Matrix(DataFrame(CSV.File(joinpath(path, "$(file_prefix)_B.csv"), header = false)))
    C = Matrix(DataFrame(CSV.File(joinpath(path, "$(file_prefix)_C.csv"), header = false)))
    D_raw =
        Matrix(DataFrame(CSV.File(joinpath(path, "$(file_prefix)_D.csv"), header = false)))
    D = MvNormal(Diagonal(map(abs2, vec(D_raw))))
    observables_raw = Matrix(
        DataFrame(
            CSV.File(joinpath(path, "$(file_prefix)_observables.csv"), header = false),
        ),
    )
    noise_raw = Matrix(
        DataFrame(CSV.File(joinpath(path, "$(file_prefix)_noise.csv"), header = false)),
    )
    observables = [observables_raw[i, :] for i = 1:size(observables_raw, 1)]
    noise = [noise_raw[i, :] for i = 1:size(noise_raw, 1)]
    u0 = zeros(size(A, 1))

    res = gradient(minimal_likelihood_test_joint_first, A, B, C, D, u0, noise, observables)

    # Some tests
    @test finite_difference_gradient(
        A -> minimal_likelihood_test_joint_first(A, B, C, D, u0, noise, observables),
        A,
    ) ≈ res[1] rtol = 1e-5
    @test finite_difference_gradient(
        B -> minimal_likelihood_test_joint_first(A, B, C, D, u0, noise, observables),
        B,
    ) ≈ res[2] rtol = 1e-5
    @test finite_difference_gradient(
        C -> minimal_likelihood_test_joint_first(A, B, C, D, u0, noise, observables),
        C,
    ) ≈ res[3] rtol = 1e-5

    # missing test for the D = MvNormal.  No finite_difference_gradient support yet.

    @test finite_difference_gradient(
        u0 -> minimal_likelihood_test_joint_first(A, B, C, D, u0, noise, observables),
        u0,
    ) ≈ res[5] rtol = 1e-7

    noise_grad = finite_difference_gradient(
        noise_mat -> minimal_likelihood_test_joint_first(
            A,
            B,
            C,
            D,
            u0,
            [noise_mat[i, :] for i = 1:size(noise_mat, 1)],
            observables,
        ),
        noise_raw,
    )
    @test [noise_grad[i, :] for i = 1:size(noise_raw, 1)] ≈ res[6] rtol = 1E-7

    observables_grad = finite_difference_gradient(
        observables_mat -> minimal_likelihood_test_joint_first(
            A,
            B,
            C,
            D,
            u0,
            noise,
            [observables_mat[i, :] for i = 1:size(observables_mat, 1)],
        ),
        observables_raw,
    )
    @test [observables_grad[i, :] for i = 1:size(observables_raw, 1)] ≈ res[7] rtol = 1E-7

    # inference
    @inferred solve(
        A,
        B,
        C,
        D,
        u0,
        (0, length(observables)),
        DifferentiableStateSpaceModels.LTILikelihood();
        noise,
        observables,
    )
end


text:
using Gadfly, RDatasets, Dates

set_default_plot_size(6inch, 3inch)

economics = dataset("HistData", "Prostitutes")
p = plot(economics, x=:Date, y=:Count, Geom.line)


text:
export GradientBanditLearner

using Flux: softmax, onehot

"""
    GradientBanditLearner(;approximator::A, optimizer::O, baseline::B)
"""
Base.@kwdef mutable struct GradientBanditLearner{A,O,B} <: AbstractLearner
    approximator::A
    optimizer::O
    baseline::B
end

(learner::GradientBanditLearner)(obs::Observation) =
    obs |> get_state |> learner.approximator |> softmax

function update!(learner::GradientBanditLearner, transitions)
    s, a, r = transitions
    probs = s |> learner.approximator |> softmax
    r̄ = learner.baseline isa Number ? learner.baseline : learner.baseline(r)
    errors = (r - r̄) .* (onehot(a, 1:length(probs)) .- probs)
    update!(learner.approximator, s => apply!(learner.optimizer, s, errors))
end

function extract_transitions(buffer::EpisodeTurnBuffer, ::GradientBanditLearner)
    if length(buffer) > 0
        state(buffer)[end-1], action(buffer)[end-1], reward(buffer)[end]
    else
        nothing
    end
end


text:
export fmtstar!

fmtstar!(P::MPProblem; kwargs...) = fmtstar!(P, length(P.V); kwargs...)
function fmtstar!{T}(P::MPProblem{T}, N::Int; rm::T = T(1),
                                              connections::Symbol = :R,
                                              k = min(ceil(Int, (2*rm)^dim(P.SS)*(e/dim(P.SS))*log(N)), N-1),
                                              r = T(0),
                                              ensure_goal_ct = 1,
                                              init_idx = 1,
                                              checkpts = true)  # TODO: bleh, prefer false
    tic()
    P.CC.count = 0

    if connections == :R
        nearF = inballF!
        nearB = inballB!
    elseif connections == :K
        nearF = mutualknnF!
        nearB = knnB!
    else
        error("Connection type must be radial (:R) or k-nearest (:K)")
    end
    r > 0 && setup_steering(P.SS, r)
    if !is_free_state(P.init, P.CC, P.SS)
        warn("Initial state is infeasible!")
        P.status = :failed
        P.solution = MPSolution(P.status, T(Inf), toq(), Dict())
        return T(Inf)
    end
    free_volume_ub = sample_free!(P, N - length(P.V), ensure_goal_ct = ensure_goal_ct)  # TODO: clean this logic up
    if checkpts
        F = trues(N)
        for i in 1:N
            F[i] = is_free_state(P.V[i], P.CC, P.SS)
        end
    end
    if r == 0
        d = dim(P.SS)
        r = T(rm*2*(1/d*free_volume_ub/(pi^(d/2)/gamma(d/2+1))*log(N)/N)^(1/d))
        setup_steering(P.SS, r)
    end

    A = zeros(Int,N)
    W = trues(N)
    H = falses(N)
    C = zeros(T,N)
    P.V.init = P.init
    if P.V[init_idx] == P.init
        W[init_idx] = false
        H[init_idx] = true
        HHeap = Collections.PriorityQueue([init_idx], T[0])
    else    # special casing the first expansion round of FMT if P.init is not in the sample set
        # HHeap = Collections.PriorityQueue(Int[], T[])
        # neighborhood = (connections == :R ? inballF(P.V, P.init, r) : knnF(P.V, P.init, r))
        # for ii in 1:length(nonzeroinds(neighborhood))
        #     x, c = nonzeroinds(neighborhood)[ii], nonzeros(neighborhood)[ii]
        #     if is_free_motion(P.init, P.V[x], P.CC, P.SS)
        #         A[x] = 0
        #         C[x] = c
        #         HHeap[x] = c
        #         H[x] = true
        #         W[x] = false
        #     end
        # end
    end
    z = Collections.dequeue!(HHeap)    # i.e. z = init_idx

    while !is_goal_pt(P.V[z], P.goal, P.SS)
        H_new = Int[]
        for x in (connections == :R ? nonzeroinds(nearF(P.V, z, r, W)) : nonzeroinds(nearF(P.V, z, k, W)))
            checkpts && !F[x] && continue
            neighborhood = (connections == :R ? nearB(P.V, x, r, H) : nearB(P.V, x, k, H))
            c_min, y_idx = findmin(C[nonzeroinds(neighborhood)] + nonzeros(neighborhood))
            y_min = nonzeroinds(neighborhood)[y_idx]
            if is_free_motion(P.V[y_min], P.V[x], P.CC, P.SS)
                A[x] = y_min
                C[x] = c_min
                HHeap[x] = c_min
                push!(H_new, x)
                W[x] = false
            end
        end
        H[H_new] = true
        H[z] = false
        if !isempty(HHeap)
            z = Collections.dequeue!(HHeap)
        else
            break
        end
    end

    sol = [z]
    costs = [C[z]]
    while sol[1] != 1
        unshift!(sol, A[sol[1]])
        if sol[1] == 0
            unshift!(costs, T(0))
            break
        end
        unshift!(costs, C[sol[1]])
    end

    P.status = is_goal_pt(P.V[z], P.goal, P.SS) ? :solved : :failed
    solution_metadata = Dict(
        "radius_multiplier" => rm,
        "collision_checks" => P.CC.count,
        "num_samples" => N,
        "cost" => C[z],
        "cumcost" => costs,
        "planner" => "FMTstar",
        "solved" => is_goal_pt(P.V[z], P.goal, P.SS),
        "tree" => A,
        "path" => sol
    )
    connections == :R && (solution_metadata["r"] = r)
    connections == :K && (solution_metadata["k"] = k)
    P.solution = MPSolution(P.status, C[z], toq(), solution_metadata)
    P.status, P.solution.cost, P.solution.elapsed
end


text:
module NeoVim

using Compat
using Docile
import MsgPack
export Nvim, NeoVimError, request, respond

immutable Val{T} end

include("nvim.jl")
include("buffer.jl")
include("window.jl")
include("tabpage.jl")
include("api.jl")

function __init__()
    n, p = Nvim(Val{:Spawn})
    global const API = api_info(n)
    close(n)
    # TODO: do this by sending :q! if possible
    kill(p)

    declare_err(API["error_types"])
    declare_type(API["types"])
    declare_func(API["functions"])
end

immutable NeoVimError <: Exception
    id::Int
    msg::UTF8String
end
@compat NeoVimError(errid, msg::Vector{UInt8}) = NeoVimError(errid, bytestring(msg))
NeoVimError(a::Vector) = NeoVimError(a[1], a[2])

const REQUEST = 0
const RESPONSE = 1
const NOTIFICATION = 2

function request(n::Nvim, func::ByteString, args...)
    requestid = n.reqid
    n.reqid += 1

    MsgPack.pack(n, Any[REQUEST, requestid, func, collect(args)])

    err, res = MsgPack.unpack(n)[3:4]
    if err !== nothing
        throw(NeoVimError(err))
    else
        return sanitize(n, res)
    end
end
function respond(n::Nvim, reqid, args...)
    MsgPack.pack(n, Any[RESPONSE, reqid, nothing, collect(args)])
end
Base.error(n::Nvim, reqid, msg) = MsgPack.pack(n, Any[RESPONSE, reqid, msg, nothing])

function eventloop(nvim::Nvim, data)
    @async while true
        msg = MsgPack.unpack(nvim)
        msgtype = msg[1]
        if msgtype == REQUEST
            reqid = msg[2]
            event = symbol(msg[3])
            args = sanitize(nvim, msg[4])
            # TODO: get logging sorted out
            try
                request_callback(Val{event}, nvim, reqid, args, data)
            catch err
                error(nvim, reqid, string(err))
            end
        elseif msgtype == RESPONSE
            # TODO: make requests asynchronous
        elseif msgtype == NOTIFICATION
            event = symbol(msg[2])
            args = sanitize(nvim, msg[3])
            try
                notification_callback(Val{event}, nvim, args, data)
            end
        else
            throw(NeoVimError(-1, "Unknown message type."))
        end
    end
end

function request_callback(::Any, nvim, reqid, args, data)
    error(nvim, reqid, "Unable to handle event.")
end

function notification_callback(::Any, nvim, args, data)
    throw(NeoVimError(-1, "You must define a notification_callback method."))
end

end # module


text:
julia> Pkg.add("CellularAutomata")
INFO: Installing CellularAutomata v0.1.2
INFO: Package database updated

julia> using CellularAutomata

julia> gameOfLife{T<:Int}(n::T, m::T, gen::T) = CA2d([3], [2,3], int(randbool(n, m)), gen)
gameOfLife (generic function with 1 method)

julia> gameOfLife(15, 30, 5)
30x15x5 Cellular Automaton

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
import ..TimePetriNets
using SparseArrays

function TimePetriNets.TPN(n::Neuron{H}, name=:neuron; M=Int) where {H}
    # collect all places and transitions to glue dendrite segments + synapses together
    P = Symbol[]
    T = Symbol[]
    m₀ = M[]
    eft = H[]
    lft = H[]
    trans_probs = Float64[]

    R = Tuple{Int,Int,M}[]
    ΔF = Tuple{Int,Int,M}[]

    function add_place!(name, initial)
        push!(P, Symbol(name))
        push!(m₀, M(initial))
        length(P)
    end

    function add_transition!(name, e, l=e, p=1.0)
        push!(T, Symbol(name))
        push!(eft, H(e))
        push!(lft, H(l))
        push!(trans_probs, p)
        length(T)
    end

    function copy_transitions!(other, names...)
        idxs = indexin(Symbol.(names), other.T)
        return add_transition!.(
            names,
            other.eft[idxs],
            other.lft[idxs],
            other.trans_probs[idxs]
        )
    end

    function copy_places!(other, names...)
        idxs = indexin(Symbol.(names), other.P)
        return add_place!.(
            names,
            other.x₀.m[idxs]
        )
    end

    sub_tpns = TPN{M,H}[]

    # go through dendritic tree
    queue = Tuple{SomaOrDendrite,Union{Nothing,Int},Union{Nothing,Int}}[(n.dendrite, nothing, nothing)]
    while !isempty(queue)
        (seg, parent_counter, parent_trans) = pop!(queue)
        has_reset = !isempty(seg.inh_inputs)

        # the dendrite is mainly a monoflop
        dendrite_tpn = TimePetriNets.Monoflop(Symbol("$(name)_$(seg.name)"),
            seg==n.dendrite ? seg.parameters.spike_duration+seg.parameters.refractory_period : seg.parameters.plateau_duration; add_reset=has_reset)
        push!(sub_tpns, dendrite_tpn)

        # need to access start, stop transition as well as the on place --> copy them
        t_names = ("$(name)_$(seg.name)_start", "$(name)_$(seg.name)_stop")
        seg_start, seg_stop = copy_transitions!(dendrite_tpn, t_names...)
        seg_reset = has_reset ? copy_transitions!(dendrite_tpn, "$(name)_$(seg.name)_reset")[1] : nothing
        (seg_on,) = copy_places!(dendrite_tpn, "$(name)_$(seg.name)_on")

        # depending on excitatory #inputs, add synaptic counter
        syn_counter_exc = length(seg.exc_inputs) == 1 ? nothing :
                          add_place!("$(name)_$(seg.name)_syn_counter_exc", 0)

        # depending on inhibitory #inputs, add synaptic counter and inverted counter
        syn_counter_inh = length(seg.inh_inputs) == 1 ? nothing :
                          add_place!("$(name)_$(seg.name)_syn_counter_inh", 0)
        syn_counter_inh_off = add_place!("$(name)_$(seg.name)_syn_counter_inh_off", length(seg.inh_inputs))

        # depending on #children, add dendritic counter
        den_counter = length(seg.children) == 1 ? nothing :
                      add_place!("$(name)_$(seg.name)_den_counter", 0)

        # if there is a parent counter, update it, too
        if !isnothing(parent_counter)
            push!(ΔF, (parent_counter, seg_start, 1), (parent_counter, seg_stop, -1))
        end

        # if there is a parent transaction, make it sensitive to this plateau
        if !isnothing(parent_trans)
            push!(R, (seg_on, parent_trans, 1))
        end

        ## Add all inhibitory synapses
        for inp in seg.inh_inputs
            # the inhibitory synapse is mainly a monoflop
            synapse_tpn = TimePetriNets.Monoflop(Symbol("$(name)_$(seg.name)_$(inp)"), seg.parameters.ipsp_duration; add_reset=false)
            push!(sub_tpns, synapse_tpn)
            # we need to access the start and stop transition and on place --> copy them
            (t_start, t_stop) = copy_transitions!(synapse_tpn, "$(name)_$(seg.name)_$(inp)_start", "$(name)_$(seg.name)_$(inp)_stop")
            (p_on,) = copy_places!(synapse_tpn, "$(name)_$(seg.name)_$(inp)_on")

            # Add trigger place
            t_trigger = add_place!("$(name)_$(seg.name)_$(inp)_trigger", 0)
            push!(ΔF, (t_trigger, t_start, -1))

            # if there is a synapse counter, update it, too
            if !isnothing(syn_counter_inh)
                push!(ΔF, (syn_counter_inh, t_start, 1), (syn_counter_inh, t_stop, -1))
            else
                syn_counter_inh = p_on
            end
            # off-counter works in opposite direction
            push!(ΔF, (syn_counter_inh_off, t_start, -1), (syn_counter_inh_off, t_stop, 1))
        end

        if has_reset
            if !isnothing(parent_counter)
                push!(ΔF, (parent_counter, seg_reset, -1))
            end

            # enable reset when at least one inhibitory synapse is on
            # enable segment if all inhibitory synapses are off
            push!(R,
                (syn_counter_inh, seg_reset, 1),
                (syn_counter_inh_off, seg_start, length(seg.inh_inputs))
            )
        end

        ## Add all excitatory synapses
        for inp in seg.exc_inputs
            # the excitatory synapse is mainly a monoflop
            synapse_tpn = TimePetriNets.Monoflop(Symbol("$(name)_$(seg.name)_$(inp)"), seg.parameters.epsp_duration; add_reset=false)
            push!(sub_tpns, synapse_tpn)
            # we need to access the start and stop transition and on place --> copy them
            (t_start, t_stop) = copy_transitions!(synapse_tpn, "$(name)_$(seg.name)_$(inp)_start", "$(name)_$(seg.name)_$(inp)_stop")
            (p_on,) = copy_places!(synapse_tpn, "$(name)_$(seg.name)_$(inp)_on")

            # Add trigger place
            t_trigger = add_place!("$(name)_$(seg.name)_$(inp)_trigger", 0)
            push!(ΔF, (t_trigger, t_start, -1))

            # if there is a synapse counter, update it, too
            if !isnothing(syn_counter_exc)
                push!(ΔF, (syn_counter_exc, t_start, 1), (syn_counter_exc, t_stop, -1))
            else
                syn_counter_exc = p_on
            end
        end

        # Make start transition sensitive to excitatory synaptic input
        push!(R, (syn_counter_exc, seg_start, seg.parameters.synaptic_threshold))

        if !isnothing(den_counter)
            # Make start transition sensitive to dendritic input
            push!(R, (den_counter, seg_start, min(length(seg.children), seg.parameters.dendritic_threshold)))

            # Go through child branches with shared dendrite counter
            for child in seg.children
                push!(queue, (child, den_counter, nothing))
            end
        else
            # Add single child branch without extra counter -> directly connect to seg_start
            push!(queue, (seg.children[1], nothing, seg_start))
        end
    end

    m = length(P)
    n = length(T)

    R = isempty(R) ? spzeros(M,m,n) : sparse(collect.(zip(R...))..., m, n)
    ΔF = isempty(ΔF) ? spzeros(M,m,n) : sparse(collect.(zip(ΔF...))..., m, n)

    glue = TimePetriNets.TPN(P, T, R, ΔF, eft, lft, trans_probs, m₀)
    |(glue, sub_tpns...)
end

function TimePetriNets.TPN(net::NeuralNetwork, args...;kwargs...)
    # construct all neurons' tpns
    neuron_tpns = [TimePetriNets.TPN(neuron,name,args...;kwargs...) for (name,neuron) in pairs(net.neurons)]
    (H,M) = eltype(neuron_tpns).parameters

    input_trans = Monoflop.(Symbol.(Ref("input_").*String.(net.inputs)), 0)

    # combine synapse and neuron tpns into one tpn
    net_tpn = |(input_trans..., neuron_tpns...)

    # set lft of all inputs to infinity
    for input in net.inputs
        idx = net_tpn.T_index[Symbol("input_$(input)_start")]
        net_tpn.eft[idx] = zero(H)
        net_tpn.lft[idx] = typemax(H)
    end

    function add_axon!(src_idx, axons)
        for axon in axons
            if isa(axon.target,Symbol)
                continue
            end

            (neuron,input) = axon.target
            neuron_name = neuron[].name[]

            for dendrite in neuron[].ports[input]
                tgt_idx = net_tpn.P_index[Symbol("$(neuron_name)_$(dendrite.name)_$(input)_trigger")]
                # add arc from transition to trigger place
                net_tpn.ΔF[tgt_idx, src_idx] += 1

                t_idx = get(net_tpn.T_index, Symbol("$(neuron_name)_$(dendrite.name)_$(input)_start"), nothing)
                if !isnothing(t_idx)
                    net_tpn.eft[t_idx] = axon.parameters.axon_delay
                    net_tpn.lft[t_idx] = axon.parameters.axon_delay
                end
                t_idx = get(net_tpn.T_index, Symbol("$(neuron_name)_$(dendrite.name)_$(input)_fail"), nothing)
                if !isnothing(t_idx)
                    net_tpn.eft[t_idx] = axon.parameters.axon_delay
                    net_tpn.lft[t_idx] = axon.parameters.axon_delay
                end
            end
        end
        nothing
    end

    # add synapses between the neurons
    for (source_name, axons) in pairs(net.axons)
        # find the spike-start transition of the source neuron or input
        source_name = if source_name ∈ net.inputs
            Symbol("input_$(source_name)")
        else
            Symbol("$(source_name)_$(net.neurons[source_name].dendrite.name)")
        end
        src_name = Symbol("$(source_name)_start")
        src_idx = net_tpn.T_index[src_name]
        add_axon!(src_idx, axons)
    end
    # update TPN
    TimePetriNets.update!(net_tpn)
    return net_tpn
end


text:
julia> collection = []
0-element Array{Any,1}

julia> push!(collection, 1,2,4,7)
4-element Array{Any,1}:
 1
 2
 4
 7


text:
using DifferentialEquations
using Plots
plotly()
# Plot a Van der Pol oscillator: \frac{d^2 x}{dt^2} - \mu(1-x^2) \frac{dx}{dt} + x = 0

function f!(dx, x, μ, t) # Can either be in-place f, in which case the first argument is dx, or out of place, in which case the first argument is x and the output is dx
    dx[1] = x[2]
    dx[2] = μ(t)*(1-x[1]^2).*x[2] - x[1]
end

X0 = [1.0, 1.0]
tspan = (0.0, 2000.0)
μ = unitStep(1000.0).*20.0 # A step from 0 to 20 at t = 1000
prob = ODEProblem(f!, X0, tspan, μ)
sol = solve(prob; tstops=sort(collect(μ.d)), saveat=0.001, reltol=1e-6)

plot(sol, vars=(1))
plot!(title="Van der Pol, parameter steps from 0 to 20")


text:
```
currentAsset: a Liquid Asset

    E.g: Cash , accountRecievable (Money from a Debitor), ...
```

mutable struct currentAsset <: Asset
    #fields
end

max_stars_repo_name: adamwillisXanax/accountingJulia

max_issues_repo_name: adamwillisXanax/accountingJulia

max_forks_repo_name: adamwillisXanax/accountingJulia


text:
using Distributed
using Distributions

@everywhere using Pkg
@everywhere Pkg.activate(".")

@everywhere using Psychometrics

using Pkg
Pkg.activate(".")
using Psychometrics

function est_pg()
    I_total = 50
    N = 3000

    # ITEM PARAMETERS AND LATENTS

    items = [
        Item(
            i,
            string("item_", i),
            ["math"],
            Parameters2PL(
                Product([LogNormal(0.0, 0.25), Normal(0, 1)]),
                [1e-5, Inf],
                [-Inf, Inf],
            ),
        ) for i = 1:I_total
    ]
    examinees = [
        Examinee(e, string("examinee_", e), Latent1D(Normal(0.0, 1.0), [-6.0, 6.0])) for
        e = 1:N
    ]
    # RESPONSES
    responses = answer(examinees, items)

    #Initial Values, need to make sure that all Variance needs to be positive
    items_est = [
        Item(
            i,
            string("item_", i),
            ["math"],
            Parameters2PL(
                Product([TruncatedNormal(1, 5, 0.0, Inf), Normal(0, 5)]),
                [1e-5, 5.0],
                [-6.0, 6.0],
            ),
        ) for i = 1:I_total
    ]
    map(i -> i.parameters.a = 1.0, items_est)
    map(i -> i.parameters.b = 0.0, items_est)
    map(i -> empty_chain!(i), items_est)

    examinees_est = [
        Examinee(e, string("examinee_", e), Latent1D(Normal(0, 1), [-6.0, 6.0])) for
        e = 1:N
    ]
    for n = 1:N
        # examinees_est[n].latent.val = examinees[n].latent.val + rand(Normal(0.0,0.3))
        examinees_est[n].latent.val = 0.0
        #examinees_est[n].latent.prior = Normal(examinees[n].latent.val, 0.3)
    end
    map(i -> empty_chain!(i), examinees_est)

    joint_estimate!(
        items_est,
        examinees_est,
        responses;
        method = "pg",
        mcmc_iter = 4000,
        max_time = 400,
        item_sampling = false,
        examinee_sampling = false,
        verbosity = 2, #default
    )

    println("a RMSE")
    println(
        sqrt(
            mean(
                map(
                    (i, i_est) -> (i.parameters.a - i_est.parameters.a)^2,
                    items,
                    items_est,
                ),
            ),
        ),
    )
    println("b RMSE")
    println(
        sqrt(
            mean(
                map(
                    (i, i_est) -> (i.parameters.b - i_est.parameters.b)^2,
                    items,
                    items_est,
                ),
            ),
        ),
    )
    println("latent RMSE")
    println(
        sqrt(
            mean(
                map(
                    (i, i_est) -> (i.latent.val - i_est.latent.val)^2,
                    examinees,
                    examinees_est,
                ),
            ),
        ),
    )
    return examinees, examinees_est, items, items_est, responses
end
examinees, examinees_est, items, items_est, responses = est_pg();


text:
using Flux
using Flux: onehotbatch, onecold, logitcrossentropy
using Flux: DataLoader
using MLDatasets, CUDA
import Dates
import BSON

Base.@kwdef struct Args
    η::Float32 = 3e-4
    batchsize::Integer = 128
    epochs::Integer = 50
    infotime::Integer = 5
    usecuda::Bool = true
    savemodel::Bool = true
    savepath::String = "runs/"
end

function dataloaders(args)
    ENV["DATADEPS_ALWAYS_ACCEPT"] = "true"

    xtrain, ytrain = MNIST.traindata(Float32)
    xtest, ytest = MNIST.testdata(Float32)

    xtrain, xtest = flatten(xtrain), flatten(xtest)

    ytrain = onehotbatch(ytrain, sort(unique(ytrain)))
    ytest = onehotbatch(ytest, sort(unique(ytest)))

    trainloader = DataLoader((xtrain, ytrain); batchsize=args.batchsize, shuffle=true)
    testloader = DataLoader((xtest, ytest); batchsize=args.batchsize)

    trainloader, testloader
end

function eval_loss_accuracy(loader, model, device)
    acc = Float32(0)
    ls = Float32(0)
    ntot = UInt64(0)

    for (x′, y′) ∈ loader
        x, y = device(x′), device(y′)
        ŷ = model(x)
        ls += logitcrossentropy(ŷ, y, agg=sum)
        acc += sum(onecold(ŷ) .== onecold(y))
        ntot += size(x, ndims(x))
    end

    round(ls/ntot, digits=4), round(100acc/ntot, digits=4)
end

function train(; kws...)
    args = Args(; kws...)

    device = (args.usecuda && CUDA.functional()) ? gpu : cpu

    model = Chain(Dense(28*28, 32, relu), Dense(32, 10), softmax) |> device

    ps = params(model)

    trainloader, testloader = dataloaders(args)

    opt = ADAM(args.η)

    function report(epoch)
        train = eval_loss_accuracy(trainloader, model, device)
        test = eval_loss_accuracy(testloader, model, device)
        @info "Epoch $epoch\n\t• Train loss/acc: $train\n\t• Test loss/acc: $test"
    end

    @info "$(args.epochs) training epochs beginning now."
    report(0)
    for epoch ∈ 1:args.epochs
        for (x′, y′) ∈ trainloader
            x, y = device(x′), device(y′)
            gs = gradient(ps) do
                logitcrossentropy(model(x), y)
            end
            Flux.Optimise.update!(opt, ps, gs)
        end
        epoch % args.infotime == 0 && report(epoch)
    end

    # Saving model
    if args.savemodel
        !ispath(args.savepath) && mkpath(args.savepath)
        modelpath = joinpath(args.savepath, "model_$(Dates.today()).bson")
        let m = cpu(model)
            BSON.@save modelpath m
        end
        @info "Model saved in $modelpath."
    end
    return cpu(model)
end

if abspath(PROGRAM_FILE) == @__FILE__
    cd(@__DIR__)
    model = @time train(; epochs=100)
end

max_stars_repo_name: benharmonics/ModelZoo

max_issues_repo_name: benharmonics/ModelZoo

max_forks_repo_name: benharmonics/ModelZoo