text:
module ContextualMacros

export @contextualmacro

using Base.Meta: isexpr

struct ContextualMacro{macroname} <: Function end

"""
    ContextualMacros.def(macroname::Symbol) -> f

Return a function callable as `f(__source__, __module__, args...)`.
"""
def(macroname::Symbol) = ContextualMacro{macroname}()

function _macroname_of(ex)
    if !(isexpr(ex, :macrocall, 2) && ex.args[1] isa Symbol)
        throw(ArgumentError("Expected a macro name (e.g., `@f`). Got:\n$ex"))
    end
    return Symbol(String(ex.args[1])[2:end])
end

"""
    @contextualmacro @macroname

Define a contextual macro `@macroname`.
"""
macro contextualmacro(ex)
    macroname = _macroname_of(ex)
    quote
        macro $macroname(args...)
            $(def(macroname))(__source__, __module__, args...)
        end
    end |> esc
end

"""
    ContextualMacros.@def @macroname

Define a contextual macro `@macroname` whose identity is solely
determined by the name of the macro.

!!! warning "Unstable"

    This is an experimental API.
"""
macro def(ex)
    macroname = _macroname_of(ex)
    withat = Symbol("@$macroname")
    esc(:(const $withat = $(QuoteNode(def(macroname)))))
end
# Experimental because the calling convention of macros is not a part
# of the stable API of Julia.

"""
    ContextualMacros.expandwith(module, expression; definitions...)

Macro-expand `expression` with contextual macro `@f` whose definition
is given by `definitions[:f]`.

Each value of `definitions` is a function that takes a single argument
with a "context" object with (at least) following properties:

- `__module__::Module`: module in which the macro is expanded
- `__source__::LineNumberNode`: line at which the macro is expanded
- `args`: arguments to the macro
"""
function expandwith(__module__::Module, ex; kwargs...)
    with(; kwargs...) do
        macroexpand(__module__, ex)
    end
end

"""
    ContextualMacros.with(f; definitions...)

Run function `f` while using the definitions of the macros specified
by `definitions`.  See also [`expandwith`](@ref).
"""
function with(f; definitions...)
    lock(DEFINITION_LOCK) do
        original = copy(DEFINITION_STACK)
        try
            push!(DEFINITION_STACK, values(definitions))
            f()
        finally
            expected = vcat(original, [values(definitions)])
            if !isequal(DEFINITION_STACK, expected)
                @error("DEFINITION_STACK is corrupted. Trying to recover.")
            end
            append!(empty!(DEFINITION_STACK), original)
        end
    end
end

# Not sure if macro expansions can run concurrently.  But let's lock
# it to play on the safe side...
const DEFINITION_LOCK = ReentrantLock()
const DEFINITION_STACK = []
const DEFINITION_CONSUMERS = Pair{String,Vector{Symbol}}[]

function (::ContextualMacro{macroname})(__source__, __module__, args...) where {macroname}
    ctx = (__module__ = __module__, __source__ = __source__, args = args)
    for def in reverse(DEFINITION_STACK)
        if (f = get(def, macroname, nothing)) !== nothing
            return f(ctx)
        end
    end
    throw(NotInContextError(macroname))
end

struct NotInContextError <: Exception
    macroname::Symbol
end

function Base.showerror(io::IO, err::NotInContextError)
    outside_context_error_message(io, err.macroname)
end

function providers_by_macroname(macroname, consumers = DEFINITION_CONSUMERS)
    providers = String[]
    for (pro, macros) in consumers
        if macroname in macros
            push!(providers, pro)
        end
    end
    return providers
end

function outside_context_error_message(io::IO, macroname)
    providers = providers_by_macroname(macroname)
    println(io, "Macro `@", macroname, "` used outside context.")
    if isempty(providers)
        print(io, "No package implements `@", macroname, "`.")
        return
    end
    print(io, "It must be used in one of the following macro(s):")
    for n in providers
        println(io)
        print(io, n)
    end
end

"""
    ContextualMacros.register(context_name::String, macronames::Symbol...)
"""
function register(context_name::String, macronames::Symbol...)
    lock(DEFINITION_LOCK) do
        push!(DEFINITION_CONSUMERS, context_name => collect(Symbol, macronames))
    end
    return
end

function _cleanup()
    empty!(DEFINITION_STACK)
    empty!(DEFINITION_CONSUMERS)
end

function __init__()
    # Make it work nicely in sysimage:
    _cleanup()
    atexit(_cleanup)
end

# Use README as the docstring of the module:
@doc let path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    replace(read(path, String), r"^```julia"m => "```jldoctest README")
end ContextualMacros

end


text:
#
# Basic definitions
#

immutable Hyper{T<:Real} <: Number
  f0::T
  f1::T
  f2::T
  f12::T
end

Hyper(x::Real, eps1::Real, eps2::Real, eps1eps2::Real) =
  Hyper(promote(x, eps1, eps2, eps1eps2)...)

Hyper(x::Real) = Hyper(x, zero(x), zero(x), zero(x))
Hyper() = Hyper(0.0, 0.0, 0.0, 0.0)

typealias Hyper256 Hyper{Float64}
Hyper256() = Hyper256(0.0, 0.0, 0.0, 0.0)
typealias Hyper128 Hyper{Float32}
Hyper128() = Hyper128(0.0, 0.0, 0.0, 0.0)

real(z::Hyper) = z.f0
eps1(z::Hyper) = z.f1
eps2(z::Hyper) = z.f2
eps1eps2(z::Hyper) = z.f12

eps(z::Hyper) = eps(real(z))
eps{T}(::Type{Hyper{T}}) = eps(T)
one(z::Hyper) = Hyper(one(real(z)))
one{T}(::Type{Hyper{T}}) = Hyper(one(T))
nan{T}(::Type{Hyper{T}}) = nan(one(T))
isnan(z::Hyper) = isnan(real(z))

convert{T<:Real}(::Type{Hyper{T}}, x::Real) =
  Hyper{T}(convert(T, x), convert(T, 0), convert(T, 0), convert(T, 0))

convert{T<:Real}(::Type{Hyper{T}}, z::Hyper{T}) = z

convert{T<:Real}(::Type{Hyper{T}}, z::Hyper) =
  Hyper{T}(convert(T, real(z)), convert(T, eps1(z)), convert(T, eps2(z)),
    convert(T, eps1eps2(z)))

convert{T<:Real}(::Type{T}, z::Hyper) =
  ((eps1(z) == 0 && eps2(z) == 0 && eps1eps2(z)) ? convert(T, real(z)) : throw(InexactError()))

promote_rule{T<:Real, S<:Real, Q<:Real, P<:Real}(::Type{Hyper{T}}, ::Type{Hyper{S}}, ::Type{Hyper{Q}}, ::Type{Hyper{P}}) =
    Hyper{promote_type(T, S, Q, P)}

promote_rule{T<:Real}(::Type{Hyper{T}}, ::Type{T}, ::Type{T}, ::Type{T}) = Hyper{T}

promote_rule{T<:Real, S<:Real, Q<:Real, P<:Real}(::Type{Hyper{T}}, ::Type{S}, ::Type{Q}, ::Type{P}) =
    Hyper{promote_type(T, S, Q, P)}

hyper(x, y, z, yz) = Hyper(x, y, z, yz)
hyper(x) = Hyper(x)
hyper() = Hyper()

@vectorize_1arg Real hyper

hyper256(s::Float64, t::Float64, q::Float64, p::Float64) = Hyper{Float64}(s, t, q, p)
hyper256(s::Real, t::Real, q::Real, p::Real) =
  hyper256(float64(s), float64(t), float64(q), float64(p))
hyper256(z) = hyper256(real(z), eps1(z), eps2(z), eps1eps2(z))

hyper128(s::Float32, t::Float32, q::Float32, p::Float32) = Hyper{Float32}(s, t, q, p)
hyper128(s::Real, t::Real, q::Real, p::Real) =
  hyper128(float32(s), float32(t), float32(q), float32(p))
hyper128(z) = hyper128(real(z), eps1(z), eps2(z), eps1eps2(z))

ishyper(x::Hyper) = true
ishyper(x::Number) = false

real_valued{T<:Real}(z::Hyper{T}) = (eps1(z) == 0 && eps2(z) == 0 && eps1eps2(z) == 0)
integer_valued(z::Hyper) = real_valued(z) && integer_valued(real(z))

isfinite(z::Hyper) = isfinite(real(z))
reim(z::Hyper) = (real(z), eps1(z), eps2(z), eps1eps2(z))

# No-op?
conjhyper(z::Hyper) = z

#
# IO definitions
#

function hyper_show(io::IO, z::Hyper, compact::Bool)
  f0, f1, f2, f12 = reim(z)
  if isnan(f0) || (isfinite(f1) || isfinite(f2) || isfinite(f12))
    compact ? showcompact(io, f0) : show(io, f0)
    if isfinite(f1)
      if signbit(f1) == 1 && !isnan(f1)
        f1 = -f1
        print(io, compact ? "-" : " - ")
      else
        print(io, compact ? "+" : " + ")
      end
      compact ? showcompact(io, f1) : show(io, f1)
      if !(isa(f1, Integer) || isa(f1, Rational) ||
           isa(f1, FloatingPoint) || isfinite(f1))
           print(io, "*")
      end
      print(io, "\u03F51")
    end
    if isfinite(f2)
      if signbit(f2) == 1 && !isnan(f2)
        f2 = -f2
        print(io, compact ? "-" : " - ")
      else
        print(io, compact ? "+" : " + ")
      end
      compact ? showcompact(io, f2) : show(io, f2)
      if !(isa(f2, Integer) || isa(f2, Rational) ||
           isa(f2, FloatingPoint) || isfinite(f2))
           print(io, "*")
      end
      print(io, "\u03F52")
    end
    if isfinite(f12)
      if signbit(f12) == 1 && !isnan(f12)
        f12 = -f12
        print(io, compact ? "-" : " - ")
      else
        print(io, compact ? "+" : " + ")
      end
      compact ? showcompact(io, f12) : show(io, f12)
      if !(isa(f12, Integer) || isa(f12, Rational) ||
           isa(f12, FloatingPoint) || isfinite(f12))
           print(io, "*")
      end
      print(io, "\u03F51\u03F52")
    end
  else
    print(io, "Hyper(", f0, ",", f1, ",", f2, ",", f12, ")")
  end
end

show(io::IO, z::Hyper) = hyper_show(io, z, false)
showcompact(io::IO, z::Hyper) = hyper_show(io, z, true)

function read{T<:Real}(s::IO, ::Type{Hyper{T}})
  f0 = read(s, T)
  f1 = read(s, T)
  f2 = read(s, T)
  f12 = read(s, T)
  Hyper{T}(f0, f1, f2, f12)
end

function write(s::IO, z::Hyper)
  write(s, real(z))
  write(s, eps1(z))
  write(s, eps2(z))
  write(s, eps1eps2(z))
end

#
# Generic function on hyperdualnumbers
#

convert(::Type{Hyper}, z::Hyper) = z
convert(::Type{Hyper}, x::Real) = hyper(x)

+(z::Hyper, w::Hyper) = hyper(real(z) + real(w), eps1(z) + eps1(w),
  eps2(z) + eps2(w), eps1eps2(z) + eps1eps2(w))
+(z::Number, w::Hyper) = hyper(z + real(w), eps1(w), eps2(w), eps1eps2(w))
+(z::Hyper, w::Number) = hyper(real(z) + w, eps1(z), eps2(z), eps1eps2(z))

-(z::Hyper) = hyper(-real(z), -eps1(z), -eps2(z), -eps1eps2(z))
-(z::Hyper, w::Hyper) = hyper(real(z) - real(w), eps1(z) - eps1(w),
  eps2(z) - eps2(w), eps1eps2(z) - eps1eps2(w))
-(z::Number, w::Hyper) = hyper(z - real(w), -eps1(w), -eps2(w), -eps1eps2(w))
-(z::Hyper, w::Number) = hyper(real(z) - w, eps1(z), eps2(z), eps1eps2(z))

*(x::Bool, z::Hyper) = ifelse(x, z, ifelse(signbit(real(z)) == 0, zero(z), -zero(z)))
*(x::Hyper, z::Bool) = z * x

*(z::Hyper, w::Hyper) = hyper(real(z) * real(w),
  real(z)*eps1(w)+eps1(z)*real(w), real(z)*eps2(w)+eps2(z)*real(w),
  real(z)*eps1eps2(w)+eps1(z)*eps2(w)+eps2(z)*eps1(w)+real(w)*eps1eps2(z))

*(z::Hyper, w::Real) = hyper(real(z) * w, eps1(z)*w, eps2(z)*w, w*eps1eps2(z))
*(z::Number, w::Hyper) = w * z

/(z::Hyper, w::Hyper) = z*(one(real(z))/w)
function /(z::Number, w::Hyper)
    invrw = one(z)/real(w)
    deriv = -z*invrw^2
    hyper(z*invrw, eps1(w)*deriv, eps2(w)*deriv,
          eps1eps2(w)*deriv-2eps1(w)*eps2(w)*deriv*invrw)
end

# Needed to prevent ambiguous warning:
#   /(Number,Complex{T<:Real}) at complex.jl:127

/(z::Hyper, w::Complex) = hyper(real(z)/w, eps1(z)/w, eps2(z)/w, eps1eps2(z)/w)
/(z::Hyper, w::Number) = hyper(real(z)/w, eps1(z)/w, eps2(z)/w, eps1eps2(z)/w)

function ^(z::Hyper, w::Rational)
  deriv = w * real(z)^(w-1)
  hyper(real(z)^w, eps1(z)*deriv, eps2(z)*deriv,
    eps1eps2(z)*deriv+w*(w-1)*eps1(z)*eps2(z)*real(z)^(w-2))
end

function ^(z::Hyper, w::Integer)
  deriv = w * real(z)^(w-1)
  hyper(real(z)^w, eps1(z)*deriv, eps2(z)*deriv,
    eps1eps2(z)*deriv+w*(w-1)*eps1(z)*eps2(z)*real(z)^(w-2))
end

