text:
using LinearAlgebra, SparseArrays, Statistics, StatsBase

# This expresses the probability the frequency of cells that carry variant a certain variant c_j is at least f.
φ(f,u,ρ,c) = (u / (u-log(1+(ρ-1)*f)))^c

# This calculates the expected number of cells in the population at a given time (measured here based on the discussion
# of the paper (see Equations 9 and 10 for details)
N(s,ρ,u) = exp((ρ/u -1/u)*s) + 2*sinh((1/u - ρ/u)*s)/ρ

function gmap(seed,u,ρ, debug = false)
    # Calculation of the approximate mapping g(c_j)
    # For details see Equation (10) of the manuscript
    if debug
        println("exp part = ",exp((ρ/u - 1/u)*seed))
        println("sinh part = ",2*sinh((1/u - ρ/u)*seed)/ρ)
    end
    return seed == 0 ? 0 : ceil(u*( exp((ρ/u - 1/u)*seed) + 2*sinh((1/u - ρ/u)*seed)/ρ))
end

function popEstimate(seed,u,ρ, debug = false)
    # This function returns an estimate of the expected number of cells in the population,
    # at the time  each of the seeding-cell-specific variants appeared. For details see the manuscript (Equations 9 and 10)
    if debug
        println("exp part = ",exp((ρ/u - 1/u)*seed))
        println("sinh part = ",2*sinh((1/u - ρ/u)*seed)/ρ)
    end
    return exp((ρ/u - 1/u)*seed) + 2*sinh((1/u - ρ/u)*seed)/ρ
end

function calculate_sampling_probabilities_populations(u, ts1, ts2, t_det, ρ, ρ1, ρ2 )
    # This refers to the 3rd population based seeding strategy
    # For context see the related discussion in Section 4.1.3 of the paper.
    p0 = zeros(t_det)
    p1 = zeros(t_det)
    p2 = zeros(t_det)
    for k=1:t_det
        if k < ts1
            N0 = N(k,ρ,u)
            p0[k] = 1
            continue
        end
        if k < ts2
            N0 = N(k,ρ,u)
            N1 = N(k-ts1,ρ1,u)
            p0[k] = N0 / (N0 + N1)
            p1[k] = N1 / (N0 + N1)
            continue
        end
        N0 = N(k,ρ,u)
        N1 = N(k-ts1,ρ1,u)
        N2 = N(k-ts2,ρ2,u)
        p0[k] = (N0) / (N0+N1+N2)
        p1[k] = (N1) / (N0+N1+N2)
        p2[k] = (N2) / (N0+N1+N2)
    end
    p0, p1, p2
end

function Bm(f,u,ρ,seed, debug = false)
    Bm = 0
    for j=0:seed
        Ptemp = 1.0
        if j != seed
            p = (1-φ(f,u,ρ,gmap(j+1,u,ρ))) - (1-φ(f,u,ρ,gmap(j,u,ρ)))
        else
            p = φ(f,u,ρ,gmap(j,u,ρ))
        end
        Bm += (seed-j)*p
        if debug
            println("P(U_j+1) = ", (1-φ(f,u,ρ,gmap(j+1,u,ρ))))
            println("P(U_j) = ", (1-φ(f,u,ρ,gmap(j,u,ρ))))
            println("Bm = ", Bm)
        end
    end
    return Bm
end

function BmGeneric(d, k_max, debug = false)
    Bm = 0
    BmPMF = zeros(k_max+1)
    for j=0:k_max
        if j != k_max
            p = (1-d[j+2]) - (1-d[j+1])
        else
            p = d[k_max]
        end
        Bm += (k_max-j)*p
        BmPMF[k_max-j+1] = p
    end
    return Bm, BmPMF
end

function BmDist(f,u,ρ,seed, debug = false)
    # Calculates the B_md distribution (See equation 6, and the related discussion in the paper)
    Bm = []
    for j=0:seed
        Ptemp = 1.0
        if j != seed
            p = (1-φ(f,u,ρ,gmap(j+1,u,ρ))) - (1-φ(f,u,ρ,gmap(j,u,ρ)))
        else
            p = φ(f,u,ρ,gmap(j,u,ρ))
        end
        append!(Bm,p)
        if debug
            println("P(U_j+1) = ", (1-φ(f,u,ρ,gmap(j+1,u,ρ))))
            println("P(U_j) = ", (1-φ(f,u,ρ,gmap(j,u,ρ))))
            println("Bm = ", Bm)
        end
    end
    return Bm
end

function BmVar(f,u,ρ,seed)
    # To measure the innate variance of the random variable Bmd that follows the above distribution BmDist
    Bm = 0
    Bm2 = 0
    for j=0:seed
        Ptemp = 1.0
        if j != seed
            p = (1-φ(f,u,ρ,gmap(j+1,u,ρ))) - (1-φ(f,u,ρ,gmap(j,u,ρ)))
        else
            p = φ(f,u,ρ,gmap(j,u,ρ))
        end
        Bm += (seed-j)*p
        Bm2 += (seed-j)^2*p
    end
    return Bm2 - Bm^2
end

# The following are simple Helper functions necessary to facilitate the production of the plots

function kmaxEstimate(population,u,ρ)
    k = 1
    while exp((ρ/u - 1/u)*k) + 2*sinh((1/u - ρ/u)*k)/ρ < population
        k += 1
    end
    return k
end

function findtd(u, ρ, population)
    s = 1
    while exp((ρ / u - 1 / u) * s) + 2 * sinh((1 / u - ρ / u) * s) / ρ < population
        s += 1
    end
    return s
end

function finds_i(u, ρ, time, population)
    x = 0.0
    ρ1 = ρ / (x + 1)
    while exp((ρ1 / u - 1 / u) * time) + 2 * sinh((1 / u - ρ1 / u) * time) / ρ1 < population
        x += 0.0005
        ρ1 = ρ / (x + 1)
    end
    return x
end

function findts_i(u, ρ, td0, population)
    t = 0
    while exp((ρ / u - 1 / u) * t) + 2 * sinh((1 / u - ρ / u) * t) / ρ < population
        t += 1
    end
    return td0 - t # returns the type origination time measured in [0..td]  (integer value, i.e., the number of mutations...)
end

max_stars_repo_name: SunPathLab/Comet

max_issues_repo_name: SunPathLab/Comet

max_forks_repo_name: SunPathLab/Comet


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    ExponentialVariogram(range=r, sill=s, nugget=n)
    ExponentialVariogram(ball; sill=s, nugget=n)

