text:
module PMIx
using  Setfield

include("api.jl")

struct PMIxException <: Exception
    status::API.pmix_status_t
end

macro check(ex)
    quote
        status = $(esc(ex))
        if status != API.PMIX_SUCCESS
            throw(PMIxException(status))
        end
    end
end

function pmix_strncopy(dst::Ptr{Cchar}, src::Ptr{Cchar}, len)
    i = 1
    while i <= len
        c = unsafe_load(src, i)
        unsafe_store!(dst, i, c)
        if Char(c) == '\0'
            break
        end
        i += 1
    end
    unsafe_store!(dst, i, '\0'%Cchar)
end

include("value.jl")

function ZeroInfo()
    # Equivalent to PMIX_INFO_CONSTRUCT
    r_info = Ref{API.pmix_info_t}()
    ccall(:memset, Ptr{Cvoid}, (Ptr{Cvoid}, Cint, Csize_t), r_info, 0, sizeof(API.pmix_info_t))
    return r_info[]
end

function Info()
    info = ZeroInfo()
    info = @set info.value = UndefValue()
    return info
end

function Info(key, value, flags = 0)
    r_key = Ref{API.pmix_key_t}()
    GC.@preserve r_key key begin
        cstr = Base.unsafe_convert(Cstring, key) # ensure null termination
        p_src = convert(Ptr{Cchar}, cstr)
        p_dst = Base.unsafe_convert(Ptr{Cchar}, r_key)
        pmix_strncopy(p_dst, p_src, API.PMIX_MAX_KEYLEN)
    end
    API.pmix_info(r_key[], flags, value)
end

function Proc()
    # Equivalent to PMIX_PROC_CONSTRUCT
    r_proc = Ref{API.pmix_proc_t}()
    ccall(:memset, Ptr{Cvoid}, (Ptr{Cvoid}, Cint, Csize_t), r_proc, 0, sizeof(API.pmix_proc_t))
    return r_proc[]
end

function Proc(nspace, rank)
    proc = Proc()
    proc = @set proc.nspace = nspace
    proc = @set proc.rank = rank
    return proc
end

function nspace(proc::API.pmix_proc_t)
    ns = Ref(proc.nspace)
    GC.@preserve ns begin
        return unsafe_string(Base.unsafe_convert(Ptr{Cchar}, ns))
    end
end

# 4. Client initialization and finalization

function initialized()
    API.PMIx_Initialized() != 0
end

function version()
    return Base.unsafe_string(API.PMIx_Get_version())
end

# TODO: Parameters
function init()
    r_proc = Ref{API.pmix_proc_t}()
    @check API.PMIx_Init(r_proc, C_NULL, 0)
    return r_proc[]
end

function finalize()
    @check API.PMIx_Finalize(C_NULL, 0)
end

function progress()
    API.PMIx_Progress()
end

# 5. Synchronization and Data Access Operations

function fence(procs)
    @check API.PMIx_Fence(procs, length(procs), C_NULL, 0)
end

## PMIx_Fence_nb

function get(proc, key, info=nothing)
    if info === nothing
        info = C_NULL
        len = 0
    else
        len= length(info)
    end
    r_ptr = Ref{Ptr{API.pmix_value_t}}()
    @check API.PMIx_Get(Ref(proc), key, info, len, r_ptr)
    ptr = r_ptr[]
    value = Base.unsafe_load(ptr)
    # TODO: Complain about pmix_free being defined inline
    # TODO: Memory rules?
    Libc.free(ptr)
    return value
end

# 5.4 Query

# 7. Process-Related Non-Reserved Keys'

function put!(scope, key, value)
    @check API.PMIx_Put(scope, key, value)
end

# PMIx_Store_internal

function commit()
    @check API.PMIx_Commit()
end

# 8. Publish/Lookup Operations
# 9. Event Notification
# 10. Data Packing and Unpacking
# 11. Process Management

end # module

max_issues_repo_name: JuliaParallel/PMIx

max_forks_repo_name: JuliaParallel/PMIx


text:
abstract type AbstractSearcher{T} end

"""
    Search object. It contains all the indexed data and related
configuration that allows for searches to be performed.
"""
struct Searcher{T<:AbstractFloat,
                E<:AbstractEmbedder{T},
                I<:AbstractIndex} <: AbstractSearcher{T}
    data::Ref
    config::NamedTuple
    input_embedder::Ref{E}                      # embeds queries
    data_embedder::Ref{E}                       # embeds data
    index::I                                    # indexed search data
    search_trees::BKTree{String}                # suggestion structure
end

# Useful methods
id(srcher::Searcher) = srcher.config.id

description(srcher::Searcher) = srcher.config.description

isenabled(srcher::Searcher) = first(srcher.config.enabled)

disable!(srcher::Searcher) = begin
    srcher.config.enabled[1] = false
    nothing
end

enable!(srcher::Searcher) = begin
    srcher.config.enabled[1] = true
    nothing
end

# push!, pushfirst!, pop!, popfirst!, delete!
Base.push!(srcher::Searcher, entry) = pushinner!(srcher, entry, :last)

Base.pushfirst!(srcher::Searcher, entry) = pushinner!(srcher, entry, :first)

pushinner!(srcher::Searcher, entry, position::Symbol) = begin
    embedded, _ = embed(srcher.data_embedder[], [entry]; fields=srcher.config.indexable_fields)
    index_operation = ifelse(position === :first, pushfirst!, push!)
    index_operation(srcher.index, firstcol(embedded))  # push first column i.e. vector
    nothing
end

Base.pop!(srcher::Searcher) = pop!(srcher.index)

Base.popfirst!(srcher::Searcher) = popfirst!(srcher.index)

Base.deleteat!(srcher::Searcher, pos) = delete_from_index!(srcher.index, pos)

# Indexing for vectors of searchers
function getindex(srchers::AbstractVector{Searcher}, an_id::String)
    idxs = Int[]
    for (i, srcher) in enumerate(srchers)
        isequal(id(srcher), an_id) && push!(idxs, i)
    end
    return srchers[idxs]
end

"""
    build_searcher(dbdata, config)

Creates a Searcher from a searcher configuration.
"""
function build_searcher(dbdata, embedders, config; id_key=DEFAULT_DB_ID_KEY)

    # Select embedders
    input_embedder = first(filter(embdr->embdr.config.id==config.input_embedder, embedders))
    data_embedder = first(filter(embdr->embdr.config.id==config.data_embedder, embedders))

    # Embed db entries
    entries = db_sorted_row_iterator(dbdata; id_key=id_key, rev=false)
    embedded, _ = embed(data_embedder, entries; fields=config.indexable_fields)

    # Build search index
    indexer= build_indexer(config.search_index,
                           config.search_index_arguments,
                           config.search_index_kwarguments)
    srchindex = indexer(embedded)

    # Build search tree (for suggestions)
    srchtree = build_bktree(dbdata, config.heuristic; id_key=id_key)

    # Build searcher
    srcher = Searcher(Ref(dbdata),
                      config,
                      Ref(input_embedder),
                      Ref(data_embedder),
                      srchindex,
                      srchtree)

    @debug "* Loaded: $srcher."
    return srcher
end

function build_bktree(dbdata, heuristic; id_key=nothing)
    if heuristic != nothing
        documents = [join(dbentry2text(dbentry, config.indexable_fields), " ") # merge sentences
                     for dbentry in db_sorted_row_iterator(dbdata; id_key=id_key, rev=false)]
        lexicon = create_lexicon(documents, 1)
        distance = get(HEURISTIC_TO_DISTANCE, heuristic, DEFAULT_DISTANCE)
        fdist = (x,y) -> evaluate(distance, x, y)
        return BKTree(fdist, collect(keys(lexicon)))
    else
        return BKTree{String}()
    end
end

# Supported indexes name to type mapping
function build_indexer(index, args, kwargs)
    default_hnsw_kwarguments = (:efConstruction=>100, :M=>16, :ef=>50)  # to ensure it works well
    default_ivfadc_kwarguments = (:kc=>2, :k=>2, :m=>1)  # to ensure it works at all
    index === :naive && return d->NaiveIndex(d, args...; kwargs...)
    index === :brutetree && return d->BruteTreeIndex(d, args...; kwargs...)
    index === :kdtree && return d->KDTreeIndex(d, args...; kwargs...)
    index === :hnsw && return d->HNSWIndex(d, args...; default_hnsw_kwarguments..., kwargs...)
    index === :ivfadc && return d->IVFIndex(d, args...; default_ivfadc_kwarguments..., kwargs...)
    index === :noop && return d->NoopIndex(d, args...; kwargs...)
end

ext: jl


text:
abstract type EoSModel end
abstract type LatticeFluidModel <: EoSModel end
abstract type SAFTModel <: LatticeFluidModel end
abstract type GCSAFTModel <: SAFTModel end
abstract type CubicModel <: EoSModel end
abstract type ABCubicModel <: CubicModel end #cubics that have an exact polynomial form to solve Z roots, this excludes CPA
abstract type ActivityModel <: EoSModel end
abstract type IdealModel <: EoSModel end
abstract type EmpiricHelmholtzModel <: EoSModel end
abstract type EoSParam end
abstract type SatPureAproximation <: EoSModel end
abstract type AlphaModel <:EoSModel end
export EoSModel
export SAFTModel,CubicModel,EmpiricHelmholtzModel
export IdealModel
export EoSParam
export AlphaModel


text:
"""
This is Tutorial_CSV.jl. It contains all the necessary commands to import a CSV file from a seismic tomography,
to convert it to the GeoData format and to export it to a Paraview format. The following steps are performed:
1. Read data from this file.
2. Put the data in a GeoData format (this is the format that is used internally in the GMG).
3. Export that data to a format readable by Paraview.

You will need to download the data from:
[https://www.seismologie.ifg.uni-kiel.de/en/research/research-data/mere2020model](https://www.seismologie.ifg.uni-kiel.de/en/research/research-data/mere2020model).
"""

using DelimitedFiles, GeophysicalModelGenerator, GeoStats

# 1. define the filename and path
filename = "./El-Sharkawy-etal-G3.2020_MeRE2020_Mediterranean.csv" # this only works if you are in the directory where the data is located

# 2. load desired data
data        =   readdlm(filename,'|',Float64,'\n', skipstart=23,header=false)
lat         =   data[:,1];
lon         =   data[:,2];
depth       =  -data[:,3];
Vs          =   data[:,4];

# Create 3D regular grid:
Depth_vec       =   unique(depth)
Lon,Lat,Depth   =   LonLatDepthGrid(-10:0.5:40,32:0.25:50,Depth_vec);

# Employ GeoStats to interpolate irregular data points to a regular grid
dLon = Lon[2,1,1]-Lon[1,1,1]
dLat = Lat[1,2,1]-Lat[1,1,1]
Cgrid = CartesianGrid((size(Lon,1),size(Lon,2)),(minimum(Lon),minimum(Lat)),(dLon,dLat))
Vs_3D = zeros(size(Depth));
for iz=1:size(Depth,3)
    println("Depth = $(Depth[1,1,iz])")
    ind   = findall(x -> x==Depth[1,1,iz], depth)
    coord = PointSet([lon[ind]'; lat[ind]'])
    Geo   = georef((Vs=Vs[ind],), coord)
    P     = EstimationProblem(Geo, Cgrid, :Vs)
    S     = IDW(:Vs => (distance=Euclidean(),neighbors=2));
    sol   = solve(P, S)
    sol_Vs= values(sol).Vs
    Vs_2D = reshape(sol_Vs, size(domain(sol)))
    Vs_3D[:,:,iz] = Vs_2D;
end

# Save data to paraview:
Data_set    =   GeophysicalModelGenerator.GeoData(Lon,Lat,Depth,(Vs_km_s=Vs_3D,))   # the GeoStats package defines its own GeoData structure, so you have to choose the correct one here
Write_Paraview(Data_set, "MeRe_ElSharkawy")


text:
__precompile__()
module VisualDL

using PyCall
const LogWriter = PyNULL()

function __init__()
    copy!(LogWriter, pywrap(pyimport("visualdl")).LogWriter)
end

include("logger.jl")

end # module


text:
# May be really useful for category theoretic constructs

struct Arrow{X,Y}
   fun::Function
end

(-->)(a,b) = Arrow{a,b}

Int --> Int
# Arrow{Int64, Int64}

f = (Int --> Int)((x) -> x*2)

# allows for currying
Base.operator_associativity(:-->) #:right

(a::Arrow{X,Y})(x) where X where Y = begin
   !(x isa X) && error("type error")
   y = a.fun(x)
   !(y isa Y) && error("type error", typeof(y))
   return y
end

f(2) # 4

dom(a::Arrow{X,Y}) where X where Y = X
cod(a::Arrow{X,Y}) where X where Y = Y

g = (Int --> Int --> Int)((x) ->
   (Int --> Int)((y) -> (x * y)))

g(3)(4)


