text:
module IJuliaWidgets

using JSON
using React
using Interact

import Interact.update_view
export mimewritable, writemime

if !isdefined(Main, :IJulia)
    error("IJuliaWidgets must be imported from inside an IJulia notebook")
end

const ijulia_js  = readall(joinpath(dirname(Base.source_path()), "ijulia.js"))

try
    display("text/html", """<script charset="utf-8">$(ijulia_js)</script>""")
catch
end

import IJulia
import IJulia: metadata, display_dict
using  IJulia.CommManager
import IJulia.CommManager: register_comm
import Base: writemime, mimewritable

const comms = Dict{Signal, Comm}()

function get_data_dict(value, mimetypes)
    dict = Dict{ASCIIString, ByteString}()
    for m in mimetypes
        if mimewritable(m, value)
            dict[m] = stringmime(m, value)
        elseif m == "text/latex" && mimewritable("application/x-latex", value)
            dict[string("text/latex")] =
                stringmime("application/x-latex", value)
        else
            warn("IPython seems to be requesting an unavailable mime type")
        end
    end
    return dict
end

function init_comm(x::Signal)
    if !haskey(comms, x)
        subscriptions = Dict{ASCIIString, Int}()
        function handle_subscriptions(msg)
            if haskey(msg.content, "data")
                action = get(msg.content["data"], "action", "")
                if action == "subscribe_mime"
                    mime = msg.content["data"]["mime"]
                    subscriptions[mime] = get(subscriptions, mime, 0) + 1
                elseif action == "unsubscribe_mime"
                    mime = msg.content["data"]["mime"]
                    subscriptions[mime] = get(subscriptions, mime, 1) - 1
                end
            end
        end
        # One Comm channel per signal object
        comm = Comm(:Signal)
        comms[x] = comm   # Backend -> Comm
        # Listen for mime type registrations
        comm.on_msg = handle_subscriptions
        # prevent resending the first time?
        function notify(value)
            mimes = keys(filter((k,v) -> v > 0, subscriptions))
            if length(mimes) > 0
                send_comm(comm, [:value =>
                                 get_data_dict(value, mimes)])
            end
        end
        lift(notify, x)
    else
        comm = comms[x]
    end

    return comm
end

function metadata(x :: Signal)
    comm = init_comm(x)
    return ["reactive"=>true,
            "comm_id"=>comm.id]
end

# Render the value of a signal.
mimewritable(m :: MIME, s :: Signal) =
    mimewritable(m, s.value)

function writemime(io:: IO, m :: MIME, s :: Signal)
    writemime(io, m, s.value)
end

function writemime(io::IO, ::MIME{symbol("text/html")},
          w::InputWidget)
    create_view(w)
end

function writemime(io::IO, ::MIME{symbol("text/html")},
                   w::Widget)
    create_view(w)
end

function writemime{T<:Widget}(io::IO, ::MIME{symbol("text/html")},
                              x::Signal{T})
    create_widget_signal(x)
end

## This is for our own widgets.
function register_comm(comm::Comm{:InputWidget}, msg)
    w_id = msg.content["data"]["widget_id"]
    comm.on_msg = (msg) -> recv(w, msg.content["data"]["value"])
end

JSON.print(io::IO, s::Signal) = JSON.print(io, s.value)

##################### IPython IPEP 23: Backbone.js Widgets #################

## ButtonView ✓
## CheckboxView ✓
## DropdownView ✓
## FloatSliderView ✓
## FloatTextView ✓
## IntSliderView ✓
## IntTextView ✓
## ProgressView
## RadioButtonsView ✓
## SelectView ✓
## TextareaView ✓
## TextView ✓
## ToggleButtonsView ✓
## ToggleButtonView ✓
## AccordionView W
## ContainerView W
## HTMLView W
## ImageView W
## LatexView W
## PopupView W
## TabView W

# Interact -> IJulia view names
view_name(w::Widget) = string(typeof(w).name, "View")
view_name{T<:Integer}(::Slider{T}) = "IntSliderView"
view_name{T<:FloatingPoint}(::Slider{T}) = "FloatSliderView"
view_name{T<:Integer}(::Textbox{T}) = "IntTextView"
view_name{T<:FloatingPoint}(::Textbox{T}) = "FloatTextView"
view_name(::Textbox) = "TextView"
view_name{view}(::Options{view}) = string(view, "View")

function metadata{T <: Widget}(x :: Signal{T})
    Dict()
end

const widget_comms = Dict{Widget, Comm}()
function update_view(w::InputWidget; src::InputWidget=w)
    msg = Dict()
    msg["method"] = "update"
    state = Dict()
    state["msg_throttle"] = 3
    state["_view_name"] = view_name(src)
    state["description"] = w.label
    state["visible"] = true
    state["disabled"] = false
    state["readout"] = true
    msg["state"] = merge(state, statedict(src))
    send_comm(widget_comms[w], msg)
end

function update_view(w::Widget; src::Widget=w)
    msg = Dict()
    msg["method"] = "update"
    state = Dict()
    state["msg_throttle"] = 3
    state["_view_name"] = view_name(src)
    state["description"] = w.label
    state["visible"] = true
    state["disabled"] = false
    msg["state"] = merge(state, statedict(src))
    send_comm(widget_comms[w], msg)
    nothing
end

function create_view(w::Widget)
    if haskey(widget_comms, w)
        comm = widget_comms[w]
    else
        comm = Comm(:WidgetModel)
        widget_comms[w] = comm
        # Send a full state update message.
        update_view(w)

        # dispatch messages to widget's handler
        comm.on_msg = msg -> handle_msg(w, msg)
        nothing # display() nothing
    end

    send_comm(comm, ["method"=>"display"])
end

function create_widget_signal(s)
    create_view(s.value)
    local target = s.value
    lift(x->update_view(target, src=x), s, init=nothing)
end

include("statedict.jl")
include("handle_msg.jl")

end


text:
module TestEnsembleMethods

using Test
using Random
using AutoMLPipeline.AbsTypes
using AutoMLPipeline.Utils
using AutoMLPipeline.EnsembleMethods
using AutoMLPipeline.DecisionTreeLearners
using DataFrames

function generateXY()
    Random.seed!(123)
    iris = getiris()
    indx = Random.shuffle(1:nrow(iris))
    features=iris[indx,1:4]
    sp = iris[indx,5] |> Vector
    (features,sp)
end

function getprediction(model,features,output)
  res = fit_transform!(model,features,output)
  sum(res .== output)/length(output)*100
end

function test_ensembles()
  tstfeatures,tstoutput = generateXY()
  models = [VoteEnsemble(),StackEnsemble(),BestLearner()]
  for model in models
    @test getprediction(model,tstfeatures,tstoutput) > 90.0
  end
end
@testset "Ensemble learners" begin
  Random.seed!(123)
  test_ensembles()
end

function test_vararg()
  rf = RandomForest()
  ada = Adaboost()
  pt = PrunedTree()

  X,Y = generateXY()
  vote = VoteEnsemble(rf,ada,pt)
  stack = StackEnsemble(rf,ada,pt)
  best = BestLearner(rf,ada,pt)
  v=fit_transform!(vote,X,Y)
  s=fit_transform!(stack,X,Y)
  p=fit_transform!(best,X,Y)
  @test score(:accuracy,v,Y) > 90.0
  @test score(:accuracy,s,Y) > 90.0
  @test score(:accuracy,p,Y) > 90.0
end # module
@testset "Vararg Ensemble Test" begin
  Random.seed!(123)
  test_vararg()
end

end

ext: jl


text:
"""
    overlap(skt::SlaterKosterTable)

Return a function equivalent to `(args...) -> overlap(skt, args...)`.

"""
overlap(skt::SlaterKosterTable) = (args...) -> overlap(skt, args...)

"""
    overlap(skt::SlaterKosterTable, r::NTuple{3,<:Real}, ϕ1, ϕ2)

Return the overlap matrix element between `ϕ2` and `ϕ1`, where `ϕ1` is
located at `r` relative to `ϕ2`, given in Angstrom.

`ϕ1`, `ϕ2` are pairs of `element_symbol => orbital_symbol`, where
`element_symbol` is a valid `element_symbol` from `elementsymbols(skt)`,
and `orbital_symbol` are tesseral spherical harmonic orbitals.

Accepted symbols for the orbitals are:

:s, :px, :py, :pz

!!! note
    The following orbitals are not yet supported: :dx2y2, :dxy, :dxz,
    :dyz, :dz2, :f...

"""
function overlap(skt::SlaterKosterTable, r, ϕ1, ϕ2)
    orbital_symbol1, orbital_symbol2 = last(ϕ1), last(ϕ2)
    if _orbitalrank(orbital_symbol2) < _orbitalrank(orbital_symbol1)
        ϕ1, ϕ2 = ϕ2, ϕ1
        r = -r
    end
    table = skt.data[first(ϕ1), first(ϕ2)]
    return overlap(table, r, last(ϕ1), last(ϕ2))
end

"""
    overlap(skt::SlaterKosterTable, ϕ)

Return the diagonal overlap matrix element corresponding to `ϕ`.
`ϕ` is a pair of `element_symbol => orbital_symbol`, where
`element_symbol` is a valid `element_symbol` from
`elementsymbols(skt)`, and `orbital_symbol` is a tesseral spherical
harmonic orbital.

Accepted symbols for the orbitals are:

:s, :px, :py, :pz

!!! note
    The following orbitals are not yet supported: :dx2y2, :dxy, :dxz,
    :dyz, :dz2, :f...

"""
overlap(skt::SlaterKosterTable, ϕ) = overlap(skt.data[first(ϕ), first(ϕ)], last(ϕ))

overlap(skt::HomoNuclearTable, orbital_symbol::Symbol) = 1

function overlap(skt::PrimitiveTable, r, orbital_symbol1::Symbol, orbital_symbol2::Symbol)
    # Convert distance to atomic units for indexing
    r = r/a0
    @assert length(r) == 3 "`length(r) = $(length(r)) != 3`: `r` is not a valid 3-dimensional vector: $r"
    min_dist, max_dist = extrema(skt.integral_table[!, :dist])
    norm_r = norm(r)
    @assert min_dist <= norm_r <= max_dist "`norm(r) = $norm_r` is outside table bounds: [$min_dist, $max_dist]"
    grid_dist = skt.grid_dist
    lo = floor(Int, norm_r/grid_dist)
    hi = norm_r == max_dist ? lo : lo + 1
    # To support automatic differentiation, we cannot use the
    # remainder operator (`%`) here, but subtract instead
    rem = norm_r/grid_dist - lo
    # Interpolate the table for the given distance
    slice = max(1, lo-1):min(skt.n_grid_points, hi+1)
    iprev, ilo = lo > 1 ? (1, 2) : (1, 1)
    ihi, inext = hi < skt.n_grid_points ? (ilo + 1, ilo + 2) : (ilo + 1, ilo + 1)
    table = DataFrame(Dict(
        name => spline_catmullrom(rem, col[iprev], col[ilo], col[ihi], col[inext]) for (name,col) in collect(pairs(eachcol(skt.integral_table[slice, :])))
    ))
    # Normalize the separation vector
    d = r ./ norm_r
    return overlap(table, d, Val(orbital_symbol1), Val(orbital_symbol2))[]
end

function overlap(df::DataFrame, d, ::Val{:s}, ::Val{:s})
    return df.Sss0
end

for (i, x) in [(1, :(:x)), (2, :(:y)), (3, :(:z))]
    @eval @inline function overlap(df::DataFrame, d, ::Val{:s}, ::Val{Symbol("p", $x)})
        return d[$i]*df.Ssp0
    end

    for (j, y) in [(1, :(:x)), (2, :(:y)), (3, :(:z))]
        if i == j
            @eval @inline function overlap(df::DataFrame, d, ::Val{Symbol("p", $x)}, ::Val{Symbol("p", $x)})
                return d[$i]^2*df.Spp0 + (1 - d[$i]^2)*df.Spp1
            end
        else
            @eval @inline function overlap(df::DataFrame, d, ::Val{Symbol("p", $x)}, ::Val{Symbol("p", $y)})
                return d[$i]*d[$j]*(df.Spp0 - df.Spp1)
            end
        end
    end
end


text:
module TestDicts

using Test
using DataFramesMeta

y = 3
d = Dict(:s => 3, :y => 44, :d => 5, :e => :(a + b))
@test @with(d, :s + :y) == d[:s] + d[:y]
@test @with(d, :s + y)  == d[:s] + y
@test @with(d, d)  == d
@test @with(d, :s + d[^(:y)])  == d[:s] + d[:y]
@test @with(d, :e.head) == d[:e].head
@test @with(Dict(:s => 3, :y => 44, :d => 5, :e => :(a + b)), :e.head) == d[:e].head

x = @with d begin
    z = y + :y - 1
    :s + z
end
@test x == y + d[:y] - 1 + d[:s]

fun = d -> @with d begin
    z = y + :y - 1
    :s + z
end
@test fun(d) == y + d[:y] - 1 + d[:s]

d2 = @transform(d, z = :y + :s)
@test d2[:z] == 47

d2 = @select(d, :y, z = :y + :s, :e)
@test d2[:z] == 47

@test DataFramesMeta.with_helper(:df, :(f.(1))) == :(f.(1))
@test DataFramesMeta.with_helper(:df, :(f.(b + c))) == :(f.(b + c))
@test DataFramesMeta.with_helper(:df, Expr(:., :a, QuoteNode(:b))) == Expr(:., :a, QuoteNode(:b))
@test DataFramesMeta.select_helper(:df, 1).args[2].args[2].args[2].args[3] == 1
@test DataFramesMeta.select_helper(:df, QuoteNode(:a)).args[2].args[2].args[2].args[2].args[2].args[2].args[3].args[1] == :a
@test DataFramesMeta.expandargs(QuoteNode(:a)) == Expr(:kw, :a, QuoteNode(:a))
@test DataFramesMeta.byrow_find_newcols(:(;), Any[]) == (Any[], Any[])

end # module


text:
module {{{PKG}}}

# Write your package code here.

end


