text:
# Hyperbolic network layer from Lensink et al. (2019)
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

export HyperbolicLayer, HyperbolicLayer3D

"""
    HyperbolicLayer(n_in, kernel, stride, pad; action=0, α=1f0, n_hidden=1)
    HyperbolicLayer(n_in, kernel, stride, pad; action=0, α=1f0, n_hidden=1, ndims=2)
    HyperbolicLayer3D(n_in, kernel, stride, pad; action=0, α=1f0, n_hidden=1)

or

    HyperbolicLayer(W, b, stride, pad; action=0, α=1f0)
    HyperbolicLayer3D(W, b, stride, pad; action=0, α=1f0)

Create an invertible hyperbolic coupling layer.

*Input*:

 - `kernel`, `stride`, `pad`: Kernel size, stride and padding of the convolutional operator

 - `action`: String that defines whether layer keeps the number of channels fixed (`0`),
    increases it by a factor of 4 (or 8 in 3D) (`1`) or decreased it by a factor of 4 (or 8) (`-1`).

 - `W`, `b`: Convolutional weight and bias. `W` has dimensions of `(kernel, kernel, n_in, n_in)`.
   `b` has dimensions of `n_in`.

 - `α`: Step size for second time derivative. Default is 1.

 - `n_hidden`: Increase the no. of channels by `n_hidden` in the forward convolution.
    After applying the transpose convolution, the dimensions are back to the input dimensions.

 - `ndims`: Number of dimension of the input (2 for 2D, 3 for 3D)

*Output*:

 - `HL`: Invertible hyperbolic coupling layer

 *Usage:*

 - Forward mode: `X_curr, X_new = HL.forward(X_prev, X_curr)`

 - Inverse mode: `X_prev, X_curr = HL.inverse(X_curr, X_new)`

 - Backward mode: `ΔX_prev, ΔX_curr, X_prev, X_curr = HL.backward(ΔX_curr, ΔX_new, X_curr, X_new)`

 *Trainable parameters:*

 - `HL.W`: Convolutional kernel

 - `HL.b`: Bias

 See also: [`get_params`](@ref), [`clear_grad!`](@ref)
"""
struct HyperbolicLayer <: NeuralNetLayer
    W::Parameter
    b::Parameter
    α::Float32
    stride::Int
    pad::Int
    action
end

@Flux.functor HyperbolicLayer

scale_a = Dict(0 => 1, 1 => 1/2, -1 => 2, "same" => 1)

# Constructor 2D
function HyperbolicLayer(n_in::Int64, kernel::Int64, stride::Int64,
                         pad::Int64; action="same", α=1f0, n_hidden=nothing, ndims=2)

    # Set ouput/hidden dimensions
    n_out = Int(n_in*scale_a[action]^ndims)
    isnothing(n_hidden) && (n_hidden = n_in)

    k = Tuple(kernel for i=1:ndims)
    W = Parameter(glorot_uniform(k..., n_out, n_hidden))
    b = Parameter(zeros(Float32, n_hidden))

    return HyperbolicLayer(W, b, α,stride, pad, action)
end

HyperbolicLayer3D(args...; kw...) =  HyperbolicLayer(args...; kw..., ndims=3)

# Constructor for given weights 2D
function HyperbolicLayer(W::AbstractArray{T, N}, b::AbstractArray{T, 1},
                         stride::Int64, pad::Int64; action=0, α=1f0) where {T, N}

    kernel, n_in, n_hidden = size(W)[N-3:N]

    # Set ouput/hidden dimensions
    n_out = Int(n_in*scale_a[action]^(N-2))
    W = Parameter(W)
    b = Parameter(b)

    return HyperbolicLayer(W, b, α, stride, pad, action)
end

HyperbolicLayer3D(W::AbstractArray{T, N}, b, stride, pad;
                  action=0, α=1f0) where {T, N} =
        HyperbolicLayer(W, b, stride, pad;action=actin, α=α)

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

# Forward pass
function forward(X_prev_in::AbstractArray{T, N}, X_curr_in::AbstractArray{T, N}, HL::HyperbolicLayer) where {T, N}

    # Change dimensions
    if HL.action == 0
        X_prev = identity(X_prev_in)
        X_curr = identity(X_curr_in)
    elseif HL.action == 1
        X_prev = invHaar_unsqueeze(X_prev_in)
        X_curr = invHaar_unsqueeze(X_curr_in)
    elseif HL.action == -1
        X_prev = Haar_squeeze(X_prev_in)
        X_curr = Haar_squeeze(X_curr_in)
    else
        throw("Specified operation not defined.")
    end

    # Symmetric convolution w/ relu activation
    cdims = DCDims(X_curr, HL.W.data; stride=HL.stride, padding=HL.pad)
    if length(size(X_curr)) == 4
        X_conv = conv(X_curr, HL.W.data, cdims) .+ reshape(HL.b.data, 1, 1, :, 1)
    else
        X_conv = conv(X_curr, HL.W.data, cdims) .+ reshape(HL.b.data, 1, 1, 1, :, 1)
    end
    X_relu = ReLU(X_conv)
    X_convT = -∇conv_data(X_relu, HL.W.data, cdims)

    # Update
    X_new = 2f0*X_curr - X_prev + HL.α*X_convT

    return X_curr, X_new
end

# Inverse pass
function inverse(X_curr::AbstractArray{T, N}, X_new::AbstractArray{T, N}, HL::HyperbolicLayer; save=false) where {T, N}
    cdims = DCDims(X_curr, HL.W.data; stride=HL.stride, padding=HL.pad)
    # Symmetric convolution w/ relu activation
    if length(size(X_curr)) == 4
        X_conv = conv(X_curr, HL.W.data, cdims) .+ reshape(HL.b.data, 1, 1, :, 1)
    else
        X_conv = conv(X_curr, HL.W.data, cdims) .+ reshape(HL.b.data, 1, 1, 1, :, 1)
    end
    X_relu = ReLU(X_conv)
    X_convT = -∇conv_data(X_relu, HL.W.data, cdims)

    # Update
    X_prev = 2*X_curr - X_new + HL.α*X_convT

    # Change dimensions
    if HL.action == 0
        X_prev_in = identity(X_prev)
        X_curr_in = identity(X_curr)
    elseif HL.action == -1
        X_prev_in = invHaar_unsqueeze(X_prev)
        X_curr_in = invHaar_unsqueeze(X_curr)
    elseif HL.action == 1
        X_prev_in = Haar_squeeze(X_prev)
        X_curr_in = Haar_squeeze(X_curr)
    else
        throw("Specified operation not defined.")
    end

    if save == false
        return X_prev_in, X_curr_in
    else
        return X_prev_in, X_curr_in, X_conv, X_relu
    end
end

# Backward pass
function backward(ΔX_curr::AbstractArray{T, N}, ΔX_new::AbstractArray{T, N}, X_curr::AbstractArray{T, N}, X_new::AbstractArray{T, N}, HL::HyperbolicLayer; set_grad::Bool=true) where {T, N}

    # Recompute forward states
    X_prev_in, X_curr_in, X_conv, X_relu = inverse(X_curr, X_new, HL; save=true)

    # Backpropagate data residual and compute gradients
    cdims = DCDims(X_curr, HL.W.data; stride=HL.stride, padding=HL.pad)
    ΔX_convT = copy(ΔX_new)
    ΔX_relu = -HL.α*conv(ΔX_convT, HL.W.data, cdims)
    ΔW = -HL.α*∇conv_filter(ΔX_convT, X_relu, cdims)

    ΔX_conv = ReLUgrad(ΔX_relu, X_conv)
    ΔX_curr += ∇conv_data(ΔX_conv, HL.W.data, cdims)
    ΔW += ∇conv_filter(X_curr, ΔX_conv, cdims)
    if length(size(X_curr)) == 4
        Δb = sum(ΔX_conv; dims=[1,2,4])[1,1,:,1]
    else
        Δb = sum(ΔX_conv; dims=[1,2,3,5])[1,1,1,:,1]
    end
    ΔX_curr += 2f0*ΔX_new
    ΔX_prev = -ΔX_new

    # Set gradients
    if set_grad
        HL.W.grad = ΔW
        HL.b.grad = Δb
    else
        Δθ = [Parameter(ΔW), Parameter(Δb)]
    end

    # Change dimensions
    if HL.action == 0
        ΔX_prev_in = identity(ΔX_prev)
        ΔX_curr_in = identity(ΔX_curr)
    elseif HL.action == -1
        ΔX_prev_in = invHaar_unsqueeze(ΔX_prev)
        ΔX_curr_in = invHaar_unsqueeze(ΔX_curr)
    elseif HL.action == 1
        ΔX_prev_in = Haar_squeeze(ΔX_prev)
        ΔX_curr_in = Haar_squeeze(ΔX_curr)
    else
        throw("Specified operation not defined.")
    end

    set_grad ? (return ΔX_prev_in, ΔX_curr_in, X_prev_in, X_curr_in) : (return ΔX_prev_in, ΔX_curr_in, Δθ, X_prev_in, X_curr_in)
end

## Jacobian utilities

# 2D
function jacobian(ΔX_prev_in::AbstractArray{T, N}, ΔX_curr_in::AbstractArray{T, N}, Δθ::Array{Parameter}, X_prev_in::AbstractArray{T, N}, X_curr_in::AbstractArray{T, N}, HL::HyperbolicLayer) where {T, N}

    # Change dimensions
    if HL.action == 0
        X_prev = identity(X_prev_in)
        X_curr = identity(X_curr_in)
        ΔX_prev = identity(ΔX_prev_in)
        ΔX_curr = identity(ΔX_curr_in)
    elseif HL.action == 1
        X_prev = invHaar_unsqueeze(X_prev_in)
        X_curr = invHaar_unsqueeze(X_curr_in)
        ΔX_prev = invHaar_unsqueeze(ΔX_prev_in)
        ΔX_curr = invHaar_unsqueeze(ΔX_curr_in)
    elseif HL.action == -1
        X_prev = Haar_squeeze(X_prev_in)
        X_curr = Haar_squeeze(X_curr_in)
        ΔX_prev = Haar_squeeze(ΔX_prev_in)
        ΔX_curr = Haar_squeeze(ΔX_curr_in)
    else
        throw("Specified operation not defined.")
    end

    cdims = DCDims(X_curr, HL.W.data; stride=HL.stride, padding=HL.pad)
    # Symmetric convolution w/ relu activation
    if length(size(X_curr)) == 4
        X_conv = conv(X_curr, HL.W.data, cdims) .+ reshape(HL.b.data, 1, 1, :, 1)
    else
        X_conv = conv(X_curr, HL.W.data, cdims) .+ reshape(HL.b.data, 1, 1, 1, :, 1)
    end
    ΔX_conv = conv(ΔX_curr, HL.W.data, cdims) .+ conv(X_curr, Δθ[1].data, cdims) .+ reshape(Δθ[2].data, 1, 1, :, 1)
    X_relu = ReLU(X_conv)
    ΔX_relu = ReLUgrad(ΔX_conv, X_conv)
    X_convT = -∇conv_data(X_relu, HL.W.data, cdims)
    ΔX_convT = -∇conv_data(ΔX_relu, HL.W.data, cdims)-∇conv_data(X_relu, Δθ[1].data, cdims)

    # Update
    X_new = 2f0*X_curr - X_prev + HL.α*X_convT
    ΔX_new = 2f0*ΔX_curr - ΔX_prev + HL.α*ΔX_convT

    return ΔX_curr, ΔX_new, X_curr, X_new
