text:
using REPL, Base64

struct JunoDisplay <: AbstractDisplay end

function plotpane_io_ctx(io::IO)
  info = @rpc plotsize()
  if info isa Array
    size = info
    dpi = 1
  else
    size = info["size"]
    dpi = info["ratio"]
  end
  size .-= 1

  IOContext(io, :juno_plotsize => size, :juno_dpi_ratio => dpi, :juno_colors => syntaxcolors())
end

const plain_mimes = [
  "image/png",
  "image/jpeg",
  "image/tiff",
  "image/bmp",
  "image/gif"
]

const plotpane_mime = "application/prs.juno.plotpane+html"
const jlpane_mime = "application/prs.juno.jlpane"
using HTTP
function displayinplotpane(x)
  # allow openening a new pane in Atom
  if showable("application/prs.juno.jlpane", x)
    msg("jlpane", x.id, x.opts)
    return true
  end

  PlotPaneEnabled[] || return false

  legacy_plotpane = false
  # Juno-specific display always takes precedence
  if showable("application/juno+plotpane", x) && !showable(plotpane_mime, x)
    legacy_plotpane = true
    m = which(show, (IO, MIME"application/juno+plotpane", typeof(x)))
    @warn("""
      The \"application/juno+plotpane\" MIME type is deprecated. Please use \"$(plotpane_mime)\" instead.
    """, maxlog=1, _id=:juno_plotpane_legacy, _file=string(m.file), _line=m.line, _module=m.module)
  end
  if showable(plotpane_mime, x) || legacy_plotpane
    try
      io = IOBuffer()
      show(plotpane_io_ctx(io), legacy_plotpane ? "application/juno+plotpane" : plotpane_mime, x)
      str = base64encode(String(take!(io)))
      startswith(str, "data:") || (str = string("data:text/html;base64,", str))
      @msg ploturl(str)
      return true
    catch err
      err isa MethodError && err.f == Base.show || rethrow(err)
    end
  end
  # xml should be in a webview as well because of the possibility of embedded JS
  if showable("image/svg+xml", x)
    try
      io = IOBuffer()

      # render(PlotPane(), HTML(string("data:image/svg+xml,", stringmime("image/svg+xml", x, context=plotpane_io_ctx(io)))))
      io = IOBuffer()
      show(plotpane_io_ctx(io), "image/svg+xml", x)
      str = base64encode(String(take!(io)))
      startswith(str, "data:") || (str = string("data:image/svg+xml;base64,", str))
      @msg ploturl(str)
      # render(PlotPane(), HTML(""))
      return true
    catch err
      err isa MethodError && err.f == Base.show || rethrow(err)
    end
  end

  # images
  for mime in plain_mimes
    if showable(mime, x)
      try
        io = IOBuffer()
        suffix = istextmime(mime) ? "," : ";base64,"
        str = stringmime(mime, x, context=plotpane_io_ctx(io))
        str = string("")
        render(PlotPane(), HTML(str))
        return true
      catch err
        err isa MethodError && err.f == Base.show || rethrow(err)
      end
    end
  end
  return false
end

using TreeViews: hastreeview, numberofnodes, treelabel, treenode, nodelabel

# called in user code and in REPL
function Base.display(d::JunoDisplay, m::MIME{T}, x) where T
  displayable(d, m) || throw(MethodError(display, (d, m, x)))
  Base.display(d::JunoDisplay, x)
end

Base.displayable(d::JunoDisplay, m::MIME"application/juno+inline") = true
Base.displayable(d::JunoDisplay, m::MIME"application/prs.juno.inline") = true

Base.displayable(d::JunoDisplay, m::MIME"application/juno+plotpane") = PlotPaneEnabled[]
Base.displayable(d::JunoDisplay, m::MIME"application/prs.juno.plotpane+html") = PlotPaneEnabled[]
Base.displayable(d::JunoDisplay, m::MIME"image/svg+xml") = PlotPaneEnabled[]
Base.displayable(d::JunoDisplay, m::MIME"application/prs.juno.jlpane") = true

function Base.displayable(d::JunoDisplay, m::MIME{T}) where T
  string(T) in plain_mimes ? PlotPaneEnabled[] : false
end

function Base.display(d::JunoDisplay, x)
  d = last(filter(x -> (x isa REPL.REPLDisplay), Base.Multimedia.displays))
  if displayinplotpane(x)
    # we shouldn't need to do this, but Plots.jl has a ugly overload on display(::REPLDisplay, ::MIME"text/plain" ::Plot)
    # which would otherwise be used
    inREPL[] && invoke(display, Tuple{typeof(d), typeof(MIME"text/plain"()), Any}, d, MIME"text/plain"(), x)
    # throw(MethodError(display, "nope"))
  else
    throw(MethodError(display, (x,)))
  end
  return nothing
end

# called when displaying results in the editor
function displayandrender(res)
  if !displayinplotpane(res)
    old = customdisplaystack()
    try
      display(res)
    catch e
    end
    restoredisplaystack(old)
  end

  if hastreeview(res)
    res = generateTreeView(res)
  end
  render′(Editor(), res)
end

function customdisplaystack()
  old = copy(Base.Multimedia.displays)
  filter!(Base.Multimedia.displays) do d
    !(d isa REPL.REPLDisplay || d isa TextDisplay ||
      d isa JunoDisplay || d isa Media.DisplayHook)
  end
  old
end

function restoredisplaystack(old)
  empty!(Base.Multimedia.displays)
  foreach(pushdisplay, old)
end

function bestlabel(f, args...)
  if applicable(f, IOBuffer(), args..., MIME"application/prs.juno.inline"())
    return HTML(sprint(f, args..., MIME"application/prs.juno.inline"()))
  elseif applicable(f, IOBuffer(), args..., MIME"application/juno+inline"()) # deprecated
    m = which(f, (IO, typeof.(args)..., MIME"application/juno+inline"));
    @warn("""
      The \"application/juno+inline\" MIME type is deprecated. Please use \"application/prs.juno.inline\" instead.
    """, maxlog=1, _id=:juno_tree_legacy, _file=string(m.file), _line=m.line, _module=m.module);
    return HTML(sprint(f, args..., MIME"application/juno+inline"()))
  elseif applicable(f, IOBuffer(), args..., MIME"text/html"())
    return HTML(sprint(f, args..., MIME"text/html"()))
  else
    return Text(sprint(f, args..., MIME"text/plain"()))
  end
end

# TreeViews.jl -> Juno.LazyTree
function generateTreeView(x)
  header = bestlabel(treelabel, x)

  numberofnodes(x) == 0 &&  return Tree(header, [])

  genchildren = function ()
    children = Any[]
    for i in 1:numberofnodes(x)
      node = treenode(x, i)
      cheader = bestlabel(nodelabel, x, i)

      if isempty(cheader.content)
        node === missing && continue
        push!(children, hastreeview(node) ? generateTreeView(node) : node)
      elseif node === missing
        push!(children, cheader)
      else
        push!(children, SubTree(Row(cheader, text" → "), hastreeview(node) ? generateTreeView(node) : node))
      end
    end
    children
  end

  return LazyTree(header, genchildren)
end


text: run(`make all`)

ext: jl


text:
# implicit scheme for advection-diffusion

using Revise
using AdFem
using ADCME
using PyPlot
using Statistics
using LinearAlgebra
using ADCMEKit

m = 40
n = 20
h = 1/n
NT = 100
Δt = 1/NT

u = 0.5*[ones(m*n);zeros(m*n)]

up_and_down = bcedge("upper|lower", m, n, h)

M = compute_fvm_mass_matrix(m, n, h)
K, rhs1 = compute_fvm_advection_matrix(u, up_and_down, zeros(size(up_and_down, 1)), m, n, h)
S, rhs2 = compute_fvm_tpfa_matrix(missing, up_and_down, zeros(size(up_and_down, 1)), m, n, h)

A = M/Δt + K - 0.01*S
A = factorize(A)

U = zeros(m*n, NT+1)
xy = fvm_nodes(m, n, h)
u0 = @. exp( - 10 * ((xy[:,1]-1.0)^2 + (xy[:,2]-0.5)^2))

function condition(i, args...)
    i <= NT
end

function body(i, u_arr)
    u = read(u_arr, i)
    u_arr = write(u_arr, i+1, A\(M*u/Δt - rhs1 + rhs2))
    return i+1, u_arr
end

i = constant(1, dtype = Int32)
u_arr = TensorArray(NT+1)
u_arr = write(u_arr, 1, u0)
_, u = while_loop(condition, body, [i, u_arr])
u = set_shape(stack(u), (NT+1, m*n))

sess = Session(); init(sess)
U = run(sess, u)
Z = zeros(NT+1, n, m)
for i = 1:NT+1
    Z[i,:,:] = reshape(U[i,:], m, n)'
end
p = visualize_scalar_on_fvm_points(Z, m, n, h)
saveanim(p, "advec.gif")

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages


text:
# Functions related to the counting bound.
# This version uses log-transformed numbers (but doesn't
# use analytic approximations to binomial coefficients).
# All logs are natural, unless otherwise noted.

module LogCountingBound

export logNumMaximalCliques

"""
  Log of the approximate number of maximal hypercliques of some size.
  Note that k < r < n .
  Also, the precision of what's returned can be set by setprecision().
  k: number of vertices per hyperedge
  r: number of vertices in the clique
  n: vertices in the larger graph
  Returns: (natural log of) the number of maximal hypercliques.
		This is approximate, as it's the expected number, and multiplies
		floating point by adding logarithms.
"""
function logNumMaximalCliques(k, r, n)
  k = BigInt(k)
  r = BigInt(r)
  n = BigInt(n)

  # log of probability that one of those is
	# not covered by a larger clique
  a = BigInt(2) ^ binomial(r, k-1)
	# ??? will this have enough significant figures?
	logP = (log(a-1) - log(a)) * (n - r)

    # number of possible cliques of this size
		(log(binomial(n, r))
    # P(clique is present ...)
		+ (- log(BigInt(2)) * binomial(r, k))
    # P( ... and not covered by a larger clique)
    + logP)
end

end

max_stars_repo_name: joshtburdick/misc

max_issues_repo_name: joshtburdick/misc

max_forks_repo_name: joshtburdick/misc


text:
import DataFrames
import Random
import StatsBase

"""
"""
function shuffle_rows!(
        dataframe::DataFrames.AbstractDataFrame,
        )
    result = shuffle_rows!(Random.GLOBAL_RNG, dataframe)
    return result
end

"""
"""
function shuffle_rows!(
        rng::Random.AbstractRNG,
        dataframe::DataFrames.AbstractDataFrame,
        )
    numrows = size(dataframe,1)
    allrows = convert(Array,1:numrows)
    rowpermutation = Random.shuffle!(rng, allrows)
    numcolumns = size(dataframe,2)
    for j = 1:numcolumns
        dataframe[:, j] = dataframe[rowpermutation, j]
    end
    return dataframe
end


text:
module Graphs

using HyperGraphs
import Graphs as _Graphs

const AbstractSimpleGraph = _Graphs.SimpleGraphs.AbstractSimpleGraph
const AbstractSimpleEdge = _Graphs.SimpleGraphs.AbstractSimpleEdge

## Graphs to HyperGraphs ##

function Base.convert(::Type{T}, g::U) where {T<:AbstractHyperGraph, U<:AbstractSimpleGraph}
    ex = collect(_Graphs.edges(g))
    etype = HyperGraphs.get_hyperedge_type(T)
    if _Graphs.is_directed(g) && isoriented(T)
        T([etype([e.src], [e.dst]) for e in ex])
    elseif !_Graphs.is_directed(g) && !isoriented(T)
        T([etype([e.src, e.dst]) for e in ex])
    else
        error("orientation of types does not match")
    end
end

## HyperGraphs to Graphs ##

# convert a HyperGraph to a Graphs.SimpleGraph
function Base.convert(::Type{_Graphs.SimpleGraph}, x::HyperGraph{T}) where {T<:Int}
    !issimple(x) && error("hypergraph is not simple")
    !isgraph(x) && error("some cardinalities != 2")
    !isequal(collect(1:maximum(vertices(x))), vertices(x)) && error("non-consecutive vertices")
    _Graphs.SimpleGraph([map((v1, v2) -> _Graphs.Edge(v1 => v2), vertices(e)...) for e in hyperedges(x)])
end

end # module


text:
mutable struct JGraphData
    graph::AbstractGraph
    layout
    node_templates
    edge_color
    edge_size
    width
    height
    frames
    numbered
end

"""
    JGraphData(graph::AbstractGraph, layout::NetworkLayout.AbstractLayout; kwargs...)

The `JGraphData` is the structure used to set all the properties of the graph animation,
the graph itself, its layout and all the options useful to animate it in the appropriate way.
It should be passed to `JGraph` to obtain an animation with the specified properties.
The keyword arguments will increase to increase in the future.

# Arguments
- `graph::AbstractGraph`: A graph as created with `LightGraphs` and affiliate packages
- `layout::NetworkLayout.AbstractLayout`: A layout as provided by `NetworkLayout` to place the graph nodes in space.

# Keywords
- `node_template` provide function that could be used as the function argument of an `Object`
- `edges_color` sets the colors of all edges
- `edges_size` sets the width of all edges
- `width` sets the width to which the graph will be adjusted
- `height` sets  the height to which the graph will be adjusted
- `frames` the frames during which the animations on the graph will last if not specified elsewhere
- `numbered` Defaults to `false` if set to `true` labels the nodes with their number as vertices in the graph.
"""
function JGraphData(
    graph::AbstractGraph,
    layout = NetworkLayout.Shell();
    node_templates = JCircle(O, 3, color=colorant"black", action=:fill),
    edge_color = colorant"black",
    edge_size = 1,
    width = 300,
    height = 300,
    frames = Javis.get_frames(Javis.CURRENT_VIDEO[1].objects[1]),
    numbered = false,
)

    return JGraphData(
        graph,
        layout,
        node_templates,
        edge_color,
        edge_size,
        width,
        height,
        frames,
        numbered,
    )
