text:
using Logging
using MPI
using JLD2
using BenchmarkTools

using Oceananigans
using Oceananigans.Distributed
using Benchmarks

Logging.global_logger(OceananigansLogger())

MPI.Init()
comm = MPI.COMM_WORLD

local_rank = MPI.Comm_rank(comm)
R = MPI.Comm_size(comm)

Nx = parse(Int, ARGS[1])
Ny = parse(Int, ARGS[2])
Nz = parse(Int, ARGS[3])

@info "Setting up distributed incompressible model with N=($Nx, $Ny, $Nz) grid points on $R rank(s)..."

topo = (Periodic, Periodic, Periodic)
distributed_grid = RegularRectilinearGrid(topology=topo, size=(Nx, Ny, Nz), extent=(1, 1, 1))
arch = MultiCPU(grid=distributed_grid, ranks=(1, R, 1))
model = DistributedIncompressibleModel(architecture=arch, grid=distributed_grid)

@info "Warming up distributed incompressible model..."

time_step!(model, 1) # warmup

@info "Benchmarking distributed incompressible model..."

trial = @benchmark begin
    @sync_gpu time_step!($model, 1)
    MPI.Barrier(comm)
end samples=10

@info "Rank $local_rank is done benchmarking!"

jldopen("strong_scaling_incompressible_model_$(R)_$local_rank.jld2", "w") do file
    file["trial"] = trial
end


text:
## Pure ######################################################################

"""
    Pure(x) :: Distribution{typeof{x}}

Wrap `x` as a distribution always yielding `x`.
`Pure(x)` is equivalent to the implicit distributions
`[x]` and `(x,)`.

# Examples
```jldoctest
julia> rand(Pure(1), 3)
3-element Array{Int64,1}:
 1
 1
 1

julia> rand([1], 3)
3-element Array{Int64,1}:
 1
 1
 1
```

!!! note
    The name `Pure` comes from the similar function of the "applicative"
    Haskell typeclass.
"""
struct Pure{T} <: Distribution{T}
    x::T
end

rand(::AbstractRNG, sp::SamplerTrivial{<:Pure}) = sp[].x

Base.show(io::IO, x::Pure) = print(io, "Pure($(x.x))")

## algebra ###################################################################

struct Op2{T,F,A,B} <: Distribution{T}
    f::F
    a::A
    b::B

    Op2{T}(f::F, a::A, b::B) where {T,F,A,B} = new{T,F,A,B}(f, a, b)
end

Op2(f::F, a::A, b::B) where {F,A,B} =
    Op2{typeof(f(one(gentype(a)), one(gentype(b))))}(f, a, b)

Sampler(RNG::Type{<:AbstractRNG}, x::Op2, n::Repetition) =
    SamplerTag{typeof(x)}((f = x.f,
                           a = Sampler(RNG, x.a, n),
                           b = Sampler(RNG, x.b, n)))

function reset!(sp::SamplerTag{<:Op2}, n...)
    reset!(sp.data.a, n...)
    reset!(sp.data.b, n...)
    sp
end

rand(rng::AbstractRNG, sp::SamplerTag{<:Op2{T}}) where {T} =
    sp.data.f(rand(rng, sp.data.a), rand(rng, sp.data.b))::T

### instances

for op = (:+, :-, :*, :/, :^)
    @eval begin
        (Base.$op)(a::Distribution, b::Distribution) = Op2($op, a,       b)
        (Base.$op)(a,               b::Distribution) = Op2($op, Pure(a), b)
        (Base.$op)(a::Distribution, b              ) = Op2($op, a,       Pure(b))
    end
end

### getindex

Op2(::typeof(getindex), a::A, b::B) where {A,B} =
    Op2{eltype(gentype(a))}(getindex, a, b)

Op2(::typeof(getindex), a::A, b::B) where {A<:Distribution{<:AbstractDict},B} =
    Op2{valtype(gentype(a))}(getindex, a, b)

"""
    getindex(X::Distribution, Y::Distribution) :: Distribution

Return a distribution yielding `x[y]` where `x <- X` and `y <- Y`.

# Examples
```julia
julia> rand(Pure('a':'z')[Uniform(1:3)])
'b': ASCII/Unicode U+0062 (category Ll: Letter, lowercase)
```
"""
Base.getindex(a::Distribution, b::Distribution) = Op2(getindex, a, b)

## Bind ######################################################################

"""
    Bind(f, X)

Generate a random value `x` from distribution `X`, and then
generate a random value from distribution specified by `f(x)`.

!!! note
    This is roughly equivalent to the `>>=` bind monadic operator
    in Haskell, but with the order of arguments reversed.

# Examples

The following expression generates 4 arrays of `Bool` (created via
`Fill(Bool, n)`), whose length `n` is chosen at random in `1:5`:
```julia
julia> rand(Bind(n -> Fill(Bool, n), 1:5), 4)
4-element Array{Array{Bool,1},1}:
 [1, 0, 1]
 [1, 1]
 [1, 0, 0, 1, 0]
 [1]
```
If multiple values must be passed to the function, `Zip` can be used
to combine distributions. For example, the following `Bind` construction
generates arrays of length yielded by `Shuffle(1:4)` (i.e. each array has
a distinct length in `1:4`), and whose elements are generated by a
distribution selected randomly in `[Bool, Normal(), 1:9]`:
```julia
julia> b = Bind(Zip([Bool, Normal(), 1:9], Shuffle(1:4))) do (t, n)
               Fill(t, n)
           end;

julia> rand(b, 4)
4-element Array{Any,1}:
 [-0.2900266833838387]
 Bool[0, 0, 1]
 Bool[0, 1, 1, 0]
 [6, 1]
```

!!! note
    The distribution specified by `f(x)` is not known until
    `rand` is called, which means that the "staged" approach
    enabled by the `Sampler` machinery cannot be fully
    exploited. So while `Bind` is a powerful general tool,
    specialized implementations of certain constructions might
    achieve better performance.

    Moreover, only `X` can have a stateful sampler (i.e. which has a
    non-trivial `reset!` method), the distribution returned by `f` can not:
    indeed, a new sampler is created every time.
    For example, `Bind(n -> Iterate(1:n), 1:9)` always yields `1`.
"""
struct Bind{F,X,T} <: Distribution{T}
    f::F
    x::X

    function Bind{T}(f::F, x) where {T,F}
        x = wrap(x)
        new{F,typeof(x),T}(f, x)
    end
end

function Bind(f::F, x) where F
    rt = Base.return_types(f, (gentype(x),))
    T = length(rt) != 1 ? Any : rt[1]
    Bind{gentype(T)}(f, x)
end

Sampler(::Type{RNG}, b::Bind, n::Repetition) where {RNG<:AbstractRNG} =
    SamplerTag{typeof(b)}((x = sampler(RNG, b.x, n),
                           f = b.f))

reset!(sp::SamplerTag{<:Bind}, n...) = (reset!(sp.data.x); sp)

function rand(rng::AbstractRNG, sp::SamplerTag{<:Bind})
    x = rand(rng, sp.data.x)
    rand(rng, sp.data.f(x))
end

Base.show(io::IO, b::Bind) =
    print(io, "Bind{", gentype(b), "}(", b.f, ", ", b.x, ')')

## Join ######################################################################

"""
    Join(x) :: Distribution

Given a distribution `x` yielding another distribution yielding `y`,
`Join(x)` is a distribution yielding `y`. This is equivalent to
`Bind(identity, x)`. When `x::Distribution{Distribution{T}}`, then
`Join(x)::Distribution{T}`; in other words, `Join` somehow "flattens"
two layers of "`Distribution` wrapping" into one.

# Examples
```julia-repl
julia> rand(Join(Pure(Normal())))
1.3938116112906205

julia> rand(Join([Normal(), Normal(10.0, 1.0)]), 5)
5-element Array{Float64,1}:
  0.8017321165720984
  9.987814723449322
 10.886780925603274
  0.17351815418302222
 11.606547172156468
```

!!! note
    * `Join` can be combined together with `Categorical` to implement
      "mixture models", e.g.
    ```julia
    MixtureModel(components, prior) = Join(Pure(components)[Categorical(prior)])
    ```

    * `Join` can be easily expressed in terms of `Bind`, but is no less
    fundamental, given `Lift`. Indeed, `Bind` could be defined as follows:
    ```julia
    Bind(f, x) = Join(Lift(f, x))
    ```
"""
struct Join{X,T} <: Distribution{T}
    x::X

    function Join(x::X) where X
        U = gentype(x)
        new{X,gentype(U)}(x)
    end
end

Sampler(::Type{RNG}, j::Join, n::Repetition) where {RNG<:AbstractRNG} =
    SamplerTag{typeof(j)}(Sampler(RNG, j.x, n))

rand(rng::AbstractRNG, sp::SamplerTag{<:Join}) =
    rand(rng, rand(rng, sp.data))

Base.show(io::IO, j::Join) = print(io, "Join(", j.x, ')')

## Thunk #####################################################################

"""
    Thunk(f) :: Distribution

Create a distribution whose each generated value is the result of calling
`d = f()` and then returning a random value drawn from distribution `d`.
This is a particular case of [`Bind`](@ref), as `Thunk(f)` is equivalent
to `Bind(_ -> f(), Pure(nothing))`.

# Examples
```julia-repl
julia> rand(let i::Int = 0
                Thunk() do
                   i += 1
                   Uniform(1:i)
                end
            end, 5)
5-element Array{Int,1}:
 1
 1
 3
 4
 2
```
"""
struct Thunk{F,T} <: Distribution{T}
    f::F

    Thunk{T}(f::F) where {T,F} = new{F,T}(f)
end

function Thunk(f::F) where F
    rt = Base.return_types(f, ())
    T = length(rt) != 1 ? Any : rt[1]
    Thunk{gentype(T)}(f)
end

rand(rng::AbstractRNG, sp::SamplerTrivial{<:Thunk}) =
    rand(rng, sp[].f())

Base.show(io::IO, t::Thunk) = print(io, "Thunk{$(gentype(t))}(", t.f, ')')

## Lift ######################################################################

"""
    Lift(f, Xs...) :: Distribution

Create a "lifted" version of `f` to the `Distribution` domain, which, given
yielded values `xs...` from `Xs...`, yields `f(xs...)`.

# Examples
```julia-repl
julia> rand(Lift(+, 0:10:20, Normal()), 5)
5-element Array{Float64,1}:
 -1.8066739854479257
 20.285455058446814
  8.674837595869976
  0.7620819803820099
  9.891731903551195
```

!!! note
    When distributions are seen as iterators, `Lift(f, Xs...)` is the iterator
    resulting from `map`ping `f` onto `Xs...`.
"""
struct Lift{T,F,D} <: Distribution{T}
    f::F
    d::D
end

Lift{T}(f::F, d...) where {T,F} = Lift{T,F,typeof(d)}(f, d)

function Lift(f::F, d...) where {F}
    rt = Base.return_types(f, map(gentype, d))
    T = length(rt) != 1 ? Any : rt[1]
    Lift{T}(f, d...)
end

function pmf(m::Lift)
    probas = Dict{gentype(m),Float64}()
    pmfs = map(pmf, m.d)
    for outcomes in Iterators.product(pmfs...)
        outcomeM = m.f((first(outcome) for outcome in outcomes)...)
        proba = prod(last.(outcomes))
        old = get(probas, outcomeM, 0.0)
        probas[outcomeM] = old + proba
    end
    PMF(m, probas)
end

function Base.show(io::IO, m::Lift)
    print(io, "Lift{$(gentype(m))}(", m.f, ", ")
    join(io, m.d, ", ")
    print(io, ')')
end

### sampling

# Repetition -> Val(1)
rand(rng::AbstractRNG, sp::SamplerTrivial{<:Lift{T}}) where {T} =
    convert(T, sp[].f((rand(rng, d) for d in sp[].d)...))

Sampler(RNG::Type{<:AbstractRNG}, m::Lift, n::Val{Inf}) =
    SamplerTag{typeof(m)}((f = m.f,
                           d = map(x -> Sampler(RNG, x, n), m.d)))

reset!(sp::SamplerTag{<:Lift}, n...) =
    (foreach(s -> reset!(s, n...), sp.data.d); sp)

rand(rng::AbstractRNG, sp::SamplerTag{<:Lift{T}}) where {T} =
    convert(T, sp.data.f((rand(rng, d) for d in sp.data.d)...))

## Keep ######################################################################

"""
    Keep(f, X) :: Distribution

Create a distribution yielding a subset of the values yielded by `X`, keeping
only those for which `f` is `true`.

# Examples
```julia-repl
julia> rand(Keep(iseven, 1:9), 5)
5-element Array{Int64,1}:
 2
 8
 2
 8
 4
```

!!! note
    `Keep` is semantically equivalent to the following construction:
    ```julia
    Keep(f, X) = Bind(X) do x
                     f(x) ? Pure(x) :
                            Keep(f, X)
                 end
    ```
"""
struct Keep{T,F,D<:Distribution{T}} <: Distribution{T}
    f::F
    d::D
end

Keep(f::F, d) where {F} = Keep(f, Uniform(d))

function pmf(k::Keep)
    dpmf = cacheall!(pmf(k.d)).pmf
    filter!(x -> k.f(first(x)), dpmf)
    p0 = sum(values(dpmf))
    replace!(x -> x[1] => x[2]/p0, dpmf)
    PMF(k, dpmf)
end

Base.show(io::IO, k::Keep) = print(io, "Keep(", k.f, ", ", k.d, ')')

### sampling

Sampler(::Type{RNG}, d::Keep, n::Repetition) where {RNG<:AbstractRNG} =
    SamplerTag{typeof(d)}((f = d.f,
                           d = Sampler(RNG, d.d, n)))

reset!(sp::SamplerTag{<:Keep}, n=0) = (reset!(sp.data.d, n); sp)

rand(rng::AbstractRNG, sp::SamplerTag{<:Keep}) =
    while true
        x = rand(rng, sp.data.d)
        sp.data.f(x) && return x
    end

## Map #######################################################################

"""
    Map(f, D...)

Given distributions `D...` yielding collections `d...`, create a distribution
yielding `map(f, d...)`.

!!! note
    `Map(f, D...)` is semantically equivalent to `Lift(map, Pure(f), D...)`.

# Examples
```julia
julia> rand(Map(+, Fill(0:10:20, 4), Fill(Normal(), 4)))
4-element Array{Float64,1}:
 20.51571632364027
 10.458305495441273
  0.24391036203770697
 18.700973042033308
```
"""
struct Map{T,F,D} <: Distribution{T}
    f::F
    d::D
end

Map{T}(f::F, d...) where {T,F} = Map{T,F,typeof(d)}(f, d)

function Map(f::F, d...) where {F}
    rt = Base.return_types(map, (F, map(gentype, d)...))
    T = length(rt) != 1 ? Any : rt[1]
    Map{T}(f, d...)
