text:
 #=

Implementations of cg and pcg that use BLAS when they can.

Started by Dan Spielman.
Contributors:

=#

"""
`cg(mat, b; tol, maxits, maxtime, verbose)` solves a symmetric linear system.
`tol` is set to 1e-6 by default,
`maxits` defaults to Inf
`maxtime` defaults to Inf.  It measures seconds.
`verbose` defaults to false
"""
function cg end

"""
`pcg(mat, b, pre; tol, maxits, maxtime, verbose)` solves a symmetric linear system using preconditioner `pre`.
`pre` can be a function or a matrix.  If a matrix, a function to solve it is created with cholFact.
`tol` is set to 1e-6 by default,
`maxits` defaults to Inf
`maxtime` defaults to Inf.  It measures seconds.
`verbose` defaults to false
"""
function pcg end

"""
`cgSolver(mat; tol, maxits, maxtime, verbose)` creates a solver for a PSD system `mat`.
The parameters are as described in cg.
"""
function cgSolver end

"""
`pcgSolver(mat, pre; tol, maxits, maxtime, verbose)` creates a solver for a PSD system using preconditioner `pre`.
The parameters are as described in pcg.
"""
function pcgSolver end

function cg(mat, b::Array{Float64,1}; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false)
    cgBLAS(mat, b, tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose)
end

function cg(mat, b::Array{Float32,1}; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false)
    cgBLAS(mat, b, tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose)
end

function cg(mat, b; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=true)
    cgSlow(mat, b, tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose)
end

function cgSolver(mat; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false)
    f(b; tol::Real=1e-6, maxits=maxits, maxtime=maxtime, verbose=verbose) = cg(mat, b, tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose)
end

function pcg(mat, b, pre::Union{AbstractArray,Matrix}; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false)
    fact = cholfact(pre)
    F = x->(fact \ x)
    pcg(mat, b, F; tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose)
end

function pcg(mat, b::Array{BigFloat,1}, pre::Function; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false)
    pcgBLAS(mat, b, pre, tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose)
end

function pcg(mat, b::Array{Float64,1}, pre::Function; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false)
    pcgBLAS(mat, b, pre, tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose)
end

function pcg(mat, b::Array{Float32,1}, pre::Function; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false)
    pcgBLAS(mat, b, pre, tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose)
end

function pcg(mat, b, pre::Function; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false)
    pcgSlow(mat, b, pre, tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose)
end

function pcgSolver(mat, pre; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false)
    f(b; tol::Real=1e-6, maxits=maxits, maxtime=maxtime, verbose=verbose) = pcg(mat, b, pre, tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose)
end

"""Create a solver that uses cg to solve Laplacian systems in mat. Specialized for the case when pre is a Laplacian matrix.  Fix the default parameters of the solver as given"""
function pcgLapSolver(mat, pre; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false)
    fact = lapChol(pre)
    f(b; tol::Real=1e-6, maxits=maxits, maxtime=maxtime, verbose=verbose) = pcg(mat, b, fact, tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose)

end

# uses BLAS.  As fast as Matlab's pcg.
# set to use similar paramaters
function cgBLAS{Tval}(mat, b::Array{Tval,1};
        tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false)

    n = size(mat,2)

    x = zeros(Tval,n)

    tol = tol * norm(b)
    r = copy(b)
    p = copy(r)

    # the following would have higher accuracy
    #       rho = sum(r.^2)
    rho = norm(r)^2

    t1 = time()

    itcnt = 0
    while itcnt < maxits
        itcnt = itcnt+1

        q = mat*p

        al = rho/dot(p, q)

        BLAS.axpy!(al,p,x)  # x = x + al * p

        BLAS.axpy!(-al,q,r)  # r -= al*q

        if norm(r) < tol #Converged?
            break
        end

        oldrho = rho
        rho = norm(r)^2

        # the following would have higher accuracy
        #       rho = sum(r.^2)

        beta = rho/oldrho

        # p = r + beta*p
        BLAS.scal!(n,beta,p,1) # p *= beta
        BLAS.axpy!(1.0,r,p) # p += r

        if (time() - t1) > maxtime
            if verbose
                println("CG stopped at maxtime.")
            end
            break
        end
      end

    if verbose
        println("CG stopped after: ", itcnt, " iterations with relative error ", (norm(r)/norm(b)), ".")
    end

    return x
end

# For when you can't use BLAS.
# set to use similar paramaters similar to MATLAB
function cgSlow{Tval}(mat, b::Array{Tval,1};
        tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false)

    n = size(mat,2)

    x = zeros(Tval,n)

    tol = tol * norm(b)
    r = copy(b)
    p = copy(r)

    # the following would have higher accuracy
    #       rho = sum(r.^2)
    rho = norm(r)^2

    t1 = time()

    itcnt = 0
    while itcnt < maxits
        itcnt = itcnt+1

        q = mat*p

        al = rho/dot(p, q)

        # x = x + al * p
        for i in 1:n
            x[i] = x[i] + al*p[i]
        end

        # r = r - al*q
        for i in 1:n
            r[i] = r[i] - al*q[i]
        end

        if norm(r) < tol
            break
        end

        oldrho = rho
        rho = norm(r)^2

        # the following would have higher accuracy
        #       rho = sum(r.^2)

        beta = rho/oldrho

        # p = r + beta*p
        for i in 1:n
            p[i] = r[i] + beta*p[i]
        end

        if (time() - t1) > maxtime
            if verbose
                println("CG stopped at maxtime.")
            end
            break
        end
      end

    if verbose
        println("CG stopped after: ", itcnt, " iterations with relative error ", (norm(r)/norm(b)), ".")
    end

    return x
end

#====================================
   PCG Code
=====================================#

# uses BLAS.  As fast as Matlab's pcg.
# set to use similar parameters "
function pcgBLAS{Tval}(mat, b::Array{Tval,1}, pre;
        tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false)

    n = size(mat,2)

    x = zeros(Tval,n)

    tol = tol * norm(b)
    r = copy(b)
    z = pre(r)
    p = copy(z)

    rho = dot(r, z) # BLAS.dot does not seem faster

    t1 = time()

    itcnt = 0
    while itcnt < maxits
        itcnt = itcnt+1

        q = mat*p

        if dot(p,q) == 0
            break
        end

        al = rho/dot(p, q)

        BLAS.axpy!(al,p,x)  # x = x + al * p

        BLAS.axpy!(-al,q,r)  # r -= al*q

        if norm(r) < tol #Converged?
            break
        end

        # here is the top of the code in numerical templates

        z = pre(r)

        oldrho = rho
        rho = dot(z, r)

        # the following would have higher accuracy
        #       rho = sum(r.^2)

        beta = rho/oldrho

        p = z + beta*p
        # BLAS.scal!(n,beta,p,1) # p *= beta
        # BLAS.axpy!(1.0,z,p) # p += z

        if (time() - t1) > maxtime
            if verbose
                println("PCG stopped at maxtime.")
            end
            break
        end

      end

    if verbose
        println("PCG stopped after: ", itcnt, " iterations with relative error ", (norm(r)/norm(b)), ".")
    end

    return x
end

function pcgSlow{Tval}(mat, b::Array{Tval,1}, pre;
        tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false)

    n = size(mat,2)

    x = zeros(Tval,n)

    tol = tol * norm(b)
    r = copy(b)
    z = pre(r)
    p = copy(z)

    rho = dot(r, z)

    t1 = time()

    itcnt = 0
    while itcnt < maxits
        itcnt = itcnt+1

        q = mat*p

        al = rho/dot(p, q)

        # x = x + al * p
        for i in 1:n
            x[i] = x[i] + al*p[i]
        end

        # r = r - al*q
        for i in 1:n
            r[i] = r[i] - al*q[i]
        end

        if norm(r) < tol #Converged?
            break
        end

        # here is the top of the code in numerical templates

        z = pre(r)

        oldrho = rho
        rho = dot(z, r)

        # the following would have higher accuracy
        #       rho = sum(r.^2)

        beta = rho/oldrho

        # p = z + beta*p
        for i in 1:n
            p[i] = z[i] + beta*p[i]
        end

        if (time() - t1) > maxtime
            if verbose
                println("PCG stopped at maxtime.")
            end
            break
        end

    end

    if verbose
        println("PCG stopped after: ", itcnt, " iterations with relative error ", (norm(r)/norm(b)), ".")
    end

    return x
end

ext: jl


text:
CONIC_OPTIMIZERS = [SCS.Optimizer, ProxSDP.Optimizer, COSMO.Optimizer]

@testset "MOI to MatOI conversion 1" begin
    # _psd1test: https://github.com/jump-dev/MathOptInterface.jl/blob/master/src/Test/contconic.jl#L2417

    for optimizer in CONIC_OPTIMIZERS
        model = MOI.instantiate(optimizer, with_bridge_type=Float64)
        δ = √(1 + (3*√2+2)*√(-116*√2+166) / 14) / 2
        ε = √((1 - 2*(√2-1)*δ^2) / (2-√2))
        y2 = 1 - ε*δ
        y1 = 1 - √2*y2
        obj = y1 + y2/2
        k = -2*δ/ε
        x2 = ((3-2obj)*(2+k^2)-4) / (4*(2+k^2)-4*√2)
        α = √(3-2obj-4x2)/2
        β = k*α

        X = MOI.add_variables(model, 6)
        x = MOI.add_variables(model, 3)

        vov = MOI.VectorOfVariables(X)

        cX = MOI.add_constraint(
            model,
            MOI.VectorAffineFunction{Float64}(vov), MOI.PositiveSemidefiniteConeTriangle(3)
        )

        cx = MOI.add_constraint(
            model,
            MOI.VectorAffineFunction{Float64}(MOI.VectorOfVariables(x)), MOI.SecondOrderCone(3)
        )

        c1 = MOI.add_constraint(
            model,
            MOI.VectorAffineFunction(
                MOI.VectorAffineTerm.(1:1, MOI.ScalarAffineTerm.([1., 1., 1., 1.], [X[1], X[3], X[end], x[1]])),
                [-1.0]
            ),
            MOI.Zeros(1)
        )

        c2 = MOI.add_constraint(
            model,
            MOI.VectorAffineFunction(
                MOI.VectorAffineTerm.(1:1, MOI.ScalarAffineTerm.([1., 2, 1, 2, 2, 1, 1, 1], [X; x[2]; x[3]])),
                [-0.5]
            ),
            MOI.Zeros(1)
        )

        objXidx = [1:3; 5:6]
        objXcoefs = 2*ones(5)
        MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([objXcoefs; 1.0], [X[objXidx]; x[1]]), 0.0))
        MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

        MatModel = MatOI.get_conic_form(Float64, model, [cX; cx; c1; c2])

        @test MatModel.c' ≈ [2. 2. 2. 0. 2. 2. 1. 0. 0.]
        @test MatModel.b' ≈ [0.   0.   0.   0.   0.   0.   0.   0.   0.  -1.  -0.5 ]
        @test MatModel.A.nzval ≈ [-1.0, -1.0, -1.0, -1.0, -2.0, -1.0, -1.0, -1.0, -1.0, -2.0, -1.0, -2.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0] atol=ATOL rtol=RTOL
    end
end

