text:
module Godunov
using LinearAlgebra
#using StaticArrays
"""
    For now, we assume a constant (fixed) viscosity "ν".
"""

TINY = 0.5 # how much to shink the timestep

α = 1.0
ν = 1.0
s = 0.1 # relaxation time
q = 3*ν/s
xv = (α/(3ν))^(1.0/3.0)
xh = (α/√q)^0.25
print("xv: $xv\n")
print("xh: $xh\n")
# Boundary conditions (Dirichlet):
Σ_lbc = 0.0
Σ_rbc = 1.0
Φ_lbc = 0.0
Φ_rbc = 0.0#-0.4

"""
A planetary particulate disk is called a "ring", but this terminology is confusing.
"""
mutable struct ring
    N::Integer              # number of cells
    Δt::Float64             # time step
    X_::Array{Float64,1}    # positions (nodes)
    ΔX_::Array{Float64,1}   # cell sizes (cell)
    U_::Array{Float64,1}    # velocity (on nodes)
    ∇U_::Array{Float64,1}   # divergence of velocity
    Σ_::Array{Float64,1}    # surface density
    S_::Array{Float64,1}    # integral of the above
    Φ_::Array{Float64,1}    # mass flux
    F_::Array{Float64,1}    # integral of the above
    r₊::Array{Float64,1}    # 2-component eigenvector, right-travelling wave (relative to fluid) [NB: independent of position]
    r₋::Array{Float64,1}    # 2-component eigenvector, left-travelling wave (relative to fluid) [NB: independent of position]
    λ₊::Array{Float64,1}    # eigenvalues, right-travelling (relative to fluid) wave
    λ₋::Array{Float64,1}    # eigenvalues, left-travelling (relative to fluid) wave
#    Λ₊R::Any                # diagonal matrix of eigenvalues, Right-travelling (absolute) right-travelling (relative)
#    Λ₊L::Any                # diagonal matrix of eigenvalues, Left-travelling (absolute) right-travelling (relative)
#    Λ₋R::Any                #
#    Λ₋L::Any                #
    α⁺::Array{Float64,1}    #
    α⁻::Array{Float64,1}    #
    t::Float64              # time
    i::Integer              # iteration counter
end
function ring(N=5;xi=(0.5*min(xv,xh)),xf=(3.0*max(xv,xh)))
    function setEs(U_::Array{Float64,1})
        n = length(U_)
        r₊ = [1.0, √q]
        r₋ = [1.0,-√q]
        λ₊ = zeros(n)
        λ₋ = zeros(n)
        for i in 1:n
            λ₊[i] = U_[i] + √q
            λ₋[i] = U_[i] - √q
        end
#=        Λ₊R =  Diagonal(λ₊ .* (λ₊ .> 0))
        Λ₊L = -Diagonal(λ₊ .* (λ₊ .< 0))
        Λ₋R =  Diagonal(λ₋ .* (λ₋ .> 0))
        Λ₋L = -Diagonal(λ₋ .* (λ₋ .< 0)) =#
        return (r₊,r₋,λ₊,λ₋)#,Λ₊R,Λ₊L,Λ₋R,Λ₋L)
    end

    x = 5.0; y = 1.0/x
    X_ = collect(range(xi^y,xf^y,length=N+1)).^x
    ΔX_ = diff(X_)
    U_ = α ./ X_.^4
    U_ .-= U_[end]
    #U_ .= -1.5
    ∇U_ = diff(U_) ./ ΔX_

    #Σ_ = ones(N); Σ_[1] = 0.0; #Σ_[N÷2:end].=1.0
    Σ_ = (X_[2:end].-xh)./xf
    Σ_ .*= (Σ_.>0)
    Φ_ = Σ_ .* (U_[1:end-1]+U_[2:end])/2.0  .* 0.0

    S_ = Σ_ .* ΔX_
    S_ ./= sum(S_)
    Σ_ = S_ ./ ΔX_
    F_ = Φ_ .* ΔX_

    (r₊,r₋,λ₊,λ₋) = setEs(U_)

    Δt = 0.5 / max( maximum(abs.(λ₊[1:end-1]./ΔX_)), maximum(abs.(λ₋[2:end]./ΔX_)), 1.0/s, maximum(abs.(∇U_))) * TINY
    print("Δt: $Δt\n")
    α⁺=zeros(N+1)
    α⁻=zeros(N+1)
    r = ring(N,Δt,X_,ΔX_,U_,∇U_,Σ_,S_,Φ_,F_,r₊,r₋,λ₊,λ₋,α⁺,α⁻,0.0,0)
    getAlphas!(r)
    return r
end

function getAlphas!(r::ring)
    N = r.N; Σ_ = r.Σ_; Φ_ = r.Φ_; α⁺ = r.α⁺; α⁻ = r.α⁻
    ΔΣ_ = cat(Σ_[1]-Σ_lbc, diff(Σ_),        0.0, dims=1) # reflective bc on rt
    ΔΦ_ = cat(Φ_[1]-Φ_lbc, diff(Φ_), -2*Φ_[end], dims=1) # reflective bc w/sign flip
    for i in 1:(N+1)
        α⁺[i] = 0.5 * [1  1/√q] ⋅ [ΔΣ_[i], ΔΦ_[i]]
        α⁻[i] = 0.5 * [1 -1/√q] ⋅ [ΔΣ_[i], ΔΦ_[i]]
    end
    r.α⁺ = α⁺
    r.α⁻ = α⁻
end#function

function step!(r::ring; Δt=r.Δt)
    N=r.N; Σ_=r.Σ_; ΔX_=r.ΔX_; S_=r.S_; Φ_=r.Φ_; F_=r.F_; α⁺=r.α⁺; α⁻=r.α⁻; r₋=r.r₋; r₊=r.r₊
    λ₊=r.λ₊; λ₋=r.λ₋
    getAlphas!(r)

    S_ = Σ_ .* ΔX_
    F_ = Φ_ .* ΔX_
    """ LHS: """
    # right-traveling.
    Σ_ -= r.α⁺[1:end-1] .* Δt .* r₊[1] .* λ₊[1:end-1] .* (λ₊[1:end-1].>0)
    Σ_ -= r.α⁺[1:end-1] .* Δt .* r₋[1] .* λ₋[1:end-1] .* (λ₋[1:end-1].>0)
    Φ_ -= r.α⁺[1:end-1] .* Δt .* r₊[2] .* λ₊[1:end-1] .* (λ₊[1:end-1].>0)
    Φ_ -= r.α⁺[1:end-1] .* Δt .* r₋[2] .* λ₋[1:end-1] .* (λ₋[1:end-1].>0)
    # left-traveling
    Σ_ -= r.α⁻[2:end]   .* Δt .* r₊[1] .* λ₊[2:end] .* (λ₊[2:end].<0)
    Σ_ -= r.α⁻[2:end]   .* Δt .* r₋[1] .* λ₋[2:end] .* (λ₋[2:end].<0)
    Φ_ -= r.α⁻[2:end]   .* Δt .* r₊[2] .* λ₊[2:end] .* (λ₊[2:end].<0)
    Φ_ -= r.α⁻[2:end]   .* Δt .* r₋[2] .* λ₋[2:end] .* (λ₋[2:end].<0)
    """ RHS: """
    # This part kinda sucks b/c it doesn't conserve very well:
    # Σ_ = S_ ./ ΔX_;  Φ_ = F_ ./ ΔX_
    Σ_[1] = Σ_lbc
    Σ_[end] = Σ_rbc
    Φ_[1] = Φ_lbc
    Φ_[end] = Φ_rbc

    #Σ_ += -Δt * r.∇U_ .* Σ_ # <- forward Euler
    Σ_ .*= (1 .- Δt*r.∇U_)
    #Σ_ ./= (1 .+ Δt*r.∇U_)  # <- backward Euler  - conserves a bit better
    #Σ_ .*= (1 .- Δt*r.∇U_)
    #Φ_ ./= (1 .+ Δt*r.∇U_ .+ Δt/s)
    Φ_ ./= (1 .+ Δt/s)
#=
    for i in 1:N
        if r.X_[i]<xh
            Σ_[i] = Φ_[i] = 0
        end
    end
=#
    #

    #S_ *= 0
    #
    # massaging.
    # 1. make sign uniform
#    Σ_ = abs.(S_); Φ_ = -abs.(Φ_)
#=    # 2. make monotonic
    for i in 2:N
        Σ_[i] = max(Σ_[i-1],Σ_[i])
#        Φ_[i] = min(Φ_[i-1],Φ_[i])
    end
#    Σ_ /= maximum(Σ_)
# =#
#    Σ_[7] = 0.0

    S_ = Σ_ .* ΔX_; F_ = Φ_ .* ΔX_
    r.S_ = S_ ./ sum(S_)
    r.F_ = F_ ./ sum(S_) # yes, sum(S_), not sum(F_)
    r.Σ_ = r.S_ ./ ΔX_
    r.Φ_ = r.F_ ./ ΔX_

    r.t += Δt
    r.i += 1
    s1 = sum(S_); s2 = sum(F_)
    #print("Sums: $s1,  $s2 \n")
    return (sum(S_),sum(F_))
end

function tstep!(r::ring,Δt)
    M = floor(Δt / r.Δt)
    Δt_last = Δt - M*r.Δt
    @assert Δt_last ≥ 0
    for i in 1:M
        step!(r)
    end
    step!(r, Δt = Δt_last)
end

end#module

max_stars_repo_name: diakosmos/razor

max_issues_repo_name: diakosmos/razor

max_forks_repo_name: diakosmos/razor


text:
using Catalyst
rn = @reaction_network begin
    α, S + I --> 2I
    β, I --> R
    S^2, R --> 0
end α β

# check can make a graph
gr = Graph(rn)

# check can save a graph
fname = Base.Filesystem.tempname()
savegraph(gr, fname, "png")

ext: jl


text:
using LTspice
using Test
using Dates

function compare_arrays(x,y)
  @test length(x) == length(y)
  for i in eachindex(x)
    if !isnan(x[i])
      @test x[i]≈y[i]
    end
  end
  return nothing
end

@testset "tests not calling LTspice.exe" begin
  include("test1.jl")
  include("test2.jl")
  include("test3.jl")
  include("test4.jl")
  include("test5.jl")
  include("test6.jl")
  include("test7.jl")
  include("test8.jl")
  include("test9.jl")
  include("test10.jl")
  include("test11.jl")
  include("test12.jl")
  include("test13.jl")
  include("test14.jl")
  include("test15.jl")
  include("testinc.jl")
end

is_ltspice_installed = (try LTspice.defaultltspiceexecutable() catch nothing end)!=nothing

if is_ltspice_installed
  @testset "tests calling LTspice.exe" begin
    include("localtests.jl")
  end
else
  println("LTspice.exe not found.  Skipping tests.")
end


