text:
export AffineSubspaceHomotopy, set_subspaces!

"""
    AffineSubspaceHomotopy(F::System, V::AffineSubspace, W::AffineSubspace)
    AffineSubspaceHomotopy(F::AbstractSystem, V::AffineSubspace, W::AffineSubspace)

Creates a homotopy ``H(x,t) = (F ∩ γ(t))(x)`` where ``γ(t)`` is a family of affine subspaces
such that ``H(x,1) = (F ∩ V)(x)`` and ``H(x,0) = (F ∩ W)(x)``.
Here ``γ(t)`` is the geodesic between `V` and `W` in the affine Grassmanian, i.e.,
it is the curve of minimal length connecting `V` and `W`.
See also [`AffineSubspace`](@ref) and [`geodesic`](@ref) and the references therein.
"""
Base.@kwdef struct AffineSubspaceHomotopy{S<:AbstractSystem} <: AbstractHomotopy
    system::S

    start::AffineSubspace{ComplexF64}
    target::AffineSubspace{ComplexF64}

    J::Matrix{ComplexF64}
    Q::Matrix{ComplexF64}
    Q_cos::Matrix{ComplexF64}
    Θ::Vector{Float64}
    U::Matrix{ComplexF64}
    γ1::Matrix{ComplexF64}

    x::Vector{ComplexF64}
    ẋ::Vector{ComplexF64}
    x_high::Vector{ComplexDF64}
    t_cache::Base.RefValue{ComplexF64}
    # For AD
    taylor_t_cache::Base.RefValue{ComplexF64}
    taylor_γ::NTuple{4,Matrix{ComplexF64}}
    v::Vector{ComplexF64}
    tx³::TaylorVector{4,ComplexF64}
    tx²::TaylorVector{3,ComplexF64}
    tx¹::TaylorVector{2,ComplexF64}
end

## Implementation details

# The computation is performed using (implicit) Stiefel coordinates. If `size(F) = m, n`
# and `V` and `W` are affine subspaces of dimension `k, then the computation is performed in
# the affine Grassmanian Graff(k,n) embedded in the Grassmanian(k+1,n+1) using the system
# ``[F(γ(t)v); (γ(t)v)[end] - 1]``. Here the `(γ(t)v)[k+1] - 1` condition ensures that
# we are in the affine Grassmanian and that `(γ(t)v)[1:k]` is the correct value in \C^n.

AffineSubspaceHomotopy(F::ModelKit.System, start::AffineSubspace, target::AffineSubspace) =
    AffineSubspaceHomotopy(ModelKitSystem(F), start, target)

function AffineSubspaceHomotopy(
    system::AbstractSystem,
    start::AffineSubspace,
    target::AffineSubspace,
)

    Q, Θ, U = geodesic_svd(target, start)
    Q_cos = target.intrinsic.Y * U
    tx³ = TaylorVector{4}(ComplexF64, size(Q, 1))

    AffineSubspaceHomotopy(
        system = system,
        start = copy(start),
        target = copy(target),
        J = zeros(ComplexF64, size(system) .+ (1, 1)),
        Q = Q,
        Q_cos = Q_cos,
        Θ = Θ,
        U = U,
        γ1 = Q_cos * LA.diagm(0 => cos.(Θ)) + Q * LA.diagm(0 => sin.(Θ)),
        x = zeros(ComplexF64, size(Q, 1)),
        ẋ = zeros(ComplexF64, size(Q, 1)),
        x_high = zeros(ComplexDF64, size(Q, 1)),
        t_cache = Ref(complex(NaN, NaN)),
        taylor_t_cache = Ref(complex(NaN, NaN)),
        taylor_γ = tuple((similar(Q) for i = 0:3)...),
        v = zeros(ComplexF64, size(Q, 2)),
        tx³ = tx³,
        tx² = TaylorVector{3}(tx³),
        tx¹ = TaylorVector{2}(tx³),
    )
end
Base.size(H::AffineSubspaceHomotopy) = (size(H.system)[1] + 1, dim(H.start) + 1)

"""
    set_subspaces!(H::AffineSubspaceHomotopy, start::AffineSubspace, target::AffineSubspace)

Update the homotopy `H` to track from the affine subspace `start` to `target`.
"""
function set_subspaces!(
    H::AffineSubspaceHomotopy,
    start::AffineSubspace,
    target::AffineSubspace,
)
    Q, Θ, U = geodesic_svd(target, start)
    copy!(H.start, start)
    copy!(H.target, target)
    H.Q .= Q
    H.Θ .= Θ
    H.U .= U
    LA.mul!(H.Q_cos, target.intrinsic.Y, U)
    LA.mul!(H.γ1, Q, LA.diagm(0 => sin.(Θ)))
    LA.mul!(H.γ1, H.Q_cos, LA.diagm(0 => cos.(Θ)), true, true)
    H
end

function γ!(H::AffineSubspaceHomotopy, t::Number)
    H.t_cache[] != t || return first(H.taylor_γ)
    if isreal(t)
        _γ!(H, real(t))
    else
        _γ!(H, t)
    end
    H.t_cache[] = t

    first(H.taylor_γ)
end
@inline function _γ!(H::AffineSubspaceHomotopy, t::Number)
    @unpack Q, Q_cos, Θ = H
    γ = first(H.taylor_γ)
    n, k = size(γ)
    @inbounds for j = 1:k
        Θⱼ = Θ[j]
        s, c = sincos(t * Θⱼ)
        for i = 1:n
            γ[i, j] = Q_cos[i, j] * c + Q[i, j] * s
        end
    end
    γ
end

γ̇!(H::AffineSubspaceHomotopy, t::Number) = isreal(t) ? _γ̇!(H, real(t)) : _γ̇!(H, t)
@inline function _γ̇!(H::AffineSubspaceHomotopy, t::Number)
    @unpack Q, Q_cos, Θ = H
    _, γ̇ = H.taylor_γ
    n, k = size(γ̇)
    @inbounds for j = 1:k
        Θⱼ = Θ[j]
        s, c = sincos(t * Θⱼ)
        ċ = -s * Θⱼ
        ṡ = c * Θⱼ
        for i = 1:n
            γ̇[i, j] = Q_cos[i, j] * ċ + Q[i, j] * ṡ
        end
    end
    γ̇
end

