text:
@testset "ressource constraints" begin
    item1 = Item(Inventory(1, 1000), Supplier(100,1), name = "item1", policy = QPolicy())
    item2 = Item(Inventory(1, 1000), Supplier(100,1), name = "item2", policy = QPolicy())
    item3 = Item(Inventory(1, 1000), Supplier(100,1), name = "item3", policy = QPolicy())
    item4 = Item(Inventory(1, 1000), Supplier(100,1), name = "item4", policy = QPolicy())
    item5 = Item(Inventory(1, 1000), Assembly(100,1, item1 => 2, item2 => 4), name = "item5", policy = QPolicy())
    item6 = Item(Inventory(1, 1000), Assembly(100,1, item3 => 10, item4 => 1), name = "item6", policy = QPolicy())
    item7 = EndProduct(
        Market(5, CVNormal{0.0}, 5, 0, 10),
        Inventory(50, 1000),
        Assembly(100,1, item5 => 1, item6 => 2),
        policy = QPolicy(),
        name = "item7"
        )

    bom = [item1,item2,item3,item4,item5,item6,item7]
    hours_cons = RessourceConstraint(1680, item5.sources[1] => 5, item6.sources[1] => 1, item7.sources[1] => 15)
    is = InventorySystem(6, bom, [hours_cons])
    is([10,10,0,0,20,0,0])
    @test state(item1) == [1000-40]
    @test state(item2) == [1000-80]
    @test state(item3) == [1000-100]
    @test state(item4) == [1000-10]
    @test state(item5) == [1000-10+20]
    @test state(item6) == [1000-20+10]
    @test state(item7.inventory) == [1000-10+10]
    hours_cons = hours_cons = RessourceConstraint(168, item5.sources[1] => 5, item6.sources[1] => 1, item7.sources[1] => 15)
    is = InventorySystem(6, bom, [hours_cons])
    reset!(is)
    is([10,10,0,0,20,0,0])
    @test state(item1) != [1000-40]
    @test state(item2) != [1000-80]
    @test state(item3) != [1000-100]
    @test state(item4) != [1000-10]
    @test state(item5) != [1000-10+20]
    @test state(item6) != [1000-20+10]
    @test state(item7.inventory) != [1000-10+10]
    r = 168/(15*10+10+5*20)
    @test state(item5) == [1000-10r+20r]
    @test state(item6) == [1000-20r+10r]
    @test state(item7.inventory) == [1000-10+10r]
end


text:
#==============================================================================#
# AWSSDB/test/runtests.jl
#
# Copyright Sam O'Connor 2016 - All rights reserved
#==============================================================================#

using AWSSDB
using Base.Test

AWSCore.set_debug_level(1)

#-------------------------------------------------------------------------------
# Load credentials...
#-------------------------------------------------------------------------------

aws = AWSCore.aws_config(region="ap-southeast-2")

#-------------------------------------------------------------------------------
# SDB tests
#-------------------------------------------------------------------------------

db = "aws_sdb_jl_test_domain"
simpledb = SimpleDB(aws, db)

if db in sdb_list_domains(aws)
    sdb_delete_domain(aws, db)
    while db in sdb_list_domains(aws)
        println("Waiting for delete...")
        sleep(1)
    end
end

sdb_create_domain(aws, db)

while !(db in sdb_list_domains(aws))
    println("Waiting for create...")
    sleep(1)
end

data = Dict("key1" => Pair["a1" => "1"],
            "key2" => Pair["a1" => "21", "a2" => "22"],
            "key3" => Pair["a1" => "31", "a2" => "32", "a3" => "33"])

for (k, a) in data
    simpledb[k] = a
end

sleep(1)
while length(collect(sdb_select(aws, "select itemName() from $db"))) < 3
    println("Waiting for put...")
    sleep(1)
end

@test Dict(sdb_get(aws, db, "key1")) == Dict(data["key1"])
@test Dict(sdb_get(aws, db, "key2")) == Dict(data["key2"])
@test Dict(sdb_get(aws, db, "key3")) == Dict(data["key3"])

@test Dict(simpledb["key1"]) == Dict(data["key1"])
@test Dict(simpledb["key2"]) == Dict(data["key2"])
@test Dict(simpledb["key3"]) == Dict(data["key3"])

sdb_put(aws, db, "key1", Pair["a1" => 2], replace=false)

sleep(1)
while length(Dict(sdb_get(aws, db, "key1"))["a1"]) < 2
    println("Waiting for put...")
    sleep(1)
end

@test sdb_get(aws, db, "key1") == Pair["a1" => ["1", "2"]]

@test Dict(sdb_select(aws, "select `a1` from $db")) ==
    Dict("key1" => Pair["a1" => ["1", "2"]],
         "key2" => Pair["a1" => "21"],
         "key3" => Pair["a1" => "31"])

@test Dict(sdb_select(aws, db, ["a1"])) ==
    Dict("key1" => Pair["a1" => ["1", "2"]],
         "key2" => Pair["a1" => "21"],
         "key3" => Pair["a1" => "31"])

@test Dict(sdb_select(aws, db, ["a1", "a2"])) ==
    Dict("key1" => Pair["a1" => ["1", "2"], "a2" => ""],
         "key2" => Pair["a1" => "21", "a2" => "22"],
         "key3" => Pair["a1" => "31", "a2" => "32"])

sdb_delete_item(aws, db, "key1")
sdb_delete_item(aws, db, "key2")
sdb_delete_item(aws, db, "key3")

while length(collect(sdb_select(aws, "select itemName() from $db"))) > 0
    println("Waiting for delete...")
    sleep(1)
end

#==============================================================================#
# End of file.
#==============================================================================#

ext: jl


text:
"""
    iterate_jkp(Fp,Fi, top_k = Val(false))

Iterable such that
```
for (j,kps) in iterate_jkp(Fp,Fi)
    for (k,ps) in kps
        for i in Fi[ps]
```
is equivalent to
```
for j = 1:n
    for k = 1:j-1 if F[k,j] != 0 && (!top_k || all(F[k+1:j-1,k] .== 0))
        for i = j:n if F[i,k] != 0
```
The body of the `j`-loop may update the sparsity pattern of `F[:,j:n]`.

Iterator creation requires `O(n)` runtime and allocates two `Vector{Int}(n)`.
Each call to `iterate()` requires `O(1)` runtime and is allocation-free.

`issorted(Fi[Fp[j]:Fp[j+1]-1])` must hold true for all `j in 1:n`. The
iterator may produce incorrect results if this property is violated.
"""
iterate_jkp(Fp,Fi, top_k = Val(false)) = Iterate_jkp(Fp,Fi, top_k)

#=
This implementation is based on

  [ESS81]
    Eisenstat, Schultz, Sherman
    Algorithms and Data Structures for Sparse Symmetric Gaussian Elimination
    SIAM Journal on Scientific and Statistical Computing (1981)

Meaning of variables:

  kl[i] with i >= j points to a linked-list representation of the set
    { k < j | F[i,k] != 0 && all(F[j+1:i-1,k] .== 0) }
  Note that this set simplifies to { k < j | F[j,k] != 0 } if i == j.
  End of list is indicated by kl[i] = 0.

  rp is a vector of length n such that
    ( rp[k] : Fp[k+1]-1 ) = nonzeros in F[j:n,k] for all k < j

The correspondences to the variables in [ESS81] are kl -> JL and rp -> IL.
=#

struct Iterate_jkp{TopK}
    Fp::Vector{Int}
    Fi::Vector{Int}
    kl::Vector{Int}
    rp::Vector{Int}
end

Base.length(jkp::Iterate_jkp) = length(jkp.Fp)-1

struct Iterate_kp{TopK}
    Fp::Vector{Int}
    Fi::Vector{Int}
    kl::Vector{Int}
    rp::Vector{Int}
    j::Int
end

Base.IteratorSize(::Type{<:Iterate_kp}) = Base.SizeUnknown()

function Iterate_jkp(Fp,Fi, top_k = Val(false))
    n = length(Fp)-1

    # Allocate workspace
    kl = fill(0,n)
    rp = Vector{Int}(undef,n)

    return Iterate_jkp{top_k}(Fp,Fi,kl,rp)
end

function Base.iterate(jkp::Iterate_jkp{T}) where {T}
    Fp = jkp.Fp
    Fi = jkp.Fi
    kl = jkp.kl
    rp = jkp.rp
    length(Fp) <= 1 && return nothing

    # Iterate over first j
    ks = Iterate_kp{T}(Fp,Fi,kl,rp,1)
    return (1,ks), 1
end

function Base.iterate(jkp::Iterate_jkp{T},j) where {T}
    Fp = jkp.Fp
    Fi = jkp.Fi
    kl = jkp.kl
    rp = jkp.rp
    length(Fp)-1 == j && return nothing

    # Update kl and rp for finished j and move to next j
    for p in @inbounds(Fp[j]:Fp[j+1]-1)
        i = Fi[p]
        if i > j
            @inbounds kl[i],kl[j] = j,kl[i]
            @inbounds rp[j] = p
            break
        end
    end
    j += 1

    # Iterate over next j
    ks = Iterate_kp{T}(Fp,Fi,kl,rp,j)
    return (j,ks),j
end

function Base.iterate(kp::Iterate_kp)
    Fp = kp.Fp
    Fi = kp.Fi
    kl = kp.kl
    rp = kp.rp
     j = kp.j

    # Iterate over first k
    @inbounds k = kl[j]
    k == 0 && return nothing
    return (k, @inbounds(rp[k]:Fp[k+1]-1)),k
end

function Base.iterate(kp::Iterate_kp{Val(false)}, k)
    Fp = kp.Fp
    Fi = kp.Fi
    kl = kp.kl
    rp = kp.rp
    j = kp.j
    n = length(Fp)-1

    # Update kl and rp for finished k and move to next k
    rpk = @inbounds( rp[k] += 1 )
    kk = @inbounds kl[k]
    if rpk < @inbounds Fp[k+1]
        i = Fi[rpk]
        @assert 1 <= i <= n
        @inbounds kl[i],kl[k] = k,kl[i]
    end
    k = kk

    # Iterate over next k
    k == 0 && return nothing
    return (k, @inbounds(rp[k]:Fp[k+1]-1)),k
end

function Base.iterate(kp::Iterate_kp{Val(true)}, k)
    Fp = kp.Fp
    Fi = kp.Fi
    kl = kp.kl
    rp = kp.rp
     j = kp.j

    # Iterate over next k
    k = @inbounds kl[k]
    k == 0 && return nothing
    return (k, @inbounds(rp[k]:Fp[k+1]-1)),k
end

max_stars_repo_name: ettersi/SparseFactorizations

max_issues_repo_name: ettersi/SparseFactorizations

max_forks_repo_name: ettersi/SparseFactorizations


text:
module DB
  using Catlab: @present
  using Catlab.Present
  using Catlab.CategoricalAlgebra.CSets
  using Catlab.CSetDataStructures: struct_acset
  export TheorySQL, generate_schema_sql, @present, get_fields, TypeToSQL, typeToSQL, @db_schema, AbstractSQL

  TypeToSQL = Dict("String" => "text",
                   "Int" => "int",
                   "Int64" => "int",
                   "Real" => "real",
                   "Bool" => "boolean")

  typeToSQL(x) = TypeToSQL[string(x)]
  @present TheorySQL(FreeSchema) begin
    Int::AttrType
    Int64::AttrType
    Real::AttrType
    String::AttrType
    Bool::AttrType
  end;

  @abstract_acset_type AbstractSQL

  # TODO: This should be replacable with a cleaner method
  macro db_schema(head)
    struct_name = gensym()
    quote
      $(esc(:eval))(struct_acset($(Meta.quot(struct_name)), AbstractSQL, $(esc(head.args[2]))))
      $(esc(head.args[1]))() = $(esc(struct_name)){$(esc(Int)), $(esc(Int)), $(esc(Real)), $(esc(String)), $(esc(Bool))}()
    end
  end

  function generate_schema_sql(schema::AbstractSQL)
    queries = map(collect(get_fields(schema))) do (name, col)
      cols = ["$n $(typeToSQL(t))" for (n,t) in col]
      "CREATE TABLE $name ($(join(cols, ", ")))"
    end
    string(join(queries, ";\n"), ";")
  end

  function get_fields(schema::AbstractSQL)
    fields = Dict{Symbol, Array{Tuple{Symbol, Type},1}}()
    for (name, table) in pairs(tables(schema))
      table_name = name

      # Get the column names and types
      col_names = propertynames(table)
      types = eltype.([schema[c] for c in col_names])
      col_names = map(x -> Symbol(split(string(x), r"_\d+_")[end]), col_names)
      fields[table_name] = map(zip(col_names,types)) do (n,t)
        (n, t)
      end
    end
    fields
  end

  # TODO: Maybe add in a way to get_fields from a Presentation?
  #       Allows more convenience when working with workflows,
  #       but doesn't allow you to take advantage of the CSet
  #       structure
end


text:
function RGP(y::REALVEC,
             x::REALMAT,
             K::Kernel,
             M::Integer = size(y),
             s2n::Real = 1)

    @assert M <= length(y) && M > 0 "Budget size, M, must be [1;T]."
    @assert s2n > 0 "Output noise, s2n, must be positive."

    T, D = size(x)

    basis = collect(1:M)

    # Run standard GP
    μ, C, _ = GP(y[basis], x[basis, :], x[basis, :], K, s2n)

    m = zeros(M)
    #C = diagm(ones(D) * 1_000)

    for t in 1:T

        # Calculate kernel
        ktt = kernel(x[1:t, :], x[1:t, :], K)
        kbt = kernel(x[basis, :], x[1:t, :], K)
        kbb = kernel(x[basis, :], x[basis, :], K)
        # Calculate gain matrix J
        J = inv(kbb) * kbt
        # Calculate mean
        #μ = m .+ J .* (μ - m)
        # calculate covariance matrices
        B = ktt - J' * kbt
        println("Size of J: $(size(J))")
        println("Size of C: $(size(C))")
        println("Size of B: $(size(B))")
        CP = B + J' * C * J
#        C = [C (C * J'); (J * C) ]

    end

    return (predictions = μ, variances = 0, loglikelihood = 0)
