text:
using SignatureGB
using Singular
P, (x1, x2, x3, x4, x5, x6, x7, x8, x9)  = Singular.PolynomialRing(Fp(65521), ["x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8", "x9"]);
gb = SignatureGB.f5([x1, x2, x3, x4, x5, x6, x7, x8, x9]);F = [729*x1^2*x2^2-2970*x1^2*x2*x5+5426*x1^2*x5^2+7056*x1^2*x6^2-1944*x1*x2^3+3960*x1*x2^2*x5-4104*x1*x2*x3*x9-1512*x1*x2*x4*x5-4368*x1*x2*x6^2+196*x1*x3*x5*x7+8360*x1*x3*x5*x9+3080*x1*x4*x5^2-12096*x1*x4*x6*x7+8568*x1*x5*x6^2+1960*x1*x5*x7^2-4018*x1*x5*x7*x9-4214*x1*x5*x8*x9+13608*x1*x6*x8*x9+1296*x2^4+5472*x2^2*x3*x9+2016*x2^2*x4*x5+676*x2^2*x6^2+3744*x2*x4*x6*x7-2652*x2*x5*x6^2-4212*x2*x6*x8*x9+4*x3^2*x7^2+5776*x3^2*x9^2+4256*x3*x4*x5*x9+80*x3*x7^3-164*x3*x7^2*x9-172*x3*x7*x8*x9+784*x4^2*x5^2+5184*x4^2*x7^2-7344*x4*x5*x6*x7-11664*x4*x7*x8*x9+2601*x5^2*x6^2+8262*x5*x6*x8*x9+400*x7^4-1640*x7^3*x9-1720*x7^2*x8*x9+1681*x7^2*x9^2+3526*x7*x8*x9^2+8410*x8^2*x9^2-98*x1*x2*x5+12264*x1*x5*x6-4*x2*x3*x7-3796*x2*x5*x6-40*x2*x7^2+82*x2*x7*x9+86*x2*x8*x9-10512*x4*x5*x7+7446*x5^2*x6+11826*x5*x8*x9+540*x1*x2-1100*x1*x5-719*x2^2-1520*x3*x9-560*x4*x5+5329*x5^2+100, 1458*x1^2*x2-2970*x1^2*x5-5832*x1*x2^2+7920*x1*x2*x5-4104*x1*x3*x9-1512*x1*x4*x5-4368*x1*x6^2+5184*x2^3+10944*x2*x3*x9+4032*x2*x4*x5+1352*x2*x6^2+3744*x4*x6*x7-2652*x5*x6^2-4212*x6*x8*x9-98*x1*x5-4*x3*x7-3796*x5*x6-40*x7^2+82*x7*x9+86*x8*x9+540*x1-1438*x2, -4104*x1*x2*x9+196*x1*x5*x7+8360*x1*x5*x9+5472*x2^2*x9+8*x3*x7^2+11552*x3*x9^2+4256*x4*x5*x9+80*x7^3-164*x7^2*x9-172*x7*x8*x9-4*x2*x7-1520*x9, -1512*x1*x2*x5+3080*x1*x5^2-12096*x1*x6*x7+2016*x2^2*x5+3744*x2*x6*x7+4256*x3*x5*x9+1568*x4*x5^2+10368*x4*x7^2-7344*x5*x6*x7-11664*x7*x8*x9-10512*x5*x7-560*x5, -2970*x1^2*x2+10852*x1^2*x5+3960*x1*x2^2-1512*x1*x2*x4+196*x1*x3*x7+8360*x1*x3*x9+6160*x1*x4*x5+8568*x1*x6^2+1960*x1*x7^2-4018*x1*x7*x9-4214*x1*x8*x9+2016*x2^2*x4-2652*x2*x6^2+4256*x3*x4*x9+1568*x4^2*x5-7344*x4*x6*x7+5202*x5*x6^2+8262*x6*x8*x9-98*x1*x2+12264*x1*x6-3796*x2*x6-10512*x4*x7+14892*x5*x6+11826*x8*x9-1100*x1-560*x4+10658*x5, 14112*x1^2*x6-8736*x1*x2*x6-12096*x1*x4*x7+17136*x1*x5*x6+13608*x1*x8*x9+1352*x2^2*x6+3744*x2*x4*x7-5304*x2*x5*x6-4212*x2*x8*x9-7344*x4*x5*x7+5202*x5^2*x6+8262*x5*x8*x9+12264*x1*x5-3796*x2*x5+7446*x5^2, 196*x1*x3*x5-12096*x1*x4*x6+3920*x1*x5*x7-4018*x1*x5*x9+3744*x2*x4*x6+8*x3^2*x7+240*x3*x7^2-328*x3*x7*x9-172*x3*x8*x9+10368*x4^2*x7-7344*x4*x5*x6-11664*x4*x8*x9+1600*x7^3-4920*x7^2*x9-3440*x7*x8*x9+3362*x7*x9^2+3526*x8*x9^2-4*x2*x3-80*x2*x7+82*x2*x9-10512*x4*x5, -4214*x1*x5*x9+13608*x1*x6*x9-4212*x2*x6*x9-172*x3*x7*x9-11664*x4*x7*x9+8262*x5*x6*x9-1720*x7^2*x9+3526*x7*x9^2+16820*x8*x9^2+86*x2*x9+11826*x5*x9, -4104*x1*x2*x3+8360*x1*x3*x5-4018*x1*x5*x7-4214*x1*x5*x8+13608*x1*x6*x8+5472*x2^2*x3-4212*x2*x6*x8+11552*x3^2*x9+4256*x3*x4*x5-164*x3*x7^2-172*x3*x7*x8-11664*x4*x7*x8+8262*x5*x6*x8-1640*x7^3-1720*x7^2*x8+3362*x7^2*x9+7052*x7*x8*x9+16820*x8^2*x9+82*x2*x7+86*x2*x8+11826*x5*x8-1520*x3];
sort!(F, by = p -> leading_monomial(p))
println(@elapsed SignatureGB.f5(F, verbose=true));exit()


text:
using QuantumCumulants
using OrdinaryDiffEq
using SteadyStateDiffEq
using ModelingToolkit
using Test

@testset "diffeq" begin

hf = FockSpace(:cavity)
ha = NLevelSpace(:atom,(:g,:e))
h = tensor(hf, ha)

@qnumbers a::Destroy(h) σ::Transition(h,:g,:e)

# Single-atom laser
@cnumbers Δ g κ γ ν

H = Δ*a'*a + g*(a'*σ + σ'*a)
J = [a,σ,σ']
he_avg = meanfield([a'*a,σ'*σ,a*σ'],H,J;rates=[κ,γ,ν])

he_exp = cumulant_expansion(he_avg,2)
@test isequal(he_exp.equations, meanfield([a'*a,σ'*σ,a*σ'],H,J;rates=[κ,γ,ν],order=2).equations)

ps = [Δ,g,κ,γ,ν]
missed = find_missing(he_exp)
@test !any(p ∈ Set(missed) for p=ps)

# Exploit phase invariance
subs = Dict([missed; QuantumCumulants._conj.(missed)] .=> 0)
he_nophase = substitute(he_exp, subs)
@test isempty(find_missing(he_nophase))

@named sys = ODESystem(he_nophase)

# Numerical solution
p0 = ps .=> [0.0,0.5,1.0,0.1,0.9]
u0 = zeros(ComplexF64,3)
tmax = 10.0

prob = ODEProblem(sys,u0,(0.0,tmax),p0,jac=true)
sol = solve(prob,RK4())
n = sol[average(a'*a)]
pe = getindex.(sol.u,2)

@test all(iszero.(imag.(n)))
@test all(iszero.(imag.(pe)))
@test all(real.(n) .>= 0.0)
@test all(1.0 .>= real.(pe) .>= 0.0)

prob_ss = SteadyStateProblem(prob)
sol_ss = solve(prob_ss, DynamicSS(Tsit5(), abstol=1e-5, reltol=1e-5), abstol=1e-6, reltol=1e-6)

@test isapprox(sol_ss[a'*a], sol[a'a][end], atol=1e-4)

end # testset

ext: jl


text:
#using PyPlot
using JuliaFEM.Testing
#using JuliaFEM.MaterialModels: stiffnessTensor, calculate_stress, State
#using JuliaFEM.MaterialModels: stiffnessTensorPlaneStress

function test_von_mises_3D_basic()

    steps = 1000
    strain_max = 0.003
    num_cycles = 3
    E = 200.0e3
    nu =  0.3
    ν = 0.3
    C = stiffnessTensor(E, ν)

    strain_tot = zeros(Float64, (steps, 6))
    strain_tot2 = zeros(Float64, (steps, 6))
    strain_tot3 = zeros(Float64, (steps, 6))

    # Adding only strain in x-axis and counting for the poisson effect
    strain_tot[:, 1] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps))
    strain_tot[:, 2] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)).*-ν
    strain_tot[:, 3] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)).*-ν
    strain_tot[:, 4] = strain_max / 10 * sin(2 * pi * linspace(0, num_cycles, steps))

    strain_last = zeros(Float64, (6))
    strain_p = zeros(Float64, (6))
    stress = zeros(Float64, (6, 1))
    stress_y =  200.0
    ss = Float64[]
    ee = Float64[]

    eig_stress = zeros(Float64, (3, 3))
    eig_vals = zeros(Float64, (steps, 3))

    function fill_tensor(a, b)
        a[1, 1] = b[1]
        a[2, 2] = b[2]
        a[3, 3] = b[3]

        a[1, 2] = b[6]
        a[1, 3] = b[5]
        a[2, 3] = b[4]

        a[2, 1] = b[6]
        a[3, 1] = b[5]
        a[3, 2] = b[4]
    end

    mat = State(C, stress_y, zeros(Float64, 6), zeros(Float64, 6))

    info("Starting calculation")
    tic()
    #=
    for i=1:steps
        strain_new = reshape(strain_tot[i, :, :], (6, 1))
        dstrain = strain_new - mat.strain
        calculate_stress!(dstrain, mat, Val{:vonMises})
        mat.strain += vec(dstrain)
        push!(ss, mat.stress[1])
        push!(ee, mat.strain[1])

        fill_tensor(eig_stress, mat.stress)
        eig_vals[i, :] = sort(eigvals(eig_stress))
    end
    =#
    stress = zeros(Float64, 6)
    strain = zeros(Float64, 6)
    for i=1:steps
        strain_new = reshape(strain_tot[i, :, :], (6, 1))
        dstrain = strain_new - strain
        calculate_stress!(dstrain, stress, C, stress_y, Val{:vonMises})
        strain = vec(strain_new)
        push!(ss, stress[1])
        push!(ee, strain[1])
        fill_tensor(eig_stress, stress)
        eig_vals[i, :] = sort(eigvals(eig_stress))
    end

    toc()
    # ================ Plotting =================== #
    n(θ, ϕ) = [sin(θ)*cos(ϕ)
               sin(θ)*sin(ϕ)
               cos(θ)]
    m(θ, ϕ, χ) = [-sin(ϕ)*cos(χ)-cos(θ)*cos(ϕ)*sin(χ)
                   cos(ϕ)*cos(χ)-cos(θ)*sin(ϕ)*sin(χ)
                   sin(θ)*sin(χ)]

    w = [sqrt(2/3) * 200 * m(54.735 * pi / 180, 45 * pi/180, x) for x=0:0.15:(2*pi+0.1)]
    base_vec =  [1 1 1] / sqrt(3)

    for i=-5:5
        tt = [w[x] + vec(base_vec) + 50 * i for x=1:length(w)]
        x = map(x->tt[x][1], collect(1:length(w)))
        y = map(x->tt[x][2], collect(1:length(w)))
        z = map(x->tt[x][3], collect(1:length(w)))
        plot3D(x, y, z, color="blue")
    end

    tt = [w[x] + vec(base_vec) + 50 * -5 for x=1:length(w)]
    x_start = map(x->tt[x][1], collect(1:length(w)))[1:5:end]
    y_start = map(x->tt[x][2], collect(1:length(w)))[1:5:end]
    z_start = map(x->tt[x][3], collect(1:length(w)))[1:5:end]

    tt = [w[x] + vec(base_vec) + 50 * 5 for x=1:length(w)]
    x_end = map(x->tt[x][1], collect(1:length(w)))[1:5:end]
    y_end = map(x->tt[x][2], collect(1:length(w)))[1:5:end]
    z_end = map(x->tt[x][3], collect(1:length(w)))[1:5:end]

    for i=1:length(x_start)
        x = [x_start[i], x_end[i]]
        y = [y_start[i], y_end[i]]
        z = [z_start[i], z_end[i]]
        plot3D(x, y, z, color="blue")
    end

    info("Calculation finished")
    #PyPlot.plot(ee, ss)
    #=
    plot3D(eig_vals[:, 1], eig_vals[:, 2], eig_vals[:, 3], color="red")
    PyPlot.title("Stress path and von Mises yield surface")
    PyPlot.xlabel("Eig Stress 1")
    PyPlot.ylabel("Eig Stress 2")
    PyPlot.zlabel("Eig Stress 3")
    PyPlot.grid()
    PyPlot.show()
    =#
end

function test_von_mises_planestress_basic()

    steps = 1000
    strain_max = 0.003
    num_cycles = 5
    E = 200.0e3
    nu =  0.3
    ν = 0.3
    C = stiffnessTensorPlaneStress(E, ν)

    strain_tot = zeros(Float64, (steps, 3))

    # Adding only strain in x-axis and counting for the poisson effect
    strain_tot[:, 1] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps))
    strain_tot[:, 2] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)).*-ν
    strain_tot[:, 3] = strain_max * sin(2 * pi * linspace(0, num_cycles, steps)).*-ν

    strain_last = zeros(Float64, (3))
    strain_p = zeros(Float64, (3))
    stress = zeros(Float64, (3, 1))
    stress_y =  200.0
    ss = Float64[]
    ee = Float64[]

    ss2 = Float64[]
    ee2 = Float64[]

    eig_stress = zeros(Float64, (3, 3))
    eig_vals = zeros(Float64, (steps, 3))
    #mat = State(C, stress_y, zeros(Float64, 6), zeros(Float64, 6))

    info("Starting calculation")
    tic()
    #=
    for i=1:steps
        strain_new = reshape(strain_tot[i, :, :], (6, 1))
        dstrain = strain_new - mat.strain
        calculate_stress!(dstrain, mat, Val{:vonMises})
        mat.strain += vec(dstrain)
        push!(ss, mat.stress[1])
        push!(ee, mat.strain[1])

        fill_tensor(eig_stress, mat.stress)
        eig_vals[i, :] = sort(eigvals(eig_stress))
    end
    =#
    stress = zeros(Float64, 3)
    strain = zeros(Float64, 3)
    for i=1:steps
        strain_new = reshape(strain_tot[i, :, :], (3, 1))
        dstrain = strain_new - strain
        stress_inc, lambda = calculate_stress(dstrain,
                                              stress,
                                              C,
                                              stress_y,
                                              Val{:vonMises},
                                              Val{:PlaneStressElasticPlasticProblem})
        stress += stress_inc
        strain = vec(strain_new)
        s1, s2, t12  = stress
        se1 = (s1 + s2)/2 + sqrt(((s1 - s2)/2)^2 + t12^2)
        se2 = (s1 + s2)/2 - sqrt(((s1 - s2)/2)^2 + t12^2)
        push!(ss, se1)
        push!(ee, se2)
    end

    toc()

    function vm_upper(a, c)
        vals = f(a[1], a[2], c)
        vm(vals[1], vals[2], 200)
    end
    vm(a,b) = sqrt(a^2 - a*b + b^2) - 200
    f(m,c) = [600*cos(c) 600*sin(c)].*m
    x_vals = []
    max_iter = 100
    y_vals = []
    for i=0:0.1:(2*pi+0.3)
        wf(x) = f(x, i)
        t = 0.01
        step = 2
        merkki = -1
        s11, s22 = wf(t)
        ii = 0
        while (abs(vm(s11, s22)) > 1e-7) && ii < max_iter
            val = vm(s11, s22)
            if sign(val) != merkki
                merkki *= -1
                step *= -0.5
            end
            t += step
            s11, s22 = wf(t)
            ii += 1
        end
        push!(x_vals, s11)
        push!(y_vals, s22)
    end
    #=
    PyPlot.plot(x_vals, y_vals)
    PyPlot.plot(ee, ss)
    PyPlot.grid()
    PyPlot.show()
    =#