text:
function test_tied_inner_product_layer(backend::Backend, n_input, T, eps)
  println("-- Testing TiedInnerProductLayer on $(typeof(backend)){$T}...")

  batch_size = abs(rand(Int)) % 10 + 5
  orig_dim = 20
  hidden_dim = 10

  inputs = Array[rand(T, orig_dim, batch_size) for i=1:n_input]
  b1 = rand(T, hidden_dim, 1)
  b2 = rand(T, orig_dim, 1)
  param_key = "shared-ip"

  println("    > Setup")
  registry_reset(backend)
  layer_data = MemoryDataLayer(tops=[symbol("i-$i") for i = 1:n_input],
      batch_size=batch_size, data=inputs)
  layer_ip = InnerProductLayer(name="ip1", param_key=param_key,
      tops=[symbol("ip1-$i") for i=1:n_input], bottoms=[symbol("i-$i") for i=1:n_input],
      output_dim=hidden_dim)
  layer_ip2 = TiedInnerProductLayer(name="ip2", tied_param_key=param_key,
      tops=[symbol("ip2-$i") for i=1:n_input], bottoms=[symbol("ip1-$i") for i=1:n_input])

  net = Net("test-tied-ip", backend, [layer_data,layer_ip,layer_ip2])
  init(net)
  copy!(net.states[2].b, b1)
  copy!(net.states[3].b, b2)
  W = to_array(net.states[2].W)

  println("    > Forward")
  forward(net)
  for i = 1:n_input
    expected_output = W * (W' * inputs[i] .+ b1) .+ b2
    got_output = to_array(net.states[3].blobs[i])
    @test all(abs(expected_output - got_output) .< eps)
  end

  println("    > Backward")
  # This network is invalid for back-propagation because there is
  # no top layer that compute loss
  @test_throws TopologyError check_bp_topology(net)

  # however, we will manually set the diffs for the top-most layer
  # to do a fake back-propagate test
  top_diffs = Array[rand(T, orig_dim, batch_size) for i=1:n_input]
  for i = 1:n_input
    copy!(net.states[3].blobs_diff[i], top_diffs[i])
  end
  backward(net)

  bias_grad = to_array(net.states[3].∇b)
  bias_grad_expected = sum([sum(top_diffs[i],2) for i = 1:n_input])
  @test all(abs(bias_grad - bias_grad_expected) .< eps)

  for i = 1:n_input
    back_grad = to_array(net.states[2].blobs_diff[i])
    back_grad_expected = W' * top_diffs[i]
    @test all(abs(back_grad - back_grad_expected) .< eps)
  end

  destroy(net)
end

function test_tied_inner_product_layer(backend::Backend)
  test_tied_inner_product_layer(backend, 3, Float32, 1e-3)
  test_tied_inner_product_layer(backend, 3, Float64, 1e-10)
end

if test_cpu
  test_tied_inner_product_layer(backend_cpu)
end
if test_gpu
  test_tied_inner_product_layer(backend_gpu)
end


text:
using Spider

include("plugins/sc_wikilinks.jl")

let
  args = (source_dir="src",
          output_dir="docs",
          clear_odir=false,
          header_file="header.html",
          footer_file="footer.html")

  run_spider(BibTexRefs(),
             ArxivLinks(),
             TOC(),
             SCWikiLinks(;toplevel_prefix="scnumerics"),
             MathJax(),
             BackLinks(),
             ;args...)

  run(`git reset -q`)
  run(`git add docs`)
  run(`git commit --allow-empty -m 'Website generated'`)
end

max_stars_repo_name: emstoudenmire/scnumerics

max_issues_repo_name: emstoudenmire/scnumerics

max_forks_repo_name: emstoudenmire/scnumerics


text:
@testset "mSGNHT" begin
    #-------------------
    # test inplace addition vs reference
    #-------------------
    function sample_reference!(s::mSGNHTState,grad::Function)
        # increment iterations counter
        s.t += 1

        # improve readability
        nparams = length(s.x)
        stepsize = s.stepsize(s.t)
        x = s.x
        p = s.p
        xi = s.xi
        D = s.D
        mu = s.mu

        # second order symmetric splitting scheme from Leimkuhler et al.
        p[:] += .5*stepsize*grad(x)
        x[:] += .5*stepsize*p
        xi[:] +=.5*stepsize/mu* (p.*p - 1)
        ids = abs(stepsize*xi) .> 1e-6
        p[ids] = exp(-stepsize*xi[ids]).*p[ids] + sqrt(.5.*(1-exp(-2.0*stepsize*xi[ids]))./xi[ids] * D) .* randn(sum(ids))
        p[!ids] += sqrt(stepsize*D)*randn(sum(!ids))
        xi[:] +=.5*stepsize/mu* (p.*p - 1)
        x[:] += .5*stepsize*p
        p[:] += .5*stepsize*grad(x)
        s
    end

    srand(1)
    s = mSGNHTState(ones(2),0.01)
    srand(1)
    ref_s = mSGNHTState(ones(2),0.01)
    msgnhtgrad(x) = x.*x

    srand(2)
    sample!(s,msgnhtgrad)

    srand(2)
    sample_reference!(ref_s,msgnhtgrad)

    @test s.x == ref_s.x

    #------------------
    # Gaussian sanity check.
    #------------------
    dm = GaussianMixtureModel([0.,],[1.],[1.])

    gradient = DataModel.getgrad(dm)
    s = mSGNHTState(zeros(1),0.1)

    x_avg = [0.]
    x2_avg = [0.]
    for i in 1:50000
        sample!(s,gradient)
        x_avg =(i-1)/i*x_avg + s.x/i
        x2_avg = (i-1)/i*x2_avg + s.x.^2/i
    end

    @test abs(x_avg-0)[1] < 0.1
    @test abs(x2_avg-1)[1] < 0.2
end


text:
"""
SEPPExpKern(μ, ϕ, γ)

SEPP with exponential kernel model type with rate λ (for all t, there is an event with probability λ(t)) viz. :
```λ(t) = μ + ϕ \\sum_{t_k < t} \\exp(γ(t-t_k))```
where `t_k` are the timestamps of all events.
"""
mutable struct SEPPExpKern <: SEPP
    data::Union{TS, Nothing}
    μ::Real
    ϕ::Real
    γ::Real
    cov_mat::Union{Array{Float64, 2}, Nothing}

    function SEPPExpKern(data::Union{TS, Nothing} = nothing; μ::Real = rand(), ϕ::Real = rand(), γ::Real = rand(), cov_mat::Union{Array{Float64, 2}, Nothing} = nothing)
        if any((μ, ϕ, γ) .< 0)
            error("paramaters must be positive or zero")
        else
            new(data,μ, ϕ, γ, cov_mat)
        end
    end

end

params(sepp::SEPPExpKern)::AbstractDict = Dict(:μ => sepp.μ, :ϕ => sepp.ϕ, :γ => sepp.γ)


text:
using Plots
using UnPack
using LinearAlgebra

using StartUpDG
using StartUpDG.ExplicitTimestepUtils
using EntropyStableEuler.Fluxes1D
import EntropyStableEuler: γ
using FluxDiffUtils

N = 3
K1D = 16
CFL = .5
T = .50

rd = init_reference_interval(N)
VX,EToV = uniform_1D_mesh(K1D)
md = init_DG_mesh(VX,EToV,rd)

@unpack M,Dr,Vq,Pq,Vf,wf,nrJ = rd
Qr = Pq'*M*Dr*Pq
Ef = Vf*Pq
Br = diagm(wf.*nrJ)
Qrh = .5*[Qr-Qr' Ef'*Br;
        -Br*Ef  Br]
Vh = [Vq;Vf]
Ph = M\transpose(Vh)
VhP = Vh*Pq

# make skew symmetric versions of the operators"
Qrhskew = .5*(Qrh-transpose(Qrh))

LX = 2
build_periodic_boundary_maps!(md,rd,LX)

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

@unpack x = md
rho = @. 2 + .5*exp(-25*(x^2))
u = @. 0*x
p = @. rho^γ

Q = primitive_to_conservative(rho,u,p)
SBP_ops = (Matrix(Qrhskew'),Vh,Ph,VhP)

# interpolate geofacs to both vol/surf nodes, pack back into md
@unpack Vq,Vf = rd
@unpack rxJ = md
rxJ = [Vq;Vf]*rxJ # interp to hybridized points
@pack! md = rxJ

# convert to rho,u,v,beta vars
function cons_to_prim_withlogs(rho,rhou,E)
    rho,u,beta = conservative_to_primitive_beta(rho,rhou,E)
    return (rho,u,beta),(log.(rho),log.(beta))
end

mxm_accum!(X,x,e) = X[:,e] .+= 2*Ph*x

function rhs(Q, SBP_ops, rd::RefElemData, md::MeshData)

    QrTr,Vh,Ph,VhP = SBP_ops
    @unpack LIFT,Vf,Vq = rd
    @unpack rxJ,J,K = md
    @unpack nxJ,sJ = md
    @unpack mapP,mapB = md

    Nh,Nq = size(VhP)

    # entropy var projection
    VU = (x->VhP*x).(v_ufun((x->Vq*x).(Q)...))
    Uh = u_vfun(VU...)
    Qh,logQh = cons_to_prim_withlogs(Uh...)
    Qh_log = (Qh...,logQh...)

    # compute face values
    QM = (x->x[Nq+1:Nh,:]).(Qh_log)
    QP = (x->x[mapP]).(QM)

    # simple lax friedrichs dissipation
    Uf =  (x->x[Nq+1:Nh,:]).(Uh)
    (rhoM,rhouM,EM) = Uf
    rhoUM_n = @. rhouM*nxJ/sJ
    lam = abs.(wavespeed_1D(rhoM,rhoUM_n,EM))
    LFc = .25*max.(lam,lam[mapP]).*sJ

    fSx = euler_fluxes(QM...,QP...)
    normal_flux(fx,u) = fx.*nxJ - LFc.*(u[mapP]-u)
    flux = map(normal_flux,fSx,Uf)
    rhsQ = (x->LIFT*x).(flux)

    # compute volume contributions using flux differencing
    rhse = zero.(getindex.(Qh,:,1))
    for e = 1:K
        QxTr = rxJ[1,e]*QrTr

        # computes sum(Ax.*Fx + Ay.*Fy,dims=2)
        # wrap QxTr, euler_fluxes in tuple to make it look "multidimensional"
        hadamard_sum!(rhse,(QxTr,),tuple∘euler_fluxes,getindex.(Qh_log,:,e);
                      skip_index=(i,j)->(i>Nq)&(j>Nq))
        mxm_accum!.(rhsQ,rhse,e)
    end

    rhstest = zero(eltype(first(rhsQ)))
    compute_rhstest = false
    if compute_rhstest
        for fld in eachindex(rhsQ)
            VUq = VU[fld][1:Nq,:]
            rhstest += sum(md.wJq.*VUq.*(Vq*rhsQ[fld]))
        end
    end

    return map(x -> -x./J,rhsQ),rhstest
end

# "Time integration coefficients"
rk4a,rk4b,rk4c = ck45()
CN = (N+1)^2/2  # estimated trace constant
dt = CFL * 2 / (CN*K1D)
Nsteps = convert(Int,ceil(T/dt))
dt = T/Nsteps

# "Perform time-stepping"
resQ = zero.(Q)
for i = 1:Nsteps
    for INTRK = 1:5
        rhsQ,rhstest = rhs(Q,SBP_ops,rd,md)
        bcopy!.(resQ, @. rk4a[INTRK]*resQ + dt*rhsQ)
        bcopy!.(Q,@. Q + rk4b[INTRK]*resQ)
    end

    if i%10==0 || i==Nsteps
        println("Time step $i out of $Nsteps")
    end
end

# plotting nodes
gr(aspect_ratio=1,legend=false)
scatter((x->rd.Vp*x).((x,Q[1])),zcolor=rd.Vp*Q[1],msw=0,ms=2,cam=(0,90))

max_stars_repo_name: yiminllin/ESDG-PosLimit

max_issues_repo_name: yiminllin/ESDG-PosLimit

max_forks_repo_name: yiminllin/ESDG-PosLimit


text:
# function angle2quat(q::Vector{Float64},ang::Vector{Float64},seq::String)
#     c = cos.(ang*0.5);
#     s = sin.(ang*0.5);
#
#     if seq=="321"
#         q[1] = c[1]*c[2]*c[3] + s[1]*s[2]*s[3]; # β0
#         q[2] = c[1]*c[2]*s[3] - s[1]*s[2]*c[3]; # β1
#         q[3] = c[1]*s[2]*c[3] + s[1]*c[2]*s[3]; # β2
#         q[4] = s[1]*c[2]*c[3] - c[1]*s[2]*s[3]; # β3
#     end
# end
#
# function angle2quat(ang::Vector{Float64},seq::String)::Vector{Float64}
#     q = Vector{Float64}(undef,4);
#     angle2quat(q,ang,seq);
#     return(q);
# end
#
# function quat2angle(ang::Vector{Float64},q::Vector{Float64},seq::String)
#     if seq == "321" # Have to program other sequences, but we may not use it.
#         # roll (x-axis rotation)
#         sinr_cosp = +2.0 * (q[1] * q[2] + q[3] * q[4]);
#         cosr_cosp = +1.0 - 2.0 * (q[2]^2 + q[3]^2);
#         ϕ = atan(sinr_cosp, cosr_cosp);
#
#         # pitch (y-axis rotation)
#         sinp = +2.0 * (q[1] * q[3] - q[4] * q[2]);
#         if (abs(sinp) >= 1)
#             θ = copysign(M_PI / 2, sinp); # use 90 degrees if out of range
#         else
#             θ = asin(sinp);
#         end
#
#         # yaw (z-axis rotation)
#         siny_cosp = +2.0 * (q[1] * q[4] + q[2] * q[3]);
#         cosy_cosp = +1.0 - 2.0 * (q[3] * q[3] + q[4] * q[4]);
#         ψ = atan(siny_cosp, cosy_cosp);
#
#         ang[1] = ψ;
#         ang[2] = θ;
#         ang[3] = ϕ;
#     end
# end
#
# function quat2angle(q::Vector{Float64},seq::String)::Vector{Float64}
#     ang = Vector{Float64}(undef,3);
#     quat2angle(ang,q,seq)
#     return(ang);
# end

function dcm2quat(q::Vector{T},C::Matrix{T})::Vector{T} where T<:Real
    # Check dimensions, etc.

    # Apply Sheppard's algorithm
    trace = tr(C);

    v0 = 0.25*(1+trace);
    v1 = 0.25*(1+2*C[1,1]-trace);
    v2 = 0.25*(1+2*C[2,2]-trace);
    v3 = 0.25*(1+2*C[3,3]-trace);

    maxV = max(v0,v1,v2,v3); # Take the maximum.
    if (v0 >= maxV)
        q[1] = sqrt(maxV);
        q[2] = (C[2,3]-C[3,2])/4.0/q[1];
        q[3] = (C[3,1]-C[1,3])/4.0/q[1];
        q[4] = (C[1,2]-C[2,1])/4.0/q[1];
    elseif (v1 >= maxV)
        q[2] = sqrt(maxV);
        q[1] = (C[2,3]-C[3,2])/4.0/q[2];
        q[3] = (C[1,2]+C[2,1])/4.0/q[2];
        q[4] = (C[3,1]+C[1,3])/4.0/q[2];
    elseif (v2 >= maxV)
        q[3] = sqrt(maxV);
        q[1] = (C[3,1]-C[1,3])/4.0/q[3];
        q[2] = (C[1,2]+C[2,1])/4.0/q[3];
        q[4] = (C[2,3]+C[3,2])/4.0/q[3];
    else
        q[4] = sqrt(maxV);
        q[1] = (C[1,2]-C[2,1])/4.0/q[4];
        q[2] = (C[3,1]+C[1,3])/4.0/q[4];
        q[3] = (C[2,3]+C[3,2])/4.0/q[4];
    end
end

function dcm2quat(C::Matrix{T})::Vector{T} where T <: Real
    q = Vector{Real}(undef,4)
    # Apply Sheppard's algorithm
    trace = tr(C);

    v0 = 0.25*(1+trace);
    v1 = 0.25*(1+2*C[1,1]-trace);
    v2 = 0.25*(1+2*C[2,2]-trace);
    v3 = 0.25*(1+2*C[3,3]-trace);

    maxV = max(v0,v1,v2,v3); # Take the maximum.
    if (v0 >= maxV)
        q[1] = sqrt(maxV);
        q[2] = (C[2,3]-C[3,2])/4.0/q[1];
        q[3] = (C[3,1]-C[1,3])/4.0/q[1];
        q[4] = (C[1,2]-C[2,1])/4.0/q[1];
    elseif (v1 >= maxV)
        q[2] = sqrt(maxV);
        q[1] = (C[2,3]-C[3,2])/4.0/q[2];
        q[3] = (C[1,2]+C[2,1])/4.0/q[2];
        q[4] = (C[3,1]+C[1,3])/4.0/q[2];
    elseif (v2 >= maxV)
        q[3] = sqrt(maxV);
        q[1] = (C[3,1]-C[1,3])/4.0/q[3];
        q[2] = (C[1,2]+C[2,1])/4.0/q[3];
        q[4] = (C[2,3]+C[3,2])/4.0/q[3];
    else
        q[4] = sqrt(maxV);
        q[1] = (C[1,2]-C[2,1])/4.0/q[4];
        q[2] = (C[3,1]+C[1,3])/4.0/q[4];
        q[3] = (C[2,3]+C[3,2])/4.0/q[4];
    end
    return q
end

function quat2dcm(C::Matrix{T},q::Vector{T})::Matrix{T} where T <: Real
    β0 = q[1];
    β1 = q[2];
    β2 = q[3];
    β3 = q[4];

    C =  [(β0^2+β1^2-β2^2-β3^2) 2*(β1*β2+β0*β3)           2*(β1*β3-β0*β2);
    2*(β1*β2-β0*β3)         (β0^2-β1^2+β2^2-β3^2)   2*(β2*β3+β0*β1);
    2*(β1*β3+β0*β2)       2*(β2*β3-β0*β1)         (β0^2-β1^2-β2^2+β3^2)];
end

function quat2dcm(q::Vector{T})::Matrix{T} where T <: Real
    # ForwardDiff has problems with inplace function definitions
    β0 = q[1];
    β1 = q[2];
    β2 = q[3];
    β3 = q[4];

    C =  [(β0^2+β1^2-β2^2-β3^2) 2*(β1*β2+β0*β3)           2*(β1*β3-β0*β2);
    2*(β1*β2-β0*β3)         (β0^2-β1^2+β2^2-β3^2)   2*(β2*β3+β0*β1);
    2*(β1*β3+β0*β2)       2*(β2*β3-β0*β1)         (β0^2-β1^2-β2^2+β3^2)];

    # C = Matrix{Float64}(undef,3,3);
    # quat2dcm(C,q);
    return(C);
end

# function ang2dcm(dcm::Matrix{Float64},ang::Vector{Float64},orientation::String)
#
#     θ1 = ang[1]; # phi
#     θ2 = ang[2]; # theta
#     θ3 = ang[3]; # psi
#
#     c1 = cos(θ1); s1 = sin(θ1);
#     c2 = cos(θ2); s2 = sin(θ2);
#     c3 = cos(θ3); s3 = sin(θ3);
#
#     M1 = [1 0 0;
#           0 c1 s1;
#           0 -s1 c1];
#
#     M2 = [c2 0 -s2;
#            0 1 0;
#           s2 0 c2];
#
#     M3 = [c3 s3 0;
#          -s3 c3 0;
#            0  1 1];
#
#     if orientation == "321"
#         C = M1*M2*M3; # Inertial to body
#     else
#         error("Uknown orientation specified")
#     end
# end

function skewX(x::Vector{T})::Matrix{T} where T<: Real
    X = [    0 -x[3]  x[2]
          x[3]     0 -x[1]
         -x[2]  x[1]    0]
    return X
end

function angVel(β::Vector{T},βdot::Vector{T})::Vector{T} where T<: Real
    # β0 = β[1];
    # β1 = β[2];
    # β2 = β[3];
    # β3 = β[4];

    # E1 = [-β1  β0  β3 -β2
    #       -β2 -β3  β0  β1
    #       -β3  β2 -β1  β0]
    E1 = genE(β)[2:4,:]
    ω = 2*E1*βdot
    return ω
end

function genE(β::Vector{T})::Matrix{T} where T<:Real
    X = [                        transpose(β)
         -β[2:4] β[1]*Matrix{Float64}(I,3,3) - skewX(β[2:4])]
    return X

end

function genEdot(βdot::Vector{T})::Matrix{T} where T<:Real
    X = [                        transpose(βdot)
         -βdot[2:4] βdot[1]*Matrix{Float64}(I,3,3) - skewX(βdot[2:4])]
    return X

end

function quaternionProduct(β1::Vector{T},β2::Vector{T})::Vector{T} where T <: Real
    x = [β1[1]*β2[1] - dot(β1[2:4],β2[2:4]);
         β1[2]*β2[2:4] + β2[1]*β1[2:4] + cross(β1[2:4],β2[2:4])]
    return x
end

function quatNormErr(β::Vector{T}) where T<:Real
    return (norm(β)-1)
end

function randomQuat()
    phi = pi/2*rand();
    e = 1/sqrt(3);
    q = [cos(phi/2);e*sin(phi/2)*ones(3)];
    # u,v,w = rand(3)
    # q = [sqrt(1-u)*sin(2*π*v);
    #      sqrt(1-u)*cos(2*π*v);
    #      sqrt(u)*sin(2*π*w);
    #      sqrt(u)*cos(2*π*w)]
    return q
end

max_stars_repo_name: isrlab/Multibody-Dynamics

max_issues_repo_name: isrlab/Multibody-Dynamics

max_forks_repo_name: isrlab/Multibody-Dynamics


text:
using Test, PolynomialBases

ufunc(x) = sinpi(x)^2

function tolerance(p, T=Float64)
    if p <= 2
        1.
    elseif p <= 4
        0.78
    elseif p <= 6
        0.02
    elseif p <= 8
        1.e-4
    elseif p <= 10
        if T==Float32
            4.e-7
        else
            1.5e-7
        end
    elseif p <= 12
        if T==Float32
            3.e-7
        else
            9.e-11
        end
    elseif p <= 14
        if T==Float32
            4.e-7
        else
            3.e-14
        end
    else
        if T==Float32
            5.e-7
        else
            5.e-16
        end
    end
end

for basis_type in subtypes(PolynomialBases.NodalBasis{PolynomialBases.Line})
    basis_type <: Union{ClosedNewtonCotes,ClenshawCurtis} && continue
    for p in 3:15, T in (Float32, Float64)
        basis = basis_type(p, T)
        u = ufunc.(basis.nodes)
        @test abs(integrate(u, basis) - 1) < tolerance(p, T)
        @test abs(integrate(u, basis.weights) - 1) < tolerance(p, T)
    end
end


text:
# Socket options

for (fset, fget, k, T) in [
    (:set_affinity,                :get_affinity,                 4, UInt64)
    (:set_type,                    :get_type,                    16,   Cint)
    (:set_linger,                  :get_linger,                  17,   Cint)
    (:set_reconnect_ivl,           :get_reconnect_ivl,           18,   Cint)
    (:set_backlog,                 :get_backlog,                 19,   Cint)
    (:set_reconnect_ivl_max,       :get_reconnect_ivl_max,       21,   Cint)
    (:set_rate,                    :get_rate,                     8,   Cint)
    (:set_recovery_ivl,            :get_recovery_ivl,             9,   Cint)
    (:set_sndbuf,                  :get_sndbuf,                  11,   Cint)
    (:set_rcvbuf,                  :get_rcvbuf,                  12,   Cint)
    (nothing,                      :_zmq_getsockopt_rcvmore,     13,   Cint)
    (nothing,                      :get_events,                  15,   Cint)
    (:set_maxmsgsize,              :get_maxmsgsize,              22,   Cint)
    (:set_sndhwm,                  :get_sndhwm,                  23,   Cint)
    (:set_rcvhwm,                  :get_rcvhwm,                  24,   Cint)
    (:set_multicast_hops,          :get_multicast_hops,          25,   Cint)
    (:set_ipv4only,                :get_ipv4only,                31,   Cint)
    (:set_tcp_keepalive,           :get_tcp_keepalive,           34,   Cint)
    (:set_tcp_keepalive_idle,      :get_tcp_keepalive_idle,      35,   Cint)
    (:set_tcp_keepalive_cnt,       :get_tcp_keepalive_cnt,       36,   Cint)
    (:set_tcp_keepalive_intvl,     :get_tcp_keepalive_intvl,     37,   Cint)
    (:set_rcvtimeo,                :get_rcvtimeo,                27,   Cint)
    (:set_sndtimeo,                :get_sndtimeo,                28,   Cint)
    (:set_immediate,               :get_immediate,               39,   Cint)
    (nothing,                      :get_fd,                      14, Sys.iswindows() ? Ptr{Cvoid} : Cint)
    ]
    if fset != nothing
        @eval function $(Symbol("_",fset))(socket::Socket, option_val::Integer)
            rc = ccall((:zmq_setsockopt, libzmq), Cint,
                       (Ptr{Cvoid}, Cint, Ref{$T}, Csize_t),
                       socket, $k, option_val, sizeof($T))
            if rc != 0
                throw(StateError(jl_zmq_error_str()))
            end
        end
        prop = QuoteNode(Symbol(String(fset)[5:end]))
        @eval @deprecate $fset(socket::Socket, option_val::Integer) setproperty!(socket, $prop, option_val)
    end
    if fget != nothing
        @eval function $(Symbol("_",fget))(socket::Socket)
            val = Ref{$T}()
            rc = ccall((:zmq_getsockopt, libzmq), Cint,
                       (Ptr{Cvoid}, Cint, Ref{$T}, Ref{Csize_t}),
                       socket, $k, val, sizeof($T))
            if rc != 0
                throw(StateError(jl_zmq_error_str()))
            end
            return Int(val[])
        end
        prop = QuoteNode(Symbol(String(fget)[5:end]))
        @eval @deprecate $fget(socket::Socket) getproperty(socket, $prop)
    end
end

# For some functions, the publicly-visible versions should require &
# return boolean:
_get_rcvmore(socket::Socket) = Bool(__zmq_getsockopt_rcvmore(socket))
@deprecate get_rcvmore(socket::Socket) getproperty(socket, :rcvmore)
@deprecate ismore(socket::Socket) getproperty(socket, :rcvmore)

# subscribe/unsubscribe options take an arbitrary byte array
for (f,k) in ((:subscribe,6), (:unsubscribe,7))
    f_ = Symbol(f, "_")
    @eval begin
        function $f_(socket::Socket, filter::Ptr{T}, len::Integer) where {T}
            rc = ccall((:zmq_setsockopt, libzmq), Cint,
                       (Ptr{Cvoid}, Cint, Ptr{T}, Csize_t),
                       socket, $k, filter, len)
            if rc != 0
                throw(StateError(jl_zmq_error_str()))
            end
        end
        $f(socket::Socket, filter::Union{Vector{UInt8},String}) =
            @preserve filter $f_(socket, pointer(filter), sizeof(filter))
        $f(socket::Socket, filter::AbstractString) = $f(socket, String(filter))
        $f(socket::Socket) = $f_(socket, C_NULL, 0)
    end
end

# string properties
for (fset, fget, k) in [
    (:_set_routing_id,           :_get_routing_id,              5)
    (nothing,                    :_get_last_endpoint,          32)
    # (:_set_tcp_accept_filter,       nothing,                     38) #  deprecated
    ]
    if fset != nothing
        @eval function ($fset)(socket::Socket, option_val::String)
            if sizeof(option_val) > 255
                throw(StateError("option value too large"))
            end
            rc = ccall((:zmq_setsockopt, libzmq), Cint,
                       (Ptr{Cvoid}, Cint, Ptr{UInt8}, Csize_t),
                       socket, $k, option_val, sizeof(option_val))
            if rc != 0
                throw(StateError(jl_zmq_error_str()))
            end
        end
        @eval ($fset)(socket::Socket, option_val::AbstractString) = $fset(socket, String(option_val))
    end
    if fget != nothing
        @eval function ($fget)(socket::Socket)
            buf = Base.StringVector(255)
            len = Ref{Csize_t}(sizeof(buf))
            rc = ccall((:zmq_getsockopt, libzmq), Cint,
                       (Ptr{Cvoid}, Cint, Ptr{UInt8}, Ref{Csize_t}),
                       socket, $k, buf, len)
            if rc != 0
                throw(StateError(jl_zmq_error_str()))
            end
            return String(resize!(buf, len[]))
        end
    end
end
@deprecate get_identity(socket::Socket) getproperty(socket, :routing_id)
@deprecate set_identity(socket::Socket, val) setproperty!(socket, :routing_id, val)

# getproperty/setproperty API for socket properties
const sockprops = (:affinity, :type, :linger, :reconnect_ivl, :backlog, :reconnect_ivl_max,
                   :rate, :recovery_ivl, :sndbuf, :rcvbuf, :rcvmore, :events, :maxmsgsize,
                   :sndhwm, :rcvhwm, :multicast_hops, :ipv4only,
                   :tcp_keepalive, :tcp_keepalive_idle, :tcp_keepalive_cnt, :tcp_keepalive_intvl,
                   :rcvtimeo, :sndtimeo, :fd, :routing_id, :last_endpoint, :immediate)

Base.propertynames(::Socket) = sockprops
@eval function Base.getproperty(value::Socket, name::Symbol)
    $(propexpression(filter!(p -> isdefined(@__MODULE__, Symbol("_get_", p)), collect(sockprops))) do p
        :($(Symbol("_get_", p))(value))
    end)
end
@eval function Base.setproperty!(value::Socket, name::Symbol, x)
    $(propexpression(filter!(p -> isdefined(@__MODULE__, Symbol("_set_", p)), collect(sockprops))) do p
        :($(Symbol("_set_", p))(value, x))
    end)
    return x
end


text:
struct Halley <: SRSolver end
struct SuperHalley <: SRSolver end
struct InverseQuadratic <: SRSolver end
struct CMethod <: SRSolver end

# 1 SuperHalley = 2 InverseQuadratic = 3 CMethod = 4
function impliedVolatilitySRHalley(
    isCall::Bool,
    price::T,
    f::T,
    strike::T,
    tte::T,
    df::T,
    ftolrel::T,
    maxEval::Int,
    solver::SRSolver,
)::T where {T}
    c, ex = normalizePrice(isCall, price, f, strike, df)
        #println(" c ",c, " 1/ex ",1/ex)
    if (c >= min(1, 1 / ex)) || (c <= 0)
        throw(DomainError(c, string("Price out of range, must be < ", min(1, 1 / ex), " and > 0")))
    end
    x = log(ex)
    guess = T(impliedVolatilitySRGuessUndiscountedCall(c, ex, x))
    b = guess
    xtolrel = 32 * eps(T) * max(1, guess)
    xtol = 0 * eps(T)
    ftolrel = max(ftolrel, eps(T))
    ftol_ = max(1, abs(c)) * ftolrel
    useLog = true
    local fb, fbOverfpb, fp2bOverfpb
    if useLog
        c = log(c)
        fb, fbOverfpb, fp2bOverfpb = objectiveLog(x, ex, b, c)
    else
        fb, fbOverfpb, fp2bOverfpb = objective(x, ex, b, c)
    end

    if abs(fb) < ftol_
        return b / sqrt(tte)
    end
    for iteration = 1:maxEval
        x0 = b
        lf = fbOverfpb * fp2bOverfpb
        x1 = x0 + srSolve(solver, fb, fbOverfpb, lf)
        a = x0
        if useLog
            fb, fbOverfpb, fp2bOverfpb = objectiveLog(x, ex, x1, c)
        else
            fb, fbOverfpb, fp2bOverfpb = objective(x, ex, b, c)
        end
        # println(iteration, " ", fb, " ", fbOverfpb, " ", fp2bOverfpb, " ", x1, " ", x0, " ", lf)
        b = x1
        xtol_ = xtol + max(1, abs(b)) * xtolrel

        if abs(b - a) <= xtol_ || abs(fb) <= ftol_
            #value = b
            break
        end
    end
    return b / sqrt(tte)
end

function srSolve(solver::Halley, fb::T, fbOverfpb::T, lf::T)::T where {T}
    return -(1 / (1 - lf / 2)) * fbOverfpb
end
function srSolve(solver::SuperHalley, fb::T, fbOverfpb::T, lf::T)::T where {T}
    return -(1 + lf / (2 * (1 - lf))) * fbOverfpb #super halley faster for med low vols/more OTM, almost similar for ATM, but unstable for very small fb
end
function srSolve(solver::InverseQuadratic, fb::T, fbOverfpb::T, lf::T)::T where {T}
    return -(1 + lf / 2) * fbOverfpb #inverse quadratic
end
function srSolve(solver::CMethod, fb::T, fbOverfpb::T, lf::T)::T where {T}
    return -(1 + lf / 2 + lf * lf / 2) * fbOverfpb  #C method seems more balanced perf
end

function objectiveLog(x::T, ex::T, v::T, logc::T) where {T}
    v = abs(v)
    h = x / v
    t = v / 2
    # naive implementation:
    # Np = normcdf(h+t)
    # Nm = normcdf(h-t)/ex
    # norm = T(1)
    # cEstimate = norm * (Np - Nm)
    # vega = normpdf(h+t)
    # volgaOverVega = (h + t) * (h - t) / v
    # logcEstimate = log(cEstimate)
    # logvega = vega/cEstimate #u'/u
    #logVolga = (volgaOverVega*logvega - logvega*logvega)#u'' u - u'2 /u^2

    # the following actually helps convergence and accuracy (up to a certain machine epsilon point):
    sqrt2 = sqrt(T(2))
    Np = erfcx(-(h + t) / sqrt2)
    Nm = erfcx(-(h - t) / sqrt2)
    eh2t2 = exp(-(h * h + t * t) / 2)
    norm = 1 / (2 * sqrt(ex)) * eh2t2
    cEstimate = norm * (Np - Nm)
    logcEstimate = log(cEstimate)
    twopi = 2 * T(pi)
    logvega = (2 /sqrt(twopi)) / (Np - Nm)  #u'/u
    volgaOverVega = (h + t) * (h - t) / v
    logvolgaOverVega = volgaOverVega - logvega
    return logcEstimate - logc, (logcEstimate - logc) / logvega, logvolgaOverVega
end
function objective(x::T, ex::T, v::T, c::T) where {T}
    v = abs(v)
    h = x / v
    t = v / 2
    # naive implementation:
    # Np = normcdf(h+t)
    # Nm = normcdf(h-t)/ex
    # norm = T(1)
    # cEstimate = norm * (Np - Nm)
    # vega = normpdf(h+t)
    # cEstimateOverVega = (cEstimate-c)/vega
    # volgaOverVega = (h + t) * (h - t) / v
    # the following actually helps convergence and accuracy (up to a certain machine epsilon point):
    sqrt2 = sqrt(T(2))
    Np = erfcx(-(h + t) / sqrt2)
    Nm = erfcx(-(h - t) / sqrt2)
    eh2t2 = exp(-(h * h + t * t) / 2)
    norm = 1 / (2 * sqrt(ex)) * eh2t2
    cEstimate = norm * (Np - Nm)
    eht = exp(-h * t)
    twopi = 2 * T(pi)
    vega = eh2t2 * eht / sqrt(twopi)    #normpdf(h + t)  #this is exp(-0.5* (h^2+t^2 +2*h*t)) = eh2t2*exp(-h*t)
    cEstimateOverVega = (Np - Nm) / (2 * eht * sqrt(ex / (twopi))) - c / vega
    volgaOverVega = (h + t) * (h - t) / v
    return cEstimate - c, cEstimateOverVega, volgaOverVega
end


text:
Base.@kwdef mutable struct information
    modelName::String = ""
    deviceSerialNumber::String = ""
    deviceFirmwareVersion::String = ""
end

Base.@kwdef mutable struct settings
    acquisitionMode::Symbol = :continuous #[:continuous,:singleFrame,:multiFrame]
    acquisitionFramerate::AbstractFloat = 30.0
    exposureMode::Symbol = :timed #[:timed,:triggerWidth]
    exposureAuto::Symbol = :off #[:off,:once,:continuous]
    exposureTime::AbstractFloat = 100.00
    gainAuto::Symbol = :off #[:off,:once,:continuous]
    gain::AbstractFloat = 0.0
    gamma::AbstractFloat = 1.25
    blackLevel::AbstractFloat = 5.0
    deviceLinkThroughputLimit::Int64 = 383328000

    width::Int64 = 100
    height::Int64 = 30
    offsetX::Int64 = 0
    offsetY::Int64 = 0
    pixelFormat::String = ""
    binningHorizontal::Int64 = 1
    binningVertical::Int64 = 1
end

Base.@kwdef mutable struct settingsLimits
    acquisitionMode::Vector{Symbol} = [:continuous,:singleFrame,:multiFrame]
    acquisitionFramerate::Tuple{AbstractFloat,AbstractFloat} = (0.0,60.0)
    exposureMode::Vector{Symbol} = [:timed,:triggerWidth]
    exposureAuto::Vector{Symbol} = [:off,:once,:continuous]
    exposureTime::Tuple{AbstractFloat,AbstractFloat} = (0.0,100.0)
    autoExposureTime::Tuple{AbstractFloat,AbstractFloat} = (0.0,100.0)
    gainAuto::Vector{Symbol} = [:off,:once,:continuous]
    gain::Tuple{AbstractFloat,AbstractFloat} = (0.0,10.0)
    gamma::Tuple{AbstractFloat,AbstractFloat} = (0.0,10.0)
    blackLevel::Tuple{AbstractFloat,AbstractFloat} = (0.0,4000.0)
    deviceLinkThroughputLimit::Tuple{Int64,Int64} = (0,383328000)

    width::Tuple{Int64,Int64} = (0,3000)
    height::Tuple{Int64,Int64} = (0,3000)
    offsetX::Tuple{Int64,Int64} = (0,3000)
    offsetY::Tuple{Int64,Int64} = (0,3000)
    binningHorizontal::Tuple{Int64,Int64} = (1,2)
    binningVertical::Tuple{Int64,Int64} = (1,2)
end

Base.@kwdef mutable struct GPIO
    triggerSelector::Symbol = :frameStart #[:frameStart,:exposureActive]
    triggerMode::Symbol = :off #[:off,:on]
    triggerSource::Symbol = :line0 #[:software,:line0,:line3,:line3]
    triggerActivation::Symbol = :fallingEdge #[:risingEdge,:fallingEdge]
    triggerOverlap::Symbol = :off #[:off,:readOut]
    triggerDelay::AbstractFloat = 0.0
    lineSelector::Symbol = :line0 #[:line0,:line1,:line2,:line3]
    lineMode::Symbol = :input #[:input]
    lineInverter::Bool = true
    lineSource::Symbol = :nothing #[:exposureActive,:externalTriggerActive,:userOutput1]
    userOutputSelector::Symbol = :UserOutputValue1 #[:UserOutputValue1,:UserOutputValue2,:UserOutputValue3]
    userOutputValue::Bool = true
end

Base.@kwdef mutable struct GPIOLimits
    triggerSelector::Vector{Symbol} = [:frameStart,:exposureActive]
    triggerMode::Vector{Symbol} = [:off,:on]
    triggerSource::Vector{Symbol} = [:software,:line0,:line3,:line3]
    triggerActivation::Vector{Symbol} = [:risingEdge,:fallingEdge]
    triggerOverlap::Vector{Symbol} = [:off,:readOut]
    triggerDelay::Tuple{AbstractFloat,AbstractFloat} = (0.0,typemax(1.0))
    lineSelector::Vector{Symbol} = [:line0,:line1,:line2,:line3]
    lineMode::Vector{Symbol} = [:input]
    lineSource::Vector{Symbol} = [:exposureActive,:externalTriggerActive,:userOutput1]
    userOutputSelector::Vector{Symbol} = [:UserOutputValue1,:UserOutputValue2,:UserOutputValue3]
end

function camSettingsRead!(cam,camSettings)
    camSettings.acquisitionFramerate = framerate(cam)
    ex,mode = exposure(cam)
    camSettings.exposureTime = ex/1000
    camSettings.gain,mode = gain(cam)
    camSettings.width,camSettings.height = imagedims(cam)
    camSettings.offsetX,camSettings.offsetY = offsetdims(cam)
end

function camSettingsLimitsRead!(cam,camSettingsLimits)
    camSettingsLimits.acquisitionFramerate = framerate_limits(cam)
    camSettingsLimits.exposureTime = exposure_limits(cam)./1000
    camSettingsLimits.autoExposureTime = autoexposure_limits(cam)./1000

    camSettingsLimits.gain = gain_limits(cam)

    # Image size
    camSettingsLimits.width,camSettingsLimits.height = imagedims_limits(cam)
    camSettingsLimits.offsetX,camSettingsLimits.offsetY = offsetdims_limits(cam)
end

function camSettingsUpdater(;timerInterval::AbstractFloat=1/10)
    global sessionStat
    global camSettings, camSettingsLimits
    global camGPIO, camGPIOLimits

    updateTimer = Timer(0,interval=timerInterval)
    lastSessionStat = deepcopy(sessionStat)
    lastCamSettings = deepcopy(camSettings)
    lastCamSettingsLimits = deepcopy(camSettingsLimits)
    lastCamGPIO = deepcopy(camGPIO)
    first_autoexposure = true
    while !sessionStat.terminate
        t_before = time()
        if isrunning(cam)

            # RECORDING

            # FRAMERATE
            if (camSettings.acquisitionFramerate != lastCamSettings.acquisitionFramerate)
                framerate!(cam,camSettings.acquisitionFramerate)
                lastCamSettings.acquisitionFramerate = camSettings.acquisitionFramerate
                camSettingsLimitsRead!(cam,camSettingsLimits)
                camSettingsRead!(cam,camSettings)
                if camSettings.exposureAuto != :off
                    camSettingsLimits.autoExposureTime = autoexposure_limits!(cam,(0,1e9)) #intentionally too small/large to find lims
                end
            end

            # EXPOSURE
            if (camSettings.exposureAuto != lastCamSettings.exposureAuto) || (camSettings.exposureTime != lastCamSettings.exposureTime)
                if camSettings.exposureAuto == :off
                    exposure!(cam,camSettings.exposureTime*1000)
                elseif camSettings.exposureAuto == :continuous
                    exposure!(cam)
                end
                lastCamSettings.exposureAuto = camSettings.exposureAuto
                lastCamSettings.exposureTime = camSettings.exposureTime
            end
            if camSettings.exposureAuto != :off
                ex,mode = exposure(cam)
                camSettings.exposureTime = ex/1000
                if first_autoexposure
                    #Catches first opportunity to set auto exposure lims, in case framerate hasn't changed first
                    #(can only be set when auto exposure is enabled)
                    camSettingsLimits.autoExposureTime = autoexposure_limits!(cam,(0,1e9)) #intentionally too small/large to find lims
                    first_autoexposure = false
                end
            end

            # GAIN
            if (camSettings.gainAuto != lastCamSettings.gainAuto) || (camSettings.gain != lastCamSettings.gain)
                if camSettings.gainAuto == :off
                    gain!(cam,camSettings.gain)
                elseif camSettings.gainAuto == :continuous
                    gain!(cam)
                end
                lastCamSettings.gainAuto = camSettings.gainAuto
                lastCamSettings.gain = camSettings.gain
            end
            if camSettings.gainAuto != :off
                camSettings.gain, _ = gain(cam)
            end

            # IMAGE OFFSET
            if (camSettings.offsetX != lastCamSettings.offsetX) || (camSettings.offsetY != lastCamSettings.offsetY)
                try
                    offsetdims!(cam,(camSettings.offsetX,camSettings.offsetY))
                    lastCamSettings.offsetX = camSettings.offsetX
                    lastCamSettings.offsetY = camSettings.offsetY
                    camSettingsLimitsRead!(cam,camSettingsLimits)
                catch e
                    #@info "Selected x or y offset not allowed"
                    camSettings.offsetX = lastCamSettings.offsetX
                    camSettings.offsetY = lastCamSettings.offsetY
                end
            end
        else
            # IMAGE SIZE (Can't be set while running)
            if (camSettings.width != lastCamSettings.width) || (camSettings.height != lastCamSettings.height)
                try
                    imagedims!(cam,(camSettings.width,camSettings.height))
                    lastCamSettings.width = camSettings.width
                    lastCamSettings.height = camSettings.height
                    camSettingsLimitsRead!(cam,camSettingsLimits)
                    sessionStat.resolutionupdate = true
                catch e
                    #@info "Selected width or height not allowed"
                    camSettings.width = lastCamSettings.width
                    camSettings.height = lastCamSettings.height
                end
            end
        end
        if time()-t_before < timerInterval
            wait(updateTimer)
        else
            yield()
        end

    end
    close(updateTimer)
    updateTimer = nothing
end

function camGPIOLimitsUpdater!(camGPIOLimits::GPIOLimits)

end


text:
using Oceananigans, Oceananigans.Units, GLMakie

grid = RectilinearGrid(CPU(), size = (64, 64, 64), halo = (3, 3, 3),
                       x = (0, 128), y = (0, 128), z = (-64, 0),
                       topology = (Periodic, Periodic, Bounded))

# Specify coriolis, viscosity / turbulence closure, and boundary condition
coriolis = FPlane(f₀=1e-4)
closure = AnisotropicMinimumDissipation()
boundary_conditions = (; b = FieldBoundaryConditions(top=FluxBoundaryCondition(1e-8)))

#=
parameters = (μ₀ = 1/day,  # surface growth rate
              λ = 5,       # sunlight attenuation length scale (m)
              m = 0.1/day) # mortality rate due to virus and zooplankton grazing
@inline growing_and_grazing(x, y, z, t, P, params) = (params.μ₀ * exp(z / params.λ) - params.m) * P
plankton_dynamics = Forcing(growing_and_grazing; field_dependencies = :P, parameters)
tracers = (:b, :P)
forcing = (; P = plankton_dynamics)
=#
tracers = :b

model = NonhydrostaticModel(; grid, coriolis, closure, boundary_conditions, tracers,
                              timestepper = :RungeKutta3, advection = WENO5(),
                              buoyancy = BuoyancyTracer())

bᵢ(x, y, z) = 1e-6 * z
ϵ(x, y, z) = 1e-6 * randn() # noise
set!(model, b=bᵢ, u=ϵ, v=ϵ, w=ϵ)

# Time-stepping
simulation = Simulation(model, Δt=10minutes, stop_time=12hours)

progress(sim) = @printf("Iter: %d, time: %s, Δt: %s, max|w|: %s\n",
                        iteration(sim), prettytime(sim), prettytime(sim.Δt), maximum(abs, model.velocities.w))

wizard = TimeStepWizard(cfl=0.8, max_change=1.1, max_Δt=2minutes)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10))

# Turbulent kinetic energy
u, v, w = model.velocities
e = Field((u^2 + v^2 + w^2) / 2)
ξ = Field(∂x(w) + ∂z(u))

simulation.output_writers[:fields] = JLD2OutputWriter(model, (; e, w, ξ, b = model.tracers.b),
                                                      schedule = TimeInterval(1hour),
                                                      prefix = "ocean_convection",
                                                      field_slicer = nothing,
                                                      force = true)

run!(simulation)

# Visualization

data_path = simulation.output_writers[:fields].filepath
e = FieldTimeSeries(data_path, "e")
b = FieldTimeSeries(data_path, "b")
x, y, z = nodes(e)
Nt = length(e.times)

fig = Figure(resolution=(1600, 800))
axe = Axis(fig[1, 1])
axξ = Axis(fig[1, 2])
slider = Slider(fig[3, :], range=1:Nt, startvalue=1)
n = slider.value

eⁿ = @lift dropdims(interior(e[$n]), dims=2)
bⁿ = @lift dropdims(interior(b[$n]), dims=2)
hm = heatmap!(axe, x, z, eⁿ; colormap=:solar, colorrange=(0, maximum(e) / 4))
Colorbar(fig[2, 1], hm, vertical=false, flipaxis=false, label="Turbulent kinetic energy (m² s⁻²)")

hm = heatmap!(axξ, x, z, bⁿ; colormap=:deep, colorrange=(minimum(b), maximum(b)))
Colorbar(fig[2, 2], hm, vertical=false, flipaxis=false, label="Buoyancy (m s⁻²)")

title = @lift "Idealized ocean convection at t = " * prettytime(e.times[$n])
Label(fig[0, :], title)

display(fig)

#record(fig, "idealized_ocean_convection.mp4", 1:Nt, framerate=12) do nn
#    @info "Rendering frame $nn of $Nt...
#    n[] = nn
#end


text:
#-------------------------------------------------------------------------------
# bearing and range available

"""
    $TYPEDEF

Bearing and Range constraint from a Pose2 to Point2 variable.
"""
mutable struct Pose2Point2BearingRange{B <: IIF.SamplableBelief, R <: IIF.SamplableBelief} <: IIF.AbstractManifoldMinimize
  bearing::B
  range::R
end

getManifold(::IIF.InstanceType{<:Pose2Point2BearingRange}) = ProductGroup(ProductManifold(SpecialOrthogonal(2), TranslationGroup(1)))

function getSample(cfo::CalcFactor{<:Pose2Point2BearingRange})
  # defaults, TODO better reuse
  M = getManifold(cfo.factor)
  e0 = ProductRepr([1 0; 0 1.], [0.])

  # vector of tangents
  smpl = hat(M, e0, [rand(cfo.factor.bearing), rand(cfo.factor.range)])

  # return IIF `::Tuple` format
  return smpl
end

function IIF.getMeasurementParametric(s::Pose2Point2BearingRange{<:Normal, <:Normal})

  meas = [mean(s.bearing), mean(s.range)]
  iΣ = [1/var(s.bearing)             0;
                      0  1/var(s.range)]

  return meas, iΣ
end

function (cfo::CalcFactor{<:Pose2Point2BearingRange})(meas, xi, lm)
  SE2 = SpecialEuclidean(2)
  #FIXME fix this factor to work correctly on manifolds
  if false

  Mt = SE2.manifold.manifolds[1]
  Mr = SE2.manifold.manifolds[2]

  # compose rotations
  rRi = xi.parts[2]
  rRo = retract(Mr, rRi, meas.parts[1])

  # new SE2 objects containing the rotation
  rTo = ProductRepr(xi.parts[1], rRo)
  oTl = ProductRepr([meas.parts[2];0], identity_element(Mr, rRi))

  # prediction landmark
  rTo = Manifolds.compose(SE2, rTo, oTl)

  # cartesian difference in predicted and estimated landmark
  δ_l = Manifolds.compose(Mt, inv(Mt, rTo.parts[1]), lm)

  # find the residuals as though in polar coordinates
  δθ = atan(δ_l[2], δ_l[1])
  δr = norm(δ_l)

  return [δθ; δr]

  # FIXME: this is close to the old BR that worked, use to ensure tests are working
  else
    Xi = vee(SE2, xi, log(SE2, identity_element(SE2, xi), xi))

    Mmeas = getManifold(cfo.factor)
    Xmeas = vee(Mmeas, Manifolds.Identity(Mmeas), meas)
    # 1-bearing
    # 2-range
    # world frame
    meas_θ = Xmeas[1]
    meas_r = Xmeas[2]

    θ = meas_θ + Xi[3]
    mx = meas_r*cos(θ)
    my = meas_r*sin(θ)

    ex = lm[1] - (mx + Xi[1])
    ey = lm[2] - (my + Xi[2])

    # res = [eθ, er]
    eθ = atan((my + Xi[2]), (mx + Xi[1])) - atan(lm[2], lm[1])  # eθ
    er= sqrt(ex^2 + ey^2) # some wasted computation here       # er

    # return hat(facM, ProductRepr([1. 0;0 1], [0.]), [eθ, er])
    return [eθ, er]
  end
end
# quick check
# pose = (0,0,0),  bear = 0.0,  range = 10.0   ==>  lm = (10,0)
# pose = (0,0,0),  bear = pi/2,  range = 10.0   ==>  lm = (0,10)
# pose = (0,0,pi/2),  bear = 0.0,  range = 10.0   ==>  lm = (0,10)
# pose = (0,0,pi/2),  bear = pi/2,  range = 10.0   ==>  lm = (-10,0)

# function (cfo::CalcFactor{<:Pose2Point2BearingRange})(measX, p, q)
#   #
#   M = SpecialEuclidean(2)
#   q_SE = ProductRepr(q, identity_element(SpecialOrthogonal(2), p.parts[2]))

#   X_se2 = log(M, identity_element(M, p), compose(M, inv(M, p), q_SE))
#   X = X_se2.parts[1]
#   # NOTE this: `X̂ = log(M, p, q_SE)` wrong for what we want
#   #go to tangent vector
#   return measX - X
# end

# Support for database based solving

passTypeThrough(d::FunctionNodeData{Pose2Point2Range}) = d

mutable struct PackedPose2Point2BearingRange <: IncrementalInference.PackedInferenceType
    bearstr::String
    rangstr::String
end

function convert( ::Type{PackedPose2Point2BearingRange}, d::Pose2Point2BearingRange )
  return PackedPose2Point2BearingRange( convert(PackedSamplableBelief, d.bearing), convert(PackedSamplableBelief, d.range) )
end

# TODO -- should not be resorting to string, consider specialized code for parametric distribution types and KDEs
function convert( ::Type{Pose2Point2BearingRange}, d::PackedPose2Point2BearingRange )
  Pose2Point2BearingRange( convert(SamplableBelief, d.bearstr), convert(SamplableBelief, d.rangstr) )
end


text:
module SpringCollab2020FlagToggleStarRotateSpinner

using ..Ahorn, Maple

@pardef FlagToggleStarRotateSpinner(x1::Integer, y1::Integer, x2::Integer=x1 + 16, y2::Integer=y1, clockwise::Bool=false, flag::String="flag_toggle_star_rotate_spinner", inverted::Bool=false) =
    Entity("SpringCollab2020/FlagToggleStarRotateSpinner", x=x1, y=y1, nodes=Tuple{Int, Int}[(x2, y2)], clockwise=clockwise, flag=flag, inverted=inverted)

function rotatingSpinnerFinalizer(entity::FlagToggleStarRotateSpinner)
    x, y = Int(entity.data["x"]), Int(entity.data["y"])
    entity.data["x"], entity.data["y"] = x + 32, y
    entity.data["nodes"] = [(x, y)]
end

const placements = Ahorn.PlacementDict(
    "Star (Rotating, Clockwise, Flag Toggle) (Spring Collab 2020)" => Ahorn.EntityPlacement(
        FlagToggleStarRotateSpinner,
        "point",
        Dict{String, Any}(
            "clockwise" => true
        ),
        rotatingSpinnerFinalizer
    ),
    "Star (Rotating, Flag Toggle) (Spring Collab 2020)" => Ahorn.EntityPlacement(
        FlagToggleStarRotateSpinner,
        "point",
        Dict{String, Any}(
            "clockwise" => false
        ),
        rotatingSpinnerFinalizer
    ),
)

Ahorn.nodeLimits(entity::FlagToggleStarRotateSpinner) = 1, 1

function Ahorn.selection(entity::FlagToggleStarRotateSpinner)
    nx, ny = Int.(entity.data["nodes"][1])
    x, y = Ahorn.position(entity)

    return [Ahorn.Rectangle(x - 8, y - 8, 16, 16), Ahorn.Rectangle(nx - 8, ny - 8, 16, 16)]
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::FlagToggleStarRotateSpinner, room::Maple.Room)
    clockwise = get(entity.data, "clockwise", false)
    dir = clockwise ? 1 : -1

    centerX, centerY = Int.(entity.data["nodes"][1])
    x, y = Ahorn.position(entity)

    radius = sqrt((centerX - x)^2 + (centerY - y)^2)

    Ahorn.drawCircle(ctx, centerX, centerY, radius, Ahorn.colors.selection_selected_fc)
    Ahorn.drawArrow(ctx, centerX + radius, centerY, centerX + radius, centerY + 0.001 * dir, Ahorn.colors.selection_selected_fc, headLength=6)
    Ahorn.drawArrow(ctx, centerX - radius, centerY, centerX - radius, centerY + 0.001 * -dir, Ahorn.colors.selection_selected_fc, headLength=6)

    Ahorn.drawSprite(ctx, "danger/starfish13", x, y)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::FlagToggleStarRotateSpinner, room::Maple.Room)
    centerX, centerY = Int.(entity.data["nodes"][1])

    Ahorn.drawSprite(ctx, "danger/starfish13", centerX, centerY)
