text:
"""
Extend this type for different `Potential` choices
"""
abstract type Potential end

"""
Extend this type for different `InitialData` choices
"""
abstract type InitialData end

"""
Extend this type for different `EvolutionEquations`.
"""
abstract type EvolutionEquations end

Base.@kwdef struct AffineNull{TP<:Potential} <: EvolutionEquations
    potential      :: TP  = ConstPotential()
end

"""
Parameters for the time evolution
"""
Base.@kwdef struct Integration{T,Tdt,S}
    dt              :: Tdt  = :auto
    tmax            :: T
    ODE_method      :: S    = AB3()
    adaptive        :: Bool = false
    # relative tolerance for adaptive integrators
    reltol          :: Float64 = 1e-6
    # filter_poststep :: Bool = true
end

"""
Parameters for Input/Output
"""
Base.@kwdef struct InOut
    # negative values suppress output
    out_bulk_every              :: Int  = -1

    out_bulk_every_t            :: Float64  = -1.0

    checkpoint_every_walltime_hours :: Float64 = -1.0

    # stop and checkpoint upon reaching this walltime
    max_walltime       :: Float64 = 1.e20

    # trigger termination
    termination_from_file :: Bool    = true
    check_file_every      :: Int     = 10
    termination_file      :: String  = "TERMINATE"

    # name of script
    _parfile           :: String  = splitext(basename(Base.source_path()))[1]

    # use name of script by default
    out_dir            :: String  = _parfile
    checkpoint_dir     :: String  = _parfile

    recover            :: Symbol  = :auto
    recover_dir        :: String  = _parfile

    # be very careful with this option! it will remove the whole folder contents
    # if set to true! use only for fast debugging runs
    remove_existing    :: Bool  = false
end

struct BulkEvolved{T,A,S} <: FlattenedVector{T,3,A}
    x :: S
end
"""
    BulkEvolved{T}(undef, Nu, Nx, Ny)

Construct a container of uninitialized Arrays to hold the bulk variable
evolved in time: phi
"""
function BulkEvolved{T}(::UndefInitializer, Nu::Int, Nx::Int, Ny::Int) where {T<:Real}
    phi = Array{T}(undef, Nu, Nx, Ny)
    x   = (phi=phi,)
    BulkEvolved{T,eltype(x),typeof(x)}(x)
end

getphi(ff::BulkEvolved) = ff.x.phi

function Base.similar(ff::BulkEvolved{T}) where {T}
    phi = similar(getphi(ff))
    x   = (phi=phi,)
    BulkEvolved{T,eltype(x),typeof(x)}(x)
end

struct BulkConstrained{T,A,S} <: FlattenedVector{T,3,A}
    x :: S
end
"""
    BulkConstrained{T}(undef, Nu, Nx, Ny)

Construct a container of uninitialized Arrays to hold all the bulk variables
that are constrained (not evolved in time): S, phid, Sd, A
"""
function BulkConstrained{T}(::UndefInitializer, Nu::Int, Nx::Int, Ny::Int) where {T<:Real}
    phid = Array{T}(undef, Nu, Nx, Ny)
    Sd   = Array{T}(undef, Nu, Nx, Ny)
    A    = Array{T}(undef, Nu, Nx, Ny)
    x    = (phid=phid, Sd=Sd, A=A)
    BulkConstrained{T,eltype(x),typeof(x)}(x)
end

getphid(ff::BulkConstrained) = ff.x.phid
getSd(ff::BulkConstrained)   = ff.x.Sd
getA(ff::BulkConstrained)    = ff.x.A

struct Boundary{T,A,S} <: FlattenedVector{T,3,A}
    x :: S
end
"""
    Boundary{T}(undef, Nx, Ny)

Construct a container of uninitialized Arrays to hold the boundary
variables a4.

These variables are automatically defined on a `(1,Nx,Ny)` grid, rather than
a `(Nx,Ny)` one.
"""
function Boundary{T}(::UndefInitializer, Nx::Int, Ny::Int) where {T<:Real}
    a4  = Array{T}(undef, 1, Nx, Ny)
    x   = (a4=a4,)
    Boundary{T,eltype(x),typeof(x)}(x)
end

geta4(ff::Boundary) = ff.x.a4

struct BulkDeriv{T}
    Du_phi  :: Array{T,3}
    Du_phid :: Array{T,3}
end

function BulkDeriv{T}(::UndefInitializer, Nu::Int, Nx::Int, Ny::Int) where {T<:Real}
    Du_phi   = Array{T}(undef, Nu, Nx, Ny)
    Du_phid  = Array{T}(undef, Nu, Nx, Ny)
    BulkDeriv{T}(Du_phi, Du_phid)
end

"""
    BulkPartition{N,A} <: AbstractPartition{N,A}

Container to store (bulk) quantities that may be spread across different grid
partitions. This is to be thought as a `Tuple` (or an `Array`) of `Bulk` objects.
"""
struct BulkPartition{N,A} <: AbstractPartition{N,A}
    x :: NTuple{N,A}
end
BulkPartition(x...) = BulkPartition(tuple(x...))

function BulkEvolved(bulks::BulkPartition{N}) where{N}
    f = ntuple(i -> BulkEvolved(bulks[i]), N)
    BulkPartition(f)
end

function BulkConstrained(bulks::BulkPartition{N}) where{N}
    f = ntuple(i -> BulkConstrained(bulks[i]), N)
    BulkPartition(f)
end

struct EvolVars{T,N,A} <: FlattenedVector{T,N,A}
    x :: NTuple{N,A}
end
EvolVars(x::NTuple{N,A}) where {N,A} = EvolVars{eltype(A),N,A}(x)

Base.similar(ff::EvolVars{T,N,S}) where {T,N,S} = EvolVars{T,N,S}(similar.(ff.x))

"""
    EvolVars(boundary::Boundary, gauge::Gauge, bulkevols::NTuple)

Build a container to store all the evolved quantities as elements of an
`NTuple`. The idea is to treat them as a single column vector for the point of
view of the time evolution routine.
"""
function EvolVars(bulkevols::BulkPartition{Nsys}) where {Nsys}
    EvolVars(bulkevols.x)
end

@inline getbulkevolvedpartition(ff::EvolVars) = BulkPartition(ff.x)


text:
# Types

const D2 = NTuple{2, Cdouble}

export Float128
struct Float128 <: AbstractFloat
    d2::D2
end

# Functions

function Base.tryparse(::Type{Float128}, str::String)
    r = tryparse(BigFloat, str)
    r === nothing && return r
    Float128(r)
end

export @f128_str
macro f128_str(str::String)
    :(Float128(BigFloat($str)))
end

# function Base.show(io::IO, x::Float128)
#     len = 40
#     buf = Vector{Cchar}(undef, len)
#     prec = 32
#     ccall((:c_dd_swrite, libqd), Cvoid, (Ref{D2}, Cint, Ptr{Cchar}, Cint),
#           x.d2, prec, buf, len)
#     print(io, unsafe_string(pointer(buf)))
# end

Base.show(io::IO, x::Float128) = show(io, BigFloat(x))

# function Float128(x::Float64)
#     r = Ref{D2}()
#     ccall((:c_dd_copy_d, libqd), Cvoid, (Cdouble, Ref{D2}), x, r)
#     Float128(r[])
# end
Float128(x::Float64) = Float128((x, 0.0))

function Float128(x::Union{Bool, Int8, Int16, Int32, UInt8, UInt16, UInt32,
                           Float16, Float32})
    Float128(Float64(x))
end

function Float128(x::Union{Int64, UInt64})
    hi = x >> 32 << 32
    lo = x - hi
    Float128(Float64(hi)) + Float128(Float64(lo))
end

function Float128(x::Rational)
    Float128(numerator(x)) / Float128(denominator(x))
end

function Float128(x::BigFloat)
    r1 = Float128(Float64(x))
    x -= BigFloat(r1)
    r2 = Float128(Float64(x))
    # x -= BigFloat(r2)
    r1 + r2
end

Base.Float64(x::Float128) = x.d2[1]
Base.BigFloat(x::Float128) = big(x.d2[1]) + big(x.d2[2])

Base.big(::Type{Float128}) = Float256

Base.promote_rule(::Type{Float128}, ::Type{Float16}) = Float128
Base.promote_rule(::Type{Float128}, ::Type{Float32}) = Float128
Base.promote_rule(::Type{Float128}, ::Type{Float64}) = Float128
Base.promote_rule(::Type{Float128}, ::Type{BigFloat}) = BigFloat
Base.promote_rule(::Type{Float128}, ::Type{I}) where {I <: Integer} = Float128
Base.promote_rule(::Type{Float128}, ::Type{BigInt}) = BigFloat

# Base.zero(::Type{Float128}) = c128_0
# Base.one(::Type{Float128}) = c128_1
# Base.eps(::Type{Float128}) = c128_eps
Base.zero(::Type{Float128}) = Float128(0.0)
Base.one(::Type{Float128}) = Float128(1.0)
Base.eps(::Type{Float128}) = Float128(2^-104)

# function Float128(::Irrational{:π})
#     r = Ref{D2}()
#     ccall((:c_dd_pi, libqd), Cvoid, (Ref{D2},), r)
#     Float128(r[])
# end
Float128(::Irrational{:ℯ}) = c128_e
Float128(::Irrational{:π}) = c128_pi

Base. +(x::Float128) = x

# function Base. -(x::Float128)
#     r = Ref{D2}()
#     ccall((:c_dd_neg, libqd), Cvoid, (Ref{D2}, Ref{D2}), x.d2, r)
#     Float128(r[])
# end

function Base. -(a::Float128)
    Float128((- a.d2[1], - a.d2[2]))
end

# function Base.abs(x::Float128)
#     r = Ref{D2}()
#     ccall((:c_dd_abs, libqd), Cvoid, (Ref{D2}, Ref{D2}), x.d2, r)
#     Float128(r[])
# end

Base.abs(a::Float128) = a.d2[1] < 0.0 ? -a : a

Base.inv(x::Float128) = one(Float64) / x

# function Base. +(x::Float128, y::Float128)
#     r = Ref{D2}()
#     ccall((:c_dd_add, libqd), Cvoid, (Ref{D2}, Ref{D2}, Ref{D2}), x.d2, y.d2, r)
#     Float128(r[])
# end
#
# function Base. +(x::Float64, y::Float128)
#     r = Ref{D2}()
#     ccall((:c_dd_add_d_dd, libqd), Cvoid, (Cdouble, Ref{D2}, Ref{D2}),
#           x, y.d2, r)
#     Float128(r[])
# end
#
# function Base. +(x::Float128, y::Float64)
#     r = Ref{D2}()
#     ccall((:c_dd_add_dd_d, libqd), Cvoid, (Ref{D2}, Cdouble, Ref{D2}),
#           x.d2, y, r)
#     Float128(r[])
# end

function Base. +(a::Float128, b::Float128)
    # This is the less accurate version ... obeys Cray-style error bound.
    s, e = two_sum(a.d2[1], b.d2[1])
    e += a.d2[2] + b.d2[2]
    s, e = quick_two_sum(s, e)
    Float128((s, e))
end

function Base. +(a::Float128, b::Float64)
    s1, s2 = two_sum(a.d2[1], b)
    s2 += a.d2[2]
    s1, s2 = quick_two_sum(s1, s2)
    Float128((s1, s2))
end

Base. +(a::Float64, b::Float128) = b + a

export wideadd
function wideadd(a::Float64, b::Float64)
    s, e = two_sum(a, b)
    Float128((s, e))
end

# function Base. -(x::Float128, y::Float128)
#     r = Ref{D2}()
#     ccall((:c_dd_sub, libqd), Cvoid, (Ref{D2}, Ref{D2}, Ref{D2}), x.d2, y.d2, r)
#     Float128(r[])
# end
#
# function Base. -(x::Float64, y::Float128)
#     r = Ref{D2}()
#     ccall((:c_dd_sub_d_dd, libqd), Cvoid, (Cdouble, Ref{D2}, Ref{D2}),
#           x, y.d2, r)
#     Float128(r[])
# end
#
# function Base. -(x::Float128, y::Float64)
#     r = Ref{D2}()
#     ccall((:c_dd_sub_dd_d, libqd), Cvoid, (Ref{D2}, Cdouble, Ref{D2}),
#           x.d2, y, r)
#     Float128(r[])
# end

function Base. -(a::Float128, b::Float128)
    # This is the less accurate version ... obeys Cray-style error bound.
    s, e = two_diff(a.d2[1], b.d2[1])
    e += a.d2[2]
    e -= b.d2[2]
    s, e = quick_two_sum(s, e)
    Float128((s, e))
end

function Base. -(a::Float128, b::Float64)
    s1, s2 = two_diff(a.d2[1], b)
    s2 += a.d2[2]
    s1, s2 = quick_two_sum(s1, s2)
    Float128((s1, s2))
end

function Base. -(a::Float64, b::Float128)
  s1, s2 = two_diff(a, b.d2[1])
    s2 -= b.d2[2]
    s1, s2 = quick_two_sum(s1, s2)
    Float128((s1, s2))
end

export widesub
function widesub(a::Float64, b::Float64)
    s, e = two_diff(a, b)
    Float128((s, e))
end

# function Base. *(x::Float128, y::Float128)
#     r = Ref{D2}()
#     ccall((:c_dd_mul, libqd), Cvoid, (Ref{D2}, Ref{D2}, Ref{D2}), x.d2, y.d2, r)
#     Float128(r[])
# end
#
# function Base. *(x::Float64, y::Float128)
#     r = Ref{D2}()
#     ccall((:c_dd_mul_d_dd, libqd), Cvoid, (Cdouble, Ref{D2}, Ref{D2}),
#           x, y.d2, r)
#     Float128(r[])
# end
#
# function Base. *(x::Float128, y::Float64)
#     r = Ref{D2}()
#     ccall((:c_dd_mul_dd_d, libqd), Cvoid, (Ref{D2}, Cdouble, Ref{D2}),
#           x.d2, y, r)
#     Float128(r[])
# end

