text:
module Optimizers
    export iterations
    using NLopt

    abstract type AbstractOptimizer end
    abstract type AbstractGradientDescent <: AbstractOptimizer end
    struct Steepest <: AbstractGradientDescent
        ∇f
        iterations::Integer
        rate::AbstractFloat
        Steepest(∇f; iterations=1, rate=1e-3) = new(∇f, iterations, rate)
    end
    (stp::Steepest)(x0; kwargs...) = optimize(stp, x0; kwargs...)

    struct Momentum <: AbstractGradientDescent
        ∇f
        iterations::Integer
        rate::AbstractFloat
        β::AbstractFloat
        Momentum(∇f; iterations=1, rate=1e-3, beta=1e-8) = new(∇f, iterations, rate, beta)
    end
    (stp::Momentum)(x0; kwargs...) = optimize(stp, x0; kwargs...)

    struct AdaGrad <: AbstractGradientDescent
        ∇f
        iterations::Integer
        rate::AbstractFloat
        ϵ::AbstractFloat
        AdaGrad(∇f; iterations=1, rate=1e-3, epsilon=1e-8) = new(∇f, iterations, rate, epsilon)
    end
    (stp::AdaGrad)(x0; kwargs...) = optimize(stp, x0; kwargs...)

    struct Adam <: AbstractGradientDescent
        ∇f
        iterations::Integer
        rate::AbstractFloat
        β1::AbstractFloat
        β2::AbstractFloat
        ϵ::AbstractFloat

        Adam(∇f; iterations=1, rate=0.001, beta1=0.9, beta2=0.999, epsilon=1e-8) = new(iterations, rate, beta1, beta2, epsilon)
    end
    (stp::Adam)(x0; kwargs...) = optimize(∇f, stp, x0; kwargs...)

    struct GlobalOpt <: AbstractOptimizer
        f
        alg::Symbol
        iterations::Integer
        xtolrel::Real
        bounds::NTuple{2}
        GlobalOpt(f; alg=:GN_MLSL_LDS, iterations=1, xtolrel=eps(), bounds=(-Float64(pi), Float64(pi))) = new(f, alg, iterations, xtolrel, bounds)
    end

    struct CRS <: AbstractOptimizer
        f
        iterations::Integer
        xtolrel::Real
        bounds::NTuple{2}
        CRS(f; iterations=1, xtolrel=eps(), bounds=(-Float64(pi), Float64(pi))) = new(f, iterations, xtolrel, bounds)
    end

    struct SGD <: AbstractGradientDescent
        ∇f
        iterations::Integer
        rate::Function
        SGD(∇f; iterations=1, rate=t->1e-3) = new(∇f, iterations, rate)
    end
    (stp::SGD)(x0; kwargs...) = optimize(stp, x0; kwargs...)

    iterations(abopt::AbstractOptimizer) = abopt.iterations

    function optimize(abopt::AbstractGradientDescent, x0::AbstractArray; isverbose::Bool=false, kwargs...)
        x = x0
        state = initialstate(abopt, x0; kwargs...)
        for k = 1:iterations(abopt)
            grad = gradient(abopt, x0)
            Δx, state = updatestep(abopt, grad, state, k; kwargs...)
            x = x - Δx

            # if norm(Δx) < eps; break
            # isverbose && println("#Iter. = $k, loss = $(f(vecpm_opt)), ||∇loss|| = $(norm(Δvecpm))")
        end
        x
    end

    gradient(abopt::AbstractOptimizer, x0) = abopt.∇f(x0)

    initialstate(::Steepest, args...; kwargs...) = ()
    initialstate(::AdaGrad, args...; kwargs...) = 0
    initialstate(adam::Adam, x0; kwargs...) = (zero(x0), 0)
    initialstate(sgd::SGD, args...; epoch=1, kwargs...) = (sgd.rate(epoch))

    updatestep(stp::Steepest, grad, args...; kwargs...) = (stp.rate .* grad, ())
    function updatestep(agd::AdaGrad, grad, v, itr; kwargs...)
        v = v + norm(grad)^2
        Δx = (agd.rate / (sqrt(v) + agd.ϵ)) * grad
        (upm, v)
    end
    function updatestep(adam::Adam, grad, (m, v), itr; kwargs...)
        m = adam.β1 * m + (1 - adam.β1) * grad
        v = adam.β2 * v + (1 - adam.β2) * norm(grad)^2
        m̂ = m / (1 - adam.β1^(itr))
        v̂ = v / (1 - adam.β2^(itr))

        Δx = adam.rate * m̂ / (sqrt(v̂) + adam.ϵ)
        (Δx, (m, v))
    end
    function updatestep(sgd::SGD, grad, rate, itr; kwargs...)
        (rate .* grad, rate)
    end

    function (gopt::GlobalOpt)(x0; kwargs...)
        opt = Opt(gopt.alg, length(x0))
        lower_bounds!(opt, gopt.bounds[1]*ones(size(x0)))
        upper_bounds!(opt, gopt.bounds[2]*ones(size(x0)))
        xtol_rel!(opt, gopt.xtolrel)
        xtol_abs!(opt, eps())
        maxeval!(opt, gopt.iterations)

        min_objective!(opt, gopt.f)
        minf, x_opt, ret = NLopt.optimize(opt, x0)

        return x_opt
    end

    function (gopt::CRS)(x0; kwargs...)
        opt = Opt(:GN_CRS2_LM, length(x0))
        lower_bounds!(opt, gopt.bounds[1]*ones(size(x0)))
        upper_bounds!(opt, gopt.bounds[2]*ones(size(x0)))
        xtol_rel!(opt, gopt.xtolrel)
        xtol_abs!(opt, eps())
        maxeval!(opt, gopt.iterations)

        min_objective!(opt, gopt.f)
        minf, x_opt, ret = NLopt.optimize(opt, x0)

        return x_opt
    end
end

max_stars_repo_name: nixir/MDCDL

max_issues_repo_name: nixir/MDCDL

max_forks_repo_name: nixir/MDCDL


