text:
# Shallow water equations using a Lax-Friedrich finite volume scheme.
# flat bottom. Equations :
# h_t + div(hU) = 0                        (1)
# (hU)_t + div(h UxU + 1/2 gh^2 Id_2) = 0  (2)
#
# Setting V = (h, hu, hv), where U = (u, v), equations write :
# V_t + F(V)_x + G(V)_y = 0
#
# Where :
# F(V) = (hu, hu^2 + 1/2 gh^2, huv)
# G(V) = (hv, huv, hv^2 + 1/2 gh^2)
#
# Lax-Friedrich for 1D equation v_t + F(v)_x = 0
# v^{n+1}_i = 1/2 (v^n_{i-1} + v^n_{i+1}) - dt / (2 dx) * (F(v^n_{i+1}) - F(v^n_{i-1}))
#
# For our problem, in 2D, with a cartesian mesh :
# V^{n+1}_{i, j} = 1/2 (V^n_{i-1, j} + V^n_{i+1, j} + V^n_{i, j-1} + V^n_{i, j+1})
#                  - dt / (2 dx) * (F(V^n_{i+1, j}) - F(V^n_{i-1, j}))
#                  - dt / (2 dy) * (G(V^n_{i, j+1}) - G(V^n_{i, j-1}))

# perf or animation ?
const INFO = false
const record_animation = false

using Printf
using StaticArrays
if record_animation
    using Makie
end

using Profile

# the type for the domain
struct Domain
    x_start::Float64
    x_end::Float64
    y_start::Float64
    y_end::Float64
end

# Initial solution
function InitSolution!(V::Array{Float64, 3}, domain::Domain)
    ny = size(V, 2)
    nx = size(V, 3)

    hx = (domain.x_end - domain.x_start) / nx
    hy = (domain.y_end - domain.y_start) / ny

    xc = 0.5  * (domain.x_start + domain.x_end) + nx/8 * hx
    yc = 0.5  * (domain.y_start + domain.y_end) - ny/8 * hy

    for i=1:nx
        x = domain.x_start + (i-0.5) * hx

        for j=1:ny
            y = domain.y_start + (j-0.5) * hy

            if (x - xc)^2 + (y - yc)^2 < 0.2^2
                V[1, j, i] = 2
            else
                V[1, j, i] = 1
            end

            V[2, j, i] = 0
            V[3, j, i] = 0
        end
    end
end

# Model
function F!(f, V::Array{Float64, 3}, j::Int64, i::Int64)
    f[1] = V[2, j, i]
    if isapprox(V[1, j, i], 0., atol=eps(Float64), rtol=0)
        f[2] = 0
        f[3] = 0
    else
        f[2] = V[2, j, i]^2 / V[1, j, i] + 0.5 * 9.81 * V[1, j, i]^2
        f[3] = V[2, j, i] * V[3, j, i] / V[1, j, i]
    end
end

function G!(g, V::Array{Float64, 3}, j::Int64, i::Int64)
    g[1] = V[3, j, i]
    if isapprox(V[1, j, i], 0., atol=eps(Float64), rtol=0)
        g[2] = 0
        g[3] = 0
    else
        g[2] = V[2, j, i] * V[3, j, i] / V[1, j, i]
        g[3] = V[3, j, i]^2 / V[1, j, i] + 0.5 * 9.81 * V[1, j, i]^2
    end
end

# Scheme
function LaxFriedrich!(V::Array{Float64, 3}, Vold::Array{Float64, 3}, lambdas::Array{Float64, 2}, dt, dx, dy)
    ny = size(V, 2)
    nx = size(V, 3)

    Cx = dt / dx
    Cy = dt / dy

    flux1 = @MVector [0., 0., 0.]
    flux2 = @MVector [0., 0., 0.]

    # moving along x axis
    for j=1:ny
        F!(flux1, Vold, j, 1)

        # left boundary
        for k=2:3
            V[k, j, 1] += Cx * (flux1[k] - lambdas[j, 1] * Vold[k, j, 1])
        end

        for i=1:nx-1
            F!(flux2, Vold, j, i+1)
            lambda = max(lambdas[j, i], lambdas[j, i+1])

            for k=1:3
                flux = 0.5 * Cx * ((flux2[k] + flux1[k]) - lambda * (Vold[k, j, i+1] - Vold[k, j, i]))
                V[k, j, i] -= flux
                V[k, j, i+1] += flux
            end
            flux1, flux2 = flux2, flux1
        end

        # right boundary
        for k=2:3
            V[k, j, nx] -= Cx * (flux1[k] + lambdas[j, nx] * Vold[k, j, nx])
        end
    end

    # moving along y-axis
    for i=1:nx
        G!(flux1, Vold, 1, i)

        # bottom boundary
        for k=2:3
            V[k, 1, i] += Cy * (flux1[k] - lambdas[1, i] * Vold[k, 1, i])
        end

        for j=1:ny-1
            G!(flux2, Vold, j+1, i)
            lambda = max(lambdas[j, i], lambdas[j+1, i])

            for k=1:3
                flux = 0.5 * Cy * ((flux2[k] + flux1[k]) - lambda * (Vold[k, j+1, i] - Vold[k, j, i]))
                V[k, j, i] -= flux
                V[k, j+1, i] += flux
            end
            flux1, flux2 = flux2, flux1
        end

        # top boundary
        for k=2:3
            V[k, ny, i] -= Cy * (flux1[k] + lambdas[ny, i] * Vold[k, ny, i])
        end
    end
end

# Computes time step satisfying the CFL condition
function UpdateCFL!(lambdas::Array{Float64, 2}, V::Array{Float64, 3}, dx::Float64, dy::Float64)
    nx = size(V, 3)
    ny = size(V, 2)

    M = 0
    for i=1:nx
        for j=1:ny
            if V[1, j, i] > 0
                lambdas[j, i] = sqrt(V[2, j, i]^2 + V[3, j, i]^2) / V[1, j, i] + sqrt(9.81 * V[1, j, i])
                M = max(M, lambdas[j, i])
            else
                lambdas[j, i] = 0.
            end
        end
    end

    return 0.5 * min(dx / M, dy / M)
end

# initialization
const domain = Domain(0., 1., 0., 1.)
const Nx = 512
const Ny = 512

# loop in time
const T = 2.0
const dx = (domain.x_end - domain.x_start) / Nx
const dy = (domain.y_end - domain.y_start) / Ny

let t=0
    V = Array{Float64, 3}(undef, 3, Nx, Ny)
    InitSolution!(V, domain)

    Vold = Array{Float64, 3}(undef, 3, Nx, Ny)
    lambdas = Array{Float64, 2}(undef, Nx, Ny)

    function update_to_time(tframe::Float64)
        while(t < tframe)
            dt = min(UpdateCFL!(lambdas, V, dx, dy), tframe - t)
            @. Vold = V
            LaxFriedrich!(V, Vold, lambdas, dt, dx, dy)
            t += dt
            INFO && @printf("\t%.5f / %.5f\n", t, tframe)
        end
    end

    if record_animation
        nframe_per_second = 1000
        nframe = convert(Int64, ceil(nframe_per_second * T)) + 1

        scene = Scene(resolution = (1920, 1080))
        x = range(domain.x_start, stop=domain.x_end, length=Nx)
        y = range(domain.y_start, stop=domain.y_end, length=Ny)

        surf = surface!(scene, x, y, view(V, 1, :, :))[end]
        record(scene, "./height2d.mp4", range(t, stop=T, length=nframe), framerate=30) do tframe
            println(tframe, " / ", T)
            update_to_time(tframe)
            surf[3] = view(V, 1, :, :)
        end
    else
        println("Initial time t = ", t)
        Profile.clear()
        @profile update_to_time(T)
        println("End of simulation, t = ", t)
    end
end

Profile.print(maxdepth=10)

max_stars_repo_name: danielcort/benchmarks-python-julia-c

max_issues_repo_name: danielcort/benchmarks-python-julia-c

max_forks_repo_name: danielcort/benchmarks-python-julia-c


text:
module _LinAlg

using GenericLinearAlgebra: svd!
using LinearAlgebra: norm, lmul!, rmul!, triu!, Givens, I, SVD, reflector!, reflectorApply!,
    QRPivoted, QRPackedQ

export tsvd, tsvd!, svd_jacobi, svd_jacobi!, rrqr, rrqr!

"""
Truncated rank-revealing QR decomposition with full column pivoting.

Decomposes a `(m, n)` matrix `A` into the product:

    A[:,piv] == Q * R

where `Q` is an `(m, k)` isometric matrix, `R` is a `(k, n)` upper
triangular matrix, `piv` is a permutation vector, and `k` is chosen such
that the relative tolerance `tol` is met in the equality above.
"""
function rrqr!(A::AbstractMatrix{T}; rtol=eps(T)) where {T<:AbstractFloat}
    # DGEQPF
    m, n = size(A)
    k = min(m, n)
    Base.require_one_based_indexing(A)

    jpvt = Vector(1:n)
    taus = Vector{T}(undef, k)
    swapcol = Vector{T}(undef, m)

    xnorms = map(norm, eachcol(A))
    pnorms = copy(xnorms)
    sqrteps = sqrt(eps(T))

    @inbounds for i in 1:k
        pvt = argmax(ii -> pnorms[ii], i:n)
        if i ≠ pvt
            jpvt[i], jpvt[pvt] = jpvt[pvt], jpvt[i]
            xnorms[pvt] = xnorms[i]
            pnorms[pvt] = pnorms[i]

            swapcol .= @view A[:, i]
            A[:, i] .= @view A[:, pvt]
            A[:, pvt] .= swapcol
        end

        tau_i = reflector!(@view A[i:end, i])
        taus[i] = tau_i
        reflectorApply!(
            (@view A[i:end, i]), tau_i, @view A[i:end, (i + 1):end]
        )

        # Lapack Working Note 176.
        for j in (i + 1):n
            temp = abs(A[i, j]) / pnorms[j]
            temp = max(zero(T), (one(T) + temp) * (one(T) - temp))
            temp2 = temp * abs2(pnorms[j] / xnorms[j])
            if temp2 < sqrteps
                recomputed = norm(@view A[(i + 1):end, j])
                pnorms[j] = recomputed
                xnorms[j] = recomputed
            else
                pnorms[j] *= sqrt(temp)
            end
        end

        # Since we did pivoting, R[i,:] is bounded by R[i,i], so we can
        # simply ignore all the other rows
        if abs(A[i, i]) < rtol * abs(A[1, 1])
            A[i:end, i:end] .= zero(T)
            taus[i:end] .= zero(T)
            k = i - 1
            break
        end
    end
    return QRPivoted{T,typeof(A)}(A, taus, jpvt), k
end

"""Truncated rank-revealing QR decomposition with full column pivoting."""
rrqr(A::AbstractMatrix{T}; rtol=eps(T)) where {T<:AbstractFloat} = rrqr!(copy(A); rtol)

"""Truncate RRQR result low-rank"""
function truncate_qr_result(qr::QRPivoted{T}, k::Integer) where {T}
    m, n = size(qr)
    0 ≤ k ≤ min(m, n) || throw(DomainError(k, "Invalid rank, must be in [0, $(min(m, n))]"))
    Qfull = QRPackedQ(view(qr.factors, :, 1:k), qr.τ[1:k])

    Q = lmul!(Qfull, Matrix{T}(I, m, k))
    R = triu!(qr.factors[1:k, :])
    return Q, R
end