function Base. *(a::Float128, b::Float128)
    p1, p2 = two_prod(a.d2[1], b.d2[1])
    p2 += (a.d2[1] * b.d2[2] + a.d2[2] * b.d2[1])
    p1, p2 = quick_two_sum(p1, p2)
    Float128((p1, p2))
end

function Base. *(a::Float128, b::Float64)
    p1, p2 = two_prod(a.d2[1], b)
    p2 += a.d2[2] * b
    p1, p2 = quick_two_sum(p1, p2)
    Float128((p1, p2))
end

Base. *(a::Float64, b::Float128) = b * a
function Base.widemul(a::Float64, b::Float64)
    p, e = two_prod(a, b)
    Float128((p, e))
end

# function Base. /(x::Float128, y::Float128)
#     r = Ref{D2}()
#     ccall((:c_dd_div, libqd), Cvoid, (Ref{D2}, Ref{D2}, Ref{D2}), x.d2, y.d2, r)
#     Float128(r[])
# end
#
# function Base. /(x::Float64, y::Float128)
#     r = Ref{D2}()
#     ccall((:c_dd_div_d_dd, libqd), Cvoid, (Cdouble, Ref{D2}, Ref{D2}),
#           x, y.d2, r)
#     Float128(r[])
# end
#
# function Base. /(x::Float128, y::Float64)
#     r = Ref{D2}()
#     ccall((:c_dd_div_dd_d, libqd), Cvoid, (Ref{D2}, Cdouble, Ref{D2}),
#           x.d2, y, r)
#     Float128(r[])
# end

function Base. /(a::Float128, b::Float128)
    # sloppy div
    q1 = a.d2[1] / b.d2[1]      # approximate quotient

    # compute  this - q1 * dd
    r = b * q1
    s1, s2 = two_diff(a.d2[1], r.d2[1])
    s2 -= r.d2[2]
    s2 += a.d2[2]

    # get next approximation
    q2 = (s1 + s2) / b.d2[1]

    # renormalize
    r_hi, r_lo = quick_two_sum(q1, q2)
    Float128((r_hi, r_lo))
end

function Base. /(a::Float128, b::Float64)
  q1 = a.d2[1] / b              # approximate quotient.

    # Compute  this - q1 * d
    p1, p2 = two_prod(q1, b)
    s, e = two_diff(a.d2[1], p1)
    e += a.d2[2]
    e -= p2

    # get next approximation.
    q2 = (s + e) / b

    # renormalize
    r_hi, r_lo = quick_two_sum(q1, q2)
    Float128((r_hi, r_lo))
end

Base. /(a::Float64, b::Float128) = Float128(a) / b

export widediv
function widediv(a::Float64, b::Float64)
    q1 = a / b

    # Compute  a - q1 * b
    p1, p2 = two_prod(q1, b)
    s, e = two_diff(a, p1)
    e -= p2

    # get next approximation
    q2 = (s + e) / b

    s, e = quick_two_sum(q1, q2)

    Float128((s, e))
end

Base. \(x::Float128, y::Float128) = y / x
Base. \(x::Float64, y::Float128) = y / x
Base. \(x::Float128, y::Float64) = y / x

function Base.cmp(x::Float128, y::Float128)
    r = Ref{Cint}()
    ccall((:c_dd_comp, libqd), Cvoid, (Ref{D2}, Ref{D2}, Ref{Cint}),
          x.d2, y.d2, r)
    Int(r[])
end

function Base.cmp(x::Float64, y::Float128)
    r = Ref{Cint}()
    ccall((:c_dd_comp_d_dd, libqd), Cvoid, (Cdouble, Ref{D2}, Ref{Cint}),
          x, y.d2, r)
    Int(r[])
end

function Base.cmp(x::Float128, y::Float64)
    r = Ref{Cint}()
    ccall((:c_dd_comp_dd_d, libqd), Cvoid, (Ref{D2}, Cdouble, Ref{Cint}),
          x.d2, y, r)
    Int(r[])
end

# Base. ==(x::Float128, y::Float128) = cmp(x, y) == 0
# Base. ==(x::Float64, y::Float128) = cmp(x, y) == 0
# Base. ==(x::Float128, y::Float64) = cmp(x, y) == 0
#
# Base. !=(x::Float128, y::Float128) = cmp(x, y) != 0
# Base. !=(x::Float64, y::Float128) = cmp(x, y) != 0
# Base. !=(x::Float128, y::Float64) = cmp(x, y) != 0
#
# Base. <(x::Float128, y::Float128) = cmp(x, y) < 0
# Base. <(x::Float64, y::Float128) = cmp(x, y) < 0
# Base. <(x::Float128, y::Float64) = cmp(x, y) < 0
#
# Base. <=(x::Float128, y::Float128) = cmp(x, y) <= 0
# Base. <=(x::Float64, y::Float128) = cmp(x, y) <= 0
# Base. <=(x::Float128, y::Float64) = cmp(x, y) <= 0
#
# Base. >(x::Float128, y::Float128) = cmp(x, y) > 0
# Base. >(x::Float64, y::Float128) = cmp(x, y) > 0
# Base. >(x::Float128, y::Float64) = cmp(x, y) > 0
#
# Base. >=(x::Float128, y::Float128) = cmp(x, y) >= 0
# Base. >=(x::Float64, y::Float128) = cmp(x, y) >= 0
# Base. >=(x::Float128, y::Float64) = cmp(x, y) >= 0

Base. ==(a::Float128, b::Float128) = a.d2[1] == b.d2[1] && a.d2[2] == b.d2[2]
Base. ==(a::Float64, b::Float128) = a == b.d2[1] && b.d2[2] == 0.0
Base. ==(a::Float128, b::Float64) = a.d2[1] == b && a.d2[2] == 0.0

Base. !=(a::Float128, b::Float128) = a.d2[1] != b.d2[1] || a.d2[2] != b.d2[2]
Base. !=(a::Float64, b::Float128) = a != b.d2[1] || b.d2[2] != 0.0
Base. !=(a::Float128, b::Float64) = a.d2[1] != b || a.d2[2] != 0.0

Base. <(a::Float128, b::Float128) =
    a.d2[1] < b.d2[1] || (a.d2[1] == b.d2[1] && a.d2[2] < b.d2[2])
Base. <(a::Float64, b::Float128) =
    a < b.d2[1] || (a == b.d2[1] && b.d2[2] > 0.0)
Base. <(a::Float128, b::Float64) =
    a.d2[1] < b || (a.d2[1] == b && a.d2[2] < 0.0)

Base. <=(a::Float128, b::Float128) =
    a.d2[1] < b.d2[1] || (a.d2[1] == b.d2[1] && a.d2[2] <= b.d2[2])
Base. <=(a::Float64, b::Float128) = b >= a
Base. <=(a::Float128, b::Float64) =
    a.d2[1] < b || (a.d2[1] == b && a.d2[2] <= 0.0)

Base. >(a::Float128, b::Float128) =
    a.d2[1] > b.d2[1] || (a.d2[1] == b.d2[1] && a.d2[2] > b.d2[2])
Base. >(a::Float64, b::Float128) =
    a > b.d2[1] || (a == b.d2[1] && b.d2[2] < 0.0)
Base. >(a::Float128, b::Float64) =
    a.d2[1] > b || (a.d2[1] == b && a.d2[2] > 0.0)

Base. >=(a::Float128, b::Float128) =
    a.d2[1] > b.d2[1] || (a.d2[1] == b.d2[1] && a.d2[2] >= b.d2[2])
Base. >=(a::Float64, b::Float128) = b <= a
Base. >=(a::Float128, b::Float64) =
    a.d2[1] > b || (a.d2[1] == b && a.d2[2] >= 0.0)

function Base.ceil(x::Float128)
    r = Ref{D2}()
    ccall((:c_dd_ceil, libqd), Cvoid, (Ref{D2}, Ref{D2}), x.d2, r)
    Float128(r[])
end

function Base.floor(x::Float128)
    r = Ref{D2}()
    ccall((:c_dd_floor, libqd), Cvoid, (Ref{D2}, Ref{D2}), x.d2, r)
    Float128(r[])
end

# Note: This breaks ties the wrong way
function Base.round(x::Float128)
    r = Ref{D2}()
    ccall((:c_dd_nint, libqd), Cvoid, (Ref{D2}, Ref{D2}), x.d2, r)
    Float128(r[])
end

function Base.trunc(x::Float128)
    r = Ref{D2}()
    ccall((:c_dd_aint, libqd), Cvoid, (Ref{D2}, Ref{D2}), x.d2, r)
    Float128(r[])
end

function Base. ^(x::Float128, i::Integer)
    r = Ref{D2}()
    ccall((:c_dd_npwr, libqd), Cvoid, (Ref{D2}, Cint, Ref{D2}), x.d2, i, r)
    Float128(r[])
end

# function Base.literal_pow(::typeof(^), x::Float128, ::Val{2})
#     r = Ref{D2}()
#     ccall((:c_dd_sqr, libqd), Cvoid, (Ref{D2}, Ref{D2}), x.d2, r)
#     Float128(r[])
# end

function Base.literal_pow(::typeof(^), a::Float128, ::Val{2})
    p1, p2 = two_sqr(a.d2[1])
    p2 += 2.0 * a.d2[1] * a.d2[2]
    p2 += a.d2[2] * a.d2[2]
    s1, s2 = quick_two_sum(p1, p2)
    Float128((s1, s2))
end

function Base.sqrt(x::Float128)
    r = Ref{D2}()
    ccall((:c_dd_sqrt, libqd), Cvoid, (Ref{D2}, Ref{D2}), x.d2, r)
    Float128(r[])
end

export nroot
function nroot(x::Float128, i::Integer)
    r = Ref{D2}()
    ccall((:c_dd_nroot, libqd), Cvoid, (Ref{D2}, Cint, Ref{D2}), x.d2, i, r)
    Float128(r[])
end

function Base.exp(x::Float128)
    r = Ref{D2}()
    ccall((:c_dd_exp, libqd), Cvoid, (Ref{D2}, Ref{D2}), x.d2, r)
    Float128(r[])
end

Base.exp10(x::Float128) = exp(c128_log10 * x)

Base.exp2(x::Float128) = exp(c128_log2 * x)

function Base.log(x::Float128)
    r = Ref{D2}()
    ccall((:c_dd_log, libqd), Cvoid, (Ref{D2}, Ref{D2}), x.d2, r)
    Float128(r[])
end

function Base.log10(x::Float128)
    r = Ref{D2}()
    ccall((:c_dd_log10, libqd), Cvoid, (Ref{D2}, Ref{D2}), x.d2, r)
    Float128(r[])
end

Base.log2(x::Float128) = c128_1_log2 * log(x)

function Base.sin(x::Float128)
    r = Ref{D2}()
    ccall((:c_dd_sin, libqd), Cvoid, (Ref{D2}, Ref{D2}), x.d2, r)
    Float128(r[])
end

function Base.cos(x::Float128)
    r = Ref{D2}()
    ccall((:c_dd_cos, libqd), Cvoid, (Ref{D2}, Ref{D2}), x.d2, r)
    Float128(r[])
end

function Base.tan(x::Float128)
    r = Ref{D2}()
    ccall((:c_dd_tan, libqd), Cvoid, (Ref{D2}, Ref{D2}), x.d2, r)
    Float128(r[])
end

function Base.asin(x::Float128)
    r = Ref{D2}()
    ccall((:c_dd_asin, libqd), Cvoid, (Ref{D2}, Ref{D2}), x.d2, r)
    Float128(r[])
end

function Base.acos(x::Float128)
    r = Ref{D2}()
    ccall((:c_dd_acos, libqd), Cvoid, (Ref{D2}, Ref{D2}), x.d2, r)
    Float128(r[])
end

function Base.atan(x::Float128)
    r = Ref{D2}()
    ccall((:c_dd_atan, libqd), Cvoid, (Ref{D2}, Ref{D2}), x.d2, r)
    Float128(r[])
end

function Base.atan(x::Float128, y::Float128)
    r = Ref{D2}()
    ccall((:c_dd_atan2, libqd), Cvoid, (Ref{D2}, Ref{D2}, Ref{D2}),
          x.d2, y.d2, r)
    Float128(r[])
end

Base.hypot(x::Float128, y::Float128) = sqrt(x^2 + y^2)

function Base.sinh(x::Float128)
    r = Ref{D2}()
    ccall((:c_dd_sinh, libqd), Cvoid, (Ref{D2}, Ref{D2}), x.d2, r)
    Float128(r[])
end

function Base.cosh(x::Float128)
    r = Ref{D2}()
    ccall((:c_dd_cosh, libqd), Cvoid, (Ref{D2}, Ref{D2}), x.d2, r)
    Float128(r[])
end

function Base.tanh(x::Float128)
    r = Ref{D2}()
    ccall((:c_dd_tanh, libqd), Cvoid, (Ref{D2}, Ref{D2}), x.d2, r)
    Float128(r[])
end

function Base.asinh(x::Float128)
    r = Ref{D2}()
    ccall((:c_dd_asinh, libqd), Cvoid, (Ref{D2}, Ref{D2}), x.d2, r)
    Float128(r[])
end

function Base.acosh(x::Float128)
    r = Ref{D2}()
    ccall((:c_dd_acosh, libqd), Cvoid, (Ref{D2}, Ref{D2}), x.d2, r)
    Float128(r[])
end

function Base.atanh(x::Float128)
    r = Ref{D2}()
    ccall((:c_dd_atanh, libqd), Cvoid, (Ref{D2}, Ref{D2}), x.d2, r)
    Float128(r[])
end

function Base.sincos(x::Float128)
    r = Ref{D2}()
    s = Ref{D2}()
    ccall((:c_dd_sincos, libqd), Cvoid, (Ref{D2}, Ref{D2}, Ref{D2}), x.d2, r, s)
    Float128(r[]), Float128(s[])
end

