text:
# an implementation of Uniform Manifold Approximation and Projection
# for Dimension Reduction, L. McInnes, J. Healy, J. Melville, 2018.

struct UMAP_{S <: Real, M <: AbstractMatrix{S}, N <: AbstractMatrix{S}, D<:AbstractMatrix{S}, K<:AbstractMatrix{<:Integer}, I<:AbstractMatrix{S}}
    graph::M
    embedding::N
    data::D
    knns::K
    dists::I

    function UMAP_{S, M, N, D, K, I}(graph, embedding, data, knns, dists) where {S<:Real,
                                                                                M<:AbstractMatrix{S},
                                                                                N<:AbstractMatrix{S},
                                                                                D<:AbstractMatrix{S},
                                                                                K<:AbstractMatrix{<:Integer},
                                                                                I<:AbstractMatrix{S}}
        issymmetric(graph) || isapprox(graph, graph') || error("UMAP_ constructor expected graph to be a symmetric matrix")
        size(knns) == size(dists) || error("UMAP_ constructor expected knns and dists to have equal size")
        new(graph, embedding, data, knns, dists)
    end
end

function UMAP_(graph::M, embedding::N, data::D, knns::K, dists::I) where {S<:Real,
                                                                          M<:AbstractMatrix{S},
                                                                          N<:AbstractMatrix{S},
                                                                          D<:AbstractMatrix{S},
                                                                          K<:AbstractMatrix{<:Integer},
                                                                          I<:AbstractMatrix{S}}
    return UMAP_{S, M, N, D, K, I}(graph, embedding, data, knns, dists)
end

const SMOOTH_K_TOLERANCE = 1e-5

"""
    umap(X::AbstractMatrix[, n_components=2]; <kwargs>) -> embedding

Embed the data `X` into a `n_components`-dimensional space. `n_neighbors` controls
how many neighbors to consider as locally connected.

See `UMAP_` for a description of keyword arguments.
"""
function umap(args...; kwargs...)
    # this is just a convenience function for now
    return UMAP_(args...; kwargs...).embedding
end

"""
    UMAP_(X::AbstractMatrix[, n_components=2]; <kwargs>) -> UMAP_ object

Create a model representing the embedding of data `X` into `n_components`-dimensional space.
The returned model has the following fields:

- `graph`: the graph representing the fuzzy simplicial set of the manifold of `X`.
- `embedding`: the `n-component`-dimensional embedding of the data `X`.
- `data`: a reference to the input data `X`.
- `knns`: a matrix of indices of `X` representing each point's nearest neighbors according to `metric`.
          `knns[j, i]` is the index of point i's jth nearest neighbor.
- `dists`: the respective distances of the above neighbors.
           `dists[j, i]` is the distance of point i's jth nearest neighbor.

# Keyword Arguments
- `n_neighbors::Integer = 15`: the number of neighbors to consider as locally connected. Larger values capture more global structure in the data, while small values capture more local structure.
- `metric::{SemiMetric, Symbol} = Euclidean()`: the metric to calculate distance in the input space. It is also possible to pass `metric = :precomputed` to treat `X` like a precomputed distance matrix.
- `n_epochs::Integer = 300`: the number of training epochs for embedding optimization
- `learning_rate::Real = 1`: the initial learning rate during optimization
- `init::Symbol = :spectral`: how to initialize the output embedding; valid options are `:spectral` and `:random`
- `min_dist::Real = 0.1`: the minimum spacing of points in the output embedding
- `spread::Real = 1`: the effective scale of embedded points. Determines how clustered embedded points are in combination with `min_dist`.
- `set_operation_ratio::Real = 1`: interpolates between fuzzy set union and fuzzy set intersection when constructing the UMAP graph (global fuzzy simplicial set). The value of this parameter should be between 1.0 and 0.0: 1.0 indicates pure fuzzy union, while 0.0 indicates pure fuzzy intersection.
- `local_connectivity::Integer = 1`: the number of nearest neighbors that should be assumed to be locally connected. The higher this value, the more connected the manifold becomes. This should not be set higher than the intrinsic dimension of the manifold.
- `repulsion_strength::Real = 1`: the weighting of negative samples during the optimization process.
- `neg_sample_rate::Integer = 5`: the number of negative samples to select for each positive sample. Higher values will increase computational cost but result in slightly more accuracy.
- `a = nothing`: this controls the embedding. By default, this is determined automatically by `min_dist` and `spread`.
- `b = nothing`: this controls the embedding. By default, this is determined automatically by `min_dist` and `spread`.
"""
function UMAP_(X::AbstractMatrix{S},
               n_components::Integer = 2;
               n_neighbors::Integer = 15,
               metric::Union{SemiMetric, Symbol} = Euclidean(),
               n_epochs::Integer = 300,
               learning_rate::Real = 1,
               init::Symbol = :spectral,
               min_dist::Real = 1//10,
               spread::Real = 1,
               set_operation_ratio::Real = 1,
               local_connectivity::Integer = 1,
               repulsion_strength::Real = 1,
               neg_sample_rate::Integer = 5,
               a::Union{Real, Nothing} = nothing,
               b::Union{Real, Nothing} = nothing
               ) where {S<:Real}
    # argument checking
    size(X, 2) > n_neighbors > 0|| throw(ArgumentError("size(X, 2) must be greater than n_neighbors and n_neighbors must be greater than 0"))
    size(X, 1) > n_components > 1 || throw(ArgumentError("size(X, 1) must be greater than n_components and n_components must be greater than 1"))
    n_epochs > 0 || throw(ArgumentError("n_epochs must be greater than 0"))
    learning_rate > 0 || throw(ArgumentError("learning_rate must be greater than 0"))
    min_dist > 0 || throw(ArgumentError("min_dist must be greater than 0"))
    0 ≤ set_operation_ratio ≤ 1 || throw(ArgumentError("set_operation_ratio must lie in [0, 1]"))
    local_connectivity > 0 || throw(ArgumentError("local_connectivity must be greater than 0"))

    # main algorithm
    knns, dists = knn_search(X, n_neighbors, metric)
    graph = fuzzy_simplicial_set(knns, dists, n_neighbors, size(X, 2), local_connectivity, set_operation_ratio)

    embedding = initialize_embedding(graph, n_components, Val(init))

    embedding = optimize_embedding(graph, embedding, embedding, n_epochs, learning_rate, min_dist, spread, repulsion_strength, neg_sample_rate, move_ref=true)
    # TODO: if target variable y is passed, then construct target graph
    #       in the same manner and do a fuzzy simpl set intersection

    return UMAP_(graph, hcat(embedding...), X, knns, dists)
end

"""
    transform(model::UMAP_, Q::AbstractMatrix; <kwargs>) -> embedding

Use the given model to embed new points into an existing embedding. `Q` is a matrix of some number of points (columns)
in the same space as `model.data`. The returned embedding is the embedding of these points in n-dimensional space, where
n is the dimensionality of `model.embedding`. This embedding is created by finding neighbors of `Q` in `model.embedding`
and optimizing cross entropy according to membership strengths according to these neighbors.

# Keyword Arguments
- `n_neighbors::Integer = 15`: the number of neighbors to consider as locally connected. Larger values capture more global structure in the data, while small values capture more local structure.
- `metric::{SemiMetric, Symbol} = Euclidean()`: the metric to calculate distance in the input space. It is also possible to pass `metric = :precomputed` to treat `X` like a precomputed distance matrix.
- `n_epochs::Integer = 300`: the number of training epochs for embedding optimization
- `learning_rate::Real = 1`: the initial learning rate during optimization
- `init::Symbol = :spectral`: how to initialize the output embedding; valid options are `:spectral` and `:random`
- `min_dist::Real = 0.1`: the minimum spacing of points in the output embedding
- `spread::Real = 1`: the effective scale of embedded points. Determines how clustered embedded points are in combination with `min_dist`.
- `set_operation_ratio::Real = 1`: interpolates between fuzzy set union and fuzzy set intersection when constructing the UMAP graph (global fuzzy simplicial set). The value of this parameter should be between 1.0 and 0.0: 1.0 indicates pure fuzzy union, while 0.0 indicates pure fuzzy intersection.
- `local_connectivity::Integer = 1`: the number of nearest neighbors that should be assumed to be locally connected. The higher this value, the more connected the manifold becomes. This should not be set higher than the intrinsic dimension of the manifold.
- `repulsion_strength::Real = 1`: the weighting of negative samples during the optimization process.
- `neg_sample_rate::Integer = 5`: the number of negative samples to select for each positive sample. Higher values will increase computational cost but result in slightly more accuracy.
- `a = nothing`: this controls the embedding. By default, this is determined automatically by `min_dist` and `spread`.
- `b = nothing`: this controls the embedding. By default, this is determined automatically by `min_dist` and `spread`.
"""
function transform(model::UMAP_,
                   Q::AbstractMatrix{S};
                   n_neighbors::Integer = 15,
                   metric::Union{SemiMetric, Symbol} = Euclidean(),
                   n_epochs::Integer = 100,
                   learning_rate::Real = 1,
                   min_dist::Real = 1//10,
                   spread::Real = 1,
                   set_operation_ratio::Real = 1,
                   local_connectivity::Integer = 1,
                   repulsion_strength::Real = 1,
                   neg_sample_rate::Integer = 5,
                   a::Union{Real, Nothing} = nothing,
                   b::Union{Real, Nothing} = nothing
                   ) where {S<:Real}
    # argument checking
    size(Q, 2) > n_neighbors > 0                     || throw(ArgumentError("size(Q, 2) must be greater than n_neighbors and n_neighbors must be greater than 0"))
    learning_rate > 0                                || throw(ArgumentError("learning_rate must be greater than 0"))
    min_dist > 0                                     || throw(ArgumentError("min_dist must be greater than 0"))
    0 ≤ set_operation_ratio ≤ 1                      || throw(ArgumentError("set_operation_ratio must lie in [0, 1]"))
    local_connectivity > 0                           || throw(ArgumentError("local_connectivity must be greater than 0"))
    size(model.data, 2) == size(model.embedding, 2)  || throw(ArgumentError("model.data must have same number of columns as model.embedding"))
    size(model.data, 1) == size(Q, 1)                || throw(ArgumentError("size(model.data, 1) must equal size(Q, 1)"))

    n_epochs = max(0, n_epochs)
    # main algorithm
    knns, dists = knn_search(model.data, Q, n_neighbors, metric, model.knns, model.dists)
    graph = fuzzy_simplicial_set(knns, dists, n_neighbors, size(model.data, 2), local_connectivity, set_operation_ratio, false)

    embedding = initialize_embedding(graph, model.embedding)
    ref_embedding = collect(eachcol(model.embedding))
    embedding = optimize_embedding(graph, embedding, ref_embedding, n_epochs, learning_rate, min_dist, spread, repulsion_strength, neg_sample_rate, a, b, move_ref=false)

    return reduce(hcat, embedding)
end

"""
    fuzzy_simplicial_set(knns, dists, n_neighbors, n_points, local_connectivity, set_op_ratio, apply_fuzzy_combine=true) -> membership_graph::SparseMatrixCSC,

Construct the local fuzzy simplicial sets of each point represented by its distances
to its `n_neighbors` nearest neighbors, stored in `knns` and `dists`, normalizing the distances
on the manifolds, and converting the metric space to a simplicial set.
`n_points` indicates the total number of points of the original data, while `knns` contains
indices of some subset of those points (ie some subset of 1:`n_points`). If `knns` represents
neighbors of the elements of some set with itself, then `knns` should have `n_points` number of
columns. Otherwise, these two values may be inequivalent.
If `apply_fuzzy_combine` is true, use intersections and unions to combine
fuzzy sets of neighbors (default true).

The returned graph will have size (`n_points`, size(knns, 2)).
"""
function fuzzy_simplicial_set(knns,
                              dists,
                              n_neighbors,
                              n_points,
                              local_connectivity,
                              set_operation_ratio,
                              apply_fuzzy_combine=true)

    σs, ρs = smooth_knn_dists(dists, n_neighbors, local_connectivity)

    rows, cols, vals = compute_membership_strengths(knns, dists, σs, ρs)
    fs_set = sparse(rows, cols, vals, n_points, size(knns, 2))

    if apply_fuzzy_combine
        res = combine_fuzzy_sets(fs_set, set_operation_ratio)
        return dropzeros(res)
    else
        return dropzeros(fs_set)
    end
end

"""
    smooth_knn_dists(dists, k, local_connectivity; <kwargs>) -> knn_dists, nn_dists

Compute the distances to the nearest neighbors for a continuous value `k`. Returns
the approximated distances to the kth nearest neighbor (`knn_dists`)
and the nearest neighbor (nn_dists) from each point.
"""
function smooth_knn_dists(knn_dists::AbstractMatrix{S},
                          k::Real,
                          local_connectivity::Real;
                          niter::Integer=64,
                          bandwidth::Real=1) where {S <: Real}

    nonzero_dists(dists) = @view dists[dists .> 0.]
    ρs = zeros(S, size(knn_dists, 2))
    σs = Array{S}(undef, size(knn_dists, 2))
    for i in 1:size(knn_dists, 2)
        nz_dists = nonzero_dists(knn_dists[:, i])
        if length(nz_dists) >= local_connectivity
            index = floor(Int, local_connectivity)
            interpolation = local_connectivity - index
            if index > 0
                ρs[i] = nz_dists[index]
                if interpolation > SMOOTH_K_TOLERANCE
                    ρs[i] += interpolation * (nz_dists[index+1] - nz_dists[index])
                end
            else
                ρs[i] = interpolation * nz_dists[1]
            end
        elseif length(nz_dists) > 0
            ρs[i] = maximum(nz_dists)
        end
        @inbounds σs[i] = smooth_knn_dist(knn_dists[:, i], ρs[i], k, bandwidth, niter)
    end

    return ρs, σs
end

# calculate sigma for an individual point
@fastmath function smooth_knn_dist(dists::AbstractVector, ρ, k, bandwidth, niter)
    target = log2(k)*bandwidth
    lo, mid, hi = 0., 1., Inf
    for n in 1:niter
        psum = sum(exp.(-max.(dists .- ρ, 0.)./mid))
        if abs(psum - target) < SMOOTH_K_TOLERANCE
            break
        end
        if psum > target
            hi = mid
            mid = (lo + hi)/2.
        else
            lo = mid
            if hi == Inf
                mid *= 2.
            else
                mid = (lo + hi) / 2.
            end
        end
    end
    # TODO: set according to min k dist scale
    return mid
end

"""
    compute_membership_strengths(knns, dists, σs, ρs) -> rows, cols, vals

Compute the membership strengths for the 1-skeleton of each fuzzy simplicial set.
"""
function compute_membership_strengths(knns::AbstractMatrix{S},
                                      dists::AbstractMatrix{T},
                                      ρs::Vector{T},
                                      σs::Vector{T}) where {S <: Integer, T}
    # set dists[i, j]
    rows = sizehint!(S[], length(knns))
    cols = sizehint!(S[], length(knns))
    vals = sizehint!(T[], length(knns))
    for i in 1:size(knns, 2), j in 1:size(knns, 1)
        @inbounds if i == knns[j, i] # dist to self
            d = 0.
        else
            @inbounds d = exp(-max(dists[j, i] - ρs[i], 0.)/σs[i])
        end
        append!(cols, i)
        append!(rows, knns[j, i])
        append!(vals, d)
    end
    return rows, cols, vals
end


text:
# Scaled dot-product attention function
function scaled_dot_product_attention(q, k, v)
    x = nGraph.bmm(q, k; transpose_b = true)
    x = x ./ size(k, 1)
    x = Flux.softmax(x)
    return nGraph.bmm(x, v)
end

# multi-head attention
#
# Arguments `q`, `k`, and `v` are the queries, keys, and values
# Arguments `wq`, `wk`, and `wv` vectors of projection matrices
function multi_head_attention(q, k, v, wq, wk, wv, wo)
    heads = map(zip(wq, wk, wv)) do x
        # Do a broadcast operation
        return scaled_dot_product_attention(x[1] * q, x[2] * k, x[3] * v)
    end

    return wo * cat(heads...; dims = 1)
end

# Point wise feed-forward networks
feed_forward(x, W::Vector, B::Vector) = W[2] * Flux.relu(W[1] * x .+ B[1]) .+ B[2]

function layer_normalization(a, g, b)
    μ = sum(a) ./ size(a, 1)
    diff = a .- μ
    σ = sqrt(sum((diff .* diff) ./ size(a, 1)))
    return ((g ./ σ) .* diff) .+ b
end

#####
##### Helpers for building the whole network
#####
struct AttentionLayer
    wq::Vector
    wk::Vector
    wv::Vector
    wo
end

function AttentionLayer(dmodel, h)
    dk = div(dmodel, h)
    wq = [param(randn(Float32, dk, dmodel)) for _ in 1:h]
    wk = [param(randn(Float32, dk, dmodel)) for _ in 1:h]
    wv = [param(randn(Float32, dk, dmodel)) for _ in 1:h]
    wo = param(randn(Float32, dmodel, dmodel))

    return AttentionLayer(wq, wk, wv, wo)
end

(A::AttentionLayer)(q, k, v) = multi_head_attention(q, k, v, A.wq, A.wk, A.wv, A.wo)

struct FeedForward
    W::Vector
    B::Vector
end

function FeedForward(dmodel::Integer, dff::Integer)
    W1 = param(randn(Float32, dff, dmodel))
    B1 = param(randn(Float32, dff))

    W2 = param(randn(Float32, dmodel, dff))
    B2 = param(randn(Float32, dmodel))
    return FeedForward([W1, W2], [B1, B2])
end

(F::FeedForward)(x) = feed_forward(x, F.W, F.B)

struct Normalize
    g
    b
end

function Normalize(dmodel::Integer)
    g = param(randn(Float32, dmodel))
    b = param(randn(Float32, dmodel))
    return Normalize(g, b)
end

(N::Normalize)(x) = layer_normalization(x, N.g, N.b)

# Encoder Layers
struct Encoder
    attentions::Vector{AttentionLayer}
    att_normalize::Vector{Normalize}
    forwards::Vector{FeedForward}
    fwd_normalize::Vector{Normalize}
end

function Encoder(N::Integer, dmodel::Integer, dff::Integer, h::Integer)
    attentions = [AttentionLayer(dmodel, h) for _ in 1:N]
    att_normalize = [Normalize(dmodel) for _ in 1:N]
    forwards = [FeedForward(dmodel, dff) for _ in 1:N]
    fwd_normalize = [Normalize(dmodel) for _ in 1:N]
    return Encoder(attentions, att_normalize, forwards, fwd_normalize)
end

function (E::Encoder)(x)
    iter = zip(E.attentions, E.att_normalize, E.forwards, E.fwd_normalize)
    for (attention, att_n, forward, fwd_n) in iter
        x = att_n(x + attention(x, x, x))
        x = fwd_n(x + forward(x))
    end
    return x
end

# Decoder Layers
struct Decoder
    first_attentions::Vector{AttentionLayer}
    first_normalize::Vector{Normalize}
    second_attentions::Vector{AttentionLayer}
    second_normalize::Vector{Normalize}
    forwards::Vector{FeedForward}
    fwd_normalize::Vector{Normalize}
end

function Decoder(N::Integer, dmodel::Integer, dff::Integer, h::Integer)
    first_attentions = [AttentionLayer(dmodel, h) for _ in 1:N]
    first_normalize = [Normalize(dmodel) for _ in 1:N]
    second_attentions = [AttentionLayer(dmodel, h) for _ in 1:N]
    second_normalize = [Normalize(dmodel) for _ in 1:N]
    forwards = [FeedForward(dmodel, dff) for _ in 1:N]
    fwd_normalize = [Normalize(dmodel) for _ in 1:N]

    return Decoder(
        first_attentions,
        first_normalize,
        second_attentions,
        second_normalize,
        forwards,
        fwd_normalize
    )
end

function (D::Decoder)(x, hidden)
    iter = zip(
        D.first_attentions,
        D.first_normalize,
        D.second_attentions,
        D.second_normalize,
        D.forwards,
        D.fwd_normalize
    )
    for (first_attention, first_n, second_attention, second_n, forward, fwd_n) in iter
        x = first_n(x + first_attention(x, x, x))
        x = second_n(x + second_attention(hidden, hidden, x))
        x = fwd_n(x + forward(x))
    end
    return x
end

# Top level entry point
struct Transformer
    encoder::Encoder
    decoder::Decoder
    embedding
    vocab_size::Int
end

function Transformer(N::Integer, dmodel::Integer, dff::Integer, h::Integer, vocab, seqlen)
    hidden_size = dmodel * seqlen
    embedding = Flux.param(randn(Float32, dmodel, vocab))

    encoder = Encoder(N, dmodel, dff, h)
    decoder = Decoder(N, dmodel, dff, h)
    return Transformer(encoder, decoder, embedding, vocab)
end

function (T::Transformer)(inputs, outputs, target)
    embedding_node = nGraph.Node(T.embedding)

    # Tell the Julia-to-nGraph converter to setup the embedding table as an inplace op.
    nGraph.__inplace(T.embedding)

    # Do the embedding table lookup
    lookup_inputs = nGraph.embedding(inputs, embedding_node)
    lookup_outputs = nGraph.embedding(outputs, embedding_node)

    hidden = T.encoder(lookup_inputs)
    y = T.decoder(lookup_outputs, hidden)
    # Reshape the first two dimensions into one dimension
    #y = transpose(reshape(y, :, size(y, 3)))
    y = transpose(y)
    y = y * embedding_node
    y = transpose(y)

    # Generate a loss
    return Flux.crossentropy(reshape(y, :), reshape(nGraph.onehot(target, T.vocab_size, 3), :))
end

function transformer_training(batchsize = 16, seq_length = 50)
    # Parameters from the Transformer Paper

    # Number of Stacked Layers
    N = 8
    # Size of the first tensor dimension
    dmodel = 512
    # Hidden size for the feed forward networks
    dff = 2048
    # Number of shards for the attention layers
    h = 8
    # Output vocab size
    vocab = 20000

    # Create inputs, outputs, and expected output layers
    selection = Int32(1):Int32(vocab)
    #X = nGraph.Node(randn(Int32, dmodel, seq_length, batchsize))
    #Y = nGraph.Node(randn(Float32, dmodel, seq_length, batchsize))
    X = nGraph.Node(rand(selection, seq_length, batchsize))
    Y = nGraph.Node(rand(selection, seq_length, batchsize))
    #expected = randn(Float32, vocab, seq_length, batchsize)
    expected = nGraph.Node(rand(selection, seq_length, batchsize))

    T = Transformer(N, dmodel, dff, h, vocab, seq_length)

    kw = (optimizer = nGraph.SGD(Float32(0.001)),)
    return T, (X, Y, expected), kw
end

ext: jl

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
struct RiskParameters{T} <: AbstractRisk{T}
  sparks::Union{Nothing, AbstractVector}
  susceptibility::Union{Nothing, AbstractVector}
  infectivity::Union{Nothing, AbstractVector}
  transmissibility::Union{Nothing, AbstractVector}
  latency::Union{Nothing, AbstractVector}
  removal::Union{Nothing, AbstractVector}

  RiskParameters{SEIR}(ϵ, θs, κ, θt, θl, θr) = new{SEIR}(ϵ, θs, κ, θt, θl, θr)
  RiskParameters{SEI}(ϵ, θs, κ, θt, θl) = new{SEI}(ϵ, θs, κ, θt, θl, nothing)
  RiskParameters{SIR}(ϵ, θs, κ, θt, θr) = new{SIR}(ϵ, θs, κ, θt, nothing, θr)
  RiskParameters{SI}(ϵ, θs, κ, θt) = new{SI}(ϵ, θs, κ, θt, nothing, nothing)
end

function Base.show(io::IO, x::RiskParameters{T}) where T <: DiseaseStateSequence
  return print(io, "$T model risk function parameters")
end

function Base.convert(::Type{Vector{Float64}},
                      x::RiskParameters{T}) where T <: DiseaseStateSequence
  return [x[i] for i = 1:length(x)]
end

function Base.convert(::Type{Array{Float64, 2}},
                      x::Vector{RiskParameters{T}}) where T <: DiseaseStateSequence
  return [x[i][j] for i = 1:length(x), j = 1:length(x[1])]
end

function _like(x::RiskParameters{T}, v::Vector{Float64}) where T <: DiseaseStateSequence
  indices = _indices(x, zeros=false)
  if indices[end] != length(v)
    @error "Incompatiable parameter vector"
  end
  if T == SEIR
    return RiskParameters{T}(v[1:(indices[1])],
                             v[(indices[1]+1):(indices[2])],
                             v[(indices[2]+1):(indices[3])],
                             v[(indices[3]+1):(indices[4])],
                             v[(indices[4]+1):(indices[5])],
                             v[(indices[5]+1):(indices[6])])
  elseif T in (SEI, SIR)
    return RiskParameters{T}(v[1:(indices[1])],
                             v[(indices[1]+1):(indices[2])],
                             v[(indices[2]+1):(indices[3])],
                             v[(indices[3]+1):(indices[4])],
                             v[(indices[4]+1):(indices[5])])
  elseif T == SI
    return RiskParameters{T}(v[1:(indices[1])],
                             v[(indices[1]+1):(indices[2])],
                             v[(indices[2]+1):(indices[3])],
                             v[(indices[3]+1):(indices[4])])
  end
end


text:
using PyCall
@pyimport pylagrit
@pyimport numpy as np

# User parameters #######################
# Discretization
nx = 61; ny = 61
# Polygon delineation filename
polygon_file = "2d_mesh.avs"
# User parameters #######################

l = pylagrit.PyLaGriT()

# Read in polygon delineation mesh
mopoly = l[:read](polygon_file)

# Create quad using polygon min and max coordinates
xs = linspace(mopoly[:xmin],mopoly[:xmax],nx)
ys = linspace(mopoly[:ymin],mopoly[:ymax],ny)
#xs = np.linspace(mopoly.xmin,mopoly.xmax,nx)
#ys = np.linspace(mopoly.ymin,mopoly.ymax,ny)
moquad = l[:gridder](xs,ys,elem_type="quad",connect=true)
mopoly[:setatt]("itetclr",1)

# Create extruded quad to tet and center at z=0
moprism = mopoly[:extrude](1.0)
moprism[:trans]((0.,0.,0.),(0.,0.,0.5))
motet = moprism[:grid2grid_prismtotet3]()

# Clean up a little
mopoly[:delete]()
moprism[:delete]()

# Add attribute to moquad to store if in polygon or not
moquad[:addatt]("if_in_poly",vtype="vint")
moquad[:setatt]("if_in_poly", 0) # Make sure attr starts as zero

# Interpolate motet's itetclr into moquad's if_in_poly
moquad[:interpolate_map]("if_in_poly",motet,"itetclr")
moquad[:intersect_elements](motet,attr_name="if_inter")

# Write out file, mapping is stored as cell attribute in if_inter
moquad[:dump]("polygon_map_julia.inp")

max_forks_repo_name: daniellivingston/LaGriT


text:
using SingleDispatchArrays
using SingleDispatchArrays: index_value_pair, add_subtype
using Test

@testset "Basic functionality" begin
    struct Point{T}
        x::T
        y::T
    end

    points = SingleDispatchVector{Point}()
    point_types = [Point{Int}, Point{Float32}, Point{Float64}]
    add_subtype.(Point, point_types)

    @test_throws ErrorException index_value_pair(points, Point{Bool}(false, false))

    push!(points, Point{Int}(3,4), Point{Float64}(12,5), Point{Float32}(1,0))

    @test_throws ErrorException push!(points, Point{Bool}(false, false))

    sum1 = 0.
    foreach(points) do p
        sum1 += hypot(p.x, p.y)
    end

    sum2 = 0.
    for p in points
        sum2 += hypot(p.x, p.y)
    end

    @test sum1 == sum2

end


text:
#ENV["CUDA_VISIBLE_DEVICES"] = "0,1,2,3,4,5,6,7"
#using Revise
using ADSeismic
using ADCME
using PyPlot
using MAT
using DelimitedFiles
matplotlib.use("Agg")
close("all")
if has_gpu()
  gpu = true
else
  gpu = false
end

output_dir = "../../../examples/nn_fwi/data/acoustic"
if !ispath(output_dir)
  mkpath(output_dir)
end

################### Generate synthetic data #####################
model_name = "../../../examples/nn_fwi/models/marmousi2-model-true.mat"
# model_name = "models/BP-model-true.mat"

## load model setting
params = load_params(model_name, vp_ref=1000, PropagatorKernel=0)
src = load_acoustic_source(model_name)
rcv = load_acoustic_receiver(model_name)
vp = constant(matread(model_name)["vp"])
## For debug, only run the first source
# src = [src[1]]
# rcv = [rcv[1]]

## assemble acoustic propagator model
model = x->AcousticPropagatorSolver(params, x, vp)

## simulated wavefield
if gpu
  Rs_ = compute_forward_GPU(model, src, rcv)
else
  [SimulatedObservation!(model(src[i]), rcv[i]) for i = 1:length(src)]
  Rs_ = [rcv[i].rcvv for i = 1:length(rcv)]
end
sess = Session(); init(sess)
Rs = run(sess, Rs_)

## save results
for i = 1:length(src)
    writedlm(joinpath(output_dir, "marmousi-r$i.txt"), Rs[i])
    # writedlm(joinpath(output_dir, "BP-r$i.txt"), Rs[i])
end
@info "Data written."

## visualize_wavefield if needed
ii = 1 # source number
u = run(sess, model(src[ii]).u)
p = visualize_wavefield(u, params)
saveanim(p, "wavefield.gif")


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

mutable struct ManagementPolicySchema <: SwaggerModel
    rules::Any # spec type: Union{ Nothing, Vector{ManagementPolicyRule} } # spec name: rules

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

const _property_map_ManagementPolicySchema = Dict{Symbol,Symbol}(Symbol("rules")=>Symbol("rules"))
const _property_types_ManagementPolicySchema = Dict{Symbol,String}(Symbol("rules")=>"Vector{ManagementPolicyRule}")
Base.propertynames(::Type{ ManagementPolicySchema }) = collect(keys(_property_map_ManagementPolicySchema))
Swagger.property_type(::Type{ ManagementPolicySchema }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ManagementPolicySchema[name]))}
Swagger.field_name(::Type{ ManagementPolicySchema }, property_name::Symbol) =  _property_map_ManagementPolicySchema[property_name]

