text:
get_test_curve1() = [CurvePt(VecSE2(0.0,0.0,0.0), 0.0),
                     CurvePt(VecSE2(1.0,0.0,0.0), 1.0),
                     CurvePt(VecSE2(3.0,0.0,0.0), 3.0)]

let
    p = lerp(CurvePt(VecSE2(0.0,0.0,0.0), 0.0), CurvePt(VecSE2(1.0,2.0,3.0), 4.0), 0.25)
    show(IOBuffer(), p)
    @test isapprox(convert(VecE2, p.pos), VecE2(0.25, 0.5))
    @test isapprox(p.pos.θ, 0.75)
    @test isapprox(p.s, 1.0)

    show(IOBuffer(), CurveIndex(1,0.0))

    @test isapprox(get_lerp_time(VecE2(0.0,0.0), VecE2(1.0,0.0), VecE2(1.0,0.0)), 1.0)
    @test isapprox(get_lerp_time(VecE2(0.0,0.0), VecE2(1.0,0.0), VecE2(0.5,0.0)), 0.5)
    @test isapprox(get_lerp_time(VecE2(0.0,0.0), VecE2(1.0,0.0), VecE2(0.5,0.5)), 0.5)
    @test isapprox(get_lerp_time(VecE2(0.0,0.0), VecE2(1.0,1.0), VecE2(0.5,0.5)), 0.5)
    @test isapprox(get_lerp_time(VecE2(1.0,0.0), VecE2(2.0,0.0), VecE2(1.5,0.5)), 0.5)
    @test isapprox(get_lerp_time(VecE2(0.0,0.0), VecE2(-1.0,0.0), VecE2(1.0,0.0)), 0.0)
    @test isapprox(get_lerp_time(VecE2(0.0,0.0), VecE2(-1.0,0.0), VecE2(-0.75,0.0)), 0.75)

    curve = get_test_curve1()
    @test index_closest_to_point(curve, VecE2(0.0,0.0)) == 1
    @test index_closest_to_point(curve, VecE2(1.0,0.0)) == 2
    @test index_closest_to_point(curve, VecE2(2.1,0.0)) == 3
    @test index_closest_to_point(curve, VecE2(0.49,0.0)) == 1
    @test index_closest_to_point(curve, VecE2(1.9,-100.0)) == 2
    @test index_closest_to_point(curve, VecSE2(1.9,-100.0,0.0)) == 2
    @test index_closest_to_point(curve, VecSE2(-1.0,0.0,0.0)) == 1
    @test index_closest_to_point([CurvePt(VecSE2(0.0,0.0,0.0),0.0)], VecE2(0.0,0.0)) == 1

    @test get_curve_index(curve, 0.0) == CurveIndex(1, 0.0)
    @test get_curve_index(curve, 1.0) == CurveIndex(2, 0.0)
    @test get_curve_index(curve, 1.5) == CurveIndex(2, 0.25)
    @test get_curve_index(curve, 3.5) == CurveIndex(2, 1.0)
    @test get_curve_index(curve,-0.5) == CurveIndex(1, 0.0)

    p = curve[CurveIndex(1, 0.75)]
    @test isapprox(convert(VecE2, p.pos), VecE2(0.75, 0.0))
    @test isapprox(p.pos.θ, 0.0)
    @test isapprox(p.s, 0.75)

    @test get_curve_index(CurveIndex(1, 0.0), curve,  0.25) == CurveIndex(1,0.25)
    @test get_curve_index(CurveIndex(1, 0.0), curve,  1.25) == CurveIndex(2,0.125)
    @test get_curve_index(CurveIndex(1, 0.5), curve,  1.50) == CurveIndex(2,0.5)
    @test get_curve_index(CurveIndex(1, 0.5), curve, -0.25) == CurveIndex(1,0.25)
    @test get_curve_index(CurveIndex(2, 0.0), curve, -0.50) == CurveIndex(1,0.50)

    res = proj(VecSE2(0.0,0.0,0.0), curve)
    show(IOBuffer(), res)
    @test res.ind == CurveIndex(1, 0.0)
    @test isapprox(res.t, 0.0)
    @test isapprox(res.ϕ, 0.0)
    res = proj(VecSE2(0.25,0.5,0.1), curve)
    @test res.ind == CurveIndex(1, 0.25)
    @test isapprox(res.t, 0.5)
    @test isapprox(res.ϕ, 0.1)
    res = proj(VecSE2(0.25,-0.5,-0.1), curve)
    @test res.ind == CurveIndex(1, 0.25)
    @test isapprox(res.t, -0.5)
    @test isapprox(res.ϕ, -0.1)
    res = proj(VecSE2(1.5,0.5,-0.1), curve)
    @test res.ind == CurveIndex(2, 0.25)
    @test isapprox(res.t,  0.5)
    @test isapprox(res.ϕ, -0.1)
    res = proj(VecSE2(10.0,0.0,0.0), curve)
    @test res.ind == CurveIndex(length(curve)-1, 1.0)
    @test isapprox(res.t, 0.0)
    @test isapprox(res.ϕ, 0.0)

    curve2 = [CurvePt(VecSE2(-1.0,-1.0,π/4), 0.0),
              CurvePt(VecSE2( 1.0, 1.0,π/4), hypot(2.0,2.0)),
              CurvePt(VecSE2( 3.0, 3.0,π/4), hypot(4.0,4.0))]

   res = proj(VecSE2(0.0,0.0,0.0), curve2)
   @test res.ind == CurveIndex(1, 0.5)
   @test isapprox(res.t, 0.0)
   @test isapprox(res.ϕ, -π/4)
end


text:
"""Demo 3x7 test case"""
y = [0.62 0.73  0.71 1.50 1.17 0.43 1.08
     0.62 1.73  0.95 1.46 1.60 1.16 0.38
     0.90 0.32 -0.48 0.95 1.08 0.02 0.40]

function write_hpp(fname)
    open(fname, "w") do io
        println(io, "// automatically generated by `julia/test/demo3x7.jl`")
        println(io, "#pragma once")
        println(io)
        let y = permutedims(y)
            print(io, "static const double demo_3x7_y[$(size(y, 1))][$(size(y, 2))] = ")
            println(
                io,
                "{\n",
                join(
                    ("  {" * join((string(y[i, j]) for j in 1:size(y, 2)), ", ") *"}"
                     for i in 1:size(y, 1)),
                    ",\n"
                ),
                "\n};"
            )
        end
    end
end

ext: jl

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
abstract type AbstractDemoCard end

struct UnmatchedCard <: AbstractDemoCard
    path::String
end

"""
    democard(path::String)::T

Constructs a concrete AbstractDemoCard instance.

The return type `T` is determined by the extension of the path
to your demofile. Currently supported types are:

* [`MarkdownDemoCard`](@ref) for markdown files
* [`JuliaDemoCard`](@ref) for julia files
* [`UnmatchedCard`](@ref) for unmatched files

"""
function democard(path::String)::AbstractDemoCard
    validate_file(path)
    _, ext = splitext(path)
    if ext in markdown_exts
        return MarkdownDemoCard(path)
    elseif ext in julia_exts
        return JuliaDemoCard(path)
    else
        return UnmatchedCard(path)
    end
end

basename(x::AbstractDemoCard) = basename(x.path)

function get_default_id(card::AbstractDemoCard)
    # drop leading coutning numbers such as 1. 1_ 1-
    m = match(r"\d*(?<name>.*)", card.title)

    id = replace(strip(m["name"]), ' ' => '-')
end

function validate_id(id::AbstractString, card::AbstractDemoCard)
    if occursin(' ', id)
        throw(ArgumentError("invalid id in $(card.path), it should not contain spaces."))
    end
end

function is_democard(file)
    try
        @suppress_err democard(file)
        return true
    catch err
        @debug err
        return false
    end
end

function load_config(card::T, key; config=Dict()) where T <: AbstractDemoCard
    isempty(config) && (config = parse(card))

    if key == "cover"
        haskey(config, key) || return nothing

        cover_path = config[key]
        if !is_remote_url(cover_path)
            cover_path = replace(cover_path, r"[/\\]" => Base.Filesystem.path_separator) # windows compatibility
        end
        return cover_path
    elseif key == "id"
        haskey(config, key) || return get_default_id(card)

        id = config[key]
        validate_id(id, card)
        return id
    elseif key == "title"
        return get(config, key, get_default_title(card))
    elseif key == "description"
        return get(config, key, card.title)
    elseif key == "hidden"
        return get(config, key, false)
    elseif key == "author"
        return get(config, key, "")
    elseif key == "date"
        return DateTime(get(config, key, DateTime(0)))
    elseif key == "julia"
        version = get(config, key, JULIA_COMPAT)
        return version isa VersionNumber ? version : VersionNumber(string(version))
    else
        throw(ArgumentError("Unrecognized key $(key) for $(T)"))
    end
end

function make_badges(card::AbstractDemoCard)
    badges = []
    if !isempty(card.author)
        for author in split(card.author, ';')
            # TODO: also split on "and"
            isempty(author) && continue

            m = match(regex_md_url, author)
            if isnothing(m)
                author_str = HTTP.escapeuri(strip(author))
                push!(badges, "![Author](https://img.shields.io/badge/Author-$(author_str)-blue)")
            else
                # when markdown url is detected, create a link for it
                # author: [Johnny Chen](https://github.com/johnnychen94)
                name, url = strip.(m.captures)
                name = HTTP.escapeuri(name)
                badge_str = "[![Author](https://img.shields.io/badge/Author-$(name)-blue)]($url)"
                push!(badges, badge_str)
            end
        end
    end
    if card.date != DateTime(0)
        date_str = string(round(Int, datetime2unix(card.date)))
        push!(badges, "![Update time](https://img.shields.io/date/$(date_str))")
    end
    isempty(badges) ? "" : join(badges, " ")
end

### load concrete implementations

include("markdown.jl")
include("julia.jl")


text: include("utility/test.jl")

max_stars_repo_name: kunyuan/DiagramCalculator

max_issues_repo_name: kunyuan/DiagramCalculator

max_forks_repo_name: kunyuan/DiagramCalculator


text:
module Articles

export Article, ArticleMetadata, title, tags, identifier, location,
       datetime, authors, LivePerformance

using Dates
using ..Tags: Tag, TagMatrix, populate!, tagobject
import ..Tags: tags, tag!

"""
Common data across all kinds of articles.
"""
struct ArticleMetadata
    "A unique, machine-friendly identifier for the media."
    identifier :: String

    "A friendly human-readable name describing the media."
    title      :: String

    "A list of tags describing the media."
    tags       :: Vector{Tag}

    "Personnel involved in the production of the media."
    authors    :: Vector{String}

    "A time indicating when the event took place."
    datetime   :: Dates.TimeType

    ArticleMetadata(identifier, title, authors, datetime) =
        new(identifier, title, Tag[], authors, datetime)
end
identifier(m::ArticleMetadata) = m.identifier
title(m::ArticleMetadata) = m.title
tags(m::ArticleMetadata) = m.tags
authors(m::ArticleMetadata) = m.authors
datetime(m::ArticleMetadata) = m.datetime
function tag!(m::ArticleMetadata, tm::TagMatrix, ts)
    tos = [tagobject(tm, t) for t in ts]
    append!(m.tags, tos)
    populate!(tm, tos)
end

"""
A particular item of publication.
"""
abstract type Article end
identifier(x::Article) = identifier(metadata(x))
title(x::Article) = title(metadata(x))
tags(x::Article) = tags(metadata(x))
authors(x::Article) = authors(metadata(x))
datetime(x::Article) = datetime(metadata(x))
tag!(x::Article, tm, t) = tag!(metadata(x), tm, t)

"""
A summary or transcript of a live performance.
"""
struct LivePerformance <: Article
    metadata   :: ArticleMetadata
    "Where the physical event took place."
    location   :: String
end
metadata(x::LivePerformance) = x.metadata
location(x::LivePerformance) = x.location

end


text:
function runNLopt!(  opt,
    p0::Vector{T},
    obj_func,
    grad_func,
    p_lbs,
    p_ubs;
    max_iters = 10000,
    xtol_rel = 1e-12,
    ftol_rel = 1e-12,
    maxtime = Inf) where T
#
    @assert length(p0) == length(p_lbs) == length(p_ubs)

    opt.maxeval = max_iters
    opt.lower_bounds = p_lbs
    opt.upper_bounds = p_ubs
    opt.xtol_rel = xtol_rel
    opt.ftol_rel = ftol_rel
    opt.maxtime = maxtime

    opt.min_objective = (xx, gg)->genericNLoptcostfunc!(xx, gg, obj_func, grad_func)

    # optimize.
    (minf, minx, ret) = NLopt.optimize(opt, p0)

    N_evals = opt.numevals

    return minf, minx, ret, N_evals
end

function genericNLoptcostfunc!(x::Vector{Float64}, df_x, f, df)::Float64

    #
    if length(df_x) > 0
    df_x[:] = df(x)
    end

    return f(x)