end

#function RGP(y::REALVEC,
#             x::REALVEC,
#             K::Kernel,
#             M::Integer = size(y),
#             s2n::Real = 1)
#   RGP(y, reshape(x, :, 1), K, s2n)
#end


text:
function slice(A::Vector{T},p::Vector{Int}) ::Vector{T} where T
    srez = Vector{Int}(undef,size(p))
    for i in 1:size(p)
        srez[i]=A[i]
    end
    return srez
end

function slice(A::Vector{T},p::Vector{Int}) ::Vector{T} where T
    return A[p]
end


text:
let p2pdevs = Set()

    global checkP2P, enableP2P

    """
    `checkP2P(did1::Int, did2::Int)::Bool` Checks for peer access between devices did1 and did2.
    Returns false if p2p access is not supported and true otherwise.
    The function checks peer access of both did1 to did2 and did2 to did1; so checkP2P(d1, d2)
    is identical to checkP2P(d2, d1).
    """
    function checkP2P(did1::Int, did2::Int)::Bool
        access = Cint[0]
        @cudart(cudaDeviceCanAccessPeer, (Ptr{Cint}, Cint, Cint), access, did1, did2)
        res1 = Bool(access[1])
        @cudart(cudaDeviceCanAccessPeer, (Ptr{Cint}, Cint, Cint), access, did2, did1)
        res2 = Bool(access[1])
        res1 && res2
    end

    """
    `enableP2P(did1::Int, did2::Int)` Enables peer access between devices did1 and did2.
    The function enables peer access of both did1 to did2 and did2 to did1; so enableP2P(d1, d2)
    is identical to enableP2P(d2, d1).
    """
    function enableP2P(did1::Int, did2::Int)
        if ((did1, did2) in p2pdevs) || ((did2, did1) in p2pdevs)
            return p2pdevs #for consistency
        end
        gpu_temp = gpu()
        gpu(did1)
        @cudart(cudaDeviceEnablePeerAccess, (Cint, Cint), did2, 0)
        gpu(did2)
        @cudart(cudaDeviceEnablePeerAccess, (Cint, Cint), did1, 0)
        gpu(gpu_temp)
        push!(p2pdevs, (did1, did2))
    end

    function all_pairs(f::Function, gpuList)
        ngpus = length(gpuList)
        for i = 1:ngpus-1
            for j = (i+1):ngpus
                f(gpuList[i], gpuList[j])
            end
        end
    end

    """
    `enableP2P(gpuList::Union{Array{Int, 1}, Nothing}=nothing)::Bool`
    Checks and enables peer access between all gpus whose ids are in the `gpuList`.
    Returns `false` if peer access does not exist between any gpu pair, `true` otherwise.
    If `gpuList =  nothing`, `[0:gpuCount()-1]` is used as the gpuList.
    """
    function enableP2P(gpuList::Union{Array{Int, 1}, Nothing}=nothing;
                       verbose::Bool=false)::Bool
        if gpuList == nothing
            gpuList = Array{Int, 1}(0:(gpuCount()-1))
        end
        # check the access status
        check = true
        all_pairs(gpuList) do d1, d2
            check = check && checkP2P(d1, d2)
            verbose && ~check && warn("No p2p access between $d1 and $d2")
        end
        if ~check; return false; end
        # enable access between peers
        all_pairs(gpuList) do d1, d2
            enableP2P(d1, d2)
        end
        return true
    end

end


text:
using Surrogates
using QuasiMonteCarlo
using QuasiMonteCarlo: KroneckerSample, SectionSample, GoldenSample
using Distributions: Cauchy, Normal
using Test

#1D
lb = 0.0
ub = 5.0
n = 5
d = 1

## Sampling methods from QuasiMonteCarlo.jl ##

# GridSample
s = Surrogates.sample(n,lb,ub,GridSample(0.1))
@test s isa Vector{Float64} && length(s) == n && all(x -> lb ≤ x ≤ ub, s)

# UniformSample
s = Surrogates.sample(n,lb,ub,UniformSample())
@test s isa Vector{Float64} && length(s) == n && all(x -> lb ≤ x ≤ ub, s)

# SobolSample
s = Surrogates.sample(n,lb,ub,SobolSample())
@test s isa Vector{Float64} && length(s) == n && all(x -> lb ≤ x ≤ ub, s)

# LatinHypercubeSample
s = Surrogates.sample(n,lb,ub,LatinHypercubeSample())
@test s isa Vector{Float64} && length(s) == n && all(x -> lb ≤ x ≤ ub, s)

# LowDiscrepancySample
s = Surrogates.sample(20,lb,ub,LowDiscrepancySample(10))
@test s isa Vector{Float64} && length(s) == 20 && all(x -> lb ≤ x ≤ ub, s)

# LatticeRuleSample (not originally in Surrogates.jl, now available through QuasiMonteCarlo.jl)
s = Surrogates.sample(20,lb,ub,LatticeRuleSample())
@test s isa Vector{Float64} && length(s) == 20 && all(x -> lb ≤ x ≤ ub, s)

# Distribution sampling (Cauchy)
s = Surrogates.sample(n,d,Cauchy())
@test s isa Vector{Float64} && length(s) == n

# Distributions sampling (Normal)
s = Surrogates.sample(n,d,Normal(0,4))
@test s isa Vector{Float64} && length(s) == n

## Sampling methods specific to Surrogates.jl ##

# KroneckerSample
s = Surrogates.sample(n,lb,ub,KroneckerSample(sqrt(2),0))
@test s isa Vector{Float64} && length(s) == n && all(x -> lb ≤ x ≤ ub, s)

# GoldenSample
s = Surrogates.sample(n,lb,ub,GoldenSample())
@test s isa Vector{Float64} && length(s) == n && all(x -> lb ≤ x ≤ ub, s)

# SectionSample
constrained_val = 1.0
s = Surrogates.sample(n, lb, ub, SectionSample([NaN64], UniformSample()))
@test s isa Vector{Float64} && length(s) == n && all(x -> lb ≤ x ≤ ub, s)

s = Surrogates.sample(n, lb, ub, SectionSample([constrained_val], UniformSample()))
@test s isa Vector{Float64} && length(s) == n && all(x -> lb ≤ x ≤ ub, s)
@test all(==(constrained_val), s)

# ND
# Now that we use QuasiMonteCarlo.jl, these tests are to make sure that we transform the output
# from a Matrix to a Vector of Tuples properly for ND problems.

lb = [0.1,-0.5]
ub = [1.0,20.0]
n = 5
d = 2

#GridSample{T}
s = Surrogates.sample(n,lb,ub,GridSample([0.1,0.5]))
@test isa(s,Array{Tuple{typeof(s[1][1]),typeof(s[1][1])},1}) == true

#UniformSample()
s = Surrogates.sample(n,lb,ub,UniformSample())
@test isa(s,Array{Tuple{typeof(s[1][1]),typeof(s[1][1])},1}) == true

#SobolSample()
s = Surrogates.sample(n,lb,ub,SobolSample())
@test isa(s,Array{Tuple{typeof(s[1][1]),typeof(s[1][1])},1}) == true

#LHS
s = Surrogates.sample(n,lb,ub,LatinHypercubeSample())
@test isa(s,Array{Tuple{typeof(s[1][1]),typeof(s[1][1])},1}) == true

#LDS
s = Surrogates.sample(n,lb,ub,LowDiscrepancySample([10,3]))
@test isa(s,Array{Tuple{typeof(s[1][1]),typeof(s[1][1])},1}) == true

#Distribution 1
s = Surrogates.sample(n,d,Cauchy())
@test isa(s,Array{Tuple{typeof(s[1][1]),typeof(s[1][1])},1}) == true

#Distribution 2
s = Surrogates.sample(n,d,Normal(3,5))
@test isa(s,Array{Tuple{typeof(s[1][1]),typeof(s[1][1])},1}) == true

#Kronecker
s = Surrogates.sample(n,lb,ub,KroneckerSample([sqrt(2),3.1415],[0,0]))
@test isa(s,Array{Tuple{typeof(s[1][1]),typeof(s[1][1])},1}) == true

#Golden
s = Surrogates.sample(n,lb,ub,GoldenSample())
@test isa(s,Array{Tuple{typeof(s[1][1]),typeof(s[1][1])},1}) == true

# SectionSample
constrained_val = 1.0
s = Surrogates.sample(n, lb, ub, SectionSample([NaN64, constrained_val], UniformSample()))
@test all(x -> x[end] == constrained_val, s)
@test isa(s,Array{Tuple{typeof(s[1][1]),typeof(s[1][1])},1}) == true
@test all(x -> lb[1] ≤ x[1] ≤ ub[1], s)

max_stars_repo_name: ludoro/Surrogates

max_issues_repo_name: ludoro/Surrogates

max_forks_repo_name: ludoro/Surrogates


text:
using ThreeOmegaMethod

include("test3omega.jl")

max_stars_repo_name: lnacquaroli/ThreeOmegaMethod

max_issues_repo_name: lnacquaroli/ThreeOmegaMethod


text:
module URDF

using RigidBodyDynamics
using LightXML
using StaticArrays
using Rotations
using DocStringExtensions
using RigidBodyDynamics.Graphs

using RigidBodyDynamics: Bounds, upper, lower
using RigidBodyDynamics: has_loops, joint_to_predecessor
using RigidBodyDynamics: DEFAULT_GRAVITATIONAL_ACCELERATION
using LinearAlgebra: ×

export
    parse_urdf,
    write_urdf

include("parse.jl")
include("write.jl")

end # module


text:
struct ScaledDiagonallyDominantVariableBridge{T} <: AbstractVariableBridge
    side_dimension::Int
    variables::Vector{NTuple{3, MOI.VariableIndex}}
    psd2x2::Vector{MOI.ConstraintIndex{MOI.VectorOfVariables,
                                       PositiveSemidefinite2x2ConeTriangle}}
end

function add_variable_bridge(
    ::Type{ScaledDiagonallyDominantVariableBridge{T}}, model::MOI.ModelLike,
    s::ScaledDiagonallyDominantConeTriangle) where {T}
    # `p.Q` is SDD iff it is the sum of psd matrices Mij that are zero except
    # for entries ii, ij and jj [Lemma 9, AM17].
    n = s.side_dimension
    @assert n > 1 # The bridges does not work with 1, `matrix_cone`
                  # should have returned `Nonnegatives` instead
    # `g[r, c]` will contain the expression `sum Mij[r, c]`.
    F = MOI.ScalarAffineFunction{T}
    g = F[zero(F) for i in 1:MOI.dimension(s)]
    num_off_diag = MOI.dimension(s) - n
    variables = Vector{NTuple{3, MOI.VariableIndex}}(undef, num_off_diag)
    psd2x2 = Vector{MOI.ConstraintIndex{
        MOI.VectorOfVariables, PositiveSemidefinite2x2ConeTriangle}}(
            undef, num_off_diag)
    diag_idx(i) = div(i * (i + 1), 2)
    k = 0
    k2x2 = 0
    for j in 1:n
        for i in 1:(j-1)
            k += 1
            k2x2 += 1
            vii, vij, vjj = MOI.add_variables(model, 3)
            variables[k2x2] = (vii, vij, vjj)
            Mii = MOI.SingleVariable(vii)
            MOIU.operate!(+, T, g[diag_idx(i)], Mii)
            Mij = MOI.SingleVariable(vij)
            MOIU.operate!(+, T, g[k], Mij)
            Mjj = MOI.SingleVariable(vjj)
            MOIU.operate!(+, T, g[diag_idx(j)], Mjj)
            # PSD constraints on 2x2 matrices are SOC representable
            psd2x2[k2x2] = MOI.add_constraint(
                model, MOI.VectorOfVariables([vii, vij, vjj]),
                PositiveSemidefinite2x2ConeTriangle())
        end
        k += 1 # diagonal entry `(i, i)`
    end
    return g, ScaledDiagonallyDominantVariableBridge{T}(n, variables, psd2x2)
end

function MOIB.added_constraint_types(
    ::Type{ScaledDiagonallyDominantVariableBridge{T}}) where {T}
    return [(MOI.VectorOfVariables, PositiveSemidefinite2x2ConeTriangle)]
end

function variable_bridge_type(::Type{ScaledDiagonallyDominantConeTriangle},
                              T::Type)
    return ScaledDiagonallyDominantVariableBridge{T}
end

# Attributes, VariableBridge acting as an model
function MOI.get(bridge::ScaledDiagonallyDominantVariableBridge,
                 ::MOI.NumberOfVariables)
    return 3 * length(bridge.variables)
end
function MOI.get(
    bridge::ScaledDiagonallyDominantVariableBridge,
    ::MOI.NumberOfConstraints{MOI.VectorOfVariables,
                              ScaledDiagonallyDominantConeTriangle})
    return length(bridge.psd2x2)
end
function MOI.get(
    bridge::ScaledDiagonallyDominantVariableBridge,
    ::MOI.ListOfConstraintIndices{MOI.VectorOfVariables,
                                  ScaledDiagonallyDominantConeTriangle})
    return bridge.psd2x2
end

# Indices
function MOI.delete(model::MOI.ModelLike,
                    bridge::ScaledDiagonallyDominantVariableBridge)
    for ci in bridge.psd2x2
        MOI.delete(model, ci)
    end
    for (vii, vij, vjj) in bridge.variables
        MOI.delete(model, vii)
        MOI.delete(model, vij)
        MOI.delete(model, vjj)
    end
end


text:
#=
Author:= Shreeda Bhat
License:= MIT
UTF-8 =#

#=
Julia  has same data structures that of python

Some data structures are
1.Tuples
2.Dictionary
3.Arrays

Dictionary and Arrays are mutable but Tuples are immutable
=#

# Let's see Dictionary

println("########################### Dictionary #################################")

my_dict = Dict("name"=>"shreeda","age"=>20,"study"=>"Julia")
println(my_dict)

#let's add a entry to existing Dict
my_dict["work"]="Astaqc Consulting"
println(my_dict)

#let's extract work from my_dict

println(my_dict["work"])

#We can delete a record from dict by using `pop!` function

deleted_entry = pop!(my_dict,"work")
println(my_dict)