function check_required(o::ManagementPolicySchema)
    (getproperty(o, Symbol("rules")) === nothing) && (return false)
    true
end

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


text:
@testset "PIPotential"  begin

#---

using ACE
using Printf, Test, LinearAlgebra, JuLIP, JuLIP.Testing, Random
using JuLIP: evaluate, evaluate_d, evaluate_ed
using ACE: combine

#---

@info("Basic test of PIPotential construction and evaluation")
maxdeg = 10
r0 = 1.0
rcut = 3.0
trans = PolyTransform(1, r0)
Pr = transformed_jacobi(maxdeg, trans, rcut; pcut = 2)
D = ACE.SparsePSHDegree()
P1 = ACE.RnYlm1pBasis(Pr; species = :X, D = D)
basis = ACE.PIBasis(P1, 2, D, maxdeg)
c = ACE.Random.randcoeffs(basis)
Vdag = combine(basis, c)
V = standardevaluator(Vdag)
Nat = 15
Rs, Zs, z0 = ACE.rand_nhd(Nat, Pr, :X)
val_basis = real(sum(c .* evaluate(basis, Rs, Zs, z0)))
val_V = evaluate(V, Rs, Zs, z0)
println(@test(val_basis ≈ val_V))
println(@test(evaluate(Vdag, Rs, Zs, z0) ≈ val_V))
J = evaluate_d(basis, Rs, Zs, z0)
grad_basis = real(sum(c[i] * J[i,:] for i = 1:length(c)))[:]
grad_V = evaluate_d(V, Rs, Zs, z0)
println(@test(grad_basis ≈ grad_V))
println(@test(evaluate_d(Vdag, Rs, Zs, z0) ≈ grad_V))

