text:
@testset "DataFrame constructor" begin
    df = DataFrame(x=1:10, y=sin.(1:10))
    p = Plot(df, x=:x, y=:y)
    @test p isa Plot
    @test length(p.data) == 1
    @test p.layout.xaxis_title == Dict(:text => :x)
    @test p.layout.yaxis_title == Dict(:text => :y)
end

@testset "overriding yaxis_title" begin
    df = DataFrame(x=1:10, y=sin.(1:10))
    p = Plot(df, x=:x, y=:y, Layout(yaxis_title="sin(x)"));
    @test p.layout.xaxis_title == Dict(:text => :x)
    @test p.layout.yaxis_title == Dict(:text => "sin(x)")
end

@testset "overriding yaxis_title_text" begin
    df = DataFrame(x=1:10, y=sin.(1:10))
    p = Plot(df, x=:x, y=:y, Layout(yaxis_title_text="sin(x)"));
    @test p.layout.xaxis_title == Dict(:text => :x)
    @test p.layout.yaxis_title == Dict(:text => "sin(x)")
end

@testset "overriding yaxis_title with attr(text=, x=)" begin
    df = DataFrame(x=1:10, y=sin.(1:10))
    p = Plot(df, x=:x, y=:y, Layout(yaxis_title=attr(x=0.5, text="sin(x)")));
    @test p.layout.xaxis_title == Dict(:text => :x)
    @test p.layout.yaxis_title_text == "sin(x)"
    @test p.layout.yaxis_title == Dict(:text => "sin(x)", :x => 0.5)
end

@testset "facet_row_wrap and facet_col_wrap" begin
    df = stack(DataFrame(x=1:10, one=1, two=2, three=3, four=4, five=5, six=6, seven=7), Not(:x))
    # _43 means 4 subplots in row 1 and 3 in row 2
    p_43 = Plot(df, x=:x, y=:value, facet_row=:variable, facet_row_wrap=2)
    @test size(p_43.layout.subplots.grid_ref) == (2, 4)

    p_2221 = Plot(df, x=:x, y=:value, facet_col=:variable, facet_col_wrap=2)
    @test size(p_2221.layout.subplots.grid_ref) == (4, 2)

    p_7x7 = @test_logs (:warn, "cannot set facet_row and facet_col_wrap -- setting facet_col_wrap to missing") Plot(df, x=:x, y=:value, facet_col=:variable, facet_col_wrap=2, facet_row=:value)
    @test size(p_7x7.layout.subplots.grid_ref) == (7, 7)

end

@testset "splom and dimensions" begin
    df = DataFrame(d1=rand(10), d2=rand(10), d3=rand(10), color_col=vcat(fill("c1", 5), fill("c2", 5)))
    p1 = Plot(df, dimensions=[:d1, :d2], kind="splom")
    @test length(p1.data) == 1
    trace = p1.data[1]
    @test !trace.diagonal_visible
    @test length(trace.dimensions) == 2
    @test trace.dimensions[1].label == :d1
    @test trace.dimensions[1].values == df[!, :d1]
    @test trace.dimensions[2].label == :d2
    @test trace.dimensions[2].values == df[!, :d2]
    @test all(x -> x.axis_matches, trace.dimensions)

    p2 = Plot(df, dimensions=[:d1, :d2], kind="splom", color=:color_col)
    @test length(p2.data) == 2
    for (trace, inds) in zip(p2.data, [1:5, 6:10])
        @test !trace.diagonal_visible
        @test length(trace.dimensions) == 2
        @test trace.dimensions[1].label == :d1
        @test trace.dimensions[1].values == df[inds, :d1]
        @test trace.dimensions[2].label == :d2
        @test trace.dimensions[2].values == df[inds, :d2]
        @test all(x -> x.axis_matches, trace.dimensions)
    end
end


text:
using OnlineStatsBase, OrderedCollections, StatsBase, Statistics, Dates, Test

import OnlineStatsBase: _fit!, _merge!

O = OnlineStatsBase

include("test_stats.jl")

#-----------------------------------------------------------------------# Random
@testset "Random Stuff" begin
string(Series(Mean(), Variance()))
string(Group(Series(Mean(), Variance()), Extrema()))
string(Mean())
string(McclainWeight())
@test Mean() != Variance()
@test !(Mean() == Variance())
@test Mean() == merge(Mean(), Mean())
@test Mean() == merge(Mean(), fit!(Variance(), 1:5))
@test_throws Exception fit!(Mean(), "abc")
@test_throws Exception OnlineStatsBase._fit!(Mean(), "a")
@test collect(OnlineStatsBase.neighbors([1,3,5])) == [(1,3), (3,5)]
@test isnan(value(fit!(Variance(), NaN)))
end

@testset "Broadcasting" begin
    o1 = Mean()
    o2 = Variance()
    @test tuple.(o1, o2, [1, 2]) == [(o1, o2, 1), (o1, o2, 2)]
    @test tuple.(o1, o2, [1, 2])[1][1] === o1
    @test tuple.(o1, o2, [1, 2])[2][1] === o1
end

#-----------------------------------------------------------------------# Weight
@testset "Weight" begin
function test_weight(w::OnlineStatsBase.Weight, f::Function)
    println("  > $w")
    @test w == copy(w)
    for i in 1:20
        @test w(i) == f(i)
    end
end
test_weight(@inferred(EqualWeight()),               i -> 1 / i)
test_weight(@inferred(ExponentialWeight(.1)),       i -> .1)
test_weight(@inferred(LearningRate(.6)),            i -> 1 / i^.6)
test_weight(@inferred(LearningRate2(.5)),           i -> 1 / (1 + .5*(i-1)))
test_weight(@inferred(HarmonicWeight(4.)),          i -> 4 / (4 + i - 1))

@test ExponentialWeight(20) == ExponentialWeight(2 / 21)

@testset "McclainWeight" begin
    w = McclainWeight(.1)
    for i in 2:100
        @test .1 < w(i) < 1
    end
end
@testset "first weight is one" begin
    for w in [EqualWeight(), LearningRate(), LearningRate2(), HarmonicWeight(), McclainWeight()]
        @test w(1) == 1
    end
end
end #Weight

ext: jl


text:
"""
    Path

Vector of tuples `(t,v)` corresponding to a timed path through a graph.
"""
const Path = Vector{Tuple{Int,Int}}  # vector of tuples (t, v)

"""
    Solution

Vector of `Path`s, one for each agent of a [`MAPF`](@ref) problem.
"""
const Solution = Vector{Path}  # one path for each agent a

"""
    Reservation

Set of tuples `(t,v)` corresponding to already-occupied vertices.
"""
const Reservation = Set{Tuple{Int,Int}}

"""
    MAPF{G,EW}
"""
struct MAPF{G<:AbstractGraph{Int}}
    graph::G
    rev_graph::G
    edge_indices::Dict{Tuple{Int,Int},Int}
    rev_edge_indices::Dict{Tuple{Int,Int},Int}
    edge_weights::Vector{Float64}
    edge_weights_mat::SparseMatrixCSC{Float64,Int64}
    sources::Vector{Int}
    destinations::Vector{Int}
    starting_times::Vector{Int}
    distances_to_destinations::Dict{Int,Vector{Float64}}
    conflict_groups::Vector{Vector{Int}}
    group_memberships::Vector{Vector{Int}}
end

function MAPF(;
    graph, sources, destinations, starting_times, conflict_groups=[[v] for v in 1:nv(graph)]
)
    @assert is_directed(graph)
    rev_graph = reverse(graph)

    edge_indices = Dict((src(ed), dst(ed)) => e for (e, ed) in enumerate(edges(graph)))
    rev_edge_indices = Dict((dst(ed), src(ed)) => e for (e, ed) in enumerate(edges(graph)))

    edge_weights_mat = Graphs.weights(graph)
    edge_weights = [edge_weights_mat[src(ed), dst(ed)] for ed in edges(graph)]

    distances_to_destinations = Dict{Int,Vector{Float64}}()
    for d in unique(destinations)
        dijkstra_state = my_dijkstra_shortest_paths(
            rev_graph, d; edge_indices=rev_edge_indices, edge_weights=edge_weights
        )
        distances_to_destinations[d] = dijkstra_state.dists
    end

    group_memberships = [Int[] for v in vertices(graph)]
    for (g, group) in enumerate(conflict_groups)
        for v in group
            push!(group_memberships[v], g)
        end
    end

    return MAPF(
        graph,
        rev_graph,
        edge_indices,
        rev_edge_indices,
        edge_weights,
        edge_weights_mat,
        sources,
        destinations,
        starting_times,
        distances_to_destinations,
        conflict_groups,
        group_memberships,
    )
end

nb_agents(mapf::MAPF) = length(mapf.sources)


text:
module ElectrochemicalKinetics

include("dos.jl")
export DOS
using .DOS: DOSData, get_dos
export DOSData, get_dos

include("kinetic_models.jl")
export ButlerVolmer, AsymptoticMarcusHushChidsey
export fermi_dirac
export integrand
export KineticModel, IntegralModel
export Marcus, MarcusHushChidsey, MarcusHushChidseyDOS

include("quantum_capacitance.jl")

include("rate_constant.jl")
export compute_k, compute_k_cq

include("fitting.jl")
include("datastructures.jl")
export fitting_params, fit_model, fit_overpotential, is_dosmodel

include("plot_fcns.jl")
export plot_models, plot_exp_and_models

end


text:
"""
Function to fetch and pool the covariances matrices.
Calls the function read_covarmat to fetch and pool_covarmats to pool the matrices.
Arguments:
myscratch::Boostscratch -> Container with data which are relevant to perform the boostinsteps
paths_in_covarmat::Array{ASCIIString, 1} -> Array with the paths which lead to the provided covariance matrices
"""
function fetch_pool_covarmat!(myscratch::Boostscratch, paths_in_covarmat::Array{ASCIIString, 1})
	covarmats = Array{Covarmat, 1}(length(paths_in_covarmat))
	for i = 1 : length(paths_in_covarmat)
		covarmats[i] = read_covarmat(paths_in_covarmat[i])
	end
	if(isempty(myscratch.pooledcovarmat.covarmat))
		myscratch.pooledcovarmat.covarmat = pool_covarmats(myscratch.pooledunibeta, covarmats)
		myscratch.pooledcovarmat.labels = myscratch.wantedlabels
	else
		pooled = pool_covarmats(myscratch.pooledunibeta, covarmats)
		merge_covarmats!(myscratch, pooled)
	end
end
"""
Function to fetch and pool the univariate estimators.
Calls the function read_unibeta to fetch and pool_unibetas to pool the estimators.
Arguments:
myscratch::Boostscratch -> Container with data which are relevant to perform the boostinsteps
paths_in_unibeta::Array{ASCIIString, 1} -> Array with the paths which lead to the provided univariate estimators
"""
function fetch_pool_unibeta!(myscratch::Boostscratch, paths_in_unibeta::Array{ASCIIString,1})
	unibetas = Array{Unibeta, 1}(length(paths_in_unibeta))
	for i = 1 : length(paths_in_unibeta)
		unibetas[i] = read_unibeta(paths_in_unibeta[i])
	end
	myscratch.pooledunibeta = pool_unibetas(unibetas)
	myscratch.actualbeta = zeros(length(myscratch.pooledunibeta.unibeta))
	myscratch.actualnom = myscratch.pooledunibeta.unibeta
end

max_stars_repo_name: danielazoeller/DistributedBoosting

max_issues_repo_name: danielazoeller/DistributedBoosting

max_forks_repo_name: danielazoeller/DistributedBoosting


text:
include("MultiFidelityABC.jl")
mkpath("figures")

using StatsPlots, Random

println("#### Repressilator")
println("# Loading data")
bm = MakeBenchmarkCloud("repressilator/output")
epsilons = (50.0,50.0)
sample_size = 10^4

println("# Fig 1")
fig1a = view_distances(bm[1:sample_size], epsilons)
fig1b = view_distances(bm[1:sample_size], epsilons, 2, L"n")
savefig(plot(fig1a, fig1b, layout=2, size=(900,360)), "figures/fig1.pdf")

println("# Fig 2 and Table 1")
fig2a = compare_efficiencies(bm, sample_size, epsilons, output="theory")
fig2b, table1, latex_table1 = compare_efficiencies(bm, sample_size, epsilons, output="plot")
savefig(plot(fig2a, fig2b, layout=2, size=(1100,440)), "figures/fig2.pdf")

println("# Table 2 and 3")
eta_tab, phi_tab = get_eta(bm, epsilons, Repressilator.F)
Random.seed!(123)
var_tab = variance_table(bm, 10^3, epsilons, eta_tab, Repressilator.F, 30.0)
Random.seed!()

println("# Supplementary")

t,x,p = simulate(Repressilator.tlm)
tt,xx = complete(Repressilator.tlm, p)

function show_plots(solutions::NTuple{N,Tuple{Times,States,String}}, j::Integer; kwargs...) where N
       fig = plot(; kwargs...);
       for (t,x,l) in solutions
       plot!(t,x[j,:]; label=l)
       end
       return fig
end

titles = ["mRNA1", "mRNA2", "mRNA3", "P1", "P2", "P3"]
figs = [show_plots(((t,x,"Tau-leap"),(tt,xx,"Gillespie")),j; title=ttl, legend=:none) for (j, ttl) in enumerate(titles)]
figs[1] = plot!(figs[1]; legend=:topleft)
repressilator_eg = plot(figs..., layout=6)
savefig(repressilator_eg, "./figures/SFig1.pdf")

println("#### Viral")
println("# Load data")

