text:
if !@isdefined(TestINet)
    struct TestINet <: ANetwork; end
    struct TestIDiNet <: ADiNetwork; end
    struct TestIIdxEdge <: AIndexedEdge; end
end

@testset "$TEST $G" begin

g = TestINet()
h = TestIDiNet()
e = TestIIdxEdge()

@test_throws ErrorException idx(e)

if G<: ANetwork

g = G(10, 20)
m = ones(10,10)
add_edge_property!(g, "rand", m)
em = eprop(g, "rand")
@test typeof(em) <: EdgeMap
@test em.data === m
@test eprop(g, "rand")[1,1] == 1

m = fill(2, 10)
add_vertex_property!(g, "rand", m)
vm = vprop(g, "rand")
@test typeof(vm) <: VertexMap
@test vm.data === m
@test vprop(g, "rand")[1] == 2

g = G(10, 20)
eprop!(g, "E", EdgeMap(g, e -> rand()))
vprop!(g, "V", VertexMap(g, v -> rand()))
eprop!(g, "E1", e -> rand())
vprop!(g, "V1", v -> rand())

e = first(edges(g))
@test eprop(g, e)["E"] == eprop(g, e, "E")
@test vprop(g, 1)["V"] == vprop(g, 1, "V")

@test eprop(g, src(e), dst(e))["E"] == eprop(g, e, "E")
@test eprop(g, src(e), dst(e))["E"] == eprop(g, src(e), dst(e), "E")

end #if

end # testset


text:
using Test
using Random: randn, rand, randn!
using AutoGrad: @gcheck, @diff, Param, Tape, value
using Knet.Ops20: elu, relu, selu, sigm, dropout, bmm, conv4, conv4w, conv4x, deconv4, mat, pool, poolx, unpool
using Knet.Ops20: batchnorm, bnmoments, bnparams, softmax, logsoftmax, logsumexp, accuracy, nll, bce, logistic, RNN

# Allow the user to override this:
if !isdefined(Main, :testparam)
    testparam(x...) = Param(randn(x...))
end

@testset "ops20" begin

    x = testparam(8,8,2,3)
    @testset "activation" begin
        @test @gcheck elu.(x)
        @test @gcheck relu.(x)
        @test @gcheck selu.(x)
        @test @gcheck sigm.(x)
        @test @gcheck dropout(x, 0.5; seed=1, drop=true)
    end

    @testset "bmm" begin
        b = testparam(8,6,2,3)
        pd(x) = permutedims(x,(2,1,3,4))
        @test @gcheck bmm(x,b)
        @test @gcheck bmm(x,pd(b),transB=true)
        @test @gcheck bmm(pd(x),b,transA=true)
        @test @gcheck bmm(pd(x),pd(b),transA=true,transB=true)
        @test (bmm(x,b) ≈ bmm(x,pd(b),transB=true) ≈ bmm(pd(x),b,transA=true) ≈ bmm(pd(x),pd(b),transA=true,transB=true))
    end

    function convtest(w,x;o...)
        c = conv4(w,x;o...)
        y = testparam(size(c)...)
        ((@gcheck conv4(w,x;o...)) &&
         (@gcheck conv4w(w,x,y;o...)) &&
         (@gcheck conv4x(w,x,y;o...)) &&
         (@gcheck deconv4(w,y;o...)))
    end

    @testset "conv" begin
        w = testparam(3,3,2,4)
        @test convtest(w, x; padding=0, stride=1, dilation=1, mode=0, alpha=1, group=1)
        @test convtest(w, x; padding=1, stride=1, dilation=1, mode=0, alpha=1, group=1)
        @test convtest(w, x; padding=0, stride=2, dilation=1, mode=0, alpha=1, group=1)
        @test convtest(w, x; padding=0, stride=1, dilation=2, mode=0, alpha=1, group=1)
        @test convtest(w, x; padding=0, stride=1, dilation=1, mode=1, alpha=1, group=1)
        @test convtest(w, x; padding=0, stride=1, dilation=1, mode=0, alpha=2, group=1)
        w1 = testparam(3,3,1,4)
        if value(w) isa Array
            # TODO: Grouped convolutions not yet implemented in NNlib, see https://github.com/JuliaGPU/CuArrays.jl/pull/523
            @test_skip convtest(w1, x; padding=0, stride=1, dilation=1, mode=0, alpha=1, group=2)
        else
            @test convtest(w1, x; padding=0, stride=1, dilation=1, mode=0, alpha=1, group=2)
        end
    end

    function pooltest(x; o...)
        y = Param(pool(x; o...))
        dy = testparam(size(y)...)
        ((@gcheck pool(x; o...)) &&
         # @gcheck poolx(x,y,dy) misleads because NNlib gets confused when x&y don't change together.
         (@gcheck poolx(value(x), value(y), dy; o...)) &&
         (@gcheck unpool(x; o...)) &&
         true)
    end

    @testset "pool" begin
        @test pooltest(x; window=2, stride=2, padding=0, mode=0, maxpoolingNanOpt=1, alpha=1)
        @test pooltest(x; window=3, stride=2, padding=0, mode=0, maxpoolingNanOpt=1, alpha=1)
        @test pooltest(x; window=2, stride=3, padding=0, mode=0, maxpoolingNanOpt=1, alpha=1)
        if value(x) isa Array
            # Pool padding is buggy: https://github.com/FluxML/NNlib.jl/issues/229
            @test_skip pooltest(x; window=2, stride=2, padding=1, mode=0, maxpoolingNanOpt=1, alpha=1)
        else
            @test pooltest(x; window=2, stride=2, padding=1, mode=0, maxpoolingNanOpt=1, alpha=1)
        end
        @test pooltest(x; window=2, stride=2, padding=0, mode=1, maxpoolingNanOpt=1, alpha=1)
        @test pooltest(x; window=2, stride=2, padding=0, mode=2, maxpoolingNanOpt=1, alpha=1)
        @test pooltest(x; window=2, stride=2, padding=0, mode=3, maxpoolingNanOpt=1, alpha=1)
        @test pooltest(x; window=2, stride=2, padding=0, mode=0, maxpoolingNanOpt=0, alpha=1)
        @test pooltest(x; window=2, stride=2, padding=0, mode=0, maxpoolingNanOpt=1, alpha=2)
    end

    @testset "mat" begin
        @test @gcheck mat(x)
        @test size(mat(x)) == (128,3)
        @test size(mat(x,dims=4)) == (384,1)
        @test size(mat(x,dims=3)) == (128,3)
        @test size(mat(x,dims=2)) == (64,6)
        @test size(mat(x,dims=1)) == (8,48)
        @test size(mat(x,dims=0)) == (1,384)
    end

    @testset "batchnorm" begin
        p = testparam(2*size(x,3))
        m = bnmoments()
        @test (@diff sum(batchnorm(x,m,p))) isa Tape
        @test @gcheck batchnorm(x,m,p; training=false)
    end

    function softtest(f,x)
        ((@gcheck f(x)) &&
         (@gcheck f(x,dims=1)) &&
         (@gcheck f(x,dims=2)))
    end

    @testset "softmax/loss" begin
        scores = testparam(10,100)
        @test softtest(softmax, scores)
        @test softtest(logsoftmax, scores)
        @test softtest(logsumexp, scores)

        labels1 = rand(1:10,100)
        labels2 = rand(1:100,10)
        @test @gcheck nll(scores, labels1, dims=1)
        @test @gcheck nll(scores, labels2, dims=2)

        value(scores)[1,:] .+= 1000
        @test accuracy(scores, labels1, dims=1) == sum(labels1 .== 1)/length(labels1)
        value(scores)[:,1] .+= 10000
        @test accuracy(scores, labels2, dims=2) == sum(labels2 .== 1)/length(labels2)

        scores = testparam(100)
        labels = rand(Bool,100)
        xlabels = labels .* 2 .- 1
        @test @gcheck bce(scores, labels)
        @test @gcheck logistic(scores, xlabels)
        @test bce(scores, labels) == logistic(scores, xlabels)
    end

    function rnntest(;ndims=1, batchSizes=nothing, o...)
        x = testparam(4:(4+ndims-1)...)
        r = RNN(4,4;atype=typeof(value(x)),o...) # need H==X for skipInput test
        @gcheck r(x; batchSizes=batchSizes)
    end

    @testset "rnn" begin
        @test rnntest(rnnType=:lstm, numLayers=1, bidirectional=false, skipInput=false, ndims=3, batchSizes=nothing)
        @test rnntest(rnnType=:gru,  numLayers=1, bidirectional=false, skipInput=false, ndims=3, batchSizes=nothing)
        @test rnntest(rnnType=:relu, numLayers=1, bidirectional=false, skipInput=false, ndims=3, batchSizes=nothing)
        @test rnntest(rnnType=:tanh, numLayers=1, bidirectional=false, skipInput=false, ndims=3, batchSizes=nothing)
        @test rnntest(rnnType=:lstm, numLayers=2, bidirectional=false, skipInput=false, ndims=3, batchSizes=nothing)
        @test rnntest(rnnType=:lstm, numLayers=1, bidirectional=true,  skipInput=false, ndims=3, batchSizes=nothing)
        @test rnntest(rnnType=:lstm, numLayers=1, bidirectional=false, skipInput=true,  ndims=3, batchSizes=nothing)
        @test rnntest(rnnType=:lstm, numLayers=1, bidirectional=false, skipInput=false, ndims=2, batchSizes=nothing)
        @test rnntest(rnnType=:lstm, numLayers=1, bidirectional=false, skipInput=false, ndims=1, batchSizes=nothing)
        if value(x) isa Array
            # TODO: batchSizes not implemented for CPU yet
            @test_skip rnntest(rnnType=:lstm, numLayers=1, bidirectional=false, skipInput=false, ndims=2, batchSizes=[2,2,1])
        else
            @test rnntest(rnnType=:lstm, numLayers=1, bidirectional=false, skipInput=false, ndims=2, batchSizes=[2,2,1])
        end
    end
end

ext: jl


text:
# Autogenerated wrapper script for librttopo_jll for i686-w64-mingw32
export librttopo

