text:
module M

function make_matrix(n, s=5.0)
    A = randn(n,n)
    for i = 1:n
        A[i,i] += s
    end
    A
end

nonsymm2symm(A) = (A + A')/2.0
make_symmat(n) = nonsymm2symm(make_matrix(n))

type MatrixDecomp
    nb::Int
    ni::Int
    Abb::Array{Float64,2}
    Aii::Array{Float64,2}
    Abi::Array{Float64,2}
    Aib::Array{Float64,2}

    function MatrixDecomp(A, nb, ni)
        Abb = A[1:nb, 1:nb]
        Abi = A[1:nb, (nb+1):end]
        Aib = A[(nb+1):end, 1:nb]
        Aii = A[(nb+1):end, (nb+1):end]

        new(nb, ni, Abb, Aii, Abi, Aib)
    end
end

using Base.LinAlg.LAPACK.getrf!
using Base.LinAlg.LAPACK.getrs!
using Base.LinAlg.BLAS.gemm!
using Base.LinAlg.BLAS.gemv!

type StaticCond
    nb::Int
    ni::Int
    B::Array{Float64,2}
    iAii::Array{Float64,2}
    X::Array{Float64,2}
    Aib::Array{Float64,2}
    ipiv::Array{Int,1}
    ipiv2::Array{Int,1}
    function StaticCond(A::MatrixDecomp)
        nb = A.nb
        ni = A.ni
        B = zeros(nb, nb)
        iAii = copy(A.Aii)
        iAii, ipiv, info = getrf!(iAii)
        X = A.Abi'
        getrs!('T', iAii, ipiv, X)

        B = copy(A.Abb)
        gemm!('T', 'N', -1.0, X, A.Aib, 1.0, B)
        B, ipiv2, info = getrf!(B)

        new(nb, ni, B, iAii, X, copy(A.Aib), ipiv, ipiv2)
    end
end
using Base.LinAlg.LAPACK.potrf!
using Base.LinAlg.LAPACK.potrs!

type ChoStaticCond
    nb::Int
    ni::Int
    B::Array{Float64,2}
    iAii::Array{Float64,2}
    X::Array{Float64,2}
    function ChoStaticCond(A::MatrixDecomp)
        nb = A.nb
        ni = A.ni
        iAii = copy(A.Aii)
        iAii, info = potrf!('L', iAii)
        X = copy(A.Aib)
        potrs!('L', iAii, X)

        B = copy(A.Abb)
        gemm!('T', 'N', -1.0, X, A.Aib, 1.0, B)

        B, info = potrf!('L', B)

        new(nb, ni, B, iAii, X)
    end
end

function mysolve(A::StaticCond, f)
    nb = A.nb
    ni = A.ni

    fb = f[1:nb]
    fi = f[(nb+1):end]

    gemv!('T', -1.0, A.X, fi, 1.0, fb)

    getrs!('N', A.B, A.ipiv2, fb)

    gemv!('N', -1.0, A.Aib, fb, 1.0, fi)

    getrs!('N', A.iAii, A.ipiv, fi)

    [fb, fi]
end

function mysolve(A::ChoStaticCond, f)
    nb = A.nb
    ni = A.ni

    fb = f[1:nb]
    fi = f[(nb+1):end]

    gemv!('T', -1.0, A.X, fi, 1.0, fb)
    potrs!('L', A.B, fb)

    potrs!('L', A.iAii, fi)
    gemv!('N', -1.0, A.X, fb, 1.0, fi)

    [fb, fi]
end

function teste(n=15, nb=10, ni=5)

    A = make_symmat(n)
    Abi = MatrixDecomp(A, nb, ni)

    iA = inv(A)
    f = linspace(1, n, n)

    x = iA * f

    S1 = StaticCond(Abi)

    x1 = mysolve(S1, f)

    S2 = ChoStaticCond(Abi)
    x2 = mysolve(S2, f)

    hcat(x, x1, x2)
end

end # Module


text:
# =========================================================================================
"""
    writexml3d_json{T}(
        stream::IOStream,
        model ::Union{Vector{Vector{T}}, Model{T, Triangle{T}}}
    )

Creates a XML3D-specific JSON file either from a given collection of nodes (representing the
latter as point cloud) or from a given surface model.

# Specification
<https://github.com/xml3d/xml3d.js/wiki/External-resources>

# Return type
`Void`

# Alias

    writexml3d_json{T}(
        fname::String,
        nodes::Union{Vector{Vector{T}}), Model{T, Triangle{T}}}
    )

Creates the JSON file by name rather than `IOStream` object.
"""
function writexml3d_json(
        stream::IOStream,
        nodes ::Vector{Vector{T}}
    ) where T
    println(stream, json(Dict(
        "format" => "xml3d-json",
        "version" => "0.4.0",
        "data" => Dict(
            "position" => Dict(
                "type" => "float3",
                "seq" => [Dict{String, Vector{Float64}}("value" => unpack(nodes))]
            )
        )
    )))
    nothing
end

function writexml3d_json(
        stream::IOStream,
        model ::Model{T, Triangle{T}}
    ) where T
    revidx = reverseindex(model.nodes)
    println(stream, json(Dict(
        "format" => "xml3d-json",
        "version" => "0.4.0",
        "data" => Dict(
            "index" => Dict(
                "type" => "int",
                "seq" => [Dict{String, Vector{Int}}(
                            "value" => [revidx[objectid(n)]-1 for n
                                in unpack([Vector{T}[o.v1, o.v2, o.v3] for o
                                in model.elements])]
                          )]
            ),
            "position" => Dict(
                "type" => "float3",
                "seq" => [Dict{String, Vector{Float64}}("value" => unpack(model.nodes))]
            ),
            "normal" => Dict(
                "type" => "float3",
                "seq" => [Dict{String, Vector{Float64}}(
                            "value" => unpack(vertexnormals(model))
                          )]
            )
        )
    )))
    nothing
end

function writexml3d_json(
        fname::String,
        model::Union{Vector{Vector{T}},Model{T, Triangle{T}}}
    ) where T
    open(fh -> writexml3d_json(fh, model), fname, "w")
end

# =========================================================================================
"""
    writexml3d_xml{T}(
        stream::IOStream,
        nodes ::Vector{Vector{T}}
    )

Creates a XML3D-specific XML file from a given collection of nodes, representing the latter
as point cloud.

# Specification
<https://github.com/xml3d/xml3d.js/wiki/External-resources>

# Return type
`Void`

# Alias

    writexml3d_xml{T}(
        fname::String,
        nodes::Vector{Vector{T}}
    )

Creates the XML file by name rather than `IOStream` object.
"""
function writexml3d_xml(
        stream::IOStream,
        nodes ::Vector{Vector{T}}
    ) where T
    xdoc = XMLDocument()
    xroot = create_root(xdoc, "xml3d")
    set_attribute(xroot, "xmlns", "http://www.xml3d.org/2009/xml3d")
    xmesh = new_child(xroot, "data")
    set_attribute(xmesh, "id", "mesh")
    xpos = new_child(xmesh, "float3")
    set_attribute(xpos, "name", "position")
    add_text(xpos, join(unpack(nodes), " "))
    println(stream, string(xdoc))
    nothing
end

function writexml3d_xml(
        fname::String,
        nodes::Vector{Vector{T}}
    ) where T
    open(fh -> writexml3d_xml(fh, nodes), fname, "w")
end

ext: jl


text:
# ------------------------------------------------------------
# Similar situation can happen with array types

function slow_positives(xs)
    ys = []   # this is actually Vector{Any}
    for x in xs
        if x > 0
            push!(ys, x)
        end
    end
    return ys
end

function fast_positives(xs)
    ys = similar(xs, 0)
    for x in xs
        if x > 0
            push!(ys, x)
        end
    end
    return ys
end

#=
julia> x = rand(Int, 100000);

julia> @btime slow_positives($x);
  1.632 ms (49730 allocations: 1.76 MiB)

julia> @btime fast_positives($x);
  1.023 ms (16 allocations: 1.00 MiB)
=#

# huge number of allocations!


text:
"""
    branch!(tree, node)

Get the branching variable with [`get_branching_variable`](@ref) and then calls [`get_branching_nodes_info`](@ref) and [`add_node!`](@ref).
"""
function branch!(tree, node)
    variable_idx = get_branching_variable(tree, tree.options.branch_strategy, node)
    # no branching variable selected => return
    variable_idx == -1 && return
    nodes_info = get_branching_nodes_info(tree, node, variable_idx)
    for node_info in nodes_info
        add_node!(tree, node_info)
    end
end

"""
    get_branching_variable(tree::BnBTree, ::FIRST, node::AbstractNode)

Return the first possible branching variable which is a branching variable based on `tree.branching_indices`
and is currently not valid based on [`is_approx_feasible`](@ref).
Return `-1` if all integer constraints are respected.
"""
function get_branching_variable(tree::BnBTree, ::FIRST, node::AbstractNode)
    values = get_relaxed_values(tree, node)
    for i in tree.branching_indices
        value = values[i]
        if !is_approx_feasible(tree, value)
            return i
        end
    end
    return -1
end

"""
    get_branching_variable(tree::BnBTree, ::MOST_INFEASIBLE, node::AbstractNode)

Return the branching variable which is furthest away from being feasible based on [`get_distance_to_feasible`](@ref)
or `-1` if all integer constraints are respected.
"""
function get_branching_variable(tree::BnBTree, ::MOST_INFEASIBLE, node::AbstractNode)
    values = get_relaxed_values(tree, node)
    best_idx = -1
    max_distance_to_feasible = 0.0
    for i in tree.branching_indices
        value = values[i]
        if !is_approx_feasible(tree, value)
            distance_to_feasible = get_distance_to_feasible(tree, value)
            if distance_to_feasible > max_distance_to_feasible
                best_idx = i
                max_distance_to_feasible = distance_to_feasible
            end
        end
    end
    return best_idx
end

"""
    get_branching_nodes_info(tree::BnBTree, node::AbstractNode, vidx::Int)

Create the information for new branching nodes based on the variable index `vidx`.
Return a list of those information as a `NamedTuple` vector.

# Example
The following would add the necessary information about a new node and return it. The necessary information are the fields required by the [`AbstractNode`](@ref).
For this examle the required fields are the lower and upper bounds of the variables as well as the status of the node.
```julia
nodes_info = NamedTuple[]
push!(nodes_info, (
    lbs = lbs,
    ubs = ubs,
    status = MOI.OPTIMIZE_NOT_CALLED,
))
return nodes_info
```
"""
function get_branching_nodes_info end


text:
using FEMAssembler, LinearAlgebra
using Test

@testset "Test: linear von mises problem" begin
  include( joinpath( "..", "examples", "vonMisesTruss", "vonMisesTruss.jl" ) )

  L0 = sqrt( d^2 + h^2 ) ;

  cosalpha = d / L0 ;
  sinalpha = h / L0 ;

  Ldef12 = norm( [ d+UG[7],h+UG[11]] - [0,0]  ) ;
  Ldef23 = norm( [2d,0] - [ d+UG[7],h+UG[11]] ) ;

  F12Analy =  +Fx/(2cosalpha) + Fy/(2sinalpha) ;
  F23Analy =  -Fx/(2cosalpha) + Fy/(2sinalpha) ;

  delta12Num = Ldef12-L0 ;
  delta23Num = Ldef23-L0 ;

  delta12Analy = F12Analy*L0 / (E*A) ;
  delta23Analy = F23Analy*L0 / (E*A) ;

  @test maximum( abs.( [ delta12Num-delta12Analy, delta12Num-delta12Analy ] ) ) < 1e-8
end


text:
"""
Provides the functions related to generating documentation stubs.
"""
module Generator

using DocStringExtensions

"""
$(SIGNATURES)

Attempts to save a file at `\$(root)/\$(filename)`. `f` will be called with file
stream (see [`open`](https://docs.julialang.org/en/latest/base/io-network/#Base.open)).

`filename` can also be a file in a subdirectory (e.g. `src/index.md`), and then
then subdirectories will be created automatically.
"""
function savefile(f, root, filename)
    filepath = joinpath(root, filename)
    if ispath(filepath) error("$(filepath) already exists") end
    @info("Generating $filename at $filepath")
    mkpath(dirname(filepath))
    open(f,filepath,"w")
end

"""
$(SIGNATURES)

Contents of the default `make.jl` file.
"""
function make(pkgname; format = :html)
    fmtpkg = format === :markdown ? ", DocumenterMarkdown" :
             format === :pdf      ? ", DocumenterLaTeX" : ""
    fmtstr = format === :html ? "Documenter.HTML()" :
             format === :markdown  ? "Markdown()" :
             format === :pdf  ? "LaTeX()" : ""

    sitename = format !== :markdown ? "\n    sitename = \"$(pkgname)\"," : ""
    """
    using Documenter$(fmtpkg)
    using $(pkgname)

    makedocs($(sitename)
        format = $(fmtstr),
        modules = [$(pkgname)]
    )

    # Documenter can also automatically deploy documentation to gh-pages.
    # See "Hosting Documentation" and deploydocs() in the Documenter manual
    # for more information.
    #=deploydocs(
        repo = "<repository url>"
    )=#
    """
end

"""
$(SIGNATURES)

Contents of the default `.gitignore` file.
"""
function gitignore()
    """
    build/
    site/
    """
end

"""
$(SIGNATURES)

Contents of the default `Project.toml` file.
"""
function project(;format = :html)
    deps = Dict("Documenter" => "e30172f5-a6a5-5a46-863b-614d45cd2de4")
    if format === :markdown
        deps["DocumenterMarkdown"] = "997ab1e6-3595-5248-9280-8efb232c3433"
    elseif format === :pdf
        deps["DocumenterLaTeX"] = "cd674d7a-5f81-5cf3-af33-235ef1834b99"
    end
    io = IOBuffer()
    print(io, "[deps]\n")
    for dep in deps
        print(io, dep.first, " = \"", dep.second, "\"\n")
    end

    return String(take!(io))
end

mkdocs_default(name, value, default) = value == nothing ? "#$name$default" : "$name$value"

"""
$(SIGNATURES)

Contents of the default `mkdocs.yml` file.
"""
function mkdocs(pkgname;
        description = nothing,
        author = nothing,
        url = nothing
    )
    s = """
    # See the mkdocs user guide for more information on these settings.
    #   http://www.mkdocs.org/user-guide/configuration/

    site_name:        $(pkgname).jl
    $(mkdocs_default("repo_url:         ", url, "https://github.com/USER_NAME/PACKAGE_NAME.jl"))
    $(mkdocs_default("site_description: ", description, "Description..."))
    $(mkdocs_default("site_author:      ", author, "USER_NAME"))

    theme: readthedocs

    extra_css:
      - assets/Documenter.css

    extra_javascript:
      - https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS_HTML
      - assets/mathjaxhelper.js

    markdown_extensions:
      - extra
      - tables
      - fenced_code
      - mdx_math

    docs_dir: 'build'

    pages:
      - Home: index.md
    """
end

"""
$(SIGNATURES)

Contents of the default `src/index.md` file.
"""
function index(pkgname)
    """
    # $(pkgname).jl

    Documentation for $(pkgname).jl
    """
end

end


text:
json_parse(r::HTTP.Response) = r.body |> String |> JSON.parse

function http_error(response)
    status_number = response.status
    status_text = response.statustext
    result_string = response.body |> String
    error("$status_number $status_text: $result_string")
end

ask(f, remote, url, body = "") =
    f(string(base_url(remote), url), headers = headers(remote), body = body)

function talk_to(f, remote, url, body = "")
    response = ask(f, remote, url, body)
    status_number = response.status
    if status_number >= 300
        http_error(response)
    end
    response
end

user_repo(r::Remote) = "$(r.username)/$(r.repo_name)"


text:
using MasterEquationSolve

rm("precompile_outdir",force=true,recursive=true)

MasterEquationSolve.mesolve_qutip("precompile_input","precompile_outdir")


text:
tol = 1e-10
#
# Test grid points.
#
Teff, L_star, logg = stellar_evolution(0.0900000036,t=1e5,Z=0.01)
@test abs(Teff   - 10^3.4680  ) < tol
@test abs(L_star - 10^(-0.733)) < tol
@test abs(logg   - 2.949) < tol

Teff, L_star, logg = stellar_evolution(0.75,t=1.41e9,Z=0.01)
@test abs(Teff   - 10^3.6942  ) < tol
@test abs(L_star - 10^(-0.606)) < tol
@test abs(logg   - 4.649) < tol

#
# Test interpolation across mass.
#
m = (0.0900000036 + 0.1000000015)/2

Teff, L_star, logg = stellar_evolution(m,t=1e5,Z=0.01)
@test abs(Teff   - 10^((3.4680 + 3.4723)/2) ) < tol
@test abs(L_star - 10^(-(0.733 + 0.658 )/2) ) < tol
@test abs(logg   - (2.949 + 2.937)/2 ) < tol

#
# Test interpolation across time.
#
m = 0.5
t = (1.00e5 + 1.12e+05) / 2

Teff, L_star, logg = stellar_evolution(m,t=t,Z=0.01)
@test abs(Teff   - 10^((3.6132 + 3.6131)/2) ) < tol
@test abs(L_star - 10^((0.557  + 0.543 )/2) ) < tol
@test abs(logg   - (2.985 + 2.999)/2 ) < tol

#
# Test interpolation across metallicity.
#
m = 0.5
t = 1e5
Z = (0.01 + 0.015) / 2

Teff, L_star, logg = stellar_evolution(m,t=t,Z=Z)
@test abs(Teff   - 10^((3.6132 + 3.6034)/2) ) < tol
@test abs(L_star - 10^((0.557  + 0.524 )/2) ) < tol
@test abs(logg   - (2.985 + 2.978)/2 ) < tol


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VpnClientConfiguration <: SwaggerModel
    vpnClientAddressPool::Any # spec type: Union{ Nothing, AddressSpace } # spec name: vpnClientAddressPool
    vpnClientRootCertificates::Any # spec type: Union{ Nothing, Vector{VpnClientRootCertificate} } # spec name: vpnClientRootCertificates
    vpnClientRevokedCertificates::Any # spec type: Union{ Nothing, Vector{VpnClientRevokedCertificate} } # spec name: vpnClientRevokedCertificates
    vpnClientProtocols::Any # spec type: Union{ Nothing, Vector{String} } # spec name: vpnClientProtocols
    vpnClientIpsecPolicies::Any # spec type: Union{ Nothing, Vector{IpsecPolicy} } # spec name: vpnClientIpsecPolicies
    radiusServerAddress::Any # spec type: Union{ Nothing, String } # spec name: radiusServerAddress
    radiusServerSecret::Any # spec type: Union{ Nothing, String } # spec name: radiusServerSecret
    radiusServers::Any # spec type: Union{ Nothing, Vector{RadiusServer} } # spec name: radiusServers
    aadTenant::Any # spec type: Union{ Nothing, String } # spec name: aadTenant
    aadAudience::Any # spec type: Union{ Nothing, String } # spec name: aadAudience
    aadIssuer::Any # spec type: Union{ Nothing, String } # spec name: aadIssuer

    function VpnClientConfiguration(;vpnClientAddressPool=nothing, vpnClientRootCertificates=nothing, vpnClientRevokedCertificates=nothing, vpnClientProtocols=nothing, vpnClientIpsecPolicies=nothing, radiusServerAddress=nothing, radiusServerSecret=nothing, radiusServers=nothing, aadTenant=nothing, aadAudience=nothing, aadIssuer=nothing)
        o = new()
        validate_property(VpnClientConfiguration, Symbol("vpnClientAddressPool"), vpnClientAddressPool)
        setfield!(o, Symbol("vpnClientAddressPool"), vpnClientAddressPool)
        validate_property(VpnClientConfiguration, Symbol("vpnClientRootCertificates"), vpnClientRootCertificates)
        setfield!(o, Symbol("vpnClientRootCertificates"), vpnClientRootCertificates)
        validate_property(VpnClientConfiguration, Symbol("vpnClientRevokedCertificates"), vpnClientRevokedCertificates)
        setfield!(o, Symbol("vpnClientRevokedCertificates"), vpnClientRevokedCertificates)
        validate_property(VpnClientConfiguration, Symbol("vpnClientProtocols"), vpnClientProtocols)
        setfield!(o, Symbol("vpnClientProtocols"), vpnClientProtocols)
        validate_property(VpnClientConfiguration, Symbol("vpnClientIpsecPolicies"), vpnClientIpsecPolicies)
        setfield!(o, Symbol("vpnClientIpsecPolicies"), vpnClientIpsecPolicies)
        validate_property(VpnClientConfiguration, Symbol("radiusServerAddress"), radiusServerAddress)
        setfield!(o, Symbol("radiusServerAddress"), radiusServerAddress)
        validate_property(VpnClientConfiguration, Symbol("radiusServerSecret"), radiusServerSecret)
        setfield!(o, Symbol("radiusServerSecret"), radiusServerSecret)
        validate_property(VpnClientConfiguration, Symbol("radiusServers"), radiusServers)
        setfield!(o, Symbol("radiusServers"), radiusServers)
        validate_property(VpnClientConfiguration, Symbol("aadTenant"), aadTenant)
        setfield!(o, Symbol("aadTenant"), aadTenant)
        validate_property(VpnClientConfiguration, Symbol("aadAudience"), aadAudience)
        setfield!(o, Symbol("aadAudience"), aadAudience)
        validate_property(VpnClientConfiguration, Symbol("aadIssuer"), aadIssuer)
        setfield!(o, Symbol("aadIssuer"), aadIssuer)
        o
    end
end # type VpnClientConfiguration

const _property_map_VpnClientConfiguration = Dict{Symbol,Symbol}(Symbol("vpnClientAddressPool")=>Symbol("vpnClientAddressPool"), Symbol("vpnClientRootCertificates")=>Symbol("vpnClientRootCertificates"), Symbol("vpnClientRevokedCertificates")=>Symbol("vpnClientRevokedCertificates"), Symbol("vpnClientProtocols")=>Symbol("vpnClientProtocols"), Symbol("vpnClientIpsecPolicies")=>Symbol("vpnClientIpsecPolicies"), Symbol("radiusServerAddress")=>Symbol("radiusServerAddress"), Symbol("radiusServerSecret")=>Symbol("radiusServerSecret"), Symbol("radiusServers")=>Symbol("radiusServers"), Symbol("aadTenant")=>Symbol("aadTenant"), Symbol("aadAudience")=>Symbol("aadAudience"), Symbol("aadIssuer")=>Symbol("aadIssuer"))
const _property_types_VpnClientConfiguration = Dict{Symbol,String}(Symbol("vpnClientAddressPool")=>"AddressSpace", Symbol("vpnClientRootCertificates")=>"Vector{VpnClientRootCertificate}", Symbol("vpnClientRevokedCertificates")=>"Vector{VpnClientRevokedCertificate}", Symbol("vpnClientProtocols")=>"Vector{String}", Symbol("vpnClientIpsecPolicies")=>"Vector{IpsecPolicy}", Symbol("radiusServerAddress")=>"String", Symbol("radiusServerSecret")=>"String", Symbol("radiusServers")=>"Vector{RadiusServer}", Symbol("aadTenant")=>"String", Symbol("aadAudience")=>"String", Symbol("aadIssuer")=>"String")
Base.propertynames(::Type{ VpnClientConfiguration }) = collect(keys(_property_map_VpnClientConfiguration))
Swagger.property_type(::Type{ VpnClientConfiguration }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VpnClientConfiguration[name]))}
Swagger.field_name(::Type{ VpnClientConfiguration }, property_name::Symbol) =  _property_map_VpnClientConfiguration[property_name]

