text:
##############################
# Global variables/constants #
##############################
const ADBACKEND = Ref(:forwarddiff)
setadbackend(backend_sym::Symbol) = setadbackend(Val(backend_sym))
function setadbackend(::Val{:forward_diff})
    Base.depwarn("`AdvancedVI.setadbackend(:forward_diff)` is deprecated. Please use `AdvancedVI.setadbackend(:forwarddiff)` to use `ForwardDiff`.", :setadbackend)
    setadbackend(Val(:forwarddiff))
end
function setadbackend(::Val{:forwarddiff})
    CHUNKSIZE[] == 0 && setchunksize(40)
    ADBACKEND[] = :forwarddiff
end

function setadbackend(::Val{:reverse_diff})
    Base.depwarn("`AdvancedVI.setadbackend(:reverse_diff)` is deprecated. Please use `AdvancedVI.setadbackend(:tracker)` to use `Tracker` or `AdvancedVI.setadbackend(:reversediff)` to use `ReverseDiff`. To use `ReverseDiff`, please make sure it is loaded separately with `using ReverseDiff`.",  :setadbackend)
    setadbackend(Val(:tracker))
end
function setadbackend(::Val{:tracker})
    ADBACKEND[] = :tracker
end

const ADSAFE = Ref(false)
function setadsafe(switch::Bool)
    @info("[AdvancedVI]: global ADSAFE is set as $switch")
    ADSAFE[] = switch
end

const CHUNKSIZE = Ref(40) # default chunksize used by AD

function setchunksize(chunk_size::Int)
    if ~(CHUNKSIZE[] == chunk_size)
        @info("[AdvancedVI]: AD chunk size is set as $chunk_size")
        CHUNKSIZE[] = chunk_size
    end
end

abstract type ADBackend end
struct ForwardDiffAD{chunk} <: ADBackend end
getchunksize(::Type{<:ForwardDiffAD{chunk}}) where chunk = chunk

struct TrackerAD <: ADBackend end

ADBackend() = ADBackend(ADBACKEND[])
ADBackend(T::Symbol) = ADBackend(Val(T))

ADBackend(::Val{:forwarddiff}) = ForwardDiffAD{CHUNKSIZE[]}
ADBackend(::Val{:tracker}) = TrackerAD
ADBackend(::Val) = error("The requested AD backend is not available. Make sure to load all required packages.")


text:
module TestReplaceComp

using Test
using Mimi
import Mimi:
    reset_compdefs

reset_compdefs()

@defcomp X begin
    x = Parameter(index = [time])
    y = Variable(index = [time])
    function run_timestep(p, v, d, t)
        v.y[t] = 1
    end
end

@defcomp X_repl begin
    x = Parameter(index = [time])
    y = Variable(index = [time])
    function run_timestep(p, v, d, t)
        v.y[t] = 2
    end
end

@defcomp bad1 begin
    x = Parameter()                 # parameter has same name but different dimensions
    y = Variable(index = [time])
end

@defcomp bad2 begin
    x = Parameter(index = [time])
    z = Variable(index = [time])    # different variable name
end

@defcomp bad3 begin
    z = Parameter()                 # different parameter name
    y = Variable(index = [time])
end

@defcomp bad4 begin
    x::Symbol = Parameter(index = [time])   # different datatype
    y = Variable()                          # different variable dimensions
end

# 1. Test scenario where the replacement works

m = Model()
set_dimension!(m, :time, 2000:2005)
add_comp!(m, X)                         # Original component X
set_param!(m, :X, :x, zeros(6))
replace_comp!(m, X_repl, :X)            # Successfully replaced by X_repl
run(m)
@test length(components(m)) == 1        # Only one component exists in the model
@test m[:X, :y] == 2 * ones(6)          # Successfully ran the run_timestep function from X_repl

# 2. Test bad internal incoming parameter

m = Model()
set_dimension!(m, :time, 2000:2005)
add_comp!(m, X, :first)                                     # Add two components
add_comp!(m, X, :second)
connect_param!(m, :second => :x, :first => :y)              # Make an internal connection with a parameter with a time dimension
@test_throws ErrorException replace_comp!(m, bad1, :second) # Cannot make reconnections because :x in bad1 has different dimensions
replace_comp!(m, bad1, :second, reconnect = false)          # Can replace without reconnecting
@test m.md.comp_defs[:second].comp_id.comp_name == :bad1    # Successfully replaced

# 3. Test bad internal outgoing variable

m = Model()
set_dimension!(m, :time, 2000:2005)
add_comp!(m, X, :first)                                     # Add two components
add_comp!(m, X, :second)
connect_param!(m, :second => :x, :first => :y)              # Make an internal connection from a variable with a time dimension
@test_throws ErrorException replace_comp!(m, bad2, :first)  # Cannot make reconnections because bad2 does not have a variable :y
replace_comp!(m, bad2, :first, reconnect = false)           # Can replace without reconnecting
@test m.md.comp_defs[:first].comp_id.comp_name == :bad2     # Successfully replaced

# 4. Test bad external parameter name

m = Model()
set_dimension!(m, :time, 2000:2005)
add_comp!(m, X)
set_param!(m, :X, :x, zeros(6))                     # Set external parameter for :x

# Replaces with bad3, but warns that there is no parameter by the same name :x
@test_logs (:warn, r".*parameter x no longer exists in component.*") replace_comp!(m, bad3, :X)

@test m.md.comp_defs[:X].comp_id.comp_name == :bad3 # The replacement was still successful
@test length(m.md.external_param_conns) == 0        # The external paramter connection was removed
@test length(m.md.external_params) == 1             # The external parameter still exists

# 5. Test bad external parameter dimensions

m = Model()
set_dimension!(m, :time, 2000:2005)
add_comp!(m, X)
set_param!(m, :X, :x, zeros(6))                         # Set external parameter for :x
@test_throws ErrorException replace_comp!(m, bad1, :X)  # Cannot reconnect external parameter, :x in bad1 has different dimensions

# 6. Test bad external parameter datatype

m = Model()
set_dimension!(m, :time, 2000:2005)
add_comp!(m, X)
set_param!(m, :X, :x, zeros(6))                         # Set external parameter for :x
@test_throws ErrorException replace_comp!(m, bad4, :X)  # Cannot reconnect external parameter, :x in bad4 has different datatype

# 7. Test component name that doesn't exist

m = Model()
set_dimension!(m, :time, 2000:2005)
add_comp!(m, X)
@test_throws ErrorException replace_comp!(m, X_repl, :Z)    # Component Z does not exist in the model, cannot be replaced

# 8. Test original postion placement functionality

m = Model()
set_dimension!(m, :time, 2000:2005)
add_comp!(m, X, :c1)
add_comp!(m, X, :c2)
add_comp!(m, X, :c3)
replace_comp!(m, X_repl, :c3)   # test replacing the last component
@test collect(values(m.md.comp_defs))[3].comp_id.comp_name == :X_repl
replace_comp!(m, X_repl, :c2)        # test replacing not the last one
@test collect(values(m.md.comp_defs))[2].comp_id.comp_name == :X_repl

end # module

ext: jl


text: include("reweighing.jl")


text:
#!/bin/env julia-1.6

module NodeCall

const lib_path = joinpath(dirname(@__DIR__), "build/lib")

using Libdl
const libjlnode = Libdl.dlopen(joinpath(lib_path, "libjlnode"))
const libjlnode = Libdl.dlopen(joinpath(lib_path, "libjlnode"))

const NapiValue = Ptr{Nothing}
const Env = Ref(Ptr{Nothing}())

struct JlnodeResult
    code::Int32
    message::Cstring
end
const err = Ref(JlnodeResult(0, C_NULL))

function _to_string(v)
    s = @ccall :libjlnode.value_to_string(err::Ptr{JlnodeResult}, Env[]::NapiValue, v::NapiValue)::NapiValue
    len = Ref{UInt}()
    ss = @ccall :libjlnode.string_to_utf8(
        err::Ptr{JlnodeResult}, Env[]::NapiValue,
        s::NapiValue,
        len::Ptr{Csize_t}
    )::Ptr{Cchar}
    @show len
    ss
end
_free_string(s) = @ccall :libjlnode.free_string(s::Ptr{Cchar})::Cvoid

test() = @ccall :libjlnode.test()::Cint
initialize(addon_path) = begin
    ret = @ccall :libjlnode.initialize(addon_path::Cstring)::Cint
    Env[] = @ccall :libjlnode.get_global_env()::Ptr{Cvoid}
    ret
end
dispose() = @ccall :libjlnode.dispose()::Cint

function run(scripts::AbstractString)
    result = Ref{NapiValue}()
    result = @ccall :libjlnode.run_script(err::Ptr{JlnodeResult}, Env[]::NapiValue, scripts::Cstring)::NapiValue
    if err[].code == 1 || err[].code == 2
        @error unsafe_string(err[].message)
    else
        s = _to_string(result)
        t = @ccall :libjlnode.value_type(
            err::Ptr{JlnodeResult},
            Env[]::NapiValue, result::NapiValue
        )::Cint
        if err[].code != 0
            @error unsafe_string(err[].message)
        end
        @info "Type: $t"
        @info unsafe_string(s)
        _free_string(s)
    end
    err[].code
end

end

@show NodeCall.test()
@show NodeCall.initialize(joinpath(NodeCall.lib_path, "jlnode_addon.node"))
@show NodeCall.Env[]

@show NodeCall.run("console.log(Math.random())")
@show NodeCall.run("Math.random()")
@show NodeCall.run("x")

if !Base.isinteractive()
    @show NodeCall.dispose()
end

max_stars_repo_name: sunoru/jlnode

max_issues_repo_name: sunoru/jlnode

max_forks_repo_name: sunoru/jlnode


text:
using DataFrames
using BSON
using RandomNumbers

include("reporting.jl")
include("problems.jl")


