text:
export revert

revert(x::Tuple) = reverse(x)
revert(x::String) = Reverse(x)

export Reverse
@auto_hash_equals struct Reverse{V}<:AbstractString
    x::V
    lastindex::Int
    Reverse(x) =
        # if lastindex(x) == 1
        #     x
        # else
        new{typeof(x)}(x,lastindex(x))
        ##end
end
revert(x::Reverse) = x.x
reverse_index(x::Reverse,i) =
    x.lastindex-i+1
reverse_index(x::AbstractString,i) =
    i

Base.SubString(x::Reverse,start::Int,stop::Int) =
    SubString(x.x, reverse_index(x, stop), reverse_index(x, start))

function set_capture(sequence::Reverse, index::Int, x)
    @warn "check"
    set_capture(sequence.x,index,x)
end

regex_string(x::Reverse) = regex_escape(x.x)
Base.ncodeunits(x::Reverse) = ncodeunits(x.x)
Base.firstindex(x::Reverse) = 1
Base.lastindex(x::Reverse) = x.lastindex
Base.getindex(x::Reverse,is::UnitRange) = getindex(x.x,reverse_index(x,is.stop):reverse_index(x,is.start))
Base.getindex(x::Reverse,i::Int) = getindex(x.x,reverse_index(x,i))
Base.nextind(x::Reverse,i::Int) = reverse_index(x,prevind(x.x,reverse_index(x,i)))
Base.prevind(x::Reverse,i::Int) = reverse_index(x,nextind(x.x,reverse_index(x,i)))
Base.nextind(x::Reverse,i::Int,n::Int) = reverse_index(x,prevind(x.x,reverse_index(x,i),n))
Base.prevind(x::Reverse,i::Int,n::Int) = reverse_index(x,nextind(x.x,reverse_index(x,i),n))
Base.iterate(x::Reverse) =
    iterate(x,1)
Base.iterate(x::Reverse,i::Int) =
    x[i], nextind(x,i)

export PositiveLookbehind
"""
    PositiveLookbehind(parser)

Parser that succeeds if and only if `parser` succeeds **before cursor**. Consumes no input.
The match is returned.
Useful for checks like "must be preceded by `parser`, don't consume its match".
"""
@auto_hash_equals struct PositiveLookbehind{T,P} <: LookAround{T}
    parser::P
    function PositiveLookbehind(p_,revert_parser=true)
        p = revert_parser ? deepmap_parser(revert,IdDict(),parser(p_)) : parser(p_)
        new{result_type(p),typeof(p)}(p)
    end
end
# result_type(p::Type{PositiveLookbehind{T}}) where T = T
regex_prefix(x::PositiveLookbehind) = "(?<="

export NegativeLookbehind
"""
    NegativeLookbehind(parser)

Parser that succeeds if and only if `parser` does not succeed **before cursor**.  Consumes no input.
`nothing` is returned as match.
Useful for checks like "must not be preceded by `parser`, don't consume its match".

```jldoctest
julia> la=NegativeLookbehind("keep")
re"(?<!keep)"

julia> parse("peek"*la,"peek")
("peek", re"(?<!keep)")
```
"""
@auto_hash_equals struct NegativeLookbehind{P} <: LookAround{NegativeLookbehind{P}}
    parser::P
    function NegativeLookbehind(p_,revert_parser=true)
        p = revert_parser ? deepmap_parser(revert,IdDict(),parser(p_)) : parser(p_)
        new{typeof(p)}(p)
    end
end
regex_prefix(x::NegativeLookbehind) = "(?<!"

export Lookbehind
"""
    Lookbehind(does_match::Bool, p)

[`PositiveLookbehind`](@ref) if `does_match==true`,
[`NegativeLookbehind`](@ref) otherwise.
"""
function Lookbehind(does_match::Bool, p)
    if does_match
        PositiveLookbehind(p)
    else
        NegativeLookbehind(p)
    end
end
@deprecate look_behind(does_match,p) Lookbehind(does_match, p)

function _iterate(t::NegativeLookbehind, str, till, posi, next_i, state::Nothing)
    rseq=revert(str)
    next_i < 1 && return next_i, MatchState()
    r = _iterate(t.parser, rseq, till,
                 reverse_index(rseq,prevind(str,next_i)), nothing)
    if r === nothing
        next_i,MatchState()
    else
        nothing
    end
end

_iterate(t::PositiveLookbehind, str, till, posi, next_i, state::MatchState) =
    nothing

function _iterate(t::PositiveLookbehind, str, till, posi, next_i, state::Nothing)
    rseq=revert(str)
    next_i < 1 && return nothing
    r = _iterate(t.parser, rseq, till,
                 reverse_index(rseq,prevind(rseq,next_i)), nothing)
    if r === nothing
        nothing
    else
        next_i,MatchState()
    end
end

regex_inner(x::Union{PositiveLookbehind,NegativeLookbehind}) =
    regex_inner(revert(x.parser))

children(x::Union{PositiveLookbehind,NegativeLookbehind}) =
    reverse(children(x.parser))

for T in [PositiveLookahead,NegativeLookahead,PositiveLookbehind,NegativeLookbehind]
    eval(quote
         deepmap_parser(f::Function,mem::AbstractDict,x_::$T,a...; kw...) =
         let x = deepmap_parser(f,mem,x_.parser,a...; kw...)
         $T(x,false)
         end
         end)
end

revert(x::Union{AnyChar,CharIn,CharNotIn,UnicodeClass,Always,Never,ConstantParser{N,Char} where N}) = x
revert(x::AtStart) = AtEnd()
revert(x::AtEnd) = AtStart()
deepmap_parser(::typeof(revert),mem::AbstractDict,x::Sequence) =
    get!(mem,x) do
        Sequence(( deepmap_parser(revert,mem,p) for p in reverse(x.parts) )...)
    end
deepmap_parser(::typeof(revert),mem::AbstractDict,x::NegativeLookbehind) =
    get!(mem,x) do
        NegativeLookahead(x.parser) ##deepmap_parser(revert,@show x.parser))
    end
deepmap_parser(::typeof(revert),mem::AbstractDict,x::NegativeLookahead) =
    get!(mem,x) do
        NegativeLookbehind(x.parser) ##deepmap_parser(revert,x.parser))
    end
deepmap_parser(::typeof(revert),mem::AbstractDict,x::PositiveLookbehind) =
    get!(mem,x) do
        PositiveLookahead(x.parser) ##deepmap_parser(revert,x.parser))
    end
deepmap_parser(::typeof(revert),mem::AbstractDict,x::PositiveLookahead) =
    get!(mem,x) do
        PositiveLookbehind(x.parser) ##deepmap_parser(revert,x.parser))
    end


text:
################################################################################
#
#  Constructor
#
################################################################################

@doc Markdown.doc"""
    NumberField(f::Vector{PolyElem{<:NumFieldElem}}, s::String="_\$", check = true)
                                              -> NumField, Vector{NumFieldElem}

Given a list $f_1, \ldots, f_n$ of univariate polynomials in $K[x]$ over
some number field $K$, constructs the extension $K[x_1, \ldots, x_n]/(f_1(x_1),
\ldots, f_n(x_n))$.
"""
NumberField(::Vector{PolyElem{<:Union{NumFieldElem, fmpq}}}, ::String, check::Bool = true)

################################################################################
#
#  Basis
#
################################################################################

@doc Markdown.doc"""
    basis(L::NonSimpleNumField) -> Vector{NumFieldElem}

Returns the canonical basis of a non-simple extension $L/K$. If
$L = K(a_1,\dotsc,a_n)$ where each $a_i$ has degree $d_i$, then the basis will
be $a_1^{i_1}\dotsm a_d^{i_d}$ with $0 \leq i_j d_j - 1$ for $1 \leq j \leq n$.
"""
basis(::NonSimpleNumField)

################################################################################
#
#  Conversion to simple extension
#
################################################################################

@doc Markdown.doc"""
    simple_extension(L::NonSimpleNumField) -> SimpleNumField, Map

Given a non-simple extension $L/K$, this function computes a simple extension
$M/K$ and a $K$-linear isomorphism $M \to L$.
"""
simple_extension(::NonSimpleNumField)

################################################################################
#
#  Check consistency
#
################################################################################

# This is used for the check in the constructor
function _check_consistency(K::NonSimpleNumField)
  QQz, z = PolynomialRing(base_field(K), "z")
  for i = 1:length(K.pol)
    v = [zero(QQz) for i in 1:length(K.pol)]
    v[i] = z
    p = evaluate(K.pol[i], v)
    if !isirreducible(p)
      return false
    end
  end
  lg = gens(K)
  el = lg[1]
  f = minpoly(el)
  d = degree(f)
  for i = 2:length(lg)
    el += lg[i]
    f = minpoly(el)
    while !issquarefree(f)
      el += lg[i]
      f = minpoly(el)
    end
    if degree(f) != prod(degree(K.pol[j], j) for j = 1:i)
      return false
    end
    if !isirreducible(f)
      return false
    end
  end
  return true
end

ext: jl


text:
import .Luxor: Drawing, finish, Point, background, sethue, fontsize, rect, text

getname(p::Dagger.OSProc) = "OS Process on worker $(p.pid)"
getname(p::Dagger.ThreadProc) = "Thread $(p.tid) on worker $(p.owner)"
getname(p) = sprint(Base.show, p)

function proclt(p1::T, p2::R) where {T,R}
    if p1.owner != p2.owner
        return p1.owner < p2.owner
    else
        return repr(T) < repr(R)
    end
end
function proclt(p1::T, p2::T) where {T}
    if p1.owner != p2.owner
        return p1.owner < p2.owner
    else
        for field in fieldnames(T)
            f1 = getfield(p1, field)
            f2 = getfield(p2, field)
            if f1 != f2
                return f1 < f2
            end
        end
    end
    false
end
proclt(p1::Dagger.OSProc, p2::Dagger.OSProc) = p1.pid < p2.pid
proclt(p1::Dagger.OSProc, p2) = p1.pid < p2.owner
proclt(p1, p2::Dagger.OSProc) = p1.owner < p2.pid

prof_to_svg(::Any, ::Any, ::Any; kwargs...) = nothing

combine_gantt_images(::Any, ::Any, ::Any, ::Any) = ("", "")

continue_rendering = Ref{Bool}(true)
render_results = Channel()

function draw_gantt(ctx, svg_path, prof_path; delay=2, width=1000, height=640, window_length=20)
    root_time = time_ns()
    window_logs = []
    if window_length !== nothing
        window_start = -window_length
    else
        window_start = 0
    end
    window_finish = 0
    procs = Set()
    for gp in ctx.procs
        for proc in Dagger.get_processors(gp)
            push!(procs, proc)
        end
    end
    if (height/length(procs)) < 50
        height = length(procs) * 50
        @warn "SVG height too small; resizing to $height pixels"
    end

    image_idx = 0
    while continue_rendering[]
        sleep(delay)
        logs = Dagger.get_logs!(ctx.log_sink, true) # get raw events
        if !isempty(logs)
            for id in keys(logs)
                append!(window_logs, map(x->(x,), filter(x->x.category==:compute, logs[id])))
            end
        end
        for idx in length(window_logs):-1:1
            log = window_logs[idx]
            if length(log) == 2
                # Clear out finished events older than window start
                log_finish_s = (log[2].timestamp-root_time)/(1000^3)
                if log_finish_s < window_start
                    @debug "Gantt: Deleted event"
                    deleteat!(window_logs, idx)
                end
            elseif log[1] isa Dagger.Event{:finish}
                # Pair finish events with start events
                sidx = findfirst(x->length(x) == 1 &&
                                    x[1] isa Dagger.Event{:start} &&
                                    x[1].id==log[1].id, window_logs)
                if sidx === nothing
                    @debug "Gantt: Removed unpaired finish"
                    deleteat!(window_logs, idx)
                    continue
                end
                window_logs[sidx] = (window_logs[sidx][1], log[1])
                @debug "Gantt: Paired event"
                deleteat!(window_logs, idx)
            end
        end
        isempty(window_logs) && continue

        # Concatenate and render profile data
        if ctx.profile
            prof_data = UInt64[]
            prof_dict = Dict{UInt64, Vector{Base.StackTraces.StackFrame}}()
            for log in filter(x->length(x)==2, window_logs)
                append!(prof_data, log[2].profiler_samples.samples)
                merge!(prof_dict, log[2].profiler_samples.lineinfo)
            end
            prof_to_svg(prof_path, prof_data, prof_dict, image_idx; width=width)
        end

        for proc in unique(map(x->x[1].timeline[2], window_logs))
            push!(procs, proc)
        end
        colors = Colors.distinguishable_colors(length(procs))
        procs_len = length(procs)
        proc_height = height/(3procs_len)
        @debug "Gantt: Start"
        if isfile(svg_path)
            Drawing(width, height, svg_path)
        else
            Drawing(width, height, joinpath(svg_path, repr(image_idx) * ".svg"))
        end
        background("white")
        for (proc_idx, proc) in enumerate(sort(collect(procs); lt=proclt))
            ypos = (proc_idx-0.5)*(height/procs_len)
            sethue("grey15")
            fontsize(round(Int,proc_height))
            text(getname(proc), Point(width/2,ypos-(proc_height/2)))
            rect(Point(1,ypos+(proc_height/3)),width-2,proc_height,:stroke)
            fontsize(8)
            text("$(window_start) s", Point(0,ypos); halign=:left)
            text("$(window_finish) s", Point(width-8,ypos); halign=:right)
            proc_color = colors[proc_idx]
            for log in filter(x->x[1].timeline[2]==proc, window_logs)
                length(log) == 1 && log[1] isa Dagger.Event{:finish} && error("Unpaired finish!")
                log_start_s = (log[1].timestamp-root_time)/(1000^3)
                log_finish_s = if length(log) == 2
                    log_finish_s = (log[2].timestamp-root_time)/(1000^3)
                else
                    window_finish+1
                end
                xstart = ((log_start_s-window_start)/(window_finish-window_start))*width
                xfinish = ((log_finish_s-window_start)/(window_finish-window_start))*width
                sethue(proc_color)
                rect(Point(xstart,ypos+(proc_height/3)+1),xfinish-xstart,proc_height-2,:fill)
                sethue("black")
                rect(Point(xstart,ypos+(proc_height/3)+1),xfinish-xstart,proc_height-2,:stroke)
            end
        end
        finish()
        @debug "Gantt: Finish"
        window_finish = (time_ns()-root_time)/(1000^3)
        if window_length !== nothing
            window_start = window_finish-window_length
        end
        image_idx += 1
    end
    if isdir(svg_path)
        final_paths = try
            combine_gantt_images(ctx, svg_path, prof_path, delay)
        catch err
            @error "Image-to-video failed" exception=err
            ("", "")
        end
        continue_rendering[] = true
        put!(render_results, final_paths)
    end
