text:
using .SnakeGames

function SnakeGameEnv(; action_style = MINIMAL_ACTION_SET, kw...)
    game = SnakeGame(; kw...)
    n_snakes = length(game.snakes)
    num_agent_style = n_snakes == 1 ? SINGLE_AGENT : MultiAgent{n_snakes}()
    SnakeGameEnv{action_style,num_agent_style,typeof(game)}(
        game,
        map(length, game.snakes),
        Vector{CartesianIndex{2}}(undef, length(game.snakes)),
        false,
    )
end

RLBase.ActionStyle(env::SnakeGameEnv{A}) where {A} = A
RLBase.NumAgentStyle(env::SnakeGameEnv{<:Any,N}) where {N} = N
RLBase.DynamicStyle(env::SnakeGameEnv{<:Any,SINGLE_AGENT}) = SEQUENTIAL
RLBase.DynamicStyle(env::SnakeGameEnv{<:Any,<:MultiAgent}) = SIMULTANEOUS

const SNAKE_GAME_ACTIONS = (
    CartesianIndex(-1, 0),
    CartesianIndex(1, 0),
    CartesianIndex(0, 1),
    CartesianIndex(0, -1),
)

function (env::SnakeGameEnv{A})(actions::Vector{CartesianIndex{2}}) where {A}
    if A === MINIMAL_ACTION_SET
        # avoid turn back
        actions = [
            a_new == -a_old ? a_old : a_new for
            (a_new, a_old) in zip(actions, env.latest_actions)
        ]
    end

    env.latest_actions .= actions
    map!(length, env.latest_snakes_length, env.game.snakes)
    env.is_terminated = !env.game(actions)
end

(env::SnakeGameEnv)(action::Int) = env([SNAKE_GAME_ACTIONS[action]])
(env::SnakeGameEnv)(actions::Vector{Int}) = env(map(a -> SNAKE_GAME_ACTIONS[a], actions))

RLBase.action_space(env::SnakeGameEnv) = 1:4
RLBase.state(env::SnakeGameEnv) = env.game.board
RLBase.state_space(env::SnakeGameEnv) = Space(fill(false..true, size(env.game.board)))
RLBase.reward(env::SnakeGameEnv{<:Any,SINGLE_AGENT}) =
    length(env.game.snakes[]) - env.latest_snakes_length[]
RLBase.reward(env::SnakeGameEnv) = length.(env.game.snakes) .- env.latest_snakes_length
RLBase.is_terminated(env::SnakeGameEnv) = env.is_terminated

RLBase.legal_action_space(env::SnakeGameEnv{FULL_ACTION_SET,SINGLE_AGENT}) =
    findall(!=(-env.latest_actions[]), SNAKE_GAME_ACTIONS)
RLBase.legal_action_space(env::SnakeGameEnv{FULL_ACTION_SET}) =
    [findall(!=(-a), SNAKE_GAME_ACTIONS) for a in env.latest_actions]

RLBase.legal_action_space_mask(env::SnakeGameEnv{FULL_ACTION_SET,SINGLE_AGENT}) =
    [a != -env.latest_actions[] for a in SNAKE_GAME_ACTIONS]
RLBase.legal_action_space_mask(env::SnakeGameEnv{FULL_ACTION_SET}) =
    [[x != -a for x in SNAKE_GAME_ACTIONS] for a in env.latest_actions]

function RLBase.reset!(env::SnakeGameEnv)
    SnakeGames.reset!(env.game)
    env.is_terminated = false
    fill!(env.latest_actions, CartesianIndex(0, 0))
    map!(length, env.latest_snakes_length, env.game.snakes)
end

Base.display(io::IO, m::MIME, env::SnakeGameEnv) = display(io, m, env.game)


text:
using DataFrames
import FileIO:save, load

# Top-level runner functions for Experimental Prediction (section 2.4)
#
# Sgm stands for Bayes Factor Surprise
# Ssh stands for Shannon Surprise

# Experimental Prediction 1
function experiment_mania_First(;
                    nsteps = [500, 1000, 2000],
                    y_hat_centers = [1.],
                    var_centers = [0.5],
                    δ_centers = collect(0.:0.1:1.),
                    Δy_hats = [0.25],
                    Δvars = [1.],
                    Δδs = [0.1],
                    seeds = collect(1:20),
                    sigma = 0.5,
                    changeprobability = 0.1,
                    foldername = "experiment",
                    filename = "experiment_")

    learnerlist = [
                    (name = "pf20", type = GParticleFilter, Npar = 20, m = 0.),
                    (name = "gnassar12", type = GNas12, Npar = 0, m = 0.)
                    ]
    for i in learnerlist
        savepath = makehomesavepath(foldername)
        mkdir(savepath)
        df = runExperimentsSingleLearner_First(nsteps = nsteps,
                        y_hat_centers = y_hat_centers,
                        var_centers = var_centers,
                        δ_centers = δ_centers,
                        Δy_hats = Δy_hats,
                        Δvars = Δvars,
                        Δδs = Δδs,
                        seeds = seeds,
                        sigma = sigma, changeprobability = changeprobability,
                        learnertype = i.type,
                        Npar = i.Npar, m = i.m)

        save(joinpath(savepath, filename * i.name * ".jld2"),
                                Dict("df" => df))
    end
end

# Experimental Prediction 2
function experiment_mania_Second(;
                    nsteps = [500, 1000, 2000],
                    p_centers = collect(0.025:0.025:0.4),
                    Δp = [0.0125],
                    seeds = collect(1:20),
                    sigma = 0.5,
                    changeprobability = 0.1,
                    foldername = "experiment_second",
                    filename = "experiment_second_")

    learnerlist = [
                    (name = "pf20", type = GParticleFilter, Npar = 20, m = 0.),
                    (name = "gnassar12", type = GNas12, Npar = 0, m = 0.)
                    ]
    for i in learnerlist
        savepath = makehomesavepath(foldername)
        mkdir(savepath)
        df = runExperimentsSingleLearner_Second(nsteps = nsteps,
                        p_centers = p_centers,
                        Δp = Δp,
                        seeds = seeds,
                        sigma = sigma, changeprobability = changeprobability,
                        learnertype = i.type,
                        Npar = i.Npar, m = i.m)

        save(joinpath(savepath, filename * i.name * ".jld2"),
                                Dict("df" => df))
    end
end

ext: jl


text:
using QContinuousMeasurement
using Test

@testset "Pretty printing" begin
    @testset "LocalDephasingModel" begin
        Nj = 2
        modelparams = ModelParameters(Nj=Nj)

        model = LocalDephasingModel(modelparams)
        initial_state = blockdiag_css(Nj)

        println(model)
        println(initial_state)
    end
end


text:
# ---
# title: 1627. Graph Connectivity With Threshold
# id: problem1627
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Math, Union Find
# link: <https://leetcode.com/problems/graph-connectivity-with-threshold/description/>
# hidden: true
# ---
#
# We have `n` cities labeled from `1` to `n`. Two different cities with labels
# `x` and `y` are directly connected by a bidirectional road if and only if `x`
# and `y` share a common divisor **strictly greater** than some `threshold`.
# More formally, cities with labels `x` and `y` have a road between them if
# there exists an integer `z` such that all of the following are true:
#
#   * `x % z == 0`,
#   * `y % z == 0`, and
#   * `z > threshold`.
#
# Given the two integers, `n` and `threshold`, and an array of `queries`, you
# must determine for each `queries[i] = [ai, bi]` if cities `ai` and `bi` are
# connected (i.e. there is some path between them).
#
# Return _an array_`answer` _, where_`answer.length == queries.length`
# _and_`answer[i]` _is_`true` _if for the_`ith` _query, there is a path
# between_`ai` _and_`bi` _, or_`answer[i]` _is_`false` _if there is no path._
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/09/ex1.jpg)
#
#
#
#     Input: n = 6, threshold = 2, queries = [[1,4],[2,5],[3,6]]
#     Output: [false,false,true]
#     Explanation: The divisors for each number:
#     1:   1
#     2:   1, 2
#     3:   1, _3_
#     4:   1, 2, _4_
#     5:   1, _5_
#     6:   1, 2, _3_ , _6_
#     Using the underlined divisors above the threshold, only cities 3 and 6 share a common divisor, so they are the
#     only ones directly connected. The result of each query:
#     [1,4]   1 is not connected to 4
#     [2,5]   2 is not connected to 5
#     [3,6]   3 is connected to 6 through path 3--6
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/10/tmp.jpg)
#
#
#
#     Input: n = 6, threshold = 0, queries = [[4,5],[3,4],[3,2],[2,6],[1,3]]
#     Output: [true,true,true,true,true]
#     Explanation: The divisors for each number are the same as the previous example. However, since the threshold is 0,
#     all divisors can be used. Since all numbers share 1 as a divisor, all cities are connected.
#
#
# **Example 3:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/17/ex3.jpg)
#
#
#
#     Input: n = 5, threshold = 1, queries = [[4,5],[4,5],[3,2],[2,3],[3,4]]
#     Output: [false,false,false,false,false]
#     Explanation: Only cities 2 and 4 share a common divisor 2 which is strictly greater than the threshold 1, so they are the only ones directly connected.
#     Please notice that there can be multiple queries for the same pair of nodes [x, y], and that the query [x, y] is equivalent to the query [y, x].
#
#
#
#
# **Constraints:**
#
#   * `2 <= n <= 104`
#   * `0 <= threshold <= n`
#   * `1 <= queries.length <= 105`
#   * `queries[i].length == 2`
#   * `1 <= ai, bi <= cities`
#   * `ai != bi`
#
#
## @lc code=start
using LeetCode

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


text:
ZygoteRules.@adjoint function ODESolution(u,args...)
  function ODESolutionAdjoint(ȳ)
    (ȳ,ntuple(_->nothing, length(args))...)
  end
  ODESolution(u,args...), ODESolutionAdjoint
end

ZygoteRules.@adjoint function ODESolution{T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11
                }(u,args...) where {T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11}

                function ODESolutionAdjoint(ȳ)
                  (ȳ,ntuple(_->nothing, length(args))...)
                end

                ODESolution{T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11}(u,args...),ODESolutionAdjoint
end

ZygoteRules.@adjoint function getindex(sol::DESolution, i)
  function DESolution_getindex_adjoint(Δ)
    Δ′ = Union{Nothing, eltype(sol.u)}[nothing for x in sol.u]
    Δ′[i] = Δ
    (Δ′,nothing)
  end
  sol[i],DESolution_getindex_adjoint
end

ZygoteRules.@adjoint function getindex(sol::DESolution, i, j...)
  function DESolution_getindex_adjoint(Δ)
    Δ′ = zero(sol)
    Δ′[i,j...] = Δ
    (Δ′, map(_ -> nothing, i)...)
  end
  sol[i,j...],DESolution_getindex_adjoint
end

ZygoteRules.@adjoint function (f::ODEFunction)(u,p,t)
  if f.vjp === nothing
    ZygoteRules._pullback(f.f,u,p,t)
  else
    f.vjp(u,p,t)
  end
end

ZygoteRules.@adjoint! function (f::ODEFunction)(du,u,p,t)
  if f.vjp === nothing
    ZygoteRules._pullback(f.f,du,u,p,t)
  else
    f.vjp(du,u,p,t)
  end
end

#=
ChainRulesCore.frule(f::ODEFunction,u,p,t)
  if f.jvp === nothing
    ChainRulesCore.frule(f.f,u,p,t)
  else
    function ode_jvp(f,du,dp,dt)
      f.jvp_u(du,u,p,t) + f.jvp_p(dp,u,p,t) + f.tgrad(u,p,t)*dt
    end
    f.f(u,p,t),ode_jvp
end
=#

function ChainRulesCore.rrule(f::ODEFunction,u,p,t)
  if f.vjp === nothing
    ChainRulesCore.rrule(f.f,u,p,t)
  else
    f.vjp(u,p,t)
  end
end

ZygoteRules.@adjoint numargs(f) = (numargs(f),df->(nothing,))
ChainRulesCore.rrule(::typeof(numargs),f) = (numargs(f),df->(nothing,))


text:
"""
Add values to extras
"""
function addextras(data, lassonumvars, betas, lambda)
    data.extras[ModelSelection.generate_extra_key(PRELIMINARYSELECTION_EXTRAKEY, data.extras)] = Dict(
        :preliminaryselection => :lasso,
        :lassonumvars => lassonumvars,
        :lassobetas => betas,
        :lassolambda => lambda,
        :nobs => data.nobs
    )
    return data
end


text:
# # Number lists (wikitext references)
#
# Wonder what parsing is useful for?
# This example presents how to read a convenient number list format, e.g.  `[1-3,9]`, of the Wikipedia wikitext format.
# This makes a simple example for what parsing into Julia types means, too.
#
# To reflect on the amazing Julia type system, the example shows
# - different ways to represent such number `Iterator`s in Julia,
1:3 == [1,2,3]
# - how you parse into any such representation
using CombinedParsers
# - inter-operation with brief regular expression (PCRE) syntax
using CombinedParsers.Regexp
dash = re" *- *"
# - inter-operation with
import TextParse
# ## Number ranges
@syntax int_range = Sequence(
    Numeric(Int), # 1
    dash,    # 2
    Numeric(Int)  # 3
) do v
    v[1]:v[3]
end;

# to match cases like "8  - 11".
# Julias number range format is `1:3`.
# The string macro `@int_range_str` is defined in `@syntax`.

int_range"1-3"
int_range"8-11"

# Julia Base.collect can be used to convert
@syntax int_vector = map(collect, int_range);
int_vector"8  - 11"

# ## Numbers
# Without `@syntax`, you can parse
int = map(x -> [x], Numeric(Int));

parse(int,"19")
int("19")

# The tree displays how a number is read and transformed to a Vector with length 1.
int

# ## Joining numbers and ranges
@syntax numbers = map(join(
    Repeat(Either(int_vector, int)),
    re" *, *"
)) do v
    vcat(v...)::Vector{Int}
end;

numbers"1-3,9"

# Prepend another parser by
(re"no *"*numbers)[end]("no 2-4,19")

# ## Inclusion in a wikitext parser
# Long and complicated texts like the Wikipedia can be parsed with `CombinedParsers.jl`.
# The parsers are less pain to write and execute at speeds comparably to PCRE implemented in C, the regular expressions industry standard.
# `CombinedParsers.jl` can inter-operate with Julia packages `TextParse.jl`.

