text:
#import Pkg
#Pkg.add("StatsBase")
using StatsBase

function find_cycle(g, chain)
    head = chain[end]
    if length(chain) > 1
        last = chain[end-1]
        neis = []
        for i in g[head]
            if i != last
                push!(neis, i)
            end
        end
    else
        neis = collect(g[head])
    end
    while true
        next = sample(neis)
        # println(next)
        if in(next, chain)
            # rec = indexin(next, chain)
            rec=0
            for i in 1:length(chain)
                if chain[i] == next
                    rec = i
                    break
                end
            end
            cycle = collect(view(chain, rec:length(chain)))
            tail  = collect(view(chain, 1:rec-1))
            return tail, cycle
        end
        push!(chain, next)
        last = chain[end-1]
        neis = []
        for i in g[next]
            if i != last
                push!(neis, i)
            end
        end
    end
end
# 書けたっぽい!

function tileByCycles(g)
    cycles = Array([])
    chain  = Array([])
    while length(g) > 0
        if length(chain) == 0
            L = collect(keys(g))
            head = sample(L)
            chain = Array([head])
        end
        # println("chain",chain)
        chain, cycle = find_cycle(g, chain)
        # println(cycle)
        push!(cycles, cycle)
        for i in 1:length(cycle)
            ii = i - 1
            if ii == 0
                ii = length(cycle)
            end
            a = cycle[ii]
            b = cycle[i]
            # println("remove edge",a,b)
            # remove edge from g
            pop!(g[a], b)
            pop!(g[b], a)
        end
        for a in cycle
            if length(g[a]) == 0
                pop!(g, a)
            end
        end
        # println(g)
    end
    return cycles
end

# どうにか書いたけど、Juliaは完全なコンパイラ言語ではなくJITだったらしい。うーむむむむ。
# このコードを、辞書を使わずに高速に書けるか???
# 今日はあきらめる。

g = Dict(1=>Set([2,3,4,5]),
         2=>Set([1,3,4,6]),
         3=>Set([1,2,5,6]),
         4=>Set([1,2,5,6]),
         5=>Set([1,3,4,6]),
         6=>Set([2,3,4,5]))

# chain = Array([1,2])
# tail, cycle = find_cycle(g, chain)
# println("tail")
# println(tail)
# println("cycle")
# println(cycle)

cycles = tileByCycles(g)
println(cycles)

max_stars_repo_name: vitroid/TileCycles

max_issues_repo_name: vitroid/TileCycles

max_forks_repo_name: vitroid/TileCycles


text:
# Copyright (c) 2021
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

using LinearAlgebra, Convex, SCS;

"""
    IRLS(A::Matrix{Float64}, b::Vector{Float64}; maxiter=100, epsilon=0.01)

Find the solution to Ax=b using Iterative Recursive Least Squares.

### Input

- `A`       -- Matrix: Ax=b
- `b`       -- Vector: Ax=b
- `maxiter` -- (optional) number of optmization iterations
- `epsilon` -- (optional) threshold to stop optimizing

### Output

Solution to Ax=b (Vector{Float64})

### Example

julia> A = randn(10, 100);
julia> b = randn(10);
julia> x = IRLS(A, b, maxiter=5, epsilon=.01);
julia> println(x)

### Algorithm

This function implements Andrew's monotone chain convex hull algorithm to
Enhancing Sparsity by Reweighted L1 Minimization

    Emmanuel J. Candes, Michael B. Wakin, and Stephen P. Boyd, "Enhancing
        Sparsity by Reweighted L1 Minimization," J Fourier Anal Appl (2008)
        14: 877–905.

"""
function IRWLS(A::Matrix{Float64},
               b::Vector{Float64};
               maxiter=100,
               epsilon=0.01)
    # get the second column dim of A
    local p, x, xhat, w;
    _, p = size(A);
    w = ones(p);
    w_old = w;
    o = ones(p);
    x = zeros(p);

    solver = () -> SCS.Optimizer(verbose=0)

    for i = 1:maxiter
        # create a diagonal matrix from the vector w
        W = diagm(w);

        # solve the convex optimization task
        xhat = Variable(p);
        prob = minimize(norm(W*xhat, 1), A*xhat==b)
        solve!(prob, solver)
        xhat = evaluate(xhat);

        # check the stop condition then move on
        w = 1.0./(epsilon*o + abs.(xhat))
        if norm(w - w_old, 2) <= 1e-6
            break;
        end
        w = w_old;
    end
    x = xhat;
    i = abs.(x) .< epsilon;
    x[i] = zeros(sum(i));
    return x;
end

# A = randn(10, 50); b = randn(10);
# x = IRWLS(A, b);
# println(x)

ext: jl


text:
export ExtensionCost_Constant,
    ExtensionCost_Zero

# -------------------------------------------
# Extension
# -------------------------------------------

@inline dflow!(::AbstractExtensionCost, args...) = nothing
@inline flowdψ(::AbstractExtensionCost, args...) = 0.0

"No extension cost"
struct ExtensionCost_Zero <: AbstractExtensionCost end
_nparm(::ExtensionCost_Zero) = 0
extensioncost(::ExtensionCost_Zero, θ) = azero(θ)
@inline function flow!(grad, x::ExtensionCost_Zero, d, obs, θ, sim, dograd::Bool)
    return extensioncost(x,θ)
end
coefnames(::ExtensionCost_Zero) = Vector{String}(undef,0)

"Constant extension cost"
struct ExtensionCost_Constant <: AbstractExtensionCost end
_nparm(::ExtensionCost_Constant) = 1
extensioncost(::ExtensionCost_Constant, θ) = θ[1]
@inline function flow!(grad, x::ExtensionCost_Constant, d, obs, θ, sim, dograd::Bool)
    sgn = _sgnext(d, obs)
    if dograd
        grad[1] = sgn
    end
    return sgn ? extensioncost(x,θ) : azero(θ)
end
coefnames(::ExtensionCost_Constant) = ["\\alpha_{ext}",]

# -------------------------------------------
# Scrap value
# -------------------------------------------

export ScrapValue_Zero, ScrapValue_Constant, ScrapValue_Price

@inline dflow!(::AbstractScrapValue, args...) = nothing
@inline flowdψ(::AbstractScrapValue, args...) = 0.0

"No scrap value"
struct ScrapValue_Zero <: AbstractScrapValue end
_nparm(::ScrapValue_Zero) = 0
@inline function flow!(grad, x::ScrapValue_Zero, d, obs, θ, sim, dograd::Bool)
    return azero(θ)
end
coefnames(::ScrapValue_Zero) = Vector{String}(undef,0)

"Constant scrap value"
struct ScrapValue_Constant <: AbstractScrapValue end
_nparm(::ScrapValue_Constant) = 1
@inline function flow!(grad, x::ScrapValue_Constant, d, obs, θ, sim, dograd::Bool)
    e = expires_today(d, obs)
    if dograd
        grad[1] = e
    end
    return e ? θ[1] : azero(θ)
end
coefnames(::ScrapValue_Constant) = ["\\alpha_{scrap}",]

"Scrap value with price"
struct ScrapValue_Price <: AbstractScrapValue end
_nparm(::ScrapValue_Price) = 2
@inline function flow!(grad, x::ScrapValue_Price, d, obs, θ, sim, dograd::Bool)
    e = expires_today(d, obs)
    p = exp(logprice(obs))
    if dograd
        grad[1] = e
        grad[2] = e*p
    end
    return e ? θ[1] + θ[2]*p : azero(θ)
end
coefnames(::ScrapValue_Price) = ["\\alpha_{scrap,0}", "\\alpha_{scrap,p}"]


text:
struct CoefNotComparable
end
Base.iszero(::CoefNotComparable) = false

@testset "Term" begin
    Mod.@polyvar x
    @test convert(Any, 1x) == 1x
    @test one(1x) == one(1.0x) == 1
    @test zero(1x) == zero(1.0x) == 0
    @test nvariables(0.0x) == 1
    @test nvariables(1x) == 1
    #@inferred one(1x)
    @inferred zero(1x)
    #@inferred one(1.0x)
    @inferred zero(1.0x)

    @test !isapproxzero((1+1e-10im)*x^2)
    @test isapproxzero(1e-10im*x^2)

    @test monic(2.0x) isa AbstractTerm{Float64}
    @test monic(2.0x) == x
    @test monic(2x^2) isa AbstractTerm{Int}
    @test monic(2x^2) == x^2

    @test leadingterm(2x^2) == 2x^2
    @test nterms(2x^2) == 1
    @test terms(2x^2) == [2x^2]
    @test nterms(0*x) == 0
    @test terms(0*x) == typeof(0*x)[]
    @test nterms(0.0x) == 0
    @test terms(0.0x) == typeof(0.0x)[]

    @test nterms(polynomial(0.0x)) == 0
    @test nterms(convert(polynomialtype(0.0x), 0.0x)) == 0

    @test term(x) isa AbstractTerm
    @test term(x^2) == x^2
    @test term(1x^2) isa AbstractTerm
    @test term(1x) == x
    @test zeroterm(1x) == 0*x

    Mod.@polyvar y
    @test degree(2x^2, x) == 2
    @test degree(2x^2, y) == 0
    @test degree(2x^2, y) == 0

    t = 3x^2*y^4
    typetests(t)
    typetests([t, 2x])
    @test (@inferred polynomial(t)) isa AbstractPolynomial{Int}
    @test (@inferred polynomial(t, Float64)) isa AbstractPolynomial{Float64}

    @test_throws InexactError push!([1], 2x)
    @test_throws InexactError push!([x^2], 2x)

    @testset "Effective variables" begin
        T = variable_union_type(x)
        @test x isa T
        @test y isa T
        @test T[x, y] == @inferred effective_variables(x * y)
        @test T[x, y] == @inferred effective_variables(y * x)
        @test T[x] == @inferred effective_variables(x * y^0)
        @test T[x] == @inferred effective_variables(y^0 * x)
        @test T[y] == @inferred effective_variables(x^0 * y)
        @test T[y] == @inferred effective_variables(y * x^0)
    end

    @testset "Compare terms" begin
        t1 = 1 * x
        t2 = 2 * x
        @test t1 < t2
        @test t1 <= t2
        @test !(t1 > t2)
        @test !(t1 >= t2)

        t1 = (1 + 1im) * x
        t2 = (2 + 2im) * x
        @test t1 < t2
        @test t1 <= t2
        @test !(t1 > t2)
        @test !(t1 >= t2)

        a = CoefNotComparable()
        b = CoefNotComparable()
        t1 = a * x
        t2 = b * y
        @test t1 > t2
        @test t1 >= t2
        @test !(t1 < t2)
        @test !(t1 <= t2)
        t1 = a * x
        t2 = b * x
        @test !(t1 > t2)
        @test t1 >= t2
        @test !(t1 < t2)
        @test t1 <= t2
    end
end


text:
# Slugify.jl -- A library that simplifies a text to an ASCII subset
# By: Emmanuel Raviart <emmanuel@raviart.com>
#
# Copyright (C) 2015 Emmanuel Raviart
# https://github.com/eraviart/Slugify.jl
#
# This file is part of Slugify.jl.
#
# The Slugify.jl package is licensed under the MIT "Expat" License.

using Base.Test

importall Slugify

# write your own tests here
@test slugify("Hello world!") == "hello-world"
@test slugify("   Hello   world!   ") == "hello-world"
@test slugify("œil, forêt, ça, où...") == "oeil-foret-ca-ou"
@test slugify("   ") == ""
@test slugify("") == ""
@test slugify("   Hello   world!   ", separator = ' ', transform = uppercase) == "HELLO WORLD"
@test slugify("   ", separator = ' ') == ""
@test slugify("", separator = ' ') == ""


text:
using Interpolations
using FFTW

"""
    Sheets.redistribute_points!(sheet, zs, Γs)

Returns the modified sheet with replacement control points at positions `zs` and strength `Γs`.

```jldoctest
julia> sheet = Vortex.Sheet(0:0.1:1, 0.0:10, 0.2)
Vortex Sheet: L ≈ 1.0, Γ = 10.0, δ = 0.2

julia> sys = (sheet,)
(Vortex Sheet: L ≈ 1.0, Γ = 10.0, δ = 0.2,)

julia> Sheets.redistribute_points!(sheet, 0:0.2:2, 0.0:0.5:5)
Vortex Sheet: L ≈ 2.0, Γ = 5.0, δ = 0.2

julia> sys[1]
Vortex Sheet: L ≈ 2.0, Γ = 5.0, δ = 0.2
```
"""
function redistribute_points!(sheet::Sheet{T}, zs, Ss) where T
    if !(sheet.Ss === Ss)
        resize!(sheet.Ss, length(Ss))
        copy!(sheet.Ss, Ss)
    end

    sheet.blobs = Blob{T}.(zs, compute_trapezoidal_weights(Ss), sheet.δ)
    sheet.zs = mappedarray(Elements.position, sheet.blobs)
    sheet