end

function show_gantt(ctx; delay=2, port=8000, width=1000, height=640, window_length=20, live=false)
    svg_prefix = "/tmp/dagger-gantt-$(String(rand(UInt8(97):UInt8(122), 16)))"
    prof_prefix = "/tmp/dagger-gantt-$(String(rand(UInt8(97):UInt8(122), 16)))"
    if live
        # Create SVG file
        # TODO: Ensure this gets cleaned up
        svg_path = svg_prefix*".svg"
        prof_path = prof_prefix*".svg"
        open(svg_path, "w") do io
            write(io, "
Waiting for Gantt data...
")
        end
        open(prof_path, "w") do io
            write(io, "
Waiting for profile data...
")
        end
    else
        svg_path = svg_prefix
        prof_path = prof_prefix
        mkpath(svg_path)
        mkpath(prof_path)
    end

    # Start drawing task
    Threads.@spawn begin
        try
            draw_gantt(ctx, svg_path, prof_path; delay=delay, width=width,
                       height=height, window_length=window_length)
        catch err
            @error exception=(err,catch_backtrace())
        end
    end

    if live
        serve_gantt(svg_path, prof_path; port=port, delay=delay)
    end
end


text:
struct TwoInts
    dim1::Int
    dim2::Int
end

Base.length(::TwoInts) = 2
Base.iterate(ti::TwoInts) = ti.dim1, false
function Base.iterate(ti::TwoInts, isdone::Bool)
    isdone && return nothing
    return ti.dim2, true
end

Base.:(-)(ti::TwoInts) = TwoInts(-ti.dim1, -ti.dim2)
Base.:(+)(v::SVector{2}, ti::TwoInts) = TwoInts(v[1] + ti.dim1, v[2] + ti.dim2)


text:
addprocs(3) # Remove addprocs if you want to use a single process only
using IncrementalInference
using KernelDensityEstimate, Gadfly # for vstack
using Cairo, Fontconfig # for drawing PNG/PDF

gt = Dict{Symbol, Array{Float64,2}}()
# HMM computed ground truth for first 3 poses only
gt[:x1]=[[-100.0; 1.96]';[0.0; 1.96]']'
gt[:x2]=[[-50.0; 3.1]';[50.0; 3.1]']'
gt[:x3]=[[100.0; 3.05]';[0.0; 3.05]']'

fg = initfg()

N=100

doors = reshape([-100.0;0.0;100.0;300.0],1,4)
cov = 3.0*ones(1,1)
# pd = kde!(doors,cov)
# pd = resample(pd,N);
# bws = getBW(pd)[:,1]
# doors2 = getPoints(pd);

v1 = addVariable!(fg,:x1,ContinuousScalar,N=N)
f1  = addFactor!(fg,[v1], Obsv2(doors, cov, [1.0]))

tem = 2.0*randn(1,N)+getVal(v1)+50.0
v2 = addVariable!(fg,:x2, ContinuousScalar, N=N)
addFactor!(fg,[v1;v2],Odo(50.0*ones(1,1),2.0*ones(1,1),[1.0]))

# monocular sighting would look something like
#addFactor!(fg, Mono, [:x3,:l1], [14.0], [1.0], [1.0])
#addFactor!(fg, Mono, [:x4,:l1], [11.0], [1.0], [1.0])

v3=addVariable!(fg,:x3,ContinuousScalar, N=N)
addFactor!(fg,[v2;v3],Odo(50.0*ones(1,1),4.0*ones(1,1),[1.0]))
f2 = addFactor!(fg,[v3], Obsv2(doors, cov', [1.0]))

if true

    v4=addVariable!(fg,:x4,ContinuousScalar, N=N)
    addFactor!(fg,[v3;v4],Odo(50.0*ones(1,1),2.0*ones(1,1),[1.0]))

    if true
        l1=addVariable!(fg, :l1, ContinuousScalar, N=N)
        addFactor!(fg, [v3,l1], Ranged([64.0],[0.5],[1.0]))
        addFactor!(fg, [v4,l1], Ranged([16.0],[0.5],[1.0]))
    end

    v5=addVariable!(fg,:x5,ContinuousScalar, N=N)
    addFactor!(fg,[v4;v5],Odo(50.0*ones(1,1),2.0*ones(1,1),[1.0]))

    v6=addVariable!(fg,:x6,ContinuousScalar, N=N)
    addFactor!(fg,[v5;v6],Odo(40.0*ones(1,1),1.25*ones(1,1),[1.0]))

    v7=addVariable!(fg,:x7,ContinuousScalar, N=N)
    addFactor!(fg,[v6;v7],Odo(60.0*ones(1,1),2.0*ones(1,1),[1.0]))

    f3 = addFactor!(fg,[v7], Obsv2(doors, cov, [1.0]))

    # HMM computed ground truth, extended for 7 poses with landmark
    gt = Dict{String, Array{Float64,2}}()
    gt[:x1]=reshape(Float64[0.0;1.97304 ],2,1) # -0.0342366
    gt[:x2]=reshape(Float64[50.0; 2.83153 ],2,1) # 49.8797
    gt[:x3]=reshape(Float64[100.0; 1.65557 ],2,1) # 99.8351
    gt[:x4]=reshape(Float64[150.0; 1.64945 ],2,1) # 148.637
    gt[:x5]=reshape(Float64[200.0; 1.77992 ],2,1) # 198.62
    gt[:x6]=reshape(Float64[240.0; 2.20466 ],2,1) # 238.492
    gt[:x7]=reshape(Float64[300.0; 2.14353 ],2,1) # 298.467
    gt[:l1]=reshape(Float64[165.0; 1.17284 ],2,1) # 164.102
end

# writeGraphPdf(fg);
tree = prepBatchTree!(fg,drawpdf=true);

# spyCliqMat(tree.cliques[1])

# list vertices in fg
@show xx,ll = ls(fg)

# from = fg.v[fg.IDs[:l1]]
# to = fg.v[fg.IDs[:x7]]
# fgs = subgraphShortestPath(fg, from=from, to=to)

# using Graphs
# el = shortest_path(fg.g, ones(19),from, to)

# do belief propagation inference over tree once
[inferOverTree!(fg, tree, N=100) for i in 1:1];

# draw all beliefs
# TODO -- Update this part of the example
using RoMEPlotting

DOYTICKS = false
xx,ll = ls(fg)
msgPlots = drawHorBeliefsList(fg, xx, gt=gt,nhor=2);
evalstr = ""
for i in 1:length(msgPlots)
    evalstr = string(evalstr, ",msgPlots[$(i)]")
end
pl = eval(parse(string("vstack(",evalstr[2:end],")")));
Gadfly.draw(PNG("4doors.png",15cm,20cm),pl) # can also do PNG

# if false
#
#   # draw upward messages
#   msgPlots=drawTreeUpwardMsgs(fg, tree, N=500); # to init memory for eval(parse(string))
#   println("Upward messages for all cliques except root")
#   # vvMsgs = vstackedDensities(fg, tree, msgPlots)
#   Gadfly.set_default_plot_size(17cm, 15cm)
#   Gadfly.draw(PNG("results/testMsgs.png",17cm,15cm),vvMsgs)
#   # vvMsgs
# end


text:
# Provide the interface as a module so that different sub-modules can reuse the symbols for dispatching
module SixelInterface

export AbstractSixelDecoder,
       AbstractSixelEncoder,
       canonical_sixel_eltype

abstract type AbstractSixelEncoder end
(enc::AbstractSixelEncoder)(io, src) = error("The encoder functor method for inputs (::$(typeof(io)), ::$(typeof(src)) is not implemented.")

abstract type AbstractSixelDecoder end
(enc::AbstractSixelDecoder)(io, src) = error("The decoder functor method for inputs (::$(typeof(io)), ::$(typeof(src)) is not implemented.")

"""
    canonical_sixel_eltype(enc, CT1) -> CT2

Given input type `CT1`, infer the expected colorant type `CT2` for encoder `enc`. This function is
used to early convert the image data into formats that backend encoder (e.g., the C libsixel) understands.
"""
function canonical_sixel_eltype(::AbstractSixelEncoder, ::DataType) end

end # module


text:
# The code to generate S4 Fig and Data for S1 Table
using PyPlot
using SurNoR_2020
using Statistics
using CSV
using DataFrames

PyPlot.svg(true)
rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
rcParams["svg.fonttype"] = "none"
rcParams["pdf.fonttype"] = 42

# A function to read the parameters' name and their optimization range given
# their index
function Func_if_VarName(num_param)
    if num_param==1
        # P_change              # The name in the codes
        VarName = "m/(1+m)"     # The name in the paper
        x = 0.00:0.0025:0.3
    elseif num_param==2
        # Eps
        VarName = "eps"
        x = -12:0.2:0.2
        x = 10 .^ x
    elseif num_param==3
        # β_R
        VarName = "beta_1"
        x = 0.1:0.2:15
    elseif num_param==4
        # β_RN1
        VarName = "beta_N1"
        x = 0.00:0.005:0.6
    elseif num_param==5
        # β_RN2
        VarName = "beta_N2"
        x = 0.00:0.005:0.6
    elseif num_param==6
        # γ_R
        VarName = "λ_R"
        x = 0.85:0.01:0.99
    elseif num_param==7
        # γ_N
        VarName = "λ_N"
        x = 0.20:0.05:0.99
    elseif num_param==8
        # λ_R
        VarName = "μ_R"
        x = 0.85:0.01:1
    elseif num_param==9
        # λ_N
        VarName = "μ_N"
        x = 0.20:0.05:1
    elseif num_param==10
        # α
        VarName = "ρ_b"
        x = 0.0:0.02:0.2
    elseif num_param==11
        # Δα
        VarName = "delta ρ"
        x = 0.0:0.05:(1-Param1.α)
    elseif num_param==12
        # Q_N_0
        VarName = "Q_N0"
        x = 0.0:0.25:12
    elseif num_param==13
        # N_PS
        VarName = "T_PS"
        x = cat(0:14,15:10:60,dims=1)
    elseif num_param==14
        # w_2
        VarName = "w_0"
        x = 0:0.05:1
    elseif num_param==15
        # w_11
        VarName = "w_11"
        x = 0:0.05:1
    elseif num_param==16
        # w_12
        VarName = "w_12"
        x = 0:0.05:1
    elseif num_param==17
        # w_scaler
        VarName = "w_scale"
        x = 0:0.1:6
    elseif num_param==18
        # β_R_2
        VarName = "beta_2"
        x = 0.1:0.2:15
    end
    return (VarName,x)
end

# 3 sets of recovered parameters
N_set = [1,2,3]

# Loading fitted parameters
Path_Params = "data/Fitted_parameters/Overall/"
Path_Params = string(Path_Params, "OverAll_SurNoR_BestParams_diff5.csv" )
Param_Data = CSV.read(Path_Params)
Row = 3

# Loading recovered parameters
Param_Data_nset_all = []
for n_set = N_set
    Path_Params_nset = "data/Fitted_parameters_simulated/Overall/"
    Path_Params_nset = string(Path_Params_nset, "nset_", string(n_set), "/", "OverAll_SurNoR_BestParams_diff5.csv" )
    Param_Data_nset = CSV.read(Path_Params_nset)
    push!(Param_Data_nset_all,Param_Data_nset)
end
Row_nset = [1,1,1]

# Parameter structure
Param1 = SurNoR_2020.Str_Param(; num_act=4, num_obs=11,
                                 count_decay=1, Q_R_0 = 0)

Param1.p_change = Param_Data[Row,1+1]
Param1.ϵ = Param_Data[Row,1+2]
Param1.β_R = Param_Data[Row,1+3]
Param1.β_ratio_N2R = Param_Data[Row,1+4]
Param1.γ_RL_R = Param_Data[Row,1+6]
Param1.γ_RL_N = Param_Data[Row,1+7]
Param1.λ_elig_R = Param_Data[Row,1+8]
Param1.λ_elig_N = Param_Data[Row,1+9]
Param1.α = Param_Data[Row,1+10]
Param1.Δα = Param_Data[Row,1+11]
Param1.Q_N_0 = Param_Data[Row,1+12]
Param1.N_prior_sweep = Param_Data[Row,1+13]
Param1.w_MF_to_MB_2 = Param_Data[Row,1+14]
Param1.w_MF_to_MB_1 = Param_Data[Row,1+15]
Param1.w_MF_to_MB_Scaler = Param_Data[Row,1+17]

Param2 = deepcopy(Param1)
Param2.β_ratio_N2R = Param_Data[Row,1+5]
Param2.w_MF_to_MB_1 = Param_Data[Row,1+16]
Param2.β_R = Param_Data[Row,1+18]

Param = [Param1, Param2]

# An array for rrrors in S1 Table
Sigma = zeros(size(Param_Data[Row,:]))

# Plotting panels of S4 Fig
Param_perm_order = 1:18
for i = 1:length(Param_perm_order)
    Path_Load = "data/LogL_Landscape/"
    num_param_i = Param_perm_order[i]
    (VarName_i, x_i) = Func_if_VarName(num_param_i)

    Path = string(Path_Load, "LogL_range_i", string(num_param_i), "_Row", string(Row), ".CSV" )

    @show VarName_i

    LogL_range = CSV.read(Path)
    LogL_range = Matrix(LogL_range)[:,2]

    LogL_range[isnan.(LogL_range)] .= -1e6
    LogL_range[isinf.(LogL_range)] .= -1e6

    y = LogL_range
    y_0 = findmax(y)[1]
    x_0_ind = findmax(y)[2]
    x_0 = x_i[x_0_ind]

    dx_res = [0.]
    # Laplace approximation
    if x_0_ind == 1
        dy = abs(y[x_0_ind+1] - y[x_0_ind])
        dx = abs(x_i[x_0_ind+1] - x_i[x_0_ind])
        c = 2 * dy / (dx^2)
        dx_res[1] = dx
    elseif x_0_ind == length(x_i)
        dy = abs(y[x_0_ind] - y[x_0_ind-1])
        dx = abs(x_i[x_0_ind] - x_i[x_0_ind-1])
        c = 2 * dy / (dx^2)
        dx_res[1] = dx
    else
        dy1 = abs(y[x_0_ind] - y[x_0_ind-1])
        dy2 = abs(y[x_0_ind+1] - y[x_0_ind])
        dx1 = abs(x_i[x_0_ind] - x_i[x_0_ind-1])
        dx2 = abs(x_i[x_0_ind+1] - x_i[x_0_ind])
        c = (dy1/dx1 + dy2/dx2) / (dx1 + dx2) * 2
        dx_res[1] = max(dx1,dx2)
    end
    Sigma[i+1] = max(1 / sqrt(c), dx_res[1])

    X_0_nset = []
    for n_set = N_set
        Param_Data_nset = Param_Data_nset_all[n_set]
        x_0_nset = Param_Data_nset[Row_nset[n_set],1+i]
        append!(X_0_nset,x_0_nset)
    end

    y = LogL_range
    y_max = findmax(y)[1]
    y_min = findmin(y)[1]

    fig = figure(); ax = gca()
    ax.plot(x_i,y,color="k")
    ax.plot([x_0,x_0],[y_min,y_max],color="#487952",alpha = 0.5)
    for n_set = N_set
        x_0_nset = X_0_nset[n_set]
        ax.plot([x_0_nset,x_0_nset],[y_min,y_max],color="#7B082E",alpha = 0.5)
    end
    title(string("SurNoR - Loss Shape - ", VarName_i))
    if VarName_i == "eps"
            xscale("log")
    end
end

# Saving data for S1 Table
A = push!(DataFrame(Param_Data[Row,:]), round.(Sigma, digits = 5))
Path_Save = "data/LogL_Landscape/"
CSV.write(string(Path_Save, "SurNoR_Params_Sigma_rounded.CSV" ),A)


text:
function parseMask(line)
    result = "XXXXXXXXXXXXXXXXXXXXXXXXXXXX" * String(split(line,"mask = ")[2])
    return result
end

struct Instruction
    pos::Int
    value::Int
end

function parseInstruction(line)
    insideSquareBrackets = r"\[([^()]*)\]"
    pos = match(insideSquareBrackets, line).captures[1]
    value = split(line, " = ")[2]
    result = Instruction(parse(Int, pos), parse(Int, value))
    return result
end

function applyCondition(x::Tuple{Char,Char})
    m = x[1]
    c = x[2]
    if m == 'X'
        return c
    elseif m == '0'
        return '0'
    elseif m == '1'
        return '1'
    else
        throw("wrong character for mask")
    end
end

function applyMask(mask::String, value::Int)
    valueBstr = bitstring(value)
    c = zip(mask, valueBstr)
    resultStr = String(map(x -> applyCondition(x), c))
    result = parse(Int, resultStr, base=2)
    return result
end

function execute(mask, instruction, memory)
    newValue = applyMask(mask, instruction.value)
    memory[instruction.pos] = newValue
end

function runProgram(inputFile)
    memory = zeros(Int, 100000)
    mask = ""
    for line in readlines(inputFile)
        if startswith(line,"mask")
            mask = parseMask(line)
        else
            instruction = parseInstruction(line)
            execute(mask, instruction, memory)
        end
    end
    println(reduce(+,memory))
end

runProgram("data.txt")


text:
function mass(u, q, w, v, x, t, mu, eta)
    m = [1.0, 1.0, 1.0];
    return m;
end
function flux(u, q, w, v, x, t, mu, eta)
    g = mu[1];
    r = u[1];  # total height: h
    ru = u[2]; # h*u
    rv = u[3]; # h*v
    r1 = 1/r;
    uv = ru*r1; # u
    vv = rv*r1; #v
    p = (0.5*g)*(r*r); # 0.5*g*h^2
    f = [ru, ru*uv+p, rv*uv, rv, ru*vv, rv*vv+p];
    # [h*u h*u^2 + 0.5*g*h^2  huv h*v h*u*v h*v^2+0.5*g*h^2]
    f = reshape(f,(3,2));

    return f;
end
function source(u, q, w, v, x, t, mu, eta)
    s = [0.0, 0.0, 0.0];
    return s;
end
function ubou(u, q, w, v, x, t, mu, eta, uhat, n, tau)
    ub = [0.0, 0.0, 0.0];
    return ub;
end
function fbou(u, q, w, v, x, t, mu, eta, uhat, n, tau)
    f = flux(u, q, w, v, x, t, mu, eta);
    fb = f[:,1]*n[1] + f[:,2]*n[2] + tau[1]*(u-uhat);
    return fb;
end
function initu(x, mu, eta)
    epsil = 0.1; # perturbation magnitude
    l = 0.5;     # Gaussian width
    k = 0.5;     # Sinusoidal wavenumber

    x1 = x[1];
    x2 = x[2];

    # The Bickley jet
    U = (1/cosh(x2))*(1/cosh(x2));

    # Slightly off-center vortical perturbations
    Psiprime = exp(-(x2 + l/10)*(x2 + l/10) / (2*(l*l))) * cos(k * x1) * cos(k * x2);

    # Vortical velocity fields (ũ, ṽ) = (-∂_y, +∂_x) ψ̃
    uprime =  Psiprime * (k * tan(k * x2) + x2 /(l*l));
    vprime = -Psiprime * k * tan(k * x1);

    u01 = 1.0; # h
    u02 = U + epsil * uprime; # h*u
    u03 = epsil * vprime;  # h*v

    u0 = [u01, u02, u03];
    return u0;
end

max_stars_repo_name: rloekvh/Exasim

max_issues_repo_name: rloekvh/Exasim

max_forks_repo_name: rloekvh/Exasim


text:
# Function to take input and return dict
function parse_bags(instructions)

    big_dict = Dict()

    # Loop over all instructions
    for instruction in instructions
        origin, contents = split(instruction, " contain ")

        contents = replace(contents, "." => "")
        contents = replace(contents, "bags" => "bag")
        origin = replace(origin, "bags" => "bag")

        # Split contents
        if ',' in contents
            big_dict[origin] = split(contents, ", ")
        else
            big_dict[origin] = [contents]
        end
    end
    return strip_number(big_dict)
end

# Function to strip the number off of text
# TODO parse number and make it a property
function strip_number(big_dict)
    for (bag, contents) in big_dict
        # println(bag)

        new_contents = []
        content_dict = Dict() # dict with "number" => nbags, "contents" => bag contents
        nbags = Dict() # nbags is dict with "bag_name" => number

        # Strip number from bag name and add to nbags
        for c in contents
            new_contents = vcat(new_contents, c[3:length(c)])
            nbags[c[3:length(c)]] = 0
            if c != "no other bag"
                nbags[c[3:length(c)]] += parse(Int, c[1])
            end
        end

        content_dict["number"] = nbags
        content_dict["contents"] = new_contents
        big_dict[bag] = content_dict
    end

    # println(big_dict)
    return big_dict
end

# Recursive function to follow bag trees
function bag_in_bag(big_dict, bag, bag_goal)

    # Loop over contents of current bag
    for c in big_dict[bag]["contents"]
        if c == " other bag"
            return false
        end

        if c == bag_goal
            return true
        elseif bag_in_bag(big_dict, c, bag_goal)
            return true
        end
    end
    return false
end

# Find if the goal is in the bag
function goal_in_bag(instructions, bag_goal)
    ins = parse_bags(instructions)

    n_goal = 0
    for (bag,_) in ins
        if bag_in_bag(ins, bag, bag_goal)
            n_goal += 1
        end
    end

    return n_goal
end

function num_bags(instructions, start_bag; orig=false)
    ins = parse_bags(instructions)
    sbag = ins[start_bag]

    nbags = 1

    # println("Checking ", sbag)
    for c in sbag["contents"]
        if c != " other bag"
            nbags += sbag["number"][c] * num_bags(instructions, c)
        end
    end

    # Return nbags-1 if this is the "original bag"
    if orig
        return nbags -1
    else
        return nbags
    end
end

instructions = readlines("./instructions.txt")

println("Part One: ", goal_in_bag(instructions, "shiny gold bag"))
println("Part Two: ", num_bags(instructions, "shiny gold bag", orig=true))


text:
"""
    dimension(::DiffusionProcess{T,DP,DW})

Return dimension of the stochastic process and driving Brownian motion.
"""
dimension(::DiffusionProcess{T,DP,DW}) where {T,DP,DW} = (
    process = DP,
    wiener = DW
)

dim_process(::DiffusionProcess{T,DP}) where {T,DP} = DP
dim_wiener(::DiffusionProcess{T,DP,DW}) where {T,DP,DW} = DW

"""
    Base.eltype(::DiffusionProcess{T}) where T = T

Return the datatype that each coordinate of the stochastic process is stored in.
"""
Base.eltype(::DiffusionProcess{T}) where T = T

"""
    state_space(::DiffusionProcess{T,DP,DW,SS})

Return the state space restrictions.
"""
state_space(::DiffusionProcess{T,DP,DW,SS}) where {T,DP,DW,SS} = SS

"""
    default_type(::DiffusionProcess{T,DP})

Allows for inference of data type that encodes the state space of a given
diffusion.
"""
default_type(::DiffusionProcess{T,DP}) where {T,DP} = SVector{DP,T}

"""
    default_wiener_type(::DiffusionProcess{T,DP,DW})

Allows for inference of data type that encodes the state space of the Brownian
motion driving a given diffusion process.
"""
default_wiener_type(::DiffusionProcess{T,DP,DW}) where {T,DP,DW} = SVector{DW,T}

"""
    parameters(P::DiffusionProcess)

Return a tuple of pairs of `parameter_name` => `parameter_value`.
"""
parameters(P::DiffusionProcess) = Dict(
    map(
        p->( p=>getfield(P, p) ),
        parameter_names(P)
    )
)

"""
    const_parameters(P::DiffusionProcess)

Return a tuple of pairs of `parameter_name` => `parameter_value`. Return only
those parameteres that are considered to be `constant`.
"""
const_parameters(P::DiffusionProcess) = Dict(
    map(
        p->( p=>getfield(P, p) ),
        const_parameter_names(P)
    )
)

"""
    var_parameters(P::DiffusionProcess)

Return a tuple of pairs of `parameter_name` => `parameter_value`. Return only
those parameteres that are considered to be `variable`.
"""
var_parameters(P::DiffusionProcess) = Dict(
    map(
        p->( p=>getfield(P, p) ),
        var_parameter_names(P)
    )
)

"""
    parameter_names(P::DiffusionProcess)

Return a tuple with the names of all paremeters.
"""
parameter_names(P::DiffusionProcess) = parameter_names(typeof(P))

"""
    parameter_names(::Type{<:DiffusionProcess})

Return a tuple with the names of all paremeters.
"""
function parameter_names(::Type{<:DiffusionProcess}) end

"""
    const_parameter_names(P::DiffusionProcess)

Return a tuple with the names of all paremeters that are considered to be
`constant`.
"""
const_parameter_names(P::DiffusionProcess) = const_parameter_names(typeof(P))

"""
    const_parameter_names(P::Type{<:DiffusionProcess})

Return a tuple with the names of all paremeters that are considered to be
`constant`.
"""
function const_parameter_names(P::Type{<:DiffusionProcess}) end

"""
    var_parameter_names(P::DiffusionProcess)

Return a tuple with the names of all paremeters that are considered to be
`variable`.
"""
var_parameter_names(P::DiffusionProcess) = var_parameter_names(typeof(P))

"""
    var_parameter_names(P::Type{<:DiffusionProcess})

Return a tuple with the names of all paremeters that are considered to be
`variable`.
"""
function var_parameter_names(P::Type{<:DiffusionProcess})
    const_pn = const_parameter_names(P)
    Tuple(filter(p->!(p in const_pn), parameter_names(P)))
end

a(t, x, P::DiffusionProcess) = σ(t, x, P) * σ(t, x, P)'

#------------------------------------------------------------------------------#
#                     Default fallbacks for linear diffusions
#------------------------------------------------------------------------------#
function b!(buffer, t, x, P::LinearDiffusion)
    B!(buffer, t, P)
    β!(buffer, t, P)
    mul!(buffer.b, buffer.B, x, true, true)
end

σ!(buffer, t, x, P::LinearDiffusion) = σ!(buffer, t, P)
a!(buffer, t, x, P::LinearDiffusion) = a!(buffer, t, P)

b(t, x, P::LinearDiffusion) = B(t, P)*x + β(t, P)
σ(t, x, P::LinearDiffusion) = σ(t, P)
a(t, x, P::LinearDiffusion) = a(t, P)
a(t, P::LinearDiffusion) = σ(t, P) * σ(t, P)'

#------------------------------------------------------------------------------#
#
#                   Default fallbacks for indexed time
#
#------------------------------------------------------------------------------#
_b((t,i)::IndexedTime, x, P::DiffusionProcess) = b(t, x, P)
_σ((t,i)::IndexedTime, x, P::DiffusionProcess) = σ(t, x, P)
_b!(buffer, (t,i)::IndexedTime, x, P::DiffusionProcess) = b!(buffer, t, x, P)
_σ!(buffer, (t,i)::IndexedTime, x, P::DiffusionProcess) = σ!(buffer, t, x, P)

"""
    nonhypo_σ((t,i)::IndexedTime, x, P::DiffusionProcess)

Return a sub-matrix of the full volatility matrix `σ` that consists of non-zero
rows of `σ`.
"""
nonhypo_σ((t,i)::IndexedTime, x, P::DiffusionProcess) = nonhypo_σ(t, x, P)

"""
    diagonaldiff(P::DiffusionProcess)

Indicator for whether the volatility coefficient is represented by a diagonal
matrix
"""
diagonaldiff(P::DiffusionProcess) = false

"""
    sparsediff(P::DiffusionProcess)

Indicator for whether the volatility coefficient is represented by a sparse
matrix
"""
sparsediff(P::DiffusionProcess) = false

"""
    diagonalBmat(P::DiffusionProcess)

Indicator for whether the B matrix (if exists) is represented by a diagonal
matrix
"""
diagonalBmat(P::DiffusionProcess) = false

"""
    sparseBmat(P::DiffusionProcess)

Indicator for whether the B matrix (if exists) is represented by a sparse
matrix
"""
sparseBmat(P::DiffusionProcess) = false

Base.zero(P::DiffusionProcess) = zero(P, Val(:process))

"""
    Base.zero(P::DiffusionProcess)

Instantiate a zero element that can represent a state of a diffusion.
"""
function Base.zero(P::DiffusionProcess, ::Val{:wiener})
    deft = default_wiener_type(P)
    zero(deft, dim_wiener(P), ismutable(deft))
end

"""
    Base.zero(P::DiffusionProcess)

Instantiate a zero element that can represent a state of a Brownian motion
driving a diffusion process.
"""
function Base.zero(P::DiffusionProcess, ::Val{:process})
    deft = default_type(P)
    zero(deft, dim_process(P), ismutable(deft))
end

"""
    Base.zero(K::Type, D, ::Val)

If `K` is a mutable type, then create `zeros` of dimension `D` and entries with
types `eltype(K)`. Otherwise, calls regular zero(`K`).
"""
Base.zero

Base.zero(K::Type, D, ::Val{true}) = zeros(eltype(K), D)
Base.zero(K::Type, D, ::Val{false}) = zero(K)

"""
    end_point_info(P::DiffusionProcess)

Return information about the end-point (works only if some information of this
kind has been passed at the time of defining a struct) TODO improve
"""
end_point_info(P::DiffusionProcess) = Dict(
    map(
        p->( p=>getfield(P, p) ),
        end_point_info_names(P)
    )
)

"""
    end_point_info_names(P::DiffusionProcess)

Return names of information pieces about the end-point (works only if some
information of this kind has been passed at the time of defining a struct)
TODO improve
"""
function end_point_info_names(P::Type{<:DiffusionProcess}) end

end_point_info_names(P::DiffusionProcess) = end_point_info_names(typeof(P))

#custom_zero(D::Integer, ::Type{K}) where K <: Array = zeros(eltype(K), D)
#custom_zero(D::Integer, ::Type{K}) where K = zero(K)

#=
function sequential_add!(save_to, x, y, a=true, b=true, c=false)
    for i in eachindex(save_to, x, y)
        @inbounds save_to[i] = a*x[i]+b*y[i]+c
    end
    return nothing
end

function sequential_mul!(save_to, x, y, a=true, b=true)
    for i in eachindex(save_to, x, y)
        @inbounds save_to[i] = a*x[i]*y[i]+b
    end
    return nothing
end

function sequential_mul_add!(save_to, X, Y, Z=Y, a=true, b=false, c=false)
    for i in eachindex(save_to, X, Y, Z)
        @inbounds save_to[i] = a * X[i] * Y[i] + b * Z[i] + c
    end
    return nothing
end
=#


text:
function bulk_mark(xgrid::ExtendableGrid{Tv,Ti}, refinement_indicators::Vector{T}, theta = 0.5; indicator_AT = ON_CELLS) where {T,Tv,Ti}

    xFaceCells = xgrid[FaceCells]
    xFaceNodes = xgrid[FaceNodes]
    nfaces = num_sources(xFaceNodes)

    # redistribute refinement_indicators to face indicators (if necessary)
    if indicator_AT == ON_CELLS
        refind4face::Array{T,1} = zeros(T,nfaces)
        cell::Int = 0
        for face = 1 : nfaces
            for k = 1 : 2
                cell = xFaceCells[k,face]
                if cell > 0
                    refind4face[face] += refinement_indicators[cell]
                end
            end
            if cell > 0
                refind4face[face] /= 2
            end
        end
    elseif indicator_AT == ON_FACES
        refind4face = refinement_indicators
    end

    # mark largest indicators until theta*total is reached
    p = Base.sortperm(refind4face[:], rev = true)
    totalsum = sum(refind4face)
    csum = 0
    j = 0
    facemarker = zeros(Bool,nfaces)
    while csum <= theta*totalsum
        j += 1
        csum += refind4face[p[j]]
        facemarker[p[j]] = true
    end

    return facemarker
end

# adaptive refinement rules fo TRIANGLE2D
# first k nodes are the CellNodes
# next m nodes are the CellFaces midpoints
# next node is the CellMidpoint (if needed)
const _no_refinement_Triangle2D = reshape([1 2 3],3,1)
const _red_refinement_Triangle2D = [1 4 6; 4 2 5; 6 5 3; 5 6 4]'
const _blueR_refinement_Triangle2D = [3 1 4; 4 2 5; 3 4 5]'
const _blueL_refinement_Triangle2D = [1 4 6; 2 3 4; 4 3 6]'
const _gree_refinement_Triangle2D = [3 1 4; 2 3 4]'

function RGB_refinement_rule(::Type{<:Triangle2D}, marked_faces::Array{Bool,1})
    if any(marked_faces) == true
        @assert marked_faces[1] = true # closure should always mark reference face
        if marked_faces[2] == true && marked_faces[3] == true
            return _red_refinement_Triangle2D
        elseif marked_faces[2] == true && marked_faces[3] == false
            return _blueR_refinement_Triangle2D
        elseif marked_faces[2] == false && marked_faces[3] == true
            return _blueL_refinement_Triangle2D
        elseif marked_faces[2] == false && marked_faces[3] == false
            return _gree_refinement_Triangle2D
        end
    else
        return _no_refinement_Triangle2D
    end
end

"""
$(TYPEDSIGNATURES)

generates a new ExtendableGrid by red-green-blue mesh refinement of triangular meshes, see e.g.

Carstensen, C.
--An Adaptive Mesh-Refining Algorithm Allowing for an H^1 Stable L^2 Projection onto Courant Finite Element Spaces--
Constr Approx 20, 549–564 (2004). https://doi.org/10.1007/s00365-003-0550-5

The bool array facemarkers determines which faces should be bisected. Note, that a closuring is performed
such that the first face in every triangle with a marked face is also refined.
"""
function RGB_refine(source_grid::ExtendableGrid{T,K}, facemarkers::Array{Bool,1}) where {T,K}

    xgrid = ExtendableGrid{T,K}()
    xgrid[CoordinateSystem]=source_grid[CoordinateSystem]

    # unpack stuff from source grid
    oldCoordinates = source_grid[Coordinates]
    oldCellGeometries = source_grid[CellGeometries]
    oldCellRegions = source_grid[CellRegions]
    EG = Base.unique(oldCellGeometries)

    # get dimension of CellGeometries
    # currently it is assumed to be the same for all cells
    dim = dim_element(EG[1])

    xCellNodes = VariableTargetAdjacency(K)
    xCellGeometries = []
    xCellRegions = zeros(K,0)
    oldCellNodes = source_grid[CellNodes]
    oldCellFaces = source_grid[CellFaces]
    oldFaceNodes = source_grid[FaceNodes]
    nfaces = num_sources(oldFaceNodes)
    ncells = num_sources(oldCellNodes)

    # closuring
    # @logmsg MoreInfo "RGB refinement with $(sum(facemarkers)) marked faces"
    is_refined = true
    closure_finished = false
    while closure_finished == false
        closure_finished = true
        for cell = 1 : ncells
            is_refined = false
            for f = 1 : 3
                if facemarkers[oldCellFaces[f,cell]] == true
                    is_refined = true
                    break;
                end
            end
            if is_refined && facemarkers[oldCellFaces[1,cell]] == false
                # mark also reference edge
                facemarkers[oldCellFaces[1,cell]] = true
                closure_finished = false
            end
        end
    end
    # @logmsg DeepInfo "marked faces after closure = $(sum(facemarkers))"

    # determine number of new vertices
    oldvertices = size(oldCoordinates,2)
    newvertices = 0
    for face = 1 : nfaces
        if facemarkers[face] == true
            newvertices += 1
        end
    end
    xCoordinates = zeros(T,size(oldCoordinates,1),oldvertices+newvertices)
    @views xCoordinates[:,1:oldvertices] = oldCoordinates

    # refine faces
    newvertex = zeros(T,size(xCoordinates,1))
    newnodenr = oldvertices
    newnode4facemidpoint = zeros(Int,nfaces)
    nnodes4item = 0
    # add face midpoints to Coordinates
    for face = 1 : nfaces
        if facemarkers[face] == true
            nnodes4item = num_targets(oldFaceNodes,face)
            fill!(newvertex,0.0)
            for k = 1 : nnodes4item, d = 1 : size(xCoordinates,1)
                newvertex[d] += xCoordinates[d,oldFaceNodes[k,face]]
            end
            newvertex ./= nnodes4item
            newnodenr += 1
            newnode4facemidpoint[face] = newnodenr
            for d = 1 : size(xCoordinates,1)
                xCoordinates[d,newnodenr] = newvertex[d]
            end
        end
    end

    # determine new cells
    nnodes4item = 3
    nfaces4item = 3
    ncells = 0
    subitemnodes = zeros(K,3)
    m = 0
    localmarked_faces = zeros(Bool,3)
    refine_rule = RGB_refinement_rule(Triangle2D,localmarked_faces)
    nnewcells::Int = 0
    for cell = 1 : num_sources(oldCellNodes)
        for f = 1 : 3
            localmarked_faces[f] = facemarkers[oldCellFaces[f,cell]]
        end
        refine_rule = RGB_refinement_rule(Triangle2D,localmarked_faces)
        nnewcells = size(refine_rule,2)
        for j = 1 : nnewcells
            for k = 1 : 3
                m = refine_rule[k,j]
                if m <= nnodes4item
                    subitemnodes[k] = oldCellNodes[m,cell]
                elseif m <= nnodes4item + nfaces4item
                    subitemnodes[k] = newnode4facemidpoint[oldCellFaces[m-nnodes4item,cell]]
                end
            end
            append!(xCellNodes,view(subitemnodes,:))
            push!(xCellGeometries,Triangle2D)
            push!(xCellRegions,oldCellRegions[cell])
        end
        ncells += nnewcells
    end

    # assign new cells to grid
    xgrid[Coordinates] = xCoordinates
    if typeof(oldCellNodes) == Array{K,2}
        nnodes4item = size(oldCellNodes,1)
        xgrid[CellNodes] = reshape(xCellNodes.colentries,nnodes4item,num_sources(xCellNodes))
    else
        xgrid[CellNodes] = xCellNodes
    end
    if typeof(oldCellRegions) <: VectorOfConstants
        xgrid[CellRegions] = VectorOfConstants{K}(1,ncells)
    else
        xgrid[CellRegions] = xCellRegions
    end
    xgrid[CellGeometries] = VectorOfConstants(Triangle2D,ncells)

    # determine new boundary faces
    oldBFaceNodes = source_grid[BFaceNodes]
    oldBFaceFaces = source_grid[BFaceFaces]
    oldBFaceRegions = source_grid[BFaceRegions]

    xBFaceRegions = zeros(K,0)
    nbfaces = num_sources(oldBFaceNodes)
    xBFaceNodes = zeros(K,0)

    refine_rule = uniform_refine_rule(Edge1D)

    newbfaces = 0
    for bface = 1 : nbfaces
        face = oldBFaceFaces[bface]
        if facemarkers[face] == true
            for j = 1 : 2
                for k = 1 : 2
                    m = refine_rule[j,k]
                    if dim == 2
                        if m <= 2
                            subitemnodes[k] = oldBFaceNodes[m,bface]
                        else
                            subitemnodes[k] = newnode4facemidpoint[face]
                        end
                    end
                end
                append!(xBFaceNodes,view(subitemnodes,1:2))
                push!(xBFaceRegions,oldBFaceRegions[bface])
            end
            newbfaces +=1
        else
            append!(xBFaceNodes,view(oldBFaceNodes,:,bface))
            push!(xBFaceRegions,oldBFaceRegions[bface])
        end
    end
    @debug "bisected bfaces = $newbfaces"

    xgrid[BFaceNodes] = reshape(xBFaceNodes,(2,nbfaces+newbfaces))
    xgrid[BFaceRegions] = xBFaceRegions
    xgrid[BFaceGeometries] = VectorOfConstants(Edge1D,nbfaces+newbfaces)

    return xgrid
end


text:
export Loss, QuadLoss, LogisticLoss, PoissonLoss, NoLoss

abstract type Loss end

###############################
# Quadratic loss (Normal data)
###############################
struct QuadLoss <: Loss
    scale::Number
end

function compute_quadloss!(Z, A, tau)
    if size(Z, 1) == 0
        loss = 0.0
    else
        loss = 0.5*sum( (Z - A).^2 .* transpose(tau))
    end
    return loss
end

function compute_quadloss_delta!(Z, A)
    Z .-= A
    return nothing
end

###############################
# Logistic loss (binary data)
###############################
struct LogisticLoss <: Loss
    scale::Number
end

function compute_logloss!(Z, A)
    if size(Z, 1) == 0
        loss = 0.0
    else
        Z .= exp.(-Z)
        Z .+= 1.0
        Z .= 1.0 ./Z
        loss = -sum( A .* log.(1e-12 .+ Z) + (1.0 .- A).*log.(1e-12 + 1.0 .- Z) )
    end
    return loss
end

function compute_logloss_delta!(Z, A)
    Z .= exp.(-Z)
    Z .+= 1.0
    Z .= 1.0 ./ Z
    Z .-= A
    return nothing
end

###############################
# Poisson loss (count data)
###############################
struct PoissonLoss <: Loss
    scale::Number
end

function compute_poissonloss!(Z, A)
    if size(Z,2) == 0
        loss = 0.0
    else
        loss = sum(1.0.*Z.*A - exp.(1.0.*Z))
    end
    return loss
end

function compute_poissonloss_delta!(Z, A)
    Z .= exp.(Z)
    Z .*= -1.0
    Z .+= A
    return nothing
end

##########################################
## Ternary loss (tailored for CNA data)
##########################################
#
#struct TernaryLoss <: Loss
#    scale::Number
#end
#
#function compute_ternaryloss!(Z, A, tau)
#
#end

##################################
# NoLoss (for latent features)
##################################
struct NoLoss <: Loss
    scale::Number
end

function compute_noloss!(Z,A)
    return 0.0
end

function compute_noloss_delta!(Z, A)
    return nothing
end

#function vprint(s; verbose=false)
#    if verbose
#        println(s)
#    end
#end

##################################
# Other functions
##################################
function compute_loss!(X, Y, Z_ql, Z_ll, Z_pl,
                       A_ql, A_ll, A_pl,
                       beta, beta_reg,
                       mu, mu_reg, theta, theta_reg,
                       tau, a_0_tau, b_0_tau,
                       X_reg_mats, Y_reg_mats;
                       verbose=false)
    loss = compute_quadloss!(Z_ql, A_ql, tau)
    loss += compute_logloss!(Z_ll, A_ll)
    loss += compute_poissonloss!(Z_pl, A_pl)

    for i=1:size(beta,1)
        loss += 0.5 * dot(beta[i,:], beta_reg*beta[i,:])
    end
    loss += compute_mat_reg_loss(X, X_reg_mats)
    loss += compute_mat_reg_loss(Y, Y_reg_mats)
    loss += compute_vec_reg_loss(mu, mu_reg)
    loss += compute_vec_reg_loss(theta, theta_reg)

    M = size(X,2)
    loss += compute_tau_loss(tau, a_0_tau, b_0_tau, M)

    return loss
end

function compute_tau_loss(tau, a_0_tau, b_0_tau, M)
    return sum((b_0_tau.*tau) .- (0.5*M + a_0_tau - 1).*log.(tau))
end

function compute_grad_delta!(Z_ql, Z_ll, Z_pl,
                             A_ql, A_ll, A_pl)
    compute_quadloss_delta!(Z_ql, A_ql)
    compute_logloss_delta!(Z_ll, A_ll)
    compute_poissonloss_delta!(Z_pl, A_pl)
end

function compute_grad_X!(grad_X, X, Y, Z,
                         Z_ql, Z_ll, Z_pl,
                         A_ql, A_ll, A_pl,
                         inst_reg_mats)

    compute_grad_delta!(Z_ql, Z_ll, Z_pl,
                        A_ql, A_ll, A_pl)

    N = size(Z,2)
    grad_X .= (Y * transpose(Z))

    add_mat_reg_grad!(grad_X, X, inst_reg_mats)
end

function compute_grad_Y!(grad_Y, X, Y, Z,
                         Z_ql, Z_ll, Z_pl,
                         A_ql, A_ll, A_pl,
                         feat_reg_mats)

    compute_grad_delta!(Z_ql, Z_ll, Z_pl,
                        A_ql, A_ll, A_pl)

    M = size(Z,1)
    grad_Y .= X * Z

    add_mat_reg_grad!(grad_Y, Y, feat_reg_mats)
end

function compute_grad_beta!(grad_beta, beta, C, Z, tau, beta_reg)
    grad_beta .= (transpose(C)*Z) .* transpose(tau)
    grad_beta .+= transpose(beta_reg*transpose(beta))
end

function compute_grad_mu!(grad_mu, mu, Z, tau, mu_reg_mat)
    grad_mu .= Z*tau
    grad_mu .+= mu_reg_mat*mu
end

function compute_grad_theta!(grad_theta, theta, Z, tau, theta_reg_mat)
   grad_theta .= (sum(Z, dims=1)[1,:] .* tau)
   grad_theta .+= theta_reg_mat*theta
end

function update_tau!(tau, Z, a_0, b_0)
    M = size(Z, 1)
    tau .= (2.0*a_0 + M)./(2.0*b_0 .+ sum(Z.^2.0, dims=1)[1,:])
end

###########################
# PRECOMPILE
gpu_mat = CuArray{Float32,2}
gpu_vec = CuArray{Float32,1}
gpu_sparse = CUDA.CUSPARSE.CuSparseMatrixCSC{Float32}

precompile(compute_quadloss!, (gpu_mat, gpu_mat, gpu_vec))
precompile(compute_logloss!, (gpu_mat, gpu_mat))
precompile(compute_poissonloss!, (gpu_mat, gpu_mat))

precompile(compute_tau_loss, (gpu_vec, Float32, Float32, Int64))

precompile(compute_loss!, (gpu_mat, gpu_mat,
                           gpu_mat, gpu_mat, gpu_mat,
                           gpu_mat, gpu_mat, gpu_mat,
                           gpu_mat, gpu_sparse,
                           gpu_mat, gpu_sparse,
                           gpu_vec, gpu_sparse, gpu_vec, gpu_sparse,
                           gpu_vec, Float32, Float32,
                           Vector{gpu_sparse}, Vector{gpu_sparse}))

precompile(compute_grad_delta!, (gpu_mat, gpu_mat, gpu_mat,
                                 gpu_mat, gpu_mat, gpu_mat))
precompile(compute_quadloss_delta!, (gpu_mat, gpu_mat))
precompile(compute_logloss_delta!, (gpu_mat, gpu_mat))
precompile(compute_poissonloss_delta!, (gpu_mat, gpu_mat))


text:
@testset "Generating fast functions" begin
	ε = 100*eps(Float64)
	t1,t2,t3 = randn(3)

	f(x1,x2,x3) = 3*x3^2 + 2*x2*x1
	@polyvar x y z
	polynomial = 3*z^2 + 2*y*x
	g = @fastfunc polynomial
	@test isapprox(f(t1,t2,t3),g(t1,t2,t3),atol=ɛ)
	# Test inside function scope
	function useInsideFunctionScope1(polynomial,t1,t2,t3)
		g = @fastfunc polynomial
		@test isapprox(f(t1,t2,t3),Base.invokelatest(g, t1,t2,t3),atol=ɛ)
	end
	useInsideFunctionScope1(polynomial,t1,t2,t3)

	h = fastfunc(polynomial)
    @test isapprox(f(t1,t2,t3),h(t1,t2,t3),atol=ɛ)
	# Test inside function scope
	function useInsideFunctionScope2(polynomial,t1,t2,t3)
        g = fastfunc(polynomial)
		@test isapprox(f(t1,t2,t3),g(t1,t2,t3),atol=ɛ)
	end
    useInsideFunctionScope2(polynomial,t1,t2,t3)
end


text:
using DiffCheckNumbers

println("Testing Constructors")

t = @elapsed include("TestConstructor.jl")

println("done (took $t seconds).")

println("Testing Promotion and Conversion")

t = @elapsed include("TestPromotionConversion.jl")

println("done (took $t seconds).")

println("Testing Arithmetic")

t = @elapsed include("TestArithmetic.jl")

println("done (took $t seconds).")


text:
module SquareSymmetries

export unit, rotate90, rotate180, rotate270, flipdiag, flipadiag, flipx, flipy, symmetries

include("symmetry.jl")

end


text:
# Declaration ... what is the clean way to do that?
parse_date() = nothing
parse_float() = nothing

## FASTA headers
const augur_header_reference = ["strain", "Strain", "STRAIN", "name", "Name", "label", "Label"]
const augur_header_date = ["date", "Date"]
const augur_header_virus = ["virus", "Virus"]
const augur_minimal_fields = [:strain, :date, :virus]
# Fieldnames below are IGNORED when reading a fasta file
const ignored_header_fields = ["?", "", '?',:?]
const special_fields = ["date", :date, "Date", :date, :fitness]
const parse_special_field = Dict("date" => parse_date, :date => parse_date, :Date => parse_date, "Date" => parse_date, :fitness =>parse_float)
## The following could be added to a BioSequences fork
# const symbol_to_alphabet = Dict(BioSequences.AminoAcid => BioSequences.AminoAcidAlphabet,
								# BioSequences.DNA => BioSequences.DNAAlphabet,
								# BioSequences.RNA => BioSequences.RNAAlphabet)
symbol_to_alphabet(::Type{AminoAcid}) = BioSequences.AminoAcidAlphabet
symbol_to_alphabet(::Type{DNA}) = BioSequences.DNAAlphabet{4}
symbol_to_alphabet(::Type{RNA}) = BioSequences.RNAAlphabet{4}
ambiguous(::Type{AminoAcid}) = AA_X
ambiguous(::Type{DNA}) = DNA_N
ambiguous(::Type{RNA}) = RNA_N
isambiguous(::T) where T<:Real = false # Artificial sequences are never ambiguous
# datatypes
const sequenceymbols = (:aa, :rna, :dna, :int8, :int64, :bool, :artificial)
function type(x::Symbol)
	if x == :aa return AminoAcidAlphabet
	elseif x == :rna return RNAAlphabet{4}
	elseif x == :dna return DNAAlphabet{4}
	elseif x == :int8 return Int8
	elseif (x == :int64 || x == :artificial) return Int64
	elseif x == :bool return Bool
	else
		@error "Possible symbols: $(sequenceymbols)"
	end
end

##
## Errors
function unknown_seqtype()
	@error "Unknown sequence type `$seqtype`-- Possible values: (`:aa`, `:dna`, `:rna`)"
end

const aa2int = Dict(AA_Gap => 1, AA_A=>2, AA_C=>3, AA_D=>4, AA_E=>5, AA_F=>6, AA_G=>7, AA_H=>8, AA_I=>9, AA_K=>10, AA_L=>11, AA_M=>12, AA_N=>13, AA_P=>14, AA_Q=>15, AA_R=>16, AA_S=>17, AA_T=>18, AA_V=>19, AA_W=>20, AA_Y=>21)
const int2aa = Dict(i => findfirst(==(i), BioTools.aa2int) for i in 1:21)
# function int2aa(a::Integer)
# 	return _aa_alphabet[a]
# end

max_stars_repo_name: PierreBarrat/FluPredictibility

max_issues_repo_name: PierreBarrat/FluPredictibility

max_forks_repo_name: PierreBarrat/FluPredictibility


text:
module ArrayAwareCSE

using Yolk
import Yolk: allow, YolkOptimizer
using InteractiveUtils
using CodeInfoTools

using Metatheory
using Metatheory.Library
using Metatheory.EGraphs

using BenchmarkTools

@metatheory_init()

#####
##### Array-aware CSE.
#####

using LinearAlgebra

function f(x1::Vector{Float64})
    x2  = x1
    x3  = x1 + x2
    x4  = x1 + x2 + x3
    x5  = x1 + x2 + x3 + x4
    #x6  = x1 + x2 + x3 + x4 + x5
    #x7  = x1 + x2 + x3 + x4 + x5 + x6
    #x8  = x1 + x2 + x3 + x4 + x5 + x6 + x7
    #x9  = x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8
    #x10 = x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9
    #x11 = x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10
    #x12 = x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8 + x9 + x10 + x11
    return x5
end

# Manual.

th = @methodtheory begin
    a + a => 2 * a
    +(a, b, c) => (a + b) + c
    +(a, b, c, d, e) => (a + b) + (c + d + e)
    a * (b * c) => (a * b) * c
    (k * a) + (n * a) => (k + n) * a
    a + (k::Number * a) => (k + 1) * a
    a::Number * b::Number |> a * b
end
display(th)

allow(::YolkOptimizer, ::typeof(*)) = true
allow(::YolkOptimizer, ::typeof(+)) = true

# Optimize.
src = opt(f, Tuple{Vector{Float64}}; ctx = YolkOptimizer(th), opt = true)
display(src)

end # module


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related to the conversion from DCM to Euler angles.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export dcm_to_angle

"""
    dcm_to_angle(dcm::DCM, rot_seq::Symbol=:ZYX)

Convert the `dcm` to Euler Angles (see [`EulerAngles`](@ref)) given a rotation
sequence `rot_seq`.

The rotation sequence is defined by a `:Symbol`. The possible values are:
`:XYX`, `XYZ`, `:XZX`, `:XZY`, `:YXY`, `:YXZ`, `:YZX`, `:YZY`, `:ZXY`, `:ZXZ`,
`:ZYX`, and `:ZYZ`. If no value is specified, then it defaults to `:ZYX`.

# Gimbal-lock and special cases

If the rotations are about three different axes, *e.g.* `:XYZ`, `:ZYX`, etc.,
then a second rotation of `±90˚` yields a gimbal-lock. This means that the
rotations between the first and third axes have the same effect. In this case,
the net rotation angle is assigned to the first rotation, and the angle of the
third rotation is set to 0.

If the rotations are about two different axes, *e.g.* `:XYX`, `:YXY`, etc., then
a rotation about the duplicated axis yields multiple representations. In this
case, the entire angle is assigned to the first rotation and the third rotation
is set to 0.

# Example

```jldoctest
julia> D = DCM([1. 0. 0.; 0. 0. -1; 0. -1 0.]);

julia> dcm_to_angle(D,:XYZ)
EulerAngles{Float64}:
  R(X) :  1.5708 rad  ( 90.0°)
  R(Y) :  0.0    rad  ( 0.0°)
  R(Z) :  0.0    rad  ( 0.0°)

julia> D = angle_to_dcm(1, -pi / 2, 2, :ZYX);

julia> dcm_to_angle(D, :ZYX)
EulerAngles{Float64}:
  R(Z) :  3.0    rad  ( 171.887°)
  R(Y) : -1.5708 rad  (-90.0°)
  R(X) :  0.0    rad  ( 0.0°)

julia> D = angle_to_dcm(1, :X) * angle_to_dcm(2, :X);

julia> dcm_to_angle(D, :XYX)
EulerAngles{Float64}:
  R(X) :  3.0 rad  ( 171.887°)
  R(Y) :  0.0 rad  ( 0.0°)
  R(X) :  0.0 rad  ( 0.0°)
```
"""
function dcm_to_angle(dcm::DCM{T}, rot_seq::Symbol=:ZYX) where T<:Number
    if rot_seq == :ZYX
        # Check for singularities.
        if !(abs(dcm[1, 3]) ≥ 1 - eps())
            return EulerAngles(
                _mod_atan(+dcm[1, 2], +dcm[1, 1]),
                asin(-dcm[1, 3]),
                _mod_atan(+dcm[2, 3], +dcm[3, 3]),
                rot_seq
            )
        else
            return EulerAngles(
                _mod_atan(-dcm[2, 1], +dcm[2, 2]),
                _mod_asin(-dcm[1, 3]),
                T(0),
                rot_seq
            )
        end
    elseif rot_seq == :XYX
        # Check for singularities.
        if !(abs(dcm[1, 1]) ≥ 1 - eps())
            return EulerAngles(
                _mod_atan(+dcm[1, 2], -dcm[1, 3]),
                acos(+dcm[1, 1]),
                _mod_atan(+dcm[2, 1], +dcm[3, 1]),
                rot_seq
            )
        else
            return EulerAngles(
                _mod_atan(+dcm[2, 3], +dcm[2, 2]),
                _mod_acos(dcm[1, 1]),
                T(0),
                rot_seq
            )
        end
    elseif rot_seq == :XYZ
        # Check for singularities.
        if !(abs(dcm[3, 1]) ≥ 1 - eps())
            return EulerAngles(
                _mod_atan(-dcm[3, 2], +dcm[3, 3]),
                asin(+dcm[3, 1]),
                _mod_atan(-dcm[2, 1], +dcm[1, 1]),
                rot_seq
            )
        else
            return EulerAngles(
                _mod_atan(+dcm[2, 3], +dcm[2, 2]),
                _mod_asin(+dcm[3, 1]),
                T(0),
                rot_seq
            )
        end
    elseif rot_seq == :XZX
        # Check for singularities.
        if !(abs(dcm[1, 1]) ≥ 1 - eps())
            return EulerAngles(
                _mod_atan(+dcm[1, 3], +dcm[1, 2]),
                acos(+dcm[1, 1]),
                _mod_atan(+dcm[3, 1], -dcm[2, 1]),
                rot_seq
            )
        else
            return EulerAngles(
                _mod_atan(-dcm[3, 2], +dcm[3, 3]),
                _mod_acos(dcm[1, 1]),
                T(0),
                rot_seq
            )
        end
    elseif rot_seq == :XZY
        # Check for singularities.
        if !(abs(dcm[2, 1]) ≥ 1 - eps())
            return EulerAngles(
                _mod_atan(+dcm[2, 3], +dcm[2, 2]),
                asin(-dcm[2, 1]),
                _mod_atan(+dcm[3, 1], +dcm[1, 1]),
                rot_seq
            )
        else
            return EulerAngles(
                _mod_atan(-dcm[3, 2], +dcm[3, 3]),
                _mod_asin(-dcm[2, 1]),
                T(0),
                rot_seq
            )
        end
    elseif rot_seq == :YXY
        # Check for singularities.
        if !(abs(dcm[2, 2]) ≥ 1 - eps())
            return EulerAngles(
                _mod_atan(+dcm[2, 1], +dcm[2, 3]),
                acos(+dcm[2, 2]),
                _mod_atan(+dcm[1, 2], -dcm[3, 2]),
                rot_seq
            )
        else
            return EulerAngles(
                _mod_atan(-dcm[1, 3], +dcm[1, 1]),
                _mod_acos(dcm[2, 2]),
                T(0),
                rot_seq
            )
        end
    elseif rot_seq == :YXZ
        if !(abs(dcm[3, 2]) ≥ 1 - eps())
            return EulerAngles(
                _mod_atan(+dcm[3, 1], +dcm[3, 3]),
                asin(-dcm[3, 2]),
                _mod_atan(+dcm[1, 2], +dcm[2, 2]),
                rot_seq
            )
        else
            return EulerAngles(
                _mod_atan(-dcm[1, 3], +dcm[1, 1]),
                _mod_asin(-dcm[3, 2]),
                T(0),
                rot_seq
            )
        end
    elseif rot_seq == :YZX
        # Check for singularities.
        if !(abs(dcm[1, 2]) ≥ 1 - eps())
            return EulerAngles(
                _mod_atan(-dcm[1, 3], +dcm[1, 1]),
                asin(+dcm[1, 2]),
                _mod_atan(-dcm[3, 2], +dcm[2, 2]),
                rot_seq
            )
        else
            return EulerAngles(
                _mod_atan(+dcm[3, 1], +dcm[3, 3]),
                _mod_asin(+dcm[1, 2]),
                T(0),
                rot_seq
            )
        end
    elseif rot_seq == :YZY
        # Check for singularities.
        if !(abs(dcm[2, 2]) ≥ 1 - eps())
            return EulerAngles(
                _mod_atan(+dcm[2, 3], -dcm[2, 1]),
                acos(+dcm[2, 2]),
                _mod_atan(+dcm[3, 2], +dcm[1, 2]),
                rot_seq
            )
        else
            return EulerAngles(
                _mod_atan(dcm[3, 1], dcm[3, 3]),
                _mod_acos(dcm[2, 2]),
                T(0),
                rot_seq
            )
        end
    elseif rot_seq == :ZXY
        # Check for singularities.
        if !(abs(dcm[2, 3]) ≥ 1 - eps())
            return EulerAngles(
                _mod_atan(-dcm[2, 1], +dcm[2, 2]),
                asin(+dcm[2, 3]),
                _mod_atan(-dcm[1, 3], +dcm[3, 3]),
                rot_seq
            )
        else
            return EulerAngles(
                _mod_atan(+dcm[1, 2], +dcm[1, 1]),
                _mod_asin(+dcm[2, 3]),
                T(0),
                rot_seq
            )
        end
    elseif rot_seq == :ZXZ
        # Check for singularities.
        if !(abs(dcm[3, 3]) ≥ 1 - eps())
            return EulerAngles(
                _mod_atan(+dcm[3, 1], -dcm[3, 2]),
                acos(+dcm[3, 3]),
                _mod_atan(+dcm[1, 3], +dcm[2, 3]),
                rot_seq
            )
        else
            return EulerAngles(
                _mod_atan(dcm[1, 2], dcm[1, 1]),
                _mod_acos(dcm[3, 3]),
                T(0),
                rot_seq
            )
        end
    elseif rot_seq == :ZYZ
        # Check for singularities.
        if !(abs(dcm[3, 3]) ≥ 1 - eps())
            return EulerAngles(
                _mod_atan(+dcm[3, 2], +dcm[3, 1]),
                acos(+dcm[3, 3]),
                _mod_atan(+dcm[2, 3], -dcm[1, 3]),
                rot_seq
            )
        else
            return EulerAngles(
                _mod_atan(-dcm[2, 1], dcm[2, 2]),
                _mod_acos(dcm[3, 3]),
                T(0),
                rot_seq
            )
        end
    else
        throw(ArgumentError("The rotation sequence :$rot_seq is not valid."))
    end
end

################################################################################
#                              Private functions
################################################################################

# This modified function computes exactly what `atan(y,x)` computes except that
# it will neglect signed zeros. Hence:
#
#   _mod_atan(0.0, -0.0) = _mod_atan(-0.0, 0.0) = 0.0
#
# The signed zero can lead to problems when converting from DCM to Euler angles.
_mod_atan(y::T, x::T) where T<:Number = atan(y + T(0), x + T(0))

# This modified function computes the `acos(x)` if `|x| <= 1` and computes
# `acos( sign(x) )`  if `|x| > 1` to avoid numerical errors when converting DCM
# to  Euler Angles.
function _mod_acos(x::T) where T<:Number
    if x > 1
        return float(T(0))
    elseif x < -1
        return float(T(π))
    else
        return acos(x)
    end
end

# This modified function computes the `asin(x)` if `|x| <= 1` and computes
# `asin( sign(x) )`  if `|x| > 1` to avoid numerical errors when converting DCM
# to  Euler Angles.
function _mod_asin(x::T) where T<:Number
    if x > 1
        return +float(T(π / 2))
    elseif x < -1
        return -float(T(π / 2))
    else
        return asin(x)
    end
end


text:
sir_petri = LabelledPetriNet([:S,:I,:R], :inf=>((:S,:I)=>(:I,:I)), :rec=>(:I=>:R))

sir = infection ⋅ recovery

@test apex(sir) == sir_petri

sir_relation = @relation (s,i,r) begin
    infection(s,i)
    recovery(i,r)
end

@test apex(oapply_epi(sir_relation)) == sir_petri

max_stars_repo_name: AlgebraicJulia/AlgebraicPetri

max_issues_repo_name: AlgebraicJulia/AlgebraicPetri

max_forks_repo_name: AlgebraicJulia/AlgebraicPetri


text:
using Makie, Colors, Observables
r = range(0, stop=5pi, length=100)
scene = lines(r, sin.(r), linewidth = 3)
lineplot = scene[end]
visible = node(:visible, false)
poprect = lift(scene.events.mouseposition) do mp
    FRect((mp .+ 5), 250, 40)
end
textpos = lift(scene.events.mouseposition) do mp
    Vec3f0((mp .+ 5 .+ (250/2, 40 / 2))..., 120)
end
popup = poly!(campixel(scene), poprect, raw = true, color = :white, strokewidth = 2, strokecolor = :black, visible = visible)
rect = popup[end]
translate!(rect, Vec3f0(0, 0, 100))
text!(popup, "( 0.000,  0.000)", textsize = 30, position = textpos, color = :darkred, align = (:center, :center), raw = true, visible = visible)
text_field = popup[end]
scene

x = Node(false)
on(scene.events.mouseposition) do event
    plot, idx = Makie.mouse_selection(scene)
    if plot == lineplot && idx > 0
        visible[] = true
        text_field[1] = sprint(io-> print(io, round.(Float64.(Tuple(lineplot[1][][idx])), digits = 3)))
    else
        visible[] = false
    end
    return
end
scene

using Makie
using LinearAlgebra
img = rand(100, 100)
scene = Scene()
heatmap!(scene, img, scale_plot = false)
clicks = Node(Point2f0[(0, 0)])
blues = Node(Point2f0[])

on(scene.events.mousebuttons) do buttons
    if ispressed(scene, Mouse.left)
        pos = to_world(scene, Point2f0(scene.events.mouseposition[]))
        found = -1
        c = clicks[]
        for i in 1:length(c)
           if norm(pos - c[i]) < 1
               found = i
           end
        end
        if found >= 1
            blues[] = push!(blues[], pos)
            deleteat!(clicks[], found)
        else
            push!(clicks[], pos)
        end
        clicks[] = clicks[]
   end
   return
end

scatter!(scene, clicks, color = :red, marker = '+', markersize = 10, raw = true)
red_clicks = scene[end]
scatter!(scene, blues, color = :blue, marker = 'o', markersize = 10, raw = true)
scene

using Makie
points = node(:poly, Point2f0[(0, 0), (0.5, 0.5), (1.0, 0.0)])

scene = poly(points, strokewidth = 2, strokecolor = :black, color = :skyblue2, show_axis = false, scale_plot = false)
scatter!(points, color = :white, strokewidth = 10, markersize = 0.05, strokecolor = :black, raw = true)
pplot = scene[end]
push!(points[], Point2f0(0.6, -0.3))
points[] = points[]

function add_move!(scene, points, pplot)
    idx = Ref(0); dragstart = Ref(false); startpos = Base.RefValue(Point2f0(0))
    on(events(scene).mousedrag) do drag
        if ispressed(scene, Mouse.left)
            if drag == Mouse.down
                plot, _idx = Makie.mouse_selection(scene)
                if plot == pplot
                    idx[] = _idx; dragstart[] = true
                    startpos[] = to_world(scene, Point2f0(scene.events.mouseposition[]))
                end
            elseif drag == Mouse.pressed && dragstart[] && checkbounds(Bool, points[], idx[])
                pos = to_world(scene, Point2f0(scene.events.mouseposition[]))
                points[][idx[]] = pos
                points[] = points[]
            end
        else
            dragstart[] = false
        end
        return
    end
end

function add_remove_add!(scene, points, pplot)
    on(events(scene).mousebuttons) do but
        if ispressed(but, Mouse.left) && ispressed(scene, Keyboard.left_control)
            pos = to_world(scene, Point2f0(events(scene).mouseposition[]))
            push!(points[], pos)
            points[] = points[]
        elseif ispressed(but, Mouse.right)
            plot, idx = Makie.mouse_selection(scene)
            if plot == pplot && checkbounds(Bool, points[], idx)
                deleteat!(points[], idx)
                points[] = points[]
            end
        end
        return
    end
end
add_move!(scene, points, pplot)
add_remove_add!(scene, points, pplot)
center!(scene)
scene

using Makie, Colors
using AbstractPlotting: transformationmatrix, textslider, colorswatch, hbox!
s1, a = textslider(0f0:50f0, "a", start = 13)
s2, b = textslider(-20f0:20f0, "b", start = 10)
s3, c = textslider(0f0:20f0, "c", start = 2)
s4, d = textslider(range(0.0, stop = 0.02, length = 100), "d", start = 0.01)
s5, scales = textslider(range(0.01, stop = 0.5, length = 100), "scale", start = 0.1)
s6, colorsw, pop = colorswatch()

function lorenz(t0, a, b, c, h)
    Point3f0(
        t0[1] + h * a * (t0[2] - t0[1]),
        t0[2] + h * (t0[1] * (b - t0[3]) - t0[2]),
        t0[3] + h * (t0[1] * t0[2] - c * t0[3]),
    )
end
# step through the `time`
function lorenz(array::Vector, a = 5.0 ,b = 2.0, c = 6.0, d = 0.01)
    t0 = Point3f0(0.1, 0, 0)
    for i = eachindex(array)
        t0 = lorenz(t0, a,b,c,d)
        array[i] = t0
    end
    array
end

n1, n2 = 18, 30
N = n1*n2
args_n = (a, b, c, d)
v0 = lorenz(zeros(Point3f0, N), to_value.(args_n)...)
positions = lift(lorenz, Node(v0), args_n...)
rotations = lift(diff, positions)
rotations = lift(x-> push!(x, x[end]), rotations)

mesh_scene = meshscatter(
    positions,
    markersize = scales, rotation = rotations,
    intensity = collect(range(0f0, stop = 1f0, length = length(positions[]))),
    color = colorsw
)

vbox(
    hbox(s1, s2, s3, s4, s5, s6),
    mesh_scene
)

function record_events(f, scene, path)
    display(scene)
    result = Vector{Pair{Float64, Pair{Symbol, Any}}}()
    for field in fieldnames(Events)
        foreach(getfield(scene.events, field)) do value
            value = isa(value, Set) ? copy(value) : value
            push!(result, time() => (field => value))
        end
    end
    f()

    open(path, "w") do io
        serialize(io, result)
    end
end
record_events(scene, "test.jls") do
    wait(Makie.global_gl_screen())
end

function replay(scene, path)
    display(scene)
    events = open(io-> deserialize(io), path)
    sort!(events, by = first)
    for i in 1:length(events)
        t1, (field, value) = events[i]
        field == :mousedrag && continue
        if field == :mousebuttons
            getfield(scene.events, field)[] = value
        else
            getfield(scene.events, field)[] = value
        end
        yield(); force_update!()
        if i < length(events)
            t2, (field, value) = events[i + 1]
            # min sleep time 0.001
            (t2 - t1 > 0.001) && sleep(t2 - t1)
        end
    end
end

replay(scene, "test.jls")
open(io-> serialize(io, Set([Mouse.left])) , "test2.jls", "w")
open(io-> deserialize(io) , "test2.jls")

scatter!(campixel(scene), lift(x-> [x], scene.events.mouseposition),color = :black, markersize = 20, raw = true)
scene


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

@doc raw"""EnvVar represents an environment variable present in a Container.

    IoK8sApiCoreV1EnvVar(;
        name=nothing,
        value=nothing,
        valueFrom=nothing,
    )

    - name::String : Name of the environment variable. Must be a C_IDENTIFIER.
    - value::String : Variable references $(VAR_NAME) are expanded using the previous defined environment variables in the container and any service environment variables. If a variable cannot be resolved, the reference in the input string will be unchanged. The $(VAR_NAME) syntax can be escaped with a double $$, ie: $$(VAR_NAME). Escaped references will never be expanded, regardless of whether the variable exists or not. Defaults to \"\".
    - valueFrom::IoK8sApiCoreV1EnvVarSource : Source for the environment variable's value. Cannot be used if value is not empty.
"""
mutable struct IoK8sApiCoreV1EnvVar <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    value::Any # spec type: Union{ Nothing, String } # spec name: value
    valueFrom::Any # spec type: Union{ Nothing, IoK8sApiCoreV1EnvVarSource } # spec name: valueFrom

    function IoK8sApiCoreV1EnvVar(;name=nothing, value=nothing, valueFrom=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1EnvVar, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiCoreV1EnvVar, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        validate_property(IoK8sApiCoreV1EnvVar, Symbol("valueFrom"), valueFrom)
        setfield!(o, Symbol("valueFrom"), valueFrom)
        o
    end
end # type IoK8sApiCoreV1EnvVar

const _property_map_IoK8sApiCoreV1EnvVar = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("value")=>Symbol("value"), Symbol("valueFrom")=>Symbol("valueFrom"))
const _property_types_IoK8sApiCoreV1EnvVar = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("value")=>"String", Symbol("valueFrom")=>"IoK8sApiCoreV1EnvVarSource")
Base.propertynames(::Type{ IoK8sApiCoreV1EnvVar }) = collect(keys(_property_map_IoK8sApiCoreV1EnvVar))
Swagger.property_type(::Type{ IoK8sApiCoreV1EnvVar }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1EnvVar[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1EnvVar }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1EnvVar[property_name]

function check_required(o::IoK8sApiCoreV1EnvVar)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1EnvVar }, name::Symbol, val)
end


text:
using FileIO

const R_FILE_PATH = r"(((\.\.){1}/)*|(/){1})?(([a-zA-Z0-9]*)/)*([a-zA-Z0-9]*)+([.dat]|[.txt]|[.dsv]|[.csv])+"
function issomepath(some_path)
    return occursin(R_FILE_PATH, some_path)
end

file_extension(path::AbstractString) = try match(r"\.[A-Za-z0-9]+$", path).match catch e "" end

function file_format(path::AbstractString)
    extension = strip(file_extension(path))
    if isempty(extension)
        return format""
    else
        @assert startswith(extension, '.')
        @assert length(extension) > 1
        return DataFormat{Symbol(uppercase(extension[2:end]))}
    end
end


text:
module MOEA

using Random, LinearAlgebra, Statistics
using TOML
using JLD
using UnPack: @unpack
using StackViews
using Plots
using Distributions
using NLSolversBase:
    NLSolversBase,
    AbstractObjective,
    ConstraintBounds,
    AbstractConstraints,
    nconstraints_x,
    nconstraints
using Infiltrator: @infiltrate
using Colors
import NLSolversBase: f_calls, value, value!
import Base:
    show,
    copy,
    minimum,
    summary,
    identity,
    getproperty,
    rand,
    getindex,
    length,
    copyto!,
    setindex!,
    replace,
    print,
    sortperm,
    sortperm!
import Distances:
    pairwise,
    Euclidean
using CImGui
using CImGui.GLFWBackend
using CImGui.OpenGLBackend
using CImGui.GLFWBackend.GLFW
using CImGui.OpenGLBackend.ModernGL
using CImGui.CSyntax
using CImGui.CSyntax.CStatic
using ImGuiGLFWBackend
using ImGuiOpenGLBackend
using LibCImGui
using GLFW
using Images
using TestImages
using ColorTypes

export optimize,
    ### TYPE
    Individual,
    Population,
    Objective,
    variables,
    objectives,
    constraints,
    ### CONSTRAINT
    BoxConstraints,
    ### METRIC
    gd,
    igd,
    FPL,
    ### TEST PROBLEM
    ZDT1,
    ZDT2,
    ZDT3,
    ZDT4,
    ZDT6,
    truepf,
    sphere,
    ackley,
    ### ALGORITHM
    GA,
    DE,
    NSGAII,
    MOEAD,
    MOEADDE,
    SMSEMOA,
    pfront,
    ### SELECTION
    ranklinear,
    uniformranking,
    roulette,
    rouletteinv,
    sus,
    susinv,
    tournament,
    truncation,
    uniformranking,
    ### CROSSOVER
    # binary
    SPX,  # single point
    TPX,  # two points
    SHFX, # shuffle
    UX,   # uniform
    BINX, # binary
    EXPX, # exponential
    BSX,  # binary subset
    # real value
    DC,   # discrete
    AX,   # average
    WAX,  # weighted average
    IC,   # intermediate
    LC,   # line
    HX,   # heuristic
    LX,   # Laplace
    MILX, # mixed integer Laplace
    SBX,  # simulated binary
    ### MUTATION
    ## GA
    # binary
    flip,
    inversion,
    # real value
    uniform,
    BGA,  # domain range
    PM,   # power
    MIPM, # mixed-integer power
    PLM,  # polynomial
    # combinatorial
    insertion,
    swap2,
    scramble,
    shifting

include("state.jl")
include("optimizer.jl")
include("individual.jl")
include("population.jl")
include("objective.jl")
include("options.jl")
include("termination.jl")
include("result.jl")
include("constraints.jl")
include("optimize.jl")
include("FPL.jl")

include("algorithms/selection.jl")
include("algorithms/crossover.jl")
include("algorithms/mutation.jl")
include("algorithms/moea.jl")

include("algorithms/GA.jl")
include("algorithms/DE.jl")
include("algorithms/NSGA-II.jl")
include("algorithms/MOEA-D.jl")
include("algorithms/MOEA-D-DE.jl")
include("algorithms/SMS-EMOA.jl")

include("problems/type.jl")
include("problems/ZDT.jl")
include("problems/DTLZ.jl")
include("problems/CEC2018.jl")
include("problems/sphere.jl")
include("problems/ackley.jl")

include("utilities.jl")

include("gui/utils.jl")
include("gui/image.jl")
include("gui/parameter.jl")
include("gui/problem.jl")
include("gui/algorithm.jl")
include("gui/playground.jl")
include("gui/experiment.jl")
include("gui/gui.jl")

end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Ninja_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Ninja")
JLLWrappers.@generate_main_file("Ninja", UUID("76642167-d241-5cee-8c94-7a494e8cb7b7"))
end  # module Ninja_jll


text:
include("./const.jl")
#include("./Preprocessing/runtests.jl")
#include("./PreliminarySelection/runtests.jl")
#include("./CrossValidation/runtests.jl")


text:
if VERSION < v"1.1"
    isnothing(x) = x === nothing
end

if VERSION < v"1.2"
    startswith(msg, x) = Base.startswith(msg, x)
    startswith(x) = msg->startswith(msg, x)
end


text:
# ---
# title: 936. Stamping The Sequence
# id: problem936
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: String, Greedy
# link: <https://leetcode.com/problems/stamping-the-sequence/description/>
# hidden: true
# ---
#
# You want to form a `target` string of **lowercase letters**.
#
# At the beginning, your sequence is `target.length` `'?'` marks.  You also have
# a `stamp` of lowercase letters.
#
# On each turn, you may place the stamp over the sequence, and replace every
# letter in the sequence with the corresponding letter from the stamp.  You can
# make up to `10 * target.length` turns.
#
# For example, if the initial sequence is "?????", and your stamp is `"abc"`,
# then you may make "abc??", "?abc?", "??abc" in the first turn.  (Note that the
# stamp must be fully contained in the boundaries of the sequence in order to
# stamp.)
#
# If the sequence is possible to stamp, then return an array of the index of the
# left-most letter being stamped at each turn.  If the sequence is not possible
# to stamp, return an empty array.
#
# For example, if the sequence is "ababc", and the stamp is `"abc"`, then we
# could return the answer `[0, 2]`, corresponding to the moves "?????" ->
# "abc??" -> "ababc".
#
# Also, if the sequence is possible to stamp, it is guaranteed it is possible to
# stamp within `10 * target.length` moves.  Any answers specifying more than
# this number of moves will not be accepted.
#
#
#
# **Example 1:**
#
#
#
#     Input: stamp = "abc", target = "ababc"
#     Output: [0,2]
#     ([1,0,2] would also be accepted as an answer, as well as some other answers.)
#
#
# **Example 2:**
#
#
#
#     Input: stamp = "abca", target = "aabcaca"
#     Output: [3,0,1]
#
#
#
#
# **Note:**
#
#   1. `1 <= stamp.length <= target.length <= 1000`
#   2. `stamp` and `target` only contain lowercase letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
function tensor_init_xeb(; enable_cuda = false)
    tn = [1.0, 0, 0]
    enable_cuda && return CuArray(tn)
    return tn
end

function tensor_single_xeb(; enable_cuda = false)
    tn = 1/2 * [1.0 1 1; 1 1 -1; 1 -1 1]
    enable_cuda && return CuArray(tn)
    return tn
end

function tensor_measure_xeb(; enable_cuda = false)
    tn = [2.0; 0; 0]
    enable_cuda && return CuArray(tn)
    return tn
end

function tensor_id_xeb(; enable_cuda = false)
    tn = [1.0 0 0; 0 1 0; 0 0 1]
    enable_cuda && return CuArray(tn)
    return tn
end
function tensor_noise_xeb(; enable_cuda = false)
    tn = [sqrt(1/2), 0, sqrt(1/2)]
    enable_cuda && return CuArray(tn)
    return tn
end

function tensor_fsim_xeb(; enable_cuda = false)
    tn = zeros(3, 3, 3, 3)
    for i = 1:3, j = 1:3
        tn[i, j, j, i] = 1.0
    end
    tn[2, 3, 3, 2] = -sqrt(3)/2
    tn[3, 2, 2, 3] = -sqrt(3)/2
    enable_cuda && return CuArray(tn)
    return tn
end

"""
    to_ein_code_xeb(layout, [s; haar = false, enable_cuda = false])

Generate the tensor network description for `layout` on a subset `s`. If the is a FSim gate out side `s`,
it will be replaced by an identity gate.
"""
function to_ein_code_xeb(layout::RQCLayout{VT, PT, SCGate}, s = collect(vertices(layout)); haar = false, enable_cuda = false) where {VT, PT}
    layout = deepcopy(layout)
    # update_id!(layout)
    vs = sort!(intersect(s, collect(vertices(layout))))
    nbits = length(vs)
    ids_size = Dict{Int, Int}(i => 3 for i in 1:nbits)
    ids_loc = Dict(i => (vs[i], 0) for i = 1:nbits)
    open_ids_dict = Dict(vs[i] => i for i = 1:nbits)
    tensor_ids = [[id] for id in 1:nbits]
    max_id = nbits
    tensors = Any[]
    for v in vs
        if gates(layout, v)[1] === Noise
            push!(tensors, sqrt(1/2)*tensor_noise_xeb(enable_cuda = enable_cuda))
        else
            push!(tensors, tensor_init_xeb(enable_cuda = enable_cuda))
        end
    end

    D = length(gates(layout, 1))
    for d = 1:D
        for v in vs
            if gates(layout, v)[d] === Single
                ts = tensor_single_xeb(; enable_cuda = enable_cuda)
                max_id += 1
                new_id_v = max_id
                push!(tensors, ts)
                push!(tensor_ids, [open_ids_dict[v], new_id_v])
                ids_size[new_id_v] = 3
                ids_loc[new_id_v] = (v, d)
                open_ids_dict[v] = new_id_v
            elseif gates(layout, v)[d] === Id
                ts = tensor_id_xeb(; enable_cuda = enable_cuda)
                max_id += 1
                new_id_v = max_id
                push!(tensors, ts)
                push!(tensor_ids, [open_ids_dict[v], new_id_v])
                ids_size[new_id_v] = 3
                ids_loc[new_id_v] = (v, d)
                open_ids_dict[v] = new_id_v
            elseif gates(layout, v)[d] === Noise
                max_id += 1
                new_id_v = max_id

                push!(tensors, tensor_noise_xeb(enable_cuda = enable_cuda), tensor_noise_xeb(enable_cuda = enable_cuda))
                push!(tensor_ids, [open_ids_dict[v]], [new_id_v])
                ids_size[new_id_v] = 3
                ids_loc[new_id_v] = (v, d)
                open_ids_dict[v] = new_id_v
            elseif gates(layout, v)[d] === FSim
                c = depth_to_cycle(d)
                u = partner(layout, v, c)
                if u < v && u in vs
                    max_id += 2
                    new_id_u = max_id - 1
                    new_id_v = max_id
                    push!(tensors, tensor_fsim_xeb(enable_cuda = enable_cuda))
                    push!(tensor_ids, [open_ids_dict[u], open_ids_dict[v], new_id_u, new_id_v])
                    ids_size[new_id_u] = 3
                    ids_size[new_id_v] = 3
                    ids_loc[new_id_u] = (u, d)
                    ids_loc[new_id_v] = (v, d)
                    open_ids_dict[u] = new_id_u
                    open_ids_dict[v] = new_id_v
                end
            end
        end
    end
    for v in vs
        push!(tensors, tensor_measure_xeb(enable_cuda = enable_cuda))
        push!(tensor_ids, [open_ids_dict[v]])
        open_ids_dict[v] = 0
    end

    open_ids = Int[]
    for v in vs
        open_ids_dict[v] > 0 && push!(open_ids, open_ids_dict[v])
    end
    # return tensor_ids, open_ids
    return EinCode(tensor_ids, open_ids), tensors, ids_size, ids_loc
end


text:
abstract type AbstractFillingObject end
abstract type AbstractOutlineObject end
abstract type AbstractWall end

"""
    Point{T}

Mutable struct that contains x and y coordinates. Each
coordinate is of type T.
"""
mutable struct Point{T}
    x::T
    y::T
end

"""
    LineWall{T}

Struct modeling a wall defined between two points. Each
point is of type T.
"""
struct LineWall{T} <: AbstractWall
    nodes::Vector{Point{T}}     # 2 points, defining start and end
    thickness::T
    n::Vector{T}                # normal unit vector of wall
    t::Vector{T}                # tangent unit vector of wall

    function LineWall(n1::Point{T},n2::Point{T},thickness::T) where T<:Real
        n = zeros(Float64,2); t = zeros(Float64,2)

        # compute wall length
        WL = sqrt((n2.x - n1.x)^2 + (n2.y - n1.y)^2)

        # comute tangent
        t[1] = (n2.x - n1.x)/WL
        t[2] = (n2.y - n1.y)/WL

        n[1] = -t[2]
        n[2] =  t[1]

        nodes = [n1,n2]

        return new{T}(nodes,thickness,n,t)
    end
    LineWall(n1,n2) = LineWall(n1,n2,0.0)
end

"""
    CircleWall{T}

Struct modeling a wall defined as a circle, with a center
and a radius. The base unit is of type T.
"""
struct CircleWall{T} <: AbstractWall
    center::Point{T}          # point that defines center
    radius::T                 # radius of circle
    thickness::T              # gives circle curve "area"

    function CircleWall(center::Point{T},radius::T) where T<:Real
        return new{T}(center,radius,zero(T))
    end
end

"""
    OutlineCircle{T}

Struct modeling an outline circle. The base unit is of type T.
"""
struct OutlineCircle{T} <: AbstractOutlineObject
    radius::T
    center::Point{T}
    wlist::Vector{CircleWall{T}}

    function OutlineCircle(radius::T,center::Point{T}) where T<:Real
        wlist = [CircleWall(center,radius)]
        return new{T}(radius,center,wlist)
    end
end

"""
    OutlineRectangle{T}

Struct modeling an outline rectangle, defined by
a lower left and upper right point. The base unit is of type T.
"""
struct OutlineRectangle{T} <: AbstractOutlineObject
    p1::Point{T}                # lower left point
    p2::Point{T}                # upper right point
    wlist::Vector{LineWall{T}}

    function OutlineRectangle(p1::Point{T},p2::Point{T}) where T<:Real
        wlist = [LineWall(Point(p2.x,p2.y),Point(p1.x,p2.y)),            # top
                 LineWall(Point(p1.x,p1.y),Point(p2.x,p1.y)),            # bottom
                 LineWall(Point(p1.x,p2.y),Point(p1.x,p1.y)),            # left
                 LineWall(Point(p2.x,p1.y),Point(p2.x,p2.y))]            # right
        return new{T}(p1,p2,wlist)
    end
end

"""
    OutlinePolygon{T}

Struct modeling an outline polygon. The polygon is
composed of an arbitrary number of points.
The base unit is of type T.
"""
struct OutlinePolygon{T} <: AbstractOutlineObject
    pList::Vector{Point{T}}
    p1_bound::Point{T}
    p2_bound::Point{T}
    wlist::Vector{LineWall{T}}

    function OutlinePolygon(pList::Vector{Point{T}}) where T<:Real
        n_objects = length(pList)
        maxX = -10_000.0; minX = -maxX
        maxY = maxX; minY = -maxX
        for ti=1:n_objects
            if pList[ti].x > maxX
                maxX = pList[ti].x
            elseif pList[ti].x < minX
                minX = pList[ti].x
            end
            if pList[ti].y > maxY
                maxY = pList[ti].y
            elseif pList[ti].y < minY
                minY = pList[ti].y
            end
        end
        p1_bound = Point(minX,minY)
        p2_bound = Point(maxX,maxY)
        wlist = Array{LineWall}(undef,n_objects)
        for ti=1:n_objects-1
            wlist[ti] = LineWall(pList[ti],pList[ti+1],zero(T))
        end
        wlist[n_objects] = LineWall(pList[n_objects],pList[1],zero(T))

        return new{T}(pList,p1_bound,p2_bound,wlist)
    end
end

struct FillingCircle{T} <: AbstractFillingObject
    radius::T
    center::Point{T}
    buffer_percent::T
end

"""
    MaterialParameters{T}

Struct that contains the expected volume fraction of the
outline, and the number of filling objects. The base unit
is of type T.
"""
struct MaterialParameters{T}
    expected_volume_fraction::T
    n_objects::Int
end

"""
    PorousStructure{T}

The main product from CircGeometry.jl. Contains a
MaterialParameters object, a list of filling objects, and
supporting array structures. The base unit is of type T.
"""
struct PorousStructure{T}
    param::MaterialParameters{T}
    olist::Vector{FillingCircle{T}}
    xArr::Vector{T}
    yArr::Vector{T}
    radiiArr::Vector{T}

    function PorousStructure(param::MaterialParameters{T},radiiArr::Vector{T},between_buffer::T) where T<:Real
        olist = [FillingCircle(radiiArr[ti],Point(zero(T),zero(T)),between_buffer) for ti=1:param.n_objects]
        xArr = zeros(T,param.n_objects)
        yArr = zeros(T,param.n_objects)
        new{T}(param, olist, xArr, yArr,radiiArr)
    end
end

mutable struct PlacementStatus
    inside::Bool
    intersection_others::Bool
    intersection_walls::Bool
    safe_placement::Bool
    marked_for_shuffling::Bool
    attempt_ind::Int

    function PlacementStatus()
        return new(false,false,false,false,false,1)
    end
end

function reset_placement!(placement::PlacementStatus)
    placement.inside = false
    placement.intersection_others = false
    placement.intersection_walls = false
    placement.safe_placement = false
    placement.marked_for_shuffling = false
    placement.attempt_ind = 1
    nothing
end

function update_attempt!(placement::PlacementStatus)
    placement.attempt_ind += 1
    nothing
end

max_stars_repo_name: pseastham/CircGeometry

max_issues_repo_name: pseastham/CircGeometry


text:
using Dolang

struct RBC end
rbc_ff = let
    # construct args
    variables = [:y, :c, :k, :i, :l, :y_l, :z]
    exogenous = [:e]
    args = [(:y, 0), (:c, 0), (:c, 1), (:k, -1), (:k, 0), (:i, 0), (:l, 0),
            (:l, 1), (:y_l, 0), (:z, -1), (:z, 0), (:z, 1), (:e, 0)]

    # params
    params = [:beta, :psi, :delta, :alpha, :rho]

    # no definitions
    # no targets

    # equations
    eqs = [
        :((1/c) = beta*(1/c(+1))*(1+alpha*(k^(alpha-1))*(exp(z(+1))*l(+1))^(1-alpha)-delta)),
        :(psi*c/(1-l) = (1-alpha)*(k(-1)^alpha)*(exp(z)^(1-alpha))*(l^(-alpha))),
        :(c+i = y),
        :(y = (k(-1)^alpha)*(exp(z)*l)^(1-alpha)),
        :(i = k-(1-delta)*k(-1)),
        :(y_l = y/l),
        :(z = rho*z(-1)+e)
    ]

    map(Dolang._filter_lines!, eqs)

    Dolang.FunctionFactory(RBC, eqs, args, params, funname=:rbc_model)
end;

print_Core.eval(ex) = (println(ex); Core.eval(ex))
print_Core.eval(make_method(rbc_ff))
print_Core.eval(make_method(Der{1}, rbc_ff))
print_Core.eval(make_method(Der{2}, rbc_ff, mutating=false))

# just some reasonable parameter values
alp = 0.33
bet = 0.99
gam = 0.003
mst = 1.011
rho = 0.7
psi = 0.787
del = 0.02
k = 6
P = 2.25
c = 0.45
e = 1
W = 4
R = 1.02
d = 0.85
n = 0.19
l = 0.86
y = 0.6
i = y-c
e = 0
z = e/(1-rho)
V = [y, c, c, k, k, i, l, l, y/l, z, z, z, e]
p = [bet, psi, del, alp, rho]
model = RBC()
levels = rbc_model(model, V, p)
jac = rbc_model(Der{1}, model, V, p)
hes = rbc_model(Der{2}, model, V, p)

max_issues_repo_name: EconForge/Dolang


text:
@enum Sides begin
    LEFT
    RIGHT
end

sign(side::Sides) = side == LEFT ? 1 : -1


text:
@testset "Modified Cholesky decomposition" begin
    @testset "SPD matrix" begin
        dim = 10
        @testset for ntest in 1:15
            a = rand(test_rng, Float64, (dim, dim))
            a = a' * a
            achol = mcholesky!(copy(a))
            # For SPD matrices, must be the same as the usual Cholesky decomposition
            @test isapprox(achol.U' * achol.U, a; rtol = 1e-4)
        end
    end

    @testset "Generic symmetric matrix" begin
        dim = 10
        @testset for ntest in 1:15
            a = rand(test_rng, Float64, (dim, dim))
            a .+= a'
            achol = mcholesky!(copy(a))
            modcholprod = achol.U' * achol.U
            # checking that the decomposition differs from the original matrix
            # only in diagonal elements
            modcholprod .-= Diagonal(modcholprod)
            a .-= Diagonal(a)
            @test isapprox(UpperTriangular(modcholprod), UpperTriangular(a); rtol = 1e-4)
        end
    end
end


text:
using BenchmarkTools
using Printf
using Random
using LinearAlgebra

# Get the number of iterations from the command line.
n, = size(ARGS)
if n < 1
    println("Usage: belief.jl N")
    println("       ---> Please specify the number of iterations.")
    exit()
end

N = parse(Int, ARGS[1])

"""
    beliefpropagation(N)

Runs the belief propagation algorithm N times.
"""
function beliefpropagation(N)
    dim = 5000
    x = ones(dim)
    A = (randn(dim,dim) .+ 1.0)/2.0

    for i = 1:N
        x = log.(A*exp.(x));
        x .-= log.(sum(exp.(x)));
    end
    x
end

println("--------------------------")
println(@sprintf "Belief calculations: %d" N)
println("--------------------------")

@btime beliefpropagation(N)