bm = MakeBenchmarkCloud("viral/output")
mf_smallBI_location = "./viral/output/mf_smallBI/"
mf_largeBI_location = "./viral/output/mf_largeBI/"

epsilons = (0.25,0.25)
eta_0 = 0.01
smallBI_size = 10^3
largeBI_size = length(bm)

function divide_cloud(c::MFABCCloud, s::Integer; stage::String)
    if stage=="bm"
        return c[1:s]
    elseif stage=="inc"
        return c[s+1:end]
    end
    error("What stage? bm or inc")
end

bm_set = Array{MFABCCloud,1}()
mf_set = Array{MFABCCloud,1}()
inc_smallBI_set = Array{MFABCCloud,1}()
inc_largeBI_set = Array{MFABCCloud,1}()

for cloud_location in mf_smallBI_location.*readdir(mf_smallBI_location)
    c = MakeMFABCCloud(cloud_location)
    push!(mf_set, c)
    push!(bm_set, divide_cloud(c, smallBI_size, stage="bm"))
    push!(inc_smallBI_set, divide_cloud(c, smallBI_size, stage="inc"))
end
for cloud_location in mf_largeBI_location.*readdir(mf_largeBI_location)
    c = MakeMFABCCloud(cloud_location)
    push!(inc_largeBI_set, divide_cloud(c, largeBI_size, stage="inc"))
end

println("# Fig 3")
savefig(view_distances(bm[1:10000], epsilons, epsilons.*2), "figures/fig3.pdf")

println("# Fig 4")
fig4 = plot_eta_estimates(bm, epsilons, (bm_set,"After burn-in"), (mf_set,"After adaptation"); method="mf", lower_eta=eta_0)
plot!(xlim=(0,0.4),ylim=(0,0.2))
savefig(fig4, "figures/fig4.pdf")

println("# Fig 5")
savefig(plot_apost_efficiencies((inc_largeBI_set,"After large burn-in"), (inc_smallBI_set, "After small burn-in"), (bm_set, "During burn-in")), "figures/fig5.pdf")

println("# Supplementary")
t,x,p = simulate(Repressilator.tlm)
tt,xx = complete(Repressilator.tlm, p)

f = [plot(t,x',label=["template" "genome" "struct" "virus"]), plot(tt,xx',legend=:none)]
savefig(plot(f..., layout=2, size=(1000,400)), "./figures/SFig2.pdf")

max_stars_repo_name: tpprescott/mf-abc

max_issues_repo_name: tpprescott/mf-abc

max_forks_repo_name: tpprescott/mf-abc


text:
# using Revise
# using Plots
using Test, PseudoArcLengthContinuation, KrylovKit, LinearAlgebra, Setfield
const PALC = PseudoArcLengthContinuation

# test for constrained problems

####################################################################################################
# we test a dummy problem with constraint

N = (x; a = 0.5, b = 0.01) -> 1 + (x + a*x^2)/(1 + b*x^2)

function F_chan(x, α, β = 0.01)
	f = similar(x)
	n = length(x)
	f[1] = x[1] - β
	f[n] = x[n] - β
	for i=2:n-1
		f[i] = (x[i-1] - 2 * x[i] + x[i+1]) * (n-1)^2 + α * N(x[i], b = β)
	end
	return f
end

n = 101
a = 3.3
ig = [(i-1)*(n-i)/n^2+0.1 for i=1:n]
optnewton = NewtonPar(verbose = true)
sol, hist, flag = newton(x ->  F_chan(x, 3.3), ig, optnewton)

_tau = BorderedArray(rand(n), 1.0)
g(x, p, tau = _tau) = dot(x, tau.u) + (p - 3.) * tau.p
pb = BorderedProblem((x, p) -> F_chan(x, p), (x, p) -> g(x, p))

# test functional with AbstractVector form
pb(vcat(sol, 3.1))
# test functional with BorderedVector form
pb(BorderedArray(sol, 3.1))

# test of newton functional, with BorderedVector / vector
optnewbd = (@set optnewton.linsolver = MatrixBLS(optnewton.linsolver))
newtonBordered(pb, BorderedArray(sol, 3.1), optnewbd )
newtonBordered(pb, vcat(sol, 3.1), optnewbd )
optnewbd = (@set optnewton.linsolver = BorderingBLS(optnewton.linsolver))
newtonBordered(pb, BorderedArray(sol, 3.1),  optnewbd)
newtonBordered(pb, vcat(sol, 3.1), optnewbd )

# test of _newtonPALC
Jchan = (x0, p0) -> PALC.finiteDifferences(x -> F_chan(x, p0), x0)
PALC._newtonPALC( (x,p) -> F_chan(x, p), Jchan, BorderedArray(sol, 3.1), _tau, BorderedArray(sol, 3.1), ContinuationPar(newtonOptions = optnewbd), PALC.DotTheta())

# test of newton functional, with BorderedVector
prob = β -> BorderedProblem((x, p) -> F_chan(x, p, β), (x, p) -> g(x, p))

PseudoArcLengthContinuation.continuationBordered(prob, BorderedArray(sol, 3.1), 0.01, ContinuationPar(newtonOptions = optnewbd, maxSteps = 4), verbosity = 0)
####################################################################################################
# problem with 2 constraints
g2 = (x, p) -> [g(x,p[1]), p[2] - 0.01]

dpF = (x0, p0) -> PALC.finiteDifferences(p ->  F_chan(x0, p[1], p[2]), p0)
dpg = (x0, p0) -> PALC.finiteDifferences(p ->  g2(x0, p), p0)

nestedpb = BorderedProblem(F = (x, p) -> F_chan(x, p[1], p[2]), g = (x, p) -> g2(x, p), dpF = dpF, dpg = dpg, npar = 2)
z0nested = BorderedArray(sol, [3.1, 0.015])
nestedpb(z0nested)

z0nested = vcat(sol, [3.1, 0.015])
nestedpb(z0nested)

optnewnested = @set optnewbd.linsolver = MatrixBLS()
newtonBordered(nestedpb, z0nested, optnewnested)


text:
#######################
# parsing AST into IR #
#######################

function is_differentiable(typ::Type)
    typ <: AbstractFloat || typ <: AbstractArray{T} where {T <: AbstractFloat}
end

include("ir.jl")

struct BasicBlockParseError <: Exception
    expr::Any
end

function parse_addr(expr::Expr)
    @assert expr.head == :macrocall && expr.args[1] == Symbol("@addr")
    line::LineNumberNode = expr.args[2]
    rest = expr.args[3:end]
    if length(rest) == 2
        change_expr = nothing
    elseif length(rest) == 3
        change_expr = rest[3]
    else
        throw(BasicBlockParseError(expr))
    end
    call = rest[1]
    if !isa(call, Expr) || call.head != :call
        throw(BasicBlockParseError(call))
    end
    generator_or_dist = Main.eval(call.args[1])
    if isa(generator_or_dist, Distribution) && change_expr != nothing
        error("Cannot pass change values to @addr for distributions")
    end
    args = call.args[2:end]
    addr = rest[2]
    if isa(addr, Expr) && length(addr.args) == 1 && isa(addr.args[1], Symbol)
        addr = addr.args[1]
    elseif isa(addr, QuoteNode)
        addr = addr.value
    end
    (addr::Symbol, generator_or_dist, args, change_expr, line)
end

function parse_change(expr::Expr)
    @assert expr.head == :macrocall && expr.args[1] == Symbol("@change")
    rest = isa(expr.args[2], LineNumberNode) ? expr.args[3:end] : expr.args[2:end]
    if length(rest) != 1
        throw(BasicBlockParseError(expr))
    end
    addr = rest[1]
    if isa(addr, Symbol)
        pass
    elseif isa(addr, Expr) && length(addr.args) == 1 && isa(addr.args[1], Symbol)
        addr = addr.args[1]
    elseif isa(addr, QuoteNode)
        addr = addr.value
    else
        throw(BasicBlockParseError(addr))
    end
    addr::Symbol
end

function parse_lhs(lhs::Expr)
    if lhs.head == :(::)
        name = lhs.args[1]
        typ = Main.eval(lhs.args[2])
        return (name, typ)
    else
        throw(BasicBlockParseError(lhs))
    end
end

function parse_lhs(lhs::Symbol)
    name = lhs
    (name, Any)
end

function parse_param(expr)
    @assert expr.head == :macrocall && expr.args[1] == Symbol("@param")
    rest = isa(expr.args[2], LineNumberNode) ? expr.args[3:end] : expr.args[2:end]
    if length(rest) != 1
        throw(BasicBlockParseError(rest))
    end
    decl = rest[1]
    if isa(decl, Symbol)
        name = decl
        typ = Any
    elseif isa(decl, Expr) && decl.head == :(::)
        name = decl.args[1]
        typ = Main.eval(decl.args[2])
    else
        throw(BasicBlockParseError(rest[1]))
    end
    (name, typ)
end

function generate_ir(args, body, output_ad, args_ad)
    ir = BasicBlockIR(output_ad, args_ad)
    if !isa(body, Expr) || body.head != :block
        throw(BasicBlockParseError(body))
    end
    for arg in args
        if isa(arg, Symbol)
            name = arg
            typ = Any
        elseif isa(arg, Expr) && arg.head == :(::)
            name = arg.args[1]
            typ = Main.eval(arg.args[2])
        else
            throw(BasicBlockParseError(body))
        end
        add_argument!(ir, name, typ)
    end
    local line::LineNumberNode
    for statement in body.args
        if isa(statement, LineNumberNode)
            line = statement
            continue
        end
        if !isa(statement, Expr)
            throw(BasicBlockParseError(statement))
        end
        if statement.head == :macrocall && statement.args[1] == Symbol("@param")
            (name, typ) = parse_param(statement)
            add_param!(ir, name, typ)
        elseif statement.head == :macrocall && statement.args[1] == Symbol("@addr")
            # an @addr statement without a left-hand-side
            (addr, dist_or_gen, args, change_expr, line) = parse_addr(statement)
            if isa(dist_or_gen, Distribution)
                @assert change_expr == nothing
                add_addr!(ir, addr, line, dist_or_gen, args)
            else
                # change_expr may be nothing, indicating nothing is known
                add_addr!(ir, addr, line, dist_or_gen, args, change_expr)
            end
        elseif statement.head == :(=)
            lhs = statement.args[1]
            rhs = statement.args[2]
            (name, typ) = parse_lhs(lhs)
            if rhs.head == :macrocall && rhs.args[1] == Symbol("@addr")
                (addr, dist_or_gen, args, change_expr, line) = parse_addr(rhs)
                if isa(dist_or_gen, Distribution)
                    @assert change_expr == nothing
                    add_addr!(ir, addr, line, dist_or_gen, args, typ, name)
                else
                    # change_expr may be nothing, indicating nothing is known
                    add_addr!(ir, addr, line, dist_or_gen, args, typ, name, change_expr)
                end
            elseif rhs.head == :macrocall && rhs.args[1] == Symbol("@argschange")
                if (length(rhs.args) == 2 && isa(rhs.args[2], LineNumberNode)) || length(rhs.args) == 1
                    add_argschange!(ir, typ, name)
                else
                    dump(rhs)
                    throw(BasicBlockParseError(statement))
                end
            elseif rhs.head == :macrocall && rhs.args[1] == Symbol("@change")
                addr = parse_change(rhs)
                add_change!(ir, addr, typ, name)
            else
                add_julia!(ir, rhs, typ, name, line)
            end
        elseif statement.head == :return
            if length(statement.args) != 1
                throw(BasicBlockParseError(statement))
            end
            expr = statement.args[1]
            if isa(expr, Symbol)
                set_return!(ir, expr)
            elseif isa(expr, Expr) && expr.head == :(::)
                expr_value = expr.args[1]
                typ = Main.eval(expr.args[2])
                set_return!(ir, expr_value, typ)
            else
                throw(BasicBlockParseError(statement))
            end
        elseif statement.head == :macrocall && statement.args[1] == Symbol("@retchange")
            if length(statement.args) != 2
                throw(BasicBlockParseError(statement))
            end
            expr = statement.args[2]
            if isa(expr, Symbol)
                set_retchange!(ir, expr)
            elseif isa(expr, Expr) && expr.head == :(::)
                expr_value = expr.args[1]
                typ = Main.eval(expr.args[2])
                set_retchange!(ir, expr_value, typ)
            else
                throw(BasicBlockParseError(statement))
            end
        else
            throw(BasicBlockParseError(statement))
        end
    end
    finish!(ir)
    ir
end

function compiled_gen_parse_inner(ast)
    @assert isa(ast, Expr) && ast.head == :macrocall && ast.args[1] == Symbol("@gen")
    ast = isa(ast.args[2], LineNumberNode) ? ast.args[3] : ast.args[2]
    if ast.head != :function
        error("syntax error in @compiled at $(ast) in $(ast.head)")
    end
    if length(ast.args) != 2
        error("syntax error in @compile at $(ast) in $(ast.args)")
    end
    signature = ast.args[1]
    body = ast.args[2]
    if signature.head != :call
        error("syntax error in @compiled at $(ast) in $(signature)")
    end
    name = signature.args[1]
    args = signature.args[2:end]
    has_argument_grads = (map(marked_for_ad, args)...,)
    args = map(strip_marked_for_ad, args)
    (name, args, body, has_argument_grads)
end

