text:
abstract type GeneralizedKroneckerProduct{T} <: AbstractMatrix{T} end

Base.eltype(K::GeneralizedKroneckerProduct{T}) where {T} = T

abstract type AbstractKroneckerProduct{T} <: GeneralizedKroneckerProduct{T} end

Base.IndexStyle(::Type{<:GeneralizedKroneckerProduct}) = IndexCartesian()

"""
    collect!(C::AbstractMatrix, K::GeneralizedKroneckerProduct)

In-place collection of `K` in `C`. If possible, specialized routines are used to
speed up the computation. The fallback is an element-wise iteration. In this case,
this function might be slow.
"""
function collect!(C::AbstractMatrix, K::GeneralizedKroneckerProduct)
    size(C) == size(K) || throw(DimensionMismatch("`K` $(size(K)) cannot be collected in `C` $(size(C))"))
    @inbounds for I in CartesianIndices(K)
        C[I] = K[I]
    end
    return C
end

"""
    collect(K::GeneralizedKroneckerProduct)

Collects a lazy instance of the `GeneralizedKroneckerProduct` type into a dense,
native matrix. Falls back to the element-wise case when not specialized method
is defined.
"""
function collect(K::GeneralizedKroneckerProduct{T}) where {T}
    C = Matrix{T}(undef, size(K)...)
    collect!(C, K)
    return C
end

"""
    KroneckerProduct{T,TA<:AbstractMatrix, TB<:AbstractMatrix} <: AbstractKroneckerProduct{T}

Concrete Kronecker product between two matrices `A` and `B`.
"""
struct KroneckerProduct{T<:Any, TA<:AbstractMatrix, TB<:AbstractMatrix} <: AbstractKroneckerProduct{T}
    A::TA
    B::TB
    function KroneckerProduct(A::AbstractMatrix{T}, B::AbstractMatrix{V}) where {T, V}
        return new{promote_type(T, V), typeof(A), typeof(B)}(A, B)
    end
end

# conversion
Base.convert(::Type{T}, K::KroneckerProduct) where {T<:KroneckerProduct} = K isa T ? K : T(K)
function KroneckerProduct{T,TA,TB}(K::KroneckerProduct) where {T,TA,TB}
    A, B = getmatrices(K)
    KroneckerProduct(convert(TA, A), convert(TB, B))
end

"""
    kronecker(A::AbstractMatrix, B::AbstractMatrix)

Construct a Kronecker product object between two arrays. Does not evaluate the
Kronecker product explictly.
"""
kronecker(A::AbstractMatrix, B::AbstractMatrix) = KroneckerProduct(A, B)

# version that have a vector as input reshape to matrices
kronecker(A::AbstractVector, B::AbstractMatrix) = KroneckerProduct(reshape(A,:,1), B)
kronecker(A::AbstractMatrix, B::AbstractVector) = KroneckerProduct(A, reshape(B,:,1))
kronecker(A::AbstractVector, B::AbstractVector) = KroneckerProduct(reshape(A,:,1), reshape(B,:,1))

"""
    kronecker(A::AbstractMatrix, B::AbstractMatrix)

Higher-order Kronecker lazy kronecker product, e.g.
```
kronecker(A, B, C, D)
```
"""
kronecker(A::AbstractVecOrMat, B::AbstractVecOrMat...) = kronecker(A,
                                                            kronecker(B...))

"""
    ⊗(A::AbstractMatrix, B::AbstractMatrix)

Binary operator for `kronecker`, computes as Lazy Kronecker product. See
`kronecker` for documentation.
"""
⊗(A::AbstractVecOrMat, B::AbstractVecOrMat) = kronecker(A, B)
⊗(A::AbstractVecOrMat...) = kronecker(A...)

"""
    getindex(K::AbstractKroneckerProduct, i1::Integer, i2::Integer)

Computes and returns the (i,j)-th element of an `AbstractKroneckerProduct` K.
Uses recursion if `K` is of an order greater than two.
"""
function getindex(K::AbstractKroneckerProduct, i1::Integer, i2::Integer)
    A, B = getmatrices(K)
    k, l = size(B)
    return (A[cld(i1, k), cld(i2, l)]::eltype(A)) * (B[(i1 - 1) % k + 1, (i2 - 1) % l + 1]::eltype(B))
end

"""
    getallfactors(K::AbstractKroneckerProduct)

Obtain all factors in an `AbstractKroneckerProduct` object.
"""
getallfactors(K::AbstractKroneckerProduct) = (getallfactors(K.A)..., getallfactors(K.B)...)
getallfactors(K::AbstractMatrix) = (K,)

"""
    getmatrices(K::AbstractKroneckerProduct)

Obtain the two matrices of an `AbstractKroneckerProduct` object.
"""
getmatrices(K::AbstractKroneckerProduct) = (K.A, K.B)

"""
    getmatrices(A::AbstractArray)

Returns a matrix itself. Needed for recursion.
"""
getmatrices(A::AbstractArray) = (A,)

firstmatrix(K::AbstractKroneckerProduct) = first(getmatrices(K))
lastmatrix(K::AbstractKroneckerProduct) = last(getmatrices(K))

"""
    size(K::AbstractKroneckerProduct)

Returns a the size of an `AbstractKroneckerProduct` instance.
"""
function size(K::AbstractKroneckerProduct)
    A, B = getmatrices(K)
    (m, n) = size(A)
    (k, l) = size(B)
    return m * k, n * l
end

"""
    size(K::GeneralizedKroneckerProduct)

Returns a the size of an `GeneralizedKroneckerProduct` instance.
"""
size(K::GeneralizedKroneckerProduct, dim::Integer) = size(K)[dim]

# CHECKS

"""
    issquare(A::AbstractMatrix)

Checks if an array is a square matrix.
"""
function issquare(A::AbstractMatrix)
    m, n = size(A)
    return m == n
end

"""
    issquare(A::Factorization)

Checks if a Factorization struct represents a square matrix.
"""
function issquare(A::Factorization)
    m, n = size(A)
    return m == n
end

"""
    issymmetric(K::AbstractKroneckerProduct)

Checks if a Kronecker product is symmetric.
"""
function LinearAlgebra.issymmetric(K::AbstractKroneckerProduct)
    A, B = getmatrices(K)
    return issymmetric(A) && issymmetric(B)
end

"""
    ishermitian(K::AbstractKroneckerProduct)

Checks if a Kronecker product is Hermitian.
"""
function LinearAlgebra.ishermitian(K::AbstractKroneckerProduct)
    A, B = getmatrices(K)
    return ishermitian(A) && ishermitian(B)
end

"""
    isposdef(K::AbstractKroneckerProduct)

Test whether a Kronecker product is positive definite (and Hermitian) by trying
to perform a Cholesky factorization of K.
"""
function LinearAlgebra.isposdef(K::AbstractKroneckerProduct)
    A, B = getmatrices(K)
    return isposdef(A) && isposdef(B)
end

# SCALAR PROPERTIES

"""
    order(M::AbstractMatrix)

Returns the order of a matrix, i.e. how many matrices are involved in the
Kronecker product (default to 1 for general matrices).
"""
order(M::AbstractMatrix) = 1
order(M::AbstractKroneckerProduct) = order(M.A) + order(M.B)

"""
    det(K::AbstractKroneckerProduct)

Compute the determinant of a Kronecker product.
"""
function LinearAlgebra.det(K::AbstractKroneckerProduct{T}) where {T}
    checksquare(K)
    A, B = getmatrices(K)
    if issquare(A) && issquare(B)
        m = size(A)[1]
        n = size(B)[1]
        return det(A)^n * det(B)^m
    else
        return zero(T)
    end
end

"""
    logdet(K::AbstractKroneckerProduct)

Compute the logarithm of the determinant of a Kronecker product.
"""
function LinearAlgebra.logdet(K::AbstractKroneckerProduct{T}) where {T}
    checksquare(K)
    A, B = getmatrices(K)
    if issquare(A) && issquare(B)
        m = size(A)[1]
        n = size(B)[1]
        return n * logdet(A) + m * logdet(B)
    else
        return real(T)(-Inf)
    end
end

"""
    tr(K::AbstractKroneckerProduct)

Compute the trace of a Kronecker product.
"""
function LinearAlgebra.tr(K::AbstractKroneckerProduct)
    checksquare(K)
    A, B = getmatrices(K)
    if issquare(A) && issquare(B)
        return tr(A) * tr(B)
    else
        return sum(diag(K))
    end
end

"""
    inv(K::AbstractKroneckerProduct)

Compute the inverse of a Kronecker product.
"""
function inv(K::AbstractKroneckerProduct)
    checksquare(K)
    A, B = getmatrices(K)
    if issquare(A) && issquare(B)
        return KroneckerProduct(inv(A), inv(B))
    else
        throw(SingularException(1))
    end
end

"""
    pinv(K::AbstractKroneckerProduct)

Compute the Moore-Penrose pseudo-inverse of a Kronecker product.
"""
function LinearAlgebra.pinv(K::AbstractKroneckerProduct)
    A, B = getmatrices(K)
    return KroneckerProduct(pinv(A), pinv(B))
end

"""
    adjoint(K::AbstractKroneckerProduct)

Compute the adjoint of a Kronecker product.
"""
function adjoint(K::AbstractKroneckerProduct)
    A, B = getmatrices(K)
    return kronecker(A', B')
end

"""
    transpose(K::AbstractKroneckerProduct)

Compute the transpose of a Kronecker product.
"""
function Base.transpose(K::AbstractKroneckerProduct)
    A, B = getmatrices(K)
    return kronecker(transpose(A), transpose(B))
end

function Base.conj(K::AbstractKroneckerProduct)
    A, B = getmatrices(K)
    return kronecker(conj(A), conj(B))
end

function LinearAlgebra.diag(K::KroneckerProduct)
    A, B = getmatrices(K)
    if issquare(A) && issquare(B)
        return kron(diag(K.A), diag(K.B))
    end
    return K[diagind(K)]
end

# COLLECTING

#=
"""
    collect(K::AbstractKroneckerProduct)

Collects a lazy instance of the `AbstractKroneckerProduct` type into a full,
native matrix. Equivalent with `Matrix(K::AbstractKroneckerProduct)`.
"""
function collect(K::AbstractKroneckerProduct)
    A, B = getmatrices(K)
    return kron(A, B)
end
=#

_maybecollect(A::GeneralizedKroneckerProduct) = collect(A)
_maybecollect(A::AbstractArray) = A

# function for in-place Kronecker product
function _kron!(C::AbstractArray, A::AbstractArray, B::AbstractArray)
    m = first(LinearIndices(C)) - 1
    @inbounds for j = axes(A,2), l = axes(B,2), i = axes(A,1)
        Aij = A[i,j]
        for k = axes(B,1)
            C[m += 1] = Aij * B[k,l]
        end
    end
    return C
end

_kron!(C::AbstractArray, A::GeneralizedKroneckerProduct, B::AbstractArray) = _kron!(C, collect(A), B)
_kron!(C::AbstractArray, A::AbstractArray, B::GeneralizedKroneckerProduct) = _kron!(C, A, collect(B))
_kron!(C::AbstractArray, A::GeneralizedKroneckerProduct, B::GeneralizedKroneckerProduct) = _kron!(C, collect(A), collect(B))

@inline function _kron!(C::AbstractArray, As::Tuple{AbstractArray, Vararg{AbstractArray}}, Bs::Tuple{AbstractArray, Vararg{AbstractArray}}, f = identity)
    m = first(LinearIndices(C)) - 1
    A1 = first(As)
    B1 = first(Bs)
    for j = axes(A1,2), l = axes(B1,2), i = axes(A1,1)
        Aijs = map(A -> A[i,j], As)
        for k = 1:size(Bs[1],1)
            Bkls = map(B -> B[k,l], Bs)
            Aijs_times_Bkls = map(*, Aijs, Bkls)
            C[m += 1] = f(Aijs_times_Bkls...)
        end
    end
    return C
end

"""
    collect!(C::AbstractMatrix, K::AbstractKroneckerProduct)

In-place collection of `K` in `C` where `K` is an `AbstractKroneckerProduct`, i.e.,
`K = A ⊗ B`. This is equivalent to the broadcasted assignment `C .= K`.

    collect!(f, C::AbstractMatrix, K1::AbstractKroneckerProduct, Ks::AbstractKroneckerProduct...)

Evaluate `f.(K1, Ks...)` and assign it in-place to `C`. This is equivalent to the broadcasted
operation `C .= f.(K1, Ks...)`.
"""
function collect!(C::AbstractMatrix, K::AbstractKroneckerProduct)
    size(C) == size(K) || throw(DimensionMismatch("`K` $(size(K)) cannot be collected in `C` $(size(C))"))
    A, B = getmatrices(K)
    return _kron!(C, A, B)
end

@inline function collect!(f, C::AbstractMatrix, K1::AbstractKroneckerProduct, Ks::AbstractKroneckerProduct...)
    @noinline throwdm(K1sz, Csz) = throw(DimensionMismatch("`K` $K1sz cannot be collected in `C` $Csz"))
    size(C) == size(K1) || throwdm(size(K1), size(C))
    for K in Ks
        size(C) == size(K) || throwdm(size(K), size(C))
    end
    Ks_all = (K1, Ks...)
    As = map(x -> first(getmatrices(x)), Ks_all)
    Bs = map(x -> last(getmatrices(x)), Ks_all)
    return _kron!(C, As, Bs, f)
end

"""
    Matrix(K::GeneralizedKroneckerProduct)

Converts a `GeneralizedKroneckerProduct` instance to a Matrix type.
"""
Base.Matrix(K::GeneralizedKroneckerProduct) = collect(K)

function Base.:+(A::AbstractKroneckerProduct, B::StridedMatrix)
    T = promote_type(eltype(A), eltype(B))
    C = similar(Array{T}, size(A))
    C .= A
    C .+= B
    return C
end
Base.:+(A::StridedMatrix, B::AbstractKroneckerProduct) = B + A

