text:
using Tokenize
function recipe_files(dir)
    ret = String[]
    for (root, dirs, files) in walkdir(dir)
        startswith(root, '.') && continue
        for file in files
            splitext(file)[2] == ".jl" || continue
            path = joinpath(root, file)
            if contains(readstring(path), "@recipe")
                push!(ret, path)
            end
        end
    end
    ret
end

immutable RecipeLocation
    startline::Int64
    stopline::Int64
    function RecipeLocation(startline, stopline)
        @assert startline <= stopline
        new(startline, stopline)
    end
end

type State
    searching_recipe_end::Bool
    current_recipe_start_line::Int
    current_line::Int
end

const OUTSIDE_LINE = -1
searching_recipe_end(s::State) = s.searching_recipe_end
searching_recipe_end!(s, b) = (s.searching_recipe_end = b; b)
function finish_recipe!(s::State)
    @assert searching_recipe_end(s)
    ret = RecipeLocation(s.current_recipe_start_line, s.current_line - 1)
    s.current_recipe_start_line = OUTSIDE_LINE
    searching_recipe_end!(s, false)
    ret
end
function start_recipe!(s::State)
    @assert searching_recipe_end(s) == false
    searching_recipe_end!(s, true)
    s.current_recipe_start_line = s.current_line
    s
end
State() = State(false, OUTSIDE_LINE, 1)

locate_recipes!(ret, ex, state) = ret
function locate_recipes!(ret, ex::Expr, state)
    if isline(ex)
        state.current_line = ex.args[1]
    end
    if searching_recipe_end(state) && isline(ex)
        push!(ret, finish_recipe!(state))
    end
    if isrecipe(ex)
        start_recipe!(state)
    end
    if !isrecipe(ex)
        for arg in ex.args
            locate_recipes!(ret, arg, state)
        end
    end
    ret
end

function locate_recipes(s::String)
    ret = RecipeLocation[]
    state = State()
    ex = parse("begin $s end")
    locate_recipes!(ret, ex, state)
    if searching_recipe_end(state)
        state.current_line = length(split(s, '\n')) + 1
        push!(ret, finish_recipe!(state))
    end
    @assert !searching_recipe_end(state)
    ret
end

function recipe_lines(s)
    lines = Int64[]
    for recipe in locate_recipes(s)
        append!(lines, recipe.startline:recipe.stopline)
    end
    lines
end

function replace_token_in_recipes(s, pattern, r)
    index = recipe_lines(s)
    lines = String[]
    for (i, line) in enumerate(String.(split(s, '\n')))
        if i in index
            line_new = ""
            for token in tokenize(line)
                if string(token) == pattern
                    line_new *= r
                else
                    line_new *= string(token)
                end
            end
            if line_new != line
                println(i)
                println(line)
                println(line_new)
            end
            push!(lines, line_new)
        else
            push!(lines, line)
        end
    end
    join(lines, '\n')
end

function replace_token_in_recipes_dir(dir, pattern, r)
    for path in recipe_files(dir)
        s = readstring(path)
        @show path
        s_new = replace_token_in_recipes(s, pattern, r)
        write(path, s_new)
    end
end

function replace_d_pkg(pkg)
    # Pkg.checkout(pkg)
    dir = Pkg.dir(pkg)
    run(Cmd(`git branch replaceD`, dir=dir))
    run(Cmd(`git checkout replaceD`, dir=dir))
    pattern = "d"
    r = "plotattributes"
    replace_token_in_recipes_dir(dir, pattern, r)
    msg = "replace d by plotattributes inside recipes"
    run(Cmd(`git commit -am $msg`, dir=dir))
end

function search_token_in_recipes(f, path)
    index = recipe_lines(readstring(path))
    for (i,line) in enumerate(readlines(path))
        if i in index
            for token in tokenize(line)
                if f(token)
                    println(" "^8, i, " ", line)
                    break
                end
            end
        end
    end
end


text:
include("naive_bayes.jl")
include("bernoulli_nb.jl")

d = GaussianNaiveBayes([0.0 10.0;
	                    0.0 2.0],
	                   [1.0 1.0;
	                    1.0 1.0],
	                   [0.25, 0.75])

X, c = rand(d, 10_000)

d = fit(GaussianNaiveBayes, X, c)

predict(d, X)

loglikelihood(d, X, c)

logpdf(d, X[:, 1], c[1])
logpdf(d, X, c)

mean(d)

fit(GaussianNaiveBayes, rand(d, 1_000)...)

