text:
using DynamicPolynomials, LinearAlgebra

import DynamicPolynomials.MultivariatePolynomials: mapcoefficientsnz

@inline function swap!(A::Matrix{Int}, i::Int, j::Int)
    tmp = Array{Int}(undef, size(A, 1))
    copyto!(tmp, view(A, :, i))
    copyto!(view(A, :, i), view(A, :, j))
    copyto!(view(A, :, j), tmp)
end

function elim!(A::Matrix{Int})
    m, n = size(A)
    h, k = one(Int), one(Int)
    while h ≤ n && k ≤ m
        x, i = findmax(A[k,h:n])
        i += h - one(h)
        if iszero(x)
            k += one(k)
        else
            swap!(A, i, h)
            for i in 1:n
                if i != h
                    f = A[k, i] ÷ A[k, h]
                    A[k, i] = zero(A[k, i])
                    @simd for u in k+1:m
                        A[u,i] ⊻= f * A[u,h]
                    end
                end
            end
            h += 1
            k += 1
        end
    end
    A
end

struct Network
    table::Vector{Tuple{UInt, Set{UInt}}}
end

Base.length(net::Network) = 2^length(net.table)

S_POMBE = Network([
    (0b000000000, Set{UInt}()),
    (0b000011100, Set([0b000000000])),
    (0b100100111, Set([0b100000000,
                       0b000000100,
                       0b100000100,
                       0b100100100,
                       0b100000110,
                       0b100000101
                      ])),
    (0b100101011, Set([0b100000000,
                       0b000001000,
                       0b100001000,
                       0b100101000,
                       0b100001010,
                       0b100001001
                      ])),
    (0b000100000, Set([0b000100000])),
    (0b011011100, Set([0b010000000])),
    (0b101000010, Set([0b100000000,
                       0b001000000,
                       0b101000000,
                       0b101000010
                      ])),
    (0b110000010, Set([0b010000000,
                       0b000000010,
                       0b010000010,
                       0b110000010
                      ])),
    (0b000010000, Set([0b000010000]))
])

NET = Network([
    (0b111, Set([0b000, 0b010, 0b100, 0b110])),
    (0b111, Set([0b010, 0b100, 0b110, 0b111])),
    (0b111, Set([0b011, 0b101, 0b110, 0b111]))
])

Base.eltype(net::Network) = UInt

Base.iterate(net::Network) = (UInt(0b0), 1)
function Base.iterate(net::Network, state::Int)
    if state ≥ length(net)
        nothing
    else
        UInt(state), state + 1
    end
end

Base.in(net::Network, state::Unsigned) = zero(state) ≤ state < length(net)

function update(net::Network, state::Unsigned)
    if !(state in net)
        error("$state is not in the network")
    end
    next = UInt(0b0)
    for (i, (mask, activations)) in enumerate(net.table)
        next ⊻= UInt((state & mask) in activations) << (i - 1)
    end
    next
end

decode(x, n) = Int[1 & (x >> (i-1)) for i in 1:n]

function transitions(net::Network)
    trans = Array{UInt}(undef, length(net))
    for state in net
        trans[state + 1] = update(net, state) + 1
    end
    trans
end

function basins(net::Network)
    trans = transitions(net)
    visited = falses(length(trans))
    basins = zeros(Int, length(trans))
    basin_number = 1
    initial_state = UInt(0b1)
    while initial_state ≤ length(trans)
        stack = UInt[]
        state = initial_state
        next_state = trans[state]
        visited[state] = true
        while !visited[next_state]
            push!(stack, state)
            state = next_state
            next_state = trans[state]
            visited[state] = true
        end
        basin = if basins[next_state] == zero(Int)
            (basin_number += 1) - 1
        else
            basin = basins[next_state]
        end

        basins[state] = basin
        while !isempty(stack)
            next_state = state
            state = pop!(stack)
            basins[state] = basin
        end

        while initial_state ≤ length(visited) && visited[initial_state]
            initial_state += 1
        end
    end
    [UInt.(findall(basins .== i) .- 1) for i in 1:basin_number-1]
end

struct Space
    n::Int
    v::Int
    Space(n) = n < 0 ? error("invalid dimension") : new(n, 2^n - 1)
end

Base.iterate(space::Space) = trues(space.n), (0, trues(space.n))

function Base.iterate(space::Space, state)
    n, s = state
    if n >= space.v - 1
        nothing
    else
        for i in 1:space.n
            if s[i]
                s[i] = zero(s[i])
                s[1:i-1] .= one(s[i])
                break
            end
        end
        copy(s), (n + 1, s)
    end
end

Base.length(space::Space) = space.v
Base.eltype(space::Space) = BitArray{1}

function invarianteq(state::Unsigned, n::Int)
    space = Space(n)
    bits = ones(Int, space.v)
    for (i, s) in enumerate(space)
        for (j, x) in enumerate(s)
            if x
                bits[i] &= (1 & state >> (j - 1))
            end
            if bits[i] == zero(bits[i])
                break
            end
        end
    end
    bits
end

function system(net::Network)
    N = length(net.table)
    A = zeros(Int, length(net) - 1, length(net))
    for (i, state) in enumerate(net)
        next = update(net, state)
        A[:,i] .= invarianteq(state, N) .⊻ invarianteq(next, N)
    end
    A
end

function prune(A::AbstractMatrix{Int})
    j = findfirst(collect(all(A[:,j] .== zero(Int)) for j in 1:size(A,2))) - 1
    @view A[:,1:j]
end