end

end


text:
⊙(s::AbstractString,::Nothing) = s
⊙(::Nothing,::Nothing) = nothing
⊙(n::Nothing,s::AbstractString) = s⊙n
⊙(l::AbstractString,r::AbstractString) = l*r

const ROMAN_NUMERALS = merge!(Base.ImmutableDict{Int,String}(),Dict(
    1 => "I",
    2 => "II",
    3 => "III",
    4 => "IV",
    5 => "V",
    6 => "VI",
    7 => "VII",
    8 => "VIII",
    9 => "IX",
    10 => "X",
    20 => "XX",
    30 => "XXX",
    40 => "XL",
    50 => "L",
    60 => "LX",
    70 => "LXX",
    80 => "LXXX",
    90 => "XC",
    100 => "C",
    200 => "CC",
    300 => "CCC",
    400 => "CD",
    500 => "D",
    600 => "DC",
    700 => "DCC",
    800 => "DCCC",
    900 => "CM",
    1000 => "M",
    2000 => "MM",
    3000 => "MMM",
))

# Using immutable dicts is a little bit of a faff, but worth it for a 100x speedup
function roman(digit)
    # The pattern is "add them up until you are 1 step away from the next digit, then do subtraction"
    # I don't fancy coding that up
    digit == 0 && return nothing
    ROMAN_NUMERALS[digit]