end

function adjointJacobian(ΔX_curr::AbstractArray{T, N}, ΔX_new::AbstractArray{T, N}, X_curr::AbstractArray{T, N}, X_new::AbstractArray{T, N}, HL::HyperbolicLayer) where {T, N}
    return backward(ΔX_curr, ΔX_new, X_curr, X_new, HL; set_grad=false)
end

## Other utils

# Clear gradients
function clear_grad!(HL::HyperbolicLayer)
    HL.W.grad = nothing
    HL.b.grad = nothing
end

# Get parameters
get_params(HL::HyperbolicLayer) = [HL.W, HL.b]


text:
@testset "Risk Measures" begin
    sample = [0:100...]

    @testset "VaR" begin
        @test VaR(sample,.9) ≈ 90
        @test VaR(sample,.9,rev=true) ≈ 10
        @test VaR(sample,.1) ≈ VaR(sample,1-.9)
        @test VaR(sample,.1) == ValueAtRisk(sample,.1)
    end

    @testset "CTE" begin
        @test CTE(sample,.9) ≈ sum(90:100) / length(90:100)
        @test CTE(sample,.1) ≈ sum(10:100) / length(10:100)
        @test CTE(sample,.15) >= CTE(sample,.1) # monotonicity
        @test CTE(sample,.9,rev=true) ≈ sum(0:10) / length(0:10)
        @test CTE(sample,.9) == ConditionalTailExpectation(sample,.9)
    end

end

ext: jl

max_stars_repo_name: alecloudenback/ActuaryUtilities

max_issues_repo_name: alecloudenback/ActuaryUtilities

max_forks_repo_name: alecloudenback/ActuaryUtilities


text:
"""
    statscholesky(xtx::Symmetric{T}, tol::Real=-1) where {T<:AbstractFloat}

Return a `CholeskyPivoted` object created from `xtx` where the pivoting scheme
retains the original order unless singularity is detected.  Columns that are
(computationally) linearly dependent on columns to their left are moved to the
right hand side in a left circular shift.
"""
function statscholesky(xtx::Symmetric{T}, tol::Real = -1) where {T<:AbstractFloat}
    n = size(xtx, 2)
    chpiv = cholesky(xtx, Val(true), tol = T(-1), check = false)
    chunp = cholesky(xtx, check = false)
    r = chpiv.rank
    piv = [1:n;]
    if r < n
        nleft = n
        while r < nleft
            k = chunp.info
            if k < nleft
                piv = piv[[1:k-1; k+1:n; k]]
                chunp = cholesky!(Symmetric(xtx[piv, piv]), check = false)
            end
            nleft -= 1
        end
    end
    for j = (r+1):n   # an MKL <-> OpenBLAS difference
        for i = (r+1):j
            chunp.factors[i, j] = zero(T)
        end
    end
    CholeskyPivoted(chunp.factors, chunp.uplo, piv, r, tol, chpiv.info)
end


text:
using LinearAlgebra: diagind

struct WideBandBath{M<:DiabaticModel,V<:AbstractVector,T} <: DiabaticFrictionModel
    model::M
    bathstates::V
    ρ::T
    function WideBandBath(model, bathstates)
        bathstates = austrip.(bathstates)
        ρ = length(bathstates) / (bathstates[end] - bathstates[begin])
        new{typeof(model),typeof(bathstates),typeof(ρ)}(model, bathstates, ρ)
    end
end

function WideBandBath(model::DiabaticModel; step, bandmin, bandmax)
    WideBandBath(model, range(bandmin, bandmax; step=step))
end

NQCModels.nstates(model::WideBandBath) = NQCModels.nstates(model.model) + length(model.bathstates) - 1
NQCModels.ndofs(model::WideBandBath) = NQCModels.ndofs(model.model)
NQCModels.nelectrons(model::WideBandBath) = fld(NQCModels.nstates(model), 2)
NQCModels.fermilevel(::WideBandBath) = 0.0

function NQCModels.potential!(model::WideBandBath, V::Hermitian, r::AbstractMatrix)

    Vsystem = NQCModels.potential(model.model, r)
    n = NQCModels.nstates(model.model)
    ϵ0 = Vsystem[1,1]

    # System states
    @views V[diagind(V)[begin:n-1]] .= Vsystem[diagind(Vsystem)[begin+1:end]] .- ϵ0

    # Bath states
    V[diagind(V)[n:end]] .= model.bathstates

    # Coupling
    for i=1:n-1
        V.data[i,n:end] .= Vsystem[i+1,1] / sqrt(model.ρ)
    end

    return V
end

function NQCModels.derivative!(model::WideBandBath, D::AbstractMatrix{<:Hermitian}, r::AbstractMatrix)

    Dsystem = NQCModels.derivative(model.model, r)
    n = NQCModels.nstates(model.model)

    for I in eachindex(Dsystem, D)
        ∂ϵ0 = Dsystem[I][1,1]

        # System states
        D_system_output = @view D[I][diagind(D[I])[begin:n-1]]
        D_system_input = @view Dsystem[I][diagind(Dsystem[I])[begin+1:end]]
        D_system_output .= D_system_input .- ∂ϵ0

        # Coupling
        for i=1:n-1
            D[I].data[i,n:end] .= Dsystem[I][i+1,1] / sqrt(model.ρ)
        end
    end

    return D
end

function NQCModels.state_independent_potential(model::WideBandBath, r)
    Vsystem = NQCModels.potential(model.model, r)
    return Vsystem[1,1]
end

function NQCModels.state_independent_derivative!(model::WideBandBath, ∂V::AbstractMatrix, r)
    Dsystem = NQCModels.derivative(model.model, r)
    for I in eachindex(∂V, Dsystem)
        ∂V[I] = Dsystem[I][1,1]
    end

    return ∂V
end


text:
using MeshCat, MeshCatMechanisms, RigidBodyDynamics
using FileIO, MeshIO, GeometryTypes, CoordinateTransformations, Rotations, Meshing
using Plots

include(joinpath(pwd(),"src/models/kuka/kuka_utils.jl"))

# Mini golf surface
t = range(-1.0, stop = 3.0, length = 100)
_k = 10.0
shift = 1.5
softplus(x) = log(1 + exp(_k * (x - shift))) / _k
dsoftplus(x) = 1.0 / (1.0 + exp(-1.0 * _k * (x - shift)))
plot(t, softplus.(t))

# ellipsoid
ellipsoid(x, y, rx, ry) = (x^2.0) / (rx^2) + (y^2) / (ry^2) - 1.0
∇ellipsoid(x,y,rx,ry) = [2.0 * x / (rx^2.0); 2.0 * y / (ry^2.0); 0.0]
get_y(x,rx,ry) = sqrt((1.0 - ((x)^2.0) / ((rx)^2.0)) * (ry^2))

ellipsoid(x,y,xc,yc,rx,ry) = ((x - xc)^2.0) / (rx^2.0) + ((y - yc)^2.0) / (ry^2.0) - 1.0
∇ellipsoid(x,y,xc,yc,rx,ry) = [2.0 * (x - xc) / (rx^2.0); 2.0*(y - yc) / (ry^2.0); 0.0]
get_y(x,xc,yc,rx,ry) = sqrt((1.0 - ((x - xc)^2.0) / ((rx)^2.0)) * (ry^2.0)) + yc

rx1 = 1.5
ry1 = 1.5

rx2 = 2.5
ry2 = 2.5

urdf_path = joinpath(pwd(), "src/models/kuka/temp/kuka.urdf")

kuka = MeshCatMechanisms.parse_urdf(urdf_path, remove_fixed_tree_joints = true)
kuka_visuals = MeshCatMechanisms.URDFVisuals(urdf_path)

state = MechanismState(kuka)
state_cache = StateCache(kuka)
result = DynamicsResult(kuka)
result_cache = DynamicsResultCache(kuka)

vis = Visualizer()
mvis = MechanismVisualizer(kuka, kuka_visuals, vis[:base])
open(vis)

nq = num_positions(kuka)
q_init = zeros(nq)
q_init[1] = 0
q_init[2] = pi/4

q_init[4] = -pi/2
q_init[5] = 0
q_init[6] = -pi/4
q_init[7] = 0
set_configuration!(state, q_init)
set_configuration!(mvis, q_init)
ee = findbody(kuka, "iiwa_link_7")
ee_point = Point3D(default_frame(ee), 0.374, 0.0, 0.0)
setelement!(mvis, ee_point, 0.01)
ee_jacobian_frame = ee_point.frame
ee_jacobian_path = path(kuka, root_body(kuka), ee)

world = root_frame(kuka)
ee_in_world = transform(state, ee_point, world).v
desired = Point3D(world, 0.5, 0.0, 0.0)
q_res1 = jacobian_transpose_ik!(state, ee, ee_point, desired,
    visualize = true, mvis = mvis)
q_res1 = Array(q_res1)
set_configuration!(mvis, q_res1)
set_configuration!(state, q_res1)
ee_in_world = transform(state, ee_point, world).v

desired = Point3D(world, 0.75, 0.0, 0.0)
q_res2 = jacobian_transpose_ik!(state, ee, ee_point, desired,
    visualize = true, mvis = mvis)