function Base.:+(A::AbstractKroneckerProduct, B1::AbstractKroneckerProduct, Brest::AbstractKroneckerProduct...)
    Bs = (B1, Brest...)
    for B in Bs
        Base.promote_shape(A, B) # check size compatibility
    end
    # special methods to handle kronecker products with singleton matrices
    # if one matrix is common to all products, we only need to add the other matrix
    if all(x -> firstmatrix(x) === firstmatrix(A), Bs)
        K1 = kronecker(firstmatrix(A), +(lastmatrix(A), map(lastmatrix, Bs)...))
        return collect(K1)
    elseif all(x -> lastmatrix(x) === lastmatrix(A), Bs)
        K2 = kronecker(+(firstmatrix(A), map(firstmatrix, Bs)...), lastmatrix(A))
        return collect(K2)
    end
    # if the sizes of the component matrices are compatible, the operation may be
    # short-circuited
    sa = map(size, getmatrices(A))
    if all(x -> map(size, getmatrices(x)) == sa, Bs)
        return broadcast(+, A, Bs...)
    end
    # collect the arrrays before adding to avoid indexing into the kronecker products
    return +(collect(A), map(collect, Bs)...)
end

function Base.:-(A::AbstractKroneckerProduct, B::StridedMatrix)
    T = promote_type(eltype(A), eltype(B))
    C = similar(Array{T}, size(A))
    C .= A
    C .-= B
    return C
end
function Base.:-(A::StridedMatrix, B::AbstractKroneckerProduct)
    T = promote_type(eltype(A), eltype(B))
    C = similar(Array{T}, size(A))
    @. C = -B
    C .+= A
    return C
end

const KronProdDiagonal = KroneckerProduct{<:Any, <:Diagonal, <:Diagonal}
for T in [:Diagonal, :UniformScaling]
    @eval Base.:+(K::KronProdDiagonal, D::$T) = Diagonal(K) + D
    @eval Base.:+(D::$T, K::KronProdDiagonal) = D + Diagonal(K)
    @eval Base.:-(K::KronProdDiagonal, D::$T) = Diagonal(K) - D
    @eval Base.:-(D::$T, K::KronProdDiagonal) = D - Diagonal(K)
end

function Base.:-(A::AbstractKroneckerProduct, B::AbstractKroneckerProduct)
    # special methods to handle kronecker products with singleton matrices
    # if one matrix is common to all products, we only need to add the other matrix
    if firstmatrix(B) === firstmatrix(A)
        K1 = kronecker(firstmatrix(A), lastmatrix(A) - lastmatrix(B))
        return collect(K1)
    elseif lastmatrix(B) === lastmatrix(A)
        K2 = kronecker(firstmatrix(A) - firstmatrix(B), lastmatrix(A))
        return collect(K2)
    end
    # if the sizes of the component matrices are compatible, the operation may be
    # short-circuited
    if map(size, getmatrices(A)) == map(size, getmatrices(B))
        return A .- B
    end
    # collect the arrrays before subtracting to avoid indexing into the kronecker products
    return collect(A) - collect(B)
end

function Base.kron(K::AbstractKroneckerProduct, C::AbstractMatrix)
    A, B = getmatrices(K)
    return kron(kron(A, B), C)
end

function Base.kron(A::AbstractMatrix, K::AbstractKroneckerProduct)
    B, C = getmatrices(K)
    return kron(A, kron(B, C))
end

Base.kron(K1::AbstractKroneckerProduct,
            K2::AbstractKroneckerProduct) = kron(collect(K1), collect(K2))

# mixed-product property
function Base.:*(K1::AbstractKroneckerProduct, K2::AbstractKroneckerProduct)
    A, B = getmatrices(K1)
    C, D = getmatrices(K2)
    # check for size
    if size(A, 2) != size(C, 1)
        throw(DimensionMismatch("Mismatch between A and C in (A ⊗ B)(C ⊗ D)"))
    end
    if size(B, 2) != size(D, 1)
        throw(DimensionMismatch("Mismatch between B and D in (A ⊗ B)(C ⊗ D)"))
    end
    return (A * C) ⊗ (B * D)
end

"""
    lmul!(a::Number, K::AbstractKroneckerProduct)

Scale an `AbstractKroneckerProduct` `K` inplace by a factor `a` by rescaling the
**left** matrix.
"""
function LinearAlgebra.lmul!(a::Number, K::AbstractKroneckerProduct)
    A, B = getmatrices(K)
    lmul!(a, A)
end

"""
    rmul!(K::AbstractKroneckerProduct, a::Number)

Scale an `AbstractKroneckerProduct` `K` inplace by a factor `a` by rescaling the
**right** matrix.
"""
function LinearAlgebra.rmul!(K::AbstractKroneckerProduct, a::Number)
    A, B = getmatrices(K)
    rmul!(B, a)
end

function Base.:*(a::Number, K::AbstractKroneckerProduct)
    A, B = getmatrices(K)
    kronecker(a * A, B)
end

function Base.:*(K::AbstractKroneckerProduct, a::Number)
    A, B = getmatrices(K)
    kronecker(A, B * a)
end

function LinearAlgebra.power_by_squaring(K::KroneckerProduct, p::Integer)
    A, B = getmatrices(K)
    kronecker(A^p, B^p)
end

function LinearAlgebra.svdvals(K::KroneckerProduct)
    A, B = getmatrices(K)
    σA = svdvals(A)
    σB = svdvals(B)
    σ = sort!(kron(σA, σB), rev = true)
    n = minimum(size(K))
    if length(σ) < n
        append!(σ, zeros(eltype(σ), n - length(σ)))
    end
    return σ
end

# Broadcasting machinery

Base.copyto!(dest::AbstractMatrix, K::AbstractKroneckerProduct) = collect!(dest, K)

struct AbsKronProdStyle <: Broadcast.AbstractArrayStyle{2} end
AbsKronProdStyle(::Val{N}) where {N} = Broadcast.DefaultArrayStyle{N}()
AbsKronProdStyle(::Val{2}) = AbsKronProdStyle()

Base.BroadcastStyle(::Type{<:AbstractKroneckerProduct}) = AbsKronProdStyle()

function Base.similar(bc::Broadcast.Broadcasted{AbsKronProdStyle}, ::Type{T}) where {T}
    similar(Array{T}, axes(bc))
end

@inline function Base.copyto!(dest::AbstractArray, bc::Broadcast.Broadcasted{AbsKronProdStyle})
    @noinline throwdm(axdest, axsrc) =
        throw(DimensionMismatch("destination axes $axdest are not compatible with source axes $axsrc"))
    axes(dest) == axes(bc) || throwdm(axes(dest), axes(bc))
    # Some common cases that may be short-circuited
    # Case 1, example: A .= B
    if bc.args isa Tuple{AbstractKroneckerProduct}
        A = bc.args[1]
        collect!(bc.f, dest, A)
        return dest
    # Case 2, example: 2 .* B
    elseif bc.args isa Tuple{Number, AbstractKroneckerProduct}
        A = last(bc.args)
        n = first(bc.args)
        collect!(let n = n; x -> bc.f(n, x); end, dest, A)
        return dest
    # Case 3, example: B .* 2
    elseif bc.args isa Tuple{AbstractKroneckerProduct, Number}
        A = first(bc.args)
        n = last(bc.args)
        collect!(let n = n; x -> bc.f(x, n); end, dest, A)
        return dest
    end
    # An operation like K1 .+ K2 may be short-circuited if the component matrices
    # have the same size.
    bcf = Broadcast.flatten(bc)
    if all(x -> x isa AbstractKroneckerProduct, bcf.args)
        sz1 = map(size, getmatrices(bcf.args[1]))
        if all(x -> map(size, getmatrices(x)) == sz1, bcf.args[2:end])
            collect!(bcf.f, dest, bcf.args...)
            return dest
        end
    elseif all(x -> x isa Union{AbstractKroneckerProduct, StridedArray}, bcf.args)
        # Performance is better if the kronecker products are collected before the
        # broadcasted operation is performed, although this incurs allocations
        broadcast!(bcf.f, dest, map(_maybecollect, bcf.args)...)
        return dest
    end
    # The general case that indexes into each array.
    # This may be slow as indexing into an AbstractKroneckerProduct is expensive
    bc′ = Broadcast.preprocess(dest, bc)
    # Performance may vary depending on whether `@inbounds` is placed outside the
    # for loop or not. (cf. https://github.com/JuliaLang/julia/issues/38086)
    @inbounds @simd for I in eachindex(bc′)
        dest[I] = bc′[I]
    end
    return dest
end


text:
# using ADCME
using PyTensorFlow
using PyCall
using LinearAlgebra
using PyPlot
using Random
Random.seed!(233)

py"""
import tensorflow as tf
libEikonal = tf.load_op_library('./build/libEikonal.so')
@tf.custom_gradient
def eikonal_lekj(f, srcx, srcy, m, n, h):
    u = libEikonal.eikonal(f, srcx, srcy, m, n, h)
    def grad(du):
        return libEikonal.eikonal_grad(du, u, f, srcx, srcy, m, n, h)
    return u, grad
"""
eikonal_ = py"eikonal_lekj"

function eikonal(f::Union{Array{Float64}, PyObject},
    srcx::Int64,srcy::Int64,h::Float64)
    n_, m_ = size(f) # m width, n depth
    n = n_-1
    m = m_-1
    # eikonal_ = load_op_and_grad("$(@__DIR__)/build/libEikonal","eikonal")
    # f,srcx,srcy,m,n,h = convert_to_tensor([f,srcx,srcy,m,n,h], [Float64,Int64,Int64,Int64,Int64,Float64])
    f = tf.cast(f, dtype=tf.float64)
    srcx = tf.cast(srcx, dtype=tf.int64)
    srcy = tf.cast(srcy, dtype=tf.int64)
    m = tf.cast(m, dtype=tf.int64)
    n = tf.cast(n, dtype=tf.int64)
    h = tf.cast(h, dtype=tf.float64)
    f = tf.reshape(f, (-1,))
    u = eikonal_(f,srcx,srcy,m,n,h)
    u.set_shape((length(f),))
    tf.reshape(u, (n_, m_))
end

# TODO: specify your input parameters
m = 60
n = 30
h = 0.1

f = ones(n+1, m+1)
for i = 1:m+1
    f[12:18, i] .= 10.
end
srcx = 30
srcy = 3
u = eikonal(f,srcx,srcy,h)
sess = Session(); init(sess)
@show run(sess, u)
pcolormesh(run(sess, u)|>Array)
axis("scaled")
colorbar()
gca().invert_yaxis()
# uncomment it for testing gradients
# error()

# TODO: change your test parameter to `m`
#       in the case of `multiple=true`, you also need to specify which component you are testings
# gradient check -- v
function scalar_function(m_)
    return sum(eikonal(m_,srcx,srcy,h)^2)
end

# TODO: change `m_` and `v_` to appropriate values
m_ = constant(rand((m+1),(n+1)))
v_ = rand((m+1),(n+1))
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)

ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 0.1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session(); init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")

ext: jl


text:
#!/usr/bin/env julia

using Pkg
Pkg.activate("../../")

using HTTP
using Gumbo
using AbstractTrees

request = HTTP.request("GET", "http://www.unicode.org/reports/tr14/tr14-34.html#PairBasedImplementation")
body = request.body
document = parsehtml(String(body))
pairTable = nothing
for element in PreOrderDFS(document.root)
  if isa(element, HTMLElement)
    if haskey(attrs(element), "class")
      #println(getattr(element, "class"))
      if getattr(element, "class") == "pair"
        global pairTable = element.children[1]
      end
    end
  end
end

function convertBreakOpportunity(opportunity::String)
  map = Dict(
    "^" => "Prohibited",
    "%" => "Indirect",
    "@" => "ProhibitedCM",
    "#" => "IndirectCM",
    "_" => "Direct"
  )

  return map[opportunity]
end

outfile = open("pairtable.h", "w")
println(outfile, """
    #ifndef PAIR_TABLE_H
    #define PAIR_TABLE_H
  """);
println(outfile, """
    typedef enum BreakOpportunity {
      Prohibited,
      Indirect,
      ProhibitedCM,
      IndirectCM,
      Direct
    } BreakOpportunity;
  """);
println(outfile, """
    #define PAIR_TABLE_SIZE 29
    static const BreakOpportunity pairTable[PAIR_TABLE_SIZE][PAIR_TABLE_SIZE] = {
  """);

for row in Iterators.drop(pairTable.children, 1)
  print(outfile, "    {")
  for cell in Iterators.drop(row.children, 1)
    element = cell
    while !isa(element, HTMLText)
      element = element.children[1]
    end
    print(outfile, string(convertBreakOpportunity(element.text)) * ", ")
  end
  println(outfile, "},")
end

println(outfile, """
    };
  """);

println(outfile, """
    #endif
  """);

close(outfile)

max_stars_repo_name: jonathanbuchanan/PhysicsEngine

max_issues_repo_name: jonathanbuchanan/PhysicsEngine

max_forks_repo_name: jonathanbuchanan/PhysicsEngine


text:
using Statistics

"""
Calculate the autocorrelation of `data` at the given distance `dist`
"""
function autocorrelation(data::Vector{<:Real}, dist::Integer)
    xᵢ = data[1:(length(data) - dist)]
    xᵢ₊ⱼ = data[(1 + dist):end]
    xᵢxᵢ₊ⱼavg = mean(xᵢ .* xᵢ₊ⱼ)
    return (xᵢxᵢ₊ⱼavg - mean(xᵢ) * mean(xᵢ₊ⱼ)) / var(data)
end

max_stars_repo_name: slhshamloo/comp-phys

max_issues_repo_name: slhshamloo/comp-phys

max_forks_repo_name: slhshamloo/comp-phys


text:
using StructuralCausalModels, Test

ProjDir = @__DIR__

amat_data = reshape([
  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,
  0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,
  0,0,0,0,1,0,1,0,1,1,0,0,0,0,0,0,
  1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,
  0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
  1,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,
  0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0
], (16,16));

vars = [Symbol("n$i") for i in 1:size(amat_data, 1)]
a = NamedArray(Int.(amat_data), (vars, vars), ("Rows", "Cols"));

dag = DAG("ag_example", a)

fn = joinpath(mktempdir(), "ag.dot")
to_graphviz(dag, fn)
#Sys.isapple() && run(`open -a GraphViz.app $(fn)`)

m = [:n3, :n5, :n6, :n15, :n16];
c = [:n4, :n7];

ag = ancestral_graph(dag; m=m, c=c)
#indx = indexin(names(ag, 1), vars)
indx = filter(x -> !isnothing(x), indexin(vars, names(ag, 1)))
ag = ag[indx, indx]