text:
mutable struct RepeatedMeasuresXMLElement <: ModelExtensionXMLElement
    el::XMLOrNothing
    precision::MatrixParameter
    precision_prior::WishartPriorXMLElement
    treeModel_el::TreeModelXMLElement
    tree_trait_ind::Int
    trait_name::String
    tp_name::String
    standardize_traits::Bool
    id::String

    function RepeatedMeasuresXMLElement(treeModel_el::TreeModelXMLElement;
                                        trait_ind::Int = 1)

        p = treeModel_el.trait_dims[trait_ind]
        trait_name = treeModel_el.node_traits[trait_ind]
        tp_name = treeModel_el.param_names[trait_ind]
        precision = MatrixParameter(Diagonal(ones(p)), bn.DEFAULT_RM_PREC_NAME)
        wishart_prior = WishartPriorXMLElement(precision)
        return new(nothing, precision, wishart_prior,
                   treeModel_el,
                   trait_ind, trait_name, tp_name, true, bn.DEFAULT_RM_NAME)
    end
end

function name(::RepeatedMeasuresXMLElement)
    return bn.REPEATED_MEASURES
end

function make_xmlelement(rm::ResidualVarianceModel, tm::TreeModelXMLElement;
                         ind::Int = 1)
    rmxml = RepeatedMeasuresXMLElement(tm, trait_ind = ind)
    set_precision(rmxml, inv(rm.Γ))
    return rmxml
end

function make_xml(rm_el::RepeatedMeasuresXMLElement)

    # repeatedMeasures xml element
    el = new_element(bn.REPEATED_MEASURES)
    attrs = [(bn.ID, rm_el.id), (bn.TRAIT_NAME, rm_el.trait_name)]
    if rm_el.standardize_traits
        set_attribute(el, bn.STANDARDIZE, bn.TRUE)
    end
    set_attributes(el, attrs)
    treeModel_el = make_xml(rm_el.treeModel_el)
    add_ref_el(el, treeModel_el)
    tp_el = new_child(el, bn.TRAIT_PARAMETER)
    add_ref_el(tp_el, bn.PARAMETER, rm_el.tp_name)
    sp_el = new_child(el, bn.SAMPLING_PRECISION)
    mp_el = make_xml(rm_el.precision)
    add_child(sp_el, mp_el)

    # make_xml(rm_el.mbd_el)

    # add_ref_el(el, rm_el.mbd_el.el)
    rm_el.el = el

    # prior xml element

    make_xml(rm_el.precision_prior)

    return el
end

function get_loggables(el::RepeatedMeasuresXMLElement)
    make_xml(el)
    return [el.el[bn.SAMPLING_PRECISION][1][bn.MATRIX_PARAMETER][1]]
end

function get_matrix_parameter(el::RepeatedMeasuresXMLElement)
    return get_loggables(el)[1]
end

function get_priors(xml::RepeatedMeasuresXMLElement)
    make_xml(xml)
    return [xml.precision_prior.el]
end

function get_precision_prior(xml::RepeatedMeasuresXMLElement)
    return xml.precision_prior.el
end

function set_precision(rm::RepeatedMeasuresXMLElement,
                       mat::AbstractArray{Float64, 2})
    check_posdef(mat)
    set_mat!(rm.precision, mat)
end


text:
reload("Rotolo")

module Try
end

module Try
using Rotolo

compile(joinpath(dirname(@__FILE__), "../test/dummy-page.jl"),
        "c:/temp/dummy.pdf")

compile(joinpath(dirname(@__FILE__), "../examples/Himmelblau.jl"),
        "c:/temp/Himmelblau.pdf")

jlsrc = normpath(joinpath(dirname(@__FILE__), "../examples/slides.jl"))
compile(jlsrc, "c:/temp/slides.pdf")

Rotolo.sessions
Rotolo.currentSession

phantomjspath = joinpath(Pkg.dir("PhantomJS"), "deps/usr/bin/phantomjs")

function renderhtml(source::IO;
                   format::String="png",
                   width::Int=1024,
                   height::Int=800,
                   clipToSelector::String="",
                   quality::Int=75,
                   paperSize::String="A4",
                   orientation::String="portrait",
                   margin::Union{String, Int}=0,
                   background::String="white")

 local result

 mktempdir() do tdir
   htmlpath = joinpath(tdir, randstring() * ".html")
   open(io -> write(io, read(source)), htmlpath, "w")
   htmlpath = replace(htmlpath, "\\", "/")

   destpath = joinpath(tdir, randstring() * ".png")
   destpath = replace(destpath, "\\", "/")

   bgjs = background == "transparent" ? "" :
            """page.evaluate(function() {
                 document.body.bgColor = '$background';
               });
            """

   clipjs = clipToSelector == "" ? "" :
              """
               var clipRect = page.evaluate(function(){
                   return document.querySelector('$clipToSelector').getBoundingClientRect();
                 });

               page.clipRect = {
                 top:    clipRect.top,
                 left:   clipRect.left,
                 width:  clipRect.width,
                 height: clipRect.height
               };
              """

   jsscript = """
     "use strict";
     var page = require('webpage').create(),
         system = require('system'),
         address, output, size, pageWidth, pageHeight;

     address = 'file:///$htmlpath';
     output = '$destpath';

     page.viewportSize = { width: $width, height: $height };

     page.paperSize = { format: '$paperSize',
                        orientation: '$orientation',
                        margin: '$margin' };

     page.onConsoleMessage = function(msg, lineNum, sourceId) {
       console.log('CONSOLE: ' + msg);
     };

     page.onError = function(msg, trace) {
       console.log('ERROR : ' + msg);
     };

     page.open(address, function (status) {
         if (status !== 'success') {
             console.log('Unable to load the address : ' + address);
             phantom.exit(1);
         } else {
             console.log('reading : ' + address);
             window.setTimeout(function () {
                 console.log('now rendering ');
                 $bgjs
                 $clipjs
                 page.render(output,
                             {format: '$format',
                              quality: '$quality'});
                 console.log('rendering done');
                 phantom.exit();
             }, 2000);
         }
     });
   """

   jspath = joinpath(tdir, randstring() * ".js")
   open(io -> write(io, jsscript), jspath, "w")

   run(`$phantomjspath $jspath`)
   result = read(destpath)
 end

 result
end

fn = "C:/Users/frtestar/AppData/Local/Temp/jl_6DE1.tmp.html"
fn = "C:/temp/test.html"
fn = "C:/temp/test2.html"
sio = open(fn, "r")
pdfbuf = renderhtml(sio, format="pdf")
open(pdfio -> write(pdfio, pdfbuf), "c:/temp/test3.pdf", "w")
close(sio)

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

reload("Rotolo")

module Try
end

module Try
using Rotolo

using Base.Markdown

@redirect Float32 Base.Markdown.MD

@session d1

Rotolo.send(Rotolo.currentSession, 0, "load",
            Dict(:assetname => "katex",
                 :assetpath => "D:/frtestar/.julia/v0.5/Rotolo/client/katex/katex.js"))

@redirect Rotolo.Katex

Katex("c = \\pm\\sqrt{a^2 + b^2}", true)

Float32(45)
"abcd" |> Rotolo.style("background-color:lightpink")

style("abcd", "background-color:lightpink")
"abcd" |> *("xyz")

theme1 = "display:flex;flex-wrap:wrap;width:300px"

@container flex2 style => "display:flex;flex-wrap:wrap;width:300px"

style("abcd", "padding:10px;")

cs = Rotolo.currentSession
function xplore(ct, level=0)
  println(" "^level, ct.name, " ($(ct.nid))")
  for c in ct.subcontainers
    xplore(c, level+2)
  end
end

xplore(cs.root_container)

@container abcd
xplore(cs.root_container)

@container cyz
xplore(cs.root_container)

@container
xplore(cs.root_container)

Rotolo.@container abcd.yo
xplore(cs.root_container)

md"test"
Float32(123)

Rotolo.@container abcd.yourk style=>"font-color:red"
xplore(cs.root_container)

md"test"

using DataFrames
@redirect DataFrame
DataFrames.head(Main.comp)

Rotolo.unfold(:(abcd.yo))

macro test(arg)
  dump(arg)
  println(Rotolo.unfold(arg))
end

@test abcd
@test abcd.youi.yo456

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

Rotolo._send(Rotolo.currentSession, 0,
             "load",
             Dict{Symbol,Any}(:assetname => "katex",
                              # :assetpath => "D:/frtestar/.julia/v0.5/Rotolo/client/katex/katex.js",
                              :assetpath => "/home/fred/.julia/v0.5/Rotolo/client/katex/katex.js")
            )

Rotolo.send("append",
     Dict(:newnid => Rotolo.getnid(),
          :compname => "katex",
          :params => Dict(:expr => "x^2+y_x=a")))

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

Rotolo._send(Rotolo.currentSession, 0,
             "load",
             Dict{Symbol,Any}(:assetname => "vegalite",
            #  :assetpath => "D:/frtestar/.julia/v0.5/Rotolo/client/vegalite/vegalite.js")
             :assetpath => "/home/fred/.julia/v0.5/Rotolo/client/VegaLite/vegalite.js")
            )

Rotolo.send("append",
     Dict(:newnid => Rotolo.getnid(),
          :compname => "vegalite",
          :params => Dict(:test => "abcd")))

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

Rotolo._send(Rotolo.currentSession, 0,
             "load",
             Dict{Symbol,Any}(:assetname => "katex",
                              # :assetpath => "D:/frtestar/.julia/v0.5/Rotolo/client/katex/katex.js",
                              :assetpath => "/home/fred/.julia/v0.5/Rotolo/client/katex/katex.js")
            )

Rotolo.send("append",
     Dict(:newnid => Rotolo.getnid(),
          :compname => "katex",
          :params => Dict(:expr => "x^2+y_x=a")))

Rotolo._send(Rotolo.currentSession, 0,
             "load", Dict{Symbol,Any}(:assetname => "vegalite",
                          :assetpath => "../vegalite.js"))

D:/frtestar/.julia/v0.5/Rotolo/client/katex/ka
D:\frtestar\devl\paper-client\dist

Rotolo.send("append",
     Dict(:newnid => Rotolo.getnid(),
          :compname => "vegalite",
          :params => Dict(:expr => "x^2+y_x=a")))

dump(:( ~test ))

456.4
Float32(456.54)

# import Media, Atom
# methods(Media.render, (Atom.Editor, Any))

md"""
  markdown Text
  # titleqsdqsd

  qsdfdf
  - qsdfgqdg
  - dfgqdrg

  ```
  code_llvm
  ```

  text *text*  qsdf  **qdfgqf**
  """

using Base.Markdown

str = md"test **test** test"
typeof(str)
methodswith(Markdown.MD)

type Abcd
  x::String
end

Abcd("test")

show(io::IO, ::MIME"text/plain",x::Abcd) = show(io, "Abcd : $(x.x)")

function show(io::IO, x::Abcd)
  show(io, "Abcd : $(x.x)")
end

show(Abcd("testtest"))

end

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

source = IOBuffer("ancd")

mktempdir() do tdir
  println(tdir)
  htmlpath = joinpath(tdir, randstring() * ".html")
  open(io -> write(io, read(source)), htmlpath, "w")
  println(read(htmlpath))
end

  # source = IOBuffer(read("c:/temp/vegalite.html"))

reload("PhantomJS")

module Try
end

module Try

using PhantomJS
using VegaLite

ts = sort(rand(10))
ys = Float64[ rand()*0.1 + cos(x) for x in ts]

v = data_values(time=ts, res=ys) +    # add the data vectors & assign to symbols 'time' and 'res'
      mark_line() +                   # mark type = line
      encoding_x_quant(:time) +       # bind x dimension to :time, quantitative scale
      encoding_y_quant(:res) ; nothing

pio = IOBuffer()
VegaLite.writehtml(pio, v)
position(pio)
seekstart(pio)

pngdata = renderhtml(pio)

open(io -> write(io, pngdata), "c:/temp/result.png", "w")

tmppath = "c:/temp/vegalite.html"
open(io -> VegaLite.writehtml(io, v), tmppath, "w")

exepath = joinpath(Pkg.dir("Wkhtmltox"), "deps/src/bin")
run(`$exepath/wkhtmltopdf c:/temp/example.html c:/temp/example.pdf`)

run(`$exepath/wkhtmltopdf -H`)

pdf_init(1)
ps = PdfSettings("out" => "c:/temp/example.pdf")

os = PdfObject("page" => "c:/temp/example.html")
os["web.printMediaType"] = "true"
os["web.printMediaType"]

conv = PdfConverter(ps, os) # pass the pdf settings and source settings
run(conv)

conv = nothing
pdf_deinit()

end


text:
#!/usr/bin/env julia

using Luxor

CAIRO_PATH_MOVE_TO = 0
CAIRO_PATH_LINE_TO = 1
CAIRO_PATH_CURVE_TO = 2
CAIRO_PATH_CLOSE_PATH = 3

using Test

using Random
Random.seed!(42)

function get_path(str)
    sethue("blue")
    setopacity(.7)
    setline(0.5)
    p = textpath(str)
    o = getpath()
    strokepath()
    sethue("red")
    x, y = 0, 0
    for e in o
        if e.element_type == CAIRO_PATH_MOVE_TO
            (x, y) = e.points
            move(x, y)
        elseif e.element_type == CAIRO_PATH_LINE_TO
            (x, y) = e.points
            # straight lines
            line(x, y)
            strokepath()
            circle(x, y, 1, :stroke)
        elseif e.element_type == CAIRO_PATH_CURVE_TO
            (x1, y1, x2, y2, x3, y3) = e.points
            # show bezier control lines
            circle(x1, y1, 1, :stroke)
            circle(x2, y2, 1, :stroke)
            circle(x3, y3, 1, :stroke)
            move(x, y)
            # Use point interface, to make sure that is also tested
            curve(Point(x1, y1), Point(x2, y2), Point(x3, y3))
            strokepath()
            (x, y) = (x3, y3) # update current point
        elseif e.element_type == CAIRO_PATH_CLOSE_PATH
            closepath()
        else
            error("unknown CairoPathEntry " * repr(e.element_type))
            error("unknown CairoPathEntry " * repr(e.points))
        end
    end
end