end

mutable struct JGraph
    data::JGraphData
    jnodes::Vector{Object}
    jedges::Dict{Pair{Int64,Int64},Object}
end

function _JGraph(g::JGraphData)

    # points = [g.scaling * GB2Luxor(point) for point in g.layout(g.graph)]
    points = _adapt_to_lims(g.layout(g.graph), g.width, g.height)

    Jnodes = if isa(g.node_templates, Function)
        [
            Object(:same, g.node_templates, point) for
            point in points
        ]

    elseif length(g.node_templates) == 1
        [
            Object(:same, g.node_templates[1], point) for
            point in points
        ]

    else
        length(g.node_templates) != length(points) &&
            error("Either a single function or one for each node should be provided as node_template")

        [
            Object(:same, node_template, point) for
            (node_template, point) in zip(g.node_templates, points)
        ]
    end

    Jedges = Dict([
        (e.src => e.dst) => Object(
            (args...; color=colorant"black") -> begin
                setline(g.edge_size)
                sethue(color)
                line(pos(Jnodes[e.src]), pos(Jnodes[e.dst]), :stroke)
            end,
        ) for e in edges(g.graph)
    ])

    if g.numbered
        for (idx, jnode) in enumerate(Jnodes)
            Object((args...) -> begin
                label(string(idx), :N, pos(jnode))
            end)
        end
    end

    return JGraph(g, Jnodes, Jedges)
end

"""
    JGraph(g::JGraphData)

The `JGraph` structure is the core structure of the package when called
with a `JGraphData` argument creates and `Object` for each node in the graph
and stores them in `jnodes` and an `Object` for each edge and stores them in `jedges`.
Note that `jnodes` is a `Vector`, whereas `jedges` is a `Dict` with pairs of nodes (integers)
as keys.

!!! note
    Since `JGraph` stores `Object`s it can only be called after a `Video`
    has been initialized look at `Javis` documentation to know all about it.
"""
function JGraph(g)
    _JGraph(g)
end

"""
    jnodes(jg::JGraph)

Accessor function to the `jnodes` stored in a `JGraph`.
"""
function jnodes(jg::JGraph)
    jg.jnodes
end

"""
    jedges(jg::JGraph)

Accessor function to the `jedges` stored in a `JGraph`.
"""
function jedges(jg::JGraph)
    jg.jedges
end

"""
    get_starting_positions(g::JGraph)

Returns the positions as obtained using the layout and the graph stored
in `g.data`.
"""
function get_starting_positions(g::JGraph)
    [node.start_pos for node in g.jnodes]
end


text:
include(joinpath(@__DIR__, "pedagogy_graph.jl"))
using Random

dims = [5, 5]
g = LG.grid(dims)
β = 1.0

hcg = PRS.HardCoreGraph(g, β)

rng = Random.MersenneTwister(123)
@time sample = PRS.generate_sample_prs(rng, hcg)

path = joinpath(@__DIR__, "output", "hard_core_graph", "hard_core_graph.pdf")
p = plot(hcg, sample, dims, path)

# @time sample = PRS.generate_sample_prs(hcg)


text:
using JuMP

type Recipe
	name::AbstractString
	time::Float64
	outs::Dict{AbstractString, Float64}
	ins::Dict{AbstractString, Float64}
	function Recipe(name, time, speed)
		recips[name] = new(name, time/speed, Dict{AbstractString, Float64}(), Dict{AbstractString, Float64}())
		recips[name]
	end
end

recips = Dict{AbstractString, Recipe}()

function production(r::Recipe, nfacs)
	@printf "%s: %d\n" r.name nfacs
	@printf "\tProduces\n"
	for (k,v) in r.outs
		@printf "\t\t%0.2f \"%s\" per min - %f per second\n" 60v * nfacs k v * nfacs
	end
	@printf "\tConsumes\n"
	for (k,v) in r.ins
		@printf "\t\t%0.2f \"%s\" per min - %f per second\n" 60v * nfacs k v * nfacs
	end
end

macro IN(iname, iqty)
	return quote
		recipe.ins[$iname] = $iqty / recipe.time
	end
end

macro OUT(iname, iqty)
	return quote
		recipe.outs[$iname] = $iqty / recipe.time
	end
end

@enum Rnum AOP HOC LOC SOL ROK

recips[AOP] = Recipe("Advanced Oil Processing", 5.0, 1.)
@IN AOP "Crude Oil" 200.
@IN AOP "Water" 100.
@OUT AOP "Heavy Oil" 10.
@OUT AOP "Light Oil" 45.
@OUT AOP "Petroleum Gas" 55.

recips[HOC] = Recipe("Heavy Oil Cracking", 3.0 , 1.25)
@IN HOC "Heavy Oil" 40.
@IN HOC "Water" 30.
@OUT HOC "Light Oil" 30.

recips[LOC] = Recipe("Light Oil Cracking", 5.0 , 1.25)
@IN LOC "Light Oil" 30.
@IN LOC "Water" 30.
@OUT LOC "Petroleum Gas" 20.

recips[SOL] = Recipe("Solid Fuel", 3., 1.25)
@IN SOL "Petroleum Gas" 20.
@OUT SOL "Solid Fuel" 1.

recips[ROK] = Recipe("Rocket Fuel", 30., 1.25)
@IN ROK "Solid Fuel" 10.
@OUT ROK "Rocket Fuel" 1.

m = Model()

@variable(m, AOPs >= 0, Int)
@variable(m, HOCs >= 0, Int)
@variable(m, LOCs >= 0, Int)
@variable(m, SOLs >= 0, Int)
@variable(m, ROKs >= 0, Int)

@constraint(m, AOPs * recips[AOP].ins["Crude Oil"] <= 2400)
@constraint(m, HOCs * recips[HOC].ins["Heavy Oil"] <= AOPs * recips[AOP].outs["Heavy Oil"])
@constraint(m, LOCs * recips[LOC].ins["Light Oil"] <= AOPs * recips[AOP].outs["Light Oil"] + HOCs * recips[HOC].outs["Light Oil"])
@constraint(m, SOLs * recips[SOL].ins["Petroleum Gas"] <= AOPs * recips[AOP].outs["Petroleum Gas"] + LOCs * recips[LOC].outs["Petroleum Gas"])
@constraint(m, ROKs * recips[ROK].ins["Solid Fuel"] <= SOLs * recips[SOL].outs["Solid Fuel"])

@objective(m, Max, ROKs)

solve(m)

production(recips[AOP], getvalue(AOPs))
production(recips[HOC], getvalue(HOCs))
production(recips[LOC], getvalue(LOCs))
production(recips[SOL], getvalue(SOLs))
production(recips[ROK], getvalue(ROKs))

max_stars_repo_name: lawless-m/factorio-utils

max_issues_repo_name: lawless-m/factorio-utils

max_forks_repo_name: lawless-m/factorio-utils


text:
using UnicodePlots
using SpikingNeuralNetworks
SNN = SpikingNeuralNetworks
using SpikeSynchrony
using Statistics
using JLD
using Plots
using UnicodePlots
using Evolutionary
using Random
##

# Override to function to include a state.
##

const Ne = 200;
const Ni = 50
const seed = 10
const k = 0.5
const this_size = 50

const Ne = 200;
const Ni = 50
const σee = 1.0
const pee = 0.5
const σei = 1.0
const pei = 0.5

global E
global spkd_ground

function make_net_from_graph_structure(arbitrary_int::Int)
    """
    # Circular ladder simulation
    Input integers affect the number of network nodes (neurons) in the circular ladder network.
    """

    arbitrary_int = Int(round(arbitrary_int))
    h = circular_ladder_graph(arbitrary_int)
    hi = circular_ladder_graph(arbitrary_int)
    E = SNN.IZ(; N = Ne, param = SNN.IZParameter(; a = 0.02, b = 0.2, c = -65, d = 8))
    I = SNN.IZ(; N = Ni, param = SNN.IZParameter(; a = 0.1, b = 0.2, c = -65, d = 2))
    EI = SNN.SpikingSynapse(E, I, :v; σ = σei, p = 1.0)
    IE = SNN.SpikingSynapse(I, E, :v; σ = -1.0, p = 1.0)
    II = SNN.SpikingSynapse(I, I, :v; σ = -1.0, p = 1.0)
    P = [E, I]
    C = [EI, IE, II]
    EE = SNN.SpikingSynapse(E, E, :v; σ = 0.5, p = 0.8)
    @inbounds for (i, j) in enumerate(h.fadjlist)
        @inbounds for k in j
            SNN.connect!(EE, i, k, 10)
        end
    end
    @inbounds for (i, j) in enumerate(hi.fadjlist)
        @inbounds for k in j
            if i < Ni && k < Ni
                SNN.connect!(EI, i, k, 10)
                SNN.connect!(IE, i, k, 10)
                SNN.connect!(II, i, k, 10)
            end
        end
    end

    P = [E, I]
    C = [EE, EI, IE, II]
    return P, C

end

function make_net_SNN(Ne, Ni; σee = 1.0, pee = 0.5, σei = 1.0, pei = 0.5)
    Ne = 200
    Ni = 50

    E = SNN.IZ(; N = Ne, param = SNN.IZParameter(; a = 0.02, b = 0.2, c = -65, d = 8))
    I = SNN.IZ(; N = Ni, param = SNN.IZParameter(; a = 0.1, b = 0.2, c = -65, d = 2))
    EE = SNN.SpikingSynapse(E, E, :v; σ = σee, p = pee)
    EI = SNN.SpikingSynapse(E, I, :v; σ = σei, p = pei)
    IE = SNN.SpikingSynapse(I, E, :v; σ = -1.0, p = 0.5)
    II = SNN.SpikingSynapse(I, I, :v; σ = -1.0, p = 0.5)
    P = [E, I]
    C = [EE, EI, IE, II]
    @show(C)

    return P, C
end
function get_trains(p)
    fire = p.records[:fire]
    x, y = Float32[], Float32[]
    for time in eachindex(fire)
        for neuron_id in findall(fire[time])
            push!(x, time)
            push!(y, neuron_id)
        end
    end
    cellsa = Array{Union{Missing,Any}}(undef, 1, Int(findmax(y)[1]))
    nac = Int(findmax(y)[1])
    for (inx, cell_id) in enumerate(1:nac)
        cellsa[inx] = []
    end
    @inbounds for cell_id in unique(y)
        @inbounds for (time, cell) in collect(zip(x, y))
            if Int(cell_id) == cell
                append!(cellsa[Int(cell_id)], time)

            end

        end
    end

    cellsa

end

function rmse(spkd)
    error = Losses(mean(spkd), spkd; agg = mean)
end

function rmse_depr(spkd)
    total = 0.0
    @inbounds for i = 1:size(spkd, 1)
        total += (spkd[i] - mean(spkd[i]))^2.0
    end
    return sqrt(total / size(spkd, 1))
end

function raster_difference(spkd0, spkd_found)
    maxi0 = size(spkd0)[2]
    maxi1 = size(spkd_found)[2]
    mini = findmin([maxi0, maxi1])[1]
    spkd = ones(mini)
    maxi = findmax([maxi0, maxi1])[1]

    if maxi > 0
        if maxi0 != maxi1
            return sum(ones(maxi))

        end
        if isempty(spkd_found[1, :])
            return sum(ones(maxi))
        end
    end
    spkd = ones(mini)
    @inbounds for i in eachindex(spkd)
        if !isempty(spkd0[i]) && !isempty(spkd_found[i])
            maxt1 = findmax(spkd0[i])[1]
            maxt2 = findmax(spkd_found[i])[1]
            maxt = findmax([maxt1, maxt2])[1]
            if maxt1 > 0.0 && maxt2 > 0.0
                t, S = SpikeSynchrony.SPIKE_distance_profile(
                    unique(sort(spkd0[i])),
                    unique(sort(spkd_found[i]));
                    t0 = 0.0,
                    tf = maxt,
                )
                spkd[i] = SpikeSynchrony.trapezoid_integral(t, S) / (t[end] - t[1])
            end
        end
    end
    spkd
end

function loss(model)
    @show(Ne, Ni)
    @show(model)

    σee = model[1]
    pee = model[2]
    σei = model[3]
    pei = model[4]
    P1, C1 = make_net_SNN(Ne, Ni, σee = σee, pee = pee, σei = σei, pei = pei)#,a=a)
    @show(C1)
    E1, I1 = P1
    SNN.monitor([E1, I1], [:fire])
    sim_length = 500
    @inbounds for t = 1:sim_length*ms
        E1.I = vec([11.5 for i = 1:sim_length])
        SNN.sim!(P1, C1, 1ms)
    end

    spkd_found = get_trains(P1[1])
    println("Ground Truth \n")
    SNN.raster([E]) |> display
    println("Best Candidate \n")
    SNN.raster([E1]) |> display

    error = raster_difference(spkd_ground, spkd_found)
    error = sum(error)
    @show(error)

    error

end

function eval_best(params)
    """
    Evaluate the best candidate network model chosen by a genetic algorithm.
    """

    σee = params[1]
    pee = params[2]
    σei = params[3]
    pei = params[4]
    P1, C1 = make_net_SNN(Ne, Ni, σee = σee, pee = pee, σei = σei, pei = pei)
    E1, I1 = P1
    SNN.monitor([E1, I1], [:fire])
    sim_length = 500
    @inbounds for t = 1:sim_length*ms
        E1.I = vec([11.5 for i = 1:sim_length])
        SNN.sim!(P1, C1, 1ms)
    end

    spkd_found = get_trains(P1[1])
    println("Ground Truth: \n")
    SNN.raster([E]) |> display
    println("candidate: \n")
    SNN.raster([E1]) |> display
    E1, spkd_found