An exponential variogram with range `r`, sill `s` and nugget `n`.
Optionally, use a custom metric `ball`.
"""
struct ExponentialVariogram{T,B} <: Variogram
  sill::T
  nugget::T
  ball::B
end

ExponentialVariogram(ball; sill=1.0, nugget=0.0) =
  ExponentialVariogram(sill, nugget, ball)

ExponentialVariogram(; range=1.0, sill=1.0, nugget=0.0) =
  ExponentialVariogram(sill, nugget, MetricBall(range))

function (γ::ExponentialVariogram)(h)
  r = radius(γ.ball)
  s = γ.sill
  n = γ.nugget
  (s - n) * (1 - exp(-3(h/r))) + (h > 0) * n
end

isstationary(::Type{<:ExponentialVariogram}) = true

ext: jl


text:
# ---
# title: 764. Largest Plus Sign
# id: problem764
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/largest-plus-sign/description/>
# hidden: true
# ---
#
# In a 2D `grid` from (0, 0) to (N-1, N-1), every cell contains a `1`, except
# those cells in the given list `mines` which are `0`. What is the largest axis-
# aligned plus sign of `1`s contained in the grid? Return the order of the plus
# sign. If there is none, return 0.
#
# An " _axis-aligned plus sign of`1`s_ of order **k** " has some center
# `grid[x][y] = 1` along with 4 arms of length `k-1` going up, down, left, and
# right, and made of `1`s. This is demonstrated in the diagrams below. Note that
# there could be `0`s or `1`s beyond the arms of the plus sign, only the
# relevant area of the plus sign is checked for 1s.
#
# **Examples of Axis-Aligned Plus Signs of Order k:**
#
#
#
#     Order 1:
#     000
#     0 **1** 0
#     000
#
#     Order 2:
#     00000
#     00 **1** 00
#     0 **111** 0
#     00 **1** 00
#     00000
#
#     Order 3:
#     0000000
#     000 **1** 000
#     000 **1** 000
#     0 **11111** 0
#     000 **1** 000
#     000 **1** 000
#     0000000
#
#
# **Example 1:**
#
#
#
#     Input: N = 5, mines = [[4, 2]]
#     Output: 2
#     Explanation:
#     11111
#     11111
#     1 **1** 111
#     **111** 11
#     1 **1** 011
#     In the above grid, the largest plus sign can only be order 2.  One of them is marked in bold.
#
#
# **Example 2:**
#
#
#
#     Input: N = 2, mines = []
#     Output: 1
#     Explanation:
#     There is no plus sign of order 2, but there is of order 1.
#
#
# **Example 3:**
#
#
#
#     Input: N = 1, mines = [[0, 0]]
#     Output: 0
#     Explanation:
#     There is no plus sign, so return 0.
#
#
# **Note:**
#
#   1. `N` will be an integer in the range `[1, 500]`.
#   2. `mines` will have length at most `5000`.
#   3. `mines[i]` will be length 2 and consist of integers in the range `[0, N-1]`.
#   4. _(Additionally, programs submitted in C, C++, or C# will be judged with a slightly smaller time limit.)_
#
#
## @lc code=start
using LeetCode

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


text:
#
# GUID and CLSID
#

immutable t_data4
    d1::Cuchar; d2::Cuchar; d3::Cuchar; d4::Cuchar; d5::Cuchar; d6::Cuchar; d7::Cuchar; d8::Cuchar
end

immutable GUID
    data1::Culong
    data2::Cushort
    data3::Cushort
    data4::t_data4
end
GUID() = GUID(0,0,0,t_data4(0,0,0,0,0,0,0,0))

immutable CLSID
    guid::GUID
    CLSID(guid) = new(guid)
end
CLSID() = CLSID(GUID())

#
# Win API aliases
#
immutable IID
    guid::GUID
    IID(guid) = new(guid)
end
IID() = IID(GUID())

typealias REFIID Ptr{IID}
typealias LPIID Ptr{GUID}
typealias REFCLSID Ptr{CLSID}
typealias LPCLSID Ptr{CLSID}
typealias LPOLESTR Ptr{Cwchar_t}
typealias LPCOLESTR Ptr{Cwchar_t}
typealias DWORD Culong

typealias LPUNKNOWN Ptr{Void} # TODO: make this better

######## Base COM API start

baremodule COINIT
    APARTMENTTHREADED  = 0x2
    MULTITHREADED      = 0x0
    DISABLE_OLE1DDE    = 0x4
    SPEED_OVER_MEMORY  = 0x8
end

typealias HResult Uint32
baremodule HRESULT
    S_OK          =   0x00000000
    S_FALSE       =   0x00000001
    E_UNEXPECTED  =   0x8000FFFF
    E_NOTIMPL     =   0x80004001
    E_OUTOFMEMORY =   0x8007000E
    E_INVALIDARG  =   0x80070057
    E_NOINTERFACE =   0x80004002
    E_POINTER     =   0x80004003
    E_HANDLE      =   0x80070006
    E_ABORT       =   0x80004004
    E_FAIL        =   0x80004005
    E_ACCESSDENIED=   0x80070005
    E_PENDING     =   0x8000000A
end # module HRESULT

baremodule REGDB
    E_FIRST           = 0x80040150
    E_LAST            = 0x8004015F
    S_FIRST           = 0x00040150
    S_LAST            = 0x0004015F
    E_READREGDB       = 0x80040150
    E_WRITEREGDB      = 0x80040151
    E_KEYMISSING      = 0x80040152
    E_INVALIDVALUE    = 0x80040153
    E_CLASSNOTREG     = 0x80040154
    E_IIDNOTREG       = 0x80040155
    E_BADTHREADINGMODEL = 0x80040156
end

baremodule CLSCTX
    import Base.|

    INPROC_SERVER           = 0x1
    INPROC_HANDLER          = 0x2
    LOCAL_SERVER            = 0x4
    INPROC_SERVER16         = 0x8
    REMOTE_SERVER           = 0x10
    INPROC_HANDLER16        = 0x20
    RESERVED1               = 0x40
    RESERVED2               = 0x80
    RESERVED3               = 0x100
    RESERVED4               = 0x200
    NO_CODE_DOWNLOAD        = 0x400
    RESERVED5               = 0x800
    NO_CUSTOM_MARSHAL       = 0x1000
    ENABLE_CODE_DOWNLOAD    = 0x2000
    NO_FAILURE_LOG          = 0x4000
    DISABLE_AAA             = 0x8000
    ENABLE_AAA              = 0x10000
    FROM_DEFAULT_CONTEXT    = 0x20000
    ACTIVATE_32_BIT_SERVER  = 0x40000
    ACTIVATE_64_BIT_SERVER  = 0x80000
    ENABLE_CLOAKING         = 0x100000
    APPCONTAINER            = 0x400000
    ACTIVATE_AAA_AS_IU      = 0x800000
    PS_DLL                  = 0x80000000

    SERVER = INPROC_SERVER | LOCAL_SERVER | REMOTE_SERVER
    ALL    = INPROC_HANDLER | SERVER
end # module CLSCTX

################################################################################

#
# CLSID
#

function CLSIDFromString(id::String)
    out = [GUID()]
    res = ccall( (:CLSIDFromString, l_ole32), Uint32,
                (LPCOLESTR, LPCLSID), utf16(id), pointer(out))
    res != HRESULT.S_OK && error("CLSIDFromString: Unable to convert $id to CLSID")
    return CLSID(out[1])
end
CLSID(id::String) = CLSIDFromString(id)

#
# Interfaces
#

abstract IUnknown
abstract IDispatch <: IUnknown
#TODO: more types

#
# Builtin IIDs
#

const BaseIIDs = {
    IUnknown => CLSID("{00000000-0000-0000-C000-000000000046}"),
    #RecordInfo => CLSID("{0000002F-0000-0000-C000-000000000046}"),
    #IRecordInfo => CLSID("{0000002F-0000-0000-C000-000000000046}"),
    IDispatch => CLSID("{00020400-0000-0000-C000-000000000046}"),
    #ITypeComp => CLSID("{00020403-0000-0000-C000-000000000046}"),
    #ITypeInfo => CLSID("{00020401-0000-0000-C000-000000000046}"),
    #ITypeInfo2 => CLSID("{00020412-0000-0000-C000-000000000046}"),
    #ITypeLib => CLSID("{00020402-0000-0000-C000-000000000046}"),
    #ITypeLib2 => CLSID("{00020411-0000-0000-C000-000000000046}"),
    None => CLSID("{00000000-0000-0000-0000-000000000000}")
    }
getindex(::Type{IID}, x::Type) = IID(BaseIIDs[x])
getindex(::Type{CLSID}, x::Type) = BaseIIDs[x]

function getindex(t::IID, x::Ptr{Void})
    k = None
    for k in keys(BaseIIDs)
        if BaseIIDs[k].guid == t break
        else continue
        end
    end
    convert(Ptr{k}, x)
end
#
# Misc
#

immutable LCID
    d::Uint32
end
convert(::Type{Uint32}, t::LCID) = t.d


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

import Test
import DensityInterface
import Documenter

Test.@testset "Package DensityInterface" begin
    include("test_interface.jl")

    # doctests
    Documenter.DocMeta.setdocmeta!(
        DensityInterface,
        :DocTestSetup,
        quote
            using DensityInterface
            object = logfuncdensity(x -> x^2)
            log_f = logdensityof(object)
            f = densityof(object)
            x = 4.2
        end;
            recursive=true,
    )
    Documenter.doctest(DensityInterface)
end # testset


text:
# Parser
# ======

# A set of parser parameters.
struct ParserParameters
    delim::UInt8
    quot::UInt8
    trim::Bool
    skip::Int
    colnames::Union{Vector{Symbol},Nothing}
    chunksize::Int

    function ParserParameters(delim::Char, quot::Char, trim::Bool, skip::Integer, header::Any, chunksize::Integer)
        if delim ∉ ALLOWED_DELIMITERS
            throw(ArgumentError("delimiter $(repr(delim)) is not allowed"))
        elseif quot ∉ ALLOWED_QUOTECHARS
            throw(ArgumentError("quotation character $(repr(quot)) is not allowed"))
        elseif delim == quot
            throw(ArgumentError("delimiter and quotation character cannot be the same character"))
        elseif delim == ' ' && trim
            throw(ArgumentError("delimiting with space and space trimming are exclusive"))
        elseif quot == ' ' && trim
            throw(ArgumentError("quoting with space and space trimming are exclusive"))
        elseif skip < 0
            throw(ArgumentError("skip cannot be negative"))
        elseif chunksize < 0
            throw(ArgumentError("chunk size cannot be negative"))
        elseif chunksize ≥ MAX_TOKEN_START
            throw(ArgumentError("chunk size must be less than $(MAX_TOKEN_START)"))
        end
        if header != nothing
            colnames = Symbol.(collect(header))
        else
            colnames = nothing
        end
        return new(
            UInt8(delim),
            UInt8(quot),
            trim,
            skip,
            colnames,
            chunksize,
        )
    end
end

# Integer
# -------

function fillcolumn!(col::Vector{Int}, nvals::Int, mem::Memory, tokens::Matrix{Token}, c::Int, quot::UInt8)
    @inbounds for i in 1:nvals
        start, length = location(tokens[c,i])
        col[end-nvals+i] = parse_integer(mem, start, length)
    end
    return col
end

function fillcolumn!(col::Vector{Union{Int,Missing}}, nvals::Int, mem::Memory, tokens::Matrix{Token}, c::Int, quot::UInt8)
    @inbounds for i in 1:nvals
        t = tokens[c,i]
        if ismissing(t)
            col[end-nvals+i] = missing
        else
            start, length = location(t)
            col[end-nvals+i] = parse_integer(mem, start, length)
        end
    end
    return col
end

const SAFE_INT_LENGTH = sizeof(string(typemax(Int))) - 1

@inline function parse_integer(mem::Memory, start::Int, length::Int)
    stop = start + length - 1
    if length > SAFE_INT_LENGTH
        # safe but slow fallback
        buf = IOBuffer()
        for i in start:stop
            write(buf, mem[i])
        end
        return parse(Int, String(take!(buf)))
    end
    i = start
    b = mem[i]
    if b == UInt8('-')
        sign = -1
        i += 1
    elseif b == UInt8('+')
        sign = +1
        i += 1
    else
        sign = +1
    end
    n::Int = 0
    while i ≤ stop
        @inbounds b = mem[i]
        n = 10n + (b - UInt8('0'))
        i += 1
    end
    return sign * n
end

# Float
# -----

function fillcolumn!(col::Vector{Float64}, nvals::Int, mem::Memory, tokens::Matrix{Token}, c::Int, quot::UInt8)
    @inbounds for i in 1:nvals
        start, length = location(tokens[c,i])
        col[end-nvals+i] = parse_float(mem, start, length)
    end
    return col
end

function fillcolumn!(col::Vector{Union{Float64,Missing}}, nvals::Int, mem::Memory, tokens::Matrix{Token}, c::Int, quot::UInt8)
    @inbounds for i in 1:nvals
        t = tokens[c,i]
        if ismissing(t)
            col[end-nvals+i] = missing
        else
            start, length = location(tokens[c,i])
            col[end-nvals+i] = parse_float(mem, start, length)
        end
    end
    return col
end

@inline function parse_float(mem::Memory, start::Int, length::Int)
    return ccall(:jl_strtod_c, Cdouble, (Ptr{UInt8}, Ptr{Cvoid}), mem.ptr + start - 1, C_NULL)
end

function fillcolumn!(col::Vector{Bool}, nvals::Int, mem::Memory, tokens::Matrix{Token}, c::Int, quot::UInt8)
    @inbounds for i in 1:nvals
        start, length = location(tokens[c,i])
        col[end-nvals+i] = parse_bool(mem, start, length)
    end
    return col
end

@inline function parse_bool(mem::Memory, start::Int, length::Int)
    c = mem[start]
    return (c == UInt8('f') || c == UInt8('F')) ? false : true
end

# String
# ------

const N_CACHED_STRINGS = 8

function fillcolumn!(col::Vector{String}, nvals::Int, mem::Memory, tokens::Matrix{Token}, c::Int, quot::UInt8)
    cache = StringCache(N_CACHED_STRINGS)
    usecache = true
    @inbounds for i in 1:nvals
        t = tokens[c,i]
        start, length = location(t)
        if kind(t) & QSTRING != 0
            s = qstring(mem, start, length, quot)
        elseif usecache
            s = allocate!(cache, mem.ptr + start - 1, length % UInt64)
        else
            s = unsafe_string(mem.ptr + start - 1, length)
        end
        col[end-nvals+i] = s
        if usecache && i % 4096 == 0 && 10 * cache.stats.hit < cache.stats.hit + cache.stats.miss
            # stop using cache because the cache hit rate is too low
            usecache = false
        end
    end
    if !usecache
        @debug "Cache is turned off"
    end
    @debug cache
    return col
end

function fillcolumn!(col::Vector{Union{String,Missing}}, nvals::Int, mem::Memory, tokens::Matrix{Token}, c::Int, quot::UInt8)
    cache = StringCache(N_CACHED_STRINGS)
    usecache = true
    @inbounds for i in 1:nvals
        t = tokens[c,i]
        if ismissing(t)
            col[end-nvals+i] = missing
        else
            start, length = location(t)
            if kind(t) & QSTRING != 0
                s = qstring(mem, start, length, quot)
            elseif usecache
                s = allocate!(cache, mem.ptr + start - 1, length % UInt64)
            else
                s = unsafe_string(mem.ptr + start - 1, length)
            end
            col[end-nvals+i] = s
            if usecache && i % 4096 == 0 && 10 * cache.stats.hit < cache.stats.hit + cache.stats.miss
                # stop using cache because the cache hit rate is too low
                usecache = false
            end
        end
    end
    if !usecache
        @debug "Cache is turned off"
    end
    @debug cache
    return col
end

function qstring(mem::Memory, start::Int, length::Int, quot::UInt8)
    i = start
    skip = false
    buf = IOBuffer()
    while i < start + length
        if skip
            skip = false
        else
            x = mem[i]
            write(buf, x)
            skip = x == quot
        end
        i += 1
    end
    return String(take!(buf))
end

# Date
# ----

const DATE_REGEX = r"^\d{4}-\d{2}-\d{2}$"

function parse_date(s::String)
    return Date(s, dateformat"y-m-d")
end

function parse_date(col::Vector{String})
    out = Vector{Date}(undef, length(col))
    for i in 1:length(col)
        out[i] = parse_date(col[i])
    end
    return out
end

function parse_date(col::Vector{Union{String,Missing}})
    out = Vector{Union{String,Missing}}(undef, length(col))
    for i in 1:length(col)
        x = col[i]
        if x === missing
            out[i] = missing
        else
            out[i] = parse_date(x)
        end
    end
    return out
end

function is_date_like(col::Vector{<:Union{String,Missing}})
    i = 1
    n = 0
    while i ≤ length(col) && n < 3
        x = col[i]
        if x isa String
            if !occursin(DATE_REGEX, x)
                return false
            end
            n += 1
        end
        i += 1
    end
    return n > 0
end

# DateTime
# --------

const DATETIME_REGEX = r"^\d{4}-\d{2}-\d{2}T\d{2}:\d{2}:\d{2}(:?\.\d+)?$"

function parse_datetime(s::String)
    return DateTime(s, dateformat"y-m-dTH:M:S.s")
end

function is_datetime_like(col::Vector{<:Union{String,Missing}})
    # Check if the first three strings (if any) are datetime-like.
    i = 1
    n = 0
    while i ≤ length(col) && n < 3
        x = col[i]
        if x isa String
            if !occursin(DATETIME_REGEX, x)
                return false
            end
            n += 1
        end
        i += 1
    end
    return n > 0
end

function parse_datetime(col::Vector{String})
    out = Vector{DateTime}(undef, length(col))
    for i in 1:length(col)
        out[i] = parse_datetime(col[i])
    end
    return out
end

function parse_datetime(col::Vector{Union{String,Missing}})
    out = Vector{Union{DateTime,Missing}}(undef, length(col))
    for i in 1:length(col)
        x = col[i]
        if x === missing
            out[i] = missing
        else
            out[i] = parse_datetime(x)
        end
    end
    return out
end


text:
using FluxTraining
using Pollen
using GraphPlot

project = Pollen.documentationproject(FluxTraining;
    refmodules = [FluxTraining, FluxTraining.Phases], watchpackage=true)

##

Pollen.serve(project)


text:
#==========================================================================================#

# INTERFACE

function fit(M::Type{<: OneStageModel}, MD::Microdata; novar::Bool = false, kwargs...)

    obj = M(MD; kwargs...)

    _fit!(obj, MD.weights)
    novar || _vcov!(obj, getcorr(obj), MD.weights)

    return obj
end

#==========================================================================================#

# ESTIMATES

coef(obj::ParObject)     = switch_stage(obj).β
vcov(obj::ParObject)     = switch_stage(obj).V
stderror(obj::ParObject) = sqrt.(diag(vcov(obj)))
tstat(obj::ParObject)    = coef(obj) ./ stderror(obj)
pval(obj::ParObject)     = 2.0 * normccdf.(abs.(tstat(obj)))

function confint(obj::ParObject, level::Real = 0.95)
    return coef(obj) .+ norminvcdf((1.0 - level) / 2.0) * stderror(obj) .* [1.0 -1.0]
end

#==========================================================================================#

# SUMMARY STATISTICS

nobs(obj::Microdata) = sum(obj.weights)
nobs(obj::AnyModel)  = nobs(switch_stage(obj).sample)

dof(obj::ParModel)          = length(coef(obj))
dof_residual(obj::ParModel) = nobs(obj) - dof(obj)

loglikelihood(obj::MLE)     = _loglikelihood(obj, getweights(obj))
nullloglikelihood(obj::MLE) = _nullloglikelihood(obj::MLE, getweights(obj))
deviance(obj::MLE)          = _deviance(obj::MLE, getweights(obj))
nulldeviance(obj::MLE)      = _nulldeviance(obj::MLE, getweights(obj))

#==========================================================================================#

# PREDICTION

linear_predictor(obj::AnyModel) = linear_predictor(obj, switch_stage(obj).sample)
predict(obj::AnyModel)          = predict(obj, switch_stage(obj).sample)
fitted(obj::AnyModel)           = predict(obj)
residuals(obj::AnyModel)        = residuals(obj, switch_stage(obj).sample)
response(obj::AnyModel)         = getvector(obj, :response)

function residuals(obj::AnyModel, MD::Microdata)
    y  = response(obj)
    r  = predict(obj, MD)
    r .= y .- r
    return r
end

#==========================================================================================#

# COEFFICIENT LABELS

coefnames(obj::ParEstimate) = obj.names

# OUTPUT

function coeftable(obj::ParObject; level::Float64 = 0.95, digits::Int = 4)

    table = formatter(hcat(coef(obj), stderror(obj), tstat(obj), pval(obj)), digits)
    label = [" Estimate", " St. Err.", "  t-stat.", "  p-value"]

    if level > 0.0
        table  = hcat(table, formatter(confint(obj, level), digits))
        label  = vcat(label, ["     C.I.", "($(format("{:.0d}", 100 * level))%)  "])
    end

    return CoefTable(table, label, coefnames(obj))
end

function Base.show(io::IO, obj::ParObject)
    isdefined(switch_stage(obj), :V) ? println(io, coeftable(obj)) : println(io, coef(obj))
end

max_stars_repo_name: lbittarello/Microeconometrics

max_issues_repo_name: lbittarello/Microeconometrics

max_forks_repo_name: lbittarello/Microeconometrics


text:
using Plots

V(x,y) = (1/2)*(x^2 + y^2 + 2*x^2*y - (2/3)*y^3)
x = y = range(-1,stop = 1, length = 1000)
V(x,y) = (1/2)*(x^2 + y^2 + 2*x^2*y - (2/3)*y^3)
x = y = range(-1,stop = 1, length = 1000)
contourf(x,y,V)
contour!(x,y,V,levels=0:1/6:1,title = "Henon-Heiles potential: Intensity map",linecolor = :white)
# savefig("potential_intensity.png")

plot(x,y,V,title = "Henon-Heiles potential: Surface map")
# savefig("potential_surface.png")

max_stars_repo_name: Ved-Mahajan/Henon-Heiles-System

max_issues_repo_name: Ved-Mahajan/Henon-Heiles-System

max_forks_repo_name: Ved-Mahajan/Henon-Heiles-System


text:
function save(net::ELM, fname::AbstractString; compress=true)
    jldopen(fname, "w", compress=compress) do file
        write(file, "machine", net)
    end
end

function load(fname::AbstractString)
    jldopen(fname, "r") do file
        read(file, "machine")
    end
end


text:
# Things that will go into PlotUtils.jl

function optimal_ticks_and_labels(axis::Axis, ticks = nothing)
    amin,amax = axis_limits(axis)

    # scale the limits
    scale = axis[:scale]
    sf = scalefunc(scale)

    # If the axis input was a Date or DateTime use a special logic to find
    # "round" Date(Time)s as ticks
    # This bypasses the rest of optimal_ticks_and_labels, because
    # optimize_datetime_ticks returns ticks AND labels: the label format (Date
    # or DateTime) is chosen based on the time span between amin and amax
    # rather than on the input format
    # TODO: maybe: non-trivial scale (:ln, :log2, :log10) for date/datetime
    if ticks == nothing && scale == :identity
        if axis[:formatter] == dateformatter
            # optimize_datetime_ticks returns ticks and labels(!) based on
            # integers/floats corresponding to the DateTime type. Thus, the axes
            # limits, which resulted from converting the Date type to integers,
            # are converted to 'DateTime integers' (actually floats) before
            # being passed to optimize_datetime_ticks.
            # (convert(Int, convert(DateTime, convert(Date, i))) == 87600000*i)
            ticks, labels = optimize_datetime_ticks(864e5 * amin, 864e5 * amax;
                k_min = 2, k_max = 4)
            # Now the ticks are converted back to floats corresponding to Dates.
            return ticks / 864e5, labels
        elseif axis[:formatter] == datetimeformatter
            return optimize_datetime_ticks(amin, amax; k_min = 2, k_max = 4)
        end
    end

    # get a list of well-laid-out ticks
    if ticks == nothing
        scaled_ticks = optimize_ticks(
            sf(amin),
            sf(amax);
            k_min = 4, # minimum number of ticks
            k_max = 8, # maximum number of ticks
        )[1]
    elseif typeof(ticks) <: Int
        scaled_ticks, viewmin, viewmax = optimize_ticks(
            sf(amin),
            sf(amax);
            k_min = ticks, # minimum number of ticks
            k_max = ticks, # maximum number of ticks
            k_ideal = ticks,
            # `strict_span = false` rewards cases where the span of the
            # chosen  ticks is not too much bigger than amin - amax:
            strict_span = false,
        )
        axis[:lims] = map(invscalefunc(scale), (viewmin, viewmax))
    else
        scaled_ticks = map(sf, (filter(t -> amin <= t <= amax, ticks)))
    end
    unscaled_ticks = map(invscalefunc(scale), scaled_ticks)

    labels = if any(isfinite, unscaled_ticks)
        formatter = axis[:formatter]
        if formatter == :auto
            # the default behavior is to make strings of the scaled values and then apply the labelfunc
            map(labelfunc(scale, backend()), Showoff.showoff(scaled_ticks, :plain))
        elseif formatter == :scientific
            Showoff.showoff(unscaled_ticks, :scientific)
        else
            # there was an override for the formatter... use that on the unscaled ticks
            map(formatter, unscaled_ticks)
        end
    else
        # no finite ticks to show...
        String[]
    end

    # @show unscaled_ticks labels
    # labels = Showoff.showoff(unscaled_ticks, scale == :log10 ? :scientific : :auto)
    unscaled_ticks, labels
end

# return (continuous_values, discrete_values) for the ticks on this axis
function get_ticks(axis::Axis)
    ticks = _transform_ticks(axis[:ticks])
    ticks in (nothing, false) && return nothing

    dvals = axis[:discrete_values]
    cv, dv = if !isempty(dvals) && ticks == :auto
        # discrete ticks...
        axis[:continuous_values], dvals
    elseif ticks == :auto
        # compute optimal ticks and labels
        optimal_ticks_and_labels(axis)
    elseif typeof(ticks) <: Union{AVec, Int}
        # override ticks, but get the labels
        optimal_ticks_and_labels(axis, ticks)
    elseif typeof(ticks) <: NTuple{2, Any}
        # assuming we're passed (ticks, labels)
        ticks
    else
        error("Unknown ticks type in get_ticks: $(typeof(ticks))")
    end
    # @show ticks dvals cv dv

    # TODO: better/smarter cutoff values for sampling ticks
    if length(cv) > 30 && ticks == :auto
        rng = Int[round(Int,i) for i in linspace(1, length(cv), 15)]
        cv[rng], dv[rng]
    else
        cv, dv
    end
end


text:
using Base.Test
using Lora

filename = joinpath(dirname(@__FILE__), "sample.txt")

println("    Testing BasicVariableIOStream constructors and close method...")

iostream = BasicVariableIOStream((), 10000, filename)

close(iostream)
rm(filename)

println("    Testing BasicVariableIOStream IO methods...")

println("      Interaction with BasicUnvVariableState...")

nstatev = Float64[1.87, -4.5, 29.55, -0.91, 0.16]
iostreamsize = ()
iostreamn = length(nstatev)

iostream = BasicVariableIOStream(iostreamsize, iostreamn, filename, "w")
for v in nstatev
  write(iostream, BasicUnvVariableState(v))
end

close(iostream)

iostream = BasicVariableIOStream(iostreamsize, iostreamn, filename, "r")
nstate = read(iostream, Float64)

@test isa(nstate, BasicUnvVariableNState{Float64})
@test nstate.value == nstatev
@test nstate.n == iostream.n

close(iostream)
rm(filename)

println("      Interaction with BasicUnvVariableNState...")

nstatev = Float32[11.5, -41.22, -5.62, 1.98, 7.16]
iostreamsize = ()
iostreamn = length(nstatev)

iostream = BasicVariableIOStream(iostreamsize, iostreamn, filename)
nstatein = BasicUnvVariableNState(nstatev)
write(iostream, nstatein)

close(iostream)

iostream = BasicVariableIOStream(iostreamsize, iostreamn, filename, "r")
nstateout = read(iostream, Float32)

@test isa(nstateout, BasicUnvVariableNState{Float32})
@test nstateout.value == nstatein.value
@test nstateout.n == nstatein.n

close(iostream)
rm(filename)

println("      Interaction with BasicMuvVariableState...")

nstatev = Float64[8.11 -0.99 -4.19 0.1; 0.01 -0.02 1.4 8.47]
iostreamsize = (size(nstatev, 1),)
iostreamn = size(nstatev, 2)

iostream = BasicVariableIOStream(iostreamsize, iostreamn, filename)
for i in 1:iostreamn
  write(iostream, BasicMuvVariableState(nstatev[:, i]))
end

close(iostream)

iostream = BasicVariableIOStream(iostreamsize, iostreamn, filename, "r")
nstate = read(iostream, Float64)

@test isa(nstate, BasicMuvVariableNState{Float64})
@test nstate.value == nstatev
@test (nstate.size,) == iostream.size
@test nstate.n == iostream.n

close(iostream)
rm(filename)

println("      Interaction with BasicMuvVariableNState...")

nstatev = Float32[-7.1 -1.19 -7.76 6.1; -3.8 4.2 3.7 2.21]
iostreamsize = (size(nstatev, 1),)
iostreamn = size(nstatev, 2)

iostream = BasicVariableIOStream(iostreamsize, iostreamn, filename)
nstatein = BasicMuvVariableNState(nstatev)
write(iostream, nstatein)

close(iostream)

iostream = BasicVariableIOStream(iostreamsize, iostreamn, filename, "r")
nstateout = read(iostream, Float32)

@test isa(nstateout, BasicMuvVariableNState{Float32})
@test nstateout.value == nstatein.value
@test nstateout.size == nstatein.size
@test nstateout.n == nstatein.n

close(iostream)
rm(filename)

println("      Interaction with BasicMavVariableState...")

nstatev = Array(Float64, 3, 4, 2)
nstatev[:, :, 1] = [
  0.765895  0.965563  0.777571  0.508304;
  0.782592  0.366614  0.736418  0.863816;
  0.691046  0.524103  0.922673  0.223199
]
nstatev[:, :, 2] = [
  0.304872  0.210775  0.757637  0.657834;
  0.46678   0.995023  0.565492  0.007096;
  0.50774   0.362511  0.616699  0.675748
]
iostreamsize = (size(nstatev, 1), size(nstatev, 2))
iostreamn = size(nstatev, 3)

iostream = BasicVariableIOStream(iostreamsize, iostreamn, filename)
for i in 1:iostreamn
  write(iostream, BasicMavVariableState(nstatev[:, :, i]))
end

close(iostream)

iostream = BasicVariableIOStream(iostreamsize, iostreamn, filename, "r")
nstate = read(iostream, Float64)

@test isa(nstate, BasicMavVariableNState{Float64})
@test nstate.value == nstatev
@test nstate.size == iostream.size
@test nstate.n == iostream.n

close(iostream)
rm(filename)

println("      Interaction with BasicMavVariableNState...")

nstatev = Array(Float32, 3, 4, 2)
nstatev[:, :, 1] = [
  0.411455  0.268814  0.793351  0.502612;
  0.302095  0.702721  0.412092  0.488215;
  0.200365  0.862643  0.263403  0.442472
]
nstatev[:, :, 2] = [
 0.183705  0.0208819  0.423192  0.90061;
 0.31084   0.115514   0.464941  0.488014;
 0.913261  0.559753   0.969542  0.931011
]
iostreamsize = (size(nstatev, 1), size(nstatev, 2))
iostreamn = size(nstatev, 3)

iostream = BasicVariableIOStream(iostreamsize, iostreamn, filename)
nstatein = BasicMavVariableNState(nstatev)
write(iostream, nstatein)

close(iostream)

iostream = BasicVariableIOStream(iostreamsize, iostreamn, filename, "r")
nstateout = read(iostream, Float32)

@test isa(nstateout, BasicMavVariableNState{Float32})
@test nstateout.value == nstatein.value
@test nstateout.size == nstatein.size
@test nstateout.n == nstatein.n

close(iostream)
rm(filename)


text:
##########################
### MultitargetKNNClassifier
##########################

"""
    MultitargetKNNClassifier(;kwargs...)