cor(X[:, c .== 1]')
cor(X[:, c .== 2]')

#

d = BernoulliNaiveBayes([0.1 0.9;
	                     0.9 0.1],
	                    [0.25, 0.75])

X, c = rand(d, 10_000)

d = fit(BernoulliNaiveBayes, X, c)

predict(d, X)

loglikelihood(d, X, c)

logpdf(d, X[:, 1], c[1])
logpdf(d, X, c)

mean(d)

fit(BernoulliNaiveBayes, rand(d, 1_000)...)

cor(X[:, c .== 1]')
cor(X[:, c .== 2]')

ext: jl


text:
using Hearthstone
using Base.Test

@test_approx_eq_eps Hearthstone.card_in_opening_hand_coin(2) 0.5098 1e-4
@test_approx_eq_eps Hearthstone.card_in_opening_hand_no_coin(2) 0.4125 1e-4
@test_approx_eq Hearthstone.card_in_next_N_turns(1, 20, 1) 1/20
@test Hearthstone.missile_lethal(10, 8, 2) == 1.0


text:
using Test
using DataFrames
using StatPlots; gr()
using ForestBiometrics
#no idea if i should use this here or in individual scripts?

@testset "EquilibriumMoistureContent" begin include("EquilibriumMoistureContent_test.jl") end

@testset "ForestStocking" begin include("ForestStocking_test.jl") end

@testset "GingrichStocking_chart" begin include("GingrichStocking_chart_test.jl") end

@testset "HeightDub" begin include("HeightDub_test.jl") end

@testset "LimitingDistance" begin include("LimitingDistance_test.jl") end

@testset "SDI_chart" begin include("SDI_chart_test.jl") end

@testset "VolumeEquations" begin include("VolumeEquations_test.jl") end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    EstimatorFitter

A method to fit a density ratio estimator.
"""
abstract type EstimatorFitter end

"""
    fit(DRE, x_nu, x_de, ranges; [optlib])

Perform hyperparameter tuning of density ratio
estimator `dre` with numerator and denominator
samples, `x_nu` and `x_de` and with hyperparameter
`ranges`. Optionally, specify the optimization
library `optlib`.

### Notes

Hyperparameter tuning is not defined for all
density ratio estimators. Therefore, this
function may not work with some estimators.
"""
fit(dre::Type{<:DensityRatioEstimator}, x_nu, x_de,
    fitter::EstimatorFitter; optlib=default_optlib(dre)) =
  _fit(dre, x_nu, x_de, fitter, optlib)

################################################
## functions to be implemented by new fitters ##
################################################

_fit(dre::Type{DensityRatioEstimator}, x_nu, x_de,
     fitter::EstimatorFitter, optlib::Type{OptimizationLibrary}) =
  @error "not implemented"


text:
export MatchState
"""
State object for a match that is defined by the triple `parser, sequence, position`.

!!! note:
    Performance tip: [`Atomic`](@ref) is masking the state of its wrapped parser with `MatchState`.
    This simplifies the state
"""
struct MatchState end
Base.show(io::IO, ::MatchState) = print(io,"∘")

"""
State object representing ncodeunits explicitely with state of match for `prevind`, `nextind` to improve performance.
    `nc::Int` and `state::S`.

See also [`MatchState`](@ref), [`prevind`](@ref), [`nextind`](@ref).

!!! note:
    `nc` as type parameter faster but slow compilation.
"""
struct NCodeunitsState{S}
    nc::Int
    state::S
end
@inline _leftof(str,i,parser::CombinedParser,x::NCodeunitsState) = i-x.nc
@inline _rightof(str,i,parser::CombinedParser,x::NCodeunitsState) = i+x.nc
"""
    NCodeunitsState(posi::Int,after::Int,state)
    NCodeunitsState{S}(posi::Int,after::Int,state)

returns `(tuple_pos, NCodeunitsState(after-posi,state))`.

!!! note:
    rename
"""
NCodeunitsState(posi::Int,after::Int,state) =
    after, NCodeunitsState(after-posi,state)
@inline NCodeunitsState{S}(posi::Int,after::Int,state) where S =
    after, NCodeunitsState{S}(after-posi,state)

"""
    tuple_pos(pos_state::Tuple)

[`_iterate`](@ref) returns a tuple `pos_state` or nothing, and
`pos_state[1]` is position after match.
"""
@inline tuple_pos(pos_state::Tuple, default...) = pos_state[1]
@inline tuple_pos(pos_state::Nothing, default) = default

"""
    tuple_state(pos_state::Tuple)

[`_iterate`](@ref) returns a tuple `pos_state` or nothing, and
`pos_state[2]` is the state of match.
"""
@inline tuple_state(pos_state::Tuple) = pos_state[2]
@inline tuple_state(pos_state::Nothing) = nothing


text:
using TensorBoardLogger, Logging, Plots
using Test

@testset "Optional: Plots.jl" begin
logger = TBLogger(test_log_dir*"Plots-jl", tb_overwrite)
p = plot(rand(5))
p2 = plot(rand(10))
with_logger(logger) do
    @info "dispatch" p
    @info "dispatch" twoplots=[p, p2]
end
log_image(logger, "explicit/p", p)

# test serialization
pb = PipeBuffer()
show(pb, MIME("image/png"), p)
data = TensorBoardLogger.preprocess("key", p, Vector())
@test length(data) == 1
@test first(data[1]) == "key"
@test last(data[1]) isa TensorBoardLogger.PngImage
@test last(data[1]).data == pb.data

# test unpacking of array of plots
data = TensorBoardLogger.preprocess("key", [p, p], Vector())
@test length(data) == 2
@test first(data[1]) == "key/1"
@test first(data[2]) == "key/2"
@test last(data[1]) isa TensorBoardLogger.PngImage
@test last(data[2]) isa TensorBoardLogger.PngImage

close.(values(logger.all_files))
end


text:
using Lindenmayer, Luxor, Colors, ColorSchemes

crystal = LSystem(Dict(
   "F" => "9F[F-]+*",
   ),
   "F")

plant = LSystem(Dict(
  "A" => "UBB8D", # initialize
  "X" => "*[-F*X*]+F*X"),
  "AX")

global x = 0
function f(t::Turtle)
      pos = Point(t.xpos, t.ypos)
      if x == 0
            # we'll just do this at the very start
            sethue("black")
            circle(O, 245, :clip)
            paint()
      end

      d = distance(pos, boxbottomcenter(BoundingBox()))
      setcolor([Luxor.julia_purple, Luxor.julia_red,  Luxor.julia_green][rand(1:end)])
      circle(pos, 10, :fill)
      global x += 1
end

drawLSystem(plant,
      forward = 70,
      turn = 17,
      iterations=4,
      #startingx=-250,
      #startingy=250,
      startingorientation = -π/2,
      startingpen = (1, 1, 1),
      width=500,
      height=500,
      asteriskfunction = f,
      filename="docs/src/assets/logo.png",
      backgroundcolor = RGBA(1, 1, 1, 0))

max_stars_repo_name: cormullion/Lindenmayer

max_issues_repo_name: cormullion/Lindenmayer

max_forks_repo_name: cormullion/Lindenmayer


text:
#___  ____       ___ _____   _   _       _ _
#|  \/  | |     |_  /  __ \ | | | |     (_| |
#| .  . | |       | | /  \/ | | | |_ __  _| |_ ___
#| |\/| | |       | | |     | | | | '_ \| | __/ _ \
#| |  | | |___/\__/ | \__/\ | |_| | | | | | || (_) |
#_______\_____\____/ \____/  _____|_|___|_____\_____ _____ _____
#| ___ \        (_)         | | \ \ / / / __  |  _  / __  |  _  |
#| |_/ _ __ ___  _  ___  ___| |_ \ V /  `' / /| |/' `' / /| |/' |
#|  __| '__/ _ \| |/ _ \/ __| __|/   \    / / |  /| | / / |  /| |
#| |  | | | (_) | |  __| (__| |_/ /^\ \ ./ /__\ |_/ ./ /__\ |_/ /
#\_|  |_|  \___/| |\___|\___|\__\/   \/ \_____/\___/\_____/\___/
#              _/ |
#             |__/
#
# This code is part of the proposal of the team "MLJC UniTo" - University of Turin
# for "ProjectX 2020" Climate Change for AI.
# The code is licensed under MIT 3.0
# Please read readme or comments for credits and further information.

# Compiler: Julia 1.5

# Short description of this file: Numerical example, model selection using
#ModelingToolkit
using ModelingToolkit, OrdinaryDiffEq

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

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

sys = ODESystem(eqs)
sys = ode_order_lowering(sys)

u0 = [D(x) => 2.0,
      x => 1.0,
      y => 0.0,
      z => 0.0]

p  = [σ => 28.0,
      ρ => 10.0,
      β => 8/3]

tspan = (0.0,100.0)
prob = ODEProblem(sys,u0,tspan,p,jac=true)
sol = solve(prob,RK4())
using Plots; plot(sol,vars=(x,y))


text:
using TextAnalysis
using Statistics:mean
using Transformers
using Transformers.BidirectionalEncoder:WordPiece
using PreSumm.Data
using PreSumm.Model

Sample = NamedTuple{(:reference, :candidate),Tuple{String,String}}

function evaluate(
    samples::Channel{SummaryPair},
    model::Translator,
    wordpiece::WordPiece,
    tokenizer::Function,
    evaluation::Function
)
    vocabulary = Vocabulary(wordpiece)
    translated = Channel{Sample}() do channel
        for sample ∈ samples
            candidate_tokens = model(
                ["[CLS]"; wordpiece(tokenizer(sample.source)); "[SEP]"],
                vocabulary
            )
            candidate = join(candidate_tokens[2:end - 1], " ")
            reference = sample.target
            push!(channel, (reference = reference, candidate = candidate))
        end
    end
    return evaluation(translated)
end

function mean_rouge_n(samples::Channel{Sample}, n::Integer)
    return mean(
        rouge_n([sample.reference], sample.candidate, n)
        for sample ∈ samples
    )
end

function mean_rouge_n(
    samples::Channel{SummaryPair},
    model::Translator,
    wordpiece::WordPiece,
    tokenizer::Function,
    n::Integer
)
    return evaluate(
        samples,
        model,
        wordpiece,
        tokenizer,
        samples -> mean_rouge_n(samples, n)
    )
end

function mean_rouge_l(samples::Channel{Sample}, β::AbstractFloat)
    return mean(
        rouge_l_summary([sample.reference], sample.candidate, β)
        for sample ∈ samples
    )
end

function mean_rouge_l(
    samples::Channel{SummaryPair},
    model::Translator,
    wordpiece::WordPiece,
    tokenizer::Function,
    β::AbstractFloat
)
    return evaluate(
        samples,
        model,
        wordpiece,
        tokenizer,
        samples -> mean_rouge_l(samples, β)
    )
end


text:
import Base: isempty

export Tank, drymass, totalmass, propellant, propellant!
export refill!, consume!
export mass

"Stores propellant for a rocket"
abstract type Tank end

"""
    totalmass(t::Tank) -> Float64
Mass of propellant tank `t` when it is full.
"""
function totalmass end

"""
    drymass(t::Tank) -> Float64
Mass of propellant tank `t` when it is empty.
"""
function drymass end

"""
    propellant(t::Tank)
Get remaining propellant in tank. Propellant is fuel plus oxidizer
"""
propellant(tank::Tank) = tank.propellant

"""
    propellant!(t::Tank, amount)
Set amount of propellant in tank to `amount`.
"""
function propellant!(tank::Tank, amount::Real)
    if 0 <= amount + drymass(tank) <= totalmass(tank)
        tank.propellant = amount
    else
        msg = "Propellant mass plus dry mass must be less than total mass"
        throw(DomainError(amount, msg))
    end
end

"""
    isempty(t::Tank)
Check if tank `t` is empty (has no propellant left).
"""
isempty(tank::Tank) = tank.propellant <= 0

"""
    mass(t::Tank)
The total mass of tank `t`, dry mass plus propellant.
"""
mass(tank::Tank) = drymass(tank) + propellant(tank)

"""
    refill!(tank::Tank) -> Tank
Fill propellant tankt to the max. Returns full tank
"""
function refill!(tank::Tank)
    propellant!(tank, totalmass(tank) - drymass(tank))
    tank
end

"""
    consume!(t::Tank, a::Real)
Comsume `a` amount of propellant in tank `t`.
"""
function consume!(tank::Tank, amount::Real)
    remaining = max(propellant(tank) - amount, 0)
    propellant!(tank, remaining)
    remaining
end


text:
"""
    ind,rse,IND = PartitionTreeMatrixDhillon_QuadTree(matrix,jmax)

Perform the quadtree partition of matrix.

### Input Arguments
* `matrix`: The matrix to be partitioned
* `jmax`: Specify the max level of partitioning the matrix

### Output Arguments
* `ind`: Vector with the same length as length(matrix[:]), indicate the order of matrix after partition.
* `rse`: Region start and end. Each pair represents the start and end of a region in the matrix, corresponding to one node in the quadtree
* `IND`: IND[i] indicate the order of matrix after iteration i.

Copyright 2018 The Regents of the University of California

Implemented by Yiqun Shao (Adviser: Dr. Naoki Saito)
"""
function PartitionTreeMatrixDhillon_QuadTree(matrix,jmax)

# 0.Preliminary

# Constants
rows,cols = size(matrix);
N = rows*cols
ind = Matrix{Int64}(reshape(1:N,rows,cols))
IND = zeros(rows,cols,jmax);
IND[:,:,1] = ind;

# region starts & ends
rse = Vector{Vector{Int64}}(jmax)
rse[1] = [1,N];

# 1. Partition the graph to yield rse, ind and jmax

j=1;
regioncount = 0;

while j<=jmax-1
    # add a column for j+1, if necessary

    rr = 1#for tracking the child regions

    rse[j+1] = zeros(Int64,4*length(rse[j]))

    for r = 1:2:(length(rse[j]) - 1)
        rs_rows,rs_cols = ind2sub((rows,cols),rse[j][r])
        re_rows,re_cols = ind2sub((rows,cols),rse[j][r+1])

        n_rows = re_rows-rs_rows+1
        n_cols = re_cols-rs_cols+1

        indrs = ind[rs_rows:re_rows,rs_cols:re_cols]

        if n_rows > 1 && n_cols > 1

            # compute the left and right singular vectors
            u_rows,v_cols = second_largest_singular_vectors(reshape(matrix[indrs[:]],n_rows,n_cols))

            # partition the current region
            pm,_ = partition_fiedler(SparseMatrixCSC{Float64,Int}(ones(1,1)),v = vcat(u_rows, v_cols))
            pm_rows = pm[1:n_rows]
            pm_cols = pm[n_rows+1:end]

            # determine the number of points in child region 1
            n1_rows = sum(pm_rows .> 0);
            n1_cols = sum(pm_cols .> 0);

            # update the indexing
            ind[rs_rows:re_rows,rs_cols:re_cols] = indrs[vcat(find(pm_rows .> 0),find(pm_rows.<=0)),vcat(find(pm_cols.>0),find(pm_cols.<=0))]

            if !(n1_rows==0 || n1_rows == n_rows) && !(n1_cols==0 ||n1_cols==n_cols)
                child_rows = [rs_rows, rs_rows+n1_rows-1, rs_rows+n1_rows, re_rows, rs_rows, rs_rows+n1_rows-1, rs_rows+n1_rows, re_rows];
                child_cols = [rs_cols, rs_cols+n1_cols-1, rs_cols, rs_cols+n1_cols-1, rs_cols+n1_cols, re_cols, rs_cols+n1_cols, re_cols];
                for k = 1:8
                    rse[j+1][rr+k-1] = sub2ind((rows,cols),child_rows[k],child_cols[k])
                end
                rr = rr + 8

            elseif (n1_rows == 0 || n1_rows == n_rows) && !(n1_cols == 0 || n1_cols == n_cols)
                child_rows = [rs_rows, re_rows, rs_rows, re_rows];
                child_cols = [rs_cols, rs_cols+n1_cols-1, rs_cols+n1_cols, re_cols];
                for k = 1:4
                    rse[j+1][rr+k-1] = sub2ind((rows,cols),child_rows[k],child_cols[k])
                end
                rr = rr + 4

            elseif (n1_cols == 0 || n1_cols == n_cols) && !(n1_rows == 0 || n1_rows == n_rows)
                child_rows = [rs_rows, rs_rows+n1_rows-1, rs_rows+n1_rows, re_rows];
                child_cols = [rs_cols, re_cols, rs_cols, re_cols];
                for k = 1:4
                    rse[j+1][rr+k-1] = sub2ind((rows,cols),child_rows[k],child_cols[k])
                end
                rr = rr + 4
            else
                rse[j+1][rr:rr+1] = rse[j][r:r+1]
                rr = rr + 2
            end

        elseif n_rows == 1 && n_cols > 1 #no bisection on rows
            pm_cols,_ = partition_fiedler(SparseMatrixCSC{Float64,Int}(ones(1,1)),v = matrix[indrs[:]])
            n1_cols = sum(pm_cols .> 0);
            ind[rs_rows:re_rows,rs_cols:re_cols] = indrs[:,vcat(find(pm_cols.>0),find(pm_cols.<=0))]

            child_rows = [rs_rows, re_rows, rs_rows, re_rows];
            child_cols = [rs_cols, rs_cols+n1_cols-1, rs_cols+n1_cols, re_cols];
            for k = 1:4
                rse[j+1][rr+k-1] = sub2ind((rows,cols),child_rows[k],child_cols[k])
            end
            rr = rr + 4

        elseif n_rows > 1 && n_cols == 1 #no bisection on cols
            pm_rows,_ = partitio_fiedler(SparseMatrixCSC{Float64,Int}(ones(1,1)),v = matrix[indrs[:]])
            n1_rows = sum(pm_rows .> 0);
            ind[rs_rows:re_rows,rs_cols:re_cols] = indrs[vcat(find(pm_rows .> 0),find(pm_rows.<=0)),:]

            child_rows = [rs_rows, rs_rows+n1_rows-1, rs_rows+n1_rows, re_rows];
            child_cols = [rs_cols, re_cols, rs_cols, re_cols];
            for k = 1:4
                rse[j+1][rr+k-1] = sub2ind((rows,cols),child_rows[k],child_cols[k])
            end
            rr = rr + 4

        else
            rse[j+1][rr:rr+1] = rse[j][r:r+1]
            rr = rr + 2
        end

    end
    j = j+1
    rse[j] = rse[j][rse[j].!=0]
    IND[:,:,j] = ind
end
return ind,rse,IND
end

"""
    u,v = second_largest_singular_vectors(A::Matrix{Float64})

Compute the second largest left and right singular vectors of matrix.

### Input Arguments
* `A`: The matrix of which the singular vectors will be computed.

### Output Arguments
* `u`: Second largest left singular vector.
* `v`: Second largest right singular vector.

Copyright 2018 The Regents of the University of California

Implemented by Yiqun Shao (Adviser: Dr. Naoki Saito)
"""
function second_largest_singular_vectors(A::Matrix{Float64})
#Compute the second largest left and right singular vectors of An = D1^-0.5*A*D2^-0.5

rows,cols = size(A)

# compute D1 and D2 and make sure there are no zeros
D1 = sum(A,2)
D1[D1 .== 0] = maximum([0.01, minimum(D1[D1 .> 0])/10])
D2 = sum(A,1)
D2[D2 .== 0] = maximum([0.01, minimum(D2[D2 .> 0])/10])

u,_,v = svd(Diagonal(D1[:].^(-0.5))*A*Diagonal(D2[:].^(-0.5)))

# extract the 2nd singular vectors and multiply by D_i^-0.5
u = Diagonal(D1[:].^-0.5)*u[:,2]
v = Diagonal(D2[:].^-0.5)*v[:,2]
return u,v
end

"""
    BBmatrix, BBmatrix_ind, tag, dmatrix = GHWT_Matrix2D_Analysis_BestBasis(matrix::Matrix{Float64},ind::Matrix{Int64} ,rse::Vector{Array{Int64,1}})

Compute the expanding coeffcients of all levels and search for the best-basis given the partition quadtree of matrix.

### Input Arguments
* `matrix`: The matrix which will be analyzed.
* `ind`: Vector with the same length as length(matrix[:]), indicate the order of matrix after partition.
* `rse`: Region start and end. Each pair represents the start and end of a region in the matrix, corresponding to one node in the quadtree

### Output Arguments
* `BBmatrix`: Coefficients of the best-basis.
* `BBmatrix_ind`: Location of the best-basis in dmatrix.
* `tag`: Tag information of the expanding coefficients.
* `dmatrix`: The expanding coefficients on all levels.

Copyright 2018 The Regents of the University of California

Implemented by Yiqun Shao (Adviser: Dr. Naoki Saito)
"""

function GHWT_Matrix2D_Analysis_BestBasis(matrix::Matrix{Float64},ind::Matrix{Int64} ,rse::Vector{Array{Int64,1}})
    # cost functional
    costfun = function(x) norm(x[:],1) end

    # constants
    rows,cols = size(matrix);
    #tol = 10^4*eps()

    # 1. Recursively partition the rows and columns of the matrix in a coupled manner
    #ind,rse,_ = PartitionTreeMatrixDhillon_QuadTree_ys(matrix)
    jmax = length(rse)

    # 2. Analyze the matrix and perform the best basis search
    dmatrix = zeros(rows,cols,jmax)
    dmatrix[:,:,jmax] = matrix[ind]
    tag = zeros(Int64,rows,cols,jmax)

    for r = 1:2:length(rse[jmax])
        if rse[jmax][r] != rse[jmax][r+1]
            rs_rows,rs_cols = ind2sub((rows,cols),rse[jmax][r])
            re_rows,re_cols = ind2sub((rows,cols),rse[jmax][r+1])

            tag[rs_rows:re_rows,rs_cols:re_cols,jmax] = reshape(0:(re_rows-rs_rows+1)*(re_cols-rs_cols+1)-1,re_rows-rs_rows+1,re_cols-rs_cols+1)
            average = sum(dmatrix[rs_rows:re_rows,rs_cols:re_cols,jmax])/sqrt((re_rows-rs_rows+1)*(re_cols-rs_cols+1))
            dmatrix[rs_rows:re_rows,rs_cols:re_cols,jmax] = 0
            dmatrix[rs_rows,rs_cols,jmax] = average
        end
    end
    #best basis
    BBmatrix = dmatrix[:,:,jmax]
    BBmatrix_ind = jmax*ones(Int64,rows,cols)

    for j = jmax-1:-1:1
        rr = 1;
        for r = 1:2:length(rse[j])
            rs_rows,rs_cols = ind2sub((rows,cols),rse[j][r])
            re_rows,re_cols = ind2sub((rows,cols),rse[j][r+1])

            # the row indices being considered
            indrs_rows = rs_rows:re_rows

            # the column indices being considered
            indrs_cols = rs_cols:re_cols

            n_rows = re_rows-rs_rows+1
            n_cols = re_cols-rs_cols+1

            #copy of children
            if rse[j][r+1] == rse[j+1][rr+1]
                dmatrix[indrs_rows,indrs_cols,j] = dmatrix[indrs_rows,indrs_cols,j+1]
                tag[indrs_rows,indrs_cols,j] = tag[indrs_rows,indrs_cols,j+1]
                rr = rr+2

                #copy of two childern
            elseif rse[j][r+1] == rse[j+1][rr+3]
                a,b = zeros(Int64,4), zeros(Int64,4)
                for k = 1:4
                    a[k],b[k] = ind2sub((rows,cols),rse[j+1][rr+k-1])
                end

                region1 = dmatrix[a[1]:a[2],b[1]:b[2],j+1]
                region2 = dmatrix[a[3]:a[4],b[3]:b[4],j+1]
                tag1 = tag[a[1]:a[2],b[1]:b[2],j+1]
                tag2 = tag[a[3]:a[4],b[3]:b[4],j+1]
                region_vec, tag_vec = twoblockcombine(region1[:],region2[:],tag1[:],tag2[:]);
                tag[indrs_rows,indrs_cols,j] = reshape(tag_vec,n_rows,n_cols);
                dmatrix[indrs_rows,indrs_cols,j] = reshape(region_vec,n_rows,n_cols)
                rr = rr + 4

            elseif rse[j][r+1] == rse[j+1][rr+7]
                a,b = zeros(Int64,8), zeros(Int64,8)
                for k = 1:8
                    a[k],b[k] = ind2sub((rows,cols),rse[j+1][rr+k-1])
                end

                region11 = dmatrix[a[1]:a[2],b[1]:b[2],j+1]
                region21 = dmatrix[a[3]:a[4],b[3]:b[4],j+1]
                tag11 = tag[a[1]:a[2],b[1]:b[2],j+1]
                tag21 = tag[a[3]:a[4],b[3]:b[4],j+1]
                region12 = dmatrix[a[5]:a[6],b[5]:b[6],j+1]
                region22 = dmatrix[a[7]:a[8],b[7]:b[8],j+1]
                tag12 = tag[a[5]:a[6],b[5]:b[6],j+1]
                tag22 = tag[a[7]:a[8],b[7]:b[8],j+1]

                region_vec_1, tag_vec_1 = twoblockcombine(region11[:],region21[:],tag11[:],tag21[:])
                region_vec_2, tag_vec_2 = twoblockcombine(region12[:],region22[:],tag12[:],tag22[:])
                region_vec, tag_vec = twoblockcombine(region_vec_1,region_vec_2,tag_vec_1,tag_vec_2)

                tag[indrs_rows,indrs_cols,j] = reshape(tag_vec, n_rows,n_cols)
                dmatrix[indrs_rows,indrs_cols,j] = reshape(region_vec,n_rows,n_cols)
                rr = rr + 8

            end
            #Obtaining BBmatrix
            if costfun(BBmatrix[indrs_rows,indrs_cols]) > costfun(dmatrix[indrs_rows,indrs_cols,j])
                BBmatrix[indrs_rows,indrs_cols] = dmatrix[indrs_rows,indrs_cols,j]
                BBmatrix_ind[indrs_rows,indrs_cols] = j;

            end
        end
    end
    return BBmatrix, BBmatrix_ind, tag, dmatrix
end

"""
    region, tag = twoblockcombine(region1::Vector{Float64}, region2::Vector{Float64}, tag1::Vector{Int64}, tag2::Vector{Int64})

Compute the coefficients of coarser basis-vectors after combining finer basis-vectors supported on two sub-regions.

### Input Arguments
* `region1`: The coeffcient of finer basis-vectors supported on region 1.
* `tag1`: The tag information of finer basis-vectors supported on region 1.
* `region2`: The coeffcient of finer basis-vectors supported on region 2.
* `tag2`: The tag information of finer basis-vectors supported on region 2.

### Output Arguments
* `region`: The coeffcient of coarser basis-vectors supported on region 1 and region 2.
* `tag`: The tag information of coarser basis-vectors supported on region 1 and region 2.

Copyright 2018 The Regents of the University of California

Implemented by Yiqun Shao (Adviser: Dr. Naoki Saito)
"""

function twoblockcombine(region1::Vector{Float64}, region2::Vector{Float64}, tag1::Vector{Int64}, tag2::Vector{Int64})
    # the number of points in the first subregion
    n1 = length(region1)
    # the number of points in the second subregion
    n2 = length(region2)

    n = n1+n2
    region = zeros(n)
    tag = zeros(Int64,n)

    ### SCALING COEFFICIENT (n > 1)
    region[1] = ( sqrt(n1)*region1[1] + sqrt(n2)*region2[1] )/sqrt(n)
    tag[1] = 0

    ### HAAR-LIKE COEFFICIENT
    region[2] = ( sqrt(n2)*region1[1] - sqrt(n1)*region2[1] )/sqrt(n);
    tag[2] = 1

    ### WALSH-LIKE COEFFICIENTS
    # sweep through the coefficients in subregion 1 and subregion 2

    # the index of the new coefficient(s) to be created on level j
    parent = 3

    # the index of the current coefficients in subregions 1 and 2
    child1 = 2
    child2 = 2

    while parent <= n
        # no matching coefficient (use subregion 1)
        if child1 <= n1 && (child2 == n2+1 || tag1[child1] < tag2[child2])
            region[parent] = region1[child1]
            tag[parent] = 2*tag1[child1]
            child1 = child1+1;
            parent = parent+1;

            # no matching coefficient (use subregion 2)
        elseif child2 <=n2 && (child1 == n1+1 || tag2[child2] < tag1[child1])
            region[parent] = region2[child2]
            tag[parent] = 2*tag2[child2]
            child2 = child2+1
            parent = parent+1

            # matching coefficients
        else
            region[parent] = ( region1[child1] + region2[child2] )/sqrt(2)
            tag[parent] = 2*tag[child1]

            region[parent+1] = ( region1[child1] - region2[child2] )/sqrt(2)
            tag[parent+1] = 2*tag[child1]+1

            child1 = child1+1
            child2 = child2+1
            parent = parent+2
        end
    end
    return region, tag
end

"""
    matrix_r = GHWT_Matrix2D_Analysis_Synthesis(BBmatrix::Matrix{Float64},BBmatrix_ind::Matrix{Int64},tag::Array{Int64,3} ,rse::Vector{Array{Int64,1}})

Synthesize the matrix given coefficients of selected basis-vectors.

### Input Arguments
* `BBmatrix`: Coefficients of the best-basis.
* `BBmatrix_ind`: Location of the best-basis in dmatrix.
* `tag`: Tag information of the expanding coefficients.
* `rse`: Region start and end. Each pair represents the start and end of a region in the matrix, corresponding to one node in the quadtree.

### Output Arguments
* `matrix_r`: The synthesized matrix.

Copyright 2018 The Regents of the University of California

Implemented by Yiqun Shao (Adviser: Dr. Naoki Saito)
"""
function GHWT_Matrix2D_Analysis_Synthesis(BBmatrix::Matrix{Float64},BBmatrix_ind::Matrix{Int64},tag::Array{Int64,3} ,rse::Vector{Array{Int64,1}})
    costfun = function(x) norm(x[:],1) end
    # constants
    rows,cols = size(BBmatrix);
    #tol = 10^4*eps()
    jmax = length(rse)

    cmatrix = zeros(rows,cols,jmax);
    for i  = 1:rows
        for j = 1:cols
            cmatrix[i,j,BBmatrix_ind[i,j]] = BBmatrix[i,j];
        end
    end

    for j = 1:1:jmax-1
        rr = 1
        for r = 1:2:length(rse[j])
            rs_rows,rs_cols = ind2sub((rows,cols),rse[j][r])
            re_rows,re_cols = ind2sub((rows,cols),rse[j][r+1])

            if rse[j][r+1] == rse[j+1][rr+1]
                if countnz(cmatrix[rs_rows:re_rows,rs_cols:re_cols,j+1]) == 0 && countnz(cmatrix[rs_rows:re_rows,rs_cols:re_cols,j]) > 0
                    cmatrix[rs_rows:re_rows,rs_cols:re_cols,j+1] = cmatrix[rs_rows:re_rows,rs_cols:re_cols,j]
                end
                rr = rr + 2

            elseif rse[j][r+1] == rse[j+1][rr+3]
                if countnz(cmatrix[rs_rows:re_rows,rs_cols:re_cols,j+1]) == 0 && countnz(cmatrix[rs_rows:re_rows,rs_cols:re_cols,j]) > 0
                    a,b = zeros(Int64,4), zeros(Int64,4)
                    for k = 1:4
                        a[k],b[k] = ind2sub((rows,cols),rse[j+1][rr+k-1])
                    end

                    region = cmatrix[rs_rows:re_rows,rs_cols:re_cols,j]
                    tag1 = tag[a[1]:a[2],b[1]:b[2],j+1]
                    tag2 = tag[a[3]:a[4],b[3]:b[4],j+1]
                    region1_vec, region2_vec = twoblocksplit(region[:],tag1[:],tag2[:]);
                    cmatrix[a[1]:a[2],b[1]:b[2],j+1] = reshape(region1_vec,a[2]-a[1]+1,b[2]-b[1]+1)
                    cmatrix[a[3]:a[4],b[3]:b[4],j+1] = reshape(region2_vec,a[4]-a[3]+1,b[4]-b[3]+1)
                end
                rr = rr + 4
            else
                if countnz(cmatrix[rs_rows:re_rows,rs_cols:re_cols,j+1]) == 0 && countnz(cmatrix[rs_rows:re_rows,rs_cols:re_cols,j]) > 0
                    a,b = zeros(Int64,8), zeros(Int64,8)
                    for k = 1:8
                        a[k],b[k] = ind2sub((rows,cols),rse[j+1][rr+k-1])
                    end
                    tag11 = tag[a[1]:a[2],b[1]:b[2],j+1]
                    tag21 = tag[a[3]:a[4],b[3]:b[4],j+1]
                    tag12 = tag[a[5]:a[6],b[5]:b[6],j+1]
                    tag22 = tag[a[7]:a[8],b[7]:b[8],j+1]
                    tag11 = tag11[:]
                    tag21 = tag21[:]
                    tag12 = tag12[:]
                    tag22 = tag22[:]
                    tag1 = twoblockcombine_tag(tag11,tag21)
                    tag2 = twoblockcombine_tag(tag12,tag22)
                    region = cmatrix[rs_rows:re_rows,rs_cols:re_cols,j]
                    region1_vec,region2_vec = twoblocksplit(region[:], tag1,tag2)
                    region11,region21 = twoblocksplit(region1_vec,tag11,tag21)
                    region12,region22 = twoblocksplit(region2_vec,tag12,tag22)
                    cmatrix[a[1]:a[2],b[1]:b[2],j+1] = reshape(region11,a[2]-a[1]+1,b[2]-b[1]+1)
                    cmatrix[a[3]:a[4],b[3]:b[4],j+1] = reshape(region21,a[4]-a[3]+1,b[4]-b[3]+1)
                    cmatrix[a[5]:a[6],b[5]:b[6],j+1] = reshape(region12,a[6]-a[5]+1,b[6]-b[5]+1)
                    cmatrix[a[7]:a[8],b[7]:b[8],j+1] = reshape(region22,a[8]-a[7]+1,b[8]-b[7]+1)
                end
                rr = rr + 8
            end
        end
    end

    matrix_r = cmatrix[:,:,end]
    for r = 1:2:length(rse[jmax])
        if rse[jmax][r] != rse[jmax][r+1]
            rs_rows,rs_cols = ind2sub((rows,cols),rse[j][r])
            re_rows,re_cols = ind2sub((rows,cols),rse[j][r+1])

            matrix_r[rs_rows:re_rows,rs_cols:re_cols] = matrix_r[rs_rows,rs_cols]/sqrt((re_rows-rs_rows+1)*(re_cols-rs_cols+1))
        end
    end
    return matrix_r
end

"""
    tag = twoblockcombine_tag(tag1::Vector{Int64},tag2::Vector{Int64})

Compute the tag information of coarser basis-vectors after combining finer basis-vectors supported on two sub-regions.

### Input Arguments
* `tag1`: The tag information of finer basis-vectors supported on region 1.
* `tag2`: The tag information of finer basis-vectors supported on region 2.

### Output Arguments
* `tag`: The tag information of coarser basis-vectors supported on region 1 and region 2.

Copyright 2018 The Regents of the University of California

Implemented by Yiqun Shao (Adviser: Dr. Naoki Saito)
"""
function twoblockcombine_tag(tag1::Vector{Int64},tag2::Vector{Int64})
    # the number of points in the first subregion
    n1 = length(tag1)
    # the number of points in the second subregion
    n2 = length(tag2)

    n = n1+n2
    tag = zeros(Int64,n)

    ### SCALING COEFFICIENT (n > 1)
    tag[1] = 0

    ### HAAR-LIKE COEFFICIENT
    tag[2] = 1

    ### WALSH-LIKE COEFFICIENTS
    # sweep through the coefficients in subregion 1 and subregion 2

    # the index of the new coefficient(s) to be created on level j
    parent = 3

    # the index of the current coefficients in subregions 1 and 2
    child1 = 2
    child2 = 2

    while parent <= n
        # no matching coefficient (use subregion 1)
        if child1 <= n1 && (child2 == n2+1 || tag1[child1] < tag2[child2])
            tag[parent] = 2*tag1[child1]
            child1 = child1+1
            parent = parent+1

            # no matching coefficient (use subregion 2)
        elseif child2 <=n2 && (child1 == n1+1 || tag2[child2] < tag1[child1])
            tag[parent] = 2*tag2[child2]
            child2 = child2+1
            parent = parent+1

            # matching coefficients
        else
            tag[parent] = 2*tag[child1]
            tag[parent+1] = 2*tag[child1]+1

            child1 = child1+1
            child2 = child2+1
            parent = parent+2
        end
    end
    return tag
end

"""
    region1,region2 = twoblocksplit(region::Matrix{Float64}, tag1::Vector{Int64}, tag2::Vector{Int64})

Compute the coefficients of finer basis-vectors supported on two sub-regions given the coeffcients of coarser basis-vectors supported on both regions.

### Input Arguments
* `region`: The coeffcient of coarser basis-vectors supported on region 1 and region 2.
* `tag1`: The tag information of finer basis-vectors supported on region 1.
* `tag2`: The tag information of finer basis-vectors supported on region 2.

### Output Arguments
* `region1`: The coeffcient of finer basis-vectors supported on region 1.
* `region2`: The coeffcient of finer basis-vectors supported on region 2.

Copyright 2018 The Regents of the University of California

Implemented by Yiqun Shao (Adviser: Dr. Naoki Saito)
"""
function twoblocksplit(region::Matrix{Float64}, tag1::Vector{Int64}, tag2::Vector{Int64})
    # the number of points in the first subregion
    n1 = length(tag1)
    # the number of points in the second subregion
    n2 = length(tag2)

    n = n1+n2
    region1 = zeros(n1)
    region2 = zeros(n2);

    ### SCALING COEFFICIENT (n > 1)
    region1[1] = ( sqrt(n1)*region[1] + sqrt(n2)*region[2] )/sqrt(n);

    ### HAAR-LIKE COEFFICIENT
    region2[1] = ( sqrt(n2)*region[1] - sqrt(n1)*region[2] )/sqrt(n);

    ### WALSH-LIKE COEFFICIENTS
    # sweep through the coefficients in subregion 1 and subregion 2

    # the index of the new coefficient(s) to be created on level j
    parent = 3

    # the index of the current coefficients in subregions 1 and 2
    child1 = 2
    child2 = 2

    while parent <= n
        # no matching coefficient (use subregion 1)
        if child1 <= n1 && (child2 == n2+1 || tag1[child1] < tag2[child2])
            region1[child1] = region[parent]
            child1 = child1+1
            parent = parent+1

            # no matching coefficient (use subregion 2)
        elseif child2 <=n2 && (child1 == n1+1 || tag2[child2] < tag1[child1])
            region2[child2] = region[parent]
            child2 = child2+1
            parent = parent+1

            # matching coefficients
        else
            region1[child1] = ( region[parent] + region[parent+1] )/sqrt(2)
            region2[child2] = ( region[parent] - region[parent+1] )/sqrt(2)
            child1 = child1+1
            child2 = child2+1
            parent = parent+2
        end

    end
    return region1,region2
end


text:
module PerformanceIndex

using LinearAlgebra
using Reexport

include("IAE.jl")
include("ISE.jl")
include("ITAE.jl")
include("ITSE.jl")

@reexport using .IAE
@reexport using .ISE
@reexport using .ITAE
@reexport using .ITSE

end # module


text:
#name: JuDecTreeTrain
#meta.mlname: JuDecTree
#meta.mlrole: train
#description: Custom Julia Train func for Decision Tree
#language: julia
#input: dataframe df
#input: string predict_column
#input: int max_depth=2 {category: Parameters}
#output: blob model

#>>> Import necessary packages
using Pkg
Pkg.add("DecisionTree")
using DecisionTree
Pkg.add("BSON")
using BSON: @save

#>>> Extract/Prepare train features and target variable
trainX=Matrix(select(df, Not(Symbol(predict_column))))
trainY=Matrix(select(df, Symbol(predict_column)))
trainY=vec(trainY)

#>>> Build/train model
DTCmodel = DecisionTreeClassifier(max_depth=max_depth)
trained_model = fit!(DTCmodel, trainX, trainY)

#>>> Save trained model
@save file=model trained_model

max_stars_repo_name: DotNetPart/public

max_issues_repo_name: DotNetPart/public

max_forks_repo_name: DotNetPart/public


text:
abstract type AbstractSource end

import Base: @kwdef

macro defsource(ex)
    nex = ex.head == :macrocall ? ex.args[end] : ex
    defsource(nex) |> esc
end

function defsource(ex)
    head = ex.head
    ismutable, name, args = ex.args
    if ismutable
        if name isa Symbol
            return quote
                Base.@kwdef mutable struct $name{TR, HS, CB} <: AbstractSource
                    $args
                    trigger::TR = 1
                    handshake::HS = 2
                    callbacks::CB = nothing
                    name::Symbol = Symbol()
                end
            end
        else
            return quote
                Base.@kwdef mutable struct $(name.args[1]){$(name.args[2:end]...), TR, HS, CB} <: AbstractSource
                    $args
                    trigger::TR = 1
                    handshake::HS = 2
                    callbacks::CB = nothing
                    name::Symbol = Symbol()
                end
            end
        end
    else
        if name isa Symbol
            return quote
                Base.@kwdef struct $name{TR, HS, CB} <: AbstractSource
                    $args
                    trigger::TR = 1
                    handshake::HS = 2
                    callbacks::CB = nothing
                    name::Symbol = Symbol()
                end
            end
        else
            return quote
                Base.@kwdef struct $(name.args[1]){$(name.args[2:end]...), TR, HS, CB} <: AbstractSource
                    $args
                    trigger::TR = 1
                    handshake::HS = 2
                    callbacks::CB = nothing
                    name::Symbol = Symbol()
                end
            end
        end
    end
end

@macroexpand Base.@kwdef struct Sinewave7
    amplitude::Int = 1
    phase::Int
    gain = 1
end


text:
# This file was generated, do not modify it. # hide
booster = Booster()

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
abstract type AbstractEmpiricalValue <: AbstractUncertainValue end

function summarise(ed::AbstractEmpiricalValue)
    _type = typeof(ed)
    disttype = typeof(ed.distribution)

    l = length(ed.values)
    summary = "$_type estimated as a $disttype from $l values"
    return summary
end

Base.show(io::IO, ed::AbstractEmpiricalValue) = println(io, summarise(ed))

@recipe function plot_empiricalval(empval::AbstractEmpiricalValue;
        nbins = 200, n_samples = 10000)
    dist = empval.distribution
    @series begin
        seriestype := :histogram
        fit(Histogram, rand(dist, n_samples), nbins = nbins)
    end
end

export
AbstractEmpiricalValue


text:
import Random: randperm
import Base: ==, <, isless

mutable struct NNTuple{R, S}
    idx::R
    dist::S
    flag::Bool
end
NNTuple(a, b) = NNTuple(a, b, true)

==(a::NNTuple, b::NNTuple) = (a.idx == b.idx) && (a.dist == b.dist) && (a.flag == b.flag)
<(a::NNTuple, b::NNTuple) = a.dist < b.dist
isless(a::NNTuple, b::NNTuple) = <(a, b)

"""
Remove the `k` nearest neighbors from each heap in `knn_heaps`.
Return two k x length(knn_heaps) arrays for the indices and
distances to each point's kNN.
"""
function deheap_knns(heaps::Vector{H}, k) where {S, T,
                                                 H <: AbstractHeap{NNTuple{S, T}}}

    ids = Array{S}(undef, (k, length(heaps)))
    dists = Array{T}(undef, (k, length(heaps)))

    for i in 1:length(heaps)
        len = length(heaps[i])
        for j in 1:len
            # NOTE: these are max heaps, so we only want the last k
            tuple = pop!(heaps[i])
            neighbor_idx = 1 + len - j
            if neighbor_idx <= k
                ids[neighbor_idx, i] = tuple.idx
                dists[neighbor_idx, i] = tuple.dist
            end
        end
    end
    return ids, dists
end

"""
Sample `n_neighbors` elements from a set of ints `1:npoints`.
The ints in `exclude` won't be sampled.
"""
function sample_neighbors(npoints::Int,
                          n_neighbors::Int,
                          sample_rate::R = 1.;
                          exclude::Vector{Int} = Vector{Int}()) where {R <: AbstractFloat}
    # num neighbors to sample = max(ρK, N)
    last = min(npoints-length(exclude), trunc(Int, sample_rate*n_neighbors))
    ids = Set{Int}()
    while length(ids) < last
        k = (abs(rand(Int)) % npoints) + 1
        if k ∉ exclude
            union!(ids, k)
        end
    end
    return collect(ids)
end

function make_knn_heaps(data::Vector{V},
                        n_neighbors::Int,
                        metric::M) where {V <: AbstractArray,
                                          M <: SemiMetric}
    np = length(data)
    D = result_type(metric, data[1], data[1])
    HeapType = BinaryMaxHeap{NNTuple{Int, D}}
    knn_heaps = HeapType[HeapType() for _ in 1:np]
    for i in 1:np
        k_idxs = sample_neighbors(np, n_neighbors, exclude=[i])
        for j in k_idxs
            d = evaluate(metric, data[i], data[j])
            _heappush!(knn_heaps[i], NNTuple(j, d), n_neighbors)
            _heappush!(knn_heaps[j], NNTuple(i, d), n_neighbors)
        end
    end
    return knn_heaps
end

"""
Get the neighbors of each point in a KNN graph, `graph`,
as an array of ids.
"""
function _neighbors(graph, sample_rate::AbstractFloat = 1.)
    old_fw_neighbors = [Vector{Int}() for _ in 1:length(graph)]
    new_fw_neighbors = [Vector{Int}() for _ in 1:length(graph)]
    old_bw_neighbors = [Vector{Int}() for _ in 1:length(graph)]
    new_bw_neighbors = [Vector{Int}() for _ in 1:length(graph)]
    for i in 1:length(graph), j in 1:length(graph[i])
        # add incoming edge ith -> jth.idx
        if rand() ≤ sample_rate
            if graph[i].valtree[j].flag
                # denote this neighbor has participated in local join
                graph[i].valtree[j].flag = false
                append!(new_fw_neighbors[i], graph[i].valtree[j].idx)
                append!(new_bw_neighbors[graph[i].valtree[j].idx], i)
            else
                append!(old_fw_neighbors[i], graph[i].valtree[j].idx)
                append!(old_bw_neighbors[graph[i].valtree[j].idx], i)
            end
        end
    end
    return old_fw_neighbors, new_fw_neighbors, old_bw_neighbors, new_bw_neighbors
end

function min_flagged(heap)
    min_tup = NNTuple(-1, typemax(typeof(heap.valtree[1].dist)))
    for t in heap.valtree
        if t < min_tup && !t.flag
            min_tup = t
        end
    end
    return min_tup
end

"""
    _heappush!(heap::BinaryHeap, tup::NNTuple, max_candidates)

Try to push a neighbor `tup` to `heap`. This will fail (return `0`) if `tup` is
already in `heap`, if `tup.dist > top(heap).dist`. Otherwise return `1`.
If `length(heap) > max_candidates` after pushing, `pop` the largest candidate.
"""
function _heappush!(heap::AbstractHeap,
                    tup::NNTuple,
                    max_candidates::Int=length(heap))

    if max_candidates == 0
        @debug "max_candidates has a size of 0"
        return 0
    # case: empty heap
    elseif length(heap) == 0
        push!(heap, tup)
        return 1
    elseif length(heap) < max_candidates || tup < top(heap)
        # check if already in heap
        for i in 1:length(heap)
            if _check_tuple(heap, i, tup)
                return 0
            end
        end
        # push and maintain size
        if length(heap) == max_candidates
            _push_or_update!(heap, tup, max_candidates)
        else
            push!(heap, tup)
        end
        return 1
    end
    return 0
end

"""
Push `tup` onto `heap` without checkout if it already exists in the heap.
"""
function _unchecked_heappush!(heap::AbstractHeap,
                              tup::NNTuple,
                              max_candidates::Int=length(heap))
    if max_candidates == 0
        @debug "max_candidates has a size of 0"
        return 0
    end
    # case: heap not full
    if length(heap) < max_candidates
        push!(heap, tup)
        return 1
    # heap full but this neighbor closer
    elseif tup < top(heap)
        # push and maintain size
        _push_or_update!(heap, tup, max_candidates)
        return 1
    end
    return 0

end

@inline function _push_or_update!(h::BinaryHeap, t, maxc)
    push!(h, t)
    if length(h) > maxc
        pop!(h)
    end
    return
end

"""
Check if a tuple exists in a heap at index `i`, and optionally update its dist.
Returns (exists::Bool, updated::Bool)
"""
@inline function _check_tuple(h::BinaryHeap, i, t)
    if h.valtree[i].idx == t.idx
        return true
    end
    return false
end


text:
using RadialBasisFunctionModels

using Test

f = x -> x^2

X = collect( LinRange(-4,4,5) )
Y = f.(X)

φ = Multiquadric()

rbf = RBFModel( X, Y, φ, 1)

Z = rbf.(X)

@test typeof(Z[1]) <: AbstractVector{<:Real}
@test length(Z[1]) == 1

@test all( isapprox(Z[i][1], Y[i]; atol = 1e-10) for i = 1 : length(Z) )

@test rbf( [ X[1], ] ) == Z[1]

rbf_scalar = RBFInterpolationModel( X, Y, φ, 1; vector_output = false)
Z_scalar = rbf_scalar.( X )

@test(#jl
Z_scalar isa Vector{Float64}
)#jl
@test(
all( Z_scalar[i] == Z[i][1] for i = 1 : length(Z) )
)

X_f0 = Float32.(X)
Y_f0 = f.(X_f0)
rbf_f0 = RBFInterpolationModel( X_f0, Y_f0, φ, 1)
@test(#jl
rbf_f0.(X_f0) isa Vector{Vector{Float32}}
)#jl

using StaticArrays
features = [ @SVector(rand(3)) for i = 1 : 5 ]
labels = [ @SVector(rand(3)) for i = 1 : 5 ]
centers = SVector{5}(features)
rbf_sized = RBFModel( features, labels; centers )

x_vec = rand(3)
x_s = SVector{3}(x_vec)
x_m = MVector{3}(x_vec)
x_sized = SizedVector{3}(x_vec)

@test(#jl
rbf_sized( x_vec ) isa Vector
)#jl
@test(#jl
rbf_sized( x_s ) isa SVector
)#jl
@test(#jl
rbf_sized( x_m ) isa MVector
)#jl
@test(#jl
rbf_sized( x_sized ) isa SizedVector
)#jl

RBFModel(features, labels, :gaussian)

RBFModel(features, labels, :multiquadric, [1.0, 1//2])

using MLJBase
X,y = @load_boston

r = RBFInterpolator(; kernel_name = :multiquadric )
R = machine(r, X, y)

MLJBase.fit!(R)
MLJBase.predict(R, X)

X = [ rand(2) for i = 1 : 10 ]
Y = [ rand(2) for i = 1 : 10 ]

R = RBFMachine(;features = X, labels = Y, kernel_name = :gaussian )
R isa RBFMachineWithKernel
RadialBasisFunctionModels.fit!(R)
@test(#jl
R( X[1] ) ≈ Y[1]
)#jl

R = RBFMachine()
add_data!(R, X, Y)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
# The "pound" symbol or "hash" symbol or "octothorp?" is the comment character
# in a julia file. So all of these commands are ignored!

##
grades = [
  30
  40
  30
  45
  50
  20
  35
  40
  50
  49
  55
]
avggrade = sum(grades)/length(grades)
##
println("The average grade is ", avggrade)

##
# Let's compute the standard deviation too!
# The formula for the standard deviation is
# sqrt(mean((x_i - mean(x)))^2)
# Julia has a built in stdev function, we'll use that in a second.

cgrades = grades .- avggrade
stdgrades = sqrt(sum(cgrades.^2)/length(cgrades))

##
using Statistics
##
mean(grades)
##
std(grades; corrected=false)


text:
using Core.Compiler: AbstractInterpreter, NativeInterpreter, InferenceState,
    OptimizationState, CodeInfo, CodeInstance, InferenceResult, WorldRange,
    IRCode, SSAValue, inlining_policy

struct InferredSource
    src::CodeInfo
    stmt_infos::Vector{Any}
    rt::Any
end

struct OptimizedSource
    ir::IRCode
    isinlineable::Bool
end

mutable struct CthulhuInterpreter <: AbstractInterpreter
    native::NativeInterpreter

    unopt::Dict{Union{MethodInstance, InferenceResult}, InferredSource}
    opt::Dict{MethodInstance, CodeInstance}

    msgs::Dict{MethodInstance, Vector{Pair{Int, String}}}
end

CthulhuInterpreter() = CthulhuInterpreter(
    NativeInterpreter(),
    Dict{MethodInstance, InferredSource}(),
    Dict{MethodInstance, CodeInstance}(),
    Dict{MethodInstance, Vector{Pair{Int, String}}}()
)

import Core.Compiler: InferenceParams, OptimizationParams, get_world_counter,
    get_inference_cache, code_cache,
    WorldView, lock_mi_inference, unlock_mi_inference, InferenceState
using Base: @invoke

Compiler.InferenceParams(interp::CthulhuInterpreter) = InferenceParams(interp.native)
Compiler.OptimizationParams(interp::CthulhuInterpreter) = OptimizationParams(interp.native)
Compiler.get_world_counter(interp::CthulhuInterpreter) = get_world_counter(interp.native)
Compiler.get_inference_cache(interp::CthulhuInterpreter) = get_inference_cache(interp.native)

# No need to do any locking since we're not putting our results into the runtime cache
Compiler.lock_mi_inference(interp::CthulhuInterpreter, mi::MethodInstance) = nothing
Compiler.unlock_mi_inference(interp::CthulhuInterpreter, mi::MethodInstance) = nothing

struct CthulhuCache
    cache::Dict{MethodInstance, CodeInstance}
end

Compiler.code_cache(interp::CthulhuInterpreter) = WorldView(CthulhuCache(interp.opt), WorldRange(get_world_counter(interp)))
Compiler.get(wvc::WorldView{CthulhuCache}, mi::MethodInstance, default) = get(wvc.cache.cache, mi, default)
Compiler.haskey(wvc::WorldView{CthulhuCache}, mi::MethodInstance) = haskey(wvc.cache.cache, mi)
Compiler.setindex!(wvc::WorldView{CthulhuCache}, ci::CodeInstance, mi::MethodInstance) = setindex!(wvc.cache.cache, ci, mi)

Compiler.may_optimize(interp::CthulhuInterpreter) = true
Compiler.may_compress(interp::CthulhuInterpreter) = false
Compiler.may_discard_trees(interp::CthulhuInterpreter) = false
Compiler.verbose_stmt_info(interp::CthulhuInterpreter) = true

function Compiler.add_remark!(interp::CthulhuInterpreter, sv::InferenceState, msg)
    push!(get!(interp.msgs, sv.linfo, Tuple{Int, String}[]),
        sv.currpc => msg)
end

function Compiler.finish(state::InferenceState, interp::CthulhuInterpreter)
    r = @invoke Compiler.finish(state::InferenceState, interp::AbstractInterpreter)
    interp.unopt[Core.Compiler.any(state.result.overridden_by_const) ? state.result : state.linfo] = InferredSource(
        copy(isa(state.src, OptimizationState) ?
            state.src.src : state.src),
        copy(state.stmt_info),
        state.result.result)
    return r
end

function Compiler.transform_result_for_cache(interp::CthulhuInterpreter, linfo::MethodInstance,
        valid_worlds::Core.Compiler.WorldRange, @nospecialize(inferred_result))
    if isa(inferred_result, OptimizationState)
        opt = inferred_result
        if isdefined(opt, :ir)
            return OptimizedSource(opt.ir::IRCode, opt.src.inlineable)
        end
    end
    return inferred_result
end

function Compiler.inlining_policy(interp::CthulhuInterpreter)
    function (src)
        @assert isa(src, OptimizedSource)
        return src.isinlineable ? src.ir : nothing
    end
end

function Compiler.finish!(interp::CthulhuInterpreter, caller::InferenceResult)
    if isa(caller.src, OptimizationState)
        opt = caller.src
        if isdefined(opt, :ir)
            caller.src = OptimizedSource(opt.ir, opt.src.inlineable)
        end
    end
end


text:
using RNGPool

include("simpleModel.jl")
include("demo.jl")

import Statistics: mean, var

setRNGs(12345)

println("Using a Simple Model with n = 100")
model, lM = setupSimpleModel(100, 1.0)
vs = CoupledConditionalSMC.couplingTimes(model, lM, 32, model.maxn, 100, :BS, false)
println(mean(vs), ", ", var(vs))

println("Using a Simple Model with n = 200")
model, lM = setupSimpleModel(200, 1.0)
vs = CoupledConditionalSMC.couplingTimes(model, lM, 32, model.maxn, 100, :BS, false)
println(mean(vs), ", ", var(vs))

println("Using a Simple Model with n = 400")
model, lM = setupSimpleModel(400, 1.0)
vs = CoupledConditionalSMC.couplingTimes(model, lM, 32, model.maxn, 100, :BS, false)
println(mean(vs), ", ", var(vs))


text:
using FunctionTabulations
using Documenter

DocMeta.setdocmeta!(FunctionTabulations, :DocTestSetup, :(using FunctionTabulations); recursive=true)

makedocs(;
    modules=[FunctionTabulations],
    authors="Aurelio Amerio <aurelio.amerio@edu.unito.it> and contributors",
    repo="https://github.com/aurelio-amerio/FunctionTabulations.jl/blob/{commit}{path}#{line}",
    sitename="FunctionTabulations.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "List of functions" => "functions.md",
    ],
)


text:
using Test
using Elemental
using LinearAlgebra

m, n = 10, 10

for T in (Float32, Float64, Complex{Float32}, Complex{Float64})
    for mat in (Elemental.Matrix, Elemental.DistMatrix)
        A = mat(T)
        Elemental.gaussian!(A, m, n)

        @test opnorm(A, 1)   ≈ opnorm(T[A[i,j] for i = 1:m, j = 1:n], 1)
        @test opnorm(A)      ≈ opnorm(T[A[i,j] for i = 1:m, j = 1:n])
        @test opnorm(A, Inf) ≈ opnorm(T[A[i,j] for i = 1:m, j = 1:n], Inf)
    end
end


text:
export step_computation

function step_computation(direction :: String,
                          h :: AbstractNLPModel,
                          t :: Float64,
						  d :: Float64,
						  gₖ :: Float64,
						  fₖ :: Float64,
						  ftestTR :: Float64,
						  gtestTR :: Float64)

  if direction == "Nwt"
    t = t + d; fₖ = ftestTR; gₖ = gtestTR; H = hess(h, [t])[1]
    return (t, fₖ, gₖ, H)
  end
  if direction == "Sec"
	tpred = t; gₖₘ₁ = gₖ; t = t + d; fₖ = ftestTR
	gₖ = gtestTR; s = t - tpred; y = gₖ - gₖₘ₁; secₖ = y / s
	return (t, fₖ, gₖ, secₖ)
  end
  if direction == "SecA"
	tpred = t; gkm1 = gₖ; fkm1 = fₖ
	t = t + d; gₖ = gtestTR; fₖ = ftestTR
 	s = t - tpred; y = gₖ - gkm1
	Γ = 3 * (gₖ + gkm1) * s - 6 * (fₖ - fkm1)
	if (y * s + Γ) < eps(Float64) * (s^2) #correction trop petite
	  yt = y
	else
	  yt = y + Γ / s
	end
	secₖ = yt / s
  	return (t, fₖ, gₖ, secₖ)
  end
end


text:
import Base.Math: @horner

const FORWARD  =  true
const BACKWARD = false

const sqrtpi = 1.772453850905516027298
const edivsqrt2pi = 1.084437551419227546612

"""
Compute a typed 0.5.
"""
half(x::Number) = oftype(x,0.5)
half(x::Integer) = half(float(x))
half(::Type{T}) where {T<:Number} = convert(T,0.5)
half(::Type{T}) where {T<:Integer} = half(AbstractFloat)

"""
Compute a typed 2.
"""
two(x::Number) = oftype(x,2)
two(::Type{T}) where {T<:Number} = convert(T,2)

"""
The Kronecker ``\\delta`` function:

```math
\\delta_{k,j} = \\left\\{\\begin{array}{ccc} 1 & {\\rm for} & k = j,\\\\ 0 & {\\rm for} & k \\ne j.\\end{array}\\right.
```
"""
δ(k::Integer,j::Integer) = k == j ? 1 : 0

"""
Pochhammer symbol ``(x)_n = \\frac{\\Gamma(x+n)}{\\Gamma(x)}`` for the rising factorial.
"""
function pochhammer(x::Number,n::Integer)
    ret = one(x)
    if n≥0
        for i=0:n-1
            ret *= x+i
        end
    else
        ret /= pochhammer(x+n,-n)
    end
    ret
end

pochhammer(x::Number,n::Number) = isinteger(n) ? pochhammer(x,Int(n)) : ogamma(x)/ogamma(x+n)

function pochhammer(x::Number,n::UnitRange{T}) where T<:Real
    ret = Vector{promote_type(typeof(x),T)}(length(n))
    ret[1] = pochhammer(x,first(n))
    for i=2:length(n)
        ret[i] = (x+n[i]-1)*ret[i-1]
    end
    ret
end

ogamma(x::Number) = (isinteger(x) && x<0) ? zero(float(x)) : inv(gamma(x))

"""
Stirling's asymptotic series for ``\\Gamma(z)``.
"""
stirlingseries(z) = gamma(z)*sqrt((z/π)/2)*exp(z)/z^z

function stirlingseries(z::Float64)
    if z ≥ 3274.12075200175       # N =  4
        @horner(inv(z),1.0,0.08333333333333333,0.003472222222222222,-0.0026813271604938273)
    elseif z ≥ 590.1021805526798  # N =  5
        @horner(inv(z),1.0,0.08333333333333333,0.003472222222222222,-0.0026813271604938273,-0.00022947209362139917)
    elseif z ≥ 195.81733962412835 # N =  6
        @horner(inv(z),1.0,0.08333333333333333,0.003472222222222222,-0.0026813271604938273,-0.00022947209362139917,0.0007840392217200666)
    elseif z ≥ 91.4692823071966   # N =  7
        @horner(inv(z),1.0,0.08333333333333333,0.003472222222222222,-0.0026813271604938273,-0.00022947209362139917,0.0007840392217200666,6.972813758365857e-5)
    elseif z ≥ 52.70218954633605  # N =  8
        @horner(inv(z),1.0,0.08333333333333333,0.003472222222222222,-0.0026813271604938273,-0.00022947209362139917,0.0007840392217200666,6.972813758365857e-5,-0.0005921664373536939)
    elseif z ≥ 34.84031591198865  # N =  9
        @horner(inv(z),1.0,0.08333333333333333,0.003472222222222222,-0.0026813271604938273,-0.00022947209362139917,0.0007840392217200666,6.972813758365857e-5,-0.0005921664373536939,-5.171790908260592e-5)
    elseif z ≥ 25.3173982783047   # N = 10
        @horner(inv(z),1.0,0.08333333333333333,0.003472222222222222,-0.0026813271604938273,-0.00022947209362139917,0.0007840392217200666,6.972813758365857e-5,-0.0005921664373536939,-5.171790908260592e-5,0.0008394987206720873)
    elseif z ≥ 19.685015283078513 # N = 11
        @horner(inv(z),1.0,0.08333333333333333,0.003472222222222222,-0.0026813271604938273,-0.00022947209362139917,0.0007840392217200666,6.972813758365857e-5,-0.0005921664373536939,-5.171790908260592e-5,0.0008394987206720873,7.204895416020011e-5)
    elseif z ≥ 16.088669099569266 # N = 12
        @horner(inv(z),1.0,0.08333333333333333,0.003472222222222222,-0.0026813271604938273,-0.00022947209362139917,0.0007840392217200666,6.972813758365857e-5,-0.0005921664373536939,-5.171790908260592e-5,0.0008394987206720873,7.204895416020011e-5,-0.0019144384985654776)
    elseif z ≥ 13.655055978888104 # N = 13
        @horner(inv(z),1.0,0.08333333333333333,0.003472222222222222,-0.0026813271604938273,-0.00022947209362139917,0.0007840392217200666,6.972813758365857e-5,-0.0005921664373536939,-5.171790908260592e-5,0.0008394987206720873,7.204895416020011e-5,-0.0019144384985654776,-0.00016251626278391583)
    elseif z ≥ 11.93238782087875  # N = 14
        @horner(inv(z),1.0,0.08333333333333333,0.003472222222222222,-0.0026813271604938273,-0.00022947209362139917,0.0007840392217200666,6.972813758365857e-5,-0.0005921664373536939,-5.171790908260592e-5,0.0008394987206720873,7.204895416020011e-5,-0.0019144384985654776,-0.00016251626278391583,0.00640336283380807)
    elseif z ≥ 10.668852439197263 # N = 15
        @horner(inv(z),1.0,0.08333333333333333,0.003472222222222222,-0.0026813271604938273,-0.00022947209362139917,0.0007840392217200666,6.972813758365857e-5,-0.0005921664373536939,-5.171790908260592e-5,0.0008394987206720873,7.204895416020011e-5,-0.0019144384985654776,-0.00016251626278391583,0.00640336283380807,0.0005401647678926045)
    elseif z ≥ 9.715358216638403  # N = 16
        @horner(inv(z),1.0,0.08333333333333333,0.003472222222222222,-0.0026813271604938273,-0.00022947209362139917,0.0007840392217200666,6.972813758365857e-5,-0.0005921664373536939,-5.171790908260592e-5,0.0008394987206720873,7.204895416020011e-5,-0.0019144384985654776,-0.00016251626278391583,0.00640336283380807,0.0005401647678926045,-0.02952788094569912)
    elseif z ≥ 8.979120323411497  # N = 17
        @horner(inv(z),1.0,0.08333333333333333,0.003472222222222222,-0.0026813271604938273,-0.00022947209362139917,0.0007840392217200666,6.972813758365857e-5,-0.0005921664373536939,-5.171790908260592e-5,0.0008394987206720873,7.204895416020011e-5,-0.0019144384985654776,-0.00016251626278391583,0.00640336283380807,0.0005401647678926045,-0.02952788094569912,-0.002481743600264998)
    else
        gamma(z)*sqrt(z/2π)*exp(z)/z^z
    end
end

stirlingremainder(z::Number,N::Int) = (1+zeta(N))*gamma(N)/((2π)^(N+1)*z^N)/stirlingseries(z)

Aratio(n::Int,α::Float64,β::Float64) = exp((n/2+α+1/4)*log1p(-β/(n+α+β+1))+(n/2+β+1/4)*log1p(-α/(n+α+β+1))+(n/2+1/4)*log1p(α/(n+1))+(n/2+1/4)*log1p(β/(n+1)))
Aratio(n::Number,α::Number,β::Number) = (1+(α+1)/n)^(n+α+1/2)*(1+(β+1)/n)^(n+β+1/2)/(1+(α+β+1)/n)^(n+α+β+1/2)/(1+(zero(α)+zero(β))/n)^(n+1/2)

Cratio(n::Int,α::Float64,β::Float64) = exp((n+α+1/2)*log1p((α-β)/(2n+α+β+2))+(n+β+1/2)*log1p((β-α)/(2n+α+β+2))-log1p((α+β+2)/2n)/2)/sqrt(n)
Cratio(n::Number,α::Number,β::Number) = n^(-1/2)*(1+(α+1)/n)^(n+α+1/2)*(1+(β+1)/n)^(n+β+1/2)/(1+(α+β+2)/2n)^(2n+α+β+3/2)

Anαβ(n::Number,α::Number,β::Number) = 2^(α+β+1)/(2n+α+β+1)*exp(lgamma(n+α+1)-lgamma(n+α+β+1)+lgamma(n+β+1)-lgamma(n+1))
function Anαβ(n::Integer,α::Number,β::Number)
    if n==0
        2^(α+β+1)*beta(α+1,β+1)
    else
        val = Anαβ(0,α,β)
        for i=1:n
            val *= (i+α)*(i+β)/(i+α+β)/i*(2i+α+β-1)/(2i+α+β+1)
        end
        val
    end
end

function Anαβ(n::Integer,α::Float64,β::Float64)
    if n+min(α,β,α+β,0) ≥ 7.979120323411497
        2 ^ (α+β+1)/(2n+α+β+1)*stirlingseries(n+α+1)*Aratio(n,α,β)/stirlingseries(n+α+β+1)*stirlingseries(n+β+1)/stirlingseries(n+one(Float64))
    else
        (n+1)*(n+α+β+1)/(n+α+1)/(n+β+1)*Anαβ(n+1,α,β)*((2n+α+β+3)/(2n+α+β+1))
    end
end

"""
The Lambda function ``\\Lambda(z) = \\frac{\\Gamma(z+\\frac{1}{2})}{\\Gamma(z+1)}`` for the ratio of gamma functions.
"""
Λ(z::Number) = exp(lgamma(z+half(z))-lgamma(z+one(z)))
"""
For 64-bit floating-point arithmetic, the Lambda function uses the asymptotic series for ``\\tau`` in Appendix B of

I. Bogaert and B. Michiels and J. Fostier, 𝒪(1) computation of Legendre polynomials and Gauss–Legendre nodes and weights for parallel computing, *SIAM J. Sci. Comput.*, **34**:C83–C101, 2012.
"""
function Λ(x::Float64)
    if x > 9.84475
        xp = x+0.25
        @horner(inv(xp^2),1.0,-1.5625e-02,2.5634765625e-03,-1.2798309326171875e-03,1.343511044979095458984375e-03,-2.432896639220416545867919921875e-03,6.7542375336415716446936130523681640625e-03)/sqrt(xp)
    else
        (x+1.0)*Λ(x+1.0)/(x+0.5)
    end
end

"""
The Lambda function ``\\Lambda(z,λ₁,λ₂) = \\frac{\\Gamma(z+\\lambda_1)}{Γ(z+\\lambda_2)}`` for the ratio of gamma functions.
"""
Λ(z::Number,λ₁::Number,λ₂::Number) = exp(lgamma(z+λ₁)-lgamma(z+λ₂))
function Λ(x::Float64,λ₁::Float64,λ₂::Float64)
    if min(x+λ₁,x+λ₂) ≥ 8.979120323411497
        exp(λ₂-λ₁+(x-.5)*log1p((λ₁-λ₂)/(x+λ₂)))*(x+λ₁)^λ₁/(x+λ₂)^λ₂*stirlingseries(x+λ₁)/stirlingseries(x+λ₂)
    else
        (x+λ₂)/(x+λ₁)*Λ(x+1.,λ₁,λ₂)
    end
end

## TODO: deprecate when Lambert-W is supported in a mainstream repository such as SpecialFunctions.jl
"""
The principal branch of the Lambert-W function, defined by ``x = W_0(x) e^{W_0(x)}``, computed using Halley's method for ``x \\in [-e^{-1},\\infty)``.
"""
function lambertw(x::AbstractFloat)
    if x < -exp(-one(x))
        return throw(DomainError())
    elseif x == -exp(-one(x))
        return -one(x)
    elseif x < 0
        w0 = ℯ*x/(1+inv(inv(sqrt(2*ℯ*x+2))+inv(ℯ-1)-inv(sqrt(2))))
    else
        log1px = log1p(x)
        w0 = log1px*(1-log1p(log1px)/(2+log1px))
    end
    expw0 = exp(w0)
    w1 = w0 - (w0*expw0 - x)/((w0 + 1)*expw0 -
        (w0 + 2) * (w0*expw0 - x)/(2w0 + 2))
    while abs(w1/w0 - 1) > 2eps(typeof(x))
        w0 = w1
        expw0 = exp(w0)
        w1 = w0 - (w0*expw0 - x)/((w0 + 1)*expw0 -
            (w0 + 2) * (w0*expw0 - x)/(2w0 + 2))
    end
    return w1
end
lambertw(x::Real) = lambertw(float(x))

Cnλ(n::Integer,λ::Float64) = 2^λ/sqrtpi*Λ(n+λ)
Cnλ(n::Integer,λ::Number) = 2^λ/sqrt(oftype(λ,π))*Λ(n+λ)
function Cnλ(n::UnitRange{T},λ::Number) where T<:Integer
    ret = Vector{typeof(λ)}(undef, length(n))
    ret[1] = Cnλ(first(n),λ)
    for i=2:length(n)
        ret[i] = (n[i]+λ-half(λ))/(n[i]+λ)*ret[i-1]
    end
    ret
end

function Cnmλ(n::Integer,m::Integer,λ::Number)
    if m == 0
        Cnλ(n,λ)
    else
        (λ+m-1)/2/m*(m-λ)/(n+λ+m)*Cnmλ(n,m-1,λ)
    end
end

function Cnαβ(n::Integer,α::Number,β::Number)
    if n==0
        2^(α+β+1)*beta(α+1,β+1)/π
    else
        val = Cnαβ(0,α,β)
        for i=1:n
            val *= (i+α)*(i+β)/(i+(α+β+1)/2)/(i+(α+β)/2)
        end
        val
    end
end

function Cnαβ(n::Integer,α::Float64,β::Float64)
    if n+min(α,β) ≥ 7.979120323411497
        stirlingseries(n+α+1)/sqrtpi/stirlingseries(2n+α+β+2)*Cratio(n,α,β)*stirlingseries(n+β+1)
    else
        (n+(α+β+3)/2)/(n+β+1)*(n+(α+β+2)/2)/(n+α+1)*Cnαβ(n+1,α,β)
    end
end

function Cnmαβ(n::Integer,m::Integer,α::Number,β::Number)
    if m == 0
        Cnαβ(n,α,β)
    else
        Cnmαβ(n,m-1,α,β)/2(2n+α+β+m+1)
    end
end

function Cnmαβ(n::Integer,m::Integer,α::AbstractArray{T},β::AbstractArray{T}) where T<:Number
    shp = promote_shape(size(α),size(β))
    reshape([ Cnmαβ(n,m,α[i],β[i]) for i in eachindex(α,β) ], shp)
end

"""
Modified Chebyshev moments of the first kind:

```math
    \\int_{-1}^{+1} T_n(x) {\\rm\\,d}x.
```
"""
function chebyshevmoments1(::Type{T}, N::Int) where T
    μ = zeros(T, N)
    for i = 0:2:N-1
        @inbounds μ[i+1] = two(T)/T(1-i^2)
    end
    μ
end

"""
Modified Chebyshev moments of the first kind with respect to the Jacobi weight:

```math
    \\int_{-1}^{+1} T_n(x) (1-x)^\\alpha(1+x)^\\beta{\\rm\\,d}x.
```
"""
function chebyshevjacobimoments1(::Type{T}, N::Int, α, β) where T
    μ = zeros(T, N)
    N > 0 && (μ[1] = 2 .^ (α+β+1)*beta(α+1,β+1))
    if N > 1
        μ[2] = μ[1]*(β-α)/(α+β+2)
        for i=1:N-2
            @inbounds μ[i+2] = (2(β-α)*μ[i+1]-(α+β-i+2)*μ[i])/(α+β+i+2)
        end
    end
    μ
end

"""
Modified Chebyshev moments of the first kind with respect to the logarithmic weight:

```math
    \\int_{-1}^{+1} T_n(x) \\log\\left(\\frac{1-x}{2}\\right){\\rm\\,d}x.
```
"""
function chebyshevlogmoments1(::Type{T}, N::Int) where T
    μ = zeros(T, N)
    N > 0 && (μ[1] = -two(T))
    if N > 1
        μ[2] = -one(T)
        for i=1:N-2
            cst = isodd(i) ? T(4)/T(i^2-4) : T(4)/T(i^2-1)
            @inbounds μ[i+2] = ((i-2)*μ[i]+cst)/(i+2)
        end
    end
    μ
end

"""
Modified Chebyshev moments of the second kind:

```math
    \\int_{-1}^{+1} U_n(x) {\\rm\\,d}x.
```
"""
function chebyshevmoments2(::Type{T}, N::Int) where T
    μ = zeros(T, N)
    for i = 0:2:N-1
        @inbounds μ[i+1] = two(T)/T(i+1)
    end
    μ
end

"""
Modified Chebyshev moments of the second kind with respect to the Jacobi weight:

```math
    \\int_{-1}^{+1} U_n(x) (1-x)^\\alpha(1+x)^\\beta{\\rm\\,d}x.
```
"""
function chebyshevjacobimoments2(::Type{T}, N::Int, α, β) where T
    μ = zeros(T, N)
    N > 0 && (μ[1] = 2 .^ (α+β+1)*beta(α+1,β+1))
    if N > 1
        μ[2] = 2μ[1]*(β-α)/(α+β+2)
        for i=1:N-2
            @inbounds μ[i+2] = (2(β-α)*μ[i+1]-(α+β-i)*μ[i])/(α+β+i+2)
        end
    end
    μ
end

"""
Modified Chebyshev moments of the second kind with respect to the logarithmic weight:

```math
    \\int_{-1}^{+1} U_n(x) \\log\\left(\\frac{1-x}{2}\\right){\\rm\\,d}x.
```
"""
function chebyshevlogmoments2(::Type{T}, N::Int) where T
    μ = chebyshevlogmoments1(T, N)
    if N > 1
        μ[2] *= two(T)
        for i=1:N-2
            @inbounds μ[i+2] = 2μ[i+2] + μ[i]
        end
    end
    μ
end

function sphrand(::Type{T}, m::Int, n::Int) where T
    A = zeros(T, m, n)
    for i = 1:m
        A[i,1] = rand(T)
    end
    for j = 1:n÷2
        for i = 1:m-j
            A[i,2j] = rand(T)
            A[i,2j+1] = rand(T)
        end
    end
    A
end

function sphrandn(::Type{T}, m::Int, n::Int) where T
    A = zeros(T, m, n)
    for i = 1:m
        A[i,1] = randn(T)
    end
    for j = 1:n÷2
        for i = 1:m-j
            A[i,2j] = randn(T)
            A[i,2j+1] = randn(T)
        end
    end
    A
end

function sphones(::Type{T}, m::Int, n::Int) where T
    A = zeros(T, m, n)
    for i = 1:m
        A[i,1] = one(T)
    end
    for j = 1:n÷2
        for i = 1:m-j
            A[i,2j] = one(T)
            A[i,2j+1] = one(T)
        end
    end
    A
end

sphzeros(::Type{T}, m::Int, n::Int) where T = zeros(T, m, n)

"""
Pointwise evaluation of real orthonormal spherical harmonic:

```math
Y_\\ell^m(\\theta,\\varphi) = (-1)^{|m|}\\sqrt{(\\ell+\\frac{1}{2})\\frac{(\\ell-|m|)!}{(\\ell+|m|)!}} P_\\ell^{|m|}(\\cos\\theta) \\sqrt{\\frac{2-\\delta_{m,0}}{2\\pi}} \\left\\{\\begin{array}{ccc} \\cos m\\varphi & {\\rm for} & m \\ge 0,\\\\ \\sin(-m\\varphi) & {\\rm for} & m < 0.\\end{array}\\right.
```
"""
sphevaluate(θ, φ, L, M) = sphevaluatepi(θ/π, φ/π, L, M)

sphevaluatepi(θ::Number, φ::Number, L::Integer, M::Integer) = sphevaluatepi(θ, L, M)*sphevaluatepi(φ, M)

function sphevaluatepi(θ::Number, L::Integer, M::Integer)
    ret = one(θ)/sqrt(two(θ))
    if M < 0 M = -M end
    c, s = cospi(θ), sinpi(θ)
    for m = 1:M
        ret *= sqrt((m+half(θ))/m)*s
    end
    tc = two(c)*c

    if L == M
        return ret
    elseif L == M+1
        return sqrt(two(θ)*M+3)*c*ret
    else
        temp = ret
        ret *= sqrt(two(θ)*M+3)*c
        for l = M+1:L-1
            ret, temp = (sqrt(l+half(θ))*tc*ret - sqrt((l-M)*(l+M)/(l-half(θ)))*temp)/sqrt((l-M+1)*(l+M+1)/(l+3half(θ))), ret
        end
        return ret
    end
end

sphevaluatepi(φ::Number, M::Integer) = sqrt((two(φ)-δ(M, 0))/(two(φ)*π))*(M ≥ 0 ? cospi(M*φ) : sinpi(-M*φ))

function sphvrand(::Type{T}, m::Int, n::Int) where T
    A = zeros(T, m, n)
    for i = 1:m-1
        A[i,1] = rand(T)
    end
    for j = 1:n÷2
        for i = 1:m-j+1
            A[i,2j] = rand(T)
            A[i,2j+1] = rand(T)
        end
    end
    A
end

function sphvrandn(::Type{T}, m::Int, n::Int) where T
    A = zeros(T, m, n)
    for i = 1:m-1
        A[i,1] = randn(T)
    end
    for j = 1:n÷2
        for i = 1:m-j+1
            A[i,2j] = randn(T)
            A[i,2j+1] = randn(T)
        end
    end
    A
end

function sphvones(::Type{T}, m::Int, n::Int) where T
    A = zeros(T, m, n)
    for i = 1:m-1
        A[i,1] = one(T)
    end
    for j = 1:n÷2
        for i = 1:m-j+1
            A[i,2j] = one(T)
            A[i,2j+1] = one(T)
        end
    end
    A
end

sphvzeros(::Type{T}, m::Int, n::Int) where T = sphzeros(T, m, n)

function diskrand(::Type{T}, m::Int, n::Int) where T
    A = zeros(T, m, n)
    for i = 1:m
        A[i,1] = rand(T)
    end
    for j = 1:n÷2
        for i = 1:m-(j+1)÷2
            A[i,2j] = rand(T)
            A[i,2j+1] = rand(T)
        end
    end
    A
end

function diskrandn(::Type{T}, m::Int, n::Int) where T
    A = zeros(T, m, n)
    for i = 1:m
        A[i,1] = randn(T)
    end
    for j = 1:n÷2
        for i = 1:m-(j+1)÷2
            A[i,2j] = randn(T)
            A[i,2j+1] = randn(T)
        end
    end
    A
end

function diskones(::Type{T}, m::Int, n::Int) where T
    A = zeros(T, m, n)
    for i = 1:m
        A[i,1] = one(T)
    end
    for j = 1:n÷2
        for i = 1:m-(j+1)÷2
            A[i,2j] = one(T)
            A[i,2j+1] = one(T)
        end
    end
    A
end

diskzeros(::Type{T}, m::Int, n::Int) where T = zeros(T, m, n)

function trirand(::Type{T}, m::Int, n::Int) where T
    A = zeros(T, m, n)
    for j = 1:n
        for i = 1:m+1-j
            A[i,j] = rand(T)
        end
    end
    A
end

function trirandn(::Type{T}, m::Int, n::Int) where T
    A = zeros(T, m, n)
    for j = 1:n
        for i = 1:m+1-j
            A[i,j] = randn(T)
        end
    end
    A
end

function triones(::Type{T}, m::Int, n::Int) where T
    A = zeros(T, m, n)
    for j = 1:n
        for i = 1:m+1-j
            A[i,j] = one(T)
        end
    end
    A
end

trizeros(::Type{T}, m::Int, n::Int) where T = zeros(T, m, n)

"""
Pointwise evaluation of triangular harmonic:

```math
\\tilde{P}_{\\ell,m}^{(\\alpha,\\beta,\\gamma)}(x,y).
```
"""
trievaluate(x, y, L, M, α, β, γ) = trievaluate(x, L, M, α, β, γ)*trievaluate(x, y, M, β, γ)

function trievaluate(x::Number, L::Integer, M::Integer, α::Number, β::Number, γ::Number)

end

function trievaluate(x::Number, y::Number, M::Integer, β::Number, γ::Number)

end

function tetrand(::Type{T}, l::Int, m::Int, n::Int) where T
    A = zeros(T, l, m, n)
    for k = 1:n
        for j = 1:m+1-k
            for i = 1:l+2-k-j
                A[i,j,k] = rand(T)
            end
        end
    end
    A
end

function tetrandn(::Type{T}, l::Int, m::Int, n::Int) where T
    A = zeros(T, l, m, n)
    for k = 1:n
        for j = 1:m+1-k
            for i = 1:l+2-k-j
                A[i,j,k] = randn(T)
            end
        end
    end
    A
end

function tetones(::Type{T}, l::Int, m::Int, n::Int) where T
    A = zeros(T, l, m, n)
    for k = 1:n
        for j = 1:m+1-k
            for i = 1:l+2-k-j
                A[i,j,k] = one(T)
            end
        end
    end
    A
end

tetzeros(::Type{T}, l::Int, m::Int, n::Int) where T = zeros(T, l, m, n)


text:
using Test
using Forecast

import DataFrames: DataFrame

@testset "p" begin

    x = [1,2,3,4,5]
    x0 = [0,1,2,3,4]
    @test p(d(x)) == x0
    @test p(d(x),[x[1]]) == x
    @test p(d(hcat(x,x,x))) == hcat(x0,x0,x0)
    @test p(d(hcat(x,x,x)), [ x[1] x[1] x[1] ] ) == hcat(x,x,x)

    x = rand(100)
    @test p(d(x),[ x[1] ]) ≈ x
    @test p(d(hcat(x,x,x)), [ x[1] x[1] x[1] ] ) ≈ hcat(x,x,x)

    x = [repeat(1:10,10) repeat(1:10,10)]
    dx = d(x,3)
    x0 = [[1,1,0] [1,1,0]]
    @test p(dx, x0) == x

    la = 11
    x = [repeat(1:la,14) repeat(1:la,14)
         1 1]
    dx = d(x,1,la)
    x0 = reshape(repeat(1:la,inner=2),1,2,la)
    @test p(dx, x0) == x

    x = repeat(1:2,30)
    dx = d(x,2,2)
    x0 = zeros(2,1,2) # order 2, 1 variable, lag 1
    x0[1,:,:] = collect(1:2)
    @test p(dx,x0) ≈ x

    tx = hcat(co2().co2, co2().co2, co2().co2)
    for col in eachcol(tx) replace!(col, missing => 0.0) end
    @test p(d(tx), [ 333.38 333.38 333.38 ] ) ≈ tx
    @test p(d(tx,2), [ 333.38 333.38 333.38;
                        -0.27  -0.27  -0.27]    ) ≈ tx

    # calculation of π
    x = 0:0.001:1
    y = sqrt.(1 .- x.^2);
    isapprox(4*p(y)[end]/1000 , π, atol = 0.01)

end


text:
module MinAtar

using PyCall
using RecipesBase
using Colors

export
    MinAtarEnv,
    act!,
    reset!,
    get_state,
    state_shape,
    num_actions,
    gamename,
    minimal_action_set

"""
    MinAtarEnv(gamename, sticky_action_prob, difficulty_ramping, random_seed)

A wrapper struct for the python MinAtar project.

"""
mutable struct MinAtarEnv
    gamename::String
    pygame::PyObject  # the pythong game
    pyact::PyObject  # the python env.step function
    pyreset::PyObject # the python env.reset function
    pystate::PyObject # the state array object referenced by the PyArray state.o
    sticky_action_prob::Float64
    difficulty_ramping::Bool
    reward::Float64
    total_reward::Float64
    done::Bool
    colors::Array{RGB{Float64},1}
    function MinAtarEnv(gamename,
                        sticky_action_prob = 0.1,
                        difficulty_ramping = true,
                        random_seed = nothing)

        minitar_env = pyimport("minatar")
        pygame = minitar_env.environment.Environment(gamename, sticky_action_prob, difficulty_ramping, random_seed)
        sns = pyimport("seaborn")
        new(gamename, pygame, pygame.act, pygame.reset,
            pygame.state, sticky_action_prob, difficulty_ramping,
            0.0, 0.0, false,
            [RGB(sns.color_palette("cubehelix", pygame.n_channels)[i]...) for i in 1:pygame.n_channels])
    end
end

act!(env::MinAtarEnv, a) = env.pyact(a)
reset!(env::MinAtarEnv) = env.pyreset()
get_state(env::MinAtarEnv) = env.pystate()
state_shape(env::MinAtarEnv) = env.pygame.state_shape()
num_actions(env::MinAtarEnv) = 6
gamename(env::MinAtarEnv) = env.gamename
minimal_action_state(env::MinAtarEnv) = env.pygame.minimal_action_set()

@recipe function f(env::MinAtarEnv)
    ticks := nothing
    foreground_color_border := nothing
    grid := false
    legend := false
    aspect_ratio := 1

    screen = fill(Colors.RGB(0.0, 0.0, 0.0), 10, 10)
    state = env.pystate()
    if eltype(state) isa Float64
        for i in 1:size(state)[end]
            screen[findall((x)->x==1.0, state[:,:, i])] .= env.colors[i]
        end
    elseif eltype(state) isa Bool
        for i in 1:size(state)[end]
            screen[findall((x)->x, state[:,:, i])] .= env.colors[i]
        end
    end
    screen
end

end # module


text:
using DataFrames, CSV, FileIO, Dates
using BenchmarkTools
using Agtor

import Agtor: run_model

import OrderedCollections: LittleDict

# Start month/day of growing season
# Represents when water allocations are announced.
const gs_start = (5, 1)

function setup_zone(data_dir::String="test/data/")
    zone_spec_dir::String = "$(data_dir)zones"
    zone_params::Dict{Symbol, Dict} = load_spec(zone_spec_dir)

    return create(zone_params[:Zone_1])
end

"""
Allocation callback for example.
"""
function allocation_callback!(zone, dt_i)

    # Resetting allocations for example run
    if monthday(dt_i) == gs_start
        # Example annual water allocations
        allocs = (surface_water=150.0, groundwater=40.0)

        update_available_water!(zone, allocs)
    end
end

data_dir = "test/data/"
rigid_zone = setup_zone(data_dir)
scen_data = joinpath(data_dir, "scenarios", "sampled_params.csv")
samples = DataFrame!(CSV.File(scen_data))

farmer = RigidManager("Rigid", 0.04)

all_results = []
tmp_zone = deepcopy(rigid_zone)
for (row_id, r) in enumerate(eachrow(samples))
    update_model!(tmp_zone, r)

    zone_results, field_results = run_model(farmer, tmp_zone; post=allocation_callback!)
    # push!(all_results, (zone_results, field_results))

    save_results!("rigid_run.jld2", row_id, (zone_results, field_results))
end

# save_results!("rigid_run.jld2", all_results)

# Loading and displaying results
saved_results = load("rigid_run.jld2")
@info saved_results

max_stars_repo_name: ConnectedSystems/agtor-jl

max_issues_repo_name: ConnectedSystems/agtor-jl

max_forks_repo_name: ConnectedSystems/agtor-jl


text:
module Interface

using Test, Random, MLJBase,
    Tables, CategoricalArrays,
    OrderedCollections,
    TypedTables, MLJModelInterface,
    StableRNGs

using ..Models

rng = StableRNG(1551234)

include("data_utils.jl")

include("model_api.jl")

end

true


text:
module MySymPy

include("walk.jl")

end # module


text:
using Test

@testset "Rcl" begin

include("unit/base.jl")
include("unit/node_base.jl")

end


text:
struct DisplayedVisualizer
    core::CoreVisualizer
end

DisplayedVisualizer(vis::Visualizer) = DisplayedVisualizer(vis.core)

url(c::DisplayedVisualizer) = url(c.core)

"""
Render a MeshCat visualizer inline in Jupyter, Juno, or VSCode.

If this is the last command in a Jupyter notebook cell, then the
visualizer should show up automatically in the corresponding output
cell.

If this is run from the Juno console, then the visualizer should show
up in the Juno plot pane. Likewise if this is run from VSCode with
the julia-vscode extension, then the visualizer should show up in the
Julia Plots pane.
"""
render(vis::Visualizer) = render(vis.core)
render(core::CoreVisualizer) = DisplayedVisualizer(core)

@deprecate IJuliaCell(v::Visualizer) render(v)

function Base.show(io::IO,
        ::Union{MIME"text/html", MIME"juliavscode/html"},
        frame::DisplayedVisualizer)
    wait_for_server(frame.core)
    print(io, """

""")
end

function Base.show(io::IO,
        ::MIME"application/prs.juno.plotpane+html",
        d::DisplayedVisualizer)
    wait_for_server(d.core)
    print(io, """

    """)
end

function _create_command(data::Vector{UInt8})
    return """
fetch("data:application/octet-binary;base64,$(base64encode(data))")
    .then(res => res.arrayBuffer())
    .then(buffer => viewer.handle_command_bytearray(new Uint8Array(buffer)));
    """
end

"""
Extract a single HTML document containing the entire MeshCat scene,
including all geometries, properties, and animations, as well as all
required javascript assets. The resulting HTML document should render
correctly after you've exited Julia, and even if you have no internet
access.
"""
static_html(vis::Visualizer) = static_html(vis.core)

function static_html(core::CoreVisualizer)
    viewer_commands = String[]

    foreach(core.tree) do node
        if node.object !== nothing
            push!(viewer_commands, _create_command(node.object));
        end
        if node.transform !== nothing
            push!(viewer_commands, _create_command(node.transform));
        end
        for data in values(node.properties)
            push!(viewer_commands, _create_command(data));
        end
        if node.animation !== nothing
            push!(viewer_commands, _create_command(node.animation))
        end
    end

    return """
        <!DOCTYPE html>

             <meta charset=utf-8> <title>MeshCat</title>

    """
end

struct StaticVisualizer
    core::CoreVisualizer
end

"""
Render a static version of the visualizer, suitable for embedding
and offline use. The embedded visualizer includes all geometries,
properties, and animations which have been added to the scene, baked
into a single HTML document. This document also includes the full
compressed MeshCat javascript source files, so it should render
properly even after you've exited Julia and even if you have no
internet access.

To get access to the raw static HTML representation, see `static_html`
"""
render_static(vis::Visualizer, args...; kw...) = render_static(vis.core)

render_static(core::CoreVisualizer) = StaticVisualizer(core)

_srcdoc_escape(x) = replace(replace(x, '&' => "&"), '\"' => """)

function static_iframe_wrapper(html)
    id = uuid1()
    return """

    """
end

function Base.show(io::IO,
        ::Union{MIME"text/html", MIME"juliavscode/html"},
        static_vis::StaticVisualizer)
    print(io, static_iframe_wrapper(static_html(static_vis.core)))
end


text:
using Sundials, Base.Test

callback_f = function (du, u, p, t)
  du[1] = u[2]
  du[2] = -9.81
end

condtion= function (u,t,integrator) # Event when event_f(u,t,k) == 0
  u[1]
end

affect! = nothing
affect_neg! = function (integrator)
  integrator.u[2] = -integrator.u[2]
end

callback = ContinuousCallback(condtion,affect!,affect_neg!)

u0 = [50.0,0.0]
tspan = (0.0,15.0)
prob = ODEProblem(callback_f,u0,tspan)

sol = solve(prob,CVODE_Adams(),callback=callback)
@test sol(4.0)[1] > 0
sol = solve(prob,CVODE_BDF(),callback=callback)
@test sol(4.0)[1] > 0

u0 = [1.,0.]
function fun2(du, u, p, t)
   du[2] = -u[1]
   du[1] = u[2]
end
tspan = (0.0,10.0)
prob = ODEProblem(fun2,u0,tspan)

function condition2(u,t,integrator)
  get_du(integrator)[1]>0
end
affect2!(integrator) = terminate!(integrator)
cb = DiscreteCallback(condition2,affect2!)
sol = solve(prob,CVODE_BDF(),callback=cb)
@test sol.t[end] < 3.5

condition3(u,t,integrator) = u[2]
affect3!(integrator) = terminate!(integrator)
cb = ContinuousCallback(condition3,affect3!)
sol = solve(prob,CVODE_Adams(),callback=cb,abstol=1e-12,reltol=1e-12)
@test sol.t[end] ≈ pi


text:
function _handle_base_kv(pwf_data::Dict, bus::Dict, dict_dgbt)
    group_identifier = bus["BASE VOLTAGE GROUP"]
    if haskey(pwf_data, "DGBT")
        if haskey(dict_dgbt, group_identifier)
            group = dict_dgbt[group_identifier]
            return group["VOLTAGE"]
        elseif length(pwf_data["DGBT"]) == 1
            return pwf_data["DGBT"]["1"]["VOLTAGE"]
        end
    end
    return 1.0
end

function _handle_vmin(pwf_data::Dict, bus::Dict, dict_dglt)
    group_identifier = bus["VOLTAGE LIMIT GROUP"]
    if haskey(pwf_data, "DGLT")
        if haskey(dict_dglt, group_identifier)
            group = dict_dglt[group_identifier]
            return group["LOWER BOUND"]
        elseif length(pwf_data["DGLT"]) == 1
            return pwf_data["DGLT"]["1"]["LOWER BOUND"]
        end
    end
    return 0.9
end

function _handle_vmax(pwf_data::Dict, bus::Dict, dict_dglt)
    group_identifier = bus["VOLTAGE LIMIT GROUP"]
    if haskey(pwf_data, "DGLT")
        if haskey(dict_dglt, group_identifier)
            group = dict_dglt[group_identifier]
            return group["UPPER BOUND"]
        elseif length(pwf_data["DGLT"]) == 1
            return pwf_data["DGLT"]["1"]["UPPER BOUND"]
        end
    end
    return 1.1
end

function _handle_bus_type(bus::Dict)
    bus_type = bus["TYPE"]
    dict_bus_type = Dict(
        0 => 1,
        3 => 1, # PQ
        1 => 2, # PV
        2 => 3 # Referência
    )
    if bus["STATUS"] == 'L'
        return dict_bus_type[bus_type]
    elseif bus["STATUS"] == 'D'
        return 4
    end
end

function _create_dict_dglt(data::Dict)
    dict_dglt = Dict{String, Any}()

    for (k,v) in data
        sub_data = Dict{String, Any}()
        sub_data["LOWER BOUND"] = v["LOWER BOUND"]
        sub_data["UPPER BOUND"] = v["UPPER BOUND"]
        sub_data["LOWER EMERGENCY BOUND"] = v["LOWER EMERGENCY BOUND"]
        sub_data["UPPER EMERGENCY BOUND"] = v["UPPER EMERGENCY BOUND"]
        dict_dglt[v["GROUP"]] = sub_data
    end
    return dict_dglt
end

function _create_dict_dgbt(data::Dict)
    dict_dgbt = Dict{String, Any}()

    for (k,v) in data
        sub_data = Dict{String, Any}("VOLTAGE" => v["VOLTAGE"])
        dict_dgbt[v["GROUP"]] = sub_data
    end
    return dict_dgbt
end

function _pwf2pm_bus!(pm_data::Dict, pwf_data::Dict, bus::Dict, dict_dgbt, dict_dglt; add_control_data::Bool=false)
    sub_data = Dict{String,Any}()

    sub_data["bus_i"] = bus["NUMBER"]
    sub_data["bus_type"] = _handle_bus_type(bus)
    sub_data["area"] = pop!(bus, "AREA")
    sub_data["vm"] = bus["VOLTAGE"]
    sub_data["va"] = pop!(bus, "ANGLE")
    sub_data["zone"] = 1
    sub_data["name"] = pop!(bus, "NAME")

    sub_data["source_id"] = ["bus", "$(bus["NUMBER"])"]
    sub_data["index"] = bus["NUMBER"]

    sub_data["base_kv"] = _handle_base_kv(pwf_data, bus, dict_dgbt)
    sub_data["vmin"] = _handle_vmin(pwf_data, bus, dict_dglt)
    sub_data["vmax"] = _handle_vmax(pwf_data, bus, dict_dglt)

    if add_control_data
        sub_data["control_data"] = Dict{String,Any}()
        sub_data["control_data"]["voltage_controlled_bus"] = bus["CONTROLLED BUS"]
        sub_data["control_data"]["vmmin"] = sub_data["vmin"]
        sub_data["control_data"]["vmmax"] = sub_data["vmax"]
    end

    idx = string(sub_data["index"])
    pm_data["bus"][idx] = sub_data
end

function _pwf2pm_bus!(pm_data::Dict, pwf_data::Dict; add_control_data::Bool=false)

    dict_dglt = haskey(pwf_data, "DGLT") ? _create_dict_dglt(pwf_data["DGLT"]) : nothing
    isa(dict_dglt, Dict) && length(dict_dglt) == 1 ? @warn("Only one limit voltage group definded, each bus will be considered as part of the group $(pwf_data["DGLT"]["1"]["GROUP"]), regardless of its defined group") : nothing
    dict_dgbt = haskey(pwf_data, "DGBT") ? _create_dict_dgbt(pwf_data["DGBT"]) : nothing
    isa(dict_dgbt, Dict) && length(dict_dgbt) == 1 ? @warn("Only one base voltage group definded, each bus will be considered as part of the group $(pwf_data["DGBT"]["1"]["GROUP"]), regardless of its defined group") : nothing

    pm_data["bus"] = Dict{String, Any}()
    if haskey(pwf_data, "DBAR")
        for (i,bus) in pwf_data["DBAR"]
            _pwf2pm_bus!(pm_data, pwf_data, bus, dict_dgbt, dict_dglt, add_control_data = add_control_data)
        end
    end
end


text:
using Estimation
using Gadfly
using DSP

#
# Create a test signal similar to Li & Stoica 96
#

# Li & Stoica fig4
spectral_lines_f = [0.04, 0.05, 0.08, 0.09, 0.195, 0.215]
spectral_lines_p = vec(repmat([pi/4], length(spectral_lines_f)))
spectral_lines_a = [10, 10, 10, 2, 4, 6]

# Simulation properties
N = 128
noise_mean = 0
noise_vari = 1
num_trials = 100
padd_times = 16

function create_example(N, spectral_lines_f, spectral_lines_p, spectral_lines_a)
    clean = vec(zeros(N, 1))
    n = 1.0:N
    for i in 1:length(spectral_lines_f)
        clean = clean + spectral_lines_a[i] * sin(2 * pi * spectral_lines_f[i] * n + spectral_lines_p[i])
    end
    return clean
end

clean_signal = create_example(N, spectral_lines_f, spectral_lines_p, spectral_lines_a)
noise = randn(N)

z = clean_signal + noise

#
# Create a plot of the true spectrum
#

# Plot true spectrum
yline = zeros(length(0.0:0.0001:0.5))
for i in 1:length(spectral_lines_f); yline[round(spectral_lines_f[i]*10000)] = spectral_lines_a[i]; end

# Save plot as a layer so you can overlay them all later
true_plot = layer(
    x=0.0:0.0001:0.5,
    y=yline,
    Geom.line(),
    Theme(default_color=color("black")))

# Plot the true spectrum
p0 = plot(true_plot,
    Guide.title("True spectrum   N = $N"),
    Guide.xlabel("Frequency (Hz)"),
    Guide.ylabel("Amplitude"),
    Scale.y_continuous(minvalue=0, maxvalue=12))

#
# Plot FFT of signal
#

p = periodogram(z, nfft=length(z)*padd_times)
fft_plot = layer(
        x=freq(p),
        y=sqrt(power(p)/(length(z)/2)),
        Theme(default_color=color("blue")),
        Geom.line);
p1 = plot(fft_plot, true_plot,
    Guide.title("FFT spectrum"),
    Guide.xlabel("Frequency (Hz)"),
    Guide.ylabel("Amplitude"),
    Scale.y_continuous(minvalue=0, maxvalue=12));

#
# Plot Kaiser window FFT
#

window = kaiser(N, 4/pi)
window_power = mean(window.^2)

p = periodogram(z, nfft=length(z)*padd_times, window=window)

kaiser_plot = layer(
    x=freq(p),
    y=sqrt(power(p)/(length(z)/2)/sqrt(window_power)),
    Theme(default_color=color("purple")),
    Geom.line)

p2 = plot(kaiser_plot, true_plot,
    Guide.title("FFT spectrum with Kaiser window"),
    Guide.xlabel("Frequency (Hz)"),
    Guide.ylabel("Amplitude"),
    Scale.y_continuous(minvalue=0, maxvalue=12))

#
# Plot APES method
#

f = 0.00:0.001:0.5
a = apes(z, f)

apes_plot = layer(
    x=f,
    y=abs(a),
Theme(default_color=color("red")),
    Geom.line)

p5 = plot(apes_plot,
    true_plot,
Guide.title("APES spectrum with M = $(int(N/2))"),
    Guide.xlabel("Frequency (Hz)"),
    Guide.ylabel("Amplitude"),
    Scale.y_continuous(minvalue=0, maxvalue=12))

#
# Connect plots and save
#

set_default_plot_size(24cm, 12cm)
Pt = Gadfly.gridstack([p0 p1; p2 p5])
# draw(PNG("Li-Stoica-Fig1.png", 28cm, 28cm), Pt)
#=draw(PDF("Li-Stoica-Fig4-Presentation.pdf", 26cm, 16cm), Pt)=#
draw(PNG("APES-example.png", 26cm, 16cm), Pt)
#=Pt=#


text:
module TestLalePipeOptimizerAG

using Random
using Test
using Lale
using Statistics
using DataFrames: DataFrame
using AutoMLPipeline.Utils

const IRIS = getiris()
const X    = IRIS[:,1:3] |> DataFrame
const XC   = IRIS[:,1:4] |> DataFrame
const YC   = IRIS[:,5] |> Vector
const Y    = IRIS[:,4] |> Vector

function pipeline_test()

   # lale ops
   pca     = LaleOp("PCA","autogen")
   rb      = LaleOp("RobustScaler","autogen")
   rfr     = LaleOp("RandomForestRegressor")
   rfc     = laleoperator("RandomForestClassifier")
   svr     = LaleOp("SVR")
   noop    = laleoperator("NoOp","lale")

   @test_throws ArgumentError LaleOp("PCAN")
   @test_throws ArgumentError LaleOp("RandomForestClassifers")

   # amlp ops
   ohe  = OneHotEncoder()
   catf = CatFeatureSelector()
   numf = NumFeatureSelector()

   pcatrained            = fit(pca,X)
   res                   = transform(pcatrained,X)
   @assert size(res,2)   == 3

   amlpipe = @pipeline  (rb + pca + noop) |> (rfr * svr)
   amlpred = fit_transform!(amlpipe,X,Y)
   amlprmse=score(:rmse,amlpred,Y)

   # regression
   lalepipe   = (rb + pca + noop) >>  (rfr | svr )
   lale_hopt  = LalePipeOptimizer(lalepipe,max_evals  = 3,cv = 3)
   lalepred   = fit_transform!(lale_hopt,X,Y)
   lalermse   = score(:rmse,lalepred,Y)

   @test amlprmse < lalermse

   laletrained = fit(lale_hopt,X,Y)
   lalepred    = transform(laletrained,X)
   lalermse    = score(:rmse,lalepred,Y)

   @test amlprmse < lalermse

   laletrained = fit(lale_hopt,X,Y)
   lalepred    = predict(laletrained,X)
   lalermse    = score(:rmse,lalepred,Y)

   @test amlprmse < lalermse

   # classification
   lalepipe  = (rb+pca+noop)  >>  rfc
   lale_hopt = LalePipeOptimizer(lalepipe,max_evals = 3,cv = 3)
   lalepred  = fit_transform!(lale_hopt,XC,YC)
   laleacc   = score(:accuracy,lalepred,YC)

   amlpipe = @pipeline  (pca + noop) |> (rfc)
   amlpred = fit_transform!(amlpipe,XC,YC)
   amlpacc = score(:accuracy,amlpred,YC)

   @test abs(amlpacc - laleacc) > 0.0

   lalepipe    = (pca  >>  rfc)
   lale_hopt   = LalePipeOptimizer(lalepipe,max_evals = 3,cv = 3)
   laletrained = fit(lale_hopt,XC,YC)
   lalepred    = predict(laletrained,XC)
   laleacc     = score(:accuracy,lalepred,YC)

   @test abs(amlpacc - laleacc) > 0.0

   plr = @pipeline (catf |> ohe) + (numf |> rb |> pca) |> rfr
   plc = @pipeline (catf |> ohe) + (numf |> rb |> pca) |> rfc
   @test crossvalidate(plr,X,Y,"mean_absolute_error",3,false).mean < 0.3
   @test crossvalidate(plc,XC,YC,"accuracy_score",3,false).mean > 0.8

end
@testset "Lale Pipeline Optimizer" begin
   pipeline_test()
end

end


text:
model = IterationControl.SquareRooter(4.0)
IterationControl.train!(model, 1)
@test model.root ≈ 2.5
@test IterationControl.loss(model) ≈ 25/4 - 4
IterationControl.train!(model, 100)
@test IterationControl.loss(model) ≈ 0
@test IterationControl.training_losses(model)[1:2] ≈
    abs.([41/20 - 5/2, 3281/1640 - 41/20])

true


text:
import Libdl

@info "This package assumes you've already installed Intel RealSense SDK 2.0 on Windows."
const WINDOWS_SDK_DIR = joinpath(homedir(), "..", "..", "Program Files (x86)") |> normpath
const WINDOWS_DEFAULT_DIR = Sys.WORD_SIZE == 64 ? joinpath(WINDOWS_SDK_DIR, "Intel RealSense SDK 2.0", "bin", "x64") : joinpath(WINDOWS_SDK_DIR, "Intel RealSense SDK 2.0", "bin", "x86")
# possible shared library locations, you could also add yours if you installed the SDK in other place
locations = [WINDOWS_DEFAULT_DIR]
# try to find the lib
lib = Libdl.find_library(["realsense2"], locations)
if lib == ""
    error("Couldn't find realsense2.dll in the default directory $WINDOWS_DEFAULT_DIR. Please install Intel RealSense SDK 2.0 and then run Pkg.build(\"LibRealSense\").")
else
    const librealsense = Libdl.dlpath(lib) |> x->replace(x, "\\"=>"\\\\")
end

open(joinpath(@__DIR__, "deps.jl"), "w") do f
    println(f, """
    import Libdl

    const LIBREALSENSE_VERSION = v"2.15.0"
    const librealsense2 = "$librealsense"

    function check_deps()
        global librealsense2
        if !isfile(librealsense2)
            error("\$(librealsense2) does not exist, Please re-run Pkg.build(\\"LibRealSense\\"), and restart Julia.")
        end

        hdl = Libdl.dlopen_e(librealsense2)
        if hdl == C_NULL
            error("\$(librealsense2) cannot be opened, Please re-run Pkg.build(\\"LibRealSense\\"), and restart Julia.")
        else
            get_api_version = Libdl.dlsym_e(hdl, :rs2_get_api_version)
            @assert get_api_version != C_NULL "Could not find `rs2_get_api_version` within \$(librealsense2)."
            err = Ref{Ptr{Cvoid}}(0)
            version = ccall(get_api_version, Cint, (Ptr{Ptr{Cvoid}},), err)
            @assert err[] == C_NULL "Something happened when calling `rs2_get_api_version`."
            major = version ÷ 10000
            minor = (version - 10000*major) ÷ 100
            patch = (version - 10000*major - 100*minor)
            if major != LIBREALSENSE_VERSION.major && minor != LIBREALSENSE_VERSION.minor
                error(\"\"\"API version mismatch:
                      The system librealsense's API version is v\$(major).\$(minor).\$(patch), but the one used in LibRealSense.jl is \$(LIBREALSENSE_VERSION).
                      Please either upgrade Intel RealSense SDK or LibRealSense.jl.\"\"\")
            end
        end
    end
    """)
end