export sincosh
function sincosh(x::Float128)
    r = Ref{D2}()
    s = Ref{D2}()
    ccall((:c_dd_sincosh, libqd), Cvoid, (Ref{D2}, Ref{D2}, Ref{D2}),
          x.d2, r, s)
    Float128(r[]), Float128(s[])
end

function Base.rand(::Type{Float128})
    r = Ref{D2}()
    ccall((:c_dd_rand, libqd), Cvoid, (Ref{D2},), r)
    Float128(r[])
end

function Base.rand(::Type{Float128}, dims::Dims)
    r = Array{Float128}(undef, dims)
    for i in LinearIndices(r)
        r[i] = rand(Float128)
    end
    r
end

Base.rand(::Type{Float128}, dims::Integer...) = rand(Float128, Dims(dims))

# Constants

# const c128_0 = Float128(0.0)
# const c128_1 = Float128(1.0)
# const c128_eps = Float128(2^-104)

const c128_pi = begin
    r = Ref{D2}()
    ccall((:c_dd_pi, libqd), Cvoid, (Ref{D2},), r)
    Float128(r[])
end
const c128_e = exp(Float128(1.0))
const c128_log10 = log(Float128(10.0))
const c128_log2 = log(Float128(2.0))
const c128_1_log2 = inv(c128_log2)

ext: jl


text:
#this file constructs the equipment lookup database for system powers, lengths
#that is stored in the ocean layout
################################################################################
################################# Cables #######################################
################################################################################
#=kbles=database["cables"]
km_mva=km_mva_set[16]
kbl66=kbles["66.0"][string(km_mva[2])][1]=#
#**
function set_cables_per_km_cost(kbles,km_mva_set)
    eps=10^-2
    ks=get_Cost_Data()
    for km_mva in km_mva_set
        for kbl66 in kbles["66.0"][string(km_mva[2])]
            c1=kbl66.costs.ttl
            c0=deepcopy(mvac_cable(kbl66.mva,kbl66.length-1,kbl66.wnd,kbles["66.0"][string(km_mva[2])],ks)).costs.ttl
            kbl66.costs.perkm_ttl=deepcopy(c1-c0)
            kbl66.mx_rng=kbl66.length+eps
        end
        for kbl220 in kbles["220.0"][string(km_mva[2])]
            c1=kbl220.costs.ttl
            c0=deepcopy(hvac_cable(kbl220.mva,kbl220.length-1,kbl220.wnd,kbles["220.0"][string(km_mva[2])],ks)).costs.ttl
            kbl220.costs.perkm_ttl=deepcopy(c1-c0)
            kbl220.mx_rng=kbl220.length+eps
        end
        for kbl400 in kbles["400.0"][string(km_mva[2])]
            c1=kbl400.costs.ttl
            c0=deepcopy(hvac_cable(kbl400.mva,kbl400.length-1,kbl400.wnd,kbles["400.0"][string(km_mva[2])],ks)).costs.ttl
            kbl400.costs.perkm_ttl=deepcopy(c1-c0)
            kbl400.mx_rng=kbl400.length+eps
        end
        for kbl300 in kbles["300.0"][string(km_mva[2])]
            c1=kbl300.costs.ttl
            c0=deepcopy(hvdc_cable(kbl300.mva,kbl300.length-1,kbl300.wnd,kbles["300.0"][string(km_mva[2])],ks)).costs.ttl
            kbl300.costs.perkm_ttl=deepcopy(c1-c0)
            kbl300.mx_rng=kbl300.length+eps
        end
    end
    return kbles
end

#**
function get_cable_table(km_mva_set,wnd)
    ks=get_Cost_Data()
    cable_tables=Array{Array{Array{cable,1},1},1}()
    #HVDC
    cbl_datas=[get_300kV_cables()]
    for cbl_data in cbl_datas
        cable_table=Array{Array{cable,1},1}()
        for km_mva in km_mva_set
            cbl=cable()
            cbl.mva=km_mva[2]
            cbl.wnd=wnd
            push!(cable_table,[optimal_hvdc_cable(cbl,cbl_data,km_mva[1],ks)])
        end
        push!(cable_tables,deepcopy(cable_table))
    end
    #HVAC
    cbl_datas=[get_220kV_cables(),get_400kV_cables()]
    for cbl_data in cbl_datas
        cable_table=Array{Array{cable,1},1}()
        for km_mva in km_mva_set
            cable_interval=Tuple
            cable_interval=(1,km_mva[1],km_mva[2])
            push!(cable_table,solve_cable_interval(cable_interval,cbl_data,wnd,ks))
        end
        push!(cable_tables,deepcopy(cable_table))
    end
    #MVAC
    cbl_data=get_66kV_cables()
    cable_table=Array{Array{cable,1},1}()
    for km_mva in km_mva_set
        cable_interval=Tuple
        cable_interval=(1,km_mva[1],km_mva[2])
        push!(cable_table,solve_cable_interval(cable_interval,cbl_data,wnd,ks))
    end
    push!(cable_tables,deepcopy(cable_table))
    cable_dictionary=make_dictionaries_cables(cable_tables)
    return cable_dictionary
end

#finds transistions in cable sizes per up to km at mva for kv
#**
function solve_cable_interval(cable_interval,cbl_data,wnd,ks)
    openQ=Array{Array{cable,1},1}()
    closedQ=Array{Array{cable,1},1}()
    cable_table=Array{cable,1}()

    #calculate cable at min of interval distance
    cbl=cable()
    km_min=cable_interval[1]
    cbl.mva=cable_interval[3]
    cbl.wnd=wnd
    cbl.length=km_min
    cbl_min=optimal_ac_cable(cbl,cbl_data,km_min,ks)

    #calculate cable at max of interval distance
    cbl=cable()
    km_max=cable_interval[2]
    cbl.mva=cable_interval[3]
    cbl.wnd=wnd
    cbl.length=km_max
    cbl_max=optimal_ac_cable(cbl,cbl_data,km_max,ks)

    if ((cbl_min.size==cbl_max.size) && (cbl_min.num==cbl_max.num))
        push!(cable_table,cbl_max)
    else
        push!(openQ,[cbl_min,cbl_max])
        while (length(openQ)>0)
            openQ=halfway_cable(openQ,cbl_data,ks)
            closedQ,openQ=update_Qs(closedQ,openQ)
        end
        cable_table_oversized=Array{cable,1}()
        for q in closedQ
            push!(cable_table_oversized,deepcopy(q[1]))
            push!(cable_table_oversized,deepcopy(q[2]))
        end
        push!(cable_table_oversized,cbl_max)
        for (index,cbl) in enumerate(cable_table_oversized[1:length(cable_table_oversized)-1])
            if ((cbl.size==cable_table_oversized[index+1].size) && (cbl.num==cable_table_oversized[index+1].num))
            else
                push!(cable_table,deepcopy(cbl))
            end
        end
        push!(cable_table,cbl_max)
    end
    return cable_table
end

#moves finished intervals form the openQ to the closedQ
#**
function update_Qs(closedQ,openQ)
    openQ_new=Array{Array{cable,1},1}()
    for q in openQ
        if ((q[1].size==q[2].size) && (q[1].num==q[2].num))
            push!(closedQ,deepcopy(q))
        elseif ((q[1].length+1>=q[2].length))
            push!(closedQ,deepcopy(q))
        else
            push!(openQ_new,deepcopy(q))
        end
    end
    return closedQ,openQ_new
end

#devides open q at mid point
#**
function halfway_cable(openQ,cbl_data,ks)
    openQ_new=Array{Array{cable,1},1}()
    for q in openQ
        cbl=cable()
        km_halfway=(q[1].length+q[2].length)/2
        cbl.mva=q[1].mva
        cbl.wnd=q[1].wnd
        cbl_halfway=optimal_ac_cable(cbl,cbl_data,km_halfway,ks)
        if ((cbl_halfway.size==q[2].size) && (cbl_halfway.num==q[2].num))
            push!(openQ_new,[q[1],cbl_halfway])
        elseif ((cbl_halfway.size==q[1].size) && (cbl_halfway.num==q[1].num))
            push!(openQ_new,[cbl_halfway,q[2]])
        else
            push!(openQ_new,[q[1],cbl_halfway])
            push!(openQ_new,[cbl_halfway,q[2]])
        end
    end
    return openQ_new
end

#finds the hvdc cable for the look up table
#**
function optimal_hvdc_cable(cbl0,cbl_data,km,ks)
    cbl=cable()
    cbl.costs.ttl=Inf
    cbl.elec.volt=cbl_data[1][1]
    cbl.length=km
    cbl.mva=cbl0.mva
    cbl0.costs.ttl=Inf

    for _cd in cbl_data
        num=2
        capacity=(_cd[1]*_cd[5])*10^-3
        while (cbl0.mva>num*capacity)
            num=num+2
        end
        fillOut_cable_struct_dc(cbl0,_cd,km,num)
        cbl0=cost_hvdc_cable(cbl0,ks)
        if (cbl0.costs.ttl<cbl.costs.ttl)
            cbl=deepcopy(cbl0)
        end
    end
    return cbl
end

############################## Change made here ########################### Temporary kana???
#returns the optimal cable for a hvac for the look up table
#**
function optimal_ac_cable(cbl0,cbl_data,km,ks)
    cbl=cable()
    cbl.costs.ttl=Inf
    cbl.elec.volt=cbl_data[1][1]
    cbl.length=km
    cbl.mva=cbl0.mva
    ##########################
    #=if (cbl0.elec.volt==66.0)
        cbl.mva=cbl0.mva*1.333
    else
        cbl.mva=cbl0.mva
    end=#
    ##########################
    cbl0.costs.ttl=Inf
    for cd in cbl_data
        #num=1
        ######################
        if (cbl0.elec.volt==66.0)
            num=6
        else
            num=1
        end
        ######################
        exception_66kv=false
        cap_at_km=get_newQ_Capacity(cbl0.elec.freq,km,cd[1],cd[4]*10^-9,cd[5])
        #max number of cables in parallel = 12
        max_in_parallel=12#if changing check economics main function mvac_cable(mva,km,wnd,cable_array,ks)
        GMAX=2000
        if (cap_at_km/cbl0.mva>1/max_in_parallel || (cbl.elec.volt==66.0 && cbl0.mva<=3000))
            while ((cbl0.mva>num*cap_at_km) && num<max_in_parallel)
                num=deepcopy(num+1)
            end
            if ((cbl.elec.volt==66.0 && cbl0.mva>num*cap_at_km))
                num=1
                exception_66kv=true
                cbl0.mva=cbl0.mva/2
                while ((cbl0.mva>num*cap_at_km) && num<max_in_parallel)
                    num=deepcopy(num+1)
                end
            end
            fillOut_cable_struct_ac(cbl0,cd,km,num)
            #cbl0=cost_hvac_cable_o2o(cbl0,ks)
            if (cbl0.elec.volt==66.0)
                cbl0=cost_mvac_cable(cbl0,ks)
                if (exception_66kv==true)
                    cbl0.num=cbl0.num*2
                    cbl0.mva=cbl0.mva*2
                    cbl0.costs.perkm_cpx=cbl0.costs.perkm_cpx*2
                    cbl0.costs.sg=cbl0.costs.sg*2
                    cbl0.costs.cpx_p=cbl0.costs.cpx_p*2
                    cbl0.costs.cpx_i=cbl0.costs.cpx_i*2
                    cbl0.costs.rlc=cbl0.costs.rlc*2
                    cbl0.costs.cm=cbl0.costs.cm*2
                    cbl0.costs.eens=cbl0.costs.eens*2
                    cbl0.costs.ttl=cbl0.costs.ttl*2
                    cbl0.costs.grand_ttl=cbl0.costs.ttl
                end
            elseif (cbl0.elec.volt==220.0 || cbl0.elec.volt==400.0)
                cbl0=cost_hvac_cable(cbl0,ks)
            else
            end
        end
        if (cbl0.costs.ttl<cbl.costs.ttl)
            cbl=deepcopy(cbl0)
        end
    end
    return cbl
end

#puts cables into dictionaries
#**
function make_dictionaries_cables(cable_table)
    all_cable_dict=Dict{String,Dict{String,Vector{cable}}}()
    for cable_sets in cable_table
        kv_cable_dict=Dict{String,Vector{cable}}()
        for cable_set in cable_sets
             push!(kv_cable_dict,(string(cable_set[1].mva)=>cable_set))
        end
        push!(all_cable_dict,(string(cable_sets[1][1].elec.volt)=>kv_cable_dict))
    end
    return all_cable_dict
end
################################################################################
######################## Transformers and converters ###########################
################################################################################
#**
function get_xfoConv_table(km_mva_set,wnd)
    ks=get_Cost_Data()
    #HVDC
    conv_table=Array{Array{converter,1},1}()
    for km_mva in km_mva_set
        conv_off=converter()
        conv_off.relia=get_offshore_conv_failure_data(conv_off)
        conv_off.wnd=wnd
        conv_off.mva=km_mva[2]
        conv_off.elec.mva=km_mva[2]
        conv_off.num=1
        conv_off=cost_hvdc_oss(conv_off,ks)
        conv_off=adjust_base_hvdc_offshore_converter(conv_off,ks)

        conv_on=converter()
        conv_on.relia=get_onshore_conv_failure_data(conv_on)
        conv_on.wnd=wnd
        conv_on.mva=km_mva[2]
        conv_on.elec.mva=km_mva[2]
        conv_on.num=1
        conv_on=cost_hvdc_pcc(conv_on,ks)
        conv_on=adjust_base_hvdc_onshore_converter(conv_on,ks)
        push!(conv_table,[conv_off,conv_on])
    end

    #HVAC
    xfo_table=Array{Array{transformer,1},1}()
    xfo_data=get_Xfo_Data()
    for km_mva in km_mva_set
        xfo_off=transformer()
        xfo_off.relia=get_offshore_xfo_failure_data(xfo_off)
        xfo_off.wnd=wnd
        xfo_off.mva=km_mva[2]
        xfo_off=xfo_oss(xfo_off,ks,xfo_data)

        xfo_on=transformer()
        xfo_on.relia=get_onshore_xfo_failure_data(xfo_on)
        xfo_on.wnd=wnd
        xfo_on.mva=km_mva[2]
        xfo_on=xfo_pcc(xfo_on,ks,xfo_data)
        push!(xfo_table,[xfo_off,xfo_on])
    end

    xfo_dictionary,conv_dictionary=make_dictionaries_xfo_conv(xfo_table,conv_table)
    return xfo_dictionary,conv_dictionary