println("########################### Tuples #################################")

tuple_ex = ("shreeda","20","Astaqc Consulting")

println(tuple_ex)
# We can retrieve the record by using this syntax variable_name[key]

println(tuple_ex[1])
#= since tuples are immutable we cant update a record
    =#

println("###########################  Arrays #################################")

Arrays_mix_ex = ["Shreeda",20,"Astaqc Consulting"] # returns type Any Array
Arrays_int_ex = [1,2,3,4,5,6,7]
println("This is mix array $Arrays_mix_ex")
println("This is integer array $Arrays_int_ex")

#=
We can delete and insert element by `pop` and `push` functions
=#
deleted_ele = pop!(Arrays_int_ex)

println("Now deleted array $Arrays_int_ex")

#push function

push_ele = push!(Arrays_mix_ex,"Julia")
println("This is updated matrix $Arrays_mix_ex")

#=

Let's do multi dimension array :yay:

=#

two_d_int_array = [[1,2,3],[4,5,6],[7,8,9]]

two_d_string_array = [["shreeda","astaqc","julia"],["python","java","C"]]

two_d_mix_array = [["shreeda",20,"astaqc"],["erin",21,"astaqc"]]

println("this is multi dimension integer array $two_d_int_array")
println("this is multi dimension string array $two_d_string_array")
println("this is multi dimension  mixed array $two_d_mix_array")


text:
using Base.Test, Whittaker

@testset "whittaker1" begin
  @testset "out-of-place" begin
    tol = 0.0001

    y = [1.0, 1.2, 0, 1.6, 1.8]
    w = [1.0, 1.0, 0, 1.0, 1.0]
    lambda = 0.0001
    z = whittaker1(y, w, lambda)
    @test_approx_eq_eps z [1.00002, 1.2, 1.4, 1.6, 1.79998] tol

    y = [1., 2, 0, 4, 5]
    w = [1, 1, 0, 1, 1]
    lambda = 0.0001
    z = whittaker1(y, w, lambda)
    @test_approx_eq_eps z [1.0001, 2.0, 3.0, 4.0, 4.9999] tol
  end

  @testset "in-place" begin
    tol = 0.0001

    y = [1.0, 1.2, 0, 1.6, 1.8]
    w = [1.0, 1.0, 0, 1.0, 1.0]
    lambda = 0.0001
    whittaker1!(y, w, lambda)
    @test_approx_eq_eps y [1.00002, 1.2, 1.4, 1.6, 1.79998] tol

    y = [1., 2, 0, 4, 5]
    w = [1, 1, 0, 1, 1]
    lambda = 0.0001
    whittaker1!(y, w, lambda)
    @test_approx_eq_eps y [1.0001, 2.0, 3.0, 4.0, 4.9999] tol
  end
end


text:
module GenericTest

using Test

using Avro.Schemas
using Avro.Io
using Avro.Generic

const TEST_RECORD_SCHEMA =
    Schemas.RecordSchema(
        Schemas.FullName("test"),
        [
            Schemas.Field("a", 0, Schemas.LONG),
            Schemas.Field("b", 1, Schemas.STRING)
        ]
    )

const RECORD_EXAMPLES =
    [
        (
            GenericRecord(TEST_RECORD_SCHEMA, [Int64(27), "foo"]),
            [0x36, 0x06, 0x66, 0x6f, 0x6f]
        )
    ]

const TEST_ENUM_SCHEMA =
    Schemas.EnumSchema(Schemas.FullName("Foo"), ["A", "B", "C", "D"])

const ENUM_EXAMPLES =
    [
        (GenericEnumSymbol(TEST_ENUM_SCHEMA, "A"), [0x00]),
        (GenericEnumSymbol(TEST_ENUM_SCHEMA, "D"), [0x06])
    ]

const TEST_FIXED_SCHEMA = Schemas.FixedSchema(Schemas.FullName("md5"), 2)

const FIXED_EXAMPLES =
    [
        (GenericFixed(TEST_FIXED_SCHEMA, [0x01, 0x02]), [0x01, 0x02])
        (GenericFixed(TEST_FIXED_SCHEMA, [0xAA, 0xBB]), [0xAA, 0xBB])
    ]

const ARRAY_EXAMPLES =
    [
        (
            Int64[3, 27],
            Schemas.ArraySchema(Schemas.LONG),
            [0x04, 0x06, 0x36, 0x00]
        )
    ]

const MAP_EXAMPLES =
    [
        (
            Dict{String, Int64}("bar" => 27, "foo" => 3),
            Schemas.MapSchema(Schemas.LONG),
            [0x04, 0x06, 0x62, 0x61, 0x72, 0x36, 0x06, 0x66, 0x6f, 0x6f, 0x06, 0x00]
        )
    ]

const UNION_EXAMPLES =
    [
        (Int64(2), Schemas.UnionSchema([Schemas.LONG, Schemas.NULL]), [0x00, 0x04]),
        (nothing, Schemas.UnionSchema([Schemas.LONG, Schemas.NULL]), [0x02])
    ]

buffer = IOBuffer()
encoder = BinaryEncoder(buffer)
decoder = BinaryDecoder(buffer)

@testset "Generic writers" begin
    @testset "Record" for (input, expected) in RECORD_EXAMPLES
        # Encode the datum
        bytes_written = write(encoder, input.schema, input)

        # Decode the datum
        seekstart(buffer)
        output = read(decoder, input.schema)

        # Inspect the contents of the buffer
        contents = take!(buffer)

        @test expected == contents
        @test input == output
        @test length(expected) == bytes_written
    end

    @testset "Enum" for (input, expected) in ENUM_EXAMPLES
        # Encode the datum
        bytes_written = write(encoder, input.schema, input)

        # Decode the datum
        seekstart(buffer)
        output = read(decoder, input.schema)

        # Inspect the contents of the buffer
        contents = take!(buffer)

        @test expected == contents
        @test input == output
        @test length(expected) == bytes_written
    end

    @testset "Fixed" for (input, expected) in FIXED_EXAMPLES
        # Encode the datum
        bytes_written = write(encoder, input.schema, input)

        # Decode the datum
        seekstart(buffer)
        output = read(decoder, input.schema)

        # Inspect the contents of the buffer
        contents = take!(buffer)

        @test expected == contents
        @test input == output
        @test length(expected) == bytes_written
    end

    @testset "Array" for (input, schema, expected) in ARRAY_EXAMPLES
        # Encode the datum
        bytes_written = write(encoder, schema, input)

        # Decode the datum
        seekstart(buffer)
        output = read(decoder, schema)

        # Inspect the contents of the buffer
        contents = take!(buffer)

        @test expected == contents
        @test input == output
        @test length(expected) == bytes_written
    end

    @testset "Map" for (input, schema, expected) in MAP_EXAMPLES
        # Encode the datum
        bytes_written = write(encoder, schema, input)

        # Decode the datum
        seekstart(buffer)
        output = read(decoder, schema)

        # Inspect the contents of the buffer
        contents = take!(buffer)

        @test expected == contents
        @test input == output
        @test length(expected) == bytes_written
    end

    @testset "Union" for (input, schema, expected) in UNION_EXAMPLES
        # Encode the datum
        bytes_written = write(encoder, schema, input)

        # Decode the datum
        seekstart(buffer)
        output = read(decoder, schema)

        # Inspect the contents of the buffer
        contents = take!(buffer)

        @test expected == contents
        @test input == output
        @test length(expected) == bytes_written
    end
end

end


text:
##### multi dimensional advection
##### For incompressible model only
##### calculate tendencies in x direction
@kernel function calc_Gcˣ_kernel!(Gc, c, u, g::AbstractGrid, ΔT)
    i, j, k = @index(Global, NTuple)
    ### offset index for halo points
    ii = i + g.Hx
    jj = j + g.Hy
    kk = k + g.Hz
    @inbounds Gc[ii, jj, kk] = adv_flux_x(ii, jj, kk, g, u, c, ΔT)
end
function calc_Gcsˣ!(Gcs, nut, u, g::AbstractGrid, ΔT, arch::Architecture)
    kernel! = calc_Gcˣ_kernel!(device(arch), (16,16), (g.Nx, g.Ny, g.Nz))
    barrier = Event(device(arch))
    events = []
    for name in nut_names
        event = kernel!(Gcs[name].data, nut[name].data, u, g, ΔT, dependencies=barrier)
        push!(events,event)
    end
    wait(device(arch), MultiEvent(Tuple(events)))
    return nothing
end

##### calculate tendencies in y direction
@kernel function calc_Gcʸ_kernel!(Gc, c, v, g::AbstractGrid, ΔT)
    i, j, k = @index(Global, NTuple)
    ### offset index for halo points
    ii = i + g.Hx
    jj = j + g.Hy
    kk = k + g.Hz
    @inbounds Gc[ii, jj, kk] = adv_flux_y(ii, jj, kk, g, v, c, ΔT)
end
function calc_Gcsʸ!(Gcs, nut, v, g::AbstractGrid, ΔT, arch::Architecture)
    kernel! = calc_Gcʸ_kernel!(device(arch), (16,16), (g.Nx, g.Ny, g.Nz))
    barrier = Event(device(arch))
    events = []
    for name in nut_names
        event = kernel!(Gcs[name].data, nut[name].data, v, g, ΔT, dependencies=barrier)
        push!(events,event)
    end
    wait(device(arch), MultiEvent(Tuple(events)))
    return nothing
end

##### calculate tendencies in z direction
@kernel function calc_Gcᶻ_kernel!(Gc, c, w, g::AbstractGrid, ΔT)
    i, j, k = @index(Global, NTuple)
    ### offset index for halo points
    ii = i + g.Hx
    jj = j + g.Hy
    kk = k + g.Hz
    @inbounds Gc[ii, jj, kk] = adv_flux_z(ii, jj, kk, g, w, c, ΔT)
end
function calc_Gcsᶻ!(Gcs, nut, w, g::AbstractGrid, ΔT, arch::Architecture)
    kernel! = calc_Gcᶻ_kernel!(device(arch), (16,16), (g.Nx, g.Ny, g.Nz))
    barrier = Event(device(arch))
    events = []
    for name in nut_names
        event = kernel!(Gcs[name].data, nut[name].data, w, g, ΔT, dependencies=barrier)
        push!(events,event)
    end
    wait(device(arch), MultiEvent(Tuple(events)))
    return nothing
end

##### apply the tendency in x direction to tracer c
@kernel function multi_dim_x_kernel!(ctemp, Gc, c, u, g::AbstractGrid, ΔT)
    i, j, k = @index(Global, NTuple)
    ### offset index for halo points
    ii = i + g.Hx
    jj = j + g.Hy
    kk = k + g.Hz
    @inbounds ctemp[ii, jj, kk] -= ΔT / volume(ii, jj, kk, g) * (δx⁺(ii, jj, kk, Gc, g) - c[ii, jj, kk] * δx⁺(ii, jj, kk, g, Trans_x, u))
end
function multi_dim_x!(nut_temp, Gcs, nut, u, g::AbstractGrid, ΔT, arch::Architecture)
    kernel! = multi_dim_x_kernel!(device(arch), (16,16), (g.Nx, g.Ny, g.Nz))
    barrier = Event(device(arch))
    events = []
    for name in nut_names
        event = kernel!(nut_temp[name].data, Gcs[name].data, nut[name].data, u, g, ΔT, dependencies=barrier)
        push!(events,event)
    end
    wait(device(arch), MultiEvent(Tuple(events)))
    return nothing
end

##### apply the tendency in y direction to tracer c
@kernel function multi_dim_y_kernel!(ctemp, Gc, c, v, g::AbstractGrid, ΔT)
    i, j, k = @index(Global, NTuple)
    ### offset index for halo points
    ii = i + g.Hx
    jj = j + g.Hy
    kk = k + g.Hz
    @inbounds ctemp[ii, jj, kk] -= ΔT / volume(ii, jj, kk, g) * (δy⁺(ii, jj, kk, Gc, g) - c[ii, jj, kk] * δy⁺(ii, jj, kk, g, Trans_y, v ))
end
function multi_dim_y!(nut_temp, Gcs, nut, v, g::AbstractGrid, ΔT, arch::Architecture)
    kernel! = multi_dim_y_kernel!(device(arch), (16,16), (g.Nx, g.Ny, g.Nz))
    barrier = Event(device(arch))
    events = []
    for name in nut_names
        event = kernel!(nut_temp[name].data, Gcs[name].data, nut[name].data, v, g, ΔT, dependencies=barrier)
        push!(events,event)
    end
    wait(device(arch), MultiEvent(Tuple(events)))
    return nothing
end

##### apply the tendency in z direction to tracer c
@kernel function multi_dim_z_kernel!(ctemp, Gc, c, w, g::AbstractGrid, ΔT)
    i, j, k = @index(Global, NTuple)
    ### offset index for halo points
    ii = i + g.Hx
    jj = j + g.Hy
    kk = k + g.Hz
    @inbounds ctemp[ii, jj, kk] -= ΔT / volume(ii, jj, kk, g) * (δz⁺(ii, jj, kk, Gc, g) - c[ii, jj, kk] * δz⁺(ii, jj, kk, g, Trans_z, w ))
end
function multi_dim_z!(nut_temp, Gcs, nut, w, g::AbstractGrid, ΔT, arch::Architecture)
    kernel! = multi_dim_z_kernel!(device(arch), (16,16), (g.Nx, g.Ny, g.Nz))
    barrier = Event(device(arch))
    events = []
    for name in nut_names
        event = kernel!(nut_temp[name].data, Gcs[name].data, nut[name].data, w, g, ΔT, dependencies=barrier)
        push!(events,event)
    end
    wait(device(arch), MultiEvent(Tuple(events)))
    return nothing
end

function calc_nut_tendency!(a, b, c, ΔT)
    for name in nut_names
        @inbounds a[name].data .= b[name].data .- c[name].data
    end
end

