text:
module RulesProbitOutTest

using Test
using ReactiveMP
using Random
using Distributions
using StatsFuns: normcdf, normccdf, normlogcdf, normlogccdf, normlogpdf, normpdf

import ReactiveMP: @test_rules

@testset "rules:Probit:out" begin

    @testset "Belief Propagation: (m_in::UnivariateNormalDistribution, )" begin

        @test_rules [ with_float_conversions = true ] Probit(:out, Marginalisation) [
            (input = (m_in = NormalMeanVariance(1, 0.5), ), output = Bernoulli(normcdf(1/sqrt(1+0.5)))),
            (input = (m_in = NormalMeanPrecision(1, 2), ), output = Bernoulli(normcdf(1/sqrt(1+0.5)))),
            (input = (m_in = NormalWeightedMeanPrecision(2, 2), ), output = Bernoulli(normcdf(1/sqrt(1+0.5)))),
            (input = (m_in = NormalMeanVariance(2, 0.25), ), output = Bernoulli(normcdf(2/sqrt(1+0.25)))),
            (input = (m_in = NormalMeanPrecision(2, 4), ), output = Bernoulli(normcdf(2/sqrt(1+0.25)))),
            (input = (m_in = NormalWeightedMeanPrecision(8, 4), ), output = Bernoulli(normcdf(2/sqrt(1+0.25)))),
        ]

    end

    @testset "Belief Propagation: (m_in::PointMass, )" begin

        @test_rules [ with_float_conversions = true ] Probit(:out, Marginalisation) [
            (input = (m_in = PointMass(1), ), output = Bernoulli(normcdf(1))),
            (input = (m_in = PointMass(2), ), output = Bernoulli(normcdf(2))),
            (input = (m_in = PointMass(3), ), output = Bernoulli(normcdf(3)))
        ]

    end

end

end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    minimumcyclebasis

function findmincycle(graph, S)
    G = edgesubgraph(graph, setdiff(edgeset(graph), S))
    U = disjointunion(G, G)
    for s in S
        (u, v) = getedge(graph, s)
        u1 = getunionnode(U, 1, u)
        v1 = getunionnode(U, 1, v)
        u2 = getunionnode(U, 2, u)
        v2 = getunionnode(U, 2, v)
        addedge!(U, u1, v2, DisjointUnionEdge(1, s))
        addedge!(U, u2, v1, DisjointUnionEdge(1, s))
    end
    paths = []
    for n in nodeset(graph)
        n1 = getunionnode(U, 1, n)
        n2 = getunionnode(U, 2, n)
        sp = shortestpathedges(U, n1, n2)
        push!(paths, sp)
    end
    minpath = sortstablemin(paths, by=length)
    return [getsourceedge(U, e).sourcekey for e in minpath]
end

"""
    minimumcyclebasis(graph::UndirectedGraph) -> Vector{Vector{Int}}

Returns minimum cycle basis represented as an array of edge sequence that make up a cycle.

# Reference

1. de Pina, J.C. Applications of shortest path methods. PhD thesis, University of Amsterdam,
   Nether-lands (1995)
1. Kavitha, T., Mehlhorn, K., Michail, D. & Paluch, K. E. An ``\\tilde{O}(m^{2}n)`` Algorithm for
   Minimum Cycle Basis of Graphs. Algorithmica 52, 333–349 (2008).
"""
@cachefirst function minimumcyclebasis(graph::UndirectedGraph)
    cycles = Vector{Int}[]
    for conn in connectedcomponents(graph)
        subg = nodesubgraph(graph, conn)
        S = [Set(e) for e in cotree_edges(subg)]
        N = length(S)  # N: circuit rank
        for k in 1:N
            minedges = findmincycle(subg, S[k])
            push!(cycles, minedges)
            for i in (k + 1):N
                if length(intersect(S[i], minedges)) % 2 == 1
                    S[i] = symdiff(S[i], S[k])
                end
            end
        end
    end
    return cycles
end

ext: jl


text:
"""
    struct TensorProductQuadrature{D,T} <: Quadrature{D,T}
      coordinates::Vector{Point{D,T}}
      weights::Vector{T}
    end

Tensor product quadrature rule (nodes and weights) on a hyper cube [0,1]^D
"""
struct TensorProductQuadrature{D,T} <: Quadrature{D,T}
  coordinates::Vector{Point{D,T}}
  weights::Vector{T}
end

get_coordinates(q::TensorProductQuadrature) = q.coordinates

get_weights(q::TensorProductQuadrature) = q.weights

"""
    TensorProductQuadrature(degrees::NTuple{D}) where D
    TensorProductQuadrature(degrees::Point{D}) where D
"""
function TensorProductQuadrature(degrees::NTuple{D,Int}) where D
  TensorProductQuadrature{D}(degrees)
end

function TensorProductQuadrature(degrees::Point{D}) where D
  TensorProductQuadrature{D}(degrees)
end

"""
    TensorProductQuadrature{D}(degree::Integer) where D
    TensorProductQuadrature{D}(degrees) where D
"""
function TensorProductQuadrature{D}(degrees) where D
    @assert D == length(degrees)
    T = Float64
    npoints = [ ceil(Int,(degrees[i]+1.0)/2.0) for i in 1:D ]
    quads = [ gauss( eltype(Point{D,T}), npoints[i] ) for i in 1:D ]
    for i in 1:D
      quads[i][1] .+= 1;
      quads[i][1] .*= 1.0/2.0
      quads[i][2] .*= 1.0/2.0
    end
    (coords, weights) = _tensor_product(Point{D,T},quads,npoints)
    TensorProductQuadrature(coords,weights)
end

function TensorProductQuadrature{D}(degree::Integer) where D
  degrees = tfill(degree,Val{D}())
  TensorProductQuadrature(degrees)
end

# Helpers

function _tensor_product(::Type{Point{D,T}},quads,npoints) where {D,T}
  @assert length(quads) == D
  @assert length(npoints) == D
  cis = CartesianIndices(tuple(npoints...))
  n = prod(npoints)
  coords = Vector{Point{D,T}}(undef,n)
  weights = Vector{T}(undef,n)
  _tensor_product!(quads,coords,weights,cis)
  (coords, weights)
end

function _tensor_product(::Type{Point{0,T}},quads,npoints) where T
  coords = [zero(Point{0,T})]
  weights = [one(T)]
  coords, weights
end

function _tensor_product!(quads,coords,weights,cis)
  p = zero(mutable(eltype(coords)))
  D = length(p)
  lis = LinearIndices(cis)
  for ci in cis
    p[:] = 0.0
    w = 1.0
    for d in 1:D
      xi = quads[d][1][ci[d]]
      wi = quads[d][2][ci[d]]
      p[d] = xi
      w *= wi
    end
    li = lis[ci]
    coords[li] = p
    weights[li] = w
  end
end


text:
# vim:set ft=julia ts=4 sw=4 sts=4 autoindent:

# Hybrid from "Dynamic Programming Algorithms for Transition-Based Dependency
#   Parsers" by Kuhlmann et al. (2011).
#
# Author:   Pontus Stenetorp    <pontus stenetorp se>
# Version:  2014-05-06

# TODO: Is there a better way to signal the transition than using the types?
# TODO: Avoid creating undo objects if they are not needed? How?

module Hybrid

# TODO: Clean up the exports.
export FeatStruct, apply!, config, featurise, islegal, isterminal, oracle
export randoracle, transitions, undo!, NUMFEATS, featurise!

require("structs.jl")
require("conllx.jl")
require("dep.jl")

using Structs
using CoNLLX
using DepGraph

immutable Config
    stack::Vector{Vertex}
    buff::Vector{Vertex}
    graph::Graph
end

import Base: isequal
function isequal(a::Config, b::Config)
    (isequal(a.stack, b.stack) && isequal(a.buff, b.buff)
        && isequal(a.graph, b.graph))
end

import Base: show
function show(io::IO, c::Config)
    print(io, "Conf(",
        '[', join(map(v -> string(v.id), c.stack), ", "), ']',
        '|',
        '[', join(map(v -> string(v.id), reverse(c.buff)), ", "), ']',
        ')')
end

function config(sent::Sentence, coder)
    graph = dgraph(sent, coder, blind=true)
    return Config(graph[1:1], reverse(graph[2:end]), graph)
end

isterminal(c::Config) = isempty(c.buff) && length(c.stack) == 1

immutable Shift
    # Intentionally left blank.
end

function apply!(conf::Config, ::Shift)
    push!(conf.stack, pop!(conf.buff))
    return shiftundo
end

islegal(c::Config, ::Shift) = !isempty(c.buff)

import Base: isequal
isequal(::Shift, ::Shift) = true

import Base: hash
hash(::Shift) = 0x4711

immutable LeftArc
    deprel::String
end

function apply!(conf::Config, trans::LeftArc)
    head = conf.buff[end]
    dependent = pop!(conf.stack)

    undo = LeftArcUndo(dependent, trans.deprel, dependent.head,
        dependent.deprel)

    if trans.deprel == NOVAL
        edge!(conf.graph, dependent, head)
    else
        edge!(conf.graph, dependent, head, trans.deprel)
    end

    return undo
end

function islegal(c::Config, ::LeftArc)
    return !isempty(c.stack) && !isempty(c.buff) && c.stack[end].id != ROOTID
end

import Base: isequal
isequal(a::LeftArc, b::LeftArc) = a.deprel == b.deprel

import Base: hash
hash(a::LeftArc) = hash(a.deprel)

immutable RightArc
    deprel::String
end

function apply!(conf::Config, trans::RightArc)
    dependent = pop!(conf.stack)
    head = conf.stack[end]

    undo = RightArcUndo(dependent, trans.deprel, dependent.head,
        dependent.deprel)

    if trans.deprel == NOVAL
        edge!(conf.graph, dependent, head)
    else
        edge!(conf.graph, dependent, head, trans.deprel)
    end

    return undo
end

islegal(c::Config, ::RightArc) = length(c.stack) > 1

import Base: isequal
isequal(a::RightArc, b::RightArc) = a.deprel == b.deprel

import Base: hash
hash(a::RightArc) = hash(a.deprel)

immutable ShiftUndo
    # Intentionally left blank.
end

const shiftundo = ShiftUndo()

immutable LeftArcUndo
    dependent::Vertex
    deprel::String
    oldhead::Uint
    olddeprel::String
end

immutable RightArcUndo
    dependent::Vertex
    deprel::String
    oldhead::Uint
    olddeprel::String
end

undo!(conf::Config, ::ShiftUndo) = push!(conf.buff, pop!(conf.stack))

function undo!(conf::Config, undo::LeftArcUndo)
    head = conf.buff[end]
    dependent = undo.dependent

    deledge!(conf.graph, dependent, head)

    # Preserve any previous gold annotations.
    dependent.head = undo.oldhead
    dependent.deprel = undo.olddeprel

    push!(conf.stack, dependent)
end

function undo!(conf::Config, undo::RightArcUndo)
    head = conf.stack[end]
    dependent = undo.dependent

    deledge!(conf.graph, dependent, head)

    # Preserve any previous gold annotations.
    dependent.head = undo.oldhead
    dependent.deprel = undo.olddeprel

    push!(conf.stack, dependent)
end

# TODO: Variant with all deprels!
function oracle(conf::Config)
    if !isempty(conf.stack) && !isempty(conf.buff)
        s0 = conf.stack[end]
        b0 = conf.buff[end]
        if s0.head == b0.id
            return LeftArc(NOVAL)
            #return LeftArc(s0.deprel)
        end
    end

    if length(conf.stack) < 2
        return Shift()
    end

    s0 = conf.stack[end]
    s1 = conf.stack[end-1]
    # Require that a head is attached and that ther are no dependents
    #   remaining in the buffer.
    if s0.head == s1.id && !any(v -> v.head == s0.id, conf.buff)
        return RightArc(NOVAL)
        #return RightArc(s0.deprel)
    end

    return Shift()
end

# TODO: Variant with all deprels!
function transitions()
    return (Shift(), LeftArc(NOVAL), RightArc(NOVAL))
end

# Return a random legal transition (useful for debugging).
function randoracle(c::Config)
    legal = filter(t -> islegal(c, t), collect(transitions()))
    return legal[rand(1:end)]
end

immutable FeatStruct
    id::Uint
    arg1::Uint
    arg2::Uint
    arg3::Uint
    arg4::Uint
    # XXX: Why is this slower than `fstruct`?
    #function FeatStruct(id, arg1, arg2=0, arg3=0, arg4=0)
    #    FeatStruct(id, arg1, arg2, arg3, arg4)
    #end
end

function fstruct(id, arg1, arg2=0, arg3=0, arg4=0)
    FeatStruct(id, arg1, arg2, arg3, arg4)
end

import Base: hash
function hash(a::FeatStruct)
    hash(a.id) $ hash(a.arg1) $ hash(a.arg2) $ hash(a.arg3) $ hash(a.arg4)
end

import Base: show
function show(io::IO, fs::FeatStruct)
    print(io, string("FeatStruct(id=$(fs.id), $(fs.arg1), $(fs.arg2), ",
        "$(fs.arg3), $(fs.arg4))"))
end

