text:
using BinaryBuilder, Pkg

name = "Goalign"
version = v"0.3.5"
hash = "9593e5fba495e67b4c15aeee53fd292547518c2c"

# Collection of sources required to build pprof
sources = [
    GitSource("https://github.com/evolbioinfo/goalign.git", hash),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/goalign/

VERSION_PACKAGE=github.com/evolbioinfo/goalign/version.Version
NAME=${bindir}/goalign${exeext}

mkdir -p ${bindir}
go build -o ${NAME} -ldflags "-X ${VERSION_PACKAGE}=${version}"
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    ExecutableProduct("goalign", :goalign),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; compilers=[:c, :go], julia_compat="1.6")

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
using BinaryProvider

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, ["libOsi"], :libOsi),
    LibraryProduct(prefix, ["libOsiCommonTests"], :libOsiCommonTests),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/JuliaBinaryWrappers/Osi_jll.jl/releases/download/Osi-v0.108.5+2"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc_any, :cxx03)) => ("$bin_prefix/Osi.v0.108.5.aarch64-linux-gnu-cxx03.tar.gz", "b22132f3806c28cda7f860de9606a2d3cf821846704cf0605fa3d271faf91231"),
    Linux(:aarch64, libc=:glibc, compiler_abi=CompilerABI(:gcc_any, :cxx11)) => ("$bin_prefix/Osi.v0.108.5.aarch64-linux-gnu-cxx11.tar.gz", "85c97895c3087ab9e70ce91d1bf8ebc974b781909b2cd889482e3ba2c5226129"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc_any, :cxx03)) => ("$bin_prefix/Osi.v0.108.5.aarch64-linux-musl-cxx03.tar.gz", "7564b897047ed37f95d8dafa1109a46851188cd822673c48739df0d1ddfadd4b"),
    Linux(:aarch64, libc=:musl, compiler_abi=CompilerABI(:gcc_any, :cxx11)) => ("$bin_prefix/Osi.v0.108.5.aarch64-linux-musl-cxx11.tar.gz", "c532560d2e1892ba5ebb8c8392b2642d519dbd4c4f669b31a5b1b4f89bc43777"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc_any, :cxx03)) => ("$bin_prefix/Osi.v0.108.5.armv7l-linux-gnueabihf-cxx03.tar.gz", "a179769ac9dff73c8a2b53cc0d56b31034cad9618b3e9fd1ec186061df2b1498"),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc_any, :cxx11)) => ("$bin_prefix/Osi.v0.108.5.armv7l-linux-gnueabihf-cxx11.tar.gz", "193e5f5c75877843db3b9c8890c9b53f9c5049bfe1fd55a345f3ecaeff200e53"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc_any, :cxx03)) => ("$bin_prefix/Osi.v0.108.5.armv7l-linux-musleabihf-cxx03.tar.gz", "ed0113954e87e6335ace16ee5d4a80768538b6a1f19594db6d4bd1b69bd3b7f2"),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf, compiler_abi=CompilerABI(:gcc_any, :cxx11)) => ("$bin_prefix/Osi.v0.108.5.armv7l-linux-musleabihf-cxx11.tar.gz", "1071c0aae46951f98997b7afe735fddb01f3cc95fa096d9f25e4dc6476273e56"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc_any, :cxx03)) => ("$bin_prefix/Osi.v0.108.5.i686-linux-gnu-cxx03.tar.gz", "84dd3e580b4a68457c0a2f2e1d52f24a3ac96254b2dd69029f1fe3adcbe78ae0"),
    Linux(:i686, libc=:glibc, compiler_abi=CompilerABI(:gcc_any, :cxx11)) => ("$bin_prefix/Osi.v0.108.5.i686-linux-gnu-cxx11.tar.gz", "7123079b46213a355190a89ba2bbc54cea49b9ea4e9a2dc7a749cace07f32ae4"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc_any, :cxx03)) => ("$bin_prefix/Osi.v0.108.5.i686-linux-musl-cxx03.tar.gz", "78231a09e6d77d8ba38519aed44cf1d2c8135bf95b5c546522e1981424a4b58d"),
    Linux(:i686, libc=:musl, compiler_abi=CompilerABI(:gcc_any, :cxx11)) => ("$bin_prefix/Osi.v0.108.5.i686-linux-musl-cxx11.tar.gz", "21ddd9a2e31609b126439e2dbdf2e85b29ed233f095322958fec8c562f2a1cca"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc_any, :cxx03)) => ("$bin_prefix/Osi.v0.108.5.i686-w64-mingw32-cxx03.tar.gz", "485d349d4631ae90f0d50ce16a2771846f5f81514101b0633abceeec26fc73f0"),
    Windows(:i686, compiler_abi=CompilerABI(:gcc_any, :cxx11)) => ("$bin_prefix/Osi.v0.108.5.i686-w64-mingw32-cxx11.tar.gz", "df32d376da4fb9e69e93cf2b8bf38132437cf2c39e65830539b85e1967c253b3"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc_any, :cxx03)) => ("$bin_prefix/Osi.v0.108.5.x86_64-apple-darwin14-cxx03.tar.gz", "e0aea3fd07013c8740a48e84bd2d7e9a5bec7755979be27911b085753f68171c"),
    MacOS(:x86_64, compiler_abi=CompilerABI(:gcc_any, :cxx11)) => ("$bin_prefix/Osi.v0.108.5.x86_64-apple-darwin14-cxx11.tar.gz", "53da6cf3d2c73f0a2bc7ee4d9855894214a1507be744571d859a610e1ee099a6"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc_any, :cxx03)) => ("$bin_prefix/Osi.v0.108.5.x86_64-linux-gnu-cxx03.tar.gz", "495a529761f2dfa8d494901dda3a1cf575d12a381cf3824a11bcf07a58646d56"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc_any, :cxx11)) => ("$bin_prefix/Osi.v0.108.5.x86_64-linux-gnu-cxx11.tar.gz", "49003732bc3ad79103ba7c3fdd463b2356d0aedd587a10589943c9e1de591edb"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc_any, :cxx03)) => ("$bin_prefix/Osi.v0.108.5.x86_64-linux-musl-cxx03.tar.gz", "bc9a1ad02fca2b93a7f4229831de7f0917a09ee00cfb934ae849e46975a6f505"),
    Linux(:x86_64, libc=:musl, compiler_abi=CompilerABI(:gcc_any, :cxx11)) => ("$bin_prefix/Osi.v0.108.5.x86_64-linux-musl-cxx11.tar.gz", "96855a8175373605401915a8b163fbe97300de1337b006e3284d31b0559480bf"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc_any, :cxx03)) => ("$bin_prefix/Osi.v0.108.5.x86_64-w64-mingw32-cxx03.tar.gz", "f3c436822340c0076c61f220669f67c3d32c6420128b6baa059ba70ea6ae4c70"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc_any, :cxx11)) => ("$bin_prefix/Osi.v0.108.5.x86_64-w64-mingw32-cxx11.tar.gz", "70c1ad0563ac300ec0cbd3c1b3420bd07d33cf06d124d746b2e349b0f0f53769"),
)

# Install unsatisfied or updated dependencies:
# We added `, isolate=true` as otherwise, it would segfault when closing `OpenBLAS32`,
# probably because it is conflicting with Julia openblas.
unsatisfied = any(!satisfied(p; verbose=verbose, isolate=true) for p in products)
dl_info = choose_download(download_info, platform_key_abi())
if dl_info === nothing && unsatisfied
    # If we don't have a compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform (\"$(Sys.MACHINE)\", parsed as \"$(triplet(platform_key_abi()))\") is not supported by this package!")
end

# If we have a download, and we are unsatisfied (or the version we're
# trying to install is not itself installed) then load it up!
if unsatisfied || !isinstalled(dl_info...; prefix=prefix)
    # Download and install binaries
    install(dl_info...; prefix=prefix, force=true, verbose=verbose)
end

ext: jl


text:
struct JacobiPreconditioner{Tv, Ti} <: AbstractExtendablePreconditioner{Tv,Ti}
    extmatrix::ExtendableSparseMatrix{Tv,Ti}
    invdiag::Array{Tv,1}
end

function update!(precon::JacobiPreconditioner)
    cscmatrix=precon.extmatrix.cscmatrix
    invdiag=precon.invdiag
    n=cscmatrix.n
    @inbounds for i=1:n
        invdiag[i]=1.0/cscmatrix[i,i]
    end
    precon
end

function JacobiPreconditioner(extmatrix::ExtendableSparseMatrix{Tv,Ti}) where {Tv,Ti}
    @assert size(extmatrix,1)==size(extmatrix,2)
    flush!(extmatrix)
    invdiag=Array{Tv,1}(undef,extmatrix.cscmatrix.n)
    precon=JacobiPreconditioner{Tv, Ti}(extmatrix,invdiag)
    update!(precon)
end

JacobiPreconditioner(cscmatrix::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}=JacobiPreconditioner(ExtendableSparseMatrix(cscmatrix))

function  LinearAlgebra.ldiv!(u::AbstractArray{T,1} where T, precon::JacobiPreconditioner, v::AbstractArray{T,1} where T)
    invdiag=precon.invdiag
    n=length(invdiag)
    @inbounds @simd for i=1:n
        u[i]=invdiag[i]*v[i]
    end
end


text:
using Main.Deldir
using Main.GeometricalPredicates

function fh(x, y)::Real
    return 1
end

function dc(p)
    return -0.3 + abs(0.7 - sqrt(sum(p .^ 2)))
end

function dca(p)
    return sqrt(sum(p .^ 2)) - 1
end

function dcb(p, x1, x2, y1, y2)
    return -min(min(min(-y1+p[2],y2-p[2], -x1+p[1],x2-p[1])))
end

struct Scaler
    scale::Any
    transx::Any
    transy::Any

    function Scaler(bbox::Array{Int,2})
        _scale = (abs(1 / (bbox[2] - bbox[1])))
        new(_scale, (0 - bbox[1] * _scale), (0 - bbox[1, 2] * _scale))
    end
end

function scale_p(p, scaler::Scaler)
    return [(p[1] * scaler.scale) + scaler.transx, (p[2] * scaler.scale) + scaler.transy]
end

function unscale_p(p, scaler::Scaler)
    return [(p[1] - scaler.transx) / scaler.scale, (p[2] - scaler.transy) / scaler.scale]
end

function scaled_point(unscaled_point::Point2D, scaler::Scaler)
    return Point(
        ((getx(unscaled_point) * scaler.scale) + scaler.transx),
        ((gety(unscaled_point) * scaler.scale) + scaler.transy),
    )
end

function unscaled_point(scaled_point::Point2D, scaler::Scaler)
    return Point(
        (getx(scaled_point) - scaler.transx) / scaler.scale,
        (gety(scaled_point) - scaler.transy) / scaler.scale,
    )
end

function unscale_y(p, scaler)
    return (p - scaler.transx) / scaler.scale
end

function unscale_x(p, scaler)
    return (p - scaler.transy) / scaler.scale
end

function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}) where {T}
    m, n = length(vy), length(vx)
    gx = reshape(repeat(vx, inner = m, outer = 1), m, n)
    gy = reshape(repeat(vy, inner = 1, outer = n), m, n)
    return gx, gy
end

function extractbox(bbox, h0::Real, h1::Real)
    v1 = bbox[1, 1]:h0:bbox[2, 1]
    v2 = bbox[1, 2]:h1:bbox[2, 2]
    return v1, v2
end

function calculateh1(h0::Real)
    return h0 * (sqrt(3) / 2)
end

function shiftevenrows!(x, h0::Real)
    x[2:2:end, :] = x[2:2:end, :] .+ h0 / 2
end

function triangles(del, summ)
    generators = Dict{Int64,Array{Int,1}}()
    for (index, row) in enumerate(eachrow(summ))
        generators[index] = []
    end
    for row in eachrow(del)
        generators[row[5]] = append!(generators[row[5]], row[6])
        generators[row[6]] = append!(generators[row[6]], row[5])
    end
    triangles = Array{Array{Int,1},1}()
    for k in keys(generators)
        for e in generators[k]
            i = findall(in(generators[e]), generators[k])
            common = generators[k][i]
            for c in common
                tri = sort([c, k, e])
                if !(tri in triangles)
                    push!(triangles, tri)
                end
            end
        end
    end
    return map(
        t -> [
            Point(summ[t[1], 1], summ[t[1], 2]),
            Point(summ[t[2], 1], summ[t[2], 2]),
            Point(summ[t[3], 1], summ[t[3], 2]),
        ],
        triangles,
    )
end

function vectorizededges(triangle)
        a = triangle[1]
        b = triangle[2]
        c = triangle[3]
        ab = [getx(a) gety(a); getx(b) gety(b)]
        ba = [getx(b) gety(b); getx(a) gety(a)]
        bc = [getx(b) gety(b); getx(c) gety(c)]
        cb = [getx(c) gety(c); getx(b) gety(b)]
        ac = [getx(a) gety(a); getx(c) gety(c)]
        ca = [getx(c) gety(c); getx(a) gety(a)]
        return ab, ba, bc, cb, ac, ca
end