@testset "MOI to MatOI conversion 2" begin
    # find equivalent diffcp program here - https://github.com/AKS1996/jump-gsoc-2020/blob/master/diffcp_sdp_3_py.ipynb

    for optimizer in CONIC_OPTIMIZERS
        model = MOI.instantiate(optimizer, with_bridge_type=Float64)

        x = MOI.add_variables(model, 7)
        @test MOI.get(model, MOI.NumberOfVariables()) == 7

        η = 10.0

        c1  = MOI.add_constraint(
            model,
            MOI.VectorAffineFunction(
                MOI.VectorAffineTerm.(1, MOI.ScalarAffineTerm.(-1.0, x[1:6])),
                [η]
            ),
            MOI.Nonnegatives(1)
        )
        c2 = MOI.add_constraint(model, MOI.VectorAffineFunction(MOI.VectorAffineTerm.(1:6, MOI.ScalarAffineTerm.(1.0, x[1:6])), zeros(6)), MOI.Nonnegatives(6))
        α = 0.8
        δ = 0.9
        c3 = MOI.add_constraint(model, MOI.VectorAffineFunction(MOI.VectorAffineTerm.([fill(1, 7); fill(2, 5);     fill(3, 6)],
                                                                MOI.ScalarAffineTerm.(
                                                                [ δ/2,       α,   δ, δ/4, δ/8,      0.0, -1.0,
                                                                    -δ/(2*√2), -δ/4, 0,     -δ/(8*√2), 0.0,
                                                                    δ/2,     δ-α,   0,      δ/8,      δ/4, -1.0],
                                                                [x[1:7];     x[1:3]; x[5:6]; x[1:3]; x[5:7]])),
                                                                zeros(3)), MOI.PositiveSemidefiniteConeTriangle(2))
        c4 = MOI.add_constraint(
            model,
            MOI.VectorAffineFunction(
                MOI.VectorAffineTerm.(1, MOI.ScalarAffineTerm.(0.0, [x[1:3]; x[5:6]])),
                [0.0]
            ),
            MOI.Zeros(1)
        )

        MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x[7])], 0.0))
        MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)

        MatModel = MatOI.get_conic_form(Float64, model, [c1; c2; c3; c4])

        @test MatModel.c ≈ [-0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -1.0]
        @test MatModel.b ≈ [10.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
        @test MatModel.A.nzval ≈ [1.0, -1.0, -0.45, 0.318198, -0.45, 1.0, -1.0, -0.8, 0.225, -0.1, 1.0, -1.0, -0.9, 1.0, -1.0, -0.225, 1.0, -1.0, -0.1125, 0.0795495, -0.1125, 1.0, -1.0, -0.225, 1.0, 1.0] atol=ATOL rtol=RTOL
    end
end

@testset "Testing minor utilities" begin
    # testing offsets
    n = rand(Int) % 100
    @test MatOI.cons_offset(MOI.Zeros(n)) == n
    @test MatOI.cons_offset(MOI.Nonnegatives(n)) == n
    @test MatOI.cons_offset(MOI.SecondOrderCone(n)) == n
    @test MatOI.cons_offset(MOI.PositiveSemidefiniteConeTriangle(n)) == ((n+1)*n)/2

    model = MOI.instantiate(SCS.Optimizer, with_bridge_type=Float64)
    cf = MatOI.get_conic_form(Float64,model,[])
    @test MOI.is_empty(cf) == false
    MOI.empty!(cf)
    @test MOI.is_empty(cf) == true
end


text:
import Automa
import Automa.RegExp: @re_str
const re = Automa.RegExp

keyword = re"break|const|continue|else|elseif|end|for|function|if|return|type|using|while"
identifier = re"[A-Za-z_][0-9A-Za-z_!]*"
operator = re"-|\+|\*|/|%|&|\||^|!|~|>|<|<<|>>|>=|<=|=>|==|==="
macrocall = re"@" * re"[A-Za-z_][0-9A-Za-z_!]*"
comment = re"#[^\r\n]*"
char = re.cat('\'', re"[ -&(-~]" | re.cat('\\', re"[ -~]"), '\'')
string = re.cat('"', re.rep(re"[ !#-~]" | re.cat("\\\"")), '"')
triplestring = re.cat("\"\"\"", (re"[ -~]*" \ re"\"\"\""), "\"\"\"")
newline = re"\r?\n"

const minijulia = Automa.compile(
    re","          => :(emit(:comma)),
    re":"          => :(emit(:colon)),
    re";"          => :(emit(:semicolon)),
    re"\."         => :(emit(:dot)),
    re"\?"         => :(emit(:question)),
    re"="          => :(emit(:equal)),
    re"\("         => :(emit(:lparen)),
    re"\)"         => :(emit(:rparen)),
    re"\["         => :(emit(:lbracket)),
    re"]"          => :(emit(:rbracket)),
    re"{"          => :(emit(:lbrace)),
    re"}"          => :(emit(:rbrace)),
    re"$"          => :(emit(:dollar)),
    re"&&"         => :(emit(:and)),
    re"\|\|"       => :(emit(:or)),
    re"::"         => :(emit(:typeannot)),
    keyword        => :(emit(:keyword)),
    identifier     => :(emit(:identifier)),
    operator       => :(emit(:operator)),
    macrocall      => :(emit(:macrocall)),
    re"[0-9]+"     => :(emit(:integer)),
    comment        => :(emit(:comment)),
    char           => :(emit(:char)),
    string         => :(emit(:string)),
    triplestring   => :(emit(:triplestring)),
    newline        => :(emit(:newline)),
    re"[\t ]+"     => :(emit(:spaces)),
)

#=
write("minijulia.dot", Automa.machine2dot(minijulia.machine))
run(`dot -Tsvg -o minijulia.svg minijulia.dot`)
=#

context = Automa.CodeGenContext()
@eval function tokenize(data)
    $(Automa.generate_init_code(context, minijulia))
    p_end = p_eof = sizeof(data)
    tokens = Tuple{Symbol,String}[]
    emit(kind) = push!(tokens, (kind, data[ts:te]))
    while p ≤ p_eof && cs > 0
        $(Automa.generate_exec_code(context, minijulia))
    end
    if cs < 0
        error("failed to tokenize")
    end
    return tokens
end

tokens = tokenize("""
quicksort(xs) = quicksort!(copy(xs))
quicksort!(xs) = quicksort!(xs, 1, length(xs))

function quicksort!(xs, lo, hi)
    if lo < hi
        p = partition(xs, lo, hi)
        quicksort!(xs, lo, p - 1)
        quicksort!(xs, p + 1, hi)
    end
    return xs
end

function partition(xs, lo, hi)
    pivot = div(lo + hi, 2)
    pvalue = xs[pivot]
    xs[pivot], xs[hi] = xs[hi], xs[pivot]
    j = lo
    @inbounds for i in lo:hi-1
        if xs[i] <= pvalue
            xs[i], xs[j] = xs[j], xs[i]
            j += 1
        end
    end
    xs[j], xs[hi] = xs[hi], xs[j]
    return j
end
""")


text:
# Constants

const AMD_QUEUE_ALIGN_BYTES = 64
const AMD_SIGNAL_ALIGN_BYTES = 64
const API_TABLE_MAJOR_VERSION = 0x01
const CORE_API_TABLE_MAJOR_VERSION = 0x01
const AMD_EXT_API_TABLE_MAJOR_VERSION = 0x01
const FINALIZER_API_TABLE_MAJOR_VERSION = 0x01
const IMAGE_API_TABLE_MAJOR_VERSION = 0x01
const AQLPROFILE_API_TABLE_MAJOR_VERSION = 0x01
const API_TABLE_STEP_VERSION = 0x00
const CORE_API_TABLE_STEP_VERSION = 0x00
const AMD_EXT_API_TABLE_STEP_VERSION = 0x00
const FINALIZER_API_TABLE_STEP_VERSION = 0x00
const IMAGE_API_TABLE_STEP_VERSION = 0x00
const AQLPROFILE_API_TABLE_STEP_VERSION = 0x00
const AMD_INTERFACE_VERSION_MAJOR = 1
const AMD_INTERFACE_VERSION_MINOR = 0
const AQLPROFILE_VERSION_MAJOR = 2
const AQLPROFILE_VERSION_MINOR = 0
const AMD_AQLPROFILE_VERSION_MAJOR = 1
#const AQLProfile_Pfn = AQLProfile_1_00_pfn
const VenAMD_loader = 1000


text:
module IntervalLinearAlgebra

using StaticArrays, Requires, Reexport
using LinearAlgebra: checksquare

import Base: *
import CommonSolve: solve
import IntervalArithmetic: mid

function  __init__()
    # @require IntervalConstraintProgramming = "138f1668-1576-5ad7-91b9-7425abbf3153" include("linear_systems/oettli_nonlinear.jl")
    @require LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043" include("linear_systems/oettli_linear.jl")

    set_multiplication_mode(config[:multiplication])
end

@reexport using LinearAlgebra, IntervalArithmetic

const IA = IntervalArithmetic

export
    set_multiplication_mode, get_multiplication_mode,
    LinearKrawczyk, Jacobi, GaussSeidel, GaussianElimination, HansenBliekRohn, NonLinearOettliPrager, LinearOettliPrager,
    NoPrecondition, InverseMidpoint, InverseDiagonalMidpoint,
    solve, enclose, epsilon_inflation,
    comparison_matrix, interval_norm, interval_isapprox, Orthants,
    is_H_matrix, is_strongly_regular, is_strictly_diagonally_dominant, is_Z_matrix, is_M_matrix,
    rref,
    eigenbox, Rohn, Hertz, verify_eigen, bound_perron_frobenius_eigenvalue

include("linear_systems/enclosures.jl")
include("linear_systems/precondition.jl")
include("linear_systems/solve.jl")
include("linear_systems/verify.jl")
include("linear_systems/oettli.jl")
include("multiplication.jl")
include("utils.jl")
include("classify.jl")
include("rref.jl")

include("eigenvalues/interval_eigenvalues.jl")
include("eigenvalues/verify_eigs.jl")
end


text:
"""
    const TargetTree{N,T}

Type alias for a `ClusterTree` of points (represented as `SVector{N,T}`)
containing no additional data. See the documentation of
`ClusterTree` in [`WavePropBase`](https://github.com/WaveProp/WavePropBase.jl)
for more information.
"""
const TargetTree{N,T} = ClusterTree{Vector{SVector{N,T}},HyperRectangle{N,T},Nothing}

"""
    initialize_target_tree(;points,splitter)

Create the tree-structure for clustering the target `points` using the splitting
strategy of `splitter`. Returns a `TargetTree`.
"""
function initialize_target_tree(;Xpoints,splitter)
    target_tree = ClusterTree(Xpoints,splitter)
    return target_tree
end

max_stars_repo_name: WaveProp/IFGF

max_issues_repo_name: WaveProp/IFGF

max_forks_repo_name: WaveProp/IFGF


text:
"""
    simulate(mdl)

    Simulation function

    # Arguments

    # Examples
"""
function simulate(mdl, data::StructVector; feedback = x -> missing, init_θ=mdl())
    initial_state = data.s[1]
    initial_hidden = data.h[1]
    state_transition(history) = data.s[min(length(history)+1, length(data))]
    hidden_transition(history) = data.h[min(length(history)+1, length(data))]
    ending_condition(history) = length(history) > length(data)

    return simulate(mdl;
            state_transition = state_transition,
            hidden_transition = hidden_transition,
            feedback = feedback,
            initial_state = initial_state,
            initial_hidden = initial_hidden,
            ending_condition = ending_condition,
            init_θ = init_θ)
end

function simulate(mdl;
    state_transition = x -> 1,
    feedback = x -> missing,
    hidden_transition = x -> missing,
    initial_state = 1,
    initial_hidden = missing,
    ending_condition = x -> length(x) > 100,
    init_θ = mdl())

    # Initialize
    θ = init_θ
    P_ = AnimalBehavior.observ(mdl, initial_state; θ...)
    a_type = Distributions.eltype(typeof(P_))
    r_type = Base.return_types(feedback, (StructVector,))[1]

    history = StructVector(s=[initial_state], a=Vector{Any}([missing]), r=Vector{Any}([missing]), h=Vector{Any}([initial_hidden]))
    latent = StructVector([deepcopy(θ)])

    while !ending_condition(history)
        # current state
        s = history[end].s

        # action
        P = AnimalBehavior.observ(mdl, s; θ...)
        a = rand(P)
        history.a[end] = a

        # feedback
        r = feedback(history)
        history.r[end] = r

        # update
        AnimalBehavior.evol!(mdl, s, a, r; θ...)

        # next state
        ns = state_transition(history)
        nh = hidden_transition(history)
        push!(history, (s=ns, a=missing, r=missing, h=nh))
        push!(latent, deepcopy(θ))
    end
    data = StructVector(s = history[1:end-1].s,
                        a = convert.(a_type,history[1:end-1].a),
                        r = convert.(r_type, history[1:end-1].r),
                        h = history[1:end-1].h)
    return Simulation(mdl.name, data, latent[1:end-1])
end

max_stars_repo_name: sqwayer/AnimalBehavior

max_issues_repo_name: sqwayer/AnimalBehavior


text:
"""
    GreaterOrEqualConstant(x::SeaPearl.AbstractIntVar, v::Int)

Inequality constraint, `x >= v`
"""
struct GreaterOrEqualConstant <: OnePropagationConstraint
    x       ::AbstractIntVar
    v       ::Int
    active  ::StateObject{Bool}
    function GreaterOrEqualConstant(x, v, trailer)
        constraint = new(x, v, StateObject(true, trailer))
        addOnDomainChange!(x, constraint)
        return constraint
    end
end

"""
    propagate!(constraint::GreaterOrEqualConstant, toPropagate::Set{Constraint}, prunedDomains::CPModification)

`GreaterOrEqualConstant` propagation function. Basically remove the values above `v`.
"""
function propagate!(constraint::GreaterOrEqualConstant, toPropagate::Set{Constraint}, prunedDomains::CPModification)
    setValue!(constraint.active, false)

    addToPrunedDomains!(prunedDomains, constraint.x, removeBelow!(constraint.x.domain, constraint.v))
    triggerDomainChange!(toPropagate, constraint.x)
    return !isempty(constraint.x.domain)
end

function Base.show(io::IO, ::MIME"text/plain", con::GreaterOrEqualConstant)
    println(io, typeof(con), ": ", con.x.id, " ≥ ", con.v, ", active = ", con.active)
    print(io, "   ", con.x)
end

function Base.show(io::IO, con::GreaterOrEqualConstant)
    print(io, typeof(con), ": ", con.x.id, " ≥ ", con.v)
end


text:
using Documenter, GigaSOM

makedocs(modules = [GigaSOM],
        clean = false,
        format = Documenter.HTML(prettyurls = !("local" in ARGS),
                canonical = "https://lcsb-biocore.github.io/GigaSOM.jl/stable/",
                assets = ["assets/gigasomlogotransp.ico"]),
        sitename = "GigaSOM.jl",
        authors = "The developers of GigaSOM.jl",
        linkcheck = !("skiplinks" in ARGS),
        pages = [
                "Home" => "index.md",
                "Background" => "background.md",
                "How to get started" => "howToGetStarted.md",
                "Functions" => "functions.md",
                "How to contribute" => "howToContribute.md",
                ],
        )

deploydocs(
    repo = "github.com/LCSB-BioCore/GigaSOM.jl.git",
    target = "build",
    branch = "gh-pages",
    devbranch = "origin/develop",
    versions = "stable" => "v^",
    )


text:
mutable struct KDensity{R<:AbstractRange}
    "Gridpoints for evaluating the density."
    x::R
    "Kernel density at corresponding gridpoints `x`."
    density::Vector{Float64}
end

function KDensity(v)
    k = KernelDensity.kde(v)
    probability = k.density .* step(k.x)
    KDensity(k.x, probability)
end

function trim_dist(vec; p_limit = 0.01)
    v = jump_missing(vec)
    k = kde(v)
    ik = InterpKDE(k)
    P = [ismissing(val) ? 1 : pdf(ik, val).*step(k.x) for val in vec]
    [x > p_limit for x in P]
end

function trim_dist(df::AbstractDataFrame,x::Symbol; p_limit = 0.01)
    v = df[:,x]
    idxs = trim_dist(v; p_limit = p_limit)
    df[idxs,:]
end

function trim_dist!(df::AbstractDataFrame,x::Symbol; p_limit = 0.01)
    v = df[:,x]
    idxs = trim_dist(v; p_limit = p_limit)
    delete!(df,Not(idxs))
end


text:
using ReachabilityAnalysis

@taylorize function lotkavolterra!(dx, x, params, t)
    local α, β, γ, δ = 1.5, 1., 3., 1.
    dx[1] = α * x[1] - β * x[1] * x[2]
    dx[2] = δ * x[1] * x[2] - γ * x[2]
    return dx
end

X0 = Hyperrectangle(low=[4.8, 1.8], high=[5.2, 2.2]);
prob = @ivp(x' = lotkavolterra!(x), dim=2, x(0) ∈ X0)
sol = solve(prob, T=8.0, alg=TMJets());
solz = overapproximate(sol, Zonotope);

using Plots

plot(solz, vars=(1, 2), alpha=0.3,lw=0., xlab="x", ylab="y", label="Flowpipe", legend=:bottomright)
plot!(X0, label="X(0)")

@taylorize function f(du, u, p, t)
    du[1] = u[3] * u[1] - u[4] * (u[1] * u[2]) - u[7] * u[1]^2
    du[2] = -u[5] * u[2] + u[6] * (u[1] * u[2])

    #encode uncertain params
    du[3] = zero(u[1]) # p[1]
    du[4] = zero(u[1]) # p[2]
    du[5] = zero(u[1]) # p[3]
    du[6] = zero(u[1]) # p[4]
    du[7] = zero(u[1]) # p[5]
end

p_int = (0.99..1.01) × (0.99..1.01) × (2.99..3.01) × (0.99..1.01) × (0.099..0.101)
U0 = Singleton([1.0, 1.0]) × convert(Hyperrectangle, p_int)
prob = @ivp(u' = f(u), dim: 7, u(0) ∈ U0);

sol = solve(prob, tspan=(0.0, 10.0));
solz = overapproximate(sol, Zonotope);

plot(solz, vars=(1, 2), lw=0.3, title="Uncertain params", lab="abstol = 1e-15", xlab="u1", ylab="u2")

u0 = Singleton([1.0, 1.0])
□(ϵ) = BallInf(zeros(2), ϵ)
U0 = (u0 ⊕ □(0.05)) × convert(Hyperrectangle, p_int)

prob = @ivp(u' = f(u), dim: 7, u(0) ∈ U0)

sol = solve(prob, tspan=(0.0, 10.0), TMJets(abstol=1e-10))
solz = overapproximate(sol, Zonotope)
plot(solz, vars=(1, 2), color=:orange, lw=0.3,
     lab="eps = 0.05", title="Uncertain u0 and uncertain params",
     xlab="u1", ylab="u2")

U0 = (u0 ⊕ □(0.01)) × convert(Hyperrectangle, p_int)
prob = @ivp(u' = f(u), dim: 7, u(0) ∈ U0)

sol = solve(prob, tspan=(0.0, 10.0), TMJets(abstol=1e-10))
solz = overapproximate(sol, Zonotope)
plot!(solz, vars=(1, 2), color=:blue, lw=0.3,
  lab="eps = 0.01", title="Uncertain u0 and uncertain params",
  xlab="u1", ylab="u2")


text:
module TestTikZ

using Test
using Catlab.Graphics.TikZ

# Pretty-print
##############

spprint(expr::Expression) = sprint(pprint, expr)

# Properties
@test spprint(Property("circle")) == "circle"
@test spprint(Property("fill","red")) == "fill=red"
@test spprint(Property("style",[Property("circle"),Property("fill","red")])) ==
  "style={circle,fill=red}"

# Node statement
@test spprint(Node("f")) == "\\node (f) {};"
@test spprint(Node("f"; props=[Property("circle")])) ==
  "\\node[circle] (f) {};"
@test spprint(Node("f"; props=[Property("circle"),Property("fill","red")])) ==
  "\\node[circle,fill=red] (f) {};"
@test spprint(Node("f"; coord=Coordinate("0","1"))) ==
  "\\node (f) at (0,1) {};"
@test spprint(Node("f"; coord=Coordinate("1","0"), content="fun")) ==
  "\\node (f) at (1,0) {fun};"
@test spprint(Node(; props=[Property("circle")])) == "\\node[circle] {};"

# Edge statement
const PathOp = PathOperation
@test spprint(Edge("f", PathOp("to"), "g")) == "\\path (f) to (g);"
@test spprint(Edge("f", PathOp("to"), "g"; props=[Property("red")])) ==
  "\\path[red] (f) to (g);"
@test spprint(Edge("f", PathOp("to"), EdgeNode(), "g")) ==
  "\\path (f) to node (g);"
@test spprint(Edge("f", PathOp("to"),
              EdgeNode(; props=[Property("circle")], content="e"), "g")) ==
  "\\path (f) to node[circle] {e} (g);"
@test spprint(Edge("f", PathOp("to"; props=[Property("out","0")]), "g")) ==
  "\\path (f) to[out=0] (g);"
@test spprint(Edge("f", PathOp("--"), "g", PathOp("--"), "h")) ==
  "\\path (f) -- (g) -- (h);"

# Picture statement
@test spprint(Picture(Node("f"), Node("g"), Edge("f",PathOp("to"),"g"))) == """
\\begin{tikzpicture}
  \\node (f) {};
  \\node (g) {};
  \\path (f) to (g);
\\end{tikzpicture}"""

@test spprint(Picture(Node("f");
  props=[Property("node distance","1cm")])) == """
\\begin{tikzpicture}[node distance=1cm]
  \\node (f) {};
\\end{tikzpicture}"""

# Document statement
@test spprint(Document(
  Picture(Edge("f",PathOp("--"),"g"; props=[Property(">-Stealth")])),
  libraries=["arrows.meta"])) == """
\\usetikzlibrary{arrows.meta}
\\begin{tikzpicture}
  \\path[>-Stealth] (f) -- (g);
\\end{tikzpicture}"""

# Scope statement
@test spprint(
Picture(
  Scope(Node("f"); props=[Property("red")]),
  Scope(Node("g"); props=[Property("blue")]),
  Edge("f",PathOp("to"),"g")
)) == """
\\begin{tikzpicture}
  \\begin{scope}[red]
    \\node (f) {};
  \\end{scope}
  \\begin{scope}[blue]
    \\node (g) {};
  \\end{scope}
  \\path (f) to (g);
\\end{tikzpicture}"""

# Graph statements
@test spprint(
Graph(
  GraphNode("f";content="fun"),
  GraphEdge("f","g")
)) == """
\\graph{
  f/\"fun\";
  f -> g;
};"""

@test spprint(
Graph(
  GraphNode("f";props=[Property("circle")]),
  GraphEdge("f","g";props=[Property("edge label","X")])
)) == """
\\graph{
  f [circle];
  f ->[edge label=X] g;
};"""

@test spprint(
Graph(
  GraphScope(
    GraphEdge("a","b"); props=[Property("same layer")]
  ),
  GraphScope(
    GraphEdge("c","d"); props=[Property("same layer")]
  )
)) == """
\\graph{
  {[same layer]
    a -> b;
  };
  {[same layer]
    c -> d;
  };
};"""

# Matrix statement
@test spprint(
MatrixNode(
  [ [[Node("f")]]              [[Edge("g1",PathOp("to"),"g2")]];
    [[Node("h1"),Node("h2")]]  [[Node("i1"),Node("i2")]] ];
  props=[Property("draw","red")]
)) == """
\\matrix[draw=red]{
  \\node (f) {}; &
  \\path (g1) to (g2); \\\\
  \\node (h1) {};
  \\node (h2) {}; &
  \\node (i1) {};
  \\node (i2) {}; \\\\
};"""

end


text:
export tabs, pages, menu, submenu, item, icon, iconbutton, toolbar
# Higher-order layouts: e.g. tabs, pages

abstract Selection <: Tile

for (fn, typ) in  [(:pages, :Pages), (:tabs, :Tabs), (:menu, :Menu)]
    @eval begin
        immutable $typ <: Selection
            tiles::AbstractArray
            multi::Bool
            selected::Any
        end
        $fn(tiles; selected=0, multi=false, name=:_tabs) =
            selectable($typ(tiles, multi, selected), name=name)
    end
end

render(ps::Pages) =
    render(ps.tiles, "core-pages") & [:selected=>ps.selected, :multi=>ps.multi]

render(tabs::Tabs) =
    Elem("paper-tabs",
        map(t -> Elem("paper-tab", render(t)), tabs.tiles),
        selected=tabs.selected, multi=tabs.multi)

render(m::Menu) =
    render(m.tiles, "core-menu") & [:selected=>m.selected, :multi=>m.multi]

immutable SubMenu <: Tile
    icon::String
    label::String
    tiles::AbstractArray
end
submenu(icon, label, tiles) = SubMenu(icon, label, tiles)
submenu(label, tiles) = SubMenu(icon, label, tiles)

render(m2::SubMenu) =
    render(m2.tiles, "core-submenu")

immutable Item <: Tile
    icon::String
    label::String
end

item(icon, label) = Item(icon, label)
item(label) = Item("", label)

render(i::Item) =
    Elem("core-item", icon=i.icon, label=i.label)

# Toolbar and Icons

@api toolbar => Toolbar <: Tile begin
    curry(tiles::AbstractArray)
end
render(t::Toolbar) =
    render(t.tiles, "core-toolbar")

@api icon => Icon <: Tile begin
    typedarg(icon::String="")
    kwarg(url::Bool=false)
end
render(i::Icon) =
    Elem("core-icon") & [(i.url ? :src : :icon) => i.icon]

@api iconbutton => IconButton <: Tile begin
    typedarg(icon::String="")
    kwarg(url::Bool=false)
end
render(i::IconButton) =
    Elem("paper-icon-button") & [(i.url ? :src : :icon) => i.icon]


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

using Markdown
using InteractiveUtils

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

# ╔═╡ ba447bcc-0094-11eb-068b-6b5e045fd602
# Packages for basic plotting
begin
	using Pkg
	Pkg.activate(".")
	#Pkg.develop("EchelleInstruments")
	#Pkg.instantiate()
	using RvSpectMLBase
	import EchelleInstruments
	import EchelleInstruments: EXPRES
	using DataFrames, Query
	using Statistics
	using Plots
	using PlutoUI
end

# ╔═╡ b29da99e-0142-11eb-3504-15cb11acd16f
# Packages for fitting GPs
begin
	using LinearAlgebra
	using PDMats
	using StaticArrays
	using Stheno, TemporalGPs

	import Stheno: AbstractGP
	using Distributions
	import Distributions: AbstractMvNormal
	const Fx_PosteriorType = Distribution{Multivariate,Continuous}
end;

# ╔═╡ 78bc57ea-0140-11eb-2f53-15d0cb1e8b7b
md"# Visualize GP fit to spectra"

# ╔═╡ 19fc60d6-009a-11eb-3d9c-0f9152ece864
md"Divide by blaze: $(@bind div_blaze CheckBox(default=true))    Show tellurics? $(@bind plt_tellurics CheckBox())  Show continuum? $(@bind plt_continuum CheckBox())"
# Log y-axis? $(@bind log_y CheckBox()) # Need to fix GP fitting to logs

# ╔═╡ 26060010-0149-11eb-0863-33bfb7908253
md"Scaling factors for GP variance: $(@bind smooth_factor_var Slider(0.05:0.05:5; default=1.0, show_value=true))   & length scale $(@bind smooth_factor_length Slider(0.05:0.05:4; default=1.0, show_value=true))"

# ╔═╡ 9d4a5640-00a3-11eb-11f7-5d5e9c00fc61
md"### Boring code below"

# ╔═╡ ce23b456-00a5-11eb-0d07-2b15dfff6257
md"### Find files to read data from"

# ╔═╡ c5b8c394-013d-11eb-2353-7da0efded8f9
begin
	expres_data_path = joinpath(homedir(),"Data/EXPRES/inputs")
	target_subdir = "101501"
 	fits_target_str = "101501"
 	output_dir = "examples/output/"
    df_files = EchelleInstruments.make_manifest(expres_data_path, target_subdir, EXPRES );
end;

# ╔═╡ c364ee3e-0095-11eb-0faa-656ede64325f
begin
	# Configuration from param.in
	linelist_for_ccf_filename = "G8.espresso.mas"
	ccf_mid_velocity = -5e3
	max_spectra_to_use = 10
	df_files_use = df_files |>
	    @filter( _.target == fits_target_str ) |>
	    @take(max_spectra_to_use) |>
	    DataFrame;
end;

# ╔═╡ cfabe62c-0098-11eb-3b6d-ed7ab2f9f0f1
md"Select the observation index to plot: $(@bind obs_idx Slider(1:size(df_files_use,1); default=1, show_value=true))"

# ╔═╡ 99e1c70c-00a5-11eb-3cc9-fb8934931026
md"### Read data files & extract each order"

# ╔═╡ 01c2dfa6-0096-11eb-1361-af95574e6f95
begin
	spectra = map(row->EXPRES.read_data(row), eachrow(df_files_use))
end;

# ╔═╡ 1a16d0ec-0095-11eb-3bac-0ddaaf387169
function extract_orders(data::AbstractVector{ST}; verbose::Bool = false, orders_to_use = RvSpectMLBase.orders_to_use_default(first(data).inst) ) where { ST<:AbstractSpectra }
      order_list_timeseries = make_order_list_timeseries(data, orders_to_use=orders_to_use)
      order_list_timeseries = filter_bad_chunks(order_list_timeseries,verbose=true)
end;

# ╔═╡ fdcf688a-0091-11eb-3fb4-9333baf67eab
 order_list_timeseries = extract_orders(spectra);

# ╔═╡ a6031a34-0098-11eb-38aa-891004881de9
md"Select a chunk ID: $(@bind chunk_idx Slider(1:num_chunks(order_list_timeseries); default=floor(Int,num_chunks(order_list_timeseries)//2), show_value=true))"

# ╔═╡ 4ab1c184-009b-11eb-3438-35502bf52add
begin
	pixel_range_min = 1
	pixel_range_max = maximum(map(chid->length(order_list_timeseries[1][chid].λ),1:num_chunks(order_list_timeseries)))
	max_width_pixels = floor(Int,(pixel_range_max)//2)
	#pixel_range_max = min(length(order_list_timeseries[obs_idx][chunk_idx].λ), length(order_list_timeseries[obs_jdx][chunk_idx].λ))
	#max_width_pixels = floor(Int,minimum(map(chid->length(order_list_timeseries[1][chid].λ),1:num_chunks(order_list_timeseries)))//2)
	md"Central pixel: $(@bind mid_pixel Slider(pixel_range_min:pixel_range_max; default=floor(Int,(pixel_range_min+pixel_range_max)//2)))  and
	log2(Half-width in pixels): $(@bind log2_halfwidth_pixels Slider(0:0.2:log(max_width_pixels)/log(2); default=round(log(max_width_pixels)/log(2)*10)/10, show_value=true))"
end

# ╔═╡ 98c79090-014a-11eb-1ea9-c1b334406a3f
md"### Transform observed data"

# ╔═╡ 5b4656e8-014a-11eb-245f-65461ed5c02c
md"### Code for fitting GPs"

# ╔═╡ 5279c31c-0143-11eb-1e2c-bf94e52c0eb1
# Code to create prior GP
begin
	#gp_param_default = [ 0.4890909216856761, 5.800274590507981e-5] .* smooth_factor   # TODO Generalize.  Values fr fit to one order of one EXPRES spectra
	#gp_param_default = [ 0.4890909216856761, 5.800274590507981e-5 * smooth_factor ]
	gp_param_default = [ 0.4890909216856761 * smooth_factor_var, 5.800274590507981e-5 * smooth_factor_length ]
	gp_σ², gp_l = gp_param_default
	gp_kernel = gp_σ² * stretch(Matern52(), 1 / gp_l)
	prior_gp_naive = GP(gp_kernel, GPC())
	#f = to_sde(f_naive)   # if develop issues with StaticArrays could revert to this
	prior_gp = to_sde(prior_gp_naive, SArrayStorage(Float64))
end;

# ╔═╡ 6e3b4a28-0150-11eb-3f5c-59a21de3d298
# GP code assumes linear in flux.
log_y = false;

# ╔═╡ 44bce5e8-0147-11eb-18be-4df6dab237de
begin
	# Make copy spectral chunks to plot, since many need to manipulate for plotting
	(pixels1, order1) = order_list_timeseries[obs_idx][chunk_idx].flux.indices
	flux1 = copy(order_list_timeseries[obs_idx][chunk_idx].flux)
	sigmasq_obs_trans = copy(order_list_timeseries[obs_idx][chunk_idx].var)
	# Either divide flux by blaze function or multiply continue by blaze
	blaze1 = order_list_timeseries.metadata[obs_idx][:blaze][pixels1,order1]
	if div_blaze
		flux1 ./= blaze1
		sigmasq_obs_trans ./= blaze1.^2
	end
	#flux1_norm = mean(flux1)
	#flux1 ./= flux1_norm
	#sigmasq_obs_trans ./= flux1_norm.^2
	if log_y
		sigmasq_obs_trans ./= flux1.^2
		flux1 = log10.(flux1)
	end
end;

# ╔═╡ d21e73e8-0147-11eb-0a9c-cd854964a3ee
# Compute GP posterior for entire order
begin
	flux1_norm = mean(flux1)
	flux_post = posterior(prior_gp(log.(order_list_timeseries[obs_idx][chunk_idx].λ), sigmasq_obs_trans./flux1_norm.^2), flux1./flux1_norm.-1)
end;

# ╔═╡ 4edfca32-0092-11eb-25d5-dbda1e798503
begin
	plt = plot(size=(1100,500))
	if log_y  ylabel!("log Flux")   else ylabel!("Flux") end
	title!("Compare spectra to GP model")

	# Choose which pixels to plot
	halfwidth_pixels = floor(Int,2^log2_halfwidth_pixels)
	min_pixel = min(max(1,mid_pixel-halfwidth_pixels),length(order_list_timeseries[obs_idx][chunk_idx].flux))
	max_pixel = max(min(mid_pixel+halfwidth_pixels, length(order_list_timeseries[obs_idx][chunk_idx].flux)),1)

	xpred = order_list_timeseries[obs_idx][chunk_idx].λ
	pred_distr = marginals(flux_post(log.(xpred)))[min_pixel:max_pixel]
	mean_ypred = (mean.(pred_distr).+1.0)*flux1_norm
	sigma_ypred = sqrt.(var.(pred_distr)).*flux1_norm
	plot!(xpred[min_pixel:max_pixel], mean_ypred, color=:blue, label=:none )
	plot!(xpred[min_pixel:max_pixel], mean_ypred.+sigma_ypred, color=:cyan, linestyle=:dash, label=:none )
	plot!(xpred[min_pixel:max_pixel], mean_ypred.-sigma_ypred, color=:cyan, linestyle=:dash, label=:none )

	if plt_continuum
		continuum1 = copy(order_list_timeseries.metadata[obs_idx][:continuum][pixels1,order1])
		#continuum1 ./= flux1_norm
		if !div_blaze
			continuum1 .*= blaze1
		end
		if log_y
			continuum1 = log.(continuum1)
		end
		plot!(plt,order_list_timeseries[obs_idx][chunk_idx].λ[min_pixel:max_pixel],continuum1[min_pixel:max_pixel], color=:grey, label=:none)
	end

	# Plot observations
	plot!(plt,order_list_timeseries[obs_idx][chunk_idx].λ[min_pixel:max_pixel],flux1[min_pixel:max_pixel], color=:black, label=:none) # string(obs_idx))

	if plt_tellurics
		tellurics1 = order_list_timeseries.metadata[obs_idx][:tellurics][pixels1,order1]
		tellurics_idx1 = findall(x->x!=one(x),tellurics1[min_pixel:max_pixel]) .+ (min_pixel-1)
		scatter!(plt,order_list_timeseries[obs_idx][chunk_idx].λ[tellurics_idx1],flux1[tellurics_idx1], markersize=1.2, color=:green, label=:none)
	end

	plt2 = plot(size=(1100,500))
	plot!(xpred[min_pixel:max_pixel], flux1[min_pixel:max_pixel].-mean_ypred, color=:blue, label=:none )
	ylabel!("Residuals")
	xlabel!("λ (Å)")
	plot(plt,plt2,layout=(2,1),size=(1100,1000))
end

# ╔═╡ Cell order:
# ╟─78bc57ea-0140-11eb-2f53-15d0cb1e8b7b
# ╟─cfabe62c-0098-11eb-3b6d-ed7ab2f9f0f1
# ╟─a6031a34-0098-11eb-38aa-891004881de9
# ╟─4ab1c184-009b-11eb-3438-35502bf52add
# ╟─19fc60d6-009a-11eb-3d9c-0f9152ece864
# ╟─26060010-0149-11eb-0863-33bfb7908253
# ╠═4edfca32-0092-11eb-25d5-dbda1e798503
# ╟─9d4a5640-00a3-11eb-11f7-5d5e9c00fc61
# ╟─ba447bcc-0094-11eb-068b-6b5e045fd602
# ╟─ce23b456-00a5-11eb-0d07-2b15dfff6257
# ╟─c5b8c394-013d-11eb-2353-7da0efded8f9
# ╟─c364ee3e-0095-11eb-0faa-656ede64325f
# ╟─99e1c70c-00a5-11eb-3cc9-fb8934931026
# ╟─01c2dfa6-0096-11eb-1361-af95574e6f95
# ╟─1a16d0ec-0095-11eb-3bac-0ddaaf387169
# ╟─fdcf688a-0091-11eb-3fb4-9333baf67eab
# ╟─98c79090-014a-11eb-1ea9-c1b334406a3f
# ╟─44bce5e8-0147-11eb-18be-4df6dab237de
# ╟─5b4656e8-014a-11eb-245f-65461ed5c02c
# ╟─b29da99e-0142-11eb-3504-15cb11acd16f
# ╟─5279c31c-0143-11eb-1e2c-bf94e52c0eb1
# ╟─6e3b4a28-0150-11eb-3f5c-59a21de3d298
# ╟─d21e73e8-0147-11eb-0a9c-cd854964a3ee


text:
module RPSAgents
using Random
import Pkg
Pkg.add("StatsBase")
Pkg.add("Distributions")
using StatsBase, Distributions
export RandomAgent, OnlyScissorAgent, OnlyRockAgent,OnlyPaperAgent,CopyAgent, FreqCountingAgent, CounterLastActionAgent, ThompsonSamplingAgent, MarkovAgent

function RandomAgent()
  # Agent chơi bằng cách random sự lựa chọn
  function run(observation)
    return rand((0, 1, 2))
  end
  () -> (run)
end

function CopyAgent()
  # Agent chơi bằng cách copy hành động của đối thủ ở lượt trước đó
  function run(observation)
    if(observation["step"] > 0)
      return observation["lastOpponentAction"]
    end
    return rand((0, 1, 2))
  end
  () -> (run)
end

function CounterLastActionAgent()
  # Agent chơi bằng cách khắc chế hành động của đối thủ ở lượt trước đó
  function run(observation)
    return (observation["lastOpponentAction"] + 1) % 3
  end
  () -> (run)
end

function OnlyRockAgent()
  function run(observation)
    return 1
  end
  () -> (run)
end

function OnlyScissorAgent()
  function run(observation)
    return 0
  end
  () -> (run)
end

function OnlyPaperAgent()
  function run(observation)
    return 2
  end
  () -> (run)
end

function FreqCountingAgent()
  # Agent chơi bằng cách khắc chế lựa chọn mà đối thủ chọn nhiều nhất
  freqs = [1, 1, 1]

  function run(observation)
    if(observation["step"] == 0)
      return rand((0, 1, 2))
    end
    pred_action = findall(action->action == maximum(freqs), freqs)[1] - 1
    return (pred_action + 1) % 3
  end
  () -> (run)
end

function MarkovAgent()
  # Agent chơi bằng cách khắc chế pattern mà đối thủ chọn nhiều nhất
  table = Dict()
  action_seq = []

  function run(observation)
    k = 2
    # In early steps we can choose randomly option
    if(observation["step"] <= 2*k +1)
      action = rand((0,1,2))
      if(observation["step"] > 0)
        append!(action_seq, [observation["lastOpponentAction"], action])
      else
        append!(action_seq, action)
      end
      return action
    end
    key = join([i for i in action_seq[1: length(action_seq)-1]], "")
    if(haskey(table, key))
      table[key][observation["lastOpponentAction"]+1] += 1
    else
      table[key] = [1, 1, 1]
    end

    # Update action sequence
    deleteat!(action_seq, [1, 2])

    if(observation["step"] < 20)
      weights = table[key] ./ sum(table[key])
      pred_action = sample([0, 1, 2], Weights(weights))
    else
    pred_action = findall(action->action == maximum(table[key]), table[key])[1] - 1
    end
    counter_action = (pred_action + 1)%3
    append!(action_seq, [observation["lastOpponentAction"], counter_action])
    return counter_action
  end
  () -> (run)
end

function ThompsonSamplingAgent()
  decay_rate = 1.05
  reward = 3
  # Lưu số lượt thắng thua của từng agent
  bandit_state = Dict()
  # Lưu những thông tin về những lượt chơi
  history = []
  # Chứa nhưng agent đại diện cho từng chiến thuật chơi
  agents = Dict(
    "random_agent" => RandomAgent(),
    "only_scissor_agent" => OnlyScissorAgent(),
    "only_rock_agent" => OnlyRockAgent(),
    "only_paper_agent" => OnlyPaperAgent(),
    "copy_agent" => CopyAgent(),
    "freq_counting_agent" => FreqCountingAgent(),
    "counter_last_action_agent" => CounterLastActionAgent()
  )
  # Đếm số lần được lựa chọn của mỗi agent
  agent_count = Dict(
    "random_agent" => 0,
    "only_scissor_agent" => 0,
    "only_rock_agent" => 0,
    "only_paper_agent" => 0,
    "copy_agent" => 0,
    "freq_counting_agent" => 0,
    "counter_last_action_agent" => 0
  )

  # Lưu những thông tin của lượt hiện tai
  function store_current_state(history, observation, agent_name)
    agent_count[agent_name] += 1;
    state = Dict("step" => observation["step"],"agent" => agent_name,"lastOpponentAction" => -1)
    append!(history, [state])
    return history
  end
  # Cập nhật lại history với lượt phản hồi của đối thủ
  function update_history(history, observation)
    history[length(history)]["lastOpponentAction"] = observation["lastOpponentAction"]
    return history
  end

  function run(observation)
    # Mặc định sẽ chọn chiến thuật random
    agent_name = "random_agent"
    if(observation["step"] < 2)
      agent_name = "random_agent"
    else
      history = update_history(history,observation)
      last_observation = history[length(history)-1]
      gt_action = history[length(history)]["lastOpponentAction"]
      best_prob = -1
      best_agent = "random_agent"
      # Chạy các agent để xem chúng có thắng lượt vừa rồi không và tính xác suất thắng của chúng ở lượt này theo phân phối Beta
      for (name, agent) in agents
        pred_action = agent.run(last_observation)
        if(haskey(bandit_state, name))
          bandit_state[name][1] = (bandit_state[name][1] .- 1) ./ decay_rate .+ 1
          bandit_state[name][2] = (bandit_state[name][2] .- 1) ./ decay_rate .+ 1
        else
          bandit_state[name] = [1.0, 1.0]
        end
        if(pred_action == gt_action)
          bandit_state[name][1] = bandit_state[name][1] .+ reward/2
          bandit_state[name][2] = bandit_state[name][2] .+ reward/2
        else
          counter_agent_action = (pred_action + 1) % 3
          if(gt_action == counter_agent_action)
            bandit_state[name][2] = bandit_state[name][2] .+ reward
          else
            bandit_state[name][1] = bandit_state[name][1] .+ reward
          end
        end
        prob = rand(Beta(bandit_state[name][1], bandit_state[name][2]), 1)[1]
        if(prob > best_prob)
          best_agent = name
          best_prob = prob
        end
      end
      agent_name = best_agent
    end
    # Dùng agent đã lựa chọn để ra quyết định
    step = agents[agent_name].run(observation)
    # Lưu thông tin về lượt hiện tại
    history = store_current_state(history, observation, agent_name)
    return step
  end
  () -> (run, agent_count)
end
end


text:
@json struct Author
    raw::String
    user::User
end
@json_struct Author

@json struct Participant
    user::User
    role::String
    approved::Bool
    state::String
    participated_on::Date
end
@json_struct Participant

@json struct Commit
    hash::String
    data::DateTime
    author::Author
    message::String
    summary::Dict{String, Any}
    parents::Vector{Commit}
    #not present in parent commit objects
    repository::Repo
    participants::Vector{Participant}
end
@json_struct Commit

endpoint(::BitbucketAPI, ::typeof(get_commit), workspace::AStr, repo::AStr, ref::AStr) =
    Endpoint(:GET, "/repositories/$workspace/$repo/commit/$ref")
@not_implemented(::BitbucketAPI, ::typeof(get_commit), ::UUID, ::String)
into(::BitbucketAPI, ::typeof(get_commit)) = Commit


text:
using MAT
using Printf

include("../core/layers.jl")
include("../core/blocks.jl")
include("../../configs.jl")

struct _InitLayer conv_bn; end

function _InitLayer(;input_dim=3, output_dim=64, dtype=Array{Float32}, pdrop=0, bias=true)
    return _InitLayer(
        ConvBn(7, 7, input_dim, output_dim, bias=bias, padding=3, stride=2, f=relu, dtype=dtype, pdrop=pdrop)
    )
end

function (il::_InitLayer)(x; train=true)
    return pool(il.conv_bn(x, train=train), window=3, stride=2, padding=1)
end

"""
Required ResNet structures are implemented there.
"""

mutable struct ResNet50 layer1; layer2; layer3; layer4; layer5; fc; dtype; classes; end

function ResNet50(;input_dim=3, dtype=Array{Float32}, pdrop=0, bias=false, include_top=true)
    layer1 = _InitLayer(input_dim=input_dim, dtype=dtype, pdrop=pdrop)
    layer2 = Chain([
        Residual_1x3x1(64,  [64, 64, 256], dtype=dtype, pdrop=pdrop, downsample=true),
        Residual_1x3x1(256, [64, 64, 256], dtype=dtype, pdrop=pdrop),
        Residual_1x3x1(256, [64, 64, 256], dtype=dtype, pdrop=pdrop)
    ])
    layer3 = Chain([
        Residual_1x3x1(256, [128, 128, 512], dtype=dtype, pdrop=pdrop, downsample=true, ds_3x3_stride=2),
        Residual_1x3x1(512, [128, 128, 512], dtype=dtype, pdrop=pdrop),
        Residual_1x3x1(512, [128, 128, 512], dtype=dtype, pdrop=pdrop),
        Residual_1x3x1(512, [128, 128, 512], dtype=dtype, pdrop=pdrop)
    ])
    layer4 = Chain([
        Residual_1x3x1(512,  [256, 256, 1024], dtype=dtype, pdrop=pdrop, downsample=true, ds_3x3_stride=2),
        Residual_1x3x1(1024, [256, 256, 1024], dtype=dtype, pdrop=pdrop),
        Residual_1x3x1(1024, [256, 256, 1024], dtype=dtype, pdrop=pdrop),
        Residual_1x3x1(1024, [256, 256, 1024], dtype=dtype, pdrop=pdrop),
        Residual_1x3x1(1024, [256, 256, 1024], dtype=dtype, pdrop=pdrop),
        Residual_1x3x1(1024, [256, 256, 1024], dtype=dtype, pdrop=pdrop)
    ])
    layer5 = Chain([
        Residual_1x3x1(1024, [512, 512, 2048], dtype=dtype, pdrop=pdrop, downsample=true, ds_3x3_stride=2),
        Residual_1x3x1(2048, [512, 512, 2048], dtype=dtype, pdrop=pdrop),
        Residual_1x3x1(2048, [512, 512, 2048], dtype=dtype, pdrop=pdrop),
    ])

    fc = nothing
    if include_top fc = Dense(2048, 1000) end

    return ResNet50(layer1, layer2, layer3, layer4, layer5, fc, dtype, nothing)
end

function (rn::ResNet50)(x; train=true, return_intermediate=true)
    c1 = rn.layer1(x, train=train)
    c2 = rn.layer2(c1, train=train)
    c3 = rn.layer3(c2, train=train)
    c4 = rn.layer4(c3, train=train)
    c5 = rn.layer5(c4, train=train)
    if rn.fc === nothing
        if return_intermediate return c2, c3, c4, c5, nothing
        else return c5
        end
    else
        h, w, c, n = size(c5)
        # assuming image is square shaped
        p = mat(pool(c5, window=h, mode=2))
        res = rn.fc(p, train=train)
        if return_intermediate return c2, c3, c4, c5, res
        else return res
        end
    end
end

function evaluate_model(rn::ResNet50, x, top=5)
    res = rn(convert(rn.dtype, x), train=false, return_intermediate=false)
    z1 = vec(Array(res))
    s1 = sortperm(z1,rev=true)
    p1 = exp.(logp(z1))
    for ind in s1[1:top]
        print("$(rn.classes[ind]): $(@sprintf("%.2f",p1[ind]*100))%\n")
    end
end

function load_mat_weights(rn::ResNet50, path; pre_weights=nothing)

    conv_w=nothing; conv_b=nothing; fc_w=nothing; fc_b=nothing;
    bn_mom=nothing; bn_b=nothing; bn_mult=nothing; rn.classes=nothing;

    if pre_weights === nothing
        conv_w, conv_b, fc_w, fc_b, bn_mom, bn_b, bn_mult, rn.classes = _get_params(path, dtype=rn.dtype)
    else
        conv_w, conv_b, fc_w, fc_b, bn_mom, bn_b, bn_mult = pre_weights
    end

    # setting fc layer parameters
    if rn.fc !== nothing && fc_w !== nothing
        rn.fc.w = convert(rn.dtype, fc_w)
        rn.fc.b = convert(rn.dtype, fc_b)
    end

    idx = 1
    mom_idx = 1

    # setting layer 1
    rn.layer1.conv_bn.conv.w = Param(conv_w[1])
    rn.layer1.conv_bn.conv.b = Param(conv_b)
    rn.layer1.conv_bn.bn.bn_moments = bnmoments(mean=bn_mom[1], var=bn_mom[2])
    rn.layer1.conv_bn.bn.bn_params = Param(vcat(vec(bn_mult[1]), vec(bn_b[1])))
    idx += 1; mom_idx += 2;

    # setting layer 2
    rn.layer2.layers[1] = load_mat_weights(
        rn.layer2.layers[1], conv_w[idx:idx+3], bn_mom[mom_idx:mom_idx+7], bn_b[idx:idx+3], bn_mult[idx:idx+3]
    )
    idx += 4; mom_idx += 8;

    for i in 2:size(rn.layer2.layers, 1)
        rn.layer2.layers[i] = load_mat_weights(
            rn.layer2.layers[i], conv_w[idx:idx+2], bn_mom[mom_idx:mom_idx+5], bn_b[idx:idx+2], bn_mult[idx:idx+2]
        )
        idx += 3; mom_idx += 6;
    end

    # setting layer 3
    rn.layer3.layers[1] = load_mat_weights(
        rn.layer3.layers[1], conv_w[idx:idx+3], bn_mom[mom_idx:mom_idx+7], bn_b[idx:idx+3], bn_mult[idx:idx+3]
    )
    idx += 4; mom_idx += 8;

    for i in 2:size(rn.layer3.layers, 1)
        rn.layer3.layers[i] = load_mat_weights(
            rn.layer3.layers[i], conv_w[idx:idx+2], bn_mom[mom_idx:mom_idx+5], bn_b[idx:idx+2], bn_mult[idx:idx+2]
        )
        idx += 3; mom_idx += 6;
    end

    # setting layer 4
    rn.layer4.layers[1] = load_mat_weights(
        rn.layer4.layers[1], conv_w[idx:idx+3], bn_mom[mom_idx:mom_idx+7], bn_b[idx:idx+3], bn_mult[idx:idx+3]
    )
    idx += 4; mom_idx += 8;

    for i in 2:size(rn.layer4.layers, 1)
        rn.layer4.layers[i] = load_mat_weights(
            rn.layer4.layers[i], conv_w[idx:idx+2], bn_mom[mom_idx:mom_idx+5], bn_b[idx:idx+2], bn_mult[idx:idx+2]
        )
        idx += 3; mom_idx += 6;
    end

    # setting layer 5
    rn.layer5.layers[1] = load_mat_weights(
        rn.layer5.layers[1], conv_w[idx:idx+3], bn_mom[mom_idx:mom_idx+7], bn_b[idx:idx+3], bn_mult[idx:idx+3]
    )
    idx += 4; mom_idx += 8;

    for i in 2:size(rn.layer5.layers, 1)
        rn.layer5.layers[i] = load_mat_weights(
            rn.layer5.layers[i], conv_w[idx:idx+2], bn_mom[mom_idx:mom_idx+5], bn_b[idx:idx+2], bn_mult[idx:idx+2]
        )
        idx += 3; mom_idx += 6;
    end
    return rn
end

function _get_params(path; dtype=Array{Float32})
    wr50 = matread(path)
    params = wr50["params"]
    total_params = size(params["value"], 2)

    classes_imagenet = wr50["meta"]["classes"]["description"]
    avg_img = wr50["meta"]["normalization"]["averageImage"]

    conv_w = []; conv_b = [];
    fc_w = nothing; fc_b = nothing;
    bn_mom = []; bn_b = []; bn_mult = [];

    for k in 1:total_params
        name = params["name"][k]
        value = convert(dtype, params["value"][k])
        if startswith(name, "bn") && endswith(name, "moments")
            push!(bn_mom, reshape(value[:,1], (1,1,size(value,1),1)))
            push!(bn_mom, (reshape(value[:,2], (1,1,size(value,1),1)).^2) .- 1e-5)
        elseif startswith(name, "bn") && endswith(name, "bias")
            push!(bn_b, reshape(value, (1,1,length(value),1)))
        elseif startswith(name, "bn")
            push!(bn_mult, reshape(value, (1,1,length(value),1)))
        elseif startswith(name, "fc") && endswith(name, "filter")
            fc_w = reshape(value,(size(value,3),size(value,4)))'
        elseif startswith(name, "fc") && endswith(name, "bias")
            fc_b = value
        elseif endswith(name, "filter")
            push!(conv_w, value)
        elseif endswith(name, "bias")
            conv_b = reshape(value, (1,1,length(value),1))
        end
    end
    return conv_w, conv_b, fc_w, fc_b, bn_mom, bn_b, bn_mult, classes_imagenet
end

max_stars_repo_name: barisbatuhan/retinaface

max_issues_repo_name: barisbatuhan/retinaface

max_forks_repo_name: barisbatuhan/retinaface


text:
export Model

"""
    Model

Data structure to store some parameters about the simulation

- rho : Background density
- g : Gravity
- buoyancy_freq_n : Background stratification
- odg_b : Amplitude of the buoyancy
"""
mutable struct Model

    rho :: Float64
    g :: Float64
    buoyancy_freq_n :: Float64
    odg_b :: Float64
    hv_order :: Int
    hv_val :: Float64

    function Model( f0 )

    rho=1e3
    g = 9.81
    buoyancy_freq_n = 3 * f0
    odg_b = 1e-3
    hv_order = 8
    hv_val = 0

    new( rho, g, buoyancy_freq_n, odg_b, hv_order, hv_val)

    end

end


text:
"""
Computes the max-flow/min-cut between source and target using
Boykov-Kolmogorov algorithm.

Returns the maximum flow in the network, the flow matrix and
the partition {S,T} in the form of a vector of 1's and 2's.
The partition vector may also contain 0's. These can be
assigned any label (1 or 2), it is a user choice.

For further details, please refer to the paper:

BOYKOV, Y.; KOLMOGOROV, V., 2004. An Experimental Comparison of
Min-Cut/Max-Flow Algorithms for Energy Minimization in Vision.

Uses a default capacity of 1 when the capacity matrix isn\'t specified.

Requires arguments:
residual_graph::ADiGraph                # the input graph
source::Int                            # the source vertex
target::Int                            # the target vertex
capacity_matrix::AbstractMatrix{T}    # edge flow capacities

Author: Júlio Hoffimann Mendes (juliohm@stanford.edu)
"""
function boykov_kolmogorov_impl(
    residual_graph::ADiGraph,               # the input graph
    source,                           # the source vertex
    target,                           # the target vertex
    capacity_matrix::AbstractMatrix{T}    # edge flow capacities
    ) where T<:Number

    n = nv(residual_graph)

    flow = 0
    flow_matrix = spzeros(T, n, n)

    TREE = zeros(Int, n)
    TREE[source] = 1
    TREE[target] = 2

    PARENT = zeros(Int, n)

    A = [source,target]
    O = Vector{Int}()

    while true
        # growth stage
        path = find_path!(residual_graph, source, target, flow_matrix, capacity_matrix, PARENT, TREE, A)

        isempty(path) && break

        # augmentation stage
        flow += augment!(path, flow_matrix, capacity_matrix, PARENT, TREE, O)

        # adoption stage
        adopt!(residual_graph, source, target, flow_matrix, capacity_matrix, PARENT, TREE, A, O)
    end

    return flow, flow_matrix, TREE
end

function find_path!(
    residual_graph::ADiGraph,               # the input graph
    source,                           # the source vertex
    target,                           # the target vertex
    flow_matrix::AbstractMatrix{T},       # the current flow matrix
    capacity_matrix::AbstractMatrix{T},   # edge flow capacities
    PARENT::Vector{Int},                   # parent table
    TREE::Vector{Int},                     # tree table
    A::Vector{Int}                         # active set
    ) where T<:Number

    tree_cap(p,q) = TREE[p] == 1 ? capacity_matrix[p,q] - flow_matrix[p,q] :
                                    capacity_matrix[q,p] - flow_matrix[q,p]
    while !isempty(A)
        p = last(A)
        for q in neighbors(residual_graph, p)
            if tree_cap(p,q) > 0
                if TREE[q] == 0
                    TREE[q] = TREE[p]
                    PARENT[q] = p
                    pushfirst!(A, q)
                end
                if TREE[q] ≠ 0 && TREE[q] ≠ TREE[p]
                    # p -> source
                    path_to_source = [p]
                    while PARENT[p] ≠ 0
                        p = PARENT[p]
                        push!(path_to_source, p)
                    end

                    # q -> target
                    path_to_target = [q]
                    while PARENT[q] ≠ 0
                        q = PARENT[q]
                        push!(path_to_target, q)
                    end

                    # source -> target
                    path = [reverse!(path_to_source); path_to_target]

                    if path[1] == source && path[end] == target
                        return path
                    elseif path[1] == target && path[end] == source
                        return reverse!(path)
                    end
                end
            end
        end
        pop!(A)
    end

    return Vector{Int}()
end

function augment!(
    path::AbstractVector,                  # path from source to target
    flow_matrix::AbstractMatrix{T},       # the current flow matrix
    capacity_matrix::AbstractMatrix{T},   # edge flow capacities
    PARENT::Vector{Int},                   # parent table
    TREE::Vector{Int},                     # tree table
    O::Vector{Int}                         # orphan set
    ) where T<:Number

    # bottleneck capacity
    Δ = Inf
    for i=1:length(path)-1
        p, q = path[i:i+1]
        cap = capacity_matrix[p,q] - flow_matrix[p,q]
        cap < Δ && (Δ = cap)
    end

    # update residual graph
    for i=1:length(path)-1
        p, q = path[i:i+1]
        flow_matrix[p,q] += Δ
        flow_matrix[q,p] -= Δ

        # create orphans
        if flow_matrix[p,q] == capacity_matrix[p,q]
            if TREE[p] == TREE[q] == 1
                PARENT[q] = 0
                pushfirst!(O, q)
            end
            if TREE[p] == TREE[q] == 2
                PARENT[p] = 0
                pushfirst!(O, p)
            end
        end
    end

    return Δ
end

function adopt!(
    residual_graph::ADiGraph,               # the input graph
    source::Int,                           # the source vertex
    target::Int,                           # the target vertex
    flow_matrix::AbstractMatrix{T},       # the current flow matrix
    capacity_matrix::AbstractMatrix{T},   # edge flow capacities
    PARENT::Vector{Int},                   # parent table
    TREE::Vector{Int},                     # tree table
    A::Vector{Int},                        # active set
    O::Vector{Int}                         # orphan set
    ) where T<:Number

    tree_cap(p,q) = TREE[p] == 1 ? capacity_matrix[p,q] - flow_matrix[p,q] :
    capacity_matrix[q,p] - flow_matrix[q,p]
    while !isempty(O)
        p = pop!(O)
        # try to find parent that is not an orphan
        parent_found = false
        for q in neighbors(residual_graph, p)
            if TREE[q] == TREE[p] && tree_cap(q,p) > 0
                # check if "origin" is either source or target
                o = q
                while PARENT[o] ≠ 0
                    o = PARENT[o]
                end
                if o == source || o == target
                    parent_found = true
                    PARENT[p] = q
                    break
                end
            end
        end

        if !parent_found
            # scan all neighbors and make the orphan a free node
            for q in neighbors(residual_graph, p)
                if TREE[q] == TREE[p]
                    if tree_cap(q,p) > 0
                        pushfirst!(A, q)
                    end
                    if PARENT[q] == p
                        PARENT[q] = 0
                        pushfirst!(O, q)
                    end
                end
            end

            TREE[p] = 0
            B = setdiff(A, p)
            resize!(A, length(B))[:] = B
        end
    end
end


text:
"""
    rectangulardomain(zbl, ztr, Δr)

Generate initial mesh node coordinates for a rectangular domain extending from the complex
bottom left corner `zbl` to the top right corner `ztr` with initial mesh step `Δr`.
"""
function rectangulardomain(zbl::Complex, ztr::Complex, Δr)
    rzbl, izbl = reim(zbl)
    rztr, iztr = reim(ztr)

    X = rztr - rzbl
    Y = iztr - izbl

    n = ceil(Int, Y/Δr)
    dy = Y/n

    ## dx = sqrt(Δr² - (dy/2)²), solved for equilateral triangle
    m = ceil(Int, X/sqrt(Δr^2 - dy^2/4))
    dx = X/m
    half_dx = dx/2

    T = promote_type(ComplexF64, typeof(zbl), typeof(ztr), typeof(Δr))
    v = Vector{T}()
    sizehint!(v, (m+1)*(n+1))

    shift = false  # we will displace every other line by dx/2
    for j = 0:n
        y = izbl + dy*j

        for i = 0:m
            x = rzbl + dx*i

            if shift && i == 0
                continue  # otherwise, we shift out of left bound
            elseif shift
                x -= half_dx
            end

            if i == m
                shift = !shift
            end

            push!(v, complex(x, y))
        end
    end

    return v
end

"""
    diskdomain(R, Δr)

Generate initial mesh coordinates for a circular disk domain of radius `R` and
center (0, 0) for ``|z| < R`` and initial mesh step `Δr`.
"""
function diskdomain(R, Δr)
    h = Δr*sqrt(3)/2
    n = 1 + round(Int, R/h)
    R_n = R/n

    f₀ = 0.0
    np = 6

    f₀step = π/(6n)

    T = promote_type(typeof(R), typeof(Δr), Float64)
    newnodes = Vector{complex(T)}(undef, 6*sum(1:n)+1)
    newnodes[1] = 0

    idx = 2
    for ii = 1:n
        # precalculate
        iiR_n = ii*R_n
        step = 2π/np

        for jj = 0:np-1
            f = f₀ + step*jj

            newnodes[idx] = iiR_n*cis(f)
            idx += 1
        end

        f₀ += f₀step
        np += 6
    end
    return newnodes
end


text:
using ClangCompiler

const CC = ClangCompiler

@testset "Status" begin
    src = joinpath(@__DIR__, "..", "code", "main.cpp")
    args = get_compiler_args()
    push!(args, "-std=c++11")
    pushfirst!(args, "clang")
    push!(args, "-I$(@__DIR__)")

    instance = CC.CompilerInstance()
    CC.create_diagnostics(instance)
    diag = CC.get_diagnostics(instance)

    invok = CC.create_compiler_invocation_from_cmd(src, args, diag)
    CC.set_invocation(instance, invok)
    CC.set_target(instance)

    CC.PrintStats(instance, CC.HeaderSearchOptions)
    CC.PrintStats(instance, CC.DiagnosticOptions)
    CC.PrintStats(instance, CC.FrontendOptions)
    CC.PrintStats(instance, CC.CodeGenOptions)
    CC.PrintStats(instance, CC.PreprocessorOptions)
    CC.PrintStats(instance, CC.TargetOptions)

    CC.create_file_manager(instance)
    CC.create_source_manager(instance)

    CC.get_file(instance, src)

    CC.create_preprocessor(instance)

    CC.PrintStats(instance, CC.HeaderSearch)

    CC.dispose(instance)
end


text:
"""
Module for working with [Keyhole Markup Language (KML) file format](https://en.wikipedia.org/wiki/Keyhole_Markup_Language).

Be aware that this module only support a part of gx:Track extension only.

Have a look at https://developers.google.com/kml/documentation for more information.
"""
module KMLTracks

using Dates
using TimeZones
using EzXML

export read_kml_file, parse_kml_string

"""
Constant defining KML Format version
"""
const KML_VERSION = "2.2"

"""
Constant dictionary defining default namespaces
"""
const KML_NS = Dict(
    "" => "http://www.opengis.net/kml/2.2",
    "gx" => "http://www.google.com/kml/ext/2.2"
)

const dt0 = ZonedDateTime(0, tz"UTC")

"""
    KMLTrackPoint(time::ZonedDateTime, longitude::Float64, latitude::Float64, altitude::Float64)

Point in a KML Track.
"""
struct KMLTrackPoint
    time::ZonedDateTime
    longitude::Float64
    latitude::Float64
    altitude::Float64
end

"""
    KMLTrack()

A KMLTrack is a collection of track points.

`gx:Track` element is an extension of the OGC KML 2.2 standard and is supported in Google Earth 5.2 and later.
"""
struct KMLTrack
    points::Vector{KMLTrackPoint}
    KMLTrack() = new(Vector{KMLTrackPoint}[])
end

"""
    KMLPlacemark()

Placemark is a special kind of Feature. It contains all of the elements that belong to Feature, and it adds some elements that are specific to the Placemark element.

See `Placemark` tag.
"""
struct KMLPlacemark
    track::KMLTrack
    KMLPlacemark() = new(KMLTrack())
end

"""
    KMLDocument()

A KML Document is a container for placemark.

Be aware that `KMLDocument` is related to `kml` tag not to `Document` tag.
"""
struct KMLDocument
    placemark::KMLPlacemark
    KMLDocument() = new(KMLPlacemark())
end

"""
    read_kml_file(fname) -> KMLDocument

Read KML file from filename `fname`.
"""
function read_kml_file(fname)
    xdoc = readxml(fname)
    return _parse_kml(xdoc)
end

"""
    parse_kml_string(s) -> KMLDocument

Parse GPX data from String `s`.
"""
function parse_kml_string(s)
    xdoc = parsexml(s)
    return _parse_kml(xdoc)
end

"""
    _parse_kml(xdoc::EzXML.Document) -> KMLDocument

Parse `EzXML.Document` and return a `KMLDocument`.
"""
function _parse_kml(xdoc::EzXML.Document)
    kmls = root(xdoc)

    kmldoc = KMLDocument()

    ns = namespaces(xdoc.root)
    d_ns = Dict(ns)
    @assert(d_ns == KML_NS, "Namespace error $d_ns different from $KML_NS")

    ns[1] = "x" => ns[1][2]  # change default namespace key

    a_when = findall("//x:Placemark/gx:Track/x:when/text()", xdoc.root, ns)
    a_when = nodecontent.(a_when)
    a_when = ZonedDateTime.(a_when, dateformat"yyyy-mm-ddTHH:MM:SSzzz")

    a_coords = findall("//x:Placemark/gx:Track/gx:coord/text()", xdoc.root, ns)
    a_coords = nodecontent.(a_coords)

    @assert(length(a_coords) == length(a_when), "a_coords/a_when length mismatch")

    for (i, dt) in enumerate(a_when)
        coord = a_coords[i]
        coord = replace(coord, "," => ".")  # fix issue with some badbly formatted kml files
        long, lat, alt = split(coord, " ")
        long, lat, alt = parse.(Float64, (long, lat, alt))
        p = KMLTrackPoint(dt, long, lat, alt)
        push!(kmldoc.placemark.track.points, p)
    end

    return kmldoc
end

end # module


text:
## test_karhunen_loeve.jl : test Karhunen-Lo\`eve expansion for GRF generation

@testset "Karhunen-Lo\u00e8ve expansion " begin

    ## 1d Mat\'ern ##
    cov = CovarianceFunction(1,Matern(0.3,1))
    pts = 0:0.01:1
    grf = GaussianRandomField(cov,KarhunenLoeve(500),pts)
    @test isa(grf,GaussianRandomField)
    @test isa(grf.cov,CovarianceFunction)
    @test isa(grf.cov.cov,Matern)
    @test ndims(grf.cov) == 1
    @test isa(grf,GaussianRandomField{KarhunenLoeve{500}})
    @test length(grf.pts[1]) == length(pts)
    @test length(sample(grf)) == length(pts)

    ## test anisotropic
    cov = CovarianceFunction(2,AnisotropicExponential(100*[1 0.2; 0.2 1]))
    pts1 = range(0,stop = 1,length = 128)
    pts2 = range(0,stop = 1,length = 128)
    grf = GaussianRandomField(cov,KarhunenLoeve(1000),pts1,pts2)
    @test isa(grf,GaussianRandomField)
    @test isa(grf.cov,CovarianceFunction)
    @test isa(grf.cov.cov,AnisotropicExponential)
    @test ndims(grf.cov) == 2
    @test isa(grf,GaussianRandomField{KarhunenLoeve{1000}})
    @test length(grf.pts) == 2
    @test length(grf.pts[1]) == length(pts1)
    @test length(grf.pts[2]) == length(pts2)
    @test size(sample(grf),1) == length(pts1)
    @test size(sample(grf),2) == length(pts2)

    ## variation of number of KL terms ##
    @test_throws DomainError KarhunenLoeve{0}()
    @test_throws DomainError KarhunenLoeve(0)

    cov = CovarianceFunction(2,Matern(0.3,1))
    pts = 0:0.01:1
    grf = GaussianRandomField(cov,KarhunenLoeve(17),pts,pts, eigensolver=EigenSolver())
    @test isa(grf,GaussianRandomField)
    @test isa(grf.cov,CovarianceFunction)
    @test isa(grf.cov.cov,Matern)
    @test ndims(grf.cov) == 2
    @test isa(grf,GaussianRandomField{KarhunenLoeve{17}})
    @test length(grf.pts[1]) == length(pts)
    @test length(grf.pts[2]) == length(pts)
    @test size(sample(grf),1) == length(pts)
    @test size(sample(grf),2) == length(pts)

    cov = CovarianceFunction(2,Matern(0.3,1))
    pts = 0:0.01:1
    grf = GaussianRandomField(cov,KarhunenLoeve(1),pts,pts, eigensolver=EigenSolver())
    @test isa(grf,GaussianRandomField)
    @test isa(grf.cov,CovarianceFunction)
    @test isa(grf.cov.cov,Matern)
    @test ndims(grf.cov) == 2
    @test isa(grf,GaussianRandomField{KarhunenLoeve{1}})
    @test length(grf.pts[1]) == length(pts)
    @test length(grf.pts[2]) == length(pts)
    @test size(sample(grf),1) == length(pts)
    @test size(sample(grf),2) == length(pts)

    mat = Matern(0.1,1.0)
    cov = CovarianceFunction(1,mat)
    pts1 = 0:0.01:1
    nterms = [1 2 5 10 20 50 100 200 500 1000]
    err = Float64[]
    for n in nterms
        grf = GaussianRandomField(cov,KarhunenLoeve(n),pts1)
        push!(err,rel_error(grf))
    end
    @test err[1]/err[end] > 1e7

    ## quadrature methods ##
    mat = Matern(0.1,1.0)
    cov = CovarianceFunction(1,mat)
    pts1 = 0:0.025:1
    grf = GaussianRandomField(cov,KarhunenLoeve(1000),pts1,quad=GaussLegendre())
    @test rel_error(grf) < 1e-6
    grf = GaussianRandomField(cov,KarhunenLoeve(1000),pts1,quad=EOLE())
    @test rel_error(grf) < 1e-6
    @suppress grf = GaussianRandomField(cov,KarhunenLoeve(1000),pts1,quad=Simpson())
    @test rel_error(grf) < 1e-6
    grf = GaussianRandomField(cov,KarhunenLoeve(1000),pts1,quad=Midpoint())
    @test rel_error(grf) < 1e-6
    grf = GaussianRandomField(cov,KarhunenLoeve(1000),pts1,quad=Trapezoidal())
    @test rel_error(grf) < 1e-6

    ## test 3D Gaussian ##
    cov = CovarianceFunction(3,Gaussian(1.0))
    pts1 = 0:0.05:1
    pts2 = 0:0.05:1
    pts3 = 0:0.05:1
    grf = GaussianRandomField(cov,KarhunenLoeve(103),pts1,pts2,pts3)
    @test isa(grf,GaussianRandomField)
    @test isa(grf.cov,CovarianceFunction)
    @test isa(grf.cov.cov,Gaussian)
    @test ndims(grf.cov) == 3
    @test isa(grf,GaussianRandomField{KarhunenLoeve{103}})
    @test length(grf.pts) == 3
    @test length(grf.pts[1]) == length(pts1)
    @test length(grf.pts[2]) == length(pts2)
    @test length(grf.pts[3]) == length(pts3)
    @test size(sample(grf),1) == length(pts1)
    @test size(sample(grf),2) == length(pts2)
    @test size(sample(grf),3) == length(pts3)

    ## 1d Whittle ##
    cov = CovarianceFunction(1,Whittle(1))
    pts = 0:0.01:1
    grf = GaussianRandomField(cov,KarhunenLoeve(10),pts)
    @test isa(grf,GaussianRandomField)
    @test isa(grf.cov,CovarianceFunction)
    @test isa(grf.cov.cov,Whittle)
    @test ndims(grf.cov) == 1
    @test isa(grf,GaussianRandomField{KarhunenLoeve{10}})
    @test length(grf.pts[1]) == length(pts)
    @test length(sample(grf)) == length(pts)

    ## 1d Linear ##
    cov = CovarianceFunction(1,Linear(3//2))
    pts = 0:0.01:1
    grf = GaussianRandomField(cov,KarhunenLoeve(25),pts)
    @test isa(grf,GaussianRandomField)
    @test isa(grf.cov,CovarianceFunction)
    @test isa(grf.cov.cov,Linear)
    @test ndims(grf.cov) == 1
    @test isa(grf,GaussianRandomField{KarhunenLoeve{25}})
    @test length(grf.pts[1]) == length(pts)
    @test length(sample(grf)) == length(pts)

    ## 1d Spherical ##
    cov = CovarianceFunction(1,Spherical(1.5))
    pts = 0:0.02:1
    grf = GaussianRandomField(cov,KarhunenLoeve(250),pts)
    @test isa(grf,GaussianRandomField)
    @test isa(grf.cov,CovarianceFunction)
    @test isa(grf.cov.cov,Spherical)
    @test ndims(grf.cov) == 1
    @test isa(grf,GaussianRandomField{KarhunenLoeve{250}})
    @test length(grf.pts[1]) == length(pts)
    @test length(sample(grf)) == length(pts)

    ## variation of number of points ##
    cov = CovarianceFunction(1,Matern(1,1))
    pts = 0:0.1:0
    @test_throws ArgumentError grf = GaussianRandomField(cov,KarhunenLoeve(300),pts)
    pts = 0:0.1:0.1
    grf = GaussianRandomField(cov,KarhunenLoeve(300),pts)
    pts = 1:-0.1:0
    grf = GaussianRandomField(cov,KarhunenLoeve(300),pts)
    grf = GaussianRandomField(cov,KarhunenLoeve(300),pts,nq=1000)
    @test_throws ArgumentError GaussianRandomField(cov,KarhunenLoeve(300),pts,nq=100)

    ## non-SPD covariance matrix ##
    cov = CovarianceFunction(2,SquaredExponential(0.5))
    @suppress grf = GaussianRandomField(cov,KarhunenLoeve(200),pts1,pts2)
    @test isa(grf,GaussianRandomField)
    @test isa(grf.cov,CovarianceFunction)
    @test isa(grf.cov.cov,SquaredExponential)
    @test ndims(grf.cov) == 2
    @test isa(grf,GaussianRandomField{<:KarhunenLoeve})
    @test length(grf.pts[1]) == length(pts1)
    @test length(grf.pts[2]) == length(pts2)
    @test size(sample(grf),1) == length(pts1)
    @test size(sample(grf),2) == length(pts2)

    ## test domain with negative sides ##
    cov = CovarianceFunction(2,Matern(1.,2.5,σ=1.,p=2))
    pts1 = collect(-2:0.1:0)
    pts2 = collect(-1:0.1:1)
    grf = GaussianRandomField(cov,KarhunenLoeve(458),pts1,pts2)
    @test isa(grf,GaussianRandomField)
    @test isa(grf.cov,CovarianceFunction)
    @test isa(grf.cov.cov,Matern)
    @test ndims(grf.cov) == 2
    @test isa(grf,GaussianRandomField{KarhunenLoeve{458}})
    @test length(grf.pts[1]) == length(pts1)
    @test length(grf.pts[2]) == length(pts2)
    @test size(sample(grf),1) == length(pts1)
    @test size(sample(grf),2) == length(pts2)
end


text:
using Test
using TreeKnit.SRG
using TreeTools

## Simple case
nwk1 = "(A:3,(B:2,C:2)BC:1)R"
nwk2 = "((A:3,B:3)AB:1,C:4)R"

t1 = parse_newick_string(nwk1);
t2 = parse_newick_string(nwk2);

MCCs = [["A","B"], ["C"]]

arg, rlm, lm1, lm2 = SRG.arg_from_trees(t1, t2, MCCs)

## Case with a lot of singletons
nwk1 = "(((((A:2,X3)AX3:1,X4)AX4:2,X5)AX5:2,B)AB:0.,X1,X2):0.R"
nwk2 = "((((A:1,X2)AX2:3,X1)AX1:3,B)AB:0.,X3,X4,X5):0.R"

t1 = parse_newick_string(nwk1);
t2 = parse_newick_string(nwk2);

MCCs = [["A","B"],["X1"],["X2"],["X3"],["X4"],["X5"]]

arg, rlm, lm1, lm2 = SRG.arg_from_trees(t1, t2, MCCs)

## Case with two different roots
nwk1 = "(((A:1,B:1)AB:1,C:2)ABC:1,D:3)R"
nwk2 = "((A:3,D:3)AD:1,(B:3,C:3)BC:1)R"

t1 = parse_newick_string(nwk1);
t2 = parse_newick_string(nwk2);

MCCs = [["A","D"], ["B"], ["C"]]

arg, rlm, lm1, lm2 = SRG.arg_from_trees(t1, t2, MCCs)

max_stars_repo_name: PierreBarrat/RecombTools

max_issues_repo_name: PierreBarrat/RecombTools

max_forks_repo_name: PierreBarrat/RecombTools


text:
# Copyright (c) 2019 Harshavardhan P K
#
# This software is released under the MIT License.
# https://opensource.org/licenses/MIT

function transform!(arr::Array{Int64,1}, noun::Int64 = 12, verb::Int64 = 2)
    start::Int16 = 1
    arr[2] = noun
    arr[3] = verb
    while arr[start] ≠ 99
        if arr[start] == 1
            arr[arr[start+3]+1] = arr[arr[start+1]+1] + arr[arr[start+2]+1]
        elseif arr[start] == 2
            arr[arr[start+3]+1] = arr[arr[start+1]+1] * arr[arr[start+2]+1]
        end
        start += 4
    end
end

function getData(fname::String)::Array{Int64,1}
    f = open(fname)
    line = read(f, String)
    line = split(strip(line), ',')
    return map(x -> parse(Int64, x), line)
end

inp = getData("input")
transform!(inp)
println(inp[1])

inp = getData("input")
for n ∈ 1:100
    for v ∈ 1:100
        i = deepcopy(inp)
        transform!(i, n, v)
        if i[1] == 19690720
            println(100 * n + v)
        end
    end
end


text:
# Install PyCall
import Pkg
Pkg.add("PyCall")


text:
module BlasIntrinsics

using Base:llvmcall

# https://software.intel.com/sites/landingpage/IntrinsicsGuide/#techs=FMA

# Plain baseline registers
const half2   = NTuple{2, VecElement{Float16}}
const half4   = NTuple{4, VecElement{Float16}}
#
const float2  = NTuple{2, VecElement{Float32}}
const float4  = NTuple{4, VecElement{Float32}}
const float8  = NTuple{8, VecElement{Float32}}
const float16 = NTuple{16,VecElement{Float32}}
#
const double2 = NTuple{2, VecElement{Float64}}
const double4 = NTuple{4, VecElement{Float64}}
const double8 = NTuple{8, VecElement{Float64}}

@inline m256dfma(a,b,c) = ccall("llvm.fma.v4f64", llvmcall, double4, (double4, double4, double4), a, b, c)
@inline m256fma(a,b,c ) = ccall("llvm.fma.v8f32", llvmcall, float8, (float8, float8, float8), a, b, c)

@inline m512dfma(a,b,c) = ccall("llvm.fma.v8f64", llvmcall, double8, (double8, double8, double8), a, b, c)
@inline m512fma(a,b,c ) = ccall("llvm.fma.v16f32", llvmcall, float16, (float16, float16, float16), a, b, c)

@inline prefetch(a, i1, i2, i3) = ccall("llvm.prefetch", llvmcall, Cvoid, (Ptr{Cchar}, Int32, Int32, Int32), a, i1, i2, i3)

@inline prefetchT0(a) = mmprefetch(a, Int32(0), Int32(3), Int(1))
@inline prefetchT1(a) = mmprefetch(a, Int32(0), Int32(2), Int(1))
@inline prefetchT2(a) = mmprefetch(a, Int32(0), Int32(1), Int(1))

@inline function mm_mul(x::T, y::T) where T <:AbstractFloat
   return x * y
end

@inline function mm_mul(x::float4, y::float4)
  (VecElement(x[1].value * y[1].value),
   VecElement(x[2].value * y[2].value),
   VecElement(x[3].value * y[3].value),
   VecElement(x[4].value * y[4].value))
end

@inline function mm_mul(x::double2, y::double2)
  (VecElement(x[1].value * y[1].value),
   VecElement(x[2].value * y[2].value))
end

@inline function mm_mul(x::float8, y::float8)
  (VecElement(x[1].value * y[1].value),
   VecElement(x[2].value * y[2].value),
   VecElement(x[3].value * y[3].value),
   VecElement(x[4].value * y[4].value),
   VecElement(x[5].value * y[5].value),
   VecElement(x[6].value * y[6].value),
   VecElement(x[7].value * y[7].value),
   VecElement(x[8].value * y[8].value))
end

@inline function mm_mul(x::double4, y::double4)
  (VecElement(x[1].value * y[1].value),
   VecElement(x[2].value * y[2].value),
   VecElement(x[3].value * y[3].value),
   VecElement(x[4].value * y[4].value))
end

@inline function mm_add(x::T, y::T) where T <:AbstractFloat
   return x + y
end

@inline function mm_add(x::float4, y::float4)
  (VecElement(x[1].value + y[1].value),
   VecElement(x[2].value + y[2].value),
   VecElement(x[3].value + y[3].value),
   VecElement(x[4].value + y[4].value))
end

@inline function mm_add(x::double2, y::double2)
  (VecElement(x[1].value + y[1].value),
   VecElement(x[2].value + y[2].value))
end

@inline function mm_add(x::float8, y::float8)
  (VecElement(x[1].value + y[1].value),
   VecElement(x[2].value + y[2].value),
   VecElement(x[3].value + y[3].value),
   VecElement(x[4].value + y[4].value),
   VecElement(x[5].value + y[5].value),
   VecElement(x[6].value + y[6].value),
   VecElement(x[7].value + y[7].value),
   VecElement(x[8].value + y[8].value))
end

@inline function mm_add(x::double4, y::double4)
  (VecElement(x[1].value + y[1].value),
   VecElement(x[2].value + y[2].value),
   VecElement(x[3].value + y[3].value),
   VecElement(x[4].value + y[4].value))
end

@inline function mm_sub(x::T, y::T) where T <:AbstractFloat
   return x - y
end

@inline function mm_sub(x::float4, y::float4)
  (VecElement(x[1].value - y[1].value),
   VecElement(x[2].value - y[2].value),
   VecElement(x[3].value - y[3].value),
   VecElement(x[4].value - y[4].value))
end

@inline function mm_sub(x::double2, y::double2)
  (VecElement(x[1].value - y[1].value),
   VecElement(x[2].value - y[2].value))
end

@inline function mm_sub(x::float8, y::float8)
  (VecElement(x[1].value - y[1].value),
   VecElement(x[2].value - y[2].value),
   VecElement(x[3].value - y[3].value),
   VecElement(x[4].value - y[4].value),
   VecElement(x[5].value - y[5].value),
   VecElement(x[6].value - y[6].value),
   VecElement(x[7].value - y[7].value),
   VecElement(x[8].value - y[8].value))
end

@inline function mm_sub(x::double4, y::double4)
  (VecElement(x[1].value - y[1].value),
   VecElement(x[2].value - y[2].value),
   VecElement(x[3].value - y[3].value),
   VecElement(x[4].value - y[4].value))
end

@inline function mm_mad(a::T, x::T, y::T) where T <:AbstractFloat
   return a * x + y
end

@inline function mm_mad(a::T, x::float4, y::float4) where T <: AbstractFloat
  (VecElement(a * x[1].value + y[1].value),
   VecElement(a * x[2].value + y[2].value),
   VecElement(a * x[3].value + y[3].value),
   VecElement(a * x[4].value + y[4].value))
end

@inline function mm_mad(a::T, x::double2, y::double2) where T <: AbstractFloat
  (VecElement(a * x[1].value + y[1].value),
   VecElement(a * x[2].value + y[2].value))
end

@inline function mm_mad(a::T, x::float8, y::float8) where T <: AbstractFloat
   z = float8((a, a, a, a, a, a, a, a))
   return m256fma(z,x,y)
end

@inline function mm_mad(a::T, x::double4, y::double4) where T <: AbstractFloat
   z = double4((a, a, a, a))
   return m256dfma(z,x,y)
end

end #QJuliaIntrinsics

max_stars_repo_name: alexstrel/solver-tuner

max_issues_repo_name: alexstrel/solver-tuner

max_forks_repo_name: alexstrel/solver-tuner


text:
"""
Bridges `CP.Count` to reification.
"""
struct Count2ReificationBridge{T, S <: MOI.AbstractScalarSet} <: MOIBC.AbstractBridge
    vars::Vector{MOI.VariableIndex}
    vars_bin::Vector{MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}}
    cons::Vector{MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, CP.Reification{S}}}
    con_sum::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}
end

function MOIBC.bridge_constraint(
    ::Type{Count2ReificationBridge{T, S}},
    model,
    f::MOI.VectorOfVariables,
    s::CP.Count{S},
) where {T, S <: MOI.AbstractScalarSet}
    return MOIBC.bridge_constraint(
        Count2ReificationBridge{T, S},
        model,
        MOI.VectorAffineFunction{T}(f),
        s,
    )
end

function MOIBC.bridge_constraint(
    ::Type{Count2ReificationBridge{T, S}},
    model,
    f::MOI.VectorAffineFunction{T},
    s::CP.Count{S},
) where {T, S <: MOI.AbstractScalarSet}
    vars, vars_bin = MOI.add_constrained_variables(model, [MOI.ZeroOne() for _ in 1:s.dimension])

    f_scalars = MOIU.scalarize(f)
    cons = MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, CP.Reification{S}}[
        MOI.add_constraint(
            model,
            MOIU.vectorize(
                MOI.ScalarAffineFunction{T}[
                    one(T) * vars[i],
                    f_scalars[1 + i],
                ]
            ),
            CP.Reification(s.set)
        )
        for i in 1:s.dimension
    ]

    con_sum = MOI.add_constraint(
        model,
        sum(one(T) .* vars) - f_scalars[1],
        MOI.EqualTo(zero(T))
    )

    return Count2ReificationBridge(vars, vars_bin, cons, con_sum)
end

function MOI.supports_constraint(
    ::Type{Count2ReificationBridge{T, S}},
    ::Union{Type{MOI.VectorOfVariables}, Type{MOI.VectorAffineFunction{T}}},
    ::Type{CP.Count{S}},
) where {T, S <: MOI.AbstractScalarSet}
    return true
end

function MOIB.added_constrained_variable_types(::Type{Count2ReificationBridge{T, S}}) where {T, S <: MOI.AbstractScalarSet}
    return [(MOI.ZeroOne,)]
end

function MOIB.added_constraint_types(::Type{Count2ReificationBridge{T, S}}) where {T, S <: MOI.AbstractScalarSet}
    return [
        (MOI.VectorAffineFunction{T}, CP.Reification{S}),
        (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}),
    ]
end

function MOI.get(b::Count2ReificationBridge{T, S}, ::MOI.NumberOfVariables) where {T, S <: MOI.AbstractScalarSet}
    return length(b.vars)
end

function MOI.get(
    b::Count2ReificationBridge{T, S},
    ::MOI.NumberOfConstraints{
        MOI.VariableIndex, MOI.ZeroOne,
    },
) where {T, S <: MOI.AbstractScalarSet}
    return length(b.vars_bin)
end

function MOI.get(
    b::Count2ReificationBridge{T, S},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{T}, CP.Reification{S},
    },
) where {T, S <: MOI.AbstractScalarSet}
    return length(b.cons)
end

function MOI.get(
    ::Count2ReificationBridge{T, S},
    ::MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{T}, MOI.EqualTo{T},
    },
) where {T, S <: MOI.AbstractScalarSet}
    return 1
end

function MOI.get(
    b::Count2ReificationBridge{T, S},
    ::MOI.ListOfVariableIndices,
) where {T, S <: MOI.AbstractScalarSet}
    return copy(b.vars)
end

function MOI.get(
    b::Count2ReificationBridge{T, S},
    ::MOI.ListOfConstraintIndices{
        MOI.VariableIndex, MOI.ZeroOne,
    },
) where {T, S <: MOI.AbstractScalarSet}
    return copy(b.vars_bin)
end

function MOI.get(
    b::Count2ReificationBridge{T, S},
    ::MOI.ListOfConstraintIndices{
        MOI.VectorAffineFunction{T}, CP.Reification{S},
    },
) where {T, S <: MOI.AbstractScalarSet}
    return copy(b.cons)
end

function MOI.get(
    b::Count2ReificationBridge{T, S},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T}, MOI.EqualTo{T},
    },
) where {T, S <: MOI.AbstractScalarSet}
    return [b.con_sum]