end

function to_roman(number::Integer)
    (number <= 0  || number >= 3999) && throw(ErrorException("That number doesn't exist :)"))
    mapreduce(t->roman(t[2]*(10^(t[1]-1))),⊙,Iterators.reverse(enumerate(digits(number,base=10))))
end

# Hat-tip to jlapeyre: https://discourse.julialang.org/t/way-to-construct-multi-key-immutabledict/10294/2#post_5
function Base.merge!(d::Base.ImmutableDict, dsrc::Dict)
    for (k,v) in dsrc
        d = Base.ImmutableDict(d,k=>v)
    end
    d
end


text:
import ..SamplingConstraints:
    StartToEnd,
    StrictlyDecreasing,
    strictly_increasing_sequence_exists,
    constrain,
    TruncateMaximum,
    TruncateRange,
    TruncateQuantiles

import IntervalArithmetic:
    interval

"""
    resample(udata::DT, sequential_constraint::StrictlyDecreasing{T};
        quantiles = [0.0001, 0.9999]) where {DT <: AbstractUncertainValueDataset, T <: StartToEnd}

Element-wise resample the uncertain values in the dataset such that each preceding value
is strictly larger in magnitude than the next one.

## Arguments:

- **`udata`**: An uncertain dataset.
- **`sequential_constraint`**: An instance of a `StrictlyDecreasing` sequential
    sampling constraint.
- **`ordered_sampling_alg`**: An instance of a `StartToEnd` ordered
    sampling constraint, indicating that the sequence of decreasing values should
    be created in one go from start to finish (as opposed to chunking the data set
    first, then gluing partial sequences together afterwards).
- **`quantiles`**: A two-element vector representing a quantile range which is used
    to truncate the supports values in the dataset before drawing the sequence of
    decreasing values. This deals with distributions with infinite support.
"""
function resample(udata::DT,
    sequential_constraint::StrictlyDecreasing{T};
    quantiles = [0.0001, 0.9999]) where {DT <: AbstractUncertainValueDataset, T <: StartToEnd}

    n_vals = length(udata)

    # If a strictly increasing sequence exists through the reversed dataset,
    # then a strictly decreasing sequence also exists.
    if !strictly_increasing_sequence_exists(udata[end:-1:1], quantiles)
        throw(ArgumentError("No strictly decreasing sequence through dataset exists."))
    end

    # Some values may have infinite support even after applying the first round of
    # constraints, so truncate to the provided (large) quantile range to ensure that all
    # values have finite support when sampling.
    c = TruncateQuantiles(quantiles...,)
    constrained_data = [udata[i] isa CertainValue ? udata[i] : constrain(udata[i], c) for i = 1:length(udata)]

    # Pre-allocate a sample
    sample = Vector{Float64}(undef, n_vals)

    # Sample the first value in a way that ensures a strictly
    # increasing sequence from indices 2:end will exist.
    hi = maximum(constrained_data[1])
    lo = minimum(maximum.(constrained_data[1:end]))

    if lo == hi
        sample[1] = resample(constrained_data[1])
    else
        sample[1] = resample(constrained_data[1], TruncateRange(lo, hi))
    end

    for i = 2:n_vals - 1
        # Find the lower and upper bounds of the support from which we can draw values
        # while still ensuring an increasing sequence of values.
        hi = sample[i - 1]
        lo = maximum([minimum(constrained_data[i]);
                        minimum.(constrained_data[(i + 1):end])])

        # Constrain the support of the furnishing distribution according
        # to the bounds and sample a value from it.
        cd = constrained_data[i]

        if cd isa CertainValue # no need to constrain if value is not uncertain
            sample[i] = resample(cd)
        else
            constrained_val = constrain(cd, TruncateRange(lo, hi))
            sample[i] = resample(constrained_val)
        end
    end

    # The last value doesn't need an upper bound, because there's no
    # value following it we need to consider.
    sample[n_vals] = resample(constrained_data[n_vals], TruncateMaximum(sample[end - 1]))

    return sample