text:
#module Fragmentation
#using ..scenario, ..Utils, ..Types, ..ExternalBallistics
#using LinearAlgebra
#export frag, fragFly!
function createFragment(z::Zone, s::Array{FragShapes,1},proj::Projectile1D,α::Float64)
    number = Int64(rand(Normal(z.number)))
    rads = rand(Uniform(0,360), number)
    axs = rand(Uniform(z.angle), number)
    f = Array{Fragment}(undef,number)
    c_bar=[cos(α), -sin(α)]
    y_bar = [0.0, 1.0]
    e_bar = cross(c_bar,y_bar)/norm(cross(c_bar,y_bar))
    d_bar = cross(e_bar,c_bar)
    for i=1:number
        velocity = rand(Normal(z.velocity))
        mass = abs(rand(Normal(z.mass)))
        shape = rand(s)
        cd = rand(Normal(shape.cd))
        γ = rand(LogNormal(s.cf))
        Aₚ = γ*(mass/z.density)^(2/3)
        #vx = velocity*cosd(axs[i])*cosd(rads[i])
        #vy = velocity*sind(axs[i])
        #vz = velocity*cosd(axs[i])*sind(rads[i])
        f_velocity = proj.velocity*[1.0,0.0] + velocity*[cosd(axs[i]), sind(axs[i])] #in the fragments axis (2D)
        f_position =proj.calibre*[1/tand(axs[i]), 1] #in the fragments axis (2D)
        f_position = proj.position+proj.calibre/2*[cosd(axs[i]), cosd(rads[i]), sind(rads[i])]
        f[i] = Fragment(rads[i], axs[i], mass, z.density,cd, Aₚ, f_position,f_velocity,0.0)
    end
    return f
end

function frag(proj::Projectile2D, zones::Array{Zone,1}, shapes::Array{FragShapes,1},cf::ShapeFactor, αₑ_bar::Array{Float64,1})
    θ  = asin(-proj.velocity[3]/norm(proj.velocity))
    Ψ = atan(proj.velocity[2]/proj.velocity[1])
    T = [cos(θ)*cos(Ψ) -sin(Ψ) sin(θ)*cos(Ψ); cos(θ)*sin(Ψ) cos(Ψ) sin(θ)*sin(Ψ); -sin(θ) 0.0 cos(θ)]
    frags = []
    for z in zones
        f = createFragment(z,shapes,cf, proj, αₑ_bar)

    #=    for i in f

            v = i.velocity#[i.vx, i.vy, i.vz]
            vrot = rotation(v, T)
            i.velocity = vrot
            #i.vy = vrot[2]
            #i.vz = vrot[3]
        end
        =#
        push!(frags,f)
    end
    return frags
end

function fragFly!(fragments, proj::Projectile2D, w_bar::Array{Float64,1}, target::Target2D, g₀::Float64, ω_bar::Array{Float64,1}, coeffAero;R=6.356766*1e6, tspan=(0.0,1000.0))
    for i=1:length(fragments)
        for j=1:length(fragments[i])

            u0=[proj.position[1], proj.position[2], proj.position[3],fragments[i][j].velocity[1], fragments[i][j].velocity[2], fragments[i][j].velocity[3]]

            p= [w_bar, fragments[i][j].Aₚ, R, g₀,ω_bar, fragments[i][j].mass,target.position,coeffAero,fragments[i][j]]
            ExternalBallistics.Pmm.trajectoryPMM!(u0, tspan, p, fragments[i][j])
        end

    end

end

#end


text:
module FilteredGroupbyMacro

import DataFrames

export @by

macro by(exp)
    # dump(exp)

    exp.head != :ref && error("You have to use [] syntax.")

    df = exp.args[1]

    # enable chained calls by applying the macro again if df is actually another nested
    # ref expression
    if df isa Expr
        df = :(@by($df))
    end

    # nargs = length(exp.args) - 1
    args = exp.args[2:end]

    # additional keyword parameters after the ;
    kwparams = if args[1] isa Expr && args[1].head == :parameters
        params = args[1].args
        deleteat!(args, 1)
        esc(params...)
    else
        # empty splatted keyword params
        esc(:(NamedTuple()...))
    end

    filterexp = args[1]
    byexp = args[2]

    action_exps = args[3:end]

    validate_byexp(byexp)

    conv_filterexp = replace_quotenodes_filterexp!(filterexp, df)

    assign_mode = is_assign_kwarg(action_exps[1])

    # for x := sum(:y) style
    if assign_mode

        conv_assign_exps = convert_assignexp.(action_exps)

        quote

            tempdf = $(esc(DataFrames.by))(
                # filter
                $(esc(df))[$(esc(conv_filterexp)), $(esc(:))],
                # groupby
                $byexp;
                # new columns from computations
                $(conv_assign_exps...),
                $kwparams
            )

            $(esc(DataFrames.join))($(esc(df)), tempdf, on = $byexp)

        end

    # for x = sum(:y) style
    else

        conv_kwexps = convert_kwexp.(action_exps)

        quote
            $(esc(DataFrames.by))(
                # filter
                $(esc(df))[$(esc(conv_filterexp)), $(esc(:))],
                # groupby
                $byexp;
                # new columns from computations
                $(conv_kwexps...),
                $kwparams
            )
        end

    end
end

is_assign_kwarg(x) = false
is_assign_kwarg(exp::Expr) = exp.head == Symbol(":=")

is_kwarg(x) = false
is_kwarg(exp::Expr) = exp.head == :kw

function validate_byexp(exp)
    is_kwarg(exp) && error("Groupby expression missing, found keyword argument $exp instead.")
end

function convert_kwexp(kwexp)
    !is_kwarg(kwexp) && error("Expected a keyword argument but received $kwexp")
    kw = kwexp.args[1]
    exp = kwexp.args[2]

    replace_quotenodes!(exp)
    quotenodes = all_quotenodes(exp)

    :($(esc(kw)) = ($(quotenodes...),) => $(esc(:subdf)) -> $(esc(exp)))
end

function convert_assignexp(assignexp)
    !is_assign_kwarg(assignexp) && error("Expected an assignment keyword argument but received $assignexp")

    kw = assignexp.args[1]
    exp = assignexp.args[2]

    replace_quotenodes!(exp)
    quotenodes = all_quotenodes(exp)

    :($(esc(kw)) = ($(quotenodes...),) => $(esc(:subdf)) -> $(esc(exp)))
end

all_quotenodes(any) = QuoteNode[]
all_quotenodes(qn::QuoteNode) = [qn]

function all_quotenodes(exp::Expr)
    quotenodes = QuoteNode[]
    for a in exp.args
        append!(quotenodes, all_quotenodes(a))
    end
    quotenodes
end

replace_quotenodes!(any) = any

function replace_quotenodes!(exp::Expr)
    for (i, a) in enumerate(exp.args)
        if a isa QuoteNode
            sym = a.value
            exp.args[i] = :(subdf.$sym)
        elseif a isa Expr
            replace_quotenodes!(a)
        end
    end
    exp
end

replace_quotenodes_filterexp!(any, dfexp) = any

function replace_quotenodes_filterexp!(exp::Expr, dfexp)
    for (i, a) in enumerate(exp.args)
        if a isa QuoteNode
            sym = a.value
            exp.args[i] = :($dfexp.$sym)
        elseif a isa Expr
            replace_quotenodes_filterexp!(a, dfexp)
        end
    end
    exp
end

end # module

max_stars_repo_name: jkrumbiegel/GroupbyIndexingMacro


text:
"""
    get_user(::Forge[, name_or_id::Union{$AStr, Integer}])

Get the currently authenticated user, or a user by name or ID.
"""
@endpoint get_user()
@endpoint get_user(name_or_id::Union{AStr, Integer})

"""
    get_users(::Forge)

Get all users.
"""
@endpoint get_users()

"""
    update_user(::Forge[, id::Integer]; kwargs...)

Update the currently authenticated user, or a user by ID.
"""
@endpoint update_user()
@endpoint update_user(id::Integer)

"""
    create_user(::Forge; kwargs...)

Create a new user.
"""
@endpoint create_user()

"""
    delete_user(::Forge, id::Integer)

Delete a user by ID.
"""
@endpoint delete_user(id::Integer)

"""
    get_user_repos(::Forge[, name_or_id::Union{$AStr, Integer}])

Get the currently authenticated user's repositories, or those of a user by name or ID.
"""
@endpoint get_user_repos()
@endpoint get_user_repos(name_or_id::Union{AStr, Integer})

"""
    get_repo(::Forge, owner::$AStr, repo::$AStr)
    get_repo(::Forge, id::Integer)

Get a repository by owner and name or ID.
"""
@endpoint get_repo(owner::AStr, repo::AStr)
@endpoint get_repo(id::Integer)

"""
    get_branch(::Forge, owner::AbstractString, repo::AbstractString, branch::AbstractString)

Get a branch from a repository.
"""
@endpoint get_branch(owner::AStr, repo::AStr, branch::AStr)

"""
    get_file_contents(
        ::Forge,
        owner::$AStr,
        repo::$AStr,
        path::$AStr,
    )
    get_file_contents(f::Forge, id::Integer, path::$AStr)

Get a file from a repository.
"""
@endpoint get_file_contents(owner::AStr, repo::AStr, path::AStr)
@endpoint get_file_contents(id::Integer, path::AStr)

"""
    get_pull_request(::Forge, owner::$AStr, repo::$AStr, number::Integer)
    get_pull_request(::Forge, project::Integer, number::Integer)

Get a specific pull request.
"""
@endpoint get_pull_request(owner::AStr, repo::AStr, number::Integer)
@endpoint get_pull_request(project::Integer, number::Integer)

