text:
using Test
using LinearAlgebra
using BQCESubroutine
using BQCESubroutine.Utilities
using YaoLocations: YaoLocations, plain, Locations, CtrlLocations, merge_locations
using BQCESubroutine: ctrl_offset
using YaoArrayRegister
using BenchmarkTools

include(joinpath(dirname(dirname(pathof(BQCESubroutine))), "benchmarks", "sets.jl"))
using .BenchmarkSets

SUITE = BenchmarkGroup()
INST = BenchmarkGroup()

SUITE["generic"] = @benchmarkset "$N-bits" for N in 5:5:20
    @benchmarkset "$(1<<M)x$(1<<M)" for M in [1,2,3,4]
        U = rand(ComplexF64, 1<<M, 1<<M)

        locs = Locations(unique([1,2,3, N-1, N-2, N-3])[1:M])
        locs = sort(locs)
        ctrl = CtrlLocations(setdiff(1:N, locs.storage)[1])
        @info "generating generic benchmark N=$N M=$M"
        @info "parameters" locs ctrl
        @benchmarkset "vector" begin
            @bench "gate" broutine!(st, $U, $locs) setup=(st=rand(ComplexF64, 1 << $N))
            @bench "ctrl" broutine!(st, $U, $locs, $ctrl) setup=(st=rand(ComplexF64, 1 << $N))
        end

        @benchmarkset "matrix" begin
            @bench "gate" broutine!(st, $U, $locs) setup=(st=rand(ComplexF64, 100, 1 << $N))
            @bench "ctrl" broutine!(st, $U, $locs, $ctrl) setup=(st=rand(ComplexF64, 100, 1 << $N))
        end
    end
end

INST["generic"] = @benchmarkset "$N-bits" for N in 5:5:20
    @benchmarkset "$(1<<M)x$(1<<M)" for M in [1,2,3,4]
        U = rand(ComplexF64, 1<<M, 1<<M)

        locs = Locations(unique([1,2,3, N-1, N-2, N-3])[1:M])
        locs = sort(locs)
        ctrl = CtrlLocations(setdiff(1:N, locs.storage)[1])
        @info "generating generic benchmark N=$N M=$M"
        @info "parameters" locs ctrl
        @benchmarkset "vector" begin
            @bench "gate" instruct!(st, $U, $(Tuple(locs))) setup=(st=rand(ComplexF64, 1 << $N))
            @bench "ctrl" instruct!(st, $U, $(Tuple(locs)), $(Tuple(ctrl.storage)), (1, )) setup=(st=rand(ComplexF64, 1 << $N))
        end

        @benchmarkset "matrix" begin
            @bench "gate" instruct!(st, $U, $(Tuple(locs))) setup=(st=transpose(rand(ComplexF64, 100, 1 << $N)))
            @bench "ctrl" instruct!(st, $U, $(Tuple(locs)), $(Tuple(ctrl.storage)), (1, )) setup=(st=transpose(rand(ComplexF64, 100, 1 << $N)))
        end
    end
end

SUITE["basic"] = @benchmarkset "$N-bits" for N in 5:5:20
    @benchmarkset "$G" for G in [:X, :Y, :Z, :H, :T, :S, :Sdag, :Tdag]
        @info "generating benchmark for $G gate"
        @benchmarkset "i=$i" for i in unique([1,2,3, N-1, N-2, N-3])
            locs = Locations(i)
            @bench "gate" broutine!(st, $(Val(G)), $locs) setup=(st=rand(ComplexF64, 1 << $N))

            if G in [:X, :Y, :Z]
                locs = sort(Locations((i, mod1(i+2, N))))
                @bench "multi" broutine!(st, $(Val(G)), $locs) setup=(st=rand(ComplexF64, 1 << $N))

                locs = Locations(i)
                ctrl = CtrlLocations(mod1(i+1, N))
                @bench "ctrl" broutine!(st, $(Val(G)), $locs, $ctrl) setup=(st=rand(ComplexF64, 1 << $N))
            end
        end
    end
end

INST["basic"] = @benchmarkset "$N-bits" for N in 5:5:20
    @benchmarkset "$G" for G in [:X, :Y, :Z, :H, :T, :S, :Sdag, :Tdag]
        @info "generating benchmark for $G gate"
        @benchmarkset "i=$i" for i in unique([1,2,3, N-1, N-2, N-3])
            locs = Locations(i)
            @bench "gate" instruct!(st, $(Val(G)), $(Tuple(locs))) setup=(st=rand(ComplexF64, 1 << $N))

            if G in [:X, :Y, :Z]
                locs = sort(Locations((i, mod1(i+2, N))))
                @bench "multi" instruct!(st, $(Val(G)), $(Tuple(locs))) setup=(st=rand(ComplexF64, 1 << $N))

                locs = Locations(i)
                ctrl = CtrlLocations(mod1(i+1, N))
                @bench "ctrl" instruct!(st, $(Val(G)), $(Tuple(locs)), $(Tuple(ctrl.storage)), (1, )) setup=(st=rand(ComplexF64, 1 << $N))
            end
        end
    end
end

SUITE["param"] = @benchmarkset "$N-bits" for N in 5:5:20
    @benchmarkset "$G" for G in [:Rx, :Ry, :Rz]
        @info "generating benchmark for $G gate"
        @benchmarkset "i=$i" for i in unique([1,2,3, N-1, N-2, N-3])
            locs = Locations(i)
            @bench "gate" broutine!(st, $(Val(G)), $locs, 0.2) setup=(st=rand(ComplexF64, 1 << $N))

            locs = Locations(i)
            ctrl = CtrlLocations(mod1(i+1, N))
            @bench "ctrl" broutine!(st, $(Val(G)), $locs, $ctrl, 0.2) setup=(st=rand(ComplexF64, 1 << $N))
        end
    end
end

INST["param"] = @benchmarkset "$N-bits" for N in 5:5:20
    @benchmarkset "$G" for G in [:Rx, :Ry, :Rz]
        @info "generating benchmark for $G gate"
        @benchmarkset "i=$i" for i in unique([1,2,3, N-1, N-2, N-3])
            locs = Locations(i)
            @bench "gate" instruct!(st, $(Val(G)), $(Tuple(locs)), 0.2) setup=(st=rand(ComplexF64, 1 << $N))

            locs = Locations(i)
            ctrl = CtrlLocations(mod1(i+1, N))
            @bench "ctrl" instruct!(st, $(Val(G)), $(Tuple(locs)), $(Tuple(ctrl.storage)), (1, ), 0.2) setup=(st=rand(ComplexF64, 1 << $N))
        end
    end
end

# tune!(SUITE;verbose=true)
# BenchmarkTools.save("params.json", params(SUITE));

results = map([SUITE, INST]) do SET
    loadparams!(SET, BenchmarkTools.load("params.json")[1], :evals, :samples)
    run(SET; verbose = true, seconds = 10)
end

function improve(brt, inst)
    d = Dict()
    for k in keys(inst)
        d[k] = improve(brt[k], inst[k])
    end
    return d
end

function improve(brt::BenchmarkTools.Trial, inst::BenchmarkTools.Trial)
    return minimum(inst).time/minimum(brt).time
end

d = improve(results[1], results[2])

xaxis = []
speedups = Float64[]
for k1 in keys(d), k2 in keys(d[k1]), k3 in keys(d[k1][k2]), k4 in keys(d[k1][k2][k3]), k5 in keys(d[k1][k2][k3][k4])
    r = d[k1][k2][k3][k4][k5]
    if r < 1
        @info "$k1-$k2-$k3-$k4-$k5"
        push!(speedups, r)
        push!(xaxis, "$k1-$k2-$k3-$k4-$k5")
    end
end

using UnicodePlots
UnicodePlots.barplot(xaxis, speedups, xscale=log10)


text:
module OpenQuantumBase

using DocStringExtensions

import LinearAlgebra: mul!, lmul!, axpy!, ishermitian, Hermitian, eigen, eigen!,
       tr, diag, Diagonal, norm, I, Bidiagonal
import StaticArrays: SMatrix, MMatrix, MVector, @MMatrix
import SparseArrays: sparse, issparse, spzeros, SparseMatrixCSC
import LinearAlgebra.BLAS: her!, gemm!
import QuadGK: quadgk!, quadgk
import TensorOperations: tensortrace

"""
$(TYPEDEF)

Suptertype for Hamiltonians with elements of type `T`. Any Hamiltonian object should implement two interfaces: `H(t)` and `H(du, u, p, t)`.
"""
abstract type AbstractHamiltonian{T <: Number} end

Base.eltype(::AbstractHamiltonian{T}) where {T} = T
Base.size(H::AbstractHamiltonian) = H.size
Base.size(H::AbstractHamiltonian, dim::T) where {T <: Integer} = H.size[dim]
get_cache(H::AbstractHamiltonian) = H.u_cache

"""
$(TYPEDEF)

Base for types defining Hamiltonians using dense matrices.
"""
abstract type AbstractDenseHamiltonian{T <: Number} <: AbstractHamiltonian{T} end

"""
$(TYPEDEF)

Base for types defining Hamiltonians using sparse matrices.
"""
abstract type AbstractSparseHamiltonian{T <: Number} <: AbstractHamiltonian{T} end

"""
$(TYPEDEF)

Base for types defining quantum annealing process.
"""
abstract type AbstractAnnealing{hType <: AbstractHamiltonian,uType <: Union{Vector,Matrix},} end

"""
$(TYPEDEF)

Base for types defining system bath coupling operators in open quantum system models.
"""
abstract type AbstractCouplings end

"""
$(TYPEDEF)

Base for types defining bath object.
"""
abstract type AbstractBath end

"""
$(TYPEDEF)
"""
abstract type AbstractLiouvillian end

include("base_util.jl")
include("unit_util.jl")
include("math_util.jl")
include("matrix_util.jl")
include("interpolation.jl")
include("dev_tools.jl")

include("coupling/coupling.jl")
include("coupling/interaction.jl")

include("bath/bath_util.jl")
include("bath/ohmic.jl")
include("bath/hybrid_ohmic.jl")
include("bath/spin_fluc.jl")
include("bath/custom.jl")

include("integration/cpvagk.jl")
include("integration/ext_util.jl")

include("hamiltonian/hamiltonian_base.jl")
include("hamiltonian/dense_hamiltonian.jl")
include("hamiltonian/sparse_hamiltonian.jl")
include("hamiltonian/adiabatic_frame_hamiltonian.jl")
include("hamiltonian/interp_hamiltonian.jl")
include("hamiltonian/custom_hamiltonian.jl")

include("annealing/annealing_type.jl")
include("annealing/displays.jl")

include("opensys/diffeq_liouvillian.jl")
include("opensys/redfield.jl")
include("opensys/davies.jl")
include("opensys/cgme.jl")
include("opensys/lindblad.jl")
include("opensys/stochastic.jl")
include("opensys/trajectory_jump.jl")

include("projection/projection.jl")

export SparseHamiltonian, DenseHamiltonian, AdiabaticFrameHamiltonian,
       InterpDenseHamiltonian, InterpSparseHamiltonian, CustomDenseHamiltonian
export eigen_decomp, EIGEN_DEFAULT

export temperature_2_β, temperature_2_freq, β_2_temperature, freq_2_temperature
export σx, σz, σy, σi, σ, ⊗, PauliVec, spσx, spσz, spσi, spσy, σ₊, σ₋,
       bloch_to_state, creation_operator, annihilation_operator, number_operator

export q_translate, standard_driver, local_field_term, two_local_term,
       single_clause, q_translate_state, collective_operator, hamming_weight_operator

export check_positivity, check_unitary, check_density_matrix, partial_trace,
       matrix_decompose, low_level_matrix, fidelity, inst_population, gibbs_state, purity, check_pure_state

export construct_interpolations, gradient, log_uniform

export hamiltonian_from_function, evaluate, issparse, get_cache, update_cache!,
       update_vectorized_cache!