const _allowed_VpnClientConfiguration_vpnClientProtocols = ["IkeV2", "SSTP", "OpenVPN"]

function check_required(o::VpnClientConfiguration)
    true
end

function validate_property(::Type{ VpnClientConfiguration }, name::Symbol, val)
    if name === Symbol("vpnClientProtocols")
        Swagger.validate_param(name, "VpnClientConfiguration", :enum, val, _allowed_VpnClientConfiguration_vpnClientProtocols)
    end
end


text:
abstract type SequentialPointProcess end
@doc raw"""
  Softcore(kernel, kernel_integral=kernel)

Define a Softcore model with kernel.

The kernel is a positive function of distance. Large values indicate repulsion.

The desity for the next point ``x_k`` is
```math
f(x_k | x_{1:k-1}) \propto \exp \left(-\sum_{i=1}^{k-1} f(\| x_i - x_k \|)\right).
```

```kernel_integral``` can be specified to
use a different kernel for computing the normalization coefficients.
"""
struct Softcore{F,F2} <: SequentialPointProcess
  kernel::F
  kernel_integral::F2
end
Softcore(f) = Softcore(f, f)
"""
  OverlappingDiscs(theta, R)

  Define a sequential model where the probability of the next point depends on
  how many discs overlap that location. Each point generates a disc `R(i)` and the
  radius of the disc may depend on the id of the point. The
  parameter `theta(noverlaps, k)` is the density corresponding the number
  of overlaps and the id of the next point.

  theta and R may be specified as numbers, vectors of length equal to the
  number of points.
"""
struct OverlappingDiscs{F, F2} <: SequentialPointProcess
  theta::F
  R::F2
  function OverlappingDiscs(theta, R)
    f, f2 = totheta(theta), to1fun(R)
    new{typeof(f), typeof(f2)}(f, f2)
  end
end
to1fun(x::AbstractVector{<:Real}) = k -> x[k]
to1fun(x::Real) = k -> x
to1fun(x) = x
totheta(x::AbstractVector{<:Real}) = (noverlaps, k) -> ifelse(noverlaps==0, 1-x[k], x[k])
totheta(x::Real) = (noverlaps, k) -> ifelse(noverlaps==0, 1-x, x)
totheta(x) = x

"""Uniform()

Uniform density.
"""
struct Uniform <: SequentialPointProcess end
"""Mixture(model1, model2, theta)

Defines a Mixture model with components ```model1``` and ```model2``` with
probibilities ```(1-theta)``` and ```theta```.
"""
struct Mixture{T1, T2} <: SequentialPointProcess
  first::T1
  second::T2
  theta::Float64
end

function normalizedlogdensities(m, xy, window, int)
  z = normconstants(m, xy[1:end-1], window, int)
  for i in eachindex(z)
    if z[i] == 0.0
      z[i] = -Inf
    else
      z[i] = logfk(m, xy[i+1], view(xy, 1:i)) - log(z[i])
    end
  end
  z
end
function normalizedlogdensities(::Uniform, xy, window, int)
  x1, x2 = window.x
  y1, y2 = window.y
  A = (x2-x1)*(y2-y1)
  fill(log(1/A), length(xy)-1)
end
function logdens(m::Mixture, xy, window, cubpar)
  0.0 <= m.theta <= 1.0 || return -Inf
  z1 = normalizedlogdensities(m.first, xy, window, cubpar)
  z2 = normalizedlogdensities(m.second, xy, window, cubpar)
  sum(zip(z1, z2)) do (z1, z2)
    log((1-m.theta) * exp(z1) + m.theta * exp(z2))
  end
end
function logdens(m, xy, window, cubpar)
  z1 = normalizedlogdensities(m, xy, window, cubpar)
  sum(z1)
end

max_stars_repo_name: mikkoku/SequentialPointProcesses

max_issues_repo_name: mikkoku/SequentialPointProcesses

max_forks_repo_name: mikkoku/SequentialPointProcesses


text:
module RulesMvNormalMeanPrecisionPrecisionTest

using Test
using ReactiveMP
using Random

import ReactiveMP: @test_rules