using GEOS_jll
JLLWrappers.@generate_wrapper_header("librttopo")
JLLWrappers.@declare_library_product(librttopo, "librttopo-1.dll")
function __init__()
    JLLWrappers.@generate_init_header(GEOS_jll)
    JLLWrappers.@init_library_product(
        librttopo,
        "bin\\librttopo-1.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# module Bukdu.Plug

struct Static <: AbstractPlug
end

struct StaticController <: ApplicationController
    conn::Conn
end

function content_type_for_file_extionsion(ext)::String
    mimes = Dict(
        ".html" => "text/html; charset=utf-8",
        ".htm"  => "text/html; charset=utf-8",
        ".wasm" => "application/wasm",
        ".css"  => "text/css",
    )
    get(mimes, ext, "application/octet-stream")
end

"""
    plug(::Type{Static}; at::String, from::String, only::Union{Vector{String},Nothing}=nothing, indexfile="index.html")
"""
function plug(::Type{Static}; at::String, from::String, only::Union{Vector{String},Nothing}=nothing, indexfile="index.html")
    @info :from1 from
    function _readfile_base(c::StaticController, f)
        reqpath = c.conn.request.target
        offset = isdirpath(at) ? 1 : 2
        targetpath = reqpath[length(at)+offset:end]
        filepath = joinpath(f(from, targetpath)...)
        @info :from from
        @info :targetpath targetpath
        @info :filepath filepath
        (_, fileext) = splitext(filepath)
        ext = lowercase(fileext)
        Render(content_type_for_file_extionsion(ext), read(filepath))
    end # function _readfile_base

    function readfile(c::StaticController)
        _readfile_base(c, (from, targetpath) -> (from, targetpath))
    end # function readfile

    function readindexfile(c::StaticController)
        _readfile_base(c, (from, targetpath) -> (from, targetpath, indexfile))
    end # function readindexfile

    has_only = only isa Vector{String} && !isempty(only)
    for (root, dirs, files) in walkdir(from)
         subpath = root[length(from)+1:end]
         for filename in files
             if has_only
                 subfilepath = normpath(subpath, filename)
                 !any(x->startswith(subfilepath, x), only) && continue # for filename
             end
             if filename == indexfile
                 reqindex = normpath(at, subpath)
                 get(reqindex, StaticController, readindexfile)
             end
             reqpath = normpath(at, subpath, filename)
             if Sys.iswindows()
                 path_arr = split(reqpath, Base.Filesystem.path_separator)
                 reqpath = join(path_arr, '/')
             end
             @info :reqpath_________ (reqpath, at, subpath, filename)
             get(reqpath, StaticController, readfile)
        end
    end
end

# module Bukdu.Plug

max_stars_repo_name: wookay/TestBukdu

max_issues_repo_name: wookay/TestBukdu

max_forks_repo_name: wookay/TestBukdu


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

"""
    module Simulate_PendulumDAE2

Simulate PendulumDAE model with x_fixed=true.
"""
module Simulate_PendulumDAE2

include(joinpath("models", "PendulumDAE.jl"))
import .PendulumDAE

import ModiaMath

model  = PendulumDAE.Model(x_fixed=true)
result = ModiaMath.simulate!(model, stopTime=2.0, log=true)

ModiaMath.plot(result, [(:x, :y), (:vx, :vy), :lambda], heading="Simulate_PendulumDAE2.jl (index3 reduced to index1)")

end


text:
# ------------ Lightweight visitor ------------ #

struct Visitor <: Selection
    addrs::Vector{Any}
    Visitor() = new([])
end

push!(vs::Visitor, addr::Address) = push!(vs.addrs, addr)

function visit!(vs::Visitor, addr)
    addr in vs.addrs && error("VisitorError (visit!): already visited address $(addr).")
    push!(vs, addr)
end

function visit!(vs::Visitor, addrs::Vector)
    for addr in addrs
        addr in vs.addrs && error("VisitorError (visit!): already visited address $(addr).")
        push!(vs, addr)
    end
end

function visit!(vs::Visitor, par::Address, addrs::Vector)
    for addr in addrs
        addr in vs.addrs && error("VisitorError (visit!): already visited address $(addr).")
        push!(vs, par => addr)
    end
end


text:
mutable struct SidePanelManager <: GtkNotebook

    handle::Ptr{Gtk.GObject}
    main_window
    panels

    function SidePanelManager(main_window)

        ntb = GtkNotebook()
        set_gtk_property!(ntb, :vexpand, true)

        # for some strange reason panel aren't type correctly, so I also have them in
        # that variable as well
        panels = Any[]
        n = new(ntb.handle, main_window, panels)
        Gtk.gobject_move_ref(n, ntb)
    end
end

function init!(sp::SidePanelManager, panels, titles)
    for (p,t) in zip(panels, titles)
        push!(sp, p)
        push!(sp.panels, p)
        Gtk.GAccessor.tab_label_text(sp, p, t)
    end
end

function give_me_a_treeview(n, rownames)

    t = ntuple(i->AbstractString, n)
    list = GtkTreeStore(t...)

    tv = GtkTreeView(GtkTreeModel(list))

    cols = GtkTreeViewColumn[]

    for i=1:n
        r1 = GtkCellRendererText()
        c1 = GtkTreeViewColumn(rownames[i], r1, Dict([("text", i-1)]))
        Gtk.G_.sort_column_id(c1, i-1)
        push!(cols, c1)
        Gtk.G_.max_width(c1, Int(200/n))
        push!(tv, c1)
    end

    return (tv, list, cols)
end


text:
"""
    IntroJuliaDemo
Este modulo define funciones de prueba para mostrar cómo funciona la
documentación.
"""
module IntroJuliaDemo

# Write your package code here.

export f

"""
    f(x, y)
Esta función calcula un plano ``2x + 3y`` con la función `f`.
"""
f(x, y) = 2x + 3y

end


text:
using DataFrames
using Distributions

USER = 1000
SKILL = 5
ITEM = SKILL*5
RECORD = 200000
INCREASE = 400 / (RECORD / USER)
DECREASE = 200 / (RECORD / USER)

srand(1)

# 項目の作成
function make_items()
  item = DataFrame(item = Int[], skill = Char[], a = Float64[], β = Float64[], c = Float64[])

  for i in 1:ITEM
    p = 0.1
    a = rand(LogNormal(0, 0.5))
    β = rand(Normal(0,2))
    c = rand(Beta(20*p+1, 20*(1-p)+1))
    push!(item, [i, Char(Int('A') + (i-1) % SKILL), a, β, c])
  end

  return item
end

# 能力値の初期値
function init_trait()
  mat = rand(Normal(-10,2), USER, SKILL)
  return mat
end

# Raschモデル
function rasch(Θ, β, a, c)
  numer = 1 - c
  denom = 1 + exp(-1.7 * a * (Θ - β))
  return c + numer / denom
end

# 忘却型Raschモデル
function decay_rasch(Θ, β, Θ̃, β̃, t, h, c)
  numer = (1 + h * t)^(-exp(Θ̃ - β̃)) - c
  denom = 1 + exp(-1.7(Θ - β))
  return c + numer / denom
end

# 回答結果
function calc_result(trait, item)
  prob = rasch(trait, item[:β][1], item[:a][1], item[:c][1])
  if prob[1] > rand()
    return 1
  else
    return 0
  end
end

# 経過時間(hour)
function passed_hour(second1, second2)
  diff = second1 - second2
  return round(Int(Dates.Second(diff)) / (60*60))
end

# 能力値を成長させる
function grow_trait(trait)
  trait += rand(Beta(10,10)) * INCREASE
  return trait
end

# 能力値を衰退させる
function decline_trait(trait, time)
  trait -= rand(Beta(10,10)) * (DECREASE * time)
  return trait
end

# 初期化
items = make_items()
trait = init_trait()

# レコードを生成
data = DataFrame(result = Int[], user_id = Int[], item_id = Int[], skill = Char[], date = DateTime[])
u_num = zeros(USER)
i_num = zeros(ITEM)

# スタート時刻
date = DateTime(2017, 5, 5, 12, 0, 0)

# 最後に回答した日時を記録
user_last_update = fill(date, USER, SKILL)

for record in 1:RECORD
  user_id = rand(1:USER)
  u_num[user_id] += 1
  item_id = rand(1:ITEM)
  i_num[item_id] += 1
  skill_id = items[:skill][item_id] - 'A' + 1

  passed_time = passed_hour(date, user_last_update[user_id, skill_id])
  trait[user_id, skill_id] = decline_trait(trait[user_id, skill_id], passed_time)

  result = calc_result(trait[user_id, skill_id], items[item_id, :])

  # ランダムな秒数だけ経過させる
  date += Dates.Second(round(rand(Uniform(1,60))))
  user_last_update[user_id, skill_id] = date

  push!(data, [result, user_id, item_id, items[:skill][item_id], date])

  trait[user_id, skill_id] = grow_trait(trait[user_id, skill_id])
end

writetable("synthetic_data.csv", data)


text:
import Base: iterate

export SegmentedSieve

function generate_siever_primes(small_sieve::SmallSieve, segment_lo)
    xs = small_sieve.xs
    sievers = Vector{Siever}(undef, vec_count_ones(xs))
    j = 0
    @inbounds for i = eachindex(xs)
        x = xs[i]
        while x != 0x00
            sievers[j += 1] = Siever(compute_prime(x, i), segment_lo)
            x &= x - 0x01
        end
    end
    return sievers
end

struct SegmentIterator{T<:AbstractUnitRange}
    range::T
    segment_length::Int
    first_byte::Int
    last_byte::Int
    sievers::Vector{Siever}
    presieve_buffer::Vector{UInt8}
    segment::Vector{UInt8}
end

struct Segment{Tr,Ts}
    range::Tr
    segment::Ts
end

function Base.show(io::IO, s::Segment)
    # compute left padding
    padding = floor(Int, log10(last(s.range))) + 1

    padding_str = " " ^ padding

    print(io, padding_str, "  ")
    for p in ps
        print(io, lpad(p, 2, "0"), "  ")
    end

    println()

    for (start, byte) in zip(s.range, s.segment)
        mask = 0b00000001
        print(lpad(start, padding, "0"), "  ")
        for i = 1 : 8
            print(io, (byte & mask) == mask ? " x  " : " .  ")
            mask <<= 1
        end
        println()
    end

    io
end

function SegmentIterator(range::T, segment_length::Integer) where {T<:AbstractUnitRange}
    from, to = first(range), last(range)
    first_byte, last_byte = cld(first(range), 30), cld(last(range), 30)
    sievers = generate_siever_primes(SmallSieve(isqrt(to)), 30 * (first_byte - 1) + 1)
    presieve_buffer = create_presieve_buffer()
    xs = zeros(UInt8, segment_length)

    return SegmentIterator{T}(range, segment_length, first_byte, last_byte, sievers, presieve_buffer, xs)
end

function iterate(iter::SegmentIterator, segment_index_start = iter.first_byte)
    @inbounds begin
        if segment_index_start ≥ iter.last_byte
            return nothing
        end

        from, to = first(iter.range), last(iter.range)

        segment_index_next = min(segment_index_start + iter.segment_length, iter.last_byte + 1)
        segment_curr_len = segment_index_next - segment_index_start

        # Presieve
        apply_presieve_buffer!(iter.segment, iter.presieve_buffer, segment_index_start, segment_index_next - 1)

        # Set the preceding so many bits before `from` to 0
        if segment_index_start == iter.first_byte
            if iter.first_byte === 1
                iter.segment[1] = 0b11111110 # just make 1 not a prime.
            end
            for i = 1 : 8
                30 * (segment_index_start - 1) + ps[i] >= from && break
                iter.segment[1] &= wheel_mask(ps[i])
            end
        end

        # Set the remaining so many bits after `to` to 0
        if segment_index_next == iter.last_byte + 1
            for i = 8 : -1 : 1
                to ≥ 30 * (segment_index_next - 2) + ps[i] && break
                iter.segment[segment_curr_len] &= wheel_mask(ps[i])
            end
        end

        # Sieve the interval, but skip the pre-sieved primes
        xs = iter.segment

        for p_idx in 5:length(iter.sievers)
            p            = iter.sievers[p_idx]
            last_idx     = 0
            n_bytes      = segment_index_next - segment_index_start
            byte_idx     = p.byte_index - segment_index_start + 1
            wheel_idx    = p.wheel_index
            increment    = p.prime_div_30
            @sieve_loop :unroll :save_on_exit
            iter.sievers[p_idx] = Siever(increment, segment_index_start + byte_idx - 1, last_idx)
        end

        segment_start = 30 * (segment_index_start - 1)
        segment_stop = 30 * (segment_index_next - 1) - 1

        segment_index_start += iter.segment_length

        return Segment(segment_start:30:segment_stop, view(xs, Base.OneTo(segment_curr_len))), segment_index_start
    end
end


text:
using HypothesisTests
using Tullio
using Bootstrap
using StatsBase

# Give two feature matrices, a test function (from HypothesisTests.jl) and calculate the test statistics for each feature common to the two feature matrices
function testFeatureDistributions(𝐟₁::AbstractVector, 𝐟₂::AbstractVector, test=TailedTestOfVariance; tail=:both)
    if tail ∈ (:firstgreater, :secondsmaller, :firstlarger)
        tail = :right
    elseif tail ∈ (:firstsmaller, :secondgreater, :secondlarger)
        tail = :left
    end
    p = pvalue(test(𝐟₁, 𝐟₂); tail=tail)
end
function testFeatureDistributions(F₁::AbstractArray, F₂::AbstractArray, test=TailedTestOfVariance; tail=:both)
    @tullio p[i] := testFeatureDistributions(F₁[i, :], F₂[i, :], test; tail=tail)
end
function testFeatureDistributions(F₁::AbstractFeatureArray, F₂::AbstractFeatureArray, test=TailedTestOfVariance; tail=:both)
    # Intersect features before testing. In this case the test statistic returned refers to the order of the first input array. NaN means F₂ did not have a feature in F₁
    F̂₁, F̂₂ = intersectFeatures(F₁, F₂)
    p = testFeatureDistributions(Array(F̂₁), Array(F̂₂), test; tail=tail)
    F̂₁, p = intersectFeatures(F̂₁, p) # label p with the features of F̂₁
    fs = Catch22.featureDims(F̂₁)
    outp = deepcopy(F̂₁[:, 1])
    outp[:] .= NaN
    outp[fs] .= p[fs]
end
export testFeatureDistributions

function TailedTestOfVariance(x::AbstractVector{T}, y::AbstractVector{T}; scorediff=abs, statistic=mean) where {T<:Real}
    zx, zy = scorediff.(x .- statistic(x)), scorediff.(y .- statistic(y))
    EqualVarianceTTest(zx, zy)
end
export TailedTestOfVariance
function RobustTailedTestOfVariance(x::AbstractVector{T}, y::AbstractVector{T}) where {T<:Real}
    zx, zy = abs.(x .- statistic(x)), scorediff.(y .- median(y))
    EqualVarianceTTest(zx, zy)
end
export RobustTailedTestOfVariance

"""
Bootstrap to get an estimate of the standard error of a sample's standard deviation
"""
function bootstrapSEσ(X::AbstractVector, n=1000, sampling=BasicSampling)
    bs = bootstrap(StatsBase.std, X, sampling(n))
    return stderror(bs)[1]
end
export bootstrapSEσ


text:
# sample sequences of ergms with different parameters' values from R package ergm
# and test the PseudoLikelihoodScoreDrivenERGM filter
using Utilities,AReg,StaticNets,DynNets , JLD,MLBase,StatsBase,CSV, RCall
 using PyCall; pygui(); using PyPlot
 # load the required packages in R
 R"library(statnet)
    library(ergm)
    library(sna)
    library(coda)
    library(network)
    sessionInfo()"

 R"set.seed(0)"

 parDgpT =  dgp_paper_SDERGM(Nsample = 100,dgpType = "sin", T = 50, N = 50,
                              Nterms = 2, Nsteps1 = 1 ,Nsteps2 = 0)

# For each t, sample the ERGM with parameters corresponding to the DGP at time t
# and run a single snapeshot estimate
@rput T; @rput parDgpT;@rput N;@rput Nsample
 #create an empty network, the formula defining ergm, sample the ensemble and
 # store the sufficient statistics and change statistics in R
 R"
 net <- network.initialize(N)
  #formula_ergm = net ~ edges + gwesp(decay = 0.25,fixed = TRUE)
  formula_ergm = net ~ edges +  gwesp(decay = 0.25,fixed = TRUE)
   sampledMat_T_R =    array(0, dim=c(N,N,T,Nsample))
   estParSS_T_R =    list()
   changeStats_T_R = list()
   stats_T_R = list()
    for(t in 1:T){
        estParSS_t_R =    list()
        changeStats_t_R = list()
        stats_t_R = list()
        print(t)
        for(n in 1:Nsample){
            print(parDgpT[,t])
             net <- simulate(formula_ergm, nsim = 1, seed = sample(1:100000000,1), coef = parDgpT[,t],control = control.simulate.formula(MCMC.burnin = 100000))
             sampledMat_T_R[,,t,n] <- as.matrix.network( net)
             tmp <- ergm(formula_ergm)#,estimate = 'MPLE')#)#
             estParSS_t_R[[n]] <- tmp[[1]]
             print(c(t,n))

             print(estParSS_t_R[[n]])
             chStat_t <- ergmMPLE(formula_ergm)
             changeStats_t_R[[n]] <- cbind(chStat_t$response, chStat_t$predictor,chStat_t$weights)
             stats_t_R[[n]] <- summary(formula_ergm)

             }
              estParSS_T_R[[t]] <- estParSS_t_R
              changeStats_T_R[[t]] <-changeStats_t_R
              stats_T_R[[t]] <- stats_t_R
         }"

 # import sampled networks in julia
 sampledMat_T = BitArray( @rget(sampledMat_T_R))
 estParSS_T = @rget(estParSS_T_R);#tmp = zeros(Nterms,T); for t=1:T tmp[:,t] = estParSS_T[t]; end ; estParSS_T = tmp
 changeStats_T = @rget changeStats_T_R;# tmp = Array{Array{Float64,2}}(T); for t=1:T tmp[t] =  changeStats_T[t];end;changeStats_T = tmp
 stats_T = @rget(stats_T_R); #tmp = zeros(Nterms,T); for t=1:T tmp[:,t] = stats_T[t];end;stats_T = tmp
 estParSS_T = Utilities.collapseArr3(estParSS_T)

# Save all Sampled Data
save_fold = "./data/estimatesTest/sdergmTest/R_MCMC_estimates/"
 @save(save_fold*"test_Nodes_$(N)_T_$(T)_Sample_$(Nsample)_Ns_" * dgpType * "_$(Nsteps1)_$(Nsteps2)_MPLE.jld",
              stats_T, changeStats_T,estParSS_T,sampledMat_T ,parDgpT,Nsample,T,N)

using Utilities,AReg,StaticNets,DynNets , JLD,MLBase,StatsBase,CSV, RCall
 using PyCall; pygui(); using PyPlot
 using JLD2,Utilities, GLM
 ## load R MCMC simulation and estimates and estimate sdergmTest

 if false
 Nsample = 100
 dgpType = "sin"
 T = 50
 N = 50
 Nterms = 2
 Nsteps1 ,Nsteps2 = 0,1
 load_fold = "./data/estimatesTest/sdergmTest/R_MCMC_estimates/"
 @load(load_fold*"test_Nodes_$(N)_T_$(T)_Sample_$(Nsample)_Ns_" * dgpType * "_$(Nsteps1)_$(Nsteps2)_MPLE.jld",
             stats_T, changeStats_T,estParSS_T,sampledMat_T ,parDgpT,Nsample)
 end

 onlyTest = false
 targetAllTv = true
 tvParFromTest = false; pValTh = 0.05
 gasParEst = fill(fill(Float64[],2), Nsample)
 startPointEst = fill(fill(0.0,2), Nsample)
 staticEst = fill(fill(0.0,2), Nsample)
 filtPar_T_Nsample = zeros(Nterms,T,Nsample)
 pVals_Nsample = zeros(Nterms,Nsample)
 scoreAutoc_Nsample = zeros(Nterms,Nsample)
 convFlag = falses(Nsample)
 tmp = Array{Array{Float64,2},2}(T,Nsample); for t=1:T,s=1:Nsample tmp[t,s] =  changeStats_T[t][s];end;#changeStats_T = tmp
 for n=1:Nsample
     @show(n)
      pVals_Nsample[:,n],tmpInfo,staticEst[n] =
         DynNets.pValStatic_SDERGM( DynNets.SdErgmDirBinGlobalPseudo(tmp[:,n],DynNets.fooGasPar,falses(Nterms),""))
      @show(tmpInfo)
     if tvParFromTest
         indTvPar =  pVals_Nsample[:,n] .< pValTh
     else
         indTvPar = BitArray([true,true])
     end
     if targetAllTv
         indTargPar = indTvPar
     else
          indTargPar =BitArray([false,false])
     end
    model = DynNets.SdErgmDirBinGlobalPseudo(tmp[:,n],DynNets.fooGasPar,indTvPar,"")

    if .!onlyTest
    estPar,convFlag[n],UM,startPoint = DynNets.estimate(model;UM =meanSq(estParSS_T[n,:,:],1),indTvPar = indTvPar, indTargPar = indTargPar)

    gasParEst[n] = estPar # store gas parameters
    startPointEst[n] = startPoint
    gasParVec = zeros(sum(indTvPar)*3); for i=0:(sum(indTvPar)-1) gasParVec[1+3i : 3(i+1)] = estPar[indTvPar][i+1]; end
    constParVec = zeros(sum(!indTvPar)); for i=1:sum(.!indTvPar) constParVec[i] = estPar[.!indTvPar][i][1]; end
    gasFiltPar , pseudolike = DynNets.score_driven_filter_or_dgp(model,gasParVec,indTvPar;
                              vConstPar = constParVec,ftot_0= startPoint )

    filtPar_T_Nsample[:,:,n] = gasFiltPar
    end
 end

 save_fold = "./data/estimatesTest/sdergmTest/"*
            "gas_MCMC_comparison_estimates/"
 @save(save_fold*"test_Nodes_$(N)_T_$(T)_Sample_$(Nsample)_Ns_" * dgpType * "_$(Nsteps1)_$(Nsteps2)_MPLE_target_$(targetAllTv).jld" ,
     stats_T, changeStats_T,estParSS_T,sampledMat_T ,parDgpT,
     Nsample,T,N,filtPar_T_Nsample,gasParEst,convFlag,pVals_Nsample,scoreAutoc_Nsample,staticEst)

if false

    sum(convFlag)
    figure()
     legTex = ["Constant ERGM";"SDERGM"; "Sequence of ERGM"; "DGP"]
     for n=1:Nsample
     indTvPar = BitArray([true,true])
     model = SdErgmDirBinGlobalPseudo(tmp[:,n],fooGasPar,indTvPar,"")
     estPar = gasParEst[n] # store gas parameters
     gasFiltPar  = filtPar_T_Nsample[:,:,n]
     parInd = 1
     subplot(1,2,1);plot(1:T,ones(T)*staticEst[n][parInd],"-b")
                    plot(1:T,gasFiltPar[parInd,:],"r")
                    plot(1:T,estParSS_T[n,:,parInd],".b")
                    plot(1:T,parDgpT[parInd,:],"k",linewidth=5)

     parInd = 2
     subplot(1,2,2);plot(1:T,ones(T)*staticEst[n][parInd],"-b")
                    plot(1:T,gasFiltPar[parInd,:],"r")
                    plot(1:T,estParSS_T[n,:,parInd],".b")
                    plot(1:T,parDgpT[parInd,:],"k",linewidth=5)
     end
     namePar1 = "Number of Links"
     namePar2 = "GWESP"
     subplot(1,2,1);    title(namePar1); legend(legTex)
     subplot(1,2,2);  title(namePar2 ); legend(legTex)

    #Plotta info relative ai tests
    figure()
         subplot(2,2,1); parInd = 1;   plt[:hist]((pVals_Nsample[parInd,:]), bins=logspace(minimum(log10(pVals_Nsample[1,:])),maximum(log10(pVals_Nsample[parInd,:])), 20))
         if maximum((pVals_Nsample[parInd,:])) >pValTh
              axvspan((pValTh),maximum((pVals_Nsample[parInd,:])),color = "r",alpha = 0.1);
         end
              xscale("log");
         subplot(2,2,2);  parInd = 2;   plt[:hist]((pVals_Nsample[parInd,:]), bins=logspace(minimum(log10(pVals_Nsample[1,:])),maximum(log10(pVals_Nsample[parInd,:])), 20))
         if maximum((pVals_Nsample[parInd,:])) >pValTh
              axvspan((pValTh),maximum((pVals_Nsample[parInd,:])),color = "r",alpha = 0.1);
         end
         subplot(2,2,3); plot((pVals_Nsample[1,:]),(pVals_Nsample[2,:]),".")
         axvline((pValTh));axhline((pValTh))
         xscale("log");
         xscale("log");yscale("log")
         xlabel("P-Value " * namePar1)
         ylabel("P-Value " * namePar2)
end

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
# Bayesian probit model with a Normal prior N(0, priorstd^2*I). The log-posterior and its gradient, as well as the
# metric tensor and its partial derivatives are passed to the model as functions

using Distributions, MCMC

# Create design matrix X and response variable y from vaso data array
vaso = readdlm(joinpath(dirname(@__FILE__),"vaso.txt"), ' ')
covariates, y = vaso[:, 1:end-1], vaso[:, end]
nsamples, npars = size(covariates)
covariates = (broadcast(-, covariates, mean(covariates, 1))./repmat(std(covariates, 1), nsamples, 1))
polynomialOrder = 1
X = ones(nsamples, npars*polynomialOrder+1)
for i = 1:polynomialOrder
  X[:, ((i-1)*npars+2):(i*npars+1)] = covariates.^i
end
npars += 1

# Define log-prior
priorstd = 10.
priorvar = priorstd^2
distribution = MvNormal(zeros(size(X, 2)), priorvar*eye(npars))
logprior(pars::Vector{Float64}) = logpdf(distribution, pars)
randprior() = rand(distribution)

# Define log-likelihood
function loglik(pars::Vector{Float64})
  XPars = X*pars
  normal = Normal()
  return (dot(logcdf(normal, XPars), y)+dot(logcdf(normal, -XPars), 1-y))
end

# Define log-posterior
log_posterior(pars::Vector{Float64}) = logprior(pars)+loglik(pars)

# Define gradient of log-posterior
function grad_log_posterior(pars::Vector{Float64})
  XPars = X*pars
  normal = Normal()
  return (X'*(y.*exp(-(XPars.^2+log(2*pi))/2-logcdf(normal, XPars))
    -(1-y).*exp(-(XPars.^2+log(2*pi))/2-logcdf(normal, -XPars)))-pars/priorvar)