end

# test_von_mises_3D_basic()

#test_von_mises_planestress_basic()


text:
using JuMP
using MathOptInterface
using GLPK
using Test
using MathOptSetDistances
const MOI = MathOptInterface
const MOD = MathOptSetDistances

include("types.jl")
include("test.jl")


text:
# Julia wrapper for header: /usr/local/include/libavutil/downmix_info.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_downmix_info_update_side_data

function av_downmix_info_update_side_data(frame)
    ccall((:av_downmix_info_update_side_data,libavutil),Ptr{AVDownmixInfo},(Ptr{AVFrame},),frame)
end


text:
module RANSAC

using LinearAlgebra
using Random
using Logging
using Logging: default_logcolor
using StaticArrays: SVector, SMatrix, MMatrix, StaticArray
using RegionTrees
#using Images: label_components, component_lengths, component_subscripts
using ExtractMacro
import JSON
import YAML

import RegionTrees: AbstractRefinery, needs_refinement, refine_data
import Base: show, length, deleteat!

# debuuugggggg
#using Infiltrator
#using Makie: scatter, vbox

export  rodriguesdeg,
        rodriguesrad

export  arbitrary_orthogonal,
        arbitrary_orthogonal2,
        isparallel,
        findAABB,
        unitdisk2square,
        smallestdistance,
        prob,
        havesameelement

export  ConfidenceInterval,
        notsoconfident,
        isoverlap,
        E

export  OctreeNode,
        OctreeRefinery,
        nomodRANSACCloud,
        RANSACCloud

export  FittedShape,
        FittedPlane,
        FittedSphere,
        FittedCylinder,
        FittedCone,
        ExtractedShape,
        ransacparameters,
        setfloattype

export  project2plane,
        refit

export  ransac

export  nosource_debuglogger,
        nosource_infologger,
        sourced_debuglogger,
        sourced_infologger,
        nosource_IterInflog,
        nosource_IterLow1log,
        nosource_IterLow2log,
        nosource_Compute1log,
        nosource_Compute2log,
        nosource_Errorlog

export  exportJSON,
        readconfig

include("utilities.jl")
include("confidenceintervals.jl")
include("octree.jl")
include("fitting.jl")
include("shapes/plane.jl")
include("shapes/sphere.jl")
include("shapes/cylinder.jl")
include("shapes/cone.jl")
#include("parameterspacebitmap.jl")
include("iterations.jl")
include("logging.jl")
include("orientedbox_.jl")
include("json-yaml.jl")

"""
`const DEFAULT_PARAMETERS =
    defaultparameters([FittedPlane, FittedCone, FittedCylinder, FittedSphere])`
"""
const DEFAULT_PARAMETERS = defaultparameters([FittedPlane, FittedCone, FittedCylinder, FittedSphere])

"""
`const DEFAULT_SHAPE_DICT = Dict("plane"=>FittedPlane,
    "cone"=>FittedCone, "cylinder"=>FittedCylinder, "sphere"=>FittedSphere)`
"""
const DEFAULT_SHAPE_DICT = Dict("plane"=>FittedPlane, "cone"=>FittedCone, "cylinder"=>FittedCylinder, "sphere"=>FittedSphere)

end #module


text:
using FFTW
using DSP
using PyPlot

include("./plotfft.jl")

# System variables
Fs = 48000
n = 1:48000

#Input signal

Fsig = sin.(2*pi*440/Fs .* n) + sin.(2*pi*16000/Fs .* n)
#Fsig = zeros(48000)

#DSP functions

#Do upsampling by inserting zero signals
function upsample(signal,factor)
    out = Float64[]
    for s = signal
        push!(out,s)
        for i = 1:factor-1
            push!(out,0.0)
        end
    end
    return out
end

#Do moving average filter y[n] = (x[n] + x[n-1])/2
function movingavg(sig)
    out = zeros(size(sig)[1])
    for i = 2:size(sig)[1]
        out[i] = (sig[i]+sig[i-1])
    end
return out
end

#Implement the final quantizer
function comparator(signal)
    if signal > 0.0001 return 1.0
    elseif signal < -0.0001 return  -1.0
    else return 0
    end
end

#2rd order delta sigma DAC

function dsdac(signal)
    outarray = zeros(size(signal)[1])
    reg1 = 0.0
    reg2 = 0.0
    reg3 = 0.0
    out = 0.0
    out2 = 0.0
    ret = 1.0
    for i = 1:size(signal)[1]
        input = signal[i]
        out = reg3
        outarray[i] = out
        ret = comparator(out)
        reg3 = reg3 + input - ret
    end
    return outarray
end

#Another delta sigma implementation that closely follows how it would be implemented on hardware
function dsdac2(signal)
    outarray = zeros(size(signal)[1])
    reg_in = [0.0,0.0,0.0,0.0]
    reg_out = [0.0,0.0,0.0,0.0]

    for i = 1:size(signal)[1]
        bitstream = comparator(reg_out[2])
        outarray[i] = reg_out[2]
        input = signal[i]
        integrator_1_out = input - bitstream + reg_out[1]
        reg_in[2] = integrator_1_out - bitstream + reg_out[2]
        #reg_in[3] = reg_out[2] - bitstream + reg_out[3]
        #reg_in[4] = reg_out[2] - bitstream + reg_out[4]
        reg_out = deepcopy(reg_in)
    end
    return outarray
end

#Test delta sigma DAC by using DSP.jl built in interpolator and upsampled signal filtered with 3x moving average filter
good_resample = comparator.(dsdac2(resample(Fsig,128)))
diy_resample = comparator.(dsdac2(movingavg(movingavg(movingavg(upsample(resample(Fsig,8),32))))))


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "oneTBB"
version = v"2021.2.1"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/oneapi-src/oneTBB.git",
              # Note: this isn't actually 2021.2.1 (there is no such versions at
              # the moment), this seems to be the revision where most platforms
              # are successful
              "9e15720bc7744f85dff611d34d65e9099e077da4"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/oneTBB/

# Adapt patch from
# https://github.com/oneapi-src/oneTBB/pull/203
atomic_patch -p1 ../patches/musl.patch