end

#puts transformers and converters into dictionaries
#**
function make_dictionaries_xfo_conv(xfo_table,conv_table)
    all_xfo_dict=Dict{String,Dict{String,transformer}}()
    for xfo_sets in xfo_table
        xfo_dict=Dict{String,transformer}()
        push!(xfo_dict,("offshore"=>xfo_sets[1]))
        push!(xfo_dict,("onshore"=>xfo_sets[2]))
        push!(all_xfo_dict,(string(xfo_sets[1].mva)=>xfo_dict))
    end
    all_conv_dict=Dict{String,Dict{String,converter}}()
    for conv_sets in conv_table
        conv_dict=Dict{String,converter}()
        push!(conv_dict,("offshore"=>conv_sets[1]))
        push!(conv_dict,("onshore"=>conv_sets[2]))
        push!(all_conv_dict,(string(conv_sets[1].mva)=>conv_dict))
    end
    return all_xfo_dict,all_conv_dict
end
################################################################################
############################# Database Bits ####################################
################################################################################
#checks if HVDC is cheapest option for longest ranges at each power for a PCC connection
#**
function check4hvdc(equipment_database,km_mva_set,wnd,kv)
    hvdc_bit=false
    for km_mva in km_mva_set
        eqps=optimal_mog2pcc_wMOG(km_mva[1],km_mva[2],kv,wnd,equipment_database)
        cheapest=argmin([eqps["300.0"]["cost"],eqps["220.0"]["cost"],eqps["400.0"]["cost"]])
        cheapest_connection_2pcc=[eqps["300.0"],eqps["220.0"],eqps["400.0"]][cheapest]
        if (haskey(cheapest_connection_2pcc,"cable"))
            if (cheapest_connection_2pcc["cable"].elec.volt==300.0)
                hvdc_bit=true
            end
        end
    end
    return hvdc_bit
end

#checks if mid point compensation is cheapest option for longest ranges at each power for a PCC connection
#**
function check4mpc_ac(equipment_database,km_mva_set,wnd,kv)
    mpc_bit=false
    for km_mva in km_mva_set
        eqps=optimal_mog2pcc_wMOG(km_mva[1],km_mva[2],kv,wnd,equipment_database)
        cheapest=argmin([eqps["300.0"]["cost"],eqps["220.0"]["cost"],eqps["400.0"]["cost"]])
        cheapest_connection_2pcc=[eqps["300.0"],eqps["220.0"],eqps["400.0"]][cheapest]
        if (haskey(cheapest_connection_2pcc,"cable"))
            if (cheapest_connection_2pcc["cable"].mpc_ac==true)
                mpc_bit=true
            end
        end
    end
    return mpc_bit
end

#checks if hvac is needed if HVDC is possible
#**
function check4hvac(equipment_database,km_mva_set,wnd,kv)
    hvac_bit=false
    for km_mva in km_mva_set
        eqps=optimal_mog2pcc_wMOG(1,km_mva[2],kv,wnd,equipment_database)
        cheapest=argmin([eqps["300.0"]["cost"],eqps["220.0"]["cost"],eqps["400.0"]["cost"]])
        cheapest_connection_2pcc=[eqps["300.0"],eqps["220.0"],eqps["400.0"]][cheapest]
        if (haskey(cheapest_connection_2pcc,"plat_aft_ac"))
            if (cheapest_connection_2pcc["plat_aft_ac"].acdc=="ac")
                hvac_bit=true
            end
        end
    end
    return hvac_bit
end


text:
"A triangle defined by three points in a plane."
abstract type Triangle end
# This is just a conceptual type that acts as a data container. No functionality is implemented
# to manipulate triangles.

"The triangle below the diagonal of a square `[a,b] × [a,b]`."
struct LowerRightTriangle{T} <: Triangle
   a  :: T
   b  :: T
end

"The triangle above the diagonal of a square `[a,b] × [a,b]`."
struct UpperRightTriangle{T} <: Triangle
   a  :: T
   b  :: T
end


text:
r = Quaternion(rand(), ℝ³(rand(3)))
t = 100ℝ³(rand(3))
dq1 = Biquaternion(r, t)

rotation = getrotation(dq1)
translation = gettranslation(dq1)

# Unit condition
@test isapprox(norm(dq1), 1)
@test isapprox(conj(conj(dq1)), dq1)
@test isapprox(rotation, normalize(r))
@test isapprox(translation, t)

qr = Quaternion(rand(), ℝ³(rand(3)))
qd = Quaternion([0; vec(100ℝ³(rand(3)))])
dq2 = Biquaternion(qr, qd)
q̂ = normalize(dq2)
dq3 = Biquaternion(2dq2)

@test isapprox(norm(q̂), 1)
@test isapprox(dq2, dq3) != true

## constructors

rotation = Quaternion(rand(4))
translation = ℝ³(rand(3))
q1 = Biquaternion(rotation)
q2 = Biquaternion(translation)
q = Biquaternion(rotation, translation)

@test isapprox(q2 * q1, q)

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
module Wip
using Simulacrum
using HyperGraphs, Symbolics, ModelingToolkit
using LinearAlgebra: diagm
import SymbolicUtils: Term, Sym

export mass_action, power_law, combinatoric_power_law, rn_assumptions, make_equations, clone_map,
AbstractModel, AbstractBioModel, TelegraphModel, BirthDeathModel, hypergraph, make_Langevin_noise

# incrementally build common reaction network assumptions on chemical hyperedges
power_law(che::ChemicalHyperEdge{Num}) = src(che) .^ src_stoich(che)
parse_stoich(s) = s <= 20 ? s : big(s) # factorial for numbers > 20 needs big
combinatoric_power_law(che::ChemicalHyperEdge{Num}) = power_law(che) ./ factorial.(parse_stoich.(src_stoich(che)))
mass_action(che::ChemicalHyperEdge{Num}; f::Function = src) = rate(che) * prod(f(che))

# assemble reaction network assumptions on a chemical hyperedge
rn_assumptions(che::ChemicalHyperEdge{Num}) = mass_action(che, f = combinatoric_power_law)
rn_assumptions(chx::ChemicalHyperGraph{Num}) = rn_assumptions.(hyperedges(chx))

# helper function that returns equations describing how the system evolves over time for the given chemical hypergraph
function make_equations(chx::ChemicalHyperGraph{Num}, f::Function; iv::Symbol = :t)
    D    = Differential(iv)
    lhss = [D(x) for x in vertices(chx)]
    rhss = incidence_matrix(chx) * f(chx)
    eqs  = Equation.(lhss, rhss)
end

# overload the ModelingToolkit.ODESystem constructor to work on chemical hypergraphs
function ModelingToolkit.ODESystem(chx::ChemicalHyperGraph, f::Function; iv = Symbolics.variable(:t), states = vertices(chx), ps = [])
    ODESystem(make_equations(chx, f), iv, states, ps)
end
function ModelingToolkit.ODESystem(chx::ChemicalHyperGraph{Num}, f::Function; iv = Symbolics.variable(:t), states = vertices(chx), ps = [])
    ODESystem(Simulacrum.get_value(make_equations(chx, f)), iv, states, ps)
end
function ModelingToolkit.ODESystem(chx::ChemicalHyperGraph{T}, f::Function; iv = Symbolics.variable(:t), ps = []) where {T<:Term}
    ODESystem(convert(ChemicalHyperGraph{Num}, chx), f, iv = iv, ps = ps)
end

# automatically convert equations in matrix form to vector (convenience)
ModelingToolkit.ODESystem(eqs::Matrix{Equation}; iv = Symbolics.variable(:t), states = get_state_variables(vec(eqs)), ps = get_parameters(vec(eqs))) = ODESystem(vec(eqs), iv, states, ps)

# helper function that takes a hypergraph that has been made by cloning a base model n times and returns a map to be used with ModelingToolkit
function clone_map(in_map::Vector{Pair{T, U}}, clones::Vector) where {T, U}
    out_map = Vector{eltype(in_map)}()
    for i in 1:length(in_map):length(clones)
        for j in 1:length(in_map)
            push!(out_map, (clones[i + j - 1] => in_map[j].second))
        end
    end
    out_map
end
function clone_map(in_map::Vector{Pair{T, U}}, x::V) where {T, U, V<:AbstractHyperGraph}
    clone_map(in_map, vertices(x))
end

# helper function that makes maps to be used with e.g. ModelingToolkit
function make_map(firsts::Vector, seconds::Vector)
    @assert length(firsts) == length(seconds)
    [(firsts[i] => seconds[i]) for i in 1:length(firsts)]
end

# dispatch susbstitute on vectors of equations (convenience)
Symbolics.substitute(eqs::Vector{Equation}, d::Dict) = [substitute(eq, d) for eq in eqs]
Symbolics.substitute(nums::Vector{Num}, d::Dict) = [substitute(num, d) for num in nums]

# helper that returns the left-hand sides of equations wrapped in Num type
get_lhs_var(eq::Equation) = Num(Symbolics.get_variables(eq.lhs)[1])
get_lhss_vars(eqs::Vector{Equation}) = get_lhs_var.(eqs)

# helper that returns equations state variables
function get_state_variables(eq::Equation)
    all_vars = Num.(union(Symbolics.get_variables(eq.rhs)))
    filter(x -> Simulacrum.get_value(x) isa Term, all_vars)
end
get_state_variables(eqs::Vector{Equation}) = union(vcat(get_state_variables.(eqs)...))

# helper that returns equations parameters
function get_parameters(eq::Equation)
    all_vars = Num.(union(Symbolics.get_variables(eq.rhs)))
    filter(x -> Simulacrum.get_value(x) isa Sym, all_vars)
end
get_parameters(eqs::Vector{Equation}) = union(vcat(get_parameters.(eqs)...))

# prototype of models.jl
const ChE = ChemicalHyperEdge
const ChX = ChemicalHyperGraph
abstract type AbstractModel end
abstract type AbstractBioModel <: AbstractModel end
abstract type AbstractPhysModel <: AbstractModel end
hypergraph(m::T) where {T<:AbstractModel} = m.X
make_equations(m::T, f::Function) where {T<:AbstractModel} = make_equations(hypergraph(m), f)
"""
    TelegraphModel

Telegraph model of mRNA production.

active and inactive variables represent time spent in either state; this means that the transition representing transcription is abstract in nature: when the promoter is active, mRNA is transcribed (but not from the active state variable, and so it must appear on either side of the transition).
"""
struct TelegraphModel <: AbstractBioModel
    name::Symbol
    X::AbstractHyperGraph
    TelegraphModel(; name = :Telegraph) = new(name, make_telegraph_x())
end
function make_telegraph_x()
    @variables t inactive(t) active(t) mRNA(t) λ μ k δ
    es = [  ChE([inactive], [active], λ),     # promoter: inactive to active
            ChE([active], [inactive], μ),     # promoter: active to inactive
            ChE([active], [mRNA, active], k), # transcription
            ChE([mRNA], Num[], δ)]            # degradation: mRNA to nothing
    ChX(es)
end
"""
    BirthDeathModel

"""
struct BirthDeathModel <: AbstractBioModel
    name::Symbol
    X::AbstractHyperGraph
    function BirthDeathModel(; name = :BirthDeathModel, var::Symbol = :X)
        @variables t λ μ
        x = Num(Symbolics.variable(var, T = Symbolics.FnType)(t))
        es = [  ChE(Num[], [x], λ),
                ChE([x], Num[], μ)]
        new(name, ChX(es))
    end
end
"""
    LorenzSystemModel

Model of the Lorenz system, known for its chaotic solutions.
"""
struct LorenzSystemModel <: AbstractPhysModel
    name::Symbol
    X::AbstractHyperGraph
    function LorenzSystemModel(; name = :LorenzSystemModel)
        @variables t x(t) y(t) z(t) β ρ σ
        es = [  # dissipation
                ChE([x], Num[], σ),
                ChE([y], Num[], 1),
                ChE([z], Num[], β),

                # one-way interactions
                ChE([x], [x, y], ρ),
                ChE([y], [y, x], σ),

                # two-way interaction
                ChE([x, y], [x, y, z], 1),

                # three-way interaction
                ChE([x, y, z], [x, z], 1/y)]
        new(name, ChX(es))
    end
end

# helper that returns an Expr that can be evaluated to expose variables into the global scope, specifically written for Num types (for now)
expose_sym(sym::Num) = :( $(sym.val.name) = $sym )
expose_term(term::Num) = :( $(term.val.f.name) = $term )
expose(num::Num) = (num.val isa Sym ? expose_sym : expose_term)(num)
expose(nums::Vector{Num}) = expose.(nums)
macro expose(arg) # only works for one symbolic variable at a time...
    expose(eval(arg)) |> esc
end

# returns a matrix of Langevin (intrinsic) noise terms
function make_Langevin_noise(x::T, propensities_f::Function; complex = false) where {T<:AbstractHyperGraph}
    !(eltype(x) <: SymTypes) && error("expected a symbolic hypergraph, got $eltype(x)")
    incidence_matrix(x) * diagm(0 => sqrt.((complex ? identity : abs).(propensities_f(x))))
end
function make_Langevin_noise(m::T, propensities_f::Function; complex = false) where {T<:AbstractModel}
    make_Langevin_noise(hypergraph(m), propensities_f, complex = complex)
end

# more methods for get_value
Simulacrum.get_value(ns::Vector{Num}) = Simulacrum.get_value.(ns)
Simulacrum.get_value(ns::Matrix{Num}) = reshape(Simulacrum.get_value(vec(ns)), size(ns))

end # module