end

# Define metric tensor
function tensor(pars::Vector{Float64})
  XPars = X*pars
  normal = Normal()
  vector = exp(-XPars.^2-logcdf(normal, XPars)-logcdf(normal, -XPars)-log(2*pi))
  return ((X'.*repmat(vector', npars, 1))*X+(eye(npars)/priorvar))
end

# Define derivatives of metric tensor
function deriv_tensor(pars::Vector{Float64})
  output = Array(Float64, npars, npars, npars)

  XPars = X*pars
  normal = Normal()
  vector01 = exp(-XPars.^2-2*logcdf(normal, XPars)-logcdf(normal, -XPars)-log(2*pi))

  for i = 1:npars
    vector02 = (vector01.*(exp(-(XPars.^2+log(2*pi))/2-logcdf(normal, -XPars))
      -2*(pdf(normal, XPars)+XPars.*cdf(normal, XPars))).*X[:, i])

    output[:, :, i] = (X'.*repmat(vector02', npars, 1))*X
  end

  return output
end

mcmodel = model(log_posterior, grad=grad_log_posterior, tensor=tensor, dtensor=deriv_tensor, init=randprior())

mcchain01 = run(mcmodel * RWM(0.5) * SerialMC(1001:10000))

acceptance(mcchain01)

mcchain02 = run(mcmodel * HMC(0.001) * SerialMC(1001:10000))

acceptance(mcchain02)

mcchain03 = run(mcmodel * RMHMC(0.5, EmpMCTuner(0.8, verbose=true)) * SerialMC(5001:10000))

acceptance(mcchain03)


text:
# FIFO queue

type Queue{S}   # S is the type of internal deque
    store::S
end

Queue{T}(ty::Type{T}) = Queue(Deque{T}())
Queue{T}(ty::Type{T}, blksize::Integer) = Queue(Deque{T}(blksize))

isempty(s::Queue) = isempty(s.store)
length(s::Queue) = length(s.store)

front(s::Queue) = front(s.store)
back(s::Queue) = back(s.store)

function enqueue!(s::Queue, x)
    push!(s.store, x)
    s
end

dequeue!(s::Queue) = shift!(s.store)

max_stars_repo_name: ambujone/DIDebug


text:
push!(LOAD_PATH, "../src/")

using Documenter, IterativeLQR

makedocs(
    modules = [IterativeLQR],
    format = Documenter.HTML(prettyurls=false),
    sitename = "IterativeLQR",
    pages = [
        ##############################################
        ## MAKE SURE TO SYNC WITH docs/src/index.md ##
        ##############################################
        "index.md",

        "faq.md",
        "api.md",
        "contributing.md",
        "citing.md"
    ]
)

deploydocs(
    repo = "github.com/thowell/IterativeLQR.jl.git",
)


text:
"A Expressive Library for Probabilistic Programming"
module Omega

import ForwardDiff
using Spec
using UnicodePlots
using DocStringExtensions: FIELDS, SIGNATURES, TYPEDEF
using Cassette
using Lens
using Reexport

import Random
import Random: GLOBAL_RNG, AbstractRNG
import Statistics: mean, var, quantile

# Util
include("util/Util.jl")
@reexport using .Util

include("ctx/ctx.jl")
@reexport using .Ctx