function get_path_test(fname)
    Drawing(500, 500, fname)
    origin()
    @layer begin
        fontsize(24)
        fontface("Times-Italic")
        setline(0.5)
        translate(-100, -100)
        for theta in range(0, length=20, stop=2π)
            @layer begin
                rotate(theta)
                translate(30, 0)
                sethue("red")
                textoutlines("Luxor", O, :fill, halign=:left, valign=:baseline)
                sethue("blue")
                textoutlines("Luxor", O, :stroke, halign=:left, valign=:baseline)
            end
        end
    end
    fontsize(80)
    get_path("Luxor")
    sethue("magenta")
    @test finish() == true
    println("...finished test: output in $(fname)")
end
get_path_test("get-path-test.pdf")


text:
__precompile__(true)
module MacroTools

using Compat

export @match, @capture

include("match.jl")
include("types.jl")
include("union.jl")
include("macro.jl")
include("utils.jl")

include("examples/destruct.jl")

end # module


text:
module Evolutionary

using Random

export Strategy, strategy, inverse, mutationwrapper,
# ES mutations
isotropic, anisotropic, isotropicSigma, anisotropicSigma,
# GA mutations
flip, domainrange, inversion, insertion, swap2, scramble, shifting,
# ES recombinations
average, marriage, averageSigma1, averageSigmaN,
# GA recombinations
singlepoint, twopoint, uniform,
discrete, waverage, intermediate, line,
pmx, ox1, cx, ox2, pos,
# GA selections
ranklinear, uniformranking, roulette, sus, tournament, #truncation
# Optimization methods
es, cmaes, ga

const Strategy = Dict{Symbol,Any}
const Individual = Union{Vector, Matrix, Function, Nothing}

# Wrapping function for strategy
function strategy(; kwargs...)
    result = Dict{Symbol,Any}()
    for (k, v) in kwargs
        result[k] = v
    end
    return result
end

# Inverse function for reversing optimization direction
function inverseFunc(f::Function)
    function fitnessFunc(x::T) where {T <: Vector}
        return 1.0/(f(x)+eps())
    end
    return fitnessFunc
end

# Obtain individual
function getIndividual(init::Individual, N::Int)
    if isa(init, Vector)
        @assert length(init) == N "Dimensionality of initial population must be $(N)"
        individual = init
    elseif isa(init, Matrix)
        @assert size(init, 1) == N "Dimensionality of initial population must be $(N)"
        populationSize = size(init, 2)
        individual = init[:, 1]
    elseif isa(init, Function) # Creation function
        individual = init(N)
    else
        individual = ones(N)
    end
    return  individual
end

# Collecting interim values
function keep(interim, v, vv, col)
    if interim
        if !haskey(col, v)
            col[v] = typeof(vv)[]
        end
        push!(col[v], vv)
    end
end

# Print a bar
function bar(n,total)
    percent = n/total*100.
    size = percent/2
    text = ""
    print("\r Progress:|")
    for a = 0:size
        text *="█"
    end
    printstyled(text, color=:red)
    text = ""
    for a = 1:(50.0 - size)
        text *="█"
    end
    printstyled(text, color=:white)
    print("| $(round(percent))%")
end

# ES & GA recombination functions
include("recombinations.jl")

# ES & GA mutation functions
include("mutations.jl")

# GA selection functions
include("selections.jl")

# Evolution Strategy
include("es.jl")
include("cmaes.jl")

# Genetic Algorithms
include("ga.jl")

end


text:
@testset "Hermite 3D" begin

    f(x, y, z) = cos(x^2 + y^2 + z^2)

    xmin, xmax, nx = -1, 1, 32
    ymin, ymax, ny = -1, 1, 32
    zmin, zmax, nz = -1, 1, 32

    ϵ = 0.1
    γ = 3.0

    hermite_x = Hermite(Chebyshev( xmin, xmax, nx ), ϵ, γ)
    hermite_y = Hermite(Chebyshev( ymin, ymax, ny ), ϵ, γ)
    hermite_z = Hermite(Chebyshev( zmin, zmax, nz ), ϵ, γ)

    xk = hermite_x.nodes
    yk = hermite_y.nodes
    zk = hermite_z.nodes

    nxe, nye, nze = 32, 64, 128

    xe = collect(range(xmin, stop=xmax, length=nxe))
    ye = collect(range(ymin, stop=ymax, length=nye))
    ze = collect(range(zmin, stop=zmax, length=nze))
    dx = abs(xe[2] - xe[1])
    dy = abs(ye[2] - ye[1])
    dz = abs(ze[2] - ze[1])

    fk = [f(x, y, z) for x in xk, y in yk, z in zk]
    fe = [f(x, y, z) for x in xe, y in ye, z in ze]

    @time s  = interpolate(hermite_x, hermite_y, hermite_z, fk, xe, ye, ze)

    max_error = maximum(abs.(s .- fe))

    l2_error  = sqrt(trapz((s - fe).^2, dx, dy, dz))

    println(" Hermite 3d max error : $max_error ")
    println(" Hermite 3d L2  error : $l2_error  ")

    @test true

end


text:
using Documenter
using WaveProp

makedocs(
    sitename = "WaveProp",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    modules = [WaveProp],
    pages = [
        "Home"      => ["index.md",]
        # "Tutorials" => ["Tutorials/soundsoft_scattering.md",]
        "Modules"   => ["Utils/Utils.md",
                        "Geometry/Geometry.md",
                        "Interpolation/Interpolation.md",
                        "Integration/Integration.md",
                        "Mesh/Mesh.md",
                        "PhysicalProblem/PhysicalProblem.md",
                        "FEM/FEM.md",
                        "Nystrom/Nystrom.md",
                        "IO/IO.md"
                        ]
        # "Examples" => "examples.md"
        # "Documentation" => "devdocs.md"
    ]
)

# Deploy to gh-pages. See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/WaveProp/WaveProp.jl.git",
    devbranch = "main"
)

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
using BinaryProvider

const platform = platform_key()
const prefix = Prefix(!isempty(ARGS) ? ARGS[1] : joinpath(@__DIR__,"usr"))

libiconv = LibraryProduct(prefix, "libiconv")

const bin_prefix = "https://github.com/Keno/IConvBuilder/releases/download/1.15"

# TODO Update hash values, only the Windows x64 hash is correct right now
const download_info = Dict(
    Linux(:aarch64, :glibc) => ("$bin_prefix/libiconv.aarch64-linux-gnu.tar.gz", "c4abb8e04b424faac414bb3abdf2bc29eed2e63a83fb5acbc99a4a548a2dc643"),
    Linux(:powerpc64le, :glibc) => ("$bin_prefix/libiconv.powerpc64le-linux-gnu.tar.gz", "e5b0474b0e77c16db21427e73452acddc08f80d6fc59656fafe383814afa4578"),
    Windows(:i686) => ("$bin_prefix/libiconv.i686-w64-mingw32.tar.gz", "0b7820ccef0d0f227fa2e9e592af87fe9cd032827f207849d015739ec2d2f0a3"),
    Linux(:x86_64, :glibc) => ("$bin_prefix/libiconv.x86_64-linux-gnu.tar.gz", "3a2f8df4ebe057d22554853df6e506de29a53ab0d2dc2f9ae96e733b529177c5"),
    MacOS() => ("$bin_prefix/libiconv.x86_64-apple-darwin14.tar.gz", "1c46e48169663bc138d0adcbee95877d0361bd8c57b8f30c3d0e0a81200d3dd1"),
    Windows(:x86_64) => ("$bin_prefix/libiconv.x86_64-w64-mingw32.tar.gz", "cf447074e801d3fb5209285580a90dd0bf7fcaaf512af3baeb6eb3ac1f5f6664"),
    Linux(:armv7l, :glibc) => ("$bin_prefix/libiconv.arm-linux-gnueabihf.tar.gz", "d855e816d781cb69eacae52f508df09b4deba7143e3c950b9ea386023f81645b"),
    Linux(:i686, :glibc) => ("$bin_prefix/libiconv.i686-linux-gnu.tar.gz", "3d9309dd256509136428072ee6b8542c40c59c5c7e11ac58a6284b4a40aac593"),
)

if platform in keys(download_info)
    # Grab the url and tarball hash for this particular platform
    url, tarball_hash = download_info[platform]

    install(url, tarball_hash; prefix=prefix, force=true, verbose=true)

    # Finaly, write out a deps file containing paths to libiconv
    @write_deps_file libiconv
else
    error("Your platform $(Sys.MACHINE) is not recognized, we cannot install libiconv.")
end


text:
module V3P

# import Base.+, Base.*, Base./, Base.norm, Base.==, Base.<, Base.>
import Base: +, *, /, ==, <, >, norm, dot
export Vec3, norm, dist

immutable Vec3{T}
    x::T
    y::T
    z::T
end

(+)(a::Vec3, b::Vec3) = Vec3(a.x+b.x, a.y+b.y, a.z+b.z)
(*)(p::Vec3, s::Real) = Vec3(p.x*s, p.y*s, p.z*s)
(*)(s::Real, p::Vec3) = p*s
(/)(p::Vec3, s::Real) = (1.0/s)*p

(==)(a::Vec3, b::Vec3) = (a.x==b.x)&&(a.y==b.y)&&(a.z==b.z) ? true : false;

dot(a::Vec3, b::Vec3) = a.x*b.x + a.y*b.y + a.z*b.z;
norm(a::Vec3) = sqrt(dot(a,a));

(<)(a::Vec3, b::Vec3) = norm(a) < norm(b) ? true : false;
(>)(a::Vec3, b::Vec3) = norm(a) > norm(b) ? true : false;

dist(a::Vec3, b::Vec3) = sqrt((a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y) + (a.z-b.z)*(a.z-b.z))

end


text:
"""
**RMSProp Optimizer**
```julia
RMSprop(; η::Real=0.001, γ::Real=0.01, ϵ::Real=1e-8)
```

Algorithm

```math
\\begin{align*}
    E[g^2]_t =& \\gamma E[g^2]_{t-1}+(1-\\gamma) g_t^2\\\\
    \\Delta x_t =& \\frac{\\eta}{\\sqrt{E[g^2]_t + \\epsilon}}g_t
\\end{align*}
```
    [Reference](http://www.cs.toronto.edu/~tijmen/csc321/slides/lecture_slides_lec6.pdf)
"""
mutable struct RMSprop <: Optimizer
    opt_type::String
    t::Int64
    ϵ::Float64
    η::Float64
    γ::Float64
    E_g²_t::AbstractArray
end

function RMSprop(; η::Real=0.001, γ::Real=0.01, ϵ::Real=1e-8)
    @assert η > 0.0 "η must be greater than 0"
    @assert γ > 0.0 "γ must be greater than 0"
    @assert ϵ > 0.0 "ϵ must be greater than 0"

    RMSprop("RMSprop", 0, ϵ, η, γ, [])
end

params(opt::RMSprop) = "ϵ=$(opt.ϵ), η=$(opt.η), γ=$(opt.γ)"

function update(opt::RMSprop, g_t::AbstractArray{T}) where {T<:Real}
    # resize accumulated and squared updates
    if opt.t == 0
        opt.E_g²_t = zero(g_t)
    end

    # accumulate gradient
    opt.E_g²_t = opt.γ * opt.E_g²_t + (one(T) - opt.γ) * (g_t .^ 2)

    # compute update
    RMS_g_t = sqrt.(opt.E_g²_t .+ opt.ϵ)
    δ = opt.η * g_t ./ RMS_g_t

    return δ
end


text:
Δx = 0.02
xlim = (-5.98,5.98)
ylim = (-5.98,5.98)
g = PhysicalGrid(xlim,ylim,Δx)

w = Nodes(Dual,size(g))
wc = Nodes(Dual,size(g),dtype=ComplexF64)

dq = Edges(Dual,w)
oc = similar(w)
oc .= 1

n = 150
radius = 1.0
body = Circle(radius,n)
X = VectorData(collect(body))
f = VectorData(X)
ϕ = ScalarData(f)
fill!(ϕ,1.0)

regop = Regularize(X,cellsize(g),I0=origin(g),weights=dlengthmid(body),ddftype=CartesianGrids.Yang3)
Hv = RegularizationMatrix(regop,f,dq)
Hs = RegularizationMatrix(regop,ϕ,w)

@testset "Double Layer" begin

  dlayer = DoubleLayer(body,Hv,weight=1.0)
  ϕ .= rand(length(ϕ))

  @test abs(dot(oc,dlayer(ϕ),g)) < 100*eps(1.0)

  dlayer2 = DoubleLayer(body,g,w)

  @test dlayer2(ϕ) == dlayer(ϕ)

  v = Edges(Primal,w)
  dvlayer = DoubleLayer(body,g,v)

  f .= rand(length(f))
  dvlayer(f)

  dlayer2 = deepcopy(dlayer)
  DoubleLayer!(dlayer2,body,g)

  @test dlayer.H == dlayer2.H

  dlayerc = DoubleLayer(body,g,wc)

  @test eltype(dlayerc(1)) == ComplexF64

end

@testset "Single Layer" begin

  slayer = SingleLayer(body,Hs) #,weight=cellsize(g)^2)
  ϕ .= 1.0

  @test abs(dot(oc,slayer(ϕ),g)-2π*radius) < 2e-3

  slayer2 = SingleLayer(body,g,w)

  @test slayer2(ϕ) == slayer(ϕ)

  slayer2 = deepcopy(slayer)
  SingleLayer!(slayer2,body,g)

  @test slayer.H == slayer2.H

end

@testset "Masks" begin

  inner = Mask(body,g,w)

  fill!(w,1)

  @test abs(dot(oc,inner(w),g) - π*radius^2) < 1e-3

  out = similar(w)
  inner(out,w)
  @test out == inner(w)

  outer = ComplementaryMask(inner)
  outer(out,w)
  @test out == outer(w)

  innerc = Mask(body,g,wc)

  @test eltype(innerc(wc)) == ComplexF64

end


text:
"""
    wrapout(x)

Wrap Symata output for printing. Most methods are defined
in Output.jl.
"""
wrapout(x) = x

## TODO: disable symwarn, so we are forced to use @symwarn below
#for f in ( :print, :println, :string, :warn, :error, :throw)
for f in ( :print, :println, :string, :error, :throw)
   @eval function $(Symbol("sym" * string(f)))(args...)
        $(f)(map(wrapout,args)...)
    end
end