"""
    get_pull_requests(::Forge, owner::$AStr, repo::$AStr)
    get_pull_requests(::Forge, project::Integer)

List a repository's pull requests.
"""
@endpoint get_pull_requests(owner::AStr, repo::AStr)
@endpoint get_pull_requests(repo::Integer)

"""
    create_pull_requests(::Forge, owner::$AStr, repo::$AStr; kwargs...)
    create_pull_requests(::Forge, project::Integer; kwargs...)

Create a pull request.
"""
@endpoint create_pull_request(owner::AStr, repo::AStr)
@endpoint create_pull_request(project::Integer)

"""
    update_pull_request(::Forge, owner::$AStr, repo::$AStr, number::Integer; kwargs...)
    update_pull_request(::Forge, project::Integer, number::Integer; kwargs...)

Update a pull request.
"""
@endpoint update_pull_request(owner::AStr, repo::AStr, number::Integer)
@endpoint update_pull_request(project::Integer, number::Integer)

"""
    get_commit(::Forge, owner::$AStr, repo::$AStr, ref::$AStr)
    get_commit(::Forge, project::Integer, ref::$AStr)

Get a commit from a repository.
"""
@endpoint get_commit(owner::AStr, repo::AStr, ref::AStr)
@endpoint get_commit(project::Integer, ref::AStr)

"""
    is_collaborator(
        ::Forge,
        owner::$AStr,
        repo::$AStr,
        name_or_id::Union{$AStr, Integer},
    )

Check whether or not a user is a collaborator on a repository.
"""
@endpoint is_collaborator(owner::AStr, repo::AStr, name_or_id::Union{AStr, Integer})

"""
    is_member(::Forge, org::$AStr, name_or_id::Union{$AStr, Integer})

Check whether or not a user is a member of an organization.
"""
@endpoint is_member(org::AStr, name_or_id::Union{AStr, Integer})

"""
    get_tags(::Forge, owner::$AStr, repo::$AStr)
    get_tags(::Forge, project::Integer)

Get a list of tags from a repository.
"""
@endpoint get_tags(owner::AStr, repo::AStr)
@endpoint get_tags(project::Integer)


text:
include("solve_method.jl")
include("centrality_corr.jl")
include("regularization.jl")
include("system_write.jl")
include("preconditioners/abstract-precond.jl")
include("solvers/augmented/augmented.jl")
include("solvers/Newton/Newton.jl")
include("solvers/normal/normal.jl")
include("solvers/Krylov_utils.jl")

function compute_α_dual(v, dir_v)
  n = length(v)
  T = eltype(v)
  if n == 0
    return one(T)
  end
  α = one(T)
  @inbounds for i = 1:n
    if dir_v[i] < zero(T)
      α_new = -v[i] * T(0.999) / dir_v[i]
      if α_new < α
        α = α_new
      end
    end
  end
  return α
end

function compute_α_primal(v, dir_v, lvar, uvar)
  n = length(v)
  T = eltype(v)
  α_l, α_u = one(T), one(T)
  @inbounds for i = 1:n
    if dir_v[i] > zero(T)
      α_u_new = (uvar[i] - v[i]) * T(0.999) / dir_v[i]
      if α_u_new < α_u
        α_u = α_u_new
      end
    elseif dir_v[i] < zero(T)
      α_l_new = (lvar[i] - v[i]) * T(0.999) / dir_v[i]
      if α_l_new < α_l
        α_l = α_l_new
      end
    end
  end
  return min(α_l, α_u)
end

@inline function compute_αs(x, s_l, s_u, lvar, uvar, Δxy, Δs_l, Δs_u, nvar)
  α_pri = @views compute_α_primal(x, Δxy[1:nvar], lvar, uvar)
  α_dual_l = compute_α_dual(s_l, Δs_l)
  α_dual_u = compute_α_dual(s_u, Δs_u)
  return α_pri, min(α_dual_l, α_dual_u)
end

@inline function compute_μ(x_m_lvar, uvar_m_x, s_l, s_u, nb_low, nb_upp)
  return (dot(s_l, x_m_lvar) + dot(s_u, uvar_m_x)) / (nb_low + nb_upp)
end

function update_pt!(x, y, s_l, s_u, α_pri, α_dual, Δxy, Δs_l, Δs_u, ncon, nvar)
  x .= @views x .+ α_pri .* Δxy[1:nvar]
  y .= @views y .+ α_dual .* Δxy[(nvar + 1):(ncon + nvar)]
  s_l .= s_l .+ α_dual .* Δs_l
  s_u .= s_u .+ α_dual .* Δs_u
end

function safe_boundary(v::T) where {T <: Real}
  if v == 0
    v = eps(T)^2
  end
  return v
end

# "security" if x is too close from lvar or uvar
function boundary_safety!(x_m_lvar, uvar_m_x)
  x_m_lvar .= safe_boundary.(x_m_lvar)
  uvar_m_x .= safe_boundary.(uvar_m_x)
end