# Allows us to generate the highly repetitive and painful featurisation code.
#
# Generates a constant and a function:
#
#   * NUMFEATS
#   * featurise!(feats::Vector{feattype}, c::Config)
#
macro featurise(feattype, feattypef, initblk, featsblk)
    exret = Expr(:block)

    featurisef = esc(quote
        function featurise!(featids, fidbyfstruct, c::Config)
            # initblk
            # feature extraction code generated from featblk
        end
    end)

    # Nasty, the block holding the instructions of the featurise function.
    featurisefblk = featurisef.args[1].args[2].args[2]
    # Inject the initialisation block into the function.
    append!(featurisefblk.args, initblk.args)

    seenfeats = Set()
    featpos = uint(0)
    for ex in featsblk.args
        # Ignore embedded line information.
        if ex.head == :line
            continue
        end
        # We assume that the non-id feature arguments are given as tuples.
        if ex.head != :tuple
            warn("ignoring expression $ex")
            continue
        end

        @assert !in(ex, seenfeats) "feature template duplicate: $ex"
        push!(seenfeats, ex)

        featpos += 1
        # Call to the feature structure function.
        fstructex = Expr(:call, feattypef, :($featpos), ex.args...)
        # Call to convert the feature structure into an id.
        idex = Expr(:call, id, :fidbyfstruct, fstructex)
        # Assign the id to its position in the output vector.
        assignex = Expr(:(=), :(featids[$featpos]), idex)
        push!(featurisefblk.args, assignex)
    end

    # Return the features.
    push!(featurisefblk.args, :(return featids))
    # Inject the function declaration into the returned expression.
    push!(exret.args, featurisef)

    # Define a constant for the total number of features.
    push!(exret.args, esc(:(const NUMFEATS = $featpos)))

    return exret
end

@featurise(FeatStruct, fstruct,
    begin
        # Short-hands.
        stack = c.stack
        buff = c.buff

        s0 = !isempty(stack) ? stack[end] : NULLVERT
        s1 = length(stack) > 1 ? stack[end-1] : NULLVERT
        s2 = length(stack) > 2 ? stack[end-2] : NULLVERT

        b0 = !isempty(buff) ? buff[end] : NULLVERT
        b1 = length(buff) > 1 ? buff[end-1] : NULLVERT

        s0_lm = lmostdep(s0)
        s1_lm = lmostdep(s1)
        b0_lm = lmostdep(b0)

        s0_rm = rmostdep(s0)
        s1_rm = rmostdep(s1)
    end,
    begin
        # Feature set used in "Training Deterministic Parsers with
        #   Non-Deterministic Oracles" by Goldberg et al. (2013) for
        #   the Hybrid model.

        # First order.
        (s0.form, )
        (s0.postag, )
        (s0.form,   s0.postag)
        (s1.form, )
        (s1.postag, )
        (s1.form,   s1.postag)
        (b0.form, )
        (b0.postag, )
        (b0.form,   b0.postag)
        (b1.form, )
        (b1.postag, )
        (b1.form,   b1.postag)

        # Second order.
        (s0.form,   s1.form)
        (s0.postag, s1.postag)
        (s0.postag, b0.postag)

        (s0.form,   s0.postag,  s1.postag)
        (s0.postag, s1.form,    s1.postag)
        (s0.form,   s1.form,    s1.postag)
        (s0.form,   s0.postag,  s1.form)
        # Goldberg et al. repeats an earlier feature template here. This is
        #   most likely an error, we correct it in order to conform with
        #   the next template group.
        #(s0.form,   s0.postag,  s1.postag)
        (s1.form,   s1.postag,  b0.postag)

        (s0.form,   b0.form)
        # Goldberg et al. repeats a template here, we omit it.
        #(s0.postag, b0.postag)
        (s0.postag, b1.postag)

        (s0.postag, b0.form,    b0.postag)
        (s0.form,   b0.form,    b0.postag)
        (s0.form,   s0.postag,  b0.form)
        (b0.form,   b0.postag,  b1.postag)

        # Third order.
        (s0.postag, b0.postag,  b1.postag)
        (s1.postag, s0.postag,  b0.postag)
        (s0.form,   b0.postag,  b1.postag)
        (s1.postag, s0.form,    b0.postag)

        # Fourth order.
        (s1.postag, s1_rm.postag,    b0.postag)
        # Goldberg et al. used the same feature ID for the feature above
        #   and below, we correct this in the macro call.
        (s1.postag, s1_rm.postag,   s0.postag)
        (s1.postag, s1_lm.postag,   s0.postag)
        (s1.postag, s0.postag,      s0_lm.postag)
        (s1.postag, s0_rm.postag,   s0.postag)
        (s0.postag, s1_lm.postag,   s0.form)
        (s1.postag, s0.form,        s0_rm.postag)
        # Goldberg et al. used the same feature ID for the next feature as for
        #   the one two steps above, we correct this in the macro call.
        (s0.postag, s1_rm.postag,   s0.form)
        # Goldberg et al. used the same feature ID for the next feature as for
        #   the one two steps above, we correct this in the macro call.
        (s1.postag, s0.form,        s0_lm.postag)

        (s0.postag, s0_rm.postag,   b0.postag)
        # Goldberg et al. used the same feature ID for the feature above
        #   and below, we correct this in the macro call.
        (s0.postag, s0_rm.postag,   b1.postag)
        (s0.postag, s0_lm.postag,   b0.postag)
        (s0.postag, b0.postag,      b0_lm.postag)
        (s0.postag, s0_lm.postag,   b0.form)
        (s0.postag, b0.form,        b0_lm.postag)

        (s0.postag, b0.postag,  b0_lm.postag,   s0_rm.postag)
        # Golberg et al. has a typo in the ID for the following feature,
        #   however, the ID still remains unique and should have no classifier
        #   performance impact.
        (s0.postag, s1.postag,  s0_lm.postag,   s1_rm.postag)

        # Fifth order.
        (s2.postag, s1.postag,  s0.postag)
        # Goldberg et al. uses the following two duplicates that should have
        #   no impact on the classifier performance, we omit them.
        #(s1.postag, s0.postag,  b0.postag)
        #(s0.postag, b0.postag,  b1.postag)
    end
)
end

max_stars_repo_name: ninjin/allen

max_issues_repo_name: ninjin/allen

max_forks_repo_name: ninjin/allen


text:
const StringLike = Union{AbstractString, Symbol}

macro ignore(args...) end

StringArray(length::Integer) = Vector{Cuchar}(undef, length)

function extof(path::AbstractString)
   m = match(r"\.\w+$", path)
   if !isnothing(m)
      return m.match
   else
      return ""
   end
end

@export macro offsetof(type, field)
   quote
      let fields = fieldnames($type)
         local index = findfirst(f -> f == $(QuoteNode(field)), fields)
         fieldoffset($type, index)
      end
   end |> esc
end

abstract type Iterable end

function Base.iterate(self::Iterable, state = 1)
   if state > fieldcount(typeof(self))
      return nothing
   end

   return getfield(self, state), state + 1
end

max_stars_repo_name: jorge-brito/Glutils

max_issues_repo_name: jorge-brito/Glutils

max_forks_repo_name: jorge-brito/Glutils


text:
using MacroTools: postwalk

"""
    @summary function ... end

Create a function to be passed to [`summarize`](@ref).

# Example

```julia
julia> @summary function energies(meta, data)
    id = "test"
    Ntw = meta.n_target
    len = nrow(data)
    μ_shift, σ_shift = mean_and_se(data.shift)
    μ_proj, _, _, σ_proj_l, σ_proj_u = med_and_errs(ratio_of_means(data.hproj, data.vproj))
end;

julia> summarize(energies, dir) # dir can be a path or a `FileTree`
10×8 DataFrame
 Row │ id      Ntw    len    μ_shift    σ_shift    μ_proj    σ_proj_l    σ_proj_u
     │ String  Int64  Int64  Float64    Float64    Float64   Float64     Float64
─────┼──────────────────────────────────────────────────────────────────────────────
   1 │ test      100  60000  -24.6576   0.145638   -2.23911  140.029     139.326
   2 │ test     1000  60000   -4.59347  0.0395217  -4.5541     0.110722    0.105786
   3 │ test      200  60000   -8.77444  0.100985   -1.08539    9.95662    23.2936
   4 │ test      300  60000   -4.4406   0.102841   -4.64745    0.325664    0.299909
   5 │ test      400  60000   -4.41512  0.0753507  -4.45474    0.239063    0.229831
   6 │ test      500  60000   -4.51584  0.0588216  -4.63888    0.190949    0.181465
   7 │ test      600  60000   -4.44399  0.0519303  -4.68629    0.147418    0.136815
   8 │ test      700  60000   -4.49641  0.0513159  -4.56       0.156817    0.146055
   9 │ test      800  60000   -4.58289  0.0455694  -4.56392    0.133013    0.128233
  10 │ test      900  60000   -4.49216  0.0429634  -4.74169    0.125266    0.119338

julia> @summary(function(meta, data)
    id = "test"
    Ntw = meta.n_target
    len = nrow(data)
    μ_shift, σ_shift = mean_and_se(data.shift)
    μ_proj, _, _, σ_proj_l, σ_proj_u = med_and_errs(ratio_of_means(data.hproj, data.vproj))
end) |> summarize(dir; skip=10_000)
10×8 DataFrame
 Row │ id      Ntw    len    μ_shift    σ_shift    μ_proj    σ_proj_l    σ_proj_u
     │ String  Int64  Int64  Float64    Float64    Float64   Float64     Float64
─────┼──────────────────────────────────────────────────────────────────────────────
   1 │ test      100  50000  -24.6273   0.154834   -4.05181  124.113     115.389
   2 │ test     1000  50000   -4.58082  0.0429057  -4.52871    0.120729    0.11253
   3 │ test      200  50000   -8.83234  0.111493   -1.77266   12.5945     22.3408
   4 │ test      300  50000   -4.39811  0.124721   -4.72357    0.351143    0.326857
   5 │ test      400  50000   -4.41533  0.0802611  -4.46048    0.269546    0.250092
   6 │ test      500  50000   -4.54364  0.0652312  -4.63162    0.19611     0.190288
   7 │ test      600  50000   -4.47056  0.0578356  -4.68352    0.165806    0.15526
   8 │ test      700  50000   -4.49484  0.0543399  -4.61254    0.152938    0.151639
   9 │ test      800  50000   -4.59333  0.0472654  -4.55456    0.148458    0.142926
  10 │ test      900  50000   -4.52037  0.047915   -4.70404    0.131165    0.125946

```
"""
macro summary(expr)
    if @capture(expr, (function fname_(args__) body_ end))
        header = Expr(:call, fname, args...)
    elseif @capture(expr, (function(args__) body_ end))
        header = Expr(:tuple, args...)
    else
        @show expr
        error("`expr` must be a `function` expression")
    end
    # Validate body
    for ex in body.args
        if ex isa LineNumberNode
            continue
        elseif ex isa Expr && ex.head == :(=)
            ex.args[2] = esc(ex.args[2])
        else
            error("only assignment expressions allowed in function body")
        end
    end
    header.args .= esc.(header.args)

    # Collect all names that appear on the left.
    left = Symbol[]
    postwalk(body) do ex
        if @capture(ex, ((names__,) = value_) | (names__ = value_))
            append!(left, names)
        else
            ex
        end
    end
    names = filter(!=(:_), (left...,))
    values = Expr(:tuple, names...)
    new_body = quote
        $body
        namedtuple($names, $values)
    end
    return Expr(:function, header, new_body)
end

"""
    summarize(f, root; skip=0)

Summarize a run by applying the function `f` to every pair of files in the `meta` and `raw`.
`f` should take a `NamedTuple` as the fist argument and `DataFrame` as the second. The first
`skip` rows of the `DataFrames` are skipped.

See also [`@summary`](@ref), a convenience macro for creating `f`.
"""
function summarize(f, root::FileTree; skip=0)
    meta = FileTrees.load(root["meta"]; lazy=true) do f
        namedtuple(JLSO.load(path(f)))
    end
    raw = FileTrees.load(root["raw"]; lazy=true) do f
        RimuIO.load_df(path(f))[skip + 1:end,:]
    end
    map(values(meta), values(raw)) do m, r
        f(exec(m), exec(r))
    end |> DataFrame
end
summarize(f, root::String; kwargs...) = summarize(f, FileTree(root); kwargs...)
summarize(root; kwargs...) = f -> summarize(f, FileTree(root); kwargs...)

export memory_use
function memory_use(df)
    if hasproperty(df, :len_before)
        return maximum(max.(df.len, df.len_before))
    else
        return maximum(df.len)
    end
end

export safe_proj_e
function safe_proj_e(df)
    try
        return med_and_errs(ratio_of_means(df.hproj, df.vproj))
    catch
        return (NaN, NaN, NaN, NaN, NaN)
    end
end