agmat = [
0 100 0  0  0  0   0  0 100;
100 0 0  0  0  0   0  0   0;
0   0 0  0  0  0   0  0   0;
0   0 1  0  0 10   0  0   0;
0   0 1  0  0  0   0  0   0;
0   1 0 10  1  0   0  0   0;
1   0 0  0  0  0   0  0 100;
0   0 0  0  1  0   1  0   0;
100 0 0  0  0  0 100  0   0;
]


text:
__precompile__()

module PostSPH

using Printf #To construct string message easily
using StaticArrays

#Add to project.toml manually: https://discourse.julialang.org/t/update-project-toml-manually/32477

include("ReadXML.jl")
include("SaveVTK.jl")

export Cat,
    readBi4Array, readBi4_NumberOfParticles, readBi4_CurrentTotalParticles, readBi4_Time

##Hardcoded enum - Cat is "category"
"""
    Cat are hardcoded enums to specify the potential arrays extractable from bi4 files.
"""
@enum Cat begin
    Points
    Idp
    Vel
    Rhop
end

##Read from binary files directly

## Make two dicts,   one for the enum -> string and one for the data type - Bi4 specific
const varNames = (:key, :offset)

const IdpSearch = transcode(UInt8, "ARRAY\x03\0\0\0Idp")[:]
const IdpOffset = 11
const IdpKey = NamedTuple{varNames}([IdpSearch, IdpOffset])

const PosSearch = transcode(UInt8, "ARRAY\x03\0\0\0Pos")[:]
const PosOffset = 11
const PosKey = NamedTuple{varNames}([PosSearch, PosOffset])

const VelSearch = transcode(UInt8, "ARRAY\x03\0\0\0Vel")[:]
const VelOffset = 11
const VelKey = NamedTuple{varNames}([VelSearch, VelOffset])

const RhopSearch = transcode(UInt8, "ARRAY\x04\0\0\0Rhop")[:]
const RhopOffset = 12
const RhopKey = NamedTuple{varNames}([RhopSearch, RhopOffset])

const searchKeyBi4 =
    Dict{Cat,NamedTuple}(Idp => IdpKey, Points => PosKey, Vel => VelKey, Rhop => RhopKey)
const catTypeBi4 =
    Dict{Cat,DataType}(Idp => Int32, Points => Float32, Vel => Float32, Rhop => Float32)
const catArrayBi4 = Dict{Cat,Int64}(Idp => 1, Points => 2, Vel => 2, Rhop => 1)
const catColBi4 = Dict{Cat,Int64}(Idp => 1, Points => 3, Vel => 3, Rhop => 1)

##Lists files in directory and only returns applicable files, ie. "Part_XXXX.bi4"
# first_file bug
function _dirFiles()
    files = readdir()
    #Operation on dirFiles instantly
    filter!(x -> occursin(r"Part_\d{4}.bi4", x), files)
    return files
end

#Command for plotting
#println(Char.(key))
#scatter(getindex.(a[1], 1), getindex.(a[1], 3))
#readBi4Array(typ::Cat,Bi4Files::Array{String,1}=_dirFiles()) = readBi4Array(typ, false, Bi4Files)
#readBi4Array(typ::Cat,Bi4Files::Array{String,1}) = readBi4Array(typ, false, Bi4Files)
readBi4Array(typ::Cat, Bi4Files::String) = readBi4Array(typ, [Bi4Files])
function readBi4Array(typ::Cat, Bi4Files::Array{String,1} = _dirFiles())

    nBi4 = size(Bi4Files)[1]

    key = searchKeyBi4[typ].key
    offset = searchKeyBi4[typ].offset
    ncol = catColBi4[typ]
    T = catTypeBi4[typ]

    # THIS BREAKS SAVE VTK
    j = Vector{Array{T,1}}(undef, nBi4)
    Threads.@threads for i = 1:nBi4
        j[i], ~ = _readBi4(Bi4Files[i], key, offset, T, ncol)
    end

    return j
end

function _readBi4(file::String, key, offset, T, ncol)

    # Import a full bi4 file as Array{UInt8,1}
    ft = open(file, read = true)
    rf = read(ft)
    close(ft)

    # Start position is found by search the file for the key and finding
    # first occurence, then adding offset
    startPos = Base._searchindex(rf, key, 1) + offset #1 byte offset

    #+1 due to hexeditor/julia?. "n id start" and "n id end"
    nid_s = startPos + 1 + sizeof(Int64)
    nid_e = nid_s - 1 + sizeof(Int32)

    # Every "Array" in bi4 mentions the number of particles, read it in
    n = reinterpret(Int32, rf[nid_s:nid_e])[1]

    #data id start
    # Multiply with 4 here since UInt8 size, times number of particles, times
    # times number of columns gives the correct indices in the rf array for
    # Float32, Int32 etc.
    did_s = nid_e + 1 + sizeof(Int32)
    did_e = did_s - 1 + 4 * n * ncol

    # Reinterpret the data as the specified data type, extract the relevant
    # snip of Array{UInt8,1} in "rf"
    data = reinterpret(T, rf[did_s:did_e])

    return data, n
end

function readBi4_NumberOfParticles(Bi4Files::Array{String,1} = _dirFiles())

    nBi4 = size(Bi4Files)[1]

    j = Vector{Array{Int32,1}}(undef, nBi4)

    ParticleString = ["CaseNp", "CaseNfixed", "CaseNmoving", "CaseNfloat", "CaseNfluid"]
    for i = 1:nBi4
        ft = open(Bi4Files[i], read = true)
        j_inner = zeros(Int32, length(ParticleString))
        for k in eachindex(j_inner)
            readuntil(ft, ParticleString[k])
            read(ft, Int32)
            j_inner[k] = read(ft, Int32)
        end
        j[i] = j_inner
        close(ft)
    end
    return ParticleString, j
end

#Npok is the current number of actual particles in the bi4 file
function readBi4_CurrentTotalParticles(Bi4Files::Array{String,1} = _dirFiles())

    nBi4 = size(Bi4Files)[1]

    j = Vector{Array{Int32,1}}(undef, nBi4)

    ParticleString = "Npok"
    j = zeros(Int32, (nBi4,))
    for i = 1:nBi4
        ft = open(Bi4Files[i], read = true)
        readuntil(ft, ParticleString)
        read(ft, Int32)
        j[i] = read(ft, Int32)
        close(ft)
    end
    return j
end

## Function to read the time at current simulation step, as in "XXXX.out"
function readBi4_Time(Bi4Files::Array{String,1} = _dirFiles())
    nBi4 = size(Bi4Files)[1]

    T = Float64
    j = Vector{T}(undef, nBi4)

    ParticleString = ["TimeStep"]
    Threads.@threads for i = 1:nBi4
        ft = open(Bi4Files[i], read = true)
        readuntil(ft, ParticleString[1])
        read(ft, Int32)
        j[i] = read(ft, T)
        close(ft)
    end
    return j
end

#Function to only find specific Idps
#for MovingSquare example "Bodies[2][1]""
#In the for loop the first index is the index of the relevant "typ" array,
#the second index is the sorting of this "typ" array corresponding to the "idp"
#array and "start:move" are the number of particles defined by the "Body", where
#0 and 1 indexing from C++ and Julia differences has been taken into account.
function readBi4Body(Body, typ)
    start = getfield(Body, :beg) + 1    #First idp

    idp_vec = readBi4Array(Idp)
    val_vec = readBi4Array(typ)

    nBi4 = length(idp_vec)

    k = []
    move = []

    if Body.bool == true
        move = readBi4Npok()
        k = collect(1:nBi4)
    else
        move = start + getfield(Body, :count) - 1  #Number of particles from first idp
        k = ones(Int, nBi4)
    end

    j = similar(val_vec)

    if typ == Idp || typ == Rhop
        for i = 1:length(j)
            j[i] = val_vec[i][sortperm(idp_vec[i])][start:move[k[i]]]
        end
    else
        for i = 1:length(j)
            id = sortperm(idp_vec[i])
            j[i] = val_vec[i][:, id][:, start:move[k[i]]]
        end
    end

    return j
end

end #PostSPH


text:
using ArgParse

function parse_cmd()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "-s"
            help = "Structured mesh"
            action = :store_true
        "--mpi"
            help = "Enable MPI"
            action = :store_true
        "N"
            help = "Number of elements in each direction"
            arg_type = Int
            default = 45
        "num_iters"
            help = "Number of iterations to run"
            arg_type = Int
            default = 10
        "--enzyme"
            help = "Run Enzyme"
            action = :store_true
    end

    return parse_args(s)
end

function getCacheCoherencePad(::T) where T
    return div(128, sizeof(T))
end

# PAD_DIV(nbytes, align) = nbytes + div(align - 1, align)
# PAD(nbytes, align) = PAD_DIV(nbytes,align) * align

@inline function calcElemVolume(
    x0, x1,x2, x3, x4, x5, x6, x7, y0, y1, y2, y3, y4, y5, y6, y7,
    z0, z1, z2, z3, z4, z5, z6, z7
)
    @inline function triple_product(x1, y1, z1, x2, y2, z2, x3, y3, z3)
        return (((x1)*((y2)*(z3) - (z2)*(y3)) + (x2)*((z1)*(y3)
                  - (y1)*(z3)) + (x3)*((y1)*(z2) - (z1)*(y2))))
    end
    #((x1)*((y2)*(z3) - (z2)*(y3)) + (x2)*((z1)*(y3) - (y1)*(z3)) + (x3)*((y1)*(z2) - (z1)*(y2)))

    twelveth = 1.0/12.0

    dx61 = x6 - x1
    dy61 = y6 - y1
    dz61 = z6 - z1

    dx70 = x7 - x0
    dy70 = y7 - y0
    dz70 = z7 - z0

    dx63 = x6 - x3
    dy63 = y6 - y3
    dz63 = z6 - z3

    dx20 = x2 - x0
    dy20 = y2 - y0
    dz20 = z2 - z0

    dx50 = x5 - x0
    dy50 = y5 - y0
    dz50 = z5 - z0

    dx64 = x6 - x4
    dy64 = y6 - y4
    dz64 = z6 - z4

    dx31 = x3 - x1
    dy31 = y3 - y1
    dz31 = z3 - z1

    dx72 = x7 - x2
    dy72 = y7 - y2
    dz72 = z7 - z2

    dx43 = x4 - x3
    dy43 = y4 - y3
    dz43 = z4 - z3

    dx57 = x5 - x7
    dy57 = y5 - y7
    dz57 = z5 - z7

    dx14 = x1 - x4
    dy14 = y1 - y4
    dz14 = z1 - z4

    dx25 = x2 - x5
    dy25 = y2 - y5
    dz25 = z2 - z5

    # 11 + 3*14
    volume = (
        triple_product(
            dx31 + dx72, dx63, dx20,
            dy31 + dy72, dy63, dy20,
            dz31 + dz72, dz63, dz20
        ) +
        triple_product(
            dx43 + dx57, dx64, dx70,
            dy43 + dy57, dy64, dy70,
            dz43 + dz57, dz64, dz70
        ) +
        triple_product(
            dx14 + dx25, dx61, dx50,
            dy14 + dy25, dy61, dy50,
            dz14 + dz25, dz61, dz50
        )
    )
    volume *= twelveth
    return volume
end

@inline calcElemVolume(x, y, z) = calcElemVolume(
    x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8],
    y[1], y[2], y[3], y[4], y[5], y[6], y[7], y[8],
    z[1], z[2], z[3], z[4], z[5], z[6], z[7], z[8]
)

function verifyAndWriteFinalOutput(elapsed_time, domain, nx, numRanks)
    # GrindTime1 only takes a single domain into account, and is thus a good way to measure
    # processor speed indepdendent of MPI parallelism.
    # GrindTime2 takes into account speedups from MPI parallelism.
    grindTime1 = ((elapsed_time*1e6)/domain.cycle)/(nx*nx*nx)
    grindTime2 = ((elapsed_time*1e6)/domain.cycle)/(nx*nx*nx*numRanks)

    @printf "Run completed:\n"
    @printf "   Problem size        =  %d\n" nx
    @printf "   MPI tasks           =  %d\n" numRanks
    @printf "   Iteration count     =  %d\n" domain.cycle
    @printf "   Final Origin Energy =  %.6e\n" domain.e[1]

    MaxAbsDiff = 0.0
    MaxRelDiff = 0.0
    TotalAbsDiff = 0.0
    for j in 1:nx
        for k in j+1:nx
            AbsDiff = abs(domain.e[(j-1)*nx+k] - domain.e[(k-1)*nx+j])
            TotalAbsDiff += AbsDiff
            if MaxAbsDiff < AbsDiff
                MaxAbsDiff = AbsDiff
            end
            RelDiff = AbsDiff / domain.e[(k-1)*nx+j]
            if MaxRelDiff < RelDiff
                MaxRelDiff = RelDiff
            end
        end
    end

    # Quick symmetry check
    @printf "   Testing Plane 0 of Energy Array on rank 0:\n"
    @printf "        MaxAbsDiff   = %.6e\n" MaxAbsDiff
    @printf "        TotalAbsDiff = %.6e\n" TotalAbsDiff
    @printf "        MaxRelDiff   = %.6e\n" MaxRelDiff

    # Timing information
    @printf "\nElapsed time         = %.2f (s)\n" elapsed_time
    @printf "Grind time (us/z/c)  =  %.8f  (per dom)  (%.8f overall)\n" grindTime1 elapsed_time
    @printf "FOM                  = %.8f (z/s)\n\n" 1000/grindTime2
end