@testset "rules:MvNormalMeanPrecision:precision" begin

    @testset "Variational: (q_out::PointMass, q_μ::MultivariateNormalDistributionsFamily)" begin

        @test_rules [ with_float_conversions = true, float32_atol = 1e-5 ] MvNormalMeanPrecision(:Λ, Marginalisation) [
            (input = (q_out = PointMass([ 1.0, 2.0 ]), q_μ = MvNormalMeanPrecision([ 3.0, 5.0 ], [ 3.0 2.0; 2.0 4.0 ])), output = Wishart(4.0, [ 75/73 -46/73; -46/73 36/73 ])),
            (input = (q_out = MvNormalMeanPrecision([1.0; 2.0], [3.0 2.0; 2.0 4.0]), q_μ = PointMass([3.0; 5.0])), output = Wishart(4.0, [ 75/73 -46/73; -46/73 36/73 ])),
            (input = (q_out = MvNormalMeanPrecision([1.0; 2.0], [3.0 2.0; 2.0 4.0]), q_μ = MvNormalMeanPrecision([3.0; 5.0], [3.0 2.0; 2.0 4.0])), output = Wishart(4.0, [ 39/74 -11/37; -11/37 10/37 ])),
            (input = (q_out = MvNormalMeanPrecision([1.0; 2.0; 3.0], [3.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 4.0]), q_μ = MvNormalMeanPrecision([3.0; 5.0; -2.0], [3.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 4.0])), output = Wishart(5.0, [ 15/11 -3/22 5/11; -3/22 19/22 5/11; 5/11 5/11 16/33 ]))
        ]
        @test_rules [ with_float_conversions = true, float32_atol = 1e-5 ] MvNormalMeanPrecision(:Λ, Marginalisation) [
            (input = (q_out = PointMass([ 1.0, 2.0 ]), q_μ = MvNormalMeanCovariance([ 3.0, 5.0 ], [ 3.0 2.0; 2.0 4.0 ])), output = Wishart(4.0, [ 13/27 -8/27; -8/27 7/27 ])),
            (input = (q_out = MvNormalMeanCovariance([1.0; 2.0], [3.0 2.0; 2.0 4.0]), q_μ = PointMass([3.0; 5.0])), output = Wishart(4.0, [ 13/27 -8/27; -8/27 7/27 ])),
            (input = (q_out = MvNormalMeanCovariance([1.0; 2.0], [3.0 2.0; 2.0 4.0]), q_μ = MvNormalMeanCovariance([3.0; 5.0], [3.0 2.0; 2.0 4.0])), output = Wishart(4.0, [ 17/70 -1/7; -1/7 1/7 ])),
            (input = (q_out = MvNormalMeanCovariance([1.0; 2.0; 3.0], [3.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 4.0]), q_μ = MvNormalMeanCovariance([3.0; 5.0; -2.0], [3.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 4.0])), output = Wishart(5.0, [ 51/338 -6/169 5/169; -6/169 115/676 45/676; 5/169 45/676 47/676 ]))
        ]
        @test_rules [ with_float_conversions = true, float32_atol = 1e-5 ] MvNormalMeanPrecision(:Λ, Marginalisation) [
            (input = (q_out = PointMass([ 1.0, 2.0 ]), q_μ = MvNormalWeightedMeanPrecision([ 3.0, 5.0 ], [ 3.0 2.0; 2.0 4.0 ])), output = Wishart(4.0, [ 73/67 -26/67; -26/67 68/67 ])),
            (input = (q_out = MvNormalWeightedMeanPrecision([1.0; 2.0], [3.0 2.0; 2.0 4.0]), q_μ = PointMass([3.0; 5.0])), output = Wishart(4.0, [ 165/163 -106/163; -106/163 76/163 ])),
            (input = (q_out = MvNormalWeightedMeanPrecision([1.0; 2.0], [3.0 2.0; 2.0 4.0]), q_μ = MvNormalWeightedMeanPrecision([3.0; 5.0], [3.0 2.0; 2.0 4.0])), output = Wishart(4.0, [ 73/70 11/35; 11/35 34/35 ])),
            (input = (q_out = MvNormalWeightedMeanPrecision([1.0; 2.0; 3.0], [3.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 4.0]), q_μ = MvNormalWeightedMeanPrecision([3.0; 5.0; -2.0], [3.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 4.0])), output = Wishart(5.0, [ 459/338 -36/169 60/169; -36/169 115/169 90/169; 60/169 90/169 188/169 ]))
        ]

    end

    @testset "Variational: (q_out_μ::MultivariateNormalDistributionsFamily)" begin

        @test_rules [ with_float_conversions = true ] MvNormalMeanPrecision(:Λ, Marginalisation) [
            (input = (q_out_μ = MvNormalMeanCovariance(ones(4), diageye(4)), ), output = Wishart(4.0, [ 0.5 0.0; 0.0 0.5 ])),
            (input = (q_out_μ = MvNormalMeanPrecision(ones(4), diageye(4)), ), output = Wishart(4.0, [ 0.5 0.0; 0.0 0.5 ])),
            (input = (q_out_μ = MvNormalWeightedMeanPrecision(ones(4), diageye(4)), ), output = Wishart(4.0, [ 0.5 0.0; 0.0 0.5 ])),
        ]

        a = [ 1.0, 2.0, -1.0, -2.0 ]
        A = [ 3.5 -0.5 -0.25 0.0; -0.5 3.0 -0.25 0.0; -0.25 -0.25 6.0 0.25; 0.0 0.0 0.25 7.0]

        @test_rules [ with_float_conversions = true, float32_atol = 1e-5 ] MvNormalMeanPrecision(:Λ, Marginalisation) [
            (input = (q_out_μ = MvNormalMeanCovariance(a, A),), output = Wishart(4.0, [ 13/150 -2/75; -2/75 7/150])),
            (input = (q_out_μ = MvNormalMeanPrecision(a, A), ), output = Wishart(4.0, [ 3751/1966  -3653/3932; -3653/3932  30259/58980 ])),
            (input = (q_out_μ = MvNormalWeightedMeanPrecision(a, A), ), output = Wishart(4.0, [ 139109613/68523766 -96827813/137047532; -96827813/137047532 62742851/68523766 ])),
        ]

    end

end

end


text:
using Polynomials
"""
    jacobi(p::Integer, α, β)
Compute the Legendre polynomial with parameters `α`, `β` of degree `p`
using the three term recursion [Karniadakis and Sherwin, Spectral/hp Element
Methods for CFD, Appendix A].
"""
function poly_jacobi(p::Integer, α, β, ::Type{T}=Float64, var=:x) where {T}
    a = Poly{T}([one(T)], var)
    b = Poly{T}([(α - β)/2, (2+α+β)/2], var)
    if p <= 0
        return a
    elseif p == 1
        return b
    end
    px = Poly{T}([zero(T), one(T)], var)
    for n in 2:p
        a1 = 2n*(n+α+β)*(2n-2+α+β)
        a2 = (2n-1+α+β)*(α+β)*(α-β)
        a3 = (2n-2+α+β)*(2n-1+α+β)*(2n+α+β)
        a4 = 2*(n-1+α)*(n-1+β)*(2n+α+β)
        a, b = b, ( (a2+a3*px)*b - a4*a ) / a1
    end
    b
end

function poly_legendre{T<:Number}(n, ::Type{T}=Float64, var=:x)
    return poly_jacobi(n,0.0,0.0,T,var)
end

function dubiner_base(j::Integer,z,w)
    S=2*w-1
    R=(2*(1+(2*z-1))/(1-S))-1
    t=-3/2+(1/2)*sqrt(1+8*j)
    n=((ceil(t)+1)*(ceil(t)+2))/2-j;
    k=2*n+1;
    m=ceil(t)-n;

    @assert ((w>=0)&&(z>=0))&&(1>=z+w) "evalued point not in domain"
    P=polyval(poly_jacobi(Int(n),0,0),R)*polyval(poly_jacobi(Int(m),k,0),S)*((1-S)/2)^n
    N=sqrt(2/((2*n+1)*(m+n+1)))
    return (2*P)/N
end

function poly_jacobi2(x, p::Integer, α, β)
    T = typeof( (2+α+β)*x / 2 )
    a = one(T)
    b = ((2+α+β)*x + α - β) / 2
    if p <= 0
        return a
    elseif p == 1
        return b
    end
    for n in 2:p
        a1 = 2n*(n+α+β)*(2n-2+α+β)
        a2 = (2n-1+α+β)*(α+β)*(α-β)
        a3 = (2n-2+α+β)*(2n-1+α+β)*(2n+α+β)
        a4 = 2*(n-1+α)*(n-1+β)*(2n+α+β)
        a, b = b, ( (a2+a3*x)*b - a4*a ) / a1
    end
    b
end

function poly_legendre2(x::T, n::Integer) where {T}
    return poly_jacobi2(x, n,zero(T),zero(T))
end

function dubiner_base2(j::Integer,z,w)
    S=2*w-1
    R=(2*(1+(2*z-1))/(1-S))-1
    t=-3/2+(1/2)*sqrt(1+8*j)
    n=((ceil(t)+1)*(ceil(t)+2))/2-j;
    k=2*n+1;
    m=ceil(t)-n;
    @assert ((w>=0)&&(z>=0))&&(1>=z+w) "evalued point not in domain"
    P=poly_jacobi2(R,Int(n),0,0)*poly_jacobi2(S, Int(m),k,0)*((1-S)/2)^n
    N=sqrt(2/((2*n+1)*(m+n+1)))
    return (2*P)/N
end

using PolynomialBases
"""
djacobi(x,n::Integer,α, β)
Evaluate the `n` order jacobi polynomial derivative
at point x∈[-1,1]
"""
function djacobi(x,n::Integer,α, β)
    T = typeof((2+α+β)*x/2)
    a = zero(T)
    if n <= 0;return a;end
    return (n+α+β+1)/2*jacobi(x,n-1,α+1,β+1)
end
"""
dubiner(x,y,n::Integer,m::Integer)
Compute the dubiner polynomial of degree `n`,`m` at point (x,y)
on the reference triangle ((0,0),(1,0),(0,1))
"""
function dubiner(x,y,n::Integer,m::Integer)
    #check domain
    @assert ((y>=0)&&(x>=0))&&(1>=x+y) "point not in domain"
    # Map to reference square
    ξ=2*x/(1-y)-1
    η=2*y-1
    k=2*n+1
    #Compute Dubiner_nm(ξ, η)
    P=jacobi(ξ,n,0,0)*jacobi(η, m,k,0)*((1-η)/2)^n
    #normalize
    N=sqrt(2/((2*n+1)*(m+n+1)))
    return (2*P)/N
end
"""
eval_dubiner(x,y,j::Integer)
Evaluate the dubiner basis `j` at point (x,y)
on the reference triangle ((0,0),(1,0),(0,1))
"""
function dubiner_basis(x,y,j::Integer)
    #Compute degrees
    t=-3/2+(1/2)*sqrt(1+8*j)
    n=((ceil(t)+1)*(ceil(t)+2))/2-j
    m=ceil(t)-n
    #Compute Dubiner_nm(ξ, η)
    return dubiner(x,y,Int(n),Int(m))
end

"""
∇dubiner(x,y,n::Integer,m::Integer)
Compute the gradient of dubiner polynomial of degree `(n,m)` at point (x,y)
on the reference triangle ((0,0),(1,0),(0,1))
"""
function ∇dubiner(x,y,n::Integer,m::Integer)
    #check domain
    @assert ((y>=0)&&(x>=0))&&(1>=x+y) "point not in domain"
    # Map to reference square
    ξ=2*x/(1-y)-1
    k=2*n+1
    η=2*y-1
    #Compute ∇Dubiner_nm(ξ, η)
    Dφn = djacobi(ξ,n,0,0)
    φn = jacobi(ξ,n,0,0)
    φm = jacobi(η,m,k,0)
    Dφm = djacobi(η,m,k,0)
    Px=2/(1-η)*Dφn*φm*((1-η)/2)^n
    N=sqrt(2/((2*n+1)*(m+n+1)))
    Py=(2*x/(1-y)^2*Dφn*((1-η)/2)^n-n*((1-η)/2)^(n-1)*φn)*φm + 2*φn*((1-η)/2)^n*Dφm
    return [(4*Px)/N,(2*Py)/N]
end

"""
∇dubiner_basis(x,y,j::Integer)
Compute the gradient of dubiner basis `j` at point (x,y)
on the reference triangle ((0,0),(1,0),(0,1))
"""
function ∇dubiner_basis(x,y,j::Integer)
    #Compute degrees
    t=-3/2+(1/2)*sqrt(1+8*j)
    n=Int(((ceil(t)+1)*(ceil(t)+2))/2-j)
    m=Int(ceil(t)-n)
    #Compute ∇Dubiner_nm(ξ, η)
    return ∇dubiner(x,y,n,m)
end

"""
eval_dubiner_ex(z,w,j::Integer)
Compute the dubiner polynomial of degree `j` at point (r,s)
on the reference triangle ((0,0),(1,0),(0,1))
Use exact functions when available
"""
function eval_dubiner_ex(r::T,s::T,j::Integer) where{T}
    if j == 0; return zero(T)
    elseif j == 1; return sqrt(2)*one(T)
    elseif j == 2; return 2*sqrt(3)*(2*r + s - 1)
    elseif j == 3; return 2*(3*s - 1)
    elseif j == 4; return sqrt(30)*(6*r^2 + 6*r*(s - 1) + s^2 - 2*s + 1)
    elseif j == 5; return 3*sqrt(2)*(5*s - 1)*(2*r + s - 1)
    elseif j == 6; return sqrt(6)*(10*s^2 - 8*s + 1)
    elseif j == 7; return 2*sqrt(14)*(2*r + s - 1)*(10*r^2 + 10*r*(s - 1) + s^2 - 2*s + 1)
    elseif j == 8; return 2*sqrt(10)*(7*s - 1)*(6*r^2 + 6*r*(s - 1) + s^2 - 2*s + 1)
    elseif j == 9; return 2*sqrt(6)*(21*s^2 - 12*s + 1)*(2*r + s - 1)
    elseif j == 10; return 2*sqrt(2)*(35*s^3 - 45*s^2 + 15*s - 1)
    elseif j == 11; return 3*sqrt(10)*(70*r^4 + 140*r^3*(s - 1) + 90*r^2*(s^2 - 2*s + 1) + 20*r*(s^3 - 3*s^2 + 3*s - 1) + s^4 - 4*s^3 + 6*s^2 - 4*s + 1)
    elseif j == 12; return sqrt(70)*(9*s - 1)*(2*r + s - 1)*(10*r^2 + 10*r*(s - 1) + s^2 - 2*s + 1)
    elseif j == 13; return 5*sqrt(2)*(36*s^2 - 16*s + 1)*(6*r^2 + 6*r*(s - 1) + s^2 - 2*s + 1)
    elseif j == 14; return sqrt(30)*(84*s^3 - 84*s^2 + 21*s - 1)*(2*r + s - 1)
    elseif j == 15; return sqrt(10)*(126*s^4 - 224*s^3 + 126*s^2 - 24*s + 1)
    else; return eval_dubiner(r,s,j)
    end
end