function ^(z::Hyper, w::Number)
  xval = real(z)
  tol = typeof(xval)==FloatingPoint ? eps(xval) : 10.0^-15
  if abs(float(xval)) < tol
    xval = ifelse(signbit(xval)==0, tol, -tol)
  end
  deriv = w * xval^(w-1)
  # Use actual value for f0, tol for deriv only
  hyper(real(z)^w, eps1(z)*deriv, eps2(z)*deriv,
    eps1eps2(z)*deriv+w*(w-1)*eps1(z)*eps2(z)*xval^(w-2))
end

^(z::Hyper, w::Hyper) = exp(w*log(z))

function exp(z::Hyper)
  deriv = exp(real(z))
  hyper(deriv, deriv*eps1(z), deriv*eps2(z), deriv*(eps1eps2(z)+eps1(z)*eps2(z)))
end

function log(z::Hyper)
  deriv1 = eps1(z)/real(z)
  deriv2 = eps2(z)/real(z)
  hyper(log(real(z)), deriv1, deriv2, eps1eps2(z)/real(z)-(deriv1*deriv2))
end

function sin(z::Hyper)
  funval = sin(real(z))
  deriv = cos(real(z))
  hyper(funval, deriv*eps1(z),deriv*eps2(z),deriv*eps1eps2(z)-funval*eps1(z)*eps2(z))
end

function cos(z::Hyper)
  funval = cos(real(z))
  deriv = -sin(real(z))
  hyper(funval, deriv*eps1(z),deriv*eps2(z),deriv*eps1eps2(z)-funval*eps1(z)*eps2(z))
end

function tan(z::Hyper)
  funval = tan(real(z))
  deriv = funval*funval+1
  hyper(funval, deriv*eps1(z),deriv*eps2(z),deriv*eps1eps2(z)+eps1(z)*eps2(z)*(2*funval*deriv))
end

function asin(z::Hyper)
  funval = asin(real(z))
  deriv1 = 1.0-real(z)*real(z)
  deriv = 1.0/sqrt(deriv1)
  hyper(funval, deriv*eps1(z),deriv*eps2(z),deriv*eps1eps2(z)+eps1(z)*eps2(z)*(real(z)/deriv1^-1.5))
end

function acos(z::Hyper)
  funval = acos(real(z))
  deriv1 = 1.0-real(z)*real(z)
  deriv = -1.0/sqrt(deriv1)
  hyper(funval, deriv*eps1(z),deriv*eps2(z),deriv*eps1eps2(z)+eps1(z)*eps2(z)*(-real(z)/deriv1^-1.5))
end

function atan(z::Hyper)
  funval = atan(real(z))
  deriv1 = 1.0+real(z)*real(z)
  deriv = 1.0/sqrt(deriv1)
  hyper(funval, deriv*eps1(z),deriv*eps2(z),
    deriv*eps1eps2(z)+eps1(z)*eps2(z)*(-2*real(z)/deriv1*deriv1))
end

sqrt(z::Hyper) = z^0.5

maximum(z::Hyper, w::Hyper) = z > w ? z : w
maximum(z::Hyper, w::Number) = z > w ? z : hyper(w)
maximum(z::Number, w::Hyper) = z > w ? hyper(z) : w

minimum(z::Hyper, w::Hyper) = z < w ? z : w
minimum(z::Hyper, w::Number) = z < w ? z : hyper(w)
minimum(z::Number, w::Hyper) = z < w ? hyper(z) : w

>(z::Hyper, w::Hyper) = real(z) > real(w)
>(z::Hyper, w::Number) = real(z) > w
>(z::Number, w::Hyper) = z > real(w)
>=(z::Hyper, w::Hyper) = real(z) >= real(w)
>=(z::Hyper, w::Number) = real(z) >= w
>=(z::Number, w::Hyper) = z >= real(w)
<(z::Hyper, w::Hyper) = real(z) < real(w)
<(z::Hyper, w::Number) = real(z) < w
<(z::Number, w::Hyper) = z < real(w)
<=(z::Hyper, w::Hyper) = real(z) <= real(w)
<=(z::Hyper, w::Number) = real(z) <= w
<=(z::Number, w::Hyper) = z <= real(w)

==(z::Hyper, w::Hyper) = real(z) == real(w) && eps1(z) == eps1(w) &&
  eps2(z) == eps2(w) && eps1eps2(z) == eps1eps2(w)

==(z::Hyper, x::Real) = real_valued(z) && real(z) == x
==(x::Real, z::Hyper) = ==(z, x)

isequal(z::Hyper, w::Hyper) = isequal(real(z), real(w)) && isequal(eps1(z), eps1(w)) &&
  isequal(eps2(z), eps2(w)) && isequal(eps1eps2(z), eps1eps2(w))

isequal(z::Hyper, x::Real) =
isequal(x::Real, z::Hyper) = ==(z, x)

ext: jl

max_stars_repo_name: hwborchers/RDualNumbers

max_issues_repo_name: hwborchers/RDualNumbers

max_forks_repo_name: hwborchers/RDualNumbers


text:
using CSV, Distributions, HypothesisTests

data1 = CSV.read("../data/machine1.csv", header=false)[:,1]
data2 = CSV.read("../data/machine2.csv", header=false)[:,1]
xBar1, xBar2 = mean(data1), mean(data2)
n1, n2 = length(data1), length(data2)
sig1, sig2 = 1.2, 1.6
alpha = 0.05
z = quantile(Normal(),1-alpha/2)

println("Calculating formula: ", (xBar1 - xBar2 - z*sqrt(sig1^2/n1+sig2^2/n2),
				  xBar1 - xBar2 + z*sqrt(sig1^2/n1+sig2^2/n2)))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
include("nearly_zero.jl")
include("smaller_vector.jl")


text:
# ---
# title: 611. Valid Triangle Number
# id: problem611
# author: Indigo
# date: 2021-06-26
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/valid-triangle-number/description/>
# hidden: true
# ---
#
# Given an array consists of non-negative integers, your task is to count the
# number of triplets chosen from the array that can make triangles if we take
# them as side lengths of a triangle.
#
# **Example 1:**
#
#
#
#     Input: [2,2,3,4]
#     Output: 3
#     Explanation:
#     Valid combinations are:
#     2,3,4 (using the first 2)
#     2,3,4 (using the second 2)
#     2,2,3
#
#
# **Note:**
#
#   1. The length of the given array won't exceed 1000.
#   2. The integers in the given array are in the range of [0, 1000].
#
#
## @lc code=start
using LeetCode

function triangle_number(nums::Vector{Int})
    res, len = 0, length(nums)
    sort!(nums)
    @inbounds for i in 1:len-2
        nums[i] == 0 && continue
        k = i + 2
        for j in i+1:len-1
            while k <= len && nums[i] + nums[j] > nums[k]; k += 1; end
            res += k - j - 1
        end
    end
    return res
end
## @lc code=end


text:
using Test
using BenchmarkTools
using Plots
using SparseArrays

model,obj = Dynamics.quadrotor
model,obj_con = Dynamics.quadrotor_3obs

##########
N = 50 # 201
integration = :rk4
opts = SolverOptions()
opts.verbose = false
opts.square_root = false
opts.cost_tolerance = 1e-5
opts.cost_tolerance_intermediate = 1e-5
opts.constraint_tolerance = 1e-4
opts.outer_loop_update_type = :feedback

# Obstacle Avoidance
model,obj_uncon = TrajectoryOptimization.Dynamics.quadrotor
r_quad = 3.0
n = model.n
m = model.m
obj_con = TrajectoryOptimization.Dynamics.quadrotor_3obs[2]
spheres = TrajectoryOptimization.Dynamics.quadrotor_3obs[3]
n_spheres = length(spheres)

solver_uncon = Solver(model,obj_uncon,integration=integration,N=N,opts=opts)
solver_con = Solver(model,obj_con,integration=integration,N=N,opts=opts)

U_hover = 0.5*9.81/4.0*ones(solver_uncon.model.m, solver_uncon.N-1)
X_hover = rollout(solver_uncon,U_hover)
# @time results_uncon, stats_uncon = solve(solver_uncon,U_hover)
# @time results_uncon_dircol, stats_uncon_dircol = TrajectoryOptimization.solve_dircol(solver_uncon, X_hover, U_hover, options=dircol_options)

@time res, stats_con = solve(solver_con,U_hover)
solver = solver_con
#########

n = 13
m = 4

# Unconstrained
Q = (1e-3)*Matrix(I,n,n)
R = (1e-2)*Matrix(I,m,m)
# R = (1e-1)*Matrix(I,m,m)
Qf = (1000.0)*Matrix(I,n,n)
tf = 1.0
dt = 0.1

# -initial state
x0 = zeros(n)
x0[1:3] = [0.; 0.; 0.]
q0 = [1.;0.;0.;0.]
x0[4:7] = q0

# -final state
xf = copy(x0)
xf[1:3] = [0.;10.;0.] # xyz position
xf[4:7] = q0
u_max = 10.0
u_min = 0.0

obj_uncon = LQRObjective(Q, R, Qf, tf, x0, xf)
obj_con = TrajectoryOptimization.ConstrainedObjective(obj_uncon,u_min=u_min,u_max=u_max)#,u_min=u_min)
solver = Solver(model,obj_con,N=101)

solver.obj
solver.opts.cost_tolerance = 1e-5
solver.opts.constraint_tolerance = 1e-4
solver.opts.penalty_max = 1e8
n,m,N = get_sizes(solver)
Random.seed!(1)
U = 0.5*9.81/4.0*ones(solver.model.m, solver.N-1)
res,stats = solve(solver,U)
plot(res.U)
plot(res.X)

cost(solver,res)
λ_update_default!(res,solver);
update_constraints!(res,solver)
max_violation(res)
@assert max_violation(res) < opts.constraint_tolerance

J_prev = cost(solver,res)

p,pI,pE = get_num_constraints(solver)
p_N,pI_N,pE_N = get_num_terminal_constraints(solver)

## Newton 2 ###############
results_new = copy(res)
newton_results = NewtonResults(solver)
newton_active_set!(newton_results,results_new,solver)
update_newton_results!(newton_results,results_new,solver)
newton_step!(results_new,newton_results,solver,1.0)
max_violation(results_new)
newton_cost(res,newton_results,solver)

# newton_solve!(res,solver)


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

mutable struct CoreDatabaseManager <: AbstractCoreDatabaseManager

end

function get_database(client::AbstractCoreClient, name::AbstractString)::CoreDatabase
    isempty(name) && throw(TypeDBClientException(CLIENT_MISSING_DB_NAME))
    if contains_database(client, name)
        return CoreDatabase(name)
    else
        throw(TypeDBClientException(CLIENT_DB_DOES_NOT_EXIST, name))
    end
end

function contains_database(client::AbstractCoreClient, name::AbstractString)
    isempty(name) && throw(TypeDBClientException(CLIENT_MISSING_DB_NAME))
    db = DatabaseManagerRequestBuilder
    req_result, status = Proto.databases_contains(client.core_stub.blockingStub, gRPCController() , db.contains_req(name))
    return grpc_result_or_error(req_result, status, result->result.contains)
end

function get_all_databases(client::AbstractCoreClient)::Vector{CoreDatabase}
    db = DatabaseManagerRequestBuilder
    req_result, status = Proto.databases_all(client.core_stub.blockingStub, gRPCController(), db.all_req())
    return grpc_result_or_error(req_result, status, result -> [CoreDatabase(db_name) for db_name in result.names])
end

function create_database(client::AbstractCoreClient, name::AbstractString)
    isempty(name) && throw(TypeDBClientException(CLIENT_MISSING_DB_NAME))
    db = DatabaseManagerRequestBuilder
    req_result, status =  Proto.databases_create(client.core_stub.blockingStub, gRPCController(), db.create_req(name))
    return grpc_result_or_error(req_result, status, result->true)
end

"""
    delete_database(client::AbstractCoreClient, name::AbstractString)

Delete the database for the given name without any question. Be carful. There is no
recovery function within this call.
"""
function delete_database(client::AbstractCoreClient, name::AbstractString)
    database = get_database(client, name)
    db = DatabaseRequestBuilder
    req_result, status =  Proto.database_delete(client.core_stub.blockingStub, gRPCController(), db.delete_req(database.name))
    return grpc_result_or_error(req_result, status, result->true)
end


text:
module SphericalHarmonics

using StaticArrays, LinearAlgebra

import ACE, ACEbase

import ACE: valtype, gradtype,
		      evaluate!, evaluate_d!, evaluate_ed!,
			   write_dict, read_dict,
				ACEBasis,
				acquire_B!, release_B!,
				acquire_dB!, release_dB!,
				acquire!, release!

import ACE: VectorPool, ArrayCache

export SHBasis

# --------------------------------------------------------
#     Coordinates
# --------------------------------------------------------

"""
`struct SphericalCoords` : a simple datatype storing spherical coordinates
of a point (x,y,z) in the format `(r, cosφ, sinφ, cosθ, sinθ)`.

Use `spher2cart` and `cart2spher` to convert between cartesian and spherical
coordinates.
"""
struct SphericalCoords{T}
	r::T
	cosφ::T
	sinφ::T
	cosθ::T
	sinθ::T
end

spher2cart(S::SphericalCoords) = S.r * SVector(S.cosφ*S.sinθ, S.sinφ*S.sinθ, S.cosθ)

function cart2spher(R::AbstractVector)
	@assert length(R) == 3
	r = norm(R)
	φ = atan(R[2], R[1])
	sinφ, cosφ = sincos(φ)
	cosθ = R[3] / r
	sinθ = sqrt(R[1]^2+R[2]^2) / r
	return SphericalCoords(r, cosφ, sinφ, cosθ, sinθ)
end

SphericalCoords(φ, θ) = SphericalCoords(1.0, cos(φ), sin(φ), cos(θ), sin(θ))
SphericalCoords(r, φ, θ) = SphericalCoords(r, cos(φ), sin(φ), cos(θ), sin(θ))