function nut_advection!(nut, nut_temp, Gcs, vel, g::AbstractGrid, ΔT, arch::Architecture)
    for name in nut_names
        @inbounds nut_temp[name].data .= nut[name].data
    end
    ##### x direction
    calc_Gcsˣ!(Gcs, nut_temp, vel.u.data, g, ΔT, arch)
    fill_halo_Gcs!(Gcs, g)
    multi_dim_x!(nut_temp, Gcs, nut, vel.u.data, g, ΔT, arch)
    fill_halo_nut!(nut_temp, g)

    ##### y direction
    calc_Gcsʸ!(Gcs, nut_temp, vel.v.data, g, ΔT, arch)
    fill_halo_Gcs!(Gcs, g)
    multi_dim_y!(nut_temp, Gcs, nut, vel.v.data, g, ΔT, arch)
    fill_halo_nut!(nut_temp, g)

    ##### z direction
    calc_Gcsᶻ!(Gcs, nut_temp, vel.w.data, g, ΔT, arch)
    fill_halo_Gcs!(Gcs, g)
    multi_dim_z!(nut_temp, Gcs, nut, vel.w.data, g, ΔT, arch)
    fill_halo_nut!(nut_temp, g)

    calc_nut_tendency!(Gcs, nut_temp, nut, ΔT)

end


text:
# Autogenerated wrapper script for Gama_jll for i686-w64-mingw32-cxx11
export gama_g3, gama_local, gama_local_gkf2yaml, gama_local_xml2sql, gama_local_xml2txt, gama_local_yaml2gkf