function compiled_gen_parse(ast)
    if ast.head == :macrocall && ast.args[1] == Symbol("@gen")
        (name, args, body, args_ad) = compiled_gen_parse_inner(ast)
        return (name, args, body, false, args_ad)
    elseif ast.head == :macrocall && ast.args[1] == Symbol("@ad")
        sub_ast = isa(ast.args[2], LineNumberNode) ? ast.args[3] : ast.args[2]
        if (isa(sub_ast, Expr)
            && sub_ast.head == :macrocall
            && sub_ast.args[1] == Symbol("@gen"))
            (name, args, body, args_ad) = compiled_gen_parse_inner(sub_ast)
            return (name, args, body, true, args_ad)
        end
    end
    error("syntax error in @compiled, expected @compiled @gen .. or @compiled @ad @gen ..")
end

###########################
# basic block trace types #
###########################

# the trace is a mutable struct with fields:
# - a field for each value node (prefixed with a gensym)
# - a field for each addr statement (either a subtrace nor a value)
# - note: there is redundancy between the value nodes and the distribution addr fields

const is_empty_field = gensym("is_empty")
const call_record_field = gensym("call_record")
const value_node_prefix = gensym("value")

function value_field(name::Symbol)
    Symbol("$(value_node_prefix)_$name")
end

function value_field(node::ValueNode)
    value_field(node.name)
end

function value_trace_ref(trace, node::ValueNode)
    fieldname = value_field(node)
    Expr(:(.), trace, QuoteNode(fieldname))
end

struct BasicBlockAssignment{T} <: Assignment
    trace::T
end

static_has_leaf_node(::BasicBlockAssignment, addr) = false
static_has_internal_node(::BasicBlockAssignment, addr) = false

get_address_schema(::Type{BasicBlockAssignment{T}}) where {T} = get_address_schema(T)
Base.isempty(trie::BasicBlockAssignment) = getproperty(trie.trace, is_empty_field)
get_leaf_node(trie::BasicBlockAssignment, key::Symbol) = static_get_leaf_node(trie, Val(key))
get_internal_node(trie::BasicBlockAssignment, key::Symbol) = static_get_internal_node(trie, Val(key))
has_leaf_node(trie::BasicBlockAssignment, key::Symbol) = static_has_leaf_node(trie, Val(key))
has_internal_node(trie::BasicBlockAssignment, key::Symbol) = static_has_internal_node(trie, Val(key))
has_leaf_node(trie::BasicBlockAssignment, addr::Pair) = _has_leaf_node(trie, addr)
get_leaf_node(trie::BasicBlockAssignment, addr::Pair) = _get_leaf_node(trie, addr)
has_internal_node(trie::BasicBlockAssignment, addr::Pair) = _has_internal_node(trie, addr)
get_internal_node(trie::BasicBlockAssignment, addr::Pair) = _get_internal_node(trie, addr)

function make_assignment_methods(trace_type, addr_dist_nodes, addr_gen_nodes)
    methods = Expr[]

    ## get_leaf_nodes ##

    leaf_addrs = map((node) -> node.address, addr_dist_nodes)
    push!(methods, quote
        function Gen.get_leaf_nodes(trie::Gen.BasicBlockAssignment{$trace_type})
            $(Expr(:tuple,
                [quote ($(QuoteNode(addr)), trie.trace.$addr) end for addr in leaf_addrs]...))
        end
    end)

    ## get_internal_nodes ##

    internal_addrs = map((node) -> node.address, addr_gen_nodes)
    push!(methods, quote
        function Gen.get_internal_nodes(trie::Gen.BasicBlockAssignment{$trace_type})
            $(Expr(:tuple,
                [quote ($(QuoteNode(addr)), get_assignment(trie.trace.$addr)) end for addr in internal_addrs]...))
        end
    end)

    for node::AddrDistNode in addr_dist_nodes
        addr = node.address

        ## static_get_leaf_node ##

        push!(methods, quote
            function Gen.static_get_leaf_node(trie::Gen.BasicBlockAssignment{$trace_type},
                                           ::Val{$(QuoteNode(addr))})
                trie.trace.$addr
            end
        end)

        ## static_has_leaf_node ##

        push!(methods, quote
            function Gen.static_has_leaf_node(trie::Gen.BasicBlockAssignment{$trace_type},
                                           ::Val{$(QuoteNode(addr))})
                true
            end
        end)

    end

    for node::AddrGeneratorNode in addr_gen_nodes
        addr = node.address

        ## static_has_internal_node ##

        push!(methods, quote
            function Gen.static_has_internal_node(trie::Gen.BasicBlockAssignment{$trace_type},
                                                  ::Val{$(QuoteNode(addr))})
                true
            end
        end)

        ## static_get_internal_node ##

        push!(methods, quote
            function Gen.static_get_internal_node(trie::Gen.BasicBlockAssignment{$trace_type},
                                                  ::Val{$(QuoteNode(addr))})
                get_assignment(trie.trace.$addr)
            end
        end)
    end

    methods
end

function generate_trace_type(ir::BasicBlockIR, name)
    trace_type_name = gensym("BasicBlockTrace_$name")
    fields = Expr[]
    for (name, node) in ir.value_nodes
        # NOTE: for now record the incremental computations
        # (ir.incremental_nodes) in the trace, but these can be removed for
        # performance optimization
        typ = get_type(node)
        push!(fields, Expr(:(::), value_field(node), QuoteNode(typ)))
    end
    for (addr, node) in ir.addr_dist_nodes
        typ::Type = get_return_type(node.dist)
        push!(fields, Expr(:(::), node.address, QuoteNode(typ)))
    end
    for (addr, node) in ir.addr_gen_nodes
        typ::Type = get_trace_type(node.gen)
        push!(fields, Expr(:(::), node.address, QuoteNode(typ)))
    end
    addresses = union(keys(ir.addr_dist_nodes), keys(ir.addr_gen_nodes))
    assignment_methods = make_assignment_methods(
        trace_type_name, values(ir.addr_dist_nodes), values(ir.addr_gen_nodes))
    retval_type = QuoteNode(ir.output_node ===  nothing ? Nothing : ir.output_node.typ)
    defn = esc(quote

        # specialized trace implementation
        mutable struct $trace_type_name
            $is_empty_field::Bool
            $call_record_field::CallRecord{$retval_type}
            $(Expr(:block, fields...))
            $trace_type_name() = new()
        end

        function Base.copy(other::$trace_type_name)
            trace = $trace_type_name()
            $(Expr(:block, [
                let fieldname = field.args[1]
                    quote trace.$fieldname = other.$fieldname end
                end for field in fields]...))
            trace
        end

        Gen.get_call_record(trace::$trace_type_name) = trace.$call_record_field
        Gen.has_choices(trace::$trace_type_name) = !trace.$is_empty_field

        # assignment view of the trace
        Gen.get_assignment(trace::$trace_type_name) = Gen.BasicBlockAssignment(trace)

        function Gen.get_address_schema(::Type{$trace_type_name})
            Gen.StaticAddressSchema(
                Set{Symbol}([$([QuoteNode(addr) for addr in keys(ir.addr_dist_nodes)]...)]),
                Set{Symbol}([$([QuoteNode(addr) for addr in keys(ir.addr_gen_nodes)]...)]))
        end
        $(Expr(:block, assignment_methods...))
    end)
    (defn, trace_type_name)
end

#########################
# basic block generator #
#########################

abstract type BasicGenFunction{T,U} <: Generator{T,U} end

# a method on the generator type that is executed during expansion of
# generator API generated functions
function get_ir end
function get_grad_fn end

function generate_generator_type(ir::BasicBlockIR, trace_type::Symbol, name::Symbol, node_to_gradient)
    generator_type = gensym("BasicBlockGenerator_$name")
    retval_type = ir.output_node === nothing ? :Nothing : ir.output_node.typ
    defn = esc(quote
        struct $generator_type <: Gen.BasicGenFunction{$retval_type, $trace_type}
            params_grad::Dict{Symbol,Any}
            params::Dict{Symbol,Any}
        end
        $generator_type() = $generator_type(Dict{Symbol,Any}(), Dict{Symbol,Any}())

        (gen::$generator_type)(args...) = get_call_record(simulate(gen, args)).retval
        Gen.get_ir(::Type{$generator_type}) = $(QuoteNode(ir))
        #Gen.render_graph(::$generator_type, fname) = Gen.render_graph(Gen.get_ir($generator_type), fname)
        Gen.get_trace_type(::Type{$generator_type}) = $trace_type
        function Gen.get_static_argument_types(::$generator_type)
            [node.typ for node in Gen.get_ir($generator_type).arg_nodes]
        end
        Gen.accepts_output_grad(::$generator_type) = $(QuoteNode(ir.output_ad))
        Gen.has_argument_grads(::$generator_type) = $(QuoteNode(ir.args_ad))
        Gen.get_grad_fn(::Type{$generator_type}, node::Gen.JuliaNode) = $(QuoteNode(node_to_gradient))[node]
    end)
    (defn, generator_type)
end

# TODO refactor and simplify:
function generate_gradient_fn(node::JuliaNode, gradient_fn::Symbol)
    if isa(node.expr_or_value, Expr) || isa(node.expr_or_value, Symbol)
        input_nodes = node.input_nodes
        inputs_do_ad = map((in_node) -> is_differentiable(get_type(in_node)), node.input_nodes)
        untracked_inputs = [gensym("untracked_$(in_node.name)") for in_node in node.input_nodes]
        maybe_tracked_inputs = [in_node.name for in_node in node.input_nodes]
        track_stmts = Expr[]
        grad_exprs = Expr[]
        grad_exprs_noop = Expr[]
        tape = gensym("tape")
        for (untracked, maybe_tracked, do_ad) in zip(untracked_inputs, maybe_tracked_inputs, inputs_do_ad)
            if do_ad
                push!(track_stmts, quote $maybe_tracked = ReverseDiff.track($untracked, $tape) end)
                push!(grad_exprs, quote ReverseDiff.deriv($maybe_tracked) end)
                push!(grad_exprs_noop, quote zero($untracked) end)
            else
                push!(track_stmts, quote $maybe_tracked = $untracked end)
                push!(grad_exprs, quote nothing end)
                push!(grad_exprs_noop, quote nothing end)
            end
        end
        output_grad = gensym("output_grad")
        given_output_value = gensym("given_output_value")
        output_value_maybe_tracked = gensym("output_value_maybe_tracked")
        err_msg = QuoteNode("julia expression was not differentiable: $(node.expr_or_value)")
        quote
            function $gradient_fn($output_grad, $given_output_value, $(untracked_inputs...))
                $tape = ReverseDiff.InstructionTape()
                $(track_stmts...)
                $output_value_maybe_tracked = $(node.expr_or_value)
                @assert isapprox(ReverseDiff.value($output_value_maybe_tracked), $given_output_value)
                if $output_grad !== nothing
                    if ReverseDiff.istracked($output_value_maybe_tracked)
                        ReverseDiff.deriv!($output_value_maybe_tracked, $output_grad)
                        ReverseDiff.reverse_pass!($tape)
                        return ($(grad_exprs...),)
                    else
                        # the output value was not tracked

                        # this could indicate that the expresssion was not
                        # differentiable, which should probably be an error

                        # or, it could indicate that the expression was a constant

                        # TODO revisit

                        return ($(grad_exprs_noop...),)
                    end
                else
                    # output_grad is nothing (i.e. not a floating point value)
                    return ($(grad_exprs_noop...),)
                end
            end
        end
    else
        # it is a constant value
        @assert length(node.input_nodes) == 0
        quote
            $gradient_fn(output_grad, given_output_value) = ()
        end
    end
end

function generate_gradient_functions(ir::BasicBlockIR)
    gradient_function_defns = Expr[]
    node_to_gradient = Dict{JuliaNode,Symbol}()
    for node::JuliaNode in filter((node) -> isa(node, JuliaNode), ir.all_nodes)
        gradient_fn = gensym("julia_grad_$(node.output.name)")
        gradient_fn_defn = esc(generate_gradient_fn(node, gradient_fn))
        push!(gradient_function_defns, gradient_fn_defn)
        node_to_gradient[node] = gradient_fn
    end
    (gradient_function_defns, node_to_gradient)
end

macro compiled(ast)

    # parse the AST
    (name, args, body, output_ad, args_ad) = compiled_gen_parse(ast)

    # geneate intermediate data-flow representation
    ir = generate_ir(args, body, output_ad, args_ad)

    # generate gradient functions
    (gradient_function_defns, node_to_gradient) = generate_gradient_functions(ir)

    # generate trace type definition
    (trace_type_defn, trace_type) = generate_trace_type(ir, name)

    # generate generator type definition
    (generator_type_defn, generator_type) = generate_generator_type(
        ir, trace_type, name, node_to_gradient)

    Expr(:block,
        trace_type_defn,
        generator_type_defn,
        gradient_function_defns...,
        quote $(esc(name)) = $(esc(generator_type))() end)
end

#####################
# static parameters #
#####################

# V1: just use a dictionary
# V2: create specialized fields.

# note that parameters will be cached (as specialized fields) in the trace;
# user will need to use assess() after changing the parameters to get a trace
# that has the new values of the parameters, before doing e.g. backprop()

# for V1, just during simulate and assess, the parameters will be read from
# dictionaries

function set_param!(gf::BasicGenFunction, name::Symbol, value)
    gf.params[name] = value
end

function get_param(gf::BasicGenFunction, name::Symbol)
    gf.params[name]
end

function get_param_grad(gf::BasicGenFunction, name::Symbol)
    gf.params_grad[name]
end

function zero_param_grad!(gf::BasicGenFunction, name::Symbol)
    gf.params_grad[name] = zero(gf.params[name])
end

function init_param!(gf::BasicGenFunction, name::Symbol, value)
    set_param!(gf, name, value)
    zero_param_grad!(gf, name)
end

######################
# change propagation #
######################

