text:
#Functions
export save_BigData!, read_BigData!, update_NeoBigDataEntries!, read_MongoData, delete_MongoData

# --- Internal utility methods ---

"""
    _saveBigDataElement!(cg, vertex, bDE)

Insert or update the actual data payload into Mongo as required. Does not update Neo4j.
"""
function _saveBigDataElement!(cg::CloudGraph, bDE::BigDataElement)::Nothing
  saveTime = string(Dates.now(Dates.UTC));

  # 1. Get the collection name, if set...
  collection = cg.mongo.cgBindataCollection
  # DB customization for multitenancy
  if(haskey(bDE.sourceParams, "database") && haskey(bDE.sourceParams, "collection"))
      collection = cg.mongo.client[bDE.sourceParams["database"]][bDE.sourceParams["collection"]]
  elseif(haskey(bDE.sourceParams, "collection"))
      collection = cg.mongo.client[_mongoDefaultDb][bDE.sourceParams["collection"]]
  end

  # 2. Check if the key exists...
  isNew = true;
  if(bDE.sourceId != "")
    numNodes = length(collection, Mongoc.BSON("""{ "cgId": "$(bDE.sourceId)" }""") )
    isNew = numNodes == 0;
  end
  if(isNew)
    # Insert the node (additional parameters for user readability, data is all that matters)
    bsonDoc = Mongoc.BSON()
    bsonDoc["cgId"] = bDE.sourceId
    bsonDoc["elementId"] = bDE.id
    bsonDoc["description"] = bDE.description
    bsonDoc["data"] = bDE.data
    bsonDoc["mimeType"] = bDE.mimeType
    bsonDoc["neoNodeId"] = bDE.neoNodeId
    bsonDoc["lastSavedTimestamp"] = saveTime
    m_oid = push!(collection, bsonDoc)
    # m_oid = insert(collection, ("cgId" =>  bDE.sourceId, "elementId" => bDE.id, "description" => bDE.description, "data" => bDE.data, "mimeType" => bDE.mimeType, "neoNodeId" => bDE.neoNodeId, "lastSavedTimestamp" => saveTime))
    @info "Inserted big data to mongo id = $(m_oid) for cgId = $(bDE.id)"
  else
    # Update the node
    selector = Mongoc.BSON("""{ "cgId": "$(bDE.sourceId)" }""")
    bsonDoc = Mongoc.BSON()
    # bsonDoc["cgId"] = bDE.sourceId # this is the selector we are using, so don't update
    bsonDoc["elementId"] = bDE.id
    bsonDoc["description"] = bDE.description
    bsonDoc["data"] = bDE.data
    bsonDoc["mimeType"] = bDE.mimeType
    bsonDoc["neoNodeId"] = bDE.neoNodeId
    bsonDoc["lastSavedTimestamp"] = saveTime
    setDoc = Mongoc.BSON()
    setDoc["\$set"] = bsonDoc
    m_oid = Mongoc.update_one(collection, selector, setDoc)
    # m_oid = update(collection, ("cgId" => bDE.sourceId), set("elementId" => bDE.id, "description" => bDE.description, "data" => bDE.data, "mimeType" => bDE.mimeType, "neoNodeId" => bDE.neoNodeId, "lastSavedTimestamp" => saveTime))
    @info "Updated big data to mongo id (result=$(m_oid)) for cgId $(bDE.id)"
  end
  return nothing
end

# ------

"""
    update_NeoBigDataEntries!(cg, vertex)

Update the bigData dictionary elements in Neo4j. Does not insert or read from Mongo.
"""
function update_NeoBigDataEntries!(cg::CloudGraph, vertex::CloudVertex)::Nothing
  savedSets = Vector{BigDataRawType}(); #Vector{savedSets = BigDataRawType}();
  for elem in vertex.bigData.dataElements
    # keep big data separate during Neo4j updates and remerge at end
    push!(savedSets, elem.data);
    elem.data = Dict{String, Any}();
  end
  vertex.bigData.isExistingOnServer = true;
  vertex.bigData.lastSavedTimestamp = string(Dates.now(Dates.UTC));

  # Get the json bigData prop.
  bdProp = json(vertex.bigData);
  # Now put the data back
  i = 0;
  for elem in vertex.bigData.dataElements
    i += 1;
    elem.data = savedSets[i];
  end

  #Update the bigdata property
  setnodeproperty(vertex.neo4jNode, "bigData", bdProp);
  return(nothing)
end

function save_BigData!(cg::CloudGraph, vertex::CloudVertex)::Nothing
  #Write to Mongo
  for bDE in vertex.bigData.dataElements
    _saveBigDataElement!(cg, bDE);
  end

  #Now update the Neo4j node.
  update_NeoBigDataEntries!(cg, vertex)
  return(nothing)
end

function read_MongoData(cg::CloudGraph, bDE::BigDataElement)::BigDataRawType
    # 1. Get the collection name, if set...
    collection = cg.mongo.cgBindataCollection
    # DB customization for multitenancy
    if(haskey(bDE.sourceParams, "database") && haskey(bDE.sourceParams, "collection"))
        # collection = Mongo.MongoCollection(cg.mongo.client, bDE.sourceParams["database"], bDE.sourceParams["collection"]);
        collection = cg.mongo.client[bDE.sourceParams["database"]][bDE.sourceParams["collection"]]
    elseif(haskey(bDE.sourceParams, "collection"))
        # collection = Mongo.MongoCollection(cg.mongo.client, _mongoDefaultDb, bDE.sourceParams["collection"]);
        collection = cg.mongo.client[_mongoDefaultDb][bDE.sourceParams["collection"]]
    end

    # 2. See if the element exists
    numNodes = length(collection, Mongoc.BSON("""{ "cgId": "$(bDE.sourceId)" }""") )
    # numNodes = count(collection, ("cgId" => bDE.sourceId));
    if(numNodes != 1)
      error("The query for data elements named $(bDE.id) with Mongo CGID $(bDE.sourceId) returned $(numNodes) values, expected 1 result for this element!");
    end
    results = Mongoc.find_one(collection, Mongoc.BSON("""{ "cgId": "$(bDE.sourceId)" }""") )
    # results = first(find(collection, ("cgId" => bDE.sourceId)))

    #Have it, now parse it until we have a native binary or dictionary datatype.
    # If new type, convert back to dictionary

    if(typeof(results) == Mongoc.BSON)
        dict = Mongoc.as_dict(results)
        return dict["data"]      # testOutput = dict(results["data"]);
        # return convert(Dict{String, Any}, testOutput) #From {Any, Any} to a more comfortable stronger type
    else
        return results["data"];
    end
end

function delete_MongoData(cg::CloudGraph, bDE::BigDataElement)::Nothing
    # 1. Get the collection name, if set...
    collection = cg.mongo.cgBindataCollection
    # DB customization for multitenancy
    if(haskey(bDE.sourceParams, "database") && haskey(bDE.sourceParams, "collection"))
        collection = cg.mongo.client[bDE.sourceParams["database"]][bDE.sourceParams["collection"]]
        # collection = Mongo.MongoCollection(cg.mongo.client, bDE.sourceParams["database"], bDE.sourceParams["collection"]);
    elseif(haskey(bDE.sourceParams, "collection"))
        collection = cg.mongo.client[_mongoDefaultDb][bDE.sourceParams["collection"]]
        # collection = Mongo.MongoCollection(cg.mongo.client, _mongoDefaultDb, bDE.sourceParams["collection"]);
    end

    # 2. See if the element exists
    selector = Mongoc.BSON("""{ "cgId": "$(bDE.sourceId)" }""")
    numNodes = length(collection, selector);

    if(numNodes != 1)
      error("The query for data elements named $(bDE.id) with Mongo CGID $(bDE.sourceId) returned $(numNodes) values, expected 1 result for this element!");
    end
    Mongoc.delete_one(collection, selector)
    @info "Deleted big data mongo id = $(bDE.sourceId) for cgId = $(bDE.id)"

    return(nothing)
end

function read_BigData!(cg::CloudGraph, vertex::CloudVertex)::BigData
  if(vertex.bigData.isExistingOnServer == false)
    error("The data does not exist on the server. 'isExistingOnServer' is false. Have you saved with set_BigData!()");
  end
  for bDE in vertex.bigData.dataElements
      #TODO: Handle the different source types.
      bDE.data = read_MongoData(cg, bDE)
  end
  vertex.bigData.isRetrieved = true
  return(vertex.bigData)
end

function delete_BigData!(cg::CloudGraph, vertex::CloudVertex)::Nothing
  if(vertex.bigData.isExistingOnServer == false)
    error("The data does not exist on the server. 'isExistingOnServer' is false. Have you saved with set_BigData!()");
  end
  # Update structure now so if it fails midway and we save again it still writes a new set of keys.
  vertex.bigData.isExistingOnServer = false
  # Delete the data.
  for bDE in vertex.bigData.dataElements
      delete_MongoData(cg, bDE)
  end
  # Delete the entries and update the indices
  empty!(vertex.bigData.dataElements)
  update_NeoBigDataEntries!(cg, vertex)
  return(nothing)
end


text:
using KernelAbstractions, CUDAKernels, Test, CUDA

if has_cuda_gpu()
    CUDA.allowscalar(false)
end

# Simple kernel for matrix multiplication
@kernel function vec_add_kernel!(a, b, c)
    i, j = @index(Global, NTuple)

    c[i,j] = a[i,j] + b[i,j]
end

# Creating a wrapper kernel for launching with error checks
function vec_add!(a, b, c)
    if size(a)[2] != size(b)[1]
        println("Matrix size mismatch!")
        return nothing
    end
    if isa(a, Array)
        kernel! = vec_add_kernel!(CPU(),4)
    else
        kernel! = vec_add_kernel!(CUDADevice(),256)
    end
    kernel!(a, b, c, ndrange=size(c))
end

@testset "CPU+GPU kernels" begin
    a = rand(256, 256)
    b = rand(256, 256)
    c = zeros(256, 256)

    # beginning CPU tests, returns event
    ev = vec_add!(a,b,c)
    wait(ev)

    @test isapprox(c, a+b)

    # beginning GPU tests
    if has_cuda_gpu()
        d_a = CuArray(a)
        d_b = CuArray(b)
        d_c = CuArray(c)

        ev = vec_add!(d_a, d_b, d_c)
        wait(ev)

        @test isapprox(Array(d_c), a+b)
    end
end

ext: jl

max_stars_repo_name: leios/simuleios


text:
@doc doc"""
    mdtable(array; latex=true, head=[], side=[], transpose=false)

Latexify the input and output a markdown-formatted table.

```julia
julia> using Latexify
juila> copy_to_clipboard(true)
julia> M = ["x/y" 1//2
            "p_m" "e^2"]
julia> mdtable(M)
```

| $\frac{x}{y}$ | $\frac{1}{2}$ |
| -------------:| -------------:|
|       $p_{m}$ |       $e^{2}$ |

```julia
julia> head = ["Column 1", "Column 2"]
julia> side = ["Row 1", "Row 2"]
julia> mdtable(M; head=head, side=side)
```

|     . |      Column 1 |      Column 2 |
| -----:| -------------:| -------------:|
| Row 1 | $\frac{x}{y}$ | $\frac{1}{2}$ |
| Row 2 |       $p_{m}$ |       $e^{2}$ |

The value in the top right corner can be set if you let the `side` vector be one element
larger than the number of rows of your input:

```julia
julia> side = ["Corner", "Row 1", "Row 2"]
julia> mdtable(M; head=head, side=side)
```

| Corner |      Column 1 |      Column 2 |
| ------:| -------------:| -------------:|
|  Row 1 | $\frac{x}{y}$ | $\frac{1}{2}$ |
|  Row 2 |       $p_{m}$ |       $e^{2}$ |

The `head` and `side` vectors are not latexifed, but you can easily do this yourself:

```julia
julia> head = ["p_1", "p_2"]
julia> mdtable(M; head=latexinline(head))
```

|       $p_{1}$ |       $p_{2}$ |
| -------------:| -------------:|
| $\frac{x}{y}$ | $\frac{1}{2}$ |
|       $p_{m}$ |       $e^{2}$ |
"""
function mdtable end

function mdtable(M::AbstractMatrix; latex::Bool=true, escape_underscores=false, head=[], side=[], transpose=false, kwargs...)
    transpose && (M = permutedims(M, [2,1]))
    if latex
        M = _latexinline.(M; kwargs...)
    elseif haskey(kwargs, :fmt)
        formatter = kwargs[:fmt] isa String ? PrintfNumberFormatter(kwargs[:fmt]) : kwargs[:fmt]
        M = map(x -> x isa Number ? formatter(x) : x, M)
    end

    if !isempty(head)
        M = vcat(safereduce(hcat, head), M)
        @assert length(head) == size(M, 2) "The length of the head does not match the shape of the input matrix."
    end
    if !isempty(side)
        length(side) == size(M, 1) - 1 && (side = [LaTeXString("∘"); side])
        @assert length(side) == size(M, 1) "The length of the side does not match the shape of the input matrix."
        M = hcat(side, M)
    end

    t = "| " * join(M[1,:], " | ") * " |\n"
    size(M, 1) > 1 && (t *= "| ---  "^(size(M,2)-1) * "| --- |\n")
    for i in 2:size(M,1)
        t *= "| " * join(M[i,:], " | ") * " |\n"
    end

    escape_underscores && (t = replace(t, "_"=>"\\_"))
    t = Markdown.parse(t)
    COPY_TO_CLIPBOARD && clipboard(t)
    return t
end

mdtable(v::AbstractArray; kwargs...) = mdtable(reshape(v, (length(v), 1)); kwargs...)
mdtable(v::AbstractArray...; kwargs...) = mdtable(safereduce(hcat, v); kwargs...)
mdtable(d::AbstractDict; kwargs...) = mdtable(collect(keys(d)), collect(values(d)); kwargs...)
mdtable(arg::Tuple; kwargs...) = mdtable(safereduce(hcat, [collect(i) for i in arg]); kwargs...)


text:
using MathLink
using Test

import MathLink: WExpr, WSymbol