end


text:
# This file was generated, do not modify it. # hide
using Gnuplot
x = -2π:0.001:2π
@gp(x, sin.(x), "w l t 'sin'", "set yrange [-1.1:1.1]", "set grid",
    x, cos.(x), "with linespoints ls 1 t 'cos' ",
    "set style line 1 lc rgb 'black' lt 1 lw 1 pt 6 pi -200 ps 1.5")
    save(term="pngcairo size 600,400", output="plt2_ex3.png") # hide


text:
function getSoftmaxProbs(q_vector, temp, bias_vector; exploration_rate=0.00001)
    # softmax policy n actions
    # adds a bias value to the q values
    biased_q = q_vector .+ bias_vector
    biased_q = biased_q .- maximum(biased_q) # Avoid exp(sth VERY large) = Inf!!!
    nom = exp.(biased_q/temp)
    # nom = checkInf.(nom)
    # if isinf(variable) #|| ratioterm >= 1e16
    denom = sum(nom)
    prob = nom ./ denom
    #avoid some extrema (ln(0)). It's reasonable to assume a very small selection chance
    # even for very big differences in q-values
    # We don't want the LL be dominated by a few extrema (e.g. explorative actions)
    # therefore, shrink each p from [0, 1] to [0.001, 0.999], then normalize to 1
    # println("exploration_rate=$exploration_rate")
    # println("prob=$prob")
    boundedVals = exploration_rate .+ (1.0 - 2*exploration_rate)*prob
    boundedVals = boundedVals ./ sum(boundedVals) # normalize to 1
    # println("boundedVals2: $boundedVals")