"""
Truncated singular value decomposition.

Decomposes an `(m, n)` matrix `A` into the product:

    A == U * (s .* VT)

where `U` is a `(m, k)` matrix with orthogonal columns, `VT` is a `(k, n)`
matrix with orthogonal rows and `s` are the singular values, a set of `k`
nonnegative numbers in non-ascending order.  The SVD is truncated in the
sense that singular values below `tol` are discarded.
"""
function tsvd!(A::AbstractMatrix{T}; rtol=eps(T)) where {T<:AbstractFloat}
    # Perform RRQR of the m x n matrix at a cost of O(m*n*k), where k is
    # the QR rank (a mild upper bound for the true rank)
    A_qr, k = rrqr!(A; rtol)
    Q, R = truncate_qr_result(A_qr, k)

    # RRQR is an excellent preconditioner for Jacobi. One should then perform
    # Jacobi on RT
    RT_svd = svd!(Matrix(R'))

    # Reconstruct A from QR
    U = Q * RT_svd.V
    V = @view RT_svd.U[invperm(A_qr.p), :]
    s = RT_svd.S
    return SVD(U, s, V')
end

"""Truncated singular value decomposition."""
tsvd(A::AbstractMatrix{T}; rtol=eps(T)) where {T<:AbstractFloat} = tsvd!(copy(A); rtol)

#################################################################################################################
### Everything below is currently not used in tsvd. (GenericLinearAlgebra.svd! is used instead of svd_jacobi) ###
#################################################################################################################

"""
Compute Givens rotation `R` matrix that satisfies:

    [  c  s ] [ f ]     [ r ]
    [ -s  c ] [ g ]  =  [ 0 ]
"""
function givens_params(f::T, g::T) where {T<:AbstractFloat}
    # ACM Trans. Math. Softw. 28(2), 206, Alg 1
    if iszero(g)
        c, s = one(T), zero(T)
        r = f
    elseif iszero(f)
        c, s = zero(T), T(signbit(g))
        r = abs(g)
    else
        r = copysign(hypot(f, g), f)
        c = f / r
        s = g / r
    end
    return (c, s), r
end

"""
Apply Givens rotation to vector:

      [ a ]  =  [  c   s ] [ x ]
      [ b ]     [ -s   c ] [ y ]
"""
function givens_lmul((c, s)::Tuple{T,T}, (x, y)) where {T}
    a = c * x + s * y
    b = c * y - s * x
    return a, b
end

"""
Perform the SVD of upper triangular two-by-two matrix:

      [ f    g   ]  =  [  cu  -su ] [ smax     0 ] [  cv   sv ]
      [ 0    h   ]     [  su   cu ] [    0  smin ] [ -sv   cv ]

Note that smax and smin can be negative.
"""
function svd2x2(f::T, g::T, h::T) where {T<:AbstractFloat}
    # Code taken from LAPACK xSLAV2:
    fa = abs(f)
    ga = abs(g)
    ha = abs(h)
    if fa < ha
        # switch h <-> f, cu <-> sv, cv <-> su
        (sv, cv), (smax, smin), (su, cu) = svd2x2(h, g, f)
    elseif iszero(ga)
        # already diagonal, fa > ha
        smax, smin = fa, ha
        cv, sv = cu, su = one(T), zero(T)
    elseif fa < eps(T) * ga
        # ga is very large
        smax = ga
        if ha > one(T)
            smin = fa / (ga / ha)
        else
            smin = (fa / ga) * ha
        end
        cv, sv = one(T), f / g
        cu, su = one(T), h / g
    else
        # normal case
        fmh = fa - ha
        d = fmh / fa
        q = g / f
        s = T(2) - d
        spq = hypot(q, s)
        dpq = hypot(d, q)
        a = (spq + dpq) / T(2)
        smax, smin = abs(fa * a), abs(ha / a)

        tmp = (q / (spq + s) + q / (dpq + d)) * (one(T) + a)
        tt = hypot(tmp, T(2))
        cv = T(2) / tt
        sv = tmp / tt
        cu = (cv + sv * q) / a
        su = ((h / f) * sv) / a
    end
    return (cu, su), (smax, smin), (cv, sv)
end

"""
Perform the SVD of an arbitrary two-by-two matrix:

      [ a11  a12 ]  =  [  cu  -su ] [ smax     0 ] [  cv   sv ]
      [ a21  a22 ]     [  su   cu ] [    0  smin ] [ -sv   cv ]

Note that smax and smin can be negative.
"""
function svd2x2(a11::T, a12::T, a21::T, a22::T) where {T}
    abs_a12 = abs(a12)
    abs_a21 = abs(a21)
    if iszero(a21)
        # upper triangular case
        (cu, su), (smax, smin), (cv, sv) = svd2x2(a11, a12, a22)
    elseif abs_a12 > abs_a21
        # closer to lower triangular - transpose matrix
        (cv, sv), (smax, smin), (cu, su) = svd2x2(a11, a21, a12, a22)
    else
        # First, we use a givens rotation  Rx
        # [  cx   sx ] [ a11  a12 ] = [ rx  a12' ]
        # [ -sx   cx ] [ a21  a22 ]   [ 0   a22' ]
        (cx, sx), rx = givens_params(a11, a21)
        a11, a21 = rx, zero(rx)
        a12, a22 = givens_lmul((cx, sx), (a12, a22))

        # Next, use the triangular routine
        # [ f  g ]  =  [  cu  -su ] [ smax     0 ] [  cv   sv ]
        # [ 0  h ]     [  su   cu ] [    0  smin ] [ -sv   cv ]
        (cu, su), (smax, smin), (cv, sv) = svd2x2(a11, a12, a22)

        # Finally, update the LHS (U) transform as follows:
        # [  cx  -sx ] [  cu  -su ] = [  cu'  -su' ]
        # [  sx   cx ] [  su   cu ]   [  su'   cu' ]
        cu, su = givens_lmul((cx, -sx), (cu, su))
    end
    return (cu, su), (smax, smin), (cv, sv)
end

function jacobi_sweep!(U::AbstractMatrix, VT::AbstractMatrix)
    ii, jj = size(U)
    ii ≥ jj || throw(ArgumentError("matrix must be 'tall'"))
    size(VT, 1) == jj || throw(ArgumentError("U and VT must be compatible"))
    Base.require_one_based_indexing(U)
    Base.require_one_based_indexing(VT)

    offd = zero(eltype(U))
    @inbounds for i in 1:ii
        for j in (i + 1):jj
            # Construct the 2x2 matrix to be diagonalized
            Hii = sum(abs2, @view U[:, i])
            Hij = sum(k -> @inbounds(U[k, i] * U[k, j]), 1:ii)
            Hjj = sum(abs2, @view U[:, j])
            offd += abs2(Hij)

            # diagonalize
            (_, _), (_, _), (cv, sv) = svd2x2(Hii, Hij, Hij, Hjj)

            # apply rotation to VT
            rot = Givens(i, j, cv, sv)
            lmul!(rot, VT)
            rmul!(U, adjoint(rot))
        end
    end
    return sqrt(offd)
end

"""Singular value decomposition using Jacobi rotations."""
function svd_jacobi!(U::AbstractMatrix{T}; rtol=eps(T), maxiter=20) where {T}
    m, n = size(U)
    m ≥ n || throw(ArgumentError("matrix must be 'tall'"))
    Base.require_one_based_indexing(U)

    VT = Matrix(one(T) * I, n, n)
    Unorm = norm(@view U[1:n, 1:n])
    for _ in 1:maxiter
        offd = jacobi_sweep!(U, VT)
        offd < rtol * Unorm && break
    end

    s = norm.(eachcol(U))
    U ./= transpose(s)
    return SVD(U, s, VT)
end

"""Singular value decomposition using Jacobi rotations."""
function svd_jacobi(U::AbstractMatrix{T}; rtol=eps(T), maxiter=20) where {T}
    return svd_jacobi!(copy(U); rtol, maxiter)
end

end # module _LinAlg

ext: jl


text:
function parseKillerJSON(json_sums)
    sums = []
    for s in json_sums
        indices = Tuple[]
        for ind in s["indices"]
            push!(indices, tuple(ind...))
        end

        push!(sums, (result = s["result"], indices = indices, color = s["color"]))
    end
    return sums
end

function solve_killer_sudoku(filename; special = false, branch_strategy=:Auto)
    json_file = joinpath(dir, "benchmark/killer_sudoku/data/$filename")
    sums = parseKillerJSON(JSON.parsefile(json_file))

    m = CS.Optimizer(logging = [], time_limit = 20,
                    branch_strategy= branch_strategy)

    x = [[MOI.add_constrained_variable(m, MOI.Integer()) for i in 1:9] for j in 1:9]
    for r in 1:9, c in 1:9
        MOI.add_constraint(m, x[r][c][1], MOI.GreaterThan(1.0))
        MOI.add_constraint(m, x[r][c][1], MOI.LessThan(9.0))
    end

    for s in sums
        saf = MOI.ScalarAffineFunction{Float64}(
            [MOI.ScalarAffineTerm(1.0, x[ind[1]][ind[2]][1]) for ind in s.indices],
            0.0,
        )
        MOI.add_constraint(m, saf, MOI.EqualTo(convert(Float64, s.result)))
        if !special
            MOI.add_constraint(
                m,
                [x[ind[1]][ind[2]][1] for ind in s.indices],
                CS.AllDifferentSetInternal(length(s.indices)),
            )
        end
    end

    # sudoku constraints
    for r in 1:9
        MOI.add_constraint(
            m,
            MOI.VectorOfVariables([x[r][c][1] for c in 1:9]),
            CS.AllDifferentSetInternal(9),
        )
    end
    for c in 1:9
        MOI.add_constraint(
            m,
            MOI.VectorOfVariables([x[r][c][1] for r in 1:9]),
            CS.AllDifferentSetInternal(9),
        )
    end
    variables = [MOI.VariableIndex(0) for _ in 1:9]
    for br in 0:2
        for bc in 0:2
            variables_i = 1
            for i in (br * 3 + 1):((br + 1) * 3), j in (bc * 3 + 1):((bc + 1) * 3)
                variables[variables_i] = x[i][j][1]
                variables_i += 1
            end
            MOI.add_constraint(m, variables, CS.AllDifferentSetInternal(9))
        end
    end

    MOI.optimize!(m)
end


text:
export musescore, musescore_drumkey

using MusicManipulations, DefaultApplication

const MUSESCORE = @static Sys.iswindows() ? "MuseScore3" : "mscore"
const MUSESCORE_EXISTS = [false]

function test_musescore()
    if !MUSESCORE_EXISTS[1]
        r = try
            r = run(`$(MUSESCORE) -v`)
        catch
            false
        end
        if r == false || ((typeof(r) == Base.Process) && r.exitcode != 0)
            throw(SystemError(
            """
            The command `$(MUSESCORE) -v` did not run, which probably means that
            MuseScore is not accessible from the command line.
			Please first install MuseScore
            on your computer and then add it to your PATH."""
            ))
        end
    end
    global MUSESCORE_EXISTS[1] = true
end

"""
    musescore(file, notes | midi; keywords...)
Use the open source software "MuseScore" to create a score and save the given
`notes` (or an entire `midi` file) to `file`, which can be either a `.pdf` or
a `.png`.If given a `.png` the actual file name will end with `-1`, `-2` etc.
for each page of the score.

MuseScore must be accessible from the command line
for this to work (add the path to `MuseScore.exe` to your PATH environment
variable (ask Google)).

The keyword `display = true` will also display the created `file`.
Keyword `rmmidi = true` deletes the MIDI file that has to be created inbetween
score convertion (use `false` to keep it).
MuseScore is run by default with commands `-n -T 10 -r 1200`, but you can
change the keyword `c` to be whatever command (enclosed in backticks) you want, see the
[command line options](https://musescore.org/en/handbook/3/command-line-options)
of MuseScore for details.
"""
function musescore(file, notes;
	display = true, rmmidi = true, c = `-n -T 20 -r 1200`
	)

    lowercase(file[end-3:end]) ∈ (".png", ".pdf") || error("file must be .pdf or .png.")
    MUSESCORE_EXISTS[1] || test_musescore()

    midiname = file[1:end-4]*".mid"
	midi = writeMIDIFile(midiname, notes)

    if file[end-3:end] == ".png"
        cmd = `$MUSESCORE $(c) -o $(file) $(midiname)`
        muspng = file[1:end-4]*"-1.png"
    else
        cmd = `$MUSESCORE $(c) -o $(file) $(midiname)`
        muspng = file
    end
	run(cmd)
    rmmidi && rm(midiname)
    display && DefaultApplication.open(muspng)
end

"""
	musescore_drumkey

A dictionary that given the drum instrument as a string it returns the
MIDI pitch that MuseScore uses. See
https://musescore.org/sites/musescore.org/files/General%20MIDI%20Standard%20Percussion%20Set%20Key%20Map.pdf
for a more complete list.

Notice that to import a MIDI file into MuseScore and the drumset to be the selected
instrument, you should write the MIDI notes to channel `9`.
"""
const musescore_drumkey = begin
	a =
	Dict(
	"Acoustic Bass Drum" => "B1",
	"Acoustic Snare" => "D2",
	"Side Stick" => "C#2",
	"Closed Hi-Hat" => "F#2",
	"Open Hi-Hat" => "A#2",
	"Pedal Hi-Hat" => "G#2",
	"Ride Cymbal" => "D#3",
	"Ride Bell" => "F3",
	"Low-Mid Tom" => "B2",
	"Cowbell" => "G#3",
	"Tambourine" => "F#3",
	"High Floor Tom" => "G2",
	"Low Floor Tom" => "F2",
	"Crash Cymbal 1" => "C#3",
	"Crash Cymbal 2" => "A3"
	)
	b = Dict((k, name_to_pitch(a[k])) for k in keys(a))
end


text:
module mfesp1
using StaticArrays
using Elfel
using Elfel.RefShapes: RefShapeTriangle, manifdim, RefShapeInterval
using Elfel.FElements: FE, refshape, FEH1_T3
using Elfel.FElements: bfun, bfungradpar
using Elfel.FESpaces: FESpace, ndofs, setebc!, nunknowns, doftype, nunknowns
using Elfel.FESpaces: numberfreedofs!, numberdatadofs!
using Elfel.FEIterators: FEIterator
using MeshCore
using MeshSteward: Mesh, load, nspacedims, baseincrel
using Test
function test()
    fe = FEH1_T3()
    mesh = load(Mesh(), "qmesh.mesh")
    fesp = FESpace(Float64, mesh, fe)

    @test doftype(fesp) == Float64

    refc = [[1, 2, 5],
    [1, 5, 4],
    [4, 5, 8],
    [4, 8, 7],
    [7, 8, 11],
    [7, 11, 10],
    [2, 3, 6],
    [2, 6, 5],
    [5, 6, 9],
    [5, 9, 8],
    [8, 9, 12],
    [8, 12, 11], ]
    it = FEIterator(fesp)
    k = 1
    for c in it
        @test isapprox(c._nodes,  refc[k])
        k = k + 1
    end
    @test length(it) == 12

    # @show summary(mesh)

    numberfreedofs!(fesp)
    @test ndofs(fesp) == 12
    @test nunknowns(fesp) == 12

    for i in [1, 4, 7, 10]
        setebc!(fesp, 0, i, 1, 0.0)
    end
    numberfreedofs!(fesp)
    @test ndofs(fesp) == 12
    @test nunknowns(fesp) == 8

    it = FEIterator(fesp)
    for el in it
         el._dofs
    end
    # sdim = nspacedims(femesh)
    # mdim = manifdim(femesh)

    # geom = geomattr(femesh)

    # el = 4
    # gradNpar = bfungradpar(fespace.fe, [1/3, 1/3])

    # conn = connectivity(femesh)
    # for c in conn
    #     J = SMatrix{sdim, mdim}(zeros(sdim, mdim))
    #     for j in 1:length(c)
    #         J += geom[c[j]] * gradNpar[j]
    #     end
    # end

    # c = iterate(conn, 1)[1]
    # J = SMatrix{sdim, mdim}(zeros(sdim, mdim))
    # for j in 1:length(c)
    #     J += geom[c[j]] * gradNpar[j]
    # end
    # @test isapprox(J, [0.5 0.5; 0.0 0.3333333333333333])
end
end
using .mfesp1
mfesp1.test()

module mfesp2
using StaticArrays
using Elfel
using Elfel.RefShapes: RefShapeTriangle, manifdim, RefShapeInterval
using Elfel.FElements: FE, refshape, FEH1_T3, shapedesc
using Elfel.FElements: FEH1_T3_BUBBLE
using Elfel.FElements: bfun, bfungradpar
using Elfel.FElements: nfeatofdim, ndofperfeat, ndofsperel
using Elfel.FESpaces: FESpace, ndofs, setebc!, nunknowns, doftype, nunknowns
using Elfel.FESpaces: numberfreedofs!, numberdatadofs!, dofnum
using Elfel.FESpaces: edofmdim, edofbfnum, edofcompnt
using Elfel.FEIterators: FEIterator
using MeshCore: ir_identity
using MeshSteward: Mesh, load, nspacedims, baseincrel, attach!
using Test
function test()
    fe = FEH1_T3()
    mesh = load(Mesh(), "qmesh.mesh")
    fesp = FESpace(Float64, mesh, fe, 3)

    @test doftype(fesp) == Float64

    @test shapedesc(fesp.fe).name == "T3"
    @test refshape(fesp.fe) == Elfel.RefShapes.RefShapeTriangle

    @test nfeatofdim(fesp.fe, 0) == 3
    @test ndofperfeat(fesp.fe, 0) == 1

    @test ndofsperel(fesp.fe) == 3

    @test isapprox(bfun(fesp.fe, [1/3, 1/3]), [0.3333333333333334, 0.3333333333333333, 0.3333333333333333])

    attach!(mesh, ir_identity(baseincrel(mesh)))
    fesp = FESpace(Float64, mesh, FEH1_T3_BUBBLE(), 1)

    @test nfeatofdim(fesp.fe, 0) == 3
    @test nfeatofdim(fesp.fe, 1) == 3
    @test nfeatofdim(fesp.fe, 2) == 1
    @test nfeatofdim(fesp.fe, 3) == 0
    @test ndofperfeat(fesp.fe, 0) == 1
    @test ndofperfeat(fesp.fe, 1) == 0
    @test ndofperfeat(fesp.fe, 2) == 1
    @test ndofperfeat(fesp.fe, 3) == 0

    @test ndofsperel(fesp.fe) == 4

    @test isapprox(bfun(fesp.fe, [1/3, 1/3]), [0.3333333333333334, 0.3333333333333333, 0.3333333333333333, 0.03703703703703704])

    fesp = FESpace(Float64, mesh, FEH1_T3_BUBBLE(), 2)

    @test nfeatofdim(fesp.fe, 0) == 3
    @test nfeatofdim(fesp.fe, 1) == 3
    @test nfeatofdim(fesp.fe, 2) == 1
    @test nfeatofdim(fesp.fe, 3) == 0
    @test ndofperfeat(fesp.fe, 0) == 1
    @test ndofperfeat(fesp.fe, 1) == 0
    @test ndofperfeat(fesp.fe, 2) == 1
    @test ndofperfeat(fesp.fe, 3) == 0

    @test ndofsperel(fesp.fe) == 4

    @test isapprox(bfun(fesp.fe, [1/3, 1/3]), [0.3333333333333334, 0.3333333333333333, 0.3333333333333333, 0.03703703703703704])

    emdim = Int64[]
    bfnum = Int64[]
    compnt = Int64[]
    bfn = 1
    for m in 0:1:3
        if ndofperfeat(fesp.fe, m) > 0
            for i in 1:nfeatofdim(fesp.fe, m)
                for k in 1:ndofperfeat(fesp.fe, m)
                    for j in 1:fesp.nfecopies
                        push!(emdim, m)
                        push!(bfnum, bfn)
                        push!(compnt, j)
                    end
                    bfn += 1
                end
            end
        end
    end
    @test edofmdim(fesp) == emdim
    @test edofbfnum(fesp) == bfnum
    @test edofcompnt(fesp) == compnt

    @test dofnum(fesp, 0, 1, 1) == 0

    true
end
end
using .mfesp2
mfesp2.test()


text:
@testset "solveMAP" begin

M  = [1 1 1; 1 1 0; 1 0 1; 1 0 0; 0 1 1; 0 1 0; 0 0 1;]
(m,n) = size(M)
R     = StrainRecon.sdiag(ones(m)*100);
atrue = [.6; 0.3; 0.1]
Data  = M*atrue

param = getStrainReconParam(Data,n,R=R)

Mmap,amap = solveMAP(param,maxTry=100,maxIterQP=10)
@test norm(amap - atrue)/norm(atrue) < 1e-8
@test norm(Mmap*amap-Data) < 1e-7

Msqp,Asqp = solveMAP(param,storeInterm=true,maxIterQP=20)
# Asqp is n-by-history of each BCD iteration-by-MaxTry
# Msqp is m-by-n-by-history of each BCD iteration-by-MaxTry

Jc = zeros(size(Asqp,3)) # Misfit for MaxTry number of iterations

for k=1:size(Asqp,3)
    # for each Asqp iteration, find column where norm(w_old - w_new) = 0, use column number - 1 to find frame for Msqp
    iter = size(Asqp,2)
    for i = 1:size(Asqp,2)
        # println(i)
        if sum(Asqp[:,i,k]) == 0
            iter = i-1
            # println("found column of all zeros ", i-1)
            break
        end
    end
    @test all(Asqp[:,iter,k].>0)
	Jc[k] = misfit(param,Msqp[:,:,iter-1,k],Asqp[:,iter,k])[1]
end
@test all(diff(Jc).>=0)
@test norm(sum(Asqp[:,1],1)-1,Inf) < 1e-7 # Asqp is n-by-history of iterations
# @test all(Asqp[:,:].>0)
end


text:
#DotsGram SepalLength
plot(iris, x=:SepalLength, y=:PetalLength, color=:Species)
#DotsGram PetalWidth
plot(iris, x=:PetalWidth, y=:PetalLength, color=:Species)
#BoxPlot
plot(iris, x=:Species, y=:PetalLength, Geom.boxplot)
#Histogram PetalLength
plot(iris, x=:PetalLength, color=:Species, Geom.histogram)
#Histogram PetalWidth
plot(iris, x=:PetalWidth, color=:Species, Geom.histogram)

CSV.write("test_data.csv", test_data)
CSV.read("test_data.csv")
Feather.write("test_data.feather", test_data)
Feather.read("test_data.feather")

test_data = iris[rand(150) .<= 0.1, [:PetalLength, :PetalWidth, :Species]]


text:
using MPIMeasurements
using Plots

#include("config.jl")

#scanner = MPIScanner(conf)
daq = getDAQ(scanner)

params = toDict(daq.params)

params["studyName"] = "Test"
params["studyDescription"] = "A very cool measurement"
params["scannerOperator"] = "Tester"
params["dfStrength"] = [1e-3]
params["dfPhase"] = [pi/2]
params["dfWaveform"] = "SINE"
params["acqNumFrames"] = 1
params["acqNumAverages"] = 1
params["controlPhase"] = false #true TODO
params["acqFFSequence"] = "None"

u = measurement(daq, params)

@info size(u)

p = plot(vec(u[:,1,1,1]))
savefig(p, joinpath(imgdir, "im1.png"))

scanner = MPIScanner(conf) #what???
daq = getDAQ(scanner)
params["dfWaveform"] = "TRIANGLE"
params["controlPhase"] = false
u = measurement(daq, params)

p = plot(vec(u[:,1,1,1]))
savefig(p, joinpath(imgdir, "im2.png"))


text:
#=
see:
https://github.com/hfawzi/cvxquad/blob/master/examples/entanglement_assisted_capacity.m
listing 2 in "Efficient optimization of the quantum relative entropy"
by H. Fawzi and O. Fawzi
=#

using SparseArrays

struct EntanglementAssisted{T <: Real} <: ExampleInstanceJuMP{T}
    nb::Int
    ne::Int
end

function build(inst::EntanglementAssisted{T}) where {T <: Float64}
    gamma = 0.2
    ampl_damp = [
        1 0;
        0 sqrt(gamma);
        0 sqrt(1 - gamma);
        0 0;
        ]
    ampl_dim = 4
    na = 2
    nb = inst.nb
    ne = inst.ne
    @assert nb * ne == ampl_dim
    rt2 = sqrt(T(2))
    sa = Cones.svec_length(ampl_dim)
    sb = Cones.svec_length(nb)

    model = JuMP.Model()
    JuMP.@variables(model, begin
        ρ[1:na, 1:na], PSD
        cond_epi
        qe_epi
    end)

    Q1 = Symmetric(ampl_damp * ρ * ampl_damp')
    Q2 = zeros(JuMP.AffExpr, nb * ne, nb * ne)
    kron!(Q2, I(nb), partial_trace(Q1, 1, [nb, ne]))
    Q3 = partial_trace(Q1, 2, [nb, ne])
    Q1_vec = Cones.smat_to_svec!(zeros(JuMP.AffExpr, sa), Q1, rt2)
    Q2_vec = Cones.smat_to_svec!(zeros(JuMP.AffExpr, sa), Q2, rt2)
    Q3_vec = Cones.smat_to_svec!(zeros(JuMP.AffExpr, sb), Q3, rt2)
    RE_cone = Hypatia.EpiTrRelEntropyTriCone{T}(1 + 2 * sa)
    E_cone = Hypatia.EpiPerSepSpectralCone{T}(Cones.NegEntropySSF(),
        Cones.MatrixCSqr{T, T}, nb)

    JuMP.@constraints(model, begin
        tr(ρ) == 1
        vcat(cond_epi, Q1_vec, Q2_vec) in RE_cone
        vcat(qe_epi, 1, Q3_vec) in E_cone
    end)

    JuMP.@objective(model, Max, (cond_epi + qe_epi) / -log(T(2)))

    return model
end

# partial trace of Q over system i given subsystem dimensions subs
function partial_trace(
    Q::Symmetric,
    i::Int,
    subs::Vector{Int},
    )
    @assert 1 <= i <= length(subs)
    @assert size(Q, 1) == prod(subs)
    return sum(partial_trace_j(j, Q, i, subs) for j in 1:subs[i])
end

function partial_trace_j(
    j::Int,
    Q::Symmetric,
    i::Int,
    subs::Vector{Int},
    )
    X1 = sparse(I, 1, 1)
    X2 = copy(X1)
    for (k, dim) in enumerate(subs)
        if k == i
            spej = spzeros(Bool, dim, 1)
            spej[j] = 1
            X1 = kron(X1, spej')
            X2 = kron(X2, spej)
        else
            spI = sparse(I, dim, dim)
            X1 = kron(X1, spI)
            X2 = kron(X2, spI)
        end
    end
    return X1 * Q * X2
end


text:
module DichotomicLocalPolytope

# """
# Module: DichotomicLocalPolytope
#
# Hard-coded computation of local polytope vertices for the CHSH scenario with
# and without auxiliary communication.
# """

using QBase: QMath
# """
# bipartite_no_signaling_strategies(;subspace="p16"):
#   Returns the set of strategies for dichotomic, bi-partite bell test. The
#   joint strategies for Alice and Bob are represented as the kronecker product
#   of their individual strategies α(x)⊗β(y).
#
# Input: subspace, string, the behavior subspace for which results are returned
#        defaults to "p16", the most general and high dimensional bipartite space.
#
# Output: An array of joint strategies which map input x⊗y → a⊗b
# """
function bipartite_no_signaling_strategies(;subspace="p16")
    # Assumption: Alice and Bob have dichotomic setups
    s = dichotomic_strategies()
    n = size(s)[1]

    no_signaling_strategies = [] # initialize empty array

    # Alice and Bob each choose a strategy from the dichotomic strategies.
    for α in s
        for β in s
            if subspace == "p16"
                push!( no_signaling_strategies, kron(α, β) )
            elseif subspace == "p12"
                push!( no_signaling_strategies, kron(α, β)[1:3,:] )
            end
        end
    end

    no_signaling_strategies
end

# """
# deterministic_no_signaling_behaviors(;subspace="p16"):
#   Computes the deterministic behaviors (vertices) of the local polytope for
#   a bi-partitie, dichotomic setup for which no-signaling is applied.
#
# Input: subspace, string, the behavior subspace for which results are returned
#        defaults to "p16", the most general and high dimensional bipartite space.
#
# Output: Column Vectors, which represent deterministic behaviors in the 16-dim
#   basis described by { P(ab|xy) } ∀ a,b,x,y.
# """
function deterministic_no_signaling_behaviors(;subspace="p16")
    # Assumptions: 2-parties with dichotomic setups
    if subspace == "p16"
        s = bipartite_no_signaling_strategies()
    elseif subspace == "p12"
        s = bipartite_no_signaling_strategies(subspace="p12")
    end

    n = size(s)[1]

    behaviors = []
    # looping through distinct strategies
    for i in 1:n
        # initializing an empty probabilistic behavior
        if subspace == "p16"
            p = zeros(Int64,16)
        elseif subspace == "p12"
            p = zeros(Int64,12)
        end

        # columns of identity matrix represent xy ∈ {00, 01, 10, 11}
        identity = QMath.id(size(s[i])[2])

        # Strategy columns represent outputs (a⊗b) for given input (x⊗y). We
        # compute the 16-dim behavior by (a⊗b)⊗(x⊗y) → P(ab|xy)
        for j in 1:size(s[i])[2]
            p = p + kron(s[i][:,j],identity[:,j])
        end

        push!(behaviors,p)
    end

    behaviors
end

# """
# Description:
#   Computes the vertices of the local polytope which bounds the system where
#   Alice sends 1-bit of communication to Bob.
#
# Input:
#   subspace: String (optional), the behavior subspace in which we represent the
#             deterministic behaviors. default, 16-dim
#
# Output:
#   an array of deterministic vertices represented in the specified subspace
# """
function deterministic_fixed_direction_signaling_behaviors(;subspace="p16")
    # Assumption: signaling goes from Alice to Bob

    ds = dichotomic_strategies()

    ts = two_input_strategies()

    binary_inputs = [
        [1;0], # zero
        [0;1]  # one
    ]

    deterministic_behaviors = []
    for α in ds
        for β in ts

            if subspace == "p16"
                behavior = zeros(Int64,16)
            elseif subspace == "p12"
                behavior = zeros(Int64,12)
            elseif subspace == "p10"
                behavior = zeros(Int64,10)
            end

            for x in binary_inputs

                for y in binary_inputs

                    xy = kron(x,y)
                    input = kron(xy,x)

                    if subspace == "p16"
                        ab = kron(α,β)*input
                        behavior += kron(ab,xy)
                    elseif subspace == "p12"
                        ab = kron(α,β)[1:3,:]*input
                        behavior += kron(ab,xy)
                    elseif subspace == "p10"
                        ab = kron(α,β)[[1,3],:]*input

                        behavior += cat(zeros(Int64,2),kron(ab,xy),dims=1)
                    end
                end
                if subspace == "p10"
                    a = α[1,:]'*x # α gets cast as a 1-D column vector
                    behavior += cat(kron(a,x),zeros(Int64,8),dims=1)
                end

            end
            push!(deterministic_behaviors, behavior)
        end
    end

    deterministic_behaviors
end

# """
# deterministic_one_way_signaling_behaviors(;subspace="p16"):
#   Computes the deterministic behaviors (vertices) of the local polytope by
#   executing a finite set of deterministic strategies for all inputs. In this
#   computation it is assumed that 1-bit of information is passed either from
#   Alice to Bob or from Bob to Alice. The case is also considered where there
#   is no-signaling between Alice and Bob.
#
# Input: subspace, string, the behavior subspace for which results are returned
#        defaults to "p16", the most general and high dimensional bipartite space.
#
# Output: An array of p12 behaviors
# """
function deterministic_one_way_signaling_behaviors(;subspace="p16")
    ds = dichotomic_strategies()

    # filter out strategies that do not depend on second input
    ts = filter(f -> (
            (f*[1;0;0;0] != f*[0;1;0;0]) | (f*[0;0;1;0] != f*[0;0;0;1])
        ),
        two_input_strategies()
    )

    binary_inputs = [[1;0],[0;1]]
    deterministic_behaviors = deterministic_no_signaling_behaviors(subspace=subspace)
    # case for when Alice receives 1-bit from Bob
    for α in ts
        for β in ds
            if subspace == "p16"
                strategy_behavior = zeros(Int64,16)
            elseif subspace == "p12"
                strategy_behavior = zeros(Int64,12)
            end

            for x in binary_inputs
                for y in binary_inputs
                    xy = kron(x,y)
                    input = kron(xy,y)

                    if subspace == "p16"
                        ab = kron(α,β)*input
                    elseif subspace == "p12"
                        ab = kron(α,β)[1:3,:]*input
                    end

                    strategy_behavior += kron(ab,xy)
                end
            end
            push!(deterministic_behaviors,strategy_behavior)
        end
    end

    # case for when Bob receives 1-bit from Alice
    for α in ds
        for β in ts
            if subspace == "p16"
                strategy_behavior = zeros(Int64,16)
            elseif subspace == "p12"
                strategy_behavior = zeros(Int64,12)
            end

            for x in binary_inputs
                for y in binary_inputs
                    xy = kron(x,y)
                    input = kron(xy,x)

                    if subspace == "p16"
                        ab = kron(α,β)*input
                    elseif subspace == "p12"
                        ab = kron(α,β)[1:3,:]*input
                    end

                    strategy_behavior += kron(ab,xy)
                end
            end
            push!(deterministic_behaviors,strategy_behavior)
        end
    end

    deterministic_behaviors
end

# """
# dichotomic_strategies():
#   Returns the set of deterministic strategies that map a binary
#   input to a binary output.
#
# Input: None.
#
# Output: An array of 2x2 matrices that represent all possible strategies.
# """
function dichotomic_strategies()
    I2 = [1 0; 0 1]
    PC0 = [1 1; 0 0]
    PC1 = [0 0; 1 1]
    P2 = [0 1; 1 0]

    [I2, PC0, PC1, P2]
end

# """
# two_input_strategies():
#   Returns the set of deterministic strategies that map two binary inputs to
#   a single binary output.
#
# Input: None.
#
# Output: An array of 2x4 matrices that represent all possible strategies
# """
function two_input_strategies()
    s = dichotomic_strategies()
    n = size(s)[1]

    two_input_strategies = []
    for i in 1:n
        for j in 1:n
            push!(two_input_strategies, cat(s[i],s[j],dims=2))
        end
    end

     two_input_strategies
end

end


text:
@reexport module WindyGridWorld

export WindyGridWorldEnv

using ReinforcementLearningBase

const NX = 7
const NY = 10
const Wind = [CartesianIndex(w, 0) for w in [0, 0, 0, -1, -1, -1, -2, -2, -1, 0]]
const StartPosition = CartesianIndex(4, 1)
const Goal = CartesianIndex(4, 8)
const Actions = [
    CartesianIndex(0, -1),  # left
    CartesianIndex(0, 1),   # right
    CartesianIndex(-1, 0),  # up
    CartesianIndex(1, 0),   # down
]

const LinearInds = LinearIndices((NX, NY))

Base.@kwdef mutable struct WindyGridWorldEnv <: AbstractEnv
    position::CartesianIndex{2} = StartPosition
end

RLBase.get_observation_space(env::WindyGridWorldEnv) = DiscreteSpace(length(LinearInds))
RLBase.get_action_space(env::WindyGridWorldEnv) = DiscreteSpace(length(Actions))

function (env::WindyGridWorldEnv)(a::Int)
    p = env.position + Wind[env.position[2]] + Actions[a]
    p = CartesianIndex(min(max(p[1], 1), NX), min(max(p[2], 1), NY))
    env.position = p
    nothing
end

RLBase.observe(env::WindyGridWorldEnv) =
    (
        state = LinearInds[env.position],
        terminal = env.position == Goal,
        reward = env.position == Goal ? 0.0 : -1.0,
    )

function RLBase.reset!(env::WindyGridWorldEnv)
    env.position = StartPosition
    nothing
end

end


text:
# PART 5 – ASSESSMENT OF LOCAL METAL LOSS

# Determine Asessment Applicability
#Determine the assessment applicability
# @doc DesignCodeCriteria
# @doc MaterialToughness
# @doc CyclicService
# @doc Part5ComponentType
print("Begin -- Assessment Applicability and Component Type Checks\n")
creep_range = CreepRangeTemperature("Carbon Steel (UTS ≤ 414MPa (60 ksi))"; design_temperature=100.0, units="lbs-in-psi")
design = DesignCodeCriteria("ASME B31.3 Piping Code")
toughness = MaterialToughness("Certain")
cyclic = CyclicService(100, "Meets Part 14")
x = Part5ComponentType("Straight Section of Piping, Elbow or Bend - No Structural Attachments", vessel_orientation="horizontal", material="Carbon and Low Alloy Steels", D=24.75,Lss=120.0,H=0.0, NPS=10.0, design_temperature=100.0, units="lbs-in-psi")
part5_applicability = Part5AsessmentApplicability(x,design,toughness,cyclic,creep_range)

# For all assessments - determine the inspection data grid
M1 = [0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300]
M2 = [0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.100, 0.220, 0.280, 0.250, 0.240, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300]
M3 = [0.300, 0.300, 0.300, 0.300, 0.300, 0.215, 0.255, 0.215, 0.145, 0.275, 0.170, 0.240, 0.250, 0.250, 0.280, 0.290, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300]
M4 = [0.300, 0.300, 0.300, 0.300, 0.300, 0.170, 0.270, 0.190, 0.190, 0.285, 0.250, 0.225, 0.275, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300]
M5 = [0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300]
M6 = [0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300, 0.300]
CTPGrid = hcat(M6,M5,M4,M3,M2,M1) # build in descending order
CTPGrid = rotl90(CTPGrid) # rotate to correct orientation

# Input variables
    annex2c_tmin_category = "Straight Pipes Subject To Internal Pressure" # ["Cylindrical Shell","Spherical Shell","Hemispherical Head","Elliptical Head","Torispherical Head","Conical Shell","Toriconical Head","Conical Transition","Nozzles Connections in Shells",
    # "Junction Reinforcement Requirements at Conical Transitions","Tubesheets","Flat head to cylinder connections","Bolted Flanges","Straight Pipes Subject To Internal Pressure","Boiler Tubes","Pipe Bends Subject To Internal Pressure",
    # "MAWP for External Pressure","Branch Connections","API 650 Storage Tanks"]
    equipment_group = "piping" # "vessel", "tank"
    flaw_location = "internal" # "External","Internal"
    pipe_code = "B31.8" # "B31.3", "B31.8"
    FCA_string = "internal"
    metal_loss_categorization = "LTA" # "LTA" or "Groove-Like Flaw"
    units = "lbs-in-psi" # "lbs-in-psi" or "nmm-mm-mpa"
    remaining_life = true
    tmm_forcing = true
    tmm = .949
    tnom = 1.0 # nominal or furnished thickness of the component adjusted for mill undertolerance as applicable.
    trd = 1.0 # uniform thickness away from the local metal loss location established by thickness measurements at the time of the assessment.
    FCAml = 0.1 # Future Corrosion Allowance applied to the region of metal loss.
    FCA = 0.0 # Future Corrosion Allowance applied to the region away from the metal loss (see Annex 2C, paragraph 2C.2.8).
    LOSS = 0.0 #the amount of uniform metal loss away from the local metal loss location at the time of the assessment.
    Do = 8.625 # Outside Diameter
    D = Do - 2*(tnom) # Inside Dia.
    P = 8139.13 # internal design pressure.
    S = 65000.0 # code allowable stress note if using B31.8 its SMYS - check design code criteria
    E = 1.0 # weld joint efficiency or quality factor from the original construction code, if unknown use 0.7.
    F = 0.6 # B31.8 design factor
    T = 1.0 # temperature factor B31.8
    MA = 0.0 # mechanical allowances (thread or groove depth); for threaded components, the nominal thread depth (dimension h of ASME B.1.20.1) shall apply.
    Yb31 = 0.4 # coefficient from ASME B31 Piping codes used for determining the pipe wall thickness, the coefficient can be determined from the following table that is valid for tmin < Do / 6 Annex 2C .
    t = trd - LOSS - FCA # thickness of the shell or pipe adjusted for mill tolerance, LOSS and FCA , or cylinder thickness at a conical transition for a junction reinforcement calculation adjusted for mill tolerance, LOSS and FCA , as applicable.
    t = 1.0
    tsl = 0.0 # supplemental thickness for mechanical loads other than pressure that result in longitudinal stress; this thickness is usually obtained from the results of a weight case in a stress analysis of the piping system (see paragraph 2C.2.7).
    spacings = 0.5 # spacings determine by visual inspection to adequately ccategorizse the corrosion -----------+ may add to CTP_Grid function for plotting purposes
    # Flaw dimensions
    s = 20.0 # longitudinal extent or length of the region of local metal loss based on future corroded thickness,
    c = 6.7  # circumferential extent or length of the region of local metal loss (see Figure 5.2 and Figure 5.10), based on future corroded thickness, tc .
    Ec = 1.0 # circumferential weld joint efficiency. note if damage on weld see # 2C.2.5 Treatment of Weld and Riveted Joint Efficiency, and Ligament Efficiency
    El = 1.0 # longitudinal weld joint efficiency. note if damage on weld see # 2C.2.5 Treatment of Weld and Riveted Joint Efficiency, and Ligament Efficiency
    RSFa = 0.9 # remaining strength factor - consult API 579 is go lower than 0.9

    # For all assessments determine far enough from structural discontinuity
    # Flaw-To-Major Structural Discontinuity Spacing
    L1msd = [12.0] # distance to the nearest major structural discontinuity.
    L2msd = [12.0] # distance to the nearest major structural discontinuity.
    L3msd = [12.0] # distance to the nearest major structural discontinuity.
    L4msd = [12.0] # distance to the nearest major structural discontinuity.
    L5msd = [12.0] # distance to the nearest major structural discontinuity.
    Lmsd = minimum([L1msd,L2msd,L3msd,L4msd,L5msd])[1]
    if (Lmsd >= (1.8*(sqrt(D*(trd - LOSS - FCA)))))
        print("Satisfied - Flaw is located far enough from structural discontinuity\n")
        lmsd_satisfied = 1
    else
        print("Not satisfied - Flaw is too close to a structural discontinuity - Conduct a level 3 assessment\n")
        lmsd_satisfied = 0
    end

    # Groove Like Flaw dimensions
    gl = .05 # length of the Groove-Like Flaw based on future corroded condition.
    gw = .4 # width of the Groove-Like Flaw based on future corroded condition.
    gr = 0.1 # radius at the base of a Groove-Like Flaw based on future corroded condition.
    β = 40.0 # see (Figure 5.4) :: orientation of the groove-like flaw with respect to the longitudinal axis or a parameter to compute an effective fracture toughness for a groove being evaluated as a crack-like flaw, as applicable.

# Perform level 1 assessment
if (part5_applicability[1] == 1 && lmsd_satisfied == 1) # begin level 1 assessment
    #let part_5_lta_output = Array{Any,2},
    part_5_lta_output = Part5LTALevel1(CTPGrid; remaining_life=remaining_life,tmm_forcing=tmm_forcing, tmm=tmm, annex2c_tmin_category=annex2c_tmin_category, pipe_code=pipe_code,equipment_group=equipment_group, flaw_location=flaw_location, FCA_string=FCA_string, metal_loss_categorization=metal_loss_categorization, units=units, Lmsd=Lmsd,tnom=tnom,
        trd=trd, FCA=FCA, FCAml=FCAml, LOSS=LOSS, Do=Do, D=D, P=P, S=S, E=E, MA=MA, Yb31=Yb31, t=t,tsl=tsl, F=F, T=T, spacings=spacings, s=s, c=c, El=El, Ec=Ec, RSFa=RSFa, gl=gl, gw=gw, gr=gr,β=β)
    #end # let end
    part_5_lta_output

elseif (part5_applicability[1] == 0 && lmsd_satisfied == 0)
    print("Level 1 Criteria Not Met - Perform Level 2 or 3 as applicable")
elseif (part5_applicability[1] == 1 && lmsd_satisfied == 0)
    print("Level 1 Criteria Not Met - Perform Level 2 or 3 as applicable")
elseif (part5_applicability[1] == 0 && lmsd_satisfied == 1)
    print("Level 1 Criteria Not Met - Perform Level 2 or 3 as applicable")
end

# Perform level 2 assessment
if (part5_applicability[2] == 1 && lmsd_satisfied == 1) # begin level 2 assessment
    #let part_5_lta_output = Array{Any,2},
    part_5_lta_output = Part5LTALevel2(CTPGrid; annex2c_tmin_category=annex2c_tmin_category, pipe_code=pipe_code, equipment_group=equipment_group, flaw_location=flaw_location, FCA_string=FCA_string, metal_loss_categorization=metal_loss_categorization, units=units, Lmsd=Lmsd, tnom=tnom,
        trd=trd, FCA=FCA, FCAml=FCAml, LOSS=LOSS, Do=Do, D=D, P=P, S=S, E=E, MA=MA, Yb31=Yb31, tsl=tsl, t=t,F=F, T=T, spacings=spacings, s=s, c=c, El=El, Ec=Ec, RSFa=RSFa, gl=gl, gw=gw, gr=gr,β=β)
    #end # let end
elseif (part5_applicability[1] == 0 && lmsd_satisfied == 0)
    print("Level 1 Criteria Not Met - Perform Level 2 or 3 as applicable")
elseif (part5_applicability[1] == 1 && lmsd_satisfied == 0)
    print("Level 1 Criteria Not Met - Perform Level 2 or 3 as applicable")
elseif (part5_applicability[1] == 0 && lmsd_satisfied == 1)
    print("Level 1 Criteria Not Met - Perform Level 2 or 3 as applicable")
end


text:
function gaussianKernel(arr, σₗ=0.5::Float, σₛ=0.5::Float)
	A = fill(0.0, size(arr,1), size(arr,1))

	for i in 1:size(A,1)
		for j in 1:size(A,1)
			if i == j
				A[i,j] = 0
			else
				# Distance in LAB space
				aₗ= exp.(-((arr[i,1] - arr[j,1])^2)/(σₗ^2))

				# Distance in pixel space
				aₛ = exp.(-euclidean(arr[i,4:5], arr[j,4:5])/(σₛ^2))

				# Combined distance metric
				A[i,j] = aₗ .* aₛ

			end
		end
	end

	return A
end

max_stars_repo_name: neuro-myoung/CellSegmenter

max_issues_repo_name: neuro-myoung/CellSegmenter

max_forks_repo_name: neuro-myoung/CellSegmenter


text:
"""
Author: Mohammad Dehghani Ashkezari <mdehghan@uw.edu>

Date: 2020-05-07

Function: HTTP requests to the CMAP RESTful API.
"""

using HTTP
using CSV
using DataFrames
using JSON

mutable struct API
    token :: String
    domain :: String
    tokenPrefix :: String

    function API(
                token::String="",
                tokenPrefix::String="Api-Key ",
                domain::String="https://simonscmap.com"
                )

        if isempty(token)
            token = get_api_key();
        end

        @assert(length(token)>0, "API key not found!")

        api = new();
        api.token = token;
        api.tokenPrefix = tokenPrefix;
        api.domain = domain;
        api
    end
end

function api_key_fname()
    return "api_key.csv"
end

"""
    set_api_key(key::String)

Saves your API key on local disk.

Examples
≡≡≡≡≡≡≡≡≡≡
set_api_key("Your API Key")

"""
function set_api_key(key::String)
    CSV.write(api_key_fname(), DataFrame(apiKey = key))
end

"""
    get_api_key()

Returns your api key, if exists.

Examples
≡≡≡≡≡≡≡≡≡≡
get_api_key()

"""
function get_api_key()
    keyFile = api_key_fname()
    if ~isfile(keyFile)
        error("API Key file not found. The following command will register the API key:\n\nset_api_key(key::String)\n\n")
    end
    return DataFrame(CSV.File(keyFile))[1, 1]
end

function atomic_request(api::API, route, payload)
    try
        baseURL = api.domain;
        tokenPrefix = api.tokenPrefix;
        token = tokenPrefix * api.token;
        queryString = "";
        if ~isempty(payload)
            queryString = HTTP.escapeuri(payload);
        end
        url = baseURL * route * queryString;
        headers = (("Authorization", token),
                   ("content-type","application/json")
                  )
        response = HTTP.request("GET", url, headers);
        status = response.status;
        respStr = String(response.body);
        df = CSV.read(IOBuffer(chomp(respStr)))
        return status, df
    catch e
        println("Error in atomic_request : $e")
        return
    end
end

"""
    query(api::API, queryStatement::String)

Takes a custom SQL query statement and returns the results in form of a dataframe.

Examples
≡≡≡≡≡≡≡≡≡≡
api = API()
status, response = query(api, "SELECT * FROM tblSensors")
println(response)

"""
function query(api::API, queryStatement::String)
    payload = Dict("query" => queryStatement);
    route = "/api/data/query?";
    return atomic_request(api, route, payload)
end


text:
function complete_lower_bound(; y=y, idx=idx, S=S, μα=μα, Σα=Σα, v²ₐ=v²ₐ, μβ=μβ, Σβ=Σβ, v²ᵦ=v²ᵦ, μf=μf, Σf=Σf, K=K)

    #
    # ∫ q(α) q(β) q(f) log𝓝(y| Af + b, S) + log𝓝(f|0,K) + log𝓝(a|0,v²ₐI) + log𝓝(b|0,v²ᵦ) dα dβ df
    #                           (1)               (2)             (3)             (4)

    Elogl = exact_expectation_inside_log_of_lower_bound(; y=y, idx=idx, S=S, μα=μα, Σα=Σα, v²ₐ=v²ₐ, μβ=μβ, Σβ=Σβ, v²ᵦ=v²ᵦ, μf=μf, Σf=Σf, K=K)

    #
    # H[q(α)]  +  H[q(β)]  +  H[q(f)]
    #

    gaussianentropy(Σ) = 0.5*logdet(2.0 * π * ℯ * Σ)

    return Elogl + gaussianentropy(Σα) + gaussianentropy(Σβ) + gaussianentropy(Σf)

end


text:
JS =  "// activate extensions only after Notebook is initialized\nrequire([\"base/js/events\"], function (events) {\n\$([IPython.events]).on(\"app_initialized.NotebookApp\", function () {\n    /* load your extension here */\n    IPython.load_extensions('hide_input');\n    });\n});"

# Test for IPython 3 (required)
const ipvers = try
    convert(VersionNumber, chomp(readall(`ipython --version`)))
catch e1
        error("IPython 3 is required for Revealables.\n", "$e1")
end

install = "install"

if ipvers < v"3.0"
    warn("IPython 3 is required for Revealables.\n" *
        "You have IPython $ipvers.\n" *
        "This installer will download hide_input.js to your nbextensions folder,\n" *
        "but you will need to manually check your custom.js file to\n" *
        "ensure that hide_input is installed correctly.\n" *
        "Instructions at https://github.com/ipython-contrib/\n" *
        "      IPython-notebook-extensions/wiki#adding-the-extension-to-customjs\n" *
        "Continue installation? yes/[no]")
    while install != "yes"
        install = chomp(readline(STDIN))
        if install == "no" || install == ""
            error("Installation cancelled.")
        end
        warn("Continue installation? yes/[no]")
    end
end

# Copy hide_input.js to correct location
println("Please enter the name of the profile you use with IJulia.
    \nLeave blank for default")
profile = chomp(readline(STDIN))
if profile == ""
    profile = "default"
end

extdir = chomp(readall(`ipython locate nbextensions`))
profiledir = chomp(readall(`ipython locate profile $profile`))

download("https://raw.githubusercontent.com/ipython-contrib/IPython-notebook-extensions/master/usability/hide_input.js", Pkg.dir(extdir,"nbextensions","hide_input.js"))

if ipvers < v"3.0"
    run(`ipython --profile $profile hide_input_setup.txt`)
else
    #open custom.js file,
    jsFilename = Pkg.dir(profiledir,"static","custom","custom.js")
    touch(jsFilename)
    f = open(jsFilename)
    lines = readlines(f)
    close(f)

    foundextensions = false
    foundhide_input = false

    tmpjsname = Pkg.dir(profiledir,"static","custom","tmp.js")
    touch(tmpjsname)
    tmpjsfile = open(tmpjsname, "w")

    for l in lines
        if foundextensions && (!foundhide_input && chomp(l) == chomp("    });\n"))
            write(tmpjsfile, "    IPython.load_extensions('hide_input');\n")
            foundhide_input = true
            println("Installed hide-input.js")
        end
        write(tmpjsfile, l)
        if chomp(l) == chomp("\$([IPython.events]).on(\"app_initialized.NotebookApp\", function () {\n")
            foundextensions = true
        elseif chomp(l) == chomp("    IPython.load_extensions('hide_input');\n")
            foundhide_input = true
            println("hide-input.js already installed; nothing to do")
        end
    end
    if !foundextensions
        write(tmpjsfile, JS)
        println("Installed hide-input.js")
    end
    close(tmpjsfile)
    mv(tmpjsname, jsFilename)
end


text:
using DataStructures
using NCDatasets
using Dates

function SaveOutput(fname)
  Dataset(fname,"c",attrib = OrderedDict("title" => "test")) do ds
    defDim(ds,"Alon",1)
    defDim(ds,"Alat",1)
    defDim(ds,"Alev",vleva)
    defDim(ds,"Olon",1)
    defDim(ds,"Olat",1)
    defDim(ds,"Olev",vlevo)
    defDim(ds,"time",1)
    defVar(ds,"time",Dates.datetime2epochms(mtime)/1000/3600/24,("time",),
      attrib = OrderedDict("units" => "days since 0000-01-01 00:00:00" ))
    defVar(ds,"Alon",lon,("Alon",),
      attrib = OrderedDict("units" => "degree"))
    defVar(ds,"Alat",lon,("Alat",),
      attrib = OrderedDict("units" => "degree"))
    defVar(ds,"Alev",ZAC,("Alev",),
      attrib = OrderedDict("units" => "m"))
    defVar(ds,"Olon",lon,("Olon",),
      attrib = OrderedDict("units" => "degree"))
    defVar(ds,"Olat",lon,("Olat",),
      attrib = OrderedDict("units" => "degree"))
    defVar(ds,"Olev",ZOC,("Olev",),
      attrib = OrderedDict("units" => "m"))
    defVar(ds,"THETA",ΘA,("Alev",),
      attrib = OrderedDict("units" => "degree Kelvin"))
    defVar(ds,"TA",TA,("Alev",),
      attrib = OrderedDict("units" => "degree Celsius"))
    defVar(ds,"qA",qA,("Alev",),
      attrib = OrderedDict("units" => "kg/kg"))
    defVar(ds,"UA",UA,("Alev",),
      attrib = OrderedDict("units" => "m/s"))
    defVar(ds,"rhoA",ρA,("Alev",),
      attrib = OrderedDict("units" => "kg/m^3"))
    defVar(ds,"TO",TO,("Olev",),
      attrib = OrderedDict("units" => "degree Celsius"))
    defVar(ds,"UO",UO,("Olev",),
      attrib = OrderedDict("units" => "m/s"))
    defVar(ds,"rhoO",ρO,("Olev",),
      attrib = OrderedDict("units" => "kg/kg"))
  end
end

function ReadOutput(fname)
  ds = NCDataset(fname)
  global lat=ds["Alat"][1]
  global lon=ds["Alon"][1]
  global time=ds["time"][1]
  global ZAC=ds["Alev"][:]
  global ZOC=ds["Olev"][:]
  global ΘA=ds["THETA"][:]
  global TA=ds["TA"][:]
  global qA=ds["qA"][:]
  global UA=ds["UA"][:]
  global ρA=ds["rhoA"][:]
  global TO=ds["TO"][:]
  global UO=ds["UO"][:]
  global ρO=ds["rhoO"][:]
end


text:
module Wright

export wright

include("main.jl")

end


text:
using Distributions
import Base.size, Base.length, Base.rand, Base.mean, Base.convert

doc"""
    MatrixNormal(M, U, V)
The *Matrix Normal distribution* with mean `M` and row and column covariances `U` and `V` has probability density function
$f(\mathbf{X}; \mathbf{M}, \mathbf{U}, \mathbf{V}) = \frac{\exp\left(
    -\frac{1}{2} \mathrm{tr}\left[
        \mathbf{V}^{-1} (\mathbf{X} - \mathbf{M})^\top \mathbf{U}^{-1}
        (\mathbf{X} - \mathbf{M})
    \right]
    \right)}{(2\pi)^{np/2} |\mathbf{V}|^{n/2} |\mathbf{U}|^{p/2}}$
when $\mathbf{X}$ is $n \times p$.
```julia
MatrixNormal(U, V)      # Matrix Normal with mean 0 and covariances U and V
MatrixNormal(M, U, V)      # Matrix Normal with mean M and covariances U and V
params(d)        # Get the parameters, i.e. (M, U, V)
mean(d)          # Get the mean, i.e., M
cov(d)           # Get the covariance, i.e., V ⊗ U
var(d)           # Get the marginal variances for individual entries
invcov(d)        # inverse of covariance
logdetcov(d)     # log determinant of covariance
logpdf(d, X)     # log of the pdf at point X
```
External links
* [Matrix Normal distribution on Wikipedia](http://en.wikipedia.org/wiki/Matrix_normal_distribution)
"""
immutable MatrixNormal{Mean<:AbstractMatrix, CovU<:AbstractPDMat, CovV<:AbstractPDMat} <: ContinuousMatrixDistribution
    M::Mean
    U::CovU
    V::CovV
end

### Constructors
typealias CovLike Union{AbstractMatrix, AbstractVector, AbstractPDMat}
normalize(σ::AbstractVector) = PDiagMat(abs2(σ))
normalize(Σ::AbstractMatrix) = PDMat(Σ)
normalize(Σ::AbstractPDMat) = Σ
normalize(d::Int, σ::Real) = PDScalMat(d, abs2(σ))

function MatrixNormal{Mean <: AbstractMatrix, CovU<:AbstractPDMat, CovV<:AbstractPDMat}(M::Mean, U::CovU, V::CovV)
    size(M) == (dim(U), dim(V)) || throw(DimensionMismatch("The dimensions of M, U, and V are inconsistent."))
    MatrixNormal{Mean, CovU, CovV}(M, U, V)
end
function MatrixNormal{Mean<:AbstractMatrix, CovU<:CovLike, CovV<:CovLike}(M::Mean, U::CovU, V::CovV)
    MatrixNormal(M, normalize(U), normalize(V))
end

MatrixNormal{CovU<:AbstractPDMat, CovV<:AbstractPDMat}(U::CovU, V::CovV) = MatrixNormal(zeros(dim(U), dim(V)), U, V)
MatrixNormal{CovU<:CovLike, CovV<:CovLike}(U::CovU, V::CovV) = MatrixNormal(normalize(U), normalize(V))

### Conversions
convert(::Type{MvNormal}, d::MatrixNormal) = MvNormal(vec(mean(d)), cov(d))

### interface
size(d::MatrixNormal) = (dim(d.U), dim(d.V))
params(d::MatrixNormal) = (d.M, d.U, d.V)
mean(d::MatrixNormal) = d.M
cov(d::MatrixNormal) = kron(full(d.V), full(d.U))
var(d::MatrixNormal) = diag(d.U) * diag(d.V)'
invcov(d::MatrixNormal) = kron(full(inv(d.V)), full(inv(d.U)))
logdetcov(d::MatrixNormal) = ((n, p) = size(d); n * logdet(d.V) + p * logdet(d.U))

rand(d::MatrixNormal) = reshape(rand(convert(MvNormal, d)), size(d))
entropy(d::MatrixNormal) = entropy(convert(MvNormal, d))

function logpdf(d::MatrixNormal, X::AbstractMatrix)
    size(X) == size(d) || error("Input matrix x has wrong size.")
    m = X - d.M
    n, p = size(d)
    invU = full(inv(d.U))
    invV = full(inv(d.V))
    lpdf = -0.5 * trace(invV * m' * invU * m)
    lpdf -= (n * p/2) * Float64(Distributions.log2π)
    lpdf -= 0.5 * logdetcov(d)
    lpdf
end


text:
te = [-0.112639 0.124715 0.124715 0.268717 0.124715 0.268717 0.268717 0.046154]
st = SymmetricTensor(reshape(te, (2,2,2)))
Random.seed!(43)
c = cumulants(randn(1000, 20), 4)

@testset "axiliary functions" begin
    @test rep((1,2,3)) == 6
    @test rep((1,2,2)) == 3
    @test rep((1,1,1)) == 1
end

@testset "norm" begin
    @test norm(st) ≈ 0.5273572868359742
    @test norm(st, 1) ≈ 1.339089
    @test norm(st, 2.5) ≈ norm(te, 2.5)
    n = norm(c[2])
    @test n ≈ norm(Array(c[2]))
    @test cnorms(c) ≈ [norm(c[3])/(n^(3/2)), norm(c[4])/(n^2)]
end


text:
#**************************************************************************************
# build.jl
# =============== part of the GeoEfficiency.jl package.
#
# this file excuted at instlation time.
#
#**************************************************************************************

# creating `.GeoEfficiency` folder at the home directory.

@info("Creating '.GeoEfficiency' data folder at '$(homedir())'.....")

@info("Creating '.GeoEfficiency' folder at '$(homedir())'.....")
try
	cp(joinpath(dirname(@__FILE__),".batch"), joinpath(homedir(),".GeoEfficiency"))
catch err
	if isdir(joinpath(homedir(),".GeoEfficiency"))
		@warn("`.GeoEfficiency` data folder already exist in your home directory.")
	else
		@error("""'.GeoEfficiency' data folder: could not be created,
		\n  the package may be unable to work in patch mode.
		\n  Try 'import Pkg; Pkg.build("GeoEfficiency")' to rebuild.""")
	end #if
end #try

max_stars_repo_name: DrKrar/GeoEfficiency

max_issues_repo_name: DrKrar/GeoEfficiency

max_forks_repo_name: DrKrar/GeoEfficiency


text:
# Wisdom & Hernandez version of Kepler solver, but with quartic
# convergence.

function calc_ds_opt(y,yp,ypp,yppp)
# Computes quartic Newton's update to equation y=0 using first through 3rd derivatives.
# Uses techniques outlined in Murray & Dermott for Kepler solver.
# Rearrange to reduce number of divisions:
num = y*yp
den1 = yp*yp-y*ypp*.5
den12 = den1*den1
den2 = yp*den12-num*.5*(ypp*den1-third*num*yppp)
return -y*den12/den2
end

function kep_elliptic!(x0::Array{Float64,1},v0::Array{Float64,1},r0::Float64,dr0dt::Float64,k::Float64,h::Float64,beta0::Float64,s0::Float64,state::Array{Float64,1})
# Solves equation (35) from Wisdom & Hernandez for the elliptic case.

r0inv = inv(r0)
beta0inv = inv(beta0)
# Now, solve for s in elliptical Kepler case:
if beta0 > 1e-15
# Initial guess (if s0 = 0):
  if s0 == 0.0
    s = h*r0inv
  else
    s = copy(s0)
  end
  s0 = copy(s)
  sqb = sqrt(beta0)
  y = 0.0; yp = 1.0
  iter = 0
  ds = Inf
  fac1 = k-r0*beta0
  fac2 = r0*dr0dt
  while iter == 0 || (abs(ds) > 1e-8 && iter < 10)
    xx = sqb*s
    sx = sqb*sin(xx)
    cx = cos(xx)
# Third derivative:
    yppp = fac1*cx - fac2*sx
# Take derivative:
    yp = (-yppp+ k)*beta0inv
# Second derivative:
    ypp = fac1*beta0inv*sx + fac2*cx
    y  = (-ypp + fac2 +k*s)*beta0inv - h  # eqn 35
# Now, compute fourth-order estimate:
    ds = calc_ds_opt(y,yp,ypp,yppp)
    s += ds
    iter +=1
  end
#  if iter > 2
#    println(iter," ",s," ",s/s0-1," ds: ",ds)
#  end
# Since we updated s, need to recompute:
  xx = 0.5*sqb*s; sx = sin(xx) ; cx = cos(xx)
# Now, compute final values:
  g1bs = 2.*sx*cx/sqb
  g2bs = 2.*sx^2*beta0inv
  f = 1.0 - k*r0inv*g2bs # eqn (25)
  g = r0*g1bs + fac2*g2bs # eqn (27)
  for j=1:3
# Position is components 2-4 of state:
    state[1+j] = x0[j]*f+v0[j]*g
  end
  r = sqrt(state[2]*state[2]+state[3]*state[3]+state[4]*state[4])
  rinv = inv(r)
  dfdt = -k*g1bs*rinv*r0inv
  dgdt = r0*(1.0-beta0*g2bs+dr0dt*g1bs)*rinv
  for j=1:3
# Velocity is components 5-7 of state:
    state[4+j] = x0[j]*dfdt+v0[j]*dgdt
  end
else
  println("Not elliptic ",beta0," x0 ",x0)
end
# recompute beta:
state[8]= r
state[9] = (state[2]*state[5]+state[3]*state[6]+state[4]*state[7])*rinv
# beta is element 10 of state:
state[10] = 2.0*k*rinv-(state[5]*state[5]+state[6]*state[6]+state[7]*state[7])
# s is element 11 of state:
state[11] = s
# ds is element 12 of state:
state[12] = ds
return iter
end

function kep_hyperbolic!(x0::Array{Float64,1},v0::Array{Float64,1},r0::Float64,dr0dt::Float64,k::Float64,h::Float64,beta0::Float64,s0::Float64,state::Array{Float64,1})
# Solves equation (35) from Wisdom & Hernandez for the hyperbolic case.

r0inv = inv(r0)
beta0inv = inv(beta0)
# Now, solve for s in hyperbolic Kepler case:
if beta0 < -1e-15
# Initial guess (if s0 = 0):
  if s0 == 0.0
    s = h*r0inv
  else
    s = copy(s0)
  end
  s0 = copy(s)
  sqb = sqrt(-beta0)
  y = 0.0; yp = 1.0
  iter = 0
  ds = Inf
  fac1 = k-r0*beta0
  fac2 = r0*dr0dt
  while iter == 0 || (abs(ds) > 1e-8 && iter < 10)
    xx = sqb*s; cx = cosh(xx); sx = sqb*(exp(xx)-cx)
# Third derivative:
    yppp = fac1*cx + fac2*sx
# Take derivative:
    yp = (-yppp+ k)*beta0inv
# Second derivative:
    ypp = -fac1*beta0inv*sx  + fac2*cx
    y  = (-ypp +fac2 +k*s)*beta0inv - h  # eqn 35
# Now, compute fourth-order estimate:
    ds = calc_ds_opt(y,yp,ypp,yppp)
    s += ds
    iter +=1
  end
#  if iter > 2
#    #println("iter: ",iter," ds/s: ",ds/s0)
#    println(iter," ",s," ",s/s0-1," ds: ",ds)
#  end
  xx = 0.5*sqb*s; cx = cosh(xx); sx = exp(xx)-cx
# Now, compute final values:
  g1bs = 2.0*sx*cx/sqb
  g2bs = -2.0*sx^2*beta0inv
  f = 1.0 - k*r0inv*g2bs # eqn (25)
  g = r0*g1bs + fac2*g2bs # eqn (27)
  for j=1:3
    state[1+j] = x0[j]*f+v0[j]*g
  end
  # r = norm(x)
  r = sqrt(state[2]*state[2]+state[3]*state[3]+state[4]*state[4])
  rinv = inv(r)
  dfdt = -k*g1bs*rinv*r0inv
  dgdt = r0*(1.0-beta0*g2bs+dr0dt*g1bs)*rinv
  for j=1:3
# Velocity is components 5-7 of state:
    state[4+j] = x0[j]*dfdt+v0[j]*dgdt
  end
else
  println("Not hyperbolic",beta0," x0 ",x0)
end
# recompute beta:
state[8]= r
state[9] = (state[2]*state[5]+state[3]*state[6]+state[4]*state[7])*rinv
# beta is element 10 of state:
state[10] = 2.0*k*rinv-(state[5]*state[5]+state[6]*state[6]+state[7]*state[7])
# s is element 11 of state:
state[11] = s
# ds is element 12 of state:
state[12] = ds
return iter
end


text:
struct FastAIDataset
    name
    subfolder
    extension
    description
    checksum
    datadepname
    subpath
    size
end

const ROOT_URL = "https://s3.amazonaws.com/fast-ai-"

function FastAIDataset(
        name, subfolder, checksum="";
        extension="tgz",
        description="",
        datadepname=name,
        subpath=name,
        size="???")
    return FastAIDataset(name, subfolder, extension, description, checksum, datadepname, subpath, size)
end

const DESCRIPTIONS = Dict(
    "imagenette" => "A subset of 10 easily classified classes from Imagenet: tench, English springer, cassette player, chain saw, church, French horn, garbage truck, gas pump, golf ball, parachute",
    "imagewoof" => "A subset of 10 harder to classify classes from Imagenet (all dog breeds): Australian terrier, Border terrier, Samoyed, beagle, Shih-Tzu, English foxhound, Rhodesian ridgeback, dingo, golden retriever, Old English sheepdog",
    "food-101" => "101 food categories, with 101,000 images; 250 test images and 750 training images per class. The training images were not cleaned. All images were rescaled to have a maximum side length of 512 pixels.",

)

const DATASETCONFIGS = [
    # imageclas
    FastAIDataset("CUB_200_2011", "imageclas", "0c685df5597a8b24909f6a7c9db6d11e008733779a671760afef78feb49bf081", size="1GiB"),
    FastAIDataset("bedroom", "imageclas", "7c95250ccb177c582f602c08f239c71f7a70512729d2e078925261cf5e349f5d", size="4.25GiB"),
    FastAIDataset("caltech_101", "imageclas", "af6ece2f339791ca20f855943d8b55dd60892c0a25105fcd631ee3d6430f9926", size="126MiB", subpath="101_ObjectCategories"),
    FastAIDataset("cifar10", "imageclas", "637c5814e11aefcb6ee76d5f59c67ddc8de7f5b5077502a195b0833d1e3e4441"),
    FastAIDataset("cifar100", "imageclas", "085ac613ceb0b3659c8072143ae553d5dd146b3c4206c3672a56ed02d0e77d28"),
    FastAIDataset("food-101", "imageclas", "abc3d6b03a9886fdea6d2a124cf88e22a99dfdb03085b2478be97de3f8e4679f", size="5.3GB", description=DESCRIPTIONS["food-101"]),
    FastAIDataset("imagenette-160", "imageclas", "1bd650bc16884ca88e4f0f537ed8569b1f8d7ae865d37eba8ecdd87d9cd9dcfa", size="1.45GiB", description=DESCRIPTIONS["imagenette"]),
    FastAIDataset("imagenette-320", "imageclas", description=DESCRIPTIONS["imagenette"]),
    FastAIDataset("imagenette", "imageclas", description=DESCRIPTIONS["imagenette"]),
    FastAIDataset("imagenette2-160", "imageclas", "64d0c4859f35a461889e0147755a999a48b49bf38a7e0f9bd27003f10db02fe5", description=DESCRIPTIONS["imagenette"]),
    FastAIDataset("imagenette2-320", "imageclas", "569b4497c98db6dd29f335d1f109cf315fe127053cedf69010d047f0188e158c", description=DESCRIPTIONS["imagenette"]),
    FastAIDataset("imagenette2", "imageclas", "6cbfac238434d89fe99e651496f0812ebc7a10fa62bd42d6874042bf01de4efd", description=DESCRIPTIONS["imagenette"]),
    FastAIDataset("imagewang-160", "imageclas", "a0d360f9d8159055b3bf2b8926a51d19b2f1ff98a1eef6034e4b891c59ca3f1a", size="182MiB"),
    FastAIDataset("imagewang-320", "imageclas", "fd53301c335aa46f0f4add68dd471cd0b8b66412382cc36f5f510d0a03fb4d9d", size="639MiB"),
    FastAIDataset("imagewang", "imageclas"),
    FastAIDataset("imagewoof-160", "imageclas", "a0d360f9d8159055b3bf2b8926a51d19b2f1ff98a1eef6034e4b891c59ca3f1a", description=DESCRIPTIONS["imagewoof"]),
    FastAIDataset("imagewoof-320", "imageclas", description=DESCRIPTIONS["imagewoof"]),
    FastAIDataset("imagewoof", "imageclas", description=DESCRIPTIONS["imagewoof"]),
    FastAIDataset("imagewoof2-160", "imageclas", "663c22f69c2802d85e2a67103c017e047096702ffddf9149a14011b7002539bf", description=DESCRIPTIONS["imagewoof"]),
    FastAIDataset("imagewoof2-320", "imageclas", "7db6120fdb9ae079e26346f89e7b00d7f184f8137791609b97fd0405d3f92305", description=DESCRIPTIONS["imagewoof"], size="313MB"),
    FastAIDataset("imagewoof2", "imageclas", description=DESCRIPTIONS["imagewoof"]),
    FastAIDataset("mnist_png", "imageclas", "9e18edaa3a08b065d8f80a019ca04329e6d9b3e391363414a9bd1ada30563672"),
    FastAIDataset("mnist_var_size_tiny", "imageclas", "8a0f6ca04c2d31810dc08e739c7fa9b612e236383f70dd9fc6e5a62e672e2283"),
    FastAIDataset("oxford-102-flowers", "imageclas"),
    FastAIDataset("oxford-iiit-pet", "imageclas"),
    FastAIDataset("stanford-cars", "imageclas"),

    # nlp
    FastAIDataset("ag_news_csv", "nlp", "9a8c300eabb45750237fcc669f61cb8a3448f3ef6f6098e1ce340e444f6872be", size="11MB"),
    FastAIDataset("amazon_review_full_csv", "nlp", "4af62eeee139d0142e0747340b68646d23483d9475c33ea0641ee9175b423443", size="600MB"),
    FastAIDataset("amazon_review_polarity_csv", "nlp", "d2a3ee7a214497a5d1b8eaed7c8d7ba2737de00ada3b0ec46243983efa100361", size="600MB"),
    FastAIDataset("dbpedia_csv", "nlp", "42db5221ddedddb673a4cabcc5f3a7d869714c878bcfe4ba94b29d14aa38e417", size="65MB"),
    FastAIDataset("giga-fren", "nlp", "11c97af99471fe641f210d8b86ccccf3b298b9199853987ee53892d709d7ca6b", size="2.4GB"),
    FastAIDataset("imdb", "nlp", "d501018afa17aee9fa1ebe8ac29859a5609980e13dc6e611aa21567cc357351f", size="140MB"),
    FastAIDataset("sogou_news_csv", "nlp", "6b77fc935561d339b82aa552d7e31ea59eff492a494920579b3ce70604efb5c2", size="360MB"),
    FastAIDataset("wikitext-103", "nlp", "27b89e94d98a9f9db74588a2e75b04378ee21569ce55d329d3e73e27d0952551", size="181MB"),
    FastAIDataset("wikitext-2", "nlp", "4e39df0e84453ae2f3d34333de2a9d8e57560a7a6e621f13e11dc21241320074", size="4MB"),
    FastAIDataset("yahoo_answers_csv", "nlp", "2d4277855faf8b35259009425fa8f7fe1888b5644b47165508942d000f4c96ae", size="305MB"),
    FastAIDataset("yelp_review_full_csv", "nlp", "56006b0a17a370f1e366504b1f2c3e3754e4a3dda17d3e718a885c552869a559", size="187MB"),
    FastAIDataset("yelp_review_polarity_csv", "nlp", "528f22e286cad085948acbc3bea7e58188416546b0e364d0ae4ca0ce666abe35", size="158MB"),

    # imagelocal
    FastAIDataset("biwi_head_pose", "imagelocal", "9cfefd53ed85f824c5908bc6eb21fc719583eec57a7df1d8141d3156645693cf", size="430MiB"),
    FastAIDataset("camvid", "imagelocal", "11db05fc3ee727fb17de7499380b20258a41beeb1002a2aee2c2244a472a4a45", size="571MB"),
    FastAIDataset("pascal-voc", "imagelocal", "10fc13a659da20fdd8302dd394d88ca7e4e60e69fd8a5212c3e3357964a58215", size="4.3GB"),
    FastAIDataset("pascal_2007", "imagelocal"),
    FastAIDataset("pascal_2012", "imagelocal"),
    FastAIDataset("siim_small", "imagelocal"),
    FastAIDataset("skin-lesion", "imagelocal"),
    FastAIDataset("tcga-small", "imagelocal"),

    # sample
    FastAIDataset("adult_sample", "sample", "47ecd1848abc976643ee82d8788b712e3006d629bbc7554efa1077a91579e99e", size="3.8MB"),
    FastAIDataset("biwi_sample", "sample"),
    FastAIDataset("camvid_tiny", "sample", "cd42a9bdd8ad3e0ce87179749beae05b4beb1ae6ab665841180b1d8022fc230b"),
    FastAIDataset("dogscats", "sample", "b79c0a5e4aa9ba7a0b83abbf61908c61e15bed0e5b236e86a0c4a080c8f70d7c", size="800MiB"),
    FastAIDataset("human_numbers", "sample"),
    FastAIDataset("imdb_sample", "sample", "8e776d995296136b3f9a3cf001796d886cb0b60e86877ce71c7abbdc3c247341", size="4KB"),
    FastAIDataset("mnist_sample", "sample", "b373a14f282298aeba0f7dd56b7cdb6c2401063d4f118c39c54982907760bd38", size="3MB"),
    FastAIDataset("mnist_tiny", "sample", "0d1fedf86243931aa3fc065d2cf4ffab339a972958d8594ae993ee32bd8e15b9", size="300KB"),
    FastAIDataset("movie_lens_sample", "sample"),
    FastAIDataset("planet_sample", "sample", "f2509212bb2dcdc147423b164564f2e63cae1d1db0b504166e5b92cfbcbb3b4c", size="14.8MB"),
    FastAIDataset("planet_tiny", "sample", "41a5fdd82db1c9fb2cff17e1a1270102414a25a34b21b770f953d28483961edb", size="1MB"),

    # coco
    FastAIDataset("coco_sample", "coco", "56960c0ac09ff35cd8588823d37e1ed0954cb88b8bfbd214a7763e72f982911c", size="3GB"),
    FastAIDataset("train2017", "coco", datadepname="coco-train2017", extension="zip"),
    FastAIDataset("val2017", "coco", datadepname="coco-val2017", extension="zip"),
    FastAIDataset("test2017", "coco", datadepname="coco-test2017", extension="zip"),
    FastAIDataset("unlabeled2017", "coco", datadepname="coco-unlabeled2017", extension="zip"),
    FastAIDataset("image_info_test2017", "coco", datadepname="coco-image_info_test2017", extension="zip"),
    FastAIDataset("image_info_unlabeled2017", "coco", datadepname="coco-image_info_unlabeled2017", extension="zip"),
    FastAIDataset("annotations_trainval2017", "coco", datadepname="coco-annotations_trainval2017", extension="zip"),
    FastAIDataset("stuff_annotations_trainval2017", "coco", datadepname="coco-stuff_annotations_trainval2017", extension="zip"),
    FastAIDataset("panoptic_annotations_trainval2017", "coco", datadepname="coco-panoptic_annotations_trainval2017", extension="zip"),
]

const DATASETS = [d.datadepname for d in DATASETCONFIGS]
const DATASETS_IMAGECLASSIFICATION = vcat(
    [d.datadepname for d in DATASETCONFIGS if d.subfolder == "imageclas"],
    ["mnist_sample", "mnist_tiny", "dogscats"],

)

function DataDeps.DataDep(d::FastAIDataset)
    return DataDep(
        "fastai-$(d.datadepname)",
        """
        "$(d.name)" from the fastai dataset repository (https://course.fast.ai/datasets)

        $(d.description)

        Download size: $(d.size)
        """,
        "$(ROOT_URL)$(d.subfolder)/$(d.name).$(d.extension)",
        d.checksum,
        post_fetch_method=DataDeps.unpack,
    )
end

function initdatadeps()
    for d in DATASETCONFIGS
        DataDeps.register(DataDep(d))
    end
end


text:
module CovariantVectors

# Re-export methods from DifferentialEquations extended here:
export init, solve, step!
import DifferentialEquations: init, solve, step!

using ..Stages: AbstractSource, AbstractStage, finish_if_not!,
    StageIterator, StageState, StagedSolver, goto!
import ..Stages: current_result, stage_index, record!

using ..LyapunovExponents: LEProblem, LESolver, dimension, is_semi_unitary,
    default_Q0, get_tangent_prob, validate_tangent_prob, time_type, ceil_if
import ..LyapunovExponents:
    phase_state, TangentRenormalizer, get_integrator
using LinearAlgebra
using Markdown

include("utils.jl")
include("clv_problem.jl")
include("core_stages.jl")
include("clv_solver.jl")
include("accessors.jl")

end

using .CovariantVectors: goto!

using .CovariantVectors:
    CLVProblem, CLVSolver, CLV, backward_dynamics!, forward_dynamics!,
    indexed_backward_dynamics!, indexed_forward_dynamics!
export
    CLVProblem, CLVSolver, CLV, backward_dynamics!, forward_dynamics!,
    indexed_backward_dynamics!, indexed_forward_dynamics!


text:
"""
Implementation of the Graphical Lasso (Glasso) for l1-penalized precision matrix estimation

Code credit: https://github.com/bnaul/GraphLasso.jl
"""

function glasso(S::Matrix{Float64}, α::Float64; tol::Float64 = 1e-5,
                    maxit::Int = 1000, penalize_diag::Bool = true, verbose::Bool = false)
    p = size(S, 1)
    adj = abs.(S) .> α
    blocks = connected_components(SimpleGraph(max.(adj, I(p))))
    if(verbose)
        print(blocks)
    end
    Θ = zeros(p, p)
    for block = blocks
        Θ[block,block] .= fit_block(S[block,block], α, tol, maxit, penalize_diag)
    end
    return Θ
end

function fit_block(S::Matrix{Float64}, α::Float64, tol::Float64,
                   maxit::Int, penalize_diag::Bool)
    p = size(S, 1)
    W = copy(S)
    if penalize_diag
        W .+= α * I(p)
    end
    if size(S, 1) == 1
        return inv(W)
    end
    Θ = zeros(p, p)
    i = 0
    β = zeros(p-1, p)
    while i < maxit
        i += 1
        W_old = copy(W)
        for j = 1:p
            inds = collect(1:p)
            splice!(inds, j)
            W11 = W[inds, inds]
            sqrtW11 = sqrt(W11)
            # β[:,j] = fit(LassoPath, sqrtW11, sqrtW11 \ S[inds,j], λ=[α/(p-1)],
            #              standardize=false, intercept=false).coefs
            β[:,j] = glmnet(sqrtW11, sqrtW11 \ S[inds, j], lambda=[α / (p-1)],
                         standardize = false, intercept = false).betas
            W[inds, j] .= W11 * β[:, j]
            W[j, inds] .= W[inds, j]
        end
        if norm(W - W_old) < tol
            break
        end
    end
    if i == maxit
        @warn "Maximum number of iterations reached, graphical lasso failed to converge"
    end
    for j in 1:p
        inds = collect(1:p)
        splice!(inds, j)
        Θ[j, j] = 1 / (W[j, j] - dot(W[inds, j], β[:, j]))
        Θ[inds, j] .= -Θ[j, j] * β[:, j]
        Θ[j, inds] .= Θ[inds, j]
    end
    return Θ
end

# # Testing
# S_iris = cov(Array(dataset("datasets", "iris")[:,1:4]))

# α = 1.0

# W_R = [1.68569351 0.000000 0.2743155 0.01969989
#        0.00000000 1.189979 0.0000000 0.00000000
#        0.27431545 0.000000 4.1162779 0.29560940
#        0.01969989 0.000000 0.2956094 1.58100626]

# Θ_R = [ 0.59973155 0.0000000 -0.03996709  0.00000000
#         0.00000000 0.8403507  0.00000000  0.00000000
#        -0.03996708 0.0000000  0.24890787 -0.04604166
#         0.00000000 0.0000000 -0.04604166  0.64111722]

# W_iris, Θ_iris = graphlasso(S_iris, α; tol=1e-6, penalize_diag=true)

# @test isapprox(W_R, W_iris, atol=1e-5)
# @test isapprox(Θ_R, Θ_iris, atol=1e-5)


text:
"""
Abstract type for the Baseline, each subtype should define the method
genBaseline(p::Baseline,,Nlon,Nlat,Ntime )
"""
abstract Baseline
export Baseline

"This function generates the Baseline 3D Array of size Nlon, Nlat, Ntime."
genBaseline(m::Baseline,Nlon,Nlat,Ntime)=error("Method genBaseline not implemented for $(typeof(m))")
export genBaseline

"""
Most trivial baseline case, constant mean `m` without seasonal cycle
"""
type ConstantBaseline <: Baseline
    m::Float64
end
export ConstantBaseline
genBaseline(c::ConstantBaseline,Ntime,Nlat,Nlon)=fill(c.m,Ntime,Nlat,Nlon)

"""
very simple periodic BaseLine returning a sine oscillation in the time domain,
where `nOsc` is the number of oscillations and `amp` is the amplitude
"""
type SineBaseline <: Baseline
	nOsc::Float64 # Number of oscillations
	amp::Float64
end
export SineBaseline
genBaseline(m::SineBaseline,Ntime,Nlat,Nlon)=[m.amp*sin(k*m.nOsc/Ntime*2*pi) for k=1:Ntime,j=1:Nlat,i=1:Nlon]

"""
A trend baseline where the mean is 0 and there is a trend in Lon, Lat and Time dimension. The range of the trend will be
`[-tlon,tlon], [-tlat,tlat], [-ttime,ttime]`
"""
type TrendBaseline <: Baseline
	tlon::Float64
	tlat::Float64
	ttime::Float64
end
export TrendBaseline
genBaseline(m::TrendBaseline,Ntime,Nlat,Nlon)=[(i-1)/(Nlon-1)*2*m.tlon-m.tlon + (j-1)/(Nlat-1)*2*m.tlat-m.tlat + (k-1)/(Ntime-1)*2*m.ttime-m.ttime for k=1:Ntime,j=1:Nlat,i=1:Nlon]

max_stars_repo_name: CAB-LAB/DataFarm

max_issues_repo_name: esa-esdl/DataFarm

max_forks_repo_name: esa-esdl/DataFarm


text:
function palindrome(s)
    len = length(s)
    for i = 1:(len/2)
        if(s[len-i+1]!=s[i])
            return false
        end
    end
    return true
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# # m2.1t.jl

# Load Julia packages (libraries) needed

using TuringModels

# ### snippet 2.8t

# Define the data

k = 6; n = 9;

# Define the model

@model globe_toss(n, k) = begin
  theta ~ Beta(1, 1) # prior
  k ~ Binomial(n, theta) # model
  return k, theta
end;

# Use Turing mcmc

chns = sample(globe_toss(n, k), NUTS(0.65), 1000)

# Look at the proper draws (in corrected chn2)

chns |> display

# Show the hpd region

hpd(chns, alpha=0.055) |> display

# End of `02/m2.1t.jl`


text:
using LightGraphicalModels
using Test

@testset "LightGraphicalModels.jl" begin
    # Write your tests here.
end


text:
# B = μ (I-σI+1/3Χ) H
#FT^-1{1/}

using Plots
using LinearAlgebra
pyplot()

Dipole(x,y) = begin
    r = sqrt.(x.^2 .+ y.^2)
    θ = atan.(x,y) #angle with respect to the main magnetic field
    r < 1e-16 ? 0 : 1/(4*π)*(1 .+ 3*cos.(2*θ))./r.^3
end

B0 = 1 #From scanner object

"""Schenck JF. The role of magnetic susceptibility in magnetic resonance
imaging: MRI magnetic compatibility of the first and second kinds.
Med Phys. 1996;23(6):815-850. doi:10.1118/1.597854"""
Sphere(x,R,Δχ) = begin
    if norm(x) <= R #Inside sphere
        2Δχ/3*B0 #2/3
    else #Outside
        Δχ/3*B0*R^3*(2x[3]^2-x[1]^2-x[2]^2)/norm(x)^5
    end
end

"""Schenck JF. The role of magnetic susceptibility in magnetic resonance
imaging: MRI magnetic compatibility of the first and second kinds.
Med Phys. 1996;23(6):815-850. doi:10.1118/1.597854"""
CylinderZ(x,R,Δχ) = begin
    if norm(x) <= R #Inside sphere
        Δχ*B0 #1
    else #Outside
        0
    end
end

"""Schenck JF. The role of magnetic susceptibility in magnetic resonance
imaging: MRI magnetic compatibility of the first and second kinds.
Med Phys. 1996;23(6):81 5-850. doi:10.1118/1.597854"""
CylinderY(x,R,Δχ) = begin
    if norm(x) <= R #Inside sphere
        Δχ/2*B0 #1/2
    else #Outside
        Δχ/2*B0*R^2*(x[3]^2-x[1]^2)/(x[3]^2+x[1]^2)^2
    end
end

# contourf([CylinderY([x,0,z],.1,.1) for z=-1:.001:1, x=-1:.001:1],aspect_ratio=:equal)
# contourf([Sphere([x,0,z],.1,.1) for z=-1:.001:1, x=-1:.001:1],aspect_ratio=:equal)


text:
module NetworkDiscovery

using LightGraphs
import LightGraphs: add_edge!, neighbors, degree
import POMDPs: POMDP
import POMDPs: reward
import POMDPs: create_state, create_transition_distribution, transition
import POMDPs: create_observation, create_observation_distribution, observation
import POMDPs: Belief, belief, create_belief
import POMDPs: Policy, action, create_action
import POMDPs: AbstractDistribution
import POMDPs: initial_belief, isterminal
import POMDPs: discount
import POMCP: sparse_actions

export CommunityNetwork, CommunityAffiliationPOMDP
export ProbeNode, GuessAffiliation
export Neighborhood

export DiscoveryHeuristic, RevealedGraph
export ProbeHighestDegree, GuessBasedOnNeighbors

export generate_network, generate_problem

# state
type CommunityNetwork <: Belief
    graph::Graph
    n_comms::Int
    labels::Vector{Int}
    budget::Int
    CommunityNetwork() = new()
end

# pomdp
type CommunityAffiliationPOMDP <: POMDP
    budget::Int
    node::Int
    reward::Float64
    true_network::CommunityNetwork
    known_nodes::Set{Int}
    known_edges::Set{Edge}
    known_labels::Dict{Int,Int}
    p_intra::Float64
    p_inter::Float64
end
Base.length(p::CommunityAffiliationPOMDP) = length(p.true_network.graph.vertices)

isterminal(::CommunityAffiliationPOMDP, s::CommunityNetwork) = s.budget < 0

# actions
create_action(pomdp::CommunityAffiliationPOMDP) = nothing
type ProbeNode
    node::Int
    ProbeNode()=new()
    ProbeNode(i)=new(i)
end
==(u::ProbeNode, v::ProbeNode) = u.node==v.node
hash(a::ProbeNode) = hash(a.node)
type GuessAffiliation
    community::Int
    GuessAffiliation()=new()
    GuessAffiliation(i)=new(i)
end
==(u::GuessAffiliation, v::GuessAffiliation) = u.community==v.community
hash(a::GuessAffiliation) = hash(a.community)

# reward
reward(p::CommunityAffiliationPOMDP, s, a::ProbeNode) = 0.0
function reward(p::CommunityAffiliationPOMDP, s, a::GuessAffiliation)
    return p.reward*(a.community==s.labels[p.node])
end

# state transitions
create_state(::CommunityAffiliationPOMDP) = CommunityNetwork()
create_transition_distribution(p::CommunityAffiliationPOMDP) = CommunityNetwork()
function transition(p::CommunityAffiliationPOMDP, s, a, d=create_transition_distribution(p))
    d.graph = s.graph
    d.n_comms = s.n_comms
    d.labels = s.labels
    d.budget = s.budget-1
    return d
end
function Base.rand!(rng::AbstractRNG, s, d::CommunityNetwork)
    s.graph = d.graph
    s.n_comms = d.n_comms
    s.labels = d.labels
    s.budget = d.budget
end

# observations
type Neighborhood <: AbstractDistribution
    neighbors::Vector{Int}
    Neighborhood() = new()
end
==(u::Neighborhood, v::Neighborhood) = u.neighbors==v.neighbors
hash(o::Neighborhood) = hash(o.neighbors)
create_observation(::CommunityAffiliationPOMDP) = Neighborhood()
create_observation_distribution(::CommunityAffiliationPOMDP) = Neighborhood()
function observation(p::CommunityAffiliationPOMDP, s, a::ProbeNode, d=create_observation_distribution(p))
    d.neighbors = neighbors(s.graph, a.node)
    return d
end
function observation(p::CommunityAffiliationPOMDP, s, a::GuessAffiliation, d=create_observation_distribution(p))
    d.neighbors = []
    return d
end
function Base.rand!(rng::AbstractRNG, o, d::Neighborhood)
    o.neighbors = d.neighbors
end

discount(::CommunityAffiliationPOMDP) = 1.0

include("heuristic_belief.jl")
include("heuristics.jl")
include("pomdp.jl")
include("problems.jl")
include("visualization.jl")

end # module