q_res2 = Array(q_res2)
set_configuration!(mvis, q_res2)
set_configuration!(state, q_res2)
ee_in_world = transform(state, ee_point, world).v

# Kuka iiwa arm parsed from URDF using RigidBodyDynamics.jl
# + 3D particle
struct KukaParticle{T} <: Model
	n::Int
	m::Int
	d::Int

	# particle
	mp::T
	rp::T

    qL::Vector{T}
    qU::Vector{T}

    state_cache1
    state_cache2
    state_cache3

    result_cache1
    result_cache2
    result_cache3

    world

	ee
    ee_point
    ee_jacobian_frame
    ee_jacobian_path

	μ_ee_p
	μ1
	μ2
	μ3

	nq
    nu
    nc
    nf
    nb

    idx_u
    idx_λ
    idx_b
    idx_ψ
    idx_η
    idx_s

end

# Dimensions
nq = 7 + 3 # configuration dim
nu = 7 # control dim
nc = 3 # number of contact points
nf = 4 # number of faces for friction cone
nb = nc * nf
ns = 1

n = 2 * nq
m_contact = nc + nb + nc + nb + 1
m = nu + m_contact
d = 0

idx_u = (1:nu)
idx_λ = nu .+ (1:nc)
idx_b = nu + nc .+ (1:nb)
idx_ψ = nu + nc + nb .+ (1:nc)
idx_η = nu + nc + nb + nc .+ (1:nb)
idx_s = nu + nc + nb + nc + nb .+ (1:ns)

mp = 0.05
rp = 0.05
μ_ee_p = 1.0
μ1 = 0.01
μ2 = 0.1
μ3 = 0.1

results_cache1 = DynamicsResultCache(kuka)
results_cache2 = DynamicsResultCache(kuka)
results_cache3 = DynamicsResultCache(kuka)

state_cache1 = StateCache(kuka)
state_cache2 = StateCache(kuka)
state_cache3 = StateCache(kuka)

function kuka_q(q)
	view(q, 1:7)
end

function particle_q(q)
	view(q, 8:10)
end

# Methods
function M_func(m::KukaParticle, q::AbstractVector{T}) where T
    state = m.state_cache3[T]
    result = m.result_cache3[T]
    set_configuration!(state, kuka_q(q))
    mass_matrix!(result.massmatrix, state)
    return cat(result.massmatrix, Diagonal(m.mp * ones(3)), dims = (1, 2))
end

function C_func(m::KukaParticle, qk, qn, h)
    if eltype(kuka_q(qk)) <: ForwardDiff.Dual
		return [dynamics_bias_qk(m, kuka_q(qk), kuka_q(qn), h)..., 0.0, 0.0, m.mp * 9.81]
	elseif eltype(kuka_q(qn)) <: ForwardDiff.Dual
		return [dynamics_bias_qn(m,kuka_q(qk), kuka_q(qn), h)..., 0.0 ,0.0, m.mp * 9.81]
	elseif typeof(h) <: ForwardDiff.Dual
		return [dynamics_bias_h(m, kuka_q(qk), kuka_q(qn), h)..., 0.0, 0.0, m.mp * 9.81]
	else
		return [_dynamics_bias(m, kuka_q(qk), kuka_q(qn), h)..., 0.0, 0.0, m.mp * 9.81]
	end
end

function _dynamics_bias(m::KukaParticle, qk, qn, h)
	T = eltype(m.qL)
	state = m.state_cache1[T]
	result = m.result_cache1[T]
	set_configuration!(state,qk)
	set_velocity!(state, (qn - qk) / h)
	dynamics_bias!(result, state)
    return result.dynamicsbias
end

function dynamics_bias_qk(m::KukaParticle, qk::AbstractVector{T}, qn, h) where T
    state = m.state_cache2[T]
    result = m.result_cache2[T]
    set_configuration!(state, qk)
    set_velocity!(state,(qn - qk) / h)
    dynamics_bias!(result, state)
    return result.dynamicsbias
end

function dynamics_bias_qn(m::KukaParticle, qk, qn::AbstractVector{T}, h) where T
    state = m.state_cache3[T]
    result = m.result_cache3[T]
    set_configuration!(state, qk)
    set_velocity!(state, (qn - qk) / h)
    dynamics_bias!(result, state)
    return result.dynamicsbias
end

function dynamics_bias_h(m::KukaParticle,qk,qn,h::T) where T
    state = m.state_cache1[T]
    result = m.result_cache1[T]
    set_configuration!(state, qk)
    set_velocity!(state,(qn - qk) / h)
    dynamics_bias!(result, state)
    return result.dynamicsbias
end

function B_func(m::KukaParticle, q)
    [Diagonal(ones(7)); zeros(3, 7)]
end

function friction_cone(model::KukaParticle, u)
	λ = u[model.idx_λ]
	b = u[model.idx_b]

    @SVector [model.μ_ee_p * λ[1] - sum(b[1:4]),#,
			  model.μ1 * λ[2] - sum(b[5:8]),
			  model.μ2 * λ[3] - sum(b[9:12])]
			  # model.μ3*λ[3] - sum(b[9:12])]
end

function maximum_dissipation(model::KukaParticle, x⁺, u, h)
	q3 = x⁺[model.nq .+ (1:model.nq)]
	q2 = x⁺[1:model.nq]

	ψ = u[model.idx_ψ]
    ψ_stack = [ψ[1] * ones(4);
			   ψ[2] * ones(4);
			   ψ[3] * ones(4)]

    η = u[model.idx_η]

    P_func(model, q3) * (q3 - q2) / h + ψ_stack - η
end

function fd(model::KukaParticle, x⁺, x, u, w, h, t)
	q3 = view(x⁺, model.nq .+ (1:model.nq))
	q2⁺ = view(x⁺, 1:model.nq)
	q2⁻ = view(x, model.nq .+ (1:model.nq))
	q1 = view(x, 1:model.nq)
	u_ctrl = view(u, model.idx_u)
	λ = view(u, model.idx_λ)
	b = view(u, model.idx_b)

    [q2⁺ - q2⁻;
    ((1.0 / h) * (M_func(model, q1) * (SVector{10}(q2⁺) - SVector{10}(q1))
    - M_func(model, q2⁺) * (SVector{10}(q3) - SVector{10}(q2⁺)))
    + B_func(model, q3) * SVector{7}(u_ctrl)
    + transpose(N_func(model, q3)) * SVector{3}(λ)
    + transpose(P_func(model, q3)) * SVector{12}(b)
    - h * C_func(model, q2⁺, q3, h))]
end

function kinematics_ee(m::KukaParticle, q::AbstractVector{T}) where T
    state = m.state_cache3[T]
    set_configuration!(state, kuka_q(q))
    return transform(state, m.ee_point, m.world).v
end

function ϕ_func(m::KukaParticle, q::AbstractVector{T}) where T
    # state = m.state_cache3[T]
    # set_configuration!(state,kuka_q(q))
	# p_ee = transform(state, m.ee_point, m.world).v

	p_ee = kinematics_ee(m, q)
	p_p = q[8:10]
	diff = (p_ee - p_p)
	d_ee_p = diff' * diff
    @SVector [d_ee_p - m.rp^2,#,
		      ellipsoid(p_p[1], p_p[2], 0.25, -1.5, 1.0, 1.0),
	          # # -1.0*ellipsoid(p_p[1],p_p[2],m.rx2,m.ry2),
	          p_p[3] - softplus(p_p[1])]
end

function N_func(m::KukaParticle, q::AbstractVector{T}) where T
    state = m.state_cache3[T]

    # N = zeros(T,m.nc,m.nq)

    set_configuration!(state, kuka_q(q))

    pj1 = PointJacobian(transform(state, m.ee_point, m.world).frame,
		zeros(T, 3, 7))

    ee_in_world = transform(state, m.ee_point, m.world)

    point_jacobian!(pj1, state, m.ee_jacobian_path,ee_in_world) #TODO confirm this is correct

    # N[1,:] = pj1.J[3,:]
    # N[2,:] = pj2.J[3,:]

	p_ee = ee_in_world.v
	p_p = q[8:10]
	diff = (p_ee - p_p)
	# d_ee_p = diff'*diff

    return [2.0 * diff' * [pj1.J[1:3,:] -1.0 * Diagonal(ones(3))];#;
			zeros(1, 7) transpose(∇ellipsoid(p_p[1], p_p[2], 0.25, -1.5, 1.0, 1.0));
			# # zeros(1,7) -1.0*transpose(∇ellipsoid(p_p[1],p_p[2],m.rx2,m.ry2));
			zeros(1, 7) -dsoftplus(p_p[1]) 0.0 1.0]
	# ϕ_tmp(y) = ϕ_func(model,y)
	# ForwardDiff.jacobian(ϕ_tmp,q)
end

# tmp(x) = ϕ_func(model,x)
# vec(ForwardDiff.jacobian(tmp,x1)) - vec(N_func(model,x1))

function P_func(m::KukaParticle, q::AbstractVector{T}) where T
    state = m.state_cache3[T]

	map = @SMatrix [1.0 0.0;
				    -1.0 0.0;
				    0.0 1.0;
				    0.0 -1.0]

    set_configuration!(state, kuka_q(q))

    pj1 = PointJacobian(transform(state, m.ee_point, m.world).frame,
		zeros(T, 3, 7))

    ee_in_world = transform(state, m.ee_point, m.world)

    point_jacobian!(pj1, state, m.ee_jacobian_path, ee_in_world) #TODO confirm this is correct

	y0 = [0.0; 1.0; 0.0]
	z0 = [0.0; 0.0; 1.0]
	p_p = particle_q(q)
    return [map * pj1.J[2:3,:] zeros(4,3);
			zeros(4,7) map * [transpose(z0); transpose(cross(z0, ∇ellipsoid(p_p[1], p_p[2], 0.25, -1.5, 1.0, 1.0)))];
			# # zeros(4,7) map*[transpose(z0); transpose(cross(z0,∇ellipsoid(p_p[1],p_p[2],m.rx2,m.ry2)))];
			zeros(4,7) map * [transpose(y0); transpose(cross(y0, [-1.0 * dsoftplus(p_p[1]); 0.0; 1.0]))]]

end

qL = -Inf * ones(nq)
qU = Inf * ones(nq)