end

"""
    Sheets.split!(sheet, n::Int)

Remove segments `0:n` from `sheet`, and return those segments as a new sheet.

# Example

```jldoctest
julia> sheet = Vortex.Sheet(0:0.1:1, 0.0:10, 0.2)
Vortex Sheet: L ≈ 1.0, Γ = 10.0, δ = 0.2

julia> sheet₋ = Sheets.split!(sheet, 5)
Vortex Sheet: L ≈ 0.4, Γ = 4.0, δ = 0.2

julia> sheet
Vortex Sheet: L ≈ 0.6, Γ = 6.0, δ = 0.2
```
"""
function split!(sheet::Sheet{T}, n::Int) where T
    @assert 2 < n < length(sheet) - 2

    blobs = splice!(sheet.blobs, 1:n-1)
    Ss = splice!(sheet.Ss, 1:n-1)
    push!(Ss, sheet.Ss[1])
    push!(blobs, Blob{T}(sheet.blobs[1].z, 0.5(Ss[end] - Ss[end-1]), sheet.δ))

    sheet.blobs[1] = Blob{T}(sheet.blobs[1].z, 0.5(sheet.Ss[2] - sheet.Ss[1]), sheet.δ)

    Sheet(blobs, Ss, sheet.δ)
end

"""
    Sheets.truncate!(sheet, n::Int)

Remove segments `0:n` from `sheet`, and return the circulation in those segments.

# Example

```jldoctest
julia> sheet = Vortex.Sheet(0:0.1:1, 0.0:10, 0.2)
Vortex Sheet: L ≈ 1.0, Γ = 10.0, δ = 0.2

julia> Sheets.truncate!(sheet, 5)
4.0
```
"""
function truncate!(sheet::Sheet{T}, n::Int) where T
    @assert 2 ≤ n ≤ length(sheet)-1
    ΔS = sheet.Ss[n] - sheet.Ss[1]
    splice!(sheet.Ss, 1:n-1)
    splice!(sheet.blobs, 1:n-1)
    sheet.blobs[1] = Blob{T}(sheet.blobs[1].z, 0.5(sheet.Ss[2] - sheet.Ss[1]), sheet.δ)
    return ΔS
end

"""
    Sheets.remesh(sheet, Δs::Float64 , params::Tuple = ())

Uniformly redistribute the control points of the sheet to have a nominal spacing of `Δs`.
Material quantities that should be redistributed along with the control points can be passed in as elements of `params`.

Returns the tuple `(z₌, Γ₌, L [, p₌])` where

- `z₌` is an array with the positions of the uniformly distributed points
- `Γ₌` is circulation interpolated onto `z₌`
- `L` is total length of the sheet
- `p₌` is a tuple containing the material quantities from `params` interpolated onto `z₌`

# Example

```jldoctest
julia> sheet = Vortex.Sheet(0:0.1:1, 0.0:10, 0.2)
Vortex Sheet: L ≈ 1.0, Γ = 10.0, δ = 0.2

julia> age = collect(10.0:-1:0);

julia> Sheets.remesh(sheet, 0.2, (age, ))
(Complex{Float64}[0.0+0.0im, 0.25+0.0im, 0.5+0.0im, 0.75+0.0im, 1.0+0.0im], [0.0, 2.5, 5.0, 7.5, 10.0], 1.0, ([10.0, 7.5, 5.0, 2.5, 0.0],))
```
"""
function remesh(sheet::Sheet{S}, Δs::Float64, params::Tuple = ()) where S
    L = arclengths(sheet)

    if L[end] < Δs
        @warn("Cannot remesh, sheet length smaller than nominal spacing")
        return Elements.position.(sheet.blobs), sheet.Ss, L[end], params
    end

    knots = (L,)
    mode = Gridded(Linear())

    zspline = interpolate(knots, sheet.zs, mode)
    Γspline = interpolate(knots, sheet.Ss, mode)

    psplines = map(params) do p
        interpolate(knots, p, mode)
    end

    N = ceil(Int, L[end]/Δs)

    L₌ = range(0, L[end], length=N)

    z₌ = zspline(L₌)
    S₌ = Γspline(L₌)

    p₌ = map(psplines) do spline
        spline(L₌)
    end

    z₌, S₌, L[end], p₌
end

"""
    Sheets.remesh!(sheet::Sheet, Δs::Float64, params::Tuple = ())

Same as [`Sheets.remesh`](@ref), except `sheet` is replaced
internally by a uniformly interpolated control points.
Returns the tuple (sheet, L, p₌) where

- `sheet` is the modified sheet
- `L` is total length of the sheet
- `p₌` is a tuple containing the material quantities from `params` interpolated onto the new control points of `sheet`

```jldoctest
julia> sheet = Vortex.Sheet(0:0.1:1, 0.0:10, 0.2)
Vortex Sheet: L ≈ 1.0, Γ = 10.0, δ = 0.2

julia> age = collect(10.0:-1:0);

julia> Sheets.remesh!(sheet, 0.2, (age,));

julia> Elements.position.(sheet.blobs)
5-element Array{Complex{Float64},1}:
  0.0 + 0.0im
 0.25 + 0.0im
  0.5 + 0.0im
 0.75 + 0.0im
  1.0 + 0.0im

julia> age
5-element Array{Float64,1}:
 10.0
  7.5
  5.0
  2.5
  0.0
```
"""
function remesh!(sheet::Sheet{S}, Δs::Float64, params::Tuple = ()) where S
    z₌, S₌, L, p₌ = remesh(sheet, Δs, params)
    if L < Δs
        return sheet, L, params
    end
    redistribute_points!(sheet, z₌, S₌)

    for i in 1:length(params)
        resize!(params[i], length(p₌[i]))
        copy!(params[i], p₌[i])
    end

    sheet, L, params
end

"""
    Sheets.arclength(s)

Compute the polygonal arc length of `s`, where `s` can be either an
vector of complex numbers or a `Vortex.Sheet`.

# Example

```jldoctest
julia> sheet = Vortex.Sheet(0:0.1:1, 0.0:10, 0.2)
Vortex Sheet: L ≈ 1.0, Γ = 10.0, δ = 0.2

julia> Sheets.arclength(sheet)
1.0
"""
function arclength(zs::AbstractVector)
    L = 0.0
    for i in 2:length(zs)
        L += abs(zs[i] - zs[i-1])
    end
    L
end
arclength(sheet::Sheet) = arclength(sheet.zs)

"""
    Sheets.arclengths(s)

Cumulative sum of the polygonal arc length of `s`, where `s` can be
either an vector of complex numbers or a `Vortex.Sheet`.

# Example

```jldoctest
julia> sheet = Vortex.Sheet(0:0.1:1, 0.0:10, 0.2)
Vortex Sheet: L ≈ 1.0, Γ = 10.0, δ = 0.2

julia> Sheets.arclengths(sheet)
11-element Array{Float64,1}:
 0.0
 0.1
 0.2
 0.3
 0.4
 0.5
 0.6
 0.7
 0.8
 0.9
 1.0
```
"""
function arclengths(zs::AbstractVector)
    N = length(zs)

    Δs = zeros(N)
    for i in 2:N
        Δs[i] = abs(zs[i] - zs[i-1])
    end
    L = accumulate(+, Δs)
end
arclengths(sheet::Sheet) = arclengths(sheet.zs)

"""
    Sheets.append_segment!(sheet::Sheet, z, Γ)

Append a new segment with circulation `Γ` extending from the end of the sheet to `z`.

# Example

```jldoctest
julia> sheet = Vortex.Sheet(0:0.1:1, 0.0:10, 0.2)
Vortex Sheet: L ≈ 1.0, Γ = 10.0, δ = 0.2

julia> sheet.blobs[end]
Vortex.Blob(1.0 + 0.0im, 0.5, 0.2)

julia> Sheets.append_segment!(sheet, 1.1, 2.0)

julia> sheet
Vortex Sheet: L ≈ 1.1, Γ = 12.0, δ = 0.2

julia> sheet.blobs[end]
Vortex.Blob(1.1 + 0.0im, 1.0, 0.2)
```
"""
function append_segment!(sheet::Sheet{T}, z, S::T) where T
    b₋ = sheet.blobs[end]
    sheet.blobs[end] = Blob{T}(b₋.z, b₋.S + 0.5S, b₋.δ)
    push!(sheet.blobs, Blob{T}(z, 0.5S, sheet.δ))
    push!(sheet.Ss, sheet.Ss[end] + S)
    nothing
end

"""
    Sheets.filter!(sheet, Δs, Δf[, params])

Redistribute and filter the control points of a vortex sheet

# Arguments

- `sheet`: the vortex sheet to be modified
- `Δs`: the nominal spacing between the uniform points
- `Δf`: the minimum length scale that the filter should allow to pass
  through
- `params`: an optional tuple of vectors containing material properties

# Returns

If `params` is passed in, then its vectors will be overwritten by
their interpolated values on the new control points, and the function
returns the tuple (sheet, params).
Otherwise, it returns (sheet, ())
"""
function filter!(sheet, Δs, Δf, params::Tuple = ())
    z₌, S₌, L, p₌ = remesh(sheet, Δs, params)

    for i in 1:length(params)
        resize!(params[i], length(p₌[i]))
        copy!(params[i], p₌[i])
    end

    if L > Δs
        filter_position!(z₌, Δf, L)
        redistribute_points!(sheet, z₌, S₌), params
    else
        @warn("Filter not applied, total sheet length smaller than nominal spacing")
        sheet, params
    end
end

"""
    Sheets.filter_position!(s, Δf, L = arclength(z₌))

Filter out any length scales in `s` that is smaller than `Δf`, storing the result back in `s`.
`s` can be either a vector of complex positions, or a `Vortex.Sheet`.
"""
function filter_position!(z₌::AbstractVector{T}, Δf, L = arclength(z₌)) where T
    Δs = abs(z₌[2] - z₌[1])

    @assert Δs < Δf

    F = plan_dct!(z₌)
    F * z₌

    cutoff = ceil(Int, 2L/Δf) + 1

    z₌[cutoff:end] .= zero(ComplexF64)

    # F \ z₌
    # Out-of-date syntax in FFTW for planned inverse transforms
    # manually construct inverse transform for now
    F⁻¹ = let X = Array{T}(undef, F.plan.sz),
             iK = FFTW.inv_kind[FFTW.REDFT10]
           FFTW.DCTPlan{T,iK,true}(FFTW.plan_r2r!(X, iK, F.region, flags=F.plan.flags),
                              F.r, F.nrm, F.region)
       end

    F⁻¹ * z₌
end

function filter_position!(sheet::Sheet, Δf, L = arclength(sheet))
    zs = Elements.position.(sheet.blobs)
    filter_position!(zs, Δf, L)

    redistribute_points!(sheet, zs, sheet.Ss)
end


text:
function get_eop_history(obj::EOPBasedTransformProvider)
    return jcall(obj, "getEOPHistory", EOPHistory, ())
end

function get_non_interpolating_provider(obj::EOPBasedTransformProvider)
    return jcall(obj, "getNonInterpolatingProvider", EOPBasedTransformProvider, ())
end

function get_transform(obj::TransformProvider, arg0::AbsoluteDate)
    return jcall(obj, "getTransform", Transform, (AbsoluteDate,), arg0)
end

function get_transform(obj::TransformProvider, arg0::FieldAbsoluteDate)
    return jcall(obj, "getTransform", FieldTransform, (FieldAbsoluteDate,), arg0)
end


text:
using Random
using LinearAlgebra
using SparseArrays
using GraphRecipes
using Plots

cd(@__DIR__)
cd("../assets")

n = 15
Random.seed!(1)
A = Float64[ rand() < 0.5 ? 0 : rand() for i=1:n, j=1:n]
for i=1:n
    A[i, 1:i-1] = A[1:i-1, i]
end
x = rand(n)
y = rand(n)
z = rand(n)
graphplot(A,
          node_weights = 1:n,
          nodecolor = range(colorant"yellow", stop=colorant"red", length=n),
          names = 1:n,
          linecolor = :darkgrey,
          layout_kw = Dict(:x => x, :y => y),
          )
savefig("random_labelled_graph.png")

graphplot(A,
          node_weights = 1:n,
          markercolor = :darkgray,
          dim = 3,
          markersize = 5,
          linecolor = :darkgrey,
          linealpha = 0.5,
          layout_kw = Dict(:x => x, :y => y, :z => z),
          )
savefig("random_3d_graph.png")

adjmat = Symmetric(sparse(rand(0:1,8,8)))

