text:
"""
Delete existing decision tree.
"""
function clear_tree!(p::BOMCPPlanner)
    p.tree = nothing
end

"""
Construct an BOMCP tree and choose the best action.
"""
POMDPs.action(p::BOMCPPlanner, b) = first(action_info(p, b))

"""
Construct an BOMCP tree and choose the best action. Also output some information.
"""
function POMDPModelTools.action_info(p::BOMCPPlanner, b; tree_in_info=false)
    local a::actiontype(p.p)
    info = Dict{Symbol, Any}()
    if p.gp isa GPLA
        reset_gp!(p.solver.action_selector.gp)
    end
    try
        B = typeof(b)
        A = actiontype(p.p)
        tree = BOMCPTree{B,A}(p.solver.n_iterations)
        p.tree = tree
        bnode = insert_belief_node!(tree, b)
        nquery = 0
        start_us = CPUtime_us()
        for i = 1:p.solver.n_iterations
            nquery += 1
            s = rand(p.rng, b)
            simulate(p, s, bnode, p.solver.depth)
            if p.gp isa GPLA
                update_gp_buffer!(p.solver.action_selector)
            end

            if CPUtime_us() - start_us >= p.solver.max_time * 1e6
                break
            end
        end
        backup_belief_node!(p.p, tree, 1)
        info[:search_time_us] = CPUtime_us() - start_us
        info[:tree_queries] = nquery
        if p.solver.tree_in_info || tree_in_info
            info[:tree] = tree
        end
        best_Q = -Inf
        banode = 0
        for child in tree.b_children[bnode]
            if tree.q[child] > best_Q
                best_Q = tree.q[child]
                banode = child
            end
        end
        a = tree.a_labels[banode] # choose action with highest approximate value
    catch ex
        a = convert(actiontype(p.p), default_action(p.solver.default_action, p.p, b, ex))
        info[:exception] = ex
    end
    return a, info
end

"""
Return the reward for one iteration of BOMCP.
    (p, s, bnode, p.solver.depth)
"""
function simulate(dpw::BOMCPPlanner, s, bnode::Int, d::Int)
    sol = dpw.solver
    tree = dpw.tree
    b = tree.b_labels[bnode]
    if isterminal(dpw.p, s)
        return 0.0
    elseif d == 0
        return estimate_value(dpw.solved_estimate, dpw.p, s, d)
    end

    # ========== Start action progressive widening ==========
    new_node = false
    if length(tree.b_children[bnode]) <= sol.k_action*tree.n_b[bnode]^sol.alpha_action
        a = next_action(sol.action_selector, dpw.p, b, BOMCPBeliefNode(tree, bnode)) # action generation step
        if !(a in tree.a_labels[tree.b_children[bnode]])
            n0 = init_N(sol.init_N, dpw.p, b, a)
            banode = insert_action_node!(tree, bnode, a, n0,
                                init_Q(sol.init_Q, dpw.p, b, a)
                               )
            tree.n_b[bnode] += n0
            new_node = true
        end
    end
    if !new_node
        best_UCB = -Inf
        banode = 0
        ltn = log(tree.n_b[bnode])
        for child in tree.b_children[bnode]
            n = tree.n_a[child]
            q = tree.q[child]
            c = sol.exploration_constant # for clarity
            UCB = q + c*sqrt(ltn/n)
            @assert !isnan(UCB) "UCB was NaN (q=$q, c=$c, ltn=$ltn, n=$n)"
            @assert !isequal(UCB, -Inf)
            if UCB > best_UCB
                best_UCB = UCB
                banode = child
                a = tree.a_labels[banode]
            end
        end
    end
    # ========== End action progressive widening ==========
    # ========== Start belief progressive widening ==========
    sp, o, r = POMDPs.gen(DDNOut(:sp, :o, :r), dpw.p, s, a, dpw.rng)
    bp = POMDPs.update(sol.belief_updater, b, a, o)
    if haskey(tree.transitions, (banode, bp))
        bpnode, n = tree.transitions[(banode, bp)]
        n += 1
        tree.transitions[(banode, bp)] = (bpnode, n)
        tree.m_b[bpnode] += 1
    elseif length(tree.a_children[banode]) <= sol.k_belief*tree.n_a[banode]^sol.alpha_belief
        bpnode = insert_belief_node!(tree, bp)
        n = 1
        tree.transitions[(banode, bp)] = (bpnode, n)
        push!(tree.a_children[banode], bpnode)
        new_node = true
    else
        children = tree.a_children[banode]
        weights = Weights(tree.m_b[children]/sum(tree.m_b[children]))
        bpnode = sample(dpw.rng, children, weights)
        bp = tree.b_labels[bpnode]
        bpnode, n = tree.transitions[(banode, bp)]
        sp = rand(bp)
        r = POMDPs.reward(dpw.p, s, a, sp)
    end
    # ========== End belief progressive widening ==========
    if new_node
        v = estimate_value(dpw.solved_estimate, dpw.p, sp, bp, d-1)
        tree.v[bpnode] = v
    else
        v = simulate(dpw, sp, bpnode, d-1)
    end

    tree.n_a[banode] += 1
    tree.n_b[bnode] += 1

    q = r + discount(dpw.p)*v
    tree.q[banode] += (q - tree.q[banode])/tree.n_a[banode]
    tree.r[banode] += (r - tree.r[banode])/tree.n_a[banode]

    tree.v[bnode] = maximum(tree.q[tree.b_children[bnode]])
    if sol.action_selector.gp isa GPLA
        update_gp!(sol.action_selector.gp, tree, bnode, new_node, sol.action_selector.belief_λ, sol.action_selector.obs_buffer) #, Val(sol.action_selector.gp.dim), Val(size(sol.action_selector.gp.x, 2)))
    else
        update_gp!(sol.action_selector.gp, tree, bnode, new_node, sol.action_selector.belief_λ)
    end
    return q
end

function backup_belief_node!(p::Union{POMDP, MDP}, tree::BOMCPTree, index::Int)
    children = tree.b_children[index]
    n_children = length(children)
    if n_children == 0
        return tree.v[index]
    else
        max_q = -Inf
        for i = 1:n_children
            child_id = children[i]
            q = backup_action_node!(p, tree, child_id)
            if q > max_q
                max_q = q
            end
        end
        tree.v[index] = max_q
        return max_q
    end
end

function backup_action_node!(p::Union{POMDP, MDP}, tree::BOMCPTree, index::Int)
    children = tree.a_children[index]
    n_children = length(children)
    q = 0.0
    m_total = 0.0
    for i = 1:n_children
        child_id = children[i]
        m = tree.m_b[child_id]
        q += backup_belief_node!(p, tree, child_id)*m
        m_total += m
    end
    q /= m_total
    q *= POMDPs.discount(p)
    q += tree.r[index]
    tree.q[index] = q
    return q
end


text:
@testset "Continuous identity system" begin
    for sd in 1:3
        s = ContinuousIdentitySystem(sd)
        @test statedim(s) == sd
        @test inputdim(s) == 0
        @test noisedim(s) == 0
        @test islinear(s) && isaffine(s) && !ispolynomial(s)
        @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s)
    end
end

@testset "Continuous constrained identity system" begin
    for sd in 1:3
        X = Singleton(ones(sd))
        s = ConstrainedContinuousIdentitySystem(sd, X)
        @test statedim(s) == sd
        @test inputdim(s) == 0
        @test noisedim(s) == 0
        @test stateset(s) == X
        @test islinear(s) && isaffine(s) && !ispolynomial(s)
        @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s)
    end
end

@testset "Continuous linear system" begin
    for sd in 1:3
        s = LinearContinuousSystem(zeros(sd, sd))
        @test statedim(s) == sd
        @test inputdim(s) == 0
        @test noisedim(s) == 0
        @test islinear(s) && isaffine(s) && !ispolynomial(s)
        @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s)
    end
end

@testset "Continuous affine system" begin
    for sd in 1:3
        s = AffineContinuousSystem(zeros(sd, sd), zeros(sd))
        @test statedim(s) == sd
        @test inputdim(s) == 0
        @test !islinear(s) && isaffine(s) && !ispolynomial(s)
        @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s)
    end
end

@testset "Continuous linear control system" begin
    for sd in 1:3
        s = LinearControlContinuousSystem(zeros(sd, sd), ones(sd, sd))
        @test statedim(s) == sd
        @test inputdim(s) == sd
        @test noisedim(s) == 0
        @test islinear(s) && isaffine(s) && !ispolynomial(s)
        @test !isnoisy(s) && iscontrolled(s) && !isconstrained(s)
    end
end

@testset "Continuous constrained linear system" begin
    A = [1. 1; 1 -1]
    X = Line([1., -1], 0.) # line x = y
    s = ConstrainedLinearContinuousSystem(A, X)
    @test statedim(s) == 2
    @test inputdim(s) == 0
    @test noisedim(s) == 0
    @test stateset(s) == X
    @test islinear(s) && isaffine(s) && !ispolynomial(s)
    @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s)
end

@testset "Continuous constrained affine system" begin
    A = [1. 1; 1 -1]
    c = [1.; 1.]
    X = Line([1., -1], 0.) # line x = y
    s = ConstrainedAffineContinuousSystem(A, c, X)
    @test statedim(s) == 2
    @test inputdim(s) == 0
    @test noisedim(s) == 0
    @test stateset(s) == X
    @test !islinear(s) && isaffine(s) && !ispolynomial(s)
    @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s)
end

@testset "Continuous affine control system with state constraints" begin
    A = zeros(2, 2)
    B = Matrix([0.5 1.5]')
    c = [0.0, 1.0]
    X = Line([1., -1], 0.) # line x = y
    U = Interval(-1.0, 1.0)  # -1 <= u <= 1
    s = ConstrainedAffineControlContinuousSystem(A, B, c, X, U)
    @test statedim(s) == 2
    @test inputdim(s) == 1
    @test noisedim(s) == 0
    @test stateset(s) == X
    @test inputset(s) == U
    @test !islinear(s) && isaffine(s) && !ispolynomial(s)
    @test !isnoisy(s) && iscontrolled(s) && isconstrained(s)
end

@testset "Continuous constrained linear control system" begin
    A = [1. 1; 1 -1]
    B = Matrix([0.5 1.5]')
    X = Line([1., -1], 0.)
    U = Interval(0.9, 1.1)
    s = ConstrainedLinearControlContinuousSystem(A, B, X, U)
    @test statedim(s) == 2
    @test inputdim(s) == 1
    @test noisedim(s) == 0
    @test stateset(s) == X
    @test inputset(s) == U
    @test islinear(s) && isaffine(s) && !ispolynomial(s)
    @test !isnoisy(s) && iscontrolled(s) && isconstrained(s)

    # initial value problem composite type
    x0 = Singleton([1.5, 2.0])
    p = InitialValueProblem(s, x0)
    @test inputset(p) == U
    @test statedim(p) == 2
    @test inputdim(p) == 1

    # check that the matrices A and B need not be of the same type
    # with A dense and B a lazy adjoint
    B = [0.5 1.5]'
    s = ConstrainedLinearControlContinuousSystem(A, B, X, U)

    # check that the matrices A and B need not be of the same type
    # with A sparse and B a lazy adjoint
    A = sparse([1], [2], [1.0], 2, 2) # sparse matrix
    s = ConstrainedLinearControlContinuousSystem(A, B, X, U)
end

@testset "Continuous linear algebraic system" begin
    for sd in 1:3
        s = LinearAlgebraicContinuousSystem(zeros(sd, sd), zeros(sd, sd))
        @test statedim(s) == sd
        @test inputdim(s) == 0
        @test noisedim(s) == 0
        @test islinear(s) && isaffine(s) && !ispolynomial(s)
        @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s)
    end
end

@testset "Continuous constrained linear algebraic system" begin
    A = [1. 1; 1 -1]
    E = [0. 1; 1 0]
    X = LinearConstraint([0, -1.], 0.) # the set y ≥ 0
    s = ConstrainedLinearAlgebraicContinuousSystem(A, E, X)
    @test statedim(s) == 2
    @test inputdim(s) == 0
    @test noisedim(s) == 0
    @test stateset(s) == X
    @test islinear(s) && isaffine(s) && !ispolynomial(s)
    @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s)
end

@testset "Initial value problem for a continuous constrained linear algebraic system" begin
    A = [1. 1; 1 -1]
    E = [0. 1; 1 0]
    X = LinearConstraint([0, -1.], 0.) # the set y ≥ 0
    s = ConstrainedLinearAlgebraicContinuousSystem(A, E, X)
    @test statedim(s) == 2
    @test inputdim(s) == 0
    @test noisedim(s) == 0
    @test stateset(s) == X
    @test islinear(s) && isaffine(s) && !ispolynomial(s)
    @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s)

    x0 = Singleton([1.5, 2.0])
    p = IVP(s, x0)
end

@testset "Polynomial system in continuous time" begin
    @polyvar x y
    p = 2x^2 - 3x + y

    # default constructor for scalar p and
    s = PolynomialContinuousSystem(p)
    @test statedim(s) == 2
    @test inputdim(s) == 0
    @test noisedim(s) == 0
    @test !islinear(s) && !isaffine(s) && ispolynomial(s)
    @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s)

    @test TypedPolynomials.nvariables(s) == 2
    @test TypedPolynomials.variables(s) == (x, y)

    # constructor for scalar p and given number of variables
    s = PolynomialContinuousSystem([p], 2)
end

@testset "Polynomial system in continuous time with state constraints" begin
    @polyvar x y
    p = 2x^2 - 3x + y
    X = BallInf(zeros(2), 0.1)

    # default constructor for scalar p and
    s = ConstrainedPolynomialContinuousSystem(p, X)
    @test statedim(s) == 2
    @test inputdim(s) == 0
    @test noisedim(s) == 0
    @test dim(stateset(s)) == dim(X)
    @test !islinear(s) && !isaffine(s) && ispolynomial(s)
    @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s)

    @test TypedPolynomials.nvariables(s) == 2
    @test TypedPolynomials.variables(s) == (x, y)