model = KukaParticle(
	n, m, d,
	mp, rp,
	qL, qU,
    state_cache1, state_cache2, state_cache3,
    results_cache1, results_cache2, results_cache3,
    world,
	ee,
	ee_point,
	ee_jacobian_frame,
	ee_jacobian_path,
	μ_ee_p,
	μ1, μ2, μ3,
	nq, nu, nc, nf, nb,
	idx_u,
	idx_λ,
	idx_b,
	idx_ψ,
	idx_η,
	idx_s)

# Visualization
function visualize!(mvis, model::KukaParticle, q;
		verbose = false, r_ball = 0.1, Δt = 0.1)

	p_hole = [0.66; -3.0; 0.0]
	f = x -> x[3] - softplus(x[1])

	sdf = SignedDistanceField(f, Rect(Vec(-5, -10, -1), Vec(10, 10, 4)))
	mesh = HomogenousMesh(sdf, MarchingTetrahedra())
	setobject!(vis["slope"], mesh,
			   MeshPhongMaterial(color=RGBA{Float32}(86/255, 125/255, 70/255, 1.0)))
	settransform!(vis["slope"], compose(Translation(0.0, 5.0, 0.0)))

	circle1 = Cylinder(Point3f0(0,0,0), Point3f0(0,0,0.25), convert(Float32,1.0 - r_ball))
		setobject!(vis["circle1"], circle1,
		MeshPhongMaterial(color = RGBA(86/255, 125/255, 20/255, 1.0)))
	settransform!(vis["circle1"], compose(Translation(0.25, -1.5, 0.0)))

	setobject!(vis["ball"], Sphere(Point3f0(0),
				convert(Float32, r_ball)),
				MeshPhongMaterial(color = RGBA(1,1,1,1.0)))

	settransform!(vis["ball"], compose(Translation(0.66,3.0,0.0)))

	hole = Cylinder(Point3f0(0,0,0), Point3f0(0,0,0.01),
		convert(Float32, r_ball * 1.5))
		setobject!(vis["hole"], hole,
		MeshPhongMaterial(color = RGBA(0,0,0,1.0)))
	settransform!(vis["hole"], compose(Translation(p_hole...)))

	club = Cylinder(Point3f0(0,0,0), Point3f0(0,0,0.374), convert(Float32,0.025))
   	setobject!(vis["club"], club, MeshPhongMaterial(color = RGBA(0,0,0,1.0)))

	state = model.state_cache1[Float64]

	ee = findbody(kuka, "iiwa_link_7")
	ee_body = Point3D(default_frame(ee), 0.0, 0.0, 0.0)
	ee_end = Point3D(default_frame(ee), 0.374, 0.0, 0.0)
	ee_body_jacobian_frame = ee_body.frame
	ee_body_jacobian_path = path(kuka, root_body(kuka), ee)
	ee_end_jacobian_frame = ee_end.frame
	ee_end_jacobian_path = path(kuka, root_body(kuka), ee)

	world = root_frame(kuka)
	ee_body_in_world = transform(state, ee_body, world).v
	ee_end_in_world = transform(state, ee_end, world).v

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

	T = length(q)
    for t = 1:T
        q_kuka = kuka_q(q[t])
		q_particle = particle_q(q[t])
		set_configuration!(state, kuka_q(q[t]))
		ee_body_in_world = transform(state, ee_body, world).v
		ee_end_in_world = transform(state, ee_end, world).v

        MeshCat.atframe(anim,t) do
			set_configuration!(mvis,q_kuka)
            settransform!(vis["ball"], compose(Translation(q_particle + [0.0;0.0;0.5 * r_ball]), LinearMap(RotZ(0))))
			settransform!(vis["club"], cable_transform(ee_body_in_world, ee_end_in_world))

			if norm(particle_q(q[t])[1:2] - p_hole[1:2]) < 1.0e-2
				setvisible!(vis["ball"], false)
			else
				setvisible!(vis["ball"], true)
			end
		end
    end

    MeshCat.setanimation!(vis, anim)
end

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
using Test, Mineos, SeisModels

@testset "Eigenfrequencies" begin
    let freqs = eigenfrequencies(LinearLayeredModel(PREM), lmax=128, nmax=10)
        @test length(freqs) == 4232
        @test freqs[0,:S,0] ≈ 0.814338 atol=1e-6
        @test freqs[0,:T,7] ≈ 1.22036 atol=1e-5
        @test freqs[8,:C,44] ≈ 40.28792 atol=1e-5
    end
end


text:
export @gen, @param, @trace

import MacroTools

const DSL_STATIC_ANNOTATION = :static
const DSL_ARG_GRAD_ANNOTATION = :grad
const DSL_RET_GRAD_ANNOTATION = :grad
const DSL_TRACK_DIFFS_ANNOTATION = :diffs
const DSL_NO_JULIA_CACHE_ANNOTATION = :nojuliacache

struct Argument
    name::Symbol
    typ::Union{Symbol,Expr}
    annotations::Set{Symbol}
    default::Union{Some{Any}, Nothing}
end

Argument(name, typ) = Argument(name, typ, Set{Symbol}(), nothing)
Argument(name, typ, annotations) = Argument(name, typ, annotations, nothing)

function parse_annotations(annotations_expr)
    annotations = Set{Symbol}()
    if isa(annotations_expr, Symbol)
        push!(annotations, annotations_expr)
    elseif isa(annotations_expr, Expr) && annotations_expr.head == :tuple
        for annotation in annotations_expr.args
            if !isa(annotation, Symbol)
                error("syntax error in annotations_expr at $annotation")
            else
                push!(annotations, annotation)
            end
        end
    else
        error("syntax error in annotations at $annotations")
    end
    annotations
end

function parse_arg(expr)
    if isa(expr, Symbol)
        # x
        arg = Argument(expr, :Any)
    elseif isa(expr, Expr) && expr.head == :(::)
        # x::Int
        arg = Argument(expr.args[1], expr.args[2])
    elseif isa(expr, Expr) && expr.head == :kw
        # x::Int=1
        sub_arg = parse_arg(expr.args[1])
        default = Some(expr.args[2])
        arg = Argument(sub_arg.name, sub_arg.typ, Set{Symbol}(), default)
    elseif isa(expr, Expr) && expr.head == :call
        # (grad,foo)(x::Int)
        annotations_expr = expr.args[1]
        sub_arg = parse_arg(expr.args[2])
        annotations = parse_annotations(annotations_expr)
        arg = Argument(sub_arg.name, sub_arg.typ, annotations, sub_arg.default)
    else
        dump(expr)
        error("syntax error in gen function argument at $expr")
    end
    arg
end

include("dynamic.jl")
include("static.jl")

function desugar_tildes(expr)
    MacroTools.postwalk(expr) do e
        # Replace tilde statements with :gentrace expressions
        if MacroTools.@capture(e, {*} ~ rhs_call)
            Expr(:gentrace, rhs, nothing)
        elseif MacroTools.@capture(e, {addr_} ~ rhs_call)
            Expr(:gentrace, rhs, Some(addr))
        elseif MacroTools.@capture(e, lhs_Symbol ~ rhs_call)
            addr = QuoteNode(lhs)
            Expr(:(=), lhs, Expr(:gentrace, rhs, Some(addr)))
        elseif MacroTools.@capture(e, lhs_ ~ rhs_call)
            error("Syntax error: Invalid left-hand side: $(e)." *
                  "Only a variable or address can appear on the left of a `~`.")
        elseif MacroTools.@capture(e, lhs_ ~ rhs_)
            error("Syntax error: Invalid right-hand side in: $(e)")
        else
            e
        end
    end
end

function extract_quoted_exprs(expr)
    quoted_exprs = []
    expr = MacroTools.prewalk(expr) do e
        if MacroTools.@capture(e, :(quoted_)) && !isa(quoted, Symbol)
            push!(quoted_exprs, e)
            Expr(:placeholder, length(quoted_exprs))
        else
            e
        end
    end
    return expr, quoted_exprs
end

function insert_quoted_exprs(expr, quoted_exprs)
    expr = MacroTools.prewalk(expr) do e
        if MacroTools.@capture(e, p_placeholder)
            idx = p.args[1]
            quoted_exprs[idx]
        else
            e
        end
    end
    return expr
end

function preprocess_body(expr, __module__)
    # Expand all macros relative to the calling module
    expr = macroexpand(__module__, expr)
    # Protect quoted expressions from pre-processing by extracting them
    expr, quoted_exprs = extract_quoted_exprs(expr)
    # Desugar tilde calls to :gentrace expressions
    expr = desugar_tildes(expr)
    # Reinsert quoted expressions after pre-processing
    expr = insert_quoted_exprs(expr, quoted_exprs)
    return expr
end

function parse_gen_function(ast, annotations, __module__)
    ast = MacroTools.longdef(ast)
    if ast.head != :function
        error("syntax error at $ast in $(ast.head)")
    end
    if length(ast.args) != 2
        error("syntax error at $ast in $(ast.args)")
    end
    signature = ast.args[1]
    if signature.head == :(::)
        (call_signature, return_type) = signature.args
    elseif signature.head == :call
        (call_signature, return_type) = (signature, :Any)
    else
        error("syntax error at $(signature)")
    end
    body = preprocess_body(ast.args[2], __module__)
    name = call_signature.args[1]
    args = map(parse_arg, call_signature.args[2:end])
    static = DSL_STATIC_ANNOTATION in annotations
    if static
        make_static_gen_function(name, args, body, return_type, annotations)
    else
        make_dynamic_gen_function(name, args, body, return_type, annotations)
    end
end

macro gen(annotations_expr, ast::Expr)
    # parse the annotations
    annotations = parse_annotations(annotations_expr)
    # parse the function definition
    parse_gen_function(ast, annotations, __module__)
end

macro gen(ast::Expr)
    parse_gen_function(ast, Set{Symbol}(), __module__)
end

macro trace(expr::Expr)
    return Expr(:gentrace, esc(expr), nothing)
end

macro trace(expr::Expr, addr)
    return Expr(:gentrace, esc(expr), Some(addr))
end

macro param(expr::Expr)
    if (expr.head != :(::)) error("Syntax in error in @param at $(expr)") end
    name, type = expr.args
    return Expr(:genparam, esc(name), esc(type))
end

macro param(sym::Symbol)
    return Expr(:genparam, esc(sym), esc(:Any))
end


text:
export draw_curve

