text:
export export_SDP

function export_SDP(sdp::SDPPrimal, path; indentedprint=true, renamemoments=true)

    # Build moment shortname dict if required
    momentdict = build_momentdict(sdp, renamemoments)

    # Collect all constraints keys
    ctr_keys = build_ctrkeysset(sdp)

    # Export blocks of constraints
    blocks_file = joinpath(path, "matrix.sdp")
    !isfile(blocks_file) || rm(blocks_file)

    fblocks = open(blocks_file, "a")
    print_blocksfile(fblocks, sdp.blocks, momentdict=momentdict, indentedprint=indentedprint)
    close(fblocks)

    # Export linear part of constraints
    lin_file = joinpath(path, "lin.sdp")
    !isfile(lin_file) || rm(lin_file)

    flin = open(lin_file, "a")
    print_linfile(flin, sdp.lin, sdp.linsym, momentdict=momentdict, indentedprint=indentedprint)
    close(flin)

    # Export constants of constraints
    cst_file = joinpath(path, "const.sdp")
    !isfile(cst_file) || rm(cst_file)

    fcst = open(cst_file, "a")
    print_cstfile(fcst, sdp.cst, momentdict=momentdict, ctr_keys=ctr_keys, indentedprint=indentedprint)
    close(fcst)

    # Export bloc types
    types_file = joinpath(path, "types.sdp")
    !isfile(types_file) || rm(types_file)

    ftypes = open(types_file, "a")
    print_typesfile(ftypes, sdp.block_to_vartype)
    close(ftypes)

    # Export moment dictionary name if required
    momname_file = joinpath(path, "names.sdp")
    !isfile(momname_file) || rm(momname_file)

    if renamemoments
        fname = open(momname_file, "a")
        print_namesfile(fname, momentdict)
        close(fname)
    end

    return nothing
end

"""
    momentdict = build_momentdict(sdp, renamemoments)

    Build a dict `momentdict::Dict{Exponent, String}` that stores short names for moments (Exponent) of `sdp`.
    NOTE: should be extended to Variable and (String, Variable) later...
"""
function build_momentdict(sdp, renamemoments::Bool)
    momentdict = Dict{Exponent, String}()

    n_moment=0
    n_matvar=0
    n_scalvar=0

    momentdict[Exponent()] = "1"
    ## Treating blocks
    for (moment, blockname, γ, δ) in keys(sdp.blocks)
        α, β = moment.conj_part, moment.expl_part
        if renamemoments
            haskey(momentdict, α) || (n_moment+=1; momentdict[α] = shortname_moment(n_moment))
            haskey(momentdict, β) || (n_moment+=1; momentdict[β] = shortname_moment(n_moment))
            haskey(momentdict, γ) || (n_moment+=1; momentdict[γ] = shortname_moment(n_moment))
            haskey(momentdict, δ) || (n_moment+=1; momentdict[δ] = shortname_moment(n_moment))
        else
            haskey(momentdict, α) || (momentdict[α] = format_string(α))
            haskey(momentdict, β) || (momentdict[β] = format_string(β))
            haskey(momentdict, γ) || (momentdict[γ] = format_string(γ))
            haskey(momentdict, δ) || (momentdict[δ] = format_string(δ))
        end
    end

    for (moment, var) in keys(sdp.lin)
        α, β = moment.conj_part, moment.expl_part
        if renamemoments
            haskey(momentdict, α) || (n_moment+=1; momentdict[α] = shortname_moment(n_moment))
            haskey(momentdict, β) || (n_moment+=1; momentdict[β] = shortname_moment(n_moment))
        else
            momentdict[α] = format_string(α)
            momentdict[β] = format_string(β)
        end
    end

    for (moment, blockname, var) in keys(sdp.linsym)
        α, β = moment.conj_part, moment.expl_part
        if renamemoments
            haskey(momentdict, α) || (n_moment+=1; momentdict[α] = shortname_moment(n_moment))
            haskey(momentdict, β) || (n_moment+=1; momentdict[β] = shortname_moment(n_moment))
        else
            momentdict[α] = format_string(α)
            momentdict[β] = format_string(β)
        end
    end

    for moment in keys(sdp.cst)
        α, β = moment.conj_part, moment.expl_part
        if renamemoments
            haskey(momentdict, α) || (n_moment+=1; momentdict[α] = shortname_moment(n_moment))
            haskey(momentdict, β) || (n_moment+=1; momentdict[β] = shortname_moment(n_moment))
        else
            momentdict[α] = format_string(α)
            momentdict[β] = format_string(β)
        end
    end

    return momentdict
end

function build_ctrkeysset(sdp::SDPPrimal{T}) where T
    ctr_keys = Set{Moment}()

    for ((moment, blockname, γ, δ), λ) in sdp.blocks
        push!(ctr_keys, moment)
    end
    for ((moment, var), λ) in sdp.lin
        push!(ctr_keys, moment)
    end
    for ((moment, blockname, var), λ) in sdp.linsym
        push!(ctr_keys, moment)
    end
    for (moment, λ) in sdp.cst
        push!(ctr_keys, moment)
    end
    return ctr_keys
end

function print_blocksfile(io::IO, sdpblocks::Dict{Tuple{Moment, String, Exponent, Exponent}, T};
                                                        momentdict::Dict{Exponent, String}=Dict{Exponent, String}(),
                                                        indentedprint=false,
                                                        print_header=true) where T
    if print_header
        println(io, "## Description of the matrices A_ji for the problem:")
        println(io, "##         max     ∑ A_0i[k,l] × Zi[k,l] + ∑ b_0[k] × x[k] + c_0")
        println(io, "##         s.t.    ∑ A_ji[k,l] × Zi[k,l] + ∑ b_j[k] × x[k] + c_j  ==  0")
        println(io, "## Constraints keys are j → (j_conj, j_expl, clique).")
        println(io, "## Objective keys are 0 → (1,1, *) for any *.")
        println(io, "#")
    end

    cstrlenα = maximum(x->length(haskey(momentdict, x[1].conj_part)?momentdict[x[1].conj_part]: string(x[1].conj_part)), keys(sdpblocks))
    cstrlenα= max(cstrlenα, length("#j_conj"))
    cstrlenβ = maximum(x->length(haskey(momentdict, x[1].expl_part)?momentdict[x[1].expl_part]: string(x[1].expl_part)), keys(sdpblocks))
    cstrlenβ= max(cstrlenβ, length("j_expl"))
    cliquelen = maximum(x->length(x[1].clique), keys(sdpblocks))
    cliquelen= max(cliquelen, length("clique"))
    blocklen = maximum(x->length(x[2]), keys(sdpblocks))
    blocklen= max(blocklen, length("Zi"))
    rowlen = maximum(x->length(haskey(momentdict, x[3])?momentdict[x[3]]:string(x[3])), keys(sdpblocks))
    rowlen = max(rowlen, length("k"))
    collen = maximum(x->length(haskey(momentdict, x[4])?momentdict[x[4]]:string(x[4])), keys(sdpblocks))
    collen = max(collen, length("l"))

    print_string(io, "#j_conj", cstrlenα, indentedprint=indentedprint)
    print_string(io, "j_expl", cstrlenβ, indentedprint=indentedprint)
    print_string(io, "clique", cliquelen, indentedprint=indentedprint)
    print_string(io, "Zi", blocklen, indentedprint=indentedprint)
    print_string(io, "k", rowlen, indentedprint=indentedprint)
    print_string(io, "l", collen, indentedprint=indentedprint)
    print_string(io, "Real(A_ij[k, l])", 23, indentedprint=indentedprint)
    print_string(io, "Imag(A_ij[k, l])", 23, indentedprint=indentedprint)
    println(io)

    for (moment, blockname, γ, δ) in sort(collect(keys(sdpblocks)))
        λ = sdpblocks[(moment, blockname, γ, δ)]

        α, β = moment.conj_part, moment.expl_part
        print_string(io, haskey(momentdict, α)?momentdict[α]: string(α), cstrlenα, indentedprint=indentedprint)
        print_string(io, haskey(momentdict, β)?momentdict[β]: string(β), cstrlenβ, indentedprint=indentedprint)
        print_string(io, moment.clique, cliquelen, indentedprint=indentedprint)
        print_string(io, blockname, blocklen, indentedprint=indentedprint)
        print_string(io, haskey(momentdict, γ)?momentdict[γ]: string(γ), rowlen, indentedprint=indentedprint)
        print_string(io, haskey(momentdict, δ)?momentdict[δ]: string(δ), collen, indentedprint=indentedprint)
        @printf(io, "% .16e % .16e\n", real(λ), imag(λ))
    end
end

function print_linfile(io::IO, sdplin::Dict{Tuple{Moment, Exponent}, T}, sdplinsym::Dict{Tuple{Moment, String, Exponent}, T};
                                                                     momentdict::Dict{Exponent, String}=Dict{Exponent, String}(),
                                                                     indentedprint=false,
                                                                     print_header=true) where T
    if print_header
        println(io, "## Description of the vectors b_j for the problem:")
        println(io, "##         max     ∑ A_0i[k,l] × Zi[k,l] + ∑ b_0[k] × x[k] + c_0")
        println(io, "##         s.t.    ∑ A_ji[k,l] × Zi[k,l] + ∑ b_j[k] × x[k] + c_j  ==  0")
        println(io, "## Constraints keys are j → (j_conj, j_expl, clique).")
        println(io, "## Objective keys are 0 → (1,1, *) for any *.")
        println(io, "#")
    end

    cstrlenα = 0
    cstrlenβ = 0
    cliquelen = 0
    varlen = varlensym = 0

    if length(union(keys(sdplin), keys(sdplinsym))) > 0
        cstrlenα = maximum(x->length(haskey(momentdict, x[1].conj_part)?momentdict[x[1].conj_part]: string(x[1].conj_part)), union(keys(sdplin), keys(sdplinsym)))
        cstrlenβ = maximum(x->length(haskey(momentdict, x[1].expl_part)?momentdict[x[1].expl_part]: string(x[1].expl_part)), union(keys(sdplin), keys(sdplinsym)))
        cliquelen = maximum(x->length(x[1].clique), union(keys(sdplin), keys(sdplinsym)))
        varlen = length(sdplin)!=0 ? maximum(x->length(format_string(x[2])), keys(sdplin)) : 0
        varlensym = length(sdplinsym)!=0 ? maximum(x->length(format_string(x[3], x[2])), keys(sdplinsym)) : 0
    end

    cstrlenα= max(cstrlenα, length("#j_conj"))
    cstrlenβ= max(cstrlenβ, length("j_expl"))
    cliquelen= max(cliquelen, length("clique"))
    varlen = max(varlen, varlensym, length("x[k]"))

    print_string(io, "#j_conj", cstrlenα, indentedprint=indentedprint)
    print_string(io, "j_expl", cstrlenβ, indentedprint=indentedprint)
    print_string(io, "clique", cliquelen, indentedprint=indentedprint)
    print_string(io, "x[k]", varlen, indentedprint=indentedprint)
    print_string(io, "Real(b_j[k])", 23, indentedprint=indentedprint)
    print_string(io, "Imag(b_j[k])", 23, indentedprint=indentedprint)
    println(io)

    if length(sdplin)!=0
        for (moment, var) in sort(collect(keys(sdplin)))
            λ = sdplin[(moment, var)]

            α, β = moment.conj_part, moment.expl_part
            print_string(io, haskey(momentdict, α)?momentdict[α]: string(α), cstrlenα, indentedprint=indentedprint)
            print_string(io, haskey(momentdict, β)?momentdict[β]: string(β), cstrlenβ, indentedprint=indentedprint)
            print_string(io, moment.clique, cliquelen, indentedprint=indentedprint)
            print_string(io, format_string(var), varlen, indentedprint=indentedprint)
            @printf(io, "% .16e % .16e\n", real(λ), imag(λ))
        end
    end
    if length(sdplinsym) != 0
        for (moment, blockname, var) in sort(collect(keys(sdplinsym)))
            λ = sdplinsym[(moment, blockname, var)]

            α, β = moment.conj_part, moment.expl_part
            print_string(io, haskey(momentdict, α)?momentdict[α]: string(α), cstrlenα, indentedprint=indentedprint)
            print_string(io, haskey(momentdict, β)?momentdict[β]: string(β), cstrlenβ, indentedprint=indentedprint)
            print_string(io, moment.clique, cliquelen, indentedprint=indentedprint)
            print_string(io, format_string(var, blockname), varlen, indentedprint=indentedprint)
            @printf(io, "% .16e % .16e\n", real(λ), imag(λ))
        end
    end
end

function print_cstfile(io::IO, sdpcst::Dict{Moment, T};
                                momentdict::Dict{Exponent, String}=Dict{Exponent, String}(),
                                ctr_keys::Set{Moment}=Set{Moment}(),
                                indentedprint=false,
                                print_header=true) where T
    if print_header
        println(io, "## Description of the scalars c_j for the problem:")
        println(io, "##         max     ∑ A_0i[k,l] × Zi[k,l] + ∑ b_0[k] × x[k] + c_0")
        println(io, "##         s.t.    ∑ A_ji[k,l] × Zi[k,l] + ∑ b_j[k] × x[k] + c_j  ==  0")
        println(io, "## Constraints keys are j → (j_conj, j_expl, clique).")
        println(io, "## Objective keys are 0 → (1,1, *) for any *.")
        println(io, "#")
    end

    cstrlenα = maximum(x->length(haskey(momentdict, x.conj_part)?momentdict[x.conj_part]: string(x.conj_part)), union(ctr_keys, keys(sdpcst)))
    cstrlenα= max(cstrlenα, length("#j_conj"))
    cstrlenβ = maximum(x->length(haskey(momentdict, x.expl_part)?momentdict[x.expl_part]: string(x.expl_part)), union(ctr_keys, keys(sdpcst)))
    cstrlenβ= max(cstrlenβ, length("j_expl"))
    cliquelen = maximum(x->length(x.clique), union(ctr_keys, keys(sdpcst)))
    cliquelen= max(cliquelen, length("clique"))

    print_string(io, "#j_conj", cstrlenα, indentedprint=indentedprint)
    print_string(io, "j_expl", cstrlenβ, indentedprint=indentedprint)
    print_string(io, "clique", cliquelen, indentedprint=indentedprint)
    print_string(io, "Real(c_j)", 23, indentedprint=indentedprint)
    print_string(io, "Imag(c_j)", 23, indentedprint=indentedprint)
    println(io)

    for moment in sort(union(ctr_keys, keys(sdpcst)))
        α, β = moment.conj_part, moment.expl_part
        λ = haskey(sdpcst, moment)?sdpcst[moment]:0
        print_string(io, haskey(momentdict, α)?momentdict[α]: string(α), cstrlenα, indentedprint=indentedprint)
        print_string(io, haskey(momentdict, β)?momentdict[β]: string(β), cstrlenβ, indentedprint=indentedprint)
        print_string(io, moment.clique, cliquelen, indentedprint=indentedprint)
        @printf(io, "% .16e % .16e\n", real(λ), imag(λ))
    end