end

function Evolutionary.trace!(
    record::Dict{String,Any},
    objfun,
    state,
    population,
    method::GA,
    options,
)
    idx = sortperm(state.fitpop)
    record["fitpop"] = state.fitpop[:]
    record["pop"] = population[:]
end


text:
function movements_to_side!(r, side :: HorizonSide)
    count = 0
    while isborder(r,side) == false
        move!(r,side)
        count = count + 1
    end
    return count
end

function ugol!(r)

    home = Int[]

    while (isborder(r, Sud) == false) || (isborder(r, West) == false)
        s = movements_to_side!(r, Sud)
        w = movements_to_side!(r, West)
        push!(home, s)
        push!(home, w)
    end
    return home
end

function count_move!(r, side :: HorizonSide, count)

    for i = 1:count
        move!(r, side)
    end
end

function go_home!(r, home)
    reverse!(home)

    for i in 1: length(home)
        if i % 2 == 0
            count_move!(r, Nord, home[i])
        else
            count_move!(r, Ost, home[i])
        end
    end
end

function corners(r)

    home = Int[]
    home = ugol!(r)

    for side in (Nord, Ost, Sud, West)
        movements_to_side!(r, side)
        putmarker!(r)
    end

    go_home!(r,home)

end

corners(r)
show(r)

max_stars_repo_name: AntonDudDm/mirea-progs

max_issues_repo_name: AntonDudDm/mirea-progs

max_forks_repo_name: AntonDudDm/mirea-progs


text:
using LVServer
using Test

@testset "LVServer.jl" begin
    # Write your tests here.

    # Here are the tests:
    @test LVServer.tmp_test(0) == 0
    @test LVServer.tmp_test(3) == 3
    @test LVServer.tmp_test(5) == 5
end


text:
function discretizecns!(system::CVSystem)
    append!(system.cns.H,zeros(system.solverparams.numsteps+1))
    append!(system.cns.Emaxlv,zeros(system.solverparams.numsteps+1))
    append!(system.cns.Emaxrv,zeros(system.solverparams.numsteps+1))
    append!(system.cns.R2L,zeros(system.solverparams.numsteps+1))
    append!(system.cns.R3L,zeros(system.solverparams.numsteps+1))
    append!(system.cns.R2U,zeros(system.solverparams.numsteps+1))
    append!(system.cns.R3U,zeros(system.solverparams.numsteps+1))
    append!(system.cns.C4L,zeros(system.solverparams.numsteps+1))
    append!(system.cns.C5L,zeros(system.solverparams.numsteps+1))
    append!(system.cns.C4U,zeros(system.solverparams.numsteps+1))
    append!(system.cns.C5U,zeros(system.solverparams.numsteps+1))
    append!(system.cns.V4L,zeros(system.solverparams.numsteps+1))
    append!(system.cns.V5L,zeros(system.solverparams.numsteps+1))
    append!(system.cns.V4U,zeros(system.solverparams.numsteps+1))
    append!(system.cns.V5U,zeros(system.solverparams.numsteps+1))
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
using Documenter
using Practice_Julia_Repositoria

makedocs(
    sitename = "Practice_Julia_Repositoria",
    format = Documenter.HTML(),
    modules = [Practice_Julia_Repositoria]
)

deploydocs(
    repo = "github.com/sandreza/Practice_Julia_Repositoria"
)


text:
# Markov Chain Monte Carlo (MCMC)
# Parallel Processing

using PyPlot
using Distributed
using Dates
using LinearAlgebra
using DelimitedFiles
using Random

struct Metropolis_Spec
    rl_impl::Function
    model_params::Array{MetropolisParam,1}
    nr_burnin::Int
    nr_samples::Int
    inter_collect_interval::Int
end

function getProposalDistribution(params)
    allVariances = [p.rwalk_step_size_std^2 for p in params]
    return Distributions.MvNormal(zeros(length(params)), diagm(allVariances))
end

function run_Metropolis(data::Array{Int,2}, mSpec, opts...)
    # Perform MCMC
    proposal_dist = getProposalDistribution(mSpec.model_params)
    states = zeros(mSpec.nr_samples, length(mSpec.model_params)+1)

    currentValues = [rand(p.prior_Distribution) for p in mSpec.model_params]

    index_alpha_0 = [i for i in 1:length(mSpec.model_params) if mSpec.model_params[i].name == "alpha_0"]
    if !isempty(index_alpha_0); index_alpha_0 = index_alpha_0[1]; end

    # ----- Make sure alpha_0 + alpha_surprise < 1
    # NOT USED eventually
    for i in 1:length(mSpec.model_params)
        if ((mSpec.model_params[i].name == "alpha_surprise")
            && !isempty(index_alpha_0) )
            while (currentValues[i] + currentValues[index_alpha_0]) > 1.
                currentValues[i] = rand(mSpec.model_params[i].prior_Distribution)
                currentValues[index_alpha_0] = rand(mSpec.model_params[index_alpha_0].prior_Distribution)
            end
        end
    end

    llPrior = sum([log(pdf(mSpec.model_params[i].prior_Distribution, currentValues[i])) for i in length(mSpec.model_params)])
    currentLL = llPrior + mSpec.rl_impl(data, currentValues, opts...)[1]

    nr_collected_samples = 0
    step_counter = 0
    while nr_collected_samples<mSpec.nr_samples
        if (step_counter > mSpec.nr_burnin && (step_counter%mSpec.inter_collect_interval ==0))
          nr_collected_samples += 1
          states[nr_collected_samples,:] = vcat(currentValues, currentLL)
        end
        proposalValues = currentValues + rand(proposal_dist)

        llPrior = 0.
        for i in 1:length(mSpec.model_params)
            if typeof(mSpec.model_params[i].prior_Distribution) == DiscreteUniform
                # Handle discrete variables as well)
                proposalValues[i] = Int(round(proposalValues[i]))
            end
            # ----- Make sure alpha_0 + alpha_surprise < 1
            # NOT USED eventually
            if ((mSpec.model_params[i].name == "alpha_surprise")
                && !isempty(index_alpha_0)
                && (proposalValues[i] > (1. - proposalValues[index_alpha_0])))
                llPrior += -Inf
            else
                llPrior += log(pdf(mSpec.model_params[i].prior_Distribution, proposalValues[i]))
            end
        end

        if llPrior > currentLL #prior can be -Inf for out of bound steps
            proposalLL = llPrior + mSpec.rl_impl(data, proposalValues, opts...)[1]
          if exp(proposalLL - currentLL) > rand()
              currentValues = proposalValues
              currentLL = proposalLL
          end
        else

        end

      step_counter+=1

    end
    println("llPrior=$llPrior")
    return states
end

function test_work_remotely(msg::String)
    println("worker proc id = $(myid()), called with echo message: $msg")
    RL_Fit.sayHi()
    a = rand(2, 3) .+ myid()
    return a
end

function run_parallel_Metropolis(data::Array{Int,2}, mSpec,rl_module_path::String,
                                opts...; nr_of_repetitions=2)
    # Run MCMC on many threads on server in parallel
    println("run_parallel_Metropolis")

    res_cat = @distributed vcat for i=1:nr_of_repetitions
        run_Metropolis(data, mSpec, opts...)
        end

    return res_cat
end

# specify CV for multiple models
struct Model_spec
    model_name::String
    learner_function_V::Function
    model_params::Array{MetropolisParam,1}
    result_folder::String
end

struct Model_spec_per_subj
    model_name::String
    learner_function_V::Function
    model_params::Array{MetropolisParam,1}
    per_param_posterior_std::Array{Float64,1}
    result_folder::String
end

# used to do a "batch" fit of multiple models
struct FitSpec_multiple_models
    all_models::Array{Any,1}  # allow  Metropolis_Spec or Metropolis_Spec_per_subj, TODO: learn inheritance in julia
    fit_max_episodes::Int
    nr_burnin::Int
    nr_samples::Int
    inter_collect_interval::Int
    top_n_fraction::Float64
    nr_repetitions::Int
    do_per_subject_fit::Bool
    log_path::String
end

function run_fit_multiple_models(description::String, data_path::String,
                                fit_spec::FitSpec_multiple_models, rl_module_path)
    fit_start_time = now()
    log_buffer = IOBuffer()
    log_print_ln!(log_buffer, "\nParameter fitting: start: $fit_start_time")
    log_print_ln!(log_buffer, "description: $description")
    plt.ioff() # don't plot, save figure

    (all_data, all_data_stats) = RL_Fit.load_SARSPEI(data_path; max_episodes_per_subj = fit_spec.fit_max_episodes)
    data_string = RL_Fit.to_short_string(all_data_stats)
    log_print_ln!(log_buffer, "Data:\n$data_string")
    nr_models = length(fit_spec.all_models)
    log_print_ln!(log_buffer, "nr of models: $nr_models")
    @show size(all_data)
    for i_model = 1:nr_models
        current_model = fit_spec.all_models[i_model]
        log_print_ln!(log_buffer, "********************************************************************************")
        log_print_ln!(log_buffer, "*********       Start parameter fitting of model $(current_model.model_name)       ************")
        log_print_ln!(log_buffer, "********************************************************************************")
        log_print_ln!(log_buffer, "Prior:")
        log_print_ln!(log_buffer, "$(current_model.model_params)")
        mSpec = RL_Fit.Metropolis_Spec(current_model.learner_function_V, current_model.model_params, fit_spec.nr_burnin, fit_spec.nr_samples, fit_spec.inter_collect_interval)
        # fit on all data (typically all subjects, max 40 episodes per subject)
        top_MAP, top_n_avgMAP, top_n_min, top_n_max, top_n_stdMAP, top_n_samples,
            mcmc_samples, LL_noPrior_MAP, LL_noPrior_avgMAP, all_signals = RL_Fit.run_fit_procedure(description, all_data, mSpec,
                                                                    fit_spec.nr_repetitions,  "TODO: mcmc_data_file",
                                                                    fit_spec.top_n_fraction, rl_module_path)
        log_print_ln!(log_buffer, "top_MAP: $top_MAP")
        log_print_ln!(log_buffer, "top_n_avgMAP: $top_n_avgMAP")
        log_print_ln!(log_buffer, "size(mcmc_samples): $(size(mcmc_samples))")
        log_print_ln!(log_buffer, "LL(top_MAP): $LL_noPrior_MAP")
        log_print_ln!(log_buffer, "LL(top_n_avgMAP): $LL_noPrior_avgMAP")

        top_n_diff = (top_n_max - top_n_min)[1:end-1]
        if maximum(top_n_diff) > 0.15
            log_print_ln!(log_buffer, "   >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>")
            log_print_ln!(log_buffer, "       !! CHECK LARGE DIFFERENCE IN TOP-N PARAM VALUES !!")
            log_print_ln!(log_buffer, "       top_n_diff=$top_n_diff")
            log_print_ln!(log_buffer, "       top_n_min: $top_n_min")
            log_print_ln!(log_buffer, "       top_n_max: $top_n_max")
            log_print_ln!(log_buffer, "   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<")
        end

        nr_params = length(current_model.model_params)
        nr_subjects = length(all_data_stats.all_subject_ids)
        AIC = 2.0 * nr_params - 2.0 * LL_noPrior_MAP
        AICc = AIC + (2*nr_params)*(nr_params+1)/(nr_subjects- nr_params - 1)

        # Store the map parameters and the log likelihood
        LL_path = "$(current_model.result_folder)maxLL_fit/"
        if !ispath(LL_path)
            println("mkpath($LL_path)")
            mkpath("$LL_path")
        end
        writedlm("$(LL_path)LogLikelihood_MAP.csv", LL_noPrior_MAP)
        writedlm("$(LL_path)LogLikelihood_avgMAP.csv", LL_noPrior_avgMAP)
        writedlm("$(LL_path)MAP_params.csv", top_MAP[1:end-1])
        writedlm("$(LL_path)MAP_n_avg_params.csv", top_n_avgMAP[1:end-1])
        writedlm("$(LL_path)AIC.csv", AIC)
        writedlm("$(LL_path)AICc.csv", AICc)
        writedlm("$(LL_path)mcmc_samples.csv", mcmc_samples)
        writedlm("$(LL_path)var_covar_mcmc_samples.csv", cov(mcmc_samples))

        # This stores the results and signals of the 1st subject:
        # open("$(LL_path)step_by_step_info.csv", "w") do newfile
        # write(newfile, to_string(all_signals[1]))
        # end

        # Plot: re-call get_MAP()
        f = RL_Fit.plot_result(mcmc_samples, current_model.model_params,
            "$(current_model.model_name), $data_string, LL=$(round(LL_noPrior_MAP,digits=2)) | AIC=$(round(AIC,digits=2))",
            top_n_fraction = 1.)
        f.savefig("$(LL_path)MCMC.png")   # save the figure to file
        plt.close("all")

        f = RL_Fit.plot_result_logLL(mcmc_samples, current_model.model_params,
            "$(current_model.model_name), $data_string, LL=$(round(LL_noPrior_MAP,digits=2)) | AIC=$(round(AIC,digits=2))",
            top_n_fraction = 1.)
        f.savefig("$(LL_path)MCMC_LLdots.png")   # save the figure to file
        plt.close("all")

        f = RL_Fit.plot_result(mcmc_samples, current_model.model_params,
            "$(current_model.model_name), $data_string, LL=$(round(LL_noPrior_MAP,digits=2)) | AIC=$(round(AIC,digits=2))",
            top_n_fraction = 1.,
            hist_density = false)
        f.savefig("$(LL_path)MCMC_notnormed.png")   # save the figure to file
        plt.close("all")

        if fit_spec.do_per_subject_fit
            log_print_ln!(log_buffer, "Per subject fit:")

            all_parameters_subj = []
            params_with_pop_prior = []
            sum_pop_LL = 0.
            sum_subj_LL = 0.
            which_top = deepcopy(top_MAP)
            for i in 1:length(current_model.model_params)
                current_param = current_model.model_params[i]
                model_param_distr = current_param.prior_Distribution

                prior_mean = which_top[i]
                prior_std = current_model.per_param_posterior_std[i]
                prior_param_distr = TruncatedNormal(prior_mean, prior_std,
                                                    minimum(model_param_distr),
                                                    maximum(model_param_distr))

                model_param_prior = MetropolisParam(
                                    current_param.name,
                                    current_param.descr,
                                    prior_param_distr,
                                    current_param.rwalk_step_size_std)

                push!(params_with_pop_prior, model_param_prior)
            end

            for s_id in all_data_stats.all_subject_ids
              log_print!(log_buffer, "fitting subject_id $s_id. ")

              subj_signals, pop_signals = RL_Fit.run_single_subj_fit_write_csv(
                                        data_path,
                                        "$(LL_path)",
                                        s_id,
                                        fit_spec.fit_max_episodes,
                                        current_model.learner_function_V,
                                        params_with_pop_prior,
                                        which_top[1:end-1],
                                        fit_spec.nr_burnin,
                                        fit_spec.nr_samples,
                                        fit_spec.inter_collect_interval,
                                        fit_spec.nr_repetitions,
                                        fit_spec.top_n_fraction,
                                        description, rl_module_path)

              push!(all_parameters_subj,subj_signals.param_vector)
              log_print!(log_buffer, "subj_signals.param_vector= $(round.(subj_signals.param_vector,digits=3))")
              log_print_ln!(log_buffer, "subj_LL= $(round(subj_signals.LL,digits=3)), pop_LL= $(round(pop_signals.LL,digits=3))")
              sum_pop_LL += pop_signals.LL
              sum_subj_LL += subj_signals.LL
            end
            log_print_ln!(log_buffer, "sum_pop_LL= $(round(sum_pop_LL,digits=3)), sum_subj_LL=$(round(sum_subj_LL,digits=3))")
            # write results
            RL_Fit.writedlm("$(LL_path)all_subj_param_vectors.csv", all_parameters_subj)
            RL_Fit.writedlm("$(LL_path)sum_pop_param_LL.csv", sum_pop_LL)
            RL_Fit.writedlm("$(LL_path)sum_subj_param_LL.csv", sum_subj_LL)

        end
        log_print_ln!(log_buffer, "END fitting of model $(current_model.model_name)")
        log_print_ln!(log_buffer, "================================================================================")
    end
    log_print_ln!(log_buffer, "END: $(now())")
    # write log file
     log_file_name = "$(fit_spec.log_path)/Model_fit_log_$(Dates.format(now(), "dd_u_yyyy_HH_MM")).txt"
     log_str = String(take!(log_buffer))
     open(log_file_name, "w") do f
         write(f, log_str)
     end