## This might take some work. We need to do wrapout on the evaluated arguments,
## which gets into the guts of @warn.
## FIXME. Probably need to write our own version of log_msg :(
"""
    @symwarn(arg)

Performs macro call `@warn` on `wrapout(arg)` so that
Julia displays objects of type `Mxpr` correctly.
"""
macro symwarn(args...)
#    exprs =  map(x -> :(wrapout($(esc(x)))), args)
    exprs =  map(x -> :(wrapout($(esc(x)))), args)
    :(@warn($(exprs...)))
#    :(@warn(wrapout($(esc(arg)))))
#    :(@warn($(wrapout(esc(wrapout(args))))))
#    :(@warn($(map(wrapout, args)...)))
#    :(@warn(map(x -> esc(wrapout(x), $args))))
end

@doc """
   symprintln(expr)

prints `expr` as a Symata expression. For instance, symbols
are printed without the leading `:`. Rational numbers are printed
like `x/y` rather than `x//y`.
""" symprintln


text:
# Autogenerated wrapper script for Clingcon_jll for x86_64-unknown-freebsd
export clingcon, libclingcon

using Clingo_jll
JLLWrappers.@generate_wrapper_header("Clingcon")
JLLWrappers.@declare_executable_product(clingcon)
JLLWrappers.@declare_library_product(libclingcon, "libclingcon.so")
function __init__()
    JLLWrappers.@generate_init_header(Clingo_jll)
    JLLWrappers.@init_executable_product(
        clingcon,
        "bin/clingcon",
    )

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

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Macro Utilities.

macro extend(x, ex)
       fields = (:($v::$t) for (v, t) in zip(fieldnames(eval(x)), fieldtypes(eval(x))))
       append!(ex.args[3].args, fields)
       esc(ex)
end

macro selfdoc(ex)

    @require ex isa Expr && ex.head in (:function, Symbol("="))

    docex = copy(ex)
    Base.remove_linenums!(docex)
    if ex.head == Symbol("=")
        body = docex.args[2]
        if body.head == :block && length(body.args) == 1
            docex.args[2] = body.args[1]
        end
    end
    call = ex.args[1]
    if call.head != :call
        call = call.args[1]
    end
    name = call.args[1]
    doc = string("```\n", docex, "\n```\n")
    esc(quote
        @doc($doc, $name)
        Base.@__doc__($ex)
    end)
end

# FIXME copy/paste from base/util.jl: https://git.io/J8vC5
"""
    @invoke f(arg::T, ...; kwargs...)
Provides a convenient way to call [`invoke`](@ref);
`@invoke f(arg1::T1, arg2::T2; kwargs...)` will be expanded into `invoke(f, Tuple{T1,T2}, arg1, arg2; kwargs...)`.
When an argument's type annotation is omitted, it's specified as `Any` argument, e.g.
`@invoke f(arg1::T, arg2)` will be expanded into `invoke(f, Tuple{T,Any}, arg1, arg2)`.
"""
macro invoke(ex)
    f, args, kwargs = destructure_callex(ex)
    arg2typs = map(args) do x
        Base.is_expr(x, :(::)) ? (x.args...,) : (x, GlobalRef(Core, :Any))
    end
    args, argtypes = first.(arg2typs), last.(arg2typs)
    return esc(:($(GlobalRef(Core, :invoke))($(f), Tuple{$(argtypes...)}, $(args...); $(kwargs...))))
end

function destructure_callex(ex)
    Base.is_expr(ex, :call) || throw(ArgumentError("a call expression f(args...; kwargs...) should be given"))

    f = first(ex.args)
    args = []
    kwargs = []
    for x in ex.args[2:end]
        if Base.is_expr(x, :parameters)
            append!(kwargs, x.args)
        elseif Base.is_expr(x, :kw)
            push!(kwargs, x)
        else
            push!(args, x)
        end
    end

    return f, args, kwargs
end

# End of file: macroutils.jl


text:
abstract type CallbackRunner end

struct LinearRunner <: CallbackRunner end

function handle(runner::LinearRunner, event::Event, phase, learner)
    idxs = Zygote.ignore() do
        topological_sort_by_dfs(learner.callbacks.graph)
    end
    for i in idxs
        _on(event, phase, learner.callbacks.cbs[i], learner)
    end
end


text:
@testset "Rewriting" begin

    Al = Alphabet(["a", "e", "b", "p"])
    KnuthBendix.set_inversion!(Al, "a", "e")
    KnuthBendix.set_inversion!(Al, "b", "p")
    lenlexord = LenLex(Al)

    a, A, b, B = [Word([i]) for i in 1:4]
    ε = one(a)

    @testset "Free Rewriting" begin
        @test KnuthBendix.rewrite_from_left(a * A, Al) == ε
        @test KnuthBendix.rewrite_from_left(a * A * b, Al) == b
        @test KnuthBendix.rewrite_from_left(a * B * b, Al) == a
        @test KnuthBendix.rewrite_from_left(a * B * b * A, Al) == ε
        @test KnuthBendix.rewrite_from_left(a * B * b * a, Al) == a * a
        @test KnuthBendix.rewrite_from_left(a * B * b * a * A, Al) == a
        @test KnuthBendix.rewrite_from_left(a * b * A * B, Al) == a * b * A * B
    end

    @testset "Rule simplification" begin
        Al = deepcopy(Al)
        push!(Al, "z")

        prefix = Word(rand(1:length(Al)-1, 100)) # all invertible
        suffix = Word(rand(1:length(Al)-1, 100)) # all invertible
        l = Word([5, 1, 2, 2])
        r = Word([5, 1, 2, 1])

        @test KnuthBendix.simplifyrule!(prefix * l, prefix * r, Al) == (l, r)
        @test KnuthBendix.simplifyrule!(l * suffix, r * suffix, Al) == (l, r)
        @test KnuthBendix.simplifyrule!(prefix * l * suffix, prefix * r * suffix, Al) ==
              (l, r)
        @test KnuthBendix.simplifyrule!(l * suffix, prefix * r * Word([5]), Al) ==
              (l * suffix, prefix * r * Word([5]))
        @test KnuthBendix.simplifyrule!(copy(l), copy(r), Al) == (l, r)

        @test KnuthBendix.simplifyrule!(a^4, a^2, Al) == (a^2, one(a))
    end

    @testset "Rewriting System operations" begin

        R = RewritingSystem(
            [a * A => ε, A * a => ε, b * B => ε, B * b => ε, a * b => b * a],
            lenlexord,
            bare = true,
        )

        Z = empty(R)

        @test R isa KnuthBendix.AbstractRewritingSystem
        @test R isa RewritingSystem

        @test R !== Z
        @test isempty(Z)
        @test !isempty(R)

        push!(Z, KnuthBendix.Rule(b * B => ε))
        @test collect(KnuthBendix.rules(Z)) == [KnuthBendix.Rule(b * B => ε)]
        @test KnuthBendix.ordering(Z) == lenlexord

        KnuthBendix.deactivate!(first(KnuthBendix.rules(Z)))
        @test isempty(collect(KnuthBendix.rules(Z)))

        @test KnuthBendix.rewrite_from_left(a * A, R) == ε
        @test KnuthBendix.rewrite_from_left(b * B, Z) == b * B

        KnuthBendix.deactivate!(first(KnuthBendix.rules(R)))
        @test KnuthBendix.rewrite_from_left(a * A, R) == a * A

        push!(Z, KnuthBendix.Rule(b * B => ε))
        @test isempty(KnuthBendix.rules(empty!(Z)))
        @test isempty(KnuthBendix.rules(empty(R)))
        @test KnuthBendix.rewrite_from_left(b * B, Z) == b * B

        @test sprint(show, Z) isa String
    end
end


text:
using Base.Cartesian

function cudims(n::Integer)
  threads = min(n, 256)
  ceil(Int, n / threads), threads
end

cudims(a::AbstractArray) = cudims(length(a))

@inline ind2sub_(a::AbstractArray{T,0}, i) where T = ()
@inline ind2sub_(a, i) = Tuple(CartesianIndices(a)[i])

macro cuindex(A)
  quote
    A = $(esc(A))
    i = (blockIdx().x-1) * blockDim().x + threadIdx().x
    i > length(A) && return
    ind2sub_(A, i)
  end
end

@generated function nindex(i::T, ls::NTuple{N,T}) where {N,T}
  na = one(i)
  quote
    Base.@_inline_meta
    $(foldr((n, els) -> :(i ≤ ls[$n] ? ($n, i) : (i -= ls[$n]; $els)), :($na, $na), one(i):i(N)))
  end
end

@inline function catindex(dim, I::NTuple{N}, shapes) where N
  @inbounds x, i = nindex(I[dim], getindex.(shapes, dim))
  x, ntuple(n -> n == dim ? i : I[n], Val{N})
end

function growdims(dim, x)
  if ndims(x) >= dim
    x
  else
    reshape(x, size.((x,), 1:dim)...)
  end
end

function _cat(dim, dest, xs...)
  function kernel(dim, dest, xs)
    I = @cuindex dest
    @inbounds n, I′ = catindex(dim, Int.(I), size.(xs))
    @inbounds dest[I...] = xs[n][I′...]
    return
  end
  xs = growdims.(dim, xs)
  blk, thr = cudims(dest)
  @cuda blocks=blk threads=thr kernel(dim, dest, xs)
  return dest
end

function Base.cat_t(dims::Integer, T::Type, x::CuArray, xs::CuArray...)
  catdims = Base.dims2cat(dims)
  shape = Base.cat_shape(catdims, (), size.((x, xs...))...)
  dest = Base.cat_similar(x, T, shape)
  _cat(dims, dest, x, xs...)
end

Base.vcat(xs::CuArray...) = cat(xs..., dims=1)
Base.hcat(xs::CuArray...) = cat(xs..., dims=2)

"""
    @sync ex

Run expression `ex` and synchronize the GPU afterwards. This is a CPU-friendly
synchronization, i.e. it performs a blocking synchronization without increasing CPU load. As
such, this operation is preferred over implicit synchronization (e.g. when performing a
memory copy) for high-performance applications.

It is also useful for timing code that executes asynchronously.
"""
macro sync(ex)
    quote
        local e = CuEvent(CUDAdrv.EVENT_BLOCKING_SYNC | CUDAdrv.EVENT_DISABLE_TIMING)
        local ret = $(esc(ex))
        CUDAdrv.record(e)
        CUDAdrv.synchronize(e)
        ret
    end
end

"""
    @argout call(..., out(output_arg), ...)

Change the behavior of a function call, returning the value of an argument instead.
A common use case is to pass a newly-created Ref and immediately dereference that output:

    @argout(some_getter(Ref{Int}()))[]

If no output argument is specified, `nothing` will be returned.
Multiple output arguments return a tuple.

"""
macro argout(ex)
    Meta.isexpr(ex, :call) || throw(ArgumentError("@argout macro should be applied to a function call"))

    block = quote end

    # look for output arguments (`out(...)`)
    output_vars = []
    args = ex.args[2:end]
    for (i,arg) in enumerate(args)
        if Meta.isexpr(arg, :call) && arg.args[1] == :out
            # allocate a variable
            @gensym output_val
            push!(block.args, :($output_val = $(ex.args[i+1].args[2]))) # strip `output(...)`
            push!(output_vars, output_val)

            # replace the argument
            ex.args[i+1] = output_val
        end
    end

    # generate a return
    push!(block.args, ex)
    if isempty(output_vars)
        push!(block.args, :(nothing))
    elseif length(output_vars) == 1
        push!(block.args, :($(output_vars[1])))
    else
        push!(block.args, :(tuple($(output_vars...))))
    end

    esc(block)
end

"""
    @workspace size=getWorkspaceSize(args...) [eltyp=UInt8] [fallback=nothing] buffer->begin
      useWorkspace(workspace, sizeof(workspace))
    end

Create a GPU workspace vector with element type `eltyp` and size in number of elements (in
the default case of an UInt8 element type this equals to the amount of bytes) determined by
calling `getWorkspaceSize`, and pass it to the closure for use in calling `useWorkspace`.
A fallback workspace size `fallback` can be specified if the regular size would lead to OOM.
Afterwards, the buffer is put back into the memory pool for reuse.

This helper protects against the rare but real issue of `getWorkspaceSize` returning
different results based on the GPU device memory pressure, which might change _after_
initial allocation of the workspace (which can cause a GC collection).

Use of this macro should be as physically close as possible to the function that actually
uses the workspace, to minimize the risk of GC interventions between the allocation and use
of the workspace.

"""
macro workspace(ex...)
    code = ex[end]
    kwargs = ex[1:end-1]

    sz = nothing
    eltyp = :UInt8
    fallback = nothing
    for kwarg in kwargs
        key,val = kwarg.args
        if key == :size
            sz = val
        elseif key == :eltyp
            eltyp = val
        elseif key == :fallback
            fallback = val
        else
            throw(ArgumentError("Unsupported keyword argument '$key'"))
        end
    end

    if sz === nothing
        throw(ArgumentError("@workspace macro needs a size argument"))
    end

    # break down the closure to a let block to prevent JuliaLang/julia#15276
    Meta.isexpr(code, :(->)) || throw(ArgumentError("@workspace macro should be applied to a closure"))
    length(code.args) == 2 || throw(ArgumentError("@workspace closure should take exactly one argument"))
    code_arg = code.args[1]
    code = code.args[2]

    return quote
        sz = $(esc(sz))
        workspace = nothing
        try
          while workspace === nothing || length(workspace) < sz
              workspace = CuArray{$(esc(eltyp))}(undef, sz)
              sz = $(esc(sz))
          end
        catch ex
            $fallback === nothing && rethrow()
            isa(ex, OutOfGPUMemoryError) || rethrow()
            workspace = CuArray{UInt8}(undef, $fallback)
        end

        let $(esc(code_arg)) = workspace
            $(esc(code))
        end
    end
end

if VERSION <= v"1.1"
  # JuliaLang/julia#30187-like functionality, but only for CuContext dicts to avoid clashes
  function Base.get!(default::Base.Callable, d::IdDict{CuContext,V}, @nospecialize(key)) where {V}
    if !haskey(d, key)
      d[key] = default()
    end
    return d[key]
  end
end


text:
using Pkg; Pkg.activate(@__DIR__)
using Comrade
using Distributions
using ComradeGalactic
using ComradeAHMC
using GalacticBBO
using Plots
using StatsBase
using GalacticOptimJL
using DistributionsAD