text:
# ---
# title: 1403. Minimum Subsequence in Non-Increasing Order
# id: problem1403
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Greedy, Sort
# link: <https://leetcode.com/problems/minimum-subsequence-in-non-increasing-order/description/>
# hidden: true
# ---
#
# Given the array `nums`, obtain a subsequence of the array whose sum of
# elements is **strictly greater** than the sum of the non included elements in
# such subsequence.
#
# If there are multiple solutions, return the subsequence with **minimum size**
# and if there still exist multiple solutions, return the subsequence with the
# **maximum total sum** of all its elements. A subsequence of an array can be
# obtained by erasing some (possibly zero) elements from the array.
#
# Note that the solution with the given constraints is guaranteed to be
# **unique**. Also return the answer sorted in **non-increasing** order.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [4,3,10,9,8]
#     Output: [10,9]
#     Explanation: The subsequences [10,9] and [10,8] are minimal such that the sum of their elements is strictly greater than the sum of elements not included, however, the subsequence [10,9] has the maximum total sum of its elements.
#
#
# **Example 2:**
#
#
#
#     Input: nums = [4,4,7,6,7]
#     Output: [7,7,6]
#     Explanation: The subsequence [7,7] has the sum of its elements equal to 14 which is not strictly greater than the sum of elements not included (14 = 4 + 4 + 6). Therefore, the subsequence [7,6,7] is the minimal satisfying the conditions. Note the subsequence has to returned in non-decreasing order.
#
#
# **Example 3:**
#
#
#
#     Input: nums = [6]
#     Output: [6]
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 500`
#   * `1 <= nums[i] <= 100`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end

ext: jl


text:
using OrdinaryDiffEq, DiffEqBase, DiffEqDevTools, Test
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
using DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear

f = (u,p,t) -> u
fun = ODEFunction(f; analytic=(u0,p,t)->u0*exp(t))
prob = ODEProblem(fun,1/2,(0.0,1.0))

sol =solve(prob,Euler();dt=1//2^(4),dense_errors=true)
sol2 =solve(prob,Vern9();dt=1//2^(10),abstol=1e-14,reltol=1e-14)

prob2 = ODEProblem(f,1/2,(0.0,1.0))
sol3 =solve(prob_ode_linear,Euler();dt=1//2^(4))

errsol1 = appxtrue(sol,sol2)

sol4 =solve(prob,Euler();dt=1//2^(4))
test_sol = TestSolution(sol2)
errsol2 = appxtrue(sol4,test_sol)

sol5 =solve(prob,Euler();dt=1//2^(4))
test_sol = TestSolution(sol2.t,sol2[end])
errsol3 = appxtrue(sol5,test_sol)

@test errsol1.errors[:L2] ≈ 0.018865798306718855
@test errsol1.errors[:L2] ≈ errsol2.errors[:L2]

@test errsol1.errors[:final] ≈ sol.errors[:final]
@test errsol1.errors[:l2] ≈ sol.errors[:l2]
@test errsol1.errors[:l∞] ≈ sol.errors[:l∞]
@test errsol1.errors[:L∞] ≈ sol.errors[:L∞]
@test errsol1.errors[:L2] ≈ sol.errors[:L2]


text:
module ChebyExp
using Base.LinAlg: BLAS.gemm!, BlasFloat, checksquare, tilebufsize, Abuf, Bbuf, Cbuf, rcswap!
import Base.A_mul_B!

export chebyexp

include("exp.jl")

# Utility functions
A_mul_B!(α::T, A::StridedVecOrMat{T}, B::StridedVecOrMat{T}, β::T, C::StridedMatrix{T}) where {T<:BlasFloat} = gemm!('N','N',α,A,B,β,C)

function A_mul_B!(α::Number, A::AbstractVecOrMat{T}, B::AbstractVecOrMat{S}, β::Number, C::AbstractVecOrMat{R}) where {T,S,R}
    mA, nA = size(A)
    mB, nB = size(B)
    if mB != nA
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA), matrix B has dimensions ($mB,$nB)"))
    end
    if size(C,1) != mA || size(C,2) != nB
        throw(DimensionMismatch("result C has dimensions $(size(C)), needs ($mA,$nB)"))
    end

    @inbounds begin
        if isbits(R) && isbits(T) && isbits(S)
            tile_size = floor(Int, sqrt(tilebufsize / max(sizeof(R), sizeof(S), sizeof(T))))
            sz = (tile_size, tile_size)
            Atile = unsafe_wrap(Array, convert(Ptr{T}, pointer(Abuf)), sz)
            Btile = unsafe_wrap(Array, convert(Ptr{S}, pointer(Bbuf)), sz)

            z1 = zero(A[1, 1]*B[1, 1] + A[1, 1]*B[1, 1])
            z = convert(promote_type(typeof(z1), R), z1)

            if mA < tile_size && nA < tile_size && nB < tile_size
                Base.copy_transpose!(Atile, 1:nA, 1:mA, A, 1:mA, 1:nA)
                copy!(Btile, 1:mB, 1:nB, B, 1:mB, 1:nB)
                for j = 1:nB
                    boff = (j-1)*tile_size
                    for i = 1:mA
                        aoff = (i-1)*tile_size
                        s = z
                        for k = 1:nA
                            s += Atile[aoff+k] * Btile[boff+k]
                        end
                        C[i,j] = α*s + β*C[i,j]
                    end
                end
            else
                Ctile = unsafe_wrap(Array, convert(Ptr{R}, pointer(Cbuf)), sz)
                for jb = 1:tile_size:nB
                    jlim = min(jb+tile_size-1,nB)
                    jlen = jlim-jb+1
                    for ib = 1:tile_size:mA
                        ilim = min(ib+tile_size-1,mA)
                        ilen = ilim-ib+1
                        fill!(Ctile, z)
                        for kb = 1:tile_size:nA
                            klim = min(kb+tile_size-1,mB)
                            klen = klim-kb+1
                            Base.copy_transpose!(Atile, 1:klen, 1:ilen, A, ib:ilim, kb:klim)
                            copy!(Btile, 1:klen, 1:jlen, B, kb:klim, jb:jlim)
                            for j = 1:jlen
                                bcoff = (j-1)*tile_size
                                for i = 1:ilen
                                    aoff = (i-1)*tile_size
                                    s = z
                                    for k = 1:klen
                                        s += Atile[aoff+k] * Btile[bcoff+k]
                                    end
                                    Ctile[bcoff+i] += s
                                end
                            end
                        end
                        for j = 1:jlen, i = 1:ilen
                            coff = (j-1)*tile_size+i
                            Ctile[coff] = α*Ctile[coff] + β*C[ib+i-1,jb+j-1]
                        end
                        copy!(C, ib:ilim, jb:jlim, Ctile, 1:ilen, 1:jlen)
                    end
                end
            end
        else # Multiplication for non-plain-data uses the naive algorithm
            for i = 1:mA, j = 1:nB
                z2 = zero(A[i, 1]*B[1, j] + A[i, 1]*B[1, j])
                Ctmp = convert(promote_type(R, typeof(z2)), z2)
                for k = 1:nA
                    Ctmp += A[i, k]*B[k, j]
                end
                C[i,j] = α*Ctmp + β*C[i,j]
            end
        end
    end # @inbounds
    return C
end

end

max_stars_repo_name: jebej/ChebyExp

max_issues_repo_name: jebej/ChebyExp

max_forks_repo_name: jebej/ChebyExp


text:
#### Extract subset
# -------------------------------------------------------------------------------------------------
#
# Subsetting, meaning extracting some portion of the 3D domain
#
# -------------------------------------------------------------------------------------------------
function extract_subset_dim1_3d(
    ncol,
    nlay,
    ngpt,
    array_in,
    colS,
    colE,
    array_out,
)#
    #    bind (C, name="extract_subset_dim1_3d")
    #    integer,                             intent(in ) :: ncol, nlay, ngpt
    #    real(FT), dimension(ncol,nlay,ngpt), intent(in ) :: array_in
    #    integer,                             intent(in ) :: colS, colE
    #    real(FT), dimension(colE-colS+1,
    #                             nlay,ngpt), intent(out) :: array_out

    #    integer :: icol, ilay, igpt
    for igpt = 1:ngpt
        for ilay = 1:nlay
            for icol = colS:colE
                array_out[icol-colS+1, ilay, igpt] = array_in[icol, ilay, igpt]
            end
        end
    end

end
# ---------------------------------
function extract_subset_dim2_4d(
    nmom,
    ncol,
    nlay,
    ngpt,
    array_in,
    colS,
    colE,
    array_out,
)#
    #    bind (C, name="extract_subset_dim2_4d")
    #    integer,                                  intent(in ) :: nmom, ncol, nlay, ngpt
    #    real(FT), dimension(nmom,ncol,nlay,ngpt), intent(in ) :: array_in
    #    integer,                                  intent(in ) :: colS, colE
    #    real(FT), dimension(nmom,colE-colS+1,
    #                                  nlay,ngpt), intent(out) :: array_out

    #    integer :: icol, ilay, igpt, imom

    for igpt = 1:ngpt
        for ilay = 1:nlay
            for icol = colS:colE
                for imom = 1:nmom
                    array_out[imom, icol-colS+1, ilay, igpt] =
                        array_in[imom, icol, ilay, igpt]
                end
            end
        end
    end

end

#
# Extract the absorption optical thickness which requires mulitplying by 1 - ssa
#
#    integer,                             intent(in ) :: ncol, nlay, ngpt
#    real(FT), dimension(ncol,nlay,ngpt), intent(in ) :: tau_in, ssa_in
#    integer,                             intent(in ) :: colS, colE
#    real(FT), dimension(colE-colS+1,
#                             nlay,ngpt), intent(out) :: tau_out
#    integer :: icol, ilay, igpt
function extract_subset_absorption_tau(
    ncol,
    nlay,
    ngpt,
    tau_in,
    ssa_in,
    colS,
    colE,
    tau_out,
)
    FT = eltype(tau_in)

    for igpt = 1:ngpt
        for ilay = 1:nlay
            for icol = colS:colE
                tau_out[icol-colS+1, ilay, igpt] =
                    tau_in[icol, ilay, igpt] *
                    (FT(1) - ssa_in[icol, ilay, igpt])
            end
        end
    end
end

# Extract a subset of n columns starting with column 'start'
function get_subset_range(
    this::GasConcs,
    start::I,
    n::I,
    subset::GasConcs,
) where {I}
    # class(GasConcs),      intent(in   ) :: this
    # integer,                  intent(in   ) :: start, n
    # class(GasConcs),      intent(inout) :: subset
    @assert !(n <= 0)
    @assert !(start < 1)
    @assert !(
        this.ncol ≠ nothing && start > this.ncol || start + n - 1 > this.ncol
    )

    this.nlay = 0
    this.ncol = 0
    # reset!(subset)
    # These two arrays should be the same length
    subset.gas_name = Array(undef, size(this.gas_name))
    subset.concs = Array(undef, size(this.concs))
    subset.nlay = this.nlay
    subset.ncol = this.ncol ≠ nothing ? n : 0
    subset.gas_name[:] .= this.gas_name[:]

    for i = 1:size(this.gas_name)
        #
        # Preserve scalar/1D/2D representation in subset,
        #   but need to ensure at least extent 1 in col dimension (ncol = 0 means no gas exploits this dimension)
        #
        allocate(subset.concs[i].conc(
            min(max(subset.ncol, 1), size(this.concs[i].conc, 1)),
            min(subset.nlay, size(this.concs[i].conc, 2)),
        ))
        if size(this.concs[i].conc, 1) > 1      # Concentration stored as 2D
            subset.concs[i].conc[:, :] .=
                this.concs[i].conc[start:(start+n-1), :]
        else
            subset.concs[i].conc[:, :] .= this.concs[i].conc[:, :]
        end
    end

end


text:
# f(x,y)=x-y'

FXMYtranspose(x::Float64,y::Float64)=Fxmy(x,y)
FXMYtranspose(x::Float64,y)=Fxmy(x,y')
FXMYtranspose(x,y::Float64)=Fxmy(x,y)
FXMYtranspose(x,y)=Fxmy(x,y')
FXMYtranspose_inplace(value,auxvalue,x,y)=Fxmy_inplace(value,auxvalue,x,y')

function DXMYtranspose(derivativeIDX,f_c,faux_c,grad_c,grad_n,x::Array,y::Array)
    if derivativeIDX==1
        if size(x)==(1,1)
            axpy!(1.0,[sum(grad_c)],grad_n)
        else
            axpy!(1.0,grad_c,grad_n)
        end
    elseif derivativeIDX==2
        if size(y)==(1,1)
            axpy!(-1.0,[sum(grad_c)],grad_n)
        else
            axpy!(-1.0,grad_c',grad_n)
        end
    end
end

if PROC=="GPU"

    function FXMYtranspose_inplace(value::CudaArray,auxvalue,x::CudaArray,y::CudaArray)
        if size(x)==(1,1)
            gfill!(value,x);
            CUBLAS.geam!('N','T',1.0,value,-1.0,y,value)
        elseif size(y)==(1,1)
            gfill!(value,y)
            CUBLAS.geam!('N','N',1.0,x,-1.0,value,value)
        else
            CUBLAS.geam!('N','T',1.0,x,-1.0,y,value)
        end
    end

    function DXMYtranspose(derivativeIDX,f_c,faux_c,grad_c,grad_n,x::CudaArray,y::CudaArray)
        if derivativeIDX==1
            if size(x)==(1,1)
                axpy!(1.0,sum(grad_c),grad_n)
            else
                tmp=CudaArray(Float64,size(grad_c))
                CUBLAS.geam!('T','N',1.0,grad_c,0.0,tmp,tmp) # tmp=grad_c'
                axpy!(1.0,tmp,grad_n); free(tmp)
            end
        elseif derivativeIDX==2
            if size(y)==(1,1)
                axpy!(-1.0,sum(grad_c),grad_n)
            else
                tmp=CudaArray(Float64,size(grad_c))
                CUBLAS.geam!('T','N',1.0,grad_c,0.0,tmp,tmp) # tmp=grad_c'
                axpy!(-1.0,tmp,grad_n); free(tmp)
            end
        end
    end

end

Derivative[FXMYtranspose]=DXMYtranspose
Inplace[FXMYtranspose]=FXMYtranspose_inplace

export FXMYtranspose

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
include("normal_cases.jl")
include("exception_cases.jl")


text:
module SimpleTasks

include("types.jl")

include("tasks/include.jl")
include("services/include.jl")
include("examples/include.jl")

end # module SimpleTasks

max_forks_repo_name: seung-lab/SimpleTasks


text:
#Kronecker delta function
#23/07/2017

function krondelt(ii,jj)
  if ii==jj
    kronecker = 1
  else
    kronecker = 0
  end
  return kronecker
end

max_stars_repo_name: CSChisholm/Rydberg

max_issues_repo_name: CSChisholm/Rydberg

max_forks_repo_name: CSChisholm/Rydberg


text:
abstract type AbstractGridWorldAction end

abstract type AbstractMoveAction <: AbstractGridWorldAction end
abstract type AbstractTurnAction <: AbstractGridWorldAction end

struct MoveForward <: AbstractMoveAction end
const MOVE_FORWARD = MoveForward()
move(::MoveForward, dir::AbstractDirection, pos::CartesianIndex{2}) = move(dir, pos)

struct MoveUp <: AbstractMoveAction end
const MOVE_UP = MoveUp()
move(::MoveUp, pos::CartesianIndex{2}) = move(UP, pos)
move(action::MoveUp, dir::AbstractDirection, pos::CartesianIndex{2}) = move(action, pos)

struct MoveDown <: AbstractMoveAction end
const MOVE_DOWN = MoveDown()
move(::MoveDown, pos::CartesianIndex{2}) = move(DOWN, pos)
move(action::MoveDown, dir::AbstractDirection, pos::CartesianIndex{2}) = move(action, pos)

struct MoveLeft <: AbstractMoveAction end
const MOVE_LEFT = MoveLeft()
move(::MoveLeft, pos::CartesianIndex{2}) = move(LEFT, pos)
move(action::MoveLeft, dir::AbstractDirection, pos::CartesianIndex{2}) = move(action, pos)

struct MoveRight <: AbstractMoveAction end
const MOVE_RIGHT = MoveRight()
move(::MoveRight, pos::CartesianIndex{2}) = move(RIGHT, pos)
move(action::MoveRight, dir::AbstractDirection, pos::CartesianIndex{2}) = move(action, pos)

struct NoMove <: AbstractMoveAction end
const NO_MOVE = NoMove()
move(::NoMove, pos::CartesianIndex{2}) = pos
move(action::NoMove, dir::AbstractDirection, pos::CartesianIndex{2}) = move(action, pos)

struct TurnLeft <: AbstractTurnAction end
const TURN_LEFT = TurnLeft()
turn(::TurnLeft, ::Up) = LEFT
turn(::TurnLeft, ::Down) = RIGHT
turn(::TurnLeft, ::Left) = DOWN
turn(::TurnLeft, ::Right) = UP

struct TurnRight <: AbstractTurnAction end
const TURN_RIGHT = TurnRight()
turn(::TurnRight, ::Up) = RIGHT
turn(::TurnRight, ::Down) = LEFT
turn(::TurnRight, ::Left) = UP
turn(::TurnRight, ::Right) = DOWN

struct PickUp <: AbstractGridWorldAction end
const PICK_UP = PickUp()

struct Drop <: AbstractGridWorldAction end
const DROP = Drop()

const DIRECTED_NAVIGATION_ACTIONS = (MOVE_FORWARD, TURN_LEFT, TURN_RIGHT)
const UNDIRECTED_NAVIGATION_ACTIONS = (MOVE_UP, MOVE_DOWN, MOVE_LEFT, MOVE_RIGHT)


text:
using DifferentialEquations

intro1 = function()
    f(u,p,t) = 1.01*u
    u0 = 1/2
    tspan = (0.0,1.0)

    prob = ODEProblem(f,u0,tspan)
    using Plots
    plot(sol,linewidth=5,title="Solution to the linear ODE with a thick line",
        xaxis="Time (t)",yaxis="u(t) (in μm)",label="My Thick Line!") # legend=false
    plot!(sol.t, t->0.5*exp(1.01t),lw=3,ls=:dash,label="True Solution!")
end

function mtk1
    using ModelingToolkit

    @variables t x(t) RHS(t)  # independent and dependent variables
    @parameters τ       # parameters
    D = Differential(t) # define an operator for the differentiation w.r.t. time

    # your first ODE, consisting of a single equation, indicated by ~
    @named fol_separate = ODESystem([ RHS  ~ (1 - x)/τ,
                                      D(x) ~ RHS ])

    using DifferentialEquations: solve
    using Plots: plot

    prob = ODEProblem(structural_simplify(fol_separate), [x => 0.0], (0.0,10.0), [τ => 3.0])
    sol = solve(prob)
    plot(sol, vars=[x,RHS])
    #sol[RHS] # get the derived variable as Array
    #hcat(sol.t, (sol[v] for v in (x,RHS))...)
end

function lorenz!(du,u,p,t)
    σ,ρ,β = p
    du[1] = σ*(u[2]-u[1])
    du[2] = u[1]*(ρ-u[3]) - u[2]
    du[3] = u[1]*u[2] - β*u[3]
end
u0 = [1.0,0.0,0.0]
p = (10,28,8/3) # we could also make this an array, or any other type!
tspan = (0.0,100.0)
prob = ODEProblem(lorenz!,u0,tspan,p)
sol = solve(prob)

function interactive()
    using Plots
    plot(sol)
end

function mtk_example()
    #using ModelingToolkit
    @parameters t σ ρ β
    @variables x(t) y(t) z(t)
    D = Differential(t)
    eqs = [D(x) ~ σ*(y-x),
        D(y) ~ x*(ρ-z)-y,
        D(z) ~ x*y - β*z]
    de = ODESystem(eqs,t,[x,y,z],[σ,ρ,β])
    sol = solve(prob)
    plot(sol)
end

function catalyst_example1()
    #using Catalyst, Latexify
    rs = @reaction_network begin
    c1, X --> 2X
    c2, X --> 0
    c3, 0 --> X
    end c1 c2 c3
    #Graph(rs)
    p = (1.0,2.0,50.) # [c1,c2,c3]
    tspan = (0.,4.)
    u0 = [5.]         # [X]
    #osys = convert(ODESystem,rs)
    #latexify(osys)
    # solve ODEs
    oprob = ODEProblem(rs, u0, tspan, p)
    osol  = solve(oprob, Tsit5())
    #plot(osol)
    equations(osys)

    # solve for Steady-States
    #ssprob = SteadyStateProblem(rs, u0, p)
    ssprob = SteadyStateProblem(osys, u0, p)
    sssol  = solve(ssprob, SSRootfind())
end

function simpmod_example()
    #using ModelingToolkit, DifferentialEquations, Plots
    @parameters t ks km kd Y HS HB
    # t in hours
    # masses per soilmass in mol/g
    @variables s(t) b(t) cr(t) r_tot(t) q(t) dec_s(t) tvr_b(t)
    D = Differential(t)
    HG = HS - Y * HB # can take out of system because does not involve t
    eqs = [
        dec_s ~ (ks*s*b)/(s + km),
        tvr_b ~ kd*b,
        D(s) ~ -dec_s + tvr_b,
        D(b) ~ Y * dec_s - tvr_b,
        r_tot ~ (1-Y) * dec_s,
        D(cr) ~ r_tot,
        q ~ -HG * dec_s,
        ]
    de = ODESystem(eqs; name=:simpmod)
    des = structural_simplify(de) # omit r,q
    p_straw = Dict(
        :parms => [ks => 0.15, km => 2.62e-6, kd => 1.35e-2, Y => 0.72,
            HB => -492, HS => -469],
        #:x0 => [s => 4.17e-6, b => 1.14e-5, cr => 0] # twutz: wrong in paper
        :x0 => [s => 4.17e-5, b => 1.14e-5, cr => 0]
    ) # see S5
    prob = ODEProblem(des,
        #[s => 1.0, b => 1.0, r => 0.0, q => 0.0],
        p_straw[:x0],
        (0.0, 48.0),
        #[ks => 0.5, km => 0.5, kd => 1/20, Y => 0.5, HG => -10]
        p_straw[:parms]
    )
    #sol = solve(prob);
    #sol = solve(prob,AutoTsit5(Rodas5()));
    #sol = solve(prob,lsoda()); # does not work with DAE
    sol = solve(prob,Rodas5());
    #sol[r_tot]
    #plot(sol)
    plot(sol,vars=[s, b, cr, dec_s, r_tot])
    # glucose taken up within 16 hours, but stored in microbial biomass
    # instead of respired
    plot(sol,vars=[q], ylim = (0, 800e-6))
    plot(sol,vars=[r_tot], ylim = (0,2e-6))
    # in the model resp and q only differ by factor (1-Y)/HG
end

function test_min()
    # using minimum function in equation works nicely
    #using ModelingToolkit, DifferentialEquations, Plots
    @parameters t ks km kd Y HG cn_s cn_b imm_N
    @variables s(t) b(t) r(t) rO(t) q(t) dec_s(t) tvr_b(t)
    @variables syn_c_pot(t) syn_n_pot(t) syn(t)
    D = Differential(t)
    HG = HS - Y * HB
    eqs = [
        dec_s ~ (ks*s*b)/(s + km),
        tvr_b ~ kd*b,
        syn_c_pot ~ Y * dec_s,
        syn_n_pot ~ (dec_s/cn_s + imm_N) * cn_b,
        syn ~ min(syn_c_pot, syn_n_pot), # note the min function
        D(s) ~ -dec_s + tvr_b,
        D(b) ~ +syn - tvr_b,
        r ~ dec_s - syn,
        rO ~ r - (1 - Y) * dec_s,
        q ~ -HG * dec_s,
        ]
    de = ODESystem(eqs)
    des = structural_simplify(de) # omit r,q
    prob = ODEProblem(des,
        #[s => 1.0, b => 1.0, r => 0.0, q => 0.0],
        [s => 100.0, b => 1.0],
        (0.0, 50.0),
        [ks => 0.5, km => 0.5, kd => 1/10, Y => 0.5, HG => -10,
        cn_b => 8, cn_s => 20,
        imm_N => 0.05];
        #jac=true # takes longer but ok here
    )
    sol = solve(prob);
    #plot(sol)
    plot(sol,vars=[r, rO, syn, syn_c_pot, syn_n_pot]) # r and q are tracked
    plot(sol,vars=[s, b, r, q, dec_s]) # r and q are tracked
end


text:
using Compat: @compat
using Test, UUIDs, Dates, Onda, Legolas, Arrow, Tables, TimeSpans, DataFrames, Random
using ConstructionBase

function has_rows(a, b)
    for name in propertynames(b)
        getproperty(a, name) == getproperty(b, name) || return false
    end
    for name in Tables.columnnames(b)
        Tables.getcolumn(a, name) == Tables.getcolumn(b, name) || return false
    end
    return true
end

include("annotations.jl")
include("signals.jl")
include("serialization.jl")
include("samples.jl")
include("deprecations.jl")
include(joinpath(dirname(@__DIR__), "examples", "flac.jl"))
include(joinpath(dirname(@__DIR__), "examples", "tour.jl"))


text:
# Path integral molecular dynamics integration.

function make_Cmat(P::Int)
    Cmat = Array{Float64}(undef, P, P)
    # This works for both even and odd P.
    M = P / 2

    for k_idx in 1:P
        k = k_idx - 1

        for j in 1:P
            if k == 0
                Cmat[j, k_idx] = 1.0
            elseif 1 <= k < M
                Cmat[j, k_idx] = sqrt(2.0) * cos(2pi * j * k / P)
            elseif k == M
                Cmat[j, k_idx] = (-1.0)^j
            elseif M < k <= P-1
                Cmat[j, k_idx] = sqrt(2.0) * sin(2pi * j * k / P)
            else
                error()
            end
        end
    end

    Cmat ./ sqrt(P)
end

MS_fict(P::Int) = MS ./ P # g / mol

struct NormalModeTransformer
    "Transformation matrix."
    Cmat::Matrix{Float64}
end

function to_normal_modes!(nms::NormalModeState{P,D,A,N}, s::State{P,D,A,N},
                          nmt::NormalModeTransformer) where {P,D,A,N}
    for n in 1:N
        for a in 1:A
            for d in 1:D
                nms.Qs[:, d, a, n] .= (s.qs[:, d, a, n]' * nmt.Cmat)'
            end
        end
    end

    for n in 1:N
        for a in 1:A
            for d in 1:D
                nms.Ps[:, d, a, n] .= (s.ps[:, d, a, n]' * nmt.Cmat)'
            end
        end
    end

    nothing
end

function from_normal_modes!(s::State{P,D,A,N}, nms::NormalModeState{P,D,A,N},
                            nmt::NormalModeTransformer) where {P,D,A,N}
    for n in 1:N
        for a in 1:A
            for d in 1:D
                s.qs[:, d, a, n] .= nmt.Cmat * nms.Qs[:, d, a, n]
            end
        end
    end

    for n in 1:N
        for a in 1:A
            for d in 1:D
                s.ps[:, d, a, n] .= nmt.Cmat * nms.Ps[:, d, a, n]
            end
        end
    end

    nothing
end

function in_normal_modes(f, nms::NormalModeState{P,D,A,N}, s::State{P,D,A,N},
                         nmt::NormalModeTransformer) where {P,D,A,N}
    to_normal_modes!(nms, s, nmt)
    f()
    from_normal_modes!(s, nms, nmt)
end

struct OscillatorPropagator
    "Propagation coefficients."
    coefs::Vector{Vector{Float64}}
end

function OscillatorPropagator(dt::Float64, omegas::AbstractVector{Float64})
    coef1 = dt .* sinc.(dt .* omegas ./ pi) # ps
    coef2 = cos.(dt .* omegas)
    coef3 = cos.(dt .* omegas)
    coef4 = -omegas .* sin.(dt .* omegas) # 1 / ps

    coefs = [coef1, coef2, coef3, coef4]

    OscillatorPropagator(coefs)
end

function propagate_oscillators!(nms::NormalModeState{P,D,A,N},
                                oscprop::OscillatorPropagator) where {P,D,A,N}
    for n in 1:N
        for a in 1:A
            m = MS_fict(P)[a] # g / mol

            for d in 1:D
                Qs_p = (oscprop.coefs[1] .* nms.Ps[:, d, a, n] ./ m
                        .+ oscprop.coefs[2] .* nms.Qs[:, d, a, n])
                Ps_p = (oscprop.coefs[3] .* nms.Ps[:, d, a, n]
                        .+ oscprop.coefs[4] .* nms.Qs[:, d, a, n] .* m)

                nms.Qs[:, d, a, n] .= Qs_p
                nms.Ps[:, d, a, n] .= Ps_p
            end
        end
    end

    nothing
end

struct Constrainer
    "Constraint distance (nm)."
    R::Float64
    "Propagation coefficients rotated to normal mode basis."
    coefs_rot::Vector{Matrix{Float64}}
end

function Constrainer(R::Float64, oscprop::OscillatorPropagator,
                     Cmat::Matrix{Float64})
    coefs_rot = [Cmat * diagm(0 => coef) * Cmat' for coef in oscprop.coefs]

    Constrainer(R, coefs_rot)
end

struct ConstraintQuantities
    "Center of mass displacement (nm)."
    disp::Vector{Float64}
    "Center of mass distance (nm)."
    dist::Float64
    "Gradient of center of mass distance."
    grad::Array{Float64,3}
end

function ConstraintQuantities(s::State)
    disp = com_disp(CN1, CN2, CJ, s)
    dist = norm(disp)
    grad = com_r_grad(CN1, CN2, CJ, s)

    ConstraintQuantities(disp, dist, grad)
end

function apply_constraint!(s::State{P,D,A,N}, cstr::Constrainer,
                           cq::ConstraintQuantities) where {P,D,A,N}
    # Reduced mass between the two molecules (g / mol).
    M_red = 0.5 * sum(MS_fict(P))

    ddisp = com_disp(CN1, CN2, CJ, s) - cq.disp # nm
    x = dot(cq.disp ./ cq.dist, ddisp) # nm
    discriminant = cstr.R^2 - (norm(ddisp)^2 - x^2) # nm^2

    discriminant >= 0.0 || error("Cannot enforce constraint.")

    lambda = -(cq.dist + x - sqrt(discriminant)) * M_red # g nm / mol

    for n in 1:N
        for a in 1:A
            for d in 1:D
                # g nm / ps mol
                kp = cq.grad[d, a, n] * lambda / cstr.coefs_rot[1][1, 1]
                kq = kp / MS_fict(P)[a] # nm / ps
                s.qs[:, d, a, n] .+= kq .* cstr.coefs_rot[1][:, 1]
                s.ps[:, d, a, n] .+= kp .* cstr.coefs_rot[3][:, 1]
            end
        end
    end

    nothing
end

function apply_velocity_constraint!(s::State{P,D,A,N}) where {P,D,A,N}
    # Reduced mass between the two molecules (g / mol).
    M_red = 0.5 * sum(MS_fict(P))

    grad = com_r_grad(CN1, CN2, CJ, s)

    lambda = 0.0 # nm / ps

    for n in 1:N
        for d in 1:D
            lambda -= dot(grad[d, :, n], s.ps[CJ, d, :, n] ./ MS_fict(P))
        end
    end

    s.ps[CJ, :, :, :] .+= M_red .* lambda .* grad

    nothing
end

struct ForceApplicator
    "Water model, containing force implementations."
    model::WaterModel
    "Additional force terms."
    extra_forces::Dict{String,Vector{AdHocForce}}
    "Time step (ps)."
    dt::Float64
    "Result of the previous force evaluation."
    current_force::ForceContainer
end

function eval_pot(forceapp::ForceApplicator, s::State{P,D,A,N}) where {P,D,A,N}
    result = Dict{String,Float64}()

    result["model"] = eval_pot(forceapp.model, s) / P

    for (label, terms) in forceapp.extra_forces
        haskey(result, label) && error("Invalid label: '$(label)'.")

        result[label] = sum(eval_pot(term, s) for term in terms)
    end

    result
end

function apply_force!(s::State{P,D,A,N}, forceapp::ForceApplicator;
                      fresh_force::Bool=true) where {P,D,A,N}
    if fresh_force || !forceapp.current_force.initialized
        fill!(forceapp.current_force.f, 0.0)
        forceapp.current_force.initialized = true

        eval_force!(forceapp.current_force, forceapp.model, s)

        forceapp.current_force.f ./= P

        for terms in values(forceapp.extra_forces)
            for term in terms
                eval_force!(forceapp.current_force, term, s)
            end
        end
    end

    s.ps .+= forceapp.dt .* forceapp.current_force.f

    nothing
end

struct Thermostat
    "Momentum dissipation coefficient (P)."
    coef1::Vector{Float64}
    "Momentum fluctuation coefficient (A, P; g nm / ps mol)."
    coef2::Matrix{Float64}
end

function Thermostat(beta::Float64, gamma0::Float64, dt::Float64,
                    omegas::AbstractVector{Float64})
    P = length(omegas)

    frictions = 2omegas # 1 / ps
    frictions[1] = gamma0

    coef1 = exp.(-dt .* frictions)
    coef2 = sqrt.((1.0 .- coef1.^2)' .* MS_fict(P) ./ beta)

    Thermostat(coef1, coef2)
end

function apply_thermostat!(nms::NormalModeState{P,D,A,N},
                           thermostat::Thermostat) where {P,D,A,N}
    for n in 1:N
        for a in 1:A
            for d in 1:D
                nms.Ps[:, d, a, n] .*= thermostat.coef1
                nms.Ps[:, d, a, n] .+= thermostat.coef2[a, :] .* randn(P)
            end
        end
    end

    nothing
end

abstract type Integrator end

const registered_integrators = Dict{String,Type{<:Integrator}}()

list_integrators() = sort(collect(keys(registered_integrators)))

get_integrator(name::String) = registered_integrators[name]

include("integrator_baoab.jl")
include("integrator_cbaoab.jl")
include("integrator_cobabo.jl")
include("integrator_obabo.jl")


text:
# using Routing
include("../src/Routing_include.jl")

# To plot
using PyPlot

function plot_solomon_solution(solomon::Solomon, sol_routes, sol_obj, duration)
    data_name = solomon.data_name
    nodes = solomon.nodes

    n_nodes = length(nodes)
    N = 1:n_nodes

    fig = figure()
    for i in N
        plot(nodes[i].cx, nodes[i].cy, marker=".", color="r", label=i)
        annotate(i, (nodes[i].cx, nodes[i].cy), textcoords="offset points", xytext=(2,2), ha="center")
    end
    for r in sol_routes
        for n in 1:length(r)-1
            i, j = r[n], r[n + 1]
            dx = nodes[j].cx - nodes[i].cx
            dy = nodes[j].cy - nodes[i].cy
            arrow(nodes[i].cx, nodes[i].cy, dx, dy,
                color="b", head_width=0.5, length_includes_head=true)
        end
    end
    sol_obj = round(sol_obj * 100) / 100
    duration = round(duration * 100) /100
    title("$(data_name): obj=$(sol_obj), duration=$(duration) s")
    savefig("$(data_name).png", dpi=1000)
    close(fig)
end

# solomon_vrptw = generate_solomon_vrptw_instance("RC102_025")

solomon_dataset_name = "R102_025"

solomon = load_solomon(solomon_dataset_name)
vrptw = generate_solomon_vrptw_instance(solomon)

start_time = time()
@time routes, obj_val = solveVRP(vrptw, pricing_method="monodirectional");
end_time = time()
duration = end_time - start_time

@show obj_val
@show routes

plot_solomon_solution(solomon, routes, obj_val, duration)
@show solomon_dataset_name, obj_val, duration

include("solomon_exact.jl")
solomon_exact(vrptw)


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

mutable struct LogicManager <: AbstractLogicManager
    transaction::AbstractCoreTransaction
end

"""
    get_rule(log_mgr::AbstractLogicManager, label::AbstractString)

The get_rule function will return the rule for the given label
"""
function get_rule(log_mgr::AbstractLogicManager, label::AbstractString)
    res = execute(log_mgr.transaction, LogicManagerRequestBuilder.get_rule_req(label))
    result = if which_oneof(res.logic_manager_res.get_rule_res, :res) == :rule
                Rule(res.logic_manager_res.get_rule_res.rule)
             else
                nothing
             end
    return result
end
"""
    get_rules(log_mgr::AbstractLogicManager)

Here we get all rules in the schema exposed by the session.
"""
function get_rules(log_mgr::AbstractLogicManager)
    dbResult = stream(log_mgr.transaction, LogicManagerRequestBuilder.get_rules_req())
    answers = [res.logic_manager_res_part.get_rules_res_part.rules for res in dbResult]
    rules = map(Rule,answers)
    result = reduce(vcat, rules)
    return result
end

"""
    put_rule(log_mgr::AbstractLogicManager, label::AbstractString, when::AbstractString, then::AbstractString)

The function gives the possibility to formulate a rule and put it in the database. The when
and then clauses will be written in TypeQL terms.
"""
function put_rule(log_mgr::AbstractLogicManager, label::AbstractString, when::AbstractString, then::AbstractString)
    res = execute(log_mgr.transaction, LogicManagerRequestBuilder.put_rule_req(label, when, then))
    result = Rule(res.logic_manager_res.get_rule_res.rule)
    return result
end


text:
abstract type AbstractGeometry end
# """The general `Geometry` type."""
# struct Geometry <: AbstractGeometry end  # commented out as it won't be used directly

abstract type AbstractGeometryCollection <: AbstractGeometry end
"""A `GeometryCollection` is a collection of `Geometry`s."""
struct GeometryCollection <: AbstractGeometryCollection end

abstract type AbstractPoint <: AbstractGeometry end
"""A simple `Point`."""
struct Point <: AbstractPoint end

abstract type AbstractCurve <: AbstractGeometry end
abstract type AbstractLineString <: AbstractCurve end
"""A `LineString` is a collection of straight lines between its `Point`s."""
struct LineString <: AbstractLineString end
"""A Line is `LineString` with just two points."""
struct Line <: AbstractLineString end
"""A LinearRing is a `LineString` with the same begin and endpoint."""
struct LinearRing <: AbstractLineString end

"""A `CircularString` is a curve, with an odd number of points.
A single segment consists of three points, where the first and last are the beginning and end,
while the second is halfway the curve."""
struct CircularString <: AbstractCurve end
"""A `CompoundCurve` is a curve that combines straight `LineString`s and curved `CircularString`s."""
struct CompoundCurve <: AbstractCurve end

abstract type AbstractSurface <: AbstractGeometry end
abstract type AbstractCurvePolygon <: AbstractSurface end
"""A `Polygon` that can contain either circular or straight curves as rings."""
struct CurvePolygon <: AbstractCurvePolygon end
abstract type AbstractPolygon <: AbstractCurvePolygon end
"""A `Polygon` with straight rings either as exterior or interior(s)."""
struct Polygon <: AbstractPolygon end
"""A `Polygon` with straight rings either as exterior or interior(s)."""
struct Triangle <: AbstractPolygon end

"""A `Polygon` that is rectangular and could be described by the minimum and maximum vertices."""
struct Rectangle <: AbstractPolygon end
"""A `Polygon` with four vertices."""
struct Quad <: AbstractPolygon end
"""A `Polygon` with five vertices."""
struct Pentagon <: AbstractPolygon end
"""A `Polygon` with six vertices."""
struct Hexagon <: AbstractPolygon end

abstract type AbstractPolyHedralSurface <: AbstractSurface end
"""A `PolyHedralSurface` is a connected surface consisting of Polygons."""
struct PolyHedralSurface <: AbstractPolyHedralSurface end
"""A `TIN` is a `PolyHedralSurface` consisting of `Triangle`s."""
struct TIN <: AbstractPolyHedralSurface end  # Surface consisting of Triangles

abstract type AbstractMultiPoint <: AbstractGeometryCollection end
"""A `MultiPoint` is a collection of `Point`s."""
struct MultiPoint <: AbstractMultiPoint end

abstract type AbstractMultiCurve <: AbstractGeometryCollection end
"""A `MultiCurve` is a collection of `Curve`s."""
struct MultiCurve <: AbstractMultiCurve end
abstract type AbstractMultiLineString <: AbstractMultiCurve end
"""A `MultiPoint` is a collection of `Point`s."""
struct MultiLineString <: AbstractMultiLineString end

abstract type AbstractMultiSurface <: AbstractGeometryCollection end
abstract type AbstractMultiPolygon <: AbstractMultiSurface end
"""A `MultiPolygon` is a collection of `Polygon`s."""
struct MultiPolygon <: AbstractMultiPolygon end


text:
using GLMakie

function equidistant_cubed_shell_warp(a, b, c, R = max(abs(a), abs(b), abs(c)))

    r = hypot(a, b, c)

    x1 = R * a / r
    x2 = R * b / r
    x3 = R * c / r

    return x1, x2, x3
end

n = 20

r = range(-1,1,length=n+1)

xa = zeros(n+1,n+1)
xb = zeros(n+1,n+1)
xc = zeros(n+1,n+1)

xa .= r
xb .= r'
xc .= 1

a = [equidistant_cubed_shell_warp(a,b,c)[1] for (a,b,c) in zip(xa,xb,xc)]
b = [equidistant_cubed_shell_warp(a,b,c)[2] for (a,b,c) in zip(xa,xb,xc)]
c = [equidistant_cubed_shell_warp(a,b,c)[3] for (a,b,c) in zip(xa,xb,xc)]

sc = Scene()

wireframe!(sc,
           a,b,c,
           show_axis=false,
           linewidth=1.2)
wireframe!(sc,
           a,b,-c,
           show_axis=false,
           linewidth=1.2)
wireframe!(sc,
           c,a,b,
           show_axis=false,
           linewidth=1.2)
wireframe!(sc,
           -c,a,b,
           show_axis=false,
           linewidth=1.2)
wireframe!(sc,
           b,c,a,
           show_axis=false,
           linewidth=1.2)
wireframe!(sc,
           b,-c,a,
           show_axis=false,
           linewidth=1.2)

sc


text:
using Pkg; Pkg.activate("GraphFusedElasticnet.jl/")
include("GraphFusedElasticNet.jl/src/ARS.jl")
using Plots

μ, σ = 1.0, 2.0
f(x) = exp(-0.5(x - μ)^2 / σ^2) / sqrt(2pi * σ^2)
support = (-10.0, 2.0)

# Build the sampler and simulate 10,000 samples
sampler = RejectionSampler(f, support, (-1.0, 1.0), max_segments = 3)
@time sim = run_sampler!(sampler, 1000);

x = range(-10.0, 10.0, length=100)
envelop = [eval_envelop(sampler.envelop, xi) for xi in x]
target = [f(xi) for xi in x]

histogram(sim, normalize = true, label = "Histogram")
plot!(x, [target envelop], width = 2, label = ["Normal(μ, σ)" "Envelop"])

max_stars_repo_name: mauriciogtec/gfen-reproduce

max_issues_repo_name: mauriciogtec/gfen-reproduce

max_forks_repo_name: mauriciogtec/gfen-reproduce


text:
include("routing.jl")
include("design/design.jl")


text:
"""
    struct CoreModel <: MetabolicModel

A "bare bones" core linear optimization problem of the form, with reaction and
metabolite names.
```
min c^T x
s.t. S x = b
      xₗ ≤ x ≤ xᵤ