println(@test(all(JuLIP.Testing.test_fio(V))))

#---

# check multi-species
maxdeg = 5
Pr = transformed_jacobi(maxdeg, trans, rcut; pcut = 2)
species = [:C, :O, :H]
P1 = ACE.RnYlm1pBasis(Pr; species = [:C, :O, :H], D = D)
basis = ACE.PIBasis(P1, 3, D, maxdeg)
c = randcoeffs(basis)
Vdag = combine(basis, c)
V = standardevaluator(Vdag)
Rs, Zs, z0 = ACE.rand_nhd(Nat, Pr, species)
AA = evaluate(basis, Rs, Zs, z0)
val_basis = real(sum(c .* evaluate(basis, Rs, Zs, z0)))
val_V = evaluate(V, Rs, Zs, z0)
println(@test(val_basis ≈ val_V))
println(@test(evaluate(Vdag, Rs, Zs, z0) ≈ val_V))
J = evaluate_d(basis, Rs, Zs, z0)
grad_basis = real(sum(c[i] * J[i,:] for i = 1:length(c)))[:]
grad_V = evaluate_d(V, Rs, Zs, z0)
println(@test(grad_basis ≈ grad_V))
println(@test(evaluate_d(Vdag, Rs, Zs, z0) ≈ grad_V))