"""
eval_Ddubiner_ex(z,w,j::Integer)
Compute the gradient of dubiner polynomial of degree `j` at point (r,s)
on the reference triangle ((0,0),(1,0),(0,1))
Use exact functions when available
"""
function eval_Ddubiner_ex(r::T,s::T,j::Integer) where{T}
    if j == 0; return [zero(T),zero(T)]
    elseif j == 1; return [zero(T),zero(T)]
    elseif j == 2; return [4*sqrt(3), 2*sqrt(3)]*one(T)
    elseif j == 3; return [zero(T), one(T)*6]
    elseif j == 4; return [12*sqrt(30)*r + 6*sqrt(30)*(s - 1),
                            6*sqrt(30)*r + 2*sqrt(30)*s - 2*sqrt(30)]
    elseif j == 5; return [6*sqrt(2)*(5*s - 1),
                            30*sqrt(2)*r + 30*sqrt(2)*s - 18*sqrt(2)]
    elseif j == 6; return [zero(T),20*sqrt(6)*s - 8*sqrt(6)]
    elseif j == 7; return [24*sqrt(14)*(5*r^2 + 5*r*(s - 1) + s^2 - 2*s + 1),
                            60*sqrt(14)*r^2 + r*(48*sqrt(14)*s - 48*sqrt(14)) + 6*sqrt(14)*s^2 - 12*sqrt(14)*s + 6*sqrt(14)]
    elseif j == 8; return [12*sqrt(10)*(7*s - 1)*(2*r + s - 1),
                            84*sqrt(10)*r^2 + r*(168*sqrt(10)*s - 96*sqrt(10)) + 42*sqrt(10)*s^2 - 60*sqrt(10)*s + 18*sqrt(10)]
    elseif j == 9; return [4*sqrt(6)*(21*s^2 - 12*s + 1),
                            r*(168*sqrt(6)*s - 48*sqrt(6)) + 126*sqrt(6)*s^2 - 132*sqrt(6)*s + 26*sqrt(6)]
    elseif j == 10; return [zero(T),210*sqrt(2)*s^2 - 180*sqrt(2)*s + 30*sqrt(2)]
    elseif j == 11; return [840*sqrt(10)*r^3 + 1260*sqrt(10)*r^2*(s - 1) + 540*sqrt(10)*r*(s^2 - 2*s + 1) + 60*sqrt(10)*(s^3 - 3*s^2 + 3*s - 1),
                            420*sqrt(10)*r^3 + r^2*(540*sqrt(10)*s - 540*sqrt(10)) + r*(180*sqrt(10)*s^2 - 360*sqrt(10)*s + 180*sqrt(10)) + 12*sqrt(10)*s^3 - 36*sqrt(10)*s^2 + 36*sqrt(10)*s - 12*sqrt(10)]
    elseif j == 12; return [12*sqrt(70)*(9*s - 1)*(5*r^2 + 5*r*(s - 1) + s^2 - 2*s + 1),
                            180*sqrt(70)*r^3 + r^2*(540*sqrt(70)*s - 300*sqrt(70)) + r*(324*sqrt(70)*s^2 - 456*sqrt(70)*s + 132*sqrt(70)) + 36*sqrt(70)*s^3 - 84*sqrt(70)*s^2 + 60*sqrt(70)*s - 12*sqrt(70)]
    elseif j == 13; return [30*sqrt(2)*(36*s^2 - 16*s + 1)*(2*r + s - 1),
                            r^2*(2160*sqrt(2)*s - 480*sqrt(2)) + r*(3240*sqrt(2)*s^2 - 3120*sqrt(2)*s + 510*sqrt(2)) + 720*sqrt(2)*s^3 - 1320*sqrt(2)*s^2 + 690*sqrt(2)*s - 90*sqrt(2)]
    elseif j == 14; return [2*sqrt(30)*(84*s^3 - 84*s^2 + 21*s - 1),
                            r*(504*sqrt(30)*s^2 - 336*sqrt(30)*s + 42*sqrt(30)) + 336*sqrt(30)*s^3 - 504*sqrt(30)*s^2 + 210*sqrt(30)*s - 22*sqrt(30)]
    elseif j == 15; return [zero(T),504*sqrt(10)*s^3 - 672*sqrt(10)*s^2 + 252*sqrt(10)*s - 24*sqrt(10)]
    else; return eval_Ddubiner(r,s,j)
    end
end

#test
const rr = 0.5
const ss = 0.5
for j = 1:15
    println(eval_dubiner(rr,ss,j) ≈ eval_dubiner_ex(rr,ss,j))
    println(eval_Ddubiner(rr,ss,j) ≈ eval_Ddubiner_ex(rr,ss,j))
end

max_stars_repo_name: Paulms/MyJuliaScripts

max_issues_repo_name: Paulms/MyJuliaScripts

max_forks_repo_name: Paulms/MyJuliaScripts


text:
import Flux: OneHotMatrix
import Mill: ArrayNode, MaybeHotMatrix
"""
	struct ExtractCategorical{V,I} <: AbstractExtractor
		keyvalemap::Dict{V,I}
		n::Int
		uniontypes::Bool
	end
	ExtractCategorical(s::Entry, uniontypes = true)
	ExtractCategorical(s::UnitRange, uniontypes = true)
	ExtractCategorical(s::Vector, uniontypes = true)

Converts a single item to a one-hot encoded vector.
Converts array of items into matrix of one-hot encoded columns.
There is always alocated an extra element for a unknown value.
If passed `missing`, if `uniontypes` is true, returns column of missing values, otherwise raises error.
If `uniontypes` is true, it allows extracting `missing` values and all extracted values will be of type
`Union{Missing, <other type>}` due to type stability reasons. Otherwise missings extraction is not allowed.

# Examples

```jldoctest
julia> using Mill: catobs

julia> e = ExtractCategorical(2:4, true);

julia> mapreduce(e, catobs, [2,3,1,4]).data
4×4 MaybeHotMatrix with eltype Union{Missing, Bool}:
  true    ⋅      ⋅      ⋅
   ⋅     true    ⋅      ⋅
   ⋅      ⋅      ⋅     true
   ⋅      ⋅     true    ⋅

julia> mapreduce(e, catobs, [1,missing,5]).data
4×3 MaybeHotMatrix with eltype Union{Missing, Bool}:
   ⋅    missing    ⋅
   ⋅    missing    ⋅
   ⋅    missing    ⋅
  true  missing   true

julia> e(4).data
4×1 MaybeHotMatrix with eltype Union{Missing, Bool}:
   ⋅
   ⋅
  true
   ⋅

julia> e(missing).data
4×1 MaybeHotMatrix with eltype Union{Missing, Bool}:
 missing
 missing
 missing
 missing

julia> e = ExtractCategorical(2:4, false);

julia> mapreduce(e, catobs, [2,3,1,4]).data
4×4 OneHotMatrix(::Vector{UInt32}) with eltype Bool:
 1  ⋅  ⋅  ⋅
 ⋅  1  ⋅  ⋅
 ⋅  ⋅  ⋅  1
 ⋅  ⋅  1  ⋅

julia> e(4).data
4×1 OneHotMatrix(::Vector{UInt32}) with eltype Bool:
 ⋅
 ⋅
 1
 ⋅
```
"""
struct ExtractCategorical{V,I} <: AbstractExtractor
	keyvalemap::Dict{V,I}
	n::UInt32
	uniontypes::Bool
end

ExtractCategorical(s::Entry, uniontypes = true) = ExtractCategorical(collect(keys(s.counts)), uniontypes)
ExtractCategorical(s::UnitRange, uniontypes = true) = ExtractCategorical(collect(s), uniontypes)
function ExtractCategorical(ks::Vector, uniontypes = true)
	if isempty(ks)
		@warn "Skipping initializing empty categorical variable does not make much sense to me"
		return nothing
	end
	ks = sort(unique(ks))
	keys_len = length(ks)
	if keys_len > typemax(UInt32)
		@error "Dictionary is too big, we support uint32 types"
	end
	# keeping int32 indices, because it's faster than int64 and int64 is too big to be usable anyway
	ExtractCategorical(Dict(zip(ks, UInt32(1):UInt32(keys_len))), UInt32(keys_len + 1), uniontypes)
end

map_val(s, ::MissingOrNothing) = s.uniontypes ? missing : error("This extractor does not support missing values")
map_val(s, v) = get(s.keyvalemap, v, s.n)
stabilize_types_categorical(s::ExtractCategorical{V,I}, x) where {V,I} = s.uniontypes ? Vector{Union{Missing, I}}(x) : x
val2idx(s::ExtractCategorical{V,I}, v::V) where {V,I} = stabilize_types_categorical(s, [map_val(s, v)])
val2idx(s::ExtractCategorical{<:Number,I}, v::Number) where {V,I} = stabilize_types_categorical(s, [map_val(s, v)])
construct(s::ExtractCategorical, x, y) = s.uniontypes ? MaybeHotMatrix(x, y) : OneHotMatrix(x, y)

make_missing_categorical(s::ExtractCategorical, v, store_input) =
    s.uniontypes ?
    _make_array_node(MaybeHotMatrix(stabilize_types_categorical(s, [missing]), s.n), [v], store_input) :
    error("This extractor does not support missing values")

(s::ExtractCategorical{V})(v::V; store_input=false) where {V<:HierarchicType} =
	_make_array_node(construct(s, val2idx(s, v), s.n), [v], store_input)
# I need to define this in order to disambiguate the method with different numbers compared to HierarchicType
(s::ExtractCategorical{V})(v::V; store_input=false) where {V<:Number} =
	_make_array_node(construct(s, val2idx(s, v), s.n), [v], store_input)

# following 2 methods are to let us extract float from int extractor and vice versa
(s::ExtractCategorical{<:Number})(v::Number; store_input=false) =
	_make_array_node(construct(s, val2idx(s, v), s.n), [v], store_input)

# following 2 methods are to let us extract numeric string from float or int extractor
# I'm trying to parse as float because integer can be parsed as float so I assume all numbers we care about
# are "floatable". Yes, this does not work for
(s::ExtractCategorical{<:Number})(v::AbstractString; store_input=false) =
	_make_array_node(construct(s, val2idx(s, tryparse(FloatType, v)), s.n), [v], store_input)
(s::ExtractCategorical)(v::MissingOrNothing; store_input=false) = make_missing_categorical(s, v, store_input)
(s::ExtractCategorical{V,I})(::ExtractEmpty; store_input=false) where {V,I} =
	ArrayNode(construct(s, s.uniontypes ? Vector{Union{Missing, I}}() : Vector{I}(), s.n))

# todo: this is more specific than V, V method and we need to fix it
(s::ExtractCategorical)(v::HierarchicType; store_input=false) = make_missing_categorical(s, v, store_input)

Base.hash(e::ExtractCategorical, h::UInt) = hash((e.keyvalemap, e.n, e.uniontypes), h)
Base.:(==)(e1::ExtractCategorical, e2::ExtractCategorical) =
	e1.keyvalemap == e2.keyvalemap && e1.n === e2.n && e1.uniontypes === e2.uniontypes


text:
module BaseOnnx

using ProtoBuf
include("onnx_pb.jl")
include("read.jl")
include("write.jl")

read(io) = readproto(io, ModelProto())

end


text:
using Convolution2D
using Test

@testset "Convolution2D" begin
    X = [
        1  4  5  4  4
        1  3  4  5  3
        2  3  5  4  4
        4  4  2  1  1
        3  3  3  2  4
    ]
    kern = [0 1 0
		1 -4 1
		0 1 0]
    out = conv2D(X, kern)
    @test out[3, 3] == -7

    X =  [5  2  5  3
        1  2  2  4
        1  1  1  4
        3  2  4  3]
    @test_broken conv2D(X, kern)
end


text:
reload("RMB")

module dev

using RMB
import RMB: inventorymodel, pidc1, pidc2, detensor, rank, pi!, rsm

dims=(5,3,10)
m2=inventorymodel(dims,5.0,20.0,[0.8,1.0,1.2])
pi!(m2,pidc1(m2))
u!(m2,[0.0,1.0])    #free stocking, expensive order
solvedp!(m2)

dd,dc,ds,dy,dr=m2.dims
dk,da=dd*dc*ds,dy*dr

ac=reshape([1.0*(ik==jk)*(ia==1) for ik=1:dk,jk=1:dk,ia=1:da],dk,dk*da)

function pi!(model::InventoryModel,eta::AbstractVector)
    pidc=RMB.z2q(eta)
    pi!(model,pidc)
end

eta=rand(dd*dc*(dd*dc-1))
j1=RMB.idrank(m2,ac,eta)

# pic1=rsm(dc,dc)
# pic2=rsm(dc,dc)
# pid1=rsm(dd,dd)
# pid2=rsm(dd,dd)

# function pi!(model::InventoryModel,eta::AbstractVector)
#     pic=eta[1]*pic1+(1-eta[1])*pic2
#     pid=eta[2]*pid1+(1-eta[2])*pid2
#     pi!(model,kron(pic,pid))
# end

# eta=[0.5,0.5]
# j2=RMB.idrank(m2,eta)

end #module


text:
using Genie.Router, Genie.Renderer.Html, Genie.Renderer.Json
import Genie.Renderer.Json.json

route("/") do
  serve_static_file("welcome.html")
end

if haskey(ENV, "TURN_ON_LOGGING")
  const GREETING = "Welcome to Genie.jl! This site is served by Julia $(Base.VERSION) with $(Threads.nthreads()) threads. Logging is turned on."
else
  const GREETING = "Welcome to Genie.jl! This site is served by Julia $(Base.VERSION) with $(Threads.nthreads()) threads. Logging is turned off."
end

route("/greeting", named = :greeting) do
  html(GREETING)
end

route("/random", named = :random) do
  p = params(:num, "1")
  num = tryparse(Int32, p)
  if (num === nothing || num < 1)
      num = 1
  end

  json(rand(Int, num))
end

max_stars_repo_name: StatisticalMice/performance-tests

max_issues_repo_name: StatisticalMice/performance-tests

max_forks_repo_name: StatisticalMice/performance-tests


text:
using ChangePointMean
using Test

@testset "ChangePointMean.jl" begin
    s0 = [1, 2, 1, 2, 1, 1, 1, 2, 2, 1, 2, 1, 2, 2, 2, 1, 1, 2, 1]
    s1 = [10, 11, 9, 9, 10, 10, 11, 11, 18, 16, 17, 17, 17, 16, 18]
    s2 = [10, 11, 9, 9, 10, 10, 11, 11, 18, 16, 17, 17, 17, 16, 18, 21, 22, 22, 23, 22, 21, 21, 23]

    @test mcpoint(s0) == 0
    @test mcplast(s0) == 0
    @test mcpall(s0) == Int[]

    @test mcpoint(s1) == 9
    @test mcplast(s1) == 9
    @test mcpall(s1) == [9]

    @test mcpoint(s2) == 9
    @test mcplast(s2) == 16
    @test mcpall(s2) == [9, 16]
end


text:
# Julia wrapper for header: /usr/local/include/libavformat/avio.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    avio_check,
    avio_alloc_context,
    avio_w8,
    avio_write,
    avio_wl64,
    avio_wb64,
    avio_wl32,
    avio_wb32,
    avio_wl24,
    avio_wb24,
    avio_wl16,
    avio_wb16,
    avio_put_str,
    avio_put_str16le,
    avio_seek,
    avio_skip,
    avio_tell,
    avio_size,
    avio_flush,
    avio_read,
    avio_r8,
    avio_rl16,
    avio_rl24,
    avio_rl32,
    avio_rl64,
    avio_rb16,
    avio_rb24,
    avio_rb32,
    avio_rb64,
    avio_get_str,
    avio_get_str16le,
    avio_get_str16be,
    avio_open,
    avio_open2,
    avio_close,
    avio_closep,
    avio_open_dyn_buf,
    avio_close_dyn_buf,
    avio_enum_protocols,
    avio_pause,
    avio_seek_time

function avio_check(url,flags::Integer)
    ccall((:avio_check,libavformat),Cint,(Ptr{UInt8},Cint),url,flags)
end

function avio_alloc_context(buffer,buffer_size::Integer,write_flag::Integer,opaque,read_packet,write_packet,seek)
    ccall((:avio_alloc_context,libavformat),Ptr{AVIOContext},(Ptr{Cuchar},Cint,Cint,Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void}),buffer,buffer_size,write_flag,opaque,read_packet,write_packet,seek)
end

function avio_w8(s,b::Integer)
    ccall((:avio_w8,libavformat),Void,(Ptr{AVIOContext},Cint),s,b)
end

function avio_write(s,buf,size::Integer)
    ccall((:avio_write,libavformat),Void,(Ptr{AVIOContext},Ptr{Cuchar},Cint),s,buf,size)
end

function avio_wl64(s,val::UInt64)
    ccall((:avio_wl64,libavformat),Void,(Ptr{AVIOContext},UInt64),s,val)
end

