text:
Base.@kwdef struct WellIndexSolver{T}
    J⁻¹::Vector{T}
    J⁻¹min::Vector{T}
    J⁻¹max::Vector{T}
    B::SparseMatrixCSC{T, Int}
    W::SparseMatrixCSC{T, Int}
    lb::BitVector
    ub::BitVector
    ΔPmin::Vector{T}
    ΔPmax::Vector{T}
    lbviews::VectorViewBool
    ubviews::VectorViewBool
end

Base.@kwdef struct PresTargetTerm{T} <: AbstractTargetTerm{T}
    Pobs::Vector{T}
    Wobs::Vector{T}
    ΔP::Vector{T}
    g::Matrix{T}
    gviews::Vector{ColumnSlice{T}}
end

Base.@kwdef struct PresBoundTerm{S, T} <: AbstractTargetTerm{T}
    Pobs::Vector{T}
    Wobs::Vector{T}
    ΔP::Vector{T}
    g::Matrix{T}
    gviews::Vector{ColumnSlice{T}}
end

Base.@kwdef struct PbhpTargetTerm{T} <: AbstractTargetTerm{T}
    Pobs::Vector{T}
    Wobs::Vector{T}
    ΔP::Vector{T}
    g::Matrix{T}
    gviews::Vector{ColumnSlice{T}}
    solver::WellIndexSolver{T}
end

Base.@kwdef struct PinjTargetTerm{T} <: AbstractTargetTerm{T}
    Pobs::Vector{T}
    Wobs::Vector{T}
    ΔP::Vector{T}
    g::Matrix{T}
    gλ::Matrix{T}
    gviews::Vector{ColumnSlice{T}}
    gλviews::Vector{ColumnSlice{T}}
    solver::WellIndexSolver{T}
    J⁻¹min_h::Vector{T}
    J⁻¹max_h::Vector{T}
    λviews::CartesianView{T}
end

Base.@kwdef struct L2TargetTerm{T} <: AbstractTargetTerm{T}
    x::Vector{T}
    g::Vector{T}
    αₓ::Vector{T}
end

Base.@kwdef struct WellIndexViews{S, T<:AbstractFloat}
    xviews::VectorView{T}
    pviews::Vector{RowRange{T}}
    yviews::CartesianView{T}
    λviews::CartesianView{T}
end

Base.@kwdef struct TargetFunction{T<:AbstractFloat}
    terms::@NamedTuple begin
        Pres::PresTargetTerm{T}
        Pbhp::PbhpTargetTerm{T}
        Pinj::PinjTargetTerm{T}
        Pmin::PresBoundTerm{:Pmin, T}
        Pmax::PresBoundTerm{:Pmax, T}
        L2::L2TargetTerm{T}
    end
    wviews::@NamedTuple begin
        Jp::WellIndexViews{:Jp, T}
        Jinj::WellIndexViews{:Jinj, T}
        Gw::WellIndexViews{:Gw, T}
    end
    prob::NonlinearProblem{T}
end