println(@test(all(JuLIP.Testing.test_fio(V))))

#---

@info("Check several properties of PIPotential")
for species in (:X, :Si, [:C, :O, :H]), N = 1:5
   local Rs, Zs, z0, V, Vdag, basis, P1, maxdeg, Nat, c, val_basis, val_V
   maxdeg = 7
   Nat = 15
   P1 = ACE.RnYlm1pBasis(Pr; species = species)
   basis = ACE.PIBasis(P1, N, D, maxdeg)
   @info("species = $species; N = $N; length = $(length(basis))")
   c = randcoeffs(basis)
   Vdag = combine(basis, c)
   V = standardevaluator(Vdag)
   @info("check (de-)serialisation")
   println(@test(all(JuLIP.Testing.test_fio(Vdag))))
   @info("Check basis and potential match")
   for ntest = 1:20
      Rs, Zs, z0 = ACE.rand_nhd(Nat, Pr, species)
      val_basis = real(sum(c .* evaluate(basis, Rs, Zs, z0)))
      val_V = evaluate(V, Rs, Zs, z0)
      print_tf(@test(val_basis ≈ val_V))
   end
   println()
   @info("Check gradients")
   for ntest = 1:20
      Rs, Zs, z0 = ACE.rand_nhd(Nat, Pr, species)
      V0 = evaluate(V, Rs, Zs, z0)
      dV0 = evaluate_d(V, Rs, Zs, z0)
      Us = [ rand(eltype(Rs)) .- 0.5 for _=1:length(Rs) ]
      dV0_dUs = sum(transpose.(dV0) .* Us)
      errs = []
      for p = 2:12
         h = 0.1^p
         V_h = evaluate(V, Rs + h * Us, Zs, z0)
         dV_h = (V_h - V0) / h
         # @show norm(dAA_h - dAA_dUs, Inf)
         push!(errs, norm(dV_h - dV0_dUs, Inf))
      end
      success = (/(extrema(errs)...) < 1e-3) || (minimum(errs) < 1e-10)
      print_tf(@test success)
   end
   println()
   @info("Check graph evaluator")
   for ntest = 1:20
      Rs, Zs, z0 = ACE.rand_nhd(Nat, Pr, species)
      v = evaluate(V, Rs, Zs, z0)
      vgr = evaluate(Vdag, Rs, Zs, z0)
      print_tf(@test(v ≈ vgr))
   end
   println()
end
println()

#---

@info("Check Correctness of SHIP.PIPotential calculators")

naive_energy(V::PIPotential, at) =
      sum( evaluate(V, Rs, at.Z[j], at.Z[i])
            for (i, j, Rs) in sites(at, cutoff(V)) )

for N = 1:5
   species = :Si
   maxdeg = 7
   Pr = transformed_jacobi(maxdeg, trans, 4.0; pcut = 2)
   P1 = ACE.RnYlm1pBasis(Pr; species = species)
   basis = ACE.PIBasis(P1, N, D, maxdeg)
   @info("N = $N; length = $(length(basis))")
   c = randcoeffs(basis)
   V = combine(basis, c)
   at = bulk(:Si) * (2,2,3)
   rattle!(at, 0.1)
   print("     energy: ")
   println(@test energy(V, at) ≈ naive_energy(V, at) )
   print("site-energy: ")
   println(@test energy(V, at) ≈ sum( site_energy(V, at, n)
                                         for n = 1:length(at) ) )
   print("     forces: ")
   println(@test JuLIP.Testing.fdtest(V, at; verbose=false))
   print("site-forces: ")
   println(@test JuLIP.Testing.fdtest( x -> site_energy(V, set_dofs!(at, x), 3),
                                       x -> mat(site_energy_d(V, set_dofs!(at, x), 3))[:],
                                       dofs(at); verbose=false ) )
end

#---

end


text:
struct ThreeDisplay <: AbstractPlotting.AbstractScreen
    session::JSServe.Session
end

JSServe.session(td::ThreeDisplay) = td.session

# We use objectid to find objects on the js side
js_uuid(object) = string(objectid(object))

function Base.insert!(td::ThreeDisplay, scene::Scene, plot::AbstractPlot)
    plot_data = serialize_plots(scene, [plot])
    WGL.insert_plot(td.session, js_uuid(scene), plot_data)
    return
end

function Base.delete!(td::ThreeDisplay, scene::Scene, plot::AbstractPlot)
    uuids = js_uuid.(AbstractPlotting.flatten_plots(plot))
    WGL.delete_plots(td.session, js_uuid(scene), uuids)
    return
end

function all_plots_scenes(scene::Scene; scene_uuids=String[], plot_uuids=String[])
    push!(scene_uuids, js_uuid(scene))
    for plot in scene.plots
        append!(plot_uuids, (js_uuid(p) for p in AbstractPlotting.flatten_plots(plot)))
    end
    for child in scene.children
        all_plots_scenes(child, plot_uuids=plot_uuids, scene_uuids=scene_uuids)
    end
    return scene_uuids, plot_uuids
end

"""
    find_plots(td::ThreeDisplay, plot::AbstractPlot)

Gets the ThreeJS object representing the plot object.
"""
function find_plots(td::ThreeDisplay, plot::AbstractPlot)
    return find_plots(JSServe.session(td), plot)
end

function find_plots(session::Session, plot::AbstractPlot)
    uuids = js_uuid.(AbstractPlotting.flatten_plots(plot))
    return WGL.find_plots(session, uuids)
end

function JSServe.print_js_code(io::IO, plot::AbstractPlot, context)
    uuids = js_uuid.(AbstractPlotting.flatten_plots(plot))
    JSServe.print_js_code(io, js"$(WGL).find_plots($(uuids))", context)
end

function three_display(session::Session, scene::Scene)
    if TEXTURE_ATLAS_CHANGED[]
        JSServe.update_cached_value!(session, AbstractPlotting.get_texture_atlas().data)
        TEXTURE_ATLAS_CHANGED[] = false
    end
    serialized = serialize_scene(scene)
    JSServe.register_resource!(session, serialized)
    window_open = scene.events.window_open

    width, height = size(scene)

    canvas = DOM.um("canvas", width=width, height=height, tabindex="0")
    wrapper = DOM.div(canvas)
    comm = Observable(Dict{String,Any}())
    push!(session, comm)

    scene_data = Observable(serialized)

    canvas_width = lift(x -> [round.(Int, widths(x))...], pixelarea(scene))

    scene_id = objectid(scene)
    setup = js"""
    function setup(scenes){
        const canvas = $(canvas)

        const scene_id = $(scene_id)
        const renderer = $(WGL).threejs_module(canvas, $comm, $width, $height)
        if ( renderer ) {
            const three_scenes = scenes.map(x=> $(WGL).deserialize_scene(x, canvas))
            const cam = new $(THREE).PerspectiveCamera(45, 1, 0, 100)
            $(WGL).start_renderloop(renderer, three_scenes, cam)
            JSServe.on_update($canvas_width, w_h => {
                renderer.setSize(w_h[0], w_h[1]);
                canvas.style.width = w_h[0];
                canvas.style.height = w_h[1];
            })
        } else {
            const warning = $(WEBGL).getWebGLErrorMessage();
            $(wrapper).removeChild(canvas)
            $(wrapper).appendChild(warning)
        }
    }
    """

    onjs(session, scene_data, setup)
    mousedrag(scene, nothing)
    scene_data[] = scene_data[]
    connect_scene_events!(scene, comm)
    mousedrag(scene, nothing)
    three = ThreeDisplay(session)

    on(session.on_close) do closed
        if closed
            scene_uuids, plot_uuids = all_plots_scenes(scene)
            WGL.delete_scenes(session, scene_uuids, plot_uuids)
            window_open[] = false
        end
    end

    return three, wrapper
end


text:
@inline function fhir_to_julia(fhirname::Symbol)
    if fhirname == :end
        return :end_fhir
    elseif fhirname == :function
            return :function_fhir
    elseif fhirname == :global
            return :global_fhir
    elseif fhirname == :import
            return :import_fhir
    elseif fhirname == :for
            return :for_fhir
    else
        return fhirname
    end
end

@inline function julia_to_fhir(julianame::Symbol)
    if julianame == :end_fhir
        return :end
    elseif julianame == :function_fhir
            return :function
    elseif julianame == :global_fhir
            return :global
    elseif julianame == :import_fhir
            return :import
    elseif julianame == :for_fhir
            return :for
    else
        return julianame
    end
end