$KNNClassifierDescription

$MultitargetKNNClassifierFields
"""
@mlj_model mutable struct MultitargetKNNClassifier <: MMI.Probabilistic
    K::Int = 5::(_ > 0)
    algorithm::Symbol = :kdtree::(_ in (:kdtree, :brutetree, :balltree))
    metric::Metric = Euclidean()
    leafsize::Int = 10::(_ ≥ 0)
    reorder::Bool = true
    weights::KNNKernel = Uniform()
    output_type::Type{<:MultiUnivariateFinite} = DictTable
end

function _predictmode_knnclassifier(weights, y, idxsvec)
    return StatsBase.mode.(_predict_knnclassifier(weights, y, idxsvec))
end

function dict_preds(::Val{:columnaccess}, func, target_table, idxsvec, weights)
    cols = Tables.columns(target_table)
    colnames = Tables.columnnames(cols)
    dict_table = Dict(
        nm => func(weights, Tables.getcolumn(cols, nm), idxsvec) for nm in colnames
    )
    return dict_table
end

function dict_preds(::Val{:noncolumnaccess}, func, target_table, idxsvec, weights)
     cols = Tables.dictcolumntable(target_table)
     colnames = Tables.columnnames(cols)
     dict_table = Dict(
        nm => func(weights, Tables.getcolumn(cols, nm), idxsvec) for nm in colnames
     )
     return dict_table
end

function dict_preds(func::F, target_table, idxsvec, weights) where {F<:Function}
    preds = if Tables.columnaccess(target_table)
        dict_preds(Val(:columnaccess), func, target_table, idxsvec, weights)
    else
        dict_preds(Val(:noncolumnaccess), func, target_table, idxsvec, weights)
    end
    return preds
end

function ntuple_preds(::Val{:columnaccess}, func, target_table, idxsvec, weights)
    cols = Tables.columns(target_table)
    colnames = Tables.columnnames(cols)
    column_table = NamedTuple{Tuple(colnames)}(
        func(weights, Tables.getcolumn(cols, nm), idxsvec) for nm in colnames
    )
    return column_table
end

function ntuple_preds(::Val{:noncolumnaccess}, func, target_table, idxsvec, weights)
    cols = Tables.columntable(target_table)
    colnames = Tables.columnnames(cols)
    column_table = NamedTuple{colnames}(
        func(weights, Tables.getcolumn(col, nm), idxsvec) for nm in colnames
    )
    return column_table
end

function ntuple_preds(func::F, target_table, idxsvec, weights) where {F <: Function}
    preds = if Tables.columnaccess(target_table)
        ntuple_preds(Val(:columnaccess), func, target_table, idxsvec, weights)
    else
        ntuple_preds(Val(:noncolumnaccess), func, target_table, idxsvec, weights)
    end
    return preds
end

function univariate_table(::Type{T}, weights, target_table, idxsvec) where {T}
    table = if T <: DictTable
        dict_preds(_predict_knnclassifier, target_table, idxsvec, weights)
    else
        ntuple_preds(_predict_knnclassifier, target_table, idxsvec, weights)
    end
    return table
end

function categorical_table(::Type{T}, weights, target_table, idxsvec) where {T}
    table = if T <: DictTable
        dict_preds(_predictmode_knnclassifier, target_table, idxsvec, weights)
    else
        ntuple_preds(_predictmode_knnclassifier, target_table, idxsvec, weights)
    end
    return table
end

function MMI.predict(m::MultitargetKNNClassifier, fitresult, X)
    err_if_given_invalid_K(m.K)
    Xmatrix = transpose(MMI.matrix(X))
    check_onebased_indexing("prediction input", Xmatrix)
    args = setup_predict_args(m, Xmatrix, fitresult)
    table_of_univariate_probs = univariate_table(m.output_type, args...)
    return table_of_univariate_probs
end

function MMI.predict_mode(m::MultitargetKNNClassifier, fitresult, X)
    err_if_given_invalid_K(m.K)
    Xmatrix = transpose(MMI.matrix(X))
    check_onebased_indexing("prediction input", Xmatrix)
    args = setup_predict_args(m, Xmatrix, fitresult)
    table_of_probs = categorical_table(m.output_type, args...)
    return table_of_probs
end

##########################
### KNNClassifier
##########################

"""
    KNNClassifier(;kwargs...)

