text:
value(x) = x
cuify(x) = error("To use LinSolveGPUFactorize, you must do `using CuArrays`")

# Piracy, should get upstreamed
LinearAlgebra.ldiv!(Y::AbstractArray, A::AbstractArray, B::AbstractArray) = (copyto!(Y,B); ldiv!(A,Y))

function __init__()
  @require ApproxFun="28f2ccd6-bb30-5033-b560-165f7b14dc2f" begin
    eval_u0(u0::ApproxFun.Fun) = false
  end

  @require Distributions="31c24e10-a181-5473-b8eb-7969acd0382f" begin
    handle_distribution_u0(_u0::Distributions.Sampleable) = rand(_u0)
  end

  @require ForwardDiff="f6369f11-7733-5829-9624-2563aa707210" begin

    value(x::Type{ForwardDiff.Dual{T,V,N}}) where {T,V,N} = V
    value(x::ForwardDiff.Dual) = value(ForwardDiff.value(x))

    @inline fastpow(x::ForwardDiff.Dual, y::ForwardDiff.Dual) = x^y

    # Support adaptive with non-dual time
    @inline function ODE_DEFAULT_NORM(u::AbstractArray{<:ForwardDiff.Dual,N},t) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline function ODE_DEFAULT_NORM(u::Array{<:ForwardDiff.Dual,N},t) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline ODE_DEFAULT_NORM(u::ForwardDiff.Dual,t) = abs(value(u))

    # When time is dual, it shouldn't drop the duals for adaptivity
    @inline function ODE_DEFAULT_NORM(u::AbstractArray{<:ForwardDiff.Dual,N},t::ForwardDiff.Dual) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((x for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline function ODE_DEFAULT_NORM(u::Array{<:ForwardDiff.Dual,N},t::ForwardDiff.Dual) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((x for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline ODE_DEFAULT_NORM(u::ForwardDiff.Dual,t::ForwardDiff.Dual) = abs(u)

    # Type piracy. Should upstream
    Base.nextfloat(d::ForwardDiff.Dual{T,V,N}) where {T,V,N} = ForwardDiff.Dual{T}(nextfloat(d.value), d.partials)
    Base.prevfloat(d::ForwardDiff.Dual{T,V,N}) where {T,V,N} = ForwardDiff.Dual{T}(prevfloat(d.value), d.partials)

    struct DiffCache{T<:AbstractArray, S<:AbstractArray}
        du::T
        dual_du::S
    end

    function DiffCache(u::AbstractArray{T}, siz, ::Type{Val{chunk_size}}) where {T, chunk_size}
        DiffCache(u, zeros(ForwardDiff.Dual{nothing,T,chunk_size}, siz...))
    end

    dualcache(u::AbstractArray, N=Val{ForwardDiff.pickchunksize(length(u))}) = DiffCache(u, size(u), N)

    get_tmp(dc::DiffCache, u::AbstractArray{T}) where T<:ForwardDiff.Dual = reinterpret(T, dc.dual_du)
    get_tmp(dc::DiffCache, u::AbstractArray) = dc.du
  end

  @require Measurements="eff96d63-e80a-5855-80a2-b1b0885c5ab7" begin

    value(x::Type{Measurements.Measurement{T}}) where {T} = T
    value(x::Measurements.Measurement) = Measurements.value(x)

    # Support adaptive steps should be errorless
    @inline function ODE_DEFAULT_NORM(u::AbstractArray{<:Measurements.Measurement,N},t) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline function ODE_DEFAULT_NORM(u::Array{<:Measurements.Measurement,N},t) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline ODE_DEFAULT_NORM(u::Measurements.Measurement,t) = abs(Measurements.value(u))
  end

  @require MonteCarloMeasurements="0987c9cc-fe09-11e8-30f0-b96dd679fdca" begin

    value(x::Type{MonteCarloMeasurements.AbstractParticles{T,N}}) where {T,N} = T
    value(x::MonteCarloMeasurements.AbstractParticles) = mean(x)

    # Support adaptive steps should be errorless
    @inline function ODE_DEFAULT_NORM(u::AbstractArray{<:MonteCarloMeasurements.AbstractParticles,N},t) where {N}
      sqrt(mean(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))))
    end
    @inline function ODE_DEFAULT_NORM(u::AbstractArray{<:MonteCarloMeasurements.AbstractParticles,N},t::AbstractArray{<:MonteCarloMeasurements.AbstractParticles,N}) where {N}
      sqrt(mean(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(value.(t)))))
    end
    @inline ODE_DEFAULT_NORM(u::MonteCarloMeasurements.AbstractParticles,t) = abs(value(u))
  end

  @require Unitful="1986cc42-f94f-5a68-af5c-568840ba703d" begin
    # Support adaptive errors should be errorless for exponentiation
    value(x::Type{Unitful.AbstractQuantity{T,D,U}}) where {T,D,U} = T
    value(x::Unitful.AbstractQuantity) = x.val
    @inline function ODE_DEFAULT_NORM(u::AbstractArray{<:Unitful.AbstractQuantity,N},t) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline function ODE_DEFAULT_NORM(u::Array{<:Unitful.AbstractQuantity,N},t) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline ODE_DEFAULT_NORM(u::Unitful.AbstractQuantity,t) = abs(value(u))
    @inline UNITLESS_ABS2(x::Unitful.AbstractQuantity) = real(abs2(x)/oneunit(x)*oneunit(x))
  end

  @require Tracker="9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" begin
    value(x::Type{Tracker.TrackedReal{T}}) where T = T
    value(x::Type{Tracker.TrackedArray{T,N,A}}) where {T,N,A} = Array{T,N}
    value(x::Tracker.TrackedReal)  = x.data
    value(x::Tracker.TrackedArray) = x.data

    # Support adaptive with non-tracked time
    @inline function ODE_DEFAULT_NORM(u::Tracker.TrackedArray,t) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline function ODE_DEFAULT_NORM(u::AbstractArray{<:Tracker.TrackedReal,N},t) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline function ODE_DEFAULT_NORM(u::Array{<:Tracker.TrackedReal,N},t) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
    end
    @inline ODE_DEFAULT_NORM(u::Tracker.TrackedReal,t) = abs(value(u))

    # Support TrackedReal time, don't drop tracking on the adaptivity there
    @inline function ODE_DEFAULT_NORM(u::Tracker.TrackedArray,t::Tracker.TrackedReal) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip(u,Iterators.repeated(t))) / length(u))
    end
    @inline function ODE_DEFAULT_NORM(u::AbstractArray{<:Tracker.TrackedReal,N},t::Tracker.TrackedReal) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip(u,Iterators.repeated(t))) / length(u))
    end
    @inline function ODE_DEFAULT_NORM(u::Array{<:Tracker.TrackedReal,N},t::Tracker.TrackedReal) where {N}
      sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip(u,Iterators.repeated(t))) / length(u))
    end
    @inline ODE_DEFAULT_NORM(u::Tracker.TrackedReal,t::Tracker.TrackedReal) = abs(u)
  end

  # Piracy, should get upstreamed
  @require CuArrays="3a865a2d-5b23-5a0f-bc46-62713ec82fae" begin
    cuify(x::AbstractArray) = CuArrays.CuArray(x)
    function LinearAlgebra.ldiv!(x::CuArrays.CuArray,_qr::CuArrays.CUSOLVER.CuQR,b::CuArrays.CuArray)
      _x = UpperTriangular(_qr.R) \ (_qr.Q' * reshape(b,length(b),1))
      x .= vec(_x)
      CuArrays.unsafe_free!(_x)
    end
    function findall_events(affect!,affect_neg!,prev_sign::CuArrays.CuArray,next_sign::CuArrays.CuArray)
      hasaffect::Bool = affect! !== nothing
      hasaffectneg::Bool = affect_neg! !== nothing
      f = (p,n)-> ((p < 0 && hasaffect) || (p > 0 && hasaffectneg)) && p*n<=0
      A = map(f,prev_sign,next_sign)
      out = findall(A)
      CuArrays.unsafe_free!(A)
      out
    end
  end
end


text:
module ArraysTests

using Test

@testset "Interfaces" begin include("InterfaceTests.jl") end

@testset "CachedArrays" begin include("CachedArraysTests.jl") end

@testset "Kernels" begin include("KernelsTests.jl") end

@testset "Apply" begin include("ApplyTests.jl") end

@testset "CompressedArrays" begin include("CompressedArraysTests.jl") end

@testset "LocalToGlobalArrays" begin include("LocalToGlobalArraysTests.jl") end

@testset "Tables" begin include("TablesTests.jl") end

end # module

ext: jl


text:
M = 10
function generate_f(N, d, k; X = rand(N, d))
    f = rand(MvNormal(zeros(N), kernelmatrix(k, X, obsdim = 1) + 1e-5I))
    return X, f
end

function tests(model1, model2, X, f, y, problem)
    train!(model1, 1)
    L = AGP.objective(model1)
    @test L < 0
    train!(model1, 5)
    @test testconv(model1, problem, X, f, y)
    @test all(proba_y(model1, X)[2] .> 0)
    train!(model2, 6)
    @test testconv(model2, problem, X, f, y)
    @test all(proba_y(model2, X)[2] .> 0)
end

function tests(model1::OnlineSVGP, model2, X, f, y, problem)
    for (X_, y_) in eachbatch((X, y), obsdim = 1, size = 10)
        train!(model1, X_, y_, iterations = 5)
    end
    L = AGP.objective(model1)
    @test testconv(model1, problem, X, f, y)
    @test all(proba_y(model1, X)[2] .> 0)
    for (X_, y_) in eachbatch((X, y), obsdim = 1, size = 10)
        train!(model2, X_, y_, iterations = 5)
    end
    @test testconv(model2, problem, X, f, y)
    @test all(proba_y(model2, X)[2] .> 0)
end

function tests(
    model1::AbstractGP{T,<:BayesianSVM},
    model2,
    X,
    f,
    y,
    problem,
) where {T}
    train!(model1, 1)
    L = AGP.objective(model1)
    train!(model1, 5)
    @test testconv(model1, problem, X, f, y)
    @test all(proba_y(model1, X)[2] .> 0)
    train!(model2, 6)
    @test testconv(model2, problem, X, f, y)
    @test all(proba_y(model2, X)[2] .> 0)
end

function tests_likelihood(
    l::AbstractLikelihood,
    ltype::Type{<:AbstractLikelihood},
    dict::Dict,
    floattypes,
    problem,
    nLatent,
    X,
    f,
    y,
    k
)
    @testset "VGP" begin
        for floattype in floattypes
            dictvgp = dict["VGP"]
            @testset "AnalyticVI" begin
                if dictvgp["AVI"]
                    model = VGP(
                        X,
                        y,
                        k,
                        l,
                        AnalyticVI(),
                        optimiser = false,
                        verbose = 0,
                    )
                    @test eltype(model) == floattype
                    @test AGP.likelihood(model) isa ltype
                    @test AGP.inference(model) isa AnalyticVI
                    @test AGP.getf(model) isa NTuple{nLatent,AGP.VarLatent}
                    @test AGP.output(model) isa AbstractVector
                    @test AGP.input(model) isa AbstractVector
                    @test AGP.nLatent(model) == nLatent
                    model_opt = VGP(
                        X,
                        y,
                        k,
                        l,
                        AnalyticVI(),
                        optimiser = true,
                        verbose = 0,
                    )
                    tests(model, model_opt, X, f, y, problem)
                else
                    @test_throws ErrorException VGP(
                        X,
                        y,
                        k,
                        l,
                        AnalyticVI(),
                        optimiser = false,
                        verbose = 0,
                    )
                end
            end  # Analytic VI
            @testset "NumericalVI" begin
                if dictvgp["QVI"]
                    model = VGP(X, y, k, l, QuadratureVI(), optimiser = false)
                    @test eltype(model) == floattype
                    @test AGP.likelihood(model) isa ltype
                    @test AGP.inference(model) isa QuadratureVI
                    @test AGP.getf(model) isa NTuple{nLatent,AGP.VarLatent}
                    @test AGP.output(model) isa AbstractVector
                    @test AGP.input(model) isa AbstractVector
                    @test AGP.nLatent(model) == nLatent
                    model_opt = VGP(
                        X,
                        y,
                        k,
                        l,
                        QuadratureVI(),
                        optimiser = true,
                        verbose = 0,
                    )
                    tests(model, model_opt, X, f, y, problem)
                else
                    @test_throws ErrorException VGP(
                        X,
                        y,
                        k,
                        l,
                        QuadratureVI(),
                        optimiser = false,
                        verbose = 0,
                    )
                end
                if dictvgp["MCVI"]
                    model =
                        VGP(X, y, k, l, MCIntegrationVI(), optimiser = false)
                        @test eltype(model) == floattype
                        @test AGP.likelihood(model) isa ltype
                        @test AGP.inference(model) isa MCIntegrationVI
                        @test AGP.getf(model) isa NTuple{nLatent,AGP.VarLatent}
                        @test AGP.output(model) isa AbstractVector
                        @test AGP.input(model) isa AbstractVector
                        @test AGP.nLatent(model) == nLatent
                    model_opt = VGP(
                        X,
                        y,
                        k,
                        l,
                        MCIntegrationVI(),
                        optimiser = true,
                        verbose = 0,
                    )
                    tests(model, model_opt, X, f, y, problem)
                else
                    @test_throws ErrorException VGP(
                        X,
                        y,
                        k,
                        l,
                        MCIntegrationVI(),
                        optimiser = false,
                        verbose = 0,
                    )
                end
            end # Loop on Numerical VI
        end # Loop on float types
    end # VGP
    @testset "OSVGP" begin
        dictosvgp = dict["OSVGP"]
        for floattype in floattypes
            @testset "AnalyticVI" begin
                if dictosvgp["AVI"]
                    global model = OnlineSVGP(
                        k,
                        l,
                        AnalyticVI(),
                        AGP.InducingPoints.OIPS(),
                        optimiser = false,
                        verbose = 0,
                    )
                    @test eltype(model) == floattype
                    @test AGP.likelihood(model) isa ltype
                    @test AGP.inference(model) isa AnalyticVI
                    @test AGP.getf(model) isa NTuple{nLatent, AGP.OnlineVarLatent}
                    @test AGP.nLatent(model) == nLatent
                    global model_opt = OnlineSVGP(
                        k,
                        l,
                        AnalyticVI(),
                        AGP.InducingPoints.OIPS(),
                        optimiser = true,
                        verbose = 0,
                    )
                    tests(model, model_opt, X, f, y, problem)
                else
                    @test_throws ErrorException OnlineSVGP(
                        k,
                        l,
                        AnalyticVI(),
                        UniGrid(10),
                        optimiser = false,
                        verbose = 0,
                    )
                end
            end  # Analytic VI
            @testset "NumericalVI" begin
                if dictosvgp["QVI"]
                    model = OnlineSVGP(
                        k,
                        l,
                        QuadratureVI(),
                        UniGrid(10),
                        optimiser = false,
                    )
                    @test eltype(model) == floattype
                    @test AGP.likelihood(model) isa ltype
                    @test AGP.inference(model) isa QuadratureVI
                    @test AGP.getf(model) isa NTuple{nLatent, AGP.OnlineVarLatent}
                    @test AGP.nLatent(model) == nLatent
                    model_opt = OnlineSVGP(
                        k,
                        l,
                        QuadratureVI(),
                        UniGrid(10),
                        optimiser = true,
                        verbose = 0,
                    )
                    tests(model, model_opt, X, f, y, problem)
                else
                    @test_throws ErrorException OnlineSVGP(
                        k,
                        l,
                        QuadratureVI(),
                        UniGrid(10),
                        optimiser = false,
                        verbose = 0,
                    )
                end
                if dictosvgp["MCVI"]
                    model = OnlineSVGP(
                        k,
                        l,
                        MCIntegrationVI(),
                        UniGrid(10),
                        optimiser = false,
                    )
                    @test eltype(model) == floattype
                    @test AGP.likelihood(model) isa ltype
                    @test AGP.inference(model) isa MCIntegrationVI
                    @test AGP.getf(model) isa NTuple{nLatent, AGP.SparseVarLatent}
                    @test AGP.nLatent(model) == nLatent
                    model_opt = OnlineSVGP(
                        k,
                        l,
                        MCIntegrationVI(),
                        UniGrid(10),
                        optimiser = true,
                        verbose = 0,
                    )
                    tests(model, model_opt, X, f, y, problem)
                else
                    @test_throws ErrorException OnlineSVGP(
                        k,
                        l,
                        MCIntegrationVI(),
                        UniGrid(10),
                        optimiser = false,
                        verbose = 0,
                    )
                end
            end # Loop on Numerical VI
        end # Loop on float types
    end #Loop on Online SVGP
    @testset "SVGP" begin
        dictsvgp = dict["SVGP"]
        for floattype in floattypes
            @testset "AnalyticVI" begin
                if dictsvgp["AVI"]
                    model = SVGP(
                        X,
                        y,
                        k,
                        l,
                        AnalyticVI(),
                        M,
                        optimiser = false,
                        verbose = 0,
                    )
                    @test eltype(model) == floattype
                    @test AGP.likelihood(model) isa ltype
                    @test AGP.inference(model) isa AnalyticVI
                    @test AGP.getf(model) isa NTuple{nLatent, AGP.SparseVarLatent}
                    @test AGP.output(model) isa AbstractVector
                    @test AGP.input(model) isa AbstractVector
                    @test AGP.nLatent(model) == nLatent
                    model_opt = SVGP(
                        X,
                        y,
                        k,
                        l,
                        AnalyticVI(),
                        M,
                        optimiser = true,
                        Zoptimiser = true,
                        verbose = 0,
                    )
                    tests(model, model_opt, X, f, y, problem)

                    model_svi = SVGP(
                        X,
                        y,
                        k,
                        l,
                        AnalyticSVI(10),
                        M,
                        optimiser = false,
                        verbose = 0,
                    )
                    tests(model_svi, model, X, f, y, problem)
                else
                    @test_throws ErrorException SVGP(
                        X,
                        y,
                        k,
                        l,
                        AnalyticVI(),
                        M,
                        optimiser = false,
                        verbose = 0,
                    )
                end
            end # Analytic VI
            @testset "NumericalVI" begin
                if dictsvgp["QVI"]
                    model =
                        SVGP(X, y, k, l, QuadratureVI(), M, optimiser = false)
                        @test eltype(model) == floattype
                        @test AGP.likelihood(model) isa ltype
                        @test AGP.inference(model) isa QuadratureVI
                        @test AGP.getf(model) isa NTuple{nLatent,AGP.SparseVarLatent}
                        @test AGP.output(model) isa AbstractVector
                        @test AGP.input(model) isa AbstractVector
                        @test AGP.nLatent(model) == nLatent
                    model_opt = SVGP(
                        X,
                        y,
                        k,
                        l,
                        QuadratureVI(),
                        M,
                        optimiser = true,
                        Zoptimiser = true,
                        verbose = 0,
                    )
                    tests(model, model_opt, X, f, y, problem)

                    model_svi = SVGP(
                        X,
                        y,
                        k,
                        l,
                        QuadratureSVI(10),
                        M,
                        optimiser = false,
                        verbose = 0,
                    )
                    tests(model_svi, model, X, f, y, problem)
                else
                    @test_throws ErrorException SVGP(
                        X,
                        y,
                        k,
                        l,
                        QuadratureVI(),
                        M,
                        optimiser = false,
                        verbose = 0,
                    )
                end
                if dictsvgp["MCVI"]
                    model = SVGP(
                        X,
                        y,
                        k,
                        l,
                        MCIntegrationVI(),
                        M,
                        optimiser = false,
                    )
                    @test eltype(model) == floattype
                    @test AGP.likelihood(model) isa ltype
                    @test AGP.inference(model) isa MCIntegrationVI
                    @test AGP.getf(model) isa NTuple{nLatent, AGP.SparseVarLatent}
                    @test AGP.output(model) isa AbstractVector
                    @test AGP.input(model) isa AbstractVector
                    @test AGP.nLatent(model) == nLatent
                    model_opt = SVGP(
                        X,
                        y,
                        k,
                        l,
                        MCIntegrationVI(),
                        M,
                        optimiser = true,
                        Zoptimiser = true,
                        verbose = 0,
                    )
                    tests(model, model_opt, X, f, y, problem)

                    model_svi = SVGP(
                        X,
                        y,
                        k,
                        l,
                        MCIntegrationSVI(10),
                        10,
                        optimiser = false,
                        verbose = 0,
                    )
                    tests(model_svi, model, X, f, y, problem)
                else
                    @test_throws ErrorException SVGP(
                        X,
                        y,
                        k,
                        l,
                        MCIntegrationVI(),
                        M,
                        optimiser = false,
                        verbose = 0,
                    )
                end
            end # Loop Numerical VI
        end # Loop on float types
    end # SVGP
    @testset "MCGP" begin
        dictmcgp = dict["MCGP"]
        for floattype in floattypes
            @testset "Gibbs Sampling" begin
                if dictmcgp["Gibbs"]
                    model = MCGP(X, y, k, l, GibbsSampling())
                    @test eltype(model) == floattype
                    @test AGP.likelihood(model) isa ltype
                    @test AGP.inference(model) isa GibbsSampling
                    @test AGP.getf(model) isa NTuple{nLatent,AGP.SampledLatent}
                    @test AGP.output(model) isa AbstractVector
                    @test AGP.input(model) isa AbstractVector
                    @test AGP.nLatent(model) == nLatent
                    samples = sample(model, 100; progress=false)
                    @test_broken samples2 =
                        sample(model, 100; cat=true, progress=false)
                else
                    @test_throws ErrorException MCGP(
                        X,
                        y,
                        k,
                        l,
                        GibbsSampling(),
                    )
                end
                if dictmcgp["HMC"]
                    model = MCGP(X, y, k, l, HMCSampling())
                    @test eltype(model) == floattype
                    @test AGP.likelihood(model) isa ltype
                    @test AGP.inference(model) isa HMCSampling
                    @test AGP.getf(model) isa NTuple{nLatent,AGP.SampledLatent}
                    @test AGP.output(model) isa AbstractVector
                    @test AGP.input(model) isa AbstractVector
                    @test AGP.nLatent(model) == nLatent
                    samples = sample(model, 20; progress=false)
                    @test_broken samples2 =
                        sample(model, 20; cat=true)
                else
                    @test_throws ErrorException MCGP(X, y, k, l, HMCSampling())
                end
            end # Gibbs Sampling
        end # Loop floattypes
    end # MCGP
end

function testconv(
    model::AbstractGP,
    problem_type::String,
    X::AbstractArray,
    f::AbstractArray,
    y::AbstractArray,
)
    μ, Σ = predict_f(model, X, cov = true, diag = false)
    μ, diagΣ = predict_f(model, X, cov = true, diag = true)
    y_pred = predict_y(model, X)
    py_pred = proba_y(model, X)
    if problem_type == "Regression"
        err = mean(abs.(y_pred - f))
        # @info "Regression Error" err
        return err < 1.5
    elseif problem_type == "Classification"
        err = mean(y_pred .!= y)
        # @info "Classification Error" err
        return err < 0.5
    elseif problem_type == "MultiClass"
        err = mean(y_pred .!= y)
        # @info "Multiclass Error" err
        return err < 0.9
    elseif problem_type == "Poisson" || problem_type == "NegBinomial"
        err = mean(abs.(y_pred - y))
        # @info "Event Error" err
        return err < 20.0
    else
        throw(ErrorException("Wrong problem description $problem_type"))
    end
end


text:
using Documenter, GenomeGraphs, Pkg

makedocs(
    modules = [GenomeGraphs, GenomeGraphs.Graphs, GenomeGraphs.GraphIndexes],
    format = Documenter.HTML(),
    sitename = "GenomeGraphs.jl",
    authors = replace(join(Pkg.TOML.parsefile("Project.toml")["authors"], ", "), r" <.*?>" => "" ),
    pages = [
        "Home" => "index.md",
        "Manual" => [
            "Installation" => "man/install.md"
            "Genome Assembly" => "man/assembly.md"
        ],
        "API" => [
            "Graphs submodule" => "api/Graphs.md",
            "GraphIndexes submodule" => "api/GraphIndexes.md"
        ]
    ],

)

deploydocs(
    repo = "github.com/BioJulia/GenomeGraphs.jl.git",
    push_preview = true,
    deps = nothing,
    make = nothing
)

max_issues_repo_name: BioJulia/BioSequenceGraphs

max_forks_repo_name: BioJulia/BioSequenceGraphs


text:
## S H E E T
# Empty sheet
Sheet(name::String, wb::Workbook) = Sheet(name, wb,
    OrderedDict{LocalPosition,Cell}() # map
)

# Getters
get(S::Sheet, p::LocalPosition, default = missing) = Base.get(S.map, p, default)
get(S::Sheet, R::LocalRangeType) = get.([S],collect(R))
Base.getindex(S::Sheet, p::LocalRef) = get(S, p)
Base.getindex(S::Sheet, p::String) = get(S, LocalRef(p))
Base.getindex(S::Sheet, i, j) = get(S, LocalPosition(i,j))

# No setters for cells: only set at the workbook level
Base.setindex!(S::Sheet, f, p::Union{LocalPosition,AbstractString}) = setindex!(S.wb, f, GlobalPosition(S.name,p))
Base.setindex!(S::Sheet, f, i,j)   = setindex!(S.wb, f, GlobalPosition(S.name,LocalPosition(i,j)))
Base.setindex!(S::Sheet, f, xs...) = setindex!(S,string(f),xs...)

## O U T P U T
const minPos = LocalPosition(1,1)
const maxPos = LocalPosition(typemax(UInt),typemax(UInt))
function hull(pp)
    isempty(pp) && return LocalPosition(0,0)
    LocalPosition(maximum(getfield.(pp,:row)),maximum(getfield.(pp,:column)))
end

function array(S::Sheet, fn = value)
    isempty(S.map.keys) && return nothing
    fn.(S[LocalPosition(1,1):hull(S.map.keys)])
end

function Base.string(S::Sheet, field = value)
    R = array(S,field)
    isnothing(R) && return "Empty sheet"
    pretty_table(String, R, colNumAsStr.(1:size(R,2)), row_names=1:size(R,1))
end
Base.show(io::IO, S::Sheet) = Base.print(io, string(S))

function Base.size(S::Sheet)
    h = hull(S.map.keys)
    (h.row,h.column)
end


text:
mutable struct Location
    lonlat::Tuple{Float64,Float64}
    placeid::String
    address::String
    locationtype::String
    tags::Vector{String}
end

function Base.show(io::IO, loc::Location)
    pprint(s::String, width::Int=50) =
        length(s) > width ? string(s[1:(width-3)], "...") : s
    println(io, "Location")
    println(io, "========")
    println(io, " *    type: ", pprint(loc.locationtype))
    println(io, " * address: ", pprint(loc.address))
    println(io, " *    tags: (",length(loc.tags),") ", pprint(join(loc.tags," | ")))
end

function getlocation(loc::Dict)
    geom = loc["geometry"]
    Location(
        (geom["location"]["lng"], geom["location"]["lat"]), loc["place_id"],
        loc["formatted_address"], geom["location_type"], loc["types"]
    )
end

function geocode(;kwargs...)
    resultjson = requestjson("geocode"; kwargs...)
    if resultjson["status"] == "OK"
        return getlocation.(resultjson["results"])
    else
        return Location[]
    end
end

"""

Geocoding is the process of converting addresses (like "1600 Amphitheatre
Parkway, Mountain View, CA") into geographic coordinates (like latitude
37.423021 and longitude -122.083739).

https://developers.google.com/maps/documentation/geocoding/intro

## Parameters
* `key`: your google maps API key
* `address`: The street address that you want to geocode, in the format used by
    the national postal service of the country concerned. Additional address
    elements such as business names and unit, suite or floor numbers should be
    avoided. Please refer to the FAQ for additional guidance.
    https://developers.google.com/maps/faq#geocoder_queryformat

## Remarks
The geocoder is designed to map street addresses to geographical coordinates.
We therefore recommend that you format geocoder requests in accordance with the
following guidelines to maximize the likelihood of a successful query:

* Specify addresses in accordance with the format used by the national postal
    service of the country concerned.
* Do not specify additional address elements such as business names, unit
    numbers, floor numbers, or suite numbers that are not included in the
    address as defined by the postal service of the country concerned. Doing so
    may result in responses with ZERO_RESULTS.
* Use the street number of a premise over the building name where possible.
* Use street number addressing over specifying cross streets where possible.
* Do not provide 'hints' such as nearby landmarks.

https://developers.google.com/maps/faq#geocoder_queryformat

"""
geocode(address::String; kwargs...) = geocode(address=address; kwargs...)

"""

The term geocoding generally refers to translating a human-readable address into
a location on a map. The process of doing the opposite, translating a location
on the map into a human-readable address, is known as reverse geocoding.

## Parameters
* `lonlat`: the longitude and latitude values specifying the location for which
    you wish to obtain the closest, human-readable address.

## Keyword Arguments
* (optional but recommended) `key`: your google maps API key

"""
geocode(lon::Float64, lat::Float64; kwargs...) =
    geocode(latlng = string(lat,",",lon); kwargs...)

geocode(lonlat::Tuple{Float64,Float64}; kwargs...) =
    geocode(lonlat...; kwargs...)


text:
# Interfaces
# ==========
#
# Interface functions for pairwise sequence alignment algorithms.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioAlignments.jl/blob/master/LICENSE.md

# include algorithms
include("alignment.jl")
include("result.jl")
include("algorithms/common.jl")
include("algorithms/needleman_wunsch.jl")
include("algorithms/banded_needleman_wunsch.jl")
include("algorithms/smith_waterman.jl")
include("algorithms/edit_distance.jl")
include("algorithms/hamming_distance.jl")

"""
    pairalign(type, seq, ref, model, [options...])

Run pairwise alignment between two sequences: `seq` and `ref`.

Available `type`s are:

* `GlobalAlignment()`
* `LocalAlignment()`
* `SemiGlobalAlignment()`
* `OverlapAlignment()`
* `EditDistance()`
* `LevenshteinDistance()`
* `HammingDistance()`

`GlobalAlignment`, `LocalAlignment`, `SemiGlobalAlignment`, and
`OverlapAlignment` are problem that maximizes alignment score between two
sequences.  Therefore, `model` should be an instance of `AbstractScoreModel`
(e.g. `AffineGapScoreModel`).

`EditDistance`, `LevenshteinDistance`, and `HammingDistance` are problem that
minimizes alignment cost between two sequences.  As for `EditDistance`, `model`
should be an instance of `AbstractCostModel` (e.g. `CostModel`).
`LevenshteinDistance` and `HammingDistance` have predefined a cost model,
so users cannot specify a cost model for these alignment types.

When you pass the `score_only=true` or `distance_only=true` option to
`pairalign`, the result of pairwise alignment holds alignment score/distance
only.  This may enable some algorithms to run faster than calculating full
alignment result.  Other available `options` are documented for each alignemnt
type.

Example
-------

    using BioSequences
    using BioAlignments

    # create affine gap scoring model
    affinegap = AffineGapScoreModel(
        match=5,
        mismatch=-4,
        gap_open=-5,
        gap_extend=-3
    )

    # run global alignment between two DNA sequences
    pairalign(GlobalAlignment(), dna"AGGTAG", dna"ATTG", affinegap)

    # run local alignment between two DNA sequences
    pairalign(LocalAlignment(), dna"AGGTAG", dna"ATTG", affinegap)

    # you cannot specify a cost model in LevenshteinDistance
    pairalign(LevenshteinDistance(), dna"AGGTAG", dna"ATTG")

See also: `AffineGapScoreModel`, `CostModel`
"""
function pairalign end

function pairalign(::GlobalAlignment, a::S1, b::S2, score::AffineGapScoreModel{T};
                   score_only::Bool=false,
                   banded::Bool=false,
                   lower_offset::Int=0,
                   upper_offset::Int=0) where {S1,S2,T}
    m = length(a)
    n = length(b)
    if banded
        if m > n
            L = m - n + lower_offset
            U = upper_offset
        else
            L = lower_offset
            U = n - m + upper_offset
        end
        bnw = BandedNeedlemanWunsch{T}(m, n, L, U)
        score = run!(bnw, a, b, score.submat, score.gap_open, score.gap_extend)
        if score_only
            return PairwiseAlignmentResult{S1,S2}(score, true)
        else
            a′ = traceback(bnw, a, b, (m, n))
            return PairwiseAlignmentResult(score, true, a′, b)
        end
    else
        nw = NeedlemanWunsch{T}(m, n)
        score = run!(nw, a, b, score.submat, score.gap_open, score.gap_extend)
        if score_only
            return PairwiseAlignmentResult{S1,S2}(score, true)
        else
            a′ = traceback(nw, a, b, (m, n))
            return PairwiseAlignmentResult(score, true, a′, b)
        end
    end
end

function pairalign(::SemiGlobalAlignment, a::S1, b::S2, score::AffineGapScoreModel{T};
                   score_only::Bool=false) where {S1,S2,T}
    m = length(a)
    n = length(b)
    nw = NeedlemanWunsch{T}(m, n)
    gap_open = score.gap_open
    gap_extend = score.gap_extend
    score = run!(nw, a, b, score.submat,
        T(0), T(0), gap_open, gap_extend, T(0), T(0),
        gap_open, gap_extend, gap_open, gap_extend, gap_open, gap_extend,
    )
    if score_only
        return PairwiseAlignmentResult{S1,S2}(score, true)
    else
        a′ = traceback(nw, a, b, (m, n))
        return PairwiseAlignmentResult(score, true, a′, b)
    end
end

function pairalign(::OverlapAlignment, a::S1, b::S2, score::AffineGapScoreModel{T};
                   score_only::Bool=false) where {S1,S2,T}
    m = length(a)
    n = length(b)
    nw = NeedlemanWunsch{T}(m, n)
    score = run!(nw, a, b, score.submat,
        T(0), T(0), score.gap_open, score.gap_extend, T(0), T(0),
        T(0), T(0), score.gap_open, score.gap_extend, T(0), T(0),
    )
    if score_only
        return PairwiseAlignmentResult{S1,S2}(score, true)
    else
        a′ = traceback(nw, a, b, (m, n))
        return PairwiseAlignmentResult(score, true, a′, b)
    end
end

function pairalign(::LocalAlignment, a::S1, b::S2, score::AffineGapScoreModel{T};
                   score_only::Bool=false) where {S1,S2,T}
    sw = SmithWaterman{T}(length(a), length(b))
    score, endpos = run!(sw, a, b, score.submat, score.gap_open, score.gap_extend)
    if score_only
        return PairwiseAlignmentResult{S1,S2}(score, true)
    else
        a′ = traceback(sw, a, b, endpos)
        return PairwiseAlignmentResult(score, true, a′, b)
    end
end

function pairalign(::EditDistance, a::S1, b::S2, cost::CostModel;
                   distance_only::Bool=false) where {S1,S2}
    dist, trace, endpos = edit_distance(a, b, cost.submat, cost.insertion, cost.deletion)
    if distance_only
        return PairwiseAlignmentResult{S1,S2}(dist, false)
    else
        a′ = edit_traceback(a, b, trace, endpos)
        return PairwiseAlignmentResult(dist, false, a′, b)
    end
end

function pairalign(::LevenshteinDistance, a::S1, b::S2;
                   distance_only::Bool=false) where {S1,S2}
    unitcost = CostModel(
        DichotomousSubstitutionMatrix{Int}(0, 1),
        insertion=1,
        deletion=1
    )
    return pairalign(EditDistance(), a, b, unitcost, distance_only=distance_only)
end

function pairalign(::HammingDistance, a::S1, b::S2;
                   distance_only::Bool=false) where {S1,S2}
    dist, anchors = hamming_distance(Int, a, b)
    if distance_only
        return PairwiseAlignmentResult{S1,S2}(dist, false)
    else
        a′ = AlignedSequence(a, anchors)
        return PairwiseAlignmentResult(dist, true, a′, b)
    end
end


text:
#=
function assign_lines_to_orders(line_list::DataFrame, order_info::DataFrame; v_center::Real=0, Δv_to_avoid_tellurics::Real = RvSpectMLBase.max_bc )
  @assert hasproperty(line_list, :lambda)
  @assert hasproperty(line_list, :lambda_lo)
  @assert hasproperty(line_list, :lambda_hi)
  ll1 = copy(line_list)
  ll2 = copy(line_list)
  ll3 = copy(line_list)
  ll1[!,:order] .= 0
  ll2[!,:order] .= 0
  ll3[!,:order] .= 0
  #if !(hasproperty(line_list,:lambda_lo) && hasproperty(line_list,:lambda_hi))
  add_line_boundaries_to_line_list(line_list; Δv_to_avoid_tellurics=Δv_to_avoid_tellurics, v_center_to_avoid_tellurics=v_center )
  #end
  for (i,line) in enumerate(eachrow(line_list))
    order_indices = findall( map(order-> order.λ_min <= line.lambda_lo <= order.λ_max && order.λ_min <= line.lambda_hi <= order.λ_max, eachrow(order_info) ) )
    if length(order_indices) == 1
      ll1[i,:order] = order_info.order[order_indices[1]]
    elseif length(order_indices) == 2
      ll1[i,:order] = order_info.order[order_indices[1]]
      ll2[i,:order] = order_info.order[order_indices[2]]
    elseif length(order_indices) == 3
      ll1[i,:order] = order_info.order[order_indices[1]]
      ll2[i,:order] = order_info.order[order_indices[2]]
      ll3[i,:order] = order_info.order[order_indices[3]]
    elseif length(order_indices) > 3
      @warn "One line appears in >3 orders.  Not implemented yet."
      println("line = ", line, "  order_indices = ", order_indices, "  line.order = ",order_info.order[order_indices])
      #=  df_tmp = DataFrame()
      for key in keys(eachcol(line_list))
        df_tmp[!,key] = fill(line[key], length(order_indices)-1 )
      end
      df_tmp[:order] .= order_indices[2:end]
      append!(line_list,df_tmp)
      =#
      #break
    end
  end
  append!(ll1,ll2)
  append!(ll1,ll3)
  ll1 |> @filter( _.order >= 1 ) |> @orderby(_.order) |> @thenby(_.lambda_lo) |> DataFrame
  #ll1 |> @filter( _.order >= 1 ) |> @orderby(_.lambda_lo) |> @thenby(_.order) |> DataFrame
end
=#

#=
function add_line_boundaries_to_line_list(line_list::DataFrame; Δv_to_avoid_tellurics::Real = 0.0, v_center_to_avoid_tellurics::Real = 0.0 )
    @assert 0.5*RvSpectMLBase.max_bc_earth_rotation <= Δv_to_avoid_tellurics <= 4*RvSpectMLBase.max_bc
    #line_list_to_search_for_tellurics = copy(line_list)
    line_list_to_search_for_tellurics = line_list
    @assert hasproperty(line_list_to_search_for_tellurics,:lambda)
    if hasproperty(line_list_to_search_for_tellurics,:lambda_lo) && hasproperty(line_list_to_search_for_tellurics,:lambda_hi)
        line_list_to_search_for_tellurics.lambda_lo = line_list_to_search_for_tellurics.lambda_lo.*(calc_doppler_factor(v_center_to_avoid_tellurics)/calc_doppler_factor(Δv_to_avoid_tellurics))
        line_list_to_search_for_tellurics.lambda_hi = line_list_to_search_for_tellurics.lambda_hi.*(calc_doppler_factor(v_center_to_avoid_tellurics)*calc_doppler_factor(Δv_to_avoid_tellurics))
    else
        line_list_to_search_for_tellurics.lambda    = line_list_to_search_for_tellurics.lambda.*(calc_doppler_factor(v_center_to_avoid_tellurics))
        line_list_to_search_for_tellurics.lambda_lo = line_list_to_search_for_tellurics.lambda.*(calc_doppler_factor(v_center_to_avoid_tellurics)/calc_doppler_factor(Δv_to_avoid_tellurics))
        line_list_to_search_for_tellurics.lambda_hi = line_list_to_search_for_tellurics.lambda.*(calc_doppler_factor(v_center_to_avoid_tellurics)*calc_doppler_factor(Δv_to_avoid_tellurics))
    end
    return line_list_to_search_for_tellurics
end
=#

function calc_snr_weights_for_lines!(line_list::DataFrame, spectra::AAS ; v_center::Real=0 ) where { AS<:AbstractSpectra2D, AAS<:AbstractVector{AS} }
  num_obs = length(spectra)
  @assert num_obs >= 1
  num_cols = size(first(spectra).λ,1)
  num_lines = size(line_list,1)
  snr_list = zeros(num_obs, num_lines)
  depth = zeros(num_obs, num_lines)
  exp_sigma_rv = zeros(num_obs, num_lines)
  line_list[!,:weight_snr] .= 0.0
  for (i,line) in enumerate(eachrow(line_list))
    order = line.order
    for obs_idx in 1:num_obs
      boost_factor = calc_doppler_factor(v_center)
      # TODO: Do something with boost_factor for each observation
      idx_min = searchsortedfirst(view(spectra[obs_idx].λ,1:num_cols,order),line.lambda_lo)
      idx_max = searchsortedlast(view(spectra[obs_idx].λ,1:num_cols,order),line.lambda_hi)
      @assert 1 <= idx_min <= num_cols
      @assert 1 <= idx_max <= num_cols
      #println("# Line ", i, " at ", line.lambda, " obs ", obs_idx, " idx = ", idx_min, " - ", idx_max)
      snr_list[obs_idx,i] = RvSpectMLBase.calc_snr( spectra[obs_idx], idx_min:idx_max, order )
      #(depth[obs_idx,i], exp_sigma_rv[obs_idx,i]) = calc_depth_and_expected_rv_precission( spectra[obs_idx], idx_min:idx_max, order )
    end
  end
  snr_weights = vec(median(snr_list,dims=1)).^2
  sum_snr_weights = sum(snr_weights)
  sum_weight_inputs = sum(line_list[!,:weight])
  snr_weights .*= sum_weight_inputs/sum_snr_weights
  line_list[!,:weight_snr] .= snr_weights
  #=
  line_list[!,:exp_sigma_rv] .= vec(median(exp_sigma_rv,dims=1))
  exp_σ_rv_weights = 1.0 ./ line_list[!,:exp_sigma_rv].^2
  sum_σ_rv_weights = sum(exp_σ_rv_weights)
  exp_σ_rv_weights ./= sum_σ_rv_weights

  line_list[!,:weight_theory] .= exp_σ_rv_weights #.* snr_weights
  #line_list[!,:weight_theory] ./= sum(line_list[!,:weight_theory])
  =#

  line_list[!,:weight_input] = deepcopy(line_list[!,:weight])
  #line_list[!,:weight] = line_list[!,:weight_theory]
  line_list[!,:weight] .*= snr_weights
  #line_list[!,:depth] .= vec(median(depth,dims=1))
  return line_list
end

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs


text:
@doc raw"""
    CostProblem{T, Manifold, TCost} <: Problem{T}

speficy a problem for solvers just based on cost functions, i.e.
gradient free ones.
# Fields

* `M`            – a manifold ``\mathcal M``
* `cost` – a function ``F: \mathcal M → ℝ`` to minimize

# Constructors

    CostProblem(M, cost; evaluation=AllocatingEvaluation())

Generate a problem. While this Problem does not have any allocating functions,
the type `T` can be set for consistency reasons with other problems.

# See also
[`NelderMead`](@ref)
"""
struct CostProblem{T,mT<:AbstractManifold,Tcost} <: Problem{T}
    M::mT
    cost::Tcost
end
function CostProblem(
    M::mT, cost::T; evaluation::AbstractEvaluationType=AllocatingEvaluation()
) where {mT<:AbstractManifold,T}
    return CostProblem{typeof(evaluation),mT,T}(M, cost)
end
@doc raw"""
    NelderMeadOptions <: Options

Describes all parameters and the state of a Nealer-Mead heuristic based
optimization algorithm.

# Fields

The naming of these parameters follows the [Wikipedia article](https://en.wikipedia.org/wiki/Nelder–Mead_method)
of the Euclidean case. The default is given in brackets, the required value range
after the description

* `population` – an `Array{`point`,1}` of ``n+1`` points ``x_i``, ``i=1,…,n+1``, where ``n`` is the
  dimension of the manifold.
* `stopping_criterion` – ([`StopAfterIteration`](@ref)`(2000)`) a [`StoppingCriterion`](@ref)
* `α` – (`1.`) reflection parameter (``α > 0``)
* `γ` – (`2.`) expansion parameter (``γ > 0``)
* `ρ` – (`1/2`) contraction parameter, ``0 < ρ ≤ \frac{1}{2}``,
* `σ` – (`1/2`) shrink coefficient, ``0 < σ ≤ 1``
* `x` – (`p[1]`) - a field to collect the current best value
* `retraction_method` – `ExponentialRetraction()` the rectraction to use, defaults to
  the exponential map
* `inverse_retraction_method` - `LogarithmicInverseRetraction` an `inverse_retraction(M,x,y)` to use.

# Constructors

    NelderMead(M[, population]; kwargs)

construct a Nelder-Mead Option with a default popultion (if not provided) of set of `dimension(M)+1` random points.

In the constructor all fields (besides the population) are keyword arguments.
"""
mutable struct NelderMeadOptions{
    T,
    S<:StoppingCriterion,
    Tα<:Real,
    Tγ<:Real,
    Tρ<:Real,
    Tσ<:Real,
    TR<:AbstractRetractionMethod,
    TI<:AbstractInverseRetractionMethod,
} <: Options
    population::Vector{T}
    stop::S
    α::Tα
    γ::Tγ
    ρ::Tρ
    σ::Tσ
    x::T
    costs::Vector{Float64}
    retraction_method::TR
    inverse_retraction_method::TI
    function NelderMeadOptions(
        M::AbstractManifold,
        population::Vector{T};
        stopping_criterion::StoppingCriterion=StopAfterIteration(2000),
        α=1.0,
        γ=2.0,
        ρ=1 / 2,
        σ=1 / 2,
        retraction_method::AbstractRetractionMethod=default_retraction_method(M),
        inverse_retraction_method::AbstractInverseRetractionMethod=default_inverse_retraction_method(
            M
        ),
    ) where {T}
        return new{
            T,
            typeof(stopping_criterion),
            typeof(α),
            typeof(γ),
            typeof(ρ),
            typeof(σ),
            typeof(retraction_method),
            typeof(inverse_retraction_method),
        }(
            population,
            stopping_criterion,
            α,
            γ,
            ρ,
            σ,
            population[1],
            [],
            retraction_method,
            inverse_retraction_method,
        )
    end
end
function NelderMeadOptions(M::AbstractManifold; kwargs...)
    p = [random_point(M) for i in 1:(manifold_dimension(M) + 1)]
    return NelderMeadOptions(M, p; kwargs...)
end
@deprecate NelderMeadOptions(
    population,
    retraction_method::AbstractRetractionMethod=ExponentialRetraction(),
    inverse_retraction_method::AbstractInverseRetractionMethod=LogarithmicInverseRetraction(),
    kwargs...,
) NelderMeadOptions(
    DefaultManifold(2),
    population,
    retraction_method=retraction_method,
    inverse_retraction_method=inverse_retraction_method,
    kwargs...,
)
get_solver_result(o::NelderMeadOptions) = o.x


text:
"""
This is the corresponding tests file for the `srcfile.jl` found in the `src`
directory. Your test files should be set up approximately like this one, i.e.:

- declare your test dependencies (Test, your package, then other deps)

- Group related tests into `@testset`s

- `@test` unit functionality

- See if a `@test_throws` some exception

Again, please note that:

- `srcfile.jl` is just a placeholder name for your actual, eventual source file
  names

- `make test` from the repo top-level will run these for you.
"""

using Test
using PkgName

@testset "bogus" begin
    @test 1 == 1
end

@testset "sharpe" begin
    @test sharpe([0.03, 0.01]) == 0.0
    @test sharpe([0.03, 0.01], rf = 0.00) == 1.4142135623730951

    # Can't pass anything other than an Array{Float64} as x
    @test_throws MethodError sharpe(1)
end

max_stars_repo_name: ryapric/boilerplate

max_issues_repo_name: ryapric/boilerplate

max_forks_repo_name: ryapric/boilerplate


text:
@testset "General Manifold functions and fallbacks" begin
M = Euclidean(3)
x = RnPoint([1.,0.,0.])
y = RnPoint([0.,1.,0.])
ξ = RnTVector([1.,0.,1.])
η = RnTVector([0.,1.,1.])
@test copy(x) == x
@test copy(ξ) == ξ
@test getValue(2*ξ) == 2*getValue(ξ)
@test getValue(ξ*2) == getValue(ξ)*2
@test getValue(2/RnTVector([1.,1.])) == 2 ./ [1.,1.]
@test getValue(ξ/2) == getValue(ξ)/2
@test getValue(ξ+η) == getValue(ξ) + getValue(η)
@test getValue(ξ-η) == getValue(ξ) - getValue(η)
@test getValue(+η) == getValue(η)
@test getValue(-η) == -getValue(η)
@test x==x
@test ξ==ξ
@test geodesic(M,x,y,3) == RnPoint.([[1.,0.,0.],[.5,.5,0.],[0.,1.,0.]])
@test geodesic(M,x,y,[0.,.5,1.]) == RnPoint.([[1.,0.,0.],[.5,.5,0.],[0.,1.,0.]])
@test reflection(M,x,y) == RnPoint(2*getValue(x)-getValue(y))
#
# Test Fallbacks – check that errors are thrown if the manifold does not fit
struct TestManifold <: Manifold end
struct TestMPoint <: MPoint end
struct TestTVector <: TVector end
M = TestManifold()
x = TestMPoint()
y = TestMPoint()
ξ = TestTVector()
η = TestTVector()
@test_throws DomainError distance(M, x, y)
@test_throws DomainError dot(M, x, ξ, η)
@test_throws ErrorException dot(Euclidean(3), ξ, η)
@test_throws DomainError exp(M, x, ξ)
@test_throws DomainError log(M, x, y)
@test_throws DomainError parallelTransport(M, x, y, ξ)
@test_throws DomainError tangentONB(M, x, ξ)
@test_throws DomainError tangentONB(M, x, y)
@test_logs (:warn,"""No valitadion for a TestMPoint on TestManifold available. Continuing without
  validation. To turn this warning off, either deactivate the validate flag
  in (one of) your extended MPoints or implement a corresponding validation.""") validateMPoint(M,x)
@test_logs (:warn,"""No valitadion for a TestMPoint and a TestTVector on TestManifold available.
  Continuing without validation. To turn this warning off, either deactivate
  the validate flag in (one of) your extended TVectors or MPoints or
  implement a corresponding validation""") validateTVector(M,x,ξ)
@test_throws DomainError zeroTVector(M,x)
#
#
@test_throws DomainError manifoldDimension(M)
@test_throws DomainError manifoldDimension(x)
@test_throws DomainError getValue(x)
@test_throws DomainError getValue(ξ)
#
@test_throws ErrorException typeofMPoint(ξ)
@test_throws ErrorException typeofTVector(x)
@test_throws ErrorException typicalDistance(M)
end


text:
using FV3Kernels

@show s = SillyType(π / 4)

# stuff with s...


text:
using DataFrames
using CSV
using Plots
#%%
bp = open("backward.dat")
fp = open("forward.dat")
cp = open("central.dat")
bp_df=CSV.read(bp, header=false, delim=" ")

plot(bp_df[1],bp_df[2])
xaxis!(:log10)
yaxis!(:log10)

max_stars_repo_name: aurelio-amerio/corso-cpp


text:
using BenchmarkTools
using JultraDark: Grids, max_phase_grad

resol = 64

grids = Grids(1.0, resol)

@benchmark max_phase_grad(grids.ψx)


text:
@safetestset "Markdown rendering" begin

  @safetestset "String markdown rendering" begin
    using Genie
    using Genie.Renderer.Html
    using Markdown

    view = raw"""
# Hello
## Welcome to Genie""" |> Markdown.parse

    @test (Html.html(view, forceparse = true).body |> String) == "<!DOCTYPE html>
Hello

Welcome to Genie
"

    view = raw"""
# Hello
## Welcome to Genie, $name""" |> Markdown.parse

    @test (Html.html(view, name = "John").body |> String) == "<!DOCTYPE html>
Hello

Welcome to Genie, John
"

    layout = raw"""

Layout header

    <% @yield %>

Layout footer

"""

    @test (Html.html(view, layout = layout, name = "John").body |> String) == "<!DOCTYPE html>

Layout header

Hello

Welcome to Genie, John

Layout footer

"
  end;

  @safetestset "Template markdown rendering" begin
    using Genie, Genie.Renderer
    using Genie.Renderer.Html

    @test Html.html(filepath("views/view.jl.md"), numbers = [1, 1, 2, 3, 5, 8, 13]).body |> String == "<!DOCTYPE html>
There are 7

-> 1      -> 1      -> 2      -> 3      -> 5      -> 8      -> 13
"
    @test Html.html(filepath("views/view.jl.md"), layout = filepath("views/layout.jl.html"), numbers = [1, 1, 2, 3, 5, 8, 13]).body |> String ==
      "<!DOCTYPE html>

Layout header

There are 7

-> 1      -> 1      -> 2      -> 3      -> 5      -> 8      -> 13

Layout footer

"
  end

  @safetestset "Markdown rendering with embedded variables" begin
    using Genie, Genie.Renderer
    using Genie.Renderer.Html

    @test Html.html(filepath("views/view-vars.jl.md")).body |> String == "<!DOCTYPE html>
There are 7

-> 1      -> 1      -> 2      -> 3      -> 5      -> 8      -> 13
"
  end;

end;


text:
"""
    AnovaResult <: AnovaEffect

A set of values for an Anova factor which has been tested

`name` - the name of this factor

`ss` - sum of squares

`df` - degrees of freedom

`ms` - mean square

`f` - the F statistic

`p` - the probability of a Type I error (imcorrect rejection of null hypothesis)
"""
struct AnovaResult <: AnovaEffect
    name::String
    ss::Float64
    df::Float64
    ms::Float64
    f::Float64
    p::Float64
end

AnovaResult(factor, f, p) = AnovaResult(factor.name, factor.ss, factor.df, factor.ms, f, p)

import Base.isapprox
isapprox(x::AnovaResult, y::AnovaResult) =
    x.name == y.name &&
    x.ss ≈ y.ss &&
    x.df == y.df &&
    x.ms ≈ y.ms &&
    x.f ≈ y.f &&
    x.p ≈ y.p

export AnovaResult


text:
module MatrixMarket2

using Sparso

# Separate "infoonly" feature out of mmread() so that mmread() returns only a
# single type: SparseMatrixCSC, instead of a union type
function mmread_infoonly(filename)
#      Reads the contents of the Matrix Market file 'filename'
#      into a matrix, which will be either sparse or dense,
#      depending on the Matrix Market format indicated by
#      'coordinate' (coordinate sparse storage), or
#      'array' (dense array storage).  The data will be duplicated
#      as appropriate if symmetry is indicated in the header. (Not yet
#      implemented).
#
#      If infoonly is true information on the size and structure is
#      returned.
        filename1 = ASCIIString(filename)
        mmfile = open(filename1,"r")
        tokens = split(chomp(readline(mmfile)))
        if length(tokens) != 5 error("Not enough words on header line") end
        if tokens[1] != "%%MatrixMarket" error("Not a valid MatrixMarket header.") end
        (head1, rep, field, symm) = map(lowercase, tokens[2:5])
        if head1 != "matrix"
            error("This seems to be a MatrixMarket $head1 file, not a MatrixMarket matrix file")
        end
        if field != "real" error("non-float fields not yet allowed") end

        ll   = readline(mmfile)         # Read through comments, ignoring them
        while length(ll) > 0 && ll[1] == '%' ll = readline(mmfile) end
        dd     = int(split(ll))         # Read dimensions
        rows   = dd[1]
        cols   = dd[2]
        entries = rep == "coordinate" ? dd[3] : rows * cols

        return rows, cols, entries, rep, field, symm
end

function mmread(filename, symmetric_only = false, force_symmetric = false)
#      Reads the contents of the Matrix Market file 'filename'
#      into a matrix, which will be either sparse or dense,
#      depending on the Matrix Market format indicated by
#      'coordinate' (coordinate sparse storage), or
#      'array' (dense array storage).  The data will be duplicated
#      as appropriate if symmetry is indicated in the header or
#      if force_symmetric is set.

    # Read into a COO array (stored statically insded mm_read)
    sizes::Vector{Cint} = [0, 0, 0, 0]
    ccall((:load_matrix_market_step1, Sparso.libcsr), Void,
        (Ptr{Uint8}, Ptr{Cint}, Bool, Bool), filename, pointer(sizes), force_symmetric, false)

    is_symmetric::Cint = sizes[1] # 0/1: true/false
    if symmetric_only && is_symmetric == 0
        throw("mmread: matrix market file is NOT symmetric. Sparse cannot handle it now, since we need a matrix to be symmetric, and thus we can use its CSC representation as CSR")
    end

    n::Cint = sizes[2]
    m::Cint = sizes[3]
    nnz::Cint = sizes[4]  #Note: if symmetric, nnz includes elements in both upper and lower triangle

    v = Array(Cdouble, nnz)
    i = Array(Cint, nnz)
    j = Array(Cint, n + 1)

    A = SparseMatrixCSC{Cdouble, Cint}(m, n, j, i, v)

    # Convert the COO array to CSR array.
    ccall((:load_matrix_market_step2, Sparso.libcsr), Void,
        (Ptr{Uint8}, Ptr{Cdouble}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Bool),
        filename, pointer(j), pointer(i), pointer(v), pointer(sizes), true)

    #distance = div(nnz, 100); # print about 100 elements to check manually
    #ccall((:CSR_PrintSomeValues, Sparso.libcsr), Void,
        #(Cint, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}, Cint, Bool),
        #m, n, pointer(j), pointer(i), pointer(v), convert(Cint, distance), true)

    A
end

function mmread_reorder(filename)
    mtx_file = open(filename, "r")
    modification_time::Float64 = mtime(mtx_file)
    close(mtx_file)

    reordered_filename = string(filename, ".reordered")
    if isfile(reordered_filename)
        # The file format is as follows:
        #   last modification time of the original mtx file
        #   m n #non_zeros Ti_num Tv_num colptr rowval nzval
        #   PermutationVector ReversePermutationVector
        reordered_file = open(reordered_filename, "r")
        modification_time1 = read(reordered_file, Float64)
        if (modification_time == modification_time1)
            # The reordered file has the up to date info
            m = read(reordered_file, Int)
            n = read(reordered_file, Int)
            nnz = read(reordered_file, Int64)

            Ti_num = read(reordered_file, Int)
            Tv_num = read(reordered_file, Int)
            dict = Dict{Int, Any}(0 => Int, 1 => Int32, 2 => Int64, 3 => Float64)
            Ti = dict[Ti_num]
            Tv = dict[Tv_num]

            colptr = read(reordered_file, Ti, n + 1)
            rowval = read(reordered_file, Ti, nnz)
            nzval  = read(reordered_file, Tv, nnz)
            A = SparseMatrixCSC(m, n, colptr, rowval, nzval)

            P      = read(reordered_file, Cint, n)
            Pprime = read(reordered_file, Cint, n)
            close(reordered_file)

            return A, P, Pprime
        else
            # The reordered file is out of date
            close(reordered_file)
            rm(reordered_filename)
        end
    end

    # Read the file to build a CSC array, and reorder it
    A          = mmread(filename)
    m::Int     = A.m
    n::Int     = A.n
    nnz::Int64 = size(A.nzval, 1)
    Tv         = eltype(A.nzval)
    Ti         = eltype(A.colptr)
    colptr     = Array(Ti, n + 1)
    rowval     = Array(Ti, nnz)
    nzval      = Array(Tv, nnz)
    newA       = SparseMatrixCSC{Tv, Ti}(m, n, colptr, rowval, nzval)

    P          = Array(Cint, n)
    Pprime     = Array(Cint, n)

    ccall((:CSR_ReorderMatrix, Sparso.libcsr), Void,
              (Cint, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble},
               Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble},
               Ptr{Cint}, Ptr{Cint}, Bool, Bool, Bool),
               A.m, A.n, pointer(A.colptr), pointer(A.rowval), pointer(A.nzval),
               pointer(newA.colptr), pointer(newA.rowval), pointer(newA.nzval),
               pointer(P), pointer(Pprime), true, true, true)

    # Store the reordered A and P, PPrime into a file
    reordered_file = open(reordered_filename, "w")

    dict = Dict{Any, Int}(Int => 0, Int32 => 1, Int64 => 2, Float64 => 3)
    write(reordered_file, modification_time)
    write(reordered_file, m)
    write(reordered_file, n)
    write(reordered_file, nnz)
    write(reordered_file, dict[Ti])
    write(reordered_file, dict[Tv])
    write(reordered_file, newA.colptr)
    write(reordered_file, newA.rowval)
    write(reordered_file, newA.nzval)
    write(reordered_file, P)
    write(reordered_file, Pprime)
    close(reordered_file)

    return newA, P, Pprime