function validedges(triangles, del, scaler, fd, geps)
    i = 1
    inside_edges = Array{Array{Float64,2},1}()
    outside_edges = Array{Array{Float64,2},1}()
    for triangle in triangles
        center = unscaled_point(centroid(Primitive(triangle[1], triangle[2], triangle[3])), scaler)
        i += 1
        a = triangle[1]
        b = triangle[2]
        c = triangle[3]
        edges = vectorizededges(triangle)
        if fd([getx(center), gety(center)]) > -geps
            push!(inside_edges, edges...)
        else
            push!(outside_edges, edges...)
        end
    end
    nedges = inside_edges
    yedges = outside_edges
    inside_edges = filter(edge -> !(edge in outside_edges), inside_edges)
    edges = [[r[1] r[2]; r[3] r[4]] for r in eachrow(del)]
    edges = filter(edge -> !(edge in inside_edges), edges)
    x = Array{Float64,1}()
    y = Array{Float64,1}()
    line_edges = Array{Line2D{Point2D},1}()
    for edge in edges
        push!(x, edge[1])
        push!(x, edge[2])
        push!(x, NaN)
        push!(y, edge[3])
        push!(y, edge[4])
        push!(y, NaN)
        push!(line_edges, Line(Point(edge[1], edge[3]), Point(edge[2], edge[4])))
    end
    return line_edges, x, y
end

function pointstoforces(edges, scaler, Fscale, pfix)
    bars = Array{Point2D,1}()
    barvec = Array{Array{Float64,1},1}()
    points_to_fvces = Dict{Point2D,Array{Float64,1}}()
    for edge in edges
        b = unscaled_point(getb(edge), scaler)
        a = unscaled_point(geta(edge), scaler)
        push!(
            bars,
            Point(getx(a) + ((getx(b) - getx(a)) / 2), gety(a) + ((gety(b) - gety(a)) / 2)),
        )
        push!(
            barvec,
            [
                getx(b) - getx(a)
                gety(b) - gety(a)
            ],
        )
        push!(points_to_fvces, geta(edge) => [0, 0])
        push!(points_to_fvces, getb(edge) => [0, 0])
    end
    L = [sqrt(sum(v_sum .^ 2)) for v_sum in barvec]
    hbars = [fh(getx(p), gety(p)) for p in bars]
    L0 = hbars * Fscale * sqrt(sum(L .^ 2) / sum(hbars .^ 2))
    Fvec = maximum.(L0 - L) ./ L .* barvec
    iterator = 1
    for edge in edges
        a = geta(edge)
        b = getb(edge)
        prev_a = points_to_fvces[a]
        prev_b = points_to_fvces[b]
        push!(points_to_fvces, a => prev_a + (-Fvec[iterator]))
        push!(points_to_fvces, b => prev_b + (Fvec[iterator]))
        iterator = iterator + 1
    end
    for p in eachrow(pfix)
        p = Point(p[1], p[2])
        delete!(points_to_fvces, p)
        push!(points_to_fvces, p => [0, 0])
    end
    return points_to_fvces
end

function finalpositions(points_to_fvces, scaler, deltat, fd, geps, deps, h0)
    new_p = Array{Point2D,1}()
    d_points = Array{Array{Float64,1},1}()
    for (point, force) in points_to_fvces
        p = unscaled_point(point, scaler)
        np = [getx(p), gety(p)] + deltat * force
        push!(new_p, Point2D(np[1], np[2]))
        d = fd(np)
        if d < -geps
            push!(d_points, force)
        end
    end
        final_p = Array{Array{Float64,1},1}()
    for p in new_p
        x = getx(p)
        y = gety(p)
        d = fd([x, y])
        if d > 0
            dgradx = (fd([x + deps, y]) - d) / deps
            dgrady = (fd([x, y + deps]) - d) / deps
            res = [x, y] - [d * dgradx, d * dgrady]
            push!(final_p, res)
        else
            push!(final_p, [x, y])
        end
    end
    final_p = unique(final_p)
    final_p = map(p -> scale_p([p[1], p[2]], scaler), final_p)
    final_p = transpose(reshape(vcat(final_p...), 2, length(final_p)))
    return final_p, move_index(d_points, deltat, h0)
end

function move_index(d_points, deltat, h0)
    d = map(row -> sum(deltat * row .^ 2), d_points)
    push!(d, 0)
    return maximum(sqrt.(d) / h0)
end

function generate(fd, fh, bbox, h0, pfix = [], ttol = 0.1, geps = 0.001 * h0, Fscale = 1.2, dptol = 0.001, deltat = 0.2, deps = sqrt(eps(Float64)) * h0)
    h1 = calculateh1(h0)
    pold = Inf
    scaler = Scaler(bbox)
    v1, v2 = extractbox(bbox, h0, h1)
    x, y = meshgrid(v1, v2)
    shiftevenrows!(x, h0)
    p = [x[:] y[:]]
    p = [vcat(row) for row in eachrow(p) if fd(row) < -geps]
    r0 = [1 ./ fh(row...) .^ 2 for row in p]
    r0_max = maximum(r0 ./ maximum(r0))
    p = [vcat(scale_p(row, scaler)) for row in p if (rand(Float64, size(p)))[1] < r0_max]
    p = transpose(reshape(vcat(p...), 2, length(p)))
    pfix = [vcat(scale_p(row, scaler)) for row in eachrow(pfix)]
    pfix = transpose(reshape(vcat(pfix...), 2, length(pfix)))
    while true
        del, vor, summ = deldir(p[:, 1], p[:, 2])
        trigs = triangles(del, summ)
        line_edges, x, y = validedges(trigs, del, scaler, fd, geps)
        points_to_fvces = pointstoforces(line_edges, scaler, Fscale, pfix)
        final_p, move_index = finalpositions(points_to_fvces, scaler, deltat, fd, geps, deps, h0)
        if move_index < dptol
            break
        else
            p = final_p
        end
    end
    x = map(x -> unscale_x(x, scaler), x)
    y = map(y -> unscale_y(y, scaler), y)
    return x, y
end

max_stars_repo_name: jstarczewski/julia-num

max_issues_repo_name: jstarczewski/julia-num

max_forks_repo_name: jstarczewski/julia-num


text:
using SDEs
using Compat, BenchmarkTools
include("../test/setup.jl")
include("../test/test_utils.jl")

θ=0.5
η=0.5

println(prob, ", dt:", t[2]-t[1], ", θ:",θ, ", η:",η)
print("EM proto  :")
@btime x, W = integrate(prob, x0, t, SDEs.EMproto(), 1)
print("EM        :")
@btime x, W = integrate(prob, x0, t, EM(), 1)
print("Diffeq EM :")
@btime x, W = integrate(prob, x0, t, DEQEM(), 1)
print("PCE proto :")
@btime x, W = integrate(prob, x0, t, SDEs.PCEproto(), 1; θ=θ, η=η)
print("PCE       :")
@btime x, W = integrate(prob, x0, t, PCE(), 1; θ=θ, η=η)

max_stars_repo_name: onoderat/SDEs

max_issues_repo_name: onoderat/SDEs

max_forks_repo_name: onoderat/SDEs


text:
using ModelKit, LinearAlgebra, Test
import SymEngine
const SE = SymEngine

@testset "ModelKit.jl" begin

    @testset "Variables" begin
        @var a b x[1:2] y[1:2, 1:3]

        @test a isa Variable
        @test b isa Variable
        @test x isa Vector{Variable}
        @test length(x) == 2
        @test y isa Matrix{Variable}
        @test size(y) == (2, 3)

        @var c, d
        @test c isa Variable
        @test d isa Variable

        let
            c2 = @unique_var c
            @test c != c2
        end

        f = a + b
        @test variables(f) == Set([a, b])
        @test nvariables(f) == 2
        g = x[1]
        @test variables([f, g]) == Set([a, b, x[1]])
        @test nvariables([f, g]) == 3
    end

    @testset "Subs" begin
        @var x y w z u
        f = x^2 * (x + y * w)

        @test subs(f, x => z) == z^2 * (z + w * y)
        @test subs([f], x => z) == [z^2 * (z + w * y)]
        @test subs(f, [x, y] => [z^2, z + 2]) == z^4 * (w * (2 + z) + z^2)
        @test subs(f, [x, y] => [z^2, z + 2], w => u) == z^4 *
                                                         (u * (2 + z) + z^2)
        @test subs(f, x => z^2, y => 3, w => u) == z^4 * (3 * u + z^2)
    end

    @testset "Evaluation" begin
        @var x y w
        f = x^2 * (x + y * w)
        @test evaluate([f, f], [x, y, w] => [2, 3, -5]) == [-52, -52]
        @test [f, 2f]([x, y, w] => [2, 3, -5]) == [-52, -104]
    end

    @testset "Linear Algebra" begin
        @var x[1:2, 1:2]
        @test det(x) == -x[2, 1] * x[1, 2] + x[2, 2] * x[1, 1]
    end

    @testset "Differentation" begin
        @var x y

        f = x^2 + y^2
        g = x^3 + y^3

        @test differentiate(f, x) == 2x
        @test differentiate(f, [x, y]) == [2x, 2y]

        @test differentiate([f, g], x) == [2x, 3 * x^2]
        @test differentiate([f, g], [x, y]) == [2x 2y; 3 * x^2 3 * y^2]
    end

    @testset "Expand" begin
        @var x y
        @test expand((x + y)^2) == 2 * x * y + x^2 + y^2
    end

    @testset "Modeling" begin
        @testset "Bottleneck" begin
            @var x y z
            f = [(0.3 * x^2 + 0.5z + 0.3x + 1.2 * y^2 - 1.1)^2 +
                 (0.7 * (y - 0.5x)^2 + y + 1.2 * z^2 - 1)^2 - 0.3]

            I = let x = [x, y, z]
                n, m = length(x), length(f)
                @unique_var y[1:n] v[1:m] w[1:m]
                J = differentiate(f, x)
                f′ = f(x => y)
                J′ = J(x => y)
                Nx = (x - y) - J' * v
                Ny = (x - y) - J′' * w
                System([f; f′; Nx; Ny], [x; y; v; w])
            end
            @test I isa System
            @test size(I) == (8, 8)
        end

        @testset "Steiner" begin
            @var x[1:2] a[1:5] c[1:6] y[1:2, 1:5]

            #tangential conics
            f = sum([a; 1] .* monomials(x, 2))
            ∇ = differentiate(f, x)
            #5 conics
            g = sum(c .* monomials(x, 2))
            ∇_2 = differentiate(g, x)
            #the general system
            #f_a_0 is tangent to g_b₀ at x₀
            function Incidence(f, a₀, g, b₀, x₀)
                fᵢ = f(x => x₀, a => a₀)
                ∇ᵢ = ∇(x => x₀, a => a₀)
                Cᵢ = g(x => x₀, c => b₀)
                ∇_Cᵢ = ∇_2(x => x₀, c => b₀)

                [fᵢ; Cᵢ; det([∇ᵢ ∇_Cᵢ])]
            end
            @var v[1:6, 1:5]
            I = vcat(map(i -> Incidence(f, a, g, v[:, i], y[:, i]), 1:5)...)
            F = System(I, [a; vec(y)], vec(v))
            @test size(F) == (15, 15)
        end

        @testset "Reach plane curve" begin
            @var x y
            f = (x^3 - x * y^2 + y + 1)^2 * (x^2 + y^2 - 1) + y^2 - 5
            ∇ = differentiate(f, [x; y]) # the gradient
            H = differentiate(∇, [x; y]) # the Hessian

            g = ∇ ⋅ ∇
            v = [-∇[2]; ∇[1]]
            h = v' * H * v
            dg = differentiate(g, [x; y])
            dh = differentiate(h, [x; y])

            ∇σ = g .* dh - ((3 / 2) * h) .* dg

            F = System([v ⋅ ∇σ; f], [x, y])
            @test size(F) == (2, 2)
        end
    end

    @testset "System" begin
        @var x y a b
        f = [(x + y)^3 + x^2 + x + 5y + 3a, 2 * x^2 + b]
        F = System(f, [x, y], [b, a])

        show_F = """
        System of length 2
         2 variables: x, y
         2 parameters: b, a

         3*a + x + 5*y + x^2 + (x + y)^3
         b + 2*x^2"""
        @test sprint(show, F) == show_F

        T = compile(F)
        F2 = ModelKit.interpret(T)
        @test F == F2
        @test sprint(show, T) == "Compiled: $show_F"
        @test size(T) == size(F) == (2, 2)
    end

    @testset "Homotopy" begin
        @var x y z t

        h = [x^2 + y + z + 2t, 4 * x^2 * z^2 * y + 4z - 6x * y * z^2]
        H = Homotopy(h, [x, y, z], t)

        show_H = """
        Homotopy in t of length 2
         3 variables: x, y, z

         2*t + y + z + x^2
         4*z - 6*x*y*z^2 + 4*x^2*y*z^2"""

        @test sprint(show, H) == show_H

        T = compile(H)
        H2 = ModelKit.interpret(T)
        @test H == H2
        @test sprint(show, T) == "Compiled: $show_H"
        @test size(T) == size(H) == (2, 3)
    end

    @testset "Codegen helpers" begin
        @test ModelKit.sqr(3 + 2im) == (3 + 2im)^2
        @test ModelKit.sqr(3) == 3^2
    end

    @testset "Instructions" begin

        @testset "Higher order pow diff" begin
            for d in [2, 5]
                @var x
                f = x^d
                list, _ = ModelKit.instruction_list([f])

                diff_map = Dict()
                diff_map[(:x, 1)] = :x1
                diff_map[(:x, 2)] = :x2
                diff_map[(:x, 3)] = :x3

                @eval ModelKit begin
                    function __diff_4_pow(x, x1, x2, x3, t)
                        $(ModelKit.to_expr(ModelKit.univariate_diff!(
                            list,
                            4,
                            diff_map,
                        )))
                    end
                end

                SE.@vars x t
                SE.@funs u

                exp1 = SE.expand(SE.subs(
                    SE.diff(u(t)^d, t, 4),
                    SE.diff(u(t), t, 4) => 0,
                ) / factorial(4),)

                u1 = SE.diff(u(t), t)
                u2 = SE.diff(u(t), t, 2) / 2
                u3 = SE.diff(u(t), t, 3) / 6
                exp2 = SE.expand(ModelKit.__diff_4_pow(u(t), u1, u2, u3, t))
                @test exp1 == exp2
            end
        end

        @testset "Higher order mul" begin
            @var x y
            f = x * y
            list, _ = ModelKit.instruction_list([f])

            diff_map = Dict()
            diff_map[(:x, 1)] = :x1
            diff_map[(:x, 2)] = :x2
            diff_map[(:x, 3)] = :x3
            diff_map[(:y, 1)] = :y1
            diff_map[(:y, 2)] = :y2
            diff_map[(:y, 3)] = :y3

            @eval ModelKit begin
                function __diff_4_mul__(x, y, t)
                    x1 = Main.SE.diff(x, t)
                    x2 = Main.SE.diff(x, t, 2) / 2
                    x3 = Main.SE.diff(x, t, 3) / 6
                    y1 = Main.SE.diff(y, t)
                    y2 = Main.SE.diff(y, t, 2) / 2
                    y3 = Main.SE.diff(y, t, 3) / 6
                    $(ModelKit.to_expr(ModelKit.univariate_diff!(
                        list,
                        4,
                        diff_map,
                    )))
                end
            end

            SE.@funs u v
            SE.@vars t

            exp1 = SE.expand(SE.subs(
                SE.diff(u(t) * v(t), t, 4),
                SE.diff(u(t), t, 4) => 0,
                SE.diff(v(t), t, 4) => 0,
            ) / factorial(4),)
            exp2 = SE.expand(ModelKit.__diff_4_mul__(u(t), v(t), t))
            @test exp1 == exp2
        end

        @testset "Higher order plus" begin
            @var x y
            f = x + y
            list, _ = ModelKit.instruction_list([f])

            diff_map = Dict()
            diff_map[(:x, 1)] = :x1
            diff_map[(:x, 2)] = :x2
            diff_map[(:x, 3)] = :x3
            diff_map[(:y, 1)] = :y1
            diff_map[(:y, 2)] = :y2
            diff_map[(:y, 3)] = :y3

            @eval ModelKit begin
                function __diff_3_plus__(x, y, t)
                    x1 = Main.SE.diff(x, t)
                    x2 = Main.SE.diff(x, t, 2) / 2
                    x3 = Main.SE.diff(x, t, 3) / 6
                    y1 = Main.SE.diff(y, t)
                    y2 = Main.SE.diff(y, t, 2) / 2
                    y3 = Main.SE.diff(y, t, 3) / 6
                    $(ModelKit.to_expr(ModelKit.univariate_diff!(
                        list,
                        3,
                        diff_map,
                    )))
                end
            end

            SE.@funs u v
            SE.@vars t

            exp1 = SE.expand(diff(u(t) + v(t), t, 3) / factorial(3))
            exp2 = SE.expand(ModelKit.__diff_3_plus__(u(t), v(t), t))
            @test exp1 == exp2
        end
    end