export ConstantCouplings, TimeDependentCoupling, TimeDependentCouplings,
       CustomCouplings, collective_coupling, Interaction, InteractionSet

export Annealing, ODEParams, set_u0!, Evolution

export Lindblad, EnsembleFluctuator, DiffEqLiouvillian

export Ohmic, OhmicBath, CustomBath, CorrelatedBath, ULEBath, HybridOhmic,
       HybridOhmicBath
export correlation, polaron_correlation, γ, S, spectrum, info_freq

export τ_B, τ_SB, coarse_grain_timescale

export InplaceUnitary, EᵨEnsemble, sample_state_vector

export ProjectedSystem, project_to_lowlevel, get_dθ, concatenate

# APIs for developement tools, may move to a different repo latter
export random_ising, alt_sec_chain

end  # module OpenQuantumBase

ext: jl


text:
function counting_sort_letters(A::Array, position::Int)
	c = zeros(Int, 26)
	for i in A
		c[Int(i[position]) - 96] += 1
	end
	for i in 2:26
		c[i] += c[i - 1]
	end
	ret = Array{String, 1}(undef, length(A))
	for i in length(A) : -1 : 1
		ret[c[Int(A[i][position]) - 96]] = A[i]
		c[Int(A[i][position]) - 96] -= 1
	end
	return ret
end

function bisect_left_(A, p, r, v)
	i = p
	if p < r
	   q = floor(Int, (p + r) / 2)
	   if v >= A[q]
		   i = bisect_left_(A, p, q, v)
	   else
		   i = bisect_left_(A, q + 1, r, v)
	   end
	end
	return i
end

function bisect_left_(A, v)
	return bisect_left_(A, 1, length(A) + 1, v)
end

function counting_sort_length(A::Array)
	max_len = 0
	for i in A
		if max_len < length(i)
			max_len = length(i)
		end
	end
	c = zeros(Int, max_len + 1)

	for i in A
		c[length(i) + 1] += 1
	end

	for i in 2 : max_len + 1
		c[i] += c[i - 1]
	end

	ret = Array{String, 1}(undef, length(A))

	for i in length(A) : -1 : 1
		ret[c[length(A[i]) + 1]] = A[i]
		c[length(A[i]) + 1] -= 1
	end
	return ret
end

function copy_into!(to, from, p, r)
	for i in p : r
		to[i] = from[i - p + 1]
	end
end

function flexradix(A::Array, max_len)
	A = counting_sort_length(A)

	num_empty = 0

	for i in A
		if length(i) == 0
			num_empty += 1
		end
		if max_len < length(i)
			max_len = length(i)
		end
	end

	if max_len == 0
		return A
	end

	c = zeros(Int, max_len + 1)

	for i in A
		c[length(i) + 1] += 1
	end

	for i in max_len : -1 : 1
		c[i] += c[i + 1]
	end

	r = length(A)

	for cur_len in (max_len + 1) : -1 : 2
		p = r - c[cur_len] + 1

		if r - p > 0
			copy_into!(A, counting_sort_letters(A[p : r], cur_len - 1), p, r)
		end
	end

	copy_into!(A, counting_sort_letters(A[(num_empty + 1) : r], 1), num_empty + 1, r)

	return A
end


text:
immutable Logistic <: ContinuousUnivariateDistribution
    μ::Float64
    β::Float64

    function Logistic(μ::Float64, β::Float64)
    	β > zero(β) || error("Logistic: scale must be positive")
    	@compat new(Float64(μ), Float64(β))
    end

    @compat Logistic(μ::Real) = new(Float64(μ), 1.0)
    Logistic() = new(0.0, 1.0)
end

@distr_support Logistic -Inf Inf

#### Parameters

location(d::Logistic) = d.μ
scale(d::Logistic) = d.β

params(d::Logistic) = (d.μ, d.β)

#### Statistics

mean(d::Logistic) = d.μ
median(d::Logistic) = d.μ
mode(d::Logistic) = d.μ

std(d::Logistic) = π * d.β / sqrt3
var(d::Logistic) = (π * d.β)^2 / 3.0
skewness(d::Logistic) = 0.0
kurtosis(d::Logistic) = 1.2

entropy(d::Logistic) = log(d.β) + 2.0

#### Evaluation

zval(d::Logistic, x::Float64) = (x - d.μ) / d.β
xval(d::Logistic, z::Float64) = d.μ + z * d.β

pdf(d::Logistic, x::Float64) = (e = exp(-zval(d, x)); e / (d.β * (1.0 + e)^2))
logpdf(d::Logistic, x::Float64) = (u = -abs(zval(d, x)); u - 2.0 * log1pexp(u) - log(d.β))

cdf(d::Logistic, x::Float64) = logistic(zval(d, x))
ccdf(d::Logistic, x::Float64) = logistic(-zval(d, x))
logcdf(d::Logistic, x::Float64) = -log1pexp(-zval(d, x))
logccdf(d::Logistic, x::Float64) = -log1pexp(zval(d, x))

quantile(d::Logistic, p::Float64) = xval(d, logit(p))
cquantile(d::Logistic, p::Float64) = xval(d, -logit(p))
invlogcdf(d::Logistic, lp::Float64) = xval(d, -logexpm1(-lp))
invlogccdf(d::Logistic, lp::Float64) = xval(d, logexpm1(-lp))

function gradlogpdf(d::Logistic, x::Float64)
    e = exp(-zval(d, x))
    ((2.0 * e) / (1.0 + e) - 1.0) / d.β
end

mgf(d::Logistic, t::Real) = exp(t * d.μ) / sinc(d.β * t)

function cf(d::Logistic, t::Real)
    a = (π * t) * d.β
    a == zero(a) ? complex(one(a)) : cis(t * d.μ) * (a / sinh(a))
end

#### Sampling

rand(d::Logistic) = quantile(d, rand())


text:
function hum_dewpoint end
function hum_massfrac end
function hum_molfrac end
function rel_hum end
function hum_wetbulb end
function hum_ratio end
function dry_volume end
function dry_enthalpy end
function dry_entropy end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* Author: Maur좩cio Carvalho Maciel, Paula B. Staudt, Rafael P. Soares
#* $Id$
#*--------------------------------------------------------------------
type splitter2
	splitter2()=begin
		new(
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet stream",
				:PosX=>0.5,
				:PosY=>0,
				:Symbol=>"_{in}"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Outlet stream 1",
				:PosX=>0.25,
				:PosY=>1,
				:Symbol=>"_{out1}"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Outlet stream 2",
				:PosX=>0.75,
				:PosY=>1,
				:Symbol=>"_{out2}"
			)),
			fill(fraction(Dict{Symbol,Any}(
				:Brief=>"Distribution of Outlets",
				:Default=>0.33,
				:Symbol=>"\\phi"
			)),(2)),
			[
				:(sum(FlowRatios) = 1),
				:(Outlet1.F = Inlet.F * FlowRatios(1)),
				:(Outlet1.F + Outlet2.F = Inlet.F),
				:(Outlet1.z = Inlet.z),
				:(Outlet2.z = Inlet.z),
				:(Outlet1.P = Inlet.P),
				:(Outlet2.P = Inlet.P),
				:(Outlet1.h = Inlet.h),
				:(Outlet2.h = Inlet.h),
				:(Outlet1.T = Inlet.T),
				:(Outlet2.T = Inlet.T),
				:(Outlet1.v = Inlet.v),
				:(Outlet2.v = Inlet.v),
			],
			[
				"Normalize Flow Ratios","Flow","","Composition","","Pressure","","Enthalpy","","Temperature","","Vapourisation Fraction","",
			],
			[:Inlet,:Outlet1,:Outlet2,:FlowRatios,]
		)
	end
	Inlet::stream
	Outlet1::stream
	Outlet2::stream
	FlowRatios::Array{fraction}
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export splitter2
function setEquationFlow(in::splitter2)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
	addEquation(9)
	addEquation(10)
	addEquation(11)
	addEquation(12)
	addEquation(13)