function update_IterData!(itd, pt, fd, id, safety)
  T = eltype(itd.x_m_lvar)
  itd.x_m_lvar .= @views pt.x[id.ilow] .- fd.lvar[id.ilow]
  itd.uvar_m_x .= @views fd.uvar[id.iupp] .- pt.x[id.iupp]
  safety && boundary_safety!(itd.x_m_lvar, itd.uvar_m_x)

  itd.μ = compute_μ(itd.x_m_lvar, itd.uvar_m_x, pt.s_l, pt.s_u, id.nlow, id.nupp)
  mul!(itd.Qx, fd.Q, pt.x)
  itd.xTQx_2 = dot(pt.x, itd.Qx) / 2
  fd.uplo == :U ? mul!(itd.ATy, fd.A, pt.y) : mul!(itd.ATy, fd.A', pt.y)
  fd.uplo == :U ? mul!(itd.Ax, fd.A', pt.x) : mul!(itd.Ax, fd.A, pt.x)
  itd.cTx = dot(fd.c, pt.x)
  itd.pri_obj = itd.xTQx_2 + itd.cTx + fd.c0
  if typeof(pt.x) <: Vector
    itd.dual_obj = @views dot(fd.b, pt.y) - itd.xTQx_2 + dot(pt.s_l, fd.lvar[id.ilow]) -
           dot(pt.s_u, fd.uvar[id.iupp]) + fd.c0
  else # views and dot not working with GPU arrays
    itd.dual_obj = dual_obj_gpu(
      fd.b,
      pt.y,
      itd.xTQx_2,
      pt.s_l,
      pt.s_u,
      fd.lvar,
      fd.uvar,
      fd.c0,
      id.ilow,
      id.iupp,
      itd.store_vdual_l,
      itd.store_vdual_u,
    )
  end
  itd.pdd = abs(itd.pri_obj - itd.dual_obj) / (one(T) + abs(itd.pri_obj))
end

function update_data!(
  pt::Point{T},
  α_pri::T,
  α_dual::T,
  itd::IterData{T},
  pad::PreallocatedData{T},
  res::AbstractResiduals{T},
  fd::QM_FloatData{T},
  id::QM_IntData,
) where {T <: Real}

  # (x, y, s_l, s_u) += α * Δ
  update_pt!(
    pt.x,
    pt.y,
    pt.s_l,
    pt.s_u,
    α_pri,
    α_dual,
    itd.Δxy,
    itd.Δs_l,
    itd.Δs_u,
    id.ncon,
    id.nvar,
  )
  update_IterData!(itd, pt, fd, id, true)

  #update Residuals
  res.rb .= itd.Ax .- fd.b
  res.rc .= itd.ATy .- itd.Qx .- fd.c
  res.rc[id.ilow] .+= pt.s_l
  res.rc[id.iupp] .-= pt.s_u
  # update stopping criterion values:
  #         rcNorm, rbNorm = norm(rc), norm(rb)
  #         xNorm = norm(x)
  #         yNorm = norm(y)
  #         optimal = pdd < ϵ_pdd && rbNorm < ϵ_rb * max(1, bNorm + ANorm * xNorm) &&
  #                     rcNorm < ϵ_rc * max(1, cNorm + QNorm * xNorm + ANorm * yNorm)
  res.rcNorm, res.rbNorm = norm(res.rc, Inf), norm(res.rb, Inf)
  typeof(res) <: ResidualsHistory && push_history_residuals!(res, itd, pad, id)
end

function iter!(
  pt::Point{T},
  itd::IterData{T},
  fd::Abstract_QM_FloatData{T},
  id::QM_IntData,
  res::AbstractResiduals{T},
  sc::StopCrit{Tc},
  dda::DescentDirectionAllocs{T},
  pad::PreallocatedData{T},
  ϵ::Tolerances{T},
  cnts::Counters,
  T0::DataType,
  display::Bool,
) where {T <: Real, Tc <: Real}
  @inbounds while cnts.k < sc.max_iter && !sc.optimal && !sc.tired
    time_fact = (cnts.kc == -1) ? time_ns() : UInt(0) # timer centrality_corr factorization
    out = update_pad!(pad, dda, pt, itd, fd, id, res, cnts, T0) # update data for the solver! function used
    time_fact = (cnts.kc == -1) ? time_ns() - time_fact : 0.0
    out == 1 && break

    time_solve = (cnts.kc == -1) ? time_ns() : 0.0 # timer centrality_corr solve

    # Solve system to find a direction of descent
    out = update_dd!(dda, pt, itd, fd, id, res, pad, cnts, T0)
    out == 1 && break

    if typeof(pt.x) <: Vector
      α_pri, α_dual =
        compute_αs(pt.x, pt.s_l, pt.s_u, fd.lvar, fd.uvar, itd.Δxy, itd.Δs_l, itd.Δs_u, id.nvar)
    else
      α_pri, α_dual = compute_αs_gpu(
        pt.x,
        pt.s_l,
        pt.s_u,
        fd.lvar,
        fd.uvar,
        itd.Δxy,
        itd.Δs_l,
        itd.Δs_u,
        id.nvar,
        itd.store_vpri,
        itd.store_vdual_l,
        itd.store_vdual_u,
      )
    end

    if cnts.kc > 0   # centrality corrections
      α_pri, α_dual =
        multi_centrality_corr!(dda, pad, pt, α_pri, α_dual, itd, fd, id, cnts, res, T0)
      ## TODO replace by centrality_corr.jl, deal with α
    end

    update_data!(pt, α_pri, α_dual, itd, pad, res, fd, id) # update point, residuals, objectives...
    time_solve = (cnts.kc == -1) ? time_ns() - time_solve : UInt(0)
    (cnts.kc == -1) && nb_corrector_steps!(cnts, time_fact, time_solve)

    sc.optimal = itd.pdd < ϵ.pdd && res.rbNorm < ϵ.tol_rb && res.rcNorm < ϵ.tol_rc
    sc.small_μ = itd.μ < ϵ.μ

    cnts.k += 1
    if T == Float32
      cnts.km += 1
    elseif T == Float64
      cnts.km += 4
    else
      cnts.km += 16
    end

    sc.Δt = time() - sc.start_time
    sc.tired = sc.Δt > sc.max_time

    if display == true
      @info log_row(
        Any[
          cnts.k,
          itd.minimize ? itd.pri_obj : -itd.pri_obj,
          itd.pdd,
          res.rbNorm,
          res.rcNorm,
          α_pri,
          α_dual,
          itd.μ,
          get_nprod!(pad),
        ],
      )
    end
  end
end


text:
include("acqr.jl")

include("acqr_simula_test.jl")
include("acqr_kfilter_test.jl")
include("acqr_ksmoother_test.jl")
include("acqr_em_starting_test.jl")
include("acqr_em_test.jl")
include("acqr_kfilter_s_test.jl")
include("acqr_ksmoother_s_test.jl")
include("acqr_em_s_test.jl")
include("acqr_kalman_test.jl")


text:
# Clique types

# this is a developmental type, will be standardized after conclusion of #1010
# TODO resolve type instability
const MsgRelativeType = Vector{NamedTuple{(:variables, :likelihood), Tuple{Vector{Symbol},DFG.AbstractRelative}}}

const MsgPriorType = Dict{Symbol, MsgPrior{<:ManifoldKernelDensity}}

"""
    $TYPEDEF

Internal development types used during consolidation.  Stores relative and prior information making up a joint likelihood
message passed upward on the Bayes tree.
"""
mutable struct _MsgJointLikelihood
  relatives::IIF.MsgRelativeType
  priors::IIF.MsgPriorType
end

"""
  $(TYPEDEF)
Belief message for message passing on the tree.  This should be considered an incomplete joint probility.

Notes:
- belief -> Dictionary of [`TreeBelief`](@ref)
- variableOrder -> Ordered variable id list of the seperators in cliqueLikelihood
- cliqueLikelihood -> marginal distribution (<: `SamplableBelief`) over clique seperators.
- Older names include: productFactor, Fnew, MsgPrior, LikelihoodMessage

DevNotes:
- Used by both nonparametric and parametric.
- Objective for parametric case: `MvNormal(μ=[:x0;:x2;:l5], Σ=[+ * *; * + *; * * +])`.
- Part of the consolidation effort, see #459.
- Better conditioning for joint structure in the works using deconvolution, see #579, #635.
  - TODO confirm why <: Singleton.

$(TYPEDFIELDS)
"""
mutable struct LikelihoodMessage{T <: MessageType} <: AbstractPrior
  sender::NamedTuple{(:id,:step),Tuple{Int,Int}}
  status::CliqStatus
  belief::Dict{Symbol, TreeBelief} # will eventually be deprecated
  variableOrder::Vector{Symbol}
  cliqueLikelihood::Union{Nothing,SamplableBelief}  # TODO drop the Union
  msgType::T
  hasPriors::Bool
  # this is different from belief[].inferdim, as the total available infer dims remaining during down msgs -- see #910
  childSolvDims::Dict{Int, Float64}
  # calc differential factors for joint in the child clique
  jointmsg::_MsgJointLikelihood
  # diffJoints::Vector{NamedTuple{(:variables, :likelihood), Tuple{Vector{Symbol},DFG.AbstractRelative}}}
end

"""
    $TYPEDEF

Cache messages being passed on the tree, one container per clique.

Notes
- See model 2 (?) on IIF #674
"""
mutable struct MessageBuffer
  # up receive message buffer (multiple children, multiple messages)
  upRx::Dict{Int, LikelihoodMessage}
  # down receive message buffer (one parent)
  downRx::Union{Nothing, LikelihoodMessage}
  # RESERVED up outgoing message buffer (one parent)
  upTx::Union{Nothing, LikelihoodMessage}
  # RESERVED down outgoing message buffer (multiple children but one message)
  downTx::Union{Nothing, LikelihoodMessage}
end
MessageBuffer() = MessageBuffer(Dict{Int, LikelihoodMessage}(), nothing, nothing, nothing)

##==============================================================================
## BayesTreeNodeData
##==============================================================================

"""
$(TYPEDEF)

Data structure for each clique in the Bayes (Junction) tree.
"""
mutable struct BayesTreeNodeData
  status::CliqStatus
  frontalIDs::Vector{Symbol}
  separatorIDs::Vector{Symbol}
  inmsgIDs::Vector{Symbol} # Int
  potIDs::Vector{Symbol} # Int # this is likely redundant TODO -- remove
  potentials::Vector{Symbol}
  partialpotential::Vector{Bool}

  dwnPotentials::Vector{Symbol}
  dwnPartialPotential::Vector{Bool}

  cliqAssocMat::Array{Bool,2}
  cliqMsgMat::Array{Bool,2}
  directvarIDs::Vector{Symbol}
  directFrtlMsgIDs::Vector{Symbol}
  msgskipIDs::Vector{Symbol}
  itervarIDs::Vector{Symbol}
  directPriorMsgIDs::Vector{Symbol}
  debug
  debugDwn

  allmarginalized::Bool
  initialized::Symbol
  upsolved::Bool
  downsolved::Bool
  isCliqReused::Bool             # holdover

  # JT Local messages saved for cache and debugging, see IIF #675
  messages::MessageBuffer
end

## Packed types for serialization

mutable struct PackedBayesTreeNodeData
  frontalIDs::Vector{Symbol}
  separatorIDs::Vector{Symbol}
  inmsgIDs::Vector{Symbol} # Int
  potIDs::Vector{Symbol} # Int # this is likely redundant TODO -- remove
  potentials::Vector{Symbol}
  partialpotential::Vector{Bool}
  dwnPotentials::Vector{Symbol}
  dwnPartialPotential::Vector{Bool}
  cliqAssocMat::Array{Bool,2}
  cliqMsgMat::Array{Bool,2}
  directvarIDs::Vector{Symbol} # Int
  directFrtlMsgIDs::Vector{Symbol} # Int
  msgskipIDs::Vector{Symbol} # Int
  itervarIDs::Vector{Symbol} # Int
  directPriorMsgIDs::Vector{Symbol} # Int
end

## Full Clique Types

struct CliqueId{T}
  value::T
end

"""
    $(TYPEDEF)
Structure to store clique data
DEV NOTES: To replace TreeClique completely
    $(FIELDS)
"""
mutable struct TreeClique
  "Interger id unique within a tree with userId, robotId, sessionId"
  id::CliqueId{Int64} # not to be confused with the underlying index used by LightGraphs.jl, see issue #540
  "Data as `BayesTreeNodeData`"
  data::BayesTreeNodeData
  "Drawing attributes"
  attributes::Dict{String, Any}
  #solveInProgress #on a clique level a "solve in progress" might be very handy
end

#


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiDiscoveryV1beta1Endpoint <: SwaggerModel
    addresses::Any # spec type: Union{ Nothing, Vector{String} } # spec name: addresses
    conditions::Any # spec type: Union{ Nothing, IoK8sApiDiscoveryV1beta1EndpointConditions } # spec name: conditions
    hostname::Any # spec type: Union{ Nothing, String } # spec name: hostname
    targetRef::Any # spec type: Union{ Nothing, IoK8sApiCoreV1ObjectReference } # spec name: targetRef
    topology::Any # spec type: Union{ Nothing, Dict{String, String} } # spec name: topology

    function IoK8sApiDiscoveryV1beta1Endpoint(;addresses=nothing, conditions=nothing, hostname=nothing, targetRef=nothing, topology=nothing)
        o = new()
        validate_property(IoK8sApiDiscoveryV1beta1Endpoint, Symbol("addresses"), addresses)
        setfield!(o, Symbol("addresses"), addresses)
        validate_property(IoK8sApiDiscoveryV1beta1Endpoint, Symbol("conditions"), conditions)
        setfield!(o, Symbol("conditions"), conditions)
        validate_property(IoK8sApiDiscoveryV1beta1Endpoint, Symbol("hostname"), hostname)
        setfield!(o, Symbol("hostname"), hostname)
        validate_property(IoK8sApiDiscoveryV1beta1Endpoint, Symbol("targetRef"), targetRef)
        setfield!(o, Symbol("targetRef"), targetRef)
        validate_property(IoK8sApiDiscoveryV1beta1Endpoint, Symbol("topology"), topology)
        setfield!(o, Symbol("topology"), topology)
        o
    end
end # type IoK8sApiDiscoveryV1beta1Endpoint

const _property_map_IoK8sApiDiscoveryV1beta1Endpoint = Dict{Symbol,Symbol}(Symbol("addresses")=>Symbol("addresses"), Symbol("conditions")=>Symbol("conditions"), Symbol("hostname")=>Symbol("hostname"), Symbol("targetRef")=>Symbol("targetRef"), Symbol("topology")=>Symbol("topology"))
const _property_types_IoK8sApiDiscoveryV1beta1Endpoint = Dict{Symbol,String}(Symbol("addresses")=>"Vector{String}", Symbol("conditions")=>"IoK8sApiDiscoveryV1beta1EndpointConditions", Symbol("hostname")=>"String", Symbol("targetRef")=>"IoK8sApiCoreV1ObjectReference", Symbol("topology")=>"Dict{String, String}")
Base.propertynames(::Type{ IoK8sApiDiscoveryV1beta1Endpoint }) = collect(keys(_property_map_IoK8sApiDiscoveryV1beta1Endpoint))
Swagger.property_type(::Type{ IoK8sApiDiscoveryV1beta1Endpoint }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiDiscoveryV1beta1Endpoint[name]))}
Swagger.field_name(::Type{ IoK8sApiDiscoveryV1beta1Endpoint }, property_name::Symbol) =  _property_map_IoK8sApiDiscoveryV1beta1Endpoint[property_name]