end


text:
# export StackedTimeSolverData
"""
    StackedTimeSolverData <: AbstractSolverData

The data structure used in the stacked time algorithm.

**TODO** Add all the details here.
"""
struct StackedTimeSolverData <: AbstractSolverData
    "Number of time periods"
    NT::Int
    "Number of variables"
    NV::Int
    "Number of shocks"
    NS::Int
    "Number of unknowns = variables + shocks + aux variables"
    NU::Int
    "Final condition type: `fcgiven` or `fclevel` or `fcslope`."
    FC::FCType
    "List of vectors for the time indexes of blocks of the matrix."
    TT::Vector{Vector{Int}}
    "List of vectors for the row indexes of each block."
    II::Vector{Vector{Int}}
    "List of vectors for the column indexes of each block."
    JJ::Vector{Vector{Int}}
    "Inverse lookup for the column indexes."
    BI::Vector{Vector{Int}}
    "Jacobian matrix cache."
    J::SparseMatrixCSC{Float64,Int}
    "Correction to J in case FC == fcslope. See description of algorithm."
    CiSc::SparseMatrixCSC{Float64,Int}   # Correction to J to account for FC in the case of fcslope
    "The steady state data, used for FC ∈ (fclevel, fcslope)."
    SS::AbstractArray{Float64,2}
    "Mask of variables with constant growth rate in the steady state"
    log_mask::AbstractVector{Bool}
    "Mask of variables with constant change in the steady state"
    lin_mask::AbstractVector{Bool}
    "The `evaldata` from the model, used to evaluate RJ and R!"
    evaldata::ModelBaseEcon.AbstractModelEvaluationData
    "Keep track of which variables are set by exogenous constraints."
    exog_mask::AbstractVector{Bool}
    "Keep track of which variables are set by final conditions."
    fc_mask::AbstractVector{Bool}
    "Keep track of which variables are \"active\" for the solver."
    solve_mask::AbstractVector{Bool}
    "Cache the LU decomposition of the active part of J"
    luJ::Ref{Any}
end

"""
    update_plan!(sd::StackedTimeSolverData, model, plan; changed=false)

Implementation for the Stacked Time algorithm. Plan must have the same range as
the original plan.

By default the data structure is updated only if an actual change in the plan is
detected. Setting the `changed` flag to `true` forces the update even if the
plan seems unchanged. This is necessary only in rare circumstances.

"""
function update_plan!(sd::StackedTimeSolverData, m::Model, p::Plan; changed=false)
    if sd.NT != length(p.range)
        error("Unable to update using a simulation plan of different length.")
    end

    # LinearIndices used for indexing the columns of the global matrix
    LI = LinearIndices((p.range, 1:sd.NU))

    sim = m.maxlag + 1:sd.NT - m.maxlead
    NTFC = m.maxlead

    # Assume initial conditions are set correctly to exogenous in the constructor
    # We only update the masks during the simulation periods
    # unknowns = m.allvars
    foo = falses(sd.NU)
    for t in sim
        # s = p[t]
        # si = indexin(s, unknowns)
        si = p[t, Val(:inds)]
        fill!(foo, false)
        foo[si] .= true
        if !all(foo .== sd.exog_mask[LI[t,:]])
            sd.exog_mask[LI[t,:]] .= foo
            changed = true
        end
    end

    # Update the solve_mask array
    if changed
        @. sd.solve_mask = !(sd.fc_mask | sd.exog_mask)
        @assert !any(sd.exog_mask .& sd.fc_mask)
        @assert sum(sd.solve_mask) == size(sd.J, 1)
    end

    # Update the Jacobian correction matrix, if exogenous plan changed
    if changed && (sd.FC == fcslope)
        # Matrix C is composed of blocks looking like this:
        #    1  0  0  0
        #   -1  1  0  0
        #    0 -1  1  0
        #    0  0 -1  1
        # The corresponding rows in Sc look like this
        #   0 0 0 0 0 -1
        #   0 0 0 0 0  0
        #   0 0 0 0 0  0
        #   0 0 0 0 0  0
        # The C^{-1} * Sc matrix then looks like this
        #   0 0 0 0 0 -1
        #   0 0 0 0 0 -1
        #   0 0 0 0 0 -1
        #   0 0 0 0 0 -1
        # The number of rows equals NTFC (the number of final condition periods)
        # The column with the -1 is the one corresponding to the last simulation period of the variable.

        # # foo is a mask array with foo[i] = true if i is a column  index for the last simulation period of a variable and foo[i] = false otherwise.
        # foo = falses(sd.NU * sd.NT)                 # all false
        # foo[LI[last(sim),1:sd.NV]] .= true          # true at last sim period of variables
        # foo[LI[last(sim),sd.NV + sd.NS + 1:sd.NU]] .= true  # true at last sim period of aux variables
        # # no final conditions for shocks - their final condition is always 0 (the expectation of the shock).

        # # The column indices - findall(foo) returns the indexes of the elements that are true.
        # #   then we repeat each NTFC times, since this is the number of rows we need
        # JJfoo = repeat(findall(foo[sd.solve_mask]), inner=NTFC)

        # The row indices
        #   the block for variable 1 gets row indexes 1 .. NTFC
        #   the block for variable 2 gets row indexes NTFC+1 .. 2*NTFC
        #   the block for variable vi gets row indexes (vi-1)*NTFC+1 .. vi*NTFC
        # IIfoo = vcat((NTFC * (vi - 1) .+ (1:NTFC) for vi in Iterators.flatten((1:sd.NV, sd.NV + sd.NS + 1:sd.NU)) if sd.solve_mask[LI[last(sim),vi]])...)
        IIfoo = Int[]
        JJfoo = Int[]
        VVfoo = Float64[]
        for vi in Iterators.flatten((1:sd.NV, (sd.NV + sd.NS + 1):sd.NU))
            colind = LI[last(sim), vi]
            if sd.solve_mask[colind]
                # if the variable is exogenous in the last period, then Sc doesn't have a column for it.
                append!(IIfoo, ((vi - 1) * NTFC + 1):(vi * NTFC))
                append!(VVfoo, fill(-1.0, NTFC))
                foo = falses(sd.NU * sd.NT)
                foo[colind] = true
                append!(JJfoo, repeat(findall(foo[sd.solve_mask]), NTFC))
            end
        end
        # Construct the sparse matrix
        sd.CiSc .= sparse(IIfoo, JJfoo, VVfoo, NTFC * sd.NU, size(sd.J, 1))
        # sd.CiSc .= sparse(IIfoo, JJfoo, fill(Float64(-1), length(JJfoo)), NTFC * sd.NU, size(sd.J, 1))
        # sd.CiSc.nzval .= -1
    end
    # Update the Jacobian correction matrix, if exogenous plan changed
    if changed && (sd.FC == fcnatural)
        # Matrix C is composed of blocks looking like this:
        #    -1  0  0  0  0
        #     2 -1  1  0  0
        #    -1  2 -1  0  0
        #     0 -1  2 -1  0
        #     0  0 -1  2 -1
        # The corresponding rows in Sc look like this
        #     0 ..  0 -1  2
        #     0 ..  0  0 -1
        #     0 ..  0  0  0
        #     0 ..  0  0  0
        #     0 ..  0  0  0
        # The C^{-1} * Sc matrix then looks like this
        #     0 ..  0  1 -2
        #     0 ..  0  2 -3
        #     0 ..  0  3 -4
        #     0 ..  0  4 -5
        #     0 ..  0  5 -6

        # The number of rows equals NTFC (the number of final condition periods)
        # The column with the -2,-3,-4... is the one corresponding to the last simulation period of the variable.
        # The column with the 1,2,3... is the one corresponding to the period before the last simulation period of the variable.

        # The row indices
        #   the block for variable 1 gets row indexes 1 .. NTFC
        #   the block for variable 2 gets row indexes NTFC+1 .. 2*NTFC
        #   the block for variable vi gets row indexes (vi-1)*NTFC+1 .. vi*NTFC
        # IIfoo = vcat((NTFC * (vi - 1) .+ (1:NTFC) for vi in Iterators.flatten((1:sd.NV, sd.NV + sd.NS + 1:sd.NU)) if sd.solve_mask[LI[last(sim),vi]])...)
        IIfoo = Int[]
        JJfoo = Int[]
        VVfoo = Float64[]
        for vi in Iterators.flatten((1:sd.NV, (sd.NV + sd.NS + 1):sd.NU))
            let col1 = LI[last(sim) - 1, vi]
                if sd.solve_mask[col1]
                    append!(IIfoo, ((vi - 1) * NTFC + 1):(vi * NTFC))
                    append!(VVfoo, 1:NTFC)
                    foo = falses(sd.NU * sd.NT)
                    foo[col1] = true
                    append!(JJfoo, repeat(findall(foo[sd.solve_mask]), NTFC))
                end
            end
            let col2 = LI[last(sim), vi]
                if sd.solve_mask[col2]
                # if the variable is exogenous in the last period, then Sc doesn't have a column for it.
                    append!(IIfoo, ((vi - 1) * NTFC + 1):(vi * NTFC))
                    append!(VVfoo, -1.0 .- (1:NTFC))   # -2, -3, ..., -NTFC-1
                    foo = falses(sd.NU * sd.NT)
                    foo[col2] = true
                    append!(JJfoo, repeat(findall(foo[sd.solve_mask]), NTFC))
                end
            end
        end
        # Construct the sparse matrix
        sd.CiSc .= sparse(IIfoo, JJfoo, VVfoo, NTFC * sd.NU, size(sd.J, 1))
    end

    if changed
        # cached lu is no longer valid, since active columns have changed
        sd.luJ[] = nothing
    end

    return sd
end

"""
    make_BI(J, II)

Prepares the `BI` array for the solver data. Called from the constructor of
`StackedTimeSolverData`.

"""
function make_BI(JMAT::SparseMatrixCSC{Float64,Int}, II::AbstractVector{<:AbstractVector{Int}})
    # Computes the set of indexes in JMAT.nzval corresponding to blocks of equations in II

    # Bar is the inverse map of II, i.e. Bar[j] = i <=> II[i] contains j
    local Bar = zeros(Int, size(JMAT, 1))
    for i in axes(II, 1)
        Bar[II[i]] .= i
    end

    # start with empty lists in BI
    local BI = [[] for _ in II]

    # iterate the non-zero elements of JMAT and record their row numbers in BI
    local rows = rowvals(JMAT)
    local n = size(JMAT, 2)
    for i = 1:n
        for j in nzrange(JMAT, i)
            row = rows[j]
            push!(BI[ Bar[row] ], j)
        end
    end
    return BI
end