function dependent(B::AbstractMatrix{Int})
    T = Array(B')
    for j in 1:size(B,2)
        i = findfirst(B[:, j] .== one(Int)) + 1
        T[j, i:end] .= zero(Int)
    end
    T
end

function independent(T::AbstractMatrix{Int})
    idx = Int[]
    for j in 1:size(T,2)
        if all(T[:,j] .== zero(Int))
            push!(idx, j)
        end
    end
    U = zeros(Int, length(idx), size(T,2))
    for (i, j) in enumerate(idx)
        U[i, j] = 1
    end
    U, idx
end

function solve(A::Matrix{Int})
    @polyvar a[1:size(A,1)]
    B = prune(A)
    T = dependent(B)
    U, idx = independent(T)
    a, a .=> (a[idx]' * U * (I - B*T))'
end

const Factors = Dict{PolyVar{true}, Polynomial{true, Int}}

function freevars(xs::AbstractVector{PolyVar{true}}, terms::AbstractVector{Polynomial{true, Int}})
    sort!(setdiff(union(effective_variables.(terms)...), xs), rev=true)
end

function groupterms(xs, terms)
    factors = Factors()
    free = freevars(xs, terms)
    for v in free
        factors[v] = zero(keytype(factors))
    end
    for term in filter(t -> t != 0, terms)
        for m in monomials(term)
            vs = intersect(free, effective_variables(m))
            @assert length(vs) == 1
            factors[vs[1]] += m
        end
    end
    for (k, v) in factors
        factors[k] = mapcoefficientsnz(c -> c % 2, v)
    end
    factors
end

const SolSub = Pair{PolyVar{true}, Polynomial{true, Int}}

struct Invariant
    x::Vector{PolyVar{true}}
    a::Vector{PolyVar{true}}
    factors::Factors
end

function invariantterms(as::AbstractVector{PolyVar{true}})
    n = Int(log2(length(as) + 1))
    space = Space(n)
    @assert length(as) == length(space)
    @polyvar x[1:n]
    terms = zeros(Polynomial{true, Int}, space.v)
    for (i, (a, t)) in enumerate(zip(as, space))
        terms[i] = a * prod(x[t])
    end
    x, terms, as
end

function Invariant(a::AbstractVector{PolyVar{true}}, ss::AbstractVector{SolSub})
    x, terms = invariantterms(a)
    factors = groupterms(x, map(t -> subs(t, ss...), terms))
    Invariant(x, a, factors)
end

function Invariant(net::Network)
    A = system(net)
    elim!(A)
    invariant(solve(A)...)
end

function evaluate(ℰ::Invariant, state)
    E = zero(keytype(ℰ.factors))
    ss = ℰ.x .=> decode(state, length(ℰ.x))
    for p in values(ℰ.factors)
        E += subs(p, ss...)
    end
    mapcoefficientsnz(c -> c % 2, E)
end

Base.eltype(::Invariant) = Invariant
Base.length(::Invariant) = 1
Base.iterate(ℰ::Invariant) = (ℰ,1)
Base.iterate(ℰ::Invariant,x) = (x == 0) ? (ℰ,1) : nothing

max_stars_repo_name: ELIFE-ASU/dglmoore

max_issues_repo_name: ELIFE-ASU/dglmoore

max_forks_repo_name: ELIFE-ASU/dglmoore


text:
@testset "Pair Interactions" begin

    # Test that each exchange gets mapped to the correct backend type
    function test_type_mapping()
        cryst = Sunny.diamond_crystal()
        latsize = (4, 4, 4)
        exchange_ints = diamond_test_exchanges()

        # Note: These tests rely on diamond_test_exchanges outputting
        #  the interactions in the order
        #  [heisenberg, on_site [heisen], diagonal, general]
        correct_types = [
            Sunny.HeisenbergCPU,
            Sunny.HeisenbergCPU,
            Sunny.DiagonalCouplingCPU,
            Sunny.GeneralCouplingCPU
        ]

        for (int, type) in zip(exchange_ints, correct_types)
            sys = SpinSystem(cryst, [int], latsize)
            ℋ = sys.hamiltonian
            if type == Sunny.HeisenbergCPU
                @test length(ℋ.heisenbergs) == 1
            elseif type == Sunny.DiagonalCouplingCPU
                @test length(ℋ.diag_coups) == 1
            else
                @test length(ℋ.gen_coups) == 1
            end
        end
    end

    test_type_mapping()

    function Base.isapprox(int::Sunny.QuadraticInteraction, int2::Sunny.QuadraticInteraction; kwargs...)
        exchange_eq = isapprox(int.J, int2.J; kwargs...)
        bond_eq = int.bond == int2.bond
        exchange_eq && bond_eq
    end

    # Test that the output of `show` on pair interactions actually creates the correct interaction again,
    #   up to truncation error in `show`. (Does not check label on interactions, which isn't printed in `show`)
    function test_quadratic_shown_info()
        io = IOBuffer()
        context = IOContext(io)
        exchange_ints = diamond_test_exchanges()
        for int in exchange_ints
            show(context, "text/plain", int)
            output = String(take!(io))
            reparsed_int = eval(Meta.parse(output))
            @test isapprox(int, reparsed_int; atol=1e-4)
        end
    end

    test_quadratic_shown_info()

    # Checks that a list of bonds contains no equivalent bonds
    function no_equivalent_bonds(bonds::Vector{Bond})
        for (i, bond) in enumerate(bonds[1:end-1])
            equiv_bond = Bond(bond.j, bond.i, -bond.n)
            rest_of_bonds = bonds[i+1:end]
            if (bond ∈ rest_of_bonds) || (equiv_bond ∈ rest_of_bonds)
                return false
            end
        end
        return true
    end

    # Test that iterating over bonds on specific sublattices both only gives
    #  bonds of that sublattice, but also that we hit all bonds if we iterate
    #  over each sublattice.
    function correct_sublat_iteration(bondtable::Sunny.BondTable)
        total_iterated = 0
        for b in 1:nbasis(bondtable)
            for (bond, _) in Sunny.sublat_bonds(bondtable, b)
                if !(bond.i == b || bond.j == b)
                    return false
                end
                total_iterated += 1
            end
        end
        return total_iterated == length(bondtable)
    end

    @testset "BondTable" begin
        cryst = Sunny.diamond_crystal()
        latsize = (4, 4, 4)
        exchange_ints = diamond_test_exchanges()
        sys = SpinSystem(cryst, exchange_ints, (4,4,4))
        ℋ = sys.hamiltonian
        pair_ints = [ℋ.heisenbergs..., ℋ.diag_coups..., ℋ.gen_coups...]

        for pair_int in pair_ints
            bondtable = pair_int.bondtable
            @test no_equivalent_bonds(bondtable.culled_bonds)
            @test correct_sublat_iteration(bondtable)
        end
    end

end

ext: jl


text:
function hfun_bar(vname)
  val = Meta.parse(vname[1])
  return round(sqrt(val), digits=2)
end

function hfun_m1fill(vname)
  var = vname[1]
  return pagevar("index", var)
end

function lx_baz(com, _)
  # keep this first line
  brace_content = Franklin.content(com.braces[1]) # input string
  # do whatever you want here
  return uppercase(brace_content)
end

function lx_google(com,_)

end

function lx_figpop(lxc::Franklin.LxCom, _)
  rpath = Franklin.stent(lxc.braces[1])
  path  = Franklin.parse_rpath(rpath; canonical=false, code=true)
  fdir, fext = splitext(path)
  alt = ""
  # there are several cases
  # A. a path with no extension --> guess extension
  # B. a path with extension --> use that
  # then in both cases there can be a relative path set but the user may mean
  # that it's in the subfolder /output/ (if generated by code) so should look
  # both in the relpath and if not found and if /output/ not already last dir
  candext = ifelse(isempty(fext),
                   (".png", ".jpeg", ".jpg", ".svg", ".gif"), (fext,))
  for ext ∈ candext
      candpath = fdir * ext
      syspath  = joinpath(Franklin.PATHS[:site], split(candpath, '/')...)
      isfile(syspath) && return Franklin.html_img(candpath, alt)
  end
  # now try in the output dir just in case (provided we weren't already
  # looking there)
  p1, p2 = splitdir(fdir)
  if splitdir(p1)[2] != "output"
      for ext ∈ candext
          candpath = p1 * "/output/" *p2 * ext
          syspath  = joinpath(Franklin.PATHS[:site], split(candpath, '/')...)
          isfile(syspath) && return "@@parent-container [\\fig{$rpath}]($candpath) @@"
      end
  end
  return Franklin.html_err("Image matching '$path' not found.")
end


text:
struct RadialVelocities{T<:AbstractFloat}
    times::Vector{T}
    rvs::Vector{T}
    s_copy::State{T}
end

function RadialVelocities(times::Vector{T}, ic::ElementsIC{T}) where T<:AbstractFloat
    rvs = zeros(T, length(times))
    s_copy = State(ic)
    return RadialVelocities(times, rvs, s_copy)
end

function (intr::Integrator)(s::State{T}, rvs::RadialVelocities{T}) where T<:AbstractFloat

    t0 = s.t[1] #initial time
    nsteps = intr.tmax/intr.h
    h = intr.h

    s_copy = rvs.s_copy #shortcut to structure(s_copy)

    tc = t0
    for (i,t) in enumerate(rvs.times)

        dt = t-tc
        while dt > h
            intr.scheme(s,h)
            tc += h
            dt = t-tc
        end

        set_state!(s_copy, s) #copying data in 's' to 's_copy'
        intr.scheme(s_copy,dt)
        rvs.rvs[i] = s_copy.v[3,1]

    end
end


text:
# MacOS tests

using Test
using Pkg
using ApplicationBuilder

builddir = mktempdir()
@assert isdir(builddir)

@testset "make_bundle_identifier Utils" begin
    @test occursin(r"""^com.[a-z0-9]+.myappnamedthisapp22$""",
                ApplicationBuilder.make_bundle_identifier("My app named this_app22")
          )
end

@testset "copy files Utils" begin
    @testset "clean_file_pattern" begin
        @test ApplicationBuilder.clean_file_pattern("/a/b/c", "") == "/a/b/c"
        @test ApplicationBuilder.clean_file_pattern(" ~/tmp/a.txt ", "") == "$(homedir())/tmp/a.txt"
        @test ApplicationBuilder.clean_file_pattern("Hello World", "") == "Hello World"
        @test_throws ArgumentError ApplicationBuilder.clean_file_pattern(" ", "-R")
        try
            ApplicationBuilder.clean_file_pattern(" ", "-R")
        catch e
            @test occursin(r".*-R.*", e.msg)
        end
    end
    @testset "copy_file_dir_or_glob" begin
        d = mkpath(joinpath(mktempdir(), "d"))
        mkpath(joinpath(d, "a"))
        for i in 1:5
            write(joinpath(d, "a", "f$i.txt"), rand(1:100))
        end
        for i in 1:5
            write(joinpath(d, "a", "f$i.bi"), rand(1:100))
        end

        mkpath(joinpath(d, "b"))
        ApplicationBuilder.copy_file_dir_or_glob("$d/a/*.txt", "$d/b")
        @test length(readdir(joinpath(d,"b"))) == 5

        rm(joinpath(d, "b"), force=true, recursive=true)
        mkpath(joinpath(d, "b"))
        ApplicationBuilder.copy_file_dir_or_glob("$(relpath(d))/a", "$d/b")
        @test length(readdir(joinpath(d,"b"))) == 1
        @test length(readdir(joinpath(d,"b","a"))) == 10
    end
end

@testset "HelloWorld.app" begin
    @test 0 == include("build_examples/hello.jl")
    @test isdir("$builddir/HelloWorld.app")
    @test success(`$builddir/HelloWorld.app/Contents/MacOS/hello`)

    # There shouldn't be a Libraries dir since none specified.
    @test !isdir("$builddir/HelloWorld.app/Contents/Libraries")

    # Ensure all dependencies on Julia libs are internal, so the app is portable.
    @testset "No external Dependencies" begin
        @test !success(pipeline(
                        `otool -l "$builddir/HelloWorld.app/Contents/MacOS/hello"`,
                        `grep 'julia'`,  # Get all julia deps
                        `grep -v '@rpath'`))  # make sure all are relative.
    end
end

@testset "commandline_app" begin
    @test 0 == include("build_examples/commandline_hello.jl")
    @test success(`open $builddir/hello.app`)
end

function testRunAndKillProgramSucceeds(cmd, timeout=10)
    p = open(cmd, "r+") # Make sure it runs correctly
    sleep(1)
    process_exited(p) && (println("Test Failed: failed to launch: \n", read(p.out, String)); return false)
    sleep(timeout)
    process_exited(p) && (println("Test Failed: Process died: \n", read(p.out, String)); return false)
    # Manually kill program after it's been running for a bit.
    kill(p); sleep(1)
    return true
end

# Test that it can run without .julia directory (Dangerous!)
function testBundledSuccessfully_macro(cmd_expr, timeout=10)
    quote
        val = false
        julia_dir = Pkg.Pkg2._pkgroot()
        mv(julia_dir, julia_dir*".bak")  # NOTE: MUST mv() THIS BACK
        try
            val = testRunAndKillProgramSucceeds($cmd_expr, $timeout)
        catch
        end
        mv(julia_dir*".bak", julia_dir)  # NOTE: MUST RUN THIS LINE IF ABOVE IS RUN
        val
    end
end
macro testBundledSuccessfully(expr...)
    testBundledSuccessfully_macro(expr...)
end

# Disabling the SDL tests since Cairo is currently broken in METADATA.
#@testset "sdl: simple example of binary dependencies" begin
#    @test 0 == include("build_examples/sdl.jl")
#    # Test that it runs correctly
#    @test testRunAndKillProgramSucceeds(`$builddir/HelloSDL2.app/Contents/MacOS/sdl`)
#    # Test that it can run without .julia directory
#    @test @testBundledSuccessfully(`$builddir/HelloSDL2.app/Contents/MacOS/sdl`, 3)
#end

# Disabling Blink Tests since Blink has changed and this no longer works.
#@testset "HelloBlink.app" begin
#    @test 0 == include("build_examples/blink.jl")
#
#    @test isdir("$builddir/HelloBlink.app")
#    # Test that it copied the correct files
#    @test isdir("$builddir/HelloBlink.app/Contents/Libraries")
#    @test isfile("$builddir/HelloBlink.app/Contents/Resources/main.js")
#    # Test that it runs correctly
#    @test testRunAndKillProgramSucceeds(`$builddir/HelloBlink.app/Contents/MacOS/blink`)
#    # Test that it can run without .julia directory
#
#    # TODO: This is broken because Blink currently can't be statically compiled
#    # https://github.com/JunoLab/Blink.jl/pull/174
#    # (It appears to work in this test, but the application does nothing because it errors.)
#    #  @test @testBundledSuccessfully(`$builddir/HelloBlink.app/Contents/MacOS/blink`, 10)
#    # Replacing with a test_broken so we remember.
#    @test_broken false
#end


text:
function OsculatingToMeanElementsConverter(arg0::SpacecraftState, arg1::jint, arg2::Propagator, arg3::jdouble)
    return OsculatingToMeanElementsConverter((SpacecraftState, jint, Propagator, jdouble), arg0, arg1, arg2, arg3)
end

function convert(obj::OsculatingToMeanElementsConverter)
    return jcall(obj, "convert", SpacecraftState, ())
end


text:
using Pkg

Pkg.add(PackageSpec(url="https://github.com/davpoolechem/JuliaChem.jl"))
Pkg.instantiate()


text:
# ---
# title: 474. Ones and Zeroes
# id: problem474
# author: Indigo
# date: 2021-06-22
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/ones-and-zeroes/description/>
# hidden: true
# ---
#
# You are given an array of binary strings `strs` and two integers `m` and `n`.
#
# Return _the size of the largest subset of`strs` such that there are **at
# most** _`m` __`0` _' s and _`n` __`1` _' s in the subset_.
#
# A set `x` is a **subset** of a set `y` if all elements of `x` are also
# elements of `y`.
#
#
#
# **Example 1:**
#
#
#
#     Input: strs = ["10","0001","111001","1","0"], m = 5, n = 3
#     Output: 4
#     Explanation: The largest subset with at most 5 0's and 3 1's is {"10", "0001", "1", "0"}, so the answer is 4.
#     Other valid but smaller subsets include {"0001", "1"} and {"10", "1", "0"}.
#     {"111001"} is an invalid subset because it contains 4 1's, greater than the maximum of 3.
#
#
# **Example 2:**
#
#
#
#     Input: strs = ["10","0","1"], m = 1, n = 1
#     Output: 2
#     Explanation: The largest subset is {"0", "1"}, so the answer is 2.
#
#
#
#
# **Constraints:**
#
#   * `1 <= strs.length <= 600`
#   * `1 <= strs[i].length <= 100`
#   * `strs[i]` consists only of digits `'0'` and `'1'`.
#   * `1 <= m, n <= 100`
#
#
## @lc code=start
using LeetCode

function find_max_form(strs::Vector{String}, m::Int, n::Int)
    dp = OffsetArray(fill(0, m + 1, n + 1), -1, -1)
    for str in strs
        o = count(==('1'), str)
        z = length(str) - o
        for j in m:-1:z, k in n:-1:o
            dp[j, k] = max(dp[j, k], dp[j - z, k - o] + 1)
        end
    end
    return dp[end]
end
## @lc code=end


text:
# HER2K / SKR2K
#

macro blis_linalg_her2k(typename, btypename, targetfunc, bliname)

    # Get method for the typed API backend.
    blifuncname = Symbol("bli_", cblas_typechar[ eval(typename) ], bliname)
    blifunc = getproperty(TypedBackend, blifuncname)

    return quote
        """
        BLIS-based HER2K/SYR2K with approximately the same interface as `LinearAlgebra.BLAS`.
        Strided arrays are directly supported without further Julia-level dispatching.
        """
        $(esc(targetfunc))(uplo::AbstractChar,
                           tAB::AbstractChar,
                           alpha::$typename,
                           A::StridedMatrix{$typename},
                           B::StridedMatrix{$typename},
                           beta::$btypename,
                           C::StridedMatrix{$typename}) = begin

            bli_uploc = char_to_uplo[uplo]
            bli_tAB = char_to_trans[tAB]

            n , k = size(A)
            n_, k_= size(B)
            m , m_= size(C)

            if bli_tAB == BLIS_TRANSPOSE || bli_tAB == BLIS_CONJ_TRANSPOSE
                n , k  = k , n
                n_, k_ = k_, n_
            end

            (m == m_) || throw(DimensionMismatch("HE(SY)R2K C storage not square."))
            (k == k_) || throw(DimensionMismatch("HE(SY)R2K contracted size mismatch."))
            (m == n &&
             n == n_) || throw(DimensionMismatch("HE(SY)R2K external size mismatch."))

            $blifunc(bli_uploc,
                     bli_tAB,
                     m, k,
                     [alpha],
                     A, strides(A)...,
                     B, strides(B)...,
                     [beta],
                     C, strides(C)...)
            C

        end
    end
end

@blis_linalg_her2k(Float32, Float32, her2k!, her2k!)
@blis_linalg_her2k(Float64, Float64, her2k!, her2k!)
@blis_linalg_her2k(ComplexF32, Float32, her2k!, her2k!)
@blis_linalg_her2k(ComplexF64, Float64, her2k!, her2k!)

@blis_linalg_her2k(Float32, Float32, syr2k!, syr2k!)
@blis_linalg_her2k(Float64, Float64, syr2k!, syr2k!)
@blis_linalg_her2k(ComplexF32, ComplexF32, syr2k!, syr2k!)
@blis_linalg_her2k(ComplexF64, ComplexF64, syr2k!, syr2k!)


text:
using ADNLPModels
using Documenter
using Printf
using StoppingInterface

makedocs(
  modules = [StoppingInterface],
  doctest = true,
  # linkcheck = true,
  strict = true,
  format = Documenter.HTML(
    assets = ["assets/style.css"],
    prettyurls = get(ENV, "CI", nothing) == "true",
  ),
  sitename = "StoppingInterface.jl",
  pages = ["Home" => "index.md", "Tutorial" => "tutorial.md", "Reference" => "reference.md"],
)

deploydocs(repo = "github.com/SolverStoppingJulia/StoppingInterface.jl.git", devbranch = "main")


text:
function assert_enoughvaluesforepochs(vname::String, v::Vector, epochs::Int)
   if length(v) < epochs
      error("Not enough `$vname` (vector of length $(length(v))) " .*
         "for training epochs ($epochs)")
   end
end

function assertinitoptimizers(optimizer::AbstractOptimizer,
      optimizers::Vector{<:AbstractOptimizer}, bm::BM,
      learningrates::Vector{Float64}, sdlearningrates::Vector{Float64},
      epochs::Int
      ) where {BM<:AbstractBM}

   if isempty(optimizers)
      if optimizer === NoOptimizer()
         # default optimization algorithm
         optimizers = map(
               i -> defaultoptimizer(bm, learningrates[i], sdlearningrates[i]),
               1:epochs)
      else
         optimizers = fill(initialized(optimizer, bm), epochs)
      end
   else
      optimizers = map(opt -> initialized(opt, bm), optimizers)
      assert_enoughvaluesforepochs("optimizers", optimizers, epochs)
   end

   optimizers
end

function defaultoptimizer(rbm::R,
      learningrate::Float64, sdlearningrate::Float64) where {R<: AbstractRBM}

   loglikelihoodoptimizer(rbm; learningrate = learningrate,
         sdlearningrate = sdlearningrate)
end

function defaultoptimizer(dbm::MultimodalDBM, learningrate::Float64,
      sdlearningrate::Float64)

   StackedOptimizer(map(
         rbm -> defaultoptimizer(rbm, learningrate, sdlearningrate),
         dbm))
end

"""
    fitrbm(x; ...)
Fits an RBM model to the data set `x`, using Stochastic Gradient Descent (SGD)
with Contrastive Divergence (CD), and returns it.

# Optional keyword arguments (ordered by importance):
* `rbmtype`: the type of the RBM that is to be trained
   This must be a subtype of `AbstractRBM` and defaults to `BernoulliRBM`.
* `nhidden`: number of hidden units for the returned RBM
* `epochs`: number of training epochs
* `learningrate`/`learningrates`: The learning rate for the weights and biases
   can be specified as single value, used throughout all epochs, or as a vector
   of `learningrates` that contains a value for each epoch. Defaults to 0.005.
* `batchsize`: number of samples that are used for making one step in the
   stochastic gradient descent optimizer algorithm. Default is 1.
* `pcd`: indicating whether Persistent Contrastive Divergence (PCD) is to
   be used (true, default) or simple CD that initializes the Gibbs Chain with
   the training sample (false)
* `cdsteps`: number of Gibbs sampling steps for (persistent)
   contrastive divergence, defaults to 1
* `monitoring`: a function that is executed after each training epoch.
   It takes an RBM and the epoch as arguments.
* `upfactor`, `downfactor`: If this function is used for pretraining a part of
   a DBM, it is necessary to multiply the weights of the RBM with factors.
* `sdlearningrate`/`sdlearningrates`: learning rate(s) for the
   standard deviation if training a `GaussianBernoulliRBM` or
   `GaussianBernoulliRBM2`. Ignored for other types of RBMs.
   It usually must be much smaller than the learning rates for
   the weights. By default it is 0.0, which means that the standard deviation
   is not learned.
* `startrbm`: start training with the parameters of the given RBM.
   If this argument is specified, `nhidden` and `rbmtype` are ignored.
* `optimizer`/`optimizers`: an object of type `AbstractOptimizer` or a vector of
   them for each epoch. If specified, the optimization is performed as implemented
   by the given optimizer type. By default, the `LoglikelihoodOptimizer`
   with the `learningrate`/`learningrates` and `sdlearningrate`/`sdlearningrates`
   is used. For other types of optimizers, the learning rates must be specified
   in the `optimizer`. For more information on how to write your own optimizer,
   see `AbstractOptimizer`.
"""
function fitrbm(x::Matrix{Float64};
      nhidden::Int = size(x,2),
      inversedispersion::AbstractArray = ones(size(x,2)) .* 0.5,
      epochs::Int = 10,
      upfactor::Float64 = 1.0,
      downfactor::Float64 = 1.0,
      learningrate::Float64 = 0.005,
      learningrates::Vector{Float64} = fill(learningrate, epochs),
      pcd::Bool = true,
      cdsteps::Int = 1,
      batchsize::Int = 1,
      fisherscoring::Int = 1,
      dispersionquantile::Float64 = 0.5,
      estimatedispersion::String = "gene",
      lambda::Float64 = 100.0,
      zeroinflation::Bool = true,
      dropoutmid::Float64 = 0.5,
      dropoutshape::Float64 = 0.5,
      limit::Int = 10,
      epsilon::Float64 = 0.0001,
      rbmtype::DataType = BernoulliRBM,
      startrbm::AbstractRBM = NoRBM(),
      monitoring::Function = nomonitoring,

      # these arguments are only relevant for GaussianBernoulliRBMs:
      sdlearningrate::Float64 = 0.0,
      sdlearningrates::Vector{Float64} = fill(sdlearningrate, epochs),
      sdinitfactor::Float64 = 0.0,

      optimizer::AbstractOptimizer = NoOptimizer(),
      optimizers::Vector{<:AbstractOptimizer} = Vector{AbstractOptimizer}(),
      sampler::AbstractSampler = (cdsteps == 1 ? NoSampler() : GibbsSampler(cdsteps - 1)))

   if startrbm === NoRBM()
      rbm = initrbm(x, nhidden, inversedispersion, fill(dropoutmid, size(x,2)), fill(dropoutshape, size(x,2)), rbmtype)
   else
      rbm = deepcopy(startrbm)
      nhidden = nhiddennodes(startrbm)
   end

   assert_enoughvaluesforepochs("learningrates", learningrates, epochs)
   assert_enoughvaluesforepochs("sdlearningrates", sdlearningrates, epochs)

   if sdinitfactor > 0 &&
         (rbmtype == GaussianBernoulliRBM || rbmtype == GaussianBernoulliRBM2)
      rbm.sd .*= sdinitfactor
   end

   optimizers = assertinitoptimizers(optimizer, optimizers, rbm,
         learningrates, sdlearningrates, epochs)

   if pcd
      chainstate = rand(batchsize, nhidden)
   else
      chainstate = Matrix{Float64}(undef, 0, 0)
   end

   # allocate space for trainrbm!
   nvisible = size(x, 2)
   nsamples = size(x,1)
   h = Matrix{Float64}(undef, batchsize, nhidden)
   hmodel = Matrix{Float64}(undef, batchsize, nhidden)
   vmodel = Matrix{Float64}(undef, batchsize, nvisible)

   for epoch = 1:epochs
      if epoch % 200 == 0
         println("Epoch = ", epoch)
      end
      # Train RBM on data set
      trainrbm!(rbm, x, cdsteps = cdsteps, chainstate = chainstate,
            upfactor = upfactor, downfactor = downfactor,
            learningrate = learningrates[epoch],
            sdlearningrate = sdlearningrates[epoch],
            fisherscoring = fisherscoring,
            zeroinflation = zeroinflation,
            optimizer = optimizers[epoch],
            sampler = sampler,
            batchsize = batchsize,
            h = h, hmodel = hmodel, vmodel = vmodel)

      # Evaluation of learning after each training epoch
      monitoring(rbm, epoch)

      # Estimate inverse dispersion using Fisher scoring algorithm
      nfisherscoring = fisherscoring
      if fisherscoring != 0 && epoch % 5 == 0

         if estimatedispersion == "gene"
            genewisedispersion = zeros(0)
            wts = zeros(0)
               for i in 1:size(x,2)
                  push!(genewisedispersion,mle_for_θ(x[:,i], estimatedmean[:,i], wts, rbm, lambda; maxIter=20))
               end

            rbm.inversedispersion .= genewisedispersion

         elseif estimatedispersion == "gene-cell"

            genewisedispersion = zeros(size(x))
            wts = zeros(0)

            for j in 1:size(x,1)
               for i in 1:size(x,2)
                  genewisedispersion[j,i] = mle_for_θ([x[j,i]], [estimatedmean[j,i]], wts, rbm, lambda; maxIter=20)
               end
            end

            global genewisedispersion = genewisedispersion
            rbm.inversedispersion .= genewisedispersion

         else
            genewisedispersion = zeros(0)
            wts = zeros(0)
            for i in 1:size(x,2)
               push!(genewisedispersion,mle_for_θ(x[:,i], estimatedmean[:,i], wts, rbm, lambda; maxIter=20))
            end

            geneorder = collect(1:size(genewisedispersion,1))
            genewisedispersion = DataFrame(hcat(genewisedispersion, geneorder))
            genewisedispersion = sort(genewisedispersion, :x1)

            lowestdispersion = Array{Int64,1}(genewisedispersion[findall(genewisedispersion[1] .> quantile(genewisedispersion[1],dispersionquantile)),2])

            globaldispersion = mle_for_θ_global(x[:,lowestdispersion], estimatedmean[:,lowestdispersion], wts, rbm, lambda; maxIter=20)
            rbm.inversedispersion .= globaldispersion
            if epoch % 30 == 0
               println("Dispersion = ", 1 ./ rbm.inversedispersion)
            end
         end

      end
   end

   rbm
end

function libsizenormalization!(x::Array{Float64,2})

   # Library size normalization (see Splatter Package Bioconductor [https://github.com/Oshlack/splatter/blob/382e22e6b22f32f2bb86a4152af369eda19aaba3/R/splat-estimate.R])
   libsize = mapslices(sum, x, dims = 2)
   libmedian = median(libsize)
   x .= (x ./ libsize) .* libmedian
   x
end

function libsizenormalization(x::Array{Float64,2})

   # Library size normalization (see Splatter Package Bioconductor [https://github.com/Oshlack/splatter/blob/382e22e6b22f32f2bb86a4152af369eda19aaba3/R/splat-estimate.R])
   libsize = mapslices(sum, x, dims = 2)
   libmedian = median(libsize)
   normalizedcounts = (x ./ libsize) .* libmedian
   normalizedcounts
end

function normalization!(x::Array{Float64,2})

   # Normalization (Eraslan et al. [https://doi.org/10.1038/s41467-018-07931-2] pp. 11-12)
   libsize = mapslices(sum, x, dims = 2)
   medianlibsize = median(libsize)
   sizefactors = libsize ./ medianlibsize
   sizefactors = Diagonal(sizefactors[:,1])

   x .= log.((inv(sizefactors) * x) .+ 1)
   x
end

function normalization(x::Array{Float64,2})

   # Normalization (Eraslan et al. [https://doi.org/10.1038/s41467-018-07931-2] pp. 11-12)
   libsize = mapslices(sum, x, dims = 2)
   medianlibsize = median(libsize)
   sizefactors = libsize ./ medianlibsize
   sizefactors = Diagonal(sizefactors[:,1])

   normalizedcounts = log.((inv(sizefactors) * x) .+ 1)
   normalizedcounts
end

function logisticfunction(x::Array{Float64,1}, x0::Array{Float64,1}, k::Array{Float64,1})
   1 ./ (1 .+ exp.(-k .* (x .- x0)))
end

function estimatedropout!(x::Matrix{Float64}, rbm::NegativeBinomialBernoulliRBM)

   normalizedcounts = x

   means = mapslices(mean, normalizedcounts, dims =1)
   logcounts = log1p.(means)
   logcounts = logcounts[1,:]
   zeroproportion = mapslices(sum, normalizedcounts .== 0, dims = 1) / size(normalizedcounts,1)
   zeroproportion = zeroproportion[1,:]

   dataset = DataFrame(hcat(zeroproportion,logcounts))
   rename!(dataset, :x1 => :zeroproportion, :x2 => :logcounts)

   tmp = fit(GeneralizedLinearModel, @formula(zeroproportion ~ logcounts), dataset, Binomial());
   rbm.dropoutmid .= coef(tmp)[1]
   rbm.dropoutshape .= coef(tmp)[2]
end

function sampledropout(x::Matrix{Float64}, rbm::NegativeBinomialBernoulliRBM)

   dropoutmask = Array{Float64,2}(undef,size(x))
   for i in 1:size(x,2)

      dropoutmask[:,i] = logisticfunction(x[:,i], rbm.dropoutmid, rbm.dropoutshape)
   end

   bernoulli!(dropoutmask)
   1 .- dropoutmask
end

function sampledropout(x::Matrix{Float64}, dbm::MultimodalDBM)

   @. model(x, p) = 1 / (1 + exp(-p[2] * (x - p[1])))
   dropoutmask = Array{Float64,2}(undef,size(x))
   for i in 1:size(x,2)

      dropoutmask[:,i] = bernoulli!(model(log1p.(x[:,i]),[dbm[1].dropoutmid[i], dbm[1].dropoutshape[i]]))
   end

   1 .- dropoutmask
end

function fittedlogistic(x::Matrix{Float64}, rbm::NegativeBinomialBernoulliRBM)

   @. model(x, p) = 1 / (1 + exp(-p[2] * (x - p[1])))
   logisticestimate = Array{Float64,2}(undef,size(x))
   for i in 1:size(x,2)

      logisticestimate[:,i] = model(log1p.(x[:,i]),[rbm.dropoutmid[i], rbm.dropoutshape[i]])
   end

   logisticestimate
end

"""
    initrbm(x, nhidden)
    initrbm(x, nhidden, rbmtype)
Creates a RBM with `nhidden` hidden units and initalizes its weights for
training on dataset `x`.
`rbmtype` can be a subtype of `AbstractRBM`, default is `BernoulliRBM`.
"""
function initrbm(x::Array{Float64,2}, nhidden::Int,
   inversedispersion::AbstractArray, dropoutmid::Array{Float64,1}, dropoutshape::Array{Float64,1},
      rbmtype::DataType = BernoulliRBM)

   nsamples, nvisible = size(x)
   weights = randn(nvisible, nhidden)/sqrt(nvisible)
   hidbias = zeros(nhidden)

   if rbmtype == BernoulliRBM
      visbias = initvisiblebias(x)
      return BernoulliRBM(weights, visbias, hidbias)

   elseif rbmtype == GaussianBernoulliRBM
      visbias = vec(mean(x, dims = 1))
      sd = vec(std(x, dims = 1))
      GaussianBernoulliRBM(weights, visbias, hidbias, sd)

   elseif rbmtype == GaussianBernoulliRBM2
      visbias = vec(mean(x, dims = 1))
      sd = vec(std(x, dims = 1))
      weights .*= sd
      return GaussianBernoulliRBM2(weights, visbias, hidbias, sd)

   elseif rbmtype == BernoulliGaussianRBM
      visbias = initvisiblebias(x)
      return BernoulliGaussianRBM(weights, visbias, hidbias)

   elseif rbmtype == Binomial2BernoulliRBM
      visbias = initvisiblebias(x/2)
      return Binomial2BernoulliRBM(weights, visbias, hidbias)

   elseif rbmtype == NegativeBinomialBernoulliRBM
      weights .= -rand(Uniform(0.0005,0.002),nvisible,nhidden)
      inversedispersion = inversedispersion
      visbias = initvisiblebiasnegativebinomial(x,inversedispersion)
      dropoutmid = [0.5]
      dropoutshape = [1.0]
      return NegativeBinomialBernoulliRBM(weights, visbias, hidbias, inversedispersion, dropoutmid, dropoutshape)

   else
      error(string("Datatype for RBM is unsupported: ", rbmtype))
   end
end

"""
    initvisiblebias(x)
Returns sensible initial values for the visible bias for training an RBM
on the data set `x`.
"""
function initvisiblebias(x::Array{Float64,2})
   nvisible = size(x,2)
   initbias = zeros(nvisible)
   for j=1:nvisible
      empprob = mean(x[:,j])
      if empprob > 0
         initbias[j] = log(empprob/(1-empprob))
      end
   end
   initbias
end

function initvisiblebiasnegativebinomial(x::Array{Float64,2}, inversedispersion::Array{Float64,1})

   nvisible = size(x,2)
   means = vec(mapslices(mean,x, dims = 1))

   theta = log.(means ./ (means .+ inversedispersion))
   theta
end

"""
    nomonitoring
Accepts a model and a number of epochs and returns nothing.
"""
function nomonitoring(bm, epoch)
end

"""
    randombatchmasks(nsamples, batchsize)
Returns BitArray-Sets for the sample indices when training on a dataset with
`nsamples` samples using minibatches of size `batchsize`.
"""
function randombatchmasks(nsamples::Int, batchsize::Int)
   if batchsize > nsamples
      error("Batchsize ($batchsize) exceeds number of samples ($nsamples).")
   end

   nfullbatches, nremainingbatch = divrem(nsamples, batchsize)
   batchsizes = fill(batchsize, nfullbatches)
   if nremainingbatch > 0
      push!(batchsizes, nremainingbatch)
   end
   randomsampleindices = randperm(nsamples)
   map(batchrange -> randomsampleindices[batchrange],
         ranges(batchsizes))
end

function sdupdateterm(gbrbm::GaussianBernoulliRBM,
         v::Matrix{Float64}, h::Matrix{Float64})
   vec(mean((v .- gbrbm.visbias').^2 ./ (gbrbm.sd' .^ 3) -
         h * gbrbm.weights' .* (v ./ (gbrbm.sd' .^ 2)), dims = 1))