$KNNClassifierDescription

$KNNFields
"""
@mlj_model mutable struct KNNClassifier <: MMI.Probabilistic
    K::Int = 5::(_ > 0)
    algorithm::Symbol = :kdtree::(_ in (:kdtree, :brutetree, :balltree))
    metric::Metric = Euclidean()
    leafsize::Int = 10::(_ ≥ 0)
    reorder::Bool = true
    weights::KNNKernel = Uniform()
end

function setup_predict_args(
    m::Union{KNNClassifier, MultitargetKNNClassifier}, Xmatrix, fitresult
)
    tree, targetvecortable, sample_weights = (
        fitresult.tree, fitresult.target, fitresult.sample_weights
    )
    kernel = m.weights
    # for each entry, get the K closest training point + their distance
    idxsvec, distsvec = NN.knn(tree, Xmatrix, m.K, sort_idxs(kernel))
    idxs_matrix = NeighIndsMatrix(idxsvec)
    dists_matrix = NeighDistsMatrix(distsvec)
    knn_weights = get_weights(kernel, dists_matrix)
    weights = combine_weights(idxs_matrix, knn_weights, sample_weights)
    return weights, targetvecortable, idxsvec
end

function _predict_knnclassifier(weights, y, idxsvec)
    # `classes` and `classes_seen` has an ordering consistent with the pool of y
    classes = @inbounds MMI.classes(y[1])
    classes_seen = filter(in(y), classes)
    nclasses = length(classes)
    nc = length(classes_seen) # Number of classes seen in `y`.

    # `labels` and `integers_seen` assigns ranks {1, 2, ..., nclasses}, to the
    # `CategoricalValue`'s in `y` according the ordering given by the pool of `y`.
    # In the case of `integers_seen` the ranks happen to sorted with the lower ranks
    # comming first and the higher rank coming last.
    # In both cases some ranks may be absent because some categorical values in `classes`
    # may be absent from `y` (in the case of `labels`) and `classes_seen`
    # (in the case of `integers_seen`).
    labels = MMI.int(y)
    integers_seen = MMI.int(classes_seen)

    # Recode `integers_seen` to be in {1,..., nc}
    # same as nc == nclasses || replace!(labels, (integers_seen .=> 1:nc)...)
    nc == nclasses || _replace!(labels, integers_seen, 1:nc)
    nsamples = length(idxsvec)
    probs = zeros(eltype(weights), nsamples, nc)

    @inbounds for i in eachindex(idxsvec)
        idxs = idxsvec[i]
        idxs_labels = @view(labels[idxs])
        @simd for j in eachindex(idxs_labels)
             @inbounds probs[i, idxs_labels[j]] += weights[i, j]
        end
    end
    return MMI.UnivariateFinite(classes_seen, scale(probs, dims=2))
end

function MMI.predict(m::KNNClassifier, fitresult, X)
    err_if_given_invalid_K(m.K)
    Xmatrix = transpose(MMI.matrix(X))
    check_onebased_indexing("prediction input", Xmatrix)
    predict_args = setup_predict_args(m, Xmatrix, fitresult)
    uni_var_probs = _predict_knnclassifier(predict_args...)
    return uni_var_probs
end

##########################
### KNNRegressor
##########################

"""
    KNNRegressoor(;kwargs...)

$KNNRegressorDescription

$KNNFields
"""
@mlj_model mutable struct KNNRegressor <: MMI.Deterministic
    K::Int = 5::(_ > 0)
    algorithm::Symbol = :kdtree::(_ in (:kdtree, :brutetree, :balltree))
    metric::Metric = Euclidean()
    leafsize::Int = 10::(_ ≥ 0)
    reorder::Bool = true
    weights::KNNKernel = Uniform()
end

function setup_predict_args(
    m::KNNRegressor, Xmatrix, fitresult
)
    tree, target, sample_weights = (
        fitresult.tree, fitresult.target, fitresult.sample_weights
    )
    kernel = m.weights
    # for each entry, get the K closest training point + their distance
    idxsvec, distsvec = NN.knn(tree, Xmatrix, m.K, sort_idxs(kernel))
    idxs_matrix = NeighIndsMatrix(idxsvec)
    dists_matrix = NeighDistsMatrix(distsvec)
    knn_weights = get_weights(kernel, dists_matrix)
    weights = combine_weights(idxs_matrix, knn_weights, sample_weights)
    return weights, target, idxs_matrix
end

function _predict_knnreg(weights, y, idxs_matrix)
    @inbounds labels = @view(y[idxs_matrix])
    preds = @inbounds(
        @view(_sum(labels .* weights, dims=2)[:]) ./ @view(_sum(weights, dims=2)[:])
    )
    return preds
end

function MMI.predict(m::KNNRegressor, fitresult, X)
    err_if_given_invalid_K(m.K)
    #Xmatrix = MMI.matrix(X, transpose=true) # NOTE: copies the data
    Xmatrix = transpose(MMI.matrix(X))
    check_onebased_indexing("prediction input", Xmatrix)
    predict_args = setup_predict_args(m, Xmatrix, fitresult)
    preds = _predict_knnreg(predict_args...)
    return preds
end

##########################
### MultitargetKNNRegressor
##########################

"""
    MultitargetKNNRegressor(;kwargs...)

$KNNRegressorDescription