adjmat = Symmetric(sparse(rand(0:1,8,8)))
plot(
     graphplot(adjmat,
               method=:chorddiagram,
               names=[text(string(i), 8) for i in 1:8],
               linecolor=:black,
               fillcolor=:lightgray),

     graphplot(adjmat,
               method=:arcdiagram,
               markersize=3,
               linecolor=:black,
               markercolor=:black)
               )
savefig("arc_chord_diagrams.png")

code = :(
function mysum(list)
    out = 0
    for value in list
        out += value
    end
    out
end
)

default(size=(1000, 1000))
plot(code, fontsize=10, shorten=0.01, axis_buffer=0.15, nodeshape=:rect)
savefig("AST_example.png")

plot(AbstractFloat, method=:tree, fontsize=10, nodeshape=:ellipse)
savefig("julia_type_tree.png")

# `AbstractTrees` example
using AbstractTrees
AbstractTrees.children(d::Dict) = [p for p in d]
AbstractTrees.children(p::Pair) = AbstractTrees.children(p[2])
function AbstractTrees.printnode(io::IO, p::Pair)
    str = isempty(AbstractTrees.children(p[2])) ? string(p[1], ": ", p[2]) : string(p[1], ": ")
    print(io, str)
end

d = Dict(:a => 2,:d => Dict(:b => 4,:c => "Hello"),:e => 5.0)

plot(TreePlot(d), method=:tree, fontsize=10, nodeshape=:ellipse)
savefig("julia_dict_tree.png")


text:
function test_euler_all()
    print(rpad("Math.Series.e...", 50, ' '))

    simple = Brainstorm._Math._Series._Euler._Taylor._Simple.euler
    inplace = Brainstorm._Math._Series._Euler._Taylor._Inplace.euler

    setprecision(4_000_000) do
        expected = [trunc(BigInt, e*(big(10)^d)) for d in (1, 10, 100, 1_000, 10_000, 100_000, 1_000_000)]

        @test all([simple(d)  for d in (1, 10, 100, 1_000, 10_000, 100_000, 1_000_000)] - expected .== 0)
        @test all([inplace(d) for d in (1, 10, 100, 1_000, 10_000, 100_000, 1_000_000)] - expected .== 0)
    end

    println("PASS")
end


text:
# indicator of nonnegative orthant

export IndNonnegative

"""
**Indicator of the nonnegative orthant**

    IndNonnegative()

Returns the indicator of the set
```math
C = \\{ x : x \\geq 0 \\}.
```
"""
struct IndNonnegative <: ProximableFunction end

is_separable(f::IndNonnegative) = true
is_convex(f::IndNonnegative) = true
is_cone(f::IndNonnegative) = true

function (f::IndNonnegative)(x::AbstractArray{R}) where R <: Real
    for k in eachindex(x)
        if x[k] < 0
            return R(Inf)
        end
    end
    return R(0)
end

function prox!(y::AbstractArray{R}, f::IndNonnegative, x::AbstractArray{R}, gamma=R(1)) where R <: Real
    for k in eachindex(x)
        if x[k] < 0
            y[k] = R(0)
        else
            y[k] = x[k]
        end
    end
    return R(0)
end

fun_name(f::IndNonnegative) = "indicator of the Nonnegative cone"
fun_dom(f::IndNonnegative) = "AbstractArray{Real}"
fun_expr(f::IndNonnegative) = "x ↦ 0 if all(0 ⩽ x), +∞ otherwise"
fun_params(f::IndNonnegative) = "none"

function prox_naive(f::IndNonnegative, x::AbstractArray{R}, gamma=R(1)) where R <: Real
    y = max.(R(0), x)
    return y, R(0)
end


text:
module Eurocode1993

export eq13_64_2006, eq13_68_2006, eq13_618_2006, eq13_627_2006, eq13_628c_2006

function eq13_64_2006(Weff, fyb, γMO)

    McRd = Weff*fyb/γMO

end

function eq13_68_2006(hw, sw, ϕ, t, fyb, E, γMO)

    λw = 0.346 * sw/t * sqrt(fyb/E)

   if λw <= 0.83
        fbv = 0.58 * fyb
   elseif (λw > 0.83) & (λw <= 1.40)
        fbv = 0.48 * fyb / λw
   elseif λw > 1.40
        fbv = 0.67 * fyb / λw^2
   end

   VbRd = (hw/sin(deg2rad(ϕ)) *t * fbv) / γMO

   return VbRd, fbv, λw

end

function eq13_618_2006(α, t, fyb, E, r, la, ϕ, γM1)

    RwRd = α * t^2 * sqrt(fyb * E) * (1 - 0.1 * sqrt(r/t)) * (0.5 + sqrt(0.02 * la / t)) * (2.4 + (ϕ/90)^2) / γM1

end

function eq13_627_2006(NEd, NRd, MyEd, MyRd, MfRd, MplRd, VEd, VwRd)

     I_N = abs(NEd/NRd)
     I_M = abs(MyEd/MyRd)
     I_V = abs((1 - MfRd/MplRd)*(2*VEd/VwRd - 1)^2)
     I_total = I_N + I_M + I_V

     return I_N, I_M, I_V, I_total
end

function eq13_628c_2006(MEd, McRd, FEd, RwRd)

     I_M = abs(MEd/McRd)
     I_F = abs(FEd/RwRd)
     I_total = I_M + I_F

     return I_M, I_F, I_total

end

end #module


text:
using Meshing
using Test
using GeometryTypes
import GeometryBasics
const GB = GeometryBasics
using ForwardDiff
using StaticArrays
using Statistics: mean
using LinearAlgebra: dot, norm

algos = (MarchingCubes, MarchingTetrahedra, NaiveSurfaceNets)

@testset "meshing algorithms" begin
    for algo in (MarchingCubes, MarchingTetrahedra, NaiveSurfaceNets)
        @test algo(5) == algo{Float64}(5.0, 0.001, true, false)
        @test algo(0x00, 0x00) == algo{UInt8}(0x00, 0x00, true, false)
        @test algo{UInt16}(5, 0x00) == algo{UInt16}(0x0005, 0x0000, true, false)
        @test algo{Float32}(1) ==  algo{Float32}(1.0f0, 0.001f0, true, false)
        @test algo{Float32}(1.0, 0x00) == algo{Float32}(1.0f0, 0.0f0, true, false)
        @test algo{Float32}(iso=1, eps=0x00, insidepositive=true, reduceverts=false) ==
            algo{Float32}(1.0f0, 0.0f0, false, true)
    end
end

@testset "type helpers" begin
    dt = Meshing._determine_types
    @test dt(HomogenousMesh) == (Point{3,Float64}, Face{3,Int64})
    @test dt(GLNormalMesh) == (Point{3,Float32}, Face{3,OffsetInteger{-1,UInt32}})
    @test dt(HomogenousMesh, Float16) == (Point{3,Float16}, Face{3,Int64})
    @test dt(HomogenousMesh{Point{3, Float64}, Face{3, UInt32}}) == (Point{3,Float64}, Face{3,UInt32})
    @test dt(HomogenousMesh{Point{3, Float32}, Face{3, UInt32}}, Float64) == (Point{3,Float32}, Face{3,UInt32})
    @test dt(HomogenousMesh, Float64, 4) == (Point{3,Float64}, Face{4,Int64})
end

@testset "surface nets" begin
    sdf_sphere = SignedDistanceField(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.))) do v
        sqrt(sum(dot(v,v))) - 1 # sphere
    end
    sdf_torus = SignedDistanceField(HyperRectangle(Vec(-2,-2,-2.),Vec(4,4,4.)), 0.05) do v
        (sqrt(v[1]^2+v[2]^2)-0.5)^2 + v[3]^2 - 0.25
    end

    snf_torus = HomogenousMesh(HyperRectangle(Vec(-2,-2,-2.),Vec(4,4,4.)), (81,81,81), NaiveSurfaceNets()) do v
        (sqrt(v[1]^2+v[2]^2)-0.5)^2 + v[3]^2 - 0.25
    end

    snf_sphere = HomogenousMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)), (21,21,21), NaiveSurfaceNets()) do v
        sqrt(sum(dot(v,v))) - 1 # sphere
    end

    # test convience constructors
    HomogenousMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)), NaiveSurfaceNets()) do v
        sqrt(sum(dot(v,v))) - 1 # sphere
    end
    HomogenousMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)), NaiveSurfaceNets(), samples=(5,5,5)) do v
        sqrt(sum(dot(v,v))) - 1 # sphere
    end

    sphere = HomogenousMesh(sdf_sphere, NaiveSurfaceNets())
    torus = HomogenousMesh(sdf_torus, NaiveSurfaceNets())
    @test length(vertices(sphere)) == 1832
    @test length(vertices(torus)) == 5532
    @test length(faces(sphere)) == 1830
    @test length(faces(torus)) == 5532

    @test length(vertices(sphere)) == length(vertices(snf_sphere))
    @test length(vertices(torus)) == length(vertices(snf_torus))
    @test length(faces(sphere)) == length(faces(snf_sphere))
    @test length(faces(torus)) == length(faces(snf_torus))
end

@testset "noisy spheres" begin
    # Produce a level set function that is a noisy version of the distance from
    # the origin (such that level sets are noisy spheres).
    #
    # The noise should exercise marching tetrahedra's ability to produce a water-
    # tight surface in all cases (unlike standard marching cubes).
    #
    N = 10
    sigma = 1.0
    distance = Float32[ sqrt(Float32(i*i+j*j+k*k)) for i = -N:N, j = -N:N, k = -N:N ]
    distance = distance + sigma*rand(2*N+1,2*N+1,2*N+1)

    # Extract an isosurface.
    #
    lambda = N-2*sigma # isovalue

    msh = HomogenousMesh(distance, MarchingTetrahedra(lambda))

    s2 = SignedDistanceField(HyperRectangle(Vec(0,0,0.),Vec(1,1,1.))) do v
        sqrt(sum(dot(v,v))) - 1 # sphere
    end

    msh = HomogenousMesh(s2, MarchingTetrahedra())
    @test length(vertices(msh)) == 973
    @test length(faces(msh)) == 1830
end

@testset "vertex interpolation" begin
    @test Meshing.vertex_interp(0, Vec(0,0,0), Vec(0,1,0), -1, 1) == Vec(0,0.5,0)
    @test Meshing.vertex_interp(-1, Vec(0,0,0), Vec(0,1,0), -1, 1) == Vec(0,0,0)
    @test Meshing.vertex_interp(1, Vec(0,0,0), Vec(0,1,0), -1, 1) == Vec(0,1,0)
end

@testset "marching cubes" begin

    algo = MarchingCubes()
    algo_pos = MarchingCubes(insidepositive=true)

    sdf = SignedDistanceField(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.))) do v
        sqrt(sum(dot(v,v))) - 1 # sphere
    end

    mf = SimpleMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)),algo, samples=(21,21,21)) do v
        sqrt(sum(dot(v,v))) - 1 # sphere
    end

    mf_pos = SimpleMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)),algo_pos, samples=(21,21,21)) do v
        -sqrt(sum(dot(v,v))) + 1 # sphere positive inside
    end

    mfrv = SimpleMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)),MarchingCubes(reduceverts=false), samples=(21,21,21)) do v
        sqrt(sum(dot(v,v))) - 1 # sphere
    end

    # convience constructors
    SimpleMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)), MarchingCubes()) do v
        sqrt(sum(dot(v,v))) - 1 # sphere
    end
    SimpleMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)), MarchingCubes(), samples=(5,6,7)) do v
        sqrt(sum(dot(v,v))) - 1 # sphere
    end

    @test length(vertices(mfrv)) == 10968
    m = SimpleMesh(sdf,algo)
    m2 = SimpleMesh(sdf,algo)
    @test length(vertices(m)) == 7320
    @test length(faces(m)) == 3656
    @test length(faces(mf)) == length(faces(mfrv))
    @test length(faces(mf)) == length(faces(mf_pos))
    @test m == m2
    @test length(vertices(m)) == length(vertices(mf))
    @test length(vertices(mf)) == length(vertices(mf_pos))
    @test length(faces(m)) == length(faces(mf))
end

@testset "marching tetrahedra" begin
    a1 = MarchingTetrahedra()
    a2 = MarchingTetrahedra(reduceverts=false)
    m1 = SimpleMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)), a1, samples=(21,21,21)) do v
        sqrt(sum(dot(v,v))) - 1 # sphere
    end
    m2 = SimpleMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)), a2, samples=(21,21,21)) do v
        sqrt(sum(dot(v,v))) - 1 # sphere
    end
    @test length(m1.vertices) == 5582
    @test length(m2.vertices) == 33480
    @test length(m1.faces) == 11160
    @test length(m2.faces) == length(m1.faces)
end