text:
@testset "explorers" begin

    @testset "EpsilonGreedyExplorer" begin
        @testset "API" begin
            explorer = EpsilonGreedyExplorer(0.1; is_break_tie = true)
            Random.seed!(explorer, 123)

            values = [0, 1, 2, -1]
            tarprob = [0.025, 0.025, 0.925, 0.025]

            # https://github.com/JuliaLang/julia/issues/10391#issuecomment-488642687
            # @test isapprox(prob(explorer, values), tarprob)
            @test isapprox(probs(prob(explorer, values)), tarprob)

            actions = [explorer(values) for _ in 1:10000]
            action_counts = countmap(actions)

            @test all(
                isapprox.(
                    [action_counts[i] for i in 1:length(values)] ./ 10000,
                    tarprob;
                    atol = 0.005,
                ),
            )
        end

        @testset "linear" begin
            explorer = EpsilonGreedyExplorer(;
                ϵ_stable = 0.1,
                ϵ_init = 0.9,
                warmup_steps = 3,
                decay_steps = 8,
                kind = :linear,
                is_break_tie = true,
            )
            xs = [0, 1, 2, -1, 2]
            mask = [true, true, false, true, false]
            E = [0.9, 0.9, range(0.9; stop = 0.1, step = -0.1)..., 0.1, 0.1, 0.1, 0.1]

            for ϵ in E
                @test RLCore.get_ϵ(explorer) ≈ ϵ
                @test isapprox(
                    probs(prob(explorer, xs)),
                    [ϵ / 5, ϵ / 5, ϵ / 5 + (1 - ϵ) / 2, ϵ / 5, ϵ / 5 + (1 - ϵ) / 2],
                )
                explorer(xs)
            end

            explorer = EpsilonGreedyExplorer(;
                ϵ_stable = 0.1,
                ϵ_init = 0.9,
                warmup_steps = 3,
                decay_steps = 8,
                kind = :linear,
                is_break_tie = true,
            )
            for ϵ in E
                @test RLCore.get_ϵ(explorer) ≈ ϵ
                @test isapprox(
                    probs(prob(explorer, xs, mask)),
                    [ϵ / 3, (1 - ϵ) + ϵ / 3, 0.0, ϵ / 3, 0.0],
                )
                explorer(xs)
            end

            explorer = EpsilonGreedyExplorer(;
                ϵ_stable = 0.1,
                ϵ_init = 0.9,
                warmup_steps = 3,
                decay_steps = 8,
                kind = :linear,
                is_break_tie = true,
            )
            for i in 1:100
                @test mask[explorer(xs, mask)]
            end
        end

        @testset "exp" begin
            explorer = EpsilonGreedyExplorer(;
                ϵ_stable = 0.1,
                ϵ_init = 0.9,
                warmup_steps = 3,
                decay_steps = 8,
                kind = :exp,
                is_break_tie = true,
            )
            xs = [0, 1, 2, -1, 2]
            mask = [true, true, false, true, false]
            for i in 1:10
                explorer(xs)
            end
            ϵ = 0.1 + (0.9 - 0.1) * exp(-1)
            @test isapprox(
                probs(prob(explorer, xs)),
                [ϵ / 5, ϵ / 5, ϵ / 5 + (1 - ϵ) / 2, ϵ / 5, ϵ / 5 + (1 - ϵ) / 2],
            )

            for i in 1:100
                explorer(xs)
            end
            ϵ = 0.1
            @test isapprox(
                probs(prob(explorer, xs)),
                [ϵ / 5, ϵ / 5, ϵ / 5 + (1 - ϵ) / 2, ϵ / 5, ϵ / 5 + (1 - ϵ) / 2];
                atol = 1e-5,
            )
        end
    end

end


text:
using StaticArrays: MVector

"""
    Spline{T}

Represents a spline function.

---

    Spline(B::AbstractBSplineBasis, coefs::AbstractVector)

Construct a spline from a B-spline basis and a vector of B-spline coefficients.

# Examples

```jldoctest; filter = r"coefficients: \\[.*\\]"
julia> B = BSplineBasis(BSplineOrder(4), -1:0.2:1);

julia> coefs = rand(length(B));

julia> S = Spline(B, coefs)
13-element Spline{Float64}:
 basis: 13-element BSplineBasis of order 4, domain [-1.0, 1.0]
 order: 4
 knots: [-1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0]
 coefficients: [0.815921, 0.076499, 0.433472, 0.672844, 0.468371, 0.348423, 0.868621, 0.0831675, 0.369734, 0.401199, 0.990734, 0.565907, 0.984855]
```

---

    Spline{T = Float64}(undef, B::AbstractBSplineBasis)

Construct a spline with uninitialised vector of coefficients.

---

    (S::Spline)(x)

Evaluate spline at coordinate `x`.
"""
struct Spline{
        T,  # type of coefficient (e.g. Float64, ComplexF64)
        Basis <: AbstractBSplineBasis,
        CoefVector <: AbstractVector{T},
    }
    basis :: Basis
    coefs :: CoefVector

    function Spline(B::AbstractBSplineBasis, coefs::AbstractVector)
        length(coefs) == length(B) ||
            throw(ArgumentError("wrong number of coefficients"))
        Basis = typeof(B)
        T = eltype(coefs)
        CoefVector = typeof(coefs)
        k = order(B)
        @assert k >= 1
        new{T, Basis, CoefVector}(B, coefs)
    end
end

Broadcast.broadcastable(S::Spline) = Ref(S)

Base.copy(S::Spline) = Spline(basis(S), copy(coefficients(S)))

function Base.show(io::IO, S::Spline)
    println(io, length(S), "-element ", nameof(typeof(S)), '{', eltype(S), '}', ':')
    print(io, " basis: ")
    summary(io, basis(S))
    println(io, "\n order: ", order(S))
    let io = IOContext(io, :compact => true, :limit => true)
        println(io, " knots: ", knots(S))
        print(io, " coefficients: ", coefficients(S))
    end
    nothing
end

Base.:(==)(P::Spline, Q::Spline) =
    basis(P) == basis(Q) && coefficients(P) == coefficients(Q)

Base.isapprox(P::Spline, Q::Spline; kwargs...) =
    basis(P) == basis(Q) &&
    isapprox(coefficients(P), coefficients(Q); kwargs...)

function Spline{T}(init, B::AbstractBSplineBasis) where {T}
    coefs = Vector{T}(init, length(B))
    Spline(B, coefs)
end

Spline(init, B::AbstractBSplineBasis) = Spline{Float64}(init, B)

# TODO deprecate?
Spline(init, B::AbstractBSplineBasis, ::Type{T}) where {T} =
    Spline{T}(init, B)

parent_spline(S::Spline) = parent_spline(basis(S), S)
parent_spline(::BSplineBasis, S::Spline) = S

"""
    coefficients(S::Spline)

Get B-spline coefficients of the spline.
"""
coefficients(S::Spline) = S.coefs

"""
    length(S::Spline)

Returns the number of coefficients in the spline.

Note that this is equal to the number of basis functions, `length(basis(S))`.
"""
Base.length(S::Spline) = length(coefficients(S))

"""
    eltype(::Type{<:Spline})
    eltype(S::Spline)

Returns type of element returned when evaluating the [`Spline`](@ref).
"""
Base.eltype(::Type{<:Spline{T}}) where {T} = T

"""
    basis(S::Spline) -> AbstractBSplineBasis

Returns the associated B-spline basis.
"""
basis(S::Spline) = S.basis

knots(S::Spline) = knots(basis(S))
order(::Type{<:Spline{T,Basis}}) where {T,Basis} = order(Basis)
order(S::Spline) = order(typeof(S))

# TODO allow evaluating derivatives at point `x` (should be much cheaper than
# constructing a new Spline for the derivative)
(S::Spline)(x) = _evaluate(basis(S), S, x)

function _evaluate(::BSplineBasis, S::Spline, x)
    T = eltype(S)
    t = knots(S)
    n = knot_interval(t, x)
    n === nothing && return zero(T)  # x is outside of knot domain
    k = order(S)
    spline_kernel(coefficients(S), t, n, x, BSplineOrder(k))
end

# Fallback, if the basis is not a regular BSplineBasis
_evaluate(::AbstractBSplineBasis, S::Spline, x) = parent_spline(S)(x)

function spline_kernel(
        c::AbstractVector{T}, t, n, x, ::BSplineOrder{k},
    ) where {T,k}
    # Algorithm adapted from https://en.wikipedia.org/wiki/De_Boor's_algorithm
    if @generated
        ex = quote
            @nexprs $k j -> d_j = @inbounds c[j + n - $k]
        end
        for r = 2:k, j = k:-1:r
            d_j = Symbol(:d_, j)
            d_p = Symbol(:d_, j - 1)
            jk = j - k
            jr = j - r
            ex = quote
                $ex
                α = @inbounds (x - t[$jk + n]) / (t[$jr + n + 1] - t[$jk + n])
                $d_j = $T((1 - α) * $d_p + α * $d_j)
            end
        end
        d_k = Symbol(:d_, k)
        quote
            $ex
            return $d_k
        end
    else
        # Similar using MVector (a bit slower than @generated version).
        spline_kernel_alt(c, t, n, x, BSplineOrder(k))
    end
end

function spline_kernel_alt(
        c::AbstractVector{T}, t, n, x, ::BSplineOrder{k},
    ) where {T, k}
    d = MVector(ntuple(j -> @inbounds(c[j + n - k]), Val(k)))
    @inbounds for r = 2:k
        dprev = d[r - 1]
        for j = r:k
            α = (x - t[j + n - k]) / (t[j + n - r + 1] - t[j + n - k])
            dtmp = dprev
            dprev = d[j]
            d[j] = (1 - α) * dtmp + α * dprev
        end
    end
    @inbounds d[k]
end

"""
    *(op::Derivative, S::Spline) -> Spline

Returns `N`-th derivative of spline `S` as a new spline.

See also [`diff`](@ref).

# Examples

```jldoctest; filter = r"coefficients: \\[.*\\]"
julia> B = BSplineBasis(BSplineOrder(4), -1:0.2:1);

julia> S = Spline(B, rand(length(B)))
13-element Spline{Float64}:
 basis: 13-element BSplineBasis of order 4, domain [-1.0, 1.0]
 order: 4
 knots: [-1.0, -1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0, 1.0]
 coefficients: [0.461501, 0.619799, 0.654451, 0.667213, 0.334672, 0.618022, 0.967496, 0.900014, 0.611195, 0.469467, 0.221618, 0.80084, 0.269533]

julia> Derivative(0) * S === S
true

julia> Derivative(1) * S
12-element Spline{Float64}:
 basis: 12-element BSplineBasis of order 3, domain [-1.0, 1.0]
 order: 3
 knots: [-1.0, -1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0, 1.0]
 coefficients: [2.37448, 0.259885, 0.0638088, -1.6627, 1.41675, 1.74737, -0.33741, -1.44409, -0.708643, -1.23925, 4.34416, -7.9696]

julia> Derivative(2) * S
11-element Spline{Float64}:
 basis: 11-element BSplineBasis of order 2, domain [-1.0, 1.0]
 order: 2
 knots: [-1.0, -1.0, -0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.0]
 coefficients: [-21.146, -0.98038, -8.63255, 15.3972, 1.65313, -10.4239, -5.53341, 3.67724, -2.65301, 27.917, -123.138]
```
"""
Base.:*(op::Derivative, S::Spline) = _diff(basis(S), S, op)

"""
    diff(S::Spline, [op::Derivative = Derivative(1)]) -> Spline

Same as `op * S`.

Returns `N`-th derivative of spline `S` as a new spline.
"""
Base.diff(S::Spline, op = Derivative(1)) = op * S

_diff(::AbstractBSplineBasis, S, etc...) = diff(parent_spline(S), etc...)

_diff(::BSplineBasis, S::Spline, ::Derivative{0}) = S

function _diff(
        ::BSplineBasis, S::Spline, ::Derivative{Ndiff} = Derivative(1),
    ) where {Ndiff}
    Ndiff :: Integer
    @assert Ndiff >= 1

    u = coefficients(S)
    t = knots(S)
    k = order(S)

    if Ndiff >= k
        throw(ArgumentError(
            "cannot differentiate order $k spline $Ndiff times!"))
    end

    Base.require_one_based_indexing(u)
    du = similar(u)
    copy!(du, u)

    @inbounds for m = 1:Ndiff, i in Iterators.Reverse(eachindex(du))
        dt = t[i + k - m] - t[i]
        if iszero(dt) || i == 1
            # In this case, the B-spline that this coefficient is
            # multiplying is zero everywhere, so we can set this to zero.
            # From de Boor (2001, p. 117): "anything times zero is zero".
            du[i] = zero(eltype(du))
        else
            du[i] = (k - m) * (du[i] - du[i - 1]) / dt
        end
    end

    # Finally, create lower-order spline with the given coefficients.
    # Note that the spline has `2 * Ndiff` fewer knots, and `Ndiff` fewer
    # B-splines.
    N = length(u)
    Nt = length(t)
    t_new = view(t, (1 + Ndiff):(Nt - Ndiff))
    B = BSplineBasis(BSplineOrder(k - Ndiff), t_new; augment = Val(false))

    Spline(B, view(du, (1 + Ndiff):N))
end

# Zeroth derivative: return S itself.
Base.diff(S::Spline, ::Derivative{0}) = S

"""
    integral(S::Spline)

Returns an antiderivative of the given spline as a new spline.

The algorithm is described in de Boor 2001, p. 127.
"""
integral(S::Spline) = _integral(basis(S), S)

_integral(::AbstractBSplineBasis, S, etc...) = integral(parent_spline(S), etc...)

function _integral(::BSplineBasis, S::Spline)
    u = coefficients(S)
    t = knots(S)
    k = order(S)
    Base.require_one_based_indexing(u)

    Nt = length(t)
    N = length(u)

    # Note that the new spline has 2 more knots and 1 more B-spline.
    t_int = similar(t, Nt + 2)
    t_int[2:(end - 1)] .= t
    t_int[1] = t_int[2]
    t_int[end] = t_int[end - 1]

    β = similar(u, N + 1)
    β[1] = zero(eltype(β))

    @inbounds for i in eachindex(u)
        m = i + 1
        β[m] = zero(eltype(β))
        for j = 1:i
            β[m] += u[j] * (t[j + k] - t[j]) / k
        end
    end

    B = BSplineBasis(BSplineOrder(k + 1), t_int; augment = Val(false))
    Spline(B, β)
end

function knot_interval(t::AbstractVector, x)
    n = searchsortedlast(t, x)  # t[n] <= x < t[n + 1]
    n == 0 && return nothing    # x < t[1]

    Nt = length(t)

    if n == Nt  # i.e. if x >= t[end]
        t_last = t[n]
        x > t_last && return nothing
        # If x is exactly on the last knot, decrease the index as necessary.
        while t[n] == t_last
            n -= one(n)
        end
    end

    n
end

max_stars_repo_name: jipolanco/BasisSplines

max_issues_repo_name: jipolanco/BasisSplines

max_forks_repo_name: jipolanco/BasisSplines


text:
# Examples from README.md