"""
Example: MaskedArgChange{Tuple{Val{:true},Val{:false}},Something}(something)
"""
struct MaskedArgChange{T <: Tuple,U}
    info::U
end

# TODO make the type parameter U part of the BasicGenFunction type parameter?
get_change_type(::BasicGenFunction) = MaskedArgChange

function mask(bits...)
    parameters = map((bit) -> Val{bit}, bits)
    MaskedArgChange{Tuple{parameters...},Nothing}(nothing)
end

export MaskedArgChange, mask

include("simulate.jl")
include("assess.jl")
include("generate.jl")
include("update.jl")
include("backprop.jl")

export @compiled

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
#Calculate spin-spin correlation <S_1 \cdot S_2> for ground state wavefunction
function Spincorr_GS(WRBM::Array{Complex{Float64},1},Nsample::Int64,i_::Int64,j_::Int64)

    #Initialize network parameters
    RBM_Par = init_RBM_par(alpha,nspins)

    #Divide number of samples among CPU workers
    nsample = partition_sample(Nsample)

    #Store number of workers for later use
    NW = Int64[]; push!(NW, NumberWorkers())

    #Update neural network parameters at the current time t
    for j in procs()
        fetch(@spawnat j WToPar(RBM_Par,WRBM))
    end

    #Each worker evaluate spin-spin correlations and sum the outcome
    SP = @distributed (+) for y = 1:NW[1]
        sp_x(RBM_Par,y,nsample,i_,j_)
    end

    spincorr = SP./NW[1]

    return spincorr

end

#Calculate spin-spin correlation <S_1 \cdot S_2> for dynamic wavefunction
function Spincorr_DYN(WRBM_t::Array{Complex{Float64},2},Nsample::Int64,i_::Int64,j_::Int64)

    #Initialize network parameters
    RBM_Par = init_RBM_par(alpha,nspins)

    #Get time interval length from size of WRBM_t
    sizeWRBMt = length(WRBM_t[:,1])

    #Divide number of samples among CPU workers
    nsample = partition_sample(Nsample)

    #Store number of workers for later use
    NW = Int64[]; push!(NW, NumberWorkers())

    #Initialize array storing spin-spin correlations
    spincorr_t = zeros(sizeWRBMt)

    save_ind = 1
    for i=1:sizeWRBMt

        #Update neural network parameters at the current time t
        for j in procs()
            fetch(@spawnat j WToPar(RBM_Par,WRBM_t[i,:]))
        end

        #Each worker evaluate spin-spin correlations and sum the outcome
        SP = @distributed (+) for y = 1:NW[1]
            sp_x(RBM_Par,y,nsample,i_,j_)
        end

        spincorr_t[i] = SP./NW[1]

        #Every 10 steps correlations are saved for back-up
        save_ind += 1
        if save_ind == 10
            writedlm("Spincorr_$(nspins)_$(alpha).jl",spincorr_t)
            save_ind = 1
        end

    end

    return spincorr_t

end

#Evaluate <S_1 \cdot S_2> with Monte Carlo sampling
function sp_x(RBMPar::RBM_par,y::Int64,_nsample::Int64,i_::Int64,j_::Int64)

    #generate rabdom seed for the Monte Carlo chain depending on the worker y
    Random.seed!(100*y)

    #Set number of spin flips for the Monte Carlo chain. Default nflips=2
    nflips=SpinFlips()
    flips = Array{Int64}(undef,nflips)

    #Initialize a random state with zero magnetization if mag0=true
    state = GenRandomState(mag0)

    #initializing look-up tables in the wave-function
    Lt = InitLt(RBMPar,state)

    #thermalization
    @inbounds for n = 1:_nsample*0.1
      @inbounds for i = 1:nspins
                 StateUpdate(RBMPar,Lt,state,flips,nflips)
      end
    end

    spinc = Float64(0)
    #sequence of sweeps
    @inbounds for n = 1:_nsample
        @inbounds for i = 1:nspins
                    StateUpdate(RBMPar,Lt,state,flips,nflips)
                end
                spinc += real(sp_bond(RBMPar,Lt,state,[i_;j_]))
        end

        return spinc/_nsample

end

max_stars_repo_name: GiamFbn/MyULTRAFAST

max_issues_repo_name: GiamFbn/MyULTRAFAST

max_forks_repo_name: GiamFbn/MyULTRAFAST


text:
include("pointCompare.jl")
include("offset.jl")
include("scatterPatch.jl")

#=
Match a patch.
A probe of a search.
The patch was built from the target.
Compare to patch around given point in corpus.

Return ProbeResult

This understands:
- iteration over points of a patch
- short circuit as soon as discover patch cannot be better

Formerly known as computeBestFit()
=#
function patchMatch(
    targetImage::MaskedImage{ValueType, DimensionCount},
    corpusImage::MaskedImage{ValueType, DimensionCount},
    # corpusPatchCenterPoint::Int64,    # framed point, using linear indexing
    corpusPatchCenterPoint::CartesianIndex{DimensionCount},
    synthPatch::ScatterPatch{ValueType, DimensionCount},
    patchDiffToBeat::Float32   # For this target patch, but only while the target is unchanged
    )::ProbeResult where {ValueType, DimensionCount}
    #=
    Iterate over patch of corpus surrounding corpusPoint,
    summing point difference from corresponding patch in source/target.

    A patch is a vector of offsets.
    Offset (can be negative) is not a subtype of Point (only positive)
    Adding an Offset to a Point yields a Point that can be out of range in its frame,
    so the comparison allows that, and compensates.

    As long as patches are not huge,
    no need to worry about loss of precision.
    We are summing many numbers each ~1.0
    =#
    differenceSum::Float32 = zero(Float32)

    # Framed point to mutate in the for loop
    #=
    patchPointCorpus = PointInMaskedImage(
        corpusPatchCenterPoint.point,          # point mutates in the loop
        corpusPatchCenterPoint.maskedImage )   # frame is constant in the loop
    =#
    # patchPointCorpus::CartesianIndex{N}   # N not defined

    #=
    Iterate over points of the patch.

    !!! The patch may be empty:
    when the context is empty, on the first point of the first pass.
    =#
    i = 0
    # @debug "Patch match count of neighbors" length(synthPatch.neighbors)
    for neighbor in synthPatch.neighbors

        #=
        Skip the first neighbor.
        As in the original code.
        TODO omit the first neighbor from synthPatch
        =#
        i += 1
        if i == 1
            continue
        end

        # CartesianIndex arithmetic to get corresponding point in corpus
        patchPointCorpus = offsetPatchPoint(corpusImage.image, corpusPatchCenterPoint, neighbor.offset)

        # is wild, might be out of bounds

        # Difference of value of patchPoint (from the target) to value of corresponding point in corpus
        # TEMP test specialized calls
        differenceSum += comparePointsTargetWithCorpus(
            neighbor, # Neighbor caches values of a framed point in target
            corpusImage, patchPointCorpus   # framed point
            )

        if differenceSum >= patchDiffToBeat
            # patch is already worse than previous best
            # @debug "Short circuit patch match"
            break
        end
    end

    #=
    # assert differenceSum >= patchDiffToBeat
    or is < and is computed over the entire patch
    =#

    betterment = classifyBetterment(differenceSum, patchDiffToBeat)

    #=
    we don't have access to the distinguished patch point in target
    corpusPatchCenterPoint is the patch distinguished point in corpus

    synthPatch.neighbors[1].targetPoint
    =#
    @debug "Patch match" corpusPatchCenterPoint betterment differenceSum

    return ProbeResult(betterment, differenceSum)
end

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
# REPL tests
using TestHelpers
import Base: REPL, LineEdit

function fake_repl()
    # Use pipes so we can easily do blocking reads
    # In the future if we want we can add a test that the right object
    # gets displayed by intercepting the display
    stdin_read,stdin_write = (Base.Pipe(C_NULL), Base.Pipe(C_NULL))
    stdout_read,stdout_write = (Base.Pipe(C_NULL), Base.Pipe(C_NULL))
    stderr_read,stderr_write = (Base.Pipe(C_NULL), Base.Pipe(C_NULL))
    Base.link_pipe(stdin_read,true,stdin_write,true)
    Base.link_pipe(stdout_read,true,stdout_write,true)
    Base.link_pipe(stderr_read,true,stderr_write,true)

    repl = Base.REPL.LineEditREPL(TestHelpers.FakeTerminal(stdin_read, stdout_write, stderr_write))
    stdin_write, stdout_read, stdout_read, repl
end

# Writing ^C to the repl will cause sigint, so let's not die on that
ccall(:jl_exit_on_sigint, Void, (Cint,), 0)
# These are integration tests. If you want to unit test test e.g. completion, or
# exact LineEdit behavior, put them in the appropriate test files.
# Furthermore since we are emulating an entire terminal, there may be control characters
# in the mix. If verification needs to be done, keep it to the bare minimum. Basically
# this should make sure nothing crashes without depending on how exactly the control
# characters are being used.
begin
stdin_write, stdout_read, stdout_read, repl = fake_repl()

repl.specialdisplay = Base.REPL.REPLDisplay(repl)
repl.no_history_file = true

repltask = @async begin
    Base.REPL.run_repl(repl)
end

sendrepl(cmd) = write(stdin_write,"inc || wait(b); r = $cmd; notify(c); r\r")

inc = false
b = Condition()
c = Condition()
sendrepl("\"Hello REPL\"")
inc=true
begin
    notify(b)
    wait(c)
end
# Latex completions
write(stdin_write, "\x32\\alpha\t")
readuntil(stdout_read, "α")
# Bracketed paste in search mode
write(stdin_write, "\e[200~paste here ;)\e[201~")
# Abort search (^C)
write(stdin_write, '\x03')
# Test basic completion in main mode
write(stdin_write, "Base.REP\t")
readuntil(stdout_read, "Base.REPL")
write(stdin_write, '\x03')
write(stdin_write, "\\alpha\t")
readuntil(stdout_read,"α")
write(stdin_write, '\x03')
# Close REPL ^D
write(stdin_write, '\x04')
wait(repltask)
end

function buffercontents(buf::IOBuffer)
    p = position(buf)
    seek(buf,0)
    c = readall(buf)
    seek(buf,p)
    c
end

# Test various history related issues
begin
    stdin_write, stdout_read, stdout_read, repl = fake_repl()
    # In the future if we want we can add a test that the right object
    # gets displayed by intercepting the display
    repl.specialdisplay = Base.REPL.REPLDisplay(repl)

    interface = REPL.setup_interface(repl)
    repl_mode = interface.modes[1]
    shell_mode = interface.modes[2]
    help_mode = interface.modes[3]
    histp = interface.modes[4]

    hp = REPL.REPLHistoryProvider((Symbol=>Any)[:julia => repl_mode,
                                           :shell => shell_mode,
                                           :help  => help_mode])
    fakehistory =
    """
    # time: 2014-06-30 17:32:49 EDT
    # mode: julia
    \tshell
    # time: 2014-06-30 17:32:59 EDT
    # mode: shell
    \tll
    # time: 2014-06-30 17:32:49 EDT
    # mode: julia
    \t1 + 1
    # time: 2014-06-30 17:35:39 EDT
    # mode: foobar
    \tbarfoo
    # time: 2014-06-30 18:44:29 EDT
    # mode: shell
    \tls
    # time: 2014-06-30 19:44:29 EDT
    # mode: foobar
    \tls
    # time: 2014-06-30 20:44:29 EDT
    # mode: julia
    \t2 + 2"""

    REPL.hist_from_file(hp, IOBuffer(fakehistory))
    REPL.history_reset_state(hp)

    histp.hp = repl_mode.hist = shell_mode.hist = help_mode.hist = hp

    # Some manual setup
    s = LineEdit.init_state(repl.t, interface)

    # Test that navigating history skips invalid modes
    # (in both directions)
    LineEdit.history_prev(s, hp)
    @test LineEdit.mode(s) == repl_mode
    @test buffercontents(LineEdit.buffer(s)) == "2 + 2"
    LineEdit.history_prev(s, hp)
    @test LineEdit.mode(s) == shell_mode
    @test buffercontents(LineEdit.buffer(s)) == "ls"
    LineEdit.history_prev(s, hp)
    @test LineEdit.mode(s) == repl_mode
    @test buffercontents(LineEdit.buffer(s)) == "1 + 1"
    LineEdit.history_next(s, hp)
    @test LineEdit.mode(s) == shell_mode
    @test buffercontents(LineEdit.buffer(s)) == "ls"
    LineEdit.history_next(s, hp)
    @test LineEdit.mode(s) == repl_mode
    @test buffercontents(LineEdit.buffer(s)) == "2 + 2"
    LineEdit.history_next(s, hp)

    # Test that searching backwards puts you into the correct mode and
    # skips invalid modes.
    LineEdit.enter_search(s, histp, true)
    ss = LineEdit.state(s, histp)
    write(ss.query_buffer, "l")
    LineEdit.update_display_buffer(ss, ss)
    LineEdit.accept_result(s, histp)
    @test LineEdit.mode(s) == shell_mode
    @test buffercontents(LineEdit.buffer(s)) == "ls"

    # Test that searching for `ll` actually matches `ll` after
    # both letters are types rather than jumping to `shell`
    LineEdit.history_prev(s, hp)
    LineEdit.enter_search(s, histp, true)
    write(ss.query_buffer, "l")
    LineEdit.update_display_buffer(ss, ss)
    write(ss.query_buffer, "l")
    LineEdit.update_display_buffer(ss, ss)
    LineEdit.accept_result(s, histp)
    @test LineEdit.mode(s) == shell_mode
    @test buffercontents(LineEdit.buffer(s)) == "ll"
end