function WellIndexSolver{T}(df::AbstractDataFrame, Qobs, Wobs, shape, α) where {T}

    # Indices of fixing intervals
    idx = @with df begin
        rng = UnitRange.(:Jstart, :Jstop)
        getindex.(Ref(LinearIndices(shape)), :Istart, rng)
    end

    # Fill in system matrix of rates (Q ⋅ J⁻¹ = ΔP)
    Q = spzeros(T, length(Qobs), nrow(df))
    @_ map(_1[_2] .= Qobs[_2], eachcol(Q), idx)

    # Diagonal weight matrix
    Dw = Diagonal(replace(Wobs, missing => zero(T)))

    # Matrix for calculating reciprocal PI (Ĵ⁻¹ = B ⋅ ΔP)
    B = dropzeros(Q' * Dw * Q)
    map!(inv, B.nzval, B.nzval)
    B = dropzeros(B * Q' * Dw)

    # Matrix for calculating weighted sum of squared residuals (β = ΔPᵀ ⋅ W ⋅ ΔP)
    lmul!(T(2) * α, Dw)
    W = I - Q * B
    W = dropzeros(W' * Dw * W)

    # Parameters for calculating PI
    params = (
        J⁻¹ = Array{T}(undef, nrow(df)),
        lb = falses(nrow(df)),
        ub = falses(nrow(df)),
        J⁻¹min = inv.(df.Max_value),
        J⁻¹max = inv.(df.Min_value),
        ΔPmin = Q * inv.(df.Max_value),
        ΔPmax = Q * inv.(df.Min_value),
    )

    # References to flags that indicate that PI (ΔP) is out of bounds [lb, ub]
    nums = Array{Int}(undef, length(Qobs))
    @_ map(nums[_1] .= _2, idx, axes(df, 1))
    lbviews = view(params.lb, nums)
    ubviews = view(params.ub, nums)

    return WellIndexSolver{T}(; B, W, lbviews, ubviews, params...)
end

function WellIndexViews{S, T}(df::AbstractDataFrame, prob::NonlinearProblem{T}, J⁻¹) where {S, T}

    # PI matrix
    J = getfield(prob.params, S)

    # References to fitting elements of the vector of reciprocal PI
    df_S = getparams(df, Val(S))
    xviews = view(J⁻¹, axes(df_S, 1)[.!df_S.Ignore])

    # References to fitting elements of the PI matrix
    df_S = getparams(df, Val(S), Val(:var))
    pviews, yviews = @with df_S begin
        view.(Ref(J), :Istart, UnitRange.(:Jstart, :Jstop)),
        view(J, CartesianIndex.(:Istart, :Jstart))
    end

    # References to injection efficiency factor (needed only for S == Jinj)
    idx = replace(df_S.Link, missing => 1)
    df_λ = view(getparams(df, Val(:λ)), idx, :)
    λviews = @with df_λ begin
        view(prob.params.λ, CartesianIndex.(:Istart, :Jstart))
    end

    return WellIndexViews{S, T}(; xviews, pviews, yviews, λviews)
end

function PresTargetTerm{T}(df::AbstractDataFrame, idx, α) where {T}
    Pobs, Wobs = @with df begin
        :Pres[vec(idx)], :Wres[vec(idx)]
    end
    params = (
        Pobs = replace(Pobs, missing => zero(T)),
        Wobs = T(2) .* α .* replace(Wobs, missing => zero(T)),
        ΔP = Array{T}(undef, length(idx)),
        g = Array{T}(undef, size(idx)),
    )
    gviews = collect(eachcol(params.g))
    PresTargetTerm{T}(; gviews, params...)
end

getpobs(df, idx, ::Val{:Pmin}) = @with(df, :Pres_min[idx])
getpobs(df, idx, ::Val{:Pmax}) = @with(df, :Pres_max[idx])

function PresBoundTerm{S, T}(df::AbstractDataFrame, idx, α) where {S, T}
    Pobs = getpobs(df, vec(idx), Val(S))
    params = @with df begin (
        Pobs = replace(Pobs, missing => zero(T)),
        Wobs = T(2) .* α .* .!ismissing.(Pobs),
        ΔP = Array{T}(undef, length(idx)),
        g = Array{T}(undef, size(idx)),
    ) end
    gviews = collect(eachcol(params.g))
    PresBoundTerm{S, T}(; gviews, params...)
end

function PbhpTargetTerm{T}(df_rates::AbstractDataFrame, df_params::AbstractDataFrame, idx, α) where {T}

    # List of all values of PI
    df_view = getparams(df_params, Val(:Gw))
    Pobs, Wobs, Qobs = @with df_rates begin
        :Pbhp_prod[vec(idx)],
        :Wbhp_prod[vec(idx)],
        (:Qliq ./ :Total_mobility)[vec(idx)]
    end

    # Solver for PI
    solver = WellIndexSolver{T}(df_view, Qobs, Wobs, size(idx), α)

    # Parameters for calculating PI
    params = (
        Pobs = replace(Pobs, missing => zero(T)),
        Wobs = T(2) .* α .* replace(Wobs, missing => zero(T)),
        ΔP = Array{T}(undef, length(idx)),
        g = Array{T}(undef, size(idx)),
    )
    gviews = collect(eachcol(params.g))

    PbhpTargetTerm{T}(; gviews, solver, params...)
end

function PinjTargetTerm{T}(df_rates::AbstractDataFrame, df_params::AbstractDataFrame, prob::NonlinearProblem{T}, idx, α) where {T}

    # List of all values of PI
    df_view = getparams(df_params, Val(:Jinj))
    Pobs, Wobs, Qobs = @with df_rates begin
        :Pbhp_inj[vec(idx)],
        :Wbhp_inj[vec(idx)],
        .-:Qinj[vec(idx)]
    end

    # Solver for PI
    solver = WellIndexSolver{T}(df_view, Qobs, Wobs, size(idx), α)

    # Parameters needed for calculating PI
    params = (
        Pobs = replace(Pobs, missing => zero(T)),
        Wobs = T(2) .* α .* replace(Wobs, missing => zero(T)),
        ΔP = Array{T}(undef, length(idx)),
        g = Array{T}(undef, size(idx)),
        gλ = Array{T}(undef, size(idx)),
        J⁻¹min_h = copy(solver.J⁻¹min),
        J⁻¹max_h = copy(solver.J⁻¹max),
    )
    gviews = collect(eachcol(params.g))
    gλviews = collect(eachcol(params.gλ))

    # References to injection efficiency factors
    df_λ = view(getparams(df_params, Val(:λ)), df_view.Link, :)
    λviews = @with df_λ begin
        view(prob.params.λ, CartesianIndex.(:Istart, :Jstart))
    end

    PinjTargetTerm{T}(; gviews, gλviews, solver, λviews, params...)
end

function L2TargetTerm{T}(df::AbstractDataFrame, x, α) where {T}
    df_view = getparams(df, Val(:tanks), Val(:var))
    params = (
        αₓ = T(2) .* α .* replace(df_view.alpha, missing => zero(T)),
        g = Array{T}(undef, length(x)),
    )
    return L2TargetTerm{T}(; x, params...)
end

function TargetFunction{T}(df_rates::AbstractDataFrame, df_params::AbstractDataFrame, prob::NonlinearProblem{T}, fset::FittingSet{T}, α) where {T}

    # Linear indexing
    Nt = size(prob.C, 2)
    Nd = length(prob.pviews)
    idx = LinearIndices((Nd, Nt))'

    # Separate terms of objective function
    terms = @with df_rates begin (
        Pres = PresTargetTerm{T}(df_rates, idx, α["alpha_resp"]),
        Pbhp = PbhpTargetTerm{T}(df_rates, df_params, idx, α["alpha_bhp"]),
        Pinj = PinjTargetTerm{T}(df_rates, df_params, prob, idx, α["alpha_inj"]),
        Pmin = PresBoundTerm{^(:Pmin), T}(df_rates, idx, α["alpha_lb"]),
        Pmax = PresBoundTerm{^(:Pmax), T}(df_rates, idx, α["alpha_ub"]),
        L2 = L2TargetTerm{T}(df_params, fset.cache.ybuf, α["alpha_l2"]),
    ) end

    # References to well parameters
    wviews = (
        Jp = WellIndexViews{:Jp, T}(df_params, prob, terms.Pbhp.solver.J⁻¹),
        Jinj = WellIndexViews{:Jinj, T}(df_params, prob, terms.Pinj.solver.J⁻¹),
        Gw = WellIndexViews{:Gw, T}(df_params, prob, terms.Pbhp.solver.J⁻¹)
    )

    TargetFunction{T}(; terms, wviews, prob)
end

function update_targ!(targ::TargetFunction)
    # FIXED: Using 'map' instead of 'for' preserves 'type-stability'
    map(term -> update_term!(term, targ.prob), values(targ.terms))
    return targ
end

function update_term!(term::PresTargetTerm, prob::NonlinearProblem)
    @unpack Pcalc = prob.params
    @unpack Wobs, Pobs, ΔP, g = term
    @turbo for i = 1:length(Pcalc)
        ΔP[i] = Pcalc[i] - Pobs[i]
        g[i] = Wobs[i] * ΔP[i]
    end
    return term
end

compare(Pcalc, Pobs, ::Val{:Pmin}) = Pcalc < Pobs
compare(Pcalc, Pobs, ::Val{:Pmax}) = Pcalc > Pobs

function update_term!(term::PresBoundTerm{S}, prob::NonlinearProblem) where {S}
    @unpack Pcalc = prob.params
    @unpack Pobs, Wobs, ΔP, g = term
    @turbo for i = 1:length(Pcalc)
        W = compare(Pcalc[i], Pobs[i], Val(S))
        ΔP[i] = Pcalc[i] - Pobs[i]
        g[i] = W * Wobs[i] * ΔP[i]
    end
    return term
end

function solve!(g, solver::WellIndexSolver, ΔP)
    @unpack J⁻¹, J⁻¹min, J⁻¹max, B, W, lb, ub = solver
    mul!(J⁻¹, B, ΔP)
    @turbo for i = 1:length(J⁻¹)
        lb[i] = J⁻¹[i] < J⁻¹min[i]
        ub[i] = J⁻¹[i] > J⁻¹max[i]
        J⁻¹[i] = clamp(J⁻¹[i], J⁻¹min[i], J⁻¹max[i])
    end
    mul!(g, W, ΔP)
    return g
end

function update_term!(term::PbhpTargetTerm, prob::NonlinearProblem)
    @unpack Pcalc = prob.params
    @unpack Pobs, Wobs, ΔP, g = term
    @unpack W, ΔPmin, ΔPmax, lbviews, ubviews = term.solver
    @turbo for i = 1:length(Pcalc)
        ΔP[i] = Pcalc[i] - Pobs[i]
    end
    solve!(vec(g), term.solver, ΔP)
    @inbounds @simd for i = 1:length(ΔP)
        if lbviews[i]
            ΔP[i] -= ΔPmin[i]
            g[i] = Wobs[i] * ΔP[i]
        elseif ubviews[i]
            ΔP[i] -= ΔPmax[i]
            g[i] = Wobs[i] * ΔP[i]
        end
    end
    return term
end

function update_term!(term::PinjTargetTerm{T}, prob::NonlinearProblem) where {T}
    @unpack Pcalc, λ = prob.params
    @unpack Pobs, Wobs, ΔP, g = term
    @unpack gλ, J⁻¹min_h, J⁻¹max_h, λviews = term
    @unpack W, ΔPmin, ΔPmax, lbviews, ubviews, J⁻¹min, J⁻¹max = term.solver
    @turbo for i = 1:length(Pcalc)
        ΔP[i] = Pcalc[i] - Pobs[i]
    end
    @inbounds @simd for i = 1:length(λviews)
        J⁻¹min[i] = λviews[i] * J⁻¹min_h[i]
        J⁻¹max[i] = λviews[i] * J⁻¹max_h[i]
    end
    solve!(vec(g), term.solver, ΔP)
    @inbounds @simd for i = 1:length(ΔP)
        if lbviews[i]
            ΔP[i] -= λ[i] * ΔPmin[i]
            g[i] = Wobs[i] * ΔP[i]
            gλ[i] = -ΔPmin[i] * g[i]
        elseif ubviews[i]
            ΔP[i] -= λ[i] * ΔPmax[i]
            g[i] = Wobs[i] * ΔP[i]
            gλ[i] = -ΔPmax[i] * g[i]
        else
            gλ[i] = zero(T)
        end
    end
    return term
end

function update_term!(term::L2TargetTerm, prob::NonlinearProblem)
    @unpack g, x, αₓ = term
    @turbo for i = 1:length(g)
        g[i] = αₓ[i] * x[i]
    end
end

function grad!(g, targ::TargetFunction{T}, n) where {T}
    fill!(g, zero(T))
    terms = @inbounds values(targ.terms)[1:end-1]
    # FIXED: Using 'map' instead of 'for' preserves 'type-stability'
    map(term -> (g .+= @inbounds term.gviews[n]), terms)
    return g
end

getvalue(targ::TargetFunction) = sum(getvalue, values(targ.terms))
getvalues(targ::TargetFunction) = map(getvalue, targ.terms)

function getvalue(term::AbstractTargetTerm{T}) where {T}
    @unpack ΔP, g = term
    val = zero(T)
    @turbo for i = 1:length(ΔP)
        val += ΔP[i] * g[i]
    end
    return T(0.5) * val
end

function getvalue(term::L2TargetTerm{T}) where {T}
    @unpack x, g = term
    val = zero(T)
    @turbo for i = 1:length(x)
        val += x[i] * g[i]
    end
    return T(0.5) * val
end

function calc_well_index!(targ::TargetFunction)
    calc_well_index!(targ.prob, targ.wviews.Gw)
    calc_well_index!(targ.prob, targ.wviews.Jinj)
end

function calc_well_index!(prob::NonlinearProblem, wviews::WellIndexViews{:Gw})
    @unpack xviews, pviews = wviews
    @unpack Pcalc, Pbhp, Jp, M, Gw, Qliq_h = prob.params
    @inbounds @simd for i = 1:length(xviews)
        fill!(pviews[i], inv(xviews[i]))
    end
    @turbo for i = 1:length(Pcalc)
        Jp[i] = Gw[i] * M[i]
        Pbhp[i] = Pcalc[i] - Qliq_h[i] / Jp[i]
    end
end

function calc_well_index!(prob::NonlinearProblem, wviews::WellIndexViews{:Jinj})
    @unpack xviews, pviews, λviews = wviews
    @unpack Pcalc, Pinj, Jinj, Qinj_h, λ = prob.params
    @inbounds @simd for i = 1:length(xviews)
        fill!(pviews[i], λviews[i] * inv(xviews[i]))
    end
    @turbo for i = 1:length(Pcalc)
        Pinj[i] = Pcalc[i] + λ[i] * Qinj_h[i] / Jinj[i]
    end
end


text:
# add all packages here
using DataFrames
using Gadfly
using Vega # for word cloud
using TextAnalysis #, DimensionalityReduction, Clustering
using Plotly
using Query

#include("clean-data.jl")

    function fetchTopNTopic(df, min_count = 100 )
        # define min count
        # filter out
        df = df |> @query(i, begin
                    @where i.count>min_count
                    @orderby descending(i.count)
                    @select {i.term, i.count}
                  end) |> DataFrame

         plot = @time Gadfly.plot(df,y="count", x="term", Geom.bar, color="term")
        # save plot
        #draw(PNG("plot/diagnosis_top_n.png", 3inch, 3inch), plot)
        # return
         return(df,plot)
     end

     function fetchNgramTopic(full_text_array, min_count = 100, ngram =3 )
         # define min count
         # filter out
         full_text_sd, full_text_tx= cleanText( join(full_text_array,"\n"))

        # create ngrams
        full_text_ng=ngrams(full_text_sd, ngram)
        # convert to dataframe
        full_text_ng_df = DataFrame(term=collect(keys(full_text_ng)), count=collect(values(full_text_ng)))

        # add new column size
        full_text_ng_df = full_text_ng_df |> @query(i, begin
                            @orderby descending(i.count)
                            @select {i.term, i.count, size=length(split(get(i.term)))}
                          end) |> DataFrame

        # filter out < min_count
        full_text_ng_df = full_text_ng_df |> @query(i, begin
                            @where i.count>min_count && i.size==ngram
                            @select {i.term, i.count,i.size }
                          end) |> DataFrame

        # plot graph
        plot = @time Gadfly.plot(full_text_ng_df,y="count", x="term",
                Geom.bar(orientation= :vertical), color="term",
                #Guide.xlabel(nothing),
                #Guide.ylabel(nothing),
                Theme(key_position = :none))
        # save plot
         #draw(PNG("plot/diagnosis_top_n.png", 3inch, 3inch), plot)
         # return
          return(full_text_ng_df,plot)
      end

    # this function generates trends for an array of genes, cells, concepts, therapy
  function generateTrends(df_full_text, terms=[], startYear = 2008, endYear =2018 )

                  # Create df
          trend_df = DataFrame(
                          term=String[],
                          count=Int64[],
                          year=Int64[])

            for yr in startYear:endYear
                # filter by year
                filtered_df = df_full_text |> @query(i, begin
                                        @where i.year==yr
                                        @select {i.fullText}
                                      end) |> DataFrame

                # clean text and convert to sd
                sd_array=[]
                for row in eachrow(filtered_df)
                    sd,tx=cleanText(row[:fullText])
                    push!(sd_array,sd)
                end

                # convert to corpus
                crp = Corpus(sd_array)
                #standardize
                standardize!(crp, StringDocument)
                #update crp
                update_lexicon!(crp)
                update_inverse_index!(crp)
                lexicon(crp)

                for term in terms
                    # query corpus
                   count=sum(crp[term])
                    # create df and push it to trend_df
                     temp_df = DataFrame(term=term, count=count, year=yr)
                     append!(trend_df, temp_df)

                end

            end

           return(trend_df)
  end

ext: jl

max_stars_repo_name: jmulloka/Breast-Cancer-Text-Mining

max_issues_repo_name: jmulloka/Breast-Cancer-Text-Mining

max_forks_repo_name: jmulloka/Breast-Cancer-Text-Mining


text:
using Test
using Spandex

@testset "cholesky_sym 1" begin
    local g = DirectedGraph{Int64}(3)
    g[1, 1] = 6
    g[2, 1] = 8
    g[2, 2] = 27
    g[3, 1] = 14
    g[3, 2] = 27
    g[3, 3] = 41
    local ata = from_graph(g, 3, 3)
    ata.layout = Spandex.lower_symmetric

    local b = cholesky_sym(ata)

    @test 6 == b.nnz

    @test contains(b, 1, 1)
    @test contains(b, 2, 1)
    @test contains(b, 2, 2)
    @test contains(b, 3, 1)
    @test contains(b, 3, 2)
    @test contains(b, 3, 3)
end

@testset "cholesky_sym 2" begin
    local g = DirectedGraph{Int64}(11)
    g[1, 1] = 1
    g[2, 2] = 2
    g[3, 2] = 3
    g[3, 3] = 3
    g[4, 4] = 4
    g[5, 5] = 5
    g[6, 1] = 6
    g[6, 4] = 6
    g[6, 6] = 6
    g[7, 1] = 7
    g[7, 7] = 7
    g[8, 2] = 8
    g[8, 5] = 8
    g[8, 8] = 8
    g[9, 6] = 9
    g[9, 9] = 9
    g[10, 3] = 10
    g[10, 4] = 10
    g[10, 6] = 10
    g[10, 8] = 10
    g[10, 10] = 10
    g[11, 3] = 11
    g[11, 5] = 11
    g[11, 7] = 11
    g[11, 8] = 11
    g[11, 10] = 11
    g[11, 11] = 11
    local ata = from_graph(g, 11, 11)
    ata.layout = Spandex.lower_symmetric

    local b = cholesky_sym(ata)

    @test 33 == b.nnz

    @test contains(b, 7, 6)
    @test contains(b, 8, 3)
    @test contains(b, 9, 7)
    @test contains(b, 10, 7)
    @test contains(b, 10, 9)
    @test contains(b, 11, 9)
end

@testset "cholesky_to! 1" begin
    local g = DirectedGraph{Float64}(3)
    g[1, 1] = 6.0
    g[2, 1] = 8.0
    g[2, 2] = 27.0
    g[3, 1] = 14.0
    g[3, 2] = 27.0
    g[3, 3] = 41.0
    local ata = from_graph(g, 3, 3)
    ata.layout = Spandex.lower_symmetric

    clear!(g)
    g[1, 1] = 6.0
    g[2, 1] = 1.0 + 1.0 / 3.0
    g[2, 2] = 16.0 + 1.0 / 3.0
    g[3, 1] = 2.0 + 1.0 / 3.0
    g[3, 2] = 0.510204081632653
    g[3, 3] = 4.08163265306122
    local e = from_graph(g, 3, 3)
    e.layout = Spandex.lower_symmetric

    local ld = cholesky_sym(ata)
    cholesky_to!(ata, ld)

    ld.values = round.(ld.values, digits = 8)
    e.values = round.(e.values, digits = 8)

    @test equals(ld, e)
end

@testset "cholesky_to! 2" begin
    local g = DirectedGraph{Float64}(5)
    g[1, 1] = 0.454154210872255
    g[2, 1] = 0.493313382040145
    g[2, 2] = 0.673117517240105
    g[3, 1] = 0.0
    g[3, 2] = 0.267770063806461
    g[3, 3] = 0.753202118586013
    g[4, 1] = 0.597779567627811
    g[4, 2] = 0.650049443659237
    g[4, 3] = 0.337426886372401
    g[4, 4] = 1.72955430012221
    g[5, 1] = 0.0939734363675742
    g[5, 2] = 0.0200830025404996
    g[5, 3] = 0.0
    g[5, 4] = 0.25140727798654
    g[5, 5] = 0.527264476708613
    local ata = from_graph(g, 5, 5)
    ata.layout = Spandex.lower_symmetric

    clear!(g)
    g[1, 1] = 0.454154210872255
    g[2, 1] = 1.08622439301549
    g[2, 2] = 0.137268488267129
    g[3, 1] = 0.0
    g[3, 2] = 1.95070308697049
    g[3, 3] = 0.230862228520464
    g[4, 1] = 1.31624799091856
    g[4, 2] = 0.00529397289085283
    g[4, 3] = 1.45545384824111
    g[4, 4] = 0.453678241855123
    g[5, 1] = 0.206919663228681
    g[5, 2] = -0.597320167013766
    g[5, 3] = 0.692813459608376
    g[5, 4] = -0.2306541683174
    g[5, 5] = 0.323895344304452
    local e = from_graph(g, 5, 5)
    e.layout = Spandex.lower_symmetric

    local ld = cholesky_sym(ata)
    cholesky_to!(ata, ld)

    ld.values = round.(ld.values, digits = 8)
    e.values = round.(e.values, digits = 8)

    @test equals(ld, e)
end

@testset "solve_to! 1" begin
    local g = DirectedGraph{Float64}(3)
    g[1, 1] = 38.0
    g[2, 1] = 22.0
    g[2, 2] = 17.0
    g[3, 1] = 14.0
    g[3, 2] = 13.0
    g[3, 3] = 11.0
    local ata = from_graph(g, 3, 3)
    ata.layout = Spandex.lower_symmetric

    local b = [1.0, 10.0, 5.0]

    local ld = cholesky_sym(ata)
    cholesky_to!(ata, ld)
    local x = zeros(Float64, 3)
    solve_to!(ld, b, x)

    @test 3 == length(x)

    @test abs(-9.5 - x[1]) < 1e-10
    @test abs(34.166666666666 - x[2]) < 1e-10
    @test abs(-27.833333333333 - x[3]) < 1e-10
end

@testset "solve_to! 2" begin
    local g = DirectedGraph{Float64}(5)
    g[1, 1] = 0.339721892599889
    g[2, 1] = -0.117511897064494
    g[2, 2] = 0.61348986691309
    g[3, 1] = 0.0221940137659884
    g[3, 2] = 0.27504299861335
    g[3, 3] = 0.355516636030647
    g[4, 1] = -0.26882517335515
    g[4, 2] = 0.0684899239330446
    g[4, 3] = 0.0333434870019433
    g[4, 4] = 0.412101013254036
    g[5, 1] = -0.159648639043336
    g[5, 2] = 0.0401190234458681
    g[5, 3] = -0.0318276646882203
    g[5, 4] = -0.0108046250911296
    g[5, 5] = 0.399329937547505
    local ata = from_graph(g, 5, 5)
    ata.layout = Spandex.lower_symmetric

    local b = [4.0, -1.0, -3.0, 4.0, -2.0]

    local ld = cholesky_sym(ata)
    cholesky_to!(ata, ld)
    local x = zeros(Float64, 5)
    solve_to!(ld, b, x)

    @test 5 == length(x)

    @test abs(84.3926371078906 - x[1]) < 1e-10
    @test abs(20.6069097491701 - x[2]) < 1e-10
    @test abs(-33.4133473479173 - x[3]) < 1e-10
    @test abs(64.7118533875248 - x[4]) < 1e-10
    @test abs(25.7485304673898 - x[5]) < 1e-10
end

@testset "solve 1" begin
    local g = DirectedGraph{Float64}(10)
    g[1, 1] = 0.360464443870286
    g[3, 2] = 0.965038079655014
    g[10, 2] = 0.806541221607173
    g[1, 3] = 0.156202523064209
    g[3, 3] = 0.70277194218269
    g[7, 3] = 0.398688926587124
    g[9, 3] = 0.158532504726658
    g[10, 3] = 0.070915819808533
    g[7, 4] = 0.552895404215196
    g[4, 5] = 0.97656582830328
    g[6, 5] = 0.362469500523493
    g[2, 6] = 0.510437505153131
    g[3, 6] = 0.473695871041683
    g[5, 7] = 0.477123911915246
    g[8, 7] = 0.582754540178946
    g[4, 8] = 0.828533162691592
    g[5, 9] = 0.612247361949774
    g[6, 10] = 0.570109021624869
    local a = from_graph(g, 10, 10)

    local solver = CholeskySolver{Float64}(10, 10)
    solver.use_permutation = false
    solver.use_normalization = false

    solve_sym(solver, a)

    local b = collect(1.0:10.0)

    local x = solve(solver, a, b)

    @test abs(4.6288540824669431 - x[1]) < 1e-10
    @test abs(9.0912268747716602 - x[2]) < 1e-10
    @test abs(-4.2799392703685015 - x[3]) < 1e-10
    @test abs(15.746856145276686 - x[4]) < 1e-10
    @test abs(-0.60521269977679459 - x[5]) < 1e-10
    @test abs(-0.59628461522440024 - x[6]) < 1e-10
    @test abs(13.727906774511695 - x[7]) < 1e-10
    @test abs(5.5411542327921595 - x[8]) < 1e-10
    @test abs(-2.5315136968936010 - x[9]) < 1e-10
    @test abs(10.909090909090908 - x[10]) < 1e-10

    solver.use_permutation = true
    solve_sym(solver, a)

    x = solve(solver, a, b)

    @test abs(4.6288540824669431 - x[1]) < 1e-10
    @test abs(9.0912268747716602 - x[2]) < 1e-10
    @test abs(-4.2799392703685015 - x[3]) < 1e-10
    @test abs(15.746856145276686 - x[4]) < 1e-10
    @test abs(x[5]) < 1e-10
    @test abs(-0.59628461522440024 - x[6]) < 1e-10
    @test abs(13.727906774511695 - x[7]) < 1e-10
    @test abs(4.8278091694066978 - x[8]) < 1e-10
    @test abs(-2.5315136968936010 - x[9]) < 1e-10
    @test abs(10.524302848075246 - x[10]) < 1e-10
end

@testset "solve 2" begin
    local g = DirectedGraph{Float64}(3)
    g[1, 2] = 1.0
    g[1, 3] = 1.0
    g[2, 1] = 2.0
    g[2, 2] = 4.0
    g[2, 3] = -2.0
    g[3, 2] = 3.0
    g[3, 3] = 15.0
    local a = from_graph(g, 3, 3)

    local solver = CholeskySolver{Float64}(3, 3)
    solver.use_permutation = false
    solver.use_normalization = false

    solve_sym(solver, a)

    local b = [17.0, 2.89, -3.3]

    local x = solve(solver, a, b)

    @test abs(-46.1300 - x[1]) < 1e-8
    @test abs(21.5250 - x[2]) < 1e-8
    @test abs(-4.5250 - x[3]) < 1e-8

    local y = [
        x[2] * get_rowwise(a, 1, 2) + x[3] * get_rowwise(a, 1, 3),
        x[1] * get_rowwise(a, 2, 1) +
        x[2] * get_rowwise(a, 2, 2) +
        x[3] * get_rowwise(a, 2, 3),
        x[2] * get_rowwise(a, 3, 2) + x[3] * get_rowwise(a, 3, 3),
    ]

    @test abs(b[1] - y[1]) < 1e-8
    @test abs(b[2] - y[2]) < 1e-8
    @test abs(b[3] - y[3]) < 1e-8

    solver.use_permutation = true
    solve_sym(solver, a)
    x = solve(solver, a, b)

    @test abs(-46.1300 - x[1]) < 1e-8
    @test abs(21.5250 - x[2]) < 1e-8
    @test abs(-4.5250 - x[3]) < 1e-8
end

@testset "solve 3" begin
    local g = DirectedGraph{Float64}(10)
    g[1, 1] = 0.360464443870286
    g[3, 2] = 0.965038079655014
    g[10, 2] = 0.806541221607173
    g[1, 3] = 0.156202523064209
    g[3, 3] = 0.70277194218269
    g[7, 3] = 0.398688926587124
    g[9, 3] = 0.158532504726658
    g[10, 3] = 0.070915819808533
    g[7, 4] = 0.552895404215196
    g[4, 5] = 0.97656582830328
    g[6, 5] = 0.362469500523493
    g[2, 6] = 0.510437505153131
    g[3, 6] = 0.473695871041683
    g[5, 7] = 0.477123911915246
    g[8, 7] = 0.582754540178946
    g[4, 8] = 0.828533162691592
    g[5, 9] = 0.612247361949774
    g[6, 10] = 0.570109021624869
    local a = from_graph(g, 10, 10)

    local solver = CholeskySolver{Float64}(10, 10)
    solver.use_permutation = false
    solver.use_normalization = false

    solve_sym(solver, a)

    local b = collect(1.0:10.0)

    local x = solve(solver, a, b)

    local y = mul(a, x)

    solver.use_normalization = true

    x = solve(solver, a, b)
    local z = mul(a, x)

    @test sum((y .- z) .^ 2) < 1e-8
end

@testset "solve 4" begin
    local g = DirectedGraph{Float64}(10)
    g[1, 1] = 0.360464443870286
    g[3, 2] = 0.965038079655014
    g[10, 2] = 0.806541221607173
    g[1, 3] = 0.156202523064209
    g[3, 3] = 0.70277194218269
    g[7, 3] = 0.398688926587124
    g[9, 3] = 0.158532504726658
    g[10, 3] = 0.070915819808533
    g[7, 4] = 0.552895404215196
    g[4, 5] = 0.97656582830328
    g[6, 5] = 0.362469500523493
    g[2, 6] = 0.510437505153131
    g[3, 6] = 0.473695871041683
    g[5, 7] = 0.477123911915246
    g[8, 7] = 0.582754540178946
    g[4, 8] = 0.828533162691592
    g[5, 9] = 0.612247361949774
    g[6, 10] = 0.570109021624869
    local a = from_graph(g, 10, 10)

    local solver = CholeskySolver{Float64}(10, 10)
    solver.use_permutation = true
    solver.use_normalization = false

    solve_sym(solver, a)

    local b = collect(1.0:10.0)

    local x = solve(solver, a, b)

    local y = mul(a, x)

    solver.use_normalization = true

    x = solve(solver, a, b)
    local z = mul(a, x)

    @test sum((y .- z) .^ 2) < 1e-8
end

@testset "update 1" begin
    local g = DirectedGraph{Float64}(3)
    g[1, 2] = 1.0
    g[1, 3] = 1.0
    g[2, 1] = 2.0
    g[2, 2] = 4.0
    g[2, 3] = -2.0
    g[3, 2] = 3.0
    g[3, 3] = 15.0
    local a = from_graph(g, 3, 3)

    local solver = CholeskySolver{Float64}(3, 3)
    solver.use_permutation = false
    solver.use_normalization = false

    solve_sym(solver, a)

    local b = [17.0, 2.89, -3.3]

    local x = solve(solver, a, b)

    local m = SparseArray{Float64}(3)
    m[1] = 7.0
    m[2] = -5.0
    m[3] = 1.0
    local u = update!(solver, m, 9.0)

    add_vertex!(g)
    g[4, 1] = m[1]
    g[4, 2] = m[2]
    g[4, 3] = m[3]
    a = from_graph(g, 4, 3)
    push!(b, 9.0)
    x = solve(solver, a, b)

    @test sum((x .- u) .^ 2) < 1e-8
end

@testset "update 2" begin
    local g = DirectedGraph{Float64}(3)
    g[1, 2] = 1.0
    g[1, 3] = 1.0
    g[2, 1] = 2.0
    g[2, 2] = 4.0
    g[2, 3] = -2.0
    g[3, 2] = 3.0
    g[3, 3] = 15.0
    local a = from_graph(g, 3, 3)

    local solver = CholeskySolver{Float64}(3, 3)
    solver.use_permutation = true
    solver.use_normalization = false

    solve_sym(solver, a)

    local b = [17.0, 2.89, -3.3]

    local x = solve(solver, a, b)

    local m = SparseArray{Float64}(3)
    m[1] = 7.0
    m[2] = -5.0
    m[3] = 1.0
    local u = update!(solver, m, 9.0)

    add_vertex!(g)
    g[4, 1] = m[1]
    g[4, 2] = m[2]
    g[4, 3] = m[3]
    a = from_graph(g, 4, 3)
    push!(b, 9.0)
    x = solve(solver, a, b)

    @test sum((x .- u) .^ 2) < 1e-8
end

@testset "update 3" begin
    local g = DirectedGraph{Float64}(3)
    g[1, 2] = 1.0
    g[1, 3] = 1.0
    g[2, 1] = 2.0
    g[2, 2] = 4.0
    g[2, 3] = -2.0
    g[3, 2] = 3.0
    g[3, 3] = 15.0
    local a = from_graph(g, 3, 3)

    local solver = CholeskySolver{Float64}(3, 3)
    solver.use_permutation = false
    solver.use_normalization = false

    solve_sym(solver, a)

    local b = [17.0, 2.89, -3.3]

    local x = solve(solver, a, b)

    local m = SparseArray{Float64}(3)
    m[3] = 1.0
    local u = update!(solver, m, 9.0)

    add_vertex!(g)
    g[4, 3] = 1.0
    a = from_graph(g, 4, 3)
    push!(b, 9.0)
    x = solve(solver, a, b)

    @test sum((x .- u) .^ 2) < 1e-8
end

@testset "update 4" begin
    local g = DirectedGraph{Float64}(10)
    g[1, 1] = 0.360464443870286
    g[3, 2] = 0.965038079655014
    g[10, 2] = 0.806541221607173
    g[1, 3] = 0.156202523064209
    g[3, 3] = 0.70277194218269
    g[7, 3] = 0.398688926587124
    g[9, 3] = 0.158532504726658
    g[10, 3] = 0.070915819808533
    g[7, 4] = 0.552895404215196
    g[4, 5] = 0.97656582830328
    g[6, 5] = 0.362469500523493
    g[2, 6] = 0.510437505153131
    g[3, 6] = 0.473695871041683
    g[5, 7] = 0.477123911915246
    g[8, 7] = 0.582754540178946
    g[4, 8] = 0.828533162691592
    g[5, 9] = 0.612247361949774
    g[6, 10] = 0.570109021624869
    local a = from_graph(g, 10, 10)

    local solver = CholeskySolver{Float64}(10, 10)
    solver.use_permutation = true
    solver.use_normalization = false

    solve_sym(solver, a)

    local b = collect(1.0:10.0)

    local x = solve(solver, a, b)

    local m = SparseArray{Float64}(10)
    m[2] = 0.5
    m[3] = 0.1
    m[6] = 0.9
    local u = update!(solver, m, 11.0)

    add_vertex!(g)
    g[11, 2] = m[2]
    g[11, 3] = m[3]
    g[11, 6] = m[6]
    a = from_graph(g, 11, 10)
    push!(b, 11.0)
    x = solve(solver, a, b)

    @test sum((x .- u) .^ 2) < 1e-8
end

@testset "downdate 1" begin
    local g = DirectedGraph{Float64}(3)
    g[1, 2] = 1.0
    g[1, 3] = 1.0
    g[2, 1] = 2.0
    g[2, 2] = 4.0
    g[2, 3] = -2.0
    g[3, 2] = 3.0
    g[3, 3] = 15.0
    local a = from_graph(g, 3, 3)

    local solver = CholeskySolver{Float64}(3, 3)
    solver.use_permutation = false
    solver.use_normalization = false

    solve_sym(solver, a)

    local b = [17.0, 2.89, -3.3]

    local x = solve(solver, a, b)

    local m = SparseArray{Float64}(3)
    m[1] = 3.0
    m[2] = 2.0
    m[3] = 1.0
    local u = update!(solver, m, 9.0)

    local d = downdate!(solver, m, 9.0)

    @test sum((x .- d) .^ 2) < 1e-8
end

@testset "downdate 2" begin
    local g = DirectedGraph{Float64}(10)
    g[1, 1] = 0.360464443870286
    g[3, 2] = 0.965038079655014
    g[10, 2] = 0.806541221607173
    g[1, 3] = 0.156202523064209
    g[3, 3] = 0.70277194218269
    g[7, 3] = 0.398688926587124
    g[9, 3] = 0.158532504726658
    g[10, 3] = 0.070915819808533
    g[7, 4] = 0.552895404215196
    g[4, 5] = 0.97656582830328
    g[6, 5] = 0.362469500523493
    g[2, 6] = 0.510437505153131
    g[3, 6] = 0.473695871041683
    g[5, 7] = 0.477123911915246
    g[8, 7] = 0.582754540178946
    g[4, 8] = 0.828533162691592
    g[5, 9] = 0.612247361949774
    g[6, 10] = 0.570109021624869
    local a = from_graph(g, 10, 10)

    local solver = CholeskySolver{Float64}(10, 10)
    solver.use_permutation = true
    solver.use_normalization = false

    solve_sym(solver, a)

    local b = collect(1.0:10.0)

    local x = solve(solver, a, b)

    local m = SparseArray{Float64}(10)
    m[2] = 0.5
    m[3] = 0.1
    m[6] = 0.9
    local u = update!(solver, m, 11.0)

    @test sum((x .- u) .^ 2) > 0

    u = downdate!(solver, m, 11.0)

    @test sum((x .- u) .^ 2) < 1e-8
end

@testset "downdate 3" begin
    local g = DirectedGraph{Float64}(10)
    g[1, 1] = 0.360464443870286
    g[3, 2] = 0.965038079655014
    g[10, 2] = 0.806541221607173
    g[1, 3] = 0.156202523064209
    g[3, 3] = 0.70277194218269
    g[7, 3] = 0.398688926587124
    g[9, 3] = 0.158532504726658
    g[10, 3] = 0.070915819808533
    g[7, 4] = 0.552895404215196
    g[4, 5] = 0.97656582830328
    g[6, 5] = 0.362469500523493
    g[2, 6] = 0.510437505153131
    g[3, 6] = 0.473695871041683
    g[5, 7] = 0.477123911915246
    g[8, 7] = 0.582754540178946
    g[4, 8] = 0.828533162691592
    g[5, 9] = 0.612247361949774
    g[6, 10] = 0.570109021624869
    local a = from_graph(g, 10, 10)

    local solver = CholeskySolver{Float64}(10, 10)
    solver.use_permutation = false
    solver.use_normalization = true

    solve_sym(solver, a)

    local b = collect(1.0:10.0)

    local x = solve(solver, a, b)

    local m = SparseArray{Float64}(10)
    m[2] = 0.5
    m[3] = 0.1
    m[6] = 0.9
    local u = update!(solver, m, 11.0)
    u = downdate!(solver, m, 11.0)

    @test sum((x .- u) .^ 2) < 1e-8
end

@testset "downdate 4" begin
    local g = DirectedGraph{Float64}(10)
    g[1, 1] = 0.360464443870286
    g[3, 2] = 0.965038079655014
    g[10, 2] = 0.806541221607173
    g[1, 3] = 0.156202523064209
    g[3, 3] = 0.70277194218269
    g[7, 3] = 0.398688926587124
    g[9, 3] = 0.158532504726658
    g[10, 3] = 0.070915819808533
    g[7, 4] = 0.552895404215196
    g[4, 5] = 0.97656582830328
    g[6, 5] = 0.362469500523493
    g[2, 6] = 0.510437505153131
    g[3, 6] = 0.473695871041683
    g[5, 7] = 0.477123911915246
    g[8, 7] = 0.582754540178946
    g[4, 8] = 0.828533162691592
    g[5, 9] = 0.612247361949774
    g[6, 10] = 0.570109021624869
    local a = from_graph(g, 10, 10)

    local solver = CholeskySolver{Float64}(10, 10)
    solver.use_permutation = true
    solver.use_normalization = true

    solve_sym(solver, a)

    local b = collect(1.0:10.0)

    local x = solve(solver, a, b)

    local m = SparseArray{Float64}(10)
    m[2] = 0.5
    m[3] = 0.1
    m[6] = 0.9
    local u = update!(solver, m, 11.0)
    u = downdate!(solver, m, 11.0)

    @test sum((x .- u) .^ 2) < 1e-8
end

@testset "example 1" begin
    local n = 3
    local a = from_csr(
        n,
        n,
        [1, 3, 6, 8],
        [2, 3, 1, 2, 3, 2, 3],
        [1.0, 1.0, 2.0, 4.0, -2.0, 3.0, 15.0],
    )

    local solver = CholeskySolver{Float64}(n, n)
    # # by default
    # solver.use_permutation = true
    # solver.use_normalization = true

    solve_sym(solver, a)

    local x = solve(solver, a, [17.0, 2.89, -3.3])

    local m = SparseArray{Float64}(n)
    m[1] = 7.0
    m[2] = -5.0
    m[3] = 1.0

    local u = update!(solver, m, 9.0)
    local v = downdate!(solver, m, 9.0)

    @assert sum(abs.(x - u)) > 0
    @assert sum(abs.(x - v)) < 1e-8
end


text:
################################################################################
#
#   DEFINITION OF LT HAMILTONIAN (STRUCT ONLY FOR COMPUTATION)
#   - label type L (of bond labels)
#   - unitcell type U<:AbstractUnitcell{S,B} where {S, B<:AbstractBond{L,N} where {N}},
#   - bond hamiltonian type HB<:AbstractBondHamiltonian{L,NS} where {NS}
#
################################################################################
struct LTHamiltonian{
        L,
        U<:AbstractUnitcell{S,B} where {S, B<:AbstractBond{L,N} where {N}},
        HB<:AbstractBondHamiltonian{L,NS} where {NS}
    }

    # the unitcell
    unitcell :: U

    # the bond Hamiltonian
    h_bond   :: HB

    # ...

end

# function to create such an object
function newLTHamiltonian(
            unitcell :: U,
            bond_hamiltonian :: HB
        ) :: LTHamiltonian{L,U,HB} where {L,U<:AbstractUnitcell{S,B} where {S, B<:AbstractBond{L,N} where {N}},HB<:AbstractBondHamiltonian{L,NS} where {NS}}

    # return such an object
    return LTHamiltonian{L,U,HB}(unitcell, bond_hamiltonian)
end

# delta function for unitcell and bond to compute the displacement vector of the bond
function delta(
        uc :: U,
        b  :: B
    ) :: Vector{Float64} where {L,N,S,B<:AbstractBond{L,N},U<:AbstractUnitcell{S,B}}

    # build up the difference vector within the unitcell
    d = point(site(uc,to(b))) .- point(site(uc,from(b)))
    # add all the bravais lattice vectors
    for i in 1:N
        d .+= wrap(b)[i] .* latticeVectors(uc)[i]
    end
    # return the vector
    return d
end

# function to obtain the Hamiltonian matrix at some k
function getMatrixAtK(
        hamiltonian :: LTHamiltonian{L,U,HB},
        k :: Vector{<:Real}
    ) :: Matrix{Complex} where {L,NS,U,HB<:AbstractBondHamiltonian{L,NS}}

    # new matrix
    h = zeros(Complex, NS*numSites(hamiltonian.unitcell),NS*numSites(hamiltonian.unitcell))

    # add all bonds to the matrix
    for b in bonds(hamiltonian.unitcell)
        # get the bond indices
        i_from  = from(b)
        i_to    = to(b)
        delta_r = delta(hamiltonian.unitcell, b)
        # get the interaction matrix and add it to the general matrix
        h[(i_from-1)*(NS)+1:(i_from)*(NS), (i_to-1)*(NS)+1:(i_to)*(NS)] .+= 0.5 * bondterm(hamiltonian.h_bond, b) .* exp(im*dot(k,delta_r))
    end

    # return the matrix
    return h
end

# function determine the energy of a k_vector for a given LT hamiltonian
function energy(
            hamiltonian :: LTHamiltonian{L,U,HB},
            k           :: Vector{Float64}
        ) :: Float64 where {L,NS,U,HB<:AbstractBondHamiltonian{L,NS}}

    # first compute the matrix at the given k
    M_k = getMatrixAtK(hamiltonian, k)

    # return the minimal eigenvalue
    return eigmin(M_k)
end


text:
__precompile__()

module Nabla

    using SpecialFunctions
    using LinearAlgebra
    using Statistics

    # Some aliases used repeatedly throughout the package.
    export ∇Scalar, ∇Array, SymOrExpr, ∇ArrayOrScalar
    const ∇Scalar = Number
    const ∇Array = AbstractArray{<:∇Scalar}
    const ∇AbstractVector = AbstractVector{<:∇Scalar}
    const ∇AbstractMatrix = AbstractMatrix{<:∇Scalar}
    const ∇ArrayOrScalar = Union{AbstractArray{<:∇Scalar}, ∇Scalar}
    const SymOrExpr = Union{Symbol, Expr}

    # ones/zeros(::AbstractArray) is deprecated in 0.7 and removed in 1.0, but it's a
    # pretty useful method, so we'll define our own for internal use
    for f in (:ones, :zeros)
        like = Symbol(f, "like")
        @eval begin
            $(like)(a::AbstractArray) = $(f)(eltype(a), size(a))
            $(like)(n::Integer) = $(f)(n)
        end
    end

    # Meta-programming utilities specific to Nabla.
    include("code_transformation/util.jl")
    include("code_transformation/differentiable.jl")

    # Functionality for constructing computational graphs.
    include("core.jl")

    # Functionality for defining new sensitivities.
    include("sensitivity.jl")

    # Finite differencing functionality - only used in tests. Would be good to move this
    # into a separate module at some point.
    include("finite_differencing.jl")

    # Sensitivities for the basics.
    include("sensitivities/indexing.jl")
    include("sensitivities/scalar.jl")
    include("sensitivities/array.jl")

    # Sensitivities for functionals.
    include("sensitivities/functional/functional.jl")
    include("sensitivities/functional/reduce.jl")
    include("sensitivities/functional/reducedim.jl")

    # Linear algebra optimisations.
    include("sensitivities/linalg/generic.jl")
    include("sensitivities/linalg/symmetric.jl")
    include("sensitivities/linalg/strided.jl")
    include("sensitivities/linalg/blas.jl")
    include("sensitivities/linalg/diagonal.jl")
    include("sensitivities/linalg/triangular.jl")
    include("sensitivities/linalg/factorization/cholesky.jl")
    include("sensitivities/linalg/factorization/svd.jl")

    # Checkpointing
    include("checkpointing.jl")

end # module Nabla


text:
##file to include in ComplexModeler.jl which contains all files associated to Matpower

include("MatpowerGenerator.jl")
include("MatpowerLoad.jl")
include("MatpowerShunt.jl")
include("MatpowerLine_π.jl")
include("MatpowerVolt.jl")
include("read_matpower.jl")
include("read_matpower_simple.jl")

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
using Calmato
using Test

@testset "Calmato.jl" begin
    db1 = search_db("Cu Zn Liang Hsiao 2015")
    print(db1)
    sys = init_system(db1)
    print(sys)
    res = solidify(sys)
    print(res)

    db2 = search_db("Cu H O S P Magnusson")
    db2 = select(db2, "Cu S", *, -1)
    print(db2)

    db3 = merge(db1, db2, 1 => 1, 2 => 4, *)
    print(db3)
    sys = init_system(db3)
    print(sys)
    X = [0.7, 0.3, 0.01]
    Ts = 600:100:1300
    res = solidify(sys, X, Ts)
    print(res)
end


text:
function varianceSD(x)

    n = length(x)

    (1/(n-1)) * sum((x .- meanSD(x)).^2)

end


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

@doc doc"""
    AdaDelta(; kwargs...)

Scale learning rates by the ratio of accumulated gradients to accumulated
updates, see [1] and notes for further description.

### Attributes
* `η`: default `1.0`, learning rate.
* `ρ`: default `0.95`, squared gradient moving average decay factor.
* `ϵ`: default `1e-6`, small value added for numerical stability.
* `clip`: default `0`, gradient clipping.
  If positive, will clip the gradient into the range `[-clip, clip]`.
* `scale`: default `0`, gradient rescaling.
  If != 0, multiply the gradient with `scale` before updating.
  Often choose to be `1.0 / batch_size`.
  If leave it default, high-level API like `fit!` will set it to
  `1.0 / batch_size`, since `fit!` knows the `batch_size`.
* `λ`: default `0.00001`, weight decay is equivalent
  to adding a global l2 regularizer for all the parameters.

### Notes
`ρ` should be between 0 and 1. A value of `ρ` close to 1 will decay the
moving average slowly and a value close to 0 will decay the moving average
fast.

`ρ = 0.95` and `ϵ = 1e-6` are suggested in the paper and reported to
work for multiple datasets (MNIST, speech). In the paper, no learning rate is
considered (so `η = 1.0`). Probably best to keep it at this value.

`ϵ` is important for the very first update (so the numerator does not become 0).

Using the step size `η` and a decay factor `ρ` the learning rate is
calculated as:

```math
\begin{align*}
  r_t &= ρ r_{t-1} + (1 - ρ) g^2 \\
  η_t &= η \frac{\sqrt{s_{t-1} + ϵ}} {\sqrt{r_t + ϵ}} \\
  s_t &= ρ s_{t-1} + (1 - ρ) (η_t \times g)^2
\end{align*}
```

### References
1. Zeiler, M. D. (2012):
   ADADELTA: An Adaptive Learning Rate Method. arXiv Preprint arXiv:1212.5701.
"""
AdaDelta

@defstruct AdaDelta <: AbstractOptimizer (
  (η      :: Real = 1.0,  η > 0),
  (ρ      :: Real = 0.95, 0 < ρ < 1 ),
  (ϵ      :: Real = 1e-6, ϵ > 0),
  (clip   :: Real = 0,    clip >= 0),
   scale  :: Real = 0,
  (λ      :: Real = 1e-5, λ >= 0),
  η_sched :: Any  = initlrsched(η)
)

mutable struct AdaDeltaState
  x  :: NDArray
  Δx :: NDArray
end

create_state(::AdaDelta, ::Int, W::NDArray) =
  AdaDeltaState(zeros(size(W), context(W)), zeros(size(W), context(W)))

function update!(ada::AdaDelta, ::Int, W::NDArray, ∇::NDArray, s::AdaDeltaState)
  η  = get(ada.η_sched)
  x  = s.x
  Δx = s.Δx
  ρ  = ada.ρ
  ϵ  = ada.ϵ

  normgrad!(ada, W, ∇)

  # Update s.acc as in RMSProp
  @inplace x .*= ρ
  @inplace x .+= (1 - ρ) .* ∇.^2

  # Compute update using the "old" Δx
  Δxₜ = ∇ .* sqrt(Δx .+ ϵ) ./ sqrt(x .+ ϵ)  # FIXME: sqrt dot-call
  @inplace W .+= -η .* Δxₜ

  # update Δx using update
  @inplace Δx .*= ρ
  @inplace Δx .+= (1 - ρ) .* Δxₜ.^2
end


text:
#=

Model definition for *SSE models

Ignacio Quintero Mächler

t(-_-t)

27 05 2019

=#

"""
    define_mod(esse_mod::String, x::Array{Float64,1}, y::Array{Float64,1}, k::Int64)

Defines ESSE model.
"""
function define_mod(esse_mod::String,
                    k       ::Int64,
                    af      ::Function,
                    md      ::Bool)

  if occursin(r"^[s|S][A-za-z]*", esse_mod)         # if speciation
    mod_ode = md ? make_esse_s(k, af, md) : make_esse_s(k, af)
    npars   = 2k + k*k
    pardic  = build_par_names(k, true)
    ws      = true
    printstyled("running speciation ESSE model \n", color=:green)

  elseif occursin(r"^[e|E][A-za-z]*", esse_mod)     # if extinction
    mod_ode = md ? make_esse_e(k, af, md) : make_esse_e(k, af)
    npars   = (2k + k*k)::Int64
    pardic  = build_par_names(k, true)
    ws      = false

    printstyled("running extinction ESSE model \n", color=:green)

  elseif occursin(r"^[t|T|r|R|q|Q][A-za-z]*", esse_mod) # if transition
    mod_ode = md ? make_esse_q(k, af, md) : make_esse_q(k, af)
    npars   = 2k*k::Int64
    pardic  = build_par_names(k, false)
    ws      = false

    printstyled("running transition ESSE model \n", color=:green)

  else
    error("esse_mod does not match any of the alternatives:
          speciation, extinction or transition")
  end

  return mod_ode, npars, pardic, md, ws

end

"""
    define_mod(cov_mod::String,
               k            ::Int64,
               h            ::Int64,
               ny           ::Int64)

Defines EGeoHiSSE model for `k` areas, `h` hidden states and `ny` covariates.
"""
function define_mod(cov_mod::NTuple{N,String},
                    k      ::Int64,
                    h      ::Int64,
                    ny     ::Int64) where {N}

  model = [false, false, false]

  for m in cov_mod
    # if speciation
    if occursin(r"^[s|S][A-za-z]*", m)
      model[1] = true
    end
    # if extinction
    if occursin(r"^[e|E][A-za-z]*", m)
      model[2] = true
    end
    # if dispersal
    if occursin(r"^[t|T|r|R|q|Q][A-za-z]*", m)
      model[3] = true
    end
  end

  mexp = "$(model[1] ? "speciation," : "")$(model[2] ? "extinction," : "")$(model[3] ? "transition," : "")"
  mexp = replace(mexp, "," => ", ")
  mexp = mexp[1:(end-2)]

  printstyled("running $mexp EGeoHiSSE model with:
    $k single areas
    $h hidden states
    $ny covariates \n", color=:green)

  return tuple(model...)
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
import Base: ndims

abstract type TestFunction{N,M} end
ndims(::TestFunction{N,M}) where{N,M} = (N,M)
constrains(::TestFunction) = ()
include("single_obj_functions.jl")
include("multi_obj_functions.jl")


text:
module GreyDecision

include("./greynumber.jl")
include("./utility.jl")
include("./stats.jl")
include("./random.jl")
include("./topsis.jl")
include("./vikor.jl")
include("./saw.jl")
include("./nds.jl")
include("./moora.jl")
include("./electre.jl")
include("./wpm.jl")
include("./marcos.jl")
include("./aras.jl")

import ..GreyNumbers: GreyNumber, kernel
import ..Topsis: topsis
import ..Vikor: vikor
import ..Saw: saw
import ..NDS: nds
import ..Moora: moora
import ..Electre: electre
import ..Wpm: wpm
import ..Marcos: marcos
import ..Aras: aras

export GreyNumbers
export Utility
export Random
export Stats

export GreyNumber
export topsis, vikor, saw, nds, moora, electre, wpm, marcos, aras

end # module GreyDecision


text:
# Autogenerated wrapper script for GEOS_jll for aarch64-apple-darwin
export libgeos, libgeos_cpp

JLLWrappers.@generate_wrapper_header("GEOS")
JLLWrappers.@declare_library_product(libgeos, "@rpath/libgeos_c.1.dylib")
JLLWrappers.@declare_library_product(libgeos_cpp, "@rpath/libgeos.3.10.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libgeos,
        "lib/libgeos_c.1.16.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgeos_cpp,
        "lib/libgeos.3.10.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
include("BUDE.jl")
using ROCKernels, CUDAKernels, KernelAbstractions, StaticArrays, CUDA, AMDGPU

@enum Backend cuda rocm cpu

function list_rocm_devices()::Vector{DeviceWithRepr}
  try
    # AMDGPU.agents()'s internal iteration order isn't stable
    sorted = sort(AMDGPU.get_agents(:gpu), by = repr)
    map(x -> (x, repr(x), rocm), sorted)
  catch
    # probably unsupported
    []
  end
end

function list_cuda_devices()::Vector{DeviceWithRepr}
  return !CUDA.functional(false) ? [] :
         map(d -> (d, "$(CUDA.name(d)) ($(repr(d)))", cuda), CUDA.devices())
end

function devices()::Vector{DeviceWithRepr}
  cudas = list_cuda_devices()
  rocms = list_rocm_devices()
  cpus = [(undef, "$(Sys.cpu_info()[1].model) ($(Threads.nthreads())T)", cpu)]
  vcat(cpus, cudas, rocms)
end

function run(params::Params, deck::Deck, device::DeviceWithRepr)

  nposes::Int64 = size(deck.poses)[2]
  blocks_size::Int64 = ceil(Int64, nposes / params.ppwi)
  nthreads::Int64 = params.wgsize

  (selected, _, backend) = device
  println("KernelAbstractions Backend: $backend")

  if backend == cpu
    println("Using CPU with max $(Threads.nthreads()) threads")
    protein = deck.protein
    ligand = deck.ligand
    forcefield = deck.forcefield
    poses = deck.poses
    etotals = Array{Float32}(undef, nposes)
    backend_impl = CPU()
  elseif backend == cuda
    CUDA.device!(selected)
    if CUDA.device() != selected
      error("Cannot select CUDA device, expecting $selected, but got $(CUDA.device())")
    end
    if !CUDA.functional(true)
      error("Non-functional CUDA configuration")
    end
    println("Using CUDA device: $(CUDA.name(selected)) ($(repr(selected)))")

    protein = CuArray{Atom}(deck.protein)
    ligand = CuArray{Atom}(deck.ligand)
    forcefield = CuArray{FFParams}(deck.forcefield)
    poses = CuArray{Float32,2}(deck.poses)
    etotals = CuArray{Float32}(undef, nposes)
    backend_impl = CUDADevice()
  elseif backend == rocm
    AMDGPU.DEFAULT_AGENT[] = selected
    if AMDGPU.get_default_agent() != selected
      error("Cannot select HSA device, expecting $selected, but got $(AMDGPU.get_default_agent())")
    end
    println("Using GPU HSA device: $(AMDGPU.get_name(selected)) ($(repr(selected)))")

    protein = ROCArray{Atom}(deck.protein)
    ligand = ROCArray{Atom}(deck.ligand)
    forcefield = ROCArray{FFParams}(deck.forcefield)
    poses = ROCArray{Float32,2}(deck.poses)
    etotals = ROCArray{Float32}(undef, nposes)
    backend_impl = ROCDevice()
  else
    error("unsupported backend $(backend)")
  end

  println("Using kernel parameters: <<<$blocks_size,$nthreads>>> 1:$nposes")

  kernel! = fasten_main(backend_impl, nthreads, blocks_size)

  wait(
    kernel!(
      nthreads,
      nposes,
      protein,
      ligand,
      forcefield,
      poses,
      etotals,
      Val(convert(Int, params.ppwi)),
      ndrange = blocks_size,
    ),
  )

  elapsed = @elapsed for i = 1:params.iterations
    wait(
      kernel!(
        nthreads,
        nposes,
        protein,
        ligand,
        forcefield,
        poses,
        etotals,
        Val(convert(Int, params.ppwi)),
        ndrange = blocks_size,
      ),
    )
  end

  (Array{Float32}(etotals), elapsed, params.ppwi)
end

@fastmath @kernel function fasten_main(
  wgsize::Int64,
  nposes::Int64,
  @Const(protein),
  @Const(ligand),
  @Const(forcefield),
  @Const(poses),
  etotals,
  ::Val{PPWI},
) where {PPWI}

  gid = @index(Group)
  lid = @index(Local)

  #  Get index of first TD
  ix = (gid - 1) * wgsize * PPWI + lid

  #  Have extra threads do the last member intead of return.
  #  A return would disable use of barriers, so not using return is better
  ix = ix <= nposes ? ix : nposes - PPWI

  etot = MArray{Tuple{PPWI},Float32}(undef)
  transform = MArray{Tuple{PPWI,3},Vec4f32}(undef)

  lsz = wgsize
  for i = 1:PPWI
    index = (ix) + (i - 1) * lsz
    @inbounds sx::Float32 = sin(poses[1, index])
    @inbounds cx::Float32 = cos(poses[1, index])
    @inbounds sy::Float32 = sin(poses[2, index])
    @inbounds cy::Float32 = cos(poses[2, index])
    @inbounds sz::Float32 = sin(poses[3, index])
    @inbounds cz::Float32 = cos(poses[3, index])
    @inbounds transform[i, 1] = #
      Vec4f32(cy * cz, sx * sy * cz - cx * sz, cx * sy * cz + sx * sz, poses[4, index])
    @inbounds transform[i, 2] = #
      Vec4f32(cy * sz, sx * sy * sz + cx * cz, cx * sy * sz - sx * cz, poses[5, index])
    @inbounds transform[i, 3] = #
      Vec4f32(-sy, sx * cy, cx * cy, poses[6, index])
    etot[i] = 0
  end

  @inbounds for l_atom::Atom in ligand

    l_params::FFParams = forcefield[l_atom.type+1]
    lhphb_ltz::Bool = l_params.hphb < Zero
    lhphb_gtz::Bool = l_params.hphb > Zero

    lpos = MArray{Tuple{PPWI},Vec3f32}(undef)

    @simd for i = 1:PPWI
      # Transform ligand atom
      @inbounds lpos[i] = Vec3f32(
        transform[i, 1].w +
        l_atom.x * transform[i, 1].x +
        l_atom.y * transform[i, 1].y +
        l_atom.z * transform[i, 1].z,
        transform[i, 2].w +
        l_atom.x * transform[i, 2].x +
        l_atom.y * transform[i, 2].y +
        l_atom.z * transform[i, 2].z,
        transform[i, 3].w +
        l_atom.x * transform[i, 3].x +
        l_atom.y * transform[i, 3].y +
        l_atom.z * transform[i, 3].z,
      )
    end

    @inbounds for p_atom::Atom in protein
      @inbounds p_params::FFParams = forcefield[p_atom.type+1]

      radij::Float32 = p_params.radius + l_params.radius
      r_radij::Float32 = One / radij

      elcdst::Float32 = (p_params.hbtype == HbtypeF && l_params.hbtype == HbtypeF) ? Four : Two
      elcdst1::Float32 = (p_params.hbtype == HbtypeF && l_params.hbtype == HbtypeF) ? Quarter : Half
      type_E::Bool = ((p_params.hbtype == HbtypeE || l_params.hbtype == HbtypeE))

      phphb_ltz::Bool = p_params.hphb < Zero
      phphb_gtz::Bool = p_params.hphb > Zero
      phphb_nz::Bool = p_params.hphb != Zero
      p_hphb::Float32 = p_params.hphb * (phphb_ltz && lhphb_gtz ? -One : One)
      l_hphb::Float32 = l_params.hphb * (phphb_gtz && lhphb_ltz ? -One : One)
      distdslv::Float32 =
        (phphb_ltz ? (lhphb_ltz ? Npnpdist : Nppdist) : (lhphb_ltz ? Nppdist : -Float32Max))
      r_distdslv::Float32 = One / distdslv

      chrg_init::Float32 = l_params.elsc * p_params.elsc
      dslv_init::Float32 = p_hphb + l_hphb

      @simd for i = 1:PPWI
        @inbounds x::Float32 = lpos[i].x - p_atom.x
        @inbounds y::Float32 = lpos[i].y - p_atom.y
        @inbounds z::Float32 = lpos[i].z - p_atom.z

        distij::Float32 = sqrt(x ^ 2 + y ^ 2 + z ^ 2)

        # Calculate the sum of the sphere radii
        distbb::Float32 = distij - radij
        zone1::Bool = (distbb < Zero)

        # Calculate steric energy
        @inbounds etot[i] += (One - (distij * r_radij)) * (zone1 ? Two * Hardness : Zero)

        # Calculate formal and dipole charge interactions
        chrg_e::Float32 =
          chrg_init * ((zone1 ? One : (One - distbb * elcdst1)) * (distbb < elcdst ? One : Zero))
        neg_chrg_e::Float32 = -abs(chrg_e)
        chrg_e = type_E ? neg_chrg_e : chrg_e
        @inbounds etot[i] += chrg_e * Cnstnt

        # Calculate the two cases for Nonpolar-Polar repulsive interactions
        coeff::Float32 = (One - (distbb * r_distdslv))
        dslv_e::Float32 = dslv_init * ((distbb < distdslv && phphb_nz) ? One : Zero)
        dslv_e *= (zone1 ? One : coeff)
        @inbounds etot[i] += dslv_e
      end
    end
  end

  td_base = (gid - 1) * wgsize * PPWI + lid
  if td_base <= nposes
    @simd for i = 1:PPWI
      @inbounds etotals[td_base+(i-1)*wgsize] = etot[i] * Half
    end
  end

  nothing
end

main()

max_stars_repo_name: wsmoses/Enzyme-BUDE

max_issues_repo_name: wsmoses/Enzyme-BUDE

max_forks_repo_name: UoB-HPC/miniBUDE


text:
using LinearAlgebra
using Distributions: InverseWishart
using Random
using HighFrequencyCovariance

twister = MersenneTwister(10)
IW_dist(n) = InverseWishart(n, Matrix(Float64.(I(n))))

# slightly non psd non diagonal version.
A = Hermitian([1.02 0.5 0.9; 0.5 0.98 0.9; 0.9 0.9 1.01])
# Diagonal weighting.
W_root = Diagonal(rand(twister, size(A)[1]))
invW = inv(W_root^2)
U_ed = project_to_U(A, invW)
all( abs.(diag(U_ed) .- 1.0) .< 100*eps())
isa(U_ed, Hermitian)
S_ed = project_to_S(A, W_root)
all(eigen(S_ed).values .> -100*eps())

# hermitian weighting.
W_root = Hermitian(rand(twister, IW_dist(3)  ))
invW = inv(W_root^2)
U_ed = project_to_U(A, invW)
all( abs.(diag(U_ed) .- 1.0) .< 100*eps())
isa(U_ed, Hermitian)
S_ed = project_to_S(A, W_root)
all(eigen(S_ed).values .> -100*eps())

# Really non psd non diagonal version.
A = Hermitian([1.02 0.5 0.9; 0.5 0.98 -0.9; 0.9 -0.9 1.01])
# Diagonal weighting.
W_root = Diagonal(rand(twister, size(A)[1]))
invW = inv(W_root^2)
U_ed = project_to_U(A, invW)
all( abs.(diag(U_ed) .- 1.0) .< 100*eps())
isa(U_ed, Hermitian)
S_ed = project_to_S(A, W_root)
all(eigen(S_ed).values .> -100*eps())

# hermitian weighting.
W_root = Hermitian(rand(twister, IW_dist(3)  ))
invW = inv(W_root^2)
U_ed = project_to_U(A, invW)
all( abs.(diag(U_ed) .- 1.0) .< 100*eps())
isa(U_ed, Hermitian)
S_ed = project_to_S(A, W_root)
all(eigen(S_ed).values .> -100000*eps())

# Diagonal A
A = Diagonal([1.02 0 0; 0 0.98 0; 0 0 1.05])
# Diagonal weighting.
W_root = Diagonal(rand(twister, size(A)[1]))
invW = inv(W_root^2)
U_ed = project_to_U(A, invW)
all( abs.(diag(U_ed) .- 1.0) .< 100*eps())
isa(U_ed, Diagonal)
S_ed = project_to_S(A, W_root)
all(eigen(S_ed).values .> -100*eps())

# hermitian weighting.
W_root = Hermitian(rand(twister, IW_dist(3)  ))
invW = inv(W_root^2)
U_ed = project_to_U(A, invW)
all( abs.(diag(U_ed) .- 1.0) .< 100*eps())
isa(U_ed, Hermitian)
S_ed = project_to_S(A, W_root)
all(eigen(S_ed).values .> -100000*eps())

# Copypasting this here to avoid exporting it
function sqrt_psd(A::Hermitian)
    eigenvalues, eigenvectors = eigen(A)
    return Hermitian(construct_matrix_from_eigen(sqrt.(eigenvalues), eigenvectors))
end
sqrt_psd(A::Diagonal) = sqrt(A)

# A identify matrix case.
N = 4
A = Hermitian(I(N))
W_root = Hermitian(rand(twister, IW_dist(N)  ))
Dykstra = Hermitian(zeros(N,N))
W_inv = inv(W_root^2)
W_inv_sqrt = sqrt_psd(W_inv)
Y = A

Y, Dykstra = iterate_higham(Y, Dykstra, W_root, W_inv, W_inv_sqrt)
all(abs.(Y .- I(N)) .< 100*eps())
Y, Dykstra = iterate_higham(Y, Dykstra, W_root, W_inv, W_inv_sqrt)
all(abs.(Y .- I(N)) .< 100*eps())

# If A were input as diagonal
A = Diagonal(I(N))
Dykstra = Diagonal(zeros(N,N))
Y = A
Y, Dykstra = iterate_higham(Y, Dykstra, W_root, W_inv, W_inv_sqrt)
all(abs.(Y .- I(N)) .< 100*eps())
Y, Dykstra = iterate_higham(Y, Dykstra, W_root, W_inv, W_inv_sqrt)
all(abs.(Y .- I(N)) .< 100*eps())

# If A and E_root were diagonal
A = Diagonal(Float64.(I(N)))
W_root = Diagonal(rand(twister, N))
Dykstra = Diagonal(zeros(N,N))
W_inv = inv(W_root^2)
W_inv_sqrt = sqrt_psd(W_inv)
Y = A
Y, Dykstra = iterate_higham(Y, Dykstra, W_root, W_inv, W_inv_sqrt)
all(abs.(Y .- I(N)) .< 100*eps())
Y, Dykstra = iterate_higham(Y, Dykstra, W_root, W_inv, W_inv_sqrt)
all(abs.(Y .- I(N)) .< 100*eps())

# A valid correlation matrix case
N = 3
A, _ = cov2cor(Hermitian(rand(twister, IW_dist(N)  )))
W_root = Hermitian(rand(twister, IW_dist(N)  ))
Dykstra = Hermitian(zeros(N,N))
W_inv = inv(W_root^2)
W_inv_sqrt = sqrt_psd(W_inv)
Y = A

Y, Dykstra = iterate_higham(Y, Dykstra, W_root, W_inv, W_inv_sqrt)
all(abs.(Y .- A) .< 100*eps())
Y, Dykstra = iterate_higham(Y, Dykstra, W_root, W_inv, W_inv_sqrt)
all(abs.(Y .- A) .< 100*eps())

# A Slightly non psd case
N = 3
A = Hermitian([1.02 0.5 0.9; 0.5 0.98 0.9; 0.9 0.9 1.01])
W = Hermitian(abs.(cov2cor(Hermitian(rand(twister, IW_dist(N)  )))[1]))
updated_matrix, counter, convergence = nearest_correlation_matrix(A, W; doDykstra = true, stop_at_first_correlation_matrix = true, max_iterates = 1000)
updated_matrix2, counter2, convergence2 = nearest_correlation_matrix(A, W; doDykstra = true, stop_at_first_correlation_matrix = false, max_iterates = 10000)
updated_matrix3, counter3, convergence3 = nearest_correlation_matrix(A, W; doDykstra = false, stop_at_first_correlation_matrix = true, max_iterates = 10000)
all(abs.(updated_matrix .- updated_matrix2) .< 100*eps()) # These should be the same as once it is psd the matrix stops updating.
any(abs.(updated_matrix .- updated_matrix3) .> 0.0001) # These should be different as without Dykstra one is biased
# Testing convergences.
minimum(eigen(updated_matrix).values) > -10*eps()
minimum(eigen(updated_matrix2).values) > -10*eps()
minimum(eigen(updated_matrix3).values) > -10*eps()

# With diagonal weighting.
N = 3
A = Hermitian([1.02 0.5 0.9; 0.5 0.98 0.9; 0.9 0.9 1.01])
W = Diagonal(rand(twister,N))
updated_matrix, counter, convergence = nearest_correlation_matrix(A, W; doDykstra = true, stop_at_first_correlation_matrix = true, max_iterates = 1000)
updated_matrix2, counter2, convergence2 = nearest_correlation_matrix(A, W; doDykstra = true, stop_at_first_correlation_matrix = false, max_iterates = 10000)
updated_matrix3, counter3, convergence3 = nearest_correlation_matrix(A, W; doDykstra = false, stop_at_first_correlation_matrix = true, max_iterates = 10000)
all(abs.(updated_matrix .- updated_matrix2) .< 100*eps()) # These should be the same as once it is psd the matrix stops updating.
any(abs.(updated_matrix .- updated_matrix3) .> 0.0001) # These should be different as without Dykstra one is biased
# Testing convergences.
minimum(eigen(updated_matrix).values) > -10*eps()
minimum(eigen(updated_matrix2).values) > -10*eps()
minimum(eigen(updated_matrix3).values) > -10*eps()

# A really non psd case
A = Hermitian([1.02 0.5 0.9; 0.5 0.98 -0.9; 0.9 -0.9 1.01])
W = Hermitian(abs.(cov2cor(Hermitian(rand(twister, IW_dist(N)  )))[1]))
updated_matrix, counter, convergence = nearest_correlation_matrix(A, W; doDykstra = true, stop_at_first_correlation_matrix = true, max_iterates = 100000)
updated_matrix2, counter2, convergence2 = nearest_correlation_matrix(A, W; doDykstra = true, stop_at_first_correlation_matrix = false, max_iterates = 10000)
updated_matrix3, counter3, convergence3 = nearest_correlation_matrix(A, W; doDykstra = false, stop_at_first_correlation_matrix = true, max_iterates = 10000)
all(abs.(updated_matrix .- updated_matrix2) .< 1000*eps()) # These should be the same as once it is psd the matrix stops updating.
any(abs.(updated_matrix .- updated_matrix3) .> 0.0001) # These should be different as without Dykstra one is biased
# Testing convergences.
minimum(eigen(updated_matrix).values) > -10*eps()
minimum(eigen(updated_matrix2).values) > -10*eps()


text:
#   Copyright 2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
Unit tests for the methods in `jltest/cli.jl`.

Notes
-----
* For the unit tests in this files, failures and errors are expected.
"""

# --- Imports

# Standard library
using Test

# External packages
using Suppressor

# Local modules
using TestTools.jltest: cli, EnhancedTestSet

# --- Tests

@testset EnhancedTestSet "jltest.cli.parse_args()" begin

    # --- Default arguments

    raw_args = Vector{String}()
    args = cli.parse_args(; raw_args=raw_args)
    expected_args = Dict(
        "fail-fast" => false,
        "no-wrapper" => false,
        "no-recursion" => false,
        "verbose" => false,
        "version" => false,
        "tests" => Vector{String}(),
    )
    @test args == expected_args

    # --- fail-fast

    # Case: raw_args = "--fail-fast"
    raw_args = ["--fail-fast"]
    args = cli.parse_args(; raw_args=raw_args)
    expected_args = Dict(
        "fail-fast" => true,
        "no-wrapper" => false,
        "no-recursion" => false,
        "verbose" => false,
        "version" => false,
        "tests" => Vector{String}(),
    )
    @test args == expected_args

    # Case: raw_args = "-x"
    raw_args = ["-x"]
    args = cli.parse_args(; raw_args=raw_args)
    expected_args = Dict(
        "fail-fast" => true,
        "no-wrapper" => false,
        "no-recursion" => false,
        "verbose" => false,
        "version" => false,
        "tests" => Vector{String}(),
    )
    @test args == expected_args

    # --- no-wrapper

    # Case: raw_args = "--no-wrapper"
    raw_args = ["--no-wrapper"]
    args = cli.parse_args(; raw_args=raw_args)
    expected_args = Dict(
        "fail-fast" => false,
        "no-wrapper" => true,
        "no-recursion" => false,
        "verbose" => false,
        "version" => false,
        "tests" => Vector{String}(),
    )
    @test args == expected_args

    # Case: raw_args = "-W"
    raw_args = ["-W"]
    args = cli.parse_args(; raw_args=raw_args)
    expected_args = Dict(
        "fail-fast" => false,
        "no-wrapper" => true,
        "no-recursion" => false,
        "verbose" => false,
        "version" => false,
        "tests" => Vector{String}(),
    )
    @test args == expected_args

    # --- no-recursion

    # Case: raw_args = "--no-recursion"
    raw_args = ["--no-recursion"]
    args = cli.parse_args(; raw_args=raw_args)
    expected_args = Dict(
        "fail-fast" => false,
        "no-wrapper" => false,
        "no-recursion" => true,
        "verbose" => false,
        "version" => false,
        "tests" => Vector{String}(),
    )
    @test args == expected_args

    # Case: raw_args = "-R"
    raw_args = ["-R"]
    args = cli.parse_args(; raw_args=raw_args)
    expected_args = Dict(
        "fail-fast" => false,
        "no-wrapper" => false,
        "no-recursion" => true,
        "verbose" => false,
        "version" => false,
        "tests" => Vector{String}(),
    )
    @test args == expected_args

    # --- verbose

    # Case: raw_args = "--verbose"
    raw_args = ["--verbose"]
    args = cli.parse_args(; raw_args=raw_args)
    expected_args = Dict(
        "fail-fast" => false,
        "no-wrapper" => false,
        "no-recursion" => false,
        "verbose" => true,
        "version" => false,
        "tests" => Vector{String}(),
    )
    @test args == expected_args

    # Case: raw_args = "-v"
    raw_args = ["-v"]
    args = cli.parse_args(; raw_args=raw_args)
    expected_args = Dict(
        "fail-fast" => false,
        "no-wrapper" => false,
        "no-recursion" => false,
        "verbose" => true,
        "version" => false,
        "tests" => Vector{String}(),
    )
    @test args == expected_args

    # --- version

    # Case: raw_args = "--version"
    raw_args = ["--version"]
    args = cli.parse_args(; raw_args=raw_args)
    expected_args = Dict(
        "fail-fast" => false,
        "no-wrapper" => false,
        "no-recursion" => false,
        "verbose" => false,
        "version" => true,
        "tests" => Vector{String}(),
    )
    @test args == expected_args

    # Case: raw_args = "-V"
    raw_args = ["-V"]
    args = cli.parse_args(; raw_args=raw_args)
    expected_args = Dict(
        "fail-fast" => false,
        "no-wrapper" => false,
        "no-recursion" => false,
        "verbose" => false,
        "version" => true,
        "tests" => Vector{String}(),
    )
    @test args == expected_args

    # --- tests

    raw_args = ["test-1", "test-2.jl"]
    args = cli.parse_args(; raw_args=raw_args)
    expected_args = Dict(
        "fail-fast" => false,
        "no-wrapper" => false,
        "no-recursion" => false,
        "verbose" => false,
        "version" => false,
        "tests" => ["test-1", "test-2.jl"],
    )
    @test args == expected_args
end

@testset EnhancedTestSet "jltest.cli.run(): basic tests" begin
    # --- Preparations

    # Construct path to test directory
    test_dir = joinpath(@__DIR__, "data-basic-tests")
    test_dir_relpath = relpath(test_dir)

    # Precompute commonly used values
    some_tests_file = joinpath(test_dir, "some_tests.jl")
    expected_output_some_tests = "$(joinpath(test_dir_relpath, "some_tests")): .."

    some_tests_no_testset_file = joinpath(test_dir, "some_tests_no_testset.jl")
    expected_output_some_tests_no_testset = "$(joinpath(test_dir_relpath, "some_tests_no_testset")): .."

    failing_tests_file = joinpath(test_dir, "failing_tests.jl")
    expected_output_failing_tests = strip(
        """
        $(joinpath(test_dir_relpath, "failing_tests")): .
        =====================================================
        failing tests: Test Failed at $(failing_tests_file):27
          Expression: 2 == 1
           Evaluated: 2 == 1

        Stacktrace:
         [1]
        """
    )

    expected_output_failing_tests_fail_fast = strip(
        """
        $(joinpath(test_dir_relpath, "failing_tests")): .
        =====================================================
        Test Failed at $(failing_tests_file):27
          Expression: 2 == 1
           Evaluated: 2 == 1

        =====================================================
        Error During Test at
        """
    )

    failing_tests_no_testset_file = joinpath(test_dir, "failing_tests_no_testset.jl")
    expected_output_failing_tests_no_testset = strip(
        """
        $(joinpath(test_dir_relpath, "failing_tests_no_testset")): .
        =====================================================
        All tests: Test Failed at $(failing_tests_no_testset_file):26
          Expression: 2 == 1
           Evaluated: 2 == 1

        Stacktrace:
          [1]
        """
    )

    more_tests_file = joinpath(test_dir, "subdir", "more_tests.jl")
    expected_output_more_tests = "$(joinpath(test_dir_relpath, "subdir", "more_tests")): .."

    # --- Tests

    # Case: normal operation
    tests = [some_tests_file, some_tests_no_testset_file]
    output = strip(@capture_out begin
        cli.run(tests)
    end)

    @test occursin(expected_output_some_tests, output)
    @test occursin(expected_output_some_tests_no_testset, output)

    # Case: fail_fast = true
    tests = [failing_tests_file, some_tests_no_testset_file]
    local error = nothing
    output = strip(@capture_out begin
        try
            cli.run(tests; fail_fast=true)
        catch error
        end
    end)

    @test error isa Test.FallbackTestSetException
    @test error.msg == "There was an error during testing"

    @test startswith(output, expected_output_failing_tests_fail_fast)
    @test !occursin("some_tests_no_testset", output)

    # Case: use_wrapper = false, fail_fast = true
    tests = [failing_tests_file, some_tests_no_testset_file]
    local error = nothing
    output = strip(@capture_out begin
        try
            cli.run(tests; fail_fast=true, use_wrapper=false)
        catch error
        end
    end)

    @test error isa Test.FallbackTestSetException
    @test error.msg == "There was an error during testing"

    @test startswith(output, expected_output_failing_tests_fail_fast)
    @test !occursin("some_tests_no_testset", output)

    # Case: use_wrapper = false, fail_fast = false
    tests = [failing_tests_file, some_tests_no_testset_file]
    local error = nothing
    output = strip(@capture_out begin
        try
            cli.run(tests; use_wrapper=false)
        catch error
        end
    end)

    @test isnothing(error)

    @test startswith(output, expected_output_failing_tests)
    @test occursin(expected_output_some_tests_no_testset, output)

    # Case: recursive = false
    tests = [test_dir]
    local error = nothing
    log_msg = strip(@capture_err begin
        output = strip(@capture_out begin
            try
                cli.run(tests; recursive=false)
            catch error
            end
        end)
    end)

    @test isnothing(error)

    @test startswith(output, expected_output_failing_tests)
    @test occursin(expected_output_failing_tests_no_testset, output)
    @test occursin(expected_output_some_tests, output)
    @test occursin(expected_output_some_tests_no_testset, output)
    @test !occursin("more_tests", output)

    # Case: `tests` is empty
    cd(test_dir)
    tests = []
    local error = nothing
    log_msg = strip(@capture_err begin
        output = strip(@capture_out begin
            try
                cli.run(tests)
            catch error
            end
        end)
    end)

    @test isnothing(error)

    expected_output_failing_tests = strip(
        """
        failing_tests: .
        =====================================================
        failing tests: Test Failed at $(failing_tests_file):27
          Expression: 2 == 1
           Evaluated: 2 == 1

        Stacktrace:
         [1]
        """
    )
    @test startswith(output, expected_output_failing_tests)

    expected_output_failing_tests_no_testset = strip(
        """
        failing_tests_no_testset: .
        =====================================================
        All tests: Test Failed at $(failing_tests_no_testset_file):26
          Expression: 2 == 1
           Evaluated: 2 == 1

        Stacktrace:
          [1]
        """
    )
    @test occursin(expected_output_failing_tests_no_testset, output)

    expected_output_some_tests = "some_tests: .."
    @test occursin(expected_output_some_tests, output)

    expected_output_some_tests_no_testset = "some_tests_no_testset: .."
    @test occursin(expected_output_some_tests_no_testset, output)

    expected_output_more_tests = "$(joinpath("subdir", "more_tests")): .."
    @test occursin(expected_output_more_tests, output)
end

@testset EnhancedTestSet "jltest.cli.run(): error cases" begin

    # --- Exercise functionality and check results

    # Case: invalid `tests` arg
    @test_throws MethodError cli.run([1, 2, 3])
end

# --- Emit message about expected failures and errors

println()
@info "For $(basename(@__FILE__)), 5 failures and 0 errors are expected."


text:
function getref(repo::GitRepo, ref::AbstractString)
    refpath = joinpath(repo.path, "refs", ref)
    if isfile(refpath)
        # look in refs subdirectory
        return SHA1(readline(refpath))
    else
        packedref = joinpath(repo.path, "packed-refs")
        if isfile(packedref)
            for str in eachline(packedref)
                # '#' are comments
                # '^' are peeled refs of the preceding lines
                if !startswith(str, ('#','^')) && ref == str[47:end]
                    return SHA1(str[1:40])
                end
            end
        end
    end
    error("Could not find $ref")
end

max_stars_repo_name: simonbyrne/GitX

max_issues_repo_name: simonbyrne/GitX

max_forks_repo_name: simonbyrne/GitX


text:
using QuasiCopula, LinearAlgebra, Random, GLM
using DataFrames, Statistics
using BenchmarkTools, Test

BLAS.set_num_threads(1)
Threads.nthreads()

p = 3   # number of fixed effects, including intercept
m = 1    # number of variance components
# true parameter values
Random.seed!(12345)
βtrue = rand(Uniform(-2, 2), p)
θtrue = [0.5]
τtrue = 100.0
σ2 = inv(τtrue)
σ = sqrt(σ2)

#simulation parameters
samplesize = 10000
ni = 25

# true Gamma
Γ = θtrue[1] * ones(ni, ni)

T = Float64
gcs = Vector{GaussianCopulaVCObs{T}}(undef, samplesize)

# for reproducibility I will simulate all the design matrices here
Random.seed!(12345)
X_samplesize = [randn(ni, p - 1) for i in 1:samplesize]

for i in 1:samplesize
    X = [ones(ni) X_samplesize[i]]
    μ = X * βtrue
    vecd = Vector{ContinuousUnivariateDistribution}(undef, ni)
    for i in 1:ni
        vecd[i] = Normal(μ[i], σ)
    end
    nonmixed_multivariate_dist = NonMixedMultivariateDistribution(vecd, Γ)
    # simuate single vector y
    y = Vector{Float64}(undef, ni)
    res = Vector{Float64}(undef, ni)
    rand(nonmixed_multivariate_dist, y, res)
    V = [ones(ni, ni)]
    gcs[i] = GaussianCopulaVCObs(y, X, V)
end

# form VarLmmModel
gcm = GaussianCopulaVCModel(gcs)
# precompile
println("precompiling Gaussian VCM fit")
gcm2 = deepcopy(gcm);
QuasiCopula.fit!(gcm2, IpoptSolver(print_level = 0, max_iter = 20));

fittime = @elapsed QuasiCopula.fit!(gcm, IpoptSolver(print_level = 5, max_iter = 100, limited_memory_max_history = 50, accept_after_max_steps = 4, tol = 10^-8, hessian_approximation = "limited-memory"))
@show fittime
@show gcm.β
@show gcm.θ
@show gcm.τ
@show gcm.∇β
@show gcm.∇θ
@show gcm.∇τ

loglikelihood!(gcm, true, true)
vcov!(gcm)
@show QuasiCopula.confint(gcm)

mseβ, mseτ, mseθ = MSE(gcm, βtrue, τtrue[1], θtrue)
@show mseβ
@show mseτ
@show mseθ

@test mseβ < 0.01
@test mseθ < 0.01
@test mseτ < 0.01


text:
using DomTerm
using Test

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


text:
@testset "curve_length" begin
    f(x) = sin(x)

    sol = 3.8201977890277120179047620821714432919099676146472747210804966567
    @test isapprox(curve_length(f,0,pi), sol; atol = 1e-10)
    @test !isapprox(curve_length(f,0,5*pi;rtol=1e-4),5*sol,atol=1e-7)
end

@testset "curve_fraction" begin
    f(x) = sin(x)

    @test isapprox(curve_fraction(f,0,pi,0.5), pi/2; atol = 1e-10)
    @test !isapprox(curve_fraction(f,0,pi,0.5;rtol=1e-2), pi/2; atol = 1e-9)
end


text:
function ob(r,Side)
    x = 0
    Flag = true
    while isborder(r,HorizonSide(Side)) && !isborder(r,HorizonSide((Side + 3) % 4))
        move!(r,HorizonSide((Side + 3) % 4))
        x += 1
    end
    if isborder(r,HorizonSide(Side)) && isborder(r,HorizonSide((Side + 3) % 4))
        Flag = false
    end
    if  !isborder(r,HorizonSide(Side))
        move!(r,HorizonSide(Side))
    end
    while !isborder(r,HorizonSide(Side)) && isborder(r,HorizonSide((Side + 1) % 4))
        move!(r,HorizonSide(Side))
    end
    for i = 1:x
        move!(r,HorizonSide((Side+1) % 4))
    end
    return Flag
end

function ob2(r,Side,y)
    x = 0
    Flag = true
    while isborder(r,HorizonSide(Side)) && !isborder(r,HorizonSide((Side + 3) % 4))
        move!(r,HorizonSide((Side + 3) % 4))
        x += 1
    end
    if isborder(r,HorizonSide(Side)) && isborder(r,HorizonSide((Side + 3) % 4))
        Flag = false
    end
    if  !isborder(r,HorizonSide(Side))
        move!(r,HorizonSide(Side))
        y+=1
    end
    while !isborder(r,HorizonSide(Side)) && isborder(r,HorizonSide((Side + 1) % 4))
        move!(r,HorizonSide(Side))
        y+= 1
    end
    for i = 1:x
        move!(r,HorizonSide((Side+1) % 4))
    end
    return Flag,y
end

function ob3(r,Side,y)
    x = 0
    Flag = true
    while isborder(r,HorizonSide(Side)) && !isborder(r,HorizonSide((Side + 3) % 4))
        move!(r,HorizonSide((Side + 3) % 4))
        x += 1
    end
    if isborder(r,HorizonSide(Side)) && isborder(r,HorizonSide((Side + 3) % 4))
        Flag = false
    end
    if  !isborder(r,HorizonSide(Side))
        move!(r,HorizonSide(Side))
        y-=1
    end
    while !isborder(r,HorizonSide(Side)) && isborder(r,HorizonSide((Side + 1) % 4))
        move!(r,HorizonSide(Side))
        y-= 1
    end
    for i = 1:x
        move!(r,HorizonSide((Side+1) % 4))
    end
    return Flag,y
end
function corner(r)
    Flag = true
    y = 0
    z = 0
    while !isborder(r,Sud) || Flag
        if !isborder(r,Sud)
            move!(r,Sud)
            y += 1
        end
        if isborder(r,Sud)
            Flag,y =  ob2(r,2,y)
        end
    end
    while !isborder(r,West)
        move!(r,West)
        z += 1
    end
    return y,z
end

function marks(r)
    Flag2 = true

    while !isborder(r,Nord) || Flag2
        Flag1 = true
        while !isborder(r,Ost) || Flag1
            putmarker!(r)
            if !isborder(r,Ost)
                move!(r,Ost)
            else
                Flag1 = ob(r,3)
            end
        end
        Flag1 = true
        while !isborder(r,West) || Flag1
            if !isborder(r,West)
                move!(r,West)
            else
                Flag1 = ob(r,1)
            end
        end
        if !isborder(r,Nord)
            move!(r,Nord)
        else
            Flag2 = ob(r,0)
        end
    end

end

function obratno(r,y,z)
    Flag = true
    while !isborder(r,Sud)
        move!(r,Sud)
    end
    while !isborder(r,West)
        move!(r,West)
    end
    for i = 1:z
        move!(r,Ost)
    end
    while y > 0
        if !isborder(r,Nord)
            move!(r,Nord)
            y -= 1
        else
            Flag, y = ob3(r,0,y)
        end
    end

end

function task_16(r)
    y = 0
    z = 0
    y, z = corner(r)
    marks(r)
    obratno(r,y,z)
end

task_16(r)


text:
module TestNotImplemented

using GroupsCore
using Test

struct SomeGroup <: Group end

struct SomeGroupElement <: GroupElement
    elts::Vector{Int} # SomeGroupElement is not isbits anymore
end

if VERSION < v"1.5.0"
    contains(haystack, needle) = occursin(needle, haystack)
end

@testset "Exceptions" begin
    @testset "InterfaceNotImplemented exception" begin
        @test GroupsCore.InterfaceNotImplemented(
            :Castle,
            "Aaaaarghhhhh....",
        ) isa Exception
        ex = GroupsCore.InterfaceNotImplemented(:Castle, "Aaargh...")
        @test ex isa GroupsCore.InterfaceNotImplemented

        @test sprint(showerror, ex) ==
              "Missing method from Castle interface: `Aaargh...`"
    end
    @testset "InfiniteOrder" begin
        G = SomeGroup()
        @test contains(
            sprint(showerror, GroupsCore.InfiniteOrder(G)),
            "isfinite",
        )
        g = SomeGroupElement(Int[1, 2, 3])
        @test contains(
            sprint(showerror, GroupsCore.InfiniteOrder(g)),
            "isfiniteorder",
        )
        @test contains(
            sprint(showerror, GroupsCore.InfiniteOrder(g, "Aaargh...")),
            "Aaargh...",
        )
    end
end

@testset "SomeGroup: No interface implemented" begin

    INI = GroupsCore.InterfaceNotImplemented
    InfO = GroupsCore.InfiniteOrder

    @testset "Group Interface" begin

        G = SomeGroup()

        # Iteration
        @test_throws INI eltype(G)
        @test_throws INI iterate(G)
        @test_throws INI iterate(G, 1)

        GroupsCore.hasgens(::SomeGroup) = false

        @test_throws ArgumentError iterate(G)
        @test_throws ArgumentError iterate(G, 1)
        @test_throws ArgumentError gens(G, 1)

        # revert to the default
        GroupsCore.hasgens(::SomeGroup) = true

        # Assumption 1: Groups are of unknown size
        @test Base.IteratorSize(G) == Base.SizeUnknown()
        @test_throws ArgumentError Base.isfinite(G)
        @test_throws ArgumentError order(G)

        Base.IteratorSize(::Type{SomeGroup}) = Base.HasShape{1}()
        @test Base.isfinite(G)
        @test_throws INI order(G)

        Base.IteratorSize(::Type{SomeGroup}) = Base.HasLength()
        @test Base.isfinite(G)
        @test_throws INI order(G)

        Base.IteratorSize(::Type{SomeGroup}) = Base.IsInfinite()
        @test !Base.isfinite(G)
        @test_throws InfO order(G)

        # return to the default:
        Base.IteratorSize(::Type{SomeGroup}) = Base.SizeUnknown()

        # Assumption 2: Groups have generators:
        @test hasgens(G)

        # Group Interface
        @test_throws INI one(G)
        @test_throws ArgumentError order(G)
        @test_throws INI gens(G)

        Base.eltype(::Type{SomeGroup}) = SomeGroupElement
        @test_throws INI rand(G, 2)

        @test_throws INI gens(G, 1)
        @test_throws INI ngens(G)

        @test_throws INI GroupsCore.rand_pseudo(G, 2, 2)
    end

    @testset "GroupElem Interface" begin

        g = SomeGroupElement(Int[])

        @test_throws INI parent(g)
        @test_throws INI g == g
        @test_throws INI isequal(g, g)

        @test_throws INI isfiniteorder(g)

        Base.IteratorSize(::SomeGroup) = Base.HasLength()
        Base.parent(::SomeGroupElement) = SomeGroup()
        @test isfiniteorder(g)
        Base.IteratorSize(::SomeGroup) = Base.SizeUnknown()

        GroupsCore.isfiniteorder(::SomeGroupElement) = false
        @test_throws InfO order(g)
        @test_throws INI deepcopy(g)

        @test_throws INI inv(g)
        @test_throws INI g * g
    end

end

end # of module TestNotImplemented


text:
using OVERTVerify
using LazySets
using Dates
using JLD2

function run_query(query_number, avoid_set, controller_name)
	controller = "nnet_files/jmlr/tora_"*controller_name*"_controller.nnet"

	println("Controller: ", controller_name)

	query = OvertQuery(
		Tora,      # problem
		controller, # network file
		Id(),    # last layer activation layer Id()=linear, or ReLU()=relu
		"MIP",     # query solver, "MIP" or "ReluPlex"
		15,        # ntime
		0.1,       # dt
		-1,        # N_overt
		)

	input_set = Hyperrectangle(low=[0.6, -0.7, -0.4, 0.5], high=[0.7, -0.6, -0.3, 0.6])
	# the other Infs are just to denote that these constraints should be ignored. So essentially the only constraint being added is:
	# x1 <= -2
	t1 = Dates.time()
	SATus, vals, stats = symbolic_satisfiability(query, input_set, avoid_set)
	t2 = Dates.time()
	dt = (t2-t1)
	println("dt is $dt")

	JLD2.@save "src/examples/jmlr/data/tora_satisfiability_"*string(controller_name)*"_controller_data_q"*string(query_number)*".jld2" query input_set avoid_set SATus vals stats dt query_number controller_name

	return SATus
end

function run_tora_satisfiability(;controller_name="smallest")

	# query 1
	avoid_set1 = HalfSpace([1., 0., 0., 0.], -2.) # checks if x1 <= -2
	SATus1 = run_query(1, avoid_set1, controller_name)

	if SATus1 == "unsat" # early stopping, potentially, if first query is sat
		# query 2
		avoid_set2 =  HalfSpace([-1., 0., 0., 0.], -2.)# checks if x1 >= 2
		SATus2 = run_query(2, avoid_set2, controller_name)

		open("src/examples/jmlr/data/tora_satisfiability_"*string(controller_name)*".txt", "w") do io
			write(io, "SATus1 = $SATus1 \n SATus2 = $SATus2")
		end;
	end
end

run_tora_satisfiability(controller_name=ARGS[1])


text:
module SteadyStateFit

export Optim, SteadyStateObjective, solve

import DiffEqBase: solve  # use CommonSolve.jl?
import Kaleido
import NLsolve
import Optim
using DiffEqBase: ODEProblem
using LinearAlgebra: lu, pinv, det
using Setfield
using ZygoteRules: @adjoint

if isdefined(Kaleido, :prettylens)
    using Kaleido: prettylens
else
    prettylens(x; kwargs...) = sprint(print, x; kwargs...)
end

include("utils.jl")
include("znlsolve.jl")
include("objectives.jl")
include("solvers.jl")

function __init__()
    # Load Zygote only at run-time:
    if ccall(:jl_generating_output, Cint, ()) != 1
        Zygote = Base.require(Base.PkgId(
            Base.UUID("e88e6eb3-aa80-5325-afca-941959d7151f"),
            "Zygote",
        ))
        @eval const pullback = $Zygote.pullback
    end
end

end # module


text:
module TestSKCrossValidator

using Test
using Random
using TSMLextra
using TSML: @pipeline

function test_skcrossvalidator()
  data=getiris()
  X=data[:,1:4]
  Y=data[:,5] |> Vector{String}
  rf = RandomForest()
  ppl1 = @pipeline rf
  @test crossvalidate(ppl1,X,Y,"accuracy_score").mean > 0.80
  @test crossvalidate(ppl1,X,Y,"balanced_accuracy_score").mean > 0.80
  @test crossvalidate(ppl1,X,Y,"cohen_kappa_score").mean > 0.80
  @test crossvalidate(ppl1,X,Y,"jaccard_score","weighted").mean > 0.80
  @test crossvalidate(ppl1,X,Y,"matthews_corrcoef").mean > 0.80
  @test crossvalidate(ppl1,X,Y,"hamming_loss").mean < 0.1
  @test crossvalidate(ppl1,X,Y,"zero_one_loss").mean < 0.1
  @test crossvalidate(ppl1,X,Y,"f1_score","weighted").mean > 0.80
  @test crossvalidate(ppl1,X,Y,"precision_score","weighted").mean > 0.80
  @test crossvalidate(ppl1,X,Y,"recall_score","weighted").mean > 0.80
end
@testset "CrossValidator" begin
  Random.seed!(123)
  test_skcrossvalidator()
end

end

max_stars_repo_name: ppalmes/TSMLextra

max_issues_repo_name: ppalmes/TSMLextra

max_forks_repo_name: ppalmes/TSMLextra


text:
# %% Compute the scalar and tensor CMB power spectrum
# Call the python wrapper for CLASS, `classy`
# Reference: https://github.com/lesgourg/class_public
classy = pyimport("classy")
# Create an instance of the CLASS wrapper
cosmo = classy.Class()
## Scalar power spectrum
# Create a dictionary of the cosmological parameters
Tcmb = 2.7255e6 # in μK
deg_ncdm = true # 3 massive neutrinos with degenerate mass?
params = Dict(
    "output" => "tCl pCl lCl",
    "modes" => "s",
    "l_max_scalars" => lmax,
    "lensing" => "yes",
    "A_s" => exp(3.0448) * 1e-10,
    "n_s" => 0.96605,
    "k_pivot" => 0.05,
    "h" => 0.6732,
    "omega_b" => 0.022383,
    "omega_cdm" => 0.12011,
    "tau_reio" => 0.0543,
    "N_ncdm" => 1,
)
if deg_ncdm
    push!(params, "m_ncdm" => 0.02, "deg_ncdm" => 3, "N_ur" => 0.00641)
    mν = params["m_ncdm"] * params["deg_ncdm"]
else
    push!(params, "m_ncdm" => 0.06, "N_ur" => 2.0328)
    mν = params["m_ncdm"]
end
# Set the parameters to the cosmological code
cosmo.set(params)
# Run the whole code. Depending on your output, it will call the
# CLASS modules more or less fast. For instance, without any
# output asked, CLASS will only compute background quantities,
# thus running almost instantaneously.
# This is equivalent to the beginning of the `main` routine of CLASS,
# with all the struct_init() methods called.
cosmo.compute()
# Access the lensed scalar cl until l=lmax. It is a dictionnary that contains the fields: ell, tt, te, ee, bb, pp, tp
cls = cosmo.lensed_cl(lmax)
# If you want to change completely the cosmology, you should also
# clean the arguments, otherwise, if you are simply running on a loop
# of different values for the same parameters, this step is not needed
cosmo.empty()
## Tensor power spectrum
rclass = 0.01
params = Dict(
    "output" => "tCl pCl",
    "modes" => "t",
    "l_max_tensors" => lmax,
    "lensing" => "no",
    "r" => rclass,
    "n_t" => 0,
    "A_s" => exp(3.0448) * 1e-10,
    "k_pivot" => 0.05,
    "h" => 0.6732,
    "omega_b" => 0.022383,
    "omega_cdm" => 0.12011,
    "tau_reio" => 0.0543,
    "N_ncdm" => 1,
)
if deg_ncdm
    push!(params, "m_ncdm" => 0.02, "deg_ncdm" => 3, "N_ur" => 0.00641)
    mν = params["m_ncdm"] * params["deg_ncdm"]
else
    push!(params, "m_ncdm" => 0.06, "N_ur" => 2.0328)
    mν = params["m_ncdm"]
end
cosmo.set(params)
cosmo.compute()
clt = cosmo.raw_cl(lmax)
cosmo.empty()


text:
module BFloat16s

    export BFloat16, LowPrecArray

    include("bfloat16.jl")
    include("lowprecarrays.jl")

end


text:
"""
    encode_motif(x, m::Int = length(x)) → s::Int

Encode the length-`m` motif `x` (a vector of indices that would sort some vector `v`
in ascending order) into its unique integer symbol ``s \\in \\{1, 2, \\ldots, m - 1 \\}``,
using Algorithm 1 in Berger et al. (2019)[^Berger2019].

## Example

```julia
v = rand(5)

# The indices that would sort `v` in ascending order. This is now a permutation
# of the index permutation (1, 2, ..., 5)
x = sortperm(v)

# Encode this permutation as an integer.
encode_motif(x)
```
[^Berger2019]: Berger, Sebastian, et al. "Teaching Ordinal Patterns to a Computer: Efficient Encoding Algorithms Based on the Lehmer Code." Entropy 21.10 (2019): 1023.
"""
function encode_motif(x, m::Int = length(x))
    n = 0
    for i = 1:m-1
        for j = i+1:m
            n += x[i] > x[j] ? 1 : 0
        end
        n = (m-i)*n
    end

    return n
end

function isless_rand(a, b)
    if a == b
        rand(Bool)
    elseif a < b
        true
    else
        false
    end
end

""" Compute probabilities of symbols `Π`, given weights `wts`. """
function symprobs(Π::AbstractVector, wts::AbstractVector; normalize = true)
    length(Π) == length(wts) || error("Need length(Π) == length(wts)")
    N = length(Π)
    idxs = sortperm(Π)
    sΠ = Π[idxs]   # sorted symbols
    sw = wts[idxs] # sorted weights

    i = 1   # symbol counter
    W = 0.0 # Initialize weight
    ps = Float64[]

    prev_sym = sΠ[1]

    while i <= length(sΠ)
        symᵢ = sΠ[i]
        wtᵢ = sw[i]
        if symᵢ == prev_sym
            W += wtᵢ
        else
            # Finished counting weights for the previous symbol, so push
            # the summed weights (normalization happens later).
            push!(ps, W)

            # We are at a new symbol, so refresh sum with the first weight
            # of the new symbol.
            W = wtᵢ
        end
        prev_sym = symᵢ
        i += 1
    end
    push!(ps, W) # last entry

    # Normalize
    Σ = sum(sw)
    if normalize
        return ps ./ Σ
    else
        return ps
    end
end


text:
using LinearAlgebra
@testset "merge" begin
    A = [3 4; 5 6.6]

    # Example which constructs a composition
    x = [3.0; 2.2; 0.1]
    # Generate a polynomial
    (graph1, cref) = graph_ps(x)
    g1_output = graph1.outputs[end]

    # Generate graph x^2
    (graph2, cref) = graph_monomial([0.0; 0; 1.0], input = :B)

    graph = merge_graphs(graph1, graph2, prefix2 = "Monomial")

    # Point the output of graph1 to the input of graph2
    add_lincomb!(graph, :MonomialB, 1.0, g1_output, 0, :I)

    # Adjust outputs
    newout = graph.outputs[end]
    empty!(graph.outputs)
    add_output!(graph, newout)

    @test (x[1] * I + x[2] * A + x[3] * A^2)^2 ≈ eval_graph(graph, A)

    A = [2.7 4; 5.3 6.6] / 7
    (graph, cref) = graph_newton_schulz(3)
    (graph_2, cref_2) = graph_newton_schulz_degopt(3)
    rename_node!(graph_2, graph_2.outputs[1], :Myout)
    (graph_3, cref_3) = graph_newton_schulz_degopt(3, input = :G2Myout)
    graph_4 = merge_graphs(
        graph_2,
        graph_3;
        prefix1 = "G2",
        prefix2 = "G3",
        input2 = :G2Myout,
    )
    @test eval_graph(graph_4, A) ≈ eval_graph(graph, eval_graph(graph, A))
end


text:
"""
Copyright Marcus Greiff 2020

The Aerial Vehicle Control package
"""
module AerialVehicleControl

export @verbose,
    recompile,
    get_matrix_and_copy,
    skew,
    SO3_hat,
    SO3_vee,
    SU2_hat,
    SU2_vee,
    quat_2_SU2,
    SU2_2_quat,
    quat_2_SO3,
    quaternion_left_product,
    attitude_dynamics,
    ntuple_2_mat,
    rand_PSD,
    allFieldsEqual,
    matrix_double_t,
    ref_state_qw_t,
    dyn_state_qw_t,
    con_state_qw_fsf_t,
    con_state_qw_fof_t,
    maximal_feasible_kc_SO3,
    maximal_feasible_kc_SU2,
    maximal_feasible_eps_SO3,
    get_M1M2W_SO3,
    get_M1M2W_SU2,
    numerical_differentiation

using LaTeXStrings
using LinearAlgebra
using Libdl
using Test

const CONT_LIB_BASE = dirname(@__DIR__)
const CONT_LIB_SRC  = joinpath(CONT_LIB_BASE, "src")
const CONT_LIB_UTIL = joinpath(CONT_LIB_BASE, "utils")
const CONT_LIB_PATH = joinpath(CONT_LIB_SRC,  "run")
const CONT_LIB_DOCS = joinpath(CONT_LIB_BASE, joinpath("docs", "images"))

# Recompilation of the C-code
include(joinpath(CONT_LIB_UTIL, "cont_compile.jl"))
# Julia types mirroring the structs in C (see cont_types.c)
include(joinpath(CONT_LIB_UTIL, "cont_types.jl"))
# The help functions used in the tests and examples
include(joinpath(CONT_LIB_UTIL, "cont_help_functions.jl"))

end