text:
@testset "CUDA: BOCcos" begin
    system = BOCcos(GPSL1(use_gpu = Val(true)), 1, 1)
    @test @inferred(Tracking.is_upcoming_integration_new_bit(system, 0xfffff00000, 50)) == true
    @test @inferred(Tracking.is_upcoming_integration_new_bit(system, 0xfffff, 10)) == false
    @test @inferred(Tracking.is_upcoming_integration_new_bit(system, 0xfffff, 40)) == true

    @test @inferred(Tracking.get_default_correlator(system, NumAnts(1))) == EarlyPromptLateCorrelator(NumAnts(1))
    @test @inferred(Tracking.get_default_correlator(system, NumAnts(3))) == EarlyPromptLateCorrelator(NumAnts(3))

    system = BOCcos(GPSL5(use_gpu = Val(true)), 1, 1)
    @test @inferred(Tracking.is_upcoming_integration_new_bit(system, 0x35, 50)) == true
    @test @inferred(Tracking.is_upcoming_integration_new_bit(system, 0x35, 5)) == false
    @test @inferred(Tracking.is_upcoming_integration_new_bit(system, 0x3ca, 10)) == true # 0x3ca == 1111001010

    @test @inferred(Tracking.get_default_correlator(system, NumAnts(1))) == EarlyPromptLateCorrelator(NumAnts(1))
    @test @inferred(Tracking.get_default_correlator(system, NumAnts(3))) == EarlyPromptLateCorrelator(NumAnts(3))

end


text:
# To-dos:
#
# - How do we handle TF/SF?  Need a capability to do subpixel smoothing only inside some
# region.
#
# - A symmetry boundary puts the mirror image of the boundary-touching object in the space
# behind the boundary.  When this happens, not only the shape but also the material of the
# object must be flipped with respect to the symmetry boundary.  Such material flip is
# currently unsupported: only the Bloch boundary condition is valid for the boundary
# touching anisotropic materials.

export smooth_param!

# Given a vector `v` of comparable objects and vector `ind` of indices that sort `v`, return
# the number `c` of different objects in `v` and the index `nₑ` of `ind` where the last new
# object occurs (i.e., v[ind[nₑ]] is the last new object).
#
# If v is uniform, nₑ = length(v) + 1 (i.e., it is out of bound).
#
# For efficiency, this function assumes `v` is a vector of objects assigned to a voxel, i.e.,
# length(v) = 8.
@inline function countdiff(ind::MArray{TP,<:Integer,K,Zᴷ}, v::MArray{TP,<:Real,K,Zᴷ}) where {TP,K,Zᴷ}
    # Zᴷ = number of entries in MArray
    c, nₑ = 1, Zᴷ+1
    @simd for n = 2:Zᴷ
        @inbounds v[ind[n-1]] ≠ v[ind[n]] && (c += 1; nₑ = n)
    end
    return c, nₑ  # c: number of changes in v; nₑ: last n where v[ind[n]] changed
end

# vectors from voxel corners to the voxel center
nout_vxl(::Val{3}) = (SVec(1.,1.,1.), SVec(-1.,1.,1.), SVec(1.,-1.,1.), SVec(-1.,-1.,1.),
                      SVec(1.,1.,-1.), SVec(-1.,1.,-1.), SVec(1.,-1.,-1.), SVec(-1.,-1.,-1.))
nout_vxl(::Val{2}) = (SVec(1.,1.), SVec(-1.,1.), SVec(1.,-1.), SVec(-1.,-1.))
nout_vxl(::Val{1}) = (SVec(1.), SVec(-1.))

# Overall smoothing algorithm
#
# Below, obj, pind, oind refers to an Object{3} instance, parameter index (integer value) that
# distinguishes different material parameters, and oind that distinguishes different objects.
#
# - Assign oind, pind, shape to arrays object-by-object (see assignment.jl).
#     - Currently, we assign only oind (in oindKd′), and pind and shape are retrieved from oind by oind2shp and oind2pind.
#     - For the locations of grid points to assign the objects to, use τlcmp (lcmp created considering BC).
# - Using oind and pind, determine voxels to perform subpixel smoothing.
# - Inside each voxel, figure out the foreground object with which subpixel smoothing is performed.
#     - Iterating voxel corners, find the object that is put latest.  That object is the foreground object.
#     - Complication occurs when the voxel corner assigned with the foreground object is outside the domain boundary.
#         - In that case, to calculate nout and rvol properly, we have to move the voxel center.  (This is effectively the same as moving the foreground object.)
#             - If the corner is outside the periodic boundary, translate the voxel center before calling surfpt_nearby (using ∆fg).
#             - If the corner is outside the symmetry boundary, zero the nout component normal to the boundary (using σvxl).
#             - ∆fg and σvxl can be obtained simply by looking at the indices and BC.

# Below, the primed quantities are for voxel corner quantities rather than voxel center
# quantities.  Compare this with the usage of the prime in assignment.el that indicates
# complementary material (i.e., magnetic for electric).
function smooth_param!(paramKd::AbsArrComplexF{K₊₂},  # parameter array to smooth at voxel centers
                       oindKd′::NTuple{Kf⏐₁,AbsArr{ObjInd,K}},  # object index array at voxel corners; does not change
                       oind2shp::AbsVec{Shape{K,K²}},  # input map from oind to shape
                       oind2pind::AbsVec{ParamInd},  # input map from oind to pind
                       pind2matprm::AbsVec{S²ComplexF{Kf,Kf²}},  # map from pind to electric (magnetic) material parameters; Kf² = Kf^2
                       gt₀::SVec{K,GridType},  # grid type of voxel corners; generated by ft2gt.(ft, boundft)
                       l::Tuple2{NTuple{K,AbsVecReal}},  # grid locations excluding ghost points (length(l[PRIM][k]) = N[k]); used for locations of field components
                       lg::Tuple2{NTuple{K,AbsVecReal}},  # grid locations including ghost points (length(lg[PRIM][k]) = N[k]+1); used for locations of voxel corners without transformation by boundary conditions
                       σ::Tuple2{NTuple{K,AbsVecBool}},  # false if on symmetry boundary; exclude ghost points (length(σ[PRIM][k]) = N[k])
                       ∆τ::Tuple2{NTuple{K,AbsVecReal}},  # amount of shift by Bloch boundary conditions; include ghost points (length(∆τ[PRIM][k]) = N[k]+1)
                       isfield˔shp::Bool=isfield_ortho_shape(Kf,K)  # true if spaces where field and shapes are defined are orthogonal complement of each other; false if they are the same
                       ) where {K,K²,K₊₂,Kf,Kf²,Kf⏐₁}
    @assert K²==K^2 && Kf²==Kf^2 && K₊₂==K+2 && (Kf⏐₁==Kf || Kf⏐₁==1)
    @assert size(paramKd,K+1)==size(paramKd,K+2)==Kf

    ci_1 = CartesianIndex(ntuple(x->1, Val(K)))
    for nw = 1:Kf⏐₁
        # Set the grid types of the x-, y-, z-locations of Fw.
        #
        # I thought I could replace the condition Kf⏐₁==1 with isfield_ortho_shape(Kf,K),
        # because gt_cmp becomes gt₀ when isfield_ortho_shape(Kf,K) = true.  The replacement
        # is incorrect, however, because we need gt_cmp to become gt₀ even when the shape
        # and field subspaces coincide (e.g., K = Kf = 3) to assign the off-diagonal entries
        # of the material parameter tensors at voxel corners.
        gt_cmp = Kf⏐₁==1 ? gt₀ : gt_w(nw, gt₀)

        # Set the grid types of the voxel corners surroundnig Fw.
        gt_cmp′ = alter.(gt_cmp)

        # Choose various vectors for Fw (which is at the centers of the voxel defined by the
        # voxel corners below).
        lcmp = t_ind(l, gt_cmp)
        σcmp = t_ind(σ, gt_cmp)

        # Choose various arrays and vectors for voxel corners.
        lcmp′ = t_ind(lg, gt_cmp′)
        ∆τcmp′ = t_ind(∆τ, gt_cmp′)

        oind_cmp′ = oindKd′[nw]

        ind_c = MArray{Tuple{ntuple(x->2,Val(K))...},Int}(undef)  # corner indices inside voxel
        pind_vxl′ = similar(ind_c)  # material parameter indices inside voxel
        oind_vxl′ = similar(ind_c)  # object indices inside voxel

        # Set various arrays for the current component.
        CI = CartesianIndices(length.(lcmp))
        for ci_cmp = CI
            # Below, ci_1 is always +1 (i.e., not -1 depending on the grid type of field
            # component), because the field component locations lcmp exclude the ghost
            # points whereas the voxel corner locations lcmp′ include the ghost points.
            ci_vxl′ = (ci_cmp, ci_cmp + ci_1)  # Tuple2{SInt{3}}

            if !is_vxl_uniform!(oind_vxl′, pind_vxl′, ci_vxl′, oind_cmp′, oind2shp, oind2pind)
                # Note that pind_vxl′ is initialized at this point.
                sort_ind!(ind_c, pind_vxl′)  # pind_vxl′[ind_c[n]] ≤ pind_vxl′[ind_c[n+1]]
                Nprm_vxl, n_diffp = countdiff(ind_c, pind_vxl′)  # n_diffp: last n where pind_vxl′[ind_c[n]] changed

                if Nprm_vxl ≥ 2
                    if Nprm_vxl == 2
                        # Attempt to apply Kottke's subpixel smoothing algorithm.
                        prm_vxl = smooth_param_vxl(ci_vxl′, ind_c, n_diffp, pind_vxl′, oind_vxl′,
                                                   oind_cmp′, oind2shp, oind2pind, pind2matprm,
                                                   lcmp, lcmp′, σcmp, ∆τcmp′, isfield˔shp)
                    else  # Nprm_vxl ≥ 3
                        # Give up Kottke's subpixel smoothing and take simple averaging.
                        # See the comments at the next usage of amean_param().
                        prm_vxl = amean_param(ci_vxl′, oind_cmp′, oind2pind, pind2matprm)::S²ComplexF{Kf,Kf²}
                    end  # if Nprm_vxl == 2

                    # Below, we have four combinations depending on the values of Kf⏐₁ and Kf:
                    # Kf⏐₁ is either ==Kf or ==1, and Kf is either ≥2 or ==1.
                    #
                    # The case of Kf⏐₁ == Kf covers Kf≥2 and Kf==1.  When Kf==1, the case
                    # corresponds to Kf⏐₁==1 and Kf==1.  Therefore, the "if" case Kf⏐₁ == Kf
                    # covers the two combinations:
                    # - Kf⏐₁==Kf and Kf≥2,
                    # - Kf⏐₁==1 and Kf==1.
                    #
                    # Therefore, the "else" case should cover the remaining two combinations:
                    # - Kf⏐₁==Kf && Kf==1,
                    # - Kf⏐₁==1 && Kf≥2.
                    # However, the first combination is equivalent to Kf⏐₁==1 && Kf==1, which
                    # was already covered in the "if" case.  Therefore, the "else" case
                    # deals with only the case of Kf⏐₁==1 && Kf≥2.
                    if Kf⏐₁ == Kf  # Kf⏐₁==Kf && Kf≥2, or Kf⏐₁==1 && Kf==1
                        # Overwrite the (nw,nw)-diagonal entry of paramKd.
                        @inbounds paramKd[ci_cmp,nw,nw] = prm_vxl[nw,nw]
                    else  # Kf⏐₁==1 && Kf≥2
                        # Overwrite the off-diagonal entires of paramKd.
                        # Below the main diagonal.
                        for c = 1:Kf, r = c+1:Kf  # column- and row-indices
                            @inbounds paramKd[ci_cmp,r,c] = prm_vxl[r,c]
                        end

                        # Above the main diagonal.
                        for c = 2:Kf, r = 1:c-1  # column- and row-indices
                            @inbounds paramKd[ci_cmp,r,c] = prm_vxl[r,c]
                        end
                    end  # if Kf⏐₁ == 1
                end  # if Nprm_vxl ≥ 2
                # If Nprm_vxl = 1, the voxel is filled with the same material, and there is
                # no need for subpixel smoothing.
            end  # if !is_vxl_uniform!
        end  # for ci_cmp = CI
    end  # for nw = 1:Kf+1

    return nothing
