text:
"""
Calculate moments of Gaussian distribution `G = (λ / π)^(D / 2) * exp[-λ(c^2 + ξ^2)]`

* internality: `gauss_moments(prim::T) where {T<:AbstractArray{<:Real,1}}`
* no internality: `gauss_moments(prim::T, inK) where {T<:AbstractArray{<:Real,1}}`

"""
function gauss_moments(prim::T) where {T<:AbstractArray{<:Real,1}}

    if eltype(prim) <: Int
        MuL = OffsetArray(similar(prim, Float64, 7), 0:6)
    else
        MuL = OffsetArray(similar(prim, 7), 0:6)
    end
    MuR = similar(MuL)
    Mu = similar(MuL)

    MuL[0] = 0.5 * SpecialFunctions.erfc(-sqrt(prim[end]) * prim[2])
    MuL[1] = prim[2] * MuL[0] + 0.5 * exp(-prim[end] * prim[2]^2) / sqrt(π * prim[end])
    MuR[0] = 0.5 * SpecialFunctions.erfc(sqrt(prim[end]) * prim[2])
    MuR[1] = prim[2] * MuR[0] - 0.5 * exp(-prim[end] * prim[2]^2) / sqrt(π * prim[end])
    for i = 2:6
        MuL[i] = prim[2] * MuL[i-1] + 0.5 * (i - 1) * MuL[i-2] / prim[end]
        MuR[i] = prim[2] * MuR[i-1] + 0.5 * (i - 1) * MuR[i-2] / prim[end]
    end
    @. Mu = MuL + MuR

    if length(prim) == 3

        return Mu, MuL, MuR

    elseif length(prim) == 4

        Mv = similar(MuL)
        Mv[0] = 1.0
        Mv[1] = prim[3]
        for i = 2:6
            Mv[i] = prim[3] * Mv[i-1] + 0.5 * (i - 1) * Mv[i-2] / prim[end]
        end

        return Mu, Mv, MuL, MuR

    elseif length(prim) == 5

        Mv = similar(MuL)
        Mv[0] = 1.0
        Mv[1] = prim[3]
        for i = 2:6
            Mv[i] = prim[3] * Mv[i-1] + 0.5 * (i - 1) * Mv[i-2] / prim[end]
        end

        Mw = similar(MuL)
        Mw[0] = 1.0
        Mw[1] = prim[4]
        for i = 2:6
            Mw[i] = prim[4] * Mw[i-1] + 0.5 * (i - 1) * Mw[i-2] / prim[end]
        end

        return Mu, Mv, Mw, MuL, MuR

    end

end

# ------------------------------------------------------------
# A more general function
# deal with absent internality by setting inK = 0
# ------------------------------------------------------------
function gauss_moments(prim::T, inK) where {T<:AbstractArray{<:Real,1}}

    if eltype(prim) <: Int
        MuL = OffsetArray(similar(prim, Float64, 7), 0:6)
    else
        MuL = OffsetArray(similar(prim, 7), 0:6)
    end
    MuR = similar(MuL)
    Mu = similar(MuL)

    MuL[0] = 0.5 * SpecialFunctions.erfc(-sqrt(prim[end]) * prim[2])
    MuL[1] = prim[2] * MuL[0] + 0.5 * exp(-prim[end] * prim[2]^2) / sqrt(π * prim[end])
    MuR[0] = 0.5 * SpecialFunctions.erfc(sqrt(prim[end]) * prim[2])
    MuR[1] = prim[2] * MuR[0] - 0.5 * exp(-prim[end] * prim[2]^2) / sqrt(π * prim[end])
    for i = 2:6
        MuL[i] = prim[2] * MuL[i-1] + 0.5 * (i - 1) * MuL[i-2] / prim[end]
        MuR[i] = prim[2] * MuR[i-1] + 0.5 * (i - 1) * MuR[i-2] / prim[end]
    end
    @. Mu = MuL + MuR

    if length(prim) == 3

        Mxi = similar(MuL, 0:2)
        Mxi[0] = 1.0
        Mxi[1] = 0.5 * inK / prim[end]
        Mxi[2] = (inK^2 + 2.0 * inK) / (4.0 * prim[end]^2)

        return Mu, Mxi, MuL, MuR

    elseif length(prim) == 4

        Mv = similar(MuL)
        Mv[0] = 1.0
        Mv[1] = prim[3]
        for i = 2:6
            Mv[i] = prim[3] * Mv[i-1] + 0.5 * (i - 1) * Mv[i-2] / prim[end]
        end

        Mxi = similar(MuL, 0:2)
        Mxi[0] = 1.0
        Mxi[1] = 0.5 * inK / prim[end]
        Mxi[2] = (inK^2 + 2.0 * inK) / (4.0 * prim[end]^2)

        return Mu, Mv, Mxi, MuL, MuR

    elseif length(prim) == 5

        Mv = similar(MuL)
        Mv[0] = 1.0
        Mv[1] = prim[3]
        for i = 2:6
            Mv[i] = prim[3] * Mv[i-1] + 0.5 * (i - 1) * Mv[i-2] / prim[end]
        end

        Mw = similar(MuL)
        Mw[0] = 1.0
        Mw[1] = prim[4]
        for i = 2:6
            Mw[i] = prim[4] * Mw[i-1] + 0.5 * (i - 1) * Mw[i-2] / prim[end]
        end

        return Mu, Mv, Mw, MuL, MuR

    end

end

"""
    moments_conserve(Mu::OffsetArray{<:AbstractFloat,1}, alpha::Integer)

    moments_conserve(Mu::OffsetArray{<:Real,1}, Mxi::OffsetArray{<:Real,1},
        alpha::Integer, delta::Integer)

    moments_conserve(Mu::OffsetArray{<:Real,1}, Mv::OffsetArray{<:Real,1},
        Mw::OffsetArray{<:Real,1}, alpha::Integer, beta::Integer, delta::Integer)

Calculate conservative moments of particle distribution

"""
moments_conserve(Mu::T, alpha::I) where {T<:OffsetArray{<:AbstractFloat,1},I<:Integer} =
    Mu[alpha]

function moments_conserve(
    Mu::T,
    Mxi::T,
    alpha::I,
    delta::I,
) where {T<:OffsetArray{<:AbstractFloat,1},I<:Integer}

    uv = similar(Mu, 3)
    uv[1] = Mu[alpha] * Mxi[delta÷2]
    uv[2] = Mu[alpha+1] * Mxi[delta÷2]
    uv[3] = 0.5 * (Mu[alpha+2] * Mxi[delta÷2] + Mu[alpha] * Mxi[(delta+2)÷2])

    return uv

end

function moments_conserve(
    Mu::T,
    Mv::T,
    Mw::T,
    alpha::I,
    beta::I,
    delta::I,
) where {T<:OffsetArray{<:AbstractFloat,1},I<:Integer}

    if length(Mw) == 3 # internal motion
        uv = similar(Mu, 4)
        uv[1] = Mu[alpha] * Mv[beta] * Mw[delta÷2]
        uv[2] = Mu[alpha+1] * Mv[beta] * Mw[delta÷2]
        uv[3] = Mu[alpha] * Mv[beta+1] * Mw[delta÷2]
        uv[4] =
            0.5 * (
                Mu[alpha+2] * Mv[beta] * Mw[delta÷2] +
                Mu[alpha] * Mv[beta+2] * Mw[delta÷2] +
                Mu[alpha] * Mv[beta] * Mw[(delta+2)÷2]
            )
    else
        uv = similar(Mu, 5)
        uv[1] = Mu[alpha] * Mv[beta] * Mw[delta]
        uv[2] = Mu[alpha+1] * Mv[beta] * Mw[delta]
        uv[3] = Mu[alpha] * Mv[beta+1] * Mw[delta]
        uv[4] = Mu[alpha] * Mv[beta] * Mw[delta+1]
        uv[5] =
            0.5 * (
                Mu[alpha+2] * Mv[beta] * Mw[delta] +
                Mu[alpha] * Mv[beta+2] * Mw[delta] +
                Mu[alpha] * Mv[beta] * Mw[delta+2]
            )
    end

    return uv

end

# ------------------------------------------------------------
# Discrete moments of conservative variables
# ------------------------------------------------------------
#--- 1F1V ---#
function moments_conserve(
    f::X,
    u::T,
    ω::T,
) where {X<:AbstractArray{<:AbstractFloat,1},T<:AbstractArray{<:AbstractFloat,1}}
    w = similar(f, 3)
    w[1] = discrete_moments(f, u, ω, 0)
    w[2] = discrete_moments(f, u, ω, 1)
    w[3] = 0.5 * discrete_moments(f, u, ω, 2)

    return w
end

#--- 2F1V ---#
function moments_conserve(
    h::X,
    b::X,
    u::T,
    ω::T,
) where {X<:AbstractArray{<:AbstractFloat,1},T<:AbstractArray{<:AbstractFloat,1}}
    w = similar(h, 3)
    w[1] = discrete_moments(h, u, ω, 0)
    w[2] = discrete_moments(h, u, ω, 1)
    w[3] = 0.5 * (discrete_moments(h, u, ω, 2) + discrete_moments(b, u, ω, 0))

    return w
end

#--- 1F2V ---#
function moments_conserve(
    f::X,
    u::T,
    v::T,
    ω::T,
) where {X<:AbstractArray{<:AbstractFloat,2},T<:AbstractArray{<:AbstractFloat,2}}
    w = similar(f, 4)
    w[1] = discrete_moments(f, u, ω, 0)
    w[2] = discrete_moments(f, u, ω, 1)
    w[3] = discrete_moments(f, v, ω, 1)
    w[4] = 0.5 * (discrete_moments(f, u, ω, 2) + discrete_moments(f, v, ω, 2))

    return w
end

#--- 2F2V ---#
function moments_conserve(
    h::X,
    b::X,
    u::T,
    v::T,
    ω::T,
) where {X<:AbstractArray{<:AbstractFloat,2},T<:AbstractArray{<:AbstractFloat,2}}
    w = similar(h, 4)
    w[1] = discrete_moments(h, u, ω, 0)
    w[2] = discrete_moments(h, u, ω, 1)
    w[3] = discrete_moments(h, v, ω, 1)
    w[4] =
        0.5 * (
            discrete_moments(h, u, ω, 2) +
            discrete_moments(h, v, ω, 2) +
            discrete_moments(b, u, ω, 0)
        )

    return w
end

#--- 3F2V ---#
function moments_conserve(
    h0::X,
    h1::X,
    h2::X,
    u::T,
    v::T,
    ω::T,
) where {X<:AbstractArray{<:AbstractFloat,2},T<:AbstractArray{<:AbstractFloat,2}}
    w = similar(h0, 5)
    w[1] = discrete_moments(h0, u, ω, 0)
    w[2] = discrete_moments(h0, u, ω, 1)
    w[3] = discrete_moments(h0, v, ω, 1)
    w[4] = discrete_moments(h1, u, ω, 0)
    w[5] =
        0.5 * (
            discrete_moments(h0, u, ω, 2) +
            discrete_moments(h0, v, ω, 2) +
            discrete_moments(h2, u, ω, 0)
        )

    return w
end

#--- 1F3V ---#
function moments_conserve(
    f::X,
    u::T,
    v::T,
    w::T,
    ω::T,
) where {X<:AbstractArray{<:AbstractFloat,3},T<:AbstractArray{<:AbstractFloat,3}}
    moments = similar(f, 5)

    moments[1] = discrete_moments(f, u, ω, 0)
    moments[2] = discrete_moments(f, u, ω, 1)
    moments[3] = discrete_moments(f, v, ω, 1)
    moments[4] = discrete_moments(f, w, ω, 1)
    moments[5] =
        0.5 * (
            discrete_moments(f, u, ω, 2) +
            discrete_moments(f, v, ω, 2) +
            discrete_moments(f, w, ω, 2)
        )

    return moments
