text:
###################################################
# Random Variable Type Hierarchy
###################################################
#=
Arguably, much of this is duplicating insupport/support from distributions,
but is more usable for our purposes of mapping between constrained and
unconstrained model parameters.
=#
abstract RVType

abstract RScalar <: RVType
abstract RVector <: RVType
abstract RMatrix <: RVType

ndims(x::RScalar) = 1
nfree(x::RScalar) = 1

storage_type(rv::RScalar, T) = T
storage_type(rv::RVector, T) = Vector{T}
storage_type(rv::RMatrix, T) = Matrix{T}
###################################################
# random variable concrete types
###################################################
"""
Real-valued scalar.
"""
immutable RReal  <: RScalar
end
constrain(rv::RReal, x::Real) = x
unconstrain(rv::RReal, x::Real) = x
logdetjac(rv::RReal, x::Real) = 0.

"""
Positively constrained variable with (optional) lower bound.
"""
immutable RPositive{T <: Real}  <: RScalar
    lb::T
end
RPositive{T <: Real}(lb::T=0.) = RPositive{T}(lb)
constrain(rv::RPositive, x::Real) = exp(x) + rv.lb
unconstrain(rv::RPositive, x::Real) = log(x - rv.lb)
logdetjac(rv::RPositive, x::Real) = x

"""
Constrained variable with (optional) upper bound.
"""
immutable RNegative{T <: Real}  <: RScalar
    ub::T
end
RNegative{T <: Real}(ub::T=0.) = RNegative{T}(ub)
constrain(rv::RNegative, x::Real) = rv.ub - exp(x)
unconstrain(rv::RNegative, x::Real) = log(rv.ub - x)
logdetjac(rv::RNegative, x::Real) = x

"""
Constrained variable with (optional) upper and lower bound.
Using logistic transformation.
"""
immutable RBounded{T <: Real}  <: RScalar
    lb::T
    ub::T
end

RBounded{T <: Real}(lb::T = 0., ub::T = 1.0) = RBounded{T}(lb, ub)
constrain(rv::RBounded, x::Real) = rv.lb + (rv.ub - rv.lb) * StatsFuns.logistic(x)
unconstrain(rv::RBounded, x::Real) = logit((x - rv.lb) / (rv.ub - rv.lb))
logdetjac(rv::RBounded, x::Real) = log(rv.ub - rv.lb) - x - 2 * StatsFuns.log1pexp(-x)

"""
Probability constrained value.
Using logistic transformation.
"""
immutable RProbability <: RScalar
end

constrain(rv::RProbability, x::Real) = StatsFuns.logistic(x)
unconstrain(rv::RProbability, x::Real) = StatsFuns.logit(x)
logdetjac(rv::RProbability, x::Real) = - x - 2log1pexp(-x)

"""
Correlation constrained value.
Using tanh(x) function.
"""
immutable RCorrelation <: RScalar
end

constrain(rv::RCorrelation, x::Real) = tanh(x)
unconstrain(rv::RCorrelation, x::Real) = atanh(x)
logdetjac(rv::RCorrelation, x::Real) = log(4) + 2x - 2 * StatsFuns.log1pexp(2x)

"""
Real-valued vector.
"""
immutable RRealVec  <: RVector
    d::Int
end
ndims(x::RRealVec) = x.d
nfree(x::RRealVec) = ndims(x)
constrain(::RRealVec, x::Vector) = x
unconstrain(::RRealVec, x::Vector) = x
logdetjac(::RRealVec, x::Vector) = 0.

"""
Unit length vector.
"""
immutable RUnitVec  <: RVector
    d::Int
end
ndims(x::RUnitVec) = x.d
nfree(x::RUnitVec) = ndims(x)
constrain(::RUnitVec, x::Vector) = x ./ sqrt(dot(x, x))
unconstrain(::RUnitVec, x::Vector, len::Real) = len * x  # need to provide the norm of original vector to free x
logdetjac(::RUnitVec, x::Vector) = - .5 * dot(x, x)

"""
Ordered Constraint.
"""
immutable ROrdered  <: RVector
    d::Int
end
ndims(x::ROrdered) = x.d
nfree(x::ROrdered) = ndims(x)
function constrain(rv::ROrdered, x::Vector)
    y = zeros(length(x))
    y[:] = x
    # Or use y = copy(x); drawback needed!
    for j in 1:(ndims(rv) - 1)
        y[j + 1] = y[j] + exp(x[j + 1]) # make sure the latter is greater than previous.
    end
    y
end

function unconstrain(rv::ROrdered, x::Vector)
    y = copy(x)
    for j in 1:(ndims(rv) - 1)
        y[j + 1] = log(x[j + 1] - x[j])
    end
    y
end

logdetjac(::ROrdered, x::Vector) = sum(x[2:end])

"""
Positive Ordered Constraint
"""
immutable RPosOrdered  <: RVector
    d::Int
end
ndims(x::RPosOrdered) = x.d
nfree(x::RPosOrdered) = ndims(x)
function constrain(rv::RPosOrdered, x::Vector)
    y = copy(x)
    y[1] = exp(x[1]) # Ensure the first is positive
    for j in 1:(ndims(rv) - 1)
        y[j + 1] = y[j] + exp(x[j + 1])
    end
    y
end

function unconstrain(rv::RPosOrdered, x::Vector)
    y = copy(x)
    for j in 1:(ndims(rv) - 1)
        y[j + 1] = log(x[j + 1] - x[j])
    end
    y[1] = log(x[1])
    y
end

logdetjac(::RPosOrdered, x::Vector) = sum(x)

"""
Simplex Constraint (entries sum to one)
"""
immutable RSimplex  <: RVector
    d::Int
end
ndims(x::RSimplex) = x.d
nfree(x::RSimplex) = ndims(x)
function constrain(rv::RSimplex, x::Vector)
    y = copy(x)
    stick_len = 1 # entries sum to 1
    for j in 1:ndims(rv)
        y[j] = stick_len * logistic(x[j] - log(ndims(rv) - (j - 1)))
        stick_len -= y[j]
    end
    y! = [y; stick_len]
end

function unconstrain(rv::RSimplex, x::Vector)
    y = zeros(ndims(rv))
    stick_len = 1
    for j in 1:ndims(rv)
        y[j] = logit(x[j] / stick_len) + log(ndims(rv) - (j - 1))
        stick_len -= x[j]
    end
    y
end

function logdetjac(rv::RSimplex, x::Vector)
    y = copy(x)
    d = length(x)
    stick_len = 1
    logdetx = 0
    for j in 1:ndims(rv)
        constan = x[j] - log(d - (j - 1))
        y[j] = stick_len * logistic(constan)
        logdetx += log(stick_len) - StatsFuns.log1pexp(-constan) - StatsFuns.log1pexp(constan)
        stick_len -= y[j]
    end
    logdetx
end

"""
Random Cholesky factor (lower triangular matrix with positive diagonal).
"""
immutable RCholFact <: RMatrix
    d::Int
end

ndims(x::RCholFact) = x.d
nfree(x::RCholFact) = (p = ndims(x); p * (p + 1) ÷ 2)

function constrain(rv::RCholFact, x::Vector)
    U = UpperTriangular(x)
    for j in 1:ndims(rv)
        U[j, j] = exp(U[j, j])  # Cholesky factor must have positive diagonals
    end
    transpose(U)
end

function unconstrain(rv::RCholFact, S::LowerTriangular)
    U = transpose(S)
    for j in 1:ndims(rv)
        U[j, j] = log(U[j, j])  # diagonal of Cholesky is positive, so take log
    end
    flatten(U)
end

logdetjac(rv::RCholFact, x::Vector) = sum(diag(UpperTriangular(x)))

"""
Cholesky Correlation Matrix (lower triangular matrix with elements between -1 and 1).
"""
immutable RCholCorr <: RMatrix
    d::Int
end

ndims(x::RCholCorr) = x.d
nfree(x::RCholCorr) = (p = ndims(x); p * (p - 1) ÷ 2)

function constrain(rv::RCholCorr, x::Vector)
    z = tanh(x)
    pos = 1
    U = eye(ndims(rv))
    for j in 1:ndims(rv)
        for i in 1:j-1
            sumsqs = 0 # faster than dot(vec(U.data[1:i, j]), vec(U.data[1:i, j]))
            for ij in 1:i
                sumsqs += U[ij, j]^2
            end
            U[i, j] = z[pos] * sqrt(1 - sumsqs)
            pos += 1
        end
        sumsqs = 0 # compute dot(U[1:j-1, j], U[1:j-1, j]) in loop
        for ij in 1:j-1
            sumsqs += U[ij, j]^2
        end
        U[j, j] = sqrt(1 - sumsqs)
    end
    transpose(UpperTriangular(U)) # faster than LowerTriangular(transpose(U))
end

function unconstrain(rv::RCholCorr, S::LowerTriangular)
    U = transpose(S)
    pos = 1
    z = zeros(Int(.5 * ndims(rv) * (ndims(rv) - 1)))
    for j in 2:ndims(rv)
        for i in 1:j-1
            sumsqs = 0  # faster than dot(vec(U.data[1:i-1, j]), vec(U.data[1:i-1, j]))
            for ij in 1:i-1
                sumsqs += U[ij, j]^2
            end
            z[pos] = U[i, j] / sqrt(1 - sumsqs)
            pos += 1
        end
    end
    atanh(z)