text:
module problem127

description = """
abc-hits
The radical of n, rad(n), is the product of distinct prime factors of n. For example, 504 = 23 × 32 × 7, so rad(504) = 2 × 3 × 7 = 42.

We shall define the triplet of positive integers (a, b, c) to be an abc-hit if:
GCD(a, b) = GCD(a, c) = GCD(b, c) = 1
a < b
a + b = c
rad(abc) < c

For example, (5, 27, 32) is an abc-hit, because:
GCD(5, 27) = GCD(5, 32) = GCD(27, 32) = 1
5 < 27
5 + 27 = 32
rad(4320) = 30 < 32

It turns out that abc-hits are quite rare and there are only thirty-one abc-hits for c < 1000, with ∑c = 12523.
Find ∑c for c < 120000.
"""

# calculate rad(x) for x=1:n using prime sieve
function radsieve(n :: Int)
  r = ones(Int, n)
  for i = 2:n
    if r[i] == 1
      for j = i:i:n
        r[j] *= i
      end
    end
  end
  return r
end

# Calculate sum(c) for c < cmax and a,b,c is an ABC-Hit
function abchits(cmax :: Int = 120000)
  # Pre-compute rad, and sort 1:cmax by rad
  r = radsieve(cmax)
  p = sortperm(r)

  csum = 0
  for c = 3:cmax
    radc = r[c]

    # rad(a*b*c) == rad(a) * rad(b) * rad(c) when a,b,c are coprime
    # rad(a)*rad(b)*rad(c) < c implies rad(a) < c / (rad(c) * rad(b))
    # using 2 as a lower bound for rad(b), max value for rad(a) is c / (2 * rad(c))
    radamax = div(c, radc * 2)
    amax = div(c, 2)

    # search the sorted list for suitable rad(a), a
    for i = 1:cmax
      a = p[i]
      rada = r[a]
      if rada > radamax break end # stop searching once rad(a) is too large
      if a  > amax continue end   # skip over a that are too large
      b = c - a
      radb = r[b]

      # gcd(a,c) == 1 implies gcd(a,b) == 1 and gcd(b,c) == 1
      # if there were some d that divides a and b, then d would also divide (a+b)
      if rada * radb * radc < c && gcd(a, c) == 1
        csum += c
      end
    end
  end

  return csum
end

end

max_stars_repo_name: mbuhot/mbuhot-euler-solutions

max_issues_repo_name: mbuhot/mbuhot-euler-solutions

max_forks_repo_name: mbuhot/mbuhot-euler-solutions


text:
#################### Posterior Statistics ####################

"""
    autocor(chains[; lags = [1, 5, 10, 50], demean = true, append_chains = true, kwargs...])

Compute the autocorrelation of each parameter for the chain.

Setting `append_chains=false` will return a vector of dataframes containing the
autocorrelations for each chain.
"""
function autocor(
    chains::Chains;
    lags::AbstractVector{<:Integer} = [1, 5, 10, 50],
    demean::Bool = true,
    append_chains = false,
    kwargs...
)
    funs = Function[]
    func_names = @. Symbol("lag ", lags)
    for i in lags
        push!(funs, x -> autocor(x, [i], demean=demean)[1])
    end

    return summarize(
        chains, funs...;
        func_names = func_names,
        append_chains = append_chains,
        name = "Autocorrelation",
        kwargs...
    )
end

"""
    cor(chains[; sections, append_chains = true, kwargs...])

Compute the Pearson correlation matrix for the chain.

Setting `append_chains=false` will return a vector of dataframes containing a correlation
matrix for each chain.
"""
function cor(
    chains::Chains;
    sections = _default_sections(chains),
    append_chains = true,
    kwargs...
)
    # Subset the chain.
    _chains = Chains(chains, _clean_sections(chains, sections))

    # Obstain names of parameters.
    names_of_params = names(_chains)

    if append_chains
        df = chaindataframe_cor("Correlation", names_of_params, to_matrix(_chains))
        return df
    else
        vector_of_df = [
            chaindataframe_cor(
                "Correlation - Chain $i", names_of_params, data
            )
            for (i, data) in enumerate(to_vector_of_matrices(_chains))
        ]
        return vector_of_df
    end
end

function chaindataframe_cor(name, names_of_params, chains::AbstractMatrix; kwargs...)
    # Compute the correlation matrix.
    cormat = cor(chains)

    # Summarize the results in a named tuple.
    nt = (; parameters = names_of_params,
          zip(names_of_params, (cormat[:, i] for i in axes(cormat, 2)))...)

    # Create a ChainDataFrame.
    return ChainDataFrame(name, nt; kwargs...)
end

"""
    changerate(chains[; sections, append_chains = true, kwargs...])

Compute the change rate for the chain.

Setting `append_chains=false` will return a vector of dataframes containing the change
rates for each chain.
"""
function changerate(
    chains::Chains{<:Real};
    sections = _default_sections(chains),
    append_chains = true,
    kwargs...
)
    # Subset the chain.
    _chains = Chains(chains, _clean_sections(chains, sections))

    # Obstain names of parameters.
    names_of_params = names(_chains)

    if append_chains
        df = chaindataframe_changerate("Change Rate", names_of_params, _chains.value.data)
        return df
    else
        vector_of_df = [
            chaindataframe_changerate(
                "Change Rate - Chain $i", names_of_params, data
            )
            for (i, data) in enumerate(to_vector_of_matrices(_chains))
        ]
        return vector_of_df
    end
end

function chaindataframe_changerate(name, names_of_params, chains; kwargs...)
    # Compute the change rates.
    changerates, mvchangerate = changerate(chains)

    # Summarize the results in a named tuple.
    nt = (; zip(names_of_params, changerates)..., multivariate = mvchangerate)

    # Create a ChainDataFrame.
    return ChainDataFrame(name, nt; kwargs...)
end

changerate(chains::AbstractMatrix{<:Real}) = changerate(reshape(chains, Val(3)))
function changerate(chains::AbstractArray{<:Real,3})
    niters, nparams, nchains = size(chains)

    changerates = zeros(nparams)
    mvchangerate = 0.0

    for chain in 1:nchains, iter in 2:niters
        isanychanged = false

        for param in 1:nparams
            # update if the sample is different from the one in the previous iteration
            if chains[iter-1, param, chain] != chains[iter, param, chain]
                changerates[param] += 1
                isanychanged = true
            end
        end

        mvchangerate += isanychanged
    end

    factor = nchains * (niters - 1)
    changerates ./= factor
    mvchangerate /= factor

    changerates, mvchangerate
end

describe(c::Chains; args...) = describe(stdout, c; args...)

"""
    describe(io, chains[;
             q = [0.025, 0.25, 0.5, 0.75, 0.975],
             etype = :bm,
             kwargs...])

Print the summary statistics and quantiles for the chain.
"""
function describe(
    io::IO,
    chains::Chains;
    q = [0.025, 0.25, 0.5, 0.75, 0.975],
    etype = :bm,
    kwargs...
)
    dfs = vcat(summarystats(chains; etype = etype, kwargs...),
               quantile(chains; q = q, kwargs...))
    return dfs
end

function _hpd(x::AbstractVector{<:Real}; alpha::Real=0.05)
    n = length(x)
    m = max(1, ceil(Int, alpha * n))

    y = sort(x)
    a = y[1:m]
    b = y[(n - m + 1):n]
    _, i = findmin(b - a)

    return [a[i], b[i]]
end

function hpd(chn::Chains; alpha::Real=0.05, kwargs...)
    labels = [:lower, :upper]
    l(x) = _hpd(x, alpha=alpha)[1]
    u(x) = _hpd(x, alpha=alpha)[2]
    return summarize(chn, l, u; name = "HPD", func_names = labels, kwargs...)
end

"""
    quantile(chains[; q = [0.025, 0.25, 0.5, 0.75, 0.975], append_chains = true, kwargs...])

Compute the quantiles for each parameter in the chain.

Setting `append_chains=false` will return a vector of dataframes containing the quantiles
for each chain.
"""
function quantile(
    chains::Chains;
    q::AbstractVector = [0.025, 0.25, 0.5, 0.75, 0.975],
    append_chains = true,
    kwargs...
)
    # compute quantiles
    funs = Function[]
    func_names = @. Symbol(100 * q, :%)
    for i in q
        push!(funs, x -> quantile(cskip(x), i))
    end

    return summarize(
        chains, funs...;
        func_names = func_names,
        append_chains = append_chains,
        name = "Quantiles",
        kwargs...
    )
end

"""
    function summarystats(
        chains;
        sections = _default_sections(chains),
        append_chains= true,
        method::AbstractESSMethod = ESSMethod(),
        maxlag = 250,
        etype = :bm,
        kwargs...
    )

Compute the mean, standard deviation, naive standard error, Monte Carlo standard error,
and effective sample size for each parameter in the chain.

Setting `append_chains=false` will return a vector of dataframes containing the summary
statistics for each chain.

When estimating the effective sample size, autocorrelations are computed for at most `maxlag` lags.
"""
function summarystats(
    chains::Chains;
    sections = _default_sections(chains),
    append_chains::Bool = true,
    method::AbstractESSMethod = ESSMethod(),
    maxlag = 250,
    etype = :bm,
    kwargs...
)
    # Store everything.
    funs = [mean∘cskip, std∘cskip, sem∘cskip, x -> mcse(cskip(x), etype; kwargs...)]
    func_names = [:mean, :std, :naive_se, :mcse]

    # Subset the chain.
    _chains = Chains(chains, _clean_sections(chains, sections))

    # Calculate ESS separately.
    ess_df = ess(_chains; sections = nothing, method = method, maxlag = maxlag)

    # Summarize.
    summary_df = summarize(
        _chains, funs...;
        func_names = func_names,
        append_chains = append_chains,
        additional_df = ess_df,
        name = "Summary Statistics",
        sections = nothing
    )

    return summary_df
end

"""
    mean(chains[, params; kwargs...])

Calculate the mean of a chain.
"""
function mean(chains::Chains; kwargs...)
    # Store everything.
    funs = [mean∘cskip]
    func_names = [:mean]

    # Summarize.
    summary_df = summarize(
        chains, funs...;
        func_names = func_names,
        name = "Mean",
        kwargs...
    )

    return summary_df
end

mean(chn::Chains, syms) = mean(chn[:, syms, :])
# resolve method ambiguity with `mean(f, ::AbstractArray)`
mean(chn::Chains, syms::AbstractVector) = mean(chn[:, syms, :])


text:
# https://discourse.julialang.org/t/write-a-rest-interface-like-flask/18538/16

# Bukdu v0.4.4
using Bukdu

struct RESTController <: ApplicationController
    conn::Conn
end

function create(c::RESTController)
    @info :payload (c.params.message, c.params.x, c.params.y)
    render(JSON, "OK")
end

routes() do
    post("/messages", RESTController, create)
    plug(Plug.Parsers, [:json])
end

Bukdu.start(8080)

#=
curl -H "Content-Type: application/json" http://127.0.0.1:8080/messages -d '{"message": "Hello Data"}'
curl -H "Content-Type: application/json" http://127.0.0.1:8080/messages -d '{"x": 1.2, "y": 2.3}'
=#


text:
# main term structures

type NullYieldTermStructure <: YieldTermStructure end

# function calculate!{T <: TermStructure}(ts::T, recalculate::Bool=false)
#   if !ts.calculated || recalculate
#     _calculate!(ts)
#   end
#
#   return ts
# end

type JumpDate
  ts_quote::Quote
  ts_date::Date
end

type JumpTime
  ts_quote::Quote
  ts_time::Float64
end

discount{T <: YieldTermStructure}(yts::T, date::Date) = discount(yts, time_from_reference(yts, date))

function discount{T <: YieldTermStructure}(yts::T, time_frac::Float64)
  disc = discount_impl(yts, time_frac)
  if isdefined(yts, :jumpTimes)
    if length(yts.jumpTimes) == 0
      return disc
    end

    jump_effect = 1.0
    for jump in yts.jumpTimes
      if jump.ts_time > 0.0 && jump.ts_time < time
        if jump.ts_quote.value > 0.0 && jump.ts_quote.value <= 1.0
          jump_effect *= jump.ts_quote.value
        end
      end
    end

    return jump_effect * disc
  else
    return disc
  end
end

function zero_rate{T <: YieldTermStructure, C <: CompoundingType, F <: Frequency}(yts::T, date::Date, dc::DayCount, comp::C, freq::F = Annual())
  if date == yts.referenceDate
    return implied_rate(1.0 / discount(yts, 0.0001), dc, comp, 0.0001, freq)
  else
    return implied_rate(1.0 / discount(yts, date), dc, comp, reference_date(yts), date, freq)
  end
end

function zero_rate{T <: YieldTermStructure, C <: CompoundingType, F <: Frequency}(yts::T, time_frac::Float64, comp::C, freq::F = Annual())
  t = time_frac == 0.0 ? 0.0001 : time_frac
  return implied_rate(1.0 / discount(yts, t), yts.dc, comp, t, freq)
end

function forward_rate{T <: YieldTermStructure, DC <: DayCount, C <: CompoundingType, F <: Frequency}(yts::T, date1::Date, date2::Date, dc::DC, comp::C, freq::F)
  if date1 == date2
    t1 = max(time_from_reference(yts, date1) - 0.0001 / 2.0, 0.0)
    t2 = t1 + 0.0001
    return implied_rate(discount(yts, t1) / discount(yts, d2), dc, comp, 0.0001, freq)
  elseif date1 < date2
    return implied_rate(discount(yts, date1) / discount(yts, date2), dc, comp, date1, date2, freq)
  else
    error("Forward start date must be before forward end date")
  end
end

forward_rate{T <: YieldTermStructure, DC <: DayCount, C <: CompoundingType, F <: Frequency}(yts::T, date::Date, period::Integer, dc::DC, comp::C, freq::F) = forward_rate(yts, date, date + Dates.Day(period), dc, comp, freq)