"""
`draw_curve(plist)` draws a closed curve through points in the
`plist` which is a 1-dimensional list of complex numbers.
Use `draw_curve(plist,false)` for an open curve.
"""
function draw_curve(pts::Array{Complex{T},1}, closed::Bool=true; opts...) where T
    x = real.(pts)
    y = imag.(pts)
    np = length(pts)

    kind = closed ? :closed : :open
    maxt = closed ? np+1 : np
    Sx = Spline(x,kind)
    Sy = Spline(y,kind)

    fx = funk(Sx)
    fy = funk(Sy)

    plot!(fx,fy,1,maxt;opts...)
end


text:
module MiniLoggersTest

using ReTest

include("test01_tokenizer.jl")
include("test02_loggerformat.jl")

end # module

MiniLoggersTest.runtests()


text:
using Test
using Dates
using Pretend
using DataFrames

using Discord:
    Client, DiscordChannel, Message, User

using HoJBot:
    IgHolding, IgPortfolio, IgUserError, PrettyView, SimpleView,
    current_date, discord_channel, discord_reply, discord_upload_file,
    recent_date_range, retrieve, retrieve_users, upload_file,
    ig_affirm_non_player, ig_affirm_player,
    ig_buy, ig_cash_entry, ig_chart, ig_count_shares, ig_execute,
    ig_file_path, ig_get_real_time_quote, ig_grouped_holdings, ig_hey,
    ig_historical_prices, ig_holdings_data_frame, ig_is_player,
    ig_load_all_portfolios, ig_load_portfolio, ig_mark_to_market!,
    ig_mark_to_market_portfolio, ig_ranking_table,
    ig_reformat_view!, ig_remove_game, ig_save_portfolio,
    ig_sell, ig_start_new_game, ig_value_all_portfolios, ig_view_table

Pretend.activate()

const USER_ID = UInt64(0)
const USER_ID2 = UInt64(1)
const USER_ID3 = UInt64(2)
const USER_NAME = "Joe"

function test_ig_cases(f, name)
    @testset "$name" begin
        try
            @test_nowarn ig_start_new_game(USER_ID)
            @test_nowarn ig_start_new_game(USER_ID2)
            f()
        finally
            @test_nowarn ig_remove_game(USER_ID)
            @test_nowarn ig_remove_game(USER_ID2)
        end
    end
end

function clean_ig_game_files()
    for id in (USER_ID, USER_ID2, USER_ID3)
        path = ig_file_path(id)
        isfile(path) && rm(path)
    end
end

# Mock functions
mock_ig_get_quote_100(symbol::AbstractString) = 100.0
mock_currrent_date() = Date(2021, 1, 11)
mock_reply(c::Client, m::Message, s::AbstractString) = nothing
mock_upload_file(c::Client, ch::DiscordChannel, filename::AbstractString; kwargs...) = nothing
mock_channel(c::Client, id::UInt64) = DiscordChannel(; id = 0x00, type = 0x00)

function mock_retrieve_users(c::Client, ids::Vector{UInt64})
    return Dict(id => User(; id=id, username="$id") for id in ids)
end

# Mock objects
mocked_client() = Client("hey")

@testset "Investment Game" begin

    # ensure everything is clean before we start
    clean_ig_game_files()

    test_ig_cases("Life cycle") do
        @test ig_is_player(USER_ID) == true
        @test ig_hey(USER_NAME, "wat") isa AbstractString
        @test_nowarn ig_affirm_player(USER_ID)
        @test_throws IgUserError ig_affirm_non_player(USER_ID)
    end

    test_ig_cases("Load/save") do
        # load
        pf = @test_nowarn ig_load_portfolio(USER_ID)
        @test pf.cash > 0
        @test isempty(pf.holdings)

        # save with holdings
        push!(pf.holdings, IgHolding("AAPL", 1, today(), 100))
        @test_nowarn ig_save_portfolio(USER_ID, pf)
        pf = @test_nowarn ig_load_portfolio(USER_ID)
        @test length(pf.holdings) == 1

        @test_nowarn ig_save_portfolio(USER_ID2, pf)
        pfs = ig_load_all_portfolios()
        @test length(pfs) >= 2   # two test pf's, plus whatever is there for others
        @test haskey(pfs, USER_ID)
        @test haskey(pfs, USER_ID2)
    end

    test_ig_cases("Buy/sell") do
        apply(ig_get_real_time_quote => mock_ig_get_quote_100) do
            # Buy something, portfolio should be updated
            executed_price = @test_nowarn ig_buy(USER_ID, "AAPL", 50)
            @test executed_price == 100.0

            # Load it again and check
            pf = ig_load_portfolio(USER_ID)
            @test length(pf.holdings) == 1
            @test pf.holdings[1].symbol == "AAPL"
            @test pf.holdings[1].shares == 50
            @test pf.cash == 1_000_000 - 50 * 100.0

            # Sell more than what you have
            @test_throws IgUserError ig_sell(USER_ID, "AAPL", 70)

            # Partial sell
            executed_price = @test_nowarn ig_sell(USER_ID, "AAPL", 30)
            @test executed_price == 100.0

            # What's left?
            pf = ig_load_portfolio(USER_ID)
            @test length(pf.holdings) == 1
            @test pf.holdings[1].symbol == "AAPL"
            @test pf.holdings[1].shares == 20

            # Sell completely
            ig_sell(USER_ID, "AAPL", 20)
            pf = ig_load_portfolio(USER_ID)
            @test length(pf.holdings) == 0
            @test pf.cash == 1_000_000

            # buy 3 lots relief
            ig_buy(USER_ID, "IBM", 50)
            ig_buy(USER_ID, "IBM", 50)
            ig_buy(USER_ID, "IBM", 50)
            ig_sell(USER_ID, "IBM", 80)
            pf = ig_load_portfolio(USER_ID)
            @test pf.cash == 1_000_000 - 70 * 100.0
            @test length(pf.holdings) == 2
            @test sort([h.shares for h in pf.holdings]) == [20, 50]
        end
    end

    test_ig_cases("Views") do
        # empty portfolio with pre-filled cash only
        df = @test_nowarn ig_mark_to_market_portfolio(USER_ID)
        @test nrow(df) == 1
        @test df.symbol[1] == "CASH:USD"

        # setup
        ig_buy(USER_ID, "AAPL", 100, 120)
        ig_buy(USER_ID, "IBM", 50, 130)
        ig_buy(USER_ID, "IBM", 50, 150)
        pf = ig_load_portfolio(USER_ID)

        # counting
        @test ig_count_shares(pf, "AAPL") == 100
        @test ig_count_shares(pf, "IBM") == 100

        # convert to holdings data frame
        df = @test_nowarn ig_holdings_data_frame(pf)
        @test df isa AbstractDataFrame
        @test sort(propertynames(df)) == sort([:symbol, :shares, :purchase_price, :purchase_date])

        # grouped with average pricing
        gdf = ig_grouped_holdings(df)
        @test hasproperty(gdf, :symbol) == true
        @test hasproperty(gdf, :shares) == true
        @test hasproperty(gdf, :purchase_price) == true
        @test gdf[gdf.symbol .== "IBM", :].purchase_price[1] == 140

        # mark to market
        df2 = copy(df)
        apply(ig_get_real_time_quote => mock_ig_get_quote_100) do
            @test_nowarn ig_mark_to_market!(df2)
            @test hasproperty(df2, :current_price)
            @test hasproperty(df2, :market_value)
            @test unique(df2.current_price) == [100]

            df3 = @test_nowarn ig_mark_to_market_portfolio(USER_ID)
            @test df3.current_price .* df3.shares == df3.market_value

            # compatible column names
            @test Set(propertynames(ig_cash_entry(pf))) == Set(propertynames(df3))

            # translation to string
            @test_nowarn ig_view_table(PrettyView(), df2)

            # SimpleView requires specific format
            @test_nowarn ig_reformat_view!(df3)
            @test_nowarn ig_view_table(SimpleView(), df3)
        end
    end

    test_ig_cases("Ranking") do
        ig_buy(USER_ID, "AAPL", 100, 60)
        ig_buy(USER_ID, "IBM", 100, 60)
        ig_buy(USER_ID2, "AAPL", 100, 120)  # bought at higher price i.e. less cash remaining
        ig_buy(USER_ID2, "IBM", 100, 120)
        apply(
            ig_get_real_time_quote => mock_ig_get_quote_100,
            retrieve_users => mock_retrieve_users,
        ) do
            valuations = @test_nowarn ig_value_all_portfolios()
            # ranking order: highest valuation is at the top
            @test valuations[1].id == USER_ID
            @test valuations[2].id == USER_ID2

            client = mocked_client()
            rt = @test_nowarn ig_ranking_table(client)
            @test rt[1, :player] == "$USER_ID"
            @test rt[2, :player] == "$USER_ID2"
        end
    end

    test_ig_cases("Executors") do
        c = Client("test")
        m = Message(; id = 0x00, channel_id = 0x00)
        u = User(; id = USER_ID3, username = USER_NAME)
        apply(
            discord_channel => mock_channel,
            discord_reply => mock_reply,
            discord_upload_file => mock_upload_file,
            ig_get_real_time_quote => mock_ig_get_quote_100,
            retrieve_users => mock_retrieve_users,
        ) do
            @test_nowarn ig_execute(c, m, u, Val(Symbol("start-game")), [])
            @test_nowarn ig_execute(c, m, u, Val(Symbol("abandon-game")), [])
            @test_nowarn ig_execute(c, m, u, Val(Symbol("abandon-game-really")), [])

            @test_nowarn ig_execute(c, m, u, Val(Symbol("start-game")), [])
            @test_nowarn ig_execute(c, m, u, Val(:view), [])
            @test_nowarn ig_execute(c, m, u, Val(:view), ["simple"])
            @test_nowarn ig_execute(c, m, u, Val(:buy), ["100", "aapl"])
            @test_nowarn ig_execute(c, m, u, Val(:sell), ["100", "aapl"])
            @test_nowarn ig_execute(c, m, u, Val(Symbol("abandon-game-really")), [])

            @test_nowarn ig_execute(c, m, u, Val(:quote), ["aapl"])
            @test_nowarn ig_execute(c, m, u, Val(:chart), ["aapl"])
            @test_nowarn ig_execute(c, m, u, Val(:chart), ["aapl", "100d"])

            @test_nowarn ig_execute(c, m, u, Val(:rank), [])
            @test_nowarn ig_execute(c, m, u, Val(:rank), ["10"])
        end
    end

    @testset "Misc" begin
        apply(current_date => mock_currrent_date) do
            r = @test_nowarn recent_date_range(Day(10))
            @test r isa Tuple{Date,Date}
            @test r[1] == Date(2021, 1, 1)
            @test r[2] == Date(2021, 1, 11)
        end
    end

    @testset "Charting" begin
        dates = collect(Date(2021, 1, 1):Day(1):Date(2021,12,31))
        values = collect(1:365)
        @test_nowarn ig_chart("AAPL", dates, values)
    end

    @testset "Integration tests" begin
        from_date, to_date = Date(2021, 1, 1), Date(2021, 1, 31)
        @test_nowarn ig_historical_prices("AAPL", from_date, to_date)
        @test_throws IgUserError ig_historical_prices("BADSYMBOL", from_date, to_date)
    end

    clean_ig_game_files()