end

#--- 4F1V ---#
function moments_conserve(
    h0::X,
    h1::X,
    h2::X,
    h3::X,
    u::T,
    ω::T,
) where {X<:AbstractArray{<:AbstractFloat,1},T<:AbstractArray{<:AbstractFloat,1}}
    moments = similar(h0, 5)

    moments[1] = discrete_moments(h0, u, ω, 0)
    moments[2] = discrete_moments(h0, u, ω, 1)
    moments[3] = discrete_moments(h1, u, ω, 0)
    moments[4] = discrete_moments(h2, u, ω, 0)
    moments[5] = 0.5 * discrete_moments(h0, u, ω, 2) + 0.5 * discrete_moments(h3, u, ω, 0)

    return moments
end

"""
Calculate conservative moments of diatomic particle distribution

- 1D: `diatomic_moments_conserve(
    h::X,
    b::X,
    r::X,
    u::T,
    ω::T,
) where {X<:AbstractArray{<:AbstractFloat,1},T<:AbstractArray{<:AbstractFloat,1}}`
- 2D: `diatomic_moments_conserve(
    h0::X,
    h1::X,
    h2::X,
    u::T,
    v::T,
    ω::T,
) where {X<:AbstractArray{<:AbstractFloat,2},T<:AbstractArray{<:AbstractFloat,2}}`

"""
function diatomic_moments_conserve(
    h::X,
    b::X,
    r::X,
    u::T,
    ω::T,
) where {X<:AbstractArray{<:AbstractFloat,1},T<:AbstractArray{<:AbstractFloat,1}}
    w = similar(h, 4)
    w[1] = discrete_moments(h, u, ω, 0)
    w[2] = discrete_moments(h, u, ω, 1)
    w[3] =
        0.5 * (
            discrete_moments(h, u, ω, 2) +
            discrete_moments(b, u, ω, 0) +
            discrete_moments(r, u, ω, 0)
        )
    w[4] = 0.5 * discrete_moments(r, u, ω, 0)

    return w
end

#--- 3F2V ---#
function diatomic_moments_conserve(
    h0::X,
    h1::X,
    h2::X,
    u::T,
    v::T,
    ω::T,
) where {X<:AbstractArray{<:AbstractFloat,2},T<:AbstractArray{<:AbstractFloat,2}}
    w = similar(h0, 5)
    w[1] = discrete_moments(h0, u, ω, 0)
    w[2] = discrete_moments(h0, u, ω, 1)
    w[3] = discrete_moments(h0, v, ω, 1)
    w[4] =
        0.5 * (
            discrete_moments(h0, u, ω, 2) +
            discrete_moments(h0, v, ω, 2) +
            discrete_moments(h1, u, ω, 0) +
            discrete_moments(h2, u, ω, 0)
        )
    w[5] = 0.5 * discrete_moments(h2, u, ω, 0)

    return w
end

"""
Calculate slope-related conservative moments
`a = a1 + u * a2 + 0.5 * u^2 * a3`

"""
moments_conserve_slope(a, Mu::T, alpha::I) where {T<:OffsetArray{<:Real,1},I<:Int} =
    a * moments_conserve(Mu, alpha)

moments_conserve_slope(
    a::X,
    Mu::Y,
    Mxi::Y,
    alpha::I,
) where {X<:AbstractArray{<:Real,1},Y<:OffsetArray{<:Real,1},I<:Int} =
    a[1] .* moments_conserve(Mu, Mxi, alpha + 0, 0) .+
    a[2] .* moments_conserve(Mu, Mxi, alpha + 1, 0) .+
    0.5 * a[3] .* moments_conserve(Mu, Mxi, alpha + 2, 0) .+
    0.5 * a[3] .* moments_conserve(Mu, Mxi, alpha + 0, 2)

function moments_conserve_slope(
    a::X,
    Mu::Y,
    Mv::Y,
    Mxi::Y,
    alpha::I,
    beta::I,
) where {X<:AbstractArray{<:Real,1},Y<:OffsetArray{<:Real,1},I<:Int}

    return a[1] .* moments_conserve(Mu, Mv, Mxi, alpha + 0, beta + 0, 0) .+
           a[2] .* moments_conserve(Mu, Mv, Mxi, alpha + 1, beta + 0, 0) .+
           a[3] .* moments_conserve(Mu, Mv, Mxi, alpha + 0, beta + 1, 0) .+
           0.5 * a[4] .* moments_conserve(Mu, Mv, Mxi, alpha + 2, beta + 0, 0) .+
           0.5 * a[4] .* moments_conserve(Mu, Mv, Mxi, alpha + 0, beta + 2, 0) .+
           0.5 * a[4] .* moments_conserve(Mu, Mv, Mxi, alpha + 0, beta + 0, 2)

end

function moments_conserve_slope(
    a::X,
    Mu::Y,
    Mv::Y,
    Mw::Y,
    alpha::I,
    beta::I,
    delta::I,
) where {X<:AbstractArray{<:Real,1},Y<:OffsetArray{<:Real,1},I<:Int}

    return a[1] .* moments_conserve(Mu, Mv, Mw, alpha + 0, beta + 0, delta + 0) .+
           a[2] .* moments_conserve(Mu, Mv, Mw, alpha + 1, beta + 0, delta + 0) .+
           a[3] .* moments_conserve(Mu, Mv, Mw, alpha + 0, beta + 1, delta + 0) .+
           a[4] .* moments_conserve(Mu, Mv, Mw, alpha + 0, beta + 0, delta + 1) .+
           0.5 * a[5] .* moments_conserve(Mu, Mv, Mw, alpha + 2, beta + 0, delta + 0) .+
           0.5 * a[5] .* moments_conserve(Mu, Mv, Mw, alpha + 0, beta + 2, delta + 0) .+
           0.5 * a[5] .* moments_conserve(Mu, Mv, Mw, alpha + 0, beta + 0, delta + 2)

end

"""
    * direct quadrature
    discrete_moments(f, ω)

    * velocity moments
    discrete_moments(f, u, ω, n)

Discrete moments of particle distribution

"""
discrete_moments(f, ω) = sum(@. ω * f)

discrete_moments(f, u, ω, n) = sum(@. ω * u^n * f)

"""
    pressure(f, prim, u, ω)
    pressure(h, b, prim, u, ω, K)
    pressure(h, b, prim, u, v, ω, K)

Calculate pressure from particle distribution function

"""
pressure(f, prim, u, ω) = sum(@. ω * (u - prim[2])^2 * f)

pressure(h, b, prim, u, ω, K) =
    (sum(@. ω * (u - prim[2])^2 * h) + sum(@. ω * b)) / (K + 1.0)

pressure(h, b, prim, u, v, ω, K) =
    (sum(@. ω * ((u - prim[2])^2 + (v - prim[3])^2) * h) + sum(@. ω * b)) / (K + 2.0)

"""
Calculate stress tensor from particle distribution function

"""
stress(f, prim, u, ω) = sum(@. ω * (u - prim[2]) * (u - prim[2]) * f)

function stress(f, prim, u, v, ω)
    P = similar(prim, 2, 2)

    P[1, 1] = sum(@. ω * (u - prim[2]) * (u - prim[2]) * f)
    P[1, 2] = sum(@. ω * (u - prim[2]) * (v - prim[3]) * f)
    P[2, 1] = P[1, 2]
    P[1, 2] = sum(@. ω * (v - prim[3]) * (v - prim[3]) * f)

    return P
end

"""
    heat_flux(f, prim, u, ω)
    heat_flux(h, b, prim, u, ω)
    heat_flux(h, b, r, prim, u, ω)
    heat_flux(f, prim, u, v, ω)
    heat_flux(h, b, prim, u, v, ω)
    heat_flux(h, b, r, prim, u, v, ω)
    heat_flux(f, prim, u, v, w, ω)

Calculate heat flux from particle distribution function

Multiple dispatch doesn't consider unstructured multi-dimensional velocity space.
In that case a new method needs to be defined.

"""
heat_flux(f, prim, u, ω) = 0.5 * sum(@. ω * (u - prim[2]) * (u - prim[2])^2 * f) # 1F1V

#--- 2F1V ---#
heat_flux(
    h::X,
    b::X,
    prim::Y,
    u::Z,
    ω::Z,
) where {
    X<:AbstractArray{<:AbstractFloat,1},
    Y<:AbstractArray{<:Real,1},
    Z<:AbstractArray{<:AbstractFloat,1},
} = 0.5 * (sum(@. ω * (u - prim[2]) * (u - prim[2])^2 * h) + sum(@. ω * (u - prim[2]) * b))

#--- 3F1V Rykov ---#
function heat_flux(
    h::X,
    b::X,
    r::X,
    prim::Y,
    u::Z,
    ω::Z,
) where {
    X<:AbstractArray{<:AbstractFloat,1},
    Y<:AbstractArray{<:Real,1},
    Z<:AbstractArray{<:AbstractFloat,1},
}

    q = similar(h, 2)

    q[1] =
        0.5 *
        (sum(@. ω * (u - prim[2]) * (u - prim[2])^2 * h) + sum(@. ω * (u - prim[2]) * b))
    q[2] = 0.5 * (sum(@. ω * (u - prim[2]) * r))

    return q

end

#--- 1F2V ---#
function heat_flux(
    h::X,
    prim::Y,
    u::Z,
    v::Z,
    ω::Z,
) where {
    X<:AbstractArray{<:AbstractFloat,2},
    Y<:AbstractArray{<:Real,1},
    Z<:AbstractArray{<:AbstractFloat,2},
}

    q = similar(h, 2)
    q[1] = 0.5 * sum(@. ω * (u - prim[2]) * ((u - prim[2])^2 + (v - prim[3])^2) * h)
    q[2] = 0.5 * sum(@. ω * (v - prim[3]) * ((u - prim[2])^2 + (v - prim[3])^2) * h)

    return q

end

#--- 2F2V ---#
function heat_flux(
    h::X,
    b::X,
    prim::Y,
    u::Z,
    v::Z,
    ω::Z,
) where {
    X<:AbstractArray{<:AbstractFloat,2},
    Y<:AbstractArray{<:Real,1},
    Z<:AbstractArray{<:AbstractFloat,2},
}

    q = similar(h, 2)

    q[1] =
        0.5 * (
            sum(@. ω * (u - prim[2]) * ((u - prim[2])^2 + (v - prim[3])^2) * h) +
            sum(@. ω * (u - prim[2]) * b)
        )
    q[2] =
        0.5 * (
            sum(@. ω * (v - prim[3]) * ((u - prim[2])^2 + (v - prim[3])^2) * h) +
            sum(@. ω * (v - prim[3]) * b)
        )

    return q

end

#--- 3F2V Rykov ---#
function heat_flux(
    h::X,
    b::X,
    r::X,
    prim::Y,
    u::Z,
    v::Z,
    ω::Z,
) where {
    X<:AbstractArray{<:AbstractFloat,2},
    Y<:AbstractArray{<:Real,1},
    Z<:AbstractArray{<:AbstractFloat,2},
}

    q = similar(h, 4)

    q[1] =
        0.5 * (
            sum(@. ω * (u - prim[2]) * ((u - prim[2])^2 + (v - prim[3])^2) * h) +
            sum(@. ω * (u - prim[2]) * b)
        )
    q[2] =
        0.5 * (
            sum(@. ω * (v - prim[3]) * ((u - prim[2])^2 + (v - prim[3])^2) * h) +
            sum(@. ω * (v - prim[3]) * b)
        )
    q[3] = 0.5 * sum(@. ω * (u - prim[2]) * r)
    q[4] = 0.5 * sum(@. ω * (v - prim[3]) * r)

    return q

end