end

function vanderpol!(x, dx)
    dx[1] = x[2]
    dx[2] = x[2] * (1-x[1]^2) - x[1]
    return dx
end

@testset "Continuous system defined by a function" begin
    s = BlackBoxContinuousSystem(vanderpol!, 2)
    f! = (x, dx) -> s.f(x, dx)
    x = [1.0, 0.0]
    dx = similar(x)
    f!(x, dx)
    @test dx ≈ [0.0, -1.0]
    @test !islinear(s) && !isaffine(s) && !ispolynomial(s)
    @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s)
end

@testset "Continuous system defined by a function with state constraints" begin
    H = HalfSpace([1.0, 0.0], 0.0) # x <= 0
    s = ConstrainedBlackBoxContinuousSystem(vanderpol!, 2, H)
    f! = (x, dx) -> s.f(x, dx)
    x = [1.0, 0.0]
    dx = similar(x)
    f!(x, dx)
    @test dx ≈ [0.0, -1.0]
    @test stateset(s) == H
    @test !islinear(s) && !isaffine(s) && !ispolynomial(s)
    @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s)
end

function vanderpol_controlled!(x, u, dx)
    dx[1] = x[2]
    dx[2] = x[2] * (1-x[1]^2) - x[1] + u[1]
    return dx
end

@testset "Continuous control system defined by a function with state constraints" begin
    H = HalfSpace([1.0, 0.0], 0.0) # x <= 0
    U = Interval(-0.1, 0.1)
    s = ConstrainedBlackBoxControlContinuousSystem(vanderpol_controlled!, 2, 1, H, U)
    f! = (x, u, dx) -> s.f(x, u, dx)
    x = [1.0, 0.0]
    u = an_element(U)
    dx = similar(x)
    f!(x, u, dx)
    @test dx ≈ [0.0, -1.0 + u[1]]
    @test statedim(s) == 2
    @test inputdim(s) == 1
    @test stateset(s) == H
    @test inputset(s) == U
    @test !islinear(s) && !isaffine(s) && !ispolynomial(s)
    @test !isnoisy(s) && iscontrolled(s) && isconstrained(s)
end

# Noisy
@testset "Noisy Continuous constrained linear system" begin
    A = [1. 1; 1 -1]
    D = [1. 2; 0 1]
    X = Line([1., -1], 0.) # line x = y
    W = BallInf(zeros(2), 2.0)
    s = NoisyConstrainedLinearContinuousSystem(A, D, X, W)
    @test s.A == A
    @test s.D == D
    @test statedim(s) == 2
    @test inputdim(s) == 0
    @test noisedim(s) == 2 == dim(W)
    @test stateset(s) == X
    @test noiseset(s) == W
    @test islinear(s) && isaffine(s) && !ispolynomial(s)
    @test isnoisy(s) && !iscontrolled(s) && isconstrained(s)
end