# load eht-imaging we use this to load eht data
load_ehtim()
# To download the data visit https://doi.org/10.25739/g85n-f134
obs = ehtim.obsdata.load_uvfits(joinpath(@__DIR__, "SR1_M87_2017_096_hi_hops_netcal_StokesI.uvfits"))
obs.add_scans()
# kill 0-baselines since we don't care about
# large scale flux and make scan-average data
obs = obs.flag_uvdist(uv_min=0.1e9).avg_coherent(0.0, scan_avg=true).add_fractional_noise(0.02)
# extract log closure amplitudes and closure phases
dlcamp = extract_lcamp(obs)
dcphase = extract_cphase(obs)
lklhd = RadioLikelihood(dlcamp, dcphase)

# Build the Model. Here we we a struct to hold some caches
# which will speed up imaging
struct ImModel{C}
    cache::C
    fov::Float64
    npix::Int
end

# For our model we will be using a rasterized image. This can be viewed as something like a
# non-parametric model. As a result of this we will need to use a `modelimage` object to
# store cache information we will need to compute the numerical FT.
function (model::ImModel)(θ)
    (;c) = θ
    #Construct the image model
    cmat = reshape(c, model.npix, model.npix)
    img = IntensityMap(cmat, model.fov, model.fov, BSplinePulse{3}())
    #Create the modelimage object that will use a cache to compute the DFT
    m = modelimage(img, cache)
end
const npix = 10
const fovxy = μas2rad(65.0)
# Now we can feed in the array information to form the cache. We will be using a DFT since
# it is efficient for so few pixels
cache = create_cache(Comrade.DFTAlg(dlcamp), IntensityMap(rand(npix,npix), fovxy, fovxy, BSplinePulse{3}()))
mms = ImModel(cache, fovxy, npix)
# We will use a Dirichlet prior to enforce that the flux sums to unity since closures are
# degenerate to total flux.
prior = (c = DistributionsAD.TuringDirichlet(npix*npix, 1.0),)

post = Posterior(lklhd, prior, mms)
tpost = asflat(post)

# Let's run an optimizer to get a nice starting location
# It turns out that gradients are really helpful here
ndim = dimension(tpost)
f = OptimizationFunction(tpost, GalacticOptim.AutoZygote())
x0 = randn(ndim)
prob = OptimizationProblem(f, x0, nothing)
sol = solve(prob, LBFGS(); maxiters=4_000, show_trace=true, show_every=50, g_tol=1e-2)

xopt = transform(tpost, sol)

# Let's see how the fit looks
residual(mms(xopt), dlcamp)
residual(mms(xopt), dcphase)
plot(mms(xopt), fovx=1.2*fovxy, fovy=1.2*fovxy)

# now we sample using hmc
metric = DenseEuclideanMetric(ndim)
hchain, stats = sample(post, AHMC(;metric, autodiff=AD.ZygoteBackend()), 3000; nadapts=2000, init_params=xopt)

foox = let dc=dc, dd=dd, pr=tpost.lpost.prior, tr=tpost.transform, plan = cache.plan, phases=cache.phases, dmat=dlcamp.config.designmat, dmatc=dcphase.config.designmat
    x->begin
        y = transform(tr, x)
        lp = logdensityof(pr, y)
        vis = conj.(plan*y.c).*phases
        logdensityof(dd, dmat*log.(abs.(vis))) + logdensityof(dc, dmatc*angle.(vis)) + lp
    end
end


text:
module Wrangling

using Compat: Compat, contains

export contains
export startswith_any, endswith_any, contains_any

include("symbol_searching.jl")

end # module


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

mutable struct CustomDomain <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    useSubDomainName::Any # spec type: Union{ Nothing, Bool } # spec name: useSubDomainName

    function CustomDomain(;name=nothing, useSubDomainName=nothing)
        o = new()
        validate_property(CustomDomain, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(CustomDomain, Symbol("useSubDomainName"), useSubDomainName)
        setfield!(o, Symbol("useSubDomainName"), useSubDomainName)
        o
    end
end # type CustomDomain

const _property_map_CustomDomain = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("useSubDomainName")=>Symbol("useSubDomainName"))
const _property_types_CustomDomain = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("useSubDomainName")=>"Bool")
Base.propertynames(::Type{ CustomDomain }) = collect(keys(_property_map_CustomDomain))
Swagger.property_type(::Type{ CustomDomain }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_CustomDomain[name]))}
Swagger.field_name(::Type{ CustomDomain }, property_name::Symbol) =  _property_map_CustomDomain[property_name]

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

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


text:
using Documenter, AutonomousMerging

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

deploydocs(
    repo = "github.com/sisl/AutonomousMerging.jl.git",
)


text:
# nicer functions for accessing the hashes from julia
module KDF

export scrypt, bcrypt, pbkdf2

end


text:
using Distributed
using Logging

# Add four workers and give them all function definitions
addprocs(4)
@everywhere using AdaptiveResonance

"""
    artscene_filter_porcelain()

Runs the artscene user-level functions on a random image.
"""
# function artscene_filter_porcelain()
@testset "ARTSCENE Filter Porcelain" begin
    # Set the logging level to Debug
    LogLevel(Logging.Debug)
    n_processes = nprocs()
    n_workers = nworkers()
    @info "Processes: $n_processes, Workers: $n_workers"

    # Random image
    raw_image = rand(3, 5, 5)

    # Process the image through the filters
    O, C = artscene_filter(raw_image)

    # Set the logging level back to Info
    LogLevel(Logging.Info)
end # @testset "ARTSCENE Filter Porcelain"

# Close the workers after testing
rmprocs(workers())


text:
using SymPy

x, y, t = @syms x y t
u1 = (x^2+y^2)*exp(-t)*0.1
u2 = (x^2-y^2)*exp(-t)*0.1
ε = [
    diff(u1, x)
    diff(u2, y)
    diff(u2, x) + diff(u1, y)
]

ν = 0.35
E = 2.0

H = zeros(3,3)
H[1,1] = E*(1. -ν)/((1+ν)*(1. -2. *ν));
H[1,2] = H[1,1]*ν/(1-ν);
H[2,1] = H[1,2];
H[2,2] = H[1,1];
H[3,3] = H[1,1]*0.5*(1. -2. *ν)/(1. -ν);

σ = H * ε
f1 = u1 - ( diff(σ[1],x) + diff(σ[3], y) )
f2 = u2 - ( diff(σ[2],y) + diff(σ[3], x) )

g1 = -σ[3]
g2 = -σ[2]

function code(s)
    s = sympy.julia_code(s)
    replace(replace(s, ".*"=>"*"), ".^"=>"^")
end

println(code(f1))
println(code(f2))
println(code(g1))
println(code(g2))


text:
using Documenter
using Lathe: models, stats, preprocess
using Lathe
Documenter.makedocs(root = "./",
       source = "src",
       build = "build",
       clean = true,
       doctest = true,
       modules = Module[Lathe,models,stats,preprocess],
       repo = "https://github.com/emmettgb/Lathe.jl/",
       highlightsig = true,
       sitename = "Lathe.jl",
       expandfirst = [],
       pages = [
               "Index" => "index.md",
               "Stats" => "stats.md",
               "Preprocess" => "preprocess.md",
               "Models" => "models.md"
               ]
       )


text:
using TrajectoryOptimization
import TrajectoryOptimization: stage_cost, cost_expansion
using TrajOptPlots
const TO = TrajectoryOptimization

# using PlanningWithAttitude
# import PlanningWithAttitude: cost_gradient, stage_cost, cost_hessian
# include("visualization.jl")
using StaticArrays
using LinearAlgebra
using BenchmarkTools
using Plots
using ForwardDiff
using Random
using MeshCat

# Dynamics
model = Dynamics.Satellite()
N = 101
tf = 5.0

if !isdefined(Main,:vis)
    vis = Visualizer(); open(vis);
    set_mesh!(vis, model)
end

# Objective
struct SatCost <: CostFunction
    Q::Diagonal{Float64,SVector{3,Float64}}
    R::Diagonal{Float64,SVector{3,Float64}}
    q::SVector{3,Float64}
    q0::SVector{4,Float64}
    Qq::Float64
    c::Float64
end

function SatCost(Q::Diagonal,R::Diagonal,q0::SVector, Qq=0.0, ω0=@SVector zeros(3))
    q = -Q*ω0
    c = 0.5*ω0'Q*ω0
    SatCost(Q,R,q,q0,Qq,c)
end

function stage_cost(cost::SatCost, x::SVector, u::SVector=(@SVector zeros(3)))
    ω = @SVector [x[1], x[2], x[3]]
    q = @SVector [x[4], x[5], x[6], x[7]]
    J = 0.5*(ω'cost.Q*ω + u'cost.R*u) + cost.q'ω + cost.c
    d = cost.q0'q
    if d ≥ 0
        J += cost.Qq*(1-d)
    else
        J += cost.Qq*(1+d)
    end
    return J
end