end

function resample(udata::DT, constraint::SCT, sequential_constraint::StrictlyDecreasing{T};
        quantiles = [0.0001, 0.9999]) where {SCT <: SamplingConstraint,
            DT <: AbstractUncertainValueDataset,
            T <: StartToEnd}

    resample(constrain(udata, constraint), sequential_constraint, quantiles = quantiles)
end

function resample(udata::DT, constraint::Vector{SCT}, sequential_constraint::StrictlyDecreasing{T};
        quantiles = [0.0001, 0.9999]) where {SCT <: SamplingConstraint,
                    DT <: AbstractUncertainValueDataset,
                    T <: StartToEnd}

    resample(constrain(udata, constraint), sequential_constraint, quantiles = quantiles)
end

"""
    resample(udata::AbstractUncertainValueDataset,
        constraint::Union{SamplingConstraint, Vector{SamplingConstraint}},
        sequential_constraint::StrictlyDecreasing{OrderedSamplingAlgorithm};
        quantiles = [0.0001, 0.9999])

Draw a sequence of values strictly decreasing in magnitude from the dataset, sampling
each of the furnishing distributions once each, after first truncating the supports
of the values in the dataset using the provided `constraint`s.

## Arguments:

- **`udata`**: An uncertain dataset.
- **`constraint`**: Sampling constraint(s) to apply to each of the values in the dataset
before drawing the sequence of decreasing values.
- **`sequential_constraint`**: An instance of a `StrictlyDecreasing` sequential
sampling constraint. For example, `StrictlyDecreasing(StartToEnd())` indicates
that a strictly decreasing sequence should be created in one go from start to
finish (as opposed to chunking the data set first, then gluing partial sequences
together afterwards).
- **`quantiles`**: A two-element vector representing a quantile range which is used
to truncate the supports values in the dataset before drawing the sequence of
decreasing values. This deals with distributions with infinite support.
"""
resample(udata::AbstractUncertainValueDataset,
    constraint::Union{SamplingConstraint, Vector{SamplingConstraint}},
    sequential_constraint::StrictlyDecreasing{OrderedSamplingAlgorithm};
    quantiles = [0.0001, 0.9999])