@testset "sign flips" begin
    for algo_type in algos
        algo = algo_type()
        algo_pos = algo_type(insidepositive=true)
        mf = SimpleMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)),algo, samples=(21,21,21)) do v
            sqrt(sum(dot(v,v))) - 1 # sphere
        end

        mf_pos = SimpleMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)),algo_pos, samples=(21,21,21)) do v
            -sqrt(sum(dot(v,v))) + 1 # sphere positive inside
        end
        @test mf.vertices == mf_pos.vertices
        @test mf.faces == mf_pos.faces
    end
end

@testset "function/array" begin

    sdf_torus = SignedDistanceField(HyperRectangle(Vec(-2,-2,-2.),Vec(4,4,4.)),0.05) do v
        (sqrt(v[1]^2+v[2]^2)-0.5)^2 + v[3]^2 - 0.25 # torus
    end

    fn_torus(v) = (sqrt(v[1]^2+v[2]^2)-0.5)^2 + v[3]^2 - 0.25 # torus

    for algo in algos
        @testset "$algo" begin
            sdf_m = HomogenousMesh(sdf_torus, algo())
            fn_m = HomogenousMesh(fn_torus, HyperRectangle(Vec(-2,-2,-2.),Vec(4,4,4.)), size(sdf_torus), algo())
            @test sdf_m.vertices == fn_m.vertices
            @test sdf_m.faces == fn_m.faces
        end
    end

end

@testset "respect origin" begin
    # verify that when we construct a mesh, that mesh:
    #   a) respects the origin of the SDF
    #   b) is correctly spaced so that symmetric functions create symmetric meshes
    f = x -> norm(x)
    bounds = HyperRectangle(Vec(-1, -1, -1), Vec(2, 2, 2))
    resolution = 0.1
    sdf = SignedDistanceField(f, bounds, resolution)

    for algorithm in (MarchingCubes(0.5),
                      MarchingTetrahedra(0.5),
                      MarchingCubes(iso=0.5,reduceverts=false),
                      MarchingTetrahedra(iso=0.5,reduceverts=false))
        @testset "$(typeof(algorithm))" begin
            mesh = @inferred GLNormalMesh(sdf, algorithm)
            # should be centered on the origin
            @test mean(vertices(mesh)) ≈ [0, 0, 0] atol=0.15*resolution
            # and should be symmetric about the origin
            @test maximum(vertices(mesh)) ≈ [0.5, 0.5, 0.5]
            @test minimum(vertices(mesh)) ≈ [-0.5, -0.5, -0.5]

            # Test functional Meshing
            fmesh = @inferred GLNormalMesh(f, bounds, size(sdf), algorithm)
            # should be centered on the origin
            @test mean(vertices(fmesh)) ≈ [0, 0, 0] atol=0.15*resolution
            # and should be symmetric about the origin
            @test maximum(vertices(fmesh)) ≈ [0.5, 0.5, 0.5]
            @test minimum(vertices(fmesh)) ≈ [-0.5, -0.5, -0.5]
        end
    end

    #TODO: SurfaceNets functional variant is bugged
    # Naive Surface Nets has no accuracy guarantee, and is a weighted sum
    # so a larger tolerance is needed for this one. In addition,
    # quad -> triangle conversion is not functioning correctly
    # see: https://github.com/JuliaGeometry/GeometryTypes.jl/issues/169
    mesh = @inferred GLNormalMesh(sdf, NaiveSurfaceNets(0.5))
    # should be centered on the origin
    @test mean(vertices(mesh)) ≈ [0, 0, 0] atol=0.15*resolution
    # and should be symmetric about the origin
    @test maximum(vertices(mesh)) ≈ [0.5, 0.5, 0.5] atol=0.2
    @test minimum(vertices(mesh)) ≈ [-0.5, -0.5, -0.5] atol=0.2
end

@testset "AbstractMeshingAlgorithm interface" begin
    f = x -> norm(x) - 0.5
    bounds = HyperRectangle(Vec(-1, -1, -1), Vec(2, 2, 2))
    resolution = 0.1
    sdf = SignedDistanceField(f, bounds, resolution)

    for algo in algos
        @testset "$algo" begin
            @test_nowarn GLNormalMesh(sdf, algo())
            @inferred GLNormalMesh(sdf, algo())
            @test_nowarn GLNormalMesh(sdf.data, algo(0.5))
            @inferred GLNormalMesh(sdf.data, algo(0.5))
        end
    end
end

@testset "mixed types" begin
    # https://github.com/JuliaGeometry/Meshing.jl/issues/9
    samples = randn(10, 10, 10)
    m1 = HomogenousMesh(SignedDistanceField(
        HyperRectangle(Vec(0,0,0), Vec(1, 1, 1)),
        samples), MarchingTetrahedra())
    m2 = HomogenousMesh(SignedDistanceField(
        HyperRectangle(Vec(0,0,0), Vec(1, 1, 1)),
        Float32.(samples)), MarchingTetrahedra())
    @test all(isapprox.(vertices(m1), vertices(m2), rtol=1e-6))
    @test all(faces(m1) == faces(m2))

    @testset "forward diff" begin
        # Demonstrate that we can even take derivatives through the meshing algorithm
        f = x -> norm(x)
        bounds = HyperRectangle(Vec(-1, -1, -1), Vec(2, 2, 2))
        resolution = 0.1
        sdf = SignedDistanceField(f, bounds, resolution)
        for algo in (MarchingCubes, MarchingTetrahedra)
            @testset "$algo" begin
                let
                    function surface_distance_from_origin(isovalue)
                        # Compute the mean distance of each vertex in the isosurface
                        # mesh from the origin. This should return a value equal to the
                        # isovalue and should have a derivative of 1.0 w.r.t. the
                        # isovalue. This function is just meant to serve as an example
                        # of some quantity you might want to differentiate in a mesh,
                        # and has the benefit for testing of having a trivial expected
                        # solution.
                        mesh = HomogenousMesh(sdf, algo(isovalue))
                        mean(norm.(vertices(mesh)))
                    end

                    isoval = 0.85
                    @test surface_distance_from_origin(isoval) ≈ isoval atol=1e-2
                    @test ForwardDiff.derivative(surface_distance_from_origin, isoval) ≈ 1 atol=1e-2
                end
            end
        end
    end

    @testset "type stability" begin
        # verify that we don't lose type stability just by mixing arguments
        # of different types
        data = randn(5, 5, 5)
        iso = 0.2
        eps = 1e-4
        for algo_ty in (MarchingTetrahedra, MarchingCubes,NaiveSurfaceNets)
            algo1 = algo_ty(iso, eps)
            algo2 = algo_ty(Float64(iso), Float16(eps))
            algo3 = algo_ty(Float32(iso), Float64(eps))
            @inferred(isosurface(data, algo1, Point{3, Float16}, Face{3,UInt32}))
            @inferred(isosurface(Float32.(data), algo2, Point{3, Float32}, Face{3,Int16}))
            @inferred(isosurface(Float64.(data), algo3, Point{3, Float64}, Face{3,UInt}))
        end
    end
    @testset "Float16" begin
        sdf_torus = SignedDistanceField(HyperRectangle(Vec{3,Float16}(-2,-2,-2.),
                                                        Vec{3,Float16}(4,4,4.)),
                                        0.1, Float16) do v
            (sqrt(v[1]^2+v[2]^2)-0.5)^2 + v[3]^2 - 0.25
        end
        @test typeof(PlainMesh{Float16,Face{4,Int}}(sdf_torus,NaiveSurfaceNets(zero(Float16)))) == PlainMesh{Float16,Face{4,Int}}
        @test typeof(PlainMesh{Float16,Face{3,Int}}(sdf_torus,MarchingTetrahedra(zero(Float16)))) == PlainMesh{Float16,Face{3,Int}}
        @test typeof(PlainMesh{Float16,Face{3,Int}}(sdf_torus,MarchingCubes(zero(Float16)))) == PlainMesh{Float16,Face{3,Int}}
    end
    @testset "Integers" begin
        A = rand(Int, 10,10,10)
        for algo in algos
            @testset "$algo" begin
                p, t = isosurface(A, algo())
                @test typeof(p) <: Vector{SVector{3,Float64}}
                @test typeof(t) <: Vector{ SVector{ Meshing.default_face_length(algo()), Int } }
            end
        end
    end
    @testset "Defaults" begin
        A = rand(10,10,10)
        @test isosurface(A) == isosurface(A, MarchingCubes())
    end
end

@testset "GeometryTypes API" begin
    @testset "vararg passing" begin
        A = rand(20,20,20)
        m = GLNormalMesh(A,MarchingTetrahedra(1.0),origin=Point(Float32(0),Float32(0),Float32(0)),widths=Point(Float32(1),Float32(1),Float32(1)))
        @test typeof(m) == GLNormalMesh
    end

end

@testset "GeometryBasics API" begin
    @testset "vararg passing" begin
        A = rand(20,20,20)
        m = GB.Mesh(A,MarchingTetrahedra(1.0),origin=Point(Float32(0),Float32(0),Float32(0)),widths=Point(Float32(1),Float32(1),Float32(1)))
        @test m isa GB.Mesh
    end
end


text:
immutable WindowsPath <: AbstractPath
    parts::Tuple{Vararg{String}}
    drive::String
    root::String
end

WindowsPath() = WindowsPath(tuple(), "", "")

WindowsPath(parts::Tuple) = WindowsPath(parts, "", "")

function WindowsPath(str::AbstractString)
    if isempty(str)
        return WindowsPath(tuple("."), "", "")
    end

    if startswith(str, "\\\\?\\")
        error("The \\\\?\\ prefix is currently not supported.")
    end

    str = replace(str, POSIX_PATH_SEPARATOR, WIN_PATH_SEPARATOR)

    if startswith(str, "\\\\")
        error("UNC paths are currently not supported.")
    elseif startswith(str, "\\")
        tokenized = split(str, WIN_PATH_SEPARATOR)

        return WindowsPath(tuple(WIN_PATH_SEPARATOR, String.(tokenized[2:end])...), "", WIN_PATH_SEPARATOR)
    elseif contains(str, ":")
        l_drive, l_path = splitdrive(str)

        tokenized = split(l_path, WIN_PATH_SEPARATOR)

        l_root = isempty(tokenized[1]) ? WIN_PATH_SEPARATOR : ""

        if isempty(tokenized[1])
            tokenized = tokenized[2:end]
        end

        if !isempty(l_drive) || !isempty(l_root)
            tokenized = tuple(string(l_drive, l_root), tokenized...)
        end

        return WindowsPath(tuple(String.(tokenized)...), l_drive, l_root)
    else
        tokenized = split(str, WIN_PATH_SEPARATOR)

        return WindowsPath(tuple(String.(tokenized)...), "", "")
    end
end

function ==(a::WindowsPath, b::WindowsPath)
    return lowercase.(parts(a)) == lowercase.(parts(b)) &&
        lowercase(drive(a)) == lowercase(drive(b)) &&
        lowercase(root(a)) == lowercase(root(b))
end
Base.String(path::WindowsPath) = joinpath(parts(path)...)
parts(path::WindowsPath) = path.parts
drive(path::WindowsPath) = path.drive
root(path::WindowsPath) = path.root

function Base.show(io::IO, path::WindowsPath)
    print(io, "p\"")
    if isabs(path)
        print(io, replace(anchor(path), "\\", "/"))
        print(io, join(parts(path)[2:end], "/"))
    else
        print(io, join(parts(path), "/"))
    end
    print(io, "\"")
end

function isabs(path::WindowsPath)
    return !isempty(drive(path)) || !isempty(root(path))
end

expanduser(path::WindowsPath) = path


text:
function create_search_tree(
    dataset::DataSet{T, I},
    progressbar::Progress;
    mincuts::I = 8,
    maxleafsize::I = 200)::SpacePartitioningTree{T, I} where {T<:AbstractFloat, I<:Integer}

    sugg_cuts = (dataset.N / maxleafsize) ^ (1 / dataset.P)
    cuts = ceil(I, max(mincuts, sugg_cuts))

    recdepth = ceil(I, log(dataset.N / maxleafsize) / log(cuts))

    #define dimension list
    dimensionlist = if cuts > mincuts
        [i for i = 1:dataset.P]
    else
        [i for i = 1:recdepth]
    end

    leafsize = ceil(I, dataset.N / cuts^recdepth)
    @debug "cuts $cuts\tleafsize $leafsize\tRec. Depth $recdepth"
    cutlist = zeros(T, 0)

    if recdepth > 0
        createleafs(dataset, progressbar, dimensionlist, cutlist, cuts, leafsize, 1)
    end

    dataset.partitioningtree = SpacePartitioningTree(cuts, leafsize, dimensionlist, length(dimensionlist), cutlist)
end