text:
@testset "Space" begin
    # basic
    @test (@inferred length(Space(3,3))) == 2
    @test (@inferred length(Space(Symmetry(3,3),3))) == 2
    @test (@inferred Tensorial.tensorsize(Space(Symmetry(3,3),3))) == (3,3,3)
    @test (@inferred Tensorial.tensororder(Space(Symmetry(3,3),3))) == 3
    @test (@inferred Tuple(Space(Symmetry(3,3),3))) == (Symmetry(3,3),3)
    @test Space(3,2)[1] == 3
    @test Space(3,2)[2] == 2
    # prohibited
    @test_throws Exception size(Space(Symmetry(3,3),3))
    @test_throws Exception ndims(Space(Symmetry(3,3),3))
    # promotion
    @test (@inferred Tensorial.promote_space(Space(3,2), Space(3,2))) == Space(3,2)
    @test (@inferred Tensorial.promote_space(Space(3,3,3), Space(Symmetry(3,3),3))) == Space(3,3,3)
    @test (@inferred Tensorial.promote_space(Space(Symmetry(3,3,3)), Space(Symmetry(3,3),3))) == Space(Symmetry(3,3),3)
    @test (@inferred Tensorial.promote_space(Space(Symmetry(3,3,3)), Space(3,Symmetry(3,3)))) == Space(3,Symmetry(3,3))
    @test (@inferred Tensorial.promote_space(Space(Symmetry(3,3,3,3)), Space(3,Symmetry(3,3),3))) == Space(3,Symmetry(3,3),3)
    # wrong promotion
    @test_throws Exception Tensorial.promote_space(Space(3,2), Space(3,1))
    @test_throws Exception Tensorial.promote_space(Space(3,2), Symmetry(Space(3,3)))
end


text:
const resolutions = Dict(2.5 => "2.5", 5.0 => "5", 10.0 => "10")

# Interface methods

function download_raster(T::Type{WorldClim{BioClim}}; layer::Integer=1, resolution::AbstractFloat=10.0)
    1 ≤ layer ≤ 19 || throw(ArgumentError("The layer must be between 1 and 19"))

    raster_path = rasterpath(T, layer, resolution)
    zip_path = zippath(T, layer, resolution)

    if !isfile(raster_path)
        _maybe_download(zipurl(T, layer, resolution), zip_path)
        mkpath(dirname(raster_path))
        raster_name = rastername(T, layer, resolution)
        zf = ZipFile.Reader(zip_path)
        write(raster_path, read(_zipfile_to_read(raster_name, zf)))
        close(zf)
    end
    return raster_path
end

# BioClim layers don't get their own folder
rasterpath(T::Type{<:WorldClim{BioClim}}, layer) = rasterpath(T)

rastername(T::Type{<:WorldClim{BioClim}}, key, res::AbstractFloat) =
    "wc2.1_$(resolutions[res])m_bio_$(key).tif"

zipname(T::Type{<:WorldClim{BioClim}}, key, res::AbstractFloat) =
    "wc2.1_$(resolutions[res])m_bio.zip"

zipurl(T::Type{<:WorldClim{BioClim}}, key, res::AbstractFloat) =
    joinpath(WORLDCLIM_URL, "base", zipname(T, key, res))

zippath(T::Type{<:WorldClim{BioClim}}, key, res::AbstractFloat) =
    joinpath(rasterpath(T), "zips", zipname(T, key, res))


text:
# This file was generated, do not modify it. # hide
@show :popfirst! => popfirst!(res)
push!(res, 7, 0.1)
@show :push! => res
@show :pop! => pop!(res)
res

max_stars_repo_name: sadit/SimilaritySearchDemos

max_issues_repo_name: sadit/SimilaritySearchDemos

max_forks_repo_name: sadit/SimilaritySearchDemos


text:
__precompile__()
module ReadGlobal

using DelimitedFiles

export readglobal, getdimsize, readpadded, readpadded!, readfield, readfield!, checkinput, getnfilter, doinchunks, read_info, readcsv, testinput

function findglobal()
    filename="global"
    if !isfile(filename)
        for i = 1:10
            filename = "../" * filename
            isfile(filename) && break
        end
    end
    isfile(filename) || throw(error("Could not find global file"))
    return readdlm(filename,String,comment_char='/')
end

function readglobal()
    text = findglobal()
    dict = Dict{Symbol,String}()
    for i in 1:size(text)[1]
        dict[Symbol(text[i,1])] = text[i,2]
    end
    return dict
end

function getdimsize()
    dict = readglobal()
    nx = parse(Int,dict[:nx])
    ny = parse(Int,dict[:ny])
    nz = parse(Int,dict[:nz])
    x  = parse(Float64,dict[:xDomainSize])
    y  = parse(Float64,dict[:yDomainSize])
    z  = parse(Float64,dict[:zDomainSize])

    return nx,ny,nz,x,y,z
end

function readpadded!(stream::IO, field::AbstractArray{T,3}) where {T}
    dim = size(field)
    nx,ny,nz = dim
    sb = sizeof(T)
    npad = iseven(nx) ? 2 : 1
    ind = LinearIndices(dim)
    for k in 1:nz
        for j in 1:ny
            unsafe_read(stream,Ref(field,ind[1,j,k]),nx*sb)
            skip(stream, npad*sb)
        end
    end
    return field
end

function readpadded!(file::AbstractString, field::AbstractArray{T,N}) where {T,N}
    open(file) do io
        dims = size(field)
        ndim = [dims...]
        ndim[1] += iseven(ndim[1]) ? 2 : 1
        filesize(file) !== prod(ndim)*sizeof(T) && error("Dimension Mismatch")
        return readpadded!(io,field)
    end
end

function readpadded(stream,T::DataType,dims)
    field = Array{T}(undef,dims)
    return readpadded!(stream,field)
end

function readpadded(stream,dims)
    field = Array{Float64}(undef,dims)
    return readpadded!(stream,field)
end

function testinput(filename::String)
    nx,ny,nz,_,_,_ = getdimsize()
    sizefile = filesize(filename)
    if sizefile == (nx+2)*ny*nz*8
        t = true
    elseif sizefile == nx*ny*nz*8
        t = true
    elseif sizefile == (nx+2)*ny*nz*4
        t = true
    elseif sizefile == nx*ny*nz*4
        t = true
    else
        t = false
    end
    return t
end

function checkinput(filename::String,nx::Int,ny::Int,nz::Int)
    sizefile = filesize(filename)
    if sizefile == (nx+2)*ny*nz*8
        dtype = Float64
        padded = true
    elseif sizefile == nx*ny*nz*8
        dtype = Float64
        padded = false
    elseif sizefile == (nx+2)*ny*nz*4
        dtype = Float32
        padded = true
    elseif sizefile == nx*ny*nz*4
        dtype = Float32
        padded = false
    else
        error("""Input file "$(filename)" not matching number of gridpoints on global file (nx:$nx,ny:$ny,nz:$nz)""")
    end
    return dtype,padded
end

function readfield(filename::String,nx::Int,ny::Int,nz::Int)
    dtype,padded = checkinput(filename,nx,ny,nz)
    return padded ? readpadded(filename,dtype,(nx,ny,nz)) : read!(filename,Array{dtype}(undef,nx,ny,nz))
end

function readfield(filename::String)
    nx,ny,nz,x,y,z = getdimsize()
    return readfield(filename,nx,ny,nz)
end

function readfield!(filename::String,field::DenseArray)
    dtype,padded = checkinput(filename,size(field)...)
    return padded ? readpadded!(filename,field) : read!(filename,field)
end

function getnfilter()
    place = split(pwd(),"/")
    if occursin(r"^Model_",place[end])
        N = place[end-1][2:end]
        Fil = place[end-2][1]
    else
        N = place[end][2:end]
        Fil = place[end-1][1]
    end
    return N, Fil
end

function best_nc(ni::Int,no::Int,nx::Int,ny::Int,nz::Int)
    nc::Int = 1
    totalmem::Int = floor(Sys.total_memory())
    use::Float64 = 0.75
    usemem = use*totalmem
    totalfields::Int = (ni+no)*nx*ny*nz*sizeof(Float64)
    if usemem < totalfields
        newsize = totalfields
        while usemem < newsize
            nc *= 2
            newsize /= 2
        end
    end
    return nc
end

init(T,L) = T(L...)
init(::Type{Array{T,N}},L) where {T,N} = Array{T,N}(undef,L...)

function doinchunks(func::Function,nc::Int=0;input::NTuple{Ni,Pair}=(),output::NTuple{No,Pair}=()) where {Ni,No}

    nx, ny, nz, xs, ys, zs = getdimsize()

    TNi = sum(length.(getfield.(input,:second)))
    TNo = sum(length.(getfield.(output,:second)))

    nc == 0 && (nc = best_nc(TNi,TNo,nx,ny,nz))

    ifnames = getfield.(input,:second)
    ofnames = getfield.(output,:second)

    @assert mod(nx*ny*nz,nc)==0
    chunk = Int(nx*ny*nz/nc)

    iarrays = init.(getfield.(input,:first),chunk)
    oarrays = init.(getfield.(output,:first),chunk)

    calculation(func,nc,iarrays,oarrays,ifnames,ofnames)

    return 0
end

mywrite(f,v) = write(f,v)
myread!(f,v) = read!(f,v)
mywrite(f::Tuple{<:Union{<:IO,<:AbstractString}},v::Array{T,N}) where {T,N} = write(f[1],v)
myread!(f::Tuple{<:Union{<:IO,<:AbstractString}},v::Array{T,N}) where {T,N} = read!(f[1],v)

function calculation(func::F,nc,iarrays,oarrays,ifn,ofn) where {F<:Function}

    ifiles = (x->open.(x,"r")).(ifn)
    ofiles = (x->open.(x,"w")).(ofn)

    for j=1:nc
        for i in Base.OneTo(length(iarrays))
            myread!(ifiles[i],iarrays[i])
        end

        func(iarrays...,oarrays...)

        for i in Base.OneTo(length(oarrays))
            mywrite(ofiles[i],oarrays[i])
        end

        #gc()

    end

    for file in ifiles
        close.(file)
    end

    for file in ofiles
        close.(file)
    end

end

function read_info(stream)
    med = parse(Float64,readline(stream)[7:end])
    std = parse(Float64,readline(stream)[5:end])
    maxv = parse(Float64,split(readline(stream))[2])
    minv = parse(Float64,split(readline(stream))[2])
    return med,std,maxv,minv
end

read_info(stream::AbstractString) = (open(stream) do f; return read_info(f);end)

function readcsv(filename::AbstractString)
    open(filename) do f
        s = Symbol.(split(readline(f),','))
        arrs = Vector{Float64}[Vector{Float64}() for i in s]
        vals = (x->parse.(Float64,x)).(split.(readlines(f),','))
        for a in vals
            for i in eachindex(a)
               push!(arrs[i],a[i])
            end
        end
        #return NamedTuple{(s...,)}((arrs...,))
        return FakeNamedTuple(s,arrs)
    end
end

struct FakeNamedTuple
    names::Vector{Symbol}
    data::Vector{Vector{Float64}}
end

function Base.getproperty(a::FakeNamedTuple,x::Symbol)
    i = findfirst(y->x==y,getfield(a,:names))
    i === nothing && return error("""No field named "$x" """)
    return getfield(a,:data)[i]
end

end # module


text:
using SimpleProbabilitySets
using Base.Test, Distributions

@testset "SimpleProbabilitySets" begin
    d = Categorical([0.1, 0.2, 0.7])
    d2 = Categorical([0.4, 0.4, 0.2])
    pb = PBox(d, d2)
    @test cdfs(pb)[2] == [0.4, 0.8, 1.0]
    @test pints(pb)[1] ≈ [0.1, 0.0, 0.2] atol = 1e-6

    pl = [0.1, 0.2, 0.5]
    pu = [0.4, 0.4, 0.7]
    uncsmall = 0.05
    unclarge = 0.6
    pi = PInterval(pl, pu)
    pi2 = PInterval(pl, uncsmall)
    pi3 = PInterval(pl, unclarge)
    @test pi.plower == [0.1, 0.2, 0.5]
    @test pi2.plower ≈ [0.05, 0.15, 0.45] atol = 1e-6
    @test pi3.pupper ≈ [0.7, 0.8, 1.0] atol = 1e-6

    p = psample(pi)
    @test sum(p) ≈ 1.0 atol = 1e-9
    @test all(pi.plower .< p .< pi.pupper)
end


text:
"""
Non-linear grid.

- `nonlinear = 1` should make it uniform
- `nonlinear > 1` concentrates around xlow
- `nonlinear < 1` concentrates around xhigh

"""
function gridpoints(xlow, xhigh, n, nonlinear)
    n < 3 && error("Use more points for irregular grid")

    if nonlinear != 1
        x = fill(0., n)
        x[1] = xlow
        for i in 2:(n-1)
            x[i]=x[i-1] + (xhigh - x[i-1])/(n-i+1)^nonlinear
        end
        x[n]=xhigh
        return x
    end

    return collect(range(xlow, xhigh, length=n))
end


text:
@testset "_simplify_angle" begin
    # test typical ouput types
    @test typeof(Augmentor._simplify_angle(Float32(0.5))) <: Float32
    @test typeof(Augmentor._simplify_angle(Float64(0.5))) <: Float64
    @test typeof(Augmentor._simplify_angle(Int32(1))) <: Float64
    @test typeof(Augmentor._simplify_angle(Int64(1))) <: Float64

    # test if returntype of +angle is the same as returntype of -angle
    for i in (Float32(0.5), Float64(0.5), Int32(1), Int64(1))
        @test typeof(Augmentor._simplify_angle(i)) <: typeof(Augmentor._simplify_angle(-i))
    end

    # test hardcoded numbers (edge chases)
    @test Augmentor._simplify_angle(0) == 0
    @test Augmentor._simplify_angle(2pi) == 0
    @test Augmentor._simplify_angle(-2pi) == 0
    @test Augmentor._simplify_angle(3pi) == 1pi
    @test Augmentor._simplify_angle(-3pi) == 1pi
    @test Augmentor._simplify_angle(4pi) == 0
    @test Augmentor._simplify_angle(-4pi) == 0

    for angle in 0:.4:4pi
        @test Augmentor._simplify_angle(angle) == angle % 2pi
        if angle > 0
            @test Augmentor._simplify_angle(-angle) == 2pi + (-angle % 2pi)
        end
    end
end

@testset "_expand_size" begin
    # test type stability
    for w in (20, 20.),
        h in (10, 10.),
        angle = (Int16(1), Int32(1), Int64(1), Float16(.5), Float32(.5), Float64(.5))
        @test typeof(Augmentor._expand_size(w, h, angle)) <: Tuple{Int64,Int64}
    end

    # test output sizes for trivial input
    for w in 10:500:1000, h in 10:500:1000
        @test Augmentor._expand_size(w, h, 0) == (w, h)
        @test Augmentor._expand_size(w, h, 1pi) == (w, h)
        @test Augmentor._expand_size(w, h, pi/2) == (h, w)
        @test Augmentor._expand_size(w, h, 3pi/2) == (h, w)
    end

    # test output sized for some known inputs
    @test Augmentor._expand_size(768, 512, deg2rad(30)) == (921, 827)
    @test Augmentor._expand_size(768, 512, deg2rad(45)) == (905, 905)
    @test Augmentor._expand_size(768, 512, deg2rad(120)) == (827, 921)
    @test Augmentor._expand_size(768, 512, deg2rad(190)) == (845, 637)
    @test Augmentor._expand_size(768, 512, deg2rad(300)) == (827, 921)
end

A = UInt8[200 150; 50 1]
img_x = grayim(A)
img_y = permutedims(img_x, [2, 1])