function check_required(o::IoK8sApiDiscoveryV1beta1Endpoint)
    (getproperty(o, Symbol("addresses")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiDiscoveryV1beta1Endpoint }, name::Symbol, val)
end


text:
using NCDatasets

config = Dict{Any, Any}(

    "DRIVER" => Dict(
        "casename"           => "Sandbox",
        "caseroot"           => joinpath(@__DIR__, "Sandbox", "caseroot"),
        "caserun"            => joinpath(@__DIR__, "Sandbox", "caserun"),
        "archive_root"       => joinpath(@__DIR__, "Sandbox", "archive"),
    ),

    "MODEL_MISC" => Dict(
        "timetype"               => "DateTimeNoLeap",
        "init_file"              => "/billevans/projects/IOM/test/test02_Qflx_finding_driver/init_ocn.jld2",
        "rpointer_file"          => "rpointer.hoom",
        "daily_record"           => [],
        "monthly_record"         => ["ESSENTIAL",],
        "enable_archive"         => true,
    ),

    "MODEL_CORE" => Dict(

        "domain_file"                  => joinpath(@__DIR__, "CESM_domains", "domain.ocn.gx3v7.120323.nc"),
        "topo_file"                    => joinpath(@__DIR__, "Nz_bot.nc"),
        "cdata_file"                   => joinpath(@__DIR__, "POP2PROFILE.g37.nc"),

#        "domain_file"                  => joinpath(@__DIR__, "CESM_domains", "domain.ocn.gx1v6.090206.nc"),
#        "topo_file"                    => joinpath(@__DIR__, "Nz_bot.nc"),
#        "cdata_file"                   => joinpath(@__DIR__, "POP2PROFILE.g16.nc"),

        "cdata_beg_time"               => "0001-01-01 00:00:00",
        "cdata_end_time"               => "0002-01-01 00:00:00",
        "cdata_align_time"             => "0001-01-01 00:00:00",

        "z_w"                          => nothing,

        "substeps"                     => 8,
        "MLD_scheme"                   => "datastream",
        "Qflx"                         => "off",
        "Qflx_finding"                 => "off",
        "convective_adjustment"        => "on",
        #"advection_scheme"             => "ekman_codron2012_partition",
        "advection_scheme"             => "ekman_AGA2020",

        "weak_restoring"               => "off",
        "τwk_TEMP"                     => 86400.0 * 365,
        "τwk_SALT"                     => 86400.0 * 365,

        "τ_frz"                        => 3600.0,
        "Ekman_layers"                 => 5,
        "Returnflow_layers"            => 25,

        "transform_vector_field"       => false,
    ),

)

Dataset(config["MODEL_CORE"]["cdata_file"], "r") do ds

    z_w_top = nomissing(ds["z_w_top"][:], NaN)
    z_w_bot = nomissing(ds["z_w_bot"][:], NaN)

    z_w      = zeros(Float64, length(z_w_top)+1)
    z_w[1:end-1]   = z_w_top
    z_w[end] = z_w_bot[end]

    global config["MODEL_CORE"]["z_w"]  = z_w
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
function parseRPNstring(rpns)
    infix = []
    rpn = split(rpns)
    for tok in rpn
        if all(isnumber, tok)
            push!(infix, parse(Int, tok))
        else
            last = pop!(infix)
            prev = pop!(infix)
            push!(infix, Expr(:call, Symbol(tok), prev, last))
            println("Current step: $infix")
        end
    end
    infix
end

unany(s) = replace(string(s), r"Any\[:\((.+)\)\]", s"\1")

println("The final infix result: ", parseRPNstring("3 4 2 * 1 5 - 2 3 ^ ^ / +") |> unany, "\n")
println("The final infix result: ", parseRPNstring("1 2 + 3 4 + ^ 5 6 + ^") |> unany)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module Model

using Distributions

# TODO: Remove this dependency
using RCall # Mclust

import LinearAlgebra # Identity matrix
import Random # shuffle
import StatsBase # wsample, counts

include("../MCMC/MCMC.jl")
import .MCMC.Util.@namedargs

include("util.jl")
include("DICparam.jl")
include("State.jl")
include("Data.jl")
include("Constants.jl")
include("Tuners.jl")
include("repFAM/repFAM.jl")
include("update.jl")
include("FeatureSelect/FeatureSelect.jl")
include("genInitialState.jl")
include("fit.jl")

#=
precompile(cytof5_fit, (State, Constants, Data, Int, Int,
                        Vector{Vector{Symbol}}, Vector{Int}, Bool, Int, Bool,
                        Bool, Bool))
=#

end # Model


text:
@testset "Solver: second-order cone program (friction)" begin
	"""
		minimize   v' b
		subject to ||b|| <= μn
	"""

    # second-order cone
    a = [2.0; 1.0]
    b = [1.0; 0.1]

    c = RoboDojo.cone_product(a, b)
    @test norm(c - [2.1; 1.2]) < 1.0e-8

    d = zeros(2)
    RoboDojo.cone_product!(d, a, b)
    @test norm(d - c) < 1.0e-8

    # number of decision variables
    n = 7
    m = 3
    idx = RoboDojo.IndicesOptimization(
        n, n,
        [Int[], Int[]],
        [Int[], Int[]],
        [[collect(1:3), collect(5:7)]],
        [[collect(1:3), collect(5:7)]],
        collect(1:4),
        collect(1:0),
        collect(5:7),
        [collect(5:7)],
        collect(5:7),
    )

    # residual
    function _r!(r, z, θ, κ)
        x = z[1:3]
        y = z[4:4]
        z = z[5:7]

        v = θ[1:2]
        μγ = θ[3]

        r[1:3] = [0.0; v] + [y; 0.0; 0.0] - z
        r[4] = x[1] - μγ
        r[5:7] = RoboDojo.cone_product(x, z) - [κ[1]; 0.0; 0.0]
        nothing
    end

    @variables r_sym[1:n]
    @variables z_sym[1:n]
    @variables θ_sym[1:m]
    @variables κ_sym[1:1]

    parallel = Symbolics.SerialForm()
    _r!(r_sym, z_sym, θ_sym, κ_sym)
    r_sym = simplify.(r_sym)
    rf! = eval(Symbolics.build_function(r_sym, z_sym, θ_sym, κ_sym,
        parallel = parallel)[2])
    rz_exp = Symbolics.jacobian(r_sym, z_sym, simplify = true)
    rθ_exp = Symbolics.jacobian(r_sym, θ_sym, simplify = true)
    rz_sp = similar(rz_exp, Float64)
    rθ_sp = similar(rθ_exp, Float64)
    rzf! = eval(Symbolics.build_function(rz_exp, z_sym, θ_sym,
        parallel = parallel)[2])
    rθf! = eval(Symbolics.build_function(rθ_exp, z_sym, θ_sym,
        parallel = parallel)[2])

    # problem setup
    v = [0.0; 0.0]
    μγ = 1.0
    θ = [v; μγ]

    z = ones(n)
    z[1] += 1.0
    z[5] += 1.0

    # solver
    opts = RoboDojo.InteriorPointOptions(diff_sol = true)
    ip = RoboDojo.interior_point(z, θ,
        idx=idx,
        r! = rf!, rz! = rzf!, rθ! = rθf!,
        rz = rz_sp,
        rθ = rθ_sp,
        opts = opts)

    # solve
    status = RoboDojo.interior_point_solve!(ip)

    # test
    @test status
    @test norm(ip.r, Inf) < opts.r_tol
    @test all([norm(ip.z[j][2:3]) <= ip.z[j][1] for i in ip.idx.socz for j in i])
    @test ip.κ[1] < opts.κ_tol
    @test norm(ip.δz, 1) != 0.0

    # problem setup
    v = [1.0; 0.0]
    μγ = 1.0
    θ = [v; μγ]

    z = ones(n)
    z[1] += 1.0
    z[5] += 1.0

    ip.z .= z
    ip.θ .= θ

    # solve
    status = RoboDojo.interior_point_solve!(ip)

    # test
    @test status
    @test norm(ip.r, Inf) < opts.r_tol
    @test all([norm(ip.z[j][2:3]) <= ip.z[j][1] for i in ip.idx.socz for j in i])
    @test ip.κ[1] < opts.κ_tol
    @test norm(ip.δz, 1) != 0.0

    # problem setup
    v = [1.0; 10.0]
    μγ = 1.0
    θ = [v; μγ]

    z = ones(n)
    z[1] += 1.0
    z[5] += 1.0

    ip.z .= z
    ip.θ .= θ

    # solve
    status = RoboDojo.interior_point_solve!(ip)

    # test
    @test status
    @test norm(ip.r, Inf) < opts.r_tol
    @test all([norm(ip.z[j][2:3]) <= ip.z[j][1] for i in ip.idx.socz for j in i])
    @test ip.κ[1] < opts.κ_tol
    @test norm(ip.δz, 1) != 0.0

    # problem setup
    v = [1.0; 10.0]
    μγ = 0.0
    θ = [v; μγ]

    z = ones(n)
    z[1] += 1.0
    z[5] += 1.0

    ip.z .= z
    ip.θ .= θ

    # solve
    status = RoboDojo.interior_point_solve!(ip)

    # test
    @test status
    @test norm(ip.r, Inf) < opts.r_tol
    @test all([norm(ip.z[j][2:3]) <= ip.z[j][1] for i in ip.idx.socz for j in i])
    @test ip.κ[1] < opts.κ_tol
    @test norm(ip.δz, 1) != 0.0

    # problem setup
    v = [0.0; 0.0]
    μγ = 0.0
    θ = [v; μγ]

    z = ones(n)
    z[1] += 1.0
    z[5] += 1.0

    ip.z .= z
    ip.θ .= θ

    # solve
    status = RoboDojo.interior_point_solve!(ip)

    # test
    @test status
    @test norm(ip.r, Inf) < opts.r_tol
    @test all([norm(ip.z[j][2:3]) <= ip.z[j][1] for i in ip.idx.socz for j in i])
    @test ip.κ[1] < opts.κ_tol
    @test norm(ip.δz, 1) != 0.0