@testset "integers" begin
    w = W"Factorial"(30)
    @test_throws MathLink.MathLinkError weval(Int, w)
    @test weval(BigInt, w) == factorial(big(30))
    @test weval(W"Factorial"(20)) === factorial(20)

    @test weval(Float64, W"N"(W"Log"(factorial(big(30))),100)) == log(Float64(factorial(big(30))))
    @test weval(BigFloat, W"N"(W"Log"(factorial(big(30))),100)) == log(factorial(big(30)))
end

@testset "floats" begin
    w = W"Sqrt"(2.0)
    @test weval(w) == sqrt(2.0)
    w = W"Sqrt"(2f0)
    @test weval(Float32, w) == sqrt(2f0)
end

@testset "BigFloats" begin
    @test weval(BigFloat, W"N"(W"Pi",200)) == big(pi)
    @test abs(weval(Float64, W"Sin"(big(pi)))) < sin(big(pi))
end

@testset "Irrationals" begin
    @test weval(W"Sin"(pi)) == 0
end

@testset "expressions" begin
    @test weval(W"Function"(W"x",W"Times"(W"x", 2))(100)) == 200

    @test weval(W"Integrate"(W"Log"(W"x"), (W"x", 1, W"E"))) == 1

    @test weval(W`Integrate[Log[x], {x,1,E}]`) == 1
end


text:
using CUDA.CUDNN: cudnnPoolingMode_t, CUDNN_POOLING_MAX,
                  CUDNN_POOLING_AVERAGE_COUNT_INCLUDE_PADDING,
                  cudnnPoolingForward!, pooldims, cudnnPoolingBackward

import NNlib: maxpool!, ∇maxpool!, meanpool!, ∇meanpool!
import CUDA.CUDNN: cudnnPoolingDescriptor

function cudnnPoolingDescriptor(pdims::PoolDims, x::DenseCuArray{T}, mode::cudnnPoolingMode_t) where T
    window, padding, stride = NNlib.kernel_size(pdims), nnlibPadding(pdims), NNlib.stride(pdims)
    nanOpt = CUDNN_NOT_PROPAGATE_NAN
    cudnnPoolingDescriptor(mode, nanOpt, Cint(ndims(x)-2), pooldims(window,size(x)), pooldims(padding,size(x)), pooldims(stride,size(x)))
end

function maxpool!(y::DenseCuArray{T}, x::DenseCuArray{T}, pdims::PoolDims) where T<:CUDNNFloat
    d = cudnnPoolingDescriptor(pdims, x, CUDNN_POOLING_MAX)
    cudnnPoolingForward!(y, x, d)
end

function ∇maxpool!(dx::DenseCuArray{T}, dy::DenseCuArray{T}, y::DenseCuArray{T}, x::DenseCuArray{T}, pdims::PoolDims) where T<:CUDNNFloat
    xDesc, yDesc = cudnnTensorDescriptor.((x, y))
    d = cudnnPoolingDescriptor(pdims, x, CUDNN_POOLING_MAX)
    alpha, beta = scalingParameter(T,1), scalingParameter(T,0)
    cudnnPoolingBackward(handle(), d, alpha, yDesc, y, yDesc, dy, xDesc, x, beta, xDesc, dx)
    return dx
end

function meanpool!(y::DenseCuArray{T}, x::DenseCuArray{T}, pdims::PoolDims) where T<:CUDNNFloat
    d = cudnnPoolingDescriptor(pdims, x, CUDNN_POOLING_AVERAGE_COUNT_INCLUDE_PADDING)
    cudnnPoolingForward!(y, x, d)
end

function ∇meanpool!(dx::DenseCuArray{T}, dy::DenseCuArray{T}, y::DenseCuArray{T}, x::DenseCuArray{T}, pdims::PoolDims) where T<:CUDNNFloat
    xDesc, yDesc = cudnnTensorDescriptor.((x, y))
    d = cudnnPoolingDescriptor(pdims, x, CUDNN_POOLING_AVERAGE_COUNT_INCLUDE_PADDING)
    alpha, beta = scalingParameter(T,1), scalingParameter(T,0)
    cudnnPoolingBackward(handle(), d, alpha, yDesc, y, yDesc, dy, xDesc, x, beta, xDesc, dx)
    return dx
end

### Since CUDA.jl does not support 1D pooling, we have to convert to 2d

add1d(x) = reshape(x, 1, size(x)...)

function fix_pooldims_1d(pdims::PoolDims{1,K,S,P,D}) where {K,S,P,D}
    PoolDims{2, K + 1, S + 1, P + 2, D + 1}((1, NNlib.input_size(pdims)...),
                                            (1, NNlib.kernel_size(pdims)...),
                                            NNlib.channels_in(pdims),
                                            (1, NNlib.stride(pdims)...),
                                            (0, 0, NNlib.padding(pdims)...),
                                            (1, NNlib.dilation(pdims)...))
end

function maxpool!(y::DenseCuArray{T,3}, x::DenseCuArray{T,3}, pdims::PoolDims) where T<:CUDNNFloat
    maxpool!(add1d(y), add1d(x), fix_pooldims_1d(pdims))
    return y
end

function meanpool!(y::DenseCuArray{T,3}, x::DenseCuArray{T,3}, pdims::PoolDims) where T<:CUDNNFloat
    meanpool!(add1d(y), add1d(x), fix_pooldims_1d(pdims))
    return y
end

function ∇maxpool!(dx::DenseCuArray{T,3}, dy::DenseCuArray{T,3}, y::DenseCuArray{T,3}, x::DenseCuArray{T,3}, pdims::PoolDims) where T<:CUDNNFloat
    ∇maxpool!(add1d(dx), add1d(dy), add1d(y), add1d(x), fix_pooldims_1d(pdims))
    return dx
end

function ∇meanpool!(dx::DenseCuArray{T,3}, dy::DenseCuArray{T,3}, y::DenseCuArray{T,3}, x::DenseCuArray{T,3}, pdims::PoolDims) where T<:CUDNNFloat
    ∇meanpool!(add1d(dx), add1d(dy), add1d(y), add1d(x), fix_pooldims_1d(pdims))
    return dx
end


text:
struct Implicit{M1<:AbstractMatrix, M2<:AbstractMatrix, A<:AbstractArray, V<:AbstractVector}
    # Constants for backwards implicit biharmonic diffusion
    dmp1::M1
    dmp1d::M1
    dmp1s::M1

    xa::M2
    xb::M2
    xc::M2
    xd::M2
    xe::M2
    xf::A
    xj::A

    tref::V
    tref1::V
    tref2::V
    tref3::V

    elz::M1
    σ_thick_x::V
end

# Initialize constants for the implicit gravity wave computation.
# It is assumed that that all implicit steps are of length 2*Δt and use
# the forward/backward parameter α. initialize_implicit has to be re-called
# whenever either of these 2.0 parameters is changed. initialize_implicit should
# be called even if the explicit option is chosen for the gravity wave
# terms (the reference state temperature tref is subtracted from some
# terms anyway to reduce roundoff error; also the constants needed for
# the biharmonic diffusion, which is assumed always to be backwards
# implicit, are defined in initialize_implicit)
function Implicit(T, geometry::Geometry, constants::Constants, params::Parameters,
                  horizontal_diffusion::HorizontalDiffusion, Δt)
    @unpack nlev, mx, nx, σ_half, σ_full, σ_thick = geometry
    @unpack Rₑ, g, akap, R, γ = constants
    @unpack α = params
    @unpack dmp, dmpd, dmps = horizontal_diffusion

    dmp1 = zeros(T, mx, nx)
    dmp1d = zeros(T, mx, nx)
    dmp1s = zeros(T, mx, nx)

    xa = zeros(T, nlev, nlev)
    xb = zeros(T, nlev, nlev)
    xc = zeros(T, nlev, nlev)
    xd = zeros(T, nlev, nlev)
    xe = zeros(T, nlev, nlev)
    xf = zeros(T, nlev, nlev, mx+nx+1)
    xj = zeros(T, nlev, nlev, mx+nx+1)

    tref = zeros(T, nlev)
    tref1 = zeros(T, nlev)
    tref2 = zeros(T, nlev)
    tref3 = zeros(T, nlev)

    elz = zeros(T, mx, nx)
    σ_thick_x = zeros(T, nlev)

    # 1. Constants for backwards implicit biharmonic diffusion
    dmp1  = 1.0./(1.0 .+ dmp*Δt)
    dmp1d = 1.0./(1.0 .+ dmpd*Δt)
    dmp1s = 1.0./(1.0 .+ dmps*Δt)

    # 1. Constants for implicit gravity wave computation
    # reference atmosphere, function of sigma only
    γ_g = γ/(1000.0*g)

    tref = 288.0max.(0.2, σ_full).^(R*γ_g)
    tref1 = R*tref
    tref2 = akap*tref
    tref3 = σ_full.*tref

    # Other constants
    xi = Δt*α
    xxi = xi/(Rₑ^2.0)
    σ_thick_x = xi*σ_thick

    for n in 1:nx
        for m in 1:mx
            elz[m,n] = T(m + n - 2)*T(m + n - 1)*xxi
        end
    end

    #T(K) = TEX(K)+YA(K,K')*D(K') + XA(K,K')*SIG(K')
    ya = zeros(T, nlev, nlev)
    for k in 1:nlev
        for k1 in 1:nlev
            ya[k,k1] = -akap*tref[k]*σ_thick[k1]
        end
    end

    for k in 2:nlev
        xa[k,k-1] = 0.5*(akap*tref[k]/σ_full[k] - (tref[k] - tref[k-1])/σ_thick[k])
    end

    for k in 1:nlev-1
        xa[k,k] = 0.5*(akap*tref[k]/σ_full[k] - (tref[k+1] - tref[k])/σ_thick[k])
    end

    #sig(k)=xb(k,k')*d(k')
    dsum = zeros(T, nlev)
    dsum[1] = σ_thick[1]
    for k in 2:nlev
        dsum[k] = dsum[k-1] + σ_thick[k]
    end

    for k in 1:nlev-1
        for k1 in 1:nlev
            xb[k,k1] = σ_thick[k1]*dsum[k]
            if k1 <= k
                xb[k,k1] = xb[k,k1] - σ_thick[k1]
            end
        end
    end

    #t(k)=tex(k)+xc(k,k')*d(k')
    for k in 1:nlev
        for k1 in 1:nlev
            xc[k,k1] = ya[k,k1]
            for k2 in 1:nlev-1
                xc[k,k1] = xc[k,k1] + xa[k,k2]*xb[k2,k1]
            end
        end
    end

    #P(K)=XD(K,K')*T(K')
    for k in 1:nlev
        for k1 in k+1:nlev
            xd[k,k1] = R*log(σ_half[k1+1]/σ_half[k1])
        end
    end
    for k in 1:nlev
        xd[k,k] = R*log(σ_half[k+1]/σ_full[k])
    end

    #P(K)=YE(K)+XE(K,K')*D(K')
    for k in 1:nlev
        for k1 in 1:nlev
            for k2 in 1:nlev
                xe[k,k1] = xe[k,k1] + xd[k,k2]*xc[k2,k1]
            end
        end
    end

    for l in 1:mx+nx+1
        xxx = (T(l)*T(l+1))/(Rₑ^2.0)
        for k in 1:nlev
            for k1 in 1:nlev
                xf[k,k1,l] = xi^2.0*xxx*(R*tref[k]*σ_thick[k1] - xe[k,k1])
            end
        end
        for k in 1:nlev
            xf[k,k,l] = xf[k,k,l] + 1.0
        end
    end

    for l in 1:mx+nx+1
        xj[:,:,l] = inv(xf[:,:,l])
    end

    for k in 1:nlev
        for k1 in 1:nlev
            xc[k,k1] = xc[k,k1]*xi
        end
    end

    Implicit(dmp1, dmp1d, dmp1s, xa, xb, xc, xd, xe, xf, xj, tref, tref1, tref2, tref3,
             elz, σ_thick_x)
end

# Correct tendencies for implicit gravity wave model
#  Input/output : D_tend  = divergence tendency
#                 Tₐ_tend = temperature tendency
#                 pₛ_tend = tendency of log(surface pressure)

"""
Compute the tendencies for implicit gravity wave model
"""
function implicit_terms!(div_tend, t_tend, ps_tend)
    ye = zeros(Complex{T}, mx, nx, nlev)
    yf = zeros(Complex{T}, mx, nx, nlev)

    for k1 in 1:nlev
        for k in 1:nlev
            ye[:,:,k] = ye[:,:,k] + xd[k,k1]*Tₐ_tend[:,:,k1]
        end
    end

    for k in 1:nlev
        ye[:,:,k] = ye[:,:,k] + tref1[k]*pₛ_tend
    end

    for k in 1:nlev
        for m in 1:mx
            for n in 1:nx
                yf[m,n,k] = D_tend[m,n,k] + elz[m,n]*ye[m,n,k]
            end
        end
    end

    D_tend *= zero

    for n in 1:nx
        for m in 1:mx
            if (m + n - 2) != 0
                for k1 in 1:nlev
                    D_tend[m,n,:] = D_tend[m,n,:] + xj[:,k1,m+n-2]*yf[m,n,k1]
                end
            end
        end
    end

    for k in 1:nlev
        pₛ_tend = pₛ_tend - D_tend[:,:,k]*σ_thick_x[k]
    end

    for k in 1:nlev
        for k1 in 1:nlev
            Tₐ_tend[:,:,k] = Tₐ_tend[:,:,k] + xc[k,k1]*D_tend[:,:,k1]
        end
    end
end


text:
using Gemini
using Test

@testset "Gemini.jl" begin
  printstyled(color=:blue, "marketdata_v2\n")
  @testset "marketdata_v2" begin
    include("marketdata_v2.jl");sleep(1)
  end
  printstyled(color=:blue, "symbols\n")
  @testset "symbols" begin
    include("symbols.jl");sleep(1)
  end
  # manual testing for now, until I can mock up an HTTP server to return dummy requests.
  # printstyled(color=:blue, "orders\n")
  # @testset "orders" begin
  #   include("orders.jl");sleep(1)
  # end
end

max_stars_repo_name: rory-linehan/Gemini

max_issues_repo_name: rory-linehan/Gemini