function set_solution!(u::Vector, H::AffineSubspaceHomotopy, x::AbstractVector, t)
    (length(x) == length(H.x) - 1) ||
        throw(ArgumentError("Cannot set solution. Expected extrinsic coordinates."))
    for i = 1:length(x)
        H.x[i] = x[i]
    end
    H.x[end] = 1

    if isone(t)
        LA.mul!(u, H.γ1', H.x)
    elseif iszero(t)
        LA.mul!(u, H.Q_cos', H.x)
    else
        LA.mul!(u, γ!(H, t), H.x)
    end
end

function get_solution(H::AffineSubspaceHomotopy, u::AbstractVector, t)
    if isone(t)
        (@view H.γ1[1:end-1, :]) * u
    elseif iszero(t)
        (@view H.Q_cos[1:end-1, :]) * u
    else
        γ = γ!(H, t)
        (@view γ[1:end-1, :]) * u
    end
end

function evaluate!(u, H::AffineSubspaceHomotopy, v::AbstractVector, t)
    γ = γ!(H, t)
    n = first(size(H.system))
    if eltype(v) isa ComplexDF64
        LA.mul!(H.x_high, γ, v)
        evaluate!(u, H.system, H.x_high)
        u[n+1] = H.x_high[end] - 1.0
    else
        LA.mul!(H.x, γ, v)
        evaluate!(u, H.system, H.x)
        u[n+1] = H.x[end] - 1.0
    end
    u
end

function evaluate_and_jacobian!(u, U, H::AffineSubspaceHomotopy, v::AbstractVector, t)
    γ = γ!(H, t)

    LA.mul!(H.x, γ, v)
    evaluate_and_jacobian!(u, H.J, H.system, H.x)
    LA.mul!(U, H.J, γ)

    n = first(size(H.system))
    u[n+1] = H.x[end] - 1

    m = length(v)
    for j = 1:m
        U[n+1, j] = γ[end, j]
    end

    nothing
end

function taylor!(u, ::Val{1}, H::AffineSubspaceHomotopy, v, t)
    γ = γ!(H, t)
    γ̇ = γ̇!(H, t)

    # apply chain rule
    #    d/dt [F(γ(t)v); (γ(t)v)[end] - 1] = [J_F(γ(t)v)* γ̇(t)*v;  (γ̇(t)v)[end]]
    LA.mul!(H.x, γ, v)
    LA.mul!(H.ẋ, γ̇, v)
    evaluate_and_jacobian!(u, H.J, H.system, H.x)
    LA.mul!(u, H.J, H.ẋ)
    M = size(H, 1)
    u[M] = H.ẋ[end]

    u
end

function _taylor_γ!(H::AffineSubspaceHomotopy, t::Number)
    @unpack Q, Q_cos, Θ, U, taylor_γ = H

    γ, γ¹, γ², γ³ = taylor_γ
    n, k = size(γ)
    @inbounds for j = 1:k
        Θⱼ = Θ[j]
        s, c = sincos(t * Θⱼ)
        c¹ = -s * Θⱼ
        s¹ = c * Θⱼ
        Θⱼ_2 = 0.5 * Θⱼ^2
        s² = -s * Θⱼ_2
        c² = -c * Θⱼ_2
        Θⱼ_3 = Θⱼ_2 * Θⱼ / 3
        s³ = -c * Θⱼ_3
        c³ = s * Θⱼ_3
        for i = 1:n
            γ[i, j] = Q_cos[i, j] * c + Q[i, j] * s
            γ¹[i, j] = Q_cos[i, j] * c¹ + Q[i, j] * s¹
            γ²[i, j] = Q_cos[i, j] * c² + Q[i, j] * s²
            γ³[i, j] = Q_cos[i, j] * c³ + Q[i, j] * s³
        end
    end

end

function taylor_γ!(H::AffineSubspaceHomotopy, t::Number)
    H.taylor_t_cache[] != t || return H.taylor_γ

    if isreal(t)
        _taylor_γ!(H, real(t))
    else
        _taylor_γ!(H, t)
    end
    H.taylor_t_cache[] = t

    H.taylor_γ
end

function taylor!(u, v::Val{2}, H::AffineSubspaceHomotopy, tv::TaylorVector, t, incr::Bool)
    γ, γ¹, γ², γ³ = taylor_γ!(H, t)
    x, x¹, x² = vectors(H.tx²)
    v, v¹ = vectors(tv)

    if !incr
        LA.mul!(x, γ, v)
        LA.mul!(x¹, γ¹, v)
    end

    H.v .= v¹
    LA.mul!(H.x, γ, v¹)
    x¹ .+= H.x

    LA.mul!(H.x, γ¹, H.v)
    H.v .= v
    LA.mul!(H.x, γ², H.v, true, true)
    x² .= H.x

    taylor!(u, Val(2), H.system, H.tx²)

    n = first(size(H.system))
    u[n+1] = x²[end]

    u
end

function taylor!(u, v::Val{3}, H::AffineSubspaceHomotopy, tv::TaylorVector, t, incr::Bool)
    γ, γ¹, γ², γ³ = taylor_γ!(H, t)
    x, x¹, x², x³ = vectors(H.tx³)
    v, v¹, v² = vectors(tv)

    if !incr
        LA.mul!(x, γ, v)
        LA.mul!(x¹, γ¹, v)
        LA.mul!(x¹, γ, v¹, true, true)
        LA.mul!(x², γ², v)
        LA.mul!(x², γ¹, v¹, true, true)
    end

    H.v .= v²
    LA.mul!(H.x, γ, H.v)
    x² .+= H.x

    LA.mul!(H.x, γ¹, H.v)
    H.v .= v
    LA.mul!(H.x, γ³, H.v, true, true)
    H.v .= v¹
    LA.mul!(H.x, γ², H.v, true, true)
    x³ .= H.x

    taylor!(u, Val(3), H.system, H.tx³)

    n = first(size(H.system))
    u[n+1] = x³[end]

    u
end


text:
"""
    Coefficients:
"""
wavenumber(i,j,d) = oftype(d.lx,(2π/d.lx)*i),oftype(d.ly,(2π/d.ly)*j)

function acoeffs(prob::BetaPlane{T,PointJet{T}}) where T<:AbstractFloat
    d,c,f = prob.d,prob.c,prob.f
    y = LinRange(-d.ly/2,d.ly/2,2d.ny-1)
    if(f.τ ≠ 1/c.μ) @warn "Relxation time τ is not equal to 1/μ" end
    u = f.τ > 0.0 ? -(f.Ξ*c.Ω/f.τ)*tanh.(y/f.Δθ) : zeros(Complex{T},2d.ny-1)
    s = convert(eltype(u),2/(2*d.ny-1))
    s*fftshift(fft(u))
end

function acoeffs(prob::BetaPlane{T,Kolmogorov{T}}) where T<:AbstractFloat
    d,f = prob.d,prob.f
    A = zeros(Complex{T},2d.ny-1)
    A[d.ny+1] = A[d.ny-1] = f.A₁
    A[d.ny+4] = A[d.ny-4] = 4*f.A₄
    A
end

acoeffs(prob) = zeros(Complex{eltype(prob)},2prob.d.ny-1)

function bcoeffs(prob)
    d,c = prob.d,prob.c
    (nx,ny) = size(d)
    B = zeros(Complex{eltype(prob)},2ny-1,nx)
    β = convert(Complex{eltype(prob)},2*c.Ω*cos(deg2rad(c.θ)))
    k₄ = norm([wavenumber(nx-1,ny-1,d)])
    α₄ = 4
    @inbounds for m = 0:nx-1
        nmin = m == 0 ? 1 : -ny+1
        @inbounds for n=nmin:ny-1
            kx,ky = wavenumber(m,n,d)
            k = (kx^2+ky^2)^0.5
            B[n+d.ny,m+1] += im*β*kx/k^2 - c.μ - c.ν*k^2 - c.ν₄*(k^2/k₄^2)^α₄
        end
    end
    B
end

function ccoeffs(prob,eqs::NL)
    d,c = prob.d,prob.c
    (nx,ny) = size(d)
    Cp = zeros(eltype(prob),2ny-1,nx,2ny-1,nx)
    Cm = zeros(eltype(prob),2ny-1,nx,2ny-1,nx)

    (c.linear == true) && return Cp,Cm

    # ++ interactions note: +0 has only (0,+n)
    @inbounds for m1=1:nx-1
        @inbounds for n1=-ny+1:ny-1
            @inbounds for m2=0:min(m1,nx-1-m1)
                n2min = m2 == 0 ? 1 : -ny+1
                @inbounds for n2=max(n2min,-ny+1-n1):min(ny-1,ny-1-n1)
                    px,py = wavenumber(m1,n1,d)
                    qx,qy = wavenumber(m2,n2,d)
                    if m1 == m2
                        Cp[n2+ny,m2+1,n1+ny,m1+1] = -(px*qy - qx*py)/(px^2 + py^2)
                    else
                        Cp[n2+ny,m2+1,n1+ny,m1+1] = -(px*qy - qx*py)*(1.0/(px^2 + py^2) - 1.0/(qx^2 + qy^2))
                    end
                end
            end
        end
    end
    # +- interactions note: - includes (0,-n) because it is conj(0,n)
    @inbounds for m1=1:nx-1
        @inbounds for n1=-ny+1:ny-1
            @inbounds for m2=0:m1
                n2min = m2 == 0 ? 1 : -ny+1
                n2max = m2 == m1 ? n1 - 1 : ny-1
                @inbounds for n2=max(n2min,n1-ny+1):1:min(n2max,n1+ny-1)
                    px,py = wavenumber(m1,n1,d)
                    qx,qy = wavenumber(m2,n2,d)
                    Cm[n2+ny,m2+1,n1+ny,m1+1] = (px*qy - qx*py)*(1.0/(px^2 + py^2) - 1.0/(qx^2 + qy^2))
                end
            end
        end
    end
    Cp,Cm
end

function ccoeffs(prob,eqs::Union{GQL,CE2,GCE2})
    d,c = prob.d,prob.c
    (nx,ny) = size(d)
    Λ = typeof(eqs) == CE2 ? 0 : eqs.Λ
    Cp = zeros(eltype(prob),2ny-1,nx,2ny-1,nx)
    Cm = zeros(eltype(prob),2ny-1,nx,2ny-1,nx)

    (c.linear == true) && return Cp,Cm

    # L + L = L
    @inbounds for m1=0:Λ
        n1min = m1 == 0 ? 1 : -ny+1
        @inbounds for n1=n1min:ny-1
            @inbounds for m2=0:min(m1,Λ-m1)
                n2min = m2 == 0 ? 1 : -ny+1
                @inbounds for n2=max(n2min,-ny+1-n1):min(ny-1,ny-1-n1)
                    px,py = wavenumber(m1,n1,d)
                    qx,qy = wavenumber(m2,n2,d)
                    if m1 == m2
                        Cp[n2+ny,m2+1,n1+ny,m1+1] = -(px*qy - qx*py)/(px^2 + py^2)
                    else
                        Cp[n2+ny,m2+1,n1+ny,m1+1] = -(px*qy - qx*py)*(1.0/(px^2 + py^2) - 1.0/(qx^2 + qy^2))
                    end
                end
            end
        end
    end
    # L - L = L
    # note: -L should always include (0,-n)
    @inbounds for m1=0:Λ
        n1min = m1 == 0 ? 1 : -ny+1
        @inbounds for n1=n1min:ny-1
            @inbounds for m2=0:m1
                n2min = m2 == 0 ? 1 : -ny+1
                n2max = m2 == m1 ? n1 - 1 : ny-1
                @inbounds for n2=max(n2min,n1-ny+1):1:min(n2max,n1+ny-1)
                    px,py = wavenumber(m1,n1,d)
                    qx,qy = wavenumber(m2,n2,d)
                    Cm[n2+ny,m2+1,n1+ny,m1+1] = (px*qy - qx*py)*(1.0/(px^2 + py^2) - 1.0/(qx^2 + qy^2))
                end
            end
        end
    end
    # H - H = L
    @inbounds for m1=Λ+1:nx-1
        @inbounds for n1=-ny+1:ny-1
            @inbounds for m2=max(Λ+1,m1-Λ):m1
                n2max = m2 == m1 ? n1 - 1 : ny-1
                @inbounds for n2=max(-ny+1,n1-ny+1):min(n2max,n1+ny-1)
                    px,py = wavenumber(m1,n1,d)
                    qx,qy = wavenumber(m2,n2,d)
                    Cm[n2+ny,m2+1,n1+ny,m1+1] = (px*qy - qx*py)*(1.0/(px^2 + py^2) - 1.0/(qx^2 + qy^2))
                end
            end
        end
    end
    # H + L = H
    @inbounds for m1=Λ+1:nx-1
        @inbounds for n1=-ny+1:ny-1
            @inbounds for m2=0:min(nx-1-m1,Λ)
                n2min = m2 == 0 ? 1 : -ny+1
                @inbounds for n2=max(n2min,-ny+1-n1):1:min(ny-1,ny-1-n1)
                    px,py = wavenumber(m1,n1,d)
                    qx,qy = wavenumber(m2,n2,d)
                    Cp[n2+ny,m2+1,n1+ny,m1+1] = -(px*qy - qx*py)*(1.0/(px^2 + py^2) - 1.0/(qx^2 + qy^2))
                end
            end
        end
    end
    # H - L = H
    # note: -L should always include (0,-n)
    @inbounds for m1=Λ+1:nx-1
        @inbounds for n1=-ny+1:ny-1
            @inbounds for m2=0:min(Λ,m1 - Λ - 1)
                n2min = m2 == 0 ? 1 : -ny+1
                @inbounds for n2=max(n2min,n1-ny+1):min(ny-1,n1+ny-1)
                    px,py = wavenumber(m1,n1,d)
                    qx,qy = wavenumber(m2,n2,d)
                    Cm[n2+ny,m2+1,n1+ny,m1+1] = (px*qy - qx*py)*(1.0/(px^2 + py^2) - 1.0/(qx^2 + qy^2))
                end
            end
        end
    end
    Cp,Cm
end

fcoeffs(prob,eqs::Union{NL,GQL}) = zeros(eltype(prob),2prob.d.ny-1,prob.d.nx)
fcoeffs(prob::BetaPlane{T,PointJet{T}},eqs::GCE2) where T = ArrayPartition(zeros(T,2prob.d.ny-1,eqs.Λ+1),zeros(T,2prob.d.ny-1,prob.d.nx-eqs.Λ,2prob.d.ny-1,prob.d.nx-eqs.Λ))
fcoeffs(prob::BetaPlane{T,Kolmogorov{T}},eqs::GCE2) where T = ArrayPartition(zeros(T,2prob.d.ny-1,eqs.Λ+1),zeros(T,2prob.d.ny-1,prob.d.nx-eqs.Λ,2prob.d.ny-1,prob.d.nx-eqs.Λ))
fcoeffs(prob::BetaPlane{T,PointJet{T}},eqs::CE2) where T = zeros(T,2prob.d.ny-1,2prob.d.ny-1,prob.d.nx-1)
fcoeffs(prob::BetaPlane{T,Kolmogorov{T}},eqs::CE2) where T = zeros(T,2prob.d.ny-1,2prob.d.ny-1,prob.d.nx-1)

stochamp(ε::T,kf::Int,N::T) where T<:AbstractFloat = convert(T,sqrt(2ε*kf^2)/sqrt(N))
stochcorr(ε::T,kf::Int,dk::Int) where T<:AbstractFloat = convert(T,2π*ε*kf/(32.0dk))
stochcorr(ε::T,kf::Int,N::T) where T<:AbstractFloat = convert(T,4*π*π*ε*kf/N/8.0)

stochamp2(ε::T,kf::Int,N::T) where T<:AbstractFloat = convert(T,sqrt(ε/N)*kf)
stochcorr2(ε::T,kf::Int,N::T) where T<:AbstractFloat = convert(T,4ε*kf/N)

function fcoeffs(prob::BetaPlane{T,Stochastic{T}},eqs::Union{NL,GQL}) where T
    d,f = prob.d,prob.f
    (nx,ny),Λ = size(d),lambda(prob,eqs)
    F = zeros(T,2ny-1,nx)
    for m=1:nx-1 # should 0 be included?
        for n=-ny+1:ny-1
            k = f.isotropic == true ? (m^2 + n^2)^0.5 : m
            if (f.kf - f.dk < k < f.kf + f.dk) F[n+ny,m+1] = one(T) end
        end
    end
    Nf = sum(F)
    # if (sum(F) ≥ 1.0) F .= F * stochamp2(f.ε,f.kf,Nf) end # this is dt unaware - dist contains dt
    if (Nf ≥ 1.0) F .= F * stochamp(f.ε,f.kf,Nf) end # this is dt unaware - dist contains dt
    return F
end

function fcoeffs(prob::BetaPlane{T,Stochastic{T}},eqs::CE2) where T
    d,f = prob.d,prob.f
    (nx,ny) = size(d)
    F = zeros(T,2ny-1,2ny-1,nx-1)
    for m=1:nx-1 # should 0 be included?
        for n=-ny+1:ny-1
            k = f.isotropic == true ? (m^2 + n^2)^0.5 : m
            if(f.kf - f.dk < k < f.kf + f.dk)
                F[n+ny,n+ny,m] = one(T)
            end
        end
    end
    Nf = sum(F)
    # F .= F * stochcorr2(f.ε,f.kf,Nf) # this is dt unaware - dist contains dt
    F .= F * stochcorr(f.ε,f.kf,Nf) # this is dt unaware - dist contains dt
    # F .= F * stochcorr(f.ε,f.kf,f.dk)
    return F
end

function fcoeffs(prob::BetaPlane{T,Stochastic{T}},eqs::GCE2) where T
    d,f = prob.d,prob.f
    (nx,ny),Λ = size(d),lambda(prob,eqs)
    F = ArrayPartition(zeros(T,2ny-1,Λ+1),zeros(T,2ny-1,nx-Λ,2ny-1,nx-Λ))
    for m=1:nx-1 # should 0 be included?
        for n=-ny+1:ny-1
            k = f.isotropic == true ? (m^2 + n^2)^0.5 : m
            if(f.kf - f.dk < k < f.kf + f.dk)
                m ≤ Λ ? F.x[1][n+ny,m+1] = one(T) : F.x[2][n+ny,m-Λ,n+ny,m-Λ] = one(T)
            end
        end
    end
    if (sum(F.x[1]) ≥ 1.0) F.x[1] .= F.x[1] * stochamp(f.ε,f.kf,sum(F.x[1])) end
    Nf = sum(F.x[2])
    F.x[2] .= F.x[2] * stochcorr(f.ε,f.kf,Nf) # this is dt unaware - dist contains dt
    return F
end

function fcoeffs2(prob::BetaPlane{T,Stochastic{T}},eqs::Union{NL,GQL}) where T
    (nx,ny) = size(prob.d)
    F0 = fcoeffs2(prob,CE2())
    F = zeros(T,2ny-1,nx)
    for m=1:nx-1
        for n=-ny+1:ny-1
            F[n+ny,m+1] = sqrt(2*F0[n+ny,n+ny,m])
        end
    end
    F
end

function fcoeffs2(prob::BetaPlane{T,Stochastic{T}},eqs::CE2) where T
    d,f = prob.d,prob.f
    (nx,ny) = size(d)
    F = zeros(T,2ny-1,2ny-1,nx-1)
    m1 = f.kf - f.dk
    m2 = f.kf + f.dk
    Nk = m2 - m1 + 1
    c  = 0.2
    for m=m1:m2
        Ck = zero(T)
        for n=-ny+1:ny-1
            kx,ky = 2π*m/d.lx,2π*n/d.ly
            F[n+ny,n+ny,m] = c^2*exp(-ky^2*c^2)
            Ck += F[n+ny,n+ny,m]/(kx^2+ky^2)
        end
        F[:,:,m] /= (2Ck*Nk)
    end
    F .* f.ε
end

function fcoeffs2(prob::BetaPlane{T,Stochastic{T}},eqs::GCE2) where T
    (nx,ny),Λ = size(prob.d),lambda(prob,eqs)
    F0 = fcoeffs2(prob,CE2())
    F = ArrayPartition(zeros(T,2ny-1,Λ+1),zeros(T,2ny-1,nx-Λ,2ny-1,nx-Λ))
    for m=1:nx-1
        for n=-ny+1:ny-1
            F.x[2][n+ny,m-Λ,n+ny,m-Λ] = F0[n+ny,n+ny,m]
        end
    end
    F
end

ext: jl


text:
function build_robot!(vis::Visualizer, model::Biped;
    r=0.035, r_contact=r * 8.0 / 7.0, color_opacity=1.0)

    r = convert(Float32, r)
    r_contact = convert(Float32, r_contact)

    body_mat = MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, color_opacity))
    contact_mat = MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0.0, color_opacity))

    default_background!(vis)

    torso = GeometryBasics.Cylinder(Point3f0(0.0), Point3f0(0.0, 0.0, model.l_torso), r)
    setobject!(vis[:robot]["torso"], torso, body_mat)

    thigh_1 = GeometryBasics.Cylinder(Point3f0(0.0), Point3f0(0.0, 0.0, model.l_thigh1), r)
    setobject!(vis[:robot]["thigh1"], thigh_1, body_mat)

    calf_1 = GeometryBasics.Cylinder(Point3f0(0.0), Point3f0(0.0, 0.0, model.l_calf1), r)
    setobject!(vis[:robot]["calf1"], calf_1, body_mat)

    foot_1 = GeometryBasics.Cylinder(Point3f0(0.0), Point3f0(0.0, 0.0, model.lt_foot1 + model.lh_foot1), r)
    setobject!(vis[:robot]["foot1"], foot_1, body_mat)

    thigh_2 = GeometryBasics.Cylinder(Point3f0(0.0), Point3f0(0.0, 0.0, model.l_thigh2), r)
    setobject!(vis[:robot]["thigh2"], thigh_2, body_mat)

    calf_2 = GeometryBasics.Cylinder(Point3f0(0.0), Point3f0(0.0, 0.0, model.l_calf2), r)
    setobject!(vis[:robot]["calf2"], calf_2, body_mat)

    foot_2 = GeometryBasics.Cylinder(Point3f0(0.0), Point3f0(0.0, 0.0, model.lt_foot2 + model.lh_foot2), r)
    setobject!(vis[:robot]["foot2"], foot_2, body_mat)

    setobject!(vis[:robot]["heel1"], GeometryBasics.Sphere(Point3f0(0.0), r_contact), contact_mat)
    setobject!(vis[:robot]["heel2"], GeometryBasics.Sphere(Point3f0(0.0), r_contact), contact_mat)
    setobject!(vis[:robot]["toe1"], GeometryBasics.Sphere(Point3f0(0.0), r_contact), contact_mat)
    setobject!(vis[:robot]["toe2"], GeometryBasics.Sphere(Point3f0(0.0), r_contact), contact_mat)
    setobject!(vis[:robot]["knee1"], GeometryBasics.Sphere(Point3f0(0.0), r), body_mat)
    setobject!(vis[:robot]["knee2"], GeometryBasics.Sphere(Point3f0(0.0), r), body_mat)
    setobject!(vis[:robot]["hip"], GeometryBasics.Sphere(Point3f0(0.0), r), body_mat)
    setobject!(vis[:robot]["torso_top"], GeometryBasics.Sphere(Point3f0(0.0), r), body_mat)

    return true
