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

__precompile__(true)

module InteractiveUtils

export apropos, edit, less, code_warntype, code_llvm, code_native, methodswith, varinfo,
    versioninfo, subtypes, peakflops, @which, @edit, @less, @functionloc, @code_warntype,
    @code_typed, @code_lowered, @code_llvm, @code_native

import Base.Docs.apropos

using Base: unwrap_unionall, rewrap_unionall, isdeprecated, Bottom, show_expr_type, show_unquoted, summarysize,
    to_tuple_type, signature_type, format_bytes

using Markdown
using LinearAlgebra  # for peakflops
import Pkg

include("editless.jl")
include("codeview.jl")
include("macros.jl")

"""
    varinfo(m::Module=Main, pattern::Regex=r"")

Return a markdown table giving information about exported global variables in a module, optionally restricted
to those matching `pattern`.

The memory consumption estimate is an approximate lower bound on the size of the internal structure of the object.
"""
function varinfo(m::Module=Main, pattern::Regex=r"")
    rows =
        Any[ let value = getfield(m, v)
                 Any[string(v),
                     (value===Base || value===Main || value===Core ? "" : format_bytes(summarysize(value))),
                     summary(value)]
             end
             for v in sort!(names(m)) if isdefined(m, v) && contains(string(v), pattern) ]

    pushfirst!(rows, Any["name", "size", "summary"])

    return Markdown.MD(Any[Markdown.Table(rows, Symbol[:l, :r, :l])])
end
varinfo(pat::Regex) = varinfo(Main, pat)