"""
    StackedTimeSolverData(model, plan, fctype)

"""
function StackedTimeSolverData(m::Model, p::Plan, fctype::FCType)

    # NOTE: we do not verify the steady state, just make sure it's been assigned
    if fctype ∈ (fclevel, fcslope) && !issssolved(m)
        throw(ArgumentError("Steady state must be solved for `$(fctype)`."))
    end

    NT = length(p.range)
    init = 1:m.maxlag
    term = NT - m.maxlead + 1:NT
    sim = m.maxlag + 1:NT - m.maxlead
    NTFC = length(term)
    NTSIM = length(sim)
    NTIC = length(init)

    if (fctype == fcnatural) && (length(sim) + length(init) < 2)
        throw(ArgumentError("Simulation must include at least 2 periods for `fcnatural`."))
    end

    # @assert NTIC+NTSIM+NTFC == NT
    # @assert isempty((collect(init)∩collect(sim))∪(collect(init)∩collect(sim))∪(collect(sim)∩collect(term)))

    unknowns = m.allvars
    nunknowns = length(unknowns)

    nvars = length(m.variables)
    nshks = length(m.shocks)
    nauxs = length(m.auxvars)

    equations = m.alleqns
    nequations = length(equations)

    # LinearIndices used for indexing the columns of the global matrix
    LI = LinearIndices((p.range, 1:nunknowns))

    # Initialize empty arrays
    TT = Vector{Int}[] # the time indexes of the block, used when updating values in global_RJ
    II = Vector{Int}[] # the row indexes
    JJ = Vector{Int}[] # the column indexes

    # Prep the Jacobian matrix
    neq = 0 # running counter of equations added to matrix
    # Model equations are the same for each sim period, just shifted according to t
    Jblock = [ti + NT * (vi - 1) for eqn in equations for (ti, vi) in eqn.vinds]
    Iblock = [i for (i, eqn) in enumerate(equations) for _ in eqn.vinds]
    Tblock = -m.maxlag:m.maxlead
    @timer for t in sim
        push!(TT, t .+ Tblock)
        push!(II, neq .+ Iblock)
        push!(JJ, t .+ Jblock)
        neq += nequations
    end
    # Construct the
    @timer begin
        I = vcat(II...)
        J = vcat(JJ...)
        JMAT = sparse(I, J, fill(NaN64, size(I)), NTSIM * nequations, NT * nunknowns)
    end

    # We no longer need the exact indexes of all non-zero entires in the Jacobian matrix.
    # We do however need the set of equation indexes for each sim period
    @timer foreach(unique! ∘ sort!, II)

    # BI holds the indexes in JMAT.nzval for each block of equations
    BI = make_BI(JMAT, II)  # same as the two lines above, but faster

    # We also need the times of the final conditions
    push!(TT, term)

    #
    exog_mask = falses(nunknowns * NT)
    # Initial conditions are set as exogenous values
    exog_mask[vec(LI[init,:])] .= true
    # The exogenous values during sim are set in update_plan!() call below.

    # Final conditions are complicated
    fc_mask = falses(nunknowns * NT)
    fc_mask[vec(LI[term,:])] .= true

    # The solve_mask is redundant. We pre-compute and store it for speed
    solve_mask = .!(fc_mask .| exog_mask)

    sd = StackedTimeSolverData(NT, nvars, nshks, nunknowns, fctype, TT, II, JJ, BI, JMAT,
                             sparse([], [], Float64[], NTFC * nunknowns, size(JMAT, 1)),
                             log_transform!(steadystatearray(m, p), m, :to_internal),
                             islog.(m.allvars), islin.(m.allvars),
                             m.evaldata, exog_mask, fc_mask, solve_mask,
                             Ref{Any}(nothing))

    return update_plan!(sd, m, p; changed=true)
end

function assign_exog_data!(x::AbstractArray{Float64,2}, exog::AbstractArray{Float64,2}, sd::StackedTimeSolverData)
    # @assert size(x,1) == size(exog,1) == sd.NT
    x[sd.exog_mask] = exog[sd.exog_mask]
    assign_final_condition!(x, exog, sd, Val(sd.FC))
    return x
end

@inline function assign_final_condition!(x::AbstractArray{Float64,2}, exog::AbstractArray{Float64,2}, sd::StackedTimeSolverData, ::Val{fcgiven})
    # We assume that exog data has been transformed for log-variables
    x[sd.TT[end],:] = exog[sd.TT[end],:]
    return x
end

@inline function assign_final_condition!(x::AbstractArray{Float64,2}, ::AbstractArray{Float64,2}, sd::StackedTimeSolverData, ::Val{fclevel})
    # We assume that SS data has been transformed for log-variables
    x[sd.TT[end],:] = sd.SS[sd.TT[end],:]
    return x
end

@inline function assign_final_condition!(x::AbstractArray{Float64,2}, ::AbstractArray{Float64,2}, sd::StackedTimeSolverData, ::Val{fcrate})
    # We assume that SS data has been transformed for log-variables
    slp_vars = sd.lin_mask .| sd.log_mask
    for t in sd.TT[end]
        x[t,:] = x[t - 1,:]
        @. x[t, slp_vars] += sd.SS[t, slp_vars] - sd.SS[t - 1, slp_vars]
    end
    return x
end

function assign_final_condition!(x::AbstractArray{Float64,2}, ::AbstractArray{Float64,2}, sd::StackedTimeSolverData, ::Val{fcnatural})
    term_Ts = sd.TT[end]
    if isempty(term_Ts)
        return x
    end
    last_T = term_Ts[1] - 1
    slp_vars = sd.lin_mask .| sd.log_mask
    SLP = zeros(size(x, 2))
    @. SLP[slp_vars] = x[last_T,slp_vars] - x[last_T - 1, slp_vars]
    for t = term_Ts
        x[t,:] .= x[t - 1,:]
        @. x[t, slp_vars] += SLP[slp_vars]
    end
    return x
end

@inline assign_final_condition!(x, e, sd::StackedTimeSolverData) = assign_final_condition!(x, e, sd, Val(sd.FC))

function global_R!(res::AbstractArray{Float64,1}, point::AbstractArray{Float64}, exog_data::AbstractArray{Float64}, sd::StackedTimeSolverData)
    point = reshape(point, sd.NT, sd.NU)
    exog_data = reshape(point, sd.NT, sd.NU)
    @timer "global_R!" begin
        for (ii, tt) in zip(sd.II, sd.TT)
            eval_R!(view(res, ii), point[tt,:], sd.evaldata)
        end
    end
    return res
end

#= disable - this is not necessary with log-transformed variables

@inline update_CiSc!(x::AbstractArray{Float64,2}, sd::StackedTimeSolverData) = any(sd.log_mask) ? update_CiSc!(x, sd, Val(sd.FC)) : nothing

@inline update_CiSc!(x, sd, ::Val{fcgiven}) = nothing
@inline update_CiSc!(x, sd, ::Val{fclevel}) = nothing
function update_CiSc!(x, sd::StackedTimeSolverData, ::Val{fcslope})
    # LinearIndices used for indexing the columns of the global matrix
    local LI = LinearIndices((1:sd.NT, 1:sd.NU))
    # The last simulation period
    local term_Ts = sd.TT[end]
    if isempty(term_Ts)
        return
    end
    local last_T = term_Ts[1] - 1
    local NTFC = length(term_Ts)

    # @info "last_T = $(last_T), NTFC = $(NTFC)"

    # Iterate over the non-zero values of CiSc and update
    # values as necessary
    local rows = rowvals(sd.CiSc)
    local vals = nonzeros(sd.CiSc)
    local ncols = size(sd.CiSc, 2)
    for col = 1:ncols
        nzr = nzrange(sd.CiSc, col)
        if isempty(nzr)
            continue
        end
        # Compute variable and time index from column index
        var, tm = divrem(LI[sd.solve_mask][col] - 1, sd.NT)  .+ 1

        # @info "var = $(var), tm = $(tm)"

        if !sd.log_mask[var]
            continue
        end
        if tm != last_T
            error("Last simulation times don't match")
        end
        for j = nzr
            row = rows[j]
            # update vals[j]
            local T = rem(row - 1, NTFC) + 1 + last_T

            # @info "Updating $((row, col)) with T=$(T), x[T, var] = $(x[T, var]), x[last_T, var] = $(x[last_T, var])"

            vals[j] = - x[T, var] / x[last_T, var]
        end
    end
    return nothing
end

function update_CiSc!(x, sd, ::Val{fcnatural})
    # LinearIndices used for indexing the columns of the global matrix
    local LI = LinearIndices((1:sd.NT, 1:sd.NU))
    # The last simulation period
    local term_Ts = sd.TT[end]
    if isempty(term_Ts)
        return
    end
    local last_T = term_Ts[1] - 1
    local NTFC = length(term_Ts)

    # @info "last_T = $(last_T), NTFC = $(NTFC)"

    # Iterate over the non-zero values of CiSc and update
    # values as necessary
    local rows = rowvals(sd.CiSc)
    local vals = nonzeros(sd.CiSc)
    local ncols = size(sd.CiSc, 2)
    for col = 1:ncols
        nzr = nzrange(sd.CiSc, col)
        if isempty(nzr)
            continue
        end
        # Compute variable and time index from column index
        var, tm = divrem(LI[sd.solve_mask][col] - 1, sd.NT)  .+ 1

        # @info "var = $(var), tm = $(tm)"

        if !sd.log_mask[var]
            continue
        end
        if tm == last_T
            s = -1
            c = 2
        elseif tm == last_T - 1
            s = +1
            c = 1
        else
            error("Last simulation times don't match")
        end
        for j = nzr
            row = rows[j]
            # update vals[j]
            local T = rem(row - 1, NTFC) + 1 + last_T

            # @info "Updating $((row, col)) with T=$(T), x[T, var] = $(x[T, var]), x[last_T, var] = $(x[last_T, var])"

            vals[j] = s * (c + T - last_T - 1 ) * x[T, var] / x[tm, var]
        end
    end
    return nothing
end =#

function global_RJ(point::AbstractArray{Float64}, exog_data::AbstractArray{Float64}, sd::StackedTimeSolverData)
    nunknowns = sd.NU
    point = reshape(point, sd.NT, nunknowns)
    exog_data = reshape(exog_data, sd.NT, nunknowns)
    JAC = sd.J
    RES = Vector{Float64}(undef, size(JAC, 1))
    # Model equations @ [1] to [end-1]
    haveJ = isa(sd.evaldata, ModelBaseEcon.LinearizedModelEvaluationData) && !any(isnan.(JAC.nzval))
    haveLU = sd.luJ[] !== nothing
    if haveJ
        # update only RES
        @timer "globalRJ/evalR!" for i = 1:length(sd.BI)
            eval_R!(view(RES, sd.II[i]), point[sd.TT[i],:], sd.evaldata)
        end
    else
        # update both RES and JAC
        @timer "globalRJ/evalRJ" for i = 1:length(sd.BI)
            R, J = eval_RJ(point[sd.TT[i],:], sd.evaldata)
            RES[sd.II[i]] .= R
            JAC.nzval[sd.BI[i]] .= J.nzval
        end
        # update_CiSc!(point, sd)
        haveLU = false
    end
    if !haveLU
        if nnz(sd.CiSc) > 0
            JJ = JAC[:, sd.solve_mask] - JAC[:, sd.fc_mask] * sd.CiSc
        else
            JJ = JAC[:, sd.solve_mask]
        end
        try
            # compute lu decomposition of the active part of J and cache it.
            @timer "LU decomposition" sd.luJ[] = lu(JJ)
        catch e
            if e isa SingularException
                @error("The system is underdetermined with the given set of equations and final conditions.")
            end
            rethrow()
            # return RES, deepcopy(JJ)
        end
    end
    return RES, sd.luJ[]
end

function assign_update_step!(x::AbstractArray{Float64}, λ::Float64, Δx::AbstractArray{Float64}, sd::StackedTimeSolverData)
    x[sd.solve_mask] .+= λ .* Δx
    if nnz(sd.CiSc) > 0
        x[sd.fc_mask] .-= λ .* (sd.CiSc * Δx)
        # if sd.FC == fcnatural && any(sd.log_mask) && !isempty(sd.TT[end])
        #     for t = sd.TT[end]
        #         @. x[t, sd.log_mask] = x[t - 1, sd.log_mask]^2 / x[t - 2,sd.log_mask]
        #     end
        # end
        # assign_final_condition!(x, zeros(0,0), sd, Val(fcslope))
    end
    return x
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 12333bc8-68ca-11eb-2e89-ab9b76b0595b
begin
	using Plots
	using PlutoUI
end

# ╔═╡ f61ae738-68c9-11eb-03c1-efb38291b3ae
md"""
### Plots
"""

# ╔═╡ 27b37bc0-68ca-11eb-0c16-914e5a5ec540
PlutoUI.TableOfContents(indent=true, depth=4, aside=true)

# ╔═╡ 2c7cfabe-68ca-11eb-1bee-338af63e7e37
begin
	f(x) = 3 * x^2 + 6 * x - 9

	with_terminal() do
		for x = -5:5
  			println("(",x, ", ", f(x), ")")
		end
	end
end

# ╔═╡ e63cdbe2-68cb-11eb-20db-4dd54383d017
# pluto interaction

@bind vx html"<input type='range' min='1' max='50' step='0.5' />"

# ╔═╡ 4807bf12-68ca-11eb-37e6-81504bbf1eae
begin
	gr()                  ## Activate the GR backend for use with Plots
	plot(f, -vx, vx)      ## plot f over [-4,3]

	plot!(zero, -vx, vx)  ## horizontal line at y = 0
end

# ╔═╡ 7109711c-68ca-11eb-2c1e-5d0feb06d4cd
with_terminal() do
	println("f(-3), is a root: $(f(-3))")
	println("f(1), is the other root: $(f(1))")
end