function avio_wb64(s,val::UInt64)
    ccall((:avio_wb64,libavformat),Void,(Ptr{AVIOContext},UInt64),s,val)
end

function avio_wl32(s,val::Integer)
    ccall((:avio_wl32,libavformat),Void,(Ptr{AVIOContext},UInt32),s,val)
end

function avio_wb32(s,val::Integer)
    ccall((:avio_wb32,libavformat),Void,(Ptr{AVIOContext},UInt32),s,val)
end

function avio_wl24(s,val::Integer)
    ccall((:avio_wl24,libavformat),Void,(Ptr{AVIOContext},UInt32),s,val)
end

function avio_wb24(s,val::Integer)
    ccall((:avio_wb24,libavformat),Void,(Ptr{AVIOContext},UInt32),s,val)
end

function avio_wl16(s,val::Integer)
    ccall((:avio_wl16,libavformat),Void,(Ptr{AVIOContext},UInt32),s,val)
end

function avio_wb16(s,val::Integer)
    ccall((:avio_wb16,libavformat),Void,(Ptr{AVIOContext},UInt32),s,val)
end

function avio_put_str(s,str)
    ccall((:avio_put_str,libavformat),Cint,(Ptr{AVIOContext},Ptr{UInt8}),s,str)
end

function avio_put_str16le(s,str)
    ccall((:avio_put_str16le,libavformat),Cint,(Ptr{AVIOContext},Ptr{UInt8}),s,str)
end

function avio_seek(s,offset::Int64,whence::Integer)
    ccall((:avio_seek,libavformat),Int64,(Ptr{AVIOContext},Int64,Cint),s,offset,whence)
end

function avio_skip(s,offset::Int64)
    ccall((:avio_skip,libavformat),Int64,(Ptr{AVIOContext},Int64),s,offset)
end

function avio_tell(s)
    ccall((:avio_tell,libavformat),Int64,(Ptr{AVIOContext},),s)
end

function avio_size(s)
    ccall((:avio_size,libavformat),Int64,(Ptr{AVIOContext},),s)
end

function avio_flush(s)
    ccall((:avio_flush,libavformat),Void,(Ptr{AVIOContext},),s)
end

function avio_read(s,buf,size::Integer)
    ccall((:avio_read,libavformat),Cint,(Ptr{AVIOContext},Ptr{Cuchar},Cint),s,buf,size)
end

function avio_r8(s)
    ccall((:avio_r8,libavformat),Cint,(Ptr{AVIOContext},),s)
end

function avio_rl16(s)
    ccall((:avio_rl16,libavformat),UInt32,(Ptr{AVIOContext},),s)
end

function avio_rl24(s)
    ccall((:avio_rl24,libavformat),UInt32,(Ptr{AVIOContext},),s)
end

function avio_rl32(s)
    ccall((:avio_rl32,libavformat),UInt32,(Ptr{AVIOContext},),s)
end

function avio_rl64(s)
    ccall((:avio_rl64,libavformat),UInt64,(Ptr{AVIOContext},),s)
end

function avio_rb16(s)
    ccall((:avio_rb16,libavformat),UInt32,(Ptr{AVIOContext},),s)
end

function avio_rb24(s)
    ccall((:avio_rb24,libavformat),UInt32,(Ptr{AVIOContext},),s)
end

function avio_rb32(s)
    ccall((:avio_rb32,libavformat),UInt32,(Ptr{AVIOContext},),s)
end

function avio_rb64(s)
    ccall((:avio_rb64,libavformat),UInt64,(Ptr{AVIOContext},),s)
end

function avio_get_str(pb,maxlen::Integer,buf,buflen::Integer)
    ccall((:avio_get_str,libavformat),Cint,(Ptr{AVIOContext},Cint,Ptr{UInt8},Cint),pb,maxlen,buf,buflen)
end

function avio_get_str16le(pb,maxlen::Integer,buf,buflen::Integer)
    ccall((:avio_get_str16le,libavformat),Cint,(Ptr{AVIOContext},Cint,Ptr{UInt8},Cint),pb,maxlen,buf,buflen)
end

function avio_get_str16be(pb,maxlen::Integer,buf,buflen::Integer)
    ccall((:avio_get_str16be,libavformat),Cint,(Ptr{AVIOContext},Cint,Ptr{UInt8},Cint),pb,maxlen,buf,buflen)
end

function avio_open(s,url,flags::Integer)
    ccall((:avio_open,libavformat),Cint,(Ptr{Ptr{AVIOContext}},Ptr{UInt8},Cint),s,url,flags)
end

function avio_open2(s,url,flags::Integer,int_cb,options)
    ccall((:avio_open2,libavformat),Cint,(Ptr{Ptr{AVIOContext}},Ptr{UInt8},Cint,Ptr{AVIOInterruptCB},Ptr{Ptr{AVDictionary}}),s,url,flags,int_cb,options)
end

function avio_close(s)
    ccall((:avio_close,libavformat),Cint,(Ptr{AVIOContext},),s)
end

function avio_closep(s)
    ccall((:avio_closep,libavformat),Cint,(Ptr{Ptr{AVIOContext}},),s)
end

function avio_open_dyn_buf(s)
    ccall((:avio_open_dyn_buf,libavformat),Cint,(Ptr{Ptr{AVIOContext}},),s)
end

function avio_close_dyn_buf(s,pbuffer)
    ccall((:avio_close_dyn_buf,libavformat),Cint,(Ptr{AVIOContext},Ptr{Ptr{UInt8}}),s,pbuffer)
end

function avio_enum_protocols(opaque,output::Integer)
    ccall((:avio_enum_protocols,libavformat),Ptr{UInt8},(Ptr{Ptr{Void}},Cint),opaque,output)
end

function avio_pause(h,pause::Integer)
    ccall((:avio_pause,libavformat),Cint,(Ptr{AVIOContext},Cint),h,pause)
end

function avio_seek_time(h,stream_index::Integer,timestamp::Int64,flags::Integer)
    ccall((:avio_seek_time,libavformat),Int64,(Ptr{AVIOContext},Cint,Int64,Cint),h,stream_index,timestamp,flags)
end

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
ackno(tr::TCPReceiver)::Union{Nothing,WrappingInt32} =
    !tr.has_syn ? nothing : wrap(UInt(stream_out(tr).bytes_written_count) + stream_out(tr).has_eof + 1, tr.isn)

window_size(tr::TCPReceiver) = tr.capacity - buffer_size(tr.reassembler |> stream_out)

unassembled_bytes(tr::TCPReceiver) = tr.reassembler.unassembled_bytes

stream_out(tr::TCPReceiver) = tr.reassembler |> stream_out

assembled_bytes(tr::TCPReceiver) = tr.reassembler.cur_idx

function segment_received!(tr::TCPReceiver, seg::TCPSegment)::Nothing
    header = seg.header
    paload = seg.payload
    if (header.syn)
        tr.has_syn && return
        tr.has_syn = true
        tr.isn = header.seqno
        tr.abs_seqno = 1
        (length_in_sequence_space(seg) == 1) && return
    elseif !tr.has_syn
        return
    else
        tr.abs_seqno = unwrap(header.seqno, tr.isn, tr.abs_seqno)
    end

    (length_in_sequence_space(seg) == 0) && return

    push_substring!(tr.reassembler, paload.storage[], tr.abs_seqno - 1, header.fin)
end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
using PyPlot, ProgressMeter;

output_dir = "pic"
if ! ispath(output_dir); mkdir(output_dir); end

function pic()
    L = 5.0 # system size
    steps = 500;
    progress = Progress(steps);
    for i in 0:steps-1
        filename = @sprintf("./data/%05d.csv", i);
        data = readcsv(filename);
        x = data[:,1]; y = data[:,2]; vx = data[:,3]; vy = data[:,4];
        plot(x,y,"b.");
        quiver(x,y,vx,vy);
        axis("square"); axis([0, L, 0, L]);
        picname = @sprintf("./pic/%05d.png", i);
        savefig(picname, dpi=60);

        clf();
        next!(progress);
    end
end
pic()
gc();

max_stars_repo_name: goropikari/vicsek-model

max_issues_repo_name: goropikari/vicsek-model

max_forks_repo_name: goropikari/vicsek-model


text:
function plot_coda(sampler ; burn = 0, thin = 1, id_par = [1:4...], smooth = false)
    chains = sampler[1]
    @rput chains
    @rput thin
    @rput burn
    @rput id_par
    @rput smooth

    reval("""
        require(coda)
        N = length(chains)
        n_row = length(chains[[1]][[1]])
        ini = floor(burn * n_row)
        if(ini == 0) ini = 1
            it = seq(ini, n_row, thin)

        plot_chains = mcmc.list(lapply(1:N, function(i){
        mcmc(sapply((id_par), function(j){
        if(j<4) {
            chains[[i]][[j]][it]
        }else{
            log(chains[[i]][[j]][it])
        }
        }))}))

        plot(plot_chains, smooth = smooth)
    """)
end

function add_to_chain!(chain, param)
    for i in 1:length(param)
        push!(chain[i], param[i])
    end
end

function initialize_ClaDS2(tree::Tree ; ini_par = [], initialize_rates = 0, ltt_steps = 10, n_chains = 3, n_trees = 10) where {N}

    new_tree = Tree(deepcopy(tree.offsprings), 0., deepcopy(tree.attributes))

    root_depth = maximum(node_depths(new_tree)) * 1.01
    live_nd = node_depths_base(new_tree)

    ltt_times = [0:ltt_steps...] * root_depth/ltt_steps
    ltt_extant = LTT(tree,ltt_times)
    if length(ini_par) == 0
        ini_par = [[[10. ^(1-i),1,0.5];fill(0.000001, tree.n_nodes); fill(0.000001, Int64((tree.n_nodes+1)/2));[0, n_extant_tips(tree)] ; ltt_extant[2]] for i in 1:n_chains]
    end

    edge_trees_s = [init_edge_tree_rates2(update_rates(new_tree, ini_par[i][4:(tree.n_nodes + 3)]), ini_par[i][4:(end-2)]) for i in 1:n_chains]

    for i in 1:n_chains
        new_tree = Tree(deepcopy(tree.offsprings), 0., deepcopy(tree.attributes))
        update_rates!(new_tree, ini_par[i][4:(tree.n_nodes + 3)])

        rates = ini_par[i][4:(tree.n_nodes + 3)]
        ε = ini_par[i][3]
        σ = ini_par[i][1]
        α = ini_par[i][2]
        edge_trees = edge_trees_s[i]

        ini_par[i][4:(tree.n_nodes + 3)] = rates
        ini_par[i][3] = ε
        edge_trees_s[i] = edge_trees
    end

    extant_branch_lengths = extract_branch_lengths(tree)
    trees = [deepcopy(update_rates(new_tree, ini_par[i][4:end])) for i in 1:n_chains]
    chains = [[[deepcopy(ini_par[i][j])] for j in 1:length(ini_par[1])] for i in 1:n_chains]
    param = deepcopy(ini_par)
    mean_rates_chains = [[time_rates(tree,ini_par[i][4:(tree.n_nodes+3)],ltt_times)] for i in 1:n_chains]
    println(" ")
    return chains, param, edge_trees_s, trees, extant_branch_lengths, ltt_times, live_nd, mean_rates_chains, Array{Tree,1}(undef,0)
end

function add_iter_ClaDS2(sampler, n_reccord::Int64; thin = 1, fs = 1., plot_tree = 0, print_state = 0,
    max_node_number = 1_000, max_try = 100_000, it_edge_tree = 1, print_all = false, it_rates = 1, n_trees = 10, ini = false, prior_ε = "uniform", logε0 = 0., sdε = 0.5)

    chain_s, param_s, edge_trees_s, tree_s, extant_branch_lengths, ltt_times, live_nd, mean_rates_chains, enhanced_trees = sampler

    n_chains = length(chain_s)
    tips_id = tips(tree_s[1])[2:end]
    lefts = n_left(tree_s[1])
    ntips = Int64((tree_s[1].n_nodes+1)/2)

    if plot_tree>0
        reval("par(mfrow=c(3,4), mar=c(5,2,2,2))")
    end

    n_par = tree_s[1].n_nodes + 3

    for k in 1:n_chains
        chain = chain_s[k]
        tree = tree_s[k]

        σ = param_s[k][1]
        α = param_s[k][2]
        ε = param_s[k][3]
        rates = param_s[k][4:(tree.n_nodes + 3)]

        edge_trees = edge_trees_s[k]

        relative_rates = Array{Float64,1}(undef,0)

        for i in 1:n_reccord
            for j in 1:thin
                for l in 1:it_edge_tree
                    update_edge_trees!(edge_trees, tree, σ, α, ε, fs, rates, extant_branch_lengths,
                        max_node_number = max_node_number, max_try = max_try, keep_if_any = true,
                        do_tips =  true)

                    relative_rates = extract_relative_rates(tree, edge_trees, rates)
                    draw_λ0_slicing!(rates, edge_trees, σ, α, ε, lefts)
                    σ = draw_σ(relative_rates, α, β0 = 0.05, α0 = 0.5)
                    α = draw_α(relative_rates, σ,  α_0 = -0.05, σ_0 = 0.1)
                    if prior_ε == "uniform"
                        ε = draw_ε_crown_priorUnif(tree, edge_trees, lefts)
                    elseif prior_ε == "uniformInf"
                        ε = draw_ε_crown(tree, edge_trees, lefts)
                    elseif prior_ε == "ClaDS0"
                        ε = 0.
                    else
                        ε = draw_ε_crown_priorln(tree, edge_trees, lefts, logε0 = logε0, sd = sdε)
                    end
                    draw_λ0!(tree::Tree, ε::Float64, rates::Array{Float64,1},
                        edge_trees::Array{EdgeTreeRates2,1})
                    tree.attributes[1] = rates[1]
                end

                if ini
                    print(".")
                    draw_λ0_slicing!(rates, edge_trees, σ, α, ε, lefts)
                    draw_λ0!(tree::Tree, ε::Float64, rates::Array{Float64,1},
                        edge_trees::Array{EdgeTreeRates2,1})
                    tree.attributes[1] = rates[1]
                end

                if it_rates > 0
                    ε, σ, α = update_edges_ETR2!(tree, edge_trees, σ, α, ε, rates,lefts,it_rates = it_rates,  prior_ε = prior_ε, logε0 = logε0, sd = sdε)
                end

                if plot_tree>0
                    if i%plot_tree == 0 && j == thin
                        enhanced_tree = graft_edge_trees(tree, edge_trees)
                        plot_ClaDS(enhanced_tree)#, extant_edges[2:end], ln = false)
                    end
                end
                if print_state>0
                    if i%print_state == 0 && (j == thin || print_all)
                        relative_rates = extract_relative_rates(tree, edge_trees)
                        println("$i : sd $(sqrt(var(relative_rates))), mean $(mean(relative_rates)), σ $σ, α $α, ε $ε, λ0 $(rates[1]), n_ext $(extract_nextinct(tree,edge_trees)) $(n_tip(tree,edge_trees)) $(n_extant_tips(tree,edge_trees))")
                    end
                end

                if length(enhanced_trees) < n_trees
                    enhanced_tree = graft_edge_trees(tree, edge_trees)
                    push!(enhanced_trees, enhanced_tree)
                else
                    u = rand()
                    if u < 0.1*n_trees/i
                        enhanced_tree = graft_edge_trees(tree, edge_trees)
                        itree = sample(1:n_trees)
                        enhanced_trees[itree] = enhanced_tree
                    end
                end

                #println(" 3 nnn")

            end
            #println("$σ ")

            tip_rates = extract_tip_rates(tree, edge_trees, tips_id, rates)
            param_s[k][1] = σ
            param_s[k][2] = α
            param_s[k][3] = ε
            param_s[k][4:n_par] = rates
            param_s[k][(n_par+1):(n_par+ntips)] = tip_rates
            param_s[k][n_par + ntips + 1] = extract_nextinct(tree,edge_trees)
            param_s[k][n_par + ntips + 2] = n_tip(tree,edge_trees) - extract_nextinct(tree,edge_trees)
            #enhanced_tree = graft_edge_trees(tree, edge_trees)
            #ltt = LTT(enhanced_tree, ltt_times)[2]
            ltt = LTT(tree, edge_trees, ltt_times)[2]
            param_s[k][(n_par + ntips + 3):end] = ltt
            edge_trees_s[k] = edge_trees
            add_to_chain!(chain_s[k], param_s[k])
            #println(time_rates(tree,edge_trees,ltt_times))
            #println(mean_rates_chains)
            for imr in 1
                #add_to_chain!(mean_rates_chains[k], time_rates(tree,edge_trees,ltt_times))
                push!(mean_rates_chains[k], time_rates(tree,edge_trees,ltt_times))
            end
            tree_s[k] = tree
        end
    end

    println(" ")
    return chain_s, param_s, edge_trees_s, tree_s, extant_branch_lengths, ltt_times, live_nd, mean_rates_chains, enhanced_trees