end


text:
function load_data(name)
    datasets = ["iris", "adult", "digits", "causal"]
    data_path = joinpath(dirname(pathof(DecisionTree)), "..", "test/data/")

    if name == "digits"
        f = open(joinpath(data_path, "digits.csv"))
        data = readlines(f)[2:end]
        data = [[parse(Float32, i)
            for i in split(row, ",")]
            for row in data]
        data = hcat(data...)
        Y = Int.(data[1, 1:end]) .+ 1
        X = convert(Matrix, transpose(data[2:end, 1:end]))
        return X, Y
    end

    if name == "iris"
        iris = DelimitedFiles.readdlm(joinpath(data_path, "iris.csv"), ',')
        X = iris[:, 1:4]
        Y = iris[:, 5]
        return X, Y
    end

    if name == "adult"
        adult = DelimitedFiles.readdlm(joinpath(data_path, "adult.csv"), ',');
        X = adult[:, 1:14];
        Y = adult[:, 15];
        return X, Y
    end

    data_path = joinpath(dirname(pathof(CausalForest)), "..", "data/")

    if name == "causal"
        data = DelimitedFiles.readdlm(joinpath(data_path, "synthetic_causal.csv"), ';', skipstart=1)
        indices = Int.(data[:, 1])
        X = data[:, 2:11]
        T = Int.(data[:, 12])
        Y = data[:, 13]
        return indices, X, T, Y
    end

    if !(name in datasets)
        throw("Available datasets are $(join(datasets,", "))")
    end
