text:
function inverse(side::HorizonSide)
    return HorizonSide(mod(Int(side) + 2, 4))
end

function f()
    marked = 0
    Temperature = 0
    side = Ost
    while !isborder(robot, Nord)
        while !isborder(robot, side)
            move!(robot, side)
            if ismarker(robot)
                Temperature += temperature(robot)
                marked += 1
            end
        end
        side = inverse(side)
        move!(robot, Nord)
    end
    while !isborder(robot, side)
        move!(robot, side)
        if ismarker(robot)
            Temperature += temperature(robot)
            marked += 1
        end
    end
    if marked
        println(0)
    else
        println(Temperature / marked)
    end
end

f()

max_stars_repo_name: TheXlebushek/mirea-progs

max_issues_repo_name: TheXlebushek/mirea-progs

max_forks_repo_name: TheXlebushek/mirea-progs


text:
module Info

project_name = "cytof-diff-dens"

awsbucket = "s3://$(project_name)"
resultsdir = joinpath(ENV["SCRATCH_DIR"], project_name)

awsbucket_datastudy = joinpath(awsbucket, "datastudy")
resultsdir_datastudy = joinpath(resultsdir, "datastudy")

awsbucket_simstudy = joinpath(awsbucket, "simstudy")
resultsdir_simstudy = joinpath(resultsdir, "simstudy")

awsbucket_compare = joinpath(awsbucket, "compare")
resultsdir_compare = joinpath(resultsdir, "compare")

end  # module

ext: jl


text:
function update_mean(μ::Real, n::Integer, x::Real; remove::Bool=false)
    mask = 2 * Int(!remove) - 1
    new_n = n + mask
    if new_n == 0
        return 0.0, 0
    elseif new_n < 0
        error("Cannot remove an element from an empty collection")
    else
        sum_x = μ * n
        new_sum_x = sum_x + mask * x
        new_μ = new_sum_x / new_n
        return new_μ, new_n
    end
end

function update_mean_std(μ::Real, σ::Real, n::Integer, x::Real; remove::Bool=false)
    new_μ, new_n = update_mean(μ, n, x; remove=remove)
    if new_n == 0
        return 0.0, 0.0, 0
    end
    σ² = σ^2
    μ_x² = σ² + μ^2
    new_μ_x², new_n = update_mean(μ_x², n, x^2; remove=remove)
    new_σ² = new_μ_x² - new_μ^2
    # if new_σ² < -0.1
    #     @warn "precision error?" new_σ²
    # end
    new_σ² = max(0, new_σ²)
    new_σ = √(new_σ²)
    return new_μ, new_σ, new_n
end

function update_mean_std_corr(μ_x::Real, μ_y::Real, σ_x::Real, σ_y::Real, ρ::Real, n::Integer, x::Real, y::Real; remove::Bool=false)
    new_μ_x, new_σ_x, new_n = update_mean_std(μ_x, σ_x, n, x; remove=remove)
    if new_n == 0
        return 0.0, 0.0, 0.0, 0.0, 0.0, 0
    end
    new_μ_y, new_σ_y, new_n = update_mean_std(μ_y, σ_y, n, y; remove=remove)
    cov = ρ * σ_x * σ_y
    μ_xy = cov + μ_x * μ_y
    new_μ_xy, new_n = update_mean(μ_xy, n, x * y; remove=remove)
    new_cov = new_μ_xy - new_μ_x * new_μ_y
    if abs(new_cov) < 1e-9
        new_ρ = 0.0
    else
        new_ρ = new_cov / (new_σ_x * new_σ_y)
    end
    # if abs(new_ρ) > 1 + 1e-9
    #     @warn "precision error?" new_ρ
    # end
    new_ρ = min(max(-1, new_ρ), 1)
    return new_μ_x, new_μ_y, new_σ_x, new_σ_y, new_ρ, new_n
end

# stats = (0.0, 0.0, 0.0, 0.0, 0.0, 0)

# stats = update_mean_std_corr(stats..., 1, 1)
# stats = update_mean_std_corr(stats..., 2, 1.4)
# stats = update_mean_std_corr(stats..., 1.4, 1.5)

function get_trace_distribution(w_trace, possible_weights)
    dist = Dict{Float64,Float64}()
    for w in possible_weights
        dist[round(w, digits=2)] = 0
    end
    for w in w_trace
        dist[round(w, digits=2)] = dist[round(w, digits=2)] + 1 / length(w_trace)
    end
    return dist
end


text:
struct NestedIntervalContext{T₁ <: AbstractControl,   T₂ <: AbstractModel,  T₃ <: AbstractDisplayProperties} <: AbstractPlotContext
    control::T₁
    model::T₂
    display_properties::T₃
end

function (context::NestedIntervalContext{<: NestedIntervalControl,   <: NestedInterval,  <: NestedIntervalDisplayProperties})()
        control = context.control
        model = context.model
        display_properties = context.display_properties
        id = get_id(display_properties)
        caption = get_caption(display_properties)
        is_new_window(display_properties) ? CImGui.Begin(caption*id,C_NULL, CImGui.ImGuiWindowFlags_NoBringToFrontOnFocus) : nothing
        isenabled(control) ? control(model, display_properties) : nothing
        is_new_window(display_properties) ? CImGui.End() : nothing
end


text:
@testset "Distances" begin
  p = P2(0, 1)
  l = Line(P2(0, 0), P2(1, 0))
  @test evaluate(Euclidean(), p, l) == T(1)
  @test evaluate(Euclidean(), l, p) == T(1)

  p1, p2 = P2(1, 0), P2(0, 1)
  @test evaluate(Chebyshev(), p1, p2) == T(1)

  p = P2(68, 259)
  l = Line(P2(68, 260), P2(69, 261))
  @test evaluate(Euclidean(), p, l) ≤ T(0.8)
end


text:
@testset "Properties" begin

  atoms = readPDB("./structure.pdb","protein")

  @test mass(atoms) ≈ 11079.704440000156
  f = formula(select(atoms,"residue < 5"))
  @test f.formula == [("H", 25), ("C", 19), ("N", 4), ("O", 7), ("S", 1)]

  @test element.(select(atoms,"residue = 1")) == ["N", "H", "H", "H", "C", "H",
                                                  "C", "H", "H", "H", "C", "O" ]
  @test element_name.(select(atoms,"residue = 1")) == [ "Nitrogen", "Hydrogen", "Hydrogen",
                                                        "Hydrogen", "Carbon", "Hydrogen",
                                                        "Carbon", "Hydrogen", "Hydrogen",
                                                        "Hydrogen", "Carbon", "Oxygen" ]
  @test atomic_number.(select(atoms,"residue = 1")) == [ 7, 1, 1, 1, 6, 1, 6, 1,
                                                         1, 1, 6, 8 ]
  residues = collect(eachresidue(atoms))
  @test length(residues) == 104
  @test name(Residue(atoms,1:12)) == "ALA"
  @test Residue(atoms,1:12).range == 1:12

  m = maxmin(atoms)
  @test m.xmin ≈ [-14.18, -17.561, -15.369]
  @test m.xmax ≈ [18.694, 14.182, 15.909]
  @test m.xlength ≈ [32.873999999999995, 31.743000000000002, 31.278]

  s = select(atoms,"residue = 3")
  @test coor(s) == [-4.383 -11.903 -6.849; -4.51 -11.263 -6.096;
                    -3.903 -11.262 -8.062; -3.731 -12.076 -8.767;
                    -4.938 -10.279 -8.612; -4.417 -9.552 -9.06;
                    -5.543 -9.911 -7.784; -5.867 -10.85 -9.684;
                    -5.451 -10.837 -10.863; -6.974 -11.289 -9.3;
                    -2.626 -10.48 -7.749; -1.94 -10.014 -8.658]

  s2 = select(atoms,"residue = 5")
  @test distance(s, s2) ≈ 3.6750402718881863
  x1 = coor(s); x2 = coor(s2);
  @test distance(x1, x2) ≈ 3.6750402718881863
  x1 = coor(s,xyz_in_cols=false); x2 = coor(s2,xyz_in_cols=false);
  @test distance(x1, x2, xyz_in_cols=false) ≈ 3.6750402718881863
  @test distance(residues[3],residues[5]) ≈ 3.6750402718881863

  r = Residue(select(atoms,"residue = 3"))
  @test coor(s) == coor(r)
  @test coor(select(atoms,"residue = 3")) == coor(residues[3])
  @test coor(residues[3])' == coor(residues[3],xyz_in_cols=false)
  @test PDBTools.same_residue(atoms[1],atoms[2]) == true
  @test PDBTools.same_residue(atoms[1],atoms[20]) == false

  @test residuename("glu") == "Glutamic acid"
  @test oneletter("Glu") == "E"
  @test threeletter("E") == "GLU"

end


text:
# This file is part of FlagSets.jl project, which is licensed under BDS-3-Clause.
# See file LICENSE.md for more information.

"""
    @symbol_flagset T[::BaseType] flag_1[=bit_1] flag_2[=bit_2] ...

Create a `FlagSet{Symbol,BaseType}` subtype with name `T` and the flags supplied.
Each instance of `T` represents a subset of `{:flag1, :flag2, ...}`.
The associated `bit_1`,`bit_2`, etc, if provided, must be distinct positive powers of 2.

The following constructors are provided:
    - `T(itr)`: construct a set of the symbols from the given iterable object. Each symbol
    must be one of :`flag_1`, :`flag_2`, etc;
    - `T()`: construct an empty set;
    - `T(sym1::Symbol, syms::Symbol...)`: equivalent to `T([sym1, syms...])`
    - `T(bitmask::Integer)`: construct a set of flags associated with each bit in `mask`;
    - `T(; kwargs...)`: each `flag_i` can be supplied as a boolean value
    (`flag_i = true` includes `:flag_i`, `flag_i = false` excludes it).

It should be safe to new constructor methods and/or overriding the methods provided
(except the default using an internal struct, i.e., `T(::FlagSets.BitMask)`).

To construct sets of non-`Symbol` constants, use the macro [`@flagset`](@ref).

# Examples
```julia
julia> @symbol_flagset FontFlags bold=1 italic=2 large=4

julia> FontFlags(1)
FontFlags with 1 element:
  :bold

julia> FontFlags(:bold, :italic)
FontFlags with 2 elements:
  :bold
  :italic

julia> FontFlags(bold=true, italic=false)
FontFlags with 1 element:
  :bold

julia> for flag in FontFlags(3); println(flag) end
bold
italic

julia> FontFlags(3) | FontFlags(4)
FontFlags with 3 elements:
  :bold
  :italic
  :large

julia> eltype(FontFlags)
Symbol

julia> typemax(FontFlags)
FontFlags with 3 elements:
  :bold
  :italic
  :large

julia> typemin(FontFlags)
FontFlags([])
```

Values can also be specified inside a `begin` block, e.g.

```julia
@symbol_flagset T begin
    flag_1 [= x]
    flag_2 [= y]
end
```

If `BaseType` is provided, it must be an `Integer` type big enough to represent each
`bit_i`. Otherwise, some convenient integer type is inferred
(`UInt32` is used if possible, to guarantee compatibility with C code).
Each instance of `T` can be converted to `BaseType`, where each bit set in the
converted `BaseType` corresponds to a flag in the flag set. `read` and `write`
perform these conversions automatically. In case the flagset is created with a non-default
`BaseType`, `Integer(flagset)` will return the integer `flagset` with the type `BaseType`.

```julia
julia> convert(Int, FontFlags(bold=true, italic=false))
1

julia> convert(Integer, FontFlags(bold=true, italic=false))
0x00000001
```

To list flags of a flag set type or instance, use `flags`, e.g.

```julia
julia> flags(FontFlags)
(:bold, :italic, :large)

julia> flags(FontFlags(3))
(:bold, :italic)
```
"""
macro symbol_flagset(typespec::Union{Symbol,Expr}, flagspecs...)
    expand_flagset(typespec, flagspecs, true, __module__)
end

"""
    @flagset T [key_1 =] [bit_1 -->] flag_1 [key_2 =] [bit_2 -->] flag_2...

Create a [`FlagSet`](@ref) subtype with name `T` and the flags supplied (evaluated).
Each instance of the created type represents a subset of `{flag_1, flag_2, ...}`.
The associated `bit_1`,`bit_2`, etc, if provided, must be distinct positive powers of 2.

The following constructors are provided:
    - `T(itr)`: construct a set of the flags from the given iterable object. Each flag
    must be one of `flag_1`, `flag_2`, etc;
    - `T()`: construct an empty set;
    - `T(bitmask::Integer)`: construct a set of flags associated with each bit in `mask`;
    - `T(; kwargs...)`: each `key_i` can be given as a keyword with boolean value
    (`key_i = true` includes `flag_i`, `key_i = false` excludes it).

It should be safe to new constructor methods and/or overriding the methods provided
(except the default which is defined for the internal struct `FlagSets.BitMask`).

# Examples
```julia
julia> @flagset RoundingFlags begin
    down = RoundDown
    up = RoundUp
    near = 16 --> RoundNearest
    64 --> RoundNearestTiesUp  # No key
end

julia> RoundingFlags([RoundDown,RoundUp])
RoundingFlags with 2 elements:
  RoundingMode{:Down}()
  RoundingMode{:Up}()

julia> RoundingFlags(near = true, up = false, RoundNearestTiesUp = true)
RoundingFlags with 2 elements:
  RoundingMode{:Nearest}()
  RoundingMode{:NearestTiesUp}()

julia> RoundingFlags(1)
RoundingFlags with 1 element:
  RoundingMode{:Down}()

julia> for flag in RoundingFlags(3); println(flag) end
RoundingMode{:Down}()
RoundingMode{:Up}()

julia> RoundingFlags(3) | RoundingFlags(16)
RoundingFlags with 3 elements:
  RoundingMode{:Down}()
  RoundingMode{:Up}()
  RoundingMode{:Nearest}()

julia> eltype(RoundingFlags)
RoundingMode

julia> typemax(RoundingFlags)
RoundingFlags with 4 elements:
  RoundingMode{:Down}()
  RoundingMode{:Up}()
  RoundingMode{:Nearest}()
  RoundingMode{:NearestTiesUp}()

julia> typemin(RoundingFlags)
RoundingFlags([])
```

If some `flag_i` is a `Symbol` and `key_i` is not explicitly given, `key_i` is set to
`flag_i`. Therefore, both macros below are equivalent:
```
@flagset FontFlags :bold :italic 8 --> :large
@flagset FontFlags begin
    bold = :bold
    italic = :italic
    large = 8 --> :large
end
```

The following syntaxes can be used to provide `FlagType` and/or `BaseType`:

```julia
@flagset T {FlagType,BaseType} ... # Use supplied types
@flagset T {FlagType} ...          # Detect BaseType
@flagset T {FlagType,_} ...        # Detect BaseType too
@flagset T {_,BaseType} ...        # Detect FlagType
@flagset T {} ...                  # Braces are ignored
```

If `FlagType` is not provided or is `_`, it is inferred from the flag values
(like array constructors do).
Otherwise, `flag_1`, `flag_2`, etc, must be of type `FlagType`.

If `BaseType` is not provided or is `_`, some convenient integer type is inferred
(`UInt32` is used if possible, for maximum compatibility with C code).
Otherwise, it must be an `Integer` type big enough to represent each `bit_i`.

Each instance of `T` can be converted to `BaseType`, where each bit set in the
converted `BaseType` corresponds to a flag in the flag set. `read` and `write`
perform these conversions automatically. In case the flagset is created with a non-default
`BaseType`, `Integer(flagset)` will return the integer `flagset` with the type `BaseType`.

```julia
julia> convert(Int, RoundingFlags(up = true, near = true))
6

julia> convert(Integer, RoundingFlags(up = true, near = true))
0x00000006
```

To list flags of a flag set type or instance, use `flags`, e.g.

```julia
julia> flags(RoundingFlags)
(RoundingMode{:Down}(), RoundingMode{:Up}(), RoundingMode{:Nearest}())

julia> flags(RoundingFlags(3))
(RoundingMode{:Down}(), RoundingMode{:Up}())
```

!!! note

    For backward compatibility, the syntax `@flagset T::BaseType ...` is equivalent to
    `@symbol_flagset T::BaseType ...`, but the former syntax deprecated
    (use [`symbol_flagset`](@ref) instead).
"""
macro flagset(typespec::Union{Symbol,Expr}, flagspecs...)
    try
        return expand_flagset(typespec, flagspecs, false, __module__)
    catch ex
        throw(ex isa UndefVarError ? undef_var_error_hint(ex) : ex)
    end