@syntax wiki_references = Sequence(2,"[",numbers,"]");

wiki_references"[1, 7-9, 2]"

# The tree displays how a bracketed comma separated sequence of numbers and number ranges is read and transformed to a Vector.
wiki_references

# ## PCRE papercuts when parsing number sequences
# The same parser as a regular expression will be tedious to understand and write (though writing the regular expression `re" *- *"` is clear).
# PCRE matching does recognize the match but makes not all required parsing parts accessible (7 is not captured).
re = "\\[([[:digit:]]+ *- *[[:digit:]]+|[[:digit:]]+)(?: *, *([[:digit:]]+ *- *[[:digit:]]+|[[:digit:]]+))*\\]"
match(Regex("^"*re*"\$"), "[1-3,7,9]")
### TODO: regex_string(wiki_references)

# To make the parsing work with regular expressions you would choose a stepwise strategy, handling `[` and `]` and stepping though `,` separated parts.
# Parsing was pain.


text:
module Day1

export day1

function day1(::Val{:a}, filename, n)
    nums = map(s -> parse(Int,s), readlines(filename)) |> Set{Int}

    for k in nums
        if n-k in nums
            return k * (n-k)
        end
    end
    nothing
end

function day1(::Val{:b}, filename, n)
    nums = map(s -> parse(Int,s), readlines(filename)) |> Set{Int}
    for i in nums
        for j in nums
            m = n-i-j
            if m in nums
                return i * j * m
            end
        end
    end
    nothing
end

end


text:
"""
    RelativeEntropyBridge{T}

The `RelativeEntropyCone` is representable with exponential cone and LP constraints, since
``u \\ge \\sum_{i=1}^n w_i \\log (\\frac{w_i}{v_i})`` if and only if there exists a vector
``y`` such that ``u \\ge \\sum_i y_i`` and ``y_i \\ge w_i \\log (\\frac{w_i}{v_i})`` or
equivalently ``v_i \\ge w_i \\exp (\\frac{-y_i}{w_i})`` or equivalently
``(-y_i, w_i, v_i) \\in ExponentialCone``, for all ``i``.
"""
struct RelativeEntropyBridge{T, F, G, H} <: AbstractBridge
    y::Vector{MOI.VariableIndex}
    ge_index::CI{F, MOI.GreaterThan{T}}
    exp_indices::Vector{CI{G, MOI.ExponentialCone}}
end
function bridge_constraint(::Type{RelativeEntropyBridge{T, F, G, H}}, model::MOI.ModelLike, f::H, s::MOI.RelativeEntropyCone) where {T, F, G, H}
    f_scalars = MOIU.eachscalar(f)
    d = MOI.dimension(s)
    v_dim = div(d - 1, 2)
    y = MOI.add_variables(model, v_dim)
    ge_index = MOIU.normalize_and_add_constraint(model, MOIU.operate(-, T, f_scalars[1], MOIU.operate(sum, T, y)), MOI.GreaterThan(zero(T)), allow_modify_function=true)
    w_start = 1 + v_dim
    exp_funcs = [MOIU.operate(vcat, T, MOIU.operate(-, T, MOI.SingleVariable(y[i])), f_scalars[w_start + i], f_scalars[1 + i]) for i in 1:v_dim]
    exp_indices = [MOI.add_constraint(model, exp_func_i, MOI.ExponentialCone()) for exp_func_i in exp_funcs]
    return RelativeEntropyBridge{T, F, G, H}(y, ge_index, exp_indices)
end

MOI.supports_constraint(::Type{RelativeEntropyBridge{T}}, ::Type{<:MOI.AbstractVectorFunction}, ::Type{MOI.RelativeEntropyCone}) where T = true
MOIB.added_constrained_variable_types(::Type{<:RelativeEntropyBridge}) = Tuple{DataType}[]
MOIB.added_constraint_types(::Type{RelativeEntropyBridge{T, F, G, H}}) where {T, F, G, H} = [(F, MOI.GreaterThan{T}), (G, MOI.ExponentialCone)]
function concrete_bridge_type(::Type{<:RelativeEntropyBridge{T}}, H::Type{<:MOI.AbstractVectorFunction}, ::Type{MOI.RelativeEntropyCone}) where T
    S = MOIU.scalar_type(H)
    F = MOIU.promote_operation(-, T, S, S)
    G = MOIU.promote_operation(vcat, T, MOIU.promote_operation(-, T, MOI.SingleVariable), S)
    return RelativeEntropyBridge{T, F, G, H}
end

# Attributes, Bridge acting as a model
MOI.get(bridge::RelativeEntropyBridge, ::MOI.NumberOfVariables) = length(bridge.y)
MOI.get(bridge::RelativeEntropyBridge, ::MOI.ListOfVariableIndices) = bridge.y
MOI.get(bridge::RelativeEntropyBridge{T, F}, ::MOI.NumberOfConstraints{F, MOI.GreaterThan{T}}) where {T, F} = 1
MOI.get(bridge::RelativeEntropyBridge{T, F, G}, ::MOI.NumberOfConstraints{G, MOI.ExponentialCone}) where {T, F, G} = length(bridge.y)
MOI.get(bridge::RelativeEntropyBridge{T, F}, ::MOI.ListOfConstraintIndices{F, MOI.GreaterThan{T}}) where {T, F} = [bridge.ge_index]
MOI.get(bridge::RelativeEntropyBridge{T, F, G}, ::MOI.ListOfConstraintIndices{G, MOI.ExponentialCone}) where {T, F, G} = bridge.exp_indices

# References
function MOI.delete(model::MOI.ModelLike, bridge::RelativeEntropyBridge)
    for exp_index_i in bridge.exp_indices
        MOI.delete(model, exp_index_i)
    end
    MOI.delete(model, bridge.ge_index)
    MOI.delete(model, bridge.y)
end

# Attributes, Bridge acting as a constraint
function MOI.get(model::MOI.ModelLike, ::MOI.ConstraintFunction, bridge::RelativeEntropyBridge{T, F, G, H}) where {T, F, G, H}
    func = MOIU.zero_with_output_dimension(G, 1 + 2 * length(bridge.y))
    MOIU.operate_output_index!(+, T, 1, func, MOI.get(model, MOI.ConstraintFunction(), bridge.ge_index))
    w_start = 1 + length(bridge.y)
    for i in eachindex(bridge.y)
        exp_func_i = MOIU.eachscalar(MOI.get(model, MOI.ConstraintFunction(), bridge.exp_indices[i]))
        MOIU.operate_output_index!(-, T, 1, func, exp_func_i[1])
        MOIU.operate_output_index!(+, T, 1 + i, func, exp_func_i[3])
        MOIU.operate_output_index!(+, T, w_start + i, func, exp_func_i[2])
    end
    return MOIU.convert_approx(H, MOIU.remove_variable(func, bridge.y))
end
MOI.get(model::MOI.ModelLike, ::MOI.ConstraintSet, bridge::RelativeEntropyBridge) = MOI.RelativeEntropyCone(1 + 2 * length(bridge.y))
MOI.supports(::MOI.ModelLike, ::Union{MOI.ConstraintPrimalStart, MOI.ConstraintDualStart}, ::Type{<:RelativeEntropyBridge}) = true
function MOI.get(model::MOI.ModelLike, attr::Union{MOI.ConstraintPrimal, MOI.ConstraintPrimalStart}, bridge::RelativeEntropyBridge{T}) where T
    primal = zeros(T, 1 + 2 * length(bridge.y))
    primal[1] = MOI.get(model, attr, bridge.ge_index)
    w_start = 1 + length(bridge.y)
    for i in eachindex(bridge.y)
        exp_primal_i = MOI.get(model, attr, bridge.exp_indices[i])
        primal[1] -= exp_primal_i[1]
        primal[1 + i] = exp_primal_i[3]
        primal[w_start + i] = exp_primal_i[2]
    end
    return primal
end
function MOI.set(model::MOI.ModelLike, attr::MOI.ConstraintPrimalStart, bridge::RelativeEntropyBridge{T}, value) where T
    v_dim = length(bridge.y)
    v_value = value[2:(v_dim + 1)]
    w_value = value[(v_dim + 2):end]
    y_value = [w_i * log(w_i / v_i) for (v_i, w_i) in zip(v_value, w_value)]
    MOI.set(model, attr, bridge.ge_index, value[1] - reduce(+, y_value, init=zero(T)))
    for i in 1:v_dim
        MOI.set(model, attr, bridge.exp_indices[i], [-y_value[i], w_value[i], v_value[i]])
        MOI.set(model, MOI.VariablePrimalStart(), bridge.y[i], y_value[i])
    end
    return
end
# Given a is dual on u - sum(y) >= 0 and (b_i, c_i, d_i) is dual on (-y_i, w_i, v_i)
# in ExponentialCone, the dual on (u, v, w) in RelativeEntropyCone is (a, d_i, c_i).
function MOI.get(model::MOI.ModelLike, attr::Union{MOI.ConstraintDual, MOI.ConstraintDualStart}, bridge::RelativeEntropyBridge{T}) where {T}
    dual = zeros(T, 1 + 2 * length(bridge.y))
    dual[1] = MOI.get(model, attr, bridge.ge_index)[1]
    w_start = 1 + length(bridge.y)
    for i in eachindex(bridge.y)
        exp_dual_i = MOI.get(model, attr, bridge.exp_indices[i])
        dual[1 + i] = exp_dual_i[3]
        dual[w_start + i] = exp_dual_i[2]
    end
    return dual
end
# Given constraint dual start of (u, v, w), constraint dual on GreaterThan is u
# and on exponential cone constraint i is (r_i, w_i, v_i), but since y_i is free,
# its dual is 0, so we have -r_i + u == 0 hence r_i = u.
# Note: alternatively, we could use the Lambert W function to calculate
# r_i = exp(W(-w_i / (-v_i * e))) * (-v_i * e), but this is more complicated.
function MOI.set(model::MOI.ModelLike, ::MOI.ConstraintDualStart, bridge::RelativeEntropyBridge, value)
    u = value[1]
    MOI.set(model, MOI.ConstraintDualStart(), bridge.ge_index, u)
    w_start = 1 + length(bridge.y)
    for i in eachindex(bridge.y)
        MOI.set(model, MOI.ConstraintDualStart(), bridge.exp_indices[i], [u, value[w_start + i], value[1 + i]])
    end
    return
end


text:
function AbstractPlotting.plot!(
        lscene::LScene, P::AbstractPlotting.PlotFunc,
        attributes::AbstractPlotting.Attributes, args...;
        kw_attributes...)

    plot = AbstractPlotting.plot!(lscene.scene, P, attributes, args...; kw_attributes...)[end]

    plot
end

protrusionnode(ls::LScene) = ls.layoutnodes.protrusions
computedsizenode(ls::LScene) = ls.layoutnodes.computedsize

function LScene(parent::Scene; bbox = nothing, scenekw = NamedTuple(), kwargs...)

    attrs = merge!(Attributes(kwargs), default_attributes(LScene, parent))

    sizeattrs = sizenode!(attrs.width, attrs.height)
    alignment = lift(tuple, attrs.halign, attrs.valign)

    suggestedbbox = create_suggested_bboxnode(bbox)

    autosizenode = Node{NTuple{2, Optional{Float32}}}((nothing, nothing))

    computedsize = computedsizenode!(sizeattrs, autosizenode)

    finalbbox = alignedbboxnode!(suggestedbbox, computedsize, alignment, sizeattrs, autosizenode)

    protrusions = Node(RectSides{Float32}(0, 0, 0, 0))

    scene = Scene(parent, lift(IRect2D, finalbbox); scenekw...)

    layoutnodes = LayoutNodes{LScene, GridLayout}(suggestedbbox, protrusions, computedsize, autosizenode, finalbbox, nothing)

    LScene(scene, attrs, layoutnodes)
end


text:
using Bigleaf

using DataFrames, Pipe, Missings
using Dates, TimeZones
using Statistics

using Latexify
using DataDeps, Suppressor
using RData
import CodecBzip2, CodecXz
#@suppress_err # error in github-actions: GitHubActionsLogger has no field stream
register(DataDep(
    "DE_Tha_Jun_2014.rda",
    "downloading exampple dataset DE_Tha_Jun_2014 from bitbucket.org/juergenknauer/bigleaf",
    "https://bitbucket.org/juergenknauer/bigleaf/raw/0ebe11626b4409305951e8add9f6436703c82584/data/DE_Tha_Jun_2014.rda",
    "395f02e1a1a2d175ac7499c200d9d48b1cb58ff4755dfd2d7fe96fd18258d73c"
))
#println(datadep"DE_Tha_Jun_2014.rda")
ENV["DATADEPS_ALWAYS_ACCEPT"]="true" # avoid question to download
DE_Tha_Jun_2014 = first(values(load(joinpath(datadep"DE_Tha_Jun_2014.rda/DE_Tha_Jun_2014.rda"))))
nothing
tha = DE_Tha_Jun_2014
set_datetime_ydh!(tha)

thaf = copy(tha);
setinvalid_qualityflag!(thaf);
setinvalid_range!(thaf,
     :PPFD => (200, Inf),
     :ustar => (0.2, Inf),
     :LE =>(0, Inf),
     :VPD => (0.01, Inf)
     );
setinvalid_nongrowingseason!(thaf, 0.4);
setinvalid_afterprecip!(thaf; min_precip=0.02, hours_after=24);

thas = subset(thaf, :valid)

function tmpf()
    Dl=0.01
    aerodynamic_conductance!(thas; Gb_model=Su2001(),
        Dl, LAI=thal.LAI, zh=thal.zh, zr=thal.zr);
    surface_conductance!(thas, InversePenmanMonteith(); G=thas.G);
end

# tha48 and thal see runtests.jl

show(thaf.wind[1:48])
show(thaf.VPD[1:48])
show(thaf.LE[1:48])
show(thaf.Rn[1:48])
show(thaf.G[1:48])

dfGPPd = @pipe tha |>
    transform(_, :datetime => ByRow(yearmonthday) => :ymd, copycols = false) |>
    groupby(_, :ymd) |>
    combine(_, :datetime => (x -> Date(first(x))) => :date, :GPP => mean => :GPPd, )
    #x -> x[!,2]