end

function logdetjac(rv::RCholCorr, x::Vector)
    z = tanh(x)
    pos = 2 # z[1] has no contribution to logdetjac
    U = eye(ndims(rv))
    logdetL = sum(log(4) + 2x - 2 * StatsFuns.log1pexp(2x))
    for j in 3:ndims(rv)
        for i in 1:j-1
            sumsqs = 0
            for ij in 1:i
                sumsqs += U[ij, j]^2
            end
            logdetL += 0.5 * log1p(-sumsqs)
            U[i, j] = z[pos] * sqrt(1 - sumsqs)
            pos += 1
        end
        sumsqs = 0 # compute dot(U[1:j-1, j], U[1:j-1, j]) in loop
        for ij in 1:j-1
            sumsqs += U[ij, j]^2
        end
        U[j, j] = sqrt(1 - sumsqs)
    end
    logdetL
end

"""
Random Correlation matrix (symmetric, positive-definite).
"""
immutable RCorrMat <: RMatrix
    d::Int
end

ndims(x::RCorrMat) = x.d
nfree(x::RCorrMat) = (p = ndims(x); p * (p - 1) ÷ 2)

function constrain(rv::RCorrMat, x::Vector)
    z = tanh(x)
    pos = 1
    L = eye(ndims(rv))
    for j in 1:ndims(rv)
        for i in j+1:ndims(rv)
            sumsqs = 0 # faster than dot(vec(L[i, 1:j]), vec(L[i, 1:j]))
            for ij in 1:j
                sumsqs += L[i, ij]^2
            end
            L[i, j] = z[pos] * sqrt(1 - sumsqs)
            pos += 1
        end
        sumsqs = 0
        for ij in 1:j-1
            sumsqs += L[j, ij]^2
        end
        L[j, j] = sqrt(1 - sumsqs)
    end
    PDMat(Base.LinAlg.Cholesky(full(L), :L))
end

function unconstrain(::RCorrMat, S::PDMat)
    L = copy(S.chol[:L])
    K = dim(S)
    z = zeros(Int(.5 * K * (K - 1)))
    z[1:K-1] = L.data[2:K, 1]
    pos = K
    for j in 2:K
        for i in j+1:K
            sumsqs = 0
            for ij in 1:j-1
                sumsqs += L[i, ij]^2
            end
            z[pos] = L[i, j] / sqrt(1 - sumsqs)
            pos += 1
        end
    end
    atanh(z)
end

function logdetjac(rv::RCorrMat, x::Vector)
    z = tanh(x)
    pos = 1
    val = zeros(length(x))
    logdetL = length(x) * log(4) + 2sum(x) - 2sum(StatsFuns.log1pexp(2x))
    for k in 1:ndims(rv)-2
        for i in k+1:ndims(rv)
            val[pos] = (ndims(rv) - k - 1) * log1p(-z[pos]^2) # from Eq (11) of LKJ paper
            pos += 1
        end
    end
    logdetL += .5 * sum(val)
end

"""
Random covariance matrix (symmetric, positive-definite).
"""
immutable RCovMat <: RMatrix
    d::Int
end

ndims(x::RCovMat) = x.d
nfree(x::RCovMat) = (p = ndims(x); p * (p + 1) ÷ 2)

function constrain(rv::RCovMat, x::Vector)
    U = UpperTriangular(x)
    for j in 1:ndims(rv)
        U[j, j] = exp(U[j, j])  # Cholesky factor must have positive diagonals
    end
    PDMat(Base.LinAlg.Cholesky(full(U), :U))
end

function unconstrain(::RCovMat, S::PDMat)
    U = copy(S.chol[:U])
    for j in 1:dim(S)
        U[j, j] = log(U[j, j])  # diagonal of Cholesky is positive, so take log
    end
    flatten(U)
end

function logdetjac(rv::RCovMat, x::Vector)
    d = ndims(rv)
    d * logtwo + (d+1:-1:2) ⋅ diag(LowerTriangular(x))
end

"""
Random LKJ covariance matrix (symmetric, positive-definite).
"""
immutable RCovLKJ <: RMatrix
    d::Int
end

ndims(x::RCovLKJ) = x.d
nfree(x::RCovLKJ) = (p = ndims(x); p * (p + 1) ÷ 2)

function constrain(rv::RCovLKJ, x::Vector)
    d = ndims(rv)
    K = Int(.5d * (d - 1))  # if K not an integer, should throw an error message
    z = tanh(x[1:K])
    diagz = exp(x[K + 1:end])
    pos = 1
    L = eye(ndims(rv))
    for j in 1:ndims(rv)
        for i in j+1:ndims(rv)
            sumsqs = 0
            for ij in 1:j
                sumsqs += L[i, ij]^2
            end
            #sumsqs = dot(vec(L[i, 1:j]), vec(L[i, 1:j]))
            L[i, j] = z[pos] * sqrt(1 - sumsqs)
            pos += 1
        end
        sumsqs = 0
        for ij in 1:j-1
            sumsqs += L[j, ij]^2
        end
        L[j, j] = sqrt(1 - sumsqs)
    end
    PDMat(Base.LinAlg.Cholesky(full(Base.LinAlg.Diagonal(diagz) * L), :L))
end

function unconstrain(::RCovLKJ, S::PDMat)
    K = dim(S)
    z = zeros(Int(.5K * (K + 1)))
    z[Int(.5K * (K - 1))+1:end] = .5log(diag(S)) # identify the σ^2 variances
    L = S.chol[:L] ./ sqrt(diag(S))
    z[1:K-1] = L[2:K, 1]
    pos = K
    for j in 2:K
        for i in j+1:K
            sumsqs = 0
            for ij in 1:j-1
                sumsqs += L[i, ij]^2
            end
            z[pos] = L[i, j] / sqrt(1 - sumsqs)
            pos += 1
        end
    end
    z[1:Int(.5K * (K - 1))] = atanh(z[1:Int(.5K * (K - 1))])
    z
end

function logdetjac(rv::RCovLKJ, x::Vector)
    d = ndims(rv)
    K = Int(.5d * (d - 1))
    z = tanh(x[1:K])
    pos = 1
    val = zeros(length(x))
    logdetL = K * log(4) + 2sum(x[1:K]) - 2sum(StatsFuns.log1pexp(2x[1:K]))
    logdetL += sum(x[K+1:end]) # log determinant for positive constraint
    logdetL += d .* (sum(x[K+1:end]) + log(2)) # correction from correlation to covariance
    for k in 1:d-2
        for i in k+1:d
            val[pos] = (d - k - 1) * log1p(-z[pos]^2)
            pos += 1
        end
    end
    logdetL += .5 * sum(val) # log determinant for correlation constraint
end

constrain(pars, d::Distribution) = map(constrain, parsupp(d), pars)
unconstrain(d::Distribution) = map(unconstrain, parsupp(d), params(d))

"""
Number of free parameters needed for (multivariate) normal approximation
to the posterior over unconstrained parameters in ADVI.
"""
function num_pars_advi(rv::RVType, full=false)
    p = nfree(rv)
    full ? p * (p + 3) ÷ 2 : 2p
end


text:
using Merly
using Test
using HTTP
using JSON
#Pkg.add("BenchmarkTools")
#using BenchmarkTools

ip = "127.0.0.1"
port = 8086

@test File("index.html") == read( open(joinpath(pwd(),"index.html")), String)

useCORS(true)
@test useCORS(true) == true
headersalways("Strict-Transport-Security","max-age=10886400; includeSubDomains; preload")
notfound("""<!DOCTYPE html>

              <title>Not found</title>

404, Not found

              """)
notfound("cosa/notfound.html")

u = "hello"

@page "/" "Hello World!"
@page "/hola/:usr>" "Hello {{usr}}!"

@page "/mifile" File("Index.html")

@route GET "/get/:data1>" begin
  "get this back: {{data1}}"
end

@route GET "/regex/(\\w+\\d+)" begin

  println("req.version ",req.version)
  println("req.headers ",req.headers)

   "datos $(req.params[1])"
end

@route POST "/post" begin
  res.body = "I did something!"
  "I did something!"
end

@route POST|PUT|DELETE "/" begin
  println("params: ",req.params)
  println("query: ",req.query)
  println("body: ",req.body)

  res.headers["Content-Type"]= "text/plain"

  "I did something!"
end

Get("/data", (req,res)->(begin
  res.headers["Content-Type"]= "text/plain"
  println("params: ",req.params)
  println("query: ",req.query)
  u*"data"
end))

Post("/data", (req,res)->(begin
  println("params: ",req.params)
  println("query: ",req.query)
  println("body: ",req.body)
  res.headers["Content-Type"]= "text/plain"
  global u="bye"
  "I did something!"
end))

@test webserverpath("cosa") == joinpath(pwd(),"cosa")
webserverfiles("*")

 @async start(host = ip, port = port, verbose = false)
 sleep(2)

 r = HTTP.get("http://$(ip):$(port)/")
 @test r.status == 200
 @test String(r.body) == "Hello World!"

