text:
"
Author: Daniel O'Connor
Email: dt.oconnor@outlook.com

"

module SpinVectorO2Dynamics
    include("..\\samplers\\tools.jl")

    using ..SamplingTools: hj_info, add_noise, add_cross_talk, rotor_aggregation
    using DifferentialEquations
    using Random
    using LinearAlgebra
    using Interpolations
    using RecursiveArrayTools

    function O2_model(spin_vector::Array{Float64, 1},
                             biases::Array{Float64, 1},
                             JM::Array{Float64, 2},
                             run_time::Float64,
                             open_system::Bool, # stochastic only
                             temp::Float64, # stochastic only
                             damping::Float64, # stochastic only
                             trajectories::Int64, # stochastic only
                             at_fn,
                             bt_fn,
                             ts,
                             abstol::Float64,
                             reltol::Float64)
        nspins = size(spin_vector, 1)
        ω_0 = zeros(Float64, nspins)
        u_0 = ArrayPartition(spin_vector, ω_0)
        t_span = (0.0, run_time)

        model! = function (du, u, p, t)
            # θ'(t) = ω(t)
            du.x[1] .= u.x[2]
            # α/I = ω'(t) = -1/I dV/dϑ = F/I
            du.x[2] .= (at_fn(t) .* cos.(u.x[1])) +
                       (bt_fn(t) .* ((biases .* sin.(u.x[1])) +
                                     (sin.(u.x[1]) .* (JM * cos.(u.x[1])))))
            nothing
        end

        model_jac! = function (J, u, p, t)
            # Jacobian J[i, j] = df_i / du_j
            J[1, 1, :] .= 0.0
            J[1, 2, :] .= 1.0
            J[2, 2, :] .= 0.0
            J[2, 1, :] .= (-at_fn(t) .* sin.(u.x[1])) +
                         (-bt_fn(t) .* ((biases .* cos.(u.x[1])) +
                                         (cos.(u.x[1]) .* (JM * cos.(u.x[1])))))
            nothing
        end

        noise_model! = function (du, u, p, t)
            # θ'(t) = ω(t)
            du.x[1] .= u.x[2]
            # ω'(t) = - dV/dϑ - λω = F/I
            du.x[2] .= (-damping .* u.x[2]) +
                       (at_fn(t) .* cos.(u.x[1])) +
                       (bt_fn(t) .* ((biases .* sin.(u.x[1])) +
                                     (sin.(u.x[1]) .* (JM * cos.(u.x[1])))))
            nothing
        end

        noise_model_jac! = function (J, u, p, t)
            # Jacobian J[i, j] = df_i / du_j
            J[1, 1, :] .= 0.0
            J[1, 2, :] .= 1.0
            J[2, 2, :] .= -damping
            J[2, 1, :] .= (-at_fn(t) .* sin.(u.x[1])) +
                         (-bt_fn(t) .* ((biases .* cos.(u.x[1])) +
                                         (cos.(u.x[1]) .* (JM * cos.(u.x[1])))))
            nothing
        end

        σ_model! = function (du, u, p, t)
            # constants infront of wiener process
            du.x[1] .= 0.0
            du.x[2] .= sqrt(2*damping*temp)
            nothing
        end

        # supplying jacobian can provide speed increases, especially for stiff
        # problems, use commented line to remove jacobian version

        if open_system
            f =  DifferentialEquations.ODEFunction(noise_model!, jac=noise_model_jac!)
            # f = noise_model!
            prob = DifferentialEquations.SDEProblem(f, σ_model!, u_0, t_span)
            ensembleprob = EnsembleProblem(prob)
            # can be upgraded to use GPU
            sol = DifferentialEquations.solve(ensembleprob, EnsembleThreads(),
                                              trajectories=trajectories,
                                              saveat=ts)
        else
            # f =  DifferentialEquations.ODEFunction(model!, jac=model_jac!)
            f = model!
            prob = DifferentialEquations.ODEProblem(f, u_0, t_span)
            sol = DifferentialEquations.solve(prob, saveat=ts, abstol=abstol,
                                              reltol=reltol)
        end

        return sol
    end

    function solve(h::Dict, J::Dict; anneal_time=100.0, run_time=0.0, # in nanoseconds
                    temp=0.1, friction_constant=1e-3, trajectories=1000,
                    # for individual qubit schedules, input a
                    # dictionary of interpolations
                    # NOTE if interpolations aren't continuous then
                    # instabilities arise in the ODE integrator
                    a_fn=nothing, b_fn=nothing, saveat=collect(0:0.01:1),
                    add_xtalk=false, noisy=false, open_system=false,
                    initial_state=Array{Float64, 1}(), abstol=1e-8,
                    reltol=1e-6)

        nspins = length(h)
        # check to see if we have individual qubit schedules
        independent_qubit_sechdules = false
        if isa(a_fn, Dict) | isa(b_fn, Dict)
            # if one input is not a dict but one if, convert to dict of lists
            if !isa(a_fn, Dict)
                interp = deepcopy(a_fn)
                a_fn = Dict([(i, interp) for i = 1:nspins])
            elseif !isa(b_fn, Dict)
                interp = deepcopy(b_fn)
                b_fn = Dict([(i, interp) for i = 1:nspins])
            end

            independent_qubit_sechdules = true
        elseif a_fn == nothing
            # linear schedules
            a_fn = LinearInterpolation(saveat, 1.0 .- saveat, extrapolation_bc=Flat())
        elseif b_fn == nothing
            b_fn = LinearInterpolation(saveat, saveat, extrapolation_bc=Flat())
        end

        if run_time == 0
            run_time = anneal_time
        end

        # define coefficient interpolation functions
        if independent_qubit_sechdules
            # vectorise the independent_schedules
            at_fn = function (t)
                a_vec = zeros(Float64, nspins)
                for i = 1:nspins
                    a_vec[i] = a_fn[i](t/anneal_time)
                end
                return a_vec
            end
            bt_fn = function (t)
                b_vec = zeros(Float64, nspins)
                for i = 1:nspins
                    b_vec[i] = b_fn[i](t/anneal_time)
                end
                return b_vec
            end
        else
            at_fn = (t) -> a_fn(t/anneal_time)
            bt_fn = (t) -> b_fn(t/anneal_time)
        end

        # prepare initial state if defined
        if length(initial_state) == nspins
            @assert ndims(initial_state) == 1
            spin_vector = acos.(deepcopy(initial_state))
        else
            # starting in the transverse field direction as default
            spin_vector = fill(0.5 * pi, nspins)
        end

        # add noise effects
        if noisy
            h, J = add_noise(h, J)
        end
        # add cross-talk effects
        if add_xtalk
            h, J = add_xtalk(h, J)
        end

        # prepare problem info
        JM, M, degs = hj_info(h, J)
        JS = deepcopy(JM)
        biases = diag(JM)
        JM -= diagm(biases)
        JM += JM'

        # run the protocol
        sol = O2_model(spin_vector, biases, JM, run_time,
                      open_system, # stochastic only
                      temp, # stochastic only
                      friction_constant, # stochastic only
                      trajectories, # stochastic only
                      at_fn, bt_fn, run_time.*saveat, abstol, reltol)

        if open_system
            spin_vectors = [sol[i].u[end].x[1] for i =1:trajectories]
            results = rotor_aggregation(spin_vectors, trajectories, collect(UpperTriangular(JS)))
        else
            final_sol = sol.u[end].x[1]
            spins = Array{Float64, 1}(sign.(cos.(final_sol)))
            energy = biases' * spins + 0.5 * (spins' * (JM * spins))
            results = [spins, energy]
        end
        return results, sol
    end

    export solve, O2_model

end

max_stars_repo_name: dtoconnor/Annealing-Samplers

max_issues_repo_name: dtoconnor/Annealing-Samplers

max_forks_repo_name: dtoconnor/Annealing-Samplers


text:
#####
##### Fields
#####

##### Auxiliary fields

# Center only
cent_aux_vars_en_2m(FT) = (;
    dissipation = FT(0),
    shear = FT(0),
    entr_gain = FT(0),
    detr_loss = FT(0),
    press = FT(0),
    buoy = FT(0),
    interdomain = FT(0),
    rain_src = FT(0),
)
cent_aux_vars_up(FT) = (;
    q_liq = FT(0),
    q_ice = FT(0),
    T = FT(0),
    RH = FT(0),
    s = FT(0),
    buoy = FT(0),
    area = FT(0),
    q_tot = FT(0),
    θ_liq_ice = FT(0),
    θ_liq_ice_tendency_precip_formation = FT(0),
    qt_tendency_precip_formation = FT(0),
    entr_sc = FT(0),
    detr_sc = FT(0),
    frac_turb_entr = FT(0),
    entr_turb_dyn = FT(0),
    detr_turb_dyn = FT(0),
    asp_ratio = FT(0),
    horiz_K_eddy = FT(0),
)
cent_aux_vars_edmf(FT, n_up) = (;
    turbconv = (;
        ϕ_temporary = FT(0),
        ψ_temporary = FT(0),
        bulk = (;
            area = FT(0),
            θ_liq_ice = FT(0),
            RH = FT(0),
            buoy = FT(0),
            q_tot = FT(0),
            q_liq = FT(0),
            q_ice = FT(0),
            T = FT(0),
            cloud_fraction = FT(0),
            θ_liq_ice_tendency_precip_formation = FT(0),
            qt_tendency_precip_formation = FT(0),
            Ri = FT(0),
        ),
        up = ntuple(i -> cent_aux_vars_up(FT), n_up),
        en = (;
            w = FT(0),
            area = FT(0),
            q_tot = FT(0),
            q_liq = FT(0),
            q_ice = FT(0),
            θ_liq_ice = FT(0),
            θ_virt = FT(0),
            θ_dry = FT(0),
            RH = FT(0),
            s = FT(0),
            T = FT(0),
            buoy = FT(0),
            cloud_fraction = FT(0),
            tke = FT(0),
            Hvar = FT(0),
            QTvar = FT(0),
            HQTcov = FT(0),
            qt_tendency_precip_formation = FT(0),
            θ_liq_ice_tendency_precip_formation = FT(0),
            unsat = (; q_tot = FT(0), θ_dry = FT(0), θ_virt = FT(0)),
            sat = (; T = FT(0), q_vap = FT(0), q_tot = FT(0), θ_dry = FT(0), θ_liq_ice = FT(0)),
            Hvar_rain_dt = FT(0),
            QTvar_rain_dt = FT(0),
            HQTcov_rain_dt = FT(0),
        ),
        θ_liq_ice_tendency_precip_sinks = FT(0),
        qt_tendency_precip_sinks = FT(0),
        qr_tendency_evap = FT(0),
        qs_tendency_melt = FT(0),
        qs_tendency_dep_sub = FT(0),
        qr_tendency_advection = FT(0),
        qs_tendency_advection = FT(0),
        en_2m = (;
            tke = cent_aux_vars_en_2m(FT),
            Hvar = cent_aux_vars_en_2m(FT),
            QTvar = cent_aux_vars_en_2m(FT),
            HQTcov = cent_aux_vars_en_2m(FT),
        ),
        term_vel_rain = FT(0),
        term_vel_snow = FT(0),
        KM = FT(0),
        KH = FT(0),
        mixing_length = FT(0),
        θ_virt = FT(0),
        massflux_tendency_h = FT(0),
        massflux_tendency_qt = FT(0),
        diffusive_tendency_h = FT(0),
        diffusive_tendency_qt = FT(0),
        prandtl_nvec = FT(0),
        # Added by Ignacio : Length scheme in use (mls), and smooth min effect (ml_ratio)
        # Variable Prandtl number initialized as neutral value.
        mls = FT(0),
        b_exch = FT(0),
        ml_ratio = FT(0),
        w_up_c = FT(0),
        w_en_c = FT(0),
        Shear² = FT(0),
        ∂θv∂z = FT(0),
        ∂qt∂z = FT(0),
        ∂θl∂z = FT(0),
        ∂θv∂z_unsat = FT(0),
        ∂qt∂z_sat = FT(0),
        ∂θl∂z_sat = FT(0),
        l_entdet = FT(0),
        ϕ_gm = FT(0), # temporary for grid-mean variables
        ϕ_gm_cov = FT(0), # temporary for grid-mean covariance variables
        ϕ_en_cov = FT(0), # temporary for environmental covariance variables
        ϕ_up_cubed = FT(0), # temporary for cubed updraft variables in grid mean 3rd moment functions
    ),
)

# Face only
face_aux_vars_up(FT) = (;
    w = FT(0),
    nh_pressure = FT(0),
    nh_pressure_b = FT(0),
    nh_pressure_adv = FT(0),
    nh_pressure_drag = FT(0),
    massflux = FT(0),
)

face_aux_vars_edmf(FT, n_up) = (;
    turbconv = (;
        bulk = (; w = FT(0)),
        ρ_ae_KM = FT(0),
        ρ_ae_KH = FT(0),
        ρ_ae_K = FT(0),
        en = (; w = FT(0)),
        up = ntuple(i -> face_aux_vars_up(FT), n_up),
        massflux_h = FT(0),
        massflux_qt = FT(0),
        ϕ_temporary = FT(0),
        diffusive_flux_h = FT(0),
        diffusive_flux_qt = FT(0),
        diffusive_flux_u = FT(0),
        diffusive_flux_v = FT(0),
    ),
)

##### Diagnostic fields

# Center only
cent_diagnostic_vars_edmf(FT, n_up) = (;
    turbconv = (;
        asp_ratio = FT(0),
        entr_sc = FT(0),
        detr_sc = FT(0),
        massflux = FT(0),
        frac_turb_entr = FT(0),
        horiz_K_eddy = FT(0),
    ),
)

# Face only
face_diagnostic_vars_edmf(FT, n_up) =
    (; turbconv = (; nh_pressure = FT(0), nh_pressure_adv = FT(0), nh_pressure_drag = FT(0), nh_pressure_b = FT(0)))