@testset "Noisy Continuous constrained control linear system" begin
    A = [1. 1; 1 -1]
    B = Matrix([0.5 1.5]')
    D = [1. 2; 0 1]
    X = Line([1., -1], 0.) # line x = y
    U = Hyperrectangle(low=[0.9], high=[1.1])
    W = BallInf(zeros(2), 2.0)
    s = NoisyConstrainedLinearControlContinuousSystem(A, B, D, X, U, W)
    @test s.A == A
    @test s.B == B
    @test s.D == D
    @test statedim(s) == 2
    @test inputdim(s) == dim(U)
    @test noisedim(s) == 2 == dim(W)
    @test stateset(s) == X
    @test inputset(s) == U
    @test noiseset(s) == W
    @test islinear(s) && isaffine(s) && !ispolynomial(s)
    @test isnoisy(s) && iscontrolled(s) && isconstrained(s)
end

@testset "Noisy Continuous constrained control affine system" begin
    A = [1. 1; 1 -1]
    B = Matrix([0.5 1.5]')
    c = [1.0, 0.5]
    D = [1. 2; 0 1]
    X = Line([1., -1], 0.) # line x = y
    U = Hyperrectangle(low=[0.9], high=[1.1])
    W = BallInf(zeros(2), 2.0)
    s = NoisyConstrainedAffineControlContinuousSystem(A, B, c, D, X, U, W)
    @test s.A == A
    @test s.B == B
    @test s.c == c
    @test s.D == D
    @test statedim(s) == 2
    @test inputdim(s) == dim(U)
    @test noisedim(s) == 2 == dim(W)
    @test stateset(s) == X
    @test inputset(s) == U
    @test noiseset(s) == W
    @test isaffine(s) && !islinear(s) && !ispolynomial(s)
    @test isnoisy(s) && iscontrolled(s) && isconstrained(s)
end

@testset "Noisy Continuous constrained control blackbox system" begin
    n = 2
    m = 1
    l = 2
    f(x,u,w) = ones(n,n)*x + ones(n,m)*u + ones(n,l)*w
    X = BallInf(zeros(n), 1.0)
    U =  BallInf(zeros(m), 1.0)
    W =  BallInf(zeros(l), 1.0)
    s = NoisyConstrainedBlackBoxControlContinuousSystem(f, n, m, l, X, U, W)
    @test s.f == f
    @test statedim(s) == n
    @test inputdim(s) == dim(U)
    @test noisedim(s) == l == dim(W)
    @test stateset(s) == X
    @test inputset(s) == U
    @test noiseset(s) == W
    @test !islinear(s) && !isaffine(s) && !ispolynomial(s)
    @test isnoisy(s) && iscontrolled(s) && isconstrained(s)
end

ext: jl


text:
import Base: size, ndims, length, show, nameof

"""
    AbstractLattice{N}

Abstract type for general lattices. `N` indicates the dimension.
For a more concrete definition please refer the following material:

- Lattice (group): https://en.wikipedia.org/wiki/Lattice_(group)
- Lattice Model: https://en.wikipedia.org/wiki/Lattice_model_(physics)
- Lattice Graph: https://en.wikipedia.org/wiki/Lattice_graph
"""
abstract type AbstractLattice{N} end

"""
    ndims(lattice) -> N

Returns number of dimensions of the lattice.
"""
ndims(::AbstractLattice{N}) where N = N

"""
    size(lattice) -> Tuple

Returns the size of this `lattice`.
"""
function size(::AbstractLattice)
    error("Not Implemented")
end

"""
    name(lattice) -> String

Returns the name of this lattice.
"""
nameof(ltc::AbstractLattice) = "Abstract Lattice"

"""
    length(lattice) -> Int

Returns the length of this lattice, or the production
of each dimension size.
"""
function length(::AbstractLattice)
    error("Not Implemented")
end


text:
include("Graph.jl")

using SparseArrays
using LinearAlgebra
import Laplacians

function RoleSim(G; delta = 1e-6, beta = 0.15, maxIt = 100)
    n = G.n
    g = adjacentList(G)
    d = degrees(G)
    Spre = zeros(n, n)
    S = zeros(n, n)
    Sp = ones(n, n)

    calc(x, y) = begin
        hx = zeros(Bool, d[x])
        hy = zeros(Bool, d[y])
        A = Vector{Tuple{Int32, Int32, Float64}}()
        for i = 1 : d[x], j = 1 : d[y]
            push!(A, (i, j, Sp[g[x][i], g[y][j]]))
        end
        sort!(A, by=x->-x[3])
        ret = 0.0
        for (u, v, w) in A
            if !hx[u] && !hy[v]
                hx[u] = true
                hy[v] = true
                ret += w
            end
        end
        return ret
    end

    ts = 0
    while true
        ts += 1
        for i = 1 : n, j = 1 : n
            S[i, j] = (1-beta)*calc(i, j)/max(d[i], d[j]) + beta
        end
        if (norm(S - Sp) < delta) || (ts > maxIt)
            break
        end
        for i = 1 : n, j = 1 : n
            Sp[i, j] = S[i, j]
        end
    end
    return S
end

function getW(G)
    IpL = zeros(G.n, G.n)
    for (u, v) in G.E
        IpL[u, u] += 1
        IpL[v, v] += 1
        IpL[u, v] -= 1
        IpL[v, u] -= 1
    end
    foreach(i -> IpL[i, i] += 1, 1 : G.n)
    return inv(IpL)
end

function ForestSim(G)
    n = G.n
    W = getW(G)
    S = zeros(n, n)
    for i = 1 : n, j = 1 : n
        S[i, j] = min(W[i, i], W[j, j]) / max(W[i, i], W[j, j])
    end
    return S
end

function getSparseIpL(G)
    d = ones(G.n)
    for (u, v) in G.E
        d[u] += 1
        d[v] += 1
    end
    Is = zeros(Int32, G.m*2+G.n)
    Js = zeros(Int32, G.m*2+G.n)
    Vs = zeros(G.m*2+G.n)
    ID = 0
    for (u, v) in G.E
        ID += 1
        Is[ID] = u
        Js[ID] = v
        Vs[ID] = -1
        Is[ID + G.m] = v
        Js[ID + G.m] = u
        Vs[ID + G.m] = -1
    end
    for i = 1 : G.n
        Is[G.m + G.m + i] = i
        Js[G.m + G.m + i] = i
        Vs[G.m + G.m + i] = d[i]
    end
    return sparse(Is, Js, Vs, G.n, G.n)
end

function getSparseB(G)
    n = G.n
    m = G.m
    Is = zeros(Int32, m*2)
    Js = zeros(Int32, m*2)
    Vs = zeros(m*2)
    ID = 0
    for (u, v) in G.E
        ID += 1
        Is[ID] = ID
        Js[ID] = u
        Vs[ID] = 1
        Is[ID + m] = ID
        Js[ID + m] = v
        Vs[ID + m] = -1
    end
    return sparse(Is, Js, Vs, m, n)
end

function approxForest(G; eps = 0.1)
    n = G.n
    m = G.m
    IpL = getSparseIpL(G)
    B = getSparseB(G)
    k = round(Int, log2(G.n) / eps^2)
    f = Laplacians.approxchol_sddm(IpL, tol=1e-12)
    ret = zeros(G.n)
    for i = 1 : k
        y1 = B' * randn(m)
        y2 = rand(n)
        z1 = f(y1)
        z2 = f(y2)
        foreach(j -> ret[j] += (z1[j]^2 + z2[j]^2), 1 : G.n)
    end
    ret ./= k
    return ret
end

function calcBins(G, maxk)
    maxbin = ceil(Int32, log2(G.n))
    n = G.n
    g = adjacentList(G)
    d = degrees(G)
    B = zeros(Int32, maxk, n, maxbin)
    dk = zeros(Int32, maxk, n)
    for s = 1 : n
        dist = zeros(Int32, n)
        fill!(dist, 1048576)
        dist[s] = 0
        Q = zeros(Int32, n)
        front = 1
        rear = 1
        Q[1] = s
        while front <= rear
            u = Q[front]
            front += 1
            if dist[u] == maxk
                continue
            end
            for v in g[u]
                if dist[u]+1 < dist[v]
                    dist[v] = dist[u] + 1
                    B[dist[v], s, floor(Int32, log2(d[v]))+1] += 1
                    dk[dist[v], s] += 1
                    rear += 1
                    Q[rear] = v
                end
            end
        end
    end
    return B, dk
end

function StructSim(G; maxk = 3, wi = 0.5)
    n = G.n
    d = degrees(G)
    B, dk = calcBins(G, maxk)
    maxb = size(B, 3)
    S = zeros(n, n)
    for i = 1 : n, j = 1 : n
        S[i, j] = min(d[i], d[j]) / max(d[i], d[j])
    end
    for k = 1 : maxk
        for i = 1 : n, j = 1 : n
            delta = 0.0
            foreach(z -> delta += min(B[k, i, z], B[k, j, z]), 1 : maxb)
            delta /= max(dk[k, i], dk[k, j])
            S[i, j] = (1.0 - wi) * S[i, j] + wi * delta
        end
    end
    return S
end


text:
mutable struct PowerSimulationsModel{T<:AbstractOperationsModel, R<:Dates.Period}
    name::String
    model::PowerOperationModel{T}
    steps::Int64
    periods::Int64
    resolution::R
    date_from::Dates.DateTime
    date_to::Dates.DateTime
    lookahead_periods::Int64
    lookahead_resolution::R
    dynamic_analysis::Bool
    timeseries::Dict{Any,Any}
end


text:
using DecFP
using DecFP: Dec128, _parse

macro define(name, definition)
    quote
        macro $(esc(name))()
            esc($(Expr(:quote, definition)))
        end
    end
end

function _dec128(x::Expr)
    y = x
    y.args .= _dec128.(y.args)
    return  y
end

_dec128(x::Int) = Dec128(x)
_dec128(x::String) = _parse(Dec128, x)

_dec128(x) = x

macro dec128(x)
    return esc(_dec128(x))
end


text:
using Documenter, StringBuilders

makedocs(
	modules = [StringBuilders],
	sitename = "StringBuilders.jl",
	pages = [
        "Home" => "index.md"
    ]
)

deploydocs(
    repo = "github.com/davidanthoff/StringBuilders.jl.git"
)


text:
using Graphs
using CustomVertexGraphs

g = set_vertices(grid((4,)), ["A", "B", "C", "D"])
@show g
@show g[["A", "B"]]
@show has_edge(g, "A" => "B")
@show add_edge!(g, "A" => "C")
@show has_edge(g, "A" => "C")


text:
"""
    module Lattice

for defining discrete Cartesian and Polar lattice grids used in Iris.
"""
module Lattices

using ...Defaults
using Formatting
using RecipesBase
using StaticArrays

export Lattice,
Cartesian,
Polar,
get_lattice_index,
LatticeError

struct Lattice
    type::Symbol

    dx::Float64
    dy::Float64
    θ::Float64
    x0::Float64
    y0::Float64

    dr::Float64
    dθ::Float64
    r0::Float64

    sinθ::Float64 # to avoid inefficient trig calls later
    cosθ::Float64 # to avoid inefficient trig calls later
    v1::SArray{Tuple{2},Float64,1,2}
    v2::SArray{Tuple{2},Float64,1,2}
    origin::SArray{Tuple{2},Float64,1,2}

    function Lattice(type,args...)
        if type==:Cartesian
            CartesianLat(args...)
        else
            PolarLat(args...)
        end
    end

    function CartesianLat(dx::Number,dy::Number,θ::Number=0,x0::Number=0,y0::Number=0)
        sinθ, cosθ = sincos(θ)
        R = @SMatrix [ cos(θ) -sin(θ); sin(θ)  cos(θ)]
        v1 = R*SVector{2}([1,0])
        v2 = R*SVector{2}([0,1])
        origin = SVector{2}([x0,y0])
        return new(:Cartesian,dx, dy, θ, x0, y0, NaN, NaN, NaN, sinθ, cosθ, v1, v2, origin)
    end

    function PolarLat(dr::Number,dθ::Number,θ::Number=0,x0::Number=0,y0::Number=0,r0::Number=dr/2)
        return new(:Polar,NaN,NaN,θ,x0,y0,dr,dθ,r0)
    end

    # pretty printing of Lattice
    function Base.show(io::IO, lat::Lattice)
        if lat.type==:Cartesian
            if !get(io, :sub, false)
                print(io, "Cartesian Lattice", ": \n",
                "\tdx=", fmt("3.2f",lat.dx), ", dy=", fmt("3.2f",lat.dy), ", ∠", fmt("3.2f",(mod2pi(lat.θ+π/2))*180/π), "°\n",
                "\torigin: (", fmt("3.2f",lat.x0), ", ", fmt("3.2f",lat.y0), ")")
            elseif !get(io, :sub1, false)
                print(io,
                "\n\tdx=", fmt("3.2f",lat.dx), ", dy=", fmt("3.2f",lat.dy), ", ∠", fmt("3.2f",(mod2pi(lat.θ+π/2))*180/π), "°\n",
                "\torigin: (", fmt("3.2f",lat.x0), ", ", fmt("3.2f",lat.y0), ")")
            else
                print(io,
                "\n\t\tdx=", fmt("3.2f",lat.dx),", dy=", fmt("3.2f",lat.dy), ", ∠", fmt("3.2f",(mod2pi(lat.θ+π))*180/π), "°\n",
                "\t\torigin: (", fmt("3.2f",lat.x0), ", ", fmt("3.2f",lat.y0), ")")
            end
        elseif lat.type==:Polar
            if !get(io, :sub, false)
                print(io, "Polar Lattice", ": \n",
                "\tdr=", fmt("3.2f",lat.dr), ", dθ=", fmt("3.2f",lat.dθ), "=π/", fmt("3.1f",π/lat.dθ),", ∠", fmt("3.2f",(mod2pi(lat.θ))*180/π), "°\n",
                "\torigin: (", fmt("3.2f",lat.x0), ", ", fmt("3.2f",lat.y0), ")")
            elseif !get(io, :sub1, false)
                print(io,
                "\n\tdr=", fmt("3.2f",lat.dr), ", dθ=", fmt("3.2f",lat.dθ), "=π/", fmt("3.1f",π/lat.dθ),", ∠", fmt("3.2f",(mod2pi(lat.θ))*180/π), "°\n",
                "\torigin: (", fmt("3.2f",lat.x0), ", ", fmt("3.2f",lat.y0), ")")
            else
                print(io,
                "\n\t\tdr=", fmt("3.2f",lat.dr), ", dθ=", fmt("3.2f",lat.dθ), "=π/", fmt("3.1f",π/lat.dθ),", ∠", fmt("3.2f",(mod2pi(lat.θ))*180/π), "°\n",
                "\t\torigin: (", fmt("3.2f",lat.x0), ", ", fmt("3.2f",lat.y0), ")")
            end
        else
            println(io, "unknown lattice type")
        end
    end
end
"""
    Lattices.Cartesian(dx,dy,θ=0,x0=0,y0=0) -> lattice

`x0,y0` is the origin, `θ` is the rotation

----
    (::Cartesian)(i,j) -> x,y

`x,y` are coordinates of the `i,j` lattice site (`i,j` continuous)
--------
    (lat::Cartesian)(dr=lat.dr,...)
"""
Cartesian(args...) = Lattice(:Cartesian,args...)

"""
    Lattices.Polar(dr,dθ,θ=0,r0=dr/2,x0=0,y0=0) -> lattice

`r0` is the minimum radius, `x0,y0` is the origin, `θ` is the zero angle

----
    (::Polar)(i,j) -> x,y

`x,y` are coordinates of the `i,j` lattice site (`i,j` continuous)

--------
    (lat::Polar)(dr=lat.dr,...)
"""
Polar(args...) = Lattice(:Polar,args...)

"""
    (::Lattice)(i,j) -> x,y
"""
function (lat::Lattice)(i::Real,j::Real)
    if lat.type==:Cartesian
        v = lat.origin + i*lat.v1*lat.dx + j*lat.v2*lat.dy
        return v[1],v[2]
    elseif lat.type==:Polar
        r = lat.r0 + i*lat.dr
        θ = lat.θ  + j*lat.dθ
        v = r.*sincos(θ)
        return lat.x0+v[2],lat.y0+v[1]
    end
end

# functor wrappers for both Cartesian and Polar lattices
(lat::Lattice)(ijk::Tuple) = lat(ijk...)
(lat::Lattice)(k::CartesianIndex) = lat(k[1],k[2])
function (lat::Lattice)(i::Real,j::Real,v::AbstractArray)
    v[:] = lat(i,j)
    return nothing
end
function (lat::Lattice)(x::AbstractArray,y::AbstractArray)
    X = Array{Float64}(undef,size(x))
    Y = Array{Float64}(undef,size(y))
    for i ∈ eachindex(X)
        X[i], Y[i] = lat(x[i],y[i])
    end
    return X,Y
end

function (lat::Lattice)(;dx=lat.dx, dy=lat.dy, θ=lat.θ, x0=lat.x0, y0=lat.y0, dr=lat.dr, dθ=lat.dθ, r0=dr/2)
    if lat.type==:Cartesian
        return Cartesian(dx,dy,θ,x0,y0)
    else
        return Polar(dr,dθ,θ,x0,y0,r0)
    end
end

struct LatticeError <: Exception
    lat::Lattice
end
function Base.showerror(io::IO,err::LatticeError)
    print(io,"LatticeError: unrecognized lattice type ")
    print(io, err.lat.type)
    print(io, ", must be one of :Cartesian or :Polar")
end

"""
    get_lattice_index(lattice,x,y) -> nx,ny

returns indices `nx`,`ny` such that `(x,y) = (x0,y0) + nx*v1*dx + ny*v2*dx`
"""
function get_lattice_index(lat::Lattice,x::AbstractArray,y::AbstractArray)
    @assert size(x)==size(y) "x and y must have the same size"
    X = Array{Float64}(undef,size(x)...)
    Y = Array{Float64}(undef,size(y)...)
    for k ∈ eachindex(X)
        X[k],Y[k] = get_lattice_index(lat,x[k],y[k])
    end
    return X,Y
end
function get_lattice_index(lat::Lattice,x::Real,y::Real)
    i,j = lat.type==:Cartesian ? get_lattice_index_cartesian(lat,x,y) : get_lattice_index_polar(lat,x,y)
    return i,j
end
function get_lattice_index_cartesian(lat::Lattice,x::Real,y::Real)
    dx, dy, v1, v2, x0, y0 = lat.dx, lat.dy, lat.v1, lat.v2, lat.x0, lat.y0
    x += -x0
    y += -y0
    rv1 = v1[1]*x + v1[2]*y
    rv2 = v2[1]*x + v2[2]*y
    return rv1/dx, rv2/dy
end
function get_lattice_index_polar(lat::Lattice,x::Real,y::Real)
    dr, dθ, x0, y0 = lat.dr, lat.dθ, lat.x0, lat.y0
    x += -x0
    y += -y0
    r, θ = hypot(y,x), atan(y,x)+π
    return (r-lat.r0)/dr, (θ-lat.θ)/dθ
end

"""
    plot(lattice, N=[4,4])

scatter plot of lattice vectors from -`N[i]`:`+N[i]` for each dimension i.

for polar lattices, plots all angles out to radial site `N[1]`.

If `N` is a scalar, uses same span for both dimensions.
"""
@recipe f(lattice::Lattice, N::Int=4) = (lattice, [N,N])
@recipe function f(lattice::Lattice, N::Array{Int,1})
    lattice.type==:Polar ? N[2]=floor(Int,2π/lattice.dθ-eps(2π/lattice.dθ)) : nothing
    seriestype --> :scatter
    legend --> false
    aspect_ratio --> 1
    @series begin
        x, y = Float64[], Float64[]
        for n ∈ 0:N[1], m ∈ 0:N[2]
            temp = lattice(n,m)
            push!(x,temp[1]); push!(y,temp[2])
        end
        (x, y)
    end
end

end # module


text:
import Base: iterate, length

"""
    ordered_edges(dag)

Iterator of edges of a dag, ordered in Chickering order:

    Perform a topological sort on the NODES
    while there are unordered EDGES in g
        Let y be the lowest ordered NODE that has an unordered EDGE incident into it
        Let x be the highest ordered NODE for which x => y is not ordered
        return x => y
    end
"""
ordered_edges(g) = OrderedEdges(g)

struct OrderedEdges
    g
    ys::Vector{Int}
    yp::Vector{Int}
    n::Int
    m::Int

    function OrderedEdges(g::DiGraph)
        ys = topological_sort_by_dfs(g)
        yp = sortperm(ys)
        new(g, ys, yp, nv(g), ne(g))
    end
end

length(iter::OrderedEdges) = iter.m

function iterate(iter::OrderedEdges, state=(0, 0, 0, Int[]))
    i, j, k, xs  = state
    k >= iter.m && return nothing
    g, ys, yp = iter.g, iter.ys, iter.yp
    while i == 0
        j = j + 1
        xs = sort(inneighbors(g, ys[j]), by=x->yp[x])
        i = length(xs)
    end
    Edge(xs[i], ys[j]), (i - 1, j, k + 1, xs)
end

# came in handy while prototyping
function chickering_order(g)
    outp = Pair{Int,Int}[]
    ys = topological_sort_by_dfs(g)
    yp = sortperm(ys)
    n = nv(g)
    i = 0 # src index
    j = 0 # dst index
    xs = Int[]
    while true
        while i == 0
            j = j + 1
            if j > n
                @goto ende
            end
            xs = sort(inneighbors(g, ys[j]), by=x->yp[x])
            i = length(xs)
        end
        x, y = xs[i], ys[j]
        push!(outp, x=>y)
        i -= 1
    end
    @label ende
    outp
end

@enum Status reversible=0 compelled=1 unknown=2
"""
    cpdag(skel::DiGraph)

Reference: M. Chickering: Learning equivalence classes of Bayesian
network structures. Journal of Machine Learning Research 2 (2002).
M. Chickering: A Transformational Characterization of Equivalent Bayesian Network Structures. (1995).

Note that the edge order defined there is already partly encoded into the
representation of a DiGraph.
"""
function cpdag(skel::DiGraph)
    g = copy(skel)
    ys = topological_sort_by_dfs(g)
    yp = sortperm(ys)
    n = nv(g)
    m = ne(g)
    label = Dict(zip(map(Pair, edges(g)), repeated(unknown)))

    i = 0 # src index
    j = 0 # dst index
    xs = Int[]
    while true
        while i == 0
            j = j + 1
            if j > n
                @goto ende
            end
            xs = sort(inneighbors(g, ys[j]), by=x->yp[x])
            i = length(xs)
        end

        x, y = xs[i], ys[j]

        if label[x => y] == unknown
            for w in inneighbors(g, x)
                label[w => x] == compelled || continue
                if !has_edge(g, w => y)
                    for z in inneighbors(g, y)
                        label[z => y] = compelled
                    end
                    @goto next
                else
                    label[w => y] = compelled
                end
            end

            lbl = reversible
            for z in inneighbors(g, y)
                if x == z || has_edge(g, z => x)
                    continue
                end
                lbl = compelled
            end
            for z in inneighbors(g, y)
                if label[z => y] == unknown
                    label[z => y] = lbl
                end
            end
        end
        @label next

        i = i - 1
    end
    @label ende
    for e in collect(edges(g))
        x, y = Tuple(e)
        if label[x=>y] == reversible
            add_edge!(g, y => x)
        end
    end
    g
end


text:
# julia1_1.3.jl Running Julia.
println("Who are you?")
s=readline()
println("Hello $s and Hello World!")
# Bad name dictionary can be as big as one wants
# and severity of baddness can be designated by the number
badnames=Dict("Dope" => 1, "Knothead" => 1,"Stupid" => 1)
# The bot wants a name
println("Please give me a name for I have none.")
name=readline()
# Is it a bad name or a good name
n=get(badnames,name,0)
if n==0
    println("$name is nice, I will take it." )
else
    println("I think that's not a very nice name.")
    name="Sky"  #could use a random name selector here
    println("I chose the name $name for myself.")
    println("$name is a strong name and I like it.")
end
# after this you could read in sentences
# have a verb and noun dictionary
# match on the verb and the noun and do stuff
# or do the psychiatrist thing and turn their
# responses back into questions.


text:
using ComputationalHomology
using BenchmarkTools
using Random

Random.seed!(2903872398473);
rows, cols = 2, 15
Xsml = rand(rows, cols)
Xmed = rand(rows, cols*10)
Xbig = rand(rows, cols*100)
l = 20
ɛ = 0.05
Csml, wsml = vietorisrips(Xsml, ɛ, maxoutdim=1)
Cmed, wmed = vietorisrips(Xmed, ɛ, maxoutdim=1)
Cbig, wbig = vietorisrips(Xbig, ɛ, maxoutdim=1)

suite = BenchmarkGroup()

suite["parts"] = BenchmarkGroup(["components"])
suite["parts"]["position-small-start"] = @benchmarkable ComputationalHomology.position(Csml, Csml.cells[0][1])
suite["parts"]["position-small-center"] = @benchmarkable ComputationalHomology.position(Csml, Csml.cells[0][length(Csml.cells[0])>>1])
suite["parts"]["position-small-end"] = @benchmarkable ComputationalHomology.position(Csml, Csml.cells[0][length(Csml.cells[0])])
suite["parts"]["position-medium-start"] = @benchmarkable ComputationalHomology.position(Cmed, Cmed.cells[0][1])
suite["parts"]["position-medium-center"] = @benchmarkable ComputationalHomology.position(Cmed, Cmed.cells[0][length(Cmed.cells[0])>>1])
suite["parts"]["position-medium-end"] = @benchmarkable ComputationalHomology.position(Cmed, Cmed.cells[0][length(Cmed.cells[0])])
suite["parts"]["position-big-start"] = @benchmarkable ComputationalHomology.position(Cbig, Cbig.cells[0][1])
suite["parts"]["position-big-center"] = @benchmarkable ComputationalHomology.position(Cbig, Cbig.cells[0][length(Cbig.cells[0])>>1])
suite["parts"]["position-big-end"] = @benchmarkable ComputationalHomology.position(Cbig, Cbig.cells[0][length(Cbig.cells[0])])

ɛ = 0.25
Asml = ComputationalHomology.adjacency_matrix(Csml, Bool)
Amed = ComputationalHomology.adjacency_matrix(Cmed, Bool)
Abig = ComputationalHomology.adjacency_matrix(Cbig, Bool)

suite["construction"] = BenchmarkGroup(["vrips", "witness"])
suite["construction"]["lowernbrs-small"] = @benchmarkable ComputationalHomology.lowernbrs(Csml, 10, Asml)
suite["construction"]["lowernbrs-medium"] = @benchmarkable ComputationalHomology.lowernbrs(Cmed, 103, Amed)
suite["construction"]["lowernbrs-big"] = @benchmarkable ComputationalHomology.lowernbrs(Cbig, 1009, Abig)
suite["construction"]["vr-dim1"] = @benchmarkable vietorisrips(Xmed, ɛ, false, maxoutdim=1)
suite["construction"]["vr-incremental"] = @benchmarkable vietorisrips(Xmed, ɛ, false, maxoutdim=2, expansion=:incremental)
suite["construction"]["vr-inductive"] = @benchmarkable vietorisrips(Xmed, ɛ, false, maxoutdim=2, expansion=:inductive)
suite["construction"]["witness-dim1"] = @benchmarkable witness(Xmed, l, ɛ, false, maxoutdim=1)
suite["construction"]["witness-incremental"] = @benchmarkable witness(Xmed, l, ɛ, false, maxoutdim=2, expansion=:incremental)
suite["construction"]["witness-inductive"] = @benchmarkable witness(Xmed, l, ɛ, false, maxoutdim=2, expansion=:inductive)

Fsml = filtration(Csml, wsml)
Fmed = filtration(Cmed, wmed)
Fbig = filtration(Cbig, wbig)

suite["persistence"] = BenchmarkGroup(["homology", "diagram"])
suite["persistence"]["standard-med"] = @benchmarkable diagram(StandardReduction, Fmed)
suite["persistence"]["standard-big"] = @benchmarkable diagram(StandardReduction, Fbig)
suite["persistence"]["twist-med"] = @benchmarkable diagram(TwistReduction, Fmed)
suite["persistence"]["twist-big"] = @benchmarkable diagram(TwistReduction, Fbig)

# parameters
loadparams!(suite, BenchmarkTools.load("bmsetup.json")[1], :evals, :samples);
# tune!(suite);
# BenchmarkTools.save("bmsetup.json", params(suite));

# benchmark
results = run(suite, verbose = true)
BenchmarkTools.save("results.json", results)


text:
using FiniteDiff
using Random
Random.seed!(1)
f(x) = x[1]^2 + cos(x[2] + 2x[6]) + x[5]*x[4] + x[3]^2
H = zeros(6,6)
x = randn(6)
cache = FiniteDiff.HessianCache(copy(x))
FiniteDiff.finite_difference_hessian!(H,f,x,cache)
H[1,1]

x = randn(6) * 1e6
FiniteDiff.finite_difference_hessian!(H,f,x,cache)
H[1,1]

cache.xmm .= x
cache.xmp .= x
cache.xpm .= x
cache.xpp .= x
FiniteDiff.finite_difference_hessian!(H,f,x,cache)
H[1,1]


text:
translate(cursor::CLDeclRefExpr) = MetaExpr(Symbol(spelling(cursor)), cursor)

function translate(cursor::CLInitListExpr)
    child_cursors = children(cursor)
    list = []
    for c in child_cursors
        literal_meta = translate(c)
        push!(list, literal_meta.expr)
    end
    return MetaExpr(Expr(:vect, list...), cursor)
end

function translate(cursor::CLArraySubscriptExpr)
    child_cursors = children(cursor)
    first_meta = translate(first(child_cursors))
    last_meta = translate(last(child_cursors))
    # Julia is 1-based
    idx = last_meta.expr isa Number ? last_meta.expr+1 : Expr(:call, :+, last_meta.expr, 1)
    return MetaExpr(Expr(:ref, first_meta.expr, idx), cursor)
end

function getranges(cursor::CLCursor)
    tu = clang_Cursor_getTranslationUnit(cursor)
    source_range = clang_getCursorExtent(cursor)
    range_start = clang_getRangeStart(source_range)
    range_end = clang_getRangeEnd(source_range)

    file = Ref{CXFile}(C_NULL)
    start_line = Ref{Cuint}(0)
    end_line = Ref{Cuint}(0)
    start_column = Ref{Cuint}(0)
    end_column = Ref{Cuint}(0)
    offset = Ref{Cuint}(0)
    clang_getExpansionLocation(range_start, file, start_line, start_column, offset)
    clang_getExpansionLocation(range_end, file, end_line, end_column, offset)

    expanded_start = clang_getLocation(tu, file[], start_line[], start_column[])
    expanded_end = clang_getLocation(tu, file[], end_line[], end_column[])
    expanded_range = clang_getRange(expanded_start, expanded_end)

    return source_range, expanded_range
end

function translate(cursor::CLParenExpr)
    source, expanded = getranges(cursor)
    if source == expanded
        return translate(first(children(cursor)))
    else
        # we do not translate function-like macros
        tu = clang_Cursor_getTranslationUnit(cursor)
        toks = Clang.TokenList(tu, expanded)
        code = mapreduce(x->x.text, *, collect(toks))
        expr = Expr(:macrocall, Symbol("@warn"), nothing, code)
        return MetaExpr(expr, cursor)
    end
end

function translate(cursor::CLMemberRefExpr)
    sym = Symbol(spelling(cursor))
    meta = translate(first(children(cursor)))
    meta.expr = Expr(:(.), meta.expr, Meta.quot(sym))
    return meta
end

function translate(cursor::CLCStyleCastExpr)
    child_cursors = children(cursor)
    jltype_sym = clang2julia(type(cursor))
    if jltype_sym == :Cvoid
        return translate(last(child_cursors))  # don't cast anything to Nothing
    else
        meta = translate(last(child_cursors))
        meta.expr = Expr(:call, jltype_sym, meta.expr)
        return meta
    end
end

function translate(cursor::CLUnexposedExpr)
    child_cursors = children(cursor)
    return translate(first(child_cursors))
end

function translate(cursor::CLUnaryExpr)
    child_cursors = children(cursor)
    return translate(first(child_cursors))
end

function translate(cursor::CLCallExpr)
    child_cursors = children(cursor)
    func_name = spelling(cursor)
    func_expr = Expr(:call, Symbol(func_name))
    has_address_op = false
    for c in child_cursors[2:end]
        meta = translate(c)
        if meta.info == "AddressOperator"
            has_address_op = true
            push!(func_expr.args, Expr(:&, meta.expr.args[2]))
        else
            push!(func_expr.args, meta.expr)
        end
    end
    if has_address_op
        return MetaExpr(Expr(:macrocall, Symbol("@c"), nothing, func_expr))
    else
        return MetaExpr(func_expr)
    end
end


text:
## this was my original day13 main, unfortunately overwrote splitX and Y

function day13()
  lines = readlines("inputs/11-15/13.txt")
  nums = map(str -> [parse(Int,str[l]) for l in findall(r"[[:digit:]]+", str)], lines)
  S = Set(nums)

  S = splitX(S, 655)
  S = splitY(S, 447)
  S = splitX(S, 327)
  S = splitY(S, 223)
  S = splitX(S, 163)
  S = splitY(S, 111)
  S = splitX(S, 81)
  S = splitY(S, 55)
  S = splitX(S, 40)
  S = splitY(S, 27)
  S = splitX(S, 13)
  S = splitY(S, 6)
  @show length(S)
end


text:
int main() {
    while(false) return 0;
}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
#==============================================================================#
# AWSCore/test/runtests.jl
#
# Copyright OC Technology Pty Ltd 2014 - All rights reserved
#==============================================================================#

using Base.Test
using AWSCore
using SymDict
using Retry
using XMLDict

using AWSCore: service_query

AWSCore.set_debug_level(1)

@testset "AWSCore" begin

aws = aws_config()

@testset "Load Credentials" begin
    user = aws_user_arn(aws)

    @test ismatch(r"^arn:aws:iam::[0-9]+:[^:]+$", user)

    println("Authenticated as: $user")

    aws[:region] = "us-east-1"

    println("Testing exceptions...")
    try
        AWSCore.Services.iam("GetFoo", Dict("ContentType" => "JSON"))
        @test false
    catch e
        println(e)
        @test ecode(e) == "InvalidAction"
    end

    try
        AWSCore.Services.iam("GetUser", Dict("UserName" => "notauser",
                                           "ContentType" => "JSON"))
        @test false
    catch e
        println(e)
        @test ecode(e) in ["AccessDenied", "NoSuchEntity"]
    end

    try
        AWSCore.Services.iam("GetUser", Dict("UserName" => "@#!%%!",
                                           "ContentType" => "JSON"))
        @test false
    catch e
        println(e)
        @test ecode(e) == "ValidationError"
    end

    try
        AWSCore.Services.iam("CreateUser", Dict("UserName" => "root",
                                              "ContentType" => "JSON"))
        @test false
    catch e
        println(e)
        @test ecode(e) in ["AccessDenied", "EntityAlreadyExists"]
    end
end

@testset "AssumeRole" begin
    mktemp() do config_file, config_io
        write(config_io, """[profile test]
                output = json
                region = us-east-1

                [profile test:dev]
                source_profile = test
                role_arn = arn:aws:iam::123456789000:role/Dev

                [profile test:sub-dev]
                source_profile = test:dev
                role_arn = arn:aws:iam::123456789000:role/SubDev

                [profile test2]
                aws_access_key_id = WRONG_ACCESS_ID
                aws_secret_access_key = WRONG_ACCESS_KEY
                output = json
                region = us-east-1

                [profile test3]
                source_profile = test:dev
                role_arn = arn:aws:iam::123456789000:role/test3

                [profile test4]
                source_profile = test:dev
                role_arn = arn:aws:iam::123456789000:role/test3
                aws_access_key_id = RIGHT_ACCESS_ID4
                aws_secret_access_key = RIGHT_ACCESS_KEY4
            """)
        close(config_io)

        mktemp() do creds_file, creds_io
            write(creds_io, """[test]
                aws_access_key_id = TEST_ACCESS_ID
                aws_secret_access_key = TEST_ACCESS_KEY

                [test2]
                aws_access_key_id = RIGHT_ACCESS_ID2
                aws_secret_access_key = RIGHT_ACCESS_KEY2

                [test3]
                aws_access_key_id = RIGHT_ACCESS_ID3
                aws_secret_access_key = RIGHT_ACCESS_KEY3
                """)
            close(creds_io)

            withenv(
                "AWS_SHARED_CREDENTIALS_FILE" => creds_file,
                "AWS_CONFIG_FILE" => config_file,
                "AWS_DEFAULT_PROFILE" => "test"
                ) do

                # Check credentials load
                config = AWSCore.aws_config()
                creds = config[:creds]

                @test creds.access_key_id == "TEST_ACCESS_ID"
                @test creds.secret_key == "TEST_ACCESS_KEY"

                # Check credential file takes precedence over config
                ENV["AWS_DEFAULT_PROFILE"] = "test2"
                config = AWSCore.aws_config()
                creds = config[:creds]

                @test creds.access_key_id == "RIGHT_ACCESS_ID2"
                @test creds.secret_key == "RIGHT_ACCESS_KEY2"

                # Check credentials take precedence over role
                ENV["AWS_DEFAULT_PROFILE"] = "test3"
                config = AWSCore.aws_config()
                creds = config[:creds]

                @test creds.access_key_id == "RIGHT_ACCESS_ID3"
                @test creds.secret_key == "RIGHT_ACCESS_KEY3"

                ENV["AWS_DEFAULT_PROFILE"] = "test4"
                config = AWSCore.aws_config()
                creds = config[:creds]

                @test creds.access_key_id == "RIGHT_ACCESS_ID4"
                @test creds.secret_key == "RIGHT_ACCESS_KEY4"

                # Check we try to assume a role
                ENV["AWS_DEFAULT_PROFILE"] = "test:dev"

                try
                    AWSCore.aws_config()
                    @test false
                catch e
                    @test e isa AWSCore.AWSException
                    @test ecode(e) == "InvalidClientTokenId"
                end

                # Check we try to assume a role
                ENV["AWS_DEFAULT_PROFILE"] = "test:sub-dev"
                let oldout = STDOUT
                    r,w = redirect_stdout()
                    try
                        AWSCore.aws_config()
                        @test false
                    catch e
                        @test e isa AWSCore.AWSException
                        @test ecode(e) == "InvalidClientTokenId"
                    end
                    redirect_stdout(oldout)
                    close(w)
                    output = convert(String, read(r))
                    contains(output, "Assuming \"test:dev\"")
                    contains(output, "Assuming \"test\"")
                    close(r)
                end
            end
        end
    end
end

@testset "XML Parsing" begin
    XML(x)=parse_xml(x)

    xml = """
    <CreateQueueResponse>
        <CreateQueueResult>
            <QueueUrl>
                http://queue.amazonaws.com/123456789012/testQueue
            </QueueUrl>
        </CreateQueueResult>
        <ResponseMetadata>
            <RequestId>
                7a62c49f-347e-4fc4-9331-6e8e7a96aa73
            </RequestId>
        </ResponseMetadata>
    </CreateQueueResponse>
    """

    @assert XML(xml)["CreateQueueResult"]["QueueUrl"] ==
          "http://queue.amazonaws.com/123456789012/testQueue"

    xml = """
    <GetUserResponse xmlns="https://iam.amazonaws.com/doc/2010-05-08/">
      <GetUserResult>
        <User>
          <PasswordLastUsed>2015-12-23T22:45:36Z</PasswordLastUsed>
          <Arn>arn:aws:iam::012541411202:root</Arn>
          <UserId>012541411202</UserId>
          <CreateDate>2015-09-15T01:07:23Z</CreateDate>
        </User>
      </GetUserResult>
      <ResponseMetadata>
        <RequestId>837446c9-abaf-11e5-9f63-65ae4344bd73</RequestId>
      </ResponseMetadata>
    </GetUserResponse>
    """

    @test XML(xml)["GetUserResult"]["User"]["Arn"] == "arn:aws:iam::012541411202:root"

    xml = """
    <GetQueueAttributesResponse>
      <GetQueueAttributesResult>
        <Attribute>
          <Name>ReceiveMessageWaitTimeSeconds</Name>
          <Value>2</Value>
        </Attribute>
        <Attribute>
          <Name>VisibilityTimeout</Name>
          <Value>30</Value>
        </Attribute>
        <Attribute>
          <Name>ApproximateNumberOfMessages</Name>
          <Value>0</Value>
        </Attribute>
        <Attribute>
          <Name>ApproximateNumberOfMessagesNotVisible</Name>
          <Value>0</Value>
        </Attribute>
        <Attribute>
          <Name>CreatedTimestamp</Name>
          <Value>1286771522</Value>
        </Attribute>
        <Attribute>
          <Name>LastModifiedTimestamp</Name>
          <Value>1286771522</Value>
        </Attribute>
        <Attribute>
          <Name>QueueArn</Name>
          <Value>arn:aws:sqs:us-east-1:123456789012:qfoo</Value>
        </Attribute>
        <Attribute>
          <Name>MaximumMessageSize</Name>
          <Value>8192</Value>
        </Attribute>
        <Attribute>
          <Name>MessageRetentionPeriod</Name>
          <Value>345600</Value>
        </Attribute>
      </GetQueueAttributesResult>
      <ResponseMetadata>
        <RequestId>1ea71be5-b5a2-4f9d-b85a-945d8d08cd0b</RequestId>
      </ResponseMetadata>
    </GetQueueAttributesResponse>
    """

    d = Dict(a["Name"] => a["Value"] for a in XML(xml)["GetQueueAttributesResult"]["Attribute"])

    @test d["MessageRetentionPeriod"] == "345600"
    @test d["CreatedTimestamp"] == "1286771522"

    xml = """
    <?xml version="1.0" encoding="UTF-8"?>
    <ListAllMyBucketsResult xmlns="http://s3.amazonaws.com/doc/2006-03-01">
      <Owner>
        <ID>bcaf1ffd86f461ca5fb16fd081034f</ID>
        <DisplayName>webfile</DisplayName>
      </Owner>
      <Buckets>
        <Bucket>
          <Name>quotes</Name>
          <CreationDate>2006-02-03T16:45:09.000Z</CreationDate>
        </Bucket>
        <Bucket>
          <Name>samples</Name>
          <CreationDate>2006-02-03T16:41:58.000Z</CreationDate>
        </Bucket>
      </Buckets>
    </ListAllMyBucketsResult>
    """

    @test map(b->b["Name"], XML(xml)["Buckets"]["Bucket"]) == ["quotes", "samples"]

    xml = """
    <ListDomainsResponse>
      <ListDomainsResult>
        <DomainName>Domain1</DomainName>
        <DomainName>Domain2</DomainName>
        <NextToken>TWV0ZXJpbmdUZXN0RG9tYWluMS0yMDA3MDYwMTE2NTY=</NextToken>
      </ListDomainsResult>
      <ResponseMetadata>
        <RequestId>eb13162f-1b95-4511-8b12-489b86acfd28</RequestId>
        <BoxUsage>0.0000219907</BoxUsage>
      </ResponseMetadata>
    </ListDomainsResponse>
    """

    @test XML(xml)["ListDomainsResult"]["DomainName"] == ["Domain1", "Domain2"]
end

@testset "AWS Signature Version 4" begin
    function aws4_request_headers_test()

        r = @SymDict(
            creds         = AWSCredentials(
                                "AKIDEXAMPLE",
                                "wJalrXUtnFEMI/K7MDENG+bPxRfiCYEXAMPLEKEY"
                            ),
            region        = "us-east-1",
            verb          = "POST",
            service       = "iam",
            url           = "http://iam.amazonaws.com/",
            content       = "Action=ListUsers&Version=2010-05-08",
            headers       = Dict(
                                "Content-Type" =>
                                "application/x-www-form-urlencoded; charset=utf-8",
                                "Host" => "iam.amazonaws.com"
                            )
        )

        AWSCore.sign!(r, DateTime("2011-09-09T23:36:00"))

        h = r[:headers]
        out = join(["$k: $(h[k])\n" for k in sort(collect(keys(h)))])

        expected = (
            "Authorization: AWS4-HMAC-SHA256 " *
            "Credential=AKIDEXAMPLE/20110909/us-east-1/iam/aws4_request, " *
            "SignedHeaders=content-md5;content-type;host;" *
            "x-amz-content-sha256;x-amz-date, " *
            "Signature=1a6db936024345449ef4507f890c5161" *
                     "bbfa2ff2490866653bb8b58b7ba1554a\n" *
            "Content-MD5: r2d9jRneykOuUqFWSFXKCg==\n" *
            "Content-Type: application/x-www-form-urlencoded; " *
                         "charset=utf-8\n" *
            "Host: iam.amazonaws.com\n" *
            "x-amz-content-sha256: b6359072c78d70ebee1e81adcbab4f01" *
                                 "bf2c23245fa365ef83fe8f1f955085e2\n" *
            "x-amz-date: 20110909T233600Z\n")

        @test out == expected
    end

    aws4_request_headers_test()
end

@testset "ARN" begin
    @test arn(aws,"s3","foo/bar") == "arn:aws:s3:::foo/bar"
    @test arn(aws,"s3","foo")     == "arn:aws:s3:::foo"
    @test arn(aws,"sqs", "au-test-queue", "ap-southeast-2", "1234") ==
          "arn:aws:sqs:ap-southeast-2:1234:au-test-queue"

    @test arn(aws,"sns","*","*",1234) == "arn:aws:sns:*:1234:*"
    @test arn(aws,"iam","role/foo-role", "", 1234) ==
          "arn:aws:iam::1234:role/foo-role"
end

@testset "Misc" begin
    @test HTTP.escapepath("invocations/function:f:PROD") ==
                          "invocations/function%3Af%3APROD"
end

end # testset "AWSCore"

#==============================================================================#
# End of file.
#==============================================================================#


text:
using Twilert
using Base.Test
using JSON

cfgfile = open(dirname(@__FILE__) * "/../.twilertrc", "r")
cfg = JSON.parse(cfgfile)
close(cfgfile)

const sid = cfg["sid"]
const token = cfg["token"]
const to = cfg["to"]
const from = cfg["from"]

const to_invalid = "+15005550001"
const from_invalid = "+15005550001"

@test Twilert.alert("Test message") == :ok
@test Twilert.alert("Test message", to) == :ok
@test Twilert.alert("Test message", to, from) == :ok

@test Twilert.alert("Test message", to_invalid) == :error
@test Twilert.alert("Test message", to, from_invalid) == :error


text:
using Controlz

# https://lpsa.swarthmore.edu/Root_Locus/RLocusExamples.html
g = (s + 1) / (s^3 + 4 * s^2 + 6 * s + 4)
g = 1 / (s * (s + 3))
g = (s + 3) / (s^2 - s - 2)
g = (s^2 + 2 * s + 2) / (s * (s^4 + 9 * s^3 + 33 * s^2 + 51 * s + 26))

root_locus(g, max_mag_Kc=3.0, nb_pts=5000)

max_stars_repo_name: SimonEnsemble/Controlz

max_issues_repo_name: SimonEnsemble/Controlz

max_forks_repo_name: SimonEnsemble/Controlz


text:
module Merly
import Base.|

using Sockets
using JSON
using HTTP
#XMLDict

include("base.jl")
include("mimetypes.jl")
include("routes.jl")
include("allformats.jl")

export App, @page, @route, GET,POST,PUT,DELETE,HEAD,OPTIONS,PATCH,Get,Post,Put,Delete

end # module


text:
using SpecialFunctions
"Abstract quadrature type"
abstract type QUADRATURE_TYPE end

"Gauss-Jacobi quadrature type"
struct GJ <: QUADRATURE_TYPE end
"Gauss-Lobatto-Jacobi quadrature type"
struct GLJ <: QUADRATURE_TYPE end
"Gauss-Radau-Jacobi quadrature type includind the left end"
struct GRJM <: QUADRATURE_TYPE end
"Gauss-Radau-Jacobi quadrature type includind the right end"
struct GRJP <: QUADRATURE_TYPE end

# Calculate a ratio of Gamma functions without overflow
function gamma_ratio(num::T, denom::T) where {T<:Number}
    if num>1 && denom>1
        exp((logabsgamma(num))[1] - (logabsgamma(denom))[1])
    else
        gamma(num) / gamma(denom)
    end
end

"""
Gauss-type quadrature

Numerical integrals in the domain [-1,1] can be computed
from knowledge of the function in a set of nodes and the
corresponding nodes, such that

```math
\\int_{-1}^1 (1-x)^a (1+x)^b (1-xf(x)\\:dx \\approx \\sum_{i=1}^N w^{a,b}_i f(x_i)
```

The parameters `a` and `b` form a famility of quadrature rules. But if one or both
of the ends are specified, other quadrature families are possible:

 * No ends are specified, Gauss-Jacobi quadrature (GJ)
 * Both ends are specified, Gauss-Lobatto-Jacobi quadrature (GLJ)
 * A single end is specified, Gauss-Radau-Jacobi quadrature (if +1 is specified, GRJP, if -1, GRJM)

To compute the nodes, the following functions are available:

 * `zgj` (Gauss-Jacobi)
 * `zglj` (Gauss-Lobatto-Jacobi)
 * `zgrjm` (Gauss-Radau-Jacobi, -1)
 * `zgrjp` (Gauss-Radau-Jacobi, +1)

All these functions have the following parameters:

 * `Q`  Number of quadrature nodes
 * `a` (a) weight
 * `b` (b) weight
 *  type Data type to be used, `Float64` is the default

To compute the weights, first the zeros (`z`) should be computed and then the weights are computed
with the following functions:

 * `wgj(z, a, b)`
 * `wglj(z, a, b)`
 * `wgrjm(z, a, b)`
 * `wgrjp(z, a, b)`

 ### Derivatives

The nodes used in the quadrature rules are convenient when using high order Lagrange interpolation
To compute derivatives. The following functions are used to compute the derivative matrix such that
`du = D*u`:

 * `dgj(z, a, b)`
 * `dglj(z, a, b)`
 * `dgrjm(z, a, b)`
 * `dgrjp(z, a, b)`

### Examples

See the notebooks availbale with the package.

"""
function zgj(Q, a, b, ::Type{T}=Float64) where {T<:Number}
    jacobi_zeros(Q, a, b, T)
end

zgj(Q) = zgj(Q, 0.0, 0.0)
zgj(Q, a) = zgj(Q, a, zero(a))

function zglj(Q, a, b, ::Type{T}=Float64) where {T<:Number}
    z = jacobi_zeros(Q-2, a+1, b+1, T)
    o = one(T)
    return [-o; z; o]
end
zglj(Q) = zglj(Q, 0.0, 0.0)
zglj(Q, a) = zglj(Q, a, zero(a))

function zgrjm(Q, a, b, ::Type{T}=Float64) where {T<:Number}
    z = jacobi_zeros(Q-1, a, b+1, T)
    return [-one(T); z]
end
zgrjm(Q) = zgrjm(Q, 0.0, 0.0)
zgrjm(Q, a) = zgrjm(Q, a, zero(a))

function zgrjp(Q, a, b, ::Type{T}=Float64) where {T<:Number}
    z = jacobi_zeros(Q-1, a+1, b, T)
    return [z; one(T)]
end
zgrjp(Q) = zgrjp(Q, 0.0, 0.0)
zgrjp(Q, a) = zgrjp(Q, a, zero(a))

function wgj(z::AbstractArray{T}, alpha=0, beta=0) where {T<:Number}
    a = convert(T, alpha)
    b = convert(T, beta)

    Q::Int = length(z)
    o = one(T)
    coef = 2^(a+b+1) * gamma_ratio(a+Q+1, Q+o) * gamma_ratio(b+Q+1, a+b+Q+1)
    w = [djacobi(zz, Q, a, b) for zz=z]

    for i = 1:Q
        ww = w[i]
        x = z[i]
        w[i] = o / (ww*ww) * coef / (o - x*x)
    end

    return w
end

function wglj(z::AbstractArray{T}, alpha=0, beta=0) where {T<:Number}
    a = convert(T, alpha)
    b = convert(T, beta)
    o = one(T)
    Q = length(z)

    coef = 2^(a+b+1) / (Q-o) * gamma_ratio(a+Q, Q*o) * gamma_ratio(b+Q, a+b+Q+1)

    w = [jacobi(zz, Q-1, a, b) for zz=z]
    w[1] = (b+1) * coef / (w[1]*w[1])
    w[Q] = (a+1) * coef / (w[Q]*w[Q])

    for i = 2:(Q-1)
        ww = w[i]
        w[i] = coef / (ww * ww)
    end

    return w
end

function wgrjm(z::AbstractArray{T}, alpha=0, beta=0) where {T<:Number}
    a = convert(T, alpha)
    b = convert(T, beta)
    o = one(T)

    Q = length(z)

    coef = 2^(a+b) / (b+Q) * gamma_ratio(a+Q, Q*o) * gamma_ratio(b+Q, a+b+Q+1)

    w = [jacobi(zz, Q-1, a, b) for zz=z]

    for i = 1:Q
        ww = w[i]
        w[i] = coef / (ww*ww) * (o - z[i])
    end

    w[1] *= (b + o)

    return w
end

function wgrjp(z::AbstractArray{T,1}, alpha=0, beta=0) where {T<:Number}
    a = convert(T, alpha)
    b = convert(T, beta)
    Q = length(z)
    o = one(T)

    coef = 2^(a+b) / (a+Q) * gamma_ratio(a+Q, Q*o) * gamma_ratio(b+Q, a+b+Q+1)

    w = [jacobi(zz, Q-1, a, b) for zz=z]

    for i = 1:Q
        ww = w[i]
        w[i] = coef / (ww*ww) * (o + z[i])
    end

    w[Q] *= (a + o)

    return w
end

function dgj(z::AbstractArray{T,1}, alpha=0, beta=0) where {T<:Number}

    Q = length(z)
    a = convert(T, alpha)
    b = convert(T, beta)
    o = one(T)

    djac = [djacobi(zz, Q, alpha, beta) for zz=z]

    D = zeros(T, Q, Q)
    for i = 1:Q
        for k = 1:Q
            if i != k
                D[i,k] = (djac[i]/djac[k]) / (z[i]-z[k])
            else
                D[i,i] = (a-b + (a + b + 2) * z[i]) / (o - z[k]^2) / 2
            end
        end
    end

    return D
end

function dglj(z::AbstractArray{T,1}, alpha=0, beta=0) where {T<:Number}

    Q = length(z)
    a = convert(T, alpha)
    b = convert(T, beta)
    o = one(T)

    djac = zeros(T,Q)
    djac[1] = (-1)^Q * 2 * gamma_ratio(Q+b, Q-o) / gamma(b+2)
    djac[Q] = -2 * gamma_ratio(Q+a, Q-o) / gamma(a+2)
    for i = 2:(Q-1)
        djac[i] = (o-z[i]*z[i]) * djacobi(z[i], Q-2, a+1, b+1)
    end

    D = zeros(T, Q, Q)
    for i = 1:Q
        for k = 1:Q
            if i != k
                D[i,k] = (djac[i]/djac[k]) / (z[i]-z[k])
            else
                D[i,i] = (a - b + (a + b)*z[i]) / (2*(o - z[i]^2))
            end
        end
    end

    D[1,1] = (a - (Q-1)*(Q+a+b)) / (2*(b + 2))
    D[Q,Q] = -(b - (Q-1)*(Q+a+b)) / (2*(a + 2))

    return D

end

function dgrjm(z::AbstractArray{T,1}, alpha=0, beta=0) where {T<:Number}

    Q = length(z)
    a = convert(T, alpha)
    b = convert(T, beta)
    o = one(T)

    djac = zeros(T,Q)
    for i = 2:Q
        djac[i] = (1+z[i]) * djacobi(z[i], Q-1, a, b+1)
    end
    djac[1] = (-1)^(Q-1) * gamma_ratio(Q+b+1, Q*o) / gamma(b+2)

    D = zeros(T, Q, Q)
    for i = 1:Q
        for k = 1:Q
            if i != k
                D[i,k] = (djac[i]/djac[k]) / (z[i]-z[k])
            else
                D[i,i] = (a - b + 1 + (a + b + 1)*z[i]) / (2*(1 - z[i]^2))
            end
        end
    end

    D[1,1] = -(Q-1)*(Q+a+b+1) / (2*(b + 2))

    return D

end

function dgrjp(z::AbstractArray{T,1}, alpha=0, beta=0) where {T<:Number}

    Q = length(z)
    a = convert(T, alpha)
    b = convert(T, beta)
    o = one(T)

    djac = zeros(T,Q)
    for i = 1:(Q-1)
        djac[i] = (1-z[i]) * djacobi(z[i], Q-1, a+1, b)
    end
    djac[Q] = - gamma_ratio(Q+a+1, Q*o) / gamma(a+2)

    D = zeros(T, Q, Q)
    for i = 1:Q
        for k = 1:Q
            if i != k
                D[i,k] = (djac[i]/djac[k]) / (z[i]-z[k])
            else
                D[i,i] = (a - b - 1 + (a + b + 1)*z[i]) / (2*(1 - z[i]^2))
            end
        end
    end

    D[Q,Q] = (Q-1)*(Q+a+b+1) / (2*(a + 2))

    return D

end

"""
Abstract interface of Gauss-type quadrature rules

Can be used for any `AbstractFloat` type data.
"""
struct Quadrature{T<:Number,QT<:QUADRATURE_TYPE}
    "Number of quadrature nodes"
    Q::Int
    "a weight"
    a::T
    "b weight"
    b::T
    "Quadrature nodes"
    z::Array{T, 1}
    "Quadrature weights"
    w::Array{T, 1}
    "Quadrature derivative matrix"
    D::Array{T, 2}
end

"""
Return the zeros of a Gauss type quadrature
"""
qzeros(::Type{GJ}, Q, a=0, b=0, ::Type{T}=Float64) where {T<:Number} = zgj(Q, a, b, T)
qzeros(::Type{GLJ}, Q, a=0, b=0, ::Type{T}=Float64) where {T<:Number} = zglj(Q, a, b, T)
qzeros(::Type{GRJM}, Q, a=0, b=0, ::Type{T}=Float64) where {T<:Number} = zgrjm(Q, a, b, T)
qzeros(::Type{GRJP}, Q, a=0, b=0, ::Type{T}=Float64) where {T<:Number} = zgrjp(Q, a, b, T)

"""
Return the weights of a Gauss type quadrature
"""
qweights(::Type{GJ}, z::AbstractArray{T}, a=0, b=0) where {T<:Number} = wgj(z, a, b)
qweights(::Type{GLJ}, z::AbstractArray{T}, a=0, b=0) where {T<:Number} = wglj(z, a, b)
qweights(::Type{GRJM}, z::AbstractArray{T}, a=0, b=0) where {T<:Number} = wgrjm(z, a, b)
qweights(::Type{GRJP}, z::AbstractArray{T}, a=0, b=0) where {T<:Number} = wgrjp(z, a, b)

"""
Return the derivative matrix of a Gauss type quadrature
"""
qdiff(::Type{GJ}, z::AbstractArray{T}, a=0, b=0) where {T<:Number} = dgj(z, a, b)
qdiff(::Type{GLJ}, z::AbstractArray{T}, a=0, b=0) where {T<:Number} = dglj(z, a, b)
qdiff(::Type{GRJM}, z::AbstractArray{T}, a=0, b=0) where {T<:Number} = dgrjm(z, a, b)
qdiff(::Type{GRJP}, z::AbstractArray{T}, a=0, b=0) where {T<:Number} = dgrjp(z, a, b)

"""
Create a `Quadrature` object given its type, order and weights.
"""
function Quadrature(::Type{QT}, Q, a=0, b=0,
                    ::Type{T}=Float64) where {T<:Number, QT<:QUADRATURE_TYPE}
    aa = convert(T, a)
    bb = convert(T, b)
    z = qzeros(QT, Q, aa, bb, T)
    w = qweights(QT, z, aa, bb)
    D = qdiff(QT, z, aa, bb)
    Quadrature{T,QT}(Q, aa, bb, z, w, D)
end

"Return quadrature type"
qtype(q::Quadrature{T,QT}) where {T,QT} = QT
"Return quadrature nodes"
qzeros(q::QT) where {QT<:Quadrature} = q.z
"Return quadrature weights"
qweights(q::QT) where {QT<:Quadrature} = q.w
"Return quadrature derivative matrix"
qdiff(q::QT) where {QT<:Quadrature} = q.D
"Return number of quadrature nodes"
num_points(q::QT) where {QT<:Quadrature} = q.Q
"Return quadrature `a` weight"
qalpha(q::QT) where {QT<:Quadrature} = q.a
"Return quadrature `b` weight"
qbeta(q::QT) where {QT<:Quadrature} = q.b

"""
Compute the Lagrange polynomial

This function computes the Lagrange polynomial at point `x` corresponding to
the `i`-th node of the set `z`

There is also a modifying version `lagrange!` used to computing the polynomials
at several points of an array `x`
"""
function lagrange(i, x, z)
    nz = length(z)

    l = one(z[1])

    for k = 1:(i-1)
        l = l * (x-z[k]) / (z[i]-z[k])
    end

    for k = (i+1):nz
        l = l * (x-z[k]) / (z[i]-z[k])
    end

    return l
end

"""
Interpolation matrix

Often it is necessary to compute the values of a function approximated using
Lagrange interpolation through a set of points `z`. If this will be repeated
often, a matrix can be computed that allows the easy computation using the simple
expression `fx = Imat * fz`
"""
function interp_mat(x::AbstractArray{T}, z::AbstractArray{T}) where {T<:Number}

    Q = length(z)
    np = length(x)

    Imat = zeros(T, np, Q)
    for i = 1:Q, k=1:np
        Imat[k,i] = lagrange(i, x[k], z)
    end

    return(Imat)

end

@doc (@doc zgj) zglj
@doc (@doc zgj) zgrjm
@doc (@doc zgj) zgrjp

@doc (@doc zgj) wgj
@doc (@doc zgj) wglj
@doc (@doc zgj) wgrjm
@doc (@doc zgj) wgrjp

@doc (@doc zgj) dgj
@doc (@doc zgj) dglj
@doc (@doc zgj) dgrjm
@doc (@doc zgj) dgrjp


text:
##High-level syntax
fn1 = tempname2()
fn2 = tempname2()
fn3 = tempname2()
fn4 = tempname2()
fn5 = tempname2()

nccreate(
    fn1,
    "v1",
    "Dim1",
    [1, 2],
    Dict("units" => "deg C"),
    "Dim2",
    collect(1:10),
    "Dim3",
    20,
    Dict("max" => 10),
    mode = NC_NETCDF4,
)
nccreate(fn1, "vstr", "Dim2", collect(1:10), t = String)
nccreate(fn1, "vchar", "DimChar", 20, "Dim2", t = NetCDF.NC_CHAR)
nccreate(
    fn2,
    "v2",
    "Dim1",
    [1, 2, 3],
    Dict("units" => "deg C"),
    "Dim2",
    collect(1:10),
    "Dim3",
    20,
    Dict("max" => 10),
    atts = Dict("a1" => "varatts"),
    gatts = Dict("Some global attributes" => 2010),
)
nccreate(fn3, "v3", "Dim1", 3)
for i = 1:length(tlist)
    nccreate(fn3, "vt$i", "Dim2", collect(1:10), t = tlist[i])
end

ncputatt(fn1, "v1", Dict("Additional String attribute" => "att"))
ncputatt(fn1, "global", Dict("Additional global attribute" => "gatt"))
numberatts = Dict{Any,Any}("Additional $t attribute" => t(20) for t in tlist)
numberatts["Additional String attribute"] = "string attribute"
arrayatts = Dict{Any,Any}("Additional $t attribute" => t[i for i = 1:20] for t in tlist)
arrayatts["Additional String array attribute"] = String[string("string attribute ", i) for i = 1:20]
ncputatt(fn1, "global", numberatts)
ncputatt(fn1, "v1", arrayatts)

#First generate the data
x1 = rand(2, 10, 20)
x2 = rand(2, 10, 20)
xt = [rand(tl, 10) for tl in tlist]
xs = [
    "a",
    "bb",
    "ccc",
    "dddd",
    "eeeee",
    "ffffff",
    "ggggggg",
    "hhhhhhhh",
    "iiiiiiiii",
    "jjjjjjjjjj",
]
#
# And write it
#
ncwrite(x1, fn1, "v1")
#Test sequential writing along one dimension
for i = 1:10
    ncwrite(x2[:, i, :], fn2, "v2", start = [1, i, 1], count = [-1, 1, -1])
end
ncwrite(xs, fn1, "vstr")
ncwrite(nc_string2char(xs), fn1, "vchar")

#Test automatic type conversion
x4 = [1, 2, 3]
ncwrite(x4, fn3, "v3")

for i = 1:length(tlist)
    ncwrite(xt[i], fn3, "vt$i")
end

# Test creating unlimited dimension file
nccreate(fn4, "myvar", "time", Inf)
nccreate(fn4, "myvar2", "dinf", Inf, "d1", 10, "d2", 20)
ncwrite(collect(1:10), fn4, "myvar")
ncwrite(collect(1:5), fn4, "myvar", start = [11;])
ncwrite(ones(5, 10, 20), fn4, "myvar2")
ncwrite(zeros(1, 10, 20), fn4, "myvar2")

# Test reading existing files
myvar2 = ones(5, 10, 20)
myvar2[1, :, :] .= 0.0
@test ncread(fn4, "myvar2") == myvar2
@test nc_char2string(ncread(fn1, "vchar")) == xs
@test_throws NetCDF.NetCDFError ncread("nonexistant file", "a var")

# Test creating 0-dimensional variable
nccreate(fn4, "scalar")
a = Array{Float64,0}(undef);
a[1] = 10.0;
ncwrite(a, fn4, "scalar")
@test ncread(fn4, "scalar")[1] == a[1]

# Test dimension variable type
nccreate(fn5, "x5_1", "dim1", collect(Int32, 1:4), Dict("units" => "m"))
nccreate(fn5, "x5_2", "dim2", collect(Float32, 1:4), Dict("units" => "m"))
@test typeof(ncread(fn5, "dim1")) == Array{Int32,1}
@test typeof(ncread(fn5, "dim2")) == Array{Float32,1}

# Open file by reading and create new variable
ncread(fn4, "myvar2")
nccreate(fn4, "myvar3", "time", Inf)

# Try opening a file that doesn't exist
@test_throws ArgumentError ncread("", "")
@test_throws ArgumentError nccreate("", "")
@test_throws NetCDF.NetCDFError ncread("DoesNotExist.nc", "")


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

mutable struct ServiceEndpointPolicyDefinitionPropertiesFormat <: SwaggerModel
    description::Any # spec type: Union{ Nothing, String } # spec name: description
    service::Any # spec type: Union{ Nothing, String } # spec name: service
    serviceResources::Any # spec type: Union{ Nothing, Vector{String} } # spec name: serviceResources
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState

    function ServiceEndpointPolicyDefinitionPropertiesFormat(;description=nothing, service=nothing, serviceResources=nothing, provisioningState=nothing)
        o = new()
        validate_property(ServiceEndpointPolicyDefinitionPropertiesFormat, Symbol("description"), description)
        setfield!(o, Symbol("description"), description)
        validate_property(ServiceEndpointPolicyDefinitionPropertiesFormat, Symbol("service"), service)
        setfield!(o, Symbol("service"), service)
        validate_property(ServiceEndpointPolicyDefinitionPropertiesFormat, Symbol("serviceResources"), serviceResources)
        setfield!(o, Symbol("serviceResources"), serviceResources)
        validate_property(ServiceEndpointPolicyDefinitionPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type ServiceEndpointPolicyDefinitionPropertiesFormat

const _property_map_ServiceEndpointPolicyDefinitionPropertiesFormat = Dict{Symbol,Symbol}(Symbol("description")=>Symbol("description"), Symbol("service")=>Symbol("service"), Symbol("serviceResources")=>Symbol("serviceResources"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_ServiceEndpointPolicyDefinitionPropertiesFormat = Dict{Symbol,String}(Symbol("description")=>"String", Symbol("service")=>"String", Symbol("serviceResources")=>"Vector{String}", Symbol("provisioningState")=>"ProvisioningState")
Base.propertynames(::Type{ ServiceEndpointPolicyDefinitionPropertiesFormat }) = collect(keys(_property_map_ServiceEndpointPolicyDefinitionPropertiesFormat))
Swagger.property_type(::Type{ ServiceEndpointPolicyDefinitionPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ServiceEndpointPolicyDefinitionPropertiesFormat[name]))}
Swagger.field_name(::Type{ ServiceEndpointPolicyDefinitionPropertiesFormat }, property_name::Symbol) =  _property_map_ServiceEndpointPolicyDefinitionPropertiesFormat[property_name]

function check_required(o::ServiceEndpointPolicyDefinitionPropertiesFormat)
    true
end

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


text:
using CImGui
using CImGui.CSyntax
using CImGui.GLFWBackend
using CImGui.OpenGLBackend
using CImGui.GLFWBackend.GLFW
using CImGui.OpenGLBackend.ModernGL
using Printf

include(joinpath(@__DIR__, "demo_window.jl"))

@static if Sys.isapple()
    # OpenGL 3.2 + GLSL 150
    const glsl_version = 150
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 2)
    GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
    GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # required on Mac
else
    # OpenGL 3.0 + GLSL 130
    const glsl_version = 130
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 0)
    # GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
    # GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # 3.0+ only
end

# setup GLFW error callback
error_callback(err::GLFW.GLFWError) = @error "GLFW ERROR: code $(err.code) msg: $(err.description)"
GLFW.SetErrorCallback(error_callback)

# create window
window = GLFW.CreateWindow(1280, 720, "Demo")
@assert window != C_NULL
GLFW.MakeContextCurrent(window)
GLFW.SwapInterval(1)  # enable vsync

# setup Dear ImGui context
ctx = CImGui.CreateContext()

# setup Dear ImGui style
CImGui.StyleColorsDark()
# CImGui.StyleColorsClassic()
# CImGui.StyleColorsLight()

# load Fonts
# - If no fonts are loaded, dear imgui will use the default font. You can also load multiple fonts and use `CImGui.PushFont/PopFont` to select them.
# - `CImGui.AddFontFromFileTTF` will return the `Ptr{ImFont}` so you can store it if you need to select the font among multiple.
# - If the file cannot be loaded, the function will return C_NULL. Please handle those errors in your application (e.g. use an assertion, or display an error and quit).
# - The fonts will be rasterized at a given size (w/ oversampling) and stored into a texture when calling `CImGui.Build()`/`GetTexDataAsXXXX()``, which `ImGui_ImplXXXX_NewFrame` below will call.
# - Read 'fonts/README.txt' for more instructions and details.
fonts_dir = joinpath(@__DIR__, "..", "fonts")
fonts = CImGui.GetIO().Fonts
# default_font = CImGui.AddFontDefault(fonts)
# CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "Cousine-Regular.ttf"), 15)
# CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "DroidSans.ttf"), 16)
# CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "Karla-Regular.ttf"), 10)
# CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "ProggyTiny.ttf"), 10)
# CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "Roboto-Medium.ttf"), 16)
CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "Recursive Mono Casual-Regular.ttf"), 16)
CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "Recursive Mono Linear-Regular.ttf"), 16)
CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "Recursive Sans Casual-Regular.ttf"), 16)
CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "Recursive Sans Linear-Regular.ttf"), 16)
# @assert default_font != C_NULL

