text:
#Test for Gadfly with different chunk options and figure formatsusing Weave
using Weave, Compat
using Base.Test

weave("documents/gadfly_formats_test.txt", "tex")
result = @compat readstring(open("documents/gadfly_formats_test.tex"))
ref = @compat readstring(open("documents/gadfly_formats_test_ref.tex"))
@test result == ref

weave("documents/gadfly_formats_test.txt", doctype="tex", fig_ext=".tex", plotlib="gadfly")
result = @compat readstring(open("documents/gadfly_formats_test.tex"))
ref = @compat readstring(open("documents/gadfly_formats_test_tikz_ref.tex"))
@test result == ref

weave("documents/gadfly_formats_test.txt", doctype="tex", fig_ext=".ps", plotlib="gadfly")
result = @compat readstring(open("documents/gadfly_formats_test.tex"))
ref = @compat readstring(open("documents/gadfly_formats_test_ps_ref.tex"))
@test result == ref

weave("documents/gadfly_formats_test.txt", doctype="pandoc", plotlib="gadfly")
result = @compat readstring(open("documents/gadfly_formats_test.md"))
ref = @compat readstring(open("documents/gadfly_formats_test_pandoc_ref.md"))
@test result == ref

weave("documents/gadfly_formats_test.txt", doctype="pandoc", plotlib="gadfly", fig_ext=".svg")
result = @compat readstring(open("documents/gadfly_formats_test.md"))
ref = @compat readstring(open("documents/gadfly_formats_test_svg_ref.md"))
@test result == ref

weave("documents/gadfly_formats_test.txt", doctype="github", plotlib="gadfly", fig_ext=".js.svg")
result = @compat readstring(open("documents/gadfly_formats_test.md"))
ref = @compat readstring(open("documents/gadfly_formats_test_jssvg_ref.md"))
@test result == ref

weave("documents/gadfly_formats_test.txt", doctype="rst", plotlib="gadfly")
result = @compat readstring(open("documents/gadfly_formats_test.rst"))
ref = @compat readstring(open("documents/gadfly_formats_test_ref.rst"))
@test result == ref

weave("documents/gadfly_formats_test.txt", doctype="multimarkdown", plotlib="gadfly")
result = @compat readstring(open("documents/gadfly_formats_test.md"))
ref = @compat readstring(open("documents/gadfly_formats_test_mmd_ref.md"))
@test result == ref

weave("documents/gadfly_formats_test.txt", doctype="asciidoc", plotlib="gadfly",
    out_path="documents/output")
result = @compat readstring(open("documents/output/gadfly_formats_test.txt"))
ref = @compat readstring(open("documents/output/gadfly_formats_test_ref.txt"))
@test result == ref

weave("documents/gadfly_markdown_test.jmd", doctype="github",plotlib="gadfly", informat="markdown")
result = @compat readstring(open("documents/gadfly_markdown_test.md"))
ref = @compat readstring(open("documents/gadfly_markdown_test_ref.md"))
@test result == ref

weave("documents/FIR_design.jl", doctype="pandoc", plotlib="gadfly", informat="script")
result = @compat readstring(open("documents/FIR_design.md"))
ref = @compat readstring(open("documents/FIR_design_ref.md"))
@test result == ref


text:
function _print_welcome_message(a::AbstractDict = ENV)::Nothing
    predictmdextra_version::VersionNumber = version()
    predictmdextra_pkgdir::String = package_directory()
    if !_suppress_welcome_message(a)
        @info(string("This is PredictMDExtra, version ",predictmdextra_version,),)
        @info(string("For help, please visit https://predictmd.net",),)
        @debug(string("PredictMDExtra package directory: ",predictmdextra_pkgdir,),)
    end
    return nothing
end

function _suppress_welcome_message(a::AbstractDict = ENV)::Bool
    return get(a, "SUPPRESS_PREDICTMDEXTRA_WELCOME_MESSAGE", "false") == "true"
end

ext: jl


text:
using OnlineStatsBase, LearnBase, Test

import OnlineStatsBase: _fit!, _merge!

#-----------------------------------------------------------------------# Stat for testing
mutable struct Counter{T} <: OnlineStat{T}
    n::Int
    Counter{T}() where {T} = new{T}(0)
end
Counter(T = Number) = Counter{T}()
_fit!(o::Counter{T}, y) where {T} = (o.n += 1)
_merge!(a::Counter{T}, b::Counter{T}) where {T} = (a.n += b.n)

@testset "abstract methods" begin
    o = Counter()
    println(o)
    @test value(o) == 0
    @test nobs(o) == 0
    fit!(o, rand())
    @test value(o) == nobs(o) == 1
    @test o == fit!(Counter(), rand())
    @test value(fit!(Counter(), rand(10))) == 10
    @test merge(fit!(Counter(), 1:5), fit!(Counter(), 1:5)) == fit!(Counter(), 1:10)
end

#-----------------------------------------------------------------------# Weight
@testset "Weight" begin
function test_weight(w::OnlineStatsBase.Weight, f::Function)
    println("  > $w")
    @test w == copy(w)
    for i in 1:20
        @test w(i) == f(i)
    end
end
test_weight(@inferred(EqualWeight()),               i -> 1 / i)
test_weight(@inferred(ExponentialWeight(.1)),       i -> ifelse(i==1, 1.0, .1))
test_weight(@inferred(LearningRate(.6)),            i -> 1 / i^.6)
test_weight(@inferred(LearningRate2(.5)),           i -> 1 / (1 + .5*(i-1)))
test_weight(@inferred(HarmonicWeight(4.)),          i -> 4 / (4 + i - 1))
test_weight(@inferred(max(EqualWeight(), .1)),      i -> max(.1, 1 / i))
test_weight(@inferred(max(LearningRate(.6), .1)),   i -> max(.1, 1 / i^.6))
test_weight(@inferred(.1 * EqualWeight()),          i -> .1 * (1 / i))
test_weight(@inferred(.1 * EqualWeight()),          i -> .1 * (1 / i))
test_weight(max(.5 * EqualWeight(), .1),            i -> max(.1, .5 / i))

@test ExponentialWeight(20) == ExponentialWeight(2 / 21)
@test max(.1, EqualWeight()) == max(EqualWeight(), .1)
@test .1 * EqualWeight() == EqualWeight() * .1
@testset "McclainWeight" begin
    w = McclainWeight(.1)
    for i in 2:100
        @test .1 < w(i) < 1
    end
end
@testset "first weight is one" begin
    for w in [EqualWeight(), ExponentialWeight(), LearningRate(), LearningRate2(),
              HarmonicWeight(), McclainWeight()]
        @test w(1) == 1
    end
end
end #Weight

#-----------------------------------------------------------------------# Iteration
@testset "Iteration" begin
    x, y = randn(100,10), randn(100)
    @test fit!(Counter(Vector), eachrow(x)).n == 100
    @test fit!(Counter(Vector), eachcol(x)).n == 10
    @test fit!(Counter(Tuple),  eachrow(x,y)).n == 100
    @test fit!(Counter(Tuple),  eachcol(x,y)).n == 10

    for (j, xj) in enumerate(eachcol(x))
        @test xj == x[:, j]
    end
    for (i, xi) in enumerate(eachrow(x))
        @test xi == x[i, :]
    end
    @inferred eachrow(x)
    @inferred eachcol(x)
    @inferred eachrow(x, y)
    @inferred eachcol(x, y)

    @test length(eachcol(x)) == 10
    @test length(eachrow(x)) == 100
    @test length(eachrow(x, y)) == 100
end


text:
# Packages that throw "not installed" error on upgrade: Compat, DSP, HDF5
using Pkg
function pkg_check(pkgs::Array{String,1})

  # why would you deprecate something this useful?
  if VERSION >= v"1.4"
    installs = Dict{String, VersionNumber}()
    for (uuid, dep) in Pkg.dependencies()
        dep.is_direct_dep || continue
        dep.version === nothing && continue
        installs[dep.name] = dep.version
    end
  else
    installs = Pkg.installed()
  end

  for p in pkgs
    if get(installs, p, nothing) == nothing
      @warn(string(p * " not found! Installing."))
      Pkg.add(p)
    else
      println(p * " found. Not installing.")
    end
  end
  return nothing
end
pkg_check(["DSP", "HDF5"])


text:
module LinearSystemsBench

import Base.A_ldiv_B!, Base.\

using BenchmarkTools
using IterativeSolvers

include("advection_diffusion.jl")

# A DiagonalMatrix that doesn't check whether it is singular in the \ op.
struct DiagonalPreconditioner{T}
    diag::Vector{T}
end

function A_ldiv_B!(y::AbstractVector{T}, A::DiagonalPreconditioner{T}, b::AbstractVector{T}) where T
    for i = 1 : length(b)
        @inbounds y[i] = A.diag[i] \ b[i]
    end
    y
end

(\)(D::DiagonalPreconditioner, b::AbstractVector) = D.diag .\ b

function posdef(n)
    A = SymTridiagonal(fill(2.01, n), fill(-1.0, n))
    b = A * ones(n)
    return A, b
end

function cg(; n = 1_000_000, tol = 1e-6, maxiter::Int = 200)
    A, b = posdef(n)
    P = DiagonalPreconditioner(collect(linspace(1.0, 2.0, n)))

    println("Symmetric positive definite matrix of size ", n)
    println("Eigenvalues in interval [0.01, 4.01]")
    println("Tolerance = ", tol, "; max #iterations = ", maxiter)

    # Dry run
    initial = rand(n)
    IterativeSolvers.cg!(copy(initial), A, b, Pl = P, maxiter = maxiter, tol = tol, log = false)

    # Actual benchmark
    @benchmark IterativeSolvers.cg!(x0, $A, $b, Pl = $P, maxiter = $maxiter, tol = $tol, log = false) setup=(x0 = copy($initial))
end