end

function set_robot!(vis::Visualizer, model::Biped, q;
   r=0.035, r_contact=r * 8.0 / 7.0)

    p_shift = [0.0; 0.0; r_contact]

    k_hip = kinematics_body(model, q, mode=:hip)
    p_hip = [k_hip[1]; 0.0; k_hip[2]] + p_shift

    k_torso = kinematics_body(model, q, mode=:ee)
    p_torso = [k_torso[1], 0.0, k_torso[2]] + p_shift

    k_thigh_1 = kinematics_thigh(model, q, leg=:leg1, mode=:ee)
    p_thigh_1 = [k_thigh_1[1], 0.0, k_thigh_1[2]] + p_shift

    k_calf_1 = kinematics_calf(model, q, leg=:leg1, mode=:ee)
    p_calf_1 = [k_calf_1[1], 0.0, k_calf_1[2]] + p_shift

    k_toe_1 = kinematics_foot(model, q, leg=:leg1, mode=:toe)
    p_toe_1 = [k_toe_1[1], 0.0, k_toe_1[2]] + p_shift

    k_heel_1 = kinematics_foot(model, q, leg=:leg1, mode=:heel)
    p_heel_1 = [k_heel_1[1], 0.0, k_heel_1[2]] + p_shift

    k_thigh_2 = kinematics_thigh(model, q, leg=:leg2, mode=:ee)
    p_thigh_2 = [k_thigh_2[1], 0.0, k_thigh_2[2]] + p_shift

    k_calf_2 = kinematics_calf(model, q, leg=:leg2, mode=:ee)
    p_calf_2 = [k_calf_2[1], 0.0, k_calf_2[2]] + p_shift

    k_toe_2 = kinematics_foot(model, q, leg=:leg2, mode=:toe)
    p_toe_2 = [k_toe_2[1], 0.0, k_toe_2[2]] + p_shift

    k_heel_2 = kinematics_foot(model, q, leg=:leg2, mode=:heel)
    p_heel_2 = [k_heel_2[1], 0.0, k_heel_2[2]] + p_shift

    settransform!(vis[:robot]["thigh1"], cable_transform(p_hip, p_thigh_1))
    settransform!(vis[:robot]["calf1"], cable_transform(p_thigh_1, p_calf_1))
    settransform!(vis[:robot]["foot1"], cable_transform(p_toe_1, p_heel_1))

    settransform!(vis[:robot]["thigh2"], cable_transform(p_hip, p_thigh_2))
    settransform!(vis[:robot]["calf2"], cable_transform(p_thigh_2, p_calf_2))
    settransform!(vis[:robot]["foot2"], cable_transform(p_toe_2, p_heel_2))

    settransform!(vis[:robot]["torso"], cable_transform(p_hip, p_torso))
    settransform!(vis[:robot]["heel1"], Translation(p_heel_1))
    settransform!(vis[:robot]["heel2"], Translation(p_heel_2))
    settransform!(vis[:robot]["toe1"], Translation(p_toe_1))
    settransform!(vis[:robot]["toe2"], Translation(p_toe_2))
    settransform!(vis[:robot]["knee1"], Translation(p_thigh_1))
    settransform!(vis[:robot]["knee2"], Translation(p_thigh_2))
    settransform!(vis[:robot]["hip"], Translation(p_hip))
    settransform!(vis[:robot]["torso_top"], Translation(p_torso))

    return true
end

function visualize!(vis, model::Biped, q;
	Δt=0.1, r=0.04, r_contact=0.04, color_opacity=1.0, fixed_camera=true)

	build_robot!(vis, model, r=r, r_contact=r_contact, color_opacity=color_opacity)

	anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

	for (t, qt) in enumerate(q)
		MeshCat.atframe(anim, t) do
			set_robot!(vis, model, qt)
		end
	end

	MeshCat.setanimation!(vis, anim)

    if fixed_camera
        settransform!(vis["/Cameras/default"],
            compose(Translation(0.0, -50.0, -1.0),LinearMap(RotZ(-pi / 2.0))))
        setprop!(vis["/Cameras/default/rotated/<object>"], "zoom", 25)
    end

    return true
end


text:
module HandlersWrapper

using JavaCall

const AbsoluteDate = @jimport org.orekit.time.AbsoluteDate
const Action = @jimport org.hipparchus.ode.events.Action
const Class = @jimport java.lang.Class
const ContinueOnEvent = @jimport org.orekit.propagation.events.handlers.ContinueOnEvent
const EventDetector = @jimport org.orekit.propagation.events.EventDetector
const EventHandler = @jimport org.orekit.propagation.events.handlers.EventHandler
const FieldAbsoluteDate = @jimport org.orekit.time.FieldAbsoluteDate
const FieldContinueOnEvent = @jimport org.orekit.propagation.events.handlers.FieldContinueOnEvent
const FieldEventDetector = @jimport org.orekit.propagation.events.FieldEventDetector
const FieldEventHandler = @jimport org.orekit.propagation.events.handlers.FieldEventHandler
const FieldRecordAndContinue = @jimport org.orekit.propagation.events.handlers.FieldRecordAndContinue
const FieldSpacecraftState = @jimport org.orekit.propagation.FieldSpacecraftState
const FieldStopOnDecreasing = @jimport org.orekit.propagation.events.handlers.FieldStopOnDecreasing
const FieldStopOnEvent = @jimport org.orekit.propagation.events.handlers.FieldStopOnEvent
const FieldStopOnIncreasing = @jimport org.orekit.propagation.events.handlers.FieldStopOnIncreasing
const JString = @jimport java.lang.String
const List = @jimport java.util.List
const Object = @jimport java.lang.Object
const RecordAndContinue = @jimport org.orekit.propagation.events.handlers.RecordAndContinue
const SpacecraftState = @jimport org.orekit.propagation.SpacecraftState
const StopOnDecreasing = @jimport org.orekit.propagation.events.handlers.StopOnDecreasing
const StopOnEvent = @jimport org.orekit.propagation.events.handlers.StopOnEvent
const StopOnIncreasing = @jimport org.orekit.propagation.events.handlers.StopOnIncreasing

export ContinueOnEvent
export EventHandler
export FieldContinueOnEvent
export FieldEventHandler
export FieldRecordAndContinue
export FieldStopOnDecreasing
export FieldStopOnEvent
export FieldStopOnIncreasing
export RecordAndContinue
export StopOnDecreasing
export StopOnEvent
export StopOnIncreasing
export clear
export equals
export event_occurred
export get_class
export get_events
export hash_code
export init
export notify
export notify_all
export reset_state
export to_string
export wait

include("continue_on_event.jl")
include("event_handler.jl")
include("field_continue_on_event.jl")
include("field_event_handler.jl")
include("field_record_and_continue.jl")
include("field_stop_on_decreasing.jl")
include("field_stop_on_event.jl")
include("field_stop_on_increasing.jl")
include("record_and_continue.jl")
include("stop_on_decreasing.jl")
include("stop_on_event.jl")
include("stop_on_increasing.jl")

end


text:
"""
    GPU <: Processor

Abstract supertype for GPU processors.
"""
abstract type GPU <: Processor end

"""
    CPUGPU <: GPU

    CPUGPU()

Uses the CUDA GPU code on CPU using KernelAbstractions, to test it.
"""
struct CPUGPU{L} <: GPU
    spinlock::L
end
CPUGPU() = CPUGPU(Base.Threads.SpinLock())
Base.Threads.lock(opt::CPUGPU) = lock(opt.spinlock)
Base.Threads.unlock(opt::CPUGPU) = unlock(opt.spinlock)

kernel_setup(::CPUGPU) = KernelAbstractions.CPU(), 1

function maprule!(
    data::AbstractSimData, proc::GPU, opt, ruletype::Val{<:Rule}, rule, rkeys, wkeys
)
    kernel! = ka_rule_kernel!(kernel_setup(proc)...)
    kernel!(data, ruletype, rule, rkeys, wkeys; ndrange=gridsize(data)) |> wait
    return nothing
end

# ka_rule_kernel!
# Runs cell_kernel! on GPU after retrieving the global index
# and setting the neighborhood buffer to a SArray window retrieved
# from the first (neighborhood) grid
@kernel function ka_rule_kernel!(data, ruletype::Val{<:NeighborhoodRule}, rule, rkeys, wkeys)
    I = @index(Global, NTuple)
    neighborhood_kernel!(data, _firstgrid(data, rkeys), ruletype, rule, rkeys, wkeys, I...)
    nothing
end
@kernel function ka_rule_kernel!(data, ruletype::Val, rule, rkeys, wkeys)
    I = @index(Global, NTuple)
    cell_kernel!(data, ruletype, rule, rkeys, wkeys, I...)
    nothing
end

### Indexing. UNSAFE / LOCKS required

# This is not safe for general use.
# It can be used where only identical transformations of a cell
# can happen from any other cell, such as setting all 1s to 2.
@propagate_inbounds function _setindex!(d::WritableGridData, opt::GPU, x, I...)
    dest(d)[I...] = x
end


text:
function disc_phantom(img, a, b, r)
    height, width = size(img)
    #find the possibility which has most surrounding pixels
    disc = zeros(Float64, (height,width))
    for x = 1:width
        for y = 1:height
            disc[y,x] = (x-a)^2 + (y-b)^2 < r^2 ? 1.0 : 0.0
        end
    end
    return disc
end

max_stars_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_issues_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_forks_repo_name: JuliaTomo/DynamicTomographyWithFlow


text:
function getcombinations(patterns::Vector{Interval}, weights::Vector{Interval})
	N = 2length(patterns)
	inners = 2 .^ collect(0:N-1)
	outers = reverse(inners)
	variables = [patterns; weights]
	combs = Float64[]
	for (inner, outer, var) in zip(inners, outers, variables)
		endpoints = [var.left, var.right]
		c = repeat(endpoints, inner=inner, outer=outer)
		if Base.isempty(combs)
			combs = c
		else
			combs = hcat(combs, c)
		end
	end
	return combs
end

function fwa(levelendpoints::Matrix{Float64})
	n = Int(size(levelendpoints)[2] / 2)
	averages = []
	for endpoints in eachrow(levelendpoints)
		numerator = 0
		denominator = 0
		for i = 1:n
			j = i + n
			x = endpoints[i]
			w = endpoints[j]
			numerator += x * w
			denominator += w
		end
		avg = denominator == 0 ? nothing : numerator / denominator
		if !isnothing(avg)
			push!(averages, avg)
		end
	end
	return Interval(minimum(averages), maximum(averages))
end

function fuzzy_weighted_average(X⃗::FuzzyVector, W⃗::FuzzyVector)
	#FIXME: more than 2-d patterns
	# pattern_vectors = [X.grades for X in X⃗.numbers]
	# weight_vectors = [W.grades for W in W⃗.numbers]

	xx = []
	ww = []
	A₁ = X⃗.numbers[1]
	A₂ = X⃗.numbers[2]
	W₁ = W⃗.numbers[1]
	W₂ = W⃗.numbers[2]
	for i = 1:length(A₁)
		push!(xx, [A₁.grades[i], A₂.grades[i]])
		push!(ww, [W₁.grades[i], W₂.grades[i]])
	end

	combinations = FuzzySets.getcombinations.(xx, ww)
	grades = FuzzySets.fwa.(combinations)
	return FuzzySets.FuzzyNumber(A₁.levels, grades)
end


text:
using ..Registrator: decodeb64

function is_pfile_parseable(c::AbstractString)
    @debug("Checking whether Project.toml is non-empty and parseable")
    if length(c) != 0
        try
            TOML.parse(c)
            return true, nothing
        catch ex
            if isa(ex, CompositeException) && isa(ex.exceptions[1], TOML.ParserError)
                err = "Error parsing project file"
                @debug(err)
                return false, err
            else
                rethrow(ex)
            end
        end
    else
        err = "Project file is empty"
        @debug(err)
        return false, err
    end