end

function is_vxl_uniform!(oind_vxl′::MArray{TP,Int,K,Zᴷ},  # scratch vector to store object indices inside voxel; TP = Tuple{2,...,2}, Zᴷ = 2^K
                         pind_vxl′::MArray{TP,Int,K,Zᴷ},  # scratch vector to store material parameter indices inside voxel; TP = Tuple{2,...,2}, Zᴷ = 2^K
                         ci_vxl′::Tuple2{CartesianIndex{K}},
                         oind_cmp′::AbsArr{ObjInd,K},  # object index array; does not change
                         oind2shp::AbsVec{Shape{K,K²}},  # input map from oind to shape
                         oind2pind::AbsVec{ParamInd}  # input map from oind to pind
                         ) where {TP,K,K²,Zᴷ}
    # Retrieve the elements assigned to voxel corners from 3D arrays.
    #
    # Unlike the rest of the code that is performed only where subpixel smoothing is
    # performed, the following for loop is performed at all grid points.  Therefore, the
    # number of assignments to perform in the following for loop can easily reach
    # several millions even for a simple 3D problem.  Simple assignment can
    # be surprisingly time-consuming when peformed such a large number of times.  The
    # use of @inbounds helps reducing assignment time significantly.
    nc = 0
    for ci = ci_vxl′[1]:ci_vxl′[2]
        nc += 1
        @inbounds oind_vxl′[nc] = oind_cmp′[ci]
    end

    is_vxl_uniform = isuniform(oind_vxl′)

    # Even if the voxel is filled with multiple objects, those objects could be made of
    # the same material.  Then, the voxel is still filled with uniform material
    if !is_vxl_uniform
        for nc = 1:Zᴷ
            @inbounds pind_vxl′[nc] = oind2pind[oind_vxl′[nc]]
        end

        # We could use Nprm_vxl calculated below instead of is_vxl_uniform: if
        # Nprm_vxl = 1, then is_vxl_uniform = true.  However, calculating Nprm_vxl
        # requires calling sort8! and countdiff, which are costly when called for all
        # voxels.  Therefore, by using is_vxl_uniform, we avoid calling those functions
        # when unnecessary.
        @inbounds is_vxl_uniform = isuniform(pind_vxl′)
    end

    # Note that when is_vxl_uniform = false at this point, pind_vxl′ is initialized.
    return is_vxl_uniform
end

# Calculate prm_vxl (averaged material parameter for Fw inside a voxel) when the voxel is
# composed of two material parameters (Nprm_vxl = 2).  This includes cases where the voxel
# is composed of two material parameters but of ≥3 shape (i.e., ≥2 shapes in the voxel have
# the same material).
#
# Below, XXX_cmp has size N, whereas XXX_cmp′ has size N+1 (and corresponds to voxel corners).
# Variable names ending with vxl store quantities within a voxel.  Primed names store values
# at voxel corners; nonprimed names store values representing a voxel.
function smooth_param_vxl(ci_vxl′::Tuple2{CartesianIndex{K}},
                          ind_c::MArray{TP,Int,K,Zᴷ},  # scratch vector to store corner indices inside voxel; TP = Tuple{2,...,2}, Zᴷ = 2^K
                          n_diffp::Int,  # last n where pind_vxl′[ind_c[n]] changed
                          pind_vxl′::MArray{TP,Int,K,Zᴷ},  # scratch vector to store material parameter indices inside voxel; TP = Tuple{2,...,2}, Zᴷ = 2^K
                          oind_vxl′::MArray{TP,Int,K,Zᴷ},  # scratch vector to store object indices inside voxel; TP = Tuple{2,...,2}, Zᴷ = 2^K
                          oind_cmp′::AbsArr{ObjInd,K},  # object index array; does not change
                          oind2shp::AbsVec{Shape{K,K²}},  # input map from oind to shape
                          oind2pind::AbsVec{ParamInd},  # input map from oind to pind
                          pind2matprm::AbsVec{S²ComplexF{Kf,Kf²}},  # map from pind to electric (magnetic) material parameters; Kf² = Kf^2
                          lcmp::NTuple{K,AbsVecReal},  # location of field components
                          lcmp′::NTuple{K,AbsVecReal},  # location of voxel corners without transformation by boundary conditions
                          σcmp::NTuple{K,AbsVecBool},  # false if on symmetry boundary
                          ∆τcmp′::NTuple{K,AbsVecReal},  # amount of shift into domain by Bloch boundary conditions
                          isfield˔shp::Bool=false  # true if spaces where material parameters and shapes are defined are orthogonal complement of each other; false if they are the same
                          ) where {TP,K,K²,Kf,Kf²,Zᴷ}
    # First, attempt to apply Kottke's subpixel smoothing algorithm.
    nout = @SVector zeros(K)
    rvol = 0.0

    with2objs = true  # will be updated, because there could be ≥3 objects in voxel even if Nprm_vxl = 2
    @inbounds ind_c1, ind_c2 = ind_c[Zᴷ], ind_c[1]::Int  # indices of corners of two different material parameters
    @inbounds oind_c1, oind_c2 = oind_vxl′[ind_c1], oind_vxl′[ind_c2]

    # The corners ind_c[n_diffp:8] are occupied with one material parameter because
    # ind_c is sorted for pind_vxl′.  However, those corners can still be occupied with
    # ≥2 objects.  In that case, the voxel is composed of ≥3 objects.
    #
    # Note that not having two objects inside the voxel can logically mean having either
    # one object or ≥3 objects, but because Nprm_vxl ≥ 2, with2objs = false only means
    # having ≥3 objects inside the voxel.
    for nc = Zᴷ-1:-1:n_diffp  # n = 2^K-1 corresponds to oind_c1 and is omitted
        @inbounds (with2objs = oind_vxl′[ind_c[nc]]==oind_c1) || break
    end

    # At this point, with2objs = false means the corners ind_c[n_diffp:8] are composed
    # of ≥2 objects.  Then, we already know the voxel has ≥3 objects because the corners
    # ind_c[1:n_diffp-1] have different materials than the corners ind_c[n_diffp:8]
    # because Nprm_vxl = 2 and therefore have a different object than those occupying
    # the corners ind_c[n_diffp:8].  Hence, if with2objs = false at this point, we don't
    # have to test further if the voxel has ≥3 objects: it already has.
    if with2objs  # single object for corners ind_c[n_diffp:8]
        # The corners ind_c[1:n_diffp-1] are occupied with one material parameter
        # because ind_c is sorted for pind_vxl′, but those corners can still be occupied
        # with more than one object.  In that case, the voxel is composed of more than
        # two objects.
        for nc = 2:n_diffp-1  # n = 1 corresponds oind_c2 and is omitted
            @inbounds (with2objs = oind_vxl′[ind_c[nc]]==oind_c2) || break
        end
    end

    # Find which of ind_c1 and ind_c2 is the index of the corner occupied by the
    # foreground object.
    #
    # When multiple objects have the same object index (because they are essentially the
    # same object across a periodic boundary), it doesn't matter which object to choose
    # as the foreground (or background) object, we translate points properly across the
    # domain (by ∆fg below) in order to evaluate the surface normal direction on the
    # correct object.  (At least that is the intention, but this needs to be tested
    # after assigning the same object index is actually implemented.)
    if oind_c1 > oind_c2  # ind_c1 is foreground corner
        ind_fg = ind_c1
        oind_fg, oind_bg = oind_c1, oind_c2
    else  # ind_c2 is foreground corner
        @assert oind_c1≠oind_c2
        ind_fg = ind_c2
        oind_fg, oind_bg = oind_c2, oind_c1
    end
    @inbounds shp_fg, shp_bg = oind2shp[oind_fg], oind2shp[oind_bg]
    @inbounds prm_fg, prm_bg = pind2matprm[oind2pind[oind_fg]], pind2matprm[oind2pind[oind_bg]]

    # Find
    # - nout (outward normal of the foreground object), and
    # - rvol (volume fraction of the foreground object inside the voxel).
    ci_cmp = ci_vxl′[1]  # CartesianIndex{K}
    σvxl = t_ind(σcmp, ci_cmp)
    if !with2objs  # two material parameters but more than two objects in voxel
        # In this case, the interface between two materials is not defined by the
        # surface of a single object, so we estimate nout simply from the locations of
        # the corners occupied by the two materials.
        nout, rvol = kottke_input_simple(ind_c, n_diffp)::Tuple{SFloat{K},Float}
    else  # two objects
        # When Nprm_vxl == Nobj_vxl == 2, different material parameters must correspond
        # to different objects.
        x₀ = t_ind(lcmp, ci_cmp)  # SFloat{3}: location of center of smoothing voxel
        @inbounds lvxl′ = (t_ind(lcmp′,ci_vxl′[1]), t_ind(lcmp′,ci_vxl′[2]))

        @inbounds ci_fg = CartesianIndices(ntuple(x->2,Val(K)))[ind_fg]  # subscritpt of corner ind_fg
        ci_1 = CartesianIndex(ntuple(x->1, Val(K)))
        ∆fg = t_ind(∆τcmp′, ci_cmp + ci_fg - ci_1)  # SFloat{3}; nonzero if corner ind_fg is outside periodic boundary

        # See "Overall smoothing algorithm" above.
        nout, rvol = kottke_input_accurate(x₀, σvxl, lvxl′, ∆fg, shp_fg, shp_bg)::Tuple{SFloat{K},Float}
    end  # if !with2objs

    if iszero(nout)
        # Give up Kottke's subpixel smoothing and take simple averaging.
        #
        # Previously this part was using either arithmetic or harmonic averaging,
        # depending on whether the field is on a primal or dual grid line.  The
        # underlying assumption was that the primal voxels were filled uniformly with
        # a single material, which is a common case when the user puts box-shaped
        # objects in the domain.  Then, the material interface coincides with a primal
        # grid plane, and the above choice of averaging scheme has some theoretical
        # support.
        #
        # However, now it is no longer possible to define each field as completely
        # primal or completely dual.  For example, the Ez-field can be on the primal
        # x-normal grid plane and dual y-normal grid plane simultaneously.  Therefore,
        # it is no longer possible to select the averaging scheme based on whether the
        # field is primal or dual, even if we assume each primal voxel is filled
        # uniformly.  Therefore, we simply perform arithmetic averaging regardless of
        # the field type.
        prm_vxl = amean_param(ci_vxl′, oind_cmp′, oind2pind, pind2matprm)::S²ComplexF{Kf,Kf²}
    else
        # Perform Kottke's subpixel smoothing (defined in material.jl).
        prm_vxl = isfield˔shp ? kottke_avg_param(prm_fg, prm_bg, rvol) :  # field is always parallel to shape boundaries
                                kottke_avg_param(prm_fg, prm_bg, nout, rvol)
    end

    return prm_vxl