end


text:
using Images, ImageFiltering
#=using Plots
const w = 5
function clean_rows(x::Matrix{Float64})
    n = size(x, 1)
    map(1:n - w) do i
        j = i + w
        x[i, j:end]
    end
end
const kernel = Kernel.gaussian((0,1))
function find_peaks(μ::Vector{Float64})
    img = centered(μ')
    s = vec(parent(imfilter(img, kernel)))
    findlocalmaxima(s, 1, false)
end
function get_parameters(row::Vector{Float64})
    m, M = extrema(row)
    m == M && return (0.0, 0.0, 0)
    Δ = M - m
    p = (row[end] - m)/Δ
    return (Δ, p, findfirst(row .> m + Δ/2))
end
function detect_nodes(file::String)
    a = readcsv(file)
    x = Matrix{Float64}(a[2:end, 2:end])
    lengths = Float64.(a[2:end, 1])
    times = Float64.(a[1, 2:end])
    rows = clean_rows(x)
    μ = mean.(rows)
    j = find_peaks(μ)
    map(j) do i
        Δ, p, ti = get_parameters(rows[i])
        ti += i[1] + w
        t = times[ti]
        l = lengths[i]
        (l, t, Δ, p)
    end
end

file = "/home/yakir/.julia/datadeps/all/9/pos1,3-dr5_pos2,4-cle44-luc-1-5/stage 2/root 1/intensities.csv"
p = detect_nodes(file)=#

#=struct Node{T <: Function}
    length_i::Int
    time_i::Int
    label::T
end=#

struct Node
    length_i::Int
    time_i::Int
    Δ::Float64
    p::Float64
end

const kernel = Kernel.gaussian((0,0.8))
function find_peaks(μ::Vector{Float64})
    img = centered(μ')
    s = vec(parent(imfilter(img, kernel)))
    findlocalmaxima(s, 1, false)
end
function get_parameters(row::Vector{Float64})
    m, M = extrema(row)
    m == M && return (0.0, 0.0, 0)
    Δ = M - m
    p = (row[end] - m)/Δ
    return (M, p, findfirst(row .> m + Δ/2))
end
const w = 15
function detect_nodes(track::Track)
    # track = tracks[2][1]
    n = length(track.intensities)
    m = minimum(minimum(i) for i in track.intensities)
    intensities = fill(m, n, n)
    for i in 1:n, j in 1:i-w
        intensities[j,i] = track.intensities[i][j]
    end
    y = vec(mean(intensities, 2))
    j = find_peaks(y)
    nodes = map(j) do i
        Δ, p, ti = get_parameters(intensities[i, :])
        # Node(i[1], ti, x -> text("($Δ, $p)", 5, RGB(1,1,1)))
        # Node(i[1], ti, x -> text("$Δ-$p", RGBA(0,1,0, Int(x > ti)), 8))
        Node(i[1], ti, Δ, p)
    end
    filter!(node -> node.Δ > 0.2, nodes)
    #=rows = [track.intensities[i][1:end-w] for i in w + 1:n]
    y = [mean(rows[i][min(j, i)] for i in 1:n-w) for j in 1:n-w]

    # μ = mean.(rows)
    j = find_peaks(y)

    map(j) do i
        Δ, p, ti = get_parameters(rows[i])
        li = i[1] + w - 1
        ti += li
        Node(li, ti, x -> text("($Δ, $p)", 5, RGB(1,1,1)))
        # Node(li, ti, x -> text("($Δ, $p)", RGBA(1,1,1, Int(x > li)), 8))
    end=#
end


text:
"""
    IOExtras

This module defines extensions to the `Base.IO` interface to support:
 - `startwrite`, `closewrite`, `startread` and `closeread` for streams
    with transactional semantics.
"""
module IOExtras

using ..Sockets
using MbedTLS: MbedException

export bytes, ByteView, nobytes, CodeUnits, IOError, isioerror,
       startwrite, closewrite, startread, closeread,
       tcpsocket, localport, peerport

"""
    bytes(s::String)

Get a `Vector{UInt8}`, a vector of bytes of a string.
"""
function bytes end
bytes(s::SubArray{UInt8}) = unsafe_wrap(Array, pointer(s), length(s))

const CodeUnits = Union{Vector{UInt8}, Base.CodeUnits}
bytes(s::Base.CodeUnits) = bytes(String(s))
bytes(s::String) = codeunits(s)
bytes(s::SubString{String}) = codeunits(s)

bytes(s::Vector{UInt8}) = s

"""
    isioerror(exception)

Is `exception` caused by a possibly recoverable IO error.
"""
isioerror(e) = false
isioerror(::Base.EOFError) = true
isioerror(::Base.IOError) = true
isioerror(e::ArgumentError) = e.msg == "stream is closed or unusable"
isioerror(::MbedException) = true

"""
    IOError <: Exception

The request terminated with due to an IO-related error.

Fields:
 - `e`, the error.
"""
struct IOError <: Exception
    e
    message
end

Base.show(io::IO, e::IOError) = print(io, "IOError(", e.e, " ", e.message, ")\n")

_doc = """
    startwrite(::IO)
    closewrite(::IO)
    startread(::IO)
    closeread(::IO)

Signal start/end of write or read operations.
"""
"$_doc"
startwrite(io) = nothing
"$_doc"
closewrite(io) = nothing
"$_doc"
startread(io) = nothing
"$_doc"
closeread(io) = nothing

using MbedTLS: SSLContext
tcpsocket(io::SSLContext)::TCPSocket = io.bio
tcpsocket(io::TCPSocket)::TCPSocket = io

localport(io) = try !isopen(tcpsocket(io)) ? 0 :
                    Sockets.getsockname(tcpsocket(io))[2]
                catch
                    0
                end

peerport(io) = try !isopen(tcpsocket(io)) ? 0 :
                  Sockets.getpeername(tcpsocket(io))[2]
               catch
                   0
               end

const ByteView = typeof(view(UInt8[], 1:0))
const nobytes = view(UInt8[], 1:0)

"""
Read from an `IO` stream until `find_delimiter(bytes)` returns non-zero.
Return view of bytes up to the delimiter.
"""
function Base.readuntil(buf::IOBuffer,
                    find_delimiter::Function #= Vector{UInt8} -> Int =#
                   )::ByteView

    l = find_delimiter(view(buf.data, buf.ptr:buf.size))
    if l == 0
        return nobytes
    end
    bytes = view(buf.data, buf.ptr:buf.ptr + l - 1)
    buf.ptr += l
    return bytes
end

end


text:
module TestSparse
    using GaussianProcesses
    using GaussianProcesses: get_params, set_params!, update_mll!, update_mll_and_dmll!, init_precompute, predictMVN, FullCovariance, SoR, FITC, DTC, FSA, SubsetOfRegsStrategy, DeterminTrainCondStrat, FullyIndepStrat, FullScaleApproxStrat, getQaa, getQab, predictMVN!
    using Test, Random
    using Distributions: Beta, Normal
    import Calculus
    using LinearAlgebra
    using Statistics
    using PDMats: PDMat

    """ The true function we will be simulating from. """
    function fstar(x::Float64)
        return abs(x-5)*cos(2*x)
    end

    σy = 10.0
    n=1000
    ntest=30

    Random.seed!(1)
    Xdistr = Beta(7,7)
    ϵdistr = Normal(0,σy)
    x = rand(Xdistr, n)*10
    Y = fstar.(x) .+ rand(ϵdistr,n)
    k = SEIso(log(0.3), log(5.0))
    gp_full = GPE(x', Y, MeanConst(mean(Y)), k, log(σy))
    init_params = get_params(gp_full; noise=true, domean=true, kern=true)
    set_params!(gp_full, init_params; noise=true, domean=true, kern=true)
    update_mll!(gp_full)

    Random.seed!(1)
    Xu    = Matrix(quantile(x, [0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.98])')
    xtest = Matrix(rand(Xdistr, ntest)'.*10)

    inearest = [argmin(abs.(xi.-Xu[1,:])) for xi in x]
    blockindices = [findall(isequal(i), inearest) for i in 1:size(Xu,2)]

    function test_pred(gp_sparse, cKPD::PDMat, covstrat::SubsetOfRegsStrategy)
        cK = gp_sparse.cK
        kernel = gp_sparse.kernel
        meanf = gp_sparse.mean
        mx = mean(meanf, xtest)
        xtrain = gp_sparse.x
        alpha = gp_sparse.alpha
        σ2 = exp(2*gp_sparse.logNoise)

        μpred, Σpred = predict_f(gp_sparse, xtest; full_cov=true)

        # see Quiñonero-Candela & Rasmussen 2005, eq. 15
        Qfx = getQab(cK, kernel, xtrain, xtest)
        # Qff = getQaa(cK, kernel, xtrain)
        Qxx = getQaa(cK, kernel, xtest)

        μ_alt, Σ_alt = predictMVN!(Qxx, cKPD, Qfx, mx, alpha)
        @test μ_alt ≈ μpred atol=1e-6
        @test Σ_alt ≈ Σpred rtol=1e-3 # should this be better?
    end
    function test_pred(gp_sparse, cKPD::PDMat, covstrat::Union{DeterminTrainCondStrat,FullyIndepStrat})
        cK = gp_sparse.cK
        kernel = gp_sparse.kernel
        meanf = gp_sparse.mean
        mx = mean(meanf, xtest)
        xtrain = gp_sparse.x
        alpha = gp_sparse.alpha
        σ2 = exp(2*gp_sparse.logNoise)

        μpred, Σpred = predict_f(gp_sparse, xtest; full_cov=true)

        # see Quiñonero-Candela & Rasmussen 2005, eq. 19 and 23
        Qfx = getQab(cK, kernel, xtrain, xtest)
        Kxx = cov(kernel, xtest)

        μ_alt, Σ_alt = predictMVN!(Kxx, cKPD, Qfx, mx, alpha)
        @test μ_alt ≈ μpred atol=1e-6
        @test Σ_alt ≈ Σpred rtol=1e-3 # should this be better?
    end
    function test_pred(gp_sparse, cKPD::PDMat, covstrat::FullScaleApproxStrat)
        cK = gp_sparse.cK
        kernel = gp_sparse.kernel
        meanf = gp_sparse.mean
        mx = mean(meanf, xtest)
        xtrain = gp_sparse.x
        alpha = gp_sparse.alpha
        σ2 = exp(2*gp_sparse.logNoise)

        iprednearest = [argmin(abs.(xi.-Xu[1,:])) for xi in vec(xtest)]
        blockindpred = [findall(isequal(i), iprednearest)
                        for i in 1:size(Xu,2)]
        ;

        μpred, Σpred = predict_f(gp_sparse, xtest, blockindpred; full_cov=true)

        Qfx = getQab(cK, kernel, xtrain, xtest)
        blockindtrain = blockindices
        Λfx = let
            nf, nx = size(xtrain, 2), size(xtest,2)
            zeros(nf, nx)
        end
        for (predblock, trainblock) in zip(blockindpred, blockindtrain)
            Xpredblock, Xtrainblock = xtest[:,predblock], xtrain[:,trainblock]
            Kfx_block = cov(kernel, Xtrainblock, Xpredblock)
            Qfx_block = getQab(cK, kernel, Xtrainblock, Xpredblock)
            Λfx[trainblock,predblock] = Kfx_block - Qfx_block
        end

        Kxx = cov(kernel, xtest)

        μ_alt, Σ_alt = predictMVN!(Kxx, cKPD, Qfx+Λfx, mx, alpha)
        @test μ_alt ≈ μpred atol=1e-6
        @test Σ_alt ≈ Σpred rtol=1e-3 # should this be better?
    end

    function test_sparse(gp_sparse, expect_mll)
        @test gp_sparse.mll ≈ gp_full.mll atol=10 # marginal loglik shouldn't be radically different
        @test gp_sparse.mll ≈ expect_mll atol=1e-6 # from previous run, check this doesn't drift

        cK = gp_sparse.cK
        cKmat = Matrix(cK)
        cKPD = PDMat(cKmat)

        gp_fromsparse = let
            gp = GPE(x', Y, MeanConst(mean(Y)), k, log(σy))
            gp.cK = cKPD
            update_mll!(gp; noise=false, kern=false, domean=true)
            gp
        end
        @test gp_sparse.mll ≈ gp_fromsparse.mll atol=1e-6

        @test tr(cK) ≈ tr(cKmat) rtol=1e-6
        @test logdet(cK) ≈ logdet(cKmat) rtol=1e-6
        xrand = randn(gp_sparse.nobs,2)
        @test cKPD\xrand ≈ cK\xrand atol=1e-6

        buf = init_precompute(gp_sparse)
        update_mll_and_dmll!(gp_sparse, buf; noise=true, domean=true, kern=true)
        grad_analytical = copy(gp_sparse.dmll)
        grad_numerical = Calculus.gradient(init_params) do params
            set_params!(gp_sparse, params; noise=true, domean=true, kern=true)
            GaussianProcesses.update_mll!(gp_sparse)
            t = gp_sparse.mll
            set_params!(gp_sparse, init_params; noise=true, domean=true, kern=true)
            t
        end
        @test grad_numerical ≈ grad_analytical  atol=1e-3
        test_pred(gp_sparse, cKPD, gp_sparse.covstrat)
    end

    @testset "Sparse Approximations" begin
        @testset "Subset of Regressors" begin
            test_sparse(SoR(x', Xu, Y, gp_full.mean, gp_full.kernel, gp_full.logNoise), -3704.0847727395367)
        end
        @testset "Deterministic Training Conditionals" begin
            test_sparse(DTC(x', Xu, Y, gp_full.mean, gp_full.kernel, gp_full.logNoise), -3704.084703493389)
        end
        @testset "Fully Independent Training Conditionals" begin
            test_sparse(FITC(x', Xu, Y, gp_full.mean, gp_full.kernel, gp_full.logNoise), -3709.601737889645)
        end
        @testset "Full Scale Approximation" begin
            test_sparse(FSA(x', Xu, blockindices, Y, gp_full.mean, gp_full.kernel, gp_full.logNoise), -3706.2892293004734)
        end
    end
end


text:
using TensegrityEquilibria, Test, HomotopyContinuation

@testset "circle animation" begin
    resolution = 75
    circleset = [[[sin(2*k*pi/resolution), cos(2*k*pi/resolution)], [sin(2*k*pi/resolution+pi), cos(2*k*pi/resolution+pi)]] for k in 0:resolution]
    animateTensegrity(circleset, [[1,2]], 30)
    sleep(10)
end

@testset "Zeeman interactive" begin
    @var p[1:4];
    display(stableEquilibria([[1.0,0],[2.0,1/4],p[1:2],p[3:4]],[[1,3,0.5]],[[2,3,1/4,1/4],[3,4,1/4,1/4]],p[1:2],p[3:4],[0.0,0.0],[],[]))
    sleep(10)
end

@testset "Zeeman animation" begin
    t=0:1/15:2*pi
    timestamps=[[0.1,0.5*sin(time)-0.25] for time in t]
    @var p[1:4];
    display(stableEquilibria([[1.0,0],[2.0,1/4],p[1:2],p[3:4]],[[1,3,0.5]],[[2,3,1/4,1/4],[3,4,1/4,1/4]],p[1:2],p[3:4],[0.0,0.0],[],[],timestamps))
    sleep(10)
end

@testset "Three points catastrophe" begin
    @var p[1:2] l;
    display(stableEquilibria([[1,2],[3,4],p],[[1,3,l]],[[2,3,1,1]],p,[l],[2.5],[],[]))
    sleep(10)
end

@testset "Bipyramid interactive" begin
    @var p[1:6] ell
    display(stableEquilibria([p[1:3],p[4:6],[0,1,0],[sin(2*pi/3),cos(2*pi/3),0], [sin(4*pi/3),cos(4*pi/3),0]],
                             [[1,2,ell]],
                             [[1,3,1,1],[1,4,1,1],[1,5,1,1],[2,3,1,1],[2,4,1,1],[2,5,1,1]],
                             p,[ell],[1.0],
                             [[3,4],[3,5],[4,5]],[])
    )
    sleep(10)
end

@testset "Bipyramid static" begin
     @var p[1:6]
    display(stableEquilibria([p[1:3],p[4:6],[0,1,0],[sin(2*pi/3),cos(2*pi/3),0], [sin(4*pi/3),cos(4*pi/3),0]],
                             [[1,2,2.0]],
                             [[1,3,1,1],[1,4,1,1],[1,5,1,1],[2,3,1,1],[2,4,1,1],[2,5,1,1]],
                             p[4:6],p[1:3],[0.0,-1.0,0.0],
                             [[3,4],[3,5],[4,5]],[])
    )
    sleep(10)
end

@testset "3-Prism interactive cable constant" begin
    @var p[1:7] c
    display(stableEquilibria([[0,1,0],[sin(2*pi/3),cos(2*pi/3),0],[sin(4*pi/3),cos(4*pi/3),0],[p[1],p[2],p[7]],[p[3],p[4],p[7]],[p[5],p[6],p[7]]],
                             [[1,4,3], [2,5,3],[3,6,3]],
                             [[1,5,2.5,c],[2,6,2.5,c],[3,4,2.5,c],[4,5,1,1.0],[5,6,1,1.0],[4,6,1,1.0]],
                             p[1:6],[c,p[7]],[2.0,2.75],
                             [],
                             [[1,2],[2,3],[1,3]])
    )
    sleep(10)
end

@testset "3-Prism static" begin
    @var p[1:9]
    display(stableEquilibria([[0,1,0],[sin(2*pi/3),cos(2*pi/3),0],[sin(4*pi/3),cos(4*pi/3),0],[p[1],p[2],p[3]],[p[4],p[5],p[6]],[p[7],p[8],p[9]]],
                             [[1,4,3], [2,5,3],[3,6,3]],
                             [[1,5,2.5,1.0],[2,6,2.5,1.0],[3,4,2.5,1.0],[4,5,1,1.0],[5,6,1,1.0],[4,6,1,1.0]],
                             p[1:9],[],[],
                             [],
                             [[1,2],[2,3],[1,3]])
    )
    sleep(10)
end

@testset "3-Prism interactive point" begin
    @var p[1:9]
    display(stableEquilibria([[0,1,0],[sin(2*pi/3),cos(2*pi/3),0],[sin(4*pi/3),cos(4*pi/3),0],[p[1],p[2],p[3]],[p[4],p[5],p[6]],[p[7],p[8],p[9]]],
                             [[1,4,3], [2,5,3],[3,6,3]],
                             [[1,5,2.5,1.0],[2,6,2.5,1.0],[3,4,2.5,1.0],[4,5,1,1.0],[5,6,1,1.0],[4,6,1,1.0]],
                             p[1:6],p[7:9],[0.6053381, 0, 2.5661428],
                             [],
                             [[1,2],[2,3],[1,3]])
    )
    sleep(10)
end

@testset "4-Prism static" begin
    @var p[1:9]
    display(stableEquilibria([[0,0,0], [1,0,0], [1,1,0], [0,1,0], p[1:3], [p[4],p[5],p[3]], [p[6],p[7],p[3]], [p[8],p[9],p[3]]],
                             [[1,7,2], [2,8,2],[3,5,2],[4,6,2]],
                             [[5,6,0.5,1], [6,7,0.5,1],[7,8,0.5,1],[5,8,0.5,1], [1,5,1.,1], [2,6,1.,1], [3,7,1.,1],[4,8,1.,1]],
                             p[1:9],[],[],
                             [],
                             [[1,2], [2,3], [3,4], [1,4]])
    )
    sleep(10)
end


text:
using Test, OMEinsum
using OMEinsum: IndexGroup, NestedEinsum, parse_nested
@testset "einsequence" begin
    @test push!(IndexGroup([],1), 'c').inds == IndexGroup(['c'], 1).inds
    @test isempty(IndexGroup([],1))

    @test_throws ArgumentError parse_nested("((ij,jk),km")
    @test_throws ArgumentError parse_nested("((ij,jk),k1)")

    a, b, c = rand(2,2), rand(2,2), rand(2,2)
    abc1 = ein"(ij,jk),km -> im"(a,b,c)
    abc2 = ein"((ij,jk),km) -> im"(a,b,c)
    abc3 = ein"ij,jk,km -> im"(a,b,c)

    @test abc1 ≈ abc2 ≈ abc3
end

@testset "macro" begin
    b, c, d = rand(2,2), rand(2,2,2), rand(2,2,2,2)
    @ein a[i,j] := b[i,k] * c[k,k,l] * d[l,m,m,j]
    @ein a2[i,j] := b[i,k] * (c[k,k,l] * d[l,m,m,j])
    @test a ≈ a2
end


text:
struct NullData end
const DEFAULT_DATA = Iterators.cycle((NullData(),))

"""
Interactive training with solution animation

"""
function vis_train(
    loss,
    _θ,
    opt,
    _data = DEFAULT_DATA;
    cb = (args...) -> false,
    maxiters = get_maxiters(data),
    progress = true,
    save_best = true,
)

    θ = copy(_θ)
    ps = Flux.params(θ)

    if _data == DEFAULT_DATA && maxiters == typemax(Int)
        error("A data iterator must be provided or the `maxiters` keyword argument must be set.")
    elseif _data == DEFAULT_DATA && maxiters != typemax(Int)
        data = Iterators.repeated((), maxiters)
    elseif maxiters != typemax(Int)
        data = Iterators.take(_data, maxiters) # an iterator that generates at most the first n elements of iter
    else
        data = _data
    end

    t0 = time()

    local x, min_err
    min_err = typemax(eltype(θ)) # dummy variables
    min_opt = 1

    anim = @animate for (i, d) in enumerate(data) # Iterator that yields (i, d) where i is a counter starting at 1,
                                                  # and d is the d-th value from dataset
        # calculate gradients of loss function
        # gradient() requires two args: f and params
        gs = Flux.Zygote.gradient(ps) do
            x = loss(θ, d...)
            first(x)
        end
        cb_call = cb(θ, x...)
        if !(typeof(cb_call) <: Bool)
            error("The callback function should return a boolean to decide whether stopping the optimization process.")
        elseif cb_call # stop
            break
        end

        DiffEqFlux.update!(opt, ps, gs)
        if save_best
            if first(x) < first(min_err) # we've found a better solution
                min_opt = opt
                min_err = x
            end
            if i == maxiters  # for last iteration, revert to best
                opt = min_opt
                cb(θ, min_err...)
            end
        end
    end

    t1 = time()

    res = Optim.MultivariateOptimizationResults(
        opt,
        _θ, # initial_x,
        θ, # pick_best_x(f_incr_pick, state),
        first(x), # pick_best_f(f_incr_pick, state, d),
        maxiters, # iteration,
        maxiters >= maxiters, # iteration == options.iterations,
        false, # x_converged,
        0.0, # T(options.x_tol),
        0.0, # T(options.x_tol),
        NaN, # x_abschange(state),
        NaN, # x_abschange(state),
        false, # f_converged,
        0.0, # T(options.f_tol),
        0.0, # T(options.f_tol),
        NaN, # f_abschange(d, state),
        NaN, # f_abschange(d, state),
        false, # g_converged,
        0.0, # T(options.g_tol),
        NaN, # g_residual(d),
        false, # f_increased,
        nothing,
        maxiters,
        maxiters,
        0,
        true,
        NaN,
        t1 - t0,
    )

    return res, anim

end


text:
using Test
import GDAL
import ArchGDAL as AG

"Test both that an ErrorException is thrown and that the message is as expected"
function eval_ogrerr(err, expected_message)
    @test (@test_throws ErrorException AG.@ogrerr err "e:").value.msg ==
          "e: ($expected_message)"
end

@testset "test_utils.jl" begin
    @testset "metadataitem" begin
        driver = AG.getdriver("DERIVED")
        @test AG.metadataitem(driver, "DMD_EXTENSIONS") == ""
        driver = AG.getdriver("GTiff")
        @test AG.metadataitem(driver, "DMD_EXTENSIONS") == "tif tiff"
    end

    @testset "OGR Errors" begin
        @test isnothing(AG.@ogrerr GDAL.OGRERR_NONE "not an error")
        eval_ogrerr(GDAL.OGRERR_NOT_ENOUGH_DATA, "Not enough data.")
        eval_ogrerr(GDAL.OGRERR_NOT_ENOUGH_MEMORY, "Not enough memory.")
        eval_ogrerr(
            GDAL.OGRERR_UNSUPPORTED_GEOMETRY_TYPE,
            "Unsupported geometry type.",
        )
        eval_ogrerr(GDAL.OGRERR_UNSUPPORTED_OPERATION, "Unsupported operation.")
        eval_ogrerr(GDAL.OGRERR_CORRUPT_DATA, "Corrupt data.")
        eval_ogrerr(GDAL.OGRERR_FAILURE, "Failure.")
        eval_ogrerr(
            GDAL.OGRERR_UNSUPPORTED_SRS,
            "Unsupported spatial reference system.",
        )
        eval_ogrerr(GDAL.OGRERR_INVALID_HANDLE, "Invalid handle.")
        eval_ogrerr(GDAL.OGRERR_NON_EXISTING_FEATURE, "Non-existing feature.")
        # OGRERR_NON_EXISTING_FEATURE is the highest error code currently in GDAL. If another one is
        # added this test will fail.
        eval_ogrerr(GDAL.OGRERR_NON_EXISTING_FEATURE + 1, "Unknown error.")
    end
end


text:
module SimpleSwitch

using SPPL

namespace = @sppl begin
    Y ~ RandInt(0, 4)
    Switch(Y, 0 : 4) do k
        X ~ Bernoulli(1 / (1 + k))
    end
end

namespace.model.logprob(namespace.X << set(1))

end # module


text:
"""
    train_test_year(df::DataFrame, testYears::Union{Int64, Vector{Int64}})

Creates a test dataset containing all chosen years, and a train dataset with the other years.

Column "Date" must be of type Dates.

# Arguments

- `df::DataFrame`: The dataframe containing the data.
- `testYears::Union{Int64, Vector{Int64}}`: The year (or a vector containing the years) to include in the test dataset.
"""

function train_test_year(df::DataFrame, testYears::Union{Int64, Vector{Int64}})
    test = filter(row -> Dates.year(row.Date) in testYears, df)
    years = unique(Dates.year.(df[:, :Date]))
    trainYears = setdiff(years, testYears)
    train = filter(row -> Dates.year(row.Date) in trainYears, df)
    return train, test
end

max_stars_repo_name: AlexLsn/CSO

max_issues_repo_name: AlexLsn/CSO

max_forks_repo_name: AlexLsn/CSO


text:
#####
# This includes short-hands that are convenient but perhaps confusing
#####

export chebyshevt, chebyshevu, legendre, ∫, ⨜, ⨍, ChebyshevWeight, 𝕀, 𝕌, 𝒟

## Constructors

Fun() = Fun(identity)
Fun(d::Domain) = Fun(identity,d)
Fun(d::Space) = Fun(identity,d)

## Chebyshev & Legendre polynomials

chebyshevt(n::Int,d::IntervalOrSegment{T}) where {T<:Number} = Fun(Chebyshev(d),[zeros(T,n);one(T)])
chebyshevu(n::Int,d::IntervalOrSegment{T}) where {T<:Number} =
    mod(n,2) == 1 ? Fun(Chebyshev(d),interlace(zeros(T,div(n+2,2)),2fill(one(T),div(n+2,2)))) :
                    Fun(Chebyshev(d),interlace(2fill(one(T),div(n+2,2)),zeros(T,div(n+2,2)))[1:n+1]-[one(T);zeros(T,n)])
legendre(n::Int,d::IntervalOrSegment{T}) where {T<:Number} = Fun(Legendre(d),[zeros(T,n);one(T)])

for poly in (:chebyshevt,:chebyshevu,:legendre)
    @eval begin
        $poly(::Type{T}, n::Int) where {T<:Number} = $poly(n,ChebyshevInterval{T}())
        $poly(n::Int) = $poly(Float64,n)
        $poly(n::AbstractRange, d::IntervalOrSegment{T}) where {T<:Number} = map(i->$poly(i,d),n)
    end
end

ChebyshevWeight(d,k)=k==0 ? JacobiWeight(-0.5,-0.5,d) : JacobiWeight(0.5,0.5,d)
ChebyshevWeight(d)=ChebyshevWeight(d,0)
ChebyshevWeight(k::Integer)=ChebyshevWeight(Segment(),k)
ChebyshevWeight()=ChebyshevWeight(0)

# shorthand for second order

ivp(d) = ivp(d,2)
bvp(d) = bvp(d,2)

## diff

# use conj(transpose(f)) for ArraySpace
adjoint(f::Fun)=differentiate(f)

∫(f::Fun)=integrate(f)
⨜(f::Fun)=cumsum(f)

for OP in (:Σ,:∮,:⨍,:⨎)
    @eval $OP(f::Fun)=sum(f)
end

∇(F::MultivariateFun) = grad(F)
dot(∇::Function,F::Vector{M}) where {M<:MultivariateFun} = div(F)
cross(∇::Function,F::Vector{M}) where {M<:MultivariateFun} = curl(F)

## Domains

const 𝕀 = ChebyshevInterval()
const ℝ = Line()
const 𝕌 = Circle()

𝒟 = Derivative()
Δ = Laplacian()


text:
module NN

using Flux
using Flux: onehotbatch, onecold
using MLDatasets
using Random
using CUDA
using BSON: @save, @load

export trainmodel
export accuracy
export getprocessedMNIST
export getprocessedFashionMNIST

function getprocessedMNIST()
	train_imgs, train_labels = MLDatasets.MNIST.traindata(Float32)
	test_imgs, test_labels = MLDatasets.MNIST.testdata(Float32)

    # Reshape from (width, height, samples) to (width, height, channels (1), samples)
    train_imgs = Flux.unsqueeze(train_imgs, ndims(train_imgs))
    test_imgs = Flux.unsqueeze(test_imgs, ndims(test_imgs))

	train_labels = Float32.(onehotbatch(train_labels, 0:9))
    test_labels = Float32.(onehotbatch(test_labels, 0:9))

    return (
        train=(
            inputs=train_imgs,
            labels=train_labels
        ),
        test=(
            inputs=test_imgs,
            labels=test_labels
        )
    )
end

function getprocessedFashionMNIST()
	train_imgs, train_labels = MLDatasets.FashionMNIST.traindata(Float32)
	test_imgs, test_labels = MLDatasets.FashionMNIST.testdata(Float32)

    # Reshape from (width, height, samples) to (width, height, channels (1), samples)
    train_imgs = Flux.unsqueeze(train_imgs, ndims(train_imgs))
    test_imgs = Flux.unsqueeze(test_imgs, ndims(test_imgs))

	train_labels = Float32.(onehotbatch(train_labels, 0:9))
    test_labels = Float32.(onehotbatch(test_labels, 0:9))

    return (
        train=(
            inputs=train_imgs,
            labels=train_labels
        ),
        test=(
            inputs=test_imgs,
            labels=test_labels
        )
    )
end

percstr(x::Real; digits=2) = "$(round(100x; digits=digits))%"
accuracy(model, x, y) = sum(onecold(cpu(model(x))) .== onecold(cpu(y))) / size(x)[end]

function trainmodel(;
    name::String,
    model,
    dataset,
    loss,
    opt,
    batchsize::Int,
    epochs::Int,
    device = gpu,
)
    savefilename = "saves/$name.bson"
    # Check for preexisting trained weights
    if (isfile(savefilename))
        @info "Found a preexisting trained model in file '$savefilename'. Loading from it..."
        @load savefilename parameters test_acc
        model = cpu(model)
        try
            Flux.loadparams!(model, parameters)
            @info "Model loaded. Accuracy: $(percstr(test_acc))"
            return device(model)
        catch e
            @warn "Current model is not compatible with the weights in file '$savefilename'"
            @info model
            @info "Training from scratch..."
            @warn "File '$savefilename' will be overwritten after training"
        end
    end
    @info "Performing experiment `$name`..."
    @info "Loading model and dataset to $((device == gpu) && has_cuda() ? "GPU" : "CPU")..."
    dataset = device(dataset)
    model = device(model)

    train_dataloader = Flux.Data.DataLoader((dataset.train.inputs, dataset.train.labels), batchsize=batchsize, shuffle=false)

    loss′(x, y) = loss(model(x), y)

    @info "Compiling model and calculating initial accuracy..."
    # Keep track of best results so far so we can roll back to it at the end of training
    best_acc = accuracy(model, dataset.test.inputs, dataset.test.labels)    # The test set is used here for debug purposes. This should be overwritten soon
    best_params = deepcopy(params(model))
    @info "Initial test accuracy: $(percstr(best_acc))"

    avgsecsperepoch = 0f0

    @info "Training..."
    # Those flushes increase responsiveness, specially if `stderr` is redirected to a file
    flush(stderr)
    for epoch_idx in 1:epochs
        # Train for a single epoch
        secsthisepoch = @elapsed Flux.train!(loss′, params(model), train_dataloader, opt)

        # The first epoch time is discarded as it includes (large) compilation times
        if epoch_idx > 1
            avgsecsperepoch += (secsthisepoch - avgsecsperepoch) / (epoch_idx - 1)
        end

        train_acc = accuracy(model, dataset.train.inputs, dataset.train.labels)

        # Log progress highlighting in red new best result
        @info "$name [$epoch_idx/$epochs] Train Accuracy: $(percstr(train_acc))\t($(round(avgsecsperepoch; digits=2)) secs/epoch)$(train_acc > best_acc ? "\t\e[91mNew best acc\e[m" : "")"

        # Register the best accuracy so far
        if train_acc > best_acc
            best_acc = train_acc
            best_params = deepcopy(params(model))
        end

        flush(stderr)
    end

    @info "Time spent training: $(round(avgsecsperepoch*epochs; digits=2)) secs"

    # Rewind weights to the best performing found during training
    Flux.loadparams!(model, best_params)
    test_acc = accuracy(model, dataset.test.inputs, dataset.test.labels)
    @save savefilename parameters=params(cpu(model)) test_acc=test_acc
    @info "Test accuracy: $(percstr(test_acc))"

    return model
end

end

max_stars_repo_name: ArthurWalraven/cnnslth

max_issues_repo_name: ArthurWalraven/cnnslth

max_forks_repo_name: ArthurWalraven/cnnslth


text:
using Distances
using Clustering
using Statistics
using LinearAlgebra
using Distances
using RGBDSegmentation: CDNImage
"""
```
function initial_centroids(plain_data,  k::Integer)
```
Return a ``9 \times k`` matrix containing the centroids of each  cluster.

Each centroid is represented as a 9-dimensional vector in which the
first three components represent the color components,
the 4:6 indices represent the spatial components and
the last three component represent the normal component.

Each slice of the array are normalized. i.e.
``norm(c[1:3, j]) = 1``

"""
function initial_centroids(plain_data,   k::Integer)
    centroids = zeros(9, k)
    indices = [1:3, 4:6, 7:9]
    for i=1:3
        model = kmeans(
            plain_data[indices[i], :],
            k,
            maxiter=500)
        labels = assignments(model)
        for j=1:k
            centroids[indices[i], j] = mean(
                plain_data[indices[i], findall(labels.==j)],
                dims=2)
	    centroids[indices[i], j] /= norm(centroids[indices[i], j])
        end
    end
    return centroids
end
function indmin(x)
    (vals, cidxs) = findmin(x, dims=2)
    idx = [cidx[2] for cidx in cidxs]
    return (vals,idx)
end
"""
```
function initial_clusters(p_data::Array{T, 2},
                          k::Integer) where T <: Number
```
Return the initial assigments using a combined distance among the patterns of the three images.
"""
function initial_clusters(p_data::Array{T, 2},
                          k::Integer) where T <: Number

    centroids = initial_centroids(p_data, k)
    indices = [1:3, 4:6, 7:9]
    diff = 10000;
    epsilon = 0.1;
    value = 100;
    iteration = 1;
    maxnumIter = 50;
    clust = nothing;
    while (diff > epsilon) && (maxnumIter>iteration)
        iteration = iteration+1;
        oldvalue = value;
        distance_to_center = zeros(size(p_data,2), k)
        distances = [Euclidean(), Euclidean(), CosineDist()]
        for i=1:3, j=1:k
            distance_to_center[:,j] += colwise(
                distances[i],
                p_data[indices[i],:],
                vec(centroids[indices[i],j])
            )
        end
        # assign points to nearest cluster (based on minimum distance)
        (distMin,clust) = indmin(distance_to_center)
        value = sum(distMin);
        # compute new cluster centroids
	for h=1:k
    	  try
	    cluster_indices = vec(clust.==h)
            centroids[1:3, h] = mean(
                p_data[1:3, cluster_indices],
                dims=2)
            centroids[4:6, h] = mean(p_data[4:6,cluster_indices],
                                    dims=2)
            centroids[7:9, h] = sum(p_data[7:9,cluster_indices],
                                    dims=2)
            centroids[7:9, h] = centroids[7:9, h]./ norm(centroids[7:9, h]);
	   catch e

                succ = 0;
         	return;
	   end
	end
	indices_invalid_centroid = unique([id[2] for id in vec(findall(isnan.(centroids)))])
    	if (!isempty(indices_invalid_centroid))
            centroids = centroids[:, setdiff(1:k, indices_invalid_centroid)]
	    k = k - length(indices_invalid_centroid);
	end
	diff = abs(value - oldvalue);
    end
    return clust
end


text:
type LinSolveFactorize{F}
  factorization::F
  A
end
LinSolveFactorize(factorization) = LinSolveFactorize(factorization,nothing)
function (p::LinSolveFactorize)(x,A,b,update_matrix=false)
  if update_matrix
    p.A = p.factorization(A)
  end
  A_ldiv_B!(x,p.A,b)
end
DEFAULT_LINSOLVE = LinSolveFactorize(lufact!)


text:
struct Cumulative{T}
    cum_curr::Vector{SharedMatrix{T}}
    max_curr::Vector{SharedMatrix{T}}
    cum_branch_curr::Vector{SharedVector{T}}
    cum_node_curr::Vector{SharedVector{T}}
end

struct GraphData{T,V}
    G::SparseMatrixCSC{T,V}
    cc::Vector{Vector{V}}
    points::Vector{V}
    user_points::Vector{V}
    exclude_pairs::Vector{Tuple{V,V}}
    nodemap::Matrix{V}
    polymap::Matrix{V}
    hbmeta::RasterMeta
    cellmap::Matrix{T}
    cum::Cumulative{T}
end

struct ComponentData{T,V}
    cc::Vector{V}
    matrix::SparseMatrixCSC{T,V}
    local_nodemap::Matrix{V}
    hbmeta::RasterMeta
    cellmap::Matrix{T}
end

struct Output{T,V}
    points::Vector{V}
    voltages::Vector{T}
    orig_pts::Tuple{V,V}
    comp_idx::Tuple{V,V}
    resistance::T
    col::V
    cum::Cumulative{T}
end

struct Shortcut{T}
    get_shortcut_resistances::Bool
    voltmatrix::Matrix{T}
    shortcut_res::Matrix{T}
end

"""
Core kernel of Circuitscape - used to solve several pairs

Input:
* data::GraphData
"""
function single_ground_all_pairs(data::GraphData{T,V}, flags, cfg, log = true) where {T,V}

    if flags.solver in AMG
        csinfo("Solver used: AMG accelerated by CG", cfg["suppress_messages"] in TRUELIST)
        amg_solver_path(data, flags, cfg, log)
    else
        csinfo("Solver used: CHOLMOD", cfg["suppress_messages"] in TRUELIST)
        bs = parse(Int, cfg["cholmod_batch_size"])
        _cholmod_solver_path(data, flags, cfg, log, bs)
    end
end

function amg_solver_path(data::GraphData{T,V}, flags, cfg, log)::Matrix{T} where {T,V}

    # Data
    a = data.G
    cc = data.cc
    points = data.points
    exclude = data.exclude_pairs
    nodemap = data.nodemap
    polymap = data.polymap
    orig_pts = data.user_points
    hbmeta = data.hbmeta
    cellmap = data.cellmap

    # Flags
    outputflags = flags.outputflags
    is_raster = flags.is_raster
    write_volt_maps = outputflags.write_volt_maps
    write_cur_maps = outputflags.write_cur_maps
    write_cum_cur_map_only = outputflags.write_cum_cur_map_only
    write_max_cur_maps = outputflags.write_max_cur_maps

    # Get number of focal points
    numpoints = size(points, 1)

    # Cumulative currents
    cum = data.cum

    csinfo("Graph has $(size(a,1)) nodes, $numpoints focal points and $(length(cc)) connected components", cfg["suppress_messages"] in TRUELIST)

    num, d = get_num_pairs(cc, points, exclude)
    log && csinfo("Total number of pair solves = $num", cfg["suppress_messages"] in TRUELIST)

    # Initialize pairwise resistance
    resistances = -1 * ones(T, numpoints, numpoints)::Matrix{T}
    voltmatrix = zeros(T, size(resistances))::Matrix{T}
    shortcut_res = deepcopy(resistances)::Matrix{T}

    # Get a vector of connected components
    comps = getindex.([a], cc, cc)

    get_shortcut_resistances = false
    if is_raster && !write_volt_maps && !write_cur_maps &&
            !write_cum_cur_map_only && !write_max_cur_maps &&
            isempty(exclude)
        get_shortcut_resistances = true
        csinfo("Triggering resistance calculation shortcut", cfg["suppress_messages"] in TRUELIST)
        num, d = get_num_pairs_shortcut(cc, points, exclude)
        csinfo("Total number of pair solves has been reduced to $num ", cfg["suppress_messages"] in TRUELIST)
    end
    shortcut = Shortcut(get_shortcut_resistances, voltmatrix, shortcut_res)

    for (cid, comp) in enumerate(cc)

        # Subset of points relevant to CC
        csub = filter(x -> x in comp, points) |> unique
        #idx = findin(c, csub)

        if isempty(csub)
            continue
        end

        # Conductance matrix corresponding to CC
        matrix = comps[cid]

        # Regularization step
        matrix.nzval .+= eps(eltype(matrix)) * norm(matrix.nzval)

        # Construct preconditioner *once* for every CC
        t1 = @elapsed P = aspreconditioner(smoothed_aggregation(matrix))
        csinfo("Time taken to construct preconditioner = $t1 seconds", cfg["suppress_messages"] in TRUELIST)

        # Get local nodemap for CC - useful for output writing
        t2 = @elapsed local_nodemap = construct_local_node_map(nodemap, comp, polymap)
        csinfo("Time taken to construct local nodemap = $t2 seconds", cfg["suppress_messages"] in TRUELIST)

        component_data = ComponentData(comp, matrix, local_nodemap, hbmeta, cellmap)

        function f(i)

            # Generate return type
            ret = Vector{Tuple{V,V,T}}()

            pi = csub[i]
            comp_i = something(findfirst(isequal(pi),comp), 0)
            comp_i = V(comp_i)
            I = findall(x -> x == pi, points)
            smash_repeats!(ret, I)

            # Preprocess matrix
            # d = matrix[comp_i, comp_i]

            # Iteration space through all possible pairs
            rng = i+1:size(csub, 1)
            if nprocs() > 1
                for j in rng
                    pj = csub[j]
                    csinfo("Scheduling pair $(d[(pi,pj)]) of $num to be solved", cfg["suppress_messages"] in TRUELIST)
                end
            end

            # Loop through all possible pairs
            for j in rng

                pj = csub[j]
                comp_j = something(findfirst(isequal(pj), comp), 0)
                comp_j = V(comp_j)
                J = findall(x -> x == pj, points)

                if pi == pj
                    continue
                end

                # Forget excluded pairs
                ex = false
                for c_i in I
                    for c_j in J
                        if (c_i, c_j) in exclude
                            continue
                        end

                        # Initialize currents
                        current = zeros(T, size(matrix, 1))
                        current[comp_i] = -1
                        current[comp_j] = 1

                        # Solve system
                        # csinfo("Solving points $pi and $pj")
                        log && csinfo("Solving pair $(d[(pi,pj)]) of $num", cfg["suppress_messages"] in TRUELIST)
                        t2 = @elapsed v = solve_linear_system(cfg, matrix, current, P)
                        csinfo("Time taken to solve linear system = $t2 seconds", cfg["suppress_messages"] in TRUELIST)
                        v .= v .- v[comp_i]

                        # Calculate resistance
                        r = v[comp_j] - v[comp_i]

                        # Return resistance value
                        push!(ret, (c_i, c_j, r))
                        if get_shortcut_resistances
                            resistances[c_i, c_j] = r
                            resistances[c_j, c_i] = r
                        end
                        output = Output(points, v, (orig_pts[c_i], orig_pts[c_j]),
                                        (comp_i, comp_j), r, V(c_j), cum)
                        postprocess(output, component_data, flags, shortcut, cfg)
                    end
                end
            end

        # matrix[comp_i, comp_i] = d

        ret
        end

        if get_shortcut_resistances
            idx = something(findfirst(isequal(csub[1]), points), 0)
            f(1)
            update_shortcut_resistances!(idx, shortcut, resistances, points, comp)
        else
            is_parallel = cfg["parallelize"] in TRUELIST
            if is_parallel
                X = pmap(x ->f(x), 1:size(csub,1))
            else
                X = map(x ->f(x), 1:size(csub,1))
            end

            # Set all resistances
            for x in X
                for i = 1:size(x, 1)
                    resistances[x[i][1], x[i][2]] = x[i][3]
                    resistances[x[i][2], x[i][1]] = x[i][3]
                end
            end
        end

    end

    if get_shortcut_resistances
        resistances = shortcut.shortcut_res
    end

    for i = 1:size(resistances,1)
        resistances[i,i] = 0
    end

    # Pad it with the user points
    r = vcat(vcat(0,orig_pts)', hcat(orig_pts, resistances))

    # Save resistances
    save_resistances(r, cfg)

    r
end

struct CholmodNode{T}
    cc_idx::Tuple{T,T}
    points_idx::Tuple{T,T}
end

function _cholmod_solver_path(data::GraphData{T,V}, flags,
                                  cfg, log, batch_size = 1000) where {T,V}

    # Data
    a = data.G
    cc = data.cc
    points = data.points
    exclude = data.exclude_pairs
    nodemap = data.nodemap
    polymap = data.polymap
    orig_pts = data.user_points
    hbmeta = data.hbmeta
    cellmap = data.cellmap

    # Flags
    outputflags = flags.outputflags
    is_raster = flags.is_raster
    write_volt_maps = outputflags.write_volt_maps
    write_cur_maps = outputflags.write_cur_maps
    write_cum_cur_map_only = outputflags.write_cum_cur_map_only
    write_max_cur_maps = outputflags.write_max_cur_maps

    # Cumulative current map
    cum = data.cum

    # CHOLMOD solver mode works only in double precision
    if eltype(a) == Float32
        cswarn("Converting single precision matrix to double")
        a = Float64.(a)
    end

    # Get number of focal points
    numpoints = size(points, 1)

    csinfo("Graph has $(size(a,1)) nodes, $numpoints focal points and $(length(cc)) connected components", cfg["suppress_messages"] in TRUELIST)

    num, d = get_num_pairs(cc, points, exclude)
    log && csinfo("Total number of pair solves = $num", cfg["suppress_messages"] in TRUELIST)

    # Initialize pairwise resistance
    resistances = -1 * ones(eltype(a), numpoints, numpoints)
    voltmatrix = zeros(eltype(a), size(resistances))
    shortcut_res = -1 * ones(eltype(a), size(resistances))

    # Get a vector of connected components
    comps = getindex.([a], cc, cc)

    get_shortcut_resistances = false
    if is_raster && !write_volt_maps && !write_cur_maps &&
            !write_cum_cur_map_only  && !write_max_cur_maps &&
            isempty(exclude)
        get_shortcut_resistances = true
        csinfo("Triggering resistance calculation shortcut", cfg["suppress_messages"] in TRUELIST)
        num, d = get_num_pairs_shortcut(cc, points, exclude)
        csinfo("Total number of pair solves has been reduced to $num ", cfg["suppress_messages"] in TRUELIST)
    end
    shortcut = Shortcut(get_shortcut_resistances, voltmatrix, shortcut_res)

    for (cid, comp) in enumerate(cc)

        # Subset of points relevant to CC
        csub = filter(x -> x in comp, points) |> unique
        #idx = findin(c, csub)

        if isempty(csub)
            continue
        end

        # Conductance matrix corresponding to CC
        matrix = comps[cid]

        t = @elapsed factor = construct_cholesky_factor(matrix)
        csinfo("Time taken to construct cholesky factor = $t", cfg["suppress_messages"] in TRUELIST)

        # Get local nodemap for CC - useful for output writing
        t2 = @elapsed local_nodemap = construct_local_node_map(nodemap, comp, polymap)
        csinfo("Time taken to construct local nodemap = $t2 seconds", cfg["suppress_messages"] in TRUELIST)

        component_data = ComponentData(comp, matrix, local_nodemap, hbmeta, cellmap)

        ret = Vector{Tuple{V,V,Float64}}()

        cholmod_batch = CholmodNode[]

        # Batched backsubstitution
        function g(i)

            pi = csub[i]
            comp_i = V(something(findfirst(isequal(pi),comp),0))
            I = findall(x -> x == pi, points)
            # smash_repeats!(ret, I)
            smash_repeats!(resistances, I)

            # Iteration space through all possible pairs
            rng = i+1:size(csub, 1)

            # Loop through all possible pairs
            for j in rng

                pj = csub[j]
                comp_j = V(something(findfirst(isequal(pj), comp),0))
                J = findall(x -> x == pj, points)

                if pi == pj
                    continue
                end

                # Forget excluded pairs
                for c_i in I, c_j in J
                    if (c_i, c_j) in exclude
                        continue
                    else
                        push!(cholmod_batch,
                          CholmodNode((comp_i, comp_j), (V(c_i), V(c_j))))
                    end
                end
            end
        end

        function f(i, rng, lhs)
            v = rng[i]
            output = Output(points, lhs[:,i],
                (orig_pts[cholmod_batch[v].points_idx[1]],
                orig_pts[cholmod_batch[v].points_idx[2]]),
                cholmod_batch[v].cc_idx,
                lhs[cholmod_batch[v].cc_idx[2], i] - lhs[cholmod_batch[v].cc_idx[1], i],
                V(cholmod_batch[v].points_idx[2]), cum)
            postprocess(output, component_data, flags, shortcut, cfg)
        end
        if get_shortcut_resistances
            idx = something(findfirst(isequal(csub[1]), points),0)
            g(1)
        else
            g.(1:size(csub, 1))
        end

        # l = min(num, length(cholmod_batch))
        l = length(cholmod_batch)

        for st in 1:batch_size:l

            rng = st + batch_size <= l ?
                            (st:(st+batch_size-1)) : (st:l)

            csinfo("Solving points $(rng.start) to $(rng.stop)", cfg["suppress_messages"] in TRUELIST)

            rhs = zeros(eltype(matrix), size(matrix, 1), length(rng))

            for (i,v) in enumerate(rng)
                node = cholmod_batch[v]
                rhs[node.cc_idx[1], i] = -1
                rhs[node.cc_idx[2], i] = 1
            end
            lhs = factor \ rhs

            # Normalisation step
            for (i,val) in enumerate(rng)
                n = cholmod_batch[val].cc_idx[1]
                v = lhs[n,i]
                for j = 1:size(matrix, 1)
                    lhs[j,i] = lhs[j,i] - v
                end
            end

            is_parallel = cfg["parallelize"] in TRUELIST
            if is_parallel
                X = pmap(x -> f(x, rng, lhs), 1:length(rng))
            else
                X = map(x -> f(x, rng, lhs), 1:length(rng))
            end

            for (i,v) in enumerate(rng)
                coords = cholmod_batch[v].points_idx
                r = lhs[cholmod_batch[v].cc_idx[2], i] -
                            lhs[cholmod_batch[v].cc_idx[1], i]
                resistances[coords...] = r
                resistances[reverse(coords)...] = r
            end
        end

        if get_shortcut_resistances
            update_shortcut_resistances!(idx, shortcut, resistances, points, comp)
        end
    end

    if get_shortcut_resistances
        resistances = shortcut.shortcut_res
    end

    for i = 1:size(resistances,1)
        resistances[i,i] = 0
    end

    # Pad it with the user points
    r = vcat(vcat(0,orig_pts)', hcat(orig_pts, resistances))

    # Save resistances
    save_resistances(r, cfg)

    r
end

function construct_cholesky_factor(matrix)
    cholesky(matrix + sparse(10eps()*I,size(matrix)...))
end

"""
Returns all possible pairs to solve.

Input:
* ccs::Vector{Vector{Int}} - vector of connected components
* fp::Vector{Int} - vector of focal points
* exclude_pairs::Vector{Tuple{Int,Int}} - vector of point pairs (tuples) to exclude

Output:
* n - total number of pairs
"""
function get_num_pairs(ccs, fp::Vector{V}, exclude_pairs) where V

    num = 0
    d = Dict{Tuple{V,V}, V}()

    for (i,cc) in enumerate(ccs)
        sub_fp = filter(x -> x in cc, fp) |> unique
        l = lastindex(sub_fp)
        for ii = 1:l
            pt1 = sub_fp[ii]
            for jj = ii+1:l
                pt2 = sub_fp[jj]
                if (pt1, pt2) in exclude_pairs
                    continue
                else
                    num += 1
                    d[(pt1, pt2)] = num
                end
            end
        end
    end
    num, d
end

function get_num_pairs_shortcut(ccs, fp::Vector{V}, exclude_pairs) where V

    num = 0
    d = Dict{Tuple{V,V}, V}()

    for (i,cc) in enumerate(ccs)
        sub_fp = filter(x -> x in cc, fp) |> unique
        l = lastindex(sub_fp)
        l == 0 && continue
        for ii = 1:1
            pt1 = sub_fp[ii]
            for jj = ii+1:l
                pt2 = sub_fp[jj]
                if (pt1, pt2) in exclude_pairs
                    continue
                else
                    num += 1
                    d[(pt1, pt2)] = num
                end
            end
        end
    end
    num, d
end
function smash_repeats!(ret, I)
    for i = 1:size(I,1)
        for j = i+1:size(I,1)
            push!(ret, (I[i], I[j], 0))
        end
    end
end

function smash_repeats!(resistances::Matrix{T}, I) where T
    for i = 1:size(I,1)
        for j = i+1:size(I,1)
            resistances[I[i], I[j]] = 0
            resistances[I[j], I[i]] = 0
        end
    end
end

"""
Calculate laplacian of the adjacency matrix of a graph
"""
function laplacian(G::SparseMatrixCSC{T,V}) where {T,V}
    n = size(G, 1)
    s = Vector{eltype(G)}(undef,n)
    for i = 1:n
        s[i] = sum_off_diag(G, i)
        for j in nzrange(G, i)
            if i == G.rowval[j]
                G.nzval[j] = 0
            else
                G.nzval[j] = -G.nzval[j]
            end
        end
    end
    r = V(1):V(n)
    S = sparse(r, r, s)
    G + S
end

function sum_off_diag(G, i)
     sum = zero(eltype(G))
     for j in nzrange(G, i)
         if G.rowval[j] != i
             sum += G.nzval[j]
         end
     end
     sum
 end

function solve_linear_system(cfg,
            G::SparseMatrixCSC{T,V},
            curr::Vector{T}, M)::Vector{T} where {T,V}
    cg(G, curr, Pl = M, tol = T(1e-6), maxiter = 100_000)
end

function postprocess(output, component_data, flags, shortcut, cfg)

    voltages = output.voltages
    matrix = component_data.matrix
    local_nodemap = component_data.local_nodemap
    hbmeta = component_data.hbmeta
    orig_pts = output.orig_pts

    # Shortcut flags and data
    get_shortcut_resistances = shortcut.get_shortcut_resistances

    if get_shortcut_resistances
        update_voltmatrix!(shortcut, output, component_data)
        return nothing
    end

    name = "_$(orig_pts[1])_$(orig_pts[2])"

    if flags.outputflags.write_volt_maps
        t = @elapsed write_volt_maps(name, output, component_data, flags, cfg)
        csinfo("Time taken to write voltage maps = $t seconds", cfg["suppress_messages"] in TRUELIST)
    end

    if flags.outputflags.write_cur_maps
        t = @elapsed write_cur_maps(name, output, component_data,
                                    [-9999.], flags, cfg)
        csinfo("Time taken to write current maps = $t seconds", cfg["suppress_messages"] in TRUELIST)
    end
    nothing
end

function update_voltmatrix!(shortcut, output, component_data)

    # Data
    voltmatrix = shortcut.voltmatrix
    c = output.points
    cc = component_data.cc
    voltages = output.voltages
    r = output.resistance
    j = output.col

    for i = 2:size(c, 1)
        ind = something(findfirst(isequal(c[i]), cc),0)
        if ind != 0
            voltageAtPoint = voltages[ind]
            voltageAtPoint = 1 - (voltageAtPoint/r)
            voltmatrix[i,j] = voltageAtPoint
        end
    end
end

function update_shortcut_resistances!(anchor, sc, resistances, points, comp)

    # Data
    voltmatrix = sc.voltmatrix
    shortcut = sc.shortcut_res

    check = map(x -> x in comp, points)
    l = size(resistances, 1)
    for pointx = 1:l
        if check[pointx]
            R1x = resistances[anchor, pointx]
            if R1x != -1
                shortcut[pointx, anchor] = shortcut[anchor, pointx] = R1x
                for point2 = pointx:l
                    if check[point2]
                        R12 = resistances[anchor, point2]
                        if R12 != -1
                            if R1x != -777
                                shortcut[anchor, point2] = shortcut[point2, anchor] = R12
                                Vx = voltmatrix[pointx, point2]
                                R2x = 2*R12*Vx + R1x - R12
                                if shortcut[point2, pointx] != -777
                                    shortcut[point2, pointx] = shortcut[pointx, point2] = R2x
                                end
                            else
                                shortcut[pointx, :] = shortcut[:, pointx] = -777
                            end
                        end
                    end
                end
            end
        end
    end
end


text:
__precompile__()

module Fire

import Base: @__doc__, doc
import Base.Meta: quot

export @main

struct FuncDef
    name::Symbol
    args::Vector
    kwargs::Vector
    func::Function
end

const BASIC_NUMBER_TYPE   = Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, Int128, UInt128, Float16, Float32, Float64, BigInt, BigFloat
const BASIC_STRING_TYPE   = String, Symbol
const BASIC_ABSTRACT_TYPE = Signed, Unsigned, Integer, AbstractFloat, Real, Number, AbstractString
const BASIC_TYPES         = (BASIC_NUMBER_TYPE..., BASIC_STRING_TYPE..., BASIC_ABSTRACT_TYPE...)

const REQUIRED  = 0x01
const OPTIONAL  = 0x02
const VARLENGTH = 0x03

const entries = FuncDef[]

function __init__()
    atexit(parse_command_line)
end

function parse_function_definition(f)
    if f.head != Symbol("function")
        throw(ArgumentError("@main can only be used on function definitions"))
    end
    name = f.args[1].args[1]
    defs = f.args[1].args[2:end]
    args, kwargs = Expr(:vect), Expr(:vect)

    for def in defs
        if isa(def, Symbol)
            push!(args.args, Expr(:tuple, quot(def), String, REQUIRED))
        elseif def.head == :(::)
            push!(args.args, Expr(:tuple, quot(def.args[1]), esc(def.args[2]), REQUIRED))
        elseif def.head == :kw
            if isa(def.args[1], Symbol)
                push!(args.args, Expr(:tuple, quot(def.args[1]), String, OPTIONAL, quot(def.args[2])))
            else
                push!(args.args, Expr(:tuple, quot(def.args[1].args[1]), esc(def.args[1].args[2]), OPTIONAL, quot(def.args[2])))
            end
        elseif def.head == :(...)
            def = def.args[1]
            if isa(def, Symbol)
                push!(args.args, Expr(:tuple, quot(def), String, VARLENGTH))
            else
                push!(args.args, Expr(:tuple, quot(def.args[1]), esc(def.args[2]), VARLENGTH))
            end
        elseif def.head == :parameters
            for def in def.args
                if def.head != :kw
                    error("Fire.jl don't support var length keyword argument")
                end
                if isa(def.args[1], Symbol)
                    push!(kwargs.args, Expr(:tuple, quot(def.args[1]), String, quot(def.args[2])))
                else
                    push!(kwargs.args, Expr(:tuple, quot(def.args[1].args[1]), esc(def.args[1].args[2]), quot(def.args[2])))
                end
            end
        else
            error("your function is too complex for Fire.jl")
        end
    end

    quot(name), args, kwargs, esc(name)
end

function parse_with_type(t, x)
    t in BASIC_NUMBER_TYPE             ? parse(t, x) :
    t in BASIC_STRING_TYPE             ? t(x) :
    t in (Signed, Integer)             ? parse(Int, x) :
    t == Unsigned                      ? parse(UInt, x) :
    t in (AbstractFloat, Real, Number) ? parse(Float64, x) :
    t == AbstractString                ? x :
    error("BUG")
end

function parse_command_line()
    isempty(entries) && return

    let ARGS = copy(ARGS)
        pargs, oargs = [], Dict{Symbol, Any}()
        entry = if length(entries) > 1
            isempty(ARGS) && return print_help_all()
            cmd = popfirst!(ARGS)
            if cmd == "--help"
                return print_help_all()
            else
                i = findfirst(x->x.name == Symbol(replace(cmd, '-' => '_')), entries)
                i == nothing && return println(stderr, "ERROR: unknown command $cmd, see --help for more avaliable commands")
                entries[i]
            end
        else
            entries[]
        end

        for kwarg in entry.kwargs
            if kwarg[2] == Bool
                oargs[kwarg[1]] = false
            end
        end

        i = 1

        while !isempty(ARGS)
            x = popfirst!(ARGS)
            if startswith(x, "--")
                if x == "--help"
                    return print_help_entry(entry)
                end

                x = Symbol(x[3:end])

                arg = let ind = findfirst(y->y[1] == Symbol(replace(String(x), '-' => '_')), entry.kwargs)
                    ind != nothing ? entry.kwargs[ind] : return println(stderr, "Unknown option --$x, see --help for avaliable options")
                end

                if arg[2] == Bool
                    oargs[arg[1]] = true
                elseif arg[2] <: Vector
                    oargs[arg[1]] = arg[2]()
                    while !isempty(ARGS)
                        x = ARGS[1]
                        startswith(x, "--") ? break : popfirst!(ARGS)
                        push!(oargs[arg[1]], parse_with_type(eltype(arg[2]), x))
                    end
                elseif arg[2] in BASIC_TYPES
                    if isempty(ARGS) || (x = popfirst!(ARGS); startswith(x, "--"))
                        return println(stderr, "Option --$(arg[1]) need an argument")
                    end

                    oargs[arg[1]] = parse_with_type(arg[2], x)
                else
                    return println(stderr, "Argument type $(arg[2]) is not supported by Fire.jl")
                end
            else
                if i > length(entry.args)
                    return println(stderr, "Too many arguments, see --help for list of arguments")
                end

                arg = entry.args[i]

                arg[2] in BASIC_TYPES || return println(stderr, "Argument type $(arg[2]) is not supported by Fire.jl")

                push!(pargs, parse_with_type(arg[2], x))

                if arg[3] != VARLENGTH
                    i += 1
                end
            end
        end

        number_of_required_arguments = count(x->x[3]==REQUIRED, entry.args)

        if length(pargs) < number_of_required_arguments
            return println(stderr, "Need $number_of_required_arguments positional arguments, see --help for what are them")
        end

        entry.func(pargs...; oargs...)
    end
end

function print_help_all()
    println("See --help of each command for usages")
    for entry in entries
        println("  ", replace(String(entry.name), '_' => '-'))
    end
    println()
end

function print_help_entry(entry)
    println(doc(entry.func))
    if !isempty(entry.args)
        println("Positional Arguments:")
        for i in entry.args
            print("  ", replace(String(i[1]), '_' => '-'), ": ", i[2])
            i[3] == REQUIRED  ? println() :
            i[3] == OPTIONAL  ? println(" (default: $(i[4]))") :
            i[3] == VARLENGTH ? println("...") :
            error("BUG")
        end
        println()
    end
    if !isempty(entry.kwargs)
        println("Optional Arguments:")
        for i in entry.kwargs
            println("  --", replace(String(i[1]), '_' => '-'), ": ", i[2], " (default: $(i[3]))")
        end
        println()
    end
end

macro main(f::Expr)
    x = parse_function_definition(f)
    quote
        @__doc__ $(esc(f))
        push!(entries, FuncDef($(x...)))
    end
end

end # module Fire


text:
### Reference state
using DocStringExtensions
using ..TemperatureProfiles
export ReferenceState, NoReferenceState, HydrostaticState
const TD = Thermodynamics
using CLIMAParameters.Planet: R_d, MSLP, cp_d, grav, T_surf_ref, T_min_ref
using ..DGMethods: fvm_balance!
using ..Mesh.Grids: polynomialorders

"""
    ReferenceState

Hydrostatic reference state, for example, used as initial
condition or for linearization.
"""
abstract type ReferenceState end

vars_state(m::ReferenceState, ::AbstractStateType, FT) = @vars()

"""
    NoReferenceState <: ReferenceState

No reference state used
"""
struct NoReferenceState <: ReferenceState end

"""
    HydrostaticState{P,T} <: ReferenceState

A hydrostatic state specified by a virtual
temperature profile and relative humidity.

By default, this is a dry hydrostatic reference
state.
"""
struct HydrostaticState{P, FT} <: ReferenceState
    virtual_temperature_profile::P
    relative_humidity::FT
    subtract_off::Bool
end
"""
    HydrostaticState(
        virtual_temperature_profile,
        relative_humidity = 0;
        subtract_off = true,
    )

Construct a `HydrostaticState` given virtual temperature profile and
relative humidity. The keyword argument `subtract_off` controls
whether the constructed state is subtracted off in the prognostic
momentum equation to remove hydrostatic contribution.
"""
function HydrostaticState(
    virtual_temperature_profile::TemperatureProfile{FT},
    relative_humidity = FT(0);
    subtract_off = true,
) where {FT}
    return HydrostaticState{typeof(virtual_temperature_profile), FT}(
        virtual_temperature_profile,
        relative_humidity,
        subtract_off,
    )
end

# TODO: change `ρe` to `energy` to support `θModel`
# TODO: Make `moisture` sub-component, so that moisture components are optional
vars_state(m::HydrostaticState, ::Auxiliary, FT) =
    @vars(ρ::FT, p::FT, T::FT, ρe::FT, ρq_tot::FT, ρq_liq::FT, ρq_ice::FT)

function ref_state_init_pρ!(
    atmos::AtmosModel,
    aux::Vars,
    tmp::Vars,
    geom::LocalGeometry,
)
    param_set = parameter_set(atmos)
    z = altitude(atmos, aux)
    ref_state = reference_state(atmos)
    T_virt, p = ref_state.virtual_temperature_profile(param_set, z)
    aux.ref_state.p = p
    aux.ref_state.ρ = p / (T_virt * R_d(param_set))
end

function ref_state_init_density_from_pressure!(
    atmos::AtmosModel,
    aux::Vars,
    tmp::Vars,
    geom::LocalGeometry,
) where {P, F}
    k = vertical_unit_vector(atmos, aux)
    ∇Φ = ∇gravitational_potential(atmos, aux)
    # density computation from pressure ρ = -1/g*dpdz
    ρ = -k' * tmp.∇p / (k' * ∇Φ)
    aux.ref_state.ρ = ρ
end

function ref_state_finalize_init!(
    atmos::AtmosModel,
    aux::Vars,
    tmp::Vars,
    geom::LocalGeometry,
)
    FT = eltype(aux)
    param_set = parameter_set(atmos)
    ρ = aux.ref_state.ρ
    p = aux.ref_state.p
    T_virt = p / (ρ * FT(R_d(param_set)))

    ref_state = reference_state(atmos)
    RH = ref_state.relative_humidity
    phase_type = PhaseEquil
    (T, q_pt) = TD.temperature_and_humidity_given_TᵥρRH(
        param_set,
        T_virt,
        ρ,
        RH,
        phase_type,
    )

    # Update temperature to be exactly consistent with
    # p, ρ, and q_pt
    if moisture_model(atmos) isa DryModel
        ts = PhaseDry_ρp(param_set, ρ, p)
    else
        ts = PhaseEquil_ρpq(param_set, ρ, p, q_pt.tot)
    end
    T = air_temperature(ts)
    q_pt = PhasePartition(ts)
    q_tot = q_pt.tot
    q_liq = q_pt.liq
    q_ice = q_pt.ice

    aux.ref_state.ρq_tot = ρ * q_tot
    aux.ref_state.ρq_liq = ρ * q_liq
    aux.ref_state.ρq_ice = ρ * q_ice
    aux.ref_state.T = T
    e_kin = FT(0)
    e_pot = gravitational_potential(atmos.orientation, aux)
    aux.ref_state.ρe = ρ * total_energy(e_kin, e_pot, ts)
end

atmos_init_aux!(::AtmosModel, ::ReferenceState, _...) = nothing
function atmos_init_aux!(
    atmos::AtmosModel,
    ::HydrostaticState,
    state_auxiliary::MPIStateArray,
    grid,
    direction,
)
    # Step 1: initialize both ρ and p
    init_state_auxiliary!(
        atmos,
        ref_state_init_pρ!,
        state_auxiliary,
        grid,
        direction,
    )

    # Step 2: correct ρ from p to satisfy discrete hydrostic balance
    vertical_fvm = polynomialorders(grid)[end] == 0

    if vertical_fvm
        # Vertical finite volume scheme
        # pᵢ - pᵢ₊₁ =  g ρᵢ₊₁ Δzᵢ₊₁/2 + g ρᵢ Δzᵢ/2
        fvm_balance!(fvm_balance_init!, atmos, state_auxiliary, grid)
    else
        ref_state = reference_state(atmos)
        ∇p = ∇reference_pressure(ref_state, state_auxiliary, grid)
        init_state_auxiliary!(
            atmos,
            ref_state_init_density_from_pressure!,
            state_auxiliary,
            grid,
            direction;
            state_temporary = ∇p,
        )
    end

    init_state_auxiliary!(
        atmos,
        ref_state_finalize_init!,
        state_auxiliary,
        grid,
        direction,
    )
end

using ..MPIStateArrays: vars
using ..DGMethods: init_ode_state
using ..DGMethods.NumericalFluxes:
    CentralNumericalFluxFirstOrder,
    CentralNumericalFluxSecondOrder,
    CentralNumericalFluxGradient

"""
    PressureGradientModel

A mini balance law that is used to take the gradient of reference
pressure. The gradient is computed as ∇ ⋅(pI) and the calculation
uses the balance law interface to be numerically consistent with
the way this gradient is computed in the dynamics.
"""
struct PressureGradientModel <: BalanceLaw end
vars_state(::PressureGradientModel, ::Auxiliary, T) = @vars(p::T)
vars_state(::PressureGradientModel, ::Prognostic, T) = @vars(∇p::SVector{3, T})
vars_state(::PressureGradientModel, ::Gradient, T) = @vars()
vars_state(::PressureGradientModel, ::GradientFlux, T) = @vars()
function init_state_auxiliary!(
    m::PressureGradientModel,
    state_auxiliary::MPIStateArray,
    grid,
    direction,
) end
function init_state_prognostic!(
    ::PressureGradientModel,
    state::Vars,
    aux::Vars,
    localgeo,
    t,
) end
function flux_first_order!(
    ::PressureGradientModel,
    flux::Grad,
    state::Vars,
    auxstate::Vars,
    t::Real,
    direction,
)
    flux.∇p -= auxstate.p * I
end
flux_second_order!(::PressureGradientModel, _...) = nothing
source!(::PressureGradientModel, _...) = nothing

boundary_conditions(::PressureGradientModel) = ntuple(i -> nothing, 6)
boundary_state!(nf, ::Nothing, ::PressureGradientModel, _...) = nothing

∇reference_pressure(::NoReferenceState, state_auxiliary, grid) = nothing
function ∇reference_pressure(::ReferenceState, state_auxiliary, grid)
    FT = eltype(state_auxiliary)
    ∇p = similar(state_auxiliary; vars = @vars(∇p::SVector{3, FT}), nstate = 3)

    grad_model = PressureGradientModel()
    # Note that the choice of numerical fluxes doesn't matter
    # for taking the gradient of a continuous field
    grad_dg = DGModel(
        grad_model,
        grid,
        CentralNumericalFluxFirstOrder(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient(),
    )

    # initialize p
    ix_p = varsindex(vars(state_auxiliary), :ref_state, :p)
    grad_dg.state_auxiliary.data .= state_auxiliary.data[:, ix_p, :]

    # FIXME: this isn't used but needs to be passed in
    gradQ = init_ode_state(grad_dg, FT(0))

    grad_dg(∇p, gradQ, nothing, FT(0))
    return ∇p
end

function fvm_balance_init!(
    atmos::AtmosModel,
    aux_top::Vars,
    aux::Vars,
    Δz::MArray{Tuple{2}, FT},
) where {FT}
    param_set = parameter_set(atmos)
    _grav::FT = grav(param_set)
    _R_d::FT = R_d(param_set)

    ref = aux.ref_state
    topref = aux_top.ref_state

    T_virt = ref.p / (_R_d * ref.ρ)
    top_T_virt = topref.p / (_R_d * topref.ρ)

    topref.ρ =
        ref.ρ * (_R_d * T_virt - _grav * Δz[1] / 2) /
        (_R_d * top_T_virt + _grav * Δz[2] / 2)
    topref.p = _R_d * topref.ρ * top_T_virt
end


text:
using DashBootstrapComponents

carousel = dbc_carousel(
    items = [
        Dict("key" => "1", "src" => "/static/images/slide1.svg"),
        Dict("key" => "2", "src" => "/static/images/slide2.svg"),
        Dict("key" => "3", "src" => "/static/images/slide3.svg"),
    ],
    controls = false,
    indicators = false,
    interval = 2000,
    ride = "carousel",
);

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
function print_rec(H, W)
    for i in 1:H
        for j in 1:W
            print("#")
        end
        println()
    end
end

while true
    H, W = parse.(Int, split(readline()))
    (H == 0 && W == 0) && break
    print_rec(H, W)
    println()
end


text:
export Role

"""
A [`User`](@ref) role.
More details [here](https://discordapp.com/developers/docs/topics/permissions#role-object).
"""
struct Role <: DiscordObject
    id::Snowflake
    name::String
    color::Optional{Int}  # These fields are missing in audit log entries.
    hoist::Optional{Bool}
    position::Optional{Int}
    permissions::Optional{String}
    managed::Optional{Bool}
    mentionable::Optional{Bool}
end
@boilerplate Role :constructors :docs :lower :merge :mock


text:
using BenchmarkTools, Plots, Plots.PlotMeasures, LaTeXStrings
using BenchmarkTools: minimum, median

SUITE = BenchmarkGroup()
SUITE["LaubLoomis"] = BenchmarkGroup()

# ==============================================================================
# Jet-based approach using Taylor Models
# ==============================================================================
include("laubloomis.jl")

# ----------------------------------------
# --- Case 1: smaller initial states ---
# ----------------------------------------
𝑃, 𝑂 = laubloomis(W=0.01, property=(t,x)->x[4] < 4.5)

𝑂₁ = Options(:abs_tol=>1e-10, :orderT=>7, :orderQ=>1, :max_steps=>1000)

# first run
sol_case_1 = solve(𝑃, 𝑂, op=TMJets(𝑂₁))

# verify that specification holds
v4 = [0.0, 1.0] # the flowpipe has been projected so we check for the second component which is x4
@assert all([ρ(v4, sol_case_1.Xk[i].X) < 4.5 for i in eachindex(sol_case_1.Xk)])

# print width of final box
H = sol_case_1.Xk[end].X
println("width of final box, case W = 0.01 : $(high(H)[2] - low(H)[2])")

# benchmark
SUITE["LaubLoomis"]["W=0.01"] = @benchmarkable solve($𝑃, $𝑂, op=TMJets($𝑂₁))

# ----------------------------------------
# --- Case 2: intermediate initial states ---
# ----------------------------------------
𝑃, 𝑂 = laubloomis(W=0.05, property=(t,x)->x[4] < 4.5)

𝑂₂ = Options(:abs_tol=>1e-10, :orderT=>7, :orderQ=>1, :max_steps=>1000)

# first run
sol_case_2 = solve(𝑃, 𝑂, op=TMJets(𝑂₂))

# verify that specification holds
@assert all([ρ(v4, sol_case_2.Xk[i].X) < 4.5 for i in eachindex(sol_case_2.Xk)])

# print width of final box
H = sol_case_2.Xk[end].X
println("width of final box, case W = 0.05 : $(high(H)[2] - low(H)[2])")

# benchmark
SUITE["LaubLoomis"]["W=0.05"] = @benchmarkable solve($𝑃, $𝑂, op=TMJets($𝑂₂))

# ----------------------------------------
# --- Case 3: larger initial states ---
# ----------------------------------------
𝑃, 𝑂 = laubloomis(W=0.1, property=(t,x)->x[4] < 5.0)

𝑂₃ = Options(:abs_tol=>1e-10, :orderT=>7, :orderQ=>1, :max_steps=>1000)

# first run
sol_case_3 = solve(𝑃, 𝑂, op=TMJets(𝑂₃))

# verify that specification holds
@assert all([ρ(v4, sol_case_3.Xk[i].X) < 5.0 for i in eachindex(sol_case_3.Xk)])

# print width of final box
H = sol_case_3.Xk[end].X
println("width of final box, case W = 0.1 : $(high(H)[2] - low(H)[2])")

# benchmark
SUITE["LaubLoomis"]["W=0.1"] = @benchmarkable solve($𝑃, $𝑂, op=TMJets($𝑂₃))

# ==============================================================================
# Execute benchmarks and save benchmark results
# ==============================================================================

# tune parameters
tune!(SUITE)

# run the benchmarks
results = run(SUITE, verbose=true)

# return the sample with the smallest time value in each test
println("minimum time for each benchmark:\n", minimum(results))

# return the median for each test
println("median time for each benchmark:\n", median(results))

# ==============================================================================
# Execute benchmarks and save benchmark results
# ==============================================================================

plot(sol_case_1,
     tickfont=font(30, "Times"), guidefontsize=45,
     xlab=L"t\raisebox{2.0mm}{\textcolor{white}{.}}",
     ylab=L"x_{4}\raisebox{1.2mm}{\textcolor{white}{.}}",
     xtick=[0., 2., 4., 6., 8., 10., 12., 14., 16., 18., 20.],
     ytick=[2, 2.5, 3, 3.5, 4, 4.5],
     xlims=(0., 20.), ylims=(1.5, 4.5),
     bottom_margin=6mm, left_margin=8mm, right_margin=4mm, top_margin=3mm,
     size=(1000, 1000), linecolor="blue")

plot!(x->x, x->4.5, 0., 20., line=2, color="red", linestyle=:dash, legend=nothing)
savefig(@relpath "laubloomis_case_1.png")

plot(sol_case_2,
     tickfont=font(30, "Times"), guidefontsize=45,
     xlab=L"t\raisebox{2.0mm}{\textcolor{white}{.}}",
     ylab=L"x_{4}\raisebox{1.2mm}{\textcolor{white}{.}}",
     xtick=[0., 2., 4., 6., 8., 10., 12., 14., 16., 18., 20.],
     ytick=[2, 2.5, 3, 3.5, 4, 4.5, 5.0],
     xlims=(0., 20.), ylims=(1.5, 5.0),
     bottom_margin=6mm, left_margin=8mm, right_margin=4mm, top_margin=3mm,
     size=(1000, 1000), linecolor="blue")

plot!(x->x, x->5.0, 0., 20., line=2, color="red", linestyle=:dash, legend=nothing)
savefig(@relpath "laubloomis_case_2.png")

plot(sol_case_3,
     tickfont=font(30, "Times"), guidefontsize=45,
     xlab=L"t\raisebox{2.0mm}{\textcolor{white}{.}}",
     ylab=L"x_{4}\raisebox{1.2mm}{\textcolor{white}{.}}",
     xtick=[0., 2., 4., 6., 8., 10., 12., 14., 16., 18., 20.],
     ytick=[2, 2.5, 3, 3.5, 4, 4.5, 5.0],
     xlims=(0., 20.), ylims=(1.5, 5.0),
     bottom_margin=6mm, left_margin=8mm, right_margin=4mm, top_margin=3mm,
     size=(1000, 1000), linecolor="blue")

plot!(x->x, x->5.0, 0., 20., line=2, color="red", linestyle=:dash, legend=nothing)
savefig(@relpath "laubloomis_case_3.png")

plot(sol_case_1, color="red")

plot!(sol_case_2, alpha=0.6, color="green")

plot!(sol_case_3, alpha=0.2,
     tickfont=font(30, "Times"), guidefontsize=45,
     xlab=L"t\raisebox{2.0mm}{\textcolor{white}{.}}",
     ylab=L"x_{4}\raisebox{1.2mm}{\textcolor{white}{.}}",
     xtick=[0., 2., 4., 6., 8., 10., 12., 14., 16., 18., 20.],
     ytick=[2, 2.5, 3, 3.5, 4, 4.5, 5.0],
     xlims=(0., 20.), ylims=(1.5, 5.0),
     bottom_margin=6mm, left_margin=8mm, right_margin=4mm, top_margin=3mm,
     size=(1000, 1000), color="blue")

plot!(x->x, x->5.0, 0., 20., line=2, color="red", linestyle=:dash, legend=nothing)
plot!(x->x, x->4.5, 0., 20., line=2, color="red", linestyle=:dash, legend=nothing)
savefig(@relpath "laubloomis_case_all.png")

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
using CSV, DataFrames, Dates, PyCall

export JuliaToPandas
export PandasToJulia
export LSTM
export train
export save_lstm
export load_lstm
export predict
export ragged_preds
export gen_news

"convert a Julia DataFrame to a Pandas DataFrame"
function JuliaToPandas(df::DataFrame)
    pd = pyimport("pandas")

    date_col_name = (eltype.(eachcol(df)) .|> string .== "Date") |> x-> names(df)[x] |> x-> x[1]
    CSV.write("tmp.csv", df)
    pd_df = pd.read_csv("tmp.csv", parse_dates=[date_col_name])
    run(`rm tmp.csv`)
    return pd_df
end

"convert a Pandas DataFrame to a Julia DataFrame"
function PandasToJulia(df::PyObject, missingstrings::Vector=["NA", "na", "nan", "NAN", "NaN", ".", ""])
    df.to_csv("tmp.csv", index=false)
    julia_df = CSV.File("tmp.csv", missingstrings=missingstrings) |> DataFrame
    run(`rm tmp.csv`)
    return julia_df
end

"""
    Primary class of the library, used for transforming data, training the model, and making predictions.
    parameters:
        data : DataFrame
            n x m+1 dataframe
        target_variable : String
            name of the target var
        n_timesteps : Int
            how many historical periods to consider when training the model. For example if the original data is monthly, n_steps=12 would consider data for the last year.
        fill_na_func : PyObject
            Python function to replace within-series NAs. Given a column, the function should return a scalar. E.g. to replace default mean with median, `LSTM(..., fill_na_func = pyimport("numpy").nanmedian)`
        fill_ragged_edges_func : PyObject
            Python function to replace NAs in ragged edges (data missing at end of series). Pass "ARMA" for ARMA filling. Not ARMA filling will be significantly slower as models have to be estimated for each variable to fill ragged edges.
        n_models : Int
            number of models to train and take the average of for more robust estimates
        train_episodes : Int
            number of epochs/episodes to train the model
        batch_size : Int
            number of observations per training batch
        decay : Float32
            learning rate decay
        n_hidden : Int
            number of hidden states in the network
        n_layers : Int
            number of LSTM layers in the network
        dropout : Float32
            dropout rate between the LSTM layers
        criterion : PyObject
            torch loss criterion, defaults to MAE. E.g. `LSTM(..., pyimport("torch").nn.L1Loss())`. Need parentheses at end of torch loss function.
        optimizer : PyObject
            torch optimizer, defaults to Adam. E.g. `LSTM(..., pyimport("torch").optim.Adam)`. Don't need parentheses at end of torch loss function.
        optimizer_parameters : Dict
            list of parameters for optimizer, including learning rate. E.g. Dict("lr" => 1e-2)
    returns: PyObject
        instantiated LSTM model
"""
function LSTM(;
    data::DataFrame = nothing,
    target_variable::String = nothing,
    n_timesteps::Int = nothing,
    fill_na_func::PyObject = pyimport("numpy").nanmean,
    fill_ragged_edges_func = pyimport("numpy").nanmean,
    n_models::Int = 1,
    train_episodes::Int = 200,
    batch_size::Int = 30,
    decay::Float64 = 0.98,
    n_hidden::Int = 20,
    n_layers::Int = 2,
    dropout::Float64 = 0.0,
    criterion::PyObject = pyimport("torch").nn.L1Loss(),
    optimizer::PyObject = pyimport("torch").optim.Adam,
    optimizer_parameters::Dict = Dict("lr" => 1e-2)
)::PyObject

    model = pyimport("nowcast_lstm.LSTM").LSTM(
        data = JuliaToPandas(data),
        target_variable = target_variable,
        n_timesteps = n_timesteps,
        fill_na_func = fill_na_func,
        fill_ragged_edges_func = fill_ragged_edges_func,
        n_models = n_models,
        train_episodes = train_episodes,
        batch_size = batch_size,
        decay = decay,
        n_hidden = n_hidden,
        n_layers = n_layers,
        dropout = dropout,
        criterion = criterion,
        optimizer = optimizer,
        optimizer_parameters = optimizer_parameters,
    )

    return model
end

"train an instantiated LSTM model. pass train(model; quiet=true). quiet=true will suppress printing of training loss at each epoch"
function train(model::PyObject; quiet::Bool = false)::PyObject
    model.train(quiet=quiet)
    return model
end

"get predictions from a trained model. pass predict(model, data; only\\_actuals\\_obs = true) to only get predictions where there is an actual value"
function predict(model::PyObject, data::DataFrame; only_actuals_obs::Bool = false)::DataFrame
    model.predict(JuliaToPandas(data), only_actuals_obs) |> PandasToJulia
end

"save a trained LSTM model to disk. pass save_lstm(model, path_to_save), filename must end in .pkl, e.g. '/Users/trained_model.pkl'"
function save_lstm(model::PyObject, filename::String)
    pyimport("dill").dump(model, py"open"(filename, mode="wb"))
end

"load a trained LSTM model. pass load_lstm(path_to_pkl)"
function load_lstm(filename::String)::PyObject
    model = pyimport("dill").load(py"open"(filename, "rb", -1))
    return model
end

"""
    Get predictions on artificial vintages
    parameters:
        model : PyObject
            trained LSTM model
        pub_lags : Vector
            list of periods back each input variable is set to missing. I.e. publication lag of the variable.
        lag : Int
            simulated periods back. E.g. -2 = simulating data as it would have been 2 months before target period, 1 = 1 month after, etc.
        data : DataFrame
            dataframe to generate the ragged datasets on, if none will calculate on training data
        start_date : String
            String in "YYYY-MM-DD" format: start date of generating ragged preds. To save calculation time, i.e. just calculating after testing date instead of all dates
        end_date : String
            String in "YYYY-MM-DD" format: end date of generating ragged preds
    returns: DataFrame
        predictions on vintages
"""
function ragged_preds(
    model::PyObject,
    pub_lags::Vector{Int64},
    lag::Int64,
    data::DataFrame;
    start_date::String = "",
    end_date::String = "",
)::DataFrame
    if start_date == ""
        start_date = nothing
    end
    if end_date == ""
        end_date = nothing
    end
    model.ragged_preds(pub_lags, lag, JuliaToPandas(data), start_date, end_date) |> PandasToJulia
end

"""
    Generate the news between two data releases using the method of holding out new data feature by feature and recording the differences in model output
    parameters:
        model : PyObject
            trained LSTM model
        target_period : String
            String in "YYYY-MM-DD", target prediction date
        old_data : DataFrame
             dataframe of previous dataset
        new_data : DataFrame
             dataframe of new dataset
    returns: Dictionary
        "news": dataframe of news contribution of each column with updated data. scaled_news is news scaled to sum to actual prediction delta.
        "old_pred": prediction on the previous dataset
        "new_pred": prediction on the new dataset
        "holdout_discrepency": difference between the sum of news via the holdout method and the actual prediction delta
"""
function gen_news(model::PyObject, target_period::String, old_data::DataFrame, new_data::DataFrame)::Dict
    news = model.gen_news(target_period, JuliaToPandas(old_data), JuliaToPandas(new_data))
    julia_news = Dict()
    julia_news["news"] = news["news"] |> PandasToJulia
    julia_news["old_pred"] = news["old_pred"]
    julia_news["new_pred"] = news["new_pred"]
    julia_news["holdout_discrepency"] = news["holdout_discrepency"]
    return julia_news
end