@test length(@bigint(genprimes(10^20, 10^20+1000))) == 24
@test length(genprimes(bi"10^20", bi"10^20+1000")) == 24
@test primepi(10^14+10^10; alg = :tabsieve) == 3205251958942
@test primepi(10^14+10^10; alg = :dr) == 3205251958942
@test primepi(10^14+10^8; alg = :dr) == 3204944853481
@test primepi(10^14+10^8; alg = :tabsieve) == 3204944853481
@test countprimes(100,1000) == 143
@test countprimes(100,tuplet=3) == 8
@test countprimes(10,1000,tuplet=6) == 1
res = countprimes("10^19+10^9")
## FIXME: bitrot
##@test res == 234057667299198865
@test typeof(res) == Int128
@test @bigint(countprimes(10^50, 10^50+1000)) == 7
@test primelookup(10^17 + 10^14 + 10^10) == (14,(2626111798288135,100100000000000000,10000000000))
primorial(n) = prod(someprimes(n))
@test @bigint(primorial(100)) == 2305567963945518424753102147331756070


text:
using LinearAlgebraicRepresentation
LAR = LinearAlgebraicRepresentation
using ViewerGL,SparseArrays
GL = ViewerGL

function svgarrangement(filename)
	V, EV = LAR.svg2lar(filename)
	GL.VIEW([ GL.GLLines(V,EV) ])
	model = Lar.Model(V)
	Lar.addModelCells!(model, 1, EV)
	model = Lar.Arrangement.planar_arrangement(model)
	bicon_comps = Lar.Arrangement.biconnected_components(model.T[1])

	# visualization of component graphs
	EV = LAR.cop2lar(model.T[1])
	V, EVs = LAR.biconnectedComponent((model.G, EV::LAR.Cells)) # 2-connected components
	VV = [[k] for k = 1 : size(V,2)]
	meshes = [ GL.numbering(.1)((V, (VV,EVs[k])), GL.COLORS[k],1) for k=1:length(EVs) ]
	GL.VIEW( cat(meshes) );

	# final solid visualizations
	FE = [SparseArrays.findnz(model.T[2][k,:])[1] for k=1:size(model.T[2],1)]
	FV = [collect(Set(cat([EV[e] for e in FE[f]]))) for f=1:length(FE)]

	EVs = Lar.FV2EVs(model)
	FVs = Lar.triangulate2D(model)
	GL.VIEW(GL.GLExplode(model.G, FVs, 1.2, 1.2, 1.2, 99, 1));
	GL.VIEW(GL.GLExplode(model.G, FVs,   1,   1,   1, 99, 1));
	GL.VIEW(GL.GLExplode(model.G, EVs,   1,   1,   1, 99, 1));
end


text:
# Autogenerated wrapper script for boostOpenMPI_jll for x86_64-linux-gnu-cxx03
export libboost_mpi, libboost_serialization, libboost_wserialization

using boost_jll
using OpenMPI_jll
JLLWrappers.@generate_wrapper_header("boostOpenMPI")
JLLWrappers.@declare_library_product(libboost_mpi, "libboost_mpi.so.1.76.0")
JLLWrappers.@declare_library_product(libboost_serialization, "libboost_serialization.so.1.76.0")
JLLWrappers.@declare_library_product(libboost_wserialization, "libboost_wserialization.so.1.76.0")
function __init__()
    JLLWrappers.@generate_init_header(boost_jll, OpenMPI_jll)
    JLLWrappers.@init_library_product(
        libboost_mpi,
        "lib/libboost_mpi.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libboost_serialization,
        "lib/libboost_serialization.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libboost_wserialization,
        "lib/libboost_wserialization.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
a = 3
b = 2
c = 1

fs = 100
T = 1/fs
f1 = 5
f2 = 10
f3 = 30

t = T:T:fs

x1 = a*sin.(2*pi*f1*T*t)+b*sin.(2*pi*f2*T*t)+c*sin.(2*pi*f3*T*t);

@testset "Peak prominence" begin
    @testset "Reciprocity" begin
        pi, pp = peakprom(x1, Maxima())
        ni, np = peakprom(-x1, Minima())

        @test pi == ni
        @test pp == np
    end

    @testset "Prominence values" begin
        i, p = peakprom(sin.(1e-5:1e-5:9*pi), Maxima())

        @test p[[1,5]] ≈ [1., 1.] atol=1e-4
        @test p[[2,3,4]] ≈ [2., 2., 2.] atol=1e-4
    end

    @testset "Minimum prominence" begin
        minprom = 1.5
        i, p = peakprom(sin.(1e-5:1e-5:9*pi), Maxima())
        mi, mp = peakprom(sin.(1e-5:1e-5:9*pi), Maxima(), 1, minprom)

        @test begin
            for j in 1:5 # number of peaks in i and p
                if i ∈ mi
                    if p[j] < minprom
                        return false
                    end
                else
                    if p[j] >= minprom
                        return false
                    end
                end
            end
            return true
        end
    end
end


text:
abstract type AbstractTwoDimCA end
abstract type AbstractGameOfLife <: AbstractTwoDimCA end

struct GameOfLife{T<:Bool} <: AbstractGameOfLife
    generations::Int
    all_runs::AbstractArray{T}
end

"""
    function GameOfLife(initial_state::AbstractArray{T}, generations::Int)

Given a starting matrix return the evolution for given generations of Conway's Game of Life, as described in [1].

[1] Gardner, Martin. “Mathematical games: The fantastic combinations of John Conway’s new solitaire game “life”.” Scientific American 223.4 (1970): 120-123.
"""
function GameOfLife(initial_state::AbstractArray{T},
        generations::Int) where T<: Bool

    height, width = size(initial_state)

    all_runs = Array{T}(undef, height, width, generations)
    all_runs[:, :, 1] = initial_state

    for g = 2:generations
        for i = 1:height, j = 1:width
            live_neighbours = 0

            for ud = (i-1):(i+1), lr = (j-1):(j+1)
                #pbc
                if ud < 1; ud = height-ud; end
                if ud > height; ud = ud-height; end
                if lr < 1; lr = width-lr; end
                if lr > width; lr = lr-width; end

                if all_runs[ud, lr, g-1] == 1; live_neighbours+=1; end
            end

            all_runs[i, j, g] = apply_standard_rules(all_runs[i, j, g-1], live_neighbours)
        end
    end
    return GameOfLife{T}(generations, all_runs)
end

function apply_standard_rules(state::T,
        live_neighbours::Int) where T <: Bool

    if state == 1
        if live_neighbours == 3 || live_neighbours == 4
            return 1
        else
            return 0
        end
    elseif state == 0
        if live_neighbours == 3
            return 1
        else
            return 0
        end
    end
end


text:
# issue 462

using Gadfly, Dates

set_default_plot_size(6inch, 3inch)

a = [unix2datetime(100)]
b = [10]

plot(x=a, y=b, Geom.point)


text:
using BenchmarkTools
using CUDA
using Oceananigans
using Oceananigans.TurbulenceClosures
using Benchmarks

# Benchmark function

function benchmark_closure(Arch, Closure)
    grid = RegularRectilinearGrid(size=(128, 128, 128), extent=(1, 1, 1))
    model = IncompressibleModel(architecture=Arch(), grid=grid, closure=Closure())

    time_step!(model, 1) # warmup

    trial = @benchmark begin
        @sync_gpu time_step!($model, 1)
    end samples=10

    return trial
end

# Benchmark parameters

Architectures = has_cuda() ? [CPU, GPU] : [CPU]

Closures = [Nothing,
            IsotropicDiffusivity,
            AnisotropicDiffusivity,
            AnisotropicBiharmonicDiffusivity,
            TwoDimensionalLeith,
            SmagorinskyLilly,
            AnisotropicMinimumDissipation]

# Run and summarize benchmarks

print_system_info()
suite = run_benchmarks(benchmark_closure; Architectures, Closures)

df = benchmarks_dataframe(suite)
sort!(df, [:Architectures, :Closures], by=(string, string))
benchmarks_pretty_table(df, title="Turbulence closure benchmarks")

if GPU in Architectures
    df_Δ = gpu_speedups_suite(suite) |> speedups_dataframe
    sort!(df_Δ, :Closures, by=string)
    benchmarks_pretty_table(df_Δ, title="Turbulence closure CPU to GPU speedup")
end

for Arch in Architectures
    suite_arch = speedups_suite(suite[@tagged Arch], base_case=(Arch, Nothing))
    df_arch = speedups_dataframe(suite_arch, slowdown=true)
    sort!(df_arch, :Closures, by=string)
    benchmarks_pretty_table(df_arch, title="Turbulence closures relative performance ($Arch)")
end


text:
using Test, Statistics, LinearAlgebra, Parameters # Test deps

# import Util Types
import BlockOpt: @lencheck, DimensionError, @contract, ContractError, @restrict, AccessError

# BlockOpt Types
import BlockOpt:
    BlockOptBackend,
    Model,
    Driver,
    DriverOptions,
    Weaver,
    EvaluationCounter,
    EvaluationTimer,
    BlockOptProfile,
    BlockOptTrace,
    Simulation

# Model imports
import BlockOpt:
    name,
    obj,
    grad!,
    grad,
    objective,
    gradient,
    initial_iterate,
    formula,
    dimension,
    directory,
    objective!,
    gradient!,
    initial_iterate!,
    formula!,
    final,
    hess_sample,
    hessAD

# DriverOptions imports
import BlockOpt:
    samples,
    Δ_max,
    δ_tol,
    ϵ_tol,
    max_iterations,
    weave_level,
    log_level,
    samples!,
    Δ_max!,
    δ_tol!,
    ϵ_tol!,
    max_iterations!,
    weave_level!,
    log_level!,
    LogLevel,
    INFO,
    DEBUG,
    WARN,
    ERROR,
    WeaveLevel,
    NONE,
    ALL

# Driver imports
import BlockOpt:
    orth,
    pflag,
    QN_update,
    SR1,
    PSB,
    S_update,
    S_update_a,
    S_update_b,
    S_update_c,
    S_update_d,
    S_update_e,
    S_update_f

# Timer imports
import BlockOpt: Δt, on!, off!

# Counter imports
import BlockOpt: evaluations, increment!

# Profile imports
import BlockOpt: profile, trs_timer, trs_counter, ghs_timer, ghs_counter

# Weaver imports
import BlockOpt: weaver, weave!, f_vals, ∇f_norms, Δ_vals, p_norms, ρ_vals

# BlockOptTrace imports
import BlockOpt:
    model, driver, profile, log_level, directory, io, info!, debug!, warn!, error!

# Backend imports
import BlockOpt:
    fₖ,
    ∇fₖ_norm,
    pₖ_norm,
    Δₖ,
    ρ,
    model,
    driver,
    terminal,
    secantQN,
    blockQN,
    update_Sₖ,
    build_trs,
    solve_trs,
    build_trial,
    update_Δₖ,
    accept_trial,
    gAD,
    gHS,
    initialize

# Simulation imports
import BlockOpt: trace, backend, optimize!, optimize

# Recipes imports
import BlockOpt:
    rhotrace,
    rhotrace!,
    steptrace,
    steptrace!,
    radiustrace,
    radiustrace!,
    objtrace,
    objtrace!,
    gradtrace,
    gradtrace!

include("setup.jl")

@testset verbose = true "Utilities (See File: test/src/util.jl)  " begin

    include("src/util.jl")
end

@testset verbose = true "Model (See File: test/src/model.jl)     " begin

    include("src/model.jl")
end

@testset verbose = true "Options (See File: test/src/options.jl) " begin

    include("src/options.jl")
end

@testset verbose = true "Driver (See File: test/src/driver.jl)   " begin

    include("src/driver.jl")
end

@testset verbose = true "Library (See Directory: test/src/lib/)  " begin

    include("src/lib/trace.jl")

    include("src/lib/backend.jl")

    include("src/simulation.jl")

    include("src/lib/recipes.jl")
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct VirtualMachineExtensionImagesApi <: SwaggerApi
    client::Swagger.Client
end

"""

Gets a virtual machine extension image.
Param: location::String (required)
Param: publisherName::String (required)
Param: type::String (required)
Param: version::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: VirtualMachineExtensionImage
"""
function virtualMachineExtensionImagesGet(_api::VirtualMachineExtensionImagesApi, location::String, publisherName::String, type::String, version::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", VirtualMachineExtensionImage, "/subscriptions/{subscriptionId}/providers/Microsoft.Compute/locations/{location}/publishers/{publisherName}/artifacttypes/vmextension/types/{type}/versions/{version}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "location", location)  # type String
    Swagger.set_param(_ctx.path, "publisherName", publisherName)  # type String
    Swagger.set_param(_ctx.path, "type", type)  # type String
    Swagger.set_param(_ctx.path, "version", version)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets a list of virtual machine extension image types.
Param: location::String (required)
Param: publisherName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Vector{VirtualMachineExtensionImage}
"""
function virtualMachineExtensionImagesListTypes(_api::VirtualMachineExtensionImagesApi, location::String, publisherName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Vector{VirtualMachineExtensionImage}, "/subscriptions/{subscriptionId}/providers/Microsoft.Compute/locations/{location}/publishers/{publisherName}/artifacttypes/vmextension/types", ["azure_auth"])
    Swagger.set_param(_ctx.path, "location", location)  # type String
    Swagger.set_param(_ctx.path, "publisherName", publisherName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets a list of virtual machine extension image versions.
Param: location::String (required)
Param: publisherName::String (required)
Param: type::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: filter::String
Param: top::Int32
Param: orderby::String
Return: Vector{VirtualMachineExtensionImage}
"""
function virtualMachineExtensionImagesListVersions(_api::VirtualMachineExtensionImagesApi, location::String, publisherName::String, type::String, api_version::String, subscriptionId::String; filter=nothing, top=nothing, orderby=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Vector{VirtualMachineExtensionImage}, "/subscriptions/{subscriptionId}/providers/Microsoft.Compute/locations/{location}/publishers/{publisherName}/artifacttypes/vmextension/types/{type}/versions", ["azure_auth"])
    Swagger.set_param(_ctx.path, "location", location)  # type String
    Swagger.set_param(_ctx.path, "publisherName", publisherName)  # type String
    Swagger.set_param(_ctx.path, "type", type)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "\$filter", filter)  # type String
    Swagger.set_param(_ctx.query, "\$top", top)  # type Int32
    Swagger.set_param(_ctx.query, "\$orderby", orderby)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