end

function assert_enoughvaluesforepochs(vname::String, v::Vector, epochs::Int)
   if length(v) < epochs
      error("Not enough `$vname` (vector of length $(length(v))) " .*
         "for training epochs ($epochs)")
   end
end

function assertinitoptimizers(optimizer::AbstractOptimizer,
      optimizers::Vector{<:AbstractOptimizer}, bm::BM,
      learningrates::Vector{Float64}, sdlearningrates::Vector{Float64},
      epochs::Int
      ) where {BM<:AbstractBM}

   if isempty(optimizers)
      if optimizer === NoOptimizer()
         # default optimization algorithm
         optimizers = map(
               i -> defaultoptimizer(bm, learningrates[i], sdlearningrates[i]),
               1:epochs)
      else
         optimizers = fill(initialized(optimizer, bm), epochs)
      end
   else
      optimizers = map(opt -> initialized(opt, bm), optimizers)
      assert_enoughvaluesforepochs("optimizers", optimizers, epochs)
   end

   optimizers
end

function defaultoptimizer(rbm::R,
      learningrate::Float64, sdlearningrate::Float64) where {R<: AbstractRBM}

   loglikelihoodoptimizer(rbm; learningrate = learningrate,
         sdlearningrate = sdlearningrate)
end

function defaultoptimizer(dbm::MultimodalDBM, learningrate::Float64,
      sdlearningrate::Float64)

   StackedOptimizer(map(
         rbm -> defaultoptimizer(rbm, learningrate, sdlearningrate),
         dbm))