"""
function resample(udata::UncertainIndexValueDataset,
    sequential_constraint::StrictlyDecreasing,
    quantiles = [0.0001, 0.9999])

Resample an uncertain index-value dataset by enforcing strictly decreasing indices.
"""
function resample(udata::UncertainIndexValueDataset,
        sequential_constraint::StrictlyDecreasing{<:OrderedSamplingAlgorithm})
    resample(udata.indices, sequential_constraint), resample(udata.values)
end

function resample(udata::UncertainIndexValueDataset,
        constraint::SamplingConstraint,
        sequential_constraint::StrictlyDecreasing{<:OrderedSamplingAlgorithm})

    inds = resample(constrain(udata.indices, constraint), sequential_constraint)
    vals = resample(constrain(udata.indices, constraint))

    inds, vals
end

function resample(udata::UncertainIndexValueDataset,
        constraint::Vector{SamplingConstraint},
        sequential_constraint::StrictlyDecreasing{<:OrderedSamplingAlgorithm})

    inds = resample(constrain(udata.indices, constraint), sequential_constraint)
    vals = resample(constrain(udata.indices, constraint))

    inds, vals
end

function resample(udata::UncertainIndexValueDataset,
        idx_constraint::SamplingConstraint,
        value_constraint::SamplingConstraint,
        sequential_constraint::StrictlyDecreasing{<:OrderedSamplingAlgorithm})

    inds = resample(constrain(udata.indices, idx_constraint), sequential_constraint)
    vals = resample(constrain(udata.indices, value_constraint))

    inds, vals
end

function resample(udata::UncertainIndexValueDataset,
        idx_constraint::Vector{<:SamplingConstraint},
        value_constraint::SamplingConstraint,
        sequential_constraint::StrictlyDecreasing{<:OrderedSamplingAlgorithm})

    inds = resample(constrain(udata.indices, idx_constraint), sequential_constraint)
    vals = resample(constrain(udata.indices, value_constraint))

    inds, vals
end

function resample(udata::UncertainIndexValueDataset,
        idx_constraint::SamplingConstraint,
        value_constraint::Vector{<:SamplingConstraint},
        sequential_constraint::StrictlyDecreasing{<:OrderedSamplingAlgorithm})

    inds = resample(constrain(udata.indices, idx_constraint), sequential_constraint)
    vals = resample(constrain(udata.indices, value_constraint))

    inds, vals
end

function resample(udata::UncertainIndexValueDataset,
        idx_constraint::Vector{<:SamplingConstraint},
        value_constraint::Vector{<:SamplingConstraint},
        sequential_constraint::StrictlyDecreasing{<:OrderedSamplingAlgorithm})

    inds = resample(constrain(udata.indices, idx_constraint), sequential_constraint)
    vals = resample(constrain(udata.indices, value_constraint))

    inds, vals
end

export resample


text:
# Example 1 : 	Oscillating masses

# Constraints :
# 				Position in box
# 				Velocity in box
#				Acceleration in box
#       		Acceleration rate in box

# problem definition
module prb

	using LinearAlgebra, StaticArrays

	const nx = 8											# state dim
	const nu = 4											# imput dim
	const N = 30											# horizon length
	const Δ = 0.25 											# sampling time
	const tvec = [(j-1)*Δ for j in 1:N]   					# time vector (t=0:N-1)

	# solution method settings
		# common
		const kmax = 500										# max iterations for main algorithms (and outer loop)
		const scale_Hg = true

		# JuMP
		const ϵ_pd = 1e-9										# solver precision
		const ϵ_gap = 1e-9										# solver precision

		# PIPG
		const ϵ_pipg = 1e-4										# termination precision
		const restrt_idx = 40									# index for restarting

		# ADMM
		const α_admm = 1.0										# step size
		const ϵ_nest_admm = 1e-10								# Nesterov's method termination precision
		const kmax_nest = 100									# max iterations for Nesterov's method

		# Chambolle and Pock (PDHG)
		const β⁰ = 2.0											# initial step size
		const γ⁰ = 1.0											# initial step size

	# system parameters
	const nxb2 = Int64(nx/2)									# position veloity dim
	const pmax = 2.0											# position
	const vmax = 2.0											# velocity
	const umax = 2.0					 						# input (forcing)
	const dumax = 0.5											# input rate

	# reference trajectory
	const xref = [SVector{nx}(cat(ones(nxb2)...,zeros(nxb2)...,dims=1)) for _ in 1:N-1]	# state reference (t=1:N-1)
	const uref = [SVector{nu}(zeros(nu)) for _ in 1:N-1]							# input reference (t=0:N-2)
	const x0 = 1 .* ones(nx)								    					# initial condition (t=0)

	# dynamics (t=0:N-2)
	const Ac = Array([zeros(nxb2,nxb2) Diagonal(ones(nxb2))
						       Tridiagonal(ones(nxb2-1),-2 .* ones(nxb2),ones(nxb2-1)) zeros(nxb2,nxb2)])
	const Ad = SMatrix{nx,nx}(exp(Ac*Δ))
	const Bc = [zeros(nxb2,nxb2)
				Diagonal(ones(nxb2,nxb2))]
	const Bd = SMatrix{nx,nu}(inv(Ac)*(Ad-Diagonal(ones(nx)))*Bc)

	const A = [Ad for _ in 1:N-1]
	const B = [Bd for _ in 1:N-1]
	const e = [SVector{nx}(zeros(nx)) for _ in 1:N-1]

	# state and control penalty (LQR weights)

	const q_tmp = [1.0,1.0]
	const r_tmp = 1.0

	# const q_tmp = [0.5,1.0]
	# const r_tmp = 1.0

	const Q = [SMatrix{nx,nx}(Array(Diagonal(cat(q_tmp[1] .* ones(nxb2), q_tmp[2] .* ones(nxb2), dims=1)))) for _ in 1:N-1] # t=1:N-1
	const R = [SMatrix{nu,nu}(r_tmp .* Array(Diagonal(ones(nu)))) for _ in 1:N-1] 											# t=0:N-2

end

include("proj_funcs.jl") 		# projection library