#--- 1F3V ---#
function heat_flux(
    f::X,
    prim::Y,
    u::Z,
    v::Z,
    w::Z,
    ω::Z,
) where {
    X<:AbstractArray{<:AbstractFloat,3},
    Y<:AbstractArray{<:Real,1},
    Z<:AbstractArray{<:AbstractFloat,3},
}

    q = similar(f, 3)

    q[1] =
        0.5 * sum(
            @. ω *
               (u - prim[2]) *
               ((u - prim[2])^2 + (v - prim[3])^2 + (w - prim[4])^2) *
               f
        )
    q[2] =
        0.5 * sum(
            @. ω *
               (v - prim[3]) *
               ((u - prim[2])^2 + (v - prim[3])^2 + (w - prim[4])^2) *
               f
        )
    q[3] =
        0.5 * sum(
            @. ω *
               (w - prim[4]) *
               ((u - prim[2])^2 + (v - prim[3])^2 + (w - prim[4])^2) *
               f
        )

    return q

end


text:
function stringcommand(str::String)

ind = findall(" ", str);
n = length(ind);
cmdstr = Array{String,1}(undef,n+1);
cmdstr[1] = str[1:(ind[1][1]-1)];
for i = 2:n
    cmdstr[i] = str[(ind[i-1][1]+1):(ind[i][1]-1)];
end
i = n+1;
cmdstr[i] = str[(ind[i-1][1]+1):end];

return Cmd(cmdstr)

end

function cmakecompile(pde, numpde=1)

cdir = pwd(); ii = findlast("Exasim", cdir);
ExasimPath = cdir[1:ii[end]];
DataPath = ".." * cdir[(ii[end]+1):end];
versiondir = ExasimPath * "/" * pde.version;
appdriverdir = versiondir  * "/Kernel/AppDriver/";

cp(appdriverdir * "opuApp.cpp", cdir * "/app/opuApp.cpp", force=true);
cp(appdriverdir * "cpuApp.cpp", cdir * "/app/cpuApp.cpp", force=true);
cp(appdriverdir * "gpuApp.cu",  cdir * "/app/gpuApp.cu", force=true);
#cp("app/*.cpp", ExasimPath * "/Applications/App", force=true);
#cp("app/*.cu",  ExasimPath * "/Applications/App", force=true);
#run(stringcommand("cp -f " * "app/*.cpp " * ExasimPath * "/Applications/App"));
#run(stringcommand("cp -f " * "app/*.cu " * ExasimPath * "/Applications/App"));

if lowercase(pde.version) == lowercase("Version0.1")
    version = " -D EXASIM_VERSION01=ON ";
    pdenum = "";
elseif lowercase(pde.version) == lowercase("Version0.2")
    version = " -D EXASIM_VERSION02=ON ";
    pdenum = string(numpde) * " ";
elseif lowercase(pde.version) == lowercase("Version0.3")
    version = " -D EXASIM_VERSION03=ON ";
    pdenum = string(numpde) * " ";
end
version = version * "-D EXASIM_APPDIR=" * DataPath * "/app ";

# create exec folder if it does not exist
bindir = "exec";
cd(ExasimPath);
if !isdir(bindir)
    mkdir(bindir);
end
cd(ExasimPath * "/" * bindir);

if isfile("CMakeCache.txt")
    rm("CMakeCache.txt");
end
if pde.platform == "gpu"
    if pde.buildexec==1
        run(stringcommand("cmake -D EXASIM_APP=ON -D EXASIM_CORES=ON -D EXASIM_EXECUTABLES=ON -D EXASIM_CUDA=ON " * version * "../Installation"));
    else
        if exist("libcpuCore.a", "file") && exist("libgpuCore.a", "file") && exist("gpuEXASIM", "file")
            run(stringcommand("cmake -D EXASIM_APP=ON -D EXASIM_CUDA=ON " * version * "../Installation"));
        else
            run(stringcommand("cmake -D EXASIM_APP=ON -D EXASIM_CORES=ON -D EXASIM_EXECUTABLES=ON -D EXASIM_CUDA=ON " * version * "../Installation"));
        end
    end
elseif pde.platform == "cpu"
    if pde.buildexec==1
        run(stringcommand("cmake -D EXASIM_APP=ON -D EXASIM_CORES=ON -D EXASIM_EXECUTABLES=ON " * version * "../Installation"));
    else
        if exist("libcpuCore.a", "file") && exist("cpuEXASIM", "file")
           run(stringcommand("cmake -D EXASIM_APP=ON " * version * "../Installation"));
        else
           run(stringcommand("cmake -D EXASIM_APP=ON -D EXASIM_CORES=ON -D EXASIM_EXECUTABLES=ON " * version * "../Installation"))s;
        end
    end
end
run(stringcommand("cmake --build ."));

mpirun = pde.mpirun;
if pde.platform == "cpu"
    if pde.mpiprocs==1
        str = "./cpuEXASIM " * pdenum * DataPath * "/datain/ " * DataPath * "/dataout/out";
    else
        str = mpirun * " -np " * string(pde.mpiprocs) * " ./cpumpiEXASIM " * pdenum * DataPath * "/datain/ " * DataPath * "/dataout/out";
    end
    run(stringcommand(str));
elseif pde.platform == "gpu"
    if pde.mpiprocs==1
        str = "./gpuEXASIM " * pdenum * DataPath * "/datain/ " * DataPath * "/dataout/out";
    else
        str = mpirun * " -np " * string(pde.mpiprocs) * " ./app/gpumpiEXASIM " * pdenum * DataPath * "/datain/ " * DataPath * "/dataout/out";
    end
    run(stringcommand(str));
end

cd(cdir);

return str;

end

ext: jl

max_stars_repo_name: rloekvh/Exasim

max_issues_repo_name: rloekvh/Exasim

max_forks_repo_name: rloekvh/Exasim


text:
##############################################################################################
module SlurmExperimentDriver
##############################################################################################
##############################################################################################
# imports and exports
using Random, Distributions
using Debugger
using Distributed
using LinearAlgebra
using JLD2
using DeSolvers
using L96
using EnsembleKalmanSchemes
using FilterExps
using SmootherExps

##############################################################################################
# Definition of globals

path = joinpath(@__DIR__, "../data/")

##############################################################################################
## Experiments to run as a single function call
##############################################################################################
# Filters
##############################################################################################

function filter_state_exp(j)
    f = load(path * "input_data/filter_state_input_args.jld2")
    args = f["experiments"][j]
    filter_state(args)
end

##############################################################################################
# Classic smoothers
##############################################################################################

function classic_smoother_state_exp(j)
    f = load(path * "input_data/classic_state_smoother_input_args.jld2")
    args = f["experiments"][j]
    classic_state(args)
end

function classic_smoother_param_exp(j)
    f = load(path * "input_data/classic_param_smoother_input_args.jld2")
    args = f["experiments"][j]
    classic_param(args)
end

##############################################################################################
# Single-iteration smoothers
##############################################################################################

function single_iteration_smoother_state_exp(j)
    f = load(path * "input_data/single_iteration_state_smoother_input_args.jld2")
    args = f["experiments"][j]
    single_iteration_state(args)
end

function single_iteration_smoother_param_exp(j)
    f = load(path * "input_data/single_iteration_param_smoother_input_args.jld2")
    args = f["experiments"][j]
    single_iteration_param(args)
end

##############################################################################################
# Iterative smoothers
##############################################################################################

function iterative_smoother_state_exp(j)
    f = load("../data/input_data/iterative_state_smoother_input_args.jld2")
    args = f["experiments"][j]
    iterative_state(args)
end

function iterative_smoother_param_exp(j)
    f = load("../data/input_data/iterative_param_smoother_input_args.jld2")
    args = f["experiments"][j]
    iterative_param(args)
end

##############################################################################################
# Run experiments
##############################################################################################

t = parse(Int64, ARGS[1])
s = ARGS[2]

if s == "filter_state"
    filter_state_exp(t)
elseif s == "classic_smoother_state"
    classic_smoother_state_exp(t)
elseif s == "classic_smoother_param"
    classic_smoother_param_exp(t)
elseif s == "single_iteration_smoother_state"
    single_iteration_smoother_state_exp(t)
elseif s == "single_iteration_smoother_param"
    single_iteration_smoother_param_exp(t)
elseif s == "iterative_smoother_state"
    iterative_smoother_state_exp(t)
elseif s == "iterative_smoother_param"
    iterative_smoother_param_exp(t)
end

##############################################################################################
# end module

end

max_stars_repo_name: cgrudz/DataAssimilationBenchmarks

max_issues_repo_name: cgrudz/DataAssimilationBenchmarks


text:
using SafeTestsets

    @time @safetestset "Constant Matrices Tests" begin include("matrices.jl") end
    @time @safetestset "SB Model tests" begin include("sbmodel.jl") end
    @time @safetestset "Ehrenfest dynamics test" begin include("ehrenfest.jl") end
    @time @safetestset "ECEID dynamics test" begin include("ECEID.jl") end
    @time @safetestset "CEID dynamics test" begin include("CEID.jl") end
    @time @safetestset "ECEID long test" begin include("ECEID_long.jl") end
    @time @safetestset "CEID long test" begin include("CEID_long.jl") end


text:
using MarketData
using Test

@testset "remote" begin
    @testset "FRED" begin
        ta = fred()
        @test ta |> timestamp |> length > 100
    end

    @testset "Yahoo" begin
        # ta = yahoo()
        # @test length(ta.timestamp) > 100
        @test_broken 1==2
    end

end


text:
using BinDeps

@BinDeps.setup

basedir = joinpath(Pkg.dir("DrakeVisualizer"), "deps")
director_version = v"0.1.0"

@static if is_linux()
    deps = [
        python = library_dependency("python", aliases=["libpython2.7.so", "libpython3.2.so", "libpython3.3.so", "libpython3.4.so", "libpython3.5.so", "libpython3.6.so", "libpython3.7.so"])
        python_vtk = library_dependency("vtkCommon", aliases=["libvtkCommon.so", "libvtkCommon.so.5.8", "libvtkCommon.so.5.10"], depends=[python],
            validate=(name, handle) -> begin
                isfile(replace(name, r"vtkCommon", "vtkCommonPythonD", 1))
            end)
        director = library_dependency("ddApp", aliases=["libddApp"], depends=[python_vtk, python])
    ]

    # The vtkPython libraries all have undeclared dependencies on libpython2.7,
    # so they cannot be dlopen()ed without missing symbol errors. As a result,
    # we can't use the regular library_dependency mechanism to look for vtk5
    # and python-vtk. Instead, we combined both dependencies into "python_vtk"
    # and make one build rule to apt-get install all the vtk-related packages.
    provides(SimpleBuild,
        () -> run(`sudo apt-get install libvtk5-qt4-dev python-vtk`),
        python_vtk)
    provides(AptGet, Dict("python2.7" => python))
    provides(BuildProcess, (@build_steps begin
        FileDownloader("http://people.csail.mit.edu/patmarion/software/director/releases/director-$(director_version)-linux.tar.gz",
                       joinpath(basedir, "downloads", "director.tar.gz"))
        CreateDirectory(joinpath(basedir, "usr"))
        (`tar xzf $(joinpath(basedir, "downloads", "director.tar.gz")) --directory=usr --strip-components=1`)
    end), director)
elseif is_apple()
    deps = [
        director = library_dependency("ddApp", aliases=["libddApp"])
    ]
    provides(BuildProcess, (@build_steps begin
        FileDownloader("http://people.csail.mit.edu/patmarion/software/director/releases/director-$(director_version)-mac.tar.gz",
                       joinpath(basedir, "downloads", "director.tar.gz"))
        CreateDirectory(joinpath(basedir, "usr"))
        (`tar xzf $(joinpath(basedir, "downloads", "director.tar.gz")) --directory=usr --strip-components=1`)
    end), director)

end

@BinDeps.install Dict(:ddApp => :libddApp)