# creat texture for image drawing
img_width, img_height = 256, 256
image_id = ImGui_ImplOpenGL3_CreateImageTexture(img_width, img_height)

# setup Platform/Renderer bindings
ImGui_ImplGlfw_InitForOpenGL(window, true)
ImGui_ImplOpenGL3_Init(glsl_version)

try
    demo_open = true
    clear_color = Cfloat[0.45, 0.55, 0.60, 1.00]
    while !GLFW.WindowShouldClose(window)
        GLFW.PollEvents()
        # start the Dear ImGui frame
        ImGui_ImplOpenGL3_NewFrame()
        ImGui_ImplGlfw_NewFrame()
        CImGui.NewFrame()

        demo_open && @c ShowDemoWindow(&demo_open)

        # show image example
        CImGui.Begin("Image Demo")
        image = rand(GLubyte, 4, img_width, img_height)
        ImGui_ImplOpenGL3_UpdateImageTexture(image_id, image, img_width, img_height)
        CImGui.Image(Ptr{Cvoid}(image_id), (img_width, img_height))
        CImGui.End()

        # rendering
        CImGui.Render()
        GLFW.MakeContextCurrent(window)
        display_w, display_h = GLFW.GetFramebufferSize(window)
        glViewport(0, 0, display_w, display_h)
        glClearColor(clear_color...)
        glClear(GL_COLOR_BUFFER_BIT)
        ImGui_ImplOpenGL3_RenderDrawData(CImGui.GetDrawData())

        GLFW.MakeContextCurrent(window)
        GLFW.SwapBuffers(window)
    end