# assemble quantities for vectorized PIPG+ implementation
module asm

	using ..prb
	using ..proj 								# projection onto simple sets
	using LinearAlgebra, StaticArrays, JuMP

	const scale_Hg = prb.scale_Hg
	const nH = 2*prb.nu*(prb.N-2)				# no. of inequality constraints

	function construct_H(nx,nu,N,A,B,Inx,Inu)
		nxnu = nx+nu
		H = zeros(nx*(N-1) + nH, nxnu*(N-1))
		H[1:nx,1:nxnu] .= [-B[1] Inx]
		for j=1:N-2
			H[j*nx+1:(j+1)*nx, nu+(j-1)*nxnu+1:nu+(j-1)*nxnu+2*nx+nu] .= [-A[j] -B[j] Inx]

			H[nx*(N-1)+(j-1)*nu+1:nx*(N-1)+j*nu, (j-1)*nxnu+1:(j-1)*nxnu+nu+nxnu] .= [Inu zeros(nu,nx) -Inu]

			H[nx*(N-1)+nu*(N-2)+(j-1)*nu+1:nx*(N-1)+nu*(N-2)+j*nu, (j-1)*nxnu+1:(j-1)*nxnu+nu+nxnu] .= [-Inu zeros(nu,nx) Inu]
		end

		norm_row_H = zeros(nx*(N-1)+nH)
		for j=1:nx*(N-1)+nH
			norm_row_H[j] = norm(H[j,1:end],Inf)
		end

		if scale_Hg
			# print("Scaling rows of H\n")
			for j=1:length(norm_row_H)
				H[j,1:end] .= H[j,1:end] ./ norm_row_H[j]
			end
		end

		return H,norm_row_H
	end

	function construct_G(nx,nu,N,Q,R)
		nxnu = nx+nu
		G = zeros(nxnu*(N-1),nxnu*(N-1))
		for j=1:N-1
			G[(j-1)*nxnu+1:j*nxnu,(j-1)*nxnu+1:j*nxnu] .= [R[j] 		zeros(nu,nx)
														   zeros(nx,nu) Q[j]]
		end
		return G
	end

	H_scl,norm_row_H = construct_H(prb.nx,prb.nu,prb.N,prb.A,prb.B,Diagonal(ones(prb.nx)),Diagonal(ones(prb.nu)))
	const H = deepcopy(H_scl)
	const HT = transpose(H)

	const zref = cat([cat(prb.uref[t],prb.xref[t],dims=1) for t=1:prb.N-1]...,dims=1)
	const h = cat([cat(-prb.R[t]*prb.uref[t], -prb.Q[t]*prb.xref[t], dims=1) for t=1:prb.N-1]...,dims=1)

	if scale_Hg
		# print("Scaling components of g\n")
		const scl_val_Hg = norm_row_H
	else
		const scl_val_Hg = 1.0
	end

	const g = cat(prb.A[1]*prb.x0+prb.e[1],cat([prb.e[t] for t=2:prb.N-1]...,dims=1), -prb.dumax .* ones(2*prb.nu*(prb.N-2)),dims=1)

	@assert length(h) == (prb.nx+prb.nu)*(prb.N-1) "Length of h is incorrect."
	@assert length(g) == 2*prb.nu*(prb.N-2)+prb.nx*(prb.N-1) "Length of g is incorrect."

	const Heq = [H -1 .* Array(Diagonal(ones(prb.nx*(prb.N-1) + nH)))]
	const HeqT = transpose(Heq)

	const G = construct_G(prb.nx,prb.nu,prb.N,prb.Q,prb.R)

	const svd_HTH = svd(HT*H)
	const svd_HeqTHeq = svd(HeqT*Heq)

	const svd_G = svd(G)

	const σ = max(svd_HTH.S...)
	const σ_eq = max(svd_HeqTHeq.S...)
	const λ = max(svd_G.S...)
	const μ = min(svd_G.S...)

	function construct_xu(z::Array{Float64,1})
		return [z[nu+(j-1)*nxnu+1:nu+(j-1)*nxnu+nx] for j=1:N-1], [z[(j-1)*nxnu+1:(j-1)*nxnu+nu] for j=1:N-1]
	end

	# useful constants
	const nx = copy(prb.nx)
	const nu = copy(prb.nu)
	const nxnu = nx+nu
	const nxb2 = Int64(nx/2)
	const N = copy(prb.N)

	const κ_x = MVector{nx}(randn(nx))
	const κ_u = MVector{nu}(randn(nu))

	# projection
	function project_U!(u)
		κ_u .= u
		proj.box!(u,κ_u,-prb.umax,prb.umax,nu)
	end
	function project_X!(x)
		κ_x .= x
		proj.box!(x[1:nxb2],κ_x[1:nxb2],-prb.pmax,prb.pmax,nxb2)
		proj.box!(x[nxb2+1:nx],κ_x[nxb2+1:nx],-prb.vmax,prb.vmax,nxb2)
	end

	# set state and input constraints in JuMP (vectorized version)
	function set_constraint_jump_vectorized!(model,zz)
		# acceleration, position and velocity in box
		@constraint(model,[j=1:N-1],-prb.umax .≤ zz[(j-1)*nxnu+1:(j-1)*nxnu+nu] .≤  prb.umax)
		@constraint(model,[j=1:N-1],-prb.pmax .≤ zz[nu+(j-1)*nxnu+1:nu+(j-1)*nxnu+nxb2] .≤  prb.pmax)
		@constraint(model,[j=1:N-1],-prb.vmax .≤ zz[nu+(j-1)*nxnu+nxb2+1:nu+(j-1)*nxnu+nx] .≤  prb.vmax)
	end

	# set state and input constraints (everything other than dynamics and initial condition) in JuMP (decomposed version)
	function set_constraint_jump_decomposed!(model,x,u)
		@constraint(model,[t=1:N-1], -prb.pmax .≤ x[1:nxb2,t+1] .≤ prb.pmax)
		@constraint(model,[t=1:N-1], -prb.vmax .≤ x[nxb2+1:nx,t+1] .≤ prb.vmax)
		@constraint(model,[t=1:N-1], -prb.umax .≤ u[1:nu,t] .≤ prb.umax)
		@constraint(model,[t=1:N-2],-prb.dumax .≤ u[:,t+1] .- u[:,t] .≤ prb.dumax)
	end

end

max_stars_repo_name: purnanandelango/pipg-demo

max_issues_repo_name: UW-ACL/pipg-demo

max_forks_repo_name: UW-ACL/pipg-demo


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

mutable struct SnapshotSku <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    tier::Any # spec type: Union{ Nothing, String } # spec name: tier

    function SnapshotSku(;name=nothing, tier=nothing)
        o = new()
        validate_property(SnapshotSku, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(SnapshotSku, Symbol("tier"), tier)
        setfield!(o, Symbol("tier"), tier)
        o
    end
end # type SnapshotSku

const _property_map_SnapshotSku = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("tier")=>Symbol("tier"))
const _property_types_SnapshotSku = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("tier")=>"String")
Base.propertynames(::Type{ SnapshotSku }) = collect(keys(_property_map_SnapshotSku))
Swagger.property_type(::Type{ SnapshotSku }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_SnapshotSku[name]))}
Swagger.field_name(::Type{ SnapshotSku }, property_name::Symbol) =  _property_map_SnapshotSku[property_name]

const _allowed_SnapshotSku_name = ["Standard_LRS", "Premium_LRS", "Standard_ZRS"]

function check_required(o::SnapshotSku)
    true
end

function validate_property(::Type{ SnapshotSku }, name::Symbol, val)
    if name === Symbol("name")
        Swagger.validate_param(name, "SnapshotSku", :enum, val, _allowed_SnapshotSku_name)
    end
end


text:
## basic
# envs
abstract type AbstractEnv end

# registered envs
mutable struct RegisteredEnvs
    __envs::Array{AbstractEnv, 1}
    __xs::Array{Any, 1}
end
__REGISTERED_ENVS = RegisteredEnvs([], [])

# generic functions
function readable end
function raw end

## zoo


text:
const TPixel = RGB{N0f8}            # Pixel type

mutable struct mjSim
    modelfile

    # visual interaction controls
    lastx::Float64
    lasty::Float64
    button_left::Bool
    button_middle::Bool
    button_right::Bool

    lastbutton::GLFW.MouseButton
    lastclicktm::Float64
    lastcmdkey::Union{GLFW.Key, Nothing}

    refreshrate::Int

    # function keys
    showhelp::Int
    showoption::Bool
    showinfo::Bool
    showdepth::Bool
    showfullscreen::Bool
    # stereo::Bool
    showsensor::Bool
    # profiler::Bool

    slowmotion::Bool
    paused::Bool
    keyreset::Int

    # Video recording
    record::Any
    vidbuf::Vector{UInt8}
    imgstack::Array{Array{TPixel,2},1}

    framecount::Float64
    #framenum::Int
    #lastframenum::Int

    # MuJoCo things
    scn::Ref{mjvScene}
    cam::Ref{mjvCamera}
    vopt::Ref{mjvOption}
    pert::Ref{mjvPerturb}
    con::Ref{mjrContext}
    figsensor::Ref{mjvFigure}
    m::jlModel
    d::jlData

    # Robot controller
    robot::Union{Robot, Nothing}

    # GLFW handle
    window::GLFW.Window
    vmode::GLFW.VidMode

    # Widgets
    # Inset camera view
    showcam::Bool
    camid::Int
    camviewport::mjrRect
    camviewid::Int
    cambuf::Vector{UInt8}

    #uistate::mjuiState
    #ui0::mjUI
    #ui1::mjUI

    function mjSim(m::jlModel, d::jlData, name::String; width=0, height=0)
        vmode = GLFW.GetVideoMode(GLFW.GetPrimaryMonitor())
        @info("monitor resolution: $(vmode.width)x$(vmode.height)")
        w = width > 0 ? width : Int(floor(vmode.width / 2))
        h = height > 0 ? height : Int(floor(vmode.height / 2))

        s = 512 # size of square camera viewport in pixels

        # grab the camera viewport image from the center of the window
        camviewport = mjrRect((w-s)/2, (h-s)/2, s, s)

        new("", 0.0, 0.0, false, false, false, GLFW.MOUSE_BUTTON_1, 0.0, nothing,
            vmode.refreshrate,
            0, false, false, false, false, false, false, true, 0,
            nothing,
            Vector{UInt8}(undef, 0),
            [],
            0.0, #0, 0,
            Ref(mjvScene()),
            Ref(mjvCamera()),
            Ref(mjvOption()),
            Ref(mjvPerturb()),
            Ref(mjrContext()),
            Ref(mjvFigure()),
            m, d,
            nothing,
            GLFW.CreateWindow(w, h, name),
            vmode,
            false, -1, camviewport, -1, Vector{UInt8}(undef, 0) # Inset camera view
            #ui1
        )
    end
end

#export mjSim


text:
# Hessian Eigenmaps (HLLE)
# ---------------------------
# Hessian eigenmaps: Locally linear embedding techniques for high-dimensional data,
# D. Donoho and C. Grimes, Proc Natl Acad Sci U S A. 2003 May 13; 100(10): 5591–5596

import Combinatorics: combinations

#### HLLE type
struct HLLE{T <: Real} <: AbstractDimensionalityReduction
    k::Int
    λ::AbstractVector{T}
    proj::Projection{T}

    HLLE{T}(k::Int, λ::AbstractVector{T}, proj::Projection{T}) where T = new(k, λ, proj)
end

## properties
outdim(R::HLLE) = size(R.proj, 1)
eigvals(R::HLLE) = R.λ
neighbors(R::HLLE) = R.k

## show
summary(io::IO, R::HLLE) = print(io, "Hessian Eigenmaps(outdim = $(outdim(R)), neighbors = $(neighbors(R)))")

## interface functions
function fit(::Type{HLLE}, X::AbstractMatrix{T}; maxoutdim::Int=2, k::Int=12) where {T<:Real}
    n = size(X, 2)

    # Identify neighbors
    D, E = find_nn(X, k)

    # Obtain tangent coordinates and develop Hessian estimator
    hs = (maxoutdim*(maxoutdim+1)) >> 1
    W = spzeros(T, hs*n, n)
    for i=1:n
        II = @view E[:,i]
        # re-center points in neighborhood
        μ = mean(X[:, II], dims=2)
        N = X[:, II] .- μ
        # calculate tangent coordinates
        tc = svd(N).V[:,1:maxoutdim]

        # Develop Hessian estimator
        Yi = [ones(T, k) tc zeros(T, k, hs)]
        for ii=1:maxoutdim
            Yi[:,maxoutdim+ii+1] = tc[:,ii].^2
        end
        yi = 2*(1+maxoutdim)
        for (ii,jj) in combinations(1:maxoutdim, 2)
            Yi[:, yi] = tc[:, ii] .* tc[:, jj]
            yi += 1
        end
        F = qr(Yi)
        H = transpose(F.Q[:,(end-(hs-1)):end])
        W[(1:hs).+(i-1)*hs, II] = H
    end

    # decomposition
    λ, V = decompose(transpose(W)*W, maxoutdim)
    return HLLE{T}(k, λ, transpose(V) .* convert(T, sqrt(n)))
end

transform(R::HLLE) = R.proj


text:
config = MOI.Test.Config()
optimize!_inf(mock) = MOIU.mock_optimize!(mock, MOI.INFEASIBLE)
optimize!(mock) = MOIU.mock_optimize!(mock,
    [1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0,
     0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
     1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0,
     0.0, 0.0, 0.0, 0.0, -0.0, 0.0,
     1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0,
     0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
     1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0,
     0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
     1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0,
     0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0,
     4.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.0, 4.0, 0.0, 0.0,
     0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
     0.0, -0.0, 0.0, -0.0, 0.0, 0.0, -0.0, 0.0])
@testset "sos $(typeof(mock))" for mock in mocks(optimize!_inf, optimize!_inf, optimize!)
    Tests.sos_horn_test(mock, config)
end
@testset "dsos $(typeof(mock))" for mock in mocks(optimize!_inf)
    Tests.dsos_horn_test(mock, config)
end
@testset "sdsos $(typeof(mock))" for mock in mocks(optimize!_inf)
    Tests.sdsos_horn_test(mock, config)
end

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
using Test, RanDAGen
const dir = joinpath(dirname(pathof(RanDAGen)), "..", "test")

@testset "RanDAGen" begin
    for f in ["DAG.jl",]
        file = joinpath(dir, f)
        println("Running $file tests...")
        if isfile(file)
            include(file)
        else
            @show readdir(dirname(file))
        end
    end
end


text:
abstract type GameObject end

function edges(this::GameObject)
    pos = this.pos
    if pos[1] > this.r + @width
        pos[1] = -this.r
    elseif this.pos[1] < -this.r
        pos[1] = this.r + @width
    end

    if pos[2] > this.r + @height
        pos[2] = -this.r
    elseif this.pos[2] < -this.r
        pos[2] = this.r + @height
    end
end

max_stars_repo_name: jorge-brito/LynxGallery

max_issues_repo_name: jorge-brito/LynxGallery

max_forks_repo_name: jorge-brito/LynxGallery


text:
function Base.show(io::IO, mime::MIME"text/plain", ruleset::Ruleset)
    indent = "  "
    ctx = IOContext(io, :compact => true, :indent => indent)
    printstyled(io, Base.nameof(typeof(ruleset)), "(\n"; color=:blue)
    println(io, indent * "rules = (")
    rulectx = IOContext(ctx, :indent => indent * indent)
    for rule in rules(ruleset)
        show(rulectx, mime, rule)
        print(io, ",\n")
    end
    print(io, indent * "),\n")
    _showsettings(ctx, mime, settings(ruleset))
    print(io, ")\n\n")
    ModelParameters.printparams(io, ruleset)
end

function Base.show(io::IO, ::MIME"text/plain", rule::T) where T<:Rule{R,W} where {R,W}
    indent = get(io, :indent, "")
    if R === :_default_ && W === :_default_
        printstyled(io, indent, Base.nameof(typeof(rule)); color=:red)
    else
        printstyled(io, indent, Base.nameof(typeof(rule)),
                    "{", sprint(show, R), ",", sprint(show, W), "}"; color=:red)
    end
    print(io, "(")
    if !get(io, :compact, false)
        if nfields(rule) > 0
            println(io)
            for fn in fieldnames(T)
                if fieldtype(T, fn) <: Union{Number,Symbol,String}
                    println(io, indent, "    ", fn, " = ", repr(getfield(rule, fn)))
                else
                    # Avoid prining arrays etc. Just show the type.
                    println(io, indent, "    ", fn, " = ", fieldtype(T, fn))
                end
            end
        end
    end
    print(io, ")")
end

function Base.show(io::IO, mime::MIME"text/plain", s::SimSettings)
    string =
"""
SimSettings(;
$(_showsettings(io, mime, s)))
"""
    print(io, string)
end