function createleafs(
    dataset::DataSet{T, I},
    progressbar::Progress,
    dimensionlist::Vector{I},
    cutlist::Vector{T},
    cuts::I,
    leafsize::I,
    StartID::I = 0) where {T<:AbstractFloat, I<:Integer}

    remainingRec::I = length(dimensionlist)

    thisLeaf::I = leafsize * cuts^(remainingRec - 1)
    bigLeaf::I = leafsize * cuts^remainingRec

    startInt::I = StartID
    stopInt::I = StartID + bigLeaf - 1
    if stopInt > dataset.N
        stopInt = dataset.N
    end

    sortID = sortperm(dataset.data[dimensionlist[1], startInt:stopInt])

    dataset.data[:, startInt:stopInt] = dataset.data[:, sortID.+startInt.-1]
    dataset.logprob[startInt:stopInt] = dataset.logprob[sortID.+startInt.-1]
    dataset.weights[startInt:stopInt] = dataset.weights[sortID.+startInt.-1]
    dataset.ids[startInt:stopInt] =     dataset.ids[sortID.+startInt.-1]
    dataset.sortids[startInt:stopInt] = dataset.sortids[sortID.+startInt.-1]

    @assert remainingRec >= 1
    start::I = 0
    stop::I = StartID - 1

    for i = 1:cuts
        if stop == dataset.N
            continue
        end

        start = stop + 1
        stop += thisLeaf
        if stop > dataset.N
            stop = dataset.N
        end

        push!(cutlist, dataset.data[dimensionlist[1], start])

        if remainingRec > 1
            createleafs(dataset, progressbar, dimensionlist[2:end], cutlist, cuts, leafsize, start)
        else
            next!(progressbar)
        end
    end
end

function search(
    dataset::DataSet{T, I},
    searchvol::HyperRectVolume{T},
    searchpoints::Bool = false)::SearchResult where {T<:AbstractFloat, I<:Integer}

    res = SearchResult(T, I)
    #searchpoints = false
    search!(res, dataset, searchvol, searchpoints)

    res
end

function search!(
    result::SearchResult{T, I},
    dataset::DataSet{T, I},
    searchvol::HyperRectVolume{T},
    searchpoints::Bool = false) where {T<:AbstractFloat, I<:Integer}

    result.points = 0
    resize!(result.pointIDs, 0)
    result.maxLogProb = -Inf
    result.minLogProb = Inf
    result.maxWeightProb = -Inf
    result.minWeightProb = Inf

    currentRecursion::I = 0
    currentDimension::I = 0
    treePos = zeros(I, 0)

    maxRecursion::I = dataset.partitioningtree.recursiondepth
    maxI::I = length(dataset.partitioningtree.cutlist)

    i::I = 1
    if maxI == 0
        #only on leaf
        for n::I = 1:dataset.N
            inVol = true
            for p::I = 1:dataset.P
                if dataset.data[p, n] < searchvol.lo[p] || dataset.data[p, n] > searchvol.hi[p]
                    inVol = false
                    break
                end
            end

            if inVol
                result.points += 1

                if searchpoints
                    push!(result.pointIDs, n)
                end
                prob = dataset.logprob[n]
                w = log(dataset.weights[n]) + prob
                result.maxLogProb = max(result.maxLogProb, prob)
                result.minLogProb = min(result.minLogProb, prob)
                result.maxWeightProb = max(result.maxWeightProb, w)
                result.minWeightProb = min(result.minWeightProb, w)
            end
        end
        return

    end
    while i <= maxI
        if currentRecursion < maxRecursion
            currentRecursion += 1
        end

        if length(treePos) < currentRecursion
            push!(treePos, 1)
        else
            treePos[currentRecursion] += 1
        end
        while treePos[currentRecursion] > dataset.partitioningtree.cuts
            deleteat!(treePos, currentRecursion)
            currentRecursion -= 1
            treePos[currentRecursion] += 1
        end
        currentDimension = dataset.partitioningtree.dimensionlist[currentRecursion]

        diff::I = 1
        for r::I = 1:(maxRecursion - currentRecursion)
            diff += dataset.partitioningtree.cuts^r
        end
        low::T = dataset.partitioningtree.cutlist[i]
        high::T = i + diff > maxI ? dataset.partitioningtree.cutlist[end] : dataset.partitioningtree.cutlist[i+diff]
        if treePos[currentDimension] == dataset.partitioningtree.cuts
            high = Inf
        end

        if searchvol.lo[currentDimension] > high || searchvol.hi[currentDimension] < low
            #skip this interval
            i += diff
            if currentDimension < maxRecursion
                currentRecursion -= 1
            end
            continue
        end

        #if on deepest recursion check for points
        if currentRecursion == maxRecursion
            startID, stopID = getDataPositions(dataset, treePos)

            searchInterval!(result, dataset, searchvol, startID, stopID, searchpoints)
        end
        i += 1
    end

end

@inline function getDataPositions(
    dataset::DataSet{T, I},
    TreePos::Vector{I}) where {T<:AbstractFloat, I<:Integer}

    maxRecursion = dataset.partitioningtree.recursiondepth
    startID = 1
    recCntr = maxRecursion
    for t in TreePos
        recCntr -= 1
        startID += dataset.partitioningtree.leafsize * dataset.partitioningtree.cuts^recCntr * (t-1)
    end
    stopID = startID + dataset.partitioningtree.leafsize - 1
    if stopID > dataset.N
        stopID = dataset.N
    end

    return startID, stopID
end

function searchInterval!(
    result::SearchResult{T, I},
    dataset::DataSet{T, I},
    searchvol::HyperRectVolume{T},
    start::I,
    stop::I,
    searchpoints::Bool) where {T<:AbstractFloat, I<:Integer}

    dimsort = dataset.partitioningtree.dimensionlist[dataset.partitioningtree.recursiondepth]

    for i = start:stop
        if dataset.data[dimsort, i] > searchvol.hi[dimsort]
            break
        end
        inVol = true
        for p = 1:dataset.P
            if dataset.data[p, i] < searchvol.lo[p] || dataset.data[p, i] > searchvol.hi[p]
                inVol = false
                break
            end
        end

        if inVol
            result.points += 1

            if searchpoints
                push!(result.pointIDs, i)
            end
            prob = dataset.logprob[i]
            w = log(dataset.weights[i]) + prob
            result.maxLogProb = max(result.maxLogProb, prob)
            result.minLogProb = min(result.minLogProb, prob)
            result.maxWeightProb = max(result.maxWeightProb, w)
            result.minWeightProb = min(result.minWeightProb, w)
        end
    end
end


text:
# "ctrl_fxns"
function count_thresh(array, thresholdcount)
    returnarray = []
    count_dict = Dict(collect(zip(array, zeros(size(array)))))
    for a in array
        count_dict[a] +=1
    end
    for (k,v) in count_dict
        if v >= thresholdcount
            push!(returnarray,k)
        end
    end
    return returnarray
end
# newlifespanfeeders = count_thresh(ar, 2)

function process_raw_output(x)
    ar = []
    for y in split(x)
        push!(ar, split(y, ",")[1])
    end
    return ar
end

function when_last_checkedin(feederarray, numberoflines2check)
    tailend = readlines(`tail -$numberoflines2check mqtt/ctimelog.csv`)
    array = []
    for x in tailend
        y = split(x, ",")
        absolutetime = parse(y[3]) + Dates.datetime2unix(DateTime(Date(y[4])))
        push!(array, [y[1] absolutetime])
    end
    array = vcat(array...)
    for feeder in feederarray
        array = array[find(array[:,1] .== feeder),:] # within time window
        if length(array) < 1
            println("Feeder $feeder not found")
            continue
        end
        println("Feeder $feeder's last checkin:")
        println(Dates.unix2datetime(maximum(array[:,2])))
    end
end
# when_last_checkedin(["148c1d00" "ef901d00" "f8921d00"],10000)

# creating lifespan commands for entering into bash of RPi3 hub

# mosquitto_pub -t deltacal -m 30aea43912b4,no,88
function debug_cal(ids)
  for id in ids
    #id_times = times+rand(range, length(times))
    run(`mosquitto_pub -t deltacal -m $id,no,88`)
    println("Setting $id for calibration next check in")
  end
end

function get_feeding_regimes()
    times7 = []
    for (i,t) in enumerate([7 9 11 13 15 17 19])
        push!(times7, t*60*60+2)
    end
    times7[end] = times7[end] - (4*60)

    times12 = []
    for (t) in 7.5:1:18.5
        push!(times12, Int(round(t*60*60+30)))
    end

    times24 = []
    for (t) in 7.5:0.46:18.5
        push!(times24, Int(round(t*60*60+30)))
    end

    morning_feed = []
    for t in 7.5:1:13
        push!(morning_feed,Int(round(t*60*60)))
    end

    evening_feed = []
    for t in 13.5:1:19
        push!(evening_feed, Int(round(t*60*60)))
    end

    one_last_feed = [Int(round(18.2*60*60))]
    return Dict("times7" => times7, "times12" => times12, "morning_feed" => morning_feed, "evening_feed" => evening_feed, "one_last_feed" => one_last_feed, "times24" => times24)
end

function tosecond(datetime)
    secondsout = Dates.Second(datetime)
    secondsout += Dates.Second(Dates.Minute(datetime))
    secondsout += Dates.Second(Dates.Hour(datetime))
    return Dates.value(secondsout)
end

function getfeedingtimesFromNow(sec_interval, numberfeedings)
    debugtimes = []
    nowsec = tosecond(now())
    for (t) in (nowsec+sec_interval):sec_interval:(nowsec+(sec_interval*numberfeedings))
        println(t)
        push!(debugtimes, t)
    end
    return debugtimes
end
# getfeedingtimesFromNow(90, 10)

function getfeedingtimesFromNow(sec_interval, numberfeedings, nowinSec)
    debugtimes = []
    nowsec = nowinSec
    for (t) in (nowsec+sec_interval):sec_interval:(nowsec+(sec_interval*numberfeedings))
        println(t)
        push!(debugtimes, t)
    end
    return debugtimes
end
# getfeedingtimesFromNow(90, 10)

function debug_commands_plusRand(ids, times)
    #array = []
  for id in ids
    #id_times = times+rand(range, length(times))
    times = [x+rand(-60:60) for x in times]
    joinedtimes = join(times, ",")
    run(`mosquitto_pub -t deltaorder -m $id:$joinedtimes`)
    #push!(array, cmdout)
    println("changed orders for $id to $joinedtimes")
  end
  #return array
end
#  x = debug_commands_plusRand(checked_in_feeders, truncated)

function debug_commands(ids, times)
  for id in ids
    #id_times = times+rand(range, length(times))
    times = [x for x in times]
    joinedtimes = join(times, ",")
    run(`mosquitto_pub -t deltaorder -m $id:$joinedtimes`)
    #push!(array, cmdout)
    println("changed orders for $id to $joinedtimes")
  end
end

function get_checked_feeders(timeInterval, repeats, filelocale::AbstractString)
    #timewindow = (tosecond(now())- timeInterval) # in seconds
    timewindow = timeInterval
    tailend = readlines(`tail -500 $filelocale`)
    array = []
    for x in tailend
        y = split(x, ",")
        push!(array, [y[1] Meta.parse(y[3])])
    end
    array = vcat(array...)
    array = array[findall(array[:,2] .> timewindow),:] # within time window
    checked_in_feeders = count_thresh(array[:,1], repeats)
    return checked_in_feeders
end
#Ft2k.get_checked_feeders(57412, 1, "/home/ubuntu/ft2k/Software/Watchman/data/ctimelog.csv")

function get_checked_feeders(timeIntervalmin, repeats)
    timewindow = (tosecond(now())- timeIntervalmin*60)
    tailend = readlines(`tail -100 mqtt/ctimelog.csv`)
    array = []
    for x in tailend
        y = split(x, ",")
        push!(array, [y[1] parse(y[3])])
    end
    array = vcat(array...)
    array = array[findall(array[:,2] .> timewindow),:] # within time window
    checked_in_feeders = count_thresh(array[:,1], repeats)
    return checked_in_feeders
end
# get_checked_feeders(100, 3)

function get_checked_feeders(starttimeInterval::DateTime, endtimeInterval::DateTime, repeats)
    timewindow = (tosecond(starttimeInterval), tosecond(endtimeInterval))# (tosecond(now())- timeIntervalmin*60)
    tailend = readlines(`tail -1000 mqtt/ctimelog.csv`)
    array = []
    for x in tailend
        y = split(x, ",")
        push!(array, [y[1] parse(y[3])])
    end
    array = vcat(array...)
    array = array[find(timewindow[1].<  array[:,2] .< timewindow[2]),:] # within time window
    checked_in_feeders = count_thresh(array[:,1], repeats)
    return checked_in_feeders
end
# starttimeInterval=  DateTime(2018,7,30,14,5)
# endtimeInterval = DateTime(2018,7,30,14,20)
# get_checked_feeders(DateTime(2018,7,30,14,10), DateTime(2018,7,30,14,20), 3)