show(dfGPPd.GPPd)

GPPd = [11.288497760271033, 13.013025930772224, 12.851774960756302, 11.996453734696843, 11.635422044472458, 11.155685572574535, 10.774393322790273, 10.181774605065584, 11.257192575993637, 12.9423423493281, 12.352468963712454, 13.402045020057509, 9.53826415212825, 12.071680051895479, 13.692589149111882, 12.845505638824156, 12.378533909407755, 11.672167064607493, 10.401156075240579, 10.705716138705611, 10.207347450816693, 11.052016352768987, 13.54435911634937, 12.060648361220956, 7.758974596237143, 9.869706534050541, 12.998054057980577, 10.627359564105669, 8.685295419767499, 10.874667977293333]
using Plots,StatsPlots
@df dfGPPd plot(:date, [:GPPd], xlab = "Date", ylab="GPP")


text:
using CSV, DataFrames

fname = only(ARGS)
df = CSV.File(fname) |> DataFrame
println(first(df, 3))

# W/o daemon:
# julia test.jl tsp_50.csv  11.65s

# time juliad test.jl tsp_50.csv
# 0.59s

max_stars_repo_name: jdurbin/sandbox

max_issues_repo_name: jdurbin/sandbox

max_forks_repo_name: jdurbin/sandbox


text:
function facepartition2(dmd,t,f,bcm,dim,elemtype,porder,nproc)

for i = 1:nproc
    fi = f[:,dmd[i].elempart[:]];
    f2t,~ = mkf2e(t[:,dmd[i].elempart[:]],elemtype,dim);

    # only on [intelem bndelem extelem]
    nelem = dmd[i].elempartpts;
    if length(nelem)>=3
        ne3 = sum(nelem[1:3]);
        ind = findall(f2t[1,:].<=ne3);
        f2t = f2t[:,ind];
    end

    # reorder so that boundary faces are last
    ina = findall(f2t[3,:].>0); # interior faces
    inb = findall(f2t[3,:].==0); # boundary faces
    #inc = sub2ind(size(fi), f2t[2,inb], f2t[1,inb]);
    inc = (f2t[1,inb].-1)*(size(fi,1)) .+ f2t[2,inb];
    fb = fi[inc]; # list of boundary indices

    fa = sort(unique(fb[:])); # boundary indices
    bcn = sort(unique(bcm[fa])); # a list of boundary conditions
    nbc = length(bcn);

    dmd[i].facepartpts = reshape([length(ina)],1,1);
    dmd[i].facepartbnd = reshape([0],1,1);
    ind = zeros(Int,1,length(fb));
    m = 1;
    for j=1:nbc # for each boundary condition bcn(j)
        bj = findall(bcm[:].==bcn[j]); # find all boundaries that have condition bcn(j)
        n = 0;
        for k = 1:length(bj) # for each boundary that has condition bcn(j)
            ii = findall(fb[:] .== bj[k]); # indices of the boundary bj(k)
            l = length(ii);
            n = n + l;
            ind[m:(m+l-1)] = ii;
            m = m + l;
        end
        dmd[i].facepartpts = [dmd[i].facepartpts n];
        dmd[i].facepartbnd = [dmd[i].facepartbnd bcn[j]];
    end
    # [interior faces, boundary faces]
    f2t = f2t[:,[ina; inb[ind[:]]]];

    dmd[i].facecon = faceconnectivity2(t[:,dmd[i].elempart[:]],f2t,dim,elemtype,porder);
end

return dmd

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
using Random: Random, AbstractRNG, Repetition, Sampler, rand!

struct SamplerRangeArbUInt{SP<:Sampler{ArbDigit}} <: Sampler{ArbUInt}
    start::ArbUInt
    gen::ArbUInt
    ndigits::Int
    nmaxdigits::Int
    highsp::SP
end

function SamplerRangeArbUInt(::Type{RNG}, r::AbstractUnitRange{ArbUInt}, N::Repetition=Val(Inf)) where {RNG <: AbstractRNG}
    isempty(r) && throw(ArgumentError("range must be non-empty"))
    gen = last(r) - first(r)
    ndigits = length(gen.data)
    hm = iszero(ndigits) ? zero(ArbDigit) : last(gen.data)
    highsp = Sampler(RNG, zero(ArbDigit):hm, N)
    nmaxdigits = max(ndigits, length(last(r).data), length(first(r).data))
    return SamplerRangeArbUInt(first(r), gen, ndigits, nmaxdigits, highsp)
end

Random.Sampler(::Type{RNG}, r::AbstractUnitRange{ArbUInt}, N::Repetition) where {RNG <: AbstractRNG} = SamplerRangeArbUInt(RNG, r, N)

Random.rand(rng::AbstractRNG, sp::SamplerRangeArbUInt) = rand!(rng, zero(ArbUInt), sp)

function Random.rand!(rng::AbstractRNG, a::ArbUInt, sp::SamplerRangeArbUInt)
    ndigits = sp.ndigits
    iszero(sp.ndigits) && return set_zero!(a)
    buf = a.data
    resize!(buf, sp.nmaxdigits)
    hm = last(sp.gen.data)
    while true
        rand!(rng, buf)
        hx = buf[ndigits] = rand(rng, sp.highsp)
        hx <= hm && break
        _less(<, buf, sp.gen.data) && break
    end
    add!(a, sp.start)
end


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

"""
    struct NewtonCG

Constructors:

* '''$(FUNCTIONNAME)(; fields...)'''

$(TYPEDFIELDS)

"""
@with_kw struct NewtonCG
    "amount of previous NewtonCG improvement guarding the lower
    bound to the improvement between consecutive cg iterations from
    the second NewtonCG step on"
    α::Float64=0.1

    "Number of total NewtonCG steps"
    steps::Int64=4

    "maximum number of cg iterations in the first NewtonCG step"
    i₀::Int64=5

    "maximum number of cg iterations from the second NewtonCG step on"
    i₁::Int64=50

    "LineSearcher that will be used after cg iterations are finished"
    linesearcher=StrongWolfe{Float64}()
end

mutable struct NewtonCGResults{O, Tx, Tf, M} <: Optim.OptimizationResults
    method::O
    initial_x::Tx
    minimizer::Tx
    initial_f::Tf
    minimum::Tf
    iterations::Int
    trace::M
    f_calls::Int
    g_calls::Int
    cg_iterations::Int
    absolute_reduction::Tf
end

function Base.show(io::IO, r::NewtonCGResults)
    println(io, "NewtonCG Results:")
    println(io, "------------------------------")
    println(io, "f value start: $(r.initial_f)")
    println(io, "f value end: $(r.minimum)")
    println(io, "absolute reduction: $(r.absolute_reduction)")
    println(io, "------------------------------")
    println(io, "f calls: $(r.f_calls)")
    println(io, "∇f calls: $(r.g_calls)")
    println(io, "------------------------------")
    f_hist, cg_hist = r.trace
    for (fₙ, cg_iter, n) in zip(f_hist, cg_hist, 0:r.iterations)
        println(io, "n: $n   fₙ: $fₙ    cg_iter: $cg_iter")
    end
end

function linesearch_args(
        f, ∇f, x::AbstractVector, Δx::AbstractVector, f_x, ∇f_x)
    # build univariate functions
    f_uni(α::Real) = f(x + α*Δx)
    df(α::Real) = dot(∇f(x + α*Δx), Δx)
    f_and_df(α::Real) = (f_uni(α), df(α))
    return (f_uni, df, f_and_df, 1.0, f_x, dot(∇f_x, Δx))
end

struct EvalCount
    f::Function
    counter::Base.Threads.Atomic{Int}
end

EvalCount(f::Function) = EvalCount(f, Base.Threads.Atomic{Int}(0))

function (F::EvalCount)(x::AbstractVector)
    Base.Threads.atomic_add!(F.counter, 1)
    F.f(x)
end

function _optimize(f::Function, ∇f!::Function, Σ̅⁻¹::Function,
        x₀::AbstractVector, optimizer::NewtonCG, optim_options)
    # fetch parameters from optimizer struct
    α = optimizer.α
    steps = optimizer.steps
    i₀ = optimizer.i₀
    i₁ = optimizer.i₁
    ls = optimizer.linesearcher
    # logging information
    cg_iterations=Int64[]
    f_history=Float64[]

    # gradient function of f
    ∇f(x) = ∇f!(similar(x), x)

    f = EvalCount(f)
    ∇f = EvalCount(∇f)

    # value of f before any optimization steps
    f⁰ = fⁿ⁻¹ = f(x₀)
    push!(f_history, f⁰)
    push!(cg_iterations, i₀)

    fⁿ = Δfⁿ = zero(f⁰)
    xₙ = x₀
    for n in 1:steps
        # preallocate/reset vector of descent
        Δx = zero(xₙ)
        ∇f_at_xₙ = ∇f(xₙ)

        # initialize cg_iterator with our mean fisher metric as the
        # hessian and our gradient of f at xₙ as our gradient
        # A = Σ̅⁻¹, b = ∇f(xₙ)
        cgiterator = cg_iterator!(Δx, Σ̅⁻¹(xₙ), ∇f_at_xₙ, initially_zero=true)
        if n == 1
            # do i₀ iterations of cg
            for (iteration, residual) in enumerate(cgiterator)
                if iteration >= i₀
                    break
                end
            end
        else
            fᵏ⁻¹ = fⁿ⁻¹
            # do at most i₁ cg iterations or move on if our improvement
            # is below α*100 percent of our previous NewtonCG step
            for (k, residual) in enumerate(cgiterator)
                fᵏ = f(xₙ - Δx)
                if k >= i₁ || abs(fᵏ - fᵏ⁻¹) < α*Δfⁿ
                    push!(cg_iterations, k)
                    break
                end
                fᵏ⁻¹ = fᵏ
            end
        end
        # finish NewtonCG step with line search in -Δx direction
        β, fⁿ = ls(linesearch_args(f, ∇f, xₙ, -Δx, fⁿ⁻¹, ∇f_at_xₙ)...)
        xₙ -= β*Δx
        Δfⁿ = abs(fⁿ - fⁿ⁻¹)
        fⁿ⁻¹ = fⁿ
        push!(f_history, fⁿ)
    end
    trace = (f_history=f_history, cg_iterations=cg_iterations)
    NewtonCGResults{typeof(optimizer), typeof(x₀), typeof(f⁰), typeof(trace)}(
        optimizer,
        x₀,
        xₙ,
        f⁰,
        fⁿ,
        steps,
        trace,
        f.counter[],
        ∇f.counter[],
        sum(cg_iterations),
        abs(fⁿ - f⁰)
    )
end

function _optimize(f::Function, ∇f!::Function, curvature::Function,
        x₀::AbstractVector, optimizer::Optim.AbstractOptimizer,
        optim_options::Optim.Options)
    optimize(f, ∇f!, x₀, optimizer, optim_options)
end


text:
using Test
using jInvVis
using jInv.Mesh
using Plots

# tests for regular mesh
domain = [0 1.1 0 1.0 ]
n      = [8 5 ]
Mr    = getRegularMesh(domain,n)
I     = randn(tuple(n...))
xc    = getCellCenteredGrid(Mr)

println("=== test viewImage2D  ===")

viewImage2D(xc[:,1],Mr)
xlabel!("intensity increases -->")
ylabel!("constant intensity")
# subplot(1,2,2)
viewImage2D(xc[:,2],Mr)
ylabel!("intensity increases -->")
xlabel!("constant intensity")
# println("==== all passed === ")


text:
using Random
using LinearAlgebra

using JuMP
import GLPK
using Plots

Random.seed!(314)

m = 12  # number of clients
n = 5  # number of facility locations

# Clients' locations
Xc = rand(m)
Yc = rand(m)

# Facilities' potential locations
Xf = rand(n)
Yf = rand(n)

# Fixed costs
f = ones(n);

# Distance
c = zeros(m, n)
for i in 1:m
    for j in 1:n
        c[i, j] = norm([Xc[i] - Xf[j], Yc[i] - Yf[j]], 2)
    end
end

# Display the data
scatter(Xc, Yc, label = "Clients", markershape=:circle, markercolor=:blue)
scatter!(Xf, Yf, label="Facility",
    markershape=:square, markercolor=:white, markersize=6,
    markerstrokecolor=:red, markerstrokewidth=2
)

# Create a JuMP model
ufl = Model(GLPK.Optimizer)

# Variables
@variable(ufl, y[1:n], Bin);
@variable(ufl, x[1:m, 1:n], Bin);

# Each client is served exactly once
@constraint(ufl, client_service[i in 1:m],
    sum(x[i, j] for j in 1:n) == 1
);

# A facility must be open to serve a client
@constraint(ufl, open_facility[i in 1:m, j in 1:n],
    x[i, j] <= y[j]
)