##### Prognostic fields

# Center only
cent_prognostic_vars_up(FT) = (; ρarea = FT(0), ρaθ_liq_ice = FT(0), ρaq_tot = FT(0))
cent_prognostic_vars_en(FT) = (; ρatke = FT(0), ρaHvar = FT(0), ρaQTvar = FT(0), ρaHQTcov = FT(0))
cent_prognostic_vars_edmf(FT, n_up) = (;
    turbconv = (;
        en = cent_prognostic_vars_en(FT),
        up = ntuple(i -> cent_prognostic_vars_up(FT), n_up),
        pr = (; q_rai = FT(0), q_sno = FT(0)),
    ),
)
# cent_prognostic_vars_edmf(FT, n_up) = (;) # could also use this for empty model

# Face only
face_prognostic_vars_up(FT) = (; ρaw = FT(0))
face_prognostic_vars_edmf(FT, n_up) = (; turbconv = (; up = ntuple(i -> face_prognostic_vars_up(FT), n_up)))

ext: jl


text:
"""
    fastrange(h, n)

Map hash to to bucket in range 1:n. See [fastrange](https://github.com/lemire/fastrange).
"""
fastrange(h, n::Int) = (((h % UInt128) * (n % UInt128)) >> 64) % Int + 1

"""
    target_capacity(curr_length)

Compute the target capacity.
"""
target_capacity(curr_length) = cld(3 * curr_length, 2)

"""
    Token

Utility struct used in finding and modifying entries in a [`Storage`](@ref).

# Fields:

* `parent`: Index of the parent in the `pairs` array. Zero if token has no parent.
* `index`: Index of the token in the `pairs` array. Zero if the key is not yet in the
  vector.
* `pointer`: The index in the `pointers` array that points to the first element of the chain
  that should contain the token.
"""
struct Token
    parent::Int
    index::Int
    pointer::Int
end

struct Storage{K,V,T<:Signed}
    pointers::Vector{T}
    pairs::Vector{Pair{K,V}}
    next::Vector{T}
    maxlength::Base.RefValue{Int}
end
function Storage{K,V,T}(; capacity=100) where {K,V,T}
    pointers = fill(-one(T), capacity)
    pairs = Pair{K,V}[]
    next = T[]
    return Storage(pointers, pairs, next, Ref(0))
end
Storage{K,V}(; kwargs...) where {K,V} = Storage{K,V,Int32}(; kwargs...)

###
### Basic properties
###
Base.length(dv::Storage) = length(dv.next)
Base.isempty(dv::Storage) = isempty(dv.next)
Rimu.storage(dv::Storage) = dv
Rimu.capacity(dv::Storage) = length(dv.pointers)
Base.keytype(::Type{<:Storage{K}}) where {K} = K
Base.keytype(::Storage{K}) where {K} = K
Base.valtype(::Type{<:Storage{<:Any,V}}) where {V} = V
Base.valtype(::Storage{<:Any,V}) where {V} = V
Base.eltype(::Type{<:Storage{K,V}}) where {K,V} = Pair{K,V}
Base.eltype(::Storage{K,V}) where {K,V} = Pair{K,V}
indextype(::Type{<:Storage{<:Any,<:Any,T}}) where {T} = T
indextype(dv) = indextype(typeof(dv))

function Base.show(io::IO, ::MIME"text/plain", dv::Storage{K,V,T}) where {K,V,T}
    print(io, "$(length(dv))-element Storage{$K,$V,$T}\n")
    Base.print_array(io, dv.pairs)
end
function Base.show(io::IO, dv::Storage{K,V,T}) where {K,V,T}
    print(io, "$(length(dv))-element Storage")
end

###
### get/set/deposit
###
function Base.getindex(dv::Storage{K,V,T}, key::K, h=hash(key)) where {K,V,T}
    @inbounds index = dv.pointers[fastrange(h, length(dv.pointers))]
    if index < 0
        return zero(V)
    else
        @inbounds while true
            pair = dv.pairs[index]
            pair.first == key && return pair.second
            index = dv.next[index]
            index < 0 && return zero(V)
        end
    end
end

function update_maxlength!(dv::Storage)
    dv.maxlength[] = max(dv.maxlength[], length(dv))
end

"""
    get_token_by_key(dv::Storage{K}, key::K)

Find a [`Token`](@ref) that corresponds to the `key`.
"""
@inline function get_token_by_key(dv::Storage{K}, key::K, h) where {K}
    pointer = fastrange(h, length(dv.pointers))
    @inbounds index = dv.pointers[pointer]
    if index < 0
        return Token(0, 0, pointer)
    else
        @inbounds pair = dv.pairs[index]
        if pair.first == key
            return Token(0, index, pointer)
        else
            @inbounds while true
                parent = index
                index = dv.next[parent]
                index < 0 && return Token(parent, 0, pointer)
                pair = dv.pairs[index]
                pair.first == key && return Token(parent, index, pointer)
            end
        end
    end
end
"""
    get_token_by_index(dv::Storage, index)

Find a [`Token`](@ref) that corresponds to the `index` in the `pairs` array.
"""
@inline function get_token_by_index(dv::Storage, index::Integer)
    T = indextype(dv)
    @inbounds pointer = dv.next[index]
    @inbounds while pointer > 0
        pointer = dv.next[pointer]
    end
    pointer = -pointer
    @inbounds curr = dv.pointers[pointer]
    parent = 0
    @inbounds while curr ≠ (index % T)
        parent = curr
        curr = dv.next[curr]
    end
    return Token(parent, index, pointer)
end
"""
    update_parent!(dv::Storage, t::Token, new::Integer)

Update the parent of token `t` to point to new. If `t` has no parent, the `pointers` array
is updated instead.
"""
@inline Base.@propagate_inbounds function update_parent!(dv::Storage, t::Token, new::Integer)
    T = indextype(dv)
    if iszero(t.parent)
        dv.pointers[t.pointer] = new % T
    else
        dv.next[t.parent] = new % T
    end
end
"""
    delete!(dv::Storage, t::Token)

Delete the entry of `dv` that corresponds to token `t`.
"""
@inline function Base.delete!(dv::Storage, t::Token)
    @inbounds if t.index == length(dv)
        # This is the last one, need only to pop it.
        pop!(dv.pairs)
        update_parent!(dv, t, pop!(dv.next))
    elseif !iszero(t.index)
        # Bridge the gap between parent and child.
        update_parent!(dv, t, dv.next[t.index])

        # Fill in the gap with the last entry of `dv`.
        l = get_token_by_index(dv, length(dv))
        update_parent!(dv, l, t.index)
        dv.pairs[t.index] = pop!(dv.pairs)
        dv.next[t.index] = pop!(dv.next)
    end
end
"""
    update!(dv::Storage, t::Token, pair)

Update the entry in `dv` that corresponds to token `t` to contain `pair`. This function
handles the decision of whether a new entry must be created, or if an existing entry must
be updated or deleted.
"""
@inline function update!(dv::Storage, t::Token, pair)
    @inbounds if iszero(pair.second)
        # Zero - delete entry.
        !isempty(dv) && delete!(dv, t)
    elseif iszero(t.index)
        # New entry - add to the end of pairs.
        push!(dv.pairs, pair)
        push!(dv.next, -t.pointer)
        update_parent!(dv, t, length(dv))
        update_maxlength!(dv)
        maybe_rehash_grow!(dv)
    else
        # Update in-place.
        dv.pairs[t.index] = pair
    end
end
"""
    maybe_swap_parent!(dv, index, parent, pair)

Swap current index and its parent if parent's value is smaller.
Not used.
"""
function maybe_swap_parent!(dv, index, parent, pair)
    if length(dv) > 1_000_000 && parent ≠ 0
        parent_pair = dv.pairs[parent]
        if abs(parent_pair.second) < abs(pair.second)
            dv.pairs[index] = parent_pair
            dv.pairs[parent] = pair
        end
    end
end

function Base.setindex!(dv::Storage{K,V}, v, key, h=hash(key)) where {K,V}
    value = V(v)
    token = get_token_by_key(dv, key, h)
    update!(dv, token, key => value)
    return value
end

function Rimu.deposit!(dv::Storage{K,V}, key, v, h=hash(key), parent=nothing) where {K,V}
    value = V(v)
    token = get_token_by_key(dv, key, h)
    if !iszero(token.index)
        oldkey, oldvalue = dv.pairs[token.index]
        value += oldvalue
    end
    update!(dv, token, key => value)
    return value
end

###
### empty!, empty, similar, rehash...
###
function Base.empty!(dv::Storage)
    # Note: empty! changes the number of pointers in dv to have enough space to comfortably
    # fit its previous length. This can be done because we empty the vector after every
    # FCIQMC step
    n_target = target_capacity(dv.maxlength[])
    empty!(dv.pairs)
    empty!(dv.next)
    if n_target > length(dv.pointers)
        resize!(dv.pointers, n_target)
    end
    fill!(dv.pointers, -1)
    dv.maxlength[] = 0
    return dv
end
function Base.empty(dv::Storage{K,V,T}) where {K,V,T}
    return Storage{K,V,T}()
end
function Base.empty(dv::Storage{K,<:Any,T}, ::Type{V}) where {K,V,T}
    return Storage{K,V,T}()
end
function Base.empty(dv::Storage{<:Any,<:Any,T}, ::Type{K}, ::Type{V}) where {K,V,T}
    return Storage{K,V,T}()
end

function rehash!(dv::Storage)
    dv.next .= .-fastrange.(hash.(first.(dv.pairs)), length(dv.pointers))
    dv.pointers .= -1
    @inbounds for i in 1:length(dv)
        pointer = -dv.next[i]
        parent = 0
        index = dv.pointers[pointer]

        while index > 0
            parent = index
            index = dv.next[index]
        end

        update_parent!(dv, Token(parent, i, pointer), i)
    end
    return dv
end

function maybe_rehash_grow!(dv::Storage)
    if length(dv.pairs) > length(dv.pointers)
        resize!(dv.pointers, target_capacity(length(dv)))
        rehash!(dv)
    end
    return dv
end

function Base.sizehint!(dv::Storage, n)
    n_target = target_capacity(n)
    if n_target > length(dv.pointers)
        resize!(dv.pointers, n_target)
        rehash!(dv)
    end
    return dv
end

function Base.copy!(dst::Storage{K,V}, src::Storage{K,V}) where {K,V}
    copy!(dst.pointers, src.pointers)
    copy!(dst.pairs, src.pairs)
    copy!(dst.next, src.next)
    return dst
end

"""
    pairs(dv::Storage)

Iterates over pairs in `dv`.
"""
struct StoragePairs{K,V,L<:Storage{K,V}} <: AbstractVector{Pair{K,V}}
    dvec::L
end
Base.size(ps::StoragePairs) = (length(ps.dvec),)
Base.getindex(ps::StoragePairs, i) = ps.dvec.pairs[i]
Base.pairs(dv::Storage) = StoragePairs(dv)
function Base.summary(io::IO, ps::StoragePairs{K,V}) where {K,V}
    print(io, "$(length(ps))-element StoragePairs{$K,$V}")
end

"""
    keys(dv::Storage)

Iterates over keys in `dv`.
"""
struct StorageKeys{K,L<:Storage{K}} <: AbstractVector{K}
    dvec::L
end
Base.size(ks::StorageKeys) = (length(ks.dvec),)
Base.getindex(ks::StorageKeys, i) = ks.dvec.pairs[i].first
Base.keys(dv::Storage) = StorageKeys(dv)
function Base.summary(io::IO, ks::StorageKeys{K}) where {K}
    print(io, "$(length(ks))-element StorageKeys{$K}")
end

"""
    values(dv::Storage)

Iterates over values in `dv`.
"""
struct StorageValues{V,L<:Storage{<:Any,V}} <: AbstractVector{V}
    dvec::L
end
Base.size(vs::StorageValues) = (length(vs.dvec),)
Base.getindex(vs::StorageValues, i) = vs.dvec.pairs[i].second
Base.values(dv::Storage) = StorageValues(dv)
function Base.summary(io::IO, vs::StorageValues{V}) where {V}
    print(io, "$(length(vs))-element StorageValues{$V}")
end

"""
    map!(f, values(dv::Storage))
    map!(f, dst::Storage, values(src::Storage))

Apply function to all values in `dv`.
When mapping over `pairs`, the function should return the value type.
"""
function Base.map!(f, vs::StorageValues)
    dv = vs.dvec
    # Iterating in reverse order makes sure replacement had f already applied when deleting.
    @inbounds for i in length(dv.pairs):-1:1
        k, v = dv.pairs[i]
        v = f(v)
        if iszero(v)
            t = get_token_by_index(dv, i)
            delete!(dv, t)
        else
            dv.pairs[i] = k => v
        end
    end
end
function Base.map!(f, dst::Storage, vs::StorageValues)
    src = vs.dvec
    if src === dst
        map!(f, vs)
    else
        empty!(dst)
        @inbounds for i in length(src.pairs):-1:1
            k, v = src.pairs[i]
            dst[k] = f(v)
        end
    end
    return dst
end

function Rimu.add!(dst::Storage, src::Storage, α=one(valtype(src)))
    for (k, v) in pairs(src)
        deposit!(dst, k, v * α)
    end
    return dst
end

function LinearAlgebra.dot(v::Storage, w::Storage)
    result = zero(promote_type(valtype(v), valtype(w)))
    for (key, val) in pairs(v)
        result += conj(val) * w[key]
    end
    return result
end

function MPI.Buffer(v::Storage)
    return MPI.Buffer(v.pairs)
end


text:
include("util.jl") # data(), addgrad!()

