text:
## Models

## A model implements the observable interface
## interface

## methods

##
## signal/slot style
## an observable has events which emit signals
## an observer can connect a slot (a method of the observer)
## connect(observed, signal, observer, slot)

function connect(o::Observable, signal::String, observer, slot)
    if !haskey(o.observers, signal)
        o.observers[signal] = Dict()
    end
    id = randstring(10)
    o.observers[signal][id] = (observer, slot)
    id
end
## A slot can be a simple function too.
function connect(o::Observable, signal::String, slot::Function)
    connect(o, signal, nothing, slot)
end

## for do ntation
connect(slot::Function, o::Observable, signal::String) = connect(o, signal, slot)

## we disconnect by id.
## XXX should include means to block
function disconnect(o::Observable, id::String)
    for (k,v) in o.observers
        if haskey(v, id)
            delete!(v, id)
            break
        end
    end
end

## when calling notify the values must be passed in
##
## XXX put in code to catch if observer is still extant? Remove if not...
function notify(o::Observable, signal, values...)
    if !haskey(o.observers, signal) return end
    obs = o.observers[signal]
    for (k, v) in obs
        observer, slot = v
        if isa(observer, Nothing)
            slot(values...)
#            try slot(value...) catch e nothing end
        else
            slot(observer, values...)
#            try slot(observer, values...) catch e nothing end
        end
    end
end

## methods
## This allows o[:value] to work.
getValue(model::Observable) = model.value
setValue(model::Observable, value::Nothing; signal::Bool=true) = nothing
function setValue(model::Observable, value; signal::Bool=true)
     if value != model.value
         model.value = value
         signal && notify(model, "valueChanged", getValue(model))
     end
     value
 end

type EventModel <: Observable
    observers::Dict
    EventModel() = new(Dict())
end

## Model has value
abstract Model <: Observable

type ItemModel <: Model
    observers::Dict
    ##
    value
    ItemModel() = new(Dict(), nothing)
    ItemModel(value) = new(Dict(), value)
end

## signals
## valueChanged, (value)

abstract AbstractArrayModel <: Model
type VectorModel <: AbstractArrayModel
    observers::Dict
    ##
    value
    items::Vector
    VectorModel() = new(Dict(), nothing, {})
    VectorModel(items) = new(Dict(), items[1], items)
    VectorModel(items, value) = new(Dict(), value, items)
end

## store index in 1:100, 1:100, getValue of widget converts
type TwoDSliderModel <: AbstractArrayModel
    observers::Dict
    ##
    value
    items1::Vector
    items2::Vector
    TwoDSliderModel() = new(Dict(), nothing, {}, {})
    TwoDSliderModel(items1, items2) = new(Dict(), [items1[1], items2[1]], items1, items2)
    TwoDSliderModel(items1, items2, value) = new(Dict(), value, items1, items2)
end

function getValue(o::TwoDSliderModel)
    i, j = o.value
    x = iround(1 + (length(o.items1) - 1)/(100 -1)*(i - 1))
    y = iround(1 + (length(o.items2) - 1)/(100 -1)*(j - 1))
    [x,y]
end

## value in items coordinates, convert to index
function setValue{T <: Number}(o::TwoDSliderModel, value::Vector{T}; signal::Bool=true)
    x, y = value
    function item2ind(item, items)
        ix = indmin(abs(items - item))
        iround( 1 + (100 - 1)/(length(items) - 1)*(ix -1))
    end
    i = item2ind(x, o.items1)
    j = item2ind(y, o.items2)
    if o.value != [i,j]
        o.value = [i,j]
        signal && notify(o, "valueChanged", getValue(o))
    end
end

type ArrayModel <: AbstractArrayModel
    observers::Dict
    ##
    value
    items::Array
    ArrayModel() = new(Dict(), nothing, {})
    ArrayModel(items) = new(Dict(), items[1], items)
    ArrayModel(value, items) = new(Dict(), value, items)
end

## o[:items] ...
getItems(model::AbstractArrayModel) = model.items
function setItems(model::AbstractArrayModel, items)
    old_items = model.items
    model.items = items
    notify(model, "itemsChanged", items, old_items)
end

### Stores
## A store is comprised of a vector of items wi
## A model item is nothing more than a composite type, but we offer
## this type in case one wants to subtype
abstract AbstractStoreItem

## May be needed by some (e.g., Tk)
to_string(m::Any, x) = string(x)
from_string(m::Any, field::Symbol, x) = convert(eltype(m.(field)), x)

abstract DataStore <: Object

type Store <: DataStore
    model::Observable
    types::Vector{DataType}
    items::Vector
    Store(types...) = new(ItemModel(), [types...],  Vector[])
end

length(s::DataStore) = length(s.items)
size(s::DataStore) = [length(s), length(s.items[1])]

getindex(s::DataStore, i::Int) = s.items[i]
getindex(s::DataStore, i::Int, j::Int) = s.items[i][j]
setindex!(s::DataStore, val, i::Int) = replace!(s, i, val)
function setindex!(s::DataStore, val, i::Int, j::Int)
    cur = s[i]
    cur[j] = val
    replace!(s, i, cur)
end

function insert!(s::DataStore, i::Int, vals)
    ## verify types
    vals = [convert(T,v) for (T, v) in zip(s.types, vals)]
    insert!(s.items, i, vals)
    notify(s.model, "rowInserted", i)
end

function splice!(s::DataStore, i::Int)
    items = splice!(s.items, i)
    notify(s.model, "rowRemoved", i)
    items
end

function replace!(s::DataStore, i::Int, item)
    if i < 1 || i > length(s.items)
        error("Index is out of bounds for replacing")
    end
    s.items[i] = item
    notify(s.model, "rowUpdated", i)
end

push!(s::DataStore, val) = insert!(s, length(s)+1, val)
prepend!(s::DataStore, val) = insert!(s, 1, val)
append!(s::DataStore, vals::Vector) = map(val -> push!(s, val), vals)

pop!(s::DataStore) = splice!(s::DataStore, length(s))

## Trees
type TreeNode <: Object
    text::String
    icon
    data                        # DataType
    index                       # from treeview (couples to tree view, don't like...)
    parent                      # nothing,
    children::Vector            # TreeNodes
end

function node_to_values(node::TreeNode)
    ds = node.data
    String[string(d) for d in ds]
end

abstract AbstractTreeStore <: DataStore

## Would like to put in data type here...
type TreeStore <: DataStore
    model::Observable
    types::Vector{DataType}
    children::Vector{TreeNode}
    attrs::Dict
    TreeStore(types...) = new(ItemModel(), [types...], TreeNode[], Dict())
end

## A store for display through `treeview`
##
## a treestore holds tree nodes. Each node holds a label, optional data [and an icon]. The
## nodes are most easily created from a label and data from insert!
##
## Nodes should have data all from the same CompositeType (as with `storeview`). It is not
## necessary, as often nodes that expand are not of the same type (directories are different from files, say).
##
## A node in a tree has a corresponding path -- a vector of
## indices. The functions `node_to_path` and `path_to_node`
## translate.
##
## Nodes can be expanded or collapsed through `expand_node` and `collapse_node`
##
##
treestore(types...) = TreeStore(types...)

## some store methods
length(store::TreeStore)  = length(store.types)

function index_of(node::TreeNode)
    parentnode = node.parent
    findfirst(parentnode.children, node)
end

## A path is 1-based index through children
## this finds path from a node
function node_to_path(node::TreeNode)
    path = Int[]
    path = [index_of(node)]
    parent = node.parent
    while !isa(parent, Union(TreeStore, Nothing))
        unshift!(path, index_of(parent))
        parent = parent.parent
    end
    path
end

## this finds node from a path
function path_to_node(store::TreeStore, path::Vector{Int})
    i = shift!(path)
    node = store.children[i]
    while length(path) > 0
        i = shift!(path)
        node = node.children[i]
    end
    node
end

function insert!(store::TreeStore, parentnode::Union(Nothing, TreeNode), i::Int, childnode::TreeNode)
    if parentnode == nothing
        parentnode = store
    end

    childnode.parent = parentnode
    insert!(parentnode.children, i, childnode)
    notify(store.model, "insertNode", parentnode, i, childnode)
end
## data an instance of a type of nothing
function insert!(store::TreeStore, parentnode::Union(Nothing, TreeNode), i::Int, text::String, data)
    childnode = TreeNode(text, nothing, data, nothing, parentnode, {})
    insert!(store, parentnode, i, childnode)
    childnode
end
function push!(tr::TreeStore, parentnode::Union(Nothing, TreeNode), node::TreeNode)
    n= parentnode == nothing ? length(tr.children) : length(parentnode.children)
    insert!(tr, parentnode, n+1, node)
end
function push!(tr::TreeStore, parentnode::Union(Nothing, TreeNode), text::String, data)
    childnode = TreeNode(text, nothing, data, nothing, parentnode, {})
    push!(tr, parentnode, childnode)
end
unshift!(tr::TreeStore, parentnode::Union(Nothing, TreeNode), node::TreeNode) = insert!(tr, parentnode, 1, node)
unshift!(tr::TreeStore, parentnode::Union(Nothing, TreeNode), text::String, data) = insert!(tr, parentnode, 1, node)

function splice!(store::TreeStore, parentnode::Union(Nothing, TreeNode), i::Int)
    if isa(parentnode, Nothing)
        parentnode = store
    end
    splice!(parentnode.children, i)
    notify(store.model, "removeNode", parentnode, i)
end

## remove node from store
function delete!(store::TreeStore, node::TreeNode)
    parentnode = node.parent
    if isa(parentnode, TreeStore)
        parentnode = nothing
    end
    ### how to get indexof node
    i = index_of(node)

    splice!(store, parentnode, i)
end

## update node properties text and/or data (XXX Icons!)
function update_node(store::TreeStore, node::TreeNode; text::Union(Nothing, String)=nothing, data::Union(Nothing, Any)=nothing)
    if !isa(text, Nothing)
        node.text = text
    end
    if !isa(data, Nothing)
        node.data = data
    end
    notify(store.model, "updatedNode", node)
end


text:
"""

`SPGBoxResult`: Data structure that contains the results of the `spgbox`
optimization method.

  `x`: Vector containing final point.

  `f`: Final objective function value.

  `gnorm`: Norm of the projected gradient on the bounds.

  `nit`: Number of iterations performed.

  `nfeval`: Number of function evaluations.

  `ierr`: Status of execution:
         0. Converged successfully to solution.
         1. Maximum number of iterations achieved.
         2. Maximum number of function evaluations achieved.

"""
struct SPGBoxResult
  x::Vector{Float64}
  f::Float64
  gnorm::Float64
  nit::Int64
  nfeval::Int64
  ierr::Int64
end

function Base.show(io::IO, R::SPGBoxResult )
  println("")
  println(" SPGBOX RESULT: ")
  println("")
  R.ierr == 0 && println(" Convergence achieved. ")
  R.ierr == 1 && println(" Maximum number of iterations (nitmax) reached.")
  R.ierr == 2 && println(" Maximum number of function evaluations (nfevalmax) reached.")
  println("")
  println(" Final objective function value = ", R.f)
  print(" Best solution found = [ ")
  for i in 1:min(length(R.x)-1,3)
    print("$(R.x[i]), ")
  end
  if length(R.x) > 4
    print("..., ")
  end
  println(R.x[length(R.x)],"]")
  println(" Projected gradient norm = ", R.gnorm)
  println("")
  println(" Number of iterations = ", R.nit)
  println(" Number of function evaluations = ", R.nfeval)
end

ext: jl


text:
function solve()
    a, b = [parse(Int, x) for x in split(readline())]
    c, d = [parse(Int, x) for x in split(readline())]
    b - c
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
module FEBasesTests

using Test
using Gridap

model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=(4,4))
order = 1
diritag = "boundary"
fespace = H1ConformingFESpace(Float64,model,order,diritag)

bh = FEBasis(fespace)

@test isa(Triangulation(bh),Triangulation)

ufun(x) = x[1]

uh = interpolate(fespace,ufun)

@test isa(bh,FEBasis)
@test isa(+bh,FEBasis)
@test isa(-bh,FEBasis)
@test isa(∇(bh),FEBasis)
@test isa(bh+uh,FEBasis)
@test isa(uh-bh,FEBasis)
@test isa(bh-1.0,FEBasis)
@test isa(1.0+bh,FEBasis)
@test isa(inner(bh,bh),CellMap{Point{2},1,Float64,3})
@test isa(inner(bh,uh),CellMap{Point{2},1,Float64,2})