end

function gelman_est(mcmc::Array{Array{Array{Float64,1},1},1}, Npar::Int64 ; thin = 1, burn = 0)
    nChains = length(mcmc)
    Nit = length(mcmc[1][1])
    ini = Int64(floor(1 + Nit * burn))

    it = [ini]
    i = ini + thin
    while i<= Nit
        push!(it,i)
        i += thin
    end
    Nit = length(it)

    # compute means
    M = [zeros(Npar) for i in 1:nChains]

    for n in 1:nChains
        for i in 1:Npar
            if ((i == 1 )|| (i == 3))
                for j in it
                    M[n][i] += sum(mcmc[n][i][j])
                end
            else
                for j in it
                    M[n][i] += sum(log(mcmc[n][i][j]))
                end
            end
        end
    end

    M ./= Nit

    m = zeros(Npar)
    for n in 1:nChains
        m .+= M[n]
    end

    m ./= nChains

    # compute vars
    Σ = [zeros(Npar) for i in 1:nChains]

    for n in 1:nChains
        for i in 1:Npar
            if ((i == 1 )|| (i == 3))
                for j in it
                    Σ[n][i] += (mcmc[n][i][j] - M[n][i])^2
                end
            else
                for j in it
                    Σ[n][i] += (log(mcmc[n][i][j]) - M[n][i])^2
                end
            end

        end
    end

    Σ ./= (Nit-1)

    # between chains variance
    B = zeros(Npar)

    for n in 1:nChains
        for i in 1:Npar
            B[i] += (M[n][i] - m[i])^2
        end
    end

    B .* (nChains + 1)/(nChains * (nChains - 1))

    # within chains variance

    W = zeros(Npar)
    for n in 1:nChains
        W .+= Σ[n]
    end
    W ./= nChains

    Wd = W .* (Nit-1)/Nit

    # pooled variance
    V = Wd .+ B

    # PSRF
    PSRF = V ./ W

    #return PSRF
    id = 0
    maxi = 0.

    for i in 1:Npar
        if PSRF[i] > maxi
            id = i
            maxi = PSRF[i]
        end
    end

    return id, sqrt(maxi)
end

function get_parent_rate(tree::Tree, edge_id::Int64, edge_trees::Array{EdgeTreeRates2,1})
    parent_edge = edge_trees[edge_id].parent_edge

    if parent_edge == 0
        return tree.attributes[1]
    else
        return edge_trees[parent_edge].tip_rate * edge_trees[parent_edge].stem_rate[1]
    end
end

function get_parent_rate(tree::Tree, edge_id::Int64, edge_trees::Array{EdgeTreeRates2,1}, rates::Array{Float64,1})
    parent_edge = edge_trees[edge_id].parent_edge

    if parent_edge == 0
        return rates[1]
    else
        return edge_trees[parent_edge].tip_rate * rates[parent_edge+1]
    end
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
export domain, accuracy, spine

"""
    domain(A::PolygonalAmoeba)

The domain `(xmin, xmax, ymin, ymax)` for which the amoeba is defined.
"""
domain(A::PolygonalAmoeba) = A.domain

"""
    accuracy(A::PolygonalAmoeba)

The accuracy  of the approximation.
"""
accuracy(A::PolygonalAmoeba) = A.error_bound

"""
    spine(A::PolygonalAmoeba)

The spine which was used to compute `A`.
"""
spine(A::PolygonalAmoeba) = A.spine

function Base.show(io::IO, A::PolygonalAmoeba)
    println(io, "PolygonalAmoeba:")
    println(io, " * accuracy: ", accuracy(A))
end

function Base.show(io::IO, mime::MIME"text/html", A::PolygonalAmoeba)
    if inIJulia() && plotsdefined()
        show(io, mime, Main.Plots.plot(A))
    else
        show(io, A)
    end
end

@recipe function plot(A::PolygonalAmoeba; spine=false)
    xmin, xmax, ymin, ymax = A.domain
    xlims --> (xmin, xmax)
    ylims --> (ymin, ymax)
    aspect_ratio --> :equal
    for (p1, p2, p3, p4) in A.polygons
        @series begin
            fillalpha--> get(plotattributes, :seriesalpha, 1.0)
            linecolor --> :dodgerblue
            fillcolor --> :dodgerblue
            linewidth --> 1.0
            linealpha --> get(plotattributes, :seriesalpha, get(plotattributes, :fillalpha, 1.0))
            legend --> false
            st := :shape
            x1, y1 = p1
            x2, y2 = p2
            x3, y3 = p3
            x4, y4 = p4
            x5, y5 = p1

            [x1, x2, x3, x4, x5], [y1, y2, y3, y4, y5]
        end
    end
    if spine
        @series begin
            A.spine
        end
    end
end


text:
"""
    genr_cell(cell_vec)

Do this will return a variable in `Data_Cell` type

# Example
```julia-repl
data_cell = genr_cell([1 2 3])
```
"""
function genr_cell(cell_vec)
    num_cells = 1
    for i in cell_vec
        num_cells *= i
    end

    cell_mat = Matrix{Int64}(undef, (num_cells, 3))
    for x in 0:cell_vec[1] - 1
        for y in 0:cell_vec[2] - 1
            for z in 0:cell_vec[3] - 1
                id = x * cell_vec[2] * cell_vec[3] + y * cell_vec[3] + z + 1
                cell_mat[id, :] = [x y z]
            end
        end
    end
    if typeof(cell_vec) == Array{Int64,2}
        cell_vec = conv(cell_vec, 1)
    end
    Data_Cell(cell_mat, cell_vec, num_cells)
end

# genr
"""
    genr(data_cell::Data_Cell, str::Str)

Do this will return a variable in `Data_Unit` type, containing all infomation needed for model of `str`.

# Example
```julia-repl
data_cell = genr_cell([5, 5, 5])
str = Tip3p()
data = genr(data_cell, str)
```
"""
function genr(data_cell::Data_Cell, str::Str)
    data = genr_atom(data_cell, str)
    str_fields = fieldnames(typeof(str))

    vec_bond = 0
    vec_angle = 0
    vec_dihedral = 0
    vec_improper = 0

    if in(:bond_mode, str_fields)
        data_bond = genr_bond(data_cell, data)
        vec_bond = data_bond.vec_bond
    end

    if in(:angle_mode, str_fields)
        data_angle = genr_angle(data_cell, data)
        vec_angle = data_angle.vec_angle
    end

    if in(:dihedral_mode, str_fields)
        data_dihedral = genr_dihedral(data_cell, data)
        vec_dihedral= data_dihedral.vec_dihedral
    end

    if in(:improper_mode, str_fields)
        data_improper = genr_improper(data_cell, data)
        vec_improper = data_improper.vec_improper
    end

    Data_Unit(data.data_basic, data.data_str, data.vec_atom,
        vec_bond, vec_angle, vec_dihedral, vec_improper)
end

# genr_atom
"""
    genr_atom(data_cell::Data_Cell, str::Str)

This will be called by `genr(...)`, genrating atomic information
"""
function genr_atom(data_cell::Data_Cell, str::Str)
    # Variables Setting
    vec_atom = Vector{Atom}(undef, data_cell.num_cells * str.num_atoms)

    # Generate Atom Data
    for cell in 1:data_cell.num_cells
        for atom in 1:str.num_atoms
            id_now = (cell - 1) * str.num_atoms + atom
            coord = conv((data_cell.cell_mat[cell, :] + str.atom_vec[atom, :]), 1) * str.cell_vec
            coord = conv(coord, 1) # Reduce Matrix to Vector
            vec_atom[id_now] = Atom(id_now, cell, str.atom_type[atom], str.atom_charge[atom], coord)
        end
    end

    # Generate Basic Data
    num_atoms = data_cell.num_cells * str.num_atoms
    num_bonds = 0
    num_angles = 0
    num_dihedrals = 0
    num_impropers = 0
    num_atom_types = str.num_atom_types
    num_bond_types = 0
    num_angle_types = 0
    num_dihedral_types = 0
    num_improper_types = 0
    box_size = Matrix{Float64}(undef, (3, 2))
    str_vec = Matrix{Float64}(undef, (3, 3))
    str_vec[:, :] = str.cell_vec[:, :]
    str_vec[2,1] = 0
    str_vec[3,1] = 0
    str_vec[3,2] = 0
    box_size[:, 1] = min(vec_atom, "coord")
    box_size[:, 2] = conv(conv(min(vec_atom, "coord"), 1) + conv(data_cell.cell_vec, 1) * str_vec, 1)
    box_tilt = Vector{Float64}(undef, 3)
    box_tilt[1] = data_cell.cell_vec[2] * str.cell_vec[2,1]
    box_tilt[2] = data_cell.cell_vec[3] * str.cell_vec[3,1]
    box_tilt[3] = data_cell.cell_vec[2] * str.cell_vec[3,2]

    data_basic = Data_Basic(num_atoms, num_bonds, num_angles, num_dihedrals, num_impropers,
            num_atom_types, num_bond_types, num_angle_types, num_dihedral_types, num_improper_types,
            box_size, box_tilt)

    # Output
    Data_Unit(data_basic, str, vec_atom, 0, 0, 0, 0)
end

# genr_bond
"""
    genr_bond!(data::Data; bond_radius=3.5)

This will genrate the bond information for `data` based on the geomarty judgement, bond will be genrated only if distance between two atoms with in `bond_radius`.

Notice:
- PBC issues has been considered.
- only structures with orthogonal unit cell are fine to use this function.
"""
function genr_bond!(data::Data; bond_radius=3.5)
	coord = get_data(data.vec_atom, "coord")
	num_atoms = data.data_basic.num_atoms
	box_diag = genr_box_diag(data)
	box_inv = inv(box_diag)
	coord_scl = coord * box_inv
	vec_bond = Vector{Bond}(undef, 1)

	for atom_01 = 1:num_atoms
		for atom_02 = atom_01 + 1:num_atoms
			r_diff = coord_scl[atom_02, :] - coord_scl[atom_01, :]
			r_diff -= round.(r_diff)
			r_diff = box_diag * r_diff
			r_diff = norm(r_diff)
			if r_diff <= bond_radius
				vec_bond = vcat(vec_bond, Bond(length(vec_bond), 1, [atom_01, atom_02]))
			end
		end
    end
    num_bonds = length(vec_bond) - 1
    if num_bonds == 0
        @printf("No bond has been generated within cutoff radius of %2.2f Angstrom", bond_radius)
        return 0
    end
	data.vec_bond = vec_bond[2:end]
	data.data_basic.num_bonds += num_bonds
    data.data_basic.num_bond_types += 1
    @printf("%4d bonds has been generated within cutoff radius of %2.2f Angstrom", data.data_basic.num_bonds, bond_radius)
end

"""
    genr_bond(data_cell::Data_Cell, data::Data)

This will be called by `genr(...)` if `bond_mode` is contained in `str`, genrating information of bonds
"""
function genr_bond(data_cell::Data_Cell, data::Data)
    # Reading Input
    data_basic = data.data_basic
    data_str = data.data_str
    num_cells = data_cell.num_cells
    num_cell_bonds = data_str.num_bonds
    num_bonds = num_cells * num_cell_bonds
    num_bond_types = data_str.num_bond_types

    data_basic.num_bonds = num_bonds
    data_basic.num_bond_types = num_bond_types

    # Generate Bonds
    vec_bond = Vector{Bond}(undef, num_bonds)
    for cell in 1:num_cells
        for bond = 1:num_cell_bonds
            id_now = (cell - 1) * num_cell_bonds + bond
            bond_now = Bond(data_str.bond_mode[bond])
            bond_now.id = id_now
            atom_tilt = (cell - 1) * data_str.num_atoms
            bond_now.atom .+= atom_tilt
            vec_bond[id_now] = bond_now
        end
    end

    # Output
    Data_Unit(data_basic, data_str, data.vec_atom, vec_bond, 0, 0, 0)
end

# genr_angle
"""
    genr_angle(data_cell::Data_Cell, data::Data)

This will be called by `genr(...)` if `angle_mode` is contained in `str`,
genrating information of angles
"""
function genr_angle(data_cell::Data_Cell, data::Data)
    # Reading Input
    data_basic = data.data_basic
    data_str = data.data_str
    num_cells = data_cell.num_cells
    num_cell_angles = data_str.num_angles
    num_angles = num_cells * num_cell_angles
    num_angle_types = data_str.num_angle_types

    data_basic.num_angles = num_angles
    data_basic.num_angle_types = num_angle_types

    # Generate Angle
    vec_angle = Vector{Angle}(undef, num_angles)
    for cell in 1:num_cells
        for angle in 1:num_cell_angles
            id_now = (cell - 1) * num_cell_angles + angle
            angle_now = Angle(data_str.angle_mode[angle])
            angle_now.id = id_now
            atom_tilt = (cell - 1) * data_str.num_atoms
            angle_now.atom .+= atom_tilt
            vec_angle[id_now] = angle_now
        end
    end

    # Output
    Data_Unit(data_basic, data_str, data.vec_atom,
            data.vec_bond, vec_angle, 0, 0)