end

function kottke_input_simple(ind_c::MArray{TP,Int,K,Zᴷ}, n_diffp::Integer) where {TP,K,Zᴷ}
    nout = @SVector zeros(K)  # nout for param_fg
    for n = n_diffp:Zᴷ  # n = Zᴷ corresponds ind_c[Zᴷ] used for param_fg
        @inbounds nout += nout_vxl(Val(K))[ind_c[n]]
    end
    rvol = (Zᴷ+1-n_diffp) / Zᴷ

    return nout, rvol
end

function kottke_input_accurate(x₀::SFloat{K}, σvxl::SBool{K}, lvxl′::Tuple2{SFloat{K}}, ∆fg::SFloat{K},
                               shp_fg::Shape{K,K²}, shp_bg::Shape{K,K²}) where {K,K²}
    r₀, nout = surfpt_nearby(x₀ + ∆fg, shp_fg)
    r₀ -= ∆fg
    nout = σvxl .* nout  # if voxel is across symmetry boundary plane, project nout to plane

    rvol = 0.0  # dummy value
    if !iszero(nout)
        rvol = volfrac(lvxl′, nout, r₀)
    end

    return nout, rvol
end

function amean_param(ci_vxl′::Tuple2{CartesianIndex{K}},
                     oind_cmp′::AbsArr{ObjInd,K},  # object index array; does not change
                     oind2pind::AbsVec{ParamInd},  # input map from oind to pind
                     pind2matprm::AbsVec{S²ComplexF{Kf,Kf²}}  # map from pind to electric (magnetic) material parameters; Kf² = Kf^2
                     ) where {K,Kf,Kf²}
    p = S²ComplexF{Kf,Kf²}(ntuple(x->0, Val(Kf²)))
    for ci = ci_vxl′[1]:ci_vxl′[2]
        @inbounds pc = pind2matprm[oind2pind[oind_cmp′[ci]]]
        p += pc
    end
    return p / 2^K
end

# function hmean_param(ci_vxl′::Tuple2{CartesianIndex{K}},
#                      oind_cmp′::AbsArr{ObjInd,K},  # object index array; does not change
#                      oind2pind::AbsVec{ParamInd},  # input map from oind to pind
#                      pind2matprm::AbsVec{S²ComplexF{Kf,Kf²}}  # map from pind to electric (magnetic) material parameters; Kf² = Kf^2
#                      ) where {K,Kf,Kf²}
#     p = S²ComplexF{Kf,Kf²}(ntuple(x->0, Val(Kf²)))
#     for ci = ci_vxl′[1]:ci_vxl′[2]
#         @inbounds pc = pind2matprm[oind2pind[oind_cmp′[ci]]]
#         p += inv(pc)
#     end
#     return inv(p / 2^K)
# end


text:
using LibPQ
using Compat: Test

import Compat: @__MODULE__
using Compat.Dates
using DataStreams
using Decimals
using Memento
using Missings
using NamedTuples
using OffsetArrays
using TestSetExtensions
using TimeZones

Memento.config("critical")

const TESTSET_TYPE = VERSION < v"0.7-" ? ExtendedTestSet : Test.DefaultTestSet

@testset TESTSET_TYPE "LibPQ" begin

@testset "ConninfoDisplay" begin
    @test parse(LibPQ.ConninfoDisplay, "") == LibPQ.Normal
    @test parse(LibPQ.ConninfoDisplay, "*") == LibPQ.Password
    @test parse(LibPQ.ConninfoDisplay, "D") == LibPQ.Debug
    @test_throws ErrorException parse(LibPQ.ConninfoDisplay, "N")
end

@testset "Version Numbers" begin
    valid_versions = [
        (LibPQ.pqv"11", v"11"),
        (LibPQ.pqv"11.80", v"11.0.80"),
        (LibPQ.pqv"10.1", v"10.0.1"),
        (LibPQ.pqv"9.1.5", v"9.1.5"),
        (LibPQ.pqv"9.2", v"9.2.0"),
        (LibPQ.pqv"8", v"8.0.0"),
    ]

    @testset "Valid Versions" for (pg_version, jl_version) in valid_versions
        @test pg_version == jl_version
    end

    invalid_versions = [
        "10.1.1",
        "10.0.1",
        "10.0.0.1",
        "9.0.0.1",
        "",
    ]

    # can't do cross-version macro testing apparently
    @testset "Invalid Versions" for pg_version_str in invalid_versions
        @test_throws ArgumentError LibPQ._pqv_str(pg_version_str)
    end
end