catch e
    @error "Error in renderloop!" exception=e
    Base.show_backtrace(stderr, catch_backtrace())
finally
    ImGui_ImplOpenGL3_Shutdown()
    ImGui_ImplGlfw_Shutdown()
    CImGui.DestroyContext(ctx)
    GLFW.HideWindow(window)
    GLFW.DestroyWindow(window)
end


text:
import ApproxFun: recA, recB, recC, recα, recβ, recγ
import ApproxFun: jacobirecA, jacobirecB, jacobirecC, jacobirecα, jacobirecβ, jacobirecγ

export JacobiQ, LegendreQ, WeightedJacobiQ

immutable JacobiQ{T,D<:Domain} <: RealUnivariateSpace{D}
    a::T
    b::T
    domain::D
end
LegendreQ(domain)=JacobiQ(0.,0.,domain)
LegendreQ()=LegendreQ(Interval())
JacobiQ(a,b,d::Domain)=JacobiQ(promote(a,b)...,d)
JacobiQ(a,b,d)=JacobiQ(a,b,Domain(d))
JacobiQ(a,b)=JacobiQ(a,b,Interval())

Base.promote_rule{T,V,D}(::Type{JacobiQ{T,D}},::Type{JacobiQ{V,D}})=JacobiQ{promote_type(T,V),D}
Base.convert{T,V,D}(::Type{JacobiQ{T,D}},J::JacobiQ{V,D})=JacobiQ{T,D}(J.a,J.b,J.domain)