end


text:
module splines

    export filter_coeffs, interpolant_cspline
    export eval_UC_spline, eval_UC_spline_G, eval_UC_multi_spline, eval_UC_multi_spline!
    include("csplines.jl")
    include("splines_filter.jl")

function interpolant_cspline(a, b, orders, V)

    coefs = filter_coeffs(a, b, orders, V)

    function fun(s::Array{Float64,2})
        return eval_UC_spline(a, b, orders, coefs, s)
    end

    function fun(p::Float64...)
        return fun([p...]')
    end

    return fun

end

end


text:
# Unless you have a reason not to precompile, you should set your package to precompile.

# Reasons to not precompile:
# 1) You have global pointers. This includes any globally defined Bigs! These
#    will cause segfaults.
# 2) You are calling `@eval` in a function that modifies your modules' global scope
#    This may have bad interactions with precompilation.
# 3) You are importing/using a package which has precompilation disabled

__precompile__()

# If you are not going to precompile your package, it's good practice to disable it:
# __precompile__(false)

module VariationalInference

# The main module file is for outlining the structure of the package.
# At the top of your package module, you should import your dependencies:

# using Compat

# Using Base modules.
using Random
# Load a plotting library.
using Plots
# Load the distributions library.
using Distributions

using StatsPlots

# Then define your abstract types:

# abstract AbstractMyType
# abstract AbstractMyType2 <: AbstractMyType
abstract type InferenceAlgorithm end

# Some people make sure all export statements are at the top of the main module as well
# Only values which are exported enter the namespace, anything else must be
# qualified by the package name (ExamplePackage.not_exported). By doing this all
# at the top file, it is easy to know what functionality is exported by a package

# Note: You should try to keep exports to a minimum. Exports are the public API:
# internal functionality should not be exported!

export  InferenceAlgorithm,
        fit,
        get_elbo

# Now include the real code
# Write your code in other files, otherwise the package outline gets muddled

include("CAVI.jl")
# include("output_string.jl")
# include("my_new_file.jl")

end # module


text:
"""
    MassInstallAction

Install or update GitHub Action workflows on repositories

API (all require qualification with `MassInstallAction`):

- Workflow creation: `Workflow`, `compat_helper`, `tag_bot`
- Workflow installation: `install`
"""
module MassInstallAction

include("types.jl")

include("public.jl")

include("default_workflows.jl")
include("repo_url.jl")
include("utils.jl")

end # module


text:
# This page was generated on DATEOFTODAY
using GeoPhyInv
using Test

# # Intro
#md # ```@docs
#md # AGeom
#md # ```

# # Examples

# Lets create a 2-D `mgrid` for the experiment.
mgrid2 = [range(0.0, stop = 10.0, step = 0.1), range(0.0, stop = 30.0, step = 0.2)];

# Then initialize an acquisition on `mgrid`, where the positions will be randomly chosen.
ageom2 = AGeom(mgrid2, SSrcs(10), Srcs(10), Recs(10));

# Similarly, we can do it for 3D too.
mgrid3 = fill(range(-10, stop=10, step=0.01),3);
ageom3 = AGeom(mgrid3, SSrcs(10), Srcs(10), Recs(10));

# You can check the dimensions.
GeoPhyInv.dim_names(ageom2), GeoPhyInv.dim_names(ageom3)

# For 2D, we can also initialize with one of the predefined acquisitions.
ageom2 = AGeom(mgrid2, :xwell, SSrcs(10), Recs(10));

# By default, the sources and receivers are randomly placed
# on the grid.
# If unsure, test it.
@test (ageom2 ∈ mgrid2)
@test (ageom3 ∈ mgrid3)

# The source and receiver positions can be updated as desired.
update!(ageom2[1], Srcs(), [0, 1], [10, 20]);
update!(ageom2[1], Recs(), [0, 0], [10, 20]);
update!(ageom2, SSrcs(), [0, 1], [10, 20]);

# It is easy to combine supersources. Now `ageom2` has 20 supersources.
ageom2_new = vcat(ageom2, ageom2);

# # Methods
#md # ```@docs
#md # update!(::AGeomss, ::Srcs)
#md # update!(::AGeomss, ::Recs)
#md # update!(::AGeom, ::SSrcs)
#md # update!(::AGeom, ::Recs)
#md # Base.in(::AGeom, ::AbstractVector{StepRangeLen})
#md # Base.isequal(::AGeom, ::AGeom)
#md # SparseArrays.SparseMatrixCSC(::AGeomss,::AbstractVector{StepRangeLen})
#md # ```


text:
# This file was generated, do not modify it. # hide
transformer_instance = Standardizer()
transformer_model = machine(transformer_instance, X)
fit!(transformer_model)
X = MLJ.transform(transformer_model, X);

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
using GeoRDD
using Test

function sim_data(n::Int, τ::Real)
    X = rand(2, n)
    f_X = sin.(X[1,:] * 2.0) .+ 3.0 .* X[2,:].^2.0
    m = -0.8 # mean offset
    treatment_radius = 0.8
    covarA = randn(n)
    covarB = randn(n)
    categ = rand(["A", "B", "C"], n) # random cateogry
    treat = sqrt.(X[1,:].^2 + X[2,:].^2) .< treatment_radius

    # outcome is: smooth surface + noise + treatment effect
    Y = m .+ # mean offset
        f_X .+ # smooth surface
        0.3*randn(n) .+ # noise
        τ .* treat .+ # treatment effect
        covarA .* 0.1 .+ # effect of real-valued covariate
        (categ.=="B") .* 0.2 # raise category B a little bit
    border_X1 = treatment_radius .* cos.(range(0,stop=π/2,length=1000))
    border_X2 = treatment_radius .* sin.(range(0,stop=π/2,length=1000))
    border_geo = LibGEOS.LineString([[border_X1[i], border_X2[i]] for i in 1:n])
    geordd_df = DataFrame(
        X1 = X[1, :],
        X2 = X[2, :],
        outcome = Y,
        region = treat, # treatment indicator
        covarA = covarA,
        covarB = covarB,
        categ = categ
        )
    categorical!(geordd_df, :categ)
    categorical!(geordd_df, :region)
    return geordd_df, border_X1, border_X2, border_geo
end
include("test_geometry.jl")
include("test_gprealisations.jl")
include("test_simulated.jl")


text:
"""
"""
module Compilation # begin submodule PredictMD.Compilation

############################################################################
# PredictMD.Compilation source files #######################################
############################################################################

# submodules/Server/buildexecutable/
include("buildexecutable/buildexecutable.jl")

# submodules/Server/packagecompiler/
include("packagecompiler/packagecompiler.jl")

end # end submodule PredictMD.Compilation


text:
"""
    common.jl

Description:
    Includes all of the unsupervised ARTMAP modules common code.
"""

"""
    train!(art::ARTMAP, x::RealMatrix, y::IntegerVector, preprocessed::Bool=false)

Train the ARTMAP model on a batch of data 'x' with supervisory labels 'y.'

# Arguments
- `art::ARTMAP`: the supervised ARTMAP model to train.
- `x::RealMatrix`: the 2-D dataset containing columns of samples with rows of features.
- `y::IntegerVector`: labels for supervisory training.
- `preprocessed::Bool=false`: flag, if the data has already been complement coded or not.
"""
function train!(art::ARTMAP, x::RealMatrix, y::IntegerVector, preprocessed::Bool=false)
    # Show a message if display is on
    art.opts.display && @info "Training $(typeof(art))"

    # Data information and setup
    n_samples = length(y)

    # Run the batch initialization procedure
    x = init_train!(x, art, preprocessed)

    # Initialize the output vector
    y_hat = zeros(Int, n_samples)
    # Learn until the stopping conditions
    art.epoch = 0
    while true
        # Increment the epoch and get the iterator
        art.epoch += 1
        iter = get_iterator(art.opts, x)
        for i = iter
            # Update the iterator if necessary
            update_iter(art, iter, i)
            # Grab the sample slice
            # sample = get_sample(x, i)
            sample = x[:, i]
            label = y[i]
            # Train upon the sample and label
            y_hat[i] = train!(art, sample, label, preprocessed=true)
        end

        # Check stopping conditions
        if stopping_conditions(art)
            break
        end
    end
    return y_hat
end # train!(art::ARTMAP, x::RealMatrix, y::IntegerVector, preprocessed::Bool=false)

# """
#     classify(art::ARTMAP, x::RealMatrix ; preprocessed::Bool=false, get_bmu::Bool=false)

# Predict categories of 'x' using the ARTMAP model.

# Returns predicted categories 'y_hat.'

# # Arguments
# - `art::ARTMAP`: supervised ARTMAP module to use for batch inference.
# - `x::RealMatrix`: the 2-D dataset containing columns of samples with rows of features.
# - `preprocessed::Bool=false`: flag, if the data has already been complement coded or not.
# - `get_bmu::Bool=false`, flag, if the model should return the best-matching-unit label in the case of total mismatch.

# # Examples
# ```julia-repl
# julia> my_DDVFA = DDVFA()
# DDVFA
#     opts: opts_DDVFA
#     ...
# julia> x, y = load_data()
# julia> train!(my_DDVFA, x)
# julia> y_hat = classify(my_DDVFA, y)
# ```
# """
# function classify(art::ARTMAP, x::RealMatrix ; preprocessed::Bool=false, get_bmu::Bool=false)
#     # Show a message if display is on
#     art.opts.display && @info "Testing $(typeof(art))"

#     # Preprocess the data
#     x = init_classify!(x, art, preprocessed)

#     # Data information and setup
#     n_samples = get_n_samples(x)

#     # Initialize the output vector
#     y_hat = zeros(Int, n_samples)

#     # Get the iterator based on the module options and data shape
#     iter = get_iterator(art.opts, x)
#     for ix = iter
#         # Update the iterator if necessary
#         update_iter(art, iter, ix)

#         # Grab the sample slice
#         sample = get_sample(x, ix)

#         # Get the classification
#         y_hat[ix] = classify(art, sample, preprocessed=true, get_bmu=get_bmu)
#     end

#     return y_hat
# end # classify(art::ARTMAP, x::RealMatrix ; preprocessed::Bool=false, get_bmu::Bool=false)


text:
function memorygame(starting, place=2020)
    seen = Dict{Int,Int}()
    for (i, x) in enumerate(starting[1:end-1])
        seen[x] = i
    end
    count = length(starting)
    previous = starting[end]
    while count < place
        if haskey(seen, previous)
            next = count - seen[previous]
            seen[previous] = count
            previous = next
        else
            seen[previous] = count
            previous = 0
        end
        count += 1
    end
    previous
end

function main()
    starting = [9, 12, 1, 4, 17, 0, 18]
    part1 = memorygame(starting, 2020)
    part2 = memorygame(starting, 30000000)
    @info "Solution" part1 part2
end

main()

max_stars_repo_name: dglmoore/advent-of-code

max_issues_repo_name: dglmoore/advent-of-code

max_forks_repo_name: dglmoore/advent-of-code


text:
using Distributed
@everywhere using MCPhylo

## Data
epil = Dict{Symbol, Any}(
  :y =>
    [ 5  3  3  3
      3  5  3  3
      2  4  0  5
      4  4  1  4
      7 18  9 21
      5  2  8  7
      6  4  0  2
     40 20 21 12
      5  6  6  5
     14 13  6  0
     26 12  6 22
     12  6  8  4
      4  4  6  2
      7  9 12 14
     16 24 10  9
     11  0  0  5
      0  0  3  3
     37 29 28 29
      3  5  2  5
      3  0  6  7
      3  4  3  4
      3  4  3  4
      2  3  3  5
      8 12  2  8
     18 24 76 25
      2  1  2  1
      3  1  4  2
     13 15 13 12
     11 14  9  8
      8  7  9  4
      0  4  3  0
      3  6  1  3
      2  6  7  4
      4  3  1  3
     22 17 19 16
      5  4  7  4
      2  4  0  4
      3  7  7  7
      4 18  2  5
      2  1  1  0
      0  2  4  0
      5  4  0  3
     11 14 25 15
     10  5  3  8
     19  7  6  7
      1  1  2  3
      6 10  8  8
      2  1  0  0
    102 65 72 63
      4  3  2  4
      8  6  5  7
      1  3  1  5
     18 11 28 13
      6  3  4  0
      3  5  4  3
      1 23 19  8
      2  3  0  1
      0  0  0  0
      1  4  3  2],
  :Trt =>
    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
     0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
     1, 1, 1, 1, 1, 1, 1, 1, 1],
  :Base =>
    [11, 11, 6, 8, 66, 27, 12, 52, 23, 10, 52, 33, 18, 42, 87, 50, 18, 111, 18,
     20, 12, 9, 17, 28, 55, 9, 10, 47, 76, 38, 19, 10, 19, 24, 31, 14, 11, 67,
     41, 7, 22, 13, 46, 36, 38, 7, 36, 11, 151, 22, 41, 32, 56, 24, 16, 22, 25,
     13, 12],
  :Age =>
    [31, 30, 25, 36, 22, 29, 31, 42, 37, 28, 36, 24, 23, 36, 26, 26, 28, 31, 32,
     21, 29, 21, 32, 25, 30, 40, 19, 22, 18, 32, 20, 30, 18, 24, 30, 35, 27, 20,
     22, 28, 23, 40, 33, 21, 35, 25, 26, 25, 22, 32, 25, 35, 21, 41, 32, 26, 21,
     36, 37],
  :V4 => [0, 0, 0, 1]
)
epil[:N] = size(epil[:y], 1)
epil[:T] = size(epil[:y], 2)