text:
"""
    SortedDataFrame(df::DataFrame, time::Symbol = :Time, grouping::Symbol = :Name, value::Symbol = :Value)

This struct wraps a `DataFrame`. In the constructor function for the dataframe
we presort the data and create a mapping dict so that it is fast to subset the
DataFrame by the group.

For the constructor pass in the dataframe, name of time column, name of grouping
 column and name of value column to the constructor.
### Inputs
* `df` - The tick data
* `time` - The column of the data representing time.
* `grouping` - The column of the data representing the asset name
* `value` - The column of the data representing price/logprice/etc.
* `time_period_per_unit` - The period that one unit (in the time column) corresponds to.

### Returns
* A `SortedDataFrame`.
"""
struct SortedDataFrame
    df::DataFrame
    time::Symbol
    grouping::Symbol
    value::Symbol
    groupingrows::Dict{Symbol,Array{Int64,1}}
    time_period_per_unit::Dates.Period
    function SortedDataFrame(df::DataFrame, time::Symbol, grouping::Symbol,
                             value::Symbol, time_period_per_unit::Dates.Period)
       df2 = sort(df, time)
       dic = Dict{Symbol,Vector{Int64}}()
       for g in unique(df2[:,grouping])
          dic[g] = findall(df2[:,grouping] .== g)
       end
       if typeof(time_period_per_unit) in [Year, Quarter, Month]
           error("time_period_per_unit cannot be a Year, Quarter or Month as it is undefined how many seconds/days/etc exist per Year, Quarter or Month.")
       end
       return new(df2, time, grouping, value, dic, time_period_per_unit)
    end
    function SortedDataFrame(df::DataFrame, timevar::Symbol, groupingvar::Symbol, valuevar::Symbol,
                             dic::Dict{Symbol,Vector{I}}, vol_unit::Dates.Period) where I<:Integer
        return new(df, timevar, groupingvar, valuevar, dic, vol_unit)
    end
end

function time_period_ratio(neww::Dates.Period, oldd::Dates.Period)
    return Nanosecond(neww) / Nanosecond(oldd)
end
function safe_multiply_period(scalar::Real, neww::Dates.Period)
    vall = Nanosecond(neww).value * scalar
    return Nanosecond(floor(vall))
end

import Base.show, Base.print

"""
    show(sdf::SortedDataFrame, number_of_rows = 10)
Show a SortedDataFrame with a set number of rows.
### Inputs
* `sdf` - The `SortedDataFrame` to show.
* `number_of_rows` - The number of rows to show.
"""
function show(sdf::SortedDataFrame, number_of_rows = 10)
    println()
    println("SortedDataFrame with " , nrow(sdf.df), " rows."  )
    show(stdout, sdf.df[1:number_of_rows,[sdf.time, sdf.grouping, sdf.value]] ; summary = false)
    println()
end
"""
    print(sdf::SortedDataFrame, number_of_rows = 10)
Show a SortedDataFrame with a set number of rows.
### Inputs
* `sdf` - The `SortedDataFrame` to print.
* `number_of_rows` - The number of rows to print.
"""
function print(sdf::SortedDataFrame, number_of_rows = 10)
    show(sdf, number_of_rows)
end

"""
    get_assets(ts::SortedDataFrame, obs_to_include::Integer = 10)

This returns a vector of all of the assets in the `SortedDataFrame` with at least
some number of observations (10 by default).
### Inputs
* `ts` - The tick data.
* `obs_to_include` - An integer for the minimum number of ticks in `ts` needed for the function to include that asset.
### Returns
* A `Vector{Symbol}` with each asset.
"""
function get_assets(ts::SortedDataFrame, obs_to_include::Integer = 10)
    all_assets = unique(ts.df[:,ts.grouping])
    assets = Array{Symbol,1}(undef,0)
    for a in all_assets
        cond1 = length(ts.groupingrows[a]) >= obs_to_include
        if sum(isnan.(ts.df[ts.groupingrows[a], ts.value])) > 0
            @warn string("There are nan values for ", a)
        end
        vals = ts.df[ts.groupingrows[a], ts.value]
        cond2 = (maximum(vals) - minimum(vals) > 1000*eps())
        if (cond1 & cond2) push!(assets, a) end
    end
    return sort(assets)
end

"""
    subset_to_tick(ts::SortedDataFrame, n::Integer)

This subsets a `SortedDataFrame` to only the first n ticks.
### Inputs
* `ts` - Tick data.
* `n` - How many ticks to subset to.
### Returns
* A (smaller) `SortedDataFrame`.
"""
function subset_to_tick(ts::SortedDataFrame, n::Integer)
    newdf = ts.df[1:n,:]
    assets = get_assets(ts, 0)
    newgroupingrows = Dict{Symbol,Vector{eltype(ts.groupingrows[assets[1]])}}()
    for a in assets
        newvec = ts.groupingrows[a][ts.groupingrows[a] .<= n]
        if length(newvec) > 0 newgroupingrows[a] = newvec end
    end
    return SortedDataFrame(newdf, ts.time, ts.grouping, ts.value, newgroupingrows, ts.time_period_per_unit)
end

"""
    subset_to_time(ts::SortedDataFrame, totime::Real)

This subsets a `SortedDataFrame` to only the first observations up until some time.
### Inputs
* `ts` - Tick data.
* `totime` - Up to what time.
### Returns
* A (smaller) `SortedDataFrame`.
"""
function subset_to_time(ts::SortedDataFrame, totime::Real)
    ind = searchsortedfirst(ts.df[:,ts.time], totime) - 1
    if ind >= nrow(ts.df) return ts end
    return subset_to_tick(ts, ind)
end

"""
    duration(ts::SortedDataFrame)

The time elapsed between the first and the last tick in a `SortedDataFrame`.
### Inputs
* `ts` - Tick data.
* `in_dates_period` - In Dates.Period format or just a number for the numeric difference between first and last tick.
### Returns
* A scalar representing this duration.
"""
function duration(ts::SortedDataFrame; in_dates_period::Bool = true)
    duration_in_time_col_units = ts.df[nrow(ts.df),ts.time] - ts.df[1,ts.time]
    if !in_dates_period return duration_in_time_col_units end
    duration_in_nanoseconds = round(Nanosecond(ts.time_period_per_unit).value * duration_in_time_col_units)
    return Nanosecond(duration_in_nanoseconds)
end

"""
    CovarianceMatrix(correlation::Hermitian{R}, volatility::Vector{R},
                     labels::Vector{Symbol}) where R<:Real

This `Struct` stores three elements. A `Hermitian` correlation matrix, a vector of volatilities
and a vector of labels. The order of the labels matches the order of the assets in
the volatility vector and correlation matrix.
The default constructor is used.
### Inputs
* `correlation` - A `Hermitian` correlation matrix.
* `volatility` - Volatilities for each asset.
* `labels` - The labels for the `correlation` and `volatility` members. The n'th entry of the `labels` vector should contain the name of the asset that has its volatility in the n'th entry of the `volatility` member and its correlations in the n'th row/column of the `correlation` member.
* `time_period_per_unit` - The period that one unit of volatility corresponds to.
### Returns
* A `CovarianceMatrix`.
"""
mutable struct CovarianceMatrix{R<:Real}
    correlation::Hermitian{R}
    volatility::Vector{R}
    labels::Vector{Symbol}
    time_period_per_unit::Dates.Period
end