text:
# ---
# title: 980. Unique Paths III
# id: problem980
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Backtracking, Depth-first Search
# link: <https://leetcode.com/problems/unique-paths-iii/description/>
# hidden: true
# ---
#
# On a 2-dimensional `grid`, there are 4 types of squares:
#
#   * `1` represents the starting square.  There is exactly one starting square.
#   * `2` represents the ending square.  There is exactly one ending square.
#   * `0` represents empty squares we can walk over.
#   * `-1` represents obstacles that we cannot walk over.
#
# Return the number of 4-directional walks from the starting square to the
# ending square, that **walk over every non-obstacle square  exactly once**.
#
#
#
# **Example 1:**
#
#
#
#     Input: [[1,0,0,0],[0,0,0,0],[0,0,2,-1]]
#     Output: 2
#     Explanation: We have the following two paths:
#     1. (0,0),(0,1),(0,2),(0,3),(1,3),(1,2),(1,1),(1,0),(2,0),(2,1),(2,2)
#     2. (0,0),(1,0),(2,0),(2,1),(1,1),(0,1),(0,2),(0,3),(1,3),(1,2),(2,2)
#
# **Example 2:**
#
#
#
#     Input: [[1,0,0,0],[0,0,0,0],[0,0,0,2]]
#     Output: 4
#     Explanation: We have the following four paths:
#     1. (0,0),(0,1),(0,2),(0,3),(1,3),(1,2),(1,1),(1,0),(2,0),(2,1),(2,2),(2,3)
#     2. (0,0),(0,1),(1,1),(1,0),(2,0),(2,1),(2,2),(1,2),(0,2),(0,3),(1,3),(2,3)
#     3. (0,0),(1,0),(2,0),(2,1),(2,2),(1,2),(1,1),(0,1),(0,2),(0,3),(1,3),(2,3)
#     4. (0,0),(1,0),(2,0),(2,1),(1,1),(0,1),(0,2),(0,3),(1,3),(1,2),(2,2),(2,3)
#
# **Example 3:**
#
#
#
#     Input: [[0,1],[2,0]]
#     Output: 0
#     Explanation:
#     There is no path that walks over every empty square exactly once.
#     Note that the starting and ending square can be anywhere in the grid.
#
#
#
#
# **Note:**
#
#   1. `1 <= grid.length * grid[0].length <= 20`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using FastTransforms, Compat
using Compat.Test, Compat.LinearAlgebra, Compat.Random
if VERSION < v"0.7"
    srand(0)
else
    using Compat.Statistics
    Random.seed!(0)
end

include("basictests.jl")

include("chebyshevtests.jl")

include("nuffttests.jl")

include("chebyshevjacobitests.jl")

include("chebyshevlegendretests.jl")

include("gaunttests.jl")

include("fftBigFloattests.jl")

include("paduatests.jl")

include("butterflytests.jl")

include("sphericalharmonics/sphericalharmonictests.jl")

include("toeplitztests.jl")


text:
module PoisonTests

using STLCutters
using Gridap
import Gridap: ∇
using GridapEmbedded
using Test

using STLCutters: compute_stl_model
using STLCutters: read_stl, merge_nodes, get_bounding_box

# Manufactured solution
u(x) = x[1] + x[2] - x[3]
∇u(x) = VectorValue( 1, 1, -1)
Δu(x) = 0
f(x) = - Δu(x)
ud(x) = u(x)
∇(::typeof(u)) = ∇u

@time X,T,N = read_stl(joinpath(@__DIR__,"data/Bunny-LowPoly.stl"))
#X,T,N = read_stl(joinpath(@__DIR__,"data/cube.stl"))
#@time X,T,N = read_stl(joinpath(@__DIR__,"data/441708_sf.obj"))
@time stl = compute_stl_model(T,X)
@time stl = merge_nodes(stl)
# writevtk(stl,"geo")
n = 20
δ = 0.2
pmin,pmax = get_bounding_box(stl)
diagonal = pmax-pmin
pmin = pmin - diagonal*δ
pmax = pmax + diagonal*δ
partition = (n,n,n)

geo = STLGeometry(stl)
bgmodel = CartesianDiscreteModel(pmin,pmax,partition)

# Cut the background model

cutter = STLCutter(threading=:threads,tolfactor=10^4)

@time cutgeo, = cut(cutter,bgmodel,geo)

# Setup integration meshes
Ω = Triangulation(cutgeo,PHYSICAL)
Γd = EmbeddedBoundary(cutgeo)

# Setup normal vectors
n_Γd = get_normal_vector(Γd)

#writevtk(Ω,"trian_O")
#writevtk(Γd,"trian_Gd",cellfields=["normal"=>n_Γd])
#writevtk(Triangulation(bgmodel),"bgtrian")

# Setup Lebesgue measures
order = 1
degree = 2*order
dΩ = Measure(Ω,degree)
dΓd = Measure(Γd,degree)

vol = sum( ∫(1)*dΩ  )
surf = sum( ∫(1)*dΓd )
@show vol,surf

# Setup FESpace
Ω_act = Triangulation(cutgeo,ACTIVE)
V = TestFESpace(Ω_act,ReferenceFE(lagrangian,Float64,order),conformity=:H1)
U = TrialFESpace(V)

# Weak form
γd = 10.0
γg = 0.1
h = (pmax - pmin)[1] / partition[1]

a(u,v) =
  ∫( ∇(v)⋅∇(u) ) * dΩ +
  ∫( (γd/h)*v*u  - v*(n_Γd⋅∇(u)) - (n_Γd⋅∇(v))*u ) * dΓd

l(v) =
  ∫( v*f ) * dΩ +
  ∫( (γd/h)*v*ud - (n_Γd⋅∇(v))*ud ) * dΓd

# FE problem
@time op = AffineFEOperator(a,l,U,V)
@time uh = solve(op)

e = u - uh

# Postprocess
l2(u) = sqrt(sum( ∫( u*u )*dΩ ))
h1(u) = sqrt(sum( ∫( u*u + ∇(u)⋅∇(u) )*dΩ ))

el2 = l2(e)
eh1 = h1(e)
ul2 = l2(uh)
uh1 = h1(uh)

#writevtk(Ω,"results",cellfields=["uh"=>uh])

@test el2/ul2 < 1.e-8
@test eh1/uh1 < 1.e-7

end # module


text:
"""
`FlipY <: ImageTransformation`

Description
============

Reverses the pixel order along the y-axis. In other words it
mirrors the Image along the x-axis (i.e. vertically).

If created using the parameter `chance`, the operation will be
lifted into a `Either`. See the documentation of `Either` for more
information.

Usage
======

    FlipY()

    FlipY(chance)

Arguments
==========

- **`chance`** : The probability of the operation to be applied to
the image if the operation is used in `transform`. Must be within
the set [0, 1].

Methods
========

- **`transform`** : Applies the transformation to the given Image
or set of images

Author(s)
==========

- Christof Stocker (Github: https://github.com/Evizero)

Examples
========

    using Images
    using Augmentor

    # Create a test image with 4 pixel.
    img = grayim(UInt8[200 150; 50 1])

    # Apply the image operation to the image
    img_new = transform(FlipY(), img)

    # Check manually if the result is as expected
    @assert img_new == grayim(UInt8[150 200; 1 50])

see also
=========

`ImageTransformation`, `Either`, `transform`
"""
immutable FlipY <: ImageTransformation
end

FlipY(chance) = ProbableOperation(FlipY(); chance = chance)

Base.showcompact(io::IO, op::FlipY) = print(io, "Flip y-axis.")
multiplier(::FlipY) = 1

function transform{T<:AbstractImage}(op::FlipY, img::T)
    result = copyproperties(img, flipy(img))::T
    _log_operation!(result, op)
end


text:
import DWARF: LEB128, LightDIERef, ULEB128, SLEB128, attr_color, DW_AT, DW_FORM, AttributeSpecification
using StructIO
import StructIO: fix_endian
using ObjFileBase
using AbstractTrees
import ObjFileBase: strtab_lookup
import AbstractTrees: printnode
using DWARF
import DWARF: readorskip
export Attribute

import Base: isequal, read, bytestring, ==, print, show, convert

# An attribute, encapsulating the name, form and value
immutable Attribute
    spec::AttributeSpecification
    value
end
convert{T}(::Type{Nullable{T}}, at::Attribute) = Nullable{T}(convert(T,at))
convert{T}(::Type{Ref{T}}, at::Attribute) = RefValue{T}(convert(T,at))
convert(::Type{Nullable{Attribute}}, at::Attribute) = Nullable{T}(at)
convert{T<:Integer}(::Type{T},at::Attribute) = convert(T, at.value)
convert(::Type{Attribute},at::Attribute) = at

function Base.show(io::IO, at::Attribute)
    indent = isa(io, IOContext) ? get(io, :indent, 0) : 0
    print(io," "^indent)
    if !haskey(DW_AT, at.spec.name)
        printfield_with_color(:red, io, string("Unknown ($(at.spec.name))"),18; align=:left)
    else
        printfield_with_color(attr_color, io, DW_AT[at.spec.name],18; align=:left)
    end
    printfield(io,string(" [", DW_FORM[at.spec.form], "] "),25, align = :left)
    if at.spec.form == DWARF.DW_FORM_addr
        print(io,"0x",hex(at.value,2sizeof(at.value)))
    elseif at.spec.form == DWARF.DW_FORM_sec_offset
        target_name = ".debug_unkown"
        if at.spec.name == DWARF.DW_AT_stmt_list
            target_name = ".debug_line"
        elseif at.spec.name == DWARF.DW_AT_location
            target_name = ".debug_loc"
        elseif at.spec.name == DWARF.DW_AT_start_scope || at.spec.name == DWARF.DW_AT_ranges
            target_name = ".debug_ranges"
        elseif at.spec.name == DWARF.DW_AT_macro_info
            target_name = ".debug_macinfo"
        end
        print(io,target_name,"[0x",hex(at.value.offset,
            2*sizeof(at.value.offset > typemax(UInt32) ?
                at.value.offset : UInt32)),"]")
    elseif at.spec.form == DWARF.DW_FORM_string
        print(io, String(at.value))
    elseif at.spec.name == DWARF.DW_AT_language && haskey(DWARF.DW_LANG, at.value)
        print(io, at.value, " (", DWARF.DW_LANG[at.value], ")")
    else
        showcompact(io, at.value)
    end
end

function Base.bytestring(at::Attribute, strtab = nothing)
    if at.spec.form == DWARF.DW_FORM_string
        return bytestring(at.value)
    elseif at.spec.form == DWARF.DW_FORM_strp
        @assert strtab !== nothing
        bytestring(at.value, strtab)
    else
        error("Not a string attribute")
    end
end

# Special form types
immutable StrTableReference
    offset::UInt64
end
StrTableReference(::Void) = nothing
Base.bytestring(ref::StrTableReference, strtab) = strtab_lookup(strtab, ref.offset)

function show(io::IO, x::StrTableReference)
    strtab = isa(io, IOContext) ? get(io, :strtab, nothing) : nothing
    if strtab === nothing
        print(io, ".debug_str[0x",
            hex(x.offset,2*sizeof(x.offset > typemax(UInt32) ?
                x.offset : UInt32)),"]")
    else
        show(io, strtab_lookup(strtab, x.offset))
    end
end

immutable SectionReference
    offset::UInt64
end
convert(::Type{UInt}, x::SectionReference) = convert(UInt, x.offset)
SectionReference(::Void) = nothing

immutable DwarfExpr{T}
    expr::Vector{UInt8}
end
(::Type{DwarfExpr{T}}){T}(::Void) = nothing
function show{T}(io::IO, expr::DwarfExpr{T})
    DWARF.Expressions.print_expression(io,T,expr.expr,:NativeEndian)
end

# Utility function for form processing
readorskip(io::IO, T::Type, endianness, ::Val{:skip}) = (skip(io, sizeof(T)); nothing)
readorskip(io::IO, T::Type, endianness, ::Val{:read}) = unpack(io, T, endianness)
readorskip{T<:LEB128}(io::IO, ::Type{T}, endianness, ::Val{:skip}) =
    while (read(io,UInt8)&0x80)!=0;end
readorskip{T<:LEB128}(io::IO, ::Type{T}, endianness, ::Val{:read}) = read(io, T)

readbytesorskip(io::IO, nbytes, ::Val{:skip}) = (skip(io, nbytes); nothing)
readbytesorskip(io::IO, nbytes, ::Val{:read}) = read(io, nbytes)

readorskipuntil(io::IO, delim::UInt8, ros) =
    (ret = readuntil(io, delim); ros == Val{:skip}() ? nothing : ret)

function make_ref(cu::LightDIERef, cuoffs)
    adjusted_offset = cu.cu.offset + cuoffs
    LightDIERef(cu.io, cu.cu, adjusted_offset)
end
make_ref(cu::LightDIERef, cuoffs::Void) = nothing

# Main function encoding how to process a given form
function readorskip(cu::LightDIERef, form, endianness, ros)
    io = cu.io
    cuheader = cu.cu.header
    if form == DWARF.DW_FORM_addr
        ret = readorskip(io, DWARF.size_to_inttype(cuheader.address_size), endianness, ros)
    elseif form == 0x2
        error("Invalid form")
    elseif form == DWARF.DW_FORM_block2
        length = read(io, UInt16)
        ret = readbytesorskip(io, length, ros)
    elseif form == DWARF.DW_FORM_block4
        length = read(io, UInt32)
        ret = readbytesorskip(io, length, ros)
    elseif form == DWARF.DW_FORM_data2
        ret = readorskip(io, UInt16, endianness, ros)
    elseif form == DWARF.DW_FORM_data4
        ret = readorskip(io, UInt32, endianness, ros)
    elseif form == DWARF.DW_FORM_data8
        ret = readorskip(io, UInt64, endianness, ros)
    elseif form == DWARF.DW_FORM_string
        ret = readorskipuntil(io, UInt8(0), ros)
    elseif form == DWARF.DW_FORM_block
        length = read(io, ULEB128)
        ret = readbytesorskip(io, UInt(length), ros)
    elseif form == DWARF.DW_FORM_exprloc
        length = read(io, ULEB128)
        ret = DwarfExpr{DWARF.size_to_inttype(cuheader.address_size)}(readbytesorskip(io, UInt(length), ros))
    elseif form == DWARF.DW_FORM_block1
        length = read(io, UInt8)
        ret = readbytesorskip(io, length, ros)
    elseif form == DWARF.DW_FORM_data1 || form == DWARF.DW_FORM_flag
        ret = readorskip(io, UInt8, endianness, ros)
    elseif form == DWARF.DW_FORM_sdata
        ret = readorskip(io, SLEB128, endianness, ros)
    elseif form == DWARF.DW_FORM_strp
        offs = readorskip(io, typeof(cuheader.unit_length), endianness, ros)
        ret = StrTableReference(offs)
    elseif form == DWARF.DW_FORM_udata
        ret = readorskip(io, ULEB128, endianness, ros)
    elseif form == DWARF.DW_FORM_ref_addr
        ret = readorskip(io, cuheader.version == 2 ?
            DWARF.size_to_inttype(cuheader.address_size) :
            typeof(cuheader.unit_length), endianness, ros)
    elseif form == DWARF.DW_FORM_ref1
        ret = make_ref(cu, readorskip(io, UInt8, endianness, ros))
    elseif form == DWARF.DW_FORM_ref2
        ret = make_ref(cu, readorskip(io, UInt16, endianness, ros))
    elseif form == DWARF.DW_FORM_ref4
        ret = make_ref(cu, readorskip(io, UInt32, endianness, ros))
    elseif form == DWARF.DW_FORM_ref8
        ret = make_ref(cu, readorskip(io, UInt64, endianness, ros))
    elseif form == DWARF.DW_FORM_ref_udata
        ret = make_ref(cu, UInt(readorskip(io, ULEB128, endianness, ros)))
    elseif form == DWARF.DW_FORM_indirect
        form = read(io, ULEB128)
        return readorskip(io, form, ros, endianness)
    elseif form == DWARF.DW_FORM_sec_offset
        offs = readorskip(io, typeof(cuheader.unit_length), endianness, ros)
        ret = SectionReference(offs)
    elseif form == DWARF.DW_FORM_flag_present
        ret = nothing
    elseif form == DWARF.DW_FORM_ref_sig8
        ret = readorskip(io, UInt64, endianness, ros)
    end
    (form, ret)