mkdir build && cd build/
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DTBB_TEST=OFF \
    -DTBB_EXAMPLES=OFF \
    ..
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms(; experimental=true))
# Windows with MinGW at the moment doesn't work, but it may change in the near
# future, watch out <https://github.com/oneapi-src/oneTBB/pull/351>. See also
# <https://stackoverflow.com/q/67572880/2442087>.
filter!(!Sys.iswindows, platforms)
# Disable platforms unlikely to work
filter!(p -> arch(p) ∉ ("armv6l", "armv7l"), platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libtbbmalloc", :libtbbmalloc),
    LibraryProduct("libtbbmalloc_proxy", :libtbbmalloc_proxy),
    LibraryProduct("libtbb", :libtbb),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: plut/Yggdrasil

max_issues_repo_name: plut/Yggdrasil

max_forks_repo_name: plut/Yggdrasil


text:
macro until(condition, expression)
    quote
        # Enquanto a condição for FALSA
        #  a expressão é executada
        while !($condition)
            $expression
        end
    end
end

max_stars_repo_name: prcastro/IntroJulia

max_issues_repo_name: prcastro/IntroJulia

max_forks_repo_name: prcastro/IntroJulia


text:
using Unitful: NoUnits
using UnitfulNuclear
using Base.Test

@testset "UnitfulNuclear.jl" begin
  @testset "unit tests" begin
    cur_dict = Dict(
      1u"u" => ( 1.66054e-27*Unitful.kg ),
      1u"n20" => ( 1e20*u"m^-3" ),
      1u"b" => ( 1e-24*u"cm^2" )
    )

    expected_value = 1.0

    for (cur_key, cur_value) in cur_dict
      actual_value = ( cur_key / cur_value )
      actual_value = actual_value |> NoUnits

      @test isapprox(expected_value, actual_value, rtol=5e-5)
    end

    @test eps_0 == Unitful.ϵ0
    @test mu_0 == Unitful.μ0
  end
end


text:
struct SubVis
    p::SubHuntPOMDP
    a::Nullable{Any}
    r::Nullable{Any}
    s::Nullable{Any}
    o::Nullable{Any}
    b::Nullable{Any}
end

SubVis(p; s=nothing, a=nothing, o=nothing, b=nothing, r=nothing) = SubVis(p, a, r, s, o, b)
SubVis(p::SubHuntPOMDP, arspobp::Tuple) = SubVis(p, arspobp...)

function Base.show(io::IO, mime::MIME"text/plain", v::SubVis)
    for y in v.p.size:-1:1
        for x in 1:v.p.size
            printed = false
            if !isnull(v.s)
                s = get(v.s)
                if Pos(x,y) == s.target
                    if s.aware
                        print(io, 'A')
                    else
                        print(io, 'U')
                    end
                    printed = true
                elseif Pos(x,y) == s.own
                    if isnull(v.a)
                        print(io, 'O')
                    else
                        print(io, ACT_LETTER[get(v.a)])
                    end
                    printed = true
                end
            end
            if !printed
                print(io, '.')
            end
            print(io, ' ')
        end
        print(io, '\n')
    end
    if !isnull(v.r)
        @printf("Reward: %8.2f\n", get(v.r))
    end
end

@recipe function f(v::SubVis)
    size = v.p.size
    xlim --> (1, size)
    ylim --> (1, size)
    if !isnull(v.b)
        b = get(v.b)
        @assert b isa ParticleCollection
        counts = zeros(size, size)
        allend = true
        for s in particles(b)
            if s != END_KILL
                allend = false
                counts[s.target[1], s.target[2]] += 1
            end
        end
        if !allend
            @series begin
                x = linspace(1, size, 200)
                y = linspace(1, size, 200)
                seriestype := :contour
                fill := :true
                seriescolor --> :Oranges
                x, y, (x,y) -> counts[round(Int,x), round(Int,y)]
            end
        end
    end
    if !isnull(v.s)
        s = get(v.s)
        @series begin
            label := "Sub"
            marker := :circle
            markercolor := :blue
            [s.own[1]], [s.own[2]]
        end
        @series begin
            label := "Target"
            marker := :x
            markerstrokecolor := :blue
            markerstrokewidth := 6
            [s.target[1]], [s.target[2]]
        end
        #=
        dir = DIR[s.goal]
        @series begin
            label := "Target"
            seriestype := :quiver
            quiver := ([dir[1]], [dir[2]])
            [s.target[1]], [s.target[2]]
        end
        =#
    end
end


text:
#   This file is part of Maxima.jl. It is licensed under the MIT license
#   Copyright (c) 2019 Nathan Smith
import REPL.LineEdit
import ReplMaker

ans = nothing

"""
    finished(s)

Examine the buffer in the repl to see if the input is complete. In the Maxima
repl an expression is terminated with ';' or '\$'. '\$' is used to supress
printing the output.
"""
function finished(s)
    str = String(take!(copy(LineEdit.buffer(s))))
    if length(str) == 0
        return false
    elseif str[end] == ';' || str[end] == '$'
        return true
    else
        return false
    end
end

"""
    parse(input)

Called on valid input from the repl. Evaluates the maxima expression and
display in terminal.
"""
function parser(input)
    # '%' is a reference to the previous output in maxima. The real underlying
    # maxima server we're talking to has this variable corrupted by many
    # Maxima.jl commands between repl inputs so we emulate the behaviour
    # manually by substituting in the last output from a global variable 'ans'
    if '%' in input
        tail = last(input)                              # Chop the tail
        input = subst("$(ans)", "%", input[1:end-1])    # Substitute
        input = string(input, tail)                     # add the tail back
    end
    if !isempty(strip(input))
        global ans = MExpr(input[1:end-1]) |> mcall
        if input[end] == ';'
            display(ans)
        end
    end
end

"""
    replinit(repl)

Start up the Maxima repl mode.
"""
function repl_init(repl)
    ReplMaker.initrepl(
        parser,
        valid_input_checker=finished,
        prompt_text="maxima> ",
        prompt_color=:cyan,
        start_key=')',
        mode_name="maxima",
        repl=repl,
    )
end


text:
#
# METHOD TEMPLATES
# basic method template: "This method has priority over the more general T:<Union{....} template"
#set_wfm(fgen::INSTR{:Agilent33500B},func="sinusoid",ch=1) = write(fgen.obj,"SOURce$ch:FUNC $func")
# more general method template, superseded by basic method template
#set_wfm(fgen::T;func="sinusoid",ch=1) where {T<:Union{F332x1,F335x2}} = write(fgen.obj,"SOURce$ch:FUNC $func")

##
# Reset
reset_instr(obj::T) where {T<:Union{E364x1,E364x2}} = write(obj,"*RST")
# Clear
clear_instr(obj::T) where {T<:Union{E364x1,E364x2}} = write(obj,"*CLS")
# IDN
get_idn(obj::T) where {T<:Union{E364x1,E364x2}} = query(obj,"*IDN?")

## OUTPUT FUNCTION
# Set voltage range
set_range(obj::T;ch::String="1",vrang="low") where {T<:E364x1} = write(obj,"VOLT:RANG $vrang")
set_range(obj::T;ch::String="1",vrang="low") where {T<:E364x2} = write(obj,"INST output$ch;VOLT:RANG $vrang")
# Set output voltage
set_volt(obj::T;ch::String="1",volt=0) where {T<:E364x1} = write(obj,"VOLT $volt")
set_volt(obj::T;ch::String="1",volt=0) where {T<:E364x2} = write(obj,"INST output$ch;VOLT $volt")
# Set compliance (current)
set_compl(obj::T;ch::String="1",crtlim=0) where {T<:E364x1} = write(obj,"CURR $crtlim")
set_compl(obj::T;ch::String="1",crtlim=0) where {T<:E364x2} = write(obj,"INST output$ch;CURR $crtlim")

## OUTPUT
# Set output on/off
set_outp(obj::T;ch::String="1",st="off") where {T<:E364x1} = write(obj,"OUTP $st")
set_outp(obj::T;ch::String="1",st="off") where {T<:E364x2} = write(obj,"INST output$ch;OUTP $st")

# high level example with multiple parameters
#function set_volt_compl(obj::T;ch::String="1",volt=0,crtlim=0.1) where {T<:Union{E364x1,E364x2}}
#	set_volt(psu,ch,volt)
#	sleep(0.1)
#	set_compl(psu,ch,crtlim)
#end

# generic, arbitray number of inputs. WIP
#function set_generic(obj::T;args...) where {T<:Union{E364x1,E364x2}}
#	args[1](psu,args[2])
#end


text:
function count_heads(n)
    c::Int = 0
    for i=1:n
        c += rand(Bool)
    end
    c
end

max_stars_repo_name: snowdj/CoursePack

max_issues_repo_name: snowdj/CoursePack

max_forks_repo_name: snowdj/CoursePack


text:
module test_octo_as

using Test
using Octo.Adapters.SQL # from to_sql

@test to_sql([as(AVG(:n), :avg_n)]) == "AVG(n) AS avg_n"

struct User
end

u = from(User)
n = as(u.name, :n)

@test to_sql([n]) == "name AS n"

end # module test_octo_as


text:
"""
`prior(parameters::ParameterVector{T}) where {T<:Number}`

Calculates log joint prior density of m.parameters.
"""
function prior(parameters::ParameterVector{T}) where {T<:Number}
    free_params = Base.filter(θ -> !θ.fixed, parameters)
    logpdfs = map(logpdf, free_params)
    return sum(logpdfs)
end

"""
```
posterior(loglikelihood::Function, parameters::ParameterVector,
                   data::AbstractArray; ϕ_smc::Float64 = 1.)
```

Calculates and returns the log of the posterior distribution for `m.parameters`:

```
log posterior  = log likelihood + log prior + const
log Pr(Θ|data) = log Pr(data|Θ) + log Pr(Θ) + const
```

### Arguments

- `loglikelihood`: Function which takes vector of parameters and data, and
        outputs the log-loglikelihood of a parameters given data.
- `parameters`: Vector of parameters whose likelihood is to be evaluated
- `data`: Matrix of data for observables

### Optional Arguments

- `φ_smc`: a tempering factor to change the relative weighting of the prior and
     the likelihood when calculating the posterior. It is used primarily in SMC.
"""
function posterior(loglikelihood::Function, parameters::ParameterVector,
                   data::AbstractArray; ϕ_smc::Float64 = 1.)
    like = loglikelihood(parameters, data)
    post = ϕ_smc * like + prior(parameters)
    return post
end

"""
```
posterior!(loglikelihood::Function, parameters::ParameterVector,
           para_draw::Vector{T}, data::Matrix{T};
           sampler::Bool = false, catch_errors::Bool = false,
           φ_smc = 1) where {T<:AbstractFloat}
```

Evaluates the log posterior density at `parameters`.

### Arguments

- `loglikelihood`: a function which takes vector of parameters and data, and
                   outputs the log-loglikelihood of a parameters given data.

- `para_draw`: New values for the model parameters
- `data`: Matrix of input data for observables

### Optional Arguments
- `sampler`: Whether metropolis_hastings or smc is the caller. If `sampler=true`,
     the log likelihood and the transition matrices for the zero-lower-bound
     period are also returned.
- `catch_errors`: Whether to catch errors of type `GensysError` or `ParamBoundsError`
     If `sampler = true`, both should always be caught.
- `φ_smc`: a tempering factor to change the relative weighting of the prior and
     the likelihood when calculating the posterior. It is used primarily in SMC.
"""
function posterior!(loglikelihood::Function, parameters::ParameterVector,
                    para_draw::Vector{T}, data::Matrix{T};
                    sampler::Bool = false, catch_errors::Bool = false,
                    ϕ_smc::Float64 = 1.) where {T<:AbstractFloat}
    catch_errors = catch_errors | sampler
    if sampler
        try
            update!(parameters, para_draw)
        catch err
            if isa(err, ParamBoundsError)
                return -Inf
            else
                throw(err)
            end
        end
    else
        update!(parameters, para_draw)
    end
    return posterior(loglikelihood, parameters, data; ϕ_smc = ϕ_smc)
end


text:
# replace all vanilla spikes with jumpthru fixer spikes (dash through spikes with "Zero Speed Only" checked)

local debugSpikeMap = Dict{String, String}(
    "spikesUp" => "NerdHelper/DashThroughSpikesUp",
    "spikesDown" => "NerdHelper/DashThroughSpikesDown",
    "spikesLeft" => "NerdHelper/DashThroughSpikesLeft",
    "spikesRight" => "NerdHelper/DashThroughSpikesRight"
)

local debugSizeKey = Dict{String, String}(
    "spikesUp" => "width",
    "spikesDown" => "width",
    "spikesLeft" => "height",
    "spikesRight" => "height"
)

for room in loadedState.side.map.rooms
	for entity in room.entities
		if entity.name in ("spikesUp", "spikesDown", "spikesLeft", "spikesRight")
            newSpike = Entity(debugSpikeMap[entity.name], Dict{String,Any}(
                "x" => entity.data["x"],
				"y" => entity.data["y"],
                debugSizeKey[entity.name] => entity.data[debugSizeKey[entity.name]],
                "type" => entity.data["type"],
                "zero_speed_only" => true,
                "red_boosters_count_as_dash" => true,
                "along" => true,
                "into" => true,
                "diag" => true,
                "invert" => false
            ), entity.id)
            push!(room.entities, newSpike)
        end
    end
    filter!(entity -> entity.name ∉ ("spikesUp", "spikesDown", "spikesLeft", "spikesRight"), room.entities)
end

max_stars_repo_name: Kalobi/NerdHelper

max_issues_repo_name: Kalobi/NerdHelper

max_forks_repo_name: Kalobi/NerdHelper


text:
abstract type AbstractIndirectTrajOptimizer end

struct IndirectTrajOptimizer{IPT,CSIT,ST} <: AbstractIndirectOptimizationProblem
    # Indirect Trajectory Optimization problem
    prob::IPT

    # Indirect Co-State Initializer
    csInit::CSIT

    # Indirect solver
    solver::ST

    # Data Output Manager
    writingData::Bool
    dataOutputManager::DataOutputManager

    # Meta data (May want to add more here)
    solMethod::Symbol
    initCost::Symbol
    initOptimizer::Symbol
    weights::Vector{Float64}
    minNeighborhoodFrac::Float64
    UBs::Vector{Float64}
    LBs::Vector{Float64}
    iUBs::Vector{Float64}
    iLBs::Vector{Float64}
    maxIters::Int
    funcTol::Float64
    maxStallIters::Float64
    maxStallTime::Float64
    maxTime::Float64
    useParallel::Bool
end

function IndirectTrajOptimizer(prob; solutionMethod = :FSS, initCostFunc = :WSS,
    initOptimizer = :PSO, numParticles = 500, numSwarms = 4, swarmInitMethod = :Uniform,
    UBs = [100, 100, 100, 50, 50, 50, 50], LBs = [-100, -100, -100, -50, -50, -50, -50],
    iUBs = nothing, iLBs = nothing, weights = [10, 10, 10, 1, 1, 1, 1], MFD = 2.8e-6, display = true,
    displayInterval = 1, maxIters = 1000, funcTol = 1e-6, maxStallIters = 25, maxStallTime = 500,
    maxTime = 1800, useParallel = true, homotopyParamVec = nothing, dataFolder = nothing,
    minNeighborhoodFraction = 0.25,writeData = false)

    # Check that homotopy parameter vector has been set if homotopy is used
    if prob.homotopy && homotopyParamVec === nothing
        throw(ArgumentError("Problem specified to use homotopy continuation but continuation parameters not provided."))
    end

    # Initialize initializer and solver
    if solutionMethod == :FSS
        # Initialize co-state initializer
        csInit = FSSCoStateInitializer(prob.initBVPFunc, prob.iConds, prob.fConds;
            costFunc = initCostFunc, optimizer = initOptimizer, numParticles = numParticles,
            numSwarms = numSwarms, initMethod = swarmInitMethod, UBs = UBs, LBs = LBs,
            iUBs = iUBs, iLBs = iLBs, weights = weights, MFD = MFD, display = display,
            displayInterval = displayInterval, maxIters = maxIters, funcTol,
            maxStallIters = maxStallIters, maxStallTime = maxStallTime, maxTime = maxTime,
            minNeighborhoodFraction = minNeighborhoodFraction, useParallel = useParallel)

        # Initialize solver
        solver = FSSSolver(zeros(7), prob.iConds, prob.fConds, prob.BVPFunc, prob.BVPWithSTMFunc;
            homotopy = prob.homotopy, homotopyParamVec = homotopyParamVec)
    else
        throw(ArgumentError("Only forward sigle shooting is implemented."))
    end

    # Initialize DataOutputManager
    if dataFolder === nothing
        dataFolder = joinpath(pwd(),"data")
    end
    dataOutputManager = DataOutputManager(dataFolder)

    if iUBs === nothing
        iUBs = Vector{Float64}(undef, 0)
    end
    if iLBs === nothing
        iLBs = Vector{Float64}(undef, 0)
    end

    IndirectTrajOptimizer{typeof(prob), typeof(csInit), typeof(solver)}(prob, csInit, solver,
        writeData, dataOutputManager, solutionMethod, initCostFunc, initOptimizer, weights,
        minNeighborhoodFraction, UBs, LBs, iUBs, iLBs, maxIters, funcTol, maxStallIters,
        maxStallTime, maxTime, useParallel)
end

function initialize!(ito::IndirectTrajOptimizer)
    # Initialize co-states
    initialize!(ito.csInit)

    # Re-initialize solvers data manager
    initializeData!(ito.solver, GetInitializedCostates(ito.csInit))

    # Write data if desired
    if ito.writingData
        writeData(ito.dataOutputManager, ito)
    end

    return nothing
end

function solve!(ito::IndirectTrajOptimizer; factor = 3.0, ftol = 1e-8, showTrace = true, convergenceAttempts = 4)

    # Solve
    solve!(ito.solver; factor = factor, ftol = ftol, showTrace = showTrace, convergenceAttempts = convergenceAttempts)

    # Write data if desired
    if ito.writingData
        writeData(ito.dataOutputManager, ito)
    end

    return nothing
end

function tSolve!(itoVec; factor = 3.0, ftol = 1e-8, showTrace = true, convergenceAttempts = 4)
    p = Progress(length(itoVec), 1, "Solving BVPs: ")
    Threads.@threads for i in 1:length(itoVec)
        solve!(itoVec[i]; factor = factor, ftol = ftol, showTrace = showTrace, convergenceAttempts = convergenceAttempts)
        next!(p)
    end
    return nothing
end


text:
module Geometry

export Point, Point2D, Point3D, Rectangle, vertices

abstract type Point{N,T} end
"""
    Point2D(x, y)

Generate a two-dimensional point.
"""
struct Point2D{T} <: Point{2,T}
    x::T
    y::T
end
function Point2D(x, y)
    T = Base.promote_typeof(x, y)
    return Point2D{T}(x, y)
end
struct Point3D{T} <: Point{3,T}
    x::T
    y::T
    z::T
end
"""
    Point3D(x, y, z)

Generate a three-dimensional point.
"""
function Point3D(x, y, z)
    T = Base.promote_typeof(x, y, z)
    return Point3D{T}(x, y, z)
end

"""
    Rectangle(lx, rx, ly, uy)

Generate a rectangle from its four corner points.
"""
struct Rectangle{T}
    lx::T
    rx::T
    ly::T
    uy::T
    function Rectangle(lx, rx, ly, uy)
        @assert lx != rx "The `lx` and `rx` arguments cannot be the same!"
        @assert ly != uy "The `ly` and `uy` arguments cannot be the same!"
        T = Base.promote_typeof(lx, rx, ly, uy)
        return new{T}(min(lx, rx), max(lx, rx), min(ly, uy), max(ly, uy))
    end
end
Rectangle(rec::Rectangle) = rec

"""
    vertices(rec::Rectangle)

Return four corner points of a `Rectangle`.
"""
function vertices(rec::Rectangle)
    lx, rx, ly, uy = rec.lx, rec.rx, rec.ly, rec.uy
    return Point2D(lx, ly), Point2D(lx, uy), Point2D(rx, ly), Point2D(rx, uy)
end

"""
    Base.in(p::Point2D, rec::Rectangle)

Test if point `p` is in rectangle `rec`.
"""
function Base.in(p::Point2D, rec::Rectangle)
    lx, rx, ly, uy = rec.lx, rec.rx, rec.ly, rec.uy
    return lx <= p.x <= rx && ly <= p.y <= uy
end

end


text:
export MOBIL

"""
    MOBIL
See Treiber & Kesting, 'Modeling Lane-Changing Decisions with MOBIL'
"""
type MOBIL <: LaneChangeModel{LaneChangeChoice}

    dir::Int
    rec::SceneRecord
    mlon::LongitudinalDriverModel
    safe_decel::Float64 # safe deceleration (positive value)
    politeness::Float64 # politeness factor (suggested p ∈ [0.2,0.5])
    advantage_threshold::Float64 # Δaₜₕ

    function MOBIL(
        action_context::IntegratedContinuous;
        rec::SceneRecord=SceneRecord(2,action_context.Δt),
        mlon::LongitudinalDriverModel=IntelligentDriverModel(),
        safe_decel::Float64=2.0, # [m/s²]
        politeness::Float64=0.35,
        advantage_threshold::Float64=0.1,
        )

        retval = new()
        retval.dir = DIR_MIDDLE
        retval.rec = rec
        retval.mlon = mlon
        retval.safe_decel = safe_decel
        retval.politeness = politeness
        retval.advantage_threshold = advantage_threshold
        retval
    end
end
get_name(::MOBIL) = "MOBIL"
function set_desired_speed!(model::MOBIL, v_des::Float64)
    set_desired_speed!(model.mlon, v_des)
    model
end
function observe!(model::MOBIL, scene::Scene, roadway::Roadway, egoid::Int)

    rec = model.rec
    update!(rec, scene)

    vehicle_index = get_index_of_first_vehicle_with_id(rec, egoid)
    veh_ego = scene[vehicle_index]
    v = veh_ego.state.v
    egostate_M = veh_ego.state

    left_lane_exists = convert(Float64, get(N_LANE_LEFT, rec, roadway, vehicle_index)) > 0
    right_lane_exists = convert(Float64, get(N_LANE_RIGHT, rec, roadway, vehicle_index)) > 0
    fore_M = get_neighbor_fore_along_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointRear(), VehicleTargetPointFront())
    rear_M = get_neighbor_rear_along_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointFront(), VehicleTargetPointRear())

    # accel if we do not make a lane change
    accel_M_orig = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, egoid))
    model.dir = DIR_MIDDLE

    advantage_threshold = model.advantage_threshold
    if right_lane_exists

        rear_R = get_neighbor_rear_along_right_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointFront(), VehicleTargetPointRear())

        # candidate position after lane change is over
        footpoint = get_footpoint(veh_ego)
        lane = roadway[veh_ego.state.posF.roadind.tag]
        lane_R = roadway[LaneTag(lane.tag.segment, lane.tag.lane - 1)]
        roadproj = proj(footpoint, lane_R, roadway)
        frenet_R = Frenet(RoadIndex(roadproj), roadway)
        egostate_R = VehicleState(frenet_R, roadway, veh_ego.state.v)

        Δaccel_n = 0.0
        passes_safety_criterion = true
        if rear_R.ind != 0
            id = scene[rear_R.ind].def.id
            accel_n_orig = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, id))
            veh_ego.state = egostate_R
            accel_n_test = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, id))

            body = inertial2body(get_rear_center(scene[rear_R.ind]), get_front_center(veh_ego)) # project target to be relative to ego
            s_gap = body.x

            veh_ego.state = egostate_M
            passes_safety_criterion = accel_n_test ≥ -model.safe_decel && s_gap ≥ 0
            Δaccel_n = accel_n_test - accel_n_orig
        end

        if passes_safety_criterion

            Δaccel_o = 0.0
            if rear_M.ind != 0
                id = scene[rear_M.ind].def.id
                accel_o_orig = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, id))
                veh_ego.state = egostate_R
                accel_o_test = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, id))
                veh_ego.state = egostate_M
                Δaccel_o = accel_o_test - accel_o_orig
            end

            veh_ego.state = egostate_R
            accel_M_test = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, egoid))
            veh_ego.state = egostate_M
            Δaccel_M = accel_M_test - accel_M_orig

            Δaₜₕ = Δaccel_M + model.politeness*(Δaccel_n + Δaccel_o)
            if Δaₜₕ > advantage_threshold
                model.dir = DIR_RIGHT
                advantage_threshold = Δaₜₕ
            end
        end
    end

    if left_lane_exists
        rear_L = get_neighbor_rear_along_left_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointFront(), VehicleTargetPointRear())

        # candidate position after lane change is over
        footpoint = get_footpoint(veh_ego)
        lane = roadway[veh_ego.state.posF.roadind.tag]
        lane_L = roadway[LaneTag(lane.tag.segment, lane.tag.lane + 1)]
        roadproj = proj(footpoint, lane_L, roadway)
        frenet_L = Frenet(RoadIndex(roadproj), roadway)
        egostate_L = VehicleState(frenet_L, roadway, veh_ego.state.v)

        Δaccel_n = 0.0
        passes_safety_criterion = true
        if rear_L.ind != 0
            id = scene[rear_L.ind].def.id
            accel_n_orig = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, id))
            veh_ego.state = egostate_L
            accel_n_test = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, id))

            body = inertial2body(get_rear_center(scene[rear_L.ind]), get_front_center(veh_ego)) # project target to be relative to ego
            s_gap = body.x

            veh_ego.state = egostate_M
            passes_safety_criterion = accel_n_test ≥ -model.safe_decel && s_gap ≥ 0
            Δaccel_n = accel_n_test - accel_n_orig
        end

        if passes_safety_criterion

            Δaccel_o = 0.0
            if rear_M.ind != 0
                id = scene[rear_M.ind].def.id
                accel_o_orig = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, id))
                veh_ego.state = egostate_L
                accel_o_test = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, id))
                veh_ego.state = egostate_M
                Δaccel_o = accel_o_test - accel_o_orig
            end

            veh_ego.state = egostate_L
            accel_M_test = rand(observe!(reset_hidden_state!(model.mlon), scene, roadway, egoid))
            veh_ego.state = egostate_M
            Δaccel_M = accel_M_test - accel_M_orig

            Δaₜₕ = Δaccel_M + model.politeness*(Δaccel_n + Δaccel_o)
            if Δaₜₕ > advantage_threshold
                model.dir = DIR_LEFT
                advantage_threshold = Δaₜₕ
            end
        end
    end

    model
