text:
using Compat
using Markdown

#Contains report global properties
mutable struct Report <: AbstractDisplay
  cwd::AbstractString
  basename::AbstractString
  formatdict::Dict{Symbol,Any}
  pending_code::AbstractString
  cur_result::AbstractString
  rich_output::AbstractString
  fignum::Int
  figures::Array{AbstractString}
  term_state::Symbol
  cur_chunk
  mimetypes::Array{AbstractString}
  first_plot::Bool
  header_script::String
  throw_errors::Bool
end

function Report(cwd, basename, formatdict, mimetypes, throw_errors)
    Report(cwd, basename, formatdict, "", "", "", 1, AbstractString[], :text, nothing,
        mimetypes, true, "", throw_errors)
end

#Default mimetypes in order, can be overridden for some inside `run method` formats
const default_mime_types = ["image/svg+xml", "image/png", "text/html", "text/plain"]
#const default_mime_types = ["image/png", "image/svg+xml", "text/html", "text/plain"]
#From IJulia as a reminder
#const supported_mime_types = [ "text/html", "text/latex", "image/svg+xml", "image/png", "image/jpeg", "text/plain", "text/markdown" ]

const mimetype_ext =
    Dict(".png" => "image/png",
         ".jpg" => "image/jpeg",
         ".jpeg" => "image/jpeg",
         ".svg" => "image/svg+xml",
         ".pdf" => "application/pdf")

function Base.display(report::Report, data)
    #Set preferred mimetypes for report based on format
    fig_ext = report.cur_chunk.options[:fig_ext]
    for m in unique([mimetype_ext[fig_ext] ; report.mimetypes])
        if showable(m, data)
            try
                if !istextmime(m)
                    Compat.invokelatest(display, report, m, data)
                elseif report.cur_chunk.options[:term]
                    Compat.invokelatest(display, report, "text/plain", data)
                else
                    Compat.invokelatest(display, report, m, data)
                end
            catch e
                throw(e)
                @warn("Failed to display data in \"$m\" format")
                continue
            end
            break
        end
    end
end

function Base.display(report::Report, m::MIME"image/png", data)
    figname = add_figure(report, data, m, ".png")
end

function Base.display(report::Report, m::MIME"image/svg+xml", data)
    figname = add_figure(report, data, m, ".svg")
end

function Base.display(report::Report, m::MIME"application/pdf", data)
    figname = add_figure(report, data, m, ".pdf")
end

#Text is written to stdout, called from "term" mode chunks
function Base.display(report::Report, m::MIME"text/plain", data)
    io = PipeBuffer()
    show(IOContext(io, :limit => true), m, data)
    flush(io)
    s = read(io, String)
    close(io)
    println(s)
end

function Base.display(report::Report, m::MIME"text/plain", data::Exception)
    println("Error: " * sprint(showerror, data))
end

function Base.display(report::Report, m::MIME"text/html", data::Exception)
    report.rich_output = sprint(show, m, data)
end