trian = Triangulation(model)
quad = CellQuadrature(trian,degree=2)

σfun(x,u,i) = x[1] + u
σ(u,i) = CellBasis(trian,σfun,u,i)
ids = ones(Int,ncells(trian))
cm = inner(bh,σ(bh,ids))
m = integrate(cm,trian,quad)
_ = collect(m)

tags = [7,6]
btrian = BoundaryTriangulation(model,tags)
bquad = CellQuadrature(btrian,degree=2)

bbh = restrict(bh,btrian)
cm = inner(bbh,bbh)

ca = integrate(cm,btrian,bquad)
@test isa(ca,CellArray{Float64,2})
_ = collect(ca)

buh = restrict(uh,btrian)
cm = inner(bbh,buh)
ca = integrate(cm,btrian,bquad)
@test isa(ca,CellArray{Float64,1})
_ = collect(ca)

strian = SkeletonTriangulation(model,"interior")
squad = CellQuadrature(strian,degree=2)

sbh = restrict(bh,strian)
cm = inner(jump(sbh),mean(sbh))
ca = integrate(cm,strian,squad)

suh = restrict(uh,strian)
cm = inner(jump(sbh*mean(suh)),mean(suh))
ca = integrate(cm,strian,squad)

model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=(4,4))
order = 2
diritag = "boundary"
T = VectorValue{2,Float64}
fespace = H1ConformingFESpace(T,model,order,diritag)
ufun(x) = VectorValue(x[2],x[1])
uh = interpolate(fespace,ufun)
bh = FEBasis(fespace)
@test isa(ε(bh),FEBasis)
@test isa(div(bh),FEBasis)
@test isa(curl(bh),FEBasis)
@test isa(inner(ε(bh),ε(bh)),CellMap{Point{2},1,Float64,3})
@test isa(inner(ε(bh),ε(uh)),CellMap{Point{2},1,Float64,2})

σfun(x,ε) = 4*tr(ε)*one(ε) + 2*ε
σ(u) = CellBasis(trian,σfun,u)
cm = inner(ε(bh),σ(ε(bh)))
m = integrate(cm,trian,quad)
_ = collect(m)

end # module


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

using Test
using Dates

using TOML: TOML, parse, tryparse, ParserError, Internals, print

function roundtrip(data)
    mktemp() do file, io
        data_parsed = TOML.parse(data)
        TOML.print(io, data_parsed)
        close(io)
        data_roundtrip = TOML.parsefile(file)
        return isequal(data_parsed, data_roundtrip)
    end
end

include("readme.jl")
include("toml_test.jl")
include("values.jl")
include("invalids.jl")
include("error_printing.jl")
include("print.jl")
include("parse.jl")

@inferred TOML.parse("foo = 3")

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
push!(_adbackends, :forwarddiff)

_gradient(f, x::Number, ::Val{:forwarddiff}) = ForwardDiff.derivative(f, x)
_gradient(f, x, ::Val{:forwarddiff}) = ForwardDiff.gradient(f, x)

_jacobian(f, x, ::Val{:forwarddiff}) = ForwardDiff.jacobian(f, x)


text:
"""
Ase fileformat
ASCII Segment File
Described as route network. Each line represents a route segment

See EUROCONTROL NEST Manual for Ase fileformat description

Ase DataFrame Column Names
#   Field           Type        Size    Comment
1   :FLIGHTCOUNT    Float64             normally number of flights using this route segment, could be a load
2   :SEGMENTPARITY  Int64               0=NO, 1=ODD, 2=EVEN, 3=ODD_LOW, 4=EVEN_LOW, 5=ODD_HIGH, 6=EVEN_HIGH
3   :SEGMENTTYPE    Int64               0=NO, 1=NORMAL, 2=ARRIVAL, 3=DEPARTURE (permanent rte segment)
                                        20=NO, 21=NORMAL, 22=ARRIVAL, 23=DEPARTURE (CDR Generic)
                                        40=NO, 41=NORMAL, 42=ARRIVAL, 43=DEPARTURE (CDR 1 )
                                        60=NO, 61=NORMAL, 62=ARRIVAL, 63=DEPARTURE (CDR 2)
                                        80=NO, 81=NORMAL, 82=ARRIVAL, 83=DEPARTURE (CDR 3)
                                        100=NO, 101=NORMAL, 102=ARRIVAL, 103=DEPARTURE (CDR 1+2)
                                        120=NO, 121=NORMAL, 122=ARRIVAL, 123=DEPARTURE (CDR 1+3)
4   :LATBEGINSEGMENT_DEG Float64
5   :LONBEGINSEGMENT_DEG Float64
6   :LATENDSEGMENT_DEG   Float64
7   :LONENDSEGMENT_DEG   Float64
8   :SEGMENTNAME    String      11      routePointNameBegin_routePointNameEnd (separator '_')

%% = TBD
%  = TBD
*  = TBD
"""
module Ase

export read

using Format
using CSV
using Dates

const fileformat = Dict(1=>Float64, 2=>Int64, 3=>Int64, 4=>Float64,
5=>Float64, 6=>Float64, 7=>Float64, 8=>String)

const header_format = ["FLIGHTCOUNT", "SEGMENTPARITY", "SEGMENTTYPE",
"LATBEGINSEGMENT_DEG", "LONBEGINSEGMENT_DEG", "LATENDSEGMENT_DEG", "LONENDSEGMENT_DEG",
"SEGMENTNAME"]

function read(file)
    df = CSV.read(file, delim=" ", types=fileformat, header=header_format,
    copycols=true)
    reformat!(df)
    return df
end

function reformat!(df)
    df[:,:LATBEGINSEGMENT_DEG] = df[:,:LATBEGINSEGMENT_DEG] / 60.0
    df[:,:LONBEGINSEGMENT_DEG] = df[:,:LONBEGINSEGMENT_DEG] / 60.0
    df[:,:LATENDSEGMENT_DEG] = df[:,:LATENDSEGMENT_DEG] / 60.0
    df[:,:LONENDSEGMENT_DEG] = df[:,:LONENDSEGMENT_DEG] / 60.0
end

end # module


text:
LCR6xxx_Measurements_Dict = Dict()
# Set Measurement
LCR6xxx_Measurements_Dict["max"] = 0
LCR6xxx_Measurements_Dict["min"] = 1
LCR6xxx_Measurements_Dict["amplitude"] = 2

# Set Measurement
function set_meas(obj::T;ch::String="1",fct::String="max",par::Int=1,vw::Int=1) where {T<:LCR6xxx}
	fct = LCR6xxx_Measurements_Dict[fct]
	write(obj,"VBS app.Measure.P$par.ParamEngine=$fct")
	write(obj,"VBS app.Measure.P$par.Source1=$ch")
	write(obj,"VBS app.Measure.P$par.View=$vw")
end
# Get Measurement
function get_meas(obj::T;par::Int=1) where {T<:LCR6xxx}
	query(obj,"VBS? Return=app.Measure.P$par.Out.Result.Value")
end


text:
# Rohan McLure, Australian National University (2019)
# Distributed implementation of the (Star) Stencil Parallel Research Kernel
using Distributed

function main()
    if length(ARGS) != 3 || nprocs() == 1
        println("Usage: julia -p <# workers> stencil.jl <# iterations> <# order> <# radius>")
        exit(1)
    end

    argv = map(x -> parse(Int, x), ARGS)
    iterations, order, r = argv
    ranks = nworkers()
    procs_x, procs_y = factor(ranks)
    out_channels = Matrix{Vector{Union{RemoteChannel,Nothing}}}(undef, procs_y, procs_x)
    dims = Matrix{Tuple{Int,Int}}(undef, procs_y, procs_x)
    std_width = Int(ceil(order / procs_x))
    std_height = Int(ceil(order / procs_y))
    last_width = order % std_width == 0 ? std_width : order % std_width
    last_height = order % std_height == 0 ? std_height : order % std_height
    for i in 1:procs_y
        for j in 1:procs_x
            v = Vector{Union{RemoteChannel, Nothing}}(undef, 4)
            v[1] = i != 1 ? RemoteChannel(() -> Channel{Matrix{Float64}}(0), (i-2)*procs_x+j) : nothing
            v[2] = j != procs_x ? RemoteChannel(() -> Channel{Matrix{Float64}}(0), (i-1)*procs_x+j+1) : nothing
            v[3] = i != procs_y ? RemoteChannel(() -> Channel{Matrix{Float64}}(0), i*procs_x+j) : nothing
            v[4] = j != 1 ? RemoteChannel(() -> Channel{Matrix{Float64}}(0), (i-1)*procs_x+j-1) : nothing
            out_channels[i,j] = v
            # Assigning coordinates
            height = i == procs_y ? last_height : std_height
            width = j == procs_x ? last_width : std_width
            dims[i,j] = (height, width)
        end
    end
    in_channels = Matrix{Vector{Union{RemoteChannel,Nothing}}}(undef, procs_y, procs_x)
    for i in 1:procs_y
        for j in 1:procs_x
            v = Vector{Union{RemoteChannel, Nothing}}(undef, 4)
            v[1] = i != 1 ? out_channels[i-1,j][3] : nothing
            v[2] = j != procs_x ? out_channels[i,j+1][4] : nothing
            v[3] = i != procs_y ? out_channels[i+1,j][1] : nothing
            v[4] = j != 1 ? out_channels[i,j-1][2] : nothing
            in_channels[i,j] = v
        end
    end

    futures = Matrix{Future}(undef, procs_y, procs_x)
    @sync for i in 1:procs_y
        for j in 1:procs_x
            wid = (i-1)*procs_x+j
            id = workers()[wid]
            futures[i,j] = @spawnat id do_stencil(iterations, order, dims[i,j], (i,j), (procs_y, procs_x), r, out_channels[i,j], in_channels[i,j])
        end
    end

    results = map(fetch, futures)

    active_points = (order - 2*r)^2
    norm = sum(map(x->x[1], results))
    norm /= active_points
    time = maximum(map(x->x[2], results))
    reference_norm = 2 * (iterations + 1)
    if abs(norm - reference_norm) > 1.e-8
        println("ERROR: L1 norm = $norm Reference L1 norm = $reference_norm")
    else
        println("Solution validates")
    end
    stencil_size = 4*r+1
    flops = (2*stencil_size+1) * active_points
    avgtime = time/iterations
    println("Rate (MFlops/s): ",1.e-6*flops/avgtime, " Avg time (s): ",avgtime)
end