# Ids
include("ids/ids.jl")
@reexport using .IDS

# Probability Spaces
include("space/space.jl")
using .Space
export Ω, SimpleΩ, LinearΩ, defΩ, defΩProj

# RandVar
include("nondet/nondet.jl" )
using .NonDet
export RandVar, MaybeRV, ciid, isconstant, elemtype, params, constant

# Conditioning
include("cond/cond.jl")
@reexport using .Cond
# export cond

# Lifted random variable operatiosn
include("lift/r.jl")
export ᵣ
include("lift/containers.jl")     # Array/Tuple primitives
export randarray, randtuple, ==ᵣ, tupleᵣ, arrayᵣ

# Lifting functions to RandVar domain
include("lift/lift.jl")
export @lift, lift, lift!

# Soft Inference
include("soft/soft.jl")
@reexport using .Soft

import .Soft: logerr, softeq, softgt, softlt, err, kf1β, kseα

Omega.lift!(:softeq, 2)
Omega.lift!(:softeq, 3)
Omega.lift!(:softgt, 2)
Omega.lift!(:softlt, 2)
Omega.lift!(:logerr, 1)
Omega.lift!(:err, 1)
Omega.lift!(:kf1β, 1)
Omega.lift!(:kseα, 1)

# Higher-Order Inference
include("higher/Higher.jl")
@reexport using .Higher

# Gradient
include("gradient/Gradient.jl")
@reexport using .Gradient

# Inference Algorithms
include("inference/Inference.jl")
@reexport using .Inference

# Causal Inference
include("causal/Causal.jl")
@reexport using .Causal

# Library
include("prim/Prim.jl")
@reexport using .Prim

# Library
include("dist/dist.jl")
@reexport using .Dist

# Neural Network Stuff
include("flux.jl")
@reexport using .OmegaFlux

# Scaling errors
include("scaling.jl")

# The Omega Lanaguage (Not Library)
include("lang/lang.jl")

# Experimental

include("symbolic/symbolic.jl")
@reexport using .Symbolic

end


text:
module tree_number_estimate

function number_estimate(df, lambda)
    fruit_sum = sum(df.Close)
    trip = findmax(df.Trip)[1]
    dfm = df[df.Close.>0, :]
    A = 49*3.16
    AT = area_estimate(dfm)
    println(AT)
    println("lalalal")
    fruit_per_area = lambda/A
    fruit_per_point = fruit_sum/size(dfm)[1]
    tree_number = round(Int, AT/A*fruit_per_point/fruit_per_area)
    return tree_number
end

function area_estimate(df)
    n = size(df)[1]
    len, wid = 108, 108
    arr = rand(len, wid)
    for i in 1:len
        for j in 1:wid
            arr[i, j] = 0
        end
    end
    for i in 1:n
        if df.Close[i]>0
            x = round(Int, df.X[i])
            y = round(Int, df.Y[i])
            arr[x, y] = 1
        end
    end
    area = sum(arr)
    return area/2
end

#test script
#df = CSV.read("H:\\julia_pro\\project\\data_proj_414.csv")
#dfm = df[df.Trip.<10, :]
#println(area_estimate(dfm))
#lambda = 65
#tree_number = number_estimate(df, lambda)
#println(tree_number)

end


text:
module Blender

using PyCall

using GeometryBasics
using CoordinateTransformations

using ColorTypes
using PrettyTables

include("bpy.jl")
include("materials.jl")
include("object.jl")
include("camera.jl")
include("render.jl")

end # module


text:
@testset "model_test_utils" begin
    storages = [
        (name="dense storage", val=ArrayStorage(Float64)),
        (name="static storage", val=SArrayStorage(Float64)),
    ]
    rng = MersenneTwister(123456)
    @testset "storage = $(storage.name)" for storage in storages
        @testset "random_vector" begin
            a = random_vector(rng, 3, storage.val)
            @test is_of_storage_type(a, storage.val)
            @test length(a) == 3
        end
        @testset "random_matrix" begin
            A = random_matrix(rng, 4, 3, storage.val)
            @test is_of_storage_type(A, storage.val)
            @test size(A) == (4, 3)
        end
        @testset "random_nice_psd_matrix" begin
            Σ = random_nice_psd_matrix(rng, 11, storage.val)
            @test all(eigvals(Σ) .> 0)
            @test all(eigvals(Σ) .< 1)
            @test is_of_storage_type(Σ, storage.val)
        end
        @testset "random_gaussian" begin
            x = random_gaussian(rng, 3, storage.val)
            @test is_of_storage_type(x, storage.val)
            @test length(x.m) == 3
            @test size(x.P) == (3, 3)
            @test all(eigvals(x.P) .> 0)
        end
        @testset "GaussMarkovModel" begin
            @testset "time-varying" begin
                x = random_tv_gmm(rng, Forward(), 3, 11, storage.val)
                @test length(x) == 11
                @test dim(x) == 3
            end
            @testset "time-invariant" begin
                x = random_ti_gmm(rng, Forward(), 2, 12, storage.val)
                @test length(x) == 12
                @test dim(x) == 2
            end
        end
        # @testset "LGSSM" begin
        #     @testset "SmallOutput -- time-varying" begin
        #         x = random_tv_lgssm(rng, Forward(), 3, 2, 11, storage.val)
        #         @test length(x) == 11

        #         # Just run `rand_tangent` without checking correctness.
        #         rand_tangent(rng, x)
        #     end
        #     @testset "SmallOutput -- time-invariant" begin
        #         x = random_ti_lgssm(rng, Forward(), 3, 2, 11, storage.val)
        #         @test length(x) == 11
        #         @test x.transitions.As isa Fill
        #         @test x.transitions.as isa Fill
        #         @test x.transitions.Qs isa Fill
        #         @test x.emissions.A isa Fill
        #         @test x.emissions.a isa Fill
        #         @test x.emissions.Q isa Fill

        #         # Just run `rand_tangent` without checking correctness.
        #         rand_tangent(rng, x)
        #     end
        #     @testset "LargeOutput -- time-varying" begin
        #     end
        # end
    end
    @testset "BottleneckLGC - Q=$(Q_type)" for Q_type in [Val(:dense), Val(:diag)]
        x = random_bottleneck_lgc(rng, 5, 2, 4, Q_type, ArrayStorage(Float64))
        @test size(x.H) == (2, 5)
        @test size(x.h) == (2,)
        @test eltype(x) == Float64
    end
end


text:
using KUparser, KUnet, HDF5, JLD, ArgParse, Compat, CUDArt
VERSION < v"0.4-" && eval(Expr(:using,:Dates))

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "datafiles"
        help = "JLD files for input: first one will be used for training"
        nargs = '+'
        required = true
        "--out"
        help = "JLD file for output trained network"
        default = ""
        "--net"
        help = "Initialize net from this file if specified"
        default = ""
        "--machinefile"
        help = "Need this to restart workers to avoid pmap leak"
        default = ""
        # "--nogpu"  # always use gpu
        # help = "Do not use gpu"
        # action = :store_true
        # "--parser"  # always use bparser
        # help = "Parsing algorithm to use: g(reedy)parser, b(eam)parser, or o(racle)parser for static training"
        # default = "oparser"
        "--arctype"
        help = "Move set to use: ArcEager{R1,13}, ArcHybrid{R1,13}"
        default = "ArcEager13"
        "--feats"
        help = "Features to use from KUparser.Flist"
        default = "tacl13hybrid"
        "--hidden"
        help = "One or more hidden layer sizes"
        nargs = '+'
        arg_type = Int
        default = [20000]
        "--epochs"
        help = "Minimum number of epochs to train"
        arg_type = Int
        default = 100
        # "--pepochs"
        # help = "Number of epochs between parsing"
        # arg_type = Int
        # default = 1
        # "--ncpu"
        # help = "Number of workers for multithreaded parsing"
        # arg_type = Int
        # default = 12
        "--pbatch"
        help = "Minibatch size for parsing"
        arg_type = Int
        default = 100
        "--tbatch"
        help = "Minibatch size for training"
        arg_type = Int
        default = 128
        "--nbeam"
        help = "Beam size for bparser"
        arg_type = Int
        default = 64
        "--adagrad"
        help = "If nonzero apply adagrad using arg as epsilon"
        arg_type = Float32
        nargs = '+'
        default = [1f-8]
        "--dropout"
        help = "Dropout probability"
        arg_type = Float32
        nargs = '+'
        default = [0.2f0, 0.7f0]
        "--l1reg"
        help = "L1 regularization parameter"
        arg_type = Float32
        nargs = '+'
        "--l2reg"
        help = "L2 regularization parameter"
        arg_type = Float32
        nargs = '+'
        "--learningRate"
        help = "Learning rate"
        arg_type = Float32
        nargs = '+'
        default = [2f-2]
        "--maxnorm"
        help = "If nonzero upper limit on weight matrix row norms"
        arg_type = Float32
        nargs = '+'
        "--momentum"
        help = "Momentum"
        arg_type = Float32
        nargs = '+'
        "--nesterov"
        help = "Apply nesterov's accelerated gradient"
        arg_type = Float32
        nargs = '+'
        "--seed"
        help = "Random number seed (-1 leaves the default seed)"
        arg_type = Int
        default = -1
    end
    args = parse_args(s)
    for (arg,val) in args
        print("$arg:$val ")
    end
    println("")
    return args
end

function pbatch(c::Corpus, n::Int)
    dist = Corpus[]
    for s1=1:n:length(c)
        s2=min(length(c), s1+n-1)
        push!(dist, c[s1:s2])
    end
    return dist
end

function evalpxy(pxy, data)
    p = map(a->a[1], pxy)
    evalp(p, data)
end

function evalp(p, d)
    pa = vcat(p...)
    da = vcat(d...)
    evalparse(pa, da)
end

function initdata(args)
    pt = eval(parse(args["arctype"]))
    data = Any[]
    ndeps = 0
    feats = eval(parse("Flist."*args["feats"]))
    net = nothing
    deprel = nothing
    for f in args["datafiles"]
        @date d=load(f)
        push!(data, pbatch(d["corpus"], args["pbatch"]))
        if deprel == nothing
            deprel = d["deprel"]
            ndeps = length(deprel)
        else
            @assert deprel == d["deprel"]
        end
    end
    if isempty(args["net"])
        s1 = data[1][1][1]
        p1 = pt(wcnt(s1),length(deprel))
        xrows=KUparser.flen(p1, s1, feats)
        yrows=p1.nmove
        net=newnet(relu, [xrows; args["hidden"]; yrows]...)
        net[end].f=KUnet.logp
        for k in [fieldnames(UpdateParam); :dropout]
            haskey(args, string(k)) || continue
            v = args[string(k)]
            if isempty(v)
                continue
            elseif length(v)==1
                setparam!(net, k, v[1])
            else
                @assert length(v)==length(net) "$k should have 1 or $(length(net)) elements"
                for i=1:length(v)
                    setparam!(net[i], k, v[i])
                end
            end
        end
        # Initialize net using oparse on first corpus
        @date (p,x,y) = oparse(pt, data[1], ndeps, feats)
        @date train(net, x, y; batch=args["tbatch"], loss=KUnet.logploss, shuffle=true)
    else
        net = newnet(args["net"])
    end
    return (pt, data, ndeps, feats, net)
end

function restart_all_workers(wlist)
    @date Base.terminate_all_workers()
    @date gc()
    @show Base.PGRP.refs
    @date addprocs(wlist)
    @date @everywhere eval(Expr(:using,:KUnet))
    @date @everywhere eval(Expr(:using,:KUparser))
end

function myparse{T<:Parser}(pt::Type{T}, ca::AbstractArray, ndeps::Integer, feats::Fvec, net::Net, args::Dict; trn::Bool=false)
    @date tnet = testnet(net)
    if trn
        @date pxy = pmap(ca) do c
            bparse(pt, c, ndeps, feats, copy(tnet,:gpu), args["nbeam"], args["pbatch"]; xy=true)
        end
        @show e = evalpxy(pxy, ca); flush(STDOUT)
        @date for (p,x,y) in pxy
            train(net, x, y; batch=args["tbatch"], loss=KUnet.logploss, shuffle=true)
        end
    else
        @date p = pmap(ca) do c
            bparse(pt, c, ndeps, feats, copy(tnet,:gpu), args["nbeam"], args["pbatch"])
        end
        @show e = evalp(p, ca); flush(STDOUT)
    end
    return e[1]
end

function rss()
    statm = split(readall("/proc/$(getpid())/statm"))
    pages = int(statm[2])
    mb = pages >> 8
end

function main()
    args = parse_commandline()
    args["seed"] >= 0 && (srand(args["seed"]); KUnet.gpuseed(args["seed"]))
    @date (pt, data, ndeps, feats, net) = initdata(args)
    @show (pt, ndeps, feats)
    @show map(size, data)
    @show net
    accuracy = Array(Float32, length(data))
    (bestscore,bestepoch,epoch)=(0,0,0)
    @show wlist = (isempty(args["machinefile"]) ? [] : split(readall(args["machinefile"])))

    while true
        @show epoch += 1; flush(STDOUT)
        for i=1:length(data)
            @show rss(); flush(STDOUT)
            @date isempty(wlist) || restart_all_workers(wlist)
            @show rss(); flush(STDOUT)
            @date gc()
            @show rss(); flush(STDOUT)
            @date accuracy[i] = myparse(pt, data[i], ndeps, feats, net, args; trn=(i==1))
            @show rss(); flush(STDOUT)
            @date gc()
            @show rss(); flush(STDOUT)
        end
        println("DATA:\t$epoch\t"*join(accuracy, '\t')); flush(STDOUT)

        # We look at accuracy[2] (dev score) by default,
        # accuracy[1] (trn score) if there is no dev.

        score = (length(accuracy) > 1 ? accuracy[2] : accuracy[1])

        # Update best score and save best net

        if (score > bestscore)
            @show (bestscore,bestepoch)=(score,epoch); flush(STDOUT);
            !isempty(args["out"]) && save(args["out"], net)
        end

        # Quit if we have the minimum number of epochs and have
        # not made progress in the last half.

        (epoch >= args["epochs"]) && (epoch >= 2*bestepoch) && break

    end # while true