function cost_expansion(cost::SatCost, model::Dynamics.Satellite, z::KnotPoint{T,N,M}, G) where {T,N,M}
    Q = cost.Q
    q0 = cost.q0
    u = control(z)
    ω = @SVector [z.z[1], z.z[2], z.z[3]]
    q = @SVector [z.z[4], z.z[5], z.z[6], z.z[7]]
    d = cost.q0'q
    Qω = Q*ω + cost.q
    if d ≥ 0
        Qq = -cost.Qq*cost.q0'Lmult(q)*Vmat()'
    else
        Qq =  cost.Qq*cost.q0'Lmult(q)*Vmat()'
    end
    Qx = [Qω; Qq']
    Qu = cost.R*u

    Qxx = Diagonal(@SVector [Q[1,1], Q[2,2], Q[3,3], 0,0,0])
    d = q0'q*cost.Qq
    if d ≤ 0
        d *= -1
    end
    Qqq = I(3)*d
    Qxx = @SMatrix [
        Q[1,1] 0 0 0 0 0;
        0 Q[2,2] 0 0 0 0;
        0 0 Q[3,3] 0 0 0;
        0 0 0 Qqq[1] Qqq[4] Qqq[7];
        0 0 0 Qqq[2] Qqq[5] Qqq[8];
        0 0 0 Qqq[3] Qqq[6] Qqq[9];
    ]
    Quu = cost.R
    Qux = @SMatrix zeros(M,N-1)
    return Qxx, Quu, Qux, Qx, Qu
end

Q_diag = @SVector [1e-3,1e-3,1e-3, 1e-2,1e-2,1e-2,1e-2]
R_diag = @SVector fill(1e-4,3)
θs = range(0,2pi,length=N)
u = normalize([1,0,0])
iω = @SVector [1,2,3]
iq = @SVector [4,5,6,7]

Xref = map(1:N) do k
    θ = θs[k] #+ deg2rad(185)
    ωf = [2pi/5,0,0]
    qf = @SVector [cos(θ/2), u[1]*sin(θ/2), u[2]*sin(θ/2), u[3]*sin(θ/2)]
    xf = @SVector [ωf[1], ωf[2], ωf[3], qf[1], qf[2], qf[3], qf[4]]
end

costs = map(1:N) do k
    ωf = Xref[k][iω]
    qf = Xref[k][iq]
    k < N ? s = 1 : s = 1
    # LQRCost(Q_diag*s, R_diag, xf)
    SatCost(Diagonal(Q_diag[iω]), Diagonal(R_diag), qf, 1.0, ωf)
    # xf = Xref[k]
    # QuatLQRCost(Diagonal(Q_diag), Diagonal(R_diag), xf)
end
obj = Objective(costs)

# Initial Condition
x0 = @SVector [0,0,0, 1,0,0,0.]
u0 = @SVector [0.3, 0, 0]
Random.seed!(1)
U0 = [u0 + randn(3)*4e-1 for k = 1:N-1]

# Final Condition
θ = θs[end]
xf = @SVector [0,0,0, cos(θ/2), u[1]*sin(θ/2), u[2]*sin(θ/2), u[3]*sin(θ/2)]

# Solver
prob = Problem(model, obj, xf, tf, x0=x0)
solver = iLQRSolver(prob)
initial_controls!(solver, U0)
rollout!(solver)
visualize!(vis, model, solver.Z)
cost(obj, solver.Z)

solver.opts.verbose = true
initial_controls!(solver, U0)

cost(solver)
TO.state_diff_jacobian!(solver.G, solver.model, solver.Z)
TO.discrete_jacobian!(solver.∇F, solver.model, solver.Z)
TO.cost_expansion!(solver.Q, solver.G, solver.obj, solver.model, solver.Z)
ΔV = TO.backwardpass!(solver)
TO.forwardpass!(solver, ΔV, cost(solver))
for k = 1:N
    solver.Z[k].z = solver.Z̄[k].z
end

initial_controls!(solver, U0)
solve!(solver)
solver.stats.iterations
states(solver)[end][4:7]'xf[4:7]
visualize!(vis, model, solver.Z)
Z = solver.Z

solver.opts.verbose = false
@btime begin
    initial_controls!($solver, $U0)
    solve!($solver)
end

############################################################################################
#                                 Intermediate Keyframes
############################################################################################

Q_nom = @SVector [1e-3,1e-3,1e-3, 1e-4,1e-4,1e-4,1e-4]
R_nom = @SVector fill(1e-4,3)

Q_key = @SVector [1e-3,1e-3,1e-3, 1e-2,1e-2,1e-2,1e-2]
R_key = @SVector fill(1e-4,3)

keyframes = Int.(range(1,N, length=3))
xf = @SVector [0,0,0, cos(θ/2), u[1]*sin(θ/2), u[2]*sin(θ/2), u[3]*sin(θ/2)]

costs = map(1:N) do k
    if k ∈ keyframes
        ωf = Xref[k][iω]
        qf = Xref[k][iq]
        SatCost(Diagonal(Q_key[iω]), Diagonal(R_key), qf, 1.0, ωf)
    else
        ωf = xf[iω]
        qf = xf[iq]
        SatCost(Diagonal(Q_nom[iω]), Diagonal(R_nom), qf, 0.0, ωf)
    end
end
obj = Objective(costs)

# Initial Condition
x0 = @SVector [0,0,0, 1,0,0,0.]
u0 = @SVector [0.3, 0, 0]
Random.seed!(1)
U0 = [u0 + randn(3)*4e-1 for k = 1:N-1]

# Solver
prob = Problem(model, obj, xf, tf, x0=x0)
solver = iLQRSolver(prob)
initial_controls!(solver, U0)
rollout!(solver)
visualize!(vis, model, solver.Z)
cost(obj, solver.Z)

initial_controls!(solver, U0)
# initial_controls!(solver, controls(Z_sol))
solver.opts.verbose = true
solve!(solver)
solver.stats.iterations
states(solver)[end][4:7]'xf[4:7]
visualize!(vis, model, solver.Z)
cost(solver)
Z = solver.Z

Zsol = deepcopy(Z)

iq = @SVector [4,5,6,7]
n =  13
function to_diag(iq,::Val{N}) where N
    D = @MVector zeros(N)
end
@btime to_diag($iq, Val($n))
@btime @MVector zeros(3)


text:
# Funciones de apoyo para generar funciones de inflación MAI

# MAI-F
InflationCoreMaiF(n::Int) = InflationCoreMai(MaiF(n))
InflationCoreMaiF(q::Vector{T}) where {T<:AbstractFloat} = InflationCoreMai(MaiF(T[0, q..., 1]))

# MAI-G
InflationCoreMaiG(n::Int) = InflationCoreMai(MaiG(n))
InflationCoreMaiG(q::Vector{T}) where {T<:AbstractFloat} = InflationCoreMai(MaiG(T[0, q..., 1]))

# MAI-FP
InflationCoreMaiFP(n::Int) = InflationCoreMai(MaiFP(n))
InflationCoreMaiFP(q::Vector{T}) where {T<:AbstractFloat} = InflationCoreMai(MaiFP(T[0, q..., 1]))


text:
#permutation routine function-
#general logic: receive array of permutation parameters, until a model more likely than the least is found:
#randomly select a model from the ensemble (the least likely having been removed by this point), then sample new models by permuting with each of hte given parameter sets until a model more likely than the current contour is found
#if none is found for the candidate model, move on to another candidate until the models_to_permute iterate is reached, after which return nothing for an error code

#four permutation modes: source (iterative random changes to sources until model lh>contour or iterate limit reached)
#						-(iterates, weight shift freq per source base, length change freq per source, weight_shift_dist (a ContinuousUnivariateDistribution)) for permute params
#						mix (iterative random changes to mix matrix as above)
#						-(iterates, unitrange of # of moves)
#						init (iteratively reinitialize sources from priors)
#						-(iterates) for init params
#						merge (iteratively copy a source + mix matrix row from another model in the ensemble until lh>contour or iterate						limit reached)
#						-(iterates) for merpge params

function run_permutation_routine(e::Bayes_IPM_ensemble, job_sets::Vector{Tuple{Vector{Tuple{String,Any}},Vector{AbstractFloat}}}, job_set_thresh::Vector{AbstractFloat}, job_limit::Integer, models_to_permute::Integer, contour::AbstractFloat)
    start=time()
    job_set,job_weights=job_sets[findlast(thresh->(contour>thresh),job_set_thresh)]
    !(length(job_weights)==length(job_set)) && throw(ArgumentError("Job set and job weight vec must be same length!"))

	for model = 1:models_to_permute
		m_record = rand(e.models)
        m = deserialize(m_record.path)

        for job in 1:job_limit
            mode, params = get_job(job_set,job_weights,m.flags,m.sources,m.source_length_limits)
            if mode == "PS"
                new_m=permute_source(m, contour, e.obs_array, e.obs_lengths, e.bg_scores, e.source_priors, params...)
            elseif mode == "PM"
                new_m=permute_mix(m, contour, e.obs_array, e.obs_lengths, e.bg_scores, e.source_priors, params...)
            elseif mode == "PSFM"
				new_m=perm_src_fit_mix(m, contour, e.obs_array, e.obs_lengths, e.bg_scores, e.source_priors, params...)
            elseif mode == "FM"
                new_m=fit_mix(m, e.obs_array, e.obs_lengths, e.bg_scores)
            elseif mode == "DM"
                new_m=distance_merge(e.models, m, contour, e.obs_array, e.obs_lengths, e.bg_scores, e.source_priors, params...)
            elseif mode == "SM"
                new_m=similarity_merge(e.models, m, contour, e.obs_array, e.obs_lengths, e.bg_scores, e.source_priors, params...)
			elseif mode == "RD"
				new_m=random_decorrelate(m, contour, e.obs_array, e.obs_lengths, e.bg_scores, e.source_priors, params...)
			elseif mode == "RI"
				new_m=reinit_src(m, contour,  e.obs_array, e.obs_lengths, e.bg_scores, e.source_priors, e.mix_prior, params...)
            elseif mode == "EM"
                new_m=erode_model(m, contour, e.obs_array, e.obs_lengths, e.bg_scores, e.source_priors, params...)
                occursin("PSFM",new_m.flags[1]) && (mode="PSFM")
            else
                @error "Malformed permute mode code! Current supported: \"PS\", \"PM\", \"PSFM\", \"FM\", \"DM\", \"SM\",\"RD\", \"RI\", \"EM\""
            end

			dupecheck(new_m,m) && new_m.log_Li > contour && return new_m, (time()-start, job, model, m.log_Li, new_m.log_Li, mode)
		end
	end
	return nothing, nothing
end

function worker_permute(e::Bayes_IPM_ensemble, job_chan::RemoteChannel, models_chan::RemoteChannel, job_sets::Vector{Tuple{Vector{Tuple{String,Any}},Vector{AbstractFloat}}}, job_set_thresh::Vector{AbstractFloat}, job_limit::Integer, models_to_permute::Integer)
	persist=true
    id=myid()
    model_ctr=1
	while persist
        wait(job_chan)

        start=time()

		models = fetch(job_chan)
        models === nothing && (persist=false) && break
		contour, ll_idx = findmin([model.log_Li for model in models])
		deleteat!(models, ll_idx)
        job_set,job_weights=job_sets[findlast(thresh->(contour>thresh),job_set_thresh)]
        !(length(job_weights)==length(job_set)) && throw(ArgumentError,"Job set and job weight vec must be same length!")

		for model=1:models_to_permute
			found::Bool=false
			m_record = rand(models)
			m = remotecall_fetch(deserialize,1,m_record.path)
            for job in 1:job_limit
                mode, params = get_job(job_set,job_weights,m.flags,m.sources,m.source_length_limits)
                if mode == "PS"
                    new_m=permute_source(m, contour, e.obs_array, e.obs_lengths, e.bg_scores, e.source_priors, params...)
                elseif mode == "PM"
                    new_m=permute_mix(m, contour, e.obs_array, e.obs_lengths, e.bg_scores, e.source_priors, params...)
                elseif mode == "PSFM"
                    new_m=perm_src_fit_mix(m, contour, e.obs_array, e.obs_lengths, e.bg_scores, e.source_priors, params...)
                elseif mode == "FM"
                    new_m=fit_mix(m, e.obs_array, e.obs_lengths, e.bg_scores)
                elseif mode == "DM"
                    new_m=distance_merge(e.models, m, contour, e.obs_array, e.obs_lengths, e.bg_scores, e.source_priors, params...)
                elseif mode == "SM"
                    new_m=similarity_merge(e.models, m, contour, e.obs_array, e.obs_lengths, e.bg_scores, e.source_priors, params...)
                elseif mode == "RD"
                    new_m=random_decorrelate(m, contour, e.obs_array, e.obs_lengths, e.bg_scores, e.source_priors, params...)
                elseif mode == "RI"
                    new_m=reinit_src(m, contour,  e.obs_array, e.obs_lengths, e.bg_scores, e.source_priors, e.mix_prior, params...)
                elseif mode == "EM"
                    new_m=erode_model(m, contour, e.obs_array, e.obs_lengths, e.bg_scores, e.source_priors, params...)
                    occursin("PSFM",new_m.flags[1]) && (mode="PSFM")
                else
                    @error "Malformed permute mode code! Current supported: \"PS\", \"PM\", \"PSFM\", \"FM\", \"DM\", \"SM\",\"RD\", \"RI\", \"EM\""
                end
				dupecheck(new_m,m) && new_m.log_Li > contour && (put!(models_chan, (new_m ,id, (time()-start, job, model_ctr, m.log_Li, new_m.log_Li, mode))); found=true; model_ctr=1; break)
			end
            found==true && break;
            model_ctr+=1
            wait(job_chan)
            fetch(job_chan)!=models && (break) #if the ensemble has changed during the search, update it
			model==models_to_permute && (put!(models_chan,nothing);persist=false)#worker to put nothing on channel if it fails to find a model more likely than contour
		end
	end
end

                function get_job(job_set, job_weights, flags, srcs, src_lims)
                    job_dist=Categorical(job_weights)
                    mode,params=job_set[rand(job_dist)]

                    if ("nofit" in flags)
                        while mode=="FM"
                            mode, params =job_set[rand(job_dist)]
                        end
                    end

                    if mode=="EM" && !any([>(size(srcs[s][1],1), src_lims[1]) for s in 1:length(srcs)])
                        while mode=="EM"
                            mode,params=job_set[rand(job_dist)]
                        end
                    end
                    return mode,params
                end

                function dupecheck(new_model, model)
                    (new_model.sources==model.sources && new_model.mix_matrix==model.mix_matrix) ? (return false) : (return true)
                end


text:
## generic `Tables.rows` and `Tables.columns` fallbacks
## if a table provides Tables.rows or Tables.columns,
## we'll provide a default implementation of the dual

# generic row iteration of columns
rowcount(cols) = length(getproperty(cols, propertynames(cols)[1]))

struct ColumnsRow{T}
    columns::T # a `Columns` object
    row::Int
end

Base.getproperty(c::ColumnsRow, ::Type{T}, col::Int, nm::Symbol) where {T} = getproperty(getfield(c, 1), T, col, nm)[getfield(c, 2)]
Base.getproperty(c::ColumnsRow, nm::Symbol) = getproperty(getfield(c, 1), nm)[getfield(c, 2)]
Base.propertynames(c::ColumnsRow) = propertynames(getfield(c, 1))

struct RowIterator{T}
    columns::T
    len::Int
end
Base.eltype(x::RowIterator{T}) where {T} = ColumnsRow{T}
Base.length(x::RowIterator) = x.len
schema(x::RowIterator) = schema(x.columns)

function Base.iterate(rows::RowIterator, st=1)
    st > length(rows) && return nothing
    return ColumnsRow(rows.columns, st), st + 1
end

function rows(x::T) where {T}
    if columnaccess(T)
        cols = columns(x)
        return RowIterator(cols, rowcount(cols))
    else
        throw(ArgumentError("no default `Tables.rows` implementation for type: $T"))
    end
end

# build columns from rows
haslength(L) = L isa Union{Base.HasShape, Base.HasLength}

"""
    Tables.allocatecolumn(::Type{T}, len) => returns a column type (usually AbstractVector) w/ size to hold `len` elements

    Custom column types can override with an appropriate "scalar" element type that should dispatch to their column allocator.
"""
allocatecolumn(T, len) = Vector{T}(undef, len)

@inline function allocatecolumns(::Schema{names, types}, len) where {names, types}
    if @generated
        vals = Tuple(:(allocatecolumn($(fieldtype(types, i)), len)) for i = 1:fieldcount(types))
        return :(NamedTuple{names}(($(vals...),)))
    else
        return NamedTuple{names}(Tuple(allocatecolumn(fieldtype(types, i), len) for i = 1:fieldcount(types)))
    end
end

# add! will push! or setindex! a value depending on if the row-iterator HasLength or not
@inline add!(val, col::Int, nm::Symbol, ::Union{Base.HasLength, Base.HasShape}, nt, row) = setindex!(nt[col], val, row)
@inline add!(val, col::Int, nm::Symbol, T, nt, row) = push!(nt[col], val)

@inline function buildcolumns(schema, rowitr::T) where {T}
    L = Base.IteratorSize(T)
    len = haslength(L) ? length(rowitr) : 0
    nt = allocatecolumns(schema, len)
    for (i, row) in enumerate(rowitr)
        eachcolumn(add!, schema, row, L, nt, i)
    end
    return nt
end

@inline add!(dest::AbstractArray, val, ::Union{Base.HasLength, Base.HasShape}, row) = setindex!(dest, val, row)
@inline add!(dest::AbstractArray, val, T, row) = push!(dest, val)

@inline function add_or_widen!(dest::AbstractArray{T}, val::S, nm::Symbol, L, row, len, updated) where {T, S}
    if S === T || promote_type(S, T) <: T
        add!(dest, val, L, row)
        return dest
    else
        new = allocatecolumn(promote_type(T, S), max(len, length(dest)))
        row > 1 && copyto!(new, 1, dest, 1, row - 1)
        add!(new, val, L, row)
        updated[] = merge(updated[], NamedTuple{(nm,)}((new,)))
        return new
    end
end

@inline function add_or_widen!(val, col, nm, L, columns, row, len, updated)
    @inbounds add_or_widen!(columns[col], val, nm, L, row, len, updated)
    return
end

# when Tables.schema(x) === nothing
function buildcolumns(::Nothing, rowitr::T) where {T}
    state = iterate(rowitr)
    state === nothing && return NamedTuple()
    row, st = state
    names = Tuple(propertynames(row))
    L = Base.IteratorSize(T)
    len = haslength(L) ? length(rowitr) : 0
    sch = Schema(names, nothing)
    columns = NamedTuple{names}(Tuple(Union{}[] for _ = 1:length(names)))
    return _buildcolumns(rowitr, row, st, sch, L, columns, 1, len, Ref{Any}(columns))
end

function _buildcolumns(rowitr, row, st, sch, L, columns, rownbr, len, updated)
    while true
        eachcolumn(add_or_widen!, sch, row, L, columns, rownbr, len, updated)
        rownbr += 1
        state = iterate(rowitr, st)
        state === nothing && break
        row, st = state
        columns !== updated[] && return _buildcolumns(rowitr, row, st, sch, L, updated[], rownbr, len, updated)
    end
    return updated[]
end

@inline function columns(x::T) where {T}
    if rowaccess(T)
        r = rows(x)
        return buildcolumns(schema(r), r)
    else
        throw(ArgumentError("no default `Tables.columns` implementation for type: $T"))
    end
end


text:
#-------------------------------------------------------------------------------
# Logger configuration
# TODO: Lots of design work to do here.  Presumably there should be a standard
# textural log system config system like the python one.

parse_level(level) = level
parse_level(level::String) = parse_level(Symbol(lowercase(level)))
function parse_level(level::Symbol)
    if      level == :belowminlevel  return  BelowMinLevel
    elseif  level == :debug          return  Debug
    elseif  level == :info           return  Info
    elseif  level == :warn           return  Warn
    elseif  level == :error          return  Error
    elseif  level == :abovemaxlevel  return  AboveMaxLevel
    else
        throw(ArgumentError("Unknown log level $level"))
    end
end

"""
    configure_logging(args...; kwargs...)

Call `configure_logging` with the current logger, and update cached log
filtering information.
"""
function configure_logging(args...; kwargs...)
    logger = configure_logging(current_logger(), args...; kwargs...)::AbstractLogger
    # FIXME: Tools for setting this should be in Base.
    if haskey(task_local_storage(), :LOGGER_STATE)
        task_local_storage()[:LOGGER_STATE] = LogState(logger)
    else
        global_logger(logger)
    end
    logger
end

configure_logging(::AbstractLogger, args...; kwargs...) =
    throw(ArgumentError("No configure_logging method matches the provided arguments."))

function configure_logging(logger::SimpleLogger, args...; min_level=Info, kwargs...)
    SimpleLogger(logger.stream, parse_level(min_level))
end

disable_logging(level) = disable_logging(parse_level(level))


text:
using ODEFilters
using Test
using LinearAlgebra

using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_fitzhughnagumo, prob_ode_vanstiff

@testset "Smoothing with small constant steps" begin
    prob = ODEFilters.remake_prob_with_jac(prob_ode_fitzhughnagumo)
    @test solve_ieks(prob, IEKS(order=4, diffusionmodel=:fixed)) isa ODEFilters.ProbODESolution
end


text:
"""
    LUtanpaP(A)
adalah fungsi yang digunakan untuk memfaktorisasi matriks `A` sembarang menjadi
dua matriks, yaitu matriks segitiga atas `U` dan bawah `L` sedemikian sehingga `LU=A`

# Examples
```jl
julia> A = [ 2 5 4 2 3
            4 1 0.4 6 5
            0.9 4 8 7 10
            5 7 0.5 2 7
            2 6 8 9 4];

julia> L,U = LUtanpaP(A);

julia> L
5×5 Array{Float64,2}:
 1.0    0.0        0.0       0.0      0.0
 2.0    1.0        0.0       0.0      0.0
 0.45  -0.194444   1.0       0.0      0.0
 2.5    0.611111  -1.02824   1.0      0.0
 1.0   -0.111111   0.668235  1.17806  1.0

julia> U
5×5 Array{Float64,2}:
 2.0   5.0   4.0      2.0        3.0
 0.0  -9.0  -7.6      2.0       -1.0
 0.0   0.0   4.72222  6.48889    8.45556
 0.0   0.0   0.0      2.44988    8.80541
 0.0   0.0   0.0      0.0      -15.1347
```
return matriks `L` dan `U`.
"""
function LUtanpaP(A);
    # Definisikan matriks L sebagai pencatat pengali.
    n,n = size(A);
    L = zeros(n,n);
    # Lakukan operasi baris dasar terhadap matriks
    Aug = copy(A);
    for p = 1:n-1
        # Jika pivot bernilai nol, maka gagal.
        if Aug[p,p]==0
            error("Pivot bernilai nol");
        end
        # Lakukan eliminasi lalu simpan pengali pada matriks L.
        for i = p+1:n
            k = Aug[i,p]/Aug[p,p];
            Aug[i,1:n] = Aug[i,1:n] - k*Aug[p,1:n];
            L[i,p] = k;
        end
    end
    U = Aug;
    L = L .+ I(n);
    return L, U
end


text:
for O in (:Adjoint, :Transpose,:Cholesky, :Hessenberg, :LDLt, :LQ, :LU, :QR)
    @eval LinearAlgebra.$O(::Type{T}) where {T} = LinearAlgebra.$O{T, Array{T,1}}
end

for O in (:BunchKaufman, :GeneralizedSVD, :Schur)
    @eval LinearAlgebra.$O(::Type{T}) where {T} = LinearAlgebra.$O{T, Array{T,2}}
end


text:
using ModelingToolkit, StaticArrays, LinearAlgebra
using DiffEqBase
using Test

# Define some variables
@parameters t σ ρ β
@variables x(t) y(t) z(t)
@derivatives D'~t

eqs = [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       0 ~ x + y + β*z]

lorenz1 = ODESystem(eqs,t,[x,y,z],[σ,ρ,β],name=:lorenz1)
lorenz2 = ODESystem(eqs,t,[x,y,z],[σ,ρ,β],name=:lorenz2)

@parameters α
@variables a(t)
connnectedeqs = [D(a) ~ a*lorenz1.x]

connected1 = ODESystem(connnectedeqs,t,[a],[α],systems=[lorenz1,lorenz2],name=:connected1)

u0map = [lorenz1.x => 1.0,
         lorenz2.x => 2.0,
         lorenz1.y => 3.0,
         lorenz2.y => 4.0,
         lorenz1.z => 5.0,
         lorenz2.z => 6.0,
         a => 7.0]

parammap = SA[lorenz1.σ => 1.0,
              lorenz1.ρ => 2.0,
              lorenz1.β => 3.0,
              lorenz2.σ => 4.0,
              lorenz2.ρ => 5.0,
              lorenz2.β => 6.0,
              α => 7.0]

prob = ODEProblem(connected1,u0map,(0.0,100.0),parammap,jac=true)

@test prob.u0 == [7.0,1.0,3.0,5.0,2.0,4.0,6.0]
@test prob.p == [7.0,1.0,2.0,3.0,4.0,5.0,6.0]
@test prob.u0 isa Vector
@test prob.p isa SVector
@test prob.f.syms == [:a,:lorenz1₊x,:lorenz1₊y,:lorenz1₊z,:lorenz2₊x,:lorenz2₊y,:lorenz2₊z]

eqs_flat = [D(a) ~ a*lorenz1.x,
            D(lorenz1.x) ~ lorenz1.σ*(lorenz1.y-lorenz1.x),
            D(lorenz1.y) ~ lorenz1.x*(lorenz1.ρ-lorenz1.z)-lorenz1.y,
            0 ~ lorenz1.x + lorenz1.y + lorenz1.β*lorenz1.z,
            D(lorenz2.x) ~ lorenz2.σ*(lorenz2.y-lorenz2.x),
            D(lorenz2.y) ~ lorenz2.x*(lorenz2.ρ-lorenz2.z)-lorenz2.y,
            0 ~ lorenz2.x + lorenz2.y + lorenz2.β*lorenz2.z]

@test states(connected1) == convert.(Variable,[a,lorenz1.x,lorenz1.y,lorenz1.z,lorenz2.x,lorenz2.y,lorenz2.z])
@test parameters(connected1) == convert.(Variable,[α,lorenz1.σ,lorenz1.ρ,lorenz1.β,lorenz2.σ,lorenz2.ρ,lorenz2.β])
@test eqs_flat == equations(connected1)

@variables lorenz1₊x(t) lorenz1₊y(t) lorenz1₊z(t) lorenz2₊x(t) lorenz2₊y(t) lorenz2₊z(t)
@parameters lorenz1₊σ lorenz1₊ρ lorenz1₊β lorenz2₊σ lorenz2₊ρ lorenz2₊β

eqs_flat2 = [D(a) ~ a*lorenz1₊x,
            D(lorenz1₊x) ~ lorenz1₊σ*(lorenz1₊y-lorenz1₊x),
            D(lorenz1₊y) ~ lorenz1₊x*(lorenz1₊ρ-lorenz1₊z)-lorenz1₊y,
            0 ~ lorenz1₊x + lorenz1₊y + lorenz1₊β*lorenz1₊z,
            D(lorenz2₊x) ~ lorenz2₊σ*(lorenz2₊y-lorenz2₊x),
            D(lorenz2₊y) ~ lorenz2₊x*(lorenz2₊ρ-lorenz2₊z)-lorenz2₊y,
            0 ~ lorenz2₊x + lorenz2₊y + lorenz2₊β*lorenz2₊z]

@test [x.name for x in states(connected1)] == [:a,:lorenz1₊x,:lorenz1₊y,:lorenz1₊z,:lorenz2₊x,:lorenz2₊y,:lorenz2₊z]
@test [x.name for x in parameters(connected1)] == [:α,:lorenz1₊σ,:lorenz1₊ρ,:lorenz1₊β,:lorenz2₊σ,:lorenz2₊ρ,:lorenz2₊β]
@test eqs_flat == equations(connected1)

connected2 = ODESystem(connnectedeqs,t,[a],[α],systems=[lorenz1,lorenz2],name=:connected2)

@parameters γ
@variables g(t)
connnectedeqs2 = [D(g) ~ g*connected1.lorenz1.x]
doublelevel = ODESystem(connnectedeqs2,t,[g],[γ],systems=[connected1,connected2],name=:doublelevel)

@test states(doublelevel) == convert.(Variable,[g,connected1.a,connected1.lorenz1.x,connected1.lorenz1.y,connected1.lorenz1.z,connected1.lorenz2.x,connected1.lorenz2.y,connected1.lorenz2.z,
                                          connected2.a,connected2.lorenz1.x,connected2.lorenz1.y,connected2.lorenz1.z,connected2.lorenz2.x,connected2.lorenz2.y,connected2.lorenz2.z])

@test parameters(doublelevel) == convert.(Variable,[γ,
                                             connected1.α,connected1.lorenz1.σ,connected1.lorenz1.ρ,connected1.lorenz1.β,connected1.lorenz2.σ,connected1.lorenz2.ρ,connected1.lorenz2.β,
                                             connected2.α,connected2.lorenz1.σ,connected2.lorenz1.ρ,connected2.lorenz1.β,connected2.lorenz2.σ,connected2.lorenz2.ρ,connected2.lorenz2.β])

eqs_flat = [D(g) ~ g*connected1.lorenz1.x,
             D(connected1.a) ~ connected1.a*connected1.lorenz1.x,
             D(connected1.lorenz1.x) ~ connected1.lorenz1.σ*(connected1.lorenz1.y-connected1.lorenz1.x),
             D(connected1.lorenz1.y) ~ connected1.lorenz1.x*(connected1.lorenz1.ρ-connected1.lorenz1.z)-connected1.lorenz1.y,
             0 ~ connected1.lorenz1.x + connected1.lorenz1.y + connected1.lorenz1.β*connected1.lorenz1.z,
             D(connected1.lorenz2.x) ~ connected1.lorenz2.σ*(connected1.lorenz2.y-connected1.lorenz2.x),
             D(connected1.lorenz2.y) ~ connected1.lorenz2.x*(connected1.lorenz2.ρ-connected1.lorenz2.z)-connected1.lorenz2.y,
             0 ~ connected1.lorenz2.x + connected1.lorenz2.y + connected1.lorenz2.β*connected1.lorenz2.z,
             D(connected2.a) ~ connected2.a*connected2.lorenz1.x,
             D(connected2.lorenz1.x) ~ connected2.lorenz1.σ*(connected2.lorenz1.y-connected2.lorenz1.x),
             D(connected2.lorenz1.y) ~ connected2.lorenz1.x*(connected2.lorenz1.ρ-connected2.lorenz1.z)-connected2.lorenz1.y,
             0 ~ connected2.lorenz1.x + connected2.lorenz1.y + connected2.lorenz1.β*connected2.lorenz1.z,
             D(connected2.lorenz2.x) ~ connected2.lorenz2.σ*(connected2.lorenz2.y-connected2.lorenz2.x),
             D(connected2.lorenz2.y) ~ connected2.lorenz2.x*(connected2.lorenz2.ρ-connected2.lorenz2.z)-connected2.lorenz2.y,
             0 ~ connected2.lorenz2.x + connected2.lorenz2.y + connected2.lorenz2.β*connected2.lorenz2.z]

@test eqs_flat == equations(doublelevel)

@test [x.name for x in states(doublelevel)] == [:g,
                                                :connected1₊a,:connected1₊lorenz1₊x,:connected1₊lorenz1₊y,:connected1₊lorenz1₊z,:connected1₊lorenz2₊x,:connected1₊lorenz2₊y,:connected1₊lorenz2₊z,
                                                :connected2₊a,:connected2₊lorenz1₊x,:connected2₊lorenz1₊y,:connected2₊lorenz1₊z,:connected2₊lorenz2₊x,:connected2₊lorenz2₊y,:connected2₊lorenz2₊z]
@test [x.name for x in parameters(doublelevel)] == [:γ,
                                                   :connected1₊α,:connected1₊lorenz1₊σ,:connected1₊lorenz1₊ρ,:connected1₊lorenz1₊β,:connected1₊lorenz2₊σ,:connected1₊lorenz2₊ρ,:connected1₊lorenz2₊β,
                                                   :connected2₊α,:connected2₊lorenz1₊σ,:connected2₊lorenz1₊ρ,:connected2₊lorenz1₊β,:connected2₊lorenz2₊σ,:connected2₊lorenz2₊ρ,:connected2₊lorenz2₊β]

@variables connected1₊a(t) connected1₊lorenz1₊x(t) connected1₊lorenz1₊y(t) connected1₊lorenz1₊z(t) connected1₊lorenz2₊x(t) connected1₊lorenz2₊y(t) connected1₊lorenz2₊z(t)
@variables connected2₊a(t) connected2₊lorenz1₊x(t) connected2₊lorenz1₊y(t) connected2₊lorenz1₊z(t) connected2₊lorenz2₊x(t) connected2₊lorenz2₊y(t) connected2₊lorenz2₊z(t)
@parameters connected1₊α connected1₊lorenz1₊σ connected1₊lorenz1₊ρ connected1₊lorenz1₊β connected1₊lorenz2₊σ connected1₊lorenz2₊ρ connected1₊lorenz2₊β
@parameters connected2₊α connected2₊lorenz1₊σ connected2₊lorenz1₊ρ connected2₊lorenz1₊β connected2₊lorenz2₊σ connected2₊lorenz2₊ρ connected2₊lorenz2₊β

eqs_flat2 = [D(g) ~ g*connected1₊lorenz1₊x,
            D(connected1₊a) ~ connected1₊a*connected1₊lorenz1₊x,
            D(connected1₊lorenz1₊x) ~ connected1₊lorenz1₊σ*(connected1₊lorenz1₊y-connected1₊lorenz1₊x),
            D(connected1₊lorenz1₊y) ~ connected1₊lorenz1₊x*(connected1₊lorenz1₊ρ-connected1₊lorenz1₊z)-connected1₊lorenz1₊y,
            0 ~ connected1₊lorenz1₊x + connected1₊lorenz1₊y + connected1₊lorenz1₊β*connected1₊lorenz1₊z,
            D(connected1₊lorenz2₊x) ~ connected1₊lorenz2₊σ*(connected1₊lorenz2₊y-connected1₊lorenz2₊x),
            D(connected1₊lorenz2₊y) ~ connected1₊lorenz2₊x*(connected1₊lorenz2₊ρ-connected1₊lorenz2₊z)-connected1₊lorenz2₊y,
            0 ~ connected1₊lorenz2₊x + connected1₊lorenz2₊y + connected1₊lorenz2₊β*connected1₊lorenz2₊z,
            D(connected2₊a) ~ connected2₊a*connected2₊lorenz1₊x,
            D(connected2₊lorenz1₊x) ~ connected2₊lorenz1₊σ*(connected2₊lorenz1₊y-connected2₊lorenz1₊x),
            D(connected2₊lorenz1₊y) ~ connected2₊lorenz1₊x*(connected2₊lorenz1₊ρ-connected2₊lorenz1₊z)-connected2₊lorenz1₊y,
            0 ~ connected2₊lorenz1₊x + connected2₊lorenz1₊y + connected2₊lorenz1₊β*connected2₊lorenz1₊z,
            D(connected2₊lorenz2₊x) ~ connected2₊lorenz2₊σ*(connected2₊lorenz2₊y-connected2₊lorenz2₊x),
            D(connected2₊lorenz2₊y) ~ connected2₊lorenz2₊x*(connected2₊lorenz2₊ρ-connected2₊lorenz2₊z)-connected2₊lorenz2₊y,
            0 ~ connected2₊lorenz2₊x + connected2₊lorenz2₊y + connected2₊lorenz2₊β*connected2₊lorenz2₊z]

@test eqs_flat2 == equations(doublelevel)

M = Array(I,15,15)
M[5,5] = false
M[8,8] = false
M[12,12] = false
M[15,15] = false
@test calculate_massmatrix(doublelevel) == M

jac = Num[connected1₊lorenz1₊x 0 g zeros(1,12)
          zeros(7,1) calculate_jacobian(connected1) zeros(7,7)
          zeros(7,8) calculate_jacobian(connected2)]

jac2 = Num[connected1₊lorenz1₊x 0 g zeros(1,12)
           zeros(7,1) ModelingToolkit.namespace_operation.(calculate_jacobian(connected1),connected1.name,:t) zeros(7,7)
           zeros(7,8) ModelingToolkit.namespace_operation.(calculate_jacobian(connected2),connected2.name,:t)]

@test all(isequal.(calculate_jacobian(doublelevel),jac2))

# Try higher hierarchy

@parameters t σ ρ β
@variables x(t) y(t) z(t)
@derivatives D'~t

eqs = [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]

lorenz1 = ODESystem(eqs,name=:lorenz1)
lorenz2 = ODESystem(eqs,name=:lorenz2)

@variables a(t)
@parameters γ
connections = [0 ~ lorenz1.x + lorenz2.y + a*γ]

# Now deviate from tutorial by introducing a 2nd level of connected blocks
# Other issue: How does scoping work? Is it ok to re-use components named lorenz1 and lorenz2, or will the hierarchy be flattened eventually
# connected1 and 2 are pairs of Lorenz oscillators. I'd like to solve multiple pairs together.
connected1 = ODESystem(connections,t,[a],[γ],systems=[lorenz1,lorenz2], name=:connected1)
connected2 = ODESystem(connections,t,[a],[γ],systems=[lorenz1,lorenz2], name=:connected2)

connections2ndLevel = Equation[] # Pairs are not connected
variables2ndLevel = [] # No extra variables
parameters2ndLevel = [] # No extra parameters

totalSystem = ODESystem(connections2ndLevel,t,variables2ndLevel, parameters2ndLevel, systems = [connected1, connected2],
                        name=:totalSystem)

connections2ndLevel = [0 ~ connected1.lorenz1.x + connected2.lorenz1.x]


text:
using Plots
#pyplot()
#plotlyjs()
#gui()

println("Plotting sine")
h = plot(sin, -2pi, pi, label = "sine function")
#h = plot(cos, -2pi, pi, label = "sine function")
#plot(tan, -2pi, pi, label = "sine function")
#show(h)


text:
timestamp(ddf::DateDataFrame) = ddf.timestamp

Base.values(ddf::DateDataFrame) = error("values is removed for DateDataFrame")

# `colnames` is removed in favor of `DataFrames.names`

names(ddf::DateDataFrame) = names(ddf.df)

function Base.diff(ddf::DateDataFrame; padding=false)
    names_df = names(ddf.df)
    df = ddf.df[2:end, :] .- ddf.df[1:end-1, :]
    if padding
        nan_pad = DataFrame(fill(NaN, 1, length(names_df)), names_df)
        df = vcat(nan_pad, df)
        ts = ddf.timestamp
    else
        ts = ddf.timestamp[2:end]
    end
    return DateDataFrame(ts, df)
end

function lead(ddf::DateDataFrame)
    return DateDataFrame(ddf.timestamp[1:end-1], ddf.df[2:end, :])
end

rename(ddf::DateDataFrame, args...) = rename(ddf.df, args...)
rename!(ddf::DateDataFrame, args...) = rename!(ddf.df, args...)

function rolling_mean(arr::AbstractVector, n::Int; padding=false)
    ret_length = length(arr) - n + 1

    ret = Vector{eltype(arr)}(undef, ret_length)
    ret[1] = sum(arr[1:n])
    for i in 2:ret_length
        ret[i] = ret[i-1] + arr[n + i - 1] - arr[i - 1]
    end
    ret = ret ./ n
    if padding
        ret = vcat(fill(NaN, n-1), ret)
    end
    return ret
end

function moving(::typeof(mean), ddf::DateDataFrame, n::Int; padding=false)
    pair_vec = [key=>rolling_mean(ddf.df[!, key], n; padding=padding) for key in names(ddf.df)]
    df = DataFrame(pair_vec)
    ts = padding ? ddf.timestamp : ddf.timestamp[n:end]
    return DateDataFrame(ts, df)
end


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

JLLWrappers.@generate_main_file_header("GAP")
JLLWrappers.@generate_main_file("GAP", UUID("5cd7a574-2c56-5be2-91dc-c8bc375b9ddf"))
end  # module GAP_jll


text:
using ITensors, Zygote
using Zygote: @adjoint
using ITensors: setinds

"""
A finite size PEPS type.
"""
mutable struct PEPS
    data::Matrix{ITensor}
end

PEPS(Nx::Int, Ny::Int) = PEPS(Matrix{ITensor}(undef, Nx, Ny))

"""
    PEPS([::Type{ElT} = Float64, sites; linkdims=1)
Construct an PEPS filled with Empty ITensors of type `ElT` from a collection of indices.
Optionally specify the link dimension with the keyword argument `linkdims`, which by default is 1.
"""
function PEPS(::Type{T}, sites::Matrix{<:Index}; linkdims::Integer = 1) where {T<:Number}
    Ny, Nx = size(sites)
    tensor_grid = Matrix{ITensor}(undef, Ny, Nx)
    # we assume the PEPS at least has size (2,2). Can generalize if necessary
    @assert(Nx >= 2 && Ny >= 2)

    lh = Matrix{Index}(undef, Ny, Nx - 1)
    for ii = 1:(Nx-1)
        for jj = 1:(Ny)
            lh[jj, ii] = Index(linkdims, "Lh,$jj,$ii")
        end
    end
    lv = Matrix{Index}(undef, Ny - 1, Nx)
    for ii = 1:(Nx)
        for jj = 1:(Ny-1)
            lv[jj, ii] = Index(linkdims, "Lv,$jj,$ii")
        end
    end

    # boundary cases
    tensor_grid[1, 1] = emptyITensor(T, lh[1, 1], lv[1, 1], sites[1, 1])
    tensor_grid[1, Nx] = emptyITensor(T, lh[1, Nx-1], lv[1, Nx], sites[1, Nx])
    tensor_grid[Ny, 1] = emptyITensor(T, lh[Ny, 1], lv[Ny-1, 1], sites[Ny, 1])
    tensor_grid[Ny, Nx] = emptyITensor(T, lh[Ny, Nx-1], lv[Ny-1, Nx], sites[Ny, Nx])
    for ii = 2:Nx-1
        tensor_grid[1, ii] =
            emptyITensor(T, lh[1, ii], lh[1, ii-1], lv[1, ii], sites[1, ii])
        tensor_grid[Ny, ii] =
            emptyITensor(T, lh[Ny, ii], lh[Ny, ii-1], lv[Ny-1, ii], sites[Ny, ii])
    end

    # inner sites
    for jj = 2:Ny-1
        tensor_grid[jj, 1] =
            emptyITensor(T, lh[jj, 1], lv[jj, 1], lv[jj-1, 1], sites[jj, 1])
        tensor_grid[jj, Nx] =
            emptyITensor(T, lh[jj, Nx-1], lv[jj, Nx], lv[jj-1, Nx], sites[jj, Nx])
        for ii = 2:Nx-1
            tensor_grid[jj, ii] = emptyITensor(
                T,
                lh[jj, ii],
                lh[jj, ii-1],
                lv[jj, ii],
                lv[jj-1, ii],
                sites[jj, ii],
            )
        end
    end

    return PEPS(tensor_grid)
end

PEPS(sites::Matrix{<:Index}, args...; kwargs...) = PEPS(Float64, sites, args...; kwargs...)

PEPS(data::Array, size::Tuple{Int64,Int64}) = PEPS(reshape(data, size[1], size[2]))

function randomizePEPS!(P::PEPS)
    dimy, dimx = size(P.data)
    for ii = 1:dimx
        for jj = 1:dimy
            randn!(P.data[jj, ii])
            normalize!(P.data[jj, ii])
        end
    end
end

function Base.:+(A::PEPS, B::PEPS)
    @assert(size(A.data) == size(B.data))
    out = [t1 + t2 for (t1, t2) in zip(vcat(A.data...), vcat(B.data...))]
    return PEPS(out, size(A.data))
end

function Base.:-(A::PEPS, B::PEPS)
    @assert(size(A.data) == size(B.data))
    out = [t1 - t2 for (t1, t2) in zip(vcat(A.data...), vcat(B.data...))]
    return PEPS(out, size(A.data))
end

function Base.:*(c::Number, A::PEPS)
    out = [c * t for t in vcat(A.data...)]
    return PEPS(out, size(A.data))
end

function Base.:*(A::PEPS, B::PEPS)
    @assert(size(A.data) == size(B.data))
    out_scalar = 0
    for (t1, t2) in zip(vcat(A.data...), vcat(B.data...))
        out_scalar += ITensors.scalar(t1 * t2)
    end
    return out_scalar
end

function prime(peps::PEPS; ham = true)
    prime_peps = []
    for tensor in vcat(peps.data...)
        indices = inds(tensor)
        if ham == true
            bonds = indices
        else
            bonds = [indices[i] for i = 1:length(indices)-1]
        end
        prime_peps = vcat(prime_peps, [ITensors.prime(tensor, bonds)])
    end
    return PEPS(prime_peps, size(peps.data))
end

@adjoint function prime(peps::PEPS; ham = true)
    dimy, dimx = size(peps.data)
    peps_vec = reshape(peps.data, dimy * dimx)
    function adjoint_pullback(dpeps_prime::PEPS)
        dpeps_prime_vec = reshape(dpeps_prime.data, dimy * dimx)
        dpeps_vec = []
        for i = 1:dimy*dimx
            indices = inds(peps_vec[i])
            indices_reorder = []
            for i_prime in inds(dpeps_prime_vec[i])
                index = findall(x -> x.id == i_prime.id, indices)
                @assert(length(index) == 1)
                push!(indices_reorder, indices[index[1]])
            end
            push!(dpeps_vec, setinds(dpeps_prime_vec[i], Tuple(indices_reorder)))
        end
        dpeps = PEPS(dpeps_vec, (dimy, dimx))
        return (dpeps, nothing)
    end
    return prime(peps; ham = ham), adjoint_pullback
end

# Get the tensor network of <peps|peps'>
function inner_network(peps::PEPS, peps_prime::PEPS)
    return vcat(vcat(peps.data...), vcat(peps_prime.data...))
end

# Get the tensor network of <peps|mpo|peps'>
# The local MPO specifies the 2-site term of the Hamiltonian
function inner_network(
    peps::PEPS,
    peps_prime::PEPS,
    peps_prime_ham::PEPS,
    mpo::MPO,
    coordinates::Array,
)
    @assert(length(mpo) == length(coordinates))
    network = vcat(peps.data...)
    dimy, dimx = size(peps.data)
    for ii = 1:dimx
        for jj = 1:dimy
            if (jj => ii) in coordinates
                index = findall(x -> x == (jj => ii), coordinates)
                @assert(length(index) == 1)
                network = vcat(network, [mpo.data[index[1]]])
                network = vcat(network, [peps_prime_ham.data[jj, ii]])
            else
                network = vcat(network, [peps_prime.data[jj, ii]])
            end
        end
    end
    return network
end

function extract_data(v::Array{<:PEPS})
    tensor_list = [vcat(peps.data...) for peps in v]
    return vcat(tensor_list...)
end

@adjoint function extract_data(v::Array{<:PEPS})
    size_list = [size(peps.data) for peps in v]
    function adjoint_pullback(dt)
        dt = [t for t in dt]
        index = 0
        dv = []
        for (dimy, dimx) in size_list
            size = dimy * dimx
            d_peps = PEPS(reshape(dt[index+1:index+size], dimy, dimx))
            index += size
            push!(dv, d_peps)
        end
        return (dv,)
    end
    return extract_data(v), adjoint_pullback
end