end


text:
function ispalindrome(n)
  s = string(n)
  s == reverse(s)
end

function radd(n)
  n + parse(BigInt, reverse(string(n)))
end

function islychrel(n)
  for i = 1:50
    n = radd(n)
    if ispalindrome(n)
      return false
    end
  end
  true
end

c=0
for i=1:10000
  global c
  if islychrel(BigInt(i))
    c += 1
  end
end

println(c)

max_stars_repo_name: pdbartlett/misc-stuff

max_issues_repo_name: pdbartlett/misc-stuff

max_forks_repo_name: pdbartlett/misc-stuff


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
Random.seed!(233)

if Sys.islinux()
py"""
import tensorflow as tf
libDirichletBD = tf.load_op_library('build/libDirichletBD.so')
@tf.custom_gradient
def dirichlet_bd(ii,jj,dof,vv):
    uu = libDirichletBD.dirichlet_bd(ii,jj,dof,vv)
    def grad(dy):
        return libDirichletBD.dirichlet_bd_grad(dy, uu, ii,jj,dof,vv)
    return uu, grad
"""
elseif Sys.isapple()
py"""
import tensorflow as tf
libDirichletBD = tf.load_op_library('build/libDirichletBD.dylib')
@tf.custom_gradient
def dirichlet_bd(ii,jj,dof,vv):
    uu = libDirichletBD.dirichlet_bd(ii,jj,dof,vv)
    def grad(dy):
        return libDirichletBD.dirichlet_bd_grad(dy, uu, ii,jj,dof,vv)
    return uu, grad
"""
elseif Sys.iswindows()
py"""
import tensorflow as tf
libDirichletBD = tf.load_op_library('build/libDirichletBD.dll')
@tf.custom_gradient
def dirichlet_bd(ii,jj,dof,vv):
    uu = libDirichletBD.dirichlet_bd(ii,jj,dof,vv)
    def grad(dy):
        return libDirichletBD.dirichlet_bd_grad(dy, uu, ii,jj,dof,vv)
    return uu, grad
"""
end

dirichlet_bd = py"dirichlet_bd"

# TODO:
ii = constant([1;2;3;2], dtype=Int32)
jj = constant([1;2;3;3], dtype=Int32)
vv = constant([2.0;3.0;4.0;4.0])
dof = constant([3], dtype=Int32)
u = dirichlet_bd(ii,jj,dof,vv)
sess = Session()
init(sess)
run(sess, u)
# error("")
# TODO:

# gradient check -- v
function scalar_function(m)
    return sum(tanh(dirichlet_bd(ii,jj,dof,m)))
end

m_ = vv
v_ = rand(4)
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session()
init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")


text:
include("compat.jl")
using HTTP

keep_open_stream = nothing

@testset "HTTP.Client" begin

using JSON

status(r) = r.status