end # main

main()


text:
# Trait to determine if an Index, Index collection, Tensor, or ITensor
# has symmetries
abstract type SymmetryStyle end

function symmetrystyle(T)
  return error("No SymmetryStyle defined for the specified object $T of type $(typeof(T))")
end

symmetrystyle(T, S, U, V...)::SymmetryStyle = (
  Base.@_inline_meta; symmetrystyle(symmetrystyle(T), symmetrystyle(S, U, V...))
)

symmetrystyle(T, S)::SymmetryStyle = symmetrystyle(symmetrystyle(T), symmetrystyle(S))

# Rules for basic collections
symmetrystyle(inds::Tuple) = symmetrystyle(inds...)
# `reduce(symmetrystyle, inds)` is not type stable for some reason
symmetrystyle(inds::AbstractVector) = mapreduce(symmetrystyle, symmetrystyle, inds)


text:
#
# oiformat1.jl --
#
# Define 1st revision of OI-FITS format.
#
#------------------------------------------------------------------------------
#
# This file is part of OIFITS.jl which is licensed under the MIT "Expat"
# License:
#
# Copyright (C) 2015-2019: Éric Thiébaut.
#
#------------------------------------------------------------------------------

# OI-FITS FORMAT DESCRIPTION TABLES
#
# The format of the OI-FITS data block is described by a vector of strings
# like:
#
#   ["KEYWORD FORMAT DESCR",
#     ...,
#     ...,
#    "---------------------------",
#    "COLUMN  FORMAT DESCR",
#     ...,
#     ...,
#     ...]
#
# where:
#
#   KEYWORD = keyword for HDU header.
#   COLUMN = column name for table (TTYPE).
#   FORMAT = (may be prefixed with a ? to indicate optional field):
#       for keywords, a single letter indicating the type,
#       for columns, a letter indicating the type followed by list of
#                    dimensions in parenthesis (letter W for a dimension means
#                    NWAVE).
#   DESCR = description/comment; units, if any are indicated at the end
#           between square brackets.
#
# There may be any number of keyword definitions and any number of column
# definitions, the two parts are separated by a dash line like
# "--------------".
#

# OI_TARGET definition (1st revision):
add_def("OI_TARGET", 1,
        ["OI_REVN    I      revision number of the table definition",
         "---------------------------------------------------------",
         "TARGET_ID  I(1)   index number",
         "TARGET     A(16)  target name",
         "RAEP0      D(1)   RA at mean equinox [deg]",
         "DECEP0     D(1)   DEC at mean equinox [deg]",
         "EQUINOX    E(1)   equinox [yr]",
         "RA_ERR     D(1)   error in RA at mean equinox [deg]",
         "DEC_ERR    D(1)   error in DEC at mean equino [deg]",
         "SYSVEL     D(1)   systemic radial velocity [m/s]",
         "VELTYP     A(8)   reference for radial velocity",
         "VELDEF     A(8)   definition of radial velocity",
         "PMRA       D(1)   proper motion in RA [deg/yr]",
         "PMDEC      D(1)   proper motion in DEC [deg/yr]",
         "PMRA_ERR   D(1)   error of proper motion in RA [deg/yr]",
         "PMDEC_ERR  D(1)   error of proper motion in DEC [deg/yr]",
         "PARALLAX   E(1)   parallax [deg]",
         "PARA_ERR   E(1)   error in parallax [deg]",
         "SPECTYP    A(16)  spectral type"])

# OI_ARRAY definition (1st revision):
add_def("OI_ARRAY", 1,
        ["OI_REVN     I     revision number of the table definition",
         "ARRNAME     A     array name for cross-referencing",
         "FRAME       A     coordinate frame",
         "ARRAYX      D     array center X-coordinate [m]",
         "ARRAYY      D     array center Y-coordinate [m]",
         "ARRAYZ      D     array center Z-coordinate [m]",
         "------------------------------------------------------------",
         "TEL_NAME   A(16)  telescope name",
         "STA_NAME   A(16)  station name",
         "STA_INDEX  I(1)   station index",
         "DIAMETER   E(1)   element diameter [m]",
         "STAXYZ     D(3)   station coordinates relative to array center [m]"])

# OI_WAVELENGTH definition (1st revision):
add_def("OI_WAVELENGTH", 1,
        ["OI_REVN     I     revision number of the table definition",
         "INSNAME     A     name of detector for cross-referencing",
         "------------------------------------------------------------",
         "EFF_WAVE   E(1)   effective wavelength of channel [m]",
         "EFF_BAND   E(1)   effective bandpass of channel [m]"])

# OI_VIS definition (1st revision):
add_def("OI_VIS", 1,
        ["OI_REVN     I    revision number of the table definition",
         "DATE-OBS    A    UTC start date of observations",
         "ARRNAME    ?A    name of corresponding array",
         "INSNAME     A    name of corresponding detector",
         "------------------------------------------------------------",
         "TARGET_ID  I(1)  target number as index into OI_TARGET table",
         "TIME       D(1)  UTC time of observation [s]",
         "MJD        D(1)  modified Julian Day [day]",
         "INT_TIME   D(1)  integration time [s]",
         "VISAMP     D(W)  visibility amplitude",
         "VISAMPERR  D(W)  error in visibility amplitude",
         "VISPHI     D(W)  visibility phase [deg]",
         "VISPHIERR  D(W)  error in visibility phase [deg]",
         "UCOORD     D(1)  U coordinate of the data [m]",
         "VCOORD     D(1)  V coordinate of the data [m]",
         "STA_INDEX  I(2)  station numbers contributing to the data",
         "FLAG       L(W)  flag"])

# OI_VIS2 definition (1st revision):
add_def("OI_VIS2", 1,
        ["OI_REVN     I    revision number of the table definition",
         "DATE-OBS    A    UTC start date of observations",
         "ARRNAME    ?A    name of corresponding array",
         "INSNAME     A    name of corresponding detector",
         "------------------------------------------------------------",
         "TARGET_ID  I(1)  target number as index into OI_TARGET table",
         "TIME       D(1)  UTC time of observation [s]",
         "MJD        D(1)  modified Julian Day [day]",
         "INT_TIME   D(1)  integration time [s]",
         "VIS2DATA   D(W)  squared visibility",
         "VIS2ERR    D(W)  error in squared visibility",
         "UCOORD     D(1)  U coordinate of the data [m]",
         "VCOORD     D(1)  V coordinate of the data [m]",
         "STA_INDEX  I(2)  station numbers contributing to the data",
         "FLAG       L(W)  flag"])

# OI_T3 definition (1st revision):
add_def("OI_T3", 1,
        ["OI_REVN     I    revision number of the table definition",
         "DATE-OBS    A    UTC start date of observations",
         "ARRNAME    ?A    name of corresponding array",
         "INSNAME     A    name of corresponding detector",
         "------------------------------------------------------------",
         "TARGET_ID  I(1)  target number as index into OI_TARGET table",
         "TIME       D(1)  UTC time of observation [s]",
         "MJD        D(1)  modified Julian Day [day]",
         "INT_TIME   D(1)  integration time [s]",
         "T3AMP      D(W)  triple product amplitude",
         "T3AMPERR   D(W)  error in triple product amplitude",
         "T3PHI      D(W)  triple product phase [deg]",
         "T3PHIERR   D(W)  error in triple product phase [deg]",
         "U1COORD    D(1)  U coordinate of baseline AB of the triangle [m]",
         "V1COORD    D(1)  V coordinate of baseline AB of the triangle [m]",
         "U2COORD    D(1)  U coordinate of baseline BC of the triangle [m]",
         "V2COORD    D(1)  V coordinate of baseline BC of the triangle [m]",
         "STA_INDEX  I(3)  station numbers contributing to the data",
         "FLAG       L(W)  flag"])


text:
# This method is called when a SemidiscretizationHyperbolic is constructed.
# It constructs the basic `cache` used throughout the simulation to compute the RHS etc.
# Additionally the cache can contain the gradient matrices for Euler equation with viscous terms

function create_cache(mesh::CurvedMesh, equations::AbstractEquations, dg::DG, ::Any, ::Type{uEltype}) where {uEltype <: Real}
    elements = init_elements(mesh, equations, dg.basis, uEltype)

    if typeof(equations) == CompressibleEulerEquations2D{Float64} && equations.viscous == true
        q1 = zeros(nvariables(equations), length(dg.basis.nodes), length(dg.basis.nodes), size(mesh, 1) * size(mesh, 2))
        q2 = zeros(nvariables(equations), length(dg.basis.nodes), length(dg.basis.nodes), size(mesh, 1) * size(mesh, 2))
        cache = (; elements, q1, q2)
    else
        cache = (; elements)
    end

    return cache
end

# Extract contravariant vector Ja^i (i = index) as SVector
@inline function get_contravariant_vector(index, contravariant_vectors, indices...)

    SVector(ntuple(dim -> contravariant_vectors[dim, index, indices...], ndims(contravariant_vectors) - 3))
end

@inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, orientation,
                                                  boundary_condition::BoundaryConditionPeriodic,
                                                  mesh::CurvedMesh, equations, dg::DG, cache,
                                                  direction, node_indices, surface_node_indices, element)
    @assert isperiodic(mesh, orientation)
end
@inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, orientation,
                                                  boundary_condition::BoundaryConditionPeriodic,
                                                  mesh::CurvedMesh, equations::AuxiliaryEquation, dg::DG, cache,
                                                  direction, node_indices, surface_node_indices, element, nabla)
    @assert isperiodic(mesh, orientation)
end

# BLZ: ersetzt und in dg_2d gepackt
# @inline function calc_boundary_flux_by_direction!(surface_flux_values, u, t, orientation,
#                                                   boundary_condition,
#                                                   mesh::CurvedMesh, equations,dg::DG, cache,
#                                                   direction, node_indices, surface_node_indices, element)
#   @unpack node_coordinates, contravariant_vectors = cache.elements
#   @unpack surface_flux = dg

#   u_inner = get_node_vars(u, equations, dg, node_indices..., element)
#   x = get_node_coords(node_coordinates, equations, dg, node_indices..., element)

#   # Contravariant vector Ja^i is the normal vector
#   normal = get_contravariant_vector(orientation, contravariant_vectors, node_indices..., element)
#   flux = boundary_condition(u_inner, normal, direction, x, t, surface_flux, equations)

#   for v in eachvariable(equations)
#     surface_flux_values[v, surface_node_indices..., direction, element] = flux[v]
#   end
# end

@inline ndofs(mesh::CurvedMesh, dg::DG, cache) = nelements(cache.elements) * nnodes(dg)^ndims(mesh)

include("containers.jl")
include("dg_1d.jl")
include("dg_2d.jl")
include("dg_3d.jl")


text:
using Test
using LogicCircuits