# NOTE: this one needs using JuliaDB, SQLite, IterableTables
function get_current_regime(feederarray)
    db = SQLite.DB(joinpath(ENV["HOME"],"Software", "Server", "db.sqlite"))
    sqlite_db = table(SQLite.query(db, "SELECT * FROM clients;"))
    x = filter(p -> in(p.client_id, feederarray), sqlite_db)
    x = pushcol(x, :feedcount, map(p -> length(split((p.orders).value,",")), x))
    println(x)
end
# get_current_regime(feederarray)

# "utils.jl"
function load_config()
    diroi = joinpath(@__DIR__, "..", "config", "config.toml")
    return Pkg.TOML.parsefile(diroi)
end


text:
#!/usr/bin/julia

# Trizen
# 09 November 2016
# https://github.com/trizen

# A simple implementation of Euler's totient function.

# See also:
#   https://www.youtube.com/watch?v=fq6SXByItUI
#   https://en.wikipedia.org/wiki/Euler%27s_totient_function

using Primes

function Φ(n::Int64)

    for p in keys(factor(n))
        n -= div(n, p)
    end

    n
end

# Φ(240) = Φ(2^4 * 3^1 * 5^1)
# Φ(240) = (2^4 - 2^3) * (3^1 - 3^0) * (5^1 - 5^0)
# Φ(240) = 64

println(Φ(240))

max_stars_repo_name: trizen/julia-scripts

max_issues_repo_name: trizen/julia-scripts

max_forks_repo_name: trizen/julia-scripts


text:
__precompile__()

module MicroLogging

# ----- Core API (in Base as of julia-0.7) -----
export
    ## Logger types
    AbstractLogger, LogLevel, NullLogger, SimpleLogger,
    # Public logger API:
    #   handle_message, shouldlog, min_enabled_level, catch_exceptions
    # (not exported, as they're not generally called by users)
    #
    ## Log creation
    @debug, @info, @warn, @error, @logmsg,
    ## Logger installation and control
    with_logger, current_logger, global_logger, disable_logging

# ----- MicroLogging & stdlib Logging API -----
export
    # TODO: configure_logging needs a big rethink (see, eg, python's logger
    # config system)
    configure_logging,
    ConsoleLogger,
    InteractiveLogger

import Base.CoreLogging:
    @debug, @info, @warn, @error, @logmsg,
    AbstractLogger, LogLevel,
    BelowMinLevel, Debug, Info, Warn, Error, AboveMaxLevel,
    disable_logging,
    handle_message, shouldlog, min_enabled_level, catch_exceptions,
    SimpleLogger, LogState,
    with_logger, current_logger, global_logger, disable_logging

include("StickyMessages.jl")

include("ConsoleLogger.jl")
include("InteractiveLogger.jl") # deprecated

include("config.jl")

end


text:
using Bilevel
using ForwardDiff
using Bilevel: svd, svd_finite
using LinearAlgebra
using BenchmarkTools

A = rand(3,3)
d = svd(A)
# d.S[3] = d.S[2]
A = d.U*Diagonal(d.S)*d.V'
z0 = A[:]
# z0 = rand(4)*200. .- 100.

function f(a)
    n = Int(sqrt(length(a)))
    A = reshape(a,n,n) .* 2.
    # if eltype(A) <: ForwardDiff.Dual
        # U,s,V = svd_finite(A)
    # else
        U,s,V = svd(A)
    # end
    return vcat(U[:],s[:],V[:])
end

# sol = f(z0)

# autodiff
J_auto = ForwardDiff.jacobian(f,z0)

# # # numerical
# ϵ = sqrt(eps(1.))
# J_num = zeros(size(J_auto))
# for i = 1:length(z0)
#     δ = zeros(length(z0))
#     δ[i] = ϵ
#     J_num[:,i] = (f(z0 + δ) .- sol)/ϵ
# end
#
# display(sol)
# println("")
# println("Auto")
# display(J_auto)
# println("")
# display(maximum(abs.(J_auto)))
# println("")
# println("")
# println("Finite diff")
# display(J_num)
# println("")
# display(maximum(abs.(J_num)))
# println("")
# println("")
#
# err = maximum(abs.(J_auto .- J_num))
# display(err)
# println("")
#
# err = sort(abs.(J_auto .- J_num)[:])
# display(err)
# println("")


text:
module ArraysTests

using Test

@testset "Interfaces" begin include("InterfaceTests.jl") end

@testset "CachedArrays" begin include("CachedArraysTests.jl") end

@testset "Maps" begin include("MapsTests.jl") end

@testset "AlgebraMaps" begin include("AlgebraMapsTests.jl") end

@testset "LazyArrays" begin include("LazyArraysTests.jl") end

@testset "CompressedArrays" begin include("CompressedArraysTests.jl") end

@testset "FilteredArraysTests" begin include("FilteredArraysTests.jl") end

@testset "Tables" begin include("TablesTests.jl") end

@testset "Reindex" begin include("ReindexTests.jl") end

@testset "PosNegReindex" begin include("PosNegReindexTests.jl") end

@testset "IdentityVectors" begin include("IdentityVectorsTests.jl") end

@testset "SubVectors" begin include("SubVectorsTests.jl") end

@testset "ArrayPairs" begin include("ArrayPairsTests.jl") end

@testset "AppendedArrays" begin include("AppendedArraysTests.jl") end

@testset "AutodiffTests" begin include("AutodiffTests.jl") end

@testset "VectorWithEntryRemovedTests" begin include("VectorWithEntryRemovedTests.jl") end

@testset "VectorWithEntryInsertedTests" begin include("VectorWithEntryInsertedTests.jl") end

@testset "PrintOpTrees" begin include("PrintOpTreesTests.jl") end

end # module


text:
using Revise
using OptiMode
using LinearAlgebra
using Statistics
using ArrayInterface
using RecursiveArrayTools
using StaticArrays
using HybridArrays
using SparseArrays
using FFTW
using LinearMaps
using GeometryPrimitives
using BenchmarkTools
using ChainRules
using Zygote
using ForwardDiff
using FiniteDifferences
using FiniteDiff
using UnicodePlots
using OhMyREPL
using Crayons.Box       # for color printing
using Zygote: @ignore, dropgrad
using Setfield: @set
using StaticArrays: Dynamic
using IterativeSolvers: bicgstabl
using Rotations: RotY, MRP
using RuntimeGeneratedFunctions
using Tullio
RuntimeGeneratedFunctions.init(@__MODULE__)

gradRM(fn,in) 			= 	Zygote.gradient(fn,in)[1]
gradFM(fn,in) 			= 	ForwardDiff.gradient(fn,in)
gradFD(fn,in;n=3)		=	FiniteDifferences.grad(central_fdm(n,1),fn,in)[1]
gradFD2(fn,in;rs=1e-2)	=	FiniteDiff.finite_difference_gradient(fn,in;relstep=rs)

derivRM(fn,in) 			= 	Zygote.gradient(fn,in)[1]
derivFM(fn,in) 			= 	ForwardDiff.gradient(fn,in)
derivFD(fn,in;n=3)		=	FiniteDifferences.grad(central_fdm(n,1),fn,in)[1]
derivFD2(fn,in;rs=1e-2)	=	FiniteDiff.finite_difference_derivative(fn,in;relstep=rs)

AD_style = BOLD*BLUE_FG #NEGATIVE*BOLD*BLUE_FG      # defined in Crayons.Box
FD_style = BOLD*RED_FG
MAN_style = BOLD*GREEN_FG
AD_style_N = NEGATIVE*BOLD*BLUE_FG #NEGATIVE*BOLD*BLUE_FG      # defined in Crayons.Box
FD_style_N = NEGATIVE*BOLD*RED_FG
MAN_style_N = NEGATIVE*BOLD*GREEN_FG

Δx,Δy,Δz,Nx,Ny,Nz = 6.0, 4.0, 0.5, 128, 128, 32;
grid = Grid(Δx,Δy,Δz,Nx,Ny,Nz)
LNx = rotate(MgO_LiNbO₃,Matrix(MRP(RotY(π/2))),name=:LiNbO₃_X);
LNxN = NumMat(LNx;expr_module=@__MODULE__())
SiO₂N = NumMat(SiO₂;expr_module=@__MODULE__())
Si₃N₄N = NumMat(Si₃N₄;expr_module=@__MODULE__())
AlOxN = NumMat(αAl₂O₃;expr_module=@__MODULE__())
##
# rwg(x) = ridge_wg_partial_etch(x[1],x[2],x[3],x[4],0.5,LNxN,SiO₂N,Δx,Δy) # partially etched ridge waveguide with dispersive materials, x[3] is partial etch fraction of top layer, x[3]*x[2] is etch depth, remaining top layer thickness = x[2]*(1-x[3]).
rwg3D(x) = ridge_wg_partial_etch3D(x[1],x[2],x[3],x[4],0.5,Si₃N₄N,SiO₂N,Δx,Δy,Δz)

p = [
       0.8,                #   top ridge width         `w_top`         [μm]
       0.4,                #   ridge thickness         `t_core`        [μm]
       0.75, #0.5,                #   ridge thickness         `t_core`        [μm]
       π / 14.0,           #   ridge sidewall angle    `θ`             [radian]
               ];
ω = inv(1.55)

rwg1 = rwg3D(p)
ε⁻¹ = smooth(ω,p,:fεs,true,rwg3D,grid)
kguess = 0.25 # k_guess(ω,ε⁻¹)
k⃗ = SVector(0.0,0.0,kguess)
M̂ = HelmholtzMap(k⃗, ε⁻¹, grid)
ms = ModeSolver(kguess, ε⁻¹, grid; nev=4, maxiter=300, tol=1e-7)
k1,Hv1 = solve_k(ω,p,rwg3D,grid;nev=1,eigind=1);
##
using KrylovKit
using IterativeSolvers
using DFTK: LOBPCG
evals,evecs,convinfo = eigsolve(x->ms.M̂*x,randn(ComplexF64,size(ms.H⃗,1)),size(ms.H⃗,2),:SR;maxiter=100,tol=1e-6,krylovdim=6,verbosity=2)
resIS = lobpcg(ms.eigs_itr; true,not_zeros=false,maxiter=100,tol=1e-6)
resDF = LOBPCG(ms.M̂,randn(ComplexF64,size(ms.H⃗)),I,ms.P̂,1e-6,300; display_progress=true)
##
function mag_mn1!(mag,mn::AbstractArray{T1,NDp2},k⃗::SVector{3,T2},g⃗) where {T1<:Real,T2<:Real,NDp2}
	local ẑ = SVector(0.,0.,1.)
	local ŷ = SVector(0.,1.,0.)
	# mv = view(mn,1:3,1,eachindex(g⃗)...)
	# nv = view(mn,1:3,2,eachindex(g⃗)...)
	# mvs = reinterpret(reshape,SVector{3,T1},mv)
	# nvs = reinterpret(reshape,SVector{3,T1},nv)
	kpg = zero(k⃗)
	@fastmath @inbounds for i ∈ eachindex(g⃗)
		@inbounds kpg = k⃗ - g⃗[i]
		@inbounds mag[i] = norm(kpg)
		@inbounds mn[1:3,2,i] .=  ( ( abs2(kpg[1]) + abs2(kpg[2]) ) > 0. ) ?  normalize( cross( ẑ, kpg ) ) : SVector(-1.,0.,0.) #[-1.,0.,0.] #ŷ
		@inbounds mn[1:3,1,i] .=  normalize( cross( mn[1:3,2,i], kpg )  )
	end
	# return mag,m,n
	return mag, mn
end

function mag_mn2!(mag,mn::AbstractArray{T1,NDp2},k⃗::SVector{3,T2},g⃗) where {T1<:Real,T2<:Real,NDp2}
	local ẑ = SVector(0.,0.,1.)
	local ŷ = SVector(0.,1.,0.)
	mv = view(mn,1:3,1,eachindex(g⃗)...)
	nv = view(mn,1:3,2,eachindex(g⃗)...)
	mvs = reinterpret(reshape,SVector{3,T1},mv)
	nvs = reinterpret(reshape,SVector{3,T1},nv)
	kpg = zero(k⃗)
	@fastmath @inbounds for i ∈ eachindex(g⃗)
		@inbounds kpg = k⃗ - g⃗[i]
		@inbounds mag[i] = norm(kpg)
		@inbounds nvs[i] =  ( ( abs2(kpg[1]) + abs2(kpg[2]) ) > 0. ) ?  normalize( cross( ẑ, kpg ) ) : SVector(-1.,0.,0.) #[-1.,0.,0.] #ŷ
		@inbounds mvs[i] =  normalize( cross( nvs[i], kpg )  )
	end
	# return mag,m,n
	return mag, mn