"""
convert a gradient with respect to spherical coordinates to a gradient
with respect to cartesian coordinates
"""
function dspher_to_dcart(S, f_φ_div_sinθ, f_θ)
	r = S.r + eps()
   return SVector( - (S.sinφ * f_φ_div_sinθ) + (S.cosφ * S.cosθ * f_θ),
			            (S.cosφ * f_φ_div_sinθ) + (S.sinφ * S.cosθ * f_θ),
			 			                                 - (   S.sinθ * f_θ) ) / r
end

# --------------------------------------------------------
#     Indexing
# --------------------------------------------------------

"""
`sizeP(maxL):`
Return the size of the set of Associated Legendre Polynomials ``P_l^m(x)`` of
degree less than or equal to the given maximum degree
"""
sizeP(maxL) = div((maxL + 1) * (maxL + 2), 2)

"""
`sizeY(maxL):`
Return the size of the set of real spherical harmonics ``Y_{l,m}(θ,φ)`` of
degree less than or equal to the given maximum degree
"""
sizeY(maxL) = (maxL + 1) * (maxL + 1)

"""
`index_p(l,m):`
Return the index into a flat array of Associated Legendre Polynomials `P_l^m`
for the given indices `(l,m)`.
`P_l^m` are stored in l-major order i.e. `[P(0,0), [P(1,0), P(1,1), P(2,0), ...]``
"""
index_p(l::Integer,m::Integer) = m + div(l*(l+1), 2) + 1

"""
`index_y(l,m):`
Return the index into a flat array of real spherical harmonics `Y_lm`
for the given indices `(l,m)`.
`Y_lm` are stored in l-major order i.e.
[Y(0,0), Y(1,-1), Y(1,0), Y(1,1), Y(2,-2), ...]
"""
index_y(l::Integer, m::Integer) = m + l + (l*l) + 1

function idx2lm(i::Integer)
	l = floor(Int, sqrt(i-1) + 1e-10)
	m = i - (l + (l*l) + 1)
	return l, m
end

# --------------------------------------------------------
#     Associated Legendre Polynomials
# --------------------------------------------------------

"""
`ALPolynomials` : an auxiliary datastructure for
evaluating the associated lagrange functions
used for the spherical harmonics

Constructor:
```julia
ALPolynomials(maxL::Integer, T::Type=Float64)
```
"""
struct ALPolynomials{T} <: ACEBasis
	L::Int
	A::Vector{T}
	B::Vector{T}
	B_pool::ArrayCache{T}
end

ALPolynomials(L::Integer, A::Vector{T}, B::Vector{T}) where {T}  =
		ALPolynomials(L, A, B, ArrayCache{T}())

Base.length(alp::ALPolynomials) = sizeP(alp.L)

import Base.==
==(B1::ALPolynomials{T}, B2::ALPolynomials{T}) where {T} =
		((B1.L == B2.L) && (B1.A ≈ B2.A) && (B1.B ≈ B2.B))

valtype(alp::ALPolynomials{T}, x::SphericalCoords{S}) where {T, S} =
			promote_type(T, S)

gradtype(alp::ALPolynomials{T}, x::SphericalCoords{S}) where {T, S} =
			promote_type(T, S)

function ALPolynomials(L::Integer, T::Type=Float64)
	# Precompute coefficients ``a_l^m`` and ``b_l^m`` for all l <= L, m <= l
	alp = ALPolynomials(L, zeros(T, sizeP(L)), zeros(T, sizeP(L)))
	for l in 2:L
		ls = l*l
		lm1s = (l-1) * (l-1)
		for m in 0:(l-2)
			ms = m * m
			alp.A[index_p(l, m)] = sqrt((4 * ls - 1.0) / (ls - ms))
			alp.B[index_p(l, m)] = -sqrt((lm1s - ms) / (4 * lm1s - 1.0))
		end
	end
	return alp
end

function evaluate(alp::ALPolynomials, S::SphericalCoords)
	P = acquire!(alp.B_pool, length(alp), valtype(alp, S))
	return evaluate!(P, alp, S)
end

function evaluate!(P, alp::ALPolynomials, S::SphericalCoords)
	L = alp.L
	A = alp.A
	B = alp.B
	@assert length(A) >= sizeP(L)
	@assert length(B) >= sizeP(L)
	@assert length(P) >= sizeP(L)

	temp = sqrt(0.5/π)
	P[index_p(0, 0)] = temp
	if L == 0; return P; end

	P[index_p(1, 0)] = S.cosθ * sqrt(3) * temp
	temp = - sqrt(1.5) * S.sinθ * temp
	P[index_p(1, 1)] = temp

	for l in 2:L
		il = ((l*(l+1)) ÷ 2) + 1
		ilm1 = il - l
		ilm2 = ilm1 - l + 1
		for m in 0:(l-2)
			@inbounds P[il+m] = A[il+m] * (     S.cosθ * P[ilm1+m]
  					                           + B[il+m] * P[ilm2+m] )
		end
		@inbounds P[il+l-1] = S.cosθ * sqrt(2 * (l - 1) + 3) * temp
		temp = -sqrt(1.0 + 0.5 / l) * S.sinθ * temp
		@inbounds P[il+l] = temp
	end

	return P
end

function _evaluate_ed(alp::ALPolynomials, S::SphericalCoords)
	VT = valtype(alp, S)
	P = acquire!(alp.B_pool, length(alp), VT)
	dP = acquire!(alp.B_pool, length(alp), VT)
	return _evaluate_ed!(P, dP, alp::ALPolynomials, S::SphericalCoords)
end

# this doesn't use the standard name because it doesn't
# technically perform the derivative w.r.t. S, but w.r.t. θ
# further, P doesn't store P but (P if m = 0) or (P * sinθ if m > 0)
# this is done for numerical stability
function _evaluate_ed!(P, dP, alp::ALPolynomials, S::SphericalCoords)
	L = alp.L
	A = alp.A
	B = alp.B
	@assert length(A) >= sizeP(L)
	@assert length(B) >= sizeP(L)
	@assert length(P) >= sizeP(L)
	@assert length(dP) >= sizeP(L)

	temp = sqrt(0.5/π)
	P[index_p(0, 0)] = temp
	temp_d = 0.0
	dP[index_p(0, 0)] = temp_d
	if L == 0; return P, dP; end

	P[index_p(1, 0)] = S.cosθ * sqrt(3) * temp
	dP[index_p(1, 0)] = -S.sinθ * sqrt(3) * temp + S.cosθ * sqrt(3) * temp_d
	temp1, temp_d = ( - sqrt(1.5) * temp,
					      - sqrt(1.5) * (S.cosθ * temp + S.sinθ * temp_d) )
	P[index_p(1, 1)] = temp1
	dP[index_p(1, 1)] = temp_d

	for l in 2:L
		m = 0
		@inbounds P[index_p(l, m)] =
				A[index_p(l, m)] * (     S.cosθ * P[index_p(l - 1, m)]
				             + B[index_p(l, m)] * P[index_p(l - 2, m)] )
		@inbounds dP[index_p(l, m)] =
			A[index_p(l, m)] * (
							- S.sinθ * P[index_p(l - 1, m)]
							+ S.cosθ * dP[index_p(l - 1, m)]
			             + B[index_p(l, m)] * dP[index_p(l - 2, m)] )

		for m in 1:(l-2)
			@inbounds P[index_p(l, m)] =
					A[index_p(l, m)] * (     S.cosθ * P[index_p(l - 1, m)]
					             + B[index_p(l, m)] * P[index_p(l - 2, m)] )
			@inbounds dP[index_p(l, m)] =
				A[index_p(l, m)] * (
								- S.sinθ^2 * P[index_p(l - 1, m)]
								+ S.cosθ * dP[index_p(l - 1, m)]
				             + B[index_p(l, m)] * dP[index_p(l - 2, m)] )
		end
		@inbounds P[index_p(l, l - 1)] = sqrt(2 * (l - 1) + 3) * S.cosθ * temp1
		@inbounds dP[index_p(l, l - 1)] = sqrt(2 * (l - 1) + 3) * (
									        -S.sinθ^2 * temp1 + S.cosθ * temp_d )

      (temp1, temp_d) = (
					-sqrt(1.0+0.5/l) * S.sinθ * temp1,
		         -sqrt(1.0+0.5/l) * (S.cosθ * temp1 * S.sinθ + S.sinθ * temp_d) )
		@inbounds P[index_p(l, l)] = temp1
		@inbounds dP[index_p(l, l)] = temp_d
	end

	return P, dP
end

# ------------------------------------------------------------------------
#                  Spherical Harmonics
# ------------------------------------------------------------------------

"""
`AbstractSHBasis`: This extra abstraction is no longer needed, but there uses
to be a real SH basis and in case this is revived, I am keeping it for now.
"""
abstract type AbstractSHBasis{T} <: ACEBasis end

"""
complex spherical harmonics
"""
struct SHBasis{T} <: AbstractSHBasis{T}
	alp::ALPolynomials{T}
	B_pool::ArrayCache{Complex{T}}
	dB_pool::ArrayCache{SVector{3, Complex{T}}}
end

SHBasis(maxL::Integer, T::Type=Float64) = SHBasis(ALPolynomials(maxL, T))

SHBasis(alp::ALPolynomials{T}) where {T} = SHBasis(alp,
		ArrayCache{Complex{T}}(), ArrayCache{SVector{3, Complex{T}}}())

Base.show(io::IO, SH::SHBasis) = print(io, "SHBasis(L=$(maxL(SH)))")

"""
max L degree for which the alp coefficients have been precomputed
"""
maxL(sh::AbstractSHBasis) = sh.alp.L

valtype(sh::SHBasis{T}) where {T} = Complex{T}

valtype(sh::SHBasis{T}, x::AbstractVector{S}) where {T, S} =
			Complex{promote_type(T, S)}

gradtype(sh::SHBasis{T}) where {T} = SVector{3, Complex{T}}

gradtype(sh::SHBasis{T}, x::AbstractVector{S})  where {T, S} =
			SVector{3, Complex{promote_type(T, S)}}

function ACE.degree(sh::SHBasis, i::Integer)
	l, m = idx2lm(i)
	return l
end

import Base.==
==(B1::AbstractSHBasis, B2::AbstractSHBasis) =
		(B1.alp == B2.alp) && (typeof(B1) == typeof(B2))

write_dict(SH::SHBasis{T}) where {T} =
		Dict("__id__" => "ACE_SHBasis",
			  "T" => write_dict(T),
			  "maxL" => maxL(SH))

read_dict(::Val{:ACE_SHBasis}, D::Dict) =
		SHBasis(D["maxL"], read_dict(D["T"]))

Base.length(S::AbstractSHBasis) = sizeY(maxL(S))

_evaluate_d!(dY, L, S, P, dP, ::SHBasis) = cYlm_d!(dY, L, S, P, dP)

_evaluate_ed!(Y, dY, L, S, P, dP, ::SHBasis) = cYlm_ed!(Y, dY, L, S, P, dP)

function ACE.evaluate(SH::SHBasis, R::AbstractVector)
	Y = acquire!(SH.B_pool, length(SH), valtype(SH, R))
	return evaluate!(Y, SH, R)
end

function evaluate!(Y, SH::AbstractSHBasis, R::AbstractVector)
	@assert length(R) == 3
	L = maxL(SH)
	S = cart2spher(R)
	P = evaluate(SH.alp, S)
	cYlm!(Y, maxL(SH), S, P)
	release!(P)
	return Y
end

function ACE.evaluate_d(SH::AbstractSHBasis, R::AbstractVector)
	B, dB = evaluate_ed(SH, R)
	release!(B)
	return dB
end

function ACE.evaluate_ed(SH::AbstractSHBasis, R::AbstractVector)
	Y = acquire!(SH.B_pool, length(SH), valtype(SH, R))
	dY = acquire!(SH.dB_pool, length(SH), gradtype(SH, R))
	return evaluate_ed!(Y, dY, SH, R)
end

function evaluate_ed!(Y, dY, SH::AbstractSHBasis, R::AbstractVector)
	@assert length(R) == 3
	L = maxL(SH)
	S = cart2spher(R)
	Y = acquire!(SH.B_pool, length(SH), valtype(SH, R))
	dY = acquire!(SH.dB_pool, length(SH), gradtype(SH, R))
	P, dP = _evaluate_ed(SH.alp, S)
	cYlm_ed!(Y, dY, maxL(SH), S, P, dP)
	release!(P)
	release!(dP)
	return Y, dY
end

"""
evaluate complex spherical harmonics
"""
function cYlm!(Y, L, S::SphericalCoords, P)
	@assert length(P) >= sizeP(L)
	@assert length(Y) >= sizeY(L)
   @assert abs(S.cosθ) <= 1.0

	ep = 1 / sqrt(2) + im * 0
	for l = 0:L
		Y[index_y(l, 0)] = P[index_p(l, 0)] * ep
	end

   sig = 1
   ep_fact = S.cosφ + im * S.sinφ
	for m in 1:L
		sig *= -1
		ep *= ep_fact            # ep =   exp(i *   m  * φ)
		em = sig * conj(ep)      # ep = ± exp(i * (-m) * φ)
		for l in m:L
			p = P[index_p(l,m)]
			@inbounds Y[index_y(l, -m)] = em * p   # (-1)^m * p * exp(-im*m*phi) / sqrt(2)
			@inbounds Y[index_y(l,  m)] = ep * p   #          p * exp( im*m*phi) / sqrt(2)
		end
	end

	return Y
end

"""
evaluate gradients of complex spherical harmonics
"""
function cYlm_ed!(Y, dY, L, S::SphericalCoords, P, dP)
	@assert length(P) >= sizeP(L)
	@assert length(Y) >= sizeY(L)
	@assert length(dY) >= sizeY(L)

	# m = 0 case
	ep = 1 / sqrt(2)
	for l = 0:L
		Y[index_y(l, 0)] = P[index_p(l, 0)] * ep
		dY[index_y(l, 0)] = dspher_to_dcart(S, 0.0, dP[index_p(l, 0)] * ep)
	end

   sig = 1
   ep_fact = S.cosφ + im * S.sinφ

	for m in 1:L
		sig *= -1
		ep *= ep_fact            # ep =   exp(i *   m  * φ)
		em = sig * conj(ep)      # ep = ± exp(i * (-m) * φ)
		dep_dφ = im *   m  * ep
		dem_dφ = im * (-m) * em
		for l in m:L
			p_div_sinθ = P[index_p(l,m)]
			@inbounds Y[index_y(l, -m)] = em * p_div_sinθ * S.sinθ
			@inbounds Y[index_y(l,  m)] = ep * p_div_sinθ * S.sinθ

			dp_dθ = dP[index_p(l,m)]
			@inbounds dY[index_y(l, -m)] = dspher_to_dcart(S, dem_dφ * p_div_sinθ, em * dp_dθ)
			@inbounds dY[index_y(l,  m)] = dspher_to_dcart(S, dep_dφ * p_div_sinθ, ep * dp_dθ)
		end
	end

	return Y, dY