@testset "Trimmed apply test" begin

    num_vars = 7
    mgr = SddMgr(num_vars, :balanced)

    x = Var(1)
    y = Var(2)

    x_c = compile(mgr, var2lit(x))
    y_c = compile(mgr, var2lit(y))

    notx = -var2lit(x)
    notx_c = compile(mgr,notx)
    true_c = compile(mgr,true)
    @test true_c isa Sdd
    @test true_c === mgr(true)
    @test true_c === compile(Sdd,mgr,true)
    @test true_c === (Sdd,mgr)(true)
    @test true_c === (mgr,Sdd)(true)
    false_c = compile(mgr,false)
    @test false_c === mgr(false)

    @test false_c & true_c == false_c
    @test false_c & notx_c == false_c
    @test false_c & x_c == false_c
    @test true_c & notx_c == notx_c
    @test x_c & notx_c == false_c
    @test true_c & true_c == true_c
    @test false_c & false_c == false_c
    @test x_c & x_c == x_c
    @test !x_c & !x_c == !x_c

    @test false_c | true_c == true_c
    @test false_c | notx_c == notx_c
    @test false_c | x_c == x_c
    @test true_c | notx_c == true_c
    @test x_c | notx_c == true_c
    @test true_c | true_c == true_c
    @test false_c | false_c == false_c
    @test x_c | x_c == x_c
    @test !x_c | !x_c == !x_c

    v1 = compile(mgr, Lit(1))
    v1_2 = compile(Sdd, mgr, Lit(1))
    @test v1 === v1_2

    v2 = compile(mgr, Lit(2))
    v3 = compile(mgr, Lit(3))
    v4 = compile(mgr, Lit(4))
    v5 = compile(mgr, Lit(5))
    v6 = compile(mgr, Lit(6))
    v7 = compile(mgr, Lit(7))

    t1 = v1 & v3
    t2 = v3 & v1

    @test t1 === t2
    @test model_count(t1,num_vars) == BigInt(2)^(num_vars-2)

    @test model_count(v1 & v2 & v6,num_vars) == BigInt(2)^(num_vars-3)

    c1 = v1 | v3
    c2 = v3 | v1

    @test c1 === c2
    @test model_count(c1,num_vars) == BigInt(2)^(num_vars-2) * 3

    f1 = (c1 & c2)

    @test f1 === (c2 & c1)
    @test f1 === (c1 & c2 & c2)
    @test f1 === (c1 & c2 & c2 & true_c)

    @test (v3 | v7) !== false_c
    @test (v3 | v7) & (!v3 | v7) !== false_c
    @test (v3 | v7) & (!v3 | v7) & (v3 | !v7) !== false_c
    @test (v3 | v7) & (!v3 | v7) & (v3 | !v7) & (!v3 | !v7) === false_c

    f2 = (c1 | c2)

    @test f2 === (c2 | c1)
    @test f2 === (c1 | c2 | c2)
    @test f2 === (c1 | c2 | c2 | false_c)

    @test (v3 & v7) !== true_c
    @test (v3 & v7) | (!v3 & v7) !== true_c
    @test (v3 & v7) | (!v3 & v7) | (v3 & !v7) !== true_c
    @test (v3 & v7) | (!v3 & v7) | (v3 & !v7) | (!v3 & !v7) === true_c

    @test (v1 & v7 & v3) === !(!v1 | !v7 | !v3)

    @test f2 & !f2 === false_c
    @test f2 | !f2 === true_c

    f3 = (v1 | !v3 | v5) & (!v2 | !v3 | v7) & (v2 | !v6 | !v5) & (v1 | v2 | v3) & (!v5 | v6 | v7) & (v2 | !v3 | !v6) & (!v1 | !v6 | !v7)

    @test f3 & !f3 === false_c
    @test f3 | !f3 === true_c
    @test f3 & !(v2 | !v6 | !v5) === false_c

    f4 = ((v1 | !v3 | v5) & (!v5 | v6 | v7| v1) & (!v5 | v6 | v7)  & (!v1 | !v6 | !v7)  & (!v2 | !v3 | v7)) & ((v2 | !v3 | !v6)  & (v2 | !v6 | !v5) & (v1 | v2 | v3) & (v1 | v2 | v3 | !v6) )

    @test f3 === f4

    @test isdeterministic(f4)
    @test respects_vtree(f4)

    io = IOBuffer()
    show(io,f3)
    show(io,t1)
    show(io,c1)
    show(io,true_c)
    show(io,false_c)
    show(io,v1)
    @test length(String(take!(io))) > 0 # see if it runs, regardless of result

    f4c = LogicCircuit(f4)
    @test f4c isa PlainLogicCircuit
    @test num_edges(f4c) == num_edges(f4)
    @test num_nodes(f4c) == num_nodes(f4)
    @test model_count(f4c) == model_count(f4)
    @test isdeterministic(f4c)

    f4c = StructLogicCircuit(mgr, f4)
    @test f4c isa PlainStructLogicCircuit
    @test num_edges(f4c) == num_edges(f4)
    @test num_nodes(f4c) == num_nodes(f4)
    @test model_count(f4c) == model_count(f4)
    @test f4c.vtree === f4.vtree
    @test isdeterministic(f4c)

    f4c = StructLogicCircuit(mgr, LogicCircuit(f4))
    @test f4c isa PlainStructLogicCircuit
    @test num_edges(f4c) == num_edges(f4)
    @test num_nodes(f4c) == num_nodes(f4)
    @test model_count(f4c) == model_count(f4)
    @test f4c.vtree === f4.vtree

    f4c = mgr(StructLogicCircuit(mgr,f4))
    @test Sdd(mgr,f4c) === f4c
    @test f4 === f4c
    @test f4c isa Sdd
    @test num_edges(f4c) == num_edges(f4)
    @test num_nodes(f4c) == num_nodes(f4)
    @test model_count(f4c) == model_count(f4)
    @test f4c.vtree === f4.vtree

    # @test_throws Exception propagate_constants(f4) # breaks bit in canonical false node, creating problems in other tests

end

@testset "Trimmed apply regression test 1" begin

    v = zoo_vtree("iscas89/s386.scan.min.vtree")
    mgr = SddMgr(v)
    v70 = compile(mgr,Lit(70))
    v71 = compile(mgr,Lit(71))
    n = (v70 & v71) | (!v70 & !v71)
    n & v70
    !n & !v70
    r = (!n & v70)
    @test "$r" == "[(70,-71),(-70,⊥)]" || "$r" == "[(-70,⊥),(70,-71)]"

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 98e662cc-8407-11ec-29f9-1dcf48b4d326
import PlutoHooks, Revise

# ╔═╡ db73b8ec-893b-4a80-8170-4223c7164e78
# PlutoLinks is defined in terms of itself
import PlutoLinks

# ╔═╡ 5d3b18ce-e0e9-4fec-b370-1a80672ea800
"""

The `@revise` macro imports a package watched by Revise. The cell from which `@revise` is called will be re-run whenever the package definition changes. This means that cells depending on this package will be updated as well.

 - The revised package should be developed in the current environment, usually requiring you to [opt out of the Pluto package manager](https://github.com/fonsp/Pluto.jl/wiki/%F0%9F%8E%81-Package-management#advanced-set-up-an-environment-with-pkgactivate).

---
Examples

```julia
# in one cell
using PlutoLinks

# in another cell
@revise using MyPackage

# import is also supported
@revise import MyPackage
```
"""
macro revise(using_)
    @assert(
		(Meta.isexpr(using_, :using, 1) || Meta.isexpr(using_, :import, 1)) &&
			Meta.isexpr(using_.args[1], :(.), 1),
		"Expected expression of from :(using Package) but got $(repr(using_))"
	)

    mod = last(using_.args[1].args)
    quote
        $(esc(using_))

		if @PlutoHooks.use_is_pluto_cell()
			# let block might be important to prevent Pluto from setting the callback
			# to nothing when jumping to a new workspace
	        let
	            _, rerun = @PlutoHooks.use_state(nothing)
	            @PlutoLinks.use_task([]) do
	                Revise.entr(() -> rerun(nothing),
						[], [$(esc(mod))]; all=false, postpone=true)
	            end
	        end
		end

		nothing
    end
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoHooks = "0ff47ea0-7a50-410d-8455-4348d5de0774"
PlutoLinks = "0ff47ea0-7a50-410d-8455-4348d5de0420"
Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"

[compat]
PlutoHooks = "~0.0.3"
PlutoLinks = "~0.1.2"
Revise = "~3.3.1"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.2"
manifest_format = "2.0"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.CodeTracking]]
deps = ["InteractiveUtils", "UUIDs"]
git-tree-sha1 = "9aa8a5ebb6b5bf469a7e0e2b5202cf6f8c291104"
uuid = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
version = "1.0.6"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.FileWatching]]
uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.JuliaInterpreter]]
deps = ["CodeTracking", "InteractiveUtils", "Random", "UUIDs"]
git-tree-sha1 = "6ca01d8e5bc75d178e8ac2d1f741d02946dc1853"
uuid = "aa1ae85d-cabe-5617-a682-6adf51b2e16a"
version = "0.9.2"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.LoweredCodeUtils]]
deps = ["JuliaInterpreter"]
git-tree-sha1 = "f46e8f4e38882b32dcc11c8d31c131d556063f39"
uuid = "6f1432cf-f94c-5a45-995e-cdbf5db27b0b"
version = "2.2.0"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PlutoHooks]]
deps = ["FileWatching", "InteractiveUtils", "Markdown", "UUIDs"]
git-tree-sha1 = "f297787f7d7507dada25f6769fe3f08f6b9b8b12"
uuid = "0ff47ea0-7a50-410d-8455-4348d5de0774"
version = "0.0.3"

[[deps.PlutoLinks]]
deps = ["FileWatching", "InteractiveUtils", "Markdown", "PlutoHooks", "Revise", "UUIDs"]
git-tree-sha1 = "c9d92822ea27c1453d3615bae44481859296fe31"
uuid = "0ff47ea0-7a50-410d-8455-4348d5de0420"
version = "0.1.2"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.3.0"

[[deps.Revise]]
deps = ["CodeTracking", "Distributed", "FileWatching", "JuliaInterpreter", "LibGit2", "LoweredCodeUtils", "OrderedCollections", "Pkg", "REPL", "Requires", "UUIDs", "Unicode"]
git-tree-sha1 = "2f9d4d6679b5f0394c52731db3794166f49d5131"
uuid = "295af30f-e4ad-537b-8983-00126c2a3abe"
version = "3.3.1"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╠═98e662cc-8407-11ec-29f9-1dcf48b4d326
# ╠═db73b8ec-893b-4a80-8170-4223c7164e78
# ╠═5d3b18ce-e0e9-4fec-b370-1a80672ea800
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
# ---
# title: 226. Invert Binary Tree
# id: problem226
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Tree
# link: <https://leetcode.com/problems/invert-binary-tree/description/>
# hidden: true
# ---
#
# Invert a binary tree.
#
# **Example:**
#
# Input:
#
#
#
#          4
#        /   \
#       2     7
#      / \   / \
#     1   3 6   9
#
# Output:
#
#
#
#          4
#        /   \
#       7     2
#      / \   / \
#     9   6 3   1
#
# **Trivia:**
# This problem was inspired by [this original
# tweet](https://twitter.com/mxcl/status/608682016205344768) by [Max
# Howell](https://twitter.com/mxcl):
#
# > Google: 90% of our engineers use the software you wrote (Homebrew), but you
# > can't invert a binary tree on a whiteboard so f*** off.
#
#
## @lc code=start
using LeetCode

function invert_tree!(root::Union{TreeNode{Int}, Nothing})::Union{TreeNode, Nothing}
    isnothing(root) && return nothing
    queue = [root]
    while !isempty(queue)
        node = popfirst!(queue)
        if !isnothing(node.left)
            push!(queue, node.left)
        end
        if !isnothing(node.right)
            push!(queue, node.right)
        end
        node.left, node.right = node.right, node.left
    end
    return root
end
## @lc code=end

## @lc test=start
@testset "226.invert-binary-tree.jl" begin
    @test invert_tree!(TreeNode{Int}([4,2,7,1,3,6,9])) == TreeNode{Int}([4,7,2,9,6,3,1])
    @test invert_tree!(TreeNode{Int}([2,1,3])) == TreeNode{Int}([2,3,1])
    @test invert_tree!(nothing) === nothing
end
## @lc test=end

## @lc info=start
# link: [solution to 226](https://leetcode-cn.com/problems/invert-binary-tree/solution/226-fan-zhuan-er-cha-shu-python-julia-by-sbu2/)
## @lc info=end


text:
module UtilTest

using Test
using Statistics, LinearAlgebra
import KLIEPInference: trimap, itrimap, pack, unpack, Ψising, rhat, KLIEP_Hessian, _find_supp
import KLIEPInference: chain, IsingSampler
import Random: rand

@testset "trimap" begin
    @test trimap(1, 2) == 1
	@test trimap(3, 4) == 6

	@test itrimap(1) == CartesianIndex(1, 2)
	@test itrimap(6) == CartesianIndex(3, 4)
end

@testset "pack" begin
	A = [1. 2.; 2. 3.]
	θ = [2.]
	@test pack(A) == θ
	@test unpack(θ) == [0. 2.; 2. 0]

	A = [1. 2. 3.; 2. 2. 0.5; 3. 0.5 3.]
	θ = [2., 3., 0.5]
	@test pack(A) == θ
	@test unpack(θ) == [0. 2. 3.; 2. 0. 0.5; 3. 0.5 0.]
end

@testset "Ising" begin
	X = Array{Bool}([1 1; 0 1; 0 0; 1 0])
	Ψ = [-1. 1; -1 -1; 1 -1; 1 -1; -1 -1; -1 1]
	@test (@test_logs (:warn, "a row of Ψ has zero variance") match_mode=:any Ψising(X)) == Ψ

	for rep=1:50
		n = 200
		p = 3
		m = 3
		θ = chain(m; lenC=3, lb=0.1, ub=0.2, sgn=0)
		spl = IsingSampler(θ; burn=5000, thin=5000)
		X = rand(spl, n)
		Ψ = Ψising(X)

		# check rhat
		_rhat = zeros(n)
		_mr = 0.
		for i = 1:n
			_rhat[i] = exp( dot(Ψ[:, i], θ) )
			_mr += _rhat[i]
		end
		_mr /= n
		for i = 1:n
			_rhat[i] /= _mr
		end
		@test rhat(θ, Ψ) ≈ _rhat atol=1e-5

		# check Hessian
		_H = zeros(3, 3)
		_hμ = zeros(3)
		for j=1:p
			for i=1:n
				_hμ[j] += Ψ[j, i] * _rhat[i]
			end
			_hμ[j] /= n
		end
		for j=1:p
			for k=1:p
				for i=1:n
					_H[j, k] += Ψ[j, i] * Ψ[k, i] * _rhat[i]
				end
				_H[j, k] /= n
			end
		end
		_H -= _hμ * transpose(_hμ)
		@test KLIEP_Hessian(θ, Ψ) ≈ _H atol=1e-5
	end
end

