text:
module PkgLicenses

export LICENSES, readlicense, license

const LICENSES = Dict(
    "MIT" => "MIT \"Expat\" License",
    "BSD" => "Simplified \"2-clause\" BSD License",
    "ISC" => "Internet Systems Consortium Licence",
    "ASL" => "Apache License, Version 2.0",
    "MPL" => "Mozilla Public License, Version 2.0",
    "GPL-2.0+" => "GNU Public License, Version 2.0+",
    "GPL-3.0+" => "GNU Public License, Version 3.0+",
    "LGPL-2.1+" => "Lesser GNU Public License, Version 2.1+",
    "LGPL-3.0+" => "Lesser GNU Public License, Version 3.0+",
    "CC0" => "Creative Commons Universal Public Domain Dedication (CC0 1.0)"
)

"Read license text from specified file and location"
function readlicense(lic::AbstractString,
                     dir::AbstractString=normpath(dirname(@__FILE__), "..", "res", "licenses"))
    return read(joinpath(dir, lic), String)
end

"""
    license([io,], [name])
Shows available licenses if no parameters specified. If a license name is specified as a parameter then the function shows the license text.
"""
function license(io::IO, lic::AbstractString="")
    if isempty(lic)
        lics =""
        for k in LICENSES |> keys |> collect |> sort
            lics *= " $k | $(LICENSES[k])\n"
        end
        maxdesclen = maximum(map(length,values(LICENSES)))
        lics = "Label| Full name\n-----|$(repeat("-", maxdesclen+2))\n$lics"
        println(io, lics)
        # Base.Markdown.term(STDOUT, lics)
    else
        try
            println(io, readlicense(lic))
        catch
            printstyled(io, "License $lic is not available.", color=:red)
        end
    end
end
license(lic::AbstractString="") = license(stdout, lic)

end # module


text:
"""
    get_initial_guess(parameters)

Returns the values of parameter vector. Useful for initial guess for optimization.
"""
function get_initial_guess(parameters)
    initial_guess = zeros(size(parameters, 1))
    for i ∈ 1:size(parameters, 1)
        initial_guess[i] = parameters[i].value
    end
    return initial_guess
end
"""
    get_lower_bounds(parameters::Array{Float64)

Returns the lower bounds of parameter vector. Useful for constrained optimization.
"""
function get_lower_bounds(parameters)
    lower_bound = zeros(size(parameters, 1))
    for i ∈ 1:size(parameters, 1)
        lower_bound[i] = parameters[i].lb
    end
    return lower_bound
end
"""
    get_upper_bounds(parameters::Array{Float64)

Returns the upper bounds of parameter vector. Useful for constrained optimization.
"""
function get_upper_bounds(parameters)
    upper_bound = zeros(size(parameters, 1))
    for i ∈ 1:size(parameters, 1)
        upper_bound[i] = parameters[i].ub
    end
    return upper_bound
end
"""
    get_bounds(name::ParameterGroup, val::Symbol)

Returns the lower and upper bound of a single parameter withing a `ParameterGroup`. Useful for constrained optimization.
"""
function get_bounds(name::ParameterGroup, val::Symbol)
    if typeof(getfield(name, val)) <: BoundedParameter
        x1 = getfield(getfield(name, val), :lb)
        x2 = getfield(getfield(name, val), :ub)
        return x1, x2
    else
        @warn "No bounds defined"
        return nothing
    end
end

ext: jl


text:
###
### LongSequence specific specializations of src/biosequence/counting.jl
###
### This file is a part of BioJulia.
### License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

# Counting GC positions
let
    @info "Compiling bit-parallel GC counter for LongSequence{<:NucleicAcidAlphabet}"
    counter = :(n += gc_bitcount(chunk, BitsPerSymbol(seq)))
    compile_bitpar(
        :count_gc_bitpar,
        arguments   = (:(seq::LongSequence{<:NucleicAcidAlphabet}),),
        init_code   = :(n = 0),
        head_code   = counter,
        body_code   = counter,
        tail_code   = counter,
        return_code = :(return n)
    ) |> eval
end

Base.count(::typeof(isGC), seq::LongSequence{<:NucleicAcidAlphabet}) = count_gc_bitpar(seq)

# Counting mismatches
let
    @info "Compiling bit-parallel mismatch counter for LongSequence{<:NucleicAcidAlphabet}"

    counter = :(count += mismatch_bitcount(x, y, A()))

    compile_2seq_bitpar(
        :count_mismatches_bitpar,
        arguments = (:(seqa::LongSequence{A}), :(seqb::LongSequence{A})),
        parameters = (:(A<:NucleicAcidAlphabet),),
        init_code = :(count = 0),
        head_code = counter,
        body_code = counter,
        tail_code = counter,
        return_code = :(return count)
    ) |> eval
end
Base.count(::typeof(!=), seqa::LongSequence{A}, seqb::LongSequence{A}) where {A<:NucleicAcidAlphabet} = count_mismatches_bitpar(seqa, seqb)
Base.count(::typeof(!=), seqa::LongNucleotideSequence, seqb::LongNucleotideSequence) = count(!=, promote(seqa, seqb)...)

# Counting matches
let
    @info "Compiling bit-parallel match counter for LongSequence{<:NucleicAcidAlphabet}"

    counter = :(count += match_bitcount(x, y, A()))

    count_empty = quote
        count += match_bitcount(x, y, A())
        nempty = div(64, bits_per_symbol(A())) - div(offs, bits_per_symbol(A()))
        count -= nempty
    end

    compile_2seq_bitpar(
        :count_matches_bitpar,
        arguments = (:(seqa::LongSequence{A}), :(seqb::LongSequence{A})),
        parameters = (:(A<:NucleicAcidAlphabet),),
        init_code = :(count = 0),
        head_code = count_empty,
        body_code = counter,
        tail_code = count_empty,
        return_code = :(return count)
    ) |> eval
end
Base.count(::typeof(==), seqa::LongSequence{A}, seqb::LongSequence{A}) where {A<:NucleicAcidAlphabet} = count_matches_bitpar(seqa, seqb)
Base.count(::typeof(==), seqa::LongNucleotideSequence, seqb::LongNucleotideSequence) = count(==, promote(seqa, seqb)...)

# Counting ambiguous sites
let
    @info "Compiling bit-parallel ambiguity counter for LongSequence{<:NucleicAcidAlphabet}"

    counter = :(count += ambiguous_bitcount(x, y, A()))

    compile_2seq_bitpar(
        :count_ambiguous_bitpar,
        arguments = (:(seqa::LongSequence{A}), :(seqb::LongSequence{A})),
        parameters = (:(A<:NucleicAcidAlphabet),),
        init_code = :(count = 0),
        head_code = counter,
        body_code = counter,
        tail_code = counter,
        return_code = :(return count)
    ) |> eval
end
Base.count(::typeof(isambiguous), seqa::LongSequence{A}, seqb::LongSequence{A}) where {A<:NucleicAcidAlphabet{4}} = count_ambiguous_bitpar(seqa, seqb)
Base.count(::typeof(isambiguous), seqa::LongSequence{<:NucleicAcidAlphabet{4}}, seqb::LongSequence{<:NucleicAcidAlphabet{2}}) = count(isambiguous, promote(seqa, seqb)...)
Base.count(::typeof(isambiguous), seqa::LongSequence{<:NucleicAcidAlphabet{2}}, seqb::LongSequence{<:NucleicAcidAlphabet{4}}) = count(isambiguous, promote(seqa, seqb)...)

# Counting certain sites
let
    @info "Compiling bit-parallel certainty counter for LongSequence{<:NucleicAcidAlphabet}"

    counter = :(count += certain_bitcount(x, y, A()))

    compile_2seq_bitpar(
        :count_certain_bitpar,
        arguments = (:(seqa::LongSequence{A}), :(seqb::LongSequence{A})),
        parameters = (:(A<:NucleicAcidAlphabet),),
        init_code = :(count = 0),
        head_code = counter,
        body_code = counter,
        tail_code = counter,
        return_code = :(return count)
    ) |> eval
end
Base.count(::typeof(iscertain), seqa::LongSequence{A}, seqb::LongSequence{A}) where {A<:NucleicAcidAlphabet{4}} = count_certain_bitpar(seqa, seqb)
Base.count(::typeof(iscertain), seqa::LongSequence{<:NucleicAcidAlphabet{4}}, seqb::LongSequence{<:NucleicAcidAlphabet{2}}) = count(iscertain, promote(seqa, seqb)...)
Base.count(::typeof(iscertain), seqa::LongSequence{<:NucleicAcidAlphabet{2}}, seqb::LongSequence{<:NucleicAcidAlphabet{4}}) = count(iscertain, promote(seqa, seqb)...)

# Counting gap sites
let
    @info "Compiling bit-parallel gap counter for LongSequence{<:NucleicAcidAlphabet}"

    counter = :(count += gap_bitcount(x, y, A()))

    count_empty = quote
        count += gap_bitcount(x, y, A())
        nempty = div(64, bits_per_symbol(A())) - div(offs, bits_per_symbol(A()))
        count -= nempty
    end

    compile_2seq_bitpar(
        :count_gap_bitpar,
        arguments = (:(seqa::LongSequence{A}), :(seqb::LongSequence{A})),
        parameters = (:(A<:NucleicAcidAlphabet),),
        init_code = :(count = 0),
        head_code = count_empty,
        body_code = counter,
        tail_code = count_empty,
        return_code = :(return count)
    ) |> eval
end
Base.count(::typeof(isgap), seqa::LongSequence{A}, seqb::LongSequence{A}) where {A<:NucleicAcidAlphabet{4}} = count_gap_bitpar(seqa, seqb)
Base.count(::typeof(isgap), seqa::LongSequence{<:NucleicAcidAlphabet{4}}, seqb::LongSequence{<:NucleicAcidAlphabet{2}}) = count(isgap, promote(seqa, seqb)...)
Base.count(::typeof(isgap), seqa::LongSequence{<:NucleicAcidAlphabet{2}}, seqb::LongSequence{<:NucleicAcidAlphabet{4}}) = count(isgap, promote(seqa, seqb)...)


text:
"""
For greedy algorithm, it will the element based on the optimal value in the array at each loop in the function

This does 0-1 (each item can be chosen only once) knapsack :
pack capacity = capacity
weight of each item = weights
value of each item = values

Each loop the function will find the highest value in the array and check if the capacity is enough to store it, if enough then
the value will be added into the total_max_value.
After that the highest current value will be deleted.

```
julia> zero_one_pack!(20,[1,3,11],[2,5,30])
37
```
"""
function ga_zero_one_pack!(capacity::N, weights::V, values::V
) where {N <: Number,V <: AbstractVector}
    total_max_value = 0
    while capacity > 0 && (length(weights)!=0 || length(values)!=0)
        max_value, index = findmax(values) # get the highest value and index in the vector
        if capacity >= weights[index]
            total_max_value += max_value
            capacity -= weights[index]
        end
        splice!(values, index)
        splice!(weights, index)
    end
    return total_max_value
end

"""
This does complete/infinite (each item can be chosen infinite times) knapsack :
pack capacity = capacity
weight of each item = weights
value of each item = values

Each loop the function will find the highest value in the array and check if the capacity is enough to store it, if enough then
the value will be added into the total_max_value until the capacity cannot hold the weight of the highest current value.
After that the highest current value will be deleted.

```
julia> complete_pack!(20,[1,2,9],[1,3,20])
43
```
"""
function ga_complete_pack!(capacity::N, weights::V, values::V
) where {N <: Number,V <: AbstractVector}
    total_max_value = 0
    while capacity > 0 && (length(weights)!=0 || length(values)!=0)
        max_value, index = findmax(values)
        while capacity >= weights[index]
            total_max_value += max_value
            capacity -= weights[index]
        end
        splice!(values, index)
        splice!(weights, index)
    end
    return total_max_value
end


text:
using Documenter, Powerlog

push!(LOAD_PATH,"../src/")

makedocs()

deploydocs(
    deps = Deps.pip("mkdocs", "python-markdown-math"),
    repo = "github.com/mwsohn/Powerlog.jl.git",
    julia = "0.6",
    osname = "linux"
)


text:
using KitePodSimulator, KiteUtils
using Test

include("plot_step_response.jl")

@testset "KitePodSimulator.jl" begin
    kcu = KCU(se())
    @test get_depower(kcu) ≈ 0.236
    @test get_steering(kcu) ≈ 0.0
    set_depower_steering(kcu, 1.0, 1.0)
    on_timer(kcu)
    @test get_depower(kcu) ≈ 0.23975
    @test get_steering(kcu) ≈ 0.01
end


text:
const Mean = let
    init() = (0, 0.0)
    op((n, sum), val) = n + 1, sum + val
    value((n, sum)) = sum / n
    (; init, op, value)
end

struct ExpectationAnalysis end

(e::ExpectationAnalysis)(input::ProcessedLayer) = groupreduce(Mean, input)