@testset "Online" begin
    const DATABASE_USER = get(ENV, "LIBPQJL_DATABASE_USER", "postgres")

    @testset "Example SELECT" begin
        conn = LibPQ.Connection("dbname=postgres user=$DATABASE_USER"; throw_error=false)
        @test conn isa LibPQ.Connection
        @test isopen(conn)
        @test status(conn) == LibPQ.libpq_c.CONNECTION_OK
        @test conn.closed == false

        text_display = sprint(show, conn)
        @test contains(text_display, "dbname = postgres")
        @test contains(text_display, "user = $DATABASE_USER")

        result = execute(
            conn,
            "SELECT typname FROM pg_type WHERE oid = 16";
            throw_error=false,
        )
        @test result isa LibPQ.Result
        @test status(result) == LibPQ.libpq_c.PGRES_TUPLES_OK
        @test result.cleared == false
        @test LibPQ.num_columns(result) == 1
        @test LibPQ.num_rows(result) == 1
        @test LibPQ.column_name(result, 1) == "typname"
        @test LibPQ.column_number(result, "typname") == 1

        data = Data.stream!(result, NamedTuple)

        @test data[:typname][1] == "bool"

        clear!(result)
        @test result.cleared == true

        # the same but with parameters
        result = execute(
            conn,
            "SELECT typname FROM pg_type WHERE oid = \$1",
            [16];
            throw_error=false,
        )
        @test result isa LibPQ.Result
        @test status(result) == LibPQ.libpq_c.PGRES_TUPLES_OK
        @test result.cleared == false
        @test LibPQ.num_columns(result) == 1
        @test LibPQ.num_rows(result) == 1
        @test LibPQ.column_name(result, 1) == "typname"

        data = Data.stream!(result, NamedTuple)

        @test data[:typname][1] == "bool"

        clear!(result)
        @test result.cleared == true

        # the same but with fetch
        data = fetch!(NamedTuple, execute(
            conn,
            "SELECT typname FROM pg_type WHERE oid = \$1",
            [16],
        ))

        @test data[:typname][1] == "bool"

        close(conn)
        @test !isopen(conn)
        @test conn.closed == true

        text_display_closed = sprint(show, conn)
        @test contains(text_display_closed, "closed")
    end

    @testset "Example INSERT and DELETE" begin
        conn = LibPQ.Connection("dbname=postgres user=$DATABASE_USER")

        result = execute(conn, """
            CREATE TEMPORARY TABLE libpqjl_test (
                no_nulls    varchar(10) PRIMARY KEY,
                yes_nulls   varchar(10)
            );
        """)
        @test status(result) == LibPQ.libpq_c.PGRES_COMMAND_OK
        clear!(result)

        # get the data from PostgreSQL and let DataStreams construct my NamedTuple
        result = execute(conn, """
            SELECT no_nulls, yes_nulls FROM (
                VALUES ('foo', 'bar'), ('baz', NULL)
            ) AS temp (no_nulls, yes_nulls)
            ORDER BY no_nulls DESC;
            """;
            throw_error=true,
        )
        @test status(result) == LibPQ.libpq_c.PGRES_TUPLES_OK
        @test LibPQ.num_rows(result) == 2
        @test LibPQ.num_columns(result) == 2

        data = Data.stream!(result, NamedTuple)

        @test data[:no_nulls] == ["foo", "baz"]
        @test data[:yes_nulls][1] == "bar"
        @test data[:yes_nulls][2] === missing

        clear!(result)

        stmt = Data.stream!(
            data,
            LibPQ.Statement,
            conn,
            "INSERT INTO libpqjl_test (no_nulls, yes_nulls) VALUES (\$1, \$2);",
        )
        @test_throws ArgumentError num_affected_rows(stmt.description)
        @test num_params(stmt) == 2
        @test num_columns(stmt) == 0  # an insert has no results
        @test LibPQ.column_number(stmt, "no_nulls") == 0
        @test LibPQ.column_names(stmt) == []

        result = execute(
            conn,
            "SELECT no_nulls, yes_nulls FROM libpqjl_test ORDER BY no_nulls DESC;";
            throw_error=true,
        )
        @test status(result) == LibPQ.libpq_c.PGRES_TUPLES_OK
        @test LibPQ.num_rows(result) == 2
        @test LibPQ.num_columns(result) == 2

        table_data = Data.stream!(result, NamedTuple)
        @test table_data[:no_nulls] == data[:no_nulls]
        @test table_data[:yes_nulls][1] == data[:yes_nulls][1]
        @test table_data[:yes_nulls][2] === missing

        clear!(result)

        result = execute(
            conn,
            "DELETE FROM libpqjl_test WHERE no_nulls = 'foo';";
            throw_error=true,
        )
        @test status(result) == LibPQ.libpq_c.PGRES_COMMAND_OK
        @test num_rows(result) == 0
        @test num_affected_rows(result) == 1

        clear!(result)

        result = execute(
            conn,
            "SELECT no_nulls, yes_nulls FROM libpqjl_test ORDER BY no_nulls DESC;";
            throw_error=true,
        )
        table_data_after_delete = Data.stream!(result, NamedTuple)
        @test table_data_after_delete[:no_nulls] == ["baz"]
        @test table_data_after_delete[:yes_nulls][1] === missing

        clear!(result)

        result = execute(
            conn,
            "INSERT INTO libpqjl_test (no_nulls, yes_nulls) VALUES (\$1, \$2), (\$3, \$4);",
            Union{String, Missing}["foo", "bar", "quz", missing],

        )
        @test status(result) == LibPQ.libpq_c.PGRES_COMMAND_OK
        @test num_rows(result) == 0
        @test num_affected_rows(result) == 2

        clear!(result)

        close(conn)
    end

    @testset "LibPQ.Connection" begin
        @testset "do" begin
            local saved_conn

            was_open = LibPQ.Connection("dbname=postgres user=$DATABASE_USER"; throw_error=true) do jl_conn
                saved_conn = jl_conn
                return isopen(jl_conn)
            end

            @test was_open
            @test !isopen(saved_conn)

            @test_throws ErrorException LibPQ.Connection("dbname=123fake"; throw_error=true) do jl_conn
                @test false
            end
        end

        @testset "Version Numbers" begin
            conn = LibPQ.Connection("dbname=postgres user=$DATABASE_USER"; throw_error=true)

            # update this test before PostgreSQL 20.0 ;)
            @test LibPQ.pqv"7" <= LibPQ.server_version(conn) <= LibPQ.pqv"20"
        end

        @testset "Encoding" begin
            conn = LibPQ.Connection("dbname=postgres user=$DATABASE_USER"; throw_error=true)

            @test LibPQ.encoding(conn) == "UTF8"

            LibPQ.set_encoding!(conn, "SQL_ASCII")
            @test LibPQ.encoding(conn) == "SQL_ASCII"
            LibPQ.reset_encoding!(conn)
            @test LibPQ.encoding(conn) == "SQL_ASCII"

            reset!(conn)
            @test LibPQ.encoding(conn) == "SQL_ASCII"
            LibPQ.set_encoding!(conn, "UTF8")
            @test LibPQ.encoding(conn) == "UTF8"
            LibPQ.reset_encoding!(conn)
            @test LibPQ.encoding(conn) == "UTF8"

            conn.encoding = "SQL_ASCII"
            LibPQ.reset_encoding!(conn)
            @test LibPQ.encoding(conn) == "SQL_ASCII"

            @test_throws ErrorException LibPQ.set_encoding!(conn, "NOT A REAL ENCODING")

            close(conn)
        end

        @testset "Bad Connection" begin
            @testset "throw_error=false" begin
                conn = LibPQ.Connection("dbname=123fake"; throw_error=false)
                @test conn isa LibPQ.Connection
                @test status(conn) == LibPQ.libpq_c.CONNECTION_BAD
                @test conn.closed == false

                reset!(conn; throw_error=false)
                @test status(conn) == LibPQ.libpq_c.CONNECTION_BAD
                @test conn.closed == false

                close(conn)
                @test !isopen(conn)
                @test conn.closed == true
                @test_throws ErrorException reset!(conn; throw_error=false)
            end

            @testset "throw_error=true" begin
                @test_throws ErrorException LibPQ.Connection("dbname=123fake"; throw_error=true)

                conn = LibPQ.Connection("dbname=123fake"; throw_error=false)
                @test conn isa LibPQ.Connection
                @test status(conn) == LibPQ.libpq_c.CONNECTION_BAD
                @test conn.closed == false

                @test_throws ErrorException reset!(conn; throw_error=true)
                @test !isopen(conn)
                @test conn.closed == true
                @test_throws ErrorException reset!(conn; throw_error=true)
            end
        end
    end

    @testset "Results" begin
        @testset "Nulls" begin
            conn = LibPQ.Connection("dbname=postgres user=$DATABASE_USER"; throw_error=true)

            result = execute(conn, "SELECT NULL"; throw_error=true)
            @test status(result) == LibPQ.libpq_c.PGRES_TUPLES_OK
            @test LibPQ.num_columns(result) == 1
            @test LibPQ.num_rows(result) == 1

            data = Data.stream!(result, NamedTuple)

            @test data[1][1] === missing

            clear!(result)

            result = execute(conn, """
                SELECT no_nulls, yes_nulls FROM (
                    VALUES ('foo', 'bar'), ('baz', NULL)
                ) AS temp (no_nulls, yes_nulls)
                ORDER BY no_nulls DESC;
                """;
                throw_error=true,
            )
            @test status(result) == LibPQ.libpq_c.PGRES_TUPLES_OK
            @test LibPQ.num_rows(result) == 2
            @test LibPQ.num_columns(result) == 2

            data = Data.stream!(result, NamedTuple)

            @test data[:no_nulls] == ["foo", "baz"]
            @test data[:yes_nulls][1] == "bar"
            @test data[:yes_nulls][2] === missing

            clear!(result)

            # NULL first this time, to check for errors that might come up with lazy
            # initialization of the output data vectors
            result = execute(conn, """
                SELECT no_nulls, yes_nulls FROM (
                    VALUES ('foo', 'bar'), ('baz', NULL)
                ) AS temp (no_nulls, yes_nulls)
                ORDER BY no_nulls ASC;
                """;
                throw_error=true,
            )
            @test result isa LibPQ.Result
            @test status(result) == LibPQ.libpq_c.PGRES_TUPLES_OK
            @test LibPQ.num_rows(result) == 2
            @test LibPQ.num_columns(result) == 2

            data = Data.stream!(result, NamedTuple)

            @test data[:no_nulls] == ["baz", "foo"]
            @test data[:yes_nulls][1] === missing
            @test data[:yes_nulls][2] == "bar"

            clear!(result)
            close(conn)
        end

        @testset "Not Nulls" begin
            conn = LibPQ.Connection("dbname=postgres user=$DATABASE_USER"; throw_error=true)

            result = execute(conn, "SELECT NULL"; not_null=[false], throw_error=true)
            @test status(result) == LibPQ.libpq_c.PGRES_TUPLES_OK
            @test LibPQ.num_columns(result) == 1
            @test LibPQ.num_rows(result) == 1

            data = Data.stream!(result, NamedTuple)

            @test data[1][1] === missing

            clear!(result)

            result = execute(conn, "SELECT NULL"; not_null=true, throw_error=true)
            @test status(result) == LibPQ.libpq_c.PGRES_TUPLES_OK
            @test LibPQ.num_columns(result) == 1
            @test LibPQ.num_rows(result) == 1

            @test_throws ErrorException Data.stream!(result, NamedTuple)

            clear!(result)

            result = execute(conn, "SELECT NULL"; not_null=[true], throw_error=true)
            @test status(result) == LibPQ.libpq_c.PGRES_TUPLES_OK
            @test LibPQ.num_columns(result) == 1
            @test LibPQ.num_rows(result) == 1

            @test_throws ErrorException Data.stream!(result, NamedTuple)

            clear!(result)

            result = execute(conn, """
                SELECT no_nulls, yes_nulls FROM (
                    VALUES ('foo', 'bar'), ('baz', NULL)
                ) AS temp (no_nulls, yes_nulls)
                ORDER BY no_nulls DESC;
                """;
                not_null=[true, false],
                throw_error=true,
            )
            @test status(result) == LibPQ.libpq_c.PGRES_TUPLES_OK
            @test LibPQ.num_rows(result) == 2
            @test LibPQ.num_columns(result) == 2

            data = Data.stream!(result, NamedTuple)

            @test data[:no_nulls] == ["foo", "baz"]
            @test data[:no_nulls] isa Vector{String}
            @test data[:yes_nulls][1] == "bar"
            @test data[:yes_nulls][2] === missing
            @test data[:yes_nulls] isa Vector{Union{String, Missing}}

            clear!(result)

            result = execute(conn, """
                SELECT no_nulls, yes_nulls FROM (
                    VALUES ('foo', 'bar'), ('baz', NULL)
                ) AS temp (no_nulls, yes_nulls)
                ORDER BY no_nulls DESC;
                """;
                not_null=false,
                throw_error=true,
            )
            @test status(result) == LibPQ.libpq_c.PGRES_TUPLES_OK
            @test LibPQ.num_rows(result) == 2
            @test LibPQ.num_columns(result) == 2

            data = Data.stream!(result, NamedTuple)

            @test data[:no_nulls] == ["foo", "baz"]
            @test data[:no_nulls] isa Vector{Union{String, Missing}}
            @test data[:yes_nulls][1] == "bar"
            @test data[:yes_nulls][2] === missing
            @test data[:yes_nulls] isa Vector{Union{String, Missing}}

            clear!(result)
        end

        @testset "Type Conversions" begin
            @testset "Automatic" begin
                conn = LibPQ.Connection("dbname=postgres user=$DATABASE_USER"; throw_error=true)

                result = execute(conn, """
                    SELECT oid, typname, typlen, typbyval, typcategory
                    FROM pg_type
                    WHERE typname IN ('bool', 'int8', 'text')
                    ORDER BY typname;
                    """;
                    throw_error=true,
                )
                @test result isa LibPQ.Result
                @test status(result) == LibPQ.libpq_c.PGRES_TUPLES_OK
                @test LibPQ.num_rows(result) == 3
                @test LibPQ.num_columns(result) == 5
                @test LibPQ.column_types(result) == [LibPQ.Oid, String, Int16, Bool, LibPQ.PQChar]

                data = Data.stream!(result, NamedTuple)

                @test map(eltype, values(data)) == map(T -> Union{T, Missing}, [LibPQ.Oid, String, Int16, Bool, LibPQ.PQChar])
                @test data[:oid] == LibPQ.Oid[LibPQ.PQ_SYSTEM_TYPES[t] for t in (:bool, :int8, :text)]
                @test data[:typname] == ["bool", "int8", "text"]
                @test data[:typlen] == [1, 8, -1]
                @test data[:typbyval] == [true, true, false]
                @test data[:typcategory] == ['B', 'N', 'S']

                clear!(result)
                close(conn)
            end

            @testset "Parsing" begin
                @testset "Default Types" begin
                    conn = LibPQ.Connection("dbname=postgres user=$DATABASE_USER"; throw_error=true)

                    test_data = [
                        ("3", Cint(3)),
                        ("3::int8", Int64(3)),
                        ("3::int4", Int32(3)),
                        ("3::int2", Int16(3)),
                        ("3::float8", Float64(3)),
                        ("3::float4", Float32(3)),
                        ("3::oid", LibPQ.Oid(3)),
                        ("3::numeric", decimal("3")),
                        ("$(BigFloat(pi))::numeric", decimal(BigFloat(pi))),
                        ("$(big"4608230166434464229556241992703")::numeric", parse(Decimal, "4608230166434464229556241992703")),
                        ("E'\\\\xDEADBEEF'::bytea", hex2bytes("DEADBEEF")),
                        ("E'\\\\000'::bytea", UInt8[0o000]),
                        ("E'\\\\047'::bytea", UInt8[0o047]),
                        ("E'\\''::bytea", UInt8[0o047]),
                        ("E'\\\\134'::bytea", UInt8[0o134]),
                        ("E'\\\\\\\\'::bytea", UInt8[0o134]),
                        ("E'\\\\001'::bytea", UInt8[0o001]),
                        ("E'\\\\176'::bytea", UInt8[0o176]),
                        ("'3'::\"char\"", LibPQ.PQChar('3')),
                        ("'t'::bool", true),
                        ("'T'::bool", true),
                        ("'true'::bool", true),
                        ("'TRUE'::bool", true),
                        ("'tRuE'::bool", true),
                        ("'y'::bool", true),
                        ("'YEs'::bool", true),
                        ("'on'::bool", true),
                        ("1::bool", true),
                        ("true", true),
                        ("'f'::bool", false),
                        ("'F'::bool", false),
                        ("'false'::bool", false),
                        ("'FALSE'::bool", false),
                        ("'fAlsE'::bool", false),
                        ("'n'::bool", false),
                        ("'nO'::bool", false),
                        ("'off'::bool", false),
                        ("0::bool", false),
                        ("false", false),
                        ("TIMESTAMP '2004-10-19 10:23:54'", DateTime(2004, 10, 19, 10, 23, 54)),
                        ("TIMESTAMP '2004-10-19 10:23:54.123'", DateTime(2004, 10, 19, 10, 23, 54,123)),
                        ("TIMESTAMP '2004-10-19 10:23:54.1234'", DateTime(2004, 10, 19, 10, 23, 54,123)),
                        ("'infinity'::timestamp", typemax(DateTime)),
                        ("'-infinity'::timestamp", typemin(DateTime)),
                        ("'epoch'::timestamp", DateTime(1970, 1, 1, 0, 0, 0)),
                        # ("TIMESTAMP WITH TIME ZONE '2004-10-19 10:23:54-00'", ZonedDateTime(2004, 10, 19, 10, 23, 54, tz"UTC")),
                        # ("TIMESTAMP WITH TIME ZONE '2004-10-19 10:23:54-02'", ZonedDateTime(2004, 10, 19, 10, 23, 54, tz"UTC-2")),
                        # ("TIMESTAMP WITH TIME ZONE '2004-10-19 10:23:54+10'", ZonedDateTime(2004, 10, 19, 10, 23, 54, tz"UTC+10")),
                        ("'infinity'::timestamptz", ZonedDateTime(typemax(DateTime), tz"UTC")),
                        ("'-infinity'::timestamptz", ZonedDateTime(typemin(DateTime), tz"UTC")),
                        # ("'epoch'::timestamptz", ZonedDateTime(1970, 1, 1, 0, 0, 0, tz"UTC")),
                        ("'{{{1,2,3},{4,5,6}}}'::int2[]", reshape(Int16[1 2 3; 4 5 6], 1, 2, 3)),
                        ("'{}'::int2[]", Int16[]),
                        ("'{{{1,2,3},{4,5,6}}}'::int4[]", reshape(Int32[1 2 3; 4 5 6], 1, 2, 3)),
                        ("'{{{1,2,3},{4,5,6}}}'::int8[]", reshape(Int64[1 2 3; 4 5 6], 1, 2, 3)),
                        ("'{{{1,2,3},{4,5,6}}}'::float4[]", reshape(Float32[1 2 3; 4 5 6], 1, 2, 3)),
                        ("'{{{1,2,3},{4,5,6}}}'::float8[]", reshape(Float64[1 2 3; 4 5 6], 1, 2, 3)),
                        ("'{{{1,2,3},{4,5,6}}}'::oid[]", reshape(LibPQ.Oid[1 2 3; 4 5 6], 1, 2, 3)),
                        ("'{{{1,2,3},{4,5,6}}}'::numeric[]", reshape(Decimal[1 2 3; 4 5 6], 1, 2, 3)),
                        ("'[1:1][-2:-1][3:5]={{{1,2,3},{4,5,6}}}'::int2[]", copy!(OffsetArray(Int16, 1:1, -2:-1, 3:5), [1 2 3; 4 5 6])),
                        ("'[1:1][-2:-1][3:5]={{{1,2,3},{4,5,6}}}'::int4[]", copy!(OffsetArray(Int32, 1:1, -2:-1, 3:5), [1 2 3; 4 5 6])),
                        ("'[1:1][-2:-1][3:5]={{{1,2,3},{4,5,6}}}'::int8[]", copy!(OffsetArray(Int64, 1:1, -2:-1, 3:5), [1 2 3; 4 5 6])),
                        ("'[1:1][-2:-1][3:5]={{{1,2,3},{4,5,6}}}'::float4[]", copy!(OffsetArray(Float32, 1:1, -2:-1, 3:5), [1 2 3; 4 5 6])),
                        ("'[1:1][-2:-1][3:5]={{{1,2,3},{4,5,6}}}'::float8[]", copy!(OffsetArray(Float64, 1:1, -2:-1, 3:5), [1 2 3; 4 5 6])),
                        ("'[1:1][-2:-1][3:5]={{{1,2,3},{4,5,6}}}'::oid[]", copy!(OffsetArray(LibPQ.Oid, 1:1, -2:-1, 3:5), [1 2 3; 4 5 6])),
                        ("'[1:1][-2:-1][3:5]={{{1,2,3},{4,5,6}}}'::numeric[]", copy!(OffsetArray(Decimal, 1:1, -2:-1, 3:5), [1 2 3; 4 5 6])),
                    ]

                    for (test_str, data) in test_data
                        result = execute(conn, "SELECT $test_str;")

                        try
                            @test LibPQ.num_rows(result) == 1
                            @test LibPQ.num_columns(result) == 1
                            @test LibPQ.column_types(result)[1] >: typeof(data)

                            oid = LibPQ.column_oids(result)[1]
                            func = result.column_funcs[1]
                            parsed = func(LibPQ.PQValue{oid}(result, 1, 1))
                            @test parsed == data
                            @test typeof(parsed) == typeof(data)
                        finally
                            clear!(result)
                        end
                    end

                    close(conn)
                end

                @testset "Specified Types" begin
                    conn = LibPQ.Connection("dbname=postgres user=$DATABASE_USER"; throw_error=true)

                    test_data = [
                        ("3", UInt, UInt(3)),
                        ("3::int8", UInt16, UInt16(3)),
                        ("3::int4", Int32, Int32(3)),
                        ("3::int2", UInt8, UInt8(3)),
                        ("3::oid", UInt32, UInt32(3)),
                        ("3::numeric", Float64, 3.0),
                        ("'3'::\"char\"", Char, '3'),
                        ("'foobar'", Symbol, :foobar),
                        ("0::int8", DateTime, DateTime(1970, 1, 1, 0)),
                        ("0::int8", ZonedDateTime, ZonedDateTime(1970, 1, 1, 0, tz"UTC")),
                    ]

                    for (test_str, typ, data) in test_data
                        result = execute(
                            conn,
                            "SELECT $test_str;",
                            column_types=Dict(1 => typ),
                        )

                        try
                            @test LibPQ.num_rows(result) == 1
                            @test LibPQ.num_columns(result) == 1
                            @test LibPQ.column_types(result)[1] >: typeof(data)

                            oid = LibPQ.column_oids(result)[1]
                            func = result.column_funcs[1]
                            parsed = func(LibPQ.PQValue{oid}(result, 1, 1))
                            @test parsed == data
                            @test typeof(parsed) == typeof(data)
                        finally
                            clear!(result)
                        end
                    end

                    close(conn)
                end
            end
        end

        @testset "Query Errors" begin
            @testset "Syntax Errors" begin
                conn = LibPQ.Connection("dbname=postgres user=$DATABASE_USER"; throw_error=true)

                result = execute(conn, "SELORCT NUUL;"; throw_error=false)
                @test status(result) == LibPQ.libpq_c.PGRES_FATAL_ERROR
                # We're expecting zeros per "33.3.2. Retrieving Query Result Information"
                # https://www.postgresql.org/docs/10/static/libpq-exec.html#LIBPQ-EXEC-SELECT-INFO
                @test LibPQ.num_rows(result) == 0
                @test LibPQ.num_columns(result) == 0
                clear!(result)
                @test result.cleared == true

                @test_throws ErrorException execute(conn, "SELORCT NUUL;"; throw_error=true)

                close(conn)
                @test !isopen(conn)
            end

            @testset "Wrong No. Parameters" begin
                conn = LibPQ.Connection("dbname=postgres user=$DATABASE_USER"; throw_error=true)

                result = execute(conn, "SELORCT \$1;", String[]; throw_error=false)
                @test status(result) == LibPQ.libpq_c.PGRES_FATAL_ERROR
                # We're expecting zeros per "33.3.2. Retrieving Query Result Information"
                # https://www.postgresql.org/docs/10/static/libpq-exec.html#LIBPQ-EXEC-SELECT-INFO
                @test LibPQ.num_rows(result) == 0
                @test LibPQ.num_columns(result) == 0
                clear!(result)
                @test result.cleared == true

                @test_throws ErrorException execute(
                    conn,
                    "SELORCT \$1;",
                    String[];
                    throw_error=true,
                )

                close(conn)
                @test !isopen(conn)
            end
        end

        @testset "Interface Errors" begin
            conn = LibPQ.Connection("dbname=postgres user=$DATABASE_USER"; throw_error=true)

            result = execute(
                conn,
                "SELECT typname FROM pg_type WHERE oid = \$1",
                [16],
            )
            clear!(result)
            @test_throws ErrorException fetch!(NamedTuple, result)

            close(conn)
            @test !isopen(conn)
        end
    end

    @testset "Statements" begin
        @testset "No Params, Output" begin
            conn = LibPQ.Connection("dbname=postgres user=$DATABASE_USER"; throw_error=true)

            stmt = prepare(conn, "SELECT oid, typname FROM pg_type")

            @test LibPQ.num_columns(stmt) == 2
            @test LibPQ.num_params(stmt) == 0
            @test LibPQ.column_names(stmt) == ["oid", "typname"]

            result = execute(stmt; throw_error=true)

            @test LibPQ.num_columns(result) == 2
            @test LibPQ.column_names(result) == ["oid", "typname"]
            @test LibPQ.column_types(result) == [LibPQ.Oid, String]
            @test LibPQ.num_rows(result) > 0

            clear!(result)

            close(conn)
        end

        @testset "Params, Output" begin
            conn = LibPQ.Connection("dbname=postgres user=$DATABASE_USER"; throw_error=true)

            stmt = prepare(conn, "SELECT oid, typname FROM pg_type WHERE oid = \$1")

            @test LibPQ.num_columns(stmt) == 2
            @test LibPQ.num_params(stmt) == 1
            @test LibPQ.column_names(stmt) == ["oid", "typname"]

            result = execute(stmt, [16]; throw_error=true)

            @test LibPQ.num_columns(result) == 2
            @test LibPQ.column_names(result) == ["oid", "typname"]
            @test LibPQ.column_types(result) == [LibPQ.Oid, String]
            @test LibPQ.num_rows(result) == 1

            data = fetch!(NamedTuple, result)
            @test data[:oid][1] == 16
            @test data[:typname][1] == "bool"

            clear!(result)

            close(conn)
        end
    end