end

end


text:
# triangle de Pascal modulo 7

mutable struct List
    head
    next
end

function get(l, i)
    if i == 0
        l.head
    else
	get(l.next, i-1)
    end
end

function set(l, i, v)
    if i == 0
        l.head = v
        return
    end
    set(l.next, i-1, v)
end

function create(n)
    if n == 0
        nothing
    else
        List(0, create(n-1))
    end
end

function print_row(r, i)
    for j = 0 : i
	if get(r, j) != 0
	    print("*")
	else
	    print(".")
	end
    end
    print("\n")
end

function compute_row(r, i)
    j = i
    while j > 0
	set(r, j, (get(r,j) + get(r,j-1)) % 7)
        j = j-1
    end
    set(r, 0, 1)
end

function pascal(n)
    r = create(n + 1)
    for i = 0 : n-1
	set(r, i, 0)
	compute_row(r, i)
	print_row(r, i)
    end
end

pascal(42)

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
## Find real zeros of a polynomial
##
## We do this two ways, depending if an exact gcd can be found (Z[x] or Q[x])
##
## For polynomials over Z[x], Q[x], we can use a modification of
## [Vincent's theorem](http://en.wikipedia.org/wiki/Vincent's_theorem)
## taken from Efficient isolation of polynomial’s real roots by
## Fabrice Rouillier; Paul Zimmermann
##

if VERSION >= v"0.5.0-"
    eval(Expr(:using, :Primes))
    eval(Expr(:import, :Primes, :factor))
else
    eval(Expr(:import, :Base, :factor))
end

###   Things for Polynomials   #####
## return "x"
## variable(p::Poly) = poly(zeros(eltype(p),1), p.var)

## make bounds work for floating point or rational.
_iszero{T <: (@compat Union{Integer, Rational})}(b::T; kwargs...) = b == 0
_iszero{T<:AbstractFloat}(b::T; xtol=1) = abs(b) <= 2*xtol*eps(T)

## return (q,r) with p(x) = (x-c)*q(x) + r using synthetic division
function Base.divrem(p::Poly, c::Number)
    ps = copy(p.a)                    # [p0, p1, ..., pn]
    qs = eltype(p)[pop!(ps)]           # take from right
    while length(ps) > 0
        unshift!(qs, c*qs[1] + pop!(ps))
    end
    r = shift!(qs)
    Poly(qs, p.var), r
end

## p(x) = q(x)*(x-c)^k for some q,k (possibly 0). Return maximal k and corresponding q
function multiplicity(p::Poly, c::Number)
    k = 0
    q, r = divrem(p,c)
    while _iszero(r)
        p = q
        q,r = divrem(p, c)
        k = k + 1
    end
    p, k
end

## Our Poly types for which we can find gcd
typealias QQ  @compat Union{Int, BigInt, Rational{Int}, Rational{BigInt}}
typealias BB  @compat Union{BigInt, Rational{BigInt}}

## Here computations are exact, as long as we return poly in Q[x]
function Base.divrem{T<:QQ, S<:QQ}(a::Poly{T}, b::Poly{S})
    degree(b) == 0 && (b[0] == 0 ? error("b is 0") : return (a/b[0], 0*b))

    lc(p::Poly) = p[degree(p)]

    a.var == b.var || DomainError() # "symbols must match"
    R = promote_type(T, S)
    x = poly(zeros(R,1), a.var)
    q, r= 0*a, a
    d,c = degree(b), lc(b)

    while degree(r) >= d
        s = lc(r)//c * x^(degree(r)- d)
        q, r = q + s, r - s*b
    end
    q,r
end

## gcd over Rational{BigInt}
function bgcd{T <: BB, S<: BB}(a::Poly{T}, b::Poly{S})
    (degree(b) == 0 && b[0] == 0) ? a : bgcd(b, divrem(a,b)[2])
end
##################################################
## Real zeros of p in Z[x] or Q[x], using VCA method
## cf. E$cient isolation of polynomial’s real roots by Fabrice Rouillier, Paul Zimmermann

## Polynomial transformations
" `R(p)` finds  `x^n p(1/x)` which is a reversal of coefficients "
R(p) = Poly(reverse(p.a), p.var)