end
Base.rand(model::MOBIL) = LaneChangeChoice(model.dir)


text:
@testset "Centering" begin
    data = (x=collect(1:10),
            y=rand(10) .+ 3,
            z=Symbol.(repeat('a':'e', 2)))

    @testset "Automatic centering" begin
        xc = concrete_term(term(:x), data, Center())
        @test xc isa CenteredTerm
        @test xc.center == mean(data.x)
        @test modelcols(xc, data) == data.x .- mean(data.x) == data.x .- xc.center

        yc = concrete_term(term(:y), data, Center())
        @test yc isa CenteredTerm
        @test yc.center == mean(data.y)
        @test modelcols(yc, data) == data.y .- mean(data.y) == data.y .- yc.center

        @testset "alternative center function" begin
            f = first
            xc = concrete_term(term(:x), data, Center(f))
            @test xc isa CenteredTerm
            @test xc.center == f(data.x)
            @test modelcols(xc, data) == data.x .- f(data.x) == data.x .- xc.center
            # why test this? well this makes sure that our tests
            # wouldn't pass if we were using the default center function
            # in other words, this tests we're actually hitting a different codepath
            @test !isapprox(mean(data.x), f(data.x))
        end
    end

    @testset "Manual centering" begin
        xc = concrete_term(term(:x), data, Center(5))
        @test xc isa CenteredTerm
        @test xc.center == 5 != mean(data.x)
        @test modelcols(xc, data) == data.x .- 5 == data.x .- xc.center

        yc = concrete_term(term(:y), data, Center(5))
        @test yc isa CenteredTerm
        @test yc.center == 5 != mean(data.y)
        @test modelcols(yc, data) == data.y .- 5 == data.y .- yc.center
    end

    @testset "Schema hints dict" begin
        sch = schema(data, Dict(:x => Center(), :y => Center(2)))
        xc = sch[term(:x)]
        @test xc isa CenteredTerm
        @test !StatsModels.needs_schema(xc)
        @test xc.center == mean(data.x)

        yc = sch[term(:y)]
        @test yc isa CenteredTerm
        @test yc.center == 2

        @test modelcols(xc, data) == data.x .- mean(data.x)
        @test modelcols(yc, data) == data.y .- 2
    end

    @testset "center function" begin
        sch = schema(data)
        x = sch[term(:x)]

        xc = center(x)
        @test xc isa CenteredTerm
        @test xc.center == x.mean

        xc2 = center(x, Center(2))
        xc22 = center(x, 2)
        @test xc2.center == xc22.center == 2

        @testset "centering vectors" begin
            x = collect(1:5)
            y = x .+ 5.0
            xc = x .- 3

            @test center(x) == xc
            @test center(y) == xc
            @test center(maximum, x) == (x .- 5)
            @test all(center(y, x) .== 5)

            y2 = copy(y)
            @test center!(y2) == xc
            @test all(y2 .!= y)
            @test y2 == xc

            # this is int on int so you can modify in place
            @test center!(maximum, copy(x)) == (x .- 5)
            # this converts exactly to int, so you can do it in place
            center!(mean, copy(x)) == center(mean, x)
            @test_throws ArgumentError center!(mean, [1, 2])
            @test_throws ArgumentError center!([1, 2])
            @test_throws MethodError center!(v -> 1, ["a", "b"])
        end
    end

    @testset "plays nicely with formula" begin
        f = @formula(0 ~ x * y)

        sch = schema(f, data)

        sch_c = schema(f, data, Dict(:x => Center(), :y => Center(2)))
        ff_c = apply_schema(f, sch_c)

        mm_c = modelcols(ff_c.rhs, data)
        @test mm_c == hcat(data.x .- mean(data.x),
                           data.y .- 2,
                           (data.x .- mean(data.x)) .* (data.y .- 2))

        @test coefnames(ff_c.rhs) ==
              ["x(centered: 5.5)", "y(centered: 2)", "x(centered: 5.5) & y(centered: 2)"]

        # round-trip schema is empty since needs_schema is false
        sch_2 = schema(ff_c, data)
        @test isempty(sch_2.schema)

        @testset "missing omit" begin
            d = (x=[1, 2, missing, 3], y=[1, missing, 2, 3], z=[missing, 1, 2, 3])
            dd, nonmissing = StatsModels.missing_omit(d, ff_c)
            @test findall(nonmissing) == [1, 4]
            @test dd.x == d.x[nonmissing]
            @test dd.y == d.y[nonmissing]
        end
    end

    @testset "printing" begin
        xc = concrete_term(term(:x), data, Center())
        @test StatsModels.termsyms(xc) == Set([:x])
        @test "$(xc)" == string_mime(MIME("text/plain"), xc) == "x(centered: 5.5)"
        @test coefnames(xc) == "x(centered: 5.5)"
    end

    @testset "categorical term" begin
        z = concrete_term(term(:z), data)
        @test_throws ArgumentError center(z)
        @test_throws ArgumentError center(z, Center())
        zc = center(z, Center(0.5))
        @test modelcols(zc, data) == modelcols(z, data) .- 0.5
        @test StatsModels.width(zc) == StatsModels.width(z)
        @test coefnames(zc) == coefnames(z) .* "(centered: 0.5)"

        zc2 = center(z, Center([1 2 3 4]))
        @test modelcols(zc2, data) == modelcols(z, data) .- [1 2 3 4]
        @test coefnames(zc2) ==
              coefnames(z) .* "(centered: " .* string.([1, 2, 3, 4]) .* ")"
    end

    # @testset "utilities" begin

    # end