function Base.show(io, m::MIME"text/html", data::Exception)
    println(io ,"
")
    println(io, Markdown.htmlesc("ERROR: " * sprint(showerror, data)))
    println(io ,"
")
end

#Catch "rich_output"
function Base.display(report::Report, m::MIME"text/html", data)
    s = repr(m, data)
    report.rich_output *= "\n" * s
end

#Catch "rich_output"
function Base.display(report::Report, m::MIME"text/markdown", data)
    s = repr(m, data)
    report.rich_output *= "\n" * s
end

function Base.display(report::Report, m::MIME"text/latex", data)
    s = repr(m, data)
    report.rich_output *= "\n" * s
end

"""Add saved figure name to results and return the name"""
function add_figure(report::Report, data, m, ext)
  chunk = report.cur_chunk
  full_name, rel_name = get_figname(report, chunk, ext = ext)

  open(full_name, "w") do io
      if ext == ".pdf"
          write(io, repr(m, data))
      else
          show(io, m, data)
      end
  end

  push!(report.figures, rel_name)
  report.fignum += 1
  return full_name
end

ext: jl


text:
using DataDeps

register(DataDep(
  "mammographic-mass",
  "https://archive.ics.uci.edu/ml/datasets/Mammographic+Mass",
  "https://archive.ics.uci.edu/ml/machine-learning-databases/mammographic-masses/mammographic_masses.data",
  "072051c896fb14f8f25c64eed30ba1333290af41e6c9ffad4e1ca573daa865a8",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=6,
        feature_indices=1:5,
        categoric_indices=3:4,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
function test_PartitionTS()
    d = 2
    lims = repeat([-1.0 1.0], d, 1)
    traj = Trajectory(d)
    retrain = [10 * 2^i for i in 1:5]
    pts = PartitionTS(1, retrain, traj, lims, 3, 100, p -> 4 * p, 2, false)
    f(x) = 5 * x[1]^3* x[2]^4 - 5 * x[2]^2 + 7 * x[1] * x[2]
    n = retrain[end]
    X, r = AdaptivePartition.gaussian_data(f, d, n, 1.0, lims)
    ContextualBandits.add_memory!(traj, n)
    for i in 1:n
        ContextualBandits.extend!(traj, X[i, :], 1, r[i], 0.0)
        update!(pts, X[i, :], 1, r[i])
    end
    m = partition_search(X, r, lims, 3, 100, p -> 4 * p, 0, false)
    @test m.scale ≈ pts.arms[1].scale
    @test choose([0.0, 0.0], pts) == 1
end

test_PartitionTS()


text:
"""
`Maker.VariableTarget`

The type created by `variable()`. Fields expected to be accessed publicly
include:

- `name::UTF8String`
- `dependencies::Vector{UTF8String}`
- `description::UTF8String`
"""
type VariableTarget <: AbstractTarget
    name::UTF8String
    dependencies::Vector{UTF8String}
    description::UTF8String
    action::Function
    funhash::UInt64
    isstale::Bool
    timestamp::DateTime
    varhash::UInt64
    m::Module
end

immutable CachedVariable <: AbstractCached
    funhash::UInt64
    timestamp::DateTime
    varhash::UInt64
end

cached(t::VariableTarget) = CachedVariable(t.funhash, t.timestamp, t.varhash)

function updatecache!(f::JLD.JldFile, t::VariableTarget)
    if t.name in names(f)
        delete!(f, t.name)
    end
    write(f, t.name, cached(t))
end

function execute(t::VariableTarget)
    x = has1arg(t) ? t.action(t) : t.action()
    eval(t.m, :($(symbol(t.name)) = $x))
    newhash = hash(x)
    if t.varhash != newhash
        t.timestamp = Dates.unix2datetime(time())
        t.varhash = newhash
    else
        ds = dependencies(t)
        if length(ds) > 0
            t.timestamp = max(t.timestamp, maximum([timestamp(d) for d in ds]))
        end
    end
    t.isstale = false
    updatecache!(t)
end

function isstale(t::VariableTarget)
    if !isdefined(t.m, symbol(t.name)) || t.isstale ||
       hash(eval(t.m, symbol(t.name))) != t.varhash
        return true
    end
    ds = dependencies(t)
    isempty(ds) ? true : timestamp(t) < maximum([timestamp(d) for d in ds])
end

timestamp(t::VariableTarget) = t.timestamp

"""
```julia
variable(action::Function, name::AbstractString, dependencies)
```
Define and register targets for Maker.jl.

- `action` is the function that operates when the target is
  executed.

- `name` refers to the name of the task or target.

- `dependencies` refers to the name (`AbstractString`) or names
  (`Vector{AbstractString}`) of targets that need to be satisfied
  for this target before running the `action`. These are also referred
  to as prerequisites.

Targets are registered globally.

`variable` targets define an action, and the result of the action will be
assigned to a global variable (within the Module where the  definition is
created) named by the argument `name`. A `variable` task keeps a timestamp
when the action runs. If the result of the action has the same hash as
the cached version of the hash, the cached timestamp is used. Using the
stored timestamp can reduce the number of unnecessary actions run.

If the `action` or `dependencies` of a target are redefined, the target will be
marked as stale, and the action will be updated at the next target check.

See also `make`, `file`, and `task`. `variable` registers a `VariableTarget` type.
"""
function variable(action::Function, name::AbstractString,
                  dependencies::AbstractArray = UTF8String[])
    dependencies = UTF8String[dependencies...]
    t = resolve(name, nothing)
    fh = funhash(action, dependencies)
    vh = 0
    isstale = true
    datetime = DateTime()
    if t === nothing
        # check the cache
        getjld() do f
            if name in names(f)
                x = read(f[utf8(name)])
                if fh == x.funhash
                    datetime = x.timestamp
                    vh = x.varhash
                    isstale = false
                end
            end
        end
    elseif fh != t.funhash || dependencies != t.dependencies
    else
        isstale = false
        datetime = t.timestamp
        vh = t.varhash
    end
    register(VariableTarget(name, dependencies, "", action, fh, isstale,
                            datetime, vh, current_module()))
end
variable(action::Function, name::AbstractString, dependencies::AbstractString) =
    variable(action, name, [utf8(dependencies)])


text:
json_schema = """
{ "irModel": {
  "irModelName": "string",
  "stateVariables": [{
    "name": "string",
    "label": "string",
    "type": "sv_type",
    "initial_value": "expression"
  }],
  "events": [{
    "name": "string",
    "label": "string",
    "rate": "expression",
    "input_predicate": {
      "enabling_condition": "expression"
    },
    "output_predicate": {
      "transition_function": ["lvalue = expression", "lvalue = expression", ...]
    }
  }],
  "constants": [{
    "name": "string",
    "value": "extern"|"expression"
  }],
  "expressions": [{
    "name": "string",
    "value": "extern"|"expression"
  }],
  "rateRewards": [{
    "name": "string",
    "variable": "string",
    "temporaltype": "instantoftime"|"intervaloftime"|"timeaveragedintervaloftime"|"steadystate",
    "samplingpoints": [{
      "time": "float"
      }]
  }],
  "impulseRewards": [{
    "name": "string",
    "event": "string",
    "temporaltype": "instantoftime"|"intervaloftime"|"timeaveragedintervaloftime"|"steadystate",
    "samplingpoints": [{
      "time": "float"
      }]
  }],
  "composedrewards": [{
    "name": "string",
    "expression": "expression"
  }]
 }}
"""

struct IRModel{T}
    name::String
    stateVariables::T
    events::E
    constants::C
    expressions::X
    rewards::R
end

struct Event{P,Q}
    name::String
    label::String
    rate::Expression
    input_predicate::P
    output_predicate::q
end

struct Input
    enabling_condition::Expression
end

struct Output
    transition::Expr
end

struct Constant{T}
    name::String
    value::T
end

struct Expression{T}
    name::String
    value::T
end

struct Rewards
    rate
    impulse
    composed
end

struct ComposedReward
    name
    expression
end

struct ImpulseReward
    name::String
    event::String
    temporaltype
    samplingpoints::Vector{Time}
end

struct RateRewards
    name
    variable
    temporaltype
    samplingpoints
end


text:
# FROM clause.

mutable struct FromClause <: AbstractSQLClause
    over::Union{SQLClause, Nothing}

    FromClause(; over = nothing) =
        new(over)
end

FromClause(over) =
    FromClause(over = over)

"""
    FROM(; over = nothing)
    FROM(over)

A `FROM` clause.

# Examples

```jldoctest
julia> c = ID(:person) |> AS(:p) |> FROM() |> SELECT((:p, :person_id));

julia> print(render(c))
SELECT "p"."person_id"
FROM "person" AS "p"
```
"""
FROM(args...; kws...) =
    FromClause(args...; kws...) |> SQLClause

dissect(scr::Symbol, ::typeof(FROM), pats::Vector{Any}) =
    dissect(scr, FromClause, pats)

function PrettyPrinting.quoteof(c::FromClause, ctx::QuoteContext)
    ex = Expr(:call, nameof(FROM))
    if c.over !== nothing
        ex = Expr(:call, :|>, quoteof(c.over, ctx), ex)
    end
    ex
end

rebase(c::FromClause, c′) =
    FromClause(over = rebase(c.over, c′))


text: include("NfMaxOrd/NfOrdCls.jl")


text:
using DependencyTrees
using Test

using DependencyTrees.TransitionParsing

const DT = DependencyTrees

using DependencyTrees: MultiWordTokenError, EmptyTokenError
using DependencyTrees: leftdeps, rightdeps, leftmostdep, rightmostdep
using DependencyTrees.TransitionParsing: stack, buffer, token, tokens

function showstr(op)
    buf = IOBuffer()
    show(buf, op)
    return String(take!(buf))
end

function test_treebank(filename)
    Treebank(joinpath(@__DIR__, "data", filename), DependencyTrees.from_conllu)
end

function test_sentence(filename)
    first(test_treebank(filename))
end

@testset "DependencyTrees" begin
    include("test_tokens.jl")
    include("test_graphs.jl")
    include("test_slp3_ch13.jl")
    include("test_nivre08.jl")
    include("test_kubleretal09.jl")
    include("test_treebanks.jl")
    include("test_conllu.jl")
    include("test_dynamic.jl")
    include("test_arc_hybrid.jl")
    include("test_arc_standard.jl")
    include("test_arc_swift.jl")
    include("test_evaluation.jl")
end


text:
"""
	gmtvector(cmd0::String="", arg1=nothing, kwargs...)

Time domain filtering of 1-D data tables.

Full option list at [`gmtvector`]($(GMTdoc)gmtvector.html)

Parameters
----------

- **A** | **single_vec** :: [Type => Str]   `Arg = m[conf]|vector`

    Specify a single, primary vector instead of reading tables.
    ($(GMTdoc)gmtvector.html#a)
- **C** | **cartesian** :: [Type => Str | []]        `Arg = [i|o]`

    Select Cartesian coordinates on input and output.
    ($(GMTdoc)gmtvector.html#c)
- **E** | **geod2geoc** :: [Type => Bool]

    Convert input geographic coordinates from geodetic to geocentric and output geographic
    coordinates from geocentric to geodetic.
    ($(GMTdoc)gmtvector.html#e)
- **N** | **normalize** :: [Type => Bool]

    Normalize the resultant vectors prior to reporting the output.
    ($(GMTdoc)gmtvector.html#n)
- **S** | **secondary_vec** :: [Type => Str | List]    `Arg = [vector]`

    Specify a single, secondary vector in the same format as the first vector.
    ($(GMTdoc)gmtvector.html#s)
- **T** | **transform** :: [Type => Str | List]     `Arg = a|d|D|paz|s|r[arg|R|x]`

    Specify the vector transformation of interest.
    ($(GMTdoc)gmtvector.html#t)
- $(GMT.opt_V)
- $(GMT.opt_write)
- $(GMT.opt_append)
- $(GMT.opt_b)
- $(GMT.opt_d)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_o)
- $(GMT.opt_swap_xy)
"""
function gmtvector(cmd0::String="", arg1=nothing; kwargs...)

	length(kwargs) == 0 && return monolitic("gmtvector", cmd0, arg1)

	d = KW(kwargs)
	cmd = parse_common_opts(d, "", [:V_params :b :d :e :f :g :h :i :o :yx])
	cmd = parse_these_opts(cmd, d, [[:A :single_vec], [:C :cartesian], [:E :geod2geoc], [:N :normalize],
				[:S :secondary_vec], [:T :transform]])

	common_grd(d, cmd0, cmd, "gmtvector ", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
gmtvector(arg1; kw...) = gmtvector("", arg1; kw...)


text:
# module Bukdu

export resources

using .Actions: routing_actions
import .Actions: index, edit, new, show, create, update, delete

const default_param_key = "id"

"""
    resources(path::String, ::Type{C}, param_types::Pair{Symbol,DataType}...; only=[], except=[]) where {C <: ApplicationController}
"""
function resources(path::String, ::Type{C}, param_types::Pair{Symbol,DataType}...; only=[], except=[]) where {C <: ApplicationController}
    param = default_param_key
    if !isempty(only)
        actions = only
    elseif !isempty(except)
        actions = setdiff(routing_actions, except)
    else
        actions = routing_actions
    end
    for (action, verb, routepath) in [(index,  get,    ""),
                                      (new,    get,    "/new"),
                                      (edit,   get,    "/:$param/edit"),
                                      (show,   get,    "/:$param"),
                                      (create, post,   ""),
                                      (delete, delete, "/:$param"),
                                      (update, patch,  "/:$param"),
                                      (update, put,    "/:$param")]
        !(action in actions) && continue
        url = string(path, routepath)
        Routing.add_route(verb, url, C, action, Dict{Symbol,DataType}(param_types...))
    end
end

# module Bukdu


text:
import Base.<
import Base.isapprox

Base.:<(x::T1, y::CertainValue{T2}) where {T1 <: Real, T2 <: Real} = x < y.value
Base.:<(x::CertainValue{T1}, y::T2) where {T1 <: Real, T2 <: Real} = x.value < y
Base.isless(x::CertainValue{T1}, y::CertainValue{T2}) where {T1 <: Real, T2 <: Real} = isless(x.value, y.value)
Base.isapprox(x::CertainValue{T1}, y::T2) where {T1 <: Real, T2 <: Real} = isapprox(x.value, y)
Base.isapprox(x::T1, y::CertainValue{T2}) where {T1 <: Real, T2 <: Real} = isapprox(x, y.value)


text:
include("../btg.jl")
include("btgBayesopt.jl")
include("../computation/finitedifference.jl")
using Plots
using Random
Random.seed!(8);
#
# trying out diff random seeds:
# failure: (3)
# success : (4, 5 works with good initialization, 6 works with init, 7 works with and w/o init.)
#use 12 for all prev experiments

#####
##### Define optimization problem, initialize GP with burn-in points, initialize btg parameters
#####

Himmelblau(x) = (x[1]^2 + x[2] -11)^2 + (x[1]+x[2]^2-7)^2 #function to optimize
y, x = sample_points(Himmelblau, [-5, -5], [5, 5]; num = 100) #burn-in points which are used as training data in GP
x = vcat(x, [5 5])
y = vcat(y, 890)
x = vcat(x, [-4.999 4.999])
y = vcat(y, 530)
x = vcat(x, [4.999 -4.999])
y = vcat(y, 610)
x = vcat(x, [-4.999 -4.999])
y = vcat(y, 249.6042059640019)
x = vcat(x, [5 2.55573])
y = vcat(y, 294.57103540062496)
x = vcat(x, [5.0 2.3475172404355718])
y = vcat(y, 279.5672977301108)

@info "y-vals", y[1:5]
### Set BTG parameters and get function handles for pdf, cdf, dpdf, etc.
Fx = linear_polynomial_basis(x)
#train = trainingData(x, Fx, y) #training data triple: location, covariates, labels
train = extensible_trainingData(x, Fx, y)
#rangeθ = reshape(select_single_theta_range(x), 1, 2) #rangetheta is auto-selected
#rangeθ = [0.0001 1000]
rangeθ = [300.0 1500]
@info "rangeθ", rangeθ
rangeλ = [-1.0 1.0]
btg1 = btg(train, rangeθ, rangeλ);
(pdf, cdf, dpdf, cdf_gradient, cdf_hessian) = solve(btg1; derivatives = true) #get function handles for pdf, cdf, derivtives

#####
##### Define acquisition function and derivatives
#####

println("defining acquisition function...")
#lx = [1, -5, -5]; ux = [3000.0, 5, 5] #box-constraints for optimization problem
lx = [1, -5, -5]; ux = [3000.0, 5, 5] #box-constraints for optimization problem
#make input an augmented vector [y, s]: label, location
cdf_fixed(v) = cdf(v[2:end], linear_polynomial_basis(v[2:end]), v[1])
cdf_gradient_fixed(v) = cdf_gradient(v[2:end], linear_polynomial_basis(v[2:end]), v[1])

d = length(lx)
fun(x...) = x[1] #x\in R^d+1
dfun(g, x...) = (g[1] = 1; g[i] = 0 for i = 2:length(g);)
#cdf_wrapper(x) = cdf_fixed(collect(x))
cdf_wrapper(x...) = cdf_fixed(collect(x))
#cdf_gradient_wrapper(g, x...) = (grad = cdf_gradient_fixed(collect(x...)); g[i] = grad[i] for i = 1:length(grad); @info grad; @info g; return nothing)
function cdf_gradient_wrapper(g, x...)
    grad = cdf_gradient_fixed(collect(x))
    @info "grad", grad
    for i = 1:length(grad)
        #@info "i", i
        #@info "grad component i", grad[i]
        global g[i] = grad[i]
    end
    return nothing
end
model = Model(Ipopt.Optimizer)
set_optimizer_attributes(model, "tol" => 1e-2, "max_iter" => 200)
#model = Model(Gurobi.Optimizer)
#NLoptSolver(algorithm=:SLSQP))
register(model, :fun, 3, fun, dfun)
register(model, :cdf_wrapper, 3, cdf_wrapper, cdf_gradient_wrapper)
#####
##### Initial guess for IPopt. Check gradient before starting optimizer
#####
#println("checking gradient at initial point...")
#####
##### Register variables and constraints in model and run optimizer.
#####
#@variable(model, lx[i] <= au[i=1:3] <= ux[i], start = initval[i])
@variable(model, lx[i] <= au[i=1:3] <= ux[i])
#@variable(model, au)
@NLobjective(model, Min, fun(au...))
@NLconstraint(model, cdf_wrapper(au...) == 0.25)

function single_optimization()
    initval = init_constrained_pt(cdf_fixed, lx, ux; quantile = 0.25)
    set_start_value(au[1], initval[1])
    set_start_value(au[2], initval[2])
    set_start_value(au[3], initval[3])
    #println("Running Optimizer...")
    #initval = [200, 4, 4]
    #@info "initval", initval
    #@variable(model, lx[i] <= au[i=1:3] <= ux[i])
    JuMP.optimize!(model)
    println(value.(au))
    vstar = value.(au)
    return (vstar, initval)
end

include("test_script.jl")
function run_func(n)
    x = [];
    vstars = [];
    init_val_cdf = []
    init_vals = []
    for i = 1:n
        Random.seed!(i+200)
        (vstar, initval) = single_optimization();
        push!(x, cdf_fixed(vstar));
        push!(vstars, vstar);
        push!(init_val_cdf, cdf_fixed(initval))
        push!(init_vals, initval)
    end
    return x, vstars, init_val_cdf, init_vals
end

if false
    (res, vstars, init_val_cdf, init_vals) = run_func(4);
    vcdfplot_sequence(vstars; res = res, upper = 2000) #shows cdfs at points we converged to
end

wait_for_key(prompt) = (print(stdout, prompt); read(stdin, 1); nothing)

if false #sanity check
    Plots.plot()
    Plots.scatter(x[:, 1], x[:, 2])
    wait_for_key("Observe the sampled points. Press any key to continue.")

    Plots.plot()
    func1 = x -> cdf([5.0 5.0], [1 5.0  5.0], x);
    plt!(func1, 1, 2000);
    println("25% quantile: ", fzero(x -> func1(x)-0.5, (1, 2000)))
    wait_for_key("CDF is large mean as expected. Press any key to continue.")

    func2 = x -> cdf([4 5], [1 4.0  5.0], x);
    Plots.plot()
    plt(func2, 1, 2000);
    println("25% quantile: ", fzero(x -> func2(x)-0.5, (1, 2000)))
    wait_for_key("CDF should be small. Press any key to continue.")

    func3 = x -> cdf([4 -2], [1 4.0  -2], x);
    Plots.plot()
    plt(func3, 1, 2000);
    println("25% quantile: ", fzero(x -> func3(x)-0.5, (1, 2000)))
    wait_for_key("CDF should be small. Press any key to continue.")

    #fzero(func)

end

if true #sanity check
    Plots.plot()
    Plots.scatter(x[:, 1], x[:, 2])

    Plots.plot()
    func1 = x -> cdf([5.0 5.0], [1 5.0  5.0], x);
    plt!(func1, 1, 2000);
    println("50% quantile: ", fzero(x -> func1(x)-0.5, (1, 2000)))
    println("actual: ", Himmelblau([5, 5]))

    func2 = x -> cdf([4 5], [1 4.0  5.0], x);
    plt(func2, 1, 2000);
    println("50% quantile: ", fzero(x -> func2(x)-0.5, (1, 2000)))
    println("actual: ", Himmelblau([4, 5]))

    func3 = x -> cdf([4 -2], [1 4.0  -2], x);
    plt(func3, 1, 2000);
    println("50% quantile: ", fzero(x -> func3(x)-0.5, (1, 2000)))
    println("actual: ", Himmelblau([4, -2]))

    func4 = x -> cdf([5 2.5], [1 5  2.5], x);
    plt(func4, 1, 2000);
    println("50% quantile: ", fzero(x -> func4(x)-0.5, (1, 2000)))
    println("actual: ", Himmelblau([5, 2.5]))
    #fzero(func)
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
# parser

function read_instance(MyFileName::String)
  path = "./" * MyFileName
  # Si le fichier path existe
  if isfile(path)
    # L’ouvrir
    myFile = open(path)
    readline(myFile)
    readline(myFile)
    readline(myFile)
    readline(myFile)

    N = parse(Int64, split(readline(myFile), " ")[3])
    Nm = parse(Int64, split(readline(myFile), " ")[3])
    Nf = parse(Int64, split(readline(myFile), " ")[3])
    C = parse(Int64, split(readline(myFile), " ")[3])
    G = parse(Int64, split(readline(myFile), " ")[3])
    A = parse(Int64, split(readline(myFile), " ")[3])
    T = parse(Int64, split(readline(myFile), " ")[3])
    init = parse(Float64, split(readline(myFile), " ")[3])

    y = Array{Int64,2}(zeros(N,A*G)) # on met a zero pour l'instant
    data = readlines(myFile)
		for datum in data
      line = split(datum, " ")
      allele = (parse(Int64,line[3]) - 1) * A + parse(Int64,line[5])
      y[parse(Int64,line[1]),allele] = y[parse(Int64,line[1]),allele] + 1
    end

    return N, Nm, Nf, C, G, A, T, init, y
  end
end

max_stars_repo_name: boumpteryx/RODD

max_issues_repo_name: boumpteryx/RODD

max_forks_repo_name: boumpteryx/RODD


text:
"""
"""
function autodiff_array_gradient(a,i_to_x,j_to_i=IdentityVector(length(i_to_x)))

  i_to_xdual = apply(i_to_x) do x
    cfg = ForwardDiff.GradientConfig(nothing, x, ForwardDiff.Chunk{length(x)}())
    xdual = cfg.duals
    xdual
  end

  j_to_f = to_array_of_functions(a,i_to_xdual,j_to_i)
  j_to_x = reindex(i_to_x,j_to_i)

  k = ForwardDiffGradientKernel()
  apply(k,j_to_f,j_to_x)

end

struct ForwardDiffGradientKernel <: Kernel end

function kernel_cache(k::ForwardDiffGradientKernel,f,x)
  cfg = ForwardDiff.GradientConfig(nothing, x, ForwardDiff.Chunk{length(x)}())
  r = copy(x)
  (r, cfg)
end

@inline function apply_kernel!(cache,k::ForwardDiffGradientKernel,f,x)
  r, cfg = cache
  @notimplementedif length(r) != length(x)
  ForwardDiff.gradient!(r,f,x,cfg)
  r
end

"""
"""
function autodiff_array_jacobian(a,i_to_x,j_to_i=IdentityVector(length(i_to_x)))

  i_to_xdual = apply(i_to_x) do x
    cfg = ForwardDiff.JacobianConfig(nothing, x, ForwardDiff.Chunk{length(x)}())
    xdual = cfg.duals
    xdual
  end

  j_to_f = to_array_of_functions(a,i_to_xdual,j_to_i)
  j_to_x = reindex(i_to_x,j_to_i)

  k = ForwardDiffJacobianKernel()
  apply(k,j_to_f,j_to_x)

end

struct ForwardDiffJacobianKernel <: Kernel end

function kernel_cache(k::ForwardDiffJacobianKernel,f,x)
  cfg = ForwardDiff.JacobianConfig(nothing, x, ForwardDiff.Chunk{length(x)}())
  n = length(x)
  j = zeros(eltype(x),n,n)
  (j, cfg)
end

@inline function apply_kernel!(cache,k::ForwardDiffJacobianKernel,f,x)
  j, cfg = cache
  @notimplementedif size(j,1) != length(x)
  @notimplementedif size(j,2) != length(x)
  ForwardDiff.jacobian!(j,f,x,cfg)
  j
end

"""
"""
function autodiff_array_hessian(a,i_to_x,j_to_i=IdentityVector(length(i_to_x)))
   agrad = i_to_y -> autodiff_array_gradient(a,i_to_y,j_to_i)
   autodiff_array_jacobian(agrad,i_to_x,j_to_i)
end

function to_array_of_functions(a,x,ids=IdentityVector(length(x)))
  k = ArrayOfFunctionsKernel(a,x)
  j = IdentityVector(length(ids))
  apply(k,j)
end

struct ArrayOfFunctionsKernel{A,X} <: Kernel
  a::A
  x::X
end

function kernel_cache(k::ArrayOfFunctionsKernel,j)
  xi = testitem(k.x)
  l = length(k.x)
  x = MutableFill(xi,l)
  ax = k.a(x)
  axc = array_cache(ax)
  (ax, x, axc)
end

@inline function apply_kernel!(cache,k::ArrayOfFunctionsKernel,j)
  ax, x, axc = cache
  @inline function f(xj)
    x.value = xj
    axj = getindex!(axc,ax,j)
  end
  f
end

mutable struct MutableFill{T} <: AbstractVector{T}
  value::T
  length::Int
end

Base.size(a::MutableFill) = (a.length,)

@inline Base.getindex(a::MutableFill,i::Integer) = a.value


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module PackageWithBuildSpecificTestDeps

f(x) = x

end

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
using Statistics
using PyPlot, CSV
using DataFrames, DataFramesMeta
import FileIO:load, save
# ---------
# Experimental Prediction 2:
# Process dataframe files (outputs of simulations)
# in order to create the csv data for the plots in Figure 11
#
# Sgm stands for Bayes Factor Surprise
# Ssh stands for Shannon Surprise
# ---------
function processAllLearners_second(;
                        filepath = "/path/to/experiment_jldfiles",
                            # foldername = "experimentPlots")
                            )

    jldfiles = searchdirendings(filepath, ".jld2")
    for i in 1:length(jldfiles)
        file = load(joinpath(filepath, jldfiles[i]))
        df = file["df"]
        plotAllExperimentsettings_second(df)
        csvAllExperimentsettings_second(df)
    end
end
# Make plots within julia: not used eventually for the paper, only for inspection
function plotAllExperimentsettings_second(df;
                                savepath = homedir()
                                )

    learnertype = unique(df[:learnertype])
    learnerparam = unique(df[:learnerparam])
    if length(learnertype) > 1
        error("There are more than one learners!")
    else
        learnername = learnertype[1]
    end
    savepathlearner = joinpath(savepath, string(learnername) * string(formatparamlist(learnerparam)) * "_exp2_plots")
    mkdir(savepathlearner)

    p_centers = df[:p_centers][1]

    dfmergedseeds = df_mergeseeds_Surprise_second(df)

    for σ in unique(df[:sigma])
        for p_c in unique(df[:changeprobability])
            for nstep in unique(df[:nsteps])
                for Δp in unique(df[:Δp])
                    dftrimmedmerged = df_experiment_trim_second(dfmergedseeds,
                                            svalue = σ,
                                            changeprobabilityvalue = p_c,
                                            nstepvalue = nstep,
                                            Δpvalue = Δp)
                    titlestring =  "σ: $(σ), pc: $(p_c), steps: $(nstep), Δp: $(Δp) "
                    figurestring = "σ$(σ)_" *
                                    "pc$(p_c)_" *
                                    "steps$(nstep)_" *
                                    "Δp$(Δp)"
                    meanSgm, stdSgm, seSgm = calc_stats_nan(dftrimmedmerged[:allmean_Sgm][1])
                    plotmeanSurpriseOverp(meanSgm, seSgm, titlestring = titlestring,
                                        ylabelstring = "Mean Bayes Factor Surprise",
                                        xdata = p_centers
                                        )
                    PyPlot.savefig(joinpath(savepathlearner,
                                    "exp_second_Sgm"* string(learnername) *
                                                        figurestring * ".png"))
                    PyPlot.close()
                    plt.ioff()
                    meanSsh, stdSsh, seSsh = calc_stats_nan(dftrimmedmerged[:allmean_Ssh][1])
                    plotmeanSurpriseOverp(meanSsh, seSsh, titlestring = titlestring,
                                        ylabelstring = "Mean Shannon Surprise",
                                        xdata = p_centers,
                                        ylims = [0., 2.]
                                        )
                    PyPlot.savefig(joinpath(savepathlearner,
                                    "exp_second_Ssh"* string(learnername) *
                                                        figurestring * ".png"))
                    PyPlot.close()
                    plt.ioff()
                end
            end
        end
    end
end
function plotmeanSurpriseOverp(mean_Surprise, se_Surprise;
                            xlabelstring = "p",
                            ylabelstring = "Mean S",
                            titlestring = "",
                            ylims = [0., 3.],
                            xdata = []
                            )
    fig = figure(); ax = gca()
    plot_shadedlines([mean_Surprise],
                    [se_Surprise],
                    [""], ax,
                    xlabelstring = xlabelstring,
                    titlestring = titlestring,
                    ylabelstring = ylabelstring,
                    ylims = ylims,
                    xdata = xdata
                    )
end
# Create the csv data for the plots for Experimental Prediction 1
function csvAllExperimentsettings_second(df;
                                savepath = homedir()
                                )

    learnertype = unique(df[:learnertype])
    learnerparam = unique(df[:learnerparam])
    if length(learnertype) > 1
        error("There are more than one learners!")
    else
        learnername = learnertype[1]
    end
    savepathlearner = joinpath(savepath, string(learnername) * string(formatparamlist(learnerparam)) * "_exp2_csv")
    mkdir(savepathlearner)

    p_centers = df[:p_centers][1]

    dfmergedseeds = df_mergeseeds_Surprise_second(df)

    for σ in unique(df[:sigma])
        for p_c in unique(df[:changeprobability])
            for nstep in unique(df[:nsteps])
                for Δp in unique(df[:Δp])
                    dfresults = DataFrame()
                    dfresults[:x] = p_centers
                    dftrimmedmerged = df_experiment_trim_second(dfmergedseeds,
                                                        svalue = σ,
                                                        changeprobabilityvalue = p_c,
                                                        nstepvalue = nstep,
                                                        Δpvalue = Δp)
                    csvfilestring = "sigma$(σ)_" *
                                    "pc$(p_c)_" *
                                    "steps$(nstep)_" *
                                    "Dp$(Δp)"

                    meanSgm, stdSgm, seSgm = calc_stats_nan(dftrimmedmerged[:allmean_Sgm][1])

                    dfresults[Symbol(string(learnername)*"_Sgm_mean")] = meanSgm
                    dfresults[Symbol(string(learnername)*"_Sgm_std")] = stdSgm
                    dfresults[Symbol(string(learnername)*"_Sgm_se")] = seSgm

                    meanSsh, stdSsh, seSsh = calc_stats_nan(dftrimmedmerged[:allmean_Ssh][1])
                    dfresults[Symbol(string(learnername)*"_Ssh_mean")] = meanSsh
                    dfresults[Symbol(string(learnername)*"_Ssh_std")] = stdSsh
                    dfresults[Symbol(string(learnername)*"_Ssh_se")] = seSsh

                    CSV.write(joinpath(savepathlearner, "exp_second_"*
                                                string(learnername) *
                                                csvfilestring *".csv"),
                                                dfresults, delim = " ")
                end
            end
        end
    end
end

# -----------------
# Handle dataframes
# -----------------
function df_mergeseeds_Surprise_second(df)
    df = by(df, [:learnertype, :learnerparam, :sigma, :changeprobability, :nsteps, :Δp],
                allmean_Sgm = :mean_Sgm => x -> [collect(reduce(hcat, x))],
                allse_Sgm = :se_Sgm => x -> [collect(reduce(hcat, x))],
                allmean_Ssh = :mean_Ssh => x -> [collect(reduce(hcat, x))],
                allse_Ssh = :se_Ssh => x -> [collect(reduce(hcat, x))],
                )
end
function df_experiment_trim_second(df; svalue = 0.5, changeprobabilityvalue = 0.05,
                    nstepvalue = 1000,
                    Δpvalue = 0.0125
                    )
    dftrimmed = @where(df, :changeprobability .== changeprobabilityvalue,
                            :sigma .== svalue,
                            :nsteps .== nstepvalue,
                            :Δp .== Δpvalue)
end


text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

function nopotentialgrad(foldername, filename, u, xij, xik, xil, xi, xj, xk, xl, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, pot)

ifile = 0;
gen = 0;
if pot==0
    stropu, strcpu, strgpu = gendensitygrad(foldername,filename, u, mu, mu, eta, kappa, gen, ifile);
elseif pot==1
    stropu, strcpu, strgpu = gensinglegrad(foldername, filename, u, xi, qi, ti, ai, mu, eta, kappa, gen, ifile);
elseif pot==2
    stropu, strcpu, strgpu = genpairgrad(foldername, filename, u, xij, qi, qj, ti, tj, ai, aj, mu, eta, kappa, gen, ifile);
elseif pot==3
    stropu, strcpu, strgpu = gentripletgrad(foldername, filename, u, xij, xik, qi, qj, qk, ti, tj, tk, ai, aj, ak, mu, eta, kappa, gen, ifile);
elseif pot==4
    stropu, strcpu, strgpu = genquadrupletgrad(foldername, filename, u, xij, xik, xil, qi, qj, qk, ql, ti, tj, tk, tl, ai, aj, ak, al, mu, eta, kappa, gen, ifile);
end
cppfiles(foldername, filename, stropu, strcpu, strgpu, 1);

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
@testset "Creation of an identity multiple" begin
    for n in [1, 2, 1000]
        In = IdentityMultiple(1.0I, n) # same as IdentityMultiple(UniformScaling(1.0), n))
        @test size(In) == (n, n)
        @test In[1, 1] == 1.0 && In[1, 2] == 0.0
    end
end

@testset "Operations between identity multiples" begin
    I2 = IdentityMultiple(UniformScaling(1.0), 2)
    I10 = IdentityMultiple(UniformScaling(1.0), 10)

    @test (I2 + I2).M == UniformScaling(2.0)
    @test (10.0 * I2).M == UniformScaling(10.0)
    @test (I2 * I2).M == I2.M

    @test_throws AssertionError I2 + I10
    @test_throws AssertionError I2 * I10
end

@testset "Create a continuous system with one matrix being a multiple of the identity" begin
    # See #40
    A, B = rand(4, 4), rand(4, 2)
    X = rand(Hyperrectangle, dim=4)
    U = rand(Ball2, dim=2)
    I4 = I(4)
    s = ConstrainedLinearControlContinuousSystem(A, I4, X, B*U);

    @test statedim(s) == 4
    @test stateset(s) == X
end

@testset "Matrix operations for identity multiple" begin
    I2 = I(2)
    @test getindex(I2, 1) == 1.
    @test getindex(I2, 2) == 0.
    @test getindex(I2, 3) == 0.
    @test getindex(I2, 4) == 1.
    x = I2 * rand(2)
    @test x isa Vector && length(x) == 2
    A = I2 * rand(2, 2)
    @test A isa Matrix && size(A) == (2, 2)
end


text:
numgrad(m, f, xs...) = FiniteDifferences.grad(m, f, xs...)
function numgrad(m, f, ps::Flux.Params)
    ccp = deepcopy(ps)
    function gp(p)
        cp = deepcopy(p)
        g = FiniteDifferences.grad(m, x -> (p .= x; f()), p) |> only
        p .= cp
        g
    end
    tuple(map(gp, ps)...)
end

# For points that are close to discontinuities or non-differentiable parts
function numgrads(f, xs...)
    bg = numgrad(backward_fdm(7, 1, max_range=1e-5), f, xs...)
    cg = numgrad(central_fdm(7, 1, max_range=1e-5), f, xs...)
    fg = numgrad(forward_fdm(7, 1, max_range=1e-5), f, xs...)
    zip(bg, cg, fg) |> collect
end

gradcheck(ag, ng; atol=1e-5, rtol=1e-5) = isnothing(ag) ? isapprox.(ng, 0; atol) : isapprox.(ng, ag; atol, rtol)
gradcheck(ag, ngs::Tuple; atol=1e-5, rtol=1e-5) = all(.|(map(ng -> gradcheck(ag, ng; atol, rtol), ngs)...))
gradcheck(ag::NotImplemented, ngs::Tuple; atol=1e-5, rtol=1e-5) = throw(NotImplementedException(ag))

function gradf(f, xs...)
    A = rand(size(f(xs...))...)
    (xs...) -> sum(A .* f(xs...))
end
function gradf(f, ps::Flux.Params)
    A = rand(size(f())...)
    () -> sum(A .* f())
end

function gradtest(f, xs...; atol=1e-5, rtol=1e-5)
    f2 = gradf(f, xs...)
    ag = Flux.gradient(f2, xs...)
    ngs = numgrads(f2, xs...)
    all(i -> gradcheck(ag[i], ngs[i]; atol, rtol), eachindex(xs))
end
function gradtest(f, ps::Flux.Params; atol=1e-5, rtol=1e-5)
    f2 = gradf(f, ps)
    ag = Flux.gradient(f2, ps)
    ngs = numgrads(f2, ps)
    all(i -> gradcheck(ag[ps[i]], ngs[i]; atol, rtol), 1:length(ps))
end


text:
import Base.Docs: Binding

# Keywords
# ––––––––

const builtins = ["abstract", "baremodule", "begin", "bitstype", "break",
                  "catch", "ccall", "const", "continue", "do", "else",
                  "elseif", "end", "export", "finally", "for", "function",
                  "global", "if", "immutable", "import", "importall", "let",
                  "local", "macro", "module", "quote", "return", "try", "type",
                  "typealias", "using", "while"]

const builtin_completions =
  [d(:text=>k, :type=>:keyword)
   for k in builtins]

# Module completions
# ––––––––––––––––––

const identifier_pattern = r"^@?[_\p{L}][\p{Xwd}′!]*+$"

_names(mod; all = false, imported = false) = filter!(x -> !Base.isdeprecated(mod, Symbol(x)), names(mod, all=all, imported=imported))

moduleusings(mod) = ccall(:jl_module_usings, Any, (Any,), mod)

filtervalid(names) = @>> names map(string) filter(x->occursin(identifier_pattern, x))

accessible(mod::Module) =
  [_names(mod, all=true, imported=true);
   map(_names, moduleusings(mod))...] |> unique |> filtervalid

Base.getindex(b::Binding) = isdefined(b.mod, b.var) ? getfield(b.mod, b.var) : nothing

completiontype(x) =
  isa(x, Module)   ? "module" :
  isa(x, DataType) ? "type"   :
  isa(x, Function) ? "λ"      :
                     "constant"

const meta_cache = Dict{Any,Any}()

function withmeta(completion::AbstractString, mod::Module)
  isdefined(mod, Symbol(completion)) || return completion
  b = Binding(mod, Symbol(completion))
  mod = b.mod
  haskey(meta_cache, (mod, completion)) && return meta_cache[(mod, completion)]
  x = b[]
  c = d(:text => completion,
        :type => startswith(completion, "@") ? "macro" : completiontype(x),
        :rightLabel => string(mod))
  c[:displayText] = signature(b)
  c[:description] = description(b)
  c
end

for name in map(string, _names(Base))
  meta_cache[(Base, name)] = withmeta(name, Base)
end

withmeta(completions::Vector, mod::Module) =
  [withmeta(completion, mod) for completion in completions]

function namecompletions_(mod::Module, qualified = false)
  if !qualified
    [withmeta(accessible(mod), mod); builtin_completions]
  else
    withmeta(filtervalid(_names(mod, all=true)), mod)
  end
end

const namecompletions = memoize_debounce(namecompletions_)

const prefix_pattern = r"(@?[_\p{L}][\p{Xwd}!]*+\.?@?)+$|@$"

function prefix(line)
  match = Base.match(prefix_pattern, line)
  match == nothing && return String[]
  split(match.match, ".")
end

# Path Completions
# ––––––––––––––––

function funcprefix(line)
  m = match(r"\b([^\s]*)\(([^\)]*)$", line) # matches `foo(bar`
  m == nothing && return
  return m.captures[1], m.captures[2]
end

function pathprefix(line)
  m = match(r"\"?([\w\\/:\-.]*)$", line)
  m == nothing && return
  return m.match, m.captures[1]
end

function stringmeta(cs, prefix)
  map(c -> Dict(:text => c, :_prefix => prefix, :type => :file), cs)
end

function pathmeta(cs, path, prefix)
  stringmeta(map(c -> replace(joinpath(path, c), r"^\./" => ""), cs), prefix)
end

function children(dir, ext = ""; depth = 0, out = String[], prefix = "")
  isdir(dir) || return []
  for f in readdir(dir)
    path = joinpath(dir, f)
    if isfile(path) && ext ≠ nothing && endswith(path, ext)
      push!(out, joinpath(prefix, f))
    elseif isdir(path)
      push!(out, joinpath(prefix, f, ""))
      if depth > 0
        children(path, ext, depth = depth-1, out = out, prefix = joinpath(prefix, f))
      else
      end
    end
  end
  return out
end

function includepath()
  path = Base.source_path()
  path == nothing && return pwd()
  return dirname(path)
end

function pathcompletions(line)
  m = funcprefix(line)
  m == nothing && return
  func, pre = m
  m = pathprefix(pre)
  m == nothing && return
  pre, path = m
  dir = dirname(path)
  try
    if func == "include"
      pathmeta(children(joinpath(includepath(), dir), ".jl", depth = 2), dir, path)
    elseif func == "readcsv"
      pathmeta(children(joinpath(pwd(), dir), ".csv", depth = 2), dir, path)
    elseif func == "cd"
      pathmeta(children(joinpath(pwd(), dir), nothing), dir, path)
    elseif func == "open"
      pathmeta(children(joinpath(pwd(), dir)), dir, path)
    end
  catch e
    []
  end
end

# Package Completions
# –––––––––––––––––––
# not sure why this errors with "ArgumentError: Module Pkg3 not found in current path."

import Pkg

pkgmeta(xs) = [Dict(:text => x, :type => "package") for x in xs]

function _stdlibs()
  ver = string("v", VERSION.major, ".", VERSION.minor)
  srcdir = joinpath(Sys.BINDIR,"..","..","stdlib",ver)
  releasedir = joinpath(Sys.BINDIR,"..","share","julia","stdlib",ver)
  return ispath(srcdir) ? readdir(srcdir) :
                          ispath(releasedir) ? readdir(releasedir) : []
end

const stdlibs = memoize(_stdlibs)

function pkgcompletions(line)
  if occursin(r"^using|^import", line)
    return pkgmeta(sort!(append!(collect(keys(Pkg.API.installed())), stdlibs())))
  end
end

# Completions
# –––––––––––

const providers = [pkgcompletions, pathcompletions]

function completions(line::AbstractString, mod::Module = Main; default = true)
  for provider in providers
    cs = provider(line)
    cs ≠ nothing && return cs
  end
  pre = prefix(line)
  mod = getthing(mod, pre[1:end-1])
  if isa(mod, Module)
    qualified = length(pre) > 1
    !(qualified || default) && return
    return namecompletions(mod, qualified)
  end
  return []
end

completions(mod::Module) = completions("", mod, default = true)


text:
abstract type Interpolator end

"""
    Batch a bunch of nonuniform interpolators
"""
struct BatchedInterpolator{F,T} <: Interpolator
    batch_size::Int
    input_dim::Int
    interpolators::F
    t_span::Tuple{T,T}
    function BatchedInterpolator(interpolators)
        batch_size = length(interpolators)

        # assume all interpolators match first, use to set bounds and datatypes
        t_span = interpolators[1].t_span
        new{typeof(Tuple(interpolators)),typeof(t_span[1])}(
            batch_size,
            input_dim,
            Tuple(interpolators),
            t_span,
        )
    end
end

function (f::BatchedInterpolator)(t)
    return reduce(hcat, [interp(t) for interp in f.interpolators])
end

"""
    times: 1d Array length L
    data: 2d array size (L, ...)
"""
struct LinearInterpolator{T} <: Interpolator
    times::Array{T,1}
    data::Array{T,2}
    t_span::Tuple{T,T}
    function LinearInterpolator(times, data)
        t_span = (times[1], times[end])
        new{typeof(times[1])}(times, data, t_span)
    end
end

function (f::LinearInterpolator)(t)
    idx = findfirst(x -> x >= t, f.times) - 1
    idx == 0 ? idx += 1 : nothing
    t_left = f.times[idx]
    t_delta = f.times[idx+1] - t_left
    θ = (t - t_left) / t_delta
    return (one(t) - θ) * f.data[:, idx] + θ * f.data[:, idx+1]
end

""" Interpolation functions """

function linear_interpolate(
    t::T,
    times::AbstractArray{T,1},
    data::AbstractArray{T,2},
) where {T}
    idx = findfirst(x -> x >= t, times) - 1
    idx == 0 ? idx += 1 : nothing
    t_left = times[idx]
    t_delta = times[idx+1] - t_left
    θ = (t - t_left) / t_delta
    return (one(t) - θ) * data[:, idx] .+ θ * data[:, idx+1]
end

function linear_interpolate(
    t::T,
    times::AbstractArray{T,1},
    data::AbstractArray{T,3},
) where {T}
    idx = findfirst(x -> x >= t, times) - 1
    idx == 0 ? idx += 1 : nothing
    t_left = times[idx]
    t_delta = times[idx+1] - t_left
    θ = (t - t_left) / t_delta
    return (one(t) - θ) * data[:, idx, :] .+ θ * data[:, idx+1, :]
end

""" TODO: add more interpolation types """
# https://github.com/PumasAI/DataInterpolations.jl/blob/master/src/interpolation_methods.jl
# https://github.com/PumasAI/DataInterpolations.jl/blob/master/src/interpolation_caches.jl

max_stars_repo_name: syanga/ISSNeuralODE

max_issues_repo_name: syanga/ISSNeuralODE

max_forks_repo_name: syanga/ISSNeuralODE


text:
using Printf

function map_filter_reduce_examples()
    println("--------------------------------------")
    println("map_filter_reduce_examples")

    println("We can declare an integer array")
    integer_array = [1,2,3,4,5,6,7,8,9,10]
    println(integer_array)

    println("We can use filter() to get the even numbers")
    even_numbers = filter(n -> ((n % 2) == 0), integer_array)
    println(even_numbers)

    println("We can use map() to double every number")
    doubled_numbers = map(n -> (n * 2), integer_array)
    println(doubled_numbers)

    println("We can use reduce() to get the sum of every number")
    sum_numbers = reduce(+, integer_array)
    println(sum_numbers)
end

max_stars_repo_name: James-P-D/JuliaDump

max_issues_repo_name: James-P-D/JuliaDump

max_forks_repo_name: James-P-D/JuliaDump


text:
using StaticArrays, ForwardDiff, BenchmarkTools, LinearAlgebra

n,m,N = 3,2,11
dt = 0.1
x = @SVector [1.,2,3,]
u = @SVector [-1.,1]
z = [x;u]
Z = [KnotPoint(x,u,dt) for k = 1:N]

xc = @SVector [0,1,2,3.]
yc = @SVector [0,2,2,0.]
r = 0.5*@SVector ones(length(xc))
vals = [@SVector zeros(length(xc)) for k = 1:N]
circlecon = CircleConstraint(n,m,xc,yc,r)
generate_jacobian(circlecon)

vals2 = [@SVector zeros(1) for k = 1:N]
normcon = NormConstraint(n,m,1.0)
generate_jacobian(normcon)

@btime evaluate($circlecon, $x, $u)
@btime evaluate($normcon, $x, $u)

@btime jacobian($circlecon, $z)
@btime jacobian($normcon, $z)

∇c = [@SMatrix zeros(length(circlecon), n+m) for k = 1:N]
@btime jacobian!($∇c, $circlecon, $Z)
∇c = [@SMatrix zeros(length(normcon), n+m) for k = 1:N]
@btime jacobian!($∇c, $normcon, $Z)

con1 = ConstraintVals(circlecon, 1:N)
con2 = ConstraintVals(normcon, 1:4)

@btime evaluate!($con1.vals, $circlecon, $Z)
@btime evaluate!($con1, $Z)
con2.∇c

jacobian!(con2.∇c, normcon, Z)

@btime evaluate!($con1, $Z)
@btime jacobian!($con1, $Z)
@btime jacobian!($con2, $Z)

@btime update_active_set!($con1, $(Val(0.01)))
@btime duals($con1, 3)

@btime max_violation!($con1)
@btime max_violation!($con2)

constraints = [con1, con2]
conSet = ConstraintSet(constraints, N)
jacobian!(con1, Z)
jacobian!(con2, Z)
for con in conSet.constraints
	jacobian!(con, Z)
end

@btime jacobian!($conSet, $Z)

evaluate!(conSet, Z)
jacobian!(conSet, Z)

max_violation!(conSet)
@btime max_violation!($conSet)

# AL Objective
xf = @SVector zeros(n)
Q = Diagonal(@SVector ones(n))
R = Diagonal(@SVector ones(m))
Qf = Diagonal(@SVector ones(n))
obj = LQRObjective(Q,R,Qf,xf,N)

alobj = ALObjective3(obj, conSet)
E = CostExpansion(n,m,N)
@btime cost_expansion($E, $obj, $Z)
cost_expansion(E,alobj,Z)
@btime cost_expansion($Q,$alobj,$Z)

cost_expansion(E, con1, Z)
@btime cost_expansion($E, $con1, $Z)
J = obj.J
@btime cost!($J, $con1, $Z)
@btime cost!($alobj, $Z)

# Custom constraints

struct CustomConstraint{C} <: AbstractConstraint{C}
	n::Int
	m::Int
	p::Int
	c::Function
	∇c::Function
end

function CustomConstraint{C}(n::Int, m::Int, p::Int, c::Function) where C
	generate_jacobian(c, n, m)
	∇c = eval(Symbol("∇" * string(c)))
	CustomConstraint{C}(n,m,p,c,∇c)
end

evaluate(con::CustomConstraint, x, u) = con.c(x,u)
evaluate(con::CustomConstraint, z) = con.c(z)
jacobian(con::CustomConstraint, x, u) = con.∇c(x, u)
jacobian(con::CustomConstraint, z) = con.∇c(z)

# Test constraints
n,m = 3,2
p1 = 3

x = @SVector [1.,2,3]
u = @SVector [-5.,5]
z = [x;u]

x_max = @SVector [3,2,Inf]
x_min = @SVector [-3,-2,-Inf]
u_max = @SVector [Inf, 0.5]
u_min = @SVector [0.0, -0.5]

bnd = BoundConstraint(n,m, x_max=x_max, x_min=x_min, u_max=u_max, u_min=u_min)

@btime evaluate($bnd, $xs, $us)
@btime evaluate($bnd, $xs)

function mycon1(x, u)
	@SVector [x[1]^2 + x[2]^2 - 5, u[1] - 1, u[2] - 1]
end
jacob_eq1(x,u) = [2x[1] 2x[2] 0 0 0;
                0     0     0 1 0;
                0     0     0 0 1];
eq2 = CustomConstraint{Equality}(n,m,3,mycon1)

evaluate(eq2, x, u) == [0,-6,4]

@btime evaluate($eq2,$x,$u)
@btime jacobian($eq2, $z)

mycon2(x,u) = @SVector [sin(x[1]), sin(x[3])]
ineq2 = CustomConstraint{Inequality}(n,m,2,mycon2)
@btime evaluate($ineq2, $x, $u)
@btime jacobian($ineq2, $z)

constraints = [bnd, eq2, ineq2]
constraints = [eq2, eq2]

out = [@SVector zeros(3) for k = 1:2]
function eval_constraints!(out, con, x, u)
	out[2] = evaluate(con, x, u)
end
constraints = (eq2, eq2)
@btime eval_constraints!($out, $eq2, $x, $u)
@btime evaluate($eq2, $x, $u)
@btime evaluate($constraints[2], $x, $u)

constraints = @SVector [eq2, eq2]
@btime evaluate($bnd, $x, $u)

# Quadrotor constraint
n,m = 13,4
dt = 0.1
xs,us = (@SVector rand(n)), (@SVector rand(m))
x,u = Array(xs), Array(us)
zs = [xs;us]
z = [x;u]
Z = KnotPoint(x,u,dt)

struct CircleConstraint{T} <: AbstractConstraint{Inequality}
	n::Int
	m::Int
	p::Int
	r0::SVector{3,T}
	radius::T
end

Base.size(con::CircleConstraint) = (con.n, con.m, con.p)

function evaluate(con::CircleConstraint, x, u)
	p = con.r0
	r = con.radius
	return @SVector [-((x[1] - p[1])^2 + (x[2]-p[2])^2  - r^2)]
end

function evaluate(con::CircleConstraint, z)
	p = con.r0
	r = con.radius
	return @SVector [-((z[1] - p[1])^2 + (z[2]-p[2])^2  - r^2)]
end

r0 = @SVector [0.5, 0.5, 0]
radius = 0.1
con = CircleConstraint(n, m, 1, r0, radius)
@btime evaluate($con, $xs, $us)

generate_jacobian(con)
jacobian(con, zs)
jacobian(con, xs, us)
@btime jacobian($con, $zs)

x_max = @SVector ones(n)
x_min = -@SVector ones(n)
u_max = 2@SVector ones(m)
u_min = -2@SVector ones(m)

bnd = BoundConstraint(n,m, x_max=x_max, x_min=x_min, u_max=u_max, u_min=u_min)

constraints = (con, con, bnd)

Z = [KnotPoint(x,u,dt) for k = 1:10]

function eval_constraints!(cval, constraint, Z::Traj)
	for k in eachindex(cval)
		cval[k] = evaluate(constraint, state(Z[k]), control(Z[k]))
	end
end

function eval_all_constraints!(cvals, constraints, Z::Traj)
	eval_constraints!(cvals[1], constraints[1], Z)
	# eval_constraints!(cvals[2], constraints[2], Z)
	# eval_constraints!(cvals[3], constraints[3], Z)
	# for i in eachindex(cvals)
	# 	eval_constraints!(cvals[i], constraints[i], Z)
	# end
end

@btime eval_constraints!($(cvals[1]), $(constraints[1]), $Z)
constraints = (bnd, bnd, con)
cvals = [[@SVector zeros(size(con)[3]) for k = 1:10] for con in constraints]
@btime eval_all_constraints!($cvals, $constraints, $Z)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "git_crypt"
version = v"0.6.0"

# Collection of sources required to complete build
sources = [
    "https://github.com/AGWA/git-crypt/archive/$(version).tar.gz" =>
    "777c0c7aadbbc758b69aff1339ca61697011ef7b92f1d1ee9518a8ee7702bb78",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/git-crypt-*/
make CPPFLAGS="-I${prefix}/include" LDFLAGS="-L${libdir} -lcrypto"
make install PREFIX=$prefix
if [ $target = "x86_64-w64-mingw32" ] || [ $target = "i686-w64-mingw32" ]; then
    mv $bindir/git-crypt $bindir/git-crypt.exe
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    ExecutableProduct("git-crypt", :git_crypt)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "OpenSSL_jll",
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"7")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using MultiDimDictionaries
using Test

@testset "MultiDimDictionaries.jl" begin
  d = MultiDimDictionary(([(1, 1), (2, 1), (1, 2)]), [1, 2, 3])

  # Get elements

  @test d[1, 1] == 1
  @test d[2, 1] == 2
  @test d[1, 2] == 3

  @test d[(1, 1)] == 1
  @test d[(2, 1)] == 2
  @test d[(1, 2)] == 3

  @test d[LinearIndex(1)] == 1
  @test d[LinearIndex(2)] == 2
  @test d[LinearIndex(3)] == 3

  # Error
  # @test d[(1, 1)] == 1
  # @test d[(2, 1)] == 2
  # @test d[(1, 2)] == 3

  # Set elements

  d[1, 1] = 4
  d[(2, 1)] = 5
  d[LinearIndex(3)] = 6

  @test d[1, 1] == 4
  @test d[2, 1] == 5
  @test d[1, 2] == 6

  insert!(d, (2, 2), 12)

  @test d[2, 2] == 12

  delete!(d, (1, 1))

  @test !haskey(d, (1, 1))
  @test haskey(d, (1, 2))
  @test haskey(d, (2, 1))
  @test haskey(d, (2, 2))

  delete!(d, 1, 2)

  @test !haskey(d, 1, 1)
  @test !haskey(d, 1, 2)
  @test haskey(d, 2, 1)
  @test haskey(d, 2, 2)

  d2 = MultiDimDictionary([(3, 2), (2, 3)], [32, 23])
  d12 = merge(d, d2)

  @test d12 isa MultiDimDictionary
  @test !haskey(d12, 1, 2)
  @test haskey(d12, 2, 1)
  @test haskey(d12, 3, 2)

  d1 = MultiDimDictionary()
  d1[1] = 1
  d1[2] = 2

  @test d1[1] == 1
  @test d1[2] == 2

  d2 = MultiDimDictionary()
  d2[1] = 3
  d2[2] = 4

  @test d2[1] == 3
  @test d2[2] == 4

  dv = [d1; d2]

  @test dv[1] == 1
  @test dv[2] == 2
  @test dv[3] == 3
  @test dv[4] == 4

  dh = hcat(d1, d2)

  @test dh[1, 1] == 1
  @test dh[2, 1] == 2
  @test dh[1, 2] == 3
  @test dh[2, 2] == 4

  dh = hcat(d1, d2; new_dim_keys=("X", "Y"))

  @test dh[1, "X"] == 1
  @test dh[2, "X"] == 2
  @test dh[1, "Y"] == 3
  @test dh[2, "Y"] == 4

  dh_X = dh[:, "X"]

  @test dh_X[1, "X"] == 1
  @test dh_X[2, "X"] == 2
  @test !isassigned(dh_X, 1, "Y")
  @test !isassigned(dh_X, 2, "Y")

  dh_Y = dh[:, "Y"]

  @test dh_Y[1, "Y"] == 3
  @test dh_Y[2, "Y"] == 4

  dh_Z = dh[:, "Z"]

  @test !isassigned(dh_Z, 1, "X")
  @test !isassigned(dh_Z, 2, "X")
  @test !isassigned(dh_Z, 1, "Y")
  @test !isassigned(dh_Z, 2, "Y")

  d = d1 ⊔ d2

  @test d[1, 1] == 1
  @test d[1, 2] == 2
  @test d[2, 1] == 3
  @test d[2, 2] == 4

  d = disjoint_union(d1, d2; new_dim_keys=("X", "Y"))

  @test d["X", 1] == 1
  @test d["X", 2] == 2
  @test d["Y", 1] == 3
  @test d["Y", 2] == 4

  d_slice = d[[("X", 1), ("Y", 2)]]

  @test d_slice["X", 1] == 1
  @test !isassigned(d_slice, "X", 2)
  @test !isassigned(d_slice, "Y", 1)
  @test d_slice["Y", 2] == 4

  d2 = (d ⊔ d)[1, :]

  @test d2[1, "X", 1] == 1
  @test d2[1, "X", 2] == 2
  @test d2[1, "Y", 1] == 3
  @test d2[1, "Y", 2] == 4
end


text:
# see examples.jl first

# Type Definition
@aml mutable struct Person "person", check_course
    age::UInt64, "~"
    field, "study-field"
    GPA::Float64 = 4.5, "~", GPAcheck
    courses::Vector{String}, "taken-courses"
    professors::UN{DataFrame} = nothing, "table"
    id::Int64, att"~"
    comment::UN{String} = nothing, txt"end"
end

@aml mutable struct University doc"university"
    name, att"university-name"
    people::Vector{Person}, "person"
end

# Value Checking Functions
GPAcheck(x) = x <= 4.5 && x >= 0

function check_course(age, field, GPA, courses, professors, id, comment)

    if field == "Mechanical Engineering"
        relevant = ["Artificial Intelligence", "Robotics", "Machine Design"]
    elseif field == "Computer Engineering"
        relevant = ["Julia", "Algorithms"]
    else
        error("study field is not known")
    end

    return any(in.(courses, Ref(relevant)))
end

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

@aml mutable struct MyGeneralXML{T} "my-general-xml"
    myfield::T, "~"
end
pxml_string = MyGeneralXML{String}(myfield = "a")
pprint(pxml_string)
pxml_vector = MyGeneralXML{Vector{String}}(myfield = ["b","c"])
pprint(pxml_vector)


text:
using CUDA

function matmul(a, b, c, d, conf;
                transform_global_to_shared_a = Transform.Elementwise(),
                transform_global_to_shared_b = Transform.Elementwise(),
                transform_global_to_shared_c = Transform.Elementwise(),
                transform_shared_to_global_d = Transform.Elementwise(),
                transform_shared_to_regs_a = Transform.Elementwise(),
                transform_shared_to_regs_b = Transform.Elementwise(),
                transform_shared_to_regs_c = Transform.Elementwise(),
                transform_regs_to_shared_d = Transform.Elementwise(),
                epilogue = Epilogue.Default(),
                kernel = Kernel.matmul_singlestage)

    args = [a, b, c, d,
            transform_global_to_shared_a, transform_global_to_shared_b, transform_global_to_shared_c, transform_shared_to_global_d,
            transform_shared_to_regs_a, transform_shared_to_regs_b, transform_shared_to_regs_c, transform_regs_to_shared_d,
            epilogue,
            conf]

    GC.@preserve args begin
        kernel_args = cudaconvert.(args)
        kernel_tt = Tuple{Core.Typeof.(kernel_args)...}
        kernel = cufunction(kernel, kernel_tt; )
        attributes(kernel.fun)[CUDA.FUNC_ATTRIBUTE_MAX_DYNAMIC_SHARED_SIZE_BYTES] = 64 * 1024
        kernel(kernel_args...; conf.launch_args...)
    end
end


text:
export initialize, getstream, setstream!, send

###### Defaults for AbstractDevice
#What to do when connecting
initialize(::AbstractDevice) = nothing
#What to do when disconnecting
close(::AbstractDevice) = nothing
#Which stream is the AbstractDevice connected to
getstream(dev::AbstractDevice) = dev.stream
#Set the stream the AbstractDevice is connected to
setstream!(dev::AbstractDevice, stream::LabStream) = dev.stream = stream

function safe_getwritecommand(dev::AbstractDevice, val)
    stream = try getstream(dev) catch
        error("Device $dev not connected to a stream")
    end
    cmd = try getwritecommand(stream, dev, val) catch
        error("Device $dev with output val $val not supported on stream $stream")
    end
    return cmd, stream
end
function safe_getreadcommand(dev::AbstractDevice)
    stream = try getstream(dev) catch
        error("Device $dev not connected to a stream")
    end
    cmd = try getreadcommand(stream, dev) catch
        error("Device $dev not supported  to read from on stream $stream")
    end
    return cmd, stream
end

function send(dev::AbstractDevice, val)
    cmd, stream = safe_getwritecommand(dev, val)
    send(stream, cmd)
    return
end
function read(dev::AbstractDevice)
    cmd, stream = safe_getreadcommand(dev)
    vals, ts = read(stream, cmd)
    return vals
end

function put!(dev::AbstractDevice, val)
    cmd, stream = safe_getwritecommand(dev, val)
    push!(stream.sendbuffer,cmd)
    return
end

function get(dev::AbstractDevice)
    cmd, stream = safe_getreadcommand(dev)
    push!(stream.readbuffer,cmd)
    return
end


text:
using HTTP.MessageRequest: body_was_streamed

struct DownloadsBackend <: AWS.AbstractBackend
    downloader::Union{Nothing,Downloads.Downloader}
end

DownloadsBackend() = DownloadsBackend(nothing)

const AWS_DOWNLOADER = Ref{Union{Nothing,Downloader}}(nothing)
const AWS_DOWNLOAD_LOCK = ReentrantLock()

# Here we mimic Download.jl's own setup for using a global downloader.
# We do this to have our own downloader (separate from Downloads.jl's global downloader)
# because we add a hook to avoid redirects in order to try to match the HTTPBackend's
# implementation, and we don't want to mutate the global downloader from Downloads.jl.
# https://github.com/JuliaLang/Downloads.jl/blob/84e948c02b8a0625552a764bf90f7d2ee97c949c/src/Downloads.jl#L293-L301
function get_downloader(downloader=nothing)
    lock(AWS_DOWNLOAD_LOCK) do
        yield() # let other downloads finish
        downloader isa Downloader && return nothing
        while true
            downloader = AWS_DOWNLOADER[]
            downloader isa Downloader && return nothing
            D = Downloader()
            D.easy_hook =
                (easy, info) -> Curl.setopt(easy, Curl.CURLOPT_FOLLOWLOCATION, false)
            AWS_DOWNLOADER[] = D
        end
    end
    return downloader
end

# https://github.com/JuliaWeb/HTTP.jl/blob/2a03ca76376162ffc3423ba7f15bd6d966edff9b/src/MessageRequest.jl#L84-L85
body_length(x::AbstractVector{UInt8}) = length(x)
body_length(x::AbstractString) = sizeof(x)

read_body(x::IOBuffer) = take!(x)
function read_body(x::IO)
    close(x)
    return read(x)
end

function _http_request(backend::DownloadsBackend, request::Request, response_stream::IO)
    # If we pass `output`, older versions of Downloads.jl will
    # expect a message body in the response. Specifically, it sets
    # <https://curl.se/libcurl/c/CURLOPT_NOBODY.html>
    # only when we do not pass the `output` argument.
    # See <https://github.com/JuliaLang/Downloads.jl/issues/131>.
    #
    # When the method is `HEAD`, the response may have a Content-Length
    # but not send any content back (which appears to be correct,
    # <https://stackoverflow.com/a/18925736/12486544>).
    #
    # Thus, if we did not set `CURLOPT_NOBODY`, and it gets a Content-Length
    # back, it will hang waiting for that body.
    #
    # Therefore, we do not pass an `output` when the `request_method` is `HEAD`.
    # (Note: this is fixed on the latest Downloads.jl, but we include this workaround
    #  for compatability).
    output = request.request_method != "HEAD" ? response_stream : nothing

    # HTTP.jl sets this header automatically.
    request.headers["Content-Length"] = string(body_length(request.content))

    # We pass an `input` only when we have content we wish to send.
    input = !isempty(request.content) ? IOBuffer(request.content) : nothing

    @repeat 4 try
        downloader = @something(backend.downloader, get_downloader())
        # set the hook so that we don't follow redirects. Only
        # need to do this on per-request downloaders, because we
        # set our global one with this hook already.
        if backend.downloader !== nothing && downloader.easy_hook === nothing
            downloader.easy_hook =
                (easy, info) -> Curl.setopt(easy, Curl.CURLOPT_FOLLOWLOCATION, false)
        end

        # We seekstart on every attempt, otherwise every attempt
        # but the first will send an empty payload.
        input !== nothing && seekstart(input)

        response = @mock Downloads.request(
            request.url;
            input=input,
            output=output,
            method=request.request_method,
            headers=request.headers,
            verbose=false,
            throw=true,
            downloader=downloader,
        )

        http_response = HTTP.Response(
            response.status, response.headers; body=body_was_streamed, request=nothing
        )

        if HTTP.iserror(http_response)
            target = HTTP.resource(HTTP.URI(request.url))
            throw(
                HTTP.StatusError(
                    http_response.status, request.request_method, target, http_response
                ),
            )
        end

        return AWS.Response(http_response, response_stream)
    catch e
        @delay_retry if (
            (isa(e, HTTP.StatusError) && AWS._http_status(e) >= 500) ||
            isa(e, Downloads.RequestError)
        )
        end
    end
end


text:
module Commons
  function backsolve(X,Y)

    row = size(X)[2]
    betacol = zeros(eltype(Y),row,1)

    iterator =row
    while iterator>0

      betacol[iterator] = (Y[iterator]- (X[iterator,:]*betacol)[1])/X[iterator,iterator]

      iterator-=1
    end

    return betacol

  end

end

max_stars_repo_name: ajermaky/CustomMLLibrary

max_issues_repo_name: ajermaky/CustomMLLibrary

max_forks_repo_name: ajermaky/CustomMLLibrary


text:
using ModelingToolkit: operation
# Method of lines discretization scheme
struct MOLFiniteDifference{T} <: DiffEqBase.AbstractDiscretization
    dxs::T
    order::Int
    MOLFiniteDifference(args...;order=2) = new{typeof(args[1])}(args[1],order)
end

# Get boundary conditions from an array
function get_bcs(bcs,tdomain,domain)
    lhs_deriv_depvars_bcs = Dict()
    no_bcs = size(bcs,1)
    for i = 1:no_bcs
        var = operation(bcs[i].lhs)
        if var isa Sym
            var = var.name
            if !haskey(lhs_deriv_depvars_bcs,var)
                lhs_deriv_depvars_bcs[var] = Array{Expr}(undef,3)
            end
            j = 0
            if isequal(bcs[i].lhs.args[1],tdomain.lower) # u(t=0,x)
                j = 1
            elseif isequal(bcs[i].lhs.args[2],domain.lower) # u(t,x=x_init)
                j = 2
            elseif isequal(bcs[i].lhs.args[2],domain.upper) # u(t,x=x_final)
                j = 3
            end
            if !(bcs[i].rhs isa ModelingToolkit.Symbolic)
                lhs_deriv_depvars_bcs[var][j] = :(var=$(bcs[i].rhs))
            else
                lhs_deriv_depvars_bcs[var][j] = toexpr(bcs[i].rhs)
            end
        end
    end
    return lhs_deriv_depvars_bcs
end

# Recursively traverses the input expression (rhs), replacing derivatives by
# finite difference schemes. It returns a time dependent expression (expr)
# that will be evaluated in the "f" ODE function (in DiffEqBase.discretize),
# Note: 'non-derived' dependent variables are inserted into the differential equations
#       E.g., Dx(u(t,x))=v(t,x)*Dx(u(t,x)), v(t,x)=t*x
#            =>  Dx(u(t,x))=t*x*Dx(u(t,x))

function discretize_2(input,deriv_order,approx_order,dx,X,len,
                      deriv_var,dep_var_idx,indep_var_idx)
    if !(input isa ModelingToolkit.Symbolic)
        return :($(input))
    else
        if input isa Sym || (input isa Term && operation(input) isa Sym)
            expr = :(0.0)
            var = nameof(input isa Sym ? input : operation(input))
            if haskey(indep_var_idx,var) # ind. var.
                if var != :(t)
                    i = indep_var_idx[var]
                    expr = :($X[$i][2:$len[$i]-1])
                else
                    expr = :(t)
                end
            else # dep. var.
                # TODO: time and cross derivatives terms
                i = indep_var_idx[deriv_var[1]]
                j = dep_var_idx[var]
                if deriv_order == 0
                    expr = :(u[:,$j])
                elseif deriv_order == 1
                    # TODO: approx_order and forward/backward should be
                    #       input parameters of each derivative
                    approx_order = 1
                    L = UpwindDifference(deriv_order,approx_order,dx[i],len[i]-2,-1)
                    expr = :(-1*($L*Q[$j]*u[:,$j]))
                elseif deriv_order == 2
                    L = CenteredDifference(deriv_order,approx_order,dx[i],len[i]-2)
                    expr = :($L*Q[$j]*u[:,$j])
                end
            end
            return expr
        elseif input isa Term && operation(input) isa Differential
            var = nameof(input.op.x)
            push!(deriv_var,var)
            return discretize_2(input.args[1],deriv_order+1,approx_order,dx,X,
                                len,deriv_var,dep_var_idx,indep_var_idx)
            pop!(deriv_var,var)
        else
            name = nameof(operation(input))
            if size(input.args,1) == 1
                aux = discretize_2(input.args[1],deriv_order,approx_order,dx,X,
                                   len,deriv_var,dep_var_idx,indep_var_idx)
                return :(broadcast($name, $aux))
            else
                aux_1 = discretize_2(input.args[1],deriv_order,approx_order,dx,X,
                                     len,deriv_var,dep_var_idx,indep_var_idx)
                aux_2 = discretize_2(input.args[2],deriv_order,approx_order,dx,X,
                                     len,deriv_var,dep_var_idx,indep_var_idx)
                return :(broadcast($name, $aux_1, $aux_2))
            end
        end
    end
end

# Convert a PDE problem into an ODE problem
function DiffEqBase.discretize(pdesys::PDESystem,discretization::MOLFiniteDifference)

    # TODO: discretize the following cases
    #
    #   1) PDE System
    #        1.a) Transient
    #                There is more than one indep. variable, including  't'
    #                E.g., du/dt = d2u/dx2 + d2u/dy2 + f(t,x,y)
    #        1.b) Stationary
    #                There is more than one indep. variable, 't' is not included
    #                E.g., 0 = d2u/dx2 + d2u/dy2 + f(x,y)
    #   2) ODE System
    #        't' is the only independent variable
    #        The ODESystem is packed inside a PDESystem
    #        E.g., du/dt = f(t)
    #
    #   Note: regarding input format, lhs must be "du/dt" or "0".
    #

    # The following code deals with 1.a case for 1-D,
    # i.e., only considering 't' and 'x'.

    ### Declare and define independent-variable data structures ###############

    tdomain = 0.0
    indep_var_idx = Dict()
    no_indep_vars = size(pdesys.domain,1)
    domain = Array{Any}(undef,no_indep_vars)
    dx = Array{Any}(undef,no_indep_vars)
    X = Array{Any}(undef,no_indep_vars)
    len = Array{Any}(undef,no_indep_vars)
    k = 0
    for i = 1:no_indep_vars
        var = nameof(pdesys.domain[i].variables)
        indep_var_idx[var] = i
        domain[i] = pdesys.domain[i].domain
        if var != :(t)
            dx[i] = discretization.dxs[i-k]
            X[i] = domain[i].lower:dx[i]:domain[i].upper
            len[i] = size(X[i],1)
        else
            dx[i] = 0.0
            X[i] = 0.0
            len[i] = 0.0
            tdomain = pdesys.domain[1].domain
            @assert tdomain isa IntervalDomain
            k = 1
        end
    end

    ### Declare and define dependent-variable data structures #################

    # TODO: specify order for each derivative
    approx_order = discretization.order

    lhs_deriv_depvars = Dict()
    dep_var_idx = Dict()
    dep_var_disc = Dict() # expressions evaluated in the ODE function (f)

    # if there is only one equation
    if pdesys.eq isa Equation
        eqs = [pdesys.eq]
    else
        eqs = pdesys.eq
    end
    no_dep_vars = size(eqs,1)
    for j = 1:no_dep_vars
        input = eqs[j].lhs
        op = operation(input)
        if op isa Sym
            var = nameof(op)
        else #var isa Differential
            var = nameof(operation(input.args[1]))
            lhs_deriv_depvars[var] = j
        end
        dep_var_idx[var] = j
    end
    for (var,j) in dep_var_idx
        aux = discretize_2( eqs[j].rhs,0,approx_order,dx,X,len,
                            [],dep_var_idx,indep_var_idx)
        # TODO: is there a better way to convert an Expr into a Function?
        dep_var_disc[var] = @eval (Q,u,t) -> $aux
    end

    ### Declare and define boundary conditions ################################

    # TODO: extend to Neumann BCs and Robin BCs
    lhs_deriv_depvars_bcs = get_bcs(pdesys.bcs,tdomain,domain[2])
    t = 0.0
    u_t0 = Array{Float64}(undef,len[2]-2,no_dep_vars)
    u_x0 = Array{Any}(undef,no_dep_vars)
    u_x1 = Array{Any}(undef,no_dep_vars)
    Q = Array{RobinBC}(undef,no_dep_vars)

    for var in keys(dep_var_idx)
        j = dep_var_idx[var]
        bcs = lhs_deriv_depvars_bcs[var]

        g = eval(:((x,t) -> $(bcs[1])))
        u_t0[:,j] = @eval $g.($(X[2][2:len[2]-1]),$t)

        u_x0[j] = @eval (x,t) -> $(bcs[2])
        u_x1[j] = @eval (x,t) -> $(bcs[3])

        a = Base.invokelatest(u_x0[j],X[2][1],0.0)
        b = Base.invokelatest(u_x1[j],last(X[2]),0.0)
        Q[j] = DirichletBC(a,b)
    end

    ### Define the discretized PDE as an ODE function #########################

    function f(du,u,p,t)

        # Boundary conditions can vary with respect to time
        for j in 1:no_dep_vars
            a = Base.invokelatest(u_x0[j],X[2][1],t)
            b = Base.invokelatest(u_x1[j],last(X[2]),t)
            Q[j] = DirichletBC(a,b)
        end

        for (var,disc) in dep_var_disc
            j = dep_var_idx[var]
            res = Base.invokelatest(disc,Q,u,t)
            if haskey(lhs_deriv_depvars,var)
                du[:,j] = res
            else
                u[:,j] .= res
            end
        end

    end

    # Return problem ##########################################################
    return PDEProblem(ODEProblem(f,u_t0,(tdomain.lower,tdomain.upper),nothing),Q,X)
end


text:
include("utils.jl")

@testset "Hyperbolic Space" begin
    M = Hyperbolic(2)
    @testset "Hyperbolic Basics" begin
        @test repr(M) == "Hyperbolic(2)"
        @test base_manifold(M) == M
        @test manifold_dimension(M) == 2
        @test typeof(get_embedding(M)) ==
              MetricManifold{ℝ,Euclidean{Tuple{3},ℝ},MinkowskiMetric}
        @test representation_size(M) == (3,)
        @test isinf(injectivity_radius(M))
        @test isinf(injectivity_radius(M, ExponentialRetraction()))
        @test isinf(injectivity_radius(M, [0.0, 0.0, 1.0]))
        @test isinf(injectivity_radius(M, [0.0, 0.0, 1.0], ExponentialRetraction()))
        @test !is_manifold_point(M, [1.0, 0.0, 0.0, 0.0])
        @test !is_tangent_vector(M, [0.0, 0.0, 1.0], [0.0, 0.0, 1.0, 0.0])
        @test_throws DomainError is_manifold_point(M, [2.0, 0.0, 0.0], true)
        @test !is_manifold_point(M, [2.0, 0.0, 0.0])
        @test !is_tangent_vector(M, [1.0, 0.0, 0.0], [1.0, 0.0, 0.0])
        @test Manifolds.default_metric_dispatch(M, MinkowskiMetric()) === Val{true}()
        @test_throws DomainError is_tangent_vector(
            M,
            [1.0, 0.0, 0.0],
            [1.0, 0.0, 0.0],
            true,
        )
        @test !is_tangent_vector(M, [0.0, 0.0, 1.0], [1.0, 0.0, 1.0])
        @test_throws DomainError is_tangent_vector(
            M,
            [0.0, 0.0, 1.0],
            [1.0, 0.0, 1.0],
            true,
        )
        @test is_default_metric(M, MinkowskiMetric())
        @test Manifolds.default_metric_dispatch(M, MinkowskiMetric()) === Val{true}()
        @test manifold_dimension(M) == 2

        for (P, T) in zip(
            [HyperboloidPoint, PoincareBallPoint, PoincareHalfSpacePoint],
            [HyperboloidTVector, PoincareBallTVector, PoincareHalfSpaceTVector],
        )
            p = convert(P, [1.0, 0.0, sqrt(2.0)])
            X = convert(T, [1.0, 0.0, sqrt(2.0)], [0.0, 1.0, 0.0])
            @test number_eltype(p) == eltype(p.value)
            @test X * 2.0 == T(X.value * 2.0)
            @test 2 \ X == T(2 \ X.value)
            @test +X == T(+X.value)
            @test Manifolds.allocate_result_type(M, log, (p, p)) == T
            @test Manifolds.allocate_result_type(M, inverse_retract, (p, p)) == T
            convert(T, p, X) == X
        end
    end
    @testset "Hyperbolic Representation Conversion I" begin
        p = [0.0, 0.0, 1.0]
        pH = HyperboloidPoint(p)
        @test minkowski_metric(pH, pH) == minkowski_metric(p, p)
        @test convert(HyperboloidPoint, p).value == pH.value
        @test convert(AbstractVector, pH) == p
        X = [1.0, 0.0, 0.0]
        XH = HyperboloidTVector(X)

        @test convert(HyperboloidTVector, X).value == XH.value
        @test convert(AbstractVector, XH) == X
        @test convert(HyperboloidPoint, p).value == pH.value
        is_manifold_point(M, pH)
        pB = convert(PoincareBallPoint, p)
        @test pB.value == convert(PoincareBallPoint, pH).value
        @test is_manifold_point(M, pB)
        @test convert(AbstractVector, pB) == p # convert back yields again p
        @test convert(HyperboloidPoint, pB).value == pH.value
        @test_throws DomainError is_manifold_point(
            M,
            PoincareBallPoint([0.9, 0.0, 0.0]),
            true,
        )
        @test_throws DomainError is_manifold_point(M, PoincareBallPoint([1.0, 0.0]), true)

        @test is_tangent_vector(M, pB, PoincareBallTVector([2.0, 2.0]))

        pS = convert(PoincareHalfSpacePoint, p)
        pS2 = convert(PoincareHalfSpacePoint, pB)
        pS3 = convert(PoincareHalfSpacePoint, pH)

        @test_throws DomainError is_manifold_point(
            M,
            PoincareHalfSpacePoint([0.0, 0.0, 1.0]),
            true,
        )
        @test_throws DomainError is_manifold_point(
            M,
            PoincareHalfSpacePoint([0.0, -1.0]),
            true,
        )

        @test pS.value == pS2.value
        @test pS.value == pS3.value
        @test convert(AbstractVector, pS) == convert(HyperboloidPoint, pS).value
        @test convert(PoincareBallPoint, pS2).value == pB.value
    end
    @testset "Hyperbolic Representation Conversion II" begin
        M = Hyperbolic(2)
        pts = [[0.0, 0.0, 1.0], [1.0, 0.0, sqrt(2.0)]]
        X = log(M, pts[2], pts[1])
        # For HyperboloidTVector we can do a plain wrap/unwrap
        X1 = convert(HyperboloidTVector, X)
        @test convert(AbstractVector, X1) == X
        # Convert to types and back to Array
        for (P, T) in zip(
            [HyperboloidPoint, PoincareBallPoint, PoincareHalfSpacePoint],
            [HyperboloidTVector, PoincareBallTVector, PoincareHalfSpaceTVector],
        )
            # convert to P,T
            p1 = convert(P, pts[2])
            X1 = convert(T, pts[2], X)
            (p2, X2) = convert(Tuple{P,T}, (pts[2], X))
            @test isapprox(M, p1, p2)
            @test isapprox(M, p1, X1, X2)
            for (P2, T2) in zip(
                [HyperboloidPoint, PoincareBallPoint, PoincareHalfSpacePoint],
                [HyperboloidTVector, PoincareBallTVector, PoincareHalfSpaceTVector],
            )
                @test isapprox(M, convert(P2, p1), convert(P2, pts[2]))
                @test convert(T, p1, X1) == convert(T, pts[2], X)
                (p3, X3) = convert(Tuple{P2,T2}, (pts[2], X))
                (p3a, X3a) = convert(Tuple{P2,T2}, (p1, X1))
                @test isapprox(M, p3, p3a)
                @test isapprox(M, p3, X3, X3a)
                @test isapprox(M, convert(P2, p2), p3)
                @test isapprox(M, pts[2], convert(AbstractVector, p3))
                @test isapprox(M, p3, convert(T2, p2, X2), X3)
                # coupled conversion
                (pT, XT) = convert(Tuple{AbstractVector,AbstractVector}, (p2, X2))
                @test isapprox(M, pts[2], pT)
                @test isapprox(M, pts[2], X, XT)
            end
        end
    end

    types = [
        Vector{Float64},
        SizedVector{3,Float64},
        HyperboloidPoint,
        PoincareBallPoint,
        PoincareHalfSpacePoint,
    ]
    TEST_FLOAT32 && push!(types, Vector{Float32})
    for T in types
        @testset "Type $T" begin
            is_plain_array =
                T ∉ [HyperboloidPoint, PoincareBallPoint, PoincareHalfSpacePoint]
            pts = [
                convert(T, [0.0, 0.0, 1.0]),
                convert(T, [1.0, 0.0, sqrt(2.0)]),
                convert(T, [0.0, 1.0, sqrt(2.0)]),
            ]
            test_manifold(
                M,
                pts,
                test_project_tangent=is_plain_array || T == HyperboloidPoint,
                test_musical_isomorphisms=is_plain_array,
                test_default_vector_transport=true,
                vector_transport_methods=[
                    ParallelTransport(),
                    SchildsLadderTransport(),
                    PoleLadderTransport(),
                ],
                is_tangent_atol_multiplier=10.0,
                exp_log_atol_multiplier=10.0,
                retraction_methods=(ExponentialRetraction(),),
                test_vee_hat=false,
                test_forward_diff=is_plain_array,
                test_reverse_diff=is_plain_array,
                test_tangent_vector_broadcasting=is_plain_array,
                test_vector_spaces=is_plain_array,
            )
        end
    end
    @testset "Embedding test" begin
        p = [0.0, 0.0, 1.0]
        X = [1.0, 1.0, 0.0]
        @test embed(M, p) == p
        q = similar(p)
        embed!(M, q, p)
        @test q == p
        @test embed(M, p, X) == X
        Y = similar(X)
        embed!(M, Y, p, X)
        @test Y == X
        p2 = HyperboloidPoint(p)
        X2 = HyperboloidTVector(X)
        q2 = HyperboloidPoint(similar(p))
        @test embed(M, p2).value == p2.value
        embed!(M, q2, p2)
        @test q2.value == p2.value
        @test embed(M, p2, X2).value == X2.value
        Y2 = HyperboloidTVector(similar(X))
        embed!(M, Y2, p2, X2)
        @test Y2.value == X2.value
    end
    @testset "Hyperbolic mean test" begin
        pts = [
            [0.0, 0.0, 1.0],
            [1.0, 0.0, sqrt(2.0)],
            [-1.0, 0.0, sqrt(2.0)],
            [0.0, 1.0, sqrt(2.0)],
            [0.0, -1.0, sqrt(2.0)],
        ]
        ws = UnitWeights{Float64}(length(pts))
        @test isapprox(M, mean(M, pts), pts[1]; atol=10^-4)
        @test isapprox(M, mean(M, pts, ws), pts[1]; atol=10^-4)
        @test_throws DimensionMismatch mean(M, pts, UnitWeights{Float64}(length(pts) + 1))
    end
    @testset "Hyperbolic ONB test" begin
        M = Hyperbolic(2)
        p = Manifolds._hyperbolize(M, [1.0, 0.0])
        B = get_basis(M, p, DefaultOrthonormalBasis())
        V = get_vectors(M, p, B)
        for v in V
            @test is_tangent_vector(M, p, v, true)
            for b in [DefaultOrthonormalBasis(), DiagonalizingOrthonormalBasis(V[1])]
                @test isapprox(M, p, v, get_vector(M, p, get_coordinates(M, p, v, b), b))
            end
        end
        for v in V, w in V
            @test inner(M, p, v, w) ≈ (v == w ? 1 : 0)
        end
        X = 0.5 * V[1] + 1.0 .* V[2]
        @test is_tangent_vector(M, p, X)
        c = get_coordinates(M, p, X, B)
        @test c ≈ [0.5, 1.0]
        B2 = DiagonalizingOrthonormalBasis(X)
        V2 = get_vectors(M, p, get_basis(M, p, B2))
        @test V2[1] ≈ X ./ norm(M, p, X)
        @test inner(M, p, V2[1], V2[2]) ≈ 0.0 atol = 5e-16
        B3 = DiagonalizingOrthonormalBasis(-V[2])
        V3 = get_vectors(M, p, get_basis(M, p, B3))
        @test isapprox(M, p, V3[1], -V[2])
    end
    @testset "Hyperboloid Minkowski Metric tests" begin
        p = Manifolds._hyperbolize(M, [0.0, 1.0])
        X = [1.0, 0.0, 0.0]
        pH = HyperboloidPoint(p)
        XH = HyperboloidTVector(X)
        @test minkowski_metric(XH, pH) == minkowski_metric(X, p)
        @test minkowski_metric(pH, XH) == minkowski_metric(p, X)
        @test minkowski_metric(XH, XH) == minkowski_metric(X, X)
    end
    @testset "PoincareBall project tangent test" begin
        p = PoincareBallPoint([0.5, 0.0])
        X = PoincareBallTVector([0.2, 0.2])
        Y = PoincareBallTVector(allocate(X.value))
        project!(M, Y, p, X)
        @test Y == X
        Z = project(M, p, X)
        @test isa(Z, PoincareBallTVector)
        @test Z == X
        p2 = PoincareHalfSpacePoint([1.1, 1.1])
        X2 = PoincareHalfSpaceTVector([0.2, 0.2])
        Y2 = PoincareHalfSpaceTVector(allocate(X2.value))
        project!(M, Y2, p2, X2)
        @test Y2 == X2
        Z2 = project(M, p2, X2)
        @test isa(Z2, PoincareHalfSpaceTVector)
        @test Z2 == X2
    end
end


text:
function astronaut_simulation(mech::Mechanism;
		tsim=1.0,
		tctrl=1.0,
		seed=0,
		ϵ=1.0e-14,
		control_amplitude=0.05)

    Random.seed!(seed)
    initialize!(mech, :humanoid)

	function ctrl!(mechanism, k)
		nu = input_dimension(mech)
		u = (k*mechanism.timestep < tctrl) * control_amplitude * [szeros(6); srand(nu-6)]
		set_input!(mech, u)
	    return
	end

    tcompute = @elapsed storage = simulate!(mech, tsim, ctrl!,
		record=true,
		opts=SolverOptions(rtol=ϵ, btol=ϵ))

    return storage, tcompute
end

function astronaut_simulation(;
	Nsim=1,
	timestep=1.0e-2,
	gravity=0.0,
	spring=0.0,
	damper=0.0,
	tsim=1.0,
	tctrl=1.0,
	seed=0,
	ϵ=1.0e-14,
	control_amplitude=0.1)

    mech = get_mechanism(:humanoid,
		timestep=timestep,
		gravity=gravity,
		spring=spring,
		damper=damper,
		contact=false)

	storage = []
	tcompute = zeros(Nsim)

	for i = 1:Nsim
		s, tcompute[i] = astronaut_simulation(mech;
			tsim=tsim,
			tctrl=tctrl,
			seed=seed+i,
			ϵ=ϵ,
			control_amplitude=control_amplitude)
		push!(storage, s)
	end
	return mech, [storage...], tcompute
end


text:
using ArgParse
using Glob
using Base.Filesystem

function getargs()
    s = ArgParseSettings()
    @add_arg_table! s begin
        "rootdir"
            help = "root directory of the github repository"
            required = true
        "outdir"
            help = "root directory of the github repository"
            required = true
    end
    args = parse_args(s)
end

getuttid(path) = splitext(basename(path))[1]
getspeaker(uttid) = split(uttid, "_")[1]

function main()
    args = getargs()

    inputdir = abspath(args["rootdir"])
    outputdir = abspath(args["outdir"])

    for x in ["train", "dev"]
        @info "preparing $x set..."

        run(`mkdir -p $(joinpath(outputdir, x))`)

        @info "creating wav.scp file"
        uttids = Set()
        open(joinpath(outputdir, x, "wav.scp"), "w") do scp
            for path in glob("*wav", joinpath(inputdir, "full_corpus_newsplit", x))
                uttid = getuttid(path)
                println(scp, "$uttid $path")
                push!(uttids, uttid)
            end
        end

        @info "creating uttids file"
        open(joinpath(outputdir, x, "uttids"), "w") do f
            for u in uttids
                println(f, "$u")
            end
        end

        @info "creating speaker mapping"
        open(joinpath(outputdir, x, "utt2spk"), "w") do utt2spk
            for uttid in uttids
                spk = getspeaker(uttid)
                println(utt2spk, "$uttid $spk")
            end
        end

        @info "preparing phonetic transcription"
        open(joinpath(outputdir, x, "trans.phn"), "w") do trans
            for path in glob("*mb", joinpath(inputdir, "full_corpus_newsplit", x))
                uttid = getuttid(path)
                text = filter(l -> !isspace(l), lowercase(readline(path)))
                phns = [l for l in text]
                println(trans, "$uttid $(join(phns, " "))")
            end
        end

        @info "preparing word transcription"
        open(joinpath(outputdir, x, "trans.wrd"), "w") do trans
            for path in glob("*mb", joinpath(inputdir, "full_corpus_newsplit", x))
                uttid = getuttid(path)
                text = lowercase(readline(path))
                words = [w for w in split(text, " ")]
                println(trans, "$uttid $(join(words, " "))")
            end
        end

        @info "prepared $(length(uttids)) utterances for $x set"
    end

    @info "preparing the lang directory..."
    run(`mkdir -p $(joinpath(outputdir, "lang"))`)

    @info "creating the word level lexicon"
    lexicon = Dict()
    phones = Set()
    for path in glob("*mb", joinpath(inputdir, "full_corpus_newsplit", "all"))
        text = lowercase(readline(path))
        for w in split(text)
            for l in w push!(phones, l) end
            lexicon[w] = [l for l in w]
        end
    end
    open(joinpath(outputdir, "lang", "lexicon.wrd"), "w") do f
        for w in sort([k for k in keys(lexicon)])
            println(f, "$w $(join(lexicon[w], " "))")
        end
    end

    @info "creating the phone level lexicon (identity map)"
    open(joinpath(outputdir, "lang", "lexicon.char"), "w") do f
        for w in sort([k for k in phones])
            println(f, "$w $w")
        end
    end

    @info "creating the unit inventory"
    open(joinpath(outputdir, "lang", "units"), "w") do f
        println(f, "sil non-speech-unit")
        for p in phones
            println(f, "$p speech-unit")
        end
    end

end

main()

max_stars_repo_name: BUTSpeechFIT/ConvHMM

max_issues_repo_name: BUTSpeechFIT/ConvHMM

max_forks_repo_name: BUTSpeechFIT/ConvHMM


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/functions/nonlinear/quadratic.jl
# Defines buffered structures to store quadratic functions:
# BufferedQuadraticIneq, BufferedQuadraticEq, as well as the
# lower_interval_bound, interval_bound, and eliminate_fixed_variables!
# functions associated with each structure.
#############################################################################

###
### Structure definitions
###

"""
$(TYPEDEF)

Stores a general quadratic inequality constraint with a buffer.
"""
mutable struct BufferedQuadraticIneq{T<:AbstractFloat} <: AbstractEAGOConstraint
    func::SQF{T}
    buffer::Dict{Int,T}
    saf::SAF{T}
    len::Int
end

"""
$(TYPEDEF)

Stores a general quadratic equality constraint with a buffer.
"""
mutable struct BufferedQuadraticEq{T<:AbstractFloat} <: AbstractEAGOConstraint
    func::SQF{T}
    minus_func::SQF{T}
    buffer::Dict{Int,T}
    saf::SAF{T}
    len::Int
end

#=
mutable struct BufferedConvexQuadratic <: AbstractEAGOConstraint
    func::SQF
    buffer::Dict{Int, Float64}
    saf::SAF
    len::Int
end
=#

###
### Constructor definitions
###

function create_buffer_dict(func::SQF{T}) where {T<:AbstractFloat}

    d = Dict{Int,T}()

    for term in func.quadratic_terms
        d[term.variable_index_1.value] = zero(T)
        d[term.variable_index_2.value] = zero(T)
    end

    for term in func.affine_terms
        d[term.variable_index.value] = zero(T)
    end

    return d
end

function BufferedQuadraticIneq{T}() where T<:AbstractFloat
    BufferedQuadraticIneq{T}(SQF{T}(SQT{T}[], SAT{T}[], zero(T)),
                             Dict{Int, T}(), SAF{T}(SAT{T}[], zero(T)), zero(T))
end

function BufferedQuadraticIneq(func::SQF{T}, set::LT{T}) where T<:AbstractFloat
    b = create_buffer_dict(func)
    saf = SAF{T}([SAT{T}(zero(T), VI(k)) for k in keys(b)], zero(T))
    cfunc = copy(func)
    cfunc.constant -= set.upper
    return BufferedQuadraticIneq{T}(cfunc, b, saf, length(b))
end

function BufferedQuadraticIneq(func::SQF{T}, set::GT{T}) where T<:AbstractFloat
    b = create_buffer_dict(func)
    saf = SAF{T}([SAT{T}(zero(T), VI(k)) for k in keys(b)], zero(T))
    cfunc = MOIU.operate(-, T, func)
    cfunc.constant += set.lower
    return BufferedQuadraticIneq(cfunc, b, saf, length(b))
end

function BufferedQuadraticEq{T}() where T<:AbstractFloat
    BufferedQuadraticEq{T}(SQF{T}(SQT{T}[], SAT{T}[], zero(T)),
                           SQF{T}(SQT{T}[], SAT{T}[], zero(T)),
                           Dict{Int,T}(), SAF{T}(SAT{T}[], zero(T)), zero(T))
end

function BufferedQuadraticEq(func::SQF{T}, set::ET{T}) where T<:AbstractFloat
    b = create_buffer_dict(func)
    saf = SAF{T}([SAT{T}(zero(T), VI(k)) for k in keys(b)], zero(T))
    cfunc1 = copy(func)
    cfunc1.constant -= set.value
    cfunc2 = MOIU.operate(-, T, func)
    cfunc2.constant += set.value
    return BufferedQuadraticEq(cfunc1, cfunc2, buffer, saf, length(b))
end

#=
function BufferedConvexQuadratic(f::BufferedQuadraticIneq)
    BufferedConvexQuadratic(copy(f.func), copy(f.buffer), copy(f.saf), f.len)
end
=#

###
### Parsing definitions
###

function eliminate_fixed_variables!(f::T, v::Vector{VariableInfo}) where T <: Union{BufferedQuadraticIneq,
                                                                                BufferedQuadraticIneq}
    deleted_count = 0
    index = 1
    while i + deleted_count <= f.len
        term = @inbounds f.sqf.terms[i]
        variable_info_1 = @inbounds v[term.variable_index_1.value]
        variable_info_2 = @inbounds v[term.variable_index_2.value]
        if variable_info_1.is_fixed && variable_index_2.is_fixed
            f.sqf.constant += coeff*variable_info_1.lower_bound*variable_index_2.lower_bound
            deleteat!(f.sqf.terms, i)
            deleted_count += 1
        else
            i += 1
        end
    end
    f.len -= deleted_count

    return nothing
end


text:
#= MATLAB Code from James Cai's PGEToolbox
%TAJIMA89D - calculates Tajima's D directly
%
% Syntax: [d,theta,pval] = tajima89d(nsam, Sn, sumk)
%
% nsam   - Number of smaples
% Sn     - Number of segregating sites or total number of mutations
% sumk   - Average number of nucleotide differences
% d      - tajima's D value
% theta  - estimated theta (4Neu) from Sn
% pval   - P-value for a double sided test assuming a normal
%          distribution for D

% Population Genetics and Evolution Toolbox (PGEToolbox)
% Author: James Cai
% Email: jcai@tamu.edu
%
% $LastChangedDate: 2013-01-06 13:39:38 -0600 (Sun, 06 Jan 2013) $
% $LastChangedRevision: 331 $
% $LastChangedBy: jcai $

if (Sn==0), d=0; theta=0; pval=1; return; end

nx=1:(nsam-1);
a1 = sum(1./nx);
a2 = sum(1./nx.^2);
b1=(nsam+1)/(3*(nsam-1));
b2=2*(nsam*nsam+nsam+3)/(9*nsam*(nsam-1));
c1=b1-1/a1;
c2=b2-(nsam+2)/(a1*nsam)+a2/(a1^2);
e1=c1/a1;
e2=c2/(a1^2+a2);

d=(sumk-Sn/a1)/sqrt(e1*Sn+e2*Sn*(Sn-1));

if nargout>1
    theta=Sn/a1;
    % P-value for a double sided test assuming a normal distribution for D
    pval = 2*normcdf(-abs(d));
    % [above] This evaluates the tail probability directly, avoiding
    % numerical cancellation issues in the upper tail.
    % better than the following:
    % p = normcdf(z,0,1);
    % p = 2*min(p,1-p);
    % SEE: Matlab Solution ID - 1-1YDCNJ
end
=#

## the julia port ##


text:
using HorizonSideRobots
include("../basics.jl")

md"""
            task_9 c перегородками
"""

"""Двигает робота на {x} шагов или останавливается на маркере"""
function super_move_for_x_or_find_marker!(r, side, x)
    for _ in 1:x
        super_move!(r, side)
        if ismarker(r)
            return false
        end
    end
    return true
end

"""Двигает робота по спирали пока он не найдет маркер"""
function super_radial_search!(r)
    x = 1
    count = 0
    side = Nord

    while super_move_for_x_or_find_marker!(r, side, x)
        count += 1
        side = next_side(side)
        if count % 2 == 0
            x += 1
        end
    end

end

 """Решение задачи 19"""
 function solve_19!(r)
    super_radial_search!(r)

    show(r)
end


text:
@testset "probabilityweightedmoment_gumbel.jl" begin
    n = 5000
    θ = [0.0 ; 0.0]

    pd = Gumbel(θ[1], exp(θ[2]))
    y = rand(pd, n)

    @testset "gumbelfitpwm(y)" begin
        # stationary model building
        fm = Extremes.gumbelfitpwm(y)

        cinterval = cint(fm)

        @test [x[1] for x in cinterval] <= θ
        @test θ <= [x[2] for x in cinterval]

    end

    @testset "gumbelfitpwm(df, datacol)" begin
        # stationary model building
        df = DataFrame(y = y)
        fm = Extremes.gumbelfitpwm(df, :y)

        cinterval = cint(fm)

        @test [x[1] for x in cinterval] <= θ
        @test θ <= [x[2] for x in cinterval]

    end

    @testset "gumbelfitpwm(model)" begin
        # non-stationary warn
        model = Extremes.BlockMaxima(Variable("y", y), locationcov = [Variable("t", collect(1:n))])

        @test_logs (:warn, "covariates cannot be included in the model when estimating the
            paramters by the probability weighted moment parameter estimation.
            The estimates for the stationary model is returned.") Extremes.gumbelfitpwm(model)

        # stationary Gumbel fit by pwm
        model = Extremes.BlockMaxima(Variable("y", y))

        fm = Extremes.gumbelfitpwm(model)

        cinterval = cint(fm)

        @test [x[1] for x in cinterval] <= θ
        @test θ <= [x[2] for x in cinterval]

    end
end