text:
include(joinpath(@__DIR__, "en", "standard_dictionary_numbers_extended.jl"))
include(joinpath(@__DIR__, "en", "large_standard_dictionary_numbers_extended.jl"))

# convert a value < 100 to English.
function small_convert_en(number::Integer; british::Bool = false, dict::Symbol = :modern)
    scale_numbers = _scale_modern # define scale type
    if number < 20
        word = _small_numbers[number + 1]

        return word
    end

    v = 0
    while v < length(_tens)
        d_cap = _tens[v + 1]
        d_number = BigInt(20 + 10 * v)

        if d_number + 10 > number
            if mod(number, 10) ≠ 0
                word = d_cap * "-" * _small_numbers[mod(number, 10) + 1]

                return word
            end

            return d_cap
        end
        v += 1
    end
end

# convert a value < 1000 to english, special cased because it is the level that excludes
# the < 100 special case.  The rest are more general.  This also allows you to get
# strings in the form of "forty-five hundred" if called directly.
function large_convert_en(number::Integer; british::Bool = false, dict::Symbol = :modern)
    scale_numbers = _scale_modern # define scale type
    word = string() # initiate empty string
    divisor = div(number, 100)
    modulus = mod(number, 100)

    if divisor > 0
        word = _small_numbers[divisor + 1] * " hundred"
        if modulus > 0
            word = word * " "
        end
    end

    if british
        if ! iszero(divisor) && ! iszero(modulus)
            word = word * "and "
        end
    end

    if modulus > 0
        word = word * small_convert_en(modulus, british=british, dict=dict)
    end

    return word
end

function spelled_out_en(number::Integer; british::Bool = false, dict::Symbol = :modern)
    scale_numbers = _scale_modern # default to :modern
    if isequal(dict, :british)
        scale_numbers = _scale_traditional_british
    elseif isequal(dict, :european)
        scale_numbers = _scale_traditional_european
    elseif isequal(dict, :modern)
    else
        throw(error("unrecognized dict value: $dict"))
    end

    number = big(number)
    isnegative = false
    if number < 0
        isnegative = true
    end

    number = abs(number)
    if number > limit - 1
        throw(error("""SpelledOut.jl does not support numbers larger than $(limit_str * " - 1").  Sorry about that!"""))
    end

    if number < 100
        word = small_convert_en(number, british=british, dict=dict)

        if isnegative
            word = "negative " * word
        end

        return word
    end

    if number < 1000
        word = large_convert_en(number, british=british, dict=dict)

        if isnegative
            word = "negative " * word
        end

        return word
    end

    v = 0
    while v ≤ length(scale_numbers)
        d_idx = v
        d_number = BigInt(round(big(1000)^v))

        if d_number > number
            modulus = BigInt(big(1000)^(d_idx - 1))
            l, r = divrem(number, modulus)
            word = large_convert_en(l, british=british, dict=dict) * " " * scale_numbers[d_idx - 1]

            if r > 0
                word = word * ", " * spelled_out_en(r, british=british, dict=dict)
            end

            if isnegative
                word = "negative " * word
            end

            return word
        end

        v += 1
    end
end

function decimal_convert_en(number::AbstractFloat; british::Bool = false, dict::Symbol = :modern)
    # decimal, whole = modf(number)
    # whole = round(BigInt, whole)
    whole, decimal = split(string(number), ".")
    word = spelled_out_en(parse(BigInt, whole), british=british, dict=dict) * string(" point")
    # word = spelled_out_en(whole, british=british, dict=dict) * string(" point")

    for i in decimal
        word = word * " " * _small_number_dictionary[i]
    end

    return word
end

function spelled_out_en(number::AbstractFloat; british::Bool = false, dict::Symbol = :modern)
    str_number = format(number)
    if occursin('.', str_number)
        _length_of_presicion = length(string(first(modf(number)))) - 2
        number = parse(BigFloat, format(number, precision = _length_of_presicion)) # convert 1.01e10 to 10100000000
    else
        # It is an integer is scientific notation, treat normally without decimal precision considerations
        number = parse(BigFloat, str_number)
    end

    if isinteger(number)
        return spelled_out_en(BigInt(number), british = british, dict = dict)
    else
        return decimal_convert_en(Dec128(number), british = british, dict = dict)
    end

    # try
    #     return spelled_out_en(parse(BigInt, number), british = british, dict = dict)
    # catch
    #     return decimal_convert_en(parse(Dec128, number), british = british, dict = dict)
    # end

    return nothing
end

function spelled_out_en(number::Complex; british::Bool = false, dict::Symbol = :modern)
    return spelled_out_en(real(number), british = british, dict = dict) * " and " * spelled_out_en(imag(number), british = british, dict=dict) * " imaginaries"
end


text:
abstract type ScikitModel end

mutable struct FactorAnalysis
    k::Int
    model
    function FactorAnalysis(k)
        if !isdefined(:FactorAnalysis)
            using ScikitLearn
            @sk_import decomposition: FactorAnalysis
        end
        new(k, FactorAnalysis(k))
    end
end

ext: jl


text:
function construct_device!(ps_m::CanonicalModel,
                            device_model::DeviceModel,
                            system_formulation::Type{S},
                            system::PSY.System,
                            time_range::UnitRange{Int64};
                            kwargs...) where {S <: PM.AbstractPowerFormulation}

    construct_device!(ps_m,
                      device_model.device,
                      device_model.formulation,
                      system_formulation,
                      system,
                      time_range;
                      kwargs...)

    return

end

function InitialCondition(ps_m::CanonicalModel,
                           device::PSY.Device,
                           value::Float64)

    return InitialCondition(device,
                             PJ.add_parameter(ps_m.JuMPmodel, value))

end


text:
EΓ(x,J) = J^2/(2pi)^(3/2)*x.*exp(-x.^2/2) + 2J/(2pi)*exp(-x.^2/2) + 1/sqrt(2pi)*(1-0.5*(1+erf(x/sqrt(2))))

function Γ(Z::AbstractMatrix{Float64},u)
    m,n = size(Z)
    F,E,P = 0,0,0
    for j in 1:n-1
        for i in 1:m-1
            if Z[i,j]>u
                P += 1
                ns1 = Z[i+1,j]>u
                ns2 = Z[i+1,j+1]>u
                ns3 = Z[i,j+1]>u
                E += ns1+ns2+ns3
                F += (ns1&ns2)+(ns2&ns3)
            end
        end
        if Z[m,j]>u
            P += 1
            E += Z[m,j+1]>u
        end
    end
    for i in 1:m-1
        if Z[i,n]>u
            P += 1
            E += Z[i+1,n]>u
        end
    end
    P += Z[m,n]>u
    F-E+P
end

function Γ(Z::AbstractMatrix{Float64},u::AbstractVector)
    γs = zeros(Int,length(u))
    for i in 1:length(u)
        γs[i] = Γ(Z,u[i])
    end
    γs
end


text:
"""
    abstract type AbstractQuadratureRule{D<:AbstractReferenceShape}

A quadrature rule for integrating a function over the domain `D`.

An instance `q` of `AbstractQuadratureRule{D}` is expected to implement the
following methods:

- `q()` : return the nodes `x` and weights `w` of the quadrature rule on the
  reference domain `D`. For performance reasons, the result shoudl depend only
  on the type of `q`.
- `q(el)` : return the nodes `x` and weights `w` of the quadrature rule on the
  elemenent `D`. This assumes that `domain(q)==domain(el)`, so that the element
  quadrature can be computed by *pushing forward* a reference quadrature to `el`.
"""
abstract type AbstractQuadratureRule{D} end

domain(q::AbstractQuadratureRule{D}) where {D} = D()

"""
    qnodes(Y)

Return the quadrature nodes associated with `Y`.
"""
qnodes(q::AbstractQuadratureRule) = q()[1]

"""
    qweights(Y)

Return the quadrature weights associated with `Y`.
"""
qweights(q::AbstractQuadratureRule) = q()[2]

"""
    qnormals(Y)

Return the normal vector at the quadrature nodes of `Y.
"""
qnormals(q::AbstractQuadratureRule) = error("abstract quadrature rule has no normal")

"""
    (q::AbstractQuadratureRule)()

Return the quadrature nodes `x` and weights `w` on the `domain(q)`.
"""
function (q::AbstractQuadratureRule)() end

"""
    integrate(f,q::AbstractQuadrature)
    integrate(f,x,w)

Integrate the function `f` using the quadrature rule `q`. This is simply
`sum(f.(x) .* w)`, where `x` and `w` are the quadrature nodes and weights, respectively.
"""
function integrate(f,q::AbstractQuadratureRule)
    x,w = q()
    if domain(q) == ReferenceLine()
        return integrate(x->f(x[1]),x,w)
    else
        return integrate(f,x,w)
    end
end

function integrate(f,x,w)
    sum(zip(x,w)) do (x,w)
        f(x)*prod(w)
    end
end

# overload quadgk for integration over reference shapes. Useful for various
# testing purposes.
"""
    integrate(f,s::AbstractReferenceShape)

Use `quadgk` to (adaptively) integrate a function over the reference shape `s`.

This is used mostly for testing purposes. It returns only the value of the
integral (i.e. without the error estimate).
"""
integrate(f,l::AbstractReferenceShape;kwargs...)     = integrate(f,typeof(l);kwargs...)

integrate(f,::Type{ReferenceLine};kwargs...) = quadgk(f,0,1;kwargs...)[1]

function integrate(f,::Type{ReferenceSquare})
    I    = x-> quadgk(y->f(SVector(x,y)),0,1)[1]
    out  = quadgk(I,0,1)[1]
end

# hacky way to compute integration over reference triangles. Only used for
# testing purposes to avoid having to compute integrals analyically when testing.
function integrate(f,::Type{ReferenceTriangle})
    I    = x -> quadgk(y->f(SVector(x,y)),0,1-x)[1]
    out  = quadgk(I,0,1)[1]
end

## Define some one-dimensional quadrature rules

"""
    struct Trapezoidal{N} <: AbstractQuadratureRule{ReferenceLine}

`N`-point trapezoidal rule for integrating a function over the interval `[0,1]`.

Note that for analytic 1-periodic functions, this rule will converge
exponentially fast.

# Examples:
```julia
q    = Trapezoidal(10)
f(x) = exp(x)*cos(x)
integrate(f,q)
```
"""
struct Trapezoidal{N} <: AbstractQuadratureRule{ReferenceLine} end

Trapezoidal(n::Int) = Trapezoidal{n}()

function (q::Trapezoidal{N})() where {N}
    h = 2/(N-1)
    x = [-1.0 + k*h for k in 0:N-1]
    w = [h for k in 1:N]
    w[1]   = h/2
    w[end] = h/2
    # convert to static arrays
    xs = svector(i->SVector(0.5*(x[i]+1)),N)
    ws = svector(i->w[i]/2,N)
    return xs,ws
end

"""
    struct TrapezoidalP{N} <: AbstractQuadratureRule{ReferenceLine}

Like [`Trapezoidal{N}`](@ref), but assumes the function is 1-periodic, and so the
nodes and weights at the endSVectors of the interval `[0,1]` need not be duplicated.
"""
struct TrapezoidalP{N} <: AbstractQuadratureRule{ReferenceLine} end

TrapezoidalP(n::Int) = TrapezoidalP{n}()

# open periodic trapezoidal rule on [0,1]
function _trapezoidalP(n)
    h = 1/n
    x = [(k-0.5)*h for k in 1:n]
    w = [h   for k in 1:n]
    return x,w
end

function (q::TrapezoidalP{N})() where {N}
    x,w = _trapezoidalP(N)
    # convert to static arrays
    xs = SVector{N}(SVector{1}.(x))
    ws = SVector{N}(w)
    return xs,ws
end

"""
    struct Fejer{N}

`N`-point Fejer's first quadrature rule for integrating a function over `[0,1]`.
Exactly integrates all polynomials of degree `≤ N-1`.
"""
struct Fejer{N} <: AbstractQuadratureRule{ReferenceLine} end

Fejer(n::Int) = Fejer{n}()

function (q::Fejer{N})() where {N}
    theta = [(2j-1)*π/(2*N) for j=1:N]
    x = -cos.(theta)
    w = zero(x)
    for j in 1:N
        tmp = 0.0
        for l in 1:floor(N/2)
            tmp += 1/(4*l^2-1) *cos(2*l*theta[j])
        end
        w[j] = 2/N * (1 - 2*tmp)
    end
    xs = svector(i->SVector(0.5*(x[i]+1)),N)
    ws = svector(i->w[i]/2,N)
    return xs,ws
end

"""
    struct GaussLegendre{N}

`N`-point Gauss-Legendre quadrature rule for integrating a function over `[0,1]`.
Exactly integrates all polynomials of degree `≤ 2N-1`.
"""
struct GaussLegendre{N} <: AbstractQuadratureRule{ReferenceLine} end

GaussLegendre(n::Int) = GaussLegendre{n}()

@generated function (q::GaussLegendre{N})() where {N}
    x,w  = gauss(N) # This is a quadgk function. Gives integral in [-1,1]. Converted to [0,1] below
    xs   = svector(i->SVector(0.5*(x[i]+1)),N)
    ws   = svector(i->w[i]/2,N)
    return :($xs,$ws)
end

function refine(q::GaussLegendre{N},k=2) where {N}
    GaussLegendre(Int(N*k))
end

"""
    struct Gauss{D,N} <: AbstractQuadratureRule{D}

Tabulated `N`-point symmetric Gauss quadrature rule for integration over `D`.

This is currently implemented *by hand* for low values of `N` on triangles and tetrahedrons.
"""
struct Gauss{D,N} <: AbstractQuadratureRule{D} end