end

function is_pfile_nuv(c)
    @debug("Checking whether Project.toml contains name, uuid and version")
    ib = IOBuffer(c)

    try
        p = Pkg.Types.read_project(copy(ib))
        if p.name === nothing || p.uuid === nothing || p.version === nothing
            err = "Project file should contain name, uuid and version"
            @debug(err)
            return false, err
        elseif !isempty(p.version.prerelease)
            err = "Pre-release version not allowed"
            @debug(err)
            return false, err
        elseif p.version == v"0"
            err = "Package version must be greater than 0.0.0"
            @debug(err)
            return false, err
        end
    catch ex
        err = "Error reading Project.toml: $(ex.msg)"
        @debug(err)
        return false, err
    end

    return true, nothing
end

function is_pfile_valid(c::AbstractString)
    for f in [is_pfile_parseable, is_pfile_nuv]
        v, err = f(c)
        v || return v, err
    end
    return true, nothing
end

function verify_projectfile_from_sha(reponame, sha; auth=GitHub.AnonymousAuth())
    projectfile_contents = nothing
    projectfile_found = false
    projectfile_valid = false
    err = nothing
    @debug("Getting gitcommit object for sha")
    gcom = gitcommit(reponame, GitCommit(Dict("sha"=>sha)); auth=auth)
    @debug("Getting tree object for sha")
    t = tree(reponame, Tree(gcom.tree); auth=auth)

    for tr in t.tree
        if tr["path"] == "Project.toml"
            projectfile_found = true
            @debug("Project file found")

            @debug("Getting projectfile blob")
            if isa(auth, GitHub.AnonymousAuth)
                a = get_user_auth()
            else
                a = auth
            end
            b = blob(reponame, Blob(tr["sha"]); auth=a)

            @debug("Decoding base64 projectfile contents")
            projectfile_contents = decodeb64(b.content)

            @debug("Checking project file validity")
            projectfile_valid, err = is_pfile_valid(projectfile_contents)
            break
        end
    end

    return projectfile_contents, t.sha, projectfile_found, projectfile_valid, err
end


text:
# !! auto-generated code, do not edit !!
module apriltag_msgs
using RobotOSData.Messages
$(Expr(:using, :($(Expr(:., :(RobotOSData.StdMsgs))))))
$(Expr(:using, :($(Expr(:., :(RobotOSData.CommonMsgs))))))
struct Point <: Readable
    x::Float64
    y::Float64
end
struct AprilTagDetection <: Readable
    family::String
    id::Int32
    hamming::Int32
    goodness::Float32
    decision_margin::Float32
    centre::Point
    corners::SVector{4, Point}
    homography::SVector{9, Float64}
end
struct AprilTagDetectionArray <: Readable
    header::std_msgs.Header
    detections::Vector{AprilTagDetection}
end
end


text:
# Copyright 2020 Google LLC

# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at

#     https://www.apache.org/licenses/LICENSE-2.0

# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
lines = readlines("day07input.txt")

function parse_requirement(line)
    color_match = match(r"^([a-z ]+) bags contain", line)
    if color_match === nothing
        println(line)
        return
    end
    color = strip(color_match.captures[1])
    children_match = match(r"contain ([a-z, 0-9]+).$", line)
    children = children_match.captures[1]
    parts = split(children, ",")
    contained_colors = []
    for part in parts
        if contains(part, "no other bags")
            break
        else
            child_match = match(r"^([ 0-9]+)([a-z ]+) bag[s]{0,1}", part)
            num = parse(Int, child_match.captures[1])
            child_color = child_match.captures[2]
            push!(contained_colors, (num, child_color))
        end
    end
    return (color, contained_colors)
end

function parse_requirements(lines)
    requirements = Dict()
    for req ∈ map(parse_requirement, lines)
        parent, children = req
        if haskey(requirements, parent)
            throw(ErrorException("Non-unique key"))
        end
        requirements[parent] = children
    end
    return requirements
end

function invert_requirements(requirements)
    parents = Dict{AbstractString, AbstractArray{AbstractString}}()
    for color ∈ keys(requirements)
        for (n, child_color) ∈ requirements[color]
            if haskey(parents, child_color)
                push!(parents[child_color], color)
            else
                parents[child_color] = [color]
            end
        end
    end
    return Dict(key => Set(value) for (key, value) ∈ parents)
end

function find_possible_containers(requirements, target::AbstractString)
    parents = invert_requirements(requirements)

    to_visit = Vector()
    push!(to_visit, target)
    visited = Set{String}()
    while length(to_visit) > 0
        color_to_find = pop!(to_visit)
        push!(visited, color_to_find)
        if haskey(parents, color_to_find)
            push!(to_visit, setdiff(parents[color_to_find], visited)...)
        end
    end
    delete!(visited, target)
    return visited
end

function part1(requirements)
    find_possible_containers(requirements, "shiny gold") |> length
end

function count_bags_inside(requirements, parent::AbstractString)
    total = 0
    for child ∈ requirements[parent]
        n, color = child
        total += n
        total += n * count_bags_inside(requirements, color)
    end
    return total
end

function part2(requirements)
    count_bags_inside(requirements, "shiny gold")
end

requirements = parse_requirements(lines)
println("Part 1: ", part1(requirements))
println("Part 2: ", part2(requirements))


text:
using Performance
using BenchmarkTools
using Test

@testset "structs" begin
    ##
    ## https://docs.julialang.org/en/v1/manual/performance-tips/#Type-declarations
    ##
    # Unfortunately, this problem is not Revise-friendly. You'll have to restart your
    # Julia session when you make substantive changes to the source code.
    #
    # Also, this problem uses an AbstractArray type as an example, but the lessons
    # apply to almost any `struct`. One specific feature, though, is that
    # `AbstractArray` has a standard
    # [interface](https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-array)
    # that you may need to learn for this problem. Not very many types in Julia are this well
    # standardized, but it's a good model for how to think about composability
    # and extensability more generally, two key goals for "good" code.
    a = DefaultArray([10, 20, 30], 0)
    @test @inferred(a[2]) === 20
    @test @inferred(a[4]) === 0
    # This one is trickier. See https://docs.julialang.org/en/v1/manual/constructors/#Parametric-Constructors
    # for help. See `?promote_type` as a hint for crafting your solution.
    a = DefaultArray(Float32[10, 20, 30], 0)
    @test @inferred(a[2]) === 20.0f0
    @test @inferred(a[4]) === 0.0f0
    a = DefaultArray(Float32[10, 20, 30], nothing)
    @test eltype(a) === Union{Nothing, Float32}
end


text:
# This file is a part of RRRMC.jl. License is MIT: http://github.com/carlobaldassi/RRRMC.jl/LICENCE.md

module Empty

using ExtractMacro
using ..Interface
using ..Common
using ..QT

export GraphEmpty

import ..Interface: energy, delta_energy, neighbors

struct GraphEmpty <: SimpleGraph{Int}
    N::Int
    GraphEmpty(N::Integer) = return new(N)
end

@doc """
    GraphEmpty(N::Integer) <: SimpleGraph{Int}

A trivial `SimpleGraph` type with `N` free, non-interacting spins.
(The energy is always `0`.)

Only useful for testing/debugging purposes.
""" -> GraphEmpty()

energy(X::GraphEmpty, C::Config) = 0
delta_energy(X::GraphEmpty, C::Config, move::Int) = 0

neighbors(X::GraphEmpty, i::Int) = ()

end # module


text:
using Pkg

function addandpin(spec)
    Pkg.add(PackageSpec(; spec...))
    Pkg.pin(spec.name)
end

# packages without external dependencies
# that have to be manually instlled

pkg1 = [(name="StatsBase", version="0.26.0"),
        (name="TimeZones", version="0.8.1"),
        (name="BSON", version="0.2.1"),
        (name="Revise", version="0.7.12"),
        (name="Distributions", version="0.16.4"),
        (name="Clp", version="0.5.0"),
        (name="HTTP", version="0.7.1"),
        (name="Gumbo", version="0.5.1"),
        (name="StringEncodings", version="0.3.1"),
        (name="ZMQ", version="1.0.0"),
        (name="CodecZlib", version="0.5.0"),
        (name="JSON", version="0.19.0"),
        (name="BenchmarkTools", version="0.4.1"),
        (name="JuliaWebAPI", version="0.5.0"),
        (name="FileIO", version="1.0.2"),
        (name="ProfileView", version="0.4.0"),
        (name="StaticArrays", version="0.8.3"),
        (name="ForwardDiff", version="0.9.0"),
        (name="Optim", version="0.17.1"),
        (name="JuMP", version="0.18.4"),
        (name="JLD2", version="0.1.2"),
        (name="XLSX", version="0.4.2"),
        (name="Cbc", version="0.4.2"),
        (name="DataFrames", version="0.14.1"),
        (name="CSV", version="0.4.3"),
        (name="DataFramesMeta", version="0.4.0"),
        (name="Feather", version="0.5.0"),
        (name="FreqTables", version="0.3.0"),
        (name="OnlineStats", version="0.19.1"),
        (name="MySQL", version="0.7.0"),
        (name="Cascadia", version="0.4.0"),
        (name="UnicodePlots", version="0.3.1"),
        (name="ParallelDataTransfer", version="0.5.0")]

foreach(addandpin, pkg1)

# packages that depend on Python Anaconda
# whose manual installation is optional

pkg2 = [(name="Conda", version="1.0.2"),
        (name="PyCall", version="1.18.4"),
        (name="PyPlot", version="2.6.3"),
        (name="Plots", version="0.20.5"),
        (name="StatsPlots", version="0.8.1")]

foreach(addandpin, pkg2)

# packages that have external dependencies
# that have to be manually installed

pkg3 = [(name="RCall", version="0.12.1"),
        (name="JDBC", version="0.4.0"),
        (name="LibPQ", version="0.5.0"),
        (name="Gurobi", version="0.5.3")]

foreach(addandpin, pkg3)