ccall(:jl_exit_on_sigint, Void, (Cint,), 1)

let exename=joinpath(JULIA_HOME,(ccall(:jl_is_debugbuild,Cint,())==0?"julia":"julia-debug"))

# Test REPL in dumb mode
@unix_only begin

const O_RDWR = Base.FS.JL_O_RDWR
const O_NOCTTY = Base.FS.JL_O_NOCTTY

fdm = ccall(:posix_openpt,Cint,(Cint,),O_RDWR|O_NOCTTY)
fdm == -1 && error("Failed to open PTY master")
rc = ccall(:grantpt,Cint,(Cint,),fdm)
rc != 0 && error("grantpt failed")
rc = ccall(:unlockpt,Cint,(Cint,),fdm)
rc != 0 && error("unlockpt")

fds = ccall(:open,Cint,(Ptr{Uint8},Cint),ccall(:ptsname,Ptr{Uint8},(Cint,),fdm), O_RDWR|O_NOCTTY)

# slave
slave   = RawFD(fds)
master = Base.TTY(RawFD(fdm); readable = true)

nENV = copy(ENV)
nENV["TERM"] = "dumb"
p = spawn(setenv(`$exename -f --quiet`,nENV),slave,slave,slave)
start_reading(master)
Base.wait_readnb(master,1)
write(master,"1\nquit()\n")

wait(p)

ccall(:close,Cint,(Cint,),fds)
output = readall(master)
@test output == "julia> 1\r\nquit()\r\n1\r\n\r\njulia> "
close(master)

end

# Test stream mode
outs, ins, p = readandwrite(`$exename -f --quiet`)
write(ins,"1\nquit()\n")
@test readall(outs) == "1\n"
end

max_stars_repo_name: Aerlinger/julia

max_issues_repo_name: Aerlinger/julia

max_forks_repo_name: Aerlinger/julia


text:
# libcublas_types.jl
#
# Initially generated with wrap_c from Clang.jl. Modified to remove anonymous
# enums and add cublasContext.
#
# Author: Nick Henderson <nwh@stanford.edu>
# Created: 2014-08-27
# License: MIT
#

# begin enum cublasStatus_t
const cublasStatus_t = UInt32
const CUBLAS_STATUS_SUCCESS = 0
const CUBLAS_STATUS_NOT_INITIALIZED = 1
const CUBLAS_STATUS_ALLOC_FAILED = 3
const CUBLAS_STATUS_INVALID_VALUE = 7
const CUBLAS_STATUS_ARCH_MISMATCH = 8
const CUBLAS_STATUS_MAPPING_ERROR = 11
const CUBLAS_STATUS_EXECUTION_FAILED = 13
const CUBLAS_STATUS_INTERNAL_ERROR = 14
const CUBLAS_STATUS_NOT_SUPPORTED = 15
const CUBLAS_STATUS_LICENSE_ERROR = 16
# end enum cublasStatus_t
# begin enum cublasFillMode_t
const cublasFillMode_t = UInt32
const CUBLAS_FILL_MODE_LOWER = 0
const CUBLAS_FILL_MODE_UPPER = 1
# end enum cublasFillMode_t
# begin enum cublasDiagType_t
const cublasDiagType_t = UInt32
const CUBLAS_DIAG_NON_UNIT = 0
const CUBLAS_DIAG_UNIT = 1
# end enum cublasDiagType_t
# begin enum cublasSideMode_t
const cublasSideMode_t = UInt32
const CUBLAS_SIDE_LEFT = 0
const CUBLAS_SIDE_RIGHT = 1
# end enum cublasSideMode_t
# begin enum cublasOperation_t
const cublasOperation_t = UInt32
const CUBLAS_OP_N = 0
const CUBLAS_OP_T = 1
const CUBLAS_OP_C = 2
# end enum cublasOperation_t
# begin enum cublasPointerMode_t
const cublasPointerMode_t = UInt32
const CUBLAS_POINTER_MODE_HOST = 0
const CUBLAS_POINTER_MODE_DEVICE = 1
# end enum cublasPointerMode_t
# begin enum cublasAtomicsMode_t
const cublasAtomicsMode_t = UInt32
const CUBLAS_ATOMICS_NOT_ALLOWED = 0
const CUBLAS_ATOMICS_ALLOWED = 1
# end enum cublasAtomicsMode_t
const cublasContext = Nothing
const cublasHandle_t = Ptr{cublasContext}
# complex numbers in cuda
const cuComplex = Complex{Float32}
const cuDoubleComplex = Complex{Float64}
# complex types from Base/linalg.jl
const CublasFloat = Union{Float64,Float32,ComplexF64,ComplexF32}
const CublasReal = Union{Float64,Float32}
const CublasComplex = Union{ComplexF64,ComplexF32}
# FP16 (cuda_fp16.h) in cuda
const __half = Float16
struct __half2
    x1::__half
    x2::__half
end

if CUDAdrv.version() >= v"0.7.5"
    # specify which GEMM algorithm to use in cublasGemmEx() (CUDA 7.5+)
    const cublasGemmAlgo_t = Int32
    const CUBLAS_GEMM_DFALT = -1
    const CUBLAS_GEMM_ALGO0 = 0
    const CUBLAS_GEMM_ALGO1 = 1
    const CUBLAS_GEMM_ALGO2 = 2
    const CUBLAS_GEMM_ALGO3 = 3
    const CUBLAS_GEMM_ALGO4 = 4
    const CUBLAS_GEMM_ALGO5 = 5
    const CUBLAS_GEMM_ALGO6 = 6
    const CUBLAS_GEMM_ALGO7 = 7
    # specify which DataType to use with cublas<t>gemmEx() and cublasGemmEx() (CUDA 7.5+) functions
    const cudaDataType_t = UInt32
    const CUDA_R_16F = UInt32(2)
    const CUDA_C_16F = UInt32(6)
    const CUDA_R_32F = UInt32(0)
    const CUDA_C_32F = UInt32(4)
    const CUDA_R_64F = UInt32(1)
    const CUDA_C_64F = UInt32(5)
    const CUDA_R_8I  = UInt32(3)
    const CUDA_C_8I  = UInt32(7)
    const CUDA_R_8U  = UInt32(8)
    const CUDA_C_8U  = UInt32(9)
    const CUDA_R_32I = UInt32(10)
    const CUDA_C_32I = UInt32(11)
    const CUDA_R_32U = UInt32(12)
    const CUDA_C_32U = UInt32(13)
end

@enum CUBLASMathMode::Cint begin
   CUBLAS_DEFAULT_MATH = 0
   CUBLAS_TENSOR_OP_MATH = 1
end


text:
@testset "Comparisons" begin
    @testset "Equality" begin
        function check_seq_kmer_equality(len)
            a = DNAMer(random_dna_kmer(len))
            b = LongDNASeq(a)
            return a == b && b == a
        end

        function check_seq_bigkmer_equality(len)
            a = BigDNAMer(random_dna_kmer(len))
            b = LongDNASeq(a)
            return a == b && b == a
        end

        for len in [1, 10, 32, 64]
            if len <= 32
                @test all(Bool[check_seq_kmer_equality(len) for _ in 1:reps])
            end
            @test all(Bool[check_seq_bigkmer_equality(len) for _ in 1:reps])
        end

        # True negatives
        @test DNAMer("ACG") != RNAMer("ACG")
        @test DNAMer("T")   != RNAMer("U")
        @test DNAMer("AC")  != DNAMer("AG")
        @test RNAMer("AC")  != RNAMer("AG")

        @test BigDNAMer("ACG") != BigRNAMer("ACG")
        @test BigDNAMer("T")   != BigRNAMer("U")
        @test BigDNAMer("AC")  != BigDNAMer("AG")
        @test BigRNAMer("AC")  != BigRNAMer("AG")

        @test DNAMer("ACG") != rna"ACG"
        @test DNAMer("T")   != rna"U"
        @test DNAMer("AC")  != dna"AG"
        @test RNAMer("AC")  != rna"AG"

        @test BigDNAMer("ACG") != rna"ACG"
        @test BigDNAMer("T")   != rna"U"
        @test BigDNAMer("AC")  != dna"AG"
        @test BigRNAMer("AC")  != rna"AG"

        @test rna"ACG" != DNAMer("ACG")
        @test rna"U"   != DNAMer("T")
        @test dna"AG"  != DNAMer("AC")
        @test rna"AG"  != RNAMer("AC")

        @test rna"ACG" != BigDNAMer("ACG")
        @test rna"U"   != BigDNAMer("T")
        @test dna"AG"  != BigDNAMer("AC")
        @test rna"AG"  != BigRNAMer("AC")
    end

    @testset "Inequality" begin
        for len in [1, 10, 32, 64]
            if len <= 32
                @test  isless(DNAMer{1}(UInt64(0)), DNAMer{1}(UInt64(1)))
                @test !isless(DNAMer{1}(UInt64(0)), DNAMer{1}(UInt64(0)))
                @test !isless(DNAMer{1}(UInt64(1)), DNAMer{1}(UInt64(0)))

                @test  isless(RNAMer{1}(UInt64(0)), RNAMer{1}(UInt64(1)))
                @test !isless(RNAMer{1}(UInt64(0)), RNAMer{1}(UInt64(0)))
                @test !isless(RNAMer{1}(UInt64(1)), RNAMer{1}(UInt64(0)))
            end

            @test  isless(BigDNAMer{1}(UInt128(0)), BigDNAMer{1}(UInt128(1)))
            @test !isless(BigDNAMer{1}(UInt128(0)), BigDNAMer{1}(UInt128(0)))
            @test !isless(BigDNAMer{1}(UInt128(1)), BigDNAMer{1}(UInt128(0)))

            @test  isless(BigRNAMer{1}(UInt128(0)), BigRNAMer{1}(UInt128(1)))
            @test !isless(BigRNAMer{1}(UInt128(0)), BigRNAMer{1}(UInt128(0)))
            @test !isless(BigRNAMer{1}(UInt128(1)), BigRNAMer{1}(UInt128(0)))
        end
    end

    @testset "Hash" begin
        kmers = map(DNAMer, ["AAAA", "AACT", "ACGT", "TGCA"])
        for x in kmers, y in kmers
            @test (x == y) == (hash(x) == hash(y))
        end

        bigkmers = map(BigDNAMer, ["AAAA", "AACT", "ACGT", "TGCA"])
        for x in bigkmers, y in bigkmers
            @test (x == y) == (hash(x) == hash(y))
        end

        kmers = map(RNAMer, ["AAAA", "AACU", "ACGU", "UGCA"])
        for x in kmers, y in kmers
            @test (x == y) == (hash(x) == hash(y))
        end

        bigkmers = map(BigRNAMer, ["AAAA", "AACU", "ACGU", "UGCA"])
        for x in bigkmers, y in bigkmers
            @test (x == y) == (hash(x) == hash(y))
        end
    end
end


text:
rules_text = readlines("day7/input.txt")

struct Bag
    adjective::String
    color::String
end

struct BagRule
    bag::Bag
    holds::Set{Tuple{Int, Bag}}
end

function parse_bag(line::String)
    adj, color = match(r"^(\w+) (\w+) bags", line).captures
    other_bags = Set{Tuple{Int, Bag}}()
    for match in eachmatch(r"(\d) (\w+) (\w+) bag[s]?", line)
        num, adj_inside, color_inside = match.captures
        push!(other_bags, (parse(Int, num), Bag(adj_inside, color_inside)))
    end
    return BagRule(Bag(adj, color), other_bags)
end

bags_rules = parse_bag.(rules_text)

function contains_bags(bag::Bag)
    global bags_rules
    bag_rule = bags_rules[findfirst(b -> b.bag == bag, bags_rules)]
    bags_inside = 0
    for (num, bag) in bag_rule.holds
        bags_inside += num*contains_bags(bag)
        bags_inside += num
    end
    return bags_inside
end

println(contains_bags(Bag("shiny", "gold")))

# that was harder, but we got through it!
# and it was fun! :D


text:
using KernelAbstractions

using KernelAbstractions.Adapt
using KernelAbstractions.CUDA
import ModelParameters.Flatten

"""
    GPU <: Processor

Abstract supertype for GPU processors.
"""
abstract type GPU <: Processor end

export GPU, CuGPU, CPUGPU

"""
    CPUGPU <: GPU

    CPUGPU()

Uses the CUDA GPU code on CPU using KernelAbstractions, to test it.
"""
struct CPUGPU <: GPU end

"""
    CuGPU <: GPU

    CuGPU()
    CuGPU{threads_per_block}()

```julia
ruleset = Ruleset(rule; proc=ThreadedCPU())
# or
output = sim!(output, rule; proc=ThreadedCPU())
```
"""
struct CuGPU{X} <: GPU end
CuGPU() = CuGPU{32}()

kernel_setup(::CPUGPU) where N = KernelAbstractions.CPU(), 1
kernel_setup(::CuGPU{N}) where N = KernelAbstractions.CUDADevice(), N

function Adapt.adapt_structure(to, x::Union{AbstractSimData,GridData,Rule})
    Flatten.modify(A -> adapt(to, A), x, Union{CuArray,Array,AbstractDimArray}, SArray)
end

function Adapt.adapt_structure(to, o::Output)
    frames = map(o.frames) do f
        if f isa NamedTuple
            a = map(g -> adapt(to, g), f)
        else
            adapt(to, f)
        end
    end
    @set o.extent = adapt(to, o.extent)
    @set o.frames = frames
end

function _proc_setup(::CuGPU, obj)
    Flatten.modify(CuArray, obj, Union{Array,BitArray}, Union{SArray,Dict})
end