for sch in ("http", "https")
    println("running $sch client tests...")

    println("simple GET, HEAD, POST, DELETE, etc.")
    @test status(HTTP.get("$sch://httpbin.org/ip")) == 200
    @test status(HTTP.head("$sch://httpbin.org/ip")) == 200
    @test status(HTTP.options("$sch://httpbin.org/ip")) == 200
    @test status(HTTP.post("$sch://httpbin.org/ip"; status_exception=false)) == 405
    @test status(HTTP.post("$sch://httpbin.org/post")) == 200
    @test status(HTTP.put("$sch://httpbin.org/put")) == 200
    @test status(HTTP.delete("$sch://httpbin.org/delete")) == 200
    @test status(HTTP.patch("$sch://httpbin.org/patch")) == 200

    # Testing within tasks, see https://github.com/JuliaWeb/HTTP.jl/issues/18
    println("async client request")
    @test status(fetch(@async HTTP.get("$sch://httpbin.org/ip"))) == 200

    @test status(HTTP.get("$sch://httpbin.org/encoding/utf8")) == 200

    println("pass query to uri")
    r = HTTP.get(merge(HTTP.URI("$sch://httpbin.org/response-headers"); query=Dict("hey"=>"dude")))
    h = Dict(r.headers)
    @test (haskey(h, "Hey") ? h["Hey"] == "dude" : h["hey"] == "dude")

    println("cookie requests")
    empty!(HTTP.CookieRequest.default_cookiejar)
    empty!(HTTP.DEFAULT_CLIENT.cookies)
    r = HTTP.get("$sch://httpbin.org/cookies", cookies=true)
    body = String(r.body)
    @test replace(replace(body, " "=>""), "\n"=>"")  == "{\"cookies\":{}}"
    r = HTTP.get("$sch://httpbin.org/cookies/set?hey=sailor&foo=bar", cookies=true)
    @test status(r) == 200
    body = String(r.body)
    @test replace(replace(body, " "=>""), "\n"=>"")  == "{\"cookies\":{\"foo\":\"bar\",\"hey\":\"sailor\"}}"

    # r = HTTP.get("$sch://httpbin.org/cookies/delete?hey")
    # @test String(take!(r)) == "{\n  \"cookies\": {\n    \"hey\": \"\"\n  }\n}\n"

    # stream
    println("client streaming tests")
    r = HTTP.post("$sch://httpbin.org/post"; body="hey")
    @test status(r) == 200
    # stream, but body is too small to actually stream
    r = HTTP.post("$sch://httpbin.org/post"; body="hey", stream=true)
    @test status(r) == 200
    r = HTTP.get("$sch://httpbin.org/stream/100")
    @test status(r) == 200
    bytes = r.body
    a = [JSON.parse(l) for l in split(chomp(String(bytes)), "\n")]
    totallen = length(bytes) # number of bytes to expect
    begin
        io = Base.BufferStream()
        r = HTTP.get("$sch://httpbin.org/stream/100"; response_stream=io)
        @test status(r) == 200

        b = [JSON.parse(l) for l in eachline(io)]
        @test a == b
    end

    # body posting: Vector{UInt8}, String, IOStream, IOBuffer, Base.BufferStream
    println("client body posting of various types")
    @test status(HTTP.post("$sch://httpbin.org/post"; body="hey")) == 200
    @test status(HTTP.post("$sch://httpbin.org/post"; body=UInt8['h','e','y'])) == 200
    io = IOBuffer("hey"); seekstart(io)
    @test status(HTTP.post("$sch://httpbin.org/post"; body=io)) == 200
    tmp = tempname()
    open(f->write(f, "hey"), tmp, "w")
    io = open(tmp)
    @test status(HTTP.post("$sch://httpbin.org/post"; body=io, enablechunked=false)) == 200
    close(io); rm(tmp)
    f = Base.BufferStream()
    write(f, "hey")
    close(f)
    @test status(HTTP.post("$sch://httpbin.org/post"; body=f, enablechunked=false)) == 200

    # chunksize
    #
    #     https://github.com/JuliaWeb/HTTP.jl/issues/60
    #     Currently httpbin.org responds with 411 status and “Length Required”
    #     message to any POST/PUT requests that are sent using chunked encoding
    #     See https://github.com/kennethreitz/httpbin/issues/340#issuecomment-330176449
    println("client transfer-encoding chunked")
    @test status(HTTP.post("$sch://httpbin.org/post"; body="hey", #=chunksize=2=#)) == 200
    @test status(HTTP.post("$sch://httpbin.org/post"; body=UInt8['h','e','y'], #=chunksize=2=#)) == 200
    io = IOBuffer("hey"); seekstart(io)
    @test status(HTTP.post("$sch://httpbin.org/post"; body=io, #=chunksize=2=#)) == 200
    tmp = tempname()
    open(f->write(f, "hey"), tmp, "w")
    io = open(tmp)
#    @test_broken status(HTTP.post("$sch://httpbin.org/post"; body=io, #=chunksize=2=#)) == 200
    close(io); rm(tmp)
    f = Base.BufferStream()
    write(f, "hey")
    close(f)
#    @test_broken status(HTTP.post("$sch://httpbin.org/post"; body=f, #=chunksize=2=#)) == 200

    # multipart
    println("client multipart body")
    r = HTTP.post("$sch://httpbin.org/post"; body=Dict("hey"=>"there"))
    @test status(r) == 200
    @test startswith(replace(replace(String(r.body), " "=>""), "\n"=>""), "{\"args\":{},\"data\":\"\",\"files\":{},\"form\":{\"hey\":\"there\"}")

    r = HTTP.post("$sch://httpbin.org/post"; body=Dict("hey"=>"there"))
    @test status(r) == 200
    @test startswith(replace(replace(String(r.body), " "=>""), "\n"=>""), "{\"args\":{},\"data\":\"\",\"files\":{},\"form\":{\"hey\":\"there\"}")

    tmp = tempname()
    open(f->write(f, "hey"), tmp, "w")
    io = open(tmp)
    r = HTTP.post("$sch://httpbin.org/post"; body=Dict("hey"=>"there", "iostream"=>io))
    close(io); rm(tmp)
    @test status(r) == 200
    str = replace(replace(String(r.body), " "=>""), "\n"=>"")
    @test startswith(str, "{\"args\":{},\"data\":\"\",\"files\":{\"iostream\":\"hey\"},\"form\":{\"hey\":\"there\"}")

    tmp = tempname()
    open(f->write(f, "hey"), tmp, "w")
    io = open(tmp)
    r = HTTP.post("$sch://httpbin.org/post"; body=Dict("hey"=>"there", "iostream"=>io))
    close(io); rm(tmp)
    @test status(r) == 200
    @test startswith(replace(replace(String(r.body), " "=>""), "\n"=>""), "{\"args\":{},\"data\":\"\",\"files\":{\"iostream\":\"hey\"},\"form\":{\"hey\":\"there\"}")

    tmp = tempname()
    open(f->write(f, "hey"), tmp, "w")
    io = open(tmp)
    m = HTTP.Multipart("mycoolfile.txt", io)
    r = HTTP.post("$sch://httpbin.org/post"; body=Dict("hey"=>"there", "multi"=>m))
    close(io); rm(tmp)
    @test status(r) == 200
    @test startswith(replace(replace(String(r.body), " "=>""), "\n"=>""), "{\"args\":{},\"data\":\"\",\"files\":{\"multi\":\"hey\"},\"form\":{\"hey\":\"there\"}")

    tmp = tempname()
    open(f->write(f, "hey"), tmp, "w")
    io = open(tmp)
    m = HTTP.Multipart("mycoolfile", io, "application/octet-stream")
    r = HTTP.post("$sch://httpbin.org/post"; body=Dict("hey"=>"there", "multi"=>m), #=chunksize=1000=#)
    close(io); rm(tmp)
    @test status(r) == 200
    @test startswith(replace(replace(String(r.body), " "=>""), "\n"=>""), "{\"args\":{},\"data\":\"\",\"files\":{\"multi\":\"hey\"},\"form\":{\"hey\":\"there\"}")

    # asynchronous
    println("asynchronous client request body")
    begin
        f = Base.BufferStream()
        write(f, "hey")
        t = @async HTTP.post("$sch://httpbin.org/post"; body=f, enablechunked=false)
        #fetch(f) # fetch for the async call to write it's first data
        write(f, " there ") # as we write to f, it triggers another chunk to be sent in our async request
        write(f, "sailor")
        close(f) # setting eof on f causes the async request to send a final chunk and return the response
        @test status(fetch(t)) == 200
    end

    # redirects
    println("client redirect following")
    for meth in ("GET","HEAD")
        @test status(HTTP.request(meth, "$sch://httpbin.org/redirect/1")) ==200
        @test status(HTTP.request(meth, "$sch://httpbin.org/redirect/1", redirect=false)) == 302
        @test status(HTTP.request(meth, "$sch://httpbin.org/redirect/6")) == 302 #over max number of redirects
        @test status(HTTP.request(meth, "$sch://httpbin.org/relative-redirect/1")) == 200
        @test status(HTTP.request(meth, "$sch://httpbin.org/absolute-redirect/1")) == 200
        @test status(HTTP.request(meth, "$sch://httpbin.org/redirect-to?url=http%3A%2F%2Fgoogle.com")) == 200
    end

    @test status(HTTP.post("$sch://httpbin.org/post"; body="√")) == 200
    println("client basic auth")
    @test status(HTTP.get("$sch://user:pwd@httpbin.org/basic-auth/user/pwd"; basic_authorization=true)) == 200
    @test status(HTTP.get("$sch://user:pwd@httpbin.org/hidden-basic-auth/user/pwd"; basic_authorization=true)) == 200
    @test status(HTTP.get("$sch://test:%40test@httpbin.org/basic-auth/test/%40test"; basic_authorization=true)) == 200

    # custom client & other high-level entries
    println("high-level client request methods")
    cli = HTTP.Client()
#=
    buf = IOBuffer()
    cli = HTTP.Client(buf)
    HTTP.get(cli, "$sch://httpbin.org/ip")
    seekstart(buf)
    @test length(String(take!(buf))) > 0
=#

    r = HTTP.request("GET", "$sch://httpbin.org/ip")
    @test status(r) == 200

    uri = HTTP.URI("$sch://httpbin.org/ip")
    r = HTTP.request("GET", uri)
    @test status(r) == 200
    r = HTTP.get(uri)
    @test status(r) == 200
    r = HTTP.get(cli, uri)
    @test status(r) == 200

    r = HTTP.request("GET", "$sch://httpbin.org/ip")
    @test status(r) == 200

    uri = HTTP.URI("$sch://httpbin.org/ip")
    r = HTTP.request("GET", uri)
    @test status(r) == 200

#= FIXME
    req = HTTP.Request(HTTP.GET, uri, HTTP.Headers(), HTTP.FIFOBuffer())
    r = HTTP.request(req)
    @test status(r) == 200
    @test HTTP.request(r) !== nothing
    @test length(take!(r)) > 0
=#

#=
    for c in HTTP.DEFAULT_CLIENT.httppool["httpbin.org"]
        HTTP.dead!(c)
    end
=#

    r = HTTP.get(cli, "$sch://httpbin.org/ip")
#    @test isempty(HTTP.cookies(r))
#    @test isempty(HTTP.history(r))

    r = HTTP.get("$sch://httpbin.org/image/png")
    @test status(r) == 200

    # ensure we can use AbstractString for requests
    r = HTTP.get(SubString("http://httpbin.org/ip",1))

    # canonicalizeheaders
    @test status(HTTP.get("$sch://httpbin.org/ip"; canonicalizeheaders=false)) == 200

    # r = HTTP.connect("http://47.89.41.164:80")
    # gzip body = "hey"
    # body = UInt8[0x1f,0x8b,0x08,0x00,0x00,0x00,0x00,0x00,0x00,0x03,0xcb,0x48,0xad,0x04,0x00,0xf0,0x15,0xd6,0x88,0x03,0x00,0x00,0x00]
    # r = HTTP.post("$sch://httpbin.org/post"; body=body, chunksize=1)

    @testset "Keep stream open" begin
        println("keep stream open")
        # Make a WebSocket upgrade which won't immediately close the socket.
        headers = Dict(
            "Upgrade" => "websocket",
            "Connection" => "Upgrade",
            "Sec-WebSocket-Key" => "dGhlIHNhbXBsZSBub25jZQ==",
            "Sec-WebSocket-Version" => "13")

        r = HTTP.open("GET", "ws://echo.websocket.org", headers; reuse_limit=0, keep_open=true) do http
            HTTP.startread(http)
            global keep_open_stream = HTTP.ConnectionPool.getrawstream(http)
        end
        @test isopen(keep_open_stream)
    end
end

end # @testset "HTTP.Client"


text:
# actually define the distribution.

struct ZeroTruncatedPoisson{T<:Real} <: DiscreteUnivariateDistribution
    λ::T
    ZeroTruncatedPoisson{T}(λ::Real) where {T <: Real} = new{T}(λ)
end

function ZeroTruncatedPoisson(λ::T; check_args=true) where {T <: Real}
    check_args && Distributions.@check_args(Poisson, λ > zero(λ))
    return ZeroTruncatedPoisson{T}(λ)
end

function Distributions.pdf(d::ZeroTruncatedPoisson, x::Real)
    pois = Poisson(d.λ)
    pois_pdf = Distributions.pdf(pois, x)
    pois_zero = 1-Distributions.pdf(pois, zero(x))
    return x > 0 ? pois_pdf/pois_zero : zero(pois_pdf)
end

function _logpdf(d::ZeroTruncatedPoisson, x::T) where {T<:Real}
    pois = Poisson(d.λ)
    pois_zero = logccdf(pois, 0)
    if x >= 1
        logpdf(pois, x) - pois_zero
    else
        TF = float(T)
        -TF(Inf)
    end
end

function Distributions.logpdf(d::ZeroTruncatedPoisson, x::T) where {T<:Real}
    isinteger(x) || return -TF(Inf)
    return _logpdf(d, x)
end

function _cdf(d::ZeroTruncatedPoisson, x::Real)
    pois = Poisson(d.λ)
    pois_pdf = Distributions.pdf(pois, zero(x))
    pois_cdf = Distributions.cdf(pois, x) - pois_pdf
    pois_zero = 1-pois_pdf
    return x > 0 ? pois_cdf/pois_zero : zero(pois_cdf)
end

Distributions.cdf(d::ZeroTruncatedPoisson, x::Real) = _cdf(d,x)
Distributions.cdf(d::ZeroTruncatedPoisson, x::Integer) = _cdf(d,x)

"""
    TruncatedPoissonSample(Z, E)

An observed sample ``Z`` drawn from a truncated Poisson distribution,

```math
Z \\sim \\text{Poisson}(\\mu \\cdot E)  \\mid Z \\geq 1.
```

The multiplying intensity ``E`` is assumed to be known (and equal to `1.0` by default), while
``\\mu`` is assumed unknown. The type above is used when the sample ``Z`` is to be used for estimation or inference of ``\\mu``.

```julia
TruncatedPoissonSample(3)
TruncatedPoissonSample(3, 1.5)
```
"""
struct TruncatedPoissonSample{T,S} <: DiscreteEBayesSample{T}
    Z::T
    E::S
end

TruncatedPoissonSample(Z) = TruncatedPoissonSample(Z, 1.0)
TruncatedPoissonSample() = TruncatedPoissonSample(missing)

response(Z::TruncatedPoissonSample) = Z.Z
nuisance_parameter(Z::TruncatedPoissonSample) = Z.E

function likelihood_distribution(Z::TruncatedPoissonSample, λ)
    ZeroTruncatedPoisson(λ * nuisance_parameter(Z))
end

summarize_by_default(Zs::Vector{<:TruncatedPoissonSample}) = skedasticity(Zs) == Homoskedastic()

function Base.show(io::IO, Z::TruncatedPoissonSample)
    resp_Z = response(Z)
    if ismissing(resp_Z) || isa(resp_Z, Interval)
        spaces_to_keep = 1
    else
        spaces_to_keep = max(3 - ndigits(response(Z)), 1)
    end
    spaces = repeat(" ", spaces_to_keep)
    print(io, "Z=", resp_Z, spaces, "| ", "E=", Z.E)
end


text:
using Church
using Distributions

mem(f::Function) = begin
    d = Dict()
    nf(args...) = begin
        if !haskey(d, args)
            d[args] = f(args...)
        end
        d[args]
    end
    nf
end

dp(concentration::Real, base_measure::Function) = begin
    sticks = mem(i::Int -> sample(Beta(1., concentration)))  #Dict{Int, Sample}
    atoms  = mem(i::Int -> base_measure())                   #Dict{Int, Sample}
    loop(i::Int) =
        @branch(sample(Bernoulli[sticks(i)]), atoms(i), loop(i+1))
    () -> loop(0)
end

dp_mixture(concentration::Real, base_measure::Function, parameter::Function) = begin
    dp_ = dp(concentration, base_measure)
    () -> parameter(dp_())
end

d = dp(1., () -> sample(Normal()))
ds = [d() for i = 1:20]

f = dp_mixture(1., () -> sample(Normal()), x -> sample(Normal[x, 0.1]))
fs = [f() for i = 1:20]


text:
using TensorNetworkAD
using Test
using Random

@testset "TensorNetworkAD.jl" begin
    @testset "autodiff" begin
        println("autodiff tests running...")
        include("autodiff.jl")
    end

    @testset "example tensors" begin
        println("exampletensors tests running...")
        include("exampletensors.jl")
    end

    @testset "fixedpoint" begin
        println("fixedpoint tests running...")
        include("fixedpoint.jl")
    end

    @testset "trg" begin
        println("trg tests running...")
        include("trg.jl")
    end

    @testset "ctmrg" begin
        println("ctmrg tests running...")
        include("ctmrg.jl")
    end

    @testset "variationalipeps" begin
        println("variationalipeps tests running...")
        include("variationalipeps.jl")
    end
end


text:
using MarketData

facts("constructors") do

  context("Ticker") do
      @fact 1 => roughly(1.0000000001)
  end

  context("CUSIP") do
      @fact 1 => roughly(1.0000000001)
  end

  context("BloombergID") do
      @fact 1 => roughly(1.0000000001)
  end

  context("ReutersID") do
      @fact 1 => roughly(1.0000000001)
  end

  context("Currency") do
      @fact 1 => roughly(1.0000000001)
  end

  context("CurrencyPair") do
      @fact 1 => roughly(1.0000000001)
  end

  context("Stock default constructor") do
      @fact 1 => roughly(1.0000000001)
  end

  context("Stock from Ticker") do
      @fact 1 => roughly(1.0000000001)
  end
end


text:
using IniFile
using Test

ini = Inifile()
read(ini, joinpath(dirname(@__FILE__),"test.ini"))

# section exists
@test has_section(ini, "not present") == false
@test has_section(ini, "section 1") == true

# get value or defaults
@test get(ini, "section 1", "strname1") == "section1name1"
@test get(ini, "section 1", "strname2") == "section 1 name 2"
@test get(ini, "section 1", "intname1") == "1"
@test get(ini, "section 1", "floatname1") == "1.1"
@test get(ini, "section 2", "floatname1") == "1.1"
@test get(ini, "section 2", "intname1") == "10"

# int, float and bool types
@test get_int(ini, "section 2", "intname1") == 10
@test get_float(ini, "section 2", "floatname1") == 1.1
@test get_bool(ini, "section 1", "boolname1") == true
@test get_bool(ini, "section 1", "boolname2") == true
@test get_bool(ini, "section 1", "boolname3") == false
@test get_bool(ini, "section 1", "boolname4") == false
@test set(ini, "section 2", "floatname1", 9.1) == 9.1
@test get(ini, "section 2", "floatname1") == 9.1
@test set(ini, "section 2", "intname1", 9) == 9
@test get(ini, "section 2", "intname1") == 9
@test set(ini, "section 2", "boolname1", true) == true
@test get(ini, "section 2", "boolname1") == true
@test set(ini, "section 2", "boolname1", false) == false
@test get(ini, "section 2", "boolname1") == false

# test delimiters
@test get(ini, "section 2", "delimtest1") == "value with :"
@test get(ini, "section 2", "delimtest2") == "name = value"

# key/section not found
@test get(ini, "section 2", "noname") == :notfound
@test get(ini, "section 3", "noname") == :notfound
@test get(ini, "section 3", "noname", nothing) == nothing
@test get(ini, "section 3", "strname1") == "val1"

# create a new section
@test set(ini, "section 3", "strname1", "new string") == "new string"
@test get(ini, "section 3", "strname1") == "new string"
@test get(ini, "section 3", "intname1") == "1"

iob = IOBuffer()
show(iob, ini)
buff = String(take!(iob))

iob = IOBuffer(buff)
ini = Inifile()
read(ini, iob)
@test get(ini, "section 3", "strname1") == "new string"
@test get(ini, "section 3", "intname1") == "1"


text:
"""
    distance_profile(dist, Q, T; kwargs...)

Calculate the distance profile corresponding to sliding `Q` over `T` and measuring distances between time windows using `dist`

The output is a vector of length `length(T) - length(Q) + 1`.
If inputs are matrices or higher-dimensional arrays, time is considered to be the last axis, and windows on the form `T[:, t:t+lastlength(Q)-1]`, where `lastlength` measures the length along the last axis, are compared.

Keyword arguments are sent to `dist`, which is called like this `dist(Q,T; kwargs...)`. `dist` can be a function or a distance from Distances.jl.

`distance_profile!(D, dist, Q, T; kwargs...)` is an in-place version.
"""
function distance_profile(dist, Q::AbstractArray{S}, T::AbstractArray{S}; kwargs...) where S
    m = lastlength(Q)
    n = lastlength(T)
    n >= m || throw(ArgumentError("Q cannot be longer than T"))
    l = n-m+1
    D = map(1:l) do i
        dist(Q, getwindow(T, m, i); kwargs...)
    end
end

"$(SIGNATURES)"
function distance_profile!(D::AbstractArray, dist, Q::AbstractArray{S}, T::AbstractArray{S}; kwargs...) where S
    m = lastlength(Q)
    n = lastlength(T)
    n >= m || throw(ArgumentError("Q cannot be longer than T"))
    l = n-m+1
    length(D) >= l || throw(ArgumentError("D is too short"))
    map!(D, 1:l) do i
        dist(Q, getwindow(T, m, i); kwargs...)
    end
end

"$(SIGNATURES)"
function distance_profile!(D::AbstractVector{S},::ZEuclidean, QT::AbstractVector{S}, μ, σ, m::Int, i::Int) where S <: Number
    @assert i <= length(D)
    mμ = m*μ[i]
    mσ = m*σ[i]
    O = zero(S)
    # @inbounds @fastmath @simd for j = eachindex(D)
    @avx unroll=4 for j = eachindex(D)
        # frac = (QT[j] - mμ*μ[j]) / (mσ*σ[j])
        # frac = (QT[j] - mμ*μ[j]) * Base.FastMath.inv_fast(mσ*σ[j])
        frac = (QT[j] - mμ*μ[j]) * LoopVectorization.VectorizationBase.vinv_fast(mσ*σ[j])
        D[j] = sqrt(max(2m*(1-frac), O))
    end
    D[i] = typemax(eltype(D))
    D
end

"$(SIGNATURES)"
distance_profile(
    ::ZEuclidean,
    QT::AbstractVector{S},
    μ::AbstractVector{S},
    σ::AbstractVector{S},
    m::Int,
) where {S<:Number} = distance_profile!(similar(μ), ZEuclidean(), QT, μ, σ, m, 1)

"""
    $(SIGNATURES)

Accepts precomputed sliding mean and std of the input arrays. `QT` is the windowed dot product and `i` is the index into `T` (the longer time series).
"""
function distance_profile!(D::AbstractVector{S},::ZEuclidean, QT::AbstractVector{S}, μA, σA, μT, σT, m::Int, i::Int) where S <: Number
    @assert i <= length(μA)
    @avx for j = eachindex(D,QT,μT,σT)
        frac = (QT[j] - m*μA[i]*μT[j]) / (m*σA[i]*σT[j])
        D[j] = sqrt(max(2m*(1-frac), 0))
    end
    D
end

"$(SIGNATURES)"
distance_profile(::ZEuclidean, QT::AbstractVector{S}, μA, σA, μT, σT, m::Int) where {S<:Number} =
    distance_profile!(similar(μT), ZEuclidean(), QT, μA, σA, μT, σT, m, 1)

"""
    distance_profile(::ZEuclidean, Q, T)

Compute the z-normalized Euclidean distance profile corresponding to sliding `Q` over `T`
"""
function distance_profile!(
    D::AbstractVector{S},
    ::ZEuclidean,
    Q::AbstractVector{S},
    T::AbstractVector{S},
) where {S<:Number}
    m = length(Q)
    μ, σ = sliding_mean_std(T, m)
    QT = window_dot(znorm(Q), T) # TODO: allocates a new znorm(Q) each time
    @avx for j in eachindex(D,QT,σ)
        frac = QT[j] / (m * σ[j])
        D[j] = sqrt(max(2m * (1 - frac), 0))
    end
    D
end

"$(SIGNATURES)"
distance_profile(::ZEuclidean, Q::AbstractArray{S}, T::AbstractArray{S}) where {S} =
    distance_profile!(similar(T, length(T) - length(Q) + 1), ZEuclidean(), Q, T)

"""
The dot product between query Q and all subsequences of the same length as Q in time series T
"""
function window_dot(Q, T)
    n   = length(T)
    m   = length(Q)
    QT  = conv(reverse(Q), T)
    return QT[m:n]
end

"""
    $(SIGNATURES)

return vectors with mean and std of sliding windows of length `m`
"""
function sliding_mean_std(x::AbstractVector{T}, m) where T
    @assert length(x) >= m
    n = length(x)-m+1
    s = ss = zero(float(T))
    μ = Vector{float(T)}(undef, n)
    σ = Vector{float(T)}(undef, n)
    @avx for i = 1:m
        s  += x[i]
        ss += x[i]^2
    end
    μ[1] = s/m
    σ[1] = sqrt(max(ss/m - μ[1]^2, 0))
    @fastmath @inbounds for i = 1:n-1 # fastmath making it more accurate here as well, but not faster
        s -= x[i]
        ss -= x[i]^2
        s += x[i+m]
        ss += x[i+m]^2
        μ[i+1] = s/m
        σ[i+1] = sqrt(max(ss/m - μ[i+1]^2, 0))
    end
    μ,σ
end

@inbounds @fastmath function sliding_mean_std(x::AbstractMatrix{T}, m) where T
    @assert lastlength(x) >= m
    n = lastlength(x)-m+1
    s  = zeros(T, size(x,1))
    ss = zeros(T, size(x,1))
    μ = Matrix{float(T)}(undef, size(x,1), n)
    σ = Matrix{float(T)}(undef, size(x,1), n)
    for i = 1:m
        @avx s  .+= view(x, :, i)
        @avx ss .+= view(x, :, i).^2
    end
    @avx μ[:,1] .= s./m
    @avx σ[:,1] .= sqrt.(max.(ss./m .- view(μ, :, 1).^2, 0))
    for i = 1:n-1 # fastmath making it more accurate here as well, but not faster
        @avx s .-= view(x, :, i)
        @avx ss .-= view(x, :, i).^2
        @avx s .+= view(x, :, i+m)
        @avx ss .+= view(x, :, i+m).^2
        @avx μ[:,i+1] .= s./m
        @avx σ[:,i+1] .= sqrt.(max.(ss./m .- view(μ, :, i+1).^2, 0))
    end
    μ,σ
end

"""
    $(SIGNATURES)

return mean of sliding windows of length `m`. Operates in-place and stores result in first argument
"""
function sliding_mean!(μ,x::AbstractVector{T}, m) where T
    n = length(x)-m+1
    @assert length(x) >= m
    @assert length(μ) >= n

    s = zero(T)
    @avx for i = 1:m
        s += x[i]
    end
    μ[1] = s/m
    @fastmath @inbounds for i = 1:n-1 # fastmath making it more accurate here as well, but not faster
        s -= x[i]
        s += x[i+m]
        μ[i+1] = s/m
    end
    μ
end

"""
$(SIGNATURES)

Returns the sliding-window entropy
"""
function sliding_entropy(x::AbstractVecOrMat{T}, m=actuallastlength(x)) where T
    N = actuallastlength(x)
    n = N-m+1
    e = Vector{T}(undef, n)
    if x isa AbstractVector
        ent = vmap(x->x*log(x), x)
    else
        ent = zeros(T, N)
        @avx for i = 1:N
            for j = 1:size(x,1)
                ent[i] += x[j,i]*log(x[j,i])
            end
        end
    end
    e0 = zero(T)
    @avx for i = 1:m
        e0 += ent[i]
    end
    e[1] = e0
    @fastmath @inbounds for i = 1:n-1
        e0 -= ent[i]
        e0 += ent[i+m]
        e[i+1] = e0
    end
    e .= .-e
end

"""
$(SIGNATURES)

Returns the sliding-window entropy, normalized to have unit sum for each feature channel.
"""
function sliding_entropy_normalized(x::AbstractArray{T}, m=actuallastlength(x)) where T
    sums = sum(x, dims=ndims(x))
    sliding_entropy(x ./ sums, m)
end

"""
    $(SIGNATURES)

Sliding-window PCA. Returns `U1 ∈ 𝐑(N)` the main principal-component loading scaled with the singular value, and `V1 ∈ 𝐑(d × N)` the main principal direction.

# Arguments:
- `X`: Signal ∈ 𝐑(d × N)
- `fs`: Sampling frequency
- `w`: Window length in seconds
- `cutoff_freq`: Frequencies below this will be removed. Set to 0 to avoid filtering.
- `pad`: if true, the return value will be padded with zeros to account for the fact that the analysis uses sliding windows and produce a slightly shorter output. If `pad = true`, the output will have the same length as the input.
- `normalize_win`: remove the mean from each window before calculating svd.
- `center`: Index into `X` at the center of the window (default) or at a position in the window proportional to how far along in the signal the window is taken.
"""
function sliding_pca(X::AbstractMatrix{A}; fs, w = 0.3, cutoff_freq = 0, pad = true, normalize_win=false, center = true) where A
    w = round(Int, w*fs) # from seconds to samples
    d,N = size(X)
    mX = mean(X, dims=2)
    X = X .- mX
    if cutoff_freq > 0
        X = filtfilt(digitalfilter(Highpass(cutoff_freq; fs=fs), Butterworth(6)), X')'
    end

    win = 1:w
    dominant_mode = Vector{A}(undef, N-w+1)
    directions = Matrix{A}(undef, d, N-w+1)
    d1 = zeros(A, d); d1[1] = 1
    wininds = range(1, stop = w, length = N-w+1)
    for i = 1:N-w+1
        v = _maindir!(X[:, win], i == 1 ? d1 : @view(directions[:, i-1]), normalize_win)
        directions[:,i] .= v
        if center
            ii = (win[1] + win[end]) ÷ 2 + 1 # integer index from middle of window
            @views dominant_mode[i] = v'X[:, ii]
        else
            ii = i + wininds[i] - 1
            ib = floor(Int, ii)
            ia = ceil(Int, ii)
            α = ia - ii
            @views dominant_mode[i] = v'*((1-α) .* X[:, ib] .+ (α) .* X[:,ia]) # interpolate
        end
        win = win .+ 1
    end
    if cutoff_freq > 0
        directions = filtfilt(digitalfilter(Lowpass(cutoff_freq; fs=fs), Butterworth(4)), directions')'
    end
    for i = axes(directions, 2)
        @views directions[:,i] ./= norm(directions[:,i])
    end
    if pad
        dominant_mode = [zeros(w÷2); dominant_mode; zeros(w÷2-1)]
        directions = [repeat(directions[:,1], 1, w÷2) directions repeat(directions[:,1], 1, w÷2-1)]
        directions
    end
    dominant_mode, directions
end

function _maindir!(x, last, normalize_win)
    if normalize_win
        x .-= mean(x, dims = 1)
    end
    s = svd!(x, alg = LinearAlgebra.QRIteration()) # alg choice important for stability of singular vectors
    V = s.U
    v = V[1,:]
    if dot(last, v) < 0
        v .*= -1
    end
    v
end

# function sliding_pca2(X::AbstractMatrix{A}; dout = 1, l=2) where A
#     d, N = size(X)
#     λ   = zeros(A, dout)
#     V   = zeros(A, d, dout, N)
#     μ   = zeros(A, d)
#     xi  = zeros(A, d)
#     v   = zeros(A, d)
#     U   = zeros(A, dout, N)
#     n   = 0

#     for x in eachcol(X)

#         n += 1
#         @. μ += (x - μ) / n
#         xi = (n > 1) ? (x .- μ) : deepcopy(x)
#         f = (one(A)+l)/n
#         @views if n > 1
#             V[:, :, n] .= V[:, :, n-1]
#         end
#         @views @inbounds for i in 1:dout
#             if i == n
#                 λ[i] = norm(xi)
#                 @. V[:, i, n] = xi / (λ[i] + eps())
#                 break
#             end
#             v .= ((1-f) * λ[i]) .* V[:, i, n] .+ (f * dot(V[:, i, n], xi)) .* xi
#             λ[i] = norm(v)
#             @. V[:, i, n] = v/(λ[i] + eps())
#             xi .= xi .- dot(V[:, i, n], xi) .* V[:, i]
#         end
#         @views mul!(U[:, n], V[:, :, n]', xi)
#         U[:, n] .*= λ
#     end
#     U, V

# end


text:
module PropDicts

export PropDict

struct PropDict{V}
    _kvs::Dict{Symbol, V}
end
PropDict(; kwargs...) = PropDict(Dict{Symbol, Any}(kwargs...))
PropDict{T}(; kwargs...) where T = PropDict(Dict{Symbol, T}(kwargs...))

function Base.getproperty(pd::PropDict, x::Symbol)
    kvs = getfield(pd, :_kvs)
    return kvs[x]
end

# Base.getproperty(pd::PropDict, ::Val{:_kvs}) = getfield(pd, :_kvs)
# Base.setproperty!(pd::PropDict, ::Val{:_kvs}, v) = setfield!(pd, :_kvs, v)

function Base.setproperty!(pd::PropDict, x::Symbol, v)
    kvs = getfield(pd, :_kvs)
    return kvs[x] = v
end

end # module PropDicts

:_kvs isa Val{:_kvs}


text:
using Oceananigans
using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBoundary
using Printf
using Plots

solid(x, y, z) = z < 0

underlying_grid = RectilinearGrid(size=128, z=(-0.5, 0.5), topology=(Flat, Flat, Bounded))
grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBoundary(solid))

closure = IsotropicDiffusivity(κ = 1.0)

model_kwargs = (tracers=:c, buoyancy=nothing, velocities=PrescribedVelocityFields())

model = HydrostaticFreeSurfaceModel(; grid=underlying_grid, closure=closure, model_kwargs...)
immersed_model = HydrostaticFreeSurfaceModel(; grid=grid, closure=closure, model_kwargs...)

initial_temperature(x, y, z) = exp(-z^2 / 0.02)
[set!(m, c=initial_temperature) for m in (model, immersed_model)]

z = znodes(model.tracers.c)
c = view(interior(model.tracers.c), 1, 1, :)
c_immersed = view(interior(immersed_model.tracers.c), 1, 1, :)
c_plot = plot(c, z, linewidth = 2, label = "t = 0", xlabel = "Tracer concentration", ylabel = "z")
plot!(c_plot, c_immersed, z, linewidth = 2, label = "t = 0, immersed model", xlabel = "Tracer concentration", ylabel = "z")

diffusion_time_scale = model.grid.Δzᵃᵃᶜ^2 / model.closure.κ.c
stop_time = 100diffusion_time_scale

simulations = [simulation = Simulation(m, Δt = 1e-1 * diffusion_time_scale, stop_time = stop_time) for m in (model, immersed_model)]
[run!(sim) for sim in simulations]

plot!(c_plot, c, z, linewidth = 2, alpha = 0.6, linestyle = :dash, label = @sprintf("Ordinary model, t = %.3e", model.clock.time))
plot!(c_plot, c_immersed, z, linewidth = 3, alpha = 0.6, label = @sprintf("Immersed model, t = %.3e", immersed_model.clock.time))

display(c_plot)


text:
@testset "savefield.jl" begin

end


text:
abstract type AbstractDiagonalSymmetricMatrix{M,T,FS,L} <: DiagonalStructure{M,T,FS,L} end

mutable struct DiagonalSymmetricMatrix{M,T,FS,L} <: AbstractDiagonalSymmetricMatrix{M,T,FS,L}
    data::NTuple{L,T}
    function DiagonalSymmetricMatrix{M,T,FS,L}(::UndefInitializer) where {M,T,FS,L}
        new()
    end
    @generated function DiagonalSymmetricMatrix{M,T}(::UndefInitializer) where {M,T}
        :(DiagonalSymmetricMatrix{$M,$T,$M,$(abs2(M))}(undef))
    end
end

struct PointerDiagonalSymmetricMatrix{M,T,FS,L}; data::Ptr{T}; end

@generated Base.length(::AbstractDiagonalSymmetricMatrix{M}) where M = abs2(M)
@inline Base.size(::AbstractDiagonalSymmetricMatrix{M}) where M = (M,M)

@inline Base.pointer(A::PointerDiagonalSymmetricMatrix) = A.data
@inline function Base.pointer(A::DiagonalSymmetricMatrix{M,T}) where {M,T}
    Base.unsafe_convert(Ptr{T}, pointer_from_objref(A))
end

@inline function Base.getindex(A::AbstractDiagonalSymmetricMatrix{M,T,FS,L}, i::Integer) where {M,T,FS,L}
    @boundscheck i > L && throw(BoundsError())
    unsafe_load(pointer(A), i)
end

@inline function Base.getindex(A::AbstractDiagonalSymmetricMatrix{M,T,FS,L}, i::Integer, j::Integer) where {M,T,FS,L}
    i == j && return A[i]
    i, j = minmax(i, j)
    diagonal_structure_getindex(A, i, j)
end

@inline function Base.setindex!(A::AbstractDiagonalSymmetricMatrix{M,T,FS,L}, v, i::Integer) where {M,T,FS,L}
    @boundscheck i > L && throw(BoundsError())
    unsafe_store!(pointer(A), convert(T,v), i)
    v
end
@inline function Base.setindex!(A::AbstractDiagonalSymmetricMatrix{M,T,FS,L}, v, i::Integer, j::Integer) where {M,T,FS,L}
    i == j && return A[i] = v
    i, j = minmax(i, j)
    diagonal_structure_setindex!(A, v, i, j)
end

@inline Base.firstindex(A::AbstractDiagonalSymmetricMatrix) = 1
@inline Base.lastindex(A::AbstractDiagonalSymmetricMatrix{M,T,FS,L}) where {M,T,FS,L}= L

@inline function Base.similar(A::AbstractDiagonalSymmetricMatrix{M,T,FS,L}) where {M,T,FS,L}
    DiagonalSymmetricMatrix{M,T,FS,L}(undef)
end

id_symbol(A, i) = Symbol(A, :_, i)
id_symbol(A, i, j) = Symbol(A, :_, i, :_, j)
# using LinearAlgebra
function gen_ip_chol_quote!(qa,::Type{T},N::Int,L,U = :U, S = :S) where T

    chunksize = 32 ÷ sizeof(T)
    # @show qa

    for i ∈ 1:N
        Ui_i = id_symbol(U, i, i)
        U_i = id_symbol(S, i, i)
        if i == 1
            push!(qa, :($Ui_i = sqrt($U_i)))
        else
            nchunks, r = divrem(i-1, chunksize)
            if r == 0
                push!(qa, :($Ui_i = $U_i - +$([:( $(id_symbol(U, j, i)) * $(id_symbol(U, j, i))) for j ∈ 1:chunksize ]...)  ))
                for chunk ∈ 2:nchunks
                    push!(qa, :($Ui_i = $Ui_i - +$([:( $(id_symbol(U, j, i)) * $(id_symbol(U, j, i))) for j ∈ 1+(chunk-1)*chunksize:chunk*chunksize ]...)  ))
                end

            else
                push!(qa, :($Ui_i = $U_i - +$([:( $(id_symbol(U, j, i)) * $(id_symbol(U, j, i))) for j ∈ 1:r ]...)  ))
                for chunk ∈ 1:nchunks
                    push!(qa, :($Ui_i = $Ui_i - +$([:( $(id_symbol(U, j, i)) * $(id_symbol(U, j, i))) for j ∈ 1+r+(chunk-1)*chunksize:chunk*chunksize+r ]...)  ))
                end
            end
            # for j ∈ 1:i - 1
            #     Uj_i = id_symbol(U, lti+j)
            #     push!(qa, :($Ui_i -= $Uj_i*$Uj_i))
            # end
            push!(qa, :($Ui_i = sqrt($Ui_i)))
        end

        for j ∈ i+1:N
            Uj_i = id_symbol(U, i, j)
            U_i = id_symbol(S, i, j)
            if i == 1
                push!(qa, :($Uj_i = $U_i / $Ui_i) )
            else

                nchunks, r = divrem(i-1, chunksize)
                if r == 0
                    push!(qa, :($Uj_i = $U_i - +$([:( $(id_symbol(U, k, j)) * $(id_symbol(U, k, i))) for k ∈ 1:chunksize ]...)  ))
                    for chunk ∈ 2:nchunks
                        push!(qa, :($Uj_i= $Uj_i - +$([:( $(id_symbol(U, k, j)) * $(id_symbol(U, k, i))) for k ∈ 1+(chunk-1)*chunksize:chunk*chunksize ]...)  ))
                    end

                else
                    push!(qa, :($Uj_i = $U_i - +$([:( $(id_symbol(U, k, j)) * $(id_symbol(U, k, i))) for k ∈ 1:r ]...)  ))
                    for chunk ∈ 1:nchunks
                        push!(qa, :($Uj_i = $Uj_i - +$([:( $(id_symbol(U, k, j)) * $(id_symbol(U, k, i))) for k ∈ 1+r+(chunk-1)*chunksize:chunk*chunksize+r ]...)  ))
                    end
                end

                # push!(qa, :($Uj_i = $U_i))
                # for k ∈ 1:i - 1
                #     Ujk, Uik = id_symbol(U, ltj+k), id_symbol(U, lti+k)
                #     push!(qa, :($Uj_i -= $Ujk * $Uik))
                # end
                push!(qa, :($Uj_i /= $Ui_i) )
            end
        end
    end

    # push!(q.args, U)
    qa
end

@generated function chol!(U::DiagonalUpperTriangular{M,T,FS,L},
                            S::AbstractDiagonalSymmetricMatrix{M,T,FS,L}) where {M,T,FS,L}
    q = quote end
    ind = 0
    for k ∈ 0:M-1, m ∈ 1:M-k
        ind += 1
        push!(q.args, :( $(Symbol(:S_, m, :_, m+k)) = S[$ind] ) )
    end

    gen_ip_chol_quote!(q.args, T, M, L, :U, :S )

    ind = 0
    for k ∈ 0:M-1, m ∈ 1:M-k
        ind += 1
        push!(q.args, :( U[$ind] = $(Symbol(:U_, m, :_, m+k)) ) )
    end
    quote
        @fastmath @inbounds begin
            $q
            nothing
        end
    end
end

S = DiagonalSymmetricMatrix{8,Float64}(undef);
for i ∈ 1:8
   for j ∈ 1:i-1
       global S[i, j] = 1.0
   end
   S[i,i] = 5.0
end
U = DiagonalUpperTriangularMatrix{8,Float64}(undef);
chol!(U, S); U


text:
"""
    Checks if the good ideal composed of the given generators
    contains in the given groebner_basis ideal
"""
function check_ideal_inclusion(groebner_basis, generators)
    Is, yoverx, basepolyring, yoverxs, basepolyrings = groebner_ideal_to_singular(groebner_basis)
    Is = Ideal(yoverxs, Is)

    for f in generators
        f = idealize_element(f, basepolyring, basepolyrings, yoverx, yoverxs)

        fs = Ideal(yoverxs, f)

        if !iszero(Singular.reduce(fs, GroebnerBasis.f4(Is)))
            return false
        end
    end

    return true
end


text:
# Pretty printing
# ---------------
# Implement pshow(io::PrettyIO, t::T) for types T to be prettyprinted.

# pshow() should print no initial or final newlines.
# It should print using pshow and pprint.

# pprint(io, a, {b, c})
# prints a, then prints b and c indented.
# (The indenting is only noticable if there is a line break within {b, c})

# pprint(io, a, io->fun(io), b)
# invokes fun(io) on the current io. I e
# pprint(io, {a, io->fun(io), b})
# invokes fun() with an indented io.

#load("utils/utils.jl")

is_expr(ex, head::Symbol) = (isa(ex, Expr) && (ex.head == head))
function is_expr(ex, head::Symbol, nargs::Int)
    is_expr(ex, head) && length(ex.args) == nargs
end

pprintln(args...) = pprint(args..., '\n')
pprint(args...) = pprint(default_pretty(), args...)
pshow(args...)  = pshow(default_pretty(), args...)

# fallback for io::IO
pprint(io::IO, args...) = print(io, args...)

# -- PrettyIO -----------------------------------------------------------------

abstract PrettyIO

pprint_nowrap(io::PrettyIO, s::String) = (for c in s; pprint_nowrap(io, c);end)
function pprint(io::PrettyIO, s::String)
    n = strlen(s)
    if !str_fits_on_line(io, n)
        pprint(io, '\n')
    end
    for c in s; pprint(io, c); end
end
pprint(io::PrettyIO, arg::Any) = pshow(io, arg)
pprint(io::PrettyIO, args...) = foreach(arg->pprint(io, arg), args)

pshow(io::PrettyIO, arg::Any) = pprint(io, sshow(arg))

pprint(io::PrettyIO, v::Vector) = pprint(indent(io), v...)
pprint(io::PrettyIO, pprinter::Function) = pprinter(io)

comment(io::PrettyIO) = PrettyChild(io, ()->"# ")

indent(io::PrettyIO) = indent(io::PrettyIO, 4)
indent(io::PrettyIO, indent::Int) = (pre=" "^indent; PrettyChild(io, ()->pre))

# -- PrettyRoot ---------------------------------------------------------------

type PrettyRoot <: PrettyIO
    parent::IO
    width::Int

    currpos::Int
    autowrap::Bool

    function PrettyRoot(parent::IO, width::Int)
        if width < 1; error("width must be >= 1, got ", width); end
        new(parent, width, 0, false)
    end
end

function str_fits_on_line(io::PrettyRoot, n::Integer)
    (!io.autowrap) || (io.currpos+n <= io.width)
end

function pprint_nowrap(io::PrettyRoot, c::Char)
    if c=='\t'
        nsp::Int = (-io.currpos)&7
        #if nsp==0; nsp=8; end
        print(io.parent, " "^nsp)
        io.currpos += nsp
    else
        print(io.parent, c)
        io.currpos += 1
    end
    if c == '\n'
        io.currpos = 0
        io.autowrap = false
        return true
    end
    return false
end
function pprint(io::PrettyRoot, c::Char)
    if 2*io.currpos < io.width; io.autowrap=true; end

    if pprint_nowrap(io, c); return true; end
    if io.autowrap && (io.currpos >= io.width)
        return pprint_nowrap(io, '\n')
    end
    return false
end

default_pretty() = PrettyRoot(OUTPUT_STREAM, 80)

# -- PrettyChild --------------------------------------------------------------

type PrettyChild <: PrettyIO
    parent::PrettyIO
    newline_hook::Function

    function PrettyChild(parent::PrettyIO, newline_hook::Function)
        new(parent, newline_hook)
    end
end

str_fits_on_line(io::PrettyChild, n::Integer) = str_fits_on_line(io.parent, n)

pprint_nowrap(io::PrettyChild, c::Char) = pprint_nowrap(io.parent, c)
function pprint(io::PrettyChild, c::Char)
    newline = pprint(io.parent, c)::Bool
    if newline
        pprint_nowrap(io.parent, io.newline_hook())
    end
    return newline
end

# == Expr prettyprinting ======================================================

const doublecolon = @eval (:(x::Int)).head

## list printing
function pshow_comma_list(io::PrettyIO, args::Vector,
                          open::String, close::String)
    pshow_delim_list(io, args, open, ", ", close)
end
function pshow_delim_list(io::PrettyIO, args::Vector, open::String,
                          delim::String, close::String)
    pprint(io, {open,
                io->pshow_list_delim(io, args, delim)},
           close)
end
function pshow_list_delim(io::PrettyIO, args::Vector, delim::String)
    for (arg, k) in enumerate(args)
        pshow(io, arg)
        if k < length(args)
            pprint(io, delim)
        end
    end
end

## show the body of a :block
pshow_mainbody(io::PrettyIO, ex) = pshow(io, ex)
function pshow_mainbody(io::PrettyIO, ex::Expr)
    if ex.head == :block
        args = ex.args
        for (arg, k) in enumerate(args)
            if !is_expr(arg, :line)
                pprint(io, "\n")
            end
            pshow(io, arg)
        end
    else
        pshow(io, ex)
    end
end

## show arguments of a block, and then body
pshow_body(io::PrettyIO, body::Expr) = pshow_body(io, {}, body)
function pshow_body(io::PrettyIO, arg, body::Expr)
    pprint(io, {arg, io->pshow_mainbody(io, body) })
end
function pshow_body(io::PrettyIO, args::Vector, body::Expr)
    pprint(io, {
            io->pshow_comma_list(io, args, "", ""),
            io->pshow_mainbody(io, body)
        })
end

## show ex as if it were quoted
function pshow_quoted_expr(io::PrettyIO, sym::Symbol)
    if !is(sym,:(:)) && !is(sym,:(==))
        pprint(io, ":$sym")
    else
        pprint(io, ":($sym)")
    end
end
function pshow_quoted_expr(io::PrettyIO, ex::Expr)
    if ex.head == :block
        pprint(io, "quote ", io->pshow_body(io, ex), "\nend")
    else
        pprint(io, "quote(", {ex}, ")")
    end
end

## show an expr
function pshow(io::PrettyIO, ex::Expr)
    const infix = {:(=)=>"=", :(.)=>".", doublecolon=>"::", :(:)=>":",
                   :(->)=>"->", :(=>)=>"=>",
                   :(&&)=>" && ", :(||)=>" || "}
    const parentypes = {:call=>("(",")"), :ref=>("[","]"), :curly=>("{","}")}

    head = ex.head
    args = ex.args
    nargs = length(args)

    if has(infix, head) && nargs==2             # infix operations
#        pprint(io, "(",{args[1], infix[head], args[2]},")")
        pprint(io, {args[1], infix[head], args[2]})
    elseif has(parentypes, head) && nargs >= 1  # :call/:ref/:curly
        pprint(io, args[1])
        pshow_comma_list(io, args[2:end], parentypes[head]...)
    elseif (head == :comparison) && (nargs>=3 && isodd(nargs)) # :comparison
        pprint("(",{args},")")
    elseif ((contains([:return, :abstract, :const] , head) && nargs==1) ||
            contains([:local, :global], head))
        pshow_comma_list(io, args, string(head)*" ", "")
    elseif head == :typealias && nargs==2
        pshow_delim_list(io, args, string(head)*" ", " ", "")
    elseif (head == :quote) && (nargs==1)       # :quote
        pshow_quoted_expr(io, args[1])
    elseif (head == :line) && (1 <= nargs <= 2) # :line
        let io=comment(io)
            if nargs == 1
                linecomment = "line "*string(args[1])*": "
            else
                @assert nargs==2
#               linecomment = "line "*string(args[1])*", "*string(args[2])*": "
                linecomment = string(args[2])*", line "*string(args[1])*": "
            end
            if str_fits_on_line(io, strlen(linecomment)+13)
                pprint(io, "\t#  ", linecomment)
            else
                pprint(io, "\n", linecomment)
            end
        end
    elseif head == :if && nargs == 3  # if/else
        pprint(io,
            "if ", io->pshow_body(io, args[1], args[2]),
            "\nelse ", io->pshow_body(io, args[3]),
            "\nend")
    elseif head == :try && nargs == 3 # try[/catch]
        pprint(io, "try ", io->pshow_body(io, args[1]))
        if !(is(args[2], false) && is_expr(args[3], :block, 0))
            pprint(io, "\ncatch ", io->pshow_body(io, args[2], args[3]))
        end
        pprint(io, "\nend")
    elseif head == :let               # :let
        pprint(io, "let ",
            io->pshow_body(io, args[2:end], args[1]), "\nend")
    elseif head == :block
        pprint(io, "begin ", io->pshow_body(io, ex), "\nend")
    elseif contains([:for, :while, :function, :if, :type], head) && nargs == 2
        pprint(io, string(head), " ",
            io->pshow_body(io, args[1], args[2]), "\nend")
    else
        pprint(io, head)
        pshow_comma_list(indent(io), args, "(", ")")
    end
end

max_stars_repo_name: toivoh/julia-kernels

max_issues_repo_name: toivoh/julia-kernels

max_forks_repo_name: toivoh/julia-kernels