end


text:
#
# This file is part of the Actors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

# those are the internal messages

"""
    Become(bhv)

An asynchronous [`Msg`](@ref) to an actor to change its
behavior to `bhv`.
"""
struct Become <: Msg
    x
end

"""
    Call(arg, from::Link)

A synchronous [`Msg`](@ref) to an actor to execute its
behavior with `arg...` and to send the result as a [`Response`](@ref)
message to `from`.
"""
struct Call <: Msg
    x
    from::Link
end
Call(from) = Call((),from)

"""
    Cast(arg)

An asynchronous [`Msg`](@ref) to an actor to execute
its behavior with `arg...` without sending a response.

If the actor is set to `state` dispatch, it updates its internal
state with the result.
"""
struct Cast <: Msg
    x
end
Cast() = Cast(())

"""
    Connect(x, remove=false)

A [`Msg`](@ref) to an actor to connect with `x`. If
`remove=true`, an existing connection gets removed.
"""
struct Connect <: Msg
    x
    remove::Bool
end
Connect(x) = Connect(x, false)

"""
    Diag(x, from::Link)

A synchronous [`Msg`](@ref) to an actor to send diagnostic
information.
"""
struct Diag <: Msg
    x
    from::Link
end

"""
    Down(from, reason, task)

A [`Msg`](@ref) to a monitor actor indicating that an error
has occurred or a [`Exit`](@ref) has been received.
"""
struct Down{T,U} <: Msg
    from::Union{Link,Task}
    reason::T
    task::U