end

### sampling

rand(rng::AbstractRNG, sp::SamplerTrivial{<:Map{T}}) where {T} =
    convert(T, map(sp[].f, (rand(rng, d) for d in sp[].d)...))

Sampler(RNG::Type{<:AbstractRNG}, m::Map, n::Val{Inf}) =
    SamplerTag{typeof(m)}((f = m.f,
                           d = map(x -> Sampler(RNG, x, n), m.d)))

reset!(sp::SamplerTag{<:Map}, n...) =
    (foreach(s -> reset!(s, n...), sp.data.d); sp)

rand(rng::AbstractRNG, sp::SamplerTag{<:Map{T}}) where {T} =
    convert(T, map(sp.data.f, (rand(rng, d) for d in sp.data.d)...))

### misc

function Base.show(io::IO, m::Map)
    print(io, "Map{", gentype(m), "}(", m.f, ", ")
    join(io, m.d, ", ")
    print(io, ')')
end

## Filter ####################################################################

"""
    Filter(   f, X) :: Distribution{eltype(X)}
    Filter{T}(f, X) :: Distribution{T}

Given distribution `X` yielding collection `x`, create a distribution
yielding `filter(f, x)`.
In the first form, `Filter(f, X)`, the collection `filter(f, x)` is
assumed to be of the same type as `x`. If this is not the case and
it has type `T` instead, use the second form `Filter{T}(f, X)`.

!!! note
    `Filter(f, X)` is semantically equivalent to `Lift(filter, Pure(f), X)`,
    which can be used as an alternative to `Filter{T}(f, X)` when `T` is not
    known.

# Examples
```julia-repl
julia>  rand(Filter(x -> x != 0, Fill(-1:1, 6)))
4-element Array{Int64,1}:
 -1
 -1
  1
 -1

julia> rand(Filter(x -> x != 0, Zip(-1:1, -1:1)))
ERROR: MethodError: Cannot `convert` an object of type
  Tuple{} to an object of type
  Tuple{Int64}
[...]

julia> rand(Filter{Tuple}(x -> x != 0, Zip(-1:1, -1:1)), 3)
3-element Array{Tuple,1}:
 (-1,)
 ()
 (1, -1)

julia> rand(Lift(filter, Pure(x -> x != 0), Zip(-1:1, -1:1)), 3)
3-element Array{Tuple{Vararg{Int64,N} where N},1}:
 (-1, 1)
 ()
 (1, -1)
```
"""
struct Filter{T,F,D} <: Distribution{T}
    f::F
    d::D
end

Filter{T}(f::F, d) where {T,F} = Filter{T,F,typeof(d)}(f, d)

Filter(f::F, d) where {F} = Filter{eltype(d)}(f, d)

### sampling

Sampler(RNG::Type{<:AbstractRNG}, r::Filter, n::Repetition) =
    SamplerTag{typeof(r)}((f = r.f,
                           d = Sampler(RNG, r.d, n)))

reset!(sp::SamplerTag{<:Filter}, n...) = (reset!(sp.data.d, n...); sp)

rand(rng::AbstractRNG, sp::SamplerTag{<:Filter{T}}) where {T} =
    convert(T, filter(sp.data.f, rand(rng, sp.data.d)))

### misc

Base.show(io::IO, f::Filter) = print(io, "Filter(", f.f, ", ", f.d, ')')

## Reduce ####################################################################

struct Reduce{T,F,D} <: Distribution{T}
    f::F
    d::D
end

Reduce{T}(f::F, d) where {T,F} = Reduce{T,F,typeof(d)}(f, d)

# we only support reduce for f(::X, ::X) -> X
# use Lift + Base.reduce for more complicated cases
Reduce(f::F, d) where {F} = Reduce{eltype(gentype(d))}(f, d)

### sampling

rand(rng::AbstractRNG, sp::SamplerTrivial{<:Reduce{T}}) where {T} =
    convert(T, reduce(sp[].f, rand(rng, sp[].d)))

Sampler(RNG::Type{<:AbstractRNG}, r::Reduce, n::Val{Inf}) =
    SamplerTag{typeof(r)}((f = r.f,
                           d = Sampler(RNG, r.d, n)))

reset!(sp::SamplerTag{<:Reduce}, n=0) = (reset!(sp.data.d, n); sp)

rand(rng::AbstractRNG, sp::SamplerTag{<:Reduce{T}}) where {T} =
    convert(T, reduce(sp.data.f, rand(rng, sp.data.d)))

## Counts

"""
    Counts(x) :: Distribution{<:Dict}

Create a distribution yielding a dictionary whose keys are the elements
of the collection yielded by distribution `x`, and whose values are the
number of times each element appeared in the collection.

# Examples
```julia
julia> rand(Counts(Fill(Categorical([1/6, 2/6, 3/6]), 600)))
Dict{Int64,Int64} with 3 entries:
  2 => 193
  3 => 306
  1 => 101
```
"""
struct Counts{T,X} <: Distribution{Dict{T,Int}}
    x::X

    Counts(x::X) where {X} = new{eltype(gentype(x)),X}(x)
end

Sampler(::Type{RNG}, c::Counts, n::Repetition) where {RNG<:AbstractRNG} =
    SamplerTag{typeof(c)}(Sampler(RNG, c.x, n))

reset!(sp::SamplerTag{<:Counts}, n...) = (reset!(sp.data, n...); sp)

function rand(rng::AbstractRNG, sp::SamplerTag{<:Counts{T}}) where T
    dict = Dict{T,Int}()
    for x in rand(rng, sp.data)
        dict[x] = get(dict, x, 0) + 1
    end
    dict
end

"""
    counts(x, [n::Integer])

Equivalent to `rand(Counts(x))`, or to `rand(Counts(Fill(x, n)))`
when `n` is specified.

!!! warning
    Experimental function.
"""
counts(x) = rand(Counts(x))
counts(x, n) = rand(Counts(Fill(x, n)))

## Unique ####################################################################

struct Unique{T,X} <: Distribution{T}
    x::X

    Unique(x::X) where {X} = new{gentype(x),X}(x)
end

Unique(::Type{X}) where {X} = Unique(Uniform(X))

Sampler(RNG::Type{<:AbstractRNG}, u::Unique, n::Val{1}) =
    Sampler(RNG, u.x, n)

Sampler(RNG::Type{<:AbstractRNG}, u::Unique, n::Val{Inf}) =
    SamplerTag{typeof(u)}((x    = Sampler(RNG, u.x, n),
                           seen = Set{gentype(u)}()))

function reset!(sp::SamplerTag{<:Unique}, n=0)
    seen = sp.data.seen
    n > length(seen) && sizehint!(seen, n)
    empty!(seen)
    sp
end

function rand(rng::AbstractRNG, sp::SamplerTag{<:Unique})
    seen = sp.data.seen
    while true
        x = rand(rng, sp.data.x)
        x in seen && continue
        push!(seen, x)
        return x
    end
end

## Iterated ##################################################################

"""
    Iterated{T}(f) :: Distribution{T}

Create a distribution representing the process of iterating `f`
an arbitrary number of times. The function `f` must have two methods,
one taking one `AbstractRNG` argument and the other taking two arguments.
When sampling an `Iterated` object, the first step consists in
calling `dist, st = f(rng)` where `rng` is the RNG which will be used.
`dist` is the distribution to use in the next iteration,
and `st` is an arbitrary object, which
can be any state that `f` needs in subsequent calls. Then, the
following steps are repeated in a loop until `dist === nothing`:
+ `x = rand(rng, dist)`
+ `dist, st = f(x, st)`
When `dist === nothing`, `st` is interpreted as the return value of the
process.

# Examples
```julia-repl
julia> geom(p) = Iterated{Int}(
            # implement the geometric distribution
            begin
                f(_) = Bernoulli(p), 0
                f(v, i) = v ? (nothing, i) : (Bernoulli(p), i+1)
            end);

julia> pmf(rand(geom(0.8), 100))
pmf for [0, 0, 0, 0, 0, 0, 0, 0, 1, 0  …  0, 0, 1, 0, 0, 0, 0, 0, 0, 0] with support of length 5:
  0 => 0.81
  1 => 0.13
  2 => 0.03
  3 => 0.01
  4 => 0.02
```

!!! note
    In the first step, i.e. `dist, st = f(rng)`, `rng` can be used within
    `f` to create a `Sampler` object, which can then be returned as `dist`
    and passed as `st` (or stored within `st`) across calls.
"""
struct Iterated{T,F} <: Distribution{T}
    f::F

    Iterated{T}(f::F) where {T,F} = new{T,F}(f)
end

function rand(rng::AbstractRNG, sp::SamplerTrivial{<:Iterated{T}})::T where T
    it = sp[]
    dist, st = it.f(rng)
    while dist !== nothing
        x = rand(rng, dist)
        dist, st = it.f(x, st)
    end
    st
end

function pmf(it::Iterated{T}; normalized::Bool=true) where T
    itpmf = Dict{T,Float64}()
    dist, st = it.f(Random.GLOBAL_RNG)
    cnt = iterated_pmf!(itpmf, it.f, 1.0, pmf(dist), st)
    PMF(it, itpmf, count=cnt)
end

function iterated_pmf!(itpmf, f, p0, distpmf, st0)
    l = 1
    for (x, p) in distpmf
        p1 = p0 * p
        p1 == 0.0 && continue
        st = deepcopy(st0)
        dist, st = f(x, st)
        if dist === nothing
            old = get(itpmf, st, 0.0)
            itpmf[st] = old + p1
        else
            l = lcm(l, iterated_pmf!(itpmf, f, p1, pmf(dist), st))
        end
    end
    l * Int(-distpmf.count)
end

ext: jl


text:
# https://zhuanlan.zhihu.com/p/39253127

function walk(L=50)
    t = 2
    x = 2
    while x != 1
        if rand() < 0.5
            x += 1
        else
            x -= 1
        end
        if x > L
            x = L
        end
        t += 1
    end
    return t
end

function run(L, N)
    times = Int64[]
    for i in 1:N
        push!(times, walk(L))
    end
    return times
end

function repeat(L, N, M)
    [mean(run(L, N)) for i in 1:M]
end

data = repeat(50, 1000, 10000);
using Plots
mean(data)
histogram(data)

max_stars_repo_name: unrealTOM/MC

max_issues_repo_name: unrealTOM/MC

max_forks_repo_name: unrealTOM/MC


text:
_mul_plotgrid(_, args) = grid(first(args))
_mul_plotgrid(::Tuple{Any,PaddedLayout}, (P,c)) = plotgrid(P[:,colsupport(c)])

function _plotgrid(lay::ExpansionLayout, P)
    args = arguments(lay,P)
    _mul_plotgrid(map(MemoryLayout,args), args)
end

_plotgrid(_, P) = grid(P)

_plotgrid(::WeightedBasisLayouts, wP) = plotgrid(unweighted(wP))
_plotgrid(::MappedBasisLayout, P) = invmap(parentindices(P)[1])[plotgrid(demap(P))]

plotgrid(g) = _plotgrid(MemoryLayout(g), g)

_split_svec(x) = (x,)
_split_svec(x::AbstractArray{<:StaticVector{2}}) = (map(first,x), map(last,x))

@recipe function f(g::AbstractQuasiVector)
    x = plotgrid(g)
    tuple(_split_svec(x)...,g[x])
end

@recipe function f(g::AbstractQuasiMatrix)
    x = plotgrid(g)
    tuple(_split_svec(x)...,g[x,:])
end


text:
# Attempt to find all zeros in an interval (a,b)

# Algorithm due to @djsegal in https://github.com/JuliaMath/Roots.jl/pull/113

# A naive approach to find zeros: split (a,b) by n points, look into each for a zero
# * k is oversampling rate for bisection. (It is relatively cheap to check for a bracket so we
#   oversample our intervals looking for brackets
# * assumes f(a) *not* a zero
function _fz(f, a, b, no_pts, k=4)
    zs = Real[]
    _fz!(zs, f, a, b, no_pts, k)
    zs
end

function _fz!(zs, f, a::T, b, no_pts, k=4) where {T}

    pts = range(a, stop=b, length=(no_pts-1)*k+1)
    n::Int = length(pts)

    fs = f.(pts)
    sfs = sign.(fs)

    u::T = first(pts)  # keep track of bigger interval
    found_bisection_zero = false

    for (i,x) in enumerate(pts[1:end])
        q,r = divrem(i-1, k)

        if i > 1 && iszero(r)
            v::T = x
            if !found_bisection_zero
                try
                    p1::T = identify_starting_point(u, v, sfs[(i-k):i])
                    rt::T = dfree(f, p1)
                    if !isnan(rt) && u < rt <= v
                        push!(zs, rt)
                    end
                catch err
                end
            end
            u = v
            found_bisection_zero = false
        end

        if i < n
            if iszero(fs[i+1])
                found_bisection_zero = true # kinda
                push!(zs, pts[i+1])
            elseif  sfs[i] * sfs[i+1] < 0
                found_bisection_zero = true
                rt = bisection(f, x, pts[i+1])
                push!(zs, rt)
            end
        end
    end

    sort!(zs)
end

# the algorithm first scans for zeros using the naive approach, then
# splits (a,b) by these zeros. This struct holds the subintervals
struct Interval{T}
a::T
b::T
depth::Int
end

Base.show(io::IO, alpha::Interval) = print(io, "($(alpha.a), $(alpha.b))")

# check if f(a) is non zero using tolerances max(atol, eps()), rtol
function _non_zero(fa, a::T, atol, rtol) where {T}
    abs(fa) >= max(atol, abs(a)*rtol*oneunit(fa)/oneunit(a), oneunit(fa)*eps(T))
end

# After splitting by zeros we have intervals (zm, zn) this is used to shrink
# to (zm+, zn-) where both are non-zeros, as defined above
function find_non_zero(f, a::T, barrier, xatol, xrtol, atol, rtol) where {T}
    nan = (0*a)/(0*a) # try to get typed NaN
    xtol = max(xatol, abs(a) * xrtol, oneunit(a) * eps(T))
    sgn = barrier > a ? 1 : -1
    ctr = 0
    x = a + 2^ctr*sgn*xtol
    while !_non_zero(f(x), x, atol, rtol)
        ctr += 1
        x += 2^ctr*sgn*xtol
        ((sgn > 0 && x > barrier) || (sgn < 0 && x < barrier)) && return nan
        ctr > 100 && return nan
    end

    x
end

# split a < z1 < z2 < ... < zn < b into intervals (a+,z1-), (z1+, z2-), ...
# where possible; push! onto ints
function make_intervals!(ints, f, a, b, zs, depth, xatol, xrtol, atol, rtol)
    pts = vcat(a, zs, b)

    for (u,v) in zip(pts[1:end-1], pts[2:end])
        ur = find_non_zero(f, u, v, xatol, xrtol, atol, rtol)
        isnan(ur) && continue

        vl = find_non_zero(f, v, u, xatol, xrtol, atol, rtol)
        isnan(vl) && continue

        push!(ints, Interval(ur, vl, depth))
    end