Gauss(ref,n) = Gauss{typeof(ref),n}()
Gauss(ref;n) = Gauss{typeof(ref),n}()

@generated function (q::Gauss{ReferenceTriangle,N})() where {N}
    if N == 1
        x = svector(i->SVector(1/3,1/3),1)
        w = svector(i->1/2,1)
    elseif N == 3
        x = SVector(SVector(1/6,1/6),
                    SVector(2/3,1/6),
                    SVector(1/6,2/3))
        w = svector(i->1/6,3)
    else
        notimplemented()
    end
    return :($x,$w)
end

@generated function (q::Gauss{ReferenceTetrahedron,N})() where {N}
    if N == 1
        x = SVector((SVector(1/4,1/4,1/4),))
        w = SVector(1/6)
    elseif N == 4
        a = (5-√5)/20
        b = (5+3*√5)/20
        x = SVector(SVector(a,a,a),
                    SVector(a,a,b),
                    SVector(a,b,a),
                    SVector(b,a,a)
                )
        w = svector(i->1/24,4)
    else
        notimplemented()
    end
    return :($x,$w)
end

"""
    TensorProductQuadrature{Q}

A tensor-product of one-dimension quadrature rules. Integrates over `[0,1]^d`,
where `d=length(quad)`.

# Examples
```julia
qx = Fejer(10)
qy = GaussLegendre(15)
q  = TensorProductQuadrature(qx,qy)
```
"""
struct TensorProductQuadrature{Q} <: AbstractQuadratureRule{ReferenceSquare}
    quad::Q
end

# FIXME: this is a workaround the need to easily construct a tensor quadrature
# based only on the the types of the quadratures. Useful in generated functions,
# but there is probably a better way
function TensorProductQuadrature{Tuple{Q1,Q2}}() where {Q1,Q2}
    TensorProductQuadrature(Q1(),Q2())
end

function TensorProductQuadrature(q...)
    TensorProductQuadrature(q)
end

# FIXME: the current implementation is rather obscure. How should we handle the
# product quadrature rules in general? Also make this into a generated function.
function (q::TensorProductQuadrature)()
    N       = length(q.quad)
    nodes   = map(q->q()[1],q.quad)
    weights = map(q->q()[2],q.quad)
    nodes_iter   = Iterators.product(nodes...)
    weights_iter = Iterators.product(weights...)
    x = map(x->vcat(x...),nodes_iter)
    w = map(w->prod(w),weights_iter)
    return SArray(x), w
end

# some helper functions

"""
    _qrule_for_reference_shape(ref,order)

Given a `ref`erence shape and a desired quadrature `order`, return
an appropiate quadrature rule.
"""
function _qrule_for_reference_shape(ref,order)
    if ref isa ReferenceLine
        n = ceil(Int,(order + 1) /  2)
        return GaussLegendre{n}()
    elseif ref isa ReferenceSquare
        n  = ceil(Int,(order + 1) ÷ 2)
        qx = GaussLegendre{n}()
        qy = qx
        return TensorProductQuadrature(qx,qy)
    elseif ref isa ReferenceTriangle
        if order <= 1
            return Gauss(ref,n=1)
        elseif order <=2
            return Gauss(ref,n=3)
        end
    elseif ref isa ReferenceTetrahedron
        if order <= 1
            return Gauss(ref;n=1)
        elseif order <=2
            return Gauss(ref;n=4)
        end
    end
    error("no appropriate quadrature rule found.")
end

"""
    _qrule_for_element(E,p)

Given an element type `E`, return an appropriate quadrature of order `p`.
"""
function _qrule_for_element(E,order)
    _qrule_for_reference_shape(domain(E),order)
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
"""
Tests the computation of the evaluation order using DFS
on the toy SPN.
"""
function test_compute_order_recursive()
    s, p1, p2, p3, s1, s2, s3, s4, i1, i2, i3, i4 = create_toy_spn()
    spn = SumProductNetwork(s)
    expected_order = [i1, i2, s1, i3, i4, s3, p1, s4, p2, s2, p3, s]
    for i in 1:length(expected_order)
        @test spn.order[i] == expected_order[i]
    end
end

function test_compute_order_stack()
    s, p1, p2, p3, s1, s2, s3, s4, i1, i2, i3, i4 = create_toy_spn()
    spn = SumProductNetwork(s, recursive=false)
    expected_order = [i4, i3, s4, i2, i1, s2, p3, s1, p2, s3, p1, s]
    for i in 1:length(expected_order)
        @test spn.order[i] == expected_order[i]
    end
end

"""
Tests wether the algorithms detect cycles correctly.
"""
function test_compute_order_recursive_cycles()

    function create_spn_with_cycle()
        s1 = ProdNode()
        s2 = ProdNode()
        s3 = ProdNode()
        i1 = IndicatorNode(1, 1.0)
        connect!(s1, s2)
        connect!(s2, s3)
        connect!(s3, s1)
        connect!(s3, i1)

        return s1
    end

    s1 = create_spn_with_cycle()
    @test_throws ErrorException SumProductNetwork(s1, recursive=true)

    s1 = create_spn_with_cycle()
    @test_throws ErrorException SumProductNetwork(s1, recursive=false)
end

function test_eval()
    s, p1, p2, p3, s1, s2, s3, s4, i1, i2, i3, i4 = create_toy_spn()
    spn = SumProductNetwork(s, recursive=false)

    x = [true, false]
    setInput!(spn, x)
    eval!(spn)

    @test i1.logval ≈ log(1)
    @test i2.logval ≈ log(0)
    @test i3.logval ≈ log(0)
    @test i4.logval ≈ log(1)

    @test s1.logval ≈ log(0.6)
    @test s2.logval ≈ log(0.9)
    @test s3.logval ≈ log(0.7)
    @test s4.logval ≈ log(0.8)

    @test p1.logval ≈ log(0.42)
    @test p2.logval ≈ log(0.48)
    @test p3.logval ≈ log(0.72)

    @test s.logval ≈ log(0.522)
end

function test_computeDerivatives()
    s, p1, p2, p3, s1, s2, s3, s4, i1, i2, i3, i4 = create_toy_spn()
    spn = SumProductNetwork(s, recursive=false)

    x = [true, false]
    setInput!(spn, x)
    eval!(spn)

    computeDerivatives!(spn)

    @test s.logdrv ≈ 0

    @test p1.logdrv ≈ log(0.5)
    @test p2.logdrv ≈ log(0.2)
    @test p3.logdrv ≈ log(0.3)

    @test s1.logdrv ≈ log(0.5*0.7 + 0.2*0.8)
    @test s2.logdrv ≈ log(0.3*0.8)
    @test s3.logdrv ≈ log(0.5*0.6)
    @test s4.logdrv ≈ log(0.2*0.6 + 0.3*0.9)

    @test i1.logdrv ≈ log((0.5*0.7 + 0.2*0.8)*0.6 + (0.3*0.8)*0.9)
    @test i2.logdrv ≈ log((0.5*0.7 + 0.2*0.8)*0.4 + (0.3*0.8)*0.1)
    @test i3.logdrv ≈ log((0.5*0.6)*0.3 + (0.2*0.6 + 0.3*0.9)*0.2)
    @test i4.logdrv ≈ log((0.5*0.6)*0.7 + (0.2*0.6 + 0.3*0.9)*0.8)
end

function test_normalize()
    s, p1, p2, p3, s1, s2, s3, s4, i1, i2, i3, i4 = create_toy_spn()
    spn = SumProductNetwork(s, recursive=false)

    s.weights *= 3
    s1.weights *= 5

    normalize!(spn)
    @test s.weights[1] ≈ 0.5
    @test s.weights[2] ≈ 0.2
    @test s.weights[3] ≈ 0.3
    @test s1.weights[1] ≈ 0.6
    @test s1.weights[2] ≈ 0.4
end

function test_number_of_nodes()
    s, p1, p2, p3, s1, s2, s3, s4, i1, i2, i3, i4 = create_toy_spn()
    spn = SumProductNetwork(s, recursive=false)

    @test numNodes(spn) == 12
    @test length(spn) == 12
    @test numSumNodes(spn) == 5
    @test numProdNodes(spn) == 3
    @test numLeafNodes(spn) == 4
    @test numNodes(spn, IndicatorNode) == 4
    @test numNodes(spn, GaussianNode) == 0
end

test_compute_order_recursive()
test_compute_order_stack()
test_compute_order_recursive_cycles()
test_eval()
test_computeDerivatives()
test_normalize()
test_number_of_nodes()


text:
module JLTestTest

using JLTest

@testcase false begin
  @casename "JLTest Private Impl Tests"

  @test begin
    tc = JLTest.TestContext("Dummy")
    expected = true
    result = JLTest.doTest(tc, "assertEquals", (a,b)->(a==b), 1, 1)
    assert(expected == result)
  end

  @test begin
    tc = JLTest.TestContext("Dummy")
    expected = false
    result = JLTest.doTest(tc, "assertEquals", (a,b)->(a==b), 1, 2)
    assert(expected == result)
  end

end

@testcase begin

  @casename "JLTest Tests"

  x = 0

  @setUp () -> (x+=1)
  @tearDown () -> (x=0)

  @test begin
    @testname "assertEqual Tests"
    @assertEqual(x,1)
    @assertEqual(x,1)
  end

  @test begin
    @testname "assertNotEqual Tests"
    @assertNotEqual(x,0)
  end

  @test begin
    @testname "assertLess Tests"
    @assertLess(-1,1)
  end

  @test begin
    @testname "assertLess Tests"
    @assertLessEqual(-1,1)
  end

  @test begin
    @testname "assertLessEqual Tests"
    @assertLessEqual(1,1)
  end

  @test begin
    @testname "assertGreater Tests"
    @assertGreater(1,-1)
  end

  @test begin
    @testname "assertGreaterEqual Tests"
    @assertGreaterEqual(-1,-1)
  end

  @test begin
    @testname "assertAproxEqual Tests"
    @assertAproxEqual(1,1.00005)
    @assertAproxEqual(1,1.01,0.1)
    @assertAproxEqual(1,1.01,"expected failure")
    @assertAproxEqual(1,1.01,0.1,"expected success")
    @assertAproxEqual([1.1, 2.1, 3.1], [1.10001, 2.09999, 3.1])
    @expectFailures(1)
  end

  myString1 = "a string"
  myString1a = myString1
  myString2 = "a string"
  myArray1 = [1,2,3]
  myArray1a = myArray1
  myArray2 = copy(myArray1)

  @test begin
    @testname "assertIs Tests"

    @assertIs(myString1, myString1a)
    @assertIs(myArray1, myArray1a)
  end

  @test begin
    @testname "assertIsNot Tests"
    @assertIsNot(myString1, myString2)
    @assertIsNot(myArray1, myArray2)
  end

  @test begin
    @assertIn(2, [1,2,3])
  end

  @test begin
    @assertNotIn(0, [1,2,3])
  end

  @test begin
     @assertItemsEqual([3,2,1], [1,2,3])
  end

  @test begin
    @assertIsA([3,2,1], Array)
  end

  @test begin
    @assertIsA([3,2,1], Array{Int64,1})
  end

   @test begin
     @assertIsNotA([3,2,1], Array{Int32,1})
  end

  @test begin
    @assertMatches(r"a[0-9]+a", "a1234a")
  end

  @test begin
    @assertNotMatches(r"a[0-3]+a", "a1234a")
  end

  @test begin
    a = 2
    b = 1
    @assertTrue(a==2*b)
  end

  @testskip begin
    @testname "uncond skip"
    @testFailed("Should have been skipped")
  end

  @testskip true begin
    @testname "skip with cond"
    @testFailed("Should have been skipped")
  end

  @test begin
    @testname "expected errors"
    @suppressAll 1 0 begin
      @testFailed("I expect this to fail")
    end
  end

  @test begin
    @testname "assertThrows specific exception"
    d = Dict()
    @assertThrows KeyError d["bogus"]
  end

  @test begin
    @testname "assertThrows one of specified exceptions"
    d = Dict()
    @assertThrows KeyError ErrorException error("ouch")
  end

   @test begin
    @testname "assertThrows any exception"
    d = Dict()
    @assertThrows error("boom")
  end

  @test begin
    @testname "assertThrows reports exception was not thrown"
    d = Dict()
    @suppressAll 1 0 begin
      @assertThrows ()->"I dont throw"()
    end
  end

  @test begin
    @testname "assertThrows reports exception was not thrown from a list of possible"
    d = Dict()
    @suppressAll 1 0 begin
      @assertThrows KeyError ErrorException ()->"I dont throw"()
    end
  end

  @test begin
    @testname "assertThrows reports exception was thrown but not in expected list"
    d = Dict()
    @suppressAll 1 0 begin
      @assertThrows KeyError error("I throw ErrorException")
    end
  end

  @test begin
    @testname "test msg"
    @assertEqual(1,2,"success if you see this msg")
    @expectFailures(1)
  end

  @test begin
    @testname "type related tests"
    @assertTypeOf(1, Int64)
    @assertSameType(1, 1000)
    @assertStrictlyEqual(4, 2 * 2)
    @suppressAll 1 0 begin
      @assertStrictlyEqual(4, 2.0 * 2)
    end
  end

  @testreport

end

type AssertEqualCheck
  numAssertEqual
  numAssertEqualSuccess
  numAssertEqualDiffTypes
end

