text:
#=
array_operation:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-10
=#

module ArrayOp
    # packages
    using LinearAlgebra

    # array
    A = randn(Float32, 3, 2)
    B = collect(reshape(1:9, 3, 3))
    C = fill(4, 3, 3)
    D = fill(2, 3, 3)

    # methods
    function info()
        println("Array A = $A")
        println("Element type of A = ", eltype(A))
        println("Element num of A = ", length(A))
        println("Dimension num of A = ", ndims(A))
        println("Size of A = ", size(A))
        println("Size of A, dim:1 = ", size(A, 1))
        println("Size of A, dim:2 = ", size(A, 2))
        println("Strides of A = ", strides(A))
        println("Stride of A, 2 = ", stride(A, 2))
        println("")
    end

    function indexing()
        println("Array B = $B")
        println("B[3, 3] = ", B[3, 3])
        println("B[2, :] = ", B[2,:])
        println("B[:, 1] = ", B[:,1])
        B[1, 2] = -10
        println("B[1, 2] = -10: ", B)
        B[2:3, 2:3] = [-5 -5; -3 -3]
        println("B[2:3, 2:3] = -10: ", B)
        println("")
    end

    function four_arithmetic()
        println("Array C = $C")
        println("Array D = $D")
        println("C + D = ", C + D)
        println("C - D = ", C - D)
        println("C * D = ", C * D)
        println("")
    end

    function sigmoid(x::Float64)
        return 1.0 / (1.0 + exp(-x))
    end

    function broad_cast()
        println("C .+ D = ", C .+ D)
        println("C .+ 10 = ", C .+ 10)
        println("C .- D = ", C .- D)
        println("C .- 10 = ", C .- 10)
        println("C .* D = ", C .* D)
        println("C .* 10 = ", C .* 10)
        println("C ./ D = ", C ./ D)
        println("C ./ 10 = ", C ./ 10)
        println("Element type of (C./D) = ", eltype(C ./ D))
        println("sigmoid.() = ", sigmoid.(rand(3, 2)))
        println("")
    end

    function map_red_fil()
        E = rand(4, 3)
        println("Array E = ", E)
        println("map(x -> x + 1.0, E) = ", map(x -> x + 1.0, E))
        println("All of mult, reduce(*, E) = ", reduce(*, E))
        println("filter(x -> x < 0.5, E) = ", filter(x -> x < 0.5, E))
        println("")
    end

    function matrix()
        F = rand(5, 5)
        println("Array F = ", F)
        println("Inversed Matrix inv(F) = ", inv(F))
        println("Transposed Matrix F' = ", F')
        println("Eigen Values eigvals(F) = ", eigvals(F))
        println("Eigen Vectors eigvecs(F) = ", eigvecs(F))
        println("")
    end

    function copy_shallow()
        println("Shallow copy")
        G = rand(5, 5)
        println("Array G = ", G)
        H = G
        println("H = G: ", H)
        H .= 0.0
        println("H .= 0.0: ", H)
        println("Array G = ", G)
    end

    function copy_deep()
        println("Deep copy")
        G = rand(5, 5)
        println("Array G = ", G)
        H = copy(G)
        println("H = copy(G): ", H)
        H .= 0.0
        println("H .= 0.0: ", H)
        println("Array G = ", G)
    end

    function main()
        info()
        indexing()
        four_arithmetic()
        broad_cast()
        map_red_fil()
        matrix()
        copy_shallow()
        copy_deep()
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .ArrayOp
    ArrayOp.main()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
type Literal <: AbstractPlanOutOp
  args::Dict{AbstractString, Any}
end

literal(args::Dict) = Literal(args)
execute(op::Literal, i::AbstractPlanOutMapper) = getArgMixed(op, "value")

type Set <: AbstractPlanOutOp
  args::Dict{AbstractString, Any}
end
set(args::Dict) = Set(args)

function execute(op::Set, i::AbstractPlanOutMapper)
  var = getArgString(op, "var")
  val = getArgMixed(op, "value")
  if var == "experiment_salt"
    setsalt!(i, value)
  end
  setindex!(i, evaluate(i, val), var)
end

type Get <: AbstractPlanOutOp
  args::Dict{AbstractString, Any}
end
get(args::Dict) = Get(args)
execute(op::Get, i::AbstractPlanOutMapper) = getindex(i, getArgString("var"))

type Seq <: AbstractPlanOutOp
  args::Dict{AbstractString, Any}
end
seq(args::Dict) = Seq(args)
execute(op::Seq, i::AbstractPlanOutMapper) = [evaluate(i, el) for el in getArgList(op, "seq")]

type Array <: AbstractPlanOutOp
  args::Dict{AbstractString, Any}
end
array(args::Dict) = Array(args)
execute(op::Array, i::AbstractPlanOutMapper) = [evaluate(i, el) for el in getArgList(op, "values")]

type Coalesce <: AbstractPlanOutOp
  args::Dict{AbstractString, Any}
end
coalesce(args::Dict) = Coalesce(args)

nullable_cond(val::Nullable) = isnull(val) ? false : get(val)
nullable_cond(val::Any) = val
function execute(op::Coalesce, i::AbstractPlanOutMapper)
  for j in getArgList(op, "values")
    eval_x = nullable_cond(evaluate(i, j))
    if eval_x != false
      return eval_x
    end
  end
  return Nullable()
end

type Cond <: AbstractPlanOutOp
  args::Dict{AbstractString, Any}
end
cond(args::Dict) = Cond(args)

function execute(op::Cond, i::AbstractPlanOutMapper)
  for i in getArgList(op, "cond")
    if_clause = getindex(i, "if")
    if evaluate(i, if_clause)
      then_clause = getindex(i, "then")
      return evaluate(i, then_clause)
    end
  end
end

type And <: AbstractPlanOutOp
  args::Dict{AbstractString, Any}
end
and(args::Dict) = And(args)

function execute(op::And, i::AbstractPlanOutMapper)
  for clause in getArgList(op, "values")
    if !Bool(evaluate(i, clause))
      return false
    end
  end
  return true
end

type Or <: AbstractPlanOutOp
  args::Dict{AbstractString, Any}
end
or(args::Dict) = Or(args)

function execute(op::Or, i::AbstractPlanOutMapper)
  for clause in getArgList(op, "values")
    if Bool(evaluate(i, clause))
      return true
    end
  end
  return false
end

type Map <: AbstractPlanOutOpSimple
  args::Dict{AbstractString, Any}
end
map(args::Dict) = Map(args)

function simpleExecute(op::Map, i::AbstractPlanOutMapper)
  c = copy(op.args)
  delete!(c, "op")
  delete!(c, "salt")
end

type Product <: AbstractPlanOutOpCommutative
  args::Dict{AbstractString, Any}
end
product(args::Dict) = Product(args)

execute(op::Product, i::AbstractPlanOutMapper, values::AbstractArray) = prod(values)

type Sum <: AbstractPlanOutOpCommutative
  args::Dict{AbstractString, Any}
end
sum(args::Dict) = Sum(args)

execute(op::Sum, i::AbstractPlanOutMapper, values::AbstractArray) = sum(values)

type Min <: AbstractPlanOutOpCommutative
  args::Dict{AbstractString, Any}
end
min(args::Dict) = Min(args)

execute(op::Min, i::AbstractPlanOutMapper, values::AbstractArray) = minimum(values)

type Max <: AbstractPlanOutOpCommutative
  args::Dict{AbstractString, Any}
end
max(args::Dict) = Max(args)

execute(op::Max, i::AbstractPlanOutMapper, values::AbstractArray) = maximum(values)

type Equals <: AbstractPlanOutOpBinary
  args::Dict{AbstractString, Any}
end
equals(args::Dict) = Equals(args)

execute(op::Equals, i::AbstractPlanOutMapper, left, right) = left == right

type GreaterThan <: AbstractPlanOutOpBinary
  args::Dict{AbstractString, Any}
end
greaterthan(args::Dict) = GreaterThan(args)

execute(op::GreaterThan, i::AbstractPlanOutMapper, left, right) = left > right

type LessThan <: AbstractPlanOutOpBinary
  args::Dict{AbstractString, Any}
end
lessthan(args::Dict) = LessThan(args)

execute(op::LessThan, i::AbstractPlanOutMapper, left, right) = left < right

type LessThanOrEqualTo <: AbstractPlanOutOpBinary
  args::Dict{AbstractString, Any}
end
lessthanorequalto(args::Dict) = LessThanOrEqualTo(args)

execute(op::LessThanOrEqualTo, i::AbstractPlanOutMapper, left, right) = left <= right

type GreaterThanOrEqualTo <: AbstractPlanOutOpBinary
  args::Dict{AbstractString, Any}
end
greaterthanorequalto(args::Dict) = GreaterThanOrEqualTo(args)

execute(op::GreaterThanOrEqualTo, i::AbstractPlanOutMapper, left, right) = left >= right

type Mod <: AbstractPlanOutOpBinary
  args::Dict{AbstractString, Any}
end
mod(args::Dict) = Mod(args)

execute(op::Mod, i::AbstractPlanOutMapper, left, right) = left % right

type Divide <: AbstractPlanOutOpBinary
  args::Dict{AbstractString, Any}
end
div(args::Dict) = Divide(args)

execute(op::Divide, i::AbstractPlanOutMapper, left, right) = float(left) / float(right)

type Round <: AbstractPlanOutOpUnary
  args::Dict{AbstractString, Any}
end
round(args::Dict) = Round(args)

execute(op::Round, i::AbstractPlanOutMapper, value::Number) = round(value)

type Not <: AbstractPlanOutOpUnary
  args::Dict{AbstractString, Any}
end
not(args::Dict) = Not(args)

execute(op::Not, i::AbstractPlanOutMapper, value) = !Bool(value)

type Negative <: AbstractPlanOutOpUnary
  args::Dict{AbstractString, Any}
end
negative(args::Dict) = Negative(args)

execute(op::Negative, i::AbstractPlanOutMapper, value::Number) = 0 - value

type Length <: AbstractPlanOutOpUnary
  args::Dict{AbstractString, Any}
end
len(args::Dict) = Length(args)

execute(op::Length, i::AbstractPlanOutMapper, value) = Base.length(value)

type Index <: AbstractPlanOutOpSimple
  args::Dict{AbstractString, Any}
end
index(args::Dict) = Index(args)

modifyIndex(a::Any) = a
modifyIndex(a::Number) = a + 1 # because other languages start with index = 0
function simpleExecute(op::Index, i::AbstractPlanOutMapper)
  base = getArgIndexish(op, "base")
  index = modifyIndex(getArgMixed(op, "index"))
  Base.get(base, index, Nullable())
end

type Return <: AbstractPlanOutOp
  args::Dict{AbstractString, Any}
end
ret(args::Dict) = Return(args)
function execute(op::Return, i::AbstractPlanOutMapper)
  value = evaluate(i, getArgMixed(op, "value"))
  in_experiment = value > 0
  throw(StopPlanOutException(in_experiment))
end

ext: jl


text:
using JSON
import Base.show

type Note
    id::Int
    tags::Array{AbstractString, 1}
    body::AbstractString
end

function grab_json(fname="")
    res = open(fname, "r") do f
        raw_json = readstring(f)
        res = JSON.parse(raw_json)
    end
    return res
end

function read_notes(raw_notes_list)
    notes_list = []
    for raw_note in raw_notes_list
        id = raw_note["id"]
        tags = raw_note["tags"]
        body = raw_note["body"]
        note = Note(id, tags, body)
        push!(notes_list, note)
    end

    return notes_list
end

function read_notes(fname::AbstractString)
    try
        raw_notes_list = grab_json(fname)
        return read_notes(raw_notes_list)
    catch e
        if e isa LoadError
            println("File '", fname, "' does not exist. Starting with no notes")
        else
            println("Error loading notes file. Starting with no notes")
        end
        return []
    end
end

function write_dict(notes_dict, fname)
    notes_list = collect(values(notes_dict))
    open(fname, "w") do f
        write(f, json(notes_list))
    end
end

function add_note(note::Note, note_dict::Dict, tag_dict::Dict; writeback=true, fname="notes.json")
    note_dict[note.id] = note
    for tag in note.tags
        if haskey(tag_dict, tag)
            push!(tag_dict[tag], note.id)
        else
            tag_dict[tag] = [note.id]
        end
    end

    if writeback
        write_dict(note_dict, fname)
    end
end

function create_dicts(notes_list::Array)
    note_dict=Dict()
    tag_dict=Dict()
    for note in notes_list
        add_note(note, note_dict, tag_dict)
    end
    return note_dict, tag_dict
end

function create_dicts(fname::AbstractString)
    return create_dicts(read_notes(fname))
end

function search_tags(tag, tag_dict)
    tag_list = get(tag_dict, tag, [])
    return tag_list
end

function search_tags(tag, tag_dict, notes_dict)
    tag_list = search_tags(tag, tag_dict)
    return [notes_dict[idx] for idx in tag_list]
end

function print_note_body(body::AbstractString)
    md = Markdown.parse(body)
    display(md)
end

function print_tags(tags)
    tag_color_str = Base.text_colors[:blue]
    tag_str = string(tag_color_str, "  ", map(t->string("#", t, " "), tags)...)
    println(tag_str, Base.text_colors[:normal])
end

function print_note(note::Note)
    print_note_body(note.body)
    println()
    print_tags(note.tags)
end

dashed_seperator() = println(string(["=" for _ in 1:40]...))

function get_notes(tags::Array, tag_dict::Dict, notes_dict::Dict; seperator_fn=dashed_seperator)
    importance = Dict()

    for tag in tags
        id_list = search_tags(tag, tag_dict)
        for id in id_list
            if haskey(importance, id)
                importance[id] += 1
            else
                importance[id] = 1
            end
        end
    end

    sorted_ids = sort(collect(keys(importance)), by=x->importance[x], rev=true)
    notes = [notes_dict[id] for id in sorted_ids]

    if isempty(notes)
        println("No notes found containing tags in [", tags..., "]")
    else
        seperator_fn()
        for note in notes
            print_note(note)
            seperator_fn()
        end
    end
    nothing
end

function get_notes(tag::AbstractString, tag_dict, notes_dict)
    get_notes([tag], tag_dict, notes_dict)
end

function verify_notes_file(fname)
    if !(basename(fname) in readdir(dirname(fname)))
        write_dict(Dict(0=>Note(0, [], "dummy")), fname)
    end
end


text:
#!/usr/bin/env julia
#
using DataFrames
showln(x) = (show(x); println())

println(" Running with Julia ", VERSION)

# This is now required to use constant 'e'
using Base.MathConstants
# A DataFrame is an in-memory database
df = DataFrame(A = [1, 2], B = [e, pi], C = ["xx", "xy"])
showln(df)
#> 2x3 DataFrame
#> |-------|---|---------|------|
#> | Row # | A | B       | C    |
#> | 1     | 1 | 2.71828 | "xx" |
#> | 2     | 2 | 3.14159 | "xy" |

# The columns of a DataFrame can be indexed using numbers or names
showln(df[!, 1])
#> [1,2]
showln(df[!, :A])
#> [1,2]

showln(df[!, 2])
#> [2.718281828459045,3.141592653589793]
showln(df[!, :B])
#> [2.718281828459045,3.141592653589793]

showln(df[!, 3])
#> String["xx","xy"]
showln(df[!, :C])
#> String["xx","xy"]

# The rows of a DataFrame can be indexed only by using numbers
showln(df[1, :])
#> 1x3 DataFrame
#> |-------|---|---------|------|
#> | Row # | A | B       | C    |
#> | 1     | 1 | 2.71828 | "xx" |
showln(df[1:2, :])
#> 2x3 DataFrame
#> |-------|---|---------|------|
#> | Row # | A | B       | C    |
#> | 1     | 1 | 2.71828 | "xx" |
#> | 2     | 2 | 3.14159 | "xy" |

# importing data into DataFrames
# ------------------------------

# DataFrames can be loaded from CSV files using readtable()
using CSV
iris = CSV.read("iris.csv")

# the iris dataset (and plenty of others) is also available from
using RData, RDatasets
iris = dataset("datasets","iris")

# you can directly import your own R .rda dataframe with
# mydf = DataFrame(read_rda("path/to/your/df.rda")["name_of_df"]), e.g.
# Warning: `Pkg.dir(pkgname, paths...)` is deprecated; instead, do `import RDatasets; joinpath(dirname(pathof(RDatasets)), "..", paths...)`.
using Pkg
diamonds = DataFrame(load(joinpath(Pkg.dir("RDatasets"),"data","ggplot2","diamonds.rda"))["diamonds"])

# showing DataFrames
# ------------------

# Check the names and element types of the columns of our new DataFrame
showln(names(iris))
#> [:SepalLength,:SepalWidth,:PetalLength,:PetalWidth,:Species]
showln(eltypes(iris))
#> Type[Float64,Float64,Float64,Float64,UTF8String]

# Subset the DataFrame to only include rows for one species
# Warning: `getindex(df::DataFrame, col_ind::ColumnIndex)` is deprecated, use `df[!, col_ind]` instead.
showln(iris[iris[:Species] .== "setosa", :])
#> 50x5 DataFrame
#> |-------|-------------|------------|-------------|------------|----------|
#> | Row # | SepalLength | SepalWidth | PetalLength | PetalWidth | Species  |
#> | 1     | 5.1         | 3.5        | 1.4         | 0.2        | "setosa" |
#> | 2     | 4.9         | 3.0        | 1.4         | 0.2        | "setosa" |
#> | 3     | 4.7         | 3.2        | 1.3         | 0.2        | "setosa" |
#> | 4     | 4.6         | 3.1        | 1.5         | 0.2        | "setosa" |
#> | 5     | 5.0         | 3.6        | 1.4         | 0.2        | "setosa" |
#> | 6     | 5.4         | 3.9        | 1.7         | 0.4        | "setosa" |
#> | 7     | 4.6         | 3.4        | 1.4         | 0.3        | "setosa" |
#> | 8     | 5.0         | 3.4        | 1.5         | 0.2        | "setosa" |
#> | 9     | 4.4         | 2.9        | 1.4         | 0.2        | "setosa" |
#> ⋮
#> | 41    | 5.0         | 3.5        | 1.3         | 0.3        | "setosa" |
#> | 42    | 4.5         | 2.3        | 1.3         | 0.3        | "setosa" |
#> | 43    | 4.4         | 3.2        | 1.3         | 0.2        | "setosa" |
#> | 44    | 5.0         | 3.5        | 1.6         | 0.6        | "setosa" |
#> | 45    | 5.1         | 3.8        | 1.9         | 0.4        | "setosa" |
#> | 46    | 4.8         | 3.0        | 1.4         | 0.3        | "setosa" |
#> | 47    | 5.1         | 3.8        | 1.6         | 0.2        | "setosa" |
#> | 48    | 4.6         | 3.2        | 1.4         | 0.2        | "setosa" |
#> | 49    | 5.3         | 3.7        | 1.5         | 0.2        | "setosa" |
#> | 50    | 5.0         | 3.3        | 1.4         | 0.2        | "setosa" |

# Count the number of rows for each species
showln(by(iris, :Species, df -> size(df, 1)))
#> 3x2 DataFrame
#> |-------|--------------|----|
#> | Row # | Species      | x1 |
#> | 1     | "setosa"     | 50 |
#> | 2     | "versicolor" | 50 |
#> | 3     | "virginica"  | 50 |

# Discretize entire columns at a time
iris[!, :SepalLength] = round.(Integer, iris[!, :SepalLength])
iris[!, :SepalWidth] = round.(Integer, iris[!, :SepalWidth])

# Tabulate data according to discretized columns to see "clusters"
tabulated = by(
    iris,
    [:Species, :SepalLength, :SepalWidth],
    df -> size(df, 1)
)
showln(tabulated)
#> 17x4 DataFrame
#> |-------|--------------|-------------|------------|----|
#> | Row # | Species      | SepalLength | SepalWidth | x1 |
#> | 1     | "setosa"     | 4           | 3          | 4  |
#> | 2     | "setosa"     | 5           | 2          | 1  |
#> | 3     | "setosa"     | 5           | 3          | 23 |
#> | 4     | "setosa"     | 5           | 4          | 17 |
#> | 5     | "setosa"     | 6           | 4          | 5  |
#> | 6     | "versicolor" | 5           | 2          | 3  |
#> | 7     | "versicolor" | 5           | 3          | 3  |
#> | 8     | "versicolor" | 6           | 2          | 6  |
#> | 9     | "versicolor" | 6           | 3          | 29 |
#> | 10    | "versicolor" | 7           | 3          | 9  |
#> | 11    | "virginica"  | 5           | 3          | 1  |
#> | 12    | "virginica"  | 6           | 2          | 1  |
#> | 13    | "virginica"  | 6           | 3          | 22 |
#> | 14    | "virginica"  | 7           | 3          | 19 |
#> | 15    | "virginica"  | 7           | 4          | 1  |
#> | 16    | "virginica"  | 8           | 3          | 4  |
#> | 17    | "virginica"  | 8           | 4          | 2  |

# you can setup a grouped dataframe like this
gdf = groupby(iris,[:Species, :SepalLength, :SepalWidth])

# and then iterate over it
for idf in gdf
	println(" idf size => ", size(idf,1))
end

# Adding/Removing columns
# -----------------------

# insert!(df::DataFrame,index::Int64,item::AbstractArray{T,1},name::Symbol)
# insert random numbers at col 5:
# insert!(iris, 5, rand(nrow(iris)), :randCol)
# insertcols!(df, col_ind, name => item; makeunique=makeunique)
insertcols!(iris, 5, :randCol => rand(nrow(iris)))

# remove it
# delete!(iris, :randCol)
select!(iris, Not(:randCol))

max_stars_repo_name: wbeebe/julia

max_issues_repo_name: wbeebe/julia

max_forks_repo_name: wbeebe/julia


text:
module Competition

using Test

include("../../example/competition.jl")

g_result = run_gillespie(MersenneTwister(1), (0.0, T))[1]
g_result_rec = run_gillespie_record(MersenneTwister(1), T)[1]
m_result = run_manually(MersenneTwister(1), T)
m_result_rec = run_manually_record(MersenneTwister(1), T)

@test g_result.v == g_result_rec.v == m_result.v == m_result_rec.v
@test g_result.m == g_result_rec.m == m_result.m == m_result_rec.m

@test all(map(==, getentries(g_result_rec.v), getentries(m_result_rec.v)))
@test all(map(==, values(getentries(g_result_rec.m)), values(getentries(m_result_rec.m))))

end


text:
"""
    blosc_cbuffer_sizes(cbuffer)::NamedTuple{(:nbytes, :cbytes, blocksize), Tuple{Int, Int, Int}}

Get information about a compressed buffer, namely the number of uncompressed bytes (`nbytes`) and compressed (`cbytes`).

It also returns the `blocksize` (which is used internally for doing the compression by blocks).

You only need to pass the first BLOSC_EXTENDED_HEADER_LENGTH bytes of a compressed buffer for this call to work.

This function should always succeed.

C signature `void blosc_cbuffer_sizes(const void *cbuffer, size_t *nbytes, size_t *cbytes, size_t *blocksize)`
"""
function blosc_cbuffer_sizes(cbuffer)::NamedTuple{(:nbytes, :cbytes, :blocksize), Tuple{Int, Int, Int}}

    nbytes = Ref{Csize_t}()
    cbytes = Ref{Csize_t}()
    blocksize = Ref{Csize_t}()
    @ccall lib.blosc_cbuffer_sizes(cbuffer::Ptr{Cvoid}, nbytes::Ptr{Csize_t}, cbytes::Ptr{Csize_t}, blocksize::Ptr{Csize_t})::Cvoid
    return (
        nbytes = nbytes[],
        cbytes = cbytes[],
        blocksize = blocksize[]
    )
end

"""
    blosc_cbuffer_metainfo(cbuffer)::NamedTuple{(:typesize, :flags), Tuple{Int, Int}}

Get information about a compressed buffer, namely the type size (`typesize`), as well as some internal `flags`.

You can use the `BLOSC_DOSHUFFLE`, `BLOSC_DOBITSHUFFLE`, `BLOSC_DODELTA` and `BLOSC_MEMCPYED` symbols for extracting the interesting bits (e.g. `flags & BLOSC_DOSHUFFLE` says whether the buffer is byte-shuffled or not).

This function should always succeed.

C signature `void blosc_cbuffer_metainfo(const void *cbuffer, size_t *typesize, int *flags)`

# `flags` bits:

- bit 0: whether the shuffle filter has been applied or not

- bit 1: whether the internal buffer is a pure memcpy or not

- bit 2: whether the bitshuffle filter has been applied or not

- bit 3: whether the delta coding filter has been applied or not

"""
function blosc_cbuffer_metainfo(cbuffer)::NamedTuple{(:typesize, :flags), Tuple{Int, Int}}
    typesize = Ref{Csize_t}()
    flags = Ref{Cint}()
    @ccall lib.blosc_cbuffer_metainfo(cbuffer::Ptr{Cvoid}, typesize::Ptr{Csize_t}, flags::Ptr{Cint})::Cvoid
    return (
        typesize = typesize[],
        flags = flags[]
    )
end

"""
    blosc_cbuffer_versions(cbuffer)::NamedTuple{(:version, :versionlz), Tuple{Int, Int}}

Get information about a compressed buffer, namely the internal Blosc format version (`version`) and the format for the internal Lempel-Ziv compressor used (`versionlz`).

This function should always succeed.

C signature `void blosc_cbuffer_versions(const void *cbuffer, int *version, int *versionlz)`
"""
function blosc_cbuffer_versions(cbuffer)::NamedTuple{(:version, :versionlz), Tuple{Int, Int}}
    version = Ref{Cint}()
    versionlz = Ref{Cint}()
    @ccall lib.blosc_cbuffer_versions(cbuffer::Ptr{Cvoid}, version::Ptr{Cint}, versionlz::Ptr{Cint})::Cvoid
    return (
        version = version[],
        versionlz = versionlz[]
    )
end

"""
    blosc_cbuffer_complib(cbuffer)

Get the compressor library/format used in a compressed buffer.

This function should always succeed.
"""
function blosc_cbuffer_complib(cbuffer)
    unsafe_string(
        @ccall lib.blosc_cbuffer_complib(cbuffer::Ptr{Cvoid})::Cstring
    )
end


text:
export RNNBlock

function RNNBlock(name::Symbol, net::Net, input::AbstractEnsemble)
    @assert ndims(input) == 1 "RNNLayer currently only supports 1-d inputs"
    x = InnerProductLayer(symbol(name, :x), net, input, length(input))
    state = FullyConnectedEnsemble(symbol(name, :state), net, length(input),
                                   length(input))
    # h = AddLayer(symbol(name, :tanh), net, x, state)
    sum = AddLayer(symbol(name, :sum), net, x, state)
    h = TanhLayer(symbol(name, :tanh), net, sum)
    add_connections(net, h, state,
                    (i) -> (tuple([Colon() for d in size(state)]... ));
                    recurrent=true)
    InnerProductLayer(name, net, h, length(input))
end


text:
module TransferEntropy

using Reexport

using Distributions
using Distances
using SpecialFunctions
using NearestNeighbors
@reexport using CausalityToolsBase
@reexport using PerronFrobenius
using StateSpaceReconstruction; export invariantize

using StatsBase

include("TEVars.jl")
include("estimators/transferentropy_kraskov.jl")
include("estimators/transferentropy_visitfreq.jl")
include("estimators/transferentropy_transferoperator.jl")
include("estimators/common_interface.jl")
include("te_embed.jl")

include("convenience_funcs_regular_te.jl")
include("convenience_funcs_conditional_te.jl")

export
transferentropy
TEVars

end # module


text:
#import Plots
import PyPlot

"""
Generates data for plotting f on [a, b]
"""
function plt_data(f, a, b, numpts=100)
    h = (b-a)/(numpts-1)
    y = zeros(1,numpts)
    x = collect(a:h:b)
    y = map(f, x)
    return (x, y)
end

"""
Displays graph of f on interval [a, b]
"""
function plt(f, a, b, numpts=100; label = "y", title="", color = nothing)
    h = (b-a)/(numpts-1)
    y = zeros(1,numpts)
    x = collect(a:h:b)
    y = map(f, x)
    if color!=nothing
        display(Plots.plot!(x, y, label = label, title = title, color = color, linewidth = 3))
    else
        display(Plots.plot!(x, y, label = label, title = title, linewidth = 3))
    end
    return (x, y)
end

function plt!(f, a, b, numpts=100; label = "y", title = "", color = nothing)
    h = (b-a)/(numpts-1)
    y = zeros(1,numpts)
    x = collect(a:h:b)
    y = map(f, x)
    if color!=nothing
        display(Plots.plot!(x, y, label = label, title = title, color = color, linewidth = 3))
    else
        display(Plots.plot!(x, y, label = label, title = title, linewidth = 3))
    end
    return(x, y)
end

function plt2D(f, arr1, arr2, numpts=100; label = "y", title="", color = nothing)
    a1 = arr1[1]
    b1 = arr1[2]
    a2 = arr2[1]
    b2 = arr2[2]

    h = (b-a)/(numpts-1)

    x1 = collect(a1:h:b1)
    x2 = collect(a2:h:b2)

    W = zeros(n, n)
    C = CartesianIndices(W)

    for I in C
        W[I] = f(x1[I[1]], x2[I[2]])
    end

    y = zeros(1,numpts)
    x = collect(a:h:b)
    y = map(f, x)
    if color!=nothing
        display(Plots.plot!(x, y, label = label, title = title, color = color, linewidth = 3))
    else
        display(Plots.plot!(x, y, label = label, title = title, linewidth = 3))
    end
    return (x, y)
end

function plt!(point::Float64; label = "", title = "")
    display(Plots.plot!([point], seriestype = :vline))
end

"""
plot pdf, median, mode and credible interval
"""
function plot_distribution(pdf, median, y0_true; CI=nothing, mytitle="Posterior Distribution")
    # plot pdf
    b = 1.2
    while pdf(b) > 0.05 && b < 10
        b += 0.5
    end
    xgrid = range(1e-2, stop=b, length=100)
    ygrid = pdf.(xgrid);
    PyPlot.plot(xgrid, ygrid, label = "pdf(y)")
    PyPlot.vlines(median, 0, pdf(median), label = "median",  colors = "b")
    PyPlot.vlines(y0_true, 0, pdf(y0_true), label = "true value",  colors = "r")
    # plot CI
    if CI != nothing
        CI_id = (xgrid .> CI[1]) .* (xgrid .< CI[2])
        CI_xrange = vcat(CI1[1], xgrid[CI_id], CI1[2])
        CI_yrange = vcat(pdf(CI[1]), ygrid[CI_id], pdf(CI[2]))
        PyPlot.fill_between(CI_xrange, 0, CI_yrange, alpha = 0.3, label = "95% confidence interval")
    end
    PyPlot.legend(fontsize=8)
    PyPlot.grid()
    PyPlot.title(mytitle, fontsize=10)
#     PyPlot.savefig("Figure/Prediction.pdf")
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
function template_positions(v::Visualization, position)

end

function template_scales(v::Visualization, )

end

function template_colors(v::Visualization)

end

function build_shader(v::Visualization)

end

function fragment_output(target::Screen)
    "fragment_output.frag"
end


text:
export visualize_pressure, visualize_displacement,
    visualize_stress,
    visualize_von_mises_stress, visualize_saturation,
    visualize_potential, visualize_scalar_on_gauss_points,
    visualize_scalar_on_fem_points, visualize_scalar_on_fvm_points

"""
    visualize_pressure(U::Array{Float64, 2}, m::Int64, n::Int64, h::Float64)

Visualizes pressure. `U` is the solution vector.
"""
function visualize_pressure(U::Array{Float64, 2}, m::Int64, n::Int64, h::Float64)
    # fig, ax = subplots()
    close("all")
    x = (1:m)*h
    y = (1:n)*h
    vmin = mean(U) - 2std(U)
    vmax = mean(U) + 2std(U)
    U = reshape(U, size(U,1), m, n)
    U = permutedims(U, [1,3,2])
    ln = pcolormesh(x, y, U[1,:,:], vmin=vmin, vmax=vmax,rasterized=true)
    t = title("snapshot = 000")
    colorbar()
    axis("scaled")
    xlabel("x")
    ylabel("y")
    levels = LinRange(vmin, vmax, 10) |> Array
    c = contour(x, y, U[1,:,:], levels, cmap="jet", vmin=vmin, vmax=vmax)
    gca().invert_yaxis()
    function update(frame)
        gca().clear()
        k = string(frame-1)
        k = repeat("0", 3-length(k))*k
        title("snapshot = $k")
        pcolormesh(x, y, U[frame,:,:], vmin=vmin, vmax=vmax,rasterized=true)
        contour(x, y, U[frame,:,:], levels, cmap="jet", vmin=vmin, vmax=vmax)
        xlabel("x")
        ylabel("y")
        gca().invert_yaxis()
    end
    anim = animate(update, 1:size(U,1))
end

function visualize_pressure(U::Array{Float64,1}, m::Int64, n::Int64, h::Float64)
    # fig, ax = subplots()
    close("all")
    x = (1:m)*h
    y = (1:n)*h
    vmin = mean(U) - 2std(U)
    vmax = mean(U) + 2std(U)
    U = reshape(U,  m, n)
    U = Array(U')
    ln = pcolormesh(x, y, U, vmin=vmin, vmax=vmax,rasterized=true)
    colorbar()
    axis("scaled")
    xlabel("x")
    ylabel("y")
    levels = LinRange(vmin, vmax, 10) |> Array
    c = contour(x, y, U[:,:], levels, cmap="jet", vmin=vmin, vmax=vmax)
    gca().invert_yaxis()
end

@doc raw"""
    visualize_stress(K::Array{Float64, 2}, U::Array{Float64, 2}, m::Int64, n::Int64, h::Float64; name::String="")

Visualizes displacement. `U` is the solution vector, `K` is the elasticity matrix ($3\times 3$).
"""
function visualize_stress(K::Array{Float64, 2}, U::Array{Float64, 1}, m::Int64, n::Int64, h::Float64; name::String="")
    close("all")
    NT = size(U,1)
    x1 = LinRange(0.5h,m*h,m)|>collect
    y1 = LinRange(0.5h,n*h,n)|>collect
    S = zeros(n, m)
    s = compute_von_mises_stress_term(K, U, m, n, h)
    for i = 1:m
        for j = 1:n
            S[j,i] = sum(s[4*(m*(j-1)+i-1)+1:4*(m*(j-1)+i)])/4.0
        end
    end
    μ = mean(S); σ = std(S)
    vmin = μ - 2σ
    vmax = μ + 2σ

    levels = LinRange(vmin, vmax, 10) |> Array
    pcolormesh(x1,y1,S[:,:], vmin=vmin,vmax=vmax,rasterized=true)
    colorbar()
    contour(x1, y1, S[:,:], levels, cmap="jet")
    axis("scaled")
    xlabel("x")
    ylabel("y")
    gca().invert_yaxis()
end

function visualize_stress(K::Array{Float64, 2}, U::Array{Float64, 1}, m::Int64, n::Int64, h::Float64; name::String="")
    close("all")
    NT = size(U,1)
    x1 = LinRange(0.5h,m*h,m)|>collect
    y1 = LinRange(0.5h,n*h,n)|>collect
    S = zeros(NT, n, m)
    for k = 1:NT
        s = compute_von_mises_stress_term(K, U, m, n, h)
        for i = 1:m
            for j = 1:n
                S[k,j,i] = sum(s[4*(m*(j-1)+i-1)+1:4*(m*(j-1)+i)])/4.0
            end
        end
    end
    μ = mean(S); σ = std(S)
    vmin = μ - 2σ
    vmax = μ + 2σ

    levels = LinRange(vmin, vmax, 10) |> Array
    pcolormesh(x1,y1,S[1,:,:], vmin=vmin,vmax=vmax,rasterized=true)
    colorbar()
    contour(x1, y1, S[1,:,:], levels, cmap="jet")
    axis("scaled")
    xlabel("x")
    ylabel("y")
    gca().invert_yaxis()
end

"""
    visualize_stress(Se::Array{Float64, 2}, m::Int64, n::Int64, h::Float64; name::String="")

Visualizes the Von Mises stress. `Se` is the Von Mises at the cell center.
"""
function visualize_stress(Se::Array{Float64, 2}, m::Int64, n::Int64, h::Float64; name::String="", kwargs...)
    NT = size(Se,2)

    x1 = LinRange(0.5h,m*h,m)|>collect
    y1 = LinRange(0.5h,n*h,n)|>collect
    X1, Y1 = np.meshgrid(x1,y1)

    if size(Se, 1)==4*m*n
        Sep = zeros(m*n, NT)
        for k = 1:NT
            Sep[:,k] = mean([Se[1:4:end, k], Se[2:4:end, k], Se[3:4:end, k], Se[4:4:end, k]])
        end
        Se = Sep
    end

    S = zeros(size(Se, 2), n, m)
    for k = 1:size(Se, 2)
        S[k,:,:] = reshape(Se[:,k], m, n)'
    end
    μ = mean(S); σ = std(S)
    vmin = μ - 2σ
    vmax = μ + 2σ

    x = (1:m)*h
    y = (1:n)*h
    close("all")
    ln = pcolormesh(x, y, S[1,:,:], vmin= vmin, vmax=vmax)

    colorbar()
    # c = contour(φ[1,:,:], 10, cmap="jet", vmin=vmin,vmax=vmax)
    t = title("t = 0")
    axis("scaled")
    xlabel("x")
    ylabel("y")
    gca().invert_yaxis()
    function update(i)
        gca().clear()
        # t.set_text("t = $(round(frame * Δt, digits=3))")
        # ln.set_array(φ[frame,:,:]'[:])
        # c.set_array(φ[frame,:,:]'[:])
        ln = gca().pcolormesh(x, y, S[i,:,:], vmin= vmin, vmax=vmax)
        c = gca().contour(x, y, S[i,:,:], 10, cmap="jet", vmin=vmin,vmax=vmax)
        xlabel("x")
        ylabel("y")

        k = string(i-1)
        k = repeat("0", 3-length(k))*k
        title("snapshot = $k")
        gca().invert_yaxis()
    end
    anim = animate(update, 1:size(S,1))
end

"""
    visualize_von_mises_stress(Se::Array{Float64, 2}, m::Int64, n::Int64, h::Float64; name::String="")

Visualizes the Von Mises stress.
"""
function visualize_von_mises_stress(Se::Array{Float64, 3}, m::Int64, n::Int64, h::Float64; name::String="", kwargs...)
    S = zeros(size(Se, 1), m*n)
    for i = 1:size(Se, 1)
        S[i,:] = compute_von_mises_stress_term(Se[i,:,:], m, n, h)
    end
    visualize_stress(S'|>Array, m, n, h; name = name, kwargs...)
end

function visualize_von_mises_stress(Se::Array{Float64, 2}, m::Int64, n::Int64, h::Float64)
    S = compute_von_mises_stress_term(Se, m, n, h)

    x1 = LinRange(0.5h,m*h,m)|>collect
    y1 = LinRange(0.5h,n*h,n)|>collect
    X1, Y1 = np.meshgrid(x1,y1)

    μ = mean(S); σ = std(S)
    vmin = μ - 2σ
    vmax = μ + 2σ

    S = reshape(S, m, n)'
    x = (1:m)*h
    y = (1:n)*h
    close("all")
    ln = pcolormesh(x, y, S, vmin= vmin, vmax=vmax, rasterized=true)

    colorbar()
    c = contour(x, y, S, cmap="jet")
    axis("scaled")
    xlabel("x")
    ylabel("y")
    gca().invert_yaxis()

end

function visualize_saturation(s2::Array{Float64,3}, m::Int64, n::Int64, h::Float64)
    # fig, ax = subplots()
    close("all")
    x = (1:m)*h
    y = (1:n)*h
    ln = pcolormesh(x, y, s2[1,:,:], vmin=0.0, vmax=1.0)
    t = title("snapshot = 000")
    colorbar()
    axis("scaled")
    xlabel("x")
    ylabel("y")
    function update(frame)
        k = string(frame-1)
        k = repeat("0", 3-length(k))*k
        t.set_text("snapshot = $k")
        ln.set_array(s2[frame,1:end-1,1:end-1]'[:])
    end
    anim = animate(update, 1:size(s2,1))
end

function visualize_saturation(s2::Array{Float64,2}, m::Int64, n::Int64, h::Float64)
    # fig, ax = subplots()
    close("all")
    x = (1:m)*h
    y = (1:n)*h
    ln = pcolormesh(x, y, s2, vmin=0.0, vmax=1.0)
    colorbar()
    axis("scaled")
    xlabel("x")
    ylabel("y")
end

@doc raw"""
    visualize_scalar_on_fvm_points(φ::Array{Float64, 3}, m::Int64, n::Int64, h::Float64;
    vmin::Union{Real, Missing} = missing, vmax::Union{Real, Missing} = missing)

Generates scattered potential animation for the potential $\phi\in \mathbb{R}^{(NT+1)\times n \times m}$.
"""
function visualize_scalar_on_fvm_points(φ::Array{Float64, 3}, m::Int64, n::Int64, h::Float64;
    vmin::Union{Real, Missing} = missing, vmax::Union{Real, Missing} = missing)
    m_ = mean(φ)
    s = std(φ)
    close("all")
    vmin = coalesce(vmin, m_ - 2s)
    vmax = coalesce(vmax, m_ + 2s)
    x = (1:m)*h .- 0.5h
    y = (1:n)*h .- 0.5h
    ln = pcolormesh(x, y, φ[1,:,:], vmin= vmin, vmax=vmax)
    colorbar()
    # c = contour(φ[1,:,:], 10, cmap="jet", vmin=vmin,vmax=vmax)
    t = title("t = 0")
    axis("scaled")
    xlabel("x")
    ylabel("y")
    gca().invert_yaxis()
    function update(i)
        gca().clear()
        # t.set_text("t = $(round(frame * Δt, digits=3))")
        # ln.set_array(φ[frame,:,:]'[:])
        # c.set_array(φ[frame,:,:]'[:])
        ln = gca().pcolormesh(x, y, φ[i,:,:], vmin= vmin, vmax=vmax)
        c = gca().contour(x, y, φ[i,:,:], 10, cmap="jet", vmin=vmin,vmax=vmax)
        xlabel("x")
        ylabel("y")
        k = string(i-1)
        k = repeat("0", 3-length(k))*k
        title("snapshot = $k")
        gca().invert_yaxis()
    end
    anim = animate(update, 1:size(φ,1))
end

function visualize_scalar_on_fvm_points(φ::Array{Float64, 2}, m::Int64, n::Int64, h::Float64;
        vmin::Union{Real, Missing} = missing, vmax::Union{Real, Missing} = missing)
    R = zeros(size(φ,1), n, m)
    for i = 1:size(φ,1)
        R[i,:,:] = reshape(φ[i,:], m, n)'
    end
    visualize_scalar_on_fvm_points(R, m, n, h, vmin = vmin, vmax=vmax)
end

function visualize_scalar_on_fvm_points(φ::Array{Float64, 1}, m::Int64, n::Int64, h::Float64;
    vmin::Union{Real, Missing} = missing, vmax::Union{Real, Missing} = missing)
    φ = reshape(φ, m, n)'|>Array
    m_ = mean(φ)
    s = std(φ)
    vmin = coalesce(vmin, m_ - 2s)
    vmax = coalesce(vmax, m_ + 2s)
    x = (1:m)*h .- 0.5h
    y = (1:n)*h .- 0.5h
    ln = pcolormesh(x, y, φ, vmin= vmin, vmax=vmax)
    colorbar()
    # c = contour(φ[1,:,:], 10, cmap="jet", vmin=vmin,vmax=vmax)
    axis("scaled")
    xlabel("x")
    ylabel("y")
    ln = gca().pcolormesh(x, y, φ, vmin= vmin, vmax=vmax)
    c = gca().contour(x, y, φ, 10, cmap="jet", vmin=vmin,vmax=vmax)
    gca().invert_yaxis()
end

visualize_potential = visualize_scalar_on_fvm_points

@doc raw"""
    visualize_displacement(u::Array{Float64, 2}, m::Int64, n::Int64, h::Float64)

Generates scattered plot animation for displacement $u\in \mathbb{R}^{(NT+1)\times 2(m+1)(n+1)}$.
"""
function visualize_displacement(u::Array{Float64, 2}, m::Int64, n::Int64, h::Float64)
    X = zeros(m+1, n+1)
    Y = zeros(m+1, n+1)
    for i = 1:m+1
        for j = 1:n+1
            X[i, j] = (i-1)*h
            Y[i, j] = (j-1)*h
        end
    end
    function disp(u)
        U1 = reshape(u[1:(m+1)*(n+1)], m+1, n+1)
        U2 = reshape(u[(m+1)*(n+1)+1:2(m+1)*(n+1)], m+1, n+1)
        U1 = X + U1
        U2 = Y + U2
        U1, U2
    end
    close("all")
    U1, U2 = disp(u[1,:])
    s = scatter(U1[:], U2[:], s=5)
    xmin = minimum(u[:,1:(m+1)*(n+1)])
    xmax = maximum(u[:,1:(m+1)*(n+1)]) + m*h
    ymin = minimum(u[:,(m+1)*(n+1)+1:2(m+1)*(n+1)])
    ymax = maximum(u[:,(m+1)*(n+1)+1:2(m+1)*(n+1)]) + n*h
    t = title("snapshot = 000")
    xlabel("x")
    ylabel("y")
    axis("equal")
    xlim(xmin.-h, xmax.+h)
    ylim(ymin.-h, ymax.+h)
    gca().invert_yaxis()
    function update(i)
        U1, U2 = disp(u[i,:])
        s.set_offsets([U1[:] U2[:]])

        k = string(i-1)
        k = repeat("0", 3-length(k))*k
        t.set_text("snapshot = $k")
        xlim(xmin.-h, xmax.+h)
        ylim(ymin.-h, ymax.+h)
        gca().invert_yaxis()
    end
    animate(update, 1:size(u,1))
end

function visualize_displacement(u::Array{Float64, 1}, m::Int64, n::Int64, h::Float64)
    X = zeros(m+1, n+1)
    Y = zeros(m+1, n+1)
    for i = 1:m+1
        for j = 1:n+1
            X[i, j] = (i-1)*h
            Y[i, j] = (j-1)*h
        end
    end
    function disp(u)
        U1 = reshape(u[1:(m+1)*(n+1)], m+1, n+1)
        U2 = reshape(u[(m+1)*(n+1)+1:2(m+1)*(n+1)], m+1, n+1)
        U1 = X + U1
        U2 = Y + U2
        U1, U2
    end
    close("all")
    U1, U2 = disp(u)
    s = scatter(U1[:], U2[:], s=5)
    xmin = minimum(u[1:(m+1)*(n+1)])
    xmax = maximum(u[1:(m+1)*(n+1)]) + m*h
    ymin = minimum(u[(m+1)*(n+1)+1:2(m+1)*(n+1)])
    ymax = maximum(u[(m+1)*(n+1)+1:2(m+1)*(n+1)]) + n*h
    xlim(xmin.-h, xmax.+h)
    ylim(ymin.-h, ymax.+h)
    gca().invert_yaxis()
    xlabel("x")
    ylabel("y")
    axis("equal")
end

@doc raw"""
    visualize_scalar_on_gauss_points(u::Array{Float64,1}, m::Int64, n::Int64, h::Float64, args...;kwargs...)

Visualizes the scalar `u` using pcolormesh. Here `u` is a length $4mn$ vector and the values are defined on the Gauss points
"""
function visualize_scalar_on_gauss_points(u::Array{Float64,1}, m::Int64, n::Int64, h::Float64, args...;kwargs...)
    # close("all")
    z = zeros(2m, 2n)
    x = zeros(2m)
    y = zeros(2n)

    for i = 1:m
        x[2*(i-1)+1] = pts[1] * h + (i-1)*h
        x[2*(i-1)+2] = pts[2] * h + (i-1)*h
    end

    for j = 1:n
        y[2*(j-1)+1] = pts[1] * h + (j-1)*h
        y[2*(j-1)+2] = pts[2] * h + (j-1)*h
    end

    for j = 1:n
    end

    for i = 1:m
        for j = 1:n
            idx = (j-1)*m+i
            z[2(i-1)+1, 2(j-1)+1] = u[4*(idx-1)+1]
            z[2(i-1)+2, 2(j-1)+1] = u[4*(idx-1)+2]
            z[2(i-1)+1, 2(j-1)+2] = u[4*(idx-1)+3]
            z[2(i-1)+2, 2(j-1)+2] = u[4*(idx-1)+4]
        end
    end

    vmin = mean(z) - 2std(z)
    vmax = mean(z) + 2std(z)
    pcolormesh(x, y, z', vmin=vmin, vmax=vmax,rasterized=true, args...; kwargs...)
    colorbar()
    axis("scaled")
    xlabel("x")
    ylabel("y")
    levels = LinRange(vmin, vmax, 10) |> Array
    c = contour(x, y, z', levels, cmap="jet", vmin=vmin, vmax=vmax)
    gca().invert_yaxis()
    return x, y, z
end

function visualize_scalar_on_gauss_points(u::Array{Float64,2}, m::Int64, n::Int64, h::Float64, args...;kwargs...)
    # close("all")
    z = zeros(2m, 2n)
    x = zeros(2m)
    y = zeros(2n)

    for i = 1:m
        x[2*(i-1)+1] = pts[1] * h + (i-1)*h
        x[2*(i-1)+2] = pts[2] * h + (i-1)*h
    end

    for j = 1:n
        y[2*(j-1)+1] = pts[1] * h + (j-1)*h
        y[2*(j-1)+2] = pts[2] * h + (j-1)*h
    end

    NT = size(u,1)-1
    z = zeros(NT+1, 2m, 2n)
    for k = 1:NT+1
        for i = 1:m
            for j = 1:n
                idx = (j-1)*m+i
                z[k, 2(i-1)+1, 2(j-1)+1] = u[k, 4*(idx-1)+1]
                z[k, 2(i-1)+2, 2(j-1)+1] = u[k, 4*(idx-1)+2]
                z[k, 2(i-1)+1, 2(j-1)+2] = u[k, 4*(idx-1)+3]
                z[k, 2(i-1)+2, 2(j-1)+2] = u[k, 4*(idx-1)+4]
            end
        end
    end

    vmin = mean(z) - 2std(z)
    vmax = mean(z) + 2std(z)
    pcolormesh(x, y, z[1,:,:]', vmin=vmin, vmax=vmax,rasterized=true, args...; kwargs...)
    t = title("snapshot = 000")
    colorbar()
    axis("scaled")
    xlabel("x")
    ylabel("y")
    levels = LinRange(vmin, vmax, 10) |> Array
    c = contour(x, y, z[1,:,:]', levels, cmap="jet", vmin=vmin, vmax=vmax)
    gca().invert_yaxis()

    function update(i)
        gca().clear()
        pcolormesh(x, y, z[i,:,:]', vmin=vmin, vmax=vmax,rasterized=true, args...; kwargs...)
        k = string(i-1)
        k = repeat("0", 3-length(k))*k
        title("snapshot = $k")
        xlim(-h, m*h+h)
        ylim(-h, n*h+h)
        c = contour(x, y, z[i,:,:]', levels, cmap="jet", vmin=vmin, vmax=vmax)
        gca().invert_yaxis()
    end
    animate(update, 2:size(u,1))
end

@doc raw"""
    visualize_scalar_on_fem_points(u::Array{Float64,1}, m::Int64, n::Int64, h::Float64, args...;kwargs...)

Visualizes the scalar `u` using pcolormesh. Here `u` is a length $(m+1)(n+1)$ vector and the values are defined on the FEM points
"""
function visualize_scalar_on_fem_points(u::Array{Float64,1}, m::Int64, n::Int64, h::Float64, args...;kwargs...)
    x = Array((0:m)*h)
    y = Array((0:n)*h)
    z = reshape(u, m+1, n+1)

    vmin = mean(z) - 2std(z)
    vmax = mean(z) + 2std(z)
    pcolormesh(x, y, z', vmin=vmin, vmax=vmax,rasterized=true, args...; kwargs...)
    colorbar()
    axis("scaled")
    xlabel("x")
    ylabel("y")
    levels = LinRange(vmin, vmax, 10) |> Array
    c = contour(x, y, z', levels, cmap="jet", vmin=vmin, vmax=vmax)
    gca().invert_yaxis()
    return x, y, z
end

@doc raw"""
    visualize_scalar_on_fem_points(u::Array{Float64,2}, m::Int64, n::Int64, h::Float64, args...;kwargs...)

Visualizes the scalar `u` using pcolormesh. Here `u` is a matrix of size $NT \times (m+1)(n+1)$ ($NT$ is the number of time steps) and the values are defined on the FEM points.
"""
function visualize_scalar_on_fem_points(u::Array{Float64,2}, m::Int64,
        n::Int64, h::Float64, args...;kwargs...)
    z = zeros(size(u,1), m+1, n+1)
    for i = 1:size(u,1)
        z[i,:,:] = reshape(u[i,:], m+1, n+1)
    end
    x = Array((0:m)*h)
    y = Array((0:n)*h)
    vmin = mean(z) - 2std(z)
    vmax = mean(z) + 2std(z)
    pcolormesh(x, y, z[1,:,:]', vmin=vmin, vmax=vmax,rasterized=true, args...; kwargs...)
    t = title("snapshot = 000")
    colorbar()
    axis("scaled")
    xlabel("x")
    ylabel("y")
    levels = LinRange(vmin, vmax, 10) |> Array
    c = contour(x, y, z[1,:,:]', levels, cmap="jet", vmin=vmin, vmax=vmax)
    gca().invert_yaxis()

    function update(i)
        gca().clear()
        pcolormesh(x, y, z[i,:,:]', vmin=vmin, vmax=vmax,rasterized=true, args...; kwargs...)
        k = string(i-1)
        k = repeat("0", 3-length(k))*k
        title("snapshot = $k")
        xlim(-h, m*h+h)
        ylim(-h, n*h+h)
        c = contour(x, y, z[i,:,:]', levels, cmap="jet", vmin=vmin, vmax=vmax)
        gca().invert_yaxis()
    end
    animate(update, 2:size(u,1))
end


text:
# syntax: proto2
using Compat
using ProtoBuf
import ProtoBuf.meta

mutable struct UserInformationProto <: ProtoType
    effectiveUser::AbstractString
    realUser::AbstractString
    UserInformationProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct UserInformationProto

mutable struct IpcConnectionContextProto <: ProtoType
    userInfo::UserInformationProto
    protocol::AbstractString
    IpcConnectionContextProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct IpcConnectionContextProto
const __fnum_IpcConnectionContextProto = Int[2,3]
meta(t::Type{IpcConnectionContextProto}) = meta(t, ProtoBuf.DEF_REQ, __fnum_IpcConnectionContextProto, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

export UserInformationProto, IpcConnectionContextProto


text:
#=
# A reference element.
# Notation follows that used in Nodal Discontinuous Galerkin Methods
# by Hesthaven and Warburton.
#  https://link.springer.com/book/10.1007/978-0-387-72067-8
#
# Can be used for 1, 2, 3, 4 dimensions
=#
include("jacobi_polynomial.jl");
include("refel_nodes.jl");
include("refel_triangle.jl");
include("refel_tet.jl");

mutable struct Refel
    dim::Int                # Dimension
    N::Int                  # Order of polynomials
    Np::Int                 # Number of nodes
    Nqp::Int                # Number of quadrature points
    Nfaces::Int             # Number of faces
    Nfp::Array{Int,1}       # Number of nodes for each face

    ######################################
    # Volume nodes and quadrature matrices
    ######################################
    r1d::Array{Float64}     # Node coordinates in 1D
    r::Array{Float64}       # dim-dim Node coordinates

    wr1d::Array{Float64}    # r1d gll Quadrature weights
    wr::Array{Float64}      # r gll Quadrature weights

    g1d::Array{Float64}     # 1D Gauss points
    wg1d::Array{Float64}    # 1D Gauss weights

    g::Array{Float64}       # dim-dim Gauss points
    wg::Array{Float64}      # dim-dim Gauss weights

    V::Array{Float64}       # basis at r
    gradV::Array{Float64}   # grad of basis at r
    invV::Array{Float64}    # Inverse V

    Vg::Array{Float64}      # basis at Gauss
    gradVg::Array{Float64}  # grad of basis at g
    invVg::Array{Float64}   # Inverse Vg

    Dr::Array{Float64}      # Differentiation matrix for r
    Ds::Array{Float64}      # Differentiation matrix for s
    Dt::Array{Float64}      # Differentiation matrix for t
    Dg::Array{Float64}      # Differentiation matrix for g

    # Useful quadrature matrices for the volume integrals
    Q1d::Array{Float64}     # 1D quadrature matrix: like Vg*invV
    Q::Array{Float64}       # dim-dim quadrature matrix
    Qr::Array{Float64}      # quad of derivative matrix: like gradVg*invV
    Qs::Array{Float64}      #
    Qt::Array{Float64}      #

    Ddr::Array{Float64}      # Derivatives at the elemental nodes, not quadrature nodes
    Dds::Array{Float64}      #
    Ddt::Array{Float64}      #

    #######################################
    # Surface nodes and quadrature matrices
    #######################################
    face2local::Array{Array{Int}}       # maps face nodes to local indices

    surf_r::Array{Array{Float64}}       # surface node coordinates
    surf_wr::Array{Array{Float64}}      # surface gll weights

    surf_g::Array{Array{Float64}}       # surface Gauss points
    surf_wg::Array{Array{Float64}}      # surface Gauss weights

    surf_V::Array{Array{Float64}}       # basis at surf_r
    surf_gradV::Array{Array{Float64}}   # grad of basis at surf_r
    surf_invV::Array{Array{Float64}}    # Inverse surf_V

    surf_Vg::Array{Array{Float64}}      # basis at surf_g
    surf_gradVg::Array{Array{Float64}}  # grad of basis at surf_g
    surf_invVg::Array{Array{Float64}}   # Inverse surf_Vg

    surf_Dr::Array{Array{Float64}}      # Differentiation matrix for surf_r
    surf_Ds::Array{Array{Float64}}      # Differentiation matrix for surf_s
    surf_Dt::Array{Array{Float64}}      # Differentiation matrix for surf_t
    surf_Dg::Array{Array{Float64}}      # Differentiation matrix for surf_g

    surf_Q::Array{Array{Float64}}       # quadrature matrix
    surf_Qr::Array{Array{Float64}}      # derivative quadrature matrix
    surf_Qs::Array{Array{Float64}}      #
    surf_Qt::Array{Array{Float64}}      #

    surf_Ddr::Array{Array{Float64}}     # Derivatives at the elemental nodes, not quadrature nodes
    surf_Dds::Array{Array{Float64}}     #
    surf_Ddt::Array{Array{Float64}}     #

    # Constructor needs at least this information
    Refel(dim, order, nnodes, nfaces, nfp) = new(
        dim,
        order,
        nnodes,
        -1,
        nfaces,
        nfp,
        [],[],
        [],[],
        [],[],
        [],[],
        [],[],[],
        [],[],[],
        [],[],[],[],
        [],[],[],[],[],
        [],[],[],
        [[]],
        [[]],[[]],
        [[]],[[]],
        [[]],[[]],[[]],
        [[]],[[]],[[]],
        [[]],[[]],[[]],[[]],
        [[]],[[]],[[]],[[]],
        [[]],[[]],[[]]
    )
end

import Base.copy
function copy(ref::Refel)
    newref = Refel(ref.dim, ref.N, ref.Np, ref.Nfaces, ref.Nfp);
    newref.Nqp = ref.Nqp;
    newref.r1d = copy(ref.r1d);
    newref.r = copy(ref.r);
    newref.wr1d = copy(ref.wr1d);
    newref.wr = copy(ref.wr);
    newref.g1d = copy(ref.g1d);
    newref.wg1d = copy(ref.wg1d);
    newref.g = copy(ref.g);
    newref.wg = copy(ref.wg);
    newref.V = copy(ref.V);
    newref.gradV = copy(ref.gradV);
    newref.invV = copy(ref.invV);
    newref.Vg = copy(ref.Vg);
    newref.gradVg = copy(ref.gradVg);
    newref.invVg = copy(ref.invVg);
    newref.Dr = copy(ref.Dr);
    newref.Ds = copy(ref.Ds);
    newref.Dt = copy(ref.Dt);
    newref.Dg = copy(ref.Dg);
    newref.Q1d = copy(ref.Q1d);
    newref.Q = copy(ref.Q);
    newref.Qr = copy(ref.Qr);
    newref.Qs = copy(ref.Qs);
    newref.Qt = copy(ref.Qt);
    newref.Ddr = copy(ref.Ddr);
    newref.Dds = copy(ref.Dds);
    newref.Ddt = copy(ref.Ddt);

    return newref;
end

function build_refel(dimension, order, nfaces, nodetype)
    # Check for errors
    if (dimension == 1 && nfaces != 2) || (dimension > 1 && nfaces < dimension+1)
        println("Error: buildRefel(dimension, order, nfaces, nodetype), check for valid parameters.");
        return nothing;
    end

    # Number of points determined by order element type
    if (dimension == 0)
        Np = 1; #point
        Nfp = []; #no faces

    elseif (dimension == 1) # line segment
        Np = order+1;
        Nfp = [1, 1]; # face points

    elseif (dimension == 2 && nfaces == 3) # triangle
        Np = (Int)((order+1)*(order+2)/2);
        Nfp = [order+1, order+1, order+1]; # face lines

    elseif (dimension == 2 && nfaces == 4) # quad
        Np = (Int)((order+1)*(order+1));
        Nfp = [order+1, order+1, order+1, order+1]; # face lines

    elseif (dimension == 3 && nfaces == 4)  # tet
        Np = (Int)((order+1)*(order+2)*(order+3)/6);
        M = (Int)((order+1)*(order+2)/2);
        Nfp = [M, M, M, M]; # face triangles

    elseif (dimension == 3 && nfaces == 6)  # hex
        Np = (Int)((order+1)*(order+1)*(order+1)); # hex
        M = (Int)((order+1)*(order+1));
        Nfp = [M, M, M, M, M, M]; # face quads

    elseif (dimension == 4)  # ??
        Np = (Int)((order+1)*(order+2)*(order+3)*(order+4)/24); # ??
        # TODO
    end

    refel = Refel(dimension, order, Np, nfaces, Nfp);

    # Get nodes on the reference element
    refel_nodes!(refel, nodetype);
    refel.Nqp = length(refel.wg);

    #  0D refels
    if dimension == 0
        refel.V = ones(1,1);
        refel.gradV = zeros(1,1);
        refel.invV = ones(1,1);
        refel.Vg = ones(1,1);
        refel.gradVg = zeros(1,1);
        refel.invVg = ones(1,1);
        refel.Dr = zeros(1,1);
        refel.Dg = zeros(1,1);
        refel.Q1d = zeros(1,1);
        refel.Q = ones(1,1);
        refel.Qr = zeros(1,1);
        refel.Qs = zeros(1,1);
        refel.Qt = zeros(1,1);
        refel.Ddr = zeros(1,1);
        refel.Dds = zeros(1,1);
        refel.Ddt = zeros(1,1);

        return refel;
    end

    if (dimension == 2 && nfaces == 3) # triangle
        refel = build_triangle_refel(refel);

    elseif (dimension == 3 && nfaces == 4) # tet
        refel = build_tetrahedron_refel(refel);

    else # line, quad, hex
        # Vandermonde matrix and grad,inv
        # Values of basis functions and derivs at points
        refel.V = zeros(order+1, order+1);
        refel.gradV = zeros(order+1, order+1);
        # Gauss versions
        refel.Vg = zeros(order+1, order+1);
        refel.gradVg = zeros(order+1, order+1);
        # surface versions
        refel.surf_V = fill!(Array{Array{Float64}}(undef, nfaces), []);
        refel.surf_gradV = fill!(Array{Array{Float64}}(undef, nfaces), []);
        refel.surf_Vg = fill!(Array{Array{Float64}}(undef, nfaces), []);
        refel.surf_gradVg = fill!(Array{Array{Float64}}(undef, nfaces), []);
        for fi=1:nfaces
            refel.surf_V[fi] = zeros(refel.Nfp[fi], Np);
            refel.surf_gradV[fi] = zeros(refel.Nfp[fi], Np);
            refel.surf_Vg[fi] = zeros(refel.Nfp[fi], Np);
            refel.surf_gradVg[fi] = zeros(refel.Nfp[fi], Np);
        end

        # nodal versions
        for i=1:refel.N+1
            refel.V[:,i] = jacobi_polynomial(refel.r1d, 0, 0, i-1);
        end
        for i=1:refel.N
            refel.gradV[:,i+1] = sqrt(i*(i+1)) .* jacobi_polynomial(refel.r1d, 1, 1, i-1);
        end
        refel.invV = inv(refel.V);

        # Gauss versions
        for i=1:refel.N+1
            refel.Vg[:,i] = jacobi_polynomial(refel.g1d, 0, 0, i-1);
        end
        for i=1:refel.N
            refel.gradVg[:,i+1] = sqrt(i*(i+1)) .* jacobi_polynomial(refel.g1d, 1, 1, i-1);
        end
        refel.invVg = inv(refel.Vg);

        # Differentiation matrices
        refel.Dr = refel.gradV*refel.invV;
        refel.Dg = refel.gradVg*refel.invV;

        refel.Q1d = refel.Vg*refel.invV;

        if dimension == 1
            # volume
            refel.Q = refel.Q1d;
            refel.Qr = refel.Dg;
            refel.Ddr = refel.Dr;
            #surface
            # r and g are the same for this case because there's only one point
            refel.surf_Q = [refel.V[[1],:] * refel.invV, refel.V[[Np],:] * refel.invV];
            refel.surf_Qr = [refel.gradV[[1],:] * refel.invV, refel.gradV[[Np],:] * refel.invV];
            refel.surf_Ddr = refel.surf_Qr;

        elseif dimension == 2
            # volume
            ident = Matrix(1.0*I,order+1,order+1);
            refel.Q = kron(refel.Q1d,refel.Q1d);
            refel.Qr = kron(refel.Q1d,refel.Dg);
            refel.Qs = kron(refel.Dg,refel.Q1d);
            refel.Ddr = kron(ident,refel.Dr);
            refel.Dds = kron(refel.Dr,ident);
            # surface
            refel.surf_Q = fill!(Array{Array{Float64}}(undef, nfaces), []);
            refel.surf_Qr = fill!(Array{Array{Float64}}(undef, nfaces), []);
            refel.surf_Qs = fill!(Array{Array{Float64}}(undef, nfaces), []);
            refel.surf_Ddr = fill!(Array{Array{Float64}}(undef, nfaces), []);
            refel.surf_Dds = fill!(Array{Array{Float64}}(undef, nfaces), []);
            ## Surfaces will use the full matrices rather than using tensor products
            fullinvV = kron(refel.invV, refel.invV);

            for fi=1:nfaces
                surf_gradVr = zeros(refel.Nfp[fi], Np);
                surf_gradVs = zeros(refel.Nfp[fi], Np);
                surf_gradVgr = zeros(refel.Nfp[fi], Np);
                surf_gradVgs = zeros(refel.Nfp[fi], Np);
                for ni=1:size(refel.surf_r[fi],1)
                    # nodal versions
                    for i=1:refel.N+1
                        for j=1:refel.N+1
                            ind = (j-1)*(refel.N+1) + i;
                            refel.surf_V[fi][ni,ind] = (jacobi_polynomial(refel.surf_r[fi][ni,1], 0, 0, i-1) .* jacobi_polynomial(refel.surf_r[fi][ni,2], 0, 0, j-1))[1];
                        end
                    end
                    for i=1:refel.N
                        for j=1:refel.N
                            ind = (j)*(refel.N) + i + 1;
                            surf_gradVr[ni,ind] = sqrt(i*(i+1)) * (jacobi_polynomial(refel.surf_r[fi][ni,1], 1, 1, i-1) .* jacobi_polynomial(refel.surf_r[fi][ni,2], 0, 0, j-1))[1];
                            surf_gradVs[ni,ind] = sqrt(j*(j+1)) * (jacobi_polynomial(refel.surf_r[fi][ni,1], 0, 0, i-1) .* jacobi_polynomial(refel.surf_r[fi][ni,2], 1, 1, j-1))[1];
                        end
                    end
                    # Gauss versions
                    for i=1:refel.N+1
                        for j=1:refel.N+1
                            ind = (j-1)*(refel.N+1) + i;
                            refel.surf_Vg[fi][ni,ind] = (jacobi_polynomial(refel.surf_g[fi][ni,1], 0, 0, i-1) .* jacobi_polynomial(refel.surf_g[fi][ni,2], 0, 0, j-1))[1];
                        end
                    end
                    for i=1:refel.N
                        for j=1:refel.N
                            ind = (j)*(refel.N) + i + 1;
                            surf_gradVgr[ni,ind] = sqrt(i*(i+1)) * (jacobi_polynomial(refel.surf_g[fi][ni,1], 1, 1, i-1) .* jacobi_polynomial(refel.surf_g[fi][ni,2], 0, 0, j-1))[1];
                            surf_gradVgs[ni,ind] = sqrt(j*(j+1)) * (jacobi_polynomial(refel.surf_g[fi][ni,1], 0, 0, i-1) .* jacobi_polynomial(refel.surf_g[fi][ni,2], 1, 1, j-1))[1];
                        end
                    end
                end

                refel.surf_Q[fi] = refel.surf_V[fi] * fullinvV;
                refel.surf_Qr[fi] = surf_gradVgr * fullinvV;
                refel.surf_Qs[fi] = surf_gradVgs * fullinvV;
                refel.surf_Ddr[fi] = surf_gradVr * fullinvV;
                refel.surf_Dds[fi] = surf_gradVs * fullinvV;
            end

        elseif dimension == 3
            ident = Matrix(1.0*I,order+1,order+1);
            refel.Q = kron(kron(refel.Q1d, refel.Q1d), refel.Q1d);
            refel.Qr = kron(kron(refel.Q1d, refel.Q1d), refel.Dg);
            refel.Qs = kron(kron(refel.Q1d, refel.Dg), refel.Q1d);
            refel.Qt = kron(kron(refel.Dg, refel.Q1d), refel.Q1d);
            refel.Ddr = kron(kron(ident, ident), refel.Dg);
            refel.Dds = kron(kron(ident, refel.Dg), ident);
            refel.Ddt = kron(kron(refel.Dg, ident), ident);

            # surface
            refel.surf_Q = fill!(Array{Array{Float64}}(undef, nfaces), []);
            refel.surf_Qr = fill!(Array{Array{Float64}}(undef, nfaces), []);
            refel.surf_Qs = fill!(Array{Array{Float64}}(undef, nfaces), []);
            refel.surf_Qt = fill!(Array{Array{Float64}}(undef, nfaces), []);
            refel.surf_Ddr = fill!(Array{Array{Float64}}(undef, nfaces), []);
            refel.surf_Dds = fill!(Array{Array{Float64}}(undef, nfaces), []);
            refel.surf_Ddt = fill!(Array{Array{Float64}}(undef, nfaces), []);
            ## Surfaces will use the full matrices rather than using tensor products
            fullinvV = kron(refel.invV, kron(refel.invV, refel.invV));

            for fi=1:nfaces
                surf_gradVr = zeros(refel.Nfp[fi], Np);
                surf_gradVs = zeros(refel.Nfp[fi], Np);
                surf_gradVt = zeros(refel.Nfp[fi], Np);
                surf_gradVgr = zeros(refel.Nfp[fi], Np);
                surf_gradVgs = zeros(refel.Nfp[fi], Np);
                surf_gradVgt = zeros(refel.Nfp[fi], Np);
                for ni=1:size(refel.surf_r[fi],1)
                    # nodal versions
                    for i=1:refel.N+1
                        for j=1:refel.N+1
                            for k=1:refel.N+1
                                ind = (k-1)*(refel.N+1)*(refel.N+1) + (j-1)*(refel.N+1) + i;
                                refel.surf_V[fi][ni,ind] = (jacobi_polynomial(refel.surf_r[fi][ni,1], 0, 0, i-1) .* jacobi_polynomial(refel.surf_r[fi][ni,2], 0, 0, j-1) .* jacobi_polynomial(refel.surf_r[fi][ni,3], 0, 0, k-1))[1];
                            end
                        end
                    end
                    for i=1:refel.N
                        for j=1:refel.N
                            for k=1:refel.N
                                ind = (k)*(refel.N)*(refel.N) + (j)*(refel.N) + i + 1;
                                surf_gradVr[ni,ind] = sqrt(i*(i+1)) * (jacobi_polynomial(refel.surf_r[fi][ni,1], 1, 1, i-1) .* jacobi_polynomial(refel.surf_r[fi][ni,2], 0, 0, j-1) .* jacobi_polynomial(refel.surf_r[fi][ni,3], 0, 0, k-1))[1];
                                surf_gradVs[ni,ind] = sqrt(j*(j+1)) * (jacobi_polynomial(refel.surf_r[fi][ni,1], 0, 0, i-1) .* jacobi_polynomial(refel.surf_r[fi][ni,2], 1, 1, j-1) .* jacobi_polynomial(refel.surf_r[fi][ni,3], 0, 0, k-1))[1];
                                surf_gradVt[ni,ind] = sqrt(k*(k+1)) * (jacobi_polynomial(refel.surf_r[fi][ni,1], 0, 0, i-1) .* jacobi_polynomial(refel.surf_r[fi][ni,2], 0, 0, j-1) .* jacobi_polynomial(refel.surf_r[fi][ni,3], 1, 1, k-1))[1];
                            end
                        end
                    end
                    # Gauss versions
                    for i=1:refel.N+1
                        for j=1:refel.N+1
                            for k=1:refel.N+1
                                ind = (k-1)*(refel.N+1)*(refel.N+1) + (j-1)*(refel.N+1) + i;
                                refel.surf_Vg[fi][ni,ind] = (jacobi_polynomial(refel.surf_g[fi][ni,1], 0, 0, i-1) .* jacobi_polynomial(refel.surf_g[fi][ni,2], 0, 0, j-1) .* jacobi_polynomial(refel.surf_g[fi][ni,3], 0, 0, k-1))[1];
                            end
                        end
                    end
                    for i=1:refel.N
                        for j=1:refel.N
                            for k=1:refel.N
                                ind = (k)*(refel.N)*(refel.N) + (j)*(refel.N) + i + 1;
                                surf_gradVgr[ni,ind] = sqrt(i*(i+1)) * (jacobi_polynomial(refel.surf_g[fi][ni,1], 1, 1, i-1) .* jacobi_polynomial(refel.surf_g[fi][ni,2], 0, 0, j-1) .* jacobi_polynomial(refel.surf_g[fi][ni,3], 0, 0, k-1))[1];
                                surf_gradVgs[ni,ind] = sqrt(j*(j+1)) * (jacobi_polynomial(refel.surf_g[fi][ni,1], 0, 0, i-1) .* jacobi_polynomial(refel.surf_g[fi][ni,2], 1, 1, j-1) .* jacobi_polynomial(refel.surf_g[fi][ni,3], 0, 0, k-1))[1];
                                surf_gradVgt[ni,ind] = sqrt(k*(k+1)) * (jacobi_polynomial(refel.surf_g[fi][ni,1], 0, 0, i-1) .* jacobi_polynomial(refel.surf_g[fi][ni,2], 0, 0, j-1) .* jacobi_polynomial(refel.surf_g[fi][ni,3], 1, 1, k-1))[1];
                            end
                        end
                    end
                end

                refel.surf_Q[fi] = refel.surf_V[fi] * fullinvV;
                refel.surf_Qr[fi] = surf_gradVgr * fullinvV;
                refel.surf_Qs[fi] = surf_gradVgs * fullinvV;
                refel.surf_Qt[fi] = surf_gradVgt * fullinvV;
                refel.surf_Ddr[fi] = surf_gradVr * fullinvV;
                refel.surf_Dds[fi] = surf_gradVs * fullinvV;
                refel.surf_Ddt[fi] = surf_gradVt * fullinvV;
            end
        end
    end

    return refel;
end

function custom_quadrature_refel(oldrefel, nodes, weights)
    # The supplied nodes will replace the gauss quadrature nodes and weights r.g and r.wg
    # NOTE: Elemental node and weight arrays will not be changed, so quadrature must be set to GAUSS.
    # NOTE: Only quadrature matrices are changed. Nothing else.
    refel = copy(oldrefel);
    refel.g = copy(nodes);
    refel.wg = copy(weights);

    # Vandermonde matrix and grad,inv for r will not change

    # Gauss versions will change
    Nn = size(nodes,2);
    refel.Vg = zeros(Nn, refel.N+1);
    refel.gradVg = zeros(Nn, refel.N+1);
    if refel.dim == 1
        for i=1:refel.N+1
            refel.Vg[:,i] = jacobi_polynomial(nodes[1,:], 0, 0, i-1);
        end
        for i=1:refel.N
            refel.gradVg[:,i+1] = sqrt(i*(i+1)) .* jacobi_polynomial(nodes[1,:], 1, 1, i-1);
        end

        # Differentiation matrices
        refel.Dg = refel.gradVg*refel.invV;
        refel.Q1d = refel.Vg*refel.invV;

        # Quadrature matrices
        refel.Q = refel.Q1d;
        refel.Qr = refel.Dg;
        refel.Ddr = refel.Dr;

    elseif refel.dim == 2
        # Now it's a little trickier, so forget the efficiency and do it carefully.
        Nn = size(nodes,2);
        tmp1 = zeros(Nn, refel.N+1);
        tmp2 = zeros(Nn, refel.N+1);
        Dtmp1 = zeros(Nn, refel.N+1);
        Dtmp2 = zeros(Nn, refel.N+1);

        for i=1:refel.N+1
            tmp1[:,i] = jacobi_polynomial(nodes[1,:], 0, 0, i-1);
            tmp2[:,i] = jacobi_polynomial(nodes[2,:], 0, 0, i-1);
        end
        for i=1:refel.N
            Dtmp1[:,i+1] = sqrt(i*(i+1)) .* jacobi_polynomial(nodes[1,:], 1, 1, i-1);
            Dtmp2[:,i+1] = sqrt(i*(i+1)) .* jacobi_polynomial(nodes[2,:], 1, 1, i-1);
        end

        # Use kron(Vg*invV, Vg,invV) = kron(Vg,Vg)*kron(invV,invV)
        VgXVg = zeros(Nn, (refel.N+1)*(refel.N+1));
        VgXDg = zeros(Nn, (refel.N+1)*(refel.N+1));
        DgXVg = zeros(Nn, (refel.N+1)*(refel.N+1));
        for i=1:Nn
            for j=1:(refel.N+1)
                for k=1:(refel.N+1)
                    ind = (k-1)*(refel.N+1) + j;
                    VgXVg[i,ind] = tmp1[i,j]*tmp2[i,k];
                    VgXDg[i,ind] = Dtmp1[i,j]*tmp2[i,k];
                    DgXVg[i,ind] = tmp1[i,j]*Dtmp2[i,k];
                end
            end
        end
        viXvi = kron(refel.invV, refel.invV);

        refel.Q = VgXVg * viXvi;
        refel.Qr = VgXDg * viXvi;
        refel.Qs = DgXVg * viXvi;

    elseif refel.dim == 3
        Nn = size(nodes,2);
        tmp1 = zeros(Nn, refel.N+1);
        tmp2 = zeros(Nn, refel.N+1);
        tmp3 = zeros(Nn, refel.N+1);
        Dtmp1 = zeros(Nn, refel.N+1);
        Dtmp2 = zeros(Nn, refel.N+1);
        Dtmp3 = zeros(Nn, refel.N+1);

        for i=1:refel.N+1
            tmp1[:,i] = jacobi_polynomial(nodes[1,:], 0, 0, i-1);
            tmp2[:,i] = jacobi_polynomial(nodes[2,:], 0, 0, i-1);
            tmp3[:,i] = jacobi_polynomial(nodes[3,:], 0, 0, i-1);
        end
        for i=1:refel.N
            Dtmp1[:,i+1] = sqrt(i*(i+1)) .* jacobi_polynomial(nodes[1,:], 1, 1, i-1);
            Dtmp2[:,i+1] = sqrt(i*(i+1)) .* jacobi_polynomial(nodes[2,:], 1, 1, i-1);
            Dtmp3[:,i+1] = sqrt(i*(i+1)) .* jacobi_polynomial(nodes[3,:], 1, 1, i-1);
        end

        viXviXvi = kron(kron(refel.invV, refel.invV), refel.invV);

        VgXVgXVg = zeros(Nn, (refel.N+1)*(refel.N+1)*(refel.N+1));
        VgXVgXDg = zeros(Nn, (refel.N+1)*(refel.N+1)*(refel.N+1));
        VgXDgXVg = zeros(Nn, (refel.N+1)*(refel.N+1)*(refel.N+1));
        DgXVgXVg = zeros(Nn, (refel.N+1)*(refel.N+1)*(refel.N+1));
        for ni=1:Nn
            for i=1:(refel.N+1)
                for j=1:(refel.N+1)
                    for k=1:(refel.N+1)
                        ind = ((k-1)*(refel.N+1) + (j-1))*(refel.N+1) + i;
                        VgXVgXVg[ni,ind] = tmp1[ni,i]*tmp2[ni,j]*tmp2[ni,k];
                        VgXVgXDg[ni,ind] = Dtmp1[ni,i]*tmp2[ni,j]*tmp2[ni,k];
                        VgXDgXVg[ni,ind] = tmp1[ni,i]*Dtmp2[ni,j]*tmp2[ni,k];
                        DgXVgXVg[ni,ind] = tmp1[ni,i]*tmp2[ni,j]*Dtmp2[ni,k];
                    end
                end
            end
        end

        refel.Q = VgXVgXVg * viXviXvi;
        refel.Qr = VgXVgXDg * viXviXvi;
        refel.Qs = VgXDgXVg * viXviXvi;
        refel.Qt = DgXVgXVg * viXviXvi;
    end

    return refel;

end

# Map 1D nodes to 2D face
# Assumes interval is scaled the same
function map_face_nodes_2d(g1d, v1, v2)
    fnodes = zeros(2,length(g1d));
    dx = v2[1] - v1[1];
    dy = v2[2] - v1[2];
    dist = sqrt(dx*dx+dy*dy);
    xmult = dx/dist;
    ymult = dy/dist;

    fnodes[1,:] = v1[1] .+ (g1d .- v1[1]) .*xmult;
    fnodes[2,:] = v1[2] .+ (g1d .- v1[2]) .*ymult;

    return fnodes;
end

# Map 2D nodes to 3D face
# Assumes g2d[1] aligns with v1 and interval is scaled the same
function map_face_nodes_3d(g2d, v1, v2, v3, v4=[])
    fnodes = zeros(3,length(g2d));
    if length(v4) > 0 # quads
        dx = 0
        dy = 0
        dz = 0
        dist = sqrt(dx*dx+dy*dy+dz*dz);
        xmult = dx/dist;
        ymult = dy/dist;
        zmult = dz/dist;

        fnodes[1,:] = v1[1] .+ (g2d .- v1[1]) .*xmult;
        fnodes[2,:] = v1[2] .+ (g2d .- v1[2]) .*ymult;
        fnodes[3,:] = v1[3] .+ (g2d .- v1[3]) .*zmult;
    else #triangles

    end

    return fnodes;
end


text:
const FLAG_SIMULATION_NEGATIVE_BINOMIAL_VARY_RANK_SMALL                 = false
const FLAG_SIMULATION_NEGATIVE_BINOMIAL_VARY_RANK_MEDIUM                = true
const FLAG_SIMULATION_NEGATIVE_BINOMIAL_VARY_RANK_LARGE                 = false

const FLAG_SIMULATION_NEGATIVE_BINOMIAL_VARY_MISSING_PERCENTAGE_SMALL   = false
const FLAG_SIMULATION_NEGATIVE_BINOMIAL_VARY_MISSING_PERCENTAGE_MEDIUM  = false
const FLAG_SIMULATION_NEGATIVE_BINOMIAL_VARY_MISSING_PERCENTAGE_LARGE   = false

FLAG_SIMULATION_NEGATIVE_BINOMIAL_VARY_RANK_SMALL ?
  include("simulation_negbin_vary_rank_small.jl")  : nothing

FLAG_SIMULATION_NEGATIVE_BINOMIAL_VARY_RANK_MEDIUM ?
  include("simulation_negbin_vary_rank_medium.jl") : nothing

FLAG_SIMULATION_NEGATIVE_BINOMIAL_VARY_RANK_LARGE ?
  include("simulation_negbin_vary_rank_large.jl") : nothing


text:
module NSGAII

import YAML
using Cambrian
import JSON

include("evolution.jl")

end


text:
using SparseArrays, Graphs

A = sprand(Bool, 20, 20, 0.5)
for i in 1:20; A[i,i] = 0; end
A = A + A'
dropzeros!(A)
g = IndexedGraph(A)

@testset "undirected graph" begin
    @testset "basics" begin
        @test !is_directed(g)
        @test length(collect(edges(g))) == ne(g)
        i = 3
        es = inedges(g, i)
        neigs = neighbors(g, i)
        @test all(
            ((src(e), dst(e)) == (i, j)) || ((src(e), dst(e)) == (j, i))
            for (e, j) in zip(es, neigs)
        )
    end

    @testset "edge indexing" begin
        @test all( e == get_edge(g, src(e), dst(e)) for e in edges(g) )
        @test all( e == get_edge(g, idx(e)) for e in edges(g) )
    end

    @testset "construct from SimpleGraph" begin
        sg = SimpleGraph(A)
        ig = IndexedGraph(sg)
        @test adjacency_matrix(sg) == adjacency_matrix(ig)
    end
end


text:
export zoo_psdd, zoo_psdd_file,
    PsddFormat, PsddVtreeFormat

struct PsddFormat <: FileFormat end

const PsddVtreeFormat = Tuple{PsddFormat,VtreeFormat}
Tuple{PsddFormat,VtreeFormat}() = (PsddFormat(),VtreeFormat())

##############################################
# Read SDDs
##############################################

zoo_psdd_file(name) =
    artifact"circuit_model_zoo" * zoo_version * "/psdds/$name"

"""
    zoo_psdd(name)

Loads PSDD file with given name from model zoo. See https://github.com/UCLA-StarAI/Circuit-Model-Zoo.
"""
zoo_psdd(name) =
    read(zoo_psdd_file(name), ProbCircuit, PsddFormat())

const psdd_grammar = raw"""
    start: _header (_NL node)+ _NL?

    _header : "psdd" (_WS INT)?

    node : "T" _WS INT _WS INT _WS INT _WS LOGPROB -> true_node
         | "L" _WS INT _WS INT _WS SIGNED_INT _WS? -> literal_node
         | "D" _WS INT _WS INT _WS INT _WS elems  -> decision_node

    elems : elem (_WS elem)*
    elem : INT _WS INT _WS LOGPROB

    %import common.INT
    %import common.SIGNED_INT
    %import common.SIGNED_NUMBER -> LOGPROB
    %import common.WS_INLINE -> _WS
    %import common.NEWLINE -> _NL
    """ * dimacs_comments

const psdd_parser = Lark(psdd_grammar)

abstract type PsddParse <: JuiceTransformer end

@rule start(t::PsddParse, x) = begin
    x[end]
end

@rule elem(t::PsddParse, x) =
    [t.nodes[x[1]], t.nodes[x[2]], Base.parse(Float64,x[3])]

@rule elems(t::PsddParse, x) =
    Array(x)

#  parse unstructured
struct PlainPsddParse <: PsddParse
    nodes::Dict{String,PlainProbCircuit}
    PlainPsddParse() = new(Dict{String,PlainProbCircuit}())
end

@rule literal_node(t::PlainPsddParse, x) =
    t.nodes[x[1]] = PlainProbLiteralNode(Base.parse(Lit,x[3]))

@rule true_node(t::PlainPsddParse, x) = begin
    litn = PlainProbLiteralNode(Base.parse(Lit, x[3]))
    log_prob = Base.parse(Float64, x[4])
    log_probs = [log_prob, log1p(-exp(log_prob))]
    t.nodes[x[1]] = PlainSumNode([litn, -litn], log_probs)
end

@rule decision_node(t::PlainPsddParse,x) = begin
    @assert length(x[4]) == Base.parse(Int,x[3])
    children = map(x[4]) do elem
        PlainMulNode(elem[1:2])
    end
    log_probs = map(e -> e[3], x[4])
    t.nodes[x[1]] = PlainSumNode(children, log_probs)
end

function Base.parse(::Type{PlainProbCircuit}, str, ::PsddFormat)
    ast = Lerche.parse(psdd_parser, str)
    Lerche.transform(PlainPsddParse(), ast)
end

Base.read(io::IO, ::Type{PlainProbCircuit}, ::PsddFormat) =
    parse(PlainProbCircuit, read(io, String), PsddFormat())

#  parse structured
struct StructPsddParse <: PsddParse
    id2vtree::Dict{String,<:Vtree}
    nodes::Dict{String,StructProbCircuit}
    StructPsddParse(id2vtree) =
        new(id2vtree,Dict{String,StructProbCircuit}())
end

@rule literal_node(t::StructPsddParse, x) = begin
    lit = Base.parse(Lit,x[3])
    vtree = t.id2vtree[x[2]]
    t.nodes[x[1]] = StructProbLiteralNode(lit, vtree)
end

@rule true_node(t::StructPsddParse, x) = begin
    vtree = t.id2vtree[x[2]]
    litn = StructProbLiteralNode(Base.parse(Lit, x[3]), vtree)
    log_prob = Base.parse(Float64, x[4])
    log_probs = [log_prob, log1p(-exp(log_prob))]
    t.nodes[x[1]] = StructSumNode([litn, -litn], log_probs, vtree)
end

@rule decision_node(t::StructPsddParse,x) = begin
    @assert length(x[4]) == Base.parse(Int,x[3])
    vtree = t.id2vtree[x[2]]
    children = map(x[4]) do elem
        StructMulNode(elem[1], elem[2], vtree)
    end
    log_probs = map(e -> e[3], x[4])
    t.nodes[x[1]] = StructSumNode(children, log_probs, vtree)
end

function Base.parse(::Type{StructProbCircuit}, str::AbstractString, ::PsddFormat, id2vtree)
    ast = Lerche.parse(psdd_parser, str)
    Lerche.transform(StructPsddParse(id2vtree), ast)
end

function Base.parse(::Type{StructProbCircuit}, strings, format::PsddVtreeFormat)
    id2vtree = parse(Dict{String,Vtree}, strings[2], format[2])
    parse(StructProbCircuit, strings[1], format[1], id2vtree)
end

Base.read(io::IO, ::Type{StructProbCircuit}, ::PsddFormat, id2vtree) =
    parse(StructProbCircuit, read(io, String), PsddFormat(), id2vtree)

function Base.read(ios::Tuple{IO,IO}, ::Type{StructProbCircuit}, ::PsddVtreeFormat)
    circuit_str = read(ios[1], String)
    vtree_str = read(ios[2], String)
    parse(StructProbCircuit, (circuit_str,vtree_str), PsddVtreeFormat())
end

##############################################
# Write PSDDs
##############################################

const PSDD_FORMAT = """c this file was saved by ProbabilisticCircuits.jl
c ids of psdd nodes start at 0
c psdd nodes appear bottom-up, children before parents
c
c file syntax:
c psdd count-of-sdd-nodes
c L id-of-literal-sdd-node id-of-vtree literal
c T id-of-trueNode-sdd-node id-of-vtree variable log(litProb)
c D id-of-decomposition-sdd-node id-of-vtree number-of-elements {id-of-prime id-of-sub log(elementProb)}*
c"""

function Base.write(io::IO, pc::ProbCircuit, ::PsddFormat, vtree2id::Function = (x -> 0))

    id = -1

    println(io, PSDD_FORMAT)
    println(io, "psdd $(sdd_num_nodes_leafs(pc))")

    f_con(_) = error("ProbCircuits have no constant nodes")

    f_lit(n) = begin
        nid = id += 1
        println(io, "L $nid $(vtree2id(n)) $(literal(n))")
        nid
    end

    f_a(n, ids) = begin
        if length(ids) != 2
            error("The PSDD file format requires multiplications/AND nodes to have exactly two inputs")
        end
        tuple(ids...)
    end

    f_o(n, ids) = begin
        nid = id += 1
        vtreeid = vtree2id(n)
        if num_children(n) == 2 && all(isliteralgate, children(n))
            pos_child = literal(children(n)[1]) > 0 ? 1 : 2
            log_prob = n.log_probs[pos_child]
            v = variable(children(n)[1])
            print(io, "T $nid $vtreeid $v $log_prob")
        else
            print(io, "D $nid $vtreeid $(length(ids))")
            for (el, log_prob) in zip(ids, n.log_probs)
                print(io, " $(el[1]) $(el[2]) $log_prob")
            end
        end
        println(io)
        nid
    end

    foldup_aggregate(pc, f_con, f_lit, f_a, f_o, Union{Int, Tuple{Int,Int}})

    nothing
end

function Base.write(ios::Tuple{IO,IO}, pc::StructProbCircuit, format::PsddVtreeFormat)
    vtree2id = write(ios[2], vtree(pc), format[2])
    write(ios[1], pc, format[1], i -> vtree2id[vtree(i)])
end


text:
using Pkg
Pkg.test("App")

max_stars_repo_name: DevOpsDemoTF/DevOpsDemo-template-Julia

max_issues_repo_name: DevOpsDemoTF/DevOpsDemo-template-Julia

max_forks_repo_name: DevOpsDemoTF/DevOpsDemo-template-Julia


text:
using Unitful
using DataStructures

t = range(0.0, stop=10.0, length=100)

const MyResult = OrderedDict{String,Tuple{Bool,Any}}  # value = (isConstant, signal)
ModiaPlot.getRawSignal(result::MyResult, name) = result[name]

result = MyResult()

result["time"]     = (false, t*u"s")
result["phi_max"]  = (true , 1.1f0*u"rad")
result["i_max"]    = (true , 2)
result["open"]     = (true , true)
result["Inertia"]  = (true , [1.1  1.2  1.3;
                              2.1  2.2  2.3;
                              3.1  3.2  3.3]u"kg*m^2")

plot(result, ["phi_max", "i_max", "open", "Inertia[2,2]", "Inertia[1,2:3]", "Inertia"], heading="Constants")

println("\n... result info:")
printResultInfo(result)


text:
const makiegallery_dir = dirname(dirname(@__DIR__))

using Base.Threads

const current_ref_version = Ref{String}("v0.2.4")
function ref_image_dir(version = string(current_ref_version[]))
    return download_reference(version)
end
"""
Downloads the reference images from ReferenceImages for a specific version
"""
function download_reference(version = string(current_ref_version[]))
    refpath_version = version[1] == 'v' && version[2] in '0':'9' ? version[2:end] : version
    download_dir = joinpath(makiegallery_dir, "testimages")
    isdir(download_dir) || mkpath(download_dir)
    tarfile = joinpath(download_dir, "gallery.zip")
    url = "https://github.com/SimonDanisch/ReferenceImages/archive/$(version).tar.gz" # TODO FIXME
    refpath = joinpath(download_dir, "ReferenceImages-$(refpath_version)", "gallery")
    if !isdir(refpath) # if not yet downloaded
        download_images() = download(url, tarfile)
        try
            @info "downloading reference images for version $version"
            download_images()
        catch e
            if isa(e, ErrorException) && occursin("Hash Mismatch", e.msg)
                rm(tarfile, force = true)
                download_images()
            else
                rethrow(e)
            end
        end
        BinaryProvider.unpack(tarfile, download_dir)
        # check again after download
        if !isdir(refpath)
            error("Something went wrong while downloading reference images. Plots can't be compared to references")
        end
    end
    refpath
end

is_image_file(path) = lowercase(splitext(path)[2]) in (".png", ".jpg", ".jpeg")

function extract_frames(video, frame_folder)
    path = joinpath(frame_folder, "frames%04d.jpg")
    FFMPEG.ffmpeg_exe(`-loglevel quiet -i $video -y $path`)
end

function compare_media(a, b; sigma = [1,1], eps = 0.02)
    file, ext = splitext(a)
    if ext in (".png", ".jpg", ".jpeg", ".JPEG", ".JPG")
        imga = load(a)
        imgb = load(b)
        return approx_difference(imga, imgb, sigma, eps)
    elseif ext in (".mp4", ".gif")
        mktempdir() do folder
            afolder = joinpath(folder, "a")
            bfolder = joinpath(folder, "b")
            mkpath(afolder); mkpath(bfolder)
            extract_frames(a, afolder)
            extract_frames(b, bfolder)
            aframes = joinpath.(afolder, readdir(afolder))
            bframes = joinpath.(bfolder, readdir(bfolder))
            if length(aframes) > 10
                # we don't want to compare too many frames since it's time costly
                # so we just compare 10 random frames if more than 10
                samples = rand(1:length(aframes), 10)
                aframes = aframes[samples]
                bframes = bframes[samples]
            end
            # test by maximum diff
            return mean(compare_media.(aframes, bframes; sigma = sigma, eps = eps))
        end
    else
        error("Unknown media extension: $ext")
    end
end

"""
Compares all media recursively in two recorded folders!
"""
function run_comparison(
        test_record_path, test_diff_path,
        reference = MakieGallery.download_reference();
        maxdiff = 0.032
    )
    @testset "Reference Image Tests" begin
        folders = joinpath.(test_record_path, readdir(test_record_path))
        l = length(folders)
        count = 1
        for folder in folders
            @debug "Running index $count" progress=count/l
            count += 1
            if isdir(folder) && isdir(joinpath(reference, basename(folder)))
                media = joinpath(folder, "media")
                ref_folder = joinpath(reference, basename(folder), "media")
                test_folder = joinpath(test_record_path, media)
                ref_media = sort(readdir(ref_folder))
                test_media = sort(readdir(test_folder))
                @testset "$(basename(folder))" begin
                    if isempty(test_media)
                        @warn("recodings are missing for $folder")
                    else
                        for test in test_media
                            ref = joinpath(ref_folder, test)
                            test = joinpath(test_folder, test)
                            diff = compare_media(ref, test)
                            if diff >= maxdiff
                                refdiff = joinpath(test_diff_path, basename(folder), string("ref_", basename(ref)))
                                testdiff = joinpath(test_diff_path, basename(folder), string("test_", basename(test)))
                                mkpath(dirname(refdiff)); mkpath(dirname(testdiff))
                                cp(ref, refdiff, force = true)
                                cp(test, testdiff, force = true)
                                @test diff <= maxdiff
                            else
                                @test diff <= maxdiff
                            end
                        end
                    end
                end
            end
        end
    end
end


text:
using Base.Test
using OptimTools

probs = Any[];

# Rosenbrock's function
	f1(x)   = (1-x[1]).^2 + 100*(x[2]- x[1].^2).^2
	df1(x)  = [
				400*x[1].^3 - 400*x[1].*x[2] + 2*x[1]-2;
				200*x[2]-200*x[1].^2
			]
	d2f1(x) = [1200*x[1].^2-400*x[2]+2 -400*x[1]; -400*x[1] 200]
xsol1   = [1,1]
push!(probs,("Rosenbrock",f1,df1,d2f1,xsol1))

# convex non-quadratic
f2(x)  = exp(x[1]+3*x[2]-0.1)+exp(x[1]-3*x[2]-0.1)+exp(-x[1]-0.1)
df2(x) = [
            exp(x[1]+3*x[2]-0.1)+exp(x[1]-3*x[2]-0.1)-exp(-x[1]-0.1);
            3*exp(x[1]+3*x[2]-0.1)-3*exp(x[1]-3*x[2]-0.1);
        ]
d2f2(x)  = [
            exp(x[1]+3*x[2]-0.1)+exp(x[1]-3*x[2]-0.1)+exp(-x[1]-0.1)   3*exp(x[1]+3*x[2]-0.1)-3*exp(x[1]-3*x[2]-0.1);
            3*exp(x[1]+3*x[2]-0.1)-3*exp(x[1]-3*x[2]-0.1)              9*exp(x[1]+3*x[2]-0.1)+9*exp(x[1]-3*x[2]-0.1);
        ]
xsol2   = [-0.34657;0.00000]
push!(probs,("ConvexNonQuad",f2,df2,d2f2,xsol2))

# convex quadratic
A      = randn(3,3)
A      = A'*A + eye(3)*3
b      = randn(3);
f3(x)   = 0.5*dot(x,A*x) - dot(b,x)
df3(x)  = A*x - b
d2f3(x) = A
xsol3   = A\b
push!(probs,("ConvexQuad",f3,df3,d2f3,xsol3))

# test test functions
for k=1:length(probs)
	prob = probs[k]

	println("testing function $(prob[1])")
	# test xsol
	@test norm(prob[3](prob[5])) < 1e-5

	# setup test function for derivative check
	function testFun(x,v=[])
		if isempty(v)
			return prob[2](x)
		else
			return prob[2](x), prob[3](x)'*v
		end
	end

	pass, = checkDerivative(testFun,prob[5]+randn(size(prob[5])),out=false)
	@test pass

	# setup test function for derivative check
	function testFun(x,v=[])
		if isempty(v)
			return prob[3](x)
		else
			return prob[3](x), prob[4](x)'*v
		end
	end
	pass, = checkDerivative(testFun,prob[5]+randn(size(prob[5])),out=false)
	@test pass
end


text:
using Documenter, ACE

makedocs(sitename="ACE.jl Documentation",
         pages = [
        "Home" => "index.md",
        "Getting Started" => "gettingstarted.md",
        "Math" => "math.md",
        "Developer Docs" => "devel.md",
        "AD" => "ad.md",
        "Types & Functions" => "docs.md"
         ])

        #  https://journals.aps.org/prb/abstract/10.1103/PhysRevB.99.014104

deploydocs(
    repo = "github.com/ACEsuit/ACE.jl.git",
    devbranch = "main"
)


text:
using LinearAlgebra
using Plots
using Distributions
using Random

function VanillaSGD(gradcalc::Function, nn_predict::Function, maxIter, nHidden,nParams, xtrain, ytrain, xvalid, yvalid,a=.0001, B=1)
    print("Vanilla SGD Running...")
    n = size(xtrain,1)
    W= randn(nParams,1)
    flist=[]
    valid=[]
    wbest=zeros(nParams)
    vallow=Inf
    for t in 1:maxIter

    	i = Base.rand(1:n,B)
    	f,g = gradcalc(W, xtrain[i,:], ytrain[i], nHidden)
    	W = W - a*g

        if (mod(t-1,round(maxIter/100)) == 0)
            push!(flist,f)
            ComputeValid(nn_predict, W, xvalid, yvalid, valid)
            if valid[end]<vallow
                wbest=W
                vallow=valid[end]
            end
        end
    end
    display(plot(1:length(valid), valid, xlabel="Number of Iterations", ylabel="Validation Error", legend=false))
    return wbest, vallow
end

function SGDBabysitter(gradcalc::Function, nn_predict::Function, maxIter, nHidden,nParams, xtrain, ytrain, xvalid, yvalid, maxbatch=500)

    n = size(xtrain,1)
    valid=[]
    flist=[]
    alist=[]
    Blist=[]

    a,B,W=BS_initialize(gradcalc, nn_predict, maxIter, nHidden,
                            nParams, xtrain, ytrain, xvalid, yvalid)
    push!(Blist,B)
    push!(alist,a)
    wbest=zeros(nParams)
    vallow=Inf
    triedDec=false
    for t in 1:maxIter

        idxshuffle = shuffle(1:n)
        i = rand(idxshuffle, B)
    	# i = rand(1:n,B)

    	f,g = gradcalc(W, xtrain[i,:], ytrain[i], nHidden)
    	W = W - a*g

        #Every few iterations, compute validation and append f,g,validation
        if (mod(t-1,round(maxIter/100)) == 0)
              push!(flist,f)
              ComputeValid(nn_predict, W, xvalid, yvalid, valid)
              if valid[end]<vallow
                  wbest=W
                  vallow = valid[end]
              end
        end
        if (mod(t-1,round(maxIter/50)) == 0)
            a,B,W,triedDec=BS_Select(validation_array=valid,f_array=flist, gradcalc=gradcalc,xtrain=xtrain,ytrain=ytrain,W=W,nHidden=nHidden, alist=alist, Blist=Blist,maxbatch=maxbatch, wbest=wbest,  maxIter=maxIter,triedDec=triedDec)

          end
    end

    display(plot(1:length(valid), valid, xlabel="Number of Iterations", ylabel="Validation Error", legend=false))

    return wbest, vallow
end

#INITIALIZATION FUNCTION--------------------------------------------------------
function BS_initialize(gradcalc::Function, nn_predict::Function, maxIter, nHidden,nParams, xtrain, ytrain, xvalid, yvalid)
switch=1
n = size(xtrain,1)
B=5
a=rand(Uniform(.095,.09999999))
j=1
W_init = randn(nParams,1)
W = zeros(nParams,1)
while switch==1 && j<10
    W=W_init
    flist=[]
    validlist=[]
    for t in 1:10
        i = rand(1:n,B)
        f,g = gradcalc(W, xtrain[i,:], ytrain[i], nHidden)
        W = W - a*g
        push!(flist,f)
        ComputeValid(nn_predict, W, xvalid, yvalid, validlist)
    end
    if validlist[end]<(4/5)*validlist[1] && flist[end]<(4/5)*flist[1]
        switch=0
    else
        a=rand(Uniform(.5,.75))*a
        print("Initializing... ")
        j+=1
    end
end
return a,B,W
end

#Selection Function-------------------------------------------------------------

ctr = 0
maxIter = 10000
decLimit = round(maxIter/(100*20))
threshold = round(maxIter/(100*10))
function BS_Select(;validation_array,f_array, gradcalc::Function,xtrain,ytrain,W,nHidden, alist, Blist, maxbatch, wbest,  maxIter,triedDec)
    global ctr, decLimit, threshold
    #Get most recent step and batch sizes
    a=alist[end]
    B=Blist[end]
    n,d=size(xtrain)

    #Get Gradient Angle
    i=rand(1:n,B)
    j=rand(1:n,B)
    _,gi=gradcalc(W, xtrain[i,:],ytrain[i],nHidden)
    _,gj=gradcalc(W,xtrain[j,:],ytrain[j], nHidden)
    costheta=dot(gi,gj)/(norm(gi)*norm(gj))

    tol = 1e-5
    if costheta > 1
        costheta = 1 - tol
    elseif costheta < -1
        costheta = -1 + tol
    end
    theta = acos(costheta)

    #Case 1: alpha decrease if validation error does not

    if length(validation_array) >= 4 && ((validation_array[end-3] < validation_array[end]||validation_array[end-2] < validation_array[end]) && (validation_array[end-1]<validation_array[end])) || length(validation_array) >= 4 && ((mean(validation_array[end-3:end-1]) <= 1.05*validation_array[end]) && (mean(validation_array[end-3:end-1]) >= 0.95*validation_array[end]))
        if triedDec==false
            a=(3/4)*a
            triedDec = true
        elseif triedDec==true
            #go back to best w
            W = wbest
            a = (4/5)*a
            # a = (4/3)*a
            triedDec = false
        end

    #---- Alterative to case 1 ----#
    #Don't delete, might use later
    #Tries increasing alpha

    # if length(validation_array) >= 4 && (((validation_array[end-3] < validation_array[end]||validation_array[end-2] < validation_array[end]) && (validation_array[end-1]<validation_array[end])) || ((mean(validation_array[end-3:end-1]) <= 1.05*validation_array[end]) && (mean(validation_array[end-3:end-1]) >= 0.95*validation_array[end])))
    #     if ctr <= decLimit
    #         a = (3/4)*a
    #         print("case 2 part 1")
    #         print(" Ctr is :",ctr, " ")
    #         ctr += 1
    #     elseif ctr > decLimit && ctr <= threshold
    #         a = (10)*a
    #         # ctr += 1
    #         ctr = 0
    #         decLimit += 1
    #         print("case 2 part 2")
    #     elseif ctr > threshold
    #         W = wbest
    #         ctr = 0
    #         print("case 2 part 3")
    #     end

    #Case 2: alpha decreases based on validation error
    # elseif length(validation_array) >= 4 && ((validation_array[end-3] < validation_array[end]||validation_array[end-2] < validation_array[end]) && (validation_array[end-1]<validation_array[end]))
    elseif length(validation_array) >= 4 && validation_array[end-2]>validation_array[end]
        a_adjust = (1 - abs(validation_array[end-2] - validation_array[end])/validation_array[end-2])
        a = a*a_adjust
        print("case 1")
    end

    #case 3 batchsize as a function of theta
    if exp(-.5/(pi-theta))==NaN
        B=maxbatch
    else
        B=convert(Int64,round(maxbatch*(1/(1+round(exp(-.5/(pi-theta+.00001)),digits=5)))))
    end

    # if B > maxbatch
    #     B = maxbatch
    # elseif B < 1
    #     B = 1
    # else
    #     if round(maxbatch*cos(theta/2)) == NaN
    #         B = maxbatch
    #     else
    #         B = convert(Int64, round(maxbatch*cos(theta/2)))
    #     end
    # end

    if a > 1
        a = 1
    end
    push!(alist,a)
    push!(Blist,B)
    print("Step Size: ", a," ","Batch Size: ", B, " ")
    return a,B,W,triedDec
end

#Compute Validation-------------------------------------------------------------
function ComputeValid(predict_func::Function, W, xvalid, yvalid, valid)
    yhat = predict_func(W,[ones(size(xvalid,1)) xvalid],nHidden)
    push!(valid,sum((yhat-yvalid).^2)/length(yvalid))
 end


text:
# Practical Machine learning
# Decision Trees based learning - Random Forests example
# Chapter 5

module RandomForests

export
    RandomForestClassifier,
    RandomForestRegressor,
    fit,
    predict,
    feature_importances,
    oob_error

using DataFrames

include("util.jl")
include("tree.jl")
include("randomforest.jl")
include("classifier.jl")
include("regressor.jl")

end # RandomForests module