" `p(λ x)`: scale x axis by λ "
Hλ(p, λ=1//2) = polyval(p, λ * variable(p))

" `p(x + λ)`: Translate polynomial left by λ "
Tλ(p, λ=1)   = polyval(p, variable(p) + λ)

" shift and scale p so that  [c,c+1]/2^k -> [0,1] "
function Pkc(p::Poly, k, c)
    n = degree(p)
    2^(k*n) * Hλ(Tλ(p, c/2^k), 1/2^k)
end

## Upper bound on size of real roots that is tighter than cauchy
## titan.princeton.edu/papers/claire/hertz-etal-99.ps
function upperbound(p::Poly)
    q, d = monic(p), degree(p)

    d == 0 && error("degree 0 is a constant")
    d == 1 && abs(q[0])

    a1 = abs(q[d-1])
    B = maximum([abs(q[i]) for i in 0:(d-2)])

    a,b,c = 1, -(1+a1), a1-B
    (-b + sqrt(b^2 - 4a*c))/2
end

"""

Use Descarte's rule of signs to compute number of *postive* roots of p

"""
¬(f::Function) = x -> !f(x)
function descartes_bound(p::Poly)
    bs = filter(¬_iszero, p.a)

    ## count terms
    ctr = 0
    bs = map(sign,bs)
    last = bs[1]
    for i in 2:length(bs)
        bs[i] == 0 && continue
        if last != bs[i]
            last = bs[i]
            ctr += 1
        end
    end
    ctr
end

## Descarte's upperbound on possible zeros in (0,1)
function DesBound(p::Poly)
    q = Tλ(R(p),1)
    descartes_bound(q)
end

## Interval with [c/2^k, (c+1)/2^k)
immutable Intervalkc
    k::Int
    c::Int
end

## convenience for two useful functions
DesBound(p,node::Intervalkc) = DesBound(Pkc(p,node))
function Pkc(p::Poly, node::Intervalkc)
    k,c = node.k, node.c
    Pkc(p, k, c)
end

## how we store the state of our algorithm to find zero in [0,1]
type State
    Internal                    # DesBound > 1
    Exact                       # c/2^k is a root
    Isol                        # DesBound == 1
    p::Poly
end
State(p) = State(Intervalkc[], Set{Rational}(), Intervalkc[], p)

## from naming convention of paper
function initTree(st::State)
    node = Intervalkc(0,0)
    ## check 0 and 1
    for x in [0//1, 1//1]
        if _iszero(polyval(st.p,x))
            push!(st.Exact, x)
            st.p,k = multiplicity(st.p, x)
        end
    end

    n = DesBound(Pkc(st.p, node))
    if n == 0
        ## nothing all done
    elseif n == 1
        push!(st.Isol, node)
    else
        push!(st.Internal, node)
    end
end

## get next node. Need to check that length(st.Internal) > 0, as it would throw error if no more Internal nodes
getNode(st::State) = shift!(st.Internal)

## add successors to [c,c+1]/2^k -> [2c, 2c+1]/2^(k+1) & [2(c+1), 2(c+1)+1]/2^(k+1) and check if c+1
function addSucc(st::State, node)
    k,c = node.k, node.c
    l = Intervalkc(k+1, 2c  )
    r = Intervalkc(k+1, 2c+1)
    p,m = multiplicity(st.p, (2c+1)//2^(k+1)); m > 0 && push!(st.Exact,  (2c+1)//2^(k+1))
    p,m = multiplicity(st.p, (2c+2)//2^(k+1)); m > 0 && push!(st.Exact,  (2c+2)//2^(k+1))
    for i in [l,r]
        n = DesBound(Pkc(p,i))
        n == 1 && push!(st.Isol, i)
        n > 1 && push!(st.Internal, i)
    end
end

## update state to resolve for roots in [0,1]
## st should have st.Interval = [], st.Exact and st.Isol
## This is basic search, not the more complicated one of the paper. To modify search
## strategy, the push!(st.Internal, i) would be changed in `addSucc`.
function find_in_01(p::Poly)
    st = State(p)

    initTree(st)
    for x in [0,1]
        p,k = multiplicity(p,x)
        k > 0 && push!(st.Exact, x)
    end
    while length(st.Internal) > 0
        node = getNode(st)
        node.k > 30 && error("too small?") # XXX arbitrary. Theory says a delta exists, this is δ > 1/2^30
        addSucc(st, node)
    end
    st
end

## get root from an Isol interval
function getIsolatedRoot(p::Poly, a::Real, b::Real)
    pa, pb = [polyval(p,x) for x in (a,b)]
    pab = sign(pa) * sign(pb)
    if pa == 0 p,k = multiplicity(p,a) end
    if pb == 0 p,k = multiplicity(p,b) end
    pa, pb = [polyval(p,x) for x in (a,b)]
    pab = sign(pa) * sign(pb)

    if pab < 0
        Roots.fzero(p, [a,b])
    else
        ## these aren't quite right, as we might get a double root this way
        abs(pa) <= 4eps() && return(a)
        abs(pb) <= 4eps() && return(b)
        error("Can't find root of $p in [$a, $b]")
    end
end

function getIsolatedRoot(p::Poly, i::Intervalkc)
    a,b = [i.c, i.c+1]//2^(i.k)
    c = (a+b)//2
    polyval(p,c) == 0 && return(c) # keeps rational if c/2^k type
    getIsolatedRoot(p, a, b)
end

## return positive roots of a square-free polynomial
function pos_roots(p::Poly)

    rts = Any[]

    degree(p) == 0 && return(rts)
    if degree(p) == 1
        c = -p[0]/p[1]
        if c > 0  return([c]) else return(rts) end
    end

    ## in case 0 is a root, we remove from p. Should be unnecessary, but ...
    p,k = multiplicity(p, 0)

    M = Roots.upperbound(p)

    for k in 0:floor(Integer,M)  ## we consider [k, k+1) for roots until [k, ∞) has just 0, 1 more roots.
        q = Tλ(p, k)      ## shift p left by k so this is p on [k, ∞)

        no_left = descartes_bound(q)
        if no_left == 0
            break
        elseif no_left == 1
            _iszero(polyval(p, k)) ? push!(rts, k) : push!(rts, getIsolatedRoot(p, k, M+1))
            break
        end

        ## more to do, [k, ∞) might have 2 or more roots...
        st = find_in_01(q)
        for i in st.Exact
            push!(rts, k + i)
            p, mult = multiplicity(p, k+i)
        end
        for i in st.Isol
            push!(rts, k + getIsolatedRoot(q,i))
        end
    end
    rts
end
neg_roots(p::Poly) = pos_roots(polyval(p, -variable(p)))

## use Rational{BigInt} to find real roots
## returns Any[] array as we might mix Rational{BigInt} and BigFloat values
function real_roots_sqfree(p::Poly)
    ## Handle simple cases
    degree(p) == 0 && error("Roots for degree 0 polynomials are either empty or all of R.")
    degree(p) == 1 && return([ -p[0]/p[1] ])

    rts = Any[]

    p,k = multiplicity(p, 0); k > 0 && push!(rts, 0) # 0
    append!(rts, pos_roots(p))                       # positive roots
    append!(rts, map(-,neg_roots(p)))                # negative roots
    rts
end

function real_roots{T <: QQ}(p::Poly{T})
    p = convert(Poly{Rational{BigInt}}, p)
    p = divrem(p, bgcd(p, Polynomials.polyder(p)))[1] # square free p/gcd(p, p')
    p = p*(1//p[degree(p)])

    ## factor first
    d = factor(p)
    d = filter((k,v) -> isa(k, Rational), d)
    rts = collect(keys(d))
    for rt in rts
        (p,k) = multiplicity(p, rt)
    end

    if degree(p) > 0
        new_rts = real_roots_sqfree(p)
        sort([rts; new_rts])
    else
        sort(rts)
    end
end

function real_roots(p::Poly)
    ## Handle simple cases
    p = convert(Poly{Float64}, p)
    if degree(p) > 1
        u, p, w, residual= agcd(p)
    end

    real_roots_sqfree(p)
end

### Some functions to find all rational roots of a poly in Z[x], Q[x]
## return set of divisors of |n|, e.g. 12 -> [1,2,3,4,6,12]
function divisors(n::Integer)
    n == 0 && return([1])
    d = factor(abs(n))
    out = [1]
    for (k,v) in d
        out = unique(out * [k^i for i in 0:v]')
    end
    sort(out)
end

# find rational roots of p ∈ Z[x]
# XXX most simple-minded approach around XXX this should be improved!
function rational_roots{T <: Integer}(p::Poly{T})
    rts = Rational{T}[]
    if p[0] == 0
        push!(rts, 0)
        p, k = multiplicity(p, 0)
    end

    ds = divisors(p[0])
    ns = divisors(p[degree(p)])
    poss = unique([d//n for d in ds, n in ns])
    for rt in poss
        for s in [-1,1]
            prt = s * rt
            polyval(p, prt) == 0 && push!(rts, prt)
        end
    end
    rts
end

# find rational roots of p ∈ Q[x] by using lcm(p)*p ∈ Z[x]
function rational_roots{T <: Integer}(p::Poly{Rational{T}})
    q = lcm(map(r -> r.den, coeffs(p))) * p
    q = convert(Poly{BigInt}, q)
    rational_roots(q)
end

## factor over the rationals, then call multroot on the remainder
## return d where d a dictionary of roots
## Okay, we should use keys which are `variable(p) - r` and not `r`...
function factor{T <: QQ}(p::Poly{T})
    rts = rational_roots(p)
    d = Dict{Number, Int}()
    for r in rts
        p,k = multiplicity(p,r)
        d[r] = k
    end
    ##
    if degree(p) > 0
        p = convert(Poly{Float64}, p)
        for (r,k) in zip(multroot(p)...)
            d[r] = k
        end
    end
    d
end

## Julian interface to multroot
function factor(p::Poly)
    degree(p) == 0 && DomainError() # degree 0
    d = Dict{Number, Int}()
    for (r,k) in zip(multroot(p)...)
        d[r] = k
    end
    d
end


text:
export HMatrix, size, type, tol, row, col, data

struct HMatrix
    m_type::UInt8
    m_size::Tuple{UInt32,UInt32}
    m_data::AbstractArray # full matrix or compressed format using ABt decomposition
    m_tol::T where {T<:Real}
    m_row::Tuple{Array{UInt32,1},Array{UInt32,1},Array{UInt32,1},Array{UInt32,1}}
    m_col::Tuple{Array{UInt32,1},Array{UInt32,1},Array{UInt32,1},Array{UInt32,1}}
    m_chd::Vector{HMatrix}
end

# accessors
Base.:size(hm::HMatrix) = hm.m_size
type(hm::HMatrix) = hm.m_type
tol(hm::HMatrix) = hm.m_tol
row(hm::HMatrix) = hm.m_row
col(hm::HMatrix) = hm.m_col
row(hm::HMatrix,i::U) where{U<:Integer} = hm.m_row[i]
col(hm::HMatrix,i::U) where{U<:Integer} = hm.m_col[i]
data(hm::HMatrix) = hm.m_data

# constructors

function HMatrix(m::T,n::T,tol::U,v) where {T<:Integer,U<:Real}
    #= Single value constructor which fills a leaf =#
    rowcol = tuple(UInt32[],UInt32[],UInt32[],UInt32[])
    return HMatrix(1,(m,n),ABtDecomposition(fill(v,(m,1)),fill(v,(1,n))),tol,rowcol,rowcol,HMatrix[])
end

max_stars_repo_name: marcbakry/HierarchicalMatrices

max_issues_repo_name: marcbakry/HierarchicalMatrices

max_forks_repo_name: marcbakry/HierarchicalMatrices


text:
module Prototype

#=
Playing around with prototype for plate method:
-----------------------------------------------
- Setup a default method signature, which dispaches solely on abstract types and throws an error (fall thru)
- result = analyticalSolution(  <: governing diff. eqn.     PlateBendingThin, PlateBendingThick, PlaneStress, PlaneStrain, etc.
                                <: analysis                 LinearElastic, TimeHistory, Modal, Buckling, SIF, etc.
                                <: approach                 Timoshenko, Irwin, etc.
                                <: dimensions               fully describe problem geometry
                                <: material                 isotropic, functionally graded, etc.
                                <: boundary condition       simply supported, clamped, spring, etc.
                                <: source                   concenetrated load, point load, line moment, etc.
                                <: result specification     result components to return
                             )
- at some point there should be some convenience factories for standard and well known numericla examples often presemted in the numerical examples section of papers.
- governing DE, analysis, dimensions, material, boundary condition and source are concrete types directly accepted from CivilEngineering.jl
- approach and result specification must be stated explicitly by the user (in the testing suite)
=#

end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
methods(factorial)

# %%
factorial(20)

# %%
factorial(21)

# %%
@which factorial(21)

# %%
Base._fact_table64

# %%
Base._fact_table64[1] = 666

# %%
factorial(1)

# %%
factorial(0.5)

# %%
using SpecialFunctions

# %%
factorial(0.5)

# %%
gamma(1.5)

# %%
@which factorial(0.5)

# %%
methods(factorial)

# %%
factorial(20)

# %%
@which factorial(20)

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
export subpopAgePlot

function subpopAgePlot( mpSim::MPsim, timeGrid::Vector{T}, ageRes::Real,
    ageType::Symbol, subpops::Subpopulation...; showPlot::Bool=true,
    savePlot::Bool=false, filename::AbstractString="",
    timeFactor::Real=12.0 ) where T <: Real

    savePlot |= filename != ""

    # Don't do anything if the plot mustn't be shown or saved.
    if !( showPlot || savePlot )
        return
    end  # if !( showPlot || savePlot )

    # Generate the report.
    ageReport = subpopulationAgeReport( mpSim, timeGrid, ageRes, ageType,
        subpops... )

    if isempty( ageReport ) || all( subpop -> isempty( ageReport[ subpop ] ),
        collect( keys( ageReport ) ) )
        return
    end  # if isempty( ageReport ) || ...

    ageTypeString = ageTypesS[ageType]

    # Get file name ready.
    if savePlot
        # If the user doesn't specify a filename, a random one is generated.
        if filename == ""
            dateStr = Dates.format( now(), "yyyymmdd HHMMSS" )
            fileroot = string( "Subpopulation ", ageTypeString, " plot (",
                dateStr, ") of " )
            extension = "svg"
        else
            extension = split( filename, "." )
            fileroot = length( extension ) == 1 ? extension[1] :
                join( extension[1:(end - 1)], "." )
            extension = length( extension ) == 1 ? "" : extension[end]
            extension = extension ∈ extensions ? extension : "svg"
        end  # if filename == ""

        if !ispath( dirname( fileroot ) )
            mkpath( dirname( fileroot ) )
        end  # if !ispath( dirname( fileroot ) )
    end  # if savePlot

    # Generate the plot.
    for subpop in filter( subpop -> haskey( ageReport, subpop.name ), subpops )
        report = ageReport[subpop.name]

        plotTitle = string( "Subpopulation ", subpop.name, ", ", ageTypeString,
            " plot" )

        timePoints = report[:, :timePoint] / timeFactor

        # Age summary plot.
        ageSummary = report[:, 2:6]

        if all( ageSummary[:, :mean] .=== missing )
            continue
        end  # if all( ageSummary[:, :mean] .=== missing )

        ymax = maximum( filter( age -> age isa Float64, ageSummary[:, :max] ) )
        ymax /= timeFactor
        agePlt = plot( timePoints,
            Matrix( ageSummary[:, [:mean, :median]] ) / timeFactor,
            title="Summary", labels=["mean" "median"], lw=2, ylim=[0, ymax] )
        plot!( agePlt, timePoints,
            Matrix( ageSummary[:, [:min, :max]] ) / timeFactor,
            labels=nothing, lw=1, color=:black )

        # Age distribution plot.
        ageDist = report[:, 7:end]
        ageLabels = parse.( Float64, string.( names( ageDist ) ) ) / timeFactor
        ageDist = Matrix( ageDist )
        zmax = maximum( ageDist )
        ageDistPlt = surface( timePoints, ageLabels, ageDist, zlims=[0, zmax],
            camera=(30, 60) )

        plt = plot( agePlt, ageDistPlt, size=(960, 1080), title=plotTitle,
            layout=grid( 2, 1 ) )

        # Show plot if needed.
        if showPlot
            gui( plt )
        end  # if showPlot

        # Save plot if needed.
        if savePlot
            plotfilename = string( fileroot, " ", subpop.name, ".", extension )
            savefig( plt, plotfilename )
        end  # if savePlot
    end  # for subpop in filter( ..., subpops )

end  # subpopAgePlot( mpSim, timeGrid, ageRes, ageType, subpops, showPlot,
     #   savePlot, fileName, timeFactor )

subpopAgePlot( mpSim::MPsim, timeRes::Real, ageRes::Real, ageType::Symbol,
    subpops::Subpopulation...; breakdownType::Symbol=:none, showPlot::Bool=true,
    savePlot::Bool=false, filename::AbstractString="", timeFactor::Real=12.0 ) =
    subpopAgePlot( mpSim, MP.generateTimeGrid( mpSim, timeRes ), ageRes,
    ageType, subpops..., showPlot=showPlot, savePlot=savePlot,
    filename=filename, timeFactor=timeFactor )


text:
"""
    data = load_data_set(name)

Load a small data set for numerical testing.
"""
function load_data_set(dataset="")
   name = lowercase(dataset)
   if name=="small5ren" || name=="small5" || name=="small5ren2"
      if name=="small5ren"
         pedlist = [  10  16  12  16   8  10  15  16  18  15  18  14  10  18  14
                       7  17   1  17  11   7   9  17  19  11  19   4   5  19   2];
         genid = [ 14,  2,  9,  5,  7,  13,  1,  3,  6]
      elseif name=="small5ren2"
         pedlist = [  10  16  12  16   8  10  15  16  16  15  17  14  10  17  14
                       7  17   1  17  11   7   9  17  16  11  16   4   5  17   2];
         genid = [ 14,  2,  9,  5,  7,  13,  1,  3,  6]
      else
         pedlist =[
            0  0  0  0  0  0  2  1  2  7  7  11  11   9  11
            0  0  0  0  0  0  5  4  3  6  4   8  10  13  10
         ]
         genid = [2,5,6,8,10,12,13,14,15]
      end
      G = [
         1.19677     0.542311     0.150994    0.457376    0.162172    0.0959141  -0.0896063   0.513584    -0.0896064
         0.542311    1.45157      0.141889    0.475572    0.535284    0.0322102  -0.126011    0.00187444  -0.126011
         0.150994    0.141889     0.941952   -0.0249475   0.373563    0.0175262   0.410923    0.511502     0.410924
         0.457376    0.475572    -0.0249475   1.09059     0.22076     0.652166    0.0564771  -0.0585542    0.0564771
         0.162172    0.535284     0.373563    0.22076     0.778136   -0.0314916   0.318484    0.331571     0.318484
         0.0959141   0.0322102    0.0175262   0.652166   -0.0314916   1.34492     0.492344   -0.118302     0.492345
        -0.0896063  -0.126011     0.410923    0.0564771   0.318484    0.492344    0.798901    0.391289     0.77654
         0.513584    0.00187444   0.511502   -0.0585542   0.331571   -0.118302    0.391289    1.03575      0.380108
        -0.0896064  -0.126011     0.410924    0.0564771   0.318484    0.492345    0.77654     0.380108     0.798901
      ]
      y =[
         -1.671506
         0.467041
        -1.949929
         0.537337
        -0.888975
        -1.672478
         0.893532
        -1.278057
         0.040003
        -1.1635
         1.364332
        -0.475337
         0.596651
        -1.775682
         0.68465
      ]
      X = ones(length(y),1)
      A = AnimalBreedingDataSetRaw(pedlist,genid,G,y,X)
   else
      error("specify one of the following dataset:\n small5, small5ren, or small5ren2")
   end
   return A
end


text:
"""
    everygram(seq::Vector{T}; min_len::Int=1, max_len::Int=-1)where { T <: AbstractString}

Return all possible ngrams generated from sequence of items, as an Array{String,1}

# Example

```julia-repl
julia> seq = ["To","be","or","not"]
julia> a = everygram(seq,min_len=1, max_len=-1)
 10-element Array{Any,1}:
  "or"
  "not"
  "To"
  "be"
  "or not"
  "be or"
  "be or not"
  "To be or"
  "To be or not"
```

"""
function everygram(seq::Vector{T}; min_len::Int=1, max_len::Int=-1) where { T <: AbstractString}
    ngram = []
    if max_len == -1
        max_len = length(seq)
    end
    for n in range(min_len, stop=max_len)
        temp = ngramizenew(seq, n)
        ngram = append!(ngram, temp)
    end
    return(ngram)
end

"""
    padding_ngram(word::Vector{T}, n=1; pad_left=false, pad_right=false, left_pad_symbol="<s>", right_pad_symbol ="</s>") where { T <: AbstractString}

padding _ngram is used to pad both left and right of sentence and out putting ngrmas of order n

   It also pad the original input Array of string

# Example
```julia-repl
julia> example = ["1","2","3","4","5"]

julia> padding_ngrams(example,2,pad_left=true,pad_right=true)
 6-element Array{Any,1}:
  "<s> 1"
  "1 2"
  "2 3"
  "3 4"
  "4 5"
  "5 </s>"
```
"""
function padding_ngram(word::Vector{T}, n=1; pad_left=false, pad_right=false, left_pad_symbol="<s>", right_pad_symbol ="</s>") where { T <: AbstractString}
    local seq
    seq = word
    if pad_left == true
        prepend!(seq, [left_pad_symbol])
    end
    if pad_right == true
        push!(seq, right_pad_symbol)
    end
    return  ngramizenew(seq, n)
end

"""
    ngramizenew( words::Vector{T}, nlist::Integer...) where { T <: AbstractString}

ngramizenew is used to out putting ngrmas in set

# Example
```julia-repl
julia> seq=["To","be","or","not","To","not","To","not"]
julia> ngramizenew(seq ,2)
 7-element Array{Any,1}:
  "To be"
  "be or"
  "or not"
  "not To"
  "To not"
  "not To"
  "To not"
```
"""
function ngramizenew(words::Vector{T}, nlist::Integer...) where { T <: AbstractString}
    n_words = length(words)

    tokens = []

    for n in nlist
        for index in 1:(n_words - n + 1)
            token = join(words[index:(index + n - 1)], " ")
            push!(tokens,token)
        end
    end
    return tokens
end


text:
"""
    MenuAdventures.Parts

A sub-module that allows things to be part of other things.

```jldoctest
julia> using MenuAdventures

julia> using MenuAdventures.Testing

julia> using MenuAdventures.Parts

julia> import MenuAdventures: ever_possible

julia> @universe struct Universe <: AbstractUniverse
        end;

julia> @noun struct Room <: AbstractRoom
            already_lit::Bool = true
        end;

julia> @noun struct Person <: Noun
        end;

julia> @noun struct StickyThing <: Noun
        end;

julia> ever_possible(::Take, ::Reachable, ::StickyThing) = true;

julia> ever_possible(::Attach, ::Inventory, ::StickyThing) = true;

julia> cd(joinpath(pkgdir(MenuAdventures), "test", "Parts")) do
            check_choices() do interface
                you = Person(
                    "Brandon",
                    grammatical_person = second_person,
                    indefinite_article = "",
                )
                room = Room("room")
                universe = Universe(you, interface = interface)
                universe[room, you] = Containing()
                universe[room, StickyThing("sticky thing")] = Containing()
                universe
            end
        end
true
```
"""
module Parts

using MenuAdventures: Action, get_object, Inventory, Reachable, Relationship, Sentence, Verb
import MenuAdventures: argument_domains, ever_possible, possible_now, print_sentence, string_relationship_to, verb_for

"""
    Incorporating()

A is `Incorporating` B means B is part of A.

A [`Relationship`](@ref)
"""
struct Incorporating <: Relationship end
export Incorporating

verb_for(::Incorporating) = Verb("incorporate")

function string_relationship_to(thing_answer, ::Incorporating, parent_thing)
    buffer = IOBuffer()
    print(buffer, thing_answer.text)
    print(buffer, ' ')
    print(buffer, "attached to ")
    show(IOContext(buffer, :is_subject => false, :subject => get_object(thing_answer)), parent_thing)
    String(take!(buffer))
end

"""
    Attach()

`Attach` something from your [`Inventory`](@ref) to something else [`Reachable`](@ref).

An [`Action`](@ref).
By default, `ever_possible(::Attach, ::Reachable, _) = true`, that is, you can attach something sticky to anything.
"""
struct Attach <: Action end

export Attach

ever_possible(::Attach, ::Reachable, _) = true

function possible_now(_, sentence::Sentence{Attach}, domain::Reachable, thing)
    ever_possible(sentence.action, domain, thing) &&
    # can't attach something to itself
    sentence.argument_answers[1].object !== thing
end

function argument_domains(::Attach)
    Inventory(), Reachable()
end

function (::Attach)(universe, thing, parent_thing)
    universe[parent_thing, thing] = Incorporating()
    return false
end

function print_sentence(io, ::Attach, thing_answer, parent_thing_answer)
    print(io, "Attach ")
    print(io, thing_answer.text)
    print(io, " to ")
    print(io, parent_thing_answer.text)
end

end


text:
@testset "MatPolynomial tests" begin
    @testset "MatPolynomial" begin
        @polyvar x y
        P = MatPolynomial{Int}((i,j) -> i + j, [x^2, x*y, y^2])
        zP = zero(typeof(P))
        @test isempty(zP.Q)
        @test zP == 0
        p = polynomial(P)
        @test coefficients(p) == [2, 6, 12, 10, 6]
        @test monomials(p) == [x^4, x^3*y, x^2*y^2, x*y^3, y^4]
        for i in 1:3
            for j in 1:3
                @test P[i, j] == i + j
            end
        end
        for P in (MatPolynomial((i,j) -> i * j, [y, x]),
                  MatPolynomial((i,j) -> (3-i) * (3-j), monovec([y, x])),
                  MatPolynomial([1 2; 2 4], [y, x]),
                  MatPolynomial([4 2; 2 1], monovec([y, x])))
            @test P.Q.Q == [4, 2, 1]
            @test P.x[1] == x
            @test P.x[2] == y
        end
        P = MatPolynomial((i,j) -> ((i,j) == (1,1) ? 2 : 0), [x*y, x^2, y^2])
        Q = MatPolynomial([0 1; 1 0], [x^2, y^2])
        @test P == Q
        p = MatPolynomial([2 3; 3 2], [x, y])
        @test polynomial(p) isa AbstractPolynomial{Int}
        @test polynomial(p, Float64) isa AbstractPolynomial{Float64}

        @test -2*x + dot(-x - x^2, 0) + MatPolynomial{Int}((i,j)->1, [1,x]) == -(-x^2 - 1)
        P = MatPolynomial{Int}((i,j) -> i + j, [x^2, x*y, y^2])
        p = polynomial(P)
        @test !iszero(P)
        @test iszero(P-P)
        @test iszero(P-p)
        @test iszero(p-P)
        @test P + P == P + p
        @test x * P == x * p
        @test 2 * P == 2 * p
        @test P * (2x) == (2x) * p

        @test differentiate(MatPolynomial{Int}((i,j)->1, [x]), y) == 0
        @test differentiate(MatPolynomial{Int}((i,j)->1, [x, y]), y, 1) == 2x + 2y
        #@inferred differentiate(MatPolynomial{Int}((i,j)->1, [x, y]), y, 0) # FIXME failing at the moment

        @test MatPolynomial([2 3; 3 2], [x, y]) == 2x^2 + 2y^2 + 6x*y

        mp = MatPolynomial{Int}((i,j) -> i+j, [x])
        @test norm(mp) == norm(mp, 2) == 2.0

        @polyvar v[1:3]
        P = [1 2 3; 2 4 5; 3 5 6]
        p = MatPolynomial(P, v)
        #@inferred p(x => ones(3))
        @test p(v => ones(3)) == 31
        #@inferred subs(p, x => ones(3))
        @test subs(p, v => ones(3)) == 31
    end
    @testset "SOSDecomposition" begin
        @polyvar x y
#       @test isempty(SOSDecomposition(typeof(x)[]))
#       ps = [1, x + y, x^2, x*y, 1 + x + x^2]
#       P = MatPolynomial(SOSDecomposition(ps))
#       P.Q == [2 0 1 0 1; 0 1 0 0 0; 1 0 2 1 1; 0 0 1 1 0; 1 0 1 0 2]
#       P.x == [x^2, x*y, x, y, 1]
#       @test P == P
#       @test isapprox(MatPolynomial(SOSDecomposition(P)), P)
        @test sprint(show, SOSDecomposition([x+y, x-y])) == "(x + y)^2 + (x - y)^2"
        @testset "SOSDecomposition equality" begin
            @polyvar x y
            @test !isapprox(SOSDecomposition([x+y, x-y]), SOSDecomposition([x+y]))
            @test !isapprox(SOSDecomposition([x+y, x-y]), SOSDecomposition([x+y, x+y]))
            @test isapprox(SOSDecomposition([x+y, x-y]), SOSDecomposition([x+y, x-y]))
            @test isapprox(SOSDecomposition([x+y, x-y]), SOSDecomposition([x-y, x+y+1e-8]), ztol=1e-7)
        end
    end
end


text:
module HScss

export hscss, printScss, compile, hscss_get, hscss_getall
#tokenize, structure, modifiers, clumps

include("Units.jl")
include("Tokenize.jl")
include("write.jl")

function findUnit(token)
  for i in 1:length(units)
    if units[i] == token.val
      return true
    end
  end
  return false
end
# ==============================================================================
# clump related tokens for easier structuring...
# ==============================================================================
function clumps(token)
  fin = length(token)
  stack = []
  t = 1
  while t <= fin
    offset = 1
    if  isa(token[t], Alpha)
        if findUnit(token[t]) # test for units.
            token[t] = Unit(token[t].val)
            offset = 1
        elseif isa(token[t+1], Dotdot)
            token[t] = Attr(token[t].val)
            offset = 2
        elseif isa(token[t+1], Open)
            token[t] = Tag(token[t].val)
            offset = 2
        elseif  isa(token[t+1], End)
            token[t] = Value(token[t].val)
            offset = 2
        end

    elseif isa(token[t], Dot)
        if isa(token[t+1], Digit)
            token[t] = Float(parse(Float64, "." * token[t+1].val))
            offset = 2
        elseif isa(token[t+1], Alpha) && isa(token[t+2], Open)
            token[t] = Class(token[t+1].val)
            offset = 3
        end

    elseif isa(token[t], Dollar)
        if isa(token[t+1], Alpha) && isa(token[t+2], Dotdot)
            token[t] = Ident(token[t+1].val)
            offset = 3
        elseif isa(token[t+1], Alpha)
            token[t] = Reference(token[t+1].val)
            offset = 2
        end

    elseif isa(token[t], Digit)
        if t<fin-1 && isa(token[t+1], Dot) && isa(token[t+2], Digit)
            token[t] = Float(parse(Float64, token[t].val * "." * token[t+2].val))
            offset = 3
        elseif isa(token[t+1], Dot)
            token[t] = Float(parse(Float64, token[t].val * ".")) # Float(token[t].val * ".")
            offset = 2
        else # this must be an integer!
            token[t] = Int_(parse(Int, token[t].val))
            offset = 1
        end

    elseif isa(token[t], Hash) && (isa(token[t+1], Digit) || isa(token[t+1], Alpha))
            token[t] = Hex(token[t+1].val)
            offset = 2
    end

    if !isa(token[t], End) && !isa(token[t], Space)
      push!(stack, token[t])
    end
    t+=offset
    end
    stack
end
# ==============================================================================
# Push attributes to
# ==============================================================================
function pushAttributes(token)
  fin = length(token)
  stack = []
  t = 1
while t < fin

    offset = 1
    if isa(token[t], Attr) || isa(token[t], Ident) # AbstractLeft
        values = []
        while t+offset <= fin && isa(token[t+offset], AbstractRight)
            # Unpackage the standard typed vars
            if isa(token[t+offset], Int_) || isa(token[t+offset], Float) || isa(token[t+offset], Hex) || isa(token[t+offset], Value) || isa(token[t+offset], Alpha)
                push!(values, token[t+offset].val)
            else
                push!(values, token[t+offset])
            end

            offset +=1
        end

        if length(values) == 1
            values = values[1]
        end
        if isa(token[t], Ident)
            token[t] = token[t] => values
        else
            token[t] = token[t].val => values
        end
    end
    push!(stack, token[t])
    t+=offset

end
    t <= fin &&  push!(stack, token[t])
    return stack
end
# ==============================================================================
function modifiers(token)
  fin = length(token)
  stack = []
  t = 1
while t <= fin
    if (isa(token[t], AbstractRight)) && isa(token[t+1], Unit) # This seems fragile!
      unit = token[t+1].val
      value = token[t].val
      if     unit == "px"
        value = Px(value)
      elseif unit == "em"
        value = Em(value)
      elseif unit == "cm"
        value = Cm(value)
      elseif unit == "in"
        value = In(value)
      elseif unit == "deg"
        value = Deg(value)
      elseif unit == "mm"
        value = Mm(value)
      elseif unit == "pt"
        value = Pt(value)
      elseif unit == "pc"
        value = Pc(value)
      elseif unit == "ex"
        value = Ex(value)
      elseif unit == "ch"
        value = CH(value)
      elseif unit == "rem"
        value = Rem(value)
      elseif unit == "vw"
        value = Vw(value)
      elseif unit == "vh"
        value = Vh(value)
      elseif unit == "vmin"
        value = Vmin(value)
      elseif unit == "vmax"
        value = Vmax(value)
      elseif unit == '%'
        value = Pcnt(value)
      end
      push!(stack, value)
          t+=2
    elseif (isa(token[t], Attr)) && ( isa(token[t+1], Value) || isa(token[t+1], String) )
        push!(stack, token[t].val => String(token[t+1].val))
        t+=2
    else
        push!(stack, token[t])
        t+=1
    end
end

    return pushAttributes(stack)
end
# ==============================================================================
function structure(token)
  fin = length(token)
  t = 1
  vars, classes, body = [], [], []

function getnode() # get everything and
        stack = []
        tags = []
        while t <= fin
            node = token[t]
            t+=1
            if isa(node, Class) #|| isa(node, Tag) # Go deeper
                push!(stack, node.val => Dict( getnode()...))
            elseif isa(node, Tag) # Go deeper
                push!(tags, node.val => Dict( getnode()...))
                #push!(stack, Dict(node.val => getnode()))
            elseif isa(node, Close) # Go back
                push!(stack, "nodes" => tags)
                return Dict(stack...)
            else # Add to current
                push!(stack, node)
            end
        end
        push!(stack, "nodes" => tags)
        return stack
end

# In theory we chould only find the three principal node types at the base level
while t <= fin
        node = token[t]
        t+=1
        if isa(node, Pair) && isa(node[1], Ident)
            push!(vars, node[1].val => node[2])
        elseif isa(node, Ident)
            push!(vars, (node.val => node))

        elseif isa(node, Class)
            push!(classes, (node.val => getnode()))
        elseif isa(node, Tag)
            push!(body, (node.val => getnode()))
        elseif isa(node, Close)
        else
            println("---ERROR: ", node," should not be here! ---")
        end
end

    return (vars, classes, body)
end

function hscss(s)
    tokens = tokenize(s)
    return structure(modifiers(clumps(tokens.tokens)))
end

# ==============================================================================
# Find ONE named element
# ==============================================================================
function hscss_get(array, name)
    for e in array
        if e[1] == name
            return e
        end
    end
    return nothing
end
# ==============================================================================
# return all elements of name
# ==============================================================================
function hscss_getall(array, name)
    elements = []
    for e in array
        if e[1] == name
            push!(elements, e)
        end
    end
    return elements
end

# maybe a future addition...
function compile(s)
    vars, classes, body = hscss(s)
    values = []
    for b in vars
        push!(values, (b[1],"var"))
    end

    for b in classes
        push!(values, (b[1],"class"))
    end

    for b in classes
        push!(values, (b[1],"class"))
    end

    # function p_pair(a, b, i)
    #     print(" "^(i*2), "$a:")
    #     if isa(b, Array)
    #         for i in b
    #             print(" ", printvalue(i))
    #         end
    #     else
    #         print(printvalue(b))
    #     end
    #     print(";\n")
    # end

    function p_dict(d)
        for (k,v) in d
            if k == "nodes"
                p_dict(v)
            elseif isa(v, Dict) #
                #println(" "^(i*2), "$k{")
                p_dict(v)
                #println(" "^(i*2), "}")
            else # All attributes...
                if k == "style"
                    println(k)
                elseif isa(v, Ident)
                    println(v)
                end
            end
        end
    end

    for b in body
        p_dict(b[2])
    end

    return values
end

end # module


text:
using CSV,DataFrames,Dates,GLM,Statistics,LinearAlgebra,Measures,Plots; pyplot()

df = CSV.read("../data/oneOnEpsilonBlogs.csv",copycols = true)
len = size(df)[1]

dow = dayofweek.(Date.(df.Day,Dates.DateFormat("m/d/y")))
dayGroups = [[7], [1,6], [2,3,4,5]]
inds = [[in(d,grp) for d in dow] for grp in dayGroups]
df2 = DataFrame(Time=1:len,Users=df.Users,
                Group1=inds[1],Group2=inds[2])

trainRange1, futureRange1 = 100:130, 130:180
trainRange2, futureRange2 = 200:300, 300:320
trainRange3, futureRange3 = 400:500, 500:600
trainRange4, futureRange4 = 560:600, 600:630

function forecast(trainRange,futureRange)
    model = lm(@formula(Users ~ Time + Group1 + Group2),df2[trainRange,:])
    pred(time) =  dot(coef(model),[1,time,inds[1][time],inds[2][time]])
    model, pred.(trainRange), pred.(futureRange)
end

function forecastPlot(train,future)
    p = plot(train[1], df.Users[train[1]], label = "Observed Users",
            xlabel = "Day", ylabel = "Daily Users", c=:blue)
        plot!(train[1], train[2], label = "Train", c=:red)
        plot!(future[1], future[2], label = "Forecast",c=:magenta)
        plot!(future[1], df.Users[future[1]], label = "Actual Users",c=:green)
    return p
end

model1, train1, fcst1 = forecast(trainRange1,futureRange1)
println(model1)

default(legend = :topleft)
p1 = forecastPlot((trainRange1, train1), (futureRange1, fcst1))

_, train2, fcst2 = forecast(trainRange2, futureRange2)
default(legend = false)
p2 = forecastPlot((trainRange2, train2), (futureRange2, fcst2))

_, train3, fcst3 = forecast(trainRange3, futureRange3)
p3 = forecastPlot((trainRange3, train3), (futureRange3, fcst3))

_, train4, fcst4 = forecast(trainRange4, futureRange4)
p4 = forecastPlot((trainRange4, train4), (futureRange4, fcst4))
plot(p1, p2, p3, p4, size = (900,600), margin = 5mm)


text:
#
# test data
#
"create a small instance (and topology) with a single pipe"
function single_pipe(;length=100.0, flow=200.0)
    nodes = [Node(0,0), Node(length,0)]
    demand = [-flow, flow]
    press = fill(Bounds(60.0, 80.0), size(nodes))
    diams = [Diameter(t...) for t in [(0.8, 1.0),(1.0, 1.2),(1.2, 1.5)]]
    inst = Instance(nodes, demand, press, diams)
    topo = Topology(nodes, [Arc(1,2)])
    inst, topo
end

"Create a grid of square cells with m by n nodes and given edge width."
function squaregrid(m::Int, n::Int, width::T; antiparallel=false) where T<:Real
    nodes = Node[]
    sizehint!(nodes, m*n)
    for j in 1:n
        x = width * (j - 1)
        for i in 1:m
            push!(nodes, Node(x, width * (m - i)))
        end
    end
    nodeidx = reshape(1:m*n, (m,n))

    arcs = Arc[]
    narcs = (m-1)*n + m*(n-1)
    sizehint!(arcs, antiparallel ? 2*narcs : narcs)
    # columns top-down
    for j in 1:n
        for i in 1:m-1
            push!(arcs, Arc(nodeidx[i,j], nodeidx[i+1,j]))
            antiparallel && push!(arcs, Arc(nodeidx[i+1,j], nodeidx[i,j]))
        end
    end
    # rows left-right
    for i in 1:m
        for j in 1:n-1
            push!(arcs, Arc(nodeidx[i,j], nodeidx[i,j+1]))
            antiparallel && push!(arcs, Arc(nodeidx[i,j+1], nodeidx[i,j]))
        end
    end

    Topology(nodes, arcs)
end


text:
using GLVisualize, GLWindow, GLAbstraction

function volume_data(N)
	volume 	= Float32[sin(x/15.0)+sin(y/15.0)+sin(z/15.0) for x=1:N, y=1:N, z=1:N]
	max 	= maximum(volume)
	min 	= minimum(volume)
	volume 	= (volume .- min) ./ (max .- min)
end

if !isdefined(:runtests)
	window = glscreen()
	timesignal = bounce(linspace(0f0,1f0,360))
    volumedata = volume_data(128)
else
    using NIfTI
    volumedata = niread(assetpath("brain.nii")).raw
end
const record_interactive = true

volume = visualize(volumedata, :mip)

_view(volume, window)

if !isdefined(:runtests)
	renderloop(window)
end


text:
# Important that oldest is removed first !
# Example showing why this is important:
# For instance, if the lattice is:
#            :Unbounded - ...
#           /
# :Unbounded
#           \
#            :Optimal - ...
# The lower node will give optimality cuts and fill the cutpruner of the root node
# Meanwhile, the upper node will not give any cut so the root not will stay unbounded
# Since the root node is :Unbounded it will not have dual variables for the stats, so
# it may happen (for instance with AvgCutPruner) that all the cuts are at tie !
# When the upper node stops being :Unbounded and gives its cut that is necessary for
# the root node to stop being :Unbounded. However, if the lower node replaces it with
# another optimality cut it will be ignored !
for algo in [AvgCutPruningAlgo(3), DecayCutPruningAlgo(3)]
    @testset "Oldest removed first for $(typeof(algo))" begin
        pruner = CutPruner{2, Int}(algo, :Max)
        addcuts!(pruner, [1 0], [0], [true])
        # add redundants cuts: this cut will not be added to pruner
        addcuts!(pruner, [1 0], [0], [true])
        @test ncuts(pruner) == 2

        # add redundants cuts: the last cut will not be added to pruner
        addcuts!(pruner, [1 1; 1 1], [0, 0], [true, true])
        @test ncuts(pruner) == 3

        addcuts!(pruner, [2 0], [0], [true])
        addcuts!(pruner, [3 0], [0], [true])
        addcuts!(pruner, [4 0], [0], [true])
        @test sort(pruner.A[:,1]) == [2, 3, 4]
        addcuts!(pruner, [5 0], [0], [true])
        @test sort(pruner.A[:,1]) == [3, 4, 5]
        # However, order should only be used to break ties
        addcuts!(pruner, [6 0], [0], [false])
        @test sort(pruner.A[:,1]) == [3, 4, 5]
        addcuts!(pruner, [7 0; 8 0], [0, 0], [true, true])
        @test sort(pruner.A[:,1]) == [5, 7, 8]
    end
end


text:
import Base.reshape, Base.vec
export reshape, vec
export sign, curvature, monotonicity, evaluate, conic_form!

struct ReshapeAtom <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr}
    size::Tuple{Int, Int}

    function ReshapeAtom(x::AbstractExpr, m::Int, n::Int)
        if m * n != length(x)
            error("Cannot reshape expression of size $(x.size) to ($(m), $(n))")
        end
        return new(:reshape, objectid(x), (x,), (m, n))
    end
end

function sign(x::ReshapeAtom)
    return sign(x.children[1])
end

function monotonicity(x::ReshapeAtom)
    return (Nondecreasing(),)
end

function curvature(x::ReshapeAtom)
    return ConstVexity()
end

function evaluate(x::ReshapeAtom)
    return reshape(evaluate(x.children[1]), x.size[1], x.size[2])
end

function conic_form!(x::ReshapeAtom, unique_conic_forms::UniqueConicForms)
    return conic_form!(x.children[1], unique_conic_forms)
end

reshape(x::AbstractExpr, m::Int, n::Int) = ReshapeAtom(x, m, n)
vec(x::AbstractExpr) = reshape(x, length(x), 1)


text:
# bash: export JULIA_NUM_THREADS=4
include("../RLAlgorithms.jl")
using .RLAlgorithms
using Plots
using Base.Threads
using Random: MersenneTwister
using Flux
using Flux: Zygote.Grads, Zygote.Params
import Flux.Optimise: update!
using Statistics: mean
using StatsBase: sample, Weights
import.Base: +

function +(x::Grads, y::Grads)
    return Grads(IdDict(p => (get(x.grads, p, nothing) + y.grads[p]) for p in y.params), y.params)
end
+(::Nothing, y) = y

function update!(opt, ps1::Params, ps2::Params, gs::Grads)
    for (p1, p2) in zip(ps1.order, ps2.order)
        gs.grads[p2] !== nothing && update!(opt, p1, gs.grads[p2])
        p2 .= p1
    end
end

struct Policy{S, A, V}
    shared::S
    action::A
    value::V
end
(m::Policy)(s) = (l = m.shared(s); (softmax(m.action(l)), m.value(l)))
Flux.@functor Policy

function episode!(policy, episode, env, maxt, rng)
    reset!(env)
    done = false
    for _ in 1:maxt-1
        s = copy(env.current)
        logits, v = policy(s)
        a = sample(rng, Weights(logits))
        _, r, done = step!(env, a)
        append!(episode, s, a, r, v[1])
        done && break
    end
    if !done
        s = env.current
        p, v = policy(s)
        append!(episode, s, argmax(p), v[1], v[1])
    end
    nothing
end

function grad(policy, ps, s, a, r, adv)
    return gradient(ps) do
        p, v = policy(s)
        logp = log.(p)
        entropy = -mean(sum(p .* logp, dims=1))
        return -mean(gather(logp, a) .* adv) + 0.5f0 * Flux.mse(vec(v), r) - 1f-2 * entropy
    end
end

function run!(sharedpolicy, opt, sharedenv; epochs=500, maxt=100, γ=99f-2)
    rewards = zeros(Float32, div(epochs, 5))
    sharedrng = [MersenneTwister(rand(1:100)) for _ in 1:nthreads()]
    sharedlockobj = ReentrantLock()
    sharedps = params(sharedpolicy)

    @threads for _ in 1:nthreads()
        rng = sharedrng[threadid()]
        env = deepcopy(sharedenv)
        policy = deepcopy(sharedpolicy)
        ps = params(policy)
        gs = Grads(IdDict(), Params())
        episode = PGEpisode{Float32, Int64, Float32}(length(env.observationspace), 1, maxt, γ=γ)

        for i in 1:epochs
            episode!(policy, episode, env, maxt, rng)
            GAEλ!(episode)
            discount!(episode)

            s, a, r, adv = batch(episode)
            gs += grad(policy, ps, s, a, r, adv)

            if i % 5 == 0
                lock(sharedlockobj)
                update!(opt, sharedps, ps, gs)
                unlock(sharedlockobj)
                gs = Grads(IdDict(), Params())
            end
            if i % 5 == 0
                ind = div(i, 5)
                if threadid() == (ind % nthreads()) + 1
                    rewards[ind] = test(s -> Flux.onecold(policy(s)[1]), env)
                end
            end
        end
    end
    return rewards
end

env = CartPole{Float32}()

sharedpolicy = Policy(
    Dense(length(env.observationspace), 64, relu),
    Dense(64, length(env.actionspace)),
    Dense(64, 1)
)

r = run!(sharedpolicy, ADAM(0.0015, (0.5, 0.99)), env)
plt = plot(r, labels="Reward")
display(plt)
#savefig(plt, "a3c")

#=
older Zygote version:
gs = Grads(IdDict())
function +(x::Grads, y::Grads)
    return Grads(IdDict(k => (get(x.grads, k, nothing) + y.grads[k]) for k in keys(y.grads)))
end
=#

#=
+(x::Tuple, y::Tuple) = Tuple([(x[k] + y[k]) for k in keys(x)])
function +(x::NamedTuple, y::NamedTuple)
    NamedTuple{keys(x)}([(x[k] === nothing ? nothing : x[k] + y[k]) for k in keys(x)])
end

function update!(opt, m1::T, m2::T, gs) where T <: AbstractArray
    update!(opt, m1, gs)
    m2 .= m1
    gs .= zero(eltype(gs))
end

function update!(opt, m1, m2, gs)
    for k in keys(gs)
        gs[k] === nothing && continue
        if typeof(k) <: Symbol
            update!(opt, getfield(m1, k), getfield(m2, k), gs[k])
        else
            update!(opt, m1[k], m2[k], gs[k])
        end
    end
end

sharedinitgs = (gradient(sharedpolicy) do p
    m = p(zeros(Float32, length(env.observationspace)))
    return Flux.mse(m[1], 0f0) + Flux.mse(m[2], 0f0)
end)[1]

gs = deepcopy(sharedinitgs)

gs += (gradient(policy) do model
    p, v = model(s)
    logp = log.(p)
    entropy = -mean(sum(p .* logp, dims=1))
    return -mean(gather(logp, a) .* adv) + 5f-1 * Flux.mse(vec(v), r) - 1f-2 * entropy
end)[1]

update!(opt, sharedpolicy, policy, gs)
=#


text:
type DB
    file::UTF8String
    handle::Ptr{Void}
    changes::Int

    function DB(f::UTF8String)
        handle = [C_NULL]
        f = isempty(f) ? f : expanduser(f)
        if @OK sqliteopen(f, handle)
            db = new(f, handle[1], 0)
            # register(db, regexp, nargs=2)
            finalizer(db, x -> sqlite3_close(handle[1]))
            return db
        else # error
            sqlite3_close(handle[1])
            sqliteerror()
        end
    end
end

function DB(f::AbstractString; actualize::Bool=false)
    res = DB(utf8(f))
    actualize && actualize!(res)
    return res
end
DB() = DB(":memory:")

Base.close(db::DB) = sqlite3_close_v2(db.handle)

type World
    actualized::Bool
    db::DB

    World() = new(false)
end

const actuality = World()

function actualize!(db::DB)
    actuality.db = db
    actuality.actualized = true;
end

type Stmt
    db::DB
    handle::Ptr{Void}
    sql::String
    table::String

    function Stmt(db::DB, sql::AbstractString, table::String="")
        handle = [C_NULL]
        sqliteprepare(db, sql, handle, [C_NULL])
        stmt = new(db, handle[1], sql, table)
        finalizer(stmt, x -> sqlite3_finalize(handle[1]))
        return stmt
    end
end

Base.close(stmt::Stmt) = sqlite3_finalize(stmt.handle)

function sqliteprepare(db, sql, stmt, null)
    @CHECK db sqlite3_prepare_v2(db.handle, utf8(sql), stmt, null)
end

function execute!(stmt::Stmt)
    r = sqlite3_step(stmt.handle)
    if r == SQLITE_DONE
        sqlite3_reset(stmt.handle)
    elseif r != SQLITE_ROW
        sqliteerror(stmt.db)
    end
    return r
end


text:
using HA_solver
using Test, Random

include("test_Ayiagari.jl")
include("test_utils.jl")
include("NGS_test.jl")


text:
using Geoo, DataStructures, Finder

curitiba = Coordinate("-25.44777778, -49.23055556")
caucuta = Coordinate("22.55694444, 88.33666667")
brasilia = Coordinate("-15.871329,-47.912457")
novayork = Coordinate("43.121948,-77.668278")

p1 = People("p1", curitiba, "trab1", true)
p2 = People("p2", caucuta, "trab2", true)
p3 = People("p3", brasilia, "trab3", true)
p4 = People("p4", novayork, "trab4", true)
D = SortedDict{Int, People}()
P = [p1 p2 p3 p4]
cc = 0

for p in P
    push!(D, cc => p)
    global cc += 1
end

function dictweightcoordinate(D)
    return  function weightcoordinate(s,t)
        return coordist(D[s].coord, D[t].coord)
    end
end


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

@doc raw"""WebhookConversion describes how to call a conversion webhook

    IoK8sApiextensionsApiserverPkgApisApiextensionsV1WebhookConversion(;
        clientConfig=nothing,
        conversionReviewVersions=nothing,
    )

    - clientConfig::IoK8sApiextensionsApiserverPkgApisApiextensionsV1WebhookClientConfig : clientConfig is the instructions for how to call the webhook if strategy is `Webhook`.
    - conversionReviewVersions::Vector{String} : conversionReviewVersions is an ordered list of preferred `ConversionReview` versions the Webhook expects. The API server will use the first version in the list which it supports. If none of the versions specified in this list are supported by API server, conversion will fail for the custom resource. If a persisted Webhook configuration specifies allowed versions and does not include any versions known to the API Server, calls to the webhook will fail.
"""
mutable struct IoK8sApiextensionsApiserverPkgApisApiextensionsV1WebhookConversion <: SwaggerModel
    clientConfig::Any # spec type: Union{ Nothing, IoK8sApiextensionsApiserverPkgApisApiextensionsV1WebhookClientConfig } # spec name: clientConfig
    conversionReviewVersions::Any # spec type: Union{ Nothing, Vector{String} } # spec name: conversionReviewVersions

    function IoK8sApiextensionsApiserverPkgApisApiextensionsV1WebhookConversion(;clientConfig=nothing, conversionReviewVersions=nothing)
        o = new()
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1WebhookConversion, Symbol("clientConfig"), clientConfig)
        setfield!(o, Symbol("clientConfig"), clientConfig)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1WebhookConversion, Symbol("conversionReviewVersions"), conversionReviewVersions)
        setfield!(o, Symbol("conversionReviewVersions"), conversionReviewVersions)
        o
    end
end # type IoK8sApiextensionsApiserverPkgApisApiextensionsV1WebhookConversion

const _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1WebhookConversion = Dict{Symbol,Symbol}(Symbol("clientConfig")=>Symbol("clientConfig"), Symbol("conversionReviewVersions")=>Symbol("conversionReviewVersions"))
const _property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1WebhookConversion = Dict{Symbol,String}(Symbol("clientConfig")=>"IoK8sApiextensionsApiserverPkgApisApiextensionsV1WebhookClientConfig", Symbol("conversionReviewVersions")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1WebhookConversion }) = collect(keys(_property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1WebhookConversion))
Swagger.property_type(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1WebhookConversion }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1WebhookConversion[name]))}
Swagger.field_name(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1WebhookConversion }, property_name::Symbol) =  _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1WebhookConversion[property_name]

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

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


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Types and structures of SGP4 model.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export SGP4_GravCte, SGP4_Structure