end


text:
const XY = collect(range(-3, 3, length=100))

function showplot()
    GLMakie.activate!()
    fig = Figure()
    buildplot!(fig[2, 1], fig[3, 1], Makie.Slider(fig[1, 1], range=XY, startvalue=0))
    fig
end

function buildplot!(pa1, pa2, slider)
    ρ = 0.7

    dist = MvNormal([0, 0], [1 ρ; ρ 1])

    grid = zeros(Float64, (100, 100))

    for (x, i) in enumerate(XY)
        for (y, j) in enumerate(XY)
            grid[x, y] = pdf(dist, [i, j])
        end
    end

    # compute shading and integral of portion above cutoff
    integrated = lift(slider.value) do ŷ
        integrated = zeros(Float64, 100)

        for (x, i) in enumerate(XY)
            for (y, j) in enumerate(XY)
                # probit likelihood function: P(ŷ + ϵ > 0), ϵ ~ N(0, 1)
                #                               = P(ϵ > -ŷ)
                if i > -ŷ
                    integrated[y] += grid[x, y]
                end
            end
        end

        integrated = integrated ./ sum(integrated)

        integrated
    end

    color = lift(slider.value) do ŷ
        color = Array{Any}(undef, 100, 100)

        for (x, i) in enumerate(XY)
            for (y, j) in enumerate(XY)
                color[x, y] = i > -ŷ ? colorant"#4B9CD3" : colorant"#AAAAAA"
            end
        end

        color
    end

    vertices = lift(slider.value) do ŷ
        top = maximum(grid)
        vertices = [
            3 -3 0;  # 1
            3 3 0;  # 2
            -ŷ 3 0; # 3
            -ŷ -3 0; # 4
            3 -3 top;  # 5
            3 3 top;  # 6
            -ŷ 3 top; # 7
            -ŷ -3 top; # 8
        ]
    end

    faces = [
        1 2 3; 4 1 3; # bottom
        5 6 7; 8 7 5; # top
        1 2 5; 5 6 2;
        2 3 6; 7 6 3;
        3 4 7; 8 7 4;
        4 1 5; 8 5 4;
    ]

    ax1 = Axis3(pa1, xlabel="Selection error term", ylabel="Outcome error term")
    surface!(ax1, XY, XY, grid)
    mesh!(ax1, vertices, faces, color=colorant"rgba(75, 156, 211, 0.5)")
    ax2 = Axis(pa2, ylims=(0, 15), yticks=[0])
    lines!(ax2, XY, integrated)
end

function wglplot(s, x...)
    WGLMakie.activate!()

    App() do session::Session
        fig = Figure(resolution=(1200, 900))
        slider = Makie.Slider(fig[1, 1], range=XY, startvalue=0)
        buildplot!(fig[2, 1], fig[3, 1], slider)
        return JSServe.record_states(session, DOM.div(slider, fig))
    end
end

function exportplot()
    JSServe.export_standalone(wglplot, joinpath(Base.source_dir(), "output"))
end


text:
A = matrixdepot("toeplitz", [1,2,3,4,5])
@test A == matrixdepot("toeplitz", 5)

B = matrixdepot("toeplitz", [1,2,3], [1,4,5])
@test B[1,3] == 5
@test B[3,1] == 3
@test B[1,2] == 4
println("'toeplitz' passed test...")

n = rand(2:10)
A = matrixdepot("prolate", n)
@test issym(A)

println("'prolate' passed test...")


text:
# taking the outputted files for the variable DataFrame and writing them to a featherfile
function write_data(vars_results::Dict{Symbol, DataFrames.DataFrame}, save_path::AbstractString; kwargs...)
    file_type = get(kwargs, :file_type, Feather)
    if file_type == Feather || file_type == CSV
        for (k,v) in vars_results
            file_path = joinpath(save_path,"$(k).$(lowercase("$file_type"))")
            file_type.write(file_path, vars_results[k])
        end
    else
        error("unsupported file type: $file_type")
    end
    return
end

# taking the outputted files for the variable DataFrame and writing them to a featherfile
function write_data(results::Dict, save_path::AbstractString, variable::String;kwargs...)
    file_type = get(kwargs, :file_type, Feather)
    if file_type == Feather || file_type == CSV

        file_path = joinpath(save_path,"$(variable).$(lowercase("$file_type"))")
        file_type.write(file_path)

    else
        error("unsupported file type: $file_type")
    end
    return
end

function write_data(data::DataFrames.DataFrame, save_path::AbstractString, file_name::String; kwargs...)
    if isfile(save_path)
        save_path = dirname(save_path)
    end
    file_type = get(kwargs, :file_type, Feather)
    if file_type == Feather || file_type == CSV
        file_path = joinpath(save_path,"$(file_name).$(lowercase("$file_type"))")
        file_type.write(file_path, data)
    else
        error("unsupported file type: $file_type")
    end
    return
end

function _write_optimizer_log(optimizer_log::Dict, save_path::AbstractString)

    optimizer_log[:dual_status] = Int(optimizer_log[:dual_status])
    optimizer_log[:termination_status] = Int(optimizer_log[:termination_status])
    optimizer_log[:primal_status] = Int(optimizer_log[:primal_status])
    optimizer_log = DataFrames.DataFrame(optimizer_log)
    file_path = joinpath(save_path,"optimizer_log.feather")
    Feather.write(file_path, optimizer_log)

    return

end

function _write_data(canonical::CanonicalModel, save_path::AbstractString; kwargs...)
    file_type = get(kwargs, :file_type, Feather)
    if file_type == Feather || file_type == CSV
        for (k, v) in vars(canonical)
            file_path = joinpath(save_path,"$(k).$(lowercase("$file_type"))")
            file_type.write(file_path, _result_dataframe_vars(v))
        end
    else
        error("unsupported file type: $file_type")
    end
    return
end

function _write_data_(canonical::CanonicalModel, save_path::AbstractString; kwargs...)
    file_type = get(kwargs, :file_type, Feather)
    if file_type == Feather || file_type == CSV
        (k, v) =  (:Price,con(canonical,:CopperPlateBalance))
        file_path = joinpath(save_path,"$(k).$(lowercase("$file_type"))")

        file_type.write(file_path, _result_dataframe_duals(v))
    else
        error("unsupported file type: $file_type")
    end
    return
end

function write_data(op_model::OperationModel, save_path::AbstractString; kwargs...)
    _write_data(op_model.canonical, save_path; kwargs...)
    return
end

function write_data(stage::_Stage, save_path::AbstractString; kwargs...)
    if stage.key == 2
        _write_data_(stage.canonical, save_path; kwargs...)
    end
    _write_data(stage.canonical, save_path; kwargs...)
    return
end

# These functions are writing directly to the feather file and skipping printing to memory.
function _export_model_result(stage::_Stage, start_time::Dates.DateTime, save_path::String)
    write_data(stage, save_path)
    write_data(get_time_stamp(stage, start_time), save_path, "time_stamp")
    return
end

function _export_optimizer_log(optimizer_log::Dict{Symbol, Any},
                               canonical_model::CanonicalModel,
                               path::String)

    optimizer_log[:obj_value] = JuMP.objective_value(canonical_model.JuMPmodel)
    optimizer_log[:termination_status] = Int(JuMP.termination_status(canonical_model.JuMPmodel))
    optimizer_log[:primal_status] = Int(JuMP.primal_status(canonical_model.JuMPmodel))
    optimizer_log[:dual_status] = Int(JuMP.dual_status(canonical_model.JuMPmodel))
    try
        optimizer_log[:solve_time] = MOI.get(canonical_model.JuMPmodel, MOI.SolveTime())
    catch
        @warn("SolveTime() property not supported by the Solver")
        optimizer_log[:solve_time] = NaN #"Not Supported by solver"
    end
    _write_optimizer_log(optimizer_log, path)
    return
end

function write_model_results(res::OperationModelResults, path::String, results::String)
    if !isdir(path)
        @error("Specified path is not valid. Run write_results to save results.")
    end
    folder_path = joinpath(path,results)
    write_data(res.variables, folder_path)
    write_data(res.optimizer_log, folder_path, "optimizer_log")
    time = DataFrames.DataFrame(Range = convert.(Dates.DateTime,res.time_stamp[!,:Range]))
    write_data(time, folder_path, "time_stamp")
    println("Files written to $folder_path folder.")
    return
end

""" Exports Operational Model Results to a path"""
function write_model_results(results::OperationModelResults, save_path::String)
    if !isdir(save_path)
        @error("Specified path is not valid. Run write_results to save results.")
    end
    _new_folder_path = replace_chars("$save_path/$(round(Dates.now(),Dates.Minute))", ":", "-")
    new_folder = mkdir(_new_folder_path)
    folder_path = new_folder
    write_data(results.variables, folder_path)
    _write_optimizer_log(results.optimizer_log, folder_path)
    write_data(results.time_stamp, folder_path, "time_stamp")
    println("Files written to $folder_path folder.")
    return