function indefinite(n)
    # Generate an indefinite "hard" matrix
    srand(1)
    A = speye(n) + sprand(n, n, 5.0 / n)
    A = (A + A') / 2
    x = ones(n)
    b = A * x

    A, b
end

function gmres(; n = 100_000, tol = 1e-5, restart::Int = 15, maxiter::Int = 210)
    A, b = indefinite(n)

    println("Matrix of size ", n, " with ~", nnz(A) / n, " nonzeros per row")
    println("Tolerance = ", tol, "; restart = ", restart, "; max #iterations = ", maxiter)

    IterativeSolvers.gmres(A, b, tol = tol, restart = restart, maxiter = maxiter, verbose = true)

    @benchmark IterativeSolvers.gmres($A, $b, tol = $tol, restart = $restart, maxiter = $maxiter)
end

function bicgstabl()
    A, b = advection_dominated()

    b1 = @benchmark IterativeSolvers.bicgstabl($A, $b, 2, max_mv_products = 1000, convex_combination = false) setup = (srand(1))
    b2 = @benchmark IterativeSolvers.bicgstabl($A, $b, 2, max_mv_products = 1000, convex_combination = true) setup = (srand(1))
    b3 = @benchmark IterativeSolvers.bicgstabl($A, $b, 4, max_mv_products = 1000, convex_combination = false) setup = (srand(1))
    b4 = @benchmark IterativeSolvers.bicgstabl($A, $b, 4, max_mv_products = 1000, convex_combination = true) setup = (srand(1))

    b1, b2, b3, b4
end

function minres(n = 100_000)
    A = SymTridiagonal(fill(2.1, n), fill(-1.0, n))
    x = ones(n)
    b = A * x

    @benchmark IterativeSolvers.minres($A, $b, maxiter = 100)
end

function sparse_stationary(n = 10_000)
    A = sprand(n, n, 5 / n)
    A += A' + 4.0I
    x = ones(n)
    b = A * x
    ω = 1.1

    b1 = @benchmark $sor!(sol, $A, $b, $ω, maxiter = 20) setup = (sol = zeros($n))
    b2 = @benchmark $ssor!(sol, $A, $b, $ω, maxiter = 20) setup = (sol = zeros($n))
    b3 = @benchmark $jacobi!(sol, $A, $b, maxiter = 20) setup = (sol = zeros($n))
    b4 = @benchmark $gauss_seidel!(sol, $A, $b, maxiter = 20) setup = (sol = zeros($n))

    b1, b2, b3, b4
end

end


text:
"""
`LinearKernel([ρ=1.0,[c=0.0]])`
The linear kernel is a Mercer kernel given by
```
    κ(x,y) = ρ²xᵀy + c
```
Where `c` is a real number
"""
struct LinearKernel{Tr, Tc<:Real} <: Kernel{Tr}
    transform::Tr
    c::Tc
end

function LinearKernel(ρ::T=1.0, c::Real=zero(T)) where {T<:Real}
    LinearKernel(ScaleTransform(ρ), c)
end

function LinearKernel(ρ::AbstractVector{T}, c::Real=zero(T)) where {T<:Real}
    LinearKernel(ARDTransform(ρ), c)
end

LinearKernel(t::Transform) = LinearKernel(t, 0.0)

params(k::LinearKernel) = (params(transform(k)),k.c)
opt_params(k::LinearKernel) = (opt_params(transform(k)),k.c)

@inline kappa(κ::LinearKernel, xᵀy::T) where {T<:Real} = xᵀy + κ.c

metric(::LinearKernel) = DotProduct()

"""
`PolynomialKernel([ρ=1.0[,d=2.0[,c=0.0]]])`
The polynomial kernel is a Mercer kernel given by
```
    κ(x,y) = (ρ²xᵀy + c)^d
```
Where `c` is a real number, and `d` is a shape parameter bigger than 1
"""
struct PolynomialKernel{Tr,Tc<:Real,Td<:Real} <: Kernel{Tr}
    transform::Tr
    d::Td
    c::Tc
    function PolynomialKernel{Tr, Tc, Td}(transform::Tr, d::Td, c::Tc) where {Tr<:Transform, Td<:Real, Tc<:Real}
        @check_args(PolynomialKernel, d, d >= one(Td), "d >= 1")
        return new{Tr, Td, Tc}(transform,d, c)
    end
end

function PolynomialKernel(ρ::Real=1.0, d::Td=2.0, c::Real=zero(Td)) where {Td<:Real}
    PolynomialKernel(ScaleTransform(ρ), d, c)
end

function PolynomialKernel(ρ::AbstractVector{T}, d::Real=2.0, c::Real=zero(T₁)) where {T<:Real}
    PolynomialKernel(ARDTransform(ρ), d, c)
end

function PolynomialKernel(t::Tr, d::Td=2.0, c::Tc=zero(eltype(Td))) where {Tr<:Transform, Td<:Real, Tc<:Real}
    PolynomialKernel{Tr, Tc, Td}(t, d, c)
end

params(k::PolynomialKernel) = (params(transform(k)),k.d,k.c)
opt_params(k::PolynomialKernel) = (opt_params(transform(k)),k.d,k.c)

@inline kappa(κ::PolynomialKernel, xᵀy::T) where {T<:Real} = (xᵀy + κ.c)^(κ.d)

metric(::PolynomialKernel) = DotProduct()


text:
include("core/environment.jl")
include("core/agent.jl")
include("core/gvf.jl")

max_stars_repo_name: mkschleg/JuliaRL

max_issues_repo_name: mkschleg/JuliaRL

max_forks_repo_name: mkschleg/JuliaRL


text:
function BicycleCar(scenario=:parallel_park, ;N=101)
    # Model
    model = RobotZoo.BicycleModel()
    n,m = size(model)

    if scenario == :parallel_park
        opts = SolverOptions(
            penalty_initial = 1e4,
            verbose = 1,
            cost_tolerance_intermediate = 1e-1
        )

        # Scenario
        tf = 3.0
        dt = tf/(N-1)
        x0 = SA_F64[0,0,0,0]
        xf = SA[0,2,deg2rad(0),0]

        # Objective
        Q = Diagonal(SA[1,1,1e-2,1e-2])
        R = Diagonal(SA[1e0,1e0])
        Qf = Diagonal(SA_F64[1,1,1,10])
        obj = LQRObjective(Q*dt,R*dt,Qf,xf,N)

        # Constraints
        cons = ConstraintList(n,m,N)
        bnd = [10,6.0,Inf,deg2rad(45)]
        bnd = BoundConstraint(n, m, x_min=-bnd, x_max=bnd)

        add_constraint!(cons, bnd, 1:N-1)
        add_constraint!(cons, GoalConstraint(xf), N)

        # Problem
        prob = Problem(model, obj, x0, tf, xf=xf, constraints=cons)
        initial_controls!(prob, SA[-0.1,0.0])
        rollout!(prob)

        return prob, opts
    elseif scenario == :turn90
        opts = SolverOptions(
            penalty_initial=1e3,
        )

        # Scenario
        tf = 3.0
        x0 = SA_F64[0,0,0,0]
        xf = SA[0,0,deg2rad(90),0]

        # Objective
        Q = Diagonal(SA[1,1,1e-2,1e-2])
        R = Diagonal(SA[1e0,1e0])
        Qf = Diagonal(SA_F64[1,1,1,1])*100
        obj = LQRObjective(Q*dt,R*dt,Qf,xf,N)

        # Constraints
        cons = ConstraintList(n,m,N)
        bnd = [1.5,1.5,Inf,deg2rad(45)]
        bnd = BoundConstraint(n, m, x_min=-bnd, x_max=bnd)

        add_constraint!(cons, bnd, 1:N-1)
        add_constraint!(cons, GoalConstraint(xf), N)

        # Problem
        prob = Problem(model, obj, x0, tf, xf=xf, constraints=cons)
        initial_controls!(prob, SA[-0.1,0.0])
        rollout!(prob)

        return prob, opts
    end
    throw(ErrorException("$scenario not a known scenario"))
end


text:
@graphs begin
	@node phase_5 = Start() -> ring_check

	@node ring_check = BinaryCondition("The Shadow has an Elven ring.") -> [n_true=ring_available, n_false=ring_not_available]
	@node ring_available = SetRingAvailable(true) -> modt_check
	@node ring_not_available = SetRingAvailable(false) -> modt_check

	@node modt_check = BinaryCondition("""
									   The Mouth of Sauron is recruited and his "Messenger of the Dark Tower" ability have not been used this turn.
									   """) -> [n_true=modt_available, n_false=modt_not_available]
	@node modt_available = SetMoDTAvailable(true) -> p5_strat
	@node modt_not_available = SetMoDTAvailable(false) -> p5_strat

	@node p5_strat = CheckStrategy("military") -> [n_true=p5_mili, n_false=p5_corr]

	@node p5_mili = JumpToGraph("select_action_mili") -> p5_end
	@node p5_corr = JumpToGraph("select_action_corr") -> p5_end

	@node p5_end = End("End of Phase") -> []
end

max_stars_repo_name: mvmorin/queller-bot

max_issues_repo_name: mvmorin/queller-bot

max_forks_repo_name: mvmorin/queller-bot


text:
using Turing, ContinuousBenchmarks

include(joinpath(ContinuousBenchmarks.STAN_DATA_DIR, "normal-mixture-stan.data.jl"))
include(joinpath(ContinuousBenchmarks.STAN_MODELS_DIR, "normal-mixture.model.jl"))
#include(joinpath(ContinuousBenchmarks.BENCH_DIR, "normal-mixture-stan.run.jl"))

# NOTE: I only run a sub-set of the data as running the whole is quite slow
@tbenchmark(Gibbs(10, HMC(1, 0.05, 1, :theta), PG(50, 1, :k), HMC(1, 0.2, 3, :mu)), nmmodel, simplenormalmixturestandata[1]["y"][1:100])
bench_res = @tbenchmark(Gibbs(1000, HMC(1, 0.05, 1, :theta), PG(50, 1, :k), HMC(1, 0.2, 3, :mu)), nmmodel, simplenormalmixturestandata[1]["y"][1:100])
logd = build_logd("Simple Gaussian Mixture Model", bench_res...)

#logd["stan"] = Dict("theta" => mean(nm_theta), "mu[1]" => mean(nm_mu_1), "mu[2]" =>mean(nm_mu_2))
#logd["time_stan"] = nm_time

print_log(logd, ["theta", "mu[1]", "mu[2]"])
send_log(logd)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct TeamApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_teamAccessLogs(_api::TeamApi; before=nothing, count=nothing, page=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/team.accessLogs", ["slackAuth"])
    Swagger.set_param(_ctx.query, "before", before)  # type String
    Swagger.set_param(_ctx.query, "count", count)  # type String
    Swagger.set_param(_ctx.query, "page", page)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Gets the access logs for the current team.
Param: before::String
Param: count::String
Param: page::String
Return: Dict{String, Any}
"""
function teamAccessLogs(_api::TeamApi; before=nothing, count=nothing, page=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_teamAccessLogs(_api; before=before, count=count, page=page, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function teamAccessLogs(_api::TeamApi, response_stream::Channel; before=nothing, count=nothing, page=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_teamAccessLogs(_api; before=before, count=count, page=page, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_teamBillableInfo(_api::TeamApi; user=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/team.billableInfo", ["slackAuth"])
    Swagger.set_param(_ctx.query, "user", user)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Gets billable users information for the current team.
Param: user::String
Return: Dict{String, Any}
"""
function teamBillableInfo(_api::TeamApi; user=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_teamBillableInfo(_api; user=user, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function teamBillableInfo(_api::TeamApi, response_stream::Channel; user=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_teamBillableInfo(_api; user=user, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_teamInfo(_api::TeamApi; team=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/team.info", ["slackAuth"])
    Swagger.set_param(_ctx.query, "team", team)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Gets information about the current team.
Param: team::String
Return: Dict{String, Any}
"""
function teamInfo(_api::TeamApi; team=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_teamInfo(_api; team=team, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function teamInfo(_api::TeamApi, response_stream::Channel; team=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_teamInfo(_api; team=team, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_teamIntegrationLogs(_api::TeamApi; app_id=nothing, change_type=nothing, count=nothing, page=nothing, service_id=nothing, user=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/team.integrationLogs", ["slackAuth"])
    Swagger.set_param(_ctx.query, "app_id", app_id)  # type String
    Swagger.set_param(_ctx.query, "change_type", change_type)  # type String
    Swagger.set_param(_ctx.query, "count", count)  # type String
    Swagger.set_param(_ctx.query, "page", page)  # type String
    Swagger.set_param(_ctx.query, "service_id", service_id)  # type String
    Swagger.set_param(_ctx.query, "user", user)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Gets the integration logs for the current team.
Param: app_id::String
Param: change_type::String
Param: count::String
Param: page::String
Param: service_id::String
Param: user::String
Return: Dict{String, Any}
"""
function teamIntegrationLogs(_api::TeamApi; app_id=nothing, change_type=nothing, count=nothing, page=nothing, service_id=nothing, user=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_teamIntegrationLogs(_api; app_id=app_id, change_type=change_type, count=count, page=page, service_id=service_id, user=user, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function teamIntegrationLogs(_api::TeamApi, response_stream::Channel; app_id=nothing, change_type=nothing, count=nothing, page=nothing, service_id=nothing, user=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_teamIntegrationLogs(_api; app_id=app_id, change_type=change_type, count=count, page=page, service_id=service_id, user=user, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_teamProfileGet(_api::TeamApi; visibility=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/team.profile.get", ["slackAuth"])
    Swagger.set_param(_ctx.query, "visibility", visibility)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Retrieve a team's profile.
Param: visibility::String
Return: Dict{String, Any}
"""
function teamProfileGet(_api::TeamApi; visibility=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_teamProfileGet(_api; visibility=visibility, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function teamProfileGet(_api::TeamApi, response_stream::Channel; visibility=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_teamProfileGet(_api; visibility=visibility, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export teamAccessLogs, teamBillableInfo, teamInfo, teamIntegrationLogs, teamProfileGet


text:
using Test
using DFTK: lda_x!, lda_c_vwn!
using Libxc

@testset "Fallback lda_x!" begin
    func = Libxc.Functional(:lda_x)

    # Create reference
    ρ = abs.(randn(100))
    Eref = similar(ρ)
    Vref = similar(ρ)
    Libxc.evaluate!(func, rho=ρ, zk=Eref, vrho=Vref)

    # Compute in DFTK in elevated precision
    ρbig = Array{BigFloat}(ρ)
    E=similar(ρbig)
    V=similar(ρbig)
    lda_x!(ρbig, E=E, Vρ=V)

    @test E ≈ Eref atol=5e-15
    @test V ≈ Vref atol=5e-15
end

@testset "Fallback lda_c_vwn!" begin
    func = Libxc.Functional(:lda_c_vwn)

    # Create reference
    ρ = abs.(randn(100))
    Eref = similar(ρ)
    Vref = similar(ρ)
    Libxc.evaluate!(func, rho=ρ, zk=Eref, vrho=Vref)

    # Compute in DFTK in elevated precision
    ρbig = Array{BigFloat}(ρ)
    E=similar(ρbig)
    V=similar(ρbig)
    lda_c_vwn!(ρbig, E=E, Vρ=V)

    @test E ≈ Eref atol=5e-15
    @test V ≈ Vref atol=5e-15
end


text:
"""
Run this script (from this directory) as

    time julia --code-coverage=tracefile-%p.info --code-coverage=user --project=. ./runtests.jl

Then, if you have `lcov` installed, you should also have `genhtml`, and you can run this

    genhtml tracefile-<your_PID>.info --output-directory coverage/ && open coverage/index.html

to view the coverage locally as HTML.  I find that this sometimes requires
removing files that aren't really there from the .info file.

It's a well-hidden fact that you can turn coverage on and off by adding certain comments around the
code you don't want to measure:

    # COV_EXCL_START
    untested_code_that_wont_show_up_in_coverage()
    # COV_EXCL_STOP

"""

using Quaternionic
using Test, Random, Symbolics, StaticArrays, ForwardDiff

@variables w x y z a b c d e  # Symbolic variables

# NOTE: `FloatTypes` and `IntTypes` must be in descending order of width
FloatTypes = [BigFloat, Float64, Float32, Float16]
IntTypes = [BigInt, Int128, Int64, Int32, Int16, Int8]
SymbolicTypes = [Num]
Types = [FloatTypes...; IntTypes...; SymbolicTypes...]
PrimitiveTypes = [T for T in Types if isbitstype(T)]

QTypes = [Quaternion, Rotor, QuatVec]

# Handy assignments for now
Base.eps(::Quaternion{T}) where {T} = eps(T)
Base.eps(T::Type{<:Integer}) = zero(T)
Base.eps(n::Num) = zero(n)
Base.:≈(a::Num, b::Num; kwargs...) = Symbolics.simplify(a-b; expand=true) == 0

enabled_tests = lowercase.(ARGS)

help = ("help" ∈ enabled_tests || "--help" ∈ enabled_tests)
helptests = []

# This block is cribbed from StaticArrays.jl/test/runtests.jl
function addtests(fname)
    key = lowercase(splitext(fname)[1])
    if help
        push!(helptests, key)
    else
        if isempty(enabled_tests) || key in enabled_tests
            println("Running $key.jl")
            Random.seed!(42)
            include(fname)
        end
    end
end

@testset verbose=true "All tests" begin
    addtests("quaternion.jl")
    addtests("basis.jl")
    addtests("base.jl")
    addtests("algebra.jl")
    addtests("math.jl")
    addtests("random.jl")
    addtests("conversion.jl")
    addtests("distance.jl")
    addtests("alignment.jl")
    addtests("interpolation.jl")
    addtests("gradients.jl")
    addtests("doctests.jl")
end

if help
    println()
    println("Pass no args to run all tests, or select one or more of the following:")
    for helptest in helptests
        println("    ", helptest)
    end
end


text:
using PyPlot

# Plots incremental sampling estimates with a confidence region.
#
# Args:
# - n_samples: a Vector of independent sets of samples.
# - estimate_fn: a lambda that takes as input a Vector of samples and returns a
#     scalar metric.
# - delta_n: compute incremental estimates every this many steps.
# - label: a String label to use for the plot.
function plot_estimates(n_samples::Vector, estimate_fn::Function,
                        delta_n::Integer, label::String)
    # Compute the x range.
    if isa(first(n_samples), Tuple)
        total_n = length(first(n_samples)[1])
    else
        total_n = length(first(n_samples))
    end
    xl = [i for i in delta_n:delta_n:total_n]
    if (last(xl) != total_n); push!(xl, total_n); end

    # Compute incremental estimates.
    n_y = [[isa(samples, Tuple) ?
        estimate_fn(samples[1][1:x], samples[2][1:x]) :
        estimate_fn(samples[1:x]) for samples in n_samples] for x in xl]
    min_y, max_y = minimum.(n_y), maximum.(n_y)
    mid_y = (max_y + min_y) ./ 2.0

    # For Plots.jl:
    #  plot(x_mc, mid_y_mc, ribbon=(max_y_mc - mid_y_mc), fillalpha=0.15,
    #       label="MC", lw=2, xlabel="no. of samples", ylabel="estimates")
    #  plot(x_is, mid_y_is, ribbon=(max_y_is - mid_y_is), fillalpha=0.15, label="IS", lw=2)
    p = plot(xl, mid_y, label=label, lw=0.5)
    fill_between(xl, min_y, max_y, alpha=0.1)
    return p
end

# Plots a histogram of samples drawn.
function plot_histogram(samples::Vector, weights::Union{Vector,Nothing}=nothing;
                        bins::Integer=10, density::Bool=true,
                        label::String=nothing)
    PyPlot.plt.hist(samples, bins, weights=weights, density=density,
                    label=label, alpha=0.3)
end


text: nothing = 42

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
function ellipseBE(x0,y0,a,b,n)

	#==
	n points defining the ellipse ((x-x0)/a)^2+((y-y0)/b)^2=1
	counter-clockwise
	==#
	dtheta=2pi/n
	coordinates=zeros(n,2)
	for i in 1:n
		coordinates[i,1]=x0+a*cos((i-1)*dtheta)
		coordinates[i,2]=y0+b*sin((i-1)*dtheta)
	end

	# corresponding edges
	edges=[i+j for i in 1:n, j in 0:1]
	edges[n,2]=1

	# and length and inner normals
	length=zeros(n)
	normalsIn=zeros(n,2)
	for i in 1:n
		tx=coordinates[edges[i,2],1]-coordinates[edges[i,1],1]
		ty=coordinates[edges[i,2],2]-coordinates[edges[i,1],2]
		length[i]=norm([tx ty])
		normalsIn[i,:]=[-ty tx]/norm([tx ty])
	end

	return (coordinates=coordinates, edges=edges, length=length, normalsIn=normalsIn)

end


text:
abstract type NoDecision <: Exception end
Base.show(io::IO, nd::NoDecision) = print(io, """
    Planner failed to choose an action because the following exception was thrown:
    $nd

    To specify an action for this case, use the default_action solver parameter.
    """)

struct AllSamplesTerminal <: NoDecision
    belief
end
Base.show(io::IO, ast::AllSamplesTerminal) = print(io, """
    Planner failed to choose an action because all states sampled from the belief were terminal.

    To see the belief, catch this exception as ex and see ex.belief.

    To specify an action for this case, use the default_action solver parameter.
    """)

# Default action stuff was moved to MCTS
# This section can be deleted once a version of MCTS with default actions is tagged
if !isdefined(MCTS, :default_action)
    struct ExceptionRethrow end

    default_action(::ExceptionRethrow, pomdp, belief, ex) = rethrow(ex)
    default_action(f::Function, pomdp, belief, ex) = f(belief, ex)
    default_action(p::POMDPs.Policy, pomdp, belief, ex) = action(p, belief)
    default_action(s::POMDPs.Solver, pomdp, belief, ex) = action(solve(s, pomdp), belief)
    default_action(a, pomdp, belief, ex) = a

    """
        ReportWhenUsed(a)

    When the planner fails, returns action `a`, but also prints the exception.
    """
    struct ReportWhenUsed{T}
        a::T
    end

    function default_action(r::ReportWhenUsed, pomdp, belief, ex)
        showerror(STDERR, ex)
        a = default_action(r.a, pomdp, belief, ex)
        @warn("Using default action $a")
        return a
    end
end


text:
"""
PropDist(args...)--> UnivariateDistribution

Fit the given data to a specific statistical distribution
Distributions.jl package is used to fit

The next table show the list of variables allowed:

|PropertyName|Args|Default|Input|Description
|---|---|---|---|---|
|Prop|Mandatory|--|Array{Number,1}| Data to fit|
|Distr|Mandatory|--|Distr=Distribution|Distribution to fit the data given|
"""
@userplot prophist

@recipe function f(h::prophist)

    x, d=h.args

    MinProp=minimum(x)
    MaxProp=maximum(x)
    PropRange=range(MinProp,stop=MaxProp,length=100)
    DisProp=Distributions.fit(d, x[:]*1)
    DisParam=params(DisProp)
    FitProp=pdf.(DisProp, PropRange);
    Quant=quantile.(DisProp,[0.1,0.5,0.9])

@series begin
        seriestype := :histogram
        normalize := true
        bins := bins=round(Int,sqrt(size(x,1)))
        seriescolor --> :red
        label --> "Histogram"
        x
    end

@series begin
        seriestype := :path
        linestyle := :dash
        linewidth := 3
        linecolor --> :black
        label --> "Probabilitic Distribution Fitted"
        PropRange, FitProp

 end
println("Percentiles 10% - %50 - 90% \n         $Quant ")
end


text:
function generate_transition(domain::Domain, state::State,
                             domain_type::Symbol, state_type::Symbol)
    transition_def = quote
        function transition(domain::$domain_type, state::$state_type,
                            term::Term)
            execute(domain, state, get_action(domain, term.name), term.args)
        end
        function transition(domain::$domain_type, state::$state_type,
                            term::Term; check::Bool=false)
            execute(domain, state, get_action(domain, term.name), term.args;
                    check=check)
        end
        function transition!(domain::$domain_type, state::$state_type,
                             term::Term)
            execute!(domain, state, get_action(domain, term.name), term.args)
        end
        function transition!(domain::$domain_type, state::$state_type,
                             term::Term; check::Bool=false)
            execute!(domain, state, get_action(domain, term.name), term.args;
                     check=check)
        end
    end
    return transition_def
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

# -------------
# DOMAIN VIEWS
# -------------

"""
    DomainView(domain, indices)

A partial view of a `domain` containing only the elements at `indices`.
"""
struct DomainView{Dim,T,D<:Domain{Dim,T},I} <: Domain{Dim,T}
  domain::D
  inds::I
end

# convenience functions
Base.view(domain::Domain, inds) = DomainView(domain, inds)

# -----------------
# DOMAIN INTERFACE
# -----------------

Base.getindex(v::DomainView, ind::Int) =
  getindex(v.domain, v.inds[ind])

nelements(v::DomainView) = length(v.inds)

coordinates!(buff, v::DomainView, ind::Int) =
  coordinates!(buff, v.domain, v.inds[ind])

# -----------
# IO METHODS
# -----------

function Base.show(io::IO, v::DomainView)
  domain  = v.domain
  nelms = length(v.inds)
  print(io, "$nelms View{$domain}")
end

# -----------
# DATA VIEWS
# -----------

# helper function for table view
function viewtable(table, rows, cols)
  t = Tables.columns(table)
  v = map(cols) do c
    col = Tables.getcolumn(t, c)
    c => view(col, rows)
  end
  (; v...)
end

"""
    DataView(data, inds, vars)

Return a view of `data` at indices `inds` and variables `vars`.
"""
struct DataView{D<:Data,I,V} <: Data
  data::D
  inds::I
  vars::V
end

# convenience functions
Base.view(data::Data, inds::AbstractVector{Int}) =
  DataView(data, inds, collect(Tables.schema(values(data)).names))
Base.view(data::Data, vars::AbstractVector{Symbol}) =
  DataView(data, 1:nelements(domain(data)), vars)
Base.view(data::Data, inds, vars) = DataView(data, inds, vars)

# ---------------
# DATA INTERFACE
# ---------------

domain(v::DataView) = view(domain(v.data), v.inds)
values(v::DataView) = viewtable(values(v.data), v.inds, v.vars)

# specialize methods for performance
==(v₁::DataView, v₂::DataView) =
  v₁.data == v₂.data && v₁.inds == v₂.inds && v₁.vars == v₂.vars

# specialize view to avoid infinite loops
Base.view(v::DataView, inds::AbstractVector{Int}) =
  DataView(v.data, v.inds[inds], v.vars)
Base.view(v::DataView, vars::AbstractVector{Symbol}) =
  DataView(v.data, v.inds, vars)
Base.view(v::DataView, inds, vars) =
  DataView(v.data, v.inds[inds], vars)

# -----------
# IO METHODS
# -----------

function Base.show(io::IO, v::DataView)
  data  = v.data
  nelms = length(v.inds)
  print(io, "$nelms View{$data}")
end


text:
using BinDeps

@BinDeps.setup

hts_engine = library_dependency("libhts_engine_API",
                                aliases=["libhts_engine_API", "hts_engine_API-1"])

const libhts_engine_API_version = "1.0.9"

github_root = "https://github.com/r9y9/hts_engine_API"
arch = Sys.WORD_SIZE == 64 ? "x86_64" : "i686"
major = libhts_engine_API_version[1]

provides(Sources,
         URI("$(github_root)/archive/v$(libhts_engine_API_version).tar.gz"),
         hts_engine,
         unpacked_dir="hts_engine_API-$(libhts_engine_API_version)")

prefix = joinpath(BinDeps.depsdir(hts_engine), "usr")
srcdir = joinpath(BinDeps.depsdir(hts_engine), "src",
                  "hts_engine_API-$(libhts_engine_API_version)", "src")

provides(SimpleBuild,
          (@build_steps begin
              GetSources(hts_engine)
              @build_steps begin
                  ChangeDirectory(srcdir)
                  `./waf configure --prefix=$prefix`
                  `./waf build`
                  `./waf install`
              end
           end), hts_engine, os = :Unix)

@BinDeps.install Dict(:libhts_engine_API => :libhts_engine_API)


text:
using Primes

function ulamspiral(ord::Int)
    # Possible directions
    dirs = [[0, 1], [-1, 0], [0, -1], [1, 0]]
    # fdir = ["→", "↑", "←", "↓"] # for debug pourpose
    cur = maxsteps = 1  # starting direction & starting max steps
    steps = n = 0       # starting steps     & starting number in cell
    pos = [ord ÷ 2 + 1, isodd(ord) ? ord ÷ 2 + 1 : ord ÷ 2] # starting position
    M = Matrix{Bool}(ord, ord) # result matrix
    while n < ord ^ 2  # main loop (stop when the matrix is filled)
        n += 1
        M[pos[1], pos[2]] = isprime(n)
        steps += 1
        # Debug print
        # @printf("M[%i, %i] = %5s (%2i), step %i/%i, nxt %s\n", pos[1], pos[2], isprime(n), n, steps, maxsteps, fdir[cur])
        pos  .+= dirs[cur] # increment position
        if steps == maxsteps # if reached max number of steps in that direction...
            steps = 0        # ...reset steps
            if iseven(cur) maxsteps += 1 end # if the current direction is even increase the number of steps
            cur  += 1        # change direction
            if cur > 4 cur -= 4 end # correct overflow
        end
    end
    return M
end

mprint(m::Matrix) = for i in 1:size(m, 1) println(join(el ? " ∙ " : "   " for el in m[i, :]), '\n') end

M = ulamspiral(9)
mprint(M)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
import MittagLefflerFunctions.Quadrature.MLQuad
using SpecialFunctions: erfcx
using PyPlot

xmin, xmax = -0.4, 0.4
ymin, ymax = -2.5, 2.5
Nx = 201
Ny = 201
x = range(xmin, xmax, length=Nx)
y = range(ymin, ymax, length=Ny)

z = Complex{Float64}[ complex(x[j], y[i]) for i=1:Ny, j=1:Nx ]

E_par = MLQuad(0.5, 1.0, 10, :parabola, 0.2)
E_hyp = MLQuad(0.5, 1.0, 10, :hyperbola, 0.2)

err_par = Float64[ abs(E_par(z[i,j])-erfcx(-z[i,j])) for i=1:Ny, j=1:Nx ]
err_hyp = Float64[ abs(E_hyp(z[i,j])-erfcx(-z[i,j])) for i=1:Ny, j=1:Nx ]

figure(1)
#mesh(x, y, err_par)
contour(x, y, err_par, 12)
xlabel(L"$\Re z$")
ylabel(L"$\Im z$")
grid(true)
colorbar()
savefig("figure06L.pdf")

figure(2)
contour(x, y, err_hyp, 12)
xlabel(L"$\Re z$")
ylabel(L"$\Im z$")
grid(true)
colorbar()
savefig("figure06R.pdf")


text:
const hotkeys = Hotkey[
    Hotkey(
        "ctrl + alt + shift + s",
        SettingsWindow.showSettings
    ),
    Hotkey(
        "ctrl + shift + s",
        showFileSaveDialog
    ),
    Hotkey(
        "ctrl + s",
        menuFileSave
    ),
    Hotkey(
        "ctrl + o",
        showFileOpenDialog
    ),
    Hotkey(
        "ctrl + n",
        createNewMap
    ),
    Hotkey(
        "ctrl + m",
        MetadataWindow.configureMetadata
    ),
    Hotkey(
        "ctrl + shift + t",
        RoomWindow.configureRoom
    ),
    Hotkey(
        "ctrl + t",
        RoomWindow.createRoom
    ),
    Hotkey(
        "ctrl + shift + z",
        History.redo!
    ),
    Hotkey(
        "ctrl + z",
        History.undo!
    ),
    Hotkey(
        "ctrl + f",
        focusFilterEntry!
    ),
    Hotkey(
        "ctrl + plus",
        zoomIn!
    ),
    Hotkey(
        "ctrl + minus",
        zoomOut!
    )
]

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn


text:
export SpaceVehicle

import Base: popfirst!, pushfirst!

mutable struct SpaceVehicle
	activestage::Rocket
end

tail(r::StagedRocket) = r.nextstage
tail(r::Rocket) = nothing

function popfirst!(ship::SpaceVehicle)
    r = tail(ship.activestage)
    if r == nothing
        throw(throw(ArgumentError("no rocket stages left")))
    else
        discarded = ship.activestage
        discarded.nextstage = emptypayload
        ship.activestage = r
    end
    discarded
end

function pushfirst!(ship::SpaceVehicle, r::StagedRocket)
    r.nextstage = ship.activestage
    ship.activestage = r
    ship
end


text:
module MGEO

import Base: convert

using DelimitedFiles
using Distributed
using LinearAlgebra
using StaticArrays
using Random

################################################################################
#                                    Files
################################################################################

include("types.jl")
include("configuration.jl")
include("misc.jl")
include("mgeo_canonical.jl")
include("mgeo_var.jl")
include("mgeo_real.jl")

end # module


text:
# Copyright 2020 Google LLC

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

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

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

    sorted = sort(input)
    sorted = vcat([0], sorted, [sorted[end]+3])
    diffs = (@view sorted[2:end]) .- (@view sorted[1:end-1])

    one_jumps = count(x->x==1, diffs)
    three_jumps = count(x->x==3, diffs)
    length(diffs)
    return one_jumps * three_jumps
end

function number_of_ways(voltages::AbstractArray{Int})
    if length(voltages) <= 2
        return 1
    elseif voltages[1] + 3 == voltages[2]
        return number_of_ways(@view voltages[2:end])
    else
        group_size = 1
        for i ∈ 2:length(voltages)
            if voltages[i-1] + 1 == voltages[i]
                group_size += 1
            else
                break
            end
        end
        # We're at an end and there's only one adjacent number
        if group_size == 1
            return number_of_ways(@view voltages[2:end])
        # Each is attached to 3-jumps on the ends
        elseif group_size == 2
            return number_of_ways(@view voltages[3:end])
        # One bit (populate or don't populate the middle one)
        elseif group_size == 3
            return 2 * number_of_ways(@view voltages[4:end])
        # Two bits (middle 2 can be present or absent)
        elseif group_size == 4
            return 4 * number_of_ways(@view voltages[5:end])
        # Would be 3 bits, except that at least one of the middle ones must be populated
        elseif group_size == 5
            return 7 * number_of_ways(@view voltages[6:end])
        end
    end
end

function kinds_of_groupings(sorted)
    consecutive_groupings = Set()
    group_num = 0
    prev = nothing
    for n ∈ sorted
        if n-1 == prev
            group_num+=1
        elseif group_num >= 1
            push!(consecutive_groupings, group_num)
            group_num = 1
        end
        prev = n
    end
    return consecutive_groupings
end

function part2(filename = "day10input.txt")
    sorted = parse.(Int, eachline(filename)) |> sort
    sorted = vcat([0], sorted, [sorted[end]])
    number_of_ways(sorted)
end

println("Part 1: ", part1())
println("Part 2: ", part2())


text:
function Triplets(
    annotations::DataFrame,
    index::Symbol,
    ntriplets::Int,
    weights::Vector{Float64} = ones(Float64, size(annotations[:, Not(index)], 2))
    )

    nannotations, nannotators = size(annotations[:, Not(index)])
    samples = sampletriplets(nannotations, ntriplets)
    counter = 0

    S = TripletEmbeddings.tripletstype(nannotations)
    triplets = Vector{Triplet{S}}(undef, 0)
    D = [distances(column) for column in eachcol(annotations[:, Not(index)])]

    for t in samples
        less_than = 0
        greater_than = 0

        for l in 1:nannotators
            if !ismissing(D[l][t[1], t[2]]) && !ismissing(D[l][t[1], t[3]])
                if D[l][t[1], t[2]] < D[l][t[1], t[3]]
                    less_than += weights[l]
                elseif D[l][t[1], t[2]] > D[l][t[1], t[3]]
                    greater_than += weights[l]
                end
            end
        end

        if less_than > greater_than
            counter += 1
            push!(triplets, Triplet(t[1], t[2], t[3]))
        elseif less_than < greater_than
            counter += 1
            push!(triplets, Triplet(t[1], t[3], t[2]))
        end
    end

    return triplets
end

function distances(annotations::AbstractArray)
    n = length(annotations)
    D = zeros(Union{Float64,Missing}, n, n)

    for i in 1:n, j in i+1:n
        if !ismissing(annotations[i]) && !ismissing(annotations[j])
            D[i,j] = sqeuclidean(annotations[i], annotations[j])
        else
            D[i,j] = missing
        end
    end

    return D + D'
end

function distances(annotations::Vector{T}) where T <: Real
    return pairwise(SqEuclidean(), annotations', dims=2)
end

function Triplets(Y::AbstractMatrix{T}, n::Int) where T <: Real
    # triplets = [Triplets(Y[:, partition]) for partition in Iterators.partition(1:size(Y, 2), n)]
    for partition in Iterators.partition(1:size(Y, 2), n)
        println(partition)
        Triplets(Y[:, partition])
        # println(partition)
    end
end


text:
using Test
using QuantumOptics
using LinearAlgebra

@testset "nlevel" begin

N = 3
b = NLevelBasis(N)

# Test basis equality
@test b == NLevelBasis(N)
@test b != NLevelBasis(N+1)
@test_throws DimensionMismatch NLevelBasis(0)

# Test transition operator
@test_throws BoundsError transition(b, 0, 1)
@test_throws BoundsError transition(b, N+1, 1)
@test_throws BoundsError transition(b, 1, 0)
@test_throws BoundsError transition(b, 1, N+1)
@test dense(transition(b, 2, 1)) == basisstate(b, 2) ⊗ dagger(basisstate(b, 1))

# Test nlevel states
@test_throws BoundsError nlevelstate(b, 0)
@test_throws BoundsError nlevelstate(b, N+1)
@test norm(nlevelstate(b, 1)) == 1.
@test norm(dagger(nlevelstate(b, 1))*nlevelstate(b, 2)) == 0.
@test norm(dagger(nlevelstate(b, 1))*transition(b, 1, 2)*nlevelstate(b, 2)) == 1.

end # testset


text:
# # Quantum Circuit Born Machine

using Yao, YaoExtensions
import QuAlgorithmZoo
include("qcbmlib.jl")

# ## DATA: Target Probability Distribution
# The gaussian probability disctribution in phase space of 2^6
nbit = 6
N = 1<<nbit
pg = gaussian_pdf(1:N, N/2-0.5, N/4);

# ## MODEL: Quantum Circuit and Loss
# Using a random differentiable circuit of depth 6 for training, the kernel function is universal RBF kernel
depth = 6
kernel = rbf_kernel(0:N-1, 0.25)
c = variational_circuit(nbit, depth, pair_ring(nbit)) |> autodiff(:QC);
dispatch!(c, :random)
qcbm = QCBM(c, kernel, pg);

# ## TRAINING: Adam Optimizer
# We probide the QCBMGo! iterative interface for training
niter = 100
optim = Adam(lr=0.1)
for info in QCBMGo!(qcbm, optim, niter)
    curr_loss = loss(qcbm, info["probs"])
    println("Step = ", info["step"], ", Loss = ", curr_loss)
end


text:
"""
    cart2radial(x, y, args...)

Compute the radial coordinate from an `n`-dimensional Cartesian.

# Examples
```julia-repl
julia> Pda.cart2radial(3,4)
5.0
julia> Pda.cart2radial(1.0, 2.0, 2.0)
3.0
```
"""
@inline function cart2radial(x, y, args...)
	r2 = x^2 + y^2
	for arg in args
		r2 += arg^2
	end
	return sqrt(r2)
end

"""
	cart2radial(x::SVector)
"""
@inline function cart2radial(x::SVector{N,T}) where{N,T <: Real}
	return cart2radial(x...)
end

"""
    body2lab(X::SVector{2,T}, xc::SVector{2,T}, theta::Real)
Convert the coordinates of a point defined in the body frame `X` to
that in the lab frame.
xc: position of the origin of the body frame
theta: orientation of the body frame x axix relative to that of the lab

# Examples
```julia-repl
julia> body2lab(SVector{2,Float64}(1,0), SVector{2,Float64}(0,0), pi/4)
2-element SVector{2, Float64} with indices SOneTo(2):
 0.7071067811865476
 0.7071067811865475
```
"""
function body2lab(X::SVector{2,T},
    xc::SVector{2,T},
    theta::Real) where T
    x1 = xc[1] + X[1] * cos(theta) - X[2] * sin(theta)
    x2 = xc[2] + X[2] * cos(theta) + X[1] * sin(theta)
    return SVector{2,T}(x1, x2)
end


text:
include("../../../opti/src/opti_module.jl")

module noisydiff
using ..opti
using LinearAlgebra

function pt1_smooth(t, input_signal; T = 1)
    # trapezoidal rule

    L = length(input_signal)
    y = zeros(L)

    y[1] = input_signal[1]

    for kk = 2:L
        dt = t[kk] - t[kk-1]
        f = dt / (2 * T + dt) # 2 * T / dt + 1
        y[kk] = (input_signal[kk] + input_signal[kk-1] - y[kk-1]) * f + y[kk-1] * 2 * T / (dt + 2 * T)
    end

    return y
end

function diff_sliding_least_squares_d1(t, y; w = 5)
# fitting a degree 1 polynomial in the window and return the derivative

    N = length(y)
    u = zeros(N-w)

    yv = view(y, 1:w)
    xv = view(t, 1:w)
    U1 = sum(yv)
    U2 = sum(yv .* xv)
    S1 = sum(xv)
    S2 = sum(xv.^2)
    jj = 0

    for ii = (w+1):N
        jj    += 1

        d     = w * S2 - S1^2
        a1    = (U2 * w - S1 * U1) / d
        u[jj] = a1

        x_old   = t[ii-w]
        x_new   = t[ii]
        y_old   = y[ii-w]
        y_new   = y[ii]
        S1      += x_new - x_old
        S2      += x_new^2 - x_old^2
        U1      += y_new - y_old
        U2      += y_new * x_new - y_old * x_old

    end
    return u
end

function diff_sliding_least_squares_d3(t, y; w = 7)
    # fitting a degree 1 polynomial in the window and return the derivative

    N = length(y)
    u = zeros(N-w)

    yv = view(y, 1:w)
    xv = view(t, 1:w)
    U1 = sum(yv)
    U2 = sum(yv .* xv)
    U3 = sum(yv .* xv.^2)
    U4 = sum(yv .* xv.^3)

    S1 = sum(xv)
    S2 = sum(xv.^2)
    S3 = sum(xv.^3)
    S4 = sum(xv.^4)
    S5 = sum(xv.^5)
    S6 = sum(xv.^6)
    M = zeros(4, 4)
    b = zeros(4)
    jj = 0
    for ii = (w+1):N
        jj    += 1

        M[1,1] = w
        M[1,2] = S1
        M[1,3] = S2
        M[1,4] = S3

        M[2,2] = S2
        M[2,3] = S3
        M[2,4] = S4
        M[3,3] = S4
        M[3,4] = S5

        M[4,4] = S6
        b[1] = U1
        b[2] = U2
        b[3] = U3
        b[4] = U4

        S = Symmetric(M, :U)
        a = S \ b # todo : replace it with inplace operations

        x_old   = t[ii-w]
        x_new   = t[ii]
        y_old   = y[ii-w]
        y_new   = y[ii]

        u[jj] = a[2] + 2 * a[3] * x_new + 3 * a[4] * x_new^2

        S1      += x_new - x_old
        S2      += x_new^2 - x_old^2
        S3      += x_new^3 - x_old^3
        S4      += x_new^4 - x_old^4
        S5      += x_new^5 - x_old^5
        S6      += x_new^6 - x_old^6
        U1      += y_new - y_old
        U2      += y_new * x_new - y_old * x_old
        U3      += y_new * x_new^2 - y_old * x_old^2
        U4      += y_new * x_new^3 - y_old * x_old^3

    end
    return u
end

# function diff_tikhonov(t, y; w = 5, lamda = 0.1)
#     Ma = zeros(w, w)
#     Md = zeros(w, w)

#     N = length(t)
#     u = zeros(N)

#     for ii = (w+1):N

#         A = Ma + lamda * Md

#     end

#     return u
# end

function diff_fw(t, y)
    dt = diff(t)
    dy = diff(y)
    dy ./= dt
    return dy
end

function cumtrapezoidal(t, y)
    dt = diff(t)
    N = length(y)
    v1 = view(y, 1:N-1)
    v2 = view(y, 2:N)
    ys = v1 .+ v2
    ints = [0.0; 0.5 .* dt .* ys]
    ints = cumsum(ints)
    return ints
end

function tvd(t, y; lamda = 0.1, tol_options = opti.ToleranceOptions())

    # u0 = rand(length(y)) # copy(y)
    t_min = t[1]
    t_max = t[end]
    T = (t_max - t_min) / 100
    y_pred = pt1_smooth(t, y; T = T)
    u0 = [0.0; diff_fw(t, y_pred)]

    obj_fcn(u) = __obj_fcn(t, y, lamda, u)
    # stat = opti.cg_pr(obj_fcn, u0; tol_options = tol_options)
    (x_sol, y_sol, iter, stopping_crit) = opti.quasi_newton(obj_fcn, u0; tol_options = tol_options)
    return (x_sol, y_sol, iter, stopping_crit)
end

function __obj_fcn(t, y, lamda, u)
    i1 = __integral_1(u)
    i2 = __integral_2(t, y, u)
    return lamda * i1 + i2
end

function __integral_1(u)
    N = length(u)
    abs_du = 0.0
    @inbounds for kk = 1:(N-1)
        du = u[kk+1] - u[kk]
        abs_du += abs(du)
    end
    return abs_du / N
end

function __integral_2(t, y, u)
    N = length(y)

    # dt = diff(t)

    # vu1 = view(u, 1:N-1)
    # vu2 = view(u, 2:N)

    # int_u =[y[1]; 0.5 .* dt .* (vu1 .+ vu2)]
    # y_pred = cumsum(int_u)
    # g = (y_pred - y).^2
    # i1 = dt[1] * g[1]
    # iN = dt[end] * g[end]

    # ibulk = 0.0
    # for kk = 2:N-1
    #     ibulk += g[kk] * (dt[kk] + dt[kk-1])
    # end
    # ibulk += i1
    # ibulk += iN

    ibulk = 0.0
    y_pred = y[1]
    dt = t[2] - t[1]

    @inbounds for kk = 2:N-1
        dt1     = t[kk] - t[kk-1]
        us      = u[kk] + u[kk-1]
        int_u   = 0.5 * dt * us
        y_pred  += int_u
        g       = (y_pred - y[kk]).^2
        ibulk   += g * (dt + dt1)
        dt      = dt1
    end

    dt = t[N] - t[N-1]
    int_u = 0.5 * dt * (u[N] + u[N-1])
    y_pred += int_u
    g = (y_pred - y[N]).^2
    iN = dt * g
    # ibulk += i1
    ibulk += iN

    return ibulk
end

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
"""
    decompose_filnam(str)

Decompose filename into its name (and, if present, extension, prefix and numerator).
#### Examples:
```
strExample = "T23.01.fits"

dict = decompose_filnam(strExample)
Dict{String,String} with 4 entries:
  "Extension" => ".fits"
  "Numerator" => "01"
  "Prefix"    => "T23."
  "Name"      => "T23.01"

get(dict,"Numerator","Numerator: Key absent")
"01"

get(dict,"Wild","Key absent")
"Key absent"
```
"""
function decompose_filnam(str::String)

    ne = Base.findlast('.',str)                # ne: first digit of extension
    nl = length(str)                           # ne: length of file name including extension
    ne == nothing ? extension = false : extension = true

    if extension
        strNam = str[1:ne-1]
        strExt = str[ne:nl]
        Base.Unicode.isdigit(str[ne-1]) ? n = ne-1 : n = nothing
        o = [("Name", strNam), ("Extension", strExt)]
    else
        strNam = str[1:ne-1]
        strExt = nothing
        Base.Unicode.isdigit(str[nl]) ? n = nl : n = nothing
        o = [("Name", strNam)]
    end

    if n != nothing
        strNum = ""
        while Base.Unicode.isdigit(str[n])
            strNum = str[n] * strNum
            n -= 1
        end
        strPre = str[1:n]
        Base.append!(o,[("Prefix", strPre),("Numerator", strNum)])
    end

    return Dict(o)

end

"""
    fits_combine(str1, str2 [; info=false])

Combine a series of .fits files into a single .fits file.
#### Example:
```
fits_combine("T01.fits", "T22.fits"; info=false)
T01-T22.FITS: file was created (for more information set info=true)
```
"""
function fits_combine(filnamFirst::String, filnamLast::String; info=false)

    _file_exists(filnamFirst) || return "$filnamFirst: file not found"
    _file_exists(filnamLast) || return "$filnamLast: file not found"

    filnamFirst = uppercase(filnamFirst)
    filnamLast = uppercase(filnamLast)

    d = decompose_filnam(filnamFirst)
    strPre = get(d,"Prefix","Error: no prefix")
    strNum = get(d,"Numerator","Error: no Numerator")
    strExt = get(d,"Extension","Error: no extension")
    valNum = parse(Int,strNum )
    numLeadingZeros = length(strNum) - length(string(valNum))

    d = decompose_filnam(filnamLast)
    strPre2 = get(d,"Prefix","Error: no prefix")
    strNum2 = get(d,"Numerator","Error: no Numerator")
    strExt2 = get(d,"Extension","Error: no extension")
    valNum2 = parse(Int,strNum2 )
    numLeadingZeros2 = length(strNum2) - length(string(valNum2))

    if strPre ≠ strPre2
        error(strPre * " ≠ " * strPre2 * " (prefixes must be identical)")
    elseif strExt ≠ strExt2
        error(strExt * " ≠ " * strExt2 * " (file extensions must be identical)")
    elseif strExt ≠ ".FITS"
        error("file extension must be '.fits'")
    end

    numFiles = 1 + valNum2 - valNum
    fileFirst = FITSIO.FITS(filnamFirst)
    metaInfo = FITSIO.read_header(fileFirst[1])
    dataFirst = Base.read(fileFirst[1])  # read an image from disk
    Base.close(fileFirst)
    t = typeof(dataFirst[1,1,1])
    s = size(dataFirst)
    dataStack =  Array{t,3}(undef, s[1], s[2] , numFiles)

    itr = valNum:valNum2
    filnamNext = filnamFirst
    for i ∈ itr
        l = length(filnamNext)
        filnamNext = strPre * "0"^numLeadingZeros * string(i) * ".fits"
        if l < length(filnamNext)
            numLeadingZeros = numLeadingZeros -1
            filnamNext = strPre * "0"^numLeadingZeros * string(i) * ".fits"
        end
        fileNext = FITS(filnamNext)
        dataNext = Base.read(fileNext[1])  # read an image from disk
        Base.close(fileNext)
        dataStack[:, :,i] = dataNext[:, :,1]
    end

    filnamOut = strPre * strNum * "-" * strPre * strNum2 * strExt
    fileOut = FITSIO.FITS(filnamOut,"w")
    Base.write(fileOut, dataStack; header=metaInfo)
    if info
        println("Output fileOut:\r\n", fileOut)
        println("\r\nOutput fileOut[1]:\r\n", fileOut[1])
        Base.close(fileOut)
        println("\r\nmetaInformation:\r\n", metaInfo)
    else
        Base.close(fileOut)
        return filnamOut * ": file was created (for more information set info=true)"
    end
end


text:
using JuMP
using Gurobi
using CSV

function get_data()
  c_ij = CSV.read("supplier_costs.csv"; header = 0);

  return c_ij
end

function optim()

end

function main()
  c_ij = get_data()

end


text:
"""
Initializes values for a particle
* `model`: model containing a likelihood function with data and priors
* `de`: differential evolution object
* `p`: a particle
* `n_iter`: the number of iterations
"""
function init_particle!(model, de, p, n_iter)
    N = n_iter - de.burnin
    p.samples = typeof(p.samples)(undef,N,length(p.Θ))
    p.accept = fill(false, N)
    p.weight = priorlike(model, p) + model.model(p.Θ)
    p.lp = fill(0.0, N)
    return nothing
end

function priorlike(model, p::Particle{Array{T,1}}) where {T<:Real}
    return sum(logpdf.(model.priors, p.Θ))
end

function priorlike(model, p)
    LL = 0.0
    for (pr,θ) in zip(model.priors, p.Θ)
        LL += get_LL(pr[1], θ)
    end
    return LL
end

# handles arbitrary distribution types. Not flexible enough otherwise
get_LL(d::MultivariateDistribution, x::Vector) = logpdf(d, x)
get_LL(d::UnivariateDistribution, x::Real) = logpdf(d, x)
get_LL(d, x) = loglikelihood(d, x)

"""
Metropolis-Hastings proposal selection
Note: assumes weights are posterior log likelihoods
* `proposal`: weight of proposal e.g. posterior log likelihood
* `current`: weight of current value e.g. posterior log likelihood
"""
function accept(proposal, current)
    p = min(1.0, exp(proposal-current))
    rand() <= p ? (return true) : (return false)
end

"""
replaces values outside of bounds with the boundary
* `b`: boundary (lowerbound,upperbound)
* `θ`: a parameter value
"""
function enforce_bounds(b, θ::Real)
    θ < b[1] ? (return b[1]) : nothing
    θ > b[2] ? (return b[2]) : nothing
    return θ
end

function enforce_bounds(b, Θ)
    map(θ->enforce_bounds(b, θ), Θ)
end

function enforce_bounds!(bounds, p)
    i = 0
    for (b,θ) in zip(bounds, p.Θ)
        i += 1
        p.Θ[i] = enforce_bounds(b, θ)
    end
    return nothing
end

"""
Returns parameters names.
* `model`: model containing a likelihood function with data and priors
* `p`: a particle
"""
function get_names(model, p)
    N = length.(p.Θ)
    parm_names = fill("", sum(N))
    cnt = 0
    for (k,n) in zip(model.names, N)
        if n > 1
            for i in 1:n
                cnt += 1
                parm_names[cnt] = string(k, "[",i,"]")
            end
        else
            cnt += 1
            parm_names[cnt] = string(k)
        end
    end
    push!(parm_names, "acceptance", "lp")
    return parm_names
end

"""
Store samples after burnin period
Selects between mutation and crossover step with probability β
* `de`: differential evolution object
* `groups`: groups of particles
"""
function store_samples!(de, groups)
    de.iter <= de.burnin ? (return) : nothing
    i = de.iter - de.burnin
    for group in groups
        for p in group
            add_sample!(p, i)
        end
    end
    return nothing
end

function add_sample!(p::Particle{T}, i) where {T<:Real}
    p.samples[i,:] = p.Θ'
end

function add_sample!(p, i)
    p.samples[i,:] = p.Θ
end

function sample_prior(priors)
    p = [rand(p...) for p in priors]
    t = findtype(p)
    return returntype(t, p)
end

findtype(p) = Union{unique(typeof.(p))...}
returntype(t, p) = t[p...]

"""
Update particle based on Metropolis-Hastings rule.
* `de`: differential evolution object
* `groups`: groups of particles
"""
function update_particle!(de, current, proposal)
    @unpack iter,burnin = de
    i = iter - burnin
    accepted = accept(proposal.weight, current.weight)
    if accepted
        current.Θ = proposal.Θ
        current.weight = proposal.weight
    end
    if iter > burnin
         current.accept[i] = accepted
         current.lp[i] = current.weight
     end
    return nothing
end

# Type-stable arithmatic operations for Union{Array{Float64,1},Float64} types (which return Any otherwise)
import Base: +, - ,*

function +(x::Particle, y::Particle)
    N = length(x.Θ)
    z = similar(x.Θ)
    for i in 1:N
        z[i] = x.Θ[i] + y.Θ[i]
    end
    return Particle(Θ=z)
end

+(x::Real, y::Particle) = +(y, x)

function +(x::Particle, y::Real)
    N = length(x.Θ)
    z = similar(x.Θ)
    for i in 1:N
        z[i] = x.Θ[i] .+ y
    end
    return Particle(Θ=z)
end

function +(x::Particle, d::Distribution)
    N = length(x.Θ)
    z = similar(x.Θ)
    for i in 1:N
        z[i] = x.Θ[i] .+′ draw(d, x.Θ[i])
    end
    return Particle(Θ=z)
end

function draw(d, v::Float64)
    return rand(d)
end

function draw(d, v)
    return rand(d, length(v))
end

function *(x::Particle, y::Particle)
    N = length(x.Θ)
    z = similar(x.Θ)
    for i in 1:N
        z[i] = x.Θ[i] .*′ y.Θ[i]
    end
    return Particle(Θ=z)
end

*(x::Real, y::Particle) = *(y, x)

function *(x::Particle, y::Real)
    N = length(x.Θ)
    z = similar(x.Θ)
    for i in 1:N
        z[i] = x.Θ[i] .*′ y
    end
    return Particle(Θ=z)
end

*(x::Array{<:Real,1}, y::Particle) = *(y, x)

function *(x::Particle, y::Array{<:Real,1})
    N = length(x.Θ)
    z = similar(x.Θ)
    for i in 1:N
        z[i] = x.Θ[i] .* y[i]
    end
    return Particle(Θ=z)
end

function -(x::Particle, y::Particle)
    N = length(x.Θ)
    z = similar(x.Θ)
    for i in 1:N
        z[i] = x.Θ[i] - y.Θ[i]
    end
    return Particle(Θ=z)
end

-(x::Real, y::Particle) = -(y, x)

function -(x::Particle, y::Real)
    N = length(x.Θ)
    z = similar(x.Θ)
    for i in 1:N
        z[i] = x.Θ[i] .- y
    end
    return Particle(Θ=z)
end

# arithmatic methods for hanlding discrete parameters
*′(x, y) = x*y
*′(x::Int64, y::Float64) = Int(round(x*y))
*′(x::Float64, y::Int64) = Int(round(x*y))
*′(x::Array{Int64,N}, y::Float64) where {N} = @. Int(round(x*y))
*′(x::Float64, y::Array{Int64,N}) where {N} = @. Int(round(x*y))
+′(x, y) = x + y
+′(x::Int64, y::Float64) = Int(round(x+y))
+′(x::Float64, y::Int64) = Int(round(x+y))
+′(x::Array{Int64,N}, y::Float64) where {N} = @. Int(round(x+y))
+′(x::Float64, y::Array{Int64,N}) where {N} = @. Int(round(x+y))


text:
module CentroidalVoronoiController

using VoronoiDelaunay
using VoronoiCells
using GeometryBasics
using GraphUtils

export
    CVController,
    voronoi_controller,
    get_vel,
    get_vels,
    set_positions!,
    update_controller!

"""
    centroid(pts::Vector{V}) where {V<:AbstractVector}

Compute the centroid of the polygon whose vertices are defined (in ccw order) by
`pts`.
"""
function centroid(pts::Vector{V}) where {V<:AbstractVector}
    p1 = pts[1]
    total_area = 0.0
    c = zero(V)
    for (p2,p3) in zip(pts[2:end-1],pts[3:end])
        a = GeometryBasics.area([p1,p2,p3])
        total_area += a
        c = ((total_area-a)/total_area)*c .+ (a/total_area) * (p1 .+ p2 .+ p3) / 3
    end
    c
end

mutable struct CVController
    tess::VoronoiCells.Tessellation
    vels::CachedElement{Vector{Point2{Float64}}}
end

function set_vels!(controller::CVController,vels)
    update_element!(controller.vels,vels)
    # controller.vels = vels
end
function compute_vels!(controller::CVController)
    tess = controller.tess
    vels = [centroid(cell) - pt for (cell,pt) in zip(tess.Cells,tess.Generators)]
    set_vels!(controller,vels)
end

function get_vels(controller::CVController)
    if !is_up_to_date(controller.vels)
        compute_vels!(controller)
    end
    get_element(controller.vels)
end
get_vel(controller::CVController,idx::Int) = get_vels(controller)[idx]

function set_positions!(controller::CVController,pts)
    controller.tess = voronoicells(pts,controller.tess.EnclosingRectangle)
    set_up_to_date!(controller.vels,false)
    controller.tess.Generators
end

function update_controller!(controller::CVController,pts)
    set_positions!(controller,pts)
    compute_vels!(controller)
    controller
end

function voronoi_controller(pts::Vector{Point2{Float64}},
        rect::Rectangle=Rectangle(Point2(0.0,0.0),Point2(1.0,1.0))
        )
    controller = CVController(
        voronoicells(pts,rect),
        CachedElement(map(p->zero(Point2{Float64}),pts))
    )
    compute_vels!(controller)
    controller
end
voronoi_controller(pts::Vector{Point2{Float64}},rect::GeometryBasics.Rect) = voronoi_controller(pts,convert(Rectangle,rect))
Base.convert(::Type{Rectangle},rect::GeometryBasics.Rect) = Rectangle(Point2((rect.origin .- rect.widths)...), Point2((rect.origin .+ rect.widths)...))

end


text:
using DelimitedFiles
using LOWDER

include("generators.jl")

function runtest_mw(
                    filename::String;
                    unconstrained_prob::Bool=false
                    )

    source_filename = "CUTEr_selected_problems.dat"

    problems = readdlm( source_filename, Int64 )
    total_prob = size( problems )[ 1 ]

    file = open( filename, "w" )

    for i = 1 : total_prob

        print("Running: $( i ) of $( total_prob ) ... ")

        # Initializes useful constants
        nprob = problems[i, 1]
        n = problems[i, 2]
        p = problems[i, 3]
        rsp = convert( Bool, problems[i, 4] )
        n_points = n + 1

        try

            # Generates the problem
            ( x, l, u, fmin ) = problem_generator_mw( nprob, n, p, rsp; unconstrained = unconstrained_prob )

            # Solves the problem using 'lowder'
            sol = LOWDER.lowder( fmin, x, l, u; m = n_points )

            # Saves info
            text = "$( n ) $( p ) $( sol.status ) $( sol.true_val ) $( sol.iter ) $( sol.nf ) $( sol.nf / p ) $( sol.stationarity ) $( sol.sample_radius ) $( sol.tr_radius ) $( sol.index ) $( sol.f )"
            print( file, text )
            for j = 1 : ( n - 1 )

                print( file, " $(sol.solution[ j ])" )

            end
            println( file, " $( sol.solution[ n ] )" )

            println("succes!")

        catch

            println( file, "execution_fail" )

            println("fail!")

        end

    end

    close( file )

end

# -----------------------------------------------
# Path to file
# -----------------------------------------------

filename = "../data_files/mw_cons_test_beta_1.dat"

# -----------------------------------------------
# Funtion call
# -----------------------------------------------

runtest_mw( filename; unconstrained_prob = false )


text:
using DifferentialEquations,ParameterizedFunctions, DiffEqParamEstim #DiffEq
using RecursiveArrayTools, StatsBase,Distributions #Vector of Arrays and stats
using StatsPlots, CSV, DataFrames,Printf,Dierckx,ProgressMeter
using AverageShiftedHistograms, DelimitedFiles
###########################
#ToDo
#Figure out how to deal with parameter boundaries
#Include priors into MH MCMC?
###########################
cd()
cd(".\\Documents\\GitHub\\6StateMCMC\\")
#Import the MCMC solver
include("main.jl")
#Import error function
include("LossFunction.jl")
include("textalert.jl")

function Model!(dy,y,par,t)
  #IFN, ODE 1 parameters
  k11=0 #PR8, RIGI is assumed antagonized
  n=3
  RIGI=1
  k12=par[1]
  k13=par[2]
  k14=par[3]
  #IFN_env, ODE 2 parameters
  k21=par[4]
  tau2=par[5]
  #STATP, ODE 3 parameters
  k31=par[6]
  k32=par[7]
  k33=par[8]
  #IRF7, ODE 4 parameters
  k41=par[9]
  k42=par[10]
  #IRF7P, ODE 5 parameters
  k51=par[11]
  #Cells, ODE 6 parameters
  k61=par[12]
  #Virus, ODE 7 parameters
  k71=par[13]
  k72=par[14]
  k73=par[15]

  #ODE System
  dy[1]=y[6]*(k11*RIGI*y[7]+(k12*y[7]^n)/(k13+y[7]^n)+k14*y[5])-k21*y[1]
  dy[2]=k21*y[1]-tau2*y[2]
  dy[3]=(k31*y[2]*y[6])/(k32+k33*y[2])-0.3*y[3]
  dy[4]=y[6]*(k41*y[3]+k42*y[5])-0.3*y[4]
  dy[5]=k51*y[4]-0.3*y[5]
  dy[6]=-k61*y[6]*y[7]
  dy[7]=(k71*y[6]*y[7])/(1+k72*(y[2]*7E-5))-k73*y[7]
end

#parNames = ["a","b","c"]
#parNum = length(parNames)
stateNames = ["IFN","IFNe","STATP","IRF7","IRF7P","Live Cells","Virus"]
parNum = 15

#Define information for ODE model
p=rand(parNum) #Parameter values
u0 = [0, 0, 0, 0.72205, 0, 1, 6.9e-8] #Shifted Initial Conditions
tspan = (0,24.0) #Time (start, end)
shift= [7.939, 0, 12.2075, 13.4271, 0, 0, 0] #Noise shift on data
#Contruct the ODE Problem
prob = ODEProblem(Model!,u0,tspan,p)
alg = Vern7()   #ODE solver

#Solve for the "True" solution
sol = solve(prob,alg)

## Generate data

data = CSV.read("PR8.csv")
t = convert(Array,data[:Time])
mock = CSV.read("Control.csv")
titer = CSV.read("ViralTiters.csv")
t_titer=convert(Array,titer[:Time])

removeCol = [:Experiment :Time :IFN_env] #Strip time and IFN_env from data
for col in removeCol
  deletecols!(data,col)
  deletecols!(mock,col)
end

removeCol = [:Experiment :Time] #Strip time from viral titers
for col in removeCol
  deletecols!(titer,col)
end

data = convert(Matrix,data)
mock = convert(Matrix,mock)
titer = convert(Matrix,titer)

#Data points being measured
measured = [1,3,4]

## Supply prior distributions
priors = fill(Uniform(0,1), parNum)
parBounds = fill([0.0, Inf], parNum) #Fill bounds with 0/inf
#Modify bounds with known values
parBounds[5]=Float64[0.0, 2.5]
parBounds[12]=Float64[0.0, 1.0]
parBounds[15]=Float64[0.0, 0.5]

lossFunc = LossLog(t,t_titer,data,measured,mock,titer,shift)

sampleNum = Int(1e4)

result = ptMCMC(prob,alg,priors,parBounds,lossFunc,sampleNum)
bestPars= dropdims(permutedims(result[1], [1, 3, 2])[argmax(result[2],dims=1),:],dims=1)
pNew=bestPars[1,:]

CSV.write("chainparams.csv",DataFrame(result[1][:,:,1]))

open("acceptRatio.csv","a") do io #Write out acceptance ratios
  writedlm(io,result[3][:,1])
end
open("energy.csv","a") do io #Write out loss function values
  writedlm(io,result[2][:,1])
end

#Pull best results from the best chain
bestPars= dropdims(permutedims(result[1], [1, 3, 2])[argmax(result[2],dims=1),:],dims=1)
bestChains = [remake(prob;p=bestPars[i,:]) for i=1:size(bestPars,1)]
chainSols = [solve(problem,alg) for problem in bestChains]

###############
ODEplots = Vector(undef,length(u0))

mockSplines = [Spline1D(t,mock[:,i];k=1) for i=1:size(mock,2)]

plotRecipe = (t,x,i,j) -> plot(t,x,label="Chain $j",title=stateNames[i],framestyle=:box,yformatter=y->"$(@sprintf("%.0e",y))",legend=false)

Lfc(c,m) = @. log2(c) - log2(m)
LfcShift(c,m,s) = @. log2(c+s) - log2(m)
#loop through states
for i in eachindex(u0)

#Make the first plot for the first chain
  if i ∈ measured
      midx = findfirst(x->x==i,measured)
      #logScaledSim = Lfc(chainSols[1][i,:], mockSplines[midx](chainSols[1].t))
      logScaledSim = LfcShift(chainSols[1][i,:], mockSplines[midx](chainSols[1].t), shift[i])
      ODEplots[i] = plotRecipe(chainSols[1].t,logScaledSim,i,1)
      scatter!(t,data[:,findfirst(x->x==i,measured)],label="Data: "*stateNames[i])
  else
      ODEplots[i] = plotRecipe(chainSols[1].t,chainSols[1][i,:],i,1)
  end

#Loop through and add the rest of the chains
  for j=2:length(chainSols)
    if i ∈ measured
      midx = findfirst(x->x==i,measured)
      #logScaledSim = Lfc(chainSols[j][i,:], mockSplines[midx](chainSols[j].t))
      logScaledSim = LfcShift(chainSols[j][i,:], mockSplines[midx](chainSols[j].t), shift[i])
      plot!(chainSols[j].t,logScaledSim,label="Chain $j")
    else
      plot!(chainSols[j].t,chainSols[j][i,:],label="Chain $j")
    end
  end

end

datafitPlot = plot(ODEplots...,size=(1920,1080))

savefig(datafitPlot,"DataFit.pdf")

#Posterior
HistData = [plot(ash(result[1][1:end,i,1]),title="Par $i") for i=1:parNum]
histPlot=plot(HistData...,legend=false,size=(1920,1080))

savefig(histPlot,"histPlot.pdf")

logPostPlot = plot(result[2],legend =false,title="logPosterior")
savefig(logPostPlot,"logPostPlot.pdf")

#acceptance Rate
AcceptPlot = plot(result[3],legend =false,title="Acceptance Rate")
savefig(AcceptPlot,"AcceptPlot.pdf")

#number=("4")
#message=("Julia has finished job")
#TextAlert(number,message)

max_stars_repo_path: Julia Monte Carlo Fitting/Example.jl

max_issues_repo_path: Julia Monte Carlo Fitting/Example.jl

max_forks_repo_path: Julia Monte Carlo Fitting/Example.jl


text:
module TSPDrone

include("main.jl")

export solve_tspd, solve_tspd_RL,
        print_summary, generate_summary,
        TSPDroneResult

end


text:
@testset "Body" begin
    @testset "Geometry" begin
    a1 = 1; b1 = 0.1
    ccoeff = ComplexF64[0.5(a1+b1),0,0.5(a1-b1)]
    b = Bodies.ConformalBody(ccoeff)
    ζ = exp(im*0)
    v = exp(im*π/4)
    @test Bodies.normal(ζ,v,b) ≈ 0.7071067811865475
    @test Bodies.tangent(ζ,v,b) ≈ 0.7071067811865475
    @test conftransform(ζ,b) ≈ 1.0
    p = Bodies.Polygon([-1.0,1.0,1.0,-1.0],[-1.0,-1.0,1.0,1.0])
    b = Bodies.ConformalBody(p)
    ζ = exp(im*0)
    v = exp(im*π/4)
    @test Bodies.normal(ζ,v,b) ≈ 0.7071067811865475
    @test Bodies.tangent(ζ,v,b) ≈ 0.7071067811865475
    @test conftransform(ζ,b) ≈ 1.0
    end

    @testset "Motion" begin
    a1 = 1; b1 = 0.1
    ccoeff = ComplexF64[0.5(a1+b1),0,0.5(a1-b1)]
    b = Bodies.ConformalBody(ccoeff)
    Δt = 1.0
    ċ = 1.0im
    motion = RigidBodyMotion(ċ, 0.0);
    b₊ = deepcopy(b)
    advect!(b₊,b,motion,Δt);
    @test b₊.c ≈ 0.0 + 1.0im
    @test b₊.α ≈ 0.0
    ċ = 0.0im
    motion = RigidBodyMotion(ċ, 1.0);
    b₊ = deepcopy(b)
    advect!(b₊,b,motion,Δt);
    @test b₊.c ≈ 0.0 + 0.0im
    @test b₊.α ≈ 1.0
    end

end


text:
using PyCall
BayesOpt = pyimport("bayes_opt")
BayesOpt.BayesianOptimization
so = pyimport("scipy.optimize")

function black_box_function(x)
    # Function with unknown internals we wish to maximize.
    return x[1] ^ 2 + x[2] ^2
end

so.minimize(black_box_function,[1.,1.],method="Nelder-Mead")

bds = Dict('x' => (2, 4), 'y' => (-3, 3))
bds = Dict("x"=> (2,2))
bds = Dict(-2,2)
rates = []
Dict(bds)
Dict(i => bd for (i,bd) in enumerate(bds))

optimizer = BayesOpt.BayesianOptimization(
    f=black_box_function,
    pbounds=bds,
    verbose=2, # verbose = 1 prints only when a maximum is observed, verbose = 0 is silent
    random_state=1)

py"""
    from bayes_opt import BayesianOptimization

    def black_box_function(x, y):
        # Function with unknown internals we wish to maximize.
        return -x ** 2 - (y - 1) ** 2 + 1

    # Bounded region of parameter space
    pbounds = {'x': (2, 4), 'y': (-3, 3)}

    optimizer = BayesianOptimization(
        f=black_box_function,
        pbounds=pbounds,
        verbose=2, # verbose = 1 prints only when a maximum is observed, verbose = 0 is silent
        random_state=1,
    )

    optimizer.maximize(
        init_points=2,
        n_iter=3,
    )
"""
optimizer.maximize(
    init_points=2,
    n_iter=3)

py"""
import numpy as np

def sinpi(x):
    return np.sin(np.pi * x)
"""

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

skopt = pyimport("skopt")
noise_level = 0.1

function f(x,ct)
    x[1]^2 + x[2]^2 + ct
end

g = (x) -> f(x,ct)
ct = 1.0
# Plot f(x) + contours
x = range(-2, stop=2, length=400) #.reshape(-1, 1)

fx = [f(x_i) for x_i in x]
plot(x, fx, label="True (unknown)")
plt.fill(np.concatenate([x, x[::-1]]),
        np.concatenate(([fx_i - 1.9600 * noise_level for fx_i in fx],
                        [fx_i + 1.9600 * noise_level for fx_i in fx[::-1]])),
        alpha=.2, fc="r", ec="None")
plt.legend()
plt.grid()
plt.show()

res = skopt.gp_minimize(g,                  # the function to minimize
                  [(-2.0, 2.0),(-2.0, 2.0)],      # the bounds on each dimension of x
                  acq_func="EI",      # the acquisition function
                  n_calls=50,         # the number of evaluations of f
                  n_random_starts=5,  # the number of random initialization points
                  noise=0.,       # the noise level (optional)
                  random_state=1234)   # the random seed

res["x"]


text:
module Ladder

export Submission,
    LadderState,
    ladder,
    score!,
    empirical_loss

include("losses.jl")

type Submission
    id::String
    team::String
    labels::Array{Float64,1}
end

Submission(a::Array{Float64,1}) = Submission("","",a)

abstract LadderState

type FreeLadder <: LadderState
    labels::Array{Float64,1}
    current::Float64
    leader::Submission
end

type FixedLadder <: LadderState
    labels::Array{Float64,1}
    step_size::Float64
    current::Float64
    leader::Submission
end

function ladder(labels)
    return FreeLadder(labels,Inf,Submission(Float64[]))
end

function ladder(labels,step_size)
    return FixedLadder(labels,step_size,Inf,Submission(Float64[]))
end

# evalulate submission on parameter free ladder instance
function score!(l::FreeLadder,sub::Submission,lossfct)
    loss = empirical_loss(l.labels,sub.labels,lossfct)
    s = 0.0
    if length(l.leader.labels) == 0
        # no submission yet
        s = std(map(lossfct,l.labels,sub.labels))
    else
        # compute paired standard deviation between previous best and given labels
        s = std(map(lossfct,l.labels,sub.labels)-map(lossfct,l.labels,l.leader.labels))
    end
    # one sided paired t-test
    if loss < l.current - s/sqrt(length(l.labels))
        l.current = loss
        l.leader = sub
    end
    return l.current
end

# evaluate submision on ladder instance with fixed step size
function score!(l::FixedLadder,sub::Submission,lossfct)
    loss = empirical_loss(l.labels,sub.labels,lossfct)
    if loss < l.current - l.step_size
        l.current = round(loss,int(log(10,length(l.labels))))
        l.leader = sub
    end
    return l.current
end

end # module


text:
#
# Update status of pacients
#

function update_status( s1 :: Int64, s2 :: Int64, r :: Float64, input :: MDInput)
  # Check contatmination
  if r < input.xsec
    if  ( s1 == 1 || s2 == 1 ) && ( s1 < 3 && s2 < 3 ) # 3 nor 4 are dead and immune
      if rand() < input.pcont
        return 1, 1, true
      end
    end
    return s1, s2, true
  end
  return s1, s2, false
end


text:
using ApproxBayes
using Distributions
using DelimitedFiles
using Distances

# Run ABC assuming selection model, with three summary statistics

# export JULIA_NUM_THREADS=32

# read from file
dataset = ARGS[1]
ncell = ARGS[2]
epsilon = ARGS[3]
s3_weight = parse(Float64, ARGS[4])

# dataset="dataset_20190409_1"
# ncell=25

# dataset="dataset_20190409_2"
# ncell=13

println(dataset)
println(ncell)
println(epsilon)

min_brate = 0.2
max_brate = 2
model = 1

# change directory accordingly
dir = "../../"
progSim=dir * "/code/simpdo"
simdir=dir * "/data/abcsmc"
frealCnvec = dir * "/data/4trees/" * dataset * "_rvec.txt"
frealTree = dir * "/data/4trees/" * dataset * "_divtree.txt"
freal_sstat = dir * "/data/4trees/" * dataset * "_sstat3.txt"

targetdataCnvec = readdlm(frealCnvec)[:, 1]
# skip header and the first three cell divisions
targetTree = readdlm(frealTree, skipstart=7)[:, 3]
tblen = sum(targetTree)*3/1440/2
push!(targetdataCnvec, tblen)
targetdata = targetdataCnvec
println(targetdata)

writedlm(freal_sstat, transpose(targetdata))

# Get the distance between simulated and real data (read from files) from copy numbers
function tumourABCneutralcn(params, constants, targetdata)
  # call simulation program to generate data
  # print(params)
  mu1=params[1]
  mu2=params[2]
  birth_rate = params[3]
  fitval = "1 " * string(params[4])
  death_rate = 0

  ncell=constants[1]

  seed=rand(UInt32, 1)[1]
  suffix=string(seed) * "_" * dataset * "_n" * string(ncell)

  skip = 3

  cmdSim=`$progSim -o $dir -e $ncell --model $model --fitness "$fitval" --chr_prob $mu1 --arm_prob $mu2 --birth_rate $birth_rate --death_rate $death_rate --skip $skip --seed $seed --suffix $suffix`
  res = chomp(read(cmdSim, String))
  arr = split(res,"\t")
  simdata = map(x->parse(Float64,x), arr)
  simdata = simdata[1:3]

  simdata = vcat((simdata[1]), (simdata[2]), s3_weight * log(simdata[3]))

  r = euclidean(targetdata, simdata)

  return r, 1
end

setup = ABCSMC(tumourABCneutralcn, #simulation function
    4, # number of parameters
    parse(Float64, epsilon), # target ϵ
    Prior([Uniform(0, 1), Uniform(0, 1), Uniform(min_brate, max_brate), Uniform(-1, 1)]), #Prior for each of the parameters
    maxiterations = 9*10^6, #Maximum number of iterations before the algorithm terminates
    nparticles = 500,
    α = 0.3,
    ϵ1 = 10000.0,
    convergence = 0.05,
    constants = [ncell],
  )

targetdata = vcat((targetdata[1]), (targetdata[2]), s3_weight * log(targetdata[3]))

smc = runabc(setup, targetdata, verbose = true, progress = true, parallel = true)

fout="smc-joint-withsel-" * dataset * "_n" * string(ncell) * "_epsilon" * string(epsilon) * "_logs3_weight" * string(s3_weight)
writeoutput(smc, dir=simdir, file=fout)


text:
using ZXCalculus
using ZXCalculus: ancilla_extraction

function gen_phase_gadget()
    zxd = ZXDiagram(2)
    push_gate!(zxd, Val(:Z), 1, 1//2)
    push_gate!(zxd, Val(:CNOT), 1, 2)
    push_gate!(zxd, Val(:Z), 1, 1//4)
    push_gate!(zxd, Val(:CNOT), 1, 2)
    push_gate!(zxd, Val(:Z), 1, 1//2)

    return zxd
end

zxd = gen_phase_gadget()
zxg = full_reduction(zxd)
anc_circ = ancilla_extraction(zxg)

# plot(anc_circ)

zxd_swap = ZXDiagram(2)
pushfirst_gate!(zxd_swap, Val(:SWAP), [1, 2])
# plot(zxd_swap)
convert_to_chain(zxd_swap)

zxg_swap = ZXGraph(zxd_swap)
zxd_anc = ancilla_extraction(zxg_swap)
# plot(zxd_anc)
@test length(ZXCalculus.convert_to_chain(zxd_anc)) == 3


text:
using DelaySSAToolkit
using Test, SafeTestsets

@time begin
    @time @safetestset "Algorithm accuracy" begin
        include("bursty_model.jl")
    end
    @time @safetestset "dep_gr_delay test" begin
        include("dep_gr_delay.jl")
    end
    @time @safetestset "save delay channel test" begin
        include("save_delay_channel.jl")
    end
    @time @safetestset "cascade of delay reaction test" begin
        include("cascade_of_delay_reaction.jl")
    end
    # @time @safetestset "Index test" begin
    #     include("check_index_reactions.jl")
    # end
end


text:
# -------------- Core functions --------------

# Packages required
using LinearAlgebra, Parameters, DifferentialEquations, RecursiveArrayTools
using Statistics: mean, std
using PyPlot

# Reference set of parameters
@with_kw mutable struct EcoPar
    In1 = 0.3
    In2 = 0.25  # used in slow node
    rn  = 0.05
    f   = 0.8
    e1  = 0.05
    e2  = 0.05
    e3  = 0.05
    a1  = 0.6
    a11 = 0.4
    a23 = 0.07
    a2  = 0.25
    h42 = 0.005   # => treated as h32 in NCDRe
    h43 = 0.005
    b1  = 0.040
    b11 = 0.047
    b2  = 0.040
    a42 = 0.250   # => treated as a32 in NCDRe
    a43 = 0.09
    dD  = 0.1
    d1  = 0.10
    d11 = 0.11
    df2 = 0.001   # used in slow node
    df  = 0.05
    d2  = 0.001
    hub_slow = false
end

# I- ODE - Edible resource

## interaction + diffusion
function eco_NDCRe!(du, u, p, t)
    @unpack In1, df = p
    #
    j = 1
    for i in 1:5
        ind = j:(j + 3)
        du[ind] = eco_NDCRe_unit!(du[ind], u[ind], p, t)
        j += 4
    end
    # diffusion and input
    du[1] += In1 - df * u[1]
    du[5] += In1 - df * u[5]
    du[9] += In1 - df * u[9]
    du[13] += In1 + df * (u[1] + u[5] + u[9])/3. - df * u[13]
    du[17] += In1 + df * u[13]

    return du
end

## interaction + diffusion - slow hub node (sh)
function eco_NDCRe_sh!(du, u, p, t)
    @unpack In1, In2, df, df2 = p
    #
    j = 1
    for i in 1:5
        ind = j:(j + 3)
        du[ind] = eco_NDCRe_unit!(du[ind], u[ind], p, t)
        j += 4
    end
    # diffusion and input
    du[1] += In1 - df * u[1]
    du[5] += In1 - df * u[5]
    du[9] += In1 - df * u[9]
    du[13] += In2 + df * (u[1] + u[5] + u[9])/3. - df2 * u[13]
    du[17] += In1 + df2 * u[13]

    return du
end

## core ODE
function eco_NDCRe_unit!(du, u, p, t)
    # u[1]: Phosphorus; u[2]: Resource; u[3]: Consumer; u[4]: Detritus
    @unpack rn, a1, a23, a2, a42, h42, b1, f, e1, e2, e3, dD, d1, d2 = p

    tmp1 = a1 * u[1] * u[2] / (b1 + u[1])
    tmp24 = a42 * u[2] * u[3] / (1. + a42 * h42 * u[2])
    du[1] = - rn * u[1] - tmp1 + dD * u[4]
    du[2] = f * tmp1 - tmp24 - (d1 + e1) * u[2]
    du[3] = f * tmp24 - (d2 + e2) * u[3]
    du[4] = d1 * u[2] + d2 * u[3] - (dD + e3) * u[4] # + d11 * u[3]

    return du
end

# II- ODE - edible and less edible resource

## interaction + diffusion
function eco_NDCReRle!(du, u, p, t)
    @unpack In1, df= p;
    #
    j = 1
    for i in 1:5
        ind = j:(j + 4)
        du[ind] = eco_NDCReRle_unit!(du[ind], u[ind], p, t)
        j += 5
    end
    # diffusion and input
    du[1] += In1 - df * u[1]
    du[6] += In1 - df * u[6]
    du[11] += In1 - df * u[11]
    du[16] += In1 + df * (u[1] + u[6] + u[11])/3. - df * u[16]
    du[21] += In1 + df * u[16]

    return du
end

## interaction + diffusion - slow hub node (sh)
function eco_NDCReRle_sh!(du, u, p, t)
    @unpack In1, In2, df, df2 = p;
    #
    j = 1
    for i in 1:5
        ind = j:(j + 4)
        du[ind] = eco_NDCReRle_unit!(du[ind], u[ind], p, t)
        j += 5
    end
    # diffusion and input
    du[1] += In1 - df * u[1]
    du[6] += In1 - df * u[6]
    du[11] += In1 - df * u[11]
    du[16] += In2 + df * (u[1] + u[6] + u[11])/3. - df2 * u[16]
    du[21] += In1 + df2 * u[16]

    return du
end

## core ODE
function eco_NDCReRle_unit!(du, u, p, t)
    # u[1]: Phosphorus; u[2]: Resource; u[3]: Resource less edible;
    # u[4]: Consumer; u[5]: Detritus
    @unpack rn, a1, a11, a23, a2, a42, a43, h42, h43, b1, b11, f, e1, e2, e3, dD, d1, d11, d2 = p

    tmp1 = a1 * u[1] * u[2] / (b1 + u[1])
    tmp11 = a11 * u[1] * u[3] / (b11 + u[1])
    K = 1. + a42 * h42 * u[2] + a43 * h43 * u[3]
    tmp24 = a42 * u[2] * u[4] / K
    tmp34 = a43 * u[3] * u[4] / K
    du[1] = -rn * u[1] - tmp1 - tmp11 + dD * u[5]
    du[2] = f * tmp1 - tmp24 - (d1 + e1) * u[2]
    du[3] = f * tmp11 - tmp34 - (d11 + e1) * u[3]
    du[4] = f * tmp24 + f * tmp34 - (d2 + e2) * u[4]
    du[5] = d1 * u[2] + d11 * u[3] + d2 * u[4] - (dD + e3) * u[5]

    return du
end

# III- Bifurcation analysis

# we use bif_analysis2 cause bif_analysis was used in a previous version
# see SI where it is still used.

function bif_analysis2(vals, FUN, u0, tspan, ind, pbif = EcoPar())
    nsp = size(u0, 1)
    mat_max = zeros(Float64, nsp, size(vals, 1))
    mat_min = zeros(Float64, nsp, size(vals, 1))
    mat_mean = zeros(Float64, nsp, size(vals, 1))
    mat_cv  = zeros(Float64, nsp, size(vals, 1))

    for (i, df) in enumerate(vals)
        printstyled("rep #", i, "  df=", df, "\n", color = :blue)
        pbif.df = df
        prob = ODEProblem(FUN, u0, tspan, pbif)
        # sol = solve(prob, Trapezoid(), reltol = 1e-16)
        sol = solve(prob, Rodas4P(), reltol = 1e-16)
        mat_max[:, i] = minimum(sol(ind), dims = 2)#[local_max(sol(ind)[k, :])[1] for k in 1:nsp]
        mat_min[:, i] = maximum(sol(ind), dims = 2)#[local_min(sol(ind)[k, :])[1] for k in 1:nsp]
        tmp = mean(sol(ind), dims = 2)
        tmp2 = std(sol(ind), dims = 2)
        mat_mean[:, i] = mean(sol(ind), dims = 2)
        mat_cv[:, i] = tmp2 ./ tmp
    end

    return (vals, mat_min, mat_max, mat_mean, mat_cv)
end

# IV- Figures

function fig_bif(bif, ind, file, xlab, ylab, ymin, ymax, col = "black", lw = 1.8)
    figure(figsize = (2.8, 9))

    subplot(311)
    plot(bif[1], bif[2][ind[1], :], color = col, linestyle = "-", linewidth = lw)
    plot(bif[1], bif[3][ind[1], :], color = col, linestyle = "-", linewidth = lw)
    xlabel(xlab, fontname = "Times New Roman", fontsize = 12)
    # ylabel(ylab)
    xlim(0., 1.0)
    ylim(ymin, ymax)

    subplot(312)
    plot(bif[1], bif[2][ind[2], :], color = col, linestyle = "-", linewidth = lw)
    plot(bif[1], bif[3][ind[2], :], color = col, linestyle = "-", linewidth = lw)
    xlabel(xlab, fontname = "Times New Roman", fontsize = 12)
    # ylabel(ylab)
    xlim(0., 1.0)
    ylim(ymin, ymax)

    subplot(313)
    plot(bif[1], bif[2][ind[3], :], color = col, linestyle = "-", linewidth = lw)
    plot(bif[1], bif[3][ind[3], :], color = col, linestyle = "-", linewidth = lw)
    xlabel(xlab, fontname = "Times New Roman", fontsize = 12)
    # ylabel(ylab)
    xlim(0., 1.0)
    ylim(ymin, ymax)

    tight_layout()
    savefig(file, dpi = 300)
    plt.close()

end

function barplotmean(bif, ind, file, ymin, ymax, col = "black")
    figure(figsize = (2.6, 3.8))
    bar(1:3, bif[4][ind, 81][1,:], color = col, align = "center")
    ylim(ymin, ymax)
    tight_layout()
    savefig(file, dpi = 300)
    plt.close()
end

function fig_bif2(bif, ind, ind2, col, file, xlab, ylab, ymin, ymax)
    figure(figsize = (2.8, 9))

    subplot(311)
    plot(bif[1], bif[4][ind[1], :], color = col[1], linestyle = "-")
    plot(bif[1], bif[4][ind2[1], :], color = col[2], linestyle = "-")
    xlim(0., 1.0)
    ylim(ymin, ymax)

    subplot(312)
    plot(bif[1], bif[4][ind[2], :], color = col[1], linestyle = "-")
    plot(bif[1], bif[4][ind2[2], :], color = col[2], linestyle = "-")
    xlim(0., 1.0)
    ylim(ymin, ymax)

    subplot(313)
    plot(bif[1], bif[4][ind[3], :], color = col[1], linestyle = "-")
    plot(bif[1], bif[4][ind2[3], :], color = col[2], linestyle = "-")
    xlabel(xlab, fontname = "Times New Roman", fontsize = 12)
    xlim(0., 1.0)
    ylim(ymin, ymax)

    tight_layout()
    savefig(file, dpi = 300)
    plt.close()

end

function fig_bif3(bif, ind, col, file, xlab, ylab, ymin, ymax)
    figure(figsize = (2.8, 9))

    subplot(311)
    for i in 1:3
        plot(bif[1], bif[2][ind[i][1], :], color = col[i], linestyle = "-")
        plot(bif[1], bif[3][ind[i][1], :], color = col[i], linestyle = "-")
    end
    xlim(0., 1.0)
    ylim(ymin, ymax)

    subplot(312)
    for i in 1:3
        plot(bif[1], bif[2][ind[i][2], :], color = col[i], linestyle = "-")
        plot(bif[1], bif[3][ind[i][2], :], color = col[i], linestyle = "-")
    end
    xlim(0., 1.0)
    ylim(ymin, ymax)

    subplot(313)
    for i in 1:3
        plot(bif[1], bif[2][ind[i][3], :], color = col[i], linestyle = "-")
        plot(bif[1], bif[3][ind[i][3], :], color = col[i], linestyle = "-")
    end
    xlim(0., 1.0)
    ylim(ymin, ymax)

    tight_layout()
    savefig(file, dpi = 300)
    plt.close()

end

# V- Print results

function show_results(sol, ind, nsp = 5, nnd = 5)

    nm = ["Nutrients" "Resources" "Resources inedible" "Consumers" "Detritus"]

    for i in 1:nsp
        printstyled(nm[i], "\n", color = :blue)
        [@show mean(sol(ind)[k, :]) for k in i:nsp:((nnd - 1)*nsp + i)]
    end

    return mean(sol(ind), dims = 2)
end

max_stars_repo_name: McCannLab/NDID

max_issues_repo_name: McCannLab/NDID

max_forks_repo_name: McCannLab/NDID


text:
import LinearAlgebra: svd,qr

function TNrSVD(A::MPT{4},R::Int,s::Int,q::Int,acc::Float64)
# Algorihtm 5.5 from Batselier 2018
#######################################

    # generate random Matrix (Lemma 5.1)
    O = genrandTTm(size(A),R+s)
    # compute orthogonal basis (Alg. 5.4)
    Q = randTTmsubspaceiter(q,A,O,acc)
    # Compute B (subsection 5.3)
    B = transpose(Q)*A
    # Compute econ. SVD (Alg. 5.3)
    W,S,V = svd(B)
    # Compute U
    U = Q
    U[1] = modenprod(U[1],W')

    return U,S,V
end

function genrandTTm(rowsizes::Vector{Tuple{Int64, Int64}},colsize::Int)
    # generate random Matrix (Lemma 5.1)
    d     = length(rowsizes)
    cores = Vector{Array{Float64,4}}(undef,d)
    cores[1] = randn(1,rowsizes[1][2],colsize,1)
    for i = 2:d
        cores[i] = randn(1,rowsizes[i][2],1,1)
    end
    return MPT(cores)
end

function randTTmsubspaceiter(q::Int,A::MPT{4},O::MPT{4},acc::Float64)
    # compute orthogonal basis (Alg. 5.4 + 5.2)
    Y = A*O
    Q = qr(Y)
    for i = 1:q
        Y = transpose(A)*Q
        Q = qr(Y)
        Q = roundTT(Q,acc)
        Y = A*Q
        Q = qr(Y)
        Q = roundTT(Q,acc)
    end
end

function qr(A::MPT{4},ϵ::Float64)
    d = order(A)
    r = rank(A)
    s = size(A)
    frobnorm    = norm(A);
    err2        = 0;
    Av = mpo2mps(A)

    δ = ϵ / sqrt(d-1) * frobnorm;
    for k = 1:d-1
        F     = svd(unfold(Av[k],[3],"right"),full=false);
        rcurr = length(F.S);
        sv2   = cumsum(reverse(F.S).^2);
        tr    = Int(findfirst(sv2 .> δ^2))-1;
        if tr > 0
            rcurr = length(F.S) - tr;
            err2 += sv2[tr];
        end
        Utr = F.U[:,1:rcurr];
        Str = F.S[1:rcurr];
        Vtr = F.Vt[1:rcurr,:];

        Av[k]   = reshape( Utr,( Int(length(Utr)/(size(Av[k],2)*size(Vtr,1))), size(Av[k],2), size(Vtr,1) ) );
        Av[k+1] = nmodeproduct(Matrix((Vtr'*Diagonal(Str))'),Av[k+1],1);
    end

    s1 = [s[i][1] for i = 1:d]
    s2 = [s[i][2] for i = 1:d]
    sz = [s1';s2']
    A = mps2mpo(Av,sz)

    cores = Vector{Array{Float64,4}}(undef,d)
    for i = d:-1:2
        Ai        = reshape(A[i],(r[i][1],s[i][1]*r[i][2]))
        F         = qr(Ai')
        cores[i]  = reshape(Array(F.Q'),(r[i][1],s[i][1],1,r[i][2]))
        A[i-1]    = nmodeproduct(Array(F.R),A[i-1],4)
    end
    A1t = permutedims(A[1],[3 1 2 4])
    A1  = reshape(A1t,(s[1][2], s[1][1]*r[1][2]))
    F   = qr(A1')
    Q1t = reshape(Array(F.Q)',(1, s[1][1], r[1][2], s[1][2]))
    cores[1] = permutedims(Q1t,[1,2,4,3])

    return MPT(cores),F.R
end

function svd(A::MPT{4})
    # Algorithm 5.3
    d = order(A)
    r = rank(A)
    s = size(A)
    cores = Vector{Array{Float64,4}}(undef,d)
    for i = d:-1:2
        Ai       = reshape(A[i],(r[i][1],s[i][1]*r[i][2]))
        F        = qr(Ai)
        cores[i] = reshape(Array(F.Q),(r[i][1],s[i][1],1,r[i][2]))
        A[i-1]   = nmodeproduct(Array(F.R),A[i-1],4)
    end
    A1t      = permutedims(A[1],[3 1 2 4])
    A1       = reshape(A1t,(s[1][2], s[1][1]*r[1][2]))
    W,S,Q1   = svd(A1)
    Q1t      = reshape(Q1,(s[1][2], 1, s[1][1], r[1][2]))
    cores[1] = permutedims(Q1t,[2 3 1 4])

    return W,S,MPT(cores)
end

max_stars_repo_name: clarazen/BigMat

max_issues_repo_name: clarazen/BigMat

max_forks_repo_name: clarazen/BigMat