end


text:
using Distributions, UncertainData, Combinatorics

# Test all combinations of different types of uncertain values
M = MixtureModel([Normal(3, 0.2), Normal(2, 1)])

r1 = UncertainValue(Normal, rand(), rand())
r2 = UncertainValue(rand(M, 10000))
r3 = UncertainValue(Normal, rand(Normal(4, 3.2), 10000))

uvals = [r1; r2; r3]

a = UncertainDataset(uvals)
b = UncertainDataset(uvals[end:-1:1])
N = length(a)

###########
# Addition
###########
@test a + b isa UncertainDataset
@test 2 + b isa UncertainDataset
@test [2 for i = 1:N] + b isa UncertainDataset
@test 2.2 + b isa UncertainDataset
@test [2.2 for i = 1:N] + b isa UncertainDataset
@test a + [2 for i = 1:N] isa UncertainDataset
@test a + [2.2 for i = 1:N] isa UncertainDataset

#############
# Subtraction
#############
@test a - b isa UncertainDataset
@test 2 - b isa UncertainDataset
@test [2 for i = 1:N] - b isa UncertainDataset
@test 2.2 - b isa UncertainDataset
@test [2.2 for i = 1:N] - b isa UncertainDataset
@test a - [2 for i = 1:N] isa UncertainDataset
@test a - [2.2 for i = 1:N] isa UncertainDataset

################
# Multiplication
################
@test a * b isa UncertainDataset
@test 2 * b isa UncertainDataset
@test [2 for i = 1:N] * b isa UncertainDataset
@test 2.2 * b isa UncertainDataset
@test [2.2 for i = 1:N] * b isa UncertainDataset
@test a * [2 for i = 1:N] isa UncertainDataset
@test a * [2.2 for i = 1:N] isa UncertainDataset

################
# Division
################
@test a * b isa UncertainDataset
@test 2 * b isa UncertainDataset
@test [2 for i = 1:N] * b isa UncertainDataset
@test 2.2 * b isa UncertainDataset
@test [2.2 for i = 1:N] * b isa UncertainDataset
@test a * [2 for i = 1:N] isa UncertainDataset
@test a * [2.2 for i = 1:N] isa UncertainDataset

################
# Exponentiation
################
@test a * b isa UncertainDataset
@test 2 * b isa UncertainDataset
@test [2 for i = 1:N] * b isa UncertainDataset
@test 2.2 * b isa UncertainDataset
@test [2.2 for i = 1:N] * b isa UncertainDataset
@test a * [2 for i = 1:N] isa UncertainDataset
@test a * [2.2 for i = 1:N] isa UncertainDataset


text:
@testset "Model manipulation" begin
    m1 = Metabolite("m1")
    m2 = Metabolite("m2")
    m3 = Metabolite("m3")
    m4 = Metabolite("m4")
    mets = [m1, m2, m3, m4]
    m5 = Metabolite("m5")
    m6 = Metabolite("m6")
    m7 = Metabolite("m7")

    g1 = Gene("g1")
    g2 = Gene("g2")
    g3 = Gene("g3")
    g4 = Gene("g4")
    genes = [g1, g2, g3, g4]
    g5 = Gene("g5")
    g6 = Gene("g6")
    g7 = Gene("g7")

    r1 = Reaction("r1", Dict(m1.id => -1.0, m2.id => 1.0), :forward)
    r2 = Reaction("r2", Dict(m2.id => -2.0, m3.id => 1.0), :bidirectional)
    r2.grr = [["g2"], ["g1", "g3"]]
    r3 = Reaction("r3", Dict(m1.id => -1.0, m4.id => 2.0), :reverse)
    r4 = Reaction("r4", Dict(m1.id => -5.0, m4.id => 2.0), :reverse)
    r5 = Reaction("r5", Dict(m1.id => -11.0, m4.id => 2.0, m3.id => 2.0), :reverse)

    rxns = [r1, r2]

    model = StandardModel()
    model.id = "model"
    model.reactions = OrderedDict(r.id => r for r in rxns)
    model.metabolites = OrderedDict(m.id => m for m in mets)
    model.genes = OrderedDict(g.id => g for g in genes)

    ### reactions
    add!(model, [r3, r4])
    @test length(model.reactions) == 4

    add!(model, r5)
    @test length(model.reactions) == 5

    rm!(Reaction, model, ["r5", "r4"])
    @test length(model.reactions) == 3

    rm!(Reaction, model, "r1")
    @test length(model.reactions) == 2

    ### metabolites
    add!(model, [m5, m6])
    @test length(model.metabolites) == 6

    add!(model, m7)
    @test length(model.metabolites) == 7

    rm!(Metabolite, model, ["m5", "m4"])
    @test length(model.metabolites) == 5

    rm!(Metabolite, model, "m1")
    @test length(model.metabolites) == 4

    ### genes
    add!(model, [g5, g6])
    @test length(model.genes) == 6

    add!(model, g7)
    @test length(model.genes) == 7

    rm!(Gene, model, ["g5", "g4"])
    @test length(model.genes) == 5

    rm!(Gene, model, "g1")
    @test length(model.genes) == 4
end


text:
#
# Jerzy Wroczyński (nr. indeksu 250075)
# Obliczenia Naukowe, Laboratorium
# Lista 3, Zadania 1.–3.
#

module IterativeMethods

FL = Float64
MAX_ITERATIONS = 10^9
NEAR_ZERO = 2^(-2.0^32)

"""
### Finds the solution to the equation `f(x) = 0` using the bisection method.

### Arguments
- `f` — the function to find a root of
- `a` — beginning of the interval
- `b` — ending of the interval
- `delta` — the precision of the solution in terms of the X axis
- `epsilon` — the precision of the solution in terms of the Y axis

### Return value: a tuple of four values
- `[0]` — the approximation of the root of the `f` function
- `[1]` — the actual value of the `f` function in `[0]`
- `[2]` — number of iterations performed
- `[3]` — error indicator, possible values:
    - `0` — the method was convergent in this case
    - `1` — the `f` function doesn’t change its sign inside the given interval (doesn’t cut through the X axis)

"""
function mbisekcji(f, a::FL, b::FL, delta::FL, epsilon::FL)

    # begin with the starting values
    u = f(a)
    v = f(b)
    # the length of the section
    e = b - a

    if sign(u) == sign(v)
        # `f` doesn’t change its sign between `a` and `b`
        return (a, u, 0, 1);
    end

    # predefine some variables so they are accessible outside the loop
    k = 0
    c = 0
    w = 0
    # iterate
    for k in 0:MAX_ITERATIONS

        # divide into two sections — divide the length of the section by two
        e = e/2
        # pick the middle point
        c = a + e
        w = f(c)

        # if the solution is sufficient with respect to the given precision parameters terminate this method and return the output values
        if abs(e) < delta || abs(w) < epsilon
            return (c, w, k, 0)
        end

        # choose that section that contains the solution
        if sign(w) != sign(u)
            b = c
            v = w
        else
            a = c
            u = w
        end

    end

    # somehow the function exceeded MAX_ITERATIONS, however it seems so unlikely that a separate error code for this edge case would be an excess
    return (c, w, k, 0)

end

"""
### Finds the solution to the equation `f(x) = 0` using the Newton’s method.

### Arguments
- `f` — the function to find a root of
- `pf` — the derivative of the `f` function
- `x0` — first approximation
- `delta` — the precision of the solution in terms of the X axis
- `epsilon` — the precision of the solution in terms of the Y axis
- `maxint` — maximum iterations this function can perform

### Return value: a tuple of four values
- `[0]` — the approximation of the root of the `f` function
- `[1]` — the actual value of the `f` function in `[0]`
- `[2]` — number of iterations performed
- `[3]` — error indicator, possible values:
    - `0` — the method was convergent in this case
    - `1` — given solution is not within the given precision parameters because of not enough iterations
    - `2` — the derivative function was near zero

"""
function mstycznych(f, pf, x0::FL, delta::FL, epsilon::FL, maxint::Int)

    v = f(x0)

    if abs(v) < epsilon
        # given x₀ is already within given precision
        return (x0, v, 0, 0)
    end

    # iterate until the MAX_ITERATIONS is reached or the precision is satisfied
    xn = x0
    xn1 = x0

    # predefine the iterator so it would be accessible outside the loop
    k = 0
    # iterate
    for k in 1:maxint

        # perform the recursive step
        xn1 = xn - (v/pf(xn))
        v = f(xn1)

        # check the derivative
        if abs(pf(xn)) < NEAR_ZERO || isinf(abs(pf(xn)))
            # the derivative function is near zero which could lead to a not very precise solution
            return (xn, f(xn), k, 2)
        end

        # if the solution is within given precision return the output values
        if abs(xn1 - xn) < delta || abs(v) < epsilon
            return (xn1, v, k, 0)
        end

        # otherwise continue iterating
        xn = xn1

    end

    # not enough iterations to satisfy given precision
    return (xn1, v, k, 1)

end