end

""" Exports the OpModel JuMP object in MathOptFormat"""
function write_op_model(op_model::OperationModel, save_path::String)
    _write_canonical_model(op_model.canonical, save_path)
    return
end

""" Exports the OpModel JuMP object in MathOptFormat"""
function _write_canonical_model(canonical_model::CanonicalModel, save_path::String)
    MOF_model = MOPFM
    MOI.copy_to(MOF_model, JuMP.backend(canonical_model.JuMPmodel))
    MOI.write_to_file(MOF_model, save_path)
    return
end


text:
#=
  Copyright (c) 2021
  https://github.com/Ismael-VC/Varvara.jl/blob/main/CONTRIBUTORS.md

  Permission to use, copy, modify, and distribute this software for any
  purpose with or without fee is hereby granted, provided that the above
  copyright notice and this permission notice appear in all copies.

  THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
  WITH REGARD TO THIS SOFTWARE.
=#

module UTOS

function main()::Int
  if length(ARGS) > 2 || !isfile(ARGS[1])
    return 1
  end

  buffer = open(ARGS[1]) do io
    read(io)
  end

  buffer .+= 0x80

  open(ARGS[2], "w") do io
    write(io, buffer)
  end

  println("\n")
  return 0
end

end  # module

if abspath(PROGRAM_FILE) == @__FILE__
  using .UTOS: main

  main()
end


text:
using BAT
using Compat.Test

using Distributions, PDMats

@testset "density_product" begin
    econtext = @inferred ExecContext()
    mvec = [-0.3, 0.3, 0.1]
    cmat = [1.0 1.5 0.0; 1.5 4.0 0.0; 0.0 0.0 1.0]
    Σ = @inferred PDMat(cmat)
    mvnorm = @inferred  MvNormal(mvec, Σ)
    mvn_density = @inferred GenericDensity(params -> logpdf(mvnorm, params), 3)

    cmat = [3.76748 0.446731 0.625418; 0.446731 3.9317 0.237361; 0.625418 0.237361 3.43867]
    mvec = [1., 2, 1.]
    Σ = @inferred PDMat(cmat)
    mvt = MvTDist(3, mvec, Σ)
    mvt_density = @inferred GenericDensity(params -> logpdf(mvt, params), 3)

    params = [0.0 -0.3; 0.0 0.3; 0.0 1.0]

    pb = BAT.HyperRectBounds([-2.0, -1.0, -0.5], [2.0, 3.0, 1.0],
                                reflective_bounds)
    pb2 = BAT.HyperRectBounds([-1.5, -2.0, -0.5], [2.0, 2.5, 1.5],
                                reflective_bounds)
    dp = @inferred BAT.DensityProduct((mvt_density, mvn_density), pb)
    dp1 = @inferred BAT.DensityProduct((mvt_density,), pb)
    dp2 = @inferred BAT.DensityProduct((mvn_density,), pb2)

    res = Array{Float64}(size(params, 2))

    @testset "DensityProduct" begin
        @test typeof(dp) <: BAT.DensityProduct{2,
            Tuple{BAT.GenericDensity{typeof(mvt_density.log_f)},
                  BAT.GenericDensity{typeof(mvn_density.log_f)}},BAT.HyperRectBounds{Float64}}

        @test parent(dp)[1] == mvt_density
        @test param_bounds(dp) == pb
        @test nparams(dp) == 3
    end

    @testset "unsafe_pod" begin
        up = @inferred BAT._unsafe_prod(mvt_density, mvn_density, pb)
        @test parent(up)[2] == mvn_density
        @test param_bounds(up) == pb

        prd = @inferred dp1*dp2
        @test typeof(prd) <: BAT.DensityProduct{2,
            Tuple{BAT.GenericDensity{typeof(mvt_density.log_f)},
                  BAT.GenericDensity{typeof(mvn_density.log_f)}},BAT.HyperRectBounds{Float64}}
        @test parent(prd)[2] == mvn_density
        prd_pb = @inferred param_bounds(prd)
        cut_pb = @inferred pb ∩ pb2
        @test prd_pb.vol.lo ≈ cut_pb.vol.lo
        @test prd_pb.vol.hi ≈ cut_pb.vol.hi

        up = @inferred BAT._unsafe_prod(dp1, mvn_density, pb)
        @test parent(up)[1] == mvt_density
        @test parent(up)[2] == mvn_density
        @test param_bounds(up) == pb

        up = @inferred BAT._unsafe_prod(mvt_density, dp2, pb)
        @test parent(up)[1] == mvt_density
        @test parent(up)[2] == mvn_density
        @test param_bounds(up) == pb

        up = @inferred BAT._unsafe_prod(dp1, dp2, pb)
        @test parent(up)[1] == mvt_density
        @test parent(up)[2] == mvn_density
        @test param_bounds(up) == pb
    end

    @testset "unsafe_density_logval" begin
        @test BAT.unsafe_density_logval(dp, params[:,1], econtext) ≈ -8.8547305

        @test_throws ArgumentError BAT.unsafe_density_logval(BAT.DensityProduct(Tuple([]), pb),
                                                             params[:,1], econtext)

        BAT.unsafe_density_logval!(res, dp, params, econtext)

        @test res ≈ [-8.8547305, -8.8634491]
    end

    @testset "ExecCapabilities" begin
        ec = @inferred BAT.exec_capabilities(BAT.unsafe_density_logval, dp, params[:, 1])
        ec_tocmp = @inferred ∩(
            BAT.exec_capabilities(BAT.unsafe_density_logval, mvt_density, params[:, 1]),
            BAT.exec_capabilities(BAT.unsafe_density_logval, mvn_density, params[:, 1]))

        @test ec.nthreads == ec_tocmp.nthreads
        @test ec.threadsafe == ec_tocmp.threadsafe
        @test ec.nprocs == ec_tocmp.nprocs
        @test ec.remotesafe == ec_tocmp.remotesafe

        ec = @inferred BAT.exec_capabilities(BAT.unsafe_density_logval!, res, dp, params)
        ec_tocmp = @inferred ∩(
            BAT.exec_capabilities(BAT.unsafe_density_logval!, res, mvt_density, params),
            BAT.exec_capabilities(BAT.unsafe_density_logval!, res, mvn_density, params))

        @test ec.nthreads == ec_tocmp.nthreads
        @test ec.threadsafe == ec_tocmp.threadsafe
        @test ec.nprocs == ec_tocmp.nprocs
        @test ec.remotesafe == ec_tocmp.remotesafe

    end
end


text:
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("colorscheme.jl"))
using DynamicalSystems, InteractiveDynamics
import GLMakie
using AbstractPlotting
using LinearAlgebra
using Statistics

# Code for 3D animation
ds = Systems.lorenz()
u0 = trajectory(ds, 10; Ttr = 100)[end]
dt = 0.01
T = 5.0
framerate = 30
systemtitle = "Lorenz63"

ds = Systems.towel()
u0 = trajectory(ds, 10; Ttr = 100)[end]
dt = 1
T = 15 # total evolove time
framerate = 2
systemtitle = "towelmap"

# %%
ε = 1e-6 # sphere radius
N = 4000 # initial conditions nmaking hypersphere

function uniform_sphere(N, u0, ε)
    V = typeof(u0)
    us = Vector{V}(undef, N)
    i = 1
    while i ≤ N
        x = rand(V) .- 0.5
        if norm(x) ≤ 1/(2*√3)
            us[i] = u0 .+ (ε .* x)
            i += 1
        end
    end
    return us
end

us = uniform_sphere(N, u0, ε)
# I must renormalize point coordinates to a sphere of radius 1 always,
# because of the internal convertion to 32-bit accuracy, small numbers are bad!
usnormed = [Point3f0(2((u0 .- u)./ε)...,) for u in us]
colors = [RGBAf0(abs.(u)..., 1) for u in usnormed]
plotus = Observable(usnormed)

# fig = Figure(resolution = (1200, 600))
fig = Figure(resolution = (800, 600))
ax3D = fig[1,1] = Axis3(fig)
scatter!(ax3D, plotus; markersize = 2000, color = colors, strokewidth = 0.5)
# ax2 = fig[1, 2] = Axis(fig; width = 400) # this axis has values of principal lengths
# ax2.title = "logarithm of principal axes"
ax3D.xticklabelsvisible = true
ax3D.yticklabelsvisible = true
ax3D.zticklabelsvisible = true
ax3D.azimuth = 4.44
ax3D.elevation = 0.38
ax3D.title = rpad(systemtitle *", t = $(0)", 40, ' ')

# function principals(v)
#     M = transpose(reshape(reinterpret(Float64, v), 3, :))
#     S = SVD(M).S # rotation matrix
#     return S
# end
# s = Observable(principals(usnormed))
# barplot!(ax2, s, color = to_color.(COLORSCHEME[1:3]))

display(fig)

# %%
ut = trajectory(ds, T, u0; dt)
evolution = [trajectory(ds, T, u; dt) for u in us]
t = 0:dt:T

# TODO: Use absolute, realsize of pointcloud

record(
        fig, joinpath(@__DIR__, "volume_growth_$(systemtitle).mp4"), 1:length(evolution[1]);
        framerate
        ) do i
    newus = [Point3f0(2((ut[i] .- e[i])./ε)...,) for e in evolution]
    plotus[] = newus
    # s[] = principals(newus)
    # autolimits!(ax2)
    # must set new limits centered on u
    maxε = maximum(norm(u) for u in newus)
    ax3D.limits = map(j -> (-maxε, +maxε), (1,2,3))
    ax3D.title = systemtitle *", t = $(t[i])"
    ax3D.elevation = 0.38
    ax3D.azimuth = 4.44