max_forks_repo_name: rory-linehan/Gemini


text:
function reverse_autodiff_gradient(
  result::DiffResults.DiffResult,
  f::Union{ReverseDiff.GradientTape, ReverseDiff.CompiledTape},
  x::Vector
)
  ReverseDiff.gradient!(result, f, x)
  return DiffResults.gradient(result)
end

function reverse_autodiff_negative_hessian(
  result::DiffResults.DiffResult,
  f::Union{ReverseDiff.HessianTape, ReverseDiff.CompiledTape},
  x::Vector
)
  ReverseDiff.hessian!(result, f, x)
  return -DiffResults.hessian(result)
end

function reverse_autodiff_upto_gradient(
  result::DiffResults.DiffResult,
  f::Union{ReverseDiff.GradientTape, ReverseDiff.CompiledTape},
  x::Vector
)
  ReverseDiff.gradient!(result, f, x)
  return DiffResults.value(result), DiffResults.gradient(result)
end

function reverse_autodiff_upto_negative_hessian(
  result::DiffResults.DiffResult,
  f::Union{ReverseDiff.HessianTape, ReverseDiff.CompiledTape},
  x::Vector
)
  ReverseDiff.hessian!(result, f, x)
  return DiffResults.value(result), DiffResults.gradient(result), -DiffResults.hessian(result)
end


text:
# p31.jl - gamma function via complex integral, trapezoid rule

N = 70; theta = -pi + (2*pi/N)*(.5:N-.5);
c = -11;                     # center of circle of integration
r = 16;                      # radius of circle of integration
x = -3.5:.1:4; y = -2.5:.1:2.5;
zz = x' .+ 1im*y; gaminv = 0*zz;
for i = 1:N
    t = c + r*exp(1im*theta[i]);
    gaminv += exp(t)*t.^(-zz)*(t-c);
end
gaminv = gaminv/N; gam = 1./gaminv;
clf(); surf(x,y,abs.(gam));
xlim(-3.5,4); ylim(-2.5,2.5); zlim(0,6);
xlabel("Re(z)"); ylabel("Im(z)");
text3D(4,-1.4,5.5,"\$|\\Gamma(z)|\$",fontsize=20);


text:
# Identification and resolution of Plant taxonomic names
# Nov 8th, 2018
# Author: Vivekanand Sharma, PhD
# Postdoctoral Research Associate
# Center for Biomedical Informatics
# Brown University
# Pre-processing involves text tokenization and noun phrase parsing.

using DataFrames
using HTTP
using JSON
using PyCall
#using HttpServer

@pyimport nltk as ntk
@pyimport textblob
@pyimport textblob.np_extractors as npx
#@pyimport textblob.taggers as ptag

include("sw_align_srs.jl")

# Escaping special characters
# replace(str2, r"[+|-|\^|\"|\~|\*|\?|\:|\|/|\!|\&|\|\(|\)]", " OR ")

const solr_host = "http://127.0.0.1:8983/solr"
#Sentence tokenization
function sentTokenization(txt, ntk)
    sent_arr = ntk.sent_tokenize(txt)
    for i in 1:length(sent_arr)
        if i < length(sent_arr)
            if occursin(r"[Nn]o\.\s[a-z0-9]", sent_arr[i] * " " * sent_arr[i+1]) || occursin(r"[A-Za-z]\.\s[a-z0-9]", sent_arr[i] * " " * sent_arr[i+1])
                sent_arr[i] = sent_arr[i] * " " * sent_arr[i+1] #* "\."
                deleteat!(sent_arr, i+1)
            end
        end
    end
    if sent_arr[length(sent_arr)] == "" || isempty(sent_arr[length(sent_arr)])
        deleteat!(sent_arr, length(sent_arr))
    end
    return sent_arr
end

# Generate query URL from herb name string
function generate_query(string::String)
    string1 = replace(string, r"[\+|\-|\^|\"|\~|\*|\?|\:|\|/|\!|\&|\|\(|\)|\{|\}|\[|\]]" => " ")
    string1 = replace(string1, r"[0-9]" => "")
    string1 = replace(string1, r"\s{2,}" => " ")
    string1 = strip(string1)
    str_arr = Any[]
    query   = ""
    for str in split(string1, " ")
        push!(str_arr, "term:$str~")
        query = join(str_arr, "%20OR%20")
    end
    url_string = solr_host * "/ubt_plant/select?fl=*,score&q=$query&sort=score%20desc"
    #url_string = "http://bcbi.brown.edu/solr/solr/ubt_plant/select?fl=*,score&q=$query&sort=score%20desc"
    return url_string
end

# Generate recursive query URL from herb name string
function generate_sub_query(string::String, not_arr)
    string1 = replace(string, r"[\+|\-|\^|\"|\~|\*|\?|\:|\|/|\!|\&|\|\(|\)|\{|\}|\[|\]\,]" => " ")
    string1 = replace(string1, r"[0-9]" => "")
    string1 = replace(string1, r"\s{2,}" => " ")
    string1 = strip(string1)
    str_arr = Any[]
    yes_query   = ""
    not_query = ""
    for str in split(string1, " ")
        push!(str_arr, "term:$str~")
        yes_query = join(str_arr, "%20OR%20")
    end
    if length(not_arr) > 0
        name_arr = String[]
        for name in not_arr
            #println("$name \=\>")
            for str in split(name, " ")
                push!(name_arr, "term:$str")
                not_query = join(name_arr, "%20OR%20")
            end
        end
        query = "($yes_query)%20NOT%20($not_query)"
        #println(query)
    else
        query = yes_query
    end

    url_string = solr_host * "/ubt_plant/select?fl=*,score&q=$query&sort=score%20desc&rows=10"
    @info "***Solr URL"
    @show url_string
    #url_string = "http://bcbi.brown.edu/solr/solr/ubt_plant/select?fl=*,score&q=$query&sort=score%20desc"
    return url_string
end

# Generate optimal matches for query string
function norm_string(url::String,herb_str::String)
    herb_name = herb_str
    norm_str = ""
    #flag     = 0
    ret_str  = HTTP.request("GET","$url";)
    json_str = JSON.parse(String(ret_str.body))
    @info "***JSON response"
    @show json_str
    numFound = json_str["response"]["numFound"]
    #println(json_str)
    if numFound > 0
        herb_str = replace(herb_str, r"\s+cf\.?\s+"i => " ")
        herb_str = replace(herb_str, r"\s+aff\.?\s+"i => " ")
        herb_str = replace(herb_str, r"\s{2,}" => " ")
        herb_str = strip(herb_str)
        res_dict = Dict{Float64,String}()
        max_scr = 0
        numFound > 10 ? n = 10 : n = numFound

        for i=1:n
            ubt_id   = json_str["response"]["docs"][i]["ubt_id"][1]
            name     = json_str["response"]["docs"][i]["term"][1]
            solr_scr = round(json_str["response"]["docs"][i]["score"][1], digits = 2)
            typ      = json_str["response"]["docs"][i]["type"][1]
            seq_b = uppercase(herb_str)
            seq_a = uppercase(name)
            # Old version:
            #seq_a = uppercase(herb_str)
            #seq_b = uppercase(name)
            if occursin(r"([A-Z][a-z\-]+\s[a-z]{0,3}\.?\s*[a-z\-]+)\.?.*", name)
                matrix, path = sw_align(seq_a,seq_b)
                algn_a, algn_b, scr, indx, indy, mstr = traceback(matrix,path,seq_a,seq_b)
                name_pd  = name
                name_len = length(name)  #length(split(herb_str, " ")[1])
                algn_len = length(algn_a)

                coverage = round((algn_len/name_len)*100, digits = 2)
                if coverage == 100.00
                    mapping = check_match(seq_a, seq_b, algn_a, indx, mstr)
                    #scr > max_scr ? max_scr = scr : max_scr = max_scr
                    if mapping == "match"
                        if scr > max_scr
                            max_scr = scr
                            filter!(x -> x!="-", algn_b)
                            res_dict[scr] = "$(uppercasefirst(lowercase(join(algn_b))))\$$ubt_id\$$name\$$solr_scr\$$scr\$$name_len\$$coverage\$$typ"
                            #println("$herb_name\$$name\$$name_pd\$$solr_scr\$$scr\$$name_len\$$coverage\$$mapping")
                        end
                    end
                else
                    name_pd  = match(r"([A-Z][a-z\-]+\s[a-z]{0,3}\.?\s*[a-z\-]+)\.?.*", name).captures[1]
                    name_pd  = String(name_pd)

                    seq_a    = uppercase(name_pd)
                    matrix, path = sw_align(seq_a,seq_b)
                    algn_a, algn_b, scr, indx, indy, mstr = traceback(matrix,path,seq_a,seq_b)
                    name_len = length(name_pd)
                    algn_len = length(algn_a)
                    coverage = round((algn_len/name_len)*100, digits = 2)
                    if coverage > 80.00
                        mapping = check_match(seq_a, seq_b, algn_a, indx, mstr)
                        if mapping == "match"
                            if scr > max_scr
                                max_scr = scr
                                filter!(x -> x!="-", algn_b)
                                res_dict[scr] = "$(uppercasefirst(lowercase(join(algn_b))))\$$ubt_id\$$name\$$solr_scr\$$scr\$$name_len\$$coverage\$$typ"
                                #println("$herb_name\$$name\$$name_pd\$$solr_scr\$$scr\$$name_len\$$coverage\$$mapping")
                            end
                        end
                    end
                end
            else
                matrix, path = sw_align(seq_a,seq_b)
                algn_a, algn_b, scr, indx, indy, mstr = traceback(matrix,path,seq_a,seq_b)
                name_len = length(name)  #length(split(herb_str, " ")[1])
                algn_len = length(algn_a)

                coverage = round((algn_len/name_len)*100, digits = 2)
                if coverage > 80.00
                    mapping = check_match(seq_a, seq_b, algn_a, indx, mstr)
                    if mapping == "match"
                        if scr > max_scr
                            max_scr = scr
                            filter!(x -> x!="-", algn_b)
                            res_dict[scr] = "$(uppercasefirst(lowercase(join(algn_b))))\$$ubt_id\$$name\$$solr_scr\$$scr\$$name_len\$$coverage\$$typ"
                            #println("$herb_name\$$name\$$name\$$solr_scr\$$scr\$$name_len\$$coverage\$$mapping")
                        end
                    end
                end
            end
        end
        return get(res_dict, max_scr, "$herb_name\$NA\$NA\$NA\$NA\$NA\$NA\$NA"), numFound
    else
        return "$herb_name\$NA\$NA\$NA\$NA\$NA\$NA\$NA", numFound
    end
end

# Check for match or no-match
function check_match(seq_a, seq_b, algn_a, indx, mstr)
    cs1 = indx-length(mstr)+1
    # ce1 = collect(rsearch(seq_b, " ", indx))
    ce1 = collect(something(findprev(" ", seq_b, indx), 0:-1))
    seq_b = "$seq_b "
    string(seq_b[indx+1]) != " " ? sp = search(seq_b, " ", indx)[end] : sp = indx
    cs1 == 1 ? cs1 = 2 : cs1 = cs1
    if string(seq_b[cs1-1]) == " " || occursin(Regex("^$(seq_b[cs1-1])"), seq_b)
        cs1 == 2 ? cs1 = 1 : cs1 = cs1
        if length(ce1) > 0
            if seq_b[cs1] == seq_a[1] && seq_b[ce1[end]+1] == split(seq_a, " ")[end][1] && abs(sp-(indx+1)) <= 2 # || occursin(Regex("$(seq_b[indx])\$"), seq_b))
                if occursin(r"\s", seq_a)
                    return "match"
                else
                    if seq_a[end] == seq_b[indx]
                        return "match"
                    else
                        return "no-match"
                    end
                end
            else
                return "no-match"
            end
        else
            if seq_b[cs1] == seq_a[1] && seq_b[1] == split(seq_a, " ")[end][1] && seq_a[end] == seq_b[indx]
                return "match"
            else
                return "no-match"
            end
        end
    else
        return "no-match"
    end
end

# Retrieve accepted name
function accepted_name(ubt_id::String)
    url_string = solr_host * "/ubt_plant/select?q=ubt_id:$ubt_id%20AND%20type:sciname"
    #url_string = "http://bcbi.brown.edu/solr/solr/ubt_plant/select?q=ubt_id:$ubt_id%20AND%20type:sciname"
    ret_str  = HTTP.request("GET","$url_string";)
    json_str = JSON.parse(String(ret_str.body))
    #ubt_id   = json_str["response"]["docs"][i]["ubt_id"][1]
    name     = json_str["response"]["docs"][1]["term"][1]
    return name
end

# Resolve plant species names
function resolve_name(herb::String)
    ret_arr  = String[]
    name_arr = String[]
    n = trunc(Int, length(split(herb, " "))/2)
    n == 0 ? n = 1 : nothing
    for i in 1:n
        herb = replace(lowercase(herb), r"(\s[A-Za-z0-9]{2}\s|\s[A-Za-z0-9]{2}$|^[A-Za-z0-9]{2}\s)" => " ")
        herb = replace(herb, r"\s+" => " ")
        println("Herb:" , herb)
        url  = generate_sub_query(herb, name_arr)
        println("URL: ", url)
        res_str, num_found = norm_string(url,herb)
        println("Resolved name: ", res_str)
        num_found == 0 ? break : nothing
        if length(split(split(res_str, "\$")[1], " ")) >= 2
            split(res_str, "\$")[2] != "NA" ? name_subm = split(res_str, "\$")[1] : name_subm  = "No Matches"
            ubt_id     = String(split(res_str, "\$")[2])
            name_match = split(res_str, "\$")[3]
            match_typ  = uppercasefirst(split(res_str, "\$")[end])
            ubt_id != "NA" ? accept_name = accepted_name(ubt_id) : accept_name = "NA"
            match_typ == "Sciname" ? match_typ = "Scientific Name" : match_typ = match_typ
            push!(name_arr, split(res_str, "\$")[3])
            push!(ret_arr, """{"NameSubmitted":"$name_subm","NameMatched":"$name_match","TaxonomicStatus":"$match_typ","uBiotaID":"$ubt_id","AcceptedName":"$accept_name"}""")
            split(res_str, "\$")[2] == "NA" ? break : nothing
        else
            push!(ret_arr, """{"NameSubmitted":"No Matches","NameMatched":"NA","TaxonomicStatus":"NA","uBiotaID":"NA","AcceptedName":"NA"}""")
        end
    end
    ret_arr = unique(ret_arr)
    length(ret_arr) != 1 ? ret_arr = ret_arr[1:end-1] : nothing
    return "[$(join(ret_arr, ","))]"