epil[:logBase4] = log.(epil[:Base] / 4)
epil[:BT] = epil[:logBase4] .* epil[:Trt]
epil[:logAge] = log.(epil[:Age])
map(key -> epil[Symbol(string(key, "bar"))] = mean(epil[key]),
    [:logBase4, :Trt, :BT, :logAge, :V4])

## Model Specification
model = Model(

  y = Stochastic(2,
    (a0, alpha_Base, logBase4, logBase4bar, alpha_Trt, Trt, Trtbar,
     alpha_BT, BT, BTbar, alpha_Age, logAge, logAgebar, alpha_V4, V4,
     V4bar, b1, b, N, T) ->
      UnivariateDistribution[(
        mu = exp(a0 + alpha_Base * (logBase4[i] - logBase4bar) +
                 alpha_Trt * (Trt[i] - Trtbar) + alpha_BT * (BT[i] - BTbar) +
                 alpha_Age * (logAge[i] - logAgebar) +
                 alpha_V4 * (V4[j] - V4bar) + b1[i] +
                 b[i, j]);
        Poisson(mu)) for i in 1:N, j in 1:T
      ],
    false
  ),

  b1 = Stochastic(1,
    s2_b1 -> Normal(0, sqrt(s2_b1)),
    false
  ),

  b = Stochastic(2,
    s2_b -> Normal(0, sqrt(s2_b)),
    false
  ),

  a0 = Stochastic(
    () -> Normal(0, 100),
    false
  ),

  alpha_Base = Stochastic(
    () -> Normal(0, 100)
  ),

  alpha_Trt = Stochastic(
    () -> Normal(0, 100)
  ),

  alpha_BT = Stochastic(
    () -> Normal(0, 100)
  ),

  alpha_Age = Stochastic(
    () -> Normal(0, 100)
  ),

  alpha_V4 = Stochastic(
    () -> Normal(0, 100)
  ),

  alpha0 = Logical(
    (a0, alpha_Base, logBase4bar, alpha_Trt, Trtbar, alpha_BT, BTbar,
     alpha_Age, logAgebar, alpha_V4, V4bar) ->
      a0 - alpha_Base * logBase4bar - alpha_Trt * Trtbar - alpha_BT * BTbar -
        alpha_Age * logAgebar - alpha_V4 * V4bar
  ),

  s2_b1 = Stochastic(
    () -> InverseGamma(0.001, 0.001)
  ),

  s2_b = Stochastic(
    () -> InverseGamma(0.001, 0.001)
  )

)

## Initial Values
inits = [
  Dict(:y => epil[:y], :a0 => 0, :alpha_Base => 0, :alpha_Trt => 0,
       :alpha_BT => 0, :alpha_Age => 0, :alpha_V4 => 0, :s2_b1 => 1,
       :s2_b => 1, :b1 => zeros(epil[:N]), :b => zeros(epil[:N], epil[:T])),
  Dict(:y => epil[:y], :a0 => 1, :alpha_Base => 1, :alpha_Trt => 1,
       :alpha_BT => 1, :alpha_Age => 1, :alpha_V4 => 1, :s2_b1 => 10,
       :s2_b => 10, :b1 => zeros(epil[:N]), :b => zeros(epil[:N], epil[:T]))
]

## Sampling Scheme
scheme = [AMWG([:a0, :alpha_Base, :alpha_Trt, :alpha_BT, :alpha_Age,
                :alpha_V4], 0.1),
          Slice(:b1, 0.5),
          Slice(:b, 0.5),
          Slice([:s2_b1, :s2_b], 1.0)]
setsamplers!(model, scheme)

## MCMC Simulations
sim = mcmc(model, epil, inits, 15000, burnin=2500, thin=2, chains=2)
describe(sim)


text:
using MomentTensors, Test

@testset "All tests" begin
    include("construction.jl")
    include("indexing.jl")
    include("arithmetic.jl")
    include("conversion.jl")
    include("planes.jl")
    include("io.jl")
    include("decompose.jl")
    include("units.jl")
end


text:
# this user supplied file creates a dictionary of arguments
# if a key is repeated, the last use of the key is used
# it is a little bit dangerous letting the user run arbitrary code
# as part of the solver
# now that this file is read inside a function, it is better encapsulated

arg_dict = Dict{String, Any}(
"physics" => "Euler",
"var1" => 1,
"var2" => "a",
"var3" => 3.5,
"var4" => [1,2,3],
"var3" => 4,
"run_type" => 5,
"jac_method" => 2,
"jac_type" => 2,
"order" => 3,
"IC_name" => "ICIsentropicVortex",
#"IC_name" => "ICFile",
#"ICfname" => "start.dat",
"Relfunc_name" => "ICRho1E2U3",
"variable_type" => :entropy,
"numBC" => 1,
"BC1" => [ 0, 1, 2, 3],
"BC1_name" => "isentropicVortexBC",
#"BC2" => [4, 10],
#"BC2_name" => "noPenetrationBC",
#"BC2_name" => "isentropicVortexBC",
"delta_t" => 0.005,
"t_max" => 500.000,
"smb_name" => "SRCMESHES/squarevortex_large.smb",
"dmg_name" => ".null",
"res_abstol" => -1e-14,
"res_reltol" => 1e-10,
"step_tol" => 1e-10,
"itermax" => 0,
"use_edgestab" => true,
"edgestab_gamma" => -0.01,
"use_dissipation" => false,
"dissipation_name" => "damp1",
"dissipation_const" => 3.00,
"calc_error" => true,
"calc_error_infname" => "solution_final.dat",
#"writeq" => true,
#"perturb_ic" => true,
#"perturb_mag" => 0.001,
#"write_sparsity" => true,
#"write_jac" => true,
#"write_edge_vertnums" => true,
#"write_face_vertnums" => true,
#"write_qic" => true,
#"writeboundary" => true,
#"write_res" => true,
#"write_finalsolution" => true,
#"write_finalresidual" => true,
"write_counts" => true,
"write_rhs" => false,
"solve" => false,
)


text:
function update_y_imputed!(i::Int, n::Int, j::Int, s::State, c::Constants,
                          d::Data, tuners::Tuners)

  function logFullCond(yinj)
    p = prob_miss(yinj, c.beta[:, i])
    k = s.lam[i][n]
    if k > 0
      z = s.Z[j, k]
      l = s.gam[i][n, j]
      mu = mus(z, l, s, c, d)
      sig = sqrt(s.sig2[i])
      logPrior = logpdf(Normal(mu, sig), yinj)
    else
      logPrior = logpdf(c.noisyDist, yinj)
    end
    return log(p) + logPrior
  end

  s.y_imputed[i][n, j] = MCMC.metropolisAdaptive(s.y_imputed[i][n, j],
                                                 logFullCond,
                                                 tuners.y_imputed[i, n, j])
end

function update_y_imputed!(s::State, c::Constants, d::Data, tuners::Tuners)
  for i in 1:d.I
    for j in 1:d.J
      for n in 1:d.N[i]
        if d.m[i][n, j] == 1
          update_y_imputed!(i, n, j, s, c, d, tuners)
        end
      end
    end
  end
end

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch

max_forks_repo_name: luiarthur/CytofResearch


text:
#mutable struct DecisionVariableRef <: RobustVariableRef
#end

struct Decision <: RobustVariableType
    stage::Int
    policy::String
    depends_on::Dict{Int64, AbstractVariable}
end

struct DecisionVariable<: RobustVariable
    info::JuMP.VariableInfo
    stage::Int
    policy::String
    depends_on::Dict{Int64, AbstractVariable}
end
export Decision, DecisionVariable
#Decision abhängig von Uncertain
# policy und depends on
# policy ist ein keyord wie affine was dann nachher einen switch
# depends on  nur im DECISIOSN= Vektor mit UNcerstain variable die definiert sein müssen


text:
write_test_output = false

path = dirname(@__FILE__)

m = AnSchorfheide()

save = normpath(joinpath(dirname(@__FILE__),"save"))
m <= Setting(:saveroot, saveroot)

data = h5read("reference/smc.h5", "data")

m <= Setting(:n_particles, 400)
m <= Setting(:n_Φ, 100)
m <= Setting(:λ, 2.0)
m <= Setting(:n_smc_blocks, 1)
m <= Setting(:use_parallel_workers, false)
m <= Setting(:step_size_smc, 0.5)
m <= Setting(:n_mh_steps_smc, 1)
m <= Setting(:resampler_smc, :multinomial)
m <= Setting(:target_accept, 0.25)

m <= Setting(:mixture_proportion, .9)
m <= Setting(:tempering_target, 0.95)
m <= Setting(:resampling_threshold, .5)
m <= Setting(:use_fixed_schedule, true)

n_parts = get_setting(m, :n_particles)

file = JLD2.jldopen("reference/mutation_inputs.jld2", "r")
old_particles = read(file, "particles")
d = read(file, "d")
blocks_free = read(file, "blocks_free")
blocks_all = read(file, "blocks_all")
ϕ_n = read(file, "ϕ_n")
ϕ_n1 = read(file, "ϕ_n1")
c = read(file, "c")
α = read(file, "α")
old_data = read(file, "old_data")
close(file)

old_part_cloud = DSGE.vector_particles_to_cloud(m, old_particles)
function my_likelihood(parameters::ParameterVector, data::Matrix{Float64})::Float64
    DSGE.update!(m, parameters)
    DSGE.likelihood(m, data; sampler = false, catch_errors = true,
                    use_chand_recursion = true, verbose = :low)
end

Random.seed!(42)
new_particles = [SMC.mutation(my_likelihood, m.parameters, data,
                              old_part_cloud.particles[j, :], d.μ, Matrix(d.Σ),
                              blocks_free, blocks_all, ϕ_n, ϕ_n1;
                              c = c, α = α, old_data = old_data) for j = 1:n_parts]

if write_test_output
    JLD2.jldopen("reference/mutation_outputs.jld2", "w") do file
        write(file, "particles", new_particles)
    end
end

saved_particles = load("reference/mutation_outputs.jld2", "particles")

@testset "Test mutation outputs, particle by particle" begin
    for i = 1:length(saved_particles)
        @test isapprox(saved_particles[i], new_particles[i], nans = true)
    end
end


text:
module Rounding
include("fenv_constants.jl")

export
    RoundingMode, RoundNearest, RoundToZero, RoundUp, RoundDown, RoundFromZero,
    get_rounding, set_rounding, with_rounding

## rounding modes ##
immutable RoundingMode
    code::Int
    RoundingMode(c::Integer) = new(c)
end
const RoundNearest = RoundingMode(0)
const RoundToZero = RoundingMode(1)
const RoundUp = RoundingMode(2)
const RoundDown = RoundingMode(3)
const RoundFromZero = RoundingMode(4)

function to_fenv(r::RoundingMode)
    if r === RoundNearest
        JL_FE_TONEAREST
    elseif r === RoundToZero
        JL_FE_TOWARDZERO
    elseif r === RoundUp
        JL_FE_UPWARD
    elseif r === RoundDown
        JL_FE_DOWNWARD
    elseif r === RoundFromZero
        error("unsupported rounding mode")
    else
        error("invalid rounding mode")
    end
end

function from_fenv(r::Integer)
    if r == JL_FE_TONEAREST
        return RoundNearest
    elseif r == JL_FE_DOWNWARD
        return RoundDown
    elseif r == JL_FE_UPWARD
        return RoundUp
    elseif r == JL_FE_TOWARDZERO
        return RoundToZero
    else
        error("invalid rounding mode code")
    end