@testcase begin @casename "Customization"

  @setUserData AssertEqualCheck(0,0,0)
  @setPreAssert (tc, assertion, args...)-> assertion == "assertEqual" && (tc.userData.numAssertEqual+=1;true)
  @setPostAssert function postStats(tc, assertion, args...)
    if assertion == "assertEqual"
      tc.userData.numAssertEqualSuccess  += 1
      if typeof(args[1]) != typeof(args[2])
        tc.userData.numAssertEqualDiffTypes+=1
      end
    end
  end

  @test begin
    @testname "test pre/post asserts"
    @assertEqual(1.0, 1)
    @assertEqual(1, 1)
    @quiet
    @assertEqual(1, 2) #will fail and not call post assert
    @expectFailures(1)
    @unquiet
    @assertNotEqual(1,0)
    aec = @getUserData()

    #Don't use assertEqual here unless you reset pre/post hooks!
    @assertTrue(aec.numAssertEqual == 3)
    @assertTrue(aec.numAssertEqualSuccess == 2)
    @assertTrue(aec.numAssertEqualDiffTypes == 1)
  end
end

end #module

max_stars_repo_name: smangano/JLTest

max_issues_repo_name: smangano/JLTest

max_forks_repo_name: smangano/JLTest


text:
#!/usr/bin/env julia

# Plot the 0th and 1st moment maps from the synthesized image.

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "--config"
    help = "a YAML configuration file"
    default = "config.yaml"
end

parsed_args = parse_args(ARGS, s)

import YAML
config = YAML.load(open(parsed_args["config"]))

using DiskJockey.constants
using DiskJockey.image
using DiskJockey.model
using HDF5

import PyPlot.plt
using LaTeXStrings
import Images

species = config["species"]
transition = config["transition"]
lam0 = lam0s[species*transition]
model = config["model"]
pars = convert_dict(config["parameters"], config["model"])

function plot_moment(mom::image.ZerothMoment, fname="zeroth_moment.png")

    # Image needs to be flipped along RA dimension
    ext = (mom.ra[end], mom.ra[1], mom.dec[1], mom.dec[end])

    fig, ax = plt[:subplots](nrows=1, figsize=(8,8))
    ax[:set_xlabel](L"$\Delta \alpha$ ('')")
    ax[:set_ylabel](L"$\Delta \delta$ ('')")

    # Flip the frame over the RA axis for Sky convention
    frame = flipdim(mom.data, 2)

    ax[:imshow](frame, extent=ext, interpolation="none", origin="lower", cmap=plt[:get_cmap]("PuBu"))

    plt[:savefig](fname)

end

# Assume that image has already been generated by synthesize_model.jl
im = imread()
skim = imToSky(im, pars.dpc)

# Do the velocity conversion here for plot labels
global nlam = length(skim.lams)
# convert wavelengths to velocities
global vels = c_kms * (skim.lams .- lam0)/lam0

moment = image.convert(image.ZerothMoment, skim)

plot_moment(moment)

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
# steady-state solution Kienle 2 layer 1998

@. D = 1/3musp

function fluence_spatial_freq(z, s, μa, μsp, l)

    D = @. 1/3*μsp
    α = @. sqrt((D*s^2 + μa)/D)
    A = 1
    zb = 2*A*D[1]
    z0 = 1/(μsp[1] + μa[1])

    ϕ = BigFloat(0.0)

    if z < l && z >= 0

        ϕ = 1/(D[1]*α[1]*cosh(α[1]*(l + zb)) + D[2]*α[2]*sinh(α[1]*(l + zb)))

        ϕ = ϕ*sinh(α[1]*(zb + z0))/(D[1]* α[1])
        ϕ = ϕ*(D[1]*α[1]*cosh(α[1]*(l - z)) + D[2]*α[2]*sinh(α[1]*(l - z)))

        if z < z0
            ϕ  = ϕ - sinh(α[1]*(z0 - z))/(D[1]*α[1])
        end

    else
        ϕ = sinh(α[1]*(zb + z0))*exp(α[2]*(l - z))
        ϕ = ϕ/(D[1]*α[1]*cosh(α[1]*(l + zb)) + D[2]*α[2]*sinh(α[1]*(l + zb)))
    end

    return ϕ
end

function fluence_steadystate(ρ, z, μa, μsp, l)

   ϕ = quadgk((s) -> compute_integral(z, s, μa, μsp, l, ρ), 0, 400)
   return ϕ[1]/2π

end

function compute_integral(z, s, μa, μsp, l, ρ)
    return fluence_spatial_freq(z, s, μa, μsp, l)*s*besselj0(ρ*s)
end

function fluence_spatial_freq(z, s, μa, μsp, l)
    # keep l under 10 because of bounds error
    D = @. 1/3μsp
    α = @. sqrt((D*s^2 + μa)/D)
    A = 1
    zb = 2*A*D[1]
    z0 = 1/(μsp[1] + μa[1])

    if 2*cosh(α[1]*(l + zb)) == Inf # cosh(715) returns inf
        ϕ = fluence_spatial_freq_large(z, α, z0, l, D)
    else
        ϕ = fluence_spatial_freq_small(z, α, zb, z0, l, D)
    end

    return ϕ
end

function fluence_spatial_freq_small(z, α, zb, z0, l, D)

    if z < l

        ϕ = (cosh(α[1]*(l - z))/(D[2]*α[2]) + sinh(α[1]*(l - z))/(D[1]*α[1]))
        ϕ = ϕ/(cosh(α[1]*(l + zb))/(D[2]*α[2]) + sinh(α[1]*(l + zb))/(D[1]*α[1]))
        ϕ = ϕ*sinh(α[1]*(zb + z0))/(D[1]*α[1])

        if z < z0
            ϕ  = ϕ - sinh(α[1]*(z0 - z))/(D[1]*α[1])
        end

    else
        ϕ = sinh(α[1]*(zb + z0))/(D[1]*α[1]*D[2]*α[2])
        ϕ = ϕ/(cosh(α[1]*(l + zb))/(D[2]*α[2]) + sinh(α[1]*(l + zb))/(D[1]*α[1]))
        ϕ = ϕ*exp(α[2]*(l - z))
    end

    if isinf(ϕ)
        return NaN
    end

    return ϕ
end

function fluence_spatial_freq_large(z, α, z0, l, D)

    if z < z0
        ϕ = 0.0

    elseif z < l && z >= z0
        ϕ = exp(α[1]*z0 - α[1]*z)/(2*D[1]*α[1])

    else
        ϕ = exp(α[1]*z0 + α[2]*l - α[2]*z - α[1]*l)/(D[1]*α[1] + D[2]*α[2])
    end

end

function ss_compare(ρ, z, μa, μsp)

    D = 1/3μsp
    A = 1
    z0 = 1/(μsp + μa)
    zb = 2*A*D
    μe = sqrt(3*μa*(μa + μsp))

    ϕ = exp(-μe*sqrt((z - z0)^2 + ρ^2))/sqrt((z - z0)^2 + ρ^2)
    ϕ = ϕ - exp(-μe*sqrt((z + z0 + 2*zb)^2 + ρ^2))/sqrt((z + z0 + 2*zb)^2 + ρ^2)
    ϕ = ϕ/(4*π*D)
end

function fluence_spatial_freq1(z, s, μa, μsp, l)

    D = @. 1/3μsp
    α = @. sqrt((D*s^2 + μa)/D)
    A = 1
    zb = 2*A*D[1]
    z0 = 1/(μsp[1] + μa[1])

    if z < l

        ϕ = D[1]*α[1]*(1 + exp(-2*α[1]*(l - z))) +  D[2]*α[2]*(1 - exp(-2*α[1]*(l - z)))
        ϕ = ϕ/(D[1]*α[1]*(1 + exp(-2*α[1]*(l + zb))) +  D[2]*α[2]*(1 - exp(-2*α[1]*(l + zb))))
        ϕ = ϕ*(1 + exp(-2*α[1]*(zb + z0)))/(2*D[1]*α[1])
        ϕ = ϕ*exp(α[1]*z0 - α[1]*z)

        if z < z0
            ϕ = ϕ - exp(α[1]*(z0 - z))*(1 - exp(-2*α[1]*(z0 - z)))/(2*D[1]*α[1])
        end

    else

        ϕ = exp(α[1]*(z0 - l) + α[2]*(l - z))
        ϕ = ϕ*(1 - exp(-2*α[1]*(zb + z0)))
        ϕ = ϕ/(D[1]*α[1]*(1 + exp(-2*α[1]*(l + zb))) + D[2]*α[2]*(1 - exp(-2*α[1]*(l + zb))))
    end

    return ϕ
end

ϕ = ϕ - exp(α[1]*(z0 - z))*(1 - exp(-2*α[1]*(z0 - z)))/(2*D[1]*α[1])

max_stars_repo_name: sammzer/LightPropagation

max_issues_repo_name: sammzer/LightPropagation

max_forks_repo_name: sammzer/LightPropagation


text:
using QMTK
using Compat.Test


text:
using Compat

if VERSION < v"0.3.3"
  function blasfunc(name)
    string(name)
  end
else
  function blasfunc(name)
    Base.blasfunc(name)
  end
end


text:
struct FINITO_DLFinito_iterable{R<:Real,C<:RealOrComplex{R},Tx<:AbstractArray{C},Tf,Tg} <: CIAO_iterable
    F::Array{Tf}            # smooth term
    g::Tg                   # nonsmooth term
    x0::Tx                  # initial point
    N::Int                    # of data points in the finite sum problem
    L::Maybe{Union{Array{R},R}}  # Lipschitz moduli of nabla f_i
    γ::Maybe{Union{Array{R},R}}  # stepsizes
    sweeping::Int8          # to only use one stepsize γ
    batch::Int              # batch size
    α::R                    # in (0, 1), e.g.: 0.99
    inseq::Int              # extra memory needed = length of inner seq (excluding z)                        # exm: 1231231234234234... as memory = 3
    rep::Int                # number of reps for each index
end

mutable struct FINITO_DLFinito_state{R<:Real,Tx}
    γ::Array{R}             # stepsize parameter
    hat_γ::R                # average γ
    av::Tx                  # the running average
    ind::Array{Array{Int}}  # running index set
    d::Int                  # number of batches
    z_M::Vector{Tx}
    # some extra placeholders
    ∇f_temp::Tx             # placeholder for gradients
    z_full::Tx
    inds::Array{Int}        # needed for shuffled only!
end

function FINITO_DLFinito_state(γ::Array{R}, hat_γ::R, av::Tx, ind, d, z_M) where {R,Tx} # why this function is needed?
    return FINITO_DLFinito_state{R,Tx}(
        γ,
        hat_γ,
        av,
        ind,
        d,
        z_M,
        zero(av),
        zero(av),
        collect(1:d),
    )
end

function Base.iterate(iter::FINITO_DLFinito_iterable{R}) where {R} # rewriting Base?
    N = iter.N
    r = iter.batch # batch size
    # create index sets
    ind = Vector{Vector{Int}}(undef, 0)
    d = Int(floor(N / r))
    for i = 1:d
        push!(ind, collect(r*(i-1)+1:i*r))
    end
    r * d < N && push!(ind, collect(r*d+1:N))
    # updating the stepsize
    if iter.γ === nothing
        if iter.L === nothing
            @warn "--> smoothness parameter absent"
            return nothing
        else
            γ = zeros(R, N)
            for i = 1:N
                isa(iter.L, R) ? (γ = fill(iter.α * R(iter.N) / iter.L, (N,))) :
                (γ[i] = iter.α * R(N) / (iter.L[i]))
            end
        end
    else
        isa(iter.γ, R) ? (γ = fill(iter.γ, (N,))) : (γ = iter.γ) # provided γ
    end
    #initializing the vectors
    hat_γ = 1 / sum(1 ./ γ)
    av = copy(iter.x0)
    for i = 1:N
        ∇f, ~ = gradient(iter.F[i], iter.x0)
        ∇f .*= hat_γ / N
        av .-= ∇f
    end
    # println(iter.inseq)
    z_M = [copy(av) for i = 1:(iter.inseq+1)]

    state = FINITO_DLFinito_state(γ, hat_γ, av, ind, cld(N, r), z_M)

    return state, state
end

function Base.iterate(
    iter::FINITO_DLFinito_iterable{R},
    state::FINITO_DLFinito_state{R},
) where {R}
    # full update
    prox!(state.z_full, iter.g, state.av, state.hat_γ)
    state.av .= (iter.N / state.hat_γ ) .* state.z_full
    for i = 1:iter.N
        gradient!(state.∇f_temp, iter.F[i], state.z_full) # update the gradient
        state.av .-= state.∇f_temp
    end
    state.av .*= state.hat_γ / iter.N
    iter.sweeping == 3 && (state.inds = randperm(state.d)) # shuffled

    # iter.sweeping == 3 && @error "not supported" # shuffled
    # println(state.inds)
    # repeat = 2 # number of repeats of each seq
    # println(iter.rep)
    Ind = Index_iterable(1, iter.inseq, iter.rep, iter.N)

    # Ind = take(Ind, repeat * iter.N)

    for (it, ind) in Iterators.enumerate(Ind)
        # println(it)
        if mod(it, iter.inseq * iter.rep) == 1
            # println("------new cycle starts-----")
            for ell in 1:(iter.inseq+1)       # for initializing the inner loop
                copyto!(state.z_M[ell], state.z_full)
            end
        end
        # println(ind)
        prox!(state.z_M[ind.m_n], iter.g, state.av, state.hat_γ)

        gradient!(state.∇f_temp, iter.F[state.inds[ind.i]], state.z_M[ind.m_o]) # update the gradient
        state.av .+= (state.hat_γ / iter.N) .* state.∇f_temp
        gradient!(state.∇f_temp, iter.F[state.inds[ind.i]], state.z_M[ind.m_n]) # update the gradient
        state.av .-= (state.hat_γ / iter.N) .* state.∇f_temp
        state.av .+= (state.hat_γ / state.γ[state.inds[ind.i]]) .* (state.z_M[ind.m_n] .- state.z_M[ind.m_o])
    end

    return state, state
end

# not exactly the last z since we should look at z_M[ind.m_n]
solution(state::FINITO_DLFinito_state) = state.z_M[end]

# taking prox to be able to have consistent output without knowing the last z_M cell updated
# function solution(state::FINITO_DLFinito_state)
#     sol = copy(state.z_full)
#     prox!(sol, iter.g, state.av, state.hat_γ)
#     return sol
# end

# count(state::FINITO_DLFinito_state) = []


text:
module bit
using HTTP, LibPQ, JSON, DataFrames, SQLStrings, JSONTables, CSV, URIs

# load preferences module
@static if VERSION >= v"1.6"
    using Preferences
end

@static if VERSION >= v"1.6"
    preference = @load_preference("preference", "default")
else
    preference = "default"
end

"""
    install_key!(bitio_pg_string)

Save a bit.io pg_string and password to a LocalPreferences.toml file

Make sure to add the LocalPreferences.toml file to .gitignore if in a version-controlled directory!

# Examples
```julia-repl
julia> bit.install_key!(ENV["bitio_pg_string"]) # if bitio_pg_string is in ENV
[ Info: replacing bit.io pg string
[ Info: bit.io pg_string and password saved to LocalPreferences.toml. Make sure to add this file to your .gitignore!
```
"""
function install_key!(bitio_pg_string::AbstractString)
    if @has_preference("bitio_pg_string")
        @info("replacing bit.io pg string")
    end

    # extract password from pg_string
    info = LibPQ.conninfo(bitio_pg_string)
    for option in info
        if option.keyword=="password"
           global pw = option.val
        end
    end

    #pw = match(r":(.*):(.*)@", bitio_pg_string)[2]

    # set preferences
    @set_preferences!("bitio_pg_string" => bitio_pg_string,
                      "bitio_key" => pw)

    @info("bit.io pg_string and password saved to LocalPreferences.toml. Make sure to add this file to your .gitignore!")
end

"""
    bit.query!(query, pg_string=missing)

Execute a query on bit.io. This function checks whether the pg\\_string was saved with `install_key!()`.
You can pass the `pg_string` as an optional argument if you do not want to save it to a config file.

# Example

```julia-repl
julia> bit.query!(raw"select count(*) from "\$username/\$repo"."\$tablename";")

```
"""
function query!(query; pg_string=missing)
    if ismissing(pg_string)
        pg_string = @load_preference("bitio_pg_string", missing)
        if ismissing(pg_string)
            throw(ErrorException("Please include a pg_string argument or install your pg_string with the install_key! method"))
        end
    end

    query = sql`$query`
    result = LibPQ.Connection(pg_string) do conn
        execute(conn, query.args[1])
    end
    return result
end

"""
    bit.download_table(username, schema, tablename, pg_string=missing)

Downloads a full table from a bit.io repository. This function checks whether the pg\\_string was saved with `install_key!()`.
You can pass the `pg_string` as an optional argument if you do not want to save it to a config file.

# Example

```julia-repl
julia> bit.download_table(myusername, myschemaname, myreponame)
DataFrame ...
```
"""
function download_table(username, schema, tablename; pg_string=missing)
    if ismissing(pg_string)
        pg_string = @load_preference("bitio_pg_string", missing)
        if ismissing(pg_string)
            throw(ErrorException("Please include a pg_string argument or install your pg_string with the install_key! method"))
        end
    end
    t = LibPQ.Connection(pg_string) do conn
        execute(conn, """SELECT * FROM "$username/$schema"."$tablename";""")
	end
	return DataFrame(t)
end

"""
    bit.import!(df, username, schema, tablename, bitio_key=missing,
                create_table_if_not_exists=true, if_exists="append")

import data to a bit.io table in the specified schema

This function checks whether the pg\\_string was saved with `install_key!()`.
You can pass the `pg_string` as an optional argument if you do not want to save it to a config file.

the `if_exists` argument specifies what should be done if the table aready exists. Options are:
- "append": append the data as new rows to the existing table
- "truncate": remove all existing data from the table but do not delete the table itself & leave the schema as-is
- "replace': completely delete the existing table and create a new one in its place.

# Example

```julia-repl
bit.import!(mydf, myusername, myschema, mytablename, if_exists="truncate")
HTTP.Messages.Response:
...
```

"""
function import!(df::DataFrame, username, schema, tablename; bitio_key = missing,
    create_table_if_not_exists = true, if_exists = "append")
    if ismissing(bitio_key)
        bitio_key = @load_preference("bitio_key", missing)
        if ismissing(bitio_key)
            throw(ErrorException("Please include your API key or install your pg_string with the install_key! method"))
        end
    end

    exists = query!("""SELECT EXISTS
(
    SELECT 1
    FROM information_schema.tables
    WHERE table_schema = '$username/$schema'
    AND table_name = '$tablename'
);"""
                    )
    exists = DataFrame(exists)[1,1]

    if exists
        if if_exists == "truncate"
            query!("""TRUNCATE "$username/$schema"."$tablename";""")
        elseif if_exists == "replace"
            query!("""DROP TABLE "$username/$schema"."$tablename";""")
        elseif if_exists != "append"
            throw(ErrorException("if_exists must be one of 'append', 'truncate', or 'replace'"))
        end
    end

    io = IOBuffer()
    CSV.write(io, df)

    headers = Dict(
        "Content-Disposition" => "attachment;filename='test.csv'",
        "Authorization" => "Bearer $(bitio_key)"
    )

    url_string = "https://import.bit.io/$(escapeuri(username))/$(escapeuri(schema))/$(escapeuri(tablename))"
    uri = URI(url_string)
    HTTP.request("POST", uri, headers, String(take!(io)))

end

"""Import method for CSV files"""
function import!(file::String, username, schema, tablename; bitio_key = missing,
    create_table_if_not_exists = true, if_exists = "append")
    if ismissing(bitio_key)
        bitio_key = @load_preference("bitio_key", missing)
        if ismissing(bitio_key)
            throw(ErrorException("Please include your API key or install your pg_string with the install_key! method"))
        end
    end

    exists = query!("""SELECT EXISTS
(
    SELECT 1
    FROM information_schema.tables
    WHERE table_schema = '$username/$schema'
    AND table_name = '$tablename'
);"""
    )
    exists = DataFrame(exists)[1, 1]

    if exists
        if if_exists == "truncate"
            query!("""TRUNCATE "$username/$schema"."$tablename";""")
        elseif if_exists == "replace"
            query!("""DROP TABLE "$username/$schema"."$tablename";""")
        elseif if_exists != "append"
            throw(ErrorException("if_exists must be one of 'append', 'truncate', or 'replace'"))
        end
    end

    headers = Dict(
        "Content-Disposition" => "attachment;filename='test.csv'",
        "Authorization" => "Bearer $(bitio_key)"
    )

    url_string = "https://import.bit.io/$(escapeuri(username))/$(escapeuri(schema))/$(escapeuri(tablename))"
    uri = URI(url_string)

    HTTP.request("POST", uri, headers, read(file, String))
end

end #module


text:
using Base.Test

using TALib
using TALib: basepath
using TALib: TA_MAType_SMA

using DataFrames

@testset "level 1" begin
    epsilon = 1e-10
    eps_price = 1e-6

    @test Initialize()

    for f in [GetVersionString, GetVersionMajor, GetVersionMinor, GetVersionBuild, GetVersionDate, GetVersionTime]
        print(string(f) * ": ")
        println(f())
    end

    @test length(GetVersionString()) > 3
    @test GetVersionMajor() == "0"

    s_xml = FunctionDescriptionXML()
    s_xml_expected_header = "<?xml version=\"1.0\" encoding=\"utf-8\" ?>"
    N = length(s_xml_expected_header)
    @test s_xml[1:N] == s_xml_expected_header

    angles = [0, pi/2, pi, 3pi/2, 0, pi/2, pi, 3pi/2]

    outReal = COS(angles)
    @test sum(outReal) ≈ 0 atol=epsilon
    @test outReal == cos.(angles)
    @test COS([float(pi)])[1] == -1
    @test ACOS([0.0])[1] == pi / 2
    outReal = ACOS(angles)

    outReal = SIN(angles)
    @test sum(outReal) ≈ 0 atol=epsilon
    @test outReal == sin.(angles)
    @test SIN([float(pi/2)])[1] == 1
    @test ASIN([1.0])[1] == pi / 2

    outReal = ASIN(angles)
    outReal = TAN(angles)
    outReal = ATAN(angles)

    filename = joinpath(basepath(), "ford_2012.csv")
    dfOHLCV = readtable(filename)
    dfOHLCV[:Date] = Date(dfOHLCV[:Date])
    dt = Array(dfOHLCV[:Date])
    opn = Array(dfOHLCV[:Open])
    hig = Array(dfOHLCV[:High])
    low = Array(dfOHLCV[:Low])
    cls = Array(dfOHLCV[:Close])
    price = Array(dfOHLCV[:Close])
    vol = Array(dfOHLCV[:Volume])

    @test price[1] == 11.13
    @test price[end] == 12.95

    indic = MA(price)
    @test indic[end] ≈ 11.546 atol=eps_price

    time_period = 10
    indic = MA(price, time_period=time_period, ma_type=TA_MAType_SMA)
    @test indic[end] ≈ 12.219 atol=eps_price
    @test sum(isnan.(indic)) == time_period - 1

    indic = BBANDS(price)
    time_period = 5
    @test indic[end, 1] ≈ 13.13191 atol=eps_price  # UpperBand
    @test indic[end, 2] ≈ 12.75400 atol=eps_price  # MiddleBand
    @test indic[end, 3] ≈ 12.37609 atol=eps_price  # LowerBand

    @test sum(isnan.(indic)) == (time_period - 1) * 3

    time_period = 10
    indic = BBANDS(price, time_period=time_period, deviations_up=2.0, deviations_down=2.0, ma_type=TA_MAType_SMA)
    @test indic[end, 1] ≈ 13.34468 atol=eps_price  # UpperBand
    @test indic[end, 2] ≈ 12.21900 atol=eps_price  # MiddleBand
    @test indic[end, 3] ≈ 11.09332 atol=eps_price  # LowerBand

    @test sum(isnan.(indic)) == (time_period - 1) * 3

    @test Shutdown()

end


text:
function compact_string(x)
    io = IOContext(IOBuffer(), :compact, true)
    print(io, x)
    return String(io.io)
end

function truncate_string(str, max_width)
    string_width = strwidth(str)
    if string_width <= max_width
        return str
    end

    newstring = ""
    string_width = 0
    for grapheme ∈ graphemes(str)
        grapheme_width = strwidth(grapheme)
        if string_width + grapheme_width >= max_width
            return newstring * "…"
        end
        newstring = newstring * grapheme
        string_width += grapheme_width
    end
    # Shouldn't get here, given that strwidth(s) = reduce(charwidth, 0, s)
end

# Two columns, for Series
function balance_widths(width1::Int, width2::Int, max_width::Int, min_width::Int = 20)
    @assert max_width >= 0
    biggest_width = max(width1, width2)

    while width1 + width2 > max_width
        biggest_width -= 1
        if biggest_width < min_width
            break
        end
        width1 = min(biggest_width, width1)
        width2 = min(biggest_width, width2)
    end

    return (width1, width2)
end

# Multiple columns, for Table
function balance_widths!(widths::Vector{Int}, ncols::Int, max_width::Int, buffer::Int = 2, min_width::Int = 20)
    @assert max_width >= 0
    biggest_width = maximum(widths)

    # First pass: shorten all the column widths
    while sum(widths) + buffer*length(widths) > max_width
        biggest_width -= 1
        if biggest_width < min_width
            break
        end
        for i = 1:length(widths)
            @inbounds widths[i] = min(biggest_width, widths[i])
        end
    end

    # Second pass: truncate columns, if necessary
    if sum(widths) + buffer*length(widths) > max_width
        total_width = 0
        for i = 1:length(widths)
            width = widths[i] + buffer
            if i == length(widths)
                ncols_shown = i-1
            else
                if total_width + width + buffer + 1 > max_width # we anticipate at least a "⋯" column to follow this one
                    ncols_shown = i-1
                    break
                end
            end
            total_width += width
        end
    else
        ncols_shown = length(widths)
    end

    return ncols_shown
end

show(io::IO, ::MIME"text/plain", s::AbstractSeries) = show_series(io, s)
show(io::IO, s::AbstractSeries) = show_series(io, s)
function show_series(io::IO, s::ANY)
    inds = collect(indices(s)[1])::Vector
    nrows = length(inds)::Int

    if nrows == 0
        print(io, "0-element Series")
        return
    end

    print(io, nrows, "-element Series:")

    max_show_rows = displaysize(io)[1] - 5

    if nrows <= max_show_rows
        ind_strings = map(compact_string, inds)::Vector{String}
        val_strings = map(i -> compact_string(s[i]), inds)::Vector{String}
    else
        ind_strings = String[compact_string(inds[i]) for i ∈ 1:max_show_rows]
        val_strings = String[compact_string(s[inds[i]]) for i ∈ 1:max_show_rows]
    end

    max_ind_len = maximum(strwidth, ind_strings)
    max_val_len = maximum(strwidth, val_strings)

    # Make sure fits horizontally on one line on screen
    max_show_width = displaysize(io)[2] - 4 # proceeding space plus 3 for seperator
    if max_ind_len + max_val_len > max_show_width
        (max_ind_len, max_val_len) = balance_widths(max_ind_len, max_val_len, max_show_width)
        map!(str -> truncate_string(str, max_ind_len), ind_strings, ind_strings)
        map!(str -> truncate_string(str, max_val_len), val_strings, val_strings)
    end

    if nrows > max_show_rows
        push!(ind_strings, (" " ^ div(max_ind_len-1, 2)) *  "⋮")
        push!(val_strings, (" " ^ div(max_val_len-1, 2)) *  "⋮")
    end

    for i ∈ 1:length(ind_strings)
        print(io, "\n ")
        @inbounds ind_str = ind_strings[i]
        print(io, ind_str)
        n_spaces = max_ind_len - strwidth(ind_str)
        if n_spaces > 0
            print(io, " " ^ n_spaces)
        end
        print(io, " │ ")
        @inbounds val_str = val_strings[i]
        print(io, val_str)
    end
end

show(io::IO, ::MIME"text/plain", t::AbstractTable) = show_table(io, t)
show(io::IO, t::AbstractTable) = show_table(io, t)
function show_table(io::IO, t::ANY)
    inds = indices(t)
    row_inds = collect(inds[1])::Vector
    col_inds = collect(inds[2])::Vector
    nrows = length(row_inds)::Int
    ncols = length(col_inds)::Int
    display_width = displaysize(io)[2]
    max_cols = display_width ÷ 3 # assuming one-width columns with two spaces in-between

    if nrows == 0 | ncols == 0
        # Perhaps we can handle one or the other being 0?
        print(io, "$nrows×$ncols Table")
        return
    end

    print(io, nrows, "×", ncols, " Table:")

    max_show_rows = displaysize(io)[1] - 7

    ncols_shown = min(ncols+1, max_cols) # First "column" shown is the indices

    strings = Vector{Vector{String}}(ncols_shown)
    for i ∈ 1:ncols_shown
        for j ∈ 1:(min(max_show_rows, nrows) + 1)
            if j == 1
                if i == 1
                    strings[i] = [""]
                else
                    strings[i] = [compact_string(col_inds[i-1])]
                end
            else
                if i == 1
                    push!(strings[i], compact_string(row_inds[j-1]))
                else
                    push!(strings[i], compact_string(t[row_inds[j-1], col_inds[i-1]]))
                end
            end
        end
    end

    max_column_widths = [max(1, maximum(strwidth, str_vec)) for str_vec ∈ strings]

    if sum(max_column_widths) + 2*length(max_column_widths) > display_width
        # Shorten each column and reduce the number of columns shown
        # If they all still don't fit, the final shown column will be filled with "⋯"

        ncols_shown = balance_widths!(max_column_widths, ncols, display_width)
        strings = strings[1:ncols_shown]
        max_column_widths = max_column_widths[1:ncols_shown]
        for i = 1:length(strings)
            map!(str -> truncate_string(str, max_column_widths[i]), strings[i], strings[i])
        end
    end

    if ncols_shown < ncols+1
        push!(strings, fill("⋯", length(strings[1])))
        push!(max_column_widths, 1)
    end

    if nrows > max_show_rows
        for i ∈ 1:length(strings)
            if i == length(strings) && strings[i][1] == "⋯"
                push!(strings[i], (" " ^ div(max_column_widths[i]-1, 2)) *  "⋱")
            else
                push!(strings[i], (" " ^ div(max_column_widths[i]-1, 2)) *  "⋮")
            end
        end
    end

    # Now produce output

    # Header: " --  Column1  Column2"
    print(io, "\n ")
    n_spaces = max_column_widths[1] + 3
    print(io, " " ^ n_spaces)
    for i ∈ 2:length(strings)
        @inbounds col_str = strings[i][1]
        print(io, col_str)
        if i != length(strings)
            n_spaces = max_column_widths[i] - strwidth(col_str) + 2
            print(io, " " ^ n_spaces)
        end
    end

    # Seperator: " ┌────────"
    print(io, "\n ")
    print(io, " " ^ max_column_widths[1])
    print(io, " ┌─")
    print(io, "─" ^ (sum(max_column_widths) - max_column_widths[1] + 2*(length(max_column_widths) - 2)))

    # Body " rowind │ val1  val2"
    for j = 2:length(strings[1])
        print(io, "\n ")
        @inbounds row_str = strings[1][j]
        print(io, row_str)
        n_spaces = max_column_widths[1] - strwidth(row_str)
        if n_spaces > 0
            print(io, " " ^ n_spaces)
        end
        print(io, " │ ")

        for i = 2:length(strings)
            @inbounds val_str = strings[i][j]
            print(io, val_str)

            if i != length(strings)
                n_spaces = max_column_widths[i] - strwidth(val_str) + 2
                print(io, " " ^ n_spaces)
            end
        end
    end
end


text:
# Copyright (c) 2014, Dreki Þórgísl <dreki@billo.systems>
#               2014, Bence Golda <bence@cursorinsight.com>
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
#  * Redistributions of source code must retain the above copyright notice,
#    this list of conditions and the following disclaimer.
#  * Redistributions in binary form must reproduce the above copyright notice,
#    this list of conditions and the following disclaimer in the documentation
#    and/or other materials provided with the distribution.
#  * Neither the name of the copyright holders nor the names of its
#    contributors may be used to endorse or promote products derived from this
#    software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
module Decode

using ErlPort.Exceptions

export decode, decodeterm, decodeatom,
decodesmallint, decodeint, decodenewfloat,
decodesmallbigint, decodelargebigint,
decodebin,
decodenil, decodestring, decodelist,
decodesmalltuple, decodelargetuple,
decompressterm,
size1unpack, size2unpack, size4unpack # for tests only XXX: do we really need this here?

include("Tags.jl")
include("Util.jl")

function decode(bytes::Array{Uint8,1})
    lencheck(bytes, 1)
    if bytes[1] != version
        throw(UnknownProtocolVersion(bytes[1]))
    end
    if length(bytes) >= 2 && bytes[2] == compressedtag
        # XXX maybe have this match the call to decode below? bytes[2:end]
        # instead of just bytes?
        return decodeterm(decompressterm(bytes))
    end
    return decodeterm(bytes[2:end])
end

function decode(unsupported)
    throw(UnsupportedType(unsupported))
end

function decodeterm(bytes::Array{Uint8,1})
    lencheck(bytes, 1)
    tag = bytes[1]
    if tag == atomtag
        return decodeatom(bytes)
    elseif tag == niltag
        return decodenil(bytes)
    elseif tag == stringtag
        return decodestring(bytes)
    elseif tag == smalltupletag
        return decodesmalltuple(bytes)
    elseif tag == largetupletag
        return decodelargetuple(bytes)
    elseif tag == listtag
        return decodelist(bytes)
    elseif tag == smallinttag
        return decodesmallint(bytes)
    elseif tag == inttag
        return decodeint(bytes)
    elseif tag == bintag
        return decodebin(bytes)
    elseif tag == newfloattag
        return decodenewfloat(bytes)
    elseif tag == smallbiginttag
        return decodesmallbigint(bytes)
    elseif tag == largebiginttag
        return decodelargebigint(bytes)
    else
        throw(UnsupportedData(bytes))
    end
end

function decodeterm(acc::Array, byte::Uint8)
    vcat(acc, decodeterm([byte]))
end

function decodeatom(bytes::Array{Uint8,1})
    len = lencheck(bytes, 3)
    unpackedlen = lencheck(len, size2unpack(bytes[2:3]) + 3, bytes)
    name = bytes[4:unpackedlen]
    if name == b"true"
        return (true, bytes[unpackedlen+1:end])
    elseif name == b"false"
        return (false, bytes[unpackedlen+1:end])
    elseif name == b"undefined"
        return (nothing, bytes[unpackedlen+1:end])
    else
        return (symbol(name), bytes[unpackedlen+1:end])
    end
end

function decodenil(bytes::Array{Uint8,1})
    lencheck(bytes, 1)
    return ([], bytes[2:end])
end

function decodestring(bytes::Array{Uint8,1})
    len = lencheck(bytes, 3)
    unpackedlen = lencheck(len, size2unpack(bytes[2:3]) + 3, bytes)
    (bytes[4:unpackedlen], bytes[unpackedlen+1:end])
end

function decodesmallint(bytes::Array{Uint8,1})
    lencheck(bytes, 2)
    (int1unpack(bytes[2]), bytes[3:end])
end

function decodeint(bytes::Array{Uint8,1})
    lencheck(bytes, 5)
    (int4unpack(bytes[2:5]), bytes[6:end])
end

function decodebin(bytes::Array{Uint8,1})
    len = lencheck(bytes, 5)
    unpackedlen = lencheck(len, size4unpack(bytes[2:5]) + 5, bytes)
    (bytes[6:unpackedlen], bytes[unpackedlen+1:end])
end

function decodenewfloat(bytes::Array{Uint8,1})
    lencheck(bytes, 9)
    (floatunpack(bytes[2:9]), bytes[10:end])
end

function decodelist(bytes::Array{Uint8,1})
    lencheck(bytes, 5)
    (results, tail) = converttoarray(size4unpack(bytes[2:5]), bytes[6:end])
    # XXX mojombo's BERT (https://github.com/mojombo/bert) does the same -- it
    # skips the improper part in lists (or throws a RuntimeError)
    (skipped, tail) = decodeterm(tail)
    (results, tail)
end

function converttoarray(len::Uint64, tail::Array{Uint8,1})
    results = map([0:1:len-1]) do i
        (term, tail) = decodeterm(tail)
        term
    end
    (results, tail)
end

function decodesmalltuple(bytes::Array{Uint8,1})
    lencheck(bytes, 2)
    converttotuple(size1unpack(bytes[2]), bytes[3:end])
end

function decodelargetuple(bytes::Array{Uint8,1})
    lencheck(bytes, 5)
    converttotuple(size4unpack(bytes[2:5]), bytes[6:end])
end

function converttotuple(len::Uint64, tail::Array{Uint8,1})
    if len < 1
        return( (), tail )
    end
    (results, tail) = converttoarray(len, tail)
    (tuple(results...), tail)
end

function decodesmallbigint(bytes::Array{Uint8,1})
    len = lencheck(bytes, 3)
    bisize = size1unpack(bytes[2])
    lencheck(len, bisize + 3, bytes)
    result = computebigint(bisize, bytes[4:bisize+3], bytes[3])
    (result, bytes[bisize+4:end])
end

function decodelargebigint(bytes::Array{Uint8,1})
    len = lencheck(bytes, 6)
    bisize = size4unpack(bytes[2:5])
    lencheck(len, bisize + 6, bytes)
    result = computebigint(bisize, bytes[7:bisize+6], bytes[6])
    (result, bytes[bisize+7:end])
end

function computebigint(len::Uint64, coefficients::Array{Uint8,1}, sign::Uint8)
    result = sum((256 .^ [0:len-1]) .* convert(Array{Int}, coefficients))
    return(sign > 0 ? -result : result)
end

end


text:
module TestTrap
using Debug
import Debug.Scope, Debug.Node, Debug.@instrument

firstline = -1
function trap(node::Node, scope::Scope)
    global firstline = (firstline == -1) ? node.loc.line : firstline
    line = node.loc.line - firstline + 1

    print(line, ":")
    if (line >  1); print("\tx = ", scope[:x]) end
    if (line == 3); print("\tk = ", scope[:k]) end
    println()
end

@instrument trap function f(n)
    x = 0       # line 1
    for k=1:n   # line 2
        x += k  # line 3
    end         # line 4
    x = x*x     # line 5
    x           # line 6
end

f(3)

errtrap(node, scope) = @assert false
@instrument trap @notrap begin; @bp; x=1; end

# #71
@debug begin
   serialize(IOBuffer(), @notrap ()->1)
end

@notrap z=5
@assert z == 5

end


text:
# parallel communication primatives

include("tag_manager.jl")  # managment of MPI tags
include("parallel_types.jl")  # datatype
include("parallel_rev1.jl")
include("parallel_rev2.jl")

"""
  This function is a thin wrapper around exchangeData().  It is used for the
  common case of sending and receiving the solution variables to other processes.
  It uses eqn.shared_data to do the parallel communication.
  eqn.shared_data *must* be passed into the corresponding finishDataExchange
  call.

  It is safe to call this function in the serial case.

  **Inputs**

   * mesh: an AbstractMesh
   * sbp: an SBP operator
   * eqn: an AbstractSolutionData
   * opts: options dictionary

  **Keyword arguments**

   * wait: wait for sends and receives to finish before exiting
"""
function startSolutionExchange(mesh::AbstractMesh, sbp::AbstractOperator,
                                  eqn::AbstractSolutionData, opts;
                                  wait=false)

  if mesh.npeers == 0
    return nothing
  end

  pdata = eqn.shared_data[1].pdata  # assume all are same
  if pdata == PARALLEL_DATA_FACE
    populate_buffer = getSendDataFace
  elseif pdata == PARALLEL_DATA_ELEMENT
    populate_buffer = getSendDataElement
  else
    throw(ErrorException("unsupported parallel_type = $(getParallelDataString(pdata))"))
  end

  exchangeData(mesh, sbp, eqn, opts, eqn.shared_data, populate_buffer, wait=wait)

  return nothing
end

"""
  This function posts the MPI sends and receives for a vector of SharedFaceData.  It works for both `PARALLEL_DATA_FACE or `PARALLEL_DATA_ELEMENT.  The only
  difference between these two cases is the populate_buffer() function.

  The previous receives using these SharedFaceData objects should have
  completed by the time this function is called.  An exception is throw
  if this is not the case.

  The previous sends are likely to have completed by the time this function
  is called, but they are waited on if not.  This function might not perform
  well if the previous sends have not completed.
  #TODO: fix this using WaitAny

  Inputs:
    mesh: an AbstractMesh
    sbp: an SBPOperator
    eqn: an AbstractSolutionData
    opts: the options dictionary
    populate_buffer: function with the signature:
                     populate_buffer(mesh, sbp, eqn, opts, data::SharedFaceData)
                     that populates data.q_send
  Inputs/Outputs:
    shared_data: vector of SharedFaceData objects representing the parallel
                 communication to be done

  Keyword Arguments:

    wait: wait for the sends and receives to finish before returning.  This
          is a debugging option only.  It will kill parallel performance.
"""
function exchangeData(mesh::AbstractMesh, sbp::AbstractOperator,
                      eqn::AbstractSolutionData, opts,
                      shared_data::Vector{SharedFaceData{T}},
                      populate_buffer::Function;
                      wait=false) where T

  npeers = length(shared_data)

  # bail out early if there is no communication to do
  # not sure if the rest of this function runs correctly if npeers == 0
  if npeers == 0
    return nothing
  end

  # this should already have happened.  If it hasn't something else has
  # gone wrong in the solver.  Throw an exception
  assertReceivesWaited(shared_data)

  # post the receives first
  for i=1:npeers
    Irecv!(shared_data[i])
  end

  # verify the sends are consistent
  assertSendsConsistent(shared_data)

  for i=1:npeers
    # wait for these in order because doing the waitany trick doesn't work
    # these should have completed long ago, so it shouldn't be a performance
    # problem

    # the waitany trick doesn't work because this loop posts new sends, reusing
    # the same array of MPI_Requests.

    # TODO: use 2 arrays for the Requests: old and new, so the WaitAny trick
    #       works

    idx = i
    data_i = shared_data[idx]
    tag = data_i.tag

    # wait on the previous send if it hasn't been waited on yet
    if !data_i.send_waited
      waitSend(data_i)
    end

    # move data to send buffer
    populate_buffer(mesh, sbp, eqn, opts, data_i)

    # post the send
    Isend(data_i)
  end

  if wait
    waitAllSends(shared_data)
    waitAllReceives(shared_data)
  end

  return nothing
end

"""
  This is the counterpart of exchangeData.  This function finishes the
  receives started in exchangeData.

  This function (efficiently) waits for a receive to finish and calls
  a function to do calculations for on that data. If `PARALLEL_DATA_FACE`
  == PARALLEL_DATA_FACE, it also permutes the data in the receive buffers to agree
  with the ordering of elementL.  For `PARALLEL_DATA_ELEMENT,
  users should call SummationByParts.interiorFaceInterpolate to interpolate
  the data to the face while ensuring proper permutation.

  Calling this function repeatedly without calling [`exchangeData`](@ref)
  in between is supported.  In this case, it only waits for the MPI
  communication to finish the first time.

  **Inputs**

   * mesh: an AbstractMesh
   * sbp: an SBPOperator
   * eqn: an AbstractSolutionData
   * opts: the options dictonary
   * calc_func: function that does calculations for a set of shared faces
               described by a single SharedFaceData.  It must have the signature
               calc_func(mesh, sbp, eqn, opts, data::SharedFaceData)

  **Inputs/Outputs**

   * shared_data: vector of SharedFaceData, one for each peer process that
                 needs to be communicated with.  By the time calc_func is
                 called, the SharedFaceData passed to it has its q_recv field
                 populated.  See note above about data permutation.
"""
function finishExchangeData(mesh, sbp, eqn, opts,
                            shared_data::Vector{SharedFaceData{T}},
                            calc_func::Function) where T

  npeers = length(shared_data)

  if mesh.npeers == 0
    return nothing
  end

  val = assertReceivesConsistent(shared_data)

  for i=1:npeers
    if val == 0  # request have not been waited on previously
      eqn.params.time.t_wait += @elapsed idx = waitAnyReceive(shared_data)
    else
      idx = i
    end

    data_idx = shared_data[idx]
    if data_idx.pdata == PARALLEL_DATA_FACE && val == 0
      # permute the received nodes to be in the elementR orientation
      permuteinterface!(mesh.sbpface, data_idx.interfaces, data_idx.q_recv)
    end

    calc_func(mesh, sbp, eqn, opts, data_idx)

  end

  return nothing
end

@doc """
### Utils.verifyCommunication

  This function checks the data provided by the Status object to verify a
  communication completed successfully.  The sender's rank and the number of
  elements is checked agains the expected sender and the buffer size

  Inputs:
    data: a SharedFaceData
"""->
function verifyReceiveCommunication(data::SharedFaceData{T}) where T
# verify a communication occured correctly by checking the fields of the
# Status object
# if the Status came from a send, then peer should be comm_rank ?

  sender = MPI.Get_source(data.recv_status)
  @assert sender == data.peernum

  ndata = MPI.Get_count(data.recv_status, T)
  @assert ndata == length(data.q_recv)

  return nothing
end

"""
  This function populates the send buffer from `eqn.q` for
  `PARALLEL_DATA_FACE`

  Inputs:
    mesh: a mesh
    sbp: an SBP operator
    eqn: an AbstractSolutionData
    opts: options dictonary

  Inputs/Outputs:
    data: a SharedFaceData.  data.q_send will be overwritten
"""
function getSendDataFace(mesh::AbstractMesh, sbp::AbstractOperator,
                         eqn::AbstractSolutionData, opts, data::SharedFaceData)

  idx = data.peeridx
  bndryfaces = data.bndries_local
  boundaryinterpolate!(mesh.sbpface, bndryfaces, eqn.q, data.q_send)

  return nothing
end

"""
  This function populates the send buffer from eqn.q for
  `PARALLEL_DATA_ELEMENT`

  Inputs:

    mesh: a mesh
    sbp: an SBP operator
    eqn: an AbstractSolutionData
    opts: options dictonary

  Inputs/Outputs:

    data: a SharedFaceData.  data.q_send will be overwritten
"""
function getSendDataElement(mesh::AbstractMesh, sbp::AbstractOperator,
                         eqn::AbstractSolutionData, opts, data::SharedFaceData)

  # copy data into send buffer
  idx = data.peeridx
  local_els = mesh.local_element_lists[idx]
  send_buff = data.q_send
  for j=1:length(local_els)
    el_j = local_els[j]
    for k=1:size(eqn.q, 2)
      for p=1:size(eqn.q, 1)
        send_buff[p, k, j] = eqn.q[p, k, el_j]
      end
    end
  end

  return nothing
end

@doc """
### Utils.mpi_master

  This macro introduces an if statement that causes the expression to be
  executed only if the variable myrank is equal to zero.  myrank must exist
  in the scope of the caller

"""->
macro mpi_master(ex)
  return quote
#    println("myrank = ", esc(myrank))
    if $(esc(:(myrank == 0)))
      $(esc(ex))
    end
  end
end

@doc """
### Utils.time_all

  This macro returns the value produced by the expression as well as
  the execution time, the GC time, and the amount of memory allocated
"""->
macro time_all(ex)
  quote
    local stats = Base.gc_num()
    local elapsedtime = time_ns()
    local val = $(esc(ex))
    elapsedtime = time_ns() - elapsedtime
    local diff = Base.GC_Diff(Base.gc_num(), stats)
    (val, elapsedtime/1e9, diff.total_time/1e9, diff.allocd)
  end
end

function print_time_all(f, t_elapsed, t_gc, alloc)
    println(f, t_elapsed, " seconds, ", t_gc, " GC seconds, ", alloc, " bytes allocated")
end

#------------------------------------------------------------------------------
# Debugging functions

"""
  Verify the contents of the receive buffers are current by doing a parallel
  communication and comparing the new contents of the buffers with the old
  ones.

  Throws an error if contents are not consistent.

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
   * f: IO, defaults to STDOUT
"""
function checkBufferConsistency(mesh, sbp, eqn::AbstractSolutionData{Tsol}, opts, f::IO=STDOUT) where Tsol

  println(f, "\nChecking buffer consistency")

  # make sure communication is finished
  assertReceivesWaited(eqn.shared_data)

  # copy buffers
  nbufs = length(eqn.shared_data)
  old_bufs = Array{Array{Tsol, 3}}(nbufs)
  for i=1:nbufs
    old_bufs[i] = copy(eqn.shared_data[i].q_recv)
    println(f, "initially, norm of buffer ", i, " = ", vecnorm(old_bufs[i]))
  end

  # do the parallel communications
  startSolutionExchange(mesh, sbp, eqn, opts)

  # finish the communication
  finishExchangeData(mesh, sbp, eqn, opts, eqn.shared_data, tmpCalcFunc)

  # compare buffers
  for i=1:nbufs
    buf1 = old_bufs[i]
    buf2 = eqn.shared_data[i].q_recv
    nrm_i = vecnorm(buf1 - buf2)
    println(f, "buffer ", i, " diffnorm = ", nrm_i)
    println(f, "buffer ", i, " diffnorm real real = ", vecnorm(real(buf1) - real(buf2)))
    println(f, "buffer ", i, " diffnorm real complex = ", vecnorm(real(buf1) - buf2))
    println(f, "buffer ", i, " diffnorm complex real = ", vecnorm(buf1 - real(buf2)))
#    @assert abs(nrm_i) < 1e-13
  end

  return nothing
end

# used for debugging
function tmpCalcFunc(mesh, sbp, eqn, opts, data::SharedFaceData)

  return nothing
end


text:
# Sample script for converting VLSV time series files to VTK image files.
#
# Usage:
#   julia -t 2 demo_convert2vti.jl
# or
#   JULIA_NUM_THREADS=2 julia demo_convert2vti.jl
#
# Hongyang Zhou, hyzhou@umich.edu

using Vlasiator, Glob

files = glob("*.vlsv")

Threads.@threads for i in 1:length(files)
   file = files[i]
   @info file, Threads.threadid()
   write_vtk(file; vti=true)
end


text:
#=
Cross-correlation plot for the paper.
Depending on the `absolute` variable, compare velocity with MTD or with
abs(MTD).

Integrates with DrWatson to not re-produce the data. Change `force` to
true to re-calculate correlations.
=#
using DrWatson
quickactivate(@__DIR__, "VelocityMTD")
using MusicManipulations, PyPlot, StatsBase, Statistics, Parameters
include(srcdir("correlations.jl"))
include(srcdir("extract_data.jl"))
include(srcdir("style.jl"))

datasets = ["tapping", "uwe", "pgmusic", "playalongs", "playalongs"]
datanames = ["tapping",  "Uwe Meile (piano)", "PG music (piano)", "drums playalongs", "drums playal. (lead)"]
subsets = ["all", "melody", "melody", "all", "lead"]
absolute = true

pcount = 0
fig = figure(figsize = (2figx, 3figx/4))
axρ = subplot(311)
axm = subplot(312)
axmn = subplot(313)
xlocs = zeros(length(datasets))
total = zeros(Int, length(datasets))
significant = zeros(Int, length(datasets))

for i in 1:length(datasets)
    # i == 3 && continue
    dataset = datasets[i]
    subset = subsets[i]
    c = @dict dataset subset absolute

    file, = produce_or_load(
        datadir("crosscor"), c, correlation_intervals;
        prefix = "crosscor",
        force = false
    )

    @unpack ρs, ρus, ρds, MIs, centers, MIus, MIds = file

    L = length(ρs)
    l = pcount + 1
    r = l + L - 1

    axρ.errorbar(l:r, ρs, yerr = hcat(ρus,ρds)',
    ecolor = "C0", color = "C$(i)", ls = "None", marker = "o")

    bot = mean(ρs) - mean(ρds)
    top = mean(ρs) + mean(ρus)
    rec = matplotlib.patches.Rectangle((l, bot), r-l, top-bot, color = "C$(i)",
    alpha = 0.2)
    # axρ.add_patch(rec)

    axρ.plot([l, r], [mean(ρs), mean(ρs)], color = "C$(i)", lw = 2, alpha = 0.75)

    # Plot default MI
    axm.plot([l, r], fill(mean(MIs),2), color = "C$(i)",
    lw = 2, alpha = 0.75)

    axm.errorbar(l:r, centers, yerr = hcat(MIus,MIds)',
    ecolor = "C0", color = "C0", ls = "None")
    axm.plot(l:r, MIs,  color = "C$(i)", ls = "None",
    marker = "o")
    y = 0.1
    # axm.axhline(0, color = "black", ls = "dashed", lw = 1.0)
    axm.set_ylim(0, 0.9)

    # plot normalized MI
    # make normalized mi
    no = zeros(length(MIs))
    for i in 1:length(no)
        mic = MIs[i] - centers[i]
        no[i] = if mic < 0
            mic/(MIds[i])
        else
            mic/(MIus[i])
        end
        no[i] = abs(no[i])
    end
    # plot it
    axmn.plot(l:r, no,  color = "C$(i)", ls = "None",
    marker = "o")
    # average
    axmn.plot([l, r], fill(mean(no), 2), color = "C$(i)",
    lw = 2, alpha = 0.75)
    axmn.axhline(1, color = "black", ls = "dashed", lw = 0.5)
    y = 3
    axmn.set_ylim(0, 4)

    axmn.text((r+l)/2, y, datanames[i], ha="center", va="center", color = "C$i")

    global pcount += L+1
    total[i] = length(no)
    significant[i] = count(n -> n > 1, no)
end

axρ.axhline(0, color = "black", ls = "dashed", lw = .5)
axρ.set_xticklabels([])
axρ.set_yticks(-0.3:0.2:0.3)
axρ.set_ylabel("\$\\rho\$")
axρ.set_title("correlation between velocity and $(absolute ? "abs(MTD)" : "MTD")")

axm.set_yticks(0:0.2:0.8)
axm.set_xticklabels([])
axm.set_ylabel("MI")

axmn.set_yticks(0:4)
axmn.set_ylabel("MI (norm.)")
axmn.set_xlabel("recording #")

#
for ax in fig.axes
    ax.minorticks_on()
    ax.tick_params(axis="y", which="minor", left=false)
    ax.set_xlim(0, pcount)
    ax.tick_params(axis="x", width =2, length=10)
    ax.tick_params(axis="x", which="minor", width = 1.5, length = 8)
end
add_identifiers!(fig; loc = (0.995, 0.975))
fig.align_ylabels()
fig.tight_layout()
fig.subplots_adjust(hspace = 0.2)
# fig.savefig(papersdir("figs", "crosscorrelations"))
fig.savefig(plotsdir("crosscorrelations", "absolute=$(absolute)"))

max_stars_repo_name: Datseris/VelocityMTD

max_issues_repo_name: Datseris/VelocityMTD

max_forks_repo_name: Datseris/VelocityMTD


text:
i2 = 2
f2 = 2.0
d2 = Double64(2.0)
b2 = BigFloat(2)
r2 = 2//1
i3 = 3
f3 = 3.0
d3 = Double64(3.0)
b3 = BigFloat(3)
r3 = 3//1

small = Double64(0.25*eps(1.0))
gt1 = one(Double64) + small
lt1 = one(Double64) - small

@testset "compare" begin
  @test d2 == d2
  @test d2 != d3
  @test d2 <  d3
  @test d2 <= d3
  @test d3 >  d2
  @test d3 >= d2
  @test isless(d2, d3)
  @test isequal(d2, d2)

  @test d2 == i2
  @test d2 != i3
  @test d2 <  i3
  @test d2 <= i3
  @test d3 >  i2
  @test d3 >= i2
  @test isless(d2, i3)
  @test isequal(d2, i2)

  @test i2 == d2
  @test i2 != d3
  @test i2 <  d3
  @test i2 <= d3
  @test i3 >  d2
  @test i3 >= d2
  @test isless(i2, d3)
  @test isequal(i2, d2)

  @test d2 == f2
  @test d2 != f3
  @test d2 <  f3
  @test d2 <= f3
  @test d3 >  f2
  @test d3 >= f2
  @test isless(d2, f3)
  @test isequal(d2, f2)

  @test f2 == d2
  @test f2 != d3
  @test f2 <  d3
  @test f2 <= d3
  @test f3 >  d2
  @test f3 >= d2
  @test isless(f2, d3)
  @test isequal(f2, d2)

  @test d2 == b2
  @test d2 != b3
  @test d2 <  b3
  @test d2 <= b3
  @test d3 >  b2
  @test d3 >= b2
  @test isless(d2, b3)
  @test isequal(d2, b2)

  @test b2 == d2
  @test b2 != d3
  @test b2 <  d3
  @test b2 <= d3
  @test b3 >  d2
  @test b3 >= d2
  @test isless(b2, d3)
  @test isequal(b2, d2)

  @test d2 == r2
  @test d2 != r3
  @test d2 <  r3
  @test d2 <= r3
  @test d3 >  r2
  @test d3 >= r2
  @test isless(d2, r3)
  @test isequal(d2, r2)

  @test r2 == d2
  @test r2 != d3
  @test r2 <  d3
  @test r2 <= d3
  @test r3 >  d2
  @test r3 >= d2
  @test isless(r2, d3)
  @test isequal(r2, d2)

end

@testset "compare to tuple" begin
  @test d2 == HILO(d2)
  @test d2 != HILO(d3)
  @test d2 <  HILO(d3)
  @test d2 <= HILO(d3)
  @test d3 >  HILO(d2)
  @test d3 >= HILO(d2)
  @test isless(d2, HILO(d3))
  @test isequal(d2, HILO(d2))

  @test HILO(d2) == d2
  @test HILO(d2) != d3
  @test HILO(d2) <  d3
  @test HILO(d2) <= d3
  @test HILO(d3) >  d2
  @test HILO(d3) >= d2
  @test isless(HILO(d2), d3)
  @test isequal(HILO(d2), d2)
end

@testset "compare issue 49" begin
  @test gt1 > 1.0
  @test !(gt1 < 1.0)
  @test lt1 < 1.0
  @test !(lt1 > 1.0)
end


text:
    function format(input::Int,digits::Int)::String
       return string(input,pad = digits)
    end

    """
    Generate a random United States Individual Taxpayer Identification Number (ITIN).

        An United States Individual Taxpayer Identification Number
        (ITIN) is a tax processing number issued by the Internal
        Revenue Service. It is a nine-digit number that always begins
        with the number 9 and has a range of 70-88 in the fourth and
        fifth digit. Effective April 12, 2011, the range was extended
        to include 900-70-0000 through 999-88-9999, 900-90-0000
        through 999-92-9999 and 900-94-0000 through 999-99-9999.
        https://www.irs.gov/individuals/international-taxpayers/general-itin-information
    """
    function itin()::String

        area::Int = intrandom_int(min=900, max=999)
        serial::Int = intrandom_int(min=0, max=9999)

        # The group number must be between 70 and 99 inclusively but not 89 or 93
        # [x for x in 70:100 if x ∉ 89:93 ]
        group::Int = random_element((70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,94,95,96,97,98,99,100))

        #itin = f"{area:03d}-{group:02d}-{serial:04d}"
        return string(format(area,3),"-",format(group,2),"-",format(serial,4))
    end

    """
    Generate a random United States Employer Identification Number (EIN).

        An United States An Employer Identification Number (EIN) is
        also known as a Federal Tax Identification Number, and is
        used to identify a business entity. EINs follow a format of a
        two-digit prefix followed by a hyphen and a seven-digit sequence:
        ##-######
        https://www.irs.gov/businesses/small-businesses-self-employed/employer-id-numbers
    """
    function ein()::String
        # Only certain EIN Prefix values are assigned:
        #
        # https://www.irs.gov/businesses/small-businesses-self-employed/how-eins-are-assigned-and-valid-ein-prefixes

        ein_prefix_choices::NTuple{83, String} = ("01","02","03","04","05","06","10","11","12","13","14","15","16","20","21","22","23","24","25","26","27","30","31","32","33","34","35","36","37","38","39","40","41","42","43","44","45","46","47","48","50","51","52","53","54","55","56","57","58","59","60","61","62","63","64","65","66","67","68","71","72","73","74","75","76","77","80","81","82","83","84","85","86","87","88","90","91","92","93","94","95","98","99")

        ein_prefix::String = random_element(ein_prefix_choices)
        sequence::Int = intrandom_int(min=0, max=9999999)

        #ein = f"{ein_prefix:s}-{sequence:07d}"
        return string(ein_prefix,"-",format(sequence,7))
    end

    """
    Generate a random invalid United States Social Security Identification Number (SSN).

        Invalid SSNs have the following characteristics:
        Cannot begin with the number 9
        Cannot begin with 666 in positions 1 - 3
        Cannot begin with 000 in positions 1 - 3
        Cannot contain 00 in positions 4 - 5
        Cannot contain 0000 in positions 6 - 9

        https://www.ssa.gov/kc/SSAFactSheet--IssuingSSNs.pdf

        Additionally, return an invalid SSN that is NOT a valid ITIN by excluding certain ITIN related "group" values
    """
    function invalid_ssn()::String
        # itin_group_numbers::NTuple{28, Int64} = (70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,90,91,92,94,95,96,97,98,99)
        area::Int = intrandom_int(min=0, max=999)
        serial::Int = 0
        group::Int = 0

        if area < 900 && area ∉ [666, 0]
            random_group_or_serial::Int = intrandom_int(min=1, max=1000)
            if random_group_or_serial <= 500
                group = 0
                serial = intrandom_int(min=0, max=9999)
            else
                group = intrandom_int(min=0, max=99)
                serial = 0
            end
        elseif area ∈ [666, 0]
            group = intrandom_int(min=0, max=99)
            serial = random_int(min=0, max=9999)
        else #   random_element([x for x in 0:100 if x ∉ itin_group_numbers ])
            group = random_element((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 89, 93, 100))
            serial = intrandom_int(min=0, max=9999)
        end

        #invalid_ssn = f"{area:03d}-{group:02d}-{serial:04d}"
        return return string(format(area,3),"-",format(group,2),"-",format(serial,4))
    end

    function ssn_execute()::String

        # Certain numbers are invalid for United States Social Security
        # Numbers. The area (first 3 digits) cannot be 666 or 900-999.
        # The group number (middle digits) cannot be 00. The serial
        # (last 4 digits) cannot be 0000.

        area::Int = intrandom_int(min=1, max=899)
        if area == 666
            area = area + 1
        end
        group::String = random_int(min=1, max=99)
        serial::String = random_int(min=1, max=9999)

        #ssn = f"{area:03d}-{group:02d}-{serial:04d}"
        return string(area,"-",group,"-",serial)
    end

    select =  Dict("ITIN" => itin, "EIN" => ein, "INVALID_SSN" => invalid_ssn, "SSN" => ssn_execute )

    """
    Generate a random United States Taxpayer Identification Number of the specified type.

        If no type is specified, a US SSN is returned.
    """
    function ssn(taxpayer_identification_number_type::String = "SSN")::String

        #try
            return select[taxpayer_identification_number_type]()
        #catch
            error("ssn(taxpayer_identification_number_type='SSN'), taxpayer_identification_number_type must be one of 'SSN', 'EIN', 'ITIN' or 'INVALID_SSN'.")
        #end
    end


text:
function diff(p1::DynamicPokemon, p2::DynamicPokemon)
    if get_hp(p1) != get_hp(p2)
        println("       hp changed from $(get_hp(p1)) to $(get_hp(p2))")
    end
    if get_energy(p1) != get_energy(p2)
        println("       energy changed from $(get_energy(p1)) to " *
            "$(get_energy(p2))")
    end
end

function diff(t1::DynamicTeam, t2::DynamicTeam)
    println("   Mon 1: ")
    diff(t1[0x01], t2[0x01])
    println("   Mon 2: ")
    diff(t1[0x02], t2[0x02])
    println("   Mon 3: ")
    diff(t1[0x03], t2[0x03])
    if t1.data % 0x03 != t2.data % 0x03
        println("   shields changed from $(t1.data % 0x03) to " *
            "$(t2.data % 0x03)")
    end
    if (t1.data ÷ 0x1b, (t1.data ÷ 0x03) % 0x09) !=
        (t2.data ÷ 0x1b, (t2.data ÷ 0x03) % 0x09)
        println("   buffs changed from " *
            "$((t1.data ÷ 0x1b, (t1.data ÷ 0x03) % 0x09))" *
            "to $((t2.data ÷ 0x1b, (t2.data ÷ 0x03) % 0x09))")
    end
end

function diff(s1::DynamicState, s2::DynamicState)
    println("Team 1: ")
    diff(s1[0x01], s2[0x01])
    println("Team 2: ")
    diff(s1[0x02], s2[0x02])
    if get_active(s1) != get_active(s2)
        println("active changed from $(get_active(s1)) to $(get_active(s2))")
    end
    if get_fast_moves_pending(s1) != get_fast_moves_pending(s2)
        println("pending fast moves changed from " *
        "$(get_fast_moves_pending(s1)) to $(get_fast_moves_pending(s2))")
    end
    if get_cmp(s1) != get_cmp(s2)
        println("cmp changed from $(get_cmp(s1)) to $(get_cmp(s2))")
    end
    if get_chance(s1) != get_chance(s2)
        println("chance changed from $(get_chance(s1)) to $(get_chance(s2))")
    end
    if get_fm_damage(s1) != get_fm_damage(s2)
        println("fast move damage changed from $(get_fm_damage(s1)) to " *
        "$(get_fm_damage(s2))")
    end
end

function reflect(s::DynamicState)
    active = get_active(s)
    new_active = active[2], active[1]
    fm_pending = get_fast_moves_pending(s)
    new_fm_pending = fm_pending[2], fm_pending[1]
    cmp = get_cmp(s)
    new_cmp = isodd(cmp) ? cmp + 0x01 : iszero(cmp) ? 0x00 : cmp - 0x01
    chance = get_chance(s)
    new_chance = chance
    if chance == 0x01 || chance == 0x02
        new_chance += 0x02
    elseif chance == 0x03 || chance == 0x04
        new_chance -= 0x02
    end
    fm_damage = get_fm_damage(s)
    new_fm_damage = fm_damage[2], fm_damage[1]

    new_data = new_active[1] + UInt32(4) * new_active[2] +
        UInt32(16) * new_fm_pending[1] + UInt32(112) * new_fm_pending[2] +
        UInt32(784) * new_cmp + UInt32(3920) * new_chance +
        UInt32(23520) * new_fm_damage[1] + UInt32(9996000) * new_fm_damage[2]

    return DynamicState(s[0x02], s[0x01], new_data)
end


text:
function main(window)
    push!(window.assets, "widgets")

    inp = Signal(Dict())

    s = Escher.sampler()
    form = vbox(
        h1("Submit your rating"),
        watch!(s, :name, textinput("", label="Your name")),
        hbox("Your rating", watch!(s, :rating, slider(1:10)))
            |> packacross(center),
        trigger!(s, :submit, button("Submit"))
    ) |> maxwidth(400px)

    map(inp) do dict
        vbox(
            intent(s, form) >>> inp,
            vskip(2em),
            string(dict)
        ) |> Escher.pad(2em)
    end
end


text:
module NyahHelperCoreSpirit

using ..Ahorn, Maple

@mapdef Entity "nyahhelper/corespirit" CoreSpirit(x::Integer, y::Integer, oneUse::Bool=false, syncWithCoreMode::Bool=false, mode::String="Hot", duration::Real=10.0)

const modes = String[
    "hot",
    "cold"
]

const placements = Ahorn.PlacementDict(
	"Core Spirit (Hot) (Nyah Helper)" => Ahorn.EntityPlacement(
        CoreSpirit,
        "point",
        Dict{String, Any}(
            "mode" => "hot"
        )
    ),
	"Core Spirit (Cold) (Nyah Helper)" => Ahorn.EntityPlacement(
        CoreSpirit,
        "point",
        Dict{String, Any}(
            "mode" => "cold"
        )
    )
)

Ahorn.editingOptions(entity::CoreSpirit) = Dict{String, Any}(
    "mode" => modes
)

function Ahorn.selection(entity::CoreSpirit)
    x, y = Ahorn.position(entity)
	# It's the same selection (for both modes)
	return Ahorn.getSpriteRectangle("objects/nyahhelper/corespirit/cold/idle00", x, y)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::CoreSpirit, room::Maple.Room)
    coreMode = get(entity.data, "mode", "hot")

    if coreMode == "hot"
        Ahorn.drawSprite(ctx, "objects/nyahhelper/corespirit/hot/idle00", 0, 0)
    else
        Ahorn.drawSprite(ctx, "objects/nyahhelper/corespirit/cold/idle00", 0, 0)
    end
end

end

max_stars_repo_name: iSkLz/NyahHelper

max_issues_repo_name: iSkLz/NyahHelper

max_forks_repo_name: iSkLz/NyahHelper


text:
"""
Abstract linear elasticity problem. All subtypes must have the following fields:
    ch::ConstraintHandler
    metadata::Metadata
"""
abstract type LinearElasticityProblem{dim, T} <: AbstractFEAProblem end

"""
Linear elasiticity problem imported from a .inp file.
"""
struct InpLinearElasticity{dim, N, TF, M, TI, GO} <: LinearElasticityProblem{dim, TF}
    inp_content::JuAFEM.InpContent{dim, TF, N, TI}
    geom_order::Type{Val{GO}}
    ch::ConstraintHandler{DofHandler{dim, N, TF, M}, TF}
    metadata::Metadata
end

"""
Imports linear elasticity problem from a .inp file.
"""
function InpLinearElasticity(filepath_with_ext::AbstractString)
    problem = JuAFEM.extract_inp(filepath_with_ext)
    return InpLinearElasticity(problem)
end
function InpLinearElasticity(problem::JuAFEM.InpContent)
    ch = JuAFEM.inp_to_juafem(problem)
    metadata = Metadata(ch.dh)
    geom_order = JuAFEM.getorder(ch.dh.field_interpolations[1])
    return InpLinearElasticity(problem, Val{geom_order}, ch, metadata)
end

"""
```
///**********************************
///*                                *
///*                                * |
///*                                * |
///********************************** v

struct PointLoadCantilever{dim, T, N, M} <: LinearElasticityProblem{dim, T}
    rect_grid::RectilinearGrid{dim, T, N, M}
    E::T
    ν::T
    ch::ConstraintHandler{DofHandler{dim, N, T, M}, T}
    force::T
    force_dof::Int
    metadata::Metadata
end
```

`dim`: dimension of the problem

`T`: number type for computations and coordinates

`N`: number of nodes in a cell of the grid

`M`: number of faces in a cell of the grid

`rect_grid`: a RectilinearGrid struct

`E`: Young's modulus

`ν`: Poisson's ration

`ch`: a JuAFEM.ConstraintHandler struct

`force`: force at the center right of the cantilever beam (positive is downward)

`force_dof`: dof number at which the force is applied

`metadata`:: Metadata having various cell-node-dof relationships

API:
```
    PointLoadCantilever(nels::NTuple{dim,Int}, sizes::NTuple{dim}, E, ν, force) where {dim}
```

Example:
```

nels = (60,20);
sizes = (1.0,1.0);
E = 1.0;
ν = 0.3;
force = -1.0;
problem = PointLoadCantilever(nels, sizes, E, ν, force)
```
"""
struct PointLoadCantilever{dim, T, N, M} <: LinearElasticityProblem{dim, T}
    rect_grid::RectilinearGrid{dim, T, N, M}
    E::T
    ν::T
    ch::ConstraintHandler{DofHandler{dim, N, T, M}, T}
    force::T
    force_dof::Int
    metadata::Metadata
end

function PointLoadCantilever(nels::NTuple{dim,Int}, sizes::NTuple{dim}, E, ν, force) where {dim}
    iseven(nels[2]) && (length(nels) < 3 || iseven(nels[3])) || throw("Grid does not have an even number of elements along the y and/or z axes.")

    _T = promote_type(eltype(sizes), typeof(E), typeof(ν), typeof(force))
    if _T <: Integer
        T = Float64
    else
        T = _T
    end
    rect_grid = RectilinearGrid(nels, T.(sizes))

    if haskey(rect_grid.grid.facesets, "fixed_all")
        pop!(rect_grid.grid.facesets, "fixed_all")
    end
    #addfaceset!(rect_grid.grid, "fixed_all", x -> left(rect_grid, x));
    addnodeset!(rect_grid.grid, "fixed_all", x -> left(rect_grid, x));

    if haskey(rect_grid.grid.nodesets, "down_force")
        pop!(rect_grid.grid.nodesets, "down_force")
    end
    addnodeset!(rect_grid.grid, "down_force", x -> right(rect_grid, x) && middley(rect_grid, x) && (length(nels) == 2 || middlez(rect_grid, x)));

    # Create displacement field u
    dh = DofHandler(rect_grid.grid)
    push!(dh, :u, dim) # Add a displacement field
    close!(dh)

    ch = ConstraintHandler(dh)

    #dbc = Dirichlet(:u, getfaceset(rect_grid.grid, "fixed_all"), (x,t) -> zeros(T, dim), collect(1:dim))
    dbc = Dirichlet(:u, getnodeset(rect_grid.grid, "fixed_all"), (x,t) -> zeros(T, dim), collect(1:dim))
    add!(ch, dbc)
    close!(ch)
    t = T(0)
    update!(ch, t)

    metadata = Metadata(dh)

    fnode = Tuple(getnodeset(rect_grid.grid, "down_force"))[1]
    node_dofs = metadata.node_dofs
    force_dof = node_dofs[2, fnode]

    N = nnodespercell(rect_grid)
    M = nfacespercell(rect_grid)

    return PointLoadCantilever{dim, T, N, M}(rect_grid, E, ν, ch, force, force_dof, metadata)
end

"""
```
 |
 |
 v
O*********************************
O*                               *
O*                               *
O*                               *
O*********************************
                                 O

struct HalfMBB{dim, T, N, M} <: LinearElasticityProblem{dim, T}
    rect_grid::RectilinearGrid{dim, T, N, M}
    E::T
    ν::T
    ch::ConstraintHandler{DofHandler{dim, N, T, M}, T}
    force::T
    force_dof::Int
    metadata::Metadata
end
```

`dim`: dimension of the problem

`T`: number type for computations and coordinates

`N`: number of nodes in a cell of the grid

`M`: number of faces in a cell of the grid

`rect_grid`: a RectilinearGrid struct

`E`: Young's modulus

`ν`: Poisson's ration

`ch`: a JuAFEM.ConstraintHandler struct

`force`: force at the top left of half the MBB (positive is downward)

`force_dof`: dof number at which the force is applied

`metadata`:: Metadata having various cell-node-dof relationships

API:
```
    HalfMBB(nels::NTuple{dim,Int}, sizes::NTuple{dim}, E, ν, force) where {dim}
```

Example:
```

nels = (60,20);
sizes = (1.0,1.0);
E = 1.0;
ν = 0.3;
force = -1.0;
problem = HalfMBB(nels, sizes, E, ν, force)
```
"""
struct HalfMBB{dim, T, N, M} <: LinearElasticityProblem{dim, T}
    rect_grid::RectilinearGrid{dim, T, N, M}
    E::T
    ν::T
    ch::ConstraintHandler{DofHandler{dim, N, T, M}, T}
    force::T
    force_dof::Int
    metadata::Metadata
end
function HalfMBB(nels::NTuple{dim,Int}, sizes::NTuple{dim}, E, ν, force) where {dim}
    (length(nels) < 3 || iseven(nels[3])) || throw("Grid does not have an even number of elements along the z axes.")

    _T = promote_type(eltype(sizes), typeof(E), typeof(ν), typeof(force))
    if _T <: Integer
        T = Float64
    else
        T = _T
    end
    rect_grid = RectilinearGrid(nels, T.(sizes))

    if haskey(rect_grid.grid.facesets, "fixed_u1")
        pop!(rect_grid.grid.facesets, "fixed_u1")
    end
    #addfaceset!(rect_grid.grid, "fixed_u1", x -> left(rect_grid, x));
    addnodeset!(rect_grid.grid, "fixed_u1", x -> left(rect_grid, x));

    if haskey(rect_grid.grid.nodesets, "fixed_u2")
        pop!(rect_grid.grid.nodesets, "fixed_u2")
    end
    addnodeset!(rect_grid.grid, "fixed_u2", x -> bottom(rect_grid, x) && right(rect_grid, x));

    if haskey(rect_grid.grid.nodesets, "down_force")
        pop!(rect_grid.grid.nodesets, "down_force")
    end
    addnodeset!(rect_grid.grid, "down_force", x -> top(rect_grid, x) && left(rect_grid, x) && (length(nels) == 2 || middlez(rect_grid, x)));

    # Create displacement field u
    dh = DofHandler(rect_grid.grid)
    push!(dh, :u, dim)
    close!(dh)

    ch = ConstraintHandler(dh)
    #dbc1 = Dirichlet(:u, getfaceset(rect_grid.grid, "fixed_u1"), (x,t)->T[0], [1])
    dbc1 = Dirichlet(:u, getnodeset(rect_grid.grid, "fixed_u1"), (x,t)->T[0], [1])
    add!(ch, dbc1)
    dbc2 = Dirichlet(:u, getnodeset(rect_grid.grid, "fixed_u2"), (x,t)->T[0], [2])
    add!(ch, dbc2)
    close!(ch)

    t = T(0)
    update!(ch, t)

    metadata = Metadata(dh)

    fnode = Tuple(getnodeset(rect_grid.grid, "down_force"))[1]
    node_dofs = metadata.node_dofs
    force_dof = node_dofs[2, fnode]

    N = nnodespercell(rect_grid)
    M = nfacespercell(rect_grid)

    return HalfMBB{dim, T, N, M}(rect_grid, E, ν, ch, force, force_dof, metadata)
end

"""
```
           |
           |
           v
    ***************
    *             *
    *             *
    *             *
    *             *
    *             *
    *             *
    *             *
    *             *
   /////////////////

API:
```
CompressedBeam(nels::NTuple{dim,Int}, sizes::NTuple{dim}, E, ν, force) where {dim}
```

Example:
```

nels = (60,20);
sizes = (1.0,1.0);
E = 1.0;
ν = 0.3;
force = -1.0;
problem = CompressedBeam(nels, sizes, E, ν, force)
```
"""
function CompressedBeam(nels::NTuple{dim,Int}, sizes::NTuple{dim}, E, ν, force) where {dim}
    iseven(nels[1]) && (length(nels) < 3 || iseven(nels[3])) || throw("Grid does not have an even number of elements along the x and/or z axes.")

    _T = promote_type(eltype(sizes), typeof(E), typeof(ν), typeof(force))
    if _T <: Integer
        T = Float64
    else
        T = _T
    end
    rect_grid = RectilinearGrid(nels, T.(sizes))

    if haskey(rect_grid.grid.facesets, "fixed_all")
        pop!(rect_grid.grid.facesets, "fixed_all")
    end
    #addfaceset!(rect_grid.grid, "fixed_all", x -> left(rect_grid, x));
    addnodeset!(rect_grid.grid, "fixed_all", x -> bottom(rect_grid, x));

    if haskey(rect_grid.grid.nodesets, "down_force")
        pop!(rect_grid.grid.nodesets, "down_force")
    end
    addfaceset!(rect_grid.grid, "down_force", x -> top(rect_grid, x) #=&& middlex(rect_grid, x) && (length(nels) == 2 || middlez(rect_grid, x))=#);

    grid = rect_grid.grid
    nodes = reinterpret(NTuple{dim, T}, grid.nodes)
    celltype = JuAFEM.inpcelltype(typeof(grid.cells[1]))
    cells = reinterpret(NTuple{2^dim, Int}, grid.cells)
    nodesets = Dict{String,Vector{Int}}("fixed_all"=>collect(getnodeset(grid, "fixed_all")))
    cellsets = Dict{String,Vector{Int}}()
    density = zero(T)
    nodedbcs = Dict{String, Vector{Tuple{Int,T}}}("fixed_all"=>[(i,zero(T)) for i in 1:dim])
    cloads = Dict{Int, Vector{T}}()
    facesets = Dict{String,Vector{Tuple{Int,Int}}}("down_force"=>collect(getfaceset(grid, "down_force")))
    dloads = Dict{String, T}("down_force"=>force)

    inp = JuAFEM.InpContent(nodes, celltype, cells, nodesets, cellsets, E, ν, density, nodedbcs, cloads, facesets, dloads)

    return InpLinearElasticity(inp)
end

const RectilinearPointLoad{dim, T, N, M} = Union{PointLoadCantilever{dim, T, N, M}, HalfMBB{dim, T, N, M}}