end
function atributes(in::splitter2,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/splitter_column"
	fields[:Brief]="Splitter with 2 outlet streams"
	fields[:Info]="== Assumptions ==
*Thermodynamics equilibrium
*Adiabatic

== Specify ==
* The inlet stream
* One FlowRatios of split of the outlet streams:

	FlowRatios(i) = (Mole Flow of the outlet stream i /
				Mole Flow of the inlet stream)
				where i = 1, 2
"
	drive!(fields,_)
	return fields
end
splitter2(_::Dict{Symbol,Any})=begin
	newModel=splitter2()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
using AdaGram

vm, dict = load_model(ARGS[1])
window = int(ARGS[2])

prepare_context(ctx) = map(int32, filter(v -> v > 0, map(w -> get(dict.word2id, lowercase(w), int32(-1)), split(ctx))))

while !eof(STDIN)
	word_full, ctx_num = split(readline())
	ctx_num = int(ctx_num)
	word = word_full[1:search(word_full, '.')-1]

	for i in 1:ctx_num
		left_ctx, right_ctx = split(readline(), '\t')

		left_ctx = prepare_context(left_ctx)
		right_ctx = prepare_context(right_ctx)

		ctx = vcat(right_ctx[1:min(length(right_ctx), window)], left_ctx[max(length(left_ctx)-window+1, 1):end])

		z = disambiguate(vm, dict.word2id[word], int32(ctx))
		s = indmax(z)
		println(word_full, " ", word_full, ".", i, " ", word_full, ".senses.", s)
	end
end


text:
function test_diophantine_solutions()
    a = diophantine_solution(1, 2)
    b = diophantine_solution(7, 9)
    c = diophantine_solution(-2, 42)
    s = [a, b, c]

    @test collect(diophantine_solutions(s)) == s
    @test eltype(diophantine_solutions(s)) == DiophantineSolutionXY{Int}

    x = diophantine_solution(big(0), big(0))
    @test eltype(diophantine_solutions([x])) == DiophantineSolutionXY{BigInt}
end

function test_diophantine_solutions_nonex_noney()
    @test collect(diophantine_nonex_noney()) == []
    @test eltype(diophantine_nonex_noney()) == DiophantineSolutionXY{Int}
    @test eltype(diophantine_nonex_noney(BigInt)) == DiophantineSolutionXY{BigInt}
end

function test_diophantine_solutions_anyx_anyy()
end

function test_diophantine_solutions_onex_anyy()
    @test collect(take(diophantine_onex_anyy(42), 5)) == [
        diophantine_solution(42, 0),
        diophantine_solution(42, 1),
        diophantine_solution(42, -1),
        diophantine_solution(42, 2),
        diophantine_solution(42, -2)
    ]
    @test eltype(diophantine_onex_anyy(42)) == DiophantineSolutionXY{Int}

    @test collect(take(diophantine_onex_anyy(big(42)), 5)) == [
        diophantine_solution(big(42), big(0)),
        diophantine_solution(big(42), big(1)),
        diophantine_solution(big(42), big(-1)),
        diophantine_solution(big(42), big(2)),
        diophantine_solution(big(42), big(-2))
    ]
    @test eltype(diophantine_onex_anyy(big(42))) == DiophantineSolutionXY{BigInt}
end

function test_diophantine_solutions_anyx_oney()
    @test collect(take(diophantine_anyx_oney(42), 5)) == [
        diophantine_solution(0, 42),
        diophantine_solution(1, 42),
        diophantine_solution(-1, 42),
        diophantine_solution(2, 42),
        diophantine_solution(-2, 42)
    ]
    @test eltype(diophantine_anyx_oney(42)) == DiophantineSolutionXY{Int}

    @test collect(take(diophantine_anyx_oney(big(42)), 5)) == [
        diophantine_solution(big(0), big(42)),
        diophantine_solution(big(1), big(42)),
        diophantine_solution(big(-1), big(42)),
        diophantine_solution(big(2), big(42)),
        diophantine_solution(big(-2), big(42))
    ]
    @test eltype(diophantine_anyx_oney(big(42))) == DiophantineSolutionXY{BigInt}
end

function test_diophantine_solutions_linearx_lineary()
    @test collect(take(diophantine_linearx_lineary(2, 0, -3, 0), 5)) == [
        diophantine_solution(0, 0),
        diophantine_solution(2, -3),
        diophantine_solution(-2, 3),
        diophantine_solution(4, -6),
        diophantine_solution(-4, 6)
    ]
    @test eltype(diophantine_linearx_lineary(2, 0, -3, 0)) == DiophantineSolutionXY{Int}

    @test collect(take(diophantine_linearx_lineary(big(2), big(0), big(-3), big(0)), 5)) == [
        diophantine_solution(big(0), big(0)),
        diophantine_solution(big(2), big(-3)),
        diophantine_solution(big(-2), big(3)),
        diophantine_solution(big(4), big(-6)),
        diophantine_solution(big(-4), big(6))
    ]
    @test eltype(diophantine_linearx_lineary(big(2), big(0), big(-3), big(0))) == DiophantineSolutionXY{BigInt}
end

function test_diophantine_solutions_quadraticx_quadraticy()
    @test collect(take(diophantine_quadraticx_quadraticy(1, 2, 3, -1, 3, 2), 5)) == [
        diophantine_solution(3, 2),
        diophantine_solution(6, 4),
        diophantine_solution(2, -2),
        diophantine_solution(11, 4),
        diophantine_solution(3, -8)
    ]
    @test eltype(diophantine_quadraticx_quadraticy(1, 2, 3, -1, 3, 2)) == DiophantineSolutionXY{Int}

    @test collect(take(diophantine_quadraticx_quadraticy(big(1), big(2), big(3), big(-1), big(3), big(2)), 5)) == [
        diophantine_solution(big(3), big(2)),
        diophantine_solution(big(6), big(4)),
        diophantine_solution(big(2), big(-2)),
        diophantine_solution(big(11), big(4)),
        diophantine_solution(big(3), big(-8))
    ]
    @test eltype(diophantine_quadraticx_quadraticy(big(1), big(2), big(3), big(-1), big(3), big(2))) == DiophantineSolutionXY{BigInt}
end


text:
using PtFEM

data = Dict(
  # Beam(ndim, nst, nxe, nip, finite_element(nod, nodof), axisymmetric)
  :struc_el => Beam(2, 1, 5, 1, :x, Line(2, 2), false),
  :properties => [
    1.924e4 0.2;
    1.924e4 0.6;
    1.924e4 1.0;
    1.924e4 1.4;
    1.924e4 1.8],
  :etype => [1, 2, 3, 4, 5],
  :x_coords => [0.0, 2.0, 4.0, 6.0, 8.0, 10.0],
  :loaded_nodes => [(1, [1.0 0.0])]
)

data |> display
println()

@time fem, dis_df, fm_df = p43(data)
println()

display(dis_df)
println()
display(fm_df)


text:
using SafeTestsets

@safetestset "Dispersion tests" begin
    include("tests_dispersion.jl")
end
@safetestset "Raytracing tests" begin
    include("tests_raytracing.jl")
end
@safetestset "LsqFit tests" begin
    include("tests_lsqfit.jl")
end


text:
module MultiFieldTestsSeq
using PartitionedArrays
include("../MultiFieldTests.jl")
prun(MultiFieldTests.main,sequential,(2,2))
end # module


text:
@testset "MerCount" begin
    v = DNAMer{4}[
        mer"AAAA", mer"AAAA", mer"AAAA", mer"ATAG",
        mer"GGGG", mer"GGGG", mer"GGGT", mer"AGGT",
        mer"AGGT", mer"AGGT", mer"AGGT"
    ]

    v2 = DNAMer{4}[
        mer"AAAA", mer"AAAA", mer"AAAA", mer"ATAG",
        mer"GGGG", mer"GGGG", mer"GGGT", mer"AGGT",
        mer"ACGT", mer"AAGT", mer"AAGT"
    ]

    function test_merge_into(a, b, a′)
        merge_into!(a, b)
        return isempty(b) && a == a′
    end

    @test collapse_into_counts!(v) == [
        MerCount{DNAMer{4}}(mer"AAAA", 3),
        MerCount{DNAMer{4}}(mer"AGGT", 4),
        MerCount{DNAMer{4}}(mer"ATAG", 1),
        MerCount{DNAMer{4}}(mer"GGGG", 2),
        MerCount{DNAMer{4}}(mer"GGGT", 1)
    ]

    @test collapse_into_counts!(v2) == [
        MerCount{DNAMer{4}}(mer"AAAA", 3),
        MerCount{DNAMer{4}}(mer"AAGT", 2),
        MerCount{DNAMer{4}}(mer"ACGT", 1),
        MerCount{DNAMer{4}}(mer"AGGT", 1),
        MerCount{DNAMer{4}}(mer"ATAG", 1),
        MerCount{DNAMer{4}}(mer"GGGG", 2),
        MerCount{DNAMer{4}}(mer"GGGT", 1)
    ]

    @test test_merge_into(collapse_into_counts!(v), collapse_into_counts!(v2), [
        MerCount{DNAMer{4}}(mer"AAAA", 6),
        MerCount{DNAMer{4}}(mer"AAGT", 2),
        MerCount{DNAMer{4}}(mer"ACGT", 1),
        MerCount{DNAMer{4}}(mer"AGGT", 5),
        MerCount{DNAMer{4}}(mer"ATAG", 2),
        MerCount{DNAMer{4}}(mer"GGGG", 4),
        MerCount{DNAMer{4}}(mer"GGGT", 2)
    ])

    seq = LongSequence{DNAAlphabet{2}}("ACTGAAGAGC")
    @test collect_mers(DNAMer{4}, NONCANONICAL, seq) == [mer"ACTG", mer"CTGA", mer"TGAA", mer"GAAG", mer"AAGA", mer"AGAG", mer"GAGC"]
    @test collect_mers(DNAMer{4}, NONCANONICAL, [seq, seq]) == repeat([mer"ACTG", mer"CTGA", mer"TGAA", mer"GAAG", mer"AAGA", mer"AGAG", mer"GAGC"], 2)

    @test freq(MerCount{DNAMer{4}}(mer"AAAA", 3)) == 3
    @test freq(MerCount{DNAMer{4}}(mer"AAAA", 3)) === 0x03
    @test freq(Int, MerCount{DNAMer{4}}(mer"AAAA", 3)) === 3
    @test freq(Float64, MerCount{DNAMer{4}}(mer"AAAA", 3)) == 3.0
end


text:
# -*- encoding: utf-8 -*-

export ADC, optimized_adc, adc_response, adc_inv_response, adc_filter

@doc raw"""
A structure representing the configuration of an Analogue-to-Digital Converter
(ADC). To simplify things, we assume that the ADC is fed with some *temperature*,
instead of a voltage. In this case we can avoid dealing with the calibration from
Volt to Kelvin, which is seldom useful in simulations.

The equation used to model the ADC is the following:

    x(V) = round(gain_k_over_adu * (V - offset_k)) + zero_point_adu

and the output is clipped within a user-defined range. The following fields are
available:

- `offset_k`
- `gain_k_over_adu`
- `zero_point_adu`
- `min_output_adu`: output values below this will be clipped
- `max_output_adu`: output values above this will be clipped
- `non_linearities_x_adu`: see below
- `non_linearities_y_adu`: see below

Non linearities must be specified using the two arrays `non_linearities_x_adu`
and `non_linearities_y_adu`. For each input voltage fed to the ADC, a non
linearity is applied to the output according to the following algorithm:

- The function computes the ideal output (i.e., without non linearities);
- It checks if the value of the ideal output is found in the array
  `non_linearities_x_adu`;
- If the value is found, add the corresponding value in the array
  `non_linearities_y_adu` (i.e., the element with the same index as the element
  in `non_linearities_x_adu`) to the ideal output;
- If the value is not found, but it falls within two consecutive values,
  use a linear interpolation;
- If the value is smaller than the first element or larger than the last element
  in `non_linearities_x_adu`, do not apply any correction to the output.

It is **fundamental** that `non_linearities_x_adu` is sorted in ascending
order, and that the order of the in elements `non_linearities_y_adu` matches
the order in `non_linearities_x_adu`.

"""
mutable struct ADC
    offset_k::Float64
    gain_k_over_adu::Float64
    zero_point_adu::Float64
    min_output_adu::Int
    max_output_adu::Int
    non_linearities_x_adu::Array{Float64,1}
    non_linearities_y_adu::Array{Float64,1}

    ADC(;
        offset_k = 0.0,
        gain_k_over_adu = 1.0,
        zero_point_adu = 0.0,
        min_output_adu = -2^19,
        max_output_adu = 2^19,
        non_linearities_x_adu = Float64[],
        non_linearities_y_adu = Float64[]) = new(offset_k, gain_k_over_adu, zero_point_adu, min_output_adu, max_output_adu, non_linearities_x_adu, non_linearities_y_adu)
end

@doc raw"""

    optimized_adc(; min_input_k = 0.0, max_input_k = 100.0, dynamic_range = 0.35, nbits = 20, non_linearities_x_adu = Float64[], non_linearities_y_adu = Float64[])

Return an object of type `ADC` that is optimized to measure temperatures between
`min_input_k` and `max_input_k`. The ADC is configured to return *signed*
numbers in the range -2^nbits…2^nbits.

The function accepts the following keywords:

- `min_input_k`: temperature that should trigger the lowest output within the
  dynamic range
- `max_input_k`: temperature that should trigger the lowest output within the
  dynamic range
- `dynamic_range`: pure number in the interval 0…1 that specifies the dynamic
  range of the output
- `nbits`: number of bits used by the ADC
- `non_linearities_x_adu`: same parameter used in the constructor for `ADC`
- `non_linearities_y_adu`: same parameter used in the constructor for `ADC`

"""
function optimized_adc(; min_input_k = 0.0, max_input_k = 100.0, dynamic_range = 0.35, nbits = 20, non_linearities_x_adu = Float64[], non_linearities_y_adu = Float64[])
    min_output_adu = -2^(nbits - 1)
    max_output_adu = 2^(nbits - 1)
    ADC(offset_k = -min_input_k,
        gain_k_over_adu = (max_input_k - min_input_k) / (2^nbits),
        zero_point_adu = min_output_adu,
        min_output_adu = min_output_adu,
        max_output_adu = max_output_adu,
        non_linearities_x_adu = non_linearities_x_adu,
        non_linearities_y_adu = non_linearities_y_adu,
    )
end

@doc raw"""
    adc_response(adc::ADC, input_k; include_non_linearities = true)

Simulate the response of an ADC on some voltage `input_k`. If
`include_non_linearities` is `true`, non-linearities specified in `adc` will be
considered; otherwise, the ADC will be assumed to be ideal.

See also [`adc_inv_response`](@ref) for the (pseudo)inverse function.
"""
function adc_response(adc::ADC, input_k; include_non_linearities = true)
    output = adc.gain_k_over_adu * (input_k - adc.offset_k) + adc.zero_point_adu

    # Let's use some shorthands here
    x = adc.non_linearities_x_adu
    y = adc.non_linearities_y_adu

    if include_non_linearities && length(x) > 0 && (x[1] ≤ output ≤ x[end])
        idx = searchsortedlast(x, output)
        if idx > 0
            δ = if idx == length(x)
                y[idx]
            else
                # Interpolate between idx and idx + 1
                y[idx] + (output - x[idx]) * (y[idx] - y[idx + 1]) / (x[idx] - x[idx + 1])
            end
            output += δ
        end
    end

    clamp(round(Int, output),
        adc.min_output_adu,
        adc.max_output_adu,
    )
end

@doc raw"""
    adc_inv_response(adc::ADC, input_adu)

Apply the inverse transformation of an ADC to get some voltage from a digital
measurement in ADU.This function assumes that the ADC is ideal, i.e., that it
does not have non-linearities.

See also [`adc_response`](@ref).
"""
function adc_inv_response(adc::ADC, input_adu)
    (input_adu - adc.zero_point_adu) / adc.gain_k_over_adu + adc.offset_k
end

@doc raw"""
    adc_filter(adc::ADC, input_k; include_non_linearities = true)

Simulate the measurement of the temperature `input_k` through the ADC `adc`. The
result is still a temperature, after it has been converted into ADUs by the ADC
and then converted back to a temperature again. If non-linearities are specified
in `adc`, they will be applied only to the K→ADU transformation.
"""
function adc_filter(adc::ADC, input_k; include_non_linearities = true)
    adu = adc_response(adc, input_k, include_non_linearities = include_non_linearities)
    adc_inv_response(adc, adu)
end


text:
@testset "test_timeevol" begin

@test ControlSystems.common_timeevol(Discrete(1.5), Discrete(1.5)) == Discrete(1.5)
@test ControlSystems.common_timeevol(Continuous(), Continuous()) == Continuous()

@test_throws ErrorException ControlSystems.common_timeevol(Continuous(), Discrete(1.5))
@test_throws ErrorException ControlSystems.common_timeevol(Discrete(1.0), Discrete(1.5))

sys = ss(-1, 1, 1, 0)
iscontinuous(sys) == true
isdiscrete(sys) == false

sysd = ss(-1, 1, 1, 0, 0.1)
iscontinuous(sysd) == false
isdiscrete(sysd) == true

end


text:
using DataFrames
using CategoricalArrays

# DataFrame
attenu = rcopy(DataFrame,reval(:attenu))
@test isa(attenu,DataFrame)
@test size(attenu) == (182,5)
@test rcopy(rcall(:dim,RObject(attenu))) == [182,5]
@test rcopy(rcall(:dim, RObject(attenu[1:2, :]))) == [2, 5]
@test rcopy(rcall(:dim, RObject(view(attenu, 1:2)))) == [2, 5]
dist = attenu[:dist]
@test isa(dist,Array{Float64})
station = attenu[:station]
@test isa(station, CategoricalArray)

# single row dataframe
@test size(rcopy(R"data.frame(a=1,b=2)")) == (1, 2)

# issue #186
df = R"""data.frame(dates = as.Date(c("2017-04-14", "2014-04-17")))"""
@test eltype(rcopy(df)[:dates]) == Date


text:
using SeisMain
using Test

param2 = Dict(:style=>"mxmyhaz", :min_imx=>10,:max_imx=>100, :min_imy=>35, :max_imy=>45,
               :min_ih=>1, :max_ih=>6, :min_iaz=>0, :max_iaz=>7);

param3 = Dict( :ix1_WL=>30,:ix1_WO=>5);

file_bin="section_bin"           # Run first Demo_Utils_1_Binning.ipynb
patch_out="patch"                # will make files patch1, patch2,....
file_final ="section_bin_final"  # file put back from patches via code SeisUnPatch

patch_out,npatch=SeisPatch(file_bin,"patch"; param2...,param3... );

SeisUnPatch(patch_out,file_final; param2...,param3...,nt=251);

db,hb,eb=SeisRead(file_bin);

df,hf,ef=SeisRead(file_final);

alpha = maximum(db-df)

@test alpha < 0.1


text:
function eval_overlapᵥ(
    xₜ::Vector{Float64},
    x::Matrix{Float64},
    fs_counts::Vector{Int64},
    ol_maxₜ::Vector{Float64},
    v::Int64,
)
    xₜ = xₜ .* fs_counts
    cons = LinearAlgebra.BLAS.gemv('T', x, xₜ) - ol_maxₜ
    cons[v] = 0
    return sum_pos(cons)::Float64
end

function eval_overlap(
    x::Matrix{Float64},
    fs_counts::Matrix{Float64},
    ol_max::Matrix{Float64},
    T::Int64,
    ol2::Vector{Float64},
)
    fs_counts = fs_counts .* x
    ol = _gemmblas(x, fs_counts, ol_max)
    for v2 = 1:T
        olₜ = ol[:, v2]
        olₜ[v2] = zero(Float64)
        ol2[v2] = sum_pos(olₜ)
    end
    return ol2::Vector{Float64}
end

function eval_exp_scoreₜ(x_Iₜ::Vector{Float64}, expected_scoreₜ::ExpectedScore) #x = I
    es = (expected_scoreₜ.val * x_Iₜ) ./ sum(x_Iₜ)
    # _gemvblasT(
    #     expected_scoreₜ.val,
    #     x_Iₜ,
    #     size(expected_scoreₜ.val, 2)
    # ) ./ sum(x_Iₜ)
    es = vcat((es .- expected_scoreₜ.max)..., (-es .+ expected_scoreₜ.min)...)
    return es
end

function eval_cons(xₜ::Vector{Float64}, Aₜ::Matrix{Float64}, bₜ::Vector{Float64})
    cons = copy(bₜ)
    cons = _gemvblas(Aₜ, xₜ, cons, size(xₜ, 1))
    return cons::Vector{Float64}
end

function check_feas(
    fs::FriendSets,
    consₜ::Constraint,
    xₜ::Vector{Float64},
    n_fs::Int64,
    n_items::Int64,
)
    es = 0
    cons = 0
    if size(consₜ.constr_A, 1) > 0
        cons = eval_cons(xₜ, consₜ.constr_A, consₜ.constr_b)
    end
    if n_items > n_fs
        x_Iₜ = fs_to_items(xₜ, n_items, fs.items)
    else
        x_Iₜ = copy(xₜ)
    end
    if size(consₜ.expected_score.val, 1) > 0
        cons = vcat(cons, eval_exp_scoreₜ(x_Iₜ, consₜ.expected_score))
    end
    return sum_pos(cons)::Float64, x_Iₜ::Vector{Float64}
end


text:
import LinearAlgebra: I, Hermitian, tr, pinv
import ForneyLab: unsafeCov, unsafeMean, unsafePrecision, VariateType,
				  collectNaiveVariationalNodeInbounds, assembleClamp!, ultimatePartner
import Zygote: gradient

export ruleVariationalNLARXOutNPPPPP,
       ruleVariationalNLARXIn1PNPPPP,
       ruleVariationalNLARXIn2PPNPPP,
       ruleVariationalNLARXIn3PPPNPP,
	   ruleVariationalNLARXIn4PPPPNP,
	   ruleVariationalNLARXIn5PPPPPN

# Autoregression order and bookkeeping matrices
order = Nothing
S = Nothing
s = Nothing

# Approximating points
approxx = 0.0
approxθ = 0.0

function defineOrder(dim::Int64)
    global order, s, S

	# Set autoregression order
    order = dim

	# Set bookkeeping matrices
    s = uvector(order)
    S = shift(order)
end

function ruleVariationalNLARXOutNPPPPP(g :: Function,
									   marg_y :: Nothing,
                                       marg_x :: ProbabilityDistribution{Multivariate},
                                       marg_θ :: ProbabilityDistribution{Multivariate},
                                       marg_η :: ProbabilityDistribution{Univariate},
                                       marg_u :: ProbabilityDistribution{Univariate},
                                       marg_γ :: ProbabilityDistribution{Univariate})

    # Expectations of incoming marginal beliefs
    mx = unsafeMean(marg_x)
    mθ = unsafeMean(marg_θ)
    mη = unsafeMean(marg_η)
    mu = unsafeMean(marg_u)
    mγ = unsafeMean(marg_γ)

    # Check order
	if order == Nothing
		defineOrder(length(mx))
	end

	# Update approximating points
	global approxx = mx
	global approxθ = mθ

    # Map transition noise to matrix
    mW = wMatrix(mγ, order)

    # Parameters of outgoing message
	Φ = mW
    ϕ = S*mx + s*g(mθ, mx) + s*mη*mu

	# Set outgoing message
	return Message(Multivariate, GaussianMeanPrecision, m=ϕ, w=Φ)
end

function ruleVariationalNLARXIn1PNPPPP(g :: Function,
									   marg_y :: ProbabilityDistribution{Multivariate},
                                       marg_x :: Nothing,
                                       marg_θ :: ProbabilityDistribution{Multivariate},
                                       marg_η :: ProbabilityDistribution{Univariate},
                                       marg_u :: ProbabilityDistribution{Univariate},
                                       marg_γ :: ProbabilityDistribution{Univariate})

    # Expectations of marginal beliefs
    my = unsafeMean(marg_y)
    mθ = unsafeMean(marg_θ)
    Vθ = unsafeCov(marg_θ)
    mη = unsafeMean(marg_η)
    mu = unsafeMean(marg_u)
    mγ = unsafeMean(marg_γ)

	# Check order
	if order == Nothing
		defineOrder(length(mx))
	end

	# Check approximating points
	if approxx == Nothing
		error("Approximating point for x not initialized")
	end

	# Gradient of supplied nonlinear function
	Jθ, Jx = gradient(g, mθ, approxx)

	# Map transition noise to matrix
	mW = wMatrix(mγ, order)

	# Parameters of outgoing message
    Φ = S'*mW*S + Jx*s'*mW*s*Jx'
    ϕ = (S + s*Jx')'*mW*(my - s*mη*mu) - Jx*(g(mθ,approxx) - Jx'*approxx)*mγ

	# Update global approximating point
	global approxx = inv(Φ + 1e-6*Matrix{Float64}(I, size(Φ)))*ϕ
	global approxθ = mθ

    return Message(Multivariate, GaussianWeightedMeanPrecision, xi=ϕ, w=Φ)
end

function ruleVariationalNLARXIn2PPNPPP(g :: Function,
									   marg_y :: ProbabilityDistribution{Multivariate},
                                       marg_x :: ProbabilityDistribution{Multivariate},
                                       marg_θ :: Nothing,
							  	       marg_η :: ProbabilityDistribution{Univariate},
								       marg_u :: ProbabilityDistribution{Univariate},
                                       marg_γ :: ProbabilityDistribution{Univariate})

	# Expectations of marginal beliefs
	my = unsafeMean(marg_y)
	mx = unsafeMean(marg_x)
	Vx = unsafeCov(marg_x)
	mη = unsafeMean(marg_η)
	mu = unsafeMean(marg_u)
	mγ = unsafeMean(marg_γ)

	# Check order
	if order == Nothing
		defineOrder(length(mx))
	end

	# Check approximating points
	if approxθ == Nothing
		error("Approximating point for θ not initialized")
	end

	# Gradient of supplied nonlinear function
	Jθ, Jx = gradient(g, approxθ, mx)

	# Map transition noise to matrix
	mW = wMatrix(mγ, order)

    # Parameters of outgoing message
    Φ = mγ*Jθ*Jθ'
    ϕ = Jθ*s'*mW*(my - s*mη*mu - s*(g(approxθ, mx) - Jθ'*approxθ))

	# Update global approximating point
	global approxx = mx
	global approxθ = inv(Φ + 1e-6*Matrix{Float64}(I, size(Φ)))*ϕ

    return Message(Multivariate, GaussianWeightedMeanPrecision, xi=ϕ, w=Φ)
end

function ruleVariationalNLARXIn3PPPNPP(g :: Function,
									   marg_y :: ProbabilityDistribution{Multivariate},
                                       marg_x :: ProbabilityDistribution{Multivariate},
                                       marg_θ :: ProbabilityDistribution{Multivariate},
								       marg_η :: Nothing,
								       marg_u :: ProbabilityDistribution{Univariate},
                                       marg_γ :: ProbabilityDistribution{Univariate})

 	# Expectations of marginal beliefs
	mθ = unsafeMean(marg_θ)
    my = unsafeMean(marg_y)
    mx = unsafeMean(marg_x)
	mu = unsafeMean(marg_u)
	vu = unsafeCov(marg_u)
	mγ = unsafeMean(marg_γ)

	# Update global approximating point
	global approxx = mx
	global approxθ = mθ

	# Check order
	if order == Nothing
		defineOrder(length(mx))
	end

    # Map transition noise to matrix
    mW = wMatrix(mγ, order)

	# Parameters of outgoing message
	Φ = mγ*(mu^2 + vu)
    ϕ = (mu*s')*mW*(my - (S*mx + s*g(mθ, mx)))

	return Message(Univariate, GaussianWeightedMeanPrecision, xi=ϕ, w=Φ)
end

function ruleVariationalNLARXIn4PPPPNP(g :: Function,
									   marg_y :: ProbabilityDistribution{Multivariate},
                                       marg_x :: ProbabilityDistribution{Multivariate},
                                       marg_θ :: ProbabilityDistribution{Multivariate},
							  	       marg_η :: ProbabilityDistribution{Univariate},
								       marg_u :: Nothing,
                                       marg_γ :: ProbabilityDistribution{Univariate})

 	# Expectations of marginal beliefs
	mθ = unsafeMean(marg_θ)
    my = unsafeMean(marg_y)
    mx = unsafeMean(marg_x)
	mη = unsafeMean(marg_η)
	vη = unsafeCov(marg_η)
	mγ = unsafeMean(marg_γ)

	# Update global approximating point
	global approxx = mx
	global approxθ = mθ

	# Check order
	if order == Nothing
		defineOrder(length(mx))
	end

    # Map transition noise to matrix
    mW = wMatrix(mγ, order)

	# Parameters of outgoing message
	Φ = mγ*(mη^2 + vη)
    ϕ = (mη*s')*mW*(my - (S*mx + s*g(mθ, mx)))

	return Message(Univariate, GaussianWeightedMeanPrecision, xi=ϕ, w=Φ)
end

function ruleVariationalNLARXIn5PPPPPN(g :: Function,
									   marg_y :: ProbabilityDistribution{Multivariate},
                                       marg_x :: ProbabilityDistribution{Multivariate},
                                       marg_θ :: ProbabilityDistribution{Multivariate},
							  	       marg_η :: ProbabilityDistribution{Univariate},
								       marg_u :: ProbabilityDistribution{Univariate},
                                       marg_γ :: Nothing)

    # Expectations of marginal beliefs
	mθ = unsafeMean(marg_θ)
    my = unsafeMean(marg_y)
    mx = unsafeMean(marg_x)
	mη = unsafeMean(marg_η)
	mu = unsafeMean(marg_u)
    Vθ = unsafeCov(marg_θ)
    Vy = unsafeCov(marg_y)
    Vx = unsafeCov(marg_x)
	vη = unsafeCov(marg_η)
	vu = unsafeCov(marg_u)

	# Update global approximating point
	global approxx = mx
	global approxθ = mθ

	# Check order
	if order == Nothing
		defineOrder(length(mx))
	end

	# Gradient of supplied nonlinear function
	Jθ, Jx = gradient(g, mθ, approxx)

	# Convenience variables
	Aθx = S*mx + s*g(mθ, mx)

	# Intermediate terms
	term1 = (my*my' + Vy)[1,1]
	term2 = -(Aθx*my')[1,1]
	term3 = -((s*mη*mu)*my')[1,1]
	term4 = -(my*Aθx')[1,1]
	term5 = (mx'*S'*S*mx)[1,1] + (S*Vx*S')[1,1] + g(mθ, mx)^2 + Jx'*Vx*Jx + Jθ'*Vθ*Jθ
	term6 = ((s*mη*mu)*Aθx')[1,1]
	term7 = -(my*(s*mη*mu)')[1,1]
	term8 = (Aθx*(s*mη*mu)')[1,1]
	term9 = (mu^2 + vu)*(mη^2 + vη)

	# Parameters of outgoing message
	a = 3/2
    B = (term1 + term2 + term3 + term4 + term5 + term6 + term7 + term8 + term9)/2

	return Message(Gamma, a=a, b=B)
end

function collectNaiveVariationalNodeInbounds(node::NLatentAutoregressiveX, entry::ScheduleEntry)
	inbounds = Any[]

	# Push function to calling signature (g needs to be defined in user scope)
	push!(inbounds, Dict{Symbol, Any}(:g => node.g, :keyword => false))

    target_to_marginal_entry = currentInferenceAlgorithm().target_to_marginal_entry

    for node_interface in entry.interface.node.interfaces
        inbound_interface = ultimatePartner(node_interface)
        if node_interface === entry.interface
            # Ignore marginal of outbound edge
            push!(inbounds, nothing)
        elseif (inbound_interface != nothing) && isa(inbound_interface.node, Clamp)
            # Hard-code marginal of constant node in schedule
            push!(inbounds, assembleClamp!(inbound_interface.node, ProbabilityDistribution))
        else
            # Collect entry from marginal schedule
            push!(inbounds, target_to_marginal_entry[node_interface.edge.variable])
        end
    end

    return inbounds
end

max_stars_repo_name: wmkouw/NLARX

max_issues_repo_name: wmkouw/NLARX

max_forks_repo_name: wmkouw/NLARX


text:
using IntArrays
using Test
using Random

Random.seed!(12345)

const Ts = (UInt8, UInt16, UInt32, UInt64)

@testset "conversion" begin
    @testset "IntVector" begin
        data = [0x00, 0x01, 0x02]
        ivec = IntArray{2}(data)
        @test typeof(ivec) == IntArray{2,UInt8,1}
        @test typeof(ivec) == IntVector{2,UInt8}

        ivec = convert(IntVector{2}, data)
        @test typeof(ivec) == IntArray{2,UInt8,1}
        @test typeof(ivec) == IntVector{2,UInt8}
    end
    @testset "IntMatrix" begin
        data = [0x00 0x01; 0x02 0x03; 0x04 0x05]
        imat = IntArray{3}(data)
        @test typeof(imat) == IntArray{3,UInt8,2}
        @test typeof(imat) == IntMatrix{3,UInt8}

        imat = convert(IntMatrix{3}, data)
        @test typeof(imat) == IntArray{3,UInt8,2}
        @test typeof(imat) == IntMatrix{3,UInt8}
    end
    @testset "three-dimensional array" begin
        data = rand(0x00:0x03, 2, 3, 4)
        iarr = IntArray{2}(data)
        @test typeof(iarr) == IntArray{2,UInt8,3}
    end
end

@testset "construction" begin
    @testset "IntVector" begin
        ivec = IntArray{9,UInt16}(10)
        @test typeof(ivec) == IntArray{9,UInt16,1}
        @test size(ivec) == (10,)
        @test length(ivec) == 10
    end
    @testset "IntMatrix" begin
        imat = IntMatrix{3,UInt8}(2, 3)
        @test typeof(imat) == IntArray{3,UInt8,2}
        @test size(imat) == (2, 3)
    end
    @testset "invalid width" begin
        @test_throws Exception IntArray{9,UInt8,1}(4)
        @test_throws Exception IntArray{20,UInt16,1}(4)
        @test_throws Exception IntArray{50,UInt32,1}(4)
    end
    @testset "mmap" begin
        @test typeof(IntArray{2,UInt8}(10, true)) == IntArray{2,UInt8,1}
    end
end

@testset "similar" begin
    @testset "IntVector" begin
        data = rand(0x00:0x02, 10)
        ivec = IntVector{2}(data)
        @test similar(ivec) !== ivec
        @test size(similar(ivec)) == (10,)
        @test typeof(similar(ivec)) == IntVector{2,UInt8}
        @test size(similar(ivec, UInt16)) == (10,)
        @test typeof(similar(ivec, UInt16)) == IntVector{2,UInt16}
        @test size(similar(ivec, UInt16, (20,))) == (20,)
        @test typeof(similar(ivec, UInt16, (20,))) == IntVector{2,UInt16}
    end
end

@testset "getindex" begin
    @testset "IntVector" begin
        data = [0x00, 0x01, 0x02, 0x03, 0x04]
        ivec = IntArray{3}(data)
        @test_throws BoundsError ivec[0]
        @test ivec[1] == 0x00
        @test ivec[2] == 0x01
        @test ivec[3] == 0x02
        @test ivec[4] == 0x03
        @test ivec[5] == 0x04
        @test_throws BoundsError ivec[6]
    end
    @testset "IntMatrix" begin
        # 2x3
        data = [0x00 0x01 0x02; 0x03 0x04 0x05]
        imat = IntArray{4}(data)
        # linear indexing
        @test_throws BoundsError imat[0]
        @test imat[1] == 0x00
        @test imat[2] == 0x03
        @test imat[3] == 0x01
        @test imat[4] == 0x04
        @test imat[5] == 0x02
        @test imat[6] == 0x05
        @test_throws BoundsError imat[7]
        # tuples
        @test_throws BoundsError imat[0,1]
        @test imat[1,1] == 0x00
        @test imat[1,2] == 0x01
        @test imat[1,3] == 0x02
        @test imat[2,1] == 0x03
        @test imat[2,2] == 0x04
        @test imat[2,3] == 0x05
        @test_throws BoundsError imat[1,4]
        @test_throws BoundsError imat[2,4]
        @test_throws BoundsError imat[3,1]
    end
end

@testset "setindex!" begin
    @testset "IntVector" begin
        data = [0x00, 0x01, 0x02]
        ivec = IntArray{2}(data)
        @test_throws BoundsError ivec[0] = 0x00
        ivec[1] = 0x01
        @test ivec[1] == 0x01
        ivec[2] = 0x02
        @test ivec[2] == 0x02
        ivec[3] = 0x03
        @test ivec[3] == 0x03
        @test_throws BoundsError ivec[4] = 0x01
    end
    @testset "IntMatrix" begin
        data = [0x00 0x00 0x00; 0x00 0x00 0x00]
        imat = IntArray{2}(data)
        # linear
        @test_throws BoundsError imat[0] = 0x00
        imat[2] = 0x01
        @test imat[2] == 0x01
        imat[2] = 0x00
        @test imat[2] == 0x00
        @test_throws BoundsError imat[7] = 0x00
        # tuple
        imat[1,1] = 0x01
        @test imat[1,1] == 0x01
        imat[1,3] = 0x03
        @test imat[1,3] == 0x03
        @test_throws BoundsError imat[1,4] = 0x00
        @test_throws BoundsError imat[3,1] = 0x01
    end
end

@testset "comparison" begin
    a = IntVector{2}([0x00, 0x01, 0x02, 0x03])
    b = IntVector{2}([0x00, 0x01, 0x02, 0x03])
    c = IntVector{4}([0x00, 0x01, 0x02, 0x03])
    d = IntVector{2}([0x00, 0x02, 0x02, 0x03])
    e = IntMatrix{2}([0x00  0x02; 0x01  0x03])
    @test a == b
    @test a == c
    @test a != d
    @test a != e
end

@testset "sizeof" begin
    @testset "smaller bits" begin
        n = 100
        data = rand(0x00:0x01, n)
        @test sizeof(IntVector{1}(data)) < sizeof(data)
        @test sizeof(IntVector{2}(data)) < sizeof(data)
        @test sizeof(IntVector{3}(data)) < sizeof(data)
        @test sizeof(IntVector{4}(data)) < sizeof(data)
    end
end

@testset "copy" begin
    @testset "same length" begin
        data = rand(0x00:0x03, 10)
        ivec = IntVector{2,UInt8}(data)
        @test copy(ivec) == ivec
        @test typeof(copy(ivec)) == typeof(ivec)
        @test size(copy(ivec)) == size(ivec)
        ivec′ = IntVector{2,UInt8}(10)
        @test copy!(ivec′, ivec) === ivec′
        @test ivec′ == ivec
    end
    @testset "larger" begin
        data = rand(0x00:0x03, 10)
        ivec = IntVector{2,UInt8}(data)
        ivec′ = IntVector{2,UInt8}(20)
        copy!(ivec′, ivec)
        @test ivec′[1:10] == ivec
    end
    @testset "smaller" begin
        data = rand(0x00:0x03, 10)
        ivec = IntVector{2,UInt8}(data)
        ivec′ = IntVector{2,UInt8}(5)
        @test_throws BoundsError copy!(ivec′, ivec)
    end
end

@testset "fill!" begin
    @testset "IntVector" begin
        n = 100
        data = rand(0x00:0x03, n)
        for x in 0x00:0x03
            ivec = IntVector{2}(data)
            @test fill!(ivec, x) === ivec
            @test ivec == ones(UInt8, n) * x
        end
    end
    @testset "IntMatrix" begin
        m, n = 10, 11
        data = rand(0x00:0x03, (m, n))
        for x in 0x00:0x03
            #imat = IntMatrix{2}(m, n)
            #@test fill!(imat, x) === imat == true
            #@test imat == ones(UInt8, (m, n)) * x
        end
    end
end

@testset "reverse" begin
    ivec = IntVector{2,UInt8}()
    @test reverse!(ivec) === ivec
    @test isempty(ivec)

    ivec = IntVector{2}([0x00])
    @test reverse!(ivec) === ivec
    @test ivec == [0x00]

    ivec = IntVector{2}([0x00, 0x01])
    @test reverse!(ivec) === ivec
    @test ivec == [0x01, 0x00]

    ivec = IntVector{2}([0x00, 0x01, 0x02])
    @test reverse!(ivec) === ivec
    @test ivec == [0x02, 0x01, 0x00]

    ivec = IntVector{2}([0x00, 0x01, 0x02, 0x03])
    @test reverse!(ivec) === ivec
    @test ivec == [0x03, 0x02, 0x01, 0x00]
end

@testset "push!/pop!" begin
    ivec = IntVector{4,UInt8}()
    @test length(ivec) == 0
    @test push!(ivec, 3) === ivec
    @test length(ivec) == 1
    @test ivec[end] == 3
    @test pop!(ivec) == 3
    @test length(ivec) == 0
end

@testset "append!" begin
    ivec = IntVector{4}([0x00])
    @test append!(ivec, [0x01, 0x02]) === ivec
    @test ivec == [0x00, 0x01, 0x02]
    append!(ivec, [3, 4])
    @test ivec == [0x00, 0x01, 0x02, 0x03, 0x04]
end

@testset "radixsort" begin
    data = UInt8[]
    ivec = IntVector{1}(data)
    @test issorted(radixsort(ivec))

    data = [0x00]
    ivec = IntVector{2}(data)
    @test issorted(radixsort(ivec))

    data = [0x01, 0x00]
    ivec = IntVector{2}(data)
    @test issorted(radixsort(ivec))

    n = 101
    data = rand(0x00:0x01, n)
    ivec = IntVector{1}(data)
    @test issorted(radixsort(ivec))
    @test radixsort!(ivec) === ivec
    @test issorted(ivec)
    data = rand(0x00:0x03, n)
    ivec = IntVector{2}(data)
    @test issorted(radixsort(ivec))
    @test radixsort!(ivec) === ivec
    @test issorted(ivec)
    data = rand(0x00:0x07, n)
    ivec = IntVector{3}(data)
    @test issorted(radixsort(ivec))
    @test radixsort!(ivec) === ivec
    @test issorted(ivec)
end

# thorough and time-consuming tests for each combination of width and element type

@testset "conversion" begin
    @testset "empty" begin
        for T in Ts, w in 1:sizeof(T)*8
            data = Vector{T}(undef, 0)
            @test typeof(IntArray{w}(data)) == IntArray{w,T,1}
            @test length(IntArray{w}(data)) == 0

            data = Matrix{T}(undef, 0, 0)
            @test typeof(IntArray{w}(data)) == IntArray{w,T,2}
            @test size(IntArray{w}(data)) == (0, 0)
        end
    end
    @testset "small" begin
        for T in Ts, w in 1:sizeof(T)*8
            data = T[0x00, 0x01]
            @test typeof(IntArray{w}(data)) == IntArray{w,T,1}
            @test length(IntArray{w}(data)) == 2

            data = T[0x00 0x01 0x00; 0x01 0x00 0x01]
            @test typeof(IntArray{w}(data)) == IntArray{w,T,2}
            @test size(IntArray{w}(data)) == (2, 3)
        end
    end
    @testset "large" begin
        n = 1000
        for T in Ts, w in 1:sizeof(T)*8
            data = rand(T, n)
            @test typeof(IntArray{w}(data)) == IntArray{w,T,1}
            @test length(IntArray{w}(data)) == n
        end
    end
end

@testset "random run" begin
    @testset "IntVector" begin
        n = 123
        for T in Ts, w in 1:sizeof(T)*8
            data = rand(T(0):T(2)^w-T(1), n)
            ivec = IntVector{w,T}(data)
            for i in 1:n
                @test ivec[i] == data[i]
            end

            # random update
            for _ in 1:100
                i = rand(1:n)
                x::T = rand(T) % w
                data[i] = x
                ivec[i] = x
                @test ivec[i] == data[i]
            end
            for i in 1:n
                @test ivec[i] == data[i]
            end
            @test_throws BoundsError ivec[0]
            @test_throws BoundsError ivec[n+1]

            # sort
            sort!(data)
            radixsort!(ivec)
            for i in 1:n
                @test ivec[i] == data[i]
            end

            # random push!/pop!
            for _ in 1:100
                if rand() < 0.5
                    x::T = rand(T) % w
                    push!(ivec, x)
                    push!(data, x)
                else
                    pop!(ivec)
                    pop!(data)
                end
                @test length(ivec) == length(data)
                @test ivec[end] == data[end]
            end
            while !isempty(ivec)
                x = pop!(ivec)
                y = pop!(data)
                @test x == y
            end
            @test isempty(ivec)
        end
    end
    @testset "IntMatrix" begin
        m, n = 11, 28
        for T in Ts, w in 1:sizeof(T)*8
            data = rand(T(0):T(2)^w-T(1), (m, n))
            imat = IntMatrix{w,T}(data)
            for i in 1:m, j in 1:n
                @test imat[i,j] == data[i,j]
            end

            # random update
            for _ in 1:100
                i = rand(1:m)
                j = rand(1:n)
                x::T = rand(T) % w
                data[i,j] = x
                imat[i,j] = x
                @test imat[i,j] == data[i,j]
            end
            for i in 1:m, j in 1:n
                @test imat[i,j] == data[i,j]
            end
        end
    end
end


text:
using SAIMFitKit

#-------------- OPTICAL PARAMETERS----------------------
opt = SAIMOptics()
opt.nB = 1.33;	                    #The refractive index of the ambient buffer / cytoplasm
opt.dOx = 1900.1;	                #The thickness of the SiO2 layer in units of nm

#First excitation laser
opt.λ_Ex_1 = 488.0; 		        #The wavelength of excitation laser #1 in units of nm
opt.nOx_1 = 1.4775;		            #The refractive index of SiO2 at excitation wavelength #1
opt.nSi_1 = 4.3707 + 0.080068im; 	#The complex refractive index of Si at excitation wavelength #1; #4.3 + 0.073im; 0.080068

#Second excitation laser
opt.λ_Ex_2 = 642.0;		        	#The wavelength of the microscope excitation laser #1 in units of nm
opt.nOx_2 = 1.4719;		            #The refractive index of SiO2
opt.nSi_2 = 3.8660 + 0.017933im;    #The complex refractive index of Si

#-------------- FIT PARAMETERS----------------------
p0 = [0.5, 1.0, 0.5, 1.0, 75.0]		            #Initial guesses for parameters [A1, B1, A2, B2, H]
lb = [0.0, 0.0, 0.0, 0.0, 0.0]					#Lower bounds for [A1, B1, A2, B2, H]
ub = [10000., 10000., 10000., 10000., 450.]		#Upper bounds for [A1, B1, A2, B2, H]

#-------------- IMAGE STACKS ---------------------
angles = range(5.0, length=32, stop=43.75)	#Incidece angle in degrees for each image frame
files = ["488_noisy_2_647_noisy_c"]				#Names of image stacks; multiple experiments can be listed to run in batch
path = "C:\\Users\\matth\\.julia\\dev\\SAIMFitKit\\example\\TestImages" #Path to image files

#-------------- CONDUCT FITS ---------------------
for f in files
    fit_SAIM(path, f, opt, angles, p0, lb, ub, show=true, glb=false, color=2)
end


text:
# Custom array type with nontraditional indices
# (we could use OffsetArrays, but that will have it's own `showarg` methods
#  and here we want to test fallbacks)
module ZArrays
using CustomUnitRanges: filename_for_zerorange
include(filename_for_zerorange)

export ZArray

struct ZArray{T,N} <: AbstractArray{T,N}
    parent::Array{T,N}
end
Base.indices(a::ZArray) = map(ZeroRange, size(a.parent))
end


text:
#
# @description Implements the 3DES (DES-EDE1, DES-EDE2, and DES-EDE3) block cipher algorithms.
#
# @author Alejandro Claro (alejandro.claro@gmail.com)
#
# Copyright 2017 All rights reserved.
# Use of this source code is governed by a MIT-style license that can be found in the LICENSE file.
#
export TripleDesCipher
export TRIPLE_DES_BLOCK_SIZE, TRIPLE_DES_LEGAL_KEY_SIZES
export block_size, key_size, encrypt, decrypt

# @description Defines the 3DES cipher data struture.
struct TripleDesCipher <: BlockCipher
  des::Array{DesCipher}
  key_size::UInt64

  # @description Initialize the 3DES cipher.
  #
  # @param {Vector{UInt8}} key The symmetric key.
  function TripleDesCipher(key::Vector{UInt8})
    des_keys = []

    if length(key) == TRIPLE_DES_LEGAL_KEY_SIZES[1]
      des_keys = [ key, key, key ]
    elseif length(key) == TRIPLE_DES_LEGAL_KEY_SIZES[2]
      des_keys = [ key[1 : 8], key[9 : 16], key[1 : 8] ]
    elseif length(key) == TRIPLE_DES_LEGAL_KEY_SIZES[3]
      des_keys = [ key[1 : 8], key[9 : 16], key[17 : 24] ]
    else
      throw(ArgumentError("Invalid key length. $(length(key) * 8) bits provided."))
    end

    return new([ DesCipher(des_keys[1]), DesCipher(des_keys[2]), DesCipher(des_keys[3]) ], length(key))
  end
end

# @description The 3DES cipher block size in bytes.
const TRIPLE_DES_BLOCK_SIZE = UInt64(8)

# @description The 3DES cipher legal key sizes in bytes.
const TRIPLE_DES_LEGAL_KEY_SIZES = UInt64[ 8, 16, 24 ]

# @description Gets the cipher block size.
#
# @param {TripleDesCipher} self The cipher data struture.
#
# @return {Integer} The block size in bytes.
function block_size(self::TripleDesCipher)
  return TRIPLE_DES_BLOCK_SIZE
end

# @description Gets the cipher key size.
#
# @param {TripleDesCipher} self The cipher data struture.
#
# @return {Integer} The key size in bytes.
function key_size(self::TripleDesCipher)
  return self.key_size
end

# @description Encrypts the given plaintext block with the key set at initialization.
#
# @param {TripleDesCipher} self      The cipher data struture.
# @param {Vector{UInt8}}   plaintext The message to encrypt.
#
# @return {Vector{UInt8}} The resulting ciphertext block.
function encrypt(self::TripleDesCipher, plaintext::Vector{UInt8})
  result = encrypt(self.des[1], plaintext)
  result = decrypt(self.des[2], result)
  result = encrypt(self.des[3], result)
  return result
end

# @description Decrypts the given ciphertext block with the key set at initialization.
#
# @param {TripleDesCipher} self       The cipher data struture.
# @param {Vector{UInt8}}   ciphertext The message to decrypt.
#
# @return {Vector{UInt8}} The resulting plaintext block.
function decrypt(self::TripleDesCipher, ciphertext::Vector{UInt8})
  result = decrypt(self.des[3], ciphertext)
  result = encrypt(self.des[2], result)
  result = decrypt(self.des[1], result)
  return result
end


text:
module test_julia_views

using Test

A = [1 2 3; 4 5 6]
first_row = @views A[1, :]
@test first_row isa SubArray
@test eltype(A) === eltype(first_row) === Int
@test ndims(first_row) == 1
@test first_row[3] == 3

first_row[3] = 0
@test A == [1 2 0; 4 5 6]

end # module test_julia_views


text:
module StreamModels

using
    StatsBase,
    DataStreams,
    CategoricalArrays,
    ArgCheck,
    Compat

using CategoricalArrays: CategoricalPool

export
    ModelBuilder,
    Terms,
    Formula,
    @formula,
    modelmatrix,
    @model,
    build,

    AbstractContrasts,
    EffectsCoding,
    DummyCoding,
    HelmertCoding,
    ContrastsCoding

function name end

const DEBUG = true

macro debug(msg)
    DEBUG ? :(println(string($(esc(msg))))) : nothing
end

include("typedefs.jl")
include("contrasts.jl")
include("utils.jl")
include("terms.jl")
include("formula.jl")
include("summarize/summarize.jl")
include("set_schema.jl")
include("modelmatrix.jl")
include("model.jl")

end # module


text:
abstract type HTTPAdapterBase <: ProtocolAdapter end

function loop_body(adp::HTTPAdapterBase, qq, poll_interval, fetch_count, response_type_func, session_key_position)
    cmd = Commands.fetchLatestMessage(fetch_count)
    while !adp.closed
        data = nothing
        try
            data = send(adp, cmd; session_key_position, log=false)::response_type_func(cmd)
            for msg in data.data
                put!(adp.output_channel, msg)
            end
        catch e
            if isnothing(data)
                put!(adp.output_channel, ExceptionAndBacktrace(e, catch_backtrace()))
            else
                put!(adp.output_channel, dat)
            end
        end
        sleep(poll_interval)
    end
    close(adp.output_channel)
    resp = post_restful("http://$(adp.server)/release", (), (; adp.sessionKey, qq))
    @info "HTTP adapter quitted"
end

struct EmbedSessionKey{T}
    sessionKey::String
    data::T
end

StructTypes.StructType(::Type{<:EmbedSessionKey}) = StructTypes.DictType()

@generated function StructTypes.keyvaluepairs(x::EmbedSessionKey{T}) where T
    fields = (:(x.data.$f) for f in fieldnames(T))
    quote
        pairs((; x.sessionKey, $(fields...)))
    end
end

@enum SessionKeyPosition begin
    SESSION_KEY_IN_HEADERS
    SESSION_KEY_IN_BODY
end

command_to_path(cmd::GeneralCommand) = replace(string(cmd.command), '_' => '/')

function send(adp::HTTPAdapterBase, cmd::GeneralCommand; session_key_position::SessionKeyPosition)
    url = "http://$(adp.server)/$(command_to_path(cmd))"
    if session_key_position == SESSION_KEY_IN_HEADERS
        headers = ["sessionKey" => adp.sessionKey]
    else
        headers = []
    end
    if session_key_position == SESSION_KEY_IN_BODY
        body = EmbedSessionKey(adp.sessionKey, cmd.content)
    else
        body = cmd.content
    end
    resp = if cmd.method == Commands.CommandMethods.GET
        get_restful(url, headers, body)
    elseif cmd.method == Commands.CommandMethods.POST
        post_restful(url, headers, body)
    elseif cmd.method == Commands.CommandMethods.UPLOAD
        upload(url, headers, body)
    end
    resp
end

function Base.close(adp::HTTPAdapterBase)
    adp.closed = true
end

function post_restful(url, headers, data)
    body = JSON3.write(data)
    json_header = "Content-Type" => "application/json"
    headers = [json_header, headers...]
    @debug "body = $body"
    r = HTTP.post(url, headers, body)
    data = JSON3.read(r.body)
    @debug "response: $(String(r.body))"
    return data
end

function get_restful(url, headers, query)
    query = StructTypes.constructfrom(Dict, query)
    map!(string, values(query))
    @debug "query = $query"
    r = HTTP.get(url, headers; query)
    data = JSON3.read(r.body)
    @debug "response: $(String(r.body))"
    return data
end

function upload(url, headers, form)
    form = StructTypes.constructfrom(Dict, form)
    map!(values(form)) do v
        v isa Union{Integer,Enum} ? string(v) : v
    end
    form = HTTP.Form(form)
    @debug "form = $form"
    r = HTTP.post(url, headers, form)
    ret = JSON3.read(r.body)
    @debug "response: $(String(r.body))"
    ret
end


text:
function interpolate(obj::FieldTimeInterpolable, arg0::FieldAbsoluteDate, arg1::Collection)
    return jcall(obj, "interpolate", FieldTimeInterpolable, (FieldAbsoluteDate, Collection), arg0, arg1)
end

function interpolate(obj::FieldTimeInterpolable, arg0::FieldAbsoluteDate, arg1::Stream)
    return jcall(obj, "interpolate", FieldTimeInterpolable, (FieldAbsoluteDate, Stream), arg0, arg1)
end


text:
### Transition distribution
using GridInterpolations

function POMDPs.transition(pomdp::SingleOCPOMDP, s::SingleOCState, a::SingleOCAction, dt::Float64 = pomdp.ΔT)
    if s.crash || s.ego.posG.x >= pomdp.x_goal || isterminal(pomdp, s)
        return Deterministic(TERMINAL_STATE)
    end
    ## Find Ego states first
    ego_states, ego_probs = ego_transition(pomdp, s.ego, a, dt)

    ### Find pedestrian states
    ped_probs, ped_states = ped_transition(pomdp, s.ped, dt)

    n_next_states = length(ego_states)*length(ped_states)

    ### Total state
    next_states = Vector{SingleOCState}(undef, n_next_states)
    next_probs = zeros(n_next_states)
    ind = 1
    for (i,ego) in enumerate(ego_states)
        for (j,ped) in enumerate(ped_states)
            collision = collision_checker(ego, ped, pomdp.ego_type, pomdp.ped_type)
            next_states[ind] = SingleOCState(collision, ego, ped)
            next_probs[ind] = ped_probs[j]*ego_probs[i]
            ind += 1
        end
    end
    normalize!(next_probs, 1)
    return SparseCat(next_states, next_probs)
end

"""
    ego_transition(pomdp::SingleOCPOMDP, ego::VehicleState, a::SingleOCAction)
Returns the distribution over the possible future state for the ego car only
"""
function ego_transition(pomdp::SingleOCPOMDP, ego::VehicleState, a::SingleOCAction, dt::Float64 = pomdp.ΔT)
    x_ = ego.posG.x + ego.v*dt + 0.5*a.acc*dt^2
    if x_ <= ego.posG.x # no backup
        x_ = ego.posG.x
    end
    v_ = ego.v + a.acc*dt # no backup
    if v_ <= 0.
        v_ = 0.
    end

    grid = RectangleGrid(get_X_grid(pomdp), get_V_grid(pomdp)) #XXX must not allocate the grid at each function call, find better implementation
    index, weight = interpolants(grid, [x_, v_])
    n_pts = length(index)

    states = Array{VehicleState}(undef, n_pts)
    probs = Array{Float64}(undef, n_pts)
    for i=1:n_pts
        xg, vg = ind2x(grid, index[i])
        states[i] = xv_to_state(pomdp, xg, vg)
        probs[i] = weight[i]
    end
    return states, probs
end

"""
    ped_transition(pomdp::SingleOCPOMDP, ped::VehicleState)
Return the distribution over the possible future state for the pedestrian only
"""
function ped_transition(pomdp::SingleOCPOMDP, ped::VehicleState, dt::Float64 = pomdp.ΔT)
    states = VehicleState[]
    probs = Float64[]
    sizehint!(states, 8)
    env = pomdp.env

    if pomdp.no_ped
        return [1.0], [get_off_the_grid(pomdp)]
    end

    if off_the_grid(pomdp, ped) # appear with random speed or stay of the grid
        p_birth = pomdp.p_birth
        V_grid = LinRange(0., env.params.ped_max_speed, Int(floor(env.params.ped_max_speed/pomdp.vel_res)) + 1)
        Y_grid = get_Y_grid(pomdp)
        for v in V_grid
            for y in Y_grid[1:div(length(Y_grid),2)- 3]
                push!(states, yv_to_state(pomdp, y, v))
            end
        end
        probs = ones(length(states) + 1)
        probs[1:end - 1] .= p_birth/length(states)
        # add the off the grid state
        push!(states, get_off_the_grid(pomdp))
        probs[end] = 1.0 - p_birth
        normalize!(probs, 1)
        return probs, states
    end

    grid = RectangleGrid(get_Y_grid(pomdp), get_V_ped_grid(pomdp)) #XXX preallocate
    y_ = ped.posG.y + ped.v*dt
    if y_ > pomdp.y_goal
        return [1.0], [get_off_the_grid(pomdp)]
    end
    for v_noise in [-pomdp.v_noise, 0., pomdp.v_noise]
        v_ = 1.0 + v_noise
        ind, weight = interpolants(grid, [y_, v_])
        for i=length(ind)
            yg, vg = ind2x(grid, ind[i])
            state = yv_to_state(pomdp, yg, vg)
            if !(state in states) # check for doublons
                push!(states, state)
                push!(probs, weight[i])
            else
                state_ind = findfirst(isequal(state), states)
                probs[state_ind] += weight[i]
            end
        end
    end
    # add roughening
    normalize!(probs, 1)
    probs .+= maximum(probs)
    normalize!(probs, 1)
    return probs, states
end

###### HELPERS #####################################################################################

"""
Helper to generate ego state
"""
function xv_to_state(pomdp::SingleOCPOMDP, x::Float64, v::Float64)
    return VehicleState(VecSE2(x, 0., 0.), pomdp.env.roadway, v)
end

"""
Helper to generate pedestrian state
"""
function yv_to_state(pomdp::SingleOCPOMDP, y::Float64, v::Float64)
    x_ped = 0.5*(pomdp.env.params.roadway_length - pomdp.env.params.crosswalk_width + 1)
    return VehicleState(VecSE2(x_ped, y, pi/2),
                        pomdp.env.roadway, v)
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

# These utility functions provide a simple interface to the Emitters package to create emitters that are commonly used in optical systems. Many optical systems by convention have their optical axis parallel to the z axis. These emitters direct the rays in the negative z direction, toward the entrance of the optical system.

"""
    pointemitter(origin::AbstractVector{T}, coneangle; λ::Length = 500nm, numrays = 100) where {T<:Real}

Creates a point source with Lambertian emission power and cone distribution of rays, emitting in the -z direction. λ is
a unitful Length quantity, e.g., 550nm.
"""
function pointemitter(origin::AbstractVector{T}, coneangle; λ::Length = 500nm, numrays = 100) where {T<:Real}
    @assert length(origin) == 3

    return Sources.Source(
        Geometry.identitytransform(T),
        Spectrum.DeltaFunction(λ),
        Origins.Point(Geometry.Vec3(origin)),
        Directions.UniformCone(-Geometry.unitZ3(T), coneangle, numrays),
        AngularPower.Lambertian()
    )
end

"""
    collimatedemitter(origin::AbstractVector{T}, halfsquaresize; λ::Length = 500nm, numrays = 100) where {T<:Real}

Creates a square collimated emitter, emitting rays in the -z direction. Rays are emitted on a square grid with
sqrt(numrays) on a side. λ can be a unitful quantity, e.g., 550nm, or a number. In the latter case the units are
implicitly microns.
"""
function collimatedemitter(origin::AbstractVector{T}, halfsquaresize; λ::Length = 500nm, numrays = 100) where {T<:Real}
    samples = Int(round(sqrt(numrays)))

    return Sources.Source(
        Geometry.translation(origin...),
        Spectrum.DeltaFunction(λ),
        Origins.RectGrid(2 * halfsquaresize, 2 * halfsquaresize, samples, samples),
        Directions.Constant(-Geometry.unitZ3(T)),
        AngularPower.Lambertian()
    )
end


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

import Test
Test.@testset "SIS3316Digitizers" begin

end # testset


text:
#unit-unum.jl

#test that the Unum pseudo-constructor bifurcates the actual construction to
#"UnumSmall" and "UnumLarge" types
small_typetest = Unum{4,6}(z64, z64, z16, z16, z16)
large_typetest = Unum{4,7}(z64, [z64, z64], z16, z16, z16)

@test isa(small_typetest, Unum{4,6})
@test isa(small_typetest, Unums.UnumSmall{4,6})
@test isa(large_typetest, Unum{4,7})
@test isa(large_typetest, Unums.UnumLarge{4,7})

#test that passing through the Unum copy constructor is valid.
small_typetest_prime = Unum{4,6}(small_typetest)
large_typetest_prime = Unum{4,7}(large_typetest)
@test isa(small_typetest_prime, Unum{4,6})
@test isa(small_typetest_prime, Unums.UnumSmall{4,6})
@test isa(large_typetest_prime, Unum{4,7})
@test isa(large_typetest_prime, Unums.UnumLarge{4,7})

#test that passing strange data types fails.
@test_throws MethodError Unum{4,6}(z16, [z64, z64], z16, z16, z16)
@test_throws MethodError Unum{4,7}(z16, z64, z16, z16, z16)
@test_throws MethodError Unum{4,8}(z16, [z64, z64], z16, z16, z16)

#cascading tests on the general unum checking procedure.
@test_throws ArgumentError Unums.__general_unum_check(7, 0, z64, z16, z16)          #ESS too big.
@test_throws ArgumentError Unums.__general_unum_check(6, 12, z64, z16, z16)         #FSS too big.
@test_throws ArgumentError Unums.__general_unum_check(0, 0, z64, z16, o16)          #fsize too big
@test_throws ArgumentError Unums.__general_unum_check(0, 11, z64, z16, UInt16(2048)) #fsize too big
@test_throws ArgumentError Unums.__general_unum_check(0, 0, z64, o16, z16)          #esize too big
@test_throws ArgumentError Unums.__general_unum_check(6, 0, z64, UInt16(64), z16)   #esize too big
@test_throws ArgumentError Unums.__general_unum_check(0, 0, UInt64(2), z16, z16)    #exponent too big

#actual tests on constructors using the safe unum constructor.
#test that fsize error gets thrown.
@test_throws ArgumentError Unum{0,0}(z64, z64, z16, z16, o16)
#test that esize error gets thrown
@test_throws ArgumentError Unum{0,0}(z64, z64, z16, o16, z16)
#test that exponent error gets thrown
@test_throws ArgumentError Unum{0,0}(UInt64(2), z64, z16, z16, z16)
#test that the fraction error gets thrown
@test_throws ArgumentError Unum{0,0}(z64, [z64, z64], z16, z16, z16)
@test_throws ArgumentError Unum{0,8}(z64, z64, z16, z16, z16)

#make sure that the constructor doesn't trim fractions that aren't too long.
warlpiri_one = Unum{0,0}(z64, t64, z16, z16, z16)
@test warlpiri_one.fraction == t64
@test warlpiri_one.fsize == z16
@test warlpiri_one.flags == z16

#test that the constructor correctly trims fractions that are too long.
warlpiri_some = Unum{0,0}(z64, f64, z16, z16, z16)
@test warlpiri_some.fraction == t64
@test warlpiri_some.fsize == z16
@test warlpiri_some.flags & Unums.UNUM_UBIT_MASK == Unums.UNUM_UBIT_MASK

#test the same thing in a really big number.
bigunum_trim1 = Unum{4,7}(z64, [z64, t64], z16, z16, UInt16(63))
@test bigunum_trim1.fraction.a == [z64, z64]
@test bigunum_trim1.fsize == 63
@test bigunum_trim1.flags & Unums.UNUM_UBIT_MASK == Unums.UNUM_UBIT_MASK

#test the unum_easy constructor
easy_warlpiri_two = unum(Unum{0,0}, 1, z64, z16)
@test easy_warlpiri_two.exponent == o64
@test easy_warlpiri_two.fraction == z64
@test easy_warlpiri_two.fsize == z16
@test easy_warlpiri_two.flags == z16

@devmode_on begin
  #test to see that these unsafe constructors fail.
  @test_throws ArgumentError Unums.UnumSmall{0,0}(      z64, z64, z16, z16, o16)
  @test_throws ArgumentError Unums.UnumSmall{0,0}(      z64, z64, z16, o16, z16)
  @test_throws ArgumentError Unums.UnumSmall{0,0}(UInt64(2), z64, z16, z16, o16)
  #we can even create an unsafe type
  @test_throws ArgumentError Unums.UnumSmall{4,8}(z64, z64, z16, z16, z16)

  #repeat the test in a high unum environment.
  @test_throws ArgumentError Unums.UnumLarge{4,8}(      z64, Unums.ArrayNum{8}([z64, z64, z64, z64]), z16, z16, UInt16(256))
  @test_throws ArgumentError Unums.UnumLarge{4,8}(      z64, Unums.ArrayNum{8}([z64, z64, z64, z64]), z16, UInt16(16), z16)
  @test_throws ArgumentError Unums.UnumLarge{4,8}(UInt64(2), Unums.ArrayNum{8}([z64, z64, z64, z64]), z16, z16, z16)
end

@devmode_off begin
  #unsafe unums can be created using the unsafe constructor.
  unsafe_fsize = Unums.UnumSmall{0,0}(z64, z64, z16, z16, o16)
  @test unsafe_fsize.fsize == o16
  unsafe_esize = Unums.UnumSmall{0,0}(z64, z64, z16, o16, z16)
  @test unsafe_esize.esize == o16
  unsafe_exponent = Unums.UnumSmall{0,0}(UInt64(2), z64, z16, o16, z16)
  @test unsafe_exponent.exponent == UInt64(2)
  #we can even create an unsafe type
  unsafe_type = Unums.UnumSmall{4,8}(z64, z64, z16, z16, z16)

  #then show that these trigger TypeErrors when passed through the safe constructor
  @test_throws ArgumentError unum(unsafe_fsize)
  @test_throws ArgumentError unum(unsafe_esize)
  @test_throws ArgumentError unum(unsafe_exponent)
  @test_throws MethodError unum(unsafe_type)

  #repeat the test in a high unum environment.
  unsafe_fsize = Unums.UnumLarge{4,8}(z64, Unums.ArrayNum{8}([z64, z64, z64, z64]), z16, z16, UInt16(256))
  @test unsafe_fsize.fsize == UInt16(256)
  unsafe_esize = Unums.UnumLarge{4,8}(z64, Unums.ArrayNum{8}([z64, z64, z64, z64]), z16, UInt16(16), z16)
  @test unsafe_esize.esize == UInt16(16)
  unsafe_fraction = Unums.UnumLarge{4,8}(z64, Unums.ArrayNum{8}([z64, z64]), z16, z16, z16)
  @test unsafe_fraction.fraction.a == [z64, z64]
  unsafe_exponent = Unums.UnumLarge{4,8}(UInt64(2), Unums.ArrayNum{8}([z64, z64, z64, z64]), z16, z16, z16)
  @test unsafe_exponent.exponent == UInt16(2)
  unsafe_type = Unums.UnumLarge{0,0}(z64, Unums.ArrayNum{0}([z64, z64, z64, z64]), z16, z16, z16)

  @test_throws ArgumentError unum(unsafe_fsize)
  @test_throws ArgumentError unum(unsafe_esize)
  @test_throws MethodError   unum(unsafe_fraction)
  @test_throws ArgumentError unum(unsafe_exponent)
  @test_throws MethodError   unum(unsafe_type)
end

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
@testset verbose=true "map2salm" begin
    # NOTE: Float16 irfft returns Float32, which leads to type conflicts, so we just don't support
    # map2salm on Float16

    function sYlm(s::Int, ell::Int, m::Int, theta::T, phi::T) where {T<:Real}
        # Eqs. (II.7) and (II.8) of https://arxiv.org/abs/0709.0093v3
        # Note their weird definition w.r.t. `-s`
        k_min = max(0, m + s)
        k_max = min(ell + m, ell + s)
        sin_half_theta, cos_half_theta = sincos(theta / 2)
        return (-1)^(-s) * sqrt((2 * ell + 1) / (4 * T(π))) *
            T(sum(
                (-1) ^ (k)
                * sqrt(factorial(big(ell + m)) * factorial(big(ell - m)) * factorial(big(ell - s)) * factorial(big(ell + s)))
                * (cos_half_theta ^ (2 * ell + m + s - 2 * k))
                * (sin_half_theta ^ (2 * k - s - m))
                / (factorial(big(ell + m - k)) * factorial(big(ell + s - k)) * factorial(big(k)) * factorial(big(k - s - m)))
                for k in k_min:k_max
            )) *
            cis(m * phi)
    end

    m_m2Y2m = [
        (2, NINJA.m2Y22),
        (1, NINJA.m2Y21),
        (0, NINJA.m2Y20),
        (-1, NINJA.m2Y2m1),
        (-2, NINJA.m2Y2m2)
    ]

    @testset "Input expressions $T" for T in [BigFloat, Float64, Float32]
        # These are just internal consistency tests of the sYlm function
        # above, against the explicit expressions `mY2.`
        s = -2
        ℓ = 2
        Nϑ = 17
        Nφ = 18
        for (m, m2Y2m) in m_m2Y2m
            f1 = mapslices(ϕθ -> sYlm(s, ℓ, m, ϕθ[2], ϕθ[1]), phi_theta(Nφ, Nϑ, T), dims=[3])
            f2 = mapslices(ϕθ -> m2Y2m(ϕθ[2], ϕθ[1]), phi_theta(Nφ, Nϑ, T), dims=[3])
            @test f1 ≈ f2 atol=10eps(T) rtol=10eps(T)
        end
    end

    @testset "map2salm $T" for T in [BigFloat, Float64, Float32]
        # These test the ability of map2salm to precisely decompose the results of `sYlm`.
        ℓmax = 7
        Nϑ = 2ℓmax + 1
        Nφ = 2ℓmax + 2
        for s in -2:2
            #for ℓmin in 0:abs(s)
            let ℓmin = 0
                for ℓ in abs(s):ℓmax
                    for m in -ℓ:ℓ
                        f = mapslices(
                            ϕθ -> sYlm(s, ℓ, m, ϕθ[2], ϕθ[1]),
                            phi_theta(Nφ, Nϑ, T),
                            dims=[3]
                        )
                        computed = map2salm(f, s, ℓmax)
                        expected = zeros(Complex{T}, size(computed))
                        expected[SphericalFunctions.Yindex(ℓ, m, ℓmin)] = one(T)
                        if ≉(computed, expected, atol=30eps(T), rtol=30eps(T))
                            @show T
                            @show ℓmax
                            @show Nϑ
                            @show Nφ
                            @show s
                            @show ℓmin
                            @show ℓ
                            @show m
                            println("computed = $computed")
                            println("expected = $expected")
                            println("max_diff = ", maximum(abs, computed .- expected), ";")
                            println()
                        end
                        @test computed ≈ expected atol=30eps(T) rtol=30eps(T)

                        plan = plan_map2salm(f, s, ℓmax)
                        computed2 = map2salm(f, plan)
                        @test array_equal(computed, computed2)
                    end
                end
            end
        end
    end

end


text:
module LibBZ2

include("lowlevel.jl")

end # module