typealias WeightedJacobiQ{T,D} JacobiQWeight{JacobiQ{T,D},D}

@compat (::Type{WeightedJacobiQ})(α,β,d::Domain)=JacobiQWeight(α,β,JacobiQ(β,α,d))
@compat (::Type{WeightedJacobiQ})(α,β)=JacobiQWeight(α,β,JacobiQ(β,α))

spacescompatible(a::JacobiQ,b::JacobiQ)=a.a==b.a && a.b==b.b

function canonicalspace(S::JacobiQ)
    #if isapproxinteger(S.a+0.5) && isapproxinteger(S.b+0.5)
    #    Chebyshev(domain(S))
    #else
        # return space with parameters in (-1,0.]
        JacobiQ(mod(S.a,-1),mod(S.b,-1),domain(S))
    #end
end

setdomain(S::JacobiQ,d::Domain)=JacobiQ(S.a,S.b,d)

for (REC,JREC) in ((:recα,:jacobirecα),(:recβ,:jacobirecβ),(:recγ,:jacobirecγ),
                   (:recA,:jacobirecA),(:recB,:jacobirecB),(:recC,:jacobirecC))
    @eval $REC{T}(::Type{T},sp::JacobiQ,k)=$JREC(T,sp.a,sp.b,k)
end

function stieltjes{T,D}(f::Fun{Jacobi{T,D}})
    g = Fun(f,Legendre(domain(f)))
    Fun(2coefficients(g),LegendreQ(domain(f)))