end

set_rounding{T<:Union(Float32,Float64)}(::Type{T},r::RoundingMode) = ccall(:fesetround, Cint, (Cint,), to_fenv(r))
get_rounding{T<:Union(Float32,Float64)}(::Type{T}) = from_fenv(ccall(:fegetround, Cint, ()))

function with_rounding{T}(f::Function, ::Type{T}, rounding::RoundingMode)
    old_rounding = get_rounding(T)
    set_rounding(T,rounding)
    try
        return f()
    finally
        set_rounding(T,old_rounding)
    end
end

end #module

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
#!/usr/bin/env julia

# Copyright (c) 2016, Jost von Hardenberg - ISAC-CNR, Italy

using RainFARM
using ArgParse
using Printf

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--nf", "-n"
            help = "Subdivisions for downscaling [default: 2]"
            arg_type = Int
            default = 2
        "--outfile", "-o", "--out"
            help = "Weights file"
            arg_type = AbstractString
            default = "grid.nc"
        "infile"
            help = "The input file to downscale"
            arg_type = AbstractString
            required = true
        "--varname", "-v"
            help = "Input variable name"
            arg_type = AbstractString
            default = ""
    end
	s.description = "Generates a sample file with a grid equal to the input file downscaled by a factor NF"
    return parse_args(s)
end

args = parse_commandline()
nf = args["nf"]
filenc = args["infile"]
outfile = args["outfile"]
varnc = args["varname"]

println("Reading file ", filenc)

(pr, lon_mat, lat_mat, varnc) = read_netcdf2d(filenc, varnc)

# Create the fine grid
nss = size(pr)
if length(nss)>=3
    pr = pr[:,:,1]
end
ns = nss[1]

(lon_f, lat_f) = lon_lat_fine(lon_mat, lat_mat,nf)

println("Output size: ", size(lon_f))

write_netcdf2d(outfile, reshape(pr, ns ,ns ,1), lon_f, lat_f, varnc, filenc)


text:
struct SingletonDict{K,V} <: AbstractDict{K,V}
    key::K
    value::V
end
SingletonDict(p::Pair{K,V}) where {K,V} = SingletonDict{K,V}(p.first, p.second)

Base.length(::SingletonDict) = 1
Base.keys(d::SingletonDict) = (d.key,)
Base.values(d::SingletonDict) = (d.value,)
Base.haskey(d::SingletonDict, key) = isequal(d.key, key)
Base.getindex(d::SingletonDict, key) = isequal(d.key, key) ? d.value : throw(KeyError(key))
Base.get(d::SingletonDict, key, default) = isequal(d.key, key) ? d.value : default

Base.iterate(d::SingletonDict, s = true) = s ? ((d.key => d.value), false) : nothing

struct VectorDict{K,V} <: AbstractDict{K,V}
    keys::Vector{K}
    values::Vector{V}
end
VectorDict{K,V}() where {K,V} = VectorDict{K,V}(Vector{K}(), Vector{V}())
function VectorDict{K,V}(kv) where {K,V}
    keys = Vector{K}()
    values = Vector{V}()
    if Base.IteratorSize(kv) !== SizeUnknown()
        sizehint!(keys, length(kv))
        sizehint!(values, length(kv))
    end
    for (k,v) in kv
        push!(keys, k)
        push!(values, v)
    end
    return VectorDict{K,V}(keys, values)
end
VectorDict(kv::Pair{K,V}...) where {K,V} = VectorDict{K,V}(kv)
VectorDict(g::Base.Generator) = VectorDict(g...)

Base.length(d::VectorDict) = length(d.keys)
Base.sizehint!(d::VectorDict, newsz) = (sizehint!(d.keys, newsz); sizehint!(d.values, newsz); return d)

@propagate_inbounds getpair(d::VectorDict, i::Integer) = d.keys[i] => d.values[i]

Base.copy(d::VectorDict) = VectorDict(copy(d.keys), copy(d.values))
Base.empty(::VectorDict, ::Type{K}, ::Type{V}) where {K, V} = VectorDict{K, V}()
Base.empty!(d::VectorDict) = (empty!(d.keys); empty!(d.values); return d)

function Base.delete!(d::VectorDict, key)
    i = findfirst(isequal(key), d.keys)
    if !(i == nothing || i == 0)
        deleteat!(d.keys  , i)
        deleteat!(d.values, i)
    end
    return d
end

Base.keys(d::VectorDict) = d.keys
Base.values(d::VectorDict) = d.values
Base.haskey(d::VectorDict, key) = key in d.keys
function Base.getindex(d::VectorDict, key)
    i = findfirst(isequal(key), d.keys)
    @inbounds begin
        return i !== nothing ? d.values[i] : throw(KeyError(key))
    end
end
function Base.setindex!(d::VectorDict, v, key)
    i = findfirst(isequal(key), d.keys)
    if i === nothing
        push!(d.keys, key)
        push!(d.values, v)
    else
        d.values[i] = v
    end
    return d
end

function Base.get(d::VectorDict, key, default)
    i = findfirst(isequal(key), d.keys)
    @inbounds begin
        return i !== nothing ? d.values[i] : default
    end
end

function Base.iterate(d::VectorDict, s = 1)
    @inbounds if s > length(d)
        return nothing
    else
        return (d.keys[s] => d.values[s]), s+1
    end
end

struct SortedVectorDict{K,V} <: AbstractDict{K,V}
    keys::Vector{K}
    values::Vector{V}
    function SortedVectorDict{K,V}(pairs::Vector{Pair{K,V}}) where {K,V}
        if !issorted(pairs, by=first)
            pairs = sort(pairs, by=first)
        end
        return new{K,V}(map(first, pairs), map(last, pairs))
    end
    SortedVectorDict{K,V}(keys::Vector{K}, values::Vector{V}) where {K,V} =
        new{K,V}(keys, values)
    SortedVectorDict{K,V}() where {K,V} = new{K,V}(Vector{K}(undef, 0), Vector{V}(undef, 0))
end
function SortedVectorDict{K,V}(kv) where {K,V}
    d = SortedVectorDict{K,V}()
    if Base.IteratorSize(kv) !== SizeUnknown()
        sizehint!(d, length(kv))
    end
    for (k,v) in kv
        push!(d, k=>v)
    end
    return d
end
SortedVectorDict(pairs::Vector{Pair{K,V}}) where {K,V} = SortedVectorDict{K,V}(pairs)
SortedVectorDict(kv::Pair{K,V}...) where {K,V} = SortedVectorDict{K,V}(kv)
SortedVectorDict(g::Base.Generator) = SortedVectorDict(g...)

Base.length(d::SortedVectorDict) = length(d.keys)
Base.sizehint!(d::SortedVectorDict, newsz) =
    (sizehint!(d.keys, newsz); sizehint!(d.values, newsz); return d)

Base.copy(d::SortedVectorDict{K,V}) where {K,V} = SortedVectorDict{K,V}(copy(d.keys), copy(d.values))
Base.empty(::SortedVectorDict, ::Type{K}, ::Type{V}) where {K, V} = SortedVectorDict{K, V}()
Base.empty!(d::SortedVectorDict) = (empty!(d.keys); empty!(d.values); return d)

# _searchsortedfirst(v::Vector, k) = searchsortedfirst(v, k)
function _searchsortedfirst(v::Vector, k)
    i = 1
    @inbounds while i <= length(v) && isless(v[i], k)
        i += 1
    end
    return i
end

function Base.delete!(d::SortedVectorDict{K}, k) where {K}
    key = convert(K, k)
    if !isequal(k, key)
        return d
    end
    i = _searchsortedfirst(d.keys, key)
    if i <= length(d) && isequal(d.keys[i], key)
        deleteat!(d.keys, i)
        deleteat!(d.values, i)
    end
    return d
end

Base.keys(d::SortedVectorDict) = d.keys
Base.values(d::SortedVectorDict) = d.values
function Base.haskey(d::SortedVectorDict{K}, k) where {K}
    key = convert(K, k)
    if !isequal(k, key)
        return false
    end
    i = _searchsortedfirst(d.keys, key)
    return (i <= length(d) && isequal(d.keys[i], key))
end
function Base.getindex(d::SortedVectorDict{K}, k) where {K}
    key = convert(K, k)
    if !isequal(k, key)
        throw(KeyError(k))
    end
    i = _searchsortedfirst(d.keys, key)
    @inbounds begin
        return (i <= length(d) && isequal(d.keys[i], key)) ? d.values[i] : throw(KeyError(key))
    end
end
function Base.setindex!(d::SortedVectorDict{K}, v, k) where {K}
    key = convert(K, k)
    if !isequal(k, key)
        throw(ArgumentError("$k is not a valid key for type $K"))
    end
    i = _searchsortedfirst(d.keys, key)
    if i <= length(d) && isequal(d.keys[i], key)
        d.values[i] = v
    else
        insert!(d.keys, i, key)
        insert!(d.values, i, v)
    end
    return d
end

function Base.get(d::SortedVectorDict{K}, k, default) where {K}
    key = convert(K, k)
    if !isequal(k, key)
        return default
    end
    i = _searchsortedfirst(d.keys, key)
    @inbounds begin
        return (i <= length(d) && isequal(d.keys[i], key)) ? d.values[i] : default
    end
end
function Base.get(f::Union{Function,Type}, d::SortedVectorDict{K}, k) where {K}
    key = convert(K, k)
    if !isequal(k, key)
        return f()
    end
    i = _searchsortedfirst(d.keys, key)
    @inbounds begin
        return (i <= length(d) && isequal(d.keys[i], key)) ? d.values[i] : f()
    end
end
function Base.iterate(d::SortedVectorDict, i = 1)
    @inbounds if i > length(d)
        return nothing
    else
        return (d.keys[i] => d.values[i]), i+1
    end
end


text:
rig_key = "ocpi-2"
push!(RIGS, rig_key)

const ocpi2_mappings = Dict("AO0"=>"axial piezo",
                      "AO1"=>"horizontal piezo",
                      "AO2"=>"galvo1",
                      "AO3"=>"galvo2",
                      "AI0"=>"axial piezo monitor",
                      "AI1"=>"horizontal piezo monitor",
                      "AI2"=>"galvo1 monitor",
                      "AI3"=>"galvo2 monitor",
                      "AI4"=>"analogin3",
                      "AI5"=>"analogin4",
                      "AI6"=>"analogin5",
                      "AI7"=>"analogin6",
                      "AI8"=>"analogin7",
                      "AI9"=>"analogin8",
                      "AI10"=>"analogin9",
                      "AI11"=>"analogin10",
                      "AI12"=>"analogin11",
                      "AI13"=>"analogin12",
                      "AI14"=>"analogin13",
                      "AI15"=>"analogin14",
                      "AI16"=>"analogin15",
                      "AI17"=>"analogin16",
                      "AI18"=>"analogin17",
                      "AI19"=>"analogin18",
                      "AI20"=>"analogin19",
                      "AI21"=>"analogin20",
                      "AI22"=>"analogin21",
                      "AI23"=>"analogin22",
                      "AI24"=>"analogin23",
                      "AI25"=>"analogin24",
                      "AI26"=>"analogin25",
                      "AI27"=>"analogin26",
                      "AI28"=>"analogin27",
                      "AI29"=>"analogin28",
                      "AI30"=>"analogin29",
                      "AI31"=>"analogin30",
                      "P0.0"=>"stimulus1",
                      "P0.1"=>"stimulus2",
                      "P0.2"=>"stimulus3",
                      "P0.3"=>"stimulus4",
                      "P0.4"=>"all lasers",
                      "P0.5"=>"camera1",
                      "P0.6"=>"camera2",
#                      "P0.7"=>"reserved",
                      "P0.8"=>"405nm laser",
                      "P0.9"=>"445nm laser",
                      "P0.10"=>"488nm laser",
                      "P0.11"=>"514nm laser",
                      "P0.12"=>"561nm laser",
                      "P0.13"=>"stimulus5",
                      "P0.14"=>"stimulus6",
                      "P0.15"=>"stimulus7",
                      "P0.16"=>"stimulus8",
                      "P0.17"=>"stimulus9",
                      "P0.18"=>"stimulus10",
                      "P0.19"=>"stimulus11",
                      "P0.20"=>"stimulus12",
                      "P0.21"=>"stimulus13",
                      "P0.22"=>"stimulus14",
                      "P0.23"=>"stimulus15",
                      "P0.24"=>"camera1 frame monitor",
                      "P0.25"=>"camera2 frame monitor",
                      "P0.26"=>"diginput1",
                      "P0.27"=>"diginput2",
                      "P0.28"=>"diginput3",
                      "P0.29"=>"diginput4",
                      "P0.30"=>"diginput5",
                      "P0.31"=>"diginput6",)

DEFAULT_DAQCHANS_TO_NAMES[rig_key] = ocpi2_mappings
DEFAULT_NAMES_TO_DAQCHANS[rig_key] = dictmap(reverse, ocpi2_mappings)