export virtualMachineExtensionImagesGet, virtualMachineExtensionImagesListTypes, virtualMachineExtensionImagesListVersions


text:
using FastGroupBy, DataBench
using Base.Test
import DataFrames.DataFrame
using SplitApplyCombine
@time addprocs(4)
@time @everywhere using FastGroupBy
#using DataFrames, CSV
#import Base.ht_keyindex
# write your own tests here

@time b = run_juliadb_bench(1_000_000, 100)
@time b1 = run_juliadb_bench_pmap(1_000_000, 100)
# @time a = R_bench(1_000_000, 100; libpath = "C:/Users/dzj/Documents/R/win-library/3.4")
# c = Dict(n => b[n]/a[n][1] for n in names(a))
#
# @test length(c) == 11

@time b1 = run_juliadb_bench_pmap()
gc()
file1 = replace("test/results/julia_pmap $(now()).csv",":","")
using CSV
CSV.write(file1, DataFrame(;b1...));

@time b = run_juliadb_bench()
gc()
file1 = replace("test/results/julia $(now()).csv",":","")
using CSV
CSV.write(file1, DataFrame(;b...));

rd = readdir("test/results/")
substr(rd,1,1)

@time a = R_bench(;libpath = "C:/Users/dzj/Documents/R/win-library/3.4")
file1 = replace("test/results/r $(now()).csv",":","")
using CSV
CSV.write(file1, a);

# compute relativities to R's data.table
c = Dict(n => b[n]/a[n][1] for n in names(a))
@test length(c) == 11

# collate the results
abc = vcat(a, DataFrame(;b...), DataFrame(;c...))
abc[:role] = ["R","Julia","Ratio"];

# files = open("test/results/"*readdir("test/results/")[2], "r")
# abc1 = deserialize( files)
#
# using DataBench
# dt = createIndexedTable(Int64(2e9/8), 100)
# @time sumby(dt, :id6, :v1)
# gc()
# @time sumby(dt, :id6, :v1)
# gc()
# @time groupreduce(x->x[1], x->(x[2],1), (x,y) ->(x[1]+y[1],x[2]+y[2]), zip(column(dt,:id6), column(dt,:v1)))
# gc()
# @time groupreduce(x->x[1], x->(x[2],1), (x,y) ->(x[1]+y[1],x[2]+y[2]), zip(column(dt,:id6), column(dt,:v1)))
# gc()


text:
# round to nearest ties to even
round(x::Posit8) = ccall((:p8_roundToInt, SoftPositPath), Posit8, (Posit8,),x)
round(x::Posit16) = ccall((:p16_roundToInt, SoftPositPath), Posit16, (Posit16,),x)
round(x::Posit32) = ccall((:p32_roundToInt, SoftPositPath), Posit32, (Posit32,),x)

# for Posits with 1 exponent bit
for T in (:Posit8_1, :Posit16_1, :Posit24_1)
    @eval begin
        round(x::$T) =  ccall((:pX1_roundToInt, SoftPositPath), $T, ($T,),x)
    end
end

# for Posits with 1 exponent bit
for T in (:Posit8_2, :Posit16_2, :Posit24_2)
    @eval begin
        round(x::$T) =  ccall((:pX2_roundToInt, SoftPositPath), $T, ($T,),x)
    end
end

function round(x::T, r::RoundingMode{:Down}) where {T<:AbstractPosit}
    xr = round(x)
    if xr > x
        return xr-one(T)
    else
        return xr
    end
end

function round(x::T, r::RoundingMode{:Up}) where {T<:AbstractPosit}
    xr = round(x)
    if xr < x
        return xr+one(T)
    else
        return xr
    end
end


text:
#=
Created on Friday 4 December 2020
Last update: Saturday 26 December 2020

@author: Michiel Stock
michielfmstock@gmail.com

@author: Bram De Jaegher
bram.de.jaegher@gmail.com

Templates heavily based on the MIT course "Computational Thinking"

https://computationalthinking.mit.edu/Fall20/installation/
=#

module DSJulia
    using Markdown
    using Markdown: MD, Admonition
    using PlutoUI

    export check_answer, still_missing, keep_working, correct, not_defined, hint, fyi
    export Question, QuestionOptional, QuestionBlock, validate
    export NoDiff, Easy, Intermediate, Hard
    export ProgressTracker, grade
    export @safe
    export @terminal

    include("styles.jl")
    include("admonition.jl")
    include("question.jl")
    include("grading.jl")
    include("flatland.jl")

    # Solutions
    export Solutions
    module Solutions
        using Colors, Images
        include("solutions.jl")
        include("solutions_flatland.jl")
    end

    # Convenience macro for terminal printing in Pluto
    macro terminal(ex)
        return quote
            PlutoUI.with_terminal() do
                $(esc(ex))
            end
        end
    end
end


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

export txt2fst
"""
`txt2fst(path_or_text,isym[,osym];W=TropicalWeight{Float32})`

Loads a WFST from a text file using [OpenFST format](http://openfst.org/twiki/bin/view/FST/FstQuickTour).

`isym` and `osym` can be either dictionaries or paths to the input/output symbol list.

If `path_or_text` is not a path, the input string is parsed directly.
```julia
julia> txt2fst("0 1 a a 3.0
       1 2 b <eps> 77.0
       2 3 c <eps> 33.0
       3 4 d d 44.0
       4 90.0
       ",Dict("a"=>1,"b"=>2,"c"=>3,"d"=>4))
WFST #states: 5, #arcs: 4, #isym: 4, #osym: 4
|1/0.0f0|
a:a/3.0f0 → (2)
(2)
b:ϵ/77.0f0 → (3)
(3)
c:ϵ/33.0f0 → (4)
(4)
d:d/44.0f0 → (5)
((5/90.0f0))

```
"""
function txt2fst(path::AbstractString,
                 isym::AbstractDict,
                 osym::AbstractDict=isym;W=TropicalWeight{Float32})
  fst = WFST(isym,osym; W=W)
  fst_txt = isfile(path) ? read(path, String) : path
  for x in split(fst_txt,"\n"; keepempty=false)
    arc = split(x)
    s = parse(Int,arc[1])+1
    if length(arc) <= 2 # final node with no weight
      N =s-length(fst)
      if N > 0
        add_states!(fst,N)
      end
      w = length(arc) == 2 ? parse(W,arc[2]) : one(W)
      final!(fst,s,w)
    elseif length(arc) >= 4
      n = parse(Int,arc[2])+1
      i = typeof(arc[3]) <: AbstractString ? String(arc[3]) : parse(I,arc[3])
      o = typeof(arc[4]) <: AbstractString ? String(arc[4]) : parse(O,arc[4])
      w = length(arc) == 5 ? parse(W,arc[5]) : one(W)
      add_arc!(fst,s,n,i,o,w.x)
    end
  end
  initial!(fst,1)
  return fst
end

txt2fst(path::AbstractString,path_sym::AbstractString; kwargs...) =
txt2fst(path,txt2sym(path_sym); kwargs...)

txt2fst(path::AbstractString,path_isym::AbstractString,path_osym::AbstractString; kwargs...) =
txt2fst(path,txt2sym(path_isym),txt2sym(path_osym); kwargs...)

export txt2sym
"""
`txt2sym(path)`

Loads a symbol list from a text file into a `Dict` using [OpenFST convention](http://openfst.org/twiki/bin/view/FST/FstQuickTour).
"""
function txt2sym(path)
  sym_txt = read(path,String)
  sym = Dict{String,Int}()
  for x in split(sym_txt,"\n"; keepempty=false)
    z = strip.(split(x))
    if length(z) > 2
      throw(ErrorException("Invalid symbol table `$x` in `$path`"))
    end
    if z[1] != "<eps>"
      sym[String(z[1])] = parse(Int,z[2])
    end
  end
  return sym
end

export fst2dot
"""
`fst2dot(fst)`

Converts the `fst` to a string in the [dot format](https://graphviz.org/doc/info/lang.html).

Useful to visualise large FSTs using Graphviz.
"""
function fst2dot(fst::WFST)
  x = "digraph FST {\nrankdir=LR;\nbgcolor=\"transparent\"\n;\ncenter=1;\n"
  int2isym = get_iisym(fst)
  int2osym = get_iosym(fst)
  for (i,s) in enumerate(fst)
    if isfinal(fst,i)
      w = get_final_weight(fst,i)
      x *= "$i [label = \"$(i)/$(w)\", shape = circle, style = bold, fontsize = 14]\n"
    elseif isinitial(fst,i)
      w = get_initial_weight(fst,i)
      x *= "$i [label = \"$(i)/$(w)\", shape = doublecircle, style = solid, fontsize = 14]\n"
    else
      x *= "$i [label = \"$(i)\", shape = circle, style = solid, fontsize = 14]\n"
    end
    for arc in s
      ilabel = idx2string(get_ilabel(arc),int2isym)
      olabel = idx2string(get_olabel(arc),int2osym)
      n = get_nextstate(arc)
      w = get_weight(arc)
      x*="\t$(i) -> $(n) [label = \"$(ilabel):$(olabel)/$(w)\", fontsize = 14];\n"
    end
  end
  x *="}"
  return x
end

export fst2pdf
"""
`fst2pdf(fst,pdffile)`

Prints the `fst` to a pdf file using [Graphviz](https://graphviz.org/).

Requires Graphviz to be intalled.
"""
function fst2pdf(fst,file)
  tmp = mktempdir()
  dot = fst2dot(fst)
  dotfile = joinpath(tmp,"fst.dot")
  write(dotfile,dot)
  write(file, read(`dot -Tpdf $dotfile`))
end


text:
"""
Parse a wire definition into a list of line segments
wire_def should look like ["R234", "U56", "L78"] etc
"""
function parse_wire_def(wire_def)
    cursor = [0 0]
    directions = Dict(
        'U' => [ 0  1],
        'D' => [ 0 -1],
        'L' => [-1  0],
        'R' => [ 1  0]
    )

    line_segs = []
    for seg in wire_def
        direction = seg[1]
        new_cursor = cursor + directions[direction] * parse(Int32, seg[2:end])
        push!(line_segs, (cursor, new_cursor))
        cursor = new_cursor
    end
    return line_segs
end

"""
Test if two line segments intersect, and if so, return their intersection point
"""
function find_intersection(seg1, seg2)
    # We only need to support vertical vs horizontal segment intersections...
    # at least I think a vert-vert or horz-horz intersection doesn't count based on task definition!
    (sx1, sy1), (ex1, ey1) = seg1
    (sx2, sy2), (ex2, ey2) = seg2

    # We need to order our segments so start x / start y are lower than end x / end y in each case
    sx1, sy1, ex1, ey1 = min(sx1, ex1), min(sy1, ey1), max(sx1, ex1), max(sy1, ey1)
    sx2, sy2, ex2, ey2 = min(sx2, ex2), min(sy2, ey2), max(sx2, ex2), max(sy2, ey2)

    # Now let's arrange as horizontal first, then vertical, if possible
    if sx1 == ex1 && sy2 == ey2
        # Vert then horz. Swap them.
        sx1, sy1, ex1, ey1, sx2, sy2, ex2, ey2 = sx2, sy2, ex2, ey2, sx1, sy1, ex1, ey1
    elseif !(sy1 == ey1 && sx2 == ex2)
        # It's not horz then vert either: no intersections
        return nothing
    end

    # OK, now sy1 == ey1 and sx2 == ey2. Is there an intersection?
    return (sy1 in sy2:ey2 && sx2 in sx1:ex1) ? [sx2, sy1] : nothing
end

# Part 1: read and parse input (should be 2 definitions, one for each wire)
w = open(f -> [parse_wire_def(split(line, ',')) for line in readlines(f)], "day03_input.txt")

# Do a stupid all-segs-vs-all-segs test to find all intersection points
isects = []
for seg1 in w[1]
    for seg2 in w[2]
        isect = find_intersection(seg1, seg2)
        if !isnothing(isect)
            push!(isects, isect)
        end
    end
end

# Find lowest Manhattan distance to origin apart from the first (0,0) intersection
filter!(isect -> isect[1] != 0 || isect[2] != 0, isects)
dists = [abs(x) + abs(y) for (x, y) in isects]
println("Minimum distance: $(minimum(dists))")

# Part 2: distance to intersection
"""
Compute lowest distance to the given intersection for the wire defined by list of segments wire_segs
"""
function distance_to_intersection(wire_segs, intersection)
    dist = 0
    for wire_seg in wire_segs
        # Is this wire segment going to hit the intersection?
        (sx, sy), (ex, ey) = wire_seg
        ix, iy = intersection
        if ix == sx == ex && iy in min(sy,ey):max(sy,ey)
            # Vertical segment reaching the intersection
            return dist + abs(iy - sy)
        elseif iy == sy == ey && ix in min(sx,ex):max(sx,ex)
            # Horizontal segment reaching the intersection
            return dist + abs(ix - sx)
        else
            dist += abs(ex - sx) + abs(ey - sy)
        end
    end
    return nothing
end