end

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

module HCatTest
using ParallelAccelerator

#ParallelAccelerator.DomainIR.set_debug_level(4)
#ParallelAccelerator.ParallelIR.set_debug_level(4)
#ParallelAccelerator.CGen.set_debug_level(4)
#ParallelAccelerator.set_debug_level(4)
#using CompilerTools
#CompilerTools.OptFramework.set_debug_level(4)

@acc function cat1(a::Array{Float64,1},b::Array{Float64,1})
    C = hcat(a,b)
    return C
end

function test1()
    return cat1([1.,2.,3.],[4.,5.,6.])
end
end

using Base.Test
println("Testing hcat...")
@test_approx_eq HCatTest.test1() hcat([1.,2.,3.],[4.,5.,6.])
println("Done testing hcat.")


text:
module ExecutableSpecifications

include("Gherkin.jl")

"Abstraction for presenting results from scenario steps."
abstract type Presenter end

"Presenting results from scenario steps as they occur."
abstract type RealTimePresenter <: Presenter end

include("stepdefinitions.jl")
include("exec_env.jl")
include("executor.jl")
include("asserts.jl")
include("presenter.jl")
include("result_accumulator.jl")
include("engine.jl")
include("runner.jl")

export @given, @when, @then, @expect, @beforescenario, @afterscenario, runspec

end # module


text:
"""
    Chisq(ν)
The *Chi squared distribution* (typically written χ²) with `ν` degrees of freedom has the
probability density function

```math
f(x; k) = \\frac{x^{k/2 - 1} e^{-x/2}}{2^{k/2} \\Gamma(k/2)}, \\quad x > 0.
```

If `ν` is an integer, then it is the distribution of the sum of squares of `ν` independent standard [`Normal`](@ref) variates.

```julia
Chisq(k)     # Chi-squared distribution with k degrees of freedom

params(d)    # Get the parameters, i.e. (k,)
dof(d)       # Get the degrees of freedom, i.e. k
```

External links

* [Chi-squared distribution on Wikipedia](http://en.wikipedia.org/wiki/Chi-squared_distribution)
"""
struct Chisq{T<:Real} <: ContinuousUnivariateDistribution
    ν::T
    Chisq{T}(ν::T) where {T} = new{T}(ν)
end

function Chisq(ν::T) where {T <: Real}
    @check_args(Chisq, ν > zero(ν))
    return Chisq{T}(ν)
end

Chisq(ν::T, ::NoArgCheck) where {T <: Real} = Chisq{T}(ν)
Chisq(ν::Integer) = Chisq(float(ν))

@distr_support Chisq 0.0 Inf

#### Parameters

dof(d::Chisq) = d.ν
params(d::Chisq) = (d.ν,)
@inline partype(d::Chisq{T}) where {T<:Real} = T

### Conversions
convert(::Type{Chisq{T}}, ν::Real) where {T<:Real} = Chisq(T(ν))
convert(::Type{Chisq{T}}, d::Chisq{S}) where {T <: Real, S <: Real} = Chisq(T(d.ν))

#### Statistics

mean(d::Chisq) = d.ν

var(d::Chisq) = 2d.ν

skewness(d::Chisq) = sqrt(8 / d.ν)

kurtosis(d::Chisq) = 12 / d.ν

mode(d::Chisq{T}) where {T<:Real} = d.ν > 2 ? d.ν - 2 : zero(T)

function median(d::Chisq; approx::Bool=false)
    if approx
        return d.ν * (1 - 2 / (9 * d.ν))^3
    else
        return quantile(d, 1//2)
    end
end

function entropy(d::Chisq)
    hν = d.ν/2
    hν + logtwo + lgamma(hν) + (1 - hν) * digamma(hν)
end

#### Evaluation

@_delegate_statsfuns Chisq chisq ν

mgf(d::Chisq, t::Real) = (1 - 2 * t)^(-d.ν/2)

cf(d::Chisq, t::Real) = (1 - 2 * im * t)^(-d.ν/2)

gradlogpdf(d::Chisq{T}, x::Real) where {T<:Real} =  x > 0 ? (d.ν/2 - 1) / x - 1//2 : zero(T)

#### Sampling

rand(rng::AbstractRNG, d::Chisq) =
    (ν = d.ν; rand(rng, Gamma(ν / 2.0, 2.0one(ν))))

sampler(d::Chisq) = (ν = d.ν; sampler(Gamma(ν / 2.0, 2.0one(ν))))


text:
# TODO Move this to GPLikelihoods.jl
@doc raw"""
    LaplaceLikelihood(β::Real)

Likelihood with a [Laplace distribution](https://en.wikipedia.org/wiki/Laplace_distribution):
```math
    p(y|f,\beta) = \frac{1}{2\beta}\exp\left(-\frac{|y-f|}{\beta}\right)
```

## Arguments
- `β::Real`, scale parameter
"""
struct LaplaceLikelihood{Tβ<:Real} <: AbstractLikelihood
    β::Tβ
end

LaplaceLikelihood() = LaplaceLikelihood(1.0)

(lik::LaplaceLikelihood)(f::Real) = Laplace(f, lik.β)

function (lik::LaplaceLikelihood)(f::AbstractVector{<:Real})
    return Product(lik.(f))
end

@inline laplace_λ(lik::LaplaceLikelihood) = inv(abs2(2 * lik.β))

aux_field(::LaplaceLikelihood, Ω) = getproperty(Ω, :ω)

function init_aux_variables(rng::AbstractRNG, ::LaplaceLikelihood, n::Int)
    return TupleVector((; ω=rand(rng, InverseGamma(), n)))
end

function init_aux_posterior(T::DataType, lik::LaplaceLikelihood, n::Int)
    λ = laplace_λ(lik)
    return For(TupleVector(; μ=zeros(T, n))) do φ
        NTDist(InverseGaussian(φ.μ, λ))
    end
end

function aux_full_conditional(lik::LaplaceLikelihood, y::Real, f::Real)
    return NTDist(InverseGaussian(inv(2 * lik.β * abs(y - f)), 2 * laplace_λ(lik)))
end

function aux_posterior!(
    qΩ, lik::LaplaceLikelihood, y::AbstractVector, qf::AbstractVector{<:Normal}
)
    φ = only(qΩ.inds)
    map!(φ.μ, y, qf) do yᵢ, qfᵢ
        inv(2 * lik.β * sqrt(second_moment(qfᵢ, yᵢ)))
    end
    return qΩ
end

function auglik_potential(::LaplaceLikelihood, Ω, y::AbstractVector)
    return (2 * Ω.ω .* y,)
end

function auglik_precision(::LaplaceLikelihood, Ω, ::AbstractVector)
    return (2 * Ω.ω,)
end

function expected_auglik_potential(::LaplaceLikelihood, qΩ, y::AbstractVector)
    return (2 * tvmean(qΩ).ω .* y,)
end

function expected_auglik_precision(::LaplaceLikelihood, qΩ, ::AbstractVector)
    return (2 * tvmean(qΩ).ω,)
end

function logtilt(
    lik::LaplaceLikelihood, Ω::TupleVector, y::AbstractVector, f::AbstractVector{<:Real}
)
    return length(y) * (loggamma(1//2) - log(sqrtπ) - log(2 * lik.β)) +
           mapreduce(+, aux_field(lik, Ω), y, f) do ωᵢ, yᵢ, fᵢ
        -abs2(yᵢ - fᵢ) * ωᵢ
    end
end

function logtilt(lik::LaplaceLikelihood, ω::Real, y::Real, f::Real)
    return loggamma(1//2) - log(sqrtπ) - log(2 * lik.β) - abs2(y - f) * ω
end

function expected_logtilt(lik::LaplaceLikelihood, qΩ, y, qf::AbstractVector{<:Normal})
    return length(y) * (loggamma(1//2) - log(sqrtπ) - log(2 * lik.β)) +
           mapreduce(+, y, qf, @ignore_derivatives marginals(qΩ)) do yᵢ, qfᵢ, qωᵢ
        -second_moment(qfᵢ, yᵢ) * ntmean(qωᵢ).ω
    end
end

function aux_prior(lik::LaplaceLikelihood, y::AbstractVector{<:Real})
    return For(length(y)) do _
        NTDist(aux_prior(lik))
    end
end

aux_prior(lik::LaplaceLikelihood) = InverseGamma(1//2, laplace_λ(lik))

function aux_kldivergence(::LaplaceLikelihood, qΩ::For, pΩ::For)
    return mapreduce(+, marginals(qΩ), marginals(pΩ)) do qωᵢ, pωᵢ
        μ = mean(dist(qωᵢ))
        λ = scale(dist(pωᵢ))
        log(2λ) / 2 - log(2π) / 2 - log(λ) / 2 + loggamma(1//2) + λ / μ
    end
end


text:
# taken from https://diffeq.sciml.ai/stable/features/noise_process/

using StochasticDiffEq,  DiffEqNoiseProcess
using Test

f1(u, p, t) = 1.01u
g1(u, p, t) = 1.01u
dt = 1//2^(4)
prob1 = SDEProblem(f1,g1,1.0,(0.0,1.0))
sol1 = solve(prob1,EM(),dt=dt,save_noise = true)

W2 = NoiseWrapper(sol1.W)
prob1 = SDEProblem(f1,g1,1.0,(0.0,1.0),noise=W2)
sol2 = solve(prob1,EM(),dt=dt)

@test sol1.u ≈ sol2.u

W3 = NoiseWrapper(sol1.W)
prob2 = SDEProblem(f1,g1,1.0,(0.0,1.0),noise=W3)

# smaller timestep
dt = 1//2^(5)
sol3 = solve(prob2,EM(),dt=dt)

using Plots
plot(sol1)
plot!(sol2)
plot!(sol3)