end

end # module

max_stars_repo_name: IntelLabs/Sparso

max_issues_repo_name: IntelLabs/Sparso

max_forks_repo_name: IntelLabs/Sparso


text:
path = "docs/src/assets/"

using ISA, Plots
𝐶₀ = AMFMtriplet(t->1,t->t^2,0.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(ψ₀; timeaxis=0.0:0.001:10.0)
png(path*"IS_exFM1.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->1,t->t^2,0.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(𝐶₀; timeaxis=0.0:0.001:10.0)
png(path*"IS_exFM2.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->1,t->t^2,0.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(𝐶₀; timeaxis=0.0:0.001:10.0, view="TF",
     left_margin=15Plots.mm, margin=5Plots.mm)
png(path*"IS_exFM3.png")

# ----------------

using ISA, Plots
𝐶₀ = AMFMtriplet(t->1,t->25+20*sin(t),0.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(ψ₀; timeaxis=0.0:0.001:10.0)
png(path*"IS_exFM4.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->1,t->25+20*sin(t),0.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(𝐶₀; timeaxis=0.0:0.001:10.0)
png(path*"IS_exFM5.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->1,t->25+20*sin(t),0.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(𝐶₀; timeaxis=0.0:0.001:10.0, view="TF",
     left_margin=15Plots.mm, margin=5Plots.mm)
png(path*"IS_exFM6.png")

#--------------------

using ISA, Plots
𝐶₀ = AMFMtriplet(t->1,t->exp(0.5t)+10*t,1.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(ψ₀; timeaxis=0.0:0.001:10.0)
png(path*"IS_exFM7.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->1,t->exp(0.5t)+10*t,1.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(𝐶₀; timeaxis=0.0:0.001:10.0)
png(path*"IS_exFM8.png")

using ISA, Plots
𝐶₀ = AMFMtriplet(t->1,t->exp(0.5t)+10*t,1.0)
ψ₀ = AMFMcomp(𝐶₀)
plot(𝐶₀; timeaxis=0.0:0.001:10.0, view="TF",
     left_margin=15Plots.mm, margin=5Plots.mm)
png(path*"IS_exFM9.png")


text:
# Copyright © 2021 Alexander L. Hayes
# Apache 2.0 License

"""
Custom Types
"""

"""RelationalDataset

The RelationalDataset is a `NamedTuple` with positive examples (pos), negative
examples (neg), and facts (facts) about a domain.

Examples
--------

Create an instance of a `RelationalDataset` type:

```julia
using RelationalDatasets

train = RelationalDatasets.RelationalDataset((
    pos=["related(a,b)."],
    neg=["related(b,c)."],
    facts=["child(a,b)."],
))
```

Instances of `RelationalDataset` are returned by `load`:

```julia
using RelationalDatasets

train, test = load("webkb")
```
"""
RelationalDataset = @NamedTuple begin
    pos::Vector{String}
    neg::Vector{String}
    facts::Vector{String}
end


text:
using Test
using Printf
using FileIO

# testset for image

# the tested code
include("../src/image.jl")

# This exercises the constructors
include("tinyImage.jl")

@testset "Image" begin

    testImage = createTestMaskedImage()

    # point out of bounds
    testPoint = PointInMaskedImage( CartesianIndex(-1, -1), testImage)
    @test ! isInBoundsAndSelected(testPoint)

    # point in bounds and in mask
    testPoint = PointInMaskedImage( CartesianIndex(3,2), testImage)

    @test isPointSelectedInMaskedImage(testPoint)
    @test isInBoundsAndSelected(testPoint)

    # point in bounds but not in mask
    testPoint = PointInMaskedImage( CartesianIndex(1,1), testImage)
    @test ! isInBoundsAndSelected(testPoint)

end

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 0d696b90-14b9-11eb-19c6-7110e9eded9e
begin
	using PotentialFlow
	using Plots
	using PlutoUI
end

# ╔═╡ 97a34c00-153a-11eb-0b35-e92e0aafac5e
md"""
### Flow over a cylinder
"""

# ╔═╡ 59e4bc40-14b9-11eb-3670-cdc8fa8f0ded
begin
	X = range(-3.5, 3.5, length=101)
	Y = range(-3.5, 3.5, length=101)
	Z = [x + im*y for y in Y, x in X]
end

# ╔═╡ ede98510-153b-11eb-1e6d-b1e9ad472350
md"""
Doublet strength: $(@bind doublet_strength Slider(0:π/10:2π,show_value=true, default = π))
"""

# ╔═╡ 3d84d310-153f-11eb-1520-17b8af3cf1d0
md"""
X location: $(@bind Δx Slider(-5.0:0.1:5.0,show_value=true,default = 0.0))
"""

# ╔═╡ 7a6a3060-153d-11eb-0208-a9fc8710ab9a
md"""
Y location: $(@bind Δy Slider(-5.0:0.1:5.0,show_value=true,default = 0.0))
"""

# ╔═╡ 2e843420-1538-11eb-14b2-d9ab38fa091b
begin

    fs = Freestreams.Freestream(1.0)
	Δx1 = convert(Float64,Δx)
	Δy1 = convert(Float64,Δy)
	location = (Δx1 + Δy1*im)

	db = Doublets.Doublet(location, doublet_strength)

	ψ_levels = collect(range(-5,5,length=31))

	streamlines(X, Y, (fs,db),
		levels = ψ_levels,
		color = cgrad([:black, :black]),
		ratio = 1,
		legend = false,
		markerstrokealpha = 0,
		grid = false,
		xlim = (-3.5, 3.5),
		ylim = (-3.5, 3.5))

	#circular cylinder
	plot!(cos.(0:0.1:2π), sin.(0:0.1:2π), linestyle = :dash, linewidth=4 )

end

# ╔═╡ Cell order:
# ╟─97a34c00-153a-11eb-0b35-e92e0aafac5e
# ╠═0d696b90-14b9-11eb-19c6-7110e9eded9e
# ╠═59e4bc40-14b9-11eb-3670-cdc8fa8f0ded
# ╠═ede98510-153b-11eb-1e6d-b1e9ad472350
# ╠═3d84d310-153f-11eb-1520-17b8af3cf1d0
# ╠═7a6a3060-153d-11eb-0208-a9fc8710ab9a
# ╠═2e843420-1538-11eb-14b2-d9ab38fa091b

max_stars_repo_path: Flow_over_a_cylinder - Doublet.jl

max_issues_repo_path: Flow_over_a_cylinder - Doublet.jl

max_forks_repo_path: Flow_over_a_cylinder - Doublet.jl


text:
using Gadfly

# Set the random seed to get consistent data
srand(1)

# Number of examples to use
n = 100

# Specify the true value of the variable
true_coeffs = [2; -2; 0.5];

# Generate data
x_data = rand(n, 1) * 5;
x_data_expanded = hcat([x_data .^ i for i in 1 : 3]...)
y_data = x_data_expanded * true_coeffs + 0.5 * rand(n, 1)

p = plot(
  x=x_data, y=y_data, Geom.point,
  Theme(panel_fill=color("white"))
)
# draw(PNG("data.png", 16cm, 12cm), p)


text:
# direct conversion from hexadecimal or binary
Posit8(x::UInt8) = reinterpret(Posit8,x)
Posit16(x::UInt16) = reinterpret(Posit16,x)
Posit32(x::UInt32) = reinterpret(Posit32,x)

Posit8_2(x::UInt32) = reinterpret(Posit8_2,x)
Posit16_2(x::UInt32) = reinterpret(Posit16_2,x)
Posit24_2(x::UInt32) = reinterpret(Posit24_2,x)

# allow also 8bit and 16bit although internally stored as 32bit
# shift bits as conversion to UInt32 adds 0 bits to the left
# but the posit formats assume 0 bits at the right for 8,16bit
Posit8_2(x::UInt8) = Posit8_2(UInt32(x) << 24)
Posit16_2(x::UInt16) = Posit16_2(UInt32(x) << 16)

# same for Posit_1
Posit8_1(x::UInt32) = reinterpret(Posit8_1,x)
Posit16_1(x::UInt32) = reinterpret(Posit16_1,x)
Posit24_1(x::UInt32) = reinterpret(Posit24_1,x)
Posit8_1(x::UInt8) = Posit8_1(UInt32(x) << 24)
Posit16_1(x::UInt16) = Posit16_1(UInt32(x) << 16)


text:
f(x) = sin(x)

function simpsons_rule(a, b, n)
    h = (b - a) / (2*n)
    res = 0.0

    for i in 0:2*n
        if (i==0 || i==2*n)
            res += f(a + i*h)
        elseif (i % 2 == 1)
            res += 4 * f(a + i*h)
        else
            res += 2 * f(a + i*h)
        end
    end

    res * h / 3
end

println("Analytical=", 2)
println("Numerical=", simpsons_rule(0, pi, 100))


text:
function halfwidth(c::Char)
        code = UInt(c)
        if 0xFF01 <= code <= 0xFF5E
                return Char(code - 0xFF01 + 0x0021)
        else
                return c
        end
end

function fullwidth(c::Char)
        code = UInt(c)
        if 0x0021 <= code <= 0x007E
                return Char(code - 0x0021 + 0xFF01)
        else
                return c
        end
end

const CH_NUMBERS = "零一二三四五六七八九十百千万亿"

function normalize_width_numeric(c::Char)
        c = halfwidth(c)
        return isnumeric(c) || any(isequal(c), CH_NUMBERS) ? 'N' : c
end

# normalize width chars & numeric chars & spaces
mask(x::AbstractString) = filter(!isspace, map(normalize_width_numeric, x))

# inverse function of mask
function demask(y::Vector{String}, x::String)
        xchrs = collect(x)
        nx = length(xchrs)
        projections = Vector{Char}(undef, nx)
        ip = 1
        for ix in 1:nx
                if isspace(xchrs[ix]) continue end
                projections[ip] = xchrs[ix]
                ip += 1
        end
        ret = Vector{String}(undef, length(y))
        ip = 1
        for (iy, word) in enumerate(y)
                nw = length(word)
                new_w = Vector{Char}(undef, nw)
                for iw in 1:nw
                        new_w[iw] = projections[ip]
                        ip += 1
                end
                ret[iy] = String(new_w)
        end
        return ret
end

#=
str = "  1999年10月5日 一八九九年 "
println(map(halfwidth, str))
println(normalize_width_numeric_space(str))
=#

isevil(c::Char) = isspace(c) || isnumeric(c) || isletter(c) || ispunct(c)
isevil(s::String) = all(isevil, s)
filter_evil(s::String) = String(filter(!isevil, s))

function shrink_evil(s::String)
        chrs = collect(s)
        nchr = length(chrs)
        ret = Char[]
        sizehint!(ret, nchr)
        i = 1
        while i <= nchr
                while i <= nchr && !isevil(chrs[i])
                        push!(ret, chrs[i])
                        i += 1
                end
                s = i
                while i <= nchr && isevil(chrs[i]) i += 1 end
                t = i
                while s < t && isspace(chrs[s]) s += 1 end
                while s < t && isspace(chrs[t - 1]) t -= 1 end
                if s < t append!(ret, 'E') end
        end
        return String(ret)
end

function expand_evil(y::Vector{String}, x::String)::Vector{String}
        xchrs = collect(x)
        nx = length(xchrs)
        projections = Vector{UnitRange{Int64}}(undef, length(xchrs))
        ix, iy = 1, 1
        while ix <= nx
                while ix <= nx && !isevil(xchrs[ix])
                        projections[iy] = ix:ix
                        ix += 1
                        iy += 1
                end
                s = ix
                while ix <= nx && isevil(xchrs[ix]) ix += 1 end
                t = ix
                while s < t && isspace(xchrs[s]) s += 1 end
                while s < t && isspace(xchrs[t - 1]) t -= 1 end
                if s < t
                        projections[iy] = s:t-1
                        iy += 1
                end
        end
        ret = Vector{String}(undef, length(y))
        offset = 0
        for (i, word) in enumerate(y)
                new_w = Char[]
                sizehint!(new_w, length(word))
                for _ in 1:length(word)
                        offset += 1
                        append!(new_w, x[projections[offset]])
                end
                ret[i] = String(new_w)
        end
        return ret
end

# traditional to simple chinese character
function simplify(c::Char)
    return c
end

simplify(s::String) = String(simplify.(collect(s)))

dir(o...) = joinpath(pathof(KongYiji), "..", "..", "data", o...)

function getid!(d::Dict{String,Ti}, k::String)
        if haskey(d, k)
                return d[k]
        else
                v = length(d) + 1
                d[k] = v
                return v
        end
end

import Base.vec
function vec(d::Dict{String, Ti})
        r = Array{String}(undef, length(d))
        for (k, v) in d r[v] = k end
        return r
end

# Assume the unique of elements of v
function dict(v::Vector{String})
        r = Dict{String, Ti}()
        for (i, k) in enumerate(v) r[k] = i end
        return r
end


text:
# ---
# title: 1562. Find Latest Group of Size M
# id: problem1562
# author: Indigo
# date: 2021-01-22
# difficulty: Medium
# categories: Binary Search
# link: <https://leetcode.com/problems/find-latest-group-of-size-m/description/>
# hidden: true
# ---
#
# Given an array `arr` that represents a permutation of numbers from `1` to `n`.
# You have a binary string of size `n` that initially has all its bits set to
# zero.
#
# At each step `i` (assuming both the binary string and `arr` are 1-indexed)
# from `1` to `n`, the bit at position `arr[i]` is set to `1`. You are given an
# integer `m` and you need to find the latest step at which there exists a group
# of ones of length `m`. A group of ones is a contiguous substring of 1s such
# that it cannot be extended in either direction.
#
# Return _the latest step at which there exists a group of ones of length
# **exactly**_  `m`. _If no such group exists, return_  `-1`.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [3,5,1,2,4], m = 1
#     Output: 4
#     Explanation: Step 1: "00 _1_ 00", groups: ["1"]
#     Step 2: "0010 _1_ ", groups: ["1", "1"]
#     Step 3: " _1_ 0101", groups: ["1", "1", "1"]
#     Step 4: "1 _1_ 101", groups: ["111", "1"]
#     Step 5: "111 _1_ 1", groups: ["11111"]
#     The latest step at which there exists a group of size 1 is step 4.
#
# **Example 2:**
#
#
#
#     Input: arr = [3,1,5,4,2], m = 2
#     Output: -1
#     Explanation: Step 1: "00 _1_ 00", groups: ["1"]
#     Step 2: " _1_ 0100", groups: ["1", "1"]
#     Step 3: "1010 _1_ ", groups: ["1", "1", "1"]
#     Step 4: "101 _1_ 1", groups: ["1", "111"]
#     Step 5: "1 _1_ 111", groups: ["11111"]
#     No group of size 2 exists during any step.
#
#
# **Example 3:**
#
#
#
#     Input: arr = [1], m = 1
#     Output: 1
#
#
# **Example 4:**
#
#
#
#     Input: arr = [2,1], m = 2
#     Output: 2
#
#
#
#
# **Constraints:**
#
#   * `n == arr.length`
#   * `1 <= n <= 10^5`
#   * `1 <= arr[i] <= n`
#   * All integers in `arr` are  **distinct**.
#   * `1 <= m <= arr.length`
#
#
## @lc code=start
using LeetCode

function find_latest_step(arr::Vector{Int}, m::Int)
    pos_dict = Dict{Int, Int}()
    neg_dict = Dict{Int, Int}()
    cnt = 0
    res = -1
    for (idx, num) in enumerate(arr)
        pos_dict[num] = num
        neg_dict[num] = num
        if haskey(pos_dict, num + 1)
            (pos_dict[num + 1] - num == m) && (cnt -= 1)
            tmp = pop!(pos_dict, num + 1)
            pos_dict[num] = tmp
            neg_dict[tmp] = num
        end
        if haskey(neg_dict, num - 1)
            (num - neg_dict[num - 1] == m) && (cnt -= 1)
            tmp = pop!(neg_dict, num - 1)
            neg_dict[num] = tmp
            pos_dict[tmp] = num
        end
        pos_dict[neg_dict[num]], neg_dict[pos_dict[num]] = pos_dict[num], neg_dict[num]
        (pos_dict[neg_dict[num]] - neg_dict[num] + 1 == m) && (cnt += 1)
        (cnt > 0) && (res = idx)
        if pos_dict[neg_dict[num]] != num
            pop!(pos_dict, num)
            pop!(neg_dict, num)
        end
    end
    return res
end

## @lc code=end


text:
using PyCall

@pyimport modeller
@pyimport modeller.automodel as am

env = modeller.environ()
a = am.automodel(env, alnfile="TvLDH-1bdmA.ali",
              knowns="1bdmA", sequence="TvLDH",
              assess_methods=(am.assess[:DOPE],
                              #soap_protein_od.Scorer(),
                              am.assess[:GA341]))
a[:starting_model] = 1
a[:ending_model] = 5
a[:make]()


text:
"""
    StatisticalModel

Abstract supertype for all statistical models.
"""
abstract type StatisticalModel end

"""
    coef(model::StatisticalModel)

Return the coefficients of the model.
"""
function coef end

"""
    coefnames(model::StatisticalModel)

Return the names of the coefficients.
"""
function coefnames end

"""
    coeftable(model::StatisticalModel; level::Real=0.95)

Return a table with coefficients and related statistics of the model.
`level` determines the level for confidence intervals (by default, 95%).

The returned `CoefTable` object implements the
[Tables.jl](https://github.com/JuliaData/Tables.jl/) interface, and can be
converted e.g. to a `DataFrame` via `using DataFrames; DataFrame(coeftable(model))`.
"""
function coeftable end

"""
    confint(model::StatisticalModel; level::Real=0.95)

Compute confidence intervals for coefficients, with confidence level `level` (by default 95%).
"""
function confint end

"""
    deviance(model::StatisticalModel)

Return the deviance of the model relative to a reference, which is usually when applicable
the saturated model. It is equal, *up to a constant*, to ``-2 \\log L``, with ``L``
the likelihood of the model.
"""
function deviance end

"""
    islinear(model::StatisticalModel)

Indicate whether the model is linear.
"""
function islinear end

"""
    nulldeviance(model::StatisticalModel)

Return the deviance of the null model, that is the one including only the intercept.
"""
function nulldeviance end

"""
    loglikelihood(model::StatisticalModel)
    loglikelihood(model::StatisticalModel, observation)

Return the log-likelihood of the model.

With an `observation` argument, return the contribution of `observation` to the
log-likelihood of `model`.

If `observation` is a `Colon`, return a vector of each observation's contribution
to the log-likelihood of the model. In other words, this is the vector of the
pointwise log-likelihood contributions.

In general, `sum(loglikehood(model, :)) == loglikelihood(model)`.
"""
function loglikelihood end

"""
    nullloglikelihood(model::StatisticalModel)

Return the log-likelihood of the null model corresponding to `model`.
This is usually the model containing only the intercept.
"""
function nullloglikelihood end

"""
    score(model::StatisticalModel)

Return the score of the model, that is the gradient of the
log-likelihood with respect to the coefficients.
"""
function score end

"""
    nobs(model::StatisticalModel)

Return the number of independent observations on which the model was fitted. Be careful
when using this information, as the definition of an independent observation may vary
depending on the model, on the format used to pass the data, on the sampling plan
(if specified), etc.
"""
function nobs end

"""
    dof(model::StatisticalModel)

Return the number of degrees of freedom consumed in the model, including
when applicable the intercept and the distribution's dispersion parameter.
"""
function dof end

"""
    mss(model::StatisticalModel)

Return the model sum of squares.
"""
function mss end

"""
    rss(model::StatisticalModel)

Return the residual sum of squares of the model.
"""
function rss end

"""
    informationmatrix(model::StatisticalModel; expected::Bool = true)

Return the information matrix of the model. By default the Fisher information matrix
is returned, while the observed information matrix can be requested with `expected = false`.
"""
function informationmatrix end

"""
    stderror(model::StatisticalModel)

Return the standard errors for the coefficients of the model.
"""
function stderror end

"""
    vcov(model::StatisticalModel)

Return the variance-covariance matrix for the coefficients of the model.
"""
function vcov end

"""
    weights(model::StatisticalModel)

Return the weights used in the model.
"""
function weights end

"""
    isfitted(model::StatisticalModel)

Indicate whether the model has been fitted.
"""
function isfitted end

"""
Fit a statistical model.
"""
function fit end

"""
Fit a statistical model in-place.
"""
function fit! end

"""
    aic(model::StatisticalModel)

Akaike's Information Criterion, defined as ``-2 \\log L + 2k``, with ``L`` the likelihood
of the model, and `k` its number of consumed degrees of freedom
(as returned by [`dof`](@ref)).
"""
aic(model::StatisticalModel) = -2loglikelihood(model) + 2dof(model)

"""
    aicc(model::StatisticalModel)

Corrected Akaike's Information Criterion for small sample sizes (Hurvich and Tsai 1989),
defined as ``-2 \\log L + 2k + 2k(k-1)/(n-k-1)``, with ``L`` the likelihood of the model,
``k`` its number of consumed degrees of freedom (as returned by [`dof`](@ref)),
and ``n`` the number of observations (as returned by [`nobs`](@ref)).
"""
function aicc(model::StatisticalModel)
    k = dof(model)
    n = nobs(model)
    -2loglikelihood(model) + 2k + 2k*(k+1)/(n-k-1)
end

"""
    bic(model::StatisticalModel)

Bayesian Information Criterion, defined as ``-2 \\log L + k \\log n``, with ``L``
the likelihood of the model,  ``k`` its number of consumed degrees of freedom
(as returned by [`dof`](@ref)), and ``n`` the number of observations
(as returned by [`nobs`](@ref)).
"""
bic(model::StatisticalModel) = -2loglikelihood(model) + dof(model)*log(nobs(model))

"""
    r2(model::StatisticalModel)
    r²(model::StatisticalModel)

Coefficient of determination (R-squared).

For a linear model, the R² is defined as ``ESS/TSS``, with ``ESS`` the explained sum of squares
and ``TSS`` the total sum of squares.
"""
function r2 end

const r² = r2

"""
    adjr2(model::StatisticalModel)
    adjr²(model::StatisticalModel)

Adjusted coefficient of determination (adjusted R-squared).

For linear models, the adjusted R² is defined as ``1 - (1 - (1-R^2)(n-1)/(n-p))``, with ``R^2``
the coefficient of determination, ``n`` the number of observations, and ``p`` the number of
coefficients (including the intercept). This definition is generally known as the Wherry Formula I.
"""
function adjr2 end

const adjr² = adjr2


text:
###### This file is part of the MendelImpute.jl package.
###### It contains relevant code to search breakpoints.
###### Currently code for double breakpoint search is broken and commented out.

"""
    continue_haplotype(X, compressed_Hunique, window, happair_prev, happair_next)

Searches the breakpoint between `happair_prev` and `happair_next`, if there is
one.

# Arguments:
- `X`: Genotype vector spanning 2 windows.
- `compressed_Hunique`: A `CompressedHaplotypes` keeping track of unique
    haplotypes for each window and some other information
- `window`: The current window being considered. It is the 2nd window of `X`.
- `happair_prev`: Haplotype pair for first window
- `happair_next`: Haplotype pair for second window

# Optional arguments:
- `phased`: Boolean indicating whether `happair_prev` and `happair_next` have
    been phased. If `false`, will try 2 different orientations.
- `search_double_bkpts`: If `true`, will search double breakpoints. If false
    the output `bkpt = (-2, -2)`

# Output
- `happair_next`: If `phase = false`, this is equal to `happair_next` for input.
    Otherwise it may be flipped
- `bkpt`: Tuple of integer indicating optimal breakpoint starting from the
    previous window. `-1` indicates no breakpoints and `-2` indicates double
    breakpoint.
"""
function continue_haplotype(
    X::AbstractVector,
    compressed_Hunique::CompressedHaplotypes,
    window::Int,
    happair_prev::Tuple,
    happair_next::Tuple;
    phased::Bool = false,
    search_double_bkpts::Bool = false
    )

    # indices for complete reference panel
    i, j = happair_prev
    k, l = happair_next

    # both strands match
    if i == k && j == l
        return (k, l), (-1, -1)
    end
    if i == l && j == k && !phased
        return (l, k), (-1, -1)
    end

    # get unique haplotypes with only typed snps
    overlap = compressed_Hunique.overlap
    if overlap
        # create views on non-overlapping regions
        rprev = nonoverlap_range(compressed_Hunique, window - 1)
        rcurr = nonoverlap_range(compressed_Hunique, window)
        Hprev = view(compressed_Hunique.CW_typed[window - 1].uniqueH, rprev, :)
        Hcurr = view(compressed_Hunique.CW_typed[window].uniqueH, rcurr, :)
    else
        Hprev = compressed_Hunique.CW_typed[window - 1].uniqueH
        Hcurr = compressed_Hunique.CW_typed[window].uniqueH
    end

    # only one strand matches
    if i == k && j ≠ l
        iu = complete_idx_to_unique_typed_idx(i, window - 1, compressed_Hunique)
        ku = complete_idx_to_unique_typed_idx(k, window, compressed_Hunique)
        ju1 = complete_idx_to_unique_typed_idx(j, window - 1, compressed_Hunique)
        ju2 = complete_idx_to_unique_typed_idx(j, window, compressed_Hunique)
        lu1 = complete_idx_to_unique_typed_idx(l, window - 1, compressed_Hunique)
        lu2 = complete_idx_to_unique_typed_idx(l, window, compressed_Hunique)

        # lazy concatenation
        s1  = ApplyArray(vcat, view(Hprev, :, iu),  view(Hcurr, :, ku))
        s21 = ApplyArray(vcat, view(Hprev, :, ju1), view(Hcurr, :, ju2))
        s22 = ApplyArray(vcat, view(Hprev, :, lu1), view(Hcurr, :, lu2))

        breakpt, errors = search_breakpoint(X, s1, s21, s22)
        return (k, l), (-1, breakpt)
    elseif i == l && j ≠ k && !phased
        iu = complete_idx_to_unique_typed_idx(i, window - 1, compressed_Hunique)
        lu = complete_idx_to_unique_typed_idx(l, window, compressed_Hunique)
        ju1 = complete_idx_to_unique_typed_idx(j, window - 1, compressed_Hunique)
        ju2 = complete_idx_to_unique_typed_idx(j, window, compressed_Hunique)
        ku1 = complete_idx_to_unique_typed_idx(k, window - 1, compressed_Hunique)
        ku2 = complete_idx_to_unique_typed_idx(k, window, compressed_Hunique)

        # lazy concatenation
        s1  = ApplyArray(vcat, view(Hprev, :, iu),  view(Hcurr, :, lu))
        s21 = ApplyArray(vcat, view(Hprev, :, ju1), view(Hcurr, :, ju2))
        s22 = ApplyArray(vcat, view(Hprev, :, ku1), view(Hcurr, :, ku2))

        breakpt, errors = search_breakpoint(X, s1, s21, s22)
        return (l, k), (-1, breakpt)
    elseif j == k && i ≠ l && !phased
        ju = complete_idx_to_unique_typed_idx(j, window - 1, compressed_Hunique)
        ku = complete_idx_to_unique_typed_idx(k, window, compressed_Hunique)
        iu1 = complete_idx_to_unique_typed_idx(i, window - 1, compressed_Hunique)
        iu2 = complete_idx_to_unique_typed_idx(i, window, compressed_Hunique)
        lu1 = complete_idx_to_unique_typed_idx(l, window - 1, compressed_Hunique)
        lu2 = complete_idx_to_unique_typed_idx(l, window, compressed_Hunique)

        # lazy concatenation
        s1  = ApplyArray(vcat, view(Hprev, :, ju),  view(Hcurr, :, ku))
        s21 = ApplyArray(vcat, view(Hprev, :, iu1), view(Hcurr, :, iu2))
        s22 = ApplyArray(vcat, view(Hprev, :, lu1), view(Hcurr, :, lu2))

        breakpt, errors = search_breakpoint(X, s1, s21, s22)
        return (l, k), (breakpt, -1)
    elseif j == l && i ≠ k
        ju = complete_idx_to_unique_typed_idx(j, window - 1, compressed_Hunique)
        lu = complete_idx_to_unique_typed_idx(l, window, compressed_Hunique)
        iu1 = complete_idx_to_unique_typed_idx(i, window - 1, compressed_Hunique)
        iu2 = complete_idx_to_unique_typed_idx(i, window, compressed_Hunique)
        ku1 = complete_idx_to_unique_typed_idx(k, window - 1, compressed_Hunique)
        ku2 = complete_idx_to_unique_typed_idx(k, window, compressed_Hunique)

        # lazy concatenation
        s1  = ApplyArray(vcat, view(Hprev, :, ju),  view(Hcurr, :, lu))
        s21 = ApplyArray(vcat, view(Hprev, :, iu1), view(Hcurr, :, iu2))
        s22 = ApplyArray(vcat, view(Hprev, :, ku1), view(Hcurr, :, ku2))

        breakpt, errors = search_breakpoint(X, s1, s21, s22)
        return (k, l), (breakpt, -1)
    end

    # both strand mismatch
    if search_double_bkpts
        iu1 = complete_idx_to_unique_typed_idx(i, window - 1, compressed_Hunique)
        iu2 = complete_idx_to_unique_typed_idx(i, window, compressed_Hunique)
        ku1 = complete_idx_to_unique_typed_idx(k, window - 1, compressed_Hunique)
        ku2 = complete_idx_to_unique_typed_idx(k, window, compressed_Hunique)
        ju1 = complete_idx_to_unique_typed_idx(j, window - 1, compressed_Hunique)
        ju2 = complete_idx_to_unique_typed_idx(j, window, compressed_Hunique)
        lu1 = complete_idx_to_unique_typed_idx(l, window - 1, compressed_Hunique)
        lu2 = complete_idx_to_unique_typed_idx(l, window, compressed_Hunique)

        # lazy concatenation
        s11 = ApplyArray(vcat, view(Hprev, :, iu1), view(Hcurr, :, iu2))
        s12 = ApplyArray(vcat, view(Hprev, :, ku1), view(Hcurr, :, ku2))
        s21 = ApplyArray(vcat, view(Hprev, :, ju1), view(Hcurr, :, ju2))
        s22 = ApplyArray(vcat, view(Hprev, :, lu1), view(Hcurr, :, lu2))

        breakpt, errors = search_breakpoint(X, s11, s12, s21, s22)
        return (k, l), breakpt

        # Always assume double bkpt searches are phased.
        # breakpt1, errors1 = search_breakpoint(X, H, (i, k), (j, l))
        # breakpt2, errors2 = search_breakpoint(X, H, (i, l), (j, k))
        # if errors1 < errors2
        #     return (k, l), breakpt1
        # else
        #     return (l, k), breakpt2
        # end
    else
        return (k, l), (-2, -2)
    end
end

"""
    search_single_breakpoint(X, s1, s21, s22)

Find the optimal break point between s21 and s22 in configuration
s1 | s21
s1 | s22
"""
function search_breakpoint(
    X::AbstractVector,
    s1::AbstractVector,
    s21::AbstractVector,
    s22::AbstractVector,
    )

    n = length(X)
    length(s1) == length(s21) == length(s22) == n ||error("search_breakpoint:" *
        " all vectors should have same length but length X = $n, s1 =" *
        " $(length(s1)), s21 = $(length(s21)), s22 = $(length(s22))")

    # count number of errors if second haplotype is all from s22
    errors = 0
    @inbounds for pos in 1:n
        if !ismissing(X[pos])
            errors += X[pos] ≠ s1[pos] + s22[pos]
        end
    end
    bkpt_optim, err_optim = 0, errors

    # quick return if perfect match
    err_optim == 0 && return 0, 0

    # extend haplotype s21 position by position
    @inbounds for bkpt in 1:n
        if !ismissing(X[bkpt]) && s21[bkpt] ≠ s22[bkpt]
            errors -= X[bkpt] ≠ s1[bkpt] + s22[bkpt]
            errors += X[bkpt] ≠ s1[bkpt] + s21[bkpt]
            if errors :: Int < err_optim
                bkpt_optim, err_optim = bkpt, errors
                # quick return if perfect match
                err_optim == 0 && return bkpt_optim, err_optim :: Int
            end
        end
    end

    return bkpt_optim, err_optim :: Int
end

"""
    search_breakpoint(X, s11, s12, s21, s22)

Find the optimal break point between s11 and s12 as well as between s21 and s22
in configuration
s11 | s21
s12 | s22
"""
function search_breakpoint(
    X::AbstractVector,
    s11::AbstractVector,
    s12::AbstractVector,
    s21::AbstractVector,
    s22::AbstractVector,
    )
    n = length(X)
    length(s11) == length(s12) == length(s21) == length(s22) == n ||
        error("search_breakpoint: all vectors should have same length but " *
        "length X = $n, s11 = $(length(s11)), s12 = $(length(s12)), s21 = " *
        "$(length(s21)), s22 = $(length(s22))")

    err_optim   = typemax(Int)
    bkpts_optim = (0, 0)

    # search over all combintations of break points in two strands
    @inbounds for bkpt1 in 0:n

        # count number of errors if second haplotype is all from s22
        errors = 0
        for pos in 1:bkpt1
            if !ismissing(X[pos])
                errors += X[pos] ≠ s11[pos] + s22[pos]
            end
        end
        for pos in (bkpt1 + 1):length(X)
            if !ismissing(X[pos])
                errors += X[pos] ≠ s12[pos] +s22[pos]
            end
        end
        if errors :: Int < err_optim
            err_optim = errors
            bkpts_optim = (bkpt1, 0)

            # quick return if perfect match
            err_optim == 0 && return bkpts_optim, err_optim :: Int
        end

        # extend haplotype s21 position by position
        for bkpt2 in 1:bkpt1
            if !ismissing(X[bkpt2]) && s21[bkpt2] ≠ s22[bkpt2]
                errors -= X[bkpt2] ≠ s11[bkpt2] + s22[bkpt2]
                errors += X[bkpt2] ≠ s11[bkpt2] + s21[bkpt2]
                if errors :: Int < err_optim
                    err_optim = errors
                    bkpts_optim = (bkpt1, bkpt2)
                end
            end
        end
        for bkpt2 in (bkpt1 + 1):length(X)
            if !ismissing(X[bkpt2]) && s21[bkpt2] ≠ s22[bkpt2]
                errors -= X[bkpt2] ≠ s12[bkpt2] + s22[bkpt2]
                errors += X[bkpt2] ≠ s12[bkpt2] + s21[bkpt2]
                if errors :: Int < err_optim
                    err_optim = errors
                    bkpts_optim = (bkpt1, bkpt2)
                    # quick return if perfect match
                    err_optim == 0 && return bkpts_optim, err_optim :: Int
                end
            end
        end
    end

    return bkpts_optim, err_optim :: Int
end


text:
export Mud,
        DrillFluid

abstract type Mud end

struct DrillFluid <: Mud
    Liquids::Vector{Fluid}
    Solids::Vector{Solid}

    nLiq::Integer
    nSol::Integer

    fL::Vector{Float64}
    sum_fL::Float64
    fS::Vector{Float64}
    sum_fS::Float64

    # wt::Vector{Float64}
    wtL::Vector{Float64}
    wtS::Vector{Float64}

    function DrillFluid(Flist::Vector{<:Fluid},Slist::Vector{<:Solid};fs::AbstractVector,normalize=false)

        nL = Integer(length(Flist))
        nS = Integer(length(Slist))

        missing_pos = findall(ismissing,fs)

        if !isempty(missing_pos)
            if length(missing_pos)>1
                @error "Too many volumetric proportions field left open. Package can complete the mixture with only one component at a time."
            end
            total_frac=sum(skipmissing(fs))
            if total_frac>1.0
                @error "Normalization only available if no missing values are provided."
            end

            fs[missing_pos[1]]=1.0-total_frac
            total_frac=1.0
        else
            total_frac=sum(fs)
            if total_frac != 1.0
                if normalize
                    @warn "Total volumetric fractions greater than 100%, normalizing volume to crrect. If that is not desired pass `normalize` as false."

                    fs=fs./total_frac
                    total_frac=1.0

                else
                    @error "Total volume fraction provided greater than 100%, if the values provided are total volumes change the keyword argument `normalize` to true and the values will be normalized."
                end
            end
        end

        fL=fs[1:length(Flist)]
        sum_fL=sum(fL)
        fS=fs[length(Flist)+1:end]
        sum_fS=sum(fS)

        Lrhos=[rho.ρ for rho ∈ Flist]
        Srhos=[rho.ρ for rho ∈ Slist]

        mLs=[rho*f for (rho,f) ∈ zip(Lrhos,fL)]
        mSs=[rho*f for (rho,f) ∈ zip(Srhos,fS)]

        mT=sum(mLs)+sum(mSs)

        wtL=mLs./mT
        wtS=mSs./mT

        new(
            Flist,
            Slist,
            nL,
            nS,
            fL,
            sum_fL,
            fS,
            sum_fS,
            wtL,
            wtS
        )

    end

    # function DrillFluid(Flist,Slist;fs=missing,wts::Vector{Float64})
    #     #TODO: completar esse trecho com a primeira parte dos modelos de mistura... implementar mais de uma versão, um com o uso de frações de massa, outro com volume. Por enquanto não usar

    # end

end

## Functions

# density of discret phases

function rho(P,T,fluid::DrillFluid)
    fluid_list=fluid.Liquids
    solids_list=fluid.Solids

    fL=fluid.fL
    fS=fluid.fS

    rho_l=[rho(P,T,f) for f ∈ fluid_list]
    rho_s=[rho(P,T,s) for s ∈ solids_list]

    rho_eq_l=sum(fL'*rho_l)
    rho_eq_s=sum(fS'*rho_s)

    f_TL=sum(fL)

    rho_T=rho_eq_l+rho_eq_s

    return rho_T
end


text:
using DFTK
using JLD2
using JSON3

# DFTK
setup_threading()

a = 10.26  # Silicon lattice constant in Bohr
lattice = a / 2 * [[0 1 1.];
                   [1 0 1.];
                   [1 1 0.]]
Si = ElementPsp(:Si, psp=load_psp("hgh/lda/Si-q4"))
atoms = [Si => [ones(3)/8, -ones(3)/8]]

model = model_LDA(lattice, atoms)
basis = PlaneWaveBasis(model, 10.0; kgrid=[1, 1, 1])

callback = DFTK.ScfSaveCheckpoints("checkpointfile.jld2") ∘ DFTK.ScfDefaultCallback()
scfres = self_consistent_field(basis, tol=1e-2, callback=callback, mixing=KerkerMixing())
compute_forces_cart(scfres)

# JSON
JSON3.read(JSON3.write(Dict("a" => 1.0)))["a"]

max_stars_repo_name: mfherbst/ase-dftk

max_issues_repo_name: mfherbst/ase-dftk

max_forks_repo_name: mfherbst/asedftk