"""
Gravitational constants for SGP4.

# Fields

* `R0`: Earth equatorial radius [km].
* `XKE`: √GM [er/s]^(3/2).
* `J2`: The second gravitational zonal harmonic of the Earth.
* `J3`: The thrid gravitational zonal harmonic of the Earth.
* `J4`: The fourth gravitational zonal harmonic of the Earth.

"""
@with_kw struct SGP4_GravCte{T<:Real}
    R0::T
    XKE::T
    J2::T
    J3::T
    J4::T
end

@with_kw mutable struct SGP4_DeepSpace{T<:Real}
    atime::T  = 0.0
    xli::T    = 0.0
    xni::T    = 0.0
    xnq::T    = 0.0
    xfact::T  = 0.0
    ssl::T    = 0.0
    ssg::T    = 0.0
    ssh::T    = 0.0
    sse::T    = 0.0
    ssi::T    = 0.0
    xlamo::T  = 0.0
    omegaq::T = 0.0
    omgdt::T  = 0.0
    gmst::T   = 0.0
    del1::T   = 0.0
    del2::T   = 0.0
    del3::T   = 0.0
    fasx2::T  = 0.0
    fasx4::T  = 0.0
    fasx6::T  = 0.0
    d2201::T  = 0.0
    d2211::T  = 0.0
    d3210::T  = 0.0
    d3222::T  = 0.0
    d4410::T  = 0.0
    d4422::T  = 0.0
    d5220::T  = 0.0
    d5232::T  = 0.0
    d5421::T  = 0.0
    d5433::T  = 0.0
    xnddt::T  = 0.0
    xndot::T  = 0.0
    xldot::T  = 0.0
    zmos::T   = 0.0
    se2::T    = 0.0
    se3::T    = 0.0
    si2::T    = 0.0
    si3::T    = 0.0
    sl2::T    = 0.0
    sl3::T    = 0.0
    sl4::T    = 0.0
    sgh2::T   = 0.0
    sgh3::T   = 0.0
    sgh4::T   = 0.0
    sh2::T    = 0.0
    sh3::T    = 0.0
    zmol::T   = 0.0
    ee2::T    = 0.0
    e3::T     = 0.0
    xi2::T    = 0.0
    xi3::T    = 0.0
    xl2::T    = 0.0
    xl3::T    = 0.0
    xl4::T    = 0.0
    xgh2::T   = 0.0
    xgh3::T   = 0.0
    xgh4::T   = 0.0
    xh2::T    = 0.0
    xh3::T    = 0.0
    pe::T     = 0.0
    pinc::T   = 0.0
    pgh::T    = 0.0
    ph::T     = 0.0
    pl::T     = 0.0
    pgh0::T   = 0.0
    ph0::T    = 0.0
    pe0::T    = 0.0
    pinc0::T  = 0.0
    pl0::T    = 0.0

    isynfl::Bool = false
    iresfl::Bool = false
    ilsz::Bool   = false