# Objective
@objective(ufl, Min, f'y + sum(c .* x));

# Solve the uncapacitated facility location problem with GLPK
optimize!(ufl)
println("Optimal value: ", objective_value(ufl))

x_ = value.(x) .> 1 - 1e-5
y_ = value.(y) .> 1 - 1e-5

# Display clients
p = scatter(Xc, Yc, markershape=:circle, markercolor=:blue, label=nothing)

# Show open facility
mc = [(y_[j] ? :red : :white) for j in 1:n]
scatter!(Xf, Yf,
    markershape=:square, markercolor=mc, markersize=6,
    markerstrokecolor=:red, markerstrokewidth=2,
    label=nothing
)

# Show client-facility assignment
for i in 1:m
    for j in 1:n
        if x_[i, j] == 1
           plot!([Xc[i], Xf[j]], [Yc[i], Yf[j]], color=:black, label=nothing)
        end
    end
end

display(p)

# Demands
a = rand(1:3, m);

# Capacities
q = rand(5:10, n);

# Display the data
scatter(Xc, Yc, label=nothing,
    markershape=:circle, markercolor=:blue, markersize= 2 .*(2 .+ a)
)

scatter!(Xf, Yf, label=nothing,
    markershape=:rect, markercolor=:white, markersize= q,
    markerstrokecolor=:red, markerstrokewidth=2
)

# Create a JuMP model
cfl = Model(GLPK.Optimizer)

# Variables
@variable(cfl, y[1:n], Bin);
@variable(cfl, x[1:m, 1:n], Bin);

# Each client is served exactly once
@constraint(cfl, client_service[i in 1:m], sum(x[i, :]) == 1)

# Capacity constraint
@constraint(cfl, capacity, x'a .<= (q .* y))

# Objective
@objective(cfl, Min, f'y + sum(c .* x));

# Solve the problem
optimize!(cfl)
println("Optimal value: ", objective_value(cfl))
@test termination_status(cfl) == MOI.OPTIMAL

x_ = value.(x) .> 1 - 1e-5
y_ = value.(y) .> 1 - 1e-5

# Display the solution
p = scatter(Xc, Yc, label=nothing,
    markershape=:circle, markercolor=:blue, markersize= 2 .*(2 .+ a)
)

mc = [(y_[j] ? :red : :white) for j in 1:n]
scatter!(Xf, Yf, label=nothing,
    markershape=:rect, markercolor=mc, markersize=q,
    markerstrokecolor=:red, markerstrokewidth=2
)

# Show client-facility assignment
for i in 1:m
    for j in 1:n
        if x_[i, j] == 1
            plot!([Xc[i], Xf[j]], [Yc[i], Yf[j]], color=:black, label=nothing)
            break
        end
    end
end

display(p)


text:
const uint64range = UnitRange(one(UInt64), typemax(UInt64))
const uint32range = UnitRange(one(UInt32), typemax(UInt32))

const mono32mask   = 0x1fffffff    # 536_870_911
const mono32offset = 0x44f619d0    # 1_156_979_152

mutable struct RomuQuadState
    wstate::UInt64
    xstate::UInt64
    ystate::UInt64
    zstate::UInt64
end

mutable struct RomuTrioState
    xstate::UInt64
    ystate::UInt64
    zstate::UInt64
end

mutable struct RomuDuoState
    xstate::UInt64
    ystate::UInt64
end

mutable struct RomuDuoJrState
    xstate::UInt64
    ystate::UInt64
end

mutable struct RomuQuad32State
    wstate::UInt32
    xstate::UInt32
    ystate::UInt32
    zstate::UInt32
end

mutable struct RomuTrio32State
    xstate::UInt32
    ystate::UInt32
    zstate::UInt32
end

mutable struct RomuMono32State
    xstate::UInt32
end

function setstate(::Type{RomuQuadState}, w::T=0, x::T=0, y::T=0, z::T=0) where {T<:UInt64}
    validate(RomuQuadState, w, x, y, z)
    return RomuQuadState(w, x, y, z)
end

function setstate(::Type{RomuTrioState}, x::T=0, y::T=0, z::T=0) where {T<:UInt64}
    validate(RomuTrioState, x, y, z)
    return RomuTrioState(x, y, z)
end

function setstate(::Type{RomuDuoState}, x::T=0, y::T=0) where {T<:UInt64}
    validate(RomuDuoState, x, y)
    return RomuDuoState(x, y)
end

function setstate(::Type{RomuDuoJrState}, x::T=0, y::T=0) where {T<:UInt64}
    validate(RomuDuoJrState, x, y)
    return RomuDuoJrState(x, y)
end

function setstate(::Type{RomuQuad32State}, w::T=0, x::T=0, y::T=0, z::T=0) where {T<:UInt32}
    validate(RomuQuad32State, w, x, y, z)
    return RomuQuad32State(w, x, y, z)
end

function setstate(::Type{RomuTrio32State}, x::T=0, y::T=0, z::T=0) where {T<:UInt32}
    validate(RomuTrio32State, x, y, z)
    return RomuTrio32State(x, y, z)
end

function setstate(::Type{RomuMono32State}, x::T=0) where {T<:UInt32}
    validate(RomuMono32State, x)
    x = x & mono32mask
    x = x + mono32offset
    return RomuMono32State(x)
end

setstate(::Type{RomuQuadState}, x::Nothing) = setstate(RomuQuadState, rand(uint64range, 4)...)
setstate(::Type{RomuTrioState}, x::Nothing) = setstate(RomuTrioState, rand(uint64range, 3)...)
setstate(::Type{RomuDuoState}, x::Nothing) = setstate(RomuDuoState, rand(uint64range, 2)...)
setstate(::Type{RomuDuoJrState}, x::Nothing) = setstate(RomuDuoJrState, rand(uint64range, 2)...)
setstate(::Type{RomuQuad32State}, x::Nothing) = setstate(RomuQuad32State, rand(uint32range, 4)...)
setstate(::Type{RomuTrio32State}, x::Nothing) = setstate(RomuTrio32State, rand(uint32range, 3)...)
setstate(::Type{RomuMono32State}, x::Nothing) = setstate(RomuMono32State, (rand(uint32range) & mono32mask) + mono32offset)

validate(::Type{RomuQuadState}, w::T=0, x::T=0, y::T=0, z::T=0) where {T<:UInt64} =
   if iszero(w) || iszero(x) || iszero(y) || iszero(z)
       throw(RomuException("All four states must be initialized to nonzero values."))
   else
       nothing
   end

validate(::Type{RomuTrioState}, x::T=0, y::T=0, z::T=0) where {T<:UInt64} =
   if iszero(x) || iszero(y) || iszero(z)
       throw(RomuException("All three states must be initialized to nonzero values."))
   else
       nothing
   end

validate(::Type{RomuDuoState}, x::T=0, y::T=0) where {T<:UInt64} =
   if iszero(x) || iszero(y)
       throw(RomuException("Both states must be initialized to nonzero values."))
   else
       nothing
   end

validate(::Type{RomuDuoJrState}, x::T=0, y::T=0) where {T<:UInt64} =
   if iszero(x) || iszero(y)
       throw(RomuException("Both states must be initialized to nonzero values."))
   else
       nothing
   end

validate(::Type{RomuQuad32State}, w::T=0, x::T=0, y::T=0, z::T=0) where {T<:UInt32} =
   if iszero(w) || iszero(x) || iszero(y) || iszero(z)
       throw(RomuException("All four states must be initialized to nonzero values."))
   else
       nothing
   end

validate(::Type{RomuTrio32State}, x::T=0, y::T=0, z::T=0) where {T<:UInt32} =
   if iszero(x) || iszero(y) || iszero(z)
       throw(RomuException("All three states must be initialized to nonzero values."))
   else
       nothing
   end

validate(::Type{RomuMono32State}, x::T=0) where {T<:UInt32} =
   if iszero(x)
       throw(RomuException("The state must be initialized to a nonzero value."))
   else
       nothing
   end


text:
function mw_h_cov_ha_i0_i1_i2(psi, param, nsim, rmatrix_strig)
#mw_h_cov_ha_i0_i1_i2 - One line description of what the function or script performs (H1 line)
# Generate nsim HA covariance mtrices for h-test for
#   (1). I(0) model.
#   (2). RW model.
#   (3). Integrated RW Model
#   Three covariance matrices are returned.
#   These are the ingredients to the covariance matrices for the local level and integrated local level models.
#
# Syntax:  [output1,output2] = function_name(input1,input2,input3)
#
# Inputs:
#   psi::Array{Float64,2} - dct transform.
#   param::Float64 - local parameters.
#   nsim::Int64 - number of simulations.
#   rmatrix_string = standard random variables.
#
# Outputs:
#   cov_ha_i0 - MATRIX version of cov_ha.
#   cov_ha_i1 - MATRIX version of cov_ha for summed process.
#   cov_ha_i2 - MATRIX version of cov_ha for double summed process.
#
#
#
#
  t = size(psi, 1)
  q = size(psi, 2)
  #qq = Int(q*(q+1)/2)
  # Get draws of variance function
  w = rmatrix_strig
  w = cumsum(w[1:t, 1:nsim], 1)/sqrt(t) # Random Walk
  w = w*(param/sqrt(q)) # Scaled
  sdpath_mat = exp.(w)
  # Construct summation matrix
  sum_mat = zeros(t, t)
  for i ∈ 1:t
    sum_mat[i, 1:i] = ones(1, i)
  end
  sum_mat2 = sum_mat*sum_mat # Double integration
  cov_ha_i0 = NaN*zeros(q, q, nsim)
  cov_ha_i1 = cov_ha_i0
  cov_ha_i2 = cov_ha_i1
  psi_0 = psi'
  psi_1 = psi'sum_mat
  psi_2 = psi'sum_mat2
  for i ∈ 1:nsim
    tmp = psi_0.*sdpath_mat[:,i]'
    tmp = tmp*tmp'
    cov_ha_i0[:, :, i] = Hermitian(tmp, :L)
    tmp = psi_1.*sdpath_mat[:,i]'
    tmp = tmp*tmp'
    cov_ha_i1[:, :, i] = Hermitian(tmp, :L)
    tmp = psi_2.*sdpath_mat[:,i]'
    tmp = tmp*tmp'
    cov_ha_i2[:, :, i] = Hermitian(tmp, :L)
  end
return cov_ha_i0, cov_ha_i1, cov_ha_i2
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

##    Error messages for Unicode / UTF support

const UTF_ERR_SHORT             = "invalid UTF-8 sequence starting at index <<1>> (0x<<2>> missing one or more continuation bytes)"
const UTF_ERR_INVALID_INDEX     = "invalid character index"

mutable struct UnicodeError <: Exception
    errmsg::AbstractString      ##< A UTF_ERR_ message
    errpos::Int32               ##< Position of invalid character
    errchr::UInt32              ##< Invalid character
end

show(io::IO, exc::UnicodeError) = print(io, replace(replace(string("UnicodeError: ",exc.errmsg),
    "<<1>>",string(exc.errpos)),"<<2>>",hex(exc.errchr)))

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia


text:
module MultiFEFunctions

using Gridap
using Gridap.Helpers

export MultiFEFunction

import Base: length
import Base: getindex
import Base: iterate
import Base: zero
import Gridap.FESpaces: free_dofs
import Gridap.FESpaces: FEFunction
import Gridap: FEFunctionForEval
import Gridap: restrict

struct MultiFEFunction
  fields::Vector{<:FEFunction}
  free_dofs_all_fields::AbstractVector
end

function MultiFEFunction(
  free_dofs_all_fields::AbstractVector, fespaces::MultiFESpace)
  fields = [
    FEFunction(U,restrict_to_field(fespaces,free_dofs_all_fields,i))
    for (i,U) in enumerate(fespaces) ]
  MultiFEFunction(fields,free_dofs_all_fields)
end

function MultiFEFunctionForEval(
  free_dofs_all_fields::AbstractVector, fespaces::MultiFESpace)
  fields = [
    FEFunctionForEval(U,restrict_to_field(fespaces,free_dofs_all_fields,i))
    for (i,U) in enumerate(fespaces) ]
  MultiFEFunction(fields,free_dofs_all_fields)
end

function FEFunction(
  fespaces::MultiFESpace, free_dofs_all_fields::AbstractVector)
  MultiFEFunction(free_dofs_all_fields,fespaces)
end

function FEFunctionForEval(
  fespaces::MultiFESpace, free_dofs_all_fields::AbstractVector)
  MultiFEFunctionForEval(free_dofs_all_fields,fespaces)
end

free_dofs(self::MultiFEFunction) = self.free_dofs_all_fields

length(self::MultiFEFunction) = length(self.fields)

getindex(self::MultiFEFunction,field::Integer) = self.fields[field]

iterate(self::MultiFEFunction) = iterate(self.fields)

iterate(self::MultiFEFunction,state) = iterate(self.fields,state)

function zero(U::MultiFESpace{E}) where E
  n = num_free_dofs(U)
  x = zeros(E,n)
  MultiFEFunction(x,U)
end

zero(U::Vector{<:FESpace}) = zero(MultiFESpace(U))

function restrict(uh::MultiFEFunction,trian::BoundaryTriangulation)
  [ restrict(ui.cellfield,trian) for ui in uh.fields ]
end

function restrict(uh::MultiFEFunction,trian::SkeletonTriangulation)
  [ restrict(ui.cellfield,trian) for ui in uh.fields ]
end

# Pretty printing

import Base: show

function show(io::IO,self::MultiFEFunction)
  s = "MultiFEFunction object with $(length(self)) fields"
  print(io,s)
end

function show(io::IO,::MIME"text/plain",self::MultiFEFunction)
  show(io,self)
  print(io,":")
  print(io,"\n nfree (all fields): $(length(self.free_dofs_all_fields))")
  for (i,U) in enumerate(self.fields)
    print(io,"\nfield $i:\n")
    show(io,"text/plain",U)
  end
end

end # module MultiFEFunctions


text:
module CellSegmentation

import AxisArrays
import ColorTypes
import ImageBinarization
import ImageFiltering
import ImageMorphology
import ImageSegmentation
import Interpolations

using Statistics: quantile, median

export
    CellImageAxes,
    CellImage,
    CellImageNoAxes,
    cell_image,
    median_filter,
    get_background_mask,
    remove_edges,
    nearest_neighbor_resize,
    segment_via_background,
    extract_segment

# Type for an axis to index a cell image: a Height x Width x Channel axes,
#   where the Height and Width dimensions are integer-indexed and Channels are
#   named by symbols
CellImageAxes = Tuple{
    AxisArrays.Axis{:height,UnitRange{Int64}},
    AxisArrays.Axis{:width,UnitRange{Int64}},
    AxisArrays.Axis{:channel,Array{Symbol,1}}
}

# Type for a cell image: an AxisArray indexed using a CellImageAxis
CellImage = (
    AxisArrays.AxisArray{T,3,S,CellImageAxes}
    where S <: Array{T,3}
    where T <: ColorTypes.Gray{R}
    where R <: Union{ColorTypes.Fractional, Bool}
)

# Certain operations may strip axes, making this a useful type alias
CellImageNoAxes = (
    Array{T,3} where T <: ColorTypes.Gray{S} where S <: ColorTypes.Fractional
)

"""
Wrap an array containing a cell image into a nicely-indexed CellImage.
"""
function cell_image(
    image::CellImageNoAxes,
    channel_names::AbstractArray{Symbol,1}
)::CellImage
    (height, width, _) = size(image)
    return AxisArrays.AxisArray(
        image,
        AxisArrays.Axis{:height}(1:height),
        AxisArrays.Axis{:width}(1:width),
        AxisArrays.Axis{:channel}(channel_names)
    )
end

"""
Smooth an image using a windowed median, with window length as a fraction of image size.
"""
function median_filter(
    image::AbstractArray;
    window_size::Real=0.01,
    border_size::Real=0.1
)::AbstractArray
    # determine window size
    height, width = size(image)
    window_px = Int(floor(window_size * (height + width) / 2))
    if window_px % 2 == 0
        window_px += 1
    end
    pad_px = Int(round(border_size * (height + width) / 2))

    # apply filter
    border = ImageFiltering.Fill(0, pad_px)
    ImageFiltering.mapwindow(median, image, (window_px, window_px); border=border)
end

"""
Finds the background of a grayscale image, resulting in a black and white image.

Binarizes based on quantile of nonzero values.

Uses median smoothing and morphological opening for improved smoothness.

An alternative to something like
    `.!ImageBinarization.binarize(image, ImageBinarization.Otsu())``
"""
function get_background_mask(
    image::AbstractArray;
    median_filter_window_size::Real=0.005,
    background_nonzero_quantile::Real=0.2,
    area_opening_window_size::Real=0.02
)::AbstractArray
    # median smoothing smooth
    image = median_filter(image; window_size=median_filter_window_size)

    # threshold
    background_threshold = quantile(image[image .> 0], background_nonzero_quantile)
    is_background = image .< background_threshold

    # close some holes
    height, width = size(image)
    hole_area_px = Int(round(area_opening_window_size * (height + width) / 2))^2
    is_background = ImageMorphology.area_opening(is_background; min_area=hole_area_px)
end

"""
Finds, isolates, and removes edges in a grayscale image.

Potentially useful to increase distance between bordering cell nuclei.

TODO: Examine if Otsu's method is best binarization for edges.
TODO: Make edge thickening (by morphological dilation) independent of image resolution.
"""
function remove_edges(
    image::AbstractArray;
    min_edge_diameter_size::Real=0.04,
    n_dilations::Integer=2
)::AbstractArray
    # use edge detection to find edges
    kern_y, kern_x = ImageFiltering.Kernel.ando5()
    grad_y = ImageFiltering.imfilter(image, kern_y)
    grad_x = ImageFiltering.imfilter(image, kern_x)
    grad = (grad_y.^2) .+ (grad_x.^2)

    # binarize edges
    edge_mask = ImageBinarization.binarize(Bool, grad, ImageBinarization.Otsu())

    # drop tiny edge pieces
    height, width = size(image)
    min_diameter_px = Int(round(min_edge_diameter_size * (height + width) / 2))
    edge_mask = ImageMorphology.diameter_opening(
        edge_mask; min_diameter=min_diameter_px
    )

    # thicken edges
    edge_mask = ImageMorphology.dilate(edge_mask)

    res = copy(image)
    res[edge_mask] .= zero(eltype(res))
    return res
end

"""
Resize using nearest-neighbor Interpolation rather than
whatever ImageTransformations.imresize uses (which averages out integers).
"""
function nearest_neighbor_resize(
    image::AbstractArray,
    new_size::Tuple{Integer, Integer}
)::AbstractArray
    h_old, w_old = size(image)
    h_new, w_new = new_size
    h_range = LinRange(1, h_old, h_new)
    w_range = LinRange(1, w_old, w_new)
    int = Interpolations.interpolate(
        image, Interpolations.BSpline(Interpolations.Constant())
    )
    result = int(h_range, w_range)
end

"""
Segment nucleui or cells using marked watershed method.

If markers are not provided, they are established using a distance from the
    background, thresholded to a quantile of foreground distances from background.

Using a different background for markers is supported as well.
"""
function segment_via_background(
    background::AbstractArray;
    marker_background::Union{AbstractArray, Nothing}=nothing,
    markers::Union{AbstractArray, Nothing}=nothing,
    watershed_marker_quantile::Real=0.4
)::ImageSegmentation.SegmentedImage
    dist = (
        background
        |> ImageMorphology.feature_transform
        |> ImageMorphology.distance_transform
    )
    if !isa(markers, Nothing)
        if !isa(marker_background, Nothing)
            throw(
                ArgumentError("Cannot specify both `marker_background` and `markers`")
            )
        end
    else
        if isa(marker_background, Nothing)
            marker_dist = dist
        else
            marker_dist = (
                marker_background
                |> ImageMorphology.feature_transform
                |> ImageMorphology.distance_transform
            )
        end
        marker_threshold = quantile(
            marker_dist[marker_dist .> 0], watershed_marker_quantile
        )
        markers = ImageMorphology.label_components(marker_dist .>= marker_threshold)
    end
    segmentation = ImageSegmentation.watershed(-dist, markers, mask=.!background)
end

"""
Extract the contents of a segment mask and place it on a black background.
"""
function extract_segment(
    segment_mask::AbstractArray,
    full_image::CellImage
)::CellImage
    segment_indices = findall(segment_mask)
    idx_min = minimum(segment_indices)
    idx_max = maximum(segment_indices)
    shifted_indices = segment_indices .- (idx_min - CartesianIndex(1, 1),)
    max_dim = maximum((idx_max - idx_min).I) + 1
    channel_names = full_image.axes[3].val
    extracted_image = zeros(
        ColorTypes.Gray{ColorTypes.N0f8},
        (max_dim, max_dim, length(channel_names))
    )
    extracted_image = cell_image(extracted_image, channel_names)
    extracted_image[shifted_indices, :] = full_image[segment_indices, :]
    return extracted_image
end

end


text:
using GLVisualize, GeometryTypes, GLAbstraction, ModernGL, FileIO, Reactive

if !isdefined(:runtests)
	window = glscreen()
	timesignal = loop(linspace(0f0, 1f0, 360))
end

let
rotation_angle = const_lift(*, timesignal, 2f0*pi)
rotation = map(rotationmatrix_z, rotation_angle)

const b = Point3f0[(rand(Point3f0)*2)-1 for i=1:64]

sprites = visualize(
	(SimpleRectangle(0f0,0f0,0.5f0, 0.5f0), b),
	billboard=true, image=loadasset("doge.png"),
	model=rotation
)

_view(sprites, window, camera=:perspective)
end

if !isdefined(:runtests)
	renderloop(window)
end


text:
export DepthMap

import ImageView

type DepthMap
  camera :: M34
  depth :: Array{Float32, 2}
  nxcorr :: Array{Float32, 2}
end

function DepthMap(view, nbrs, voi, w = 3)
  cam = view.camera
  im = view.image
  (nc, nx, ny) = size(im)
  mn = LibAminda.mean_and_inverse_deviation(im, w)

  # determine depth range, resolution
  bnds = bounds(cam, voi)
  near = bnds[1][3]
  far = bnds[2][3]
  dz = near * camera_resolution(cam)
  nz = ceil(Int, (far - near) / dz)
  dz = (far - near) / nz

  max_nxcorr = LibAminda.fill(nx, ny, -1.0)
  depth = LibAminda.fill(nx, ny, near)
  for k in 1:nz

    z = near + dz*(k-0.5)
    nxcorr = LibAminda.fill(nx, ny, -1.0)
    for nbr in nbrs
      cam2 = nbr.camera
      hom = Array(homography(cam, cam2, z))
      im2 = LibAminda.map_to_plane(nbr.image, hom, nx, ny)
      mn2 = LibAminda.mean_and_inverse_deviation(im2, w)
      nxc = LibAminda.normalized_cross_correlation(im, mn, im2, mn2, w)
      nxcorr = LibAminda.maximum(nxcorr, nxc)
    end

    LibAminda.update_depth(nxcorr, max_nxcorr, depth, z)
  end

  return DepthMap(cam, depth, max_nxcorr)
end

max_stars_repo_name: twadleigh/aminda

max_issues_repo_name: twadleigh/aminda

max_forks_repo_name: twadleigh/aminda


text:
"""
    mapmany(f, a...)

Like `map`, but `f(x)` for each `x ∈ a` may return an arbitrary number of values to insert
into the output.

# Example

```jldoctest
julia> mapmany(x -> 1:x, [1,2,3])
6-element Array{Int64,1}:
 1
 1
 2
 1
 2
 3
```
"""
function mapmany(f, a)
    T = eltype(promote_op(f, eltype(a)))
    out = T[]
    for x ∈ a
        append!(out, f(x))
    end
    return out
end

mapmany(f, a, b, c...) = mapmany((x,y,z...)->f(x,y,z...), zip(a, b, c...))

"""
    flatten(a)

Takes a collection of collections `a` and returns a collection containing all the elements
of the subcollecitons of `a`. Equivalent to `mapmany(idenity, a)`.

# Example

```jldoctest
julia> flatten([1:1, 1:2, 1:3])
6-element Array{Int64,1}:
 1
 1
 2
 1
 2
 3
```
"""
flatten(x) = mapmany(identity, x)

# mapview for lazy map?


text:
## Horseshoe methods
#==========================================================================================#

"""
Compute the quarter point between two points in the x-z plane.
"""
quarter_point(p1, p2) = weighted_vector(p1, p2, SVector(1/4, 0, 1/4))

"""
Compute the 3-quarter point between two points in the x-z plane.
"""
three_quarter_point(p1, p2) = weighted_vector(p1, p2, SVector(3/4, 0, 3/4))

collocation_point(p1, p2, p3, p4) = ( three_quarter_point(p1, p2) + three_quarter_point(p4, p3) ) / 2
bound_leg(p1, p2, p3, p4) = (quarter_point(p1, p2), quarter_point(p4, p3))

"""
    bound_leg(panel :: Panel3D)

Compute the bound leg for a `Panel3D`, for horseshoes/vortex rings.
"""
bound_leg(panel :: Panel3D) = bound_leg(panel.p1, panel.p2, panel.p3, panel.p4)

"""
    collocation_point(panel :: Panel3D)

Compute the collocation point of a `Panel3D` for horseshoes/vortex rings, which is the 3-quarter point on each side in the ``x``-``z`` plane.
"""
horseshoe_point(panel :: Panel3D) = collocation_point(panel.p1, panel.p2, panel.p3, panel.p4)

"""
    Line(r1 :: SVector{3,<: Real}, r2 :: SVector{3,<: Real})

A composite type consisting of two vectors to define a line.
"""
struct Line{T <: Real}
    r1 :: SVector{3,T}
    r2 :: SVector{3,T}
end

Line(r1 :: FieldVector{3,T}, r2 :: FieldVector{3,T}) where T <: Real = Line{T}(r1, r2)

r1(line :: Line) = line.r1
r2(line :: Line) = line.r2
vector(line :: Line) = r2(line) - r1(line)
center(line :: Line) = (r1(line) + r2(line)) / 2

points(lines :: Vector{<: Line}) = [ r1.(lines); [(r2 ∘ last)(lines)] ]

transform(line :: Line, rotation, translation) = let trans = Translation(translation) ∘ LinearMap(rotation); Line((trans ∘ r1)(line), (trans ∘ r2)(line)) end

r1(r, line :: Line) = r - r1(line)
r2(r, line :: Line) = r - r2(line)

bound_leg_velocity(a, b, Γ)    = Γ/4π * (1/norm(a) + 1/norm(b)) * a × b / (norm(a) * norm(b) + dot(a, b))
trailing_leg_velocity(r, Γ, u) = Γ/4π * normalize(r) × u / (norm(r) - dot(r, u))

trailing_legs_velocities(a, b, Γ, u) = trailing_leg_velocity(a, Γ, u) - trailing_leg_velocity(b, Γ, u)
total_horseshoe_velocity(a, b, Γ, u) = bound_leg_velocity(a, b, Γ) + trailing_legs_velocities(a, b, Γ, u)

horseshoe_velocity(r, line :: Line, Γ, direction) = total_horseshoe_velocity(r1(r, line), r2(r, line), Γ, direction)

## Arrays of vortex lines
#==========================================================================================#

"""
Placeholder. Vortex rings and horseshoes basically have the same methods, and are arrays of vortex lines.
"""
abstract type AbstractVortexArray end

"""
A horseshoe type consisting of a bound leg of type Line represening a vortex line.
"""
struct Horseshoe{T <: Real} <: AbstractVortexArray
    bound_leg         :: Line{T}
    collocation_point :: SVector{3,T}
    chord             :: T
end

"""
    bound_leg(horseshoe :: Horseshoe)

Getter for bound leg field of a `Horseshoe`.
"""
bound_leg(horseshoe :: Horseshoe) = horseshoe.bound_leg
collocation_point(horseshoe :: Horseshoe) = horseshoe.collocation_point

r1(r, horseshoe :: Horseshoe) = r1(r, bound_leg(horseshoe))
r2(r, horseshoe :: Horseshoe) = r2(r, bound_leg(horseshoe))

"""
Return a `Horseshoe` bound leg corresponding to a `Panel3D`.
"""
horseshoe_line(panel :: Panel3D, drift = SVector(0., 0., 0.)) = let (r1, r2) = bound_leg(panel);
    Horseshoe(Line(r1, r2), horseshoe_point(panel) .+ drift, (norm ∘ average_chord)(panel)) end

"""
Compute the midpoint of the bound leg of a `Horseshoe`.
"""
bound_leg_center(horseshoe :: Horseshoe) = (center ∘ bound_leg)(horseshoe)

"""
Compute the direction vector of the bound leg of a `Horseshoe`.
"""
bound_leg_vector(horseshoe :: Horseshoe) = (vector ∘ bound_leg)(horseshoe)

"""
    velocity(r, horseshoe, Γ, V_hat)

Compute the induced velocities at a point ``r`` of a given Horseshoe with constant strength ``Γ`` and trailing legs pointing in a given direction ``\\hat V``.
"""
function velocity(r, horseshoe :: Horseshoe, Γ :: Real, V_hat, finite_core = false)
    if finite_core
        width = (norm ∘ bound_leg_vector)(horseshoe)
        ε = max(horseshoe.chord, width) # Wrong core size? Consider options...
        horseshoe_velocity(r, bound_leg(horseshoe), Γ, V_hat, ε)
    else
        horseshoe_velocity(r, bound_leg(horseshoe), Γ, V_hat)
    end
end

max_forks_repo_name: GodotMisogi/AeroMDAO


text:
"""
SamplerState for Hamiltonian Monte Carlo

Fields:
    - x: the current state
    - p: the current momentum
    - stepsize
    - niters: number of iterations
    - mass: mass matrix
"""
type HMCState <: SamplerState
    x::Array{Float64}
    p::Array{Float64} # included in sampler state to allow autoregressive updates.
    stepsize::Float64
    niters::Int64
    mass
    function HMCState(x::Array{Float64},stepsize::Float64;p=randn(length(x)),niters=10,mass=1.0)
        if isa(mass,Number)
          mass = mass * ones(length(x))
        end
        new(x,p,stepsize,niters,mass)
    end
end

"""
sample! performs one HMC update
"""
function sample!(s::HMCState,llik,grad)
  # hamiltonian monte carlo (radford neal's version)
  nparams = length(s.x)
  mass = s.mass
  stepsize = s.stepsize
  niters = s.niters

  s.p = sqrt(mass).*randn(nparams)
  curx = s.x
  curp = s.p
  s.p += .5*stepsize * grad(s.x)
  for iter = 1:niters
    s.x += stepsize * s.p./mass
    s.p += (iter<niters ? stepsize : .5*stepsize) * grad(s.x) # two leapfrog steps rolled in one unless at the end.
  end

  logaccratio = llik(s.x) - llik(curx) -.5*sum((s.p.*s.p - curp.*curp)./mass)[1]
  if 0.0 > logaccratio - log(rand())
      #reject
      s.x = curx
      s.p = curp
  else
      #accept
      #negate momentum for symmetric Metropolis proposal
      s.p = -s.p
  end
  return s
end


text:
@deprecate WassersteinExponentialKernel() ExponentialKernel(; metric=Wasserstein())
@deprecate MixtureWassersteinExponentialKernel() ExponentialKernel(;
    metric=MixtureWasserstein()
)


text:
module Sparse23raw
export PreLapl2, PreLapl3
using SparseArrays

function PreLapl2(size::Int64)
    size2 = size * size

    nc = 5*(size - 2)^2 + 16*(size - 2) + 12
    colptr = Array{Int64}(undef, size2+1)
    rowval = Array{Int64}(undef, nc)
    nzval = Array{Float64}(undef, nc)

    h = 1. / (size - 1)
    h2 = h * h
    cd = -4. / h2
    hd = 1. / h2

    colptr[1] = 1
    count = 1
    for j=1:size
        for i=1:size
            l = (j-1) * size + i

            rowval[count] = l
            nzval[count] = cd
            numc = 1
            count += 1

            for i1 = -1:2:1
                if i+i1 > 0 && i+i1 <= size
                    rowval[count] = l + i1
                    nzval[count] = hd
                    numc += 1
                    count += 1
                end

                if j+i1 > 0 && j+i1 <= size
                    rowval[count] = l + i1 * size
                    nzval[count] = hd
                    numc += 1
                    count += 1
                end
            end

            colptr[l+1] = colptr[l] + numc
        end
    end

    M = SparseMatrixCSC{Float64, Int64}(size2, size2, colptr, rowval, nzval)
    return M, size2, count-1
end

function PreLapl3(size)
    size2 = size * size
    size3 = size2 * size

    nc = 7*(size - 2)^3 + 36*(size - 2)^2 + 60*(size - 2) + 8*4
    colptr = Array{Int64}(undef, size3+1)
    rowval = Array{Int64}(undef, nc)
    nzval = Array{Float64}(undef, nc)

    h = 1. / (size - 1)
    h2 = h * h
    cd = -6. / h2
    hd = 1. / h2

    colptr[1] = 1
    count = 1
    for k=1:size
        for j=1:size
            for i=1:size
                l = (k - 1) * size2 + (j - 1) * size + i

                rowval[count] = l
                nzval[count] = cd
                numc = 1
                count += 1

                for i1 = -1:2:1
                    if i+i1 > 0 && i+i1 <= size
                        rowval[count] = l + i1
                        nzval[count] = hd
                        numc += 1
                        count += 1
                    end

                    if j+i1 > 0 && j+i1 <= size
                        rowval[count] = l + i1 * size
                        nzval[count] = hd
                        numc += 1
                        count += 1
                    end

                    if k+i1 > 0 && k+i1 <= size
                        rowval[count] = l + i1 * size2
                        nzval[count] = hd
                        numc += 1
                        count += 1
                    end
                end

                colptr[l+1] = colptr[l] + numc
            end
        end
    end

    M = SparseMatrixCSC{Float64, Int64}(size3, size3, colptr, rowval, nzval)
    return M, size^3, count-1
end
end

max_stars_repo_name: danielcort/benchmarks-python-julia-c

max_issues_repo_name: danielcort/benchmarks-python-julia-c

max_forks_repo_name: danielcort/benchmarks-python-julia-c


text:
@testset "simulation" begin
    @test WhittleLikelihoodInference.GaussianProcess(OU(1.0,1.0), 10, 1.0).X isa WhittleLikelihoodInference.Distributions.MvNormal
    @test WhittleLikelihoodInference.GaussianProcess(CorrelatedOU(1.0,1.0,0.5), 10, 1.0).X isa WhittleLikelihoodInference.Distributions.MatrixReshaped
    @test_throws ArgumentError simulate_gp(OU(1.0,1.0), -10, 1.0, 1)
    @test_throws ArgumentError simulate_gp(OU(1.0,1.0), 10, -1.0, 1)
    @test_throws ArgumentError simulate_gp(OU(1.0,1.0), 10, 1.0, -1)
    ts1 = simulate_gp(OU(1.0,1.0), 10, 1.0, 1)[1]
    ts2 = simulate_gp(CorrelatedOU(1.0,1.0,0.5), 10, 1.0, 1)[1]
    @test ndims(ts1) == 1
    @test length(ts1) == 10
    @test ndims(ts2) == 2
    @test length(ts2) == 10
end


text:
export addressing, @addressing

macro addressing(expr)
    expr |> addressing |> esc
end

function addressing(expr::Expr)
    wrapfunc = isfuncdef(expr) ? identity : awrap
    expr |> at |> wrapfunc
end

const addressarg = symbol("##address")

awrap(expr) = :(($addressarg -> $expr)(list(0)))

function at(expr::Expr)
    if expr.head == :function
        @assert length(expr.args) == 2
        @assert expr.args[1].head == :call
        name = expr.args[1].args[1]
        args = expr.args[1].args[2:end]
        body = expr.args[2]
        :(function $name($addressarg, $(args...)); $(at(body)); end)
    elseif expr.head == :->
        args = procargs(expr.args[1])
        body = expr.args[2]
        :(($addressarg, $(args...)) -> $(at(body)))
    elseif expr.head == :call && expr.args[1] in primitives
        Expr(:call, expr.args[1], map(at, expr.args[2:end])...)
    elseif expr.head == :call
        func = at(expr.args[1]) # Recurse for e.g. dot syntax.
        args = expr.args[2:end]
        rest = map(at, args)
        id = genaddr()
        :($(func)(cons($id, $addressarg), $(rest...)))
    elseif expr.head == :local
        @assert expr.args[1].head == symbol("=")
        @assert length(expr.args[1].args) == 2
        var = expr.args[1].args[1]
        val = expr.args[1].args[2]
        :(local $var = $(at(val)))
    elseif expr.head in [:comparison, :block, :if, :&&, :||, :vcat, :ref, :., :...]
        expr.head == :comparion && (@assert expr.args[2] in compops)
        Expr(expr.head, map(at, expr.args)...)
    elseif expr.head in [:line, :quote]
        expr
    else
        error("unhandled $(expr.head) expression in addressing transform")
    end
end

at(expr::Symbol) = expr
at(expr::LineNumberNode) = expr
at(expr::Number) = expr # Bool subclasses number apparantly.
at(expr::QuoteNode) = expr
at(expr::String) = expr

block(exprs...) = Expr(:block, exprs...)

# I don't really like this, but threading a counter through at() seems
# tedious. I think I have to maintain some state anyway, in order to
# ensure that addresses are unique across all @pp blocks.

_addr = 0

function genaddr()
    global _addr
    _addr += 1
end

function resetaddressingcounter()
    global _addr
    _addr = 0
end


text:
module TestTypeDef
using Test
using CategoricalArrays
using CategoricalArrays: DefaultRefType, level,  reftype, leveltype

@testset "CategoricalPool, a b c order" begin
    pool = CategoricalPool(
        [
            "a",
            "b",
            "c"
        ],
        Dict(
            "a" => DefaultRefType(1),
            "b" => DefaultRefType(2),
            "c" => DefaultRefType(3),
        )
    )

    @test isa(pool, CategoricalPool)

    @test isa(pool.index, Vector)
    @test length(pool.index) == 3
    @test pool.index[1] == "a"
    @test pool.index[2] == "b"
    @test pool.index[3] == "c"

    @test isa(pool.invindex, Dict)
    @test length(pool.invindex) == 3
    @test pool.invindex["a"] === DefaultRefType(1)
    @test pool.invindex["b"] === DefaultRefType(2)
    @test pool.invindex["c"] === DefaultRefType(3)

    @test isa(pool.order, Vector{DefaultRefType})
    @test length(pool.order) == 3
    @test pool.order[1] === DefaultRefType(1)
    @test pool.order[2] === DefaultRefType(2)
    @test pool.order[3] === DefaultRefType(3)

    @test leveltype("abc") === String
    @test leveltype(String) === String
    @test leveltype(1.0) === Float64
    @test leveltype(Float64) === Float64

    for i in 1:3
        x = CategoricalValue(i, pool)

        @test leveltype(x) === String
        @test leveltype(typeof(x)) === String
        @test reftype(x) === DefaultRefType
        @test reftype(typeof(x)) === DefaultRefType
        @test x isa CategoricalValue{String, DefaultRefType}

        @test isa(level(x), DefaultRefType)
        @test level(x) === DefaultRefType(i)

        @test isa(CategoricalArrays.pool(x), CategoricalPool)
        @test CategoricalArrays.pool(x) === pool

        @test typeof(x)(x) === x

        @test CategoricalValue(UInt8(i), pool) == x
    end
end

@testset "CategoricalPool, c b a order" begin
    pool = CategoricalPool(
        [
            "a",
            "b",
            "c"
        ],
        Dict(
            "a" => DefaultRefType(1),
            "b" => DefaultRefType(2),
            "c" => DefaultRefType(3),
        ),
        [
            DefaultRefType(3),
            DefaultRefType(2),
            DefaultRefType(1),
        ]
    )

    @test isa(pool, CategoricalPool)

    @test isa(pool.index, Vector)
    @test length(pool.index) == 3
    @test pool.index[1] == "a"
    @test pool.index[2] == "b"
    @test pool.index[3] == "c"

    @test isa(pool.invindex, Dict)
    @test length(pool.invindex) == 3
    @test pool.invindex["a"] === DefaultRefType(1)
    @test pool.invindex["b"] === DefaultRefType(2)
    @test pool.invindex["c"] === DefaultRefType(3)

    @test isa(pool.order, Vector{DefaultRefType})
    @test length(pool.order) == 3
    @test pool.order[1] === DefaultRefType(3)
    @test pool.order[2] === DefaultRefType(2)
    @test pool.order[3] === DefaultRefType(1)

    for i in 1:3
        y = CategoricalValue(i, pool)

        @test isa(level(y), DefaultRefType)
        @test level(y) === DefaultRefType(i)

        @test isa(CategoricalArrays.pool(y), CategoricalPool)
        @test CategoricalArrays.pool(y) === pool

        @test typeof(y)(y) === y

        @test CategoricalValue(UInt8(i), pool) == y
    end
end

end


text:
y1 = zeros(20);
y2 = zeros(20);
z1 = zeros(10,10);
z2 = zeros(10,10);

xpoints1 = [rand(10) for i=1:3]
xpoints2 = [rand(10).+1 for i=1:3]
ppoints = [rand(10) for i=1:3]

function force_lower_tri(z)
    for i=1:size(z,1)
        for j=1:size(z,2)
            i<j && (z[i,j]=0)
        end
    end
end

get_gorig!(forig) = (y,x,p) -> (y.=0;ForwardDiff.gradient!(y,x->forig(x,p),x))
get_sgorig!(forig,gsparsity) = (y,x,p) ->  (y.=0;y[1:length(gsparsity)] .= ForwardDiff.gradient(x->forig(x,p),x)[gsparsity])
get_horig!(forig) = (z,x,p) -> (z.=0;ForwardDiff.hessian!(z,x->forig(x,p),x); force_lower_tri(z))
get_shorig!(forig,hsparsity) = (z,x,p) ->  (z.=0;hess=ForwardDiff.hessian(x->forig(x,p),x); z[1:length(hsparsity)] .= [hess[i,j] for (i,j) in hsparsity])
get_jorig!(forig) = (z,x,p) -> (z.=0;ForwardDiff.jacobian!(z,(y,x)->forig(y,x,p),zeros(size(z,1)),x))
get_sjorig!(forig,jsparsity) = (y,x,p) ->  (y.=0;n=length(y);jac=ForwardDiff.jacobian((y,x)->forig(y,x,p),zeros(n),x); y[1:length(jsparsity)] .= [jac[i,j] for (i,j) in jsparsity])

const test_list_1 = [
    ("function-test-1-1",(x,p)->beta(erf(x[1]/x[2]/p[3])+p[3]*x[2],erf(x[9])^2)),
    ("function-test-1-2",(x,p)->0*x[1]),
    ("function-test-1-3",(x,p)->beta(cos(log(inv(inv(x[1])))),erfc(tanh(0*x[1])))),
    ("function-test-1-4",(x,p)->(0*x[1]^x[3]^p[1]+x[1])/x[9]/x[10]),
    ("function-test-1-5",(x,p)->(x[1]+1.)^x[2]*log(x[3])/tanh(x[2])),
    ("function-test-1-6",(x,p)->beta(2*logbeta(x[1],x[5]),beta(x[2],x[3]))),
    ("function-test-1-7",(x,p)->besselj0(exp(erf(-x[1])))),
    ("function-test-1-8",(x,p)->erfc((x[1]^2/x[2])^x[9]/x[10])),
    ("function-test-1-9",(x,p)->erfc(x[1])^erf(2.5x[2])),
    ("function-test-1-10",(x,p)->sin(1/x[1])),
    ("function-test-1-11",(x,p)->sum(isodd(i) ? sin(x[i])/cos(x[i+1]) : csc(x[i+1]/p[2]+x[i-1]/p[3]) for i=2:9)),
    ("function-test-1-12",(x,p)->sum(isodd(i) ? 4+x[i] : sin(x[1]) for i=1:4)*erf(x[1]/p[1])/p[10]/p[9]*p[5]),
    ("function-test-1-13",(x,p)->exp(x[2])/cos(x[1])^2+sin(x[1]^2)),
    ("function-test-1-14",(x,p)->airyai(exp(x[1]+x[2]*p[2]^8))),
    ("function-test-1-15",(x,p)->prod(x[1]+x[2] for i=1:4)),
    ("function-test-1-16",(x,p)->sin(x[9]inv(x[1])-x[8]inv(x[2]))),
    ("function-test-1-17",(x,p)->-(x[2]-x[1])*(x[1]*x[2])*(x[1]*x[2])/sum(x[i]+p[j] for i=1:3 for j=4:8)),
    ("function-test-1-18",(x,p)->sum(sum(isodd(i) ? x[1]/x[2] : cos(x[2]) for i =1:3) / sum(isodd(i) ? x[i]^x[j] : x[i]^2 for i=1:4) for j=1:3)),
    ("function-test-1-19",(x,p)->x[1]/log(x[2]^2+9.)),
    ("function-test-1-20",(x,p)->sum(cos(x[i]) for i=1:4)),
    ("function-test-1-21",(x,p)->beta(beta(tan(beta(x[1],1)+p[2]),cos(sin(x[2]))),x[3])),
    ("function-test-1-22",(x,p)->sum(f(x[1]) for f in [sin, cos, exp])),
    ("function-test-1-23",(x,p)->sum(1 +beta(sin(x[1]+erf(p[2])/p[1])^2,cos(p[1]*x[4])/x[5]^2-x[3])+sin(x[i]) for i=1:7)),
    ("function-test-1-24",(x,p)->beta(cos(beta(beta(x[1]^9,x[2]),x[2]*x[3])),sin(x[2]*x[3]/p[2])/p[1]))
]

const test_list_2 = [
    ("function-test-2-1",(y,x,p)->begin
     y[1] = x[1]/log(x[2]^2+9.);
     y[2] = sin(1/x[1]);
     y[3] = sin(exp(x[1]+x[2]*p[2]^8));
     end),
    ("function-test-2-2",(y,x,p)->begin
     y[1] = (0*x[1]^x[3]^p[1]+x[1])/x[9]/x[10];
     y[2] = beta(2*logbeta(x[1],x[5]),beta(x[2],x[3]));
     y[3] = besselj0(exp(erf(-x[1])));
     y[4] = erfc((x[1]^2/x[2])^x[9]/x[10]);
     end),
    ("function-test-2-3",(y,x,p)-> begin
     y[1] = erfc(x[1])^erf(2.5x[2]);
     y[2] = prod(x[1]+x[2] for i=1:4);
     y[3] = sin(x[9]inv(x[1])-x[8]inv(x[2]));
     y[4] = -(x[2]-x[1])*(x[1]*x[2])*(x[1]*x[2])/sum(x[i]+p[j] for i=1:3 for j=4:8);
     y[5] = sum(sum(x[1]/x[2] for i =1:3) / sum(x[i]^x[j] for i=1:4) for j=1:3)
     end),
 ]

function test_function_1(forig,name,xpoints,ppoints)
    ff = forig(Variable(),Parameter())
    f  = function_evaluator(ff)

    g! = gradient_evaluator(ff)
    sg!,gsparsity = sparse_gradient_evaluator(ff)

    h! = hessian_evaluator(ff)
    sh!,hsparsity = sparse_hessian_evaluator(ff)

    gorig!  = get_gorig!(forig)
    sgorig! = get_sgorig!(forig,gsparsity)
    horig!  = get_horig!(forig)
    shorig!  = get_shorig!(forig,hsparsity)

    @testset "$name" begin
        @test compare(forig,f,xpoints,ppoints)
        @test compare(gorig!,g!,y1,y2,xpoints,ppoints)
        @test compare(sgorig!,sg!,y1,y2,xpoints,ppoints)
        @test compare(horig!,h!,z1,z2,xpoints,ppoints)
        @test compare(shorig!,sh!,z1,z2,xpoints,ppoints)
    end
end

function test_function_2(forig,name,xpoints,ppoints)
    ff = Field()
    forig(ff,Variable(),Parameter())

    f  = field_evaluator(ff)

    j! = jacobian_evaluator(ff)
    sj!,jsparsity = sparse_jacobian_evaluator(ff)

    jorig!  = get_jorig!(forig)
    sjorig! = get_sjorig!(forig,jsparsity)

    @testset "$name" begin
        @test compare(forig,f,y1,y2,xpoints,ppoints)
        @test compare(jorig!,j!,z1,z2,xpoints,ppoints)
        @test compare(sjorig!,sj!,y1,y2,xpoints,ppoints)
    end
end

for (f,~,~) in SimpleNL.f_nargs_1
    forig = @eval (x,p)->$f(x[1])
    name = "basic-function-test-$f"
    test_function_1(forig,name, f != :acoth ? xpoints1 : xpoints2, ppoints)
end

for (name,forig) in test_list_1
    test_function_1(forig,name,xpoints1, ppoints)
end

for (name,forig) in test_list_2
    test_function_2(forig,name,xpoints1, ppoints)
end


text:
"""
    @until condition begin
        body
    end

Repeatedly evaluate a conditional expression, and continue evaluating the
body of the loop as long as the condition remains false.

# Examples
```jldoctest
julia> i = 1
1

julia> @until i > 4 begin
           println(i)
           i += 1
       end
1
2
3
4
```
"""
macro until(condition, body)
    quote
        while ! $condition
            $body
        end
    end |> esc
end


text:
"""
    pb, init_point = import_from_dat(instancepath::String, precondcstrspath::String)

Build the polynomial optimization problem described by the `instancepath` file,
with the dispensory preconditionning descritpion `precondcstrspath` along with
the iniitial point possibly provided in the file (defaults value is null).
"""
function import_from_dat(instancepath::String, precondcstrspath::String)
    init_point = Point()
    variables = Dict{String, Variable}()
    exponents = Dict{String, Exponent}()
    pb = Problem()

    instance_str = open(instancepath)
    l = jump_comments!(instance_str)

    ## Collect and define variables
    line = matchall(r"\S+", l)
    while line[1] == "VAR_TYPE" && !eof(instance_str)
        if line[2] == "REAL"
            variables[line[3]] = Variable(line[3], Real)
        elseif line[2] == "BOOL"
            variables[line[3]] = Variable(line[3], Bool)
        elseif line[2] == "CPLX"
            variables[line[3]] = Variable(line[3], Complex)
        else
            error("import_to_dat(): Unknown variable type $(line[2]) for variable $(line[3]).")
        end

        init_point[variables[line[3]]] = parse(line[5]) + im*parse(line[6])

        l = readline(instance_str)
        line = matchall(r"\S+", l)
    end

    ## Mark where objective defintion begins
    mark(instance_str)

    ## Move forward to the monomial definition section
    while line[1] != "MONO_DEF" && !eof(instance_str)
        l = readline(instance_str)
        line = matchall(r"\S+", l)
    end

    while line[1] == "MONO_DEF" && !eof(instance_str)
        exponame = line[2]
        var = variables[line[3]]
        if !haskey(exponents, exponame)
            exponents[exponame] = Exponent()
        end
        exponents[exponame] = product(exponents[exponame], Exponent(Dict(var=>Degree(parse(line[5]), parse(line[6])))))
        l = readline(instance_str)
        line = matchall(r"\S+", l)
    end

    ## Reset stream to the objective defintion
    reset(instance_str)
    l = readline(instance_str)
    line = matchall(r"\S+", l)

    ## Build polynomial objective
    p = Polynomial()
    while line[2] == "OBJ"
        λ = parse(line[5]) + im*parse(line[6])
        var1, var2 = line[3:4]
        if line[1] == "MONO"
            p += λ * exponents[var1]
        else
            p += λ * (var1!="NONE" ? variables[var1] : 1) * (var2!="NONE" ? variables[var2] : 1)
        end
        l = readline(instance_str)
        line = matchall(r"\S+", l)
    end
    set_objective!(pb, p)

    while line[1] != "MONO_DEF" && !eof(instance_str)
        cstrname = line[2]
        lb = -Inf
        ub = +Inf
        p = Polynomial()
        while line[2] == cstrname
            λ = parse(line[5]) + im*parse(line[6])
            var1, var2 = line[3:4]
            if line[1] == "MONO"
                p += λ * exponents[var1]
            elseif line[1] ∈ Set(["CONST", "LIN", "QUAD"])
                p += λ * (var1!="NONE" ? variables[var1] : 1) * (var2!="NONE" ? variables[var2] : 1)
            elseif line[1] == "UB"
                ub = real(λ)
            elseif line[1] == "LB"
                lb = real(λ)
            else
                error("import_to_dat(): Unknown variable type $(line[1]) for constraint $(line[2]).")
            end
            l = readline(instance_str)
            line = matchall(r"\S+", l)
        end
        add_constraint!(pb, String(cstrname), lb << p << ub)
        l = readline(instance_str)
        line = matchall(r"\S+", l)
    end

    precond_str = open(precondcstrspath)

    l = jump_comments!(precond_str)
    line = matchall(r"\S+", l)
    while !eof(precond_str)
        if line[2] == "SQRT"
            pb.constraints[line[1]].precond = :sqrt
        else
            error("import_from_dat(): Unknown preconditioning $(line[2]) for constraint $(line[1]).")
        end
        l = readline(precond_str)
        line = matchall(r"\S+", l)
    end
    return pb, init_point
end

"""
    l = jump_comments!(io::IOStream)

Jump comments, i.e. lines starting by '#', in the `io` stream, and return the
first non commented line.
"""
function jump_comments!(io::IOStream)
    l = ""
    while !eof(io)
        l = readline(io)
        ismatch(r"\s*#", l) || break
    end
    return l
end

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
export opt_gauss_newton!

"""
    (iter,resnorm)=opt_gauss_newton!(
        graph,
        objfun,
        discr;
        maxit = 100,
        logger = 0,
        errtype = :abserr,
        stoptol = 1e-6,
        cref = get_all_cref(graph),
        input = :A,
        γ0 = 1.0,
        linlsqr = :backslash,
        droptol = 0,
    )

Applies the Gauss–Newton algorithm to solve the nonlinear least squares problem:
Fit the output of the `graph` to the values of `objfun`, in the points `discr`.

The variable `graph` is modified during the iterations. The function returns the
iteration `iter` where it terminated, and the corresponding residual norm
`resnorm`.

The kwargs are as follows:

  - `maxit` determines the maximum number of iterations. If `logger` has a
    value >0, then the intermediate results are printed.
  - `errtype` determines the error type measured, see `adjust_for_errtype!`.
  - `stoptol` is the corresponding stopping tolerance.
  - `cref` is a `Vector{Tuple{Symbol,Int}}` that determines which coefficients
    of `graph` that are considered free variables and optimized.
  - `input` is the label corresponding to the input node of the graph.
  - The stepsize can be scaled with `γ0`.
  - `linlsqr` and `droptol` determines how the inner linear least squares
    problem is solved; see [`solve_linlsqr!`](@ref).
"""
function opt_gauss_newton!(
    graph,
    objfun,
    discr;
    maxit = 100,
    logger = 0,
    errtype = :abserr,
    stoptol = 1e-6,
    cref = get_all_cref(graph),
    input = :A,
    γ0 = 1.0,
    linlsqr = :backslash,
    droptol = 0,
)
    resnorm = Inf
    iter = maxit
    vals = init_vals_eval_graph!(graph, discr, nothing, input)
    objfun_vals = objfun.(discr)

    for j = 1:maxit
        F = eval_graph(graph, discr, vals = vals, input = input)
        res = F - objfun_vals
        Jac = eval_jac(graph, discr, cref, vals = vals, input = input)
        adjust_for_errtype!(Jac, res, objfun_vals, errtype)

        resnorm = norm(res)
        if (logger > 0)
            @show j, resnorm
        end
        if (resnorm < stoptol)
            iter = j
            break
        end

        d = solve_linlsqr!(Jac, res, linlsqr, droptol)
        x = get_coeffs(graph, cref)
        x -= γ0 * d
        set_coeffs!(graph, x, cref)
    end
    return (iter, resnorm)
end


text:
using Plots,Statistics
dim=10

function NeighborReturner(network_, i, j)
    neighbors=[]
    if network_[i,j-1]!=0
        push!(neighbors, network_[i,j-1]) # the left neighbor--> neighbor[1]
    end
    if i!=1 && network_[i-1,j]!=0 #if true --> the neighbors[2] is the upper neighbor
        push!(neighbors,network_[i-1,j])
    end
    return neighbors
end
function InitialLabelFinder(a, L)
    while L[a] != L[L[a]]
        a = L[a]
    end
    return L[a]
end
function percolation(dim,p)
    network_=zeros(Int, dim, dim)
    L=[1]
    S=[dim]
    counter=2
    for i in 1:dim
        network_[i,1]=1
    end
    for col in 2:dim
        for row in 1:dim
            if p>rand()
                neighbors=NeighborReturner(network_, row, col)
                #println(row, col,neighbors)
                if length(neighbors)==0
                    network_[row,col]=counter
                    push!(L, counter)
                    push!(S, 1)
                    counter+=1
                elseif length(neighbors)==1
                    network_[row,col]=InitialLabelFinder(neighbors[1], L)
                    S[InitialLabelFinder(neighbors[1], L)]+=1
                elseif length(neighbors)!=0 && InitialLabelFinder(neighbors[1],L)==InitialLabelFinder(neighbors[2],L)
                    network_[row,col]=InitialLabelFinder(neighbors[1], L)
                    S[InitialLabelFinder(neighbors[1], L)]+=1
                else
                    network_[row,col]=InitialLabelFinder(neighbors[1], L)
                    L[InitialLabelFinder(neighbors[2], L)]=InitialLabelFinder(neighbors[1], L)
                    S[InitialLabelFinder(neighbors[1], L)]+=1+S[InitialLabelFinder(neighbors[2], L)]
                    S[InitialLabelFinder(neighbors[2], L)]=0

                end
            end
        end
    end
    return network_,L
end
function percolation_check(network_,L,dim)
    first_col=[]
    last_col=[]
    for i in 1:dim
        if network_[i,1]!=0
            push!(first_col,InitialLabelFinder(network_[i,1],L))
        end
    end
    for i in 1:dim
        if network_[i,dim]!=0
            push!(last_col,InitialLabelFinder(network_[i,dim],L))
        end
    end
    col=rand(1:dim) #choosing a random entry of the matrix
    row=rand(1:dim)
    if network_[row,col]!=0 && InitialLabelFinder(network_[row,col],L) in intersect(first_col,last_col)
        return 1
    else
        return 0
    end
end
#network_=percolation(dim, p)
Meanlist = []
STDlist = []
for p in 0.25:0.01:1
    check=[]
    for run_num in 1:1000
        network_,L=percolation(dim,p)
        push!(check,percolation_check(network_,L,dim))
    end
    push!(STDlist, std(check))
    push!(Meanlist, mean(check))
end
Meanlist
scatter(0.25:0.01:1, Meanlist, yerr=STDlist, xlabel="P", ylabel="Average infinite Q", legend=false)
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter4\\Fig\\4.4_L=10.png")


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

"""
    module Modia3D.Signals

All supported signals are created and computed here.

# Main developers
Andrea Neumayr and Martin Otter, [DLR - Institute of System Dynamics and Control](https://www.dlr.de/sr/en)
"""
module Signals

export @signal

using  StaticArrays
import DataStructures
import ModiaMath

import Base
import Modia3D
import Modia3D.Basics

include("signals.jl")
end


text:
using Test
import ArchGDAL;
const AG = ArchGDAL;
import ImageCore
import ColorTypes

@testset "test_images.jl" begin
    @testset "Test Gray colors" begin
        @test eltype(AG.imread("data/utmsmall.tif")) ==
              ColorTypes.Gray{ImageCore.N0f8}
    end

    AG.read("gdalworkshop/world.tif") do dataset
        @testset "Test RGB colors" begin
            @test eltype(AG.imread(dataset)) == ColorTypes.RGB{ImageCore.N0f8}
        end
    end

    AG.create(
        AG.getdriver("MEM"),
        width = 2,
        height = 2,
        nbands = 4,
        dtype = UInt8,
    ) do dataset
        @testset "imview" begin
            @test eltype(AG.imview(AG.GCI_GrayIndex, ones(UInt8, 2, 2))) ==
                  ColorTypes.Gray{ImageCore.N0f8}
            @test eltype(AG.imview(AG.GCI_Undefined, ones(UInt8, 2, 2))) ==
                  ColorTypes.Gray{ImageCore.N0f8}
            @test eltype(AG.imview(AG.GCI_RedBand, ones(UInt8, 2, 2))) ==
                  ColorTypes.RGB{ImageCore.N0f8}
            @test eltype(AG.imview(AG.GCI_GreenBand, ones(UInt8, 2, 2))) ==
                  ColorTypes.RGB{ImageCore.N0f8}
            @test eltype(AG.imview(AG.GCI_BlueBand, ones(UInt8, 2, 2))) ==
                  ColorTypes.RGB{ImageCore.N0f8}
            @test_throws ErrorException AG.imview(AG.GCI_Max, ones(UInt8, 2, 2))
            @test eltype(AG.imview(AG.GPI_Gray, ones(UInt8, 2, 2))) ==
                  ColorTypes.Gray{ImageCore.N0f8}
            @test eltype(
                AG.imview(
                    AG.GPI_Gray,
                    ones(UInt8, 2, 2),
                    ones(UInt8, 2, 2),
                    ones(UInt8, 2, 2),
                ),
            ) == ColorTypes.Gray{ImageCore.N0f8}
            @test_throws ErrorException AG.imview(AG.GPI_HLS, ones(UInt8, 2, 2))
            @test_throws ErrorException AG.imview(
                AG.GPI_HLS,
                ones(UInt8, 2, 2),
                ones(UInt8, 2, 2),
                ones(UInt8, 2, 2),
            )
            @test_throws ErrorException AG.imview(
                AG.GPI_HLS,
                ones(UInt8, 2, 2),
                ones(UInt8, 2, 2),
                ones(UInt8, 2, 2),
                ones(UInt8, 2, 2),
            )
        end

        @testset "imread(color, ..., xoffset, yoffset, xsize, ysize)" begin
            @test eltype(
                AG.imread(AG.GCI_GrayIndex, AG.getband(dataset, 1), 0, 0, 2, 2),
            ) == ColorTypes.Gray{ImageCore.N0f8}
            @test eltype(
                AG.imread(AG.GPI_Gray, AG.getband(dataset, 1), 0, 0, 2, 2),
            ) == ColorTypes.Gray{ImageCore.N0f8}
            @test eltype(AG.imread(AG.GCI_GrayIndex, dataset, 1, 0, 0, 2, 2)) ==
                  ColorTypes.Gray{ImageCore.N0f8}
            @test eltype(AG.imread(AG.GPI_Gray, dataset, 1, 0, 0, 2, 2)) ==
                  ColorTypes.Gray{ImageCore.N0f8}
        end

        @testset "imread(color, ..., rows::UnitRange, cols::UnitRange)" begin
            @test eltype(
                AG.imread(AG.GCI_GrayIndex, AG.getband(dataset, 1), 1:2, 1:2),
            ) == ColorTypes.Gray{ImageCore.N0f8}
            @test eltype(
                AG.imread(AG.GPI_Gray, AG.getband(dataset, 1), 1:2, 1:2),
            ) == ColorTypes.Gray{ImageCore.N0f8}
            @test eltype(AG.imread(AG.GCI_GrayIndex, dataset, 1, 1:2, 1:2)) ==
                  ColorTypes.Gray{ImageCore.N0f8}
            @test eltype(AG.imread(AG.GPI_Gray, dataset, 1, 1:2, 1:2)) ==
                  ColorTypes.Gray{ImageCore.N0f8}
        end

        @testset "imread(color, ...)" begin
            @test eltype(AG.imread(AG.GCI_GrayIndex, AG.getband(dataset, 1))) ==
                  ColorTypes.Gray{ImageCore.N0f8}
            @test eltype(AG.imread(AG.GPI_Gray, AG.getband(dataset, 1))) ==
                  ColorTypes.Gray{ImageCore.N0f8}
            @test eltype(AG.imread(AG.GCI_GrayIndex, dataset, 1)) ==
                  ColorTypes.Gray{ImageCore.N0f8}
            @test eltype(AG.imread(AG.GPI_Gray, dataset, 1)) ==
                  ColorTypes.Gray{ImageCore.N0f8}
        end

        @testset "Test RGBA colors" begin
            AG.setcolorinterp!(AG.getband(dataset, 1), AG.GCI_RedBand)
            AG.setcolorinterp!(AG.getband(dataset, 2), AG.GCI_GreenBand)
            AG.setcolorinterp!(AG.getband(dataset, 3), AG.GCI_BlueBand)
            AG.setcolorinterp!(AG.getband(dataset, 4), AG.GCI_AlphaBand)
            @test eltype(AG.imread(AG.getband(dataset, 1))) ==
                  ColorTypes.RGB{ImageCore.N0f8}
            @test eltype(AG.imread(dataset, 1)) ==
                  ColorTypes.RGB{ImageCore.N0f8}
            @test eltype(AG.imread(AG.getband(dataset, 1), 0, 0, 2, 2)) ==
                  ColorTypes.RGB{ImageCore.N0f8}
            @test eltype(AG.imread(dataset, 1, 0, 0, 2, 2)) ==
                  ColorTypes.RGB{ImageCore.N0f8}
            @test eltype(AG.imread(AG.getband(dataset, 1), 1:2, 1:2)) ==
                  ColorTypes.RGB{ImageCore.N0f8}
            @test eltype(AG.imread(dataset, 1, 1:2, 1:2)) ==
                  ColorTypes.RGB{ImageCore.N0f8}
            @test eltype(AG.imread(AG.getband(dataset, 1))) ==
                  ColorTypes.RGB{ImageCore.N0f8}
            @test eltype(AG.imread(dataset, 1)) ==
                  ColorTypes.RGB{ImageCore.N0f8}
            @test eltype(AG.imread(dataset, 1:4, 0, 0, 2, 2)) ==
                  ColorTypes.RGBA{ImageCore.N0f8}
            @test eltype(AG.imread(dataset, 1:4, 1:2, 1:2)) ==
                  ColorTypes.RGBA{ImageCore.N0f8}
            @test eltype(AG.imread(dataset)) == ColorTypes.RGBA{ImageCore.N0f8}
        end

        @testset "Test HSL colors" begin
            AG.setcolorinterp!(AG.getband(dataset, 1), AG.GCI_HueBand)
            AG.setcolorinterp!(AG.getband(dataset, 2), AG.GCI_SaturationBand)
            AG.setcolorinterp!(AG.getband(dataset, 3), AG.GCI_LightnessBand)
            @test_throws ErrorException AG.imread(dataset, 1:3)
        end

        @testset "Test ColorTable colors" begin
            AG.setcolorinterp!(AG.getband(dataset, 1), AG.GCI_PaletteIndex)
            AG.setcolorinterp!(AG.getband(dataset, 2), AG.GCI_PaletteIndex)
            AG.setcolorinterp!(AG.getband(dataset, 3), AG.GCI_PaletteIndex)
            AG.setcolorinterp!(AG.getband(dataset, 4), AG.GCI_PaletteIndex)

            AG.createcolortable(AG.GPI_RGB) do ct
                AG.setcolorentry!(
                    ct,
                    typemax(UInt8),
                    GDAL.GDALColorEntry(0, 0, 0, 0),
                )
                AG.setcolortable!(AG.getband(dataset, 1), ct)
                AG.setcolortable!(AG.getband(dataset, 2), ct)
                AG.setcolortable!(AG.getband(dataset, 3), ct)
                return AG.setcolortable!(AG.getband(dataset, 4), ct)
            end
            @test eltype(AG.imread(dataset)) == ColorTypes.RGBA{ImageCore.N0f8}

            AG.createcolortable(AG.GPI_Gray) do ct
                AG.setcolorentry!(
                    ct,
                    typemax(UInt8),
                    GDAL.GDALColorEntry(0, 0, 0, 0),
                )
                return AG.setcolortable!(AG.getband(dataset, 4), ct)
            end
            @test eltype(AG.imread(dataset, 4)) ==
                  ColorTypes.Gray{ImageCore.N0f8}

            AG.createcolortable(AG.GPI_CMYK) do ct # CMYK not supported yet
                AG.setcolortable!(AG.getband(dataset, 1), ct)
                AG.setcolortable!(AG.getband(dataset, 2), ct)
                return AG.setcolortable!(AG.getband(dataset, 3), ct)
            end
            @test_throws ErrorException AG.imread(dataset, 1:3)
        end
    end
end


text:
module TestUtils

using HDF5
const API = HDF5
#const API = HDF5.API

export test_filter

function test_filter_init(; cd_values = Cuint[], data = ones(UInt8, 1024))
    flags = Cuint(0)
    nbytes = sizeof(data)
    buf_size = Ref(Csize_t(sizeof(data)))
    databuf = Libc.malloc(sizeof(data))
    data = reinterpret(UInt8, data)
    unsafe_copyto!(Ptr{UInt8}(databuf), pointer(data), sizeof(data))
    buf = Ref(Ptr{Cvoid}(databuf))
    return flags, cd_values, nbytes, buf_size, buf
end

function test_filter_compress!(filter_func, flags::Cuint, cd_values::Vector{Cuint}, nbytes::Integer, buf_size::Ref{Csize_t}, buf::Ref{Ptr{Cvoid}})
    nbytes = Csize_t(nbytes)
    cd_nelmts = Csize_t(length(cd_values))
    GC.@preserve flags cd_nelmts cd_values nbytes buf_size buf begin
        ret_code = filter_func(
            flags,
            cd_nelmts,
            pointer(cd_values),
            Csize_t(nbytes),
            Base.unsafe_convert(Ptr{Csize_t}, buf_size),
            Base.unsafe_convert(Ptr{Ptr{Cvoid}}, buf)
        )
        @info "Compression:" ret_code buf_size[]
        if ret_code <= 0
            error("Test compression failed: $ret_code.")
        end
    end
    return ret_code
end

function test_filter_decompress!(filter_func, flags::Cuint, cd_values::Vector{Cuint}, nbytes::Integer, buf_size::Ref{Csize_t}, buf::Ref{Ptr{Cvoid}})
    nbytes = Csize_t(nbytes)
    cd_nelmts = Csize_t(length(cd_values))
    flags |= UInt32(API.H5Z_FLAG_REVERSE)
    GC.@preserve flags cd_nelmts cd_values nbytes buf_size buf begin
        ret_code = filter_func(
            flags,
            cd_nelmts,
            pointer(cd_values),
            Csize_t(nbytes),
            Base.unsafe_convert(Ptr{Csize_t},buf_size),
            Base.unsafe_convert(Ptr{Ptr{Cvoid}}, buf)
        )
        @info "Decompression:" ret_code buf_size[]
    end
    return ret_code
end

function test_filter_cleanup!(buf::Ref{Ptr{Cvoid}})
    Libc.free(buf[])
end

function test_filter(filter_func; cd_values::Vector{Cuint} = Cuint[], data = ones(UInt8, 1024))
    flags, cd_values, nbytes, buf_size, buf = test_filter_init(; cd_values = cd_values, data = data)
    nbytes_compressed, nbytes_decompressed = 0, 0
    try
        nbytes_compressed = test_filter_compress!(filter_func, flags, cd_values, nbytes, buf_size, buf)
        nbytes_decompressed = test_filter_decompress!(filter_func, flags, cd_values, nbytes_compressed, buf_size, buf)
        if nbytes_decompressed > 0
            # ret_code is the number of bytes out
            round_trip_data = unsafe_wrap(Array,Ptr{UInt8}(buf[]), nbytes_decompressed)
            @info "Is the data the same after a roundtrip?" data == round_trip_data
        end
    catch err
        rethrow(err)
    finally
        test_filter_cleanup!(buf)
    end
    @info "Compression Ratio" nbytes_compressed / nbytes_decompressed
    return nbytes_compressed, nbytes_decompressed
end

end