# ╔═╡ a2dc225c-68ca-11eb-14b9-3be05028c029
with_terminal() do
	println("f(-1), is the min: $(f(-1))")
end

# ╔═╡ b652dc2c-68ca-11eb-1388-4bff4ac93f75
data = [1.6800483  -1.641695388;
        0.501309281 -0.977697538;
        1.528012113 0.52771122;
        1.70012253 1.711524991;
        1.992493625 1.891000015;
        2.706075824 -0.463427794;
        2.994931927 -0.443566619;
        3.491852811 -1.275179133;
        3.501191722 -0.690499597;
        4.459924502 -5.516130799;
        4.936965851 -6.001703074;
        5.023289852 -8.36416901;
        5.04233698 -7.924477517;
        5.50739285 -10.77482371;
        5.568665171 -10.9171878]

# ╔═╡ c21d59ec-68ca-11eb-04e6-d3453387da46
begin
	gr()           # Activate the GR backend for use with Plots

	x, y = data[:, 1], data[:, 2]
	plot(x, y, linetype=:scatter, leg=false)
	# scatter(x, y) # this is an alternative method, but does make a legend
end

# ╔═╡ d8bb398a-68ca-11eb-2cb7-c799ea584ad5
begin
	gr()           # Activate the GR backend for use with Plots

	# x, y = data[:, 1], data[:, 2]
	scatter(x, y)  # this is an alternative method, but does make a legend
end

# ╔═╡ fd1df772-68ca-11eb-3d84-2306e2f8b2ff
begin
	n = size(data, 1)
	view(data, n-3:n, :)  # displays 4 last lines
	#                     # no copy made, contrary to data[end-3:end, :]
end

# ╔═╡ 5bb58dae-68cb-11eb-0f14-15a9b6a27a93
begin
	m = 20
	nx = sort(rand(m)); ny = rand(m)

	Plots.scatter(nx, ny)
end

# ╔═╡ 84c33e12-68cb-11eb-2c3f-91955e670357
plot!(nx, ny, leg=false, title="A sample plot")

# ╔═╡ 903feeaa-68cb-11eb-1544-c18b07df0652
plot(nx, ny, leg=false, title="Another sample plot")

# ╔═╡ 223e2023-c974-45ac-920a-489e20bfc3f3
html"""
<style>
  main {
        max-width: calc(800px + 25px + 6px);
  }
  .plutoui-toc.aside {
    background: linen;
  }
  h3, h4 {
        background: wheat;
  }
</style>
"""

# ╔═╡ Cell order:
# ╟─f61ae738-68c9-11eb-03c1-efb38291b3ae
# ╠═12333bc8-68ca-11eb-2e89-ab9b76b0595b
# ╠═27b37bc0-68ca-11eb-0c16-914e5a5ec540
# ╠═2c7cfabe-68ca-11eb-1bee-338af63e7e37
# ╠═e63cdbe2-68cb-11eb-20db-4dd54383d017
# ╠═4807bf12-68ca-11eb-37e6-81504bbf1eae
# ╠═7109711c-68ca-11eb-2c1e-5d0feb06d4cd
# ╠═a2dc225c-68ca-11eb-14b9-3be05028c029
# ╠═b652dc2c-68ca-11eb-1388-4bff4ac93f75
# ╠═c21d59ec-68ca-11eb-04e6-d3453387da46
# ╠═d8bb398a-68ca-11eb-2cb7-c799ea584ad5
# ╠═fd1df772-68ca-11eb-3d84-2306e2f8b2ff
# ╠═5bb58dae-68cb-11eb-0f14-15a9b6a27a93
# ╠═84c33e12-68cb-11eb-2c3f-91955e670357
# ╠═903feeaa-68cb-11eb-1544-c18b07df0652
# ╟─223e2023-c974-45ac-920a-489e20bfc3f3

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks


text:
import Base.convert
import Base.==

@compat abstract type InferenceType end
@compat abstract type PackedInferenceType end

@compat abstract type FunctorInferenceType <: Function end

abstract type InferenceVariable end

# been replaced by Functor types, but may be reused for non-numerical cases
@compat abstract type Pairwise <: InferenceType end
@compat abstract type Singleton <: InferenceType end

@compat abstract type FunctorSingleton <: FunctorInferenceType end
# @compat abstract type FunctorPartialSingleton <: FunctorInferenceType end
@compat abstract type FunctorSingletonNH <: FunctorSingleton end

@compat abstract type FunctorPairwise <: FunctorInferenceType end
@compat abstract type FunctorPairwiseMinimize <: FunctorInferenceType end
@compat abstract type FunctorPairwiseNH <: FunctorPairwise end
# @compat abstract type FunctorPairwiseNHMinimize <: FunctorPairwiseMinimize end # TODO

struct ContinuousScalar <: InferenceVariable
  dims::Int
  ContinuousScalar() = new(1)
end
struct ContinuousMultivariate <:InferenceVariable
  dims::Int
  ContinuousMultivariate() = new()
  ContinuousMultivariate(x) = new(x)
end

const FGG = Graphs.GenericIncidenceList{Graphs.ExVertex,Graphs.Edge{Graphs.ExVertex},Array{Graphs.ExVertex,1},Array{Array{Graphs.Edge{Graphs.ExVertex},1},1}}
const FGGdict = Graphs.GenericIncidenceList{Graphs.ExVertex,Graphs.Edge{Graphs.ExVertex},Dict{Int,Graphs.ExVertex},Dict{Int,Array{Graphs.Edge{Graphs.ExVertex},1}}}

# Condensed representation of KernelDensityEstimate, by saving points and bandwidth
mutable struct EasyMessage
  pts::Array{Float64,2}
  bws::Array{Float64,1}
  EasyMessage() = new()
  EasyMessage(a::Array{Float64,2}, b::Array{Float64,1}) = new(a,b)
  EasyMessage(p::BallTreeDensity) = new(getPoints(p), getBW(p)[:,1])
end

mutable struct FactorGraph
  g::FGGdict
  bn
  IDs::Dict{Symbol,Int}
  fIDs::Dict{Symbol,Int}
  id::Int
  nodeIDs::Array{Int,1} # TODO -- ordering seems improved to use adj permutation -- pending merge JuliaArchive/Graphs.jl/#225
  factorIDs::Array{Int,1}
  bnverts::Dict{Int,Graphs.ExVertex} # TODO -- not sure if this is still used, remove
  bnid::Int # TODO -- not sure if this is still used
  dimID::Int
  cg
  cgIDs::Dict{Int,Int} # cgIDs[exvid] = neoid
  sessionname::String
  registeredModuleFunctions::VoidUnion{Dict{Symbol, Function}}
  reference::VoidUnion{Dict{Symbol, Tuple{Symbol, Vector{Float64}}}}
  stateless::Bool
  FactorGraph() = new()
  FactorGraph(
    x1,
    x2,
    x3,
    x4,
    x5,
    x6,
    x7,
    x8,
    x9,
    x10,
    x11,
    x12,
    x13,
    x14,
    x15
   ) = new(
    x1,
    x2,
    x3,
    x4,
    x5,
    x6,
    x7,
    x8,
    x9,
    x10,
    x11,
    x12,
    x13,
    x14,
    x15,
    false )
end

function emptyFactorGraph(;reference::VoidUnion{Dict{Symbol, Tuple{Symbol, Vector{Float64}}}}=nothing)
    fg = FactorGraph(Graphs.incdict(Graphs.ExVertex,is_directed=false),
                     Graphs.incdict(Graphs.ExVertex,is_directed=true),
                    #  Dict{Int,Graphs.ExVertex}(),
                    #  Dict{Int,Graphs.ExVertex}(),
                     Dict{Symbol,Int}(),
                     Dict{Symbol,Int}(),
                     0,
                     [],
                     [],
                     Dict{Int,Graphs.ExVertex}(),
                     0,
                     0,
                     nothing,
                     Dict{Int,Int}(),
                     "",
                     Dict{Symbol, Function}(:IncrementalInference=>IncrementalInference.getSample), # TODO likely to be removed
                     reference  ) #evalPotential
    return fg
end

mutable struct VariableNodeData
  initval::Array{Float64,2}
  initstdev::Array{Float64,2}
  val::Array{Float64,2}
  bw::Array{Float64,2}
  BayesNetOutVertIDs::Array{Int,1}
  dimIDs::Array{Int,1}
  dims::Int
  eliminated::Bool
  BayesNetVertID::Int
  separator::Array{Int,1}
  groundtruth::VoidUnion{ Dict{ Tuple{Symbol, Vector{Float64}} } } # not packed yet
  softtype
  initialized::Bool
  VariableNodeData() = new()
  function VariableNodeData(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11)
    warn("Deprecated use of VariableNodeData(11 param), use 13 parameters instead")
    new(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11, nothing, true) # TODO ensure this is initialized true is working for most cases
  end
  VariableNodeData(x1::Array{Float64,2},
                   x2::Array{Float64,2},
                   x3::Array{Float64,2},
                   x4::Array{Float64,2},
                   x5::Vector{Int},
                   x6::Vector{Int},
                   x7::Int,
                   x8::Bool,
                   x9::Int,
                   x10::Vector{Int},
                   x11::VoidUnion{ Dict{ Tuple{Symbol, Vector{Float64}} } },
                   x12,
                   x13::Bool) =
    new(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13)
end

mutable struct PackedVariableNodeData
  vecinitval::Array{Float64,1}
  diminitval::Int
  vecinitstdev::Array{Float64,1}
  diminitdev::Int
  vecval::Array{Float64,1}
  dimval::Int
  vecbw::Array{Float64,1}
  dimbw::Int
  BayesNetOutVertIDs::Array{Int,1}
  dimIDs::Array{Int,1}
  dims::Int
  eliminated::Bool
  BayesNetVertID::Int
  separator::Array{Int,1}
  # groundtruth::VoidUnion{ Dict{ Tuple{Symbol, Vector{Float64}} } }
  softtype::String
  initialized::Bool
  PackedVariableNodeData() = new()
  PackedVariableNodeData(x1::Vector{Float64},
                         x2::Int,
                         x3::Vector{Float64},
                         x4::Int,
                         x5::Vector{Float64},
                         x6::Int,
                         x7::Vector{Float64},
                         x8::Int,
                         x9::Vector{Int},
                         x10::Vector{Int},
                         x11::Int,
                         x12::Bool,
                         x13::Int,
                         x14::Vector{Int},
                         x15::String,
                         x16::Bool) = new(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16)
end

mutable struct GenericWrapParam{T} <: FunctorInferenceType
  usrfnc!::T
  params::Vector{Array{Float64,2}}
  varidx::Int
  particleidx::Int
  measurement::Tuple #Array{Float64,2}
  samplerfnc::Function # TODO -- remove, since no required. Direct multiple dispatch at solve
  specialzDim::Bool
  partial::Bool
  GenericWrapParam{T}() where {T} = new()
  GenericWrapParam{T}(fnc::T, t::Vector{Array{Float64,2}}) where {T} = new(fnc, t, 1,1, (zeros(0,1),) , +, false, false)
  GenericWrapParam{T}(fnc::T, t::Vector{Array{Float64,2}}, i::Int, j::Int) where {T} = new(fnc, t, i, j, (zeros(0,1),) , +, false, false)
  GenericWrapParam{T}(fnc::T, t::Vector{Array{Float64,2}}, i::Int, j::Int, meas::Tuple, smpl::Function) where {T} = new(fnc, t, i, j, meas, smpl, false, false)
  GenericWrapParam{T}(fnc::T, t::Vector{Array{Float64,2}}, i::Int, j::Int, meas::Tuple, smpl::Function, szd::Bool) where {T} = new(fnc, t, i, j, meas, smpl, szd, false)
  GenericWrapParam{T}(fnc::T, t::Vector{Array{Float64,2}}, i::Int, j::Int, meas::Tuple, smpl::Function, szd::Bool, partial::Bool) where {T} = new(fnc, t, i, j, meas, smpl, szd, partial)
end

mutable struct FastRootGenericWrapParam{T} <: Function
  p::Vector{Int}
  perturb::Vector{Float64}
  X::Array{Float64,2}
  Y::Vector{Float64}
  xDim::Int
  zDim::Int
  gwp::GenericWrapParam{T}
  FastRootGenericWrapParam{T}(xArr::Array{Float64,2}, zDim::Int, residfnc::GenericWrapParam{T}) where {T} =
      new(collect(1:size(xArr,1)), zeros(zDim), xArr, zeros(size(xArr,1)), size(xArr,1), zDim, residfnc)
end

mutable struct GenericFunctionNodeData{T, S}
  fncargvID::Array{Int,1}
  eliminated::Bool
  potentialused::Bool
  edgeIDs::Array{Int,1}
  frommodule::S #Union{Symbol, AbstractString}
  fnc::T
  GenericFunctionNodeData{T, S}() where {T, S} = new()
  GenericFunctionNodeData{T, S}(x1, x2, x3, x4, x5, x6) where {T, S} = new(x1, x2, x3, x4, x5, x6)
end

FunctionNodeData{T <: Union{InferenceType, FunctorInferenceType}} = GenericFunctionNodeData{T, Symbol}
FunctionNodeData() = GenericFunctionNodeData{T, Symbol}()
FunctionNodeData(x1, x2, x3, x4, x5, x6) = GenericFunctionNodeData{T, Symbol}(x1, x2, x3, x4, x5, x6)

# typealias PackedFunctionNodeData{T <: PackedInferenceType} GenericFunctionNodeData{T, AbstractString}
PackedFunctionNodeData{T <: PackedInferenceType} = GenericFunctionNodeData{T, AbstractString}
PackedFunctionNodeData() = GenericFunctionNodeData{T, AbstractString}()
PackedFunctionNodeData(x1, x2, x3, x4, x5, x6) = GenericFunctionNodeData{T, AbstractString}(x1, x2, x3, x4, x5, x6)