end

"""
Low level SGP4 structure.
"""
@with_kw mutable struct SGP4_Structure{T<:Real}
    # TLE parameters.
    epoch::T
    n_0::T
    e_0::T
    i_0::T
    Ω_0::T
    ω_0::T
    M_0::T
    bstar::T
    # Propagation time from epoch.
    Δt::T
    # Current mean orbital elements.
    a_k::T
    e_k::T
    i_k::T
    Ω_k::T
    ω_k::T
    M_k::T
    n_k::T
    # Parameters related with the orbit.
    all_0::T
    nll_0::T
    # Useful constants to decrease the computational burden.
    AE::T
    QOMS2T::T
    β_0::T
    ξ::T
    η::T
    sin_i_0::T
    θ::T
    θ²::T
    A_30::T
    k_2::T
    k_4::T
    C1::T
    C3::T
    C4::T
    C5::T
    D2::T
    D3::T
    D4::T
    dotM::T
    dotω::T
    dotΩ1::T
    dotΩ::T
    # Selected algorithm.
    algorithm::Symbol
    # SGP4 gravitational constants.
    sgp4_gc::SGP4_GravCte{T}
    # SGP4 deep space structure.
    sgp4_ds::SGP4_DeepSpace{T}
end


text:
module FiniteDifferences