"""
    versioninfo(io::IO=stdout; verbose::Bool=false, packages::Bool=false)

Print information about the version of Julia in use. The output is
controlled with boolean keyword arguments:

- `packages`: print information about installed packages
- `verbose`: print all additional information
"""
function versioninfo(io::IO=stdout; verbose::Bool=false, packages::Bool=false)
    println(io, "Julia Version $VERSION")
    if !isempty(Base.GIT_VERSION_INFO.commit_short)
        println(io, "Commit $(Base.GIT_VERSION_INFO.commit_short) ($(Base.GIT_VERSION_INFO.date_string))")
    end
    if ccall(:jl_is_debugbuild, Cint, ())!=0
        println(io, "DEBUG build")
    end
    println(io, "Platform Info:")
    println(io, "  OS: ", Sys.iswindows() ? "Windows" : Sys.isapple() ?
        "macOS" : Sys.KERNEL, " (", Sys.MACHINE, ")")

    if verbose
        lsb = ""
        if Sys.islinux()
            try lsb = readchomp(pipeline(`lsb_release -ds`, stderr=devnull)) end
        end
        if Sys.iswindows()
            try lsb = strip(read(`$(ENV["COMSPEC"]) /c ver`, String)) end
        end
        if !isempty(lsb)
            println(io, "      ", lsb)
        end
        if Sys.isunix()
            println(io, "  uname: ", readchomp(`uname -mprsv`))
        end
    end

    if verbose
        cpuio = IOBuffer() # print cpu_summary with correct alignment
        Sys.cpu_summary(cpuio)
        for (i, line) in enumerate(split(String(take!(cpuio)), "\n"))
            prefix = i == 1 ? "  CPU: " : "       "
            println(io, prefix, line)
        end
    else
        cpu = Sys.cpu_info()
        println(io, "  CPU: ", cpu[1].model)
    end

    if verbose
        println(io, "  Memory: $(Sys.total_memory()/2^30) GB ($(Sys.free_memory()/2^20) MB free)")
        try println(io, "  Uptime: $(Sys.uptime()) sec") end
        print(io, "  Load Avg: ")
        Base.print_matrix(io, Sys.loadavg()')
        println(io)
    end
    println(io, "  WORD_SIZE: ", Sys.WORD_SIZE)
    println(io, "  LIBM: ",Base.libm_name)
    println(io, "  LLVM: libLLVM-",Base.libllvm_version," (", Sys.JIT, ", ", Sys.CPU_NAME, ")")

    println(io, "Environment:")
    for (k,v) in ENV
        if contains(String(k), r"JULIA")
            println(io, "  $(k) = $(v)")
        end
    end
    if verbose
        for (k,v) in ENV
            if contains(String(k), r"PATH|FLAG|^TERM$|HOME")
                println(io, "  $(k) = $(v)")
            end
        end
    end
    if packages || verbose
        println(io, "Packages:")
        println(io, "  Package Directory: ", Pkg.dir())
        print(io, "  Package Status:")
        if isdir(Pkg.dir())
            println(io, "")
            Pkg.status(io)
        else
            println(io, " no packages installed")
        end
    end
end

function type_close_enough(@nospecialize(x), @nospecialize(t))
    x == t && return true
    # TODO: handle UnionAll properly
    return (isa(x, DataType) && isa(t, DataType) && x.name === t.name && x <: t) ||
           (isa(x, Union) && isa(t, DataType) && (type_close_enough(x.a, t) || type_close_enough(x.b, t)))
end

# `methodswith` -- shows a list of methods using the type given
"""
    methodswith(typ[, module or function]; supertypes::Bool=false])

Return an array of methods with an argument of type `typ`.

The optional second argument restricts the search to a particular module or function
(the default is all top-level modules).

If keyword `supertypes` is `true`, also return arguments with a parent type of `typ`,
excluding type `Any`.
"""
function methodswith(t::Type, f::Function, meths = Method[]; supertypes::Bool=false)
    for d in methods(f)
        if any(function (x)
                   let x = rewrap_unionall(x, d.sig)
                       (type_close_enough(x, t) ||
                        (supertypes ? (t <: x && (!isa(x,TypeVar) || x.ub != Any)) :
                         (isa(x,TypeVar) && x.ub != Any && t == x.ub)) &&
                        x != Any)
                   end
               end,
               unwrap_unionall(d.sig).parameters)
            push!(meths, d)
        end
    end
    return meths
end

function _methodswith(t::Type, m::Module, supertypes::Bool)
    meths = Method[]
    for nm in names(m)
        if isdefined(m, nm)
            f = getfield(m, nm)
            if isa(f, Function)
                methodswith(t, f, meths; supertypes = supertypes)
            end
        end
    end
    return unique(meths)
end

methodswith(t::Type, m::Module; supertypes::Bool=false) = _methodswith(t, m, supertypes)

function methodswith(t::Type; supertypes::Bool=false)
    meths = Method[]
    for mod in Base.loaded_modules_array()
        append!(meths, _methodswith(t, mod, supertypes))
    end
    return unique(meths)
end

# subtypes
function _subtypes(m::Module, x::Type, sts=Set{Any}(), visited=Set{Module}())
    push!(visited, m)
    xt = unwrap_unionall(x)
    if !isa(xt, DataType)
        return sts
    end
    xt = xt::DataType
    for s in names(m, all = true)
        if isdefined(m, s) && !isdeprecated(m, s)
            t = getfield(m, s)
            if isa(t, DataType)
                t = t::DataType
                if t.name.name === s && supertype(t).name == xt.name
                    ti = typeintersect(t, x)
                    ti != Bottom && push!(sts, ti)
                end
            elseif isa(t, UnionAll)
                t = t::UnionAll
                tt = unwrap_unionall(t)
                isa(tt, DataType) || continue
                tt = tt::DataType
                if tt.name.name === s && supertype(tt).name == xt.name
                    ti = typeintersect(t, x)
                    ti != Bottom && push!(sts, ti)
                end
            elseif isa(t, Module)
                t = t::Module
                in(t, visited) || _subtypes(t, x, sts, visited)
            end
        end
    end
    return sts
end

function _subtypes_in(mods::Array, x::Type)
    if !isabstracttype(x)
        # Fast path
        return Type[]
    end
    sts = Set{Any}()
    visited = Set{Module}()
    for m in mods
        _subtypes(m, x, sts, visited)
    end
    return sort!(collect(sts), by=string)
end

subtypes(m::Module, x::Type) = _subtypes_in([m], x)

"""
    subtypes(T::DataType)

Return a list of immediate subtypes of DataType `T`. Note that all currently loaded subtypes
are included, including those not visible in the current module.

# Examples
```jldoctest
julia> subtypes(Integer)
3-element Array{Any,1}:
 Bool
 Signed
 Unsigned
```
"""
subtypes(x::Type) = _subtypes_in(Base.loaded_modules_array(), x)

# dumptype is for displaying abstract type hierarchies,
# based on Jameson Nash's typetree.jl in https://github.com/JuliaArchive/Examples
function dumptype(io::IO, @nospecialize(x), n::Int, indent)
    print(io, x)
    n == 0 && return  # too deeply nested
    isa(x, DataType) && x.abstract && dumpsubtypes(io, x, Main, n, indent)
    nothing
end

directsubtype(a::DataType, b::DataType) = supertype(a).name === b.name
directsubtype(a::UnionAll, b::DataType) = directsubtype(a.body, b)
directsubtype(a::Union, b::DataType) = directsubtype(a.a, b) || directsubtype(a.b, b)
# Fallback to handle TypeVar's
directsubtype(a, b::DataType) = false
function dumpsubtypes(io::IO, x::DataType, m::Module, n::Int, indent)
    for s in names(m, all = true)
        if isdefined(m, s) && !isdeprecated(m, s)
            t = getfield(m, s)
            if t === x || t === m
                continue
            elseif isa(t, Module) && nameof(t) === s && parentmodule(t) === m
                # recurse into primary module bindings
                dumpsubtypes(io, x, t, n, indent)
            elseif isa(t, UnionAll) && directsubtype(t::UnionAll, x)
                dt = unwrap_unionall(t)
                println(io)
                if isa(dt, DataType) && dt.name.wrapper === t
                    # primary type binding
                    print(io, indent, "  ")
                    dumptype(io, dt, n - 1, string(indent, "  "))
                else
                    # aliases to types
                    print(io, indent, "  ", m, ".", s, "{")
                    tvar_io::IOContext = io
                    tp = t
                    while true
                        show(tvar_io, tp.var)
                        tvar_io = IOContext(tvar_io, :unionall_env => tp.var)
                        tp = tp.body
                        if isa(tp, UnionAll)
                            print(io, ", ")
                        else
                            print(io, "} = ")
                            break
                        end
                    end
                    show(tvar_io, tp)
                end
            elseif isa(t, Union) && directsubtype(t::Union, x)
                println(io)
                print(io, indent, "  ", m, ".", s, " = ", t)
            elseif isa(t, DataType) && directsubtype(t::DataType, x)
                println(io)
                if t.name.module !== m || t.name.name != s
                    # aliases to types
                    print(io, indent, "  ", m, ".", s, " = ")
                    show(io, t)
                else
                    # primary type binding
                    print(io, indent, "  ")
                    dumptype(io, t, n - 1, string(indent, "  "))
                end
            end
        end
    end
    nothing
end

const Distributed_modref = Ref{Module}()

"""
    peakflops(n::Integer=2000; parallel::Bool=false)

`peakflops` computes the peak flop rate of the computer by using double precision
[`gemm!`](@ref LinearAlgebra.BLAS.gemm!). By default, if no arguments are specified, it
multiplies a matrix of size `n x n`, where `n = 2000`. If the underlying BLAS is using
multiple threads, higher flop rates are realized. The number of BLAS threads can be set with
[`BLAS.set_num_threads(n)`](@ref).

If the keyword argument `parallel` is set to `true`, `peakflops` is run in parallel on all
the worker processors. The flop rate of the entire parallel computer is returned. When
running in parallel, only 1 BLAS thread is used. The argument `n` still refers to the size
of the problem that is solved on each processor.
"""
function peakflops(n::Integer=2000; parallel::Bool=false)
    a = fill(1.,100,100)
    t = @elapsed a2 = a*a
    a = fill(1.,n,n)
    t = @elapsed a2 = a*a
    @assert a2[1,1] == n
    if parallel
        if !isassigned(Distributed_modref)
            Distributed_modref[] = Base.require(Base, :Distributed)
        end
        Dist = Distributed_modref[]
        sum(Dist.pmap(peakflops, fill(n, Dist.nworkers())))
    else
        2*Float64(n)^3 / t
    end
end

@deprecate methodswith(typ, supertypes) methodswith(typ, supertypes = supertypes)
@deprecate whos(io::IO, m::Module, pat::Regex) show(io, varinfo(m, pat))
@deprecate whos(io::IO, m::Module)             show(io, varinfo(m))
@deprecate whos(io::IO)                        show(io, varinfo())
@deprecate whos(m::Module, pat::Regex)         varinfo(m, pat)
@deprecate whos(m::Module)                     varinfo(m)
@deprecate whos(pat::Regex)                    varinfo(pat)
@deprecate whos()                              varinfo()
@deprecate code_native(io, f, types, syntax) code_native(io, f, types, syntax = syntax)
@deprecate code_native(f, types, syntax) code_native(f, types, syntax = syntax)
# PR #21974
@deprecate versioninfo(verbose::Bool) versioninfo(verbose=verbose)
@deprecate versioninfo(io::IO, verbose::Bool) versioninfo(io, verbose=verbose)

end

max_stars_repo_name: blegat/julia

max_issues_repo_name: blegat/julia

max_forks_repo_name: blegat/julia


text:
@testset "Recombinations" begin

    compare(iter1, iter2) = length(setdiff(Set(iter1), Set(iter2))) == 0

    pop = [
        [1.0, 0.0],
        [0.0, 1.0],
        [0.0, 0.0]
    ]
    @testset "ES" begin
        @test average(pop) ≈ [1/3.,1/3.]
        @test compare( (marriage(pop) for i in 1:100), ([p1,p2] for p1 in [0,1.], p2 in [0,1.]))
    end

    @testset "ES Strategies" begin
        ss = [IsotropicStrategy{Float64}(1.0,1,0),
              IsotropicStrategy{Float64}(2.0,0,0), IsotropicStrategy{Float64}(3.0,0,0)]
        @test average(ss).σ == ss[2].σ
        @test gaussian(ss[1]).σ == ss[1].σ
        @test ss[1].σ != 1.0
        ss = [AnisotropicStrategy{Float64}(fill(1.0,3),1,1),
              AnisotropicStrategy{Float64}(fill(2.0,3),0,0),
              AnisotropicStrategy{Float64}(fill(3.0,3),0,0)]
        @test average(ss).σ == ss[2].σ
        @test gaussian(ss[1]).σ == ss[1].σ
        @test ss[1].σ != fill(1.0,3)
    end

    @testset "GA" begin
        @test identity(pop[1], pop[2]) == (pop[1], pop[2])
        @test compare((singlepoint(pop[1], pop[2]) for i in 1:100),
                      [([0.0, 1.0],[1.0, 0.0]),([1.0, 1.0],[0.0, 0.0])])
        @test compare((twopoint(pop[1], pop[2]) for i in 1:100),
                      [([0.0, 1.0],[1.0, 0.0]),([1.0, 1.0],[0.0, 0.0]),([0.0, 0.0], [1.0, 1.0])])
        @test compare((uniform(pop[1], pop[2]) for i in 1:100),
                      [([0.0, 1.0],[1.0, 0.0]),([1.0, 1.0],[0.0, 0.0]),([0.0, 0.0], [1.0, 1.0]),([1.0, 0.0], [0.0, 1.0])])
        @test compare((discrete(pop[1], pop[2]) for i in 1:100),
                      ([p1,p2],[p3,p4]) for p1 in [0,1.], p2 in [0,1.], p3 in [0,1.], p4 in [0,1.])

        @test waverage([1.0,0.5])(pop[1], pop[2]) == ([1.0, 2.0], [1.0, 2.0])

        @test map(sum, HX(pop[1], pop[2])) |> collect ≈ [1.0, 1.0]
        @test map(abs∘first∘diff, LX()(pop[1], pop[2])) |> collect ≈ [1.0, 1.0]
        @test mapslices(diff, hcat(MILX(1., 0., 1.)(Real[0, 0.0], Real[1, 0.0])...), dims=2) |> vec ≈ [1.0, 0.0]

        @test PMX(pop[1], pop[2]) == ([0.0, 1.0], [1.0, 0.0])
        @test sum(PMX(collect(1:9), collect(9:-1:1))) == fill(10,9)
        @test sum(OX1(collect(1:9), collect(9:-1:1))) == fill(10,9)
        @test sum(OX2(collect(1:9), collect(9:-1:1))) == fill(10,9)
        @test sum(CX(collect(1:9), collect(9:-1:1))) == fill(10,9)
        @test sum(POS(collect(1:9), collect(9:-1:1))) == fill(10,9)
        @test CX([8,4,7,3,6,2,5,1,9,0],collect(0:9)) == ([8,1,2,3,4,5,6,7,9,0],[0,4,7,3,6,2,5,1,8,9])
    end

    @testset "DE" begin
        v1 = collect(1:10)
        v2 = fill(0,10)

        xvr = uniformbin(0.0)
        @test first( xvr(v1, v2) )== v1
        xvr = uniformbin(1.0)
        @test last( xvr(v1, v2) )== v1
        xvr = uniformbin(0.5)
        m1, m2 = xvr(v1, v2)
        @test m1.+m2 == v1

        xvr = exponential(0.0)
        @test sum(first( xvr(v1, v2) )) <= 10
        xvr = exponential(1.0)
        @test first( xvr(v1, v2) )== v1
        xvr = exponential(0.05)
        @test sum(first( xvr(v1, v2) )) <= 19
    end

end

ext: jl


text:
@safetestset hilbert_decode_zero = "hilbert decode zero" begin
  using BijectiveHilbert
  table2 = [
    15 12 11 10
    14 13  8  9
    1  2  7  6
    0  3  4  5
  ]
  table2 = table2[end:-1:1, :]

  table2ind = CartesianIndices(table2)
  gg = Simple2D(Int)
  X = zeros(Int, 2)
  for z in 0:15
    BijectiveHilbert.decode_hilbert_zero!(gg, X, z)
    found = table2ind[table2 .== z][1]
    @test found[2] - 1 == X[1]
    @test found[1] - 1 == X[2]
  end
end

@safetestset simple2d_own_inverse = "simple2d is its own inverse" begin
  using BijectiveHilbert: Simple2D
  using ..HilbertTestSuite: check_own_inverse
  gg = Simple2D(Int)
  for b in 2:7
    @test check_own_inverse(gg, b, 2)
  end
end

@safetestset simple2d_complete_set = "simple2d is a complete set" begin
  using BijectiveHilbert: Simple2D
  using ..HilbertTestSuite: check_complete_set
  gg = Simple2D(Int)
  for b in 2:7
    @test check_complete_set(gg, b, 2)
  end
end

@safetestset simple2d_type_interactions = "Simple2D type interactions" begin
    using BijectiveHilbert
    using UnitTestDesign
    using Random
    rng = Random.MersenneTwister(9790323)
    for retrial in 1:5
        AxisTypes = shuffle(rng, [Int8, Int, UInt, Int128, UInt8, UInt128])
        IndexTypes = shuffle(rng, [Int8, UInt8, Int, UInt, Int128, UInt128])
        Count= shuffle(rng, [0, 1])
        Dims = shuffle(rng, [2, 3, 4])
        Bits = shuffle(rng, [2, 3, 4, 5])
        test_set = all_pairs(
            AxisTypes, IndexTypes, Count, Dims, Bits;
        )
        for (A, I, C, D, B) in test_set
            gg = Simple2D(I)
            if B * D > log2(typemax(I))
                continue
            end
            # Add this because these values aren't in Hilbert order because
            # It's an asymmetrical space.
            if B * D > log2(typemax(A))
              continue
            end
          last = (one(I) << (B * D)) - one(I) + I(C)
            mid = one(I) << (B * D - 1)
            few = 5
            X = zeros(A, D)
            hlarr = vcat(C:min(mid, few), max(mid + 1, last - few):last)
            for hl in hlarr
                hli = I(hl)
                if C == 0
                    decode_hilbert_zero!(gg, X, hli)
                    hl2 = encode_hilbert_zero(gg, X)
                    if hl2 != hli
                        @show A, I, C, D, B, X
                        @test hl2 == hli
                    end
                    @test typeof(hl2) == typeof(hli)
                else
                    decode_hilbert!(gg, X, hli)
                    hl3 = encode_hilbert(gg, X)
                    @test hl3 == hli
                    @test typeof(hl3) == typeof(hli)
                end
            end
        end
    end
end


text:
# Note: figs4-7 involve memory-intensive calculations
# All simulations were done in 3 machines:
# AMD Ryzen 9 5950X 16-Core Processor, 64Gb
# Combining distributed parallelization and MKL multithreading
using MBSinBLG
fig2run()
fig3run()
fig4run()
fig5run()
fig6run()
fig7run()

max_stars_repo_name: fernandopenaranda/MBSinBLG

max_issues_repo_name: fernandopenaranda/MBSinBLG

max_forks_repo_name: fernandopenaranda/MBSinBLG


text:
# TODO may want to test this but it's very slow on azure
@otestset "Examples" begin
    @test_all_no_arg_functions Examples
end # testset Examples


text:
#==
    Note on implementation of factorisition types:
    The general strategy is to create a named factorization type which wraps other factorizations.
    This allows us to define `getproperty` methods which return `NamedDimsArray`s around the
    corresponding properties of the parent factorization types.

    NOTE: Maybe we could just have a `Named` type that works with arrays and factorizations
    kind of like `Adjoint`?
==#

struct NamedFactorization{L, T, F<:Factorization{T}} <: Factorization{T}
    fact::F
end

function NamedFactorization{L}(fact::F) where {L, T, F <: Factorization{T}}
    return NamedFactorization{L,T,F}(fact)
end

# Need parent to explicitly call `getfield` to avoid infinitely calling `getproperty`
Base.parent(named::NamedFactorization) = getfield(named, :fact)
Base.size(named::NamedFactorization) = size(parent(named))
Base.propertynames(named::NamedFactorization; kwargs...) = propertynames(parent(named))

# Factorization type specific initial iterate calls
Base.iterate(named::NamedFactorization{L, T, <:LU}) where {L, T} = (named.L, Val(:U))
Base.iterate(named::NamedFactorization{L, T, <:LQ}) where {L, T} = (named.L, Val(:Q))
Base.iterate(named::NamedFactorization{L, T, <:SVD}) where {L, T} = (named.U, Val(:S))
function Base.iterate(
    named::NamedFactorization{L, T, <:Union{QR, LinearAlgebra.QRCompactWY, QRPivoted}}
) where {L, T}
    return (named.Q, Val(:R))
end

# Generic iterate follow up calls
function Base.iterate(named::NamedFactorization, st::Val{D}) where D
    r = iterate(parent(named), st)
    r === nothing && return nothing
    return (getproperty(named, D), last(r))
end

# Convenience constructors
for func in (:lu, :lu!, :lq, :lq!, :svd, :svd!, :qr, :qr!)
    @eval begin
        function LinearAlgebra.$func(nda::NamedDimsArray{L, T}, args...; kwargs...) where {L, T}
            return NamedFactorization{L}($func(parent(nda), args...; kwargs...))
        end
    end
end

# `getproperty` wrappers
## LU
function Base.getproperty(fact::NamedFactorization{L, T, <:LU}, d::Symbol) where {L, T}
    inner = getproperty(parent(fact), d)
    n1, n2 = L
    if d === :L
        return NamedDimsArray{(n1, :_)}(inner)
    elseif d === :U
        return NamedDimsArray{(:_, n2)}(inner)
    elseif d === :P
        perm_matrix_labels = (first(L), first(L))
        return NamedDimsArray{perm_matrix_labels}(inner)
    elseif d === :p
        perm_vector_labels = (first(L),)
        return NamedDimsArray{perm_vector_labels}(inner)
    else
        return inner
    end
end

## LQ
function Base.getproperty(fact::NamedFactorization{L, T, <:LQ}, d::Symbol) where {L, T}
    inner = getproperty(parent(fact), d)
    n1, n2 = L
    if d === :L
        return NamedDimsArray{(n1, :_)}(inner)
    elseif d === :Q
        return NamedDimsArray{(:_, n2)}(inner)
    else
        return inner
    end
end

## svd
function Base.getproperty(fact::NamedFactorization{L, T, <:SVD}, d::Symbol) where {L, T}
    inner = getproperty(parent(fact), d)
    n1, n2 = L
    # Naming based off the SVD visualization on wikipedia
    # # https://en.wikipedia.org/wiki/File:Singular_value_decomposition_visualisation.svg
    if d === :U
        return NamedDimsArray{(n1, :_)}(inner)
    elseif d === :V
        return NamedDimsArray{(n2, :_)}(inner)
    elseif d === :Vt
        return NamedDimsArray{(:_, n2)}(inner)
    else # :S
        return inner
    end
end

## qr
function Base.getproperty(
    fact::NamedFactorization{L, T, F},
    d::Symbol
) where {L, T, F<:Union{QR, LinearAlgebra.QRCompactWY, QRPivoted}}
    inner = getproperty(parent(fact), d)
    n1, n2 = L
    if d === :Q
        return NamedDimsArray{(n1,:_)}(inner)
    elseif d === :R
        return NamedDimsArray{(:_, n2)}(inner)
    elseif F <: QRPivoted && d === :P
        return NamedDimsArray{(n2, n2)}(inner)
    elseif F <: QRPivoted && d === :p
        return NamedDimsArray{(n2,)}(inner)
    else
        return inner
    end
end


text:
using JLD
using GaussianMixtures

setmem(0.1)
include("data.jl")
include("bayes.jl")
include("train.jl")
include("scikitlearn.jl")


text:
export AnEnKF

"""
    EnKF( np )

Ensemble Kalman Filters
"""
struct AnEnKF <: MonteCarloMethod

    np::Int64

end

"""
    forecast( da, yo, mc; progress=true)

- da: DataAssimilation
- yo: Observations
- mc: Monte Carlo method
- progress: display the progress bar.

Apply stochastic and sequential data assimilation technics using
model forecasting or analog forecasting.
"""
function forecast(da::DataAssimilation, yo::TimeSeries, mc::AnEnKF; progress = true)

    # dimensions
    nt = yo.nt        # number of observations
    np = mc.np        # number of particles
    nv = yo.nv        # number of variables (dimensions of problem)

    # initialization
    x̂ = TimeSeries(nt, nv)

    m_xa_part = [zeros(Float64, (nv, np)) for i = 1:nt]
    part = [zeros(Float64, (nv, np)) for i = 1:nt]
    pf = [zeros(Float64, (nv, nv)) for i = 1:nt]
    xf = zeros(Float64, (nv, np))
    xf_mean = similar(xf)
    ef = similar(xf)
    Ks = zeros(Float64, (3, 3))

    if progress
        p = Progress(nt)
    end

    for k = 1:nt

        if progress
            next!(p)
        end

        # update step (compute forecasts)
        if k == 1
            xf .= rand(MvNormal(da.xb, da.B), np)
        else
            xf, xf_mean = da.m(part[k-1])
            m_xa_part[k] .= xf_mean
        end

        part[k] .= xf

        ef .= xf * (Matrix(I, np, np) .- 1 / np)
        pf[k] .= (ef * ef') ./ (np - 1)

        # analysis step (correct forecasts with observations)
        ivar_obs = findall(.!isnan.(yo.u[k]))
        n = length(ivar_obs)

        if n > 0
            μ = zeros(Float64, n)
            σ = da.R[ivar_obs, ivar_obs]
            eps = rand(MvNormal(μ, σ), np)
            yf = da.H[ivar_obs, :] * xf
            Σ = (da.H[ivar_obs, :] * pf[k]) * da.H[ivar_obs, :]'
            Σ .+= da.R[ivar_obs, ivar_obs]
            invΣ = inv(Σ)
            K = (pf[k] * da.H[ivar_obs, :]') * invΣ
            d = yo.u[k][ivar_obs] .- yf .+ eps
            part[k] .= xf .+ K * d
            # compute likelihood
            innov_ll = mean(yo.u[k][ivar_obs] .- yf, dims = 2)
        else
            part[k] .= xf
        end

        x̂.u[k] .= vec(sum(part[k] ./ np, dims = 2))

    end

    x̂

end


text:
@testset "restructure" begin
    x = rand(Float32, 2, 2)
    y = rand(4)
    yr = ArrayInterface.restructure(x, y)
    @test yr isa Matrix{Float64}
    @test size(yr) == (2,2)
    @test vec(yr) == vec(y)

    @testset "views" begin
        x = @view rand(4)[1:2]
        y = rand(2)
        yr = ArrayInterface.restructure(x, y)
        @test yr isa Vector{Float64}
        @test size(yr) == (2,)
        @test yr == y

        x = @view rand(4,4)[1:2,1:2]
        y = rand(2,2)
        yr = ArrayInterface.restructure(x, y)
        @test yr isa Matrix{Float64}
        @test size(yr) == (2,2)
        @test yr == y

        x = @view rand(4,4)[1]
        y = @view rand(2,2)[1]
        yr = ArrayInterface.restructure(x, y)
        @test yr isa Array{Float64,0}
        @test size(yr) == ()
        @test yr == y
    end
end


text:
const ModelSolution{T,NPT_FULL,NP_SAVED} = Union{ODECompositeSolution{T,NPT_FULL,<:AbstractArray{<:AbstractArray{T,NP_SAVED}}},ODESolution{T,NPT_FULL,<:AbstractArray{<:AbstractArray{T,NP_SAVED}}}}#,RODESolution{T,NPT_FULL,<:AbstractArray{<:AbstractArray{T,NP_SAVED}}}}

# FIXME not real dispatch, since it's just an alias
@inline population(A::AbstractArray{T,N}, i) where {T,N} = view_slice_last(A, i)
function population_coordinates(coordinates::AbstractArray{<:CartesianIndex,N}, P) where N
    cat([[CartesianIndex(coord,i) for coord in coordinates] for i in 1:P]...; dims=N+1)
end
population_repeat(arr::AbstractArray{T,N}, P) where {T,N} = repeat(arr, outer=([1 for _ in 1:N]..., P))
"""
    population_timepoint(solution, pop_dx, time_dx)

Return spatial frame for a given population `pop_dx` and time `time_dx`.
"""
@generated function population_timepoint(solution::ModelSolution{T,NPT_FULL,NP_SAVED}, pop_dx::Int, time_dx::Int) where {T,NPT_FULL,NP_SAVED}
    N = NP_SAVED - 1 # N + pops
    colons = [:(:) for i in 1:N]
    :(solution(solution.t[time_dx])[$(colons...), pop_dx])
end

# Going forward, make macro to construct actions and interactions between
# desired numbers of populations

# Essentially: Slow, bad type inference at construction level, but good, strong type
# inference at the solver-level.
# e.g.: use length of vector input to choose which type to construct, but then when the type
# is given to the solver, it's set and so all the inference is set for the hard part.

type_name(T::DataType) = T.name.wrapper
type_name(T::UnionAll) = type_name(T.body)
function (ap::Type{<:AbstractParameter})(pop::POP) where {POP <: AbstractPopulationP}
	type_name(POP)((ap(param) for param in getfield.(Ref(pop), fieldnames(POP)))...)
end

abstract type AbstractPOneD{NP,AP} <: AbstractPopulationP{NP,AP} end
abstract type AbstractPTwoD{NP,AP} <: AbstractPopulationP{NP,AP} end
abstract type AbstractPopulationActionsParameters{NP,AP} <: AbstractPOneD{NP,AP} end
abstract type AbstractPopulationInteractionsParameters{NP,AP} <: AbstractPTwoD{NP,AP} end
abstract type AbstractPopulationActions{NP,AP} <: AbstractPOneD{NP,AP} end
abstract type AbstractPopulationInteractions{NP,AP} <: AbstractPTwoD{NP,AP} end

const PopActParam{P,T} = Union{Array{<:AbstractPopulationActionsParameters,1},AbstractPopulationActionsParameters{P,<:T}}
const PopInteractParam{P,T} = Union{Array{<:AbstractPopulationInteractionsParameters,1},AbstractPopulationInteractionsParameters{P,T}}
const PopAct{P,T} = Union{Array{<:AbstractPopulationActions,1},AbstractPopulationActions{P,T}}
const PopInteract{P,T} = Union{Array{<:AbstractPopulationInteractions,1},AbstractPopulationInteractions{P,T}}

function array(p1d::AbstractPOneD{2})
    [p1d.p1, p1d.p2]
end

function array(p2d::AbstractPTwoD{2})
    [p2d.p11 p2d.p12;
     p2d.p21 p2d.p22]
end

struct PopulationActionsParameters2{P,P1<:P,P2<:P} <: AbstractPopulationActionsParameters{2,P}
    p1::P1
    p2::P2
    PopulationActionsParameters2(p1::P1,p2::P2) where {P1,P2} = new{typejoin(P1,P2),P1,P2}(p1,p2)
    PopulationActionsParameters2{P,P1,P2}(p1::P1,p2::P2) where {P,P1<:P,P2<:P} = new(p1,p2)
end
(pp::PopulationActionsParameters2)(args...) = PopulationActions2(pp.p1(args...), pp.p2(args...))

struct PopulationActions2{P,P1<:P,P2<:P} <: AbstractPopulationActions{2,P}
    p1::P1
    p2::P2
    PopulationActions2(p1::P1,p2::P2) where {P1,P2} = new{typejoin(P1,P2),P1,P2}(p1,p2)
end
function (pop_actions::PopulationActions2)(inplace::AbstractArray{T,N}, source::AbstractArray{T,N}, t) where {T,N}
    pop_actions.p1(population(inplace, 1), population(source, 1), t)
    pop_actions.p2(population(inplace, 2), population(source, 2), t)
end

struct PopulationInteractionsParameters2{P,P11<:P,P12<:P,P21<:P,P22<:P} <: AbstractPopulationInteractionsParameters{2,P}
    p11::P11
	p21::P21
    p12::P12
    p22::P22
    PopulationInteractionsParameters2(p11::P11,p21::P21,p12::P12,p22::P22) where {P11,P12,P21,P22} = new{typejoin(P11,P12,P21,P22),P11,P12,P21,P22}(p11,p21,p12,p22)
end
struct PopulationInteractions2{P,P11<:P,P12<:P,P21<:P,P22<:P} <: AbstractPopulationInteractions{2,P}
    p11::P11
    p21::P21
    p12::P12
    p22::P22
    PopulationInteractions2(p11::P11,p21::P21,p12::P12,p22::P22) where {P11,P12,P21,P22} = new{typejoin(P11,P12,P21,P22),P11,P12,P21,P22}(p11,p21,p12,p22)
end
(pp::PopulationInteractionsParameters2)(args...) = PopulationInteractions2(pp.p11(args...), pp.p21(args...),
                                                                           pp.p12(args...), pp.p22(args...))
function (pop_interactions::PopulationInteractions2)(inplace::ARR1, source::ARR2, t) where {T,N,ARR1<:AbstractArray{T,N},ARR2<:AbstractArray{T,N}}
	pop_interactions.p11(population(inplace, 1), population(source, 1), t)
    pop_interactions.p12(population(inplace, 1), population(source, 2), t)
    pop_interactions.p21(population(inplace, 2), population(source, 1), t)
    pop_interactions.p22(population(inplace, 2), population(source, 2), t)
end
export PopulationInteractions2

function (arr::Array{<:AbstractPopulationP})(args...)
	[p(args...) for p in arr]
end

PopulationActionsParameters(p1::AbstractParameter, p2::AbstractParameter) = PopulationActionsParameters2(p1, p2)
PopulationActionsParameters(p1::AbstractAction, p2::AbstractAction) = PopulationActions2(p1, p2)
PopulationInteractionsParameters(p11::AbstractParameter,p21::AbstractParameter,p12::AbstractParameter,p22::AbstractParameter) = PopulationInteractionsParameters2(p11,p21,p12,p22)
PopulationInteractionsParameters(p11::AbstractAction,p21::AbstractAction,p12::AbstractAction,p22::AbstractAction) = PopulationInteractions2(p11,p21,p12,p22)

const MissingAP = Union{Missing, <:AbstractParameter}
const MissingAA = Union{Missing, <:AbstractAction}

PopulationActionsParameters(p1::MissingAP, p2::MissingAP) = missing
PopulationInteractionsParameters(p11::MissingAP, p21::MissingAP, p12::MissingAP, p22::MissingAP) = missing

function PopulationActionsParameters(p1::AbstractAction, p2::Missing)
    missing
end
function PopulationActionsParameters(p1::Missing, p2::AbstractAction)
    missing
end
function PopulationActionsParameters(p1::Missing, p2::Missing)
    missing
end
PopulationInteractionsParameters(p11::MissingAA, p21::MissingAA, p12::MissingAA, p22::MissingAA) = missing

function pops(t::Type{T}; kwargs...) where T
    syms = keys(kwargs)
    args = zip(values(kwargs)...)
    new_kwargs = (zip(syms, arg) for arg in args)
    params = map((single_kwargs) -> t(;single_kwargs...), new_kwargs)
	if ndims(params) == 1
		return PopulationActionsParameters(params...)
	elseif ndims(params) == 2
		return PopulationInteractionsParameters(params...)
	else
		error("Weirdly shaped parameters.")
	end
end

function pops(ps::Vector)# where {P<:AbstractParameter}
    PopulationActionsParameters(ps...)
end
function pops(ps::Matrix)# where {P<:AbstractParameter}
    PopulationInteractionsParameters(ps...)
end

function pops(ps::Vector{A}) where {A <: AbstractPopulationInteractionsParameters}
    PopulationInteractionsParameters([CompositeParameter([pop_params...]) for pop_params in zip(array.(ps)...)]...)
end


text:
struct RestrictLowDistance
    max_width::Int
end

function Diderot.process(
    restrict::RestrictLowDistance, layer::Layer{S,D,V}
) where {S,D,V}
    if length(layer) <= restrict.max_width
        return layer
    end

    candidates = collect(layer)
    sort!(candidates, by=tup -> tup.second.distance, rev=true)
    return Layer{S,D,V}(Dict(candidates[1:restrict.max_width]), false)
end


text:
using LinearAlgebra
### Solver Options ###
opts = TrajectoryOptimization.SolverOptions()
opts.verbose = true
######################

### Set up model, objective, solver ###
# Model
dt = 0.1
traj_folder = joinpath(dirname(pathof(TrajectoryOptimization)),"..")
urdf_folder = joinpath(traj_folder, "dynamics/urdf")
urdf_cartpole = joinpath(urdf_folder, "cartpole.urdf")

model_urdf = Model(urdf_cartpole,[1.;0.])
n = model.n
m = model.m
# Objective
Q = 0.01*Diagonal(I,n)
Qf = 1000.0*Diagonal(I,n)
R = 0.01*Diagonal(I,m)

x0 = [0.;pi;0.;0.]
xf = [0.;0.;0.;0.]

tf = 5.0

obj_uncon = LQRObjective(Q, R, Qf, tf, x0, xf)

# -Constraints
u_min = -20
u_max = 1000
x_min = [-1000; -1000; -1000; -1000]
x_max = [1000; 1000; 1000; 1000]

obj_con = TrajectoryOptimization.ConstrainedObjective(obj_uncon, u_min=u_min, u_max=u_max, x_min=x_min, x_max=x_max) # constrained objective

# Solver (foh & zoh)
# solver_foh = Solver(model, obj_con, integration=:rk3_foh, dt=dt, opts=opts)
solver_zoh = Solver(model, obj_uncon, integration=:rk3, dt=dt, opts=opts)

# -Initial control and state trajectories
U = ones(solver_zoh.model.m,solver_zoh.N)
X_interp = line_trajectory(solver_zoh)
X_rollout = TrajectoryOptimization.rollout(solver,U)
# X_interp = ones(solver_zoh.model.n,solver.N)
#######################################

### Solve ###
# @time sol_foh, = TrajectoryOptimization.solve(solver_foh,X_interp,U)
solver_zoh.opts.verbose = false
@time sol_zoh, = TrajectoryOptimization.solve(solver_zoh,U)
solve_dircol(solver_zoh,X_rollout,U,method=:hermite_simpson)
#############

# ### Results ###
# if opts.verbose
#     # println("Final state (foh): $(sol_foh.X[:,end])")
#     println("Final state (zoh): $(sol_zoh.X[:,end])")
#
#     # println("Termination index\n foh: $(sol_foh.termination_index)\n zoh: $(sol_foh.termination_index)")
#
#     # println("Final cost (foh): $(sol_foh.cost[sol_foh.termination_index])")
#     println("Final cost (zoh): $(sol_zoh.cost[sol_zoh.termination_index])")
#
#     # plot((sol_foh.cost[1:sol_foh.termination_index]))
#     plot!((sol_zoh.cost[1:sol_zoh.termination_index]))
#
#     # plot(sol_foh.U')
#     plot!(sol_zoh.U')
#
#     # plot(sol_foh.X[1:2,:]')
#     plot!(sol_zoh.X[1:2,:]')
# end
###############


text:
# ---
# title: 206. Reverse Linked List
# id: problem206
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Linked List
# link: <https://leetcode.com/problems/reverse-linked-list/description/>
# hidden: true
# ---
#
# Reverse a singly linked list.
#
# **Example:**
#
#
#
#     Input: 1->2->3->4->5->NULL
#     Output: 5->4->3->2->1->NULL
#
#
# **Follow up:**
#
# A linked list can be reversed either iteratively or recursively. Could you
# implement both?
#
#
## @lc code=start
using LeetCode

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


text:
using Test, Waveforms

x = 0.4
y = 3π/7
θ = 0.05

@testset "square waves" begin
    @testset "squarewave" begin
        @test squarewave(x) == 1.0
        @test squarewave(y) == 1.0
        @test squarewave(x, θ) == -1.0
        @test squarewave(y, θ) == -1.0
        @test_throws DomainError squarewave(x, -1.0)
        @test_throws DomainError squarewave(x, 2.0)
    end
    @testset "squarewave1" begin
        @test squarewave1(x) == 1.0
        @test squarewave1(y) == 1.0
        @test squarewave1(x, θ) == -1.0
        @test squarewave1(y, θ) == -1.0
        @test_throws DomainError squarewave1(x, -1.0)
        @test_throws DomainError squarewave1(x, 2.0)
    end
end

@testset "triangle waves" begin
    @testset "trianglewave" begin
        @test trianglewave(x) == 0.2546479089470326
        @test trianglewave(y) == 0.8571428571428572
    end
    @testset "trianglewave1" begin
        @test trianglewave1(x) == 0.3999999999999999
        @test trianglewave1(y) == 0.6144125938460689
    end
end

@testset "sawtooth waves" begin
    @testset "sawtoothwave" begin
        @test sawtoothwave(x) == 0.12732395447351627
        @test sawtoothwave(y) == 0.42857142857142855
    end
    @testset "sawtoothwave1" begin
        @test sawtoothwave1(x) == 0.8
        @test sawtoothwave1(y) == 0.6927937030769655
    end
end


text:
using Agents, Random, AgentsPlots, Plots, ProgressMeter

mutable struct FMP_Agent <: AbstractAgent
    id::Int
    pos::NTuple{2, Float64}
    vel::NTuple{2, Float64}
    tau::NTuple{2, Float64}
    color::String
    type::Symbol
    radius::Float64
    SSdims::NTuple{2, Float64}  # include this for plotting
    # NOTE: if you change anything in this you need to restart the REPL
    # (I think it is the precompilation step)
end

"""
Initialization function for FMP simulation. Contains all model parameters.
"""
function FMP_Model(simtype;
                   rho = 7.5e6,
                   rho_obstacle = 7.5e6,
                   step_inc = 2,
                   dt = 0.01,
                   num_agents = 10,
                   SS_dims = (1, 1),  # x,y should be equal for proper plot scaling
                   num_steps = 300,
                   terminal_max_dis = 0.01,
                   c1 = 10,
                   c2 = 10,
                   vmax = 0.1,
                   d = 0.02, # distance from centroid to centroid
                   r = (3*vmax^2/(2*rho))^(1/3)+d,
                   obstacle_list = [],
                  )

    # define AgentBasedModel (ABM)
    properties = Dict(:rho=>rho,
                      :rho_obstacle=>rho_obstacle,
                      :step_inc=>step_inc,
                      :r=>r,
                      :d=>d,
                      :dt=>dt,
                      :num_agents=>num_agents,
                      :num_steps=>num_steps,
                      :terminal_max_dis=>terminal_max_dis,
                      :c1=>c1,
                      :c2=>c2,
                      :vmax=>vmax,
                      :obstacle_list=>obstacle_list,
                     )

    space2d = ContinuousSpace(SS_dims; periodic=true)
    model = ABM(FMP_Agent, space2d, properties=properties)
    AgentPositionInit(model, num_agents; type=simtype)

    # append obstacles into obstacle_list
    for agent in allagents(model)
        if agent.type == :O
            append!(model.obstacle_list, agent.id)
        end
    end

    return model

end

"""
Simulation wrapper for FMP simulations.

Initializes model based on "type" parameter which dictates type of simulation to perform. Different simulation descriptions can be found in `/multiagent/simulation_init.jl`.

Possible inputs include:
- `circle` have agents move/swap places around perimeter of a circle
- `circle_object` have agents move/swap places around perimeter of a circle with an object in the middle
- `line` have agents move left to right in vertical line
- `centered_line_object` have agents remain stationary in vertical line and an object move through them
- `moving_line` have agents move left to right in vertical line past an object
- `random` have agents start in random positions with random velocities

Next it loops through the number of simulation steps (specified in model params) and create simulation display using `plotabm()`.

Finished by saving at the location specified by `outputpath` variable.
"""
function FMP_Simulation(simtype::String; outputpath = "output/simresult.gif")
    gr()
    cd(@__DIR__)

    # init model
    model = FMP_Model(simtype)
    agent_step!(agent, model) = move_agent!(agent, model, model.dt)

    # init state space
    e = model.space.extent
    step_range = 1:model.step_inc:model.num_steps

    mean_norms = Array{Float64}(undef,1,)

    # setup progress meter counter
    p = Progress(round(Int,model.num_steps/model.step_inc))
    anim = @animate for i in step_range

        # step model including plot stuff
        FMP(model)
        p1 = AgentsPlots.plotabm(
            model,
            as = PlotABM_RadiusUtil,
            ac = PlotABM_ColorUtil,
            am = PlotABM_ShapeUtil,
            #showaxis = false,
            grid = false,
            xlims = (0, e[1]),
            ylims = (0, e[2]),
            aspect_ratio=:equal,
            scheduler = PlotABM_Scheduler,
        )
        title!(p1, "FMP Simulation (step $(i))")

        # step model and progress counter
        step!(model, agent_step!, model.step_inc)
        next!(p)
    end
    gif(anim, outputpath, fps = 100)

end

max_stars_repo_name: rmcsqrd/antfarm

max_issues_repo_name: rmcsqrd/antfarm

max_forks_repo_name: rmcsqrd/antfarm


text:
"""
Service Income

any retruns  (in the form) of service conducted

online , or offline:

1.paid Classes
2.Online Teaching
3.Coaching

"""

mutable struct serviceIncome :> Income end

max_stars_repo_name: adamwillisXanax/accountingJulia

max_issues_repo_name: adamwillisXanax/accountingJulia

max_forks_repo_name: adamwillisXanax/accountingJulia


text:
# ---
# title: 229. Majority Element II
# id: problem229
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/majority-element-ii/description/>
# hidden: true
# ---
#
# Given an integer array of size `n`, find all elements that appear more than `⌊
# n/3 ⌋` times.
#
# **Follow-up:** Could you solve the problem in linear time and in O(1) space?
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [3,2,3]
#     Output: [3]
#
#
# **Example 2:**
#
#
#
#     Input: nums = [1]
#     Output: [1]
#
#
# **Example 3:**
#
#
#
#     Input: nums = [1,2]
#     Output: [1,2]
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 5 * 104`
#   * `-109 <= nums[i] <= 109`
#
#
## @lc code=start
using LeetCode

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


text:
using SimulatedNeuralMoments, Flux, MCMCChains, StatsPlots, DelimitedFiles, Optim
using BSON:@save
using BSON:@load
using DelimitedFiles
include("JDlib.jl")
lb, ub = PriorSupport()
# fill in the structure that defines the model
model = SNMmodel("SP500 estimation", lb, ub, InSupport, Prior, PriorDraw, auxstat)
# train the net, and save it and the transformation info
nnmodel, nninfo = MakeNeuralMoments(model)
@save "neuralmodel.bson" nnmodel nninfo  # use this line to save the trained neural net
@load "neuralmodel.bson" nnmodel nninfo # use this to load a trained net
data = readdlm("sp500.txt")
m = NeuralMoments(auxstat(data), model, nnmodel, nninfo)
@show m
m = min.(max.(lb,m),ub)
@show m
chain, P, tuning = MCMC(m, 10000, model, nnmodel, nninfo, verbosity=true, do_cue = true, tuningloops=2, nthreads=10, burnin=40, tuning=0.75)
# save visualize results
writedlm("chain", chain)
writedlm("P", P)
writedlm("tuning", tuning)

max_stars_repo_name: mcreel/ILSNM

max_issues_repo_name: mcreel/ILSNM

max_forks_repo_name: mcreel/ILSNM


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Private functions.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

################################################################################
#                                Preprocessing
################################################################################

# Those functions apply a preprocessing to the data that will be printed
# depending on its type.

function _preprocess_vec_or_mat(
    data::AbstractVecOrMat,
    header::Union{Nothing, AbstractVector, Tuple}
)
    if header === nothing
        pheader = (["Col. " * string(i) for i = 1:size(data,2)],)
    elseif header isa AbstractVector
        pheader = (header,)
    else
        pheader = header
    end

    return data, pheader
end

function _preprocess_dict(
    dict::AbstractDict{K, V};
    sortkeys::Bool = false
) where {K, V}
    pheader = (
        ["Keys", "Values"],
        [compact_type_str(K), compact_type_str(V)]
    )

    k = collect(keys(dict))
    v = collect(values(dict))

    if sortkeys
        ind = sortperm(collect(keys(dict)))
        vk  = k[ind]
        vv  = v[ind]
    else
        vk = k
        vv = v
    end

    pdata = hcat(vk, vv)

    return pdata, pheader
end

function _preprocess_Tables_column(
    data::Any,
    header::Union{Nothing, AbstractVector, Tuple}
)
    # Access the table using the columns.
    table = Tables.columns(data)

    # Get the column names.
    names = collect(Symbol, Tables.columnnames(table))

    # Compute the table size and get the column types.
    size_j::Int = length(names)
    size_i::Int = Tables.rowcount(table)

    pdata = ColumnTable(data, table, names, (size_i, size_j))

    # For the header, we have the following priority:
    #
    #     1. If the user passed a vector `header`, then use it.
    #     2. Otherwise, check if the table defines a schema to create the
    #        header.
    #     3. If the table does not have a schema, then build a default header
    #        based on the column name and type.
    if header === nothing
        sch = Tables.schema(data)

        if sch !== nothing
            types::Vector{String} = compact_type_str.([sch.types...])
            pheader = (names, types)
        else
            pheader = (pdata.column_names,)
        end
    elseif header isa AbstractVector
        pheader = (header,)
    else
        pheader = header
    end

    return pdata, pheader
end

function _preprocess_Tables_row(
    data::Any,
    header::Union{Nothing, AbstractVector, Tuple}
)
    # Access the table using the rows.
    table = Tables.rows(data)

    # We need to fetch the first row to get information about the columns.
    row₁, ~ = iterate(table, 1)

    # Get the column names.
    names = collect(Symbol, Tables.columnnames(row₁))

    # Compute the table size.
    size_i::Int = length(table)
    size_j::Int = length(names)

    pdata = RowTable(data, table, names, (size_i, size_j))

    # For the header, we have the following priority:
    #
    #     1. If the user passed a vector `header`, then use it.
    #     2. Otherwise, check if the table defines a schema to create the
    #        header.
    #     3. If the table does not have a schema, then build a default header
    #        based on the column name and type.
    if header === nothing
        sch = Tables.schema(data)

        if sch !== nothing
            types::Vector{String} = compact_type_str.([sch.types...])
            pheader = (names, types)
        else
            pheader = (pdata.column_names,)
        end
    elseif header isa AbstractVector
        pheader = (header,)
    else
        pheader = header
    end

    return pdata, pheader
end

################################################################################
#                              Print information
################################################################################

# This function creates the structure that holds the global print information.
function _print_info(
    data::Any;
    alignment::Union{Symbol, Vector{Symbol}} = :r,
    cell_alignment::Union{
        Nothing,
        Dict{Tuple{Int, Int}, Symbol},
        Function,
        Tuple
    } = nothing,
    cell_first_line_only::Bool = false,
    compact_printing::Bool = true,
    filters_row::Union{Nothing, Tuple} = nothing,
    filters_col::Union{Nothing, Tuple} = nothing,
    formatters::Union{Nothing, Function, Tuple} = nothing,
    header::Union{Nothing, AbstractVector, Tuple} = nothing,
    header_alignment::Union{Symbol, Vector{Symbol}} = :s,
    header_cell_alignment::Union{
        Nothing,
        Dict{Tuple{Int, Int}, Symbol},
        Function,
        Tuple
    } = nothing,
    limit_printing::Bool = true,
    renderer::Symbol = :print,
    row_names::Union{Nothing, AbstractVector} = nothing,
    row_name_alignment::Symbol = :r,
    row_name_column_title::AbstractString = "",
    row_number_column_title::AbstractString = "Row",
    show_row_number::Bool = false,
    title::AbstractString = "",
    title_alignment::Symbol = :l
)
    # Get information about the table we have to print based on the format of
    # `data`, which must be an `AbstractMatrix` or an `AbstractVector`.
    dims     = size(data)
    num_dims = length(dims)

    if num_dims == 1
        num_rows = dims[1]
        num_cols = 1
    elseif num_dims == 2
        num_rows = dims[1]
        num_cols = dims[2]
    else
        throw(ArgumentError("`data` must not have more than 2 dimensions."))
    end

    _header = header isa Tuple ? header : (header,)

    header_num_cols = length(first(_header))
    header_num_rows = length(_header)

    if alignment isa Symbol
        alignment = fill(alignment, num_cols)
    else
        if length(alignment) != num_cols
            error("The length of `alignment` must be the same as the number of rows.")
        end
    end

    if header_alignment isa Symbol
        header_alignment = fill(header_alignment, num_cols)
    else
        if length(header_alignment) != num_cols
            error("The length of `header_alignment` must be the same as the number of rows.")
        end
    end

    # If there is a vector of row names, then it must have the same size of the
    # number of rows.
    if row_names !== nothing
        if length(row_names) != num_rows
            error("The number of lines in `row_names` must match the number of lines in the matrix.")
        end

        show_row_names = true
    else
        show_row_names = false
    end

    # If the user wants to filter the data, then check which columns and rows
    # must be printed. Notice that if a data is filtered, then it means that it
    # passed the filter and must be printed.
    filtered_rows = ones(Bool, num_rows)
    filtered_cols = ones(Bool, num_cols)

    if filters_row !== nothing
        @inbounds for i = 1:num_rows
            filtered_i = true

            for filter in filters_row
                !filter(_getdata(data), i) && (filtered_i = false) && break
            end

            filtered_rows[i] = filtered_i
        end
    end

    if filters_col !== nothing
        @inbounds for i = 1:num_cols
            filtered_i = true

            for filter in filters_col
                !filter(_getdata(data), i) && (filtered_i = false) && break
            end

            filtered_cols[i] = filtered_i
        end
    end

    # `id_cols` and `id_rows` contains the indices of the data array that will
    # be printed.
    id_cols          = findall(filtered_cols)
    id_rows          = findall(filtered_rows)
    num_printed_cols = length(id_cols)
    num_printed_rows = length(id_rows)

    # Make sure that `cell_alignment` is a tuple.
    if cell_alignment === nothing
        cell_alignment = ()
    elseif typeof(cell_alignment) <: Dict
        # If it is a `Dict`, then `cell_alignment[(i,j)]` contains the desired
        # alignment for the cell `(i,j)`. Thus, we need to create a wrapper
        # function.
        cell_alignment_dict = copy(cell_alignment)
        cell_alignment = ((data, i, j) -> begin
            if haskey(cell_alignment_dict, (i, j))
                return cell_alignment_dict[(i, j)]
            else
                return nothing
            end
        end,)
    elseif typeof(cell_alignment) <: Function
        cell_alignment = (cell_alignment,)
    end

    # Make sure that `header_cell_alignment` is a tuple.
    if header_cell_alignment === nothing
        header_cell_alignment = ()
    elseif typeof(header_cell_alignment) <: Dict
        # If it is a `Dict`, then `header_cell_alignment[(i,j)]` contains the
        # desired alignment for the cell `(i,j)`. Thus, we need to create a
        # wrapper function.
        header_cell_alignment_dict = copy(header_cell_alignment)
        header_cell_alignment = ((data, i, j) -> begin
            if haskey(header_cell_alignment_dict, (i, j))
                return header_cell_alignment_dict[(i, j)]
            else
                return nothing
            end
        end,)
    elseif typeof(header_cell_alignment) <: Function
        header_cell_alignment = (header_cell_alignment,)
    end

    # Make sure that `formatters` is a tuple.
    formatters === nothing  && (formatters = ())
    typeof(formatters) <: Function && (formatters = (formatters,))

    # Render.
    renderer_val = renderer == :show ? Val(:show) : Val(:print)

    # Create the structure that stores the print information.
    pinfo = PrintInfo(
        data,
        _header,
        id_cols,
        id_rows,
        num_rows,
        num_cols,
        num_printed_cols,
        num_printed_rows,
        header_num_rows,
        header_num_cols,
        show_row_number,
        row_number_column_title,
        show_row_names,
        row_names,
        row_name_alignment,
        row_name_column_title,
        alignment,
        cell_alignment,
        formatters,
        compact_printing,
        title,
        title_alignment,
        header_alignment,
        header_cell_alignment,
        cell_first_line_only,
        renderer_val,
        limit_printing
    )

    return pinfo
end

################################################################################
#                                   Printing
################################################################################

# This is a middleware function to apply the preprocess step to the data that
# will be printed.
function _pretty_table(
    (@nospecialize io::IO),
    data::Any;
    header::Union{Nothing, AbstractVector, Tuple} = nothing,
    kwargs...
)

    if Tables.istable(data)
        if Tables.columnaccess(data)
            pdata, pheader = _preprocess_Tables_column(data, header)
        elseif Tables.rowaccess(data)
            pdata, pheader = _preprocess_Tables_row(data, header)
        else
            error("The object does not have a valid Tables.jl implementation.")
        end

    elseif typeof(data) <: AbstractVecOrMat
        pdata, pheader = _preprocess_vec_or_mat(data, header)
    elseif typeof(data) <: AbstractDict
        sortkeys = get(kwargs, :sortkeys, false)
        pdata, pheader = _preprocess_dict(data; sortkeys = sortkeys)
    else
        error("The type $(typeof(data)) is not supported.")
    end

    return _pt(io, pdata; header = pheader, kwargs...)
end

# This is the low level function that prints the table. In this case, `data`
# must be accessed by `[i,j]` and the size of the `header` must be equal to the
# number of columns in `data`.
function _pt(
    (@nospecialize io::IO),
    data::Any;
    alignment::Union{Symbol, Vector{Symbol}} = :r,
    backend::Union{Symbol, T_BACKENDS} = Val(:auto),
    cell_alignment::Union{
        Nothing,
        Dict{Tuple{Int, Int}, Symbol},
        Function,
        Tuple
    } = nothing,
    cell_first_line_only::Bool = false,
    compact_printing::Bool = true,
    filters_row::Union{Nothing, Tuple} = nothing,
    filters_col::Union{Nothing, Tuple} = nothing,
    formatters::Union{Nothing, Function, Tuple} = nothing,
    header::Union{Nothing, AbstractVector, Tuple} = nothing,
    header_alignment::Union{Symbol, Vector{Symbol}} = :s,
    header_cell_alignment::Union{
        Nothing,
        Dict{Tuple{Int, Int}, Symbol},
        Function,
        Tuple
    } = nothing,
    limit_printing::Bool = true,
    renderer::Symbol = :print,
    row_names::Union{Nothing, AbstractVector} = nothing,
    row_name_alignment::Symbol = :r,
    row_name_column_title::AbstractString = "",
    row_number_column_title::AbstractString = "Row",
    show_row_number::Bool = false,
    title::AbstractString = "",
    title_alignment::Symbol = :l,
    kwargs...
)
    # NOTE: Inform that using `backend` as `Symbol` is deprecated. This will be
    # removed in the next PrettyTables.jl version.
    if backend isa Symbol
      _backend = Val(backend)
      Base.depwarn(
          "Using `backend` as a symbol is deprecated. Use `Val(:$backend)` instead.",
          nothing
      )
    else
      _backend = backend
    end

    if _backend === Val(:auto)
        # In this case, if we do not have the `tf` keyword, then we just
        # fallback to the text _backend. Otherwise, check if the type of `tf`.
        if haskey(kwargs, :tf)
            tf = kwargs[:tf]

            if tf isa TextFormat
                _backend = Val(:text)
            elseif tf isa HTMLTableFormat
                _backend = Val(:html)
            elseif tf isa LatexTableFormat
                _backend = Val(:latex)
            else
                throw(
                    TypeError(
                        :_pt,
                        Union{TextFormat, HTMLTableFormat, LatexTableFormat},
                        typeof(tf)
                    )
                )
            end
        else
            _backend = Val(:text)
        end
    end

    # Create the structure that stores the print information.
    pinfo = _print_info(
        data;
        alignment               = alignment,
        cell_alignment          = cell_alignment,
        cell_first_line_only    = cell_first_line_only,
        compact_printing        = compact_printing,
        filters_row             = filters_row,
        filters_col             = filters_col,
        formatters              = formatters,
        header                  = header,
        header_alignment        = header_alignment,
        header_cell_alignment   = header_cell_alignment,
        limit_printing          = limit_printing,
        renderer                = renderer,
        row_names               = row_names,
        row_name_alignment      = row_name_alignment,
        row_name_column_title   = row_name_column_title,
        row_number_column_title = row_number_column_title,
        show_row_number         = show_row_number,
        title                   = title,
        title_alignment         = title_alignment
    )

    # Select the appropriate backend.
    if _backend === Val(:text)
        _pt_text(Ref{Any}(io), pinfo; kwargs...)
    elseif _backend === Val(:html)
        _pt_html(Ref{Any}(io), pinfo; kwargs...)
    elseif _backend === Val(:latex)
        _pt_latex(Ref{Any}(io), pinfo; kwargs...)
    end

    return nothing
end


text:
# Copyright (c) 2021, Joseph M. Rutherford
# All rights reserved.
#
# Code provided under the license contained in the LICENSE file.
#
# This file depends upon sudoku_common.jl to be included first.

"""
    valid_subarray(rank,test)

Traverse an length rank*rank array of integers to report valid symbols 1:rank*rank.
Rank is required here because we could be getting Vector or Matrix arguments.
"""
function valid_subarray(rank::Integer,test::Array)
    if rank < 1 || Unsigned(rank) != rank
        throw(DomainError("Invalid rank"))
    end
    rank_squared = Unsigned(rank)*Unsigned(rank)
    if length(test) != rank_squared
        throw(DomainError("Array length mismatch to rank-squared"))
    end
    add_entries = Unsigned(sum(test)) == sum(1:rank_squared)
    multiply_entries = prod(test) == factorial(rank_squared)
    return add_entries && multiply_entries
end

"""
    valid_puzzle(puzzle)

Confirm that every row, column, and block of entries in 1:rank*rank.
"""
function valid_puzzle(puzzle::Matrix)
    rank=get_rank(puzzle)
    rank_squared = rank*rank
    result = true
    # use logical AND throughout loops
    for rowcol = 1:rank_squared
        # check full row and col
        if result
            result =
            valid_subarray(rank, puzzle[rowcol,:]) &&
            valid_subarray(rank, puzzle[:,rowcol])
        else
            break # don't bother checking more
        end
    end
    # Check rank x rank blocks
    for row_block = 1:rank
        if result
            start_row = rank*(row_block-1)+1
            stop_row = rank*(row_block-1)+rank
            for col_block = 1:rank
                if result
                    start_col = rank*(col_block-1)+1
                    stop_col = rank*(col_block-1)+rank
                    result = valid_subarray(rank, puzzle[start_row:stop_row,start_col:stop_col])
                else
                    break # don't bother checking more col blocks
                end
            end
        else
            break # don't bother checking more row or col blocks
        end
    end
    return result
end

"""
   satisfies(puzzle,solution)

Verify that a valid solution is a candidate for the puzzle.
"""
function satisfies(puzzle::Matrix,solution::Matrix)
    if size(solution) != size(puzzle)
        throw(DomainError("Solution cannot be checked against a different size puzzle"))
    end
    if !(valid_puzzle(solution))
        throw(DomainError("Invalid solution cannot be checked against a puzzle"))
    end
    satisfies = true
    for i = 1:length(puzzle)
        if puzzle[i] == 0
            continue
        elseif puzzle[i] != solution[i]
            satisfies = false
            break
        end
    end
    return satisfies
end

max_stars_repo_name: joseph-m-rutherford/sudoku-julia

max_issues_repo_name: joseph-m-rutherford/sudoku-julia

max_forks_repo_name: joseph-m-rutherford/sudoku-julia


text:
@discord_object struct IntegrationAccount
    id::String
    name::String
end


text:
@test Nemo.accuracy_bits(GibbsTypePriors.Vnk_NGG(950, 50, 0.5, 0.2, 200)) ≤ 200
@test Nemo.accuracy_bits(GibbsTypePriors.Vnk_NGG(950, 50, 0.5, 0.2, 2000)) ≥ 200
@test Nemo.accuracy_bits(GibbsTypePriors.Vnk_2PD(7, 6, 0.5, 0.01, 200)) ≤ 200
@test Nemo.accuracy_bits(GibbsTypePriors.Vnk_2PD(7, 6, 0.5, 0.01, 2000)) ≥ 200

max_stars_repo_name: konkam/GibbsTypePriors

max_issues_repo_name: konkam/GibbsTypePriors

max_forks_repo_name: konkam/GibbsTypePriors


text:
mutable struct DirectedGraph{V, E} <: AbstractGraph{V, E}
    vertices::Vector{V}
    edges::Vector{E}
    sources::Vector{V}
    targets::Vector{V}
    inedges::Vector{Vector{E}}
    outedges::Vector{Vector{E}}
end

DirectedGraph{V, E}() where {V, E} = DirectedGraph{V, E}(V[], E[], V[], V[], Set{E}[], Set{E}[])

function DirectedGraph(vertexfun::Base.Callable, edgefun::Base.Callable, other::DirectedGraph)
    vertices = map(vertexfun, other.vertices)
    edges = map(edgefun, other.edges)
    vertexmap = Dict(zip(other.vertices, vertices))
    edgemap = Dict(zip(other.edges, edges))
    sources = [vertexmap[v] for v in other.sources]
    targets = [vertexmap[v] for v in other.targets]
    inedges = [[edgemap[e] for e in vec] for vec in other.inedges]
    outedges = [[edgemap[e] for e in vec] for vec in other.outedges]
    DirectedGraph(vertices, edges, sources, targets, inedges, outedges)
end

# AbstractGraph interface
vertices(g::DirectedGraph) = g.vertices
edges(g::DirectedGraph) = g.edges
source(edge, g::DirectedGraph{V, E}) where {V, E} = g.sources[edge_index(edge)]
target(edge, g::DirectedGraph{V, E}) where {V, E} = g.targets[edge_index(edge)]
in_edges(vertex::V, g::DirectedGraph{V, E}) where {V, E} = g.inedges[vertex_index(vertex)]
out_edges(vertex::V, g::DirectedGraph{V, E}) where {V, E} = g.outedges[vertex_index(vertex)]

Base.show(io::IO, ::DirectedGraph{V, E}) where {V, E} = print(io, "DirectedGraph{$V, $E}(…)")

function Base.empty!(g::DirectedGraph)
    empty!(g.vertices)
    empty!(g.edges)
    empty!(g.sources)
    empty!(g.targets)
    empty!(g.inedges)
    empty!(g.outedges)
    g
end

function add_vertex!(g::DirectedGraph{V, E}, vertex::V) where {V, E}
    @assert vertex ∉ vertices(g)
    vertex_index!(vertex, num_vertices(g) + 1)
    push!(g.vertices, vertex)
    push!(g.outedges, E[])
    push!(g.inedges, E[])
    g
end

function add_edge!(g::DirectedGraph{V, E}, source::V, target::V, edge::E) where {V, E}
    @assert edge ∉ edges(g)
    source ∈ vertices(g) || add_vertex!(g, source)
    target ∈ vertices(g) || add_vertex!(g, target)
    edge_index!(edge, num_edges(g) + 1)
    push!(g.edges, edge)
    push!(g.sources, source)
    push!(g.targets, target)
    push!(out_edges(source, g), edge)
    push!(in_edges(target, g), edge)
    g
end

function remove_vertex!(g::DirectedGraph{V, E}, vertex::V) where {V, E}
    disconnected = isempty(in_edges(vertex, g)) && isempty(out_edges(vertex, g))
    disconnected || error("Vertex must be disconnected from the rest of the graph before it can be removed.")
    index = vertex_index(vertex)
    deleteat!(g.vertices, index)
    deleteat!(g.inedges, index)
    deleteat!(g.outedges, index)
    for i = index : num_vertices(g)
        vertex_index!(g.vertices[i], i)
    end
    vertex_index!(vertex, -1)
    g
end

function remove_edge!(g::DirectedGraph{V, E}, edge::E) where {V, E}
    target_inedges = in_edges(target(edge, g), g)
    deleteat!(target_inedges, findfirst(isequal(edge), target_inedges))
    source_outedges = out_edges(source(edge, g), g)
    deleteat!(source_outedges, findfirst(isequal(edge), source_outedges))
    index = edge_index(edge)
    deleteat!(g.edges, index)
    deleteat!(g.sources, index)
    deleteat!(g.targets, index)
    for i = index : num_edges(g)
        edge_index!(g.edges[i], i)
    end
    edge_index!(edge, -1)
    g
end

function rewire!(g::DirectedGraph{V, E}, edge::E, newsource::V, newtarget::V) where {V, E}
    oldsource = source(edge, g)
    oldtarget = target(edge, g)

    g.sources[edge_index(edge)] = newsource
    g.targets[edge_index(edge)] = newtarget

    oldtarget_inedges = in_edges(oldtarget, g)
    deleteat!(oldtarget_inedges, findfirst(isequal(edge), oldtarget_inedges))
    oldsource_outedges = out_edges(oldsource, g)
    deleteat!(oldsource_outedges, findfirst(isequal(edge), oldsource_outedges))

    push!(out_edges(newsource, g), edge)
    push!(in_edges(newtarget, g), edge)

    g
end

function replace_edge!(g::DirectedGraph{V, E}, old_edge::E, new_edge::E) where {V, E}
    if new_edge !== old_edge
        src = source(old_edge, g)
        dest = target(old_edge, g)
        index = edge_index(old_edge)
        edge_index!(new_edge, index)
        g.edges[index] = new_edge
        out_edge_index = findfirst(isequal(old_edge), out_edges(src, g))
        out_edges(src, g)[out_edge_index] = new_edge
        in_edge_index = findfirst(isequal(old_edge), in_edges(dest, g))
        in_edges(dest, g)[in_edge_index] = new_edge
        edge_index!(old_edge, -1)
    end
    g
end

function reindex!(g::DirectedGraph{V, E}, vertices_in_order, edges_in_order) where {V, E}
    @assert isempty(setdiff(vertices(g), vertices_in_order))
    @assert isempty(setdiff(edges(g), edges_in_order))
    sources = source.(edges_in_order, Ref(g))
    targets = target.(edges_in_order, Ref(g))
    empty!(g)
    for v in vertices_in_order
        add_vertex!(g, v)
    end
    for (source, target, e) in zip(sources, targets, edges_in_order)
        add_edge!(g, source, target, e)
    end
    g
end


text:
function get_blocks_list(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/blocks/list.json", options))
end

function get_blocks_ids(auth::TwAuth, options::Dict)
  Requests.json(get_oauth(auth, "https://api.twitter.com/1.1/blocks/ids.json", options))
end

max_stars_repo_name: nick-paul/JuliaTwitter

max_issues_repo_name: nick-paul/JuliaTwitter

max_forks_repo_name: nick-paul/JuliaTwitter


text:
using Plots

function branin(x, a=1, b=5.1/(4π^2), c=5/π, r=6, s=10, t=1/(8π))
    return a * (x[2] - b * x[1]^2 + c * x[1] - r)^2 + s*(1-t) * cos(x[1]) + s;
end

# plot

contour(-50:50, -50:50, (x, y)->branin([x, y]))


text:
# Ported from PRScs/gigrnd.py

function psi(x, alpha, lam)
    f = -alpha*(cosh(x)-1)-lam*(exp(x)-x-1)
    return f
end

function dpsi(x, alpha, lam)
    f = -alpha*sinh(x)-lam*(exp(x)-1)
    return f
end

function g(x, sd, td, f1, f2)
    if (x >= -sd) && (x <= td)
        f = 1
    elseif x > td
        f = f1
    elseif x < -sd
        f = f2
    end

    return f
end

function gigrnd(p, a, b)
    # setup -- sample from the two-parameter version gig(lam,omega)
    p = float(p); a = float(a); b = float(b)
    lam = p
    omega = sqrt(a*b)

    if lam < 0
        lam = -lam
        swap = true
    else
        swap = false
    end

    alpha = sqrt(omega^2+lam^2)-lam

    # find t
    x = -psi(1, alpha, lam)
    if (x >= 1/2) && (x <= 2)
        t = 1
    elseif x > 2
        t = sqrt(2/(alpha+lam))
    elseif x < 1/2
        t = log(4/(alpha+2*lam))
    end

    # find s
    x = -psi(-1, alpha, lam)
    if (x >= 1/2) && (x <= 2)
        s = 1
    elseif x > 2
        s = sqrt(4/(alpha*cosh(1)+lam))
    elseif x < 1/2
        if alpha == 0
            s = 1/lam
        else
            s = min(1/lam, log(1+1/alpha+sqrt(1/alpha^2+2/alpha)))
        end
    end

    # find auxiliary parameters
    eta = -psi(t, alpha, lam)
    zeta = -dpsi(t, alpha, lam)
    theta = -psi(-s, alpha, lam)
    xi = dpsi(-s, alpha, lam)

    p = 1/xi
    r = 1/zeta

    td = t-r*eta
    sd = s-p*theta
    q = td+sd

    # random variate generation
    while true
        U = rand()
        V = rand()
        W = rand()
        if U < q/(p+q+r)
            rnd = -sd+q*V
        elseif U < (q+r)/(p+q+r)
            rnd = td-r*log(V)
        else
            rnd = -sd+p*log(V)
        end

        f1 = exp(-eta-zeta*(rnd-t))
        f2 = exp(-theta+xi*(rnd+s))
        if W*g(rnd, sd, td, f1, f2) <= exp(psi(rnd, alpha, lam))
            break
        end
    end

    # transform back to the three-parameter version gig(p,a,b)
    rnd = exp(rnd)*(lam/omega+sqrt(1+lam^2/omega^2))
    if swap
        rnd = 1/rnd
    end

    rnd = rnd/sqrt(a/b)
    return rnd
end


text:
function test_time_window(
    times, values, windows, op; op_bang=nothing, approximate_equality::Bool=false
)
    # Work out what types we should be using for times and values
    Value = typeof(first(values))
    Time = typeof(first(times))

    for window in windows
        state = if isnothing(op_bang)
            TimeWindowAssociativeOp{Value,op,Time}(window)
        else
            TimeWindowAssociativeOp{Value,op,op_bang,Time}(window)
        end
        for (i, (time, value)) in enumerate(zip(times, values))
            @test update_state!(state, time, value) == state

            # Passing the same time in twice is illegal, since it is not monotonically
            # increasing.
            @test_throws ArgumentError update_state!(state, time, value)

            received = window_value(state)
            # This is the expected range of the input values over which we are reducing.
            i_start = searchsortedlast(times, time - window) + 1
            i_stop = i
            @test window_size(state) == i_stop - i_start + 1
            expected = reduce(op, values[i_start:i_stop])
            if approximate_equality
                @test received ≈ expected
            else
                @test received == expected
            end

            @test window_full(state) == (first(times) <= time - window)
        end
    end
end

@testset "TimeWindowAssociativeOp" begin
    @testset "trivial" begin
        T = Int64
        op = +
        window = 2
        state = TimeWindowAssociativeOp{T,op,typeof(window)}(window)

        @test update_state!(state, 0, 3) == state
        @test !window_full(state)
        @test window_size(state) == 1
        @test window_value(state) == 3

        @test update_state!(state, 1, 4) == state
        @test !window_full(state)
        @test window_size(state) == window
        @test window_value(state) == 7

        @test update_state!(state, 2, 4) == state
        @test window_full(state)
        @test window_size(state) == window
        @test window_value(state) == 8

        @test update_state!(state, 3, -5) == state
        @test window_full(state)
        @test window_size(state) == window
        @test window_value(state) == -1
    end

    @testset "incompatible time type" begin
        # Time is of type Int64, but window is given as a Float64 -> this should throw.
        @test_throws ArgumentError TimeWindowAssociativeOp{Float64,+,Int64}(10.0)
    end

    @testset "integral" begin
        for op in (+, *, max, min)
            test_time_window(1:20, 1:20, 1:20, op)
            test_time_window(
                [Dates.DateTime(2000, 1, x) for x in 1:20],
                1:20,
                [Dates.Day(x) for x in 1:30],
                op,
            )
        end
    end

    @testset "random matrix multiply" begin
        times = [DateTime(2000, 1, x) for x in 1:20]
        values = [rand(-5:5, (2, 2)) for _ in 1:20]
        windows = [Dates.Day(x) for x in 1:30]
        test_time_window(times, values, windows, *)
    end

    @testset "mutating" begin
        times = [DateTime(2000, 1, x) for x in 1:20]
        values = [Data(1, rand(-5:5, (2, 2))) for _ in 1:20]
        windows = [Dates.Day(x) for x in 1:30]
        test_time_window(times, values, windows, merge)
        test_time_window(times, values, windows, merge; op_bang=merge!)
    end
end


text:
# script to fit model to data
#
# The script is called as
#
# fitmodel.jl [model] [objective] [[rat]]
#
# and writes data to ../fitdata/fits/fit_[model]_[objective]_[rat]. Here, [rat]
# is optional and indicates the rat identifier (1-4). If not provided, the data
# for all rats is collapsed into a single rat.

using Random, CSV

# include required scripts
include("common.jl")
include("models.jl")
include("simmodel.jl")
include("behaviorstats.jl")
include("modelfitting.jl")

# settings
const blockreps = 10
const simtrials = 10^5
const burnin = 1000
const nrats = 1 # (do not divide dataset by 4 rats, nrats = 4)
const datapath = "../fitdata/fits"
const paraminiobjfn = "psychchronseq1"

# returns model for given name
function createmodel(modelname)
    modeldict = Dict("opt" => OptimalModel,
                     "optlapse" => OptimalLapseModel,
                     "optcollapse" => OptimalColLapseModel,
                     "optbiaslapse" => OptimalBiasLapseModel,
                     "optcolbiaslapse" => OptimalColBiasLapseModel,
                     "adflapse" => ADFLapseModel,
                     "adfcollapse" => ADFColLapseModel,
                     "adfbiaslapse" => ADFBiasLapseModel,
                     "adfcolbiaslapse" => ADFColBiasLapseModel,
                     "randomcollapse" => RandomColLapseModel,
                     "randomcolbiaslapse" => RandomColBiasLapseModel,
                     "gammalapse" => GammaLapseModel,
                     "gammacollapse" => GammaColLapseModel,
                     "gammabiaslapse" => GammaBiasLapseModel,
                     "gammacolbiaslapse" => GammaColBiasLapseModel,
                     "deltalapse" => DeltaLapseModel,
                     "deltacollapse" => DeltaColLapseModel,
                     "deltabiaslapse" => DeltaBiasLapseModel,
                     "deltacolbiaslapse" => DeltaColBiasLapseModel,
                     "deltaexpcolbiaslapse" => DeltaExpColBiasLapseModel,
                     "lcalapse" => LCALapseModel,
                     "lcacollapse" => LCAColLapseModel)
    if modelname in keys(modeldict)
        return modeldict[modelname]()
    else
        error("Uknown model identifier $modelname")
    end
end

# returns fitting objective for given name and model
function objectivefn(objfnname, m, ratid::Int=0)
    # load data for individual / all rats
    d_ident = load_data(:ident, ratid)
    d_categ = load_data(:categ, ratid)
    d_itlvd = load_data(:itlvd, ratid)
    # simulated trials
    Random.seed!(0)
    evidence, ids = catevidence(ident_cs, categ_cs, 10, 10^5)
    itlvdevidence, itlvdids = catitlvdevidence(itlvd_cs, 10, 10^5)
    # create objective function
    if objfnname == "psychchron1"
        o = PsychChron1ObjFn(m, evidence, ids, burnin, d_ident, d_categ, nrats)
    elseif objfnname == "psychchron1ident"
        o = PsychChron1IdentObjFn(m, evidence, ids, burnin, d_ident, d_categ, nrats)
    elseif objfnname == "psychchron1categ"
        o = PsychChron1CategObjFn(m, evidence, ids, burnin, d_ident, d_categ, nrats)
    elseif objfnname == "psychchron2"
        o = PsychChron2ObjFn(m, evidence, ids, burnin, d_ident, d_categ, nrats)
    elseif objfnname == "psychchron3"
        o = PsychChron3ObjFn(m, evidence, ids, burnin, d_ident, d_categ, nrats)
    elseif objfnname == "psychchronseq1"
        o = PsychChronSeq1ObjFn(m, evidence, ids, burnin, d_ident, d_categ, nrats)
    elseif objfnname == "psychchronseq2"
        o = PsychChronSeq2ObjFn(m, evidence, ids, burnin, d_ident, d_categ, nrats)
    elseif objfnname == "psychchronseq3"
        o = PsychChronSeq3ObjFn(m, evidence, ids, burnin, d_ident, d_categ, nrats)
    elseif objfnname == "psychchron1itlvd"
        o = PsychChron1ItlvdObjFn(m, itlvdevidence, itlvdids, burnin, d_itlvd, nrats)
    elseif objfnname == "psychchron1itlvdedge"
        # only fit stimuli that also appear in ident/categ condition
        o = PsychChron1ItlvdObjFn(m, itlvdevidence, itlvdids, burnin, d_itlvd, nrats,
                                  [1, 2, 3, 4, 5, 6, 7, 8, 9, 16, 17, 24, 25, 32])
    else
        error("Unknown object function identifier $objfnname")
    end
end

# process arguments
if length(ARGS) ∉ [2,3]
    error("Expected two or three arguments")
end
const modelname = lowercase(ARGS[1])
const objfnname = lowercase(ARGS[2])
ratid = length(ARGS) >= 3 ? parse(Int, ARGS[3]) : 0  # ratid = 0 for all rats
m = createmodel(modelname)
o = objectivefn(objfnname, m, ratid)

# check if model has already been fit with parameter initialization objective fn
datafile(modelname, objfnname, ratid) = ratid == 0 ?
    "$(datapath)/fit_$(modelname)_$(objfnname).csv" :
    "$(datapath)/fit_$(modelname)_$(objfnname)_rat$(ratid).csv"
paraminifilename = datafile(modelname, paraminiobjfn, ratid)
if ratid != 0 && !isfile(paraminifilename)
    # default to collective rat fits if they don't exist for individual rats
    paraminifilename = datafile(modelname, paraminiobjfn, 0)
end
if isfile(paraminifilename)
    # loading initial parameters from previous model fit with $paraminiobjfn
    println("Loading initial parameters from $paraminifilename")
    d = CSV.read(paraminifilename)
    (maxllh, maxi) = findmax(d[:llh])
    if maxi != size(d,1) && maxllh > d[end,end]
        warn("Maximum llh $maxllh found in line $maxi instead of last line")
    end
    ϕini = Vector(d[maxi,1:end-1])
else
    # use initial parameters are determined by model
    ϕini = initialϕ(m)
end

# fit model
if ratid == 0
    println("Fitting $modelname model with objective $objfnname across all rats")
else
    println("Fitting $modelname model with objective $objfnname for rat $ratid")
end
d, ret = fitmodel(m, o, ϕini)
println("Optimizer stopped with return value $ret")
println("Minimum found at")
println(d[end,:])

if !(ret ∈ (:SUCCESS, :XTOL_REACHED, :STOPVAL_REACHED,
            :FTOL_REACHED, :MAXEVAL_REACHED, :MAXTIME_REACHED))
    println("Not writing to file, as optimization didn't complete successfully")
else
    # write to file
    outfile = datafile(modelname, objfnname, ratid)
    println("Writing data to $(outfile)")
    CSV.write(outfile, d)
end


text:
const DEFAULT_VERSION = 0

"Converts between bitstypes, and little endian to host."
function convert_array(T::Type{<:Base.BitInteger64}, src::Vector{<:Base.BitInteger64})
    nbytes = length(src) * sizeof(eltype(src))
    dst = Vector{T}(undef, nbytes ÷ sizeof(T))
    GC.@preserve dst src begin
        srcptr = Ptr{UInt}(pointer(src))
        dstptr = Ptr{UInt}(pointer(dst))
        for i in 1:nbytes ÷ sizeof(UInt)
            unsafe_store!(dstptr, htol(unsafe_load(srcptr, i)), i)
        end
    end
    return dst
end

convert_hashfunction(sk::KmerHashes{K}) where K = KmerHashes{K,hash}(sk.bases, sk.sketch)
convert_hashfunction(sk::Vector{<:KmerHashes}) = [convert_hashfunction(s) for s in sk]

@noinline function validate_hashfunction(F)
    F === hash || error("""Kash currently only supports serializing with hashfunction
    Base.hash, see `convert_hashfunction`""")
end

function bson_dict(v::Val{0}, x::KmerHashes)
    d = Dict{Symbol,Any}()
    d[:r] = Int64(x.sketch.requested)
    d[:b] = Int64(x.bases)
    d[:h] = convert_array(UInt8, x.sketch.hashes)
    return d
end

function bson_dict(v::Val{0}, sketches::Vector{KmerHashes{K,F}}) where {K,F}
    validate_hashfunction(F)
    d = Dict{Symbol,Any}()
    d[:m] = Int32(1752392011) # UTF-8 for "Kash", magic header
    d[:v] = Int32(0)
    d[:k] = Int32(K)

    array = Vector{Any}() # must be a Vector{Any} for BSON
    for sketch in sketches
        push!(array, bson_dict(Val(0), sketch))
    end
    d[:h] = array
    return d
end

# This function should only verify magic header and check version.
# Then it dispatches to the version-specific parser
function load_sketches_bson(dict::Dict{Symbol,Any})
    if !haskey(dict, :m) || (dict[:m] != Int32(1752392011))
        throw(ValueError("Magic number 1752392011 not found in BSON"))
    end
    v = Int(dict[:v])
    if !(v in (0,))
        throw(ValueError("This version of Kash cannot read Kash BSON file v. $v"))
    end
    return load_sketches_bson(Val{v}(), dict)
end

function load_sketches_bson(v::Val{0}, dict::Dict{Symbol,Any})
    K = Int(dict[:k])
    T = KmerHashes{K,hash}
    array = dict[:h]
    sketches = Vector{T}(undef, length(array))
    for (i, sketchdict) in enumerate(array)
        sketch = load_sketch_bson(Val(0), T, sketchdict)
        sketches[i] = sketch
    end
    return sketches
end

function load_sketch_bson(v::Val{0}, T::Type{<:KmerHashes}, dict::Dict{Symbol,Any})
    r = Int(dict[:r])
    b = Int(dict[:b])
    h = convert_array(UInt, dict[:h])
    sketch = MinHash.MinHashSketch(hash, r, h)
    return T(b, sketch)
end

function load_sketches(io::IO)
    d = BSON.parse(io)
    load_sketches_bson(d)
end

load_sketches(path::String) = open(load_sketches, path)

function save_sketches(io::IO, sketches::Vector{KmerHashes{K,F}}) where {K,F}
    d = bson_dict(Val(DEFAULT_VERSION), sketches)
    bson(io, d)
end

function save_sketches(path::String, sketches::Vector{KmerHashes{K,F}}) where {K,F}
    open(io -> save_sketches(io, sketches), path, "w")
end


text:
"""
    AlternateSelector <: AbstractActionSelector

Used to ensure that all actions are selected alternatively.

    AlternateSelector(n::Int)

`n::Int` means the optional actions are `1:n`.
"""
mutable struct AlternateSelector <: AbstractActionSelector
    n::Int
    cur::Int
    AlternateSelector(n::Int) = new(n, 0)
end

function (s::AlternateSelector)()
    s.cur = s.cur + 1 > s.n ? 1 : s.cur + 1
    s.cur
end

"""
    (s::AlternateSelector)(values::Any)

Ignore the action `values`, generate an action alternatively.

## Example

```julia
julia> selector = AlternateSelector(3)
AlternateSelector(3, 0)

julia> any_state = 0 # for AlternateSelector, state can be anything

julia> [selector(any_state) for i in 1:10]  # iterate through all actions
10-element Array{Int64,1}:
 1
 2
 3
 1
 2
 3
 1
 2
 3
 1
```
"""
(s::AlternateSelector)(values::Any) = s()


text:
using Revise
	using PseudoArcLengthContinuation, LinearAlgebra, Plots, SparseArrays, Setfield, Parameters
	const PALC = PseudoArcLengthContinuation

f1(u, v) = u^2 * v
norminf = x -> norm(x, Inf)

function plotsol(x; kwargs...)
	N = div(length(x), 2)
	plot!(x[1:N], label="u"; kwargs...)
	plot!(x[N+1:2N], label="v"; kwargs...)
end

function Fbru!(f, x, p)
	@unpack α, β, D1, D2, l = p
	n = div(length(x), 2)
	h = 1.0 / (n+0); h2 = h*h
	c1 = D1 / l^2 / h2
	c2 = D2 / l^2 / h2

	u = @view x[1:n]
	v = @view x[n+1:2n]

	# Dirichlet boundary conditions
	f[1]   = c1 * (α	  - 2u[1] + u[2] ) + α - (β + 1) * u[1] + f1(u[1], v[1])
	f[end] = c2 * (v[n-1] - 2v[n] + β / α)			 + β * u[n] - f1(u[n], v[n])

	f[n]   = c1 * (u[n-1] - 2u[n] +  α   ) + α - (β + 1) * u[n] + f1(u[n], v[n])
	f[n+1] = c2 * (β / α  - 2v[1] + v[2])			 + β * u[1] - f1(u[1], v[1])

	for i=2:n-1
		  f[i] = c1 * (u[i-1] - 2u[i] + u[i+1]) + α - (β + 1) * u[i] + f1(u[i], v[i])
		f[n+i] = c2 * (v[i-1] - 2v[i] + v[i+1])			  + β * u[i] - f1(u[i], v[i])
	end
	return f
end

function Fbru(x, p)
	f = similar(x)
	Fbru!(f, x, p)
end

function Jbru_sp(x, p)
	@unpack α, β, D1, D2, l = p
	# compute the Jacobian using a sparse representation
	n = div(length(x), 2)
	h = 1.0 / n; h2 = h*h

	c1 = D1 / p.l^2 / h2
	c2 = D2 / p.l^2 / h2

	u = @view x[1:n]
	v = @view x[n+1:2n]

	diag   = zeros(eltype(x), 2n)
	diagp1 = zeros(eltype(x), 2n-1)
	diagm1 = zeros(eltype(x), 2n-1)

	diagpn = zeros(eltype(x), n)
	diagmn = zeros(eltype(x), n)

	@. diagmn = β - 2 * u * v
	@. diagm1[1:n-1] = c1
	@. diagm1[n+1:end] = c2

	@. diag[1:n]    = -2c1 - (β + 1) + 2 * u * v
	@. diag[n+1:2n] = -2c2 - u * u

	@. diagp1[1:n-1]   = c1
	@. diagp1[n+1:end] = c2

	@. diagpn = u * u
	J = spdiagm(0 => diag, 1 => diagp1, -1 => diagm1, n => diagpn, -n => diagmn)
	return J
end

function finalise_solution(z, tau, step, contResult)
	n = div(length(z.u), 2)
	printstyled(color=:red, "--> Solution constant = ", norm(diff(z.u[1:n])), " - ", norm(diff(z.u[n+1:2n])), "\n")
	return true
end

n = 100
####################################################################################################
# test for the Jacobian expression
# using ForwardDiff
# sol0 = rand(2n)
# J0 = ForwardDiff.jacobian(x-> Fbru(x, par_bru), sol0) |> sparse
# J1 = Jbru_sp(sol0, par_bru)
# J0 - J1
####################################################################################################
# different parameters to define the Brusselator model and guess for the stationary solution
par_bru = (α = 2., β = 5.45, D1 = 0.008, D2 = 0.004, l = 0.3)
	sol0 = vcat(par_bru.α * ones(n), par_bru.β/par_bru.α * ones(n))

# par_bru = (α = 2., β = 4.6, D1 = 0.0016, D2 = 0.008, l = 0.061)
# 	xspace = LinRange(0, par_bru.l, n)
# 	sol0 = vcat(		par_bru.α .+ 2 .* sin.(pi*xspace/par_bru.l),
# 			par_bru.β/par_bru.α .- 0.5 .* sin.(pi*xspace/par_bru.l))

# eigls = EigArpack(1.1, :LM)
# 	opt_newton = PALC.NewtonPar(eigsolver = eigls)
# 	out, hist, flag = @time PALC.newton(
# 		x ->  Fbru(x, par_bru),
# 		x -> Jbru_sp(x, par_bru),
# 		sol0, opt_newton, normN = norminf)
#
# 		plot();plotsol(out);plotsol(sol0, label = "sol0",line=:dash)
####################################################################################################
eigls = EigArpack(1.1, :LM)
opts_br_eq = ContinuationPar(dsmin = 0.001, dsmax = 0.02, ds = 0.01, pMax = 1.9, detectBifurcation = 2, nev = 21, plotEveryNsteps = 50, newtonOptions = NewtonPar(eigsolver = eigls, tol = 1e-9), maxSteps = 1060)

	br, _ = @time PALC.continuation(
		(x, p) ->    Fbru(x, @set par_bru.l = p),
		(x, p) -> Jbru_sp(x, @set par_bru.l = p),
		sol0, par_bru.l,
		opts_br_eq, verbosity = 0,
		plot = false,
		printSolution = (x, p) -> x[div(n,2)], normC = norminf)
#################################################################################################### Continuation of the Hopf Point using Jacobian expression
ind_hopf = 1
	# hopfpt = PALC.HopfPoint(br, ind_hopf)

	hopfpoint, _, flag = @time PALC.newtonHopf(
		(x, p) ->  Fbru(x, @set par_bru.l = p),
		(x, p) -> Jbru_sp(x, @set par_bru.l = p),
		br, ind_hopf,
		opts_br_eq.newtonOptions, normN = norminf)
	flag && printstyled(color=:red, "--> We found a Hopf Point at l = ", hopfpoint.p[1], ", ω = ", hopfpoint.p[2], ", from l = ", br.bifpoint[ind_hopf].param, "\n")
####################################################################################################Continuation of Periodic Orbit
M = 10

l_hopf, Th, orbitguess2, hopfpt, vec_hopf = PALC.guessFromHopf(br, ind_hopf, opts_br_eq.newtonOptions.eigsolver, M, 22*0.05)
#
orbitguess_f2 = reduce(hcat, orbitguess2)
orbitguess_f = vcat(vec(orbitguess_f2), Th) |> vec
####################################################################################################
# essai Shooting
using DifferentialEquations, DiffEqOperators, ForwardDiff

FOde(f, x, p, t) = Fbru!(f, x, p)
Jbru(x, dx, p) = ForwardDiff.derivative(t -> Fbru(x .+ t .* dx, p), 0.)
JacOde(J, x, p, t) = copyto!(J, Jbru_sp(x, p))

u0 = vec(sol0) .+ 0.01 .* rand(2n)
par_hopf = (@set par_bru.l = l_hopf + 0.01)

jac_prototype = Jbru_sp(ones(sol0 |> length), @set par_bru.β = 0)
	jac_prototype.nzval .= ones(length(jac_prototype.nzval))

ff = ODEFunction(FOde, jac_prototype = JacVecOperator{Float64}(FOde, u0, par_hopf))
probsundials = ODEProblem(FOde, u0, (0., 5200.), par_hopf)

vf = ODEFunction((u,p,t)->Fbru(u,p); jac = (u,p,t) -> Jbru_sp(u,p))
	prob = ODEProblem(vf,  u0, (0.0, .1), par_hopf)

# heatmap(sol[:,:], color = :viridis)
####################################################################################################
# M = 10
dM = 3
orbitsection = Array(orbitguess_f2[:, 1:dM:M])

initpo = vcat(vec(orbitsection), 3.0)

PALC.plotPeriodicShooting(initpo[1:end-1], length(1:dM:M));title!("")

# PALC.sectionShooting(initpo, Array(orbitguess_f2[:,1:dM:M]), par_hopf)

probSh = p -> PALC.ShootingProblem(u -> Fbru(u, p), p, prob, Rodas4P(),
		length(1:dM:M), x -> PALC.sectionShooting(x, Array(orbitguess_f2[:,1:dM:M]), p, Fbru); atol = 1e-10, rtol = 1e-8)

res = @time probSh(par_hopf)(initpo)
norminf(res)
res = probSh(par_hopf)(initpo, initpo)
norminf(res)

ls = GMRESIterativeSolvers(tol = 1e-7, N = length(initpo), maxiter = 100, verbose = false)
	# ls = GMRESKrylovKit{Float64}(verbose = 0, dim = 200, atol = 1e-9, rtol = 1e-5)
	optn_po = NewtonPar(verbose = true, tol = 1e-9,  maxIter = 25, linsolver = ls)
	# deflationOp = PALC.DeflationOperator(2.0, (x,y) -> dot(x[1:end-1], y[1:end-1]),1.0, [outpo])
	outpo ,_ = @time PALC.newton(probSh(par_hopf),
			initpo, optn_po;
			normN = norminf)
	plot(initpo[1:end-1], label = "Init guess")
	plot!(outpo[1:end-1], label = "sol")

opts_po_cont = ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds= 0.01, pMax = 1.5, maxSteps = 500, newtonOptions = (@set optn_po.tol = 1e-7), nev = 25, precisionStability = 1e-8, detectBifurcation = 0)
	br_po, _, _= @time PALC.continuationPOShooting(
		p -> probSh(@set par_hopf.l = p),
		outpo, par_hopf.l,
		opts_po_cont; verbosity = 2,
		plot = true,
		plotSolution = (x; kwargs...) -> PALC.plotPeriodicShooting!(x[1:end-1], length(1:dM:M); kwargs...),
		printSolution = (u, p) -> u[end], normC = norminf)

# simplified call
eig = EigKrylovKit(tol= 1e-12, x₀ = rand(2n), verbose = 0, dim = 40)
# eig = DefaultEig()
opts_po_cont_floquet = @set opts_po_cont.newtonOptions = NewtonPar(linsolver = ls, eigsolver = eig, tol = 1e-7, verbose = true)
opts_po_cont_floquet = setproperties(opts_po_cont_floquet; nev = 3, precisionStability = 1e-2, detectBifurcation = 2, maxSteps = 5000)

br_po, _ , _ = @time PALC.continuationPOShooting(
		p -> probSh(@set par_hopf.l = p),
		outpo, par_hopf.l,
		opts_po_cont_floquet; verbosity = 2,
		plot = true,
		# callbackN = cb_ss,
		plotSolution = (x; kwargs...) -> PALC.plotPeriodicShooting!(x[1:end-1], length(1:dM:M); kwargs...),
		printSolution = (u, p) -> u[end], normC = norminf)
####################################################################################################
# Simple Poincare Shooting
function sectionP(x, po::AbstractMatrix, p)
	res = 1.0
	N, M = size(po)
	# @show size(x,2)
	for ii = 1:size(x, 2)
		res *= dot(x .- po[:, ii], Fbru(po[:, ii], p))
	end
	# @show p res
	res
end

M = 1
orbitsection = Array(sol[:,1:M])

initpo = vec(orbitsection)

sectionP(initpo[1:2n], Array(sol[:,1:M]), par_hopf)

probPSh = p -> PALC.PoincareShootingProblem(u -> Fbru(u, p), p, probsundials, Rodas4P(), M, x -> sectionP(x, Array(sol[:,1:M]), p), atol = 1e-9)

probPSh(par_hopf)

res = @time probPSh(par_hopf)(initpo)
norminf(res)
res = probPSh(par_hopf)(initpo, initpo)
norminf(res)

ls = GMRESKrylovKit(verbose = 0, dim = 200, atol = 1e-9, rtol = 1e-5)
	ls = GMRESIterativeSolvers(tol = 1e-4, N = length(initpo), maxiter=50, verbose = false)
	optn = NewtonPar(verbose = true, tol = 1e-9,  maxIter = 65, linsolver = ls)
	# deflationOp = PALC.DeflationOperator(2.0, (x,y) -> dot(x[1:end-1], y[1:end-1]),1.0, [outpo])
	outpo, _ = @time PALC.newton(x -> probPSh(par_hopf)(x),
			x -> (dx -> probPSh(par_hopf)(x, dx)),
			initpo,
			optn;
			# callback = cb_ss,
			normN = norminf)
	plot(initpo[1:end-1], label = "Init guess")
	plot!(outpo[1:end-1], label = "sol")

getPeriod(probPSh(par_hopf), outpo)

opts_poP_cont = ContinuationPar(dsmin = 0.0001, dsmax = 0.02, ds= 0.02, pMax = 4.0, maxSteps = 3, newtonOptions = @set optn.tol = 1e-8)
	eig = EigKrylovKit(tol= 1e-12, x₀ = rand(2n), verbose = 2, dim = 40)
	# eig = DefaultEig()
	opts_poP_cont_floquet = setproperties(opts_poP_cont.newtonOptions; eigsolver = eig,maxIter = 25)
	opts_poP_cont_floquet = setproperties(opts_poP_cont; computeEigenValues = true, nev = 10, precisionStability = 5e-3, detectBifurcation = 1, maxSteps = 100)

br_pok2, upo , _= @time PALC.continuationPOShooting(
		p -> probPSh(@set par_hopf.l = p),
		outpo, par_hopf.l,
		opts_poP_cont_floquet;
		verbosity = 2, plot = true,
		# callbackN = cb_ss,
		plotSolution = (x;kwargs...) -> plot!(x; label = "", kwargs...),
		# printSolution = (x,p)->norm(x),
		normC = norminf)
####################################################################################################
# Multiple Poincare Shooting
function sectionP!(out, x, po::AbstractMatrix, p)
	res = 1.0
	N, M = size(po)
	# @show N M
	for ii = 1:M
		out[ii] = dot(x .- po[:, ii], Fbru(po[:, ii], p))
	end
	out
end

M = 1
orbitsection = Array(sol[:,1:M])

initpo = vec(orbitsection)

sectionP(initpo[1:2n], Array(sol[:,1:M]), par_hopf)

probPSh = p -> PALC.PoincareShootingProblem(u -> Fbru(u, p), p, probsundials, Rodas4P(), M, x -> sectionP!(x, Array(sol[:,1:M]), p), atol = 1e-9)

probPSh(par_hopf)

res = @time probPSh(par_hopf)(initpo)
norminf(res)
res = probPSh(par_hopf)(initpo, initpo)
norminf(res)

ls = GMRESIterativeSolvers(tol = 1e-4, N = length(initpo), maxiter=50, verbose = false)
	optn = NewtonPar(verbose = true, tol = 1e-9,  maxIter = 65, linsolver = ls)
	# deflationOp = PALC.DeflationOperator(2.0, (x,y) -> dot(x[1:end-1], y[1:end-1]),1.0, [outpo])
	outpo, _ = @time PALC.newton(x -> probPSh(par_hopf)(x),
			x -> (dx -> probPSh(par_hopf)(x, dx)),
			initpo,
			optn;
			callback = cb_ss,
			normN = norminf)
	plot(initpo[1:end-1], label = "Init guess")
	plot!(outpo[1:end-1], label = "sol")

getPeriod(probPSh(par_hopf), outpo)
####################################################################################################
# Multiple Poincare Shooting with Hyperplane parametrization
M = size(sol, 2)
dM = 15
normals = [Fbru(sol[:,ii], par_hopf)/(norm(Fbru(sol[:,ii], par_hopf)))^2 for ii = 1:dM:M]
	centers = [sol[:,ii] for ii = 1:dM:M]

probHPsh = p -> PALC.PoincareShootingProblem(u -> Fbru(u, p), p, probmatrixfree, ImplicitEuler(), normals, centers; atol = 1e-10, rtol = 1e-9, dt = 0.1)

# probHPsh = p -> PALC.PoincareShootingProblem(u -> Fbru(u, p), p, probsundials, ImplicitEuler(), normals, centers; atol = 1e-10, rtol = 1e-9, dt = 0.1)

hyper = probHPsh(par_hopf).psh.section
initpo_bar = zeros(size(sol,1)-1, length(normals))
	for ii=1:length(normals)
		initpo_bar[:, ii] .= PALC.R(hyper, centers[ii], ii)
	end

probHPsh(par_hopf)(vec(initpo_bar))
probHPsh(par_hopf)(vec(initpo_bar)) |> norminf

probHPsh(par_hopf)(vec(initpo_bar), vec(initpo_bar) )

ls = GMRESIterativeSolvers(tol = 1e-5, N = length(vec(initpo_bar)), maxiter = 500, verbose = true)
	# ls = GMRESKrylovKit(rtol = 1e-5, verbose = 2)
	deflationOp = PALC.DeflationOperator(2.0, (x,y) -> dot(x,y), 1.0, [sol0])
	optn = NewtonPar(verbose = true, tol = 1e-7,  maxIter = 140, linsolver = ls)
	outpo_psh, _ = @time PALC.newton(x -> probHPsh(par_hopf)(x; verbose = false),
			x -> (dx -> probHPsh(par_hopf)(x, dx; δ = 1e-7)),
			vec(initpo_bar), optn,
			# deflationOp,
			; normN = norminf)

plot(outpo_psh)

####################################################################################################
using DiffEqBase

function sectionHyper(out, x, normals = normals, centers = centers)
	for ii = 1:length(normals)
		out[ii] = dot(normals[ii], x .- centers[ii])
	end
	out
endw

orbitsection = vec(sol[:,1:1])

pSection(out, u, t, integrator) = sectionHyper(out, u) * (integrator.iter > 1)
affect!(integrator, idx) = terminate!(integrator)
cb = VectorContinuousCallback(pSection, affect!, length(normals); interp_points = 150, affect_neg! = nothing)
# _probsundials = DiffEqBase.remake(probsundials; u0 = vec(sol[:,1:1]), tspan = (0, Inf64), p = par_hopf)

PALC.flow(vec(sol[:,1:1]), Inf64, probsundials; alg = Rodas4(), callback = cb)

PALC.dflow_fd(vec(sol[:,1:1]), vec(sol[:,1:1]), Inf64, probsundials; alg = Rodas4(), callback = cb)