end
function stieltjes{D}(f::Fun{Chebyshev{D}})
    g = Fun(f,Legendre(domain(f)))
    Fun(2coefficients(g),LegendreQ(domain(f)))
end

function stieltjes{S,D}(f::Fun{JacobiWeight{S,D}})
    # Jacobi parameters need to transform to:
    α,β = f.space.α,f.space.β
    g = Fun(f,WeightedJacobi(α,β,domain(f)))
    Fun(2coefficients(g),WeightedJacobiQ(α,β,domain(f)))
end

evaluate{T,D<:Interval}(f::AbstractVector,S::JacobiQ{T,D},x) = stieltjesintervalrecurrence(S,f,tocanonical(S,x))./2jacobiQweight(S.b,S.a,tocanonical(S,x))
function evaluate{T,D<:Curve}(f::AbstractVector,S::JacobiQ{T,D},z::Number)
    sum(evaluate(f,setcanonicaldomain(S),complexroots(domain(S).curve-z)))
end
evaluate{T,D<:Curve}(f::AbstractVector,S::JacobiQ{T,D},z::AbstractArray) = reshape(promote_type(eltype(f),T,eltype(z))[ evaluate(f,S,z[i]) for i in eachindex(z) ], size(z))


text:
module Plots

if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@optlevel"))
    @eval Base.Experimental.@optlevel 1
end

const _current_plots_version = VersionNumber(
    split(
        first(
            filter(
                line -> occursin("version", line),
                readlines(normpath(@__DIR__, "..", "Project.toml")),
            ),
        ),
        "\"",
    )[2],
)

using Reexport

import GeometryBasics
using Dates, Printf, Statistics, Base64, LinearAlgebra, Random
using SparseArrays

using FFMPEG

@reexport using RecipesBase
import RecipesBase: plot, plot!, animate, is_explicit, grid
using Base.Meta
@reexport using PlotUtils
@reexport using PlotThemes
import UnicodeFun
import StatsBase
import Downloads
import Showoff
import JSON

using Requires

#! format: off
export
    grid,
    bbox,
    plotarea,
    KW,

    wrap,
    theme,

    plot,
    plot!,
    attr!,

    current,
    default,
    with,
    twinx,

    pie,
    pie!,
    plot3d,
    plot3d!,

    title!,
    xlabel!,
    ylabel!,
    xlims!,
    ylims!,
    zlims!,
    xticks!,
    yticks!,
    annotate!,
    xflip!,
    yflip!,
    xaxis!,
    yaxis!,
    xgrid!,
    ygrid!,

    xlims,
    ylims,
    zlims,

    savefig,
    png,
    gui,
    inline,
    closeall,

    backend,
    backends,
    backend_name,
    backend_object,
    aliases,

    Shape,
    text,
    font,
    stroke,
    brush,
    Surface,
    OHLC,
    arrow,
    Segments,
    Formatted,

    Animation,
    frame,
    gif,
    mov,
    mp4,
    webm,
    animate,
    @animate,
    @gif,

    test_examples,
    iter_segments,
    coords,

    translate,
    translate!,
    rotate,
    rotate!,
    center,
    BezierCurve,

    plotattr,
    scalefontsize,
    scalefontsizes,
    resetfontsizes
#! format: on
# ---------------------------------------------------------

import NaNMath # define functions that ignores NaNs. To overcome the destructive effects of https://github.com/JuliaLang/julia/pull/12563
ignorenan_minimum(x::AbstractArray{F}) where {F<:AbstractFloat} = NaNMath.minimum(x)
ignorenan_minimum(x) = Base.minimum(x)
ignorenan_maximum(x::AbstractArray{F}) where {F<:AbstractFloat} = NaNMath.maximum(x)
ignorenan_maximum(x) = Base.maximum(x)
ignorenan_mean(x::AbstractArray{F}) where {F<:AbstractFloat} = NaNMath.mean(x)
ignorenan_mean(x) = Statistics.mean(x)
ignorenan_extrema(x::AbstractArray{F}) where {F<:AbstractFloat} = NaNMath.extrema(x)
ignorenan_extrema(x) = Base.extrema(x)