end

function print_typesfile(io::IO, block_to_vartype::Dict{String, Symbol})
    println(io, "## Description of the matrix and scalar variables Zi and x[k] for the problem:")
    println(io, "##         max     ∑ A_0i[k,l] × Zi[k,l] + ∑ b_0[k] × x[k] + c_0")
    println(io, "##         s.t.    ∑ A_ji[k,l] × Zi[k,l] + ∑ b_j[k] × x[k] + c_j  ==  0")
    println(io, "## Matrix variable types Zi are \"SDPC\" or \"SDP\".")
    println(io, "## By default, scalar variables x[k] are assumed to be real, free.")
    println(io, "#")

    cstrlen = maximum(x->length(x), keys(block_to_vartype))
    cstrlen = max(cstrlen, length("#Zi"))

    print_string(io, "#Zi", cstrlen, alignright=false); println(io, " type")

    for blockname in sort(collect(keys(block_to_vartype)))
        vartype = block_to_vartype[blockname]
        if vartype in Set([:SDP, :SDPC])
            print_string(io, blockname, cstrlen)
            println(io, " $(string(vartype))")
        end
        # warn(LOGGER, "Ignoring variable $blockname of type $vartype") ## NOTE: better logging system.
    end
end

function print_namesfile(io::IO, momentdict::Dict{Exponent, String})
    println(io, "## Description of the scalars c_j for the problem:")
    println(io, "##         max     ∑ A_0i[k,l] × Zi[k,l] + ∑ b_0[k] × x[k] + c_0")
    println(io, "##         s.t.    ∑ A_ji[k,l] × Zi[k,l] + ∑ b_j[k] × x[k] + c_j  ==  0")
    println(io, "## Constraints keys are j → (j_conj, j_expl, clique).")
    println(io, "## Objective keys are 0 → (1,1, *) for any *.")
    println(io, "#")
    println(io, "#shortname  Explicit_name")

    for α in sort(collect(keys(momentdict)))
        shortname = momentdict[α]
        println(io, "$shortname $(format_string(α))")
    end
end


text:
using ExpFamily
using Base.Test

x = randn(10)

ssNP  = suffstats(GaussianNatParam, x)
ssMP  = suffstats(GaussianMeanParam, x)
ssdNP = suffstats(DiagGaussianNatParam, x)
ssdMP = suffstats(DiagGaussianMeanParam, x)

ssNP2   = 2*ssNP
ssdNP2  = ssdNP*2
ssNP2b  = ssNP + ssNP
ssdNP2b = ssdNP + ssdNP