end

"""
    Exec(func::Bhv, from::Link)

A synchronous [`Msg`](@ref) to an actor to execute `func`
and to send a `Response` message with the return value to `from`.
"""
struct Exec <: Msg
    func
    from::Link
end
Exec(t::Tuple, from::Link) = Exec(first(t), from)

"""
    Exit(reason, from, task, state)

A [`Msg`](@ref) to an actor causing it to terminate.
`Exit` messages are sent to [`connect`](@ref)ed actors
if an error has occurred and then are propagated further.
They are not propagated by `:sticky` actors, see
[`trapExit`](@ref).
"""
struct Exit{T,U,V} <: Msg
    reason::T
    from
    task::U
    state::V
end

"""
    Init(f::Bhv)

A [`Msg`](@ref) to an actor to execute the given
[`Bhv`](@ref) and to store it in the [`_ACT`](@ref)
variable.
"""
struct Init <: Msg
    x::Bhv
end

"""
    Query(s::Symbol, from::Link)

A [`Msg`](@ref) to an actor to send a
`Response` message with one of its internal state
variables `s` to `from`.

- `s::Symbol` can be one of `:sta`, `:res`, `:bhv`, `:dsp`.
"""
struct Query <: Msg
    x::Symbol
    from::Link
end
Query(t::Tuple, from::Link) = Query(first(t), from)