function _copyto_output!(outgrid, grid, proc::GPU)
    src = adapt(Array, source(grid))
    copyto!(outgrid, CartesianIndices(outgrid), src, CartesianIndices(outgrid))
    return nothing
end
function _copyto_output!(outgrid::CuArray, grid, proc::CuGPU)
    gridindices = CartesianIndices(map(a -> a .+ radius(grid), axes(outgrid)))
    copyto!(outgrid, CartesianIndices(outgrid), parent(source(grid)), gridindices)
    return nothing
end

function maprule!(
    wgrids::Union{<:GridData{Y,X,R},Tuple{<:GridData{Y,X,R},Vararg}},
    simdata, proc::GPU, opt, ruletype::Val{<:Rule}, rule, args...
) where {Y,X,R}
    kernel! = cu_cell_kernel!(kernel_setup(proc)...)
    kernel!(wgrids, simdata, ruletype, rule, args...; ndrange=gridsize(simdata)) |> wait
end
function maprule!(
    wgrids::Union{<:GridData{Y,X,R},Tuple{<:GridData{Y,X,R},Vararg}},
    simdata, proc::GPU, opt, ruletype::Val{<:NeighborhoodRule}, rule, args...
) where {Y,X,R}
    grid = simdata[neighborhoodkey(rule)]
    kernel! = cu_neighborhood_kernel!(kernel_setup(proc)...)
    # n = _indtoblock.(gridsize(simdata), 8) .- 1
    n = gridsize(simdata)
    kernel!(wgrids, simdata, grid, opt, ruletype, rule, args..., ndrange=n) |> wait
    return nothing
end

@kernel function cu_neighborhood_kernel!(
    wgrids, data, grid::GridData{Y,X,R}, opt, ruletype::Val{<:NeighborhoodRule}, rule, args...
) where {Y,X,R}
    I, J = @index(Global, NTuple)
    src = parent(source(grid))
    @inbounds buf = view(src, I:I+2R, J:J+2R)
    bufrule = _setbuffer(rule, buf)
    cell_kernel!(wgrids, data, ruletype, bufrule, args..., I, J)
    nothing
end

@kernel function cu_cell_kernel!(wgrids, simdata, ruletype::Val, rule, rkeys, rgrids, wkeys)
    i, j = @index(Global, NTuple)
    cell_kernel!(wgrids, simdata, ruletype, rule, rkeys, rgrids, wkeys, i, j)
    nothing
end

# @kernel function cu_cell_kernel!(
#     simdata::SimData, grid::GridData{Y,X,1}, rule::NeighborhoodRule,
#     rkeys, rgrids, wkeys, wgrids
# ) where {Y,X}
#     gi, gj = @index(Group, NTuple)
#     i, j = @index(Local,  NTuple)
#     N = @uniform groupsize()[1]
#     M = @uniform groupsize()[2]
#     # +1 to avoid bank conflicts on shared memory
#     tile = @localmem eltype(src) 8, 8
#     for n in 1:4
#         b = ((i - 1) & 1)
#         k = 4b + n
#         l = 2j + b - 1
#         @inbounds tile[k, l] = src[(gi-1) * N + k, (gi-1) * N + l]
#     end
#     @synchronize
#     # Manually calculate global indexes
#     I = (gi-1) * N + i
#     J = (gj-1) * M + j
#     @inbounds buf = view(tile, i:i+2, j:j+2)
#     bufrule = setbuffer(rule, buf)
#     readval = _readgrids(rkeys, rgrids, I, J)
#     writeval = applyrule(simdata, bufrule, readval, (I, J))
#     _writegrids!(wgrids, writeval, I, J)
#     nothing
# end

for (f, op) in atomic_ops
    atomic_f = Symbol(:atomic_, f)
    @eval begin
        @propagate_inbounds function ($f)(::CuGPU, d::WritableGridData{Y,X,R}, x, I...) where {Y,X,R}
            A = parent(dest(d))
            i = Base._to_linear_index(A, (I .+ R)...)
            ($f)(A, x, i)
        end
        @propagate_inbounds function ($f)(A::CUDA.CuDeviceArray, x, i) where {Y,X,R}
            (CUDA.$atomic_f)(pointer(A, i), x)
        end
    end
end


text:
Base.HOME_PROJECT[] = abspath(Base.HOME_PROJECT[]) # JuliaLang/julia/pull/28625

using CLIMA, Documenter

# TODO: Add generated examples back
# include("generate.jl")

GENERATED_BL_EXAMPLES = [
    joinpath("examples", "DGmethods_old", "generated", f)
    for
    f in ("ex_001_periodic_advection.md", "ex_002_solid_body_rotation.md")
]
GENERATED_BL_EXAMPLES = filter!(x -> isfile(x), GENERATED_BL_EXAMPLES)

pages = Any[
    "Home" => "index.md",
    "Common" => Any["MoistThermodynamics" => "Common/MoistThermodynamics.md",],
    "Utilites" => Any[],
    "Atmos" => Any[
        "Atmos/SurfaceFluxes.md",
        "Atmos/EDMFEquations.md",
        "Microphysics" => "Atmos/Microphysics.md",
    ],
    "ODESolvers" => "ODESolvers.md",
    "LinearSolvers" => "LinearSolvers.md",
    "Mesh" => "Mesh.md",
    "Arrays" => "Arrays.md",
    "DGmethods_old" => "DGmethods_old.md",
    "InputOutput.md",
    "Developer docs" => Any[
        "CodingConventions.md",
        "AcceptableUnicode.md",
        "VariableList.md",
    ],
]

if !isempty(GENERATED_BL_EXAMPLES)
    push!(
        pages,
        "Balance Law Examples" =>
            ["BalanceLawOverview.md", GENERATED_BL_EXAMPLES...],
    )
end

makedocs(
    sitename = "CLIMA",
    doctest = false,
    strict = false,
    linkcheck = false,
    checkdocs = :exports,
    # checkdocs = :all,
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        mathengine = MathJax(Dict(
            :TeX => Dict(
                :equationNumbers => Dict(:autoNumber => "AMS"),
                :Macros => Dict(),
            ),
        )),
        # prettyurls = !("local" in ARGS),
        # canonical = "https://climate-machine.github.io/CLIMA/stable/",
    ),
    clean = true,
    modules = [Documenter, CLIMA],
    pages = pages,
)

# make sure there are no *.vtu files left around from the build
p = joinpath(@__DIR__, "build", "examples", "DGmethods_old", "generated")
if ispath(p)
    cd(p) do
        foreach(file -> endswith(file, ".vtu") && rm(file), readdir())
    end
end

deploydocs(repo = "github.com/climate-machine/CLIMA.git", target = "build")

max_stars_repo_name: trontrytel/CLIMA

max_issues_repo_name: trontrytel/CLIMA

max_forks_repo_name: trontrytel/CLIMA


text:
"Holds the tableau of a stochastic explicit Runge-Kutta method."
struct TableauSERK{T} <: AbstractTableauERK{T}
    name::Symbol
    s::Int

    qdrift::CoefficientsRK{T}
    qdiff::CoefficientsRK{T}
    qdiff2::CoefficientsRK{T}

    # Order of the tableau is not included, because unlike in the deterministic
    # setting, it depends on the properties of the noise (e.g., the dimension of
    # the Wiener process and the commutativity properties of the diffusion matrix)
    #
    # Orders stored in qdrift, qdiff and qdiff2 are understood as the classical orders of these methods.

    function TableauSERK{T}(name, qdrift, qdiff, qdiff2) where {T}
        @assert qdrift.s == qdiff.s == qdiff2.s
        @assert qdrift.c[1] == qdiff.c[1] == qdiff2.c[1] == 0.
        @assert istrilstrict(qdrift.a)
        @assert istrilstrict(qdiff.a)
        @assert istrilstrict(qdiff2.a)
        @assert !(qdrift.s==1. && qdrift.a[1,1] ≠ 0.)
        @assert !(qdiff.s==1. && qdiff.a[1,1] ≠ 0.)
        @assert !(qdiff2.s==1. && qdiff2.a[1,1] ≠ 0.)
        new(name, qdrift.s, qdrift, qdiff, qdiff2)
    end
end

function TableauSERK(name, qdrift::CoefficientsRK{T}, qdiff::CoefficientsRK{T}, qdiff2::CoefficientsRK{T}) where {T}
    TableauSERK{T}(name, qdrift, qdiff, qdiff2)
end

function TableauSERK(name, qdrift::CoefficientsRK{T}, qdiff::CoefficientsRK{T}) where {T}
    TableauSERK{T}(name, qdrift, qdiff, CoefficientsRK(T, :NULL, 0, zero(qdrift.a), zero(qdrift.b), zero(qdrift.c)) )
end

function TableauSERK(name::Symbol, order_drift::Int, a_drift::Matrix{T}, b_drift::Vector{T}, c_drift::Vector{T},
                                    order_diff::Int, a_diff::Matrix{T}, b_diff::Vector{T}, c_diff::Vector{T},
                                    order_diff2::Int, a_diff2::Matrix{T}, b_diff2::Vector{T}, c_diff2::Vector{T}) where {T}
    TableauSERK{T}(name, CoefficientsRK(name, order_drift, a_drift, b_drift, c_drift), CoefficientsRK(name, order_diff, a_diff, b_diff, c_diff),
                    CoefficientsRK(name, order_diff2, a_diff2, b_diff2, c_diff2))
end

function TableauSERK(name::Symbol, order_drift::Int, a_drift::Matrix{T}, b_drift::Vector{T}, c_drift::Vector{T},
                                    order_diff::Int, a_diff::Matrix{T}, b_diff::Vector{T}, c_diff::Vector{T}) where {T}
    TableauSERK{T}(name, CoefficientsRK(name, order_drift, a_drift, b_drift, c_drift), CoefficientsRK(name, order_diff, a_diff, b_diff, c_diff),
                    CoefficientsRK(:NULL, 0, zero(a_drift), zero(b_drift), zero(c_drift)))
end

"Stochastic Explicit Runge-Kutta integrator."
struct IntegratorSERK{DT, TT, ET <: SDE{DT,TT}} <: StochasticIntegrator{DT,TT}
    equation::ET
    tableau::TableauSERK{TT}
    Δt::TT

    Δy::Vector{DT}

    Q::Vector{Vector{DT}}     # Q[j][k] - the k-th component of the j-th internal stage
    V::Vector{Vector{DT}}     # V[j][k] - the k-th component of v(Q[j])
    B::Vector{Matrix{DT}}     # B[j]    - the diffusion matrix B(Q[j])

    function IntegratorSERK{DT,TT}(equation::ET, tableau, Δt::TT) where {DT, TT, ET <: SDE{DT,TT}}
        D = equation.d
        M = equation.m
        NS= max(equation.ns,equation.ni)
        S = tableau.s

        # create internal stage vectors
        Q = create_internal_stage_vector(DT, D, S)
        V = create_internal_stage_vector(DT, D, S)
        B = create_internal_stage_vector(DT, D, M, S)

        new{DT,TT,ET}(equation, tableau, Δt, zeros(DT,M), Q, V, B)
    end
end

function IntegratorSERK(equation::SDE{DT,TT}, tableau::TableauSERK{TT}, Δt::TT) where {DT,TT}
    IntegratorSERK{DT,TT}(equation, tableau, Δt)
end

@inline equation(integrator::IntegratorSERK) = integrator.equation
@inline timestep(integrator::IntegratorSERK) = integrator.Δt
@inline tableau(integrator::IntegratorSERK)  = integrator.tableau
@inline nstages(integrator::IntegratorSERK)  = nstages(tableau(integrator))
@inline eachstage(integrator::IntegratorSERK) = 1:nstages(integrator)

"""
Integrate SDE with explicit Runge-Kutta integrator.
  Calculating the n-th time step of the explicit integrator for the sample path m
"""
function integrate_step!(int::IntegratorSERK{DT,TT}, sol::AtomicSolutionSDE{DT,TT}) where {DT,TT}
    local tᵢ::TT
    local ydrift::DT

    # reset cache
    reset!(sol, timestep(int))

    # calculates v(t,tQ) and assigns to the i-th column of V
    int.equation.v(sol.t̅, sol.q̅, int.V[1])
    # calculates B(t,tQ) and assigns to the matrix BQ[1][:,:]
    int.equation.B(sol.t̅, sol.q̅, int.B[1])

    @inbounds for i in 2:int.tableau.s
        for k in eachindex(int.Q[i])
            # contribution from the drift part
            ydrift = 0
            for j = 1:i-1
                ydrift += int.tableau.qdrift.a[i,j] * int.V[j][k]
            end

            # ΔW contribution from the diffusion part
            int.Δy .= 0
            for j = 1:i-1
                for l = 1:equation(int).m
                    int.Δy[l] += int.tableau.qdiff.a[i,j] * int.B[j][k,l]
                end
            end

            int.Q[i][k] = sol.q̅[k] + int.Δt * ydrift + dot(int.Δy,sol.ΔW)

            # ΔZ contribution from the diffusion part
            if int.tableau.qdiff2.name ≠ :NULL
                int.Δy .= 0
                for j = 1:i-1
                    for l = 1:equation(int).m
                        int.Δy[l] += int.tableau.qdiff2.a[i,j] * int.B[j][k,l]
                    end
                end

                int.Q[i][k] += dot(int.Δy,sol.ΔZ)/int.Δt
            end

        end
        tᵢ = sol.t̅ + int.Δt * int.tableau.qdrift.c[i]
        int.equation.v(tᵢ, int.Q[i], int.V[i])
        int.equation.B(tᵢ, int.Q[i], int.B[i])
    end

    # compute final update
    if int.tableau.qdiff2.name == :NULL
        update_solution!(sol.q, int.V, int.B, int.tableau.qdrift.b, int.tableau.qdiff.b, int.Δt, sol.ΔW, int.Δy)
    else
        update_solution!(sol.q, int.V, int.B, int.tableau.qdrift.b, int.tableau.qdiff.b, int.tableau.qdiff2.b, int.Δt, sol.ΔW, sol.ΔZ, int.Δy)
    end