using LinearAlgebra
using SparseArrays

import Base: eltype, show

struct Basis{T}
    j::UnitRange{Integer}
    ρ::T
    Z::T
    δβ₁::T # Correction used for bare Coulomb potentials, Eq. (22) Schafer2009
end

function Basis(n::I, ρ::T, j₀::I, Z::T=one(T)) where {I<:Integer, T}
    j₀ < zero(I) && throw(ArgumentError("Illegal starting point of grid"))
    j = (1:n) .+ j₀
    Basis{T}(j, ρ, Z, j₀ == 0 ? (Z*ρ/8 * (one(T) + Z*ρ)) : zero(T))
end

Basis(n::I, ρ::T, Z::T=one(T), r₀::T=zero(T)) where {I<:Integer, T} =
    Basis(n, ρ, floor(Int, r₀/ρ), Z)

# Basis(rₘₐₓ::T, n::I, Z::T=one(T)) where {I<:Integer,T} =
#     Basis{T}(n, rₘₐₓ/(n-1/2), Z)

eltype(::Basis{T}) where T = T

locs(basis::Basis{T}) where T = (basis.j .- 1/2)*basis.ρ

function (basis::Basis)(D::Diagonal)
    n = length(basis.j)
    @assert size(D) == (n,n)
    D
end
(basis::Basis)(f::Function) = Diagonal(f.(locs(basis)))
(basis::Basis)(::UniformScaling) = I

# Not really useful, only implemented to be compatible with FEDVR and
# BSplines.
function (basis::Basis)(r::AbstractVector)
    B = spzeros(eltype(r), length(r), length(basis.j))
    for j = basis.j[1:end-1]
        a = (j-1/2) * basis.ρ
        b = (j+1/2) * basis.ρ
        B[:,j] = (a .<= r .< b)
    end
    B
end

# Variationally derived coefficients for three-point stencil
α(j::Integer) = j^2/(j^2 - 1/4)
β(j::Integer) = (j^2 - j + 1/2)/(j^2 - j + 1/4)

function β(basis::Basis)
    b = β.(basis.j)
    b[1] += basis.δβ₁
    b
end

function derop(::Type{U}, basis::Basis{T}, o::Integer) where {U,T}
    o ∉ 0:2 && error("Unsupported derivative order $(o) ∉ 0:2")
    if o == 0
        basis(I)
    elseif o == 1
        a = α.(basis.j[1:end-1])
        Tridiagonal{U}(-a, zeros(T, length(basis.j)), a)/2basis.ρ
    else
        SymTridiagonal{U}(-2β(basis), α.(basis.j[1:end-1]))/basis.ρ^2
    end
end
derop(basis::Basis{T}, o::Integer) where T = derop(T, basis, o)

show(io::IO, basis::Basis) =
    write(io, "Finite differences basis with $(length(basis.j)) points spaced by ρ = $(basis.ρ)")

end # module


text:
# Generate random seeds a priori for reproducibility

using Distributions

# random seeds for Monte Carlo estimation of optimality gap
srand(9518)
const randomSeeds_MC = ceil.(Int64,rand(Uniform(0,10000),maxNumMCReplicates))

# random seeds for constructing different model replicates
srand(8543)
const randomSeeds_models = ceil.(Int64,rand(Uniform(0,10000),maxNumModelReplicates))

# random seeds for constructing different data replicates for a given model
srand(7654)
const randomSeeds_data = ceil.(Int64,rand(Uniform(0,10000),maxNumDataReplicates))

max_stars_repo_name: rohitkannan/DD-SAA

max_issues_repo_name: rohitkannan/DD-SAA

max_forks_repo_name: rohitkannan/DD-SAA


text:
# https://adventofcode.com/2021/day/9
using AdventOfCode
using StatsBase
CI=CartesianIndex
input = split.(readlines("2021/data/day_9.txt"),"") .|> (x->parse.(Int,x)) |> Base.Fix1(reduce, hcat)

function part_1(input)
    lows=@views min.(input[[begin+1:end;end-1],begin:end],
        input[[begin+1;begin:end-1],begin:end],
        input[begin:end,[begin+1:end;end-1]],
        input[begin:end,[begin+1;begin:end-1]]) .> input
    sum((@view input[lows]).+1)
end
function findbasin(point ,input, lowpoints)
    if input[point]==9 || (point in lowpoints) return(point) end
    CIs=CartesianIndices(input)
    fI, lI = first(CIs), last(CIs)
    neighbors=min.(Ref(lI),max.(Ref(fI), CI.([(1,0),(-1,0),(0,1),(0,-1)]).+ Ref(point)))
    filter!(!in((point,)),neighbors) # get neighbors not including self
    lowneighbor=neighbors[findmin(input[neighbors])[2]]
    findbasin(lowneighbor,input, lowpoints)
end
function part_2(input)
    lows=@views min.(input[[begin+1:end;end-1],begin:end],
        input[[begin+1;begin:end-1],begin:end],
        input[begin:end,[begin+1:end;end-1]],
        input[begin:end,[begin+1;begin:end-1]]) .> input
    lowpoints=findall(lows)
    cm=countmap(findbasin(x,input,lowpoints) for x in CartesianIndices(input))
    collect(sort(cm, byvalue=true))[end-2:end].|>last|>prod
end


text:
using DataFrames
datavector = @data(["A", "A", "A","B", "B", "B"])
pooleddatavector = @pdata(["A", "A", "A","B", "B", "B"])

levels(pooleddatavector)

pooleddatavector = compact(pooleddatavector)

dataframe_notpooled = DataFrame(A = [10, 10, 10, 20, 20, 20], B = ["X", "X", "X", "Y", "Y", "Y"])

pooleddf = pool!(dataframe_notpooled, [:A, :B])

max_stars_repo_path: Chapter 2/9.jl

max_stars_repo_name: PacktPublishing/Julia-for-data-science

max_issues_repo_path: Chapter 2/9.jl

max_issues_repo_name: PacktPublishing/Julia-for-data-science

max_forks_repo_path: Chapter 2/9.jl

max_forks_repo_name: PacktPublishing/Julia-for-data-science