"""
    Term(x::Bhv)

A [`Msg`](@ref) to an actor to save the given [`Bhv`](@ref)
and to execute it upon termination.
"""
struct Term <: Msg
    x::Bhv
end

"""
    Timeout()

A return value to signal that a timeout has occurred.
"""
struct Timeout <: Msg end

"""
    Update(s::Symbol, x)

An asynchronous [`Msg`](@ref) to an actor to update its
internal state `s` to `x`.

- `s::Symbol` can be one of `:arg`, `:self`, `:sta`, `:usr`.
"""
struct Update <: Msg
    s::Symbol
    x
end


text:
using StatsBase, Distributions, Plots; pyplot()

function rouletteSpins(p)
    x = 0
    while true
        x += 1
        if rand() < p
            return x
        end
    end
end

p, xGrid, N = 18/37, 1:7, 10^6

mcEstimate = counts([rouletteSpins(p) for _ in 1:N],xGrid)/N

gDist = Geometric(p)
gPmf = [pdf(gDist,x-1) for x in xGrid]

plot(xGrid, mcEstimate, line=:stem, marker=:circle,
	c=:blue, ms=10, msw=0, lw=4, label="MC estimate")
plot!( xGrid, gPmf, line=:stem, marker=:xcross,
	c=:red, ms=6, msw=0, lw=2, label="PMF",
	ylims=(0,0.5), xlabel="x", ylabel="Probability")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using Documenter, Stheno