end

# Take input as text for resolving names
function extract_plants(text, extractor)
    ret_arr = Any[]
    txt = replace(text, r"\n+" => " ")
    txt = replace(txt, r"\*\s" => "")
    txt = replace(txt, r"\%" => "")
    txt = strip(txt)
    println("text after stripping: ", txt)
    sent_arr = sentTokenization(txt, ntk)
    println("text after tokenization: ", sent_arr)

    for sent in sent_arr
        nphrs = textblob.TextBlob(sent, np_extractor=extractor)

        for phrs in nphrs["noun_phrases"]
            println("Noun phrase: ", phrs)
            try
                norm_str = resolve_name(phrs)
                println("Resolved name: ", norm_str)
                json_str = JSON.parse(String(norm_str))
                for i in 1:length(json_str)
                    if json_str[i]["NameSubmitted"] != "No Matches" && json_str[i]["NameSubmitted"] != "Unknown"
                        name_subm  = json_str[i]["NameSubmitted"]
                        name_match = json_str[i]["NameMatched"]
                        match_typ  = json_str[i]["TaxonomicStatus"]
                        ubt_id     = json_str[i]["uBiotaID"]
                        accept_name = json_str[i]["AcceptedName"]
                        #println("$(json_str[i]["NameSubmitted"])\$$(json_str[i]["NameMatched"])\$$(json_str[i]["AcceptedName"])")
                        push!(ret_arr, """{"NameSubmitted":"$name_subm","NameMatched":"$name_match","TaxonomicStatus":"$match_typ","uBiotaID":$ubt_id,"AcceptedName":"$accept_name"}""")
                        println("Ret Array: ", ret_arr)
                    end
                end
            catch err
                println(err)
                continue
            end
        end
    end
    return "[$(join(ret_arr, ","))]"
end

# function run_server()

#     query_dict = Dict()

#     HTTP.listen() do request::HTTP.Request

#         println("******************")
#         @show request
#         println("******************")

#         uri = parse(HTTP.URI, request.target)
#         query_dict = HTTP.queryparams(uri)

#         headers = Dict{AbstractString,AbstractString}(
#             "Server"            => "Julia/$VERSION",
#             "Content-Type"      => "text/html; charset=utf-8",
#             "Content-Language"  => "en",
#             "Date"              => Dates.format(now(Dates.UTC), Dates.RFC1123Format),
#             "Access-Control-Allow-Origin" => "*" )

#         return HTTP.Response(200, HTTP.Headers(collect(headers)), body = String(extract_plants(query_dict["text"])))
#     end

# end

# run_server()

#############
### TESTS ###
#############
#herb = "Strychnos nux-vomica"
herb = "This sentence contains Raulfia serpentina, Mangifera indica and Arabidoopsis thaliana of plantae, glycine and fabaceae family in it."
#herb = "This sentence has no analysis plantae ginkgo biloba." # names in it."
#herb = "Arabidopsis"
#herb = "This study focussed on the effect of increasing nitrogen (N) supply on root uptake and root-to-shoot translocation of zinc (Zn) as well as retranslocation of foliar-applied Zn in durum wheat (Triticum durum)."
#herb = "plantae"
#herb = "Apical buds of Norway spruce (Picea abies) trees at high and low elevation were heated during 2006 and 2007."
#herb_arr = ["Strychnos nux-vomica",
#    "This sentence contains Raulfia serpentina, Mangifera indica and Arabidoopsis thaliana of plantae and fabaceae family in it.",
#    "This sentence has no plantae names in it.",
#    "Arabidopsis",
#    "This sentence contains Raulfia serpentina, Mangifera indica and Arabidoopsis thaliana of plantae and fabaceae family in it.",
#    "This sentence contains Raulfia serpentina, Mangifera indica and Arabidoopsis thaliana of plantae and fabaceae family in it."]
#herb = "Arabidoopsis thaliane Pigment stripes associated arabidoopsis thaliane cycle with veins (venation) is a common flower colour pattern."
#herb = "Apical hydraulic conductivity (k) was estimated from anatomical data."
#herb = "Here we quantified intraspecific variation and covariation of leaf mass per area (LMA) and wood density (WD) in monospecific forests of the widespread tree species Nothofagus pumilio to determine its magnitude and whether it is related to environmental conditions and ontogeny."
#herb = "arabidopsis thaliana"
#herb = "Arabidopsis thaliana"
#herb = "This sentence has no plant names in it."
#startTime = time()
#for herb in herb_arr
#    println(resolve_name(herb))
#    println("\=\=\=\=\=")
#end
#endTime = time()
#println("########################")
#println("Time taken: $(round(endTime-startTime, 2)) sec.")
#println("########################")

# extractor = npx.ConllExtractor()
# text = herb
# pt_spec = extract_plants(text, extractor)
# println(pt_spec)

max_stars_repo_name: bcbi/SolrPlantAPI

max_issues_repo_name: bcbi/SolrPlantAPI

max_forks_repo_name: bcbi/SolrPlantAPI


text:
function get_erf_integral(JUPL_PATH)
    data = readdlm("$JUPL_PATH/assets/erf_integral.dat")
    erf_table = data[:,2]
end

function write_sim_parameters()
    filename = "output/parameters.csv"
    writedlm(filename, [@sprintf("%.5e", dtm) @sprintf("%d",nz) @sprintf("%.5e", dt_print) @sprintf("%.5e", tm_quiet) @sprintf("%.5e", L)], ",")
end

function open_files()
    	c1::Int = 0
    	c2::Int = 0
    	c3::Int = 0

    	for i = 1:(nspec+1)*3+2
    		if i==1
    			filename = "r.csv"
    		elseif i==(nspec+1)*3+2
    			filename = "efield.csv"
    		elseif (i>1) & (i<=(nspec+1)+1)
    			c1 += 1
    			filename = "vel" * string(c1) * ".csv"
    		elseif (i>(nspec+1)+1) & (i<=2*(nspec+1)+1)
                c2 += 1
    			filename = "den" * string(c2) * ".csv"
    		elseif (i>2*(nspec+1)+1) & (i<=3*(nspec+1)+1)
    			c3 += 1
    			filename = "temp" * string(c3) * ".csv"
    		end
            f = open("output/" * filename, "a")
    		@eval $(Symbol("f_$i")) = $f
    	end
end

function write_data(hydro)
	c1::Int = 0
	c2::Int = 0
	c3::Int = 0
    var = Array{Float64}(nz)

	for i = 1:(nspec+1)*3+2
		if i==1
			var = r
            filename = "r.csv"
		elseif i==(nspec+1)*3+2
			var = hydro.Efield
            filename = "efield.csv"
		elseif (i>1) & (i<=(nspec+1)+1)
			c1 += 1
			var = hydro.u[:,c1]
            filename = "vel" * string(c1) * ".csv"
		elseif (i>(nspec+1)+1) & (i<=2*(nspec+1)+1)
			c2 += 1
			if c2<=nspec
				var = 1.e-6 * hydro.rho[:,c2]	/ mi[c2]
			else #electrons
				var = 1.e-6 * hydro.rho[:,c2]	/ me
			end
            filename = "den" * string(c2) * ".csv"
		elseif (i>2*(nspec+1)+1) & (i<=3*(nspec+1)+1)
			c3 += 1
			var = hydro.T[:,c3] / qe_C
            filename = "temp" * string(c3) * ".csv"
		end
        f = @eval $(Symbol("f_$i"))
		for k = 1:nz
			writedlm(f, [@sprintf("%.24e", var[k])])
		end
	end
end

function close_files()
    for i = 1:(nspec+1)*3+2
        f = @eval $(Symbol("f_$i"))
        close(f)
    end
end

#---------------------------------------------------------
function write_all_data(hydro)
	c1::Int = 0
	c2::Int = 0
	c3::Int = 0
    var = Array{Float64}(nz)

    filename = "./output/U1D.csv"
    var = hydro.U1D
    writedlm(filename, var, ",")

    filename = "./output/U1D_p.csv"
    var = hydro.U1D_p
    writedlm(filename, var, ",")

    filename = "./output/U1D_c.csv"
    var = hydro.U1D_c
    writedlm(filename, var, ",")

    filename = "./output/F1D.csv"
    var = hydro.F1D
    writedlm(filename, var, ",")

    filename = "./output/G1D.csv"
    var = hydro.G1D
    writedlm(filename, var, ",")

    filename = "./output/R1D.csv"
    var = hydro.R1D
    writedlm(filename, var, ",")

    filename = "./output/C1D.csv"
    var = hydro.C1D
    writedlm(filename, var, ",")

    filename = "./output/T.csv"
    var = hydro.T
    writedlm(filename, var, ",")

end

max_stars_repo_name: cbellei/JUPL

max_issues_repo_name: cbellei/JUPL

max_forks_repo_name: cbellei/JUPL


text:
# test some utility functions specific to the python interface cosmo-python
using COSMO, Test, LinearAlgebra, Random, SparseArrays

rng = Random.MersenneTwister(41)

@testset "cosmo-python related utility function" begin

    # create a set of convex sets from a dictionary
    cone = Dict("f" => 2, "l" => 3, "q" => [3; 4], "s" => [3; 6; 10], "ep" => 2, "ed" => 1, "p" => [0.3; -0.4], "b" => 2)
    convex_sets = COSMO.convex_sets_from_dict(cone, -rand(rng, 2), rand(rng, 2))
    @test convex_sets[3] isa COSMO.SecondOrderCone
    @test convex_sets[3].dim == 3
    @test convex_sets[6] isa COSMO.PsdConeTriangle
    @test convex_sets[11] isa COSMO.PowerCone
    @test convex_sets[11].α == 0.3
    @test convex_sets[13] isa COSMO.Box

    # solve QP using set! function with sparse matrix input in (rowval, colptr, nzval)-format
    q = [1; 1.];
    P = sparse([4. 1; 1 2]);
    A1 = [1. 1; 1 0; 0 1];
    A = sparse([A1; -A1])

    lo = [1.; 0; 0];
    uo = [1; 0.7; 0.7];
    b = [uo; -lo]
    cone = Dict("l" => 6)
    m, n = size(A)
    model = COSMO.Model();
    COSMO.set!(model, P.rowval, P.colptr, P.nzval, q, A.rowval, A.colptr, A.nzval, b, cone, nothing, nothing, m, n);
    res = COSMO.optimize!(model);
end


text:
include("AlgAssOrd/AlgAssOrd.jl")
include("AlgAssOrd/Ideal.jl")
include("AlgAssOrd/FracIdeal.jl")
include("AlgAssOrd/Elem.jl")
include("AlgAssOrd/PicardGroup.jl")
include("AlgAssOrd/ResidueRingMultGrp.jl")
include("AlgAssOrd/LocallyFreeClassGroup.jl")


text:
function is_sum_of_pair(arr, num)
    n = length(arr)
    for i=1:n, j=1:n
        if i == j
            continue
        end
        if arr[i] + arr[j] == num
            return true
        end
    end
    false
end

#=
is_sum_of_pair(1:25, 49)
=#

function find_first_invalid(arr, preamble_size)
    for i=preamble_size+1:length(arr)
        check_arr = @view arr[i-preamble_size:i]
        if !is_sum_of_pair(check_arr, arr[i])
            return arr[i]
        end
    end
    nothing
end

function contiguous_set(arr, num)
    for i=1:length(arr)
        s = arr[i]
        for j=i+1:length(arr)
            s += arr[j]
            if s == num
                return @view arr[i:j]
            end
            if s > num
                break
            end
        end
    end
    nothing
end

function day09part1(arr)
    find_first_invalid(arr, 25)
end

function day09part2(arr)
    num = find_first_invalid(arr, 25)
    contig = contiguous_set(arr, num)
    minimum(contig) + maximum(contig)
end

function day09(file)
    s = read(file, String)
    arr = map(x->parse(Int, x), split(s, "\n", keepempty=false))
    (day09part1(arr), day09part2(arr))
end

#=
day09("input.txt")
=#


text:
@testset "Rules left of" begin
	grammar = @cgrammar begin
		Real = 1 | 2 | 3
		Real = Real + Real
		Real = Real * Real
	end

	iter = ContextExpressionIterator(grammar, 3, :Real)
	elems = collect(iter)

	focus_node = elems[50]

	path = Vector{Int}([1])
	@test rulesonleft(focus_node, path) == Set([4])

	path = Vector{Int}([2,1])
	@test rulesonleft(focus_node, path) == Set([4,3])

	path = Vector{Int}([2,2])
	@test rulesonleft(focus_node, path) == Set([4,3,1])

end


text:
@testset "Integers" begin
    let table = [0, 1, 2, 2, 3, 3, 4, 3, 4, 4,
                 5, 4, 5, 5, 5, 4, 5, 5, 6, 5,
                 6, 6, 6, 5, 6, 6, 6, 6, 7, 6,
                 7, 5, 6, 6, 7, 6, 7, 7, 7, 6,
                 7, 7, 7, 7, 7, 7, 8, 6, 7, 7]

        for (n, expected) in enumerate(table)
            try
                @test assembly(n) == expected
            catch
                @error "Test failed on" n
                rethrow
            end
        end
    end
end

@testset "Strings" begin
    try
        for n in 2:8
            elapsed = []
            for str in outer(join, 'a':'b', n)
                try
                    @test assembly(str) == stringpa(str)
                catch
                    @error "Test failed on" str
                    rethrow
                end
            end
        end
    finally
        rm.(["test.txt", "test_pathway.txt"])
    end
end


text:
###
  Name     : c5_31_CAPM.jl
  Book     : Hands-on Data Science with Anaconda )
  Publisher: Packt Publishing Ltd.
  Author   : Yuxing Yan and James Yan
  Date     : 5/30/2018
  email    : yany@canisius.edu
             paulyxy@hotmail.com