end

# adjust what we mean by x1 ~ x2 for purposes of adding a new zero
function approx_close(z1, z2, xatol, xrtol)
    tol = max(sqrt(xatol), max(abs(z1), abs(z2)) * sqrt(xrtol))
    abs(z1 - z2) < tol
end

# is proposed not near xs? (false and we add proposed)
function not_near(proposed, xs, xatol, xrtol)
    n = length(xs)
    n <= 1 && return true
    ind = n + 1
    for (i,rt) in enumerate(xs)
        if proposed < rt
            ind = i
            break
        end
    end
    if ind > 1 # xs[ind-1] <= propose < xs[ind]
        rt = xs[ind-1]
        approx_close(proposed, rt, xatol, xrtol) && return false
    end
    if ind <= n # value to right
        rt = xs[ind]
        approx_close(proposed, rt, xatol, xrtol) && return false
    end
    return true
end

"""
    find_zeros(f, a, b; [no_pts=12, k=8, naive=false, xatol, xrtol, atol, rtol])

Search for zeros of `f` in the interval `[a,b]`.

# Examples

```jldoctest find_zeros
julia> using Roots

julia> find_zeros(x -> exp(x) - x^4, -5, 20)        # a few well-spaced zeros
3-element Vector{Float64}:
 -0.8155534188089606
  1.4296118247255556
  8.613169456441398

julia> find_zeros(x -> sin(x^2) + cos(x)^2, 0, 2pi)  # many zeros
12-element Vector{Float64}:
 1.78518032659534
 2.391345462376604
 3.2852368649448853
 3.3625557095737544
 4.016412952618305
 4.325091924521049
 4.68952781386834
 5.00494459113514
 5.35145266881871
 5.552319796014526
 5.974560835055425
 6.039177477770888

julia> find_zeros(x -> cos(x) + cos(2x), 0, 4pi)    # mix of simple, non-simple zeros
6-element Vector{Float64}:
  1.0471975511965976
  3.141592653589943
  5.235987755982988
  7.330382858376184
  9.424777960769228
 11.519173063162574

julia> f(x) = (x-0.5) * (x-0.5001) * (x-1)          # nearby zeros
f (generic function with 1 method)

julia> find_zeros(f, 0, 2)
3-element Vector{Float64}:
 0.5
 0.5001
 1.0

julia> f(x) = (x-0.5) * (x-0.5001) * (x-4) * (x-4.001) * (x-4.2)
f (generic function with 1 method)

julia> find_zeros(f, 0, 10)
3-element Vector{Float64}:
 0.5
 0.5001
 4.2

julia> f(x) = (x-0.5)^2 * (x-0.5001)^3 * (x-4) * (x-4.001) * (x-4.2)^2  # hard to identify
f (generic function with 1 method)

julia> find_zeros(f, 0, 10, no_pts=21)                # too hard for default
5-element Vector{Float64}:
 0.49999999999999994
 0.5001
 4.0
 4.001
 4.200000000000001
```

!!! note
    There are two cases where the number of zeros may be underreported:
    * if the initial interval, [a,b], is too wide
    * if there are zeros  that are very nearby

----

The basic algorithm checks for zeros among the endpoints, and then
divides the interval `(a,b)` into `no_pts-1` subintervals and then
proceeds to look for zeros through bisection or a derivative-free
method.  As checking for a bracketing interval is relatively cheap and
bisection is guaranteed to converge, each interval has `k` pairs of
intermediate points checked for a bracket.

If any zeros are found, the algorithm uses these to partition `(a,b)`
into subintervals. Each subinterval is shrunk so that the endpoints
are not zeros and the process is repeated on the subinterval. If the
initial interval is too large, then the naive scanning for zeros may
be fruitless and no zeros will be reported. If there are nearby zeros,
the shrinking of the interval may jump over them, though as seen in
the examples, nearby roots can be identified correctly, though for
really nearby points, or very flat functions, it may help to increase
`no_pts`.

The tolerances are used to shrink the intervals, but not to find zeros
within a search. For searches, bisection is guaranteed to converge
with no specified tolerance. For the derivative free search, a
modification of the `Order0` method is used, which at worst case
compares `|f(x)| <= 8*eps(x)` to identify a zero. The algorithm might
identify more than one value for a zero, due to floating point
approximations. If a potential pair of zeros satisfy
`isapprox(a,b,atol=sqrt(xatol), rtol=sqrt(xrtol))` then they are
consolidated.

The algorithm can make many function calls. When zeros are found in an
interval, the naive search is carried out on each subinterval. To cut
down on function calls, though with some increased chance of missing
some zeros, the adaptive nature can be skipped with the argument
`naive=true` or the number of points stepped down.

The algorithm is derived from one in a
[PR](https://github.com/JuliaMath/Roots.jl/pull/113) by @djsegal.

!!! Note
The `IntervalRootFinding` package provides a rigorous alternative to this heuristic one.
That package use interval arithmetic, so can compute bounds on the size of the image of
an interval under a function `f`. If this image includes `0`, then it can look for the zero.
Bisection, on the other hand, only will look for a zero if the two endpoints have different signs,
a much more rigid condition for a potential zero.

For example, this function (due to `@truculentmath`) is particularly tricky, as it is positive at every floating point number, but has two zeros (the vertical asymptote at `15//11` is negative between adjacent floating point values):

```jldoctest interval_root_finding
julia> using IntervalArithmetic, IntervalRootFinding, Roots

julia> g(x) = x^2 + 1 +log(abs( 11*x-15 ))/99
g (generic function with 1 method)

julia> find_zeros(g, -3, 3)
Float64[]

julia> IntervalRootFinding.roots(g, -3..3, IntervalRootFinding.Bisection)
1-element Vector{Root{Interval{Float64}}}:
 Root([1.36363, 1.36364], :unknown)
```

A less extreme usage might be the following, where `unique` indicates Bisection could be useful and indeed `find_zeros` will identify these values:

```jldoctest interval_root_finding
julia> g(x) = exp(x) - x^5
g (generic function with 1 method)

julia> rts = IntervalRootFinding.roots(g, -20..20)
2-element Vector{Root{Interval{Float64}}}:
 Root([12.7132, 12.7133], :unique)
 Root([1.29585, 1.29586], :unique)

julia> find_zeros(g, -20, 20)
2-element Vector{Float64}:
  1.2958555090953687
 12.713206788867632
```

"""
function find_zeros(f, a, b; no_pts = 12, k=8,
                    naive = false,
                    kwargs...
                    )

    a0, b0 = promote(float(a), float(b))
    a0 = isinf(a0) ? nextfloat(a0) : a0
    b0 = isinf(b0) ? prevfloat(b0) : b0

    # set tolerances if not specified
    fa0, fb0 = f(a0), f(b0)
    d = Dict(kwargs)
    T, S = eltype(a0), eltype(fa0)
    xatol::T = get(d, :xatol, eps(one(T))^(4//5) * oneunit(T))
    xrtol = get(d, :xrtol, eps(one(T)) * one(T))
    atol::S  = get(d, :atol,  eps(float(S)) * oneunit(S))
    rtol  = get(d, :rtol,  eps(float(S)) * one(S))

    zs = T[]  # collect zeros

    # check endpoints for exact zeros, then narrow
    abs(fa0) * oneunit(T) / oneunit(S) <= 8 * eps(a0) && push!(zs, a0)
    abs(fb0) * oneunit(T) / oneunit(S) <= 8 * eps(b0) && push!(zs, b0)
    a0 = find_non_zero(f, a0, b0, xatol, xrtol, atol, rtol)
    b0 = find_non_zero(f, b0, a0, xatol, xrtol, atol, rtol)

    _fz!(zs, f, a0, b0, no_pts,k)  # initial zeros

    ints = Interval{T}[] # collect subintervals
    !naive && !isempty(zs) &&  make_intervals!(ints, f, a0, b0, zs, 1, xatol, xrtol, atol, rtol)

    nzs = T[]
    cnt = 0

    while !naive && !isempty(ints)
        cnt += 1
        i = pop!(ints)
        # this is fussy. Ideally, we would want to explore intervals
        # with fewer points (each interval is already probed for k=4
        # bisections) but how many fewer?  We already had ~ (i.b-i.a) *
        # no_pts / (b-a) points. Would want an increased density so n
        # > (i.b - i.a) / (b - a) * no_pts but this doesn't perform as
        # well as we might expect

        # sub_no_pts = ceil(Int, (i.b - i.a) / (b-a) * no_pts * 2^(i.depth))
        #sub_no_pts <= 2 && continue  # stop on depth, always divide if roots
        #sub_no_pts = max(3, floor(Int, no_pts  / (2.0)^(i.depth)))
        sub_no_pts = floor(Int, no_pts  / (2.0)^(i.depth))

        empty!(nzs)
        if sub_no_pts >= 2
            _fz!(nzs, f, i.a, i.b, sub_no_pts, k)
        end

        if !isempty(nzs)
            azs = filter(rt->not_near(rt, zs, xatol, xrtol), nzs) # trim out nearby roots
            length(azs) == 0 && continue
            append!(zs, azs)
            sort!(zs)
            i.depth > 4 && continue
            make_intervals!(ints, f, i.a, i.b, azs, i.depth+1, xatol, xrtol, atol, rtol)
        end
    end

    length(zs) <= 1 && return zs
    sort!(zs)

    # may identify same zero with nearby values along the way
    # this trims out with a relaxed tolerance on how close
    # nearby roots can be. Default is epsilon^(2/5) ≈ 5e-7
    inds = Int[1]
    z1 = first(zs)
    for i in 2:length(zs)
        z2 = zs[i]
        if !approx_close(z1, z2, xatol, xrtol)
            push!(inds, i)
        end
        z1=z2
    end

    return zs[inds]

end


text:
const l_exp = 1:4
n_exp = length(l_exp)

l_train = []
l_val = []
for i = 1:n_exp
    j = l_exp[i]
    if !(j in [3])
        push!(l_train, i)
    else
        push!(l_val, i)
    end
end

l_exp_data = [];
l_exp_info = zeros(Float64, length(l_exp), 3);
l_exp_info[:, 2] = [5.0, 10.0, 20.0, 30.0];
l_exp_info[:, 3] .= -100.0;

l_fname = ["1.Beechwood-05.xlsx", "2.Beechwood-10.xlsx",
           "3.Beechwood-20.xlsx", "4.Beechwood-30.xlsx"]

for (i_exp, value) in enumerate(l_exp)
    fname = l_fname[i_exp]
    xf = XLSX.readxlsx("./results/$fuel_name/exp_data/$fname")
    exp_data = Float64.(xf[1][:])
    index = indexin(unique(exp_data[:, 1]), exp_data[:, 1])
    exp_data = exp_data[index, :]
    exp_data[:, 1] .+= 273.15

    l_exp_info[i_exp, 1] = exp_data[1, 1] # initial temperature, K

    ind = findmax(exp_data[:, 1])[2]
    exp_data_ = zeros(ind, 3)
    exp_data_[:, 1] .= (exp_data[1:ind, 1] .- exp_data[1, 1]) ./ l_exp_info[i_exp, 2] .* 60.0
    exp_data_[:, 2] .= exp_data[1:ind, 1]
    exp_data_[:, 3] .= exp_data[1:ind, 2] ./ exp_data[1, 2]

    push!(l_exp_data, exp_data_)
end

max_stars_repo_name: DENG-MIT/CRNN-Pyrolysis

max_issues_repo_name: DENG-MIT/CRNN-Pyrolysis

max_forks_repo_name: DENG-MIT/CRNN-Pyrolysis


text:
#module Dist
    import Distances: evaluate, SemiMetric
    export evaluate, SemiMetric  # reexporting

    include("bits.jl")
    include("sets.jl")
    include("strings.jl")
    include("vectors.jl")
    include("cos.jl")
    include("cloud.jl")
#end

max_stars_repo_name: sadit/NNS

max_issues_repo_name: sadit/NNS

max_forks_repo_name: sadit/NNS


text:
"""
Generic done method, some steppers may implement their own versions.
"""

function done(s::Solver, state::AbstractState)
    if state.step.t >= s.options.tstop
        return true
    elseif state.dt < s.options.minstep
        warn("minstep reached.")
        return true
    elseif state.iters >= s.options.maxiters
        warn("Maximum number of iterations ($(Int(s.options.maxiters))) reached, consider setting a larger maxiter.")
        return true
    end
    return false
end


text:
module Normalization

using Statistics
using JuliennedArrays
import LinearAlgebra:   normalize,
                        normalize!

export  fit,
        fit!,
        normalize!,
        normalize,
        denormalize!,
        denormalize,
        nansafe,
        ZScore,
        RobustZScore,
        Sigmoid,
        RobustSigmoid,
        MinMax

abstract type AbstractNormalization end
(𝒯::Type{<:AbstractNormalization})(dims) = 𝒯(;dims)
function (𝒯::Type{<:AbstractNormalization})(dims, p)
    isnothing(p) || (all(x->x==p[1], length.(p)) && error("Inconsistent parameter dimensions"))
    𝒯(;dims, p)
end
(T::AbstractNormalization)(dims) = dims == () || (T.dims = dims)

macro _Normalization(name, 𝑝, 𝑓, 𝑓⁻¹)
    :(mutable struct $(esc(name)) <: AbstractNormalization
        dims
        p::Union{Nothing, NTuple{length($𝑝), AbstractArray}}
        𝑝::NTuple{length($𝑝), Function}
        𝑓::Function
        𝑓⁻¹::Function
     end;
     ($(esc(name)))(; dims = nothing,
                         p = nothing,
                         𝑝 = $𝑝,
                         𝑓 = $𝑓,
                         𝑓⁻¹ = $𝑓⁻¹) = $(esc(name))(dims, p, 𝑝, 𝑓, 𝑓⁻¹)
     )
end

# * Common normalizations
@_Normalization ZScore (mean, std)         (x, 𝜇, 𝜎) -> x .= (x .- 𝜇)./𝜎  #=
                                        =# (y, 𝜇, 𝜎) -> y .= y.*𝜎 .+ 𝜇
@_Normalization Sigmoid (mean, std)        (x, 𝜇, 𝜎)->x.=1.0./(1 .+exp.(.-(x.-𝜇)./𝜎)) #=
                                        =# (y, 𝜇, 𝜎) -> y .= .-𝜎.*log.(1.0./y .- 1) .+ 𝜇
@_Normalization MinMax (minimum, maximum)  (x, l, u) -> x .= (x.-l)./(u-l) #=
                                        =# (y, l, u) -> y .= (u-l).*y .+ l
@_Normalization Center (mean,)             (x, 𝜇) -> x .= x .- 𝜇     (y, 𝜇) -> y .= y .+ 𝜇
@_Normalization RobustCenter (median,)     Centre().𝑓   Centre().𝑓⁻¹

# * Robust versions of typical 2-parameter normalizations
common_norms = [:ZScore, :Sigmoid,]
_iqr = x -> (quantile(x[:], 0.75) - quantile(x[:], 0.25))/1.35 # ? Divide by 1.35 so that std(x) ≈ _iqr(x) when x contains normally distributed values
_robustNorm(N::Symbol; name="Robust"*string(N)|>Symbol) = eval(:(@_Normalization $name (median, _iqr) ($N)().𝑓 ($N)().𝑓⁻¹))
_robustNorm.(common_norms)

# * Mixed versions of typical 2-parameter normalizations
mixedcenter(x) = (_iqr(x) == 0) ? mean(x) : median(x)
mixedscale(x) = (𝜎 = _iqr(x); 𝜎 == 0 ? std(x) : 𝜎)
_mixedNorm(N::Symbol; name="Mixed"*string(N)|>Symbol) = eval(:(@_Normalization $name (mixedcenter, mixedscale) ($N)().𝑓 ($N)().𝑓⁻¹))
_mixedNorm.(common_norms)

# * NaN-safe versions
_nansafe(p) = x -> p(filter(!isnan, x))
nansafe!(T::AbstractNormalization) = (T.𝑝=_nansafe.(T.𝑝); ())
nansafe(T::AbstractNormalization) = (N = deepcopy(T); nansafe!(N); N)
nansafe(𝒯::Type{<:AbstractNormalization}; dims=nothing) = dims |> 𝒯 |> nansafe

Base.reshape(x::Number, dims...) = reshape([x], dims...)
function fit!(T::AbstractNormalization, X::AbstractArray; dims=())
    T(dims)
    dims = isnothing(T.dims) ? (1:ndims(X)) : T.dims
    psz = size(X) |> collect
    psz[[dims...]] .= 1
    T.p = reshape.(map.(T.𝑝, (Slices(X, dims...),)), psz...)
end
fit(T::AbstractNormalization, X::AbstractArray; kw...)=(T=deepcopy(T); fit!(T, X; kw...); T)
fit(𝒯::Type{<:AbstractNormalization}, X::AbstractArray; dims=nothing) = (T = 𝒯(dims); fit!(T, X); T)

function normalize!(X::AbstractArray, T::AbstractNormalization)
    isnothing(T.p) && fit!(T, X)
    mapdims!(T.𝑓, X, T.p...; T.dims)
end
normalize!(X::AbstractArray, 𝒯::Type{<:AbstractNormalization}; dims=nothing) = normalize!(X, fit(𝒯, X; dims))
normalize(X::AbstractArray, T::Union{AbstractNormalization, Type{<:AbstractNormalization}}; kwargs...) = (Y=copy(X); normalize!(Y, T; kwargs...); Y)

function denormalize!(X::AbstractArray, T::AbstractNormalization)
    isnothing(T.p) && error("Cannot denormalize with an unfit normalization")
    mapdims!(T.𝑓⁻¹, X, T.p...; T.dims)
end
denormalize(X::AbstractArray, args...) = (Y=copy(X); denormalize!(Y, args...); Y)

"""
Map the function `f` over the `dims` of all of the arguments. `f` should accept the same number of arguments as there are variables in `x...`. The first element of `x` is the considered as the reference array, and all other arguments must have sizes consistent with the reference array, or equal to 1.
"""
function mapdims!(f, x...; dims)
    n = ndims(x[1])
    isnothing(dims) && (dims = 1:n)
    dims = sort([dims...])
    @assert max(dims...) <= n
    @assert unique(dims) == dims
    length(dims) == n && return f(x...) # Shortcut for global normalisation
    negdims = Base._negdims(n, dims)
    @assert all(all(size.(x[2:end], i) .== 1) for i ∈ dims)
    @assert all(all(size(x[1], i) .== size.(x, i)) for i ∈ negdims)
    idxs = Base.compute_itspace(x[1], (negdims...,)|>Val)
    Threads.@threads for i ∈ idxs # map(f!, Slices.(x, dims...)...)
        selectslice = x -> view(x, i...)
        f(selectslice.(x)...)
    end
end

end


text:
##################################################################
# Source code taken from: https://github.com/jbytecode/MCCGA
# Credits to corresponding author.
##################################################################

function initialprobs(
    lower::Array{T,1},
    upper::Array{T,1};
    maxsamples = 10000,
) where {T<:Number}
    p = length(lower)
    @assert p == length(upper)

    bitlen = p * 32
    mat = Array{Int8,2}(undef, maxsamples, bitlen)
    for tries = 1:maxsamples
        randvalues = map((x, y) -> x + rand() * (y - x), lower, upper)
        mat[tries, :] .= bits(randvalues)
    end
    return colmeans(mat)
end

function sample(probs::Array{T,1}) where {T<:Number}
    return map(x -> if rand() < x
        1
    else
        0
    end, probs)
end

function colmeans(mat::Array{T,2})::Array{Float64,1} where {T<:Number}
    _, p = size(mat)
    means = Array{Float64,1}(undef, p)
    for i = 1:p
        means[i] = mean(mat[:, i])
    end
    return means
end

"""
        bits(value)

Returns bits of the value. The value is supposed to be 32-bits integer or float.
If the value has larger bits, it is converted to its 32-bits counterpart.

# Arguments:
 - `value::Number`: An integer or float

# Output
- `::Array{Int8, 1}`: Array of bits of size 32

# Example
```julia
julia> bits(3.14159265)
32-element Vector{Int8}:
 0
 1
 0
 0
 0
 0
 0
 0
 0
 1
 0
 ⋮
 1
 1
 1
 1
 1
 0
 1
 1
 0
 1
 1
````

# References
Satman, Mehmet Hakan, and Emre Akadal. "Machine Coded Compact Genetic Algorithms
for Real Parameter Optimization Problems." Alphanumeric Journal 8.1 (2020): 43-58.

Satman, Mehmet Hakan, and Emre Akadal. "Makine Kodlu Hibrit Kompakt Genetik Algoritmalar
Optimizasyon Yöntemi", TR Patent 2018-GE-510,239
"""
function bits(f::T)::Array{Int8,1} where {T<:Number}
    newf = Float32(f)
    casted = reinterpret(Int32, newf)
    return reverse(map(i-> casted >> i & 1, 0:31))
end

"""
        bits(values)

Returns bits of the values. The values is supposed to be an array of 32-bits integer or float.
If the values have larger bits, they are converted to their 32-bits counterpart.

# Arguments:
 - `values::Array{Number, 1}`: An array of integers or floats

# Output
- `::Array{Int8, 1}`: Array of bits of size 32 for each item in the values.

# Example
```
julia> mybits = bits([3.14159265, exp(1.0)]);

julia> println(mybits)
Int8[0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0]
```

# References
Satman, Mehmet Hakan, and Emre Akadal. "Machine Coded Compact Genetic Algorithms
for Real Parameter Optimization Problems." Alphanumeric Journal 8.1 (2020): 43-58.

Satman, Mehmet Hakan, and Emre Akadal. "Makine Kodlu Hibrit Kompakt Genetik Algoritmalar
Optimizasyon Yöntemi", TR Patent 2018-GE-510,239
"""
function bits(fs::Array{T,1})::Array{Int8,1} where {T<:Number}
    return map(x -> bits(x), fs) |> x -> vcat(x...)
end

"""
        floats(bitarray)

Extract float values using the bits in bitarray.

# Arguments:
 - `bitarray::Array{Int, 1}`: An array of integers that are either 0s or 1s.

# Output
- `::Array{Float64, 1}`: Float values that are constructed using the bitarray.

# Example

```

julia> mybits = [0,0,1,1,1,1,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];

julia> floats(mybits)
1-element Vector{Float64}:
 0.25

 julia> mybits = bits([3.14, 2.71]);

julia> floats(mybits)
2-element Vector{Float64}:
 3.140000104904175
 2.7100000381469727

```

# References
Satman, Mehmet Hakan, and Emre Akadal. "Machine Coded Compact Genetic Algorithms
for Real Parameter Optimization Problems." Alphanumeric Journal 8.1 (2020): 43-58.

Satman, Mehmet Hakan, and Emre Akadal. "Makine Kodlu Hibrit Kompakt Genetik Algoritmalar
Optimizasyon Yöntemi", TR Patent 2018-GE-510,239
"""
function floats(bitarray::Array{T,1})::Array{Float64,1} where {T<:Integer}
    numbits = length(bitarray)
    values = Array{Float64,1}(undef, 0)
    starting = 1
    while true
        ending = starting + 32 - 1
        if ending > numbits
            break
        end
        bitpart = bitarray[starting:ending]
        bitsasstring = join(map(x -> string(x), bitpart))
        numasint = parse(UInt32, bitsasstring, base = 2)
        numasfloat = reinterpret(Float32, numasint)
        push!(values, Float64(numasfloat))
        starting = ending + 1
    end
    return values
end


text:
function calculate_constraints(state_array, parameters, constraint_index, syn_data,time_index,tstep)
  alpha = 1.1 # looseness/tightness of constraints
  beta = 10
  n_species = length(state_array[:,1])
  syn_upper = syn_data["upper"]
  syn_lower = syn_data["lower"]
  syn_mean = syn_data["mean"]
  state_array = max(0,state_array)

  species_drift_amplitude = data_dictionary["species_drift_amplitude"]
  species_drift_decay_rate = data_dictionary["species_drift_decay_rate"]
  species_drift = species_drift_amplitude*exp(-time_index*tstep*species_drift_decay_rate)

  lower_bound = 0
  upper_bound = 100000000
  species_lower_bound = max(min(state_array*(-species_drift)/tstep,-0.00000000001/tstep),(lower_bound-state_array)/tstep);
  species_upper_bound = min(max(state_array*(species_drift)/tstep,0.00000000001/tstep),(upper_bound-state_array)/tstep);

  # constrain the species to synthetic data
  for c_index in constraint_index
      species_lower_bound[c_index] = (syn_lower[c_index,time_index+1]-state_array[c_index]) / tstep
      species_upper_bound[c_index] = (syn_upper[c_index,time_index+1]-state_array[c_index]) / tstep
  end
  # oxygen
  species_lower_bound[101] = -10000
  species_upper_bound[101] = 10000
  #h2o
  species_lower_bound[72] = -10000
  species_upper_bound[72] = 10000

  species_constraints = [species_lower_bound species_upper_bound]
  #flux_constraints = [flux_lower_bound flux_upper_bound]

  return species_constraints#, flux_constraints
end


text:
_lowest_eltype(nt::NamedTuple) = _lowest_eltype(values(nt))
_lowest_eltype(nt::Tuple) = _lowest_eltype.(nt)
_lowest_eltype(nt::T) where {T <: Real} = T
_lowest_eltype(nt::AbstractArray{T}) where {T <: Real} = T

_promote_type(nt::NamedTuple) = _promote_type(_lowest_eltype(nt)...)
_promote_type(::Type{T}, t::Tuple) where {T <: Real} = promote_type(T, t...)
_promote_type(t::Tuple, ::Type{T}) where {T <: Real} = promote_type(T, t...)
_promote_type(::Type{X}, ::Type{Y}, args... ) where {X <: Real, Y <: Real} =
    _promote_type(promote_type(X, Y), args...)

@inline _par_in_input_form(p::RT) where {RT} = p
@inline _par_in_input_form(p::T) where {T <: ShapedAsNT} = p[]

function _remove_const_shapes_from_nt(nt::NamedTuple, shape::ValueShapes.AbstractValueShape)
    ks = keys(nt)
    ds = filter(iks -> typeof(shape[iks]) <: ValueAccessor{<:ConstValueShape}, eachindex(ks))
    return !isempty(ds) ? delete(nt, ks[ds]) : nt
end

function subhist(h::Histogram{<:Any, 1}, r::Tuple{<:Real,<:Real})
    first_bin, last_bin = (StatsBase.binindex(h, r[1]), StatsBase.binindex(h, r[2]))
    if first_bin < 1 first_bin = 1 end
    if (last_bin > length(h.weights)) last_bin = length(h.weights) end
    Histogram(h.edges[1][first_bin:last_bin+1], h.weights[first_bin:last_bin])
end
subhist(h, i::Interval) = subhist(h, (i.left, i.right))


text:
"""
Macro that inputs an expression corresponding to a reaction network and outputs
a `ReactionNetwork` that can be used as input to generation of ODE, SDE, and
Jump problems.

Most arrows accepted (both right, left, and bi-drectional arrows). Note that
while --> is a correct arrow, neither <-- nor <--> works. Using non-filled
arrows (⇐, ⟽, ⇒, ⟾, ⇔, ⟺) will disable mass kinetics and let you cutomize
reaction rates yourself. Use 0 or ∅ for degradation/creation to/from nothing.
Example systems: ### Basic Usage ### rn = @reaction_network begin           #
Creates a ReactionSystem. 2.0, X + Y --> XY                  # This will have
reaction rate corresponding to 2.0*[X][Y] 2.0, XY ← X + Y                    #
Identical to 2.0, X + Y --> XY end

    ### Manipulating Reaction Rates ###
    rn = @reaction_network begin
        2.0, X + Y ⟾ XY                   # Ignores mass kinetics. This will have reaction rate corresponding to 2.0.
        2.0X, X + Y --> XY                 # Reaction rate needs not be constant. This will have reaction rate corresponding to 2.0*[X]*[X]*[Y].
        XY+log(X)^2, X + Y --> XY          # Reaction rate accepts quite complicated expressions.
        hill(XY,2,2,2), X + Y --> XY       # Reaction inis activated by XY according to a hill function. hill(x,v,K,N).
        mm(XY,2,2), X + Y --> XY           # Reaction inis activated by XY according to a michaelis menten function. mm(x,v,K).
    end

    ### Multiple Reactions on a Single Line ###
    rn = @reaction_network begin
        (2.0,1.0), X + Y ↔ XY              # Identical to reactions (2.0, X + Y --> XY) and (1.0, XY --> X + Y).
        2.0, (X,Y) --> 0                   # This corresponds to both X and Y degrading at rate 2.0.
        (2.0, 1.0), (X,Y) --> 0            # This corresponds to X and Y degrading at rates 2.0 and 1.0, respectively.
        2.0, (X1,Y1) --> (X2,Y2)           # X1 and Y1 becomes X2 and Y2, respectively, at rate 2.0.
    end

    ### Adding Parameters ###
    kB = 2.0; kD = 1.0
    p = [kB, kD]
    p = []
    rn = @reaction_network begin
        (kB, kD), X + Y ↔ XY               # Lets you define parameters outside on network. Parameters can be changed without recalling the network.
    end kB, kD

    ### Defining New Functions ###
    my_hill_repression(x, v, k, n) = v*k^n/(k^n+x^n)

    # may be necessary to
    # @register my_hill_repression(x, v, k, n)
    # see https://mtk.sciml.ai/stable/tutorials/symbolic_functions/#Registering-Functions-1

    r = @reaction_network MyReactionType begin
        my_hill_repression(x, v_x, k_x, n_x), 0 --> x
    end v_x k_x n_x

    ### Simulating Reaction Networks ###
    probODE = ODEProblem(rn, args...; kwargs...)        # Using multiple dispatch the reaction network can be used as input to create ODE, SDE and Jump problems.
    probSDE = SDEProblem(rn, args...; kwargs...)
    probJump = JumpProblem(prob,aggregator::Direct,rn)
"""

# Declare various arrow types symbols used for the empty set (also 0).
empty_set = Set{Symbol}([:∅])
fwd_arrows = Set{Symbol}([:>, :→, :↣, :↦, :⇾, :⟶, :⟼, :⥟, :⥟, :⇀, :⇁, :⇒, :⟾])
bwd_arrows = Set{Symbol}([:<, :←, :↢, :↤, :⇽, :⟵, :⟻, :⥚, :⥞, :↼, :↽, :⇐, :⟽])
double_arrows = Set{Symbol}([:↔, :⟷, :⇄, :⇆, :⇌, :⇋, :⇔, :⟺])
pure_rate_arrows = Set{Symbol}([:⇐, :⟽, :⇒, :⟾, :⇔, :⟺])

# unfortunately the following doesn't seem to work on 1.5, so supporting these
# operators seems to require us to only support 1.6 and up...
@static if VERSION >= v"1.6.0"
    push!(bwd_arrows, Symbol("<--"))
    push!(double_arrows, Symbol("<-->"))
end

# Declares symbols which may neither be used as parameters not varriables.
forbidden_symbols = [:t, :π, :pi, :ℯ, :im, :nothing, :∅]

### The main macro, takes reaction network notation and returns a ReactionSystem. ###
"""
    @reaction_network

Generates a [`ReactionSystem`](@ref) that encodes a chemical reaction
network.

See the [The Reaction DSL](@ref) documentation for details on
parameters to the macro.

Examples:
```julia
# a basic SIR model, with name SIR
sir_model = @reaction_network SIR begin
    c1, s + i --> 2i
    c2, i --> r
end c1 c2

# a basic SIR model, with random generated name
sir_model = @reaction_network begin
    c1, s + i --> 2i
    c2, i --> r
end c1 c2

# an empty network with name empty
emptyrn = @reaction_network empty

# an empty network with random generated name
emptyrn = @reaction_network
```
"""
macro reaction_network(name::Symbol, ex::Expr, parameters...)
    make_reaction_system(MacroTools.striplines(ex), parameters; name=name)
end

macro reaction_network(ex::Expr, parameters...)
    make_reaction_system(MacroTools.striplines(ex), parameters)
end

#Returns a empty network (with, or without, a declared name)
macro reaction_network(name::Symbol=gensym(:ReactionSystem))
    return Expr(:block,:(@parameters t),
                :(ReactionSystem(Reaction[],
                                 t,
                                 [],
                                 [];
                                 name=$(QuoteNode(name)))))
end

### Macros used for manipulating, and successively builing up, reaction systems. ###

"""
    @add_reactions

Adds the reactions declared to a preexisting [`ReactionSystem`](@ref). All
parameters used in the added reactions need to be declared after the
reactions.

See the [Catalyst.jl for Reaction Models](@ref) documentation for details on
parameters to the macro.
"""
macro add_reactions(rn::Symbol, ex::Expr, parameters...)
    :(merge!($(esc(rn)),$(make_reaction_system(MacroTools.striplines(ex), parameters))))
end

### Sturctures containing information about a reactant, and a reaction, respectively.

#Structure containing information about one reactant in one reaction.
struct ReactantStruct
    reactant::Symbol
    stoichiometry::Number
end
#Structure containing information about one Reaction. Contain all its substrates and products as well as its rate. Contains a specialized constructor.
struct ReactionStruct
    substrates::Vector{ReactantStruct}
    products::Vector{ReactantStruct}
    rate::ExprValues
    only_use_rate::Bool

    function ReactionStruct(sub_line::ExprValues, prod_line::ExprValues, rate::ExprValues, only_use_rate::Bool)
        sub = recursive_find_reactants!(sub_line,1,Vector{ReactantStruct}(undef,0))
        prod = recursive_find_reactants!(prod_line,1,Vector{ReactantStruct}(undef,0))
        new(sub, prod, rate, only_use_rate)
    end
end

### Functions that process the input and rephrase it as a reaction system ###

# Takes the reactions, and rephrases it as a "ReactionSystem" call, as designated by the ModelingToolkit IR.
function make_reaction_system(ex::Expr, parameters; name=gensym(:ReactionSystem))
    reactions = get_reactions(ex)
    reactants = get_reactants(reactions)
    !isempty(intersect(forbidden_symbols,union(reactants,parameters))) && error("The following symbol(s) are used as reactants or parameters: "*((map(s -> "'"*string(s)*"', ",intersect(forbidden_symbols,union(reactants,parameters)))...))*"this is not permited.")
    network_code = Expr(:block,:(@parameters t),:(@variables), :(ReactionSystem([],t,[],[]; name=$(QuoteNode(name)))))
    foreach(parameter-> push!(network_code.args[1].args, parameter), parameters)
    foreach(reactant -> push!(network_code.args[2].args, Expr(:call,reactant,:t)), reactants)
    foreach(parameter-> push!(network_code.args[3].args[6].args, parameter), parameters)
    foreach(reactant -> push!(network_code.args[3].args[5].args, reactant), reactants)
    for reaction in reactions
        subs_init = isempty(reaction.substrates) ? nothing : :([]); subs_stoich_init = deepcopy(subs_init)
        prod_init = isempty(reaction.products) ? nothing : :([]); prod_stoich_init = deepcopy(prod_init)
        reaction_func = :(Reaction($(recursive_expand_functions!(reaction.rate)), $subs_init, $prod_init, $subs_stoich_init, $prod_stoich_init, only_use_rate=$(reaction.only_use_rate)))
        for sub in reaction.substrates
            push!(reaction_func.args[3].args, sub.reactant)
            push!(reaction_func.args[5].args, sub.stoichiometry)
        end
        for prod in reaction.products
            push!(reaction_func.args[4].args, prod.reactant)
            push!(reaction_func.args[6].args, prod.stoichiometry)
        end
        push!(network_code.args[3].args[3].args,reaction_func)
    end
    return network_code
end

#Generates a vector containing a number of reaction structures, each containing the information about one reaction.
function get_reactions(ex::Expr, reactions = Vector{ReactionStruct}(undef,0))
    for line in ex.args
        (line.head != :tuple) && (continue)
        (rate,r_line) = line.args
        (r_line.head  == :-->) && (r_line = Expr(:call,:→,r_line.args[1],r_line.args[2]))

        arrow = r_line.args[1]
        if in(arrow,double_arrows)
            (typeof(rate) == Expr && rate.head == :tuple) || error("Error: Must provide a tuple of reaction rates when declaring a bi-directional reaction.")
            push_reactions!(reactions, r_line.args[2], r_line.args[3], rate.args[1], in(arrow,pure_rate_arrows))
            push_reactions!(reactions, r_line.args[3], r_line.args[2], rate.args[2], in(arrow,pure_rate_arrows))
        elseif in(arrow,fwd_arrows)
            push_reactions!(reactions, r_line.args[2], r_line.args[3], rate, in(arrow,pure_rate_arrows))
        elseif in(arrow,bwd_arrows)
            push_reactions!(reactions, r_line.args[3], r_line.args[2], rate, in(arrow,pure_rate_arrows))
        else
            throw("malformed reaction")
        end
    end
    return reactions
end

#Takes a reaction line and creates reactions from it and pushes those to the reaction array. Used to create multiple reactions from, for instance, 1.0, (X,Y) --> 0.
function push_reactions!(reactions::Vector{ReactionStruct}, sub_line::ExprValues, prod_line::ExprValues, rate::ExprValues, only_use_rate::Bool)
    lengs = [tup_leng(sub_line), tup_leng(prod_line), tup_leng(rate)]
    for i = 1:maximum(lengs)
        (count(lengs.==1) + count(lengs.==maximum(lengs)) < 3) && (throw("malformed reaction"))
        push!(reactions, ReactionStruct(get_tup_arg(sub_line,i), get_tup_arg(prod_line,i), get_tup_arg(rate,i), only_use_rate))
    end
end

#Recursive function that loops through the reaction line and finds the reactants and their stoichiometry. Recursion makes it able to handle weird cases like 2(X+Y+3(Z+XY)).
function recursive_find_reactants!(ex::ExprValues, mult::Int, reactants::Vector{ReactantStruct})
    if typeof(ex)!=Expr
        (ex == 0 || in(ex,empty_set)) && (return reactants)
        if in(ex, getfield.(reactants,:reactant))
            idx = findall(x -> x==ex ,getfield.(reactants,:reactant))[1]
            reactants[idx] = ReactantStruct(ex,mult+reactants[idx].stoichiometry)
        else
            push!(reactants, ReactantStruct(ex,mult))
        end
    elseif ex.args[1] == :*
        recursive_find_reactants!(ex.args[3],mult*ex.args[2],reactants)
    elseif ex.args[1] == :+
        for i = 2:length(ex.args)
            recursive_find_reactants!(ex.args[i],mult,reactants)
        end
    else
        throw("malformed reaction")
    end
    return reactants
end

# Extract the reactants from the set of reactions.
function get_reactants(reactions::Vector{ReactionStruct})
    reactants = Vector{Symbol}()
    for reaction in reactions, reactant in union(reaction.substrates,reaction.products)
        !in(reactant.reactant,reactants) && push!(reactants,reactant.reactant)
    end
    return reactants
end

### Functionality for expanding function call to custom and specific functions ###

#Recursively traverses an expression and replaces special function call like "hill(...)" with the actual corresponding expression.
function recursive_expand_functions!(expr::ExprValues)
    (typeof(expr)!=Expr) && (return expr)
    foreach(i -> expr.args[i] = recursive_expand_functions!(expr.args[i]), 1:length(expr.args))
    if expr.head == :call
        !isdefined(Catalyst,expr.args[1]) && (expr.args[1] = esc(expr.args[1]))
    end
    return expr
end


text:
using LinearAlgebra: mul!

export Layer, forward!, backward!

"A fully-connected layer"
mutable struct Layer
    weights::Matrix{Float64}  # W
    biases::Matrix{Float64}  # b
    activation_fn::Function  # σ
    grad_fn::Function  # ∇σ, or dy/dz

    # Varying
    inputs::Matrix{Float64}  # x
    weighted_sums::Matrix{Float64}  # z = W * x
    activations::Matrix{Float64}  # y = σ(z)
    _z::Union{Nothing, Matrix{Float64}}  # pre-allocated z
end

# Weight initialization:
#
# - The mean of the activations should be zero.
# - The variance of the activations should stay the same across every layer.
function Layer(
    sz::Tuple{Int, Int};
    activation_fn::Function = identity,
    grad_fn::Function = (_, _) -> 1.0,
)
    (n_out, n_in) = sz
    weights = randn(n_out, n_in) .* sqrt(2.0 / (n_out + n_in))  # Xavier initialization
    biases = zeros(n_out, 1)

    inputs, weighted_sums, activations = map(1:3) do _
        Matrix{Float64}(undef, 0, 0)
    end
    _z = nothing

    Layer(weights, biases, activation_fn, grad_fn,
        inputs, weighted_sums, activations, _z)
end

"""
Does the forward pass.

x: n_in × batch_size
Return value: n_out × batch_size
"""
function forward!(l::Layer, x::AbstractMatrix{Float64}; training::Bool=true)::Matrix{Float64}
    (d, n) = size(x)

    W, b, σ = l.weights, l.biases, l.activation_fn

    # Writes to pre-allocated matrix if possible.
    # (This is actually slightly slower on the current toy datasets.)

    # The last minibatch may not have the same size with the previous ones.
    z = if isnothing(l._z) || size(l._z) != (size(W, 1), size(x, 2))
        l._z = W * x .+ b
    else
        l._z .= b
        mul!(l._z, W, x, 1.0, 1.0)
    end
    y = σ.(z)

    if training
        l.inputs = x
        l.weighted_sums = z
        l.activations = y
    end

    y
end

# Learning
# --------
#
# Computes gradients by back-propagation, optimizes with (stochastic) gradient descent.
#
# - Connection: Node (i) -> Node (j)
# - e: error, computed by loss/cost function
# - z: weighted sum, y: activation. y = σ(z). z is logit if σ is sigmoid
# - ∇σ: a function to compute dy/dz that takes arguments of both z and y
#
# Computes `∂e/∂y_i`:
#
#     ∂e/∂y_i = Σ_j w_{ji} ∂e/∂z_j
#     ∂e/∂z_j = ∂e/∂y_j ⋅ ∂y_j/∂z_j
#
#     ⟹
#
#     ∂e/∂y_i = Σ_j w_{ji} ∂e/∂y_j ⋅ ∇σ(z_j, y_j)
#
# Computes `Δw_{ji}` after obtaining `∂e/∂y_i`:
#
#     Δw_{ji} = -η (∂e/∂w_{ji})
#             = -η (∂e/∂z_j) y_i
#     ∂e/∂z_j = ∂e/∂y_j ⋅ ∂y_j/∂z_j
#
#     ⟹
#
#     Δw_{ji} = -η (∂e/∂z_j) y_i
#             = -η (∂e/∂y_j ⋅ ∇σ(z_j, y_j)) y_i
#
# Here the weights are updated layer by layer. In Flux, they are handled in a holistic manner so the
# code is much more concise there.
#
# - - - - - - - -
#
# Within each minibatch, W stays the same across all the inputs. The computational graph is:
#
#         +--(*x1)--> z1 --+
#         |                |
#      W -+--(*x2)--> z2 --+--(+)--> Σe --(/n)--> e
#         |                |
#         +--(*x3)--> z3 --+
#
# The weights (W & b) are computed as the mean rather than sum of all data so it is consistent with
# how the error is computed, and both the error and the magnitude of the gradient are independent of
# batch size.

"""
Given ∂e/∂y of the next layer, does back-propagation, and returns ∂e/∂y of this layer."

∂e_over_∂y: n_out × batch_size
Return value: n_in × batch_size
"""
function backward!(l::Layer, ∂e_over_∂y::AbstractMatrix{Float64}, η::Float64)::Matrix{Float64}
    # Architecture: Layer i -> Layer j. BP: Layer j ⋯> Layer i
    W, b = l.weights, l.biases
    y_i, z_j, y_j = l.inputs, l.weighted_sums, l.activations
    ∇σ, ∇y_j = l.grad_fn, ∂e_over_∂y
    n = size(y_i, 2)

    ∇σ_j = ∇σ.(z_j, y_j)
    ∇z_j = ∇y_j .* ∇σ_j
    ∇y_i = W' * ∇z_j
    ∇W = (∇z_j * y_i') / n  # sum of all data in the batch
    ∇b = sum(∇z_j, dims=2) / n
    W .-= η * ∇W
    b .-= η * ∇b
    ∇y_i
end

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
using Distributed
using CUDAnative

const ndevices = length(CUDAnative.devices())

# add one worker per device
addprocs(ndevices)

using CuArrays
using DistributedArrays

# Set device to be used on each worker
asyncmap(collect(zip(workers(), CUDAnative.devices()))) do (p, d)
    remotecall_wait(() -> CUDAnative.device!(d), p)
    nothing
end

A = DArray((400, 400)) do I
    dims = map(length, I)
    reshape(CuArrays.CURAND.curand(Float32, prod(dims)), dims)
end;

sum(A)

# simple broadcast
B = A .^2

C = A .+ A .* sin.(B)


text:
##### Beginning of file

import Test # stdlib

Test.@test( isdir(PredictMDExtra.package_directory()) )

Test.@test( isdir(PredictMDExtra.package_directory("ci")) )

Test.@test( isdir(PredictMDExtra.package_directory("ci", "travis")) )

Test.@test( isdir(PredictMDExtra.package_directory(PredictMDExtraTestModuleA)) )

Test.@test( isdir(PredictMDExtra.package_directory(PredictMDExtraTestModuleB)) )

Test.@test( isdir( PredictMDExtra.package_directory(
            PredictMDExtraTestModuleB, "directory2",
            ) ) )

Test.@test( isdir( PredictMDExtra.package_directory(
            PredictMDExtraTestModuleB, "directory2", "directory3",
            ) ) )

Test.@test_throws(
    ErrorException,
    PredictMDExtra.package_directory(PredictMDExtraTestModuleC),
    )

##### End of file


text:
const has_bokeh_png_deps = false

function initialize_bokeh()
    ispynull(bokeh) || return nothing
    try
        copy!(bokeh, _import_dependency("bokeh", "bokeh"; channel="conda-forge"))
    catch err
        copy!(bokeh, PyNULL())
        throw(err)
    end
    return nothing
end

# install dependencies for saving PNGs if using conda
function initialize_bokeh_png_deps()
    has_bokeh_png_deps && return nothing
    try
        _import_dependency("selenium", "selenium"; channel="conda-forge")
        has_bokeh_png_deps = true
    catch err
        has_bokeh_png_deps = false
        throw(err)
    end
    return nothing
end

load_backend(::Val{:bokeh}) = initialize_bokeh()

convert_result(f, axis, ::Val{:bokeh}) = BokehPlot(axis)
function convert_result(f, axes::AbstractArray, ::Val{:bokeh})
    return BokehPlot(arviz.plots.backends.create_layout(axes))
end

"""
    BokehPlot(::PyObject)

Loose wrapper around a Bokeh figure, mostly used for dispatch.

In most cases, use one of the plotting functions with `backend=:bokeh` to create a
`BokehPlot` instead of using a constructor.
"""
struct BokehPlot
    o::PyObject
end

BokehPlot(plot::BokehPlot) = plot

@inline PyObject(plot::BokehPlot) = getfield(plot, :o)

Base.convert(::Type{BokehPlot}, o::PyObject) = BokehPlot(o)

Base.hash(plot::BokehPlot) = hash(PyObject(plot))

Base.propertynames(plot::BokehPlot) = propertynames(PyObject(plot))

function Base.getproperty(plot::BokehPlot, name::Symbol)
    o = PyObject(plot)
    name === :o && return o
    return getproperty(o, name)
end

function render_html(plot::BokehPlot, name=nothing)
    obj = PyObject(plot)
    return bokeh.embed.file_html(obj, bokeh.resources.CDN, name)
end

Base.display(::REPL.REPLDisplay, plot::BokehPlot) = bokeh.plotting.show(plot)

Base.show(io::IO, ::MIME"text/html", plot::BokehPlot) = print(io, render_html(plot))
function Base.show(
    io::IO,
    ::Union{MIME"application/juno+plotpane",MIME"application/prs.juno.plotpane+html"},
    plot::BokehPlot,
)
    return print(io, render_html(plot))
end
function Base.show(io::IO, ::MIME"juliavscode/html", plot::BokehPlot)
    return print(io, render_html(plot))
end
function Base.show(io::IO, ::MIME"image/png", plot::BokehPlot)
    initialize_bokeh_png_deps()
    image = bokeh.io.export.get_screenshot_as_png(plot)
    print(io, image._repr_png_())
    return nothing
end

"""
    write(io::IO, plot::BokehPlot)
    write(filename::AbstractString, plot::BokehPlot)

Write the HTML representation of the Bokeh plot to the I/O stream or file.
"""
Base.write(io::IO, plot::BokehPlot) = print(io, render_html(plot))


text:
##########  Configuration   #########

have_imagemagick = false
@unix_only which_cmd = "which"
@windows_only which_cmd = "where"
_, p = readsfrom(`$which_cmd convert`)
wait(p)
if !success(p)
    warn("ImageMagick utilities not found. Install for more file format support.")
else
    have_imagemagick = true
end


text:
"""
Automatic Differentiation
=========================

## References

- https://en.wikipedia.org/wiki/Automatic_differentiation
- https://github.com/MikeInnes/diff-zoo
- https://rufflewind.com/2016-12-30/reverse-mode-automatic-differentiation
"""

module AutomaticDifferentiation

include("forward.jl")
include("backward.jl")

#=

∂y/∂x can be expressed as (supposing the computation forms a polytree):

$$
{\partial{y} \over \partial{x}} =
  \prod_{i=1}^{n} {\partial{w_i} \over \partial{w_{i-1}}} =
  {\partial{y} \over \partial{w_{n-1}}} \cdot
  {\partial{w_{n-1}} \over \partial{w_{n-2}}} \cdot \ldots \cdot
  {\partial{w_1} \over \partial{x}}
$$

- Forward: Evaluates ∂y/∂x with right association
- Backward: Evaluates ∂y/∂x with left association

=#

end # module


text:
"""
returns sigma for  j
"""
function sigma(j,as,bs,jmax)
    rho=(j-as)/bs
1/(1+exp(-rho))
end

"""
returns u-function for j
"""
function ushape(j,a1m,a2m,b1m,b2m,cm,d1m,d2m,jmax)
        rho1=(a1m-j)/b1m
        rho2=(j-a2m)/b2m
        cm+
    if j<a1m rho1^d1m
            elseif a1m<=j<=a2m 0
        elseif a2m<=j<=jmax rho2^d2m
    end
end

"""
returns wieght for t,j (according to the ovsyannikov data)
"""
function w(t,j,jmax)
   wv=[0.012 0.052 0.11 0.179 0.253 0.323 0.483 0.677 0.908 1.116 1.297 1.541 1.758 1.971 1.997 2.358 2.475];
    return wv[j]
end

"""
wtj
"""
function gwtj(w,tmax,jmax)
  wtj= fill(1.0,tmax,jmax);
    for t in [1:tmax;]
        for j in [1:jmax;] wtj[t,j]=w(t,j,jmax)
            end
    end
    return wtj
end

"""
fills H with 1's
"""
function h(t,j,jmax)
    hv=fill(1.0,jmax)';
    return hv[j]
end

"""
mtj
"""
function ghtj(h,tmax,jmax)
   htj= fill(0.0,tmax,jmax);
    for t in [1:tmax;]
        for j in [1:jmax;] htj[t,j]=h(t,j,jmax)
            end
    end
    return htj
end

"""
mtj
"""

function gmtj(a1m,a2m,b1m,b2m,cm,d1m,d2m,htj,tmax,jmax)
   mtj= fill(1.0,tmax,jmax);
    for t in [1:tmax;]
        for j in [1:jmax;] mtj[t,j]= htj[t,j]*ushape(j,a1m,a2m,b1m,b2m,cm,d1m,d2m,jmax)
            end
    end
    return mtj
end

"""
ft=0.2
"""
function gft(tmax,jmax)
    ft= fill(0.2,tmax);
    return ft
end

"""
ftj
"""
function gftj(as,bs,ft,tmax,jmax)
    ftj= fill(1.0,tmax,jmax);
    for t in [1:tmax;]
        for j in [1:jmax;] ftj[t,j]=ft[t]*sigma(j,as,bs,jmax)
            end
    end
    return ftj
end

"""
x
"""

function gx(x,t0,t1,tmax,j0,j1,jmax,M,F)
    for t in [t0:t1;]

            for i in [0:j1-j0-1;]
            #println(t+i+1," ",j0+i+1," ",x[t+i,j0+i])
                x[t+i+1,j0+i+1]=x[t+i,j0+i]*exp(-M[t+i,j0+i]-F[t+i,j0+i])

            end

    end
    return x
end

"""
y
"""
function gy(y,x,t0,t1,tmax,j0,j1,jmax,M,F,wtj,kx,kw)
    for t in [t0:t1;]
        for i in [0:j1-j0;]
      y[t+i,j0+i]=(kw*wtj[t+i,j0+i])*(kx*x[t+i,j0+i])*F[t+i,j0+i]/(M[t+i,j0+i]+F[t+i,j0+i])*(1- exp(-M[t+i,j0+i]-F[t+i,j0+i]))

        end
    end
    return y
end

"""
x,y
"""
function gm_burago(tmax,jmax,M,F,W,kx,kw,x1j)
#0
    x=fill(0.0,tmax,jmax);

    for j in [1:jmax;] x[1,j]=x1j[j] end

    for t in 1:tmax x[t,1]=x1j[1] end

 #1
        for t in [1:tmax-jmax+1;]
           for i in [0:jmax-2;]
            #println(t+i+1," ",i+2,"<-", t+i," ",i+1)
            x[t+i+1,i+1+1]=x[t+i,i+1]*exp(-M[t+i,i+1]-F[t+i,i+1])
           end
        end

 #2
    for t in [tmax-jmax+2:tmax-1;]
           for i in [0:tmax-t-1;]
            #println(t+i+1," ",i+2,"<-", t+i," ",i+1)
            x[t+i+1,i+1+1]=x[t+i,i+1]*exp(-M[t+i,i+1]-F[t+i,i+1])
           end
        end
   #3

    for j in [2:jmax;]
           for i in [0:jmax-j-1;]
           #println(i+1+1,"! ",j+i+1,"<-", i+1," ",j+i)
            x[i+1+1,j+i+1]=x[i+1,j+i]*exp(-M[i+1,j+i]-F[i+1,j+i])
           end
        end

 #0
     y=fill(0.0,tmax,jmax);

#1
   for t in [1:tmax-jmax+1;]
        for i in [0:jmax-1;]
            #println(t+i," ",i+1)
            Z=M[t+i,i+1]+F[t+i,i+1]
            y[t+i,i+1]=
            (kw*W[t+i,i+1])*(kx*x[t+i,i+1])*F[t+i,i+1]/Z*(1- exp(-Z))

        end
    end

   #2
     for t in [tmax-jmax+2:tmax;]
        for i in [0:tmax-t;]
            #println(t+i," ",i+1)
            Z=M[t+i,i+1]+F[t+i,i+1]
            y[t+i,i+1]=
            (kw*W[t+i,i+1])*(kx*x[t+i,i+1])*F[t+i,i+1]/Z*(1- exp(-Z))

        end
    end

    #3

    for j in [2:jmax;]
           for i in [0:jmax-j;]

     Z=M[i+1,j+i]+F[i+1,j+i]
            y[i+1,j+i]=
            (kw*W[i+1,j+i])*(kx*x[i+1,j+i])*F[i+1,j+i]/Z*(1- exp(-Z))

        end
    end

    #gy(y,x,1,tmax-jmax+1,tmax,1,jmax,jmax,M,F,W,kx,kw)

    return x,y
end

"""
x,y
"""

function gm(tmax,jmax,M,F,W,kx,kw,x1j,xt1)
#0
    x=fill(0.0,tmax,jmax);

    for j in [1:jmax;] x[1,j]=x1j[j] end

    for t in 1:tmax x[t,1]=xt1[t] end

 #1
        for t in [1:tmax-jmax+1;]
           for i in [0:jmax-2;]
            #println(t+i+1," ",i+2,"<-", t+i," ",i+1)
            x[t+i+1,i+1+1]=x[t+i,i+1]*exp(-M[t+i,i+1]-F[t+i,i+1])
           end
        end

 #2
    for t in [tmax-jmax+2:tmax-1;]
           for i in [0:tmax-t-1;]
            #println(t+i+1," ",i+2,"<-", t+i," ",i+1)
            x[t+i+1,i+1+1]=x[t+i,i+1]*exp(-M[t+i,i+1]-F[t+i,i+1])
           end
        end
   #3

    for j in [2:jmax;]
           for i in [0:jmax-j-1;]
           #println(i+1+1,"! ",j+i+1,"<-", i+1," ",j+i)
            x[i+1+1,j+i+1]=x[i+1,j+i]*exp(-M[i+1,j+i]-F[i+1,j+i])
           end
        end

 #0
     y=fill(0.0,tmax,jmax);

#1
   for t in [1:tmax-jmax+1;]
        for i in [0:jmax-1;]
            #println(t+i," ",i+1)
            Z=M[t+i,i+1]+F[t+i,i+1]
            y[t+i,i+1]=
            (kw*W[t+i,i+1])*(kx*x[t+i,i+1])*F[t+i,i+1]/Z*(1- exp(-Z))
        end
    end

   #2
     for t in [tmax-jmax+2:tmax;]
        for i in [0:tmax-t;]
            #println(t+i," ",i+1)
            Z=M[t+i,i+1]+F[t+i,i+1]
            y[t+i,i+1]=
            (kw*W[t+i,i+1])*(kx*x[t+i,i+1])*F[t+i,i+1]/Z*(1- exp(-Z))
        end
    end

    #3

    for j in [2:jmax;]
           for i in [0:jmax-j;]

     Z=M[i+1,j+i]+F[i+1,j+i]
            y[i+1,j+i]=
            (kw*W[i+1,j+i])*(kx*x[i+1,j+i])*F[i+1,j+i]/Z*(1- exp(-Z))
        end
    end

    #gy(y,x,1,tmax-jmax+1,tmax,1,jmax,jmax,M,F,W,kx,kw)

    return x,y
end

"""
adds some noise with d(m,s) (laplace, normal) to a matrix
"""
function noise2ts_burago(ts,dn,dc)
    (tmax,jmax)=size(ts.x)
      x1=fill(0.0,tmax,jmax)
    y1=fill(0.0,tmax,jmax)
    for t in 1:tmax
        for j in 1:jmax
            x1[t,j]=ts.x[t,j]*exp(rand(dn))
            y1[t,j]=ts.y[t,j]*exp(rand(dc))
        end
    end
    return cohorts_t(ts.year_b,ts.year_e,x1,y1,ts.w)

end

"""
constraint
"""
function buragoconstraint1(x::Vector, jmax)
    tmax=div((length(x)-9),4)
        x[4*tmax+1]-x[4*tmax+2]
end
"""
constraint
"""
function buragoconstraint2(x::Vector, jmax)
     tmax=div((length(x)-9),4)
    x[4*tmax+2]-jmax
end
"""
constraint
"""
 function buragoconstraint3(x::Vector, jmax)
     tmax=div((length(x)-9),4)
    x[4*tmax+8]-jmax
end
"""
constraint
"""
function buragoconstraint4(x::Vector, kappa)
     tmax=div((length(x)-9),4)
    p=true
    for i in 1:3*tmax
        p=p&(1/kappa<=x[i]<=kappa)
    end
    if p
        -1.
    else 1.
    end
end
"""
constraint
"""
function buragoconstraint5(x::Vector,i, kappa1)
    kappa1-x[i]
    end
"""
constraint
"""
function buragoconstraint6(x::Vector,i, kappa2)
    x[i]-kappa2
    end
"""
criterion (Laplace) evaluation for cohorts of cst type/Lp, p=1
"""
function buragofunc(x::Vector, grad::Vector)

    global cst
    (tmax,jmax)=size(cst.x)

    local Hl=fill(1.0,tmax,jmax)
    local ftl=fill(1,tmax)

    global kx,kw
    global count

    if length(grad) > 0
        grad[1] = 0
        grad[2] = 0.5/sqrt(x[2])
    end

    count::Int += 1

for t in 1:tmax
        Hl[t,1:3]=x[3*(t-1)+1:3*t]
    end

    local a1ml,a2ml,b1ml,b2ml,cml,d1ml,d2ml,asl,bsl
    ftl=x[3*tmax+1:4*tmax]
    a1ml,a2ml,b1ml,b2ml,cml,d1ml,d2ml=x[4*tmax+1:4*tmax+7]
    asl,bsl=x[4*tmax+8:4*tmax+9]

    #println(cst,psl,Hl,ftl)
    return Lp(cst,a1ml,a2ml,b1ml,b2ml,cml,d1ml,d2ml,Hl,asl,bsl,ftl,kx,kw,1)

    #println("f_$count($x)=$y")

end

"""
criterion (Gauss) evaluation for cohorts of cst type/Lp, p=2
"""
#p=2
function buragofuncp2(x::Vector, grad::Vector)

    global cst
    (tmax,jmax)=size(cst.x)

    local Hl=fill(1.0,tmax,jmax)
    local ftl=fill(1,tmax)

    global kx,kw
    global count

    if length(grad) > 0
        grad[1] = 0
        grad[2] = 0.5/sqrt(x[2])
    end

    count::Int += 1

for t in 1:tmax
        Hl[t,1:3]=x[3*(t-1)+1:3*t]
    end

    local a1ml,a2ml,b1ml,b2ml,cml,d1ml,d2ml,asl,bsl
    ftl=x[3*tmax+1:4*tmax]
    a1ml,a2ml,b1ml,b2ml,cml,d1ml,d2ml=x[4*tmax+1:4*tmax+7]
    asl,bsl=x[4*tmax+8:4*tmax+9]

    #println(cst,psl,Hl,ftl)
    return Lp(cst,a1ml,a2ml,b1ml,b2ml,cml,d1ml,d2ml,Hl,asl,bsl,ftl,kx,kw,2)

    #println("f_$count($x)=$y")

end

"""
criterion (normal, deviation) evaluation for cohorts of cst type/Lmine2, p=2
"""
function minefuncp2(x::Vector, grad::Vector)

    global cst
    (tmax,jmax)=size(cst.x)

    local Hl=fill(1.0,tmax,jmax)
    local ftl=fill(1,tmax)

    global kx,kw
    global count

    if length(grad) > 0
        grad[1] = 0
        grad[2] = 0.5/sqrt(x[2])
    end

    count::Int += 1

for t in 1:tmax
        Hl[t,1:3]=x[3*(t-1)+1:3*t]
    end

    local a1ml,a2ml,b1ml,b2ml,cml,d1ml,d2ml,asl,bsl
    ftl=x[3*tmax+1:4*tmax]
    a1ml,a2ml,b1ml,b2ml,cml,d1ml,d2ml=x[4*tmax+1:4*tmax+7]
    asl,bsl=x[4*tmax+8:4*tmax+9]

    #println(cst,psl,Hl,ftl)
    return Lmine2(cst,a1ml,a2ml,b1ml,b2ml,cml,d1ml,d2ml,Hl,asl,bsl,ftl,kx,kw,2)

    #println("f_$count($x)=$y")

end

"""
criterion (normal, deviation) evaluation for cohorts of cst type/Lmine2max, p=2
"""
function minefuncp2max(x::Vector, grad::Vector)

    global cst
    (tmax,jmax)=size(cst.x)

    local Hl=fill(1.0,tmax,jmax)
    local ftl=fill(1,tmax)

    global kx,kw
    global count

    if length(grad) > 0
        grad[1] = 0
        grad[2] = 0.5/sqrt(x[2])
    end

    count::Int += 1

for t in 1:tmax
        Hl[t,1:3]=x[3*(t-1)+1:3*t]
    end

    local a1ml,a2ml,b1ml,b2ml,cml,d1ml,d2ml,asl,bsl
    ftl=x[3*tmax+1:4*tmax]
    a1ml,a2ml,b1ml,b2ml,cml,d1ml,d2ml=x[4*tmax+1:4*tmax+7]
    asl,bsl=x[4*tmax+8:4*tmax+9]

    #println(cst,psl,Hl,ftl)
    return Lmine2max(cst,a1ml,a2ml,b1ml,b2ml,cml,d1ml,d2ml,Hl,asl,bsl,ftl,kx,kw,2)

    #println("f_$count($x)=$y")

end

"""
criterion (normal, deviation) evaluation for cohorts of cst type/Lp12, p=2
"""
function buragofunc12(x::Vector, grad::Vector)

    global cst
    (tmax,jmax)=size(cst.x)

    local Hl=fill(1.0,tmax,jmax)
    local ftl=fill(1,tmax)

    global kx,kw
    global count

    if length(grad) > 0
        grad[1] = 0
        grad[2] = 0.5/sqrt(x[2])
    end

    count::Int += 1

for t in 1:tmax
        Hl[t,1:3]=x[3*(t-1)+1:3*t]
    end

    local a1ml,a2ml,b1ml,b2ml,cml,d1ml,d2ml,asl,bsl
    ftl=x[3*tmax+1:4*tmax]
    a1ml,a2ml,b1ml,b2ml,cml,d1ml,d2ml=x[4*tmax+1:4*tmax+7]
    asl,bsl=x[4*tmax+8:4*tmax+9]

    #println(cst,psl,Hl,ftl)
    return Lp12(cst,a1ml,a2ml,b1ml,b2ml,cml,d1ml,d2ml,Hl,asl,bsl,ftl,kx,kw,1)

    #println("f_$count($x)=$y")

end

"""
list of optimization algorithms
"""
ALGS=(:LN_NELDERMEAD,:LN_SBPLX,:LN_COBYLA,:LN_NEWUOA,:LN_PRAXIS,:LN_BOBYQA,:LN_AUGLAG,:LN_AUGLAG_EQ)

#

"""
optimizes func with 2 algs (a -main and a1 - local), (burago)func (Lp, Lp12, etc)  and 4 buragoconstraints
"""
function optimizeLf(x,func,tmax,jmax,a,a1,lb,lb1,xtol,xtol1,bc4,ctol,is1,is2,is3,pop,pop1,maxe)

global count = 0# keeps track of # of function evaluations
opt = Opt(ALGS[a], 4*tmax+9)
l_opt= Opt(ALGS[a1], 4*tmax+9)

    println(ALGS[a])
    println(ALGS[a1])

    lower_bounds!(opt, fill(lb,4*tmax+9))
    lower_bounds!(l_opt, fill(lb1,4*tmax+9))
    xtol_rel!(opt,xtol)
    xtol_rel!(l_opt,xtol1)

min_objective!(opt, func)
    inequality_constraint!(opt, (x,g) -> buragoconstraint1(x,jmax), ctol)
    inequality_constraint!(opt, (x,g) -> buragoconstraint2(x,jmax), ctol)
    inequality_constraint!(opt, (x,g) -> buragoconstraint3(x,jmax),ctol )
    inequality_constraint!(opt, (x,g) -> buragoconstraint4(x,bc4), ctol)
#for i in 1:3*tmax
#    inequality_constraint!(opt, (x,g) -> buragoconstraint5(x,i,0.1), 1e-8)
#     inequality_constraint!(opt, (x,g) -> buragoconstraint6(x,i,5.), 1e-8)
#end

local_optimizer!(opt,l_opt)

is=fill(0.4,4*tmax+9)
    is[1:3*tmax]=fill(is1,3*tmax)#h
    is[3*tmax+1:4*tmax]=fill(is2,tmax)#f
    is[4*tmax+1:4*tmax+9]=fill(is3,9)#pm&pf

initial_step!(opt,is)
population!(opt,pop)
population!(l_opt,pop1)
NLopt.srand_time()
    maxeval!(opt, maxe)

(minf,minx,ret) = NLopt.optimize(opt,x)

end

"""
optimizes Lp with 2 algs (a -main and a1 - local), buragofunc   and 4 buragoconstraints
"""
#
function optimizeL(x,tmax,jmax,a,a1,lb,lb1,xtol,xtol1,bc4,ctol,is1,is2,is3,pop,pop1,maxe)

global count = 0# keep track of # of function evaluations
opt = Opt(ALGS[a], 4*tmax+9)
l_opt= Opt(ALGS[a1], 4*tmax+9)

   # println(ALGS[a])
    # println(ALGS[a1])

    lower_bounds!(opt, fill(lb,4*tmax+9))
    lower_bounds!(l_opt, fill(lb1,4*tmax+9))
    xtol_rel!(opt,xtol)
    xtol_rel!(l_opt,xtol1)

min_objective!(opt, buragofunc)
    inequality_constraint!(opt, (x,g) -> buragoconstraint1(x,jmax), ctol)
    inequality_constraint!(opt, (x,g) -> buragoconstraint2(x,jmax), ctol)
    inequality_constraint!(opt, (x,g) -> buragoconstraint3(x,jmax),ctol )
    inequality_constraint!(opt, (x,g) -> buragoconstraint4(x,bc4), ctol)
#for i in 1:3*tmax
#    inequality_constraint!(opt, (x,g) -> buragoconstraint5(x,i,0.1), 1e-8)
#     inequality_constraint!(opt, (x,g) -> buragoconstraint6(x,i,5.), 1e-8)
#end

local_optimizer!(opt,l_opt)

is=fill(0.4,4*tmax+9)
    is[1:3*tmax]=fill(is1,3*tmax)#h
    is[3*tmax+1:4*tmax]=fill(is2,tmax)#f
    is[4*tmax+1:4*tmax+9]=fill(is3,9)#pm&pf

initial_step!(opt,is)
population!(opt,pop)
population!(l_opt,pop1)
NLopt.srand_time()
    maxeval!(opt, maxe)

(minf,minx,ret) = NLopt.optimize(opt,x)

end

"""
optimizes Lp with 1 algs (a -main), buragofunc   and 4 buragoconstraints, no lower_bounds, no x_tol, no local_optimizer?
no population
"""
function optimizeLs(x,tmax,jmax,a,lb,lb1,xtol,xtol1,bc4,ctol,is1,is2,is3,pop,pop1,maxe)

   global count = 0
opt = Opt(ALGS[a], 4*tmax+9)
#l_opt= Opt(ALGS[a1], 4*tmax+9)
    lower_bounds!(opt, fill(lb,4*tmax+9))
    #lower_bounds!(l_opt, fill(lb1,4*tmax+9))
    xtol_rel!(opt,xtol)
    #xtol_rel!(l_opt,xtol1)

min_objective!(opt, buragofunc)
    inequality_constraint!(opt, (x,g) -> buragoconstraint1(x,jmax), ctol)
    inequality_constraint!(opt, (x,g) -> buragoconstraint2(x,jmax), ctol)
    inequality_constraint!(opt, (x,g) -> buragoconstraint3(x,jmax),ctol )
    inequality_constraint!(opt, (x,g) -> buragoconstraint4(x,bc4), ctol)
#for i in 1:3*tmax
#    inequality_constraint!(opt, (x,g) -> buragoconstraint5(x,i,0.1), 1e-8)
#     inequality_constraint!(opt, (x,g) -> buragoconstraint6(x,i,5.), 1e-8)
#end

#local_optimizer!(opt,l_opt)

is=fill(0.4,4*tmax+9)
    is[1:3*tmax]=fill(is1,3*tmax)#h
    is[3*tmax+1:4*tmax]=fill(is2,tmax)#f
    is[4*tmax+1:4*tmax+9]=fill(is3,9)#pm&pf

initial_step!(opt,is)
population!(opt,pop)
#population!(l_opt,pop1)
NLopt.srand_time()
    maxeval!(opt, maxe)

(minf,minx,ret) = NLopt.optimize(opt,x)

end

"""
optimizes Lp with 2 algs (a -main and a1 - local), buragofunc_v   and 4 buragoconstraints
"""
function optimizeLp_v2(x,tmax,jmax,a,a1,lb,lb1,xtol,xtol1,bc4,ctol,is1,is2,is3,pop,pop1,maxe)

global count = 0# keep track of # of function evaluations
opt = Opt(ALGS[a], 4*tmax+9)
l_opt= Opt(ALGS[a1], 4*tmax+9)

    println(ALGS[a])
    println(ALGS[a1])

    lower_bounds!(opt, fill(lb,4*tmax+9))
    lower_bounds!(l_opt, fill(lb1,4*tmax+9))
    xtol_rel!(opt,xtol)
    xtol_rel!(l_opt,xtol1)

min_objective!(opt, buragofunc_v)
    inequality_constraint!(opt, (x,g) -> buragoconstraint1(x,jmax), ctol)
    inequality_constraint!(opt, (x,g) -> buragoconstraint2(x,jmax), ctol)
    inequality_constraint!(opt, (x,g) -> buragoconstraint3(x,jmax),ctol )
    inequality_constraint!(opt, (x,g) -> buragoconstraint4(x,bc4), ctol)
#for i in 1:3*tmax
#    inequality_constraint!(opt, (x,g) -> buragoconstraint5(x,i,0.1), 1e-8)
#     inequality_constraint!(opt, (x,g) -> buragoconstraint6(x,i,5.), 1e-8)
#end

local_optimizer!(opt,l_opt)

is=fill(0.4,4*tmax+9)
    is[1:3*tmax]=fill(is1,3*tmax)#h
    is[3*tmax+1:4*tmax]=fill(is2,tmax)#f
    is[4*tmax+1:4*tmax+9]=fill(is3,9)#pm&pf

initial_step!(opt,is)
population!(opt,pop)
population!(l_opt,pop1)
NLopt.srand_time()
    maxeval!(opt, maxe)

(minf,minx,ret) = NLopt.optimize(opt,x)

end

"""
evaluates vector of unknown pars with use of given
 ALGS=(:LN_NELDERMEAD, :LN_SBPLX, :LN_COBYLA, :LN_NEWUOA, :LN_PRAXIS, :LN_BOBYQA, :LN_AUGLAG, :LN_AUGLAG_EQ)
from NLopt and LN_AUGLAG (7),
xx - initial approximation for pars,
k - number of iterations,
algs - list of algs
"""
function opt_v(xx,k,algs)
@eval @everywhere  xxx=$xx
flag=true
ni=0
result=0.0
while flag
    ni+=1
    println(ni)
    pm= pmap(i->optimizeLp_v2(xxx,tmax,jmax,7,i,0.0,0.0,1e-20,1e-20,2.,1e-20,1.,1.,0.1,5,5,50000),algs)

    pm1=filter(x->typeof(x)==Tuple{Float64,Array{Float64,1},Symbol},pm)
   minfs=map(i->pm1[i][1],1:length(pm1))

    println(minfs)
    result=minfs

    if minimum(minfs)==maximum(minfs) || ni>k
         println( minimum(minfs)," ",maximum(minfs));
            flag=false
        else
        println( minimum(minfs)," ",maximum(minfs));
        mm=minimum(minfs);
        flag=true;
            @eval @everywhere  xxx=$pm[findfirst(x->$mm==x,$minfs)][2]
    end
end

    return(xxx)
end

max_stars_repo_name: picestcode/Fishmetica

max_issues_repo_name: picestcode/Fishmetica

max_forks_repo_name: picestcode/Fishmetica


text:
using cart_pole_example
using Test

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


text:
module Entropy
using Docile
@docstrings
using HDF5, JLD
include("types.jl")
include("nsb2.jl")
include("utils.jl")
include("bias.jl")
include("hist.jl")
using .Histograms
include("binless.jl")
include("binned.jl")
include("algorithms.jl")
include("process.jl")
import GUICheck
if GUICheck.hasgui()
	include("plot.jl")
end

export Hist2d

end # module


text:
using Vann
using Base.Test

include("model_tests.jl")
include("calib_tests.jl")

max_stars_repo_name: jmgnve/Vann

max_forks_repo_name: jmgnve/Vann


text:
using AcuteML
@aml mutable struct Body "~"
    h1, "~"
    p::Vector{String}, "~"
end

@aml mutable struct Page doc"html"
    body::Body, "~"
end

b = Body(h1 = "My heading", p = ["Paragraph1", "Paragraph2"])
d = Page(body = b)
pprint(d)


text:
module Augmentor

using MappedArrays
using ImageCore
using ImageTransformations
using ImageFiltering
using CoordinateTransformations
using Rotations
using Interpolations
using StaticArrays
using IdentityRanges
using MLDataPattern
using ComputationalResources
using FileIO
using Base.PermutedDimsArrays: PermutedDimsArray
using LinearAlgebra
using OffsetArrays

# axes(::OffsetArray) changes from Base.Slice to Base.IdentityUnitRange in julia 1.1
# https://github.com/JuliaArrays/OffsetArrays.jl/pull/62
# TODO: switch to Base.IdentityUnitRange when we decide to drop 1.0 compatibility
using OffsetArrays: IdentityUnitRange

using InteractiveUtils: methodswith

export

    CPU1,
    CPUThreads,

    Gray,
    RGB,

    SplitChannels,
    CombineChannels,
    PermuteDims,
    Reshape,

    ConvertEltype,
    MapFun,
    AggregateThenMapFun,

    Rotate90,
    Rotate180,
    Rotate270,
    Rotate,

    ShearX,
    ShearY,

    FlipX,
    FlipY,

    Crop,
    CropNative,
    CropSize,
    CropRatio,
    RCropRatio,

    Resize,

    Scale,
    Zoom,

    CacheImage,
    NoOp,
    Either,

    ElasticDistortion,

    augment,
    augment!,
    augmentbatch!,

    testpattern

include("compat.jl")
include("utils.jl")
include("types.jl")
include("operation.jl")

include("operations/channels.jl")
include("operations/dims.jl")
include("operations/convert.jl")
include("operations/mapfun.jl")

include("operations/noop.jl")
include("operations/cache.jl")
include("operations/rotation.jl")
include("operations/shear.jl")
include("operations/flip.jl")
include("operations/crop.jl")
include("operations/resize.jl")
include("operations/scale.jl")
include("operations/zoom.jl")
include("operations/either.jl")

include("distortionfields.jl")
include("distortedview.jl")
include("operations/distortion.jl")

include("pipeline.jl")
include("codegen.jl")
include("augment.jl")
include("augmentbatch.jl")

function __init__()
    if VERSION < v"1.3"
        # see compat.jl
        rand_mutex[] = Threads.Mutex()
    end
end

end # module


text:
function oper(+)
    @eval import Base.$+
    @eval function $+(r::FinMCDict{num1, num2}, d::FinMCDict{num1, num3}) where {num1 <: Number, num2 <: Number, num3 <: Number}
        T_r = collect(keys_(r))
        T_d = collect(keys_(d))
        if (length(T_d) > length(T_r))
            return $+(d, r)
        end
        d_complete = complete_(T_r, d)
        #key_type=typeof(zero(num1));
        val_type = typeof(zero(num2) + zero(num3))
        out = FinMCDict{num1, val_type}()
        for t in T_r
            insert!(out, t, $+(r[t], d_complete[t]))
        end

        return out
    end
    @eval function $+(r::FinMCDict{num1, num2}, d::num3) where {num1 <: Number, num2 <: Number, num3 <: Number}
        return $+.(r, d)
    end
end
function oper2(+)
    @eval import Base.$+
    @eval function $+(d::num3, r::FinMCDict{num1, num2}) where {num1 <: Number, num2 <: Number, num3 <: Number}
        return $+(r, d)
    end
end

import Base.-
function -(d::num3, r::FinMCDict{num1, num2}) where {num1 <: Number, num2 <: Number, num3 <: Number}
    return +(-1 * r, d)
end

oper(Symbol(+))
oper(Symbol(-))
oper(Symbol(/))
oper(Symbol(*))

oper2(Symbol(*))
oper2(Symbol(+))

#keys__(x::Dictionary)=x.indices
#keys__(x)=keys_(x)
#
#values__(x::Dictionary)=x.values
#values__(x)=values_(x)

function complete_(T::Array{num}, d::FinMCDict{num1, num2}) where {num <: Number, num1 <: Number, num2 <: Number}
    T_d = collect(keys_(d))
    idx_ = sortperm(T_d)
    T_d = T_d[idx_]
    d_val = collect(values_(d))
    d_val = d_val[idx_]
    out = FinMCDict{num1, num2}()
    itp = LinearInterpolation(T_d, d_val)
    for t in T
        insert!(out, t, itp(t))
    end
    return out
end


text:
using Test
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL

@testset "GLVertexArray.jl" begin

   # ffunction GLVertexArray()
   @testset "GLVertexArray" begin
      @test
      @test
      @test
      @test
   end

   # function releaseGpuResources(array::GLVertexArray)
   @testset "releaseGpuResources" begin
      @test
      @test
      @test
      @test
   end

   # function enableVertexArray(array::GLVertexArray)
   @testset "enableVertexArray" begin
      @test
      @test
      @test
      @test
   end

   # function disableVertexArray(array::GLVertexArray)
   @testset "disableVertexArray" begin
      @test
      @test
      @test
      @test
   end

end


text:
@defcomp MarketDamages begin
    region = Index()
    y_year = Parameter(index=[time], unit="year")

    #incoming parameters from Climate
    rtl_realizedtemperature = Parameter(index=[time, region], unit="degreeC")

    #tolerability variables
    atl_adjustedtolerableleveloftemprise = Parameter(index=[time,region], unit="degreeC")
    imp_actualreduction = Parameter(index=[time, region], unit= "%")
    i_regionalimpact = Variable(index=[time, region], unit="degreeC")

    #impact Parameters
    rcons_per_cap_SLRRemainConsumption = Parameter(index=[time, region], unit = "\$/person")
    rgdp_per_cap_SLRRemainGDP = Parameter(index=[time, region], unit = "\$/person")

    save_savingsrate = Parameter(unit= "%", default=15.)
    wincf_weightsfactor_market =Parameter(index=[region], unit="")
    W_MarketImpactsatCalibrationTemp =Parameter(unit="%GDP", default=0.6)
    ipow_MarketIncomeFxnExponent =Parameter(default=-0.13333333333333333)
    iben_MarketInitialBenefit=Parameter(default=.1333333333333)
    tcal_CalibrationTemp = Parameter(default=3.)
    GDP_per_cap_focus_0_FocusRegionEU = Parameter(unit="\$/person", default=34298.93698672955)

    #impact variables
    isatg_impactfxnsaturation = Parameter(unit="unitless")
    rcons_per_cap_MarketRemainConsumption = Variable(index=[time, region], unit = "\$/person")
    rgdp_per_cap_MarketRemainGDP = Variable(index=[time, region], unit = "\$/person")
    iref_ImpactatReferenceGDPperCap=Variable(index=[time, region])
    igdp_ImpactatActualGDPperCap=Variable(index=[time, region])
    impmax_maxtempriseforadaptpolicyM = Parameter(index=[region], unit= "degreeC")

    isat_ImpactinclSaturationandAdaptation= Variable(index=[time,region])
    isat_per_cap_ImpactperCapinclSaturationandAdaptation = Variable(index=[time,region])
    pow_MarketImpactExponent=Parameter(unit="", default=2.1666666666666665)

    function run_timestep(p, v, d, t)

        for r in d.region

            #calculate tolerability
            if (p.rtl_realizedtemperature[t,r]-p.atl_adjustedtolerableleveloftemprise[t,r]) < 0
                v.i_regionalimpact[t,r] = 0
            else
                v.i_regionalimpact[t,r] = p.rtl_realizedtemperature[t,r]-p.atl_adjustedtolerableleveloftemprise[t,r]
            end

            v.iref_ImpactatReferenceGDPperCap[t,r]= p.wincf_weightsfactor_market[r]*((p.W_MarketImpactsatCalibrationTemp + p.iben_MarketInitialBenefit * p.tcal_CalibrationTemp)*
                (v.i_regionalimpact[t,r]/p.tcal_CalibrationTemp)^p.pow_MarketImpactExponent - v.i_regionalimpact[t,r] * p.iben_MarketInitialBenefit)

            v.igdp_ImpactatActualGDPperCap[t,r]= v.iref_ImpactatReferenceGDPperCap[t,r]*
                (p.rgdp_per_cap_SLRRemainGDP[t,r]/p.GDP_per_cap_focus_0_FocusRegionEU)^p.ipow_MarketIncomeFxnExponent

            if v.igdp_ImpactatActualGDPperCap[t,r] < p.isatg_impactfxnsaturation
                v.isat_ImpactinclSaturationandAdaptation[t,r] = v.igdp_ImpactatActualGDPperCap[t,r]
            else
                v.isat_ImpactinclSaturationandAdaptation[t,r] = p.isatg_impactfxnsaturation+
                    ((100-p.save_savingsrate)-p.isatg_impactfxnsaturation)*
                    ((v.igdp_ImpactatActualGDPperCap[t,r]-p.isatg_impactfxnsaturation)/
                    (((100-p.save_savingsrate)-p.isatg_impactfxnsaturation)+
                    (v.igdp_ImpactatActualGDPperCap[t,r]-
                    p.isatg_impactfxnsaturation)))
                end

            if v.i_regionalimpact[t,r] < p.impmax_maxtempriseforadaptpolicyM[r]
                v.isat_ImpactinclSaturationandAdaptation[t,r]=v.isat_ImpactinclSaturationandAdaptation[t,r]*(1-p.imp_actualreduction[t,r]/100)
            else
                v.isat_ImpactinclSaturationandAdaptation[t,r] = v.isat_ImpactinclSaturationandAdaptation[t,r] *
                    (1-(p.imp_actualreduction[t,r]/100)* p.impmax_maxtempriseforadaptpolicyM[r] /
                    v.i_regionalimpact[t,r])
            end

            v.isat_per_cap_ImpactperCapinclSaturationandAdaptation[t,r] = (v.isat_ImpactinclSaturationandAdaptation[t,r]/100)*p.rgdp_per_cap_SLRRemainGDP[t,r]
            v.rcons_per_cap_MarketRemainConsumption[t,r] = p.rcons_per_cap_SLRRemainConsumption[t,r] - v.isat_per_cap_ImpactperCapinclSaturationandAdaptation[t,r]
            v.rgdp_per_cap_MarketRemainGDP[t,r] = v.rcons_per_cap_MarketRemainConsumption[t,r]/(1-p.save_savingsrate/100)
        end

    end
end

# Still need this function in order to set the parameters than depend on
# readpagedata, which takes model as an input. These cannot be set using
# the default keyword arg for now.

function addmarketdamages(model::Model)
    marketdamagescomp = add_comp!(model, MarketDamages)
    marketdamagescomp[:impmax_maxtempriseforadaptpolicyM] = readpagedata(model, "data/impmax_economic.csv")

    # fix the current bug which implements the regional weights from SLR and discontinuity also for market and non-market damages (where weights should be uniformly one)
    marketdamagescomp[:wincf_weightsfactor_market] = readpagedata(model, "data/wincf_weightsfactor_market.csv")

    return marketdamagescomp
end

max_stars_repo_name: lrennels/PAGE-ICE

max_issues_repo_name: lrennels/PAGE-ICE

max_forks_repo_name: lrennels/PAGE-ICE


text:
module SequenceInference

# Write your package code here.
greet()=print("hello world")

end

max_stars_repo_name: guchchcug/SequenceInference

max_issues_repo_name: guchchcug/SequenceInference

max_forks_repo_name: guchchcug/SequenceInference


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

JLLWrappers.@generate_main_file_header("UCX")
JLLWrappers.@generate_main_file("UCX", UUID("16e4e860-d6b8-5056-a518-93e88b6392ae"))
end  # module UCX_jll