"""
### Finds the solution to the equation `f(x) = 0` using the secant method.

### Arguments
- `f` — the function to find a root of
- `x0` — the first approximation
- `x1` — the second approximation
- `delta` — the precision of the solution in terms of the X axis
- `epsilon` — the precision of the solution in terms of the Y axis
- `maxint` — maximum iterations this function can perform

### Return value: a tuple of four values
- `[0]` — the approximation of the root of the `f` function
- `[1]` — the actual value of the `f` function in `[0]`
- `[2]` — number of iterations performed
- `[3]` — error indicator, possible values:
    - `0` — the method was convergent in this case
    - `1` — given solution is not within the given precision parameters because of not enough iterations

"""
function msiecznych(f, x0::FL, x1::FL, delta::FL, epsilon::FL, maxint::Int)

    # calc the first values
    v = f(x0)
    w = f(x1)
    a = x0
    b = x1

    # predefine the iterator so it would be accessible outside the loop
    k = 0
    #iterate
    for k in 0:maxint

        if abs(v) > abs(w)
            # swap the values, if necessary to maintain the order such that the latter value (`f(b)`) is equal or bigger than former (`f(a)`)
            # because of that the absolute value of the function will not rise — and that’s what we want since we’re calculating the root of the `f` function
            a,b = b,a
            v,w = w,v
        end

        # calculate the pseudo-derivative
        d = (b-a)/(w-v)
        # shift the values
        b = a
        w = v
        # calc the next value
        a = a - d * v
        v = f(a)

        if abs(v) < epsilon || abs(b-a) < delta
            # given solution is within given precision
            return (a, v, k, 0)
        end

    end

    # not enough iterations to give an accurate solution (not within given precision)
    return (a, v, k, 1)

end

end

max_stars_repo_name: jerry-sky/academic-notebook

max_issues_repo_name: jerry-sky/academic-notebook

max_forks_repo_name: jerry-sky/academic-notebook


text:
@testset "Takacs B&B vertical single large segment" begin

pressurecorrelation = BeggsAndBrill
p_initial = 346.6 - 14.65
dp_est = 100
t_avg = 107.2
md_initial = 0
md_end = 3700
tvd_initial = 0
tvd_end = 3700
inclination = 0
id = 2.441
roughness = 0.0003
q_o = 375
q_w = 0
GLR = 480
APIoil = 41.06
sg_water = 1
sg_gas = 0.916
molFracCO2 = 0
molFracH2S = 0
pseudocrit_pressure_correlation = HankinsonWithWichertPseudoCriticalPressure
pseudocrit_temp_correlation = HankinsonWithWichertPseudoCriticalTemp
Z_correlation = PapayZFactor
gas_viscosity_correlation = LeeGasViscosity
solutionGORcorrelation = StandingSolutionGOR
bubblepoint = 3000
oilVolumeFactor_correlation = StandingOilVolumeFactor
waterVolumeFactor_correlation = GouldWaterVolumeFactor
dead_oil_viscosity_correlation = GlasoDeadOilViscosity
live_oil_viscosity_correlation = ChewAndConnallySaturatedOilViscosity
frictionfactor = SerghideFrictionFactor
error_tolerance = 1.0 #psi

P_pc = pseudocrit_pressure_correlation(sg_gas, molFracCO2, molFracH2S)
    _, T_pc, _ = pseudocrit_temp_correlation(sg_gas, molFracCO2, molFracH2S)

ΔP_est = PressureDrop.calculate_pressuresegment(pressurecorrelation, p_initial, dp_est, t_avg,
                                md_end - md_initial, tvd_end - tvd_initial, inclination, true, id, roughness,
                                q_o, q_w, GLR, GLR, APIoil, sg_water, sg_gas,
                                Z_correlation, P_pc, T_pc,
                                gas_viscosity_correlation, solutionGORcorrelation, bubblepoint, oilVolumeFactor_correlation, waterVolumeFactor_correlation,
                                dead_oil_viscosity_correlation, live_oil_viscosity_correlation, frictionfactor, error_tolerance)

@test ΔP_est ≈ 3770 * (0.170 + 0.312) / 2 atol = 100 #order of magnitude test based on average between example points in Takacs (50).

end #testset Takacs B&B vertical

@testset "IHS Cleveland 6 - B&B full wellbore" begin

#%% end to end test

testpath = joinpath(dirname(dirname(pathof(PressureDrop))), "test/testdata/Cleveland_6/Test_survey_Cleveland_6.csv")

testwell = read_survey(path = testpath, id_included = false, maxdepth = 10000, id = 2.441)
test_temp = collect(range(85, stop = 160, length = length(testwell.md)))

pressure_values = traverse_topdown(wellbore = testwell, roughness = 0.0006, temperatureprofile = test_temp,
                                    pressurecorrelation = BeggsAndBrill, dp_est = 10, error_tolerance = 0.1,
                                    q_o = 400, q_w = 500, GLR = 2000, APIoil = 36, sg_water = 1.05, sg_gas = 0.75,
                                    WHP = 150 - 14.65)

ihs_data = joinpath(dirname(dirname(pathof(PressureDrop))), "test/testdata/Cleveland_6/Perform_results_cleveland_6_long.csv")
ihs_pressures = [parse.(Float64, split(line, ',', keepempty = false)) for line in readlines(ihs_data)[2:end]] |>
                        x -> hcat(x...)'

@test (ihs_pressures[end, 2] - 14.65) ≈ pressure_values[end] atol = 25

model = WellModel(wellbore = testwell, roughness = 0.0006, temperatureprofile = test_temp,
                    pressurecorrelation = BeggsAndBrill, dp_est = 10, error_tolerance = 0.1,
                    q_o = 400, q_w = 500, GLR = 2000, APIoil = 36, sg_water = 1.05, sg_gas = 0.75,
                    WHP = 150 - 14.65)

@test (ihs_pressures[end, 2] - 14.65) ≈ traverse_topdown(model)[end] atol = 25

#= view results
ihs_temps = "C:/pressuredrop.git/test/testdata/Cleveland_6/Perform_temps_cleveland_6_long.csv"
ihs_temps = [parse.(Float64, split(line, ',', keepempty = false)) for line in readlines(ihs_temps)[2:end]] |>
                        x -> hcat(x...)' ;

using Gadfly
set_default_plot_size(8.5inch, 11inch)

plot(   layer(x = pressure_values, y = testwell.md, Geom.line),
        layer(x = test_temp, y = testwell.md, Geom.line, Theme(default_color = "red")),
        layer(x = ihs_pressures[:,2], y = ihs_pressures[:,1], Geom.line, Theme(default_color = "green")),
        layer(x = ihs_temps[:,2], y = ihs_temps[:,1], Geom.line, Theme(default_color = "orange")),
        Coord.cartesian(yflip = true)) #matches

[pressure_values[i] - pressure_values[i-1] for i in 2:length(pressure_values)] #negative drops only occur where inclination is negative
=#

end #testset IHS Cleveland 6 - B&B


text:
function harmonic(MX, θ0, tgα0, R, N=100)
    tgα0 = -tgα0
    ω=sqrt(1.0-3*MX*G/C/C/R)

    dudθ0 = dudθ(θ0, tgα0, R)
    u0 = 1.0/R

    ϕ = atan(-dudθ0/ω/u0)-ω*θ0
    A = u0/cos(ω*θ0+ϕ)

    ϕm=π/2-ω*ϕ*2.0

    ϕl = linspace(θ0, ϕm, N)
    u = A*cos(ω*ϕl+ϕ)
    du = -A*ω*sin(ω*ϕl+ϕ)

    ixs = find(u.>=u0)
    u=u[ixs]
    du = du[ixs]
    ϕl = ϕl[ixs]

    x = [v for v in -cos(ϕl)./u]
    y = [v for v in sin(ϕl)./u]

    out_angle =  ray_angle( ϕl[end], 1/u[end], du[end])
    in_angle =  atan(tgα0)

    x, y, u, du, -in_angle, out_angle, -out_angle-in_angle, ϕl
end

harmonic200b(z, M200b, Θ0, tgα0, N) =
    harmonic(M200b, Θ0, tgα0, R200b(z, M200b), N)
harmonic200c(z, M200c, Θ0, tgα0, N) =
    harmonic(M200c, Θ0, tgα0, R200c(z, M200c), N)
harmonicVir(z, MVir, Θ0, tgα0, N) =
    harmonic(MVir, Θ0, tgα0, RVIR(z, MVir), N)

function harmonic_e(MX, rmin, tgα, R, N=100)
    tgα =- tgα
    m0 = MX*G/C/C
    ω=sqrt(1.0-3*m0/R)

    ϕ = linspace(-π/2, π/2, N)
    α = ω.*ϕ
    cα = cos(α)
    sα = sin(α)

    u = +cα/rmin + m0/R/rmin/ω/ω * (α.*sα+cα.*log(cα))

    ixs = find(u.>((1./R)))
    u=u[ixs]
    ϕ=ϕ[ixs]

    a1 = atan(tgα)
    dx = -cos(ϕ[2]+π/2)./u[2] + cos(ϕ[1]+π/2)./u[1]
    dy =  sin(ϕ[2]+π/2)./u[2] - sin(ϕ[1]+π/2)./u[1]
    a2 = -atan2(dy,dx)

    @show dα = π/2+a1-a2
    θ = ϕ+dα
    @show dα-π/2

    x = [v for v in -cos(θ)./u]
    y = [v for v in sin(θ)./u]

    # TODO: du, outgoing_angle, deflection_angle
    in_angle = atan(tgα)
    x, y, u, θ, -in_angle
end

harmonic_e200b(z, M200b, Θ0, tgα0, N) =
    harmonic_e(M200b, Θ0, tgα0, R200b(z, M200b), N)
harmonic_e200c(z, M200c, Θ0, tgα0, N) =
    harmonic_e(M200c, Θ0, tgα0, R200c(z, M200c), N)
harmonic_eVir(z, MVir, Θ0, tgα0, N) =
    harmonic_e(MVir, Θ0, tgα0, RVIR(z, MVir), N)

max_stars_repo_name: skariel/Lensing

max_issues_repo_name: skariel/Lensing

max_forks_repo_name: skariel/Lensing