end

##
gr = grid
g⃗s = g⃗(gr)
mag, m⃗, n⃗ = mag_m_n(k⃗,g⃗s)
mag, mn = mag_mn(k⃗,g⃗s)
mag, mn = mag_mn(k⃗,grid)

@btime mag_m_n($k⃗,$g⃗s)
@btime mag_mn($k⃗,$g⃗s)
@btime mag_mn($k⃗,$grid)

##

g⃗s = g⃗(gr)
mag, m⃗, n⃗ = mag_m_n(k⃗,g⃗s)
d0 = randn(Complex{Float64}, (3,size(gr)...))
d1 = randn(Complex{Float64}, (3,128,128))
fftax = _fftaxes(gr)
Fp1 = plan_fft!(d0,fftax,flags=FFTW.PATIENT)
Fp2 = plan_bfft!(d0,fftax,flags=FFTW.PATIENT)
Fp3 = plan_fft(d0,fftax,flags=FFTW.PATIENT)
Fp4 = plan_bfft(d0,fftax,flags=FFTW.PATIENT)

fftax1 = (2:3)
Fp11 = plan_fft!(d1,fftax1,flags=FFTW.PATIENT)
Fp21 = plan_bfft!(d1,fftax1,flags=FFTW.PATIENT)
Fp31 = plan_fft(d1,fftax1,flags=FFTW.PATIENT)
Fp41 = plan_bfft(d1,fftax1,flags=FFTW.PATIENT)

# return HelmholtzMap{3,T}(
HelmholtzMap(
                SVector{3,Float64}(k⃗),
                gr.Nx,
                gr.Ny,
                gr.Nz,
                N(gr),
                1. / N(gr),
                g⃗s,
                mag,
                m⃗,
                n⃗,
                HybridArray{Tuple{3,Dynamic(),Dynamic(),Dynamic()},Float64}(reinterpret(reshape,Float64,m⃗)),
                HybridArray{Tuple{3,Dynamic(),Dynamic(),Dynamic()},Float64}(reinterpret(reshape,Float64,n⃗)),
            HybridArray{Tuple{3,Dynamic(),Dynamic(),Dynamic()},Complex{Float64}}(d0),# (Array{Float64}(undef,(Nx,Ny,Nz,3))),
            HybridArray{Tuple{3,Dynamic(),Dynamic(),Dynamic()},Complex{Float64}}(d0),# (Array{Float64}(undef,(Nx,Ny,Nz,3))),
                plan_fft!(d0,fftax,flags=FFTW.PATIENT), # planned in-place FFT operator 𝓕!
                plan_bfft!(d0,fftax,flags=FFTW.PATIENT), # planned in-place iFFT operator 𝓕⁻¹!
                plan_fft(d0,fftax,flags=FFTW.PATIENT), # planned in-place FFT operator 𝓕!
                plan_bfft(d0,fftax,flags=FFTW.PATIENT), # planned in-place iFFT operator 𝓕⁻¹!
                ε⁻¹,
                [ 3. * inv(sum(diag(ε⁻¹[:,:,I]))) for I in eachindex(gr)], #[ 3. * inv(sum(diag(einv))) for einv in ε⁻¹],
                [ inv(mm) for mm in mag ], # inverse |k⃗+g⃗| magnitudes for precond. ops
                0.0,
        )

##
# g⃗s = g⃗(grid)
# mag, m⃗, n⃗ = mag_m_n(k⃗,g⃗s)
# d0 = randn(Complex{Float64}, (3,size(grid)...))
# fftax = _fftaxes(grid)
# [ 3. * inv(sum(diag(ε⁻¹[:,:,ix,iy,iz]))) for ix=1:grid.Nx,iy=1:grid.Ny,iz=1:grid.Nz]
# epsave = [ 3. * inv(sum(diag(ε⁻¹[:,:,I]))) for I in CartesianIndices(size(grid))]
# ε⁻¹[:,:,gridinds[64,64,64]]

##

##

##


text:
module Triangulate
using DocStringExtensions
using Triangle_jll

include("ctriangulateio.jl")
include("triangulateio.jl")
include("plot.jl")

export triangulate
export triunsuitable
export TriangulateIO
export numberofpoints
export numberofsegments
export numberoftriangles
export TriangulateError
export plot_triangulateio, plot_in_out, isplots, ispyplot

end


text:
mutable struct Point{X<:FR,Y<:FR} <: Object
    x::X
    y::Y
    opts::Option{LittleDict{Symbol,Any}}
end

point(x::FR, y::FR; kw...) = Point(x, y, dict(;kw...))
point(x::FR, y::FR, d::LittleDict{Symbol,Any}) = Point(x, y, d)

function str(p::Point, b::Board)
    xs, xss, xrp = strf(p.x, "FX", b)
    ys, yss, yrp = strf(p.y, "FY", b)
    opts = get_opts(p)
    jss = js".create('point', [$xss, $yss], $opts);"
    return xs * ys * b.name * replacefn(jss.s, xrp, yrp)
end

mutable struct Segment <: Object
    a::String
    b::String
    opts::Option{LittleDict{Symbol,Any}}
end

segment(a::String, b::String; kw...) = Segment(a, b, dict(;kw...))
segment(a::String, b::String, d::LittleDict{Symbol,Any}) = Segment(a, b, d)

function str(s::Segment, b::Board)
    opts = get_opts(s)
    sa, sb = s.a, s.b
    jss = js".create('segment', [$sa, $sb], $opts);"
    arp = "INSERT_$(sa)_$(sb)" => opts.name
    return b.name * replacefn(jss.s, arp)
end


text:
using Plots
include("../src/Types.jl")
@recipe function f(::Type{Vertex2D}, V::Vertex2D, add_marker = false)
    linecolor --> :blue
    seriestype := :scatter
    markershape --> (add_marker ? :circle : :none)
    delete!(plotattributes, :add_marker)
    vec.(V)
end
V = Vertex2D.([rand([-1, -2, 0, -3], 2) for i in 1:5])
plt = plot(; framestyle = :origin,
    aspect_ratio = 1)
# plot!(plt, V)
vec.(V)
plot!(plt, vec.(V), seriestype = :scatter)


text:
mldatasets = ["MNIST", "FashionMNIST", "CIFAR10", "SVHN2"]

"""
	load_mldatasets_data(dataset::String[; anomaly_class_ind=1])

Loads MNIST, FMNIST, SVHN2 and CIFAR10 datasets.
"""
function load_mldatasets_data(dataset::String; anomaly_class_ind::Int=1)
	(dataset in mldatasets) ? nothing : error("$dataset not available in MLDatasets.jl")
	# since the functions for MLDatasets.MNIST, MLDatasets.CIFAR10 are the same
	sublib = getfield(MLDatasets, Symbol(dataset))

	# do we need to download it?
	isdir(joinpath(first(DataDeps.standard_loadpath), dataset)) ?
		nothing : sublib.download()

	# now get the data
	tr_x, tst_x = sublib.traintensor(Float32), sublib.testtensor(Float32)
	if ndims(tr_x) == 3 # some datasets are 3D tensors :(
		tr_x = reshape(tr_x, size(tr_x,1), size(tr_x,2), 1, :)
		tst_x = reshape(tst_x, size(tst_x,1), size(tst_x,2), 1, :)
	end
	data = cat(tr_x, tst_x, dims=4)
	labels = cat(sublib.trainlabels(), sublib.testlabels(), dims=1)

	# return the normal and anomalous data
	label_list = unique(labels)
	aclass_inds = labels .== label_list[anomaly_class_ind] # indices of anomalous data
	data[:,:,:,.!aclass_inds], data[:,:,:,aclass_inds]
end

"""
	mnist_c_categories()

List of the available MNIST-C categories.
"""
mnist_c_categories() = readdir(joinpath(datadep"MNIST-C", "mnist_c"))

"""
	load_mnist_c_data(;category::String="brightness")

Loads the corrupted MNIST dataset for given category. Returns normal (non-corrupted) and anomalous data. For
a list of available corruption categories, run `GenerativeAD.Datasets.mnist_c_categories()`.
"""
function load_mnist_c_data(;category::String="brightness", kwargs...)
	dp = joinpath(datadep"MNIST-C", "mnist_c")
	available_categories = mnist_c_categories()
	!(category in available_categories) ? error("Requested category $category not found, $(available_categories) available.") : nothing
	tr_x_a, tst_x_a = map(x->Float32.(permutedims(npzread(joinpath(dp, category, x))/255, (2,3,4,1))),
		["train_images.npy", "test_images.npy"])
	tr_x_n, tst_x_n = map(x->Float32.(permutedims(npzread(joinpath(dp, "identity", x))/255, (2,3,4,1))),
		["train_images.npy", "test_images.npy"])
	return cat(tr_x_n, tst_x_n, dims=4), cat(tr_x_a, tst_x_a, dims=4)
end

"""
	mvtec_ad_categories()

List of the available MVTec-AD categories.
"""
mvtec_ad_categories() = readdir(datadep"MVTec-AD")

"""
	load_mvtec_ad_data(;category::String="bottle", img_size=nothing)

Loads the corrupted MNIST dataset for given category. Returns normal (non-corrupted) and anomalous data. For
a list of available corruption categories, run `GenerativeAD.Datasets.mnist_c_categories()`.
"""
function load_mvtec_ad_data(;category::String="bottle", img_size=nothing, kwargs...)
	# return the downscaled version if needed
	isnothing(img_size) ? nothing :
		return load_mvtec_ad_data_downscaled(category=category, img_size=img_size)
	# else get the original
	dp = datadep"MVTec-AD"
	available_categories = mvtec_ad_categories()
	!(category in available_categories) ? error("Requested category $category not found, $(available_categories) available.") : nothing
	cdp = joinpath(dp, category)
	# training and testing normal images
	tr_imgs_n = load_images_rgb(joinpath(cdp, "train/good"))
	tst_imgs_n = load_images_rgb(joinpath(cdp, "test/good"))
	# anomalous images
	tst_imgs_a = map(x->load_images_rgb(joinpath(cdp, "test", x)), filter(x->x!="good", readdir(joinpath(cdp, "test"))))

	return cat(tr_imgs_n, tst_imgs_n, dims=4), cat(tst_imgs_a..., dims=4)
end

"""
	load_mvtec_ad_data_downscaled(;img_size = 256, category::String="wood")

Returns the normal and anomalous MVTec-AD data for a given category.
"""
function load_mvtec_ad_data_downscaled(;img_size = 256, category::String="bottle", kwargs...)
	inpath = datadir("mvtec_ad/downscaled_data")
	!ispath(inpath) ? error("downscaled data not available at $inpath") : nothing
	data = load(joinpath(inpath, "$(category)_$(img_size).bson"))
	return data[:normal], data[:anomalous]
end

"""
	img_to_array_rgb(img)

Convert RGB img to a 3D tensor.
"""
img_to_array_rgb(img) = permutedims(Float32.(channelview(RGB.(img))), (2,3,1))

"""
	load_images_rgb(path)

Load all images in a path.
"""
function load_images_rgb(path)
	imgs = []
	for f in readdir(path, join=true)
		try # some files are corrupted
			img = img_to_array_rgb(load(f))
			push!(imgs, img)
		catch
			nothing
		end
	end
	cat(imgs..., dims=4)
end


text:
function displayMainWindow()
    initConfigs()
    initLoadedState()

    windowTitle = baseTitle
    if loadedState.map !== nothing
        windowTitle = "$baseTitle - $(Maple.getSideName(loadedState.side))"
    end

    # Everything else should be ready, safe to make the window
    global window = Window(windowTitle, 1280, 720, true, true, gravity=GdkGravity.GDK_GRAVITY_CENTER)

    disableLoadingScreen = get(debug.config, "DISABLE_LOADING_SCREEN", false)
    progressDialog = nothing

    if !disableLoadingScreen
        progressDialog = getProgressDialog(
            title="Starting Ahorn",
            description="Getting ready...",
            stylesheet="",
            parent=window
        )

        showall(progressDialog)
    end

    configured = configureCelesteDir()
    if !configured
        error("Celeste installation not configured")
    end

    @progress extractGamedata(storageDirectory, get(debug.config, "ALWAYS_FORCE_GAME_EXTRACTION", false)) "Extracting game assets..." progressDialog

    @progress initSprites() "Getting sprites ready..." progressDialog
    @progress initTilerMetas() "Getting tiles ready..." progressDialog
    @progress initCamera() "Getting camera ready..." progressDialog
    @progress initLangdata() "Getting language files ready..." progressDialog
    @progress initMenubar() "Getting menubar ready..." progressDialog

    @progress Backup.initBackup(persistence) "Starting backup handler..." progressDialog

    @progress updateToolDisplayNames!(loadedTools) "Getting tools ready..." progressDialog
    @progress updateToolList!(toolList) "Getting tools ready..." progressDialog

    @progress initExternalModules() "Loading external plugins..." progressDialog

    @progress updateTreeView!(roomList, getTreeData(loadedState.map)) "Populating room list..." progressDialog

    roomListVisiblityFunctions[get(persistence, "room_list_column_visibility", "all")]()

    selectLoadedRoom!(roomList)
    selectRow!(toolList, 1)

    global box = Box(:v)
    global grid = generateMainGrid()

    push!(box, grid)
    push!(window, box)

    @progress initSignals(window, canvas) "Hooking up signals..." progressDialog

    # If the window was previously maximized, maximize again
    if get(persistence, "start_maximized", false)
        maximize(window)
    end

    @progress draw(canvas) "Warming up drawing..." progressDialog
    @progress showall(window) "Waiting for main window..." progressDialog

    visible(window, true)

    if progressDialog !== nothing
        Gtk.destroy(progressDialog)
    end

    interactiveWorkaround(window)