"""
    show(cm::CovarianceMatrix)
    show(cm::CovarianceMatrix, decimal_places_volatility::Integer, decimal_places_correlation::Integer)

This prints the `CovarianceMatrix` in a nice format.
### Inputs
* `cm` - The `CovarianceMatrix` you want to print.
* `decimal_places_volatility`  - The number of digits you want to show for volatilities. If not input then all digits are shown.
* `decimal_places_correlation` - The number of digits you want to show. If not input then all digits are shown.
"""
function show(cm::CovarianceMatrix)
    println()
    println("Volatilities per time interval of ", cm.time_period_per_unit)
    flat_labels = reshape( cm.labels , (1, length(cm.labels)))
    Base.print_matrix(stdout, vcat(flat_labels, cm.volatility'))
    println("\n")
    println("Correlations")
    corr = cm.correlation
    corr = vcat(flat_labels, corr)
    corr = hcat([:___, cm.labels...], corr)
    Base.print_matrix(stdout, corr)
    println("\n")
end
function show(cm::CovarianceMatrix, sig_figs_volatility::Integer, decimal_places_correlation::Integer)
    cm2 = deepcopy(cm)
    cm2.volatility  = round.(cm2.volatility , sigdigits = sig_figs_volatility)
    cm2.correlation = Hermitian(round.(cm2.correlation, digits = decimal_places_correlation))
    show(cm2)
end

"""
    convert_vol(vol::Union{Vector{<:Real},Real}, vol_period::Dates.Period, new_vol_period::Dates.Period)

### Inputs
* `Vol` - The volatility you want to convert.
* `vol_period` - The period of the volatility you have.
* `new_vol_period` - The period that you want to convert the volatility to.
"""
function convert_vol(vol::Union{Vector{<:Real},Real}, vol_period::Dates.Period, new_vol_period::Dates.Period)
    return vol .* sqrt( time_period_ratio(new_vol_period, vol_period) )
end
convert_vol(vol::Missing, vol_period::Dates.Period, new_vol_period::Dates.Period) = missing

"""
    print(cm::CovarianceMatrix)

This prints the `CovarianceMatrix` in a nice format.
### Inputs
* `cm` - The `CovarianceMatrix` you want to print.
* `decimal_places_volatility`  - The number of digits you want to show for volatilities. If not input then all digits are shown.
* `decimal_places_correlation` - The number of digits you want to show. If not input then all digits are shown.
"""
function print(cm::CovarianceMatrix)
    show(cm)
end
function print(cm::CovarianceMatrix, decimal_places_volatility::Integer, decimal_places_correlation::Integer)
    show(cm, decimal_places_volatility, decimal_places_correlation)
end

"""
    make_nan_covariance_matrix(labels::Vector{Symbol}, time_period_per_unit::Dates.Period)

This makes an empty `CovarianceMatrix` struct with all volatilities and correlations being NaNs.
### Inputs
* `labels` - The names of the asset names for this (empty) `CovarianceMatrix`.
* `time_period_per_unit` - The time interval the volatilities will be for.
### Returns
* An (empty) `CovarianceMatrix`
"""
function make_nan_covariance_matrix(labels::Vector{Symbol}, time_period_per_unit::Dates.Period)
    d = length(labels)
    correlation = ones(d,d)
    correlation .= NaN
    correlation[diagind(correlation)] .= 1
    vols = ones(d)
    vols .= NaN
    return CovarianceMatrix(Hermitian(correlation), vols, labels, time_period_per_unit)
end

"""
    calculate_mean_abs_distance(cov1::CovarianceMatrix, cov2::CovarianceMatrix, decimal_places::Integer = 8)

Calculates the mean absolute distance (elementwise in L1 norm) between two `CovarianceMatrix`s.
Undefined if any labels differ between the two `CovarianceMatrix`s.
### Inputs
* `cov1` - The first `CovarianceMatrix`
* `cov2` - The second `CovarianceMatrix`
* `decimal_places` - How many decimal places to show the result to.
### Returns
* An `Tuple` with the distance for correlations in first entry and distance for volatilities in the second.

    calculate_mean_abs_distance(d1::Dict{Symbol,<:Real}, d2::Dict{Symbol,<:Real})

Calculates the mean absolute distance (elementwise in L1 norm) between two `CovarianceMatrix`s.
### Inputs
* `d1` - The first `Dict`
* `d2` - The second `Dict`
### Returns
* A scalar with the mean distance between matching elements.
"""
function calculate_mean_abs_distance(cov1::CovarianceMatrix, cov2::CovarianceMatrix, decimal_places::Integer = 8)
    if length(symdiff(cov1.labels, cov2.labels)) != 0 return NaN, NaN end
    N = length(cov1.labels)
    cov11 = rearrange(cov1, cov2.labels)
    cor_error = round(sum(abs.(cov11.correlation .- cov2.correlation)) / ((N^2-N)/2), digits = decimal_places)
    vol_error = round(mean(abs.(cov11.volatility  .- cov2.volatility)), digits  = decimal_places)
    return (Correlation_error = cor_error, Volatility_error = vol_error)
end
function calculate_mean_abs_distance(d1::Dict{Symbol,<:Real}, d2::Dict{Symbol,<:Real})
    dist = 0.0
    common_labels = intersect(collect(keys(d1)), collect(keys(d2)))
    for label in common_labels
        dist += abs(d1[label] - d2[label])
    end
    return dist / length(common_labels)
end

"""
    get_correlation(covar::CovarianceMatrix, asset1::Symbol, asset2::Symbol)

Extract the correlation between two assets stored in a CovarianceMatrix.
### Inputs
* `covar` - A `CovarianceMatrix`
* `asset1` - A `Symbol` representing an asset.
* `asset2` - A `Symbol` representing an asset.
### Returns
* A Scalar (the correlation coefficient).
"""
function get_correlation(covar::CovarianceMatrix, asset1::Symbol, asset2::Symbol)
    index1 = findfirst(asset1 .== covar.labels)
    index2 = findfirst(asset2 .== covar.labels)
    if isnothing(index1) | isnothing(index2) return missing end
    return covar.correlation[index1, index2]
end

"""
    get_volatility(covar::CovarianceMatrix, asset1::Symbol, time_period_per_unit::Dates.Period = Second(1))

Get the volatility for a stock from a `CovarianceMatrix`.
### Inputs
* `covar` - A `CovarianceMatrix`
* `asset1` - A `Symbol` representing an asset.
* `time_period_per_unit` - The time interval the volatilities will be for.
### Returns
* A Scalar (the volatility).
"""
function get_volatility(covar::CovarianceMatrix, asset1::Symbol, time_period_per_unit::Dates.Period = covar.time_period_per_unit)
    index1 = findfirst(asset1 .== covar.labels)
    if isnothing(index1) return missing end
    return convert_vol(covar.volatility[index1], covar.time_period_per_unit, time_period_per_unit)
end

"""
    is_psd_matrix(mat::Hermitian, min_eigen_threshold::Real = 0.0)
    is_psd_matrix(covar::CovarianceMatrix, min_eigen_threshold::Real = 0.0)

Test if a matrix is psd (Positive Semi-Definite). This is done by seeing if all eigenvalues are positive.
If a `Hermitian` is input then it will be tested. If a `CovarianceMatrix` is input then its correlation matrix will be tested.
### Inputs
* `mat` - A `Hermitian` matrix or a `CovarianceMatrix`
* `min_eigen_threshold` - How big does the smallest eigenvalue have to be.
### Returns
* A `Bool` that is true if mat is psd and false if not.
"""
function is_psd_matrix(mat::Hermitian, min_eigen_threshold::Real = 0.0)
    eig = eigen(mat).values
    if length(eig) == 0 return false end # There is no eigenvalue decomposition.
    return minimum(eig) > 0 # is it PSD
end
is_psd_matrix(covar::CovarianceMatrix) = is_psd_matrix(covar.correlation)

"""
    valid_correlation_matrix(mat::Hermitian)
    valid_correlation_matrix(covar::CovarianceMatrix)

Test if a `Hermitian` matrix is a valid correlation matrix. This is done by testing if it is psd, if it has a unit diagonal and if all other elements are less than one.
If a `Hermitian` is input then it will be tested. If a `CovarianceMatrix` is input then its correlation matrix will be tested.
### Inputs
* `mat` - A `Hermitian` matrix or a `CovarianceMatrix`
* `min_eigen_threshold` - How big does the smallest eigenvalue have to be.
### Returns
* A `Bool` that is true if mat is a valid correlation matrix and false if not.
"""
function valid_correlation_matrix(mat::Hermitian, min_eigen_threshold::Real = 0.0)
    A = is_psd_matrix(mat, min_eigen_threshold)
    B = all(abs.(diag(mat) .- 1) .< 10*eps()) # does it have a unit diagonal
    C = all(abs.(mat) .<= 1 + 10*eps())       # all all off diagonals less than one in absolute value
    return all([A,B,C])
end
valid_correlation_matrix(covar::CovarianceMatrix, min_eigen_threshold::Real = 0.0) = valid_correlation_matrix(covar.correlation, min_eigen_threshold)

"""
    ticks_per_asset(ts::SortedDataFrame, assets::Vector{Symbol} = get_assets(ts))

Count the number of observations for each asset.
### Inputs
* `ts` - The tick data
* `assets` - A vector with asset `Symbol`s.
### Returns
* A `Dict` with the number of observations for each input asset.
"""
function ticks_per_asset(ts::SortedDataFrame, assets::Vector{Symbol} = get_assets(ts))
    ticks_per_asset = map(a -> length(ts.groupingrows[a]), assets)
    return Dict{Symbol,eltype(ticks_per_asset)}(assets .=> ticks_per_asset)
end


text:
### A Pluto.jl notebook ###
# v0.14.7

using Markdown
using InteractiveUtils

# ╔═╡ ab5e05cd-840c-4133-8808-b883c4f3a69e
begin
    import Pkg
    Pkg.activate(mktempdir())
    Pkg.add([
        Pkg.PackageSpec(name="PlutoUI", version="0.7"),
    ])
    using PlutoUI
end

# ╔═╡ ad530bb1-78bb-40a3-93c6-adec42207417
PlutoUI.TableOfContents(aside=true, title="Contenidos")

# ╔═╡ 74499a68-c58c-11eb-278e-478d05fe686e
md"""
# Scratch y sintaxis de Julia

**Un aspecto de los lenguajes basados en texto con el muchos estudiantes kidian es la comprención de la sintaxis especícifica (las reglas del lenguaje) requerida. Muchas veces los errores en un programa se deben a que no se siguene estas reglas, y estos erres son denominados errores de sintaxis o _sintax errors_.
Por lo tanto, es útil mostrar los paralelismos y las diferencias entre un idioma que un alumno ya ha dominado y el que está tratando de aprender.**

En los siguientes pasos, les mostraré algunos bloques de Scratch y su código equivalente en Julia. La lista está lejos de ser exhaustiva y está pensada como una guía de referencia más que como un ejercicio sobre el que trabajar.

"""

# ╔═╡ ee6d8355-542e-4e08-9279-4c05132d7302
md"""## Asignación de variables

- En scratch una variable necesita ser creada antes de que se le pueda asignar un valor, mientras que en Julia, una variable es creada cuando se le asigna un valor.
"""

# ╔═╡ 71935899-12a3-4ad4-9084-b8b8f7a3c720
md"""
![](assets/img/crearVariables.png)"""

# ╔═╡ 0af4e606-09c5-4a6c-b4be-9c202058e192
md"""
![](assets/img/puntosBienvendia.png)
"""

# ╔═╡ 7987cc38-f28d-4ee4-8a1c-738ecc5f81da
md"""
!!! hint "Si querés interactuar con esta página"
	Para usar la notebook hacé click en el boton arriba a la derecha
	![](https://alvarmaciel.github.io/paginas-pluto/assets/img/edit.png)

	y luego en el botón Binder

	![](https://alvarmaciel.github.io/paginas-pluto/assets/img/binder.png)
	>Toma un rato en cargar todo, pero con estos dos pasos se peude uasr la notebook
"""

# ╔═╡ 190473c5-abf5-488c-a81c-63b33d7c91e2
md"""
En Julia es necesario rodear los Strings o cualquier texto con comillas smples (') o dobles (")."""

# ╔═╡ c694ee4f-bcb0-4fd5-9fe7-c011523fd864
md"""
Corregí la variable `bienvenida` para que sea un string y poné un mensaje de bienvenida
"""

# ╔═╡ 0bbbe9f5-4870-4d11-886f-84aae2cc2cbd
bienvenida = 42

# ╔═╡ c34c3202-0d9e-4fca-bbdb-65217aa8bbef
if !(@isdefined bienvenida)
	md"""
!!! danger "Atencion"
	no cambies el nombre de la variable es `bienvenida`
	"""
end

# ╔═╡ bc4a39ea-74b0-45d1-92d2-de89bb8b1c3d
md"""
!!! hint "Ayuda"
    Los strings o cadenas de textos van entre comillas simples o dobles
"""

# ╔═╡ 4b014a11-c362-42dd-8c3e-7e30cb564fa2
if (@isdefined bienvenida) && (bienvenida isa String)
	md"""
!!! correct "Muy Bien"
    Bien hecho
"""

else
	md"""
!!! warning "Atencion"
Todavía falta. No cambiaste la variable
"""
end

# ╔═╡ 0a6df781-6bef-4955-bb8d-75a99fd5f09c
md"""
## Incrementar una variable

- En Scratch, se puede aumentar o reducir el valor de una variable. Para hacer esto, se una un bloque _cambiar_ con la variable y el número a incrementar, y usamos un número negativo para reducir el valor.
- Para aumentar o o reducir el valor de una variable en Julia, tienen que calcular el nuevo monto usando el valor guardado en la variable, y asignandole el nuevo vvalor a la variable.
![](assets/img/sumarA.png)
``` julia
puntos = puntos + 1
```
"""

# ╔═╡ 553e2e22-9b56-42ce-8f5a-4676c9faf024
md"""
## Salida simple

- En Scratch, hacemos que un objeto hable para generar una salida al usuario del programa.
- Julia usa las sentencias ```print``` y ```println``` para generar una salida de texto, números o simbolos.
- En Julia, nececitan usar comillas simples o dobles para imprimer texto.
- En las notebooks interactivas como esta (PlutoJl) la salida de las sentencias ```print``` y ```println``` se visualizan en el REPL. Para visaulizar una salida en PlutoJl con solo enunciar la variable, la celda mostrará su valor.

![](assets/img/puntajePuntos.png)

En julia se escribe:
``` julia
print("Tu puntaje es asombroso!!!")
print(puntos)
```
"""

# ╔═╡ fb03bd9c-ed6c-48a5-a3cd-08a4af5b405a
md"""
En Pluto simplemente enunciamos las variables, si son más de dos sentencias en un bloque. además podemos incluir texto en formato markdown como este que estás viendo.
```julia
begin

end
```
"""

# ╔═╡ 1e9e7329-bcf0-4b72-a4c6-d13276d9577f
begin
	puntos = 8
	md"""
	En esta celda pueden ver el código.
	- Todo lo que está entre ```md``` y las comillas es texto.
	- Para incorporar el valor de la variable al texto se usa el símbolo ```$```

	$bienvenida
	Tu puntaje es $puntos"""
end

# ╔═╡ 77c8a723-499e-478d-a519-49b7274ad179
md"""
## Estrucutras repetitivas condicionales

- La estructura de repetición condicional en Scratch es repetir **hasta** que cieta condición es Verdadera.
- La repetición condicional en Julia, repite **mientras**  cierta condición sea Verdadera.
- Pueden observar que el código dentro de la repetición está identado, La identación no es necesaria en Julia pero se recomienda. Se usan 4 espacios o un tab. Esto puede compararse a la forma en la que la los bloques de repeteción de Scratch encapsulan el código dentro de él.
![](assets/img/repetirHasta.png)
``` julia
while puntos <= 10
	print("jugando")
end
```
"""

# ╔═╡ a33d895b-34ba-468f-8366-40b069849041
md"""
## Repeticiones infinitas (en proceso)
"""

# ╔═╡ 0cba104c-627a-457c-8aeb-dfa938c5b467
md"""
## Estrucuturas Condicionales  (en proceso)
"""

# ╔═╡ 31a22a9a-5ac0-4304-8e2a-104cea9407f2
md"""
## Pruebas de igualdad  (en proceso)
"""

# ╔═╡ ca296a81-ce11-40d1-a5db-fd10329acfe1
md"""
# Desafío  (en proceso)
"""

# ╔═╡ f8e19d35-cb0c-4188-bd0f-ec298a493c1b
md"""
# Creditos y Licencia
>Estos textos/imagenes/videos fueron creados en base al material de Raspberry Pi Foundation y están bajo la licencia [Open Government Licence v3.0.](http://www.nationalarchives.gov.uk/doc/open-government-licence/version/3/).

>This text/image/video was created by the Raspberry Pi Foundation and is licensed under the [Open Government Licence v3.0.](http://www.nationalarchives.gov.uk/doc/open-government-licence/version/3/).
"""

# ╔═╡ cc1818e5-4608-4d06-aff9-76bf8dc41a54
md"""
!!! warning "Fuente:"
- [Scratch and Python syntax en ](https://www.futurelearn.com/courses/block-to-text-based-programming/13/steps/1129535) en [Scratch to Python: Moving from Block- to Text-based Programming](https://www.futurelearn.com/courses/block-to-text-based-programming)
"""

# ╔═╡ Cell order:
# ╟─ad530bb1-78bb-40a3-93c6-adec42207417
# ╟─ab5e05cd-840c-4133-8808-b883c4f3a69e
# ╟─74499a68-c58c-11eb-278e-478d05fe686e
# ╟─ee6d8355-542e-4e08-9279-4c05132d7302
# ╟─71935899-12a3-4ad4-9084-b8b8f7a3c720
# ╟─0af4e606-09c5-4a6c-b4be-9c202058e192
# ╟─7987cc38-f28d-4ee4-8a1c-738ecc5f81da
# ╟─190473c5-abf5-488c-a81c-63b33d7c91e2
# ╟─c694ee4f-bcb0-4fd5-9fe7-c011523fd864
# ╠═0bbbe9f5-4870-4d11-886f-84aae2cc2cbd
# ╟─c34c3202-0d9e-4fca-bbdb-65217aa8bbef
# ╟─bc4a39ea-74b0-45d1-92d2-de89bb8b1c3d
# ╟─4b014a11-c362-42dd-8c3e-7e30cb564fa2
# ╟─0a6df781-6bef-4955-bb8d-75a99fd5f09c
# ╟─553e2e22-9b56-42ce-8f5a-4676c9faf024
# ╟─fb03bd9c-ed6c-48a5-a3cd-08a4af5b405a
# ╠═1e9e7329-bcf0-4b72-a4c6-d13276d9577f
# ╟─77c8a723-499e-478d-a519-49b7274ad179
# ╟─a33d895b-34ba-468f-8366-40b069849041
# ╟─0cba104c-627a-457c-8aeb-dfa938c5b467
# ╟─31a22a9a-5ac0-4304-8e2a-104cea9407f2
# ╟─ca296a81-ce11-40d1-a5db-fd10329acfe1
# ╟─f8e19d35-cb0c-4188-bd0f-ec298a493c1b
# ╟─cc1818e5-4608-4d06-aff9-76bf8dc41a54

max_stars_repo_name: alvarmaciel/paginas-pluto

max_issues_repo_name: alvarmaciel/paginas-pluto

max_forks_repo_name: alvarmaciel/paginas-pluto


text:
n = 9
p = 4
nd= 9
x0 = [0,0,0,0,0,0,0,0,0]
d(i,j,N) = (j==1 ? 1*sin(2*pi/N*i) : 0) + (j==3 ? 2*sin(4*pi/N*i) : 0) + (j==5 ? 2*i/N : 0)
dt = .01
Q = [1,0,1,0,1,0,1,1,1]
Qf= [1,0,1,0,1,0,1,1,1]/dt
R = ones(4)/10

function nlm_ocp(;N=1,optimizer=Ipopt.Optimizer,opt...)
    m = SimpleNL.Model(optimizer;opt...)

    x=[variable(m,start = 0) for i=1:N+1,j=1:n]
    u=[variable(m,start = 0) for i=1:N,j=1:p]
    [constraint(m,x[1,i]-x0[i]) for i=1:n]
    for i=1:N
        constraint(m, -x[i+1,1] + x[i,1] + (x[i,2])*dt)
        constraint(m, -x[i+1,2] + x[i,2] + (u[i,1]*cos(x[i,7])*sin(x[i,8])*cos(x[i,9])+u[i,1]*sin(x[i,7])*sin(x[i,9]))*dt)
        constraint(m, -x[i+1,3] + x[i,3] + (x[i,4])*dt)
        constraint(m, -x[i+1,4] + x[i,4] + (u[i,1]*cos(x[i,7])*sin(x[i,8])*sin(x[i,9])-u[i,1]*sin(x[i,7])*cos(x[i,9]))*dt)
        constraint(m, -x[i+1,5] + x[i,5] + (x[i,6])*dt)
        constraint(m, -x[i+1,6] + x[i,6] + (u[i,1]*cos(x[i,7])*cos(x[i,8])-9.8)*dt)
        constraint(m, -x[i+1,7] + x[i,7] + (u[i,2]*cos(x[i,7])/cos(x[i,8])+u[i,3]*sin(x[i,7])/cos(x[i,8]))*dt)
        constraint(m, -x[i+1,8] + x[i,8] + (-u[i,2]*sin(x[i,7])+u[i,3]*cos(x[i,7]))*dt)
        constraint(m, -x[i+1,9] + x[i,9] + (u[i,2]*cos(x[i,7])*tan(x[i,8])+u[i,3]*sin(x[i,7])*tan(x[i,8])+u[i,4])*dt)
    end
    objective(m,.5*sum(Q[j]*(x[i,j]-d(i,j,N))^2 for i=1:N for j=1:n) + .5*sum(R[j]*(u[i,j]^2) for i=1:N for j=1:p)
              + .5*sum(Qf[j]*(x[N+1,j]-d(N+1,j,N))^2 for j=1:n))

    return m
end

function jump_ocp(;N=1,optimizer=Ipopt.Optimizer)
    m = JuMP.Model(optimizer)
    @variable(m,x[1:N+1,1:n],start = 0)
    @variable(m,u[1:N,1:p],start = 0)
    @constraint(m,[i=1:n],x[1,i]==x0[i])
    @constraint(m,[i=1:N],x[i+1,1] == x[i,1] + (x[i,2])*dt)
    @NLconstraint(m,[i=1:N], x[i+1,2] == x[i,2] + (u[i,1]*cos(x[i,7])*sin(x[i,8])*cos(x[i,9])+u[i,1]*sin(x[i,7])*sin(x[i,9]))*dt)
    @constraint(m,[i=1:N], x[i+1,3] == x[i,3] + (x[i,4])*dt)
    @NLconstraint(m,[i=1:N], x[i+1,4] == x[i,4] + (u[i,1]*cos(x[i,7])*sin(x[i,8])*sin(x[i,9])-u[i,1]*sin(x[i,7])*cos(x[i,9]))*dt)
    @constraint(m,[i=1:N], x[i+1,5] == x[i,5] + (x[i,6])*dt)
    @NLconstraint(m,[i=1:N], x[i+1,6] == x[i,6] + (u[i,1]*cos(x[i,7])*cos(x[i,8])-9.8)*dt)
    @NLconstraint(m,[i=1:N], x[i+1,7] == x[i,7] + (u[i,2]*cos(x[i,7])/cos(x[i,8])+u[i,3]*sin(x[i,7])/cos(x[i,8]))*dt)
    @NLconstraint(m,[i=1:N], x[i+1,8] == x[i,8] + (-u[i,2]*sin(x[i,7])+u[i,3]*cos(x[i,7]))*dt)
    @NLconstraint(m,[i=1:N], x[i+1,9] == x[i,9] + (u[i,2]*cos(x[i,7])*tan(x[i,8])+u[i,3]*sin(x[i,7])*tan(x[i,8])+u[i,4])*dt)
    @objective(m,Min, .5*sum(Q[j]*(x[i,j]-d(i,j,N))^2 for i=1:N for j=1:n) + .5*sum(R[j]*(u[i,j]^2) for i=1:N for j=1:p)
                 + .5*sum(Qf[j]*(x[N+1,j]-d(N+1,j,N))^2 for j=1:n))
    return m
end

function nlm_luksan_vlcek_501(;N=1,optimizer=Ipopt.Optimizer,opt...)
    m = SimpleNL.Model(optimizer;opt...)

    x = [SimpleNL.variable(m;start=mod(i,2)==1 ? -1.2 : 1.) for i=1:N];
    for i=1:N-2
        SimpleNL.constraint(m,3x[i+1]^3+2*x[i+2]-5+sin(x[i+1]-x[i+2])sin(x[i+1]+x[i+2])+4x[i+1]-x[i]exp(x[i]-x[i+1])-3)
    end
    SimpleNL.objective(m,sum(100(x[i-1]^2-x[i])^2+(x[i-1]-1)^2 for i=2:N))

    return m
end

function jump_luksan_vlcek_501(;N=1,optimizer=Ipopt.Optimizer)
    m=JuMP.Model(optimizer)
    @variable(m,x[i=1:N], start= mod(i,2)==1 ? -1.2 : 1.)
    @NLconstraint(m,[i=1:N-2], 3x[i+1]^3+2x[i+2]-5+sin(x[i+1]-x[i+2])sin(x[i+1]+x[i+2])+4x[i+1]-x[i]exp(x[i]-x[i+1])-3==0.)
    @NLobjective(m,Min,sum(100(x[i-1]^2-x[i])^2+(x[i-1]-1)^2 for i=2:N))
    return m
end

function parse_ipopt_output(filename)
    output = read(filename,String);
    i = match(r"Total CPU secs in IPOPT",output).offset
    j = match(r"Total CPU secs in NLP function evaluations",output).offset;
    k = match(r"EXIT:",output).offset

    cpu_time_ipopt = parse(Float64,output[i+54:j-3])
    cpu_time_nlp_function = parse(Float64,output[j+54:k-3])
    cpu_time_total = cpu_time_ipopt + cpu_time_nlp_function

    return cpu_time_total, cpu_time_nlp_function
end


text:
using TSXPlayground.TSX

function f()
    x = 0
    if (status = TSX.xbegin()) == TSX.XBEGIN_STARTED
        x = 1
        TSX.xend()
    end
    return (x, status)
end


text:
function mockline(; a=a, b=b, t_range=t_range, σ=0.1)

    observations = zeros(3, length(t_range))

    for (index, t) in enumerate(t_range)

        for i in 1:3

            observations[i, index] = a[i]*t + b[i] + randn()*σ

        end

    end

    @printf("Intersection point is at\n") ; display(b) ; @printf("\n")

    figure()

    plot3D(observations[1, :], observations[2, :], observations[3, :], "bo", label="observations")

    plotvector(b, "k-", "vector b")

    legend()

    return observations, σ*ones(size(observations))

end


text:
@testset "Basic" begin
    @testset "Basic: prefix_sum" begin
        @test prefix_sum([1, 1, 1]) == [1, 2, 3]
        @test prefix_sum([1, 2, 3]) == [1, 3, 6]
    end

    @testset "Basic: DifferenceArray" begin
        arr = [1, 3, 6, 10, 15, 21]
        diff_arr = DifferenceArray.create_diff_arr(arr)
        @test diff_arr == [1, 2, 3, 4, 5, 6]

        #Increase elements with index from [3, 5] by 1
        DifferenceArray.add_to_arr(diff_arr, 3, 5, 1)
        @test DifferenceArray.calculate_arr(diff_arr) == [1, 3, 7, 11, 16, 21]

        #Increase all elements by 1
        DifferenceArray.add_to_arr(diff_arr, 1, 6, 1)
        @test DifferenceArray.calculate_arr(diff_arr) == [2, 4, 8, 12, 17, 22]
    end

    @testset "Basic: The Tower of Hanoi" begin
        # a, b, c represents the names of the first, second and third rod respectively
        @test Hanoi.solve('a', 'b', 'c', 1) == ['a' => 'c'] # One ring, just move it from the first to the third rod
        @test Hanoi.solve('a', 'b', 'c', 2) == ['a' => 'b', 'a' => 'c', 'b' => 'c']
        @test Hanoi.solve('a', 'b', 'c', 3) == ['a' => 'c', 'a' => 'b', 'c' => 'b', 'a' => 'c', 'b' => 'a', 'b' => 'c', 'a' => 'c']
    end
end

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia


text:
using LeetCode
using Test

@testset "LeetCode.jl" begin
    include("common.jl")
    include("problems/problems.jl")
end


text:
# Autogenerated wrapper script for libavtp_jll for armv7l-linux-gnueabihf
export libavtp

JLLWrappers.@generate_wrapper_header("libavtp")
JLLWrappers.@declare_library_product(libavtp, "libavtp.so.0")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libavtp,
        "lib/libavtp.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# this script generates the moveout plot for all cross-corrlations
# in the benchmarking test

using SeisNoise, Glob, Plots, Statistics, Plots.PlotMeasures
# plot all correlations
CORRDIR = expanduser("~/SeisNoise_paper/SINGLECORR")
files = glob("*",CORRDIR)
freqmin = 1 / 50
freqmax = 1 / 30
max_length = 1000.

# load first file and allocate
C = load_corr(files[10000],"ZZ")
t = -C.maxlag:C.maxlag
corr = Array{eltype(C.corr)}(undef,size(C.corr,1),length(files))
dist = zeros(Float64,length(files))
azi = zeros(Float64,length(files))
baz = zeros(Float64,length(files))
lat = zeros(Float64,length(files))
lon = zeros(Float64,length(files))

# loop through and plot
for ii = 1:length(files)
    C = load_corr(files[ii],"ZZ")
    clean_up!(C,freqmin,freqmax,max_length=max_length)
    println(ii," ",C.name)
    corr[:,ii] .= C.corr[:,1]
    dist[ii] = C.dist
    azi[ii] = C.azi
    baz[ii] = C.baz
    lon[ii] = C.loc.lon
    lat[ii] = C.loc.lat
end

# sort by cross-correlation by distance
distind = sortperm(dist)
corr .= corr[:,distind]
dist .= dist[distind]
azi .= azi[distind]
baz .= baz[distind]
lon .= lon[distind]
lat .= lat[distind]

### plotting ###
# just find the highest snr traces
maxamp = maximum(abs.(corr),dims=1)
maxmad = mad(corr)
snr = (maxamp ./ maxmad)[:]
nonan = findall(.!isnan.(snr[:]))
snr = snr[nonan]
corr = corr[:,nonan]
dist = dist[nonan]
azi = azi[nonan]
baz = baz[nonan]
lon = lon[nonan]
lat = lat[nonan]

scaling = 10000 # scale waveforms by 1000
lw = 0.5  # linewidth
al = 0.1  # alpha
mindist = 100. # minimum distance in km for plot
maxdist = 10000.
distave = 10.
dists = collect(mindist:distave:maxdist)
maxt = 3000
tind = findall(abs.(t) .< maxt)
distplot = zeros(eltype(C.corr),length(tind),length(dists))
for ii = 1:length(dists)
    cind = findall((dist .> dists[ii]) .& (dist .< dists[ii] + distave))
    c = mean(corr[tind,cind] ./ (maximum(abs.(corr[tind,cind]),dims=1) .+ 1),dims=2)
    distplot[:,ii] = c ./ maximum(abs.(c))
end
heatmap(dists,t[tind],distplot,color=:RdBu, xlabel="Interstation Distance [km]",ylabel="Lag [s]")

p = plot()
for ii = 1:length(dists)
    cind = findall((dist .> dists[ii]) .& (dist .< dists[ii] + distave))
    c = mean(corr[tind,cind] ./ (maximum(abs.(corr[tind,cind]),dims=1) .+ 1),dims=2)
    # c = corr[tind,snrind[ii]] ./ maximum(abs.(corr[tind,snrind[ii]]),dims=1)
    plot!(p,t[tind],dists[ii] .+ c .* scaling,color=:black,linewidth=lw,alpha=al,
         label="",xlabel = "Lag [s]",ylabel="Interstation Distance [km]",
         xlims=(-maxt,maxt),ylims=(mindist,maxdist),dpi=500,left_margin=-60px,
                bottom_margin=-20px,right_margin=25px)
end
p
png(expanduser("~/SeisNoise_paper/LHZ-moveout"))


text:
# Adapted from mnist.jl

using AutoGrad
using ArrayFire
using GZip
using Compat

function timeit(f=train0, w = weights(64; seed=1); epochs=10)
    isdefined(:dtrn) || loaddata()
    for i=1:3
        sleep(2)
        #gc_enable(false)
        @time f(w; epochs=epochs)
        #gc_enable(true)
    end
end

function profit(f=train0, w = weights(64; seed=1); epochs=10)
    timeit(f, w; epochs=epochs)
    sleep(2)
    #gc_enable(false)
    @profile f(w; epochs=epochs)
    #gc_enable(true)
end

function predict(w, x)
    i = 1
    while i+2 < length(w)
        x = max(0, w[i]*x .+ w[i+1])
        i += 2
    end
    return w[i]*x .+ w[i+1]
end

function loss(w,x,ygold)
    ypred = predict(w, x)
    sum((ypred-ygold).^2)
end

function loss1(w,x,ygold)
    ypred = predict(w, x)
    ynorm = ypred .- log(sum(exp(ypred),1))
    -sum(ygold .* ynorm) / size(ygold,2)
end

function train0(w=weights(64); lr=.1, epochs=1)
    gradfun = grad(loss)
    for epoch=1:epochs
        for (x,y) in dtrn
            g = gradfun(w, x, y)
            for i in 1:length(w)
                w[i] -= lr*g[i]
            end
        end
    end
    return w
end

function train1(w=weights(64); lr=.1, epochs=1)
    gradfun = grad(loss)
    for epoch=1:epochs
        for (x,y) in dtrn
            g = gradfun(w, x, y)
        end
    end
end

function train2(w=weights(64); lr=.1, epochs=1)
    for epoch=1:epochs
        for (x,y) in dtrn
            z = loss(w, x, y)
        end
    end
end

function train3(w=weights(64); lr=.1, epochs=1)
    for epoch=1:epochs
        for (x,ygold) = dtrn
            y = predict(w,x)
        end
    end
end

function weights(h...; seed=nothing)
    seed==nothing || srand(seed)
    w = Any[]
    x = 28*28
    for y in [h..., 10]
        push!(w, convert(Array{Float32}, 0.1*randn(y,x)))
        push!(w, zeros(Float32,y))
        x = y
    end
    return map(AFArray,w)
end

function loaddata()
    info("Loading data...")
    global xtrn, xtst, ytrn, ytst, dtrn
    xshape(a)=reshape(a./255f0,784,div(length(a),784))
    yshape(a)=(a[a.==0]=10; full(sparse(convert(Vector{Int},a),1:length(a),1f0)))
    xtrn = xshape(gzload("train-images-idx3-ubyte.gz")[17:end])
    xtst = xshape(gzload("t10k-images-idx3-ubyte.gz")[17:end])
    ytrn = yshape(gzload("train-labels-idx1-ubyte.gz")[9:end])
    ytst = yshape(gzload("t10k-labels-idx1-ubyte.gz")[9:end])
    dtrn = map(a->map(AFArray,a), minibatch(xtrn, ytrn, 100))
    # dtrn = [(CudaArray(xtst),CudaArray(ytst))]
    info("Loading done...")
end

function gzload(file; path=Knet.dir("data",file), url="http://yann.lecun.com/exdb/mnist/$file")
    isfile(path) || download(url, path)
    f = gzopen(path)
    a = @compat read(f)
    close(f)
    return(a)
end

function minibatch(x, y, batchsize)
    data = Any[]
    nx = size(x,2)
    for i=1:batchsize:nx
        j=min(i+batchsize-1,nx)
        push!(data, (x[:,i:j], y[:,i:j]))
    end
    return data
end

isdefined(:dtrn) || loaddata()

using AutoGrad: id2
import AutoGrad: _dbg, sum_helper
_dbg(x::AFArray)=Symbol("AF$(join([id2(x),size(x)...,8*sizeof(eltype(x))],'_'))")
sum_helper{T}(a::AFArray{T},b::AFArray{T},c::AFArray{T}...) = +(a,b,c...)

w1 = weights();
w2 = weights(64);
(x,y) = first(dtrn);
g = grad(loss)

# timeit(train0,w1); # => 6.713386 seconds (7.34 M allocations: 276.245 MB, 4.68% gc time) # compare to 1.977053 seconds (2.91 M allocations: 173.705 MB)  # compare to 1.8276 in Knet7
# timeit(train1,w1); # => 5.613988 seconds (6.75 M allocations: 258.575 MB, 4.26% gc time) # compare to 1.898212 seconds (2.88 M allocations: 173.156 MB)
# timeit(train2,w1); # => 1.213765 seconds (668.72 k allocations: 25.752 MB) # compare to 1.025943 seconds (631.80 k allocations: 74.277 MB)
# timeit(train3,w1); # => 0.464109 seconds (259.80 k allocations: 11.838 MB) # compare to 0.608363 seconds (601.80 k allocations: 27.127 MB)

# timeit(train0,w2); # => 11.899407 seconds (11.37 M allocations: 421.264 MB, 4.53% gc time) # compare to 3.659655 seconds (6.97 M allocations: 344.909 MB)  # compare to 2.9201 in Knet7
# timeit(train1,w2); # => 9.333509 seconds (10.21 M allocations: 386.199 MB, 4.29% gc time) # compare to 3.505348 seconds (6.90 M allocations: 343.810 MB)
# timeit(train2,w2); # => 1.874707 seconds (907.80 k allocations: 33.353 MB) # compare to 1.391519 seconds (1.36 M allocations: 107.053 MB)
# timeit(train3,w2); # => 0.856982 seconds (511.80 k allocations: 20.078 MB) # compare to 0.985798 seconds (1.33 M allocations: 59.903 MB)

# Plain forward is faster.  relu and softmax are much slower.  Also probably too much alloc.

# Replacing softloss with quadloss:
# timeit(train0,w1); # => 4.318921 seconds (4.92 M allocations: 181.396 MB, 5.20% gc time)
# timeit(train1,w1); # => 3.246138 seconds (4.33 M allocations: 163.726 MB, 5.24% gc time)
# timeit(train2,w1); # => 0.847536 seconds (373.80 k allocations: 12.113 MB)
# timeit(train3,w1); # => 0.464719 seconds (109.80 k allocations: 3.690 MB)

# timeit(train0,w2); # => 8.889076 seconds (8.94 M allocations: 324.584 MB, 4.87% gc time)
# timeit(train1,w2); # => 6.513208 seconds (7.78 M allocations: 289.520 MB, 4.31% gc time)
# timeit(train2,w2); # => 1.375004 seconds (625.80 k allocations: 20.352 MB, 3.89% gc time)
# timeit(train3,w2); # => 0.749930 seconds (361.80 k allocations: 11.929 MB)

# There is still a big slowdown in the backward_pass.  What operation is responsible?
# The number of allocations is smaller in forward, larger in backward compared to old impl.
# Is it related to not being able to fuse loops?
# Is there some ungetindex applied to a AFArray?


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: detective
using AWS.Compat
using AWS.UUIDs

"""
    AcceptInvitation()

Accepts an invitation for the member account to contribute data to a behavior graph. This operation can only be called by an invited member account.  The request provides the ARN of behavior graph. The member account status in the graph must be INVITED.

# Required Parameters
- `GraphArn`: The ARN of the behavior graph that the member account is accepting the invitation for. The member account status in the behavior graph must be INVITED.

"""
accept_invitation(GraphArn; aws_config::AWSConfig=global_aws_config()) = detective("PUT", "/invitation", Dict{String, Any}("GraphArn"=>GraphArn); aws_config=aws_config)
accept_invitation(GraphArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = detective("PUT", "/invitation", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("GraphArn"=>GraphArn), args)); aws_config=aws_config)

"""
    CreateGraph()

Creates a new behavior graph for the calling account, and sets that account as the master account. This operation is called by the account that is enabling Detective. Before you try to enable Detective, make sure that your account has been enrolled in Amazon GuardDuty for at least 48 hours. If you do not meet this requirement, you cannot enable Detective. If you do meet the GuardDuty prerequisite, then when you make the request to enable Detective, it checks whether your data volume is within the Detective quota. If it exceeds the quota, then you cannot enable Detective.  The operation also enables Detective for the calling account in the currently selected Region. It returns the ARN of the new behavior graph.  CreateGraph triggers a process to create the corresponding data tables for the new behavior graph. An account can only be the master account for one behavior graph within a Region. If the same account calls CreateGraph with the same master account, it always returns the same behavior graph ARN. It does not create a new behavior graph.

"""
create_graph(; aws_config::AWSConfig=global_aws_config()) = detective("POST", "/graph"; aws_config=aws_config)
create_graph(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = detective("POST", "/graph", args; aws_config=aws_config)

"""
    CreateMembers()

Sends a request to invite the specified AWS accounts to be member accounts in the behavior graph. This operation can only be called by the master account for a behavior graph.   CreateMembers verifies the accounts and then sends invitations to the verified accounts. The request provides the behavior graph ARN and the list of accounts to invite. The response separates the requested accounts into two lists:   The accounts that CreateMembers was able to start the verification for. This list includes member accounts that are being verified, that have passed verification and are being sent an invitation, and that have failed verification.   The accounts that CreateMembers was unable to process. This list includes accounts that were already invited to be member accounts in the behavior graph.

# Required Parameters
- `Accounts`: The list of AWS accounts to invite to become member accounts in the behavior graph. For each invited account, the account list contains the account identifier and the AWS account root user email address.
- `GraphArn`: The ARN of the behavior graph to invite the member accounts to contribute their data to.

# Optional Parameters
- `Message`: Customized message text to include in the invitation email message to the invited member accounts.
"""
create_members(Accounts, GraphArn; aws_config::AWSConfig=global_aws_config()) = detective("POST", "/graph/members", Dict{String, Any}("Accounts"=>Accounts, "GraphArn"=>GraphArn); aws_config=aws_config)
create_members(Accounts, GraphArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = detective("POST", "/graph/members", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Accounts"=>Accounts, "GraphArn"=>GraphArn), args)); aws_config=aws_config)

"""
    DeleteGraph()

Disables the specified behavior graph and queues it to be deleted. This operation removes the graph from each member account's list of behavior graphs.  DeleteGraph can only be called by the master account for a behavior graph.

# Required Parameters
- `GraphArn`: The ARN of the behavior graph to disable.

"""
delete_graph(GraphArn; aws_config::AWSConfig=global_aws_config()) = detective("POST", "/graph/removal", Dict{String, Any}("GraphArn"=>GraphArn); aws_config=aws_config)
delete_graph(GraphArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = detective("POST", "/graph/removal", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("GraphArn"=>GraphArn), args)); aws_config=aws_config)

"""
    DeleteMembers()

Deletes one or more member accounts from the master account behavior graph. This operation can only be called by a Detective master account. That account cannot use DeleteMembers to delete their own account from the behavior graph. To disable a behavior graph, the master account uses the DeleteGraph API method.

# Required Parameters
- `AccountIds`: The list of AWS account identifiers for the member accounts to delete from the behavior graph.
- `GraphArn`: The ARN of the behavior graph to delete members from.

"""
delete_members(AccountIds, GraphArn; aws_config::AWSConfig=global_aws_config()) = detective("POST", "/graph/members/removal", Dict{String, Any}("AccountIds"=>AccountIds, "GraphArn"=>GraphArn); aws_config=aws_config)
delete_members(AccountIds, GraphArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = detective("POST", "/graph/members/removal", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AccountIds"=>AccountIds, "GraphArn"=>GraphArn), args)); aws_config=aws_config)

"""
    DisassociateMembership()

Removes the member account from the specified behavior graph. This operation can only be called by a member account that has the ENABLED status.

# Required Parameters
- `GraphArn`: The ARN of the behavior graph to remove the member account from. The member account's member status in the behavior graph must be ENABLED.

"""
disassociate_membership(GraphArn; aws_config::AWSConfig=global_aws_config()) = detective("POST", "/membership/removal", Dict{String, Any}("GraphArn"=>GraphArn); aws_config=aws_config)
disassociate_membership(GraphArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = detective("POST", "/membership/removal", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("GraphArn"=>GraphArn), args)); aws_config=aws_config)

"""
    GetMembers()

Returns the membership details for specified member accounts for a behavior graph.

# Required Parameters
- `AccountIds`: The list of AWS account identifiers for the member account for which to return member details. You cannot use GetMembers to retrieve information about member accounts that were removed from the behavior graph.
- `GraphArn`: The ARN of the behavior graph for which to request the member details.

"""
get_members(AccountIds, GraphArn; aws_config::AWSConfig=global_aws_config()) = detective("POST", "/graph/members/get", Dict{String, Any}("AccountIds"=>AccountIds, "GraphArn"=>GraphArn); aws_config=aws_config)
get_members(AccountIds, GraphArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = detective("POST", "/graph/members/get", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AccountIds"=>AccountIds, "GraphArn"=>GraphArn), args)); aws_config=aws_config)

"""
    ListGraphs()

Returns the list of behavior graphs that the calling account is a master of. This operation can only be called by a master account. Because an account can currently only be the master of one behavior graph within a Region, the results always contain a single graph.

# Optional Parameters
- `MaxResults`: The maximum number of graphs to return at a time. The total must be less than the overall limit on the number of results to return, which is currently 200.
- `NextToken`: For requests to get the next page of results, the pagination token that was returned with the previous set of results. The initial request does not include a pagination token.
"""
list_graphs(; aws_config::AWSConfig=global_aws_config()) = detective("POST", "/graphs/list"; aws_config=aws_config)
list_graphs(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = detective("POST", "/graphs/list", args; aws_config=aws_config)

"""
    ListInvitations()

Retrieves the list of open and accepted behavior graph invitations for the member account. This operation can only be called by a member account. Open invitations are invitations that the member account has not responded to. The results do not include behavior graphs for which the member account declined the invitation. The results also do not include behavior graphs that the member account resigned from or was removed from.

# Optional Parameters
- `MaxResults`: The maximum number of behavior graph invitations to return in the response. The total must be less than the overall limit on the number of results to return, which is currently 200.
- `NextToken`: For requests to retrieve the next page of results, the pagination token that was returned with the previous page of results. The initial request does not include a pagination token.
"""
list_invitations(; aws_config::AWSConfig=global_aws_config()) = detective("POST", "/invitations/list"; aws_config=aws_config)
list_invitations(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = detective("POST", "/invitations/list", args; aws_config=aws_config)

"""
    ListMembers()

Retrieves the list of member accounts for a behavior graph. Does not return member accounts that were removed from the behavior graph.

# Required Parameters
- `GraphArn`: The ARN of the behavior graph for which to retrieve the list of member accounts.

# Optional Parameters
- `MaxResults`: The maximum number of member accounts to include in the response. The total must be less than the overall limit on the number of results to return, which is currently 200.
- `NextToken`: For requests to retrieve the next page of member account results, the pagination token that was returned with the previous page of results. The initial request does not include a pagination token.
"""
list_members(GraphArn; aws_config::AWSConfig=global_aws_config()) = detective("POST", "/graph/members/list", Dict{String, Any}("GraphArn"=>GraphArn); aws_config=aws_config)
list_members(GraphArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = detective("POST", "/graph/members/list", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("GraphArn"=>GraphArn), args)); aws_config=aws_config)

"""
    RejectInvitation()

Rejects an invitation to contribute the account data to a behavior graph. This operation must be called by a member account that has the INVITED status.

# Required Parameters
- `GraphArn`: The ARN of the behavior graph to reject the invitation to. The member account's current member status in the behavior graph must be INVITED.

"""
reject_invitation(GraphArn; aws_config::AWSConfig=global_aws_config()) = detective("POST", "/invitation/removal", Dict{String, Any}("GraphArn"=>GraphArn); aws_config=aws_config)
reject_invitation(GraphArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = detective("POST", "/invitation/removal", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("GraphArn"=>GraphArn), args)); aws_config=aws_config)

"""
    StartMonitoringMember()

Sends a request to enable data ingest for a member account that has a status of ACCEPTED_BUT_DISABLED. For valid member accounts, the status is updated as follows.   If Detective enabled the member account, then the new status is ENABLED.   If Detective cannot enable the member account, the status remains ACCEPTED_BUT_DISABLED.

# Required Parameters
- `AccountId`: The account ID of the member account to try to enable. The account must be an invited member account with a status of ACCEPTED_BUT_DISABLED.
- `GraphArn`: The ARN of the behavior graph.

"""
start_monitoring_member(AccountId, GraphArn; aws_config::AWSConfig=global_aws_config()) = detective("POST", "/graph/member/monitoringstate", Dict{String, Any}("AccountId"=>AccountId, "GraphArn"=>GraphArn); aws_config=aws_config)
start_monitoring_member(AccountId, GraphArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = detective("POST", "/graph/member/monitoringstate", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("AccountId"=>AccountId, "GraphArn"=>GraphArn), args)); aws_config=aws_config)


text:
## Exercise 8-6
## Run the program on paper, changing the values of i and j during each iteration. Find and fix the second error in this function.
println("Ans: ")

function isreverse(word1, word2)
    if length(word1) != length(word2)
        return false
    end

    i = firstindex(word1)
    j = lastindex(word2)

    while j > 0
        if word1[i] != word2[j]
            return false
        end
        j = prevind(word2, j)
        i = nextind(word1, i)
    end
    true
end

println(isreverse("haha", "ahah"))
println(isreverse("I cant think of any other string", "maybe"))
println(isreverse("manye", "maybe"))

println("End.")


text:
using IndexableBitVectors

srand(12345)

const step = 16

function bench_getindex(bv, n)
    range = 1:step:endof(bv)
    # warming up
    for i in range
        bv[i]
    end
    t = @elapsed for _ in 1:n
        for i in range
            bv[i]
        end
    end
    return t / (n * length(range))
end

function bench_rank1(bv, n)
    range = 1:step:endof(bv)
    for i in range
        rank1(bv, i)
    end
    t = @elapsed for _ in 1:n
        for i in range
            rank1(bv, i)
        end
    end
    return t / (n * length(range))
end

function bench_select1(bv, n)
    range = 1:step:endof(bv)
    for i in range
        select1(bv, i)
    end
    t = @elapsed for _ in 1:n
        for i in range
            select1(bv, i)
        end
    end
    return t / (n * length(range))
end

let
    name = shift!(ARGS)
    T = name == "CompactBitVector" ? CompactBitVector :
        name == "SucVector" ? SucVector :
        name == "RRR" ? RRR :
        name == "LargeRRR" ? LargeRRR : error("unknown type name: ", name)
    columns = ["type", "length", "r", "bench_getindex", "bench_rank1", "bench_select1"]
    println(join(columns, '\t'))
    for r in [0.01, 0.1, 0.5], p in 10:2:30
        #info("$r $p")
        len = 2^p
        vec = rand(Float32, len) .< r
        bv = convert(T, vec)
        print(T, '\t', len, '\t', r, '\t')
        n = max(1, div(2^26, 2^p))
        gc()
        print(bench_getindex(bv, n), '\t')
        gc()
        print(bench_rank1(bv, n), '\t')
        gc()
        print(bench_select1(bv, n), '\n')
    end
end