Note #1: go to http://finance.yahoo.com to dowload the last 5-year monthly
          data for IBM and S&P500. The ticker for S&P500 is ^GSPC
      #2: delete the first observation which contains 'null'
      #3: check if the last line is a blank line. If it is, delete it.
###

using DataFrames, CSV
using GLM, StatsModels
#
function f(x)
  p=x[Symbol("Adj Close")]
  n=length(p)
  a=p[2:n]./p[1:(n-1)]-1.0
  return convert(Array{Float64,1}, a)
end
#
x = CSV.read("c:/temp/ibmMonthly5years.txt")
stockRet=f(x)
y = CSV.read("c:/temp/sp500Monthly5years.txt")
mktRet=f(y)
n3=min(length(stockRet),length(mktRet))
data = DataFrame(X=mktRet[1:n3], Y=stockRet[1:n3])
OLS = glm(@formula(Y ~ X), data, Normal(), IdentityLink())


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule blis_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("blis")
JLLWrappers.@generate_main_file("blis", UUID("6136c539-28a5-5bf0-87cc-b183200dce32"))
end  # module blis_jll


text: abstract type MultiPtsDriver <: DriverModel{AccelSteeringDirection} end

max_stars_repo_name: sisl/UrbanDrivingSimulation

max_issues_repo_name: sisl/UrbanDrivingSimulation

max_forks_repo_name: sisl/UrbanDrivingSimulation


text:
__precompile__(true)

module NLSolversBase

using DiffEqDiffTools, ForwardDiff, Calculus, DiffResults
import Distributed: clear!
export AbstractObjective,
       NonDifferentiable,
       OnceDifferentiable,
       TwiceDifferentiable,
       TwiceDifferentiableHV,
       value,
       value!,
       value_gradient!,
       value_jacobian!,
       gradient,
       gradient!,
       jacobian,
       jacobian!,
       hessian,
       hessian!,
       value!!,
       value_gradient!!,
       value_jacobian!!,
       hessian!!,
       hv_product,
       hv_product!,
       only_fg!,
       only_fgh!,
       only_fj!,
       only_fg,
       only_fj,
       only_g_and_fg,
       only_j_and_fj,
       clear!,
       f_calls,
       g_calls,
       h_calls,
       hv_calls

export AbstractConstraints, OnceDifferentiableConstraints,
    TwiceDifferentiableConstraints, ConstraintBounds

function diffeqdiff_fdtype(autodiff)
    if autodiff == :finiteforward
        fdtype = Val{:forward}
    elseif autodiff == :finitecomplex
        fdtype = Val{:complex}
    elseif any(autodiff .== (:finite, :central, :finitecentral))
        fdtype = Val{:central}
    end
    fdtype
end

is_finitediff(autodiff) = autodiff ∈ (:central, :finite, :finiteforward, :finitecomplex)
is_forwarddiff(autodiff) = autodiff ∈ (:forward, :forwarddiff, true)

function x_of_nans(x)
    x_out = fill!(similar(x), NaN)
    x_out
end

include("objective_types/inplace_factory.jl")
include("objective_types/abstract.jl")
include("objective_types/nondifferentiable.jl")
include("objective_types/oncedifferentiable.jl")
include("objective_types/twicedifferentiable.jl")
include("objective_types/twicedifferentiablehv.jl")
include("objective_types/incomplete.jl")
include("objective_types/constraints.jl")
include("interface.jl")

end # module


text:
using PyCall

# Some utility functions to help with testing

export switch_endianness, get_statevector_using_picoquant, get_statevector_using_qiskit

"""
    function switch_endianness(vec)

Reverses the qubit ordering for a given state vector
"""
function switch_endianness(vec)
    n = convert(Int, log2(length(vec)))
    vec = reshape(vec, Tuple([2 for _ in 1:n]))
    vec = permutedims(vec, [i for i in n:-1:1])
    reshape(vec, 2^n)
end

"""
    function get_statevector_using_qiskit(circ; big_endian=true)

Uses qiskit to calculate the state vector for a given circuit
"""
function get_statevector_using_qiskit(circ; big_endian=true)
    qiskit = pyimport("qiskit")
    backend = qiskit.Aer.get_backend("statevector_simulator")
    job = qiskit.execute(circ, backend)
    result = job.result()
    vec = result.get_statevector()
    if !big_endian
        return switch_endianness(vec)
    else
        return vec
    end
end

"""
    function get_statevector_using_picoquant(circ; big_endian=false)

Uses picoquant to calculate statevector resulting from given circuit
"""
function get_statevector_using_picoquant(circ; big_endian=false)
    tn = convert_qiskit_circ_to_network(circ, InteractiveBackend{Array{ComplexF64}}(), decompose=false, transpile=false)
    qubits = circ.n_qubits
    add_input!(tn, "0"^qubits)
    full_wavefunction_contraction!(tn)
    result_tensor = load_tensor_data(tn, :result)
    dims = size(result_tensor)
    if !big_endian
        index_order = length(dims):-1:1
        result_tensor = permutedims(result_tensor, index_order)
    end
    reshape(result_tensor, prod(dims))
end


text:
@inline function logbtl(choices::Array{<:Int, 2},
                        goodness::Array{<:Real},
                        scale::Real)
    goodness_matrix = goodness[choices]
    scaled  = goodness_matrix / scale
    logpref = scaled[:,1] - logsumexp(scaled, dims=2)[:,1]
    sum(logpref)
end

@inline function logbtl_matrix(choices::Array{<:Int, 2},
                               goodness::Array{<:Real},
                               scale::Real)
    # goodness ∈ R^{npoints, ncandidates}
    # Compute the preferences of all entries
    goodness_matrix = goodness[choices]
    scaled  = goodness_matrix / scale
    Z       = logsumexp(scaled, dims=2)[:,1]
    @simd for i = 1:size(scaled, 1)
        @inbounds scaled[i,:] = scaled[i,:] .- Z[i]
    end
    return scaled
end

@inline function ∇logbtl(logbtl_matrix::Array{<:Real, 2},
                         latent::Array{<:Real, 1},
                         choices::Array{<:Int, 2},
                         scale::Real)
    pref         = exp.(logbtl_matrix)
    partial_grad = zeros(size(logbtl_matrix))

    @inbounds begin
        # first column derivative
        pick  = pref[:,1]
        ∇pick = (1 .- pick) / scale
        partial_grad[:,1] = ∇pick

        # other column derivative
        comps  = pref[:,2:end]
        ∇comp  = comps / -scale
        partial_grad[:,2:end] = ∇comp
    end

    partial_grad = reshape(partial_grad, :)
    choices      = reshape(choices, :)
    full_grad    = zeros(length(latent))
    @simd for i = 1:length(partial_grad)
        @inbounds full_grad[choices[i]] += partial_grad[i]
    end
    return full_grad
end

function ∇²logbtl(logbtl_matrix::Array{<:Real, 2},
                  latent::Array{<:Real},
                  choices::Array{<:Int, 2},
                  scale::Real)
    pref    = exp.(logbtl_matrix)

    # hessian
    nlatent = length(latent)
    ndata   = size(logbtl_matrix, 1)
    ncand   = size(logbtl_matrix, 2)
    result  = zeros(nlatent, nlatent)

    @inbounds for entry_j = 1:ncand
        @inbounds for entry_i = 1:ncand
            @simd for block = 1:ndata
                H_i_j = pref[block, entry_j] / (scale*scale)
                if(entry_i == entry_j)
                    H_i_j *= pref[block, entry_j] - 1
                else
                    H_i_j *= pref[block, entry_i]
                end
                latent_i = choices[block, entry_i]
                latent_j = choices[block, entry_j]
                result[latent_i, latent_j] += H_i_j
            end
        end
    end
    return Symmetric(result)
end

@inline function approx_marginallike(a::Vector,
                                     B::LU,
                                     choices::Array{<:Int, 2},
                                     latent::Array{<:Real},
                                     scale::Real)
    # GPML 3.32
    loglike = logbtl(choices, latent, scale)
    t1 = loglike
    t2 = dot(latent, a) / -2
    t3 = logdet(B) / -2
    t1 + t2 + t3
end

@inline function approx_marginallike(K::Matrix,
                                     choices::Array{<:Int, 2},
                                     latent::Array{<:Real},
                                     scale::Real)
    # GPML 3.32
    btl_mat = logbtl_matrix(choices, latent, scale)
    ∇ll     = ∇logbtl(btl_mat, latent, choices, scale)
    W       = -∇²logbtl(btl_mat, latent, choices, scale)

    f   = reshape(latent, :)
    WK  = W*K
    b   = W*f + ∇ll
    B   = I + WK
    Blu = lu(B)
    a   = (b - Blu \ (WK*b))

    loglike = sum(btl_mat[:,1])
    t1 = loglike
    t2 = dot(f, a) / -2
    t3 = logdet(Blu) / -2

    # GPML 3.32
    t1 + t2 + t3
end

@inline function logjointlike(K::PDMats.PDMat,
                              choices::Array{<:Int, 2},
                              latent::Array{<:Real},
                              scale::Real)
    # GPML 3.12
    t3 = logdet(K) / -2
    t4 = size(K, 1) * log(2*π) / -2
    loglike = logbtl(choices, latent, scale)

    t1 = loglike
    t2 = PDMats.invquad(K, latent) / -2

    t1 + t2 + t3 + t4
end


text:
# Important packages
println("== Compiling useful packages...")
using JLD

# Importing the words file
println("== Importing dictionary...")
f = open("preprocessed_data/words.txt")
raw_words = readlines(f)
close(f)

# Preprocessing words
println("== Preprocessing dictionary...")
words = Array(AbstractString, 0)
for raw_word in raw_words
	push!(words,strip(raw_word, ['\n']))
end

# Importing the sents file
println("== Importing sentences...")
f = open("preprocessed_data/sents.txt")
sents = readlines(f)
close(f)

# Checking array sizes
println("- Dictionary size: ", length(words))
println("- Number of sentences: ", length(sents))

# Initializing the count array
println("== Initializing the count array...")
count = zeros(Int32, length(words), length(words))

# Counting word matches (symmetric)
println("== Counting...")
w = 0
for sent in sents
	for i = 1:length(words)
		if (contains(sent, words[i]))
			for j = i:length(words)
				if (contains(sent, words[j]))
					count[i,j] += 1
					count[j,i] += 1
				end
			end
		end
	end
	w+=1
	println(w, "/", length(sents), " sents processed")
end

# Saving the array
println("== Exporting count array...")
save("backup_matrices/counts.jld", "counts", count)


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

mutable struct TroubleshootingRecommendedActions <: SwaggerModel
    actionId::Any # spec type: Union{ Nothing, String } # spec name: actionId
    actionText::Any # spec type: Union{ Nothing, String } # spec name: actionText
    actionUri::Any # spec type: Union{ Nothing, String } # spec name: actionUri
    actionUriText::Any # spec type: Union{ Nothing, String } # spec name: actionUriText

    function TroubleshootingRecommendedActions(;actionId=nothing, actionText=nothing, actionUri=nothing, actionUriText=nothing)
        o = new()
        validate_property(TroubleshootingRecommendedActions, Symbol("actionId"), actionId)
        setfield!(o, Symbol("actionId"), actionId)
        validate_property(TroubleshootingRecommendedActions, Symbol("actionText"), actionText)
        setfield!(o, Symbol("actionText"), actionText)
        validate_property(TroubleshootingRecommendedActions, Symbol("actionUri"), actionUri)
        setfield!(o, Symbol("actionUri"), actionUri)
        validate_property(TroubleshootingRecommendedActions, Symbol("actionUriText"), actionUriText)
        setfield!(o, Symbol("actionUriText"), actionUriText)
        o
    end
end # type TroubleshootingRecommendedActions

const _property_map_TroubleshootingRecommendedActions = Dict{Symbol,Symbol}(Symbol("actionId")=>Symbol("actionId"), Symbol("actionText")=>Symbol("actionText"), Symbol("actionUri")=>Symbol("actionUri"), Symbol("actionUriText")=>Symbol("actionUriText"))
const _property_types_TroubleshootingRecommendedActions = Dict{Symbol,String}(Symbol("actionId")=>"String", Symbol("actionText")=>"String", Symbol("actionUri")=>"String", Symbol("actionUriText")=>"String")
Base.propertynames(::Type{ TroubleshootingRecommendedActions }) = collect(keys(_property_map_TroubleshootingRecommendedActions))
Swagger.property_type(::Type{ TroubleshootingRecommendedActions }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_TroubleshootingRecommendedActions[name]))}
Swagger.field_name(::Type{ TroubleshootingRecommendedActions }, property_name::Symbol) =  _property_map_TroubleshootingRecommendedActions[property_name]

function check_required(o::TroubleshootingRecommendedActions)
    true
end

function validate_property(::Type{ TroubleshootingRecommendedActions }, name::Symbol, val)
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 ProtocolConfiguration <: SwaggerModel
    HTTPConfiguration::Any # spec type: Union{ Nothing, HTTPConfiguration } # spec name: HTTPConfiguration

    function ProtocolConfiguration(;HTTPConfiguration=nothing)
        o = new()
        validate_property(ProtocolConfiguration, Symbol("HTTPConfiguration"), HTTPConfiguration)
        setfield!(o, Symbol("HTTPConfiguration"), HTTPConfiguration)
        o
    end
end # type ProtocolConfiguration

const _property_map_ProtocolConfiguration = Dict{Symbol,Symbol}(Symbol("HTTPConfiguration")=>Symbol("HTTPConfiguration"))
const _property_types_ProtocolConfiguration = Dict{Symbol,String}(Symbol("HTTPConfiguration")=>"HTTPConfiguration")
Base.propertynames(::Type{ ProtocolConfiguration }) = collect(keys(_property_map_ProtocolConfiguration))
Swagger.property_type(::Type{ ProtocolConfiguration }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ProtocolConfiguration[name]))}
Swagger.field_name(::Type{ ProtocolConfiguration }, property_name::Symbol) =  _property_map_ProtocolConfiguration[property_name]