end


text:
function _ffill(weights_df::AbstractDataFrame, date_df::AbstractDataFrame)::AbstractDataFrame
    """
    Description:
    - Implementation of LOCF from: https://hongtaoh.com/en/2021/06/27/julia-ffill/
    Input:
    - weights_df: A dataframe of the weights over the limited period
    - date_df: A dataframe over the full period
    """
    res_df = leftjoin(date_df, weights_df, on = :Date)
    sort!(res_df, :Date)
    ffill(v) = v[accumulate(max, [i*!ismissing(v[i]) for i in 1:length(v)], init=1)]
    return DataFrame([ffill(res_df[!, c]) for c in names(res_df)], names(res_df))
end

function _realised_returns(returns_df::AbstractDataFrame,
        ffill_df::AbstractDataFrame,
        date_df::AbstractDataFrame)::AbstractDataFrame
    cols = ncol(returns_df)
    matrix_prod = sum(Matrix(returns_df[:, 2:cols]) .* Matrix(ffill_df[:, 2:cols]), dims = 2)
    temp_cumsum = cumsum(matrix_prod, dims = 1)
    cumsum_df = DataFrame(Realised = vec(temp_cumsum))
    return hcat(date_df, cumsum_df)
end

"""
    walk_forward
The *Arcsine distribution* has probability density function
"""
function walk_forward(returns::AbstractDataFrame, func, rebalancing::Int)
    """
    Input:
    - Returns over full period
    - A function to calculate the weights
    - How often to rebalance the portfolio
    Output:
    """
    period, cols  = size(returns) # The length of the time series number of assets + date
    @assert rebalancing < period
    num_assets = cols - 1
    returns_assets = returns[:, 2:cols] # Filter out the dates
    index = 1 # Indexing in the while loop
    boolean = true
    # Create empty dataframe to populate with dates
    dates_df = select(returns, :Date)[1:rebalancing:end, :]
    # Preallocate the weights
    tot_weights = zeros(nrow(dates_df), num_assets)
    while boolean
        #date = returns.Dates[rebalancing * index]
        if index * rebalancing > period
            returns_index = Matrix(returns_assets[1:end, 1:num_assets])
            calc_weights = equal_weights(returns_index)
            tot_weights[index, :] = calc_weights'
            boolean = false # End the while statement
        else
            returns_index = Matrix(returns_assets[1:index * rebalancing, 1:num_assets]) # index * rebalancing - (rebalancing-1)
            calc_weights = func(returns_index)
            tot_weights[index, :] = calc_weights'
        end
        index = index + 1
    end
    weights_df = DataFrame(tot_weights, names(returns)[2:end])
    weights_df = hcat(dates_df, weights_df)
    # Fill all the dates where there is no weights!
    ffill_df = _ffill(weights_df, select(returns, :Date))
    # Calculate the total returns!
    realised_returns = _realised_returns(returns, ffill_df, select(returns, :Date))
    # Plot the results
    _plot(realised_returns)
end


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

export backward!

"""
`backward!(beta,Nt,A,c,y)`

In-place version of [`backward`](@ref).
"""
function backward!(beta::AbstractMatrix{T},Nt::Integer,
                   A::AbstractMatrix{T},c::Vector{T},y::AbstractMatrix{T}) where {T<:AbstractFloat}
  Ns = size(A,1)
  fill!(beta,zero(T))
  beta[end,:] .= 1
  for t = Nt-1:-1:1
    for j = 1:Ns
      for k = 1:Ns
        aij = A[j,k]
        if  aij != T(0)
          beta[t,j] += beta[t+1,k]*aij*y[t+1,k]
        end
      end
    end
    for j = 1:Ns
      beta[t,j] /= c[t+1]
    end
  end
  return beta
end

"""
`backward!(beta,Nt,A,c,t2tr,y)`

In-place version [`backward`](@ref) with constrined path.
"""
function backward!(beta::AbstractMatrix{T},
                   Nt::Integer,
                   A::AbstractMatrix{T},
                   c::Vector{T},
                   t2IJ::Vector{Pair{Vector{D},Vector{D}}},
                   y::AbstractMatrix{T}) where {D<:Integer,T<:AbstractFloat}
  Ns = size(A,1)
  fill!(beta,zero(T))
  @views beta[end,:] .= 1

  for t = Nt-1:-1:1
    I,J = t2IJ[t]
    for k in eachindex(I)
      i,j = I[k],J[k]
      beta[t,i] += beta[t+1,j]*A[i,j]*y[t+1,j]
    end
    for j = 1:Ns
      beta[t,j] /= c[t+1]
    end
  end
  return beta
end

export backward

"""
`backward(Nt,A,c,y)`

Computes the scaled backward probabilities (``\\bar{\\beta}``) of an Hidden Markov Model (HMM) of `Ns` states with transition matrix `A` (must be a `Ns` × `Ns` matrix), scale coefficients `c` (must be a `Nt`-long vector which is produced by [`forward`](@ref)) and observation likelihoods `y` (must be of size `Nt2` × `Ns` with `Nt2` ≥ `Nt`).
"""
function backward(Nt::Integer,A::AbstractArray{T},c::Vector{T},y::AbstractMatrix{T}) where {T <: AbstractFloat}
  Ns = size(y,2)
  beta = zeros(T,Nt,Ns)
  backward!(beta,Nt,A,c,y)
end

"""
`backward(Nt,A,c,t2tr,y)`

Computes the scaled backward probabilities with constrained path indicated by the `t2tr` vector. See [`forward`](@ref) for a description of the structure and meaning of `t2tr`.
"""
function backward(Nt::Integer,
                  A::AbstractMatrix{T},
                  c::Vector{T},
                  t2tr,
                  y::AbstractMatrix{T}) where {D <: Integer,
                                               T <: AbstractFloat}
  Ns = size(A,1)
  beta = zeros(T,Nt,Ns)
  backward!(beta,Nt,A,c,t2tr,y)
end


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

"""
    ConjArray(array)

A lazy-view wrapper of an `AbstractArray`, taking the elementwise complex conjugate. This
type is usually constructed (and unwrapped) via the [`conj`](@ref) function (or related
[`ctranspose`](@ref)), but currently this is the default behavior for `RowVector` only. For
other arrays, the `ConjArray` constructor can be used directly.

# Examples
```jldoctest
julia> [1+im, 1-im]'
1×2 RowVector{Complex{Int64},ConjArray{Complex{Int64},1,Array{Complex{Int64},1}}}:
 1-1im  1+1im

julia> ConjArray([1+im 0; 0 1-im])
2×2 ConjArray{Complex{Int64},2,Array{Complex{Int64},2}}:
 1-1im  0+0im
 0+0im  1+1im
```
"""
struct ConjArray{T,N,A<:AbstractArray} <: AbstractArray{T,N}
    parent::A
end

@inline ConjArray(a::AbstractArray{T,N}) where {T,N} = ConjArray{conj_type(T),N,typeof(a)}(a)

const ConjVector{T,V<:AbstractVector} = ConjArray{T,1,V}
@inline ConjVector(v::AbstractVector{T}) where {T} = ConjArray{conj_type(T),1,typeof(v)}(v)

const ConjMatrix{T,M<:AbstractMatrix} = ConjArray{T,2,M}
@inline ConjMatrix(m::AbstractMatrix{T}) where {T} = ConjArray{conj_type(T),2,typeof(m)}(m)

# This type can cause the element type to change under conjugation - e.g. an array of complex arrays.
@inline conj_type(x) = conj_type(typeof(x))
@inline conj_type(::Type{T}) where {T} = promote_op(conj, T)

@inline parent(c::ConjArray) = c.parent
@inline parent_type(c::ConjArray) = parent_type(typeof(c))
@inline parent_type(::Type{ConjArray{T,N,A}}) where {T,N,A} = A

@inline size(a::ConjArray) = size(a.parent)
IndexStyle(::CA) where {CA<:ConjArray} = IndexStyle(parent_type(CA))
IndexStyle(::Type{CA}) where {CA<:ConjArray} = IndexStyle(parent_type(CA))

@propagate_inbounds getindex(a::ConjArray{T,N}, i::Int) where {T,N} = conj(getindex(a.parent, i))
@propagate_inbounds getindex(a::ConjArray{T,N}, i::Vararg{Int,N}) where {T,N} = conj(getindex(a.parent, i...))
@propagate_inbounds setindex!(a::ConjArray{T,N}, v, i::Int) where {T,N} = setindex!(a.parent, conj(v), i)
@propagate_inbounds setindex!(a::ConjArray{T,N}, v, i::Vararg{Int,N}) where {T,N} = setindex!(a.parent, conj(v), i...)

@inline similar(a::ConjArray, ::Type{T}, dims::Dims{N}) where {T,N} = similar(parent(a), T, dims)

# Currently, this is default behavior for RowVector only
@inline conj(a::ConjArray) = parent(a)

# Helper functions, currently used by RowVector
@inline _conj(a::AbstractArray) = ConjArray(a)
@inline _conj(a::AbstractArray{T}) where {T<:Real} = a
@inline _conj(a::ConjArray) = parent(a)
@inline _conj(a::ConjArray{T}) where {T<:Real} = parent(a)

max_stars_repo_name: Mikewl/julia

max_issues_repo_name: Mikewl/julia

max_forks_repo_name: Mikewl/julia


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

using Markdown
using InteractiveUtils

# ╔═╡ f5cc0f04-5d99-11eb-3abe-bf3fccdac9e6
begin
	import Pkg
	Pkg.activate(Base.current_project())
	using ReinforcementLearning
	using Flux
	using Statistics
	using Plots
	using Distributions
end

# ╔═╡ 28cb401c-5d9a-11eb-2098-1939ea55d918
begin
	const pₕ = 0.4
	const WinCapital = 100

	decode_state(s::Int) = s - 1
	encode_state(s::Int) = s + 1

	function nextstep(s::Int, a::Int)
		s = decode_state(s)
		a = min(s, a)
		if s == WinCapital || s==0
			[(0., false,encode_state(s))=>1.0]
		else
			[
				((s+a >= WinCapital ? 1.0 : 0.), false, encode_state(min(s+a, WinCapital))) => pₕ,
			 	(0., false, encode_state(max(s-a, 0))) => 1-pₕ
			]
		end
	end

	struct GamblerProblemEnvModel <: AbstractEnvironmentModel
		cache
	end

	function GamblerProblemEnvModel()
		GamblerProblemEnvModel(
			Dict(
				(s,a) => nextstep(s,a)
				for s in 1:(WinCapital+1) for a in 1:WinCapital
			)
		)
	end

	RLBase.state_space(m::GamblerProblemEnvModel) = Base.OneTo(WinCapital+1)
	RLBase.action_space(m::GamblerProblemEnvModel) = Base.OneTo(WinCapital)

	(m::GamblerProblemEnvModel)(s, a) = m.cache[(s,a)]
end

# ╔═╡ 7cfd0a76-5d9c-11eb-35dc-176affb90ad3
V = TabularVApproximator(;n_state=1+WinCapital,opt=Descent(1))

# ╔═╡ bcede678-5d9c-11eb-05d0-fb25818336cf
RLZoo.value_iteration!(V=V, model=GamblerProblemEnvModel(), γ=1.0, max_iter=1000)

# ╔═╡ 005c156e-5d9e-11eb-11be-1b69be12b13b
plot(V.table[2:end-1])

# ╔═╡ Cell order:
# ╠═f5cc0f04-5d99-11eb-3abe-bf3fccdac9e6
# ╠═28cb401c-5d9a-11eb-2098-1939ea55d918
# ╠═7cfd0a76-5d9c-11eb-35dc-176affb90ad3
# ╠═bcede678-5d9c-11eb-05d0-fb25818336cf
# ╠═005c156e-5d9e-11eb-11be-1b69be12b13b