const ocpi2_aochans = map(x->"AO$(x)", 0:3)
const ocpi2_aichans = map(x->"AI$(x)", 0:31)
const ocpi2_dochans = map(x->"P0.$(x)", vcat([0:6...], [8:23...]))
const ocpi2_dichans = map(x->"P0.$(x)", 24:31)
const ocpi2_pos_ctrl_chans = map(x->DEFAULT_NAMES_TO_DAQCHANS[rig_key][x], ["axial piezo"; "horizontal piezo"])
const ocpi2_pos_mon_chans = map(x->DEFAULT_NAMES_TO_DAQCHANS[rig_key][x], ["axial piezo monitor"; "horizontal piezo monitor"])
const ocpi2_cam_ctrl_chans = map(x->DEFAULT_NAMES_TO_DAQCHANS[rig_key][x], ["camera1"; "camera2"])
const ocpi2_cam_mon_chans = map(x->DEFAULT_NAMES_TO_DAQCHANS[rig_key][x], ["camera1 frame monitor"; "camera2 frame monitor"])
const ocpi2_laschans = map(x->DEFAULT_NAMES_TO_DAQCHANS[rig_key][x], ["405nm laser"; "445nm laser"; "488nm laser"; "514nm laser"; "561nm laser"; "all lasers"])
const ocpi2_stimchans = map(x->DEFAULT_NAMES_TO_DAQCHANS[rig_key][x], ["stimulus$x" for x = 1:15])
const ocpi2_galvo_ctrl_chans = map(x->DEFAULT_NAMES_TO_DAQCHANS[rig_key][x], ["galvo1"; "galvo2"])
const ocpi2_galvo_mon_chans = map(x->DEFAULT_NAMES_TO_DAQCHANS[rig_key][x], ["galvo1 monitor"; "galvo2 monitor"])
const ocpi2_fixed_names = ["axial piezo", "axial piezo monitor", "horizontal piezo", "horizontal piezo monitor", "camera1", "camera1 frame monitor", "camera2", "camera2 frame monitor", "405nm laser", "445nm laser", "488nm laser", "514nm laser", "561nm laser", "galvo1", "galvo2", "galvo1 monitor", "galvo2 monitor"]
const ocpi2_fixed_daqchans = map(x->DEFAULT_NAMES_TO_DAQCHANS[rig_key][x], ocpi2_fixed_names)

AO_CHANS[rig_key] = OrderedSet(ocpi2_aochans)
AI_CHANS[rig_key] = OrderedSet(ocpi2_aichans)
DO_CHANS[rig_key] = OrderedSet(ocpi2_dochans)
DI_CHANS[rig_key] = OrderedSet(ocpi2_dichans)
POS_CONTROL_CHANS[rig_key] = ocpi2_pos_ctrl_chans
POS_MONITOR_CHANS[rig_key] = ocpi2_pos_mon_chans
CAM_CONTROL_CHANS[rig_key] = ocpi2_cam_ctrl_chans
CAM_MONITOR_CHANS[rig_key] = ocpi2_cam_mon_chans
LAS_CONTROL_CHANS[rig_key] = ocpi2_laschans
STIM_CHANS[rig_key] =OrderedSet(ocpi2_stimchans)
GALVO_CONTROL_CHANS[rig_key] = ocpi2_galvo_ctrl_chans
GALVO_MONITOR_CHANS[rig_key] = ocpi2_galvo_mon_chans
FIXED_NAMES[rig_key] = ocpi2_fixed_names
FIXED_DAQ_CHANS[rig_key] = ocpi2_fixed_daqchans

RIG_CHIP_SIZES[rig_key] = PCO_EDGE_4_2_CHIP_SIZE
RIG_FRAMERATE_FUNCS[rig_key] = PCO_EDGE_4_2_FRAMERATE_FUNC
RIG_ROI_FUNCS[rig_key] = PCO_EDGE_4_2_ROI_FUNC

#TODO: measure these
LASER_ON_TIME[rig_key] = 9 * Unitful.μs
LASER_OFF_TIME[rig_key] = 9 * Unitful.μs

CAMERA_ON_TIME[rig_key] = 19526.0 * Unitful.ns #This is the worst jitter measured with Edge 5.5 and 4.2 cameras, see ImagineInterface issue #18
CAMERA_OFF_TIME[rig_key] = 19526.0 * Unitful.ns #This is the worst jitter measured with Edge 5.5 and 4.2 cameras, see ImagineInterface issue #18

PIEZO_RANGES[rig_key] = (0.0μm .. 800.0μm, 0.0V .. 10.0V)
PIEZO_MAX_SPEED[rig_key] = 2000*Unitful.μm / Unitful.s
AO_RANGE[rig_key] = -10.0V .. 10.0V
AI_RANGE[rig_key] = AO_RANGE[rig_key] #TODO: make sure this is true.  (true if we are recording -10..10V on analog inputs)

max_forks_repo_name: HolyLab/ImagineInterface


text:
#########################
# Guido Dipietro - 2020 #
#########################

using BenchmarkTools

cp = (1,2,3,4,5,6,7,8) # like,speffz
co = (0,0,0,0,0,0,0,0) # 1 = clockwise
cube = (cp, co)

# Atoms
Rp((a,b,c,d,e,f,g,h)) = (a,c,f,d,e,g,b,h)
Ro(a) = @. (a + (0,1,2,0,0,1,2,0)) % 3
Up((a,b,c,d,e,f,g,h)) = (d,a,b,c,e,f,g,h)
xp((a,b,c,d,e,f,g,h)) = (d,c,f,e,h,g,b,a)
xo(a) = @. (a + (2,1,2,1,2,1,2,1)) % 3

# Movedefs
R((p,o)) = Rp(p), (Ro∘Rp)(o)
U((p,o)) = Up(p), Up(o)
x((p,o)) = xp(p), (xo∘xp)(o)
D = x ∘ x ∘ U ∘ x ∘ x
y = U ∘ D ∘ D ∘ D
F = x ∘ x ∘ x ∘ U ∘ x
B = y ∘ y ∘ F ∘ y ∘ y
z = F ∘ B ∘ B ∘ B
L = y ∘ y ∘ R ∘ y ∘ y

# Move (string) to Move (func)
function m_to_mf(m, funcs)
    func = getfield(Main,Symbol(m[1]))
    length(m)==2 ? (m[2]=='2' ? func∘func : func∘func∘func) : func
end

# String of moves to single func
str_to_func(string) = foldl(∘, map(m_to_mf, reverse(split(string))))

############

# Running
scramble = "R' U' F D2 B2 U2 L2 B2 U2 F2 R' B2 R' U' F' R D2 L D L' U2 F R' U' F"
solution = "U F2 B U B D R2 L2 F2 U' F2 L2 D F2 D' R2 D' B D L B'"

# One single function that represents this cube state! It is the solved state here.
# pos = str_to_func(scramble*" "*solution)

# Benchmarking:
@btime _pos(cube) setup=(_pos=str_to_func(scramble*" "*solution))
# 19.999 μs (1 allocation: 144 bytes)
# ((1, 2, 3, 4, 5, 6, 7, 8), (0, 0, 0, 0, 0, 0, 0, 0))
# 😲😲😲😲😲😲😲😲😲😲😲😲😲😲


text:
# binned memory pool allocator
#
# the core design is a pretty simple:
# - bin allocations into multiple pools according to their size (see `poolidx`)
# - when requested memory, check the pool for unused memory, or allocate dynamically
# - conversely, when released memory, put it in the appropriate pool for future use

## tunables

const MAX_POOL = 2^27 # 128 MiB

const USAGE_WINDOW = 5

# min and max time between successive background task iterations.
# when the pool usages don't change, scan less regularly.
#
# together with USAGE_WINDOW, this determines how long it takes for objects to get reclaimed
const MIN_DELAY = 1.0
const MAX_DELAY = 5.0

## infrastructure

Base.@kwdef struct BinnedPool <: AbstractPool
  stream_ordered::Bool

  # TODO: npools
  lock::ReentrantLock = ReentrantLock()
  active::Vector{Set{Block}} = Vector{Set{Block}}()
  cache::Vector{Vector{Block}} = Vector{Vector{Block}}()

  freed_lock::NonReentrantLock = NonReentrantLock()
  freed::Vector{Block} = Vector{Block}()
end

poolidx(n) = Base.ceil(Int, Base.log2(n))+1
poolsize(idx) = 2^(idx-1)

@assert poolsize(poolidx(MAX_POOL)) <= MAX_POOL "MAX_POOL cutoff should close a pool"

function create_pools(pool::BinnedPool, idx)
  if length(pool.cache) >= idx
    # fast-path without taking a lock
    return
  end

  @lock pool.lock begin
    while length(pool.cache) < idx
      push!(pool.active, Set{Block}())
      push!(pool.cache, Vector{Block}())
    end
  end
end

## pooling

# reclaim unused buffers
function reclaim(pool::BinnedPool, target_bytes::Int=typemax(Int))
  pool_repopulate(pool)

  @lock pool.lock begin
    freed_bytes = 0

    # process pools in reverse, to discard largest buffers first
    for pid in reverse(1:length(pool.cache))
      bytes = poolsize(pid)
      cache = pool.cache[pid]

      bufcount = length(cache)
      for i in 1:bufcount
        block = pop!(cache)

        actual_free(block; pool.stream_ordered)

        freed_bytes += bytes
        if freed_bytes >= target_bytes
          return freed_bytes
        end
      end
    end

    return freed_bytes
  end
end

# repopulate the "available" pools from the list of freed blocks
function pool_repopulate(pool::BinnedPool)
  blocks = @lock pool.freed_lock begin
    isempty(pool.freed) && return
    blocks = Set(pool.freed)
    empty!(pool.freed)
    blocks
  end

  @lock pool.lock begin
    for block in blocks
      pid = poolidx(sizeof(block))

      @inbounds active = pool.active[pid]
      @inbounds cache = pool.cache[pid]

      # mark the buffer as available
      delete!(active, block)
      push!(cache, block)
    end
  end

  return
end

function alloc(pool::BinnedPool, bytes; stream::CuStream)
  if bytes <= MAX_POOL
    pid = poolidx(bytes)
    create_pools(pool, pid)
    bytes = poolsize(pid)
  else
    pid = -1
  end

  block = nothing

  # NOTE: checking the pool is really fast, and not included in the timings
  @lock pool.lock begin
    if pid != -1 && !isempty(pool.cache[pid])
      block = pop!(pool.cache[pid])
    end
  end

  if block === nothing
    pool_repopulate(pool)

    @lock pool.lock begin
      if pid != -1 && !isempty(pool.cache[pid])
        block = pop!(pool.cache[pid])
      end
    end
  end

  if block === nothing
    block = actual_alloc(bytes; pool.stream_ordered)
  end

  if block === nothing
    GC.gc(false)
    pool_repopulate(pool)

    @lock pool.lock begin
      if pid != -1 && !isempty(pool.cache[pid])
        block = pop!(pool.cache[pid])
      end
    end
  end

  # TODO: we could return a larger allocation here, but that increases memory pressure and
  #       would require proper block splitting + compaction to be any efficient.

  if block === nothing
    reclaim(pool, bytes)
    block = actual_alloc(bytes; pool.stream_ordered)
  end

  if block === nothing
    GC.gc(true)
    pool_repopulate(pool)

    @lock pool.lock begin
      if pid != -1 && !isempty(pool.cache[pid])
        block = pop!(pool.cache[pid])
      end
    end
  end

  if block === nothing
    reclaim(pool, bytes)
    block = actual_alloc(bytes; pool.stream_ordered)
  end

  if block === nothing
    reclaim(pool, typemax(Int))
    block = actual_alloc(bytes, true; pool.stream_ordered)
  end

  if block !== nothing && pid != -1
    @lock pool.lock begin
      @inbounds active = pool.active[pid]
      @inbounds cache = pool.cache[pid]

      # mark the buffer as active
      push!(active, block)
    end
  end

  return block
end

function free(pool::BinnedPool, block; stream::CuStream)
  # was this a pooled buffer?
  bytes = sizeof(block)
  if bytes > MAX_POOL
    actual_free(block; pool.stream_ordered)
    return
  end

  # we don't do any work here to reduce pressure on the GC (spending time in finalizers)
  # and to simplify locking (preventing concurrent access during GC interventions)
  @spinlock pool.freed_lock begin
    push!(pool.freed, block)
  end
end

function cached_memory(pool::BinnedPool)
  sz = @lock pool.freed_lock mapreduce(sizeof, +, pool.freed; init=0)
  @lock pool.lock for (pid, pl) in enumerate(pool.cache)
    bytes = poolsize(pid)
    sz += bytes * length(pl)
  end
  return sz
end


text:
# TODO Add Explanation

mutable struct Online_r2 <: Scheme
    #check is last stored checkpoint
	check::Int
	#capo is the temporary fine
	capo::Int
    tail::Int
    acp::Int
    cstart::Int
    cend::Int
    numfwd::Int
    numcmd::Int
    numstore::Int
    rwcp::Int
    prevcend::Int
    firstuturned::Bool
    verbose::Int
    fstore::Union{Function,Nothing}
    frestore::Union{Function,Nothing}
    ch::Vector{Int}
    ord_ch::Vector{Int}
    num_rep::Vector{Int}
end

function Online_r2(
    #steps::Int,
    checkpoints::Int,
    fstore::Union{Function,Nothing} = nothing,
    frestore::Union{Function,Nothing} = nothing,
    anActionInstance::Union{Nothing,Action} = nothing,
    #bundle_::Union{Nothing,Int} = nothing,
    verbose::Int = 0
)
    if !isa(anActionInstance, Nothing)
        # same as default init above
        anActionInstance.actionflag = 0
        anActionInstance.iteration  = 0
        anActionInstance.cpNum      = 0
    end
    if steps<0
       error("Online_r2: negative steps")
    elseif checkpoints < 0
       error("Online_r2: negative checkpoints")
    end
    cstart = 0
    tail   = 1
    cend            = steps
    acp             = checkpoints
    numfwd          = 0
    numcmd          = 0
    numstore        = 0
    rwcp            = -1
    prevcend        = 0
    firstuturned    = false
    check = -1
    capo = -1
    ch = Vector{Int}(undef, acp)
    ord_ch = Vector{Int}(undef, acp)
    num_rep = Vector{Int}(undef, acp)
    revolve = Online_r2(check, capo, tail, acp, cstart, cend, numfwd, numcmd, numstore, rwcp, prevcend, firstuturned, verbose, fstore, frestore, ch, ord_ch, num_rep)
    return revolve
end