function convert(::Type{PackedVariableNodeData}, d::VariableNodeData)
  return PackedVariableNodeData(d.initval[:],size(d.initval,1),
                              d.initstdev[:],size(d.initstdev,1),
                              d.val[:],size(d.val,1),
                              d.bw[:], size(d.bw,1),
                              d.BayesNetOutVertIDs,
                              d.dimIDs, d.dims, d.eliminated,
                              d.BayesNetVertID, d.separator,
                              string(d.softtype), d.initialized)
end
function convert(::Type{VariableNodeData}, d::PackedVariableNodeData)

  r1 = d.diminitval
  c1 = r1 > 0 ? floor(Int,length(d.vecinitval)/r1) : 0
  M1 = reshape(d.vecinitval,r1,c1)

  r2 = d.diminitdev
  c2 = r2 > 0 ? floor(Int,length(d.vecinitstdev)/r2) : 0
  M2 = reshape(d.vecinitstdev,r2,c2)

  r3 = d.dimval
  c3 = r3 > 0 ? floor(Int,length(d.vecval)/r3) : 0
  M3 = reshape(d.vecval,r3,c3)

  r4 = d.dimbw
  c4 = r4 > 0 ? floor(Int,length(d.vecbw)/r4) : 0
  M4 = reshape(d.vecbw,r4,c4)

  # TODO -- allow out of module type allocation (future feature, not currently in use)
  st = IncrementalInference.ContinuousMultivariate # eval(parse(d.softtype))

  return VariableNodeData(M1,M2,M3,M4, d.BayesNetOutVertIDs,
    d.dimIDs, d.dims, d.eliminated, d.BayesNetVertID, d.separator,
    nothing, st, d.initialized )
end
function VNDencoder(P::Type{PackedVariableNodeData}, d::VariableNodeData)
  return convert(P, d) #PackedVariableNodeData
end
function VNDdecoder(T::Type{VariableNodeData}, d::PackedVariableNodeData)
  return convert(T, d) #VariableNodeData
end

function compare(a::VariableNodeData,b::VariableNodeData)
    TP = true
    TP = TP && a.initval == b.initval
    TP = TP && a.initstdev == b.initstdev
    TP = TP && a.val == b.val
    TP = TP && a.bw == b.bw
    TP = TP && a.BayesNetOutVertIDs == b.BayesNetOutVertIDs
    TP = TP && a.dimIDs == b.dimIDs
    TP = TP && a.dims == b.dims
    TP = TP && a.eliminated == b.eliminated
    TP = TP && a.BayesNetVertID == b.BayesNetVertID
    TP = TP && a.separator == b.separator
    return TP
end

function ==(a::VariableNodeData,b::VariableNodeData, nt::Symbol=:var)
  return IncrementalInference.compare(a,b)
end

# heavy use of multiple dispatch for converting between packed and original data types during DB usage
function convert{T <: InferenceType, P <: PackedInferenceType}(::Type{FunctionNodeData{T}}, d::PackedFunctionNodeData{P})
  return FunctionNodeData{T}(d.fncargvID, d.eliminated, d.potentialused, d.edgeIDs,
          Symbol(d.frommodule), convert(T, d.fnc))
end
function convert{P <: PackedInferenceType, T <: InferenceType}(::Type{PackedFunctionNodeData{P}}, d::FunctionNodeData{T})
  return PackedFunctionNodeData{P}(d.fncargvID, d.eliminated, d.potentialused, d.edgeIDs,
          string(d.frommodule), convert(P, d.fnc))
end

# Functor version -- TODO, abstraction can be improved here
function convert{T <: FunctorInferenceType, P <: PackedInferenceType}(::Type{FunctionNodeData{GenericWrapParam{T}}}, d::PackedFunctionNodeData{P})
  usrfnc = convert(T, d.fnc)
  gwpf = prepgenericwrapper(Graphs.ExVertex[], usrfnc, getSample)
  return FunctionNodeData{GenericWrapParam{T}}(d.fncargvID, d.eliminated, d.potentialused, d.edgeIDs,
          Symbol(d.frommodule), gwpf) #{T}
end
function convert{P <: PackedInferenceType, T <: FunctorInferenceType}(::Type{PackedFunctionNodeData{P}}, d::FunctionNodeData{T})
  return PackedFunctionNodeData{P}(d.fncargvID, d.eliminated, d.potentialused, d.edgeIDs,
          string(d.frommodule), convert(P, d.fnc.usrfnc!))
end

function FNDencode{T <: FunctorInferenceType, P <: PackedInferenceType}(::Type{PackedFunctionNodeData{P}}, d::FunctionNodeData{T})
  return convert(PackedFunctionNodeData{P}, d) #PackedFunctionNodeData{P}
end
function FNDdecode{T <: FunctorInferenceType, P <: PackedInferenceType}(::Type{FunctionNodeData{T}}, d::PackedFunctionNodeData{P})
  return convert(FunctionNodeData{T}, d) #FunctionNodeData{T}
end

function FNDencode{T <: InferenceType, P <: PackedInferenceType}(::Type{PackedFunctionNodeData{P}}, d::FunctionNodeData{T})
  return convert(PackedFunctionNodeData{P}, d) #PackedFunctionNodeData{P}
end
function FNDdecode{T <: InferenceType, P <: PackedInferenceType}(::Type{FunctionNodeData{T}}, d::PackedFunctionNodeData{P})
  return convert(FunctionNodeData{T}, d) #FunctionNodeData{T}
end

# Compare FunctionNodeData
function compare{T,S}(a::GenericFunctionNodeData{T,S},b::GenericFunctionNodeData{T,S})
  # TODO -- beef up this comparison to include the gwp
  TP = true
  TP = TP && a.fncargvID == b.fncargvID
  TP = TP && a.eliminated == b.eliminated
  TP = TP && a.potentialused == b.potentialused
  TP = TP && a.edgeIDs == b.edgeIDs
  TP = TP && a.frommodule == b.frommodule
  TP = TP && typeof(a.fnc) == typeof(b.fnc)
  return TP
end

function addGraphsVert!(fgl::FactorGraph,
            exvert::Graphs.ExVertex;
            labels::Vector{<:AbstractString}=String[])
  #
  Graphs.add_vertex!(fgl.g, exvert)
end

function getVertNode(fgl::FactorGraph, id::Int; nt::Symbol=:var, bigData::Bool=false)
  return fgl.g.vertices[id] # check equivalence between fgl.v/f[i] and fgl.g.vertices[i]
  # return nt == :var ? fgl.v[id] : fgl.f[id]
end
function getVertNode(fgl::FactorGraph, lbl::Symbol; nt::Symbol=:var, bigData::Bool=false)
  return getVertNode(fgl, (nt == :var ? fgl.IDs[lbl] : fgl.fIDs[lbl]), nt=nt, bigData=bigData)
end
getVertNode{T <: AbstractString}(fgl::FactorGraph, lbl::T; nt::Symbol=:var, bigData::Bool=false) = getVertNode(fgl, Symbol(lbl), nt=nt, bigData=bigData)

# excessive function, needs refactoring
function updateFullVertData!(fgl::FactorGraph,
    nv::Graphs.ExVertex;
    updateMAPest::Bool=false )
  #

  # not required, since we using reference -- placeholder function CloudGraphs interface
  # getVertNode(fgl, nv.index).attributes["data"] = nv.attributes["data"]
  nothing
end

function makeAddEdge!(fgl::FactorGraph, v1::Graphs.ExVertex, v2::Graphs.ExVertex; saveedgeID::Bool=true)
  edge = Graphs.make_edge(fgl.g, v1, v2)
  Graphs.add_edge!(fgl.g, edge)
  if saveedgeID push!(getData(v2).edgeIDs,edge.index) end #.attributes["data"]
  edge
end

function graphsOutNeighbors(fgl::FactorGraph, vert::Graphs.ExVertex; ready::Int=1,backendset::Int=1, needdata::Bool=false)
  Graphs.out_neighbors(vert, fgl.g)
end
function graphsOutNeighbors(fgl::FactorGraph, exVertId::Int; ready::Int=1,backendset::Int=1, needdata::Bool=false)
  graphsOutNeighbors(fgl.g, getVert(fgl,exVertId), ready=ready, backendset=backendset, needdata=needdata)
end

function graphsGetEdge(fgl::FactorGraph, id::Int)
  nothing
end

function graphsDeleteVertex!(fgl::FactorGraph, vert::Graphs.ExVertex)
  warn("graphsDeleteVertex! -- not deleting Graphs.jl vertex id=$(vert.index)")
  nothing
end

#


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

using Markdown
using InteractiveUtils

# ╔═╡ ecf99396-5e16-11ec-1735-c5307d7bf94a
data = [(rand(), rand(), rand()) for _ in 1:200]

# ╔═╡ 026cc2c9-059b-4d59-a95d-e0d884629b84
begin
	using Ripserer
	ripserer(data)
end

# ╔═╡ fca029d7-51f9-4ae2-96eb-d0aa5f086304
begin
	using Plots
	result=ripserer(data,dim_max=2)
	plot(plot(result),barcode(result))
end

# ╔═╡ Cell order:
# ╠═ecf99396-5e16-11ec-1735-c5307d7bf94a
# ╠═026cc2c9-059b-4d59-a95d-e0d884629b84
# ╠═fca029d7-51f9-4ae2-96eb-d0aa5f086304

max_stars_repo_name: BioMedAI/Pluto-Note

max_issues_repo_name: BioMedAI/Pluto-Note

max_forks_repo_name: BioMedAI/Pluto-Note


text:
function test_strictly_lexicographicallylessthan_vectorofvariables(
    model::MOI.ModelLike,
    config::MOIT.Config{T},
) where {T <: Real}
    @MOIT.requires MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer) # x1, x2, x3, x4
    @MOIT.requires MOI.supports_constraint(model, MOI.VariableIndex, MOI.EqualTo{T}) # c1, c3, c4
    @MOIT.requires MOI.supports_constraint(model, MOI.VariableIndex, MOI.GreaterThan{T}) # c2
    @MOIT.requires MOI.supports_constraint(model, MOI.VectorOfVariables, CP.Strictly{CP.LexicographicallyLessThan, T}) # c5

    x1, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x2, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x3, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x4, _ = MOI.add_constrained_variable(model, MOI.Integer())

    c1 = MOI.add_constraint(model, x1, MOI.EqualTo(2))
    c2 = MOI.add_constraint(model, x2, MOI.GreaterThan(1))
    c3 = MOI.add_constraint(model, x3, MOI.EqualTo(2))
    c4 = MOI.add_constraint(model, x4, MOI.EqualTo(2))

    c5 = MOI.add_constraint(
        model,
        MOI.VectorOfVariables([x1, x2, x3, x4]),
        CP.Strictly(CP.LexicographicallyLessThan(2)),
    )

    @test MOI.is_valid(model, x1)
    @test MOI.is_valid(model, x2)
    @test MOI.is_valid(model, x3)
    @test MOI.is_valid(model, x4)
    @test MOI.is_valid(model, c1)
    @test MOI.is_valid(model, c2)
    @test MOI.is_valid(model, c3)
    @test MOI.is_valid(model, c4)
    @test MOI.is_valid(model, c5)

    if MOIT._supports(config, MOI.optimize!)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.VariablePrimal(), x1) == 2
        @test MOI.get(model, MOI.VariablePrimal(), x2) in [1, 2]
        @test MOI.get(model, MOI.VariablePrimal(), x3) == 2
        @test MOI.get(model, MOI.VariablePrimal(), x4) == 2
    end
end

function MOIT.setup_test(
    ::typeof(test_strictly_lexicographicallylessthan_vectorofvariables),
    mock::MOIU.MockOptimizer,
    ::Config{T},
) where {T <: Real}
    MOIU.set_mock_optimize!(
        mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, T[2, 2, 2, 2]),
    )
    return
end

function test_strictly_lexicographicallylessthan_vectoraffinefunction(
    model::MOI.ModelLike,
    config::MOIT.Config{T},
) where {T <: Real}
    @MOIT.requires MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer) # x1, x2, x3, x4
    @MOIT.requires MOI.supports_constraint(model, MOI.VariableIndex, MOI.EqualTo{T}) # c1, c3, c4
    @MOIT.requires MOI.supports_constraint(model, MOI.VariableIndex, MOI.GreaterThan{T}) # c2
    @MOIT.requires MOI.supports_constraint(model, MOI.VectorOfVariables, CP.Strictly{CP.LexicographicallyLessThan, T}) # c5

    x1, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x2, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x3, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x4, _ = MOI.add_constrained_variable(model, MOI.Integer())

    c1 = MOI.add_constraint(model, one(T) * x1, MOI.EqualTo(2))
    c2 = MOI.add_constraint(model, one(T) * x2, MOI.GreaterThan(1))
    c3 = MOI.add_constraint(model, one(T) * x3, MOI.EqualTo(2))
    c4 = MOI.add_constraint(model, one(T) * x4, MOI.EqualTo(2))

    c5 = MOI.add_constraint(
        model,
        MOIU.vectorize(one(T) .* [x1, x2, x3, x4]),
        CP.Strictly(CP.LexicographicallyLessThan(2)),
    )

    @test MOI.is_valid(model, x1)
    @test MOI.is_valid(model, x2)
    @test MOI.is_valid(model, x3)
    @test MOI.is_valid(model, x4)
    @test MOI.is_valid(model, c1)
    @test MOI.is_valid(model, c2)
    @test MOI.is_valid(model, c3)
    @test MOI.is_valid(model, c4)
    @test MOI.is_valid(model, c5)

    if MOIT._supports(config, MOI.optimize!)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.VariablePrimal(), x1) == 2
        @test MOI.get(model, MOI.VariablePrimal(), x2) == 1
        @test MOI.get(model, MOI.VariablePrimal(), x3) == 2
        @test MOI.get(model, MOI.VariablePrimal(), x4) == 2
    end