@test isapprox(ssNP[1],   x)
@test isapprox(ssMP[2],   (x*x')/2)
@test isapprox(ssMP[1],   x)
@test isapprox(ssMP[2],   (x*x')/2)
@test isapprox(ssdNP[1],  x)
@test isapprox(ssdNP[2],  (x.^2)/2)
@test isapprox(ssdMP[1],  x)
@test isapprox(ssdMP[2],  (x.^2)/2)
@test isapprox(ssNP2[1],  2x)
@test isapprox(ssNP2[2],  (x*x'))
@test isapprox(ssdNP2[1], 2x)
@test isapprox(ssdNP2[2], x.^2)
@test isapprox(ssNP2[1],  ssNP2b[1])
@test isapprox(ssNP2[2],  ssNP2b[2])
@test isapprox(ssdNP2[1], ssdNP2b[1])
@test isapprox(ssdNP2[2], ssdNP2b[2])

@test isapprox( vec(ssNP),              vcat(x,vec(x*x'/2)))
@test isapprox( norm(ssNP),             norm(vcat(x,vec(x*x'/2))))
@test isapprox( vec(2*ssNP-ssNP),       vec(ssNP))
@test isapprox( vec((2*ssNP)/2),        vec(ssNP))
@test isapprox( vec((5*ssdNP-ssdNP)/4), vec(ssdNP))

ext: jl


text:
# AUTO GENERATED FILE - DO NOT EDIT

export dashbio_igv

"""
    dashbio_igv(;kwargs...)

An Igv component.
The Igv component is an interactive genome visualization component
developed by the Integrative Genomics Viewer (IGV) team. It uses an
example integration of igv.js and React (https://www.npmjs.com/package/igv).
Keyword arguments:
- `id` (String; optional): The ID of this component, used to identify dash components in callbacks.
The ID needs to be unique across all of the components in an app.
- `className` (String; optional): className of the component div.
- `genome` (String; optional): String identifier defining genome (e.g. "hg19"). See https://github.com/igvteam/igv.js/wiki/Reference-Genome
    for details and list of supported identifiers. Note: One (but only one) of
    either genome or reference properties must be set. If both are set,
    the genome property will be ignored.
- `locus` (String; optional): Initial genomic location(s). Either a string or an array of strings.
    If an array a viewport is created for each location.
- `minimumBases` (Real; optional): Minimum window size in base pairs when zooming in
- `reference` (Dict; optional): Object defining reference genome. see https://github.com/igvteam/igv.js/wiki/Reference-Genome
    Note: One (but only one) of either genome or reference properties must be set. If both are set,
    the genome property will be ignored.
- `style` (Dict; optional): Generic style overrides on the plot div
- `tracks` (Array; optional): Array of configuration objects defining tracks initially displayed when app launches.
    see https://github.com/igvteam/igv.js/wiki/Tracks-2.0
"""
function dashbio_igv(; kwargs...)
        available_props = Symbol[:id, :className, :genome, :locus, :minimumBases, :reference, :style, :tracks]
        wild_props = Symbol[]
        return Component("dashbio_igv", "Igv", "dash_bio", available_props, wild_props; kwargs...)
end


text:
@testset "Knapsack" begin

    @testset "Knapsack: ZeroOnePack!" begin
        dp = zeros(Int, 30)
        @test zero_one_pack!(20, [1,3,11], [2,5,30], dp) == 37
        dp = zeros(Int, 30)
        @test zero_one_pack!(10, [1,3,11], [20,5,80], dp) == 25
    end

    @testset "Knapsack: CompletePack!" begin
        dp = zeros(Int, 30)
        @test complete_pack!(20, [1,2,9], [1,3,20], dp) == 43
        dp = zeros(Int, 30)
        @test complete_pack!(10, [1,3,11], [20,5,80], dp) == 200
    end

end

max_stars_repo_name: ashwani-rathee/Julia

max_issues_repo_name: ashwani-rathee/Julia

max_forks_repo_name: ashwani-rathee/Julia


text:
module FeatureTransforms

using Dates: TimeType, Period, Day, hour
using NamedDims: dim
using Statistics: mean, std
using Tables

export Transform, transform, transform!
export HoD, LinearCombination, OneHotEncoding, Periodic, Power
export AbstractScaling, IdentityScaling, MeanStdScaling

include("utils.jl")
include("traits.jl")
include("transform.jl")
include("apply.jl")

# Transform implementations
include("linear_combination.jl")
include("one_hot_encoding.jl")
include("periodic.jl")
include("power.jl")
include("scaling.jl")
include("temporal.jl")

include("test_utils.jl")

# TODO: remove in v0.4 https://github.com/invenia/FeatureTransforms.jl/issues/82
Base.@deprecate_binding is_transformable TestUtils.is_transformable

end


text:
# Represents a pattern, a unique sequence of stops visited by a transit vehicle

struct TripPattern
    stops::Vector{Int64}
    service::Int64
end


text:
x = 2

macro tenTimesX()
  :(y = 10; x * y)
end

(function ()
  y = 3                      # Will this y be changed?
  z = @tenTimesX()           # Changed to 10 perhaps?
  @assert y == 3 && z == 20  # No, y will not change
end)()

max_stars_repo_name: idrougge/ple

max_issues_repo_name: idrougge/ple

max_forks_repo_name: idrougge/ple


text:
"""
    KMeansTreeOptions <: TreeOptions

Is the datatype that discribes which tree the [`create_tree`](@ref) function creates.

# Fields
- `iterations`: number of iterations on each level, default is one iterations
- `nchildren`: defines the number of children of each node, default is two
- `nmin`: defines the minimum amount of datapoints which are needed in a
    cluster so that it gets split up in subclusters, default is 1
- `maxlevel`: defines the maximum amount of levels, default is 100.
"""
struct KMeansTreeOptions <: TreeOptions
    iterations
    nchildren
    nmin
    maxlevel
end

function KMeansTreeOptions(;
    iterations=1,
    nchildren=2,
    nmin=1,
    maxlevel=100
)
    return KMeansTreeOptions(iterations, nchildren, nmin, maxlevel)
end

"""
    KMeansTreeNode{T} <: AbstractNode

Is the datatype of each node in the [K-Means Clustering Tree](@ref).

# Fields
- `parent::Union{KMeansTreeNode{T},Nothing}`: is the superordinate cluster of
    of the represented cluster
- `children::Union{Vector{KMeansTreeNode{T}}, Nothing}`: all directly
    subordinated clusters of the represented cluster
- `level::Integer`: the level of the represented cluster
- `center`: the center by which the cluster is defined
- `radius`: the euclidian distance between the center and the farthest away
    point
- `data::T`: array containig the indices of the points in this cluster
"""
mutable struct KMeansTreeNode{T} <: AbstractNode
    parent::Union{KMeansTreeNode{T},Nothing}
    children::Union{Vector{KMeansTreeNode{T}}, Nothing}
    level::Integer
    center
    radius
    data::T
end

function KMeansTreeNode{T}(iterations, data::T) where T
    return KMeansTreeNode{T}(nothing, nothing, 0, nothing, 0.0, data)
end

function KMeansTreeNode{T}(level, center::Vector{Float64}, radius,data::T) where T
    return KMeansTreeNode{T}(nothing, nothing, level, center, radius, data)
end

function add_child!(parent::KMeansTreeNode, center, radius, data::T) where T
    childnode = KMeansTreeNode(nothing, nothing, parent.level + 1, center, radius, data)
    childnode.parent = parent
    if parent.children === nothing
        parent.children = [childnode]
    else
        push!(parent.children, childnode)
    end
end

"""
    create_tree(points::Array{SVector{D, T}, 1}; treeoptions)

Creates an algebraic tree for an givn set of datapoints. The returned
datastructure is the foundation for the algorithms in this package.

# Arguments
- `points::Array{SVector{D, T}, 1}`: is an array of
    [SVector](https://juliaarrays.github.io/StaticArrays.jl/stable/pages/api/#SVector-1).
    Each
    [SVector](https://juliaarrays.github.io/StaticArrays.jl/stable/pages/api/#SVector-1)
    contains in general two or three float values, which discribe the position
    in the space.

# Keywords
- `treeoptions::TreeOptions`: this keyword defines by which tree is build.
    `TreeOptions` is an abstract type which either can be `BoxTreeOptions` or
    [`KMeansTreeOptions`](@ref). Default type is `BoxTreeOptions`.

# Returns
- `AbstractNode`: the root of the created tree. AbstractNode is an abstract type
    which either can be `BoxTreeNode` or [`KMeansTreeNode`](@ref), depending on the keyword.
"""
function create_tree(
    points::Vector{SVector{D,T}},
    treeoptions::KMeansTreeOptions
) where {D,T}
    root = KMeansTreeNode( nothing, nothing,  0, nothing, 0.0, Vector(1:length(points)))
    fill_tree!(
        root,
        points,
        nmin=treeoptions.nmin,
        maxlevel=treeoptions.maxlevel,
        iterations=treeoptions.iterations,
        nchildren=treeoptions.nchildren
    )

    return root
end

function whichcenter(center, point, nchildren)
    centerindex = argmin(norm.([center[i]-point for i = 1:nchildren]))

    return centerindex
end

function fill_tree!(
    node::KMeansTreeNode, points::Vector{SVector{D,T}};
    nmin = 1,
    maxlevel=log2(eps(eltype(points))),
    iterations,
    nchildren
) where {D,T}
    if length(node.data) <= nmin|| node.level >= maxlevel - 1 || length(node.data) <= nchildren
        return
    end

    center = points[node.data[1:nchildren]]
    sorted_points = zeros(eltype(node.data), length(node.data)+1, nchildren)

    for iter = 1:iterations
        sorted_points .= 0
        for pindex in node.data
            id = whichcenter(center, points[pindex], nchildren)
            sorted_points[1, id] +=  1
            sorted_points[1 + sorted_points[1, id], id] = pindex
        end
        for k=1:nchildren
            center[k] = sum(
                points[sorted_points[2:(sorted_points[1, k]+1), k]]
            ) ./ sorted_points[1,k]
        end
    end

    for i = 1:nchildren
        if sorted_points[1, i] > 0
            childpointindices = sorted_points[2:(sorted_points[1, i]+1), i]
            radius = maximum(norm.([points[childpointindices[j]].-center[i] for j = 1:length(childpointindices)]))
            add_child!(node, center[i], radius, childpointindices)
            fill_tree!(
                node.children[end],
                points,
                nmin=nmin,
                maxlevel=maxlevel,
                iterations=iterations,
                nchildren=nchildren
            )
        end
    end
end

"""
    iscompressable(sourcenode::AbstractNode, testnode::AbstractNode)

Determins whether two nodes of a tree are comressable. The criteria differs
between the [Box Tree](@ref) and the [K-Means Clustering Tree](@ref).
For the [K-Means Clustering Tree](@ref) two nodes can be compressed, if the distance
between the centers of two clusters is greater than the sum of their radius
multiplied by a factor of 1.5.
For the [Box Tree](@ref) two nodes can be compressed, if the distance between the centers
of the two boxes is greater than the sum of the distances of each box's center to
one of its corners multiplied by a factor of 1.1.

# Arguments
- `sourcenode::AbstractNode`: the node which is observed
- `testnode::AbstractNode`: the node which is tested for compression

# Returns
- `true`: if the input nodes are compressable
- `false`: if the input nodes are not compressable
"""
function iscompressable(sourcenode::KMeansTreeNode, testnode::KMeansTreeNode)
    if sourcenode.level > 0 && testnode.level > 0
        dist = norm(sourcenode.center - testnode.center)
        factor = 1.5
        if factor * (sourcenode.radius + testnode.radius) < dist
            return true
        else
            return false
        end
    else
        return false
    end
end

max_stars_repo_name: sbadrian/FastBEAST

max_issues_repo_name: sbadrian/FastBEAST

max_forks_repo_name: sbadrian/FastBEAST


text:
#
# @description Implementation of conversion helper functions.
#
# @author Alejandro Claro (alejandro.claro@gmail.com)
#
# Copyright 2017 All rights reserved.
# Use of this source code is governed by a MIT-style license that can be found in the LICENSE file.
#

# @description Unpacks a sequence of bytes into a array of bits.
#
# @param {Vector{UInt8}} data The sequence of bytes.
#
# @return {Vector{UInt8}} The array of bits.
function unpack_bits(data::Vector{UInt8})
  result = zeros(UInt8, length(data) * 8)
  index  = 1

  for x in data
  bit = 7

  while bit >= 0
      mask = (1 << bit)
      result[index] = (x & mask) != 0 ? 1 : 0
  index += 1
  bit -= 1
    end
  end

  return result
end

# @description Packs a sequence of bits into a array of bytes.
#
# @param {Vector{UInt8}} data The sequence of bits.
#
# @return {Vector{UInt8}} The array of bytes.
function pack_bits(data::Vector{UInt8})
  if length(data) % 8 != 0
    throw(ArgumentError("Could not pack bits. The bit array must be divisible by 8."))
  end

  result = zeros(UInt8, div(length(data) + 7, 8))

  for index in 1 : length(result)
    byte  = (index - 1) * 8
    value = 0

    for bit in  1 :8
      value += data[byte + bit] << (8 - bit)
    end

    result[index] = value
  end

  return result
end

# @description Permutates the given block using the specified table.
#
# @param {Vector{UInt8}} block The block to transform.
# @param {Vector{UInt8}} table The permutation table.
#
# @return {Vector{UInt8}} The permutation.
function permutate(block::Vector{UInt8}, table::Vector{UInt8})
  return [block[index] for index in table]
end

# @description Permutates the given block of bytes using the specified table of bits.
#
# @param {Vector{UInt8}} block The block to transform.
# @param {Vector{UInt8}} table The bits permutation table.
#
# @return {Vector{UInt8}} The permuted block of bytes.
function permutate_bits(block::Vector{UInt8}, table::Vector{UInt8})
  bits = unpack_bits(block)
  bits = [bits[index] for index in table]
  return pack_bits(bits)
end


text:
# This file was generated, do not modify it. # hide
fitted_params(m).best_model.max_depth


text:
# no main

void f()
{
   return;
}

max_stars_repo_name: lukaszcz/javalette

max_issues_repo_name: lukaszcz/javalette

max_forks_repo_name: lukaszcz/javalette


text:
"""
	dynamics(y, env, x, u, w)

	evaluates an environment's dynamics (in-place)

	y: state at next time step
	env: Environment
	x: state at current time step
	u: input
	w: system parameters
"""
function dynamics(y, env::Environment, x, u, w;
	gradients=false,
	attitude_decompress=false)
	step(env, x, u,
		gradients=gradients,
		attitude_decompress=attitude_decompress)[1]
    y .= env.state
end

"""
	dynamics_jacobian_state(dx, env, x, u, w)

	evaluates an environment's dynamics Jacobian wrt state (in-place)

	dx: state Jacobian at next time step
	env: Environment
	x: state at current time step
	u: input
	w: system parameters
"""
function dynamics_jacobian_state(dx, env::Environment, x, u, w;
	attitude_decompress=false)
	step(env, x, u,
		gradients=true,
		attitude_decompress=attitude_decompress)
    dx .= env.dynamics_jacobian_state
end

"""
	dynamics_jacobian_input(du, env, x, u, w)

	evaluates an environment's dynamics Jacobian wrt input (in-place)

	du: input Jacobian at next time step
	env: Environment
	x: state at current time step
	u: input
	w: system parameters
"""
function dynamics_jacobian_input(du, env::Environment, x, u, w;
	attitude_decompress=false)
	# step(env, x, u, diff=true) # this is run in dynamics_jacobian_state
	du .= env.dynamics_jacobian_input
end


text:
try
    run(`which mvn`)
catch
    error("Cannot find maven. Is it installed?")
end

cd("../jvm/sparkjl")
run(`mvn clean package`)


text:
"""
    wagner_liftingline_model(; C1=0.165, C2=0.335, eps1 = 0.0455, eps2 = 0.3)

Construct a model by coupling an unsteady aerodynamic model based on Wagner's function
(see [`wagner_model`](@ref)) and a lifting line section model (see
[`liftingline_section_model`](@ref)).
"""
function wagner_liftingline_model(; C1=0.165, C2=0.335, eps1 = 0.0455, eps2 = 0.3)

    # aerodynamic model
    aero = wagner_model(; C1, C2, eps1, eps2)

    # structural model
    stru = liftingline_section_model()

    # submodels
    submodels = (aero, stru)

    # construct coupling
    coupling = wagner_liftingline_coupling(aero, stru)

    # return the coupled model
    return CoupledModel(submodels, coupling)
end

# --- Internal Methods for this Coupling --- #

# coupling definition
function wagner_liftingline_coupling(aero, stru)

    # model constants
    C1 = aero.constants.C1
    C2 = aero.constants.C2

    # coupling function
    g = (dx, x, p, t) -> wagner_liftingline_inputs(dx, x, p, t; C1, C2)

    # number of states, inputs, and parameters (use Val(N) to use inferrable dimensions)
    nx = Val(8)
    ny = Val(9)
    np = Val(8)

    # number of parameters introduced by the coupling (use Val(N) to use inferrable dimensions)
    npc = Val(0)

    # jacobian definitions
    ratejac = Nonlinear() # TODO: define rate jacobian function
    statejac = Nonlinear() # TODO: define state jacobian function
    paramjac = Nonlinear() # TODO: define parameter jacobian function
    tgrad = Zeros()

    # return resulting coupling
    return Coupling{false}(g, nx, ny, np, npc;
        ratejac = ratejac,
        statejac = statejac,
        paramjac = paramjac,
        tgrad = tgrad)
end

# coupling function
function wagner_liftingline_inputs(dx, x, p, t; C1, C2)
    # rate variables
    dλ1, dλ2, dvx, dvy, dvz, dωx, dωy, dωz = dx
    # extract state variables
    λ1, λ2, vx, vy, vz, ωx, ωy, ωz = x
    # extract parameters
    a, b, a0, α0, cd0, cm0, ρ, c = p
    # freestream velocity components
    u, v, ω = liftingline_section_velocities(vx, vz, ωy)
    udot, vdot, ωdot = liftingline_section_accelerations(dvx, dvz, dωy)
    # calculate aerodynamic loads
    N, A, M = wagner_loads(a, b, ρ, c, a0, α0, cd0, cm0, C1, C2, u, v, ω, vdot, ωdot, λ1, λ2)
    # forces and moments per unit span
    f = SVector(A, 0, N)
    m = SVector(0, M, 0)
    # return inputs
    return vcat(u, v, ω, f, m)
end


text:
function parsedatetime(dt::String)
    res =  match(r"(\d{4})-(\d{2})-(\d{2})T(\d{2}):(\d{2}):(\d{2})\s(\w{3})",dt).captures
    datetime(int(res[1]),int(res[2]),int(res[3]),int(res[4]),int(res[5]),int(res[6]),0,eval(parse(res[7])))
end

function parsedatetime1(dt::String)
    res =  match(r"(\d{4})-(\d{2})-(\d{2})\s(\d{2}):(\d{2}):(\d{2})",dt).captures
    datetime(int(res[1]),int(res[2]),int(res[3]),int(res[4]),int(res[5]),int(res[6]),1)
end

function parsedatetime_from_TOS(dt::String)
    res =  match(r"(\d{1,2})/(\d{1,2})/(\d{2})\s(\d{2}):(\d{2}):(\d{2})",dt).captures
    datetime(2000+int(res[3]),int(res[1]),int(res[2]),int(res[4]),int(res[5]),int(res[6]))
end

function makedatetime(a::Array{Int64,1})
    datetime(a[1],a[2],a[3],a[4],a[5],a[6],a[7])
end

function datetolastsecond(datearray::Vector{Date})
#function datetolastsecond(datearray::Array{Date{ISOCalendar,UTC},1})
    res = Date{ISOCalendar, UTC}[]
    for d in 1:length(datearray)
        dtarray = [int(split(string(datearray[d]),"-")),23,59,59,59]
        push!(res, makedatetime(dtarray))
    end
    res
end

function discretesignal(ba::TimeArray{Bool,1})
    iba = float64(ba.values)
    idx = [1]
    for b in 2:length(ba)
        if ba.values[b] !== ba.values[b-1]
            push!(idx, b)
        end
    end
    TimeArray(ba.timestamp[idx], iba[idx], ["signal"])
end


text:
function backTrackingLS(f,g,maxLSiter,c1,θ,∇,fval)
    α = 1
    lineSearchIter = 0
    for i in 1:maxLSiter
        f_xk_alpha_k = f( θ + α * -∇)
        f_xk_plus_grad = fval + ( c1 * α * ∇' * -∇)[1]
		#@printf("%f \t < \t %f\n",f_xk_alpha_k,f_xk_plus_grad)
        if  ! (f_xk_alpha_k < f_xk_plus_grad[1])
            α = α * 0.5
            lineSearchIter += 1
        else
            break
        end
    end
    return (α,lineSearchIter)
end

# function / gradient-function / direction / ls iter / zoom iter / c1 / v2 / theta
function wolfeLineSearch(f,g,pk,maxLSiter,maxZoomIter,c1,c2,xk,grad,fval)
    # NEED ITERATION COUNTER
    a_0 = 0.0
    a_iminus1 = a_0
    a_i = 0.1
    a_max = 65536.0

    phi_0 = fval
    phi_a_iminus1 = phi_0
    phi_a_i = NaN

    # phi'(α) = innerProduct(pk, grad(x))
    phi_prime_0 = (pk'*grad)[1]
    phi_prime_a_i = NaN

    for i in 1:maxLSiter

        x_new = xk + a_i * pk
        phi_a_i = f(x_new)

        # check wolfe condition
        if phi_a_i > phi_0 + c1 * a_i * phi_prime_0 ||
            (i > 1 && phi_a_i >= phi_a_iminus1)
            (a_star,iter) = zoom(a_iminus1,a_i,maxZoomIter,phi_0,phi_prime_0,f,g,pk,xk,c1,c2)
            return (a_star, i , iter)
        end

        phi_prime_a_i = (pk'*g(x_new))[1]
        # check condition
        if abs(phi_prime_a_i) <= -c2*phi_prime_0
            return (a_i,i,0)
        end

        if phi_prime_a_i >= 0
            (a_star,iter) = zoom(a_i,a_iminus1,maxZoomIter,phi_0,phi_prime_0,f,g,pk,xk,c1,c2)
            return (a_star, i , iter)
        end

        # choose alphai+1 in (a_i,a_max)
        a_iminus1 = a_i
        a_i = a_i * 2

        phi_a_iminus1 = phi_a_i
    end
    return (a_i,maxLSiter,0)
end

function zoom(alpha_l, alpha_u,maxzoomiter,phi_0,phi_prime_0,f,g,pk,xk,c1,c2)

   phi_prime_al = (pk'*g(xk + alpha_l * pk))[1]
   phi_prime_au = (pk'*g(xk + alpha_u * pk))[1]

   for i in 1:maxzoomiter
        phi_al = f(xk + alpha_l*pk)
        phi_au = f(xk + alpha_u*pk)

        d1 = phi_prime_al + phi_prime_au - 3 * (phi_al - phi_au) / (alpha_l - alpha_u)
        d2 = sign(alpha_u-alpha_l)*sqrt(d1*d1-phi_prime_al*phi_prime_au)

        alpha_j = alpha_u - (alpha_u-alpha_l) * (phi_prime_au + d2 - d1) / (phi_prime_au - phi_prime_al + 2*d2)

        phi_alpha_j = f( xk + alpha_j * pk)
        # check condition
        if phi_alpha_j > phi_0 + c1 * alpha_j * phi_prime_0 || phi_alpha_j >= phi_al
            alpha_u = alpha_j
        else
            phi_prime_aj = (pk' * g(xk + alpha_j * pk))[1]
            if abs(phi_prime_aj) <= -c2*phi_prime_0
                return (alpha_j,i)
            end
            if phi_prime_aj* (alpha_u-alpha_l) >= 0
                alpha_u = alpha_l
            end
            alpha_l = alpha_j
        end

        if i == maxzoomiter
            return (alpha_j,maxzoomiter)
        end
        phi_prime_al = (pk'*g(xk + alpha_l))[1]
        phi_prime_au = (pk'*g(xk + alpha_u))[1]

    end
end

max_stars_repo_name: zagazao/num-opt

max_issues_repo_name: zagazao/num-opt

max_forks_repo_name: zagazao/num-opt


text:
title = "Home"
weight = 1
active = true


text:
viz_global_config = VegaLite.@vlfrag(legend = {orient = "top", padding = 0})
viz_defaults = (; width = 500, height = 300, frame = [-40, 0], scatter_opacity = 0.2)

estimator_color_encoding = VegaLite.@vlfrag(
    field = "estimator_name",
    type = "nominal",
    title = "Estimator",
    scale = {
        domain = ["Ours Full", "Ours Partial", "Baseline Full", "Baseline Partial"],
        # Paired color scheme:
        # range = ["#a6cee3", "#1f78b4", "#bdbdbd", "#636363"],
        range = ["#4c78a8", "#72b7b2", "#e45756", "#f58518"],
    },
)

function visualize_paramerr(;
    scatter_opacity = viz_defaults.scatter_opacity,
    width = viz_defaults.width,
    height = viz_defaults.height,
    frame = viz_defaults.frame,
    y_label = "Mean Parameter Cosine Error",
    round_x_axis = true,
)
    @vlplot(
        config = viz_global_config,
        height = height,
        width = width,
        color = estimator_color_encoding,
        x = {
            "position_observation_error:q",
            title = "Mean Absolute Postion Observation Error [m]",
            scale = {nice = round_x_axis},
        },
        transform = [
            {
                window = [
                    {field = "parameter_estimation_error", op = "median", as = "error_average"},
                    {field = "parameter_estimation_error", op = "q1", as = "error_band_lower"},
                    {field = "parameter_estimation_error", op = "q3", as = "error_band_upper"},
                ],
                groupby = ["estimator_name"],
                frame = frame,
            },
        ]
    ) +
    @vlplot(
        mark = {"point", tooltip = {content = "data"}, opacity = scatter_opacity, filled = true},
        y = {"parameter_estimation_error:q", title = y_label},
    ) +
    @vlplot(mark = "line", y = "error_average:q",) +
    @vlplot(
        mark = {"errorband", extent = "iqr"},
        y = {"error_band_lower:q", title = y_label},
        y2 = "error_band_upper:q"
    )
end

function visualize_poserr(;
    scatter_opacity = viz_defaults.scatter_opacity,
    width = viz_defaults.width,
    height = viz_defaults.height,
    frame = viz_defaults.frame,
    y_label = "Mean Absolute Position Prediction Error [m]",
    round_x_axis = true,
)
    @vlplot(
        config = viz_global_config,
        height = height,
        width = width,
        color = estimator_color_encoding,
        x = {
            "position_observation_error:q",
            title = "Mean Absolute Postion Observation Error [m]",
            scale = {nice = round_x_axis},
        },
        transform = [
            {
                window = [
                    {field = "position_estimation_error", op = "median", as = "error_average"},
                    {field = "position_estimation_error", op = "q1", as = "error_band_lower"},
                    {field = "position_estimation_error", op = "q3", as = "error_band_upper"},
                ],
                groupby = ["estimator_name"],
                frame = frame,
            },
        ]
    ) +
    @vlplot(
        mark = {"point", tooltip = {content = "data"}, opacity = scatter_opacity, filled = true},
        y = {
            "position_estimation_error:q",
            title = y_label,
            scale = {type = "symlog", constant = 0.01},
        },
        shape = {
            "converged:n",
            title = "Trajectory Reconstructable",
            legend = nothing,
            scale = {domain = [true, false], range = ["circle", "triangle-down"]},
        },
    ) +
    @vlplot(mark = {"line"}, y = "error_average:q",) +
    @vlplot(
        mark = {"errorband", extent = "iqr"},
        y = {"error_band_lower:q", title = y_label},
        y2 = "error_band_upper:q",
    )
end


text:
# missing methods (in Cairo)
function define_missing_cairo_methods()
    if !isdefined(ComplexVisual, :get_line_width)
        @eval ComplexVisual begin
            get_line_width(ctx::CairoContext) = ccall(
                (:cairo_get_line_width, Cairo.libcairo),
                Float64, (Ptr{Nothing},), ctx.ptr)
        end
    end
end
define_missing_cairo_methods()

# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4:


text:
# This file is a part of SimilaritySearch.jl

using SimilaritySearch
using Distances
using LinearAlgebra
using Test

function test_seq(db, queries, dist::SemiMetric, ksearch, valid_lower=1e-3)
    seq = ExhaustiveSearch(dist, db)
    id, dist = searchbatch(seq, queries, 10)

    for c in eachcol(dist)
        @test c[1] < valid_lower
    end
end

@testset "indexing vectors with ExhaustiveSearch" begin
    # NOTE: The following algorithms are complex enough to say we are testing it doesn't have syntax errors, a more grained test functions are required
    ksearch = 3
    db = MatrixDatabase(rand(4, 1000))
    queries = rand(db, 100)
    @info typeof(db), typeof(queries)
    for (recall_lower_bound, dist) in [
        (1.0, L2Distance()), # 1.0 -> metric, < 1.0 if dist is not a metric
        (1.0, L1Distance()),
        (1.0, LInftyDistance()),
        (1.0, SqL2Distance()),
        (1.0, LpDistance(3.0)),
        (0.1, LpDistance(0.5)),
        (1.0, AngleDistance()),
        (1.0, CosineDistance())
    ]
        test_seq(db, queries, dist, ksearch)
    end
end

@testset "indexing sequences with ExhaustiveSearch" begin
    # NOTE: The following algorithms are complex enough to say we are testing it doesn't have syntax errors, a more grained test functions are required
    ksearch = 10
    db = VectorDatabase([create_sequence(5, false) for i in 1:1000])
    queries = rand(db, 100)
    @info typeof(db), typeof(queries)

    # metric distances should achieve recall=1 (perhaps lesser because of numerical inestability)
    for dist in [
        CommonPrefixDissimilarity(),
        LevenshteinDistance(),
        LcsDistance(),
        StringHammingDistance()
    ]
        test_seq(db, queries, dist, ksearch)
    end
end

@testset "indexing sets with ExhaustiveSearch" begin
    # NOTE: The following algorithms are complex enough to say we are testing it doesn't have syntax errors, a more grained test functions are required
    ksearch = 10
    db = VectorDatabase([create_sequence(5, true) for i in 1:1000])
    queries = rand(db, 100)
    @info typeof(db), typeof(queries)

    # metric distances should achieve recall=1 (perhaps lesser because of numerical inestability)
    for dist in [
        JaccardDistance(),
        DiceDistance(),
        IntersectionDissimilarity()
    ]
        test_seq(db, queries, dist, ksearch)
    end
end

@testset "Normalized Cosine and Normalized Angle distances" begin
    # cosine and angle distance
    ksearch = 10
    X = MatrixDatabase(rand(4, 1000))
    queries = rand(X, 100)
    normalize!.(X)

    test_seq(X, queries, NormalizedAngleDistance(), ksearch)
    test_seq(X, queries, NormalizedCosineDistance(), ksearch)
end

@testset "Binary hamming distance" begin
    ksearch = 10
    db = MatrixDatabase(rand(UInt64, 8, 1000))
    queries = rand(db, 100)
    test_seq(db, queries, BinaryHammingDistance(), ksearch)
end

max_stars_repo_name: sadit/SimilaritySearch

max_issues_repo_name: sadit/NNS

max_forks_repo_name: sadit/SimilaritySearch


text:
mutable struct SumFun <: Function# OscLaurent{D<:SPeriodicLine,R}?
    funs::Vector{Fun}
end

evaluate(f::SumFun,x) = sum(map(y -> y(x),f.funs))

SumFun(f::Fun{T}) where T = SumFun([f])

conj(s::SumFun) = SumFun(map(conj,s.funs))
conj(s::Array{T}) where T<:SumFun = map(conj,s)

sum(s::SumFun) = sum(map(sum,s.funs))
sum(s::Array{T}) where T<:SumFun = map(sum,s)

combine!(s::SumFun) = SumFun(combine!(s.funs))
combine!(s::Vector{SumFun}) = map(combine!,s)
function combine!(f::Vector{T}) where T <: Fun
    i = 0
    while i < length(f)
        i += 1
        j = i
        while j < length(f)
            j += 1
            if spacescompatible(f[i],f[j])
                f[i] += f[j]
                deleteat!(f,j)
                j -= 1
            end
        end
    end
    f
end

+(f::SumFun,g::SumFun) = combine!(SumFun(vcat(g.funs,f.funs)))
-(f::SumFun,g::SumFun) = combine!(SumFun(vcat(-g.funs,f.funs)))

+(f::SumFun,g::Fun) = f + SumFun(g)
+(g::Fun,f::SumFun) = f + SumFun(g)

-(f::SumFun,g::Fun) = f - SumFun(g)
-(g::Fun,f::SumFun) = -f + SumFun(g)

*(f::Fun,g::SumFun) = SumFun(map(y -> f*y,g.funs))
*(f::SumFun,g::Fun) = SumFun(map(y -> g*y,f.funs))

*(f::SumFun,g::SumFun) = combine!(+(map(y -> y*g,f.funs)...))

*(c::Float64,g::SumFun) = SumFun(map(y -> c*y,g.funs))
*(c::Complex{Float64},g::SumFun) = SumFun(map(y -> c*y,g.funs))

(f::SumFun)(x) = evaluate(f,x)
+(f::SumFun) = f

(f::Array)(x) = map(y -> evaluate(y,x),f)
chop!(f::SumFun) = SumFun(map(chop!,f.funs))
chop!(f::Array) = map(chop!,f)

chop!(f::SumFun,tol::Float64) = SumFun(map(x -> chop!(x,tol),f.funs))
chop!(f::Array,tol::Float64) = map(x -> chop!(x,tol),f)

SumFun(f::Fun{T}) where T<: SumSpace = SumFun(components(f))
copy(f::SumFun) = SumFun(copy(f.funs))


text:
include("population.jl")

function geneticAlgorithm!(data, elitePercentage, crossoverPoint, mutationPercentage, repeatSize)
     calculatePopulationFitness!(data)
     populationSize = length(data)
     popGen = 0
     bestFit = data[1].fitness
     repeatNum = 1
     eliteSize = Int(trunc(elitePercentage*populationSize))
      # Ako nam za ukrstanje ostaje neparan broj jedinki, povecamo eliteSize za 1
     eliteSize = eliteSize+(populationSize-eliteSize)%2
     while (data[1].fitness>0) && (repeatNum<repeatSize)
        elite = deepcopy(selectPopulation(data, eliteSize))
        data = selectPopulation(data, populationSize-eliteSize)
        data = crossoverPopulation(data, crossoverPoint)
        mutatePopulation!(data, mutationPercentage)
        data = [data; elite]
        calculatePopulationFitness!(data)
        popGen = popGen+1
        print("Generacija $popGen bestFit $bestFit broj ponavljanja $repeatNum \n")
        if abs(bestFit-data[1].fitness)<0.01
           repeatNum = repeatNum + 1
        else
           bestFit = data[1].fitness
           repeatNum = 1
	end
     end
     return popGen, repeatNum, data
end

max_stars_repo_name: pikacxe/Primenjeni-Algoritmi

max_issues_repo_name: pikacxe/Primenjeni-Algoritmi

max_forks_repo_name: pikacxe/Primenjeni-Algoritmi


text:
module GymWrappers

using OpenAIGym
import OpenAIGym: AbstractGymEnv, actionset
import Reinforce: reward, total_reward

abstract type AbstractGymWrapper <: AbstractGymEnv end

export AbstractGymEnv, AbstractGymWrapper, gymenv, getwrapper, maybe_wrap
# other exports in included files

# Defaults
for fn in [:state, :reset!, :finished, :actions, :step!, :reward, :total_reward]
    @eval Reinforce.$fn(wrapped_env::AbstractGymWrapper, args...; kwargs...) =
            Reinforce.$fn(wrapped_env.env, args...; kwargs...)
end

# for fn in [:render]
#     @eval OpenAIGym.$fn(wrapped_env::AbstractGymWrapper; kwargs...) =
#             OpenAIGym.$fn(wrapped_env.env; kwargs...)
# end

# With multiple wraps, eventually this will get to the actual GymEnv
getwrapper{WrapType <: AbstractGymEnv}(::Type{WrapType}, wrpenv::AbstractGymWrapper) =
    getwrapper(WrapType, wrpenv.env) # recurses down the wrappers

# If the env is of the right type, return it
getwrapper{WrapType <: AbstractGymWrapper}(::Type{WrapType}, env::WrapType) = env

# Special case when looking for a GymEnv, to win the method specialisation race
getwrapper(::Type{GymEnv}, env::GymEnv) = env

# If we're here we looked through all the wrappers and there weren't any `WrapType` wrappers
getwrapper{WrapType <: AbstractGymEnv}(::Type{WrapType}, env::GymEnv) =
    error("env is not wrapped with a wrapper of type $WrapType")

# short syntax for getting the gym env
gymenv(env::AbstractGymEnv) = getwrapper(GymEnv, env)

# The Business
"""
`maybe_wrap(env::AbstractGymEnv; wrapper_specs=default_wrapper_specs, wrap_deepmind=true)`

Wrap an environment using the wrapper_specs in `wrapper_specs[gamename]`, where
`gamename` is "Breakout" for `env.name == "BreakoutNoFrameskip-v4"`.

If `!haskey(wrapper_specs, gamename)` then the `env` is returned unwrapped.
Example wrapper spec:

```
default_wrapper_specs = Dict(
    "Pong"=>[ActionSetWrapper=>([2, 3],), GreyMeanWrapper=>(false,)]
)
```

Which will call `ActionSetWrapper(env, [2,3])`, `GreyMeanWrapper(env, false)`
"""
function maybe_wrap(env::AbstractGymEnv; wrapper_specs=default_wrapper_specs,
                    wrap_deepmind=true)
    gamename, gamever = env.name, env.ver
    replace(gamename, r"(NoFrameskip|Deterministic)","")
    wrapped_env = env
    all_wrappers = get(wrapper_specs, gamename, [])
    wrap_deepmind && append!(all_wrappers, wrapper_specs["deepmind_defaults"])
    for wrapper_spec in all_wrappers
        wrapper, args = wrapper_spec
        wrapped_env = wrapper(wrapped_env, args...)
    end
    wrapped_env
end

include("misc_utils.jl")
include("reward_wrappers.jl")
include("observation_wrappers.jl")
include("action_wrappers.jl")
include("episode_wrappers.jl")

# --------------------------------
# Default Wrappers
# --------------------------------
default_wrapper_specs = Dict{String, Vector{Pair{Type{T} where T <: AbstractGymWrapper, Tuple}}}(
    "deepmind_defaults"=>[
        # wrappers applied from inner most (top) -> outer most (bottom)
        EpisodicLifeWrapper=>(),
        GreyMeanWrapper=>(),
        # GreyChanWrapper=>(2,),
        DownsizeWrapper=>(84, 84), # (210, 160) -> (84, 84)
        MaxAndSkipWrapper=>(4,), # repeat action for 4 frames
        MultiFrameWrapper=>(4,), # stack last 4 frames
    ]
)

# --------------------------------
# Tests
# --------------------------------
#---
# envname = "Pong-v0"
# env = GymEnv(envname) |> maybe_wrap
# #---
# gryimg = OpenAIGym.render(env; mode="rgb_array")
# #---
# colorview(Gray, (env.env.env.state[:,:,3] ./ 255)) == gryimg
# #---
# imshow(gryimg)
# #---
# OpenAIGym.render(env) # should call imshow automatically
# #---
# actions(env)
#---
end


text:
"""
A `FunctionSpace{S,T}` is the supertype of all function spaces. The space
contains functions that map type `S` to type `T`.
Thus, `S` is the domaintype and `T` is the codomain type.
"""
abstract type FunctionSpace{S,T}
end

domaintype(::Type{FunctionSpace{S,T}}) where {S,T} = S
domaintype(::Type{S}) where {S <: FunctionSpace} = domaintype(supertype(S))
domaintype(::FunctionSpace{S,T}) where {S,T} = S

codomaintype(::Type{FunctionSpace{S,T}}) where {S,T} = T
codomaintype(::Type{S}) where {S <: FunctionSpace} = codomaintype(supertype(S))
codomaintype(::FunctionSpace{S,T}) where {S,T} = T

length(s::FunctionSpace) = Inf

show(io::IO, s::FunctionSpace) = print(io, "Space : $(domaintype(s)) → $(codomaintype(s))")
Display.object_parentheses(s::FunctionSpace) = true

"""
`GenericFunctionSpace` is a generic function space for functions about which we
have no additional information apart from their domain and codomain types.
"""
struct GenericFunctionSpace{S,T} <: FunctionSpace{S,T}
end

# The zero element of a generic function space
zero(space::GenericFunctionSpace{S,T}) where {S,T} = x::S -> zero(T)

# The multiplicative identity element of a generic function space
one(space::GenericFunctionSpace{S,T}) where {S,T} = x::S -> one(T)

"""
A function space equipped with a measure.

The measure `μ` induces a bilinear form `(f,g) = int(f, g, dμ)` and possibly a
(semi)norm or inner product. This depends on the properties of the measure.
"""
struct MeasureSpace{M <: Measure,S,T} <: FunctionSpace{S,T}
    measure     ::  M

    # Ensure that the domain types of the space and the measure are the same
    MeasureSpace{M,S,T}(measure::Weight{S}) where {M,S,T} = new(measure)
end

MeasureSpace(measure::Weight{T}) where {T} =
    MeasureSpace{typeof(measure),T,codomaintype(measure)}(measure)

measure(s::MeasureSpace) = s.measure

domain(s::MeasureSpace) = support(measure(s))

space(m::Weight) = MeasureSpace(m)

show(io::IO, s::MeasureSpace) = print(io, "Function space with measure: $(measure(s))")
Display.object_parentheses(s::MeasureSpace) = true

const FourierSpace{T} = MeasureSpace{FourierWeight{T},T,T}

FourierSpace() = FourierSpace{Float64}()
FourierSpace{T}() where {T} = MeasureSpace(FourierWeight{T}())
show(io::IO, s::FourierSpace) = print(io, "Fourier space : $(domaintype(s)) → $(codomaintype(s))")

const ChebyshevTSpace{T} = MeasureSpace{ChebyshevTWeight{T},T,T}

ChebyshevTSpace() = ChebyshevTSpace{Float64}()
ChebyshevTSpace{T}() where {T} = MeasureSpace(ChebyshevTWeight{T}())

struct L2{S,T} <: FunctionSpace{S,T}
    domain  ::  Domain{S}
end

L2() = L2{Float64}()
L2{T}() where {T} = L2{T}(DomainSets.FullSpace{T}())

L2(domain::Domain{S}) where {S} = L2{subeltype(S)}(domain)
L2{T}(domain::Domain{S}) where {S,T} = L2{S,T}(domain)

domain(space::L2) = space.domain

abstract type WeightedL2{S,T} <: FunctionSpace{S,T} end

struct ChebyshevSpace{S,T} <: WeightedL2{S,T}
end

ChebyshevSpace() = ChebyshevSpace{Float64}()
ChebyshevSpace{T}() where {T} = ChebyshevSpace{T,T}()

domain(space::ChebyshevSpace{S,T}) where {S,T} = ChebyshevInterval{S}()

struct ProductSpace{S,T} <: FunctionSpace{S,T}
    spaces
end

function tensorproduct(spaces::FunctionSpace...)
    S = Tuple{map(domaintype, spaces)...}
    T = mapreduce(codomaintype, promote_type, spaces)
    ProductSpace{S,T}(spaces)
end

components(space::ProductSpace) = space.spaces


text:
# test local projection stabilization

function test_lps()

  opts = Dict{String, Any}(
    "physics" => "Euler",
    "operator_type" => "SBPDiagonalE",
    "dimensions" => 2,
    "run_type" => 5,
    "jac_method" => 2,
    "jac_type" => 2,
    "order" => 3,
    "IC_name" => "ICIsentropicVortex",
    "use_DG" => true,
    "volume_integral_type" => 2,
    "Volume_flux_name" => "IRFlux",
    "Flux_name" => "IRSLFFlux",
    "numBC" => 3,
    "BC1" => [0],
    "BC1_name" => "isentropicVortexBC",  # outlet
    "BC2" => [2],
    "BC2_name" => "isentropicVortexBC", # inlet
    "BC3" => [1, 3],
    "BC3_name" => "noPenetrationBC",  # was noPenetrationBC
    "aoa" => 0.0,
    "smb_name" => "SRCMESHES/vortex_3x3_.smb",
    "dmg_name" => ".null",
    "solve" => false,
    )

  mesh, sbp, eqn, opts = solvePDE(opts)

  @testset "LP Stabiization" begin
    # test the reference element
    degree = sbp.degree
    q = zeros(mesh.numNodesPerElement)
    coords = calcnodes(sbp, sbp.vtx)

    # the projection operator should not modify polynomials up to degree
    # p, and should decrease the energy of polynomials of higher degrees

    for d=0:(degree + 1)
      for i=1:mesh.numNodesPerElement
        x = coords[1, i]; y = coords[2, i]
        q[i] = x^(d) + 2*y^(d)
      end

      q2 = eqn.params.lps_data.P*q
      e1 = dot(q, sbp.w .* q)
      e2 = dot(q2, sbp.w .* q2)  # energy of higher order modes

      if d <= degree
        @test abs(e2) < 1e-15
      else
        @test e2 < 0.8*e1  # 0,8 is arbitrary factor
      end
    end   # end d

    # now test on non-reference elements
    for d=0:(degree + 1)
      for i=1:mesh.numEl
        for j=1:mesh.numNodesPerElement
          x = mesh.coords[1, j, i]; y = mesh.coords[2, j, i]
          q[j] = x^(d) + 2*y^(d)
        end

        q2 = eqn.params.lps_data.P*q
        e1 = dot(q, sbp.w .* q)
        e2 = dot(q2, sbp.w .* q2)  # energy of higher order modes

        if d <= degree
          @test abs(e2) < 1e-15
        else
          @test e2 < 0.8*e1  # 0,8 is arbitrary factor
        end
      end  # end i
    end  # end d

    # check entropy stability
    eqn.q .+= 0.01*rand(size(eqn.q))
    w = zeros(mesh.numDofPerNode, mesh.numNodesPerElement, mesh.numEl)
    entropy_vars = eqn.params.lps_data.entropy_vars
    for i=1:mesh.numEl
      for j=1:mesh.numNodesPerElement
        q_i = ro_sview(eqn.q, :, j, i)
        w_i = sview(w, :, j, i)
        EulerEquationMod.convertToEntropy(entropy_vars, eqn.params, q_i, w_i)
      end
    end

    fill!(eqn.res, 0)
    EulerEquationMod.applyLPStab(mesh, sbp, eqn, opts)
    val = sum(w .* eqn.res)
    println("val = ", val)
    @test val < 0
  end  # end testset

  return nothing
end

add_func1!(EulerTests, test_lps, [TAG_SHORTTEST])


text:
"""
    enclose(f::Function, dom::Interval_or_IntervalBox,
            solver::Symbol=:IntervalArithmetic; [kwargs]...)

Return a range enclosure of a univariate or multivariate function on the given
domain.

### Input

- `f`      -- function
- `dom`    -- hyperrectangular domain, either a unidimensional  `Interval` or
              a multidimensional `IntervalBox`
- `solver` -- (optional, default: `IntervalArithmetic`) choose one among the
              available solvers; see `RangeEnclosures.available_solvers` for the
              full list
- `kwargs` -- optional keyword arguments passed to the solver; for available
              options see the documentation of each solver

### Output

An interval representing the range enclosure (minimum and maximum) of `f` over
its domain `dom`.

### Examples

```jldoctest enclose_examples
julia> using RangeEnclosures

julia> enclose(x -> 1 - x^4 + x^5, 0..1) # use default solver
[0, 2]

julia> enclose(x -> 1 - x^4 + x^5, 0..1, :IntervalArithmetic)
[0, 2]

julia> enclose(x -> 1 - x^4 + x^5, 0..1, :TaylorModels, order=4)
[0.78125, 1.125]

julia> enclose(x -> 1 - x^4 + x^5, 0..1, :TaylorModels, order=10)
[0.8125, 1.09375]

julia> enclose(x -> 1 - x^4 + x^5, 0..1, :IntervalOptimisation)
[0.916034, 1.00213]
```
You can also try other solvers such as `SumOfSquares` and `AffineArithmetic`.

A vector of solvers can be passed in the `solver` options. Then, the result is
obtained by intersecting the range enclosure of each solver.
In the previous example,

```jldoctest enclose_examples
julia> using RangeEnclosures

julia> enclose(x -> 1 - x^4 + x^5, 0..1, [:TaylorModels, :IntervalArithmetic])
[0.8125, 1.09375]
```
"""
function enclose(f::Function, dom::Interval_or_IntervalBox,
                 solver::Symbol=:IntervalArithmetic; kwargs...)

    𝑂 = Dict(kwargs)
    numvars = length(dom)

    if solver == :IntervalArithmetic
        # solver
        R = enclose_IntervalArithmetic(f, dom)

    elseif solver == :BranchandBound
        tol =  :tol ∈ keys(𝑂) ? 𝑂[:tol] : 0.6
        order = :order ∈ keys(𝑂) ? 𝑂[:order] : 10
        #solver
        R = enclose_BranchandBound(f, dom, order=order, tol=tol)

    elseif solver == :IntervalOptimisation
        # unpack options or set defaults
        structure =  :structure ∈ keys(𝑂) ? 𝑂[:structure] : HeapedVector
        tol =  :tol ∈ keys(𝑂) ? 𝑂[:tol] : 1e-3

        # solver
        R = enclose_IntervalOptimisation(f, dom, structure=structure, tol=tol)

    elseif solver == :TaylorModels
        # unpack options or set defaults
        normalize =  :normalize ∈ keys(𝑂) ? 𝑂[:normalize] : true
        order = :order ∈ keys(𝑂) ? 𝑂[:order] : 10

        # solver
        R = enclose_TaylorModels(f, dom; order=order, normalize=normalize)

    elseif solver == :SumOfSquares
        # unpack options or set defaults
        if :order ∈ keys(𝑂)
            order = 𝑂[:order]
            pop!(𝑂, :order)
        else
            order = 5
        end

        if :backend ∈ keys(𝑂)
            backend = 𝑂[:backend]
            pop!(𝑂, :backend)
        else
            backend = SDPA.Optimizer
        end

        R = enclose_SumOfSquares(f, dom; order=order, backend=backend, 𝑂...)

    elseif solver == :AffineArithmetic
        # requires affine arithmetic to be loaded
        @assert isdefined(@__MODULE__, :AffineArithmetic) "package `AffineArithmetic` " *
        "not loaded but it is required for executing `enclose`"

        # solver
        R = enclose_AffineArithmetic(f, dom)

    else
        error("algorithm $algorithm unknown")
    end

    return R
end

function enclose(p::AbstractPolynomialLike, dom::Interval_or_IntervalBox,
                 solver=:IntervalArithmetic; kwargs...)
    f(x...) = p(variables(p) => x)
    return enclose(f, dom, solver; kwargs...)
end

function enclose(f::Function, dom::Interval_or_IntervalBox,
                 method::Vector{Symbol}; kwargs...)
   return mapreduce(ξ -> enclose(f, dom, ξ, kwargs...), ∩, method)
end

"""
    relative_precision(x::Interval{N}, xref::Interval{N}) where {N}

Return the relative precision of an interval representing a range enclosure
given a reference interval.

### Input

- `x`    -- test interval
- `xref` -- reference interval

### Output

An interval containing the left (resp. right) percentages that describe the
relative precision of `x` with respect to the reference interval `xref`.

### Examples

Suppose that the reference interval is ``[-1.2, 4.6]``, and let ``[-1.25, 7.45]``
be the overapproximation of the reference interval obtained with some global
optimisation method that returns a range enclosure.
The relative precision is then:

```jldoctest relative_precision
julia> using RangeEnclosures: relative_precision

julia> xref = Interval(-1.2, 4.6)
[-1.2, 4.6]

julia> x = Interval(-1.25, 7.45)
[-1.25, 7.45001]

julia> relative_precision(x, xref)
[0.862068, 49.138]
```
Notice how the percentage of relative error of the maximum is big, close to
``50%``, while the approximation of the minimum is much smaller, close to ``1%``
of the reference value.

### Algorithm

This function measures the relative precision of the result in a more informative
way than taking the scalar overestimation because it evaluates the precision of
the lower and the upper range bounds separately (cf. Eq. (20) in [1]).

[1] Althoff, Matthias, Dmitry Grebenyuk, and Niklas Kochdumper.
   *Implementation of Taylor models in CORA 2018.*
   Proc. of the 5th International Workshop on Applied Verification for
   Continuous and Hybrid Systems. 2018.

### Notes

It is assumed that the test, or paragon interval, is an overapproximation of the
reference interval; otherwise the left or right percentages will be negative.
"""
function relative_precision(x::Interval{N}, xref::Interval{N}) where {N}
    x⁻, x⁺ = inf(x), sup(x)
    xref⁻, xref⁺ = inf(xref), sup(xref)

    Δref = xref⁺ - xref⁻
    @assert Δref != zero(N) "the width of the reference interval is $Δref, but " *
    "it should have a non-zero width"

    rel⁻ = -(x⁻ - xref⁻) / Δref
    rel⁺ = (x⁺ - xref⁺) / Δref
    return 100 * Interval(rel⁻, rel⁺)
end


text:
export McmcInitialiser, InitMode, InitRandEdit, InitRandIns, get_init
export perturb, perturb_ins

# Initilisers
# -----------

# These can be passed to mcmc samplers to determine the defualt initialisation scheme.

"""
Abstract type representing initialisation schemes for SIS model samplers.
"""
abstract type McmcInitialiser end

"""
`InitMode <: McmcInitialiser` - this is a MCMC initialisation scheme for SIS model samplers which starts the MCMC chain at the model mode by default.
"""
struct InitMode <: McmcInitialiser
    function InitMode()
        return new()
    end
end

function get_init(
    initiliaser::InitMode,
    model::Union{SIS,SIM},
    )
    return model.mode
end

struct InitRandEdit <: McmcInitialiser
    δ::Int
    function InitRandEdit(δ::Int)
        return new(δ)
    end
end

function perturb(
    S::InteractionSequence,
    V::AbstractArray,
    K_inner::DimensionRange,
    δ::Int
    )

    S_init = deepcopy(S)
    N = length(S)

    ind_del = zeros(Int, δ)
    ind_add = zeros(Int, δ)
    vals = zeros(Int, δ)

    rem_edits = δ

    for i in 1:N
        if i == N
            δ_tmp = rem_edits
        else
            p = 1/(N-i+1)
            δ_tmp = rand(Binomial(rem_edits, p))
        end

        if δ_tmp == 0
            continue
        else

            n = length(S[i])
            d = rand(max(0, ceil(Int, (n + δ_tmp + - K_inner.u)/2)):min(n-K_inner.l, δ_tmp))
            m = n + δ_tmp - 2*d

            ind_del_v = view(ind_del, 1:d)
            ind_add_v = view(ind_add, 1:(δ_tmp-d))
            vals_v = view(vals, 1:(δ_tmp-d))

            StatsBase.seqsample_a!(1:n, ind_del_v)
            StatsBase.seqsample_a!(1:m, ind_add_v)
            sample!(V, vals)

            delete_insert!(S_init[i], ind_del_v, ind_add_v, vals_v)

        end

        rem_edits -= δ_tmp

        if rem_edits == 0
            break
        end
    end

    return S_init

end

function get_init(
    initialiser::InitRandEdit,
    model::Union{SIS,SIM}
    )
    return perturb(model.mode, model.V, model.K_inner, initialiser.δ)
end

struct InitRandIns <: McmcInitialiser
    δ::Int
    function InitRandIns(δ::Int)
        return new(δ)
    end
end

function perturb_ins(
    S::InteractionSequence,
    V::AbstractArray,
    K_inner::DimensionRange,
    δ::Int
    )

    S_init = deepcopy(S)
    N = length(S)

    ind_add = zeros(Int, δ)

    rem_edits = δ

    for i in 1:N
        if i == N
            δ_tmp = rem_edits
        else
            p = 1/(N-i+1)
            δ_tmp = rand(Binomial(rem_edits, p))
        end

        if δ_tmp == 0
            continue
        else

            n = length(S[i])
            m = n + δ_tmp
            if m > K_inner.u
                δ_tmp = (K_inner.u-n)
                m = n + δ_tmp
            end
            ind_add_v = view(ind_add, 1:δ_tmp)
            StatsBase.seqsample_a!(1:m, ind_add_v)

            for j in ind_add_v
                insert!(S_init[i], j, rand(V))
            end

        end

        rem_edits -= δ_tmp

        if rem_edits == 0
            break
        end
    end

    return S_init

end

function get_init(
    initialiser::InitRandIns,
    model::Union{SIS,SIM}
    )
    return perturb_ins(model.mode, model.V, model.K_inner, initialiser.δ)
end


text:
using UnalignedVectors
using Base.Test

ENDIAN_BOM == 0x04030201 || error("tests implemented only for little endian machines")

@testset "UnalignedVectors" begin
    @test_throws DimensionMismatch UnalignedVector{UInt16}([0x01, 0x02, 0x03])
    @test_throws ErrorException UnalignedVector{String}([0x01, 0x02])
    @test_throws MethodError UnalignedVector{UInt16}([1, 2])
    @test UnalignedVector{UInt16}([0x01, 0x02]) == [0x0201]

    a = collect(0x01:0x14)
    v = UnalignedVector{UInt16}(a)
    @test IndexStyle(v) == IndexLinear()
    @test length(v) == 10
    @test size(v) == (10,)
    for i = 1:10
        @test v[i] == (2i)<<8 + 2i - 1
    end
    @test_throws BoundsError v[0]
    @test_throws BoundsError v[11]
    for i = 1:10
        v[i] = i
    end
    @test v[1:10] == [v[i] for i in 1:10] == [i for i in 1:10]
    @test v == collect(0x0001:0x000a)
    @test a == vec(vcat((0x01:0x0a)', zeros(UInt8,1,10)))
    @test_throws InexactError v[1] = 0xffffffff

    @test unaligned_reinterpret(UInt8, a) === a
    v = unaligned_reinterpret(Int16, a)
    @test v == collect(1:10)
    r = reinterpret(Int32, v, (1, 5))
    @test size(r) == (1, 5)
    for i = 1:5
        @test r[1,i] == 2i*2^16 + 2i - 1
    end
    r2 = reinterpret(UInt16, r, (10,))
    @test isa(r2, UnalignedVector{UInt16})
    @test isa(reinterpret(UInt8, r, (20,)), Vector{UInt8})
end

@testset "Mmap" begin
    fn = tempname()
    A = [1.0 2.0;
         3.0 4.0]
    open(fn, "w") do io
        write(io, "BIGARRAY\n")
        write(io, ndims(A))
        for s in size(A)
            write(io, s)
        end
        write(io, A)
    end
    function reader(io)
        String(read(io, 9)) == "BIGARRAY\n" || error("file is not a BIGARRAY file")
        n = read(io, Int)
        sz = (read(io, Int, n)...)
        a = Mmap.mmap(io, Vector{UInt8}, sizeof(Float64)*prod(sz), position(io))
        v = UnalignedVector{Float64}(a)
        return reshape(v, sz), a
    end
    B, mmapped = open(fn) do io
        reader(io)
    end
    @test B == A
    # Clean up. Some platforms (e.g., Windows) may have a hard time
    # with this, since the file was mmapped.
    finalize(mmapped)
    gc(true)
    sleep(1.0)
    rm(fn)
end

nothing


text:
using Test, Statistics, Random
using CUDA: CUDA, functional
using Knet.KnetArrays: KnetArray
using AutoGrad: Param, @gcheck

# This is buggy in julia (missing dims arg) as of Sep 16, 2019 so fix it here:
# https://github.com/JuliaLang/julia/issues/33280
Statistics.stdm(A::AbstractArray, m; corrected::Bool=true, dims=:) =
    ((dims === :) ? sqrt.(varm(A, m; corrected=corrected)) : sqrt.(varm(A, m; corrected=corrected, dims=dims)))

if CUDA.functional(); @testset "statistics" begin
    for T in (Float32,Float64)
        a = randn(T,3,4)
        k = KnetArray(a)
        p = Param(k)
        @test mean(a) ≈ mean(k)
        @test @gcheck mean(p)
        @test mean(a,dims=1) ≈ mean(k,dims=1)
        @test @gcheck mean(p,dims=1)
        @test mean(a,dims=2) ≈ mean(k,dims=2)
        @test @gcheck mean(p,dims=2)
        @test mean(abs,a) ≈ mean(abs,k)
        @test @gcheck mean(abs,p)
        @test mean(abs2,a) ≈ mean(abs2,k)
        @test @gcheck mean(abs2,p)
        @test std(a) ≈ std(k)
        @test @gcheck std(p)
        @test std(a,dims=1) ≈ std(k,dims=1)
        @test @gcheck std(p,dims=1)
        @test std(a,dims=2) ≈ std(k,dims=2)
        @test @gcheck std(p,dims=2)
        @test stdm(a,mean(a)) ≈ stdm(k,mean(k))
        @test @gcheck stdm(p,mean(p))
        @test stdm(a,mean(a,dims=1),dims=1) ≈ stdm(k,mean(k,dims=1),dims=1)
        @test @gcheck stdm(p,mean(p,dims=1),dims=1)
        @test stdm(a,mean(a,dims=2),dims=2) ≈ stdm(k,mean(k,dims=2),dims=2)
        @test @gcheck stdm(p,mean(p,dims=2),dims=2)
        @test var(a) ≈ var(k)
        @test @gcheck var(p)
        @test var(a,dims=1) ≈ var(k,dims=1)
        @test @gcheck var(p,dims=1)
        @test var(a,dims=2) ≈ var(k,dims=2)
        @test @gcheck var(p,dims=2)
        @test varm(a,mean(a)) ≈ varm(k,mean(k))
        @test @gcheck varm(p,mean(p))
        @test varm(a,mean(a,dims=1),dims=1) ≈ varm(k,mean(k,dims=1),dims=1)
        @test @gcheck varm(p,mean(p,dims=1),dims=1)
        @test varm(a,mean(a,dims=2),dims=2) ≈ varm(k,mean(k,dims=2),dims=2)
        @test @gcheck varm(p,mean(p,dims=2),dims=2)
    end
end; end


text:
# MIT License

#Copyright (c) 2021 hondoRandale <jules.rasetaharison@tutanota.com> and contributors

#Permission is hereby granted, free of charge, to any person obtaining a copy
#of this software and associated documentation files (the "Software"), to deal
#in the Software without restriction, including without limitation the rights
#to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
#copies of the Software, and to permit persons to whom the Software is
#furnished to do so, subject to the following conditions:

#The above copyright notice and this permission notice shall be included in all
#copies or substantial portions of the Software.

#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
#IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
#FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
#AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
#LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
#OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
#SOFTWARE.

include( "benchmarkOptimProblems.jl" )

using SlimeMoldOptim
using Test

@testset "SlimeMoldOptim.jl" begin
    # 2-D objective functions
    @test Beale_2D( [3.0f0,0.5f0]; kwargs=() )                     ≈ 0.0f0
    @test typeof( Beale_2D( [3.0f0,0.5f0]; kwargs=() ) )          == Float32

    @test Himmelblau_2D( [ 3.0f0, 2.0f0 ]; kwargs=() )             ≈ 0.0f0
    @test typeof( Himmelblau_2D( [ 3.0f0, 2.0f0 ]; kwargs=() ) )  == Float32
    @test Himmelblau_2D( [ -2.805118f0,  3.131312f0 ]; kwargs=() ) ≈ 0.0f0 atol=1f-3
    @test Himmelblau_2D( [ -3.779310f0, -3.283186f0 ]; kwargs=() ) ≈ 0.0f0 atol=1f-10
    @test Himmelblau_2D( [ 3.584428f0, -1.848126f0 ]; kwargs=() )  ≈ 0.0f0 atol=1f-5

    @test Levi_no13_2D( [ 1.0f0, 1.0f0 ]; kwargs=() )                      ≈ 0.0f0 atol=1f-15
    @test typeof( Levi_no13_2D( [ 1.0f0, 1.0f0 ]; kwargs=() ) )           == Float32

    @test Ackley_2D( [ 0.0f0, 0.0f0 ]; kwargs=() )                         ≈ 0.0f0
    @test typeof( Ackley_2D( [ 0.0f0, 0.0f0 ]; kwargs=() ) )              == Float32

    @test Booth_2D(  [ 1.0f0, 3.0f0 ]; kwargs=() )                         ≈ 0.0f0
    @test typeof( Booth_2D(  [ 1.0f0, 3.0f0 ]; kwargs=() ) )              == Float32

    @test Rosenbrock_2D( [ 1.0f0, 1.0f0 ]; kwargs=() )                     ≈ 0.0f0
    @test typeof( Rosenbrock_2D( [ 1.0f0, 1.0f0 ]; kwargs=() ) )          == Float32

    @test Matyas_2D( [ 0.0f0, 0.0f0 ]; kwargs=() )                         ≈ 0.0f0
    @test typeof( Matyas_2D( [ 0.0f0, 0.0f0 ]; kwargs=() ) )              == Float32

    @test Three_hump_camel_2D( [ 0.0f0, 0.0f0 ]; kwargs=() )               ≈ 0.0f0
    @test typeof( Three_hump_camel_2D( [ 0.0f0, 0.0f0 ]; kwargs=() ) )    == Float32

    @test Easom_2D( [ Float32(π), Float32(π) ]; kwargs=() )                ≈ -1.0f0
    @test typeof( Easom_2D( [ Float32(π), Float32(π) ]; kwargs=() ) )     == Float32

    @test Sphere_2D( [0.0f0, 0.0f0]; kwargs=() )                           ≈ 0.0f0
    @test typeof( Sphere_2D( [0.0f0, 0.0f0]; kwargs=() ) )                == Float32

    @test Schaffer_function_no2_2D( [0.0f0, 0.0f0]; kwargs=() )            ≈ 0.0f0
    @test typeof( Schaffer_function_no2_2D( [0.0f0, 0.0f0]; kwargs=() ) ) == Float32

    @test Bukin_function_no6_2D( [-10.0f0, 1.0f0 ]; kwargs=() )            ≈ 0.0f0
    @test typeof( Bukin_function_no6_2D( [-10.0f0, 1.0f0 ]; kwargs=() ) ) == Float32

    # N-D objective functions
    @test Michalewicz_ND( [ 2.20f0, 1.57f0 ]; kwargs=() )                          ≈ -1.8081267f0
    @test typeof( Michalewicz_ND( [ 2.20f0, 1.57f0 ]; kwargs=() ) )               == Float32

    @test Zakharov_ND( [ 0.0f0, 0.0f0 ]; kwargs=() )                               ≈ 0.0f0
    @test typeof( Zakharov_ND( [ 0.0f0, 0.0f0 ]; kwargs=() ) )                    == Float32

    @test AxisParallelHyperEll_ND( [ 0.0f0, 0.0f0 ]; kwargs=() )                   ≈ 0.0f0
    @test typeof( AxisParallelHyperEll_ND( [ 0.0f0, 0.0f0 ]; kwargs=() ) )        == Float32

    @test Sphere_ND( [ 0.0f0, 0.0f0 ]; kwargs=() )                                 ≈ 0.0f0
    @test typeof( Sphere_ND( [ 0.0f0, 0.0f0 ]; kwargs=() ) )                      == Float32

    @test Rastrigin_ND( [ 0.0f0, 0.0f0 ]; kwargs=() )                              ≈ 0.0f0
    @test typeof( Rastrigin_ND( [ 0.0f0, 0.0f0 ]; kwargs=() ) )                   == Float32

    @test Styblinski_Tang_ND( [ -2.903534f0, -2.903534f0 ]; kwargs=() )            ≈ -78.33233f0
    @test typeof( Styblinski_Tang_ND( [ -2.903534f0, -2.903534f0 ]; kwargs=() ) ) == Float32

    @test Exponential_ND( [ 0.0f0, 0.0f0 ]; kwargs=() )                            ≈ -1.0f0
    @test typeof( Exponential_ND( [ 0.0f0, 0.0f0 ]; kwargs=() ) )                 == Float32

    @test Griewank_ND( [ 0.0f0, 0.0f0 ]; kwargs=() )                               ≈ 0.0f0
    @test typeof( Griewank_ND( [ 0.0f0, 0.0f0 ]; kwargs=() ) )                    == Float32

    lower    = [ -10.0f0, -10.0f0 ];
    upper    = [ 10.0f0, 10.0f0 ];

    sma( objFunction ) = slimeMoldAlgo( 30, 1000, ();
                                        lower       = lower,
                                        upper       = upper,
                                        objFunction = objFunction,
                                        maxiter     = 1000,
                                        ϵ_conv      = 1f-20,
                                        z           = 0.03f0 )

   @test sum( [ sma( fObj ).minFGlobal for fObj ∈ obj_functions_2D ] ) < 0.0f0
   @test sum( [ sma( fObj ).minFGlobal for fObj ∈ obj_functions_Nd ] ) < 750.0f0
end


text:
using Test

testcases = filter(f -> endswith(f, ".jl") & (f != "runtests.jl"), readdir())

for testcase in testcases
    @testset "$(testcase)" begin
        include(joinpath(testcase))
    end
end

max_stars_repo_name: crinders/BayesianSumProductNetworks

max_issues_repo_name: crinders/BayesianSumProductNetworks

max_forks_repo_name: crinders/BayesianSumProductNetworks


text:
println("\n\nmake_d_datestamp: start")

prepend!(LOAD_PATH, ["Project.toml"])
import LibPQ

DISASTER_PASSWORD = read("./etc/disaster-pass", String)

postgres = LibPQ.Connection("host=127.0.0.1 port=5432 dbname=disaster user=disaster password=$DISASTER_PASSWORD")

sql =
"""
drop table if exists public.d_datestamp;
create table if not exists public.d_datestamp (
    d_datestamp_id  date primary key,
    year          integer not null,
    month         integer not null,
    day           integer not null,
    weekofyear    integer not null
);

do \$body\$
declare
    _date  date;
begin
    for _date in select generate_series('1992-01-01'::date, '2015-12-31', '1 day') loop
        -- raise notice '_date=%', _date;
        insert into public.d_datestamp (
            d_datestamp_id,
            year,
            month,
            day,
            weekofyear
        ) values (
            _date,
            date_part('year', _date),
            date_part('month', _date),
            date_part('day', _date),
            date_part('week', _date)
        ) on conflict do nothing;
    end loop;

end; \$body\$ language plpgsql;
"""
LibPQ.execute(postgres, sql)
LibPQ.close(postgres)

println("make_d_datestamp: stop")


text:
# Extend Impute.jl methods to work across multiple components at once
# Since Impute.jl only supports single dimension operations and handling subselection
# filtering would be tricky we aren't going to handle merging `dims` with other
# pattern keys for now.

"""
    Impute.validate(ds::KeyedDataset, validator::Validator; dims=:)

Apply the validator to components in the [`KeyedDataset`](@ref) with the specified `dims`.
"""
function Impute.validate(ds::KeyedDataset, validator::Validator; dims=:)
    patterns = dims === Colon() ? Pattern[(:__, )] : Pattern[(:__, dims)]
    return map(a -> Impute.validate(a, validator; dims=dims), ds, patterns...)
end

"""
    Impute.apply(ds::KeyedDataset, imp::DeclareMissings)

Declare missing values across all components in the [`KeyedDataset`](@ref).
"""
Impute.apply(ds::KeyedDataset, imp::DeclareMissings) = Impute.apply!(deepcopy(ds), imp)
function Impute.apply!(ds::KeyedDataset, imp::DeclareMissings)
    return map!(a -> Impute.apply!(a, imp), ds, ds)
end

"""
    Impute.impute(ds, imp; dims)

Apply the imputation algorithm `imp` along the `dims` for all components of the
[`KeyedDataset`](@ref) with that dimension.

# Example
```jldoctest
julia> using AxisKeys, Impute; using AxisSets: KeyedDataset, flatten;

julia> ds = KeyedDataset(
           flatten([
               :train => [
                   :temp => KeyedArray([1.0 1.1; missing 2.2; 3.0 3.3]; time=1:3, id=[:a, :b]),
                   :load => KeyedArray([7.0 7.7; 8.0 missing; 9.0 9.9]; time=1:3, loc=[:x, :y]),
                ],
                :predict => [
                   :temp => KeyedArray([1.0 missing; 2.0 2.2; 3.0 3.3]; time=1:3, id=[:a, :b]),
                   :load => KeyedArray([7.0 7.7; 8.1 missing; 9.0 9.9]; time=1:3, loc=[:x, :y]),
                ]
            ])...
       );

julia> [k => parent(parent(v)) for (k, v) in Impute.substitute(ds; dims=:time).data]  # KeyedArray printing isn't consistent in jldoctests
4-element Vector{Pair{Tuple{Symbol, Symbol}, Matrix{Union{Missing, Float64}}}}:
   (:train, :temp) => [1.0 1.1; 2.2 2.2; 3.0 3.3]
   (:train, :load) => [7.0 7.7; 8.0 8.0; 9.0 9.9]
 (:predict, :temp) => [1.0 1.0; 2.0 2.2; 3.0 3.3]
 (:predict, :load) => [7.0 7.7; 8.1 8.1; 9.0 9.9]

julia> [k => parent(parent(v)) for (k, v) in Impute.substitute(ds; dims=:loc).data]
4-element Vector{Pair{Tuple{Symbol, Symbol}, Matrix{Union{Missing, Float64}}}}:
   (:train, :temp) => [1.0 1.1; missing 2.2; 3.0 3.3]
   (:train, :load) => [7.0 7.7; 8.0 8.8; 9.0 9.9]
 (:predict, :temp) => [1.0 missing; 2.0 2.2; 3.0 3.3]
 (:predict, :load) => [7.0 7.7; 8.1 8.8; 9.0 9.9]
```
"""
function Impute.impute(ds::KeyedDataset, imp::Imputor; dims, kwargs...)
    return Impute.impute!(deepcopy(ds), imp; dims=dims, kwargs...)
end

function Impute.impute!(ds::KeyedDataset, imp::Imputor; dims, kwargs...)
    # NOTE: We don't use `map!` here because we know we can mutate the underlying data,
    # so we don't need to assign the value back into the `ds.data`
    selection = filter(p -> dims in dimnames(last(p)), ds.data)

    for (k, v) in selection
        Impute.impute!(v, imp; dims=dims, kwargs...)
    end

    # Just to be safe, we check that our imputation method didn't break our constraints.
    validate(ds)
    return ds
end

"""
    Impute.apply(ds, filter; dims)

Filter out missing data along the `dims` for each component in the [`KeyedDataset`](@ref)
with that dimension.

# Example
```jldoctest
julia> using AxisKeys, Impute; using AxisSets: KeyedDataset, Pattern, flatten;

julia> ds = KeyedDataset(
           flatten([
               :train => [
                   :temp => KeyedArray([1.0 1.1; missing 2.2; 3.0 3.3]; time=1:3, id=[:a, :b]),
                   :load => KeyedArray([7.0 7.7; 8.0 missing; 9.0 9.9]; time=1:3, loc=[:x, :y]),
                ],
                :predict => [
                   :temp => KeyedArray([1.0 missing; 2.0 2.2; 3.0 3.3]; time=1:3, id=[:a, :b]),
                   :load => KeyedArray([7.0 7.7; 8.1 missing; 9.0 9.9]; time=1:3, loc=[:x, :y]),
                ]
            ])...
       );

julia> [k => parent(parent(v)) for (k, v) in Impute.filter(ds; dims=:time).data]  # KeyedArray printing isn't consistent in jldoctests
4-element Vector{Pair{Tuple{Symbol, Symbol}, Matrix{Union{Missing, Float64}}}}:
   (:train, :temp) => [3.0 3.3]
   (:train, :load) => [9.0 9.9]
 (:predict, :temp) => [3.0 3.3]
 (:predict, :load) => [9.0 9.9]

julia> [k => parent(parent(v)) for (k, v) in Impute.filter(ds; dims=Pattern(:train, :__, :time)).data]
4-element Vector{Pair{Tuple{Symbol, Symbol}, Matrix{Union{Missing, Float64}}}}:
   (:train, :temp) => [1.0 1.1; 3.0 3.3]
   (:train, :load) => [7.0 7.7; 9.0 9.9]
 (:predict, :temp) => [1.0 missing; 3.0 3.3]
 (:predict, :load) => [7.0 7.7; 9.0 9.9]

julia> [k => parent(parent(v)) for (k, v) in Impute.filter(ds; dims=:loc).data]
4-element Vector{Pair{Tuple{Symbol, Symbol}, Matrix{Union{Missing, Float64}}}}:
   (:train, :temp) => [1.0 1.1; missing 2.2; 3.0 3.3]
   (:train, :load) => [7.0; 8.0; 9.0]
 (:predict, :temp) => [1.0 missing; 2.0 2.2; 3.0 3.3]
 (:predict, :load) => [7.0; 8.1; 9.0]
```
"""
Impute.apply(ds::KeyedDataset, f::Filter; dims) = Impute.apply!(deepcopy(ds), f; dims=dims)

function Impute.apply!(ds::KeyedDataset, f::Filter; dims)
    pattern = _pattern(dims)
    dim = pattern.segments[end]

    dim in (:_, :__) && throw(ArgumentError(
        "$pattern points to an ambiguous dimension in the dataset. " *
        "Kwarg `dims` must end with a dimname."
    ))
    filter_paths = dimpaths(ds, pattern)

    # Limit our constraint map to paths containing the supplied dim
    cmap = filter(constraintmap(ds)) do (constraint, paths)
        any(p -> p in filter_paths, paths)
    end

    # Add remaining unconstrained components
    constraint_paths = union(values(cmap))
    for p in filter(!in(constraint_paths), filter_paths)
        cmap[p] =  Set([p])
    end

    # Extract component paths
    filter_paths = [p[1:end-1] for p in filter_paths]

    # Apply our shared filter mask for each set of paths
    for (constraint, paths) in cmap
        @debug "$constraint => $paths"
        # We're assuming this dataset has already been validated so all dimpaths are
        # already equal
        mask = trues(length(axiskeys(ds, first(paths))))

        # Pre-extract our component keys and values
        selection = [p[1:end-1] => ds.data[p[1:end-1]] for p in paths]

        # First pass to determine our shared key mask
        for (k, v) in selection
            # Only filter components defined by dims
            if k in filter_paths
                for (i, s) in enumerate(eachslice(v; dims=dim))
                    mask[i] &= f.func(s)
                end
            end
        end

        # Second pass to use selectdim on each component with our mask
        for (k, v) in selection
            # copy is so we don't change the data element type to a view
            ds.data[k] = copy(selectdim(v, NamedDims.dim(dimnames(v), dim), mask))
        end
    end

    return ds
end


text:
abstract type AbstractKoopmanAlgorithm end;

"""
    DMDPINV()

Approximates the Koopman operator `K` based on

```julia
K = Y / X
```

where `Y` and `X` are data matrices.

"""
mutable struct DMDPINV <: AbstractKoopmanAlgorithm end;

# Fast but more allocations
(x::DMDPINV)(X::AbstractArray, Y::AbstractArray) = Y / X

"""
    DMDSVD(rtol)

Approximates the Koopman operator `K` based on the singular value decomposition
of `X` such that:

```julia
K = Y*V*Σ*U'
```

where `Y` and `X = U*Σ*V'` are data matrices.
If `rtol` ∈ (0, 1) is given, the singular value decomposition is reduced to include only
entries bigger than `rtol*maximum(Σ)`. If `rtol` is an integer, the reduced SVD up to `rtol` is used
for computation.
"""
mutable struct DMDSVD{T} <: AbstractKoopmanAlgorithm
    rtol::T
end;

DMDSVD() = DMDSVD(0.0)

# Slower but fewer allocations
function (x::DMDSVD)(X::AbstractArray, Y::AbstractArray)
    U, S, V = svd(X)
    if typeof(x.rtol) <: Int && x.rtol < size(X, 1)
        V = V[:, x.rtol]
    elseif zero(x.rtol) < x.rtol < one(x.rtol)
        r = vec(S .> x.rtol*maximum(S))
        U = U[:, r]
        S = S[r]
        V = V[:, r]
    end

    return Y*V*Diagonal(one(eltype(X)) ./ S)*U'
end

"""
    TOTALDMD(rtol, alg)

Approximates the Koopman operator `K` with the algorithm `alg` over the rank-reduced data
matrices `Xᵣ = X Qᵣ` and `Yᵣ = Y Qᵣ`, where `Qᵣ` originates from the singular value decomposition of
the joint data `Z = [X; Y]`. Based on [this paper](http://cwrowley.princeton.edu/papers/Hemati-2017a.pdf).

If `rtol` ∈ (0, 1) is given, the singular value decomposition is reduced to include only
entries bigger than `rtol*maximum(Σ)`. If `rtol` is an integer, the reduced SVD up to `rtol` is used
for computation.
"""
mutable struct TOTALDMD{R, A} <: AbstractKoopmanAlgorithm
    rtol::R
    alg::A
end

TOTALDMD() = TOTALDMD(0.0, DMDPINV())

function (x::TOTALDMD)(X::AbstractArray, Y::AbstractArray)
    _ , S, Q = svd([X; Y])
    if typeof(x.rtol) <: Int && x.rtol < size(X, 1)
        Q = Q[:, 1:x.rtol]
    elseif zero(x.rtol) < x.rtol < one(x.rtol)
        r = vec(S .> x.rtol*maximum(S))
        Q = Q[:, r]
    end

    return x.alg(X*Q, Y*Q)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ d9f9bd70-508d-11ec-1d5d-b3592d653992
# ╠═╡ show_logs = false
begin
	let
		using Pkg
		Pkg.activate(mktempdir())
		Pkg.Registry.update()
		Pkg.add("PlutoUI")
		Pkg.add(url="https://github.com/Dale-Black/DistanceTransforms.jl")
		Pkg.add(url="https://github.com/Dale-Black/Losers.jl")
	end

	using PlutoUI
	using DistanceTransforms
	using Losers
end

# ╔═╡ 5398d829-ad81-4971-94b5-f87eeaaa381e
md"""
## Import packages
"""

# ╔═╡ 44119707-2397-4e7f-83f7-f362eae6004f
TableOfContents()

# ╔═╡ 6e03547c-d9b6-4972-9574-7cb4f30948f9
md"""
## Distance transforms + deep learning
One area of recent interest in deep learning has been distance transform-based loss functions to directly reduce Hausdorff loss and other boundary-specific loss metrics [D. Karimi and S. E. Salcudean](http://arxiv.org/abs/1904.10030).

This package can be utilized directly in various training loops, along with [Losers.jl](https://github.com/Dale-Black/Losers.jl) to train segmentation networks, like the common UNet architecture, to directly reduce Hausdorff distance using the `hausdorff()` loss function.

In practice, it is common to combine `hausdorff` with the `dice` loss function to help stabalize during training.

Below we will show how one might set this up, without actually training on anything.
"""

# ╔═╡ 2810025c-9940-4da2-a587-938c425f511e
y1 = [1 1 1 0; 1 1 1 0; 1 1 1 0; 1 1 1 0]

# ╔═╡ dfc81bb5-29fa-46f2-a335-66f541205d22
y = cat(y1, y1; dims=3)

# ╔═╡ 89347efa-ee13-4c61-a4e6-3d4d1a77955d
ŷ = copy(y)

# ╔═╡ 96d6cb3f-a500-490b-9bc9-500731af2651
md"""
Given some fake three dimensional array, lets apply a distance transform to this array for use in the `hausdorff` loss function
"""

# ╔═╡ 646d3c5b-e72b-49ac-b1e1-74d2d89c57c2
y_dtm = euclidean(y)

# ╔═╡ 10f64f57-e052-4895-9fe3-58e1c4fb2d88
ŷ_dtm = euclidean(ŷ)

# ╔═╡ 6e52eae8-e1cd-4d07-a000-d591b7bfa3ff
hausdorff(ŷ, y, ŷ_dtm, y_dtm)

# ╔═╡ b7c539d1-0e65-4716-ad23-cf4400bf4bf7
md"""
As one would expect, the difference between two identical arrays is zero
"""

# ╔═╡ 2780485d-ce45-4d5f-a88c-d4823bcf4273
begin
	y2 = rand([0, 1], 10, 10, 10)
	ŷ2 = rand([0, 1], 10, 10, 10)

	y2_dtm = euclidean(y2)
	ŷ2_dtm = euclidean(ŷ2)
end

# ╔═╡ 423704a8-926e-4fd3-98c4-7da8374a35a5
hausdorff(ŷ2, y2, ŷ2_dtm, y2_dtm)

# ╔═╡ 5d1ac21a-8073-4817-869b-5500649e486f
md"""
## GPU considerations
When using distance transform-based loss functions for deep learning, one must consider GPUs for many real-world deep learning tasks. The benefit of DistanceTransforms.jl is that GPU compatible distance transform algorithms are ready to use, out-of-the-box.

We will look at how one would go about this, without actually calling a CUDA array from the CUDA.jl library.
"""

# ╔═╡ cfd9c96d-1ce7-4edc-954c-4d844aff6629
begin
	y3 = rand([0, 1], 10, 10, 10)
	ŷ3 = rand([0, 1], 10, 10, 10)
end

# ╔═╡ 6345df57-dba7-4608-9e08-9c9d981810e5
begin
	y3_bool = boolean_indicator(y3)
	ŷ3_bool = boolean_indicator(ŷ3)
end

# ╔═╡ dde2505d-1a69-4201-a75a-e29a62197163
begin
	tfm1 = SquaredEuclidean(y3_bool)
	tfm2 = SquaredEuclidean(ŷ3_bool)
end

# ╔═╡ d2a987a3-f44c-4239-9e96-d83ae16c265e
begin
	y3_dtm = transform(y3_bool, tfm1)
	ŷ3_dtm = transform(ŷ3_bool, tfm1)
end

# ╔═╡ 90467d7a-785d-44d3-92bc-9d05bcb5f5b4
hausdorff(ŷ3, y3, ŷ3_dtm, y3_dtm)

# ╔═╡ 43e694a4-bbd5-4072-a4de-0ea65b16c906
md"""
If we were to be running an actual training loop it might look something like this...

```julia
# Begin training loop
for epoch in 1:max_epochs
	step = 0
    α = ...
	@info epoch

	# Loop through training data
	for (xs, ys) in train_loader
		step += 1
		@info step

		# Send data to GPU
		xs, ys = xs |> gpu, ys |> gpu
		gs = Flux.gradient(ps) do
			ŷs = model(xs)

			# Apply distance transform using GPU compatible `SquaredEuclidean`
			# Data will usually be 4D or 5D [x, y, (z), channel, batch]

			ys_dtm = CuArray{Float32}(undef, size(ys))
			ŷs_dtm = CuArray{Float32}(undef, size(ŷs))
			for b in size(ys, 5)
				for c in size(ys, 4)
					for z in size(ys, 3)
						bool_arr_gt = boolean_indicator(ys[:, :, z, c, b])
						bool_arr_pred = boolean_indicator(ŷs[:, :, z, c, b])
						tfm_gt = SquaredEuclidean(bool_arr_gt)
						tfm_pred = SquaredEuclidean(bool_arr_pred)
						ys_dtm[:, :, z, c, b] = transform!(bool_arr_gt, tfm_gt, threads)
						ŷs_dtm[:, :, z, c, b] = transform!(bool_arr_pred, tfm_pred, threads)
					end
				end
			end
			hd_loss = hausdorff(ŷs, ys, ŷs_dtm, ys_dtm)
			dice_loss = dice(ŷs, ys)
			loss = α*dice_loss + (1-α)*hausdorff_loss
			return loss
		end
		Flux.update!(optimizer, ps, gs)
	end
end
```
"""

# ╔═╡ Cell order:
# ╟─5398d829-ad81-4971-94b5-f87eeaaa381e
# ╠═d9f9bd70-508d-11ec-1d5d-b3592d653992
# ╠═44119707-2397-4e7f-83f7-f362eae6004f
# ╟─6e03547c-d9b6-4972-9574-7cb4f30948f9
# ╠═2810025c-9940-4da2-a587-938c425f511e
# ╠═dfc81bb5-29fa-46f2-a335-66f541205d22
# ╠═89347efa-ee13-4c61-a4e6-3d4d1a77955d
# ╟─96d6cb3f-a500-490b-9bc9-500731af2651
# ╠═646d3c5b-e72b-49ac-b1e1-74d2d89c57c2
# ╠═10f64f57-e052-4895-9fe3-58e1c4fb2d88
# ╠═6e52eae8-e1cd-4d07-a000-d591b7bfa3ff
# ╟─b7c539d1-0e65-4716-ad23-cf4400bf4bf7
# ╠═2780485d-ce45-4d5f-a88c-d4823bcf4273
# ╠═423704a8-926e-4fd3-98c4-7da8374a35a5
# ╟─5d1ac21a-8073-4817-869b-5500649e486f
# ╠═cfd9c96d-1ce7-4edc-954c-4d844aff6629
# ╠═6345df57-dba7-4608-9e08-9c9d981810e5
# ╠═dde2505d-1a69-4201-a75a-e29a62197163
# ╠═d2a987a3-f44c-4239-9e96-d83ae16c265e
# ╠═90467d7a-785d-44d3-92bc-9d05bcb5f5b4
# ╟─43e694a4-bbd5-4072-a4de-0ea65b16c906

max_stars_repo_name: Dale-Black/DistanceTransformsDocs

max_issues_repo_name: Dale-Black/DistanceTransformsDocs

max_forks_repo_name: Dale-Black/DistanceTransformsDocs


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

"""
	CVortex

A GPU accelerated vortex particle and vortex filament library.
"""
module CVortex

	export 	particle_induced_velocity,
			filament_induced_velocity,
			particle_induced_dvort,
			filament_induced_dvort,
			particle_visc_induced_dvort,
			particle_field_vorticity,
			filament_induced_velocity_influence_matrix,
			RegularisationFunction,
			singular_regularisation,
			planetary_regularisation,
			gaussian_regularisation,
			winckelmans_regularisation,
			redistribute_particles_on_grid,
			particle_pedrizzetti_relaxation,
			RedistributionFunction,
			lambda0_redistribution,
			lambda1_redistribution,
			lambda2_redistribution,
			lambda3_redistribution,
			m4p_redistribution,
			number_of_accelerators,
			number_of_enabled_accelerators,
			accelerator_name,
			accelerator_enabled,
			accelerator_enable,
			accelerator_disable
	#-------------------------------------------------------------------------
	# Loading shared library binary (cvortex)

	import Libdl: dlopen

	const libcvortex = joinpath(
		dirname(dirname(@__FILE__)), "deps/libcvortex")
	function __init__()
		try
			dlopen(libcvortex)
		catch
			error("$(libcvortex) cannot be opened. Consider "*
				"rebuilding the package.")
		end
		# Inialisation is automatic:
		ccall(
			("cvtx_initialise", libcvortex),
			Cvoid, ())
	end

	#-------------------------------------------------------------------------
	include("ConvertionChecks.jl")
	include("Vec3f.jl")
	include("Vec2f.jl")
	include("Accelerators.jl")
	include("RegularisationFunction.jl")
	include("RedistributionFunction.jl")
	include("VortexParticle3D.jl")
	include("VortexParticle2D.jl")
	include("VortexFilament.jl")
	include("CVortexInfo.jl")

end #module