function forward_rate{T <: YieldTermStructure, C <: CompoundingType, F <: Frequency}(yts::T, time1::Float64, time2::Float64, comp::C, freq::F)
  if time1 == time2
    t1 = max(time1 - 0.0001 / 2.0, 0.0)
    t2 = t1 + 0.0001
    interval, compound = (t2 - t1, discount(yts, t1) / discount(yts, t2))
  else
    interval, compound = (time2 - time1, discount(yts, time1) / discount(yts, time2))
  end

  return implied_rate(compound, yts.dc, comp, interval, freq)
end

## FlatForwardTermStructure
type FlatForwardTermStructure{B <: BusinessCalendar, DC <: DayCount, C <: CompoundingType, F <: Frequency} <: YieldTermStructure
  settlementDays::Int
  referenceDate::Date
  calendar::B
  forward::Quote
  dc::DC
  comp::C
  freq::F
  rate::InterestRate
  jumpTimes::Vector{JumpTime}
  jumpDates::Vector{JumpDate}
end

function FlatForwardTermStructure{B <: BusinessCalendar, DC <: DayCount, C <: CompoundingType, F <: Frequency}(settlement_days::Int, referenceDate::Date, calendar::B, forward::Quote, dc::DC,
                                  comp::C = ContinuousCompounding(), freq::F = QuantLib.Time.Annual())
  rate = InterestRate(forward.value, dc, comp, freq)
  FlatForwardTermStructure{B, DC, C, F}(settlement_days, referenceDate, calendar, forward, dc, comp, freq, rate, Vector{JumpTime}(0), Vector{JumpDate}(0))
end

FlatForwardTermStructure{B <: BusinessCalendar, DC <: DayCount, C <: CompoundingType, F <: Frequency}(referenceDate::Date, calendar::B, forward::Quote, dc::DC, comp::C = ContinuousCompounding(), freq::F = QuantLib.Time.Annual()) =
                        FlatForwardTermStructure(0, referenceDate, calendar, forward, dc, comp, freq)

FlatForwardTermStructure{B <: BusinessCalendar, DC <: DayCount, C <: CompoundingType, F <: Frequency}(settlementDays::Int, calendar::B, forward::Quote, dc::DC, comp::C = ContinuousCompounding(), freq::F = QuantLib.Time.Annual()) =
                        FlatForwardTermStructure(settlementDays, Date(), calendar, forward, dc, comp, freq)

FlatForwardTermStructure{DC <: DayCount}(referenceDate::Date, forward::Float64, dc::DC) =
                        FlatForwardTermStructure(0, referenceDate, TargetCalendar(), Quote(forward), dc, ContinuousCompounding(), Annual())

FlatForwardTermStructure{DC <: DayCount, C <: CompoundingType, F <: Frequency}(referenceDate::Date, forward::Float64, dc::DC, compounding::C, freq::F) =
                        FlatForwardTermStructure(0, referenceDate, TargetCalendar(), Quote(forward), dc, compounding, freq)

discount_impl(ffts::FlatForwardTermStructure, time_frac::Float64) = discount_factor(ffts.rate, time_frac)

# discount_impl(ffts::FlatForwardTermStructure, time_frac::Float64) = discount_factor(ffts.rate, time_frac)


text:
function pr_comment_is_blocking(c::GitHub.Comment)
    c_body = body(c)
    if occursin("[noblock]", c_body)
        return false
    else
        return true
    end
end

function pr_has_no_blocking_comments(registry::GitHub.Repo,
                                     pr::GitHub.PullRequest;
                                     auth::GitHub.Authorization)
    all_pr_comments = get_all_pull_request_comments(registry, pr; auth=auth)
    if isempty(all_pr_comments)
        return true
    else
        num_comments = length(all_pr_comments)
        comment_is_blocking = BitVector(undef, num_comments)
        for i = 1:num_comments
            comment_is_blocking[i] = pr_comment_is_blocking(all_pr_comments[i])
        end
        if any(comment_is_blocking)
            return false
        else
            return true
        end
    end
end

function pr_is_old_enough(pr_type::Symbol,
                          pr_age::Dates.Period;
                          new_package_waiting_period::Dates.Period,
                          new_version_waiting_period::Dates.Period)
    if pr_type == :NewPackage
        return pr_age > new_package_waiting_period
    elseif pr_type == :NewVersion
        return pr_age > new_version_waiting_period
    else
        throw(ArgumentError("pr_type must be either :NewPackage or :NewVersion"))
    end
end

@inline my_isless(r1::GitHub.Review, r2::GitHub.Review) = isless(r1.id, r2.id)

function pr_was_approved_by_me(repo::GitHub.Repo,
                               pr::GitHub.PullRequest;
                               auth::GitHub.Authorization,
                               whoami::String)
    all_pr_reviews = get_all_pull_request_reviews(repo, pr; auth=auth)
    my_pr_reviews = Vector{GitHub.Review}(undef, 0)
    i_left_at_least_one_approving_review = false
    i_did_not_leave_any_nonapproving_reviews = true
    for rev in all_pr_reviews
        if reviewer_login(rev) == whoami
            if review_state(rev) == "APPROVED"
                i_left_at_least_one_approving_review = true
            else
                i_did_not_leave_any_nonapproving_reviews = false
            end
            push!(my_pr_reviews, rev)
        end
    end
    i_approved_this_pr = i_left_at_least_one_approving_review && i_did_not_leave_any_nonapproving_reviews
    if i_approved_this_pr
        if length(my_pr_reviews) == 1
            my_review = my_pr_reviews[1]
        else
            @warn("I did not leave exactly one review. Using the most recent review")
            sort!(my_pr_reviews; lt = my_isless)
            my_review = my_pr_reviews[end]
        end
        my_review_body = my_review.body
        approving_review_regex = r"approved_pr_head_commit_sha=\"([A-Za-z0-9]*?)\""
        if occursin(approving_review_regex, my_review_body)
            m = match(approving_review_regex, my_review_body)
            approved_pr_head_sha = m[1]
            return true, approved_pr_head_sha
        else
            @error("The review I left does not match the regex", my_review_body)
            return false, ""
        end
    else
        return false, ""
    end
end

function cron_or_api_build(registry::GitHub.Repo;
                           auth::GitHub.Authorization,
                           authorized_authors::Vector{String},
                           merge_new_packages::Bool,
                           merge_new_versions::Bool,
                           new_package_waiting_period,
                           new_version_waiting_period,
                           whoami::String)
    # first, get a list of ALL open pull requests on this repository
    # then, loop through each of them.
    all_currently_open_pull_requests = my_retry(() -> get_all_pull_requests(registry, "open"; auth=auth))
    reverse!(all_currently_open_pull_requests)
    at_least_one_exception_was_thrown = false
    num_retries = 0
    if isempty(all_currently_open_pull_requests)
        @info("There are no open pull requests.")
    else
        for pr in all_currently_open_pull_requests
            try
                my_retry(() -> cron_or_api_build(pr,
                                                 registry::GitHub.Repo;
                                                 auth = auth,
                                                 authorized_authors = authorized_authors,
                                                 merge_new_packages = merge_new_packages,
                                                 merge_new_versions = merge_new_versions,
                                                 new_package_waiting_period = new_package_waiting_period,
                                                 new_version_waiting_period = new_version_waiting_period,
                                                 whoami = whoami),
                         num_retries)
            catch ex
                at_least_one_exception_was_thrown = true
                showerror(stderr, ex)
                Base.show_backtrace(stderr, catch_backtrace())
                println(stderr)
            end
        end
        if at_least_one_exception_was_thrown
            error("At least one exception was thrown. Check the logs for details.")
        end
    end
    return nothing
end

function cron_or_api_build(pr::GitHub.PullRequest,
                           registry::GitHub.Repo;
                           auth::GitHub.Authorization,
                           authorized_authors::Vector{String},
                           merge_new_packages::Bool,
                           merge_new_versions::Bool,
                           new_package_waiting_period,
                           new_version_waiting_period,
                           whoami::String)
    #       first, see if the author is an approved author. if not, then skip.
    #       next, see if the title matches either the "New Version" regex or
    #               the "New Package regex". if it is not either a new
    #               package or a new version, skip.
    #       next, see if it is old enough. if it is not old enough, then skip.
    #       then, get all of the reviews. make sure that (1) I left at least one
    #               review, and (2) all of my reviews are approving. if this criterion
    #               is not met, skip
    #       then, get all of the pull request comments. if there is any comment that is
    #               (1) not by me, and (2) does not contain the text [noblock], then skip
    #       if all of the above criteria were met, then merge the pull request
    pr_number = number(pr)
    @info("Now examining pull request $(pr_number)")
    pr_author = author_login(pr)
    if pr_author in authorized_authors
        if is_new_package(pr) || is_new_version(pr)
            if is_new_package(pr) # it is a new package
                pr_type = :NewPackage
                pkg, version = parse_pull_request_title(NewPackage(), pr)
            else # it is a new version
                pr_type = :NewVersion
                pkg, version = parse_pull_request_title(NewVersion(), pr)
            end
            pr_age = time_since_pr_creation(pr)
            _pr_is_old_enough = pr_is_old_enough(pr_type,
                                                 pr_age;
                                                 new_package_waiting_period = new_package_waiting_period,
                                                 new_version_waiting_period = new_version_waiting_period)
            if _pr_is_old_enough
                i_approved_this_pr, approved_pr_head_sha = pr_was_approved_by_me(registry,
                                                                                pr;
                                                                                auth = auth,
                                                                                whoami = whoami)
                if i_approved_this_pr
                    if pr_has_no_blocking_comments(registry, pr; auth = auth)
                        "Pull request: $(pr_number). "
                        "Type: $(pr_type). "
                        "Decision: merge. "
                        if pr_type == :NewPackage # it is a new package
                            if merge_new_packages
                                my_comment = comment_text_merge_now()
                                @info(string("Pull request: $(pr_number). ",
                                             "Type: $(pr_type). ",
                                             "Decision: merge now."))
                                # my_retry(() -> post_comment!(registry, pr, my_comment; auth = auth))
                                my_retry(() -> merge!(registry, pr, approved_pr_head_sha; auth = auth))
                            else
                                @info(string("Pull request: $(pr_number). ",
                                             "Type: $(pr_type). ",
                                             "Decision: do not merge. ",
                                             "Reason: ",
                                             "This is a new package pull request. ",
                                             "All of the criteria for automerging ",
                                             "were met. ",
                                             "However, merge_new_packages is false, ",
                                             "so I will not merge. ",
                                             "If merge_new_packages had been set to ",
                                             "true, I would have merged this ",
                                             "pull request right now."))
                            end
                        else # it is a new version
                            if merge_new_versions
                                my_comment = comment_text_merge_now()
                                @info(string("Pull request: $(pr_number). ",
                                             "Type: $(pr_type). ",
                                             "Decision: merge now."))
                                # my_retry(() -> post_comment!(registry, pr, my_comment; auth = auth))
                                my_retry(() -> merge!(registry, pr, approved_pr_head_sha; auth = auth))
                            else
                                @info(string("Pull request: $(pr_number). ",
                                             "Type: $(pr_type). ",
                                             "Decision: do not merge. ",
                                             "Reason: merge_new_versions is false",
                                             "This is a new version pull request. ",
                                             "All of the criteria for automerging ",
                                             "were met. ",
                                             "However, merge_new_versions is false, ",
                                             "so I will not merge. ",
                                             "If merge_new_versions had been set to ",
                                             "true, I would have merged this ",
                                             "pull request right now."))
                            end
                        end
                    else
                        @info(string("Pull request: $(pr_number). ",
                                     "Type: $(pr_type). ",
                                     "Decision: do not merge. ",
                                     "Reason: pull request has one or more blocking comments."))
                    end
                else
                    @info(string("Pull request: $(pr_number). ",
                                 "Type: $(pr_type). ",
                                 "Decision: do not merge. ",
                                 "Reason: it is not the case that ",
                                 "both of the following conditions ",
                                 "are true: ",
                                 "(1) I left at least one approving ",
                                 "review. ",
                                 "(2) I did not leave any ",
                                 "non-approving reviews."),
                          whoami)
                end
            else
                @info(string("Pull request: $(pr_number). ",
                             "Type: $(pr_type). ",
                             "Decision: do not merge. ",
                             "Reason: mandatory waiting period has not elapsed."),
                      pr_type,
                      pr_age,
                      new_package_waiting_period,
                      new_version_waiting_period)
            end
        else
            @info(string("Pull request: $(pr_number). ",
                         "Decision: do not merge. ",
                         "Reason: pull request is neither a new package nor a new version."),
                  title(pr))
        end
    else
        @info(string("Pull request: $(pr_number). ",
                     "Decision: do not merge. ",
                     "Reason: pull request author is not authorized to automerge."),
              pr_author,
              authorized_authors)
    end
    return nothing
end


text:
using Test

include("interp.jl")

@test interp(numC(5), top_level_env) == numV(5)
@test interp(numC(0), top_level_env) == numV(0)

@test interp(ifC(idC("true"), numC(-1), numC(1)), top_level_env) == numV(-1)
@test interp(ifC(idC("false"), numC(-1), numC(1)), top_level_env) == numV(1)
@test_throws ErrorException("Invalid TULI5 if statement") interp(ifC(numC(324), numC(-1), numC(1)), top_level_env)

@test interp(appC(idC("<="), [numC(3), numC(2)]), top_level_env) == boolV(false)
@test interp(appC(idC("<="), [numC(3), numC(3)]), top_level_env) == boolV(true)

@test interp(appC(idC("+"), [numC(3), numC(5)]), top_level_env) == numV(8)
@test interp(appC(idC("-"), [numC(15), numC(4.5)]), top_level_env) == numV(10.5)
@test interp(appC(idC("*"), [numC(3), numC(31)]), top_level_env) == numV(93)
@test interp(appC(idC("/"), [numC(10), numC(5)]), top_level_env) == numV(2.0)
@test interp(appC(idC("equal?"), [numC(10), numC(5)]), top_level_env) == boolV(false)
@test interp(appC(idC("equal?"), [numC(10), numC(10)]), top_level_env) == boolV(true)