end

"""
    trainrbm!(rbm, x)
Trains the given `rbm` for one epoch using the data set `x`.
(See also function `fitrbm`.)

# Optional keyword arguments:
* `learningrate`, `cdsteps`, `sdlearningrate`, `upfactor`, `downfactor`,
   `optimizer`:
   See documentation of function `fitrbm`.
* `chainstate`: a matrix for holding the states of the RBM's hidden nodes. If
   it is specified, PCD is used.
"""
function trainrbm!(rbm::AbstractRBM, x::Array{Float64,2};
      chainstate::Matrix{Float64} = Matrix{Float64}(undef, 0, 0),
      upfactor::Float64 = 1.0,
      downfactor::Float64 = 1.0,
      learningrate::Float64 = 0.005,
      cdsteps::Int = 1,
      batchsize::Int = 1,
      sdlearningrate::Float64 = 0.0,
      fisherscoring::Int = 1,
      zeroinflation::Bool = true,
      limit::Int = 10,
      epsilon::Float64 = 0.0001,
      optimizer::AbstractOptimizer = LoglikelihoodOptimizer(rbm;
            learningrate = learningrate, sdlearningrate = sdlearningrate),
      sampler::AbstractSampler = (cdsteps == 1 ? NoSampler() : GibbsSampler(cdsteps - 1)),

      # write-only arguments for reusing allocated space:
      v::Matrix{Float64} = Matrix{Float64}(undef, batchsize, length(rbm.visbias)),
      h::Matrix{Float64} = Matrix{Float64}(undef, batchsize, length(rbm.hidbias)),
      hmodel::Matrix{Float64} = Matrix{Float64}(undef, batchsize, length(rbm.hidbias)),
      vmodel::Matrix{Float64} = Matrix{Float64}(undef, batchsize, length(rbm.visbias)))

   nsamples = size(x, 1)

   # perform PCD if a chain state is provided as parameter
   pcd = !isempty(chainstate)

   batchmasks = randombatchmasks(nsamples, batchsize)
   nbatchmasks = length(batchmasks)

   normalbatchsize = true

   # Zero-Inflation
   if zeroinflation && typeof(rbm) == NegativeBinomialBernoulliRBM

      estimatedropout!(x,rbm)
   end

   global estimatedmean = zeros(0,size(x, 2))
   # go through all samples or thorugh all batches of samples
   for batchindex in eachindex(batchmasks)
      batchmask = batchmasks[batchindex]

      normalbatchsize = (batchindex < nbatchmasks || nsamples % batchsize == 0)

      if normalbatchsize
         v .= view(x, batchmask, :)
      else
         v = x[batchmask, :]
         thisbatchsize = nsamples % batchsize
         h = Matrix{Float64}(undef,thisbatchsize, nhiddennodes(rbm))
         if !pcd
            vmodel = Matrix{Float64}(undef,size(v))
            hmodel = Matrix{Float64}(undef,size(h))
         end # in case of pcd, vmodel and hmodel are not downsized
      end

      # Calculate potential induced by visible nodes, used for update
      hiddenpotential!(h, rbm, v, upfactor)

      # In case of CD, start Gibbs chain with the hidden state induced by the
      # sample. In case of PCD, start Gibbs chain with
      # previous state of the Gibbs chain.
      if pcd
         hmodel = chainstate # note: state of chain will be visible by the caller
      else
         copyto!(hmodel, h)
      end

      samplehiddenpotential!(hmodel, rbm)

      # Additional sampling steps, may be customized
      sample!(vmodel, hmodel, sampler, rbm, upfactor, downfactor)

      # Do not sample in last step to avoid unnecessary sampling noise
      visiblepotential!(vmodel, rbm, hmodel, downfactor)

      estimatedmean = vcat(estimatedmean, vmodel)
      global estimatedmean = estimatedmean

      hiddenpotential!(hmodel, rbm, vmodel, upfactor)

      if !normalbatchsize && pcd
         # remove additional samples
         # that are unnecessary for computing the gradient
         vmodel = vmodel[1:thisbatchsize, :]
         hmodel = hmodel[1:thisbatchsize, :]
      end
      computegradient!(optimizer, v, vmodel, h, hmodel, rbm)
      updateparameters!(rbm, optimizer)
   end

   rbm
end

function linkfunction(v::M, rbm::NegativeBinomialBernoulliRBM,
   ) where{M <: AbstractArray{Float64,1}}

   tmp = zeros(size(v))
   tmp .= log1p.(v ./ (rbm.inversedispersion .+ v))
   tmp
end

function linkfunction(v::M, rbm::NegativeBinomialBernoulliRBM,
   ) where{M <: AbstractArray{Float64,2}}

   tmp = zeros(size(v))
   tmp .= log1p.(v ./ (rbm.inversedispersion .+ v))
   tmp
end

function inverselinkfunction(v::M, rbm::NegativeBinomialBernoulliRBM,
   ) where{M <: AbstractArray{Float64,1}}

   tmp = zeros(size(v))
   tmp = (exp.(v) .* rbm.inversedispersion) ./ (1 .- exp.(v))
   tmp
end