function next_action!(revolve::Online_r2)::Action
    # Default values for next action
    actionflag     = none
    info(" check = ", check,  " ch[check] ", ch[check], " capo ", capo)
    numcmd+=1
    cpnum = 0

    if ((check == -1) || ((ch[check] != capo) && (capo <= snaps-1)))
    #condition for takeshot for r=1
        oldcapo_o = capo
        check += 1
        ch[check] = capo
        t = 0
        if (snaps < 4)
            for i in 0:snaps-1
            #for(int i=0;i<snaps;i++)
      	        num_rep[i] = 2
            end
            incr = 2
            iter = 1
            oldind = snaps-1
        else
            iter = 1
            incr = 1
            oldind = 1
            for i in 0:snaps-1
            #for(int i=0;i<snaps;i++)
      	        num_rep[i] = 1
      	        ord_ch[i] = i
            end
            offset = snaps-1
        end
        if (capo == snaps-1)
            ind = 2
            old_f = 1
        end
        # Increase the number of takeshots and the corresponding checkpoint
        numstore+=1
        #checkpoint->number_of_writes[check]++
        return Action(store, check, -1, -1)
        #return ACTION::store
    elseif (capo < snaps-1)
    #condition for advance for r=1
        capo = oldcapo_o+1
        nunmfwd+=1
        return Action(forward, capo, oldcapo_o, cpnum)
        #return ACTION::advance
    else
    #Online_r2-Checkpointing for r=2
        if (ch[check] == capo)
        # condition for advance for r=2
            if (snaps == 1)
                capo = MAXINT-1
      		    nunmfwd+=1
                return Action(forward, capo, oldcapo_o, cpnum)
      		    #return ACTION::advance
            elseif (snaps == 2)
                capo = ch[1]+incr
      		    nunmfwd+=1
                return Action(forward, capo, oldcapo_o, cpnum)
      		    #return ACTION::advance
            elseif (snaps == 3)
                nunmfwd+=incr
      		    if (iter == 0)
      			    capo = ch[oldind]
                    for i=0:(t+1)/2
      			        #for(int i=0;i<=(t+1)/2;i++)
      				    capo += incr
      				    incr = incr + 1
      				    iter = iter + 1
                    end
      		    else
      			    capo = ch[ind]+incr
      			    incr = incr + 1
      				iter = iter + 1
                end
                actionflag = forward
                return Action(forward, capo, oldcapo_o, cpnum)
      		    #return ACTION::advance
      	    else
                if (capo == snaps-1)
      			    capo = capo+2
      		    	ind=snaps-1
      			    nunmfwd+=2
                    return Action(forward, capo, oldcapo_o, cpnum)
      			    #return ACTION::advance
                end
      		    if (output)
      			    info(" iter ", iter,  " incr ", incr, "offset", offset)
                end
      		    if (t == 0)
      			    if (iter < offset)
      				    capo = capo+1
      				    nunmfwd+=1
      			    else
      				    capo = capo+2
      				    nunmfwd+=2
                    end
      			    if (offset == 1)
      				    t+=1
                    end
                    return Action(forward, capo, oldcapo_o, cpnum)
      			    #return ACTION::advance
                end
      		    if (output)
                    info(" iter ", iter, "incr ", incr)
                end
                error(" not implemented yet")
                return Action(done, capo, oldcapo_o, cpnum)
      		    #return ACTION::error
            end
        else
            #takeshot for r=2
            if (snaps == 2)
                ch[1] = capo
      		    incr+=1
      		    #Increase the number of takeshots and the corresponding checkpoint
      		    numstore+=1
      		    #checkpoint->number_of_writes[1]++
                return Action(store, check, -1, -1)
      		    #return ACTION::takeshot
            elseif (snaps == 3)
                ch[ind] = capo
      		    check = ind
      		    info(" iter ", iter, " num_rep[1] ", num_rep[1])
      		    if (iter == num_rep[1])
      			    iter = 0
      			    t+=1
      		    	oldind = ind
      			    num_rep[1]+=1
      			    ind = 2 - num_rep[1]%2
      			    incr=1
                  end
      		    #Increase the number of takeshots and the corresponding checkpoint
      		    numstore+=1
      		    #checkpoint->number_of_writes[check]++
                return Action(store, check, -1, -1)
      		    #return ACTION::takeshot
            else
                if (capo < snaps+2)
      			    ch[ind] = capo
      			    check = ind
      			    if (capo == snaps+1)
      				    oldind = ord_ch[snaps-1]
      				    ind = ch[ord_ch[snaps-1]]
      				    if (output)
                            info(" oldind ", oldind, " ind ", ind)
                        end
                        for k=snaps-1:-1:2
      				    #for(int k=snaps-1;k>1;k--)
      					    ord_ch[k]=ord_ch[k-1]
      					    ch[ord_ch[k]] = ch[ord_ch[k-1]]
                        end
      				    ord_ch[1] = oldind
      				    ch[ord_ch[1]] = ind
      				    incr=2
      				    ind = 2
      				    if (output)
      				    	info(" ind ", ind, " incr ", incr, " iter ", iter)
                            for j=0:snaps-1
      				    	#for(int j=0;j<snaps;j++)
      				    	    info(" j ", j, " ord_ch ", ord_ch[j], " ch ", ch[ord_ch[j]], " rep ", num_rep[ord_ch[j]])
                            end
                        end
                    end
      			    #Increase the number of takeshots and the corresponding checkpoint
      			    numstore+=1
      			    #checkpoint->number_of_writes[check]++
                    return Action(store, check, -1, -1)
      			    #return ACTION::takeshot
                end
            end
      		if (t == 0)
      			if (output)
                    info(" ind ", ind, " incr ",  incr, " iter ", iter, " offset ", offset)
                end
                if (iter == offset)
      				offset=offset-1
      				iter = 1
      				check = ord_ch[snaps-1]
      				ch[ord_ch[snaps-1]] = capo
      				oldind = ord_ch[snaps-1]
      				ind = ch[ord_ch[snaps-1]]
      				if (output)
                        info(" oldind " , oldind , " ind " , ind)
                    end
                    for k=snaps-1:-1:incr+1
      				#for(int k=snaps-1;k>incr;k--)
      					ord_ch[k]=ord_ch[k-1]
      					ch[ord_ch[k]] = ch[ord_ch[k-1]]
                    end
      				ord_ch[incr] = oldind
      				ch[ord_ch[incr]] = ind
      				incr+=1
      				ind=incr
      				if (output)
                        info(" ind ", ind, " incr ", incr, " iter ", iter)
                        for j=0:snaps-1
      					#for(int j=0;j<snaps;j++)
                            info(" j ", j << " ord_ch ", ord_ch[j], " ch ", ch[ord_ch[j]], " rep ", num_rep[ord_ch[j]])
                        end
                    end
      			else
      				ch[ord_ch[ind]] = capo
      				check = ord_ch[ind]
      				iter+=1
      				ind+=1
      				if (output)
                        info(" xx ind ", ind, " incr ", incr, " iter ", iter)
                    end
                end
      			#Increase the number of takeshots and the corresponding checkpoint
      			numstore=numstore+1
      			#checkpoint->number_of_writes[check]++
                return Action(store, check, -1, -1)
      			#return ACTION::takeshot
      		end
        end
    end
    actionflag = done
    return Action(done, capo, oldcapo_o, cpnum)
    #return ACTION::terminate # This means that the end of Online_r2 Checkpointing for r=2 is reached and
    #  another Online_r2 Checkpointing class must be started
end

function Revolve(
    online::Online_r2,
    checkpoints::Int,
    fstore::Union{Function,Nothing} = nothing,
    frestore::Union{Function,Nothing} = nothing;
)
    if !isa(anActionInstance, Nothing)
        # same as default init above
        anActionInstance.actionflag = 0
        anActionInstance.iteration  = 0
        anActionInstance.cpNum      = 0
    end
    cstart = online_r2.cstart
    tail   = 1
    cend            = online_r2.steps
    acp             = online_r2.checkpoints
    numfwd          = online_r2.numfwd
    numinv          = online_r2.numinv
    numstore        = online_r2.numstore
    rwcp            = online_r2.rwcp
    prevcend        = online_r2.prevcend
    firstuturned    = online_r2.firstuturned #true
    stepof = Vector{Int}(undef, acp+1)

    revolve = Revolve(steps, bundle, tail, acp, cstart, cend, numfwd, numinv, numstore, rwcp, prevcend, firstuturned, stepof, verbose, fstore, frestore)

    if verbose > 0
        predfwdcnt = forwardcount(revolve)
        if predfwdcnt == -1
            error("Revolve: error returned by  revolve::forwardcount")
        else
            @info "prediction:"
            @info " overhead forward steps : $predfwdcnt"
            @info " overhead factor        : $(predfwdcnt/steps)"
        end
    end
    return revolve
end

function checkpoint_struct(body::Function, alg::Online_r2, model_input::MT, shadowmodel::MT) where {MT}
    model = deepcopy(model_input)
    storemap = Dict{Int32,Int32}()
    check = 0
    model_check = Array{MT}(undef, alg.acp)
    copyto!(model_input, deepcopy(model))
    while true
        next_action = next_action!(alg)
        if (next_action.actionflag == Checkpointing.store)
            check = check+1
            storemap[next_action.iteration-1]=check
            model_check[check] = deepcopy(model)
        elseif (next_action.actionflag == Checkpointing.forward)
            for j= next_action.startiteration:(next_action.iteration - 1)
                body(model)
            end
        elseif (next_action.actionflag == Checkpointing.firstuturn)
            error("Unexpected firstuturn")
        elseif (next_action.actionflag == Checkpointing.uturn)
            error("Unexpected uturn")
        elseif (next_action.actionflag == Checkpointing.restore)
            error("Unexpected restore")
        elseif next_action.actionflag == Checkpointing.done
            info("Done with online phase")
            break
        end
    end
end


text:
# Dependencies
using Distributed
@everywhere push!(LOAD_PATH, "./")
using CpelNano

# IO
dataDir = dirname(dirname(pathof(CpelNano)))
nano = "$(dataDir)/examples/targeted_example/nanopolish/targeted_example_noise2.0_sorted_methylation_calls.sorted.tsv"
fasta = "$(dataDir)/examples/targeted_example/reference/hg38_chr17_42879379_43290399.fa"
bed = "$(dataDir)/examples/targeted_example/regions_of_interest/brca1.bed"
outdir = "$(dataDir)/examples/targeted_example/cpelnano/"

# Output files
mml_path = "$(outdir)/targeted_example_mml.txt"
nme_path = "$(outdir)/targeted_example_nme.txt"
theta_path = "$(outdir)/targeted_example_theta.txt"

# Remove output in outdir
isfile(mml_path) && rm(mml_path)
isfile(nme_path) && rm(nme_path)
isfile(theta_path) && rm(theta_path)

# Show region of interest
CpelNano.print_log("Regions of interest")
print(String(read(bed)))

# Configuration
min_cov = 10.0; max_size_subreg = 350; size_est_reg = 3000; max_em_init = 5; max_em_iters = 50;
config = CpelNano.CpelNanoConfig(min_cov, max_size_subreg, size_est_reg, max_em_init, max_em_iters);
config.out_dir = outdir; config.out_prefix = "targeted_example"; config.bed_reg = bed; config.verbose = false;

# Analyze each region in chr
CpelNano.analyze_nano(nano,fasta,config)

# Print files
println("*******************************MML*******************************")
println(String(read(mml_path)))
println("*****************************************************************")
readline()
println("*******************************NME*******************************")
println(String(read(nme_path)))
println("*****************************************************************")
readline()


text:
"""
    opnorm(A::IntervalMatrix, p::Real=Inf)

The matrix norm of an interval matrix.

### Input

- `A` -- interval matrix
- `p` -- (optional, default: `Inf`) the class of `p`-norm

### Notes

The matrix ``p``-norm of an interval matrix ``A`` is defined as

```math
    ‖A‖_p := ‖\\max(|\\text{inf}(A)|, |\\text{sup}(A)|)‖_p
```

where ``\\max`` and ``|·|`` are taken elementwise.
"""
function LinearAlgebra.opnorm(A::IntervalMatrix, p::Real=Inf)
    if p == Inf
        return _opnorm_inf(A)
    elseif p == 1
        return _opnorm_1(A)
    else
        error("the interval matrix norm for this value of p=$p is not implemented")
    end
end

# The operator norm in the infinity norm corresponds to the
# maximum absolute value row-sum.
function _opnorm_inf(A::IntervalMatrix{T}) where {T}
    m, n = size(A)
    res = zero(T)
    @inbounds @simd for i in 1:m
        acc = zero(T)
        for j in 1:n
            x = A[i, j]
            acc += max(abs(inf(x)), abs(sup(x)))
        end
        if acc > res
            res = acc
        end
    end
    return res
end

# The operator norm in the 1-norm corresponds to the
# maximum absolute value column-sum.
function _opnorm_1(A::IntervalMatrix{T}) where {T}
    m, n = size(A)
    res = zero(T)
    @inbounds @simd for j in 1:n
        acc = zero(T)
        for i in 1:m
            x = A[i, j]
            acc += max(abs(inf(x)), abs(sup(x)))
        end
        if acc > res
            res = acc
        end
    end
    return res
end

"""
    diam_norm(A::IntervalMatrix, p=Inf)

Return the diameter norm of the interval matrix.

### Input

- `A` -- interval matrix
- `p` -- (optional, default: `Inf`) the `p`-norm used; valid options are:
         `1`, `2`, `Inf`

### Output

The operator norm, in the `p`-norm, of the scalar matrix obtained
by taking the element-wise `diam` function, where `diam(x) := sup(x) - inf(x)`
for an interval `x`.

### Notes

This function gives a measure of the *width* of the interval matrix.
"""
function diam_norm(A::IntervalMatrix, p=Inf)
    return opnorm(diam.(A), p)
end