```
"""
mutable struct CoreModel <: MetabolicModel
    S::SparseMat
    b::SparseVec
    c::SparseVec
    xl::Vector{Float64}
    xu::Vector{Float64}
    rxns::Vector{String}
    mets::Vector{String}
    grrs::Vector{Maybe{GeneAssociation}}

    function CoreModel(
        S::MatType,
        b::VecType,
        c::VecType,
        xl::VecType,
        xu::VecType,
        rxns::StringVecType,
        mets::StringVecType,
        grrs::Vector{Maybe{GeneAssociation}} = Vector{Maybe{GeneAssociation}}(
            nothing,
            length(rxns),
        ),
    )
        all([length(b), length(mets)] .== size(S, 1)) ||
            throw(DimensionMismatch("inconsistent number of metabolites"))

        all(
            [length(c), length(xl), length(xu), length(rxns), length(grrs)] .== size(S, 2),
        ) || throw(DimensionMismatch("inconsistent number of reactions"))

        new(sparse(S), sparse(b), sparse(c), collect(xl), collect(xu), rxns, mets, grrs)
    end
end

"""
    reactions(a::CoreModel)::Vector{String}

Get the reactions in a `CoreModel`.
"""
reactions(a::CoreModel)::Vector{String} = a.rxns

"""
    metabolites(a::CoreModel)::Vector{String}

Metabolites in a `CoreModel`.
"""
metabolites(a::CoreModel)::Vector{String} = a.mets

"""
    stoichiometry(a::CoreModel)::SparseMat

`CoreModel` stoichiometry matrix.
"""
stoichiometry(a::CoreModel)::SparseMat = a.S

"""
    bounds(a::CoreModel)::Tuple{Vector{Float64},Vector{Float64}}

`CoreModel` flux bounds.
"""
bounds(a::CoreModel)::Tuple{Vector{Float64},Vector{Float64}} = (a.xl, a.xu)

"""
    balance(a::CoreModel)::SparseVec

`CoreModel` target flux balance.
"""
balance(a::CoreModel)::SparseVec = a.b

"""
    objective(a::CoreModel)::SparseVec

`CoreModel` objective vector.
"""
objective(a::CoreModel)::SparseVec = a.c

"""
    reaction_stoichiometry(model::CoreModel, rid::String)::Dict{String, Float64}

Return the stoichiometry of reaction with ID `rid`.
"""
reaction_stoichiometry(m::CoreModel, rid::String)::Dict{String,Float64} =
    Dict(m.mets[k] => v for (k, v) in zip(findnz(m.S[:, first(indexin([rid], m.rxns))])...))

"""
    reaction_stoichiometry(model::CoreModel, ridx)::Dict{String, Float64}

Return the stoichiometry of reaction at index `ridx`.
"""
reaction_stoichiometry(m::CoreModel, ridx)::Dict{String,Float64} =
    Dict(m.mets[k] => v for (k, v) in zip(findnz(m.S[:, ridx])...))

"""
    reaction_gene_association_vec(model::CoreModel)::Vector{Maybe{GeneAssociation}}

Retrieve a vector of all gene associations in a [`CoreModel`](@ref), in the
same order as `reactions(model)`.
"""
reaction_gene_association_vec(model::CoreModel)::Vector{Maybe{GeneAssociation}} = model.grrs

"""
    reaction_gene_association(model::CoreModel, ridx::Int)::Maybe{GeneAssociation}

Retrieve the [`GeneAssociation`](@ref) from [`CoreModel`](@ref) by reaction
index.
"""
reaction_gene_association(model::CoreModel, ridx::Int)::Maybe{GeneAssociation} =
    model.grrs[ridx]

"""
    reaction_gene_association(model::CoreModel, rid::String)::Maybe{GeneAssociation}

Retrieve the [`GeneAssociation`](@ref) from [`CoreModel`](@ref) by reaction ID.
"""
reaction_gene_association(model::CoreModel, rid::String)::Maybe{GeneAssociation} =
    model.grrs[first(indexin([rid], model.rxns))]

"""
    Base.convert(::Type{CoreModel}, m::M) where {M <: MetabolicModel}

Make a `CoreModel` out of any compatible model type.
"""
function Base.convert(::Type{CoreModel}, m::M) where {M<:MetabolicModel}
    if typeof(m) == CoreModel
        return m
    end

    (xl, xu) = bounds(m)
    CoreModel(
        stoichiometry(m),
        balance(m),
        objective(m),
        xl,
        xu,
        reactions(m),
        metabolites(m),
        Vector{Maybe{GeneAssociation}}([
            reaction_gene_association(m, id) for id in reactions(m)
        ]),
    )
end


text:
using StarWarsArrays
using Test

import StarWarsArrays: StarWarsError, order

@testset "OriginalOrder" begin
    v = StarWarsArray(collect(1:9))
    @test order(v) == OriginalOrder
    @test size(v) == size(parent(v))
    @test v[5] == 2
    @test v[3] == 6
    @test v[7] == 7
    @test_throws BoundsError v[0]
    @test join(split(repr("text/plain", v), '\n')[begin+1:end], '\n') ==
        " 4\n 5\n 6\n 1\n 2\n 3\n 7\n 8\n 9"
    v[4] = -42
    @test v[4] == -42
    @test axes(v) == ([4, 5, 6, 1, 2, 3, 7, 8, 9],)
    @test eachindex(v) == [4, 5, 6, 1, 2, 3, 7, 8, 9]
    @test length(v) == 9
    @test size(v) == (9,)
    for (idx, x) in enumerate(v)
        # Explanation: iterating over `v` should be equivalent to iterating over its parent,
        # because the iteration protocol should follow the order 4, 5, 6, 1, 2, 3, ...
        @test x == parent(v)[idx]
    end

    m = StarWarsArray(reshape(collect(1:81), 9, 9), OriginalOrder)
    @test order(m) == OriginalOrder
    @test size(m) == size(parent(m))
    @test m[5] == 2
    @test m[3,4] == 6
    @test m[7,2] == 43
    @test m[4,4] == 1
    @test_throws BoundsError m[0, 12]
    @test join(split(repr("text/plain", m), '\n')[begin+1:begin+3], '\n') ==
        " 31  40  49  4  13  22  58  67  76\n 32  41  50  5  14  23  59  68  77\n 33  42  51  6  15  24  60  69  78"
    m[4,4] = -42
    @test m[4] == -42
    m[1] = 0
    @test m[1,4] == 0
    @test axes(m) == ([4, 5, 6, 1, 2, 3, 7, 8, 9], [4, 5, 6, 1, 2, 3, 7, 8, 9])
    @test length(m) == 81
    @test size(m) == (9, 9)
    for (idx, x) in enumerate(m)
        @test x == parent(m)[idx]
    end

    @test_throws ArgumentError StarWarsArray(rand(2), OriginalOrder)
    @test_throws ArgumentError StarWarsArray(rand(8, 2), OriginalOrder)
end

@testset "MacheteOrder" begin
    v = StarWarsArray(collect(1:9), MacheteOrder)
    @test order(v) == MacheteOrder
    @test size(v) == (8,)
    @test v[5] == 2
    @test v[3] == 4
    @test v[6] == 5
    @test v[7] == 6
    @test_throws StarWarsError v[1]
    @test_throws BoundsError v[0]
    @test join(split(repr("text/plain", v), '\n')[begin+1:end], '\n') ==
        " 3\n 4\n 1\n 2\n 5\n 6\n 7\n 8"
    v[4] = -42
    @test v[4] == -42
    @test axes(v) == ([3, 4, 1, 2, 5, 6, 7, 8],)
    @test eachindex(v) == [3, 4, 1, 2, 5, 6, 7, 8]
    @test length(v) == 8
    @test size(v) == (8,)
    for (idx, x) in enumerate(v)
        @test x == parent(v)[idx]
    end

    m = StarWarsArray(reshape(collect(1:81), 9, 9), MacheteOrder)
    @test order(m) == MacheteOrder
    @test size(m) == (8, 8)
    @test m[5] == 2
    @test m[3,4] == 4
    @test m[7,2] == 24
    @test m[4,4] == 1
    @test m[9,4] == 8
    @test_throws BoundsError m[0, 12]
    @test_throws BoundsError m[82]
    @test_throws StarWarsError m[1]
    @test_throws StarWarsError m[3, 1]
    rpr = split(repr("text/plain", m), '\n')
    @test rpr[2] == " 21  30  3  12  39  48  57  66"
    @test rpr[6] == " 23  32  5  14  41  50  59  68"
    m[4,4] = -42
    @test m[4] == -42
    m[2] = 0
    @test m[2,4] == 0
    @test axes(m) == ([3, 4, 1, 2, 5, 6, 7, 8], [3, 4, 1, 2, 5, 6, 7, 8])
    @test length(m) == 64
    @test size(m) == (8, 8)
    for (idx, x) in enumerate(m)
        @test x == parent(m)[idx]
    end

    @test_throws ArgumentError StarWarsArray(rand(4), MacheteOrder)
    @test_throws ArgumentError StarWarsArray(rand(6, 2, 9), MacheteOrder)

    @test sprint(showerror, StarWarsArrays.StarWarsError(1, MacheteOrder)) ==
        "StarWarsError: there is no episode 1 in MacheteOrder"
end

max_stars_repo_name: giordano/StarWarsArrays

max_issues_repo_name: giordano/StarWarsArrays

max_forks_repo_name: giordano/StarWarsArrays


text:
module RocketTeardownTest

using Test

import Rocket
import Rocket: TeardownLogic, UnsubscribableTeardownLogic, CallableTeardownLogic, VoidTeardownLogic, InvalidTeardownLogic
import Rocket: Teardown, as_teardown
import Rocket: unsubscribe!, teardown!, on_unsubscribe!

import Rocket: InvalidTeardownLogicTraitUsageError, MissingOnUnsubscribeImplementationError

@testset "Teardown" begin

    struct DummyType end

    struct AnotherDummyType end
    Rocket.as_teardown(::Type{<:AnotherDummyType}) = VoidTeardownLogic()

    struct DummySubscription end
    Rocket.as_teardown(::Type{<:DummySubscription}) = UnsubscribableTeardownLogic()

    struct ImplementedSubscription end
    Rocket.as_teardown(::Type{<:ImplementedSubscription}) = UnsubscribableTeardownLogic()
    Rocket.on_unsubscribe!(::ImplementedSubscription)     = "unsubscribed"

    @testset "as_teardown" begin
        # Check if arbitrary dummy type has undefined teardown logic
        @test as_teardown(DummyType) === InvalidTeardownLogic()

        # Check if as_teardown returns specified teardown logic
        @test as_teardown(AnotherDummyType) === VoidTeardownLogic()
        @test as_teardown(DummySubscription) === UnsubscribableTeardownLogic()
        @test as_teardown(ImplementedSubscription) === UnsubscribableTeardownLogic()

        # Check if as_teardown returns CallableTeardownLogic for Function object
        @test as_teardown(Function) === CallableTeardownLogic()
    end

    @testset "unsubscribe!" begin
        # Check if arbitrary dummy type throws an error in unsubscribe!
        @test_throws InvalidTeardownLogicTraitUsageError unsubscribe!(DummyType())

        # Check if void teardown object does nothing
        @test unsubscribe!(AnotherDummyType()) === nothing

        # Check if function object calls itself in unsubscribe!
        @test unsubscribe!(() -> return 1) === 1

        #Check if dummy subscription throws an error in unusubscribe!
        @test_throws MissingOnUnsubscribeImplementationError unsubscribe!(DummySubscription())

        #Check if implemented subscription calls on_unsubscribe!
        @test unsubscribe!(ImplementedSubscription()) === "unsubscribed"
    end

end

end


text:
module BSplines

using RecipesBase

include("knot_sets.jl")
include("quadrature.jl")
include("basis.jl")
include("splines.jl")

end