DocMeta.setdocmeta!(
    Stheno,
    :DocTestSetup,
    :(using Stheno, Random, LinearAlgebra);
    recursive=true,
)

makedocs(
	modules = [Stheno],
    format = Documenter.HTML(),
    sitename = "Stheno.jl",
    pages = [
        "Home" => "index.md",
        "Getting Started" => "getting_started.md",
        "GP API" => "gp_api.md",
        "CompositeGP API" => "composite_gp_api.md",
        "Input Types" => "input_types.md",
        "Kernel Design" => "kernel_design.md",
        "Internals" => "internals.md",
    ],
)

deploydocs(repo="github.com/willtebbutt/Stheno.jl.git")


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

Random.seed!(100)
## Convergence Testing
dts = (1/2) .^ (7:-1:4)
testTol = 0.2

prob = prob_ode_linear
sol = solve(prob,OwrenZen3())
@test length(sol) < 20
sol = solve(prob,OwrenZen4())
@test length(sol) < 20
sol = solve(prob,OwrenZen5())
@test length(sol) < 20

sim = test_convergence(dts,prob,OwrenZen3(),dense_errors=true)
@test abs(sim.𝒪est[:final]-3) < testTol
@test abs(sim.𝒪est[:L2]-3) < testTol
sim = test_convergence(dts,prob,OwrenZen4(),dense_errors=true)
@test abs(sim.𝒪est[:final]-4) < testTol
@test abs(sim.𝒪est[:L2]-4) < testTol
sim = test_convergence(dts,prob,OwrenZen5(),dense_errors=true)
@test abs(sim.𝒪est[:final]-5) < testTol
@test abs(sim.𝒪est[:L2]-5) < testTol

prob = prob_ode_2Dlinear
sol = solve(prob,OwrenZen3())
@test length(sol) < 20
sol = solve(prob,OwrenZen4())
@test length(sol) < 20
sol = solve(prob,OwrenZen5())
@test length(sol) < 20

sim = test_convergence(dts,prob,OwrenZen3(),dense_errors=true)
@test abs(sim.𝒪est[:final]-3) < testTol
@test abs(sim.𝒪est[:L2]-3) < testTol
sim = test_convergence(dts,prob,OwrenZen4(),dense_errors=true)
@test abs(sim.𝒪est[:final]-4) < testTol
@test abs(sim.𝒪est[:L2]-4) < testTol
sim = test_convergence(dts,prob,OwrenZen5(),dense_errors=true)
@test abs(sim.𝒪est[:final]-5) < testTol
@test abs(sim.𝒪est[:L2]-5) < testTol


text:
function cnsactivations!(system::CVSystem,n::Int64)
    # activation based on Gompertz's sigmoidal function
    push!(system.cns.np,system.cns.a*exp(
        -system.cns.b*exp(-system.cns.c*
        (system.cns.Paverage[system.solverparams.numbeats]/
        system.cns.Ptarget - 1))));
    push!(system.cns.ns,system.cns.a*exp(
        -system.cns.b*exp(system.cns.c*
        (system.cns.Paverage[system.solverparams.numbeats]/
        system.cns.Ptarget - 1))));
    # # hard-code activations
    # nst = 0.25;
    # r = 1/system.cns.c*log(log(nst^(-1/system.cns.b)))+1;
    # ns = exp(-system.cns.b*exp(system.cns.c*(r-1)));
    # np = exp(-system.cns.b*exp(-system.cns.c*(r-1)));
    # push!(system.cns.ns,ns);
    # push!(system.cns.np,np);
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling