text:
export ARC_generic_ls

"""
A one dimensionnal ARC algorithm designed for Line Search purposes.

In this generic version we assume a second order polynomial enriched with a cubic
regularization term is used to approximate our h function.
We use c(d) = h(θ) + h'(θ)⋆d + 0.5 * h''(θ) ⋆ d² + Δ/3 |d|³.
Also h''(θ) is sometimes approximated by a secant (Sec) or an improved secant (SecA).
That approximation
is decided with hyper parameter "direction"

Some specific parameters for theses algorithms:
  - Δ :: The cubic regularization parameter
Some parameters are used for TR and ARC methods:
  -eps1 and eps2 :: treshold to determine if q is a bad or a good++ approximation
    of h
  -red and aug :: Controls the speed at which we reduce or augment the size of
    the trust region depending on the previously mentionned treshold.

All the parameters mentionned above can be found in the LS_Function_Meta
"""
function ARC_generic_ls(h         :: LineModel, stop_ls   :: LS_Stopping;
                        f_meta = LS_Function_Meta(),
                        φ_dφ      :: Function = (x, y) -> phi_dphi(x, y),
                        verboseLS :: Bool = false, kwargs...)

    state = stop_ls.current_state

    update!(state, x = 1.0, ht = obj(h, 1.0), gt = grad(h, 1.0), tmps = time())

    # Check if 1.0 is an admissible step size.
    OK = start!(stop_ls)

    # If 1.0 isn't an admissible step size we start our line search algorithm
    iter = 0
    # If 1.0 is a good step size we don't want to compute  unnecessary function
    # evaluations
    H = NaN; φt = NaN; dφt = NaN                   # only for the verbose
    if OK == false
      φt, dφt = φ_dφ(h, state)

      tprec = NaN; φtestTR = NaN; dφtestTR = NaN; φtprec = NaN; dφtprec = NaN;

      # H will denote the approximation to φ'' hereafter
      if f_meta.dir == :Nwt
        H = hess(h, state.x)
      elseif f_meta.dir == :Sec || f_meta.dir == :SecA
        H = 1.0
      end
    end # if OK == false

    verboseLS && @printf("   iter   t       φt        dφt        Δ")
    verboseLS && @printf("   Successful        d        ratio\n")
    verboseLS && @printf("%4d %9.2e  %9.2e  %9.2e %9.2e\n",
                         iter, state.x, φt, dφt, f_meta.Δ);

    while !OK # admissible: respecte armijo et wolfe,
                                # tired: nb d'itérations

        d = ARC_step_computation(H, dφt, f_meta.Δ; kwargs...)

        # We see where the direction d gets us
        candidate_state = copy(state)
        update!(candidate_state, x = state.x + d)
        φtestTR, dφtestTR = φ_dφ(h, candidate_state; kwargs...)

        # depending on the approximation of the second derivate we need different
        # information
        if f_meta.dir == :Sec || f_meta.dir == :SecA
          tprec = state.x
          dφtprec = dφt
        end
        if f_meta.dir == :SecA φtprec = φt end

        # test d'arrêt sur dφ
        (pred, ared, ratio) =
               pred_ared_computation(dφt, φt, H, d, φtestTR, dφtestTR)

        if (state.x + d < 0.0) || ratio < f_meta.eps1  # Unsuccessful
          f_meta.Δ = f_meta.red * f_meta.Δ
          iter += 1
          verboseLS &&
            @printf("%4d %9.2e %9.2e  %9.2e  %9.2e %9.2e %7.2e\n",
                    iter, state.x, φt, dφt, f_meta.Δ, 0, d);
        else             # Successful
            stop_ls.current_state = copy(candidate_state)
            state = stop_ls.current_state
            φt, dφt = φtestTR, dφtestTR
            H = update_H(f_meta.dir, h, state.x, tprec, φt, dφt, φtprec, dφtprec)

          if ratio > f_meta.eps2
            f_meta.Δ = f_meta.aug * f_meta.Δ
          end
          iter += 1
          verboseLS &&
            @printf("%4d %9.2e %9.2e  %9.2e  %9.2e  %9.2e %7.2e \n",
                    iter, state.x, φt, dφt, f_meta.Δ, 1, d);
        end;
        OK = stop!(stop_ls)
    end;

    return (state, stop_ls.meta.optimal)

end


text:
#=
equivalent code would normally be executed by Travis
make sure that this is run from pwd = GpABC.jl/docs
using Pkg
Pkg.develop(PackageSpec(path=".."))
Pkg.instantiate()
local_build = true
=#
using Documenter, GpABC

local_build = "local" in ARGS

gpabc_makedocs() = makedocs(
    modules = [GpABC],
    format = Documenter.HTML(
        prettyurls = !local_build,
        canonical = "https://tanhevg.github.io/GpABC.jl/stable/",
    ),
    sitename = "GpABC",
    authors = "Tara Abdul Hameed, Fei He, Jonathan Ish-Horowitz, Istvan Klein, Michael Stumpf, Evgeny Tankhilevich",
    pages = [
        "Home" => "index.md",
        "Package Overview" => [
            "ABC Parameter Inference" => "overview-abc.md",
            "ABC Model Selection" => "overview-ms.md",
            "LNA" => "overview-lna.md",
            "Gaussian Process Regression" => "overview-gp.md",
            "summary_stats.md"
        ],
        "Examples" => [
            "ABC Parameter Inference" => "example-abc.md",
            "ABC Model Selection" => "example-ms.md",
            "Stochastic inference (LNA)" => "example-lna.md",
            "Gaussian Processes" => "example-gp.md",
        ],
        "Reference" => [
            "ABC Basic" => "ref-abc.md",
            "ABC Advanced" => "ref-abc-advanced.md",
            "Stochastic inference (LNA)" => "ref-lna.md",
            "Model Selection" => "ref-ms.md",
            "Gaussian Processes" => "ref-gp.md",
            "Kernels" => "ref-kernels.md"
        ]
    ]
)

gpabc_makedocs()

deploydocs(
    repo = "github.com/tanhevg/GpABC.jl.git"
)

ext: jl


text:
using XDiff
import ReverseDiffSource

function ann(w1, w2, w3, x1)
    x2 = w1 * x1
    # x2 = log.(1. + exp.(_x2))   # soft RELU unit
    x3 = w2 * x2
    # x3 = log.(1. + exp.(_x3))   # soft RELU unit
    x4 = sum.(w3 * x3)
    # return 1. ./ (1. + exp.(-x4))  # sigmoid output
    return x4
end

# function ann(w1, w2, w3, x1)
#     _x2 = w1 * x1
#     x2 = log.(1. + exp.(_x2))   # soft RELU unit
#     _x3 = w2 * x2
#     x3 = log.(1. + exp.(_x3))   # soft RELU unit
#     x4 = sum.(w3 * x3)
#     return 1. ./ (1. + exp.(-x4))  # sigmoid output
# end

# ANN input parameter types
types = (Matrix{Float64}, Matrix{Float64}, Matrix{Float64}, Vector{Float64})

# generate example input
w1, w2, w3 = randn(10,10), randn(10,10), randn(1,10)
x1 = randn(10)

dann = ReverseDiffSource.rdiff(ann, types)
_, rds_dw1_val, _, _, _ = dann(w1, w2, w3, x1)

dw1, _, _, _ = fdiff(ann, types)
dw1_val = dw1(w1, w2, w3, x1)

if isapprox(dw1_val, rds_dw1_val)
    println("Equal!")
else
    error("Not equal!")
end


text:
using BumpFuncs
using Test

@testset "BumpFuncs.jl" begin
    @test bump(0.0) ≈ 1.0 / ℯ
    @test bump([0.0; 0.0]) ≈ 1.0 / ℯ
    @test bump( 3.5, 3.5, 1.0 ) ≈ 1.0 / ℯ
    @test bump( [3.5; 3.5], [3.5, 3.5], [1.0,1.0] ) ≈ 1.0 / ℯ
    @test bump(1.0) == 0.0
    @test dbump(0.5) ≈ -0.46861713442795866
    @test ddbump(0.1) ≈ -0.758014676692494
    @test dbump(0.5,0.0,1.0) ≈ dbump(0.5)
    @test ddbump(0.5,0.0,1.0) ≈ ddbump(0.5)
    @test dbump([0.5],[0.2],[2.0])[1] ≈ dbump(0.5,0.2,2.0)
    @test ddbump([0.5],[0.2],[2.0])[1,1] ≈ ddbump(0.5,0.2,2.0)
    @test dbump([0.3])[1] ≈ dbump(0.3)
    @test ddbump([0.3])[1,1] ≈ ddbump(0.3)
end


text:
#=

A call to Resynthesizer using tiny 1D text

Usage:

>julia
julia> include("test/1D/tiny.jl")
julia> @time test(text, mask)
=#

# The tested item
include("../../src/apps/inpaint.jl")

# mock test data
# global so it persists and out of the test

#=
A String is not an AbstractArray.
But a byte-array string literal is, but is immutable.
Converts to a Vector{UInt8}, which is mutable.
=#
textLiteral = b"Now is the time for all good men to come to the aid of their country."
text = Vector{UInt8}(textLiteral)

# small mask
mask = falses(size(text))
mask[3:4] .= true

function test(text, mask)
    inpaint(text, mask)
end

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

"""
    struct BAT.LogUniform{T<:Real} <: Distributions.Distribution{Univariate,Continuous}

*Experimental feature, not part of stable public API.*

The log-uniform distribution (reciprocal distribution) over an interval ``[a, b]``.

Constructors:

* ```LogUniform(a::Real, b::Real)```

See also [Reciprocal distribution on Wikipedia](https://en.wikipedia.org/wiki/Reciprocal_distribution).
"""
struct LogUniform{T<:Real} <: Distributions.Distribution{Univariate,Continuous}
    a::T
    b::T
    _normf::T
end

function LogUniform{T}(a::Real, b::Real) where {T<:Real}
    @argcheck 0 < T(a) < T(b)
    LogUniform{T}(a, b, inv(log(b / a)))
end

function LogUniform(a::Real, b::Real)
    T = promote_type(typeof(a), typeof(b))
    LogUniform{T}(a, b)
end

LogUniform() = LogUniform{Float64}(1, 2)

Base.minimum(d::LogUniform{T}) where T = d.a
Base.maximum(d::LogUniform{T}) where T = d.b

StatsBase.params(d::LogUniform{T}) where T = (d.a, d.b)
@inline Distributions.partype(::LogUniform{T}) where T = T

Distributions.location(d::LogUniform) = d.a
Distributions.scale(d::LogUniform) = d.b - d.a

Base.eltype(::Type{LogUniform{T}}) where T = T

Statistics.mean(d::LogUniform{T}) where T = (d.b - d.a) * d._normf
StatsBase.median(d::LogUniform{T}) where T = quantile(d, 1//2)
StatsBase.mode(d::LogUniform{T}) where T = d.a
StatsBase.modes(d::LogUniform{T}) where T = Fill(mode(d), 0)

Statistics.var(d::LogUniform{T}) where T = (d.b^2 - d.a^2)/2 * d._normf - mean(d)^2
StatsBase.std(d::LogUniform{T}) where T = sqrt(var(d))
# StatsBase.skewness(d::LogUniform{T}) where T = ...
# StatsBase.kurtosis(d::LogUniform{T}) where T = ...

# StatsBase.entropy(d::LogUniform{T}) where T = ...

function Distributions.logpdf(d::LogUniform{T}, x::U) where {T,U<:Real}
    log(pdf(d, x))
end

function Distributions.pdf(d::LogUniform{T}, x::U) where {T,U<:Real}
    R = promote_type(T,U)
    y = inv(x) * d._normf
    ifelse(insupport(d, x), R(y), zero(R))
end

Distributions.logcdf(d::LogUniform, x::Real) = log(cdf(d, x))

function Distributions.cdf(d::LogUniform{T}, x::U) where {T,U<:Real}
    R = promote_type(T,U)
    p = log(x / d.a) * d._normf
    ifelse(insupport(d, x), R(p), zero(R))
end

Distributions.logccdf(d::LogUniform, x::Real) = log(ccdf(d, x))

function Distributions.ccdf(d::LogUniform{T}, x::U) where {T,U<:Real}
    y = cdf(d, x)
    one(y) - y
end

function Distributions.quantile(d::LogUniform{T}, p::U) where {T,U<:Real}
   R = promote_type(T,U)
   x = exp(p / d._normf) * d.a
   convert(float(R), x)
end

function Distributions.cquantile(d::LogUniform{T}, p::U) where {T,U<:Real}
    y = quantile(d, p)
    one(y) - y
end

# Distributions.mgf(d::LogUniform, t::Real) = ...
# Distributions.cf(d::LogUniform, t::Real) = ...

# # Distributions.gradlogpdf(d::LogUniform, x::Real) = ...

# Implemented implicitly via quantile:
# Base.rand(rng::AbstractRNG, d::LogUniform)

Distributions.truncated(d::LogUniform, l::Real, u::Real) = LogUniform(promote(max(l, d.a), min(u, d.b))...)
Distributions.truncated(d::LogUniform, l::Integer, u::Integer) = LogUniform(promote(max(l, d.a), min(u, d.b))...)


text:
export arithmetic_mean, geometric_mean, harmonic_mean, generalized_means
export pref_attach3, common3, jaccard3, adamic_adar3
export common_nbr_set

"""
arithmetic_mean
---------------

Returns the arithmetic mean of the weights of the edges of a list of triangles.

arithmetic_mean(triangles::Vector{NTuple{3,Int64}}, B::SpIntMat)

Input parameters:
- triangles::Vector{NTuple{3,Int64}}: The vector of triangles upon which to compute scores.
- B::SpIntMat: Projected graph as a Sparse integer matrix, where B[i, j] is the number of times that i and j co-appear in a simplex.
"""
function arithmetic_mean(triangles::Vector{NTuple{3,Int64}}, B::SpIntMat)
    scores = zeros(Float64, length(triangles))
    Threads.@threads for ind = 1:length(triangles)
        i, j, k = triangles[ind]
        scores[ind] = (B[i, j] + B[j, k] + B[i, k]) / 3
    end
    return scores
end

"""
geometric_mean
---------------

Returns the geometric mean of the weights of the edges of a list of triangles.

geometric_mean(triangles::Vector{NTuple{3,Int64}}, B::SpIntMat)

Input parameters:
- triangles::Vector{NTuple{3,Int64}}: The vector of triangles upon which to compute scores.
- B::SpIntMat: Projected graph as a Sparse integer matrix, where B[i, j] is the number of times that i and j co-appear in a simplex.
"""
function geometric_mean(triangles::Vector{NTuple{3,Int64}}, B::SpIntMat)
    scores = zeros(Float64, length(triangles))
    Threads.@threads for ind = 1:length(triangles)
        i, j, k = triangles[ind]
        scores[ind] = (B[i, j] * B[j, k] * B[i, k])^(1 / 3)
    end
    return scores
end

"""
harmonic_mean
---------------

Returns the harmonic mean of the weights of the edges of a list of triangles.

harmonic_mean(triangles::Vector{NTuple{3,Int64}}, B::SpIntMat)

Input parameters:
- triangles::Vector{NTuple{3,Int64}}: The vector of triangles upon which to compute scores.
- B::SpIntMat: Projected graph as a Sparse integer matrix, where B[i, j] is the number of times that i and j co-appear in a simplex.
"""
function harmonic_mean(triangles::Vector{NTuple{3,Int64}}, B::SpIntMat)
    scores = zeros(Float64, length(triangles))
    Threads.@threads for ind = 1:length(triangles)
        i, j, k = triangles[ind]
        scores[ind] = 3.0 / (1.0 / B[i, j] + 1.0 / B[j, k] + 1.0 / B[i, k])
    end
    return scores
end

"""
generalized_means
----------------

Computes the generalized p-means of the weights of the edges of a list of triangles.
The generalized mean of 3 values is

```math
M_p(x, y, z) = ((x^p + y^p + z^p) / 3)^{1/p}
```

generalized_means(triangles::Vector{NTuple{3,Int64}}, B::SpIntMat, ps::Float64=[-Inf; collect(-4:0.25:4); Inf])

Input parameters:
- triangles::Vector{NTuple{3,Int64}}: The vector of triangles upon which to compute scores.
- B::SpIntMat: Projected graph as a Sparse integer matrix, where B[i, j] is the number of times that i and j co-appear in a simplex.
- ps::Vector{Float64}=[-Inf; collect(-4:0.25:4); Inf]: the values of p for which to compute the means

Returns a matrix of size length(triangles) x length(ps) of the scores for the various generalized means.
"""
function generalized_means(triangles::Vector{NTuple{3,Int64}}, B::SpIntMat, ps::Vector{Float64}=[-Inf; collect(-4:0.25:4); Inf])
    scores = zeros(Float64, length(triangles), length(ps))
    Threads.@threads for ind = 1:length(triangles)
        i, j, k = triangles[ind]
        Bij, Bjk, Bik = B[i, j], B[j, k], B[i, k]
        for (r, p) in enumerate(ps)
            if     p == -Inf; scores[ind, r] = min(Bij, Bjk, Bik)
            elseif p == Inf;  scores[ind, r] = max(Bij, Bjk, Bik)
            elseif p == 0;    scores[ind, r] = (Bij * Bjk * Bik)^(1.0 / 3)
            else              scores[ind, r] = ((Bij^p + Bjk^p + Bik^p) / 3)^(1.0 / p)
            end
        end
    end
    return scores
end

"""
pref_attach3
------------

Returns the preferential attachment score for a degree vector

pref_attach3(triangles::Vector{NTuple{3,Int64}}, degrees::Vector{Int64})

The score of triangle (i, j, k) is degrees[i] * degrees[j] * degrees[k]

Input parameters:
- triangles::Vector{NTuple{3,Int64}}: The vector of triangles upon which to compute scores.
- degrees::Vector{Int64}: the degree vector
"""
function pref_attach3(triangles::Vector{NTuple{3,Int64}}, degrees::Vector{Int64})
    scores = zeros(Float64, length(triangles))
    Threads.@threads for ind = 1:length(triangles)
        i, j, k = triangles[ind]
        scores[ind] = 1.0 * degrees[i] * degrees[j] * degrees[k]
    end
    return scores
end

""" Return common neighbors of two nodes u and v. """
function common_nbr_set(common_nbrs::NbrSetMap, u::Int64, v::Int64)
    min_node, max_node = min(u, v), max(u, v)
    if !haskey(common_nbrs, (min_node, max_node))
        return Set{Int64}()
    end
    return common_nbrs[(min_node, max_node)]
end

"""
common3
------------

Returns the number of common 4th neighbors of a list of triangles.

common3(triangles::Vector{NTuple{3,Int64}}, common_nbrs::NbrSetMap)

Input parameters:
- triangles::Vector{NTuple{3,Int64}}: The vector of triangles upon which to compute scores.
- common_nbrs::NbrSetMap: the common neighbors map attained from common_neighbors_map()
"""
function common3(triangles::Vector{NTuple{3,Int64}}, common_nbrs::NbrSetMap)
    scores = zeros(Float64, length(triangles))
    Threads.@threads for ind = 1:length(triangles)
        i, j, k = triangles[ind]
        common_ij = common_nbr_set(common_nbrs, i, j)
        common_ik = common_nbr_set(common_nbrs, i, k)
        common_jk = common_nbr_set(common_nbrs, j, k)
        common_ijk = intersect(common_ij, common_ik, common_jk)
        scores[ind] = length(common_ijk)
    end
    return scores
end

"""
jaccard3
------------

Returns the 3-way Jaccard index of the neighbor lists of nodes in triangle.

jaccard3(triangles::Vector{NTuple{3,Int64}}, common_nbrs::NbrSetMap, degrees::Vector{Int64})

Input parameters:
- triangles::Vector{NTuple{3,Int64}}: The vector of triangles upon which to compute scores.
- common_nbrs::NbrSetMap: the common neighbors map attained from common_neighbors_map()
- degrees::Vector{Int64}: degree of each node in the projected graph
"""
function jaccard3(triangles::Vector{NTuple{3,Int64}}, common_nbrs::NbrSetMap,
                  degrees::Vector{Int64})
    scores = zeros(Float64, length(triangles))
    Threads.@threads for ind = 1:length(triangles)
        i, j, k = triangles[ind]
        common_ij = common_nbr_set(common_nbrs, i, j)
        common_ik = common_nbr_set(common_nbrs, i, k)
        common_jk = common_nbr_set(common_nbrs, j, k)
        common_ijk = intersect(common_ij, common_ik, common_jk)
        nij, nik, njk = length(common_ij), length(common_ik), length(common_jk)
        nijk = length(common_ijk)
        di, dj, dk = degrees[[i, j, k]]
        scores[ind] = nijk / (di + dj + dk - nij - nik - njk + nijk)
    end
    return scores
end

"""
adamic_adar3
------------

Returns the 3-way Adamic-Adar score for a list of triangles.

adamic_adar3(triangles::Vector{NTuple{3,Int64}}, common_nbrs::NbrSetMap, degrees::Vector{Int64})

Input parameters:
- triangles::Vector{NTuple{3,Int64}}: The vector of triangles upon which to compute scores.
- common_nbrs::NbrSetMap: the common neighbors map attained from common_neighbors_map()
- degrees::Vector{Int64}: degree of each node in the projected graph
"""
function adamic_adar3(triangles::Vector{NTuple{3,Int64}}, common_nbrs::NbrSetMap,
                      degrees::Vector{Int64})
    scores = zeros(Float64, length(triangles))
    Threads.@threads for ind = 1:length(triangles)
        i, j, k = triangles[ind]
        common_ij = common_nbr_set(common_nbrs, i, j)
        common_ik = common_nbr_set(common_nbrs, i, k)
        common_jk = common_nbr_set(common_nbrs, j, k)
        common_ijk = intersect(common_ij, common_ik, common_jk)
        scores[ind] = sum([1.0 / log(degrees[z]) for z in common_ijk])
    end
    return scores
end


text:
using LinearAlgebra
using Statistics
"""
    InterpretMe.partial_dependence(data::Array{<:Any,2},
                                   features::Union{<:Integer,Vector{<:Integer}},
                                   estimate_fn, grid_size::Int,
                                   return_ice::Bool)
Calculate the partial dependance of a machine learning model.

`estimate_fn` is a function from the model that generates a vector of predictions from a matrix of inputs.
`feature_names` should be the feature names corresponding to the features in `data`
`return_ice` whether to return the ICE lines.
Returns
`pd` the partial dependancy results
`grid_values` the grid values for all the selected features
"""
function partial_dependence(data::Array{<:Any,2},
                            features::Union{<:Integer,Vector{<:Integer}},
                            estimate_fn, grid_size::Int, return_ice::Bool)
    # Mainly inspired by Sk-learn's function.
    # https://scikit-learn.org/stable/modules/generated/sklearn.inspection.partial_dependence.html
    #Step 1 Necessary checking
    if grid_size <=1
         raise(error("grid_size should be greater than 1, $grid_size was passed"))
    end
    #Step 2 Get_grid values
    grid, grid_values = get_grid(data[:,features], grid_size)
    #Step 3 evaluate partial_dependance
    probe = estimate_fn(data[1:2,:])
    pd = typeof(probe)(undef, size(grid,1), size(probe,2))
    # pd = []
    for grid_value_index in 1:size(grid,1)
        temp_data = copy(data)
        # (n,length(features)) .= (1,length(features))
        temp_data[:,features] .= grid[grid_value_index,:]
        predictions = estimate_fn(temp_data)
        pd[grid_value_index,:] = mean(predictions, dims=1)
    end
    pd, grid_values
end

function get_grid(data_features, grid_size)
    values = Nothing
    for feature in 1:size(data_features,2)
        unique_vals = unique(data_features[:,feature])
        value =
            if length(unique_vals) <= grid_size
                sort(unique_vals)
            else
                LinRange(minimum(data_features[:,feature]),
                 maximum(data_features[:,feature]), grid_size)
            end
        if values == Nothing
            values = [value]
        else
            push!(values, value)
        end
    end
    cartesian_prod_looped(values), values
end

function cartesian_prod(values)
    if length(values) == 1
        values[1]
    else
        cartesian_prod([
        hcat(repeat(values[1],inner=(size(values[2],1),
                                    [1 for i in 1:ndims(values[1])-1]...))
            ,repeat(values[2],outer=(size(values[1],1),
                                    [1 for i in 1:ndims(values[2])-1]...))
            )
         , values[3:end]...])
     end
 end

 function cartesian_prod_looped(values)
     while length(values) > 1
         values = [
         hcat(repeat(values[1],inner=(size(values[2],1),
                                     [1 for i in 1:ndims(values[1])-1]...))
             ,repeat(values[2],outer=(size(values[1],1),
                                     [1 for i in 1:ndims(values[2])-1]...))
             )
          , values[3:end]...]
      end
      values[1]
      # To return transposed.
      # res_values = typeof(values[1])(undef,size(transpose(values[1]))...)
      # using LinearAlgebra
      # transpose!(res_values, values[1])
      # or
      # typeof(values[1])(transpose(values[1]))
  end


text:
module DirichletFESpacesTests

using Test
using Gridap.Algebra
using Gridap.Geometry
using Gridap.FESpaces
using Gridap.CellData
using Gridap.ReferenceFEs

n = 2
mesh = (n,n)
domain = (0,1,0,1)
order = 1
dirichlet_tags = [1,3,4,6,7]
model = CartesianDiscreteModel(domain, mesh)
labels = get_face_labeling(model)
add_tag_from_tags!(labels,"dirichlet",dirichlet_tags)

reffe = ReferenceFE(lagrangian,Float64,order)
Vf = TestFESpace( model,reffe; conformity=:H1, labels=labels, dirichlet_tags="dirichlet")

Vd = DirichletFESpace(Vf)
matvecdata = ([],[],[])
matdata = ([],[],[])
vecdata = ([],[])
test_single_field_fe_space(Vd,matvecdata,matdata,vecdata)

Uf = TrialFESpace(Vf,1)

Ud = TrialFESpace(Vd)
test_single_field_fe_space(Ud,matvecdata,matdata,vecdata)

trian = Triangulation(model)
degree = 2*order
dΩ = Measure(trian,degree)

a(u,v) = ∫( u*v )*dΩ
l(v) = ∫( v*0 )*dΩ

op_ff = AffineFEOperator(a,l,Uf,Vf)
op_fd = AffineFEOperator(a,l,Ud,Vf)
op_df = AffineFEOperator(a,l,Uf,Vd)

xd = get_dirichlet_values(Uf)

@test maximum(abs.(get_vector(op_ff) + get_matrix(op_fd)*xd)) < 1.0e-10

end # module


text:
# Вычисление начального приближения с использованием асимптотических методов

function front_velocity(f1_data, Tₘ, M)

    df1dt = zero(f1_data);
    τ = Tₘ[2] - Tₘ[1];

    df1dt[1] = (-3 * f1_data[1] + 2 * f1_data[2] - f1_data[3]) / (2 * τ);
    for m in 2:M-2
        df1dt[m] = (f1_data[m+1] - f1_data[m-1]) / (2 * τ);
    end
    df1dt[M+1] = (f1_data[M+1 - 2] - 2 * f1_data[M+1 - 1] + 3 * f1_data[M+1]) / (2 * τ);

    return df1dt;
end

# Найдем минимальный и максимальный номер узлов сетки Xₙ, которые пересекает
# переходный слой
# .     Узлы сетки
# x     f1(t) положение переходного слоя
#                     k                                   K
# .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .
#                   x x x x x x x x x x x x x x x x x x x x x
function get_nodes_range_where_fron_passed(x, Xₙ)
    k = findfirst(z -> z > minimum(x), Xₙ)
    K = findfirst(z -> z > maximum(x), Xₙ) - 1

    # Первые и последние точки, часто оказываются плохими
    # Отрежем на всякий случай
    return k, K
end

function get_A_matrix(k, K, Xₙ, N)
    h = Xₙ[2] - Xₙ[1];

    Al = LowerTriangular(h * ones(K-k+1, K-k+1));
    Ar = UpperTriangular(-h * ones(K-k+1, K-k+1));
    Al[:, 1]    .= (k-1)*h;
    Ar[:, end]  .= (N+1-K)*h;
    # Повторим чертеж для объяснения двух прошлых строк
    # N = 17, тогда нумерация массива будет следующей
    # k = 6, K = 16
    # Интервалов лежащих слева от k -- 5, (k-1)
    # Интервалов лежащих справа от K -- 2, (N+1 - K)
    # 1   2   3   4   5   6   7   8   9   11  12  13  14  15  16  17  18
    #                     k                                   K
    # .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .
    #                   x x x x x x x x x x x x x x x x x x x x x

    return A = -(Al + Ar) / 2.0;
end

function get_B_vector(k, K, Vn, ulₘ, urₘ)
    # `B` составим из ненулевых v_n
    B = Vn[k:K];
    B = B .- ( -(ulₘ[1] + ulₘ[1])  / 2.0);

    return B;
end

function solve_system_for_initial_guess(k, K, N, A, B, α = 0.001)

    qx = ( transpose(A) * A + α*I) \ transpose(A)*B

    l = [qx[2] for i in 1:k-1];     # Дополняем начальное приближение константами
    r = [qx[end-1] for i in K+1:N+1]; # Дополняем Справа
    qx[end] = qx[end-1];
    qx[1] = qx[2];

    q_guess = [l; qx; r];
end

function initial_guess(f1_data, Xₙ, N, Tₘ, M, ulₘ, urₘ, α = 0.001)

    df1dt = front_velocity(f1_data, Tₘ, M);

    # Последние точки некрасивые, обрежем их
    x = f1_data[2:end-3];
    y = df1dt[2:end-3];
    tm = Tₘ[2:end-3];

    k, K = get_nodes_range_where_fron_passed(x, Xₙ);
    @info "Фронт прошел точки $((k, K)), $(div((100*(K-k)),(N+1)))% кадра"

    # Найдем скорость фронта как функцию координат
    # Иначе говоря, найдем скорость, которой обладал фронт проходя узел сетки
    vspl = Spline1D(x, y);
    Vn = zero(Xₙ);
    for n in k:K
            Vn[n] = vspl(Xₙ[n])
    end

    A = get_A_matrix(k, K, Xₙ, N);
    B = get_B_vector(k, K, Vn, ulₘ, urₘ);
    q_guess = solve_system_for_initial_guess(k, K, N, A, B, α);

    return q_guess;
end


text:
"""
    emptybatch(model::T) where {T<:Model}

Create an empty batched input array on the GPU if available.
"""
function emptybatch(model::T) where {T<:Model}
    modelInputSize = getModelInputSize(model)
    gpu(Array{Float32}(undef, modelInputSize...))
end

"""
    flipdict(dict::Dict)

Flip the key=>value pair for each entry in a dict.
"""
flipdict(dict::Dict) = Dict(map(x->(dict[x],x),collect(keys(dict))))

"""
    createcountdict(dict::Dict)

Create a dict copy of namesdict, for counting the occurances of each named object.
"""
createcountdict(dict::Dict) = Dict(map(x->(x,0),collect(keys(dict))))

"""
    drawBoxes(img::Array, results)

Draw boxes on image for each BBOX result.
"""
function drawBoxes(img::Array, results)
    w = size(img,1)
    h = size(img,2)
    imgCopy = copy(img)
    size(results) == (0,) && return imgCopy
    for i in 1:size(results,2)
        bbox = results[1:4,i]
        class = results[end-1,i]
        conf = results[end-2,i]
        p = Point(round(Int, bbox[1]*w), round(Int, bbox[2]*h))
        q = Point(round(Int, bbox[1]*w), round(Int, bbox[4]*h))
        r = Point(round(Int, bbox[3]*w), round(Int, bbox[2]*h))
        s = Point(round(Int, bbox[3]*w), round(Int, bbox[4]*h))
        draw!(imgCopy, Polygon([p,q,s,r]), zero(eltype(imgCopy)))
    end
    return imgCopy
end

"""
    benchmark(;select = [1,3,4,5,6], reverseAfter:Bool=false)

Convenient benchmarking
"""
function benchmark(;select = [1,3,4,5,6], reverseAfter::Bool = false)
    pretrained_list = [
                        YOLO.v2_tiny_416_COCO,
                        YOLO.v2_608_COCO,
                        YOLO.v3_tiny_416_COCO,
                        YOLO.v3_320_COCO,
                        YOLO.v3_416_COCO,
                        YOLO.v3_608_COCO,
                        YOLO.v3_608_spp_COCO
                        ][select]
    reverseAfter && (pretrained_list = vcat(pretrained_list, reverse(pretrained_list)))
    IMG = rand(RGB,416,416)

    header = ["Model" "loaded?" "load time (s)" "ran?" "run time (s)" "run time (fps)"]
    table = Array{Any}(undef, length(pretrained_list), 6)
    for (i, pretrained) in pairs(pretrained_list)
        modelname = string(pretrained)
        @info "Loading and running $modelname"
        table[i,:] = [modelname false "-" "-" "-" "-"]

        t_load = @elapsed begin
            mod = pretrained(silent=true)
        end
        table[i, 2] = true
        table[i, 3] = round(t_load, digits=3)

        batch = emptybatch(mod)
        batch[:,:,:,1] .= gpu(resizePadImage(IMG, mod))

        res = mod(batch) #run once
        t_run = @belapsed $mod($batch);
        table[i, 4] = true
        table[i, 5] = round(t_run, digits=4)
        table[i, 6] = round(1/t_run, digits=1)

        mod = nothing
        batch = nothing
        GC.gc()
    end
    pretty_table(table, header)
end


text:
"""
    PolynomialCalibrationModel{T,F}(k,β,c::Calibration{T,F})

A CalibrationModel representing a polynomial regression.
"""
struct PolynomialCalibrationModel{T<:Quantity, F<:Quantity} <: CalibrationModel{T,F}
    k::Int
    β::Vector{Float64}
    c::Calibration{T,F}
end

StatsBase.coef(m::PolynomialCalibrationModel) = m.β

"""
    design_polynomial(Q,k)

Create the design matrix for a polynomial regression
of highest order k
"""
function design_polynomial(Q,k)
    X = ones(length(Q),k+1)
    for i in 1:k+1
        X[:,i] = Q.^(i-1)
    end
    X
end

function design_polynomial(Q::Quantity,k)
    design_polynomial(quantity(Q),k)
end

function design_polynomial(c::Calibration,k)
    qi = interpolatecal(c)
    design_polynomial(qi,k)
end

function StatsBase.fit(::Type{PolynomialCalibrationModel},c::Calibration,k)
    X = design_polynomial(c,k)
    y = quantity(to_quantity(c))
    β = X\y
    PolynomialCalibrationModel(k,β,c)
end

function StatsBase.predict{T,F}(m::PolynomialCalibrationModel{T,F},q::F)
    X = design_polynomial(q,m.k)
    T(times(q),X*m.β)
end


text:
using Test
using TropicalTensors, TropicalTensors.Reversible

@testset "reversible square" begin
    include("square.jl")
end

@testset "reversible chimera" begin
    include("chimera.jl")
end

@testset "reversible second_neighbor" begin
    include("second_neighbor.jl")
end


text:
julia> sizeof(Int8)
1

julia> t = 1
1

julia> sizeof(t)
8

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
Base.broadcastable(layer::SimpleSDMResponse) = layer
Base.BroadcastStyle(::Type{<:SimpleSDMResponse}) = Broadcast.Style{SimpleSDMResponse}()

Base.broadcastable(layer::SimpleSDMPredictor) = similar(layer)
Base.BroadcastStyle(::Type{<:SimpleSDMPredictor}) = Broadcast.Style{SimpleSDMPredictor}()

function Base.Broadcast.broadcast(f, L::LT) where {LT <: SimpleSDMLayer}
    newgrid = Array{Any}(nothing, size(L))
    N = SimpleSDMResponse(newgrid, L)
    v = filter(!isnothing, L.grid)
    fv = f.(v)
    N.grid[findall(!isnothing, L.grid)] .= fv

    internal_types = unique(typeof.(N.grid))

    RT = LT <: SimpleSDMResponse ? SimpleSDMResponse : SimpleSDMPredictor
    return RT(convert(Matrix{Union{internal_types...}}, N.grid), N)
end


text:
"""
This is a file where code to be used for your module can appear. Files like
these are *not* modules themselves, so don't wrap them in one, or your
`include()`s won't run!

The example function, `sharpe()`, also demonstrates how to use docstrings in
Julia. However, refer to the official docs to read about best practices, like
not documenting kwargs if they're self-explanatory.
"""

using Statistics

"""
    sharpe(x[, rf])

Calculate the Sharpe Ratio of a financial returns series `x`, using `rf` as the
risk-free rate. If `rf` is missing, use 0.02.

# Arguments
- `x::Array{Float64}`: An array of returns data
- `rf::Float64`: The risk-free rate

# Examples
```jldoctest
using PkgName

returns = [0.03, 0.04, 0.01]
sharpe(returns, rf = 0.02)

# output

0.4364357804719848
```
"""
function sharpe(x::Array{Float64}; rf::Float64 = 0.02)::Float64
    return ((mean(x) - rf) / std(x))
end

max_stars_repo_name: ryapric/boilerplate

max_issues_repo_name: ryapric/boilerplate

max_forks_repo_name: ryapric/boilerplate


text:
const TSYSTEM_TABLE = Dict{UInt,Vector{System}}()

"""
    CompiledSystem <: AbstractSystem

An [`AbstractSystem`](@ref) which automatically compiles a straight line program for the
fast evaluation of a system `F` and its Jacobian.
For large systems the compilation can take some time and require a large amount of memory.
If this is a problem consider [`InterpretedSystem`](@ref).

    CompiledSystem(F::System; optimizations = true)

Construct a `CompiledSystem` from the given [`System`](@ref) `F`. If `optimizations = true`
then [`optimize`](@ref) is called on `F` before compiling.
"""
struct CompiledSystem{HI} <: AbstractSystem
    nexpressions::Int
    nvariables::Int
    nparameters::Int
    system::System
end

function CompiledSystem(F::System; optimizations::Bool = true)
    F = optimizations ? optimize(F) : F
    n = length(F)
    nvars = nvariables(F)
    nparams = nparameters(F)

    # We substitute all variable and parameter names away such that two systems
    # compile to the same code independently of the names
    @var _x_[1:nvars] _p_[1:nparams]
    D = ExpressionMap()
    for (x, y) in zip(variables(F), _x_)
        D[x] = y
    end
    for (x, y) in zip(parameters(F), _p_)
        D[x] = y
    end
    cleard_exprs = subs(expressions(F), D)
    sys = System(cleard_exprs, _x_, _p_)
    h = hash(cleard_exprs)
    k = 0
    if haskey(TSYSTEM_TABLE, h)
        # check that it is identical
        for (i, vi) in enumerate(TSYSTEM_TABLE[h])
            if vi == sys
                k = i
                break
            end
        end
        if k == 0
            push!(TSYSTEM_TABLE[h], sys)
            k = length(TSYSTEM_TABLE[h])
        end
    else
        k = 1
        TSYSTEM_TABLE[h] = [sys]
    end
    return CompiledSystem{(h, k)}(n, nvars, nparams, F)
end
function Base.show(io::IO, TS::CompiledSystem)
    print(io, "Compiled: ")
    show(io, TS.system)
end

(F::CompiledSystem)(x, p = nothing) = F.system(x, p)

interpret(TS::CompiledSystem) = TS.system
interpret(::Type{CompiledSystem{HI}}) where {HI} = TSYSTEM_TABLE[first(HI)][last(HI)]

Base.size(CS::CompiledSystem) = (CS.nexpressions, CS.nvariables)
Base.size(CS::CompiledSystem, i::Integer) = size(CS)[i]
Base.length(CS::CompiledSystem) = CS.nexpressions
nparameters(CS::CompiledSystem) = CS.nparameters
nvariables(CS::CompiledSystem) = CS.nvariables

variables(F::CompiledSystem) = variables(F.system)
parameters(F::CompiledSystem) = parameters(F.system)
variable_groups(F::CompiledSystem) = variable_groups(F.system)
System(F::CompiledSystem) = F.system

Base.:(==)(F::CompiledSystem{A}, G::CompiledSystem{B}) where {A,B} = A == B

######################
## CompiledHomotopy ##
######################

const THOMOTOPY_TABLE = Dict{UInt,Vector{Homotopy}}()

"""
    CompiledHomotopy <: AbstractHomotopy

An [`AbstractHomotopy`](@ref) which automatically compiles a straight line program for the
fast evaluation of a homotopy `H` and its Jacobian.
For large homotopies the compilation can take some time and require a large amount of
memory. If this is a problem consider [`InterpretedHomotopy`](@ref).

    CompiledSystem(F::System; optimizations = true)

Construct a `CompiledSystem` from the given [`System`](@ref) `F`. If `optimizations = true`
then [`optimize`](@ref) is called on `F` before compiling.
"""
struct CompiledHomotopy{HI} <: AbstractHomotopy
    nexpressions::Int
    nvariables::Int
    nparameters::Int
    homotopy::Homotopy
end

function CompiledHomotopy(H::Homotopy)
    n = length(H)
    nvars = nvariables(H)
    nparams = nparameters(H)

    # We substitute all variable and parameter names away such that two homotopies
    # compile to the same code independently of the names
    @var _x_[1:nvars] _t_ _p_[1:nparams]
    D = ExpressionMap()
    for (x, y) in zip(variables(H), _x_)
        D[x] = y
    end
    for (x, y) in zip(parameters(H), _p_)
        D[x] = y
    end
    D[H.t] = _t_
    cleard_exprs = subs(expressions(H), D)
    homotopy = Homotopy(cleard_exprs, _x_, _t_, _p_)
    h = hash(cleard_exprs)

    k = 0
    if haskey(THOMOTOPY_TABLE, h)
        # check that it is identical
        for (i, vi) in enumerate(THOMOTOPY_TABLE[h])
            if vi == homotopy
                k = 1
                break
            end
        end
        if k == 0
            push!(THOMOTOPY_TABLE[h], homotopy)
            k = 1
        end
    else
        k = 1
        THOMOTOPY_TABLE[h] = [homotopy]
    end
    return CompiledHomotopy{(h, 1)}(n, nvars, nparams, H)
end

Base.size(CH::CompiledHomotopy) = (CH.nexpressions, CH.nvariables)
Base.size(CH::CompiledHomotopy, i::Integer) = size(CH)[i]
Base.length(CH::CompiledHomotopy) = CH.nexpressions
nparameters(CH::CompiledHomotopy) = CH.nparameters
nvariables(CH::CompiledHomotopy) = CH.nvariables

function Base.show(io::IO, TH::CompiledHomotopy)
    print(io, "Compiled: ")
    show(io, TH.homotopy)
end

interpret(CH::CompiledHomotopy) = CH.homotopy
interpret(::Type{CompiledHomotopy{HI}}) where {HI} = THOMOTOPY_TABLE[first(HI)][last(HI)]

################
## EVALUATION ##
################
function _evaluate!_impl(T)
    I = interpret(T)
    list, out = instruction_list(expressions(I))
    generate_evaluate_impl(
        list,
        out;
        variables = Symbol.(variables(I)),
        parameters = Symbol.(parameters(I)),
        t = I isa Homotopy ? Symbol(I.t) : nothing,
    )
end

"""
    evaluate!(u, T::CompiledSystem, x, p = nothing)

Evaluate `T` for variables `x` and parameters `p` and store result in `u`.
"""
@generated function evaluate!(u, T::CompiledSystem, x, p = nothing)
    _evaluate!_impl(T)
end

"""
    evaluate!(u, T::CompiledHomotopy, x, t, p = nothing)

Evaluate `T` for variables `x`, `t` and parameters `p` and store result in `u`.
"""
@generated function evaluate!(u, T::CompiledHomotopy, x, t, p = nothing)
    _evaluate!_impl(T)
end

function _evaluate_and_jacobian!_impl(T)
    I = interpret(T)
    list, out = instruction_list(expressions(I))
    generate_evaluate_and_jacobian_impl(
        list,
        out;
        variables = Symbol.(variables(I)),
        parameters = Symbol.(parameters(I)),
        t = I isa Homotopy ? Symbol(I.t) : nothing,
    )
end

"""
    evaluate_and_jacobian!(u, U, T::CompiledSystem, x, p = nothing)

Evaluate `T` and its Jacobian for variables `x` and parameters `p` and
store result in `u`.
"""
@generated function evaluate_and_jacobian!(u, U, T::CompiledSystem, x, p = nothing)
    _evaluate_and_jacobian!_impl(T)
end
#
"""
    evaluate_and_jacobian!(u, U, T::CompiledHomotopy, x, t, p = nothing)

Evaluate `T` and its Jacobian for variables `x`, `t` and parameters `p` and
store result in `u`.
"""
@generated function evaluate_and_jacobian!(u, U, T::CompiledHomotopy, x, t, p = nothing)
    _evaluate_and_jacobian!_impl(T)
end

#
"""
    jacobian!(U, T::CompiledSystem, x, p = nothing)

Evaluate the Jacobian of `T` for variables `x`, `t` and parameters `p`
and store result in `u`.
"""
function jacobian!(U, T::CompiledSystem, x, p = nothing)
    evaluate_and_jacobian!(nothing, U, T, x, p)
end

"""
    jacobian!(U, T::CompiledHomotopy, x, t, p = nothing)

Evaluate the Jacobian of `T` for variables `x`, `t` and parameters `p` and
store result in `u`.
"""
function jacobian!(U, T::CompiledHomotopy, x, t, p = nothing)
    evaluate_and_jacobian!(nothing, U, T, x, t, p)
end

function _taylor!_impl(T)
    I = interpret(T)
    list, out = instruction_list(expressions(I))
    generate_taylor_impl(
        list,
        out;
        variables = Symbol.(variables(I)),
        parameters = Symbol.(parameters(I)),
        t = I isa Homotopy ? Symbol(I.t) : nothing,
    )
end

"""
    taylor!(
        u,
        ::Val{K},
        F::CompiledSystem,
        x,
        p = nothing,
    )

Compute the Taylor series order order `K` of ``u = F(x,p)``.
"""
@generated function taylor!(
    u::AbstractArray,
    Order::Val{K},
    T::CompiledSystem,
    x::AbstractArray,
    p::Union{Nothing,AbstractArray} = nothing,
) where {K}
    _taylor!_impl(T)
end

"""
    taylor!(
        u,
        ::Val{K},
        H::CompiledHomotopy,
        x,
        p = nothing,
    )

Compute the Taylor series order order `K` of ``u = H(x,t,p)``.
"""
@generated function taylor!(
    u::AbstractArray,
    Order::Val{K},
    T::CompiledHomotopy,
    x::AbstractArray,
    t,
    p::Union{Nothing,AbstractArray} = nothing,
) where {K}
    _taylor!_impl(T)
end


text:
module SemanticNames

using Static
using Static: static, True, False, StaticSymbol

export @defname

struct IsName{name} <: Function
    IsName{name}() where {name} = new{name::Symbol}()
end

(isn::IsName{name})(x::Symbol) where {name} = isn(static(x))
(isn::IsName{name})(x::Tuple{Vararg{Symbol}}) where {name} = isn(static(x))
(isn::IsName{name})(x::Tuple{Vararg{StaticSymbol}}) where {name} = isn(typeof(x))
(isn::IsName{name})(x::StaticSymbol) where {name} = isn(typeof(x))
(isn::IsName{name})(::Type{StaticSymbol{sym}}) where {name,sym} = False()
(isn::IsName{name})(::Type{StaticSymbol{name}}) where {name} = True()
(isn::IsName{name})(::Type{IsName{sym}}) where {name,sym} = isn(static(sym))

macro name_str(s) :(SemanticName{$(Expr(:quote, Symbol(s)))}) end

macro defname(name)
    is_name = esc(Symbol(:is_, name))
    is_not_name = esc(Symbol(:is_not_, name))

    quote
        const $is_name = IsName{$(QuoteNode(name))}()

        macro $is_name(name)
            new_is_name($(QuoteNode(name)), name)
        end

        macro $is_name(name, ns...)
            new_is_name($(QuoteNode(name)), name, ns...)
        end

        macro $is_not_name(name)
            new_is_not_name($(QuoteNode(name)), name)
        end

        macro $is_not_name(name, ns...)
            new_is_not_name($(QuoteNode(name)), name, ns...)
        end

        nothing
    end
end

function new_is_name(is_name, name)
    @eval begin
        (::IsName{$(QuoteNode(is_name))})(::Type{StaticSymbol{$name}}) = True()
        nothing
    end
end

function new_is_name(is_name, name, ns...)
    T = Expr(:curly, :Tuple)
    push!(T.args, Expr(:curly, :StaticSymbol, name))
    for n in ns
        push!(T.args, Expr(:curly, :StaticSymbol, n))
    end
    @eval begin
        (::IsName{$(QuoteNode(is_name))})(::Type{$T}) = True()
        nothing
    end
end

function new_is_not_name(is_name, name)
    @eval begin
        (::IsName{$(QuoteNode(is_name))})(::Type{StaticSymbol{$name}}) = False()
        nothing
    end
end

function new_is_note_name(is_name, name, ns...)
    T = Expr(:curly, :Tuple)
    push!(T.args, Expr(:curly, :StaticSymbol, name))
    for n in ns
        push!(T.args, Expr(:curly, :StaticSymbol, n))
    end
    @eval begin
        (::IsName{$(QuoteNode(is_name))})(::Type{$T}) = False()
        nothing
    end
end

function Base.show(io::IO, ::MIME"text/plain", isn::IsName{name}) where {name}
    fxnname = Symbol(:is_, name)
    nmethods = length(methods(isn))
    if nmethods == 1
        print(io, "$fxnname (generic function with 1 method)")
    else
        print(io, "$fxnname (generic function with $nmethods methods)")
    end
end

end


text:
#=

  Euler #17 in Julia.

  Problem 17:
  """
  If the numbers 1 to 5 are written out in words: one, two, three, four, five,
  then there are 3 + 3 + 5 + 4 + 4 = 19 letters used in total.

  If all the numbers from 1 to 1000 (one thousand) inclusive were written out in
  words, how many letters would be used?

  NOTE: Do not count spaces or hyphens. For example, 342 (three hundred and forty-two)
  contains 23 letters and 115 (one hundred and fifteen) contains 20 letters. The use of
  "and" when writing out numbers is in compliance with British usage.
  """

  This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

include("Euler.jl")

function english(n)
    divs  =      [1_000_000_000, 1_000_000,  1_000,       100]
    divnames  =  ["billion", "million", "thousand", "hundred"]
    prefixes  =  ["0", "twen", "thir", "for", "fif", "six", "seven", "eigh", "nine"]
    _ordinals  = ["first", "second", "third", "fourth", "fifth", "sixth", "seventh",
                      "eighth", "ninth", "tenth", "eleventh", "twelfth", "thirteenth",
                      "fourteenth","fifteenth", "sixteenth", "seventeenth",
                      "eighteenth", "nineteenth"]
    cardinals =  ["one", "two", "three", "four", "five", "six", "seven",
                      "eight", "nine", "ten", "eleven", "twelve", "thirteen", "fourteen",
                      "fifteen", "sixteen", "seventeen", "eighteen", "nineteen"];

    s = ""
    printed = 0
    if n < 0
        s = "minus" * s;
        n = -n
    end
    d = 0
    for i in 1:length(divs)
        d = floor(Int,n / divs[i])
        n %= divs[i]
        if d != 0
            s *= english(d) * divnames[i];
            printed = 1
        end
    end

    if n > 0 && printed == 1
        s *= "and";
    end
    if n == 0
        # dummy
    elseif n > 19
        d = floor(Int,n / 10)
        n %= 10
        s *= prefixes[d] * "ty" * english(n)
    else
        s *= cardinals[n]
    end

    return s
end

# 0.00138457s
function euler17a()
    return sum((1:1000).|>i->length(english(i)))
end

run_euler(euler17a)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
module AbstractOperator
# basic operator types to do inversion with
abstract type Operator end
abstract type Operator1D <: Operator end
abstract type Operator2D <: Operator end

# any physics function must be a subtype of the
# above operators and needs to provide a misfit
function get_misfit end

# for gradient based inversion
function getresidual end

# many geophysics data types require a sounding type
abstract type Sounding end
# these soundings often need to be fed into a function
# to make a physics operator with the relevant sounding data
function makeoperator end
# for nuisances this may change from sounding to sounding
function make_tdgp_opt end

export Operator, Operator1D, Operator2D
end

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
function rsu7(in_ch=3, mid_ch=12, out_ch=3)
    rebnconvin = Chain(rebnconv(in_ch,out_ch,1)...)

    rebnconv1 = rebnconv_pool(out_ch, mid_ch, 1)
    rebnconv2 = rebnconv_pool(mid_ch, mid_ch, 1)
    rebnconv3 = rebnconv_pool(mid_ch, mid_ch, 1)
    rebnconv4 = rebnconv_pool(mid_ch, mid_ch, 1)
    rebnconv5 = rebnconv_pool(mid_ch, mid_ch, 1)

    rebnconv6 = Chain(rebnconv(mid_ch,mid_ch,1)...)
    rebnconv7 = Chain(rebnconv(mid_ch,mid_ch,2)...)

    rebnconv6d = rebnconv_upsample(mid_ch*2,mid_ch,1)
    rebnconv5d = rebnconv_upsample(mid_ch*2,mid_ch,1)
    rebnconv4d = rebnconv_upsample(mid_ch*2,mid_ch,1)
    rebnconv3d = rebnconv_upsample(mid_ch*2,mid_ch,1)
    rebnconv2d = rebnconv_upsample(mid_ch*2,mid_ch,1)

    rebnconv1d = Chain(rebnconv(mid_ch*2,out_ch,1)...)

    layer = x -> begin
        hxin = rebnconvin(x)

        hx1, hx = rebnconv1(hxin)
        hx2, hx = rebnconv2(hx)
        hx3, hx = rebnconv3(hx)
        hx4, hx = rebnconv4(hx)
        hx5, hx = rebnconv5(hx)

        hx6 = rebnconv6(hx)
        hx7 = rebnconv7(hx6)

        hx6dup = rebnconv6d(hx5, hx6, hx7)
        hx5dup = rebnconv5d(hx4, hx5, hx6dup)
        hx4dup = rebnconv4d(hx3, hx4, hx5dup)
        hx3dup = rebnconv3d(hx2, hx3, hx4dup)
        hx2dup = rebnconv6d(hx1, hx2, hx3dup)

        hx1d = rebnconv1d(cat(hx2dup, hx1; dims=3))

        return hx1d + hxin
    end
end

function rsu6(in_ch=3, mid_ch=12, out_ch=3)
    rebnconvin = Chain(rebnconv(in_ch,out_ch,1)...)

    rebnconv1 = rebnconv_pool(out_ch, mid_ch, 1)
    rebnconv2 = rebnconv_pool(mid_ch, mid_ch, 1)
    rebnconv3 = rebnconv_pool(mid_ch, mid_ch, 1)
    rebnconv4 = rebnconv_pool(mid_ch, mid_ch, 1)

    rebnconv5 = Chain(rebnconv(mid_ch,mid_ch,1)...)
    rebnconv6 = Chain(rebnconv(mid_ch,mid_ch,2)...)

    rebnconv5d = rebnconv_upsample(mid_ch*2,mid_ch,1)
    rebnconv4d = rebnconv_upsample(mid_ch*2,mid_ch,1)
    rebnconv3d = rebnconv_upsample(mid_ch*2,mid_ch,1)
    rebnconv2d = rebnconv_upsample(mid_ch*2,mid_ch,1)

    rebnconv1d = Chain(rebnconv(mid_ch*2,out_ch,1)...)

    layer = x -> begin
        hxin = rebnconvin(x)

        hx1, hx = rebnconv1(hxin)
        hx2, hx = rebnconv2(hx)
        hx3, hx = rebnconv3(hx)
        hx4, hx = rebnconv4(hx)

        hx5 = rebnconv5(hx)
        hx6 = rebnconv6(hx5)

        hx5dup = rebnconv5d(hx4, hx5, hx6)
        hx4dup = rebnconv4d(hx3, hx4, hx5dup)
        hx3dup = rebnconv3d(hx2, hx3, hx4dup)
        hx2dup = rebnconv2d(hx1, hx2, hx3dup)

        hx1d = rebnconv1d(cat(hx2dup, hx1; dims=3))

        return hx1d + hxin
    end
end

function rsu5(in_ch=3, mid_ch=12, out_ch=3)
    rebnconvin = Chain(rebnconv(in_ch,out_ch,1)...)

    rebnconv1 = rebnconv_pool(out_ch, mid_ch, 1)
    rebnconv2 = rebnconv_pool(mid_ch, mid_ch, 1)
    rebnconv3 = rebnconv_pool(mid_ch, mid_ch, 1)

    rebnconv4 = Chain(rebnconv(mid_ch,mid_ch,1)...)
    rebnconv5 = Chain(rebnconv(mid_ch,mid_ch,2)...)

    rebnconv4d = rebnconv_upsample(mid_ch*2,mid_ch,1)
    rebnconv3d = rebnconv_upsample(mid_ch*2,mid_ch,1)
    rebnconv2d = rebnconv_upsample(mid_ch*2,mid_ch,1)

    rebnconv1d = Chain(rebnconv(mid_ch*2,out_ch,1)...)

    layer = x -> begin
        hxin = rebnconvin(x)

        hx1, hx = rebnconv1(hxin)
        hx2, hx = rebnconv2(hx)
        hx3, hx = rebnconv3(hx)

        hx4 = rebnconv4(hx)
        hx5 = rebnconv5(hx4)

        hx4dup = rebnconv4d(hx3, hx4, hx5)
        hx3dup = rebnconv3d(hx2, hx3, hx4dup)
        hx2dup = rebnconv2d(hx1, hx2, hx3dup)

        hx1d = rebnconv1d(cat(hx2dup, hx1; dims=3))

        return hx1d + hxin
    end
end

function rsu4(in_ch=3, mid_ch=12, out_ch=3)
    rebnconvin = Chain(rebnconv(in_ch,out_ch,1)...)

    rebnconv1 = rebnconv_pool(out_ch, mid_ch, 1)
    rebnconv2 = rebnconv_pool(mid_ch, mid_ch, 1)

    rebnconv3 = Chain(rebnconv(mid_ch,mid_ch,1)...)
    rebnconv4 = Chain(rebnconv(mid_ch,mid_ch,2)...)

    rebnconv3d = rebnconv_upsample(mid_ch*2,mid_ch,1)
    rebnconv2d = rebnconv_upsample(mid_ch*2,mid_ch,1)

    rebnconv1d = Chain(rebnconv(mid_ch*2,out_ch,1)...)

    layer = x -> begin
        hxin = rebnconvin(x)

        hx1, hx = rebnconv1(hxin)
        hx2, hx = rebnconv2(hx)

        hx3 = rebnconv3(hx)
        hx4 = rebnconv4(hx3)

        hx3dup = rebnconv3d(hx2, hx3, hx4)
        hx2dup = rebnconv2d(hx1, hx2, hx3dup)

        hx1d = rebnconv1d(cat(hx2dup, hx1; dims=3))

        return hx1d + hxin
    end
end

function rsu4f(in_ch=3, mid_ch=12, out_ch=3)
    rebnconvin = Chain(rebnconv(in_ch,out_ch,1)...)

    rebnconv1 = Chain(rebnconv(out_ch,mid_ch,1)...)
    rebnconv2 = Chain(rebnconv(mid_ch,mid_ch,2)...)
    rebnconv3 = Chain(rebnconv(mid_ch,mid_ch,4)...)

    rebnconv4 = Chain(rebnconv(mid_ch,mid_ch,8)...)

    rebnconv3d = Chain(rebnconv(mid_ch*2,mid_ch,4)...)
    rebnconv2d = Chain(rebnconv(mid_ch*2,mid_ch,2)...)

    rebnconv1d = Chain(rebnconv(mid_ch*2,out_ch,1)...)

    layer = x -> begin
        hxin = rebnconvin(x)

        hx1 = rebnconv1(hxin)
        hx2 = rebnconv2(hx1)
        hx3 = rebnconv3(hx2)

        hx4 = rebnconv4(hx3)

        hx3d = rebnconv3d(cat(hx4, hx3; dims=3))
        hx2d = rebnconv2d(cat(hx3d, hx2; dims=3))
        hx1d = rebnconv1d(cat(hx2d, hx1; dims=3))

        return hx1d + hxin
    end
end


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

export get_uttID2file
function get_uttID2file(dataset_path,folder::String)
  uttID2file = Dict{String,String}()
  for (root,dir,files) in walkdir(joinpath(dataset_path,folder);follow_symlinks=true)
    wavs = files[findall(contains.(files,".wav"))]
    folders = split(root, "/")
    spkID = folders[end]
    type = folders[end-1]
    for f in wavs
      sentenceID = split(f, "."; limit=2)[1]
      uttID = "$(spkID)-$(type)-$(sentenceID)"
      uttID2file[uttID] = joinpath(root,f)
    end
  end
  return uttID2file
end

export get_uttID2text
function get_uttID2text(uttID2file::Dict)
  d = Dict(
           'z' => "ZERO", '3' => "THREE", '7' => "SEVEN",
           'o' => "OH",   '4' => "FOUR",  '8' => "EIGHT",
           '1' => "ONE",  '5' => "FIVE",  '9' => "NINE",
           '2' => "TWO",  '6' => "SIX",   'a' =>"",  'b'=>"")
  uttID2text = Dict{String,String}()
  for uttID in keys(uttID2file)
    text = split(uttID,"-")[3]
    try
      uttID2text[uttID] = strip(prod([d[t] for t in text].*" "))
    catch
      error("$text is an invalid filename, invalid dataset!")
    end
  end
  return uttID2text
end

function get_uttID2wav(uttID2file::Dict; T=Float32)
  uttID2wav = Dict(uttID => T.(load(uttID2file[uttID]).data)[:]
                   for uttID in keys(uttID2file))
end


text:
using DataFormatter
using DataFrames
using OrdinaryDiffEq
using NamedArrays

# functions defined to allow simulation of DataSets from DataFormatter.jl

function simulate_ponctual(p, data_params, target)
    if target == [:Δσ_peak]
        df_sol = get_peak_values(p, data_params)
    elseif target == [:ϵ_peak, :Δσ_peak]
        df_sol = get_peak_ϵ_and_Δσ(p, data_params)
    end
    return df_sol
end

function simulate_timeseries(p, data_params, target, time_vec=nothing)
    if any(target .∈ Ref((:Δσ,)))
        df_sol = integrate_csr(p, data_params ; time_vec)
    elseif any(target .∈ Ref((:ϵ̇_dev,)))
        #model = build_model(mp, data.pc ; control = ConstantStress(data.Δσ), geom = p.geom)
    end
    return df_sol
end

function get_peak_values(p, data)
    df = integrate_csr(p, data ; stop_at_peak=true)
    row_peak = df[argmin(df.Δσ_sim),:]
    return DataFrame(["Δσ_peak_sim", "D_peak_sim"] .=> [row_peak."Δσ_sim",row_peak."D_sim"])
end

function get_peak_ϵ_and_Δσ(p, data)
    df = integrate_csr(p, data ; stop_at_peak=true)
    row_peak = df[end,:] #argmax(abs.(df.Δσ_sim))
    ϵ_peak = row_peak.t * data.ϵ̇
    Δσ_peak = row_peak.Δσ_sim
    D_peak = row_peak.D_sim
    return DataFrame(["ϵ_peak_sim", "Δσ_peak_sim", "D_peak_sim"] .=> [ϵ_peak, Δσ_peak, D_peak])
end

function integrate_csr(p, data ; time_vec=Float64[], stop_at_peak=false)
    # define model
    model = build_model(p, data)

    #define derivative update functions
    update_func!(du,u,p,t) = update_derivatives!(du,u,p,t,model)

    # get initial conditions and tspan
    Dᵢ = (:Dᵢ ∈ names(p.mp)[1]) ? p.mp[:Dᵢ] : getp(:D₀,p.mp,p.mp_default)
    u0 = init_variables(model, Dᵢ)
    tspan = (0.0, p.model.ϵmax/data.ϵ̇_dev)

    # create callbacks
    conditionD(u,t,i) = u[2] - p.model.Dmax # if damage overshoots 1
    affectD!(integrator) = (integrator.u[2] = p.model.Dmax) # set it to 1
    cbD = ContinuousCallback(conditionD, affectD!, save_positions=(false,false)) # check that continuously

    condition_last_time(u,t,i) = t > tspan[2]
    cb_last_time = DiscreteCallback(condition_last_time, terminate! ; save_positions=(false,false))

    if stop_at_peak
        conditionpeak(u,t,i) = (abs(i.u[1]) < abs(i.uprev[1])) # if axial stress decreases
        cbpeak = DiscreteCallback(conditionpeak, terminate! ; save_positions=(true,true))
        cb = CallbackSet(cbD, cbpeak, cb_last_time)
    else
        cb = CallbackSet(cbD, cb_last_time)
    end

    # create ODEProblem and solve
    prob = ODEProblem(update_func!, convert.(eltype(p.mp), u0), tspan, p.model)
    sol = solve(prob, AutoTsit5(Rosenbrock23());#Rosenbrock23();#Tsit5();
            abstol=p.solver.abstol,
            reltol=p.solver.reltol,
            callback = cb,
            maxiters=p.solver.maxiters,
            saveat=time_vec)

    # find s_peak and convert to Δσ
    s_vec = sol[1,:]
    Δσ_vec = sₐₓ2Δσ.(Ref(model),s_vec)
    return DataFrame(["t","Δσ_sim", "D_sim"].=> [sol.t, Δσ_vec, sol[2,:]])
end

function build_model(p, data)
    mp = p.mp
    mpd = p.mp_default
    geom = p.model.geom
    DamType = p.model.damage_type
    pc = data.pc
    ϵ̇ = data.ϵ̇_dev

    r = Rheology(; μ=getp(:μ,mp,mpd), ψ=getp(:ψ,mp,mpd), a=getp(:a,mp,mpd), D₀=getp(:D₀,mp,mpd), n=getp(:n,mp,mpd), K₁c=getp(:K₁c,mp,mpd), l̇₀=getp(:l̇₀,mp,mpd))
    plasticity = Plasticity(MinViscosityThreshold(), CoulombYieldStress(μ=getp(:μ,mp,mpd),C=0))
    cm = ConstitutiveModel(;weakening = LinearWeakening(γ=getp(:γ,mp,mpd)),
                       damage = DamType(r),
                       elasticity = IncompressibleElasticity(G=getp(:G,mp,mpd)), #3.146e9
                       plasticity)
    setup = TriaxialSetup(;geom,
                           control = ConstantStrainRate(ϵ̇),
                           pc,
                           )
    return Model(;setup, cm)
end

getp(sym::Symbol,A::NamedArray,B::NamedArray) = haskey(A.dicts[1],sym) ? A[sym] : B[sym]

max_stars_repo_name: Leooop/SimpleDamage

max_issues_repo_name: Leooop/SimpleDamage

max_forks_repo_name: Leooop/SimpleDamage


text:
module test

using Random
using Dates

using Backtest.Ids: AssetId
using Backtest.Engine: CalcLattice, addfields!, newbar!
using Backtest.ConcreteFields: Open, High, Low, Close, Volume, SMA, ZScore, Rank, Returns, LogReturns
using Backtest.Events: FieldCompletedProcessingEvent, AbstractOrderEvent
using Backtest.DataReaders: InMemoryDataReader, PerFileDataReader, AbstractDataReader
using Backtest.Orders: MarketOrder, LimitOrder
using Backtest: Strategy, StrategyOptions, run, order!, log, INFO, TRANSACTIONS, NOVERBOSITY
using Backtest.Utils: crossover, crossunder, writejson

#
# function enginetest()
#
#   # Make a backtest object
#   nbarstostore = 5
#   assetids = ["AAPL", "MSFT", "TSLA"]
#   lattice = CalcLattice(nbarstostore, assetids)
#
#   function rand_(base::Number)
#     SEED = 1234
#     rng = MersenneTwister() # Optional: include SEED for reproducibility
#     result = base*(1+randn(rng))
#     return result
#   end
#
#   # Make data
#   nbarstorun = max(nbarstostore, 120_960)
#   allbars = [
#     Dict(
#       "AAPL" => Dict("Open" => rand_(1.0*i), "Close" => rand_(1.2*i)),
#       "MSFT" => Dict("Open" => rand_(2.0*i), "Close" => rand_(2.5*i)),
#       "TSLA" => Dict("Open" => rand_(3.0*i), "Close" => rand_(2.5*i))
#     ) for i in 1:nbarstorun
#   ]
#   # println(allbars)
#   # println()
#   # println(allbars[6][AssetId("TSLA")][FieldId("Open")])
#
#   # Make fields
#   fieldoperations = Vector([
#     Open("Open"),
#     Close("Close"),
#     SMA("SMA2-Open", "Open", 2),
#     SMA("SMA3-Close", "Close", 3),
#     ZScore("ZScore-Open", "Open")
#   ])
#   # print("Fields: ")
#   # println(fields)
#
#   addfields!(lattice, fieldoperations)
#
#   # print("Dict afterward: ")
#   # println(lattice.windowdependentfields)
#
#
#   # Feed bars into engine
#   # (for benchmarking purposes) add first bar separately so that compile time is not included in the benchmark
#   newbar!(lattice, allbars[1]);
#
#   @time begin
#   for barindex in 2:length(allbars)
#     newbar!(lattice, allbars[barindex])
#   end end
#
#   nbarstoprint = 0
#   nbars = length(lattice.recentbars)
#   for i in (nbars-nbarstoprint+1):nbars
#     println(lattice.recentbars[i])
#   end
# end
#
function getexamplecryptodatareader(symbol::String)
  # basepath = "../fakedata/stocks/minute/$symbol"
  basepath = "../fakedata/crypto/$symbol"
  daysofmonth = [lpad(i, 2, "0") for i=1:10]
  # daysofmonth = ["06", "07", "08", "09", "10", "13", "14", "15", "16", "17", "18", "19", ""]
  dayofmonthtopath = (dom -> "$basepath/2015-07-$(dom).csv")
  sources = [dayofmonthtopath(dom) for dom in daysofmonth]
  # return InMemoryDataReader(symbol, sources, datetimecol="Datetime")
  return PerFileDataReader(symbol, sources, datetimecol="Datetime")
end

function datareadertest()
  # Test that this doesn't error out.
  getexamplecryptodatareader("BTCUSDT")
  getexamplecryptodatareader("LTCUSDT")
  getexamplecryptodatareader("ETHUSDT")
end

function bttest()
  ## Make the data readers ##
  assetids = ["A", "B", "C"]
  datareaders = Dict{AssetId, AbstractDataReader}()
  for assetid in assetids
    datareaders[assetid] = getexamplecryptodatareader(assetid)
  end

  ## Make the field operations ##
  (dt, o, h, l, c, v) = ("Datetime", "Open", "High", "Low", "Close", "Volume")
  fieldoperations = [
    SMA("SMA30-Close", c, 30),
    SMA("SMA60-Close", c, 60),
    ZScore("ZScore-Open", o),
    Rank("Rank-Close", c),
    Returns("Returns-Close", c),
    Returns("Returns-Close-3", c, 3),
    LogReturns("LogReturns-Close", c),
    LogReturns("LogReturns-Close-3", c, 3),
  ]

  # function ondata(strat::Strategy, event::FieldCompletedProcessingEvent)
  #   # order!(strat, MarketOrder("LTCUSDT", 1))
  #   # order!(strat, LimitOrder("LTCUSDT", -1, 10))
  #   fast, slow = "SMA30-Close", "SMA60-Close"
  #   if crossover(strat, "A", fast, slow)
  #     order!(strat, MarketOrder("A", .1))
  #   elseif crossunder(strat, "A", fast, slow)
  #     order!(strat, MarketOrder("A", -.1))
  #   end
  # end
  #
  # function onorder(strat::Strategy, event::ET) where {ET<:AbstractOrderEvent}
  # end
  #
  # stratoptions = StrategyOptions(
  #   datareaders=datareaders,
  #   fieldoperations=fieldoperations,
  #   numlookbackbars=-1,
  #   start=Dates.DateTime(2015, 7, 1, 0, 0),
  #   endtime=Dates.DateTime(2015, 7, 10, 0, 0),
  #   tradinginterval=Dates.Minute(1),
  #   verbosity=NOVERBOSITY,
  #   datadelay=Dates.Second(5),
  #   messagelatency=Dates.Second(3),
  #   fieldoptimeout=Dates.Second(2),
  #   datetimecol=dt,
  #   opencol=o,
  #   highcol=h,
  #   lowcol=l,
  #   closecol=c,
  #   volumecol=v,
  #   ondataevent=ondata,
  #   onorderevent=onorder,
  #   principal=1_000
  # )
  #
  # # Run the backtest ##
  # @time begin
  #   run(stratoptions)
  # end

  # Test running a data-only job
  jsonfile = "jsondata.json"
  @time writejson(jsonfile;
    datareaders=datareaders,
    fieldoperations=fieldoperations,
    start=Dates.DateTime(2015, 7, 8, 10, 0),
    endtime=Dates.DateTime(2015, 7, 8, 10, 2),
    tradinginterval=Dates.Minute(1),
    datetimecol=dt,
    opencol=o,
    highcol=h,
    lowcol=l,
    closecol=c,
    volumecol=v,
  )
end
end

max_stars_repo_name: maxholloway/backtest

max_issues_repo_name: maxholloway/backtest

max_forks_repo_name: maxholloway/backtest


text:
#=
Check that the Poisson solver in the PIC framework works as expected.
It is used in ``PIC.compute_phi``
=#

using Test

@testset "Poisson solve cos" begin
    using LinearAlgebra, Roots
    using VlasovSolvers
    # using Plots

    include("PIC.jl")

    function samples_verif2(nsamples)
        x0 = Array{Float64}(undef, nsamples)
        seq = SobolSeq(1)
        for i = 1:nsamples
            r = Sobol.next!(seq)
            resx = find_zero(x-> x + 0.99 / kx * sin(kx*x) - r[1]*2π/kx, 0.0)
            x0[i] = resx
        end
        v0 = rand(Uniform(-10, 10), nsamples)
        wei = 2π/kx/nsamples .* ones(Float64, nsamples)
        return x0, v0, wei
    end

    kx = 0.5
    L = 2π / kx
    mesh1 = OneDGrid(CPU(), nx, 0, L)

    ρ = cos.(kx .* mesh1.points)
    ρ_tot = 0
    ϕ = compute_phi(ρ, ρ_tot, mesh1)

    # plot(mesh1.points, ρ, label="ρ")
    # plot!(mesh1.points, ϕ .* kx^2, label="ϕ")
    @test maximum(abs.(ϕ  .- ρ ./ kx^2)) < 10^-3
end

@testset "Poisson solve sin" begin
    using LinearAlgebra, Roots
    using VlasovSolvers
    # using Plots

    include("PIC.jl")

    function samples_verif2(nsamples)
        x0 = Array{Float64}(undef, nsamples)
        seq = SobolSeq(1)
        for i = 1:nsamples
            r = Sobol.next!(seq)
            resx = find_zero(x-> x + 0.99 / kx * sin(kx*x) - r[1]*2π/kx, 0.0)
            x0[i] = resx
        end
        v0 = rand(Uniform(-10, 10), nsamples)
        wei = 2π/kx/nsamples .* ones(Float64, nsamples)
        return x0, v0, wei
    end

    nx = 256
    kx = 0.4
    L = 2π / kx
    mesh1 = OneDGrid(CPU(), nx, 0, L)

    ρ = sin.(5*kx .* mesh1.points)
    ρ_tot = 0
    ϕ = compute_phi(ρ, ρ_tot, mesh1)

    # plot(mesh1.points, ρ, label="ρ")
    # plot!(mesh1.points, ϕ .* 25kx^2, label="ϕ")
    @test maximum(abs.(ϕ  .- ρ ./ (25*kx^2))) < 10^-3
end


text:
# ---
# title: 1013. Partition Array Into Three Parts With Equal Sum
# id: problem1013
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/partition-array-into-three-parts-with-equal-sum/description/>
# hidden: true
# ---
#
# Given an array `A` of integers, return `true` if and only if we can partition
# the array into three **non-empty** parts with equal sums.
#
# Formally, we can partition the array if we can find indexes `i+1 < j` with
# `(A[0] + A[1] + ... + A[i] == A[i+1] + A[i+2] + ... + A[j-1] == A[j] + A[j-1]
# + ... + A[A.length - 1])`
#
#
#
# **Example 1:**
#
#
#
#     Input: A = [0,2,1,-6,6,-7,9,1,2,0,1]
#     Output: true
#     Explanation: 0 + 2 + 1 = -6 + 6 - 7 + 9 + 1 = 2 + 0 + 1
#
#
# **Example 2:**
#
#
#
#     Input: A = [0,2,1,-6,6,7,9,-1,2,0,1]
#     Output: false
#
#
# **Example 3:**
#
#
#
#     Input: A = [3,3,6,5,-2,2,5,1,-9,4]
#     Output: true
#     Explanation: 3 + 3 = 6 = 5 - 2 + 2 + 5 + 1 - 9 + 4
#
#
#
#
# **Constraints:**
#
#   * `3 <= A.length <= 50000`
#   * `-10^4 <= A[i] <= 10^4`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


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

using Markdown
using InteractiveUtils

# ╔═╡ d09d95a0-6b17-11eb-11dc-37663c2e26ed
using ChemEquations

# ╔═╡ bfdf0fbe-6b17-11eb-1c12-c758f0464228
import Pkg

# ╔═╡ c48abd6a-6b17-11eb-36f5-5f3ba22cc688
Pkg.add("ChemEquations")

# ╔═╡ 3cf182e8-6b13-11eb-0207-39526e08850f
md" ## _Solubility Equilibria_
The purpose of this notebook is to explore solubility, specifically the reaction system and how changes will affect the overall system. The topics that will be covered are:
1. *Ksp*
2. *Qsp*
3. *Saturation*
4. *Predicting the formation of a precipitiate*
5. *Common Ion Effect*
6. *ΔGᴼ*

"

# ╔═╡ 33092d3c-6b16-11eb-0701-91cce981323e
md" ## Ksp
**Solubility**: Amount of salt of salt that will dissolve to form a saturated solution at a given temperature.

Depending on the type of salt, the value for solubility changes. Salts that have low solubility are known as *insoluble salts*. This includes salts like silver chloride or calcium sulfate. Meanwhile, salts with high solubility, like Sodium Chloride or Potassium suflate are know as *soluble salts*.

How do you calculate *Ksp*, the solubility product constant?

Lets use this equation an example!

$\ce{PbCl2(s) -> Pb+2(aq) + 2Cl-(aq) }$
"

# ╔═╡ a3c5dbe6-6b17-11eb-3f3a-290e677b30d0
equation = ce" PbCl2 = Pb + Cl"

# ╔═╡ 0ce2aee2-6b18-11eb-1d80-3526c27ffe80
balance(equation)

# ╔═╡ 9265a3ce-6fe1-11eb-325f-1bf2c77e36f6
md" Ksp is another equilibrium expression, specifically for soluble salts. Therefore, the equation for Ksp is
$\ce{Ksp =[Pb+2][Cl-]^2}$

In a normal equilibrium expression, we include the reactants, however, since the reactant is a solid, it doesn't get included into the expression, hence the reason the expression only has the products of the reaction

"

# ╔═╡ 8decce0e-6fe2-11eb-0226-2f4ca2954e94
md" Here is a Ksp Calculator. Plug in the concentrations of the products (the aqueous components) and the exponents"

# ╔═╡ bef3f120-6fe3-11eb-3ce0-7d2de07c388c
Ksp(c1,e1,c2,e2) = (c1^e1)*(c2^e2)

# ╔═╡ ea82fa7e-6fe4-11eb-2648-01ecd4757558
Ksp(0.08967,3,0.6789,2)

# ╔═╡ 691f1508-6b18-11eb-2f11-8354277afae2
md" ## Qsp
Qsp, or Q, is the number of ions in solution at any given moment. It is calculated the same way as Ksp.
"

# ╔═╡ 3f96aba4-6fe6-11eb-2c80-d988d9209ea7
Qsp(c,e,C,E)=(c^e)*(C^E)

# ╔═╡ 543b6202-6fe6-11eb-0969-8d26d5c50f6c
Qsp(0.067,1, 0.00324,2)

# ╔═╡ 67bad146-6fe6-11eb-2742-63398669caa5
md" The difference between Ksp and Qsp is that Ksp is the number of ions in solution at the saturation point whereas Qsp is the number of ions at any time.
"

# ╔═╡ 5260d00c-6fe1-11eb-2dd1-39e92ef68a07
md" ## Saturation
"

# ╔═╡ d0fed526-6fe8-11eb-23df-65103a653151
precipitate(Q,Ksp1)
if Qsp>Ksp1
	md"The solution is supersaturated. This means there are too many ions so by Le Chatliers principle, the reaction shifts to the left(towards the solid), creating more precipitate in the solution. Currently in the beaker, there is solid at the bottom of the beaker."
elseif Qsp < Ksp1
		md"The solution is not saturated. The solid can still dissolve. The only thing in the beaker is ions"
else
		md"Since Q and K are equal,the solution is saturated. This means that the rate of dissolution(the forward reaction) and the reverse reaction(formation of the precipipate) are equal. Currently in the beaker, there is a small amount of crystallization."
end

# ╔═╡ 1147b294-6fea-11eb-12a3-15c40863664e
precipitate(0.01,0.02)

# ╔═╡ 59151304-6fe1-11eb-09ee-0f3da19ba7cb
md" ## Formation of a precipitate
"

# ╔═╡ 6322a014-6fe1-11eb-1cad-6d6aefa75fde
md" ## Common Ion Effect
"

# ╔═╡ Cell order:
# ╠═bfdf0fbe-6b17-11eb-1c12-c758f0464228
# ╠═c48abd6a-6b17-11eb-36f5-5f3ba22cc688
# ╠═d09d95a0-6b17-11eb-11dc-37663c2e26ed
# ╟─3cf182e8-6b13-11eb-0207-39526e08850f
# ╟─33092d3c-6b16-11eb-0701-91cce981323e
# ╠═a3c5dbe6-6b17-11eb-3f3a-290e677b30d0
# ╠═0ce2aee2-6b18-11eb-1d80-3526c27ffe80
# ╟─9265a3ce-6fe1-11eb-325f-1bf2c77e36f6
# ╟─8decce0e-6fe2-11eb-0226-2f4ca2954e94
# ╠═bef3f120-6fe3-11eb-3ce0-7d2de07c388c
# ╠═ea82fa7e-6fe4-11eb-2648-01ecd4757558
# ╟─691f1508-6b18-11eb-2f11-8354277afae2
# ╠═3f96aba4-6fe6-11eb-2c80-d988d9209ea7
# ╠═543b6202-6fe6-11eb-0969-8d26d5c50f6c
# ╠═67bad146-6fe6-11eb-2742-63398669caa5
# ╠═5260d00c-6fe1-11eb-2dd1-39e92ef68a07
# ╠═d0fed526-6fe8-11eb-23df-65103a653151
# ╠═1147b294-6fea-11eb-12a3-15c40863664e
# ╠═59151304-6fe1-11eb-09ee-0f3da19ba7cb
# ╠═6322a014-6fe1-11eb-1cad-6d6aefa75fde


text:
function Nc(T,S::SemiconductorMaterial)
    T_0 = S.T0
    S.Nc300*(T/T_0)^(3/2)
end

function Nv(T,S::SemiconductorMaterial)
    T_0 = S.T0
    S.Nv300*(T/T_0)^(3/2)
end

function ni(T,f::DistributionFunc,S::SemiconductorMaterial)
    sqrt(n(T,f,S)*p(T,f,S))
end

function N(T,f::DistributionFunc,S::SemiconductorMaterial)
    return N(T,Ef(T,f,S),f,S)
end

function N(T,E_f,f::DistributionFunc,S::SemiconductorMaterial)
    E_v = Ev(T,S)
    return n(T,E_f,f,S)+Na(S.Dopping,T,E_f,0,E_v)
end

function Pt(T,f::DistributionFunc,S::SemiconductorMaterial)
    return Pt(T,Ef(T,f,S),f,S)
end
function Pt(T,E_f,f::DistributionFunc,S::SemiconductorMaterial)
    E_c = Ec(T,S)
    return p(T,E_f,f,S)+Na(S.Dopping,T,E_f,E_c,0.0)
end

function Nt(T,f::DistributionFunc,S::SemiconductorMaterial)
    return Nt(T,Ef(T,f,S),f,S)
end
function Nt(T,E_f,f::DistributionFunc,S::SemiconductorMaterial)
    E_c = Ec(T,S)
    return n(T,E_f,f,S)+Nd(S.Dopping,T,E_f,E_c,0)
end

function n(T,f::DistributionFunc,S::SemiconductorMaterial)
    return n(T,Ef(T,f,S),Ec(T,S),f,S)
end

function n(T,E_f,f::DistributionFunc,S::SemiconductorMaterial)
    return n(T,E_f,Ec(T,S),f,S)
end

function n(T,E_f,E_C,f::DistributionFunc,S::SemiconductorMaterial)
    ΔE = E_f - E_C
    return Nc(T,S)*f(ΔE,T)
end

function P(T,f::DistributionFunc,S::SemiconductorMaterial)
    η = Ef(T,f,S)
    return P(T,η,f,S)
end

function P(T,E_f,f::DistributionFunc,S::SemiconductorMaterial)
    E_c = Ec(T,S)
    return p(T,E_f,f,S)+Nd(S.Dopping,T,E_f,E_c,0)
end

function p(T,f::DistributionFunc,S::SemiconductorMaterial)
    η = Ef(T,f,S)
    return p(T,η,f,S)
end

function p(T,E_f,f::DistributionFunc,S::SemiconductorMaterial)
    p(T,E_f,Ev(T,S),f,S)
end

function p(T,E_f,E_V,f::DistributionFunc,S::SemiconductorMaterial)
    ΔE = E_V-E_f
    return Nv(T,S)*f(ΔE,T)
end


text:
module SineFit

import LsqFit
import FFTW
import Statistics

export estimate_wave_parameters, calculate_wave_shape, sin

struct WaveFitParams
    amplitude::Float64
    frequency::Float64
    phaseshft::Float64
    vertoffst::Float64
end

"""
    estimate_wave_parameters(xdata::Vector{T}, ydate::Vector{T};
        sample_rate=nothing)::WaveFitParams where {T<:Number}

Use Fourier analysis and statistics to find the best monochromatic sine wave parameters for data.
`xsym` and `ysym` are the column symbols for the x and y columns within DataFrame `df`, respectively.
Returns a `WaveFitParams` struct with frequency and phase shift as circular frequencies resp. radians.

If `sample_rate` is not specified, the rate is determined automatically from the `xdata` argument.
"""
function estimate_wave_parameters(xdata::Vector{T}, ydata::Vector{T};
        sample_rate=nothing)::WaveFitParams where {T<:Number}
    if sample_rate === nothing
        sample_rate = 1 / Statistics.mean(xdata[2:end] - xdata[1:end-1])
    end

    offset = (maximum(ydata)+minimum(ydata))/2
    ampl = Statistics.mean(Statistics.quantile(abs.(ydata .- offset), 99/100))
    fft = FFTW.rfft(ydata .- offset)
    fftfreq = FFTW.rfftfreq(length(ydata), sample_rate)

    # Determine dominant frequency and calculate phase shift.
    maxfreqloc = argmax(abs.(fft))
    maxfreq = fftfreq[maxfreqloc]
    shift = mod(atan(imag(fft[maxfreqloc]), real(fft[maxfreqloc])) + pi/2, 2pi)
    freq = maxfreq * 2pi

    WaveFitParams(ampl, freq, shift, offset)
end

function sin_model(xs, params)
    params[1] .* Base.sin.(params[2] .* xs .+ params[3]) .+ params[4]
end

"""
    calculate_wave_shape(xdata::Vector{T}, ydate::Vector{T};
        sample_rate=nothing)::WaveFitParams where {T<:Number}

Fit a sine function (`sin_model`, i.e. `a sin(b x + c) + d`) to the given data `xdata`/`ydata`.

This function first estimates the parameters using normal statistics and fourier analysis using `estimate_wave_parameters()`.
Then, a least-squares fit (`LsqFit`) is used to refine the parameters -- for example, to adjust the frequency, as a FFT has
a very limited frequency resolution. This helps a lot with accuracy.

If `sample_rate` is not specified, the rate is determined automatically from the `xdata` argument.
"""
function calculate_wave_shape(xdata::Vector{T}, ydata::Vector{T};
        sample_rate=nothing)::WaveFitParams where {T<:Number}

    estimate = estimate_wave_parameters(xdata, ydata, sample_rate=sample_rate)
    initial = [estimate.amplitude, estimate.frequency, estimate.phaseshft, estimate.vertoffst]

    fit = LsqFit.curve_fit(sin_model, xdata, ydata, initial)
    params = fit.param

    if params[1] < 0
        params[1] *= -1
        params[3] -= 180
    end
    params[3] = mod(params[3], 360)
    if params[3] > 180
        params[3] = params[3] - 360
    elseif params[3] < -180
        params[3] = params[3] + 360
    end

    WaveFitParams(params...)
end

"""
    sin(x::Vector{T}, wp::WaveFitParams)::Vector{T} where {T<:Number}

Return sin for `x` with the parameters given in `wp`. Useful in order to plot the fit function, for example.

```
wp = calculate_wave_shape(xs, ys)
Plots.plot(xs, ys)
Plots.plot!(xs, sin.(xs, wp))
```

"""
function sin(xs::Vector{T}, wp::WaveFitParams)::Vector{T} where {T<:Number}
    wp.amplitude .* Base.sin.(wp.frequency .* xs .+ wp.phaseshft) .+ wp.vertoffst
end

end # module


text:
function ellipse!(μ,Σ;kw...)
    for i=1:3
        plot!([tuple((i*sqrt(Σ)*[cos(t),sin(t)]+μ)...) for t in range(0,stop=2π,length=100)];
              seriestype=:shape,linealpha=0.5,linecolor=:gray,fillalpha=[0.4,0.2,0.1][i],kw...)
    end
    current()
end

function mixtureplot(X₁,X₂,μ₀,Σ₀,μ₁,Σ₁,Π=nothing)
    scatter(X₁,X₂;fillalpha=0.5,markerstrokewidth=0.5,
        marker_z=Π,mc=ColorGradient([:lightblue,:orange]),
        colorbar=:false)
    ellipse!(μ₀,Σ₀,fillcolor=:lightblue)
    ellipse!(μ₁,Σ₁,fillcolor=:orange)
    plot!(;bg=:transparent,xlims=(-4.5,8.5),ylims=(-3.5,10.5),
        leg=false,ticks=:none,ratio=:equal,bgcolor=:white)
end

max_stars_repo_name: data-gymnasia/julia-binder

max_issues_repo_name: data-gymnasia/julia-binder

max_forks_repo_name: data-gymnasia/julia-binder


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    DRAW_SETTING,
    atomcolor, isatomvisible, sdfcoords2d, coords2d,
    chargesign, atommarkup, atomhtml,
    draw2d!, drawatomindex!, sethighlight!

"""
    DRAW_SETTING

Default setting parameters of the molecule drawing canvas.

# Required fields

- `:display_terminal_carbon`(Bool) whether to display terminal C atom or not
- `:double_bond_notation`(Symbol)
    - `:alongside`: all double bonds are represented as a carbon skeleton and
                    a segment alongside it.
    - `:dual`: all double bonds are represented as two equal length parallel
               segments.
    - `:chain`: `:dual` for chain bonds and `:alongside` for ring bonds
    - `:terminal`: `:dual` for terminal bonds (adjacent to degree=1 node) and
                   `:alongside` for others (default)
- `:atomcolor`(Dict{Symbol,Color}) atom symbol and bond colors for organic atoms
- `:defaul_atom_color`(Dict{Symbol,Color}) colors for other atoms
"""
const DRAW_SETTING = Dict(
    :display_terminal_carbon => false,
    :double_bond_style => :terminal,  # :alongside, :dual, :chain, :terminal
    :atomcolor => Dict(
        :H => Color(0, 0, 0),
        :B => Color(128, 0, 0),
        :C => Color(0, 0, 0),
        :N => Color(0, 0, 255),
        :O => Color(255, 0, 0),
        :F => Color(0, 255, 0),
        :Si => Color(128, 64, 192),
        :P => Color(192, 0, 192),
        :S => Color(192, 192, 0),
        :Cl => Color(64, 192, 64),
        :As => Color(128, 0, 128),
        :Se => Color(128, 128, 0),
        :Br => Color(0, 192, 0),
        :I => Color(0, 128, 0)
    ),
    :default_atom_color => Color(0, 192, 192)
)

"""
    atomcolor(mol::GraphMol; setting=DRAW_SETTING) -> Vector{Color}

Return atom colors for molecule 2D drawing
"""
@cachefirst function atomcolor(
        mol::GraphMol; setting=DRAW_SETTING, kwargs...)
    atomc = setting[:atomcolor]
    dfc =  setting[:default_atom_color]
    return [get(atomc, sym, dfc) for sym in atomsymbol(mol)]
end

atomcolor(view::SubgraphView; kwargs...) = atomcolor(view.graph; kwargs...)

"""
    isatomvisible(mol::GraphMol; setting=DRAW_SETTING) -> Vector{Bool}

Return whether the atom is visible in the 2D drawing.
"""
@cachefirst function isatomvisible(
        mol::GraphMol; setting=DRAW_SETTING, kwargs...)
    termc = setting[:display_terminal_carbon]
    vec = Bool[]
    deg_ = nodedegree(mol)
    sym_ = atomsymbol(mol)
    for i in 1:nodecount(mol)
        isvisible = deg_[i] == 0 || sym_[i] !== :C || (termc && deg_[i] == 1)
        push!(vec, isvisible)
    end
    return vec
end

isatomvisible(view::SubgraphView; kwargs...
    ) = isatomvisible(view.graph; kwargs...)

@cachefirst function sdfcoords2d(mol::SDFile)
    coords = zeros(Float64, nodecount(mol), 2)
    for (i, node) in enumerate(nodeattrs(mol))
        coords[i, :] = node.coords[1:2]
    end
    return coords
end

sdfcoords2d(view::SubgraphView) = sdfcoords2d(view.graph)

@cachefirst function coords2d(
        mol::GraphMol; forcecoordgen=false, setting=DRAW_SETTING, kwargs...)
    if nodeattrtype(mol) === SDFileAtom && !forcecoordgen
        coords = sdfcoords2d(mol)
        style = getproperty.(edgeattrs(mol), :notation)
    else
        coords, style = coordgen(mol)
    end
    bondorder_ = bondorder(mol)

    # All double bonds to be "="
    if setting[:double_bond_style] == :dual
        for b in findall(bondorder_ .== 2)
            style[b] = 2
        end
        return
    end

    # Or only non-ring bonds to be "="
    if setting[:double_bond_style] == :terminal
        for i in findall(nodedegree(mol) .== 1)
            b = iterate(incidences(mol, i))[1]
            if bondorder_[b] == 2
                style[b] = 2
            end
        end
    elseif setting[:double_bond_style] == :chain
        for b in findall(.!isringbond(mol) .* (bondorder_ .== 2))
            style[b] = 2
        end
    end

    # Align double bonds alongside the ring
    for ring in sort(sssr(mol), by=length, rev=true)
        cw = isclockwise(toarray(coords, ring))
        cw === nothing && continue
        ordered = cw ? ring : reverse(ring)
        rr = vcat(ordered, ordered)
        for i in 1:length(ordered)
            e = findedgekey(mol, rr[i], rr[i + 1])
            (u, v) = getedge(mol, e)
            if bondorder_[e] == 2 && u != rr[i]
                style[e] = 1
            end
        end
    end

    return coords, style
end

coords2d(view::SubgraphView; kwargs...) = coords2d(view.graph; kwargs...)

"""
    chargesign(charge::Int) -> String

Get a charge sign.
"""
function chargesign(charge::Int)
    charge == 0 && return ""
    sign = charge > 0 ? "+" : "–" # en dash, not hyphen-minus
    num = abs(charge)
    return num > 1 ? string(num, sign) : sign
end

function atommarkup(symbol, charge, hydrogenconnected, direction,
                    substart, subend, supstart, supend)
    if hydrogenconnected == 1
        hs = "H"
    elseif hydrogenconnected > 1
        hs = string("H", substart, hydrogenconnected, subend)
    else
        hs = ""
    end
    chg = charge == 0 ? "" : string(supstart, chargesign(charge), supend)
    txt = direction == :left ? [chg, hs, symbol] : [symbol, hs, chg]
    return string(txt...)
end

atomhtml(
    symbol, charge, hydrogenconnected, direction
) = atommarkup(
    symbol, charge, hydrogenconnected, direction, "", "", "", "")

"""
    draw2d!(canvas::Canvas, mol::UndirectedGraph; kwargs...)

Draw molecular image to the canvas.
"""
function draw2d!(canvas::Canvas, mol::UndirectedGraph; kwargs...)
    # Canvas settings
    coords_, bondstyles_ = coords2d(mol; kwargs...)
    initcanvas!(canvas, coords_, boundary(mol, coords_))
    canvas.valid || return

    # Properties
    atomsymbol_ = atomsymbol(mol)
    atomcolor_ = atomcolor(mol; kwargs...)
    charge_ = charge(mol)
    isatomvisible_ = isatomvisible(mol; kwargs...)
    bondorder_ = bondorder(mol)
    implicithconnected_ = implicithconnected(mol)

    # Draw bonds
    for i in edgeset(mol)
        (u, v) = getedge(mol, i)
        if bondorder_[i] == 1 && bondstyles_[i] in (2, 7)
            # cordgen reverse bonds
            (t, s) = (u, v)
            bs = bondstyles_[i] - 1
        else
            (s, t) = (u, v)
            bs = bondstyles_[i]
        end
        setbond!(
            canvas, bondorder_[i], bs,
            Segment{Point2D}(canvas.coords, s, t),
            atomcolor_[u], atomcolor_[v],
            isatomvisible_[u], isatomvisible_[v]
        )
    end

    # Draw atoms
    for i in nodeset(mol)
        isatomvisible_[i] || continue
        pos = Point2D(canvas.coords, i)
        # Determine text direction
        if implicithconnected_[i] > 0
            cosnbrs = []
            hrzn = pos + (1.0, 0.0)
            for adj in adjacencies(mol, i)
                posnbr = Point2D(canvas.coords, adj)
                dist = Geometry.distance(pos, posnbr)
                if dist > 0
                    dp = dot(hrzn - pos, posnbr - pos)
                    push!(cosnbrs, dp / dist)
                end
            end
            if isempty(cosnbrs) || minimum(cosnbrs) > 0
                # [atom]< or isolated node(ex. H2O, HCl)
                setatomright!(
                    canvas, pos, atomsymbol_[i], atomcolor_[i],
                    implicithconnected_[i], charge_[i]
                )
                continue
            elseif maximum(cosnbrs) < 0
                # >[atom]
                setatomleft!(
                    canvas, pos, atomsymbol_[i], atomcolor_[i],
                    implicithconnected_[i], charge_[i]
                )
                continue
            end
        end
        # -[atom]- or no hydrogens
        setatomcenter!(
            canvas, pos, atomsymbol_[i], atomcolor_[i],
            implicithconnected_[i], charge_[i]
        )
    end
    return
end

function drawatomindex!(canvas::Canvas, mol::UndirectedGraph;
                        color=Color(0, 0, 0), bgcolor=Color(240, 240, 255))
    isatomvisible_ = isatomvisible(mol)
    for i in nodeset(mol)
        offset = isatomvisible_[i] ? (0.0, canvas.fontsize/2.0) : (0.0, 0.0)
        pos = Point2D(canvas.coords, i) + offset
        setatomnote!(canvas, pos, string(i), color, bgcolor)
    end
    return
end

function sethighlight!(
        canvas::Canvas, substr::UndirectedGraph; color=Color(253, 216, 53))
    isatomvisible_ = isatomvisible(substr)
    for i in edgeset(substr)
        (u, v) = getedge(substr, i)
        setbondhighlight!(canvas, Segment{Point2D}(canvas.coords, u, v), color)
    end
    for i in nodeset(substr)
        isatomvisible_[i] || continue
        pos = Point2D(canvas.coords, i)
        setatomhighlight!(canvas, pos, color)
    end
    return
end


text:
s = Dict("output" => Dict("branch_flows" => true), "conv_losses_mp" => true)
@testset "test ac polar pf" begin
    @testset "5-bus ac dc case" begin
        result = run_acdcpf("../test/data/case5_acdc.m", ACPPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 0; atol = 1e-2)

        @test isapprox(result["solution"]["gen"]["1"]["pg"], 1.3494; atol = 2e-3)
        @test isapprox(result["solution"]["gen"]["2"]["pg"], 0.40; atol = 2e-3)

        @test isapprox(result["solution"]["bus"]["1"]["vm"], 1.06; atol = 2e-3)
        @test isapprox(result["solution"]["bus"]["1"]["va"], 0.00000; atol = 2e-3)
        @test isapprox(result["solution"]["bus"]["2"]["vm"], 1.00; atol = 2e-3)
        @test isapprox(result["solution"]["bus"]["3"]["vm"], 0.995; atol = 2e-3)

        @test isapprox(result["solution"]["convdc"]["2"]["pgrid"], -0.1954; atol = 2e-3)
        @test isapprox(result["solution"]["convdc"]["3"]["pdc"], 0.36421; atol = 2e-3)
        @test isapprox(result["solution"]["busdc"]["1"]["vm"], 1.008; atol = 2e-3)

    end
    @testset "5-bus ac dc case with 2 seperate ac grids" begin
        result = run_acdcpf("../test/data/case5_2grids.m", ACPPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 0; atol = 1e-2)

        @test isapprox(result["solution"]["gen"]["1"]["pg"], 1.9326; atol = 2e-3)
        @test isapprox(result["solution"]["gen"]["4"]["pg"], 0.40; atol = 2e-3)

        @test isapprox(result["solution"]["bus"]["1"]["vm"], 1.06; atol = 2e-3)
        @test isapprox(result["solution"]["bus"]["1"]["va"], 0.00000; atol = 2e-3)
        @test isapprox(result["solution"]["bus"]["3"]["vm"], 0.987; atol = 2e-3)
        @test isapprox(result["solution"]["bus"]["7"]["va"], -0.0065; atol = 2e-3)
        @test isapprox(result["solution"]["bus"]["10"]["vm"], 0.972; atol = 2e-3)

        @test isapprox(result["solution"]["convdc"]["1"]["pgrid"], 0.6; atol = 2e-3)
        @test isapprox(result["solution"]["convdc"]["2"]["pdc"], 0.56872; atol = 2e-3)
        @test isapprox(result["solution"]["busdc"]["1"]["vm"], 1.015; atol = 2e-3)
    end
    # REMOVED for TRAVIS, otherwise case ok
    # @testset "24-bus rts ac dc case with three zones" begin
    #     result = run_acdcpf("../test/data/case24_3zones_acdc.m", ACPPowerModel, ipopt_solver; setting = s)
    #
    #     @test result["termination_status"] == LOCALLY_SOLVED
    #     @test isapprox(result["objective"], 0; atol = 1e-2)
    #
    #     @test isapprox(result["solution"]["gen"]["65"]["pg"], 1.419; atol = 2e-3)
    #     @test isapprox(result["solution"]["gen"]["65"]["qg"], -1.29648; atol = 2e-3)
    #
    #
    #     @test isapprox(result["solution"]["bus"]["101"]["vm"], 1.035; atol = 2e-3)
    #     @test isapprox(result["solution"]["bus"]["101"]["va"], -0.1389; atol = 2e-3)
    #     @test isapprox(result["solution"]["bus"]["205"]["vm"], 1.0318; atol = 2e-3)
    #     @test isapprox(result["solution"]["bus"]["301"]["vm"], 1.026266; atol = 2e-3)
    #
    #     @test isapprox(result["solution"]["convdc"]["2"]["pgrid"], -0.753; atol = 2e-3)
    #     @test isapprox(result["solution"]["convdc"]["3"]["pdc"], -1.37301; atol = 2e-3)
    #     @test isapprox(result["solution"]["busdc"]["5"]["vm"], 1.01731; atol = 2e-3)
    # end
end

@testset "test dc pf" begin
    @testset "5-bus ac dc case" begin
        result = run_acdcpf("../test/data/case5_acdc.m", DCPPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 0; atol = 1e-2)

        @test isapprox(result["solution"]["gen"]["1"]["pg"], 1.28313; atol = 2e-3)
        @test isapprox(result["solution"]["bus"]["5"]["va"], -0.09668; atol = 2e-3)
        @test isapprox(result["solution"]["bus"]["3"]["va"], -0.08918; atol = 2e-3)
        @test isapprox(result["solution"]["convdc"]["2"]["pgrid"], -0.21686; atol = 2e-3)
        @test isapprox(result["solution"]["convdc"]["3"]["pdc"], 0.36051; atol = 2e-3)
    end
    @testset "24-bus rts ac dc case with three zones" begin
        result = run_acdcpf("../test/data/case24_3zones_acdc.m", DCPPowerModel, ipopt_solver; setting = s)

        @test result["termination_status"] == LOCALLY_SOLVED
        @test isapprox(result["objective"], 0; atol = 1e-2)
        @test isapprox(result["solution"]["gen"]["65"]["pg"], 1.419; atol = 2e-3)

        @test isapprox(result["solution"]["bus"]["119"]["va"], 0.17208; atol = 2e-3)
        @test isapprox(result["solution"]["bus"]["224"]["va"], 0.07803; atol = 2e-3)

        @test isapprox(result["solution"]["convdc"]["2"]["pgrid"], -0.753; atol = 2e-3)
        @test isapprox(result["solution"]["convdc"]["3"]["pdc"], -1.39167; atol = 2e-3)
    end
end


text:
using LinearAlgebra
using PolynomialBasis
using ImplicitDomainQuadrature
using Revise
using CutCell
include("../../../test/useful_routines.jl")

function bulk_modulus(l, m)
    return l + 2m / 3
end

function lame_lambda(k, m)
    return k - 2m / 3
end

function solve_for_displacement(
    basis,
    cellquads,
    interfacequads,
    facequads,
    cutmesh,
    levelset,
    levelsetcoeffs,
    stiffness,
    theta0,
    penalty,
)

    bilinearforms = CutCell.BilinearForms(basis, cellquads, stiffness, cutmesh)
    interfacecondition = CutCell.incoherent_interface_condition(
        basis,
        interfacequads,
        stiffness,
        cutmesh,
        penalty,
    )

    sysmatrix = CutCell.SystemMatrix()
    sysrhs = CutCell.SystemRHS()

    CutCell.assemble_bilinear_form!(sysmatrix, bilinearforms, cutmesh)
    CutCell.assemble_interface_condition!(sysmatrix, interfacecondition, cutmesh)
    CutCell.assemble_bulk_transformation_linear_form!(
        sysrhs,
        transfstress,
        basis,
        cellquads,
        cutmesh,
    )
    CutCell.assemble_incoherent_interface_transformation_rhs!(
        sysrhs,
        transfstress,
        basis,
        interfacequads,
        cutmesh,
    )

    matrix = CutCell.make_sparse(sysmatrix, cutmesh)
    rhs = CutCell.rhs(sysrhs, cutmesh)

    topleftnodeid = CutCell.nodes_per_mesh_side(cutmesh.mesh)[2]
    CutCell.apply_dirichlet_bc!(matrix, rhs, [1, topleftnodeid], 1, 0.0, 2)
    CutCell.apply_dirichlet_bc!(matrix, rhs, [1], 2, 0.0, 2)

    nodaldisplacement = matrix \ rhs

    return nodaldisplacement
end


text:
#!/bin/env julia
using InteractiveUtils

const PLACEHOLDER = "|"

function main()
    clip = clipboard()
    clip_ = String(clip)

    if !isempty(ARGS)

        d = []

        for i in readlines(ARGS[1])
            s = split(i)
            push!(d, [String(s[1]),String(s[2])])
        end

        sort!(d, rev=true, lt=(x,y)->isless(length(x[1]),length(y[1])))

        i = 1
        while i <= length(d)
            pos = findfirst(d[i][1], clip_)
            if pos == nothing
                i += 1
            else
                clip_ = replace(clip_, d[i][1] => repeat(PLACEHOLDER, sizeof(d[i][2])), count=1)
                c = ""
                if pos[1] > 1
                    c *= clip[1:prevind(clip, pos[1])]
                end
                c *= d[i][2]
                replace_end = nextind(clip, last(pos))
                if replace_end <= lastindex(clip)
                    c *= clip[replace_end:lastindex(clip)]
                end
                clip = c
            end
        end

    end

    clipboard(clip)

end

main()

max_stars_repo_name: gstavrinos/rclip

max_issues_repo_name: gstavrinos/rclip

max_forks_repo_name: gstavrinos/rclip


text:
using .Curves: ECPoint, gx, gy
using .Specs: PRG

abstract type Group end

Base.broadcasted(f::Function, x::Group, y::AbstractVector{<:Integer}) = f.((x for i in 1:length(y)), y)

Base.broadcasted(::typeof(*), x::G, y::Vector{G}) where G <: Group = (x for i in 1:length(y)) .* y
Base.broadcasted(::typeof(*), x::Vector{G}, y::G) where G <: Group =  x .* (y for i in 1:length(x))

order(x::G) where G <: Group = order(G)

import Base./
/(x::G, y::G) where G <: Group = x * inv(y)

name(x::G) where G <: Group = name(G)

Base.convert(::Type{Vector{G}}, x::Vector) where G <: Group = G[ G <| i for i in x]

Base.rand(prg::PRG, ::Type{G}, N::Integer; nr::Integer = 0) where G <: Group = Vector{G} <| rand(prg, spec(G), N; nr)

struct ECGroup{P<:ECPoint} <: Group
    x::P
end

Base.convert(::Type{ECGroup{P}}, x) where P <: ECPoint = ECGroup{P}(P <| x)
Base.convert(::Type{G}, x::G) where G <: ECGroup = x

Base.:*(x::G, y::G) where G <: ECGroup = G(x.x + y.x)

function Base.:^(x::G, n::Integer) where G <: ECGroup
    n_mod = mod(n, order(G))

    @assert n_mod != 0 "A bad exponent"

    return G(n_mod * x.x)
end

Base.inv(g::G) where G <: ECGroup = g^(order(G) - 1)

order(::Type{ECGroup{P}}) where P = order(P)

Base.isvalid(x::ECGroup) = isvalid(x.x)

Base.:(==)(x::G, y::G) where G <: ECGroup = x.x == y.x

modulus(::Type{ECGroup{P}}) where P <: ECPoint = modulus(P)

name(::Type{ECGroup{P}}) where P <: ECPoint = name(P)

Base.isless(x::G, y::G) where G <: ECGroup = isless(x.x, y.x)

Curves.gx(g::ECGroup) = gx(g.x)
Curves.gy(g::ECGroup) = gy(g.x)

function Base.show(io::IO, g::G) where G <: ECGroup
    show(io, G)
    print(io, " <| (")
    show(io, gx(g))
    print(io, ", ")
    show(io, gy(g))
    print(io, ")")
end

struct PGroup{S} <: Group
    g::BigInt

    PGroup(p::Integer, q::Integer; name=nothing) = PGroup{static(; p, q, name)}

    function PGroup{S}(x::BigInt) where S
        #@assert 1 < x < S.p "Not in range"
        # Relaxing a little so that group products could happen propeerly
        # But is a presence of 1 in multiplicative products a vulnerability?
        @assert 0 < x < S.p "Not in range"
        new{S}(x)
    end

    PGroup{S}(x::Integer) where S = PGroup{S}(BigInt(x))

end

modulus(::Type{PGroup{S}}) where S = BigInt(S.p)
modulus(::G) where G <: PGroup = modulus(G)

order(::Type{PGroup{S}}) where S = BigInt(S.q)

name(::Type{PGroup}) = nothing

name(::Type{PGroup{S}}) where S =  isnothing(S.name) ? nothing : convert(Symbol, S.name)

Base.show(io::IO, g::PGroup) = print(io, value(g))

Base.show(io::IO, ::Type{PGroup}) = print(io, "PGroup")

#dublicate at Fields.
function trimnumber(x::String)
    if length(x) < 30
        return x
    else
        return x[1:10] * "..." * x[end-10:end]
    end
end

trimnumber(x::Integer) = trimnumber(string(x))

groupstr(m) = "𝐙/($(trimnumber(m)))"

function Base.show(io::IO, ::Type{G}) where G <: PGroup

    if name(G) == nothing
        print(io, groupstr(modulus(G)))
    else
        print(io, name(G))
    end
end

function Base.display(::Type{G}) where G <: PGroup
    show(G)
    print(" (order = $(order(G)))")
end

value(g::PGroup) = g.g

Base.convert(::Type{P}, x::Integer) where P <: PGroup = P(BigInt(x))

Base.isvalid(g::G) where G <: PGroup = value(g) != 1 && powermod(value(g), order(G), modulus(G)) == 1

import Base.*
*(x::G, y::G) where G <: PGroup = G(mod(value(x) * value(y), modulus(G)))

import Base.^
function ^(x::G, n::Integer) where G <: PGroup

    n_mod = mod(n, order(G))
    @assert n_mod != 0 "A bad exponent"
    @assert value(x) != 1 "A value 1 is not part of prime group"

    return G(powermod(value(x), n_mod, modulus(G)))
end

import Base.inv
inv(x::G) where G <: PGroup = G(modinv(value(x), modulus(G)))

Base.:(==)(x::G, y::G) where G <: PGroup = x.g == y.g

Base.isless(x::G, y::G) where G <: PGroup = value(x) < value(y)

#Base.convert(::Type{G}, x::Integer) where G <: PGroup = G(x)

# function Base.prod(x::Vector{G}) where G <: PGroup

#     p = modulus(G)1
#     s = value(x[1])

#     for i in x[2:end]
#         s = mod(s * value(i), p)
#     end

#     return G(s)
# end


text:
function EigenDecomposition(arg0::RealMatrix)
    return EigenDecomposition((RealMatrix,), arg0)
end

function EigenDecomposition(arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return EigenDecomposition((Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function get_d(obj::EigenDecomposition)
    return jcall(obj, "getD", RealMatrix, ())
end

function get_determinant(obj::EigenDecomposition)
    return jcall(obj, "getDeterminant", jdouble, ())
end

function get_eigenvector(obj::EigenDecomposition, arg0::jint)
    return jcall(obj, "getEigenvector", RealVector, (jint,), arg0)
end

function get_imag_eigenvalue(obj::EigenDecomposition, arg0::jint)
    return jcall(obj, "getImagEigenvalue", jdouble, (jint,), arg0)
end

function get_imag_eigenvalues(obj::EigenDecomposition)
    return jcall(obj, "getImagEigenvalues", Vector{jdouble}, ())
end

function get_real_eigenvalue(obj::EigenDecomposition, arg0::jint)
    return jcall(obj, "getRealEigenvalue", jdouble, (jint,), arg0)
end

function get_real_eigenvalues(obj::EigenDecomposition)
    return jcall(obj, "getRealEigenvalues", Vector{jdouble}, ())
end

function get_solver(obj::EigenDecomposition)
    return jcall(obj, "getSolver", DecompositionSolver, ())
end

function get_square_root(obj::EigenDecomposition)
    return jcall(obj, "getSquareRoot", RealMatrix, ())
end

function get_v(obj::EigenDecomposition)
    return jcall(obj, "getV", RealMatrix, ())
end

function get_vt(obj::EigenDecomposition)
    return jcall(obj, "getVT", RealMatrix, ())
end

function has_complex_eigenvalues(obj::EigenDecomposition)
    return jcall(obj, "hasComplexEigenvalues", jboolean, ())
end