end

## Error creation

function undef_var_error_hint(ex)
    ErrorException(
        "An $ex has been caught during macroexpansion of @flagset. If you are using " *
        "the old syntax (before version 0.3), use @symbol_flagset or the following syntax:" *
        "\n @flagset T [bit_1 -->] :flag_1 [bit_2 -->] :flag_2 ..."
    )
end

function deprecated_syntax_message(typename, BaseType)
    ArgumentError(
        "Deprecated syntax for macro @flagset. Use @symbol_flagset or use the syntax:" *
        "\n @flagset $typename {Symbol,$BaseType)} [bit_1 -->] :flag_1 [bit_2 -->] :flag_2 ..."
    )
end

function invalid_bit_error(typename, flagspec)
    ArgumentError("invalid bit for FlagSet $typename: $flagspec; should be an integer positive power of 2")
end

function invalid_flagspec_error(typename, flagspec)
    ArgumentError("invalid flag argument for FlagSet $typename: $flagspec")
end

function overflow_error(typename, flagspec)
    ArgumentError("overflow in bit of flag $flagspec for FlagSet $typename")
end

function not_unique_error(field::Symbol, typename, value)
    ArgumentError("$field for FlagSet $typename is not unique: $value")
end

function invalid_type_error(kind::Symbol, typename, type)
    ArgumentError("invalid $kind type for FlagSet $typename: $type")
end

function check_flag_type(typename, FlagType)
    (isnothing(FlagType) || FlagType isa Type) ||
        throw(invalid_type_error(:flag, typename, FlagType))
end

function check_base_type(typename, BaseType)
    (isnothing(BaseType) || BaseType isa Type && BaseType <: Integer) ||
        throw(invalid_type_error(:base, typename, BaseType))
end

## Helper functions

function _to_nothing(expr, __module__)
    expr == :_ ? nothing : Core.eval(__module__, expr)
end

function parse_flag_set_type(typespec, of_type, symflags::Bool, __module__)
    BaseType = nothing
    FlagType = symflags ? Symbol : nothing
    typename = nothing
    if !isnothing(of_type)
        if typespec isa Symbol && 0 ≤ length(of_type.args) ≤ 2
            typename = typespec
            length(of_type.args) ≥ 1 && (FlagType = _to_nothing(of_type.args[1], __module__))
            length(of_type.args) ≥ 2 && (BaseType = _to_nothing(of_type.args[2], __module__))
            # else typename = nothing # indicate error
        end
    elseif (
        isa(typespec, Expr) &&
        typespec.head === :(::) &&
        length(typespec.args) == 2 &&
        isa(typespec.args[1], Symbol)
    )
        symflags || @warn(deprecated_syntax_message(typespec.args[1], typespec.args[2]), maxlog = 1)
        FlagType = Symbol
        typename = typespec.args[1]
        BaseType = Core.eval(__module__, typespec.args[2])
    elseif isa(typespec, Symbol)
        typename = typespec
        # else typename = nothing
    end
    if isnothing(typename)
        throw(ArgumentError("invalid type expression for FlagSet: $typespec"))
    end
    check_flag_type(typename, FlagType)
    check_base_type(typename, BaseType)
    (typename, FlagType, BaseType, symflags)
end

function parse_flag_spec(typename, FlagType, flagspec, symflags::Bool, __module__)
    key, flag, bit = nothing, nothing, nothing
    if (
        isa(flagspec, Expr) &&
        (flagspec.head === :(=) || flagspec.head === :kw) &&
        length(flagspec.args) == 2
    )
        key = flagspec.args[1]
        flag_value = flagspec.args[2]
    elseif FlagType == Symbol && flagspec isa Symbol
        key = flag = flagspec
        return (key, flag, nothing)
    elseif symflags
        flag_value = nothing # indicates an error
    else
        flag_value = flagspec
    end
    if (isa(flag_value, Expr) && (flag_value.head === :-->) && length(flag_value.args) == 2)
        bit = Core.eval(__module__, flag_value.args[1]) # allow consts and exprs, e.g. uint128"1"
        flag_expr = flag_value.args[2]
        flag = Core.eval(__module__, flag_expr)
    else
        flag_expr = flag_value
        val = Core.eval(__module__, flag_value)
        flag, bit = FlagType == Symbol && !(val isa Symbol) ? (key, val) : (val, nothing)
    end
    if isnothing(key)
        key = flag isa Symbol && Base.isidentifier(flag) && flag != :missing ? flag : nothing
    end
    if isnothing(flag_value) || !(key isa Union{Symbol,Nothing})
        throw(invalid_flagspec_error(typename, flagspec))
    end
    if !isnothing(bit)
        (bit isa Real && bit >= 1 && ispow2(bit)) || throw(invalid_bit_error(typename, flagspec))
    end
    (key, flag, bit)
end

function parse_flag_specs(typename, FlagType, BaseType, flagspecs, symflags::Bool, __module__)
    if length(flagspecs) == 1 && flagspecs[1] isa Expr && flagspecs[1].head === :block
        flagspecs = flagspecs[1].args
    end
    result = Any[]
    bit = one(BaseType)
    mask = zero(BaseType)
    flags = Set([])
    keys = Set([])
    for flagspec in flagspecs
        flagspec isa LineNumberNode && continue
        (key, flag, next_bit) = parse_flag_spec(typename, FlagType, flagspec, symflags, __module__)
        bit::BaseType = something(next_bit, bit)
        bit <= zero(BaseType) && throw(overflow_error(typename, flagspec))
        flag::FlagType = flag
        (bit & mask) != 0 && throw(not_unique_error(:bit, typename, bit))
        flag in flags && throw(not_unique_error(:flag, typename, flag))
        key in keys && throw(not_unique_error(:key, typename, key))
        push!(flags, flag)
        isnothing(key) || push!(keys, key)
        mask |= bit
        index = trailing_zeros(bit) + 1
        push!(result, (key, flag, bit, index))
        bit <<= 1
    end
    (sort!(result, by = t -> t[4]), mask)
end

function infer_types(FlagType, BaseType, parsed, mask)
    if isnothing(FlagType)
        # Construct an array from flags and let Julia tell us what is the best flag type
        FlagType = eltype(collect(p[2] for p ∈ parsed))
    end
    if isnothing(BaseType)
        BaseType = if mask ≤ typemax(UInt32)
            UInt32
        elseif mask ≤ typemax(UInt64)
            UInt64
        elseif mask ≤ typemax(UInt128)
            UInt128
        else
            BigInt
        end
    end
    (FlagType, BaseType)
end

function expand_flagset(typespec, flagspecs, symflags::Bool, __module__)
    of_type = nothing
    if !isempty(flagspecs) && first(flagspecs) isa Expr && first(flagspecs).head == :braces
        of_type = first(flagspecs)
        flagspecs = flagspecs[2:end]
    end
    if isempty(flagspecs)
        throw(ArgumentError("no arguments given for FlagSet $typespec"))
    end
    (typename, FlagType, BaseType, symflags) =
        parse_flag_set_type(typespec, of_type, symflags, __module__)

    FT = something(FlagType, Any)
    BT = something(BaseType, BigInt)
    parsed, mask = parse_flag_specs(typename, FT, BT, flagspecs, symflags, __module__)
    (FlagType, BaseType) = infer_types(FlagType, BaseType, parsed, mask)

    len = last(parsed)[4] # last index
    key_vector = Vector{Union{Symbol,Nothing}}(undef, len)
    fill!(key_vector, nothing) # Avoid undefined reference errors
    flag_vector = Vector{Union{FlagType,Nothing}}(undef, len)
    fill!(flag_vector, nothing) # Avoid undefined reference errors
    indices = BaseType[]
    flag_bit_map = Base.ImmutableDict{FlagType,BaseType}()

    for (key, flag, bit, index) ∈ parsed
        key::Union{Symbol,Nothing} = key
        flag::FlagType = flag
        bit::BaseType = bit
        key_vector[index] = key
        flag_vector[index] = flag
        flag_bit_map = Base.ImmutableDict(flag_bit_map, flag => bit)
        push!(indices, index)
    end

    flagset_flags = Tuple(flag_vector)
    flags = Tuple(flag_vector[idx] for idx in indices)
    keys = Tuple(key_vector[idx] for idx in indices)
    keys_flags = ((key, flag) for (key, flag) ∈ zip(keys, flags) if !isnothing(key))

    blk = quote
        # flagset definition
        Base.@__doc__(struct $(esc(typename)) <: FlagSet{$FlagType,$BaseType}
                bitflags::$BaseType
                function $(esc(typename))(bitmask::FlagSets.BitMask)
                    bitmask = bitmask.bitmask
                    bitmask & $(mask) == bitmask || flagset_argument_error($(esc(typename)), bitmask)
                    return new(convert($BaseType, bitmask))
                end
            end)
        function $(esc(typename))(; $((Expr(:kw, :($key::Bool), false) for (key, _) ∈ keys_flags)...))
            bitmask::$BaseType = zero($BaseType)
            $((:($key && (bitmask |= $(flag_bit_map[flag]))) for (key, flag) ∈ keys_flags)...)
            $(esc(typename))(FlagSets.BitMask(bitmask))
        end

        FlagSets.flagset_flags(::Type{$(esc(typename))}) = $(esc(flagset_flags))
        FlagSets.flags(::Type{$(esc(typename))}) = $(esc(flags))
        FlagSets.flag_bit_map(::Type{$(esc(typename))}) = $(esc(flag_bit_map))
        FlagSets.flagkeys(::Type{$(esc(typename))}) =
            $(keys === flags ? :(FlagSets.flags($(esc(typename)))) : esc(keys))
        Base.typemin(x::Type{$(esc(typename))}) = $(esc(typename))(FlagSets.BitMask($(zero(BaseType))))
        Base.typemax(x::Type{$(esc(typename))}) = $(esc(typename))(FlagSets.BitMask($mask))
    end
    push!(blk.args, :nothing)
    blk.head = :toplevel
    return blk
end


text:
# The Bio.Phylo module
# ====================
#
# Types and methods for phylogenetic trees.
#
# Part of the Bio.Phylo module.
#
# This file is a part of BioJulia. License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

module Phylo

using Compat

using LightGraphs

using Bio.Tokenize

using Bio.Indexers

export

    # Types
    Phylogeny,

    # Phylogeny Methods

    ## Roots & Clades
    isrooted,
    isrerootable,
    root,
    isroot,
    clades,
    isclade,

    ## Leaves / Taxa
    leaves,
    isleaf,

    ## Structure manipulation
    hasparent,
    parent,
    haschildren,
    nchildren,
    haschild,
    children,
    isredundant,
    ispreterminal,
    issemipreterminal,

    ## Branches
    branchlength,
    branchlength!,
    parent_branch,
    child_branches,
    rem_branch!,
    add_branch!,

    ## Misc
    n_possible_rooted,
    n_possible_unrooted

include("phylogeny.jl")
include("node_basics.jl")
include("branch_basics.jl")
include("manipulation.jl")
include("dating.jl")

end # module Phylo