@everywhere function do_stencil(
    iterations::Int,
    order::Int,
    local_dims::Tuple{Int,Int},
    proc_dims::Tuple{Int,Int},
    proc_setup::Tuple{Int,Int},
    r::Int,
    out_channels::Vector{Union{RemoteChannel,Nothing}},
    in_channels::Vector{Union{RemoteChannel,Nothing}}
)
    ub, rb, db, lb = out_channels .!== nothing
    h, w = local_dims
    id_y, id_x = proc_dims
    procs_y, procs_x = proc_setup
    std_height, std_width = Int(ceil(order/procs_y)), Int(ceil(order/procs_x))
    i_start, j_start = (id_y-1)*std_height, (id_x-1)*std_width
    i_end, j_end = i_start + h - 1, j_start + w - 1
    i_o_start, j_o_start = (ub ? 1 : r+1), (lb ? 1 : r+1)
    i_o_end, j_o_end = (db ? h : h-r), (rb ? w : w-r)
    A = zeros(Float64, (local_dims .+ 2*r)...)
    for i in 1:h
        for j in 1:w
            A[i+r,j+r] = (i_start + i - 1) + (j_start + j - 1)
        end
    end
    B = zeros(Float64, local_dims...)
    W = zeros(Float64, 2*r+1, 2*r+1)
    for i in 1:r
        W[r+1, r+i+1] = 1.0/(2*i*r)
        W[r+i+1, r+1] = 1.0/(2*i*r)
        W[r+1, r-i+1] = -1.0/(2*i*r)
        W[r-i+1, r+1] = -1.0/(2*i*r)
    end

    u_out, r_out, d_out, l_out = [Matrix{Float64}(undef, (i % 2 == 0 ? (h,r) : (r,w))...) for i in 1:4]
    u_out_v, r_out_v, d_out_v, l_out_v = [
        view(A, r+1:2*r, r+1:r+w),
        view(A, r+1:r+h, w+1:r+w),
        view(A, h+1:r+h, r+1:r+w),
        view(A, r+1:r+h, r+1:2*r)
    ]
    u_in_v, r_in_v, d_in_v, l_in_v = [
        view(A, 1:r, r+1:r+w),
        view(A, r+1:r+h, r+w+1:2*r+w),
        view(A, r+h+1:2*r+h, r+1:r+w),
        view(A, r+1:r+h, 1:r)
    ]

    # Channels
    u_out_c, r_out_c, d_out_c, l_out_c = out_channels
    u_in_c, r_in_c, d_in_c, l_in_c = in_channels

    local t0, t1

    for k in 0 : iterations
        if k == 1
            t0 = time_ns()
        end
        # Send
        @sync begin
            if ub
                @async begin
                    copy!(u_out, u_out_v)
                    put!(u_out_c, u_out)
                end
                @async begin
                    u_in = take!(u_in_c)
                    copy!(u_in_v, u_in)
                end
            end
            if rb
                @async begin
                    copy!(r_out, r_out_v)
                    put!(r_out_c, r_out)
                end
                @async begin
                    r_in = take!(r_in_c)
                    copy!(r_in_v, r_in)
                end
            end
            if db
                @async begin
                    copy!(d_out, d_out_v)
                    put!(d_out_c, d_out)
                end
                @async begin
                    d_in = take!(d_in_c)
                    copy!(d_in_v, d_in)
                end
            end
            if lb
                @async begin
                    copy!(l_out, l_out_v)
                    put!(l_out_c, l_out)
                end
                @async begin
                    l_in = take!(l_in_c)
                    copy!(l_in_v, l_in)
                end
            end
        end

        for j in j_o_start : j_o_end
            for i in i_o_start : i_o_end
                tmp = 0.0
                for ii in -r:r
                    @inbounds tmp += W[r+ii+1, r+1] * A[i+ii+r, j+r]
                end
                for jj in -r:-1
                    @inbounds tmp += W[r+1, r+jj+1] * A[i+r, j+jj+r]
                end
                for jj in 1:r
                    @inbounds tmp += W[r+1, r+jj+1] * A[i+r, j+jj+r]
                end
                @inbounds B[i,j] += tmp
            end
        end

        for j in 1:w
            for i in 1:h
                @inbounds A[i+r, j+r] += 1.0
            end
        end
    end

    t1 = time_ns()

    local_norm = 0.0
    for i in i_o_start : i_o_end
        for j in j_o_start : j_o_end
            local_norm += abs(B[i,j])
        end
    end

    return (local_norm, (t1-t0)*1.e-9)
end

function factor(x::Int)
    s = Int(floor(sqrt(x)))
    for p in s:-1:1
        if x % p == 0
            return (p, x ÷ p)
        end
    end
end

main()


text:
using Arrow, DataFrameMacros, DataFrames, MixedModels
    # read in the data
axb_data = DataFrame(Arrow.Table("./instructors/Patricia_Amillos/data/axb_data.arrow"))

    # for a binary response it can help to summarize the number of true/false responses
acctbl = combine(groupby(axb_data, [:ID, :accuracy]), nrow)
unstack(acctbl, :ID, :accuracy, :nrow)

    # similarly for item
unstack(combine(groupby(axb_data, [:item, :accuracy]), nrow), :item, :accuracy, :nrow)

    # set the contrasts
contr = Dict(
    :ID => Grouping(),
    :item => Grouping(),
    :Task => EffectsCoding(),
    :training => EffectsCoding(),
)

    # fit an initial model with scalar random effects to the hindi contrast group
m1 = fit(
    MixedModel,
    @formula(accuracy ~ 1 + Task * trained + (1|ID) + (1|item)),
    @subset(axb_data, :contrast == "hindi"),
    Bernoulli();
    contrasts=contr,
)


text:
include("imports.jl")
include("dice_metric.jl")
include("hausdorff_metric.jl")
include("utils.jl")


text:
@testset "135.candy.jl" begin
    @test candy([1, 0, 2]) == 5
    @test candy([1, 2, 2]) == 4
end


text:
using SupervisedLearning

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

msg("TrainingHistory: Basic functions")

_history = TrainingHistory()

function f(i, b; muh=10)
  @test b == "yo"
  @test muh == .3
  i
end

function f2(i, b; muh=10)
  @test b == "my"
  @test muh == .1
  i
end

@test_throws ArgumentError push!(_history, -1, :myf, f(10, "yo", muh = .3))

numbers = collect(0:2:100)
for i = numbers
  @test push!(_history, i, :myf, f(i + 1, "yo", muh = .3)) == i + 1
  @test @push!(_history, i, f2(i - 1, "my", muh = .1)) == i - 1
end

for (i, v) in enumerate(_history, :myf)
  @test in(i, numbers)
  @test i + 1 == v
end

for (i, v) in @enumerate(_history, f2(i - 1, "my", muh = .1))
  @test in(i, numbers)
  @test i - 1 == v
end

a1, a2 = get(_history, :myf)
@test typeof(a1) <: Vector && typeof(a2) <: Vector
@test length(a1) == length(a2) == length(numbers) == length(_history, :myf)
@test a1 + 1 == a2

a1, a2 = @get(_history, f2(i - 1, "my", muh = .1))
@test typeof(a1) <: Vector && typeof(a2) <: Vector
@test length(a1) == length(a2) == length(numbers) == @length(_history, f2(i - 1, "my", muh = .1))
@test a1 - 1 == a2

@test_throws ArgumentError push!(_history, 10, :myf, f(10, "yo", muh = .3))
@test_throws KeyError enumerate(_history, :sign)
@test_throws KeyError length(_history, :sign)

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

msg("TrainingHistory: Storing arbitrary types")

_history = TrainingHistory(Uint8)

for i = 1:100
  @test push!(_history, i % Uint8, :mystring, string("i=", i + 1)) == string("i=", i+1)
  @test push!(_history, i % Uint8, :myfloat, float(i + 1)) == float(i+1)
end

a1, a2 = get(_history, :mystring)
@test typeof(a1) <: Vector{Uint8}
for entry in a2
  @test typeof(entry) <: String
end

a1, a2 = get(_history, :myfloat)
@test typeof(a1) <: Vector{Uint8}
for entry in a2
  @test typeof(entry) <: FloatingPoint
end


text:
push!(LOAD_PATH,"../src/")
using DisjunctiveProgramming
using Documenter
makedocs(
         sitename = "DisjunctiveProgramming.jl",
         modules  = [DisjunctiveProgramming],
         pages=[
                "Home" => "index.md"
               ])
deploydocs(;
    repo="github.com/hdavid16/DisjunctiveProgramming.jl",
)


text:
struct KrylovKit_IRAM{Tv} <: AbstractEigenMethod{Tv}
    precond::Vector{Tv}
end

(::Type{<:KrylovKit_IRAM})(h::AbstractMatrix{Tv}) where {Tv} = KrylovKit_IRAM(rand(Tv, size(h, 2)))

function (d::Diagonalizer{<:KrylovKit_IRAM, Tv})(nev::Integer; precond = true, kw...) where {Tv}
    if isfinite(d.point)
        λs, ϕv, _ = KrylovKit.eigsolve(x -> d.lmap * x, d.method.precond, nev; kw...)
        λs .= 1 ./ λs .+ d.point
    else
        λs, ϕv, _ = KrylovKit.eigsolve(d.matrix, d.method.precond, nev,
                                       d.point > 0 ? :LR : :SR, Lanczos(kw...))
    end
    if precond
        d.method.precond .= zero(Tv)
        foreach(ϕ -> (d.method.precond .+= ϕ), ϕv)
    end
    return Eigen(real(λs), hcat(ϕv))
end


text:
push!(LOAD_PATH, "./src")
using DelimitedFiles
using Plots
using SetupDG
using UnPack
using UniformTriMesh
using PlotThemes
using ToeplitzMatrices
using Basis2DTri

macro Name(arg)
   string(arg)
end

gr(aspect_ratio=1,legend=false,colorbar=:right,framestyle=:none,xlims=(-1,1),ylims=(-1,1),
   markerstrokewidth=0,markersize=1,c=:haline)

N = 4
Nq_P = countlines("src/QuadratureData/quad_nodes_tri_N$(2*N).txt")
Nd = 5
K1D = 128
K = K1D*K1D*2
Np_F = 8
Nplot = 15
Mach = ".7"
delta = ".05"
extra_attribute = "varwidthsquared"

#T = 2.5
#T_arr = [0.5;1.0;1.5;2.0;2.5]
T_arr = [2.0]
for T in T_arr
    folder_name = "K$(K1D)NPF$(Np_F)N$(N)delta$(delta)Mach$(Mach)$(extra_attribute)"
    file_name = "array_end"
    if T == 0.5
        file_name = "array_1"
    elseif T == 1.0
        file_name = "array_2"
    elseif T == 1.5
        file_name = "array_3"
    elseif T == 2.0
        file_name = "array_4"
    elseif T == 2.5
        file_name = "array_end"
    end
    Q = readdlm("./examples/ES_Fourier-DG/$folder_name/$file_name.txt",'\t',Float64,'\n')
    xq = readdlm("./examples/ES_Fourier-DG/$folder_name/xq.txt",'\t',Float64,'\n')
    yq = readdlm("./examples/ES_Fourier-DG/$folder_name/yq.txt",'\t',Float64,'\n')

    xq = reshape(xq,Nq_P,Np_F,K)
    yq = reshape(yq,Nq_P,Np_F,K)

    h       = 2*pi/Np_F
    column  = [0; .5*(-1).^(1:Np_F-1).*cot.((1:Np_F-1)*h/2)]
    Dt      = Array{Float64,2}(Toeplitz(column,column[[1;Np_F:-1:2]]))

    rd = init_reference_tri(N);
    @unpack fv,Nfaces,r,s,VDM,V1,Vp,Dr,Ds,rf,sf,wf,nrJ,nsJ,rq,sq,wq,Vq,M,Pq,Vf,LIFT = rd
    rp, sp = Basis2DTri.equi_nodes_2D(Nplot)
    Vp = Basis2DTri.vandermonde_2D(N,rp,sp)/VDM
    Np = length(rp)

    "Mesh related variables"
    # Initialize 2D triangular mesh
    VX,VY,EToV = uniform_tri_mesh(K1D,K1D)
    @. VX = VX
    @. VY = VY
    md   = init_mesh((VX,VY),EToV,rd)
    # Intialize triangular prism
    VX   = repeat(VX,2)
    VY   = repeat(VY,2)
    VZ   = [2/Np_F*ones((K1D+1)*(K1D+1),1); ones((K1D+1)*(K1D+1),1)]
    @unpack rxJ,sxJ,ryJ,syJ = md

    JP = 1/K1D^2
    JF = 1/pi
    J = JF*JP
    Dr = JF*Dr
    Ds = JF*Ds
    Dt = JP*Dt
    VpPq = Vp*Pq

    Q = reshape(Q,Nq_P,Np_F,Nd,K)

    for f = [1]
        xq_f = zeros(Nq_P,K)
        yq_f = zeros(Nq_P,K)
        for k = 1:K
            @. xq_f[:,k] = xq[:,f,k]
            @. yq_f[:,k] = yq[:,f,k]
        end

        #rho_f = zeros(Nq_P,K)
        u_f = zeros(Nq_P,K)
        v_f = zeros(Nq_P,K)
        #w_f = zeros(Nq_P,K)
        vort_z_f = zeros(Nq_P,K)
        for k = 1:K
            # @. rho_f[:,k] = Q[:,f,1,k]
            @. u_f[:,k] = Q[:,f,2,k]./Q[:,f,1,k]
            @. v_f[:,k] = Q[:,f,3,k]./Q[:,f,1,k]
            #@. w_f[:,k] = Q[:,f,4,k]./Q[:,f,1,k]
        end

        v_f_modal = Pq*v_f
        u_f_modal = Pq*u_f
        #w_f_modal = Pq*w_f

        vort_z_f = Vq*(rxJ.*(Dr*v_f_modal)+sxJ.*(Ds*v_f_modal)) - Vq*(ryJ.*(Dr*u_f_modal)+syJ.*(Ds*u_f_modal))

        xq_f = VpPq*xq_f
        yq_f = VpPq*yq_f
        # rho_f = VpPq*rho_f
        # u_f = VpPq*u_f
        # v_f = VpPq*v_f
        # w_f = VpPq*w_f
        vort_z_f = VpPq*vort_z_f

        # rho_f = rho_f[:]
        # u_f = u_f[:]
        # v_f = v_f[:]
        # w_f = w_f[:]
        xq_f = xq_f[:]
        yq_f = yq_f[:]
        vort_z_f = vort_z_f[:]
        # absu_f = @. u_f^2+v_f^2+w_f^2

        # scatter(xq_f,yq_f,rho_f,zcolor=rho_f,camera=(0,90),axis=nothing)
        # png("./examples/ES_Fourier-DG/$folder_name/rho_f$(f)_T=$T.png")
        # scatter(xq_f,yq_f,u_f,zcolor=u_f,camera=(0,90),axis=nothing)
        # png("./examples/ES_Fourier-DG/$folder_name/u_f$(f)_T=$T.png")
        # scatter(xq_f,yq_f,v_f,zcolor=v_f,camera=(0,90),axis=nothing)
        # png("./examples/ES_Fourier-DG/$folder_name/v_f$(f)_T=$T.png")
        # scatter(xq_f,yq_f,w_f,zcolor=w_f,camera=(0,90),axis=nothing)
        # png("./examples/ES_Fourier-DG/$folder_name/w_f$(f)_T=$T.png")
        # scatter(xq_f,yq_f,absu_f,zcolor=absu_f,camera=(0,90),axis=nothing)
        # png("./examples/ES_Fourier-DG/$folder_name/speed_f$(f)_T=$T.png")
        scatter(xq_f,yq_f,vort_z_f,zcolor=vort_z_f,camera=(0,90),clims=(-0.001,0.001),axis=nothing)
        png("./examples/ES_Fourier-DG/$folder_name/vort_f$(f)_T=$T.png")
    end