@test_throws ErrorException("TULI5 invalid + call") interp(appC(idC("+"), [boolV(true), numC(5)]), top_level_env)
@test_throws ErrorException("TULI5 invalid - call") interp(appC(idC("-"), [boolV(true), numC(5)]), top_level_env)
@test_throws ErrorException("TULI5 invalid * call") interp(appC(idC("*"), [boolV(true), numC(5)]), top_level_env)
@test_throws ErrorException("TULI5 invalid <= call") interp(appC(idC("<="), [boolV(true), numC(5)]), top_level_env)
@test_throws ErrorException("TULI5 invalid / call") interp(appC(idC("/"), [boolV(true), numC(5)]), top_level_env)
@test_throws ErrorException("TULI5 division by 0") interp(appC(idC("/"), [numC(6), numC(0)]), top_level_env)

@test interp(appC(lamC(["a", "b"], appC(idC("+"), [idC("a"), idC("b")])), [numC(3), numC(7)]), top_level_env) == numV(10)
@test interp(appC(lamC(["a", "b"], appC(idC("a"), [idC("b")])), [lamC(["x"], appC(idC("*"), [idC("x"), idC("x")])), numC(10)]), top_level_env) == numV(100)

max_stars_repo_name: nhausmna/julia-hello-world

max_issues_repo_name: nhausmna/julia-hello-world

max_forks_repo_name: nhausmna/julia-hello-world


text:
# using CuArrays
# using CUDAnative
import CUDA.CUDNN

# CuArrays doesn't provide a dedicated version of conv() for CuArrays, and so call to
# NNlib.conv(CuArray(), CuArray()) leads to CPU-based algorithm on GPU arrays, which is terribly slow.
# Here we implement our wrapper conv2d() for CuArrays via in-place version CuArrays.conv!().
# The same applies to maxpool() which we implement using maxpool!().
# Note that later we will have in-place rules and rewrite them on the tape, but we still
# need to make all operations works without any reference to gradients.

## conv2d

# function conv2d(x::CuArray{T,4}, w::CuArray{T,4}; stride=1, padding=0, dilation=1) where T
#     cdims = DenseConvDims(x, w; stride=stride, padding=padding, dilation=dilation)
#     y = similar(x, NNlib.output_size(cdims)..., NNlib.channels_out(cdims), size(x, 4))
#     return CUDNN.conv!(y, x, w, cdims)
# end

# function ∇conv2d_w(dy::CuArray{T,4}, x::CuArray{T,4}, w::CuArray{T,4}; stride=1, padding=0, dilation=1) where T
#     cdims = DenseConvDims(x, w; stride=stride, padding=padding, dilation=dilation)
#     dw = similar(w)
#     return CUDNN.∇conv_filter!(dw, x, dy, cdims)
# end

# function ∇conv2d_x(dy::CuArray{T,4}, x::CuArray{T,4}, w::CuArray{T,4}; stride=1, padding=0, dilation=1) where T
#     cdims = DenseConvDims(x, w; stride=stride, padding=padding, dilation=dilation)
#     dx = similar(x)
#     return CUDNN.∇conv_data!(dx, dy, w, cdims)
# end

# ## maxpool2d

# function maxpool2d(x::CuArray, kernel_size; stride=kernel_size, padding=0, dilation=1)
#     pdims = PoolDims(x, kernel_size; stride=stride, padding=padding, dilation=dilation)
#     y = similar(x, NNlib.output_size(pdims)..., NNlib.channels_out(pdims), size(x, 4))
#     return CUDNN.maxpool!(y, x, pdims)
# end

# function ∇maxpool2d_x(dy::CuArray, y::CuArray, x::CuArray, kernel_size; stride=kernel_size, padding=0, dilation=1)
#     pdims = PoolDims(x, kernel_size; stride=stride, padding=padding, dilation=dilation)
#     dx = similar(x)
#     return NNlib.∇maxpool!(dx, dy, y, x, pdims)
# end

# ## activations

# culogistic(x) = one(x) / (one(x) + CUDA.exp(-x))
# CUDA.cufunction(::typeof(logistic)) = culogistic

# ∇culogistic(dy, x) = culogistic(x) * (one(x) - culogistic(x)) * dy
# CUDA.cufunction(::typeof(∇logistic)) = ∇culogistic

# cusoftplus(x) = CUDA.log(CUDA.exp(x) + one(x))
# CUDA.cufunction(::typeof(softplus)) = cusoftplus

# ∇cusoftplus(dy, x) = culogistic(x) * dy
# CUDA.cufunction(::typeof(∇softplus)) = ∇cusoftplus

# ∇cuelu(dy, x::Real, alpha) = ifelse(x >= 0, x/1, alpha * CUDA.exp(x))
# CUDA.cufunction(::typeof(∇elu)) = ∇cuelu

## batchnorm

function batchnorm_impl(gamma::CuArray, beta::CuArray, x::CuArray,
                        mu::CuArray, sigma2::CuArray, momentum::Real; eps, training)
    NNlibCUDA.batchnorm(gamma, beta, x, mu, sigma2, momentum; eps=eps, training=training)
end

function ∇batchnorm_impl(gamma::CuArray, beta::CuArray, x::CuArray, dy::CuArray,
                         mu::CuArray, sigma2::CuArray, momentum; eps, training)
    NNlibCUDA.∇batchnorm(gamma, beta, x, dy, mu, sigma2, momentum; eps=eps, training=training)
end


text:
function preprocess_points(
    pts,
    margin;
    ORIGIN=nothing,
    vy_min_plus_max_input=0
    )
    (vx_min, vx_max, vy_min, vy_max) = find_bounding_box(pts)
    # flip within the box
    sx(x) = x + (ORIGIN===nothing ? (margin-vx_min) : ORIGIN[1])
    sy(y) = y + (ORIGIN===nothing ? (margin-vy_min) : ORIGIN[2])
    # parameter may be inherited, when ORIGIN is given
    # (vy_min+vy_max) must be given
    vy_min_plus_max = (ORIGIN===nothing ? (vy_min+vy_max) : vy_min_plus_max_input)
    # transform the point into a region with padding
    invY(v) = (sx(v[1]),sy(vy_min_plus_max-v[2]))
    pts_refl = invY.(pts)
    origin = invY((0,0))
    svg_head = "<svg width=\"$(real2str(vx_max-vx_min+2margin))\" height=\"$(real2str(vy_max-vy_min+2margin))\">"
    svg_tail = "</svg>"
    return (svg_head,svg_tail), (vx_min, vx_max, vy_min, vy_max), origin, pts_refl
end

function plot_vectors(
    vectors;
    headsize=0.1,
    line_style=("blue",4.0,:solid),
    margin=20.0
    )
    # ---------------------
    ϵ0 = 1e-8.*vectors[1]
    B = [v.+ϵ0 for v ∈ vectors]
    # ---------------------
    # find bounding box and revert y-axis within bbox
    ( (svg_head, svg_tail),
      (vx_min, vx_max, vy_min, vy_max),
      origin, pts_refl ) = preprocess_points( B, margin )
    lines = [(origin,v) for v ∈ pts_refl]
    # ---------------------
    (dashed,dotted) = LineStyle[line_style[3]]
    arrows = svg_arrows( lines,
                         line_style[1], line_style[2];
                         headsize=headsize,
                         dashed=dashed, dotted=dotted, dasharray=nothing )
    # ---------------------
    return (svg_head * "\n" * arrows * "\n" * svg_tail)
end

#NOTE This function is not responsible for any coordinate transformations !!
function plot_graph(
    V::Vector{Tuple{Float64,Float64}},
    E::Vector{Tuple{Int64,Int64}};
    MARGIN = 20.0,
    EDGECOLOR="red", EDGEWIDTH=1.0, EDGESTYLE=:solid,
    VERTEXCOLOR="black", VERTEXRAD=3.0, VERTEXSTYLE=:dot,
    HIGHLIGHT=[]
    )

    lines_svg = []

    # edges
    EC = (typeof(EDGECOLOR)<:Vector) ? EDGECOLOR : fill(EDGECOLOR,length(E))
    EW = (typeof(EDGEWIDTH)<:Vector) ? EDGEWIDTH : fill(EDGEWIDTH,length(E))
    ET = (typeof(EDGESTYLE)<:Vector) ? EDGESTYLE : fill(EDGESTYLE,length(E))
    EINFO = collect(zip(EC,EW,ET))

    # find the distinct pairs of (color, width)
    edges = Dict(p=>Vector{Tuple}() for p ∈ unique(EINFO))
    for (c,w,t) ∈ keys(edges)
        for (l,e) ∈ enumerate(E)
            if EINFO[l] == (c,w,t)
                e1 = V[e[1]]
                e2 = V[e[2]]
                push!( edges[(c,w,t)], ([e1[1],e1[2]],[e2[1],e2[2]]) )
            end
        end
    end

    # draw lines
    for (c,w,t) ∈ keys(edges)
        push!( lines_svg, svg_lines( edges[(c,w,t)], c, w,
                                     dashed=LineStyle[t][1],
                                     dotted=LineStyle[t][2] ) ) # (dashed,dotted)
    end

    # vertices
    VC = (typeof(VERTEXCOLOR)<:Vector) ? VERTEXCOLOR : fill(VERTEXCOLOR,length(V))
    VR = (typeof(VERTEXRAD  )<:Vector) ? VERTEXRAD   : fill(VERTEXRAD,  length(V))
    VT = (typeof(VERTEXSTYLE)<:Vector) ? VERTEXSTYLE : fill(VERTEXSTYLE,length(V))
    VINFO = collect(zip(VC,VR,VT))

    # find the distinct pairs of (color, width)
    vertices = Dict(p=>Tuple[] for p ∈ unique(VINFO))
    for (c,r,t) ∈ keys(vertices) # for each combination of c(color) r(radius) t(type)
        for (l,v) ∈ enumerate(V)
            if VINFO[l] == (c,r,t)
                push!( vertices[(c,r,t)], (v[1],v[2]) )
            end
        end
    end
    for (c,r,t) ∈ keys(vertices)
        push!( lines_svg, CircFuncDict[t](vertices[(c,r,t)], r, c) )
    end

    # highlights
    minVR = 0.618*minimum(VR)
    for (pos,col) ∈ HIGHLIGHT
        push!( lines_svg, svg_points(pos,minVR,col) )
    end

    # finalize
    (vx_min, vx_max, vy_min, vy_max) = find_bounding_box(V)
    W = real2str(vx_max-vx_min+2MARGIN)
    H = real2str(vy_max-vy_min+2MARGIN)
    svg_head = "<svg width=\"$W\" height=\"$H\">"
    svg_tail = "</svg>"

    return (svg_head * "\n" * join(lines_svg,"\n") * "\n" * svg_tail)
end

max_stars_repo_name: algorithmx/PlotSVG

max_issues_repo_name: algorithmx/PlotSVG

max_forks_repo_name: algorithmx/PlotSVG