end


text:
###############################################################################
#
# Create LeafBios
#
###############################################################################
"""
    create_leaf_bios(FT, rt_dim::RTDimensions)

Create a [`LeafBios`](@ref) type struct, given
- `FT` Floating number type
- `rt_dim` [`RTDimensions`](@ref) type struct
"""
function create_leaf_bios(FT, rt_dim::RTDimensions)
    @unpack nWL, nWLE, nWLF = rt_dim;

    return LeafBios{FT}(nWL  = nWL ,
                        nWLE = nWLE,
                        nWLF = nWLF)
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
import IntervalArithmetic: (⊆), (≺), contains_zero, isempty, isnan,
    precedes, ⊂, isinterior, isdisjoint

export (⊆), (≺), contains_zero, isempty, isnan, precedes, ⊂, isinterior,
    isdisjoint

#NOTE missing: interior, disjoint
for f in (:(⊆), :(≺), :precedes, :⊂, :isinterior, :isdisjoint)
    @eval $f(a::NumberInterval, b::NumberInterval) = $f(Interval(a), Interval(b))
end

for f in (:contains_zero, :isempty, :isnan)
    @eval $f(a::NumberInterval) = $f(Interval(a))
end

# these fit better in IntervalArithmetic so we don't export them here
issingleton(a::NumberInterval) = iszero(radius(a))


text:
"""
Print the sequence of outbound nodes of a graph starting from the given index.
pystyle implies reversing the final order and substracting 1 from the vertex number.
"""
function prettyp(G::AbstractGraph, strtidx::Int; pystyle=false)
    arrow = G isa LightGraphs.SimpleDiGraph ? "->" : "--"
    v = unroll(G, strtidx)
    println(join(string.(pystyle ? reverse(v .- 1) : v), arrow))
end

"""
Walk the given graph, starting from the given vertex,
returning the list of all outbound vertices that are visited.
"""
function unroll(G::AbstractGraph, idx::Int; visited=Int[])
    idx in visited && return visited
    push!(visited, idx)
    N = outneighbors(G, idx)
    length(N) < 1 && length(visited) == nv(G) ? (return visited) : nothing
    for n in N
        unroll(G, n, visited=visited)
    end
    return reverse(visited)
end

"""
Return the first and last `len` elements of the vector as a string.

Default is to print 3 elements at head and tail.

```julia
julia> v = collect(1:10);
julia> prettyp(v)

"1,2,3...8,9,10"
```

Example w/5 elements visible

```julia
julia> v = collect(1:10);
julia> prettyp(v, 5)

"1,2,3,4,5...6,7,8,9,10"
```

"""
function prettyp(v::AbstractVector, len::Int = 3)
    length(v) < 7 && return join(string.(v), ",")
    head = join(string.(v[1:len]),",")
    tail = join(string.(v[(end-(len-1)):end]),",")
    return join([head,tail],"...")
end


text:
module ImagePartitioning

abstract type PartitionScheme end
struct AllowContraction <: PartitionScheme end

include("partition.jl")

export
	# main functions
    partition_image,
	AllowContraction

end # module


text:
module TestRandomMC

using QuantEcon
using Base.Test
using FactCheck

facts("Testing random_mc.jl") do
    context("Test random_markov_chain") do
        n, k = 5, 3
        mc_dicts = (Dict("P" => random_markov_chain(n).p, "k" => n),
                    Dict("P" => random_markov_chain(n, k).p, "k" => k))
        for d in mc_dicts
            P = d["P"]
            @fact size(P) --> (n, n)
            @fact all(x->(countnz(x)==d["k"]),
                      [P[i, :] for i in 1:size(P)[1]]) --> true
        end
    end

    context("Test random_stochastic_matrix") do
        n, k = 5, 3
        Ps = (random_stochastic_matrix(n), random_stochastic_matrix(n, k))
        for P in Ps
            @fact all(P .>= 0) --> true
            @fact all(x->isapprox(sum(x), 1),
                      [P[i, :] for i in 1:size(P)[1]]) --> true
        end
    end

    context("Test errors properly thrown") do
        # n <= 0
        @fact_throws random_markov_chain(0)

        # k <= 0
        @fact_throws random_markov_chain(2, 0)

        # k > n
        @fact_throws random_markov_chain(2, 3)
    end
end  # facts

end  # module


text:
export iauGmst82
"""
Universal Time to Greenwich mean sidereal time (IAU 1982 model).

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  canonical model.

Given:
   dj1,dj2    double    UT1 Julian Date (see note)

Returned (function value):
               double    Greenwich mean sidereal time (radians)

Notes:

   1. The UT1 date dj1+dj2 is a Julian Date, apportioned in any
   convenient way between the arguments dj1 and dj2.  For example,
   JD(UT1)=2450123.7 could be expressed in any of these ways,
   among others:

            dj1            dj2

      2450123.7          0          (JD method)
         2451545        -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5         0.2         (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 and MJD methods are good compromises
   between resolution and convenience.  The date & time method is
   best matched to the algorithm used:  maximum accuracy (or, at
   least, minimum noise) is delivered when the dj1 argument is for
   0hrs UT1 on the day in question and the dj2 argument lies in the
   range 0 to 1, or vice versa.

   2. The algorithm is based on the IAU 1982 expression.  This is
   always described as giving the GMST at 0 hours UT1.  In fact, it
   gives the difference between the GMST and the UT, the steady
   4-minutes-per-day drawing-ahead of ST with respect to UT.  When
   whole days are ignored, the expression happens to equal the GMST
   at 0 hours UT1 each day.

   3. In this function, the entire UT1 (the sum of the two arguments
   dj1 and dj2) is used directly as the argument for the standard
   formula, the constant term of which is adjusted by 12 hours to
   take account of the noon phasing of Julian Date.  The UT1 is then
   added, but omitting whole days to conserve accuracy.

Called:
   iauAnp       normalize angle into range 0 to 2pi

References:

   Transactions of the International Astronomical Union,
   XVIII B, 67 (1983).

   Aoki et al., Astron.Astrophys., 105, 359-361 (1982).

This revision:  2017 October 12

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# double iauGmst82(double dj1, double dj2)

function iauGmst82(dj1::Real, dj2::Real)
   return ccall((:iauGmst82, libsofa_c), Cdouble,
               (Cdouble, Cdouble),
               convert(Float64, dj1), convert(Float64, dj2))
end


text:
using RAFF
using DelimitedFiles
using Printf

# Set Debug for Logging
using Logging
using Base.CoreLogging

function run_lmlovo(p::Int, initguess=nothing; model_str="logistic")

    global_logger(ConsoleLogger(stdout, Logging.Error))

    n, model, modelstr = RAFF.model_list[model_str]

    open("/tmp/output.txt") do fp

        global N = parse(Int, readline(fp))

        global data = readdlm(fp)

    end

    if initguess == nothing

        initguess = zeros(Float64, n)

    end

    rsol = lmlovo(model, initguess, data[:, 1:2], n, p)

    @printf("Solution found:
            fbest = %f
            p     = %d\n", rsol.f, rsol.p)
    println(rsol.solution)

    return rsol

end


text:
abstract AssetSpread <: FinancialAsset

type CoveredCall <: AssetSpread
    long::LongStock
    short::ShortCall
    risk::Float64

    CoveredCall(long, short) = new(long, short, (long.basis - short.basis))
end

type LongPutVertical <: AssetSpread
    long::LongPut
    short::ShortPut
    risk::Float64

    function LongPutVertical(long, short)
        if typeof(long) == LongPut && typeof(short) == ShortPut
            risk = long.basis - short.basis
            new(long, short, risk)
        elseif typeof(long) == ShortPut && typeof(short) == LongPut
            risk = short.basis - long.basis
            new(short, long, risk) # flip arguments if they come in the wrong order
        else
            error("elements must include one ShortPut and one LongPut") # enforce one short and one long and switch them if order is incorrect
        end
    end
end

type LongCallVertical <: AssetSpread
    long::LongCall
    short::ShortCall
    risk::Float64

    function LongCallVertical(long, short)
        _checkcall(long,short)
        risk = short.basis - long.basis
        new(long, short, risk)
    end
end

type ShortPutVertical <: AssetSpread
    long::LongPut
    short::ShortPut
    risk::Float64

    function ShortPutVertical(long, short)
        _checkput(long,short)
        risk = abs(long.strike - short.strike) - (short.basis - long.basis)
        new(long, short, risk)
    end
end

type ShortCallVertical <: AssetSpread
    long::LongCall
    short::ShortCall
    risk::Float64

    function ShortCallVertical(long, short)
        _checkcall(long,short)
        risk = abs(long.strike - short.strike) - (short.basis - long.basis)
        new(long, short, risk)
    end
end

# pull out invariant code to avoid duplication

function _checkput(long, short)

    if long.ticker.id != short.ticker.id
        error("Spreads with different underlyings not supported") # enforce same underlying
    else
        nothing
    end

    if typeof(long) == ShortPut && typeof(short) == LongPut
        return(temp = short; short = long; long = temp) # flip arguments if they come in the wrong order
    elseif typeof(long) == LongPut && typeof(short) == ShortPut
        nothing
    else
        error("elements must include one ShortPut and one LongPut") # enforce one short and one long and switch them if order is incorrect
    end
end

function _checkcall(long, short)

    if long.ticker.id != short.ticker.id
        error("Spreads with different underlyings not supported") # enforce same underlying
    else
        nothing
    end

    if typeof(long) == ShortCall && typeof(short) == LongCall
        return(short = long ; long  = short) # flip arguments if they come in the wrong order
    elseif typeof(long) == LongCall && typeof(short) == ShortCall
        nothing
    else
        error("elements must include one ShortPut and one LongPut") # enforce one short and one long and switch them if order is incorrect
    end
end

# show
function show(io::IO, s::AssetSpread)
    typ = split(string(typeof(s)), ".")[2] # FinancialAssets.LongPut -> LongPut (e.g.)
    println(io, @sprintf("spread:         %s%s", typ, " Spread"))
    println(io, @sprintf("ticker:         %s", s.long.ticker))
    println(io, @sprintf("risk:           %s", trunc(s.risk,2))) # only show two decimals
end


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

using Markdown
using InteractiveUtils

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

# ╔═╡ f27d484a-0993-11eb-08c9-030d9f88b754
using Plots, Distributions, PlutoUI, DSP, Statistics, CardioModels

# ╔═╡ 60c79388-12b7-11eb-0493-51151fa697ab
md"""
### About the model

A main feature of the model is that besides parasympathetic action through the baroreflex, sympathetic action is added. Sympathetic activity is triggered if the diastolic pressure falls below a reference value and influences heart rate, inotropy, and systemic resistance with some delay. In this way, a 0.1Hz oscillation is induced.

#### The overall model is formulated as follows

##### Sympathetic activation

$SympD_n = 0 \hspace{5mm} \text{if} \hspace{22mm} D_{n} > D_{\text{ref}}$
$SympD_n = D_{\text{ref}} - D_n \hspace{5mm} \text{if} \hspace{5mm} D_n \leq D_{\text{ref}}$
$Symp_n = \text{atan} \left(G_{SA}(z) \cdot SympD_{n}\right)$

##### Baroreflex

$I_n^{\star} = a_1 S_{n-1} + I_0 + w_I + A_{\rho}^I\sin(2\pi f_{\rho}t + \phi) \hspace{5mm} \text{if} \hspace{5mm}I_{n-1} > 700\ ms$

$I_n^{\star} = a_2 S_{n-2} + I_0 + w_I + A_{\rho}^I\sin(2\pi f_{\rho}t + \phi) \hspace{5mm} \text{if} \hspace{5mm}I_{n-1} \leq 700\ ms$

$I_n = \frac{I_n^{\star}}{1 + \beta Symp_n}$

with $\beta$ the sympathetic 'force' on $I$

##### Windkessel

$D_n = S_{n} \exp\left(\frac{-I_n }{R C (1 + \delta Symp_n)}\right)$

with $RC$ the diastolic runoff and $\delta$ the sympathetic 'force' on $D$

##### Starling/Restitution

$P_n = (\gamma I_n + P_0)(1+ \varepsilon Symp_n) + w_P + A_{\rho}^P\sin(2\pi f_{\rho}t + \phi)$
with $\varepsilon$ the sympathetic 'force' on $P$ (inotropy by symathetics and respiratory BP modulation)
\
The next systole is determined by $S_{n+1} = D_n + P_n$
"""

# ╔═╡ 004f4ad1-4160-4e52-8813-70b9c73e3a29
md"# Create a model:"

# ╔═╡ 7efa1b0e-0a5c-11eb-14c9-032c18080c9b
md"""
Diastolic reference value:
\
$(@bind Dref NumberField(30:160, default = 80))
mmHg
"""

# ╔═╡ 25b183da-0998-11eb-0a1f-03534abcdc2e
md"""
Diastolic runoff:
\
$(@bind t_runoff NumberField(1:5000, default = 2000))
ms
"""

# ╔═╡ cdba0c98-120b-11eb-2aa3-d9e84e989773
md"""
Offset:
\
$(@bind offset NumberField(1:100, default = 20))
mmHg
"""

# ╔═╡ ac0c0478-0998-11eb-31bc-e3e4193eb47f
md"""
Baro strength:
\
a1:
$(@bind a1 NumberField(0:0.01:30, default = 9))
mmHg/ms
"""

# ╔═╡ f7a5c9ca-0998-11eb-1a7f-c57908a3cee5
md"""
a2:
$(@bind a2 NumberField(0:0.01:30, default = 9))
mmHg/ms
"""

# ╔═╡ 918ea734-120a-11eb-33e7-4340ad35c9f1
md"""
Sympathetic Force on D:
\
$(@bind b_D NumberField(0:0.001:1, default = 0.25))