end

max_stars_repo_name: yiminllin/ESDG-Fourier

max_issues_repo_name: yiminllin/ESDG-Fourier

max_forks_repo_name: yiminllin/ESDG-Fourier


text:
using Unitful
using YAML
using PyCall
using StaticArrays

module Calc
    include("Calculators/RateUncertainty.jl")
    include("Calculators/ThermoUncertainty.jl")
    include("Calculators/Thermo.jl")
    include("Calculators/Diffusion.jl")
    include("Calculators/Rate.jl")
    include("Calculators/Viscosity.jl")
    include("Calculators/Thermovec.jl")
    include("Calculators/Ratevec.jl")
end
module Spc
    include("Calculators/ThermoUncertainty.jl")
    include("Calculators/Thermo.jl")
    include("Calculators/Diffusion.jl")
    include("Species.jl")
end
module Rxn
    include("Calculators/RateUncertainty.jl")
    include("Calculators/ThermoUncertainty.jl")
    include("Calculators/Thermo.jl")
    include("Calculators/Diffusion.jl")
    include("Calculators/Rate.jl")
    include("Species.jl")
    include("Reaction.jl")
end
module Solv
    include("Calculators/Viscosity.jl")
    include("Solvent.jl")
end

const unitsdict = Dict()
const elementdict = Dict([1=>"H",6=>"C",8=>"O",7=>"N",17=>"Cl",16=>"S",18=>"Ar",10=>"Ne",2=>"He",
        15=>"P",9=>"F",35=>"Br",53=>"I",289=>"Fl"])

const allowedfcnlist = vcat(names(Calc),names(Spc),names(Rxn),names(Solv))

"""
Zhao et al 2003
"""
const mcgowanvolumes = Dict(["H"=> upreferred(8.71u"mL/mol").val, "He"=> upreferred(6.75u"mL/mol").val,
        "C"=> upreferred(16.35u"mL/mol").val, "N"=> upreferred(14.39u"mL/mol").val,
         "O"=> upreferred(12.43u"mL/mol").val, "F"=> upreferred(10.47u"mL/mol").val,
         "Ne"=> upreferred(8.51u"mL/mol").val,"Si"=> upreferred(26.83u"mL/mol").val, "P"=> upreferred(24.87u"mL/mol").val,
        "S"=> upreferred(22.91u"mL/mol").val, "Cl"=> upreferred(20.95u"mL/mol").val,
         "Ar"=> upreferred(18.99u"mL/mol").val,"Br"=> upreferred(26.21u"mL/mol").val,])