function _showsettings(io::IO, ::MIME"text/plain", s::SimSettings)
    indent = get(io, :indent, "")
    settings = """
    $(indent)boundary = $(s.boundary),
    $(indent)proc = $(s.proc),
    $(indent)opt = $(s.opt),
    $(indent)cellsize = $(s.cellsize),
    $(indent)timestep = $(s.timestep),
    """
    print(io, settings)
end

function Base.show(io::IO, mime::MIME"text/plain", chain::Chain{R,W}) where {R,W}
    indent = get(io, :indent, "")
    ctx = IOContext(io, :compact => true, :indent => indent * "  ")
    printstyled(io, indent, string("Chain{", sprint(show, R), ",", sprint(show, W), "}"); color=:green)
    print(io, "(")
    for rule in rules(chain)
        println(io)
        show(ctx, mime, rule)
        print(io, ",")
    end
    print(io, "\n)")
end


text:
"""
	Nihilist

Performs encryption, decryption, and chart generation based on the VIC Soviet cipher;
a straddling checkerboard of the nihilist family of ciphers.
"""

module Nihilist

using DataFrames
using CSV

#Type union for interpreting loading a chart from file in the encode/decode stages
DataFrameOrPath=Union{DataFrame, String}

"""
	chart_maker(; file_loc::String="", file_out::String="")

Creates a chart.
EX:
┌───┬┬───┬───┬───┬───┬───┬───┬───┬───┬───┬───┐
│   ││ 3 │ 4 │ 9 │ 0 │ 7 │ 6 │ 8 │ 2 │ 5 │ 1 │
├───┼┼───┼───┼───┼───┼───┼───┼───┼───┼───┼───┤
│   ││   │ i │ r │ e │ a │ t │ s │ n │ o │   │
├───┼┼───┼───┼───┼───┼───┼───┼───┼───┼───┼───┤
│ 1 ││ v │ k │ y │ j │ d │ m │ z │ g │ c │ w │
├───┼┼───┼───┼───┼───┼───┼───┼───┼───┼───┼───┤
│ 3 ││ u │ b │ l │ q │ f │ h │ p │ . │ / │ x │
└───┴┴───┴───┴───┴───┴───┴───┴───┴───┴───┴───┘

The most commonly used letters are placed at the top row, and encoded as a single digit. Two gaps are left
for the second and third rows to double encode leading with the row digit. The codes
are placed inline as the output text, and commonly added against a key (modulo 10) of
some psuedo-random data, pre-selected and shared.

##Parameters:
* `file_loc` - location to import chart
* `file_out` - location to output chart

##Returns:
DataFrame of a chart
"""

function chart_maker(; file_loc::String="", file_out::String="")
	if length(file_loc)>0 #load from file
		try
			return CSV.read(file_loc)
		catch
			error("File could not be read:\n$file_loc")
		end

	else #make a chart
		#hardcode the numbers just to grab them randomly easier,
		#and hardcode the alphabets rather then rulefind for "aeinorst"
		number_pool="0123456789"
		top_letters="aeinorst" #most commonly used letters; placed at top to encode as single digits
		bottom_letters="bcdfghjklmpquvwxyz./" #/ is special character to denote numbers
		row_numbers=[" "] #The top row has no 'row' (single digit vs double digit encoding)
		col_numbers=[]

		#build column by filling randomly 0-9, and removing that choice from list (avoid duplicates)
		while length(number_pool)>0
			rand_spot=1+(rand(UInt8)%length(number_pool))
			push!(col_numbers, Symbol(string(number_pool[rand_spot])))
			number_pool=join(split(number_pool, number_pool[rand_spot]))
		end
		#set row numbers- any two unequal numbers 0-9
		while length(row_numbers)<3
			digit=string(rand(UInt8)%10)
			if !in(digit, row_numbers)
				push!(row_numbers, digit)
			end
		end

		#get avoid POS's - find position in col_nums where row_nums are, and use as mask
		avoid=[]
		for i in row_numbers[2:end]
			push!(avoid, findfirst(col_numbers, Symbol(i)))
		end

		#Get alphabet array ready here; set row/col indexers, and fill alphabet with ""
		row_pos=1; col_pos=1
		chart_alph=fill(" ", length(row_numbers), length(col_numbers))
		#just randomly 'pull' letters from list, just like we did the numbers, but dont pull if its an avoid
		#Note on 'avoid': The cipherchart double encodes letters using a leading row digit. To avoid
		#confusion with the single encoded top row when decoding, those two digits are never used alone.
		while length(top_letters)>0
			if !in(col_pos, avoid)
				rand_spot=1+(rand(UInt8)%length(top_letters))
				chart_alph[row_pos]=string(top_letters[rand_spot])
				top_letters=join(split(top_letters, top_letters[rand_spot]))
			end
			row_pos+=length(row_numbers); col_pos+=1
		end
		#row_pos is shifted for the second row
		row_pos=row_pos%length(chart_alph)+1
		while length(bottom_letters)>0
			rand_spot=1+(rand(UInt8)%length(bottom_letters))
			chart_alph[row_pos]=string(bottom_letters[rand_spot])
			bottom_letters=join(split(bottom_letters, bottom_letters[rand_spot]))
			row_pos+=length(row_numbers)
			if row_pos>length(chart_alph) #hit limit, move to 'next' row
				row_pos=row_pos%length(chart_alph)+1
			end
		end

		#Pass alphabet into DataFrame, and rename according to the col/rows generated.
		new_chart=DataFrame(chart_alph)
		rename!(new_chart, f=>t for (f,t)=zip(names(new_chart), col_numbers))
		new_chart[:CODE]=row_numbers

		#Save to file
		if length(file_out)>0
			try
				#Write a df to a file
				CSV.write(file_out, new_chart)
			catch
				error("File could not be read:\n$file_out")
			end
		end

		return new_chart
	end
end

"""
	encode(message::String, chart::DataFrameOrPath; key::Integer=0, spacing_enable::Bool=true, group_length::Integer=5,
			group_per_line::Integer=8, key_pos::Integer=1, num_repeat::Integer=3, file_out::String="", padding::Bool=true)

Handles encoding of a plaintext message into code according to the provided chart Dataframe. If a key is provided,
it will be added to the message modulo 10.

##Parameters:
* `message` - The message to be encoded - pass as a filepath to load from a file
* `chart` - The chart Dataframe to use in encoding
* `key` - If a key is passed, added to the message %10
* `spacing_enable` - flag to enable group spacing
* `group_length` - length of each group of codes
* `group_per_line` - how many groups to include on a line
* `key_pos` - position of keygroup in output message
* `num_repeat` - How many times to repeat a number in code
* `file_out` - location to output cipher
* `padding` - Flag to pad message to length or not.

##Returns:
The encoded message, formatted accordingly.
"""

function encode(message::String, chart::DataFrameOrPath; key::Integer=0, spacing_enable::Bool=true, group_length::Integer=5,
				group_per_line::Integer=8, key_pos::Integer=1, num_repeat::Integer=3, file_out::String="", padding::Bool=true)
	#file-load; check for a filepath, and throw a catch
	if length(message)<45 && isfile(message)
		try
			message=open(message)
		catch
			error("File could not be read:\n$file_out")
		end
	end

	#message handling - remove special chars and spaces
	plaintext=join(split(lowercase(message)))

	#check if we load chart or use passed chart
	chart=(isequal(typeof(chart), String)?chart_maker(file_loc=chart):chart)

	#find the slash in the chart
	slashcode=""
	for col in eachcol(chart)
		row_pos=1
		for i in col[2]
			if isequal(i, "/")
				if !isequal(chart[:CODE][row_pos], " ") #make sure its not a single encoding row
					slashcode*=string(chart[:CODE][row_pos])
				end
				slashcode*=string(col[1])
			end
			row_pos+=1
		end
	end

	#key group handling - BROKEN/NOT IMPLEMENTED
	#=
	key_group=""
	if key>0
		for i in 1:5-length(string(key))
			key_group*="0"
		end
		key_group*=string(key)
		key=factorial(key) #too big too fast
	end
	=#

	#Generate the actual ciphertext
	ciphertext=""; output=""
	for letter in plaintext
		#check if we got a real, live number
		if isdigit(letter)
			ciphertext*=slashcode
			for i in 1:num_repeat
				ciphertext*=letter
			end
		else
			#set up the codegroup for the given letter
			code=""
			for col in eachcol(chart)
				row_pos=1
				for i in col[2]
					if isequal(i, string(letter))
						#match to a letter
						if !isequal(chart[:CODE][row_pos], " ")
							#row encoding IF in bottom 2 rows
							code*=string(chart[:CODE][row_pos])
						end
						#col encoding
						code*=string(col[1])
					end
					row_pos+=1
				end
			end
			ciphertext*=code
		end
	end

	#Key handling - could handle as you make the output, could handle here(slower) - tbd

	if padding
		#Pad current group
		if length(ciphertext)%group_length>0
			for i in 1:(group_length-(length(ciphertext)%group_length))
				ciphertext*=string(rand(UInt8)%10)
			end
		end
	end

	#Break cipher into char groups, and handle newlines
	group_space=0; group_line=0; group_count=0
	for char in ciphertext
		output*=char
		group_space+=1
		if group_space == group_length
			output*=(spacing_enable?" ":"")
			group_space=0
			group_line+=1
			group_count+=1
		end
		if group_line == group_per_line
			output*="\n"
			group_line=0
		end
	end

	if padding
		#Pad entire cipher to square off
		if group_count%group_per_line>0
			for i in 1:(group_per_line-(group_count%group_per_line))
				for j in 1:group_length
					output*=string(rand(UInt8)%10)
				end
				output*=(spacing_enable?" ":"")
			end
		end
	end

	#if file_out given, try to output to that file
	if length(file_out)>0
		try
			open(file_out, "w") do file
				write(file, output)
			end
		catch
			error("File could not be read:\n$file_out")
		end
	end

	return output
end

"""
	decode(message::String, chart::DataFrameOrPath; key_pos::Integer=1, num_repeat::Integer=3, file_out::String="")

Handles decoding of message. If a key is provided, will be subtracted from the message, modulo 10

##Parameters:
* `message` - message to decode
* `chart` - chart to use in decoding
* `key_pos` - position of keygroup - strip off leading 0's
* `num_repeat` - used to interpret numbers with proper formatting
* `file_out` - location to output plaintext message

##Returns:
Decoded message.
"""

function decode(message::String, chart::DataFrameOrPath; key_pos::Integer=1, num_repeat::Integer=3, file_out::String="")
	#check if a file, if so, load as message
	if length(message)<45 && isfile(message)
		try
			message=read(message, String)
		catch
			error("File could not be read:\n$file_loc")
		end
	end

	ciphertext=join(split(message))

	#check if we load chart or use passed chart
	chart=(isequal(typeof(chart),String)?chart_maker(file_loc=chart):chart)

	#key handling- grab key from keygroup, strip it, and run it

	#build plaintext
	plaintext=""
	index=1
	while index<length(ciphertext)
		#check if its in a row- if so, use next letter as the col, and move forward
		if in(string(ciphertext[index]), chart[:CODE]) #double encoded
			#use match in ROWCODE to determine which letter
			letter=chart[Symbol(ciphertext[index+1])][findfirst(chart[:CODE],string(ciphertext[index]))]
			index+=2
		else #single encoded
			#grab by the col symbol
			letter=chart[Symbol(ciphertext[index])][1]
			index+=1
		end
		if isequal(letter, "/") && index+1<length(ciphertext) #watch for bounds errors
			#numbers next, handle accordingly
			plaintext*=ciphertext[index+1] #just pass it as is
			index+=num_repeat
		else
			plaintext*=letter
		end
	end

	#if file_out given, try to output to that file
	if length(file_out)>0
		try
			open(file_out, "w") do file
				write(file, plaintext)
			end
		catch
			error("File could not be read:\n$file_out")
		end
	end

	return plaintext
end

end


text:
# A utility script for running a comparison benchmark.

using PkgBenchmark

import BenchmarkTools:
    BenchmarkGroup,
    TrialJudgement,
    prettydiff,
    memory,
    ratio,
    isregression

function flatten(group::BenchmarkGroup; prefix="")
    results = Pair[]

    for (name, val) in group.data
        if isa(val, BenchmarkGroup)
            append!(results, flatten(val; prefix="$prefix/$name"))
        elseif isa(val, TrialJudgement)
            push!(results, "$prefix/$name" => val)
        end
    end

    return results
end

@testset "Benchmarks" begin
    curr_commit = get(ENV, "MEMENTO_CURR_COMMIT", "16b620c")
    base_commit = get(ENV, "MEMENTO_BASE_COMMIT", "be07e98")
    result = judge(
        "Memento",
        curr_commit,
        base_commit;
        saveresults=false,
        promptsave=false,
        f=median,
        judgekwargs=Dict(
            :time_tolerance => 0.2,
            :memory_tolerance => 0.2,
        ),
    )

    flat_result = flatten(result["Memento"])

    println("Benchmark, Time, Memory")
    for element in flat_result
        println(string(
            element.first,
            ", ",
            prettydiff(time(ratio(element.second))),
            " ( ",
            time(element.second),
            " )",
            ", ",
            prettydiff(memory(ratio(element.second))),
            " ( ",
            memory(element.second),
            " )",
        ))
    end

    for element in flat_result
        @test !isregression(element.second)
    end
end


text:
using Distributions
using SplineHazard
using SplineHazard.Spline

function test_data(n=200, tmax=36)
    ## 25 patients per month
    max_entry = n / 25

    entry = rand(Uniform(0, max_entry), n)

    ## 12 months median survival time
    T = -log.(1 .-rand(Uniform(0, 1), n)) ./ (log(2)/12)

    ## approx: 10% admin cens, 10% drop-out
    C = min.(rand(Uniform(0, 3*36), n), tmax)

    time = min.(T, C)
    status = T .< C
    id = 1:n

    (id=id, entry=entry, time=time, status=status)
end

nonpar_loglik(s, data) = sum(log.(hazard(data.time[data.status], s))) - sum(cumhaz(data.time, s))

function nonpar(data, M=10000)
    s = setup_sampler(M, data.time[data.status], (sp,par) -> nonpar_loglik(sp, data))
    SplineHazard.sample!(s, M, SplineHazard.extract(s, 1))
    summary(s, data.time)
end

function nonpar_example()
    data = test_data()
    nonpar(data, 10000)
end