$KNNFields
"""
@mlj_model mutable struct MultitargetKNNRegressor <: MMI.Deterministic
    K::Int = 5::(_ > 0)
    algorithm::Symbol = :kdtree::(_ in (:kdtree, :brutetree, :balltree))
    metric::Metric = Euclidean()
    leafsize::Int = 10::(_ ≥ 0)
    reorder::Bool = true
    weights::KNNKernel = Uniform()
end

function setup_predict_args(
    m::MultitargetKNNRegressor, Xmatrix, fitresult
)
    tree, names, target, sample_weights = (
        fitresult.tree,
        fitresult.target.names,
        fitresult.target.target,
        fitresult.sample_weights
    )
    kernel = m.weights
    # for each entry, get the K closest training point + their distance
    idxsvec, distsvec = NN.knn(tree, Xmatrix, m.K, sort_idxs(kernel))
    idxs_matrix = NeighIndsMatrix(idxsvec)
    dists_matrix = NeighDistsMatrix(distsvec)
    knn_weights = get_weights(kernel, dists_matrix)
    weights = combine_weights(idxs_matrix, knn_weights, sample_weights)
    return names, (weights, target, idxs_matrix)
end

function mul_label_by_weight(idxs_matrix, weights, Ymatrix, ind)
    i, j, k = ind.I
    @inbounds(label = Ymatrix[idxs_matrix[i, k], j])
    @inbounds(weight = weights[i, k])
    return label * weight
end

function _predict_multiknnreg(weights, Ymatrix, idxs_matrix)
    cartesian_inds = CartesianIndices(
        (axes(idxs_matrix, 1), axes(Ymatrix, 2), axes(idxs_matrix, 2))
    )
    weighted_labels = _sum(
        mul_label_by_weight.(
            Ref(idxs_matrix),
            Ref(weights),
            Ref(Ymatrix),
            cartesian_inds
        );
        dims = 3
   )
   # Normalizing the weighted labels gives us our final preds
   preds = @inbounds(@view(weighted_labels[:, :, 1])) ./ _sum(weights, dims=2)
   return preds
end

function MMI.predict(m::MultitargetKNNRegressor, fitresult, X)
    err_if_given_invalid_K(m.K)
    #Xmatrix = MMI.matrix(X, transpose=true) # NOTE: copies the data
    Xmatrix = transpose(MMI.matrix(X))
    check_onebased_indexing("prediction input", Xmatrix)
    names, args = setup_predict_args(m, Xmatrix, fitresult)
    preds = _predict_multiknnreg(args...)
    # A matrix table was intentionally outputed.
    # Users can coerce to their desired Table type.
    return MMI.table(preds, names = names)
end

##################
### FIT MODELS
##################

struct MultiKNNRegressorTarget{S, M<:AbstractMatrix}
    names::S
    target::M
end

_getcolnames(::Nothing) = nothing
_getcolnames(sch) = sch.names

function get_target(::MultitargetKNNRegressor, y)
    target = MMI.matrix(y)
    check_onebased_indexing("input_target", target)
    names = _getcolnames(MMI.schema(y))
    return MultiKNNRegressorTarget(names, target)
end

function get_target(::KNNRegressor, y)
   check_onebased_indexing("input_target", y)
   return y
end

# `CategoricalArrays` and `Table(Finite)` are one-based indexed arrays and tables.
get_target(::Union{KNNClassifier, MultitargetKNNClassifier}, y) = y

struct KNNResult{T, U, W}
    tree::T
    target::U
    sample_weights::W
end

const KNN = Union{
    KNNClassifier, KNNRegressor, MultitargetKNNRegressor, MultitargetKNNClassifier
}

function MMI.fit(m::KNN, verbosity::Int, X, y, w::Union{Nothing, Vec{<:Real}}=nothing)
    Xmatrix = transpose(MMI.matrix(X))
    target = get_target(m, y)
    if w !== nothing
        check_onebased_indexing("sample_weights", w)
        length(w) == size(Xmatrix, 2) || throw(
            DimensionMismatch(
                "Differing number of observations in sample_weights and target vector."
            )
        )
    end
    if m.algorithm == :kdtree
        tree = NN.KDTree(Xmatrix; leafsize=m.leafsize, reorder=m.reorder)
    elseif m.algorithm == :balltree
        tree = NN.BallTree(Xmatrix; leafsize=m.leafsize, reorder=m.reorder)
    elseif m.algorithm == :brutetree
        tree = NN.BruteTree(Xmatrix; leafsize=m.leafsize, reorder=m.reorder)
    end
    report = NamedTuple{}()
    return KNNResult(tree, target, w), nothing, report
end

MMI.fitted_params(model::KNN, fitresult) = (tree=fitresult.tree,)

########################
### Models Metadata
########################

metadata_model(
    KNNRegressor,
    input = Table(Continuous),
    target = Vec{Continuous},
    weights = true,
    descr = KNNRegressorDescription,
    path = "$(PKG).KNNRegressor"
)

metadata_model(
    KNNClassifier,
    input = Table(Continuous),
    target = Vec{<:Finite},
    weights = true,
    descr = KNNClassifierDescription,
    path = "$(PKG).KNNClassifier"
)

metadata_model(
    MultitargetKNNClassifier,
    input = Table(Continuous),
    target = Table(Finite),
    weights = true,
    descr = KNNClassifierDescription,
    path = "$(PKG).MultitargetKNNClassifier"
)

metadata_model(
    MultitargetKNNRegressor,
    input = Table(Continuous),
    target = Table(Continuous),
    weights = true,
    descr = KNNRegressorDescription,
    path = "$(PKG).MultitargetKNNRegressor"
)


text:
using Test
using LogicCircuits
using ProbabilisticCircuits
using DiscriminativeCircuits
using DataFrames: DataFrame
using CUDA

# This tests are supposed to test queries on the circuits
@testset "Logistic Circuit Query and Parameter Tests" begin
    # Uses a Logistic Circuit with 4 variables, and tests 3 of the configurations to
    # match with python version.

    CLASSES = 2

    logistic_circuit = zoo_dc("little_4var.circuit")
    @test logistic_circuit isa LogisticCircuit

    # check probabilities for binary samples
    data_raw = @. Bool([0 0 0 0; 0 1 1 0; 0 0 1 1])
    data = DataFrame(data_raw, :auto)
    # true_weight_func = [3.43147972 4.66740416;
    #                     4.27595352 2.83503504;
    #                     3.67415087 4.93793472]
    true_prob = [0.9686740008311808 0.9906908445371728;
                 0.9862917392724188 0.9445399509069984;
                 0.9752568185086389 0.9928816444223209]

    class_prob = class_likelihood_per_instance(logistic_circuit, data)
    for i = 1:size(true_prob)[1]
        for j = 1:CLASSES
            @test true_prob[i,j] ≈ class_prob[i,j]
        end
    end

    # check probabilities for float samples
    data = DataFrame(Float32.(data_raw), :auto)
    class_prob = class_likelihood_per_instance(logistic_circuit, data)
    for i = 1:size(true_prob)[1]
        for j = 1:CLASSES
            @test true_prob[i,j] ≈ class_prob[i,j]
        end
    end

    # check predicted_classes
    true_labels = [2, 1, 2]
    predicted_classes = predict_class(logistic_circuit, data)
    @test all(predicted_classes .== true_labels)

    # check accuracy
    @test accuracy(logistic_circuit, data, true_labels) == 1.0

    # check parameter updates
    original_literal_parameters = Dict{Int, Vector{Float64}}()
    foreach(logistic_circuit) do ln
        if ln isa Logistic⋁Node
            foreach(ln.children, eachrow(ln.thetas)) do c, theta
                if c isa LogisticLiteralNode
                    original_literal_parameters[c.literal] = copy(theta)
                end
            end
        end
    end

    one_hot_labels = [0.0 1.0;
                      1.0 0.0;
                      0.0 1.0]
    one_hot_labels = Float32.(one_hot_labels)
    true_error = true_prob .- one_hot_labels
    step_size = 0.1
    learn_parameters(logistic_circuit, data, true_labels; num_epochs=1, step_size=step_size)

    if CUDA.functional()
        data_gpu = to_gpu(data)
        labels_gpu = to_gpu(true_labels)
        learn_parameters(logistic_circuit, data_gpu, labels_gpu; num_epochs=1, step_size=step_size)
    end

    #####
    ### TODO; not valid code anymore `c.data`, not exactly sure what this was doing but seems was when we had upflow circuits
    ####
    # foreach(logistic_circuit) do ln
    #     if ln isa Logistic⋁Node
    #         foreach(ln.children, eachrow(ln.thetas)) do c, theta
    #             if c isa LogisticLiteralNode
    #                 for class = 1:CLASSES
    #                     true_update_amount = -step_size * sum(c.data.upflow .* true_error[:, class]) / size(true_error)[1]
    #                     updated_amount = theta[class] - original_literal_parameters[c.literal][class]
    #                     @test updated_amount ≈ true_update_amount atol=1e-7
    #                 end
    #             end
    #         end
    #     end
    # end

end

max_stars_repo_name: Juice-jl/DiscriminativeCircuits

max_issues_repo_name: Juice-jl/DiscriminativeCircuits

max_forks_repo_name: Juice-jl/DiscriminativeCircuits


text:
rig_key = "ocpi-lsk"
push!(RIGS, rig_key)

#Much like OCPI2, but with missing channels for missing hardware
const ocpi_lsk_mappings = Dict("AO0"=>"axial piezo",
#                      "AO1"=>"reserved",
                      "AO2"=>"analogout3",
                      "AO3"=>"analogout4",
                      "AI0"=>"axial piezo monitor",
#                      "AI1"=>"reserved",
                      "AI2"=>"analogin1",
                      "AI3"=>"analogin2",
                      "AI4"=>"analogin3",
                      "AI5"=>"analogin4",
                      "AI6"=>"analogin5",
                      "AI7"=>"analogin6",
                      "AI8"=>"analogin7",
                      "AI9"=>"analogin8",
                      "AI10"=>"analogin9",
                      "AI11"=>"analogin10",
                      "AI12"=>"analogin11",
                      "AI13"=>"analogin12",
                      "AI14"=>"analogin13",
                      "AI15"=>"analogin14",
                      "AI16"=>"analogin15",
                      "AI17"=>"analogin16",
                      "AI18"=>"analogin17",
                      "AI19"=>"analogin18",
                      "AI20"=>"analogin19",
                      "AI21"=>"analogin20",
                      "AI22"=>"analogin21",
                      "AI23"=>"analogin22",
                      "AI24"=>"analogin23",
                      "AI25"=>"analogin24",
                      "AI26"=>"analogin25",
                      "AI27"=>"analogin26",
                      "AI28"=>"analogin27",
                      "AI29"=>"analogin28",
                      "AI30"=>"analogin29",
                      "AI31"=>"analogin30",
                      "P0.0"=>"stimulus1",
                      "P0.1"=>"stimulus2",
                      "P0.2"=>"stimulus3",
                      "P0.3"=>"stimulus4",
                      "P0.4"=>"all lasers",
                      "P0.5"=>"camera1",
#                      "P0.6"=>"reserved",
#                      "P0.7"=>"reserved",
#                      "P0.8"=>"reserved",
#                      "P0.9"=>"reserved",
#                      "P0.10"=>"reserved",
#                      "P0.11"=>"reserved",
#                      "P0.12"=>"reserved",
                      "P0.13"=>"stimulus5",
                      "P0.14"=>"stimulus6",
                      "P0.15"=>"stimulus7",
                      "P0.16"=>"stimulus8",
                      "P0.17"=>"stimulus9",
                      "P0.18"=>"stimulus10",
                      "P0.19"=>"stimulus11",
                      "P0.20"=>"stimulus12",
                      "P0.21"=>"stimulus13",
                      "P0.22"=>"stimulus14",
                      "P0.23"=>"stimulus15",
                      "P0.24"=>"camera1 frame monitor",
#                      "P0.25"=>"reserved",
                      "P0.26"=>"diginput1",
                      "P0.27"=>"diginput2",
                      "P0.28"=>"diginput3",
                      "P0.29"=>"diginput4",
                      "P0.30"=>"diginput5",
                      "P0.31"=>"diginput6",)

DEFAULT_DAQCHANS_TO_NAMES[rig_key] = ocpi_lsk_mappings
DEFAULT_NAMES_TO_DAQCHANS[rig_key] = dictmap(reverse, ocpi_lsk_mappings)

const ocpi_lsk_pos_ctrl_chans = map(x->DEFAULT_NAMES_TO_DAQCHANS[rig_key][x], ["axial piezo"])
const ocpi_lsk_aochans = map(x->"AO$(x)", [0;2;3])
const ocpi_lsk_aichans = map(x->"AI$(x)", vcat([0;], [2:15...])) #currently AI1 is unused
const ocpi_lsk_dochans = map(x->"P0.$(x)", vcat([0:5...], [13:23...]))
const ocpi_lsk_dichans = map(x->"P0.$(x)", vcat([24;], [26:31...]))
const ocpi_lsk_pos_mon_chans = map(x->DEFAULT_NAMES_TO_DAQCHANS[rig_key][x], ["axial piezo monitor"])
const ocpi_lsk_cam_ctrl_chans = map(x->DEFAULT_NAMES_TO_DAQCHANS[rig_key][x], ["camera1"])
const ocpi_lsk_cam_mon_chans = map(x->DEFAULT_NAMES_TO_DAQCHANS[rig_key][x], ["camera1 frame monitor"])
const ocpi_lsk_laschans = map(x->DEFAULT_NAMES_TO_DAQCHANS[rig_key][x], ["all lasers"])
const ocpi_lsk_stimchans = map(x->DEFAULT_NAMES_TO_DAQCHANS[rig_key][x], ["stimulus$x" for x = 1:15])
const ocpi_lsk_fixed_names = ["axial piezo", "axial piezo monitor", "all lasers", "camera1", "camera1 frame monitor"]
const ocpi_lsk_fixed_daqchans = map(x->DEFAULT_NAMES_TO_DAQCHANS[rig_key][x], ocpi_lsk_fixed_names)

AO_CHANS[rig_key] = OrderedSet(ocpi_lsk_aochans)
AI_CHANS[rig_key] = OrderedSet(ocpi_lsk_aichans)
DO_CHANS[rig_key] = OrderedSet(ocpi_lsk_dochans)
DI_CHANS[rig_key] = OrderedSet(ocpi_lsk_dichans)
POS_CONTROL_CHANS[rig_key] = ocpi_lsk_pos_ctrl_chans
POS_MONITOR_CHANS[rig_key] = ocpi_lsk_pos_mon_chans
CAM_CONTROL_CHANS[rig_key] = ocpi_lsk_cam_ctrl_chans
CAM_MONITOR_CHANS[rig_key] = ocpi_lsk_cam_mon_chans
LAS_CONTROL_CHANS[rig_key] = ocpi_lsk_laschans
STIM_CHANS[rig_key] =OrderedSet(ocpi_lsk_stimchans)
GALVO_CONTROL_CHANS[rig_key] = String[]
GALVO_MONITOR_CHANS[rig_key] = String[]

FIXED_NAMES[rig_key] = ocpi_lsk_fixed_names
FIXED_DAQ_CHANS[rig_key] = ocpi_lsk_fixed_daqchans

#TODO: check these
#RIG_CHIP_SIZES[rig_key] = PCO_EDGE_5_5_CHIP_SIZE
#RIG_FRAMERATE_FUNCS[rig_key] = PCO_EDGE_5_5_FRAMERATE_FUNC
#RIG_ROI_FUNCS[rig_key] = PCO_EDGE_5_5_ROI_FUNC

#TODO: measure these
LASER_ON_TIME[rig_key] = 100 * Unitful.μs
LASER_OFF_TIME[rig_key] = 100 * Unitful.μs

CAMERA_ON_TIME[rig_key] = 19526.0 * Unitful.ns #This is the worst jitter measured with Edge 5.5 and 4.2 cameras, see ImagineInterface issue #18
CAMERA_OFF_TIME[rig_key] = 19526.0 * Unitful.ns #This is the worst jitter measured with Edge 5.5 and 4.2 cameras, see ImagineInterface issue #18

PIEZO_RANGES[rig_key] = (0.0μm .. 400.0μm, 0.0V .. 10.0V)
PIEZO_MAX_SPEED[rig_key] = 2000*Unitful.μm / Unitful.s
AO_RANGE[rig_key] = -10.0V .. 10.0V
AI_RANGE[rig_key] = AO_RANGE[rig_key] #TODO: make sure this is true.  (true if we are recording -10..10V on analog inputs)

max_forks_repo_name: HolyLab/ImagineInterface


text:
using Test
using MPI
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.MPIStateArrays

function main()
  FT = Float64
  Nhorz = 3
  Nstack = 5
  N = 4
  DA = Array

  MPI.Initialized() || MPI.Init()

  comm = MPI.COMM_WORLD
  crank = MPI.Comm_rank(comm)
  csize = MPI.Comm_size(comm)

  Rrange = FT.(accumulate(+,1:Nstack+1))
  topology = StackedCubedSphereTopology(MPI.COMM_SELF, Nhorz, Rrange;
                                        boundary=(1,2))
  grid = DiscontinuousSpectralElementGrid(topology; FloatType=FT,
                                          DeviceArray=DA, polynomialorder=N,
                                          meshwarp=Topologies.cubedshellwarp)

  #=
  @show elems       = topology.elems
  @show realelems   = topology.realelems
  @show ghostelems  = topology.ghostelems
  @show sendelems   = topology.sendelems
  @show elemtocoord = topology.elemtocoord
  @show elemtoelem  = topology.elemtoelem
  @show elemtoface  = topology.elemtoface
  @show elemtoordr  = topology.elemtoordr
  @show elemtobndy  = topology.elemtobndy
  @show nabrtorank  = topology.nabrtorank
  @show nabrtorecv  = topology.nabrtorecv
  @show nabrtosend  = topology.nabrtosend
  =#

  # Check x1x2x3 matches before comm
  x1 = @view grid.vgeo[:, Grids._x1, :]
  x2 = @view grid.vgeo[:, Grids._x2, :]
  x3 = @view grid.vgeo[:, Grids._x3, :]

  @test x1[grid.vmap⁻] ≈ x1[grid.vmap⁺]
  @test x2[grid.vmap⁻] ≈ x2[grid.vmap⁺]
  @test x3[grid.vmap⁻] ≈ x3[grid.vmap⁺]

  Np = (N+1)^3
  x1x2x3 = MPIStateArray{FT}(topology.mpicomm, DA, Np, 3,
                             length(topology.elems),
                             realelems=topology.realelems,
                             ghostelems=topology.ghostelems,
                             vmaprecv=grid.vmaprecv,
                             vmapsend=grid.vmapsend,
                             nabrtorank=topology.nabrtorank,
                             nabrtovmaprecv=grid.nabrtovmaprecv,
                             nabrtovmapsend=grid.nabrtovmapsend)
  x1x2x3.data[:,:,topology.realelems] .=
        @view grid.vgeo[:, [Grids._x1, Grids._x2, Grids._x3], topology.realelems]
  MPIStateArrays.start_ghost_exchange!(x1x2x3)
  MPIStateArrays.finish_ghost_exchange!(x1x2x3)

  # Check x1x2x3 matches after
  x1 = @view x1x2x3.data[:, 1, :]
  x2 = @view x1x2x3.data[:, 2, :]
  x3 = @view x1x2x3.data[:, 3, :]

  @test x1[grid.vmap⁻] ≈ x1[grid.vmap⁺]
  @test x2[grid.vmap⁻] ≈ x2[grid.vmap⁺]
  @test x3[grid.vmap⁻] ≈ x3[grid.vmap⁺]

  nothing
end
isinteractive() || main()

max_stars_repo_name: trontrytel/CLIMA

max_issues_repo_name: trontrytel/CLIMA

max_forks_repo_name: trontrytel/CLIMA


text:
function build_robot!(vis::Visualizer, model::Quadruped;
	r_body=0.02, r_thigh=0.0175, r_calf=0.0125, r_hip=0.02, r_knee=0.02, r_foot=0.02,
	color_opacity=1.0)

	body_mat = MeshPhongMaterial(color=RGBA(0.0, 0.0, 0.0, color_opacity))
	contact_mat = MeshPhongMaterial(color=RGBA(1.0, 165.0 / 255.0, 0.0, color_opacity))

	default_background!(vis)

	torso = GeometryBasics.Cylinder(Point3f0(0.0, 0.0, -model.l_torso2),
		Point3f0(0.0, 0.0, model.l_torso1),
		convert(Float32, r_body))
	setobject!(vis[:robot]["torso"], torso, body_mat)

	thigh_1 = GeometryBasics.Cylinder(Point3f0(0),
		Point3f0(0.0, 0.0, model.l_thigh1),
		convert(Float32, r_thigh))
	setobject!(vis[:robot]["thigh1"], thigh_1, body_mat)

	calf_1 = GeometryBasics.Cylinder(Point3f0(0),
		Point3f0(0.0, 0.0, model.l_calf1),
		convert(Float32, r_calf))
	setobject!(vis[:robot]["calf1"], calf_1, body_mat)

	thigh_2 = GeometryBasics.Cylinder(Point3f0(0),
		Point3f0(0.0, 0.0, model.l_thigh2),
		convert(Float32, r_thigh))
	setobject!(vis[:robot]["thigh2"], thigh_2, body_mat)

	calf_2 = GeometryBasics.Cylinder(Point3f0(0),
		Point3f0(0.0, 0.0, model.l_calf2),
		convert(Float32, r_calf))
	setobject!(vis[:robot]["calf2"], calf_2, body_mat)

	thigh_3 = GeometryBasics.Cylinder(Point3f0(0),
		Point3f0(0.0, 0.0, model.l_thigh3),
		convert(Float32, r_thigh))
	setobject!(vis[:robot]["thigh3"], thigh_3, body_mat)

	calf_3 = GeometryBasics.Cylinder(Point3f0(0),
		Point3f0(0.0, 0.0, model.l_calf3),
		convert(Float32, r_calf))
	setobject!(vis[:robot]["calf3"], calf_3, body_mat)

	thigh_4 = GeometryBasics.Cylinder(Point3f0(0),
		Point3f0(0.0, 0.0, model.l_thigh4),
		convert(Float32, r_thigh))
	setobject!(vis[:robot]["thigh4"], thigh_4, body_mat)

	calf_4 = GeometryBasics.Cylinder(Point3f0(0),
		Point3f0(0.0, 0.0, model.l_calf4),
		convert(Float32, r_calf))
	setobject!(vis[:robot]["calf4"], calf_4, body_mat)

	hip1 = setobject!(vis[:robot]["hip1"],
		GeometryBasics.Sphere(Point3f0(0), convert(Float32, r_hip)), body_mat)
	hip2 = setobject!(vis[:robot]["hip2"],
		GeometryBasics.Sphere(Point3f0(0), convert(Float32, r_hip)), body_mat)
	knee1 = setobject!(vis[:robot]["knee1"],
		GeometryBasics.Sphere(Point3f0(0), r_knee), body_mat)
	knee2 = setobject!(vis[:robot]["knee2"],
		GeometryBasics.Sphere(Point3f0(0), r_knee), body_mat)
	knee3 = setobject!(vis[:robot]["knee3"],
		GeometryBasics.Sphere(Point3f0(0), r_knee), body_mat)
	knee4 = setobject!(vis[:robot]["knee4"],
		GeometryBasics.Sphere(Point3f0(0), r_knee), body_mat)

	feet1 = setobject!(vis[:robot]["feet1"],
		GeometryBasics.Sphere(Point3f0(0), r_foot), contact_mat)
	feet2 = setobject!(vis[:robot]["feet2"],
		GeometryBasics.Sphere(Point3f0(0), r_foot), contact_mat)
	feet3 = setobject!(vis[:robot]["feet3"],
		GeometryBasics.Sphere(Point3f0(0), r_foot), contact_mat)
	feet4 = setobject!(vis[:robot]["feet4"],
		GeometryBasics.Sphere(Point3f0(0), r_foot), contact_mat)

	return true
end

function set_robot!(vis::Visualizer, model::Quadruped, q::AbstractVector;
		r_foot=0.02, y_offset=0.0)

	p_shift = [0.0; y_offset; r_foot]

	k_torso = q[1:2]
	p_torso = [k_torso[1], 0.0, k_torso[2]] + p_shift

	k_hip1 = kinematics_hip(model, q, hip=:hip1)
	p_hip1 = [k_hip1[1], 0.0, k_hip1[2]] + p_shift

	k_hip2 = kinematics_hip(model, q, hip=:hip2)
	p_hip2 = [k_hip2[1], 0.0, k_hip2[2]] + p_shift

	k_thigh_1 = kinematics_thigh(model, q, leg=:leg1, mode=:ee)
	p_thigh_1 = [k_thigh_1[1], 0.0, k_thigh_1[2]] + p_shift

	k_calf_1 = kinematics_calf(model, q, leg=:leg1, mode=:ee)
	p_calf_1 = [k_calf_1[1], 0.0, k_calf_1[2]] + p_shift

	k_thigh_2 = kinematics_thigh(model, q, leg=:leg2, mode=:ee)
	p_thigh_2 = [k_thigh_2[1], 0.0, k_thigh_2[2]] + p_shift

	k_calf_2 = kinematics_calf(model, q, leg=:leg2, mode=:ee)
	p_calf_2 = [k_calf_2[1], 0.0, k_calf_2[2]] + p_shift

	k_thigh_3 = kinematics_thigh(model, q, leg=:leg3, mode=:ee)
	p_thigh_3 = [k_thigh_3[1], 0.0, k_thigh_3[2]] + p_shift

	k_calf_3 = kinematics_calf(model, q, leg=:leg3, mode=:ee)
	p_calf_3 = [k_calf_3[1], 0.0, k_calf_3[2]] + p_shift

	k_thigh_4 = kinematics_thigh(model, q, leg=:leg4, mode=:ee)
	p_thigh_4 = [k_thigh_4[1], 0.0, k_thigh_4[2]] + p_shift

	k_calf_4 = kinematics_calf(model, q, leg=:leg4, mode=:ee)
	p_calf_4 = [k_calf_4[1], 0.0, k_calf_4[2]] + p_shift

	settransform!(vis[:robot]["thigh1"], cable_transform(p_hip1, p_thigh_1))
	settransform!(vis[:robot]["calf1"], cable_transform(p_thigh_1, p_calf_1))
	settransform!(vis[:robot]["thigh2"], cable_transform(p_hip1, p_thigh_2))
	settransform!(vis[:robot]["calf2"], cable_transform(p_thigh_2, p_calf_2))
	settransform!(vis[:robot]["thigh3"], cable_transform(p_hip2, p_thigh_3))
	settransform!(vis[:robot]["calf3"], cable_transform(p_thigh_3, p_calf_3))
	settransform!(vis[:robot]["thigh4"], cable_transform(p_hip2, p_thigh_4))
	settransform!(vis[:robot]["calf4"], cable_transform(p_thigh_4, p_calf_4))
	settransform!(vis[:robot]["torso"],
		compose(
			MeshCat.Translation(p_torso[1], p_torso[2], p_torso[3]),
			MeshCat.LinearMap(RotY(-q[3] + 0.5 * π))))
	settransform!(vis[:robot]["hip1"], MeshCat.Translation(p_hip1))
	settransform!(vis[:robot]["hip2"], MeshCat.Translation(p_hip2))
	settransform!(vis[:robot]["knee1"], MeshCat.Translation(p_thigh_1))
	settransform!(vis[:robot]["knee2"], MeshCat.Translation(p_thigh_2))
	settransform!(vis[:robot]["knee3"], MeshCat.Translation(p_thigh_3))
	settransform!(vis[:robot]["knee4"], MeshCat.Translation(p_thigh_4))
	settransform!(vis[:robot]["feet1"], MeshCat.Translation(p_calf_1))
	settransform!(vis[:robot]["feet2"], MeshCat.Translation(p_calf_2))
	settransform!(vis[:robot]["feet3"], MeshCat.Translation(p_calf_3))
	settransform!(vis[:robot]["feet4"], MeshCat.Translation(p_calf_4))

	return true
end

function visualize!(vis, model, q;
	Δt=0.1,
	r_body=0.02, r_thigh=0.0175, r_calf=0.0125, r_hip=0.02, r_knee=0.02, r_foot=0.02,
	color_opacity=1.0, fixed_camera=true)

	build_robot!(vis, model,
		r_body=r_body, r_thigh=r_thigh, r_calf=r_calf, r_hip=r_hip, r_knee=r_knee, r_foot=r_foot,
		color_opacity=color_opacity)

	anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

	for (t, qt) in enumerate(q)
		MeshCat.atframe(anim, t) do
			set_robot!(vis, model, qt)
		end
	end

	MeshCat.setanimation!(vis, anim)

	if fixed_camera
		settransform!(vis["/Cameras/default"],
        compose(Translation(0.0, -50.0, -1.0),LinearMap(RotZ(-pi / 2.0))))
        setprop!(vis["/Cameras/default/rotated/<object>"], "zoom", 25)
	end

	return true
end


text:
"""
    test_1layerqg_rossbywave(stepper, dt, nsteps, dev; kwargs...)