end

function MOIT.setup_test(
    ::typeof(test_strictly_lexicographicallylessthan_vectoraffinefunction),
    mock::MOIU.MockOptimizer,
    ::Config{T},
) where {T <: Real}
    MOIU.set_mock_optimize!(
        mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, T[2, 1, 2, 2]),
    )
    return
end


text:
# To run this script, `cd` to the `./test/fixtures` directory and then, from the Julia terminal, `include("./runner.jl")`.

import JSON

# Small values:
x = linspace( 1e-200, 1e-308, 1007 )

y = Array( Any, length(x) )
for i in eachindex(x)
	y[i] = bits( x[i] )
end

data = Dict([
	("x", x),
	("expected", y)
])

outfile = open("bits_1e-200_1e-308.json", "w")
write( outfile, JSON.json(data) )
close( outfile );

# Medium values:
x = linspace( -1e3, 1e3, 1007 )

y = Array( Any, length(x) )
for i in eachindex(x)
	y[i] = bits( x[i] )
end

data = Dict([
	("x", x),
	("expected", y)
])

outfile = open("bits_-1e3_1e3.json", "w")
write( outfile, JSON.json(data) )
close( outfile );

# Large values:
x = linspace( 1e200, 1e308, 1007 )

y = Array( Any, length(x) )
for i in eachindex(x)
	y[i] = bits( x[i] )
end

data = Dict([
	("x", x),
	("expected", y)
])

outfile = open("bits_1e200_1e308.json", "w")
write( outfile, JSON.json(data) )
close( outfile );

# Subnormal values:
x = linspace( 1.12e-308, 5e-324, 1007 )

y = Array( Any, length(x) )
for i in eachindex(x)
	y[i] = bits( x[i] )
end

data = Dict([
	("x", x),
	("expected", y)
])

outfile = open("bits_1e-308_5e-324.json", "w")
write( outfile, JSON.json(data) )
close( outfile );

max_issues_repo_name: compute-io/bits

max_forks_repo_name: compute-io/bits


text:
using Lazy

"""
    Input a tree for display as a fringed structure.
"""
function fringe(tree)
    fringey(node::Pair) = [fringey(i) for i in node]
    fringey(leaf::Int) = leaf
    fringey(tree)
end

"""
    equalsfringe() uses a reduction to a lazy 1D list via
    getleaflist() for its "equality" of fringes
"""
getleaflist(tree::Int) = [tree]
getleaflist(tree::Pair) = vcat(getleaflist(seq(tree[1])), getleaflist(seq(tree[2])))
getleaflist(tree::Lazy.LazyList) = vcat(getleaflist(tree[1]), getleaflist(tree[2]))
getleaflist(tree::Void) = []
equalsfringe(t1, t2) = (getleaflist(t1) == getleaflist(t2))

a = 1 => 2 => 3 => 4 => 5 => 6 => 7 => 8
b = 1 => (( 2 => 3 ) => (4 => (5 => ((6 => 7) => 8))))
c = (((1 => 2) => 3) => 4) => 5 => 6 => 7 => 8

x = 1 => 2 => 3 => 4 => 5 => 6 => 7 => 8 => 9
y = 0 => 2 => 3 => 4 => 5 => 6 => 7 => 8
z = 1 => 2 => (4 => 3) => 5 => 6 => 7 => 8

prettyprint(s) = println(replace("$s", r"\{Any,1\}|Any|Array\{T,1\}\swhere\sT|Array|", ""))
prettyprint(fringe(a))
prettyprint(fringe(b))
prettyprint(fringe(c))
prettyprint(fringe(x))
prettyprint(fringe(y))
prettyprint(fringe(z))

prettyprint(getleaflist(a))
prettyprint(getleaflist(b))
prettyprint(getleaflist(c))

println(equalsfringe(a, a))
println(equalsfringe(a, b))
println(equalsfringe(a, c))
println(equalsfringe(b, c))
println(equalsfringe(a, x) == false)
println(equalsfringe(a, y) == false)
println(equalsfringe(a, z) == false)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{Core.kwftype(typeof(MusicXML.Type)),NamedTuple{(:channel, :program, :volume, :pan, :id), Tuple{Int64, Int64, Int64, Int64, String}},Type{MidiInstrument}})
    Base.precompile(Tuple{Core.kwftype(typeof(MusicXML.Type)),NamedTuple{(:name, :midiinstrument, :id), Tuple{String, MidiInstrument, String}},Type{ScorePart}})
    Base.precompile(Tuple{Core.kwftype(typeof(MusicXML.Type)),NamedTuple{(:partlist, :parts), Tuple{PartList, Vector{Part}}},Type{ScorePartwise}})
    Base.precompile(Tuple{Core.kwftype(typeof(MusicXML.Type)),NamedTuple{(:time, :divisions, :clef, :staves, :key), Tuple{Time, Int64, Vector{Clef}, Int64, Key}},Type{Attributes}})
    Base.precompile(Tuple{typeof(println),Base.PipeEndpoint,Pitch})
    Base.precompile(Tuple{typeof(println),Pitch})
    Base.precompile(Tuple{typeof(readmusicxml),String})
    Base.precompile(Tuple{typeof(vcat),Measure,Measure})
    Base.precompile(Tuple{typeof(writemusicxml),String,Vararg{Any, N} where N})
end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

for variables X in R^{n, m} and Y in S^n:
maximize    tr(C*X) :
subject to  Y - X*X' in S^n_+
            Y_ij = P_ij for (i, j) in Omega
where Omega is a set of fixed indices and P is a random PSD matrix