function inverselinkfunction(v::M, rbm::NegativeBinomialBernoulliRBM,
   ) where{M <: AbstractArray{Float64,2}}

   v = (exp.(v) .* rbm.inversedispersion') ./ (1 .- exp.(v))
   v
end

function mu(rbm::NegativeBinomialBernoulliRBM)
   tmp = (exp.(rbm.visbias) .* rbm.inversedispersion) ./ (1 .- exp.(rbm.visbias))
   tmp
end

function mu(dbm::MultimodalDBM)
   tmp = (exp.(dbm[1].visbias) .* dbm[1].inversedispersion) ./ (1 .- exp.(dbm[1].visbias))
   tmp
end

function nbvariance(rbm::NegativeBinomialBernoulliRBM)
   nbmean = mu(rbm)
   nbvariance = nbmean .+ (nbmean.^2 ./ rbm.inversedispersion)

   nbvariance
end

function nbvariance(dbm::MultimodalDBM)
   nbmean = mu(dbm)
   nbvariance = nbmean .+ (nbmean.^2 ./ dbm[1].inversedispersion)

   nbvariance
end

function successprob(rbm::NegativeBinomialBernoulliRBM)
   tmp = mu(rbm) ./ (mu(rbm) .+ rbm.inversedispersion')
   tmp
end

function successprob(dbm::MultimodalDBM)
   tmp = dbm[1].inversedispersion' ./ (mu(dbm) .+ dbm[1].inversedispersion')
   tmp
end

function negativebinomialparameters(rbm::NegativeBinomialBernoulliRBM)
   nbmean,successprobability = mu(rbm),successprob(rbm)
   tmp = hcat(nbmean,successprobability)
   estimtatedparameters = DataFrame(tmp)
   rename!(estimtatedparameters, :x1 => :Mu, :x2 => :SuccessProb)
   estimtatedparameters
end

function negativebinomialparameters(dbm::MultimodalDBM)
   nbmean,successprobability = mu(dbm),successprob(dbm)
   tmp = hcat(nbmean,successprobability)
   estimtatedparameters = DataFrame(tmp)
   rename!(estimtatedparameters, :x1 => :Mu, :x2 => :SuccessProb)
   estimtatedparameters
end

function mle_for_θ(y::AbstractVector, μ::AbstractVector, wts::AbstractVector, rbm::NegativeBinomialBernoulliRBM, lambda::Float64;
   maxIter=3, convTol=1.e-6)

   function first_derivative(θ::Real)
      tmp(yi, μi) = -((yi+θ)/(μi+θ) + log(μi+θ) - 1 - log(θ) - digamma(θ+yi) + digamma(θ))
      unit_weights ? sum(tmp(yi, μi) for (yi, μi) in zip(y, μ)) :
      sum(wti * tmp(yi, μi) for (wti, yi, μi) in zip(wts, y, μ))
   end

   function second_derivative(θ::Real)
      tmp(yi, μi) = -(yi+θ)/(μi+θ)^2 + 2/(μi+θ) - 1/θ - trigamma(θ+yi) + trigamma(θ)
      unit_weights ? sum(tmp(yi, μi) for (yi, μi) in zip(y, μ)) :
      sum(wti * tmp(yi, μi) for (wti, yi, μi) in zip(wts, y, μ))
   end

   unit_weights = length(wts) == 0
   if unit_weights
      n = length(y)
      θ = n / sum((yi/μi - 1)^2 for (yi, μi) in zip(y, μ))
   else
      n = sum(wts)
      θ = n / sum(wti * (yi/μi - 1)^2 for (wti, yi, μi) in zip(wts, y, μ))
   end
   δ, converged = one(θ), false

   for t = 1:maxIter
      θ = max(θ,0.01)

      # Add regularization term lambda*(1/theta^2) to the loglikelihood
      # of Fisher scoring algorithm. In this case lambda = 10 and the
      # first derivative is lambda*2/theta^3. The second derivative is
      # lambda*6/theta^4.
      score = first_derivative(θ)
      fisher = second_derivative(θ)
      δ = (score + lambda*2/θ^3) / (fisher + lambda*6/θ^4)
      if abs(δ) <= convTol
         converged = true
      break
   end
   θ = θ + δ

   end
   converged
   θ
end

function mle_for_θ_global(y::Array{Float64,2}, μ::Array{Float64,2}, wts::AbstractVector, rbm::NegativeBinomialBernoulliRBM;
   maxIter=10, convTol=1.e-16)

   function first_derivative(θ::Real)
      tmp(yi, μi) = (yi+θ)/(μi+θ) + log(μi+θ) - 1 - log(θ) - digamma(θ+yi) + digamma(θ)
      unit_weights ? sum(tmp(yi, μi) for (yi, μi) in zip(y, μ)) :
      sum(wti * tmp(yi, μi) for (wti, yi, μi) in zip(wts, y, μ))
   end

   function second_derivative(θ::Real)
      tmp(yi, μi) = -(yi+θ)/(μi+θ)^2 + 2/(μi+θ) - 1/θ - trigamma(θ+yi) + trigamma(θ)
      unit_weights ? sum(tmp(yi, μi) for (yi, μi) in zip(y, μ)) :
      sum(wti * tmp(yi, μi) for (wti, yi, μi) in zip(wts, y, μ))
   end

   unit_weights = length(wts) == 0
   if unit_weights
      n = length(y)
      θ = n / sum((yi/μi - 1)^2 for (yi, μi) in zip(y, μ))
   else
      n = sum(wts)
      θ = n / sum(wti * (yi/μi - 1)^2 for (wti, yi, μi) in zip(wts, y, μ))
   end
   δ, converged = one(θ), false

   for t = 1:maxIter
      θ = abs(θ)
      δ = first_derivative(θ) / second_derivative(θ)

      if abs(δ) <= convTol
         converged = true
      break
   end
   θ = θ - δ

   end
   converged
   θ
end

function mle_for_θ_start(y::AbstractVector, μ::AbstractVector, wts::AbstractVector, rbm::NegativeBinomialBernoulliRBM;
   maxIter=30, convTol=1.e-6)

   function first_derivative(θ::Real)
      tmp(yi, μi) = (yi+θ)/(μi+θ) + log(μi+θ) - 1 - log(θ) - digamma(θ+yi) + digamma(θ)
      unit_weights ? sum(tmp(yi, μi) for (yi, μi) in zip(y, μ)) :
      sum(wti * tmp(yi, μi) for (wti, yi, μi) in zip(wts, y, μ))
   end

   function second_derivative(θ::Real)
      tmp(yi, μi) = -(yi+θ)/(μi+θ)^2 + 2/(μi+θ) - 1/θ - trigamma(θ+yi) + trigamma(θ)
      unit_weights ? sum(tmp(yi, μi) for (yi, μi) in zip(y, μ)) :
      sum(wti * tmp(yi, μi) for (wti, yi, μi) in zip(wts, y, μ))
   end

   unit_weights = length(wts) == 0
   if unit_weights
      n = length(y)

      θ = n / sum((yi/μi - 1)^2 for (yi, μi) in zip(y, μ))
   else
      n = sum(wts)
      θ = n / sum(wti * (yi/μi - 1)^2 for (wti, yi, μi) in zip(wts, y, μ))
   end
   δ, converged = one(θ), false

   for t = 1:maxIter
      θ = abs(θ)
      δ = first_derivative(θ) / second_derivative(θ)
      if abs(δ) <= convTol
         converged = true
      break
   end
   θ = θ - δ
   end
   converged || throw(ConvergenceException(maxIter))
   θ
end

"""
    ConvergenceException(iters::Int, lastchange::Real=NaN, tol::Real=NaN)
The fitting procedure failed to converge in `iters` number of iterations,
i.e. the `lastchange` between the cost of the final and penultimate iteration was greater than
specified tolerance `tol`.
"""
struct ConvergenceException{T<:Real} <: Exception
    iters::Int
    lastchange::T
    tol::T
    function ConvergenceException{T}(iters, lastchange::T, tol::T) where T<:Real
        if tol > lastchange
            throw(ArgumentError("Change must be greater than tol."))
        else
            new(iters, lastchange, tol)
        end
    end
end

ConvergenceException(iters, lastchange::T=NaN, tol::T=NaN) where {T<:Real} =
    ConvergenceException{T}(iters, lastchange, tol)

function Base.showerror(io::IO, ce::ConvergenceException)
    print(io, "failure to converge after $(ce.iters) iterations.")
    if !isnan(ce.lastchange)
        print(io, " Last change ($(ce.lastchange)) was greater than tolerance ($(ce.tol)).")
    end
end

"""Standardizes a data matrix x """
function standardize!(x;standardize::Bool=true)
   xmean = mean(x,dims=1)
   xstd = std(x;dims=1)

   for j in 1:size(x)[2]
      for i in 1:size(x)[1]
         if standardize
            x[i,j] = (x[i,j] - xmean[j])/xstd[j]
         else
            x[i,j] = (x[i,j] - xmean[j])
         end
      end
   end
   nothing
end

""" Calculates univariate betas, given an explanatory matrix x and a response y """
function calcunibeta(x,y)
   n, p = size(x)
   unibeta = zeros(p)

   for i=1:p
      denom = 0.0
       for j=1:n
         unibeta[i] += x[j,i]*y[j]
         denom += x[j,i]*x[j,i]
      end
      unibeta[i] /= denom
   end
   unibeta
end

""" Calculates univariate betas, between a selected variable (candidate) in x and all other variables in x. Betas that have been calculated in advance are stored in the matrix covcache."""
function calcunibeta(x,candidate,covcache)
   n, p = size(x)
   unibeta = zeros(p-1)

   for j=1:(p-1)
      targetj = j < candidate ? j : j+1

      if isnan(covcache[candidate,targetj])
          unibeta[j] = sum(x[:,targetj].*x[:,candidate])
          denom = sum(x[:,targetj].*x[:,targetj])

         covcache[candidate,targetj] = covcache[targetj,candidate] = unibeta[j]
      else
          unibeta[j] = covcache[candidate,targetj]
          denom = sum(x[:,targetj].*x[:,targetj])
      end
      unibeta[j] /= denom
   end
   unibeta
end

""" For each variable "candidate" in inmat, the variables of inmat are selected that best explain the state of "candidate", given stepno steps."""
function partboost(inmat;
                   stepno::Int=10,nu::Float64=0.1,csf::Float64=0.9)
   n, inp = size(inmat)
   p = inp - 1
   betamat = zeros(Float64,inp,inp)

   covcache = fill(NaN,inp,inp)

   indenom = var(inmat;dims=1).*(n-1)

   for candidate=1:inp

       y = view(inmat,:,candidate)
       x = view(inmat,:,[1:(candidate-1); (candidate+1):inp])
      denom = view(indenom,[1:(candidate-1); (candidate+1):inp])

      nuvec = fill(nu,p)
      penvec = fill(n*(1/nu-1),p)

      actualnom = calcunibeta(inmat,candidate,covcache)
      beta = view(betamat,candidate,[1:(candidate-1); (candidate+1):inp])

      actualupdate = 0.0
      actualsel = -1

      for step in 1:stepno
          if step > 1
            for j=1:p
               inj = j < candidate ? j : j+1
               insel = actualsel < candidate ? actualsel : actualsel+1

               if isnan(covcache[inj,insel])
                   covval = sum(x[:,j] .* x[:,actualsel])
                   covcache[inj,insel] = covcache[insel,inj] = covval
               else
                  covval = covcache[inj,insel]
               end
               actualnom[j] -= actualupdate*covval/denom[j]
            end
         end

          actualsel = argmax(((actualnom.*denom) ./ (denom.+penvec)).^2)
         actualupdate = nuvec[actualsel]*actualnom[actualsel]
         beta[actualsel] += actualupdate
         nuvec[actualsel] = 1-(1-nuvec[actualsel])^csf
         penvec[actualsel] = n*(1/nuvec[actualsel]-1)

      end
   end

   betamat
end

""" Aggregates distances in distmat using aggfun."""
function clusterdist(distmat,indices1,indices2,aggfun=mean)
    aggfun(view(distmat,indices1,indices2))
end

""" Hierarchically clusters variables based on their pairwise dissimilarity stored in distmat. Minimum and maximum cluster size is given by minsize and maxsize respectively. By default, clustering is performed based on average linkage (aggfun=mean). Returns a dictionary of clusters, where each cluster is represented by a vector containing the cluster mebmbers."""
function clusterbeta(distmat,aggfun=mean;minsize::Int=2 ,maxsize::Int=size(distmat)[2])
   p = size(distmat)[2]
   maxdist = maximum(distmat) + 1.0

   cldist = copy(distmat)
   cldict = Dict{Int,Vector{Int}}()
   for i=1:p
      cldict[i] = [i]
   end

   global test = cldict
   while true
      cllen = map(key -> length(cldict[key]),collect(keys(cldict))) # added collect() function for julia 1.0.3

      global test = cllen
      if minimum(cllen) >= minsize
         break
      end

       _keys = sort(collect(keys(cldict)))
       source, target, actualmin = findmindist(cldist,_keys)

      if (length(cldict[source]) + length(cldict[target]) > maxsize &&

          length(cldict[source]) >= minsize &&
          length(cldict[target]) >= minsize)

          cldist[source,target] = cldist[target,source] = maxdist
         continue
      end

      append!(cldict[target],cldict[source])
      delete!(cldict,source)
      for key in keys(cldict)
         if key != target
            cldist[target,key] = cldist[key,target] =
               clusterdist(distmat,cldict[target],cldict[key],aggfun)
         end
      end
   end

    cldict
end

"""Makes an asymmetric similarity matrix symmetric. The distance (D) between i and j or j and i is replaced by max(|D_ij|,|D_ji|). """
function calcbetawish(betamat)
   p = size(betamat)[2]
   wishmat = zeros(p,p)
   for i=1:(p-1)
      for j=(i+1):p
         wishmat[i,j] = wishmat[j,i] = max(abs(betamat[i,j]),abs(betamat[j,i]))
      end
   end
   wishmat
end

""" Identifies the variables besti and bestjj with the lowest distance curmin given in the distance matrix distmat."""
function findmindist(distmat,indices=1:size(distmat,2))

   curmin = -1
   besti = bestj = -1

   for i=1:(length(indices)-1)
      for j=(i+1):length(indices)
         if (i == 1 && j == 2) || distmat[indices[i],indices[j]] < curmin
            besti = indices[i]
            bestj = indices[j]
            curmin = distmat[indices[i],indices[j]]
         end
      end
   end

   besti, bestj, curmin
end

""" Identifies clusters of related SNPs using stagewise regression. SNP clusters are allowed to have a maximum of maxsize SNPs while at least are required to have a minimum of minsize SNPs. stepno indicates the number of steps in stagewise regression."""
function findpartition(mat;minsize::Int=5,maxsize::Int=20,stepno::Int=100)
    dat = deepcopy(mat)
    standardize!(dat,standardize=true)
    betamat = partboost(dat,stepno=stepno)
    betawish = calcbetawish(betamat)
    cldict = clusterbeta(-betawish,minsize=minsize,maxsize=maxsize)
    mvisibleindex = Array{Array{Int,1}}(undef,length(cldict))
    clkeys = collect(keys(cldict))
    for i=1:length(clkeys)
        mvisibleindex[i] = cldict[clkeys[i]]
    end
    mvisibleindex
end

""" Identifies adequate size of hidden units for each separate DBM"""
function partitionhiddens(partitions,nhiddens)
    nparts = length(partitions)
    p= length(cat(partitions..., dims=1))
    nhiddensmat = zeros(Int,nparts,length(nhiddens))
    for i=1:(nparts-1)
        nhiddensmat[i,:] = floor.(nhiddens ./ (p/length(partitions[i])))
    end
    nhiddensmat[nparts,:] = nhiddens .- vec(sum(nhiddensmat,dims=1))
    nhiddensmat
end

""" Given the partitions, identified by findpartition, a separate DBM is trained on each partition. The overall DBM architecture is given by nhiddens. The units in nhiddens are distributed across the partitions, such that each partition has at least p units in hidden layer 1 and at least a single unit in hidden layer 2."""
function trainpartdbm(mat;partitions::Array{Array{Int64,1},1}=[collect(1:size(mat)[2])],
                      nhiddens::Array{Int64,1}=[size(mat)[2],div(size(mat)[2],10)],
                      epochs::Int64=20,
                      nparticles::Int64=100,
                      learningrate::Float64=0.005)

    partDBMS = Vector{BasicDBM}()
    nhiddensmat = partitionhiddens(partitions,nhiddens)
    println(nhiddensmat)
     for i=1:length(partitions)
        partx = mat[:,partitions[i]]
        pdbm = fitdbm(partx,
                                        nhiddens=nhiddensmat[i,:],
                                        epochs =epochs,
                                        nparticles=nparticles,
                                        learningrate=learningrate
                                        )
        push!(partDBMS,pdbm)
    end
    jointdbm = joindbms(partDBMS,partitions)
    jointdbm,partDBMS
end


text: # Julia code goes in this folder

max_stars_repo_name: TheProtaganist/CodeJ

max_issues_repo_name: TheProtaganist/CodeJ

max_forks_repo_name: TheProtaganist/CodeJ


text:
include("./types.jl")
include("./tree.jl")
include("./boost.jl")
# using Juno
# using Gallium

# Train a deepboost model on the given examples, using
# numIter iterations (which not necessarily means numIter trees)
function Train(train_examples::Vector{Example}, model::Model, tree_depth::Int64,
        num_iter::Int64, beta::Float64, lambda::Float64, loss_type::String, verbose::Bool)
# model: Model, need return
# train_examples: pointer, since it is a Vector, change is in-place, no need to return
	# Train the model
	for iter in 1:num_iter
		(train_examples, model) = AddTreeToModel(train_examples, model, loss_type,
                                             beta, lambda, tree_depth)
    # println("successfully run for ", iter, "iterations.")
    #Juno.@step gcd(12, 4)
    # @enter gcd(5, 20)
    println(model[1].weight)
		if verbose
   		(error, avg_tree_size, num_trees) = EvaluateModel(train_examples, model)
			println("Iteration: ", iter, ", error: ", error,
			        ", avg tree size: ", avg_tree_size,
			        ", num trees: ", num_trees, "\n")

		end
	end

  return model
end

# Classify examples using model
function Predict(examples::Vector{Example}, model::Model)
	# TODO::initiate labels
	labels  = zeros(Int64, length(examples)) # Vector{Label::Int64}
	for i in 1:length(examples)
		labels[i] = ClassifyExampleModel(examples[i], model)
  end
	return labels
end

# Compute the error of model on examples. Also compute the number of trees in
# model and their average size.
function Evaluate(examples::Vector{Example}, model::Model)
# error:
# avg_tree_size
# num_trees
  (error, avg_tree_size, num_trees) = EvaluateModel(examples, model)
  return (error, avg_tree_size, num_trees)
end


text:
## ExoplanetsSysSim/src/eval_model.jl
## (c) 2015 Eric B. Ford

# This merely demonstrates how to setup your own function to evaluate a simulated model against an observed catalog
# This code need not be used for any actual calculations
module TestEvalModel
export test_eval_model, evaluate_model_scalar_ret
using ExoplanetsSysSim

sim_param_closure = SimParam()
cat_phys_try_closure = KeplerPhysicalCatalog([])
cat_obs_try_closure = KeplerObsCatalog([])
summary_stat_try_closure =  CatalogSummaryStatistics()
summary_stat_ref_closure =  CatalogSummaryStatistics()

include("eval_model_macro.jl")
@make_evaluate_model(sim_param_closure,cat_phys_try_closure,cat_obs_try_closure,summary_stat_try_closure,summary_stat_ref_closure)

function test_eval_model()
  global sim_param_closure = setup_sim_param_demo()
  cat_phys = generate_kepler_physical_catalog(sim_param_closure)
  cat_obs = observe_kepler_targets_single_obs(cat_phys,sim_param_closure)
  global summary_stat_ref_closure = calc_summary_stats_obs_demo(cat_obs,sim_param_closure)
  global cat_phys_try_closure  = generate_kepler_physical_catalog(sim_param_closure)
  global cat_obs_try_closure  = observe_kepler_targets_single_obs(cat_phys_try_closure,sim_param_closure)
  global summary_stat_try_closure  = calc_summary_stats_sim_pass_one_demo(cat_obs_try_closure,cat_phys_try_closure,sim_param_closure)
  summary_stat_try_closure   = calc_summary_stats_sim_pass_two_demo(cat_obs_try_closure,cat_phys_try_closure,summary_stat_try_closure,sim_param_closure)
  param_guess = make_vector_of_sim_param(sim_param_closure)
  evaluate_model_scalar_ret( param_guess)
end

function test_eval_model_vs_sim_data()
  global sim_param_closure = setup_sim_param_demo()
  cat_phys = generate_kepler_physical_catalog(sim_param_closure)
  # cat_obs = observe_kepler_targets_single_obs(cat_phys,sim_param_closure)
  cat_obs = simulated_read_kepler_observations(sim_param_closure)
  global summary_stat_ref_closure = calc_summary_stats_obs_demo(cat_obs,sim_param_closure)
  global cat_phys_try_closure  = generate_kepler_physical_catalog(sim_param_closure)
  global cat_obs_try_closure  = observe_kepler_targets_single_obs(cat_phys_try_closure,sim_param_closure)
  global summary_stat_try_closure  = calc_summary_stats_sim_pass_one_demo(cat_obs_try_closure,cat_phys_try_closure,sim_param_closure)
  summary_stat_try_closure   = calc_summary_stats_sim_pass_two_demo(cat_obs_try_closure,cat_phys_try_closure,summary_stat_try_closure,sim_param_closure)
  param_guess = make_vector_of_sim_param(sim_param_closure)
  evaluate_model_scalar_ret( param_guess)
end

end


text:
"""
$(TYPEDEF)

An equality relationship between two expressions.

# Fields
$(FIELDS)
"""
struct Equation
    """The expression on the left-hand side of the equation."""
    lhs
    """The expression on the right-hand side of the equation."""
    rhs
end
Base.:(==)(a::Equation, b::Equation) = all(isequal.((a.lhs, a.rhs), (b.lhs, b.rhs)))
Base.hash(a::Equation, salt::UInt) = hash(a.lhs, hash(a.rhs, salt))

Base.show(io::IO, eq::Equation) = print(io, eq.lhs, " ~ ", eq.rhs)

SymbolicUtils.simplify(x::Equation; kw...) = simplify(x.lhs; kw...) ~ simplify(x.rhs; kw...)

lhss(xs) = map(x->x.lhs, xs)
rhss(xs) = map(x->x.rhs, xs)

"""
$(TYPEDSIGNATURES)

Create an [`Equation`](@ref) out of two [`Num`](@ref) instances, or an
`Num` and a `Number`.

# Examples

```jldoctest
julia> using Symbolics

julia> @variables x y;

julia> @variables A[1:3, 1:3] B[1:3, 1:3];

julia> x ~ y
x ~ y

julia> x - y ~ 0
x - y ~ 0

julia> A .~ B
3×3 Array{Equation,2}:
 A₁ˏ₁ ~ B₁ˏ₁  A₁ˏ₂ ~ B₁ˏ₂  A₁ˏ₃ ~ B₁ˏ₃
 A₂ˏ₁ ~ B₂ˏ₁  A₂ˏ₂ ~ B₂ˏ₂  A₂ˏ₃ ~ B₂ˏ₃
 A₃ˏ₁ ~ B₃ˏ₁  A₃ˏ₂ ~ B₃ˏ₂  A₃ˏ₃ ~ B₃ˏ₃

julia> A .~ 3x
3×3 Array{Equation,2}:
 A₁ˏ₁ ~ 3x  A₁ˏ₂ ~ 3x  A₁ˏ₃ ~ 3x
 A₂ˏ₁ ~ 3x  A₂ˏ₂ ~ 3x  A₂ˏ₃ ~ 3x
 A₃ˏ₁ ~ 3x  A₃ˏ₂ ~ 3x  A₃ˏ₃ ~ 3x
```
"""
Base.:~(lhs::Num, rhs::Num) = Equation(value(lhs), value(rhs))
Base.:~(lhs::Num, rhs::Number    ) = Equation(value(lhs), value(rhs))
Base.:~(lhs::Number    , rhs::Num) = Equation(value(lhs), value(rhs))
Base.:~(lhs::Symbolic, rhs::Symbolic) = Equation(value(lhs), value(rhs))
Base.:~(lhs::Symbolic, rhs::Any    ) = Equation(value(lhs), value(rhs))
Base.:~(lhs::Any, rhs::Symbolic    ) = Equation(value(lhs), value(rhs))

struct ConstrainedEquation
  constraints
  eq
end

function expand_derivatives(eq::Equation, simplify=false)
    return Equation(expand_derivatives(eq.lhs, simplify), expand_derivatives(eq.rhs, simplify))
end


text:
# ---
# title: 1300. Sum of Mutated Array Closest to Target
# id: problem1300
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Binary Search
# link: <https://leetcode.com/problems/sum-of-mutated-array-closest-to-target/description/>
# hidden: true
# ---
#
# Given an integer array `arr` and a target value `target`, return the integer
# `value` such that when we change all the integers larger than `value` in the
# given array to be equal to `value`, the sum of the array gets as close as
# possible (in absolute difference) to `target`.
#
# In case of a tie, return the minimum such integer.
#
# Notice that the answer is not neccesarilly a number from `arr`.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [4,9,3], target = 10
#     Output: 3
#     Explanation: When using 3 arr converts to [3, 3, 3] which sums 9 and that's the optimal answer.
#
#
# **Example 2:**
#
#
#
#     Input: arr = [2,3,5], target = 10
#     Output: 5
#
#
# **Example 3:**
#
#
#
#     Input: arr = [60864,25176,27249,21296,20204], target = 56803
#     Output: 11361
#
#
#
#
# **Constraints:**
#
#   * `1 <= arr.length <= 10^4`
#   * `1 <= arr[i], target <= 10^5`
#
#
## @lc code=start
using LeetCode

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


text:
"""
Arborist.jl
🆒🌳📦
An arborist performs surgery on tree.
"""
module Arborist
using CodeTracking
using Revise: Revise  # Revise must be loaded for CodeTracking to work
using MacroTools

include("hack_revise.jl")

include("reflection.jl")
include("ast_tools.jl")

include("grafters.jl")
# must be last as generate function can't call things defined afteer them
include("core.jl")
end # module


text:
"""
    WebsocketServer([; options...])

Constructs a new WebsocketServer, overriding [`serverConfig`](@ref) with the passed options.

# Example
```julia
using Websocket
server = WebsocketServer([; options...])
```
"""
struct WebsocketServer
    config::NamedTuple
    callbacks::Dict{Symbol, Union{Bool, Function}}
    flags::Dict{Symbol, Bool}
    server::Dict{Symbol, Union{Sockets.TCPServer, Array{WebsocketConnection, 1}, Nothing}}

    function WebsocketServer(; config...)
        @debug "WebsocketClient"
        config = merge(serverConfig, (; config...))
        config = merge(config, (; maskOutgoingPackets = false, type = "server",))
        self = new(
            config,
            Dict{Symbol, Union{Bool, Function}}(
                :client => false,
                :connectError => false,
                :closed => false,
                :listening => false
            ),
            Dict{Symbol, Bool}(
                :isopen => false
            ),
            Dict{Symbol, Union{Sockets.TCPServer, Array{WebsocketConnection, 1}, Nothing}}(
                :clients => Array{WebsocketConnection, 1}(),
                :socket => nothing
            )
        )
    end

end
"""
    listen(callback::Function, server::Websocket.WebsocketServer, event::Symbol)
Register event callbacks onto a server. The callback must be a function with exactly one argument.

Valid events are:
- :listening
- :client
- :connectError
- :closed

# Example
```julia
listen(server, :client) do client
    #...
end
```
"""
function listen(
    cb::Function,
    self::WebsocketServer,
    event::Symbol
)
    if !haskey(self.callbacks, event)
        return @warn "WebsocketServer has no listener for :$event."
    end
    if haskey(self.callbacks, event) && !(self.callbacks[event] isa Function)
        self.callbacks[event] = data -> (
            @async try
                cb(data)
            catch err
                err = CallbackError(err, catch_backtrace())
                err.log()
                exit()
            end
        )
    end
end

function validateUpgrade(headers::HTTP.Messages.Request)
    if !HTTP.hasheader(headers, "Upgrade", "websocket")
        throw(error("""did not receive "Upgrade: websocket" """))
    end
    if !HTTP.hasheader(headers, "Connection", "Upgrade") || HTTP.hasheader(headers, "Connection", "keep-alive upgrade")
        throw(error("""did not receive "Connection: Upgrade" """))
    end
    if !HTTP.hasheader(headers, "Sec-WebSocket-Version", "13") && !HTTP.hasheader(headers, "Sec-WebSocket-Version", "8")
        throw(error("""did not receive "Sec-WebSocket-Version: [13 or 8]" """))
    end
    if !HTTP.hasheader(headers, "Sec-WebSocket-Key")
        throw(error("""did not receive "Sec-WebSocket-Key" header."""))
    end
end

"""
    serve(server::Websocket.WebsocketServer, [port::Int, host::String; options...])
Opens up a TCP connection listener. Blocks while the server is listening.

Defaults:
- port: 8080
- host: "localhost"

`; options...` are passed to the underlying [HTTP.servers.listen](https://juliaweb.github.io/HTTP.jl/stable/public_interface/#Server-/-Handlers-1)

# Example
```julia
closed = Condition()
#...
@async serve(server, 8081, "localhost")
#...
wait(closed)
```
"""
function serve(self::WebsocketServer, port::Int = 8080, host::String = "localhost"; options...)
    @debug "WebsocketServer.listen"
    config = self.config
    options = merge(serverOptions, (; options...))

    try
        host = getaddrinfo(host)
        if config.ssl
            tlsconfig = HTTP.Servers.SSLConfig(config.sslcert, config.sslkey)
            options = merge(options, (; sslconfig = tlsconfig))
        end
        callback = self.callbacks[:client]
        callback === false && throw(error("tried to bind the server before registering \":client\" handler"))
        self.server[:server] = Sockets.listen(host, port)

        VERSION >= v"1.3" && Sockets.nagle(self.server[:server], config.useNagleAlgorithm) #Sockets.nagle needs Julia >= 1.3

        self.flags[:isopen] = true
        if self.callbacks[:listening] isa Function
            self.callbacks[:listening]((; port = port, host = host))
        end
        HTTP.listen(; server = self.server[:server],  options...) do io
            try
                headers = io.message
                validateUpgrade(headers)
                HTTP.setstatus(io, 101)
                key = string(HTTP.header(headers, "Sec-WebSocket-Key"))
                HTTP.setheader(io, "Sec-WebSocket-Accept" => acceptHash(key))
                HTTP.setheader(io, "Upgrade" => "websocket")
                HTTP.setheader(io, "Connection" => "Upgrade")

                startwrite(io)

                client = WebsocketConnection(io.stream, config, self.server[:clients])
                push!(self.server[:clients], client)
                callback(client)
                if HTTP.hasheader(headers, "Sec-WebSocket-Extensions")
                    close(client, CLOSE_REASON_EXTENSION_REQUIRED)
                else
                    startConnection(client, io)
                end
            catch err
                @error err exception = (err, catch_backtrace())
                HTTP.setstatus(io, 400)
                startwrite(io)
            end
        end
    catch err
        self.flags[:isopen] = false
        if typeof(err) === Base.IOError && occursin("software caused connection abort", err.msg)
            callback = self.callbacks[:closed]
            if callback isa Function
                callback((; host = host, port = port))
            else
                @info "The websocket server was closed cleanly:" host = host port = port
            end
            return
        end
        err = ConnectError(err, catch_backtrace())
        callback = self.callbacks[:connectError]
        if callback isa Function
            callback(err)
        else
            err.log()
        end
    end
end

"""
    emit(server::Websocket.WebsocketServer, data::Union{Array{UInt8,1}, String, Number})
Sends the given `data` as a message to all clients connected to the server.

# Example
```julia
emit(server, "Hello everybody from your loving server.")
```
"""
function emit(self::WebsocketServer, data::Union{Array{UInt8,1}, String, Number})
    for client in self.server[:clients]
        send(client, data)
    end
end
"""
    close(server::Websocket.WebsocketServer)
Gracefully disconnects all connected clients, then closes the TCP socket listener.
"""
function Base.close(self::WebsocketServer)
    @sync begin
        for client in self.server[:clients]
            close(client, CLOSE_REASON_GOING_AWAY)
            @async wait(client.closed)
        end
    end
    close(self.server[:server])
    while isopen(self.server[:server]) || isopen(self)
        sleep(0.1)
    end
end
"""
    isopen(server::Websocket.WebsocketServer)::Bool
Returns a Bool indication if the server TCP listener is open.
"""
function Base.isopen(server::WebsocketServer)
    server.flags[:isopen]
end
"""
    length(server::WebsocketServer)::Int
Returns the number of clients connected to the server.
"""
function Base.length(self::WebsocketServer)
    length(self.server[:clients])
end

max_stars_repo_name: citkane/Websocket

max_issues_repo_name: citkane/Websocket

max_forks_repo_name: citkane/Websocket


text:
using Crayons

function Base.show(io::IO, ::MIME"text/markdown", env::AbstractGridWorld)
    p, d = get_agent_pos(env), get_agent_dir(env)
    w = convert(GridWorldBase, env)

    println(io, "World:")
    for i in 1:size(w, 2)
        for j in 1:size(w, 3)
            if CartesianIndex(i, j) ∈ get_agent_view_inds(env)
                bg = :dark_gray
            else
                bg = :black
            end
            if i == p[1] && j == p[2]
                agent = get_agent(env)
                print(io, Crayon(background=bg, foreground=get_color(agent),reset=true), convert(Char, agent))
            else
                o = get_object(env)[findfirst(w.grid[:, i, j])]
                print(io, Crayon(background=bg, foreground=get_color(o),reset=true),  convert(Char, o))
            end
        end
        println(io, Crayon(reset=true))
    end
    println(io)

    println(io, "Agent's view:")
    v = get_agent_view(env)
    for i in 1:size(v, 2)
        for j in 1:size(v, 3)
            if i == 1 && j == size(v, 3) ÷ 2 + 1
                print(io, Agent(dir=DOWN))
            else
                x = findfirst(v[:, i, j])
                if isnothing(x)
                    print(io, '_')
                else
                    print(io, get_object(env)[x])
                end
            end
        end
        println(io)
    end
    println(io)
end


text:
using BinaryBuilder

# Collection of sources required to build Nettle
name = "ntl"
version = v"11.4.3"
sources = [
    ArchiveSource("https://www.shoup.net/ntl/ntl-$(version).tar.gz",
                  "b7c1ccdc64840e6a24351eb4a1e68887d29974f03073a1941c906562c0b83ad2"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/ntl-*/src
./configure PREFIX="${prefix}" GMP_PREFIX="${prefix}" NATIVE=off SHARED=on

make -j${nproc}
make install

install_license ../doc/copying.txt
"""

# Bootstrapping problem; only do natively-runnable platforms
platforms = [
    Linux(:x86_64),
    Linux(:i686),
    Linux(:x86_64; libc=:musl),
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libntl", :libntl),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("GMP_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
;; ---------------------------------------------------------------------
;;	text_create.jl
;;
;;					Jun/06/2013
;;
;; --------------------------------------------------------------------
(require 'tables)
(load-file "/var/www/data_base/common/rep_common/text_manipulate.jl")
;; --------------------------------------------------------------------

;; --------------------------------------------------------------------
(defun data_prepare_proc ()
	(setq dict_aa (make-table string-hash string=))
;
	(dict_append_proc dict_aa "t2381" "名古屋" "21386" "1930-2-14")
	(dict_append_proc dict_aa "t2382" "豊橋" "69724" "1930-10-21")
	(dict_append_proc dict_aa "t2383" "岡崎" "34513" "1930-5-16")
	(dict_append_proc dict_aa "t2384" "一宮" "52371" "1930-9-23")
	(dict_append_proc dict_aa "t2385" "蒲郡" "41852" "1930-8-9")
	(dict_append_proc dict_aa "t2386" "常滑" "82573" "1930-5-19")
	(dict_append_proc dict_aa "t2387" "大府" "51248" "1930-9-12")
	(dict_append_proc dict_aa "t2388" "瀬戸" "27453" "1930-1-18")
	(dict_append_proc dict_aa "t2389" "犬山" "51749" "1930-2-11")
)

;; --------------------------------------------------------------------
(format t "*** 開始 ***")
(setq file_out (cadr command-line-args))
(format t file_out)

(setq dict_aa (data_prepare_proc))
(table-walk dict_display_proc_single dict_aa)

(text_write_proc file_out dict_aa)

(format t "*** 終了 ***")

;; ---------------------------------------------------------------------

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
immutable DiffCache{T, S}
    du::Vector{T}
    dual_du::Vector{S}
end

Base.@pure function DiffCache{chunk_size}(T, length, ::Type{Val{chunk_size}})
    DiffCache(zeros(T, length), zeros(Dual{chunk_size, T}, length))
end

Base.@pure DiffCache(u::AbstractArray) = DiffCache(eltype(u),length(u),Val{ForwardDiff.pickchunksize(length(u))})
Base.@pure DiffCache(u::AbstractArray,alg) = DiffCache(eltype(u),length(u),Val{get_chunksize(alg)})
Base.@pure DiffCache{CS}(u::AbstractArray,T::Type{Val{CS}}) = DiffCache(eltype(u),length(u),T)

get_du{T<:Dual}(dc::DiffCache, ::Type{T}) = dc.dual_du
get_du(dc::DiffCache, T) = dc.du

type MutableReference{T}
val::T
end

Base.@pure function autodiff_setup(f!, initial_x::Vector,alg)
  autodiff_setup(f!, initial_x,Val{determine_chunksize(initial_x,alg)})
end

Base.@pure function determine_chunksize(u,alg)
  if get_chunksize(alg) != 0
    return get_chunksize(alg)
  else
    return ForwardDiff.pickchunksize(length(u))
  end
end

Base.@pure function autodiff_setup{CS}(f!, initial_x::Vector,chunk_size::Type{Val{CS}})

    permf! = (fx, x) -> f!(x, fx)

    fx2 = copy(initial_x)
    jac_cfg = ForwardDiff.JacobianConfig{CS}(initial_x, initial_x)
    g! = (x, gx) -> ForwardDiff.jacobian!(gx, permf!, fx2, x, jac_cfg)

    fg! = (x, fx, gx) -> begin
        jac_res = DiffBase.DiffResult(fx, gx)
        ForwardDiff.jacobian!(jac_res, permf!, fx2, x, jac_cfg)
        DiffBase.value(jac_res)
    end

    return DifferentiableMultivariateFunction(f!, g!, fg!)
end

#Base.getindex(m::MutableReference)=m.val
#Base.setindex!(m::MutableReference{T}, v::{T},I...) = m.val=v

realtype{T}(::Type{T}) = T
realtype{T}(::Type{Complex{T}}) = T


text:
function collatz_steps(n)
    if n < 1
        throw(DomainError("number needs to be >= 1"))
    end

    i = 0
    while n > 1
        n = iseven(n) ? div(n, 2) : 3n + 1
        i += 1
    end
    return i
end

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
# v0.6

function shellsort!(a::Array{Int})::Array{Int}
    incr = div(length(a), 2)
    while incr > 0
        for i in incr+1:length(a)
            j = i
            tmp = a[i]
            while j > incr && a[j - incr] > tmp
                a[j] = a[j-incr]
                j -= incr
            end
            a[j] = tmp
        end
        if incr == 2
            incr = 1
        else
            incr = floor(Int, incr * 5.0 / 11)
        end
    end
    return a
end

x = rand(1:10, 10)
@show x shellsort!(x)
@assert issorted(x)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module FixedEffectsBinarization

using StatsFuns.logistic
using DataFrames
using Optim: optimize, minimizer, LBFGS, NewtonTrustRegion
using Distributions
using Plots

include("strucs.jl")
include("datageneration.jl")
include("datahandling.jl")
include("estimation.jl")
include("postestimation.jl")

end


text:
function _make_features(layers, filters, batch_norm)
    featurizer = Sequential()
    for (filter, num) in zip(filters, layers)
        for _ = 1:num
            add(featurizer, Conv2D(filter, (3, 3), padding = 1))
            if batch_norm
                add(featurizer, BatchNorm())
            end
            add(featurizer, Activation(relu))
        end
        add(featurizer, MaxPool2D())
    end
    return featurizer
end

struct VGG <: Photon.Layer
    features
    output

    function VGG(layers, filters; classes = 1000, batch_norm = false)
        features = _make_features(layers, filters, batch_norm)
        add(features, Flatten())
        add(features, Dense(4096, relu))
        add(features, Dropout(0.5))
        add(features, Dense(4096, relu))
        add(features, Dropout(0.5))
        output = Dense(classes)
        new(features, output)
    end
end

function (model::VGG)(x)
    x = model.features(x)
    model.output(x)
end

vgg_spec = Dict(
    11 => ([1, 1, 2, 2, 2], [64, 128, 256, 512, 512]),
    13 => ([2, 2, 2, 2, 2], [64, 128, 256, 512, 512]),
    16 => ([2, 2, 3, 3, 3], [64, 128, 256, 512, 512]),
    19 => ([2, 2, 4, 4, 4], [64, 128, 256, 512, 512]),
)

VGG11(; classes = 1000, batch_norm = true) =
    VGG(vgg_spec[11]..., classes = classes, batch_norm = batch_norm)
VGG13(; classes = 1000, batch_norm = true) =
    VGG(vgg_spec[13]..., classes = classes, batch_norm = batch_norm)
VGG16(; classes = 1000, batch_norm = true) =
    VGG(vgg_spec[16]..., classes = classes, batch_norm = batch_norm)
VGG19(; classes = 1000, batch_norm = true) =
    VGG(vgg_spec[19]..., classes = classes, batch_norm = batch_norm)


text:
export divides

"""
    divides(t1::AbstractTermLike, t2::AbstractTermLike)

Returns whether the monomial of t1 divides the monomial of t2.

### Examples

Calling `divides(2x^2y, 3xy)` should return false because `x^2y` does not divide `xy` since `x` has a degree 2 in `x^2y` which is greater than the degree of `x` on `xy`.
However, calling `divides(3xy, 2x^2y)` should return true.
"""
function divides(t1::AbstractTermLike, t2::AbstractTermLike)
    divides(monomial(t1), monomial(t2))
end
divides(t1::AbstractVariable, t2::AbstractVariable) = t1 == t2

Base.gcd(m1::AbstractMonomialLike, m2::AbstractMonomialLike) = mapexponents(min, m1, m2)
Base.lcm(m1::AbstractMonomialLike, m2::AbstractMonomialLike) = mapexponents(max, m1, m2)

# _div(a, b) assumes that b divides a
_div(m1::AbstractMonomialLike, m2::AbstractMonomialLike) = mapexponents(-, m1, m2)
function _div(t::AbstractTerm, m::AbstractMonomial)
    coefficient(t) * _div(monomial(t), m)
end
function _div(t1::AbstractTermLike, t2::AbstractTermLike)
    (coefficient(t1) / coefficient(t2)) * _div(monomial(t1), monomial(t2))
end

Base.div(f::APL, g::Union{APL, AbstractVector{<:APL}}; kwargs...) = divrem(f, g; kwargs...)[1]
Base.rem(f::APL, g::Union{APL, AbstractVector{<:APL}}; kwargs...) = divrem(f, g; kwargs...)[2]

proddiff(x, y) = x/y - x/y
function Base.divrem(f::APL{T}, g::APL{S}; kwargs...) where {T, S}
    rf = convert(polynomialtype(f, Base.promote_op(proddiff, T, S)), f)
    q = r = zero(rf)
    lt = leadingterm(g)
    rg = removeleadingterm(g)
    lm = monomial(lt)
    while !iszero(rf)
        ltf = leadingterm(rf)
        if isapproxzero(ltf; kwargs...)
            rf = removeleadingterm(rf)
        elseif divides(lm, ltf)
            qt = _div(ltf, lt)
            q += qt
            rf = removeleadingterm(rf) - qt * rg
        elseif lm > monomial(ltf)
            # Since the monomials are sorted in decreasing order,
            # lm is larger than all of them hence it cannot divide any of them
            r += rf
            break
        else
            r += ltf
            rf = removeleadingterm(rf)
        end
    end
    q, r
end
function Base.divrem(f::APL{T}, g::AbstractVector{<:APL{S}}; kwargs...) where {T, S}
    rf = convert(polynomialtype(f, Base.promote_op(proddiff, T, S)), f)
    r = zero(rf)
    q = similar(g, typeof(rf))
    for i in eachindex(q)
        q[i] = zero(rf)
    end
    lt = leadingterm.(g)
    rg = removeleadingterm.(g)
    lm = monomial.(lt)
    useful = BitSet(eachindex(g))
    while !iszero(rf)
        ltf = leadingterm(rf)
        if isapproxzero(ltf; kwargs...)
            rf = removeleadingterm(rf)
            continue
        end
        divisionoccured = false
        for i in useful
            if divides(lm[i], ltf)
                qt = _div(ltf, lt[i])
                q[i] += qt
                rf = removeleadingterm(rf) - qt * rg[i]
                divisionoccured = true
                break
            elseif lm[i] > monomial(ltf)
                # Since the monomials are sorted in decreasing order,
                # lm is larger than all of them hence it cannot divide any of them
                delete!(useful, i)
            end
        end
        if !divisionoccured
            if isempty(useful)
                r += rf
                break
            else
                r += ltf
                rf = removeleadingterm(rf)
            end
        end
    end
    q, r
end

function Base.gcd(p::AbstractPolynomialLike{T}, q::AbstractPolynomialLike{S}) where {T, S}
    if isapproxzero(q)
        convert(polynomialtype(p, Base.promote_op(proddiff, T, S)), p)
    else
        gcd(q, rem(p, q))
    end
end
Base.lcm(p::AbstractPolynomialLike, q::AbstractPolynomialLike) = p * div(q, gcd(p, q))


text:
import Robotics: skew

using LinearAlgebra

puma560_dhparams(q) =
         [[  pi/2,        0,         0,   q[1]],
          [     0,   0.4318,         0,   q[2]],
          [ -pi/2,   0.0203,   0.15005,   q[3]],
          [  pi/2,        0,    0.4318,   q[4]],
          [ -pi/2,        0,         0,   q[5]],
          [     0,        0,         0,   q[6]]]

dof = 6

m1 = 0
m2 = 17.4
m3 = 4.8
m4 = 0.82
m5 = 0.34
m6 = 0.09
m = [m1, m2, m3, m4, m5, m6]

r1 = [0., 0., 0.]
r2 = [ -.3638,  .006,    .2275]
r3 = [ -.0203,  -.0141,  .070]
r4 = [ 0,    .019,    0]
r5 = [0.,   0.,   0. ]
r6 = [0.,   0.,   .032]
r = [r1, r2, r3, r4, r5, r6]

# %        Ixx     Iyy      Izz    -Ixy     -Iyz     -Ixz
cI1 = [0.,      0.35,    0.,      0.,   0.,   0.]
cI2 = [ .13,     .524,    .539,   0.,   0.,   0.]
cI3 = [ .066,    .086,    .0125,  0.,   0.,   0.]
cI4 = [1.8e-3,  1.3e-3,  1.8e-3,  0.,   0.,   0.]
cI5 = [ .3e-3,   .4e-3,   .3e-3,  0.,   0.,   0.]
cI6 = [ .15e-3,  .15e-3,  .04e-3, 0.,   0.,   0.]
cI = [cI1, cI2, cI3, cI4, cI5, cI6]

cI2sI(cI) = [cI[1], -cI[4], -cI[6], cI[2], -cI[5], cI[3]]
If = [vectosymm(cI2sI(ci)) for ci in cI]

l = Vector{Float64}[m[i]*r[i] for i in 1:dof]

symmetrize! = x->LinearAlgebra.copytri!(x, 'U')
L = Matrix{Float64}[symmetrize!(If[i] + m[i]*skew(r[i])'*skew(r[i])) for i in 1:dof]

### Dynamic parameters ###

dynp, format = dynp_init( L=(L,:symm3), l=(l,:vec3), m=(m,:scalar) )
dp = dynp_vect2dict(dynp, dof, format)

@test dp[:L] == L
@test dp[:l] == l
@test dp[:m] == m

@test dynp_dict2vect(dp, format) == dynp

dynparms_ok = [0.0, -0, -0, 0.34999999999999998, -0, 0.0, 0.0, 0.0, 0.0, 0, 1.03118515, 0.037980719999999996, 1.4401022999999999,
                    3.7274564059999999, -0.023751000000000001, 2.8425240559999998, -6.33012, 0.10439999999999999, 3.9585, 17.4,
                    0.090474288000000014, -0.0013739039999999998, 0.0068208000000000001, 0.111498032, 0.0047375999999999998,
                    0.015432319999999999, -0.09743999999999998, -0.06767999999999999, 0.336, 4.8, 0.0020960200000000001, -0, -0,
                    0.0012999999999999999, -0, 0.0020960200000000001, 0.0, 0.015579999999999998, 0.0, 0.82, 0.00029999999999999997, -0, -0,
                    0.00040000000000000002, -0, 0.00029999999999999997, 0.0, 0.0, 0.0, 0.34, 0.00024216, -0, -0, 0.00024216, -0,
                    4.0000000000000003e-05, 0.0, 0.0, 0.0028799999999999997, 0.09]

@test dynp ≈ dynparms_ok

### Geometry ###

q_test = [0.7504516826728697, 0.8395156106908136, 0.16851233582594916, 0.3849629637427072, 0.5252993946810777, 0.6701207256444748]
dq_test = [0.24721855939629367, 0.9805915670454258, 0.9895299755642817, 0.7861135739668947, 0.273842245476577, 0.17182358900767503]
ddq_test = [0.707405815485141, 0.25295715193420953, 0.9763909835998361, 0.8412822676113918, 0.4867768296473465, 0.11480270540937143]

Ti_test = interlink_transfs_dhstd(puma560_dhparams(q_test));
Ti_inv_test = [inv_transfmat(T) for T in Ti_test];
T_test = transfs_tobase(Ti_test);
S_test = interlink_lie_S_rotjoint_dhstd(puma560_dhparams(q_test));

T_ok = [-0.655113870655343 -0.474277925274361 -0.588120962109346 0.0540498757911011
         0.524340309498464  0.275038163391149 -0.805866768463298 -0.154761559833806
         0.543960528264717 -0.836310025215453 0.0685017183295358  0.568944734102513
                       0.0                0.0                0.0                1.0]

@test T_test[end] ≈ T_ok

### Dynamics ###

g = [0, 0, -9.81]

tau_test = rne_park(q_test, dq_test, ddq_test, Ti_inv_test, S_test, g, dp[:L], dp[:l], dp[:m]);

tau_ok = [  0.986185688341252
             16.1055550633721
            -6.54839494827661
          0.00626452648190803
          -0.0208972707132494
          5.59805923645945e-5]

@test tau_test ≈ tau_ok


text:
include("euler/euler.jl")
using .Math: fast_max
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute(number::Int)::Int
    numbers_below_10 = length.(["one", "two", "three", "four", "five", "six", "seven", "eight",
                                "nine"])
    numbers_below_20 = length.(["ten", "eleven", "twelve", "thirteen", "fourteen", "fifteen",
                                "sixteen", "seventeen", "eighteen", "nineteen"])
    tenth = length.(["twenty", "thirty", "forty", "fifty", "sixty", "seventy", "eighty", "ninety"])
    sum_numbers_below_10, sum_numbers_below_20 = sum(numbers_below_10), sum(numbers_below_20)
    sum_tenths = sum(tenth)
    sum_hundreds = length("hundred")
    sum_and = length("and")
    sum_thousands = length("thousand")
    thousands, hundreds = number ÷ 1_000, number % 1_000 ÷ 100
    tenths, digits = number % 100 ÷ 10, number % 10
    sum_numbers_below_10 = 190thousands * sum_numbers_below_10 +
        (thousands > 0 ? 1_000sum(numbers_below_10[1:thousands - 1]) +
            (100hundreds + 10tenths + digits + 1)numbers_below_10[thousands] : 0) +
        (hundreds > 0 ? 100sum(numbers_below_10[1:hundreds - 1]) + (10tenths + digits + 1) *
            numbers_below_10[hundreds] : 0) +
        9hundreds * sum_numbers_below_10 +
        (tenths ≠ 1 ? fast_max((tenths - 1)sum_numbers_below_10, 0) +
            sum(numbers_below_10[1:digits]) : sum_numbers_below_10)
    sum_numbers_below_20 = 10thousands * sum_numbers_below_20 + hundreds * sum_numbers_below_20 +
        (tenths > 1 ? sum_numbers_below_20 : 0) +
        (tenths == 1 ? sum(numbers_below_20[1:digits]) : 0)
    sum_tenths = 100thousands * sum_tenths + 10hundreds * sum_tenths +
        (tenths > 2 ? 10sum(tenth[1:tenths - 2]) : 0) +
        (tenths == 2 || tenths == 9 ? (digits + 1)tenth[tenths - 1] : 0)
    sum_hundreds = 900thousands * sum_hundreds + fast_max(100 * (hundreds - 1)sum_hundreds, 0) +
        (hundreds > 0 ? (10tenths + digits + 1)sum_hundreds : 0)
    sum_and = 891thousands * sum_and + fast_max(99 * (hundreds - 1)sum_and, 0) +
        (hundreds > 0 ? (10tenths + digits)sum_and : 0)
    sum_thousands = fast_max(1_000 * (thousands - 1), 0) +
        (thousands > 0 ? (100hundreds + 10tenths + digits + 1)sum_thousands : 0)
    return sum_numbers_below_10 + sum_numbers_below_20 + sum_tenths + sum_hundreds + sum_and +
        sum_thousands
end

compute(5)

compute(1_000)

@benchmark compute(1_000)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
export time_window2index_window, time2seconds, times_of_flow,
    plot_res, quantile_filter!, lowpass_filter2,
    lowpass_filter, median_filter

"Plot results"
function plot_res end

time2seconds(time,ex) = time2seconds(time, ex.experiment_start)
time2seconds(time,experiment_start::DateTime) =
    Dates.value.(Dates.Nanosecond.((time - Time(experiment_start))))/1e9

"""
    time_window2index_window(times, window)

Return the indices of all t∈times for which window[1]<=t<window[2].
"""
time_window2index_window(time, window) =
    findfirst(t->t>=window[1],time):findfirst(t->t>window[end], time)

"""
    times_of_flow(time, ex)
    times_of_flow(time, experiment_running)

Return times (in sec since start) and indices when water was flowing.
"""
times_of_flow(time, ex) = times_of_flow(time, ex.experiment_start, ex.experiment_running)
function times_of_flow(time, experiment_start::DateTime,
                       experiment_running::AbstractVector)
    irun = vcat([time_window2index_window(time, w) for w in experiment_running]...)
    t = time2seconds(time, experiment_start)
    dt = [0;diff(t)]
    trun = cumsum(dt[irun])
    return trun, irun
end

###############
# Filtering

"Image median filter"
median_filter(img, window=(3,3)) = Images.mapwindow(median!, img, window)

"""
Lowpass filter with a cutoff at `wavelength`,
expressed as ratio of wavelength/length(lines).
"""
function lowpass_filter(lines, wavelength::Float64=0.03)
    wavelength_pix = round(Int, wavelength*size(lines,1))
    frq = 2/wavelength_pix
    responsetype = DSP.Lowpass(frq)
    designmethod = DSP.Butterworth(4)
    round.(Int, DSP.filtfilt(DSP.digitalfilter(responsetype, designmethod), lines))
end

"""
Lowpass filter with a cutoff at `wavelength`,
expressed as ratio of wavelength/length(lines).
"""
function lowpass_filter2(v, wavelength_samples)
    frq = 2/wavelength_samples
    responsetype = DSP.Lowpass(frq)
    designmethod = DSP.Butterworth(4)
    DSP.filtfilt(DSP.digitalfilter(responsetype, designmethod), v)
end

"Filters outliers above certain quantiles"
function quantile_filter!(vec,
                          quant_cut=0.1, quant_offset=0, quant_factor=1,
                          pre_window=length(vec)÷10, post_window = pre_window+1,
                          niter = 10,
                          verbose=false)
    lenu = lenl = -1
    for i=1:niter
        upper, medu = find_spikes_quantile(vec, >,
                                           quant_cut, quant_offset, quant_factor,
                                           pre_window, post_window;
                                           verbose=verbose)
        lower, medl = find_spikes_quantile(vec, <,
                                           1-quant_cut, quant_offset, quant_factor,
                                           pre_window, post_window;
                                           verbose=verbose)

        lenu==length(upper) && break
        lenl==length(lower) && break
        lenu = length(upper)
        lenl = length(lower)
        vec[upper] = medu
        vec[lower] = medl
    end
    vec
end