end

# genr_dihedral
"""
    genr_dihedral(data_cell::Data_Cell, data::Data)

This will be called by `genr(...)` if `dihedral_mode` is contained in `str`,
genrating information of dihedral
"""
function genr_dihedral(data_cell::Data_Cell, data::Data)
    # Reading Input
    data_basic = data.data_basic
    data_str = data.data_str
    num_cells = data_cell.num_cells
    num_cell_dihedrals = data_str.num_dihedrals
    num_dihedrals = num_cells * num_cell_dihedrals
    num_dihedral_types = data_str.num_dihedral_types

    data_basic.num_dihedrals = num_dihedrals
    data_basic.num_dihedral_types = num_dihedral_types

    # Generate Dihedral
    vec_dihedral = Vector{Dihedral}(undef, num_dihedrals)
    for cell in 1:num_cells
        for dihedral in 1:num_cell_dihedrals
            id_now = (cell - 1) * num_cell_dihedrals + dihedral
            dihedral_now = dihedral(data_str.dihedral_mode[dihedral])
            dihedral_now.id = id_now
            atom_tilt = (cell - 1) * data_str.num_atoms
            dihedral_now.atom .+= atom_tilt
            vec_dihedral[id_now] = dihedral_now
        end
    end

    # Output
    Data_Unit(data_basic, data_str, data.vec_atom,
            data.vec_bond, vec_angle, vec_dihedral, 0)
end

# genr_improper
"""
    genr_improper(data_cell::Data_Cell, data::Data)

This will be called by `genr(...)` if `improper_mode` is contained in `str`,
genrating information of improper
"""
function genr_improper(data_cell::Data_Cell, data::Data)
    # Reading Input
    data_basic = data.data_basic
    data_str = data.data_str
    num_cells = data_cell.num_cells
    num_cell_impropers = data_str.num_impropers
    num_impropers = num_cells * num_cell_impropers
    num_improper_types = data_str.num_improper_types

    data_basic.num_impropers = num_impropers
    data_basic.num_improper_types = num_improper_types

    # Generate Improper
    vec_improper = Vector{Improper}(undef, num_impropers)
    for cell in 1:num_cells
        for improper in 1:num_cell_impropers
            id_now = (cell - 1) * num_cell_impropers + improper
            improper_now = Improper(data_str.improper_mode[improper])
            improper_now.id = id_now
            atom_tilt = (cell - 1) * data_str.num_atoms
            improper_now.atom .+= atom_tilt
            vec_improper[id_now] = improper_now
        end
    end

    # Output
    Data_Unit(data_basic, data_str, data.vec_atom,
            data.vec_bond, vec_angle, vec_dihedral, vec_improper)
end


text:
using DelimitedFiles
using Test

"""
    get_seat_ids(data)

Get all seat ids. Map to binary strings, parse as Int-base2, add.
"""
function get_seat_ids(data)
    rowbin = Dict('F' => '0', 'B' => '1')
    colbin = Dict('L' => '0', 'R' => '1')
    rows = [parse(Int, join(rowbin[c] for c in row[1:7]), base=2) for row in data]
    cols = [parse(Int, join(colbin[c] for c in row[8:10]), base=2) for row in data]
    return [8*row + col for (row, col) in zip(rows, cols)]
    end

test = readdlm("./data/day05_test.txt")
@test all(get_seat_ids(test) .== [567, 119, 820])

data = readdlm("./data/day05.txt")
seats = get_seat_ids(data)
println("Part 1 max id: $(maximum(seats))")

"""
    get_my-seat(data)

Find missing id for which id-1 and id+1 are present.
Just going to brute force this, there aren't that many.
"""
function get_my_seat(seats)
    for s in 1:maximum(seats)
        if s ∉ seats && s-1 ∈ seats && s+1 ∈ seats
            return s
            end
        end
    end

println("Part 2 my seat: $(get_my_seat(seats))")


text:
using Pkg

Pkg.add(url="https://github.com/snotskie/EpistemicNetworkAnalysis.jl")
Pkg.add("Plots")
Pkg.add("LinearAlgebra")
Pkg.add("DataFrames")
Pkg.add("CSV")
Pkg.add("Statistics")
Pkg.add("HypothesisTests")
Pkg.add("GLM")


text:
# common operations specific to time series data structures: lagging/shifting, differencing, etc.

import Base: diff

function rowdx!(dx::A, x::A, n::Int, r::Int=size(x,1)) where {A<:AbstractArray}
    idx = n > 0 ? (n+1:r) : (1:r+n)
    @inbounds for i=idx
        dx[i,:] .= x[i,:] - x[i-n,:]
    end
    nothing
end

function coldx!(dx::A, x::A, n::Int, c::Int=size(x,2)) where {A<:AbstractArray}
    idx = n > 0 ? (n+1:c) : (1:c+n)
    @inbounds for j=idx
        dx[:,j] .= x[:,j] - x[:,j-n]
    end
    nothing
end

function diffn(x::A, dim::Int=1, n::Int=1) where {A<:AbstractArray}
    @assert dim == 1 || dim == 2 "Argument `dim` must be 1 (rows) or 2 (columns)."
    @assert abs(n) < size(x,dim) "Argument `n` out of bounds."
    if n == 0
        return x
    end
    dx = zeros(eltype(x), size(x))
    if dim == 1
        rowdx!(dx, x, n)
    else
        coldx!(dx, x, n)
    end
    return dx
end

function diff(x::TS{V}, n::Int=1; dim::Int=1, pad::Bool=false, padval::V=zero(eltype(x))) where {V}
    @assert dim == 1 || dim == 2 "Argument dim must be either 1 (rows) or 2 (columns)."
    r = size(x, 1)
    c = size(x, 2)
    dx = diffn(x.values, dim, n)
    if dim == 1
        if pad
            idx = n>0 ? (1:n) : (r+n+1:r)
            dx[idx,:] .= padval
            return TS(dx, x.index, x.fields)
        else
            idx = n > 0 ? (n+1:r) : (1:r+n)
            return TS(dx[idx,:], x.index[idx], x.fields)
        end
    else
        if pad
            idx = n > 0 ? (1:c) : (c+1+1:c)
            dx[:,idx] .= padval
            return TS(dx, x.index, x.fields[idx])
        else
            idx = n > 0 ? (n+1:c) : (1:c+n)
            return TS(dx[:,idx], x.index, x.fields[idx])
        end
    end
end

function lag(x::TS{V}, n::Int=1; pad::Bool=false, padval::V=zero(eltype(x))) where {V}
    @assert abs(n) < size(x,1) "Argument `n` out of bounds."
    if n == 0
        return x
    elseif n > 0
        out = zeros(eltype(x), size(x))
        out[n+1:end,:] .= x.values[1:end-n,:]
    elseif n < 0
        out = zeros(eltype(x), size(x))
        out[1:end+n,:] .= x.values[1-n:end,:]
    end
    r = size(x, 1)
    c = size(x, 2)
    if pad
        idx = n>0 ? (1:n) : (r+n+1:r)
        out[idx,:] .= padval
        return TS(out, x.index, x.fields)
    else
        idx = n > 0 ? (n+1:r) : (1:r+n)
        return TS(out[idx,:], x.index[idx], x.fields)
    end
end

const shift = lag

function pct_change(x::TS{V}, n::Int=1; continuous::Bool=true, pad::Bool=false, padval::V=zero(eltype(x))) where {V}
    if continuous
        return diff(log.(x), n; pad=pad, padval=padval)
    else
        if pad
            return diff(x, n, pad=pad, padval=padval) / x
        else
            return diff(x, n) / x[n+1:end,:]
        end
    end
end


text:
# derived_sequence.jl

"""
A DerivedSequence inherits functionality from an underlying sequence that is
stored. Yet, it is its own type, and hence specific functionality can be added.

Any type that inherits from DerivedSequence has a full interface of a Sequence.
Without any other definitions, the new type is functionally equivalent to the
underlying one.
"""
abstract DerivedSequence{T} <: Sequence{T}
end

"Return the underlying sequence of the derived sequence."
supersequence(s::DerivedSequence) = s.supersequence
# We assume that the underlying sequence is stored in a field called supersequence.

for op in (:iscompact, :nzrange, )
    @eval $op(s::DerivedSequence) = $op(supersequence(s))
end

"An example of a concrete derived sequence."
immutable ConcreteDerivedSequence{T} <: Sequence{T}
    supersequence   ::  Sequence{T}
end


text:
# ----- SampleArray ---------------------

export SampleArray, getduration, toindexdelta

struct SampleArray{T} <: AbstractSampleArray{T}
    # data::AbstractMatrix{T} # MUCH SLOWER! non-interleaving frames x channels
    data::Matrix{T} # non-interleaving frames x channels
    rate::Float64 # in Hz
    names::Vector{Symbol}

    function SampleArray{T}(data::Matrix{T}, rate::Float64, names::Vector{Symbol}) where T
        _check_channel_names(names)
        length(names) == size(data, 2) || throw(DimensionMismatch("the number of names ($(length(names))) does not match the number of channels ($(size(data, 2)))!"))
        n = new(data, rate, copy(names))
    end
end

SampleArray(x::AbstractMatrix{T}, rate::Frequency, names::Union{Nothing,Vector{Symbol}}=nothing) where T =
    SampleArray{T}(x, toHz(rate), isnothing(names) ? _default_channel_names(size(x, 2)) : names)

SampleArray(x::AbstractVector{T}, rate::Frequency, names::Union{Nothing,Vector{Symbol}}=nothing) where T =
    SampleArray{T}(reshape(x, :, 1), toHz(rate), isnothing(names) ? _default_channel_names(size(x, 2)) : names)

domain(x::SampleArray) = range(0, ((nframes(x)-1)/rate(x)); length=nframes(x))

getduration(x::SampleArray) = (nframes(x)/rate(x)) * s

toindex(x::SampleArray{T}, t::Time) where T = round(Int, convert(Float64, tos(t) * rate(x))) + 1

toindexdelta(x::SampleArray, t::Int) = t # delta frames
toindexdelta(x::SampleArray, t::Time) = round(Int, convert(Float64, tos(t) * rate(x))) # delta time

@inline toframeidx(::SampleArray{T}, ::R) where {T, R <: ClosedInterval} = throw(ArgumentError("only Time intervals allowed!"))
@inline toframeidx(x::SampleArray{T}, ti::R) where {T, R <: ClosedInterval{<:Time}} = toindex(x, tos(minimum(ti))):toindex(x, tos(maximum(ti)))

# redefined so channel names are treated
function Base.getindex(x::SampleArray{T}, I::R, J::S) where {T, R <: FrameIndex, S <: ChannelIndex}
    xnew = invoke(Base.getindex, Tuple{AbstractSampleArray{T}, R, S}, x, I, J)
    # println("GETINDEX")
    names!(xnew, names(x)[tochannelidx(x, J)])
    xnew
end

function Base.similar(x::SampleArray, t::Type{T}, dims::Dims) where T
    # tries to copy channel names
    # if there are fever names in the source array use default ones
    # TODO: move to abstracttypes.jl?
    # println("SIMILAR $(dims[2]) $(nchannels(x))")
    ns = dims[2] ≤ nchannels(x) ? x.names[1:dims[2]] : _default_channel_names(dims[2])
    SampleArray{T}(similar(data(x), t, dims), rate(x), ns)
end

function slice(X::SampleArray, t::Time)
    error("missing TEST! is this used anywhere?")
    data(X)[toindex(X, t), :]
end

function Base.show(io::IO, ::MIME"text/plain", x::SampleArray{T}) where T
    dur = nframes(x) / rate(x)
    dom = domain(x)
    step = dom[2] - dom[1]
    clist = join([":$(n)" for n in names(x)], ", ")
    print(io, "SampleArray{$T}: $(nchannels(x)) channels: $(clist)\n $(nframes(x)) frames ≈ $(dur) s, step ≈ $(step) s sampled at $(rate(x)) Hz:\n ")
    print(io, data(x))
end

function DSP.resample(x::SampleArray, to::Frequency)
    r = toHz(to) / rate(x)
    SampleArray(mapslices(c -> resample(c, r), data(x); dims=1), to)
end


text:
using SwitchTimeOpt
using MathOptInterface
const MOI = MathOptInterface
using Ipopt

# Define Solver options
maxiter = 100
maxtime = 1000.0
verbose = 5
tolerance = 1e-06

#Define solver
solver = Ipopt.Optimizer()
MOI.set(solver, MOI.RawOptimizerAttribute("tol"), tolerance)
MOI.set(solver, MOI.RawOptimizerAttribute("print_level"), verbose)
MOI.set(solver, MOI.RawOptimizerAttribute("max_cpu_time"), maxtime)
MOI.set(solver, MOI.RawOptimizerAttribute("max_iter"), maxiter)

### Define system parameters
# Time vector
t0 = 0.0
tf = 12.0