r = HTTP.get("http://$(ip):$(port)/hola/usuario")
@test r.status == 200
@test String(r.body) == "Hello usuario!"

r = HTTP.get("http://$(ip):$(port)/get/testdata")
@test r.status == 200
@test String(r.body) == "get this back: testdata"

r= HTTP.get("http://$(ip):$(port)/regex/test1")
@test String(r.body) == "datos test1"

r = HTTP.get("http://$(ip):$(port)/data?hola=1")
@test r.status == 200
@test String(r.body) == "hellodata"

myjson = Dict("query"=>"data")
my_headers = HTTP.mkheaders(["Accept" => "application/json","Content-Type" => "application/xml"])
r = HTTP.post("http://$(ip):$(port)/data",my_headers,JSON.json(myjson))
@test r.status == 200
@test String(r.body) == "I did something!"

r = HTTP.get("http://$(ip):$(port)/data")
@test r.status == 200
@test String(r.body) == "byedata"

my_headers = HTTP.mkheaders(["Accept" => "application/json"])
r = HTTP.post("http://$(ip):$(port)/post",my_headers,JSON.json(myjson))
@test r.status == 200
@test String(r.body) == "I did something!"

my_headers = HTTP.mkheaders(["Accept" => "application/json","Content-Type" => "application/json"])
r = HTTP.post("http://$(ip):$(port)/",my_headers,JSON.json(myjson))
@test r.status == 200
@test String(r.body) == "I did something!"

r = HTTP.put("http://$(ip):$(port)/",my_headers,JSON.json(myjson))
@test r.status == 200
@test String(r.body) == "I did something!"

r = HTTP.delete("http://$(ip):$(port)/")
@test r.status == 200
@test String(r.body) == "I did something!"
@test r.headers == Pair{SubString{String},SubString{String}}["Content-Type"=>"text/plain", "Access-Control-Allow-Origin"=>"*", "Access-Control-Allow-Methods"=>"POST,GET,OPTIONS", "Strict-Transport-Security"=>"max-age=10886400; includeSubDomains; preload", "Transfer-Encoding"=>"chunked"]

try
  r = HTTP.get("http://$(ip):$(port)/nada")
catch e
  @test String(e.response.body) == File("notfound.html")
end

r = HTTP.get("http://$(ip):$(port)/prueba.txt")
@test r.status == 200
@test String(r.body) == "probando webserver"
@test r.headers == Pair{SubString{String},SubString{String}}["Content-Type"=>"text/plain", "Access-Control-Allow-Origin"=>"*", "Access-Control-Allow-Methods"=>"POST,GET,OPTIONS", "Strict-Transport-Security"=>"max-age=10886400; includeSubDomains; preload", "Transfer-Encoding"=>"chunked"]

r= HTTP.get("http://$(ip):$(port)/algomas/ja.txt")
@test String(r.body) == "jajajajaj"

#@btime HTTP.get("http://$(ip):$(port)/?hola=5")
# 3.864 ms (8304 allocations: 381.20 KiB)
# 453.803 μs (748 allocations: 30.92 KiB)
#@btime HTTP.get("http://$(ip):$(port)/hola/usuario")
# 4.211 ms (7685 allocations: 353.44 KiB)
# 483.861 μs (743 allocations: 30.41 KiB)
#@benchmark r= HTTP.get("http://$(ip):$(port)/get/testdata")
# 3.906 ms (7693 allocations: 353.78 KiB)
# 484.960 μs (744 allocations: 30.44 KiB)
#=
 minimum time:     740.087 μs (0.00% GC)
  median time:      843.457 μs (0.00% GC)
  mean time:        1.099 ms (1.45% GC)
  maximum time:     367.611 ms (0.00% GC)
=#

ext: jl


text:
# metropolis_intro.jl

include("../src/intro.jl")
include("../src/fast_LP.jl")
include("new_LP_test.jl")

max_stars_repo_name: sebmart/NYC-timing

max_issues_repo_name: sebmart/NYC-timing

max_forks_repo_name: sebmart/NYC-timing


text:
export BrickletDualButtonLEDState
struct BrickletDualButtonLEDState
    led_l::Integer
    led_r::Integer
end

export BrickletDualButtonButtonState
struct BrickletDualButtonButtonState
    button_l::Integer
    button_r::Integer
end