using Expat_jll
using yaml_cpp_jll
using SQLite_jll
JLLWrappers.@generate_wrapper_header("Gama")
JLLWrappers.@declare_executable_product(gama_g3)
JLLWrappers.@declare_executable_product(gama_local)
JLLWrappers.@declare_executable_product(gama_local_gkf2yaml)
JLLWrappers.@declare_executable_product(gama_local_xml2sql)
JLLWrappers.@declare_executable_product(gama_local_xml2txt)
JLLWrappers.@declare_executable_product(gama_local_yaml2gkf)
function __init__()
    JLLWrappers.@generate_init_header(Expat_jll, yaml_cpp_jll, SQLite_jll)
    JLLWrappers.@init_executable_product(
        gama_g3,
        "bin\\gama-g3.exe",
    )

    JLLWrappers.@init_executable_product(
        gama_local,
        "bin\\gama-local.exe",
    )

    JLLWrappers.@init_executable_product(
        gama_local_gkf2yaml,
        "bin\\gama-local-gkf2yaml.exe",
    )

    JLLWrappers.@init_executable_product(
        gama_local_xml2sql,
        "bin\\gama-local-xml2sql.exe",
    )

    JLLWrappers.@init_executable_product(
        gama_local_xml2txt,
        "bin\\gama-local-xml2txt.exe",
    )

    JLLWrappers.@init_executable_product(
        gama_local_yaml2gkf,
        "bin\\gama-local-yaml2gkf.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function IodGibbs(arg0::jdouble)
    return IodGibbs((jdouble,), arg0)
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function estimate(obj::IodGibbs, arg0::Frame, arg1::PV, arg2::PV, arg3::PV)
    return jcall(obj, "estimate", KeplerianOrbit, (Frame, PV, PV, PV), arg0, arg1, arg2, arg3)
end

function estimate(obj::IodGibbs, arg0::Frame, arg1::Vector3D, arg2::AbsoluteDate, arg3::Vector3D, arg4::AbsoluteDate, arg5::Vector3D, arg6::AbsoluteDate)
    return jcall(obj, "estimate", KeplerianOrbit, (Frame, Vector3D, AbsoluteDate, Vector3D, AbsoluteDate, Vector3D, AbsoluteDate), arg0, arg1, arg2, arg3, arg4, arg5, arg6)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
# page constructor

#

max_stars_repo_name: Clpr/HtmlConstructor

max_issues_repo_name: Clpr/HtmlConstructor

max_forks_repo_name: Clpr/HtmlConstructor


text:
using NISTStatisticalReferenceData
using Test

@testset "NISTStatisticalReferenceData.jl" begin
    @testset "Load Every Dataset" begin
        dsets = NISTStatisticalReferenceData.datasets()
        for (k,v) in dsets
            for vi in v
                println("Loading $vi...")
                t = NISTStatisticalReferenceData.TestCase(vi)
                # check that at least one certified value was found in every dat file.
                @test !isempty(t.certified_values)
            end
        end
    end
end


text:
import DataFrames
import StatsBase

function get_countmap(itr; skip_missings::Bool = true)
    if skip_missings
        result = get_countmap_skip_missings(itr)
    else
        result = get_countmap_include_missings(itr)
    end
    return result
end

function get_countmap_skip_missings(itr)
    result = StatsBase.countmap(collect(skipmissing(itr)))
    return result
end

function get_countmap_include_missings(itr)
    result = StatsBase.countmap(collect(itr))
    return result
end

function get_unique_values(itr; skip_missings::Bool = true)
    if skip_missings
        result = get_unique_values_skip_missings(itr)
    else
        result = get_unique_values_include_missings(itr)
    end
    return result
end

function get_unique_values_skip_missings(itr)
    result = keys(get_countmap_skip_missings(itr))
    return result
end

function get_unique_values_include_missings(itr)
    result = keys(get_countmap_include_missings(itr))
    return result
end

function get_number_of_unique_values(itr; skip_missings::Bool = true)::Int
    if skip_missings
        result = get_number_of_unique_values_skip_missings(itr)
    else
        result = get_number_of_unique_values_include_missings(itr)
    end
    return result
end

function get_number_of_unique_values_skip_missings(itr)::Int
    result::Int = length(get_unique_values_skip_missings(itr))
    return result
end

function get_number_of_unique_values_include_missings(itr)::Int
    result::Int = length(get_unique_values_include_missings(itr))
    return result
end

function find_constant_columns(
        df::DataFrames.AbstractDataFrame,
        )::Vector{Symbol}
    list_of_constant_column_names = Symbol[]
    for x in DataFrames.names(df)
        if get_number_of_unique_values(df[x]) < 2
            push!(list_of_constant_column_names, x)
        end
    end
    return list_of_constant_column_names
end

function check_no_constant_columns(
        df::DataFrames.AbstractDataFrame,
        )::Bool
    list_of_constant_column_names = find_constant_columns(df)
    if length(list_of_constant_column_names) > 0
        error(
            string(
                "Data frame contains the following constant columns: ",
                list_of_constant_column_names,
                "\"",
                )
            )
    else
        return true
    end
end


text:
__precompile__()

module BinDeps

using Compat
using Compat: String, KERNEL

export @make_run, @build_steps, find_library, download_cmd, unpack_cmd,
        Choice, Choices, CCompile, FileDownloader, FileRule,
        ChangeDirectory, FileDownloader, FileUnpacker, prepare_src,
        autotools_install, CreateDirectory, MakeTargets, SystemLibInstall

const dlext = Libdl.dlext
const shlib_ext = dlext # compatibility with older packages (e.g. ZMQ)

function find_library(pkg,libname,files)
    Base.warn_once("BinDeps.find_library is deprecated, use Base.find_library instead.")
    dl = C_NULL
    for filename in files
        dl = Libdl.dlopen_e(joinpath(Pkg.dir(),pkg,"deps","usr","lib",filename))
        if dl != C_NULL
            ccall(:add_library_mapping,Cint,(Ptr{Cchar},Ptr{Void}),libname,dl)
            return true
        end

        dl = Libdl.dlopen_e(filename)
        if dl != C_NULL
            ccall(:add_library_mapping,Cint,(Ptr{Cchar},Ptr{Void}),libname,dl)
            return true
        end
    end

    dl = Libdl.dlopen_e(libname)
    dl != C_NULL ? true : false
end

macro make_rule(condition,command)
    quote
        if(!$(esc(condition)))
            $(esc(command))
            @assert $(esc(condition))
        end
    end
end

@compat abstract type BuildStep end

downloadcmd = nothing
function download_cmd(url::AbstractString, filename::AbstractString)
    global downloadcmd
    if downloadcmd === nothing
        for download_engine in (is_windows() ? ("C:\\Windows\\System32\\WindowsPowerShell\\v1.0\\powershell",
                :powershell, :curl, :wget, :fetch) : (:curl, :wget, :fetch))
            if endswith(string(download_engine), "powershell")
                checkcmd = `$download_engine -NoProfile -Command ""`
            else
                checkcmd = `$download_engine --help`
            end
            try
                if success(checkcmd)
                    downloadcmd = download_engine
                    break
                end
            catch
                continue # don't bail if one of these fails
            end
        end
    end
    if downloadcmd == :wget
        return `$downloadcmd -O $filename $url`
    elseif downloadcmd == :curl
        return `$downloadcmd -f -o $filename -L $url`
    elseif downloadcmd == :fetch
        return `$downloadcmd -f $filename $url`
    elseif endswith(string(downloadcmd), "powershell")
        return `$downloadcmd -NoProfile -Command "(new-object net.webclient).DownloadFile(\"$url\", \"$filename\")"`
    else
        extraerr = is_windows() ? "check if powershell is on your path or " : ""
        error("No download agent available; $(extraerr)install curl, wget, or fetch.")
    end
end

if is_unix()
    function unpack_cmd(file,directory,extension,secondary_extension)
        if ((extension == ".gz" || extension == ".Z") && secondary_extension == ".tar") || extension == ".tgz"
            return (`tar xzf $file --directory=$directory`)
        elseif (extension == ".bz2" && secondary_extension == ".tar") || extension == ".tbz"
            return (`tar xjf $file --directory=$directory`)
        elseif extension == ".xz" && secondary_extension == ".tar"
            return pipeline(`unxz -c $file `, `tar xv --directory=$directory`)
        elseif extension == ".tar"
            return (`tar xf $file --directory=$directory`)
        elseif extension == ".zip"
            return (`unzip -x $file -d $directory`)
        elseif extension == ".gz"
            return pipeline(`mkdir $directory`, `cp $file $directory`, `gzip -d $directory/$file`)
        end
        error("I don't know how to unpack $file")
    end
end

if is_windows()
    function unpack_cmd(file,directory,extension,secondary_extension)
        if ((extension == ".Z" || extension == ".gz" || extension == ".xz" || extension == ".bz2") &&
                secondary_extension == ".tar") || extension == ".tgz" || extension == ".tbz"
            return pipeline(`7z x $file -y -so`, `7z x -si -y -ttar -o$directory`)
        elseif (extension == ".zip" || extension == ".7z" || extension == ".tar" ||
                (extension == ".exe" && secondary_extension == ".7z"))
            return (`7z x $file -y -o$directory`)
        end
        error("I don't know how to unpack $file")
    end
end

type SynchronousStepCollection
    steps::Vector{Any}
    cwd::AbstractString
    oldcwd::AbstractString
    SynchronousStepCollection(cwd) = new(Any[],cwd,cwd)
    SynchronousStepCollection() = new(Any[],"","")
end

import Base.push!, Base.run, Base.(|)
push!(a::SynchronousStepCollection,args...) = push!(a.steps,args...)

type ChangeDirectory <: BuildStep
    dir::AbstractString
end

type CreateDirectory <: BuildStep
    dest::AbstractString
    mayexist::Bool
    CreateDirectory(dest, me) = new(dest,me)
    CreateDirectory(dest) = new(dest,true)
end

immutable RemoveDirectory <: BuildStep
    dest::AbstractString
end

type FileDownloader <: BuildStep
    src::AbstractString     #url
    dest::AbstractString    #local_file
end

type ChecksumValidator <: BuildStep
    sha::AbstractString
    path::AbstractString
end

type FileUnpacker <: BuildStep
    src::AbstractString     #file
    dest::AbstractString    #directory
    target::AbstractString  #file or directory inside the archive to test
                    #for existence (or blank to check for a.tgz => a/)
end

type MakeTargets <: BuildStep
    dir::AbstractString
    targets::Vector{String}
    env::Dict
    MakeTargets(dir,target;env = Dict{AbstractString,AbstractString}()) = new(dir,target,env)
    MakeTargets{S<:AbstractString}(target::Vector{S};env = Dict{AbstractString,AbstractString}()) = new("",target,env)
    MakeTargets(target::String;env = Dict{AbstractString,AbstractString}()) = new("",[target],env)
    MakeTargets(;env = Dict{AbstractString,AbstractString}()) = new("",String[],env)
end

type AutotoolsDependency <: BuildStep
    src::AbstractString     #src direcory
    prefix::AbstractString
    builddir::AbstractString
    configure_options::Vector{AbstractString}
    libtarget::Vector{AbstractString}
    include_dirs::Vector{AbstractString}
    lib_dirs::Vector{AbstractString}
    rpath_dirs::Vector{AbstractString}
    installed_libpath::Vector{String} # The library is considered installed if any of these paths exist
    config_status_dir::AbstractString
    force_rebuild::Bool
    env
    AutotoolsDependency(;srcdir::AbstractString = "", prefix = "", builddir = "", configure_options=AbstractString[], libtarget = AbstractString[], include_dirs=AbstractString[], lib_dirs=AbstractString[], rpath_dirs=AbstractString[], installed_libpath = String[], force_rebuild=false, config_status_dir = "", env = Dict{String,String}()) =
        new(srcdir,prefix,builddir,configure_options,isa(libtarget,Vector)?libtarget:AbstractString[libtarget],include_dirs,lib_dirs,rpath_dirs,installed_libpath,config_status_dir,force_rebuild,env)
end

### Choices

type Choice
    name::Symbol
    description::AbstractString
    step::SynchronousStepCollection
    Choice(name,description,step) = (s=SynchronousStepCollection();lower(step,s);new(name,description,s))
end

type Choices <: BuildStep
    choices::Vector{Choice}
    Choices() = new(Choice[])
    Choices(choices::Vector{Choice}) = new(choices)
end

push!(c::Choices, args...) = push!(c.choices, args...)

function run(c::Choices)
    println()
    info("There are multiple options available for installing this dependency:")
    while true
        for x in c.choices
            println("- "*string(x.name)*": "*x.description)
        end
        while true
            print("Plese select desired method: ")
            method = Symbol(chomp(readline(STDIN)))
            for x in c.choices
                if(method == x.name)
                    return run(x.step)
                end
            end
            warn("Invalid Method")
        end
    end
end

type CCompile <: BuildStep
    srcFile::AbstractString
    destFile::AbstractString
    options::Vector{String}
    libs::Vector{String}
end

lower(cc::CCompile,c) = lower(FileRule(cc.destFile,`gcc $(cc.options) $(cc.srcFile) $(cc.libs) -o $(cc.destFile)`),c)
##

type DirectoryRule <: BuildStep
    dir::AbstractString
    step
end

type PathRule <: BuildStep
    path::AbstractString
    step
end

function meta_lower(a::Expr,blk::Expr,collection)
    if(a.head == :block || a.head == :tuple)
        for x in a.args
            if(isa(x,Expr))
                if(x.head == :block)
                    new_collection = gensym()
                    push!(blk.args,quote
                        $(esc(new_collection)) = SynchronousStepCollection($(esc(collection)).cwd)
                        push!($(esc(collection)),$(esc(new_collection)))
                    end)
                    meta_lower(x,blk,new_collection)
                 elseif(x.head != :line)
                     push!(blk.args,quote
                         lower($(esc(x)), $(esc(collection)))
                     end)
                 end
            elseif(!isa(x,LineNumberNode))
                meta_lower(x,blk,collection)
            end
        end
    else
        push!(blk.args,quote
            $(esc(collection)),lower($(esc(a)), $(esc(collection)))
        end)
    end
end

function meta_lower(a::Tuple,blk::Expr,collection)
    for x in a
        meta_lower(a,blk,collection)
    end
end

function meta_lower(a,blk::Expr,collection)
    push!(blk.args,quote
        $(esc(collection)), lower($(esc(a)), $(esc(collection)))
    end)
end

macro dependent_steps(steps)
    blk = Expr(:block)
    meta_lower(steps,blk,:collection)
    blk
end

mypwd() = chomp(readall(`pwd`))

macro build_steps(steps)
    collection = gensym()
    blk = Expr(:block)
    push!(blk.args,quote
        $(esc(collection)) = SynchronousStepCollection()
    end)
    meta_lower(steps,blk,collection)
    push!(blk.args, quote; $(esc(collection)); end)
    blk
end

src(b::BuildStep) = b.src
dest(b::BuildStep) = b.dest

(|)(a::BuildStep,b::BuildStep) = SynchronousStepCollection()
function (|)(a::SynchronousStepCollection,b::SynchronousStepCollection)
	if a.cwd == b.cwd
		append!(a.steps,b.steps)
	else
		push!(a.steps,b)
	end
	a
end
(|)(a::SynchronousStepCollection,b::Function) = (lower(b,a);a)
(|)(a::SynchronousStepCollection,b) = (lower(b,a);a)

(|)(b::Function,a::SynchronousStepCollection) = (c=SynchronousStepCollection(); ((c|b)|a))
(|)(b,a::SynchronousStepCollection) = (c=SynchronousStepCollection(); ((c|b)|a))

# Create any of these files
type FileRule <: BuildStep
    file::Array{AbstractString}
    step
    FileRule(file::AbstractString,step) = FileRule(AbstractString[file],step)
    function FileRule(files::Vector{AbstractString},step)
        new(files,@build_steps (step,) )
    end
end
FileRule{T<:AbstractString}(files::Vector{T},step) = FileRule(AbstractString[f for f in files],step)

function lower(s::ChangeDirectory,collection)
    if !isempty(collection.steps)
        error("Change of Directory must be the first instruction")
    end
    collection.cwd = s.dir
end
lower(s::Void,collection) = nothing
lower(s::Function,collection) = push!(collection,s)
lower(s::CreateDirectory,collection) = @dependent_steps ( DirectoryRule(s.dest,()->(mkpath(s.dest))), )
lower(s::RemoveDirectory,collection) = @dependent_steps ( `rm -rf $(s.dest)` )
lower(s::BuildStep,collection) = push!(collection,s)
lower(s::Base.AbstractCmd,collection) = push!(collection,s)
lower(s::FileDownloader,collection) = @dependent_steps ( CreateDirectory(dirname(s.dest),true), ()->info("Downloading file $(s.src)"), FileRule(s.dest,download_cmd(s.src,s.dest)), ()->info("Done downloading file $(s.src)") )
lower(s::ChecksumValidator,collection) = isempty(s.sha) || @dependent_steps ()->sha_check(s.path, s.sha)
function splittarpath(path)
    path,extension = splitext(path)
    base_filename,secondary_extension = splitext(path)
    if extension == ".tgz" || extension == ".tbz" || extension == ".zip" && !isempty(secondary_extension)
        base_filename *= secondary_extension
        secondary_extension = ""
    end
    (base_filename,extension,secondary_extension)
end
function lower(s::FileUnpacker,collection)
    base_filename,extension,secondary_extension = splittarpath(s.src)
    target = !isempty(s.target) ? s.target : basename(base_filename)
    @dependent_steps begin
        CreateDirectory(dirname(s.dest),true)
        PathRule(joinpath(s.dest,target),unpack_cmd(s.src,s.dest,extension,secondary_extension))
    end
end

function adjust_env(env)
    ret = similar(env)
    merge!(ret,ENV)
    merge!(ret,env) #s.env overrides ENV
    ret
end

if is_unix()
    function lower(a::MakeTargets,collection)
        cmd = `make -j8`
        if !isempty(a.dir)
            cmd = `$cmd -C $(a.dir)`
        end
        if !isempty(a.targets)
            cmd = `$cmd $(a.targets)`
        end
        @dependent_steps ( setenv(cmd, adjust_env(a.env)), )
    end
end
is_windows() && (lower(a::MakeTargets,collection) = @dependent_steps ( setenv(`make $(!isempty(a.dir)?"-C "*a.dir:"") $(a.targets)`, adjust_env(a.env)), ))
lower(s::SynchronousStepCollection,collection) = (collection|=s)

lower(s) = (c=SynchronousStepCollection();lower(s,c);c)

#run(s::MakeTargets) = run(@make_steps (s,))

function lower(s::AutotoolsDependency,collection)
    prefix = s.prefix
    if is_windows()
        prefix = replace(replace(s.prefix,"\\","/"),"C:/","/c/")
    end
	cmdstring = "pwd && ./configure --prefix=$(prefix) "*join(s.configure_options," ")

    env = adjust_env(s.env)

    for path in s.include_dirs
        if !haskey(env,"CPPFLAGS")
            env["CPPFLAGS"] = ""
        end
        env["CPPFLAGS"]*=" -I$path"
    end

    for path in s.lib_dirs
        if !haskey(env,"LDFLAGS")
            env["LDFLAGS"] = ""
        end
        env["LDFLAGS"]*=" -L$path"
    end

    for path in s.rpath_dirs
        if !haskey(env,"LDFLAGS")
            env["LDFLAGS"] = ""
        end
        env["LDFLAGS"]*=" -Wl,-rpath -Wl,$path"
    end

    if s.force_rebuild
        @dependent_steps begin
            RemoveDirectory(s.builddir)
        end
    end

    @static if is_unix()
        @dependent_steps begin
            CreateDirectory(s.builddir)
            begin
                ChangeDirectory(s.builddir)
                FileRule(isempty(s.config_status_dir)?"config.status":joinpath(s.config_status_dir,"config.status"), setenv(`$(s.src)/configure $(s.configure_options) --prefix=$(prefix)`,env))
                FileRule(s.libtarget,MakeTargets(;env=s.env))
                MakeTargets("install";env=env)
            end
        end
    end

    @static if is_windows()
        @dependent_steps begin
            ChangeDirectory(s.src)
            FileRule(isempty(s.config_status_dir)?"config.status":joinpath(s.config_status_dir,"config.status"),setenv(`sh -c $cmdstring`,env))
            FileRule(s.libtarget,MakeTargets())
            MakeTargets("install")
        end
    end
end

function run(f::Function)
	f()
end

function run(s::FileRule)
    if !any(map(isfile,s.file))
        run(s.step)
		if !any(map(isfile,s.file))
			error("File $(s.file) was not created successfully (Tried to run $(s.step) )")
		end
    end
end
function run(s::DirectoryRule)
	info("Attempting to Create directory $(s.dir)")
    if !isdir(s.dir)
        run(s.step)
		if !isdir(s.dir)
			error("Directory $(s.dir) was not created successfully (Tried to run $(s.step) )")
		end
	else
		info("Directory $(s.dir) already created")
    end
end

function run(s::PathRule)
    if !ispath(s.path)
        run(s.step)
        if !ispath(s.path)
            error("Path $(s.path) was not created successfully (Tried to run $(s.step) )")
        end
    else
        info("Path $(s.path) already created")
    end
end

function run(s::BuildStep)
    error("Unimplemented BuildStep: $(typeof(s))")
end
function run(s::SynchronousStepCollection)
    for x in s.steps
		if !isempty(s.cwd)
			info("Changing Directory to $(s.cwd)")
			cd(s.cwd)
		end
        run(x)
        if !isempty(s.oldcwd)
			info("Changing Directory to $(s.oldcwd)")
			cd(s.oldcwd)
		end
    end
end

is_unix() && (make_command = `make -j8`)
is_windows() && (make_command = `make`)

function prepare_src(depsdir,url, downloaded_file, directory_name)
    local_file = joinpath(joinpath(depsdir,"downloads"),downloaded_file)
	@build_steps begin
        FileDownloader(url,local_file)
        FileUnpacker(local_file,joinpath(depsdir,"src"),directory_name)
	end
end

function autotools_install(depsdir,url, downloaded_file, configure_opts, directory_name, directory, libname, installed_libname, confstatusdir)
    prefix = joinpath(depsdir,"usr")
    libdir = joinpath(prefix,"lib")
    srcdir = joinpath(depsdir,"src",directory)
    dir = joinpath(joinpath(depsdir,"builds"),directory)
    prepare_src(depsdir,url, downloaded_file,directory_name) |
	@build_steps begin
        AutotoolsDependency(srcdir=srcdir,prefix=prefix,builddir=dir,configure_options=configure_opts,libtarget=libname,installed_libpath=[joinpath(libdir,installed_libname)],config_status_dir=confstatusdir)
    end
end
autotools_install(depsdir,url, downloaded_file, configure_opts, directory_name, directory, libname, installed_libname) = autotools_install(depsdir,url, downloaded_file, configure_opts, directory_name, directory, libname, installed_libname, "")
autotools_install(depsdir,url, downloaded_file, configure_opts, directory, libname)=autotools_install(depsdir,url,downloaded_file,configure_opts,directory,directory,libname,libname)

autotools_install(args...) = error("autotools_install has been removed")

include("dependencies.jl")
include("debug.jl")
include("show.jl")
end


text:
# This file includes investigates the 2D interpolation mean square error(MSE) with respect to inteprolation point number.

using FractalTools
using Makie

# Construct interpolation data
f(x, y) = x^2 + y^2 + 1
vtx = [
    BigFloat.([0.0, 0.0]),
    BigFloat.([1.0, 0.0]),
    BigFloat.([0.5, 1.0])
    ]
freevar = 0.001
npts    = 50 : 5 : 150
ntpts   = 2 * npts[end]

# Construct test data
tpts = getdata(vtx, ntpts)

# Compute errors
mse = map(npts) do npt
    @info npt
    # Construct interpolation data
    pts = getdata(f, vtx, npt)

    # Construct interpolant
    interp = interpolate(pts, Interp2D(freevar))

    # Compute error
    fvals = map(pnt -> f(pnt...), tpts)
    ivals = map(pnt -> interp(pnt...), tpts)
    sum((fvals - ivals).^2) / length(tpts)
end

# Plot mse
fig = Figure()
ax = fig[1, 1] = Axis(fig, xlabel="Number of Points", ylabel="MSE", title="2D Interpolation MSE")
stem!(ax, npts, mse, color=:black)
save(joinpath(@__DIR__, "interp2d_mse.png"), fig)
display(fig)


text:
using Pkg
Pkg.activate(@__DIR__)
cd(@__DIR__)
Pkg.update()

using Documenter, Agents, DataFrames, Random, Statistics, SQLite
using AgentsPlots
using Literate
using Plots, StatsPlots
const CI = get(ENV, "CI", nothing) == "true"
CI && (ENV["GKSwstype"] = "100")

ENV["GKS_ENCODING"]="utf-8"
println("Agents status:")
Pkg.status("Agents")

println("AgentsPlots status:")
Pkg.status("AgentsPlots")

# %% Literate convertion
indir = joinpath(@__DIR__, "..", "examples")
outdir = joinpath(@__DIR__, "src", "examples")
mkpath(outdir)
toskip = ()
for file in readdir(indir)
    file ∈ toskip && continue
    Literate.markdown(joinpath(indir, file), outdir; credit = false)
end

# %%
# download the themes
using DocumenterTools: Themes
for file in ("juliadynamics-lightdefs.scss", "juliadynamics-darkdefs.scss", "juliadynamics-style.scss")
    download("https://raw.githubusercontent.com/JuliaDynamics/doctheme/master/$file", joinpath(@__DIR__, file))
end
# create the themes
for w in ("light", "dark")
    header = read(joinpath(@__DIR__, "juliadynamics-style.scss"), String)
    theme = read(joinpath(@__DIR__, "juliadynamics-$(w)defs.scss"), String)
    write(joinpath(@__DIR__, "juliadynamics-$(w).scss"), header*"\n"*theme)
end
# compile the themes
Themes.compile(joinpath(@__DIR__, "juliadynamics-light.scss"), joinpath(@__DIR__, "src/assets/themes/documenter-light.css"))
Themes.compile(joinpath(@__DIR__, "juliadynamics-dark.scss"), joinpath(@__DIR__, "src/assets/themes/documenter-dark.css"))

# %%
cd(@__DIR__)
ENV["JULIA_DEBUG"] = "Documenter"
makedocs(modules = [Agents,AgentsPlots],
sitename= "Agents.jl",
authors = "Ali R. Vahdati, George Datseris, Tim DuBois and contributors.",
doctest = false,
format = Documenter.HTML(
    prettyurls = CI,
    assets = [
        asset("https://fonts.googleapis.com/css?family=Montserrat|Source+Code+Pro&display=swap", class=:css),
    ],
    ),
pages = [
    "Introduction" => "index.md",
    "Tutorial" => "tutorial.md",
    "Examples" => [
        "Schelling's segregation model" => "examples/schelling.md",
        "SIR model for the spread of COVID-19" => "examples/sir.md",
        "Continuous space social distancing for COVID-19" => "examples/social_distancing.md",
        "Wealth distribution" => "examples/wealth_distribution.md",
        "Forest fire" => "examples/forest_fire.md",
        "Conway's game of life" => "examples/game_of_life_2D_CA.md",
        "Wright-Fisher model of evolution" => "examples/wright-fisher.md",
        "Hegselmann-Krause opinion dynamics" => "examples/HK.md",
        "Flocking" => "examples/flock.md",
        "Daisyworld" => "examples/daisyworld.md",
        "Predator-Prey" => "examples/predator_prey.md",
        "Bacteria Growth" => "examples/growing_bacteria.md"
        ],
    "API" => "api.md",
    "Comparison against Mesa (Python)" => "mesa.md"
    ],
)

if CI
    deploydocs(
        repo = "github.com/JuliaDynamics/Agents.jl.git",
        target = "build",
        push_preview = true
    )
end

println("done")

max_stars_repo_name: werzum/AgentsJLFork

max_issues_repo_name: werzum/AgentsJLFork

max_forks_repo_name: werzum/AgentsJLFork


text:
# setup project paths -
const _PATH_TO_SRC = dirname(pathof(@__MODULE__))
const _PATH_TO_BASE = joinpath(_PATH_TO_SRC, "base")
const _PATH_TO_EXPIRATION = joinpath(_PATH_TO_SRC, "expiration")
const _PATH_TO_LATTICE = joinpath(_PATH_TO_SRC, "lattice")
const _PATH_TO_MONTECARLO = joinpath(_PATH_TO_SRC, "montecarlo")

# import packages that we will use -
using JSON
using Dates
using DataFrames
using CSV
using Plots

# include our codes -
include(joinpath(_PATH_TO_BASE, "VLTypes.jl"))
include(joinpath(_PATH_TO_BASE, "VLBase.jl"))
include(joinpath(_PATH_TO_BASE, "VLFactory.jl"))
include(joinpath(_PATH_TO_BASE, "VLIntrinsic.jl"))
include(joinpath(_PATH_TO_EXPIRATION, "VLExpiration.jl"))
include(joinpath(_PATH_TO_LATTICE, "VLBinomialLattice.jl"))
include(joinpath(_PATH_TO_LATTICE, "VLLatticeGreeks.jl"))
include(joinpath(_PATH_TO_MONTECARLO, "VLLongstaffMC.jl"))


text:
using CudaPOMDPs
using Test

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


text:
const c=299792458.

function IC(l,p,h,X,Y,Z,tilt,azimut,phase,amplitude,order)
  const Rx=vcat([1 0 0],[0 -1 0],[0 0 -1])
  const Ry=vcat([-1 0 0],[0 1 0],[0 0 -1])
  const Rz=vcat([-1 0 0],[0 -1 0],[0 0 1])
  const i0=amplitude*[sin(tilt).*cos(azimut), sin(tilt).*sin(azimut), cos(tilt)]
  POS=zeros((order+1)^3*8,7)

  count=1
  for i=0:order
    perc=round(i/order*1000)/10
    println("$perc %")
    for j=0:order
      for k=0:order
        n = int(i+j+k);
        ip = (Rx^i * Ry^j * Rz^k) * i0;
        x = (i+mod(i,2))*l+(-1)^i*X;
        y = (j+mod(j,2))*p+(-1)^j*Y;
        z = (k+mod(k,2))*h+(-1)^k*Z;
        POS[count,:]=[ip' x  y  z n ]
        POS[count+1,:]=[(Rz*ip)' x  y  -z n+1 ]
        POS[count+2,:]=[(Ry*ip)' x  -y  z n+1 ]
        POS[count+3,:]=[(Ry*Rz*ip)' x  -y  -z n+2 ]
        POS[count+4,:]=[(Rx*ip)' -x  y  z n+1 ]
        POS[count+5,:]=[(Rx*Rz*ip)' -x  y  -z n+2 ]
        POS[count+6,:]=[(Rx*Ry*ip)' -x  -y  z n+2 ]
        POS[count+7,:]=[(Rx*Ry*Rz*ip)' -x  -y  -z n+3 ]
        count+=8
      end
    end
  end
  perm=sortperm(POS[:,7])
  return POS[perm,:]
end

max_stars_repo_name: manuamador/Z

max_issues_repo_name: manuamador/Z

max_forks_repo_name: manuamador/Z


text:
using Test
using DirectedAcyclicGraphs

module TestNodes

    using Test
    using ..DirectedAcyclicGraphs

    abstract type TestNode <: Tree end
    mutable struct TestINode <: TestNode
        id::Int
        children::Vector{Tree}
        parent
        TestINode(i,c) = begin
            x = new(i,c, nothing)
            for n in c
                n.parent = x
            end
            x
        end
    end

    mutable struct TestLNode <: TestNode
        id::Int
        parent
        TestLNode(i) = new(i, nothing)
    end

    DirectedAcyclicGraphs.NodeType(::Type{<:TestINode}) = Inner()
    DirectedAcyclicGraphs.NodeType(::Type{<:TestLNode}) = Leaf()
    DirectedAcyclicGraphs.children(n::TestINode) = n.children

    @testset "Tree utils for TestNodes" begin

        l1 = TestLNode(1)
        l2 = TestLNode(2)

        @test !has_children(l1)
        @test num_children(l1) == 0
        @test isleaf(l1)
        @test !isinner(l1)

        i12 = TestINode(3,[l1,l2])

        @test has_children(i12)

        @test num_children(i12) == 2

        l3 = TestLNode(4)
        l4 = TestLNode(5)

        i34 = TestINode(6,[l3,l4])

        r = TestINode(7,[i12,i34])

        @test has_children(r)
        @test num_children(r) == 2

        @test has_parent(l1)
        @test has_parent(i12)
        @test !has_parent(r)
        @test parent(l1) === i12
        @test parent(i12) === r

        @test isroot(r)
        @test !isroot(i12)
        @test root(l4) === r

        foreach(r) do n
            n.id += 1
        end
        @test l1.id == 2
        @test l2.id == 3
        @test i12.id == 4
        @test l3.id == 5
        @test l4.id == 6
        @test i34.id == 7
        @test r.id == 8

        foreach(r, l -> l.id += 1, i -> i.id -= 1)

        @test l1.id == 2 + 1
        @test l2.id == 3 + 1
        @test i12.id == 4 - 1
        @test l3.id == 5 + 1
        @test l4.id == 6 + 1
        @test i34.id == 7 - 1
        @test r.id == 8 - 1

        @test filter(n -> iseven(n.id), r) == [l2,l3,i34]

        lastvisited = nothing
        foreach(n -> lastvisited=n,r)
        @test lastvisited === r

        lastvisited = nothing
        foreach_down(n -> lastvisited=n, r)
        @assert isleaf(lastvisited) "$lastvisited"
        @test isleaf(lastvisited)

        @test num_nodes(r) == 7
        @test num_edges(r) == 6

        @test isempty(innernodes(l1))
        @test leafnodes(l1) == [l1]

        @test issetequal(innernodes(r), [i12,i34,r])
        @test issetequal(leafnodes(r), [l1,l2,l3,l4])

        @test tree_num_edges(r) == 6

        @test linearize(r)[end] == r
        @test isleaf(linearize(r)[1])
        @test linearize(l2) == [l2]
        @test length(linearize(i12)) == 3

        @test eltype(linearize(r)) == TestNode
        @test eltype(linearize(l1)) == TestLNode
        @test eltype(linearize(r, Any)) == Any

        @test left_most_descendent(r) == l1
        @test right_most_descendent(r) == l4

        np = num_parents(r)
        @test np[r] == 0
        @test np[l1] == 1
        @test np[i34] == 1

        f_l(n) = 1
        f_i1(n,cs) = sum(cs) + 1
        f_i2(n,call) = mapreduce(call, +, children(n)) + 1

        @test foldup(r, f_l, f_i2, Int) == foldup_aggregate(r, f_l, f_i1, Int)

        function df(n,m)
            v = if m === i12
                n === l1 || n === l2
            elseif m === i34
                n === l3 || n === l4
            elseif m === r
                true
            else
                @assert isleaf(m)
                false
            end
            v || (n === m)
        end

        DirectedAcyclicGraphs.lca(x::TestNode, y::TestNode) =
            lca(x,y,df)

        @test lca(l1,l2) == i12
        @test lca(l1,l2,i34) == r
        @test lca(l1) == l1
        @test_throws ErrorException lca(l3, TestLNode(2))
        @test lca(nothing, l1, df) == l1
        @test lca(l1, nothing, df) == l1
        @test lca(nothing, nothing, df) === nothing

        @test find_inode(l1,l2,df) == i12
        @test find_inode(l3,l4,df) == i34
        @test find_inode(l1,l4,df) == r
        @test find_inode(l2,nothing,df) == r
        @test find_inode(nothing,l3,df) == r
        @test find_inode(l2,l1,df) === nothing

        @test find_leaf(r, x -> true) == l1
        @test_throws ErrorException find_leaf(r, x -> false)

        @test depth(r, x -> true) == 2
        @test_throws ErrorException depth(r, x -> false)

        b = IOBuffer()
        print_tree(r,b)
        @test !isempty(String(take!(b)))

        k1 = TestLNode(3)
        k2 = TestLNode(4)
        j12 = TestINode(3,[k1,k2])

        @test k1 != l1
        @test j12 != i12

        DirectedAcyclicGraphs.isequal_local(x::TestNode,y::TestNode) =
            x.id == y.id

        @test k1 != l2
        @test k1 == l1
        @test j12 != i34
        @test j12 == i12

    end

end


text:
########## Computational Harmonic Analysis and Wavelet transform ##########

"""
    native_conv(u::AbstractVector{<:Number}, v::AbstractVector{<:Number})

Convolution between two vectors computed in the native way.
"""
function native_conv(u::AbstractVector{<:Number}, v::AbstractVector{<:Number})
    m = length(u)
    n = length(v)
    w = zeros(promote_type(eltype(u),eltype(v)), m+n-1)

    @inbounds begin
        for j=1:m, k=1:n
            w[j+k-1] += u[j] * v[k]
        end
    end
    return w
end

# Normalizations of Wavelets.jl for different families are not coherent:
# Assuming the sqrt(2) factor in the cascade algorithm (see wavefun), then the filters of the following
# family have to be rescaled by the corresponding factor
# - Daubechies: 1
# - Coiflet: 1
# - Symlet: 1/sqrt(2)
# - Battle-Lemarie: sqrt(2)
# - Beylkin: 1
# - Vaidyanathan: 1

"""
    convolution_matrix(h::AbstractVector{<:Number}, N::Integer)

Construct the matrix representation of the full convolution between a kernel `h` and a vector of length `N`.
"""
function convolution_matrix(h::AbstractVector{<:Number}, N::Integer)
    @assert N>0
    M = length(h)

    K = zeros(eltype(h), (M+N-1, M+N-1))
    hz = vcat(h[end:-1:1], zeros(eltype(h), N-1)) # zero padding of the kernel
    # K = hcat([circshift(hz,n) for n=0:M+N-2]...)
    for n=0:M+N-2
        K[:,n+1] = circshift(hz,n)
    end
    return K[M:end,:]'

    # Another more efficient implementation
    # for r=1:M
    #     for d=1:M+N-r
    #         K[d+r-1,d] = h[r]
    #     end
    # end
    # return K[:,1:N]
end

"""
    fullpastfilt(X::AbstractVector{<:Number}, filters::AbstractVector)

Filter a signal `X` of length `N` by `N` filters of increasing length from 1 to `N`.
"""
function fullpastfilt(X::AbstractVector{<:Number}, filters::AbstractVector, n::Integer=2, method::Symbol=:recursive)
    @assert length(X) >= length(filters)
    @assert all(length(f)==n for (n,f) in enumerate(filters))

#     Φ = V -> vcat(V[1], [sum(V[1:t].*filters[t]) for t=1:length(filters)])
    Φ = V -> [sum(V[1:t].*filters[t]) for t=1:length(filters)]
    Xϕ = Φ(X)

    # Experimental: phase-lag correction
    X1 = copy(X)
    for t=2:n
        if method==:recursive
            X1 += X - Xϕ
        else
            X1 += X - circshift(X, 1)
        end
        Xϕ = Φ(X1)
    end

    return Xϕ
end

"""
    conv_operator(X::AbstractVector{<:Real}, A::AbstractVector{<:Real}, b::Real=0.; nan::Symbol=:zero, mode::Symbol=:causal)

Apply the convolution kernel `A` and bias `b` to the input vector `X` and yield `A * X + b`.
"""
function conv_operator(X::AbstractVector{<:Real}, A::AbstractVector{<:Real}, b::Real=0.; nan::Symbol=:zero, mode::Symbol=:causal)
    X = copy(X)
    if nan==:zero
        X[isnan.(X)] .= 0
    end
    Y = native_conv(A, X) .+ b
    return (mode==:causal) ? Y[1:length(X)] : Y[end-length(X)+1:end]
end

"""
Compute the scaling and the wavelet function using the cascade algorithm.

The implementation follows the reference 2, but with a modified initialization.

# Returns
- ϕ, ψ, g: scaling, wavelet function and the associated sampling grid (computed for the Daubechies wavelet)

# References
- https://en.wikipedia.org/wiki/Daubechies_wavelet
- https://en.wikipedia.org/wiki/Cascade_algorithm
- http://cnx.org/content/m10486/latest/
"""
function wavefunc(lo::AbstractVector{<:Real}, hi::AbstractVector{<:Real}; level::Integer=10, nflag::Bool=true)
    @assert length(lo)==length(hi)

    # Initialization of the cascade algorithm
    # Method 1: using constant 1, this gives the best results (criteria of orthogonality etc.)
    ϕ = [1]
#
#     # Method 2: using the original filter
#     ϕ = copy(lo)
#
#     # Method 3: take one specific eigen vector of the decimated convolution matrix, see Reference 2.
#     K = convolution_matrix(lo, length(lo))[1:2:end, :]
#     μ, V = eig(K)  # such that (K * V) - V * diagm(μ) = 0
#     idx = find(abs.(μ - 1/sqrt(2)) .< 1e-3)[1]
#     ϕ = real(V[:, idx])

    # Normalization: this is necessary to get the correct numerical range
    ϕ ./= sum(ϕ)
    ψ = Float64[]

    fct = sqrt(2)  # scaling factor

    # Iteration of the cascade algorithm
    for n = 1:level
        # up-sampling of low-pass filter
        s = 2^(n-1)
        l = (length(lo)-1) * s + 1
        lo_up = zeros(Float64, l)
        lo_up[1:s:end] = lo

        if n==level
            # Last iteration only
            # up-sampling of high-pass filter
            hi_up = zeros(Float64, l)
            hi_up[1:s:end] = hi
            ψ = conv(hi_up, ϕ) * fct
        end
        ϕ = conv(lo_up, ϕ) * fct
    end

    # sampling grid:
    # the Daubechies wavelet of N vanishing moments has support [0, 2N-1] and its qmf filter has length 2N
    g = (length(lo)-1) * collect(0:(length(ϕ)-1))/length(ϕ)

    if nflag # force unit norm
        δ = g[2]-g[1]  # step of the sampling grid
        ϕ ./= sqrt(ϕ' * ϕ * δ)
        ψ ./= sqrt(ψ' * ψ * δ)
    end

    return ϕ, ψ, g
end

"""
Compute the masks for convolution for a mode of truncation.

For a input signal `x` and a kernel `h`, the full convolution x*h has length `length(x)+length(h)-1`. This function computes two masks:
- `kmask`: corresponding to the left/center/right part in x*h having the length of `x`
- `vmask`: corresponding to the valid coefficients (boundary free) in x*h

Note that `vmask` does not depend on `mode` and `vmask[kmask]` gives the mask of length of `x` corresponding to the valid coefficients in `kmask`.

# Args
- nx: length of input signal
- nh: length of kernel
- mode: `:left` (or `:causal`), `:right` (or `:anticausal`), `:center` (or `:valid`)

# Returns
- kmask, vmask

# Examples
```julia-repl
julia> y = conv(x, h)
julia> kmask, vmask = convmask(length(x), length(h), :center)
julia> tmask = vmask[kmask]
julia> y[kmask] # center part, same size as x
julia> y[vmask] # valide part, same as y[kmask][dmask]
```
"""
function convmask(nx::Integer, nh::Integer, mode::Symbol)
    kmask = zeros(Bool, nx+nh-1)

    if mode == :left || mode == :causal
        # kernel's support: 0:nh-1
        kmask[1:nx] .= true
    elseif mode == :right || mode == :anticausal
        # kernel's support: 1-nh:0
        kmask[nh:end] .= true  # or mask0[end-nx+1:end] = true
    elseif mode == :center || mode == :valid
        # kernel's support: 1-nh÷2 : nh-nh÷2
        m = max(1, floor(Int, nh/2))  # or nh÷2
        kmask[m:m+nx-1] .= true
    else
        error("Unknown mode: $mode")
    end

    vmask = zeros(Bool, nx+nh-1); vmask[nh:nx] .= true
    return kmask, vmask
end

"""
Down-sampling operator.
"""
function downsampling(x::AbstractVector, s::Integer=2)
    return x[1:s:end]
end

"""
Up-sampling operator.
"""
function upsampling(x::AbstractVector, s::Integer=2; tight::Bool=true)
    y = zeros(length(x)*s)
    y[1:s:end] = x
    return tight ? y[1:(end-(s-1))] : y
end

↑ = upsampling  # \uparrow
↓ = downsampling  # \downarrow

# ∗(x,y) = conv(x,y[end:-1:1])  # correlation, \ast
∗(x,y) = conv(x,y)  # convolution, \ast
⊛(x,y) = ↓(x ∗ ↑(y, 2, tight=true), 2)  # up-down convolution, \circledast

"""
Quadratic mirrored filter.

# Notes
- IMPORTANT: convention of qmf is (-1).^(0:length(h)-1) NOT (-1).^(1:length(h))
"""
function qmf(h::AbstractVector{<:Number})
    reverse(h .* (-1).^(0:length(h)-1))  # reverse(x) = x[end:-1:1]
end

"""
Dyadic scale stationary wavelet transform using python library `pywt`.

# Args
- x0: input signal
- wvl: name of wavelet, see `pywt.swt` documentation
- maxlevel: levels of decomposition
- mode: mode of zero-padding

# Returns
- (ac0, dc0, mask): ac: matrix of approximation coefficients with increasing scale index in column direction; dc: matrix of detail coefficients; mask: mask of zero-padding
- (ac, dc): arrays of approximation and detail coefficients without boundary elements

# Notes
- Original signal is zero-padded since `pywt.swt` requires signals of length 2^N
- The mode of convolution used by `pywt.swt` is not clear.
"""
function swt(x0::AbstractVector{<:Number}, wvl::String, maxlevel::Integer; mode::Symbol=:left)
    # @PyCall.pyimport pywt
    wavelet = pywt[:Wavelet](wvl)
    dec_lo, dec_hi, rec_lo, rec_hi = wavelet[:inverse_filter_bank]
    nf = length(dec_lo)  # filter length

    # zero padding
    n0 = length(x0)
    nx = 2^(ceil(Int, log2(n0)))  # length of zero-padded signal
    nlvl = min(nx, maxlevel)  # number of levels
    x = zeros(T, nx)
    # starting position for zero padding
    s::Int = if mode == :right
            1
        elseif mode == :center
            max((nx-n0)÷2,1)
        elseif mode == :left
            nx-n0+1
        else
            error("Unknown mode: $mode")
        end
    x[s:s+n0-1] = x0

    zmask = zeros(Bool, nx); zmask[s:s+n0-1] .= true
    # mask = ones(Bool, nx)

    w0 = pywt[:swt](x, wvl, level=nlvl)

    # ac = hcat([c[1][mask] for c in w0]...)
    # dc = hcat([c[2][mask] for c in w0]...)
    ac0 = hcat([c[1] for c in w0]...)
    dc0 = hcat([c[2] for c in w0]...)

    vmasks = []
    # ac, dc = [], []
    for s = 1:nlvl
        mask = zeros(Bool, nx)
        # nb = max(1, floor(Int, nf*2^(nlvl-s-1.)))  # number of boundary elements
        nb = nf*2^(nlvl-s)  # number of boundary elements
        idx = findall(zmask)[nb:end-nb+1]
        # idx = findall(zmask)[nb:end]
        mask[idx] .= true
        push!(vmasks, mask)
        # # truncation both sides
        # push!(ac, ac0[mask0,s][nb:end-nb+1])
        # push!(dc, dc0[mask0,s][nb:end-nb+1])
    end

    # return (ac0, dc0, mask0), (ac, dc, mask), nf
    return ac0, dc0, zmask, vmasks
end

"""
Inverse stationary wavelet transform using python library `pywt`.

# Args
- ac, dc: see outputs of `swt`
- wvl: name of wavelet, see `pywt.swt` documentation
- zmask: mask of zero-padding
"""
function iswt(ac::AbstractVector{<:Number}, dc::AbstractMatrix{<:Number}, wvl::String, zmask::AbstractVector{Bool})
    # @PyCall.pyimport pywt
    level = size(dc,2)
    w = [(ac, dc[:,1])]
    for n=2:level
        push!(w, (fill(0, size(ac)), dc[:,n]))
    end
    # w = [(emb(ac, mask), emb(dc[:,1], zmask))]
    # for n=2:level
    #     push!(w, (fill(0, size(ac)), emb(dc[:,n], zmask)))
    # end
    return pywt[:iswt](w, wvl)[findall(zmask)]
end

"""
Continuous wavelet transform based on quadrature.

# Args
- x: input signal
- wfunc: function for evaluation of wavelet at integer scales
- sclrng: range of integer scales
- mode: mode of convolution {:left, :center, :right} or {:causal, :valid, :anticausal}
"""
function cwt_quad(x::AbstractVector{<:Number}, wfunc::Function, sclrng::AbstractVector{<:Integer}, mode::Symbol)
    Ns = length(sclrng)
    Nx = length(x)

    dc = zeros((Nx, Ns))
    mc = zeros(Bool, (Nx, Ns))

    for (n,k) in enumerate(sclrng)
        f = wfunc(k) # get the function ψ(./k)
        km, vm = convmask(Nx, length(f), mode)

        # correlation via time-reversed convolution
        Y = conv(x, f[end:-1:1])
        dc[:,n] = Y[km] / √k
        mc[:,n] = vm[km]
    end
    return dc, mc
end

"""
Evaluate the wavelet function at integer scales by looking-up table.

# Args
- k: scale
- ψ: values of compact wavelet function on a grid
- Sψ: support of ψ
- v: desired number of vanishing moments of ψ, if >0 apply orthogonal projection.

# Returns
- f: the vector (ψ(n/k))_n such that n/k lies in Sψ.

# Note
For accuracy, increase the density of grid for pre-evaluation of ψ.
"""
function intscale_wavelet_filter(k::Integer, ψ::AbstractVector{<:Real}, Sψ::Tuple{Real,Real}, v::Integer=0)
    # @assert k > 0
    # @assert Sψ[2] > Sψ[1]

    Nψ = length(ψ)
    dh = (Sψ[2]-Sψ[1])/Nψ  # sampling step
    # @assert k < 1/dh  # upper bound of scale range

    Imin, Imax = ceil(Int, k*Sψ[1]), floor(Int, k*Sψ[2])
    idx = [max(1, min(Nψ, floor(Int, n/k/dh))) for n in Imin:Imax]
    f::Vector{Float64} = ψ[idx]

    # Forcing vanishing moments: necessary to avoid inhomogenity due to sampling of ψ
    # Projection onto the kernel of a under-determined Vandermonde matrix:
    if v>0
        V = vandermonde((length(f), v))'
        f -= V\(V*f)
    end

    return f
end

"""
Evaluate the wavelet function at integer scales.
"""
function intscale_wavelet_filter(k::Integer, ψ::Function, Sψ::Tuple{Real,Real}, v::Integer=0)
    # @assert k > 0
    # @assert Sψ[2] > Sψ[1]

    Imin, Imax = ceil(Int, k*Sψ[1]), floor(Int, k*Sψ[2])
    f::Vector{Float64} = ψ.((Imin:Imax)/k)

    # Forcing vanishing moments
    if v>0
        V = vandermonde((length(f), v))'
        f -= V\(V*f)
    end

    return f
end

mexhat(t::Real) = -exp(-t^2) * (4t^2-2t) / (2*sqrt(2π))

function intscale_mexhat_filter(k::Integer)
    return intscale_wavelet_filter(k, mexhat, (-5.,5.), 2)  # Mexhat has two vanishing moments
end

"""
Continous Mexican hat transform
"""
function cwt_mexhat(x::AbstractVector{<:Real}, sclrng::AbstractVector{<:Integer}, mode::Symbol=:center)
    return cwt_quad(x, intscale_mexhat_filter, sclrng, mode)
end

"""
Integer (even) scale Haar filter.

The original Haar wavelet takes value 1 on [0,1/2) and -1 on [1/2, 1) and 0 elsewhere.
"""
function intscale_haar_filter(s::Integer)
    @assert s > 0 && iseven(s)

    return vcat(ones(Float64, s÷2), -ones(Float64, s÷2))
end

"""
Continous Haar transform.
"""
function cwt_haar(x::AbstractVector{<:Real}, sclrng::AbstractVector{<:Integer}, mode::Symbol)
    all(iseven.(sclrng)) || error("Only even integer scale is admitted.")

    return cwt_quad(x, intscale_haar_filter, sclrng, mode)
end

# """
# Compute B-Spline filters of `v` vanishing moments at scale `2k`.
# """
# function bspline_filters(k::Int, v::Int)
#     @assert v>0
#     lo = vcat(ones(k), ones(k))
#     hi = vcat(ones(k),-ones(k))
#     # return col_normalize(wpt_filter(lo, hi, v-1)) * sqrt(k)
#     return col_normalize(biconv_filter(lo, hi, v-1)) * sqrt(k)
# end

"""
    intscale_bspline_filter(s::Int, v::Int)

Integer (even) scale B-Spline filter, which is defined as the auto-convolution of Haar filter.
"""
function intscale_bspline_filter(s::Integer, v::Integer)
    @assert s > 0 && iseven(s)
    @assert v > 0

    hi = vcat(ones(Float64, s÷2), -ones(Float64, s÷2))  # Haar filter at scale s
    ψv = reduce(∗, [hi for n=1:v]) / (s)^(v-1)

    return ψv
end

"""
    cwt_bspline(x::Vector{Float64}, sclrng::AbstractVector{Int}, v::Int, mode::Symbol)

Continous B-Spline transform at integer (even) scales.

# Args
- x: Input vector
- sclrng: vector of integer scales, all numbers must be even
- v: vanishing moment
- mode: mode of convolution

# TODO: parallelization
"""
function cwt_bspline(x::AbstractVector{<:Real}, sclrng::AbstractVector{<:Integer}, v::Integer, mode::Symbol)
    all(iseven.(sclrng)) || error("Only even integer scale is admitted.")

    # bsfilter = s->normalize(intscale_bspline_filter(s, v))
    bsfilter = s->intscale_bspline_filter(s, v)
    return cwt_quad(x, bsfilter, sclrng, mode)
end

"""
Evaluate the Fourier transform of a centered B-Spline wavelet.

# Notes
- Convention of Fourier transform:
    F(ω) = 1/√2π \\int f(t)e^{iωt} dt
- Convention of Fourier transform:
    f(t) = 1/√2π \\int F(ω)e^{-iωt} dω
"""
function _bspline_ft(ω::Real, v::Integer)
    #     @assert v>0  # check vanishing moment
    # # Version 1: non centered ψ. This works with convolution mode `:left`` and produces artefacts of radiancy straight lines.
    # return (ω==0) ? 0 : (2π)^((v-1)/2) * (-(1-exp(1im*ω/2))^2/(sqrt(2π)*1im*ω))^(v)  # non-centered bspline: supported on [0, v]
    # Version 2: centered ψ. This works with convolution mode `:center` which is non-causal add produces artefacts of radial circles.
    return (ω==0) ? 0 : (2π)^((v-1)/2) * (-(1-exp(1im*ω/2))^2/(sqrt(2π)*1im*ω) * exp(-1im*ω/2))^(v)  # centered bspline: supported on [-v/2, v/2]
end

#### B-Spline wavelet analysis ####

"""
    maxscale_bspline(N::Int, v::Int)

Maximum scale of B-Spline wavelet transform of `v` vanishing moments for a signal of length `N`.
"""
maxscale_bspline(N::Integer, v::Integer) = floor(Int, (N+1)/v/2)

"""
The integrand function of C^ψ_ρ(τ, H) with a centered B-spline wavelet.
"""
function Cψρ_bspline_integrand_center(τ::Real, ω::Real, ρ::Real, H::Real, v::Integer)
    return 1/(16^v * π) * (ω^2)^(v-(H+1/2)) * (sinc(ω*√ρ/4π)*sinc(ω/√ρ/4π))^(2v) * cos(ω*τ)
end

"""
Evaluate the integrand function of C^ψ_ρ(τ, H)

# Args
- τ, ω, ρ, H: see definition of C^ψ_ρ(τ, H).
- v: vanishing moments of the Bspline wavelet ψ, e.g. v=1 for Haar wavelet.

# Note
- We use the fact that C^ψ_ρ(τ, H) is a real function to simplify the implementation.
- In Julia the sinc function is defined as `sinc(x)=sin(πx)/(πx)`.
"""
function Cψρ_bspline_integrand(τ::Real, ω::Real, ρ::Real, H::Real, v::Integer, mode::Symbol)
    # The integrand is, by definition
    # _bspline_ft(ω*√ρ, v) * conj(_bspline_ft(ω/√ρ, v)) / abs(ω)^(2H+1) * exp(-1im*ω*τ)
    # this should be modulated by
    #    exp(-1im*ω*v*(√ρ-1/√ρ)/2), if the convolution mode is :left, i.e. causal
    # or exp(+1im*ω*v*(√ρ-1/√ρ)/2), if the convolution mode is :right, i.e. anti-causal
    # However, such implementation is numerically unstable due to singularities. We rewrite the function in an equivalent form using the sinc function which is numerically stable.
    if mode == :center || mode == :valid
        return Cψρ_bspline_integrand_center(τ, ω, ρ, H, v)
    elseif mode == :left || mode == :causal
        return Cψρ_bspline_integrand_center(τ+v*(√ρ-1/√ρ)/2, ω, ρ, H, v)
    elseif mode == :right || mode == :anticausal
        return Cψρ_bspline_integrand_center(τ-v*(√ρ-1/√ρ)/2, ω, ρ, H, v)
    else
        throw(UndefRefError("Unknown mode: $(mode)"))
    end
end

"""
Evaluate the function `C^ψ_ρ(τ,H)` by numerical integration.
"""
function Cψρ_bspline(τ::Real, ρ::Real, H::Real, v::Integer, mode::Symbol)
    @assert ρ>0
    @assert 1>H>0

    f = ω -> Cψρ_bspline_integrand(τ, ω, ρ, H, v, mode)

    # quadgk cannot handle the case v<=H+1/2: it may raise `DomainError with 0.0`
    return QuadGK.quadgk(f, (v>H+1/2) ? 0 : 1e-8, 50)[1]
end

# """
# Derivative w.r.t. H
# """
# function diff_Cψρ_bspline(τ::Real, ρ::Real, H::Real, v::Integer, mode::Symbol; rng::Tuple{Real, Real}=(0, 50))
#     f = ω -> ((ω==0) ? 0 : (-log(ω^2) * Cψρ_bspline_integrand(τ, ω, ρ, H, v, mode)))
#     return QuadGK.quadgk(f, rng...)[1]
# end

"""
Evaluate the function `A^ψ_ρ(τ, H)`

# Args
- τ, ρ, H: see definition
- v: vanishing moments of the wavelet ψ
- mode: {:left, :center, :right} for causal, centered, anti-causal ψ
"""
Aψρ_bspline(τ::Real, ρ::Real, H::Real, v::Int, mode::Symbol) = gamma(2H+1) * sin(π*H) * Cψρ_bspline(τ, ρ, H, v, mode)

"""
Evaluate the matrix of `A^ψ_ρ(H, τ)` for varing ρ with the B-Spline wavelet.

# Notes
- The true scale is two times the scale index due to the special implementation of B-Spline wavelet, see also `intscale_bspline_filter()`.
- Consider parallelization.
"""
function Amat_bspline(H::Real, v::Integer, lag::Real, sclrng::AbstractVector{<:Integer}, mode::Symbol)
    # all(iseven.(sclrng)) || error("Only even integer scale is admitted.")
    return [Aψρ_bspline(lag/sqrt(i*j), j/i, H, v, mode) for i in sclrng, j in sclrng]
    # return gamma(2H+1) * sin(π*H) * [Cψρ_bspline(lag/sqrt(i*j), j/i, H, v, mode) for i in sclrng, j in sclrng]
end

max_stars_repo_name: yanncalec/FracFin

max_issues_repo_name: yanncalec/FracFin

max_forks_repo_name: yanncalec/FracFin


text:
examples = joinpath(@__DIR__, "../examples")

@testset "Simulator tests" begin

    @testset "Add two numbers" begin
        program = load(joinpath(examples, "add-two-inputs.machine"))

        output = simulate!(copy(program), [3, 4])
        @test output[1] == 3 + 4

        output = simulate!(copy(program), [12, 8])
        @test output[1] == 12 + 8
    end

    @testset "Count down" begin
        program = load(joinpath(examples, "count-down.machine"))

        output = simulate!(copy(program), [3])
        @test output == [3, 2, 1, 0]

        output = simulate!(copy(program), [6])
        @test output == [6, 5, 4, 3, 2, 1, 0]
    end

    @testset "Asm and Sim - add two numbers" begin
        path = joinpath(examples, "add-two-inputs.lmc")
        program = assemble(path)
        output = simulate!(copy(program), [3, 4])
        @test output[1] == 3 + 4
    end

    @testset "Asm and Sim - count down" begin
        path = joinpath(examples, "count-down.lmc")
        program = assemble(path)
        output = simulate!(copy(program), [8])
        @test output == [8, 7, 6, 5, 4, 3, 2, 1, 0]
    end

    @testset "Maximizer" begin
        path = joinpath(examples, "maximizer.lmc")
        program = assemble(path)
        output = simulate!(copy(program), [8, 3])
        @test output == [8]

        output = simulate!(copy(program), [2,3, 4,5, 10,1, 8,5, 4,5])
        @test output == [3, 5, 10, 8, 5]
    end

    @testset "Multiply" begin
        path = joinpath(examples, "multiply.lmc")
        program = assemble(path)
        output = simulate!(copy(program), [3, 4])
        @test output == [12]

        output = simulate!(copy(program), [4, 5])
        @test output == [20]
    end

    # @testset "Factorial" begin
    #     program = load(joinpath(examples, "factorial.machine"))
    #
    #     output = simulate!(copy(program), [3])
    #     @test_broken  output == [6]
    #
    #     output = simulate!(copy(program), [4])
    #     @test_broken  output == [24]
    # end

end


text:
using SpgpGpu
using Test

@testset "SpgpGpu.jl" begin
    x = 2
    y = 2
    @test SpgpGpu.sum_values(x, y) == 4
end


text:
# This file includes the main test set of Causal.
# To include new tests, write your tests in files and save them in directories under `test` directory.

using Test
using Causal

using Logging
using Random
using JLD2, FileIO
using UUIDs
using Statistics
using LightGraphs
import Causal.process

# --------------------------- Deprecated --------------------------

function prepare(comp, kickoff::Bool=true)
    @warn "`prepare` function has been deprecated in favor of `equip`"
    equip(comp, kickoff)
end

# ---------------------------------- Include all test files ---------------------

# Construct the file tree in `test` directory.
filetree = walkdir(@__DIR__)
take!(filetree) # Pop the root directory `test` in which `runtests.jl` is.

# Include all test files under `test`
@time @testset "CausalTestSet" begin
    for (root, dirs, files) in filetree
        foreach(file -> include(joinpath(root, file)), files)
    end
end


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

mutable struct IoK8sApiCoreV1AttachedVolume <: SwaggerModel
    devicePath::Any # spec type: Union{ Nothing, String } # spec name: devicePath
    name::Any # spec type: Union{ Nothing, String } # spec name: name

    function IoK8sApiCoreV1AttachedVolume(;devicePath=nothing, name=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1AttachedVolume, Symbol("devicePath"), devicePath)
        setfield!(o, Symbol("devicePath"), devicePath)
        validate_property(IoK8sApiCoreV1AttachedVolume, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        o
    end
end # type IoK8sApiCoreV1AttachedVolume

const _property_map_IoK8sApiCoreV1AttachedVolume = Dict{Symbol,Symbol}(Symbol("devicePath")=>Symbol("devicePath"), Symbol("name")=>Symbol("name"))
const _property_types_IoK8sApiCoreV1AttachedVolume = Dict{Symbol,String}(Symbol("devicePath")=>"String", Symbol("name")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1AttachedVolume }) = collect(keys(_property_map_IoK8sApiCoreV1AttachedVolume))
Swagger.property_type(::Type{ IoK8sApiCoreV1AttachedVolume }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1AttachedVolume[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1AttachedVolume }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1AttachedVolume[property_name]

function check_required(o::IoK8sApiCoreV1AttachedVolume)
    (getproperty(o, Symbol("devicePath")) === nothing) && (return false)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1AttachedVolume }, name::Symbol, val)
end


text:
# backward compatibilty

@static if VERSION < v"1.1"
    function fieldtypes(T::Type)
        ntuple(fieldcount(T)) do idx
            fieldtype(T, idx)
        end
    end
end

# NOTE: for 1.0 compat
@static if !@isdefined(hasfield)
    function hasfield(T::Type, name::Symbol)
        return fieldindex(T, name, false) > 0
    end
end


text:
module DiscreteMarkovChains
import LinearAlgebra

include("Utils.jl")

include("MC.jl")
include("DMC.jl")
include("CMC.jl")

export DiscreteMarkovChain, state_space, transition_matrix, probability_matrix,
communication_classes, periodicities, decompose, canonical_form,
is_regular, is_ergodic, is_absorbing, is_reversible,
stationary_distribution, fundamental_matrix,
mean_time_to_absorption, exit_probabilities, first_passage_probabilities,
mean_recurrence_time, mean_first_passage_time

export ContinuousMarkovChain, embedded

end  # Module


text:
using Suppressor
using TypeClasses
using Test
using DataTypesBasic

state = Ref(0)
mycontext = @ContextManager cont -> begin
  state[] += 1
  returnvalue = cont(state[])
  returnvalue
end

three_times = mapn(mycontext, mycontext, mycontext) do x, y, z
  (x, y, z)
end

@test three_times(x -> x) == (1,2,3)
@test three_times(x -> x) == (4,5,6)
@test state[] == 6

cm(i=4, prefix="") = @ContextManager function(cont)
  println("$(prefix)before $i")
  r = cont(i)
  println("$(prefix)after $i")
  r
end

cm2 = flatten(map(x -> cm(x*x), cm(4)))
@test @suppress(cm2(x -> x)) == 16
@test splitln(@capture_out cm2(x -> x)) == [
  "before 4",
  "before 16",
  "after 16",
  "after 4"
]

nestedcm = @syntax_flatmap begin
  c1 = cm(1, "c1_")
  c2 = cm(c1 + 10, "c2_")
  @pure h = c1 + c2
  c3 = cm(h + 100, "c3_")
  @pure (c1, c2, c3)
end
@test @suppress(nestedcm(x -> x)) == (1, 11, 112)
@test splitln(@capture_out nestedcm(x -> x)) == [
  "c1_before 1",
  "c2_before 11",
  "c3_before 112",
  "c3_after 112",
  "c2_after 11",
  "c1_after 1"
]

@test pure(ContextManager, 3)(x -> x) == 3


text:
#= test/snorm.jl
=#

@testset "snorm" begin
  m = 20
  n = 10

  @time for T in (Float32, Float64, ComplexF32, ComplexF64)
    let A
      println("  $T")

      rtol = 100*eps(real(T))

      A = rand(T, m, n)
      nrm = opnorm(A)
      @test ≈(nrm, snorm(A); atol = rtol*nrm)

      A = rand(T, n, n)
      A += A'
      nrm = opnorm(A)
      @test ≈(nrm, snorm(A); atol = rtol*nrm)

      B = rand(T, size(A))
      nrm = opnorm(A - B)
      @test ≈(nrm, snormdiff(A, B); atol = rtol*nrm)
    end
  end
end


text:
# Execute POETS on EMT model
# DB 5/28/16

# Load POETS
using POETs

# process flags for cross validation
outdir = ARGS[1]
global HOLD=0
if length(matchall(r"fold_pop_obj",outdir)) == 1
global HOLD = parse(Int,split(outdir,"_")[6])
@show "Hold out", HOLD
end

# INITIALIZE MODEL
# FORMAT DATA AND INIT PARAMETER ARRAY
# INITIALIZE OBJECTIVE FUNCTION
# DECOMPOSE THIS INTO TWO FILES, ONE FOR INIT AND ANOTHER FOR FUNCS
fname = "../poets/emt_grn_model_objective_function.jl"
include(fname)

# Set parameters and pointers

initial_state = Array{Float64,1}(parameter_array)

const RANK_CUTOFF = 5.0
const MAX_NUM_ITER= 20

# Execute the beast
println(now())
(EC,PC,RA) = estimate_ensemble(objective_function,neighbor_function,acceptance_probability_function,cooling_function,initial_state;rank_cutoff=RANK_CUTOFF,maximum_number_of_iterations=MAX_NUM_ITER)
println(now())
# Save output

writedlm("$outdir/EC.dat", EC)
writedlm("$outdir/PC.dat", PC)
writedlm("$outdir/RA.dat", RA)