text:
function jacobian(m::RealBond)

    # Load in endogenous state and eq cond indices
    endo = augment_model_states(m.endogenous_states_unnormalized,
                                n_model_states_unnormalized(m))
    eq   = m.equilibrium_conditions

    # Load in parameters, steady-state parameters, and grids
    γ::Float64     = m[:γ].value
    abar::Float64  = m[:abar].value
    R::Float64     = m[:R].value
    ν::Float64     = m[:ν].value
    ρz::Float64    = m[:ρ_z].value
    ρmon::Float64  = m[:ρmon].value
    κ::Float64     = m[:κ].value
    phipi::Float64 = m[:phipi].value
    aborrow::Float64 = abar/R

    ell::Vector{Float64}  = m[:lstar].value
    c::Vector{Float64}    = m[:cstar].value
    μ::Vector{Float64}    = m[:μstar].value
    η::Vector{Float64}    = m[:ηstar].value
    χss::Vector{Float64}  = m[:χstar].value
    β::Float64            = m[:βstar].value

    xgrid::Vector{Float64} = m.grids[:xgrid].points
    xwts::Vector{Float64}  = m.grids[:xgrid].weights
    sgrid::Vector{Float64} = m.grids[:sgrid].points
    swts::Vector{Float64}  = m.grids[:sgrid].weights
    ggrid::Vector{Float64} = m.grids[:ggrid]
    xgrid_total::Vector{Float64} = m.grids[:xgrid_total]
    sgrid_total::Vector{Float64} = m.grids[:sgrid_total]
    weights_total::Vector{Float64} = m.grids[:weights_total]

    elo::Float64 = get_setting(m, :elo)
    ehi::Float64 = get_setting(m, :ehi)

    nx::Int = get_setting(m, :nx)
    ns::Int = get_setting(m, :ns)

    # Mollifiers
    qfunction(x) = mollifier_realbond(x,ehi,elo)
    qp(x) = dmollifier_realbond(x,ehi,elo)

    # Jacobian stuff
    chipW, chipR, chipX = construct_chip_realbond(xgrid_total, γ, ν, aborrow, abar, R, χss)

    # Create auxiliary variables
    # c = min.(ell.^(-1.0/γ),χss) # Consumption Decision Rule
    # η = (sgrid_total.^(1.0/ν)).*(c.^(-γ/ν))

    # EE
    dF1_dELLP, dF1_dRRP, dF1_dWWP, dF1_dTTP, dF1_dELL, dF1_dWW, dF1_dRR, dF1_dTT, unc, a =
    euler_equation_realbond(nx, ns, qp, qfunction, xgrid, sgrid, ggrid, sgrid_total, xwts, swts,
                            chipR, chipW, chipX, R, γ, β, ν,
                            c, η, ell, χss)

    # KF
    dF2_dMU, dF2_dELL, dF2_dRR, dF2_dWW, dF2_dTT =
    kolmogorov_fwd_realbond(nx, ns, qfunction, qp, xgrid, sgrid, ggrid, xwts, swts,
                            unc, a, R, ν, γ, μ, η, c)

    # MKT clearing
    dF3_dMU, dF3_dELL, dF3_dWW, dF3_dRR, dF3_dTT, GDP =
    market_clearing_realbond(nx, ns, xgrid, sgrid, xwts, swts,
                             η, c, μ, ell, γ, ν, chipW, chipR, chipX, unc)

    # Make the Jacobian
    JJ = zeros(2*nx*ns+7,4*nx*ns+14)

    # Euler equation
    JJ[eq[:eq_euler], endo[:l′_t]]  = dF1_dELLP
    JJ[eq[:eq_euler], endo[:R′_t]]  = dF1_dRRP
    JJ[eq[:eq_euler], endo[:w′_t]]  = dF1_dWWP
    JJ[eq[:eq_euler], endo[:t′_t]]  = dF1_dTTP
    JJ[eq[:eq_euler], endo[:l_t]]   = dF1_dELL
    JJ[eq[:eq_euler], endo[:R_t]]   = dF1_dRR
    JJ[eq[:eq_euler], endo[:w_t]]   = dF1_dWW
    JJ[eq[:eq_euler], endo[:t_t]]   = dF1_dTT

    # KF
    JJ[eq[:eq_kolmogorov_fwd], endo[:μ′_t]] = -eye(nx*ns)
    JJ[eq[:eq_kolmogorov_fwd], endo[:μ_t]]  = dF2_dMU
    JJ[eq[:eq_kolmogorov_fwd], endo[:l_t]]  = dF2_dELL
    JJ[eq[:eq_kolmogorov_fwd], endo[:w_t]]  = dF2_dWW
    JJ[eq[:eq_kolmogorov_fwd], endo[:R_t]]  = dF2_dRR
    JJ[eq[:eq_kolmogorov_fwd], endo[:t_t]]  = dF2_dTT

    # mkt ckr
    JJ[eq[:eq_market_clearing], endo[:μ_t]]  = dF3_dMU
    JJ[eq[:eq_market_clearing], endo[:z_t]]  .= GDP
    JJ[eq[:eq_market_clearing], endo[:l_t]]  = dF3_dELL
    JJ[eq[:eq_market_clearing], endo[:w_t]]  .= dF3_dWW
    JJ[eq[:eq_market_clearing], endo[:R_t]]  .= dF3_dRR
    JJ[eq[:eq_market_clearing], endo[:t_t]]  .= dF3_dTT

    # TFP
    JJ[eq[:eq_TFP], endo[:z′_t]]  .= 1.0
    JJ[eq[:eq_TFP], endo[:z_t]]   .= -ρz

    # Phillips
    JJ[eq[:eq_phillips], endo[:π′_t]] .= -1.0/R
    JJ[eq[:eq_phillips], endo[:z_t]]  .= κ
    JJ[eq[:eq_phillips], endo[:w_t]]  .= -κ
    JJ[eq[:eq_phillips], endo[:π_t]]  .= 1.0

    # taylor
    JJ[eq[:eq_taylor], endo[:i_t]]   .= 1.0
    JJ[eq[:eq_taylor], endo[:π_t]]   .= -phipi*R
    JJ[eq[:eq_taylor], endo[:mon_t]] .= -1.0

    # fisher
    JJ[eq[:eq_fisher], endo[:i_t]]  .= 1.0
    JJ[eq[:eq_fisher], endo[:π′_t]] .= -R
    JJ[eq[:eq_fisher], endo[:R_t]]  .= -1.0

    # transfers
    JJ[eq[:eq_transfers], endo[:t_t]]  .= 1.0
    JJ[eq[:eq_transfers], endo[:z_t]]  .= -GDP
    JJ[eq[:eq_transfers], endo[:w_t]]  .= GDP

    # MP
    JJ[eq[:eq_monetary_policy], endo[:mon′_t]] .= 1.0
    JJ[eq[:eq_monetary_policy], endo[:mon_t]]  .= -ρmon

    if !m.testing && get_setting(m, :normalize_distr_variables)
        JJ = normalize(m, JJ)
    end

    return JJ
end

function construct_chip_realbond(xgrid_total::Vector{Float64},
                                 γ::Float64, ν::Float64, aborrow::Float64,
                                 abar::Float64, R::Float64,
                                 χss::Vector{Float64})
    chipW = (1.0+ν)./(ν./(aborrow .+ χss - xgrid_total) + γ./χss)
    chipR = (1.0/(R^2))*(ν*abar./(aborrow .+ χss - xgrid_total))./(ν./(aborrow .+ χss - xgrid_total) + γ./χss)
    chipX = (ν*χss)./(ν*χss + γ*(aborrow .+ χss - xgrid_total))

    return chipW, chipR, chipX
end

function euler_equation_realbond(nx::Int, ns::Int,
                                 qp::Function, qfunction::Function,
                                 xgrid::Vector{Float64}, sgrid::Vector{Float64},
                                 ggrid::Vector{Float64}, sgrid_total::Vector{Float64},
                                 xwts::Vector{Float64}, swts::Vector{Float64},
                                 chipR::Vector{Float64}, chipW::Vector{Float64}, chipX::Vector{Float64},
                                 R::Float64, γ::Float64, β::Float64, ν::Float64,
                                 c::Vector{Float64}, η::Vector{Float64}, ell::Vector{Float64}, χss::Vector{Float64})
    ξ   = zeros(nx*ns,nx*ns)
    Γ   = zeros(nx*ns)
    Ξ   = zeros(nx*ns,nx*ns)
    dF1_dELLP = zeros(nx*ns,nx*ns)
    dF1_dRRP  = zeros(nx*ns)
    dF1_dWWP  = zeros(nx*ns)
    dF1_dTTP  = zeros(nx*ns)
    dF1_dELL  = zeros(nx*ns,nx*ns)
    dF1_dWW   = zeros(nx*ns)
    dF1_dRR   = zeros(nx*ns)
    dF1_dTT   = zeros(nx*ns)
    unc = zeros(nx*ns)
    a   = zeros(nx*ns)

    # EE jacobian terms
    for ix=1:nx
        for is =1:ns
            i = ix + nx*(is-1)
            sumΓ   = 0.0
            sumRRP = 0.0
            sumWWP = 0.0
            sumTTP = 0.0
            sumWW  = 0.0
            sumTT  = 0.0
            sumξ   = 0.0
            sumELL = 0.0
            a[i] = R*(sgrid[is]*η[i]+xgrid[ix] - c[i])
            for ixp=1:nx
                for isp = 1:ns
                    ip = ixp + nx*(isp-1)
                    unc[ip] = ((ell[ip]^(-1/γ))<=χss[ip]) # =1 if unconstrained
                    mf = R*(sgrid[is]*η[i]+xgrid[ix] - c[i]) - xgrid[ixp]
                    ξ[i,ip] = β*R*(c[ip]^(-γ))*qp(mf)*ggrid[isp]*swts[isp]*xwts[ixp]
                    Ξ[i,ip] = β*R*(c[ip]^(-γ-1))*qfunction(mf)*ggrid[isp]*swts[isp]*xwts[ixp]
                    dF1_dELLP[i,ip] = Ξ[i,ip]*(ell[ip]^(-1- 1/γ))*unc[ip]
                    sumRRP += -Ξ[i,ip]*chipR[ip]*unc[ip]
                    sumWWP += -Ξ[i,ip]*chipW[ip]*unc[ip]
                    sumTTP += -Ξ[i,ip]*chipX[ip]*unc[ip]
                    sumELL += ξ[i,ip]*R*unc[i]*((1/ν) + sgrid[is]*η[i]/(ν*c[i]))*ell[i]^(-(1/γ) - 1)
                    sumWW += ξ[i,ip]*R*unc[i]*(1.0+1.0/ν)*sgrid[is]*η[i]
                    sumTT += ξ[i,ip]*R*unc[i]
                    sumξ += ξ[i,ip]
                end
            end
            Γ[i] = sumΓ
            dF1_dELL[i,i] = -1.0 + sumELL
            dF1_dRRP[i] = sumRRP
            dF1_dWWP[i] = sumWWP
            dF1_dTTP[i] = sumTTP
            dF1_dRR[i] = (ell[i] + sumξ*a[i])/R
            dF1_dWW[i] = sumξ*R*unc[i]*(1.0+1.0/ν)*sgrid_total[i]*η[i]
            dF1_dTT[i] = sumξ*R*unc[i]
        end
    end
    return dF1_dELLP, dF1_dRRP, dF1_dWWP, dF1_dTTP, dF1_dELL, dF1_dWW, dF1_dRR, dF1_dTT, unc, a
end

function kolmogorov_fwd_realbond(nx::Int, ns::Int,
                                 qfunction::Function, qp::Function,
                                 xgrid::Vector{Float64}, sgrid::Vector{Float64},
                                 ggrid::Vector{Float64}, xwts::Vector{Float64}, swts::Vector{Float64},
                                 unc::Vector{Float64}, a::Vector{Float64},
                                 R::Float64, ν::Float64, γ::Float64,
                                 μ::Vector{Float64}, η::Vector{Float64}, c::Vector{Float64})
    dF2_dMU   = zeros(nx*ns,nx*ns)
    dF2_dELL  = zeros(nx*ns,nx*ns)
    dF2_dRR = zeros(nx*ns)
    dF2_dWW = zeros(nx*ns)
    dF2_dTT = zeros(nx*ns)

    # KF jacobian terms
    for ixp = 1:nx
        for isp = 1:ns
            ip = ixp + nx*(isp-1)
            sumRR = 0.0
            sumWW = 0.0
            sumTT = 0.0
            for ix = 1:nx
                for is = 1:ns
                    i = ix + nx*(is-1)
                    mf = R*(sgrid[is]*η[i]+xgrid[ix] - c[i]) - xgrid[ixp]
                    dF2_dMU[ip,i]  = qfunction(mf)*ggrid[isp]*xwts[ix]*swts[is]*swts[isp]
                    dF2_dELL[ip,i] = qp(mf)*ggrid[isp]*xwts[ix]*swts[is]*swts[isp]*μ[ix]*ggrid[isp]*R*unc[i]*((1/γ) + sgrid[is]*η[i]/(ν*c[i]))
                    sumRR += qp(mf)*μ[ix]*ggrid[isp]*xwts[ix]*swts[is]*swts[isp]*unc[i]*a[i]/R
                    sumWW += qp(mf)*μ[ix]*ggrid[isp]*xwts[ix]*swts[is]*swts[isp]*unc[i]*R*(1+(1/ν))*sgrid[is]*η[i]
                    sumTT += qp(mf)*μ[ix]*ggrid[isp]*xwts[ix]*swts[is]*swts[isp]*unc[i]*R
                end
            end
            dF2_dRR[ip] = sumRR
            dF2_dWW[ip] = sumWW
            dF2_dTT[ip] = sumTT
        end
    end
    return dF2_dMU, dF2_dELL, dF2_dRR, dF2_dWW, dF2_dTT
end

function market_clearing_realbond(nx::Int, ns::Int,
                                  xgrid::Vector{Float64}, sgrid::Vector{Float64},
                                  xwts::Vector{Float64}, swts::Vector{Float64},
                                  η::Vector{Float64}, c::Vector{Float64},
                                  μ::Vector{Float64}, ell::Vector{Float64},
                                  γ::Float64, ν::Float64,
                                  chipW::Vector{Float64}, chipR::Vector{Float64},
                                  chipX::Vector{Float64}, unc::Vector{Float64})
    # MKT clearing
    dF3_dMU  = zeros(nx*ns)
    dF3_dELL = zeros(nx*ns)
    dF3_dWW  = 0.0
    dF3_dRR  = 0.0
    dF3_dTT  = 0.0
    GDP      = 0.0

    for ix =1:nx
        for is=1:ns
            i = ix + nx*(is-1)
            dF3_dMU[i]  = xwts[ix]*swts[is]*(sgrid[is]*η[i] - c[i])
            dF3_dELL[i] = xwts[ix]*swts[is]*μ[ix]*(1+sgrid[is]*γ*η[i]/(ν*c[i]))*(1/γ)*unc[i]*ell[i]^(-(1/γ)-1)
            dF3_dWW += swts[is]*xwts[ix]*μ[i]*((sgrid[is]*η[i])/ν - (1+sgrid[is]*γ*η[i]/(ν*c[i]))*(1-unc[i])*chipW[i])
            dF3_dRR += -swts[is]*xwts[ix]*μ[i]*(1+sgrid[is]*γ*η[i]/(ν*c[i]))*(1-unc[i])*chipR[i]
            dF3_dTT += -swts[is]*xwts[ix]*μ[i]*(1+sgrid[is]*γ*η[i]/(ν*c[i]))*(1-unc[i])*chipX[i]
            GDP += swts[is]*xwts[ix]*sgrid[is]*η[i]*μ[i]
        end
    end

    return dF3_dMU, dF3_dELL, dF3_dWW, dF3_dRR, dF3_dTT, GDP
end

function normalize(m::RealBond, JJ::Matrix{Float64})

    Qx, _, Qleft, Qright = compose_normalization_matrices(m)

    m <= Setting(:n_predetermined_variables, size(Qx, 1))

	Jac1 = Qleft*JJ*Qright

    return Jac1
end

function compose_normalization_matrices(m::RealBond)
    nx::Int = get_setting(m, :nx)
    ns::Int = get_setting(m, :ns)

    # Create PPP matrix
    P1 = kron(eye(ns),ones(nx,1))
    P2 = kron(eye(ns), eye(nx)[:, 2:end])
    P  = hcat(P1, P2)

    (Q,R) = qr(P)

    S         = Q[:, ns+1:end]'
    Qleft     = cat(eye(nx*ns),S,[1],[1],[1],[1],[1],[1],[1], dims = [1 2])
    Qx        = cat(S,[1],[1], dims = [1 2])
    Qy        = cat(eye(nx*ns),[1],[1],[1],[1],[1], dims = [1 2])
    Qright    = cat(Qx',Qy',Qx',Qy', dims = [1 2])

    return Qx, Qy, Qleft, Qright
end