the nonlinear constraint Y - X*X' in S^n_+ is equivalent to
the conic constraint (Y, 0.5, X) in MatrixEpiPerSquareCone(),
and also to the larger conic constraint [I X'; X Y] in S^{n + m}_+

simple case of nonlinfear matrix inequality example from Lectures on Modern
Convex Optimization by Aharon Ben-Tal and Arkadi Nemirovski, pg 154
=#

using SparseArrays

struct MatrixQuadraticJuMP{T <: Real} <: ExampleInstanceJuMP{T}
    num_rows::Int
    num_cols::Int
    use_matrixepipersquare::Bool # use matrixepipersquare cone, else PSD cone
end

function build(inst::MatrixQuadraticJuMP{T}) where {T <: Float64}
    (num_rows, num_cols) = (inst.num_rows, inst.num_cols)
    C = randn(num_cols, num_rows)
    P = randn(num_rows, num_rows)
    P = Symmetric(P * P')
    (row_idxs, col_idxs, _) = findnz(tril!(sprand(Bool, num_rows,
        num_rows, inv(sqrt(num_rows)))) + I)

    model = JuMP.Model()
    JuMP.@variable(model, X[1:num_rows, 1:num_cols])
    JuMP.@variable(model, Y[1:num_rows, 1:num_rows], Symmetric)
    JuMP.@objective(model, Max, tr(C * X))
    JuMP.@constraint(model, [(row, col) in zip(row_idxs, col_idxs)],
        Y[row, col] == P[row, col])

    if inst.use_matrixepipersquare
        U_svec = zeros(JuMP.GenericAffExpr{T, JuMP.VariableRef},
            Cones.svec_length(num_rows))
        U_svec = Cones.smat_to_svec!(U_svec, 1.0 * Y, sqrt(2))
        JuMP.@constraint(model, vcat(U_svec, 0.5, vec(X)) in
            Hypatia.MatrixEpiPerSquareCone{T, T}(num_rows, num_cols))
    else
        JuMP.@constraint(model, Symmetric(
            [Matrix(I, num_cols, num_cols) X'; X Y]) in JuMP.PSDCone())
    end

    return model
end


text:
function no_slip_bc(grid)
   v_bcs = VVelocityBoundaryConditions(grid,
              top = BoundaryCondition(Value, 0.0),
           bottom = BoundaryCondition(Value, 0.0),
            north = BoundaryCondition(NormalFlow, 0.0),
            south = BoundaryCondition(NormalFlow, 0.0))
   w_bcs = WVelocityBoundaryConditions(grid,
           north = BoundaryCondition(Value, 0.0),
           south = BoundaryCondition(Value, 0.0),
             top = BoundaryCondition(NormalFlow, 0.0),
          bottom = BoundaryCondition(NormalFlow, 0.0))
   return v_bcs, w_bcs
end


text:
using Distributed
addprocs(16)
@everywhere include("plot_les.jl")
average_les()

max_stars_repo_name: nishaChandramoorthy/energies

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies


text:
module TriPMatrices

# package code goes here
include("tri_per.jl")
module Lapack
include("lapack.jl")
end
include("lu.jl")

end # module


text:
## covariance_functions.jl : utilities for Gaussian random field covariance functions

## CovarianceStructure ##
abstract type CovarianceStructure{T<:Real} end
abstract type IsotropicCovarianceStructure{T} <: CovarianceStructure{T} end
abstract type AnisotropicCovarianceStructure{T} <: CovarianceStructure{T} end

## CovarianceFunction ##
abstract type AbstractCovarianceFunction{d} end

# return number of dimension
Base.ndims(::AbstractCovarianceFunction{d}) where d = d

struct CovarianceFunction{d,C<:CovarianceStructure} <: AbstractCovarianceFunction{d}
    cov::C

    function CovarianceFunction{d,C}(cov::C) where {d,C}
        d > 0 || throw(DomainError(d, "dimension must be positive, got $(d)"))
        new{d,C}(cov)
    end
end

"""
	CovarianceFunction(d, cov)

Create a covariance function in `d` dimensions with covariance structure `cov`.

# Examples
```jldoctest
julia> m = Matern(0.1,1.0)
Matérn (λ=0.1, ν=1.0, σ=1.0, p=2.0)

julia> c = CovarianceFunction(2,m)
2d Matérn covariance function (λ=0.1, ν=1.0, σ=1.0, p=2.0)

```
"""
CovarianceFunction(d::Integer, cov::CovarianceStructure) =
    CovarianceFunction{d,typeof(cov)}(cov)

# return standard deviation of the Gaussian random field
Statistics.std(cov::CovarianceFunction) = cov.cov.σ

# evaluate the covariance function
apply(cov::CovarianceFunction, x, y) = apply(cov.cov, x, y)

# element type of covariance
Base.eltype(::CovarianceStructure{T}) where T = T
Base.eltype(::CovarianceFunction{d,<:CovarianceStructure{T}}) where {d,T} = T

# apply for isotropic random fields
apply(cov::IsotropicCovarianceStructure, dx::Vector{<:Real}) = apply(cov, norm(dx, cov.p))

# evaluate when pts is given as a kron product of 1d points
function apply(cov::CovarianceStructure, x::NTuple{d,AbstractVector},
               y::NTuple{d,AbstractVector}) where d
    C = zeros(eltype(cov), prod(length, x), prod(length, y))
    if  size(C, 1) == size(C, 2)
        return apply_symmetric!(C, cov, x, y)
    else
        return apply_non_symmetric!(C, cov, x, y)
    end
end

function apply_symmetric!(C::Matrix, cov::CovarianceStructure, x::NTuple{d,AbstractVector},
                          y::NTuple{d,AbstractVector}) where d
    z = Vector{eltype(C)}(undef, d)
    xiterator = enumerate(Iterators.product(x...))
    for (j, idy) in enumerate(Iterators.product(y...))
        for (i, idx) in Iterators.take(xiterator, j)
            @. z = idx - idy
	        @inbounds C[i, j] = apply(cov, z)
        end
    end
    Symmetric(C, :U)
end

function apply_non_symmetric!(C::Matrix, cov::CovarianceStructure,
                              x::NTuple{d,AbstractVector},
                              y::NTuple{d,AbstractVector}) where d
    z = Vector{eltype(C)}(undef, d)
    xiterator = enumerate(Iterators.product(x...))
    for (j, idy) in enumerate(Iterators.product(y...))
        for (i, idx) in xiterator
            @. z = idx - idy
            @inbounds C[i, j] = apply(cov, z)
        end
    end
    C
end

# evaluate when pts is given as a Finite Element mesh
function apply(cov::CovarianceStructure, tx::NTuple{2,AbstractMatrix},
               ty::NTuple{2,AbstractMatrix})
    x, y = first(tx), first(ty) # select FE nodes
    size(x) == size(y) || throw(DimensionMismatch())
    d, n = size(x)

    C = zeros(eltype(cov), n, n)
    z = Vector{eltype(C)}(undef, d)
    @inbounds for j in 1:n
        yj = view(y, :, j)
        for i in 1:j
            z .= view(x, :, i) .- yj
            C[i, j] = apply(cov, z)
        end
    end
    Symmetric(C, :U)
end

# evaluate for KL eigenfunctions
function apply(cov::CovarianceStructure, tx::NTuple{2,AbstractMatrix}, y::Tuple)
    x = first(tx) # select FE nodes
    d, nx = size(x)
    d == size(y, 1) || throw(DimensionMismatch())

    C = zeros(eltype(cov), nx, prod(length, y))
    z = Vector{eltype(C)}(undef, d)
    for (j, idy) in enumerate(Iterators.product(y...))
        for i in 1:nx
            z .= view(x, :, i) .- idy
            @inbounds C[i, j] = apply(cov, z)
        end
    end
    C
end

function Base.show(io::IO, cov::CovarianceFunction)
    print(io, ndims(cov), "d ", shortname(cov.cov), " covariance function")
    showparams(io, cov.cov)
end

function Base.show(io::IO, cov::CovarianceStructure)
    print(io, shortname(cov))
    showparams(io, cov)
end

function showparams(io::IO, cov::CovarianceStructure)
    n = nfields(cov)
    if n > 0
        names = fieldnames(typeof(cov))
        firstname = names[1]
        print(io, " (", firstname, "=", getfield(cov, firstname))
        for i in 2:n
            name = names[i]
            print(io, ", ", name, "=", getfield(cov, name))
        end
        print(io, ")")
    end
end


text:
# Returns all of the edges that are source or destination of a vertex
function edges_of(graph::AbstractGraph, vertex::Integer)
    Ev = Vector{AbstractEdge}()
    E = edges(graph)
    for e = E
        if vertex in (src(e),dst(e)); push!(Ev,e) end
    end
    return Ev
end

function unweighted_complement(graph::AbstractGraph)
    return complement(SimpleGraph(graph))
end


text:
macro exported_enum(name, args...)
    esc(quote
        @enum($name, $(args...))
        export $name
        $([:(export $arg) for arg in args]...)
    end)
end

@exported_enum SortBy relevance lastUpdatedDate submittedDate
@exported_enum SortOrder ascending descending
@exported_enum Field title author abstract comment jour_ref report_num id_list all_fields

mutable struct BibInfo
    url::String
    year::String
    authors::Vector{String}
    affiliations::Vector{String} # ! This field is currently not used in bib generation.
    title::String
    key::String
    journal::String
    BibInfo() = new()
end

n2f = Dict{Field,String}(
    title => "ti",
    author => "au",
    abstract => "abs",
    comment => "co",
    jour_ref => "jr",
    report_num => "rn",
    id_list => "id",
    all_fields => "all",
)


text:
@info "Installing MNE-Python"
using PyCall
pip = pyimport("pip")
flags = split(get(ENV, "PIPFLAGS", ""))
ver = get(ENV, "MNEVERSION", "")
packages = ["""mne$(isempty(ver) ? "" : "==")$(ver)"""]

@info "Package requirements:" packages
@info "Flags for pip install:" flags
ver = isempty(ver) ? "latest" : ver
@info "MNE version:" ver
pip.main(["install"; flags; packages])


text:
# Pose3Pose3 evaluation functions

"""
$(TYPEDEF)

Direct observation information of `Pose3` variable type.
"""
mutable struct PriorPose3 <: IncrementalInference.FunctorSingleton
    Zi::Distribution
    PriorPose3() = new()
    PriorPose3(st::FloatInt, sr::Float64) = new( MvNormal(zeros(6), [[st*Matrix{Float64}(LinearAlgebra.I, 3,3);zeros(3,3)];[zeros(3);sr*Matrix{Float64}(LinearAlgebra.I, 3,3)]] )  )
    PriorPose3(s::Distribution) = new(s)
end
function getSample(p3::PriorPose3, N::Int=1)
  # mv = Distributions.MvNormal(veeEuler(p3.Zi), p3.Cov)
  return (rand(p3.Zi, N),)
end

"""
$(TYPEDEF)

Serialization type for PriorPose3.
"""
mutable struct PackedPriorPose3  <: IncrementalInference.PackedInferenceType
    vecZi::Array{Float64,1} # 0rotations, 1translation in each column
    vecCov::Array{Float64,1}
    dimc::Int
    PackedPriorPose3() = new()
    PackedPriorPose3(x...) = new(x[1], x[2], x[3])
end
function convert(::Type{PriorPose3}, d::PackedPriorPose3)
  Zi = SE3(d.vecZi[1:3], Quaternion(d.vecZi[4],d.vecZi[5:7]))
  Cov = reshapeVec2Mat(d.vecCov, d.dimc)
  return PriorPose3( MvNormal( veeEuler(Zi), Cov) )
end
function convert(::Type{PackedPriorPose3}, d::PriorPose3)
  tf = SE3(d.Zi.μ[1:3], Euler(d.Zi.μ[4:6]...) )
  v1 = veeQuaternion(tf)
  v2 = d.Zi.Σ.mat[:];
  return PackedPriorPose3(v1, v2, size(d.Zi.Σ.mat,1))
end

# ------------------------------------

"""
$(TYPEDEF)
"""
mutable struct PP3REUSE
  wTi::SE3
  wTj::SE3
  iTi::SE3
  PP3REUSE() = new(SE3(0),SE3(0),SE3(0))
end

function fastpose3pose3residual!(reusethrid::PP3REUSE,
                                 res::Array{Float64},
                                 idx::Int,
                                 meas::Tuple,
                                 wXi::Array{Float64,2},
                                 wXj::Array{Float64,2}  )
  #
  reusethrid.wTi.t[1:3] = wXi[1:3,idx]
  TransformUtils.convert!(reusethrid.wTi.R, Euler(wXi[4,idx],wXi[5,idx],wXi[6,idx]))
  reusethrid.wTj.t[1:3] = wXj[1:3,idx]
  TransformUtils.convert!(reusethrid.wTj.R, Euler(wXj[4,idx],wXj[5,idx],wXj[6,idx]))

  # TODO -- convert to in place convert! functions, many speed-ups possible here
  jTi = SE3( matrix(reusethrid.wTj)\matrix(reusethrid.wTi) )
  # also wasted memory here, should operate directly on iTi and not be assigning new memory
  reusethrid.iTi = (SE3(meas[1][1:3,idx],Euler(meas[1][4:6,idx]...)) * jTi)
  res[:] = veeEuler(reusethrid.iTi)
  nothing
end

"""
$(TYPEDEF)

Rigid transform factor between two Pose3 compliant variables.
"""
mutable struct Pose3Pose3 <: FunctorPairwise
    Zij::Distribution
    reuse::Vector{PP3REUSE}
    Pose3Pose3() = new()
    Pose3Pose3(s) = new(s, PP3REUSE[PP3REUSE() for i in 1:Threads.nthreads()]  )
end
function getSample(pp3::Pose3Pose3, N::Int=1)
  return (rand(pp3.Zij, N), )
end
function (pp3::Pose3Pose3)(res::Array{Float64},
                           userdata::FactorMetadata,
                           idx::Int,
                           meas::Tuple,
                           wXi::Array{Float64,2},
                           wXj::Array{Float64,2}  )
  #
  reusethrid = pp3.reuse[Threads.threadid()]
  fastpose3pose3residual!(reusethrid, res, idx, meas, wXi, wXj)
  nothing
end

"""
$(TYPEDEF)

Serialization type for `Pose3Pose3`.
"""
mutable struct PackedPose3Pose3 <: IncrementalInference.PackedInferenceType
  vecZij::Array{Float64,1} # 3translations, 3rotation
  vecCov::Array{Float64,1}
  dimc::Int
  PackedPose3Pose3() = new()
  PackedPose3Pose3(x...) = new(x[1], x[2], x[3])
end
function convert(::Type{Pose3Pose3}, d::PackedPose3Pose3)
  qu = Quaternion(d.vecZij[4], d.vecZij[5:7])
  se3val = SE3(d.vecZij[1:3], qu)
  cov = reshapeVec2Mat(d.vecCov, d.dimc)
  return Pose3Pose3( MvNormal(veeEuler(se3val), cov) )
end
function convert(::Type{PackedPose3Pose3}, d::Pose3Pose3)
  val = d.Zij.μ
  se3val = SE3(val[1:3], Euler(val[4:6]...))
  v1 = veeQuaternion(se3val)
  v2 = d.Zij.Σ.mat
  return PackedPose3Pose3(v1[:], v2[:], size(v2,1) )
end

# -----------------------

"""
$(TYPEDEF)

Obsolete, see issue https://github.com/JuliaRobotics/IncrementalInference.jl/issues/237.
"""
mutable struct Pose3Pose3NH <: IncrementalInference.FunctorPairwiseNH
    Zij::Distribution
    nullhypothesis::Distributions.Categorical
    reuse::Vector{PP3REUSE}
    Pose3Pose3NH() = new()
    Pose3Pose3NH(s::Distribution, vh::Vector{Float64}) = new(s, Distributions.Categorical(vh), fill(PP3REUSE(), Threads.nthreads() )  )
    # Pose3Pose3NH(s::SE3, c::Array{Float64,2}, vh::Float64) = new(s,c, Distributions.Categorical([(1.0-vh);vh]),SE3(0),SE3(0),SE3(0))
    # Pose3Pose3NH(st::FloatInt, sr::Float64;vh::Float64=1.0) = new(SE3(0), [[st*Matrix{Float64}(LinearAlgebra.I, 3,3);zeros(3,3)];[zeros(3);sr*Matrix{Float64}(LinearAlgebra.I, 3,3)]], Distributions.Categorical([(1.0-vh);vh]),SE3(0),SE3(0),SE3(0))
end
function getSample(pp3::Pose3Pose3NH, N::Int=1)
  return (rand(pp3.Zij, N), )
end
function (pp3::Pose3Pose3NH)(res::Array{Float64},
            userdata,
            idx::Int,
            meas::Tuple,
            wXi::Array{Float64,2},
            wXj::Array{Float64,2}  )
  #
  reusethrid = pp3.reuse[Threads.threadid()]
  fastpose3pose3residual!(reusethrid, res, idx, meas, wXi, wXj)
  nothing
end

"""
$(TYPEDEF)

Obsolete, see issue https://github.com/JuliaRobotics/IncrementalInference.jl/issues/237.
"""
mutable struct PackedPose3Pose3NH <: IncrementalInference.PackedInferenceType
  vecZij::Vector{Float64} # 3translations, 3rotation
  vecCov::Vector{Float64}
  dimc::Int
  nullhypothesis::Vector{Float64}
  PackedPose3Pose3NH() = new()
  PackedPose3Pose3NH(x1::Vector{Float64},x2::Vector{Float64},x3::Int,x4::Vector{Float64}) = new(x1, x2, x3, x4)
end

function convert(::Type{Pose3Pose3NH}, d::PackedPose3Pose3NH)
  qu = Quaternion(d.vecZij[4], d.vecZij[5:7])
  se3val = SE3(d.vecZij[1:3], qu)
  cov = reshapeVec2Mat(d.vecCov, d.dimc)
  return Pose3Pose3NH( MvNormal(veeEuler(se3val), cov), d.nullhypothesis )
end
function convert(::Type{PackedPose3Pose3NH}, d::Pose3Pose3NH)
  val = d.Zij.μ
  se3val = SE3(val[1:3], Euler(val[4:6]...))
  v1 = veeQuaternion(se3val)
  v2 = d.Zij.Σ.mat
  return PackedPose3Pose3NH(v1[:], v2[:], size(v2,1), d.nullhypothesis.p )
end

# TODO -- stronger type safety required here
# Project all particles (columns) Xval with Z, that is for all  SE3(Xval[:,i])*Z
function projectParticles(Xval::Array{Float64,2}, Z::Distribution)
  # TODO optimize for more speed with threads and better memory management
  r,c = size(Xval)
  RES = zeros(r,c) #*cz

  ent, x = SE3(0), SE3(0)
  ENT = rand( Z, c )
  # ENT = rand( MvNormal(zeros(6), Cov), c )
  j=1
  # for j in 1:cz
  for i in 1:c
    x.R = TransformUtils.convert(SO3,Euler(Xval[4,i],Xval[5,i],Xval[6,i]))
    x.t = Xval[1:3,i]
    ent.R =  TransformUtils.convert(SO3, Euler(ENT[4:6,i]...)) # so3
    ent.t = ENT[1:3,i]
    # newval = Z*ent
    # res = x*newval
    res = x*ent
    RES[1:r,i*j] = veeEuler(res)
  end
  # end
  #
  return RES
end

⊕(Xpts::Array{Float64,2}, z::Pose3Pose3) = projectParticles(Xpts, z.Zij)
⊕(Xvert::Graphs.ExVertex, z::Pose3Pose3) = ⊕(getVal(Xvert), z)

#


text:
using Test

input = joinpath(@__DIR__, "input")
in = readlines(input)

adapters = map(x -> parse(Int, x), in)

function joltage(adapters)
    append!(adapters, 0)
    append!(adapters, maximum(adapters) + 3)
    sort!(adapters)

    deltas = diff(adapters)
    one = length(filter(x -> x == 1, deltas))
    three = length(filter(x -> x == 3, deltas))
    return one * three
end

function arrangements(adapters)
    append!(adapters, 0)
    append!(adapters, maximum(adapters) + 3)
    sort!(adapters)

    return arrangements(adapters, 1, Dict{Int, Int}())
end

function arrangements(adapters, idx, cache)
    haskey(cache, idx) && return cache[idx]
    idx == length(adapters) && return 1
    one = adapters[idx + 1] - adapters[idx] <= 3 ? arrangements(adapters, idx + 1, cache) : 0
    two = idx + 2 <= length(adapters) && adapters[idx + 2] - adapters[idx] <= 3 ? arrangements(adapters, idx + 2, cache) : 0
    three = idx + 3 <= length(adapters) && adapters[idx + 3] - adapters[idx] <= 3 ? arrangements(adapters, idx + 3, cache) : 0
    cache[idx] = one + two + three
    return one + two + three
end

sample = [
    16,
    10,
    15,
    5,
    1,
    11,
    7,
    19,
    6,
    12,
    4
]

@test joltage(copy(sample)) == 7 * 5
@test arrangements(copy(sample)) == 8

p1 = joltage(copy(adapters))
p2 = arrangements(copy(adapters))

@assert(p1 == 1984)
@assert(p2 == 3543369523456)

println("-----------------------------------------------------------------------")
println("adapter array -- part one :: $p1")
println("adapter array -- part two :: $p2")
println("-----------------------------------------------------------------------")