"""
    expectation()

Compute the expected value of the last argument conditioned on the preceding ones.
"""
expectation() = transformation(ExpectationAnalysis())


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -3.97969642451919 <= q <= 6.896074133539049)

                     add_NL_constraint(m, :(softplus(0.518279586272822 + -0.9811103721569414*$(x[1]) + -0.6250390763797746*$(x[2]) + -0.8382325906616188*$(x[3]) + 0.6951086592760602*$(x[4])) + softplus(0.9399092682371073 + -0.6199086749993752*$(x[1]) + -0.7631209251662368*$(x[2]) + -0.009236754359899013*$(x[3]) + -0.9061282260292134*$(x[4])) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
"""
    @neutral(SET, NEUT)

Create functions to make a lazy set operation commutative with a given neutral
element set type.

### Input

- `SET`  -- lazy set operation type
- `NEUT` -- set type for neutral element

### Output

Nothing.

### Notes

This macro generates four functions (possibly two more if `@absorbing` has been
used in advance) (possibly two or four more if `@declare_array_version` has been
used in advance).

### Examples

`@neutral(MinkowskiSum, N)` creates at least the following functions:
* `neutral(::MinkowskiSum) = N`
* `MinkowskiSum(X, N) = X`
* `MinkowskiSum(N, X) = X`
* `MinkowskiSum(N, N) = N`
"""
macro neutral(SET, NEUT)
    @eval begin
        # create function to obtain the neutral element
        function neutral(::Type{$SET})
            return $NEUT
        end

        # create functions to declare the neutral element
        function $SET(X::LazySet{N}, ::$NEUT{N}) where {N<:Real}
            return X
        end
        function $SET(::$NEUT{N}, X::LazySet{N}) where {N<:Real}
            return X
        end
        function $SET(Y::$NEUT{N}, ::$NEUT{N}) where {N<:Real}
            return Y
        end

        # if the absorbing element has already been defined, create combinations
        if isdefined(@__MODULE__, :absorbing) && hasmethod(absorbing, (Type{$SET},))
            @eval(@neutral_absorbing($(esc(SET)),
                                     $(esc(NEUT)),
                                     absorbing($(esc(SET)))))
        end

        # if the array set type has already been defined, create combinations
        if isdefined(@__MODULE__, :array_constructor) &&
                hasmethod(array_constructor, (Type{$SET},))
            @eval(@array_neutral($(esc(SET)),
                                 $(esc(NEUT)),
                                 array_constructor($(esc(SET)))))
        elseif isdefined(@__MODULE__, :is_array_constructor) &&
                hasmethod(is_array_constructor, (Type{$SET},))
            @eval(@array_neutral($(esc(SET)),
                                 $(esc(NEUT)),
                                 $(esc(SET))))
        end
    end
    return nothing
end

"""
    @absorbing(SET, ABS)

Create functions to make a lazy set operation commutative with a given absorbing
element set type.

### Input

- `SET` -- lazy set operation type
- `ABS` -- set type for absorbing element

### Output

Nothing.

### Notes

This macro generates four functions (possibly two more if `@neutral` has been
used in advance) (possibly two or four more if `@declare_array_version` has been
used in advance).

### Examples

`@absorbing(MinkowskiSum, A)` creates at least the following functions:
* `absorbing(::MinkowskiSum) = A`
* `MinkowskiSum(X, A) = A`
* `MinkowskiSum(A, X) = A`
* `MinkowskiSum(A, A) = A`
"""
macro absorbing(SET, ABS)
    @eval begin
        # create function to obtain the absorbing element
        function absorbing(::Type{$SET})
            return $ABS
        end

        # create functions to declare the absorbing element
        function $SET(::LazySet{N}, Y::$ABS{N}) where {N<:Real}
            return Y
        end
        function $SET(Y::$ABS{N}, ::LazySet{N}) where {N<:Real}
            return Y
        end
        function $SET(Y::$ABS{N}, ::$ABS{N}) where {N<:Real}
            return Y
        end

        # if the neutral element has already been defined, create combinations
        if isdefined(@__MODULE__, :neutral) && hasmethod(neutral, (Type{$SET},))
            @eval(@neutral_absorbing($(esc(SET)),
                                     neutral($(esc(SET))),
                                     $(esc(ABS))))
        end

        # if the array set type has already been defined, create combinations
        if isdefined(@__MODULE__, :array_constructor) &&
                hasmethod(array_constructor, (Type{$SET},))
            @eval(@array_absorbing($(esc(SET)),
                                   $(esc(ABS)),
                                   array_constructor($(esc(SET)))))
        elseif isdefined(@__MODULE__, :is_array_constructor) &&
                hasmethod(is_array_constructor, (Type{$SET},))
            @eval(@array_absorbing($(esc(SET)),
                                   $(esc(ABS)),
                                   $(esc(SET))))
        end
    end
    return nothing
end

"""
    @declare_array_version(SET, SETARR)

Create functions to connect a lazy set operation with its array set type.

### Input

- `SET`    -- lazy set operation type
- `SETARR` -- array set type

### Output

Nothing.

### Notes

This macro generates six functions (and possibly up to eight more if
`@neutral`/`@absorbing` has been used in advance for the base and/or array set
type).

### Examples

`@declare_array_version(MinkowskiSum, MinkowskiSumArray)` creates at least the
following functions:
* `array_constructor(::MinkowskiSum) = MinkowskiSumArray`
* `is_array_constructor(::MinkowskiSumArray) = true`
* `MinkowskiSum!(X, Y)`
* `MinkowskiSum!(X, arr)`
* `MinkowskiSum!(arr, X)`
* `MinkowskiSum!(arr1, arr2)`
"""
macro declare_array_version(SET, SETARR)
    _SET! = Symbol(string(SET), '!')
    @eval begin
        # create function to obtain the array version
        function array_constructor(::Type{$SET})
            return $SETARR
        end

        # create function to check that this is an array version
        function is_array_constructor(::Type{$SETARR})
            return true
        end

        # create in-place modification functions for array version
        function $_SET!(X::LazySet{N}, Y::LazySet{N}) where {N<:Real}
            # no array type: just use the lazy operation
            return $SET(X, Y)
        end
        function $_SET!(X::LazySet{N}, arr::$SETARR{N}) where {N<:Real}
            push!(array(arr), X)
            return arr
        end
        function $_SET!(arr::$SETARR{N}, X::LazySet{N}) where {N<:Real}
            push!(array(arr), X)
            return arr
        end
        function $_SET!(arr1::$SETARR{N}, arr2::$SETARR{N}) where {N<:Real}
            append!(array(arr1), array(arr2))
            return arr1
        end
        # handle method ambiguities with neutral elements
        if isdefined(@__MODULE__, :neutral) && hasmethod(neutral, (Type{$SET},))
            @eval(@array_neutral($(esc(SET)),
                                 neutral($(esc(SET))),
                                 $(esc(SETARR))))
        end
        if isdefined(@__MODULE__, :neutral) && hasmethod(neutral, (Type{$SETARR},))
            @eval(@array_neutral($(esc(SETARR)),
                                 neutral($(esc(SETARR))),
                                 $(esc(SETARR))))
        end
        # handle method ambiguities with absorbing elements
        if isdefined(@__MODULE__, :absorbing) && hasmethod(absorbing, (Type{$SET},))
            @eval(@array_absorbing($(esc(SET)),
                                   absorbing($(esc(SET))),
                                   $(esc(SETARR))))
        end
        if isdefined(@__MODULE__, :absorbing) && hasmethod(absorbing, (Type{$SETARR},))
            @eval(@array_absorbing($(esc(SETARR)),
                                   absorbing($(esc(SETARR))),
                                   $(esc(SETARR))))
        end
    end
    return nothing
end

"""
    @neutral_absorbing(SET, NEUT, ABS)

Create two functions to avoid method ambiguties for a lazy set operation with
respect to neutral and absorbing element set types.

### Input

- `SET`  -- lazy set operation type
- `NEUT` -- set type for neutral element
- `ABS`  -- set type for absorbing element

### Output

A quoted expression containing the function definitions.

### Examples

`@neutral_absorbing(MinkowskiSum, N, A)` creates the following functions as
quoted expressions:
* `MinkowskiSum(N, A) = A`
* `MinkowskiSum(A, N) = A`
"""
macro neutral_absorbing(SET, NEUT, ABS)
    return quote
        function $SET(::$NEUT{N}, Y::$ABS{N}) where {N<:Real}
            return Y
        end
        function $SET(Y::$ABS{N}, ::$NEUT{N}) where {N<:Real}
            return Y
        end
    end
end

"""
    @array_neutral(FUN, NEUT, SETARR)

Create two functions to avoid method ambiguities for a lazy set operation with
respect to the neutral element set type and the array set type.

### Input

- `FUN`     -- function name
- `NEUT`    -- set type for neutral element
- `SETARR`  -- array set type

### Output

A quoted expression containing the function definitions.

### Examples

`@array_neutral(MinkowskiSum, N, ARR)` creates the following functions as
quoted expressions:
* `MinkowskiSum(N, ARR) = ARR`
* `MinkowskiSum(ARR, N) = ARR`
"""
macro array_neutral(FUN, NEUT, SETARR)
    return quote
        function $FUN(::$NEUT{N}, X::$SETARR{N}) where {N<:Real}
            return X
        end
        function $FUN(X::$SETARR{N}, ::$NEUT{N}) where {N<:Real}
            return X
        end
    end
end

"""
    @array_absorbing(FUN, ABS, SETARR)

Create two functions to avoid method ambiguities for a lazy set operation with
respect to the absorbing element set type and the array set type.

### Input

- `FUN`     -- function name
- `ABS`     -- set type for absorbing element
- `SETARR`  -- array set type

### Output

A quoted expression containing the function definitions.

### Examples

`@array_absorbing(MinkowskiSum, ABS, ARR)` creates the following functions as
quoted expressions:
* `MinkowskiSum(ABS, ARR) = ABS`
* `MinkowskiSum(ARR, ABS) = ABS`
"""
macro array_absorbing(FUN, ABS, SETARR)
    return quote
        function $FUN(Y::$ABS{N}, ::$SETARR{N}) where {N<:Real}
            return Y
        end
        function $FUN(::$SETARR{N}, Y::$ABS{N}) where {N<:Real}
            return Y
        end
    end
end


text:
# This is adapted from Base.isapprox
# https://github.com/JuliaLang/julia/blob/381693d3dfc9b7072707f6d544f82f6637fc5e7c/base/floatfuncs.jl#L222-L291
function isapprox_le(x::Number, y::Number; atol::Real=0, rtol::Real=Base.rtoldefault(x,y,atol))
    x <= y || (isfinite(x) && isfinite(y) && abs(x-y) <= max(atol, rtol*max(abs(x), abs(y))))
end


text:
JSONRPC.parse_params(::Type{Val{Symbol("workspace/didChangeWatchedFiles")}}, params) = DidChangeWatchedFilesParams(params)
function process(r::JSONRPC.Request{Val{Symbol("workspace/didChangeWatchedFiles")},DidChangeWatchedFilesParams}, server)
    for change in r.params.changes
        uri = change.uri

        startswith(uri, "file:") || continue

        if change.type == FileChangeTypes["Created"] || change.type == FileChangeTypes["Changed"]
            if hasdocument(server, URI2(uri))
                doc = getdocument(server, URI2(uri))

                # Currently managed by the client, we don't do anything
                if get_open_in_editor(doc)
                    continue
                else
                    filepath = uri2filepath(uri)
                    content = try
                        s = read(filepath, String)
                        if !isvalid(s)
                            deletedocument!(server, URI2(uri))
                            continue
                        end
                        s
                    catch err
                        isa(err, Base.IOError) || isa(err, Base.SystemError) || rethrow()
                        deletedocument!(server, URI2(uri))
                        continue
                    end

                    set_text!(doc, content)
                    set_is_workspace_file(doc, true)
                    parse_all(doc, server)
                end
            else
                filepath = uri2filepath(uri)
                content = try
                    s = read(filepath, String)
                    isvalid(s) || continue
                    s
                catch err
                    isa(err, Base.IOError) || isa(err, Base.SystemError) || rethrow()
                    continue
                end

                doc = Document(uri, content, true, server)
                setdocument!(server, URI2(uri), doc)
                parse_all(doc, server)
            end
        elseif change.type == FileChangeTypes["Deleted"]
            if hasdocument(server, URI2(uri))
                doc = getdocument(server, URI2(uri))

                # We only handle if currently not managed by client
                if !get_open_in_editor(doc)
                    deletedocument!(server, URI2(uri))

                    publishDiagnosticsParams = PublishDiagnosticsParams(uri, missing, Diagnostic[])
                    JSONRPCEndpoints.send_notification(server.jr_endpoint, "textDocument/publishDiagnostics", publishDiagnosticsParams)
                else
                    # TODO replace with accessor function once the other PR
                    # that introduces the accessor is merged
                    doc._workspace_file = false
                end
            end
        else
            error("Unknown change type.")
        end
    end
end

JSONRPC.parse_params(::Type{Val{Symbol("workspace/didChangeConfiguration")}}, params) = params
function process(r::JSONRPC.Request{Val{Symbol("workspace/didChangeConfiguration")},Dict{String,Any}}, server::LanguageServerInstance)
    request_julia_config(server)
end

function request_julia_config(server)
    response = JSONRPCEndpoints.send_request(server.jr_endpoint, "workspace/configuration", ConfigurationParams([
        (ConfigurationItem(missing, "julia.format.$opt") for opt in fieldnames(DocumentFormat.FormatOptions))...;
        ConfigurationItem(missing, "julia.lint.run");
        (ConfigurationItem(missing, "julia.lint.$opt") for opt in fieldnames(StaticLint.LintOptions))...
        ]))

    # TODO Make sure update_julia_config can deal with the response
    if length(response) == length(fieldnames(DocumentFormat.FormatOptions)) + 1 + length(fieldnames(StaticLint.LintOptions))
        server.format_options = DocumentFormat.FormatOptions(
            response[1]===nothing ? 0 : response[1],
            response[2]===nothing ? false : response[2],
            response[3]===nothing ? false : response[3],
            response[4]===nothing ? false : response[4],
            response[5]===nothing ? false : response[5],
            response[6]===nothing ? false : response[6],
            response[7]===nothing ? false : response[7],
            response[8]===nothing ? false : response[8],
            response[9]===nothing ? false : response[9],
            response[10]===nothing ? false : response[10],
            response[11]===nothing ? false : response[11])

        N = length(fieldnames(DocumentFormat.FormatOptions)) + 1
        x = response[N]
        new_lint_opts = StaticLint.LintOptions(
            response[N + 1]===nothing ? false : response[N + 1],
            response[N + 2]===nothing ? false : response[N + 2],
            response[N + 3]===nothing ? false : response[N + 3],
            response[N + 4]===nothing ? false : response[N + 4],
            response[N + 5]===nothing ? false : response[N + 5],
            response[N + 6]===nothing ? false : response[N + 6],
            response[N + 7]===nothing ? false : response[N + 7],
            response[N + 8]===nothing ? false : response[N + 8],
            response[N + 9]===nothing ? false : response[N + 9],
        )

        new_run_lint_value = x===nothing ? false : true
        if new_run_lint_value != server.runlinter || any(getfield(new_lint_opts, n) != getfield(server.lint_options, n) for n in fieldnames(StaticLint.LintOptions))
            server.lint_options = new_lint_opts
            server.runlinter = new_run_lint_value
            for doc in getdocuments_value(server)
                StaticLint.check_all(getcst(doc), server.lint_options, server)
                empty!(doc.diagnostics)
                mark_errors(doc, doc.diagnostics)
                publish_diagnostics(doc, server)
            end
        end
    end
end

JSONRPC.parse_params(::Type{Val{Symbol("workspace/didChangeWorkspaceFolders")}}, params) = DidChangeWorkspaceFoldersParams(params)
function process(r::JSONRPC.Request{Val{Symbol("workspace/didChangeWorkspaceFolders")}}, server)
    for wksp in r.params.event.added
        push!(server.workspaceFolders, uri2filepath(wksp.uri))
        load_folder(wksp, server)
    end
    for wksp in r.params.event.removed
        delete!(server.workspaceFolders, uri2filepath(wksp.uri))
        remove_workspace_files(wksp, server)
    end
end

JSONRPC.parse_params(::Type{Val{Symbol("workspace/symbol")}}, params) = WorkspaceSymbolParams(params)
function process(r::JSONRPC.Request{Val{Symbol("workspace/symbol")},WorkspaceSymbolParams}, server)
    syms = SymbolInformation[]
    for doc in getdocuments_value(server)
        bs = collect_toplevel_bindings_w_loc(getcst(doc), query = r.params.query)
        for x in bs
            p, b = x[1], x[2]
            push!(syms, SymbolInformation(valof(b.name), 1, false, Location(doc._uri, Range(doc, p)), nothing))
        end
    end

    return syms
end


text:
using DecisionMakingEnvironments
import DecisionMakingEnvironments: create_simple_discrete_bandit, create_simple_discrete_contextualbandit
using Test
using Distributions

@testset "Simple Bandit Tests" begin
    b = create_simple_discrete_bandit(3; noise=2.0)
    @test length(b.A) == 3
    @test eltype(b.A) <: Int
    @test b.r(1) isa Float64
    @test meta_information(b)[:minreward] isa Float64

    b = create_simple_discrete_contextualbandit(3; noise_x=1.0, noise_r=0.5)
    @test size(b.X) == (1,2)
    @test all(b.X[1, :] .== [-2.0, 2.0])
    @test length(b.A) == 3
    @test eltype(b.A) <: Int
    x = b.d()
    @test x isa Float64
    @test x ≥ b.X[1, 1]
    @test x ≤ b.X[1, 2]
    @test b.r(x, 1) isa Float64
    @test meta_information(b)[:maxreward] isa Float64
end

@testset "Bandit Interface Tests" begin
    prob = create_simple_discrete_bandit(3; noise=0.0)
    r = sample(prob, 2)
    @test r == 2

    π = (x=0)->Categorical(ones(length(prob.A)) ./ length(prob.A))
    J = sample_objective(prob, π)
    @test J isa Float64
    @test J ≥ meta_information(prob)[:minreward]
    @test J ≤ meta_information(prob)[:maxreward]

    prob = create_simple_discrete_contextualbandit(3; noise_x=1.0, noise_r=0.0)
    x = sample(prob)
    @test x isa Float64
    @test x ≥ prob.X[1, 1]
    @test x ≤ prob.X[1, 2]
    x, r = sample(prob, x, 1)
    @test x isa Float64
    @test x ≥ prob.X[1, 1]
    @test x ≤ prob.X[1, 2]
    @test r ≥ meta_information(prob)[:minreward]
    @test r ≤ meta_information(prob)[:maxreward]
    J = sample_objective(prob, π)
    @test J isa Float64
    @test J ≥ meta_information(prob)[:minreward]
    @test J ≤ meta_information(prob)[:maxreward]
end


text:
module Utils

export spikes2chain

"""
    spikes2chain(ts, inds, net, dt; start=0, stop=nothing)

Convert the spike sequence of a network into the underlying Markov chain.

# Arguments
- `ts::Vector{Float64}` - vector of spike times
- `inds::Vector{Int}` - vector of indices of neurons that emitted the spikes
- `net` - the network that generated the spikes
- `dt::Float64` - time step of network simulation / sampling process
- `start` - time point at which to start sampling (0 by default)
- `stop` - time point at which to stop sampling (if nothing is passed,
           end point is determined by last spike)
"""
function spikes2chain end
# here: only inferface definition;
# actual implememtation differs per network type and is provided in
# respective modules

"The sigmoid / logistic function: σ(x) = 1 / (1 - exp(-x))"
function σ(x)
    if x >= 0
        return 1 / (1 - exp(-x))
    else
        return exp(x) / (1 + exp(x)) # numerically stable version for negative numbers
    end
end

end # module Utils


text:
module SailingPolars

using DelimitedFiles

export SailingPolar, SurfaceParameters, getpolardata, deg2rad, rad2deg, cartesian2polar
export polar2cartesian, haversine, inverse_haversine, boatspeed, bestvectorspeed
export sailingspeed, sailsegmenttime

"""
    Structure to represent a polar CSV file's data.
Contains a matrix, speeds, of sailing speeds indexed by wind velocity and angle of boat to wind
winds is a list of wind speeds
degrees is a list of angles in degrees of direction relative to the wind
Note 0.0 degrees is directly into the wind, 180 degrees is directly downwind.
"""
struct SailingPolar
    winds::Vector{Float32}
    degrees::Vector{Float32}
    speeds::Matrix{Float32} # speeds[wind direction degrees, windspeed knots]
end

"""
    struct SurfaceParameters
Structure that represents wind and surface current direction and velocity for a given position
Angles in degrees, velocities in knots
"""
struct SurfaceParameters
    winddeg::Float32
    windkts::Float32
    currentdeg::Float32
    currentkts::Float32
end

"""
function getpolardata(filename)
Read a sailing polar CSV file and return a SailingPolar containing the file data.
A sailing polar file is a CSV file, with ';' used as the comma separator instead of a comma.
The first line of file contains labels for the wind velocities that make up columns, and
the first entry of each row makes up a column of angle of sailing direction from wind in degrees
"""
function getpolardata(filename)
    datacells, headercells = readdlm(filename, ';', header=true)
    winds = map(x -> parse(Float32, x), headercells[2:end])
    degrees = datacells[:, 1]
    speeds = datacells[:, 2:end]
    return SailingPolar(winds, degrees, speeds)
end

const R = 6372800  # Earth's approximate radius in meters

"""
    deg2rad(deg)
Convert degrees to radians
"""
deg2rad(deg) = (deg * π / 180.0 + 2π) % 2π

"""
    rad2deg(rad)
Convert radians to degrees
"""
rad2deg(rad) = (rad * (180.0 / π) + 360.0) % 360.0

"""
    cartesian2polard(x, y)
Convert x, y coordinates to polar coordinates with angle in degrees
"""
cartesian2polard(x, y) = sqrt(x * x + y * y), atand(x, y)

"""
    polard2cartesian(r, deg)
Convert polar coordinates in degrees to cartesian x, y coordinates
"""
polard2cartesian(r, deg) = r .* sincosd(deg)

"""
    function haversine(lat1, lon1, lat2, lon2)
Calculate the haversine function for two points on the Earth's surface.
Given two latitude, longitude pairs in degrees for a point on the Earth,
get distance in meters and the initial direction of travel in degrees for
movement from point 1 to point 2.
"""
function haversine(lat1, lon1, lat2, lon2)
    dlat = lat2 - lat1
    dlon = lon2 - lon1
    a = sind(dlat / 2)^2 + cosd(lat1) * cosd(lat2) * sind(dlon / 2)^2
    c = 2.0 * asind(sqrt(a))
    theta = atand(sind(dlon) * cosd(lat2),
        cosd(lat1) * sind(lat2) - sind(lat1) * cosd(lat2) * cosd(dlon))
    theta = (theta + 360) % 360
    return R * c * 0.5399565, theta
end

"""
    function inverse_haversine(lat1, lon1, distance, direction)
Calculate an inverse haversine function.
Takes the point of origin in degrees (latitude, longitude), distance in meters, and
initial direction in degrees, and returns the latitude and longitude of the endpoint
in degrees after traveling the specified distance.
"""
function inverse_haversine(lat1, lon1, distance, direction)
    lat2 = asind(sind(lat1) * cos(distance / R) + cosd(lat1) * sin(distance / R) * cosd(direction))
    lon2 = lon1 + atand(sind(direction) * sin(distance / R) * cosd(lat1),
                       cos(distance / R) - sind(lat1) * sind(lat2))
    return lat2, lon2
end

"""
    function boatspeed(sp::SailingPolar, pointofsail, windspeed)
Calculate the expected sailing speed in a specified direction in knots,
given sailing polar data, a desired point of sail in degrees, and wind speed in knots
"""
function boatspeed(sp::SailingPolar, pointofsail, windspeed)
    winds, degrees, speeds = sp.winds, sp.degrees, sp.speeds
    udeg = findlast(t -> t <= pointofsail, degrees)
    odeg = findfirst(t -> t >= pointofsail, degrees)
    uvel = findlast(t -> t <= windspeed, winds)
    ovel = findfirst(t -> t >= windspeed, winds)
    if any(t -> t == nothing, [udeg, odeg, uvel, ovel])
        return -1.0
    end
    frac = (odeg == udeg && uvel == ovel) ? 1.0 :
            (odeg == udeg) ? (windspeed - winds[uvel]) / (winds[ovel] - winds[uvel]) :
            (uvel == ovel) ? (pointofsail - degrees[udeg]) / (degrees[odeg] - degrees[udeg]) :
            ((pointofsail - degrees[udeg]) / (degrees[odeg] - degrees[udeg]) +
            (windspeed - winds[uvel]) / (winds[ovel] - winds[uvel])) / 2
    return speeds[udeg, uvel] + frac * (speeds[odeg, ovel] - speeds[udeg, uvel])
end

"""
    sailingspeed(sp::SailingPolar, azimuth, dir, ws)
Calculate the expected net boat speed in a desired direction (azimuth).
This is generally different from the actual boat speed in its actual direction.
Directions are in degrees (dir is wind direction), velocity of wind (ws) is in knots.
"""
sailingspeed(sp, azimuth, dir, ws) = boatspeed(sp, dir, ws) * cosd(abs(dir - azimuth))

"""
    function bestvectorspeed(sp::SailingPolar, dirtravel, dirwind, windspeed, dircur, velcur)
Calculate the net direction and velocity of a sailing ship.
Arguments are sailing polar data, direction of travel in degrees from north, wind direction in
degrees from north, wind velocity in knots, surface current direction in degrees, and
current velocity in knots.
"""
function bestvectorspeed(sp::SailingPolar, dirtravel, dirwind, windspeed, dircur, velcur)
    pointofsail = (dirtravel - dirwind) % 360.0
    pointofsail = pointofsail < 0 ? pointofsail + 360.0 : pointofsail
    pointofsail = pointofsail > 180.0 ? 360.0 - pointofsail : pointofsail
    VMG = boatspeed(sp, pointofsail, windspeed)
    other, idx = findmax([sailingspeed(sp, pointofsail, x, windspeed) for x in sp.degrees])
    if other > VMG
        pointofsail = sp.degrees[idx]
        VMG = other
    end
    dirchosen = deg2rad(dirwind + pointofsail)
    wx, wy = VMG * sin(dirchosen), VMG * cos(dirchosen)
    curx, cury = velcur * sin(deg2rad(dircur)), velcur * cos(deg2rad(dircur))
    return rad2deg(atan(wy + cury, wx + curx)), sqrt((wx + curx)^2 + (wy + cury)^2)
end

"""
    function sailsegmenttime(sp::SailingPolar, p::SurfaceParameters, lat1, lon1, lat2, lon2)
Calculate the trip time in minutes from (lat1, lon1) to the destination (lat2, lon2).
Uses the data in SurfaceParameters for wind and current velocity and direction.
"""
function sailsegmenttime(sp::SailingPolar, p::SurfaceParameters, lat1, lon1, lat2, lon2)
    distance, dir = haversine(lat1, lon1, lat2, lon2)
    dir2, vel = bestvectorspeed(sp, dir, p.winddeg, p.windkts, p.currentdeg, p.currentkts)
    # endlat2, endlon2 = inverse_haversine(lat1, lon1, distance, dir2)
    # minutes/s * m / (knots * (m/s / knot)) = minutes
    return (1 / 60) * distance / (vel * 1.94384)
end

end # module


text:
using ALifeBenchmark
using Printf
using Serialization
using Dates
import Random
using Profile

id = string(time_ns())

ANCESTOR = LARGE_ANCESTOR
const NUM_ITERATIONS = 1e10

create_model(logger) = TierraModel(ANCESTOR, logger)
function print_function(trial, step, model)
    if step % 1e7 == 0
        @info "$id \t $step \t $(step / NUM_ITERATIONS)"
        @printf "%i %i (%s): Memory in use: %2i %% \t Num of organisms: %i \t Num of free blocks: %i\n" trial step Dates.format(now(), "HH:MM:SS")  round(Int64, 100*model.used_memory / model.memory_size) length(model.organism_keys) length(model.free_blocks)
    end
end

const NUM_TRIALS = 1

collect_distribution(create_model, execute_slice!, print_function, NUM_ITERATIONS, NUM_TRIALS, SLICE_SIZE)

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
using Documenter, Kaleido

makedocs(;
    modules=[Kaleido],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        hide("internals.md"),
    ],
    repo="https://github.com/tkf/Kaleido.jl/blob/{commit}{path}#L{line}",
    sitename="Kaleido.jl",
    authors="Takafumi Arakaki <aka.tkf@gmail.com>",
    strict=true,
)

deploydocs(;
    repo="github.com/tkf/Kaleido.jl",
)


text:
module actuator_examples
using FinEtools
using PyCall

function actuator2()
    # MEMS actuator.   Thermal analysis.
    x0 =  0.0*phun("micro*m");
    x1 = x0+5.0*phun("micro*m");
    x2 = x1+10.0*phun("micro*m");
    x3 = x2+10.0*phun("micro*m");
    x4 = x3+10.0*phun("micro*m");
    y0 = 0.0*phun("micro*m");
    y4 = 250.0*phun("micro*m");
    y3 = y4-10.0*phun("micro*m");
    y2 = y3-10.0*phun("micro*m");
    y1 = y2-10.0*phun("micro*m");
    t = 2.0*phun("micro*m");
    h = 0.1*phun("micro*m");
    z0 = 0.0*phun("micro*m");
    z3 = 2*t+h;
    z2 = z3-t;
    z1 = z2-h;
    m1 = 2*2;
    m2 = 2*2;
    m3 = 2*2;
    m4 = 3*2;
    n1 = 20*2;
    n2 = 4*2;
    n3 = 2*2;
    n4 = 2*2;
    n5 = 7*2;
    p1 = 1*2;
    p2 = 1*2;
    p3 = 1*2;
    kappa = 157*eye(3, 3)*phun("W/m/K"); # W/m/K, conductivity matrix
    DV = 5*phun("V"); # voltage drop in volt
    l  = 2*(y1+y2)/2+2*(x1+x2)/2; # length of the conductor
    resistivity  =  1.1e-5*phun("Ohm*m"); # Ohm m
    Q = DV^2/resistivity/l^2; # rate of Joule heating, W/m^3
    T_substrate = 293; # substrate temperature in degrees Kelvin

    fens,fes =  H8hexahedron([x1 y0 z0; x2 y1 z1],m2,n1,p1);
    fens1,fes1  =  H8hexahedron([x1 y1 z0;x2 y2 z1],m2,n2,p1);
    fens,fes1,fes2  =  mergemeshes(fens1, fes1, fens, fes, eps(h));
    fes =  cat(fes1,fes2);
    fens1,fes1  =  H8hexahedron([x0 y1 z0;x1 y2 z1],m1,n2,p1);
    fens,fes1,fes2  =  mergemeshes(fens1, fes1, fens, fes, eps(h));
    fes =  cat(fes1,fes2);
    fens1,fes1  =  H8hexahedron([x0 y1 z1;x1 y2 z2], m1,n2,p2);
    fens,fes1,fes2  =  mergemeshes(fens1, fes1, fens, fes, eps(h));
    fes =  cat(fes1,fes2);
    fens1,fes1  =  H8hexahedron([x0 y1 z2;x1 y2 z3],m1,n2,p3);
    fens,fes1,fes2  =  mergemeshes(fens1, fes1, fens, fes, eps(h));
    fes =  cat(fes1,fes2);
    fens1,fes1  =  H8hexahedron([x0 y2 z2;x1 y3 z3],m1,n3,p3);
    fens,fes1,fes2  =  mergemeshes(fens1, fes1, fens, fes, eps(h));
    fes =  cat(fes1,fes2);
    fens1,fes1  =  H8hexahedron([x0 y3 z2;x1 y4 z3], m1,n4,p3);
    fens,fes1,fes2  =  mergemeshes(fens1, fes1, fens, fes, eps(h));
    fes =  cat(fes1,fes2);
    fens1,fes1  =  H8hexahedron([x1 y3 z2;x3 y4 z3],m4,n4,p3);
    fens,fes1,fes2  =  mergemeshes(fens1, fes1, fens, fes, eps(h));
    fes =  cat(fes1,fes2);
    fens1,fes1  =  H8hexahedron([x3 y3 z2;x4 y4 z3],m3,n4,p3);
    fens,fes1,fes2  =  mergemeshes(fens1, fes1, fens, fes, eps(h));
    fes =  cat(fes1,fes2);
    fens1,fes1  =  H8hexahedron([x3 y0 z2;x4 y3 z3], m3,n5,p3);
    fens,fes1,fes2  =  mergemeshes(fens1, fes1, fens, fes, eps(h));
    fes =  cat(fes1,fes2);

    fens,fes  =  H8toH20(fens,fes);

    hotmater = MatHeatDiff(kappa)
    coldmater = MatHeatDiff(kappa)
    cl =  selectelem(fens, fes, box=[x0,x2,y0,y2,z0,z1],inflate = t/100);

    hotfemm  =  FEMMHeatDiff(IntegData(subset(fes,cl), GaussRule(3, 3)), hotmater)
    coldfemm  = FEMMHeatDiff(IntegData(subset(fes,setdiff(collect(1:count(fes)), cl)),  GaussRule(3, 3)), coldmater)
    geom = NodalField(fens.xyz)
    Temp = NodalField(zeros(size(fens.xyz,1),1))
    fenids = selectnode(fens, box=[x0,x4,y0,y0,z0,z3], inflate=t/1000) ; # fixed temperature on substrate
    setebc!(Temp, fenids, true, 1, T_substrate)
    applyebc!(Temp)
    numberdofs!(Temp)

    K = conductivity(hotfemm, geom, Temp) +
    conductivity(coldfemm, geom, Temp)
    fi = ForceIntensity(FFlt[Q]);
    F = distribloads(hotfemm, geom, Temp, fi, 3) +
    nzebcloadsconductivity(hotfemm, geom, Temp) +
    nzebcloadsconductivity(coldfemm, geom, Temp)

    U = K\F
    scattersysvec!(Temp,U[:])

    @pyimport matplotlib.pyplot as plt
    plt.style[:use]("seaborn-whitegrid")
    fig = plt.figure()
    ax = plt.axes()
    nList = selectnode(fens, box=[x1,x1,y0,y1,z1,z1], inflate=t/1000)
    y_i = geom.values[nList, 2]
    T_i = Temp.values[nList, 1]
    ix = sortperm(y_i)
    ax[:plot](y_i[ix], T_i[ix], color=:red, label= "hot leg")
    println("maximum(T_i) = $(maximum(T_i))")

    nList = selectnode(fens, box=[x3,x3,y0,y3,z2,z2], inflate=t/1000)
    y_o = geom.values[nList, 2]
    T_o = Temp.values[nList, 1]
    ix = sortperm(y_o)
    ax[:plot](y_o[ix], T_o[ix], color=:blue, label= "cold leg")
    plt.legend()
    ax[:set_xlabel]("Distance")
    ax[:set_ylabel]("Temperature")
    plt.show()
    # plt.xlim(0.0, 1.0)
    # plt.axis("equal")
end # actuator2

function allrun()
    println("#####################################################")
    println("# actuator2 ")
    actuator2()
    return true
end # function allrun

end # module actuator_examples


text:
using ChaosTools
using Test
using DynamicalSystemsBase #to get Systems
using DelayEmbeddings #to get columns() and embed()
@testset "DyCA Tests" begin
    println("Running dyca() tests...")

    @testset "Eigenvalues of Roessler system" begin
        # Initial conditions from `dyca()` paper
        eigenthresold = collect(0.10:0.01:0.995)
        ds = Systems.roessler(rand(3); a = 0.15, b = 0.2, c = 10);
        ts = trajectory(ds, 1000.0, dt = 0.05);
        for thresold in eigenthresold
            eigenvalues,_ = dyca(Matrix(ts),thresold)
            @test sum([round(i) for i in eigenvalues]) == 2.0
        end
    end

    @testset "Eigenvalues of Embeded Roessler system" begin
        # Initial conditions from `dyca()` paper
        eigenthresold = collect(0.90:0.01:0.995)
        ds = Systems.roessler(rand(3); a = 0.15, b = 0.2, c = 10);
        ts = trajectory(ds, 1000.0, dt = 0.05);
        cols = columns(ts);
        for thresold in eigenthresold
            for i in 1:3
                Embedded_system = Matrix(embed(cols[i], 25, 1));
                eigenvals,_ = dyca(Embedded_system,thresold)
                # eigenvalues satisfy `dyca` condtiion and have imaginary part = 0.0
                wanted_eigenvalues = eigenvals[vec(0.999 .< abs.(eigenvals) .< 1.0) .& vec(imag(eigenvals) .== 0)]
                # due to floating point errors more than 2 eigenvalues might satisy the above condition
                # might be better to use `length(psum) in [1,2,3]`
                @test length(wanted_eigenvalues) in 0:4
            end
        end
    end
end

#Test using EEG seizure data, based on Datseris's python notebook available at https://gist.github.com/Datseris/8121e5019071fe9e3cb9b1a0811ac26a
# using PyPlot, PyCall
# pyimport("mpl_toolkits.mplot3d").Axes3D
# using DataFrames, CSV
# eeg_data = DataFrame(CSV.File("input/seizure.csv")) #data available in https://github.com/JuliaDynamics/ChaosTools.jl/files/3065748/dyca.zip

# figure()
# for i=2:size(eeg_data,2)
# 	plot( eeg_data[:,i], "-")
# end
# savefig("eeg_data.png")

# eigenvalues, proj_mat, proj_eeg = dyca(Matrix(eeg_data[:, 2:end]), 0.9; norm_eigenvectors=true)
# #the eigendecomposition function used in the python version normalizes the eigenvectors. To compare the figures more easily, I also normalized them here.

# figure()
# plot(proj_eeg)
# savefig("projected_eeg_data_time_series.png")

# fig = figure(figsize=(15,14))
# ax = fig.add_subplot(2, 2, 1, projection="3d")
# ax.plot(proj_eeg[:,1], proj_eeg[:,2], proj_eeg[:,3])
# ax.view_init(elev=51, azim=86)
# ax.set_xlabel(L"x_1"); ax.set_ylabel(L"x_2"); ax.set_zlabel(L"x_3");
# ax = fig.add_subplot(2, 2, 2, projection="3d")
# ax.plot(proj_eeg[:,1], proj_eeg[:,2], proj_eeg[:,4])
# ax.view_init(elev=37, azim=173)
# ax.set_xlabel(L"x_1"); ax.set_ylabel(L"x_2"); ax.set_zlabel(L"x_3");
# ax = fig.add_subplot(2, 2, 3, projection="3d")
# ax.plot(proj_eeg[:,2], proj_eeg[:,3], proj_eeg[:,4])
# ax.view_init(elev=17, azim=36)
# ax.set_xlabel(L"x_1"); ax.set_ylabel(L"x_2"); ax.set_zlabel(L"x_3");
# ax = fig.add_subplot(2, 2, 4, projection="3d")
# ax.plot(proj_eeg[:,1], proj_eeg[:,3], proj_eeg[:,4])
# ax.view_init(elev=43, azim=148)
# ax.set_xlabel(L"x_1"); ax.set_ylabel(L"x_2"); ax.set_zlabel(L"x_3");
# savefig("projected_eeg_data_phasespace.png")


text:
# Run examples in BoltzmannMachines.jl with plotting
import BoltzmannMachines
include(joinpath(dirname(pathof(BoltzmannMachines)), "..", "test", "examples.jl"))

using BoltzmannMachinesPlots
function test_scatterhidden()
   hidden = rand(100, 2)
   labels = rand(["1", "2", "3"], 100)
   scatterhidden(hidden,
         opacity = 0.5, labels = labels)
end
test_scatterhidden()


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

module SphericalHarmonics

using StaticArrays, LinearAlgebra

import ACE

import JuLIP.MLIPs: IPBasis, alloc_B, alloc_dB, fltype
import JuLIP: alloc_temp, alloc_temp_d, evaluate!, evaluate_d!,
			     write_dict, read_dict

const JVec = SVector{3}

export SHBasis, RSHBasis

# --------------------------------------------------------
#     Coordinates
# --------------------------------------------------------

"""
`struct PseudoSpherical` : a simple datatype storing spherical coordinates
of a point (x,y,z) in the format (r, cosφ, sinφ, cosθ, sinθ).

Use `spher2cart` and `cart2spher` to convert between cartesian and spherical
coordinates.
"""
struct PseudoSpherical{T}
	r::T
	cosφ::T
	sinφ::T
	cosθ::T
	sinθ::T
end

spher2cart(S::PseudoSpherical) = S.r * JVec(S.cosφ*S.sinθ, S.sinφ*S.sinθ, S.cosθ)

function cart2spher(R::AbstractVector)
	@assert length(R) == 3
	r = norm(R)
	φ = atan(R[2], R[1])
	sinφ, cosφ = sincos(φ)
	cosθ = R[3] / r
	sinθ = sqrt(R[1]^2+R[2]^2) / r
	return PseudoSpherical(r, cosφ, sinφ, cosθ, sinθ)
end

PseudoSpherical(φ, θ) = PseudoSpherical(1.0, cos(φ), sin(φ), cos(θ), sin(θ))

"""
convert a gradient with respect to spherical coordinates to a gradient
with respect to cartesian coordinates
"""
function dspher_to_dcart(S, f_φ_div_sinθ, f_θ)
	r = S.r + eps()
   return SVector( - (S.sinφ * f_φ_div_sinθ) + (S.cosφ * S.cosθ * f_θ),
			            (S.cosφ * f_φ_div_sinθ) + (S.sinφ * S.cosθ * f_θ),
			 			                                 - (   S.sinθ * f_θ) ) / r
end

# --------------------------------------------------------
#     Indexing
# --------------------------------------------------------

"""
	sizeP(maxDegree)

Return the size of the set of Associated Legendre Polynomials ``P_l^m(x)`` of
degree less than or equal to the given maximum degree
"""
sizeP(maxDegree) = div((maxDegree + 1) * (maxDegree + 2), 2)

"""
	sizeY(maxDegree)

Return the size of the set of real spherical harmonics ``Y_{l,m}(θ,φ)`` of
degree less than or equal to the given maximum degree
"""
sizeY(maxDegree) = (maxDegree + 1) * (maxDegree + 1)

"""
	index_p(l,m)

Return the index into a flat array of Associated Legendre Polynomials ``P_l^m``
for the given indices ``(l,m)``.
``P_l^m`` are stored in l-major order i.e. [P(0,0), [P(1,0), P(1,1), P(2,0), ...]
"""
# Base.@pure index_p(l,m) = m + div(l*(l+1), 2) + 1
@inline index_p(l::Int,m::Int) = m + div(l*(l+1), 2) + 1
index_p(l::Integer, m::Integer) = index_y(Int(l), Int(m))

"""
	index_y(l,m)

Return the index into a flat array of real spherical harmonics ``Y_{l,m}``
for the given indices ``(l,m)``.
``Y_{l,m}`` are stored in l-major order i.e.
[Y(0,0), [Y(1,-1), Y(1,0), Y(1,1), Y(2,-2), ...]
"""
Base.@pure index_y(l::Int, m::Int) = m + l + (l*l) + 1
index_y(l::Integer, m::Integer) = index_y(Int(l), Int(m))

# # i = m + l + l^2 + 1
# # m = i - l - l^2 - 1
# _li2m(l, i) = i - l - l^2 - 1

# --------------------------------------------------------
#     Associated Legendre Polynomials
#     TODO: rewrite within general interface?
#           - alloc_B, alloc_dB, alloc_temp, ...
# --------------------------------------------------------

"""
`ALPCoefficients` : an auxiliary datastructure for
evaluating the associated lagrange functions
used for the spherical harmonics
"""
struct ALPCoefficients{T}
	A::Vector{T}
	B::Vector{T}
end

ALPCoefficients(maxDegree::Integer, T::Type=Float64) =
	ALPCoefficients( Vector{T}(undef, sizeP(maxDegree)),
						  Vector{T}(undef, sizeP(maxDegree)) )

"""
	compute_coefficients(L)

Precompute coefficients ``a_l^m`` and ``b_l^m`` for all l <= L, m <= l
"""
function compute_coefficients(L::Integer, T::Type=Float64)
	coeff = ALPCoefficients(L, T)
	for l in 2:L
		ls = l*l
		lm1s = (l-1) * (l-1)
		for m in 0:(l-2)
			ms = m * m
			coeff.A[index_p(l, m)] = sqrt((4 * ls - 1.0) / (ls - ms))
			coeff.B[index_p(l, m)] = -sqrt((lm1s - ms) / (4 * lm1s - 1.0))
		end
	end
	return coeff
end

"""
	compute_coefficients(L)

Create an array large enough to store an entire set of Associated Legendre
Polynomials ``P_l^m(x)`` of maximum degree L.
"""
allocate_p(L::Int) = Array{Float64}(undef, sizeP(L))

"""
	compute_p(L, x, coeff, P)

Compute an entire set of Associated Legendre Polynomials ``P_l^m(x)``
using the given coefficients, and store in the array P.
"""
function compute_p!(L::Integer, S::PseudoSpherical{T}, coeff::ALPCoefficients{T},
					     P::Array{T,1}) where {T}
	@assert length(coeff.A) >= sizeP(L)
	@assert length(coeff.B) >= sizeP(L)
	@assert length(P) >= sizeP(L)

	temp = sqrt(0.5/π)
	P[index_p(0, 0)] = temp
	if L == 0; return P; end

	P[index_p(1, 0)] = S.cosθ * sqrt(3) * temp
	temp = - sqrt(1.5) * S.sinθ * temp
	P[index_p(1, 1)] = temp

	for l in 2:L
		il = ((l*(l+1)) ÷ 2) + 1
		ilm1 = il - l
		ilm2 = ilm1 - l + 1
		for m in 0:(l-2)
			# FOR DEBUGGING TURN ON THESE ASSERTS!
			# @assert il+m == index_p(l, m)
			# @assert ilm1+m == index_p(l-1,m)
			# @assert ilm2+m == index_p(l-2,m)
			@inbounds P[il+m] = coeff.A[il+m] * (     S.cosθ * P[ilm1+m]
					                           + coeff.B[il+m] * P[ilm2+m] )
		end
		@inbounds P[il+l-1] = S.cosθ * sqrt(2 * (l - 1) + 3) * temp
		temp = -sqrt(1.0 + 0.5 / l) * S.sinθ * temp
		@inbounds P[il+l] = temp
	end

	return P
end

function compute_dp!(L::Integer, S::PseudoSpherical{T}, coeff::ALPCoefficients{T},
					     P::Array{T,1}, dP::Array{T,1}) where T
	@assert length(coeff.A) >= sizeP(L)
	@assert length(coeff.B) >= sizeP(L)
	@assert length(P) >= sizeP(L)

	temp = sqrt(0.5/π)
	P[index_p(0, 0)] = temp
	temp_d = 0.0
	dP[index_p(0, 0)] = temp_d
	if L == 0; return nothing; end

	P[index_p(1, 0)] = S.cosθ * sqrt(3) * temp
	dP[index_p(1, 0)] = -S.sinθ * sqrt(3) * temp + S.cosθ * sqrt(3) * temp_d
	temp1, temp_d = ( - sqrt(1.5) * temp,
					      - sqrt(1.5) * (S.cosθ * temp + S.sinθ * temp_d) )
	P[index_p(1, 1)] = temp1
	dP[index_p(1, 1)] = temp_d

	for l in 2:L
		m = 0
		@inbounds P[index_p(l, m)] =
				coeff.A[index_p(l, m)] * (     S.cosθ * P[index_p(l - 1, m)]
				             + coeff.B[index_p(l, m)] * P[index_p(l - 2, m)] )
		@inbounds dP[index_p(l, m)] =
			coeff.A[index_p(l, m)] * (
							- S.sinθ * P[index_p(l - 1, m)]
							+ S.cosθ * dP[index_p(l - 1, m)]
			             + coeff.B[index_p(l, m)] * dP[index_p(l - 2, m)] )

		for m in 1:(l-2)
			@inbounds P[index_p(l, m)] =
					coeff.A[index_p(l, m)] * (     S.cosθ * P[index_p(l - 1, m)]
					             + coeff.B[index_p(l, m)] * P[index_p(l - 2, m)] )
			@inbounds dP[index_p(l, m)] =
				coeff.A[index_p(l, m)] * (
								- S.sinθ^2 * P[index_p(l - 1, m)]
								+ S.cosθ * dP[index_p(l - 1, m)]
				             + coeff.B[index_p(l, m)] * dP[index_p(l - 2, m)] )
		end
		@inbounds P[index_p(l, l - 1)] = sqrt(2 * (l - 1) + 3) * S.cosθ * temp1
		@inbounds dP[index_p(l, l - 1)] = sqrt(2 * (l - 1) + 3) * (
									        -S.sinθ^2 * temp1 + S.cosθ * temp_d )

      (temp1, temp_d) = (
					-sqrt(1.0+0.5/l) * S.sinθ * temp1,
		         -sqrt(1.0+0.5/l) * (S.cosθ * temp1 * S.sinθ + S.sinθ * temp_d) )
		@inbounds P[index_p(l, l)] = temp1
		@inbounds dP[index_p(l, l)] = temp_d
	end
	return nothing
end

"""
	compute_p(L, x)

Compute an entire set of Associated Legendre Polynomials ``P_l^m(x)`` where
``0 ≤ l ≤ L`` and ``0 ≤ m ≤ l``. Assumes ``|x| ≤ 1``.
"""
function compute_p(L::Integer, S::PseudoSpherical{T}) where {T}
	P = Array{T}(undef, sizeP(L))
	coeff = compute_coefficients(L)
	compute_p!(L, S, coeff, P)
	return P
end

function compute_dp(L::Integer, S::PseudoSpherical{T}) where {T}
	P = Array{T}(undef, sizeP(L))
	dP = Array{T}(undef, sizeP(L))
	coeff = compute_coefficients(L)
	compute_dp!(L, S, coeff, P, dP)
	return P, dP
end

compute_p(L::Integer, θ::Real) =
	compute_p(L, PseudoSpherical(0.0, 0.0, 0.0, cos(θ), sin(θ)))

compute_dp(L::Integer, θ::Real) =
	compute_dp(L, PseudoSpherical(0.0, 0.0, 0.0, cos(θ), sin(θ)))

# ------------------------------------------------------------------------
#                  Spherical Harmonics
# ------------------------------------------------------------------------

"""
evaluate complex spherical harmonics
"""
function cYlm!(Y, L, S::PseudoSpherical, P)
	@assert length(P) >= sizeP(L)
	@assert length(Y) >= sizeY(L)
   @assert abs(S.cosθ) <= 1.0

	ep = 1 / sqrt(2) + im * 0
	for l = 0:L
		Y[index_y(l, 0)] = P[index_p(l, 0)] * ep
	end

   sig = 1
   ep_fact = S.cosφ + im * S.sinφ
	for m in 1:L
		sig *= -1
		ep *= ep_fact            # ep =   exp(i *   m  * φ)
		em = sig * conj(ep)      # ep = ± exp(i * (-m) * φ)
		for l in m:L
			p = P[index_p(l,m)]
			@inbounds Y[index_y(l, -m)] = em * p   # (-1)^m * p * exp(-im*m*phi) / sqrt(2)
			@inbounds Y[index_y(l,  m)] = ep * p   #          p * exp( im*m*phi) / sqrt(2)
		end
	end

	return Y
end

"""
evaluate gradients of complex spherical harmonics
"""
function cYlm_d!(Y, dY, L, S::PseudoSpherical, P, dP)
	@assert length(P) >= sizeP(L)
	@assert length(Y) >= sizeY(L)
	@assert length(dY) >= sizeY(L)
   # @assert abs(S.cosθ) < 1.0

	# m = 0 case
	ep = 1 / sqrt(2)
	for l = 0:L
		Y[index_y(l, 0)] = P[index_p(l, 0)] * ep
		dY[index_y(l, 0)] = dspher_to_dcart(S, 0.0, dP[index_p(l, 0)] * ep)
	end

   sig = 1
   ep_fact = S.cosφ + im * S.sinφ

	for m in 1:L
		sig *= -1
		ep *= ep_fact            # ep =   exp(i *   m  * φ)
		em = sig * conj(ep)      # ep = ± exp(i * (-m) * φ)
		dep_dφ = im *   m  * ep
		dem_dφ = im * (-m) * em
		for l in m:L
			p_div_sinθ = P[index_p(l,m)]
			@inbounds Y[index_y(l, -m)] = em * p_div_sinθ * S.sinθ
			@inbounds Y[index_y(l,  m)] = ep * p_div_sinθ * S.sinθ

			dp_dθ = dP[index_p(l,m)]
			@inbounds dY[index_y(l, -m)] = dspher_to_dcart(S, dem_dφ * p_div_sinθ, em * dp_dθ)
			@inbounds dY[index_y(l,  m)] = dspher_to_dcart(S, dep_dφ * p_div_sinθ, ep * dp_dθ)
		end
	end
	# return Y, dY
end

"""
evaluate real spherical harmonics
"""
function rYlm!(Y::AbstractVector{T}, L, S::PseudoSpherical, P) where {T <: Real}
	@assert length(P) >= sizeP(L)
	@assert length(Y) >= sizeY(L)
   @assert abs(S.cosθ) <= 1.0

   oort2 = 1 / sqrt(2)
	for l = 0:L
		Y[index_y(l, 0)] = P[index_p(l, 0)] * oort2
	end

   sig = 1
	ec = 1.0 + 0 * im
   ec_fact = S.cosφ + im * S.sinφ
	for m in 1:L
		sig *= -1                # sig = (-1)^m
		ec *= ec_fact            # ec = exp(i * m  * φ) / sqrt(2)
		# cYlm = p * ec,    (also cYl{-m} = sig * p * conj(ec)), but not needed)
		# rYlm    =  Re(cYlm)
		# rYl{-m} = -Im(cYlm)
		for l in m:L
			p = P[index_p(l,m)]
			@inbounds Y[index_y(l, -m)] = -p * imag(ec)
			@inbounds Y[index_y(l,  m)] =  p * real(ec)
		end
	end

	return Y
end

"""
evaluate gradients of real spherical harmonics
"""
function rYlm_d!(Y::AbstractVector{T}, dY,
					  L, S::PseudoSpherical, P, dP) where {T <: Real}
	@assert length(P) >= sizeP(L)
	@assert length(Y) >= sizeY(L)
   @assert abs(S.cosθ) <= 1.0

   oort2 = 1 / sqrt(2)
	for l = 0:L
		Y[index_y(l, 0)] = P[index_p(l, 0)] * oort2
		dY[index_y(l, 0)] = dspher_to_dcart(S, 0.0, dP[index_p(l, 0)] * oort2)
	end

   sig = 1
	ec = 1.0 + 0 * im
   ec_fact = S.cosφ + im * S.sinφ
	drec_dφ = 0.0
	for m in 1:L
		sig *= -1                # sig = (-1)^m
		ec *= ec_fact            # ec = exp(i * m  * φ) / sqrt(2)
		dec_dφ = im * m * ec

		# cYlm = p * ec,    (also cYl{-m} = sig * p * conj(ec)), but not needed)
		# rYlm    =  Re(cYlm)
		# rYl{-m} = -Im(cYlm)
		for l in m:L
			p_div_sinθ = P[index_p(l,m)]
			p = p_div_sinθ * S.sinθ
			Y[index_y(l, -m)] = -p * imag(ec)
			Y[index_y(l,  m)] =  p * real(ec)

			dp_dθ = dP[index_p(l,m)]
			dY[index_y(l, -m)] = dspher_to_dcart(S, - imag(dec_dφ) * p_div_sinθ,
															    - imag(ec) * dp_dθ)
			dY[index_y(l,  m)] = dspher_to_dcart(S,   real(dec_dφ) * p_div_sinθ,
															      real(ec) * dp_dθ)
		end
	end

	return dY
end

# ---------------------------------------------
#      The nice basis interface
# ---------------------------------------------

abstract type AbstractSHBasis{T} <: IPBasis end

"""
complex spherical harmonics
"""
struct SHBasis{T} <: AbstractSHBasis{T}
	maxL::Int
	coeff::ALPCoefficients{T}
end

"""
real spherical harmonics
"""
struct RSHBasis{T} <: AbstractSHBasis{T}
	maxL::Int
	coeff::ALPCoefficients{T}
end

import Base.==
==(B1::AbstractSHBasis, B2::AbstractSHBasis) =
	( (B1.maxL == B1.maxL) &&
	  (typeof(B1) == typeof(B2)) )

write_dict(SH::SHBasis{T}) where {T} =
		Dict("__id__" => "ACE_SHBasis",
			  "T" => write_dict(T),
			  "maxL" => SH.maxL)

read_dict(::Val{:SHIPs_SHBasis}, D::Dict) =
	read_dict(Val{:ACE_SHBasis}(), D)

read_dict(::Val{:ACE_SHBasis}, D::Dict) =
		SHBasis(D["maxL"], read_dict(D["T"]))

SHBasis(maxL::Integer, T::Type=Float64) =
		SHBasis(Int(maxL), compute_coefficients(maxL, T))

RSHBasis(maxL::Integer, T::Type=Float64) =
		RSHBasis(Int(maxL), compute_coefficients(maxL, T))

rfltype(SH::AbstractSHBasis{T}) where {T} = T
Base.length(S::AbstractSHBasis) = sizeY(S.maxL)

fltype(::SHBasis{T}) where {T} = Complex{T}
fltype(::RSHBasis{T}) where {T} = T

alloc_B( S::SHBasis{T}, args...) where {T} =
		Vector{Complex{T}}(undef, length(S))

alloc_B( S::RSHBasis{T}, args...) where {T} =
		Vector{T}(undef, length(S))

alloc_dB(S::SHBasis{T}) where {T} =
		Vector{JVec{Complex{T}}}(undef, length(S))

alloc_dB(S::RSHBasis{T}) where {T} =
		Vector{JVec{T}}(undef, length(S))

alloc_dB(S::AbstractSHBasis, N::Integer) = alloc_dB(S)

alloc_temp(SH::AbstractSHBasis{T}, args...) where {T} = (
		P = Vector{T}(undef, sizeP(SH.maxL)), )

alloc_temp_d(SH::AbstractSHBasis{T}, args...) where {T} = (
		 P = Vector{T}(undef, sizeP(SH.maxL)),
		dP = Vector{T}(undef, sizeP(SH.maxL)) )

_evaluate!(Y, L, S, P, ::SHBasis) = cYlm!(Y, L, S, P)
_evaluate!(Y, L, S, P, ::RSHBasis) = rYlm!(Y, L, S, P)

_evaluate_d!(Y, dY, L, S, P, dP, ::SHBasis) = cYlm_d!(Y, dY, L, S, P, dP)
_evaluate_d!(Y, dY, L, S, P, dP, ::RSHBasis) = rYlm_d!(Y, dY, L, S, P, dP)

function evaluate!(Y, tmp, SH::AbstractSHBasis, R::JVec)
	L=SH.maxL
	@assert 0 <= L <= SH.maxL
	@assert length(Y) >= sizeY(L)
	S = cart2spher(R)
	compute_p!(L, S, SH.coeff, tmp.P)
	_evaluate!(Y, L, S, tmp.P, SH)
	return Y
end

function evaluate_d!(Y, dY, tmp, SH::AbstractSHBasis, R::JVec)
	L=SH.maxL
	@assert 0 <= L <= SH.maxL
	@assert length(Y) >= sizeY(L)
	# if R[1]^2+R[2]^2 < 1e-20 * R[3]^2
	# 	R = JVec(R[1]+1e-9, R[2], R[3])
	# end
	S = cart2spher(R)
	compute_dp!(L, S, SH.coeff, tmp.P, tmp.dP)
	_evaluate_d!(Y, dY, L, S, tmp.P, tmp.dP, SH)
	# return Y, dY
end

end


text:
function create_subgraph(vertex_set::AbstractArray,edge_list::AbstractArray,directed::Bool)
  new_graph = graph(Int64[],ExEdge{Int64}[],is_directed = directed)
  vertex_set_d = Dict()
  for vertex in vertex_set
    add_vertex!(new_graph,vertex)
    vertex_set_d[vertex] = 1
  end
  for edge in edge_list
    if (haskey(vertex_set_d,edge.source)) && (haskey(vertex_set_d,edge.target))
      new_edge = ExEdge(edge.index,edge.source,edge.target)
      for (k,v) in edge.attributes
        new_edge.attributes[k] = v
      end
        add_edge!(new_graph,new_edge)
    end
  end
  return new_graph
end

function all_neighbors(v::Int64,g::GenericGraph)
  if is_directed(g)
    neighbr = out_neighbors(v,g)
    for n in in_neighbors(v,g)
      append!(neighbr,n)
    end
    return neighbr
  else
    return out_neighbors(v,g)
  end
end

function all_degree(v::Int64,g::GenericGraph)
  if is_directed(g)
    return out_degree(v,g) + in_degree(v,g)
  else
    return out_degree(v,g)
  end
end

function occurrence_of_degree(g::GenericGraph)
  g_vertices = vertices(g)
  degrees = zeros(Int64,0)
  max = 0
  d = 0
  for i = 1:length(g_vertices)
    d = all_degree(g_vertices[i],g)+1
    if max < d
      for j=1:(d-max)
        push!(degrees,0)
      end
      degrees[d] +=1
      max = d
    else
      degrees[d] +=1
    end
  end
  return degrees
end

max_stars_repo_name: johnybx/GraphsExt

max_issues_repo_name: johnybx/GraphsExt

max_forks_repo_name: johnybx/GraphsExt


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

JLLWrappers.@generate_main_file_header("GAP_pkg_browse")
JLLWrappers.@generate_main_file("GAP_pkg_browse", UUID("6e540885-1cbc-538c-85e7-b800ca60d645"))
end  # module GAP_pkg_browse_jll


text:
## perform the Cauchy radial fit on the ocs
##

rootdir= ENV["GAIA_ROOT"]
push!(LOAD_PATH,"$rootdir/run/src")
using GaiaClustering
using Printf, Glob, Statistics

import PyPlot , CSV , DataFrames

## directory
wdir    = "$rootdir/e2e_products"
plotdir = "$wdir/plots-individuals"
ocdir=    "$wdir/oc"

cd(wdir)

let

## fit parameters##########################################################
       NTEST= 100 ; NITER= 30000 ; NBURNIN= 5000 ##MCMC parameters for the fit
       OUTFILE= "votlist-CauchyFit.csv"
###########################################################################
    @printf("## Compute the Cauchy radial fit on the OCs\n")
    @printf("##\n")
############################functions #####################################
function plotFitDensityFig(filelist, θlist, ifig, nrow, ncol, ocdir)
    fig= PyPlot.figure(figsize=(15,10))

    for i in 1:length(filelist)
        ocname= filelist[i][1:end-4]
        θ= θlist[i]
        ocfile= "$ocdir/$ocname-oc.csv"
        plotDensity(ocfile, θ,nrow, ncol, i, ocname)
    end

    figname= "densityOC-$ifig.png"
    PyPlot.savefig(figname)
    PyPlot.close(fig)
end

function plotDensity(ocfile, θ, irow, icol, nplot, title)
    oc= CSV.read(ocfile, delim= ";")
    # binned density
    nstar= length(oc.X)
    nbin= min(trunc(Int, nstar/10),20)
    r2d,ρ2d,err2d= density2D(oc.Y, oc.Z, nbin)
    ρ2dfit= model_rad(r2d, θ, fdens1)

    ax= PyPlot.subplot(irow, icol, nplot)
    ax.set_xlabel("r [pc]") ; ax.set_ylabel("ρ")
    ax.set_yscale("log") ; ax.set_xscale("log")
    ax.set_xlim(r2d[1]*0.9, r2d[end]*1.1)
    ax.set_ylim(minimum(ρ2d[ρ2d .> 0])*0.3,maximum(ρ2d)*1.5)
    PyPlot.grid("on")
    PyPlot.scatter(r2d, ρ2d , s=4, facecolor="blue" )
    PyPlot.errorbar(r2d, ρ2d, yerr=2 .* err2d, linewidth=0.5)
    PyPlot.plot(r2d, ρ2dfit, "k-", linewidth=1)

    ## TEXT
    ax.text(0.02, 0.01, title, verticalalignment= "bottom" , horizontalalignment="left",
        transform=ax.transAxes, color="green", fontsize=15)
    val=θ.C ; txt= @sprintf("C= %3.3f",val)
    ax.text(0.02, 0.2, txt, verticalalignment= "bottom" , horizontalalignment="left",
        transform=ax.transAxes, color="red", fontsize=12)
    val=θ.s ; txt= @sprintf("s= %3.3f [pc]",val)
    ax.text(0.02, 0.3, txt, verticalalignment= "bottom" , horizontalalignment="left",
        transform=ax.transAxes, color="red", fontsize=12)
    val=θ.m ; txt= @sprintf("m= %3.3f",val)
    ax.text(0.02, 0.4, txt, verticalalignment= "bottom" , horizontalalignment="left",
        transform=ax.transAxes, color="red", fontsize=12)
    val= nstar ; txt= @sprintf("Nstar= %d",val)
    ax.text(0.02, 0.52, txt, verticalalignment= "bottom" , horizontalalignment="left",
        transform=ax.transAxes, color="blue", fontsize=12)
end

## main loop
    cd(ocdir)
    files= glob("*-oc.csv")
    cd(wdir)

    df = DataFrames.DataFrame(votname= String[], C= Float64[], Cerr= Float64[], s= Float64[], serr= Float64[], m= Float64[],
merr= Float64[])
    i= 0 ; noc= 0 ; nfail= 0
    for f in files
        i+= 1
        ocfile= "$ocdir/$f"
        oc= CSV.read("$ocdir/$f" , delim= ";")
        println("### ocfile: $f")
        println("### fitting the radial profile....")
        votname= @sprintf("%s.vot", f[1:end-7])

        nstar= length(oc.X)
        nbin= min(trunc(Int, nstar/10),20)

        verbose= false
        θfit, θfiterr, fitfound= spatialParameter(ocfile, NTEST, nbin, NITER, verbose, NBURNIN)
        println("### Fit:")
        println(θfit)
        println(θfiterr)

        push!(df, [votname,θfit.C , θfiterr.C, θfit.s , θfiterr.s, θfit.m , θfiterr.m])

        noc += 1
        if !fitfound
            println("################################ Failing...########### ")
            println("###################################################### ")
            nfail += 1
        end
        println("## OCs analyzed: $noc")
        println("## Fit failed: $nfail \n\n")
    end
    CSV.write(OUTFILE, df, delim=";")

    ######################## plots
    println("## Plotting the individual fits...")
    cd(plotdir)
    rowperfig= 4 ; colperfig= 3 ; plotperfig= rowperfig*colperfig

    ifig= 1 ; plotaccum=0
    filevot= [] ; θlist= [] ; θerrlist= []
    for item in 1:length(df.votname)
        push!(filevot, df.votname[item])
        θ= GaiaClustering.modelCauchy(df.C[item], df.s[item], df.m[item])
        push!(θlist, θ)
        θerr= GaiaClustering.modelCauchy(df.Cerr[item], df.serr[item], df.merr[item])
        push!(θerrlist, θerr)
        plotaccum += 1

        if plotaccum == plotperfig
            plotFitDensityFig(filevot, θlist, ifig,rowperfig, colperfig , ocdir)
            println("## Fig $ifig done.")
            filevot= [] ; θlist= [] ; θerrlist= []
            plotaccum= 0 ; ifig += 1
        end
    end
    if plotaccum > 0
        plotFitDensityFig(filevot, θlist, ifig,rowperfig, colperfig , ocdir)
        println("## Fig $ifig done.")
    end
    cd(wdir)
end

max_stars_repo_name: bosscha/gaia-shock

max_issues_repo_name: bosscha/gaia-shock

max_forks_repo_name: bosscha/gaia-shock


text:
export gridData,
       ROIData,
       generateInputData,
       GridData,
       createGrid,
       vtkUnsGridPlot,
       vtkStrGridPlot,
       computeFEM2GridInterpMatrix


text:
# this file will build Pumi if it cannot be located by pkg-config

# install PkgFix if not present
if !isdir(joinpath(Pkg.dir(), "PkgFix"))
  Pkg.clone("https://github.com/OptimalDesignLab/PkgFix.jl.git")
end
start_dir=pwd()
cd(Pkg.dir("PkgFix"))
run(`git checkout upgrade_0.6`)
cd(start_dir)

using PkgFix  # from now on, use PkgFix instead of Pkg for everything

include("./stringmatch.jl")
include("install_pumi.jl")

# package URLs and versions

global const ARRAYVIEWS_URL = "https://github.com/JaredCrean2/ArrayViews.jl.git"
global const ARRAYVIEWS_VER = "work"

global const MPI_URL = "https://github.com/JuliaParallel/MPI.jl.git"
global const MPI_VER = "v0.5.0"

global const SBP_URL = "https://github.com/OptimalDesignLab/SummationByParts.jl.git"
global const SBP_VER = "jcwork"

start_dir = pwd()

run(`./cleanup.sh`)

install_pumi = !haskey(ENV, "SCOREC_PREFIX") || try run(`cmake --find-package -DNAME=SCOREC -DCOMPILER_ID=GNU -DLANGUAGE=CXX -DMODE=EXIST`)
          false
          catch
	  true
	end
println("install_pumi = ", install_pumi)

install_mpi = try run(`which mpicxx`)
          false
          catch
	  true
	end
println("install_mpi = ", install_mpi)

# install MPI itself
if install_mpi
  cmd_string = "./travis-install-mpi.sh"
  arg_str = "mpich3"
  run(`$cmd_string $arg_str`)
end

#=
# install MPI.jl if needed
if PkgFix.installed("MPI") == nothing
  PkgFix.add(MPI_URL, branch_ish=MPI_VER)
end
=#

if install_pumi  # did not find pumi
  installPumi()
end

# install non-metadata dependencies
pkg_dict = PkgFix.installed()  # get dictionary of installed package names to version numbers

if !haskey(pkg_dict, "SummationByParts")
  PkgFix.add(SBP_URL, branch_ish=SBP_VER)
  # SBP installs ODLCommonTools
end

if !haskey(pkg_dict, "ArrayViews")
  PkgFix.add(ARRAYVIEWS_URL, branch_ish=ARRAYVIEWS_VER)
end

# now that all dependencies exist, install this package
cd("../src")
run(`./config.sh`)
run(`./makeinstall.sh`)
cd(start_dir)

#ldd ./libfuncs1.so | grep libmds.so | awk -F"=>" '{print $2}' | awk '{print $1}' |  sed 's%/[^/]*$%/%'


text:
using Test

@testset "Number conversations" begin
    include("utils.jl")
end

@testset "Binary tree parser tests" begin
    include("parser.jl")
end

@testset "Random number generators and oracles" begin
    include("rho.jl")
end

@testset "Groups, generators" begin
    include("generators.jl")
end

@testset "Independent basis vectors" begin
    include("crs.jl")
end

@testset "ElGamal test" begin
    include("ciphertexts.jl")
end

@testset "Bare Verifier" begin
    include("proofparser.jl")
end


text:
# imresize.jl

using ImageGeoms: downsample, downsample!
using ImageGeoms: upsample, upsample!
#include("../src/imresize.jl")

using Test: @test, @testset, @test_throws, @inferred

@testset "downsample" begin
    x = trues(6, 8)
    y = @inferred downsample(x, (3,2))
    @test y == trues(2, 4)

    x = ones(Bool, 6, 8)
    y = @inferred downsample(x, (3,2))
    @test y == trues(2, 4)

    x = trues(6, 9)
    z = @view x[1:6,1:8]
    y = @inferred downsample(z, (3,2))
    @test y == trues(2, 4)

    x = trues(6, 8, 10, 5)
    y = @inferred downsample(x, (3,2,5,5))
    @test y == trues(2, 4, 2, 1)

    x = trues(6, 7)
    y = @inferred downsample(x, (2,2))
    @test y == trues(3, 4)

    x = ones(Bool, 6, 7)
    y = @inferred downsample(x, (2,2))
    @test y == trues(3, 4)

    x = Bool[0 0 0 0 0; 0 1 1 0 0; 0 1 0 0 0; 0 0 0 0 0]
    y = @inferred downsample(x, (2,2))
    @test y == Bool[1 1 0; 1 0 0]

    x = ones(Bool, 6, 7)
    y = falses(2,2)
    @test_throws DimensionMismatch downsample!(y, x, (2,2))
end

@testset "upsample" begin
    x = trues(2, 3)
    y = @inferred upsample(x, (5,4))
    @test y == trues(10, 12)

    x = Bool[1 0 1; 0 1 0]
    y = @inferred upsample(x, (2,2))
    @test y == Bool[1 1 0 0 1 1; 1 1 0 0 1 1; 0 0 1 1 0 0; 0 0 1 1 0 0]

    x = ones(Bool, 2, 3)
    y = @inferred upsample(x, (4,1))
    @test y == trues(8, 3)

    @test_throws DimensionMismatch upsample!(y, x, (2,2))
end


text:
function class_group_proof(clg::ClassGrpCtx, lb::fmpz, ub::fmpz; extra :: fmpz=fmpz(0), prec::Int = 100, do_it=1:ub)
  #for all prime ideals P with lb <= norm <= ub, find a relation
  #tying that prime to the factor base
  # if extra is useful, assume that the function was already run for all primes
  # up to norm extra

  if extra==0
    extra = norm(clg.FB.ideals[1])
  end
  lb = max(lb, norm(clg.FB.ideals[1]))
  #println("expect to need ", Int(floor(li(ub*1.0) - li(lb*1.0))), " primes")
  O = order(clg.FB.ideals[1])
  n = degree(O)
  p = next_prime(root(lb, n))
  np = Int(floor(log(abs(discriminant(O)))/log(2)/2))
  no_primes = 0
  no_ideals = 0
  if do_it.start > 1
    p = fmpz(next_prime(do_it.start))
  end
  r = fmpz()
  #gc_enable(false)
  while p < do_it.stop
    no_primes += 1
    if no_primes % 100 == 0
      println("did $no_primes prime numbers so far, now $p, need to reach $ub")
    end
    deg_lim = Int(floor(log(ub)/log(p)))
    low_lim = Int(floor(log(lb)/log(p)))
    fac = prime_decomposition(O, Int(p), deg_lim, low_lim)
    for _k in fac
      k = _k[1]
      if norm(k) <= lb
        continue
      end
      no_ideals += 1
      #if no_ideals % 10 == 0
      #  println("done $no_ideals ideals so far...")
      #  x = Base.gc_num()
      #  println("alloc $(x.malloc)   free $(x.freed)  diff: $(x.malloc - x.freed)")
      #  gc_enable(true)
      #  gc()
      #  gc_enable(false)
      #  x = Base.gc_num()
      #  println("alloc $(x.malloc)   free $(x.freed)  diff: $(x.malloc - x.freed)")
      #end
      #println("to be more precise: $k")
      E = class_group_small_real_elements_relation_start(clg, k, limit=10, prec=prec)
      while true
        sucess = false
        a = class_group_small_real_elements_relation_next(E)
        n = norm_div(a, norm(k), np)
        if gcd(numerator(n), p) > extra
#          println("a: $a, $(norm(a)), $(norm(k)), $n")
#          println("contains too many conjugates, bad")
          continue
        end
        f, r = issmooth!(clg.FB.fb_int, numerator(n))
        if f
          M = SMat{Int}()
          fl = _factor!(clg.FB, a, false, n)[1]
          if fl
            break
          else
#            println("not smooth, ideal")
          end
        else
#          println("not smooth, int")
        end
      end
    end
    p = next_prime(p)
  end
  #println("success: used $no_primes numbers and $no_ideals ideals")
  #gc_enable(true)
end


text:
function test_projection(order)
    if order == 1
        grid = JuAFEM.generate_grid(Quadrilateral, (1, 1), Vec((0.,0.)), Vec((1.,1.)))
    elseif order == 2
        # grid = JuAFEM.generate_grid(QuadraticQuadrilateral, (1, 1), Vec((0.,0.)), Vec((1.,1.)))
        grid = JuAFEM.generate_grid(Quadrilateral, (1, 1), Vec((0.,0.)), Vec((1.,1.)))
    end

    dim = 2
    ip = Lagrange{dim, RefCube, order}()
    ip_geom = Lagrange{dim, RefCube, 1}()
    qr = QuadratureRule{dim, RefCube}(order+1)
    cv = CellScalarValues(qr, ip, ip_geom)

    # Create node values for the cell
    f(x) = Tensor{1,1,Float64}((1 + x[1]^2 + (2x[2])^2, ))
    xe = getcoordinates(grid, 1)
    # analytical values
    qp_values = [[f(spatial_coordinate(cv, qp, xe)) for qp in 1:getnquadpoints(cv)]]

    # Now recover the nodal values using a L2 projection. Since f is quadratic and the interpolation as well, we should recover the exact nodal values
    projector = L2Projector(cv, ip, grid)

    point_vars = JuAFEM.project(qp_values, projector)

    ae = compute_vertex_values(grid, f)
    # The projection gives the values in node order -> reorder ae
    # @test point_vars[1:4] ≈ [ae[1], ae[2], ae[4], ae[3]]
    # return point_vars[1:4], [ae[1], ae[2], ae[4], ae[3]]
    return point_vars[1:4], [ae[1], ae[2], ae[3], ae[4]]
end

function test_projection_mixedgrid()
    # generate a mesh with 1 quadrilateral and 2 triangular elements
    dim = 2
    nodes = Node{dim, Float64}[]
    push!(nodes, Node((0.0, 0.0)))
    push!(nodes, Node((1.0, 0.0)))
    push!(nodes, Node((2.0, 0.0)))
    push!(nodes, Node((0.0, 1.0)))
    push!(nodes, Node((1.0, 1.0)))
    push!(nodes, Node((2.0, 1.0)))

    cells = JuAFEM.AbstractCell[]
    push!(cells, Quadrilateral((1,2,5,4)))
    push!(cells, Triangle((2,3,6)))
    push!(cells, Triangle((2,6,5)))

    mesh = Grid(cells, nodes)

    order = 2
    ip = Lagrange{dim, RefCube, order}()
    ip_geom = Lagrange{dim, RefCube, 1}()
    qr = QuadratureRule{dim, RefCube}(order+1)
    cv = CellScalarValues(qr, ip, ip_geom)

    # Create node values for the 1st cell
    # use a SymmetricTensor here for testing the symmetric version of project
    f(x) = SymmetricTensor{2,2,Float64}((1 + x[1]^2, 2x[2]^2, x[1]*x[2]))
    xe = getcoordinates(mesh, 1)
    ae = compute_vertex_values(mesh, f)
    # analytical values
    qp_values = [[f(spatial_coordinate(cv, qp, xe)) for qp in 1:getnquadpoints(cv)]]

    # Now recover the nodal values using a L2 projection.
    # Assume f would only exist on the first cell, we project it to the nodes of the
    # 1st cell while ignoring the rest of the domain. NaNs should be stored in all
    # nodes that do not belong to the 1st cell
    projector = L2Projector(cv, ip, mesh, 1:1)
    point_vars = JuAFEM.project(qp_values, projector)

    # In the nodes of the 1st cell we should recover the field
    for node in mesh.cells[1].nodes
        @test ae[node] ≈ point_vars[node]
    end

    # in all other nodes we should have NaNs
    for node in setdiff(1:getnnodes(mesh), mesh.cells[1].nodes)
        for d1 = 1:dim, d2 = 1:dim
             @test isnan(point_vars[node][d1, d2])
         end
    end
end

function test_projection_newprojector(order)
    if order == 1
        grid = JuAFEM.generate_grid(Quadrilateral, (1, 1), Vec((0.,0.)), Vec((1.,1.)))
    elseif order == 2
        # grid = JuAFEM.generate_grid(QuadraticQuadrilateral, (1, 1), Vec((0.,0.)), Vec((1.,1.)))
        grid = JuAFEM.generate_grid(Quadrilateral, (1, 1), Vec((0.,0.)), Vec((1.,1.)))
    end

    dim = 2
    ip = Lagrange{dim, RefCube, order}()
    ip_geom = Lagrange{dim, RefCube, 1}()
    qr = QuadratureRule{dim, RefCube}(order+1)
    cv = CellScalarValues(qr, ip, ip_geom)

    # Create node values for the cell
    f(x) = Tensor{1,1,Float64}((1 + x[1]^2 + (2x[2])^2, ))
    xe = getcoordinates(grid, 1)
    # analytical values
    qp_values = [[f(spatial_coordinate(cv, qp, xe)) for qp in 1:getnquadpoints(cv)]]

    # Now recover the nodal values using a L2 projection. Since f is quadratic and the interpolation as well, we should recover the exact nodal values
    projector = L2Projector(qr, ip, grid)

    point_vars = JuAFEM.project(qp_values, projector)

    ae = compute_vertex_values(grid, f)
    # The projection gives the values in node order -> reorder ae
    # @test point_vars[1:4] ≈ [ae[1], ae[2], ae[4], ae[3]]
    # return point_vars[1:4], [ae[1], ae[2], ae[4], ae[3]]
    return point_vars[1:4], [ae[1], ae[2], ae[3], ae[4]]
end

@testset "Test L2-Projection" begin
    # Tests a L2-projection of integration point values (to nodal values), determined from the function y = 1 + x[1]^2 + (2x[2])^2

    # A linear approximation can not recover a quadratic solution, so projected values will be different from the analytical ones
    projected_vars, analytical_vars = test_projection(1)
    @test projected_vars ≉  analytical_vars

    # For a quadratic approximation the analytical solution is recovered
    projected_vars, analytical_vars = test_projection(2)
    @test projected_vars ≈ analytical_vars

    # Test a mixed grid, where only a subset of the cells contains a field
    test_projection_mixedgrid()

end


text:
using Game2048
using Serialization: serialize
using Game2048: value, play_via_monte_carlo_w_recorder, play_via_monte_carlo_wo_recorder

# @time x = [play_via_monte_carlo_wo_recorder(randompolicy; n=1000) for _ in 1:10]

maximum.(x)

while true
    results = play_n_game(randompolicy;  ntries_per_turn=1000)
    objhash=hash(results)
    serialize("c:/data/game2048/raw/game2048_$objhash.jls", results)
end

# find the score

scores = map(x->value(x.final_state), results)
fnlmax = map(x->maximum(x.final_state), results)

using DataFrames, Chain, DataFrameMacros, Statistics

@chain DataFrame(;scores, fnlmax) begin
    groupby(:fnlmax)
    @combine(mean(:scores))
end

y = map(results) do x
    maximum(x.final_state)
end

countmap(y)

@time results_cnn = play_n_game(policy, 100)

@time play_n_game(policy, 1)

y_cnn = map(results_cnn) do x
    maximum(x.final_state)
end

countmap(y_cnn)

using Alert
alert()


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__()

module Docile

"""
!!summary(Documentation extraction package for the Julia Language.)

Docile supports several styles of docstring. Namely:

The original ``@doc`` macro, which is now also available in Julia ``0.4``:

    using Docile

    @doc \"\"\"
    ...
    \"\"\" ->
    f(x) = x

Plain strings using the ``@document`` macro. Requires importing ``Docile`` into
the module:

    using Docile
    @document

    \"\"\"
    ...
    \"\"\"
    f(x) = x

Automatic extraction for plain docstrings from packages that have been
``import``ed:

    \"\"\"
    ...
    \"\"\"
    f(x) = x

The third style is the newest, and preferred, way to document packages in a
Docile-compatible way. ``@doc`` and ``@document`` docstring usage will remain
part of Docile for the foreseeable future to allow for backward compatibility.
"""
Docile

include("Utilities.jl")                          # Code useful across submodules.

include(joinpath("Formats", "Formats.jl"))       # Format-dispatch for parsing docstrings.

include(joinpath("Runner", "Runner.jl"))         # Execute expressions.
include(joinpath("Collector", "Collector.jl"))   # Collect information from modules.

include(joinpath("Cache", "Cache.jl"))           # Store collected information.

include(joinpath("Legacy", "Legacy.jl"))         # Compatibility with base and docile.
include(joinpath("Interface", "Interface.jl"))   # Public interface to the package.

include(joinpath("Extensions", "Extensions.jl")) # Additional formatters and parsers.

## Direct Exports. ##

import .Legacy:

    @doc, @comment,

    @document, @docstrings,

    @doc_str, @doc_mstr, @file_str, @md_str, @md_mstr,

    meta

export

    @doc, @comment,

    @document, @docstrings,

    @doc_str, @doc_mstr, @file_str, @md_str, @md_mstr,

    meta,

    Interface

end


text:
# following Asymptotic Gaussian law for
# noninteracting indistinguishable
# particles in random networks
# Valery S. Shchesnovich

"""
    partition_expectation_values(partition_size_vector, partition_counts)
    partition_expectation_values(part_occ::PartitionOccupancy)

Return the Haar averaged probability of photon number count in binned outputs
for [`Distinguishable`](@ref) and [`Bosonic`](@ref) particles.

!!! note "Reference"
    [https://www.nature.com/articles/s41598-017-00044-8.pdf](https://www.nature.com/articles/s41598-017-00044-8.pdf)
"""
function partition_expectation_values(partition_size_vector, partition_counts)

    m = sum(partition_size_vector)
    number_partitions = length(partition_size_vector)
    partition_size_ratio = partition_size_vector ./ m # q
    n = sum(partition_counts)
    @test all(partition_counts .>= 0)

    # factorials can quickly get pretty big so this is required
    if any(partition_counts .> 20) || n > 20
        n = big(n)
        partition_counts = big.(partition_counts)
    end

    proba_dist = factorial(n) / vector_factorial(partition_counts) * prod(partition_size_ratio .^ partition_counts)

    proba_bosonic = proba_dist * prod([partition_counts[i] > 1 ? prod([1+l/partition_size_vector[i] for l in 1:partition_counts[i]-1]) : 1 for i in 1:length(partition_size_vector)]) /prod([1+l/m for l in 1:n-1])

    proba_dist, proba_bosonic

end

partition_expectation_values(part_occ::PartitionOccupancy) = partition_expectation_values(partition_occupancy_to_partition_size_vector_and_counts(part_occ)...)

"""
    subset_expectation_value(subset_size, k, n, m)

Return the Haar averaged probability to find `k` from `n` photons inside a subset of
binned output modes of length `size_subset` among `m` modes for [`Distinguishable`](@ref)
and [`Bosonic`](@ref) cases.
"""
function subset_expectation_value(subset_size, k,n,m)

    partition_size_vector = [subset_size, m-subset_size]
    partition_counts = [k,n-k]

    partition_expectation_values(partition_size_vector, partition_counts)

end

"""
    subset_relative_distance_of_averages(subset_size, n, m)
"""
function subset_relative_distance_of_averages(subset_size,n,m)

    @warn "check tvd conventions"
    0.5*abs(sum([abs(subset_expectation_value(subset_size, k,n,m)[1] - subset_expectation_value(subset_size, k,n,m)[2]) for k in 0:n]))

end

"""
    choose_best_average_subset(;m, n, distance=tvd)

Return the ideal subset size on average and its total variance distance.
"""
function choose_best_average_subset(;m,n, distance = tvd)

    """returns the ideal subset size on average and its TVD,
    to compare with the full bunching test of shsnovitch"""

    function distance_this_subset_size(subset_size)

        proba_dist = [subset_expectation_value(subset_size,k,n,m)[1] for k in 0:n]
        proba_bos = [subset_expectation_value(subset_size,k,n,m)[2] for k in 0:n]

        distance(proba_dist, proba_bos)

    end

    max_distance = 0
    subset_size_max_ratio = nothing

    for subset_size in 1:m-1

        if distance_this_subset_size(subset_size) > max_distance
            max_distance = distance_this_subset_size(subset_size)
            subset_size_max_ratio = subset_size
        end
    end

    subset_size_max_ratio, distance_this_subset_size(subset_size_max_ratio)

end

"""
    best_partition_size(;m, n, n_subsets, distance=tvd)

Return the ideal `partition_size_vector` for a given number of subsets `n_subsets`
and the Haar averaged TVD in second parameter.

!!! note
    For a single subset, `n_subsets`=2 as we need a complete partition, occupying all modes.
"""
function best_partition_size(;m,n, n_subsets, distance = tvd)

    """return the ideal partition_size_vector for a given number
    of subsets n_subsets as well as the haar averaged tvd in second parameter

    for a single subset, n_subsets = 2 as we need a complete partition, occupying all modes"""

    @argcheck n_subsets >= 2 "we need a complete partition, occupying all modes"
    @argcheck n_subsets <= m "more partition bins than output modes"

    part_list = all_mode_configurations(m,n_subsets, only_photon_number_conserving = true)

    remove_trivial_partitions!(part_list)

    max_distance = 0
    part_max_ratio = nothing

    for part in ranked_partition_list(part_list)

        events = all_mode_configurations(n,n_subsets, only_photon_number_conserving = true)

        pdf = [partition_expectation_values(part, event) for event in events]

        pdf_dist = hcat(collect.(pdf)...)[1,:]
        pdf_bos = hcat(collect.(pdf)...)[2,:]

        this_distance = distance(pdf_bos,pdf_dist)

        println(part, " : ", this_distance)

        if this_distance > max_distance
            max_distance = this_distance
            part_max_ratio = part
        end
    end

    part_max_ratio, max_distance

end


text:
using Test, Polyhedra

# Taken from Base: test/show.jl
replstr(x, lim=true) = sprint((io,x) -> show(IOContext(io, :limit => lim, :displaysize => (24, 80)), MIME("text/plain"), x), x)
showstr(x) = sprint((io,x) -> show(IOContext(io, :limit => true, :displaysize => (24, 80)), x), x)

macro test_show(expr)
    esc(quote
            @test string($expr) == "$($expr)"
    end)
end

@testset "Show" begin
    @testset "Simple show" begin
        @test_show convexhull([1, 0]) + convexhull(Line([1, 1]))
        @test_show HyperPlane([1, 0], 0) ∩ HyperPlane([0, 1], 0) ∩ HalfSpace([1, 1], 1)
    end
    vshort = vrep([[i, -i] for i in 1:3])
    vlong = vrep([[i, -i] for i in 1:30])
    PT = Polyhedra.pointtype(vshort)
    RT = Polyhedra.raytype(vshort)
    @testset "Iterator" begin
        @test replstr(points(vshort))        == "3-element iterator of $PT:\n [1, -1]\n [2, -2]\n [3, -3]"
        @test replstr(points(vshort), false) == "3-element iterator of $PT:\n [1, -1]\n [2, -2]\n [3, -3]"
        @test replstr(points(vlong)) == "30-element iterator of $PT:\n [1, -1]\n [2, -2]\n [3, -3]\n [4, -4]\n [5, -5]\n [6, -6]\n [7, -7]\n [8, -8]\n [9, -9]\n [10, -10]\n [11, -11]\n [12, -12]\n [13, -13]\n [14, -14]\n [15, -15]\n [16, -16]\n [17, -17]\n [18, -18]\n [19, -19]\n [20, -20]\n  ⋮"
        @test replstr(points(vlong), false) == "30-element iterator of $PT:\n [1, -1]\n [2, -2]\n [3, -3]\n [4, -4]\n [5, -5]\n [6, -6]\n [7, -7]\n [8, -8]\n [9, -9]\n [10, -10]\n [11, -11]\n [12, -12]\n [13, -13]\n [14, -14]\n [15, -15]\n [16, -16]\n [17, -17]\n [18, -18]\n [19, -19]\n [20, -20]\n [21, -21]\n [22, -22]\n [23, -23]\n [24, -24]\n [25, -25]\n [26, -26]\n [27, -27]\n [28, -28]\n [29, -29]\n [30, -30]"
        @test showstr(points(vshort)) == string([[1, -1], [2, -2], [3, -3]])
        @test showstr(points(vlong)) == string([[i, -i] for i in 1:20])[1:end-1] * " … ]"
    end
    v = vshort + Ray([1, 0])
    p = polyhedron(v)
    PPT = Polyhedra.pointtype(p)
    PRT = Polyhedra.raytype(p)
    @testset "Polyhedron without H-representation" begin
        @test replstr(p)        == "Polyhedron $(typeof(p)):\n3-element iterator of $PPT:\n Rational{BigInt}[1//1, -1//1]\n Rational{BigInt}[2//1, -2//1]\n Rational{BigInt}[3//1, -3//1],\n1-element iterator of $PRT:\n Ray(Rational{BigInt}[1//1, 0//1])"
        @test replstr(p, false) == "Polyhedron $(typeof(p)):\n3-element iterator of $PPT:\n Rational{BigInt}[1//1, -1//1]\n Rational{BigInt}[2//1, -2//1]\n Rational{BigInt}[3//1, -3//1],\n1-element iterator of $PRT:\n Ray(Rational{BigInt}[1//1, 0//1])"
        @test showstr(p) == "convexhull([1//1, -1//1], [2//1, -2//1], [3//1, -3//1]) + convexhull(Ray(Rational{BigInt}[1//1, 0//1]))"
    end
    h = hrep(p)
    HS = Polyhedra.halfspacetype(p)
    @testset "V-Representation" begin
        @test replstr(v)        == "V-representation $(typeof(v)):\n3-element iterator of $PT:\n [1, -1]\n [2, -2]\n [3, -3],\n1-element iterator of $RT:\n Ray([1, 0])"
        @test replstr(v, false) == "V-representation $(typeof(v)):\n3-element iterator of $PT:\n [1, -1]\n [2, -2]\n [3, -3],\n1-element iterator of $RT:\n Ray([1, 0])"
        @test showstr(v) == "convexhull([1, -1], [2, -2], [3, -3]) + convexhull(Ray([1, 0]))"
    end
    it_expected = """
3-element iterator of $HS:
 HalfSpace(Rational{BigInt}[0//1, -1//1], 3//1)
 HalfSpace(Rational{BigInt}[0//1, 1//3], -1//3)
 HalfSpace(Rational{BigInt}[-1//1, -1//1], 0//1)"""
    expected = """
H-representation $(typeof(h)):
$it_expected"""
    short = "HalfSpace(Rational{BigInt}[0//1, -1//1], 3//1) ∩ HalfSpace(Rational{BigInt}[0//1, 1//3], -1//3) ∩ HalfSpace(Rational{BigInt}[-1//1, -1//1], 0//1)"
    @testset "H-Representation" begin
        @test replstr(h)        == expected
        @test replstr(h, false) == expected
        @test showstr(h) == short
    end
    @testset "Polyhedron" begin
        expected = """
Polyhedron $(typeof(p)):
$it_expected:
3-element iterator of $PPT:
 Rational{BigInt}[1//1, -1//1]
 Rational{BigInt}[2//1, -2//1]
 Rational{BigInt}[3//1, -3//1],
1-element iterator of $PRT:
 Ray(Rational{BigInt}[1//1, 0//1])"""
        @test replstr(p)        == expected
        @test replstr(p, false) == expected
        @test showstr(p) == "$short : convexhull([1//1, -1//1], [2//1, -2//1], [3//1, -3//1]) + convexhull(Ray(Rational{BigInt}[1//1, 0//1]))"
    end
    p = polyhedron(h)
    @testset "Polyhedron without V-representation" begin
        @test replstr(p)        == "Polyhedron $(typeof(p)):\n$it_expected"
        @test replstr(p, false) == "Polyhedron $(typeof(p)):\n$it_expected"
        @test showstr(p) == short
    end
end

max_stars_repo_name: blegat/Polyhedra

max_issues_repo_name: blegat/Polyhedra

max_forks_repo_name: blegat/Polyhedra


text:
module AutodiffTests

using Gridap.Arrays
using BlockArrays

function user_cell_energy(cell_u)
  function f(u)
    e = zero(eltype(u))
    for ui in u
      e += ui^2
    end
    e
  end
  lazy_map(f,cell_u)
end

function user_cell_residual(cell_u)
  function f(u)
    r = 2*u
    r
    #r = copy(u)
    #r .= zero(eltype(u))
    #for (i,ui) in enumerate(u)
    #  r[i] = 2*ui
    #end
    #r
  end
  lazy_map(f,cell_u)
end

function user_cell_jacobian(cell_u)
  function f(u)
    n = length(u)
    j = zeros(n,n)
    for i in 1:n
      j[i,i] = 2
    end
    j
  end
  lazy_map(f,cell_u)
end

L = 10
l = 24 # Do not use a number <13 (too easy for ForwardDiff)

blocksids = [(1,),(2,)]
axs = (blockedrange([4,4]),)
cell_u = [ BlockArrayCoo(axs,blocksids,[rand(4),rand(4)]) for i in 1:L ]

cell_e = user_cell_energy(cell_u)
cell_r = user_cell_residual(cell_u)
cell_j = user_cell_jacobian(cell_u)
cell_h = cell_j

cell_r_auto = autodiff_array_gradient(user_cell_energy,cell_u)
cell_j_auto = autodiff_array_jacobian(user_cell_residual,cell_u)
cell_h_auto = autodiff_array_hessian(user_cell_energy,cell_u)

test_array(cell_r_auto,cell_r)
test_array(cell_j_auto,cell_j)
test_array(cell_h_auto,cell_h)

cell_u = [ rand(l) for i in 1:L ]

cell_e = user_cell_energy(cell_u)
cell_r = user_cell_residual(cell_u)
cell_j = user_cell_jacobian(cell_u)
cell_h = cell_j

cell_r_auto = autodiff_array_gradient(user_cell_energy,cell_u)
cell_j_auto = autodiff_array_jacobian(user_cell_residual,cell_u)
cell_h_auto = autodiff_array_hessian(user_cell_energy,cell_u)

test_array(cell_r_auto,cell_r)
test_array(cell_j_auto,cell_j)
test_array(cell_h_auto,cell_h)

ids = [3,4,1,2]

function user_cell_energy_Γ(cell_u)
  cell_u_Γ = lazy_map(Reindex(cell_u),ids)
  user_cell_energy(cell_u_Γ)
end

function user_cell_residual_Γ(cell_u)
  cell_u_Γ = lazy_map(Reindex(cell_u),ids)
  user_cell_residual(cell_u_Γ)
end

function user_cell_jacobian_Γ(cell_u)
  cell_u_Γ = lazy_map(Reindex(cell_u),ids)
  user_cell_jacobian(cell_u_Γ)
end

cell_e_Γ = user_cell_energy_Γ(cell_u)
cell_r_Γ = user_cell_residual_Γ(cell_u)
cell_j_Γ = user_cell_jacobian_Γ(cell_u)
cell_h_Γ = cell_j_Γ

cell_r_Γ_auto = autodiff_array_gradient(user_cell_energy_Γ,cell_u,ids)
cell_j_Γ_auto = autodiff_array_jacobian(user_cell_residual_Γ,cell_u,ids)
cell_h_Γ_auto = autodiff_array_hessian(user_cell_energy_Γ,cell_u,ids)

test_array(cell_r_Γ_auto,cell_r_Γ)
test_array(cell_j_Γ_auto,cell_j_Γ)
test_array(cell_h_Γ_auto,cell_h_Γ)

#cell_u = [ rand(l) for i in 1:L ]

end # module