Evolves a Rossby wave and compares with the analytic solution.
"""
function test_1layerqg_rossbywave(stepper, dt, nsteps, dev::Device=CPU(); deformation_radius=Inf)
  nx = 64
  Lx = 2π
   β = 2.0
   μ = 0.0
   ν = 0.0
   T = Float64

  # the following if statement is called so that all the cases of
  # Problem() fuction are tested
  if stepper == "ForwardEuler"
    eta = zeros(dev, T, (nx, nx))
  else
    eta(x, y) = 0 * x
  end

  prob = SingleLayerQG.Problem(dev; nx=nx, Lx=Lx, eta=eta, deformation_radius=deformation_radius, β=β, μ=μ, ν=ν, stepper=stepper, dt=dt)
  sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid

  x, y = gridpoints(grid)

  # the Rossby wave initial condition
   ampl = 1e-2
  kwave = 3 * 2π/grid.Lx
  lwave = 2 * 2π/grid.Ly
      ω = -params.β * kwave / (kwave^2 + lwave^2 + 1 / deformation_radius^2)
    q₀  = @. ampl * cos(kwave * x) * cos(lwave * y)
    q₀h = rfft(q₀)

  SingleLayerQG.set_q!(prob, q₀)

  stepforward!(prob, nsteps)
  dealias!(sol, grid)
  SingleLayerQG.updatevars!(prob)

  q_theory = @. ampl * cos(kwave * (x - ω / kwave * clock.t)) * cos(lwave * y)

  return isapprox(q_theory, vars.q, rtol=grid.nx * grid.ny * nsteps * 1e-12)
end

"""
    test_1layerqg_stochasticforcing_energy_budgets(dev; kwargs...)

Tests if the energy budget is closed for a SingleLayerQG problem with stochastic forcing.
"""
function test_1layerqg_stochasticforcing_energybudget(dev::Device=CPU(); n=256, dt=0.01, L=2π, ν=1e-7, nν=2, μ=1e-1, T=Float64)
  n, L  = 256, 2π
  ν, nν = 1e-7, 2
  μ = 1e-1
  dt, tf = 0.005, 0.1/μ
  nt = round(Int, tf/dt)

  grid = TwoDGrid(dev, n, L)
  x, y = gridpoints(grid)
  K  = @. sqrt(grid.Krsq)                      # a 2D array with the total wavenumber

  # Forcing with spectrum ~ exp(-(k-kf)²/(2 dkf²))
  kf, dkf = 12.0, 2.0
  ε = 0.1

  forcing_spectrum = zeros(dev, T, (grid.nkr, grid.nl))
  @. forcing_spectrum = exp(-(K - kf)^2 / (2 * dkf^2))
  @. forcing_spectrum = ifelse(grid.Krsq < 2^2,  0, forcing_spectrum)   # no power at low wavenumbers
  @. forcing_spectrum = ifelse(grid.Krsq > 20^2, 0, forcing_spectrum)   # no power at high wavenumbers
  ε0 = parsevalsum(forcing_spectrum .* grid.invKrsq / 2, grid) / (grid.Lx * grid.Ly)
  @. forcing_spectrum = ε / ε0 * forcing_spectrum

  Random.seed!(1234)

  function calcF!(Fh, sol, t, clock, vars, params, grid)
    eta = ArrayType(dev)(exp.(2π * im * rand(T, size(sol))) / sqrt(clock.dt))
    CUDA.@allowscalar eta[1, 1] = 0
    @. Fh = eta * sqrt(forcing_spectrum)

    return nothing
  end

  prob = SingleLayerQG.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, μ=μ, dt=dt,
   stepper="RK4", calcF=calcF!, stochastic=true)

  SingleLayerQG.set_q!(prob, 0*x)

  E = Diagnostic(SingleLayerQG.kinetic_energy,     prob, nsteps=nt)
  D = Diagnostic(SingleLayerQG.energy_dissipation, prob, nsteps=nt)
  R = Diagnostic(SingleLayerQG.energy_drag,        prob, nsteps=nt)
  W = Diagnostic(SingleLayerQG.energy_work,        prob, nsteps=nt)
  diags = [E, D, W, R]

  stepforward!(prob, diags, nt)

  SingleLayerQG.updatevars!(prob)

  dEdt_numerical = (E[2:E.i] - E[1:E.i-1]) / prob.clock.dt

  dEdt_computed = W[2:E.i] - D[1:E.i-1] - R[1:E.i-1]

  return isapprox(dEdt_numerical, dEdt_computed, rtol=1e-3)
end

"""
    test_1layerqg_deterministicforcing_energy_budgets(dev; kwargs...)

Tests if the energy budget is closed for a SingleLayerQG problem with deterministic forcing.
"""
function test_1layerqg_deterministicforcing_energybudget(dev::Device=CPU(); n=256, dt=0.01, L=2π, ν=1e-7, nν=2, μ=1e-1)
  n, L  = 256, 2π
  ν, nν = 1e-7, 2
  μ = 1e-1
  dt, tf = 0.005, 0.1/μ
  nt = round(Int, tf/dt)

  grid = TwoDGrid(dev, n, L)
  x, y = gridpoints(grid)
  k₀, l₀ = 2π / grid.Lx, 2π/grid.Ly

  # Forcing = 0.01cos(4x)cos(5y)cos(2t)
  f = @. 0.01 * cos(4k₀*x) * cos(5l₀*y)
  fh = rfft(f)

  function calcF!(Fh, sol, t, clock, vars, params, grid)
    @. Fh = fh * cos(2t)

    return nothing
  end

  prob = SingleLayerQG.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, μ=μ, dt=dt,
   stepper="RK4", calcF=calcF!, stochastic=false)

  SingleLayerQG.set_q!(prob, 0*x)

  E = Diagnostic(SingleLayerQG.kinetic_energy,     prob, nsteps=nt)
  D = Diagnostic(SingleLayerQG.energy_dissipation, prob, nsteps=nt)
  R = Diagnostic(SingleLayerQG.energy_drag,        prob, nsteps=nt)
  W = Diagnostic(SingleLayerQG.energy_work,        prob, nsteps=nt)
  diags = [E, D, W, R]

  stepforward!(prob, diags, round(Int, nt))

  SingleLayerQG.updatevars!(prob)

  dEdt_numerical = (E[3:E.i] - E[1:E.i-2]) / (2 * prob.clock.dt)
  dEdt_computed  = W[2:E.i-1] - D[2:E.i-1] - R[2:E.i-1]

  return isapprox(dEdt_numerical, dEdt_computed, rtol=1e-4)
end

"""
    test_1layerqg_stochasticforcing_enstrophy_budgets(dev; kwargs...)