@testset "rotate_expand" begin
    # test types stability
    @test typeof(rotate_expand(A, deg2rad(0))) <: typeof(A)
    @test typeof(rotate_expand(img_x, deg2rad(0))) <: typeof(img_x)
    @test typeof(rotate_expand(A, deg2rad(90))) <: typeof(A)
    @test typeof(rotate_expand(img_x, deg2rad(90))) <: typeof(img_x)
    @test typeof(rotate_expand(A, deg2rad(180))) <: typeof(A)
    @test typeof(rotate_expand(img_x, deg2rad(180))) <: typeof(img_x)
    @test typeof(rotate_expand(A, deg2rad(270))) <: typeof(A)
    @test typeof(rotate_expand(img_x, deg2rad(270))) <: typeof(img_x)

    # test values for some sample matrices
    @test rotate_expand(A, deg2rad(0)) == A
    @test rotate_expand(A, deg2rad(90)) == rotr90(A)
    @test rotate_expand(A, deg2rad(180)) == rot180(A)
    @test rotate_expand(A, deg2rad(270)) == rotl90(A)

    # test values for some sample images
    @test raw(rotate_expand(img_x, deg2rad(0))) == A
    @test raw(rotate_expand(img_y, deg2rad(0))) == A'
    @test raw(rotate_expand(img_x, deg2rad(90))) == rotr90(A)
    @test raw(rotate_expand(img_y, deg2rad(90))) == rotl90(A')
    @test raw(rotate_expand(img_x, deg2rad(270))) == rotl90(A)
    @test raw(rotate_expand(img_y, deg2rad(270))) == rotr90(A')
    @test raw(rotate_expand(img_x, deg2rad(180))) == rot180(A)
    @test raw(rotate_expand(img_y, deg2rad(180))) == rot180(A')

    @imagetest "rotate_expand_45" rotate_expand(testimg, deg2rad(45))
    @imagetest "rotate_expand_120" rotate_expand(testimg, deg2rad(120))
    @imagetest "rotate_expand_220" rotate_expand(testimg, deg2rad(220))
    @imagetest "rotate_expand_320" rotate_expand(testimg, deg2rad(320))
end

@testset "rotate_crop" begin
    #pending "test types for some sample matrices"
    @test typeof(rotate_crop(A, deg2rad(0))) <: typeof(A)
    @test typeof(rotate_crop(img_x, deg2rad(0))) <: typeof(img_x)
    @test typeof(rotate_crop(A, deg2rad(180))) <: typeof(A)
    @test typeof(rotate_crop(img_x, deg2rad(180))) <: typeof(img_x)

    #pending "test values for some sample matrices"
    @test rotate_crop(A, deg2rad(0)) == A
    @test rotate_crop(A, deg2rad(180)) == rot180(A)

    #pending "test values for some sample images"
    @test raw(rotate_crop(img_x, deg2rad(0))) == A
    @test raw(rotate_crop(img_y, deg2rad(0))) == A'
    @test raw(rotate_crop(img_x, deg2rad(180))) == rot180(A)
    @test raw(rotate_crop(img_y, deg2rad(180))) == rot180(A')

    @imagetest "rotate_crop_45" rotate_crop(testimg, deg2rad(45))
    @imagetest "rotate_crop_120" rotate_crop(testimg, deg2rad(120))
    @imagetest "rotate_crop_220" rotate_crop(testimg, deg2rad(220))
    @imagetest "rotate_crop_320" rotate_crop(testimg, deg2rad(320))
end


text:
#=
Using a function rand5() that returns an integer from 1 to 5 (inclusive) with uniform probability, implement a function rand7() that returns an integer from 1 to 7 (inclusive).
=#

using Test
include("Solutions/problem45_rand7_from_rand5.jl")

@test rand7() in collect(1:7)

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
# Utilities to rewrite UIDs in .squashfs files.

"""
    getuid()

Wrapper around libc's `getuid()` function
"""
function getuid()
    return ccall(:getuid, Cint, ())
end

"""
    getgid()

Wrapper around libc's `getgid()` function
"""
function getgid()
    return ccall(:getgid, Cint, ())
end

# Note: produce these values by #including squashfs_fs.h from linux in Cxx.jl
# and running the indicated command
const offsetof_id_table_start = 0x30    # offsetof(struct suqashfs_super_block, id_table_start)
const offsetof_no_ids = 0x1a            # offsetof(struct suqashfs_super_block, no_ids)

# From squashfs_fs.h
const SQUASHFS_COMPRESSED_BIT = UInt16(1) << 15
const SQUASHFS_MAGIC = 0x73717368

"""
    generate_per_uid_squashfs(cs, new_uid = getuid())

In order for the sandbox to work well, we need to have the uids of the squashfs
images match the uid of the current unpriviledged user.  Unfortunately there is
no mount-time option to do this for us.  Fortunately, squashfs is simple enough
that if the ID table is uncompressed, we can just manually patch the uids to be
what we need.  This function performs this operation, by rewriting all UIDs and
GIDs to the given `new_uid` (which defaults to the current user's UID)..
"""
function generate_per_uid_squashfs(cs, new_uid = getuid(); verbose::Bool = false)
    # Because the source .squashfs file comes as an immutable artifact, we will create
    # a new artifact that is based off of the progenitor .squashfs artifact, but with
    # rewritten UIDs.  We do this by maintaining a separate `MutableArtifacts.toml`
    # file that is explicitly `.gitignore`'d.
    artifacts_toml = joinpath(dirname(@__DIR__), "Artifacts.toml")
    name = artifact_name(cs)
    progenitor_hash = artifact_hash(name, artifacts_toml; platform=cs.host)
    if progenitor_hash == nothing
        error("Compiler shard $(name) not found in $(artifacts_toml)")
    end

    # We now generate a name that will be bound within `MutableArtifacts.toml`.
    # This name encodes not only which artifact we're using, but also the source hash
    # and the desired UID.  This allows us to intelligently cache output from a name
    # that is tied to a specific version of the input, and all build options.
    cache_name = "$(name)_$(bytes2hex(progenitor_hash.bytes[end-8:end]))_$(new_uid)"

    # If this name isn't bound to anything, then create it by remapping the UIDs
    if get_mutable_artifact_path(cache_name) === nothing
        # Grab the progenitor path, downloading it if necessary
        progenitor_path = joinpath(artifact_path(progenitor_hash), name)

        if !isfile(progenitor_path)
            error("Compiler shard $(name) missing from disk at $(progenitor_path)")
        end

        # Create a new artifact that has the modified `.squashfs` file within it
        create_and_bind_mutable_artifact!(cache_name) do dir
            # Copy .squashfs file over to our local directory, make it writable
            cp(progenitor_path, joinpath(dir, name))
            chmod(joinpath(dir, name), 0o644)

            squash_path = joinpath(dir, name)
            open(squash_path, "r+") do file
               # Check magic
               if read(file, UInt32) != SQUASHFS_MAGIC
                   error("`$squash_path` is not a squashfs file")
               end
               # Check that the image contains only one id (which we will rewrite)
               seek(file, offsetof_no_ids)
               if read(file, UInt16) != 1
                   error("`$squash_path` uses more than one uid/gid")
               end
               # Find the index table
               seek(file, offsetof_id_table_start)
               offset = read(file, UInt64)
               seek(file, offset)
               # Find the correct metdata block
               index = read(file, UInt64)
               seek(file, index)
               # Read the metadata block
               size = read(file, UInt16)
               # Make sure it's uncompressed (yes, I know that flag is terribly
               # named - it indicates that the data is uncompressed)
               if ((size & SQUASHFS_COMPRESSED_BIT) == 0)
                   error("Metadata block is compressed")
               end
               p = position(file)
               uid = read(file, UInt32)
               if uid == new_uid
                   return
               end
               if verbose
                   @info("Rewriting $(basename(squash_path)) from UID $(uid) -> $(new_uid)")
               end
               seek(file, p)
               write(file, UInt32(new_uid))
            end
        end
    end

    # Finally, get the path to that cached artifact
    return joinpath(get_mutable_artifact_path(cache_name), name)
end


text:
# coding: utf-8

# In[1]:

#radix sort is fun
#there are two ways to do it
#the first approach is we start from the right
#the second one is to start from the left
#in this context, i choose the first approach
#they are equivalent
#we get the first digit of all elements
#sort these elements based upon the first digit
#for instance, 215 and 19
#9 for 19 is larger than 5 for 215
#for the first round, we get 215,19
#we move on to the second digit of all elements from the right
#we apply the same technique of sorting based upon the second digit
#until we reach to the very left digit of the maximum value in the list
#the sort is completed
#for instance, the second round would be 215,19
#since the second digit of both numbers are 1
#the final round is 19,215
#as 0 for 19 is smaller than 2 for 215
#period

#to implement this method in python
#we can apply fld(x,10)%10 technique
#for the first round, fld(x,1)%10 gets the first digit from the left
#for the second round, fld(x,10)%10 gets the second digit from the left
#and so on
#we can use a nested list to store elements in its rightful place
#i use digits its index indicates its digit+1
#becuz julia starts from 1
#there are other elegant solutions, plz check the link at the very end

# In[2]:

function radix_sort(arr)

    maxval=maximum(arr)
    digits=0

    #the first part is to find the largest element in the list
    #next, we find out how many digits it contains
    #so we can determine how many times we need to divide in the later stage
    while maxval>=1
        maxval=maxval/10
        digits+=1
    end

    for digit in 1:digits

        #each time, we clear the nested list
        digit_sort=[[] for _ in 1:10]
        radix=10^(digit-1)

        #this is a traversal on all elements in the list
        for i in arr

            #fld(i,radix)%10 is the technique to get digits
            #when we get it, we use i+1 as index
            #so that we can store the element in its rightful place
            ind=fld(i,radix)%10

            push!(digit_sort[ind+1],i)

        end

        #we have to empty the list after we put all elements into a nested list
        arr=[]

        #this loop is unnecessary
        #we can use append(digit_sort[1],...digit_sort[10])
        #again, i am a lazy person
        #cut me some slack, s’il vous plaît
        for j in 1:10
            append!(arr,digit_sort[j])
        end

    end

    return arr

end

# In[3]:

#surprisingly, the radix sort is really fast
#the time complexity of radix sort is o(k*n)
#k depends on the digits of the maximum value in the list
#well, its harder to get a number like 10^1000 than getting a list with 10^1000 elements
#when we get a list with smaller n and larger k
#radix sort may be a terrible choice

#there are other ways of writing radix sort
#i believe my way has too much space complexity
#this one is rather smart and elegant
# https://www.geeksforgeeks.org/radix-sort/

# In[4]:

for _ in 1:100

    test_arr=rand(1:1000,100)

    if !(radix_sort(test_arr)==sort(test_arr))
        printstyled("Erreur",color=:red)
    end

end

max_stars_repo_path: radix sort.jl

max_stars_repo_name: je-suis-tm/search-and-sort

max_issues_repo_path: radix sort.jl

max_issues_repo_name: HafizSabia/search-and-sort

max_forks_repo_path: radix sort.jl

max_forks_repo_name: HafizSabia/search-and-sort


text:
module Plot

using Gadfly
using Query
using DataFrames
using GeCKOReader

export plot_all_qualities,
       plot_hit_qualities,
       plot_guide_trajectories

function _plot(s::Matrix{Int})
    rows = Int[]
    trim_result = []
    for row in 1:size(s)[1]
        if sum(s[row, :] .!= 0) > 0
            push!(trim_result, s[row, :])
            push!(rows, row)
        end
    end
    trim_result = transpose(hcat(trim_result...))
    spy(trim_result, Scale.y_discrete(labels = i ->rows[i]),
        Guide.xticks(ticks=collect(1:3:size(trim_result)[2])))
end

"""
    plot_all_qualities(s::ScreenSample)

Generate a 2D histogram of read qualities from the FASTQ file per position in
the read using **all** of the read info.
"""
plot_all_qualities(s::GeCKOReader.ScreenSample) = _plot(s.all_qualities)

"""
    plot_hit_qualities(s::ScreenSample)

Generate a 2D histogram of read qualities from the FASTQ file per position in
the read using **only reads that were matched**.
"""
plot_hit_qualities(s::GeCKOReader.ScreenSample) = _plot(s.hit_qualities)

"""
    plot_guide_trajectories(lines::DataFrame, gene::String)

Plot the trajectories of the guides for a specific `gene` given the data in
`lines`. The `lines` dataframe should be formatted like the output of
[`fit_lines`](@ref Main.GeCKOReader.fit_lines).
"""
function plot_guide_trajectories(lines::DataFrame, gene::String)

    (!([:t0_rel_freqs, :t1_rel_freqs, :t2_rel_freqs,
        :intercept, :slope, :id, :gene] ⊆ names(lines))) && error("DataFrame malformed missing required columns")

    df = @from i in lines begin
        @where i.gene == gene
        @select i
        @collect DataFrame
    end

    lens = size(df, 1)
    tmp = DataFrame(x=repeat([0,1,2], inner=lens),
                    y=vcat(df[:t0_rel_freqs], df[:t1_rel_freqs], df[:t2_rel_freqs]),
                    guides=repeat(df[:id], outer=3))

    xs = repeat([0, 2], inner=lens)
    tmp2 = DataFrame(x= xs,
                     y=vcat(df[:intercept], df[:intercept] + df[:slope].*2),
                     guides=repeat(df[:id], outer=2))

    plot(
        layer(tmp, x=:x, y=:y, color=:guides, Geom.point),
        layer(tmp2, x=:x, y=:y, color=:guides, Geom.line),
        Guide.xlabel("# of osmoshock treatments"),
        Guide.ylabel("Relative frequencies of guide"),
        Guide.title(gene),
        Guide.xticks(ticks=[0, 1, 2])
        )
end

end #module


text:
function common_subexpressions(n)

  reps = parse(Int64, ENV["LAMP_REPS"])

  A = randn(n,n)
  B = randn(n,n)

  Benchmarker.add_data(csv, "subexpr_nai", Benchmarker.measure(reps, implicit_subexpression, A, B))
  Benchmarker.add_data(csv, "subexpr_rec", Benchmarker.measure(reps, explicit_subexpression, A, B))

  @test isapprox(implicit_subexpression(A, B)[2], explicit_subexpression(A, B)[2])
end

function implicit_subexpression(A::DenseArray{Float64, 2}, B::DenseArray{Float64, 2})
  K = zeros(size(A))
  start::Float64 = 0.0
  finish::Float64 = 0.0
  Benchmarker.cachescrub()
  Benchmarker.gcscrub()
  start = time_ns()

  K = A * B + A * B

  finish = time_ns()
  return ((finish-start)*1e-9, K)
end

function explicit_subexpression(A::DenseArray{Float64, 2}, B::DenseArray{Float64, 2})
  K = zeros(size(A))
  start::Float64 = 0.0
  finish::Float64 = 0.0
  Benchmarker.cachescrub()
  Benchmarker.gcscrub()
  start = time_ns()

  K = A * B
  K = K + K

  finish = time_ns()
  return ((finish-start)*1e-9, K)
end


text:
abstract type Component end

include("objects.jl")
include("elements.jl")

include("cache.jl")

mutable struct OrgDoc
    settings::Dict{AbstractString, Any}
    contents::Vector{Union{Heading, Section}}
    cache::OrgCache
    function OrgDoc(settings::Dict, contents::Vector{<:Union{Heading, Section}})
        o = new(settings, contents, OrgCache())
        setfield!(o.cache, :doc, o)
        o
    end
end

OrgDoc(contents::Vector{<:Union{Heading, Section}}) =
    OrgDoc(Dict{AbstractString, Any}(), contents)
OrgDoc() = OrgDoc(Union{Heading, Section}[])


text:
using StochasticDiffEq, DiffEqProblemLibrary, Test, Random
using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
import DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_linear_stratonovich, prob_sde_2Dlinear_stratonovich
Random.seed!(100)
dts = 1 ./2 .^(10:-1:2) #14->7 good plot

prob = prob_sde_linear_stratonovich
sim  = test_convergence(dts,prob,EulerHeun(),numMonte=Int(5e2))
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim  = test_convergence(dts,prob,LambaEulerHeun(),numMonte=Int(5e2))
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim  = test_convergence(dts,prob,ISSEulerHeun(),numMonte=Int(5e2))
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim  = test_convergence(dts,prob,ImplicitEulerHeun(),numMonte=Int(5e2))
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim  = test_convergence(dts,prob,ImplicitEulerHeun(theta=1),numMonte=Int(5e2))
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim  = test_convergence(dts,prob,ImplicitEulerHeun(symplectic=true),numMonte=Int(5e2))
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim = test_convergence(dts,prob,RKMil(interpretation=:Stratonovich),numMonte=Int(5e2))
@test abs(sim.𝒪est[:l2]-1) < 0.2

sim  = test_convergence(dts,prob,ImplicitRKMil(interpretation=:Stratonovich),numMonte=Int(5e2))
@test abs(sim.𝒪est[:l2]-1) < 0.1

println("Now 2D")

prob = prob_sde_2Dlinear_stratonovich

sim  = test_convergence(dts,prob,EulerHeun(),numMonte=Int(5e1))
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim  = test_convergence(dts,prob,LambaEulerHeun(),numMonte=Int(5e1))
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim  = test_convergence(dts,prob,ISSEulerHeun(),numMonte=Int(5e1))
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim  = test_convergence(dts,prob,ImplicitEulerHeun(),numMonte=Int(5e1))
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim  = test_convergence(dts,prob,ImplicitEulerHeun(theta=1),numMonte=Int(5e1))
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim  = test_convergence(dts,prob,ImplicitEulerHeun(symplectic=true),numMonte=Int(5e1))
@test abs(sim.𝒪est[:l2]-1) < 0.1

println("RKMils")

sim = test_convergence(dts,prob,RKMil(interpretation=:Stratonovich),numMonte=Int(1e2))
@test abs(sim.𝒪est[:l2]-1) < 0.2

sim  = test_convergence(dts,prob,ImplicitRKMil(interpretation=:Stratonovich),
                        numMonte=Int(1e2))
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim  = test_convergence(dts,prob,ImplicitRKMil(theta=1,interpretation=:Stratonovich),
                        numMonte=Int(1e2))
@test abs(sim.𝒪est[:l2]-1) < 0.1

sim  = test_convergence(dts,prob,ImplicitRKMil(symplectic=true,interpretation=:Stratonovich),
                        numMonte=Int(1e2))
@test abs(sim.𝒪est[:l2]-1) < 0.1

Random.seed!(200)
sol = solve(prob,EulerHeun(),dt=1/4)


text:
abstract type StepSize end

struct UserStepSize{F} <: StepSize
    f::F
    function UserStepSize(f::F) where {F}
        @assert isa(f(1), Float64) "Invalid return type, must be Float64"
        return new{F}(f)
    end
end
(s::UserStepSize)(η::Int) = s.f(max(η, 1))

struct SpecifiedStepSize <: StepSize
    s::Vector{Float64}
    function SpecifiedStepSize(a::A) where {T<:Real,A<:AbstractVector{T}}
        return new(convert(Vector{Float64}, a))
    end
end
(s::SpecifiedStepSize)(η::Int) = s.s[min(max(η, 1), length(s.s))]

struct ExpStepSize <: StepSize
    ds0::Float64
    ratio::Float64
end
(s::ExpStepSize)(η::Int) = s.ds0 * s.ratio^(max(η, 1) - 1)

struct ConstStepSize <: StepSize
    ds::Float64
end
(s::ConstStepSize)(::Int) = s.ds

struct FilterStepSize <: StepSize
    η::Int
    sa::Float64
    sb::Float64
end
(s::FilterStepSize)(η::Int) = η <= s.η ? s.sa : s.sb

struct MinStepSize <: StepSize
    children::Vector{StepSize}
    MinStepSize(a...) =  new(StepSize[StepSize(f) for f in a])
end
(s::MinStepSize)(η::Int) = minimum(map(ss->ss(η), s.children))

StepSize(a::Real) = ConstStepSize(a)
StepSize(a::A) where {A<:StepSize} = a
StepSize(a::Function) = UserStepSize(a)
StepSize(a::AbstractVector) = SpecifiedStepSize(a)

abstract type BoundaryCondition end

struct FreeBoundary <: BoundaryCondition end

struct SymmetryBoundary <: BoundaryCondition end

struct LineBoundary <: BoundaryCondition
    dir::SVector{2,Float64}
    LineBoundary(x, y) = new(normalize(SVector{2,Float64}(x, y)))
end

ConstantX() = LineBoundary(0, 1)
ConstantY() = LineBoundary(1, 0)
ConstantAngle(α) = LineBoundary(cosd(α), sind(α))

one_step(::FreeBoundary, h::Float64, d0::MVector{2,Float64}) = normalize(d0) * h
function one_step(b::LineBoundary, h::Float64, d0::MVector{2,Float64})
    return dot(b.dir, d0) < 0 ? -h * b.dir : h * b.dir
end

struct WakeInfo
    len::Float64
    ratio::Float64
    dsmax::Float64
    aoa::Float64
end


text:
using Pumas.NCA

using Pumas, CSV

file = Pumas.example_nmtran_data("nca_test_data/dapa_IV")
data = CSV.read(file);

first(data, 6) # take first 6 rows

NCA.auc(data[:CObs][1:16], data[:TIME][1:16])

NCA.auc(data[:CObs][1:16], data[:TIME][1:16], method=:linuplogdown)

timeu = u"hr"
concu = u"mg/L"
amtu  = u"mg"
data.id = data.ID
data.time = data.TIME
data.conc = data.CObs
data.amt = data.AMT_IV
data.route = "iv"
pop = read_nca(data, llq=0concu, timeu=timeu, concu=concu, amtu=amtu)

NCA.auc(pop)

NCA.auc(pop[2], auctype=:last)

NCA.auc(pop, auctype=:last)

NCA.auc(pop[1], interval=(10,Inf).*timeu)

NCA.auc(pop[1], interval=[(10,Inf).*timeu, (10, 15).*timeu])

NCA.auc_extrap_percent(pop[1])

NCA.aumc_extrap_percent(pop[1])
NCA.aumc(pop[1])

NCA.lambdaz(pop[1])

NCA.lambdazr2(pop)
NCA.lambdazadjr2(pop)
NCA.lambdazintercept(pop)
NCA.lambdaztimefirst(pop)
NCA.lambdaznpoints(pop)

NCA.lambdaz(pop[1], threshold=3)

NCA.lambdaz(pop[1], idxs=[10, 15, 16])

NCA.lambdaz(pop[1], slopetimes=[1,2,3].*timeu)

NCA.tmax(pop[1])
NCA.cmax(pop[1])
NCA.cmax(pop[1], interval=(20, 24).*timeu)
NCA.cmax(pop[1], interval=[(20, 24).*timeu, (10, 15).*timeu])

NCA.tlast(pop[1])
NCA.clast(pop[1])

NCA.thalf(pop[1])

NCA.interpextrapconc(pop[1], 12timeu, method=:linear)

using Plots # load the plotting library
plot(pop)

plot(pop, loglinear=false)

plot(pop, linear=false)

report = NCAReport(pop)

reportdf = NCA.to_dataframe(report)
first(reportdf,6) # Print only the first 6 rows

multiple_doses_file = Pumas.example_nmtran_data("nca_test_data/dapa_IV_ORAL")
mdata = CSV.read(multiple_doses_file)

timeu = u"hr"
concu = u"mg/L"
amtu  = u"mg"
mdata.id = mdata.ID
mdata.time = mdata.TIME
mdata.conc = mdata.COBS
mdata.amt = mdata.AMT
mdata.route = replace(mdata.FORMULATION, "IV"=>"iv", "ORAL"=>"ev")
mdata.occasion = mdata.OCC
mpop = read_nca(mdata, timeu=timeu, concu=concu, amtu=amtu)

plot(mpop)

df = NCA.auc(mpop)
first(df,6) # Print only the first 6 rows

rep = NCAReport(mpop, ithdose=1)
reportdf = NCA.to_dataframe(rep)
first(reportdf,6) # Print only the first 6 rows

using PumasTutorials
PumasTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
s1 = "The quick brown fox jumps over the lazy dog α,β,γ"

# [search](https://docs.julialang.org/en/v1/base/strings/#Base.findfirst-Tuple{AbstractString,AbstractString}) returns the first index of a char
i = findfirst(isequal('b'), s1)
println(i)
#> 11
# the second argument is equivalent to the second argument of split, see below

# or a range if called with another string
r = findfirst("brown", s1)
println(r)
#> 11:15

# string [replace](https://docs.julialang.org/en/v1/base/strings/#Base.replace-Tuple{AbstractString,Pair}) is done thus:
r = replace(s1, "brown" => "red")
show(r); println()
#> "The quick red fox jumps over the lazy dog α,β,γ"

# search and replace can also take a regular expressions by preceding the string with 'r'
r = findfirst(r"b[\w]*n", s1)
println(r)
#> 11:15

# again with a regular expression
r = replace(s1, r"b[\w]*n" => "red")
show(r); println()
#> "The quick red fox jumps over the lazy dog α,β,γ"

# there are also functions for regular expressions that return RegexMatch types
# [match](https://docs.julialang.org/en/v1/base/strings/#Base.match) scans left to right for the first match (specified starting index optional)
r = match(r"b[\w]*n", s1)
println(r)
#> RegexMatch("brown")

# RegexMatch types have a property match that holds the matched string
show(r.match); println()
#> "brown"

# [eachmatch](https://docs.julialang.org/en/v1/base/strings/#Base.eachmatch) returns an iterator over all the matches
r = eachmatch(r"[\w]{4,}", s1)
for i in r print("\"$(i.match)\" ") end
#> "quick" "brown" "jumps" "over" "lazy"
println()

r = collect(m.match for m = eachmatch(r"[\w]{4,}", s1))
println(r)
#> SubString{String}["quick", "brown", "jumps", "over", "lazy"]

# a string can be repeated using the [repeat](http://julia.readthedocs.org/en/latest/manual/strings/#common-operations) function,
# or more succinctly with the [^ syntax](https://docs.julialang.org/en/v1/base/strings/#Base.:^-Tuple{AbstractString,Integer}):
r = "hello "^3
show(r); println() #> "hello hello hello "

# the [strip](https://docs.julialang.org/en/v1/base/strings/#Base.strip) function works the same as python:
# e.g., with one argument it strips the outer whitespace
r = strip("hello ")
show(r); println() #> "hello"
# or with a second argument of an array of chars it strips any of them;
r = strip("hello ", ['h', ' '])
show(r); println() #> "ello"
# (note the array is of chars and not strings)

# similarly [split](https://docs.julialang.org/en/v1/base/strings/#Base.split) works in basically the same way as python:
r = split("hello, there,bob", ',')
show(r); println() #> SubString{String}["hello", " there", "bob"]
r = split("hello, there,bob", ", ")
show(r); println() #> SubString{String}["hello", "there,bob"]
r = split("hello, there,bob", [',', ' '], limit=0, keepempty=false)
show(r); println() #> SubString{String}["hello", "there", "bob"]
# (the last two arguements are limit and include_empty, see docs)

# the opposite of split: [join](https://docs.julialang.org/en/v1/base/strings/#Base.join) is simply
r = join(collect(1:10), ", ")
println(r) #> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10

max_stars_repo_name: talgatomarov/JuliaByExample

max_issues_repo_name: talgatomarov/JuliaByExample

max_forks_repo_name: talgatomarov/JuliaByExample