text:
for N in [Float64, Rational{Int}, Float32]
    # random hyperrectangle
    rand(Hyperrectangle)

    # constructor with mixed vectors
    Hyperrectangle(sparsevec([1], N[1], 1), N[1])
    Hyperrectangle(N[1], sparsevec([1], N[1], 1))

    # center/radius/high/low/generators
    c = N[0, 0]
    r = N[1, 1]
    h = Hyperrectangle(c, r)
    @test center(h) == c
    @test radius_hyperrectangle(h) == r
    @test high(h) == r
    @test low(h) == -r
    for i in 1:2
        @test radius_hyperrectangle(h, i) == r[i]
        @test high(h, i) == r[i]
        @test low(h, i) == -r[i]
    end
    gens = [Vector(g) for g in generators(h)]
    @test ispermutation(gens, [N[1, 0], N[0, 1]])
    @test genmat(h) ∈ [N[1 0; 0 1], N[0 1; 1 0]]
    @test ngens(h) == 2
    h_flat = Hyperrectangle(N[1, 2, 3, 4, 5], N[1, 0, 2, 0, 3])
    @test collect(generators(h_flat)) ==
        [SingleEntryVector(1, 5, N(1)), SingleEntryVector(3, 5, N(2)),
         SingleEntryVector(5, 5, N(3))]
    @test genmat(h_flat) == N[1 0 0; 0 0 0; 0 2 0; 0 0 0; 0 0 3]
    @test ngens(h_flat) == 3

    # 1D Hyperrectangle
    h = Hyperrectangle(N[0], N[1])
    # Test Dimension
    @test dim(h) == 1
    # Test Support Vector
    d = N[1]
    @test σ(d, h) == N[1]
    d = N[-1]
    @test σ(d, h) == N[-1]

    # 2D Hyperrectangle
    h = Hyperrectangle(N[0, 0], N[1, 1])
    # Test Dimension
    @test dim(h) == 2
    # Test Support Vector
    d = N[1, 1]
    @test σ(d, h) == N[1, 1]
    d = N[-1, 1]
    @test σ(d, h) == N[-1, 1]
    d = N[-1, -1]
    @test σ(d, h) == N[-1, -1]
    d = N[1, -1]
    @test σ(d, h) == N[1, -1]

    # 2D Hyperrectangle not 0-centered
    h = Hyperrectangle(N[1, 2], N[1, 1])
    # Test Dimension
    @test dim(h) == 2
    # Test Support Vector
    d = N[1, 1]
    @test σ(d, h) == N[2, 3]
    d = N[-1, 1]
    @test σ(d, h) == N[0, 3]
    d = N[-1, -1]
    @test σ(d, h) == N[0, 1]
    d = N[0, -1]
    @test σ(d, h) == N[2, 1]

    # 2D Hyperrectangle not same radius in each direction
    h = Hyperrectangle(N[0, 0], N[1, 2])
    # Test Dimension
    @test dim(h) == 2
    # Test Support Vector
    d = N[1, 1]
    @test σ(d, h) == N[1, 2]
    d = N[-1, 1]
    @test σ(d, h) == N[-1, 2]
    d = N[-1, -1]
    @test σ(d, h) == N[-1, -2]
    d = N[1, -1]
    @test σ(d, h) == N[1, -2]

    # 2D Hyperrectangle not centered, not same radius, for vertex representation,
    # radius, and diameter
    h = Hyperrectangle(N[3, 2], N[2, 1])
    vl = vertices_list(h)
    # Test Vertices
    @test ispermutation(vl, [N[1, 1], N[1, 3], N[5, 1], N[5, 3]])
    # norm
    @test norm(h) == norm(N[5, 3], Inf)
    # radius
    @test radius(h) == norm(N[2, 1], Inf)
    # diameter
    @test diameter(h) == norm(N[5, 3] - N[1, 1], Inf)

    # alternative constructor
    c = ones(N, 2)
    r = N[2, 3]
    l = N[-1, -2]
    h = N[3, 4]
    H1 = Hyperrectangle(c, r)
    H2 = Hyperrectangle(low=l, high=h)
    @test H1.center == H2.center
    @test H1.radius == H2.radius
    @test_throws AssertionError Hyperrectangle(low=h, high=l)

    # constructor without bounds check
    @test_throws AssertionError Hyperrectangle(low=N[1], high=N[0]) # default: true
    @test_throws AssertionError Hyperrectangle(low=N[1], high=N[0], check_bounds=true)
    @test Hyperrectangle(low=N[1], high=N[0], check_bounds=false) isa Hyperrectangle

    # Test low and high methods for a hyperrectangle
    H = Hyperrectangle(low=l, high=h)
    @test low(H) == l
    @test high(H) == h

    # boundedness
    @test isbounded(H)

    # isempty
    @test !isempty(H)

    # isuniversal
    answer, w = isuniversal(H, true)
    @test !isuniversal(H) && !answer && w ∉ H

    # membership
    H = Hyperrectangle(N[1, 1], N[2, 3])
    @test N[-1.1, 4.1] ∉ H
    @test N[-1, 4] ∈ H

    # "robust" membership of the support vector
    c = N[1.68, 0.73, 0.64]
    r = N[0.46, 0.24, 1.38]
    H = Hyperrectangle(c, r)
    @test σ(ones(N, 3), H) ∈ H

    # an_element function
    H = Hyperrectangle(N[1, 2], N[3, 4])
    @test an_element(H) ∈ H

    # isflat
    @test isflat(Hyperrectangle(N[1, 2], N[3, 0])) &&
          !isflat(Hyperrectangle(N[1, 2], N[3, 4]))

    # split
    @test ispermutation(split(H, [2, 2]),
        [Hyperrectangle(N[-0.5, 0], N[1.5, 2]),
         Hyperrectangle(N[2.5, 0], N[1.5, 2]),
         Hyperrectangle(N[-0.5, 4], N[1.5, 2]),
         Hyperrectangle(N[2.5, 4], N[1.5, 2])])

    # subset
    H1 = Hyperrectangle(N[1, 3], N[0.5, 0.5])
    H2 = Hyperrectangle(N[2, 2.5], N[0.5, 0.5])
    H3 = Hyperrectangle(N[2, 2], N[2, 3])
    B1 = BallInf(N[2, 2.5], N(0.5))
    B2 = BallInf(N[2, 2], N(1))
    @test H1 ⊈ H2 && H1 ⊆ H3 && H2 ⊆ H3
    subset, point = ⊆(H1, H2, true)
    @test !subset && point ∈ H1 && point ∉ H2
    subset, point = ⊆(H2, H1, true)
    @test !subset && point ∈ H2 && point ∉ H1
    subset, point = ⊆(H1, H3, true)
    @test subset
    @test H2 ⊆ B1 && B1 ⊆ H2
    @test B1 ⊆ B2 && B2 ⊈ B1

    # intersection & intersection emptiness
    H1 = Hyperrectangle(N[1, 1], N[2, 2])
    H2 = Hyperrectangle(N[3, 3], N[2, 2])
    B1 = BallInf(N[2, 4], N(0.5))
    B2 = BallInf(N[3, 3], N(0.5))
    for (X1, X2) in [(H1, H2), (H2, H1)]
        intersection_empty, point = is_intersection_empty(X1, X2, true)
        cap = intersection(X1, X2)
        @test cap isa Hyperrectangle{N} && center(cap) == N[2, 2] &&
              radius_hyperrectangle(cap) == N[1, 1]
        @test !is_intersection_empty(X1, X2) &&
              !intersection_empty && point ∈ X1 && point ∈ X2
    end
    cap = intersection(H1, B1)
    @test cap isa EmptySet{N}
    @test is_intersection_empty(H1, B1) && is_intersection_empty(H1, B1, true)[1]
    cap = intersection(H1, B2)
    @test cap == Hyperrectangle(N[2.75, 2.75], N[0.25, 0.25])
    intersection_empty, point = is_intersection_empty(H1, B2, true)
    @test !is_intersection_empty(H1, B2) && !intersection_empty && point ∈ H1 &&
          point ∈ B2

    # linear map (concrete)
    P = linear_map(N[1 0; 0 2], H1)
    @test P isa Zonotope
    P = linear_map(Diagonal(N[1, 2, 3, 4]), overapproximate(H1 * H1))
    @test P isa Zonotope

    # check that vertices_list is computed correctly if the hyperrectangle
    # is "degenerate"/flat, i.e., its radius contains zeros
    # these tests would crash if all 2^100 vertices were computed (see #92)
    H = Hyperrectangle(ones(N, 100), zeros(N, 100))
    vl = vertices_list(H)
    @test vl == [H.center]
    r = zeros(N, 100); r[1] = N(1)
    H = Hyperrectangle(fill(N(1), 100), r)
    vl = vertices_list(H)
    @test ispermutation(vl, [H.center + r, H.center - r])

    # transform hyperrectangle into a polygon
    H1pol = convert(HPolygon, H1)
    vlist = vertices_list(H1pol)
    @test ispermutation(vlist, [N[3, 3], N[3, -1], N[-1, -1], N[-1, 3]])

    # test that we can produce the list of constraints
    clist = constraints_list(H1)
    @test ispermutation(clist,
                        [HalfSpace(N[1, 0], N(3)), HalfSpace(N[0, 1], N(3)),
                         HalfSpace(N[-1, 0], N(1)), HalfSpace(N[0, -1], N(1))])

    # conversion to and from IntervalArithmetic's IntervalBox type
    B = IntervalBox(IntervalArithmetic.Interval(0, 1), IntervalArithmetic.Interval(0, 1))
    H = convert(Hyperrectangle, B)
    @test convert(IntervalBox, H) == B

    # conversion from other hyperrectangular sets
    H = Hyperrectangle(N[1], N[1])
    @test convert(Hyperrectangle, BallInf(N[1], N(1))) == H
    if VERSION >= v"1.1" || N == Float64
        # in pre-v1.1 versions, IntervalArithmetic always returned a Float64
        # vector in the `center(·)` implementation, leading to an error
        @test convert(Hyperrectangle, Interval(N(0), N(2))) == H
    end

    @test convert(Hyperrectangle, SymmetricIntervalHull(Singleton(N[1, 1]))) ==
          Hyperrectangle(N[0, 0], N[1, 1])

    # conversion to a zonotope
    H = Hyperrectangle(low=N[5//10, 6//10], high=N[124//10, 2355//10])
    Hz = convert(Zonotope, H)
    @test Hz == Zonotope(N[129//20, 2361//20], N[119//20 0//1; 0//1 2349//20])

    # conversion of a hyperrectangle with static array components to a zonotope
    # the specialized method for 2D static arrays is also tested
    H = Hyperrectangle(low=SA[N(5//10), N(6//10)], high=N[N(124//10), N(2355//10)])
    Hz = convert(Zonotope, H)
    Hz_sp = LazySets._convert_2D_static(Zonotope, H)
    Z = Zonotope(SA[N(129//20), N(2361//20)], SA[N(119//20) N(0//1); N(0//1) N(2349//20)])
    @test Hz == Z && Hz_sp == Z

    # rectification
    H = Hyperrectangle(N[-1, 2], N[4, 5])
    Hrect = rectify(H)
    @test Hrect.center == N[1.5, 3.5] &&  Hrect.radius == [1.5, 3.5]

    # Minkowski sum
    H1 = Hyperrectangle(N[0, 1], N[2, 3])
    H2 = Hyperrectangle(N[3, 2], N[1, 0])
    @test minkowski_sum(H1, H2) == Hyperrectangle(N[3, 3], N[3, 3])

    # set difference
    h = Hyperrectangle(low=N[0], high=N[1])
    q = Hyperrectangle(low=N[0], high=N[0.5])
    @test convert(Interval, difference(h, q).array[1]) == Interval(N(0.5), N(1))

    # another set difference test in higher-dimensions
    b = BallInf(N[0, 0, 0], N(1))
    @test isempty(difference(b, b))

    # volume of a hyperrectangular set
    b = BallInf(N[0, 0, 0], N(1))
    @test volume(b) ≈ N(8)

    # concete projection
    @test project(BallInf(zeros(N, 4), N(1)), [1, 2]) == Hyperrectangle(zeros(N, 2), ones(N, 2))
end


text:
@reexport module LeftRight

export LeftRightEnv

using ReinforcementLearningBase

using StatsBase

mutable struct LeftRightEnv <: AbstractEnv
    transitions::Array{Float64,3}
    current_state::Int
end

RLBase.get_observation_space(env::LeftRightEnv) = DiscreteSpace(2)
RLBase.get_action_space(env::LeftRightEnv) = DiscreteSpace(2)

function LeftRightEnv()
    t = zeros(2, 2, 2)
    t[1, :, :] = [0.9 0.1; 0.0 1.0]
    t[2, :, :] = [0.0 1.0; 0.0 1.0]
    LeftRightEnv(t, rand(1:2))
end

function (env::LeftRightEnv)(a::Int)
    env.current_state = sample(Weights(
        @view(env.transitions[env.current_state, a, :]),
        1.0,
    ))
    nothing
end

function RLBase.reset!(env::LeftRightEnv)
    env.current_state = 1
    nothing
end

RLBase.observe(env::LeftRightEnv) =
    (
        reward = Float64(env.current_state == 2),
        terminal = env.current_state == 2,
        state = env.current_state,
    )

end


text:
nindexes(N::Int) = N == 1 ? "1 index" : "$N indexes"

type FilledExtrapolation{T,N,ITP<:AbstractInterpolation,IT,GT,FT} <: AbstractExtrapolation{T,N,ITP,IT,GT}
    itp::ITP
    fillvalue::FT
end

function FilledExtrapolation{T,N,IT,GT}(itp::AbstractInterpolation{T,N,IT,GT}, fillvalue)
    FilledExtrapolation{T,N,typeof(itp),IT,GT,typeof(fillvalue)}(itp, fillvalue)
end

"""
`extrapolate(itp, fillvalue)` creates an extrapolation object that returns the `fillvalue` any time the indexes in `itp[x1,x2,...]` are out-of-bounds.
"""
extrapolate{T,N,IT,GT}(itp::AbstractInterpolation{T,N,IT,GT}, fillvalue) = FilledExtrapolation(itp, convert(eltype(itp), fillvalue))

@generated function getindex{T,N,ITP,IT,GT,FT}(fitp::FilledExtrapolation{T,N,ITP,IT,GT,FT}, args::Number...)
    n = length(args)
    n == N || return error("Must index $(N)-dimensional interpolation objects with $(nindexes(N))")
    Tret = FT<:Number ? getindex_return_type(ITP, args) : FT
    meta = Expr(:meta, :inline)
    quote
        $meta
        # Check to see if we're in the extrapolation region, i.e.,
        # out-of-bounds in an index
        @nexprs $N d->((args[d] < lbound(fitp,d) || args[d] > ubound(fitp, d)) && return convert($Tret, fitp.fillvalue))
        # In the interpolation region
        return getindex(fitp.itp,args...)
    end
end

getindex{T}(fitp::FilledExtrapolation{T,1}, x::Number, y::Int) = y == 1 ? fitp[x] : throw(BoundsError())

lbound(etp::FilledExtrapolation, d) = lbound(etp.itp, d)
ubound(etp::FilledExtrapolation, d) = ubound(etp.itp, d)


text:
module TestTools

using Test
using Mimi

import Mimi:
    getproperty, pretty_string

#utils: pretty_string
@test pretty_string("camelCaseBasic") == pretty_string(:camelCaseBasic) == "Camel Case Basic"
@test pretty_string("camelWithAOneLetterWord") == pretty_string(:camelWithAOneLetterWord) == "Camel With A One Letter Word"
@test pretty_string("snake_case_basic") == pretty_string(:snake_case_basic) == "Snake Case Basic"
@test pretty_string("_snake__case__weird_") == pretty_string(:_snake__case__weird_) == "Snake Case Weird"

#utils: interpolate
stepsize = 2
final = 10
ts = 10
@test Mimi.interpolate(collect(0:stepsize:final), ts) == collect(0:stepsize/ts:final)

end #module


text:
using MatrixNetworks
import MatrixNetworks: bipartite_matching,edge_list
using SparseArrays
using LinearAlgebra
using Statistics

include("../lowrank_spectral/decomposeX.jl")
include("../lowrank_spectral/newbound_methods.jl")
include("../lowrank_spectral/newbound_methods_relaxed.jl")
include("../lowrank_spectral/helper_functions.jl")
include("../lowrank_spectral/greedy_lowrank_EigenAlign.jl")
include("../lowrank_spectral/updated_MN_bipartite_matching.jl")
                #assuming these are being called from src/

struct LowRankEigenAlign_Return{T} <: returnType
  matchScore::Union{Int,Vector{Int}}
  motifCounts::NTuple{2,Union{Int,Vector{Int}}}
  matching::Union{Dict{Int,Int},Vector{Int}}
  embedding::Tuple{Matrix{T},Matrix{T}}
  profile::Union{Nothing,Float64}
end

function align_networks_eigenalign(A,B,iters,method,bmatch; default_params = true)
  D = 0
  s1,s2,s3 = find_parameters(A,B)
  if default_params == false
    s1 += 100
    s2 += 10
    s3 += 5
    @show s1
    @show s2
    @show s3
  end

  c1 = s1+s2-2s3; c2 = s3-s2; c3 = s2
  Uk,Vk,Wk,W1,W2 = decomposeX_balance_allfactors(A,B,iters+1,c1,c2,c3)
  Un,Vn = split_balanced_decomposition(Uk,Wk,Vk)
  timematching = 0

  nA = size(A,1)
  nB = size(B,1)

  if method == "full_multiplication"
    Xsample = 1e6*Uk*Wk*Vk'
    ma,mb = edge_list_local(bipartite_matching_local(sparse(Xsample)))
  elseif method == "lowrank_unbalanced_union"
    U = copy(Uk)
    V = Vk*Wk'
    P = newbound_rounding(U,V)
    Xsample = spones(P).*Xn
    ma,mb = edge_list_local(bipartite_matching_local(Xsample))
  elseif method == "lowrank_unbalanced_best"
    U = copy(Uk)
    V = Vk*Wk'
    ma,mb,D = newbound_rounding_max_approx(U,V)
  elseif method == "lowrank_balanced_union"
    U = Un
    V = Vn
    Xn = Un*Vn'
    P = newbound_rounding(U,V,bmatch)
    Xsample = spones(P).*Xn
    ma,mb = edge_list_local(bipartite_matching_local(Xsample))
    D = nnz(P)/prod(size(P))
  elseif method == "lowrank_balanced_union_lowrank_evaluation"
    U = Un
    V = Vn
    X = newbound_rounding_lowrank_evaluation_relaxed(U,V,bmatch).*(10^8) #bmatch
    D = nnz(X)/prod(size(X))
    ma,mb = edge_list_local(bipartite_matching_local(X))
  elseif method == "lowrank_balanced_best"
    U = Un
    V = Vn
    ma,mb,D = newbound_rounding_max_approx(U,V)
  elseif method == "lowrank_balanced_best_eff"
    U = Un
    V = Vn
    ma,mb,dval = newbound_rounding_max_approx_eff(U,V)
    X = U*V'
    X = X*100
    ei,ej = edge_list_local(bipartite_matching_local(sparse(X)))
    v0 = sum(sparse(ei,ej,1,size(X,1),size(X,2)).*X)
    M = sparse(ma,mb,1,size(X,1),size(X,2))
    v1 = sum(M.*X)
    row,col,M = greedy_match(X)
    v2 = sum(M.*X)
    D = [dval,v0/v1,v0/v2]
  elseif method == "lowrank_Wkdecomposed_union"
    U = Uk*W1
    V = Vk*W2'
    P = newbound_rounding(U,V)
    Xn = Un*Vn'
    Xsample = spones(P).*Xn
    ma,mb = edge_list_local(bipartite_matching_local(Xsample))
  elseif method == "lowrank_Wkdecomposed_best"
    U = Uk*W1
    V = Vk*W2'
    ma,mb,D = newbound_rounding_max_approx(U,V)
  elseif method == "lowrank_svd_union"

    U,S,V = svd(Wk)
    U1 = Uk*U*diagm(sqrt.(S)); V1 = Vk*V*diagm(sqrt.(S));
    X = newbound_rounding_lowrank_evaluation_relaxed(U1,V1,bmatch).*(10^8); #bmatch
    avgdeg = map(x->sum(X[x,:].!=0),1:size(X,1))
    avgdeg = mean(avgdeg)
    (ma,mb),timematching = @timed edge_list_local(bipartite_matching_local(X;normalize=false))
    D = avgdeg;#nnz(X)/prod(size(X))
  elseif method == "lowrank_lu_union"
    L,U = lu(Wk,NoPivot())
    U1 = Uk*L;
    V1 = Vk*U';
    X = newbound_rounding_lowrank_evaluation_relaxed(U1,V1,bmatch).*(10^8) #bmatch
    ma,mb = edge_list_local(bipartite_matching_local(X))

  elseif method == "greedy"
    U,S,V = svd(Wk)
    U1 = Uk*U*diagm(sqrt.(S)); V1 = Vk*V*diagm(sqrt.(S));
    ma,mb = greedy_lowrank(U1,V1)
  else
    error("method should be one of the following: (1)eigenalign,
          (2)lowrank_unbalanced_best, (3)lowrank_unbalanced_union,
          (4)lowrank_balanced_best, (5)lowrank_balanced_union,
          (6)lowrank_Wkdecomposed_best, (7)lowrank_Wkdecomposed_union")
  end

  return U1,V1,ma,mb,D,timematching
end

function find_parameters(A,B)
  nB = size(B,1)
  nA = size(A,1)
  nmatches = sum(A)*sum(B)
  nmismatches = sum(A)*(nB^2 - sum(B)) + sum(B)*(nA^2 - sum(A))
  mygamma = nmatches/nmismatches
  myalpha = (1/mygamma) + 1
  myeps = 0.001
  s1 = myalpha + myeps
  s2 = 1+myeps
  s3 = myeps
  return s1,s2,s3
end

function balance_decomposition(Uk,Wk,Vk)
  Du = diagm(vec(maximum(Uk,1)))
  Dv = diagm(vec(maximum(Vk,1)))
  Ukt = Uk*diagm(vec(1.0 ./maximum(Uk,1)))
  Vkt = Vk*diagm(vec(1.0 ./maximum(Vk,1)))
  rho = (maximum(Du)*maximum(Dv))
  Wkt2 = Du*Wk*Dv./rho
  rho *= maximum(Wkt2)
  Wkt2 = Wkt2/maximum(Wkt2)
  L,U = lu(Wkt2,NoPivot())
  Ud = diagm(sqrt.(abs.(diag(U)))) ##here
  L2 = L*Ud
  U2 = diagm(ones(size(U,1))./sqrt.(diag(U)))*U
  Un = Ukt*L2
  Vn = Vkt*U2'
  Xn = Ukt*L2*U2*Vkt'
  return Un,Vn,Xn
end

# can use sqrt here (like above)
function split_balanced_decomposition(Uk,Wk,Vk)
  L,U = lu(Wk,NoPivot())
  Ud = Diagonal(sqrt.(abs.(diag(U))))
  L2 = L*Ud
  U2 = Diagonal(ones(size(U,1))./sqrt.(diag(U)))*U
  Un = Uk*L2
  Vn = Vk*U2'
  return Un,Vn
end

function split_svd(Uk,Wk,Vk)
  U,S,V = svd(Wk)
  D = Diagonal(sqrt.(S))
  Unew = Uk*U*D
  Vnew = Vk*V*D
  return Unew,Vnew
end

max_stars_repo_name: chuckcol/LambdaTAME

max_issues_repo_name: chuckcol/LambdaTAME

max_forks_repo_name: chuckcol/LambdaTAME


text:
module ml2gen1
using StaticArrays
using MeshCore: nshapes
using MeshSteward: vtkwrite
using MeshSteward: L2blockx
using Test
function test()
    connectivity = L2blockx([0.0, 1.0, 3.0])
    @test (nshapes(connectivity.right), nshapes(connectivity.left)) == (3, 2)
    vtkwrite("ml2gen1a", connectivity)
    try rm("ml2gen1a.vtu"); catch end
    true
end
end
using .ml2gen1
ml2gen1.test()

module ml2gen1b
using StaticArrays
using MeshCore: nshapes
using MeshSteward: vtkwrite
using MeshSteward: L2blockx
using Test
function test()
    connectivity = L2blockx([0.0, 1.0, 3.0])
    @test (nshapes(connectivity.right), nshapes(connectivity.left)) == (3, 2)
    vtkwrite("ml2gen1b", connectivity)
    try rm("ml2gen1b.vtu"); catch end
    true
end
end
using .ml2gen1b
ml2gen1b.test()

module ml2gen2
using StaticArrays
using MeshCore: nshapes
using MeshSteward: vtkwrite
using MeshSteward: L2block
using Test
function test()
    connectivity = L2block(1.0, 7)
    # @show (nshapes(connectivity.right), nshapes(connectivity.left))
    @test (nshapes(connectivity.right), nshapes(connectivity.left)) == (8, 7)
    vtkwrite("ml2gen2", connectivity)
    try rm("ml2gen2.vtu"); catch end
    true
end
end
using .ml2gen2
ml2gen2.test()

module ml2gen3
using StaticArrays
using MeshCore: nshapes
using MeshSteward: vtkwrite
using MeshSteward: export_MESH
using MeshSteward: L2blockx, gradedspace
using Test
function test()
    connectivity = L2blockx(gradedspace(0.0, 5.0, 7, 2))
    # @show (nshapes(connectivity.right), nshapes(connectivity.left))
    @test (nshapes(connectivity.right), nshapes(connectivity.left)) == (7, 6)
    vtkwrite("ml2gen3", connectivity)
    try rm("ml2gen3.vtu"); catch end
    export_MESH("ml2gen3", connectivity)
    try rm("ml2gen3.mesh"); catch end
    try rm("ml2gen3-xyz.dat"); catch end
    try rm("ml2gen3-conn.dat"); catch end
    true
end
end
using .ml2gen3
ml2gen3.test()


text:
"Problem 11 in the Hock-Schittkowski suite"
function hs11()

  nlp = Model()

  @variable(nlp, x[i=1:2])
  setvalue(x[1], 4.9)
  setvalue(x[2], 0.1)

  @NLobjective(
    nlp,
    Min,
    (x[1] - 5)^2 + x[2]^2 - 25
  )

  @NLconstraint(
    nlp,
    -x[1]^2 + x[2] >= 0
  )

  return nlp
end


text:
mutable struct Xoshiro256p <: AbstractStreamableRNG
    Cg::Vector{UInt64}  # the current state of the RNG
    Bg::Vector{UInt64}  # the start point of the current substream
    Ig::Vector{UInt64}  # the start point of the current stream

    function Xoshiro256p(x::Vector{UInt})
        new(copy(x),copy(x),copy(x))
    end
    function Xoshiro256p(x::Vector{UInt}, y::Vector{UInt}, z::Vector{UInt})
        new(copy(x),copy(y),copy(z))
    end
end
rolt(x::UInt64, k::Int64) = (x<<k) | (x >> (64 - k))

function copy(Xos::Xoshiro256p)
    return Xoshiro256p(Xos.Cg, Xos.Bg, Xos.Ig)
end

function next(xos::Xoshiro256p)
    result_plus = xos.Cg[1] +xos.Cg[4]

    t = xos.Cg[2] << 17

    xos.Cg[3] = xor(xos.Cg[3], xos.Cg[1])
    xos.Cg[4] = xor(xos.Cg[4], xos.Cg[2])
    xos.Cg[2] = xor(xos.Cg[2], xos.Cg[3])
    xos.Cg[1] = xor(xos.Cg[1], xos.Cg[4])

    xos.Cg[3] = xor(xos.Cg[3], t)

    xos.Cg[4] = rolt(xos.Cg[4], 45)

    return result_plus
end

function short_jump(xos::Xoshiro256p)
    JUMP = [0x180ec6d33cfd0aba, 0xd5a61266f0c9392c, 0xa9582618e03fc9aa, 0x39abdc4529b1661c]
    xos.Cg[:] = xos.Bg[:]

    s = [UInt64(0) for i in 1:4]
    for i in 1:4
        for b in 0:63
            if (JUMP[i] & (UInt64(1) << b)) != 0
                #println("-")
                s[:] = xor.(s[:], xos.Cg[:])
            end
            next(xos)
        end
    end
    xos.Cg[:] = s[:]
    xos.Bg[:] = s[:]
end

function long_jump(xos::Xoshiro256p)
    LONG_JUMP = [0x76e15d3efefdcbbf, 0xc5004e441c522fb3, 0x77710069854ee241, 0x39109bb02acbe635]
    xos.Cg[:] = xos.Ig[:]
    s = [UInt64(0) for i in 1:4]
    for i in 1:4
        for b in 0:63
            if (LONG_JUMP[i] & (UInt64(1) << b)) != 0
                s[:] = xor.(s[:], xos.Cg[:])
            end
            next(xos)
        end
    end
    xos.Cg[:] = s[:]
    xos.Bg[:] = s[:]
    xos.Ig[:] = s[:]
end


text:
using AverageShiftedHistograms
using StatPlots
using Plots
using Distributions
plotly()
using JLD2

function Bias(A::Array{<:Real,2},B::Array{<:Real,2})
    return  Bias =    (A.-B)./std(A .-B)
end
Bias(A::SharedArray{<:Real,2},B::SharedArray{<:Real,2}) = BiasRel( convert(Array{Real},A),convert(Array{Real},B))

##

N = 10
GBA = 1
GW = 10
N_est = 900

Tvals= [500 1000 5000 10000] #   [100 250 500 1000]#
useStartVal = false
useStartVal? startValFlag = "T" : startValFlag = "F";# false #false #
folderName = "./data/estimatesTest/"
file_name = folderName * "EstTest_$(N)_$(GW)_$(GBA)_$(Tvals)_$(N_est)_$(useStartVal).jld"
estDataT = @load( file_name,simPar_1,simPar_2,simPar_3,estPar_1,estPar_2,estPar_3,est_times)
#

T  = Tvals[1]
ind = find(x->x==T,Tvals)

plot_tot =Vector{Plots.Plot{Plots.PlotlyBackend}}

rng_space = -3:.1:3
pW1 =0
pB1 =0
pA1 =0
for ind = 1:4
#ind = 1

errA = Bias(simPar_3[ind,:,:],estPar_3[ind,:,:])
errB = Bias(simPar_2[ind,:,:],estPar_2[ind,:,:])
errW =  (simPar_1[ind,:,:].- estPar_1[ind,:,:])./(repmat(std(estPar_1[ind,:,:],1),N_est,1))

err= errW
tmpAsh = ash(err;rng =rng_space )
pW1 = plot(tmpAsh,
    title ="  T= $(Tvals[ind]),
    W   mu =  $(round(mean(err[1,:]),2)) ")
plot!(Normal(0,1),linewidth = 5)

err= errB
tmpAsh = ash(err;rng = rng_space)
pB1 = plot(tmpAsh,
title =    " $(N_est) Est   B   mu =  $(round(mean(err[1,:]),2)) ")
plot!(Normal(0,1),linewidth = 5)

err= errA
tmpAsh = ash(err;rng = rng_space)
pA1 = plot(tmpAsh,
title = startValFlag * "  $(GW) / $(GBA) / $(N) A  mu =  $(round(mean(err[1,:]),2)) ")
plot!(Normal(0,1),linewidth = 5)

exp = "pl_T_$(ind)=  plot(pW1,pB1,pA1,layout=(3,1) )"
eval(parse(exp))
end

exp2 = " " ; for i = 1:ind exp2 *=  "pl_T_$(i) , " end

exp3 = "plot( " * exp2 * " layout = (1,ind),size=(1350,600),legend =:none)"
eval(parse(exp3))
#plot(pl_T_1 , pl_T_2,pl_T_3 , pl_T_4, layout = (1,4),size=(1350,600),legend =:none)

##

plot(errW[1,:])

##


text:
function get_constant_normalized_provider(obj::GravityFields, arg0::jint, arg1::jint)
    return jcall(obj, "getConstantNormalizedProvider", NormalizedSphericalHarmonicsProvider, (jint, jint), arg0, arg1)
end

function get_constant_unnormalized_provider(obj::GravityFields, arg0::jint, arg1::jint)
    return jcall(obj, "getConstantUnnormalizedProvider", UnnormalizedSphericalHarmonicsProvider, (jint, jint), arg0, arg1)
end

function get_normalized_provider(obj::GravityFields, arg0::jint, arg1::jint)
    return jcall(obj, "getNormalizedProvider", NormalizedSphericalHarmonicsProvider, (jint, jint), arg0, arg1)
end

function get_ocean_tides_waves(obj::GravityFields, arg0::jint, arg1::jint)
    return jcall(obj, "getOceanTidesWaves", List, (jint, jint), arg0, arg1)
end

function get_unnormalized_provider(obj::GravityFields, arg0::jint, arg1::jint)
    return jcall(obj, "getUnnormalizedProvider", UnnormalizedSphericalHarmonicsProvider, (jint, jint), arg0, arg1)
end


text:
using TestImages, Base64

test_name, test_image = first(
    first(splitext(img_name)) => testimage(img_name)
    for img_name
    in TestImages.remotefiles
    if endswith(img_name, ".png")
)

@testset "IOBuffer" begin
    @test test_image == let
        b = IOBuffer()
        PNGFiles.save(b, test_image)
        seek(b, 0)
        PNGFiles.load(b)
    end
end

@testset "Base64EncodePipe" begin
    @test test_image == let
        io = Base64EncodePipe(IOBuffer())
        PNGFiles.save(io, test_image)
        close(io)
        decoded = IOBuffer(base64decode(take!(io.io)))
        PNGFiles.load(decoded)
    end
end


text:
# Autogenerated wrapper script for NOVAS_jll for aarch64-linux-gnu
export libnovas

JLLWrappers.@generate_wrapper_header("NOVAS")
JLLWrappers.@declare_library_product(libnovas, "libnovas.so.3")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libnovas,
        "lib/libnovas.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# This file is a part of LegendDataTypes.jl, licensed under the MIT License (MIT).

using LegendDataTypes
using Test

using EncodedArrays

@testset "radware_sigcompress" begin
    function test_enc_dec(data::Vector{<:Real}, ref_codeunits::Vector{UInt8})
        codec = @inferred RadwareSigcompress(eltype(data))

        data_enc = EncodedArray{eltype(data)}(codec, size(data), ref_codeunits)

        @test @inferred(data |> codec) == data_enc
        @test @inferred(Array(data_enc)) == data
        @test typeof(Array(data_enc)) == typeof(data)
        @test @inferred(Array(data |> codec))  == data
    end

    test_enc_dec(
        [5, -7, -14, 17, -21, 0, -10, -2, -17, -14, 22, -5, -7, 7, 14, -2, 1, 0, -7, -21, -5, -15, -5, 11, 2, -24, 18, 2, -9, -3, 4, 7, -14, 12, 8, -23, 9, -8, 6, -2, 3, 9, -4, -10, -20, 5, 9, 12, 21, 3, 95, 189, 294, 405, 487, 593, 702, 781, 896, 1034, 986, 1004, 1011, 996, 990, 1005, 1003, 987, 1003, 994, 1008, 996, 989, 995, 997, 996, 985, 1010, 1006, 995, 988, 993, 1013, 1005, 985, 1025, 1008, 995, 1004, 1003],
        UInt8[0x00, 0x5a, 0x00, 0x32, 0x00, 0x06, 0xff, 0xe8, 0x75, 0x12, 0xa9, 0x0d, 0x83, 0x96, 0x1c, 0xab, 0x93, 0x45, 0xf9, 0x96, 0x65, 0x84, 0x43, 0x4c, 0x94, 0xe3, 0x68, 0x0a, 0x9a, 0x3d, 0x57, 0x1f, 0x2a, 0x48, 0x01, 0x85, 0x07, 0x96, 0x6e, 0x15, 0x0e, 0x11, 0xd8, 0x64, 0xb5, 0xb0, 0x00, 0x28, 0x00, 0x28, 0x00, 0x5f, 0xff, 0xd0, 0x8e, 0x99, 0x9f, 0x82, 0x9a, 0x9d, 0x7f, 0xa3, 0xba, 0x00, 0x42, 0x37, 0x21, 0x2a, 0x3f, 0x2e, 0x20, 0x40, 0x27, 0x3e, 0x24, 0x29, 0x36, 0x32, 0x2f, 0x25, 0x49, 0x2c, 0x25, 0x29, 0x35, 0x44, 0x28, 0x1c, 0x58, 0x1f, 0x23, 0x39, 0x2f, 0x00, 0x00, 0x00]
    )

    test_enc_dec(
        [107, 105, 113, 112, 105, 91, 119, 126, 110, 117, 105, 98, 129, 91, 112, 102, -33, 213, -54, 312, 107, 97, 107, 123, 114, 88, 130, 114, 103, 109, 116, 119, 98, 124, 120, 89, 121, 104, 118, 110, 115, 121, 108, 102, 92, 117, 121, 124, 133, 115, 207, 301, 406, 517, 599, 705, 814, 893, 1008, 1146, 1098, 1116, 1123, 1108, 1102, 1117, 1115, 1099, 1115, 1106, 1120, 1108, 1101, 1107, 1109, 1108, 1097, 1122, 1118, 1107, 1100, 1105, 1125, 1117, 1097, 1137, 1120, 1107, 1116, 1115],
        UInt8[0x00, 0x5a, 0x00, 0x35, 0x00, 0x09, 0xff, 0xca, 0x50, 0xa7, 0xd4, 0xea, 0x64, 0xfa, 0x45, 0x5a, 0xb4, 0x52, 0x2a, 0xd3, 0xe9, 0x85, 0xba, 0x45, 0x4c, 0x9c, 0x0a, 0xc2, 0xc0, 0x16, 0xe5, 0x0a, 0x5d, 0x42, 0xb1, 0x54, 0x23, 0x97, 0x0a, 0x84, 0xea, 0x8d, 0x54, 0xad, 0x4c, 0x2c, 0x95, 0xc8, 0xf5, 0x7a, 0x79, 0x58, 0xa4, 0x54, 0xab, 0xd4, 0x49, 0xc4, 0x92, 0xad, 0x5e, 0xb2, 0x5d, 0xaa, 0x60, 0xb6, 0x3e, 0x60, 0x00, 0x25, 0x00, 0x28, 0x02, 0x05, 0xff, 0xd0, 0x82, 0x9a, 0x9d, 0x7f, 0xa3, 0xba, 0x00, 0x42, 0x37, 0x21, 0x2a, 0x3f, 0x2e, 0x20, 0x40, 0x27, 0x3e, 0x24, 0x29, 0x36, 0x32, 0x2f, 0x25, 0x49, 0x2c, 0x25, 0x29, 0x35, 0x44, 0x28, 0x1c, 0x58, 0x1f, 0x23, 0x39, 0x2f]
    )
end


text:
import .Filters: AbstractCTFilter, CTFilter, RamLak

struct IsIRadonDiag end
struct IsIRadonSquare end

include("fbp_fft.jl")
include("fbp_square.jl")
include("fbpa.jl")
include("fbpa_square.jl")

@inline function _iradon(
    ::IsIRadonDiag,
    a::A,
    sinog::AbstractMatrix{T};
    rows::Optional{I} = nothing,
    cols::Optional{J} = nothing,
    α::Real = 360,
    α₀::Real = 0,
    kwargs...
) where {
    A <: AbstractIRadonAlgorithm,
    T <: Real,
    I <: Integer,
    J <: Integer,
}
    nd, nϕ = size(sinog)
    rows = maybe(cols, rows)
    rows = maybe(round(Int, nd / √2), rows)
    cols = maybe(rows, cols)
    ϕs = deg2rad(α₀)..(deg2rad(α₀) + deg2rad(α))
    sθ, cθ = sincos(atan(rows, cols))
    x₀, y₀ = cθ, sθ
    xs = linspace(eltype(sinog), -x₀..x₀, cols)
    ys = linspace(eltype(sinog), -y₀..y₀, rows)
    a(sinog, xs, ys; ϕs, kwargs...)
end

@inline function _iradon(
    ::IsIRadonDiag,
    a::A,
    sinog::AbstractMatrix{T},
    xsi::Interval,
    ysi::Interval,
    coo::Cartesian = Cartesian();
    rows::Optional{I} = nothing,
    cols::Optional{J} = nothing,
    kwargs...
) where {
    A <: AbstractIRadonAlgorithm,
    T <: Real,
    I <: Integer,
    J <: Integer,
}
    nd, nϕ = size(sinog)
    rows = maybe(cols, rows)
    rows = maybe(round(Int, nd / √2), rows)
    cols = maybe(rows, cols)
    xs = linspace(eltype(sinog), xsi, cols)
    ys = linspace(eltype(sinog), ysi, rows)
    a(sinog, xs, ys, coo; kwargs...)
end

@inline function _iradon(
    ::IsIRadonSquare,
    a::A,
    sinog::AbstractMatrix{T};
    rows::Optional{I} = nothing,
    cols::Optional{J} = nothing,
    α::Real = 360,
    α₀::Real = 0,
    kwargs...
) where {
    A <: AbstractIRadonAlgorithm,
    T <: Real,
    I <: Integer,
    J <: Integer,
}
    nd, nϕ = size(sinog)
    rows = maybe(cols, rows)
    rows = maybe(round(Int, nd), rows)
    cols = maybe(rows, cols)
    ϕs = deg2rad(α₀)..(deg2rad(α₀) + deg2rad(α))
    sθ, cθ = sincos(atan(rows, cols))
    x₀, y₀ = cθ, sθ
    xs = linspace(eltype(sinog), -x₀..x₀, cols)
    ys = linspace(eltype(sinog), -y₀..y₀, rows)
    a(sinog, xs, ys; ϕs, kwargs...)
end

@inline function _iradon(
    ::IsIRadonSquare,
    a::A,
    sinog::AbstractMatrix{T},
    xs::Interval,
    ys::Interval,
    coo::Cartesian = Cartesian();
    rows::Optional{I} = nothing,
    cols::Optional{J} = nothing,
    kwargs...
) where {
    A <: AbstractIRadonAlgorithm,
    T <: Real,
    I <: Integer,
    J <: Integer,
}
    nd, nϕ = size(sinog)
    rows = maybe(cols, rows)
    rows = maybe(round(Int, nd), rows)
    cols = maybe(rows, cols)
    xs = linspace(eltype(sinog), xsi, cols)
    ys = linspace(eltype(sinog), ysi, rows)
    a(sinog, xs, ys, coo; kwargs...)
end

@inline function iradon(sinog::AbstractMatrix; kwargs...)
    iradon(sinog, FBP(); kwargs...)
end

@inline function iradon(
    sinog::AbstractMatrix,
    xs::Union{AbstractVector,Interval},
    ys::Union{AbstractVector,Interval};
    kwargs...
)
    iradon(sinog, xs, ys, FBP(); kwargs...)
end

@inline function iradon(
    sinog::AbstractMatrix,
    geometry::AbstractGeometry;
    kwargs...
)
    iradon(sinog, geometry, FBP(); kwargs...)
end

@inline function iradon(g::AbstractGeometry; kwargs...)
    iradon(g, FBP(); kwargs...)
end


text:
function gaussian_data(d, n, f; σ=1)
    X = rand(d, n)
    y = mapslices(f, X; dims=1) + rand(Normal(0, σ), 1, n)
    return X, y
end


text:
# Bagging Algorithm
function bagged_tree(dtmtable::DTMTable, sett::ModelSettings)
    trnidx = dtmtable.trnidx
    validx = dtmtable.validx
    features = dtmtable.features
    weight = dtmtable.weight
    actualNumerator = dtmtable.numerator
    denominator = dtmtable.denominator
    mappings = dtmtable.mappings
    candMatWOMaxValues = dtmtable.candMatWOMaxValues
    obstrn, obsval, trn_meanobservedvalue, val_meanobservedvalue, trn_numtot, val_numtot, trn_denomtot, val_denomtot, empty_rows_after_iteration_stats, showTimeUsedByEachIteration, chosen_apply_tree_fn, startAtMean, adaptiveLearningRate, moderationvector, iterations, nameOflistofpredictorsSheet, nameOfpredictorsSheet, nameOfModelStatisticsSheet, nameOfScoresSheet, nameOfOverviewSheet, nameOfSettingsSheet, nameOfValidationSheet, xlData, showProgressBar_time, prroduceEstAndLeafMatrices = initSettingsWhichAreTheSameForBoostingAndBagging(trnidx, validx, actualNumerator, denominator, sett)
    obs = obstrn + obsval
    current_error = 0.0
    moderationfactor = moderationvector[1]

    if (size(moderationvector, 1) != iterations && size(moderationvector, 1) != 1)
	error("DTM: This should not have happend: size(mdf_vector) is neither 1 nor #iterations")
    end
    if size(moderationvector, 1) == 1
	actual_moderationvector = Array{Float64}(undef, iterations)
	fill!(actual_moderationvector, deepcopy(moderationvector[1]))
    else
	actual_moderationvector = deepcopy(moderationvector)
    end
    if startAtMean
	estimatedRatio = ones(length(weight)) .* trn_meanobservedvalue
    else
	error("DTM: This is currently not working: please set startAtMean=true)")
    end

# obs,obsVAL,trn_meanobservedvalue,empty_rows_after_iteration_stats,showTimeUsedByEachIteration,chosen_apply_tree_fn,startAtMean,adaptiveLearningRate,moderationvector,iterations,nameOflistofpredictorsSheet,nameOfpredictorsSheet,nameOfModelStatisticsSheet,nameOfScoresSheet,nameOfOverviewSheet,nameOfSettingsSheet,nameOfValidationSheet,xlData,showProgressBar_time=initSettingsWhichAreTheSameForBoostingAndBagging(actualNumerator,actualNumeratorVAL,denominator,sett)
    iterationsPerCore = calcIterationsPerCore(sett.iterations, Distributed.Distributed.nprocs())
    sampleSizeForEachTreeCanBeNEG = convert(Int, round(sett.subsampling_prop * length(trnidx)))
    if abs(sampleSizeForEachTreeCanBeNEG) < 1
	sampleSizeForEachTreeCanBeNEG = convert(Int, sign(sett.subsampling_prop))
    end

    vecTreesWErrs = Distributed.@distributed (vcat) for ii = 1:Distributed.Distributed.nprocs()
	# create_bagged_trees(iterationsPerCore[ii],sampleSizeForEachTreeCanBeNEG,trn_meanobservedvalue,candMatWOMaxValues,mappings,sett,actualNumerator,denominator,weight,numfeatures,charfeatures)
	# sample_data_and_build_tree!(trnidx,validx,indicatedRelativityForApplyTree_reused,candMatWOMaxValues,mappings,deepcopy(sett),actualNumerator,estimatedNumerator,weight,features,sampleSizeCanBeNEGATIVE,abssampleSize,sampleVector,T_Uint8_or_UInt16)
	create_bagged_trees(iterationsPerCore[ii], trnidx, sampleSizeForEachTreeCanBeNEG, trn_meanobservedvalue, candMatWOMaxValues, mappings, sett, actualNumerator, denominator, weight, features)
    end

	# Derive cumulative estimates for each iteration for trn and val data
		# treesOnly=[x.tree for x in vecTreesWErrs]
		# warn("think about the choice of goodness of fit criterium here!")
		weightPerTree = defineWeightPerTree(vecTreesWErrs, sett.baggingWeightTreesError)
		# this header needs to be in line with the output of the function createTrnValStatsForThisIteration
			statsPerIteration = global_statsperiter_header
			scoreBandLabels = createScorebandsUTF8List(sett.scorebandsstartingpoints, sett.nScores, addtotal=true)
		# this row is only here such that the variables are initialized outside the for loop
		local	obsPerScore, vectorWeightPerScore, numPerScore, denomPerScore, scores, scoresVAL, cumulativeStatsPerScoreBand, relativities, relativitiesVAL, maxRawRelativityPerScoreSorted, rawObservedRatioPerScore, MAPPINGSmoothedEstimatePerScore

		estimatedRatio = ones(obs) .* trn_meanobservedvalue
		currentEstimateOfIndividualTree = zeros(obs)
		leafNumbersThisTree = zeros(Int, obs)
		currentRelativity = ones(obs)
		sortvec_reused_trn_only = zeros(Int, length(trnidx))
		scores = zeros(Int, obs)
		estFromScores = zeros(obs)

		leaves_of_tree = Array{Array{Leaf{T},1}}(undef, iterations)

		if prroduceEstAndLeafMatrices
			est_matrix = Array{Float64}(undef, obs, iterations + 1)
			est_matrixFromScores = copy(est_matrix)
			MatrixOfLeafNumbers = Array{Int}(undef, obs, iterations + 1)
			MatrixOfLeafNumbers[:,1] = 0
			est_matrix[:,1] = copy(transpose(estimatedRatio))
			est_matrixFromScores[:,1] = copy(transpose(estimatedRatio))
		else
			est_matrix = Array{Float64}(undef, 0, 0)
			est_matrixFromScores = copy(est_matrix)
			MatrixOfLeafNumbers = Array{Int}(undef, 0, 0)
		end
    vectorOfLeafArrays = Array{Array{Leaf{T},1}}(undef, iterations + 1)
		vectorOfLeafArrays[1] = Array{Leaf{T}}(undef, 0) # the first entry is not defined
		p = ProgressMeter.Progress(iterations, 5, "Progress of Bagging Model:") # minimum update interval: 5 second
		cumulativeWeight = 0.0
		local fristRowOfThisTable
	for iter = 1:iterations
	error("BK: TODO, need to port the functionality regarding sett.fitForStatsAndCharts to bagging")
		# fit current tree
			vectorOfLeafArrays[iter + 1] = create_leaves_array(vecTreesWErrs[iter].tree)
			# note the individual trees were constructed on a part of the training data, here we apply them to the full training data set
			apply_tree_by_leaf!(currentEstimateOfIndividualTree, leafNumbersThisTree, trnidx, vecTreesWErrs[iter].tree.rootnode, features)
			apply_tree_by_leaf!(currentEstimateOfIndividualTree, leafNumbersThisTree, validx, vecTreesWErrs[iter].tree.rootnode, features)
			# currentEstimateOfIndividualTree[:],leafNumbersThisTree[:]=Core.eval(parse(sett.chosen_apply_tree_fn))(vectorOfRulePathsToLeavesArrays[iter+1],vecTreesWErrs[iter].tree.rootnode,numfeatures,charfeatures)
			# vectorOfLeafNumbersTrn[:,iter+1]=leafNumbersThisTree
		# update weights
			wi = weightPerTree[iter]
			cumulativeWeight += wi
		# derive current cumulative estimate
			estimatedRatio .= ((estimatedRatio .* cumulativeWeight) .+ (currentEstimateOfIndividualTree .* wi)) ./ (cumulativeWeight + wi)
			currentRelativity .= estimatedRatio ./ trn_meanobservedvalue

			if prroduceEstAndLeafMatrices
				est_matrix[:,iter + 1] = estimatedRatio
				MatrixOfLeafNumbers[:,iter + 1] = leafNumbersThisTree # leaf_numbers(vectorOfLeafArrays[1+iter],obs) #TODO / TBD adjust this for subsampling. this may not work properly, as each tree will use a different amount of data AND DIFFERENT OBSERVATIONS! thus we cannot rely on the *.idx field of the leaves!
			end
		# Derive Scores for this iteration
		maxRawRelativityPerScoreSorted, MAPPINGSmoothedEstimatePerScore, vectorWeightPerScore, obsPerScore, rawObservedRatioPerScore, numPerScore, denomPerScore, nscoresPotentiallyReduced				= constructANDderiveScores!(trnidx, validx, sortvec_reused_trn_only, estimatedRatio, currentRelativity, actualNumerator, denominator, weight, trn_meanobservedvalue, iter, sett)
		update_and_derive_scores!(scores, maxRawRelativityPerScoreSorted, currentRelativity)
	# Define estimated Numerator
		sett.smoothEstimates == "0" ? referenceForNumEstimates = rawObservedRatioPerScore : referenceForNumEstimates = MAPPINGSmoothedEstimatePerScore
		update_mapped_estFromScores!(estFromScores, referenceForNumEstimates, scores)
		# estFromScores=map(x->referenceForNumEstimates[x],scores)
		if prroduceEstAndLeafMatrices
			write_column!(est_matrixFromScores, iter + 1, estFromScores)
		end
		#	maxRawRelativityPerScoreSorted,MAPPINGSmoothedEstimatePerScore,vectorWeightPerScore,obsPerScore,rawObservedRatioPerScore,numPerScore,denomPerScore,scores,uniqueRelativitiesSorted,nscoresPotentiallyReduced=constructANDderiveScores(estimatedRatio,currentRelativity,sett.nScores,actualNumerator,denominator,weight,trn_meanobservedvalue,iter,sett.smoothEstimates,sett.print_details)
		# Apply this iteration to the validation data set
			# currentEstimateOfIndividualTreeVAL[:],leafNumbersThisTreeVAL[:]=Core.eval(parse(sett.chosen_apply_tree_fn))(vectorOfLeafArrays[iter+1],vecTreesWErrs.tree[iter.rootnode,numfeaturesVAL,charfeaturesVAL)
		# derive current cumulative estimate
			# estimatedRatioVAL[:]=(estimatedRatioVAL.*cumulativeWeight+currentEstimateOfIndividualTreeVAL.*wi)./(cumulativeWeight+wi)
			# currentRelativityVAL[:]=estimatedRatioVAL./trn_meanobservedvalue
			# est_matrixVAL[:,iter+1]=estimatedRatioVAL
			# scoresVAL=derive_scores(maxRawRelativityPerScoreSorted,currentRelativityVAL)
		# Define estimated Numerator
		# Derive Cumulative Statistics
			 # todo tbd there is potential for optimization here
			statsThisIteration, singleRowWithKeyMetrics, columnOfRelativityTrn = createTrnValStatsForThisIteration(scoreBandLabels, iter, sett.scorebandsstartingpoints, view(actualNumerator, trnidx), view(denominator, trnidx), view(weight, trnidx), view(estFromScores, trnidx), view(scores, trnidx), view(actualNumerator, validx), view(denominator, validx), view(weight, validx), view(estFromScores, validx), view(scores, validx), sett)
			statsPerIteration = vcat(statsPerIteration, singleRowWithKeyMetrics)
			if iter == 1
				fristRowOfThisTable = sett.iterations + 2 + empty_rows_after_iteration_stats
				cumulativeStatsPerScoreBand = copy(statsThisIteration)
			else
				fristRowOfThisTable += size(statsThisIteration, 1)
				cumulativeStatsPerScoreBand = vcat(cumulativeStatsPerScoreBand, statsThisIteration)
			end
			# Add the respective chart
				loc = string(excelLetter(size(cumulativeStatsPerScoreBand, 2) + 2), fristRowOfThisTable)
				thischart = defineRelativityChart(nameOfModelStatisticsSheet, nameOfModelStatisticsSheet, loc, size(statsThisIteration, 1) - 2, columnOfRelativityTrn, fristRowOfThisTable, xscale=2.0)
				push!(xlData.charts, deepcopy(thischart))

			# sett.smoothEstimates=="0" ? referenceForNumEstimates=rawObservedRatioPerScore : referenceForNumEstimates=MAPPINGSmoothedEstimatePerScore
			# estimatedNumeratorDerivedFromScores=[referenceForNumEstimates[x] for x in scores]
			# estimatedNumeratorDerivedFromScoresVAL=[referenceForNumEstimates[x] for x in scoresVAL]
			# est_matrixFromScores[:,iter+1]=estimatedNumeratorDerivedFromScores[:]
			# est_matrixFromScoresVAL[:,iter+1]=estimatedNumeratorDerivedFromScoresVAL[:]
			if showProgressBar_time # TimeUsedByEachIteration
				ProgressMeter.next!(p)
			end
	end
	# todo/tbd: a large part of this code is redundant to boosting, this should be calucalted in a separate function somewhere..... (such that we do not need to maintain two functions with the same code!)

	estimateUnsmoothed = map_these_values(rawObservedRatioPerScore, scores)
	estimateSmoothed = map_these_values(MAPPINGSmoothedEstimatePerScore, scores)
	estimateFromRelativities = currentRelativity .* trn_meanobservedvalue
	finalEstimateForCharts = sett.smoothEstimates == "0" ? estimateSmoothed : estimateUnsmoothed
	# add empty lines
	statsPerIteration = vcat(statsPerIteration, repeat([""], empty_rows_after_iteration_stats, size(statsPerIteration, 2)))
	# Add Charts which show stats for all iterations
		add_iteration_charts!(xlData, sett, 1)
		size(statsPerIteration, 2) > size(cumulativeStatsPerScoreBand, 2) ? cumulativeStatsPerScoreBand = hcat(cumulativeStatsPerScoreBand, repeat([""], size(cumulativeStatsPerScoreBand, 1), size(statsPerIteration, 2) - size(cumulativeStatsPerScoreBand, 2))) : statsPerIteration = hcat(statsPerIteration, repeat([""], size(statsPerIteration, 1), -size(statsPerIteration, 2) + size(cumulativeStatsPerScoreBand, 2)))
		stats = vcat(statsPerIteration, cumulativeStatsPerScoreBand)
# Create Score List
	scoreheader = ["Score" "Row Count" "Weight" "Numerator" "Denominator" "Maximal Raw Relativity" "Observed Ratio" "Smoothed Fit"]
	scoreMatrix = vcat(scoreheader, hcat(collect(1:size(maxRawRelativityPerScoreSorted, 1)), obsPerScore, vectorWeightPerScore, numPerScore, denomPerScore, maxRawRelativityPerScoreSorted, rawObservedRatioPerScore, MAPPINGSmoothedEstimatePerScore))
	scoreMatrix = hcat(scoreMatrix, repeat([""], size(scoreMatrix, 1), size(stats, 2) - size(scoreMatrix, 2)))
# Add Score Chart
	thischart = defineScoreChart(nameOfScoresSheet, nameOfScoresSheet, "B6", size(scoreMatrix, 1) - 1, 1, 3, 7, 8)
	push!(xlData.charts, deepcopy(thischart))
# Create univariate graphs
		predictorsData, predictorCharts = createPredictorData(trnidx, nameOfpredictorsSheet, mappings, candMatWOMaxValues, sett, scores, actualNumerator, denominator, weight, finalEstimateForCharts, features)
	# push predictorCharts into Excelobject
		push!(xlData.charts, predictorCharts...)
		predictorsSheet = ExcelSheet(nameOfpredictorsSheet, DataFrame(predictorsData,:auto))
	# Populate ExcelSheet object with data
		statsSheet = ExcelSheet(nameOfModelStatisticsSheet, DataFrame(stats,:auto))
		overviewSheet = ExcelSheet(nameOfOverviewSheet, DataFrame())
		scoreMatrixDF = DataFrame( scoreMatrix,:auto)
		scoresSheet = ExcelSheet(nameOfScoresSheet, scoreMatrixDF)
		modelsettingsSheet = ExcelSheet(nameOfSettingsSheet, DataFrame( writeAllFieldsToArray(sett),:auto))
# Create two way validation charts
	twoWayValidation, validationCharts = createTwoWayValidationCharts(trnidx, validx, nameOfValidationSheet, scoreBandLabels, mappings, candMatWOMaxValues, sett, scores, actualNumerator, denominator, weight, finalEstimateForCharts, features)
	validationSheet = ExcelSheet(nameOfValidationSheet, DataFrame(twoWayValidation,:auto))
	push!(xlData.charts, validationCharts...)
# Define Excel
	xlData.sheets = [modelsettingsSheet,overviewSheet,statsSheet,scoresSheet,predictorsSheet,validationSheet]

	z = deepcopy(stats[1:sett.iterations + 1,:])
	modelstats = DataFrame(convert(Array{Float64,2}, z[2:size(z, 1),:]), :auto)
	DataFrames.names!(modelstats, Symbol[Symbol(x) for x in view(z, 1, :)])
# resulting Bagged Tree, construct Bagged Tree
	intVarsUsed = [x.tree.intVarsUsed for x in vecTreesWErrs]
	inds_considered = [x.tree.inds_considered for x in vecTreesWErrs]
	nodesOnly = [x.tree.rootnode for x in vecTreesWErrs]
	errs = [x.err for x in vecTreesWErrs]

	# create trnidx
	@assert issorted(trnidx)
	trnidx_one_zero_full_length = map(x->UInt8(length(searchsorted(trnidx, x))), 1:length(scores))
	fp = get_feature_pools(features)
	resultingBaggedTree = BaggedTree(nodesOnly, weightPerTree, errs, sett, intVarsUsed, candMatWOMaxValues, mappings, inds_considered, scores, currentRelativity, maxRawRelativityPerScoreSorted, trn_meanobservedvalue, startAtMean, MAPPINGSmoothedEstimatePerScore, rawObservedRatioPerScorexx, est_matrix, modelstats, trnidx_one_zero_full_length, fp)

# resulting Bagged Tree
    return xlData, estimatedRatio, MatrixOfLeafNumbers, vectorOfLeafArrays, rawObservedRatioPerScore, est_matrixFromScores, stats, estimateUnsmoothed, estimateSmoothed, estimateFromRelativities, resultingBaggedTree
end

function create_bagged_trees(itr::Int, trnidx_which_should_only_be_used_once_here::Vector{Int}, sampleSizeCanBeNEGATIVE::Int, trn_meanobservedvalue::Float64, candMatWOMaxValues::Array{Array{Float64,1},1}, mappings::Array{Array{String,1},1}, sett::ModelSettings, numerator::Array{Float64,1}, denominator::Array{Float64,1}, weight::Array{Float64,1}, features::DataFrame)
# function create_bagged_trees(iterationsPerCore[ii],sampleSizeForEachTreeCanBeNEG,trn_meanobservedvalue,candMatWOMaxValues,mappings,sett,actualNumerator,denominator,weight,features)
	res = Array{TreeWithErrorStats}(itr)
	resFitted = Array{Array{Float64,1}}(undef, 1)

	abssampleSize, sampleVector = initBootstrapSample(sampleSizeCanBeNEGATIVE)
	# abssampleSize,sampleVector,num,denom,w,numf,charf,ooBagnum,ooBagdenom,ooBagw,ooBagnumf,ooBagcharf,ooBagsize=initBoostrapSample(sampleSizeCanBeNEGATIVE,numerator,denominator,weight,charfeatures,numfeatures)

	fitted_values_all_data_this_vector_is_modified_by_build_tree = zeros(length(numerator))
	reused_fitted_leafnr_vector = zeros(Int, length(numerator))
	estimatedNumFromRelativities = zeros(length(numerator))
	sortvec_reused_trn_only = zeros(Int, length(sampleVector))
	# modelling loop
	for i = 1:itr
		# sample data
		unusedSamplePart = sampleData!(trnidx_which_should_only_be_used_once_here, sampleSizeCanBeNEGATIVE, sampleVector)
	# build tree	on subsample of data
			thistree = build_tree!(sampleVector, unusedSamplePart, candMatWOMaxValues, mappings, deepcopy(sett), numerator, denominator, weight, features, fitted_values_all_data_this_vector_is_modified_by_build_tree, T_Uint8_or_UInt16)
	# apply tree to ooBagSample
			leaves_of_tree = create_leaves_array(thistree.rootnode)
			apply_tree_by_leaf!(fitted_values_all_data_this_vector_is_modified_by_build_tree, reused_fitted_leafnr_vector, sampleVector, thistree.rootnode, features)
			if length(unusedSamplePart) > 0
				apply_tree_by_leaf!(fitted_values_all_data_this_vector_is_modified_by_build_tree, reused_fitted_leafnr_vector, unusedSamplePart, thistree.rootnode, features)
			end
		# Determine Goodness of Fit
			# I think it is best to consider the goodness of fit per leaf here
			# out of bag performance should be measured here, however this is not coded yet (see also bagging!)
				# ooBagEstNumerator=ooBagEstimates.*ooBagdenom
				# fittedPerLeaf=Float64[x.fitted for x in leaves_of_tree]
				# ooBagcntperLeaf,ooBagsumnumeratorEST,ooBagsumnumerator,ooBagsumdenominator,ooBagsumweight=aggregateByLeafNr(leafNrooBag,ooBagEstNumerator,ooBagnum,ooBagdenom,ooBagw)
				relativities = fitted_values_all_data_this_vector_is_modified_by_build_tree ./ trn_meanobservedvalue
				estimatedRatio = fitted_values_all_data_this_vector_is_modified_by_build_tree
				resFitted[1] = fitted_values_all_data_this_vector_is_modified_by_build_tree
			# derive scores
			maxRawRelativityPerScoreSorted, MAPPINGSmoothedEstimatePerScore, vectorWeightPerScore, obsPerScore, rawObservedRatioPerScore, numPerScore, denomPerScore, nscoresPotentiallyReduced				= constructANDderiveScores!(sampleVector,
			unusedSamplePart,sortvec_reused_trn_only,estimatedRatio,relativities,numerator,denominator,weight,trn_meanobservedvalue,itr,sett)
			# old call constructScores(estimatedRatio,relativities,uniqueRelativitiesSorted,num,denom,w,trn_meanobservedvalue,ns,sett.nScores,sett.print_details)
			# uniqueRelativitiesSorted=unique(relativities)
			# sort!(uniqueRelativitiesSorted,alg=QuickSort)
			# ns=min(sett.nScores,size(uniqueRelativitiesSorted,1))
			# ns<nScores&&warn("There are only $(ns) distinct relativities. Number of scores will be limited to $(ns) instead of $(nScores).")
			# aggregatedModelledRatioPerScore,maxRawRelativityPerScoreSorted,MAPPINGSmoothedEstimatePerScore,vectorWeightPerScore,obsPerScore,rawObservedRatioPerScore,numPerScore,denomPerScore=	constructScores(estimatedRatio,relativities,uniqueRelativitiesSorted,num,denom,w,trn_meanobservedvalue,ns,sett.nScores,sett.print_details)

		# apply tree to ooBagSample
			if size(unusedSamplePart, 1) > 0
				# subsampling is performed
				# warn("todo: PICK a good function to determine the goodness of FIT / determine out of bag performance here! which will determine the weight")
				# @show sampleSizeCanBeNEGATIVE
				# if sampleSizeCanBeNEGATIVE<0
				#	resize!(leafNrooBag,size(unusedSamplePart,1))
				#	resize!(ooBagEstimates,size(unusedSamplePart,1))
				# end
				# leaves_of_tree=create_leaves_array(thistree.rootnode)
				# if sett.prroduceEstAndLeafMatrices
				#	error("this has not yet been updated for bagging")
					# todo/tbd this could possibly be optimized. We are calling apply_tree_by_leaf! twice (if prroduceEstAndLeafMatrices==true), once for val and once for trn
					# apply_tree_by_leaf!(fitted_values_all_data_this_vector_is_modified_by_build_tree,reused_fitted_leafnr_vector,sampleVector,thistree.rootnode,features)
					# MatrixOfLeafNumbers[:,iter+1]=copy(reused_fitted_leafnr_vector) #leaf_numbers(vectorOfLeafArrays[1+iter],obs) #TODO / TBD adjust this for subsampling. this may not work properly, as each tree will use a different amount of data AND DIFFERENT OBSERVATIONS! thus we cannot rely on the *.idx field of the leaves!
				# end
                # Determine Goodness of Fit
					# I think it is best to consider the goodness of fit per leaf here
					# ooBagEstNumerator=ooBagEstimates.*ooBagdenom

				# relativitiesVAL=ooBagEstimates./trn_meanobservedvalue
				# scoresVAL=derive_scores(maxRawRelativityPerScoreSorted,relativitiesVAL)
				# todo tbd decide which estimate we should use here.....!
				estimatedNumFromRelativities .= indicatedRelativityForApplyTree_reused .* denominator
				# sett.smoothEstimates=="0" ? referenceForNumEstimates=rawObservedRatioPerScore : referenceForNumEstimates=MAPPINGSmoothedEstimatePerScore
				# estimatedNumeratorDerivedFromScoresVAL=[referenceForNumEstimates[x] for x in scoresVAL]
				# would this make sense? aggregate_per_score(estimatedNumFromRelativitiesVal,ooBagnum,ooBagw,ooBagdenom)
				thisError = goodnessOfFit(unusedSamplePart, indicatedRelativityForApplyTree_reused, numerator, weight, denominator)
			else
			# no subsampling is performed
				thisError = ErrorStats() # Empty/Default Error Stats
			end
			res[i] = TreeWithErrorStats(deepcopy(thistree), deepcopy(thisError))
	end
	# return resFitted,res
	return res # todo tbd check this: currently we do not return fitted, but I think it is probably not needed anyway
end


text:
using Blink,WebSockets,HttpServer

function gui()
    wsh = WebSocketHandler() do req,client
        while true
            msg = read(client)
            println(String(msg))
        end
    end
    server = Server(wsh)
    hadsgewsid=rand(1:9999)
    @async run(server,hadsgewsid)

    W = Window()
    Blink.addstyle!(W,"""
    body {background-color:#2d2d2d;font-size: 10;color: white}
    .tooltip {
        position: relative;display: inline-block;border-bottom: 1px dotted black;}
    .tooltip .tooltiptext {visibility: hidden;width: 120px;background-color: black;color: #fff;text-align:center;padding: 5px 0;border-radius:6px;position:absolute;z-index:1;}
    .tooltip:hover .tooltiptext {visibility: visible;}
    .field-tip {position:relative;cursor:default}
    .field-tip .tip-content {position:absolute;top:-10px;right:9999px;width:200px;margin-right:-220px;padding:10px;color:#fff;background:#333;opacity:0;
        -webkit-transition:opacity 250ms ease-out;
           -moz-transition:opacity 250ms ease-out;
            -ms-transition:opacity 250ms ease-out;
             -o-transition:opacity 250ms ease-out;
                transition:opacity 250ms ease-out;}
    .field-tip .tip-content:before {
        content:' ';position:absolute;top:50%;
        left:-16px;width:0;height:0;margin-top:-8px;border:8px solid transparent;border-right-color:#333;}
    .field-tip:hover .tip-content {right:-20px;opacity:1;}
    """)
    body!(W,"""
    <form>

Equations

    <textarea id="equations" style="color:white;background-color:#4d4d4d;width:80%;height:50%">
    equations1
    </textarea>
    </form>
    <button id="submit" onclick="hadsgews.send('00')">submit</button>
    """)
    @js W resizeTo(683,400)
    @js W moveTo(0,5)

    Blink.addscript!(W,"""
    hadsgews = new WebSocket("ws://localhost:$hadsgewsid")
    """)

    return W,server
end

w,s=gui()

Blink.addscript!(w,"""
function addEq()
{
    var eqs = document.getElementById('eqs');
    var eq = document.createElement("textarea");
    eq.id = "eq"+eqs.childElementCount;
    eqs.appendChild(eq);
}
function addRow()
{
    var eqs = document.getElementById('eqs').children[0];
    var eq = document.createElement("tr");
    var n = eqs.childElementCount+1;
    eq.innerHTML = '<textarea id="eq'+n+'"style="width:100%">Eq.'+n+'</textarea>        	          <input size="5" type="text" name="polname'+n+'">
	        <input size="2" type="text" name="pollb'+n+'">
	        <input size="2" type="text" name="polub'+n+'">
	        <input size="10" type="text" name="polinit'+n+'">
	'
    eqs.appendChild(eq);
}

function remLastEq()
{
    var eqs = document.getElementById('eqs');
    eqs.removeChild(eqs.lastChild);
}
""")

body!(w,"""

Equations

        <textarea id="eq1" style="width:100%">Eq.1</textarea>

      <input size="5" type="text" name="polname1">

    <input size="2" type="text" name="pollb1">

    <input size="2" type="text" name="polub1">

    <input size="10" type="text" name="polinit1">

<button id="addeq" onclick="addEq()">add equation</button>
<button id="addr" onclick="addRow()">+</button>
<button id="addeq" onclick="remLastEq()">X</button>
<button id="submit" onclick="hadsgews.send('00')">submit</button>
""")

@js w document.getElementById("equations").value

if false
    close(w)
    close(s)
end

max_stars_repo_name: ZacLN/HADSGE

max_issues_repo_name: ZacLN/HADSGE

max_forks_repo_name: ZacLN/HADSGE


text:
using Plots
using MeshCat
using GeometryTypes
using CoordinateTransformations
using FileIO
using MeshIO
using Random

##########
## Maze ##
##########
integration = :rk4
r_quad = 3.0 # based on size of mesh file
model,obj_uncon = TrajectoryOptimization.Dynamics.quadrotor
n = model.n
m = model.m
N = 101
tf = 5.0
q0 = [1.;0.;0.;0.]

# -initial state
x0 = zeros(model.n)
x0[1:3] = [0.; 0.; 10.]
q0 = [1.;0.;0.;0.]
x0[4:7] = q0

# -final state
xf = copy(x0)
xf[1:3] = [0.;60.; 10.] # xyz position
xf[4:7] = q0

# -control, state limits
u_min = 0.0
u_max = 50.0
x_max = Inf*ones(model.n)
x_min = -Inf*ones(model.n)

x_max[1:3] = [25.0; Inf; 20]
x_min[1:3] = [-25.0; -Inf; 0.]

Q = (1.0)*Matrix(I,model.n,model.n)
R = (1.0)*Matrix(I,model.m,model.m)
Qf = (1000.0)*Matrix(I,model.n,model.n)

# obstacle constraints
r_cylinder = 2.
cylinders = []
zh = 3
l1 = 5
l2 = 4
l3 = 5
l4 = 10

for i = range(-25,stop=-10,length=l1)
    push!(cylinders,(i, 10,r_cylinder))
end

for i = range(10,stop=25,length=l1)
    push!(cylinders,(i, 10, r_cylinder))
end

for i = range(-7.5,stop=7.5,length=l3)
    push!(cylinders,(i, 30, r_cylinder))
end

for i = range(-25,stop=-10,length=l1)
    push!(cylinders,(i, 50, r_cylinder))
end

for i = range(10,stop=25,length=l1)
    push!(cylinders,(i, 50, r_cylinder))
end

for i = range(10+2*r_cylinder,stop=50-2*r_cylinder,length=l4)
    push!(cylinders,(-25, i, r_cylinder))
end

for i = range(10+2*r_cylinder,stop=50-2*r_cylinder,length=l4)
    push!(cylinders,(25, i, r_cylinder))
end

n_cylinders = length(cylinders)

function cI(c,x,u)
    for i = 1:n_cylinders
        c[i] = circle_constraint(x,cylinders[i][1],cylinders[i][2],cylinders[i][3]+r_quad)
    end
    c
end

opts = SolverOptions()
obj_uncon_maze = LQRObjective(Q, R, Qf, tf, x0, xf)
obj_con = TrajectoryOptimization.ConstrainedObjective(obj_uncon_maze,x_min=x_min,x_max=x_max,u_max=u_max,u_min=u_min,cI=cI)
solver_uncon = Solver(model,obj_uncon_maze,integration=integration,N=N,opts=opts)
solver_con = Solver(model,obj_con,integration=integration,N=N,opts=opts)
solver_con.opts.square_root = true
solver_con.opts.R_infeasible = 0.01
solver_con.opts.resolve_feasible = true
solver_con.opts.cost_tolerance = 1e-4
solver_con.opts.cost_tolerance_intermediate = 1e-4
solver_con.opts.constraint_tolerance = 1e-4
solver_con.opts.constraint_tolerance_intermediate = 0.01

solver_con.opts.penalty_scaling = 10.0
solver_con.opts.penalty_initial = 1.0
solver_con.opts.outer_loop_update_type = :feedback
solver_con.opts.iterations_outerloop = 25
solver_con.opts.iterations = 500
solver_con.opts.iterations_innerloop = 300
solver_con.opts.use_penalty_burnin = false
solver_con.opts.verbose = false
solver_con.opts.live_plotting = false

# Initial control trajectory
U_hover = 0.5*9.81/4.0*ones(solver_uncon.model.m, solver_uncon.N-1)

# Initial state trajectory
X_guess = zeros(model.n,7)
X_guess[:,1] = x0
X_guess[:,7] = xf
X_guess[1:3,2:6] .= [0 -12.5 -20 -12.5 0 ;15 20 30 40 45 ;10 10 10 10 10]

X_guess[4:7,:] .= q0
X0 = TrajectoryOptimization.interp_rows(N,solver_uncon.obj.tf,X_guess)

plot(X_guess[1:3,:]')

# Unconstrained solve
@time results_uncon, stats_uncon = solve(solver_uncon,U_hover)

# Constrained solve
@time results_con, stats_con = solve(solver_con,X0,U_hover)
@btime solve($solver_con,$X0,$U_hover)

# Dircol solve
dircol_options = Dict("tol"=>solver_con.opts.cost_tolerance,"constr_viol_tol"=>solver_con.opts.constraint_tolerance,"max_iter"=>6000)
@time results_dircol, stats_dircol = TrajectoryOptimization.solve_dircol(solver_con, to_array(results_con.X), [to_array(results_con.U) zeros(m)], options=dircol_options)

stats_con["iterations"]
stats_con["runtime"]
evals(solver_con,:f) / stats_con["iterations"]

stats_dircol["iterations"]
evals(solver,:f)/stats_dircol["iterations"]
stats_dircol["runtime"]

# Trajectory Plots
plot(to_array(results_uncon.U)',title="Quadrotor Unconstrained",xlabel="time",ylabel="control",labels="")
plot(to_array(results_uncon.X)[1:3,:]',title="Quadrotor Unconstrained",xlabel="time",ylabel="position",labels=["x";"y";"z"],legend=:topleft)

t_array = range(0,stop=solver_uncon.obj.tf,length=solver_uncon.N)
plot(t_array[1:N-1],to_array(results_uncon.U)',title="Quadrotor Maze",xlabel="time",ylabel="control",labels="")
plot(t_array,to_array(results_uncon.X)[1:3,:]',title="Quadrotor Maze",xlabel="time",ylabel="position",labels=["x";"y";"z"],legend=:topleft)

t_array = range(0,stop=solver_con.obj.tf,length=solver_con.N)
plot(t_array[1:end-1],to_array(results_con.U)',title="Quadrotor Maze",xlabel="time",ylabel="control",labels="")
plot(t_array,to_array(results_con.X)[1:3,:]',title="Quadrotor Maze",xlabel="time",ylabel="position",labels=["x";"y";"z"],legend=:topleft)
@assert max_violation(results_con) <= opts.constraint_tolerance
stats_con["iterations"]
# Constraint convergence plot
plot(stats_con["c_max"],yscale=:log10,title="Quadrotor Maze",xlabel="iteration",ylabel="log(max constraint violation)",label="sqrt",legend=:bottomleft)

# Constrained results
@show stats_con["iterations"]
@show stats_con["outer loop iterations"]
@show stats_con["c_max"][end]
@show stats_con["cost"][end]

#################
# Visualization #
#################

vis = Visualizer()
open(vis)

# Import quadrotor obj file
traj_folder = joinpath(dirname(pathof(TrajectoryOptimization)),"..")
urdf_folder = joinpath(traj_folder, "dynamics","urdf")
obj = joinpath(urdf_folder, "quadrotor_base.obj")

# color options
green_ = MeshPhongMaterial(color=RGBA(0, 1, 0, 1.0))
green_transparent = MeshPhongMaterial(color=RGBA(0, 1, 0, 0.1))
red_ = MeshPhongMaterial(color=RGBA(1, 0, 0, 1.0))
red_transparent = MeshPhongMaterial(color=RGBA(1, 0, 0, 0.1))
blue_ = MeshPhongMaterial(color=RGBA(0, 0, 1, 1.0))
blue_transparent = MeshPhongMaterial(color=RGBA(0, 0, 1, 0.1))
blue_semi = MeshPhongMaterial(color=RGBA(0, 0, 1, 0.5))
yellow_ = MeshPhongMaterial(color=RGBA(1, 1, 0, 1.0))
yellow_transparent = MeshPhongMaterial(color=RGBA(1, 1, 0, 0.75))

orange_ = MeshPhongMaterial(color=RGBA(233/255, 164/255, 16/255, 1.0))
orange_transparent = MeshPhongMaterial(color=RGBA(233/255, 164/255, 16/255, 0.1))
black_ = MeshPhongMaterial(color=RGBA(0, 0, 0, 1.0))
black_transparent = MeshPhongMaterial(color=RGBA(0, 0, 0, 0.1))
black_semi = MeshPhongMaterial(color=RGBA(0, 0, 0, 0.5))

function plot_cylinder(c1,c2,radius,mat,name="")
    geom = Cylinder(Point3f0(c1),Point3f0(c2),convert(Float32,radius))
    setobject!(vis["cyl"][name],geom,red_)
end

function addcylinders!(vis,cylinders,height=1.5)
    for (i,cyl) in enumerate(cylinders)
        plot_cylinder([cyl[1],cyl[2],0],[cyl[1],cyl[2],height],cyl[3],blue_,"cyl_$i")
        # plot_cylinder([cyl[1],cyl[2],0],[0,0,height],cyl[3],blue_,"cyl_$i")
    end
end

# geometries
robot_obj = FileIO.load(obj)

sphere_small = HyperSphere(Point3f0(0), convert(Float32,0.1*r_quad)) # trajectory points
sphere_medium = HyperSphere(Point3f0(0), convert(Float32,0.25*r_quad))

obstacles = vis["obs"]
traj = vis["traj"]
traj_x0 = vis["traj_x0"]
traj_uncon = vis["traj_uncon"]
target = vis["target"]
robot = vis["robot"]
robot_uncon = vis["robot_uncon"]

# Set camera location
settransform!(vis["/Cameras/default"], compose(Translation(0., 72., 60.),LinearMap(RotX(pi/7.5)*RotZ(pi/2))))
# settransform!(vis["/Cameras/default"], compose(Translation(0., 75., 50.),LinearMap(RotX(pi/10)*RotZ(pi/2))))
# settransform!(vis["/Cameras/default"], compose(Translation(0., 35., 65.),LinearMap(RotX(pi/3)*RotZ(pi/2))))

# Create and place obstacles
addcylinders!(vis,cylinders,16.0)

# Create and place trajectory
# for i = 1:N
#     setobject!(vis["traj_uncon"]["t$i"],sphere_small,blue_)
#     settransform!(vis["traj_uncon"]["t$i"], Translation(results_uncon.X[i][1], results_uncon.X[i][2], results_uncon.X[i][3]))
# end

# for i = 1:N
#     setobject!(vis["traj_x0"]["t$i"],sphere_small,blue_)
#     settransform!(vis["traj_x0"]["t$i"], Translation(X0[1,i], X0[2,i], X0[3,i]))
# end
for i = 1:size(X_guess,2)
    setobject!(vis["traj_x0"]["t$i"],sphere_medium,yellow_transparent)
    settransform!(vis["traj_x0"]["t$i"], Translation(X_guess[1:3,i]...))
end

for i = 1:N
    setobject!(vis["traj"]["t$i"],sphere_small,green_)
    settransform!(vis["traj"]["t$i"], Translation(results_con.X[i][1:3]...))
end

# setobject!(vis["robot_uncon"]["ball"],sphere_medium,orange_transparent)
# setobject!(vis["robot_uncon"]["quad"],robot_obj,black_)

# setobject!(vis["robot"]["ball"],sphere_medium,green_transparent)
setobject!(vis["robot"]["quad"],robot_obj,black_)

# Animate quadrotor
# for i = 1:N
#     settransform!(vis["robot_uncon"], compose(Translation(results_uncon.X[i][1], results_uncon.X[i][2], results_uncon.X[i][3]),LinearMap(Quat(results_uncon.X[i][4:7]...))))
#     sleep(solver_uncon.dt)
# end

for i = 1:N
    settransform!(vis["robot"], compose(Translation(results_con.X[i][1], results_con.X[i][2], results_con.X[i][3]),LinearMap(Quat(results_con.X[i][4:7]...))))
    sleep(solver_con.dt)
end

# Ghose quadrotor scene
traj_idx = [1;12;20;30;40;50;N]
n_robots = length(traj_idx)
for i = 1:n_robots
    robot = vis["robot_$i"]
    setobject!(vis["robot_$i"]["quad"],robot_obj,black_semi)
    settransform!(vis["robot_$i"], compose(Translation(results_con.X[traj_idx[i]][1], results_con.X[traj_idx[i]][2], results_con.X[traj_idx[i]][3]),LinearMap(quat2rot(results_con.X[traj_idx[i]][4:7]))))
end

# Animation
# (N-1)/tf
# anim = MeshCat.Animation(20)
# for i = 1:N
#     MeshCat.atframe(anim,vis,i) do frame
#         settransform!(frame["robot"], compose(Translation(results_con.X[i][1:3]...),LinearMap(Quat(results_con.X[i][4:7]...))))
#     end
# end
# MeshCat.setanimation!(vis,anim)


text:
# The HermiteFSE basis consists of the functions

struct HermiteFSE{T<:Real} <: Space{Line{false,T},T}
    t :: T
end

HermiteFSE{T}() where T = HermiteFSE{T}(zero(T))
HermiteFSE() = HermiteFSE{Float64}()

domain(::HermiteFSE{T}) where T = Line{false,T}()
canonicaldomain(::HermiteFSE{T}) where T = Line{false,T}()
normalization(::Type{T}, ::HermiteFSE, ::Int) where T = one(T)

spacescompatible(a::HermiteFSE,b::HermiteFSE) = a.t == b.t
canonicalspace(::HermiteFSE{T}) where T = HermiteFSE(zero(T))

function Derivative(sp::HermiteFSE,order::Integer)
    ConcreteDerivative(sp,order)
end

space(D::ConcreteDerivative{HFSE}) where HFSE <: HermiteFSE = D.space
rangespace(D::ConcreteDerivative{HFSE}) where HFSE <: HermiteFSE = domainspace(D)
bandwidths(D::ConcreteDerivative{HFSE}) where HFSE <: HermiteFSE = D.order,D.order
Base.stride(D::ConcreteDerivative{HFSE}) where HFSE <: HermiteFSE = D.order

function getindex(D::ConcreteDerivative{HermiteFSE{T},K,KK},k::Integer,j::Integer) where {T,K,KK}
    # differentiation matrix does not depend on HermiteFSE.t (see Kropielnicka, Iserles, Schratz, Webb 2021)

    m = D.order
    bandind = j-k

    if m == 1
        if bandind == 1
            convert(T,sqrt(k/2))
        elseif bandind == -1
            convert(T,-sqrt(j/2))
        else
            zero(T)
        end
    elseif m == 2
        if bandind == 0
            convert(T,.5-k)
        elseif bandind == 2
            convert(T,sqrt(k*(k+1))/2)
        elseif bandind == -2
            convert(T,sqrt(j*(j+1))/2)
        else
            zero(T)
        end
    else
        error("Higher order HermiteFSE derivatives are not supported yet")
    end
end

function evaluate(cfs::AbstractVector,Hspace::HermiteFSE,x)
    # Clenshaw's algorithm (working implementation given here for Hspace.t = 0) can have overflow/underflow problems:
    # T = eltype(x)
    # n = length(cfs)
    # bk2 = zero(T)
    # bk1 = convert(T,cfs[n])
    # for k = n-1:-1:1
    #     bk1,bk2 = convert(T,cfs[k]) + x*sqrt(2/k)*bk1 - sqrt(k/(k+1)) * bk2 , bk1
    # end
    # return bk1 * exp(-x^2 / 2) * convert(T,π)^(-1/4)

    # Rescaled forward recurrence to avoid overflow/underflow
    T = eltype(x)
    t = Hspace.t
    x = x/sqrt(1+4*t^2)
    n = length(cfs)-1
    cfs = cfs .* ((1+2*one(T)im*t)/(1-2*one(T)im*t)).^((0:n)*one(T)/2)
    ret = zero(T)
    if (n ≥ 0) ret += cfs[1] * exp(-x^2 / 2) * convert(T,π)^(-1/4) end
    if (n ≥ 1) ret += cfs[2] * x * exp.(-x^2 / 2) * (4/convert(T,π))^(1/4) end

    if (n ≥ 2)
        hkm1 = convert(T,π)^(-1/4)       #h_0(x)
        hk = sqrt(one(T)*2) * x * hkm1  #h_1(x)
        sum_log_scale = zero(T)
        for k = 2:n
            # hk = h_k(x), hkm1 = h_{k-1}(x) (actually, recaled versions thereof)
            hkm1, hk = hk, sqrt(one(T)*2/k)*x*hk - sqrt((k-one(T))/k)*hkm1
            # rescale values if necessary
            scale = (s->(s<one(T)*100) ? one(T) : inv(s))(abs(hk))
            hk *= scale
            hkm1 *= scale
            # keep track of final rescaling factor
            sum_log_scale += log(scale)
            ret += cfs[k+1] * hk * exp(-x^2 / 2 - sum_log_scale)
        end
    end
    ret *= exp(-one(T)im*t*x^2)/sqrt(1-2*one(T)im*t)
    return ret
end

eval_Hermite_function(n::Integer,x::T, t=0.0) where T <: Number = evaluate([zeros(T,n);one(T)],HermiteFSE(t),x)

## Transforms for HermiteFSE space
hasfasttransform(::HermiteFSE) = true
# Gauss-Hermite points
points(Hspace::HermiteFSE,n::Int) = gausshermite(n)[1] * sqrt(1+4*(Hspace.t)^2)

function plan_Hermite_transform(n::Integer)
   # builds the orthogonal matrix Q and weight vector 'valweights' such that
   # the val2coeffs transform is Q * (valweights .* vals) and
   # the coeffs2vals transform is (Q' * cfs) ./ valweights.

    Q = zeros(n+1,n+1)
    x = points(HermiteFSE(),n+1)

    hkm1 = ones(n+1)*π^(-1/4)       #h_0(x)
    hk = sqrt(2) * x .* hkm1  #h_1(x)
    cum_log_scale = zeros(n+1,n+1)

    Q[1,:] = hkm1
    if (n ≥ 1) Q[2,:] = hk end

    for k = 2:n
        # hk = h_k(x), hkm1 = h_{k-1}(x) (actually, recaled versions thereof)
        hkm1, hk = hk, sqrt(2/k)* x .* hk - sqrt((k-1)/k)*hkm1
        # rescale values if necessary
        scale = (x->(x<100) ? 1 : inv(x)).(abs.(hk))
        hk .*= scale
        hkm1 .*=  scale
        Q[k+1,:] = hk
        # keep track of rescaling factors
        cum_log_scale[k+1,:] = cum_log_scale[k,:] + log.(scale)
    end

    # the function values should be weighted by these
    valweights = exp.(cum_log_scale[n+1,:] + x.^2 / 2) ./ (sqrt(n+1)*abs.(Q[n+1,:]))

    # Q contains rescaled Hermite polynomials. Need to convert to unscaled ψ_k(x_j)/(sqrt(n+1)ψ_n(x_j)),
    # where ψ_n is the degree n Hermite function
    for k = 1:n+1
        Q[k,:] = Q[k,:]./abs.(Q[n+1,:]) .* exp.(cum_log_scale[n+1,:]-cum_log_scale[k,:] .- .5*log(n+1))
    end

    return valweights, Q
end

function plan_Hemite_transform_coeff_scaling(n::Integer, t = 0.0)
    # the coefficients should be weighted by these
    if t == 0.0
        coeffweights = ones(n+1)
    else
        x = points(HermiteFSE(t),n+1)
        coeffweights = ((1+2.0im*t)/(1-2.0im*t)).^(0:.5:n/2) .* exp.(-1.0im*t*x.^2)/sqrt(1-2.0im*t)
    end
    return coeffweights
end

function plan_transform(S::HermiteFSE,vals::AbstractVector)
    valweights, Q = plan_Hermite_transform(length(vals)-1)
    coeffweights = plan_Hemite_transform_coeff_scaling(length(vals)-1, S.t)
    TransformPlan(S,(valweights,Q,coeffweights),Val{false})
end
function plan_itransform(S::HermiteFSE,cfs::AbstractVector)
    valweights, Q = plan_Hermite_transform(length(cfs)-1)
    coeffweights = plan_Hemite_transform_coeff_scaling(length(cfs)-1, S.t)
    ITransformPlan(S,(valweights,Q,coeffweights),Val{false})
end

function *(P::TransformPlan{T,S,false},vals::AbstractVector) where {T,S<:HermiteFSE}
    valweights, Q, coeffweights = P.plan
    (Q*(valweights.*vals)) ./ coeffweights
end
function *(P::ITransformPlan{T,S,false},cfs::AbstractVector) where {T,S<:HermiteFSE}
    valweights, Q, coeffweights = P.plan
    (Q' * (coeffweights.*cfs)) ./ valweights
end