Tests if the enstrophy budget is closed for a SingleLayerQG problem with stochastic forcing.
"""
function test_1layerqg_stochasticforcing_enstrophybudget(dev::Device=CPU(); n=256, dt=0.01, L=2π, ν=1e-7, nν=2, μ=1e-1, T=Float64)
  n, L  = 256, 2π
  ν, nν = 1e-7, 2
  μ = 1e-1
  dt, tf = 0.005, 0.1/μ
  nt = round(Int, tf/dt)

  grid = TwoDGrid(dev, n, L)
  x, y = gridpoints(grid)
  K  = @. sqrt(grid.Krsq)                      # a 2D array with the total wavenumber

  # Forcing with spectrum ~ exp(-(k-kf)²/(2 dkf²))
  kf, dkf = 12.0, 2.0
  εᶻ = 0.1

  forcing_spectrum = zeros(dev, T, (grid.nkr, grid.nl))
  @. forcing_spectrum = exp(-(K - kf)^2 / (2 * dkf^2))
  @. forcing_spectrum = ifelse(grid.Krsq < 2^2,  0, forcing_spectrum)   # no power at low wavenumbers
  @. forcing_spectrum = ifelse(grid.Krsq > 20^2, 0, forcing_spectrum)   # no power at high wavenumbers
  εᶻ0 = parsevalsum(forcing_spectrum / 2, grid) / (grid.Lx * grid.Ly)
  forcing_spectrum .= εᶻ / εᶻ0 * forcing_spectrum

  Random.seed!(1234)

  function calcF!(Fh, sol, t, clock, vars, params, grid)
    eta = ArrayType(dev)(exp.(2π * im * rand(T, size(sol))) / sqrt(clock.dt))
    CUDA.@allowscalar eta[1, 1] = 0
    @. Fh = eta * sqrt(forcing_spectrum)

    return nothing
  end

  prob = SingleLayerQG.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, μ=μ, dt=dt,
   stepper="RK4", calcF=calcF!, stochastic=true)

  SingleLayerQG.set_q!(prob, 0*x)

  Z = Diagnostic(SingleLayerQG.enstrophy,             prob, nsteps=nt)
  D = Diagnostic(SingleLayerQG.enstrophy_dissipation, prob, nsteps=nt)
  R = Diagnostic(SingleLayerQG.enstrophy_drag,        prob, nsteps=nt)
  W = Diagnostic(SingleLayerQG.enstrophy_work,        prob, nsteps=nt)
  diags = [Z, D, W, R]

  stepforward!(prob, diags, nt)

  SingleLayerQG.updatevars!(prob)

  dZdt_numerical = (Z[2:Z.i] - Z[1:Z.i-1]) / prob.clock.dt

  dZdt_computed = W[2:Z.i] - D[1:Z.i-1] - R[1:Z.i-1]

  return isapprox(dZdt_numerical, dZdt_computed, rtol=1e-3)
end

"""
    test_1layerqg_deterministicforcing_enstrophy_budgets(dev; kwargs...)

Tests if the enstrophy budget is closed for a SingleLayerQG problem with deterministic forcing.
"""
function test_1layerqg_deterministicforcing_enstrophybudget(dev::Device=CPU(); n=256, dt=0.01, L=2π, ν=1e-7, nν=2, μ=1e-1)
  n, L  = 256, 2π
  ν, nν = 1e-7, 2
  μ = 1e-1
  dt, tf = 0.005, 0.1/μ
  nt = round(Int, tf/dt)

  grid = TwoDGrid(dev, n, L)
  x, y = gridpoints(grid)
  k₀, l₀ = 2π/grid.Lx, 2π/grid.Ly

  # Forcing = 0.01cos(4x)cos(5y)cos(2t)
  f = @. 0.01 * cos(4k₀ * x) * cos(5l₀ * y)
  fh = rfft(f)

  function calcF!(Fh, sol, t, clock, vars, params, grid)
    @. Fh = fh * cos(2t)
    return nothing
  end

  prob = SingleLayerQG.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, μ=μ, dt=dt,
   stepper="RK4", calcF=calcF!, stochastic=false)

  SingleLayerQG.set_q!(prob, 0*x)

  Z = Diagnostic(SingleLayerQG.enstrophy,             prob, nsteps=nt)
  D = Diagnostic(SingleLayerQG.enstrophy_dissipation, prob, nsteps=nt)
  R = Diagnostic(SingleLayerQG.enstrophy_drag,        prob, nsteps=nt)
  W = Diagnostic(SingleLayerQG.enstrophy_work,        prob, nsteps=nt)
  diags = [Z, D, W, R]

  stepforward!(prob, diags, round(Int, nt))

  SingleLayerQG.updatevars!(prob)

  dZdt_numerical = (Z[3:Z.i] - Z[1:Z.i-2]) / (2 * prob.clock.dt)
  dZdt_computed  = W[2:Z.i-1] - D[2:Z.i-1] - R[2:Z.i-1]

  return isapprox(dZdt_numerical, dZdt_computed, rtol=1e-4)
end

"""
    test_1layerqg_nonlinearadvection(dt, stepper, dev; kwargs...)

Tests the advection term in the SingleLayerQG module by timestepping a
test problem with timestep dt and timestepper identified by the string stepper.
The test problem is derived by picking a solution ζf (with associated
streamfunction ψf) for which the advection term J(ψf, ζf) is non-zero. Next, a
forcing Ff is derived according to Ff = ∂ζf/∂t + J(ψf, ζf) - ν∇²ζf. One solution
to the vorticity equation forced by this Ff is then ζf. (This solution may not
be realized, at least at long times, if it is unstable.)
"""
function test_1layerqg_advection(dt, stepper, dev::Device=CPU(); n=128, L=2π, ν=1e-2, nν=1, μ=0.0)
  n, L  = 128, 2π
  ν, nν = 1e-2, 1
   μ = 0.0
  tf = 1.0
  nt = round(Int, tf/dt)

  grid = TwoDGrid(dev, n, L)
  x, y = gridpoints(grid)

  ψf = @. sin(2x) * cos(2y) + 2sin(x) * cos(3y)
  qf = @. -8sin(2x) * cos(2y) - 20sin(x) * cos(3y)

  Ff = @. -(
    ν*( 64sin(2x) * cos(2y) + 200sin(x) * cos(3y) )
    + 8 * ( cos(x) * cos(3y) * sin(2x) * sin(2y) - 3cos(2x) * cos(2y) * sin(x) * sin(3y) )
    )

  Ffh = rfft(Ff)

  function calcF!(Fh, sol, t, clock, vars, params, grid)
    Fh .= Ffh
    return nothing
  end

  prob = SingleLayerQG.Problem(dev; nx=n, Lx=L, ν=ν, nν=nν, μ=μ, dt=dt, stepper=stepper, calcF=calcF!)

  SingleLayerQG.set_q!(prob, qf)

  stepforward!(prob, round(Int, nt))

  SingleLayerQG.updatevars!(prob)

  return isapprox(prob.vars.q, qf, rtol=rtol_singlelayerqg)
end

"""
    test_1layerqg_energyenstrophy_BarotropicQG(dev)

Tests the energy and enstrophy function for a SingleLayerQG problem.
"""
function test_1layerqg_energyenstrophy_BarotropicQG(dev::Device=CPU())
  nx, Lx  = 64, 2π
  ny, Ly  = 64, 3π
  grid = TwoDGrid(dev, nx, Lx, ny, Ly)
  k₀, l₀ = 2π/Lx, 2π/Ly # fundamental wavenumbers
  x, y = gridpoints(grid)

  energy_calc = 29/9
  enstrophy_calc = 10885/648

  η  = @. cos(10k₀ * x) * cos(10l₀ * y)
  ψ₀ = @. sin(2k₀ * x) * cos(2l₀ * y) + 2sin(k₀ * x) * cos(3l₀ * y)
  q₀ = @. - ((2k₀)^2 + (2l₀)^2) * sin(2k₀ * x) * cos(2l₀ * y) - (k₀^2 + (3l₀)^2) * 2sin(k₀ * x) * cos(3l₀*y)

  prob = SingleLayerQG.Problem(dev; nx=nx, Lx=Lx, ny=ny, Ly=Ly, eta=η, stepper="ForwardEuler")
  sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
  SingleLayerQG.set_q!(prob, q₀)
  SingleLayerQG.updatevars!(prob)

  energyq₀ = SingleLayerQG.energy(prob)
  enstrophyq₀ = SingleLayerQG.enstrophy(prob)

  return isapprox(energyq₀, energy_calc, rtol=rtol_singlelayerqg) && isapprox(enstrophyq₀, enstrophy_calc, rtol=rtol_singlelayerqg) && SingleLayerQG.potential_energy(prob)==0 &&
  SingleLayerQG.addforcing!(prob.timestepper.N, sol, clock.t, clock, vars, params, grid) == nothing
end

"""
    test_1layerqg_energies_EquivalentBarotropicQG(dev)

Tests the kinetic and potential energy for an equivalent barotropic SingleLayerQG problem.
"""
function test_1layerqg_energies_EquivalentBarotropicQG(dev; deformation_radius=1.23)
  nx, Lx  = 64, 2π
  ny, Ly  = 64, 3π
  grid = TwoDGrid(dev, nx, Lx, ny, Ly)
  k₀, l₀ = 2π/Lx, 2π/Ly # fundamental wavenumbers
  x, y = gridpoints(grid)

  kinetic_energy_calc = 29/9
  potential_energy_calc = 5/(8*deformation_radius^2)
  energy_calc = kinetic_energy_calc + potential_energy_calc

  η  = @. cos(10k₀ * x) * cos(10l₀ * y)
  ψ₀ = @. sin(2k₀ * x) * cos(2l₀ * y) + 2sin(k₀ * x) * cos(3l₀ * y)
  q₀ = @. - ((2k₀)^2 + (2l₀)^2) * sin(2k₀ * x) * cos(2l₀ * y) - (k₀^2 + (3l₀)^2) * 2sin(k₀ * x) * cos(3l₀*y) - 1/deformation_radius^2 * ψ₀

  prob = SingleLayerQG.Problem(dev; nx=nx, Lx=Lx, ny=ny, Ly=Ly, eta=η, deformation_radius=deformation_radius, stepper="ForwardEuler")
  sol, clock, vars, params, grid = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
  SingleLayerQG.set_q!(prob, q₀)
  SingleLayerQG.updatevars!(prob)

  kinetic_energyq₀ = SingleLayerQG.kinetic_energy(prob)
  potential_energyq₀ = SingleLayerQG.potential_energy(prob)
  energyq₀ = SingleLayerQG.energy(prob)

  return isapprox(kinetic_energyq₀, kinetic_energy_calc, rtol=rtol_singlelayerqg) && isapprox(potential_energyq₀, potential_energy_calc, rtol=rtol_singlelayerqg) && isapprox(energyq₀, energy_calc, rtol=rtol_singlelayerqg) &&
  SingleLayerQG.addforcing!(prob.timestepper.N, sol, clock.t, clock, vars, params, grid) == nothing
end

"""
    test_1layerqg_problemtype(dev, T; deformation_radius=Inf)