# Integer input
uvec = [repeat([0.0; 1.0], 4, 1); 0.0]
uvec = copy(uvec')

# Number of switching times
N = size(uvec,2) - 1

# Cost funcction matrix
C = [1.0 0.0 -1.0;
     0.0 1.0 -1.0]
Q = C'*C

# Define initial state
x0 = [0.5; 0.7; 1]

### Define system dynamics
function nldyn(x::Array{Float64,1}, u::Array{Float64,1})
  n = length(x)
  f = zeros(n)
  f[1] = x[1] - x[1]*x[2] - 0.4*x[1]*u[1]
  f[2] = -x[2] + x[1]*x[2] - 0.2*x[2]*u[1]
  return f
end

function nldyn_deriv(x::Array{Float64,1}, u::Array{Float64,1})
  df = [1.0-x[2]-0.4*u[1]       -x[1]                   0;
        x[2]                     -1+x[1]-0.2*u[1]       0;
        0                        0                      0]
  return df
end

# Generate and solve problem
ngrid = 300

# Initialize the problem
m = stoproblem(
  x0,                 # Initial state
  nldyn,              # Nonlinear dynamics
  nldyn_deriv,        # Nonlinear dynamics derivative
  uvec,               # Vector of integer inputs
  ngrid = ngrid,      # Number of grid points
  t0 = t0,            # InitialtTime
  tf = tf,            # Final time
  Q = Q,              # Cost matrix
  solver = solver)

solve!(m)
objlin = getobjval(m)
tauopt = gettau(m)
status = getstat(m)
xsim, xpts, objode45, t = simulate(m, tauopt)

println("Status: ", status)
println("Objective value: ", objlin)

showPlots = true
savePlots = false
plotSolution(m, "Lotka-Volterra", "Lotka-Volterra", showPlots, savePlots)


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

JLLWrappers.@generate_main_file_header("ONNXRuntime")
JLLWrappers.@generate_main_file("ONNXRuntime", UUID("09e6dd1b-8208-5c7e-a336-6e9061773d0b"))
end  # module ONNXRuntime_jll


text:
module Plotting

export plotLogReg,plotArray,plotDoubleArray,plotTripleArray,plotQuadArray

using Gadfly
using DataFrames
using RDatasets

iris = dataset("datasets", "iris")

function plotDoubleArray(array1,array2,filename)
    df = DataFrame(iter = [identity(x) for x in 1:size(array1,1)], obj = array1,obj2 = array2)
    p = plot(layer(df,x=:iter,y=:obj,Geom.line,Theme(default_color=colorant"red")),
            layer(df,x=:iter,y=:obj2,Geom.line,Theme(default_color=colorant"blue")))
    img = SVG(filename, 10inch, 10inch)
    draw(img, p)
end

function plotTripleArray(array1,array2,array3,filename)
    df = DataFrame(iter = [identity(x) for x in 1:size(array1,1)], obj = array1,obj2 = array2,obj3 = array3)
    p = plot(layer(df,x=:iter,y=:obj,Geom.line,Theme(default_color=colorant"red")),
            layer(df,x=:iter,y=:obj2,Geom.line,Theme(default_color=colorant"blue")),
            layer(df,x=:iter,y=:obj3,Geom.line,Theme(default_color=colorant"orange")))
    img = SVG(filename, 10inch, 10inch)
    draw(img, p)
end

function plotQuadArray(array1,array2,array3,array4,filename)
    df = DataFrame(iter = [identity(x) for x in 1:size(array1,1)], obj = array1,obj2 = array2,obj3 = array3,obj4 = array4)
    p = plot(layer(df,x=:iter,y=:obj,Geom.line,Theme(default_color=colorant"red")),
            layer(df,x=:iter,y=:obj2,Geom.line,Theme(default_color=colorant"blue")),
            layer(df,x=:iter,y=:obj3,Geom.line,Theme(default_color=colorant"orange")),
            layer(df,x=:iter,y=:obj4,Geom.line,Theme(default_color=colorant"green"))
            )
    img = SVG(filename, 6inch, 4inch)
    draw(img, p)
end

function plotArray(array,filename)
    df = DataFrame(iter = [identity(x) for x in 1:size(array,1)], obj = array)
    p = plot(layer(df,x=:iter,y=:obj,Geom.line,Theme(default_color=colorant"red")))
    img = SVG(filename, 6inch, 4inch)
    draw(img, p)
end

function plotLogReg(X)
    p = plot(iris, x=:X, y=:y, Geom.line)
    img = SVG("iris_plot.svg", 6inch, 4inch)
    draw(img, p)
end

end

max_stars_repo_name: zagazao/num-opt

max_issues_repo_name: zagazao/num-opt

max_forks_repo_name: zagazao/num-opt


text:
using Test
import Pkg

@testset "Chemostat.jl" begin

    for pkg in ["MetNets", "MetLP", "MetEP"]
        println("-"^60)
        println("Tesing: ", pkg)
        Pkg.test(pkg)
        println()
        @test true
    end
end

max_stars_repo_name: josePereiro/Chemostat

max_issues_repo_name: josePereiro/Chemostat

max_forks_repo_name: josePereiro/Chemostat


text:
module PhotonTest

using Photon

function test()

  model = Sequential(
    Conv2D(16, (3,3), relu),
    Conv2D(16, (3,3), relu),
    BatchNorm(),
    MaxPool2D(4),
    Conv2D(64, (3,3), relu),
    Conv2D(64, (3,3), relu),
    BatchNorm(),
    MaxPool2D(4),
    Conv2D(256, (3,3), relu),
    Conv2D(256, (3,3), relu),
    BatchNorm(),
    MaxPool2D(4),
    Dense(64, relu),
    Dense(10, relu)
  )

  data = randn(Float32,224,224,3,16)

  model(KorA(data))

  @time for i in 1:1000
    X = KorA(data)
    model(X)
  end
end

test()

end


text:
"""
	getmodel(bbin::String)
	getmodel(bbin::Array{Float64,1})

Return LAR model (V,EV,FV) of a box, aligned or not to axes.
"""
function Common.getmodel(bbin::String)
	# path of voume json
	return json2LARvolume(bbin)
end

function Common.getmodel(bbin::Array{Float64,1})
	# input format -> bbin = [x_min y_min z_min x_max y_max z_max ]
	bb = AABB(bbin[4], bbin[1], bbin[5], bbin[2], bbin[6], bbin[3])
	return Common.getmodel(bb)
end

"""
	searchfile(path::String,key::String)

In recursive mode, search all files in `path` with `key` in filename.
"""
function searchfile(path::String,key::String)
	files = String[]
	for (root, _, _) in walkdir(path)
		thisfiles = filter(x->occursin(key,x), readdir(root,join=true))
		union!(files,thisfiles)
	end
	return files
end

"""
	mkdir_project(folder::String, project_name::String)

Create a new folder `project_name` in `folder`.
"""
function mkdir_project(folder::String, project_name::String)
	@assert isdir(folder) "$folder not an existing folder"
	proj_folder = joinpath(folder,project_name)

	mkdir_if(proj_folder)

	return proj_folder
end

"""
	mkdir_if(folder::String)

Create `folder` if not exist.
"""
function mkdir_if(folder::String)
	if !isdir(folder)
		mkdir(folder)
	end
end

"""
	clearfolder(folder::String)

Clear the given `folder`.
"""
function clearfolder(folder::String)
	root, dirs, files = first(walkdir(folder))
	for dir in dirs
		rm(joinpath(root,dir),recursive=true)
	end
	for file in files
		rm(joinpath(root,file))
	end
	return 1
end

"""
	source2pc(source::String, lod::Union{Nothing,Int64})

Get point cloud from source.
"""
function source2pc(source::String, lod=0::Int64)

	if isdir(source) # se source è un potree
		println("Read data from: Potree struct")
		cloud_metadata = CloudMetadata(source)

		if lod == -1
			trie = potree2trie(source)
			max_level = FileManager.max_depth(trie)
			all_files = FileManager.get_all_values(trie)
			PC = FileManager.las2pointcloud(all_files...)
			return PC
		else
			all_files = FileManager.get_files_in_potree_folder(source,lod)
			PC = FileManager.las2pointcloud(all_files...)
			return PC
		end

	elseif isfile(source) # se source è un file
		println("Read data from: Single file")
		PC = FileManager.las2pointcloud(source)
		return PC
	end

end

"""
	successful(test::Bool,folder::String; message=""::String)
"""
function successful(test::Bool,folder::String; message=""::String, filename = "execution.probe")
	if test
		io = open(joinpath(folder,filename),"w")
		write(io, message)
		close(io)
	end
end

# """
# Read file line by line.
# """
# function get_directories(filepath::String)
# 	return readlines(filepath)
# end

# """
# 	boxmodelfromjson(volume::String)
#
# Read volume model from a file JSON.
# """
# function boxmodel_from_json(volume::String)
# 	V,EV,FV = json2LARvolume(volume)
# 	return V,EV,FV
# end


text:
using Test

import ReinforcementLearning as RL
import ReinforcementLearning.FiniteMarkovDecisionProcesses as MDP

include("FiniteMarkovDecisionProcesses/FiniteMarkovDecisionProcesses.jl")


text:
using Documenter, PlotDocs, Plots

# Set matplotlib gui backend
ENV["MPLBACKEND"]="agg"

# Initialize backends
pyplot()
pgfplots()
pgfplotsx()
plotly()
gr()

const PAGES = Any[
    "Home" => "index.md",
    "Getting Started" => [
        "Installation" => "install.md",
        "Basics" => "basics.md",
        "Tutorial" => "tutorial.md",
    ],
    "Manual" => [
        "Input Data" => "input_data.md",
        "Output" => "output.md",
        "Attributes" => [
            "Overview" => "attributes.md",
            "Series Attributes" => "generated/attributes_series.md",
            "Plot Attributes" => "generated/attributes_plot.md",
            "Subplot Attributes" => "generated/attributes_subplot.md",
            "Axis Attributes" => "generated/attributes_axis.md",
        ],
        "Layouts" => "layouts.md",
        "Recipes" => "recipes.md",
        "Colors" => "colors.md",
        "Animations" => "animations.md",
        "Backends" => "backends.md",
        "Supported Attributes" => "generated/supported.md",
    ],
    "Learning" => "learning.md",
    "Contributing" => "contributing.md", # TODO: testing
    "Ecosystem" => [
        "Overview" => "ecosystem.md",
        "GraphRecipes" => [
            "Introduction" => "graphrecipes/introduction.md",
            "Examples" => "graphrecipes/examples.md",
            "Attributes" => "generated/graph_attributes.md",
        ],
    ],
    "Advanced Topics" => [
        "Internals" => "pipeline.md",
    ],
]

const EXAMPLES = "Examples" => [
    "GR" => "generated/gr.md",
    "Plotly" => "generated/plotly.md",
    "PyPlot" => "generated/pyplot.md",
    "PGFPlots" => "generated/pgfplots.md",
    "PGFPlotsX" => "generated/pgfplotsx.md",
    "UnicodePlots" => "examples/unicodeplots.md",
    "InspectDR" => "examples/inspectdr.md",
]

function builddocs(; examples=true)
    generate_attr_markdown()
    generate_supported_markdown()
    generate_graph_attr_markdown()
    # save_html(PlotDocs.GRAPH_ATTR_DF, "graph_attr.html")
    for be in (:gr, :plotly, :pyplot, :pgfplots, :pgfplotsx)
        if examples
            generate_markdown(be)
        else
            rm(joinpath(GENDIR, "$be.md"), force=true)
        end
    end
    pages = examples ? [PAGES; EXAMPLES] : PAGES
    makedocs(
        format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
        sitename = "Plots",
        authors = "Thomas Breloff",
        pages = pages,
    )
end

builddocs()

deploydocs(
    repo = "github.com/JuliaPlots/PlotDocs.jl.git",
    push_preview = true,
)


text:
# Author: Lukasz Bratos
include("methods.jl")

# Różnica naszych funkcji
f(x) = exp(x) - 3 * x
g(x) = 3 * x - exp(x)
delta = 10^-4
epsilon = 10^-4

println("Przedział [0.0, 1.0: ", Methods.mbisekcji(f, 0.0, 1.0, delta, epsilon))
println("Przedział [1.0, 2.0: ", Methods.mbisekcji(f, 1.0, 2.0, delta, epsilon))

println("Przedział [0.0, 1.0: ", Methods.mbisekcji(g, 0.0, 1.0, delta, epsilon))
println("Przedział [1.0, 2.0: ", Methods.mbisekcji(g, 1.0, 2.0, delta, epsilon))


text:
using DataFrames

function lagrangianpolynomial(np::Int, x::Vector{Float64}, y::Vector{Float64}, xi::Vector{Float64})
  local yi = Vector{Float64}(undef, length(xi))
  m = sort(hcat(x, y), dims=1)
  for l in 1:length(xi)
    for i in 1:np
      term = 1.0
      for j in 1:np
        j !== i && (term *= (xi[l] - m[j, 1])/(m[i, 1] - m[j, 1]))
      end
      yi[l] += term * m[i, 2]
    end
  end
  columns=[]
  push!(columns, m[:,1])
  push!(columns, m[:,2])
  df = DataFrame(columns, [:x, :y])

  columns=[]
  push!(columns, xi)
  push!(columns, yi)
  dfxi = DataFrame(columns, [:xi, :yi])

  (df, dfxi)
end


text:
# Return log density for population of states
function eval_population( theta::Array{Float64,2}, calc_log_pdf::Function )
   popsize = size(theta,2)
   logpdf = Array{Float64}( popsize )
   for i in 1:popsize
      logpdf[i] = calc_log_pdf(theta[:,i])
   end
   return logpdf
end

## Differential Evolution Sampler
#    based on ter Braak (2006) and Nelson et al. (2014)

# Generate ensemble of trial states according to DEMCMC
function generate_trial_states_demcmc( theta::Array{Float64,2}; gamma_o::Float64 = 2.38/sqrt(2*size(theta,1)), epsilon = 0.01  )
   num_param = size(theta,1)
   num_pop = size(theta,2)
   theta_trial = similar(theta)

   for i in 1:num_pop
      # Choose head (k) and tail (j) for perturbation vector
  j = rand(1:(num_pop-1))
  if( j >= i ) j = j+1 end
  k = rand(1:(num_pop-2))
  if( k >= i ) k = k+1 end
  if( k >= j ) k = k+1 end

  # Choose scale factor
  scale = (1.0 + epsilon*randn()) * gamma_o

  theta_trial[:,i] = theta[:,i] + scale * (theta[:,k] - theta[:,j])
   end
   return theta_trial
end

# Evolve population of states (theta_init) with target density
function run_demcmc( theta_init::Array{Float64,2}, calc_log_pdf::Function; num_gen::Integer = 100, epsilon::Float64 = 0.01 )
   @assert(num_gen >=1 )
   num_param = size(theta_init,1)
   num_pop = size(theta_init,2)
   @assert(num_pop > num_param)

   # Allocate arrays before starting loop
   pop = Array{Float64}(num_param, num_pop, num_gen)
   poplogpdf = Array{Float64}(num_pop, num_gen )
   accepts_chain = zeros(Int64,num_pop )
   rejects_chain = zeros(Int64,num_pop )
   accepts_generation = zeros(Int64,num_gen )
   rejects_generation = zeros(Int64,num_gen )

   # Make a working copy of the current population of parameter values
   theta = copy(theta_init)
   logpdf = eval_population( theta, calc_log_pdf )

   for g in 1:num_gen # first(pop.range):last(pop.range)  # loop over generations
      # Generate population of trial sets of model parameters
      gamma_o = (mod(g,10)!=0) ? 2.38/sqrt(2*num_param) : 1.0 # every 10th generation try full-size steps
      theta_trial = generate_trial_states_demcmc( theta, gamma_o=gamma_o, epsilon=epsilon )

      # Evaluate model for each set of trial parameters
      logpdf_trial = eval_population( theta_trial, calc_log_pdf )

      # For each member of population
      for i in 1:num_pop
        log_pdf_ratio = logpdf_trial[i] - logpdf[i]
        if( (log_pdf_ratio>0) || (log_pdf_ratio>log(rand())) )     # Decide whether to Accept
            theta[:,i] = theta_trial[:,i]
            logpdf[i] = logpdf_trial[i]
            accepts_chain[i] = accepts_chain[i]+1
            accepts_generation[g] = accepts_generation[g]+1
        else
            rejects_chain[i] = rejects_chain[i]+1
            rejects_generation[g] = rejects_generation[g]+1
        end
      end

      # Log results
      if(true)  # (in(g,pop.range))
        pop[:,:,g] = copy(theta)
        poplogpdf[:,g] = copy(logpdf)
      end
   end
   return Dict("theta_last"=>theta, "logpdf_last"=>logpdf, "theta_all"=>pop, "logpdf_all"=>poplogpdf,
           "accepts_chain"=>accepts_chain, "rejects_chain"=>rejects_chain, "accepts_generation"=>accepts_generation, "rejects_generation"=>rejects_generation )
end