# ---------------------------------------------------------

# to cater for block matrices, Base.transpose is recursive.
# This makes it impossible to create row vectors of String and Symbol with the transpose operator.
# This solves this issue, internally in Plots at least.

# commented out on the insistence of the METADATA maintainers

#Base.transpose(x::Symbol) = x
#Base.transpose(x::String) = x

# ---------------------------------------------------------

import Measures
module PlotMeasures
import Measures
import Measures:
    Length, AbsoluteLength, Measure, BoundingBox, mm, cm, inch, pt, width, height, w, h
const BBox = Measures.Absolute2DBox

# allow pixels and percentages
const px = AbsoluteLength(0.254)
const pct = Length{:pct,Float64}(1.0)

Base.convert(::Type{<:Measure}, x::Float64) = x * pct

Base.:*(m1::AbsoluteLength, m2::Length{:pct}) = AbsoluteLength(m1.value * m2.value)
Base.:*(m1::Length{:pct}, m2::AbsoluteLength) = AbsoluteLength(m2.value * m1.value)
Base.:/(m1::AbsoluteLength, m2::Length{:pct}) = AbsoluteLength(m1.value / m2.value)
Base.:/(m1::Length{:pct}, m2::AbsoluteLength) = AbsoluteLength(m2.value / m1.value)

export BBox, BoundingBox, mm, cm, inch, px, pct, pt, w, h
end

using .PlotMeasures
import .PlotMeasures: Length, AbsoluteLength, Measure, width, height
# ---------------------------------------------------------

import RecipesPipeline
import RecipesPipeline:
    SliceIt,
    DefaultsDict,
    Formatted,
    AbstractSurface,
    Surface,
    Volume,
    is3d,
    is_surface,
    needs_3d_axes,
    group_as_matrix, # for StatsPlots
    reset_kw!,
    pop_kw!,
    scale_func,
    inverse_scale_func,
    unzip,
    dateformatter,
    datetimeformatter,
    timeformatter

# Use fixed version of Plotly instead of the latest one for stable dependency
# Ref: https://github.com/JuliaPlots/Plots.jl/pull/2779
const _plotly_min_js_filename = "plotly-2.6.3.min.js"

include("types.jl")
include("utils.jl")
include("colorbars.jl")
include("axes.jl")
include("args.jl")
include("components.jl")
include("consts.jl")
include("themes.jl")
include("plot.jl")
include("pipeline.jl")
include("layouts.jl")
include("subplots.jl")
include("recipes.jl")
include("animation.jl")
include("examples.jl")
include("arg_desc.jl")
include("plotattr.jl")
include("backends.jl")
include("output.jl")
include("ijulia.jl")
include("fileio.jl")
include("init.jl")
include("legend.jl")

include("backends/plotly.jl")
include("backends/gr.jl")
include("backends/web.jl")

include("shorthands.jl")

let PlotOrSubplot = Union{Plot,Subplot}
    global title!(plt::PlotOrSubplot, s::AbstractString; kw...)                                                  = plot!(plt; title = s, kw...)
    global xlabel!(plt::PlotOrSubplot, s::AbstractString; kw...)                                                 = plot!(plt; xlabel = s, kw...)
    global ylabel!(plt::PlotOrSubplot, s::AbstractString; kw...)                                                 = plot!(plt; ylabel = s, kw...)
    global xlims!(plt::PlotOrSubplot, lims::Tuple{T,S}; kw...) where {T<:Real,S<:Real}                           = plot!(plt; xlims = lims, kw...)
    global ylims!(plt::PlotOrSubplot, lims::Tuple{T,S}; kw...) where {T<:Real,S<:Real}                           = plot!(plt; ylims = lims, kw...)
    global zlims!(plt::PlotOrSubplot, lims::Tuple{T,S}; kw...) where {T<:Real,S<:Real}                           = plot!(plt; zlims = lims, kw...)
    global xlims!(plt::PlotOrSubplot, xmin::Real, xmax::Real; kw...)                                             = plot!(plt; xlims = (xmin, xmax), kw...)
    global ylims!(plt::PlotOrSubplot, ymin::Real, ymax::Real; kw...)                                             = plot!(plt; ylims = (ymin, ymax), kw...)
    global zlims!(plt::PlotOrSubplot, zmin::Real, zmax::Real; kw...)                                             = plot!(plt; zlims = (zmin, zmax), kw...)
    global xticks!(plt::PlotOrSubplot, ticks::TicksArgs; kw...)                                                  = plot!(plt; xticks = ticks, kw...)
    global yticks!(plt::PlotOrSubplot, ticks::TicksArgs; kw...)                                                  = plot!(plt; yticks = ticks, kw...)
    global xticks!(plt::PlotOrSubplot, ticks::AVec{T}, labels::AVec{S}; kw...) where {T<:Real,S<:AbstractString} = plot!(plt; xticks = (ticks, labels), kw...)
    global yticks!(plt::PlotOrSubplot, ticks::AVec{T}, labels::AVec{S}; kw...) where {T<:Real,S<:AbstractString} = plot!(plt; yticks = (ticks, labels), kw...)
    global xgrid!(plt::PlotOrSubplot, args...; kw...)                                                            = plot!(plt; xgrid = args, kw...)
    global ygrid!(plt::PlotOrSubplot, args...; kw...)                                                            = plot!(plt; ygrid = args, kw...)
    global annotate!(plt::PlotOrSubplot, anns...; kw...)                                                         = plot!(plt; annotation = anns, kw...)
    global annotate!(plt::PlotOrSubplot, anns::AVec{T}; kw...) where {T<:Tuple}                                  = plot!(plt; annotation = anns, kw...)
    global xflip!(plt::PlotOrSubplot, flip::Bool = true; kw...)                                                  = plot!(plt; xflip = flip, kw...)
    global yflip!(plt::PlotOrSubplot, flip::Bool = true; kw...)                                                  = plot!(plt; yflip = flip, kw...)
    global xaxis!(plt::PlotOrSubplot, args...; kw...)                                                            = plot!(plt; xaxis = args, kw...)
    global yaxis!(plt::PlotOrSubplot, args...; kw...)                                                            = plot!(plt; yaxis = args, kw...)
end

# ---------------------------------------------------------

const CURRENT_BACKEND = CurrentBackend(:none)
const PLOTS_SEED = 1234

include("precompile_includer.jl")

end # module


text:
@testset "Vertex" begin

    @testset "SourceVertex" begin
        using ONNXNaiveNASflux.NaiveNASflux
        import ONNXNaiveNASflux: SourceVertex

        data = randn(3,4,5)
        sv = SourceVertex(data, "sv")

        @test name(sv) == "sv"
        @test sv() == data

        @test nout(sv) == 3
        @test nin(sv) == []
    end
end


text:
using Test
using Random
using SparseArrays: issparse, sparse
using Statistics: mean
using Transducers
using Transducers: Transducer, simple_transduce, Reduced, isexpansive,
    TeeZip, GetIndex, SetIndex, Inject, @~, IdentityTransducer,
    EmptyResultError, IdentityNotDefinedError, MissingInitError, AbortIf, @next
using InitialValues: Init
using Logging: NullLogger, with_logger

"""
    ==ₜ(x, y)

Check that _type_ and value of `x` and `y` are equal.
"""
==ₜ(_, _) = false
==ₜ(x::T, y::T) where T = x == y

prettytypeof(x) = prettytypeof(typeof(x))
prettytypeof(x::Type) = nameof(x)

inc(x) = x + oneunit(x)

isfiniteiter(iter) =
    Base.IteratorSize(iter) isa Union{Base.HasLength, Base.HasShape}

function iterator_variants(xs)
    iters = Any[xs]
    if !(xs isa Array) && isfiniteiter(xs)
        push!(iters, collect(xs))
    end
    push!(iters, Base.Generator(identity, xs))
    if !issparse(xs) && !any(ismissing, xs)
        # Currently not optimal at all, but it should work.
        if VERSION ≤ v"1.3-"
            push!(iters, sparse(xs))
        else
            # If `SparseArrays` does not support it (e.g.,
            # `sparse([""])` is not supported in Julia 1.3), don't
            # bother.  However, don't stop testing it if it is
            # re-supported again at some point.  So at least try:
            try
                push!(iters, sparse(xs))
            catch
            end
        end
    end
    return iters
end

macro expressions(ex)
    @assert ex isa Expr
    @assert ex.head == :block
    :([$([QuoteNode(a) for a in ex.args if a isa Expr]...)])
end

macro test_error(ex)
    quote
        let err = nothing
            try
                $(esc(ex))
            catch err
            end
            @test err isa Exception
            err
        end
    end
end

macro test_broken_if(cond, ex)
    ex = quote
        if $cond
            $Test.@test_broken $ex
        else
            $Test.@test $ex
        end
    end
    esc(ex)
end

macro test_inferred(ex)
    ex = quote
        $Test.@test (($Test.@inferred $ex); true)
    end
    esc(ex)
end

macro test_broken_inferred(ex)
    ex = quote
        $Test.@test_broken (($Test.@inferred $ex); true)
    end
    esc(ex)
end

function slow_test(f, title, limit)
    if lowercase(get(ENV, "TRANSDUCERS_JL_TEST_SKIP_SLOW", "false")) == "true"
        @info "Skipping $title"
        return
    end

    @info "Running $title. Limit: $limit min."

    code = """
    sleep(60)
    println("1 minute passed...")
    flush(stdout)
    for i in 2:$limit
        sleep(60)
        println("\$i minutes passed...")
        flush(stdout)
    end
    """
    # TODO: maybe terminate the parent process

    cmd = `$(Base.julia_cmd()) --startup-file=no -e $code`
    proc = run(pipeline(cmd, stdout=stdout), wait=false)

    try
        @time f()
    finally
        ch = Channel(ch -> wait(proc); csize=1)
        for _ in 1:10
            isready(ch) && break
            kill(proc)
            sleep(0.1)
        end
    end
    return
end


text:
export rescale, rescale!

@doc """
  rescale!(S::GphysData, g::Float64; chans=CC)
  rescale(S, g; chans=CC)

Rescale all channels of S to gain `g`. By default, all channels are rescaled.

  rescale!(S::GphysData; c=c, chans=CC)
  rescale(S; c=c, chans=CC)

Change `S` to use `S.gain[c]` for all channels. Rescales data as needed. By default, `c=1`.

  rescale!(St::GphysData, Ss::GphysData)
  rescale(St, Ss)

Rescale data in `St.x` to `Ss.gain` using channel ID matching; also changes `St.gain`.

  rescale!(C::GphysChannel, g::Float64)
  rescale(C, g)

Rescale `C.x` to gain `g` and set `C.gain = g`.

  rescale!(Ct::GphysChannel, Cs::GphysChannel)
  rescale(Ct, Cs)

Rescale data in `Ct.x` to `Cs.gain` and change `Ct.gain`. ID fields must match.

"""
function rescale!(S::GphysData, gain::Float64; chans::ChanSpec=Int64[])
  CC    = mkchans(chans, S, keepirr=false)
  GAIN  = getfield(S, :gain)
  X     = getfield(S, :x)
  for i in CC
    T = eltype(X[i])
    scalefac = T(gain / getindex(GAIN, i))
    if scalefac != one(T)
      rmul!(getindex(X, i), scalefac)
    end
    GAIN[i] = gain
  end
  return nothing
end
function rescale!(S_targ::GphysData, S_src::GphysData)
  N   = getfield(S_targ, :n)
  IT  = getfield(S_targ, :id)
  IS  = getfield(S_src,  :id)
  GT  = getfield(S_targ, :gain)
  GS  = getfield(S_src,  :gain)
  X   = getfield(S_targ, :x)
  for i in 1:N
    id = getindex(IT, i)
    j  = findid(id, IS)
    (j == 0) && continue
    g_new = getindex(GS, j)
    g_old = getindex(GT, i)
    if isapprox(g_new, g_old) == false
      x = getindex(X, i)
      T = eltype(x)
      rmul!(x, T(g_new/g_old))
      GT[i] = g_new
    end
  end
  return nothing
end
rescale!(S::GphysData; c::Int=1, chans::ChanSpec=Int64[]) = rescale!(S, S.gain[c], chans=chans)
function rescale!(Ct::GphysChannel, Cs::GphysChannel)
  (Ct.id == Cs.id) || error("ID mismatch!")
  gt = Ct.gain
  gs = Cs.gain
  if gt != gs
    T = eltype(Ct.x)
    rmul!(Ct.x, T(gs/gt))
    Ct.gain = gs
  end
  return nothing
end
function rescale!(C::GphysChannel, gt::Float64)
  gs = C.gain
  if gs != gt
    T = eltype(C.x)
    rmul!(C.x, T(gs/gt))
    C.gain = gt
  end
  return nothing
end

@doc (@doc rescale!)
function rescale(S::GphysData, gain::Float64; chans::ChanSpec=Int64[])
  U = deepcopy(S)
  rescale!(U, gain, chans=chans)
  return U
end
function rescale(S_targ::GphysData, S_src::GphysData)
  U = deepcopy(S_targ)
  rescale!(U, S_src)
  return U
end
function rescale(S::GphysData; c::Int=1, chans::ChanSpec=Int64[])
  U = deepcopy(S)
  rescale!(U, U.gain[c], chans=chans)
  return U
end
function rescale(Ct::GphysChannel, Cs::GphysChannel)
  U = deepcopy(Ct)
  rescale!(U, Cs)
  return U
end
rescale(C::GphysChannel, g::Float64) = (U = deepcopy(C); rescale!(U, g); return U)


text:
#File: minVar.jl
#Author:
#Date: 29-June-2020

#OBJECTIVE FUNCTION: MININIMISE VARIANCE
function minVar(W,VCV)
###########################################################################
#INPUTS:
#W      - weight vector
#VCV    - VCV matrix
###########################################################################

    fval = W'*VCV*W;
    return fval
end