# Compute distance to each intersection along wire
wdists = [distance_to_intersection(w[1], isect) + distance_to_intersection(w[2], isect) for isect in isects]
println("Lowest combined distance along wire: $(minimum(wdists))")

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
## Solve the PBCC objective for a range of parameter settings on a suite of small graphs
using MAT
using SparseArrays
include("../src/PBCC.jl")
graphs = ["Cities_H", "Newsgroups100", "Zoo_H", "Amazon_Appliances", "Amazon_Amazon_Fashion"]

##
C = [1 0 0;
0 .75 0;
0 0 1;
0 .5 .5;
.5 0 .5;
1 .5 0;
0 0 0;
.75 .75 0;
0.1 0.5 0.6;
0 .75 .75;
.4 .5 .3]

mu = 0.0
s1 = 300
s2 = 250
ms = 5
lw = 2
y_label = "Ratio to LP Lower Bound"
using LaTeXStrings
x_label = L"\mu"
##

using Plots
plot()
mus = collect(0.00:0.01:0.2)
beta = 0.5
Plotdata = zeros(5,length(mus))
for ii = 1:length(graphs)
    graph = graphs[ii]
    M = matread("../data/"*graph*".mat")
    H = M["H"]
    y1 = zeros(length(mus))
    bestDs = zeros(length(mus))
    for i = 1:length(mus)
        mu = mus[i]
        outputmat = "Output/"*graph*"/beta_$beta"*"_mu_$mu.mat"

        mat = matread(outputmat)
        X = mat["X"]
        LPbound = mat["LPbound"]

        C, Labels, obj, bestD = RangeRoundX(H,X,mu,beta,10)
        x,y = Biclustering(H,Labels)
        ratio = obj/LPbound
        y1[i] = ratio
        bestDs[i] = bestD

    end
    Plotdata[ii,:] = y1
    @show bestDs
    plot!(mus,y1,linewidth = lw, markershape =:circle, grid = false,
    color = RGBA(C[ii,1],C[ii,2],C[ii,3],1), label = graph, legend = :bottomright,
    size = (s1,s2), markerstrokewidth = 0, markersize = ms,
    xlabel = x_label, ylabel = y_label)
end

##
matwrite("Plots/Plotdata_for_smaller_mu.mat", Dict("Plotdata"=>Plotdata))
plot!(legend = false)
# ii = 1
# annotate!(.65, 1.5, text("Cities",font(9,RGBA(C[ii,1],C[ii,2],C[ii,3],1))))
# ii = 2
# annotate!(.2, 2.3, text("News100",font(9,RGBA(C[ii,1],C[ii,2],C[ii,3],1))))
# ii = 3
# annotate!(.53, 2.5, text("Zoo",font(9,RGBA(C[ii,1],C[ii,2],C[ii,3],1))))
# ii = 4
# annotate!(.85, 1.9, text("Appliances",font(9,RGBA(C[ii,1],C[ii,2],C[ii,3],1))))
# ii = 5
# annotate!(.8, 1.1, text("Fashion",font(9,RGBA(C[ii,1],C[ii,2],C[ii,3],1))))

#["Cities_H", "Newsgroups100", "Zoo_H", "Amazon_Appliances", "Amazon_Amazon_Fashion"]

savefig("Plots/Mu_smaller_bestround.pdf")

max_stars_repo_name: nveldt/ParamCC

max_issues_repo_name: nveldt/ParamCC

max_forks_repo_name: nveldt/ParamCC


text:
mutable struct Plot{T<:Sampling}
    cascade::Cascade{T}
    legend::Vector{Pair{Handle,Any}}
    axes::Vector
    title::Label{Vector{String}}
    path::String
end

function Plot( ; cascade, legend, xaxis, yaxis, title, path)
    return Plot(cascade, legend, [xaxis, yaxis], title, path)
end


text:
using ColorVectorSpace: AbstractGray, TransparentGray, TransparentRGB

# Specialize the low-level broadcasting machinery for ImageMeta.
# This make it work for all operators, rather than specializing
# each operator individually.
Base.BroadcastStyle(::Type{<:ImageMeta}) = Broadcast.ArrayStyle{ImageMeta}()
function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{ImageMeta}}, ::Type{ElType}) where ElType
    Mt = imagemeta(bc.args...)
    if length(Mt) > 0
        M = Mt[1]
        if length(Mt) > 1
            for i = 2:length(Mt)
                Mt[i] == M || return ambigop(:Broadcast)
            end
        end
        # Use the properties field of img to create the output
        ret = ImageMeta(similar(Array{ElType}, axes(bc)), properties(M))
    else
        ret = ImageMeta(similar(Array{ElType}, axes(bc)))
    end
    ret
end
# Select all the ImageMeta arrays
@inline imagemeta(As...) = _imagemeta((), As...)
_imagemeta(out) = out
@inline _imagemeta(out, A::ImageMeta, As...) = _imagemeta((out..., A), As...)
@inline _imagemeta(out, A, As...) = _imagemeta(out, As...)

(-)(img::ImageMeta) = shareproperties(img, -arraydata(img))

batch2 = (:+, :-)

for op in batch2
    @eval begin
        ($op)(img::ImageMeta, B::BitArray) = shareproperties(img, ($op)(arraydata(img), B))
        ($op)(B::BitArray, img::ImageMeta) = shareproperties(img, ($op)(B, arraydata(img)))
        ($op)(img::ImageMeta, B::ImageMeta) = ambigop(Symbol($op))
        ($op)(img::ImageMeta, B::AbstractArray) = shareproperties(img, ($op)(arraydata(img), B))
        ($op)(B::AbstractArray, img::ImageMeta) = shareproperties(img, ($op)(B, arraydata(img)))
    end

    for CV in (:AbstractGray, :TransparentGray, :AbstractRGB, :TransparentRGB)
        @eval begin
            ($op)(img::ImageMeta{CV}, n::$CV) where {CV<:$CV} =
                shareproperties(img, ($op)(arraydata(img), n))
            ($op)(n::$CV, img::ImageMeta{CV}) where {CV<:$CV} =
                shareproperties(img, ($op)(n, arraydata(img)))
        end
    end
end

(/)(img::ImageMeta, n::Number) = shareproperties(img, arraydata(img)/n)

ambigop(s::Symbol) = error("$s with two ImageMeta arrays: dictionary choice is ambiguous")


text:
export dipole_pocklington
# This file includes some testing functions

# ----------------------------------- Univariate test functions ----------------------------------- #
function dipole_pocklington(zp, z, ra=0.005, k=2π)
    R = sqrt.(ra^2 .+ (z .- zp).^2)
    exp.(-im*k .*R ) .*  ((1 .+ im*k .*R) .* (2 .* R.^2 .- 3*ra^2) .+ (k*ra .*R).^2) ./ (4π .*R.^5)
end

parabol(x) = 1 - (2x - 1)^2

sinusoid(x, f=1) = sin(2π * f * x)

weierstrass(x, a=0.5, b=3, N=500) = sum([a^n * cos(b^n *  π * x) for n = 1 : N])

cellerier(x, a=2, N=500) = sum([1 / (a^k) * sin(a^k * x) for k = 1 : N])

riemann(x, N=500) = sum([1 / k^2 * sin(k^2 * x) for k = 1 : N])

darboux(x, N=10) = sum([1 / factorial(k) * sin(factorial(k + 1) * x) for k = 1 : N])

function p(x)
    xi = abs(x) % 2
    0       ≤   xi  ≤ 1 / 3 ?   0           :
    1 / 3   ≤   xi  ≤ 2 / 3 ?   3 * xi - 1  :
    2 / 3   ≤   xi  ≤ 4 / 3 ?   1           :
    4 / 3   ≤   xi  ≤ 5 / 3 ?   5 - 3 * xi  :
    5 / 3   ≤   xi  ≤ 2     ?   0           :
    throw(DomainError("$xi is not in the expected domian"))
end
schoenberg1(x, N=500) = sum([(1 / 2)^k * p(3^(2k) * x) for k = 0 : N])
schoenberg2(x, N=500) = sum([(1 / 2)^k * p(3^(2k + 1) * x) for k = 0 : N])

function g(x)
    xi = abs(x) % 4
    0   ≤   xi  ≤ 2 ?  1 - xi :
    2   ≤   xi  ≤ 4 ? -3 + xi :
    throw(DomainError("$xi is not in the expected domian"))
end
mccarthy(x, N=500) = sum([(1 / 2)^k * g(4^k * x) for k = 1 : N])

wen(x, N=1000) = prod([1 + (1 / 2)^n * sin(6^n * π * x) for n = 1 : N])

# --------------------------------------- Bivariate test functions --------------------------------------- #

paraboloid(x, y) = x^2 + y^2

crossintray(x, y) = -1e-4 * (abs(sin(x) * sin(y) * exp(abs(100 - sqrt(x^2 + y^2) / π))) + 1)^(0.1)

dropwave(x, y) = -(1 + cos(12 * sqrt(x^2 + y^2))) / (0.5 * (x^2 + y^2) + 2)

griewank(x, y) = (x^2 + y^2) / 4000 - cos(x) * cos(y / sqrt(2)) + 1

holder(x, y) = -abs(sin(x) * cos(y) * exp(abs(1 - sqrt(x^2 + y^2) / π)))

schewefel(x, y) = 837.9658 - x * sin(sqrt(abs(x))) - y * sin(sqrt(abs(y)))

rastrigin(x, y) = 20 + (x^2 - 10 * cos(2π * x)) + (y^2 - 10 * cos(2π * y))

levy(x, y) = sin(3π * x)^2 + (x - 1)^2 * (1 + sin(3π * y)^2) + (y - 1)^2 * (1 + sin(2π * y)^2)

eggholder(x, y) = -(y + 47) * sin(sqrt(abs(y + x / 2 + 47))) - x * sin(sqrt(abs(x - (y + 47))))

ackley(x, y, a=20, b=0.2, c=2π) =
    -a * exp(-b * sqrt(1 / 2 * (x^2 + y^2))) - exp(1 / 2 * (cos(c * x) + cos(x * y))) + a + MathConstants.e

schaffer1(x, y) = 0.5 + (sin(x^2 - y^2)^2 - 0.5) / (1 + 0.001 * (x^2 + y^2))^2

schaffer2(x, y) = 0.5 + (cos(sin(x^2 - y^2))^2 - 0.5) / (1 + 0.001 * (x^2 + y^2))^2

schubert(x, y) = sum([i * cos((i + 1) * x + i) for i in 1 : 5]) * sum([i * cos((i + 1) * y + i) for i in 1 : 5])


text:
include("src/sf_synthetic_likelihood.jl")
Qamhi_polymer(200, 2, 100, 20000, 0.03, [0.,0,0], "features.txt", "saw_conf_data20.txt")


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

mutable struct InvasiveLinkedList{T}
    # Invasive list requires that T have a field `.next >: U{T, Nothing}` and `.queue >: U{ILL{T}, Nothing}`
    head::Union{T, Nothing}
    tail::Union{T, Nothing}
    InvasiveLinkedList{T}() where {T} = new{T}(nothing, nothing)
end

#const list_append!! = append!
#const list_deletefirst! = delete!

eltype(::Type{<:InvasiveLinkedList{T}}) where {T} = @isdefined(T) ? T : Any

iterate(q::InvasiveLinkedList) = (h = q.head; h === nothing ? nothing : (h, h))
iterate(q::InvasiveLinkedList{T}, v::T) where {T} = (h = v.next; h === nothing ? nothing : (h, h))

isempty(q::InvasiveLinkedList) = (q.head === nothing)

function length(q::InvasiveLinkedList)
    i = 0
    head = q.head
    while head !== nothing
        i += 1
        head = head.next
    end
    return i
end

function list_append!!(q::InvasiveLinkedList{T}, q2::InvasiveLinkedList{T}) where T
    q === q2 && error("can't append list to itself")
    head2 = q2.head
    if head2 !== nothing
        tail2 = q2.tail::T
        q2.head = nothing
        q2.tail = nothing
        tail = q.tail
        q.tail = tail2
        if tail === nothing
            q.head = head2
        else
            tail.next = head2
        end
        while head2 !== nothing
            head2.queue = q
            head2 = head2.next
        end
    end
    return q
end

function push!(q::InvasiveLinkedList{T}, val::T) where T
    val.queue === nothing || error("val already in a list")
    val.queue = q
    tail = q.tail
    if tail === nothing
        q.head = q.tail = val
    else
        tail.next = val
        q.tail = val
    end
    return q
end

function pushfirst!(q::InvasiveLinkedList{T}, val::T) where T
    val.queue === nothing || error("val already in a list")
    val.queue = q
    head = q.head
    if head === nothing
        q.head = q.tail = val
    else
        val.next = head
        q.head = val
    end
    return q
end

function pop!(q::InvasiveLinkedList{T}) where {T}
    val = q.tail::T
    list_deletefirst!(q, val) # expensive!
    return val
end

function popfirst!(q::InvasiveLinkedList{T}) where {T}
    val = q.head::T
    list_deletefirst!(q, val) # cheap
    return val
end

function list_deletefirst!(q::InvasiveLinkedList{T}, val::T) where T
    val.queue === q || return
    head = q.head::T
    if head === val
        if q.tail::T === val
            q.head = q.tail = nothing
        else
            q.head = val.next::T
        end
    else
        head_next = head.next
        while head_next !== val
            head = head_next
            head_next = head.next::Union{T, Nothing}
        end
        if q.tail::T === val
            head.next = nothing
            q.tail = head
        else
            head.next = val.next::T
        end
    end
    val.next = nothing
    val.queue = nothing
    return q
end

#function list_deletefirst!(q::Array{T}, val::T) where T
#    i = findfirst(isequal(val), q)
#    i === nothing || deleteat!(q, i)
#    return q
#end