export BrickletDualButtonIdentity
struct BrickletDualButtonIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletDualButton
"""
Two tactile buttons with built-in blue LEDs
"""
mutable struct BrickletDualButton <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletDualButton(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_dual_button")
        deviceInternal = package.BrickletDualButton(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export set_led_state
"""
    $(SIGNATURES)

Sets the state of the LEDs. Possible states are:

* 0 = AutoToggleOn: Enables auto toggle with initially enabled LED.
* 1 = AutoToggleOff: Activates auto toggle with initially disabled LED.
* 2 = On: Enables LED (auto toggle is disabled).
* 3 = Off: Disables LED (auto toggle is disabled).

In auto toggle mode the LED is toggled automatically at each press of a button.

If you just want to set one of the LEDs and don't know the current state
of the other LED, you can get the state with :func:`Get LED State` or you
can use :func:`Set Selected LED State`.
"""
function set_led_state(device::BrickletDualButton, led_l, led_r)
    device.deviceInternal.set_led_state(led_l, led_r)
end

export get_led_state
"""
    $(SIGNATURES)

Returns the current state of the LEDs, as set by :func:`Set LED State`.
"""
function get_led_state(device::BrickletDualButton)
    return device.deviceInternal.get_led_state()
end

export get_button_state
"""
    $(SIGNATURES)

Returns the current state for both buttons. Possible states are:

* 0 = pressed
* 1 = released
"""
function get_button_state(device::BrickletDualButton)
    return device.deviceInternal.get_button_state()
end

export set_selected_led_state
"""
    $(SIGNATURES)

Sets the state of the selected LED (0 or 1).

The other LED remains untouched.
"""
function set_selected_led_state(device::BrickletDualButton, led, state)
    device.deviceInternal.set_selected_led_state(led, state)
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletDualButton)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletDualButton, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
function init_constraint!(
    com::CS.CoM,
    constraint::XNorConstraint,
    fct,
    set::XNorSet;
)
    !init_lhs_and_rhs!(com, constraint, fct, set) && return false

    if constraint.complement_lhs !== nothing
        init_constraint!(com, constraint.complement_lhs, constraint.complement_lhs.fct, constraint.complement_lhs.set)
    end
    if constraint.complement_rhs !== nothing
        init_constraint!(com, constraint.complement_rhs, constraint.complement_rhs.fct, constraint.complement_rhs.set)
    end
    return true
end

"""
    function is_constraint_violated(
        com::CoM,
        constraint::BoolConstraint,
        fct,
        set::XNorSet,
    )

Check whether one side is solved and other is violated
"""
function is_constraint_violated(
    com::CoM,
    constraint::BoolConstraint,
    fct,
    set::XNorSet,
)
    lhs_solved = is_constraint_solved(com, constraint.lhs, constraint.lhs.fct, constraint.lhs.set)
    rhs_solved = is_constraint_solved(com, constraint.rhs, constraint.rhs.fct, constraint.rhs.set)
    # neither of them is solved => it's not violated yet
    if !lhs_solved && !rhs_solved
        return false
    end
    if lhs_solved && rhs_solved
        return false
    end

    lhs_violated = is_lhs_constraint_violated(com, constraint)
    rhs_violated = is_rhs_constraint_violated(com, constraint)
    if lhs_solved && rhs_violated
        return true
    end
    if rhs_solved && lhs_violated
        return true
    end

    return false
end

"""
    still_feasible(com::CoM, constraint::XNorConstraint, fct, set::XNorSet, vidx::Int, value::Int)

Return whether the constraint can be still fulfilled when setting a variable with index `vidx` to `value`.
**Attention:** This assumes that it isn't violated before.
"""
function still_feasible(
    com::CoM,
    constraint::XNorConstraint,
    fct,
    set::XNorSet,
    vidx::Int,
    value::Int,
)
    lhs_violated = is_constraint_violated(com, constraint.lhs, constraint.lhs.fct, constraint.lhs.set)
    lhs_feasible = !lhs_violated
    lhs_solved = false
    if lhs_feasible
        lhs_solved = is_constraint_solved_when_fixed(com, constraint.lhs, constraint.lhs.fct, constraint.lhs.set, vidx, value)
    end
    if lhs_feasible && !lhs_solved
        lhs_indices = constraint.lhs.indices
        for i in 1:length(lhs_indices)
            if lhs_indices[i] == vidx
                lhs_feasible = still_feasible(com, constraint.lhs, constraint.lhs.fct, constraint.lhs.set, vidx, value)
                break
            end
        end
    end
    rhs_violated = is_constraint_violated(com, constraint.rhs, constraint.rhs.fct, constraint.rhs.set)
    rhs_feasible = !rhs_violated
    rhs_solved = false
    if rhs_feasible
        rhs_solved = is_constraint_solved_when_fixed(com, constraint.rhs, constraint.rhs.fct, constraint.rhs.set, vidx, value)
    end
    # both must be violated or both must be solved in the end
    if !lhs_feasible && rhs_solved
        return false
    end
    if rhs_feasible && !rhs_solved
        rhs_indices = constraint.rhs.indices
        for i in 1:length(rhs_indices)
            if rhs_indices[i] == vidx
                rhs_feasible = still_feasible(com, constraint.rhs, constraint.rhs.fct, constraint.rhs.set, vidx, value)
                break
            end
        end
    end

    if !rhs_feasible && lhs_solved
        return false
    end
    return true
end

"""
    prune_constraint!(com::CS.CoM, constraint::XNorConstraint, fct, set::XNorSet; logs = true)

Reduce the number of possibilities given the `XNorConstraint` by pruning both parts
Return whether still feasible
"""
function prune_constraint!(
    com::CS.CoM,
    constraint::XNorConstraint,
    fct,
    set::XNorSet;
    logs = true,
)
    lhs_violated = is_constraint_violated(com, constraint.lhs, constraint.lhs.fct, constraint.lhs.set)
    rhs_violated = is_constraint_violated(com, constraint.rhs, constraint.rhs.fct, constraint.rhs.set)
    if lhs_violated && rhs_violated
        return true
    end
    # check if one is already solved
    lhs_solved = is_constraint_solved(com, constraint.lhs, constraint.lhs.fct, constraint.lhs.set)
    rhs_solved = is_constraint_solved(com, constraint.rhs, constraint.rhs.fct, constraint.rhs.set)
    if lhs_solved && rhs_solved
        return true
    end

    # if one is solved => prune the other
    if lhs_solved
        activate_rhs!(com, constraint)
        return prune_constraint!(com, constraint.rhs, constraint.rhs.fct, constraint.rhs.set; logs=logs)
    end
    if rhs_solved
        activate_lhs!(com, constraint)
        return prune_constraint!(com, constraint.lhs, constraint.lhs.fct, constraint.lhs.set; logs=logs)
    end

    # if one is violated complement prune the other
    # Todo implement for activated complement constraints
    if lhs_violated && constraint.complement_rhs !== nothing && !implements_activate(typeof(constraint.complement_rhs), typeof(constraint.complement_rhs.fct), typeof(constraint.complement_rhs.set))
        return prune_constraint!(com, constraint.complement_rhs, constraint.complement_rhs.fct, constraint.complement_rhs.set; logs=logs)
    end
    if rhs_violated && constraint.complement_lhs !== nothing && !implements_activate(typeof(constraint.complement_lhs), typeof(constraint.complement_lhs.fct), typeof(constraint.complement_lhs.set))
        return prune_constraint!(com, constraint.complement_lhs, constraint.complement_lhs.fct, constraint.complement_lhs.set; logs=logs)
    end

    return true
end


text:
#Code written by Scott Hanna 6/27/2019
#Last Update: 8/14/2019
#Contact: shanna7@jhu.edu
#
#The Modulus type is here defined as the residue class of a modulo n where
#[a] = a + nZ := {a +nz | z ∈ Z}
#the quotient set with respect to Rn, is often denoted by Z/nZ.
#Thus, Z/nZ := {[0], [1], . . . , [n − 1]}.
#
#The set of all classes modulo n is a field only in the case that n is prime.
#i.e. it only guarantees an inverse element if n is prime.  If n is not prime
#it is possible that a particular value, a,has an inverse (that is a*a^-1 = 1),
#but the set is still not considered a field.
module ModulusN

using LinearAlgebra, Random
import Base: convert,promote_rule,+,-,*,/,zero,one,abs,abs2,^,<,<=,>,>=,==,isequal,inv,isnan
import LinearAlgebra: conj
import Random: rand

export getVal,getMod,Modulus,hasInverse

struct Modulus{N} <: Number
    val::Real
    function Modulus{N}(val) where N
        @assert N isa Real
        @assert N != 0
        new{N}(mod(val,N))
    end
    Modulus{N}(x::Modulus{N}) where N = Modulus{N}(x.val)
end

convert(::Modulus{N}, x::M) where {N,M<:Real} = Modulus{N}(x)
promote_rule(::Type{Modulus{N}}, ::Type{M}) where {N,M<:Real} = Modulus{N}
promote_rule(::Type{M}, ::Type{Modulus{N}}) where {N,M<:Real} = Modulus{N}

getVal(x::Modulus{N}) where N = return x.val
getMod(x::Modulus{N}) where N = return N

+(x::Modulus{N},y::Modulus{N}) where N = Modulus{N}(x.val + y.val)
-(x::Modulus{N},y::Modulus{N}) where N = Modulus{N}(x.val - y.val)
*(x::Modulus{N},y::Modulus{N}) where N = Modulus{N}(x.val * y.val)
/(x::Modulus{N},y::Modulus{N}) where N = inv(y) * x
-(x::Modulus{N}) where N = Modulus{N}(-x.val)
zero(x::Modulus{N}) where N = Modulus{N}(0)
one(x::Modulus{N}) where N = Modulus{N}(1)
abs(x::Modulus{N}) where N = Modulus{N}(abs(x.val))
^(x::Modulus{N}, y::Int) where N = Modulus{N}(x.val^y)
abs2(x::Modulus{N}) where N = abs(x)^2

<(x::Modulus{N}, y::Modulus{N}) where N = x.val < y.val
<=(x::Modulus{N}, y::Modulus{N}) where N = x.val <= y.val
>(x::Modulus{N}, y::Modulus{N}) where N = x.val > y.val
>=(x::Modulus{N}, y::Modulus{N}) where N = x.val >= y.val
==(x::Modulus{N}, y::Modulus{M}) where {N,M} = x.val == y.val && N==M
isequal(x::Modulus{N}, y::Modulus{M}) where {N,M} = x.val == y.val && N==M

isnan(x::Modulus{N}) where N = isnan(x.val)
conj(x::Modulus{N}) where N = abs(x)

function inv(x::Modulus{N}) where N
    divs,v = gcdx(x.val,N)
    if divs != 1 || N == 1
        error("No Inverse Exists")
    end
    return Modulus{N}(v)
end

function hasInverse(x::Modulus{N}) where N
    return gcd(x.val,N)==1 && N != 1
end

rand(rng::AbstractRNG, ::Random.SamplerType{Modulus{N}}) where N = Modulus{N}(rand(rng,Int))

end #end module ModulusN

max_stars_repo_name: sphanna/ModulusN

max_issues_repo_name: sphanna/ModulusN

max_forks_repo_name: sphanna/ModulusN


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

# Thanks @tlienart

using Franklin, JSON

using Markdown, Dates

include("youtube_videos.jl")

const DATEFMT = dateformat"yyyy-mm-dd HH:MMp"
const TZ = "America/New_York"

function hfun_doc(params)
    fname = join(params[1:max(1, length(params)-2)], " ")
    head = params[end-1]
    type = params[end]
    doc = eval(Meta.parse("@doc $fname"))
    txt = Markdown.plain(doc)
    # possibly further processing here
    body = Franklin.fd2html(txt, internal=true)
    return """

            $head

$type

$body

    """
end

function hfun_youtube_placeholder(params)
    id = params[1]
    return """

    """
end

function hfun_yt_tsp(params)
    start = params[1]
    final = params[2]
    s = Time(0) + Second(start)
    if second(s) < 10
      display = "$(minute(s)):0$(second(s))"
    elseif hour(s) < 1
      display = "$(minute(s)):$(second(s))"
    else
      display = display = "$(hour(s)):0$(minute(s)):$(second(s))"
    end
    title = join(params[3:end], " ")
    return """

     $display $title
    """
end

function hfun_youtube(params)
    id = params[1]
    return """

    """
end

function hfun_youtube_start(params)
  id = params[1]
  start = params[2]
  return """

  """
end

function hfun_showtime(params)
    id = params[1]
    str = locvar(id)
    if isnothing(str)
        @warn "Unknown datetime variable $str"
        return ""
    end
    try
        DateTime(str, DATEFMT)
    catch err
        @warn "There was an error parsing date $str, the format is yyyy-mm-dd HH:MMp (see ?DateFormat)"
        rethrow(err)
    end
end

function parse_duration(str)
    str = replace(str, r"^PT"=>"")
    hrex, mrex, srex = Regex.(string.("^([0-9]+)", ["H","M","S"]))

    t = 0
    hmatch = match(hrex, str)
    if !isnothing(hmatch)
        h = parse(Int, hmatch[1])
        t += 60*60*h
        str = replace(str, hrex=>"")
    end

    mmatch = match(mrex, str)
    if !isnothing(mmatch)
        m = parse(Int, mmatch[1])
        t += 60*m
        str = replace(str, mrex=>"")
    end

    smatch = match(srex, str)
    if !isnothing(smatch)
        s = parse(Int, smatch[1])
        t += s
        str = replace(str, srex=>"")
    end

    t
end

function hfun_go_live()
    seq = locvar("sequence")
    airtime = locvar("airtime")

    if isnothing(seq)
        @warn "airtime set, but no `sequence` variable not defined." *
        "sequence is an array of video IDs to play in order on this page"
    end

    vid_ids = [get(videos, s, s) for s in seq]

    f = tempname()
    # Get the duration of each video
    download("https://www.googleapis.com/youtube/v3/videos?id=$(join(vid_ids, ","))&part=contentDetails&key=AIzaSyDZhbWHc2PTEFTx173MaTgddnWCGPqdbB8", f)
    dict = JSON.parse(String(read(f)))

    durations = [parse_duration(video["contentDetails"]["duration"])
                 for video in dict["items"]]

    jrepr(x) = sprint(io->JSON.print(io, x))
    """

    """
end

max_stars_repo_name: thanktua-spp/website

max_issues_repo_name: thanktua-spp/website

max_forks_repo_name: thanktua-spp/website


text:
using Shrike
using Shrike: traverse_tree, getdepth, leafsize, _allknn, get_dists
using Test

@testset "Constructor" begin
    X = rand(250, 16)
    shi = ShrikeIndex(X, 3, 3)
    # Check for self consistency
    @test all((size(shi.splits) .== (shi.ntrees, 2^shi.depth - 1)))
    @test all(size(shi.indexes) .== (2^shi.depth, shi.ntrees))
    @test all(size(shi.random_vectors) .== (shi.ndims, shi.ntrees * shi.depth))

    shi = ShrikeIndex(X, depth=5, ntrees=3)
    # Check for self consistency
    @test all((size(shi.splits) .== (shi.ntrees, 2^shi.depth - 1)))
    @test all(size(shi.indexes) .== (2^shi.depth, shi.ntrees))
    @test all(size(shi.random_vectors) .== (shi.ndims, shi.ntrees * shi.depth))

    shi = ShrikeIndex(X, 3)
    # Check for self consistency
    @test all((size(shi.splits) .== (shi.ntrees, 2^shi.depth - 1)))
    @test all(size(shi.indexes) .== (2^shi.depth, shi.ntrees))
    @test all(size(shi.random_vectors) .== (shi.ndims, shi.ntrees * shi.depth))

    shi = ShrikeIndex(X, 3, depth=3, ntrees=10)
    # Check for self consistency
    @test all((size(shi.splits) .== (shi.ntrees, 2^shi.depth - 1)))
    @test all(size(shi.indexes) .== (2^shi.depth, shi.ntrees))
    @test all(size(shi.random_vectors) .== (shi.ndims, shi.ntrees * shi.depth))

end

# Check that for each data point used to build the tree
# traverse_tree returns leaf nodes containing the given datapoint
@testset "Traverse" begin
    X = rand(250, 16)
    shi = ShrikeIndex(X, 3, 3)
    for j =1:size(X)[2]
        @test all(map(idxs -> j in idxs, traverse_tree(shi, reshape(X[:, j], :, 1))))
    end
end

@testset "Nearest Neighbors" begin
    k = 10
    npoints = 1000
    X = randn(300, npoints)
    shi = ShrikeIndex(X)
    @test length(ann(shi, X[:, 1:1], k; vote_cutoff=1)) == k
    @test all(size(allknn(shi, k, ne_iters=1)) .== (npoints, k))
    @test (knngraph(shi, k, ne_iters=1)).ne == k * npoints
    shi100 = ShrikeIndex(X, ntrees=100)
    @test length(ann(shi, X[:, 1:1], k; vote_cutoff=2)) == k
    @test all(size(allknn(shi, k, ne_iters=2)) .== (npoints, k))
    @test (knngraph(shi, k, ne_iters=2)).ne == k * npoints
    @test (knngraph(shi, k, ne_iters=2, vote_cutoff=2)).ne == k * npoints
end

@testset "Utils" begin
    l = 10
    k = 5
    npoints = 1010
    X = randn(10, npoints)
    @test leafsize(ShrikeIndex(X, depth=getdepth(l, 100))) >= l
    ne = _allknn(ShrikeIndex(X, ntrees=2), k)
    @test length(get_dists(ne[1])) == k
end


text:
module MLKernelsMissing

export
    kPCA, fit, transform,
    jacobian,
    Kernel,
        MercerKernel,
            AbstractExponentialKernel,
                ExponentialKernel,
                LaplacianKernel,
                SquaredExponentialKernel,
                GaussianKernel,
                RadialBasisKernel,
                GammaExponentialKernel,
            AbstractRationalQuadraticKernel,
                RationalQuadraticKernel,
                GammaRationalQuadraticKernel,
            MaternKernel,
            LinearKernel,
            PolynomialKernel,
            ExponentiatedKernel,
            PeriodicKernel,
        NegativeDefiniteKernel,
            PowerKernel,
            LogKernel,
        SigmoidKernel

import MLKernels
import MLKernels:
    BaseFunction,
    Orientation,
    kappa,
    kernelmatrix, kernelmatrix!,
    basematrix!,
    basefunction,
    checkdimensions, subvector,
    base_initiate,
    base_aggregate,
    base_return,
    unsafe_base_evaluate,

# reexported
    Kernel,

# reexported, so we don't have to say `using MLKernels`
        MercerKernel,
            AbstractExponentialKernel,
                ExponentialKernel,
                LaplacianKernel,
                SquaredExponentialKernel,
                GaussianKernel,
                RadialBasisKernel,
                GammaExponentialKernel,
            AbstractRationalQuadraticKernel,
                RationalQuadraticKernel,
                GammaRationalQuadraticKernel,
            MaternKernel,
            LinearKernel,
            PolynomialKernel,
            ExponentiatedKernel,
            PeriodicKernel,
        NegativeDefiniteKernel,
            PowerKernel,
            LogKernel,
        SigmoidKernel

import LinearAlgebra: Hermitian, eigen, copytri!
import Statistics: mean

# reexported
import StatsBase: fit
import MultivariateStats: transform

include("util.jl")
include("basematrix.jl")
include("kernelmatrix.jl")
include("kpca.jl")
include("jacobian.jl")
# include("datasets.jl")

end # module MLKernels


text:
x = 0.0
for i in 1:1000000
    global x += 0.1
end

println("result=", x)


text:
#### EDMF model
using DocStringExtensions
using CLIMAParameters: AbstractEarthParameterSet
using CLIMAParameters.Atmos.EDMF
using CLIMAParameters.SubgridScale

"""
    EntrainmentDetrainment

An Entrainment-Detrainment model for EDMF, containing
all related model and free parameters.

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct EntrainmentDetrainment{FT <: AbstractFloat}
    "Entrainment TKE scale"
    c_λ::FT
    "Entrainment factor"
    c_ε::FT
    "Detrainment factor"
    c_δ::FT
    "Turbulent Entrainment factor"
    c_t::FT
    "Detrainment RH power"
    β::FT
    "Logistic function scale ‵[1/s]‵"
    μ_0::FT
    "Updraft mixing fraction"
    χ::FT
    "Minimum updraft velocity"
    w_min::FT
    "Exponential area limiter scale"
    lim_ϵ::FT
    "Exponential area limiter amplitude"
    lim_amp::FT
end

"""
    EntrainmentDetrainment{FT}(param_set) where {FT}

Constructor for `EntrainmentDetrainment` for EDMF, given:
 - `param_set`, an AbstractEarthParameterSet
"""
function EntrainmentDetrainment{FT}(
    param_set::AbstractEarthParameterSet,
) where {FT}
    c_λ_ = c_λ(param_set)
    c_ε_ = c_ε(param_set)
    c_δ_ = c_δ(param_set)
    c_t_ = c_t(param_set)
    β_ = β(param_set)
    μ_0_ = μ_0(param_set)
    χ_ = χ(param_set)
    w_min_ = w_min(param_set)
    lim_ϵ_ = lim_ϵ(param_set)
    lim_amp_ = lim_amp(param_set)

    args = (c_λ_, c_ε_, c_δ_, c_t_, β_, μ_0_, χ_, w_min_, lim_ϵ_, lim_amp_)

    return EntrainmentDetrainment{FT}(args...)
end

"""
    SubdomainModel

A subdomain model for EDMF, containing
all related model and free parameters.

TODO: `a_max` is valid for all subdomains,
    but it is insufficient to ensure `a_en`
    is not negative. Limits can be imposed
    for updrafts, but this is a limit
    is dictated by 1 - Σᵢ aᵢ, which must be
    somehow satisfied by regularizing prognostic
    source terms.

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct SubdomainModel{FT <: AbstractFloat}
    "Minimum area fraction for any subdomain"
    a_min::FT
    "Maximum area fraction for any subdomain"
    a_max::FT
end

function SubdomainModel(
    ::Type{FT},
    N_up;
    a_min::FT = FT(0),
    a_max::FT = 1 - N_up * a_min,
) where {FT}
    return SubdomainModel(; a_min = a_min, a_max = a_max)
end

"""
    SurfaceModel

A surface model for EDMF, containing all boundary
values and parameters needed by the model.

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct SurfaceModel{FT <: AbstractFloat, SV}
    "Area"
    a::FT
    "Surface covariance stability coefficient"
    ψϕ_stab::FT
    "Square ratio of rms turbulent velocity to friction velocity"
    κ_star²::FT
    "Updraft normalized standard deviation at the surface"
    upd_surface_std::SV
    # The following will be deleted after SurfaceFlux coupling
    "Liquid water potential temperature ‵[k]‵"
    θ_liq::FT = 299.1
    "Specific humidity ‵[kg/kg]‵"
    q_tot::FT = 22.45e-3
    "Sensible heat flux ‵[w/m^2]‵"
    shf::FT = 9.5
    "Latent heat flux ‵[w/m^2]‵"
    lhf::FT = 147.2
    "Friction velocity"
    ustar::FT = 0.28
    "Monin - Obukhov length"
    obukhov_length::FT = 0
    "Height of the lowest level"
    zLL::FT = 60
end

"""
    SurfaceModel{FT}(N_up, param_set) where {FT}

Constructor for `SurfaceModel` for EDMF, given:
 - `N_up`, the number of updrafts
 - `param_set`, an AbstractEarthParameterSet
"""
function SurfaceModel{FT}(N_up, param_set::AbstractEarthParameterSet) where {FT}
    a_surf_ = a_surf(param_set)
    κ_star²_ = κ_star²(param_set)
    ψϕ_stab_ = ψϕ_stab(param_set)

    if a_surf_ > FT(0)
        upd_surface_std = SVector(
            ntuple(N_up) do i
                percentile_bounds_mean_norm(
                    1 - a_surf_ + (i - 1) * FT(a_surf_ / N_up),
                    1 - a_surf_ + i * FT(a_surf_ / N_up),
                    1000,
                )
            end,
        )
    else
        upd_surface_std = SVector(ntuple(i -> FT(0), N_up))
    end
    SV = typeof(upd_surface_std)
    return SurfaceModel{FT, SV}(;
        upd_surface_std = upd_surface_std,
        a = a_surf_,
        κ_star² = κ_star²_,
        ψϕ_stab = ψϕ_stab_,
    )
end

"""
    NeutralDrySurfaceModel

A surface model for EDMF simulations in a
dry, neutral environment, containing all boundary
values and parameters needed by the model.

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct NeutralDrySurfaceModel{FT <: AbstractFloat}
    "Area"
    a::FT
    "Square ratio of rms turbulent velocity to friction velocity"
    κ_star²::FT
    "Friction velocity"
    ustar::FT = 0.3
    "Height of the lowest level"
    zLL::FT = 60
    "Monin - Obukhov length"
    obukhov_length::FT = 0
end

"""
    NeutralDrySurfaceModel{FT}(N_up, param_set) where {FT}

Constructor for `NeutralDrySurfaceModel` for EDMF, given:
 - `param_set`, an AbstractEarthParameterSet
"""
function NeutralDrySurfaceModel{FT}(
    param_set::AbstractEarthParameterSet,
) where {FT}
    a_surf_ = a_surf(param_set)
    κ_star²_ = κ_star²(param_set)
    return NeutralDrySurfaceModel{FT}(; a = a_surf_, κ_star² = κ_star²_)
end

"""
    PressureModel

A pressure model for EDMF, containing
all related model and free parameters.

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct PressureModel{FT <: AbstractFloat}
    "Pressure drag"
    α_d::FT
    "Pressure advection"
    α_a::FT
    "Pressure buoyancy"
    α_b::FT
    "Minimum diagnostic updraft height for closures"
    H_up_min::FT
end

"""
    PressureModel{FT}(param_set) where {FT}

Constructor for `PressureModel` for EDMF, given:
 - `param_set`, an AbstractEarthParameterSet
"""
function PressureModel{FT}(param_set::AbstractEarthParameterSet) where {FT}
    α_d_ = α_d(param_set)
    α_a_ = α_a(param_set)
    α_b_ = α_b(param_set)
    H_up_min_ = H_up_min(param_set)

    args = (α_d_, α_a_, α_b_, H_up_min_)

    return PressureModel{FT}(args...)
end

"""
    MixingLengthModel

A mixing length model for EDMF, containing
all related model and free parameters.

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct MixingLengthModel{FT <: AbstractFloat}
    "dissipation coefficient"
    c_d::FT
    "Eddy Viscosity"
    c_m::FT
    "Static Stability coefficient"
    c_b::FT
    "Empirical stability function coefficient"
    a1::FT
    "Empirical stability function coefficient"
    a2::FT
    "Von Karmen constant"
    κ::FT
    "Prandtl number empirical coefficient"
    ω_pr::FT
    "Prandtl number scale"
    Pr_n::FT
    "Critical Richardson number"
    Ri_c::FT
    "smooth minimum's fractional upper bound"
    smin_ub::FT
    "smooth minimum's regularization minimum"
    smin_rm::FT
    "Maximum mixing length"
    max_length::FT
    "Random small number variable that should be addressed"
    random_minval::FT
end

"""
    MixingLengthModel{FT}(param_set) where {FT}

Constructor for `MixingLengthModel` for EDMF, given:
 - `param_set`, an AbstractEarthParameterSet
"""
function MixingLengthModel{FT}(param_set::AbstractEarthParameterSet) where {FT}
    c_d_ = c_d(param_set)
    c_m_ = c_m(param_set)
    c_b_ = c_b(param_set)
    a1_ = a1(param_set)
    a2_ = a2(param_set)
    κ = von_karman_const(param_set)
    ω_pr_ = ω_pr(param_set)
    Pr_n_ = Pr_n(param_set)
    Ri_c_ = Ri_c(param_set)
    smin_ub_ = smin_ub(param_set)
    smin_rm_ = smin_rm(param_set)
    max_length = 1e6
    random_minval = 1e-9

    args = (
        c_d_,
        c_m_,
        c_b_,
        a1_,
        a2_,
        κ,
        ω_pr_,
        Pr_n_,
        Ri_c_,
        smin_ub_,
        smin_rm_,
        max_length,
        random_minval,
    )

    return MixingLengthModel{FT}(args...)
end

abstract type AbstractStatisticalModel end
struct SubdomainMean <: AbstractStatisticalModel end
struct GaussQuad <: AbstractStatisticalModel end
struct LogNormalQuad <: AbstractStatisticalModel end

"""
    MicrophysicsModel

A microphysics model for EDMF, containing
all related model and free parameters and
assumed subdomain distributions.

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct MicrophysicsModel{FT <: AbstractFloat, SM}
    "Subdomain statistical model"
    statistical_model::SM
end

"""
    MicrophysicsModel(
        FT;
        statistical_model = SubdomainMean()
    )

Constructor for `MicrophysicsModel` for EDMF, given:
 - `FT`, the float type used
 - `statistical_model`, the assumed environmental distribution
            of thermodynamic variables.
"""
function MicrophysicsModel(FT; statistical_model = SubdomainMean())
    args = (statistical_model,)
    return MicrophysicsModel{FT, typeof(statistical_model)}(args...)
end

"""
    Environment <: BalanceLaw
A `BalanceLaw` for the environment subdomain arising in EDMF.
"""
Base.@kwdef struct Environment{FT <: AbstractFloat, N_quad} <: BalanceLaw end

"""
    Updraft <: BalanceLaw
A `BalanceLaw` for the updraft subdomains arising in EDMF.
"""
Base.@kwdef struct Updraft{FT <: AbstractFloat} <: BalanceLaw end

abstract type Coupling end

"""
    Decoupled <: Coupling

Dispatch on decoupled model (default)

 - The EDMF SGS tendencies do not modify the grid-mean equations.
"""
struct Decoupled <: Coupling end

"""
    Coupled <: Coupling

Dispatch on coupled model

 - The EDMF SGS tendencies modify the grid-mean equations.
"""
struct Coupled <: Coupling end

"""
    EDMF <: TurbulenceConvectionModel

A turbulence convection model for the EDMF
scheme, containing all closure models and
free parameters.

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct EDMF{
    FT <: AbstractFloat,
    N,
    UP,
    EN,
    ED,
    P,
    S,
    MP,
    ML,
    SD,
    C,
} <: TurbulenceConvectionModel
    "Updrafts"
    updraft::UP
    "Environment"
    environment::EN
    "Entrainment-Detrainment model"
    entr_detr::ED
    "Pressure model"
    pressure::P
    "Surface model"
    surface::S
    "Microphysics model"
    micro_phys::MP
    "Mixing length model"
    mix_len::ML
    "Subdomain model"
    subdomains::SD
    "Coupling mode"
    coupling::C
end

"""
    EDMF(
        FT, N_up, N_quad, param_set;
        updraft = ntuple(i -> Updraft{FT}(), N_up),
        environment = Environment{FT, N_quad}(),
        entr_detr = EntrainmentDetrainment{FT}(param_set),
        pressure = PressureModel{FT}(param_set),
        surface = SurfaceModel{FT}(N_up, param_set),
        micro_phys = MicrophysicsModel(FT),
        mix_len = MixingLengthModel{FT}(param_set),
        subdomain = SubdomainModel(FT, N_up),
    )
Constructor for `EDMF` subgrid-scale scheme, given:
 - `FT`, the AbstractFloat type used
 - `N_up`, the number of updrafts
 - `N_quad`, the quadrature order. `N_quad^2` is
        the total number of quadrature points
        used for environmental distributions.
 - `updraft`, a tuple containing N_up updraft BalanceLaws
 - `environment`, the environment BalanceLaw
 - `entr_detr`, an `EntrainmentDetrainment` model
 - `pressure`, a `PressureModel`
 - `surface`, a `SurfaceModel`
 - `micro_phys`, a `MicrophysicsModel`
 - `mix_len`, a `MixingLengthModel`
 - `subdomain`, a `SubdomainModel`
 - `coupling`, a coupling type
"""
function EDMF(
    FT,
    N_up,
    N_quad,
    param_set;
    updraft = ntuple(i -> Updraft{FT}(), N_up),
    environment = Environment{FT, N_quad}(),
    entr_detr = EntrainmentDetrainment{FT}(param_set),
    pressure = PressureModel{FT}(param_set),
    surface = SurfaceModel{FT}(N_up, param_set),
    micro_phys = MicrophysicsModel(FT),
    mix_len = MixingLengthModel{FT}(param_set),
    subdomain = SubdomainModel(FT, N_up),
    coupling = Decoupled(),
)
    args = (
        updraft,
        environment,
        entr_detr,
        pressure,
        surface,
        micro_phys,
        mix_len,
        subdomain,
        coupling,
    )
    return EDMF{FT, N_up, typeof.(args)...}(args...)
end

import ClimateMachine.TurbulenceConvection: turbconv_sources, turbconv_bcs

"""
    EDMFTopBC <: TurbConvBC

Boundary conditions for the top of the EDMF.
"""
struct EDMFTopBC <: TurbConvBC end

"""
    EDMFBottomBC <: TurbConvBC

Boundary conditions for the bottom of the EDMF.
"""
struct EDMFBottomBC <: TurbConvBC end

n_updrafts(m::EDMF{FT, N_up}) where {FT, N_up} = N_up
n_updrafts(m::TurbulenceConvectionModel) = 0
turbconv_filters(m::TurbulenceConvectionModel) = ()
turbconv_filters(m::EDMF) = (
    "turbconv.environment.ρatke",
    "turbconv.environment.ρaθ_liq_cv",
    "turbconv.environment.ρaq_tot_cv",
    "turbconv.updraft",
)
n_quad_points(m::Environment{FT, N_quad}) where {FT, N_quad} = N_quad
turbconv_sources(m::EDMF) = ()
turbconv_bcs(::EDMF) = (EDMFBottomBC(), EDMFTopBC())


text:
using DataFrames
using Serialization

include("../utils.jl")
include("simulations.jl")
include("../../parallel/dataframes.jl")

function generateP1dataset(config::Dict)
    dfs = [DataFrame(df = Vector{DataFrame}(), Δm = Vector{Float64}(), ΔJ = Vector{Float64}(), friction = Vector{Float64}()) for _ in 1:4]
    nsteps = 2*Int(1/config["Δtsim"])  # Equivalent to 2 seconds
    threadlock = ReentrantLock()
    for i in 1:config["Ndfs"]
        @info ("Working on $(i)/$(config["Ndfs"]) P1 datasets")
        friction = rand()
        ΔJ = 1 + config["Σ"]["ΔJ"]*2(rand()-0.5)
        Δm = 1 + config["Σ"]["Δm"]*2(rand()-0.5)
        exp1 = () -> simplependulum2D(nsteps, Δt=config["Δtsim"], θstart=(rand() - 0.5) * π, ωstart = 2(rand()-0.5), Δm = Δm, ΔJ = ΔJ, friction=friction, threadlock = threadlock)[1]
        exp2 = () -> simplependulum2D(nsteps, Δt=config["Δtsim"], θstart=((rand()/2 + 0.5)*rand((-1,1))) * π, ωstart = 2(rand()-0.5), Δm = Δm, ΔJ = ΔJ, friction=friction, threadlock = threadlock)[1]  # [-π:-π/2; π/2:π] [-π:π]
        exptest = () -> simplependulum2D(nsteps, Δt=config["Δtsim"], θstart=(rand() - 0.5)*2π, ωstart = 2(rand()-0.5), Δm = Δm, ΔJ = ΔJ, friction=friction, threadlock = threadlock)[1]
        _dfs = generate_dataframes(config, exp1, exp2, exptest; generate_uniform = true)
        for (df, _df) in zip(dfs, _dfs)
            push!(df, (_df, Δm, ΔJ, friction))
        end
    end
    return dfs
end

function generateP2dataset(config::Dict)
    dfs = [DataFrame(df = Vector{DataFrame}(), Δm = Vector{Vector{Float64}}(), ΔJ = Vector{Vector{Float64}}(), friction = Vector{Vector{Float64}}()) for _ in 1:4]
    nsteps = 2*Int(1/config["Δtsim"])  # Equivalent to 2 seconds
    threadlock = ReentrantLock()
    for i in 1:config["Ndfs"]
        @info ("Working on $(i)/$(config["Ndfs"]) P2 datasets")
        friction = rand(2) .* [2, 0.5]
        ΔJ = ones(2) + config["Σ"]["ΔJ"].*2(rand(2).-0.5)
        Δm = ones(2) + config["Σ"]["Δm"].*2(rand(2).-0.5)
        exp1 = () -> doublependulum2D(nsteps, Δt=config["Δtsim"], θstart=(rand(2).-0.5) .* [π, 2π], Δm = Δm, ΔJ = ΔJ, friction=friction, threadlock = threadlock)[1]
        exp2 = () -> doublependulum2D(nsteps, Δt=config["Δtsim"], θstart=[(rand()/2 + 0.5)*rand([-1,1]), 2(rand()-0.5)] .* π, Δm = Δm, ΔJ = ΔJ, friction=friction, threadlock = threadlock)[1]    # [-π:-π/2; π/2:π] [-π:π]
        exptest = () -> doublependulum2D(nsteps, Δt=config["Δtsim"], θstart=(rand(2).-0.5).*2π, Δm = Δm, ΔJ = ΔJ, friction=friction, threadlock = threadlock)[1]
        _dfs = generate_dataframes(config, exp1, exp2, exptest; generate_uniform = true)
        for (df, _df) in zip(dfs, _dfs)
            push!(df, (_df, Δm, ΔJ, friction))
        end
    end
    return dfs
end

function generateCPdataset(config::Dict)
    dfs = [DataFrame(df = Vector{DataFrame}(), Δm = Vector{Vector{Float64}}(), ΔJ = Vector{Float64}(), friction = Vector{Vector{Float64}}()) for _ in 1:4]
    nsteps = 2*Int(1/config["Δtsim"])  # Equivalent to 2 seconds
    threadlock = ReentrantLock()
    for i in 1:config["Ndfs"]
        @info ("Working on $(i)/$(config["Ndfs"]) CP datasets")
        ΔJ = 1 + config["Σ"]["ΔJ"].*2(rand()-0.5)  # Only change inertia of the pendulum
        Δm = ones(2) + config["Σ"]["Δm"]*2(rand(2).-0.5)
        friction = rand(2) .* [1., 0.5]
        nsteps = 2*Int(1/config["Δtsim"])  # Equivalent to 2 seconds
        exp1 = () -> cartpole(nsteps, Δt=config["Δtsim"], xstart=rand()-0.5, θstart=(rand()-0.5)π, vstart=2(rand()-0.5), ωstart=2(rand()-0.5), Δm = Δm, ΔJ = ΔJ, friction=friction, threadlock = threadlock)[1]
        exp2 = () -> cartpole(nsteps, Δt=config["Δtsim"], xstart=rand()-0.5, θstart=(rand()/2+0.5)*rand([-1,1])π, vstart=2(rand()-0.5), ωstart=2(rand()-0.5), Δm = Δm, ΔJ = ΔJ, friction=friction, threadlock = threadlock)[1]
        exptest = () -> cartpole(nsteps, Δt=config["Δtsim"], xstart=rand()-0.5, θstart=2π*(rand()-0.5), vstart=2(rand()-0.5), ωstart=2(rand()-0.5), Δm = Δm, ΔJ = ΔJ, friction=friction, threadlock = threadlock)[1]
        _dfs = generate_dataframes(config, exp1, exp2, exptest; generate_uniform = true)
        for (df, _df) in zip(dfs, _dfs)
            push!(df, (_df, Δm, ΔJ, friction))
        end
    end
    return dfs
end

function generateFBdataset(config::Dict)
    dfs = [DataFrame(df = Vector{DataFrame}(), Δm = Vector{Vector{Float64}}(), ΔJ = Vector{Vector{Float64}}(), friction = Vector{Vector{Float64}}()) for _ in 1:4]
    nsteps = 2*Int(1/config["Δtsim"])  # Equivalent to 2 seconds
    threadlock = ReentrantLock()
    for i in 1:config["Ndfs"]
        @info ("Working on $(i)/$(config["Ndfs"]) FB datasets")
        ΔJ = ones(4) + config["Σ"]["ΔJ"]*2(rand(4).-0.5)
        Δm = ones(4) + config["Σ"]["Δm"]*2(rand(4).-0.5)
        friction = rand(4) .* [2., 0.5, 2., 0.5]
        exp1 = () -> fourbar(nsteps, Δt=config["Δtsim"], θstart=(rand(2).-0.5)π, Δm = Δm, ΔJ = ΔJ, friction=friction, threadlock = threadlock)[1]
        exp2 = () -> fourbar(nsteps, Δt=config["Δtsim"], θstart=(rand(2).-0.5)π, Δm = Δm, ΔJ = ΔJ, friction=friction, threadlock = threadlock)[1]
        exptest = () -> fourbar(nsteps, Δt=config["Δtsim"], θstart=(rand(2).-0.5)π, Δm = Δm, ΔJ = ΔJ, friction=friction, threadlock = threadlock)[1]
        _dfs = generate_dataframes(config, exp1, exp2, exptest; generate_uniform = true)
        for (df, _df) in zip(dfs, _dfs)
            push!(df, (_df, Δm, ΔJ, friction))
        end
    end
    return dfs
end

function savedatasets(id, dfs, dfnames)
    root = joinpath(dirname(dirname(dirname(@__FILE__))), "datasets")
    if !Base.Filesystem.isdir(root)
        Base.Filesystem.mkpath(root)
    end
    for (df, dfname) in zip(dfs, dfnames)
        path = joinpath(root, id*"_"*dfname*".jls")
        serialize(path, df)
    end
end

function loaddatasets(id)
    root = joinpath(dirname(dirname(dirname(@__FILE__))), "datasets")
    path = joinpath(root, id*"_trainset.jls")
    traindf = deserialize(path)
    path = joinpath(root, id*"_testset.jls")
    testdf = deserialize(path)
    return traindf, testdf
end

function getconfig()
    Σ = Dict("Δm" => 0.1, "ΔJ" => 0.1)
    config = Dict("Δtsim" => 0.0001, "Ndfs" => 100, "trainsamples" => 2048, "testsamples" => 100, "simsteps" => 20, "Σ" => Σ)
    return config
end

function generatedataset(id::String, config::Dict)
    "P1" == id && return generateP1dataset(config)
    "P2" == id && return generateP2dataset(config)
    "CP" == id && return generateCPdataset(config)
    "FB" == id && return generateFBdataset(config)
    throw(ArgumentError("Unsupported experiment ID: $id"))
end

function main()
    config = getconfig()
    for id in ["FB"]  # "P1", "P2", "CP", "FB"
        dfs = generatedataset(id, config)
        savedatasets(id, dfs, ["trainset", "testset", "trainset_uniform", "testset_uniform"])
        @info ("Completed dataset $id")
    end
end

# main()


text:
#=
   Logistic regresiion
   Inspired by WebPPL's examples/logisticRegression.wppl

   Summary Statistics (for PG(15))
      parameters       mean       std   naive_se      mcse         ess      rhat   ess_per_sec
          Symbol    Float64   Float64    Float64   Float64     Float64   Float64       Float64

               m     1.0844    0.5870     0.0059    0.0077   6361.4279    1.0009      153.4390
               b     0.1315    0.9292     0.0093    0.0111   6431.4331    0.9999      155.1275
           sigma     0.8779    0.8513     0.0085    0.0107   6452.7250    0.9999      155.6411
            y[1]   -10.7253    6.0611     0.0606    0.0801   6408.9867    1.0008      154.5861
            y[2]    -5.3599    3.2295     0.0323    0.0439   6465.7579    1.0007      155.9555
            y[3]     2.5165    1.7809     0.0178    0.0200   6562.8112    1.0010      158.2964
            y[4]     6.6892    3.7534     0.0375    0.0454   6516.1303    1.0008      157.1705
            y[5]    11.0158    5.9711     0.0597    0.0746   6263.6413    1.0010      151.0804
  labels_post[1]     0.0163    0.1266     0.0013    0.0013   8339.5233    0.9999      201.1511
  labels_post[2]     0.0552    0.2284     0.0023    0.0026   8486.3592    1.0001      204.6928
  labels_post[3]     0.8459    0.3611     0.0036    0.0041   8013.9274    1.0001      193.2977
  labels_post[4]     0.9707    0.1687     0.0017    0.0018   8543.2209    1.0002      206.0643
  labels_post[5]     0.9873    0.1120     0.0011    0.0012   9349.4013    0.9999      225.5096

  The labels_post are the posterior predictive and are fairly close to labels values: [0, 0, 1, 1, 1]

=#

using Turing, StatsPlots, DataFrames

include("jl_utils.jl")

@model function logistic_regression(x, labels)
    m ~ Normal(0, 1)
    b ~ Normal(0, 1)
    sigma ~ Gamma(1, 1)

    n = length(x)
    y = tzeros(n)
    sigmoid = tzeros(n)
    for i in 1:n
        y[i] ~ Normal(m * x[i] + b, sigma)
        sigmoid[i] = 1 / (1 + exp(-1 * y[i]))
        labels[i] ~ Bernoulli(sigmoid[i])
    end

    # Post predictive of labels
    labels_post = tzeros(n)
    for i in 1:n
        labels_post[i] ~ Bernoulli(sigmoid[i])
    end
end

xs = [-10, -5, 2, 6, 10]
labels = [0, 0, 1, 1, 1]
model = logistic_regression(xs,labels)

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(15), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
const rgbTupleType = NTuple{3, Float64}
const rgbaTupleType = NTuple{4, Float64}
const colorTupleType = Union{rgbTupleType, rgbaTupleType}

ARGB(r, g, b, a) = reinterpret(Cairo.ARGB32, UInt32(a)<<24 | UInt32(b)<<16 | UInt32(g)<<8 | UInt32(r))

function argb32ToRGBATuple(n)
    return ((n >> 16) & 255, (n >> 8) & 255, n & 255, (n >> 24) & 255)
end

function normRGBA32Tuple2ARGB32(c::rgbaTupleType)
    bytes = round.(UInt, c .* 255)

    return bytes[1] << 16 | bytes[2] << 8 | bytes[3] | bytes[4] << 24
end

max_stars_repo_name: asdfsdfadfasfasd/Ahorn

max_issues_repo_name: asdfsdfadfasfasd/Ahorn

max_forks_repo_name: asdfsdfadfasfasd/Ahorn


text:
# ---
# title: 315. Count of Smaller Numbers After Self
# id: problem315
# author: zhwang
# date: 2022-03-26
# difficulty: Hard
# categories: Binary Search, Divide and Conquer, Sort, Binary Indexed Tree, Segment Tree
# link: <https://leetcode.com/problems/count-of-smaller-numbers-after-self/description/>
# hidden: true
# ---
#
# You are given an integer array _nums_ and you have to return a new _counts_
# array. The _counts_ array has the property where `counts[i]` is the number of
# smaller elements to the right of `nums[i]`.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [5,2,6,1]
#     Output: [2,1,1,0]
#     Explanation:
#     To the right of 5 there are **2** smaller elements (2 and 1).
#     To the right of 2 there is only **1** smaller element (1).
#     To the right of 6 there is **1** smaller element (1).
#     To the right of 1 there is **0** smaller element.
#
#
#
#
# **Constraints:**
#
#   * `0 <= nums.length <= 10^5`
#   * `-10^4 <= nums[i] <= 10^4`
#
#
## @lc code=start
using LeetCode

## method 1: using merge sort
function count_smaller_method1(arr::Vector{Int})::Vector{Int}
    res = fill(0, length(arr))
    function merge_sort(arr::AbstractVector)::AbstractVector
        ## arr = [(idx1, num1), ...]
        (n = length(arr)) == 1 && return arr
        mid = n >> 1
        lpart = merge_sort(@view(arr[1:mid]))
        rpart = merge_sort(@view(arr[mid + 1:end]))
        return merge_sorted(lpart, rpart)
    end
    function merge_sorted(lpart::AbstractVector, rpart::AbstractVector)::Vector{Tuple{Int,Int}}
        l1, l2 = length(lpart), length(rpart)
        combined = Vector{Tuple{Int,Int}}(undef, l1 + l2)
        p1 = p2 = 1
        for pos in eachindex(combined)
            if p2 > l2 || p1 <= l1 && lpart[p1][2] <= rpart[p2][2]
                combined[pos] = lpart[p1]
                res[lpart[p1][1]] += (p2 - 1)
                p1 += 1
            else
                combined[pos] = rpart[p2]
                p2 += 1
            end
        end
        combined
    end

    merge_sort(collect(enumerate(arr)))
    return res
end

## method 2: using binary search
function search_left_border(left::Int, right::Int, key::Function)
    ## similar to findfirst(key, left:right) but using bs method
    while left <= right
        mid = left + (right - left) >> 1
        if key(mid) ## lies in the right part
            right = mid - 1
        else
            left = mid + 1
        end
    end
    left
end
function count_smaller_method2(arr::Vector{Int})::Vector{Int}
    queue, res = Int[], Array{Int}(undef, length(reverse!(arr)))
    for (i, num) in enumerate(arr)
        loc = search_left_border(1, length(queue), x -> queue[x] >= num)
        res[i] = loc - 1
        insert!(queue, loc, num)
    end
    reverse!(arr)
    return reverse!(res)
end
## @lc code=end

## @lc test=start
@testset "315.count-of-smaller-numbers-after-self.jl" begin
    f1 = count_smaller_method1
    f2 = count_smaller_method2
    @test f1([5,2,6,1]) == f2([5,2,6,1]) == [2, 1, 1, 0]
    @test f1([2,0,1]) == f2([2,0,1]) == [2, 0, 0]
    @test f1([-1]) == f2([-1]) == [0]
    @test f1([-1, -1]) == f2([-1, -1]) == [0, 0]
    @test f1([7,6,5,4,3,2,1]) == f2([7,6,5,4,3,2,1]) == [6, 5, 4, 3, 2, 1, 0]
end
## @lc test=end

## @lc info=start
# link: [solution to 315]()
## @lc info=end