Tests the SingleLayerQG problem constructor for different DataType `T`.
"""
function test_1layerqg_problemtype(dev, T; deformation_radius=Inf)
  prob = SingleLayerQG.Problem(dev; T=T, deformation_radius=deformation_radius)

  A = ArrayType(dev)

  return (typeof(prob.sol)<:A{Complex{T}, 2} && typeof(prob.grid.Lx)==T && eltype(prob.grid.x)==T && typeof(prob.vars.u)<:A{T, 2})
end

function test_streamfunctionfrompv(dev; deformation_radius=1.23)
  prob_barotropicQG = SingleLayerQG.Problem(dev; nx=64, deformation_radius=Inf)
  prob_equivalentbarotropicQG = SingleLayerQG.Problem(dev; nx=64, deformation_radius=deformation_radius)

  grid = prob_barotropicQG.grid
  k₀, l₀ = 2π/grid.Lx, 2π/grid.Ly # fundamental wavenumbers
  x, y = gridpoints(grid)

  ψ = @. sin(2k₀ * x) * cos(3l₀ * y) + 0*sin(3k₀ * x)

  q_barotropic = @. -((2k₀)^2 + (3l₀)^2) * sin(2k₀ * x) * cos(3l₀ * y) - 0*(3k₀)^2 * sin(3k₀ * x)
  q_equivalentbarotropic = @. -((2k₀)^2 + (3l₀)^2 + 1/deformation_radius^2) * sin(2k₀ * x) * cos(3l₀ * y) - 0*((3k₀)^2+ 1/deformation_radius^2) * sin(3k₀ * x)

  SingleLayerQG.set_q!(prob_barotropicQG, q_barotropic)
  SingleLayerQG.set_q!(prob_equivalentbarotropicQG, q_equivalentbarotropic)

  SingleLayerQG.streamfunctionfrompv!(prob_barotropicQG.vars.ψh, prob_barotropicQG.vars.qh, prob_barotropicQG.params, prob_barotropicQG.grid)

  SingleLayerQG.streamfunctionfrompv!(prob_equivalentbarotropicQG.vars.ψh, prob_equivalentbarotropicQG.vars.qh, prob_equivalentbarotropicQG.params, prob_equivalentbarotropicQG.grid)

  return (prob_barotropicQG.vars.ψ ≈ ψ && prob_equivalentbarotropicQG.vars.ψ ≈ ψ)
end

function test_1layerqg_energy_dissipation(dev; deformation_radius=2.23)
  prob = SingleLayerQG.Problem(dev; nx=16, deformation_radius=deformation_radius)
  SingleLayerQG.energy_dissipation(prob)
  return nothing
end

function test_1layerqg_enstrophy_dissipation(dev; deformation_radius=2.23)
  prob = SingleLayerQG.Problem(dev; nx=16, deformation_radius=deformation_radius)
  SingleLayerQG.enstrophy_dissipation(prob)
  return nothing
end

function test_1layerqg_energy_work(dev; deformation_radius=2.23)
  prob = SingleLayerQG.Problem(dev; nx=16, deformation_radius=deformation_radius)
  SingleLayerQG.energy_work(prob)
  return nothing
end

function test_1layerqg_enstrophy_work(dev; deformation_radius=2.23)
  prob = SingleLayerQG.Problem(dev; nx=16, deformation_radius=deformation_radius)
  SingleLayerQG.enstrophy_work(prob)
  return nothing
end

function test_1layerqg_energy_drag(dev; deformation_radius=2.23)
  prob = SingleLayerQG.Problem(dev; nx=16, deformation_radius=deformation_radius)
  SingleLayerQG.energy_drag(prob)
  return nothing
end

function test_1layerqg_enstrophy_drag(dev; deformation_radius=2.23)
  prob = SingleLayerQG.Problem(dev; nx=16, deformation_radius=deformation_radius)
  SingleLayerQG.enstrophy_drag(prob)
  return nothing
end


text:
export append_barrier!, device_barrier

append_barrier!(list::ZeCommandList, signal_event=nothing, wait_events::ZeEvent...) =
    zeCommandListAppendBarrier(list, something(signal_event, C_NULL),
                               length(wait_events), [wait_events...])

device_barrier(dev::ZeDevice) = zeDeviceSystemBarrier(dev)


text:
@testset "4.4.1.2 (d csc)^m (a+b cot)^n" begin
    (a, b, n, x, ) = @variables a b n x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*of*the*form*(d*csc(e+f*x))^m*(a+b*cot(e+f*x))^n=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(d*csc(e+f*x))^m*(a+b*cot(e+f*x))^n*with*a^2+b^2=0=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*csc(e+f*x)^m*(a+a*im*cot(e+f*x))^n=#

    #= ::Subsubsection:: =#
    #=n>0=#

    #= ::Subsubsection::Closed:: =#
    #=n<0=#

    @test_int [sin(x)^4/(im + cot(x)), x, 4, -((5*im*x)/16) + 1/(32*(im - cot(x))^2) - im/(8*(im - cot(x))) - im/(24*(im + cot(x))^3) - 3/(32*(im + cot(x))^2) + (3*im)/(16*(im + cot(x)))]
    @test_int [sin(x)^3/(im + cot(x)), x, 3, (4/5)*im*cos(x) - (4/15)*im*cos(x)^3 + (im*sin(x)^3)/(5*(im + cot(x)))]
    @test_int [sin(x)^2/(im + cot(x)), x, 4, -((3*im*x)/8) - im/(8*(im - cot(x))) - 1/(8*(im + cot(x))^2) + im/(4*(im + cot(x)))]
    @test_int [sin(x)^1/(im + cot(x)), x, 2, (2/3)*im*cos(x) + (im*sin(x))/(3*(im + cot(x)))]
    @test_int [csc(x)^1/(im + cot(x)), x, 1, (im*csc(x))/(im + cot(x))]
    @test_int [csc(x)^2/(im + cot(x)), x, 2, (-im)*x + log(sin(x))]
    @test_int [csc(x)^3/(im + cot(x)), x, 2, im*atanh(cos(x)) - csc(x)]
    @test_int [csc(x)^4/(im + cot(x)), x, 2, im*cot(x) - cot(x)^2/2]
    @test_int [csc(x)^5/(im + cot(x)), x, 3, (1/2)*im*atanh(cos(x)) + (1/2)*im*cot(x)*csc(x) - csc(x)^3/3]
    @test_int [csc(x)^6/(im + cot(x)), x, 3, im*cot(x) - cot(x)^2/2 + (1/3)*im*cot(x)^3 - cot(x)^4/4]
    @test_int [csc(x)^7/(im + cot(x)), x, 4, (3/8)*im*atanh(cos(x)) + (3/8)*im*cot(x)*csc(x) + (1/4)*im*cot(x)*csc(x)^3 - csc(x)^5/5]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(d*csc(e+f*x))^m*(a+b*cot(e+f*x))^n=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*csc(e+f*x)^m*(a+b*cot(e+f*x))^n=#

    #= ::Subsubsection:: =#
    #=n>0=#

    #= ::Subsubsection::Closed:: =#
    #=n<0=#

    @test_int [csc(x)^6/(a + b*cot(x)), x, 3, (a*(a^2 + 2*b^2)*cot(x))/b^4 - ((a^2 + 2*b^2)*cot(x)^2)/(2*b^3) + (a*cot(x)^3)/(3*b^2) - cot(x)^4/(4*b) - ((a^2 + b^2)^2*log(a + b*cot(x)))/b^5]
    @test_int [csc(x)^4/(a + b*cot(x)), x, 3, (a*cot(x))/b^2 - cot(x)^2/(2*b) - ((a^2 + b^2)*log(a + b*cot(x)))/b^3]
    @test_int [csc(x)^2/(a + b*cot(x)), x, 2, -(log(a + b*cot(x))/b)]
    @test_int [sin(x)^2/(a + b*cot(x)), x, 7, (a*(a^2 + 3*b^2)*x)/(2*(a^2 + b^2)^2) - (b^3*log(b*cos(x) + a*sin(x)))/(a^2 + b^2)^2 - ((b + a*cot(x))*sin(x)^2)/(2*(a^2 + b^2))]
    @test_int [sin(x)^4/(a + b*cot(x)), x, 8, (a*(3*a^4 + 10*a^2*b^2 + 15*b^4)*x)/(8*(a^2 + b^2)^3) - (b^5*log(b*cos(x) + a*sin(x)))/(a^2 + b^2)^3 - ((4*b^3 + a*(3*a^2 + 7*b^2)*cot(x))*sin(x)^2)/(8*(a^2 + b^2)^2) - ((b + a*cot(x))*sin(x)^4)/(4*(a^2 + b^2))]

    @test_int [csc(x)^5/(a + b*cot(x)), x, 9, (a*atanh(cos(x)))/(2*b^2) + (a*(a^2 + b^2)*atanh(cos(x)))/b^4 + ((a^2 + b^2)^(3/2)*atanh(((b - a*cot(x))*sin(x))/sqrt(a^2 + b^2)))/b^4 - ((a^2 + b^2)*csc(x))/b^3 + (a*cot(x)*csc(x))/(2*b^2) - csc(x)^3/(3*b)]
    @test_int [csc(x)^3/(a + b*cot(x)), x, 5, (a*atanh(cos(x)))/b^2 + (sqrt(a^2 + b^2)*atanh(((b - a*cot(x))*sin(x))/sqrt(a^2 + b^2)))/b^2 - csc(x)/b]
    @test_int [csc(x)^1/(a + b*cot(x)), x, 2, -(atanh(((-b + a*cot(x))*sin(x))/sqrt(a^2 + b^2))/sqrt(a^2 + b^2))]
    @test_int [sin(x)^1/(a + b*cot(x)), x, 5, (b^2*atanh(((b - a*cot(x))*sin(x))/sqrt(a^2 + b^2)))/(a^2 + b^2)^(3/2) - (a*cos(x))/(a^2 + b^2) - (b*sin(x))/(a^2 + b^2)]
    @test_int [sin(x)^3/(a + b*cot(x)), x, 9, (b^4*atanh(((b - a*cot(x))*sin(x))/sqrt(a^2 + b^2)))/(a^2 + b^2)^(5/2) - (a*b^2*cos(x))/(a^2 + b^2)^2 - (a*cos(x))/(a^2 + b^2) + (a*cos(x)^3)/(3*(a^2 + b^2)) - (b^3*sin(x))/(a^2 + b^2)^2 - (b*sin(x)^3)/(3*(a^2 + b^2))]

    @test_int [csc(x)^2/(a + b*cot(x))^2, x, 2, 1/(b*(a + b*cot(x)))]

    #= ::Subsection:: =#
    #=Integrands*of*the*form*csc(e+f*x)^m*(a+b*cot(e+f*x))^(n/2)=#

    #= ::Subsection:: =#
    #=Integrands*of*the*form*csc(e+f*x)^(m/2)*(a+b*cot(e+f*x))^n=#

    #= ::Subsection:: =#
    #=Integrands*of*the*form*csc(e+f*x)^(m/2)*(a+b*cot(e+f*x))^(n/2)=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(d*csc(e+f*x))^m*(a+b*cot(e+f*x))^n*with*n*symbolic=#
    #==#

    @test_int [csc(x)^2*(a + b*cot(x))^n, x, 2, -((a + b*cot(x))^(1 + n)/(b*(1 + n)))]

    #= ::Subsection:: =#
    #=Integrands*of*the*form*sin(e+f*x)^(m/2)*(a+b*cot(e+f*x))^n=#

    #= ::Subsection:: =#
    #=Integrands*of*the*form*sin(e+f*x)^(m/2)*(a+b*cot(e+f*x))^(n/2)=#
    end

max_stars_repo_path: test/4 Trig functions/4.4 Cotangent/4.4.1.2 (d csc)^m (a+b cot)^n.jl


text:
# ---
# title: 1578. Minimum Deletion Cost to Avoid Repeating Letters
# id: problem1578
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Greedy
# link: <https://leetcode.com/problems/minimum-deletion-cost-to-avoid-repeating-letters/description/>
# hidden: true
# ---
#
# Given a string `s` and an array of integers `cost` where `cost[i]` is the cost
# of deleting the `ith` character in `s`.
#
# Return the minimum cost of deletions such that there are no two identical
# letters next to each other.
#
# Notice that you will delete the chosen characters at the same time, in other
# words, after deleting a character, the costs of deleting other characters will
# not change.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "abaac", cost = [1,2,3,4,5]
#     Output: 3
#     Explanation: Delete the letter "a" with cost 3 to get "abac" (String without two identical letters next to each other).
#
#
# **Example 2:**
#
#
#
#     Input: s = "abc", cost = [1,2,3]
#     Output: 0
#     Explanation: You don't need to delete any character because there are no identical letters next to each other.
#
#
# **Example 3:**
#
#
#
#     Input: s = "aabaa", cost = [1,2,3,4,1]
#     Output: 2
#     Explanation: Delete the first and the last character, getting the string ("aba").
#
#
#
#
# **Constraints:**
#
#   * `s.length == cost.length`
#   * `1 <= s.length, cost.length <= 10^5`
#   * `1 <= cost[i] <= 10^4`
#   * `s` contains only lowercase English letters.
#
#
## @lc code=start
using LeetCode

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


text:
# This file was generated, do not modify it. # hide
#hideall
using Markdown

mdC_code = Markdown.htmlesc(raw"""#include <stdio.h>
#include <stdint.h>
#include <time.h>

uint8_t cnt(double ReZ,double ImZ,double ReC,double ImC){
    uint8_t k=0;
    while(k < 255){
        double new_ReZ = ReZ * ReZ - ImZ * ImZ + ReC ;
        double new_ImZ = 2 * ReZ * ImZ + ImC;
        if (new_ReZ * new_ReZ+new_ImZ * new_ImZ>4){
            return k;
        }
        ReZ = new_ReZ;
        ImZ = new_ImZ;
        k+=1;
    }
    return k;
}

int main()
{
    clock_t start,end;
    start = clock();
    int M = 2500;
    int N = 2500;
    uint8_t grid[M][N];
    double init_Re = 0.0;
    double init_Im = 0.0;
    for (int j = 0; j < M; j++)
    {
        for (int i = 0; i < N; i++)
        {
            double ReC = 4.0 / (N - 1) * i - 2.0;
            double ImC = 4.0 / (M - 1) * j - 2.0;
            double ReZ = init_Re;
            double ImZ = init_Im;
            grid[j][i] = cnt(ReZ,ImZ,ReC,ImC);
        }
    }
    end = clock();
    double cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC;
    printf("Elapsed time %f [sec]\n", cpu_time_used);
    FILE *fp;
    fp = fopen("/tmp/benchC.txt", "w");
    for (int j = 0; j < M; j++)
    {
        for (int i = 0; i < N; i++)
        {
            fprintf(fp, "%d", grid[j][i]);
            if (i != N - 1)
            {
                fprintf(fp, ",");
            }
            else
            {
                fprintf(fp, "\n");
            }
        }
    }
    fclose(fp);
    return 0;
}""")
mdfile = joinpath(dirname(@OUTPUT), "mandC.md")
open(mdfile,"w") do f
    print(f, mdC_code)
end

C_code=raw"""
#include <stdio.h>
#include <stdint.h>
#include <time.h>

uint8_t cnt(double ReZ,double ImZ,double ReC,double ImC){
    uint8_t k=0;
    while(k < 255){
        double new_ReZ = ReZ * ReZ - ImZ * ImZ + ReC ;
        double new_ImZ = 2 * ReZ * ImZ + ImC;
        if (new_ReZ * new_ReZ+new_ImZ * new_ImZ>4){
            return k;
        }
        ReZ = new_ReZ;
        ImZ = new_ImZ;
        k+=1;
    }
    return k;
}

int main()
{
    clock_t start,end;
    start = clock();
    int M = 2500;
    int N = 2500;
    uint8_t grid[M][N];
    double init_Re = 0.0;
    double init_Im = 0.0;
    for (int j = 0; j < M; j++)
    {
        for (int i = 0; i < N; i++)
        {
            double ReC = 4.0 / (N - 1) * i - 2.0;
            double ImC = 4.0 / (M - 1) * j - 2.0;
            double ReZ = init_Re;
            double ImZ = init_Im;
            grid[j][i] = cnt(ReZ,ImZ,ReC,ImC);
        }
    }
    end = clock();
    double cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC;
    printf("Elapsed time %f [sec]\n", cpu_time_used);
    FILE *fp;
    fp = fopen("/tmp/benchC.txt", "w");
    for (int j = 0; j < M; j++)
    {
        for (int i = 0; i < N; i++)
        {
            fprintf(fp, "%d", grid[j][i]);
            if (i != N - 1)
            {
                fprintf(fp, ",");
            }
            else
            {
                fprintf(fp, "\n");
            }
        }
    }
    fclose(fp);
    return 0;
}
"""

exefile = tempname()

#=
This trick is taken from

https://discourse.julialang.org/t/how-to-make-a-c-function-compiled-by-myself-available-to-ccall/7972/26
=#

open(`gcc -Wall -O2 -march=native -xc -o $exefile -`, "w") do f
    print(f, C_code)
end

run(`$exefile`)