mutable struct LinkedListItem{T}
    # Adapter class to use any `T` in a LinkedList
    next::Union{LinkedListItem{T}, Nothing}
    queue::Union{InvasiveLinkedList{LinkedListItem{T}}, Nothing}
    value::T
    LinkedListItem{T}(value::T) where {T} = new{T}(nothing, nothing, value)
end
const LinkedList{T} = InvasiveLinkedList{LinkedListItem{T}}

# delegate methods, as needed
eltype(::Type{<:LinkedList{T}}) where {T} = @isdefined(T) ? T : Any
iterate(q::LinkedList) = (h = q.head; h === nothing ? nothing : (h.value, h))
iterate(q::InvasiveLinkedList{LLT}, v::LLT) where {LLT<:LinkedListItem} = (h = v.next; h === nothing ? nothing : (h.value, h))
push!(q::LinkedList{T}, val::T) where {T} = push!(q, LinkedListItem{T}(val))
pushfirst!(q::LinkedList{T}, val::T) where {T} = pushfirst!(q, LinkedListItem{T}(val))
pop!(q::LinkedList) = invoke(pop!, Tuple{InvasiveLinkedList,}, q).value
popfirst!(q::LinkedList) = invoke(popfirst!, Tuple{InvasiveLinkedList,}, q).value
function list_deletefirst!(q::LinkedList{T}, val::T) where T
    h = q.head
    while h !== nothing
        if isequal(h.value, val)
            list_deletefirst!(q, h)
            break
        end
        h = h.next
    end
    return q
end

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
function fire!(network, excitatory_inputs)
    for (input_node, excitatory) in zip(network.inputs, excitatory_inputs)
        input_node.io_value = excitatory
    end

    apply_to_all(fire!, network)
    apply_to_all(apply_temp_activations!, network)
end

function fire!(node)
    excitatory_output = 0
    inhibitory_output = 0

    inhibitory_activation = sum([n.inhibitory_output for n in node.in_inhibitory])
    excitatory_activation = sum([n.excitatory_output for n in node.in_excitatory])

    if inhibitory_activation <= 0.0
        excitatory_output = excitatory_activation

        if node.type == InputNode()
            excitatory_output += node.io_value
        end

        noise = 2.0*NOISE_LEVEL_NODES*rand() - NOISE_LEVEL_NODES
        excitatory_output += noise

        excitatory_output = max(excitatory_output, EXCITATORY_MIN)
        excitatory_output = min(excitatory_output, EXCITATORY_MAX)

        excitatory_output -= EXCITATORY_MIN
        excitatory_output /= (EXCITATORY_MAX - EXCITATORY_MIN)
    end

    if INHIBITORY_THRESHOLD <= excitatory_activation
        inhibitory_output = 1
    end

    if node.type == OutputNode()
        node.io_value = excitatory_output
    end

    node.temp_inhibitory_output = inhibitory_output
    node.temp_excitatory_output = excitatory_output
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
export Queue, QueueRequest, createByFamily

mutable struct QueueRequest
    name::String
    priority::Float32
    flags::vk.VkQueueFlags
    index::UInt32
end
const QUEUE_REQUEST_INDEX_NO_INDEX = typemax(UInt32)

function createGraphics(::Type{QueueRequest}, name::String)
    return QueueRequest(name, 1.0, vk.VK_QUEUE_GRAPHICS_BIT, QUEUE_REQUEST_INDEX_NO_INDEX)
end

function createTransfer(::Type{QueueRequest}, name::String)
    return QueueRequest(name, 1.0, vk.VK_QUEUE_TRANSFER_BIT, QUEUE_REQUEST_INDEX_NO_INDEX)
end

function createCompute(::Type{QueueRequest}, name::String)
    return QueueRequest(name, 1.0, vk.VK_QUEUE_COMPUTE_BIT, QUEUE_REQUEST_INDEX_NO_INDEX)
end

function createByFlags(::Type{QueueRequest}, name::String, flags::vk.VkQueueFlags, priority::Float32)
    return QueueRequest(name, priority, flags, QUEUE_REQUEST_INDEX_NO_INDEX)
end

function createByFamily(::Type{QueueRequest}, name::String, index::Integer, priority::Float32)
    return QueueRequest(name, priority, 0, index)
end

mutable struct Queue
    mFamilyIndex::UInt32
    mQueue::vk.VkQueue
    mPool::vk.VkCommandPool
    mVkDevice::vk.VkDevice

    mFencePool::Vector{vk.VkFence}
    mSubmissionFences::Vector{vk.VkFence}
    mSubmissionBuffers

    function Queue(family::UInt32, queue::vk.VkQueue, pool::vk.VkCommandPool, vkDevice::vk.VkDevice)
        this = new()
        this.mFamilyIndex = family
        this.mQueue = queue
        this.mPool = pool
        this.mVkDevice = vkDevice

        this.mFencePool = Vector{vk.VkFence}()
        this.mSubmissionFences = Vector{vk.VkFence}()
        this.mSubmissionBuffers = []
        return this
    end
end

@class Queue

function handle(this::Queue)::vk.VkQueue
    return this.mQueue
end

function family(this::Queue)::UInt32
    return this.mFamilyIndex
end

function getPool(this::Queue)::vk.VkCommandPool
    return this.mPool
end

function getLogicalDevice(this::Queue)::vk.VkDevice
    return this.mVkDevice
end

function catchUp(this::Queue, inflightBuffers::Int32)
    while(length(this.mSubmissionFences) > inflightBuffers)
        VkExt.waitForFences(this.mVkDevice, this.mSubmissionFences, VkExt.VK_FALSE, typemax(UInt64))
        gc(this)
    end
end

function activeBuffers(this::Queue)::UInt32
    this.gc()
    return length(this.mSubmissionFences)
end

function gc(this::Queue)
    # Garbage collection of pool
    len = length(this.mSubmissionFences)
    i = 1
    while i <= len
        fence = this.mSubmissionFences[i]
        if vk.vkGetFenceStatus(this.mVkDevice, fence) == vk.VK_SUCCESS
            push!(this.mFencePool, fence)
            deleteat!(this.mSubmissionFences, i)
            freeCommandBuffer(this.mSubmissionBuffers[i])
            deleteat!(this.mSubmissionBuffers, i)
            len -= 1
        else
            i += 1
        end
    end
end

function submit(this::Queue, cmd,
                  waitSemaphores::Vector{vk.VkSemaphore},
                      waitStages::Vector{vk.VkPipelineStageFlags},
                signalSemaphores::Vector{vk.VkSemaphore})
    fence = findFreeFence(this)
    push!(this.mSubmissionFences, fence)
    push!(this.mSubmissionBuffers, cmd)
    cmds = [handle(cmd)]
    info = Ref(vk.VkSubmitInfo(
        vk.VK_STRUCTURE_TYPE_SUBMIT_INFO, #sType::VkStructureType
        C_NULL, #pNext::Ptr{Cvoid}
        length(waitSemaphores), #waitSemaphoreCount::UInt32
        pointer(waitSemaphores), #pWaitSemaphores::Ptr{VkSemaphore}
        pointer(waitStages), #pWaitDstStageMask::Ptr{VkPipelineStageFlags}
        1, #commandBufferCount::UInt32
        pointer(cmds), #pCommandBuffers::Ptr{VkCommandBuffer}
        length(signalSemaphores), #signalSemaphoreCount::UInt32
        pointer(signalSemaphores) #pSignalSemaphores::Ptr{VkSemaphore}
    ))

    GC.@preserve cmds info begin
        if vk.vkQueueSubmit(this.mQueue, 1, info, fence) != vk.VK_SUCCESS
            error("Failed to submit draw command buffer!")
        end
    end
end

 function findFreeFence(this::Queue)::vk.VkFence
    if !isempty(this.mFencePool)
        ret = pop!(this.mFencePool)
        vk.vkResetFences(this.mVkDevice, 1, Ref(ret))
        return ret
    else
        ret = VkExt.createFence(this.mVkDevice)
        return ret
    end
 end


text:
## Literals

# Literals look like 3, i.e. x₃ or -25, i.e. ¬x₂₅
# truth_index(literal) = Int(literal > 0)  # positive is true, negative is false

# assignments are like [1, -2, 3, -4]

debug(kw) = haskey(kw, :debug) && kw[:debug]

const unassigned = 0

index(literal) = abs(literal)   # which variable a literal is: index(-3) == 3
is_unassigned(literal) = literal == unassigned
is_assigned(literal) = literal != unassigned

is_positive(literal) = literal > 0

make_literal(variable, sign) = copysign(variable, sign)

num_assigned(assignments) = count(is_assigned, assignments)

struct Action
    action::Symbol
    literal::Int
end

Base.show(io::IO, action::Action) =
    print(io, "$(String(action.action)[1]):$(action.literal)")

const actions = Action[]
sizehint!(actions, 1000)

const assignments = Int[]

const unit_literals = Int[]

"""
Process a clause to check sat or find next unassigned variable
`assignments` has -1 if unassigned, 0 or 1 if assigned and false/true respectively

Clause looks like [1, 3, -25]

Returns :sat if the clause is not falsified, :unsat if it is falsified,
:unit if it found a single unassigned literal (i.e. a unit literal)

This assumes that there are no clauses with both x and ¬x in

TODO: Deal with x ∨ ¬x clauses in pre-processing
"""
function process(clause)

    unit_literal = 0

    for literal in clause

        variable = index(literal)
        current = assignments[variable]

        if is_unassigned(current)

            if unit_literal != 0   # already found a previous unit literal in the clause
                return :sat, -1  # > 1 literal unassigned in clause, so not falsified
            end

            unit_literal = literal

              # return :unassigned, literal
         end

        # if any literal in clause has correct value, then clause is satisfiable
        if current == literal
            return :sat, -1
        end
    end

    # exactly one is unassigned
    if unit_literal != 0
        return :unit, unit_literal
    end

    # no literals are satisfiable, hence the clause is unsatisfiable
    return :unsat, -1

end

indent(level) = print(" " ^ level)

function check_clauses(p, variable; kw...)
    for clause_number in p.clause_list[variable]
        clause = p.clauses[clause_number]

        status, literal = process(clause)

        # if variable == 12
        #     @show clause, status, literal
        # end

        if status == :unit
            push!(unit_literals, literal)
        end

        if status == :unsat
            return :unsat, clause
        end

    end

    return :sat, Int[]
end

"Unassign actions back to given literal"
function unassign!(original_literal)
    while true
        action = actions[end]
        literal = action.literal

        if literal == original_literal
            # finished
            return
        end

        variable = index(literal)
        assignments[variable] = unassigned

        pop!(actions)  # remove the action we just processed

    end

end

"""Do unit propagation for the literals implied by `original_literal`.
These were previously collected by `check_clauses`.
"""
function unit_propagation(p, original_literal; kw...)

    if isempty(unit_literals)
        return
    end

    status = :sat

    if debug(kw)
        @info "Doing unit prop"
        @show assignments
        @show original_literal
    end

    while !isempty(unit_literals)

        literal = popfirst!(unit_literals)
        variable = index(literal)
        push!(actions, Action(:unitprop, literal))

        assignments[variable] = literal

        status, clause = check_clauses(p, variable; kw...)

        if status == :unsat
            if debug(kw)
                @info "Unsat: ", clause
            end
            break
        end

        # no clauses were falsified

    end

    if status == :unsat
        unassign!(original_literal)
    end

    empty!(unit_literals)

    return status
end

function assign!(p, literal; kw...)

    variable = index(literal)

    push!(actions, Action(:assign, literal))
    assignments[variable] = literal

    status, clause = check_clauses(p, variable; kw...)

    if status == :unsat

        if debug(kw)
            @info "Unsat clause: ", clause
        end

        empty!(unit_literals)
        return :unsat
    end

    status = unit_propagation(p, literal; kw...)  # remove this line to remove unit prop

    return status

end

function select_variable(assignments)

    ## random choice:
    # possible = (1:length(assignments))[is_unassigned.(assignments)]
    # return rand(possible)

    ## first available:
    variable = findfirst(is_unassigned, assignments)
end

iterative_solve(p; kw...) =
    iterative_solve(StructuredSATProblem(p); kw...)

"""Solve a SAT problem by tree search.
- `assignments` specifies if each variable is unassigned (-1) or assigned with index false (0) or true (1)
- `starting_clause` says which clauses can be skipped since they are satisfied by
the current set of assignments
"""
function iterative_solve(p::StructuredSATProblem; kw...)

    empty!(actions)
    empty!(assignments)

    append!(assignments, fill(unassigned, p.num_variables))

    backtrack = false

    while true

        if debug(kw)
            # decisions = filter(x -> x.action == :assign, actions)
            # @show decisions

            @show actions

            assigned = filter(is_assigned, assignments)
            @show assigned
        end

        if !backtrack

            if num_assigned(assignments) == p.num_variables
                # finished!
                return :sat, assignments
            end

            # choose and assign a new variable:
            variable = select_variable(assignments)
            literal = make_literal(variable, 1)

            status = assign!(p, literal; kw...)

            if status == :unsat
                backtrack = true
            end

        else # backtrack
            if debug(kw)
                @info "Backtracking..."
            end

            if isempty(actions)
                # @error "Action list empty -- should not happen"
                break
            end

            # undo the last action
            action = pop!(actions)

            if action.action == :unitprop   # keep undoing unit props
                variable = index(action.literal)
                assignments[variable] = unassigned
                continue
            end

            literal = action.literal

            variable = index(literal)

            if is_positive(literal)
                # switch positive to negative and try again
                literal = make_literal(variable, -1)

                status = assign!(p, literal; kw...)

                if status == :sat
                   backtrack = false  # stop backtracking
                end

            else
                # neither assignment to variable worked
                # keep backtracking
                assignments[variable] = unassigned
            end

        end

    end

    return :unsat, assignments