function check_required(o::ProtocolConfiguration)
    true
end

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


text:
using Plots
using Random
Random.seed!(9)

function square(c,side)
  x = [ c[1]-side/2, c[1]+side/2, c[1]+side/2, c[1]-side/2, c[1]-side/2]
  y = [ c[2]-side/2, c[2]-side/2, c[2]+side/2, c[2]+side/2, c[2]-side/2]
  return x, y
end

plot()

x,y=square([5,5],2)
plot!(x,y,seriestype=[:shape],
      linewidth=2,fillalpha=0.05,color="green",label="")

x,y=square([5,5],6)
plot!(x,y,seriestype=[:shape],
      linewidth=2,fillalpha=0.05,color="orange",label="")

lines = collect(2:2:8)
vline!(lines,color="gray",label="",style=:dash)
hline!(lines,color="gray",label="",style=:dash)

px = [ 0.1 + 9.8*rand() for i in 1:100 ]
py = [ 0.1 + 9.8*rand() for i in 1:100 ]
scatter!(px,py,label="",alpha=0.05,color="blue")

using LaTeXStrings

fontsize=8
annotate!(3,3,text("(i-1,j-1)",fontsize,:Courier))
annotate!(5,3,text("(i-1,j)",fontsize,:Courier))
annotate!(7,3,text("(i-1,j+1)",fontsize,:Courier))

annotate!(3,5,text("(i,j-1)",fontsize,:Courier))
annotate!(5,5,text("(i,j)",fontsize,:Courier))
annotate!(7,5,text("(i,j+1)",fontsize,:Courier))

annotate!(3,7,text("(i+1,j-1)",fontsize,:Courier))
annotate!(5,7,text("(i+1,j)",fontsize,:Courier))
annotate!(7,7,text("(i+1,j+1)",fontsize,:Courier))

plot!(size=(400,400),
      xlim=(1.3,8.7),xticks=:none,
      ylim=(1.3,8.7),yticks=:none,
      framestyle=:box,
      xlabel="x",ylabel="y",grid=false)
savefig("lcells2.pdf")


text:
using Statistics

"""
```
mlr_beta(y::Array{T}; n::Int64=10, x::Array{T}=collect(1.0:n))::Matrix{Float64} where {T<:Real}
```

Moving linear regression intercept (column 1) and slope (column 2)
"""
function mlr_beta(y::Array{T}; n::Int64=10, x::Array{T}=collect(1.0:n))::Matrix{Float64} where {T<:Real}
    @assert n<length(y) && n>0 "Argument n out of bounds."
    @assert size(y,2) == 1
    @assert size(x,1) == n || size(x,1) == size(y,1)
    const_x = size(x,1) == n
    out = zeros(T, (length(y),2))
    out[1:n-1,:] .= NaN
    xbar = mean(x)
    ybar = runmean(y, n=n, cumulative=false)
    @inbounds for i = n:length(y)
        yi = y[i-n+1:i]
        xi = const_x ? x : x[i-n+1:i]
        out[i,2] = cov(xi,yi) / var(xi)
        out[i,1] = ybar[i] - out[i,2]*xbar
    end
    return out
end

"""
```
mlr_slope(y::Array{T}; n::Int64=10, x::Array{T}=collect(1.0:n))::Array{Float64} where {T<:Real}
```

Moving linear regression slope
"""
function mlr_slope(y::Array{T}; n::Int64=10, x::Array{T}=collect(1.0:n))::Array{Float64} where {T<:Real}
    @assert n<length(y) && n>0 "Argument n out of bounds."
    @assert size(y,2) == 1
    @assert size(x,1) == n || size(x,1) == size(y,1)
    const_x = size(x,1) == n
    out = zeros(size(y))
    out[1:n-1] .= NaN
    @inbounds for i = n:length(y)
        yi = y[i-n+1:i]
        xi = const_x ? x : x[i-n+1:i]
        out[i] = cov(xi,yi) / var(xi)
    end
    return out
end

"""
```
mlr_intercept(y::Array{T}; n::Int64=10, x::Array{T}=collect(1.0:n))::Array{Float64} where {T<:Real}
```

Moving linear regression y-intercept
"""
function mlr_intercept(y::Array{T}; n::Int64=10, x::Array{T}=collect(1.0:n))::Array{Float64} where {T<:Real}
    @assert n<length(y) && n>0 "Argument n out of bounds."
    @assert size(y,2) == 1
    @assert size(x,1) == n || size(x,1) == size(y,1)
    const_x = size(x,1) == n
    out = zeros(size(y))
    out[1:n-1] .= NaN
    xbar = mean(x)
    ybar = runmean(y, n=n, cumulative=false)
    @inbounds for i = n:length(y)
        yi = y[i-n+1:i]
        xi = const_x ? x : x[i-n+1:i]
        out[i] = ybar[i] - xbar*(cov(xi,yi)/var(xi))
    end
    return out
end

"""
```
mlr(y::Array{T}; n::Int64=10)::Array{Float64} where {T<:Real}
```

Moving linear regression predictions
"""
function mlr(y::Array{T}; n::Int64=10)::Array{Float64} where {T<:Real}
    b = mlr_beta(y, n=n)
    return b[:,1] + b[:,2]*float(n)
end

"""
```
mlr_se(y::Array{T}; n::Int64=10)::Array{Float64} where {T<:Real}
```

Moving linear regression standard errors
"""
function mlr_se(y::Array{T}; n::Int64=10)::Array{Float64} where {T<:Real}
    yhat = mlr(y, n=n)
    r = zeros(T, n)
    out = zeros(size(y))
    out[1:n-1] .= NaN
    nf = float(n)
    @inbounds for i = n:length(y)
        r = y[i-n+1:i] .- yhat[i]
        out[i] = sqrt(sum(r.^2)/nf)
    end
    return out
end

"""
```
mlr_ub(y::Array{T}; n::Int64=10, se::T=2.0)::Array{Float64} where {T<:Real}
```

Moving linear regression upper bound
"""
function mlr_ub(y::Array{T}; n::Int64=10, se::T=2.0)::Array{Float64} where {T<:Real}
    return y + se*mlr_se(y, n=n)
end

"""
```
mlr_lb(y::Array{T}; n::Int64=10, se::T=2.0)::Array{Float64} where {T<:Real}
```

Moving linear regression lower bound
"""
function mlr_lb(y::Array{T}; n::Int64=10, se::T=2.0)::Array{Float64} where {T<:Real}
    return y - se*mlr_se(y, n=n)
end

"""
```
mlr_bands(y::Array{T}; n::Int64=10, se::T=2.0)::Matrix{Float64} where {T<:Real}
```

Moving linear regression bands

*Output:*

Column 1: Lower bound

Column 2: Regression estimate

Column 3: Upper bound
"""
function mlr_bands(y::Array{T}; n::Int64=10, se::T=2.0)::Matrix{Float64} where {T<:Real}
    out = zeros(T, (length(y),3))
    out[1:n-1,:] .= NaN
    out[:,2] = mlr(y, n=n)
    out[:,1] = mlr_lb(y, n=n, se=se)
    out[:,3] = mlr_ub(y, n=n, se=se)
    return out
end

"""
```
mlr_rsq(y::Array{T}; n::Int64=10, adjusted::Bool=false)::Array{Float64} where {T<:Real}
```

Moving linear regression R-squared or adjusted R-squared
"""
function mlr_rsq(y::Array{T}; n::Int64=10, adjusted::Bool=false)::Array{Float64} where {T<:Real}
    yhat = mlr(y, n=n)
    rsq = runcor(y, yhat, n=n, cumulative=false) .^ 2.0
    if adjusted
        return rsq .- (1.0.-rsq)*(1.0/(float(n).-2.0))
    else
        return rsq
    end
end


text:
using DelimitedFiles
using Arpack
import Statistics: mean, median

writecsv(filename::AbstractString, data) = writedlm(filename, data, ',')
readcsv(filename::AbstractString) = readdlm(filename, ',')
readcsv(filename::AbstractString, ::Type{T}) where {T} = readdlm(filename, ',', T)

input = ARGS[1]
output1 = ARGS[2]
output2 = ARGS[3]
dim = ARGS[4]

dim = parse(Int64, dim)
println(input)
println(output1)
println(output2)
println(dim)