end

end


text:
function KlobucharIonoModel(arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return KlobucharIonoModel((Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end

function KlobucharIonoModel(arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::TimeScale)
    return KlobucharIonoModel((Vector{jdouble}, Vector{jdouble}, TimeScale), arg0, arg1, arg2)
end

function get_parameters_drivers(obj::KlobucharIonoModel)
    return jcall(obj, "getParametersDrivers", List, ())
end

function path_delay(obj::KlobucharIonoModel, arg0::AbsoluteDate, arg1::GeodeticPoint, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::Vector{jdouble})
    return jcall(obj, "pathDelay", jdouble, (AbsoluteDate, GeodeticPoint, jdouble, jdouble, jdouble, Vector{jdouble}), arg0, arg1, arg2, arg3, arg4, arg5)
end

function path_delay(obj::KlobucharIonoModel, arg0::FieldAbsoluteDate, arg1::FieldGeodeticPoint, arg2::RealFieldElement, arg3::RealFieldElement, arg4::jdouble, arg5::Vector{RealFieldElement})
    return jcall(obj, "pathDelay", RealFieldElement, (FieldAbsoluteDate, FieldGeodeticPoint, RealFieldElement, RealFieldElement, jdouble, Vector{RealFieldElement}), arg0, arg1, arg2, arg3, arg4, arg5)
end

function path_delay(obj::KlobucharIonoModel, arg0::FieldSpacecraftState, arg1::TopocentricFrame, arg2::jdouble, arg3::Vector{RealFieldElement})
    return jcall(obj, "pathDelay", RealFieldElement, (FieldSpacecraftState, TopocentricFrame, jdouble, Vector{RealFieldElement}), arg0, arg1, arg2, arg3)
end

function path_delay(obj::KlobucharIonoModel, arg0::SpacecraftState, arg1::TopocentricFrame, arg2::jdouble, arg3::Vector{jdouble})
    return jcall(obj, "pathDelay", jdouble, (SpacecraftState, TopocentricFrame, jdouble, Vector{jdouble}), arg0, arg1, arg2, arg3)
end


text:
using Test

using SeqLoggers
using LoggingExtras

@testset "Unit Tests SeqLoggers" begin include("unit_tests_loggers.jl") end
@testset "Load from Configs" begin include("unit_load_from_config.jl") end

# only for local tests (fail on server which is not windows)
if !(haskey(ENV, "LOCAL_TEST") && ENV["LOCAL_TEST"] != "TRUE")
@info "Running local tests"
@testset "Advanced file logger" begin include("unit_advanced_file_logger.jl") end
@testset "Integration Tests SeqLoggers" begin include("integration_tests_loggers.jl") end
end


text:
@testset "OIPS" begin
    seed!(42)
    N = 30
    D = 3
    nInd = 20
    kernel = SqExponentialKernel()
    X = ColVecs(rand(D, N) * 10)
    ρ_accept = 0.8
    ρ_remove = 0.9
    alg = OIPS(ρ_accept; ρ_remove=ρ_remove)
    @test repr(alg) ==
        "Online Inducing Point Selection (ρ_accept : $(ρ_accept), ρ_remove : $(ρ_remove), kmax : Inf)"
    Z = inducingpoints(alg, X; kernel=kernel)
    updateZ!(Z, alg, X; kernel=kernel)
    alg = OIPS(nInd)
    Z = inducingpoints(alg, X; kernel=kernel)
    @test length(Z) <= nInd

    @test_throws ArgumentError OIPS(2.0)
    @test_throws ArgumentError OIPS(ρ_accept; η=2.0)
    @test_throws ArgumentError OIPS(ρ_accept; ρ_remove=2.0)

    @test_throws ArgumentError OIPS(-1)
    @test_throws ArgumentError OIPS(nInd, 2.0)

    test_Zalg(OIPS(ρ_accept; ρ_remove=ρ_remove); kernel=kernel)
    test_Zalg(OIPS(nInd); kernel=kernel)
end


text:
#= Raytracer in Julia =#

using Random
using SharedArrays
using Distributed

#==============================================================
Vectors
==============================================================#

struct Vec3
  x::Float64
  y::Float64
  z::Float64
end

Vec3Zero()::Vec3 = Vec3(0,0,0)
vec_sqrt(v::Vec3)::Vec3 = Vec3(sqrt(v.x), sqrt(v.y), sqrt(v.z))
float_eq(a::Float64, b::Float64)::Bool = abs(a - b) < 0.0001
vec_eq(a::Vec3, b::Vec3)::Bool = float_eq(a.x, b.x) && float_eq(a.y, b.y) && float_eq(a.z, b.z)
dot(a::Vec3, b::Vec3)::Float64 = a.x*b.x + a.y*b.y + a.z*b.z
cross(a::Vec3, b::Vec3)::Vec3 = Vec3(a.y*b.z-a.z*b.y, a.z*b.x-a.x*b.z, a.x*b.y-a.y*b.x)
lengthSquared(a::Vec3)::Float64 = a.x*a.x + a.y*a.y + a.z*a.z
length(a::Vec3)::Float64 = sqrt(lengthSquared(a))
add(a::Vec3, b::Vec3)::Vec3 = Vec3(a.x+b.x, a.y+b.y, a.z+b.z)
add(a::Vec3, b::Float64)::Vec3 = Vec3(a.x+b, a.y+b, a.z+b)
subtract(a::Vec3, b::Vec3)::Vec3 = Vec3(a.x-b.x, a.y-b.y, a.z-b.z)
mul(a::Vec3, b::Vec3)::Vec3 = Vec3(a.x*b.x, a.y*b.y, a.z*b.z)
mul(a::Vec3, b::Float64)::Vec3 = Vec3(a.x*b, a.y*b, a.z*b)
mul(a::Float64, b::Vec3)::Vec3 = Vec3(a*b.x, a*b.y, a*b.z)
div(a::Vec3, b::Vec3)::Vec3 = Vec3(a.x/b.x, a.y/b.y, a.z/b.z)
div(a::Vec3, b::Float64)::Vec3 = Vec3(a.x/b, a.y/b, a.z/b)
eq(a::Vec3, b::Vec3)::Bool = a.x==b.x && a.y==b.x && a.z==b.z
negate(a::Vec3)::Vec3 = Vec3(-a.x, -a.y, -a.z)
unit_vector(a::Vec3)::Vec3 = div(a, length(a))

function normalized(a::Vec3)::Vec3
  l_squared = lengthSquared(a)
  if float_eq(l_squared, 0) || float_eq(l_squared, 1)
    return a
  end
  return div(a, sqrt(l_squared))
end

#==============================================================
Rays
==============================================================#
struct Ray
  origin::Vec3
  direction::Vec3
end
pointAtParameter(ray::Ray, t::Float64)::Vec3 = add(ray.origin, mul(ray.direction, t))

#==============================================================
Utils
==============================================================#

function randomInUnitSphere()::Vec3
  p = Vec3(typemax(Float64), typemax(Float64), typemax(Float64))
  while lengthSquared(p) >= 1
    p = subtract(mul(2.0, Vec3(rand(), rand(), rand())), Vec3(1,1,1))
  end; p
end

function reflect(v::Vec3, n::Vec3)::Vec3
  subtract(v, mul(2.0*dot(v, n),n))
end

# This function either returns true and a refracted vector
# or false and Vec3Zero
function refract(v::Vec3, n::Vec3, ni_over_nt::Float64)::Tuple{Bool, Vec3}
  uv = unit_vector(v)
  dt = dot(uv, n)
  discriminant = 1.0 - ni_over_nt * ni_over_nt * (1.0 - dt * dt)
  if discriminant > 0
    refracted = subtract(mul(ni_over_nt, subtract(uv, mul(n, dt))), mul(n, sqrt(discriminant)))
    return (true, refracted)
  end
  return (false, Vec3Zero())
end

# Polynomial approximation
function schlick(cosine::Float64, reflectiveIndex::Float64)::Float64
  r = ((1-reflectiveIndex) / (1+reflectiveIndex))^2
  r+(1-r)*(1-cosine)^5
end

#==============================================================
Image output
==============================================================#

function writePixelArrayToFile(pixels::Array{Vec3, 2})
  open("output.ppm", "w") do f
    width = size(pixels, 1)
    height = size(pixels, 2)
    write(f, "P3\n$width $height\n255\n")
    for y = reverse(1:height), x = 1:width
        ir = convert(Int64, round(255 * pixels[x, y].x))
        ig = convert(Int64, round(255 * pixels[x, y].y))
        ib = convert(Int64, round(255 * pixels[x, y].z))
        write(f, "$ir $ig $ib\n")
    end
  end
end

#==============================================================
Hitable & HitRecord
==============================================================#
struct HitRecord
  t::Float64
  p::Vec3
  normal::Vec3
  material::Any
end
abstract type Hitable end

#==============================================================
Materials
==============================================================#
abstract type Material end
struct Lambertian <: Material
  albedo::Vec3
end

struct Metal <: Material
  albedo::Vec3
  fuzz::Float64
  Metal(albedo::Vec3, fuzz::Float64) = new(albedo, fuzz < 1 ? fuzz : 1)
end

struct Dialetric <: Material
  reflectiveIndex::Float64
end

function scatter(material::Lambertian, ray::Ray, hitRecord::HitRecord)::Tuple{Bool, Vec3, Ray}
  target = add(add(hitRecord.p, hitRecord.normal), randomInUnitSphere())
  attenuation = material.albedo
  scattered = Ray(hitRecord.p, subtract(target, hitRecord.p))
  return (true, attenuation, scattered)
end

function scatter(material::Metal, ray::Ray, hitRecord::HitRecord)::Tuple{Bool, Vec3, Ray}
  reflected = reflect(unit_vector(ray.direction), hitRecord.normal)
  scattered = Ray(hitRecord.p, add(reflected, mul(material.fuzz, randomInUnitSphere())))
  attenuation = material.albedo
  result = dot(scattered.direction, hitRecord.normal) > 0
  return (result, attenuation, scattered)
end

function scatter(material::Dialetric, ray::Ray, hitRecord::HitRecord)::Tuple{Bool, Vec3, Ray}
  ni_over_nt, reflectProb, cosine = 0.0, 0.0, 0.0
  outwardNormal, refracted = Vec3Zero(), Vec3Zero()
  reflected = reflect(ray.direction, hitRecord.normal)
  attenuation = Vec3(1,1,1)
  scattered = Ray(Vec3Zero(), Vec3Zero())

  if dot(ray.direction, hitRecord.normal) > 0
    outwardNormal = negate(hitRecord.normal)
    ni_over_nt = material.reflectiveIndex
    cosine = material.reflectiveIndex * dot(ray.direction, hitRecord.normal) / length(ray.direction)
  else
    outwardNormal = hitRecord.normal
    ni_over_nt = 1.0 / material.reflectiveIndex
    cosine = -dot(ray.direction, hitRecord.normal) / length(ray.direction)
  end
  result, refracted = refract(ray.direction, outwardNormal, ni_over_nt)
  reflectProb = result ? schlick(cosine, material.reflectiveIndex) : 1
  scattered = (rand() < reflectProb) ? Ray(hitRecord.p, reflected) : Ray(hitRecord.p, refracted)
  return (true, attenuation, scattered)
end

#==============================================================
Geometry
==============================================================#
struct Sphere <: Hitable
  center::Vec3
  radius::Float64
  material::Material
end

# Specialized Sphere hitable
# Returns result, hitRecord
function hit(hitable::Sphere, ray::Ray, t_min::Float64, t_max::Float64, hitRecord)::Tuple{Bool, Any}
  oc = subtract(ray.origin, hitable.center)         # Vec3
  a = dot(ray.direction, ray.direction)             # Float64
  b = dot(oc, ray.direction)                        # Float64
  c = dot(oc, oc) - hitable.radius * hitable.radius # Float64
  discriminant = b*b - a*c                          # Float64
  if discriminant > 0
    temp = (-b - sqrt(b*b-a*c)) / a
    if temp < t_max && temp > t_min
      t = temp
      p = pointAtParameter(ray, t)
      normal = div(subtract(p, hitable.center), hitable.radius)
      return (true, HitRecord(t, p, normal, hitable.material))
    end

    temp = (-b + sqrt(b*b-a*c)) / a
    if temp < t_max && temp > t_min
      t = temp
      p = pointAtParameter(ray, t)
      normal = div(subtract(p, hitable.center), hitable.radius)
      return (true, HitRecord(t, p, normal, hitable.material))
    end
  end
  return (false, nothing)
end

struct HitableList <: Hitable
  list::Array{Hitable}
end

# Specialized hitable list hit function
# Returns result, hitRecord
function hit(hitable::HitableList, ray::Ray, t_min::Float64, t_max::Float64, hitRecord)::Tuple{Bool, Any}
    hitAnything = false
    closestSoFar = t_max
    tempHitRecord, hitRecordResult = nothing, nothing
    for element in hitable.list
      result, tempHitRecord = hit(element, ray, t_min, closestSoFar, hitRecord)
      if result
        hitAnything = true
        closestSoFar = tempHitRecord.t
        hitRecordResult = tempHitRecord
      end
    end
    return (hitAnything, hitRecordResult)
end

#==============================================================
Camera
==============================================================#

mutable struct Camera
  origin::Vec3
  lowerLeftCorner::Vec3
  horizontal::Vec3
  vertical::Vec3
  u::Vec3
  v::Vec3
  w::Vec3
  lensRadius::Float64

  Camera(lookFrom::Vec3, lookAt::Vec3, viewUp::Vec3, verticalFOV::Float64, aspect::Float64, aperture::Float64, focusDistance::Float64) = new(
    origin,
    lowerLeftCorner,
    horizontal,
    vertical,
    u, v, w, lensRadius)

  # Convenience constructor that takes the tuple produced by _Camera method
  Camera(t::Tuple{Vec3, Vec3, Vec3, Vec3, Vec3, Vec3, Vec3, Float64}) = new(t[1], t[2], t[3], t[4], t[5], t[6], t[7], t[8])
end

# Helper methods for Camera constructor
halfHeight(vfov::Float64)::Float64 = tan((vfov*pi/180)/2)
halfWidth(aspect::Float64, half_height::Float64)::Float64 = aspect * half_height
function _Camera(lookFrom::Vec3, lookAt::Vec3, vup::Vec3, vfov::Float64, aspect::Float64, aperture::Float64, focusDistance::Float64)
    lensRadius = aperture / 2
    half_height = halfHeight(vfov)
    half_width = halfWidth(aspect, half_height)
    origin = lookFrom
    w = unit_vector(subtract(lookFrom, lookAt))
    u = unit_vector(cross(vup, w))
    v = cross(w, u)
    a = mul((half_width*focusDistance), u)
    b = mul((half_height*focusDistance), v)
    c = mul(focusDistance, w)
    lowerLeftCorner = subtract(subtract(subtract(origin, a), b), c)
    horizontal = mul(2.0,mul(half_width,mul(focusDistance,u)))
    vertical = mul(2.0,mul(half_height,mul(focusDistance, v)))
    return origin, lowerLeftCorner, horizontal, vertical, u, v, w, lensRadius
end

function randomInUnitDisk()::Vec3
  p = Vec3(typemax(Float64), typemax(Float64), typemax(Float64))
  while dot(p, p) >= 1
    p = subtract(mul(2.0, Vec3(rand(), rand(), 0)), Vec3(1,1,0))
  end
  p
end

function getRay(camera::Camera, s::Float64, t::Float64)::Ray
  rd = mul(camera.lensRadius, randomInUnitDisk())
  offset = mul(camera.u, rd.x * rd.y)
  direction = subtract(subtract(add(add(camera.lowerLeftCorner, mul(s, camera.horizontal)), mul(t, camera.vertical)), camera.origin), offset)
  Ray(add(camera.origin, offset), direction)
end

#==============================================================
Rendering
==============================================================#

function colorFromRay(ray::Ray, world::Hitable, depth::Float64)::Vec3
  hitRecord = nothing
  result, hitRecord = hit(world, ray, 0.001, typemax(Float64), hitRecord)
  if result
    scatterResult, attenuation, scattered = scatter(hitRecord.material, ray, hitRecord)
    if depth < 50 && scatterResult
      return mul(attenuation, colorFromRay(scattered, world, depth+1))
    else
      return Vec3Zero()
    end
  else
    unit_direction = unit_vector(ray.direction)
    t = 0.5 * (unit_direction.y + 1)
    return add(mul((1.0 - t), Vec3(1,1,1)), mul(t, Vec3(0.5, 0.7, 1.0)))
  end
end

function scene()::HitableList
  list = Array{Hitable,1}()
  push!(list, Sphere(Vec3(0,-1000,0), 1000, Lambertian(Vec3(0.5, 0.5, 0.5))))

  for a in -11:10
    for b in -11:10
      chooseMat = rand()
      center = Vec3(Float64(a)+0.9*rand(), 0.2, Float64(b)+0.9*rand())
      if length(subtract(center, Vec3(4,0.2,0))) > 0.9
        # diffuse
        if chooseMat < 0.8
          push!(list, Sphere(center, 0.2, Lambertian(Vec3(rand()*rand(), rand()*rand(), rand()*rand()))))
        # metal
        elseif chooseMat < 0.95
            push!(list, Sphere(center, 0.2, Metal(Vec3(0.5*(1+rand()), 0.5*(1+rand()), 0.5*(1+rand())), 0.5*rand())))
        # glass
        else
          push!(list, Sphere(center, 0.2, Dialetric(1.5)))
        end
      end
    end
  end
  push!(list, Sphere(Vec3(0,1,0), 1.0, Dialetric(1.5)))
  push!(list, Sphere(Vec3(-4,1,0), 1.0, Lambertian(Vec3(0.4,0.2,0.1))))
  push!(list, Sphere(Vec3(4,1,0), 1.0, Metal(Vec3(0.7,0.6,0.5),0.0)))
  HitableList(list)
end

function main()

  # Seed random generator
  Random.seed!(0)

  width = 800
  height = 600
  samples = 50
  pixelArray = fill(Vec3Zero(), width, height)
  world = scene()

  lookFrom = Vec3(13,2,3)
  lookAt = Vec3(0,0,0)
  distToFocus = Float64(10)
  aperture = Float64(0.1)
  camera = Camera(_Camera(lookFrom, lookAt, Vec3(0,1,0), Float64(20), Float64(width) / Float64(height), aperture, distToFocus))

  Threads.@threads for j = reverse(1:height)
   for i = 1:width
    color = Vec3Zero()
    for sample = 1:samples
      u = Float64(i + rand()) / Float64(width)
      v = Float64(j + rand()) / Float64(height)
      ray = getRay(camera, u, v)
      color = add(color, colorFromRay(ray, world, 0.0))
    end
    color = div(color, Float64(samples))
    color = vec_sqrt(color) # Gamma correction
    pixelArray[i, j] = color
   end
  end
  writePixelArrayToFile(pixelArray)
end
@time main()

max_stars_repo_name: daniel-beard/DBRayTracerJulia

max_issues_repo_name: daniel-beard/DBRayTracerJulia

max_forks_repo_name: daniel-beard/DBRayTracerJulia


text:
"""
    leave_one_out_split(table; col_user = :userid, col_time = :timestamp)

Leave-one-out split for the input `table`. For each user, whose column is specifed by `col_user`, the items are sorted by `col_time`, and the last one is split into the test set. The others remain in the train set.

# Returns
- `train_table`
- `test_table`
"""
function leave_one_out_split(table; col_user = :userid, col_time = :timestamp)
    last_action = Dict()
    for row in Tables.rows(table)
        user = row[col_user]
        timestamp = row[col_time]
        if !(user in keys(last_action))
            last_action[user] = timestamp
        else
            last_action[user] = max(timestamp, last_action[user])
        end
    end

    train_table =
        table |> TableOperations.filter(row -> row[col_time] < last_action[row[col_user]])
    test_table =
        table |> TableOperations.filter(row -> row[col_time] == last_action[row[col_user]])

    m = Tables.materializer(table)

    return m(train_table), m(test_table)
end

"""
    ratio_split(table, train_ratio = 0.7)

Split the `table` randomly, with the train set ratio specifed by `train_ratio` argument. Current implementaion assumes `table` object that can be converted to `DataFrame`.

# Returns
- `train_table`
- `test_table`
"""
function ratio_split(table, train_ratio = 0.7)
    (train_ratio < 0 || train_ratio > 1) &&
        throw(ArgumentError("train_ratio must be between 0 and 1."))
    # currently it requires table that is converted to DataFrame
    df = DataFrame(table)
    n = nrow(df)
    n_train = round(Int, n * train_ratio)
    idx = randperm(n)
    idx_train = idx[1:n_train]
    idx_test = idx[(n_train+1):end]
    df_train = df[idx_train, :]
    df_test = df[idx_test, :]

    m = Tables.materializer(table)

    return m(df_train), m(df_test)
end


text:
## target

struct CUDACompilerTarget <: CompositeCompilerTarget
    parent::PTXCompilerTarget

    CUDACompilerTarget(cap::VersionNumber) = new(PTXCompilerTarget(cap))
end

Base.parent(target::CUDACompilerTarget) = target.parent

# filter out functions from libdevice and cudadevrt
GPUCompiler.isintrinsic(target::CUDACompilerTarget, fn::String) =
    GPUCompiler.isintrinsic(target.parent, fn) ||
    fn == "__nvvm_reflect" || startswith(fn, "cuda")

GPUCompiler.runtime_module(::CUDACompilerTarget) = CUDAnative

## job

struct CUDACompilerJob <: CompositeCompilerJob
    parent::PTXCompilerJob
end

CUDACompilerJob(target::AbstractCompilerTarget, source::FunctionSpec; kwargs...) =
    CUDACompilerJob(PTXCompilerJob(target, source; kwargs...))

Base.similar(job::CUDACompilerJob, source::FunctionSpec; kwargs...) =
    CUDACompilerJob(similar(job.parent, source; kwargs...))

Base.parent(job::CUDACompilerJob) = job.parent

function GPUCompiler.process_module!(job::CUDACompilerJob, mod::LLVM.Module)
    GPUCompiler.process_module!(job.parent, mod)

    emit_exception_flag!(mod)
end

function GPUCompiler.link_libraries!(job::CUDACompilerJob, mod::LLVM.Module,
                                     undefined_fns::Vector{String})
    GPUCompiler.link_libraries!(job.parent, mod, undefined_fns)

    # FIXME: this is ugly
    link_libdevice!(mod, job.parent.target.parent.cap, undefined_fns)
end


text:
"""Value-at-risk."""
function value_at_risk(x::Vector{Float64}, f::Vector{Float64}, α::Float64)
    i = findfirst(p -> p≥α, cumsum(f))
    if i === nothing
        return x[end]
    else
        return x[i]
    end
end

"""Conditional value-at-risk."""
function conditional_value_at_risk(x::Vector{Float64}, f::Vector{Float64}, α::Float64)
    x_α = value_at_risk(x, f, α)
    if iszero(α)
        return x_α
    else
        tail = x .≤ x_α
        return (sum(x[tail] .* f[tail]) - (sum(f[tail]) - α) * x_α) / α
    end
end

max_stars_repo_name: kafisatz/ConditionalValueAtRisk

max_issues_repo_name: kafisatz/ConditionalValueAtRisk

max_forks_repo_name: kafisatz/ConditionalValueAtRisk


text:
using BasicInterpolators
using Documenter

DocMeta.setdocmeta!(BasicInterpolators, :DocTestSetup, :(using BasicInterpolators); recursive=true)

makedocs(;
    modules=[BasicInterpolators],
    authors="Mark Baum <markmbaum@protonmail.com>",
    repo="https://github.com/markmbaum/BasicInterpolators.jl/blob/{commit}{path}#{line}",
    sitename="BasicInterpolators.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://markmbaum.github.io/BasicInterpolators.jl",
        assets=String[],
    ),
    pages=[
        "Home"=>"index.md",
        "Tutorial"=>"tutorial.md",
        "One Dimension"=>"1d.md",
        "Two Dimensions"=>"2d.md",
        "Scattered, N Dimensions"=>"scattered.md",
        "Boundary Behavior"=>"boundaries.md",
        "Miscellaneous"=>"misc.md"
    ],
)

deploydocs(;
    repo="github.com/markmbaum/BasicInterpolators.jl",
    devbranch="main",
    versions=["stable"=>"v#.#"]
)