text:
function fib150()
#
#  test polynomial suggested by Goedecker
#
    n = 150;
    p = fib(n)
    z = [-.9927187559253492;
         -.9918530882118113-.4147913261323188e-1*im;
         -.9918530882118113+.4147913261323188e-1*im;
         -.9892575836821723-.8288617318453414e-1*im;
         -.9892575836821723+.8288617318453414e-1*im;
         -.9849367355771876-.1241491541448754*im;
         -.9849367355771876+.1241491541448754*im;
         -.9788980239954073-.1651963566541725*im;
         -.9788980239954073+.1651963566541725*im;
         -.9711519029549269-.2059564344240104*im;
         -.9711519029549269+.2059564344240104*im;
         -.9617117823114389-.2463585368912788*im;
         -.9617117823114389+.2463585368912788*im;
         -.9505940045641273-.2863324315280673*im;
         -.9505940045641273+.2863324315280673*im;
         -.9378178165898683-.3258086250746164*im;
         -.9378178165898683+.3258086250746164*im;
         -.9234053363550660-.3647184834839322*im;
         -.9234053363550660+.3647184834839322*im;
         -.9073815146632301-.4029943503688435*im;
         -.9073815146632301+.4029943503688435*im;
         -.8897740920050947-.4405696637447960*im;
         -.8897740920050947+.4405696637447960*im;
         -.8706135505866639-.4773790708645437*im;
         -.8706135505866639+.4773790708645437*im;
         -.8499330616190360-.5133585409440989*im;
         -.8499330616190360+.5133585409440989*im;
         -.8277684279621956-.5484454755828406*im;
         -.8277684279621956+.5484454755828406*im;
         -.8041580222231623-.5825788166845422*im;
         -.8041580222231623+.5825788166845422*im;
         -.7791427204169229-.6156991516902623*im;
         -.7791427204169229+.6156991516902623*im;
         -.7527658313064592-.6477488159385299*im;
         -.7527658313064592+.6477488159385299*im;
         -.7250730215458899-.6786719919730538*im;
         -.7250730215458899+.6786719919730538*im;
         -.6961122367582722-.7084148056232615*im;
         -.6961122367582722+.7084148056232615*im;
         -.6659336186869521-.7369254186883471*im;
         -.6659336186869521+.7369254186883471*im;
         -.6345894185664953-.7641541180611353*im;
         -.6345894185664953+.7641541180611353*im;
         -.6021339068661845-.7900534011339715*im;
         -.6021339068661845+.7900534011339715*im;
         -.5686232795658151-.8145780573349891*im;
         -.5686232795658151+.8145780573349891*im;
         -.5341155611300741-.8376852456494883*im;
         -.5341155611300741+.8376852456494883*im;
         -.4986705043541354-.8593345679877696*im;
         -.4986705043541354+.8593345679877696*im;
         -.4623494872592649-.8794881382675888*im;
         -.4623494872592649+.8794881382675888*im;
         -.4252154072232028-.8981106470864277*im;
         -.4252154072232028+.8981106470864277*im;
         -.3873325725358873-.9151694218659986*im;
         -.3873325725358873+.9151694218659986*im;
         -.3487665915767304-.9306344823588107*im;
         -.3487665915767304+.9306344823588107*im;
         -.3095842598151608-.9444785914142191*im;
         -.3095842598151608+.9444785914142191*im;
         -.2698534448415497-.9566773009091463*im;
         -.2698534448415497+.9566773009091463*im;
         -.2296429696409571-.9672089927566195*im;
         -.2296429696409571+.9672089927566195*im;
         -.1890224943274280-.9760549149134032*im;
         -.1890224943274280+.9760549149134032*im;
         -.1480623965618780-.9831992123163535*im;
         -.1480623965618780+.9831992123163535*im;
         -.1068336508819937-.9886289526856903*im;
         -.1068336508819937+.9886289526856903*im;
         -.6540770717811572e-1-.9923341471422300*im;
         -.6540770717811572e-1+.9923341471422300*im;
         -.2385636855484892e-1-.9943077655947946*im;
         -.2385636855484892e-1+.9943077655947946*im;
         .1774833117573815e-1-.9945457468636005*im;
         .1774833117573815e-1+.9945457468636005*im;
         .5933425011689494e-1-.9930470035155510*im;
         .5933425011689494e-1+.9930470035155510*im;
         .1008292607663695-.9898134213981652*im;
         .1008292607663695+.9898134213981652*im;
         .1421613720845280-.9848498538705906*im;
         .1421613720845280+.9848498538705906*im;
         .1832588508751938-.9781641107430547*im;
         .1832588508751938+.9781641107430547*im;
         .2240503421931883-.9697669419505880*im;
         .2240503421931883+.9697669419505880*im;
         .2644649884815679-.9596720160034037*im;
         .2644649884815679+.9596720160034037*im;
         .3044325471292468-.9478958932755882*im;
         .3044325471292468+.9478958932755882*im;
         .3438835061260140-.9344579942165983*im;
         .3438835061260140+.9344579942165983*im;
         .3827491974769262-.9193805625975399*im;
         .3827491974769262+.9193805625975399*im;
         .4209619080219071-.9026886239377113*im;
         .4209619080219071+.9026886239377113*im;
         .4584549872896030-.8844099392981356*im;
         .4584549872896030+.8844099392981356*im;
         .4951629519980681-.8645749546799259*im;
         .4951629519980681+.8645749546799259*im;
         .5310215868002686-.8432167463288874*im;
         .5310215868002686+.8432167463288874*im;
         .5659680408622439-.8203709623267845*im;
         .5659680408622439+.8203709623267845*im;
         .5999409198599586-.7960757609476053*im;
         .5999409198599586+.7960757609476053*im;
         .6328803729932399-.7703717463775087*im;
         .6328803729932399+.7703717463775087*im;
         .6647281746501417-.7433019025433064*im;
         .6647281746501417+.7433019025433064*im;
         .6954278004244162-.7149115259686257*im;
         .6954278004244162+.7149115259686257*im;
         .7249244973086152-.6852481587793065*im;
         .7249244973086152+.6852481587793065*im;
         .7531653480767721-.6543615232056414*im;
         .7531653480767721+.6543615232056414*im;
         .7800993301600754-.6223034591666788*im;
         .7800993301600754+.6223034591666788*im;
         .8056773697373049-.5891278667468168*im;
         .8056773697373049+.5891278667468168*im;
         .8298523923415625-.5548906555452699*im;
         .8298523923415625+.5548906555452699*im;
         .8525793720535097-.5196497029280600*im;
         .8525793720535097+.5196497029280600*im;
         .8738153823196024-.4834648230425889*im;
         .8738153823196024+.4834648230425889*im;
         .8935196525762829-.4463977479358605*im;
         .8935196525762829+.4463977479358605*im;
         .9116536360890624-.4085121210914591*im;
         .9116536360890624+.4085121210914591*im;
         .9281810955442824-.3698735020064003*im;
         .9281810955442824+.3698735020064003*im;
         .9430682136510916-.3305493779521499*im;
         .9430682136510916+.3305493779521499*im;
         .9562837358786358-.2906091758178016*im;
         .9562837358786358+.2906091758178016*im;
         .9677991509344817-.2501242631678211*im;
         .9677991509344817+.2501242631678211*im;
         .9775889111839516-.2091679239550551*im;
         .9775889111839516+.2091679239550551*im;
         .9856306896705420-.1678152916898196*im;
         .9856306896705420+.1678152916898196*im;
         .9919056630215224-.1261432225318404*im;
         .9919056630215224+.1261432225318404*im;
         .9963988013915535-.8423009394365190e-1*im;
         .9963988013915535+.8423009394365190e-1*im;
         .9990991395939178-.4215552183469491e-1*im;
         .9990991395939178+.4215552183469491e-1*im;
         2.000000000000000];

    z = [z ones(n)];
    p, PolyZeros(z)
end


text:
module multimaterial_nas_examples
using FinEtools
using FinEtools.MeshExportModule: MESH
using FinEtoolsDeforLinear.AlgoDeforLinearModule
using FinEtoolsDeforLinear
using LinearAlgebra
using SparseArrays
using Arpack

function multimaterial_nas()
    println("""
    Vibration modes of block composed of two materials. NASTRAN input file.
    """)

    materials = Dict(
        1=>(name = "steel", E = 205000*phun("MPa"), nu = 0.3, rho = 7850*phun("KG*M^-3")),
        2=>(name = "aluminum", E = 70000*phun("MPa"), nu = 0.34, rho = 2700*phun("KG*M^-3"))
        )
    OmegaShift = (2*pi*100) ^ 2; # to resolve rigid body modes
    neigvs = 20;

    # No need to change anything below this line ##########
    MR = DeforModelRed3D
    output = import_NASTRAN("twoblocks.nas")
    fens, fesets, pids = output["fens"], output["fesets"], output["property_ids"]

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),3)) # displacement field

    numberdofs!(u)

    K = spzeros(u.nfreedofs, u.nfreedofs)
    M = spzeros(u.nfreedofs, u.nfreedofs)
    allfes = nothing
    for i in 1:length(fesets)
        pid = pids[i]
        @show E, nu, rho = materials[pid].E, materials[pid].nu, materials[pid].rho
        material = MatDeforElastIso(MR, rho, E, nu, 0.0)
        femm = FEMMDeforLinearESNICET4(MR, IntegDomain(fesets[i], NodalSimplexRule(3)), material)
        femm = associategeometry!(femm, geom)
        K += stiffness(femm, geom, u)
        M += mass(femm, geom, u)
        if allfes == nothing
            allfes = fesets[i]
        else
            allfes = cat(allfes, fesets[i])
        end
    end

    # eigs returns the nev requested eigenvalues in d, the corresponding Ritz vectors
    # v (only if ritzvec=true), the number of converged eigenvalues nconv, the number
    # of iterations niter and the number of matrix vector multiplications nmult, as
    # well as the final residual vector resid.

    if true
        d,v,nev,nconv = eigs(K+OmegaShift*M, M; nev=neigvs, which=:SM)
        d = d .- OmegaShift;
        fs = real(sqrt.(complex(d)))/(2*pi)
        println("Eigenvalues: $fs [Hz]")
        mode = 7
        scattersysvec!(u, v[:,mode])
        File =  "multimaterial_nas.vtk"
        vtkexportmesh(File, fens, allfes; vectors=[("mode$mode", u.values)])
        @async run(`"paraview.exe" $File`)
    end

    # Extract the boundary
    bfes = meshboundary(allfes)
    bconn = connasarray(bfes)

    MESH.write_MESH("multimaterial_nas.mesh", fens, bfes)

    true

end # multimaterial

function allrun()
    println("#####################################################")
    println("# multimaterial_nas ")
    multimaterial_nas()
    return true
end # function allrun

end # module multimaterial_nas_examples


text:
#=
Fundamental result of the entire algorithm:
for each point in the masked region of the target image,
the coordinates (in the corpus frame)
of the best matching patch.

When the algorithm is "with replacement"
the target image is also mutated as the algorithm proceeds.

It is not the same thing to execute without replacement,
and only then mutate (set the color of) the target image
with the color from the best match in the corpus.
Since replacing as you go along gives an ever changing comparand for searches.

But you could also say that the mutated target is a result.
Thats the most used result, say for the "retouching" application.

After the algorithm,
1) hasValue is all true
2) and mapFromTargetToCorpusPoints is valid CartesianIndex in the synth region
and uninitialized elsewhere.
=#

#=
The arrays are the same shape as the target.

The arrays are ordinary.
The mapFromTargetToCorpusPoints is sparse:
most values not initialized, only values in the synth region ever get values.
But we don't use an implementation that uses less storage when sparse.

The frame of the arrays is the same as the frame of the target image.

TODO use a subarray with the same frame, but less storage???
i.e. an array that encloses the synthesized region of the target
=#

# Depends on image.jl because methods take instance of MaskedImage
include("../image.jl")  # MaskedImage

# OLD array of CartesianIndex
# MapElementType =
# MapElementType = Int64

struct SynthResult{DimensionCount}
    #=
    Both the array, and its elements, have same DimensionCount.
    !!! The elements of the array are CartesianIndex.
    The eltype() is (CartesianIndex{DimensionCount})

    NOT AbstractArray.
    Resynthesizer chooses an implementation.
    =#

    # WAS
    mapFromTargetToCorpusPoints::Array{CartesianIndex{DimensionCount}, DimensionCount}

    # NEW array of linear index i.e. Int64
    # mapFromTargetToCorpusPoints::Array{MapElementType, DimensionCount}

    hasValue::BitArray{DimensionCount}
end

#=
Return the initial SynthResult

No points in the synth region have been synthesized.
The coordinates of the best match for points in the synth region is unknown.

In original C, prepare_target_sources()
=#
function initialSynthResult(
        targetImage::MaskedImage{ValueType, DimensionCount}
        ) where {ValueType, DimensionCount}
    @debug "SynthResult initializer called"

    #=
    Uninitialized array of CartesianIndex.
    Whose elements have same concrete type as the indices of targetImage
    e.g. CartesianIndex{N}

    TODO Really the indices of the corpusImage.
    For more generality, allow targetImage and corpusImage to have differing dimensions.
    =#
    # elementType = CartesianIndex{ndims(targetImage.image)}
    elementType = CartesianIndex{DimensionCount}

    # TODO this is doing one allocation for each element of the image
    # Why isn't it just one huge allocation?
    mapFromTargetToCorpusPoints = similar(targetImage.image, elementType) # {DimensionCount})

    # TODO understand what an undef is
    # mapFromTargetToCorpusPoints = Array{CartesianIndex{DimensionCount}}(undef, size(targetImage.image))
    #mapFromTargetToCorpusPoints = Array{MapElementType}(undef, size(targetImage.image))

    # boolean array false only in the synth region (the selected region give by the mask)
    hasValue = trues(size(targetImage.image))
    hasValue[targetImage.mask] .= false
    @debug "Initial hasValue map" hasValue

    result = SynthResult(mapFromTargetToCorpusPoints, hasValue)
    @debug "typeof SynthResult" typeof(result)
    # Useless  to assert isconcretetype(result) since all structs are concrete
    return result
end

#=
Setter method of the SynthResult object.
=#
function setBetterSearchResult(
        synthResult::SynthResult{DimensionCount},
        targetIndex::CartesianIndex{DimensionCount},
        corpusIndex::CartesianIndex{DimensionCount}
        )::Nothing where {DimensionCount}
    synthResult.hasValue[targetIndex] = true
    synthResult.mapFromTargetToCorpusPoints[targetIndex] = corpusIndex
    return nothing
end

#=
Does the point have a value?
I.E. has been synthesized already?
Getter method of the SynthResult object.

All points in target outside the synth region have a value
(from initial image, never mutated.)
Otherwise only some points in the synth region of the target have a value
(have been synthesized i.e. mutated as search progresses and improves.)
This is only pertinent for the first pass, since the algorithm starts
by nullify or blacken the values of point in the synth region.
After the first pass, they have a value, albeit a crude value.
=#
function doesPointHaveValue(synthResult, patchPoint )::Bool
    # assert patchPoint and synthResult in same frame, the target
    return synthResult.hasValue[patchPoint]
end

#=
Is the targetPoint in the synth region (not the context)
AND has been synthesized (has a best match in the corpus.)

In original C code: has_source_neighbor()

!!! This is not the original implementation.
The original relies on sourceOf field of Neighbor,
i.e. on caching SynthResult.mapFromTargetToCorpusPoints[targetIndex] in the Neighbor.
Which probably gives better cpu cache coherence.
Since the implementation below:
1) accesses the original mask
2) accesses an inverted copy of the mask (SynthResult.hasValue)
both of which are large and both might not fit in the cpu cache.
=#
function isInSynthAndWasSynthesized(
        synthResult,
        targetMaskedImage, patchPoint   # framed point ::PointInMaskedImage)
        )::Bool
        # assert patchPoint is in target frame

    #=
    Order of the terms is important for performance:
    hasValue is  general, for any target point, and is usually true when target is much larger than synth region.
    =#
    return ( isPointSelectedInMaskedImage(targetMaskedImage, patchPoint)  # is in synth region
            && doesPointHaveValue(synthResult, patchPoint)
            )
end

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
# For generalized eigenvalue problems, and application of linear
# operators, we need to also apply the metric (or its inverse) to the
# coefficients, to stay in the correct space.
operator_metric(B) = B'B

# However, for uniform finite-differences, the coefficients coincide
# with the function value at the nodes, hence we should /not/ apply
# the metric.
operator_metric(::Uniform,     ::FiniteDifferencesOrRestricted) = I
# For non-uniform finite-differences, we should, however.
operator_metric(::NonUniform, B::FiniteDifferencesOrRestricted) = I

operator_metric(B::FiniteDifferencesOrRestricted) =
    operator_metric(distribution(B), B)

operator_metric(B::FEDVROrRestricted) = I

# * Linear operator
"""
    LinearOperator(A, B⁻¹, temp)

A helper object used to apply the action of a linear operator, whose
matrix representation is given by `A`, in the basis whose metric
inverse is `B⁻¹`. For orthogonal bases (such as
e.g. [`FiniteDifferences`](@ref) and [`FEDVR`](@ref)), only
multiplication by `A` is necessary, but non-orthonal bases (such as
[`BSpline`](@ref)) necessitates the application of the metric inverse as well.
"""
struct LinearOperator{TA,TB,TT}
    A::TA
    B⁻¹::TB
    temp::TT
end

Base.size(S::LinearOperator, args...) = size(S.A, args...)
Base.eltype(S::LinearOperator) = eltype(S.A)

function LinearAlgebra.mul!(y, M::LinearOperator, x,
                            α::Number=true, β::Number=false)
    mul!(M.temp, M.A, x, α, false)
    if iszero(β)
        ldiv!(y, M.B⁻¹, M.temp)
    else
        y .*= β
        ldiv!(M.B⁻¹, M.temp)
        y .+= M.temp
    end
end

LinearAlgebra.mul!(y, M::LinearOperator{<:Any,Nothing}, x,
                   α::Number=true, β::Number=false) =
    mul!(y, M.A, x, α, β)

Base.:(*)(M::LinearOperator, x) =
    LinearAlgebra.mul!(similar(x), M, x)

LinearOperator(A, B) =
    LinearOperator(A, factorize(B), Vector{eltype(A)}(undef, size(A,1)))

LinearOperator(A, ::UniformScaling) =
    LinearOperator(A, nothing, nothing)

"""
    LinearOperator(A, R)

Construct the linear operator whose matrix representation in the basis
`R` is `A`, automatically deducing if application of the metric
inverse is also necessary.
"""
LinearOperator(A, R::BasisOrRestricted) =
    LinearOperator(A, operator_metric(R))

# * Diagonal operator

"""
    DiagonalOperator

Helper structure that when acting on a function ``f(x)`` produces the
product ``f(x)g(x)``, where ``g(x)`` is e.g. a potential. This is
similar to `R'*QuasiDiagonal(g.(x))*R` in intent, but simplifies the
case when ``g(x)`` needs to be updated, without computing a lot
overlap integrals. This is accomplished by computing the
[`FunctionProduct`](@ref) of the two functions directly.
"""
struct DiagonalOperator{Diag,FunctionProduct,B}
    diag::Diag
    ρ::FunctionProduct
    R::B
end

"""
    DiagonalOperator(f)

Construct [`DiagonalOperator`](@ref) from a function expansion `f`.
"""
function DiagonalOperator(f; kwargs...)
    T = eltype(f)
    R,c = f.args

    DiagonalOperator(c, FunctionProduct{false}(R, R, T; kwargs...), R)
end

Base.size(S::DiagonalOperator) = (length(S.ρ.ρ),length(S.ρ.ρ))
Base.size(S::DiagonalOperator, i) = size(S)[i]
Base.eltype(S::DiagonalOperator) = eltype(S.ρ)

function copyto_if_different!(dst, src)
    dst === src && return dst
    copyto!(dst, src)
end

"""
    copyto!(o::DiagonalOperator, diag::AbstractVector)

Update the [`DiagonalOperator`](@ref) to represent multiplication by
the function whose expansion coefficients are `diag`.
"""
function Base.copyto!(o::DiagonalOperator, diag::AbstractVector)
    copyto_if_different!(o.diag, diag)
    o
end

"""
    mul!(y, L::DiagonalOperator, x[, α=1, β=0])

Compute the action of the [`DiagonalOperator`](@ref) `L` on `x` and
store the result in `y`.
"""
function LinearAlgebra.mul!(y, L::DiagonalOperator, x,
                            α::Number=true, β::Number=false)
    copyto!(L.ρ, L.diag, x)
    if iszero(β)
        y .= false
    elseif !isone(β)
        lmul!(β, y)
    end
    y .+= α .* L.ρ.ρ
end

Base.:(*)(L::DiagonalOperator, x) =
    LinearAlgebra.mul!(similar(x), L, x)

function Base.copyto!(M::AbstractMatrix, L::DiagonalOperator,
                      x = L.R \ (r -> one(eltype(L.R))).(axes(L.R,1)))
    copyto!(L.ρ, L.diag, x)
    copyto!(M, L.ρ)
    M
end

Base.similar(::Type{<:Diagonal}, L::DiagonalOperator) =
    Diagonal(Vector{eltype(L)}(undef, size(L.R,2)))

Base.similar(::Type{<:BandedMatrix}, L::DiagonalOperator) =
    Matrix(undef, L.R, eltype(L))

Base.similar(L::DiagonalOperator) =
    similar(metric_shape(L.R), L)

"""
    Matrix(L::DiagonalOperator[, x])

Compute the matrix representation of the [`DiagonalOperator`](@ref)
`L`, optionally weighted by the function whose expansion coefficients
are given by `x`.

Note that this is allocating and potentially inefficient (especially
in the case of non-orthogonal bases), and should thus not be used in
compute-intense code. Furthermore, if you have a Julia function `f`
corresponding to the diagonal operator, it is most likely more
reasonable to take the `R'*QuasiDiagonal(f.(x))*R` approach instead.
"""
Matrix(L::DiagonalOperator, args...) =
    copyto!(similar(L), L, args...)

# * Shift-and-invert

"""
    ShiftAndInvert(A⁻¹, B, temp)

Represents a shifted-and-inverted operator, suitable for Krylov
iterations when targetting an interior point of the eigenspectrum.
`A⁻¹` is a factorization of the shifted operator and `B` is the metric
matrix.
"""
struct ShiftAndInvert{TA,TB,TT}
    A⁻¹::TA
    B::TB
    temp::TT

@doc raw"""
    ShiftAndInvert(A, B[, σ=0])

Construct the shifted-and-inverted operator corresponding to the
eigenproblem ``(\mat{A}-σ\mat{B})\vec{x} = \lambda\mat{B}\vec{x}``.
"""
    function ShiftAndInvert(A, B, σ::Number=0)
        A⁻¹ = factorize(A-σ*B)
        temp = Vector{eltype(A)}(undef, size(A,1))

        new{typeof(A⁻¹),typeof(B),typeof(temp)}(A⁻¹,B,temp)
    end

@doc raw"""
    ShiftAndInvert(A, ::UniformScaling[, σ=0])

Construct the shifted-and-inverted operator corresponding to the
eigenproblem ``(\mat{A}-σ\mat{I})\vec{x} = \lambda\vec{x}``.
"""
    function ShiftAndInvert(A, ::UniformScaling, σ::Number=0)
        A⁻¹ = factorize(A-σ*I)
        new{typeof(A⁻¹), Nothing, Nothing}(A⁻¹, nothing, nothing)
    end
end

@doc raw"""
    ShiftAndInvert(A, R[, σ=0])

Construct the shifted-and-inverted operator corresponding to the
eigenproblem ``(\mat{A}-σ\mat{B})\vec{x} = \lambda\mat{B}\vec{x}``,
where ``\mat{B}`` is the suitable operator metric, depending on the
`AbstractQuasiMatrix` `R` employed.
"""
ShiftAndInvert(A, R::BasisOrRestricted, σ::Number=0) =
    ShiftAndInvert(A, operator_metric(R), σ)

Base.size(S::ShiftAndInvert, args...) = size(S.A⁻¹, args...)
Base.eltype(S::ShiftAndInvert) = eltype(S.A⁻¹)

function LinearAlgebra.mul!(y,M::ShiftAndInvert,x)
    mul!(M.temp, M.B, x)
    ldiv!(y, M.A⁻¹, M.temp)
end

LinearAlgebra.mul!(y,M::ShiftAndInvert{<:Any,Nothing},x) =
    ldiv!(y, M.A⁻¹, x)


text:
import Base: size, axpy!
import Base.LinAlg: scale!
import Base.BLAS: gemm!
import Base: shmem_rand, shmem_randn

export ShareGLRM, share

### GLRM TYPE
type ShareGLRM{L<:Loss, R<:Regularizer}<:AbstractGLRM
    A::SharedArray               # The data table transformed into a coded array
    losses::Array{L,1}           # array of loss functions
    rx::Regularizer              # The regularization to be applied to each row of Xᵀ (column of X)
    ry::Array{R,1}               # Array of regularizers to be applied to each column of Y
    k::Int                       # Desired rank
    observed_features::ObsArray  # for each example, an array telling which features were observed
    observed_examples::ObsArray  # for each feature, an array telling in which examples the feature was observed
    X::SharedArray{Float64,2}    # Representation of data in low-rank space. A ≈ X'Y
    Y::SharedArray{Float64,2}    # Representation of features in low-rank space. A ≈ X'Y
end

function share(glrm::GLRM)
    isa(glrm.A, SharedArray) ? A = glrm.A : A = convert(SharedArray,glrm.A)
    isa(glrm.X, SharedArray) ? X = glrm.X : X = convert(SharedArray, glrm.X)
    isa(glrm.Y, SharedArray) ? Y = glrm.Y : Y = convert(SharedArray, glrm.Y)
    return ShareGLRM(A, glrm.losses, glrm.rx, glrm.ry, glrm.k,
                     glrm.observed_features, glrm.observed_examples,
                     X, Y)
end

### todo: define objective for shared arrays so it's evaluated (safely) in parallel


text:
using MetaGraphs
using LightGraphs
using JSON

"""
Creates a node in the graph, given its properties.

# Arguments
- `id` :  stands for the key name of the index of the graph
"""
function createNode!(g::MetaDiGraph, properties::Dict=Dict(),
    id=nothing)
    # First, we add the vertex to the graph
    if add_vertex!(g) & !isempty(properties)
        # The added index vertex must be the last one, i.e.
        # the number of vertices in the graph
        vertexIdx = nv(g)

        # Add the key if exists
        if (id ∈ keys(properties))
            set_prop!(g, vertexIdx, id, string(properties[id]))
            delete!(properties, id)
        end
        # all the other properties
        set_props!(g, vertexIdx, properties)
    end
end

"""
Updates a Node properties accessing via:
+ idxNode
+ idKey + idValue
"""
function updateNode!(g::MetaDiGraph, newProperties;
    idxNode=nothing, idKey=nothing, idValue=nothing, idProperty::Symbol=:id)
    if (idKey !== nothing) & (idValue !== nothing)
        idxNode = g[idValue, idKey]
    end
    if idxNode !== nothing
        currentProperties = props(g, idxNode)
        newProperties = merge(newProperties,
            currentProperties)
        delete!(newProperties, idProperty)
        set_props!(g, idxNode, newProperties)
    end
end

"""
Creates a edge given a graph g and the nodes u and v,
both expressed by one of the following ways:
+ idxNode
+ idKey + idValue
"""
function createEdge!(g::MetaDiGraph; properties::Dict = Dict(),
    idxNodeU=nothing, idKeyU=nothing, idValueU=nothing,
    idxNodeV=nothing, idKeyV=nothing, idValueV=nothing,
    undirected::Bool=false)
    # get the indices on the graph
    if (idKeyU !== nothing) & (idValueU !== nothing)
        idxNodeU = g[idValueU, idKeyU]
    end
    if (idKeyV !== nothing) & (idValueV !== nothing)
        idxNodeV = g[idValueV, idKeyV]
    end
    # check if they exist, and if they don't, create them # TO-DO
    if (idxNodeU !== nothing) & (idxNodeV !== nothing)
        add_edge!(g, idxNodeU, idxNodeV)
        set_props!(g, idxNodeU, idxNodeV, properties)
    end

    if undirected
        createEdge!(g; properties = properties,
            idxNodeU = idxNodeV, idKeyU = idKeyV, idValueU = idValueV,
            idxNodeV = idxNodeU, idKeyV = idKeyU, idValueV = idValueU)
    end
end

function jolt_parse(g::MetaDiGraph, splitted::Vector)
    idx1, idx2 = 1, 2
    for line in splitted
        if line != ""
            parsed = JSON.parse(line)
            try
            if "data" in keys(parsed) # other will be headers and commits
                response = parsed["data"]
                # the response will be a bunch of entities
                for entityDict in response
                    # each entity can be "()", "->" or "<-"
                    entityType = collect(keys(entityDict))[1]
                    entity = entityDict[entityType]

                    if entityType == "()"
                        # {
                        #   "()":
                        #       [ node_id, [ node_labels],
                        #       {"prop1": "value1", "prop2": "value2"}]
                        # }

                        # exists?
                        try
                            idxNode = g[string(entity[1]), :id]
                            properties = Dict([(Symbol(x), y) for (x, y) ∈ entity[3]])
                            updateNode!(g, properties; idxNode = idxNode)
                        catch err
                            # does not exist yet
                            delete!(entity[3], "id") # just to make sure there won't be problems
                            properties = merge(Dict(:type => string(entity[2][1])),
                                entity[3], Dict(:id => string(entity[1])))
                            properties = Dict([(Symbol(x), y) for (x, y) ∈ properties])
                            createNode!(g, properties, :id)
                        end

                    elseif entityType == "->"
                        #  {"->":
                            # [
                                    # rel_id, start_node_id,
                                    # rel_type, end_node_id,
                                    # {properties}
                            # ]}
                        delete!(entity[5], "id") # just to make sure there won't be problems

                        idx1 = try
                            g[string(entity[2]), :id]
                        catch err
                            createNode!(g, Dict(:id => entity[2]), :id)
                            nv(g)
                        end

                        idx2 = try
                            g[string(entity[4]), :id]
                        catch err
                            createNode!(g, Dict(:id => entity[4]), :id)
                            nv(g)
                        end

                        createEdge!(g;
                            idxNodeU = idx1,
                            idxNodeV = idx2,
                            properties = Dict(
                                :type => string(entity[3]),
                                :rel_id => string(entity[1]),
                                [Symbol(k) =>  v for (k, v) in entity[5]]...
                        ))

                    elseif entityType == "<-"
                        #  {"->":
                            # [
                                    # rel_id, end_node_id,
                                    # rel_type, start_node_id,
                                    # {properties}
                            # ]}
                        delete!(entity[5], "id") # just to make sure there won't be problems

                        idx2 = try
                            g[string(entity[2]), :id]
                        catch err
                            createNode!(g, Dict(:id => entity[2]), :id)
                            nv(g)
                        end

                        idx1 = try
                            g[string(entity[4]), :id]
                        catch err
                            createNode!(g, Dict(:id => entity[4]), :id)
                            nv(g)
                        end

                        createEdge!(g;
                            idxNodeU = idx1,
                            idxNodeV = idx2,
                            properties = Dict(
                                :type => string(entity[3]),
                                :rel_id => string(entity[1]),
                                [Symbol(k) =>  v for (k, v) in entity[5]]...
                        ))
                    end
                end
            end
        catch err
        end
        end
    end
    return g
end


text:
# Tuple constructor that either keeps
# it as a Tuple or turns it into a Tuple.
tuple_convert(t::Tuple) = t
tuple_convert(x) = tuple(x)
tuple_convert(x1, x2, xs...) = (x1, x2, xs...)


text:
###############################################################################
"""
    plt = plotsfgmax(fg::VisClaw.FixedGrid, fgmax::VisClaw.FGmax, var::Symbol; kwargs...)

    plotsfgmax!(plt::Plots.Plot, fg::VisClaw.FixedGrid, fgmax::VisClaw.FGmax, var::Symbol; kwargs...)

"""
function plotsfgmax!(plt, fg::VisClaw.FixedGrid, fgmax::VisClaw.FGmax, var::Symbol=:D; kwargs...)
    # parse keyword args
    kwdict = KWARG(kwargs)

	# get var
	val = copy(getfield(fgmax, var))
	# check
	isempty(val) && error("Empty: $var")

    # vector
	if fg.style == 0 || fg.style == 1 || fg.style == 3

		# correct
	    (var==:D) && (val = val + fgmax.topo)
		(var==:Dmin) && (val = val - fgmax.topo)

		# plot
	    plt = Plots.plot!(plt, fg.x, fg.y, val; kwdict...)

    elseif fg.style == 2
		seriestype, kwdict = VisClaw.kwarg_default(kwdict, VisClaw.parse_seriestype, :heatmap)

        x = collect(Float64, LinRange(fg.xlims[1], fg.xlims[end], fg.nx))
        y = collect(Float64, LinRange(fg.ylims[1], fg.ylims[end], fg.ny))

		# wet cells
		wet = fgmax.D .!= 0.0
		land = fgmax.topo .> 0.0

	    # correct
	    (var==:D) && (val[wet] = val[wet] + fgmax.topo[wet])
		(var==:Dmin) && (val[wet] = val[wet] - fgmax.topo[wet])
		(var==:Dmin) && (val[land] .= NaN)
		isa(val, AbstractArray{Dates.DateTime}) || (val[.!wet] .= NaN)
	    #val[.!wet] .= NaN

	    # plot
	    plt = Plots.plot!(plt, x, y, val; ratio=:equal, xlims=fg.xlims, ylims=fg.ylims, seriestype=seriestype, kwdict...)

	elseif fg.style == 4
		seriestype, kwdict = VisClaw.kwarg_default(kwdict, VisClaw.parse_seriestype, :heatmap)

		x = collect(Float64, LinRange(fg.xlims[1], fg.xlims[end], fg.nx))
        y = collect(Float64, LinRange(fg.ylims[1], fg.ylims[end], fg.ny))

		# correct
	    (var==:D) && (val = val + fgmax.topo)
		(var==:Dmin) && (val = val - fgmax.topo)

		var_vec = copy(val)
		if !isa(val, AbstractArray{Dates.DateTime})
		    val = NaN*zeros(Float64 ,(fg.ny, fg.nx))
		    val[fg.flag] = var_vec
		end
		val = reverse(val, dims=1)

		# plot
	    plt = Plots.plot!(plt, x, y, val; ratio=:equal, xlims=fg.xlims, ylims=fg.ylims, seriestype=seriestype, kwdict...)
    end

    # return
    return plt
end
###############################################################################
"""
$(@doc plotsfgmax!)
"""
plotsfgmax(fg::VisClaw.FixedGrid, fgmax::VisClaw.FGmax, var::Symbol=:D; kwargs...) =
plotsfgmax!(Plots.plot(), fg, fgmax, var; kwargs...)
###############################################################################


text:
# NUM_MAX_ITERATIONS=100000
NUM_MAX_ITERATIONS=10000

function get_pixel_color(fracgen::FractalGenerator,px::Int64,py::Int64)
   r=0.0
   g=0.0
   b=0.0

   for i in 1:fracgen.imagespec.num_supersamples
      dx=Float64(i-1)/fracgen.imagespec.num_supersamples
      x=((Float64(px)+dx)/fracgen.imagespec.width-0.5)*fracgen.cameraspec.zoom+
         fracgen.cameraspec.center.re

      for j in 1:fracgen.imagespec.num_supersamples
         dy=Float64(j-1)/fracgen.imagespec.num_supersamples
         y=((Float64(py)+dy)/fracgen.imagespec.height-0.5)*fracgen.cameraspec.zoom+
            fracgen.cameraspec.center.im
         prop=get_num_iters(fracgen.fracspec.fracfunc_closure(Complex(x,y))...)
         if fracgen.cameraspec.affine
            prop=max(0.0,prop)
         else
            prop=max(0.0,prop/NUM_MAX_ITERATIONS)
         end
         prop=prop^fracgen.cameraspec.pow
         color=blend(fracgen.cameraspec.gradient,prop)
         r+=color.r
         g+=color.g
         b+=color.b
      end
   end

   nss=fracgen.imagespec.num_supersamples*fracgen.imagespec.num_supersamples
   r/=nss
   g/=nss
   b/=nss
   return RGB8bit(r,g,b)
end

function get_num_iters(f::Function,zinit::Complex=zero(Complex))
   i=0
   THRESHOLD=2^16
   z=zinit
   while i < NUM_MAX_ITERATIONS
      # tempreal=zx
      # zx=zx*zx-zy*zy-cx
      # zy=2*tempreal*zy-cy
      z=f(z)
      if real(z*z') >= THRESHOLD
         break
      end
      i+=1
   end
   if i==NUM_MAX_ITERATIONS
      return 0.0
   end
   i=Float64(i)
   i+=(1.0-log2(log2(real(z*z'))/2.0))
   return i
end
# get_num_iters(x::Float64,y::Float64)=julia(x,y,0.0,0.0)

# inverse plane
# 1/(a+b*i) = (a-b*i)/(a^2+b^2)
# @inline function get_num_iters(c::Complex)
#    # xysqrd=x*x+y*y
#    julia(zero(Complex),1.0/c-0.25)
# end
# @inline function get_num_iters(c::Complex)
#    get_num_iters_inverse(c+0.25)
# end
# function julia(z::Complex,c::Complex)
#    i=0
#    THRESHOLD=2^16
#    while i < NUM_MAX_ITERATIONS
#       # tempreal=zx
#       # zx=zx*zx-zy*zy-cx
#       # zy=2*tempreal*zy-cy
#       z=z*z-c
#       if real(z*z') >= THRESHOLD
#          break
#       end
#       i+=1
#    end
#    if i==NUM_MAX_ITERATIONS
#       return 0.0
#    end
#    i=Float64(i)
#    i+=(1.0-log2(log2(real(z*z'))/2.0))
#    return i
# end
function blend(gradient::Gradient,prop::Float64)
   if prop == 0.0
      return gradient.nilcolor
   end
   return blend(gradient.startcolor,gradient.endcolor,prop)
end
function blend(startcolor::RGB8bit,endcolor::RGB8bit,prop::Float64)
   bl=(s,e,p)->p*s+e*(1.0-p) # blend between s and e
   function blnd(s,e,p)
      # have color values over 1.0 wrap around
      o=bl(s,e,p)%1.0
      return o > 0 ? o : 1.0+o
   end
   return RGB8bit(
      blnd(startcolor.r,endcolor.r,prop),
      blnd(startcolor.g,endcolor.g,prop),
      blnd(startcolor.b,endcolor.b,prop))
end

max_stars_repo_name: jblondin/Fractal

max_issues_repo_name: jblondin/Fractal

max_forks_repo_name: jblondin/Fractal


text:
using Impute
using Test
using DataFrames
using RDatasets
using Statistics

@testset "Impute" begin
    a = Vector{Union{Float64, Missing}}(1.0:1.0:20.0)
    a[[2, 3, 7]] .= missing
    mask = map(!ismissing, a)

    @testset "Drop" begin
        result = impute(a, :drop; limit=0.2)
        expected = copy(a)
        deleteat!(expected, [2, 3, 7])

        @test result == expected
    end

    @testset "Interpolate" begin
        result = impute(a, :interp; limit=0.2)
        @test result == collect(1.0:1.0:20)
        @test result == interp(a)
    end

    @testset "Fill" begin
        @testset "Value" begin
            fill_val = -1.0
            result = impute(a, :fill, fill_val; limit=0.2)
            expected = copy(a)
            expected[[2, 3, 7]] .= fill_val

            @test result == expected
        end

        @testset "Mean" begin
            result = impute(a, :fill; limit=0.2)
            expected = copy(a)
            expected[[2, 3, 7]] .= mean(a[mask])

            @test result == expected
        end
    end

    @testset "LOCF" begin
        result = impute(a, :locf; limit=0.2)
        expected = copy(a)
        expected[2] = 1.0
        expected[3] = 1.0
        expected[7] = 6.0

        @test result == expected
    end

    @testset "NOCB" begin
        result = impute(a, :nocb; limit=0.2)
        expected = copy(a)
        expected[2] = 4.0
        expected[3] = 4.0
        expected[7] = 8.0

        @test result == expected
    end

    @testset "DataFrame" begin
        data = dataset("boot", "neuro")
        df = impute(data, :interp; limit=1.0)
    end

    @testset "Matrix" begin
        data = Matrix(dataset("boot", "neuro"))

        @testset "Drop" begin
            result = Iterators.drop(data)
            @test size(result, 1) == 4
        end

        @testset "Fill" begin
            result = impute(data, :fill, 0.0; limit=1.0)
            @test size(result) == size(data)
        end
    end

    @testset "Not enough data" begin
        @test_throws ImputeError impute(a, :drop)
    end

    @testset "Chain" begin
        data = Matrix(dataset("boot", "neuro"))
        result = chain(
            data,
            Impute.Interpolate(),
            Impute.LOCF(),
            Impute.NOCB();
            limit=1.0
        )

        @test size(result) == size(data)
        # Confirm that we don't have any more missing values
        @test !any(ismissing, result)
    end

    @testset "Alternate missing functions" begin
        data1 = dataset("boot", "neuro")                # Missing values with `missing`
        data2 = impute(data1, :fill, NaN; limit=1.0)     # Missing values with `NaN`

        @test impute(data1, :drop; limit=1.0) == dropmissing(data1)

        result1 = chain(data1, Impute.Interpolate(), Impute.Drop(); limit=1.0)
        result2 = chain(data2, isnan, Impute.Interpolate(), Impute.Drop(); limit=1.0)
        @test result1 == result2
    end
end


text:
MCMC_Driver.getchi2forall(opt)
savefig("convergence.png", dpi=300)

irange=1:10:101
Tools3D.savepngs(opt, irange, dz=dz, extendfrac=extendfrac, slicesx=slicesx, slicesy=slicesy, slicesz=slicesz)

M  = Tools3D.assembleTat1(opt;burninfrac=0.33)
Mmean = mean(M, dims=2)
Msd = sqrt.(var(M, dims=2))

Tools3D.slicemodel(vec(Mmean),0,[0. 0.],[0.], opt, slicesx=slicesx, slicesy=slicesy, slicesz=slicesz, dz=dz, extendfrac=extendfrac)
savefig("meanmodel.png", dpi=300)
Tools3D.slicemodel(vec(Msd),0,[0. 0.],[0.], opt, slicesx=slicesx, slicesy=slicesy, slicesz=slicesz, dz=dz, extendfrac=extendfrac)
savefig("sdevmodel.png")


text:
using WaveformCommunications
using Test

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


text:
using Unitful: Unitful, Units, Quantity, @u_str
export @u_str

unitfy(::Nothing, u) = nothing
unitfy(::Nothing, ::Nothing) = nothing
unitfy(::Missing, u) = missing
unitfy(::Missing, ::Nothing) = missing
unitfy(v, ::Nothing) = v
unitfy(v::Number, u::Units) = Quantity(v, u)
unitfy(v::Array, u::Units) = unitfy.(v, u)
unitfy(v::Tuple, u::Units) = unitfy.(v, u)
unitfy(v::Quantity, u::Units) = Unitful.uconvert(u, v)
unitfy(v::Array{<:Union{Quantity,Missing}}, u::Units) = unitfy.(v, u)
unitfy(v::Tuple{Vararg{Union{Quantity,Missing}}}, u::Units) = unitfy.(v, u)
unitfy(v, u) = u(v)
unitfy(v::V, ::Type{V}) where V = v

deunitfy(v) = v
deunitfy(v::Quantity) = Unitful.ustrip(v)
deunitfy(v::Array) = deunitfy.(v)
deunitfy(v::Tuple) = deunitfy.(v)
deunitfy(v, u) = deunitfy(unitfy(v, u))

promoteunit(u...) = Unitful.promote_unit(filter(!isnothing, u)...)
promoteunit(::Nothing) = nothing
promoteunit() = nothing

hasunit(v::Units) = !Unitful.isunitless(v)
hasunit(::Nothing) = false

using DataFrames: DataFrame, DataFrames
for f in (:unitfy, :deunitfy)
    @eval $f(df::DataFrame, U::Vector) = begin
        r = DataFrame()
        for (n, c, u) in zip(propertynames(df), eachcol(df), U)
            r[!, n] = $f.(c, u)
        end
        r
    end
end

import Dates
unitfy(df::DataFrame; kw...) = begin
    #HACK: default constructor for common types to avoid scope binding issue
    D = merge(Dict(
        :Date => Dates.Date,
    ), Dict(kw))
    p = r"(.+)\(([^\(\)]+)\)$"
    M = match.(p, names(df))
    n(m::RegexMatch) = m.match => strip(m.captures[1])
    n(m) = nothing
    N = filter(!isnothing, n.(M))
    isempty(N) && return df
    u(m::RegexMatch) = begin
        s = m.captures[2]
        #HACK: assume type constructor if the label starts with `:`
        e = startswith(s, ":") ? Symbol(s[2:end]) : :(@u_str($s))
        #HACK: use Main scope for type constructor evaluation
        #TODO: remove fallback eval in favor of explict constructor mapping
        haskey(D, e) ? D[e] : Main.eval(e)
    end
    u(m) = nothing
    U = u.(M)
    DataFrames.rename(unitfy(df, U), N...)
end
deunitfy(df::DataFrame) = deunitfy(df, repeat([nothing], DataFrames.ncol(df)))

export unitfy, deunitfy


text:
info("Testing ", MOD.RefVector)
for T in IntegerTypes
    ref = T[1; 1; 2; 3; 3]
    k = convert(T, 3)
    y = MOD.RefVector(ref,k)

    @test all(y.ref .== ref)
    @test y.k == k

    @test_throws ErrorException MOD.RefVector(ref[ref .!= 2], k)
    @test_throws ErrorException MOD.RefVector(ref[ref .!= 2], k + one(T))
    @test_throws ErrorException MOD.RefVector(ref[ref .!= 2], k - one(T))
end

# Class Functions

n_k = [80; 100; 120]
k = length(n_k)
n = sum(n_k)
p = 3

y, X = sampledata(n_k, p)
Z = vcat([sum(X[y.ref .== i,:],1) for i = 1:k]...)
M = Z ./ n_k
H = X .- M[y, :]
Σ = H'H/(n-1)

#=
y = MOD.RefVector(vcat([Int64[i for j = 1:n_k[i]] for i = 1:k]...), k)
X = vcat([rand(n_k[i], p) .+ (10rand(1,p) .- 5) for i = 1:k]...)
σ = sortperm(rand(sum(n_k)))
y = MOD.RefVector(y[σ])
X = X[σ,:]
=#

info("Testing ", MOD.classcounts)
for U in IntegerTypes
    @test all(n_k .== MOD.classcounts(convert(MOD.RefVector{U}, y)))
end

info("Testing ", MOD.classtotals)
for T in FloatingPointTypes
    for U in IntegerTypes
        X_tmp = convert(Array{T}, X)
        y_tmp = convert(MOD.RefVector{U}, y)
        Z_tmp = convert(Array{T}, Z)
        @test_approx_eq MOD.classtotals(Val{:row}, X_tmp, y_tmp)  Z_tmp
        @test_approx_eq MOD.classtotals(Val{:col}, X_tmp', y_tmp) Z_tmp'
        @test_throws DimensionMismatch MOD.classtotals(Val{:row}, X_tmp', y_tmp)
        @test_throws DimensionMismatch MOD.classtotals(Val{:col}, X_tmp, y_tmp)
    end
end

info("Testing ", MOD.classmeans)
for T in FloatingPointTypes
    for U in IntegerTypes
        X_tmp = convert(Array{T}, X)
        y_tmp = convert(MOD.RefVector{U}, y)
        M_tmp = convert(Array{T}, M)
        @test_approx_eq MOD.classmeans(Val{:row}, X_tmp,  y_tmp) M
        @test_approx_eq MOD.classmeans(Val{:col}, X_tmp', y_tmp) M'
    end
end

info("Testing ", MOD.centerclasses!)
for T in FloatingPointTypes
    for U in IntegerTypes
        X_tmp = copy(convert(Array{T}, X))
        y_tmp = convert(MOD.RefVector{U}, y)
        M_tmp = convert(Array{T}, M)
        H_tmp = convert(Array{T}, H)
        @test_approx_eq MOD.centerclasses!(Val{:row}, copy(X_tmp), M_tmp, y_tmp) H
        @test_approx_eq MOD.centerclasses!(Val{:col}, copy(X_tmp)', M_tmp', y_tmp) H'

        @test_throws DimensionMismatch MOD.centerclasses!(Val{:row}, X_tmp, Array(T,k,p+1), y_tmp)
        @test_throws DimensionMismatch MOD.centerclasses!(Val{:col}, X_tmp', Array(T,p+1,k), y_tmp)

        @test_throws DimensionMismatch MOD.centerclasses!(Val{:row}, X_tmp, Array(T,k+1,p), y_tmp)
        @test_throws DimensionMismatch MOD.centerclasses!(Val{:col}, X_tmp', Array(T,p,k+1), y_tmp)
    end
end

info("Testing ", MOD.regularize!)
for T in FloatingPointTypes
    S1 = T[1 2 3;
           4 5 6;
           7 8 9]

    n = size(S1,1)

    S2 = T[7 8 9;
           8 3 6;
           1 7 4]

    B = T[1 2;
          3 4;
          5 6]

    b = T[1;
          3]

    @test_approx_eq MOD.regularize!(copy(S1), zero(T),  S2) S1
    @test_approx_eq MOD.regularize!(copy(S1), one(T),   S2) S2
    @test_approx_eq MOD.regularize!(copy(S1), one(T)/2, S2) (1-one(T)/2)*S1 + (one(T)/2)*S2
    @test_throws ErrorException MOD.regularize!(copy(S1), -one(T),  S2)
    @test_throws ErrorException MOD.regularize!(copy(S1), 2*one(T), S2)
    @test_throws DimensionMismatch MOD.regularize!(S1, one(T), B)

    @test_approx_eq MOD.regularize!(copy(S1), zero(T)) S1
    @test_approx_eq MOD.regularize!(copy(S1), one(T)) (trace(S1)/n)*eye(n)
    @test_approx_eq MOD.regularize!(copy(S1), one(T)/2) (1-one(T)/2)*S1 + (one(T)/2)*(trace(S1)/n)*eye(n)
end

info("Testing ", MOD.dotvectors)
for T in FloatingPointTypes
    X_tmp = convert(Array{T}, X)

    @test_approx_eq MOD.dotvectors(Val{:row}, X_tmp)  sum(X_tmp .* X_tmp,2)
    @test_approx_eq MOD.dotvectors(Val{:col}, X_tmp') sum(X_tmp .* X_tmp,2)

    @test_throws DimensionMismatch MOD.dotvectors!(Val{:row}, X_tmp,  Array(T, n+1))
    @test_throws DimensionMismatch MOD.dotvectors!(Val{:col}, X_tmp', Array(T, n+1))
end

info("Testing ", MOD.whitendata_svd!)
for T in FloatingPointTypes
    H_tmp = copy(convert(Array{T}, H))
    Σ_tmp = copy(convert(Array{T}, Σ))

    # Test full rank case
    for λ in (zero(T), convert(T, 0.25), convert(T, 0.5), convert(T, 0.75), one(T))
        W_tmp = MOD.whitendata_svd!(copy(H_tmp), λ)
        S_tmp = (1-λ)*Σ_tmp + (λ*trace(Σ_tmp)/p)*I
        @test_approx_eq eye(T,3) W_tmp*S_tmp*(W_tmp')
    end

    # Test degenerate case
    H_tmp = eye(T,3) .- mean(eye(T,3))
    @test_throws ErrorException MOD.whitendata_svd!(copy(H_tmp), zero(T))
end

info("Testing ", MOD.whitendata_qr!)
for T in FloatingPointTypes
    H_tmp = copy(convert(Array{T}, H))
    Σ_tmp = copy(convert(Array{T}, Σ))

    # Test full rank case
    W_tmp = MOD.whitendata_qr!(copy(H_tmp))
    @test_approx_eq eye(T,3) (W_tmp')*Σ_tmp*W_tmp

    # Test degenerate case
    H_tmp = eye(T,3)
    @test_throws ErrorException MOD.whitendata_qr!(copy(H_tmp))

    H_tmp = ones(T,4,3)
    @test_throws ErrorException MOD.whitendata_qr!(copy(H_tmp))
end

info("Testing ", MOD.whitencov_chol!)
for T in FloatingPointTypes
    H_tmp = copy(convert(Array{T}, H))
    Σ_tmp = copy(convert(Array{T}, Σ))

    # Test unregularized, full rank case
    W_tmp = MOD.whitencov_chol!(copy(Σ_tmp), Nullable{T}())
    @test_approx_eq eye(T,3) (W_tmp')*Σ_tmp*W_tmp

    # Test regularized cases
    for λ in (zero(T), convert(T, 0.25), convert(T, 0.5), convert(T, 0.75), one(T))
        W_tmp = MOD.whitencov_chol!(copy(Σ_tmp), Nullable(λ))
        S_tmp = (1-λ)*Σ_tmp + (λ*trace(Σ_tmp)/p)*I
        @test_approx_eq eye(T,3) (W_tmp')*S_tmp*W_tmp
    end

    # Test degenerate case
    @test_throws ErrorException MOD.whitencov_chol!(diagm([one(T); zero(T)]), Nullable{T}())
    @test_throws ErrorException MOD.whitencov_chol!(diagm([one(T); zero(T)]), Nullable(zero(T)))
end


text:
import Base: ==

function ==(a::UnivalueHistory, b::UnivalueHistory)
    return (a.values == b.values) && (a.lastiter == b.lastiter) && (a.iterations == b.iterations)
end

function ValueHistories.History(DT::Type{V}, lastiter::I, iterations::Vector{I}, values::Vector{V}) where {I,V}
    h = History(DT)
    h.lastiter = lastiter
    h.iterations = iterations
    h.values = values
    return h
end

function ValueHistories.MVHistory(x::Dict{Symbol, UnivalueHistory})
    mvh = MVHistory()
    merge!(mvh.storage, x)
    return mvh
end

function Unmarshal.unmarshal(DT::Type{MVHistory}, parsedJson::AbstractDict, verbose::Bool = false, verboseLvl::Int = 0)
    mvh = Dict{Symbol, UnivalueHistory}()
    for (k,v) in parsedJson["mvhistory"]
        itype = eval(Symbol(parsedJson["mv_itypes"][k]))
        vtype = Symbol(parsedJson["mv_vtypes"][k])
        vtype = vtype == Symbol(OneClassActiveLearning.ConfusionMatrix) ? eval(OneClassActiveLearning.ConfusionMatrix) : eval(vtype)
        push!(mvh, Symbol(k) => History(vtype, v["lastiter"], convert(Vector{itype}, v["iterations"]), convert(Array{vtype}, v["values"])))
    end
    return MVHistory(mvh)
end

function JSON.json(x::MVHistory)
    itypes = Dict(k => typeof(x[k]).parameters[1] for (k,v) in x.storage)
    vtypes = Dict(k => typeof(x[k]).parameters[2] for (k,v) in x.storage)
    return JSON.json(Dict(:mvhistory => x.storage, :mv_itypes => itypes, :mv_vtypes => vtypes))
end

function JSON.lower(x::T) where T <: OneClassActiveLearning.SplitStrategy
    return JSON.lower(typeof(x))
end

function JSON.lower(x::T) where T <: MLKernels.Kernel
    return JSON.lower(sprint(print, x))
end

function JSON.lower(x::T) where T <: SVDD.InitializationStrategy
    return JSON.lower(sprint(print, x))
end

function JSON.lower(x::T) where T <: OneClassActiveLearning.Oracle
    return JSON.lower(typeof(x))
end

function JSON.json(res::Result)
    jsonstring = "{"
    jsonstring *= "\"id\":" * JSON.json(res.id) * ","
    jsonstring *= "\"experiment\":" * JSON.json(res.experiment) * ","
    jsonstring *= "\"al_history\":" * JSON.json(res.al_history) * ","
    jsonstring *= "\"worker_info\":" * JSON.json(res.worker_info) * ","
    jsonstring *= "\"data_stats\":" * JSON.json(res.data_stats) * ","
    jsonstring *= "\"al_summary\":" * JSON.json(res.al_summary) * ","
    jsonstring *= "\"status\":" * JSON.json(res.status)
    jsonstring *= "}"
    return jsonstring
end

function write_result_to_file(output_file, r::OneClassActiveLearning.Result)
    r_reparsed = JSON.parse(JSON.json(r))
    d = OrderedDict()
    for k in ["id", "experiment", "worker_info", "data_stats", "al_history", "al_summary", "status"]
        d[k] = r_reparsed[k]
    end
    open(output_file, "w") do f
        JSON.print(f, d, 2)
    end
    return nothing
end

function unmarshal_al_summary(parsedJson::AbstractDict)
    al_summary = Dict()
    for (k_metric, v_metric) in parsedJson
        d = Dict{Symbol, Any}()
        for (k_summary, v_summary) in v_metric
            cur_res = nothing
            try
                if typeof(v_summary) <: Vector
                    cur_res = unmarshal(Vector{Float64}, v_summary)
                else
                    cur_res = unmarshal(Float64, v_summary)
                end
            catch ArgumentError e
                cur_res = typeof(v_summary) <: Vector ? zeros(length(v_summary)) : 0.0
            end
            d[Symbol(k_summary)] = cur_res
        end
        push!(al_summary, Symbol(k_metric) => d)
    end
    return al_summary
end

function convert_to_symbol_keys(x)
    helper(x) = x
    helper(d::Dict) = Dict(Symbol(k) => helper(v) for (k, v) in d)
    return helper(x)
end

function Unmarshal.unmarshal(DT::Type{OneClassActiveLearning.Result}, parsedJson::AbstractDict)
    id = parsedJson["id"]
    al_history = Unmarshal.unmarshal(MVHistory, parsedJson["al_history"])
    experiment = convert_to_symbol_keys(parsedJson["experiment"])
    worker_info = convert_to_symbol_keys(parsedJson["worker_info"])
    data_stats = Unmarshal.unmarshal(DataStats, parsedJson["data_stats"])
    al_summary = unmarshal_al_summary(parsedJson["al_summary"])
    status = Dict(:exit_code => Symbol(parsedJson["status"]["exit_code"]))
    return OneClassActiveLearning.Result(id, experiment, worker_info, data_stats, al_history, al_summary, status)
end