"""

# ╔═╡ ee91b962-12c3-11eb-24aa-ff59fbe0e53e
md"""
Sympathetic Force on P:
\
$(@bind b_P NumberField(0:0.001:1, default = 0.25))
"""

# ╔═╡ ed8b4472-12c3-11eb-0d55-47162f88dc88
md"""
Sympathetic Force on I:
\
$(@bind b_I NumberField(0:0.001:1, default = 0.25))
"""

# ╔═╡ eb927d16-120a-11eb-2a23-277da77b03f6
md"""
Starling Effect:
\
$(@bind starling NumberField(0:0.001:1, default = 0.02))
mmHg/ms
"""

# ╔═╡ 78a242de-120b-11eb-1662-a948f6a40cf5
w = [0.000,0.00, 0.015, 0.030, 0.02, 0.010];

# ╔═╡ a10eb0f2-50f1-11eb-0de4-f3e26e2ed897
md"""
Breathing period T:
\
$(@bind T_resp NumberField(0:50:10000, default = 4000))
ms
"""

# ╔═╡ 3d310136-0991-11eb-05ea-4351c72e2260
model = KaremakerModel(a1 = a1, a2 = a2, t_runoff = t_runoff, Dref = Dref, w = w, b_D = b_D, b_I = b_I, b_P = b_P, starling = starling, offset = offset, T_resp = T_resp, A_resp = 60, magn_pulse = 5, A_BP = 10);
#model = KaremakerModel();

# ╔═╡ 2bc64aaf-c07a-4e75-84e7-ce71f743f280
md"## Generate data:"

# ╔═╡ 23cc206d-78a3-4be0-b653-a8c5123d8dd9
md"""
N:
$(@bind N NumberField(1:1000, default = 100))
"""

# ╔═╡ 4935f00e-0991-11eb-2017-210b39b3aa56
S, D, P, I, Symp, ρ = predict(model, N, burnIn = 200);

# ╔═╡ b01c33d6-0997-11eb-222c-d34da1f9aa4b
md"""
Schow spectral density:
$(@bind p_psd CheckBox())
"""

# ╔═╡ e3c9adde-0993-11eb-191c-b501deb4a46f
begin
	if !p_psd
		pi = plot(I, title = "RR [ms]", lab = "")
		pbp = plot(S, title = "BP [mmHg]", lab = "", color = :red)
		plot!(pbp, D, lab = "", color = :green)
		plot(pi, pbp, layout = (2,1), dpi = 500)
	else
		psdi = psdplot(I, title = "RR [ms]", method = "welch", lab = "",  fs = 1000/ (mean(I)))
		psdbp = psdplot(S, title = "BP [mmHg]", lab = "", color = :red, fs = 1000/ (mean(I)))
		psdplot!(psdbp, D, lab = "", color = :green, fs = 1000/ (mean(I)))
		plot(psdi, psdbp, layout = (2,1), size = (700, 500), fs = 1000/ (mean(I)))
		vline!([1/(0.001T_resp)], lab = "f_resp")
	end
end

# ╔═╡ 0f3b19aa-120c-11eb-3c2a-6f5e3e1e42fa
md"""
**Experiment on LF / HF ratio as measure for sympatho-vagal balance:**
\
Diastolic runoff as a measure of overall systemic resistence is varied resulting in a change in sympathetic drive, which is activated only over the diastolic pressure in this model.
"""

# ╔═╡ 4eb77de8-120d-11eb-3d3b-f1414ea32fef
begin
	model2 = KaremakerModel(b_I = 0.5)
	runoffs = range(1250, 2850, length = 100)
	rr,symp_drive, lfhf, sb, db = vagal_balance(55, runoffs, model = model2)
	pv = plot(symp_drive, lfhf, ylab = "LF/HF ratio", lab = "", title = "Sympathovagal Balance")
	pr = plot(symp_drive, rr, ylab = "RR [ms]", lab = "")
	p_BP = plot(symp_drive, sb, lab = "SB")
	plot!(symp_drive, db, lab = "DB", xlab = "sympathetic drive")
	#plot(pv, pr, p_BP, layout = (3,1))
	plot(pv, pr, layout = (2,1))
end

# ╔═╡ Cell order:
# ╠═f27d484a-0993-11eb-08c9-030d9f88b754
# ╟─60c79388-12b7-11eb-0493-51151fa697ab
# ╟─004f4ad1-4160-4e52-8813-70b9c73e3a29
# ╠═3d310136-0991-11eb-05ea-4351c72e2260
# ╟─7efa1b0e-0a5c-11eb-14c9-032c18080c9b
# ╟─25b183da-0998-11eb-0a1f-03534abcdc2e
# ╟─cdba0c98-120b-11eb-2aa3-d9e84e989773
# ╟─ac0c0478-0998-11eb-31bc-e3e4193eb47f
# ╟─f7a5c9ca-0998-11eb-1a7f-c57908a3cee5
# ╟─918ea734-120a-11eb-33e7-4340ad35c9f1
# ╟─ee91b962-12c3-11eb-24aa-ff59fbe0e53e
# ╟─ed8b4472-12c3-11eb-0d55-47162f88dc88
# ╟─eb927d16-120a-11eb-2a23-277da77b03f6
# ╠═78a242de-120b-11eb-1662-a948f6a40cf5
# ╟─a10eb0f2-50f1-11eb-0de4-f3e26e2ed897
# ╟─2bc64aaf-c07a-4e75-84e7-ce71f743f280
# ╟─23cc206d-78a3-4be0-b653-a8c5123d8dd9
# ╠═4935f00e-0991-11eb-2017-210b39b3aa56
# ╟─b01c33d6-0997-11eb-222c-d34da1f9aa4b
# ╟─e3c9adde-0993-11eb-191c-b501deb4a46f
# ╟─0f3b19aa-120c-11eb-3c2a-6f5e3e1e42fa
# ╠═4eb77de8-120d-11eb-3d3b-f1414ea32fef


text:
#=
pt=normpath(joinpath(pwd(),"testing","iterator.jl"))
@assert isfile(pt)
include(pt) =#

import Base: length,eltype,==,hash,iterate,append!,isless,resize!,convert

struct MyGrayCodeSubsetsHALF # <: DTSubsets
		size_of_set::Int
		stopat::Int # (1 << it.size_of_set)-2
	end

	eltype(it::MyGrayCodeSubsetsHALF) = Int # Array{eltype(1),1}
	length(it::MyGrayCodeSubsetsHALF) = it.stopat + 1 # 1 << it.size_of_set
	start(it::MyGrayCodeSubsetsHALF) = zero(1) # state = simply count through the subsets
	done(it::MyGrayCodeSubsetsHALF, state) = state > it.stopat  # if we have equality here, then the full set will be shifted / this is not needed for a binary split, hence we stop beforehand

	function next(it::MyGrayCodeSubsetsHALF, state)
	# returns the element which needs to be flipped
	  bit_to_flip = 1 + trailing_ones(state) # in the function print_graycode bit_to_flip is defined without the plus 1
	  state += one(1)
	  bit_to_flip, state
	end
@inline iterate(x::MyGrayCodeSubsetsHALF) = next(x, start(x))
@inline iterate(x::MyGrayCodeSubsetsHALF, i) = done(x, i) ? nothing : next(x, i)

	bitflip_graycode_subsetsHALF(xs) = MyGrayCodeSubsetsHALF(length(xs), (1 << (length(xs) - 1)) - 2)

subs = bitflip_graycode_subsetsHALF(["ab","cd","mixa"])

for i in subs
    @show i
end


text:
module DocumenterCitations

using Documenter
using Documenter.Anchors
using Documenter.Builder
using Documenter.Documents
using Documenter.Selectors
using Documenter.Utilities
using Documenter.Expanders

using Markdown
using Bibliography
using Bibliography: xnames, xyear, xlink, xtitle, xin
using Unicode

export CitationBibliography
struct CitationBibliography <: Documenter.Plugin
    bib::Dict
end
function CitationBibliography(filename::AbstractString="")
    filename == "" && return CitationBibliography(Dict())
    bf = import_bibtex(filename)
    return CitationBibliography(bf)
end

"""
    Example

Foo.

[Mayer2012](@cite)
"""
struct Example end

include("citations.jl")
include("bibliography.jl")

end


text:
const GENERATOR_MAPPING_FILE =
    joinpath(dirname(pathof(PowerSystems)), "parsers", "generator_mapping.yaml")

const STRING2FUEL = Dict((string(e) => e) for e in instances(ThermalFuels.ThermalFuel))
merge!(
    STRING2FUEL,
    Dict(
        "NG" => ThermalFuels.NATURAL_GAS,
        "NUC" => ThermalFuels.NUCLEAR,
        "GAS" => ThermalFuels.NATURAL_GAS,
        "OIL" => ThermalFuels.DISTILLATE_FUEL_OIL,
        "DFO" => ThermalFuels.DISTILLATE_FUEL_OIL,
        "SYNC_COND" => ThermalFuels.OTHER,
    ),
)

const STRING2PRIMEMOVER = Dict((string(e) => e) for e in instances(PrimeMovers.PrimeMover))
merge!(
    STRING2PRIMEMOVER,
    Dict(
        "W2" => PrimeMovers.WT,
        "WIND" => PrimeMovers.WT,
        "PV" => PrimeMovers.PVe,
        "RTPV" => PrimeMovers.PVe,
        "NB" => PrimeMovers.ST,
        "STEAM" => PrimeMovers.ST,
        "HYDRO" => PrimeMovers.HY,
        "ROR" => PrimeMovers.HY,
        "NUCLEAR" => PrimeMovers.ST,
        "SYNC_COND" => PrimeMovers.OT,
        "CSP" => PrimeMovers.CP,
        "UN" => PrimeMovers.OT,
        "STORAGE" => PrimeMovers.BA,
    ),
)

"""Return a dict where keys are a tuple of input parameters (fuel, unit_type) and values are
generator types."""
function get_generator_mapping(filename = nothing)
    if isnothing(filename)
        filename = GENERATOR_MAPPING_FILE
    end
    genmap = open(filename) do file
        YAML.load(file)
    end

    mappings = Dict{NamedTuple, DataType}()
    for (gen_type, vals) in genmap
        gen = getfield(PowerSystems, Symbol(gen_type))
        for val in vals
            key = (fuel = val["fuel"], unit_type = val["type"])
            if haskey(mappings, key)
                error("duplicate generator mappings: $gen $(key.fuel) $(key.unit_type)")
            end
            mappings[key] = gen
        end
    end

    return mappings
end

"""Return the PowerSystems generator type for this fuel and unit_type."""
function get_generator_type(fuel, unit_type, mappings::Dict{NamedTuple, DataType})
    fuel = uppercase(fuel)
    unit_type = uppercase(unit_type)
    generator = nothing

    # Try to match the unit_type if it's defined. If it's nothing then just match on fuel.
    for ut in (unit_type, nothing)
        key = (fuel = fuel, unit_type = ut)
        if haskey(mappings, key)
            generator = mappings[key]
            break
        end
    end

    if isnothing(generator)
        @error "No mapping for generator fuel=$fuel unit_type=$unit_type"
    end

    return generator
end

function get_branch_type(tap::Float64, alpha::Float64)
    if tap <= 0.0
        branch_type = Line
    elseif tap == 1.0
        branch_type = Transformer2W
    else
        if alpha == 0.0
            branch_type = TapTransformer
        else
            branch_type = PhaseShiftingTransformer
        end
    end

    return branch_type
end

function calculate_rating(active_power_max::Float64, reactive_power_max::Float64)
    return sqrt(active_power_max^2 + reactive_power_max^2)
end

function convert_units!(
    value::Float64,
    unit_conversion::NamedTuple{(:From, :To), Tuple{String, String}},
)

    if unit_conversion.From == "degree" && unit_conversion.To == "radian"
        value = deg2rad(value)
    elseif unit_conversion.From == "radian" && unit_conversion.To == "degree"
        value = rad2deg(value)
    else
        throw(DataFormatError("Unit conversion from $(unit_conversion.From) to $(unit_conversion.To) not supported"))
    end
    return value
end

function Base.convert(::Type{ThermalFuels.ThermalFuel}, fuel::AbstractString)
    return STRING2FUEL[uppercase(fuel)]
end

function Base.convert(::Type{ThermalFuels.ThermalFuel}, fuel::Symbol)
    return convert(ThermalFuels.ThermalFuel, string(fuel))
end

function Base.convert(::Type{PrimeMovers.PrimeMover}, primemover::AbstractString)
    return STRING2PRIMEMOVER[uppercase(primemover)]
end

function Base.convert(::Type{PrimeMovers.PrimeMover}, primemover::Symbol)
    return convert(PrimeMovers.PrimeMover, string(primemover))
end


text:
#Julia Gadfly Histogram

using Gadfly
plot(x=randn(113), Geom.histogram(bincount=10))

max_stars_repo_name: PacktPublishing/Learning-Jupyter

max_issues_repo_name: PacktPublishing/Learning-Jupyter

max_forks_repo_name: PacktPublishing/Learning-Jupyter