end

max_stars_repo_name: dpsanders/SimpleSATSolver

max_issues_repo_name: dpsanders/SimpleSATSolver

max_forks_repo_name: dpsanders/SimpleSATSolver


text:
"Load a `NULL`-terminated list of strings"
function loadstringlist(pstringlist::Ptr{Cstring})
    stringlist = Vector{ASCIIString}()
    (pstringlist == C_NULL) && return stringlist
    i = 1
    item = unsafe_load(pstringlist, i)
    while Ptr{UInt8}(item) != C_NULL
        push!(stringlist, bytestring(item))
        i += 1
        item = unsafe_load(pstringlist, i)
    end
    stringlist
end

"Fetch list of (non-empty) metadata domains. (Since: GDAL 1.11)"
getmetadatadomainlist{T <: GDAL.GDALMajorObjectH}(obj::Ptr{T}) =
    loadstringlist(GDAL.C.GDALGetMetadataDomainList(obj))

"""
Fetch metadata. Use `""` or `NULL` for the default domain.

The returned string list is owned by the object, and may change at any time.
It is formated as a "Name=value" list with the last pointer value being `NULL`.
Note that relatively few formats return any metadata at this time.
"""
getmetadata{T <: GDAL.GDALMajorObjectH}(obj::Ptr{T}; domain::AbstractString = "") =
    loadstringlist(GDAL.C.GDALGetMetadata(obj, pointer(domain)))

# """
#     GDALSetMetadata(GDALMajorObjectH hObject,
#                     char ** papszMD,
#                     const char * pszDomain) -> CPLErr
# Set metadata.
# """
# function setmetadata{T <: GDALMajorObjectH}(arg1::Ptr{T},arg2,arg3)
#     ccall((:GDALSetMetadata,libgdal),CPLErr,(Ptr{GDALMajorObjectH},Ptr{Cstring},Cstring),arg1,arg2,arg3)
# end

# """
#     GDALGetMetadataItem(GDALMajorObjectH hObject,
#                         const char * pszName,
#                         const char * pszDomain) -> const char *
# Fetch single metadata item.
# """
# function getmetadataitem{T <: GDALMajorObjectH}(arg1::Ptr{T},arg2,arg3)
#     bytestring(ccall((:GDALGetMetadataItem,libgdal),Cstring,(Ptr{GDALMajorObjectH},Cstring,Cstring),arg1,arg2,arg3))
# end

# """
#     GDALSetMetadataItem(GDALMajorObjectH hObject,
#                         const char * pszName,
#                         const char * pszValue,
#                         const char * pszDomain) -> CPLErr
# Set single metadata item.
# """
# function setmetadataitem{T <: GDALMajorObjectH}(arg1::Ptr{T},arg2,arg3,arg4)
#     ccall((:GDALSetMetadataItem,libgdal),CPLErr,(Ptr{GDALMajorObjectH},Cstring,Cstring,Cstring),arg1,arg2,arg3,arg4)
# end

# """
#     GDALGetDescription(GDALMajorObjectH hObject) -> const char *
# Fetch object description.
# """
# function getdescription{T <: GDALMajorObjectH}(arg1::Ptr{T})
#     bytestring(ccall((:GDALGetDescription,libgdal),Cstring,(Ptr{GDALMajorObjectH},),arg1))
# end

# """
#     GDALSetDescription(GDALMajorObjectH hObject,
#                        const char * pszNewDesc) -> void
# Set object description.
# """
# function setdescription{T <: GDALMajorObjectH}(arg1::Ptr{T},arg2)
#     ccall((:GDALSetDescription,libgdal),Void,(Ptr{GDALMajorObjectH},Cstring),arg1,arg2)
# end


text:
#=

This module provides code to compute Ramsey equilibria in a LQ economy with
distortionary taxation.  The program computes allocations (consumption,
leisure), tax rates, revenues, the net present value of the debt and other
related quantities.

Functions for plotting the results are also provided below.

@author : Spencer Lyon <spencer.lyon@nyu.edu>

@date: 2014-08-21

References
----------

Simple port of the file examples/lqramsey.py

http://quant-econ.net/lqramsey.html

=#
using QuantEcon
using PyPlot

abstract AbstractStochProcess

type ContStochProcess <: AbstractStochProcess
    A::Matrix
    C::Matrix
end

type DiscreteStochProcess <: AbstractStochProcess
    P::Matrix
    x_vals::Array
end

type Economy{SP <: AbstractStochProcess}
    bet::Real
    Sg::Matrix
    Sd::Matrix
    Sb::Matrix
    Ss::Matrix
    is_discrete::Bool
    proc::SP
end

type Path
    g
    d
    b
    s
    c
    l
    p
    tau
    rvn
    B
    R
    pi
    Pi
    xi
end

function compute_exog_sequences(econ::Economy, x)
    # Compute exogenous variable sequences
    Sg, Sd, Sb, Ss = econ.Sg, econ.Sd, econ.Sb, econ.Ss
    g, d, b, s = [squeeze(S * x, 1) for S in (Sg, Sd, Sb, Ss)]

    #= Solve for Lagrange multiplier in the govt budget constraint
    In fact we solve for nu = lambda / (1 + 2*lambda).  Here nu is the
    solution to a quadratic equation a(nu^2 - nu) + b = 0 where
    a and b are expected discounted sums of quadratic forms of the state. =#
    Sm = Sb - Sd - Ss

    return g, d, b, s, Sm
end

function compute_allocation(econ::Economy, Sm, nu, x, b)
    Sg, Sd, Sb, Ss = econ.Sg, econ.Sd, econ.Sb, econ.Ss

    # Solve for the allocation given nu and x
    Sc = 0.5 .* (Sb + Sd - Sg - nu .* Sm)
    Sl = 0.5 .* (Sb - Sd + Sg - nu .* Sm)
    c = squeeze(Sc * x, 1)
    l = squeeze(Sl * x, 1)
    p = squeeze((Sb - Sc) * x, 1)  # Price without normalization
    tau = 1 .- l ./ (b .- c)
    rvn = l .* tau

    return Sc, Sl, c, l, p, tau, rvn
end

function compute_nu(a0, b0)
    disc = a0^2 - 4a0*b0

    if disc >= 0
        nu = 0.5 *(a0 - sqrt(disc)) / a0
    else
        println("There is no Ramsey equilibrium for these parameters.")
        error("Government spending (economy.g) too low")
    end

    # Test that the Lagrange multiplier has the right sign
    if nu * (0.5 - nu) < 0
        print("Negative multiplier on the government budget constraint.")
        error("Government spending (economy.g) too low")
    end

    return nu
end

function compute_Pi(B, R, rvn, g, xi)
    pi = B[2:end] - R[1:end-1] .* B[1:end-1] - rvn[1:end-1] + g[1:end-1]
    Pi = cumsum(pi .* xi)
    return pi, Pi
end

function compute_paths(econ::Economy{DiscreteStochProcess}, T)
    # simplify notation
    bet, Sg, Sd, Sb, Ss = econ.bet, econ.Sg, econ.Sd, econ.Sb, econ.Ss
    P, x_vals = econ.proc.P, econ.proc.x_vals

    state = mc_sample_path(P, 1, T)
    x = x_vals[:, state]

    # Compute exogenous sequence
    g, d, b, s, Sm = compute_exog_sequences(econ, x)

    # compute a0, b0
    ns = size(P, 1)
    F = eye(ns) - bet.*P
    a0 = (F \ ((Sm * x_vals)'.^2))[1] ./ 2
    H = ((Sb - Sd + Sg) * x_vals) .* ((Sg - Ss)*x_vals)
    b0 = (F \ H')[1] ./ 2

    # compute lagrange multiplier
    nu = compute_nu(a0, b0)

    # Solve for the allocation given nu and x
    Sc, Sl, c, l, p, tau, rvn = compute_allocation(econ, Sm, nu, x, b)

    # compute remaining variables
    H = ((Sb - Sc)*x_vals) .* ((Sl - Sg)*x_vals) - (Sl*x_vals).^2
    temp = squeeze(F*H', 2)
    B = temp[state] ./ p
    H = squeeze(P[state, :] * ((Sb - Sc)*x_vals)', 2)
    R = p ./ (bet .* H)
    temp = squeeze(P[state, :] *((Sb - Sc) * x_vals)', 2)
    xi = p[2:end] ./ temp[1:end-1]

    # compute pi
    pi, Pi = compute_Pi(B, R, rvn, g, xi)

    Path(g, d, b, s, c, l, p, tau, rvn, B, R, pi, Pi, xi)
end

function compute_paths(econ::Economy{ContStochProcess}, T)
    # simplify notation
    bet, Sg, Sd, Sb, Ss = econ.bet, econ.Sg, econ.Sd, econ.Sb, econ.Ss
    A, C = econ.proc.A, econ.proc.C

    # Generate an initial condition x0 satisfying x0 = A x0
    nx, nx = size(A)
    x0 = null((eye(nx) - A))
    x0 = x0[end] < 0 ? -x0 : x0
    x0 = x0 ./ x0[end]
    x0 = squeeze(x0, 2)

    # Generate a time series x of length T starting from x0
    nx, nw = size(C)
    x = zeros(nx, T)
    w = randn(nw, T)
    x[:, 1] = x0
    for t=2:T
        x[:, t] = A *x[:, t-1] + C * w[:, t]
    end

    # Compute exogenous sequence
    g, d, b, s, Sm = compute_exog_sequences(econ, x)

    # compute a0 and b0
    H = Sm'Sm
    a0 = 0.5 * var_quadratic_sum(A, C, H, bet, x0)
    H = (Sb - Sd + Sg)'*(Sg + Ss)
    b0 = 0.5 * var_quadratic_sum(A, C, H, bet, x0)

    # compute lagrange multiplier
    nu = compute_nu(a0, b0)

    # Solve for the allocation given nu and x
    Sc, Sl, c, l, p, tau, rvn = compute_allocation(econ, Sm, nu, x, b)

    # compute remaining variables
    H = Sl'Sl - (Sb - Sc)' *(Sl - Sg)
    L = Array(Float64, T)
    for t=1:T
        L[t] = var_quadratic_sum(A, C, H, bet, x[:, t])
    end
    B = L ./ p
    Rinv = squeeze(bet .* (Sb- Sc)*A*x, 1) ./ p
    R = 1 ./ Rinv
    AF1 = (Sb - Sc) * x[:, 2:end]
    AF2 = (Sb - Sc) * A * x[:, 1:end-1]
    xi =  AF1 ./ AF2
    xi = squeeze(xi, 1)

    # compute pi
    pi, Pi = compute_Pi(B, R, rvn, g, xi)

    Path(g, d, b, s, c, l, p, tau, rvn, B, R, pi, Pi, xi)
end

function gen_fig_1(path::Path)
    T = length(path.c)

    num_rows, num_cols = 2, 2
    fig, axes = subplots(num_rows, num_cols, figsize=(14, 10))
    plt.subplots_adjust(hspace=0.4)
    for i=1:num_rows
        for j=1:num_cols
            axes[i, j][:grid]()
            axes[i, j][:set_xlabel]("Time")
        end
    end

    bbox = (0., 1.02, 1., .102)
    legend_args = {:bbox_to_anchor => bbox, :loc => 3, :mode => :expand}
    p_args = {:lw => 2, :alpha => 0.7}

    # Plot consumption, govt expenditure and revenue
    ax = axes[1, 1]
    ax[:plot](path.rvn, label=L"$\tau_t \ell_t$"; p_args...)
    ax[:plot](path.g, label=L"$g_t$"; p_args...)
    ax[:plot](path.c, label=L"$c_t$"; p_args...)
    ax[:legend](ncol=3; legend_args...)

    # Plot govt expenditure and debt
    ax = axes[1, 2]
    ax[:plot](1:T, path.rvn, label=L"$\tau_t \ell_t$"; p_args...)
    ax[:plot](1:T, path.g, label=L"$g_t$"; p_args...)
    ax[:plot](1:T-1, path.B[2:end], label=L"$B_{t+1}$"; p_args...)
    ax[:legend](ncol=3; legend_args...)

    # Plot risk free return
    ax = axes[2, 1]
    ax[:plot](1:T, path.R - 1, label=L"$R_{t - 1}$"; p_args...)
    ax[:legend](ncol=1; legend_args...)

    # Plot revenue, expenditure and risk free rate
    ax = axes[2, 2]
    ax[:plot](1:T, path.rvn, label=L"$\tau_t \ell_t$"; p_args...)
    ax[:plot](1:T, path.g, label=L"$g_t$"; p_args...)
    ax[:plot](1:T-1, path.pi, label=L"$\pi_{t+1}$"; p_args...)
    ax[:legend](ncol=3; legend_args...)

    plt.show()

end

function gen_fig_2(path::Path)
    T = length(path.c)

    # Prepare axes
    num_rows, num_cols = 2, 1
    fig, axes = subplots(num_rows, num_cols, figsize=(10, 10))
    plt.subplots_adjust(hspace=0.5)
    bbox = (0., 1.02, 1., .102)
    legend_args = {:bbox_to_anchor => bbox, :loc => 3, :mode => :expand}
    p_args = {:lw => 2, :alpha => 0.7}

    # Plot adjustment factor
    ax = axes[1]
    ax[:plot](2:T, path.xi, label=L"$\xi_t$"; p_args...)
    ax[:grid]()
    ax[:set_xlabel]("Time")
    ax[:legend](ncol=1; legend_args...)

    # Plot adjusted cumulative return
    ax = axes[2]
    ax[:plot](2:T, path.Pi, label=L"$\Pi_t$"; p_args...)
    ax[:grid]()
    ax[:set_xlabel]("Time")
    ax[:legend](ncol=1; legend_args...)

    plt.show()
end