end


text:
using Test
using QuantumOpticsBase

@testset "printing" begin

@test sprint(show, GenericBasis([2, 3])) == "Basis(shape=[2,3])"
@test sprint(show, GenericBasis(2)) == "Basis(dim=2)"
@test sprint(show, SpinBasis(1//1)) == "Spin(1)"
@test sprint(show, SpinBasis(3//2)) == "Spin(3/2)"
@test sprint(show, FockBasis(1)) == "Fock(cutoff=1)"
@test sprint(show, NLevelBasis(2)) == "NLevel(N=2)"

@test sprint(show, PositionBasis(-4, 4, 10)) == "Position(xmin=-4, xmax=4, N=10)"
@test sprint(show, MomentumBasis(-4, 4, 10)) == "Momentum(pmin=-4, pmax=4, N=10)"

b_fock = FockBasis(4)
states = [fockstate(b_fock, 2), coherentstate(b_fock, 0.4)]
@test sprint(show, SubspaceBasis(states)) == "Subspace(superbasis=Fock(cutoff=4), states:2)"

b_mb = ManyBodyBasis(b_fock, fermionstates(b_fock, 2))
@test sprint(show, b_mb) == "ManyBody(onebodybasis=Fock(cutoff=4), states:10)"

state = fockstate(FockBasis(2), 2)
@test sprint(show, state) == "Ket(dim=3)\n  basis: Fock(cutoff=2)\n 0.0 + 0.0im\n 0.0 + 0.0im\n 1.0 + 0.0im"
state = dagger(state)
@test sprint(show, state) == "Bra(dim=3)\n  basis: Fock(cutoff=2)\n 0.0 - 0.0im\n 0.0 - 0.0im\n 1.0 - 0.0im"

op = DenseOperator(FockBasis(1))
@test sprint(show, op) == "Operator(dim=2x2)
  basis: Fock(cutoff=1)
 0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im"

op = DenseOperator(b_fock, b_fock ⊗ SpinBasis(1//2))
@test sprint(show, op) == "Operator(dim=5x10)
  basis left:  Fock(cutoff=4)
  basis right: [Fock(cutoff=4) ⊗ Spin(1/2)]
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im"

op = SparseOperator(b_fock, b_fock ⊗ SpinBasis(1//2))
if VERSION < v"1.6.0-beta1"
@test sprint(show, op) == "Operator(dim=5x10)
  basis left:  Fock(cutoff=4)
  basis right: [Fock(cutoff=4) ⊗ Spin(1/2)]
    []"
else
@test sprint(show, op) == "Operator(dim=5x10)\n  basis left:  Fock(cutoff=4)\n  basis right: [Fock(cutoff=4) ⊗ Spin(1/2)]\n     ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅    \n     ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅    \n     ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅    \n     ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅    \n     ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅          ⋅    "
end

op = SparseOperator(b_fock)
op.data[2,2] = 1
if VERSION < v"1.6.0-beta1"
@test replace(sprint(show, op), "\t" => "  ") == "Operator(dim=5x5)
  basis: Fock(cutoff=4)\n  [2, 2]  =  1.0+0.0im"
else
@test replace(sprint(show, op), "\t" => "  ") == "Operator(dim=5x5)\n  basis: Fock(cutoff=4)\n     ⋅          ⋅          ⋅          ⋅          ⋅    \n     ⋅      1.0+0.0im      ⋅          ⋅          ⋅    \n     ⋅          ⋅          ⋅          ⋅          ⋅    \n     ⋅          ⋅          ⋅          ⋅          ⋅    \n     ⋅          ⋅          ⋅          ⋅          ⋅    "
end

op = LazySum(SparseOperator(b_fock), DenseOperator(b_fock))
@test sprint(show, op) == "LazySum(dim=5x5)
  basis: Fock(cutoff=4)
  operators: 2"

op = LazyProduct(SparseOperator(b_fock), DenseOperator(b_fock))
@test sprint(show, op) == "LazyProduct(dim=5x5)
  basis: Fock(cutoff=4)
  operators: 2"

b_fock = FockBasis(2)
b_spin = SpinBasis(1//2)
b_mb = ManyBodyBasis(b_spin, fermionstates(b_spin, 1))
op = LazyTensor(b_fock ⊗ b_mb ⊗ b_spin, [1, 3], (SparseOperator(b_fock), DenseOperator(b_spin)))
@test sprint(show, op) == "LazyTensor(dim=12x12)
  basis: [Fock(cutoff=2) ⊗ ManyBody(onebodybasis=Spin(1/2), states:2) ⊗ Spin(1/2)]
  operators: 2
  indices: [1,3]"

bx = PositionBasis(-2, 2, 4)
bp = MomentumBasis(bx)
Tpx = transform(bp, bx)
@test sprint(show, Tpx) == "FFTOperators(dim=4x4)
  basis left:  Momentum(pmin=-3.141592653589793, pmax=3.141592653589793, N=4)
  basis right: Position(xmin=-2, xmax=2, N=4)"

# Inversed tensor product ordering
QuantumOpticsBase.set_printing(standard_order=true)

n = fockstate(b_fock, 1)
@test sprint(show, n) == "Ket(dim=3)\n  basis: Fock(cutoff=2)\n 0.0 + 0.0im\n 1.0 + 0.0im\n 0.0 + 0.0im"

spin1 = spindown(b_spin)
spin2 = spinup(b_spin)
state = n ⊗ spin1 ⊗ spin2
state_data = kron(n.data, spin1.data, spin2.data)
state_data_str = sprint(Base.print_array, state_data)
@test sprint(show, state) == "Ket(dim=12)
  basis: [Fock(cutoff=2) ⊗ Spin(1/2) ⊗ Spin(1/2)]\n"*state_data_str

state_data_str = sprint(Base.print_array, conj.(state_data))
@test sprint(show, dagger(state)) == "Bra(dim=12)
  basis: [Fock(cutoff=2) ⊗ Spin(1/2) ⊗ Spin(1/2)]\n"*state_data_str

op = dm(state)
op_data = state_data * state_data'
type_len = if VERSION < v"1.6.0-beta1"
    length("Complex{Float64}")
else
    length("ComplexF64")
end
op_data_str1 = split(sprint(show, op_data)[type_len+2:end-1], ";")
for i=1:length(op_data_str1)
    op_data_str1[i] = join(split(op_data_str1[i]), "  ")
    op_data_str1[i] = replace(op_data_str1[i], "  +  "=>"+")
end
op_data_str = join(op_data_str1, "\n ")
@test sprint(show, op) == "Operator(dim=12x12)
  basis: [Fock(cutoff=2) ⊗ Spin(1/2) ⊗ Spin(1/2)]\n "*op_data_str

op = sparse(op)
op_data = sparse(op_data)
op_data_str = sprint(show, op_data)[4:end]
# @test sprint(show, op) == "Operator(dim=12x12)
#   basis: [Fock(cutoff=2) ⊗ Spin(1/2) ⊗ Spin(1/2)]\n  "*op_data_str

paulix, pauliy = sigmax(b_spin), sigmay(b_spin)
pauli = paulix ⊗ pauliy
@test sprint(show, pauli) == "Operator(dim=4x4)\n  basis: [Spin(1/2) ⊗ Spin(1/2)]\n  [4, 1]  =  0.0+1.0im\n  [3, 2]  =  0.0-1.0im\n  [2, 3]  =  0.0+1.0im\n  [1, 4]  =  0.0-1.0im"
@test sprint((io, x) -> show(IOContext(io, :limit=>true, :displaysize=>(20,80)), x), pauli ⊗ pauli) == "Operator(dim=16x16)\n  basis: [Spin(1/2) ⊗ Spin(1/2) ⊗ Spin(1/2) ⊗ Spin(1/2)]\n  [16,  1]  =  -1.0+0.0im\n  [15,  2]  =  1.0+0.0im\n  [14,  3]  =  -1.0+0.0im\n  [13,  4]  =  1.0+0.0im\n  [12,  5]  =  1.0+0.0im\n  ⋮\n  [6 , 11]  =  -1.0+0.0im\n  [5 , 12]  =  1.0+0.0im\n  [4 , 13]  =  1.0+0.0im\n  [3 , 14]  =  -1.0-0.0im\n  [2 , 15]  =  1.0+0.0im\n  [1 , 16]  =  -1.0-0.0im"
@test sprint(show, dense(pauli)) == "Operator(dim=4x4)
  basis: [Spin(1/2) ⊗ Spin(1/2)]
 0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0-1.0im
 0.0+0.0im  0.0+0.0im  0.0+1.0im  0.0+0.0im
 0.0+0.0im  0.0-1.0im  0.0+0.0im  0.0+0.0im
 0.0+1.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im"

hadamard = DenseOperator(b_spin, 1/sqrt(2) * Complex.([1 1; 1 -1]))
@test sprint(show, sigmax(b_spin) ⊗ hadamard) == "Operator(dim=4x4)\n  basis: [Spin(1/2) ⊗ Spin(1/2)]\n  [3, 1]  =  0.707107+0.0im\n  [4, 1]  =  0.707107+0.0im\n  [3, 2]  =  0.707107+0.0im\n  [4, 2]  =  -0.707107+0.0im\n  [1, 3]  =  0.707107+0.0im\n  [2, 3]  =  0.707107+0.0im\n  [1, 4]  =  0.707107+0.0im\n  [2, 4]  =  -0.707107+0.0im"
@test sprint((io, x) -> show(IOContext(io, :limit=>true, :displaysize=>(20,80)), x), dense(sigmax(b_spin) ⊗ hadamard ⊗ hadamard ⊗ hadamard)) ==
"Operator(dim=16x16)\n  basis: [Spin(1/2) ⊗ Spin(1/2) ⊗ Spin(1/2) ⊗ Spin(1/2)]\n      0.0+0.0im        0.0+0.0im  …   0.353553+0.0im   0.353553+0.0im\n      0.0+0.0im        0.0+0.0im      0.353553+0.0im  -0.353553+0.0im\n      0.0+0.0im        0.0+0.0im     -0.353553+0.0im  -0.353553+0.0im\n      0.0+0.0im        0.0+0.0im     -0.353553+0.0im   0.353553-0.0im\n      0.0+0.0im        0.0+0.0im     -0.353553+0.0im  -0.353553+0.0im\n      0.0+0.0im        0.0+0.0im  …  -0.353553+0.0im   0.353553-0.0im\n      0.0+0.0im        0.0+0.0im      0.353553+0.0im   0.353553+0.0im\n      0.0+0.0im        0.0+0.0im      0.353553+0.0im  -0.353553+0.0im\n 0.353553+0.0im   0.353553+0.0im           0.0+0.0im        0.0+0.0im\n 0.353553+0.0im  -0.353553+0.0im           0.0+0.0im        0.0+0.0im\n 0.353553+0.0im   0.353553+0.0im  …        0.0+0.0im        0.0+0.0im\n 0.353553+0.0im  -0.353553+0.0im           0.0+0.0im        0.0+0.0im\n 0.353553+0.0im   0.353553+0.0im           0.0+0.0im        0.0+0.0im\n 0.353553+0.0im  -0.353553+0.0im           0.0+0.0im        0.0+0.0im\n 0.353553+0.0im   0.353553+0.0im           0.0+0.0im        0.0+0.0im\n 0.353553+0.0im  -0.353553+0.0im  …        0.0+0.0im        0.0+0.0im"

@test sprint((io, x) -> show(IOContext(io, :limit=>true, :displaysize=>(16,80)), x), dense(sigmax(b_spin) ⊗ hadamard ⊗ hadamard ⊗ hadamard)) ==
"Operator(dim=16x16)\n  basis: [Spin(1/2) ⊗ Spin(1/2) ⊗ Spin(1/2) ⊗ Spin(1/2)]\n      0.0+0.0im        0.0+0.0im  …   0.353553+0.0im   0.353553+0.0im\n      0.0+0.0im        0.0+0.0im      0.353553+0.0im  -0.353553+0.0im\n      0.0+0.0im        0.0+0.0im     -0.353553+0.0im  -0.353553+0.0im\n      0.0+0.0im        0.0+0.0im     -0.353553+0.0im   0.353553-0.0im\n      0.0+0.0im        0.0+0.0im     -0.353553+0.0im  -0.353553+0.0im\n      0.0+0.0im        0.0+0.0im  …  -0.353553+0.0im   0.353553-0.0im\n         ⋮                        ⋱                            ⋮     \n 0.353553+0.0im  -0.353553+0.0im           0.0+0.0im        0.0+0.0im\n 0.353553+0.0im   0.353553+0.0im           0.0+0.0im        0.0+0.0im\n 0.353553+0.0im  -0.353553+0.0im           0.0+0.0im        0.0+0.0im\n 0.353553+0.0im   0.353553+0.0im           0.0+0.0im        0.0+0.0im\n 0.353553+0.0im  -0.353553+0.0im  …        0.0+0.0im        0.0+0.0im"

nlevel1 = basisstate(NLevelBasis(4), 1)
nlevel2 = basisstate(NLevelBasis(4), 2)
nlevel3 = basisstate(NLevelBasis(3), 3)
@test sprint((io, x) -> show(IOContext(io, :limit=>true, :displaysize=>(16,80)), x), (nlevel1 ⊗ nlevel2) ⊗ dagger(nlevel3)) ==
"Operator(dim=16x3)\n  basis left:  [NLevel(N=4) ⊗ NLevel(N=4)]\n  basis right: NLevel(N=3)\n 0.0+0.0im  0.0+0.0im  0.0+0.0im\n 0.0+0.0im  0.0+0.0im  1.0+0.0im\n 0.0+0.0im  0.0+0.0im  0.0+0.0im\n 0.0+0.0im  0.0+0.0im  0.0+0.0im\n 0.0+0.0im  0.0+0.0im  0.0+0.0im\n 0.0+0.0im  0.0+0.0im  0.0+0.0im\n    ⋮                           \n 0.0+0.0im  0.0+0.0im  0.0+0.0im\n 0.0+0.0im  0.0+0.0im  0.0+0.0im\n 0.0+0.0im  0.0+0.0im  0.0+0.0im\n 0.0+0.0im  0.0+0.0im  0.0+0.0im\n 0.0+0.0im  0.0+0.0im  0.0+0.0im"

# Test switching back
QuantumOpticsBase.set_printing(standard_order=false)
state_data = kron(spin2.data, spin1.data, n.data)
state_data_str = sprint(Base.print_array, state_data)
@test sprint(show, state) == "Ket(dim=12)
  basis: [Fock(cutoff=2) ⊗ Spin(1/2) ⊗ Spin(1/2)]\n"*state_data_str

end # testset


text:
using Test

@testset "Core" begin
  include("MarkovCategories.jl")
  include("StochMapsExpr.jl")
end


text:
export EGSParticle, EGSParticleZ, Latch, BitRegions
export @BitRegions

include("particle.jl")
include("io.jl")


text:
# Test linearity of sources
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#
# Mathias Louboutin, mlouboutin3@gatech.edu
# Updated July 2020

### Model
model, model0, dm = setup_model(tti, viscoacoustic, nlayer)
q1, srcGeometry1, recGeometry, f0 = setup_geom(model)
srcGeometry2 = deepcopy(srcGeometry1)
srcGeometry2.xloc[:] .= .9*srcGeometry2.xloc[:]
srcGeometry2.zloc[:] .= .9*srcGeometry2.zloc[:]
dt = srcGeometry1.dt[1]

opt = Options(free_surface=fs, f0=f0)
ftol = 5f-5

####################### Modeling operators ##########################################

@testset "Linearity test with $(nlayer) layers and tti $(tti) and viscoacoustic $(viscoacoustic) and freesurface $(fs)" begin
    @timeit TIMEROUTPUT "Linearity" begin
        # Modeling operators
        Pr = judiProjection(recGeometry)
        Ps1 = judiProjection(srcGeometry1)
        Ps2 = judiProjection(srcGeometry2)
        F = judiModeling(model; options=opt)
        q2 = judiVector(srcGeometry2,q1.data[1])

        A1 = Pr*F*adjoint(Ps1)
        A2 = Pr*F*adjoint(Ps2)
        J = judiJacobian(A1, q1)

        d1 = A1*q1
        d2 = A2*q2
        d3 = Pr*F*(adjoint(Ps1)*q1 + adjoint(Ps2)*q2)
        d4 = Pr*F*(adjoint(Ps1)*q1 - adjoint(Ps2)*q2)
        d5 = A1 * (2f0 * q1)
        d6 = (2f0 * A1) * q1

        q3 = adjoint(A1) * d1
        q4 = adjoint(A1) * (2f0 * d1)
        q5 = (2f0 * adjoint(A1)) * d1

        # Test linearity F * (a + b) == F * a + F * b

        println("Test linearity of F: F * (a + b) == F * a + F * b")
        nd1 = norm(d3)
        nd2 = norm(d1 + d2)
        nd3 = norm(d3 - d1 - d2)/norm(d3)
        @printf(" F * (a + b): %2.5e, F * a + F * b : %2.5e, relative error : %2.5e \n", nd1, nd2, nd3)
        nd1 = norm(d4)
        nd2 = norm(d1 - d2)
        nd3 = norm(d4 - d1 + d2)/norm(d4)
        @printf(" F * (a - b): %2.5e, F * a - F * b : %2.5e, relative error : %2.5e \n", nd1, nd2, nd3)

        @test isapprox(d3, d1 + d2, rtol=ftol)
        @test isapprox(d4, d1 - d2, rtol=ftol)

        # Test linearity F a x == a F x
        println("Test linearity of F: F * (a * b) == a * F * b")
        nd1 = norm(2f0 * d1)
        nd1_b = norm(d6)
        nd2 = norm(d5)
        nd3 = norm(2f0 * d1 - d5)/norm(d5)
        nd4 = norm(d6 - d5)/norm(d5)

        nm1 = norm(2f0 * q3)
        nm1_b = norm(q5)
        nm2 = norm(q4)
        nm3 = norm(2f0 * q3 - q4)/norm(q4)
        nm4 = norm(q5 - q4)/norm(q4)
        @printf(" a (F x): %2.5e, F a x : %2.5e, relative error : %2.5e \n", nd1, nd2, nd3)
        @printf(" (a F) x: %2.5e, F a x : %2.5e, relative error : %2.5e \n", nd1_b, nd2, nd4)
        @printf(" a (F' x): %2.5e, F' a x : %2.5e, relative error : %2.5e \n", nm1, nm2, nm3)
        @printf(" (a F') x: %2.5e, F' a x : %2.5e, relative error : %2.5e \n", nm1_b, nm2, nm4)
        @test isapprox(2f0 * d1, d5, rtol=ftol)
        @test isapprox(2f0 * d1, d6, rtol=ftol)
        @test isapprox(2f0 * q3, q4, rtol=ftol)
        @test isapprox(2f0 * q3, q5, rtol=ftol)

        # Test linearity J * (a + b) == J * a + J * b
        dm2 = .5f0 * dm
        dm2.data .= circshift(dm2.data, (0, 20))
        lind =  J * dm
        lind2 = J * (2f0 .* dm)
        lind3 = J * dm2
        lind4 = J * (dm + dm2)
        lind5 = J * (dm - dm2)
        lind6 = (2f0 * J) * dm

        println("Test linearity of J: J * (a + b) == J * a + J * b")
        nd1 = norm(lind4)
        nd2 = norm(lind + lind3)
        nd3 = norm(lind4 - lind - lind3)/norm(lind4)
        @printf(" J * (a + b): %2.5e, J * a + J * b : %2.5e, relative error : %2.5e \n", nd1, nd2, nd3)
        nd1 = norm(lind5)
        nd2 = norm(lind - lind3)
        nd3 = norm(lind5 - lind + lind3)/norm(lind5)
        @printf(" J * (a - b): %2.5e, J * a - J * b : %2.5e, relative error : %2.5e \n", nd1, nd2, nd3)

        @test isapprox(lind4, lind + lind3, rtol=ftol)
        @test isapprox(lind5, lind - lind3, rtol=ftol)

        # Test linearity J a x == a J x
        println("Test linearity of J: J * (a * b) == a * J * b")
        nd1 = norm(2f0 * lind)
        nd1_b = norm(lind6)
        nd2 = norm(lind2)
        nd3 = norm(2f0 * lind - lind2)/norm(lind2)
        nd4 = norm(lind6 - lind2)/norm(lind2)

        # Adjoint J
        dma = adjoint(J) * d1
        dmb = adjoint(J) * (2f0 * d1)
        dmc = adjoint(J) * d2
        dmd = adjoint(J) * (d1 + d2)
        dme = adjoint(2f0 * J) * d1

        nm1 = norm(2f0 * dma)
        nm1_b = norm(dme)
        nm2 = norm(dmb)
        nm3 = norm(2f0*dma - dmb)/norm(dmb)
        nm4 = norm(dme - dmb)/norm(dmb)

        @printf(" a (J x): %2.5e, J a x : %2.5e, relative error : %2.5e \n", nd1, nd2, nd3)
        @printf(" (a J) x: %2.5e, J a x : %2.5e, relative error : %2.5e \n", nd1_b, nd2, nd4)

        @printf(" a (J' x): %2.5e, J' a x : %2.5e, relative error : %2.5e \n", nm1, nm2, nm3)
        @printf(" (a J') x: %2.5e, J' a x : %2.5e, relative error : %2.5e \n", nm1_b, nm2, nm4)

        @test isapprox(2f0 * lind, lind2, rtol=ftol)
        @test isapprox(2f0 * lind, lind6, rtol=ftol)
        @test isapprox(2f0 * dma, dmb, rtol=ftol)
        @test isapprox(2f0 * dma, dme, rtol=ftol)
    end
end

####################### Extended source operators ##########################################
if parsed_args["tti"]
    ftol = 5f-4
end

@testset "Extended source linearity test with $(nlayer) layers and tti $(tti) and freesurface $(fs)" begin
    @timeit TIMEROUTPUT "Extended source Linearity" begin
        # Modeling operators
        Pr = judiProjection(recGeometry)
        F = judiModeling(model; options=opt)
        Pw = judiLRWF(q1.geometry.dt[1], q1.data[1])

        A = Pr*F*adjoint(Pw)
        Aa = adjoint(A)
        # Extended source weights
        w = randn(Float32, model0.n)
        x = randn(Float32, model0.n)
        w[:, 1] .= 0f0; w = vec(w)
        x[:, 1] .= 0f0; x = vec(x)

        J = judiJacobian(A, w)

        d1 = A*w
        d2 = A*x
        d3 = A*(w+x)
        d4 = A*(w-x)
        d5 = A*(2f0 * w)
        d6 = (2f0 * A) * w

        q3 = Aa * d1
        q4 = Aa * (2f0 * d1)
        q5 = (2f0 * Aa) * d1

        dm2 =  .5f0 * dm
        dm2.data .= circshift(dm2.data, (0, 20))
        lind =  J * dm
        lind2 = J * (2f0 .* dm)
        lind3 = J * dm2
        lind4 = J * (dm + dm2)
        lind5 = J * (dm - dm2)
        lind6 = (2f0 * J) * dm

        dma = adjoint(J) * d1
        dmb = adjoint(J) * (2f0 * d1)
        dmc = adjoint(J) * d2
        dmd = adjoint(J) * (d1 + d2)
        dme = adjoint(2f0 * J) * d1

        # Test linearity F * (a + b) == F * a + F * b

        println("Test linearity of F: F * (a + b) == F * a + F * b")
        nd1 = norm(d3)
        nd2 = norm(d1 + d2)
        nd3 = norm(d3 - d1 - d2)/norm(d3)
        @printf(" F * (a + b): %2.5e, F * a + F * b : %2.5e, relative error : %2.5e \n", nd1, nd2, nd3)
        nd1 = norm(d4)
        nd2 = norm(d1 - d2)
        nd3 = norm(d4 - d1 + d2)/norm(d4)
        @printf(" F * (a - b): %2.5e, F * a - F * b : %2.5e, relative error : %2.5e \n", nd1, nd2, nd3)

        @test isapprox(d3, d1 + d2, rtol=ftol)
        @test isapprox(d4, d1 - d2, rtol=ftol)

        # Test linearity F a x == a F x

        println("Test linearity of F: F * (a * b) == a * F * b")
        nd1 = norm(2f0 * d1)
        nd1_b = norm(d6)
        nd2 = norm(d5)
        nd3 = norm(2f0 * d1 - d5)/norm(d5)
        nd4 = norm(d6 - d5)/norm(d5)

        nm1 = norm(2f0 * q3)
        nm1_b = norm(q5)
        nm2 = norm(q4)
        nm3 = norm(2f0 * q3 - q4)/norm(q4)
        nm4 = norm(q5 - q4)/norm(q4)
        @printf(" a (F x): %2.5e, F a x : %2.5e, relative error : %2.5e \n", nd1, nd2, nd3)
        @printf(" (a F) x: %2.5e, F a x : %2.5e, relative error : %2.5e \n", nd1_b, nd2, nd4)
        @printf(" a (F' x): %2.5e, F' a x : %2.5e, relative error : %2.5e \n", nm1, nm2, nm3)
        @printf(" (a F') x: %2.5e, F' a x : %2.5e, relative error : %2.5e \n", nm1_b, nm2, nm4)
        @test isapprox(2f0 * d1, d5, rtol=ftol)
        @test isapprox(2f0 * d1, d6, rtol=ftol)
        @test isapprox(2f0 * q3, q4, rtol=ftol)
        @test isapprox(2f0 * q3, q5, rtol=ftol)

        # Test linearity J * (a + b) == J * a + J * b

        println("Test linearity of J: J * (a + b) == J * a + J * b")
        nd1 = norm(lind4)
        nd2 = norm(lind + lind3)
        nd3 = norm(lind4 - lind - lind3)/norm(lind4)
        @printf(" J * (a + b): %2.5e, J * a + J * b : %2.5e, relative error : %2.5e \n", nd1, nd2, nd3)
        nd1 = norm(lind5)
        nd2 = norm(lind - lind3)
        nd3 = norm(lind5 - lind + lind3)/norm(lind5)
        @printf(" J * (a - b): %2.5e, J * a - J * b : %2.5e, relative error : %2.5e \n", nd1, nd2, nd3)

        @test isapprox(lind4, lind + lind3, rtol=ftol)
        @test isapprox(lind5, lind - lind3, rtol=ftol)

        # Test linearity J a x == a J x
        println("Test linearity of J: J * (a * b) == a * J * b")
        nd1 = norm(2f0 * lind)
        nd1_b = norm(lind6)
        nd2 = norm(lind2)
        nd3 = norm(2f0 * lind - lind2)/norm(lind2)
        nd4 = norm(lind6 - lind2)/norm(lind2)

        nm1 = norm(2f0 * dma)
        nm1_b = norm(dme)
        nm2 = norm(dmb)
        nm3 = norm(2f0*dma - dmb)/norm(dmb)
        nm4 = norm(dme - dmb)/norm(dmb)

        @printf(" a (J x): %2.5e, J a x : %2.5e, relative error : %2.5e \n", nd1, nd2, nd3)
        @printf(" (a J) x: %2.5e, J a x : %2.5e, relative error : %2.5e \n", nd1_b, nd2, nd4)

        @printf(" a (J' x): %2.5e, J' a x : %2.5e, relative error : %2.5e \n", nm1, nm2, nm3)
        @printf(" (a J') x: %2.5e, J' a x : %2.5e, relative error : %2.5e \n", nm1_b, nm2, nm4)

        @test isapprox(2f0 * lind, lind2, rtol=ftol)
        @test isapprox(2f0 * lind, lind6, rtol=ftol)
        @test isapprox(2f0 * dma, dmb, rtol=ftol)
        @test isapprox(2f0 * dma, dme, rtol=ftol)
    end
end


text:
# # Using MLJ to classifiy the MNIST image dataset

using Pkg
const DIR = @__DIR__
Pkg.activate(DIR)
Pkg.instantiate()

# **Julia version** is assumed to be 1.6.*

using MLJ
using Flux
import MLJFlux
using Random
Random.seed!(123)

MLJ.color_off()

using Plots
pyplot(size=(600, 300*(sqrt(5)-1)));

# ## Basic training

# Downloading the MNIST image dataset:

import MLDatasets: MNIST

ENV["DATADEPS_ALWAYS_ACCEPT"] = true
images, labels = MNIST.traindata();

# In MLJ, integers cannot be used for encoding categorical data, so we
# must force the labels to have the `Multiclass` [scientific
# type](https://alan-turing-institute.github.io/MLJScientificTypes.jl/dev/). For
# more on this, see [Working with Categorical
# Data](https://alan-turing-institute.github.io/MLJ.jl/dev/working_with_categorical_data/).

labels = coerce(labels, Multiclass);
images = coerce(images, GrayImage);

# Checking scientific types:

@assert scitype(images) <: AbstractVector{<:Image}
@assert scitype(labels) <: AbstractVector{<:Finite}

# Looks good.

# For general instructions on coercing image data, see [Type coercion
# for image
# data](https://alan-turing-institute.github.io/MLJScientificTypes.jl/dev/#Type-coercion-for-image-data-1)

images[1]

# We start by defining a suitable `Builder` object. This is a recipe
# for building the neural network. Our builder will work for images of
# any (constant) size, whether they be color or black and white (ie,
# single or multi-channel).  The architecture always consists of six
# alternating convolution and max-pool layers, and a final dense
# layer; the filter size and the number of channels after each
# convolution layer is customisable.

import MLJFlux
struct MyConvBuilder
    filter_size::Int
    channels1::Int
    channels2::Int
    channels3::Int
end

flatten(x::AbstractArray) = reshape(x, :, size(x)[end])
half(x) = div(x, 2)

function MLJFlux.build(b::MyConvBuilder, n_in, n_out, n_channels)

    k, c1, c2, c3 = b.filter_size, b.channels1, b.channels2, b.channels3

    mod(k, 2) == 1 || error("`filter_size` must be odd. ")

    p = div(k - 1, 2) # padding to preserve image size on convolution:

    h = n_in[1] |> half |> half |> half # final "image" height
    w = n_in[2] |> half |> half |> half # final "image" width

    return Chain(
        Conv((k, k), n_channels => c1, pad=(p, p), relu),
        MaxPool((2, 2)),
        Conv((k, k), c1 => c2, pad=(p, p), relu),
        MaxPool((2, 2)),
        Conv((k, k), c2 => c3, pad=(p, p), relu),
        MaxPool((2 ,2)),
        flatten,
        Dense(h*w*c3, n_out))
end

# **Note.** There is no final `softmax` here, as this is applied by
# default in all MLJFLux classifiers. Customisation of this behaviour
# is controlled using using the `finaliser` hyperparameter of the
# classifier.

# We now define the MLJ model. If you have a GPU, substitute
# `acceleration=CUDALibs()` below:

ImageClassifier = @load ImageClassifier
clf = ImageClassifier(builder=MyConvBuilder(3, 16, 32, 32),
                      acceleration=CPU1(),
                      batch_size=50,
                      epochs=10)

# You can add Flux options `optimiser=...` and `loss=...` here. At
# present, `loss` must be a Flux-compatible loss, not an MLJ measure.

# Binding the model with data in an MLJ machine:
mach = machine(clf, images, labels);

# Training for 10 epochs on the first 500 images:

fit!(mach, rows=1:500, verbosity=2);

# Inspecting:

report(mach)

#-

chain = fitted_params(mach)

#-

Flux.params(chain)[2]

#-

# Adding 20 more epochs:

clf.epochs = clf.epochs + 20
fit!(mach, rows=1:500);

# Computing an out-of-sample estimate of the loss:

predicted_labels = predict(mach, rows=501:1000);
cross_entropy(predicted_labels, labels[501:1000]) |> mean

# Or, in one line (after resetting the RNG seed to ensure the same
# result):

Random.seed!(123)
evaluate!(mach,
          resampling=Holdout(fraction_train=0.5),
          measure=cross_entropy,
          rows=1:1000,
          verbosity=0)

# ## Wrapping the MLJFlux model with iteration controls

# Any iterative MLJ model implementing the warm restart functionality
# illustrated above for `ImageClassifier` can be wrapped in *iteration
# controls*, as we demonstrate next. For more on MLJ's
# `IteratedModel` wrapper, see the [MLJ
# documentation](https://alan-turing-institute.github.io/MLJ.jl/dev/controlling_iterative_models/).

# The "self-iterating" model, called `imodel` below, is for iterating the
# image classifier defined above until one of the following stopping
# criterion apply:

# - `Patience(3)` (3 consecutive increases in the loss)

# - `InvalidValue()` (an out-of-sample loss, or a training loss,
#   is `NaN`, `Inf`, or `-Inf`)

# - `TimeLimit(t=1/60)` (training time has exceeded one minute)

# Additionally, training a machine bound to `imodel` will:

# - save a snapshot of the machine every three epochs

# - record the out-of-sample loss and training losses for plotting

# For a complete list of controls, see [this
# table](https://alan-turing-institute.github.io/MLJ.jl/dev/controlling_iterative_models/#Controls-provided).

losses = []
training_losses = [];

add_loss(loss) = push!(losses, loss)
add_training_loss(losses) = push!(training_losses, losses[end])

imodel = IteratedModel(model=clf,
                       controls=[Step(1),      # train one epoch at a time
                                 Patience(2),
                                 InvalidValue(),
                                 TimeLimit(0.5),
                                 Save(joinpath(DIR, "mnist_machine.jlso")),
                                 WithLossDo(), # for logging to `Info`
                                 WithLossDo(add_loss),
                                 WithTrainingLossesDo(add_training_loss)],
                       resampling=Holdout(fraction_train=0.7),
                       measure=log_loss,
                       retrain=false)

# Binding our self-iterating model to data:

mach = machine(imodel, images, labels)

# And training on the first 500 images:

fit!(mach, rows=1:500)

# A comparison of the training and out-of-sample losses:

plot(losses,
     title="Cross Entropy",
     xlab = "epoch",
     label="out-of-sample")
plot!(training_losses, label="training")

# Retrieving a snapshot for a prediction:

mach2 = machine(joinpath(DIR, "mnist_machine5.jlso"))
predict_mode(mach2, images[501:503])

#-

using Literate #src
Literate.markdown(@__FILE__, @__DIR__, execute=false) #src
Literate.notebook(@__FILE__, @__DIR__, execute=true) #src


text:
using Flux, Flux.Data.CMUDict
using Flux: onehot, batchseq
using Base.Iterators: partition

dict = cmudict()
alphabet = [:end, CMUDict.alphabet()...]
phones = [:start, :end, CMUDict.symbols()...]

tokenise(s, α) = [onehot(c, α) for c in s]

# Turn a word into a sequence of vectors
tokenise("PHYLOGENY", alphabet)
# Same for phoneme lists
tokenise(dict["PHYLOGENY"], phones)

words = sort(collect(keys(dict)), by = length)

# Finally, create iterators for our inputs and outputs.
batches(xs, p) = [batchseq(b, p) for b in partition(xs, 50)]

Xs = batches([tokenise(word, alphabet) for word in words],
             onehot(:end, alphabet))

Ys = batches([tokenise([dict[word]..., :end], phones) for word in words],
             onehot(:end, phones))

Yo = batches([tokenise([:start, dict[word]...], phones) for word in words],
             onehot(:end, phones))

data = collect(zip(Xs, Yo, Ys))

max_stars_repo_name: jstrube/model-zoo

max_issues_repo_name: jstrube/model-zoo

max_forks_repo_name: jstrube/model-zoo


text:
using Shannon
using Base.Test

v = collect(0:10) ./ 10

@test [1,1,2,3,4,5,6,7,8,9,10] == bin_vector(v, 0.0, 1.0, 10)

r = zeros((10,5))
i = 1
for j=1:10
  for k=1:5
    r[j,k]=i
    i = i+1
  end
end
v = r ./ 50
r = convert(Matrix{Int64}, r)

@test r == bin_matrix(v, 0.0, 1.0, 50)

c = [[1 0 0 0]; [1 2 0 0]; [1 2 3 0]; [1 2 3 4]]
r = [1; 1 + 2 * 4^1; 1 + 2*4^1 + 3*4^2; 1 + 2*4^1 + 3*4^2 + 4*4^3]

@test 1    == combine_binned_vector(squeeze(c[1,:],1), 10)
@test 21   == combine_binned_vector(squeeze(c[2,:],1), 10)
@test 321  == combine_binned_vector(squeeze(c[3,:],1), 10)
@test 4321 == combine_binned_vector(squeeze(c[4,:],1), 10)

@test r[1] == combine_binned_vector(squeeze(c[1,:],1), 4)
@test r[2] == combine_binned_vector(squeeze(c[2,:],1), 4)
@test r[3] == combine_binned_vector(squeeze(c[3,:],1), 4)
@test r[4] == combine_binned_vector(squeeze(c[4,:],1), 4)

@test r    == combine_binned_matrix(c)

# entropy tests
@test 0           == entropy([1,1,1,1,1])
@test abs(log(2,6) - entropy([1,2,3,4,5,6])) < 0.0000001
@test abs(log(4,6) - entropy([1,2,3,4,5,6], base=4)) < 0.0000001

#= @test abs(5.540741 - entropy([1,2,3,4,5,6], mode="ChaoShen")) < 0.00001 =#
#= @test abs(5.540741 - entropy([1,2,3,4,5,6], mode="ChaoShen", base=2)) < 0.00001 =#
#= @test abs(1.667929 - entropy([1,2,3,4,5,6], mode="ChaoShen", base=10)) < 0.00001 =#
#= @test abs(3.840549 - entropy([1,2,3,4,5,6], mode="ChaoShen", base=e)) < 0.00001 =#

# test basic functions - bin a value
@test 1  == bin_value(-1.0, 10)
@test 1  == bin_value(-0.9, 10)
@test 1  == bin_value(-0.8, 10)
@test 2  == bin_value(-0.7, 10)
@test 2  == bin_value(-0.6, 10)
@test 3  == bin_value(-0.5, 10)
@test 3  == bin_value(-0.4, 10)
@test 4  == bin_value(-0.3, 10)
@test 4  == bin_value(-0.2, 10)
@test 5  == bin_value(-0.1, 10)
@test 5  == bin_value(0.0,  10)
@test 6  == bin_value(0.1,  10)
@test 6  == bin_value(0.2,  10)
@test 7  == bin_value(0.3,  10)
@test 7  == bin_value(0.4,  10)
@test 8  == bin_value(0.5,  10)
@test 8  == bin_value(0.6,  10)
@test 9  == bin_value(0.7,  10)
@test 9  == bin_value(0.8,  10)
@test 10 == bin_value(0.9,  10)
@test 10 == bin_value(1.0,  10)

# test basic functions - bin a vector
@test [1, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10] == bin_vector([-1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0], -1.0, 1.0, 10)
@test [0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10] != bin_vector([-1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0], -1.0, 1.0, 10)

# test basic functions - bin a matrix
@test [[1, 2, 3] [4, 5, 6] [7, 8, 9]] == bin_matrix([[-1.0, -0.7, -0.5] [-0.3, -0.1, 0.1] [0.3, 0.5, 0.7]], -1.0, 1.0, 10)
@test [[2, 2, 3] [4, 5, 6] [7, 8, 9]] != bin_matrix([[-1.0, -0.7, -0.5] [-0.3, -0.1, 0.1] [0.3, 0.5, 0.7]], -1.0, 1.0, 10)

# test inverse binning
@test abs(-1.0 - unbin_value(1,  10, -1.0, 1.0, mode="lower")) < 0.0001
@test abs(-0.8 - unbin_value(2,  10, -1.0, 1.0, mode="lower")) < 0.0001
@test abs(-0.6 - unbin_value(3,  10, -1.0, 1.0, mode="lower")) < 0.0001
@test abs(-0.4 - unbin_value(4,  10, -1.0, 1.0, mode="lower")) < 0.0001
@test abs(-0.2 - unbin_value(5,  10, -1.0, 1.0, mode="lower")) < 0.0001
@test abs( 0.0 - unbin_value(6,  10, -1.0, 1.0, mode="lower")) < 0.0001
@test abs( 0.2 - unbin_value(7,  10, -1.0, 1.0, mode="lower")) < 0.0001
@test abs( 0.4 - unbin_value(8,  10, -1.0, 1.0, mode="lower")) < 0.0001
@test abs( 0.6 - unbin_value(9,  10, -1.0, 1.0, mode="lower")) < 0.0001
@test abs( 0.8 - unbin_value(10, 10, -1.0, 1.0, mode="lower")) < 0.0001

@test abs(-0.8 - unbin_value(1,  10, -1.0, 1.0, mode="upper")) < 0.0001
@test abs(-0.6 - unbin_value(2,  10, -1.0, 1.0, mode="upper")) < 0.0001
@test abs(-0.4 - unbin_value(3,  10, -1.0, 1.0, mode="upper")) < 0.0001
@test abs(-0.2 - unbin_value(4,  10, -1.0, 1.0, mode="upper")) < 0.0001
@test abs(-0.0 - unbin_value(5,  10, -1.0, 1.0, mode="upper")) < 0.0001
@test abs( 0.2 - unbin_value(6,  10, -1.0, 1.0, mode="upper")) < 0.0001
@test abs( 0.4 - unbin_value(7,  10, -1.0, 1.0, mode="upper")) < 0.0001
@test abs( 0.6 - unbin_value(8,  10, -1.0, 1.0, mode="upper")) < 0.0001
@test abs( 0.8 - unbin_value(9,  10, -1.0, 1.0, mode="upper")) < 0.0001
@test abs( 1.0 - unbin_value(10, 10, -1.0, 1.0, mode="upper")) < 0.0001

@test abs(-0.9 - unbin_value(1,  10, -1.0, 1.0, mode="centre")) < 0.0001
@test abs(-0.7 - unbin_value(2,  10, -1.0, 1.0, mode="centre")) < 0.0001
@test abs(-0.5 - unbin_value(3,  10, -1.0, 1.0, mode="centre")) < 0.0001
@test abs(-0.3 - unbin_value(4,  10, -1.0, 1.0, mode="centre")) < 0.0001
@test abs(-0.1 - unbin_value(5,  10, -1.0, 1.0, mode="centre")) < 0.0001
@test abs( 0.1 - unbin_value(6,  10, -1.0, 1.0, mode="centre")) < 0.0001
@test abs( 0.3 - unbin_value(7,  10, -1.0, 1.0, mode="centre")) < 0.0001
@test abs( 0.5 - unbin_value(8,  10, -1.0, 1.0, mode="centre")) < 0.0001
@test abs( 0.7 - unbin_value(9,  10, -1.0, 1.0, mode="centre")) < 0.0001
@test abs( 0.9 - unbin_value(10, 10, -1.0, 1.0, mode="centre")) < 0.0001


text:
# query_string = "SELECT a.* FROM appuser a
#                     WHERE a.login = \$1"
function retrieve_entity(filter_object::Union{IEntity,Missing},
                         data_type::DataType,
                         table_name::String,
                         columns_selection_and_mapping::Dict,
                         retrieve_complex_props::Bool,
                         dbconn::LibPQ.Connection)

    query_string = "SELECT o.* FROM " * table_name * " o "

    # Initialize the properties variables so that we can just test it's size
    #   later on in the function
    props = Dict{Symbol,Any}()

    # Instantiate an empty instance of the data type just to get the module
    orm_module = PostgresORM.get_orm(data_type())

    if !ismissing(filter_object)
        # Loop over the properties of the object and build the appropriate where clause
        props = util_get_entity_props_for_db_actions(filter_object,
                                                     dbconn,
                                                     false # do not include missing
                                                     )

        if length(props) > 0
            n = 0
            query_string *= " WHERE "
            for (propertyname,value) in props
                n += 1
                # Add the 'AND' if it's not the first restriction
                if n > 1
                    query_string *= " AND "
                end

                # NOTE: The mapping may not have the properties because some of the
                #          properties are the columns names of the complex objects
                #          which are added at runtime
                columnname = if haskey(columns_selection_and_mapping, propertyname)
                                columns_selection_and_mapping[propertyname]
                             else
                                 propertyname
                             end

                # Add the restriction
                query_string *= " o.$columnname = \$$n " # eg. 'a.login = \$1'

            end
        end
    end

    # prepared_query = LibPQ.prepare(dbconn,
    #                                query_string)

    # Prepare the query aruments
    # NOTE: It would be more elegant to convert the keys in an array of String
    #         but I didn't find a way to do it
    query_args = collect(values(props))

    result = execute_query_and_handle_result(query_string,
                                             data_type,
                                             query_args,
                                             columns_selection_and_mapping,
                                             retrieve_complex_props,
                                             dbconn)

    return(result)

end

function execute_query_and_handle_result(query_string::String,
                                         data_type::DataType,
                                         query_args::Union{Vector,Missing},
                                         retrieve_complex_props::Bool,
                                         dbconn::LibPQ.Connection)

    # Create a dummy object to get the ORM module
    orm_module = PostgresORM.get_orm(data_type())

    columns_selection_and_mapping = util_get_columns_selection_and_mapping(orm_module)

    execute_query_and_handle_result(query_string,
                                    data_type,
                                    query_args,
                                    columns_selection_and_mapping,
                                    retrieve_complex_props,
                                    dbconn)

end

function execute_plain_query(query_string::String,
                             query_args::Union{Vector,Missing},
                             dbconn::LibPQ.Connection)

   prepared_query = LibPQ.prepare(dbconn,
                                  query_string)

   if ismissing(query_args)
       query_args = []
   end

   query_result = execute(prepared_query,
                          query_args
                          ;throw_error=true)

    # Get the result as a dataframe
    data = DataFrame(query_result)

    # Close the result
    close(query_result)

    return(data)

end

function execute_query_and_handle_result(query_string::String,
                                         data_type::DataType,
                                         query_args::Union{Vector,Missing},
                                         columns_selection_and_mapping::Dict,
                                         retrieve_complex_props::Bool,
                                         dbconn::LibPQ.Connection)

     prepared_query = LibPQ.prepare(dbconn,
                                    query_string)

     if ismissing(query_args)
         query_args = []
     end

     query_result = execute(prepared_query,
                            query_args
                            ;throw_error=true)

      # Get the result as a dataframe
      data = DataFrame(query_result)

      # Close the result
      close(query_result)

      # Format the result to be mapped to the data model
      result = util_clean_and_rename_query_result(data,columns_selection_and_mapping)

      result = dataframe2vector_of_namedtuples(result)

      result =
         util_convert_namedtuple_to_object.(result,data_type,
                                        retrieve_complex_props,
                                        dbconn)
      return result
end

function retrieve_entity(filter_object::Union{IEntity,Missing},
                         orm_module::Module,
                         retrieve_complex_props::Bool,
                         dbconn::LibPQ.Connection)

     retrieve_entity(filter_object,
                     orm_module.data_type,
                     util_get_table_with_schema_name(orm_module),
                     util_get_columns_selection_and_mapping(orm_module),
                     retrieve_complex_props,
                     dbconn)
end

function retrieve_entity(filter_object::IEntity,
                         retrieve_complex_props::Bool,
                         dbconn::LibPQ.Connection)

         retrieve_entity(filter_object,
                         PostgresORM.get_orm(filter_object),
                         retrieve_complex_props,
                         dbconn)
end

function retrieve_one_entity(filter_object::T,
                             retrieve_complex_props::Bool,
                             dbconn::LibPQ.Connection) where T <: PostgresORM.IEntity

      results = retrieve_entity(filter_object,
                                retrieve_complex_props,
                                dbconn)

      if length(results) > 1
          error("Too many results for filter_object[$filter_object]")
      end
      if length(results) == 0
          return missing
      end
      result = results[1]
      return result
end


text:
using RecipesBase

"""
    VlasovSolution1D1V

Data structure that stores the solution of the Vlasov problem of one dimension
in both physical space and phase space.

"""
struct VlasovSolution1D1V

    values :: Vector{Array{Float64, 2}}
    times :: Vector{Float64}
    energy :: Vector{Float64}
    saveat :: Float64

    function VlasovSolution1D1V( )

        values = Array{Float64, 2}[]
        times = Float64[]
        energy = Float64[]
        new( values, times, energy, 0.0)

    end

    function VlasovSolution1D1V( saveat :: Float64 )

        values = Array{Float64, 2}[]
        energy = Float64[]
        times = Float64[]
        new( values, times, energy, saveat)

    end

end

@recipe function plot(sol::VlasovSolution1D1V)

    title --> "Vlasov problem 1D1V"
    xlabel --> "time"
    ylabel --> "energy"

    sol.times, sol.energy

end

max_stars_repo_name: JuliaVlasov/VlasovSolvers

max_issues_repo_name: JuliaVlasov/VlasovSolvers

max_forks_repo_name: JuliaVlasov/VlasovSolvers


text:
using Documenter, Xfoil

makedocs(;
    modules = [Xfoil],
    pages = [
        "Home" => "index.md",
        "Examples" => "examples.md",
        "Library" => "library.md"
    ],
    sitename = "Xfoil.jl",
    authors = "Taylor McDonnell <taylormcd@byu.edu>"
)

deploydocs(
    repo = "github.com/byuflowlab/Xfoil.jl.git"
)


text:
using SnoopCompile

### Log the compiles
# This only needs to be run once (to generate "/tmp/images_compiles.csv")

SnoopCompile.@snoop "http_compiles.csv" begin
    include(Pkg.dir("HTTP", "test","runtests.jl"))
end

### Parse the compiles and generate precompilation scripts
# This can be run repeatedly to tweak the scripts

# IMPORTANT: we must have the module(s) defined for the parcelation
# step, otherwise we will get no precompiles for the Images module
using HTTP

data = SnoopCompile.read("http_compiles.csv")

# Use these two lines if you want to create precompile functions for
# individual packages
pc, discards = SnoopCompile.parcel(data[end:-1:1,2])
SnoopCompile.write("src/precompile", pc)


text:
### Fetch required packages and reaction networks ###
using Catalyst, OrdinaryDiffEq, Random, SteadyStateDiffEq, Test
using ModelingToolkit: get_states, get_ps
include("test_networks.jl")

using StableRNGs
rng = StableRNG(12345)

### Compares to netowork with know steady state ###
steady_state_network_1 = @reaction_network begin
    (k1,k2), ∅ ↔ X1
    (k3,k4), ∅ ↔ 3X2
    (k5,k6), ∅ ↔ X3 + X4
end k1 k2 k3 k4 k5 k6

for factor in [1e-1, 1e0, 1e1], repeat = 1:3
    u0 = factor*rand(rng,length(get_states(steady_state_network_1))); u0[4] = u0[3];
    p = 0.01 .+ factor*rand(rng,length(get_ps(steady_state_network_1)))
    prob = SteadyStateProblem(steady_state_network_1,u0,p)
    sol = solve(prob,SSRootfind()).u
    (minimum(sol[1:1]) > 1e-2) && (@test abs.(sol[1] - p[1]/p[2]) < 0.01)
    (minimum(sol[2:2]) > 1e-2) && (@test abs.(sol[2]^3/factorial(3) - p[3]/p[4]) < 0.01)
    (minimum(sol[3:4]) > 1e-2) && (@test abs.(sol[3]*sol[4] - p[5]/p[6]) < 0.01)
end

steady_state_network_2 = @reaction_network begin
    v/10+hill(X,v,K,n), ∅ → X
    d, X → ∅
end v K n d

for factor in [1e-1, 1e1, 1e1], repeat = 1:3
    u0_small = factor*rand(rng,length(get_states(steady_state_network_2)))/100
    u0_large = factor*rand(rng,length(get_states(steady_state_network_2)))*100
    p = factor*rand(rng,length(get_ps(steady_state_network_2)))
    p[3] = round(p[3])+1
    sol1 = solve(SteadyStateProblem(steady_state_network_2,u0_small,p),SSRootfind()).u[1]
    sol2 = solve(SteadyStateProblem(steady_state_network_2,u0_large,p),SSRootfind()).u[1]
    diff1 = abs(p[1]/10 + p[1]*(sol1^p[3])/(sol1^p[3]+p[2]^p[3]) - p[4]*sol1)
    diff2 = abs(p[1]/10 + p[1]*(sol2^p[3])/(sol2^p[3]+p[2]^p[3]) - p[4]*sol2)
    @test (diff1 < 1e-8) || (diff2 < 1e-8)
end

### For a couple of networks, test that the steady state solution is identical to the long term ODE solution. ###
steady_state_test_networks = [reaction_networks_standard[8], reaction_networks_standard[10], reaction_networks_weird[1]]
for network in steady_state_test_networks, factor in [1e-1, 1e0, 1e1]
    u0 = factor*rand(rng,length(get_states(network)))
    p = factor*rand(rng,length(get_ps(network)))
    sol_ode = solve(ODEProblem(network,u0,(0.,1000000),p),Rosenbrock23())
    sol_ss = solve(SteadyStateProblem(network,u0,p),SSRootfind())

    @test all(abs.(sol_ode.u[end] .- sol_ss.u) .< 1e-4)
end

### No parameter test ###

no_param_network = @reaction_network begin
    (0.6,3.2), ∅ ↔ X
end
for factor in [1e0, 1e1, 1e2]
    u0 = factor*rand(rng,length(get_states(no_param_network)))
    sol_ss = solve(SteadyStateProblem(no_param_network,u0),SSRootfind(),abstol=1e-11)
    @test abs.(sol_ss.u[1]-0.6/3.2) < 1e-8
end