function parseentry(q,units=nothing)
    """
    parses single YAML lines into strings, numbers and arrays of Numbers in SI units
    """
    if isa(q,String)
        q = parsestring(q)
    end

    if isa(q,String)
        return q
    elseif isa(q,Quantity)
        return upreferred(q).val
    elseif isa(q,Number)
        if units == nothing
            return q
        else
            return upreferred(Quantity(q,units)).val
        end
    elseif isa(q,AbstractArray)
        if units == nothing
            if typeof(q[1]) <: AbstractArray
                return convert(Array,hcat(q...)')
            else
                return q
            end
        else
            if typeof(q[1]) <: AbstractArray
                return convert(Array,hcat(upreferred(Quantity(q,units)).val...)')
            else
                return upreferred(Quantity(q,units)).val
            end
        end
    elseif isa(q,Dict)
        return q
    else
        tq = typeof(q)
        throw(error("unable to parse $q of type $tq"))
    end

end
export parseentry

function parsestring(s)
    """
    Detects if a given string is actually a quantity input using metaprogramming
    the string is parsed by julia's parser and then analyzed to check if it is in
    the format of a quantity input
    """
    try
        ex = Meta.parse(s)
        if !isa(ex,Symbol) && !isa(ex,String) && length(ex.args) == 3 && ex.args[1] == :* &&
            (isa(ex.args[2],Number) || (ex.args[2].head == :vec && all([isa(x,Number) for x in ex.args[2].args]))) &&
            ex.args[3].args[1] == Symbol("@u_str") &&
            isa(ex.args[3].args[2],LineNumberNode) && ex.args[3].args[2].line == 1 &&
            ex.args[3].args[2].file == :none && isa(ex.args[3].args[3],String)
            return eval(ex)
        else
            return s
        end
    catch
        return s
    end
end
export parsestring

function getatomdictfromrdkit(mol)
    """
    retrives the number of each type of atom and the number of bonds of an rdkit molecule
    """
    atmD = Dict{String,Int64}()
    for atm in mol.GetAtoms()
        v = elementdict[atm.GetAtomicNum()]
        if v in keys(atmD)
            atmD[v] += 1
        else
            atmD[v] = 1
        end
    end
    nbonds = length(mol.GetBonds())
    return atmD,nbonds
end
export getatomdictfromrdkit

getatomdictsmiles(smiles) = getatomdictfromrdkit(Chem.AddHs(Chem.MolFromSmiles(smiles)))
export getatomdictsmiles
getatomdictinchi(inchi) = getatomdictfromrdkit(Chem.AddHs(Chem.MolFromInchi(inchi)))
export getatomdictinchi

function getspeciesradius(atomdict::Dict{String,Int64},nbonds::Int64)
    """
    estimates the McGowan radius by calculating the McGowan Volume
    from the number of each type of atom and the number of bonds
    as described in Zhao et al. 2003
    """
    Vtot = 0.0
    for (atm,n) in atomdict
        Vtot += n*mcgowanvolumes[atm]
    end
    Vtot -= nbonds*upreferred(6.56u"mL/mol").val
    V = Vtot/Na
    r = (0.75/Base.pi*V)^(1.0/3.0)
    return r
end
export getspeciesradius

function fcndict2obj(d::T,ymlunitsdict::Q) where {T,Q<:Any}
    """
    constructs an object from a dictionary by recursively constructing
    the objects in its fields
    the dictionary entry for "type" denotes the object being constructed
    which must be in the allowedfcnlist
    """
    strfcn = d["type"]
    fcn = Symbol(strfcn)
    @assert fcn in allowedfcnlist fcn
    kwexprs = Array{Expr,1}()
    for (key,val) in d
        if key == "type"
            continue
        end
        if isa(val,AbstractDict) && "type" in keys(val)
            val = fcndict2obj(val,ymlunitsdict)
        elseif isa(val,AbstractArray) && typeof(val).parameters[1] <: AbstractDict
            val = [fcndict2obj(x,ymlunitsdict) for x in val]
        else
            if (strfcn,key) in keys(unitsdict)
                if unitsdict[(strfcn,key)] in keys(ymlunitsdict)
                    val = parseentry(val,ymlunitsdict[unitsdict[(strfcn,key)]])
                else
                    val = parseentry(val)
                end
            else
                #@warn("assuming values for $strfcn,$key are in SI units")
                val = parseentry(val)
            end
        end
        kwexpr = Expr(:kw,Symbol(key),val)
        push!(kwexprs,kwexpr)
    end
    ex = Expr(:call,fcn,Expr(:parameters,kwexprs...))
    return eval(ex)
end
export fcndict2obj

"""
examines the input file and parses as appropriate
for chemkin (.inp files) with or without species dictionaries
calls rmgpy to convert it to yml and then parses the yml file
for rms (.rms, .yml, .rmg) it parses it directly
"""
function readinput(fname::String;spcdict::String="",output::String="chem.rms")
    extension = split(fname,".")[end]
    if extension in ["rms","yml","rmg"]
        return readinputyml(fname)
    elseif extension in ["inp"]
        if spcdict == ""
            convertchemkin2yml(fname,output=output)
        else
            convertchemkin2yml(fname,spcdictpath=spcdict,output=output)
        end
        return readinputyml(output)
    else
        throw(error("file extension $extension not understood, use .inp for chemkin files or .rms for rms yaml files"))
    end
end
export readinput

"""
parses a YAML input file into a dictionary containing
partitions of Species and Reaction objects
"""
function readinputyml(fname::String)
    D = YAML.load(open(fname))
    outdict = Dict()

    #Units
    ymlunitsdict=D["Units"] #for now don't use

    #Solvents
    if "Solvents" in keys(D)
        s = D["Solvents"]
        for sol in s
            sol["type"] = "Solvent"
        end
        solvents = [fcndict2obj(q,ymlunitsdict) for q in s]
        outdict["Solvents"] = solvents
    end

    #phases
    spcindex = 1
    spcdict = Dict()
    for p in D["Phases"]
        name = p["name"]
        spclist = Array{Species,1}()
        for d in p["Species"]
            d["index"] = spcindex
            spcindex += 1
            spcname = d["name"]
            #attempt to generate molecular information from rdkit if possible
            if !("atomnums" in keys(d)) || !("bondnum" in keys(d))
                if "smiles" in keys(d)
                    try
                        d["atomnums"],d["bondnum"] = getatomdictsmiles(d["smiles"])
                    catch
                        @warn("failed to generate molecular information from smiles for species $spcname")
                    end
                elseif "inchi" in keys(d)
                    try
                        d["atomnums"],d["bondnum"] = getatomdictinchi(d["inchi"])
                    catch
                        @warn("failed to generate molecular information from inchi for species $spcname")
                    end
                end
            end
            #if diffusion model is not present and we can calculate molecular radius and stokesdiffusivity
            if !("diffusion" in keys(d)) && "atomnums" in keys(d) && "bondnum" in keys(d)
                try
                    diffD = Dict()
                    diffD["type"] = "StokesDiffusivity"
                    diffD["r"] = getspeciesradius(d["atomnums"],d["bondnum"])
                    d["diffusion"] = diffD
                    if !("radius" in keys(d))
                        d["radius"] = diffD["r"]
                    end
                catch
                    @warn("failed to generate StokesDiffusivity model for species $spcname")
                end
            end

            spc = fcndict2obj(d,ymlunitsdict)
            push!(spclist,spc)
            spcdict[spc.name] = (spc,p["name"])
        end
        outdict[name] = Dict()
        outdict[name]["Species"] = spclist
        outdict[name]["Reactions"] = Array{ElementaryReaction,1}()
    end

    #reactions
    rxnindex = 1
    for rxn in D["Reactions"]
        phs = Array{String,1}()
        rxn["reactants"] = convert(Array{Any},rxn["reactants"])
        rxn["index"] = rxnindex
        rxnindex += 1
        for (i,item) in enumerate(rxn["reactants"])
            spc,ph = spcdict[item]
            push!(phs,ph)
            rxn["reactants"][i] = spc
        end
        rxn["products"] = convert(Array{Any},rxn["products"])
        for (i,item) in enumerate(rxn["products"])
            spc,ph = spcdict[item]
            push!(phs,ph)
            rxn["products"][i] = spc
        end
        rxn["reactants"] = [r for r in rxn["reactants"]]
        rxn["reactants"] = SVector(rxn["reactants"]...)
        rxn["products"] = [r for r in rxn["products"]]
        rxn["products"] = SVector(rxn["products"]...)
        rxn["reactantinds"] =  [r.index for r in rxn["reactants"]]
        rxn["reactantinds"] = SVector(rxn["reactantinds"]...)
        rxn["productinds"] =  [r.index for r in rxn["products"]]
        rxn["productinds"] = SVector(rxn["productinds"]...)
        if "efficiencies" in keys(rxn["kinetics"])
            rxn["kinetics"]["efficiencies"] = Dict{Int64,Float64}([spcdict[name][1].index=>val-1.0 for (name,val) in rxn["kinetics"]["efficiencies"]]) #in RMS we correct [M] rather than calculate it so we subtract 1
        end
        r = fcndict2obj(rxn,ymlunitsdict)
        unique!(phs)
        if length(phs) == 1
            push!(outdict[phs[1]]["Reactions"],r)
        else
            if Set(phs) in keys(outDict)
                push!(outdict[Set(phs)],r)
            else
                outdict[Set(phs)]=[r]
            end
        end
    end

    return outdict
end

export readinputyml


text:
event_register = Dict{AbstractString, Any}()

register_event(name::AbstractString, t::Any) = event_register[name] = t

find_event(name::AbstractString) = haskey(event_register, name) ? event_register[name] : nothing


text:
"""
    RootedSpanningForest{T<:LG.SimpleDiGraph{Int64}} <: AbstractGraphPointProcess{T}

Concrete type reprensenting a point process defined on the edges of a **connected** `graph` characterizing the uniform distribution on the directed [spanning forests](https://en.wikipedia.org/wiki/Spanning_tree) of `graph` rooted at `roots`.

It can be viewed as a the product distribution of the uniform distribution on the set of neighbors of each vertex conditioned on forming no cycles.

The object has two fields:

- `graph::LG.SimpleGraph{Int64}`
- `roots::Vector{Int64}`

**See also**

- Section 4.2 of [GuJeLi19](@cite).

# Example

A realization from a ``5\\times 5`` grid graph with `roots=[13]`.

![assets/rooted_spanning_tree.png](assets/rooted_spanning_tree.png)
"""
struct RootedSpanningForest{T<:LG.SimpleDiGraph{Int64}} <: AbstractGraphPointProcess{T}
    "Graph"
    graph::LG.SimpleGraph{Int64}
    "Roots"
    roots::Vector{Int64}
end

function Base.show(io::IO, pp::RootedSpanningForest{T}) where {T}
    print(io, "RootedSpanningForest{$T}\n- graph = $(pp.graph)\n- roots = $(pp.roots)")
end

"""
    RootedSpanningForest(
        graph::LG.SimpleGraph{T},
        roots::Union{Nothing,T,AbstractVector{T},AbstractSet{T}}=nothing
    ) where {T<:Int}

Construct a [`PRS.RootedSpanningForest`](@ref) model on a **connected** `graph`, rooted at `roots`.
If `roots === nothing`, a random vertex is selected uniformly at random among `LG.vertices(g)` and considered as `roots`.

```jldoctest; output = true
using PartialRejectionSampling
using LightGraphs; const LG = LightGraphs

g, roots = LG.grid([5, 5]), [1, 2, 3]

rsf = PRS.RootedSpanningForest(g, roots)

# output

RootedSpanningForest{LightGraphs.SimpleGraphs.SimpleDiGraph{Int64}}
- graph = {25, 40} undirected simple Int64 graph
- roots = [1, 2, 3]
```
"""
function RootedSpanningForest(
    graph::LG.SimpleGraph{T},
    roots::Union{Nothing,T,AbstractVector{T},AbstractSet{T}}=nothing
) where {T<:Int}
    roots_ = unique(isnothing(roots) ? rand(LG.vertices(graph)) : roots)
    if issubset(roots_, LG.vertices(graph))
        if LG.is_connected(graph)
            return RootedSpanningForest{LG.SimpleDiGraph{T}}(graph, roots_)
        else
            throw(DomainError(graph, "The graph must be connected"))
        end
    else
        throw(DomainError(roots, "roots must be contained in vertices(graph)"))
    end
end

"""
    generate_sample(
        [rng::Random.AbstractRNG,]
        pp::RootedSpanningForest
    )

Generate an exact sample from the [`PRS.RootedSpanningForest`](@ref).

Default sampler is [`PRS.generate_sample_prs`](@ref).
"""
function generate_sample(
    rng::Random.AbstractRNG,
    pp::RootedSpanningForest
)
    return generate_sample_prs(rng, pp)
end

# generate_sample(pp::RootedSpanningForest) = generate_sample(Random.default_rng(), pp)

"""
    generate_sample_prs(
        [rng::Random.AbstractRNG,]
        pp::RootedSpanningForest{T}
    )::T where {T<:LG.SimpleDiGraph{Int64}}

Generate a rooted spanning forest of `pp.`graph` with prescribed `pp.roots`, uniformly at random among all rooted spanning forests rooted at `pp.roots`, using Partial Rejection Sampling (PRS).

**See also**
- Section 4.2 of [GuJeLi19](@cite).

# Example

An illustration of the procedure on a ``5\\times 5`` grid graph, with `roots=[13]`

- red: variables involved in constraints violation (edges forming a cycle)
- orange: variables resampled (edges originating from a red/orange node)

![assets/rooted_spanning_tree_prs.gif](assets/rooted_spanning_tree_prs.gif)
"""
function generate_sample_prs(
    rng::Random.AbstractRNG,
    pp::RootedSpanningForest{T}
)::T where {T<:LG.SimpleDiGraph{Int64}}
    return _generate_sample_rooted_spanning_forest_prs(rng, pp.graph, pp.roots)
end

# generate_sample_prs(pp::RootedSpanningForest) = generate_sample_prs(Random.default_rng(), pp)

"""
    _generate_sample_rooted_spanning_forest_prs(
        rng::Random.AbstractRNG,
        graph::LG.SimpleGraph{T},
        roots
    )::LG.SimpleDiGraph{T} where {T}

Generate a rooted spanning forest from a **connected** `graph` with prescribed `roots`, uniformly at random among all rooted spanning forests rooted at `roots`, using Partial Rejection Sampling.

**See also**
- Section 4.2 of [GuJeLi19](@cite).
"""
function _generate_sample_rooted_spanning_forest_prs(
    rng::Random.AbstractRNG,
    graph::LG.SimpleGraph{T},
    roots
)::LG.SimpleDiGraph{T} where {T}
    g = random_neighbor_assignment(rng, graph, roots)
    while true
        vertices_in_cycles = Set(Iterators.flatten(LG.simplecycles(g)))
        isempty(vertices_in_cycles) && break
        # Resample the successor of vertices involved in cycles
        for v in vertices_in_cycles
            # Remove current edge (v, w)
            w = LG.neighbors(g, v)[1]
            LG.rem_edge!(g, v, w)
            # Resample the successor w of v and add edge (v, w)
            w = rand(rng, LG.neighbors(graph, v))
            LG.add_edge!(g, v, w)
        end
    end
    return g
end

function _generate_sample_rooted_spanning_forest_prs(
    graph::LG.SimpleGraph{T},
    roots
)::LG.SimpleDiGraph{T} where {T}
    return _generate_sample_rooted_spanning_forest_prs(Random.default_rng(), graph, roots)
end

@doc raw"""
    _generate_sample_rooted_spanning_forest_wilson(
        [rng::Random.AbstractRNG,]
        graph::LG.SimpleGraph{T},
        roots
    )

Generate a spanning forest of a **connected** `graph` with prescribed `roots`, uniformly at random among all rooted spanning forests rooted at `roots`, using [Wilson's algorithm](https://en.wikipedia.org/wiki/Loop-erased_random_walk#Uniform_spanning_tree).

**See also**
- [RandomForests.jl implementation](https://gricad-gitlab.univ-grenoble-alpes.fr/barthesi/RandomForests.jl/-/blob/master/src/random_spanning_tree.jl).
"""
function _generate_sample_rooted_spanning_forest_wilson(
    rng::Random.AbstractRNG,
    graph::LG.SimpleGraph{T},
    roots
)::LG.SimpleDiGraph{T} where {T}
    in_tree = falses(LG.nv(graph))
    for r in roots
        @inbounds in_tree[r] = true
    end

    successors = zeros(T, LG.nv(graph))
    for i in LG.vertices(graph)
        # Run a natural random walk on g from i until the walk hits a vertex in tree
        u = i
        while !in_tree[u]
            successors[u] = rand(rng, LG.neighbors(graph, u))
            u = successors[u]
        end
        # Erase loops formed during the walk
        u = i
        while !in_tree[u]
            in_tree[u] = true
            u = successors[u]
        end
    end

    return directed_forest_from_successors(successors)
end

function _generate_sample_rooted_spanning_forest_wilson(
    graph::LG.SimpleGraph{T},
    roots
)::LG.SimpleDiGraph{T} where {T}
    return _generate_sample_rooted_spanning_forest_wilson(Random.default_rng(), graph, roots)
end

"""
    directed_forest_from_successors(successors)::LG.SimpleDiGraph

Return a directed graph of size `length(successors)` with edges `i => successors[i]`.
"""
function directed_forest_from_successors(successors)::LG.SimpleDiGraph
    graph = LG.SimpleDiGraph(length(successors))
    for (i, j) in enumerate(successors)
        LG.add_edge!(graph, i, j)
    end
    return graph
end


text:
module DynamicPPL

using AbstractMCMC: AbstractSampler, AbstractChains, AbstractModel
using Distributions
using Bijectors
using MacroTools
import ZygoteRules

import Base: string,
             Symbol,
             ==,
             hash,
             in,
             getindex,
             setindex!,
             push!,
             show,
             isempty,
             empty!,
             getproperty,
             setproperty!,
             keys,
             haskey

# VarInfo
export  VarName,
        AbstractVarInfo,
        VarInfo,
        UntypedVarInfo,
        getlogp,
        setlogp!,
        acclogp!,
        resetlogp!,
        get_num_produce,
        set_num_produce!,
        reset_num_produce!,
        increment_num_produce!,
        set_retained_vns_del_by_spl!,
        is_flagged,
        unset_flag!,
        setgid!,
        updategid!,
        setorder!,
        istrans,
        link!,
        invlink!,
        tonamedtuple,
# Compiler
        ModelGen,
        @model,
        @varname,
        @varinfo,
        @logpdf,
        @sampler,
# Utilities
        vectorize,
        reconstruct,
        reconstruct!,
        Sample,
        Chain,
        init,
        vectorize,
        set_resume!,
# Model
        ModelGen,
        Model,
        getmissings,
        getargnames,
        getdefaults,
        getgenerator,
        runmodel!,
# Samplers
        Sampler,
        SampleFromPrior,
        SampleFromUniform,
# Contexts
        DefaultContext,
        LikelihoodContext,
        PriorContext,
        MiniBatchContext,
        assume,
        dot_assume,
        observer,
        dot_observe,
        tilde,
        dot_tilde,
# Pseudo distributions
        NamedDist,
        NoDist,
# Prob macros
        @prob_str,
        @logprob_str

const DEBUG = Bool(parse(Int, get(ENV, "DEBUG_DYNAMICPPL", "0")))

# Used here and overloaded in Turing
function getspace end

# Necessary forward declarations
abstract type AbstractVarInfo end
abstract type AbstractContext end

include("utils.jl")
include("selector.jl")
include("model.jl")
include("sampler.jl")
include("varname.jl")
include("distribution_wrappers.jl")
include("contexts.jl")
include("varinfo.jl")
include("context_implementations.jl")
include("compiler.jl")
include("prob_macro.jl")
include("ad.jl")

end # module


text:
include("../examples/paperNumerics/example4_2.jl");
include("../examples/paperNumerics/example4_3.jl");
include("../examples/paperNumerics/example4_4.jl");
include("../examples/paperNumerics/example4_5.jl");
include("../examples/paperNumerics/toymodel.jl");

include("../examples/meNumerics/example4_2.jl");
include("../examples/meNumerics/example4_3.jl");
include("../examples/meNumerics/example4_5.jl");
include("../examples/meNumerics/discontinuities.jl");
include("../examples/meNumerics/oneDwave.jl");

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
import Base.Broadcast: Broadcasted, AbstractArrayStyle, dotview, materialize

# """
#     mutable struct Scalar{T} <: AbstractArray{T,0}
#
# A 0-dimensional (scalar) object, much like a `RefValue{T}`, but one which is an
# abstract array and assignable under broadcasting operations.
# """
mutable struct Scalar{T} <: AbstractArray{T,0}
    x::T
    Scalar{T}()  where {T} = new()
    Scalar{T}(x) where {T} = new(x)
end
Scalar(s::T) where {T} = Scalar{T}(s)
Base.isassigned(s::Scalar) = isdefined(s, :x)

# AbstractArray interfaces
Base.size(s::Scalar) = ()
@propagate_inbounds Base.getindex(s::Scalar) = s.x
@propagate_inbounds Base.getindex(s::Scalar, ::CartesianIndex{0}) = s.x
@propagate_inbounds Base.setindex!(s::Scalar, x) = (s.x = x; s)
@propagate_inbounds Base.setindex!(s::Scalar, x, ::CartesianIndex{0}) = (s.x = x; s)
Base.similar(s::Scalar{T}) where {T} = Scalar{T}()
Base.similar(s::Scalar, ::Type{T}) where {T} = Scalar{T}()
Base.similar(s::Scalar{T}, ::Tuple{}) where {T} = Scalar{T}()
Base.similar(s::Scalar, ::Type{T}, ::Tuple{}) where {T} = Scalar{T}()

# Iteration on a scalar
Base.iterate(s::Scalar) = (s.x, nothing)
Base.iterate(s::Scalar, ::Nothing) = nothing

# Broadcasting extensions
@propagate_inbounds dotview(A::Scalar, ::CartesianIndex{0}) = A
@propagate_inbounds dotview(A::Scalar, ::CartesianIndices{0,Tuple{}}) = A
@propagate_inbounds Base.copyto!(dest::Scalar, bc::Broadcasted{<:AbstractArrayStyle{0}}) =
    dest[] = materialize(bc)


text:
module CVMAbstractModel

require("CommonTools.jl")

using DataTypes
using Corpora.Token
using CommonTools.CVMIntermediates
using CommonTools: mult_by_0!, sigmoid, loopdot, get_rand
using CommonTools: tanh, square!, l2sq, l2sq_distance

export AbstractModel
export AbstractParameters
export accumulate_gradients_batch!
export accumulate_closeness_gradients_batch!
export objective_gradients_batch
export regularize_l2!

abstract AbstractModel;
abstract AbstractParameters;

# Whenever we create a new model, these two methods need to be implemented
# (compose! and store_gradients!)

# Given the gradient of the error function w.r.t. the phrase representation,
# backpropagate this gradient to the parameters of the model
function store_gradients!(model::AbstractModel,
                    pre_grad::Vector{Float},
                    s::Vector{Token}, intermediate::CVMIntermediates;
                    intermediates_precomputed::Bool=false)
    error("This is an abstract method");
end

# Given a list of tokens, compose a phrase vector representation
function compose!(model::AbstractModel,
                                    s::Vector{Token},
                                    intermediate::CVMIntermediates)
    error("This is an abstract method");
end

# Adds l2 regularization to the given parameters
function regularize_l2!(model::AbstractModel, params::AbstractParameters;
                            sparse_updates::Bool=false)
    acc_reg = const_zero;
    if model.regularizer > const_zero
        if sparse_updates
            # Only add regularization to the gradients that are not 0
            for i=1:length(params.wordvecs)
                if params.gradients[i] != 0
                    params.gradients[i] += params.wordvecs[i] * model.regularizer
                    acc_reg += model.regularizer * params.wordvecs[i]^2
                end
            end
        else
            # Add l2 regularization
            Base.LinAlg.BLAS.axpy!(model.regularizer, params.wordvecs, params.gradients)
            # Compute l2 error
            acc_reg += l2sq(params.wordvecs) * model.regularizer
        end
    end
    return acc_reg*0.5f0;
end

# Accumulates the gradients of |s1-s2|^2
function accumulate_closeness_gradients_batch!(model::AbstractModel,
                s1::Vector{Token}, s2::Vector{Token},
                intermediates::Vector{CVMIntermediates}; scaling=1.0)
    err = 0.0
    # Compose vector representation sentence1 and sentence2
    o1 = compose!(model, s1, intermediates[1])
    o2 = compose!(model, s2, intermediates[2])

    # Cache distance of the pair of aligned sentences (valid pair)
    dist_o1_o2 = l2sq_distance(o1, o2)
    # Update error
    err += scaling * dist_o1_o2
    # Calculate the gradient of |o1-o2|^2 w.r.t. o1
    pre_grad_closeness = (o1-o2)
    # Propagate that gradient to the individual word representations
    # of sentence1
    Base.LinAlg.BLAS.scal!(model.wordvecdim,
                                    scaling,
                                    pre_grad_closeness, 1);
    store_gradients!(model, pre_grad_closeness, s1, intermediates[1])
    # The gradient of |o1-o2|^2 w.r.t. o2 is just the negative o1 gradient
    # => invert the sign
    Base.LinAlg.BLAS.scal!(model.wordvecdim,
                                    -1.0,
                                    pre_grad_closeness, 1);
    # Propagage the gradient down to the word representations of sentence 2
    store_gradients!(model, pre_grad_closeness, s2, intermediates[2])

    return err;
end

# Calculates the error of the hinge loss and the closeness
function objective_gradients_batch(model::AbstractModel,
                s1::Vector{Token}, s2::Vector{Token},
                sns::Vector{Vector{Token}},
                intermediates::Vector{CVMIntermediates};
                scaling::Float=const_one)

    # Calculate the representations of sentence1 and sentence2
    # In case of the inclusion criterion, sentence1 is the outer window
    # and sentence 2 is the inner window.
    o1 = compose!(model, s1, intermediates[1])
    o2 = compose!(model, s2, intermediates[2])

    # Cache distance
    dist_o1_o2 = l2sq_distance(o1, o2)
    # Update error
    err_closeness = dist_o1_o2 * 0.5

    # Calculate the contrastive error (hinge loss) for each noise sample
    err_contrastive = 0.0
    inv_len_sns = 1.0/size(sns,1)
    cn_scaling = scaling * model.gamma * inv_len_sns
    for sn in sns
        # Compose the vector representation of the noise sample
        on = compose!(model, sn, intermediates[3])

        # Compute the contastive part of the hinge loss
        contrastive = 0.5*(dist_o1_o2 - l2sq_distance(o1, on))

        # Calculate the hinge loss
        contrastive_err = max(0.0, model.margin + contrastive)

        # If hinge loss greater than 0: Update the accumulative error
        if  contrastive_err > 0
            err_contrastive += contrastive_err * cn_scaling;
        end
    end

    return err_contrastive + err_closeness;
end

# Accumulates the gradients of the hinge loss and the closeness function
function accumulate_gradients_batch!(model::AbstractModel,
                s1::Vector{Token}, s2::Vector{Token},
                sns::Vector{Vector{Token}},
                intermediates::Vector{CVMIntermediates};
                scaling::Float=const_one)

    # Calculate the representations of sentence1 and sentence2
    # In case of the inclusion criterion, sentence1 is the outer window
    # and sentence 2 is the inner window.
    o1 = compose!(model, s1, intermediates[1])
    o2 = compose!(model, s2, intermediates[2])

    # Cache distance of the pair of aligned sentences (valid pair)
    dist_o1_o2 = l2sq_distance(o1, o2)
    err_closeness = dist_o1_o2 * 0.5

    # Calculate the gradient of |o1-o2|^2 w.r.t. o1
    pre_grad_closeness = (o1-o2)
    Base.LinAlg.BLAS.scal!(model.wordvecdim,
                                    scaling,
                                    pre_grad_closeness, 1);
    # Propagage the gradient down to the word representations of sentence 2
    store_gradients!(model, pre_grad_closeness, s1, intermediates[1])
    # The gradient of |o1-o2|^2 w.r.t. o2 is just the negative o1 gradient
    # => invert the sign
    Base.LinAlg.BLAS.scal!(model.wordvecdim,
                                    -1.0,
                                    pre_grad_closeness, 1);
    # Propagage the gradient down to the word representations of sentence 2
    store_gradients!(model, pre_grad_closeness, s2, intermediates[2])

    err_contrastive = 0.0
    # Scale the number contribution of each individual noise sample
    # by 1/their number
    inv_len_sns = 1.0/size(sns,1)
    cn_scaling = scaling * model.gamma * inv_len_sns
    # Apply the scaling to the |o1-o2|^2 w.r.t. o2 gradient
    Base.LinAlg.BLAS.scal!(model.wordvecdim,
                                    model.gamma * inv_len_sns,
                                    pre_grad_closeness, 1);
    for sn in sns
        on = compose!(model, sn, intermediates[3])
        contrastive = 0.5*(dist_o1_o2 - l2sq_distance(o1, on))

        # Calculate the hinge loss
        contrastive_err = max(0.0, model.margin + contrastive)

        # If hinge loss greater than 0: Minimize the contrastive distances
        if  contrastive_err > 0

            # Compute gradients of |o1,o2|^2-|o1,on|^2 w.r.t. on
            pre_gradn = o1-on
            Base.LinAlg.BLAS.scal!(model.wordvecdim,
                                    cn_scaling,
                                    pre_gradn, 1);
            ### Dangerous! Overwriting "on" and reuising it is potentially
            # dangerous!!
            # Original, clean version
            #pre_grad1 = on-o2
            # hacky, high-speed version reusing "on"
            pre_grad1 = on
            Base.LinAlg.BLAS.axpy!(-1.0, o2, pre_grad1)
            Base.LinAlg.BLAS.scal!(model.wordvecdim,
                                    cn_scaling,
                                    pre_grad1, 1);

            # Add the gradients up in the gradient field of each token
            ## Hinge objective
            store_gradients!(model, pre_grad1, s1, intermediates[1],
                                    intermediates_precomputed=true)
            # Could optimize this by pulling it out of the loop
            store_gradients!(model, pre_grad_closeness, s2, intermediates[2],
                                    intermediates_precomputed=true)
            # Propagage the gradients down to the word representations
            # of the noise sentence
            store_gradients!(model, pre_gradn, sn, intermediates[3])

            # Update accumulative error
            err_contrastive += contrastive_err * model.gamma;
        end
    end

    return err_closeness, err_contrastive;
end

end

max_stars_repo_name: ogh/binclusion

max_issues_repo_name: ogh/binclusion

max_forks_repo_name: ogh/binclusion


text:
module Js

import Logging, HTTP, Reexport

Reexport.@reexport using Genie
Reexport.@reexport using Genie.Renderer

const JS_FILE_EXT   = ".jl"
const TEMPLATE_EXT  = ".jl.js"

const SUPPORTED_JS_OUTPUT_FILE_FORMATS = [TEMPLATE_EXT]

const JSString = String

const NBSP_REPLACEMENT = (" "=>"!!nbsp;;")

export js

function get_template(path::String; context::Module = @__MODULE__, vars...) :: Function
  orig_path = path

  path, extension = Genie.Renderer.view_file_info(path, SUPPORTED_JS_OUTPUT_FILE_FORMATS)

  if ! isfile(path)
    error_message = length(SUPPORTED_JS_OUTPUT_FILE_FORMATS) == 1 ?
                    """JS file "$orig_path$(SUPPORTED_JS_OUTPUT_FILE_FORMATS[1])" does not exist""" :
                    """JS file "$orig_path" with extensions $SUPPORTED_JS_OUTPUT_FILE_FORMATS does not exist"""
    error(error_message)
  end

  f_name = Genie.Renderer.function_name(path) |> Symbol
  mod_name = Genie.Renderer.m_name(path) * ".jl"
  f_path = joinpath(Genie.config.path_build, Genie.Renderer.BUILD_NAME, mod_name)
  f_stale = Genie.Renderer.build_is_stale(path, f_path)

  if f_stale || ! isdefined(context, f_name)
    f_stale && Genie.Renderer.build_module(to_js(read(path, String), extension = extension), path, mod_name)

    return Base.include(context, joinpath(Genie.config.path_build, Genie.Renderer.BUILD_NAME, mod_name))
  end

  getfield(context, f_name)
end

function to_js(data::String; prepend = "\n", extension = TEMPLATE_EXT) :: String
  output = string("function $(Genie.Renderer.function_name(data))($(Genie.Renderer.injectkwvars())) :: String \n", prepend)

  output *= if extension == TEMPLATE_EXT
          "\"\"\"
          $data
          \"\"\""
  elseif extension == JS_FILE_EXT
    data
  else
    error("Unsuported template extension $extension")
  end

  string(output, "\nend \n")
end

function render(data::String; context::Module = @__MODULE__, vars...) :: Function
  Genie.Renderer.registervars(; context = context, vars...)

  data_hash = hash(data)
  path = "Genie_" * string(data_hash)

  func_name = Genie.Renderer.function_name(string(data_hash)) |> Symbol
  mod_name = Genie.Renderer.m_name(path) * ".jl"
  f_path = joinpath(Genie.config.path_build, Genie.Renderer.BUILD_NAME, mod_name)
  f_stale = Genie.Renderer.build_is_stale(f_path, f_path)

  if f_stale || ! isdefined(context, func_name)
    f_stale && Genie.Renderer.build_module(to_js(data), path, mod_name)

    return Base.include(context, joinpath(Genie.config.path_build, Genie.Renderer.BUILD_NAME, mod_name))
  end

  getfield(context, func_name)
end

function render(viewfile::Genie.Renderer.FilePath; context::Module = @__MODULE__, vars...) :: Function
  Genie.Renderer.registervars(; context = context, vars...)

  get_template(string(viewfile); partial = false, context = context, vars...)
end

function render(::Type{MIME"application/javascript"}, data::String; context::Module = @__MODULE__, vars...) :: Genie.Renderer.WebRenderable
  try
    Genie.Renderer.WebRenderable(render(data; context = context, vars...), :javascript)
  catch ex
    isa(ex, KeyError) && Genie.Renderer.changebuilds() # it's a view error so don't reuse them
    rethrow(ex)
  end
end

function render(::Type{MIME"application/javascript"}, viewfile::Genie.Renderer.FilePath; context::Module = @__MODULE__, vars...) :: Genie.Renderer.WebRenderable
  try
    Genie.Renderer.WebRenderable(render(viewfile; context = context, vars...), :javascript)
  catch ex
    isa(ex, KeyError) && Genie.Renderer.changebuilds() # it's a view error so don't reuse them
    rethrow(ex)
  end
end

function js(data::String; context::Module = @__MODULE__, status::Int = 200,
            headers::Genie.Renderer.HTTPHeaders = Genie.Renderer.HTTPHeaders("Content-Type" => Genie.Renderer.CONTENT_TYPES[:javascript]),
            forceparse::Bool = false, noparse::Bool = false, vars...) :: Genie.Renderer.HTTP.Response
  if (occursin(raw"$", data) || occursin("<%", data) || forceparse) && ! noparse
    Genie.Renderer.WebRenderable(render(MIME"application/javascript", data; context = context, vars...), :javascript, status, headers) |> Genie.Renderer.respond
  else
    js!(data; status, headers) |> Genie.Renderer.respond
  end
end

function js!(data::S; status::Int = 200,
              headers::Genie.Renderer.HTTPHeaders = Genie.Renderer.HTTPHeaders("Content-Type" => Genie.Renderer.CONTENT_TYPES[:javascript]))::Genie.Renderer.HTTP.Response where {S<:AbstractString}
  Genie.Renderer.WebRenderable(data, :javascript, status, headers) |> Genie.Renderer.respond
end

function js(viewfile::Genie.Renderer.FilePath; context::Module = @__MODULE__, status::Int = 200,
            headers::Genie.Renderer.HTTPHeaders = Genie.Renderer.HTTPHeaders("Content-Type" => Genie.Renderer.CONTENT_TYPES[:javascript]), vars...) :: Genie.Renderer.HTTP.Response
  Genie.Renderer.WebRenderable(render(MIME"application/javascript", viewfile; context = context, vars...), :javascript, status, headers) |> Genie.Renderer.respond
end

### === ###
### EXCEPTIONS ###

function Genie.Router.error(error_message::String, ::Type{MIME"application/javascript"}, ::Val{500}; error_info::String = "") :: HTTP.Response
  HTTP.Response(Dict("error" => "500 Internal Error - $error_message", "info" => error_info), status = 500) |> js
end

function Genie.Router.error(error_message::String, ::Type{MIME"application/javascript"}, ::Val{404}; error_info::String = "") :: HTTP.Response
  HTTP.Response(Dict("error" => "404 Not Found - $error_message", "info" => error_info), status = 404) |> js
end

function Genie.Router.error(error_code::Int, error_message::String, ::Type{MIME"application/javascript"}; error_info::String = "") :: HTTP.Response
  HTTP.Response(Dict("error" => "$error_code Error - $error_message", "info" => error_info), status = error_code) |> js
end

end


text:
export verify_facet

"""
    verify_facet(
        d :: Int64,
        facet :: BellScenario.BellGame,
        aff_ind_vertices :: Vector{Matrix{Int64}}
    ) :: Bool

Returns `true` if the affinely independent set of vertices
`aff_ind_vertices` verifies that the `facet` is indeed a tight Bell inequality
of the ``\\matcal{C}_d^{X \\to Y}`` signaling polytope.
Input `d` must align with the considered facet.
This method checks that all provided vertices are:
* extreme points of the ``\\mathcal{C}_d^{X \\to Y}`` signaling polytope
* satisfy the facet inequality with equality
* affinely independent

!!! note "Note"
    This method does not check that all signaling polytope vertices satisfy the
    facet inequality. For the cataloged Bell inequalities, this requirement can be
    analytically shown to hold in general. To verify the tightness of the
    Bell inequality, it therefore remains to find a sufficient number of affinely
    independent vertices that saturate the inequality.
"""
function verify_facet(d :: Int64, facet :: BellScenario.BellGame, aff_ind_vertices :: Vector{Matrix{Int64}}) :: Bool
    (Y, X) = size(facet)
    polytope_dim = X*(Y-1)

    _verify_vertex(v) = (
        rank(v) == d
        && facet[:]'*v[:] == facet.β
        && is_deterministic(v)
    )

    all_valid_vertices = all(_verify_vertex, aff_ind_vertices)

    affine_independence = (
        length(aff_ind_vertices) == polytope_dim
        && LocalPolytope.dimension(aff_ind_vertices) == polytope_dim - 1
    )

    (all_valid_vertices && affine_independence)
end


text:
export jetcolormap, asimage, asimagesc, asimagescrgb, blocksvisu, pad, image2array, poly2mask, embedvisu
export jetcolors, jetcolorants
export aslogimage

function jetcolormap(n)
    step(m, i) = i>=m ? 1.0 : 0.0
    grad(m, i) = (i-m)/(n/4) * step(m, i)

    a = n/8
    f(m, i) = grad(m+3a, i) - grad(m+a, i) - grad(m-a, i) + grad(m-3a, i)

    r = zeros(3,n)
    r[1, :] = [f(3n/4, i) for i in 1:n]
    r[2, :] = [f(n/2, i) for i in 1:n]
    r[3, :] = [f(n/4, i) for i in 1:n]
    clamp(r, 0f0, 1f0)
end
function jetcolors(a,mi,ma, n = 101)
    j = jetcolormap(n)
    @p collect a | minus mi | divby (ma-mi) | clamp 0. 1. | times 99.9 | plus 1 | round Int _ | part j _
end
jetcolors(a, args...) = jetcolors(a, minimum(a), maximum(a)+0.01)
jetcolorants(a...) = mapvec(jetcolors(a...),x->RGB(x[1],x[2],x[3]))

function asimagescrgb(a, norm = true)
    a = asfloat32(a)
    cm = jetcolormap(256)
    r = Array(Float32, size(a, 1), size(a, 2), 3)
    normf = norm ? norm01 : identity
    b = round(Int, normf(a)*255) .+ 1
    r[:,:,1] = cm[1, b[:]]
    r[:,:,2] = cm[2, b[:]]
    r[:,:,3] = cm[3, b[:]]
    r
end

if isinstalled("Images")
    function asimage{T}(a::Array{T,2})
        grayim(a')
    end

    function asimage{T}(a::Array{T,3})
        colorim(permutedims(a,[3,2,1]))
    end
    asimagesc = asimage*asimagescrgb
end

function blocksvisu(a, padding = 0)
    a = unstack(a)
    n = len(a)
    typ = eltype(fst(a))
    paddingf(a...) = ones(typ,a...)*padding
    padsize(a) = ceil(Int,a/4)
    try
        a = @p map a reshape | unstack
    end
    @p mapvec a size | uniq | len | isequal 1 | assert
    if ndims(fst(a))==3
        return @p map (1:3) (i->@p map a at i | blocksvisu padding) | stack
    end
    z = @p paddingf padsize(size(fst(a),1)) size(fst(a),2)
    a = @p partsoflen a ceil(Int,sqrt(n)) | map riffle z | map col | map flatten
    z = @p paddingf size(fst(a),1) padsize(size(fst(a),2))
    a[end] = @p pad a[end] siz(fst(a)) padding
    r = @p riffle a z | row | flatten
end

function embedvisu{T<:Number,N}(a::Matrix, patches::DenseArray{T,N}, s = 2000)
    mi = minimum(a,2)
    ma = maximum(a,2)
    rm = rangem(patches)
    rn = rangen(patches)
    dim = ndims(patches) == 3 ? 1 : 3
    t(a) = round(Int, (a-mi)./(ma-mi) * s)
    r = zeros(Float32, s+sizem(patches),s+sizen(patches),dim)
    for i = 1:len(patches)
        x = @p at a i | t
        p = @p at patches i
        r[x[1]+rm, x[2]+rn,:] = p
    end
    dim == 1 ? squeeze(r,3) : r
end

function pad(a, siz, value = 0)
    r = value * onessiz(siz, eltype(a))
    if ndims(r) == 1
        r[1:length(a)] = a
    elseif ndims(r) == 2
        r[1:size(a,1),1:size(a,2)] = a
    elseif ndims(r) == 3
        r[1:size(a,1),1:size(a,2),1:size(a,3)] = a
    elseif ndims(r) == 4
        r[1:size(a,1),1:size(a,2),1:size(a,3),1:size(a,4)] = a
    elseif ndims(r) == 5
        r[1:size(a,1),1:size(a,2),1:size(a,3),1:size(a,4),1:size(a,5)] = a
    else
        error("ndims not supperted")
    end
    r
end

image2array(img) = @p map Any[:r,:g,:b] (i->map(y->y.(i), img.data)) | stack

aslogimage(a) = @p clamp a 0 Inf | plus 1 | log | asimagesc


text:
include("mvrandn_tests.jl")
include("timeseries_distributions_tests.jl")
include("inverse_logic_tests.jl")
include("constrained_distributions_tests.jl")
include("optimization_tests.jl")


text:
@testset "call_at combinator" begin

    @gen (grad) function foo((grad)(x::Float64))
        return x + @trace(normal(x, 1), :y)
    end

    at = call_at(foo, Int)

    @testset "assess" begin
        y = 1.234
        choices = choicemap()
        choices[3 => :y] = y
        (weight, value) = assess(at, (0.4, 3), choices)
        @test isapprox(weight, logpdf(normal, y, 0.4, 1))
        @test isapprox(value, 0.4 + y)
    end

    @testset "propose" begin
        (choices, weight, value) = propose(at, (0.4, 3))
        y = choices[3 => :y]
        @test isapprox(weight, logpdf(normal, y, 0.4, 1))
        @test length(collect(get_values_shallow(choices))) == 0
        @test length(collect(get_submaps_shallow(choices))) == 1
    end

    @testset "generate" begin

        # without constraint
        (trace, weight) = generate(at, (0.4, 3), EmptyChoiceMap())
        choices = get_choices(trace)
        @test weight == 0.
        y = choices[3 => :y]
        @test get_retval(trace) == 0.4 + y
        @test length(collect(get_values_shallow(choices))) == 0
        @test length(collect(get_submaps_shallow(choices))) == 1

        # with constraints
        y = 1.234
        constraints = choicemap()
        constraints[3 => :y] = y
        (trace, weight) = generate(at, (0.4, 3), constraints)
        choices = get_choices(trace)
        @test choices[3 => :y] == y
        @test get_retval(trace) == 0.4 + y
        @test isapprox(weight, logpdf(normal, y, 0.4, 1.))
        @test length(collect(get_values_shallow(choices))) == 0
        @test length(collect(get_submaps_shallow(choices))) == 1
    end

    function get_trace()
        y = 1.234
        constraints = choicemap()
        constraints[3 => :y] = y
        (trace, _) = generate(at, (0.4, 3), constraints)
        (trace, y)
    end

    @testset "project" begin
        (trace, y) = get_trace()
        @test isapprox(project(trace, EmptySelection()), 0.)
        selection = select(3 => :y)
        @test isapprox(project(trace, selection), logpdf(normal, y, 0.4, 1))
    end

    @testset "force update" begin
        (trace, y) = get_trace()

        # change kernel_args, same key, no constraint
        (new_trace, weight, retdiff, discard) = update(trace,
            (0.2, 3), (UnknownChange(), UnknownChange()), EmptyChoiceMap())
        choices = get_choices(new_trace)
        @test choices[3 => :y] == y
        @test length(collect(get_values_shallow(choices))) == 0
        @test length(collect(get_submaps_shallow(choices))) == 1
        @test isapprox(weight, logpdf(normal, y, 0.2, 1) - logpdf(normal, y, 0.4, 1))
        @test get_retval(new_trace) == 0.2 + y
        @test isempty(discard)
        @test isapprox(get_score(new_trace), logpdf(normal, y, 0.2, 1))

        # change kernel_args, same key, with constraint
        constraints = choicemap()
        y_new = 2.345
        constraints[3 => :y] = y_new
        (new_trace, weight, retdiff, discard) = update(trace,
            (0.2, 3), (UnknownChange(), UnknownChange()), constraints)
        choices = get_choices(new_trace)
        @test choices[3 => :y] == y_new
        @test length(collect(get_values_shallow(choices))) == 0
        @test length(collect(get_submaps_shallow(choices))) == 1
        @test isapprox(weight, logpdf(normal, y_new, 0.2, 1) - logpdf(normal, y, 0.4, 1))
        @test get_retval(new_trace) == 0.2 + y_new
        @test discard[3 => :y] == y
        @test length(collect(get_values_shallow(discard))) == 0
        @test length(collect(get_submaps_shallow(discard))) == 1
        @test isapprox(get_score(new_trace), logpdf(normal, y_new, 0.2, 1))

        # change kernel_args, different key, with constraint
        y_new = 2.345
        constraints = choicemap()
        constraints[4 => :y] = y_new
        (new_trace, weight, retdiff, discard) = update(trace,
            (0.2, 4), (UnknownChange(), UnknownChange()), constraints)
        choices = get_choices(new_trace)
        @test choices[4 => :y] == y_new
        @test length(collect(get_values_shallow(choices))) == 0
        @test length(collect(get_submaps_shallow(choices))) == 1
        @test isapprox(weight, logpdf(normal, y_new, 0.2, 1) - logpdf(normal, y, 0.4, 1))
        @test get_retval(new_trace) == 0.2 + y_new
        @test discard[3 => :y] == y
        @test length(collect(get_values_shallow(discard))) == 0
        @test length(collect(get_submaps_shallow(discard))) == 1
        @test isapprox(get_score(new_trace), logpdf(normal, y_new, 0.2, 1))
    end

    @testset "free update" begin
        (trace, y) = get_trace()

        # change kernel_args, same key, not selected
        (new_trace, weight, retdiff) = regenerate(trace,
            (0.2, 3), (UnknownChange(), UnknownChange()), EmptySelection())
        choices = get_choices(new_trace)
        @test choices[3 => :y] == y
        @test length(collect(get_values_shallow(choices))) == 0
        @test length(collect(get_submaps_shallow(choices))) == 1
        @test isapprox(weight, logpdf(normal, y, 0.2, 1) - logpdf(normal, y, 0.4, 1))
        @test get_retval(new_trace) == 0.2 + y
        @test isapprox(get_score(new_trace), logpdf(normal, y, 0.2, 1))

        # change kernel_args, same key, selected
        selection = select(3 => :y)
        (new_trace, weight, retdiff) = regenerate(trace,
            (0.2, 3), (UnknownChange(), UnknownChange()), selection)
        choices = get_choices(new_trace)
        y_new = choices[3 => :y]
        @test length(collect(get_values_shallow(choices))) == 0
        @test length(collect(get_submaps_shallow(choices))) == 1
        @test weight == 0.
        @test get_retval(new_trace) == 0.2 + y_new
        @test isapprox(get_score(new_trace), logpdf(normal, y_new, 0.2, 1))

        # change kernel_args, different key, not selected
        (new_trace, weight, retdiff) = regenerate(trace,
            (0.2, 4), (UnknownChange(), UnknownChange()), EmptySelection())
        choices = get_choices(new_trace)
        y_new = choices[4 => :y]
        @test length(collect(get_values_shallow(choices))) == 0
        @test length(collect(get_submaps_shallow(choices))) == 1
        @test weight == 0.
        @test get_retval(new_trace) == 0.2 + y_new
        @test isapprox(get_score(new_trace), logpdf(normal, y_new, 0.2, 1))
    end

    @testset "choice_gradients" begin
        (trace, y) = get_trace()

        # not selected
        retval_grad = 1.234
        (input_grads, choices, gradients) = choice_gradients(
            trace, EmptySelection(), retval_grad)
        @test isempty(choices)
        @test isempty(gradients)
        @test length(input_grads) == 2
        @test isapprox(input_grads[1], logpdf_grad(normal, y, 0.4, 1.0)[2] + retval_grad)
        @test input_grads[2] == nothing # the key has no gradient

        # selected with retval_grad
        retval_grad = 1.234
        selection = select(3 => :y)
        (input_grads, choices, gradients) = choice_gradients(
            trace, selection, retval_grad)
        @test choices[3 => :y] == y
        @test isapprox(gradients[3 => :y], logpdf_grad(normal, y, 0.4, 1.0)[1] + retval_grad)
        @test length(collect(get_values_shallow(gradients))) == 0
        @test length(collect(get_submaps_shallow(gradients))) == 1
        @test length(collect(get_values_shallow(choices))) == 0
        @test length(collect(get_submaps_shallow(choices))) == 1
        @test length(input_grads) == 2
        @test isapprox(input_grads[1], logpdf_grad(normal, y, 0.4, 1.0)[2] + retval_grad)
        @test input_grads[2] == nothing # the key has no gradient
    end

    @testset "accumulate_param_gradients!" begin
        (trace, y) = get_trace()
        retval_grad = 1.234
        input_grads = accumulate_param_gradients!(trace, retval_grad)
        @test length(input_grads) == 2
        @test isapprox(input_grads[1], logpdf_grad(normal, y, 0.4, 1)[2] + retval_grad)
        @test input_grads[2] == nothing # the key has no gradient
    end
end


text:
export bloobs, bars

function bloobs(;centers = 3, samples = 100, randomize = true)

  μ = reduce(hcat, [rand(Uniform(-10, 10), centers) for i in 1:2])
  Σ = reduce(hcat, [rand(Uniform(0.5, 4), centers) for i in 1:2])

  samplespcenter = ones(Int, centers) * round(Int, samples / centers)

  for i = 1:(centers % samples)
      samplespcenter[i] += 1
  end

  X = reduce(hcat, [rand(MvNormal(ones(2) .* μ[i], Matrix(I, 2, 2) .* Σ[i]), samplespcenter[i]) for i in 1:centers])'
  Y = reduce(vcat, [ones(Int, samplespcenter[i]) * i for i in 1:centers])

  ids = collect(1:size(X, 1))
  if randomize
    shuffle!(ids)
  end

  X = X[ids,:]
  Y = Y[ids]

  return (X, Y)
end

function bars(;img_size = 6, noise_level = 1e-4, num_per_mixture = [0 ones(1,3) / 3], num_group = 20, num_data = 50)

  numdim   = img_size^2;
  numbars  = img_size*2

  mix_theta = zeros(numdim, numbars)

  # add horizontal bars
  for i in 1:img_size
    img = ones(img_size, img_size) .* noise_level / img_size^2;
    img[:,i] = img[:,i] + 1/img_size;
    mix_theta[:,i] = img[:] / sum(img[:]);
  end

  # add vertical bars
  for i = 1:img_size
    img = ones(img_size, img_size) .* noise_level / img_size^2;
    img[i,:] = img[i,:] + 1/img_size;
    mix_theta[:,img_size+i] = img[:]/sum(img[:]);
  end

  # generate samples
  # Note: strange behavior of cumsum!!!
  cumbar = cumsum(num_per_mixture, 2)

  samples = Vector{Vector{Int}}(num_group)

  for g in 1:num_group

      # random number of mixture components
      n_bars = 1 + sum(rand() .> cumbar)
      k = randperm(numbars)[1:n_bars]

      # get weigths
      theta = mean(mix_theta[:,k], 2)

      # get samples (e.g. words)
      samples[g] = vec( 1 + sum(repmat(rand(num_data)', numdim, 1) .> repmat(cumsum(theta, 1), 1, num_data), 1) )
  end

  # return samples and bars
  return (samples, mix_theta)

end


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

export
    edgemincycles, edgemincyclemembership, mincycles, mincyclemembership,
    canonicalcycle

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

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

Calculate minimum cycle basis of the unweighted undirected graph by de Pina algorithm re-interpreted by Kavitha et al. This returns an array of edge sets that consist of the cycles.
"""
@cachefirst function edgemincycles(graph::UndirectedGraph)
    cycles = Vector{Int}[]
    for conn in connectedcomponents(graph)
        subg = nodesubgraph(graph, conn)
        S = [Set(e) for e in cotree_edges(subg)]
        N = length(S)  # N: circuit rank
        for k in 1:N
            minedges = findmincycle(subg, S[k])
            push!(cycles, minedges)
            for i in (k + 1):N
                if length(intersect(S[i], minedges)) % 2 == 1
                    S[i] = symdiff(S[i], S[k])
                end
            end
        end
    end
    return cycles
end

@cachefirst function edgemincyclemembership(graph::OrderedGraph)
    edges = [Set{Int}() for n in 1:edgecount(graph)]
    for (i, cyc) in enumerate(edgemincycles(graph))
        for e in cyc
            push!(edges[e], i)
        end
    end
    return edges
end

edgemincyclemembership(view::SubgraphView) = edgemincyclemembership(view.graph)

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

Return an array of node sets that consist of the minimum cycle basis.
"""
@cachefirst function mincycles(graph::UndirectedGraph)
    mincycs = Vector{Int}[]
    for cy in edgemincycles(graph)
        cycy = vcat(cy, cy)
        nodes = Int[]
        for i in 1:length(cy)
            a = getedge(graph, cycy[i])
            b = getedge(graph, cycy[i + 1])
            push!(nodes, iterate(intersect(a, b))[1])
        end
        push!(mincycs, nodes)
    end
    return mincycs
end

mincycles(view::SubgraphView) = mincycles(view.graph)

@cachefirst function mincyclemembership(graph::OrderedGraph)
    nodes = [Set{Int}() for n in 1:nodecount(graph)]
    for (i, cyc) in enumerate(mincycles(graph))
        for n in cyc
            push!(nodes[n], i)
        end
    end
    return nodes
end

mincyclemembership(view::SubgraphView) = mincyclemembership(view.graph)

function canonicalcycle(nodes::Vector{Int})
    """Align cycle indices to start from lowest index and following one of
    neighbors that have the lower index
    """
    # TODO: not used
    (fst, fstidx) = findmin(nodes)
    succidx = fstidx == lastindex(nodes) ? 1 : fstidx + 1
    succ = nodes[succidx]
    predidx = fstidx == 1 ? lastindex(nodes) : fstidx - 1
    pred = nodes[predidx]
    cp = succ < pred ? copy(nodes) : reverse(copy(nodes))
    while cp[1] != fst
        pushfirst!(cp, pop!(cp))
    end
    return cp
end