text:
# #Base
# Base.parse(::T, ::Missing, kwargs...) where T = missing
#
# # Dates.jl
# Dates.Date(::Missing, kwargs...) = missing
# Dates.DateTime(::Missing, kwargs...) = missing
# Dates.year(::Missing) = missing
# Dates.month(::Missing) = missing
# Dates.day(::Missing) = missing
# Dates.canonicalize(::Missing) = missing
# Dates.yearmonthday(::Missing) = missing
# Dates.yearmonth(::Missing) = missing
# Dates.monthday(::Missing) = missing
# Dates.week(::Missing) = missing
# Dates.hour(::Missing) = missing
# Dates.minute(::Missing) = missing
# Dates.second(::Missing) = missing
# Dates.millisecond(::Missing) = missing
# Dates.dayofmonth(::Missing) = missing
# Dates.microsecond(::Missing) = missing
# Dates.nanosecond(::Missing) = missing
# Dates.dayofweek(::Missing) = missing
# Dates.isleapyear(::Missing) = missing
# Dates.daysinmonth(::Missing) = missing
# Dates.daysinyear(::Missing) = missing
# Dates.dayofyear(::Missing) = missing
# Dates.dayname(::Missing) = missing
# Dates.dayabbr(::Missing) = missing
# Dates.dayofweekofmonth(::Missing) = missing
# Dates.daysofweekinmonth(::Missing) = missing
# Dates.monthname(::Missing) = missing
# Dates.monthabbr(::Missing) = missing
# Dates.quarterofyear(::Missing) = missing
# Dates.dayofquarter(::Missing) = missing
# Dates.unix2datetime(::Missing) = missing
# Dates.datetime2unix(::Missing) = missing
# Dates.firstdayofweek(::Missing) = missing
# Dates.lastdayofweek(::Missing) = missing
# Dates.firstdayofmonth(::Missing) = missing
# Dates.lastdayofmonth(::Missing) = missing
# Dates.firstdayofyear(::Missing) = missing
# Dates.lastdayofyear(::Missing) = missing
# Dates.firstdayofquarter(::Missing) = missing
# Dates.lastdayofquarter(::Missing) = missing


text:
"""
Calculates the average of the elements in a vector ``x``.
"""
function sample_mean(x::Vector)
    n = length(x)
    x̅ = 1 / n * sum(x)
    return x̅
end

"""
Calculates the sample variance of the elements in a vector ``x``.
"""
function sample_variance(x::Vector)
    n = length(x)
    s² = 1 / (n - 1) * (sum(x.^2) - 1 / n * sum(x)^2)
    return s²
end

"""
Calculates the covariance of vector ``x``.
"""
function sample_covariance(x::Vector)
    n = length(x)
    cᵪᵪ = sum_of_squares(x) / (n - 1)
    return cᵪᵪ

end

"""
Calculates the covariance of two vectors ``x`` and ``y``.
"""
function two_sample_covariance(x::Vector, y::Vector)
    n = length(x)
    cᵪᵧ = 1 / (n -1) * (sum(x .* y) - 1 / n * sum(x) * sum(y))
    return cᵪᵧ
end

"""
Calculates the correlation coefficient of two vectors ``x`` and ``y``.
See also: [`sample_variance`](@ref), [`two_sample_covariance`](@ref).
"""
function two_sample_correlation_coefficient(x::Vector, y::Vector)
    sᵪ = √(sample_variance(x))
    sᵧ = √(sample_variance(y))
    cᵪᵧ = two_sample_covariance(x, y)
    rᵪᵧ = cᵪᵧ / (sᵪ * sᵧ)
    return rᵪᵧ
end

"""
Calculates the sample variance of two vectors ``x`` and ``y``.
Used mainly for calculating reference variables for significance tests when the variance is unknown.
See also: [`sample_variance`](@ref).
"""
function two_sample_variance(x::Vector, y::Vector)
    nᵪ = length(x)
    nᵧ = length(y)
    s²ᵪ = sample_variance(x)
    s²ᵧ = sample_variance(y)
    s²ᵪᵧ = ((nᵪ - 1) * s²ᵪ + (nᵧ - 1) * s²ᵧ) / (nᵪ + nᵧ - 2)
    return s²ᵪᵧ
end

"""
This function calculates the sum of squares.
This is used mainly for estimating the parameters in simple linear regression.
See also: [`not_implemented`](@ref).
"""
function sum_of_squares(x::Vector, y=x::Vector)
    n = length(x)
    Sᵪᵧ = (sum(x .* y) - 1 / n * sum(x) * sum(y))
    return Sᵪᵧ
end

max_issues_repo_name: persson-io/Stochastics

max_forks_repo_name: persson-io/Stochastics


text:
export Associative, Nonassociative, associativity

"""
Trait to indicate that a binary operation • is associative over set S. Only
methods of • with the signature •(x::S, y::S) are to be considered.

That is to say, ``\forall x, y, z \in S. (x \bullet y) \bullet z = x \bullet (y
\bullet z)``. Note that in practice, this is not so absolute. Minor deviations
such as `Float64` arithmetic not being perfectly associative are allowed for
practical reasons.
"""
abstract Associative

"""
Trait to indicate that a binary operation • is not necessarily associative over
set S. Despite the name, this trait does not imply the existence of a
counterexample to associativity.
"""
abstract Nonassociative

"""
Return the associativity trait of the set S under the binary operation •.
"""
associativity(::Type, ::Function) = Nonassociative

# complex numbers are associative under + and *
# note that this is not in general true for arbitrary Numbers
associativity{N<:AbstractComplex}(::Type{N}, ::typeof(+)) = Associative
associativity{N<:AbstractComplex}(::Type{N}, ::typeof(*)) = Associative

# Strings are associative under *
associativity(::Type{String}, ::typeof(*)) = Associative


text:
struct SimulationState{S1, S2, TF<:AbstractFilterState}
    hidden_state::S1
    obs::S2
    filter_state::TF
end

"""
    simulate!(filtering_algorithm, filtering_problem, no_of_timesteps, dt)

Runs a simulation of the hidden state, observation, and filtering algorithm for a duration of `no_of_timesteps`.
"""
function SimulationState(filt_prob, filt_algo)
    hidden_state = deepcopy(initialize(state_model(filt_prob)))
    obs          = obs_model(filt_prob)(hidden_state, zero(eltype(hidden_state)))
    filter_state = deepcopy(initialize(filt_algo))
    return SimulationState(hidden_state, obs, filter_state)
end

hidden_state(st::SimulationState) = st.hidden_state
obs(st::SimulationState)          = st.obs
filter_state(st::SimulationState) = st.filter_state
cond_mean(x::SimulationState)     = mean(filter_state(x))
cond_cov(x::SimulationState)      = cov(filter_state(x))
cond_var(x::SimulationState)      = var(filter_state(x))

"""
    propagate!(sfs, filtering_problem, filtering_algorithm; dt) --> sfs

Propagates the system and filter states for one time-step according to the specified filtering problem and algorithm.
"""
function propagate!(
    sfs::SimulationState,
    filt_prob::AbstractFilteringProblem{S1, S2, ContinuousTime, ContinuousTime},
    filt_algo::AbstractFilteringAlgorithm{ContinuousTime, ContinuousTime},
    dt) where {S1, S2}

    hidden_state = sfs.hidden_state
    obs          = sfs.obs
    filter_state = sfs.filter_state
    state_model  = filt_prob.state_model
    obs_model    = filt_prob.obs_model

    propagate!(hidden_state, state_model, dt)
    emit!(obs, hidden_state, obs_model, dt)
    update!(filter_state, filt_algo, obs, dt)
end

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
using DataStructures

function countcard(str::String)::Accumulator{Int, Int}
    counter([parse(Int, x) for x in str])
end

function score(cards::Accumulator{Int, Int})::Int64
    ret = 0
    for i=1:9
        ret += i*10^cards[i]
    end
    ret
end

function solve()
    k = parse(Int, readline())
    s = readline()
    t = readline()

    s = countcard(s[1:end-1])
    t = countcard(t[1:end-1])

    total = 0
    win = 0
    for a=1:9, b=1:9
        n = if a == b
            (k - (s[a] + t[a])) * (k - (s[b] + t[b] + 1))
        else
            (k - (s[a] + t[a])) * (k - (s[b] + t[b]))
        end
        total += n

        s[a] += 1
        t[b] += 1
        if score(s) > score(t)
            win += n
        end
        s[a] -= 1
        t[b] -= 1
    end
    win / total
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
### RandomPolicy ###
# maintained by @zsunberg
"""
a generic policy that uses the actions function to create a list of actions and then randomly samples an action from it.
"""
mutable struct RandomPolicy{RNG<:AbstractRNG, P<:Union{POMDP,MDP}, U<:Updater} <: Policy
    rng::RNG
    problem::P
    updater::U # set this to use a custom updater, by default it will be a void updater
end
# The constructor below should be used to create the policy so that the action space is initialized correctly
RandomPolicy(problem::Union{POMDP,MDP};
             rng=Base.GLOBAL_RNG,
             updater=VoidUpdater()) = RandomPolicy(rng, problem, updater)

## policy execution ##
function action(policy::RandomPolicy, s)
    return rand(policy.rng, actions(policy.problem, s))
end

function action(policy::RandomPolicy, b::Void)
    return rand(policy.rng, actions(policy.problem))
end

## convenience functions ##
updater(policy::RandomPolicy) = policy.updater

"""
solver that produces a random policy
"""
mutable struct RandomSolver <: Solver
    rng::AbstractRNG
end
RandomSolver(;rng=Base.GLOBAL_RNG) = RandomSolver(rng)
solve(solver::RandomSolver, problem::Union{POMDP,MDP}) = RandomPolicy(solver.rng, problem, VoidUpdater())


text:
# Ensure that we can run all files in the examples directory with no errors.

exampledir = joinpath(dirname(@__DIR__), "examples")
examples = readdir(exampledir)

mktemp() do filename, io
    redirect_stdout(io) do
        for ex in examples
            haskey(ENV, "CI") && Sys.isapple() && ex == "fstree.jl" && continue
            include(joinpath(exampledir, ex))
        end
    end
end


text:
using Test
using Shark

#=
    Depth     Nodes         Performance      Result
    -----   ---------   ------------------  --------
      0   →         1     0.000000 seconds     ✓
      1   →        20     0.001304 seconds     ✓
      2   →       400     0.027820 seconds     ✓
      3   →     8,902     0.591598 seconds     ✓
      4   →   197,281    14.038067 seconds     ✓
      5   → 4,865,609   352.093975 seconds     ✓
=#
function perft(s::State, depth::Int)::Int
    nodes = 0
    if depth == 0 return 1; end
    for move in generate_legal_moves(s)
        x = make_move(s, move)
        nodes += perft(x, depth - 1)
    end
    return nodes;
end

@testset begin
    s = State()
    @test perft(s, 0) === 1
    @test perft(s, 1) === 20
    @test perft(s, 2) === 400
    @test perft(s, 3) === 8902
    @test perft(s, 4) === 197281
    # @test perft(s, 5) === 4865609
end

max_stars_repo_name: Joshua-Robison/Shark

max_issues_repo_name: Joshua-Robison/Shark

max_forks_repo_name: Joshua-Robison/Shark