@testset "_find_supp" begin
	x = [0., 1., 0., 3., 0.]
	@test _find_supp(1, x) == [2, 4, 1]
	@test _find_supp(2, x) == [4, 2]
	@test _find_supp(3, x) == [2, 4, 3]
	@test _find_supp(4, x) == [2, 4]
	@test _find_supp(5, x) == [2, 4, 5]

	x = [0., 0., 0., 3., 0.]
	y = [0., 1., 1., 3., 0.]
	@test _find_supp(1, x, y) == [4, 2, 3, 1]
	@test _find_supp(2, x, y) == [4, 3, 2]
	@test _find_supp(3, x, y) == [4, 2, 3]
	@test _find_supp(4, x, y) == [3, 2, 4]
	@test _find_supp(5, x, y) == [4, 2, 3, 5]
end

end


text:
import Tables
import Base: show, setproperty!
using PrettyTables

"""
    CleanTable <: Tables.AbstractColumns

A Tables.jl implementation that stores column names and columns for Cleaner.jl internal use.

The default behavior of this type is to try to copy the columns of the original Tables
implementation a.k.a: the source, but the user can call the second constructor specifiying
`copycols=false` to override this behavior and try to use the original columns directly.
If the source column type is not mutable, it might end up in errors.

# Constructors
```julia
CleanTable(names::Vector{Symbol}, cols; copycols::Bool=false)
CleanTable(table; copycols::Bool=true)
CleanTable(table::CleanTable; copycols::Bool=true)
```
"""
mutable struct CleanTable <: Tables.AbstractColumns
    names::Vector{Symbol}
    cols::Vector{AbstractVector}

    function CleanTable(names::Vector{Symbol}, cols; copycols::Bool=false)
        length(names) != length(cols) && error("Inconsistent length between names given and
        amount of columns")

        nrow = length(cols[1])
        for col in cols
            nrow != length(col) && error("All columns must be of the same length")
        end

        if copycols
            if Threads.nthreads() > 1 && length(names) > 1 && nrow >= 1_000_000
                copied_cols = Vector{AbstractVector}(undef, length(names))

                Threads.@threads for i in 1:length(cols)
                    copied_cols[i] = copy(cols[i])
                end

                return new(copy(names), copied_cols)
            end

            return new(copy(names), [copy(col) for col in cols])
        else
            return new(names, cols)
        end
    end
end

_getvector(x::AbstractVector) = x
_getvector(x) = collect(x)

function CleanTable(table; copycols::Bool=true)
    columns = Tables.columns(table)

    names = [Symbol(name) for name in Tables.columnnames(columns)]
    cols = Vector[_getvector(Tables.getcolumn(columns, name)) for name in names]

    return CleanTable(names, cols; copycols=copycols)
end

function CleanTable(table::CleanTable; copycols::Bool=true)
    return CleanTable(names(table), cols(table); copycols=copycols)
end

Tables.istable(::Type{<:CleanTable}) = true

names(ct::CleanTable) = getfield(ct, :names)
cols(ct::CleanTable) = getfield(ct, :cols)

Tables.schema(ct::CleanTable) = Tables.Schema(names(ct), [eltype(col) for col in cols(ct)])

Tables.columnaccess(::Type{<:CleanTable}) = true
Tables.columns(ct::CleanTable) = ct
Tables.getcolumn(ct::CleanTable, i::Int) = getindex(cols(ct), i)
function Tables.getcolumn(ct::CleanTable, name::Symbol)
    return getindex(cols(ct), findfirst(isequal(name), names(ct)))
end
function Tables.getcolumn(ct::CleanTable, ::Type{T}, i::Int, nm::Symbol) where {T}
    return getindex(cols(ct), i)
end
Tables.columnnames(ct::CleanTable) = names(ct)

Tables.materializer(::CleanTable) = CleanTable

Base.show(io::IO, ct::CleanTable) = pretty_table(io, ct)

"""
    size(table::CleanTable)

Returns a tuple containing the number of rows and columns of the given `CleanTable`.
"""
function size(table::CleanTable)
    return (length(cols(table)[1]), length(names(table)))
end

"""
    rename!(ct::CleanTable, names::Vector{Symbol})

Changes in-place the column names of a `CleanTable` to be `names`.
"""
function rename!(ct::CleanTable, names::Vector{Symbol})
    if length(cols(ct)) != length(names)
        error("Inconsistent length between names given and amount of columns")
    end

    setfield!(ct, :names, names)

    return ct
end

"""
    rename(table, names::Vector{Symbol})

Creates a `CleanTable` with copied columns and changes its column names to be `names`.
"""
function rename(table, names::Vector{Symbol})
    return rename!(CleanTable(table), names)
end

function Base.setproperty!(ct::CleanTable, s::Symbol, x::AbstractVector)
    if s in names(ct)
        if length(cols(ct)[1]) == length(x) || length(cols(ct)) == 1
            index = findfirst(x -> isequal(x, s), names(ct))
            cols(ct)[index] = x

            return nothing
        else
            error(
                "Length of the AbstractVector being assigned to column '$s' must be consistent with the length of the rest of columns",
            )
        end
    elseif s == :names
        error(
            "Changing property 'names' directly is unsupported, please use the `rename!` function",
        )
    elseif s == :cols
        error("Property 'cols' cannot be changed")
    else
        error("Property '$s' does not exist")
    end
end

function Base.setproperty!(::CleanTable, ::Symbol, x)
    return error("Value being assigned must be an AbstractVector")
end

max_stars_repo_name: TheRoniOne/Cleaner

max_issues_repo_name: TheRoniOne/Cleaner

max_forks_repo_name: TheRoniOne/Cleaner


text:
struct PosteriorCredibleInterval{TP<:Prior,TD<:AbstractDesign,TI<:Integer,TR<:Real} <: IntervalEstimator
    prior::TP
    design::TD
    XX::Array{TI,2}
    bounds::Array{TR,2}
    alpha::TR
end

string(ci::PosteriorCredibleInterval) = @sprintf "PosteriorCredibleInterval<%s>" string(ci.prior)

function PosteriorCredibleInterval(prior::Prior, design::AbstractDesign, α::Real; epsilon = 1e-4)

    ϵ = epsilon
    XX = sample_space(design)
    nn = size(XX, 1)

    bounds = zeros(Float64, nn, 2)
    for i in 1:nn
        x1, x2    = XX[i,:]
        posterior = update(prior, x1 + x2, n(design, x1))
        phat      = (x1 + x2) / n(design, x1)
        quantile  = prob -> Roots.find_zero(p -> cdf(p, posterior) - prob, (ϵ, 1 - ϵ))
        l = try
                quantile(α)
            catch e
                isa(e, Roots.ConvergenceFailed) ? 0.0 : e
        end
        u = try
                quantile(1 - α)
            catch e
                isa(e, Roots.ConvergenceFailed) ? 1.0 : e
        end
        u <= l ? error("credible interval collapsed") : nothing
        bounds[i,:] .= [l, u]
    end
    return PosteriorCredibleInterval{typeof(prior),typeof(design),eltype(XX),Float64}(prior, design, XX, bounds, α)
end


text:
using AbstractGPs, IntervalArithmetic

mutable struct ZeroOrderGPSurrogate <: Function
    f::Function
    gps::Vector{<:GP}
    X::Vector{Vector{Float64}}
    y::Union{Vector{Vector{Float64}}, Vector{Float64}}
    noise::Float64
    std_multiple::Float64
    mode::Symbol
    N::Int
    every::Int
    skip::Int
    last::Union{Symbol, Int}
    fit_prior::Bool
    fit_noise::Bool
    counter::Int
end
function ZeroOrderGPSurrogate(
    f, x; kernel = SqExponentialKernel(), X = [x], y = [f(x)], noise = 1e-8,
    std_multiple = 3.0, mode = :interval, every = 5, skip = 10, last = :all,
    fit_prior = true, fit_noise = false,
)
    @assert noise > 0
    gps = [GP(AbstractGPs.ConstMean(0.0), kernel) for _ in 1:length(y[1])]
    N = length(y[1])
    return ZeroOrderGPSurrogate(
        f, gps, X, y, noise, std_multiple, mode, N, every, skip, last,
        fit_prior, fit_noise, 0,
    )
end

function get_x_lb_ub(k::Kernel)
    x, un = flatten(k)
    lb = fill(-Inf, length(x))
    ub = fill(Inf, length(x))
    return x, lb, ub, un
end

function fit_mle!(s, gps, X, y, noise, last, fit_noise)
    xs_uns_1 = flatten.(getproperty.(gps, :mean))
    xs_lbs_ubs_uns_2 = get_x_lb_ub.(getproperty.(gps, :kernel))

    x1 = mapreduce(vcat, xs_uns_1) do x_un
        x_un[1]
    end
    x2 = mapreduce(vcat, xs_lbs_ubs_uns_2) do x_lb_ub_un
        x_lb_ub_un[1]
    end
    lb2 = mapreduce(vcat, xs_lbs_ubs_uns_2) do x_lb_ub_un
        x_lb_ub_un[2]
    end
    ub2 = mapreduce(vcat, xs_lbs_ubs_uns_2) do x_lb_ub_un
        x_lb_ub_un[3]
    end
    if fit_noise
        x = [x1; x2; noise]
        lb = [fill(-Inf, length(x1)); lb2; 1e-8]
        ub = [fill(Inf, length(x1)); ub2; Inf]
    else
        x = [x1; x2]
        lb = [fill(-Inf, length(x1)); lb2]
        ub = [fill(Inf, length(x1)); ub2]
    end
    if last == :all
        _X = X
        _y = y
    else
        _X = X[max(end-last+1, 1):end]
        _y = y[max(end-last+1, 1):end]
    end
    obj(θ) = begin
        if _y isa Vector{<:Vector}
            offset = 0
            return -sum(map(1:length(_y[1])) do i
                l1 = length(xs_uns_1[i][1])
                un1 = xs_uns_1[i][2]
                l2 = length(xs_lbs_ubs_uns_2[i][1])
                un2 = xs_lbs_ubs_uns_2[i][4]
                _gp = GP(un1(θ[offset+1:offset+l1]), un2(θ[offset+l1+1:offset+l1+l2]))
                offset += l1 + l2
                if fit_noise
                    return logpdf(_gp(_X, θ[end]), getindex.(_y, i))
                else
                    return logpdf(_gp(_X, noise), getindex.(_y, i))
                end
            end)
        else
            l1 = length(xs_uns_1[1][1])
            un1 = xs_uns_1[1][2]
            l2 = length(xs_lbs_ubs_uns_2[1][1])
            un2 = xs_lbs_ubs_uns_2[1][4]
            _gp = GP(un1(θ[1:l1]), un2(θ[l1+1:l1+l2]))
            if fit_noise
                return -logpdf(_gp(_X, θ[end]), _y)
            else
                return -logpdf(_gp(_X, noise), _y)
            end
        end
    end
    m = Model(obj)
    addvar!(m, lb, ub)
    options = IpoptOptions(max_iter = 20, print_level = 0)
    res = optimize(m, IpoptAlg(), x, options = options)
    x = res.minimizer
    offset = 0
    gps, noise = map(1:length(_y[1])) do i
        l1 = length(xs_uns_1[i][1])
        un1 = xs_uns_1[i][2]
        l2 = length(xs_lbs_ubs_uns_2[i][1])
        un2 = xs_lbs_ubs_uns_2[i][4]
        out = GP(un1(x[offset+1:offset+l1]), un2(x[offset+l1+1:offset+l1+l2]))
        offset += l1 + l2
        return out
    end, fit_noise ? x[end] : noise
    s.gps = gps
    s.noise = noise
    return s
end
function ChainRulesCore.rrule(::typeof(fit_mle!), s, gp, X, y, noise, last, fit_noise)
    return fit_mle!(s, gp, X, y, noise, last, fit_noise), _ -> begin
        return ntuple(_ -> NoTangent(), Val(8))
    end
end

function (s::ZeroOrderGPSurrogate)(x)
    if s.mode == :exact
        y = s.f(x)
        s.X = vcat(s.X, [x])
        s.y = vcat(s.y, [y])
        s.counter += 1
        if s.fit_prior && s.counter > s.skip * s.every && (s.counter % s.every) == 0
            fit_mle!(s, s.gps, s.X, s.y, s.noise, s.last, s.fit_noise)
        end
  	    return Interval.(y, y)
    else
        if eltype(s.y) <: Real
            _m, _v = mean_and_var(posterior(
                s.gps[1](s.X, s.noise), s.y,
            ), [x])
            m, v = _m[1], _v[1]
        else
            ms_vs = map(1:s.N) do i
                _gp = s.gps[i](s.X, s.noise)
                mean_and_var(posterior(_gp, getindex.(s.y, i)), [x])
            end
            m = reduce(vcat, getindex.(ms_vs, 1))
            v = reduce(vcat, getindex.(ms_vs, 2))
        end
        r = s.std_multiple .* sqrt.(v)
        return Interval.(m .- r, m .+ r)
    end
end

function surrogate(f, x; kwargs...)
    s = ZeroOrderGPSurrogate(f, x; mode = :exact, kwargs...)
    s(x)
    s.mode = :interval
    return s
end

_lower_f(s) = x -> getproperty.(s(x), :lo)
_equality_f(s) = x -> begin
    t = s(x)
    return [getproperty.(t, :lo); .- getproperty(t, :hi)]
end

function surrogate_model(vecmodel::VecModel; kwargs...)
    surrogates = Function[]
    x0 = getinit(vecmodel)
    if :expensive in vecmodel.objective.flags
        s = surrogate(vecmodel.objective.f, copy(x0); kwargs...)
        push!(surrogates, s)
        obj = Objective(
            _lower_f(s),
            vecmodel.objective.multiple,
            vecmodel.objective.flags,
        )
    else
        obj = vecmodel.objective
    end
    expensive_eq_constraints = filter(vecmodel.eq_constraints.fs) do c
        :expensive in c.flags
    end
    cheap_eq_constraints = filter(vecmodel.eq_constraints.fs) do c
        @assert c isa EqConstraint
        !(:expensive in c.flags)
    end
    eq_constraints = VectorOfFunctions(cheap_eq_constraints)
    ineq_constraints1 = mapreduce(vcat, expensive_eq_constraints; init = Union{}[]) do c
        @assert c isa EqConstraint
        s = surrogate(c, copy(x0); kwargs...)
        push!(surrogates, s)
        return IneqConstraint(
            _equality_f(s), [zero.(c.rhs); zero.(c.rhs)], c.dim * 2, c.flags,
        )
    end
    ineq_constraints2 = map(vecmodel.ineq_constraints.fs) do c
        @assert c isa IneqConstraint
        if :expensive in c.flags
            s = surrogate(c, copy(x0); kwargs...)
            push!(surrogates, s)
            return IneqConstraint(
                _lower_f(s), zero.(c.rhs), c.dim, c.flags,
            )
        else
            return c
        end
    end
    ineq_constraints = VectorOfFunctions(
        vcat(ineq_constraints1, ineq_constraints2),
    )
    return VecModel(
        obj, eq_constraints, ineq_constraints, vecmodel.box_min,
        vecmodel.box_max, vecmodel.init, vecmodel.integer,
    ), surrogates
end

function update_surrogates!(model, surrogates, x)
    foreach(surrogates) do s
        s.mode = :exact
    end
    o = model.objective(x)
    i = model.ineq_constraints(x)
    e = model.eq_constraints(x)
    foreach(surrogates) do s
        s.mode = :interval
    end
    return o, i, e
end

@params struct BayesOptOptions
    sub_options
    maxiter::Int
    initialize::Bool
    ninit::Int
    ctol::Float64
    ftol::Float64
    postoptimize::Bool
    nt::NamedTuple
end
function BayesOptOptions(;
    sub_options = IpoptOptions(),
    initialize = true,
    ninit = 10,
    maxiter = 100,
    std_multiple = 3.0,
    kernel = SqExponentialKernel(),
    noise = 1e-8,
    ctol = 1e-4,
    ftol = 1e-4,
    postoptimize = false,
    fit_prior = true,
    fit_noise = false,
    every = 2,
    skip = 2,
    last = :all,
)
    return BayesOptOptions(
        sub_options,
        maxiter,
        initialize,
        ninit,
        ctol,
        ftol,
        postoptimize,
        (
            std_multiple = std_multiple,
            kernel = kernel,
            noise = noise,
            every = every,
            skip = skip,
            last = last,
            fit_prior = fit_prior,
            fit_noise = fit_noise,
        ),
    )
end

@params mutable struct BayesOptWorkspace <: Workspace
    model::VecModel
    sub_workspace
    x0::AbstractVector
    options::BayesOptOptions
    surrogates::AbstractVector
end

@params struct BayesOptResult <: AbstractResult
    minimizer
    minimum
    niters::Int
    sub_result::AbstractResult
    surrogates
end

# ScaledSobolIterator adapted from BayesianOptimization.jl
struct ScaledSobolIterator{T,D}
    lowerbounds::Vector{T}
    upperbounds::Vector{T}
    N::Int
    seq::SobolSeq{D}
end

"""
    ScaledSobolIterator(lowerbounds, upperbounds, N;
                        seq = SobolSeq(length(lowerbounds)))
Returns an iterator over `N` elements of a Sobol sequence between `lowerbounds`
and `upperbounds`. The first `N` elements of the Sobol sequence are skipped for
better uniformity (see https://github.com/stevengj/Sobol.jl)
"""
function ScaledSobolIterator(lowerbounds, upperbounds, N;
                             seq = SobolSeq(length(lowerbounds)))
    N > 0 && skip(seq, N)
    ScaledSobolIterator(lowerbounds, upperbounds, N, seq)
end
Base.length(it::ScaledSobolIterator) = it.N
@inline function Base.iterate(it::ScaledSobolIterator, s = 1)
    s == it.N + 1 && return nothing
    Sobol.next!(it.seq, it.lowerbounds, it.upperbounds), s + 1
end

function optimize!(workspace::BayesOptWorkspace)
    @unpack sub_workspace, surrogates = workspace
    @unpack options, x0 = workspace
    @unpack maxiter, ctol, ftol, postoptimize = options
    @unpack ninit, initialize = options
    smodel = sub_workspace.model
    minval = Inf
    minimizer = copy(x0)
    lb, ub = getmin(smodel), getmax(smodel)
    if initialize && !(all(isfinite, lb) && all(isfinite, ub))
        @warn "Skipping initialization because the variable bounds are not all finite."
    elseif initialize
        initializer = ScaledSobolIterator(lb, ub, ninit)
        x2 = x0
        for x1 in initializer
            sub_workspace.x0 .= x1
            r = optimize!(sub_workspace)
            if r.minimizer != x2
                ob, ineq, eq = update_surrogates!(smodel, surrogates, r.minimizer)
                feasible = all(ineq .<= ctol) && all(abs.(eq) .<= ctol)
                if feasible && ob <= minval
                    minval = ob
                    minimizer = copy(r.minimizer)
                end
            end
            x2 = r.minimizer
        end
    end
    sub_workspace.x0 .= minimizer
    r = optimize!(sub_workspace)
    x = copy(r.minimizer)
    prevm = minval
    iter = 1
    while iter < maxiter
        iter += 1
        ob, ineq, eq = update_surrogates!(smodel, surrogates, x)
        feasible = all(ineq .<= ctol) && all(abs.(eq) .<= ctol)
        if feasible && ob <= minval
            prevm = minval
            minval = ob
            minimizer = copy(x)
        end
        converged = feasible && abs((ob - prevm) / (prevm + ftol)) <= ftol
        converged && break
        r = optimize!(sub_workspace)
        m = ob
        x = copy(r.minimizer)
    end
    ob, ineq, eq = update_surrogates!(smodel, surrogates, x)
    feasible = all(ineq .<= ctol) && all(abs.(eq) .<= ctol)
    if feasible && ob <= minval
        minval = ob
        minimizer = r.minimizer
    end
    if postoptimize
        foreach(surrogates) do s
            s.mode = :exact
        end
    end
    sub_workspace.x0 .= minimizer
    r = optimize!(sub_workspace)
    if postoptimize
        if r.minimum <= minval
            minval = r.minimum
            minimizer = copy(r.minimizer)
        end
        foreach(surrogates) do s
            s.mode = :interval
        end
    end
    return BayesOptResult(
        minimizer, minval, iter, r, surrogates,
    )
end

struct BayesOptAlg{A} <: AbstractOptimizer
    sub_alg::A
end

function Workspace(
    model::VecModel, optimizer::BayesOptAlg, x0::AbstractVector;
    options::BayesOptOptions, surrogates = nothing, kwargs...,
)
    if surrogates !== nothing
        smodel = model
    else
        smodel, surrogates = surrogate_model(model; options.nt...)
    end
    sub_workspace = Workspace(
        smodel, optimizer.sub_alg, x0;
        options = options.sub_options, kwargs...,
    )
    return BayesOptWorkspace(
        model,
        sub_workspace,
        x0,
        options,
        surrogates,
    )
end


text:
# Test using the Random data
@testset "Random" begin

    @testset "Negative Moran" begin
        geom = reggeomlattice(10, 10)
        W = polyneigh(geom)
        X = rand(StableRNG(1234567), 100)
        Wmat = Matrix(W)
        rho = -0.9

        Y = X + rho * Wmat * X
        mnegative = moran(Y, W, rng = StableRNG(1234567))

        @test score(mnegative) ≈ -0.197453 atol = 1e-5
        @test mean(mnegative) ≈ -0.010786 atol = 1e-5
        @test std(mnegative) ≈ 0.053380 atol = 1e-5
        @test zscore(mnegative) ≈ -3.496951 atol = 1e-5
        @test pvalue(mnegative) ≈ 0.0001 atol = 1e-5
    end

end


text:
abstract type AbstractEvent end
abstract type Environment end

@enum EVENT_STATE idle=0 scheduled=1 processed=2

struct EventProcessed <: Exception
  ev :: AbstractEvent
end

struct EventNotIdle <: Exception
  ev :: AbstractEvent
end

mutable struct BaseEvent
  env :: Environment
  id :: UInt
  callbacks :: Vector{Function}
  state :: EVENT_STATE
  value :: Any
  function BaseEvent(env::Environment)
    new(env, env.eid+=one(UInt), Vector{Function}(), idle, nothing)
  end
end

function show(io::IO, ev::AbstractEvent)
  print(io, "$(typeof(ev)) $(ev.bev.id)")
end

function show(io::IO, env::Environment)
  if env.active_proc == nothing
    print(io, "$(typeof(env)) time: $(now(env)) active_process: nothing")
  else
    print(io, "$(typeof(env)) time: $(now(env)) active_process: $(active_proc)")
  end
end

function environment(ev::AbstractEvent) :: Environment
  ev.bev.env
end

function value(ev::AbstractEvent) :: Any
  ev.bev.value
end

function state(ev::AbstractEvent) :: EVENT_STATE
  ev.bev.state
end

function append_callback(func::Function, ev::AbstractEvent, args::Any...) :: Function
  ev.bev.state == processed && throw(EventProcessed(ev))
  cb = ()->func(ev, args...)
  push!(ev.bev.callbacks, cb)
  cb
end

macro callback(expr::Expr)
  expr.head != :call && error("Expression is not a function call!")
  esc(:(SimHPC.append_callback($(expr.args...))))
end

function remove_callback(cb::Function, ev::AbstractEvent)
  i = indexin(ev.bev.callbacks, [cb])[1]
  i != 0 && deleteat!(ev.bev.callbacks, i)
end

function schedule(ev::AbstractEvent, delay::Number=zero(Float64); priority::Int=0, value::Any=nothing)
  state(ev) == processed && throw(EventProcessed(ev))
  env = environment(ev)
  bev = ev.bev
  bev.value = value
  env.heap[bev] = EventKey(now(env) + delay, priority, env.sid+=one(UInt))
  bev.state = scheduled
  ev
end

struct StopSimulation <: Exception
  value :: Any
  function StopSimulation(value::Any=nothing)
    new(value)
  end
end

function stop_simulation(ev::AbstractEvent)
  throw(StopSimulation(value(ev)))
end

function run(env::Environment, until::AbstractEvent)
  @callback stop_simulation(until)
  try
    while true
      step(env)
    end
  catch exc
    if isa(exc, StopSimulation)
      return exc.value
    else
      rethrow(exc)
    end
  end
end

max_stars_repo_name: zhangliye/SimHPC

max_issues_repo_name: zhangliye/SimHPC

max_forks_repo_name: zhangliye/SimHPC


text:
module WrappedR

using RCall

include("sprintr.jl")

export fitsprintr, predictsprintr

end


text:
#
# Protein - Water (no plots, to compute external timming)
#

using MDDF
using PDBTools

# Here we use the PDBTools package to read the pdb file (from http://github.com/m3g/PDBTools)
atoms = PDBTools.readPDB("./structure.pdb")

# The solute is a single protein molecule (infinte dilution case). In this case,
# use the option nmols=1
solute_indexes = PDBTools.select(atoms,"protein")
solute = MDDF.Selection( solute_indexes, nmols=1 )

# The solvent is Water, which has 3 atoms. Use the natomspermol to indicate how many
# atoms each molecule has, such that there is no ambiguity on how to split the coordinates
# of the selection into individual molecules.
solvent_indexes = PDBTools.select(atoms,"water")
solvent = MDDF.Selection( solvent_indexes, natomspermol=3 )

# Initialize trajectroy data structure and open input stream
trajectory = MDDF.Trajectory("./trajectory.dcd",solute,solvent)

# Input options for the calcualtion
options = MDDF.Options(binstep=0.2)

# Run MDDF calculation, and get the resutls in the R structure
@time R = MDDF.mddf_linkedcells_parallel(trajectory,options)


text:
using Sim, Printf

struct Guy
    name
end

abstract type Encounter end
struct Meet <: Encounter
    someone
end
struct Greet <: Encounter
    num
    from
end
struct Response <: Encounter
    num
    from
end

comm = ("Nice to meet you!", "How are you?", "Have a nice day!", "bye bye")
say(name, n) =  @printf("%5.2f s, %s: %s\n", tau(), name, comm[n])

function step!(me::Guy, σ::Meet)
    event!(𝐶, SimFunction(step!, σ.someone, Greet(1, me)), after, 2*rand())
    say(me.name, 1)
end

function step!(me::Guy, σ::Greet)
    if σ.num < 3
        event!(𝐶, SimFunction(step!, σ.from, Response(σ.num, me)), after, 2*rand())
        say(me.name, σ.num)
    else
        say(me.name, 4)
    end
end

function step!(me::Guy, σ::Response)
    event!(𝐶, SimFunction(step!, σ.from, Greet(σ.num+1, me)), after, 2*rand())
    say(me.name, σ.num+1)
end

foo = Guy("Foo")
bar = Guy("Bar")

event!(𝐶, SimFunction(step!, foo, Meet(bar)), at, 10*rand())
run!(𝐶, 20)

reset!(𝐶)