# Data loading
out = readcsv(input)
# Library size
libsize = sum(out, dims=1)
# CPMED
out = median(libsize) .* (out ./ libsize)
# Log transformation
out = log10.(out .+ 1)
# Centering
out = out .- vec(mean(out, dims=2))
# PCA
res = svds(out, nsv=dim)[1]
# Eigenvectors
writecsv(output1, res.Vt')
# Eigenvalues
writecsv(output2, res.S.*res.S/size(res.Vt)[2])

max_stars_repo_name: rikenbit/onlinePCA-experiments

max_issues_repo_name: rikenbit/onlinePCA-experiments

max_forks_repo_name: rikenbit/onlinePCA-experiments


text:
module PokerBots

import Logging
import TexasHoldem
const TH = TexasHoldem

include("bots.jl")

function battle!(n_tournaments, bots...)
    n_players = length(bots)

    initial_players = ntuple(n_players) do i
        Player(bots[i], i)
    end

    FT = Float64

    stats = Dict{Symbol,Any}()

    stats[:player_bank_rolls] = map(initial_players) do player
        (Vector{FT}(), TH.seat_number(player), TH.life_form(player))
    end

    stats[:tournament_winnings] = zeros(length(initial_players))

    # Assumes all BRs start out equal
    tournament_winnings = sum(map(initial_players) do player
        TH.bank_roll(player)
    end) - TH.bank_roll(first(initial_players))

    n_games = 0
    local winners
    for i in 1:n_tournaments
        players = deepcopy(initial_players)

        # @info "******************************** Playing tournament!"
        # TODO: ensure button moves continuously through multiple tournaments
        game = TH.Game(players)
        table = game.table
        players = TH.players_at_table(table)
        while length(players) > 1

            Logging.with_logger(Logging.NullLogger()) do
                winners = TH.play!(game)
            end

            n_games+=1
            for (j, player) in enumerate(players)
                @assert TH.seat_number(player) == j
                sn = TH.seat_number(player)
                br = TH.bank_roll(player)
                ptw = stats[:tournament_winnings][sn]
                push!(stats[:player_bank_rolls][sn][1], br+ptw)
            end

            n_players_remaining = count(map(x->!(TH.bank_roll(x) ≈ 0), players))
            if n_players_remaining ≤ 1
                if !(length(winners.players) == 1) # 1 victor per tournament
                    # @show winners.players
                    # @show n_players_remaining
                    if !all(TH.seat_number.(winners.players) .== Ref(TH.seat_number(first(winners.players))))
                        @show winners.players
                        @show TH.seat_number.(winners.players)
                        error("Nooooo")
                    end
                    # TODO: fix this bug in TH.
                    winners.players = (winners.players[1],)
                end
                for winner in winners.players
                    sn = TH.seat_number(winner)
                    stats[:tournament_winnings][sn]+=tournament_winnings
                end
                # println("Victor emerges!")
                break
            end
            TH.reset_game!(game)
            TH.move_buttons!(game)
        end
        # @info "******************************** Finished tournament!"
        if n_tournaments > 1000
            if mod(i, 100) == 0
                @info "progress = $(i/n_tournaments*100)%"
            end
        else
            @info "progress = $(i/n_tournaments*100)%"
        end
    end
    return stats, n_games
end

end # module


text:
# # Linear Program
#
# We want to solve the following linear program with decision variable `x`:
# $$
# \begin{array}{ll} \text{minimize} &  c^\top x\\
# \text{subject to} &  A x \leq b \\
#                   &  x \geq 1 \\
#                   &  x_2 \geq 5 \\
#                   &  x_1 + x_3 \geq 4.
# \end{array}
# $$
# The problem can be solved with `COSMO` in the following way:

using COSMO, LinearAlgebra, SparseArrays, Test
#-
##Define problem data:
c = [1; 2; 3; 4.];
A = Matrix(1.0I, 4, 4);
b = [10.; 10; 10; 10];
n = 4;

# Create the constraints $Ax + b \in \mathcal{K}$:

## Ax <= b
c1 = COSMO.Constraint(-A, b, COSMO.Nonnegatives);
## x >= 1
c2 = COSMO.Constraint(Matrix(1.0I, n, n), -ones(n), COSMO.Nonnegatives);
## x2 >= 5
c3 = COSMO.Constraint(1, -5, COSMO.Nonnegatives, n, 2:2);
## x1 + x3 >= 4
c4 = COSMO.Constraint([1 0 1 0], -4, COSMO.Nonnegatives);

# Define matrix $P$ and vector $q$ for the objective function:
P = spzeros(4, 4);
q = c;

# Create, assemble and solve the model:
settings = COSMO.Settings(verbose=true, eps_abs = 1e-4, eps_rel = 1e-5);
model = COSMO.Model();
assemble!(model, P, q, [c1; c2; c3; c4], settings = settings);
res = COSMO.optimize!(model)

# Compare the result to the known solution:
@test isapprox(res.x[1:4], [3; 5; 1; 1], atol=1e-2, norm = (x -> norm(x, Inf)))
#-
@test isapprox(res.obj_val, 20.0, atol=1e-2)


text:
using PkgMassInstaller
using Test

function try_remove_install_pkg(pkgname)
    for depot in DEPOT_PATH
        pkg_dir = joinpath(depot, "packages", pkgname)
        rm(pkg_dir; recursive = true, force = true)
    end
end

function count_installed_versions(pkgname)
    count = 0
    for depot in DEPOT_PATH
        pkg_dir = joinpath(depot, "packages", pkgname)
        !isdir(pkg_dir) && continue
        for dir in readdir(pkg_dir)
            ver_dir = joinpath(pkg_dir, dir)
            isdir(ver_dir) && (count += 1)
        end
    end
    count
end

@testset "PkgMassInstaller.jl" begin

    proj_dir = @__DIR__

    pkgname = "Example"
    try_remove_install_pkg(pkgname)

    println("\n" ^ 3)
    @info("Real Run, -deep = 3")
    deep = 3
    dry_run = false
    PkgMassInstaller.install_projs(proj_dir, deep, dry_run)
    count = count_installed_versions(pkgname)
    @test count >= deep

end


text:
# ==============================================================================
# Xavier Gandibleux - November 2021
#   Implemented in Julia 1.6

# ==============================================================================

println("""\nProject MOMH 2021" --------------------------------\n""")

const verbose = false
const graphic = true

include("parserMomkpZL.jl")
include("parserMomkpPG.jl")
include("vOptMomkp.jl")
include("displayGraphic.jl")
include("VNS.jl")
include("MOKP.jl")

# ==============================================================================
# Datastructure of a generic bi-objective 0/1 IP where all coefficients are integer
#
#   Max  sum{j=1,...,n} C(k,j) x(j)                k=1,...,p
#   st   sum{j=1,...,n} A(i,j) x(j) <= b_{i}       i=1,...,m
#                       x(j) = 0 or 1

# ==============================================================================

# Example ----------------------------------------------------------------------
fname = "instancesZL/knapsack.250.3"

# Read and load an instance of MO-01MKP from the collection of E. Zitzler / M. Laumanns
momkpZL = readInstanceMOMKPformatZL(verbose, fname)

# Reduce the MO-MKP instance to the two first objectives and store it as a generic Bi-01IP
dat1 = _bi01IP(momkpZL.P[1:2, :], momkpZL.W, momkpZL.ω)

# Example ----------------------------------------------------------------------
#fname = "instancesPG/set1/ZL28.DAT"
#fname = "instancesPG/set2/kp28W-Perm.DAT"
#fname = "instancesPG/set3/W7BI-rnd1-1800.DAT"

# Read and load an instance of Bi-01BKP from the collection of O. Perederieieva / X. Gandibleux
#momkpPG = readInstanceMOMKPformatPG(verbose, fname)

# Store it as a generic bi-01IP
#dat2 = _bi01IP(momkpPG.P, momkpPG.W, momkpPG.ω)

# Example ----------------------------------------------------------------------
#solverSelected = GLPK.Optimizer
#YN, XE = vSolveBi01IP(solverSelected, dat1.C, dat1.A, dat1.b)
#YN, XE = vSolveBi01IP(solverSelected, dat2.C, dat2.A, dat2.b)

# ---- Displaying the results (XE and YN)
# for n = 1:length(YN)
#     X = value.(XE, n)
#     print(findall(elt -> elt ≈ 1, X))
#     println("| z = ", YN[n])
# end

# for n = 1:length(YN)
#     println(YN[n])
# end

function TOOL_read_YN(prob::_MOMKP)

    f = open(string("YN/", prob.name, ".DAT"), "r")
    lines = readlines(f)
    close(f)

    len = length(lines) - 1

    YN = Vector{Vector{Int}}(undef, len)

    for i = 1:len
        YN[i] = parse.(Int64, split(lines[i+1]))
    end

    return YN
end
# Example ----------------------------------------------------------------------

function pointsExtremes(prob::_bi01IP, obj::Int)
    model = Model(with_optimizer(GLPK.Optimizer))

    @variable(model, x[1:size(prob.C, 2)], Bin)

    @constraint(model, [i = 1:length(prob.b)], sum((prob.A[i, j] * x[j]) for j = 1:size(prob.C, 2)) <= prob.b[i])

    if obj == 1
        @objective(model, Max, sum((prob.C[1, i]) * x[i] for i = 1:size(prob.C, 2)))
    else
        @objective(model, Max, sum((prob.C[2, i]) * x[i] for i = 1:size(prob.C, 2)))
    end
    optimize!(model)
    return value.(x)
end

pop = initPop(100, dat1)
x1 = solution([], pointsExtremes(dat1, 1), [])
x2 = solution([], pointsExtremes(dat1, 2), [])
verification(dat1, x1)
verification(dat1, x2)
push!(pop, x1)
push!(pop, x2)
sols = GVNS(pop, 2, 30, 1, dat1)
YN_VNS = map(x -> -x.val_objectif, sols)
displayGraphics(fname, TOOL_read_YN(momkpZL), "black")
displayGraphics(fname, YN_VNS, "green")


text:
using OrthogonalPolynomialsQuasi, BandedMatrices, ArrayLayouts, Test
import OrthogonalPolynomialsQuasi: recurrencecoefficients, PaddedLayout

@testset "Lanczos" begin
    @testset "Legendre" begin
        P = Legendre();
        w = P * [1; zeros(∞)];
        Q = LanczosPolynomial(w);
        @test Q.data.W[1:10,1:10] isa BandedMatrix

        Q̃ = Normalized(P);
        A,B,C = recurrencecoefficients(Q);
        Ã,B̃,C̃ = recurrencecoefficients(Q̃);
        @test @inferred(A[1:10]) ≈ Ã[1:10] ≈ [A[k] for k=1:10]
        @test @inferred(B[1:10]) ≈ B̃[1:10] ≈ [B[k] for k=1:10]
        @test @inferred(C[2:10]) ≈ C̃[2:10] ≈ [C[k] for k=2:10]

        @test A[1:10] isa Vector{Float64}
        @test B[1:10] isa Vector{Float64}
        @test C[1:10] isa Vector{Float64}

        @test Q[0.1,1:10] ≈ Q̃[0.1,1:10]

        R = P \ Q
        @test R[1:10,1:10] ≈ (P \ Q̃)[1:10,1:10]
        @test (Q'Q)[1:10,1:10] ≈ I

        # Q'Q == I => Q*sqrt(M) = P

        x = axes(P,1)
        X = Q' * (x .* Q)
        X̃ = Q̃' * (x .* Q̃)
        @test X[1:10,1:10] ≈ X̃[1:10,1:10]
    end

    @testset "other weight" begin
        P = Normalized(Legendre())
        x = axes(P,1)

        Q = LanczosPolynomial(exp.(x))
        R = P \ Q
        @test MemoryLayout(R[:,2]) isa PaddedLayout

        A,B,C = recurrencecoefficients(Q)
        @test A[1] ≈ 1.903680130866564 # emperical from Mathematica
        @test B[1] ≈ -0.5959190532652192
        @test A[2] ≈ 1.9150612001588696
        @test B[2] ≈ 0.0845629033308663
        @test C[2] ≈ 1.005978456731134

        @test Q[0.1,1] ≈ (Q * [1; zeros(∞)])[0.1] ≈ P[0.1,:]'*R[:,1] ≈ 0.6522722316024658
        @test Q[0.1,2] ≈ (Q * [0; 1; zeros(∞)])[0.1] ≈ P[0.1,:]'*R[:,2] ≈ -0.26452968200597243
        @test Q[0.1,3] ≈ (Q * [zeros(2); 1; zeros(∞)])[0.1] ≈ P[0.1,:]'*R[:,3] ≈ -0.7292002638736375
        @test Q[0.1,5] ≈ 0.7576999562707534 # emperical

        Q = LanczosPolynomial(  1 ./ (2 .+ x));
        R = P \ Q
        @test norm(R[1,3:10]) ≤ 1E-14

        Q = LanczosPolynomial(  1 ./ (2 .+ x).^2);
        R = P \ Q
        @test norm(R[1,4:10]) ≤ 1E-14

        # polys
        Q = LanczosPolynomial( 2 .+ x);
        R = P \ Q;
        @test norm(inv(R)[1,3:10]) ≤ 1E-14

        w = P * (P \ (1 .+ x))
        Q = LanczosPolynomial(w)
        @test Q[0.1,5] ≈ Normalized(Jacobi(0,1))[0.1,5] ≈ 0.742799258138176

        Q = LanczosPolynomial( 1 .+ x.^2);
        R = P \ Q;
        @test norm(inv(R)[1,4:10]) ≤ 1E-14
    end

    @testset "Expansion" begin
        P = Legendre();
        x = axes(P,1)
        w = P * [1; zeros(∞)];
        Q = LanczosPolynomial(w);
        R = Normalized(P) \ Q
        @test R * [1; 2; 3; zeros(∞)] ≈ [R[1:3,1:3] * [1,2,3]; zeros(∞)]
        @test R \ [1; 2; 3; zeros(∞)] ≈ [1; 2; 3; zeros(∞)]
        @test (Q * (Q \ (1 .- x.^2)))[0.1] ≈ (1-0.1^2)
    end

    @testset "Jacobi via Lanczos" begin
        P = Legendre(); x = axes(P,1)
        w = P * (P \ (1 .- x.^2))
        Q = LanczosPolynomial(w)
        A,B,C = recurrencecoefficients(Q)

        @test @inferred(Q[0.1,1]) ≈ sqrt(3)/sqrt(4)
        @test Q[0.1,2] ≈ 2*0.1 * sqrt(15)/sqrt(16)
    end

    @testset "Singularity" begin
        T = Chebyshev(); wT = WeightedChebyshev()
        x = axes(T,1)

        w = wT * [1; zeros(∞)];
        Q = LanczosPolynomial(w)
        @test Q[0.1,1:10] ≈ Normalized(T)[0.1,1:10]

        @test (Q'*(ChebyshevTWeight() .* Q))[1:10,1:10] ≈ I
        @test (Q'*(w .* Q))[1:10,1:10] ≈ I
    end

    @testset "BigFloat" begin
        P = Legendre{BigFloat}()
        x = axes(P,1)
        w = P * (P \ exp.(x))
        W = P \ (w .* P)
        v = [[1,2,3]; zeros(BigFloat,∞)];
        Q = LanczosPolynomial(w)

        x̃ = BigFloat(1)/10
        @test Q[x̃,1] ≈ 0.652272231602465791008015756161075576539994569266308567422055126278763683344388252
        @test Q[x̃,2] ≈ -0.26452968200597244253463861322599173806126678155361307561211048667577270734771616
        @test Q[x̃,3] ≈ -0.72920026387366053084159259908849371062183891778315602761397748592062615496583854

        X = Q \ (x .* Q)
        # empirical test
        @test X[5,5] ≈ -0.001489975039238321407179828331585356464766466154894764141171294038822525312179884

        @test (Q*[1; 2; 3; zeros(BigFloat,∞)])[0.1] ≈ -2.0643879240304606865860392675563890314480557471903856666440983048346601962485597
        @test 0.1*(Q*[1; 2; 3; zeros(BigFloat,∞)])[0.1] ≈ (Q * (X * [1; 2; 3; zeros(BigFloat,∞)]))[0.1]
    end

    @testset "Mixed Jacobi" begin
        P = Jacobi(1/2,0)
        x = axes(P,1)

        w = @. sqrt(1-x)
        Q = LanczosPolynomial(w, P)
        @test Q[0.1,1:10] ≈ Normalized(P)[0.1,1:10]

        w = @. exp(x) * sqrt(1-x)
        Q = LanczosPolynomial(w, P)
        # emperical from Julia
        @test Q[0.1,10] ≈ 0.5947384257847858
    end

    @testset "Mapped" begin
        x = Inclusion(0..1)
        w = @. sqrt(1 - x^2)
        Q = LanczosPolynomial(w, jacobi(1/2,0,0..1))
        # emperical from Julia
        @test Q[0.1,10] ≈ -0.936819626414421

        @test (Q * (Q \ exp.(x)))[0.1] ≈ exp(0.1)
        @test (Q[:,1:20] \ exp.(x)) ≈ (Q \ exp.(x))[1:20]

        @testset "Mapped Conversion" begin
            P₊ = jacobi(0,1/2,0..1)
            x = axes(P₊,1)
            y = @.(sqrt(x)*sqrt(2-x))
            U = LanczosPolynomial(y, P₊)
            @test P₊ ≠ U
            R = P₊ \ U
            @test U[0.1,5] ≈ (P₊ * R * [zeros(4); 1; zeros(∞)])[0.1]
        end
    end

    @testset "orthogonality (#68)" begin
        α = -0.5
        β = -0.5
        w = JacobiWeight(α,β)
        p = LanczosPolynomial(w)
        x = axes(w,1)
        b = 2
        wP = WeightedJacobi(α+1,β+1)
        ϕw = wP * (Jacobi(α+1,β+1) \ (b .- x))
        pϕ = LanczosPolynomial(ϕw)
        @test (pϕ.P' * (ϕw .* pϕ.P))[1:3,1:3] ≈ [2 -0.5 0; -0.5 2 -0.5; 0 -0.5 2]
        @test (pϕ' * (ϕw .* pϕ))[1:5,1:5] ≈ I

        ϕw = JacobiWeight(α+1,β+1) .* (b .- x)
        pϕ = LanczosPolynomial(ϕw)
        @test (pϕ.P' * (ϕw .* pϕ.P))[1:3,1:3] ≈ [2 -0.5 0; -0.5 2 -0.5; 0 -0.5 2]
        @test (pϕ' * (ϕw .* pϕ))[1:5,1:5] ≈ I
    end
end


text:
using AutoCorrelationResampling
using XUnit
using FFTW

@testset "Test resampling" begin include("resampling.jl") end


text:
using LightGraphs
using LightGraphs.SimpleGraphs
using LinearAlgebra
using DifferentialEquations
using ProgressMeter

include("utils\\randutils.jl")
include("utils\\categorical_tree.jl")
include("utils\\utils.jl")

struct Metaplex{T1<:AbstractSimpleGraph,T2<:AbstractSimpleGraph}
    g::T1 # global relation network between individuals
    h::T2 # relation network of components
end

abstract type EpidemicEvent end
struct InfectionEvent <: EpidemicEvent
    infectee::Int
    infector::Int
end

struct MigrationEvent <: EpidemicEvent
    individual::Int
    from::Int
    to::Int
end

function metaplex_gillespie(mp::Metaplex{T1,T2}, Xi0::BitVector, Xμ0::Vector{Int}, β, D::Vector; nmax=length(Xi0), tmax=100.0, sampling_method=:tree) where {T1,T2}
    N = nv(mp.g) # number of individuals
    g = mp.g
    M = nv(mp.h) # number of populations
    h = mp.h
    #od = map(k-> k>0 ? 1/float(k) : 0, outdegree(h))
    #od = outdegree(h)
    od = outdegree(h) .> 0
    Xi = copy(Xi0) # epidemic state of each individual
    Xμ = copy(Xμ0) # location of each individual
    popcounts = zeros(Int,2,M) # number of of individuals susceptible and infected for each component
    for k in 1:N
        if !Xi[k]
            #push!(pops_s[Xμ[k]],k)
            popcounts[1,Xμ[k]] += 1
        else
            #push!(pops_i[Xμ[k]],k)
            popcounts[2,Xμ[k]] += 1
        end
    end
    a = zeros(N) # aggregate reactions for each individual becoming infected and or migrating
    for k in 1:N
        if !Xi[k]
            l = length(filter(i->Xi[i] && Xμ[i]==Xμ[k], inneighbors(g,k)))
            a[k] = β*l # probability of k recieving infection is β time its number of infected incoming neighbors in the same population
            a[k] += D[1]
        else
            a[k] = D[2]
        end
    end

    if sampling_method == :tree
        a = CategoricalTree(a)
    end
    a0 = sum(a)
    ts = zeros(nmax)
    pcs = zeros(Int,nmax,2,M)
    pcs[1,:,:] .= popcounts
    t = 0.0
    n = 1
    while n < nmax && t < tmax && a0 > 0
        τ = log(1/rand())/a0
        k = rand_categorical(a, a0) # pick a random node
        μ = Xμ[k]
        if !Xi[k]
            if rand()*a[k] < D[1] # node k moves to another cluster
                ν = rand(outneighbors(h,μ))
                Xμ[k] = ν
                popcounts[1,μ] -= 1 # length(pops_s[μ])
                popcounts[1,ν] += 1 #length(pops_s[ν])
                l = length(filter(i->Xi[i] && Xμ[i]==ν, inneighbors(g,k)))
                a[k] = β*l + D[1]
            else # k gets infected
                Xi[k] = true
                popcounts[1,μ] -= 1 #length(pops_s[μ])
                popcounts[2,μ] += 1 #length(pops_i[μ])
                a[k] = D[2]
                for i in Iterators.filter(j->!Xi[j] && Xμ[j]==μ, outneighbors(g,k))
                    a[i] += β
                    #a[i] = D[1] + β*length(filter(i->Xi[i] && Xμ[i]==μ, inneighbors(g,i)))
                end
            end
        else
            ν = rand(outneighbors(h,μ))
            Xμ[k] = ν
            popcounts[2,μ] -= 1
            popcounts[2,ν] += 1
            for i in Iterators.filter(l->!Xi[l], outneighbors(g,k))
                if Xμ[i] == μ
                    a[i] -= β
                    #a[i] = D[1] + β*length(filter(i->Xi[i] && Xμ[i]==μ, inneighbors(g,i)))
                elseif Xμ[i] == ν
                    a[i] += β
                    #a[i] = D[1] + β*length(filter(i->Xi[i] && Xμ[i]==ν, inneighbors(g,i)))
                end
            end
        end
        a0 = sum(a)
        t += τ
        ts[n+1] = t
        pcs[n+1,:,:] .= popcounts
        n += 1
    end
    return ts[1:n], pcs[1:n,:,:]
end

function metaplex_ode(mp::Metaplex, Xi0::BitVector, Xμ0::Vector{Int}, β, D::Vector; tmax=100.0, saveat=[])
    N = nv(mp.g)
    M = nv(mp.h)
    A = adjacency_matrix(mp.g)
    L = normalized_laplacian(mp.h)
    u0 = zeros(2,N,M)
    for i in 1:N
        if !Xi0[i]
            u0[1,i,Xμ0[i]] = 1.0
        else
            u0[2,i,Xμ0[i]] = 1.0
        end
    end
    f! = function(du,u,p,t)
        for i in 1:N
            du[1,i,:] .= .-D[1]*L'*u[1,i,:]
            du[2,i,:] .= .-D[2]*L'*u[2,i,:]
        end
        for μ in 1:M
            infection = β*(u[1,:,μ].*(A*u[2,:,μ]) )# .+ u[1,:,μ].*u[2,:,μ])
            du[1,:,μ] .-= infection
            du[2,:,μ] .+= infection
        end
    end
    prob = ODEProblem(f!, u0, (0.0,tmax), N)
    sol = solve(prob, Tsit5(), saveat=saveat, abstol=1e-09)
    s = zeros(length(sol.t), N, M)
    i = zeros(length(sol.t), N, M)
    for t in 1:length(sol.t)
        s[t,:,:] = sol.u[t][1,:,:]
        i[t,:,:] = sol.u[t][2,:,:]
    end
    sμ = reshape(sum(s, dims=2), length(sol.t), M)
    iμ = reshape(sum(i, dims=2), length(sol.t), M)
    return sol.t, s, i, sμ, iμ, sol
end

function metaplex_montecarlo(mp::Metaplex, Xi0::BitVector, Xμ0::Vector{Int}, β, D::Vector; nmax=length(Xi0), tmax=100.0, nsims=1000, nbins=400)
    ts = LinRange(0.0, tmax, nbins)
    M = nv(mp.h)
    Xs_sum = zeros(nbins, M)
    Ms_sum = zeros(nbins, M)
    Xi_sum = zeros(nbins, M)
    Mi_sum = zeros(nbins, M) # running SSE
    Xn = Vector(undef, M)
    Δn = Vector(undef, M)
    p = Progress(nsims, dt=1.0)
    for n in 1:nsims
        t, pcs = metaplex_gillespie(mp, Xi0, Xμ0, β, D, nmax=nmax, tmax=tmax)
        l = 1
        for k in 1:nbins
            while l < length(t) && t[l+1] < ts[k]
                l += 1
            end
            Xn .= pcs[l,1,:]
            Δn .= Xn .- Xs_sum[k]/(n>1 ? n-1 : 1)
            Xs_sum[k,:] .+= Xn
            Ms_sum[k,:] .+= Δn.*(Xn.-Xs_sum[k]/n)
            Xn .= pcs[l,2,:]
            Δn .= Xn .- Xi_sum[k]/(n>1 ? n-1 : 1)
            Xi_sum[k,:] .+= Xn
            Mi_sum[k,:] .+= Δn.*(Xn.-Xi_sum[k]/n)
        end
        ProgressMeter.next!(p; showvalues= [(:n,n)])
    end
    return ts, Xs_sum/nsims, Xi_sum/nsims, sqrt.(abs.(Ms_sum/(nsims-1))), sqrt.(abs.(Mi_sum/(nsims-1)))
end

max_stars_repo_name: csimal/metaplex-networks

max_issues_repo_name: csimal/metaplex-networks

max_forks_repo_name: csimal/metaplex-networks


text:
using RealSense

include("rs-utils.jl")
include("gl-utils.jl")

# initialize GLFW window
glfwWidth = 1280
glfwHeight = 720
window = create_window_context(glfwWidth, glfwHeight)

# prepare drawing via OpenGL
window_frame = GLfloat[-1.0,-1.0,   1.0,-1.0,   1.0, 1.0,
                        1.0, 1.0,  -1.0, 1.0,  -1.0,-1.0]
vertex_texcoords = GLfloat[0.0,1.0,  1.0,1.0,  1.0,0.0,
                           1.0,0.0,  0.0,0.0,  0.0,1.0]
# create Vertex Buffer Objects
vbo = Ref{GLuint}(0)
glGenBuffers(1, vbo)
glBindBuffer(GL_ARRAY_BUFFER, vbo[])
glBufferData(GL_ARRAY_BUFFER, sizeof(window_frame), window_frame, GL_STATIC_DRAW)
texcoordsVBO = Ref{GLuint}(0)
glGenBuffers(1, texcoordsVBO)
glBindBuffer(GL_ARRAY_BUFFER, texcoordsVBO[])
glBufferData(GL_ARRAY_BUFFER, sizeof(vertex_texcoords), vertex_texcoords, GL_STATIC_DRAW)
# create Vertex Array Objects
vao = Ref{GLuint}(0)
glGenVertexArrays(1, vao)
glBindVertexArray(vao[])
glBindBuffer(GL_ARRAY_BUFFER, vao[])
glVertexAttribPointer(0, 2, GL_FLOAT, GL_FALSE, 0, C_NULL)
glBindBuffer(GL_ARRAY_BUFFER, texcoordsVBO[])
glVertexAttribPointer(1, 2, GL_FLOAT, GL_FALSE, 0, C_NULL)
glEnableVertexAttribArray(0)
glEnableVertexAttribArray(1)
# create shader program
vertSource= joinpath(@__DIR__, "rp.vert")
fragSource = joinpath(@__DIR__, "rp.frag")
shaderProgramID = createprogram(vertSource, fragSource)
# generate texture handle
tex = Ref{GLuint}(0)
glGenTextures(1, tex)

# initialize RealSense camera
err = Ref{Ptr{rs2_error}}(0)
ctx = rs2_create_context(RS2_API_VERSION, err)
checkerror(err)
device_list = rs2_query_devices(ctx, err)
checkerror(err)
dev_count = rs2_get_device_count(device_list, err)
checkerror(err)
@info "There are $dev_count connected RealSense devices."
_pipeline = rs2_create_pipeline(ctx, err)
checkerror(err)
# start with a playback config
config = rs2_create_config(err)
checkerror(err)
rs2_config_enable_device_from_file(config, joinpath(@__DIR__, "a.bag"), err)
checkerror(err)
# start the pipeline streaming
@static if is_apple()
    # we sleep 1s on MacOS, otherwise it will hit issue #1586:
    # https://github.com/IntelRealSense/librealsense/issues/1586
    sleep(1)
end
rs2_pipeline_start_with_config(_pipeline, config, err)
checkerror(err)
@assert err[] == C_NULL "Failed to start playback pipeline!"
pipeline_profile = rs2_pipeline_get_active_profile(_pipeline, err)
checkerror(err)
dev = rs2_pipeline_profile_get_device(pipeline_profile, err)
checkerror(err)

# create Depth-Colorizer processing block that can be used to quickly visualize the depth data
color_map = rs2_create_colorizer(err)
checkerror(err)
# create target output queue for the results
frame_queue = rs2_create_frame_queue(1, err)
checkerror(err)
# set the output of the processing block to be the queue
rs2_start_processing_queue(color_map, frame_queue, err)
checkerror(err)
@info "Playing back..."
# OpenGL rendering loop
while !GLFW.WindowShouldClose(window)
    updatefps!(window)
    # clear drawing surface
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glViewport(0, 0, glfwWidth, glfwHeight)
    glUseProgram(shaderProgramID)

    # poll next frame set and extract depth frame
    framesRef = Ref{Ptr{rs2_frame}}(C_NULL)
    has_frames = rs2_pipeline_poll_for_frames(_pipeline, framesRef, err)
    checkerror(err)
    frames = framesRef[]
    has_frames == 0 && (rs2_release_frame(frames); continue)
    depth_frame = get_depth_frame(frames)
    # select stream
    if depth_frame != C_NULL
        # frame reference is being "moved" into processing, so unless you
        # don't need frame anymore you must add an extra rs2_frame_add_ref
        rs2_frame_add_ref(depth_frame, err)
        checkerror(err)
        # invoke processing
        rs2_process_frame(color_map, depth_frame, err)
        checkerror(err)
        colorized_frame = rs2_wait_for_frame(frame_queue, 5000, err)
        checkerror(err)
        # get data
        depth_data = Ptr{UInt8}(rs2_get_frame_data(colorized_frame, err))
        checkerror(err)
        width = rs2_get_frame_width(depth_frame, err)
        height = rs2_get_frame_height(depth_frame, err)
        glBindTexture(GL_TEXTURE_2D, tex[])
        glTexImage2D(GL_TEXTURE_2D, 0, GL_RGB, width, height, 0, GL_RGB, GL_UNSIGNED_BYTE, depth_data)
        mipmap()
        glBindVertexArray(vao[])
        glDrawArrays(GL_TRIANGLES, 0, 6)
        # release resource
        rs2_release_frame(colorized_frame)
    end
    # release resource
    rs2_release_frame(depth_frame)
    rs2_release_frame(frames)
    # check and call events
    GLFW.PollEvents()
    # swap the buffers
    GLFW.SwapBuffers(window)
end

# release rs2 resources
rs2_delete_frame_queue(frame_queue)
rs2_delete_processing_block(color_map)
rs2_delete_pipeline_profile(pipeline_profile)
rs2_delete_config(config)
rs2_delete_pipeline(_pipeline)
rs2_delete_device(dev)
rs2_delete_device_list(device_list)
rs2_delete_context(ctx)

# release OpenGL resources
GLFW.DestroyWindow(window)


text:
import Escher: Elem, Tile, render
import Base: convert
import JSON: json

type VisPlot <: Tile
    nodes
    edges
end

nodes_json(v::VisJSGraph) = json([Dict(:id=>i, :label=>"Node $i") for i in 1:v.n])

function edges_json(v::VisJSGraph)
    edges = Vector{Dict}()
    for item in v.e
        push!(edges, Dict(:from=>item[1],:to=>item[2]))
    end
    json(edges)
end

convert(::Type{Tile}, v::VisJSGraph) = VisPlot(nodes_json(v), edges_json(v))

render(plot::VisPlot, state) =
    Elem(:"vis-plot", attributes = Dict(:nodes=>plot.nodes, :edges=>plot.edges))


text:
include("odehelper.jl")

function lorenz!(du,u,p,t)
    du[1] = 10.0*(u[2]-u[1])
    du[2] = u[1]*(28.0-u[3]) - u[2]
    du[3] = u[1]*u[2] - (8/3)*u[3]
end
solve_plot_print(lorenz!,[1.0;0.0;0.0],(0.0,100.0))
#http://www.math.utah.edu/~gustafso/2250systems-de.pdf
function brinetankvascade!(du,u,p,t)
    du[1] = -(1/2)*u[1]
    du[2] = (1/2) * u[1] - (1/4)*u[2]
    du[3] = (1/4)*u[2] - (1/6)*u[3]
end

solve_plot_print(lorenz!,[20,40,60],(0.0,10.0))

max_stars_repo_name: bkgsur/DifferentialEquationsWithJulia

max_issues_repo_name: bkgsur/DifferentialEquationsWithJulia

max_forks_repo_name: bkgsur/DifferentialEquationsWithJulia