# Gradient calculation methods for non broadcasted ops
# Disptach on op type
function calculategrad!(::Type{typeof(*)},
                        t::Tensor, u::Union{Tensor,AbstractArray}, res::Tensor)
    t.grad .+= res.grad * data(u)'
    addgrad!(u, t.data' * res.grad)
    return nothing
end

# Gradient calculation methods for broadcasted ops
# Dispatch on op type
import Base.Broadcast: BroadcastFunction

# Generate methods to calculate gradients for every order of arguments
for (T, S) in ((Tensor, Union{AbstractArray,Number}),
                (Union{AbstractArray,Number}, Tensor),
                (Tensor, Tensor)) # Needed to fix ambiguity
    @eval begin
        function calculategrad!(::Type{BroadcastFunction{typeof(+)}},
                                t::$T, u::$S, res::Tensor)
            addgrad!(t, res.grad)
            addgrad!(u, res.grad)
            return nothing
        end

        function calculategrad!(::Type{BroadcastFunction{typeof(-)}},
                                t::$T, u::$S, res::Tensor)
            addgrad!(t, res.grad)
            addgrad!(u, -res.grad)
            return nothing
        end

        function calculategrad!(::Type{BroadcastFunction{typeof(*)}},
                                t::$T, u::$S, res::Tensor)
            addgrad!(t, res.grad .* data(u))
            addgrad!(u, res.grad .* data(t))
            return nothing
        end
    end
end

max_stars_repo_name: leonvol/magrad

max_issues_repo_name: leonvol/magrad

max_forks_repo_name: leonvol/magrad


text:
export ad, mad, mad_core

"""
`mad_model(G)` returns the solved linear program whose optimum
value is `mad(G)`.
"""
function mad_model(G::SimpleGraph)
    EE = elist(G)
    VV = vlist(G)

    n = length(VV)
    m = length(EE)

    MOD = Model(get_solver())

    # variables
    @variable(MOD, z >= 0)
    @variable(MOD, x[i = 1:n, j = 1:m; in(VV[i], EE[j])] >= 0)

    # vertex constraints
    for i = 1:n
        v = VV[i]
        # get all edges that have v as an end point
        j_iterator = (j for j = 1:m if EE[j][1] == v || EE[j][2] == v)
        j_list = collect(j_iterator)
        @constraint(MOD, sum(x[i, j] for j in j_list) <= z)
    end

    # edge constraints
    for j = 1:m
        v, w = EE[j]
        i1 = first(findall(x -> x == v, VV))
        i2 = first(findall(x -> x == w, VV))
        @constraint(MOD, x[i1, j] + x[i2, j] == 2)
    end

    @objective(MOD, Min, z)
    optimize!(MOD)
    return MOD, x, VV, EE
end

"""
`mad(G)` computes the maximum average degree of `G`.
See also `mad_core`.
"""
function mad(G::SimpleGraph)
    if cache_check(G, :mad)
        return cache_recall(G, :mad)
    end
    MOD, x, VV, EE = mad_model(G)
    mval = objective_value(MOD)
    cache_save(G, :mad, mval)
    return mval
end

"""
`ad(G)` is the average degree of a vertex in `G`.
"""
ad(G::SimpleGraph) = 2 * NE(G) / NV(G)

"""
`mad_core(G)` returns a subgraph `H` of `G` for which
`ad(H)==mad(G)`.
"""
function mad_core(G::SimpleGraph)
    if cache_check(G, :mad_core)
        return cache_recall(G, :mad_core)
    end
    T = eltype(G)
    GG = deepcopy(G)
    while true
        if NV(GG) == 0
            error("Bad graph or something went wrong")
        end
        # solve the LP for GG
        MOD, x, VV, EE = mad_model(GG)
        n = length(VV)
        m = length(EE)
        err = 0.1 / n
        mval = objective_value(MOD)

        # if balanced, return
        if abs(ad(GG) - mval) <= err
            cache_save(G, :mad_core, GG)
            return GG
        end

        # otherwise, find and destroy defective vertices
        for i = 1:n
            v = VV[i]
            total = 0.0
            for j = 1:m
                e = EE[j]
                if in(v, e)
                    val = value(x[i, j])
                    total += val
                end
            end
            if total < mval - err   # this is a defective vertex
                delete!(GG, v)
            end
        end
    end # end while
    error("This can't happen (but it did)")
    cache_save(G, :mad_core, GG)
    return GG
end

max_stars_repo_name: scheinerman/SimpleGraphAlgorithms

max_issues_repo_name: scheinerman/SimpleGraphAlgorithms

max_forks_repo_name: scheinerman/SimpleGraphAlgorithms


text:
# Number differentiation

function gradient(f::Number)
  @inline function grad_f(x::Point)
    zero(return_type(outer,x,f))
  end
end

function ∇∇(f::Number)
  @inline function hess_f(x::Point)
    g = gradient(f)(x)
    gradient(g)(x)
  end
end

# Automatic differentiation of functions

function gradient(f::Function)
  @inline function _gradient(x)
    gradient(f,x)
  end
end

function symmetric_gradient(f::Function)
  @inline function _symmetric_gradient(x)
    symmetric_gradient(f,x)
  end
end

function divergence(f::Function)
  @inline function _divergence(x)
    divergence(f,x)
  end
end

function curl(f::Function)
  @inline function _curl(x)
    curl(f,x)
  end
end

function laplacian(f::Function)
  @inline function _laplacian(x)
    laplacian(f,x)
  end
end

# Specialization when x::Point

function gradient(f::Function,x::Point)
  gradient(f,x,return_value(f,x))
end

@inline function gradient(f::Function,x::Point,fx::Number)
  VectorValue(ForwardDiff.gradient(f,get_array(x)))
end

@inline function gradient(f::Function,x::Point,fx::VectorValue)
  TensorValue(transpose(ForwardDiff.jacobian(y->get_array(f(y)),get_array(x))))
end

function gradient(f::Function,x::Point,fx::MultiValue)
  @notimplemented
end

function symmetric_gradient(f::Function,x::Point)
  symmetric_part(gradient(f,x))
end

function divergence(f::Function,x::Point)
  divergence(f,x,return_value(f,x))
end

function divergence(f::Function,x::Point,fx)
  tr(gradient(f,x,fx))
end

function divergence(f::Function,x::Point,fx::TensorValue{2,2})
  g(x) = SVector(f(x).data)
  a = ForwardDiff.jacobian(g,get_array(x))
  VectorValue(
    a[1,1]+a[2,2],
    a[3,1]+a[4,2],
  )
end

function divergence(f::Function,x::Point,fx::TensorValue{3,3})
  g(x) = SVector(f(x).data)
  a = ForwardDiff.jacobian(g,get_array(x))
  VectorValue(
    a[1,1]+a[2,2]+a[3,3],
    a[4,1]+a[5,2]+a[6,3],
    a[7,1]+a[8,2]+a[9,3],
   )
end

function divergence(f::Function,x::Point,fx::SymTensorValue{2})
  g(x) = SVector(f(x).data)
  a = ForwardDiff.jacobian(g,get_array(x))
  VectorValue(
    a[1,1]+a[2,2],
    a[2,1]+a[3,2],
  )
end

function divergence(f::Function,x::Point,fx::SymTensorValue{3})
  g(x) = SVector(f(x).data)
  a = ForwardDiff.jacobian(g,get_array(x))
  VectorValue(
    a[1,1]+a[2,2]+a[3,3],
    a[2,1]+a[4,2]+a[5,3],
    a[3,1]+a[5,2]+a[6,3],
   )
end

function curl(f::Function,x::Point)
  grad2curl(gradient(f,x))
end

function laplacian(f::Function,x::Point)
  laplacian(f,x,return_value(f,x))
end

@inline function laplacian(f::Function,x::Point,fx::Number)
  tr(ForwardDiff.jacobian(y->ForwardDiff.gradient(f,y), get_array(x)))
end

@inline function laplacian(f::Function,x::Point,fx::VectorValue)
  A = length(x)
  B = length(fx)
  a = ForwardDiff.jacobian(y->transpose(ForwardDiff.jacobian(z->get_array(f(z)),y)), get_array(x))
  tr(ThirdOrderTensorValue{A,A,B}(Tuple(transpose(a))))
end

function laplacian(f::Function,x::Point,fx::MultiValue)
  @notimplemented
end

# Specialization when x::StaticVector

@inline gradient(f::Function,x::SVector) = gradient(f,Point(x))
@inline symmetric_gradient(f::Function,x::SVector) = symmetric_gradient(f,Point(x))
@inline divergence(f::Function,x::SVector) = divergence(f,Point(x))
@inline curl(f::Function,x::SVector) = curl(f,Point(x))
@inline laplacian(f::Function,x::SVector) = laplacian(f,Point(x))


text:
function expect_molecule_compartments!(assignment_probs::Matrix{Float64}, adjacent_points::Vector{Vector{Int}}, adjacent_weights::Vector{Vector{Float64}}, is_locked::BitVector)
    for i in 1:length(adjacent_points)
        if is_locked[i]
            continue
        end
        cur_weights = adjacent_weights[i]
        cur_points = adjacent_points[i]

        if length(cur_points) == 0
            continue
        end

        dense_sum = 0.0
        for ri in 1:size(assignment_probs, 1)
            c_d = 0.0
            for j in 1:length(cur_points)
                c_d += cur_weights[j] * assignment_probs[ri, cur_points[j]]
            end

            dense_sum += assignment_probs[ri, i] = exp(c_d)
        end

        assignment_probs[:, i] ./= dense_sum
    end
end

init_nuclei_cyto_compartments(args...; nuclei_genes::T, cyto_genes::T, kwargs...) where T<:Union{Vector{String}, Dict{String, Float64}} =
    init_molecule_compartments(args..., [nuclei_genes, cyto_genes]; kwargs...)

init_molecule_compartments(pos_data::Matrix{Float64}, genes::Vector{Int}, comp_genes::Vector{Vector{String}}; kwargs...) =
    init_molecule_compartments(pos_data, genes, [Dict(g => 1.0 for g in gs) for gs in comp_genes]; kwargs...)

function init_molecule_compartments(pos_data::Matrix{Float64}, genes::Vector{Int}, comp_genes::Vector{Dict{String, Float64}};
        gene_names::Vector{String}, scale::Float64, required_comps::Union{Vector{Int}, Nothing}=nothing)

    id_per_gene = Dict(g => i for (i,g) in enumerate(gene_names))
    comp_genes = [Dict(id_per_gene[g] => v for (g,v) in gsd) for gsd in comp_genes]

    default_probs = vcat(ones(length(comp_genes)) ./ length(comp_genes), 0.0)
    assignment_probs = zeros(length(comp_genes) + 1, length(genes))
    assignment_probs .= default_probs

    is_assigned = [falses(size(pos_data, 2)) for c in comp_genes]
    for (i,gs) in enumerate(comp_genes)
        for (g,p) in gs
            mask = (genes .== g)
            gene_vec = setindex!(zeros(length(comp_genes) + 1), p, i)
            assignment_probs[:, mask] .= gene_vec .+ (1 - p) .* default_probs
            if (required_comps === nothing) || (i in required_comps)
                is_assigned[i] .|= mask
            end
        end
    end

    in_compartment_probs = Float64[]
    is_locked = any(hcat(is_assigned...), dims=2)[:]
    dists = zeros(size(pos_data, 2))

    for ia in is_assigned
        dists = max.(dists, maximum.(knn(KDTree(pos_data[:,ia]), pos_data, 2)[2]))
    end
    in_compartment_probs = [1.0 - max(min((d / scale - 1.) / 2., 1.), 0.) for d in dists]

    is_locked .&= (in_compartment_probs .> 0.99)
    is_locked .|= (in_compartment_probs .< 0.01)

    assignment_probs[end, :] .= 1 .- in_compartment_probs;
    assignment_probs[1:(end-1), :] .*= in_compartment_probs';

    return assignment_probs, is_locked
end

function segment_molecule_compartments(assignment_probs::Matrix{Float64}, is_locked::BitVector,
        adjacent_points::Vector{Vector{Int}}, adjacent_weights::Vector{Vector{Float64}}, confidence::Vector{Float64};
        n_iters_without_update=20, weights_pre_adjusted::Bool=false, weight_mult::Float64=1.0, tol::Float64=0.01,
        max_iter::Int=500, verbose::Bool=true)

    if !weights_pre_adjusted
        adjacent_weights = [weight_mult .* adjacent_weights[i] .* confidence[adjacent_points[i]] for i in 1:length(adjacent_weights)] # instead of multiplying each time in expect
    end

    max_diffs, change_fracs = Float64[], Float64[]
    assignment_probs_prev = deepcopy(assignment_probs)
    progress = verbose ? Progress(max_iter) : nothing

    # Run
    for i in 1:max_iter
        assignment_probs_prev .= assignment_probs
        # Iteration
        expect_molecule_compartments!(assignment_probs, adjacent_points, adjacent_weights, is_locked)

        # Stop criterion
        md, cf = estimate_difference_l0(assignment_probs, assignment_probs_prev, col_weights=confidence)
        push!(max_diffs, md)
        push!(change_fracs, cf)

        prog_vals = [("Iteration", i), ("Max. difference", md), ("Fraction of probs changed", cf)]
        if progress !== nothing
            next!(progress, showvalues=prog_vals)
        end
        if (i > n_iters_without_update) && (maximum(max_diffs[(end - n_iters_without_update):end]) < tol)
            if verbose
                finish!(progress, showvalues=prog_vals)
            end
            break
        end
    end

    # Wrap results
    assignment = vec(mapslices(x -> findmax(x)[2], assignment_probs, dims=1));
    return (assignment=assignment, diffs=max_diffs, assignment_probs=assignment_probs, change_fracs=change_fracs)
end

adjust_mrf_with_compartments(adjacent_points::Vector{Vector{Int}}, adjacent_weights::Vector{Vector{Float64}}, args...; kwargs...) =
    adjust_mrf_with_compartments!(deepcopy(adjacent_points), deepcopy(adjacent_weights), args...; kwargs...)

function adjust_mrf_with_compartments!(adjacent_points::Vector{Vector{Int}}, adjacent_weights::Vector{Vector{Float64}}, nuclei_probs::Vector{Float64}, cyto_probs::Vector{Float64}; min_nuc_prob::Float64=0.25, min_weight::Float64=0.01)
    for (m1,ids) in enumerate(adjacent_points)
        nuc_prob = nuclei_probs[m1]
        (min_nuc_prob > 0.25) || continue
        weights = adjacent_weights[m1]
        for (i2,m2) in enumerate(ids)
            weights[i2] *= 1 - cyto_probs[m2] * nuc_prob
        end

        mask = (weights .> min_weight)
        any(mask) || continue

        adjacent_points[m1] = adjacent_points[m1][mask]
        adjacent_weights[m1] = adjacent_weights[m1][mask]
    end
    return adjacent_points, adjacent_weights
end


text:
# isdebugging() = isdefined(Atom, :Debugger) && Debugger.isdebugging()

module Debugger
using ..Atom, MacroTools, ASTInterpreter2, Lazy, Hiccup

include("stepper.jl")

end # module


text:
using HPAT
using MPI

@acc @inline function cndf2( In::Array{Float64,1} )
    Out = 0.5 .+ 0.5 .* erf(0.707106781 .* In);
    return Out;
end

@acc hpat function blackscholes(iterations)

    sptprice   = Float64[ 42.0 for i=1:iterations ]
    strike     = Float64[ 40.0 + (i / iterations) for i=1:iterations ]
    rate       = Float64[ 0.5 for i=1:iterations ]
    volatility = Float64[ 0.2 for i=1:iterations ]
    time       = Float64[ 0.5 for i=1:iterations ]

    logterm = log10(sptprice ./ strike);
    powterm = .5 .* volatility .* volatility;
    den = volatility .* sqrt(time);
    d1 = (((rate .+ powterm) .* time) .+ logterm) ./ den;
    d2 = d1 .- den;
    NofXd1 = cndf2(d1);
    NofXd2 = cndf2(d2);
    futureValue = strike .* exp(- rate .* time);
    c1 = futureValue .* NofXd2
    Call = sptprice .* NofXd1 .- c1;
    Put  = Call .- futureValue .+ sptprice;
    sum(Put)
end

n = 10^8

if (length(ARGS) > 0)
	n = parse(Int, ARGS[1])
end

checksum = blackscholes(n)

if MPI.Comm_rank(MPI.COMM_WORLD)==0	println("checksum: ", checksum) end


text:
using LightXML
using PyPlot

lattice_file = ARGS[1]

lfile = parse_file(lattice_file)
graph = LightXML.root(lfile)
N = parse(Int, attribute(graph, "vertices"; required=true))
L = Int(sqrt(N))

figure()

# # vertices
# xmlvertices = get_elements_by_tagname(graph, "VERTEX")
# x = Array{Float64}(length(xmlvertices))
# y = Array{Float64}(length(xmlvertices))
# for (i,v) in enumerate(xmlvertices)
#   cords = map(x->parse(Float64,x), split(content(find_element(v, "COORDINATE")), " "))
#   x[i] = cords[1]
#   y[i] = cords[2]
# end
# scatter(x,y, c=1:length(x), cmap="Blues")

# edges
xmledges = get_elements_by_tagname(graph, "EDGE")
x = Array{Float64}(length(xmledges))
y = Array{Float64}(length(xmledges))
sql = reshape(1:N, L, L)
switch = true
types = Int[]
labels = String[]

for (i,edge) in enumerate(xmledges)
    attrs = attributes_dict(edge)
    src = parse(Int, attrs["source"])
    trg = parse(Int, attrs["target"])
    typ = parse(Int, attrs["type"])
    if !(typ in types) push!(types, typ); push!(labels, "Type $(length(types))"); end
    # if typ != 2
    #     continue
    # end

    v1 = [reverse(ind2sub(sql,src))...]
    v2 = [reverse(ind2sub(sql,trg))...]

    if src > trg # pbc bond
      continue # leave out pbs bonds
      if (v1[2] == L && v1[1] != L) || (v1[2] == L && switch)
        v2 = copy(v1)
        v2[2] += 1
        if v1 == [L, L]
          switch = false
        end
      elseif (v1[2] != L && v1[1] == L) || (v1[2] == L && !switch)
        v2 = copy(v1)
        v2[1] += 1
      end
    end
    plot(zip(v1-1,v2-1)..., "-", color="C$typ")
    # plot(zip(v1-1,v2-1)..., "-")
    plot(zip(v1-1,v2-1)..., "o", color="black")
end

#Create custom artists
handles = [plt[:Line2D]((0,1),(0,0), color="C$(i-1)") for i in 1:length(types)]

#Create legend from custom artist/label lists
legend(handles,labels,loc="upper center", bbox_to_anchor=(0.5, 1.12), ncol=4, fancybox=true)
xlabel("x")
ylabel("y")

axis("equal")
savefig("$(lattice_file)_edges.pdf")

close("all")

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
struct Zero{T}
    _val::T
    function Zero{T}() where {T}
        return new(zero(T))
    end
end

function Base.getindex(obj::Zero, idx::Int)
    return obj._val
end

function Base.getindex(obj::Zero, I::Vararg{Int, N}) where {N}
    return ntuple( i -> obj._val, N)
end

struct One{T}
    _val::T
    function One{T}() where {T}
        return new(one(T))
    end
end

function Base.getindex(obj::Zero, idx::Int)
    return obj._val
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
abstract type MCDMResult end
abstract type SCDMResult end
abstract type MCDMMethod end

"""
    struct MCDMSetting
        df::DataFrame
        weights::Array{Float64, 1}
        fns::Array{Function, 1}
    end

    Immutable data structure for a MCDM setting.

# Arguments
- `df::DataFrame`: The decision matrix in type of DataFrame.
- `weights::Array{Float64,1}`: Array of weights for each criterion.
- `fns::Array{Function, 1}`: Array of functions. The elements are either minimum or maximum.

# Description
Many methods including Topsis, Electre, Waspas, etc., use a decision matrix, weights, and directions
of optimizations in types of DataFrame, Vector, and Vector, respectively. The type MCDMSetting simply
holds these information to pass them into methods easly. Once a MCDMSetting object is created, the problem
can be passed into several methods like topsis(setting), electre(setting), waspas(setting), etc.

# Examples

```julia-repl
julia> df = DataFrame();
julia> df[:, :x] = Float64[9, 8, 7];
julia> df[:, :y] = Float64[7, 7, 8];
julia> df[:, :z] = Float64[6, 9, 6];
julia> df[:, :q] = Float64[7, 6, 6];

julia> w = Float64[4, 2, 6, 8];

julia> fns = makeminmax([maximum, maximum, maximum, maximum]);

julia> setting = MCDMSetting(df, w, fns)

julia> result = topsis(setting);
julia> # Same result can be obtained using
julia> result2 = mcdm(setting, TopsisMethod())
```

"""
struct MCDMSetting
    df::DataFrame
    weights::Array{Float64, 1}
    fns::Array{Function, 1}
end

struct TopsisResult <: MCDMResult
    decisionMatrix::DataFrame
    weights::Array{Float64,1}
    normalizedDecisionMatrix::DataFrame
    normalizedWeightedDecisionMatrix::DataFrame
    bestIndex::Int64
    scores::Array{Float64,1}
end

struct VikorResult <: MCDMResult
    decisionMatrix::DataFrame
    weights::Array{Float64,1}
    weightedDecisionMatrix::DataFrame
    bestIndex::Int64
    scores::Array{Float64,1}
end

struct ElectreResult <: MCDMResult
    decisionMatrix::DataFrame
    weights::Array{Float64,1}
    weightedDecisionMatrix::DataFrame
    fitnessTable::Array{Dict,1}
    nonfitnessTable::Array{Dict,1}
    fitnessMatrix::Array{Float64,2}
    nonfitnessMatrix::Array{Float64,2}
    C::Array{Float64,1}
    D::Array{Float64,1}
    bestIndex::Tuple
end

struct MooraResult <: MCDMResult
    mooraType::Symbol
    decisionMatrix::DataFrame
    weights::Array{Float64,1}
    weightedDecisionMatrix::DataFrame
    referenceMatrix::Union{DataFrame, Nothing}
    scores::Array{Float64,1}
    bestIndex::Int64
end

struct DematelResult <: MCDMResult
    comparisonMatrix::DataFrame
    NDMatrix::Array{Float64,2}
    relationShipMatrix::Array{Float64,2}
    c::Array{Float64,1}
    r::Array{Float64,1}
    c_plus_r::Array{Float64,1}
    c_minus_r::Array{Float64,1}
    threshold::Float64
    influenceMatrix::Array{Float64,2}
    weights::Array{Float64,1}
end

struct AHPConsistencyResult <: MCDMResult
    comparisonMatrix::DataFrame
    normalizedComparisonMatrix::DataFrame
    consistencyVector::Array{Float64,1}
    priority::Array{Float64,1}
    pc::Array{Float64,1}
    lambda_max::Float64
    CI::Float64
    RI::Float64
    CR::Float64
    isConsistent::Bool
end

struct AHPResult <: MCDMResult
    comparisonMatrixList::Array{DataFrame,1}
    criteriaComparisonMatrix::DataFrame
    criteriaConsistency::AHPConsistencyResult
    decisionMatrix::DataFrame
    scores::Array{Float64,1}
    weights::Array{Float64,1}
    bestIndex::Int64
end

struct NDSResult <: MCDMResult
    ranks::Array{Int64,1}
    bestIndex::Int64
end

struct LaplaceResult <: SCDMResult
    expected_values::Array{Float64,1}
    bestIndex::Int64
end

struct MaximinResult <: SCDMResult
    rowmins::Array{Float64,1}
    bestIndex::Int64
end

struct MaximaxResult <: SCDMResult
    rowmaxs::Array{Float64,1}
    bestIndex::Int64
end

struct MinimaxResult <: SCDMResult
    rowmaxs::Array{Float64,1}
    bestIndex::Int64
end

struct MiniminResult <: SCDMResult
    rowmins::Array{Float64,1}
    bestIndex::Int64
end

struct SavageResult <: SCDMResult
    regretMatrix::DataFrame
    scores::Array{Float64,1}
    bestIndex::Int64
end

struct HurwiczResult <: SCDMResult
    scores::Array{Float64,1}
    bestIndex::Int64
end

struct MLEResult <: SCDMResult
    scores::Array{Float64,1}
    bestIndex::Int64
end

struct ExpectedRegretResult <: SCDMResult
    scores::Array{Float64,1}
    bestIndex::Int64
end

struct GameResult <: MCDMResult
    row_player_probabilities::Array{Float64,1}
    value::Float64
end

struct DataEnvelopResult <: MCDMResult
    efficiencies::Array{Float64,1}
    references::DataFrame
    orderedcases::Array{Symbol,1}
end

struct GreyResult <: MCDMResult
    referenceRow::Array{Float64,1}
    normalizedMat::Matrix
    absoluteValueMat::Matrix
    greyTable::Matrix
    scores::Array{Float64,1}
    ordering::Array{Int64,1}
    bestIndex::Int64
end

struct SawResult <: MCDMResult
    decisionMatrix::DataFrame
    normalizedDecisionMatrix::DataFrame
    weights::Array{Float64,1}
    scores::Array{Float64,1}
    ranking::Array{Int64,1}
    bestIndex::Int64
end

struct ARASResult <: MCDMResult
    referenceRow::Array{Float64,1}
    extendMat::Array{Float64,2}
    normalizedMat::Array{Float64,2}
    optimalitydegrees::Array{Float64,1}
    scores::Array{Float64,1}
    orderings::Array{Int64,1}
    bestIndex::Int64
end

struct WPMResult <: MCDMResult
    decisionMatrix::DataFrame
    normalizedDecisionMatrix::DataFrame
    weights::Array{Float64,1}
    scores::Array{Float64,1}
    ranking::Array{Int64,1}
    bestIndex::Int64
end

struct WASPASResult <: MCDMResult
    decisionMatrix::DataFrame
    normalizedDecisionMatrix::DataFrame
    weights::Array{Float64,1}
    scores::Array{Float64,1}
    ranking::Array{Int64,1}
    bestIndex::Int64
end

struct EDASResult <: MCDMResult
    decisionMatrix::DataFrame
    weights::Array{Float64,1}
    scores::Array{Float64,1}
    ranking::Array{Int64,1}
    bestIndex::Int64
end

struct MARCOSResult <: MCDMResult
    decisionMatrix::DataFrame
    weights::Array{Float64,1}
    scores::Array{Float64,1}
    ranking::Array{Int64,1}
    bestIndex::Int64
end

struct MABACResult <: MCDMResult
    decisionMatrix::DataFrame
    weights::Array{Float64,1}
    scores::Array{Float64,1}
    ranking::Array{Int64,1}
    bestIndex::Int64
end

struct MAIRCAResult <: MCDMResult
    decisionMatrix::DataFrame
    weights::Array{Float64,1}
    scores::Array{Float64,1}
    ranking::Array{Int64,1}
    bestIndex::Int64
end

struct COPRASResult <: MCDMResult
    decisionMatrix::DataFrame
    weights::Array{Float64,1}
    scores::Array{Float64,1}
    ranking::Array{Int64,1}
    bestIndex::Int64
end

struct PrometheeResult <: MCDMResult
    decisionMatrix::DataFrame
    weights::Array{Float64,1}
    scores::Array{Float64,1}
    ranking::Array{Int64,1}
    bestIndex::Int64
end

struct CoCoSoResult <: MCDMResult
    decisionMatrix::DataFrame
    weights::Array{Float64,1}
    scores::Array{Float64,1}
    ranking::Array{Int64,1}
    bestIndex::Int64
end

struct CRITICResult <: MCDMResult
    decisionMatrix::DataFrame
    w::Array{Float64,1}
    ranking::Array{Int64,1}
    bestIndex::Int64
end

struct EntropyResult <: MCDMResult
    decisionMatrix::DataFrame
    w::Array{Float64,1}
    ranking::Array{Int64,1}
    bestIndex::Int64
end

struct CODASResult <: MCDMResult
    decisionMatrix::DataFrame
    w::Array{Float64,1}
    scores::Array{Float64,1}
    ranking::Array{Int64,1}
    bestIndex::Int64
end

struct TopsisMethod <: MCDMMethod end

struct ElectreMethod <: MCDMMethod
end

struct ArasMethod <: MCDMMethod
end

struct CocosoMethod <: MCDMMethod
    lambda::Float64
end

CocosoMethod()::CocosoMethod = CocosoMethod(0.5)

struct CodasMethod <: MCDMMethod
    tau::Float64
end

CodasMethod()::CodasMethod = CodasMethod(0.02)

struct CoprasMethod <: MCDMMethod
end

struct CriticMethod <: MCDMMethod
end

struct EdasMethod <: MCDMMethod
end

struct GreyMethod <: MCDMMethod
    zeta::Float64
end

GreyMethod() :: GreyMethod = GreyMethod(0.5)

struct MabacMethod <: MCDMMethod
end

struct MaircaMethod <: MCDMMethod
end

struct MooraMethod <: MCDMMethod
    method::Symbol
end

MooraMethod() :: MooraMethod = MooraMethod(:reference)

struct PrometheeMethod <: MCDMMethod
    pref::Array{Function, 1}
    qs::Array{Union{Nothing, Float64}, 1}
    ps::Array{Union{Nothing, Float64}, 1}
end

PrometheeMethod(pref::Array{Function, 1}, qs::Array{Float64,1}, ps::Array{Float64, 1}) :: PrometheeMethod = PrometheeMethod(pref, qs, ps)

struct SawMethod <: MCDMMethod
end

struct VikorMethod <: MCDMMethod
    v::Float64
end

VikorMethod()::VikorMethod = VikorMethod(0.5)

struct WaspasMethod <: MCDMMethod
    lambda::Float64
end

WaspasMethod() :: WaspasMethod = WaspasMethod(0.5)

struct WPMMethod <: MCDMMethod
end

struct MarcosMethod <: MCDMMethod
end

max_stars_repo_name: jbytecode/jMCDM

max_issues_repo_name: jbytecode/jMCDM

max_forks_repo_name: jbytecode/jMCDM


text:
using GalerkinSparseGrids
using Test
using HCubature
using StaticArrays
using SparseArrays

#---------------------------------------
# Testing vector hier DG reconstruction
#---------------------------------------

@info "Testing full DG reconstruction 1-D with vector coefficients... "

@testset "vhier_DC.jl" begin
    @info "Testing full DG reconstruction with vector coefficients in 1D"
    @testset "Full DG reconstruction 1D with vector coefficients" begin
        for k in 1:5
            for l in 1:6
                vect = vcoeffs_DG(1, k, l, x->sin(4*x[1]); scheme="full")
                dict = V2D(1, k, l, vect; scheme="full")
                err = x->(reconstruct_DG(dict,[x[1]])-sin(4*x[1]))^2
                @test hquadrature(err,0,1; rtol=1.0e-9, atol=1.0e-12)[1] < 1/(1<<((l+k-1)))
            end
        end
    end

    # Testing D2V and V2D

    # Full case:

    @info "Testing D2V and V2D Full Case 2D... "
    @testset "D2V and V2D full case 2D" begin
        for k in 1:4
            for l in 1:3
                dict = coeffs_DG(2, k, l, x->sin(4*x[1]+x[2]); scheme="full")
                vect = D2V(2, k, l, dict; scheme="full")
                @test V2D(2, k, l, vect; scheme="full") == dict
            end
        end

        for k in 1:4
            for l in 1:3
                vect = vcoeffs_DG(2, k, l, x->sin(4*x[1]+x[2]); scheme="full")
                dict = V2D(2, k, l, vect; scheme="full")
                @test D2V(2, k, l, dict; scheme="full") == vect
            end
        end
    end

    # Sparse case:

    @info "Testing D2V and V2D Sparse Case 2D... "
    @testset "D2V and V2D sparse case 2D" begin
        for k in 1:4
            for l in 1:3
                dict = coeffs_DG(2, k, l, x->sin(4*x[1]+x[2]); scheme="sparse")
                vect = D2V(2, k, l, dict; scheme="sparse")
                @test V2D(2, k, l, vect; scheme="sparse")==dict
            end
        end

        for k in 1:4
            for l in 1:3
                vect = vcoeffs_DG(2, k, l, x->sin(4*x[1]+x[2]); scheme="sparse")
                dict= V2D(2, k, l, vect; scheme="sparse")
                @test D2V(2, k, l, dict; scheme="sparse")==vect
            end
        end
    end
end


text:
function create_input_data(X, y, labels, nntr::Int, nnval::Int, natr::Int, naval::Int, ntst::Int;
		seed=nothing)
    # create anomalous and normal splits
    nX = X[:,:,:,labels]
    aX = X[:,:,:,.!labels];
    labels = Int.(labels);

    # fix the seed
    isnothing(seed) ? nothing : Random.seed!(seed)

    # sample normal data
    nn = size(nX,4)
    randinds = sample(1:nn, nn, replace=false);
    tr_nX = nX[:,:,:,randinds[1:nntr]]
    val_nX = nX[:,:,:,randinds[nntr+1:nntr+nnval]]
    tst_nX = nX[:,:,:,randinds[nntr+nnval+1:end]];

    # sample anomalous data
    an = size(aX, 4)
    randinds = sample(1:an, an, replace=false);
    tr_aX = aX[:,:,:,randinds[1:natr]]
    val_aX = aX[:,:,:,randinds[natr+1:natr+naval]]
    tst_aX = aX[:,:,:,randinds[natr+naval+1:end]];

    # mixed data
    tr_X = cat(tr_nX, tr_aX, dims=4)
    tr_y = cat(zeros(nntr), ones(natr), dims=1)
    val_X = cat(val_nX, val_aX, dims=4)
    val_y = cat(zeros(nnval), ones(naval), dims=1)
    tst_X = cat(sample_tensor(tst_nX, ntst), sample_tensor(tst_aX, ntst), dims=4)
    tst_y = cat(zeros(ntst), ones(ntst), dims=1)

    # restart seed
    Random.seed!()

    return tr_nX, (tr_X, tr_y), (val_X, val_y), (tst_X, tst_y)
end

function construct_conv_classifier(ks, activation, strd, ncs, datasize)
	af = eval(Meta.parse(activation))
    inncs = vcat([datasize[3]], ncs[1:end-1])
    classifier_conv_part = Chain([Conv(k, inc => outc, af) for (k, inc, outc) in zip(ks, inncs, ncs)]...)
    odims = outdims(classifier_conv_part, datasize)
    odim = reduce(*, odims[1:3])
    gpu(Chain(classifier_conv_part..., x->reshape(x, odim, size(x,4)),  Dense(odim, 2)))
end

function train_conv_classifier!(classifier, tr_X, tr_y, val_X, val_y;
        lr = 0.001f0,
        λ = 0f0,
        batchsize = 64,
        nepochs = 1000,
        patience = 50,
        verb = true
    )

    tr_classifier = deepcopy(classifier)
    opt = ADAM(lr)
    ps = params(tr_classifier)
    loss(x,y) = Flux.Losses.logitcrossentropy(tr_classifier(gpu(x)), gpu(y)) + λ*sum(sqnorm, ps)
    loss(x) = loss(x...)
    oh_tr_y = Flux.onehotbatch(Bool.(tr_y), 0:1)
    oh_val_y = Flux.onehotbatch(Bool.(val_y), 0:1)

    local data_itr
    @suppress begin
        data_itr = Flux.Data.DataLoader((tr_X, oh_tr_y), batchsize=batchsize)
    end

    _val_auc = 0
    _patience = 1
    history = MVHistory()
    @info "Starting convolutional classifier training..."
    for epoch in 1:nepochs
        for (x, y) in data_itr
            param_update!(loss, ps, (x, y), opt)
        end
        local tr_l, val_l, tr_auc, val_auc
        @suppress begin
			tr_l = mean(map(x->cpu(loss(x)), Flux.Data.DataLoader((tr_X, oh_tr_y), batchsize=batchsize)))
	        val_l = mean(map(x->cpu(loss(x)), Flux.Data.DataLoader((val_X, oh_val_y), batchsize=batchsize)))
            tr_auc = auc_val(tr_y, cpu(classifier_score(tr_classifier, tr_X, batchsize)))
	        val_auc = auc_val(val_y, cpu(classifier_score(tr_classifier, val_X, batchsize)))
	    end
        if verb
            println("Epoch $epoch: loss = $(tr_l), tr AUC = $(tr_auc), val AUC = $(val_auc)")
        end
        push!(history, :tr_loss, epoch, tr_l)
        push!(history, :val_loss, epoch, val_l)
        push!(history, :tr_auc, epoch, tr_auc)
        push!(history, :val_auc, epoch, val_auc)

        if (val_auc <= _val_auc)
            if _patience >= patience
                if verb
                    @info "Stopping early after $epoch epochs."
                end
                return classifier, history, opt
            else
                _patience += 1
            end
        else
        	classifier = deepcopy(tr_classifier)
            _val_auc = val_auc
        end
    end
    @info "Finished."

    return classifier, history, opt
end

max_stars_repo_name: vitskvara/HierarchicalAD

max_issues_repo_name: vitskvara/HierarchicalAD

max_forks_repo_name: vitskvara/HierarchicalAD


text:
"""
    rsgd(;input::AbstractString="", outdir::Union{Nothing,AbstractString}=nothing, scale::AbstractString="ftt", pseudocount::Number=1f0, rowmeanlist::AbstractString="", rowvarlist::AbstractString="",colsumlist::AbstractString="", dim::Number=3, stepsize::Number=0.1, numbatch::Number=100, numepoch::Number=3, scheduling::AbstractString="robbins-monro", g::Number=0.9, epsilon::Number=1.0e-8, lower::Number=0, upper::Number=1.0f+38, expvar::Number=0.1f0, evalfreq::Number=5000, offsetFull::Number=1f-20, offsetStoch::Number=1f-6, initW::Union{Nothing,AbstractString}=nothing, initV::Union{Nothing,AbstractString}=nothing, logdir::Union{Nothing,AbstractString}=nothing, perm::Bool=false, cper::Number=1f0)

Riemannian stochastic gradient descent method.

Input Arguments
---------
- `input` : Julia Binary file generated by `OnlinePCA.csv2bin` function.
- `outdir` : The directory specified the directory you want to save the result.
- `scale` : {log,ftt,raw}-scaling of the value.
- `pseudocount` : The number specified to avoid NaN by log10(0) and used when `Feature_LogMeans.csv` <log10(mean+pseudocount) value of each feature> is generated.
- `rowmeanlist` : The mean of each row of matrix. The CSV file is generated by `OnlinePCA.sumr` functions.
- `rowvarlist` : The variance of each row of matrix. The CSV file is generated by `OnlinePCA.sumr` functions.
- `colsumlist` : The sum of counts of each columns of matrix. The CSV file is generated by `OnlinePCA.sumr` functions.
- `dim` : The number of dimension of PCA.
- `stepsize` : The parameter used in every iteration.
- `numbatch` : The number of batch size.
- `numepoch` : The number of epoch.
- `scheduling` : Learning parameter scheduling. `robbins-monro`, `momentum`, `nag`, and `adagrad` are available.
- `g` : The parameter that is used when scheduling is specified as nag.
- `epsilon` : The parameter that is used when scheduling is specified as adagrad.
- `lower` : Stopping Criteria (When the relative change of error is below this value, the calculation is terminated)
- `upper` : Stopping Criteria (When the relative change of error is above this value, the calculation is terminated)
- `evalfreq` : Evaluation Frequency of Reconstruction Error
- `offsetFull` : Off set value for avoding overflow when calculating full gradient
- `offsetStoch` : Off set value for avoding overflow when calculating stochastic gradient
- `initW` : The CSV file saving the initial values of eigenvectors.
- `initV` : The CSV file saving the initial values of loadings.
- `logdir` : The directory where intermediate files are saved, in every evalfreq (e.g. 5000) iteration.
- `perm` : Whether the data matrix is shuffled at random.
- `cper` : Count per X (e.g. CPM: Count per million <1e+6>)

Output Arguments
---------
- `W` : Eigen vectors of covariance matrix (No. columns of the data matrix × dim)
- `λ` : Eigen values (dim × dim)
- `V` : Loading vectors of covariance matrix (No. rows of the data matrix × dim)
- `Scores` : Principal component scores
- `ExpVar` : Explained variance by the eigenvectors
- `TotalVar` : Total variance of the data matrix
- stop : Whether the calculation is converged

Reference
---------
- RSGD-PCA : [Silvere Bonnabel, 2013](https://arxiv.org/abs/1111.5280)
"""
function rsgd(;input::AbstractString="", outdir::Union{Nothing,AbstractString}=nothing, scale::AbstractString="ftt", pseudocount::Number=1f0, rowmeanlist::AbstractString="", rowvarlist::AbstractString="",colsumlist::AbstractString="", dim::Number=3, stepsize::Number=0.1, numbatch::Number=100, numepoch::Number=3, scheduling::AbstractString="robbins-monro", g::Number=0.9, epsilon::Number=1.0e-8, lower::Number=0, upper::Number=1.0f+38, expvar::Number=0.1f0, evalfreq::Number=5000, offsetFull::Number=1f-20, offsetStoch::Number=1f-6, initW::Union{Nothing,AbstractString}=nothing, initV::Union{Nothing,AbstractString}=nothing, logdir::Union{Nothing,AbstractString}=nothing, perm::Bool=false, cper::Number=1f0)
    # Initial Setting
    pca = RSGD()
    if scheduling == "robbins-monro"
        scheduling = ROBBINS_MONRO()
    elseif scheduling == "momentum"
        scheduling = MOMENTUM()
    elseif scheduling == "nag"
        scheduling = NAG()
    elseif scheduling == "adagrad"
        scheduling = ADAGRAD()
    else
        error("Specify the scheduling as robbins-monro, momentum, nag or adagrad")
    end
    pseudocount, stepsize, g, epsilon, W, v, D, rowmeanvec, rowvarvec, colsumvec, N, M, TotalVar, lower, upper, evalfreq, offsetFull, offsetStoch = init(input, pseudocount, stepsize, g, epsilon, dim, rowmeanlist, rowvarlist, colsumlist, initW, initV, logdir, pca, lower, upper, evalfreq, offsetFull, offsetStoch, cper, scale)
    # Perform PCA
    out = rsgd(input, outdir, scale, pseudocount, rowmeanlist, rowvarlist, colsumlist, dim, stepsize, numbatch, numepoch, scheduling, g, epsilon, logdir, pca, W, v, D, rowmeanvec, rowvarvec, colsumvec, N, M, TotalVar, lower, upper, evalfreq, offsetStoch, perm, cper)
    if outdir isa String
        output(outdir, out, expvar)
    end
    return out
end

function rsgd(input, outdir, scale, pseudocount, rowmeanlist, rowvarlist, colsumlist, dim, stepsize, numbatch, numepoch, scheduling, g, epsilon, logdir, pca, W, v, D, rowmeanvec, rowvarvec, colsumvec, N, M, TotalVar, lower, upper, evalfreq, offsetStoch, perm, cper)
    N, M = nm(input)
    tmpN = zeros(UInt32, 1)
    tmpM = zeros(UInt32, 1)
    x = zeros(UInt32, M)
    normx = zeros(Float32, M)
    batchgrad = zeros(Float32, M, dim)
    # If not 0 the calculation is converged
    stop = 0
    s = 1
    n = 1
    batchidx = 1
    batchcount = 1
    # Each epoch s
    progress = Progress(numepoch*N)
    while(stop == 0 && s <= numepoch)
        open(input) do file
            stream = ZstdDecompressorStream(file)
            read!(stream, tmpN)
            read!(stream, tmpM)
            # Each step n
            while(stop == 0 && n <= N)
                next!(progress)
                # Row vector of data matrix
                read!(stream, x)
                normx = normalizex(x, n, stream, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, cper)
                if perm
                    normx .= normx[randperm(length(normx))]
                end
                # Update Eigen vector
                W, v, batchgrad, batchidx, batchcount = rsgdupdate(scheduling, stepsize, g, epsilon, D, N, M, dim, W, v, batchgrad, batchidx, batchcount, numbatch, normx, s, n, offsetStoch)
                # NaN
                checkNaN(N, s, n, W, evalfreq, pca)
                # Retraction
                W .= Array(qr!(W).Q)
                # save log file
                if logdir isa String
                    stop = outputlog(N, s, n, input, dim, logdir, W, pca, TotalVar, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, lower, upper, stop, evalfreq, cper)
                end
                n += 1
            end
            close(stream)
        end
        # save log file
        if logdir isa String
            stop = outputlog(s, input, dim, logdir, W, GD(), TotalVar, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, lower, upper, stop, cper)
        end
        s += 1
        if n == N + 1
            n = 1
        end
    end

    # Return, W, λ, V
    out = WλV(W, input, dim, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, TotalVar, cper)
    return (out[1], out[2], out[3], out[4], out[5], out[6], stop)
end

# RSGD × Robbins-Monro
function rsgdupdate(scheduling::ROBBINS_MONRO, stepsize, g, epsilon, D, N, M, dim, W, v, batchgrad, batchidx, batchcount, numbatch, normx, s, n, offsetStoch)
    if mod(batchidx, numbatch) == 0 && n != 1
        W .= W .+ Pw(batchgrad, W)
        batchgrad = zeros(Float32, M, dim)
        batchidx = 1
        batchcount += 1
    else
        batchgrad .+= ∇fn(W, normx, D , M, stepsize/batchcount, offsetStoch)
        batchidx += 1
    end
    v = nothing
    return W, v, batchgrad, batchidx, batchcount
end

# RSGD × Momentum
function rsgdupdate(scheduling::MOMENTUM, stepsize, g, epsilon, D, N, M, dim, W, v, batchgrad, batchidx, batchcount, numbatch, normx, s, n, offsetStoch)
    if mod(batchidx, numbatch) == 0 && n != 1
        v .= g .* v .+ Pw(batchgrad, W)
        W .= W .+ v
        batchgrad = zeros(Float32, M, dim)
        batchidx = 1
        batchcount += 1
    else
        batchgrad .+= ∇fn(W, normx, D, M, stepsize, offsetStoch)
        batchidx += 1
    end
    return W, v, batchgrad, batchidx, batchcount
end

# RSGD × NAG
function rsgdupdate(scheduling::NAG, stepsize, g, epsilon, D, N, M, dim, W, v, batchgrad, batchidx, batchcount, numbatch, normx, s, n, offsetStoch)
    if mod(batchidx, numbatch) == 0 && n != 1
        v = g .* v .+ Pw(batchgrad, W)
        W .= W .+ v
        batchgrad = zeros(Float32, M, dim)
        batchidx = 1
        batchcount += 1
    else
        batchgrad .+= ∇fn(W - g .* v, normx, D, M, stepsize, offsetStoch)
        batchidx += 1
    end
    return W, v, batchgrad, batchidx, batchcount
end

# RSGD × Adagrad
function rsgdupdate(scheduling::ADAGRAD, stepsize, g, epsilon, D, N, M, dim, W, v, batchgrad, batchidx, batchcount, numbatch, normx, s, n, offsetStoch)
    if mod(batchidx, numbatch) == 0 && n != 1
        grad = Pw(batchgrad, W)
        grad = grad / stepsize
        v .= v .+ grad .* grad
        W .= W .+ stepsize ./ (sqrt.(v) .+ epsilon) .* grad
        batchgrad = zeros(Float32, M, dim)
        batchidx = 1
        batchcount += 1
    else
        batchgrad .+= ∇fn(W, normx, D, M, stepsize, offsetStoch)
        batchidx += 1
    end
    return W, v, batchgrad, batchidx, batchcount
end


text:
module TestLale
using Test

# suppress warnings
@info "suppressing PyCall warnings"
using PyCall
warnings = pyimport("warnings")
warnings.filterwarnings("ignore")

# test modules
include("test_sklearn.jl")
include("test_autogen.jl")
include("test_lalepreprocessing.jl")
include("test_laleoptimizer_sk.jl")
include("test_laleoptimizer_ag.jl")
include("test_imblearn.jl")
include("test_project.jl")

end


text:
@test e1.first == src(e1) == 1
@test e1.second == dst(e1) == 2
@test reverse(e1) == re1

e2 = Edge(1,3)
e3 = Edge(1,4)
e4 = Edge(2,5)
e5 = Edge(3,5)

@test LightGraphs.is_ordered(e5)
@test !LightGraphs.is_ordered(reverse(e5))

g = Graph(5)
@test add_edge!(g, 1, 2)
@test add_edge!(g, e2)
@test add_edge!(g, e3)
@test add_edge!(g, e4)
@test add_edge!(g, e5)

h = DiGraph(5)
@test add_edge!(h, 1, 2)
@test add_edge!(h, e2)
@test add_edge!(h, e3)
@test add_edge!(h, e4)
@test add_edge!(h, e5)

@test vertices(g) == 1:5
i = 0
for e in edges(g)
    i+=1
end
@test i == 5
@test has_edge(g,3,5)
# @test edges(g) == Set([e1, e2, e3, e4, e5])
# @test Set{Edge}(edges(g)) == Set([e1, e2, e3, e4, e5])
# fadj, badj, and adj tested in graphdigraph.jl

@test degree(g) == [3, 2, 2, 1, 2]
@test indegree(g) == [3, 2, 2, 1, 2]
@test indegree(g,1) == 3
@test outdegree(g) == [3, 2, 2, 1, 2]
@test outdegree(g,1) == 3
@test degree(h) == [3, 2, 2, 1, 2]
@test indegree(h) == [0, 1, 1, 1, 2]
@test indegree(h,1) == 0
@test outdegree(h) == [3, 1, 1, 0, 0]
@test outdegree(h,1) == 3
@test in_neighbors(h,5) == LightGraphs.badj(h)[5] == LightGraphs.badj(h,5) == [2, 3]
@test out_neighbors(h,1) == LightGraphs.fadj(h)[1] == LightGraphs.fadj(h,1) == [2, 3, 4]

@test p1 == g2
@test issubset(h2, h1)

@test has_edge(g, 1, 2)
@test in_edges(g, 2) == [e1, reverse(e4)]
@test out_edges(g, 1) == [e1, e2, e3]

@test add_vertex!(g) && nv(g) == 6
@test add_vertices!(g,5) && nv(g) == 11
@test has_vertex(g, 11)
@test ne(g) == 5
@test !is_directed(g)
@test is_directed(h)

@test δ(g) == δin(g) == δout(g) == 0
@test Δ(g) == Δout(g) == 3
@test Δin(h) == 2
@test δ(h) == 1
@test δin(h) == 0
@test δout(h) == 0
@test CompleteGraph(4) == CompleteGraph(4)
@test CompleteGraph(4) != PathGraph(4)
@test CompleteDiGraph(4) != PathDiGraph(4)
@test CompleteDiGraph(4) == CompleteDiGraph(4)

@test degree_histogram(CompleteDiGraph(10)).weights == [10]
@test degree_histogram(CompleteGraph(10)).weights == [10]

@test neighbors(g, 1) == [2, 3, 4]
@test common_neighbors(g, 2, 3) == [1, 5]
@test common_neighbors(h, 2, 3) == [5]

@test add_edge!(g, 1, 1)
@test !add_edge!(g, 1, 1)
@test rem_edge!(g, 1, 1)
@test !rem_edge!(g, 1, 1)
@test ne(g) == 5
@test rem_edge!(g, 1, 2)
@test ne(g) == 4
@test !rem_edge!(g, 2, 1)
add_edge!(g, 1, 2)
@test ne(g) == 5

@test has_edge(g,2,1)
@test has_edge(g,1,2)
@test rem_edge!(g, 2, 1)
@test add_edge!(h, 1, 1)
@test rem_edge!(h, 1, 1)
@test rem_edge!(h, 1, 2)
@test !rem_edge!(h, 1, 2)

function test_rem_edge(g, srcv)
    srcv = 2
    for dstv in collect(neighbors(g, srcv))
        rem_edge!(g, srcv, dstv)
    end
    @test length(neighbors(g,srcv)) == 0
end
for v in vertices(g)
    test_rem_edge(copy(g),v)
end
for v in vertices(h)
    test_rem_edge(copy(h),v)
end

@test g == copy(g)
@test !(g === copy(g))
g10 = CompleteGraph(5)
@test rem_vertex!(g10, 1)
@test g10 == CompleteGraph(4)
@test rem_vertex!(g10, 4)
@test g10 == CompleteGraph(3)
@test !rem_vertex!(g10, 9)

g10 = CompleteDiGraph(5)
@test rem_vertex!(g10, 1)
@test g10 == CompleteDiGraph(4)
rem_vertex!(g10, 4)
@test g10 == CompleteDiGraph(3)
@test !rem_vertex!(g10, 9)
g10 = PathGraph(5)
@test rem_vertex!(g10, 5)
@test g10 == PathGraph(4)
@test rem_vertex!(g10, 4)
@test g10 == PathGraph(3)

g10 = PathDiGraph(5)
@test rem_vertex!(g10, 5)
@test g10 == PathDiGraph(4)
@test rem_vertex!(g10, 4)
@test g10 == PathDiGraph(3)

g10 = PathDiGraph(5)
@test rem_vertex!(g10, 1)
h10 = PathDiGraph(6)
@test rem_vertex!(h10, 1)
@test rem_vertex!(h10, 1)
@test g10 == h10

g10 = CycleGraph(5)
@test rem_vertex!(g10, 5)
@test g10 == PathGraph(4)

g10 = PathGraph(3)
@test rem_vertex!(g10, 2)
@test g10 == Graph(2)

g10 = PathGraph(4)
@test rem_vertex!(g10, 3)
h10 =Graph(3)
@test add_edge!(h10,1,2)
@test g10 == h10

g10 = CompleteGraph(5)
@test rem_vertex!(g10, 3)
@test g10 == CompleteGraph(4)


text:
#######################################################################################################################################################################################################
#
# Changes to this structure
# General
#     2022-Jan-14: refactor the Leaf structure within BIO, PRC, PSM as fields
#     2022-Jan-24: add p_CO₂_s to the structure
#     2022-Jan-24: fix documentation
#     2022-Feb-07: moved FLM to PRC
# Bug fixes:
#     2022-Jan-24: add FT control to p_CO₂_i
# To do
#     TODO: add leaf hydraulics as a field as well
#     TODO: link leaf water content to BIO_PHYSICS.l_H₂O
#
#######################################################################################################################################################################################################
"""

$(TYPEDEF)

Structure to save leaf parameters

# Fields

$(TYPEDFIELDS)

"""
mutable struct Leaf{FT<:AbstractFloat}
    # parameters that do not change with time
    "[`LeafBiophysics`](@ref) type leaf biophysical parameters"
    BIO::LeafBiophysics{FT}
    "[`AbstractReactionCenter`](@ref) type photosynthesis reaction center"
    PRC::AbstractReactionCenter{FT}
    "[`AbstractPhotosynthesisModel`](@ref) type photosynthesis model"
    PSM::AbstractPhotosynthesisModel{FT}

    # prognostic variables that change with time
    "Absorbed photosynthetically active radiation `[μmol m⁻² s⁻¹]`"
    apar::FT
    "Stomatal conductance to water vapor `[mol m⁻² s⁻¹]`"
    g_H₂O_s::FT
    "Current leaf temperature"
    t::FT

    # dignostic variables that change with time
    "Total leaf diffusive conductance to CO₂ `[mol m⁻² s⁻¹]`"
    g_CO₂::FT
    "Boundary leaf diffusive conductance to CO₂ `[mol m⁻² s⁻¹]`"
    g_CO₂_b::FT
    "Leaf internal CO₂ partial pressure `[Pa]`"
    p_CO₂_i::FT
    "Leaf surface CO₂ partial pressure `[Pa]`"
    p_CO₂_s::FT
    "Saturation H₂O vapor pressure, need to update with temperature and leaf water pressure `[Pa]`"
    p_H₂O_sat::FT

    # caches to speed up calculations
    "Last leaf temperature. If different from t, then make temperature correction"
    _t::FT
end

#######################################################################################################################################################################################################
#
# Changes to this constructor
# General
#     2022-Jan-14: add C3 and C4 constructors
#     2022-Jan-24: add C3Cytochrome constructor
#     2022-Jan-24: add p_CO₂_s to the constructor
#     2022-Jan-24: add documentation
#     2022-Feb-07: remove fluorescence model from Leaf struct
#     2022-Feb-11: set default APAR = 1000
#     2022-Feb-11: add colimit option in constructor to enable quick deployment of quadratic colimitation
#
#######################################################################################################################################################################################################
"""

    Leaf{FT}(psm::String, wls::WaveLengthSet{FT} = WaveLengthSet{FT}(); colimit::Bool = false) where {FT<:AbstractFloat}

Constructor for `Leaf`, given
- `psm` Photosynthesis model type, must be `C3`, `C3Cytochrome`, or `C4`
- `wls` [`WaveLengthSet`](@ref) type structure that determines the dimensions of leaf parameters

---
# Examples
```julia
leaf_c3 = Leaf{Float64}("C3");
leaf_c4 = Leaf{Float64}("C4");
leaf_cy = Leaf{Float64}("C3Cytochrome");
leaf_c3 = Leaf{Float64}("C3"; colimit = true);
leaf_c4 = Leaf{Float64}("C4"; colimit = true);
leaf_cy = Leaf{Float64}("C3Cytochrome"; colimit = true);
wls = WaveLengthSet{FT}(collect(400:10:2500));
leaf_c3 = Leaf{Float64}("C3", wls);
leaf_c4 = Leaf{Float64}("C4", wls);
leaf_cy = Leaf{Float64}("C3Cytochrome", wls);
```
"""
Leaf{FT}(psm::String, wls::WaveLengthSet{FT} = WaveLengthSet{FT}(); colimit::Bool = false) where {FT<:AbstractFloat} = (
    @assert psm in ["C3", "C3Cytochrome", "C4"] "Photosynthesis model ID must be C3, C4, or C3Cytochrome!";

    if psm == "C3"
        _prc = VJPReactionCenter{FT}();
        _psm = C3VJPModel{FT}(colimit = colimit);
    elseif psm == "C3Cytochrome"
        _prc = CytochromeReactionCenter{FT}();
        _psm = C3CytochromeModel{FT}(colimit = colimit);
    elseif psm == "C4"
        _prc = VJPReactionCenter{FT}();
        _psm = C4VJPModel{FT}(colimit = colimit);
    end;

    return Leaf{FT}(
                LeafBiophysics{FT}(wls),            # BIO
                _prc,                               # PRC
                _psm,                               # PSM
                1000,                               # apar
                0.01,                               # g_H₂O_s
                T_25(),                             # t
                0.01,                               # g_CO₂
                3.0,                                # g_CO₂_b
                20,                                 # p_CO₂_i
                40,                                 # p_CO₂_s
                saturation_vapor_pressure(T_25()),  # p_H₂O_sat
                0)                                  # _t
);


text:
## R code 15.3
dlist <- list(
    D_obs = standardize( d$Divorce ),
    D_sd = d$Divorce.SE / sd( d$Divorce ),
    M = standardize( d$Marriage ),
    A = standardize( d$MedianAgeMarriage ),
    N = nrow(d)
)

m15.1 <- ulam(
    alist(
        D_obs ~ dnorm( D_true , D_sd ),
        vector[N]:D_true ~ dnorm( mu , sigma ),
        mu <- a + bA*A + bM*M,
        a ~ dnorm(0,0.2),
        bA ~ dnorm(0,0.5),
        bM ~ dnorm(0,0.5),
        sigma ~ dexp(1)
    ) , data=dlist , chains=4 , cores=4 )

## R code 15.4
precis( m15.1 , depth=2 )

## R code 15.5
dlist <- list(
    D_obs = standardize( d$Divorce ),
    D_sd = d$Divorce.SE / sd( d$Divorce ),
    M_obs = standardize( d$Marriage ),
    M_sd = d$Marriage.SE / sd( d$Marriage ),
    A = standardize( d$MedianAgeMarriage ),
    N = nrow(d)
)

m15.2 <- ulam(
    alist(
        D_obs ~ dnorm( D_est , D_sd ),
        vector[N]:D_est ~ dnorm( mu , sigma ),
        mu <- a + bA*A + bM*M_est[i],
        M_obs ~ dnorm( M_est , M_sd ),
        vector[N]:M_est ~ dnorm( 0 , 1 ),
        a ~ dnorm(0,0.2),
        bA ~ dnorm(0,0.5),
        bM ~ dnorm(0,0.5),
        sigma ~ dexp( 1 )
    ) , data=dlist , chains=4 , cores=4 )

## R code 15.7
N <- 500
A <- rnorm(N)
M <- rnorm(N,-A)
D <- rnorm(N,A)
A_obs <- rnorm(N,A)

## R code 15.8
N <- 100
S <- rnorm( N )
H <- rbinom( N , size=10 , inv_logit(S) )

## R code 15.9
D <- rbern( N ) # dogs completely random
Hm <- H
Hm[D==1] <- NA

## R code 15.10
D <- ifelse( S > 0 , 1 , 0 )
Hm <- H
Hm[D==1] <- NA

## R code 15.11
set.seed(501)
N <- 1000
X <- rnorm(N)
S <- rnorm(N)
H <- rbinom( N , size=10 , inv_logit( 2 + S - 2*X ) )
D <- ifelse( X > 1 , 1 , 0 )
Hm <- H
Hm[D==1] <- NA

## R code 15.12
dat_list <- list(
    H = H,
    S = S )

m15.3 <- ulam(
    alist(
        H ~ binomial( 10 , p ),
        logit(p) <- a + bS*S,
        a ~ normal( 0 , 1 ),
        bS ~ normal( 0 , 0.5 )
    ), data=dat_list , chains=4 )

## R code 15.13
dat_list0 <- list(
    H = H[D==0],
    S = S[D==0] )

m15.4 <- ulam(
    alist(
        H ~ binomial( 10 , p ),
        logit(p) <- a + bS*S,
        a ~ normal( 0 , 1 ),
        bS ~ normal( 0 , 0.5 )
    ), data=dat_list0 , chains=4 )

## R code 15.17
dat_list <- list(
    K = standardize( d$kcal.per.g ),
    B = standardize( d$neocortex.prop ),
    M = standardize( d$logmass )
)

m15.3 <- ulam(
    alist(
        K ~ dnorm( mu , sigma ),
        mu <- a + bB*B + bM*M,
        B ~ dnorm( nu , sigma_B ),
        c(a,nu) ~ dnorm( 0 , 0.5 ),
        c(bB,bM) ~ dnorm( 0, 0.5 ),
        sigma_B ~ dexp( 1 ),
        sigma ~ dexp( 1 )
    ) , data=dat_list , chains=4 , cores=4 )

## R code 15.18
precis( m15.3 , depth=2 )

## R code 15.19
obs_idx <- which( !is.na(d$neocortex.prop) )
dat_list_obs <- list(
    K = dat_list$K[obs_idx],
    B = dat_list$B[obs_idx],
    M = dat_list$M[obs_idx]
)
m15.4 <- ulam(
    alist(
        K ~ dnorm( mu , sigma ),
        mu <- a + bB*B + bM*M,
        B ~ dnorm( nu , sigma_B ),
        c(a,nu) ~ dnorm( 0 , 0.5 ),
        c(bB,bM) ~ dnorm( 0, 0.5 ),
        sigma_B ~ dexp( 1 ),
        sigma ~ dexp( 1 )
    ) , data=dat_list_obs , chains=4 , cores=4 )

## R code 15.22
m15.5 <- ulam(
    alist(
       # K as function of B and M
        K ~ dnorm( mu , sigma ),
        mu <- a + bB*B_merge + bM*M,

       # M and B correlation
        MB ~ multi_normal( c(muM,muB) , Rho_BM , Sigma_BM ),
        matrix[29,2]:MB <<- append_col( M , B_merge ),

       # define B_merge as mix of observed and imputed values
        vector[29]:B_merge <- merge_missing( B , B_impute ),

       # priors
        c(a,muB,muM) ~ dnorm( 0 , 0.5 ),
        c(bB,bM) ~ dnorm( 0, 0.5 ),
        sigma ~ dexp( 1 ),
        Rho_BM ~ lkj_corr(2),
        Sigma_BM ~ exponential(1)
    ) , data=dat_list , chains=4 , cores=4 )

m15.6 <- ulam(
    alist(
        # singing bird model
        ## cat known present/absent:
        notes|RC==0 ~ poisson( lambda ),
        log(lambda) <- a + b*cat,
        ## cat NA:
        notes|RC==1 ~ custom( log_sum_exp(
                log(k) + poisson_lpmf( notes | exp(a + b) ),
                log(1-k) + poisson_lpmf( notes | exp(a) )
            ) ),

        # priors
        a ~ normal(0,1),
        b ~ normal(0,0.5),

        # sneaking cat model
        cat|RC==0 ~ bernoulli(k),
        k ~ beta(2,2)
    ), data=dat , chains=4 , cores=4 )

## R code 15.31
m15.7 <- ulam(
    alist(
        # singing bird model
        notes|RC==0 ~ poisson( lambda ),
        notes|RC==1 ~ custom( log_sum_exp(
                log(k) + poisson_lpmf( notes | exp(a + b) ),
                log(1-k) + poisson_lpmf( notes | exp(a) )
            ) ),
        log(lambda) <- a + b*cat,
        a ~ normal(0,1),
        b ~ normal(0,0.5),

        # sneaking cat model
        cat|RC==0 ~ bernoulli(k),
        k ~ beta(2,2),

        # imputed values
        gq> vector[N]:PrC1 <- exp(lpC1)/(exp(lpC1)+exp(lpC0)),
        gq> vector[N]:lpC1 <- log(k) + poisson_lpmf( notes[i] | exp(a+b) ),
        gq> vector[N]:lpC0 <- log(1-k) + poisson_lpmf( notes[i] | exp(a) )
    ), data=dat , chains=4 , cores=4 )


text:
function vis(KS::KineticSolver, psi::Array{Float64,3})

    rho = computeMoments(KS, psi)
    #println("Final mass = $(sum(rho[:]))." )
    nicewrite(string(KS.outputfolderprefix, "data/rhofinal.txt"), rho)

    if KS.ictype == "linesource"
        fig, ax = subplots(figsize = (10.5, 8), dpi = 100)
        minVal = 0.0
        maxVal = 0.5
        #	title("$(KS.Q.quadtype), nq=$(KS.nquadpoints), cs=$(KS.convolutionmagnitude), rot=$(KS.rotationmagnitude)", fontsize=20)
        pcolormesh(rho, cmap = "RdBu_r", vmin = minVal, vmax = maxVal)
        colorbar()
        PyPlot.savefig(string(KS.outputfolderprefix, "rho.png"))

        rhoana = readdlm("../vis/exactLineSource.txt", ',')
        fig, ax = subplots(figsize = (10.5, 8), dpi = 100)
        plot(rhoana[:, 1], rhoana[:, 2])
        rhong = rho[3:end-2, 3:end-2] # remove ghosts
        cuthori = rhong[Int64(KS.ny / 2), :]
        x = range(-1.5, stop = 1.5, length = KS.nx - 4)
        plot(x, cuthori)
        cutverti = rhong[:, Int64(KS.nx / 2)]
        plot(x, cutverti)
        cutdiag = [rhong[i, i] for i = 1:KS.nx-4]
        plot(sqrt(2) * x, cutdiag)
        cutdiag = [rhong[KS.ny-4+1-i, i] for i = 1:KS.nx-4]
        plot(sqrt(2) * x, cutdiag)
        title(
            "$(KS.Q.quadtype), nq=$(KS.nquadpoints), convolution=$(KS.convolutionmagnitude) ",
            fontsize = 30,
        )
        PyPlot.savefig(string(KS.outputfolderprefix, "rhocut.png"))

    elseif KS.ictype == "checkerboard"
        fig, ax = subplots(figsize = (10.5, 8), dpi = 100)
        minVal = -7.0
        maxVal = 0.0
        title(
            "$(KS.Q.quadtype), nq=$(KS.nquadpoints), cs=$(KS.convolutionmagnitude)",
            fontsize = 20,
        )
        pcolormesh(log10.(abs.(rho)), cmap = "RdBu_r", vmin = minVal, vmax = maxVal)
        colorbar()
        PyPlot.savefig(string(KS.outputfolderprefix, "rho_log10.png"))

    elseif KS.ictype == "beam"

        fig, ax = subplots(figsize = (10.5, 8), dpi = 100)
        minVal = -7.0
        maxVal = 0.0
        title(
            "$(KS.Q.quadtype), nq=$(KS.nquadpoints), cs=$(KS.convolutionmagnitude)",
            fontsize = 20,
        )
        pcolormesh(log10.(abs.(rho[3:end-2, 3:end-2])), cmap = "RdBu_r")#,vmin=minVal, vmax=maxVal)
        colorbar()
        PyPlot.savefig(string(KS.outputfolderprefix, "rho_log10.png"))

        fig, ax = subplots(figsize = (10.5, 8), dpi = 100)
        minVal = -7.0
        maxVal = 0.0
        title(
            "$(KS.Q.quadtype), nq=$(KS.nquadpoints), cs=$(KS.convolutionmagnitude)",
            fontsize = 20,
        )
        pcolormesh(((rho[3:end-2, 3:end-2])), cmap = "RdBu_r")#,vmin=minVal, vmax=maxVal)
        colorbar()
        PyPlot.savefig(string(KS.outputfolderprefix, "rho.png"))

    end

    if KS.convolutionmagnitude > 0
        nicewrite(
            string(KS.outputfolderprefix, "data/convolutionmat.txt"),
            Matrix(KS.ScatteringKernelConvolution),
        )
        fig, ax = subplots(figsize = (10.5, 8), dpi = 100)
        title("Convolutionmatrix", fontsize = 20)
        pcolormesh((Matrix(KS.ScatteringKernelConvolution)), cmap = "RdBu_r")
        colorbar()
        PyPlot.savefig(string(KS.outputfolderprefix, "convmat.png"))
    end
end

function visMoments(KS::KineticSolver, phi, name::String)

    rho = phi[1, :, :]
    #println("Final mass = $(sum(rho[:]))." )
    #nicewrite(string(KS.outputfolderprefix,"data/rhofinalMoments",name,".txt"),rho)

    if KS.ictype == "linesource"
        fig, ax = subplots(figsize = (10.5, 8), dpi = 100)
        minVal = 0.0
        maxVal = 0.5
        #	title("$(KS.Q.quadtype), nq=$(KS.nquadpoints), cs=$(KS.convolutionmagnitude), rot=$(KS.rotationmagnitude)", fontsize=20)
        pcolormesh(rho, cmap = "RdBu_r", vmin = minVal, vmax = maxVal)
        colorbar()
        PyPlot.savefig(string(KS.outputfolderprefix, "rho", name, ".png"))

        rhoana = readdlm("../vis/exactLineSource.txt", ',')
        fig, ax = subplots(figsize = (10.5, 8), dpi = 100)
        plot(rhoana[:, 1], rhoana[:, 2])
        rhong = rho[3:end-2, 3:end-2] # remove ghosts
        cuthori = rhong[Int64(KS.ny / 2), :]
        x = range(-1.5, stop = 1.5, length = KS.nx - 4)
        plot(x, cuthori)
        cutverti = rhong[:, Int64(KS.nx / 2)]
        plot(x, cutverti)
        cutdiag = [rhong[i, i] for i = 1:KS.nx-4]
        plot(sqrt(2) * x, cutdiag)
        cutdiag = [rhong[KS.ny-4+1-i, i] for i = 1:KS.nx-4]
        plot(sqrt(2) * x, cutdiag)
        title(
            "$(KS.Q.quadtype), nq=$(KS.nquadpoints), convolution=$(KS.convolutionmagnitude) ",
            fontsize = 30,
        )
        PyPlot.savefig(string(KS.outputfolderprefix, "rhocut", name, ".png"))

    elseif KS.ictype == "checkerboard"
        fig, ax = subplots(figsize = (10.5, 8), dpi = 100)
        minVal = -7.0
        maxVal = 0.0
        title(
            "$(KS.Q.quadtype), nq=$(KS.nquadpoints), cs=$(KS.convolutionmagnitude)",
            fontsize = 20,
        )
        pcolormesh(log10.(abs.(rho)), cmap = "RdBu_r", vmin = minVal, vmax = maxVal)
        colorbar()
        PyPlot.savefig(string(KS.outputfolderprefix, "rho_log10.png"))

    elseif KS.ictype == "beam"

        fig, ax = subplots(figsize = (10.5, 8), dpi = 100)
        minVal = -7.0
        maxVal = 0.0
        title(
            "$(KS.Q.quadtype), nq=$(KS.nquadpoints), cs=$(KS.convolutionmagnitude)",
            fontsize = 20,
        )
        pcolormesh(log10.(abs.(rho[3:end-2, 3:end-2])), cmap = "RdBu_r")#,vmin=minVal, vmax=maxVal)
        colorbar()
        PyPlot.savefig(string(KS.outputfolderprefix, "rho_log10.png"))

        fig, ax = subplots(figsize = (10.5, 8), dpi = 100)
        minVal = -7.0
        maxVal = 0.0
        title(
            "$(KS.Q.quadtype), nq=$(KS.nquadpoints), cs=$(KS.convolutionmagnitude)",
            fontsize = 20,
        )
        pcolormesh(((rho[3:end-2, 3:end-2])), cmap = "RdBu_r")#,vmin=minVal, vmax=maxVal)
        colorbar()
        PyPlot.savefig(string(KS.outputfolderprefix, "rho.png"))

    end

    if KS.convolutionmagnitude > 0
        nicewrite(
            string(KS.outputfolderprefix, "data/convolutionmat.txt"),
            Matrix(KS.ScatteringKernelConvolution),
        )
        fig, ax = subplots(figsize = (10.5, 8), dpi = 100)
        title("Convolutionmatrix", fontsize = 20)
        pcolormesh((Matrix(KS.ScatteringKernelConvolution)), cmap = "RdBu_r")
        colorbar()
        PyPlot.savefig(string(KS.outputfolderprefix, "convmat.png"))
    end
end

function summarizeAsPlot(KS::KineticSolver, psi::Array{Float64,3})
    figure()
    subplot(221)
    pcolormesh(KS.SigmaS)
    subplot(222)
    pcolormesh(KS.SigmaT - KS.SigmaS)
    subplot(223)
    pcolormesh(computeMoments(KS, KS.Source))
    subplot(224)
    pcolormesh(computeMoments(KS, KS.psi0))
    PyPlot.savefig(string(KS.outputfolderprefix, "testcase.png"))

end

max_stars_repo_name: vavrines/CSN

max_issues_repo_name: vavrines/CSN

max_forks_repo_name: vavrines/CSN


text:
module TestVectors

using UnionArrays
using Test

function test_homogeneous_size()
    xs = UnionVector(Any[1, 2.0])
    @test xs[1] === 1
    @test xs[2] === 2.0
    @test (xs[1] = 3.0) isa Any
    @test xs[1] === 3.0
end

function test_heterogeneous_size()
    xs = UnionVector(Any[UInt8(1), 2.0, (a=1, b=2)])
    @test xs[1] === UInt8(1)
    @test xs[2] === 2.0
    @test xs[3] === (a=1, b=2)
    @test (xs[1] = 3.0) isa Any
    @test xs[1] === 3.0
end

end  # module


text:
module MonoMath

using Pipe
using EzXML
using DataFrames
using DelimitedFiles
using Mustache
using JSON3

include("samples.jl")
include("fonts.jl")

# Find which font is actually used to render a glyph (to find fallback font)
function pango_view(font, text)
  withenv("FC_DEBUG"=>4) do
    name = @pipe `pango-view --font=$font -q -t $text` |>
      eachline |>
      Iterators.filter(contains("family:"), _) |>
      foldl((x,y)->y, _) |> # Keep only last element
      split(_, '"')[2]
  end
end

# Get fc-match data as a Dict
function fc_match(family)
    props = Dict{String,String}()
    for line in eachline(`fc-match -v $family`)
        if ':' in line
            key, value = strip.(split(line, ':'))
            push!(props, key => value)
        end
    end
    return props
end

# Convert an option=>value pair to an hb-view flag (without leading --)
hbflag(option) = "$(replace(String(option[1]), '_'=>'-'))=$(option[2])"

function hbview(file, text, outfile; options...)
  flags = ["--"*hbflag(opt) for opt in options]
  run(`hb-view -o $outfile $flags $file $text`)
end

function default_shaper_list()
  if Sys.islinux()
    ["ot"]
  elseif Sys.iswindows()
    ["uniscribe", "ot"]
  elseif Sys.isapple()
    ["coretext", "ot"]
  else
    throw("Unsupported kernel: $(Sys.KERNEL)")
  end
end

# Get list of codepoints in a font
function font_codepoints(file)
  cmap = parsexml(read(pipeline(`ttx -y 0 -i -t cmap -o /dev/stdout $file`)))
  return [parse(UInt32, m["code"]) for m in findall("/ttFont/cmap//map", cmap)]
end

# Get versions stored in the head and name tables of the given font file
function font_versions(file)
  error_stream = IOBuffer()
  doc = parsexml(read(pipeline(`ttx -y 0 -i -t head -t name -o /dev/stdout $file`, stderr=error_stream)))

  # Print unexpected lines from error stream
  seekstart(error_stream)
  for line in eachline(error_stream)
    if !startswith(line, "Dumping ") # We expect lines starting with "Dumping "
      println(line)
    end
  end

  head_version = nodecontent(findfirst("/ttFont/head/fontRevision/@value", doc))
  name_version = nodecontent(findfirst("/ttFont/name/namerecord[@nameID=5]", doc)) |> strip
  return (head=head_version, name=name_version)
end

# Filename of a sample
sample_file(index, format) = lpad(index, 2, '0') * ".$format"

# Directory and filename of a sample, relative to the public root
function sample_path(; fontsize, os, shaper_list, format, fontname, index)
  shapers = join(shaper_list, '+')
  dir = "samples/$fontsize/$os/$shapers/$format/$fontname"
  file = sample_file(index, format)
  return (dir, file)
end

function make_sample(; fontname, fontfile, text, index, shaper_list, fontsize, format)
  (dir, file) = sample_path(; fontname, index, shaper_list, fontsize, format, os=Sys.KERNEL)
  fulldir = joinpath("public", dir)
  mkpath(fulldir)

  options = [
             :line_space => 2
             :margin => fontsize ÷ 2
             :background => "311d2f"
             :foreground => "e6e4e5"
             :language => "en"
             :shapers => join(shaper_list, ",")
             :font_size => fontsize
            ];

  # Write image
  outfile = joinpath(fulldir, file)
  println("Writing $outfile...")
  hbview(fontfile, text, outfile; options...)
end

function make_font_samples(; fontname, fontfile, options...)
  for (index, (desc, text)) in enumerate(all_samples)
    make_sample(; fontname, fontfile, text, index, options...)
  end
end

function make_samples(; format=:svg, fontsize=14, shaper_list=default_shaper_list())
  # Make sample text files
  for (index, (desc, text)) in enumerate(all_samples)
    write(joinpath("public", "samples", sample_file(index, :txt)), text)
  end

  for (name, file) in eachrow(fonts)
    make_font_samples(; fontname=name, fontfile=file, format, fontsize, shaper_list)
  end
end

function make_html(; sample_indices=eachindex(text_samples), fontsize=14, os=Sys.KERNEL,
                     shaper_list=default_shaper_list(), format=:svg, fontname="JuliaMono", filename="index.html")
  args = (; fontsize, os, shaper_list, format, fontname)
  img = [joinpath(sample_path(; args..., index)...) for index in sample_indices]
  txt = [joinpath("samples", sample_file(index, :txt)) for index in sample_indices]
  desc = first.(all_samples[sample_indices])
  samples = DataFrame(; img, txt, desc)
  tok = Mustache.load("resources/$filename.template")
  str = render(tok; fonts, samples)
  write("public/$filename", str)
  return str
end

rounded_version(v::Number) = round(v, digits=4)
rounded_version(v::AbstractString) = rounded_version(parse(Float64, v))

# Convert maj.min.patch version to a single number, as used by Iosevka
function parse_maj_min_patch(str)
    println("@@ parsing $str...")
    result = match(r"^(\d+)\.(\d+)\.(\d+)(.*)$", str)
    isnothing(result) && return nothing
    maj, min, patch, suffix = result
    isempty(suffix) || println("Version suffix ignored in $str")
    return parse(Int, maj) + parse(Int, min)/100 + parse(Int, patch)/1000
end

function make_font_data()
  data = Dict{String,String}()

  for (name, file) in eachrow(fonts)
    versions = font_versions(file)

    # First version in semicolon-spearated list
    name_first_version = first(split(versions.name, ';', limit=2))

    # Remove Version prefix if present
    name_value = last(split(name_first_version, "Version ", limit=2))

    head_rounded = rounded_version(versions.head)

    # Name table version: first try to parse as simple number
    name_rounded = try
      rounded_version(name_value)
    catch
      # Try to parse as x.y.z
      v = parse_maj_min_patch(name_value)
      isnothing(v) ? nothing : rounded_version(v)
    end
    println("Versions for $name: head $head_rounded, name $name_rounded")

    same = head_rounded == name_rounded
    if same || versions.head == "0.0"
      summary = name_value
    else
      summary  = name_value * " (head table: $(versions.head))"
    end

    data[name] = summary
  end

  write("public/font-data.js", "font_data = ", JSON3.write(data))
end

function make_all()
  make_html()
  make_font_data()
  make_samples()

  extra_indices = length(text_samples) .+ eachindex(extra_samples)
  make_html(sample_indices=extra_indices, filename="extras.html")
end

end

max_stars_repo_name: waldyrious/mono-math

max_issues_repo_name: waldyrious/mono-math

max_forks_repo_name: waldyrious/mono-math


text:
################################################################################

# read file
file = "data/signal.tsv"
v = readdlm(file)

# setup hidden Markov model
hmm = setup(v)

################################################################################

max_stars_repo_name: DanielRivasMD/HiddenMarkovModelReaders

max_issues_repo_name: DanielRivasMD/HiddenMarkovModelReaders