text:
using FITSIO

"""
    readfits("abc.fits", [hdu=:], [slices...=:])

Load a FITS file from the given file name or in-memory Byte array (Vector{UInt8}).

By default, this will return a vector of DirectImage for each data extension in the file.
You can also specify `hdu=1` for example, to only load the first image.

```julia
images = readfits("multiextension.fits", :)
one_img = readfits("multiextension.fits", 1)
```

By default, this will load all the data from each extension. You can also specify a subset
of data to load:

```julia
readfits("image.fits", 1, 200:300, 400:900, ...)
```
"""
function readfits(
    fname::Union{String,Vector{UInt8}},
    hdu::Union{Integer,AbstractRange,typeof(:)}=1,
    slices...
)

    return FITS(fname, "r") do fits

        try

            if hdu == (:)
                hdu = 1:length(fits)
            end

            # Handle loading multiple hdus if specified
            map(hdu) do h_no
                data = read(fits[h_no], slices...)
                headers = read_header(fits[h_no])

                img = DirectImage(data)
                for key in eachindex(headers)
                    img[Symbol(key)] = headers[key]
                    img[Symbol(key),/] = get_comment(headers, key)
                end

                if !haskey(headers, "FNAME")
                    img[:FNAME] = basename(fname)
                    img[:FNAME,/] = "file name, as loaded into Julia"
                end
                return img
            end

        catch err
            if !(typeof(err) <: InterruptException)
                @error "Error reading FITS file" exception=err
            end
            rethrow(err)
        end
    end
end
export readfits

"""
Write an array to a FITS file. Works for any array. Does not add headers
"""
function writefits(fname, images::DirectImage...)
    return FITS(fname, "w") do file
        try
            for img in images
                props = propertynames(img)
                headers = FITSHeader(
                    [String(prop) for prop in props],
                    Any[img[prop] for prop in props],
                    String[img[prop,/] for prop in props],
                )
                write(file, collect(arraydata(img)), header=headers)
            end
        catch exp
            println(stderr,exp)
            Base.show_backtrace(stderr)
            rethrow(exp)
        end
    end
end

"""
Write an array to a FITS file. Works for any array. Does not add headers
"""
function writefits(fname, images::AbstractVector{<:DirectImage})
    return FITS(fname, "w") do file
        try
            for img in images
                props = propertynames(img)
                headers = FITSHeader(
                    [String(prop) for prop in props],
                    Any[img[prop] for prop in props],
                    String[img[prop,/] for prop in props],
                )
                write(file, collect(arraydata(img)), header=headers)
            end
        catch exp
            println(stderr,exp)
            Base.show_backtrace(stderr)
            rethrow(exp)
        end
    end
end

"""
Write an array to a FITS file. Works for any array. Does not add headers
"""
function writefits(fname, images::AbstractArray{<:Number}...)
    return FITS(fname, "w") do fits
        try
            for img in images
                write(fits, collect(img))
            end
        catch exp
            println(stderr,exp)
            Base.show_backtrace(stderr)
            rethrow(exp)
        end
    end
end
export writefits

"""
Write an array to a FITS file. Works for any array. Does not add headers
"""
function writefits(fname, images::AbstractVector{<:AbstractArray})
    return FITS(fname, "w") do fits
        try
            for img in images
                write(fits, collect(img))
            end
        catch exp
            println(stderr,exp)
            Base.show_backtrace(stderr)
            rethrow(exp)
        end
    end
end
export writefits

"""
    readheader("info.fits")

Read only the headers from a FITS file.
Returns an empty DirectImage with the headers populated and accessible as properties.

Note: readfits returns data combined with headers. There is no need to run both functions.

```julia
headers = readheader("info.fits")
headers.DATE_OBS == "2020-01-01"
```
"""
function readheader(
    fname::Union{String,Vector{UInt8}},
)

    images = FITS(fname, "r") do fits

        data = Int[]
        headers = read_header(first(fits))

        # Handle loading multiple hdus if specified
        map(hdu) do h_no
            data = read(fits[h_no], slices...)
            headers = read_header(fits[h_no])

            img = DirectImage(data)
            for key in eachindex(headers)
                img[Symbol(key)] = headers[key]
                img[Symbol(key),/] = get_comment(headers, key)
            end

            return hdu
        end
    end

    if typeof(hdu) <: Integer
        return images[1]
    else
        return images
    end
end
export readheader

ext: jl

max_stars_repo_name: sefffal/SpaceImages


text:
using Chain

const Num = Int64

function subject(base::Num, exponent::Num)
    modulus::Num = 20201227
    @assert (modulus-1)^2 >= 0

    re = 1
    base %= modulus
    while exponent > 0
        if exponent % 2 == 1
            re = (re * base) % modulus
        end
        exponent >>= 1
        base = base^2 % modulus
    end
    re
end

@assert subject(7, 8) == 5764801
@assert subject(7, 11) == 17807724
@assert subject(5764801, 11) == 14897079
@assert subject(17807724, 8) == 14897079

function inverse(base::Num, exponent::Num)
    re = 0
    while true
        subject(base, re) == exponent && return re
        re += 1
    end
end

function ex1(keys::Array{Num})
    log = @chain keys begin
        map(i -> inverse(7, i), _)
    end
    subject(keys[1], log[2])
end

problem = [8184785, 5293040]

println("Encryption key: ", ex1(problem))

using Test

@testset "Combo Breaker" begin
    @testset "example 1" begin
    end

    @testset "results" begin
    end
end

max_stars_repo_name: CmdQ/AoC


text:
export InitState

const LocalState = Union{Int,String}

struct InitState
  sts::Vector{LocalState}
  sites::SiteSet
  InitState(sites::SiteSet) = new(fill("",length(sites)),sites)
  InitState(sites::SiteSet,default::LocalState) = new(fill(default,length(sites)),sites)
  InitState(sites::SiteSet,def_vec::Vector{LocalState}) = new(def_vec, sites)
end

getindex(is::InitState,n::Integer) = getindex(is.sts,n)
setindex!(is::InitState,i::LocalState,n::Integer) = setindex!(is.sts,i,n)

sites(is::InitState)::SiteSet = is.sites
site(is::InitState,n::Integer)::Site = site(is.sites,n)
state(is::InitState,n::Integer)::IndexVal = state(site(is,n),is.sts[n])

length(is::InitState) = length(is.sts)

function show(io::IO,
              is::InitState)
  print(io,"InitState")
  (length(is) > 0) && print(io,"\n")
  for n=1:length(is)
    print(io,"  $(site(is,n)) ")
    if typeof(is[n])==Int
      print(io,"$(is[n])\n")
    else
      print(io,"\"$(is[n])\"\n")
    end
  end
end


text:
using PointProcessWavelet
using Test

@testset "PointProcessWavelet.jl" begin
    x0 = rand(2, 100)
    wp = WaveletParams(0.5, 128, 3, 8, 2, 1 / 128)

    vs = v_λ_k_all(x0, wp)
    k0 = K_all(x0, wp, vs)
    w0 = W_all(x0, wp, vs)

    k1 = K_all(x0, x0, wp, vs, vs)
    w1 = W_all(x0, x0, wp, vs, vs)

    @test k0 ≈ k1
    @test w0 ≈ w1
end


text:
### A Pluto.jl notebook ###
# v0.18.4

using Markdown
using InteractiveUtils

# ╔═╡ 2658c9f2-af72-11ec-1bd4-6d9866631eca
using HypertextLiteral

# ╔═╡ e1bb1bd1-a270-42d9-9a15-fd99a821da08
import Deno_jll: deno

# ╔═╡ c5c2c93f-e81b-464b-ae44-318d3d2e800a
bundle_path = joinpath(@__DIR__, "js/index.js")

# ╔═╡ 1be81c04-3e61-46a4-ac52-b3c613d286c2
bundle_code = begin
	error_file = tempname()

	try
		read(pipeline(`$(deno()) bundle $(bundle_path)`, stderr=error_file), String)
	catch e
		if e isa ProcessFailedException
			@htl("""

❌ Bundle failed

$(read(error_file, String))

			""")
		else
			rethrow(e)
		end
	end
end

# ╔═╡ 13e1cd4f-7f83-4d0e-98e7-d9419a9d2f0e
Text(bundle_code)

# ╔═╡ 72856f1a-9784-4463-802e-24e201b2dc3c
function import_local_js(code::AbstractString)

	code_js =
		try
		Main.PlutoRunner.publish_to_js(code)
	catch
		repr(code)
	end

	HypertextLiteral.JavaScript(
		"""
		await (() => {

		window.created_imports = window.created_imports ?? new Map()

		let code = $(code_js)

		if(created_imports.has(code)){
			return created_imports.get(code)
		} else {
			let blob_promise = new Promise((r) => {
        		const reader = new FileReader()
        		reader.onload = async () => r(await import(reader.result))
        		reader.readAsDataURL(
				new Blob([code], {type : "text/javascript"}))
    		})
			created_imports.set(code, blob_promise)
			return blob_promise
		}
		})()
		"""
	)
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Deno_jll = "04572ae6-984a-583e-9378-9577a1c2574d"
HypertextLiteral = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"

[compat]
Deno_jll = "~1.20.4"
HypertextLiteral = "~0.9.3"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.2"
manifest_format = "2.0"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.Deno_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "970da1e64a94f13b51c81691c376a1d5a83a0b3c"
uuid = "04572ae6-984a-583e-9378-9577a1c2574d"
version = "1.20.4+0"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.4.1"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.Preferences]]
deps = ["TOML"]
git-tree-sha1 = "d3538e7f8a790dc8903519090857ef8e1283eecd"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.5"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╠═2658c9f2-af72-11ec-1bd4-6d9866631eca
# ╠═e1bb1bd1-a270-42d9-9a15-fd99a821da08
# ╠═c5c2c93f-e81b-464b-ae44-318d3d2e800a
# ╠═1be81c04-3e61-46a4-ac52-b3c613d286c2
# ╠═13e1cd4f-7f83-4d0e-98e7-d9419a9d2f0e
# ╠═72856f1a-9784-4463-802e-24e201b2dc3c
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
"""
    JopSincRegular(dom, rng, dtdom, dtrng)

Return a sinc interpolation filter applied to a signal in `dom::JotSpace`, regularly resampling
from dtdom to dtrng. It is assumed that the forward operation maps from coarse to dense sampling
along the fast dimension of arrays.
"""
# TODO remove constraint number of dimensions == 2
function JopSincRegular(dom::JetSpace{T,2}, rng::JetSpace{T,2}, dtdom, dtrng, nthreads=12) where {T}
	n1dom,n2dom = size(dom)
	n1rng,n2rng = size(rng)
    @assert n1dom < n1rng   # coarse to dense
    @assert n2dom == n2rng  # same size slow dimension
    filters = WaveFD.interpfilters(T(dtrng), T(dtdom), 0, WaveFD.LangC(), nthreads)
	JopLn(dom = dom, rng = rng, df! = JopSincRegular_df!, df′! = JopSincRegular_df′!, s = (filters=filters, nthreads=nthreads))
end

export JopSincRegular

function JopSincRegular_df!(d::AbstractArray{T,2}, m::AbstractArray{T,2}; filters, nthreads, kwargs...) where {T}
    WaveFD.interpadjoint!(filters, d, m) # from coarse to dense
    d
end

function JopSincRegular_df′!(m::AbstractArray{T,2}, d::AbstractArray{T,2}; filters, nthreads, kwargs...) where {T}
    WaveFD.interpforward!(filters, m, d) # from dense to coarse
    m
end

max_stars_repo_name: ChevronETC/JetPackWaveFD

max_issues_repo_name: ChevronETC/JetPackWaveFD

max_forks_repo_name: ChevronETC/JetPackWaveFD


text:
using JuliaWeBWorK
using SymPy

qs = JuliaWeBWorK.QUESTIONS()

## --------------------------------------------------

meta=(Project="Some sample calculus questions",
     Question="3"
      )

intro = jmt"""

[Details for students](https://webwork.maa.org/wiki/Student_Information)

----

# Sample test

Justify answers and show all work for full credit. No graphing calculators. No internet.

(Wait, that is **joke**).

"""

## --------------------------------------------------

numericq(jmt"""
# Problem 1
Let

```math
f(x) = \frac{4}{5} x^5 + 2x^4 - 20x^3 +  9
```

(a) Find the critical points of  ``f(x)``
""",
         () -> begin
         @syms x
         ex = 4//5 * x^5 + 2x^4  -  20x^3 + 9
         cps = solve(diff(ex,x))
         List(cps)
         end,
         ()) |> qs

## --------------------------------------------------

numericq(jmt"""
(b) Specify all intervals where \(f(x)\) is *increasing*.
""",
         () -> begin
         @syms x
         ex = 4//5 * x^5 + 2x^4  -  20x^3 + 9
         dex = diff(ex,x)
         cps = sort(solve(dex))
         delta = minimum(diff(cps))/2
         l = -Inf
         Ints = []
         for  r in cps
         if dex(r-delta) > 0
         push!(Ints, Interval(l,r))
         end
         l = r
         end
         if dex(l+delta) > 0
         push!(Ints, Interval(l, Inf))
         end
         List(Ints)
         end,
         ()) |> qs

## --------------------------------------------------

numericq(jmt"""

(c) For which critical points does the first derivative change sign from positive to negative?"

""",
         () -> begin
         @syms x
         ex = 4//5 * x^5 + 2x^4  -  20x^3 + 9
         dex = diff(ex,x)
         cps = sort(solve(dex))
         delta = minimum(diff(cps))/2
         xs =  []
         for r in cps
         if (dex(r-delta) > 0) && (dex(r+delta) < 0)
         push!(xs, r)
         end
         end
         List(xs)
         end,
         ()) |> qs

## --------------------------------------------------

r  = randomizer(8:15) |> qs
function q4_helper(a)
    @vars x
    ex = x^4 - 12x^2 + a
    dp = diff(ex,x)
    dpp = diff(ex,(x,2))
    cps =  (sort∘solve)(dp)
    ips =  (sort∘solve)(dpp)
    (ex, cps, ips)
end

numericq(jmt"""

# Problem 2
Let

```math
f(x) = x^4 -  12x^2  + {{:a1}}
```

(a) Find the critical points  of ``f(x)``.
""",
         (a) -> begin
         ex,cps,ips = q4_helper(a)
         List(cps)
         end,
         r) |> qs

numericq(jmt"""

(b) Find  the inflection points of  ``f(x)``.

""",
         (a) -> begin
         ex,cps,ips = q4_helper(a)
         List(ips)
         end,
         r) |> qs

numericq(jmt"""

(c) Find the intervals where  ``f(x)``  is concave  down.

""",
         (a) -> begin
         ex,cps,ips = q4_helper(a)
         x = free_symbols(ex)[1]
         dex = diff(ex,x)
         ddex = diff(ex,x,x)

         δ = minimum(diff(ips))/2
         Ints =  []
         l =  -Inf
         for r in ips
         if  ddex(r-δ) <  0
         push!(Ints, Interval(l,r))
         end
         l = r
         end
         if ddex(l + δ) < 0
         push!(Ints,  Interval(l,Inf))
         end
         List(Ints)
         end,
         r) |> qs

numericq(jmt"""

(d) Which critical points have a negative second derivative?

""",
         (a)  -> begin
         ex,cps,ips = q4_helper(a)
         x = free_symbols(ex)[1]
         xs =  []
         for c  in  cps
         diff(ex,(x,2))(c) < 0 && push!(xs,  c)
         end
         List(xs)
         end,
         r) |> qs

##
## --------------------------------------------------
##
p  =  Page(intro, qs; context="Interval", meta...)  #  Interval here  is needed


text:
using BandedMatrices, BlockBandedMatrices, BlockArrays, LinearAlgebra, ArrayLayouts, Test

@testset "broadcasting" begin
    @testset "general" begin
        N = 10
        A = BlockBandedMatrix{Float64}(undef, 1:N,1:N, (1,1))
        A.data .= randn.()
        n = size(A,1)
        B = Matrix{Float64}(undef, n,n)
        B .= exp.(A)
        @test B == exp.(Matrix(A)) == exp.(A)
        @test exp.(A) isa PseudoBlockMatrix
        @test A .+ 1 isa PseudoBlockMatrix

        A = BandedBlockBandedMatrix{Float64}(undef, 1:N,1:N, (1,1), (1,1))
            A.data .= randn.()
            n = size(A,1)
        B = Matrix{Float64}(undef, n,n)
        B .= exp.(A)
        @test B == exp.(Matrix(A)) == exp.(A)
        @test exp.(A) isa PseudoBlockMatrix
        @test A .+ 1 isa PseudoBlockMatrix
    end

    @testset "lmul!/rmul!" begin
        N = 10
        A = BlockBandedMatrix{Float64}(undef, 1:N,1:N, (1,1))
            A.data .= randn.()
        B = BlockBandedMatrix{Float64}(undef, 1:N,1:N, (2,2))
        B .= (-).(A)
        @test similar(A) isa typeof(A)
        @test similar(A,Float64) isa typeof(A)
        @test -A isa typeof(A)
        @test (-).(A) isa typeof(A)
        @test blockbandwidths(A) == blockbandwidths(-A) == blockbandwidths((-).(A))
        @test B == -A == (-).(A)
        @test A-I isa typeof(A)
        @test I-A isa typeof(A)
        @test blockbandwidths(A) == blockbandwidths(A-I) == blockbandwidths(I-A)

        B .= 2.0.*A

        @test B ==  2A == 2.0.*A
        @test 2A isa typeof(A)
        @test 2.0.*A isa typeof(A)
        @test blockbandwidths(2A) == blockbandwidths(2.0.*A) == blockbandwidths(A)

        A .= 2.0.*A
        @test A == B

        B .= A.*2.0

        @test B ==  A*2 == A.*2.0
        @test A*2 isa typeof(A)
        @test A .* 2.0 isa typeof(A)
        @test blockbandwidths(A*2) == blockbandwidths(A.*2.0) == blockbandwidths(A)
        A .= A.*2.0
        @test A == B

        B .= A ./ 2.0

        @test B == A/2 == A ./ 2.0
        @test A/2 isa typeof(A)
        @test A ./ 2.0 isa typeof(A)
        @test blockbandwidths(A/2) == blockbandwidths(A ./ 2.0) == blockbandwidths(A)

        B .= 2.0 .\ A

        @test B == A/2 == A ./ 2.0
        @test 2\A isa typeof(A)
        @test 2.0 .\ A isa typeof(A)
        @test blockbandwidths(2\A) == blockbandwidths(2.0 .\ A) == blockbandwidths(A)

        A = BandedBlockBandedMatrix{Float64}(undef, 1:N,1:N, (1,1),(1,1))
            A.data .= randn.()
        B = BandedBlockBandedMatrix{Float64}(undef, 1:N,1:N, (2,2),(2,2))
        B .= (-).(A)
        @test similar(A) isa typeof(A)
        @test similar(A,Float64) isa typeof(A)
        @test -A isa typeof(A)
        @test (-).(A) isa typeof(A)
        @test blockbandwidths(A) == blockbandwidths(-A) == blockbandwidths((-).(A))
        @test B == -A == (-).(A)
        @test A-I isa typeof(A)
        @test I-A isa typeof(A)
        @test bandwidths(A) == bandwidths(A-I) == bandwidths(I-A)

        B .= 2.0.*A

        @test B ==  2A == 2.0.*A
        @test 2A isa typeof(A)
        @test 2.0.*A isa typeof(A)
        @test blockbandwidths(2A) == blockbandwidths(2.0.*A) == blockbandwidths(A)
        @test subblockbandwidths(2A) == subblockbandwidths(2.0.*A) == subblockbandwidths(A)

        A .= 2.0.*A
        @test A == B

        B .= A.*2.0

        @test B ==  A*2 == A.*2.0
        @test A*2 isa typeof(A)
        @test A .* 2.0 isa typeof(A)
        @test blockbandwidths(A*2) == blockbandwidths(A.*2.0) == blockbandwidths(A)
        @test subblockbandwidths(A*2) == subblockbandwidths(A.*2.0) == subblockbandwidths(A)
        A .= A.*2.0
        @test A == B

        B .= A ./ 2.0

        @test B == A/2 == A ./ 2.0
        @test A/2 isa typeof(A)
        @test A ./ 2.0 isa typeof(A)
        @test blockbandwidths(A/2) == blockbandwidths(A ./ 2.0) == blockbandwidths(A)
        @test subblockbandwidths(A/2) == subblockbandwidths(A ./ 2.0) == subblockbandwidths(A)

        B .= 2.0 .\ A

        @test B == A/2 == A ./ 2.0
        @test 2\A isa typeof(A)
        @test 2.0 .\ A isa typeof(A)
        @test blockbandwidths(2\A) == blockbandwidths(2.0 .\ A) == blockbandwidths(A)
        @test subblockbandwidths(2\A) == subblockbandwidths(2.0 .\ A) == subblockbandwidths(A)
    end

    @testset "axpy!" begin
        N = 10
        A = BlockBandedMatrix{Float64}(undef, 1:N,1:N, (1,1))
            A.data .= randn.()
        B = BlockBandedMatrix{Float64}(undef, 1:N,1:N, (2,2))
            B.data .= randn.()
        C = BlockBandedMatrix{Float64}(undef, 1:N,1:N, (3,3))
        @time C .= A .+ B
        @test C == A + B == A .+ B

        @test A + B isa typeof(A)
        @test A .+ B isa typeof(A)
        @test blockbandwidths(A+B) == blockbandwidths(A.+B) == (2,2)
        @time B .= A .+ B
        @test B == C

        C .= 2.0 .* A .+ B
        @test C == 2A+B == 2.0.*A .+ B
        @test 2A + B isa typeof(A)
        @test 2.0.*A .+ B isa typeof(A)
        @test blockbandwidths(2A+B) == blockbandwidths(2.0.*A .+ B) == (2,2)
        B .= 2.0 .* A .+ B
        @test B == C

        N = 10
        A = BandedBlockBandedMatrix{Float64}(undef, 1:N,1:N, (1,1), (1,1))
            A.data .= randn.()
        B = BandedBlockBandedMatrix{Float64}(undef, 1:N,1:N, (2,2), (2,2))
            B.data .= randn.()
        C = BandedBlockBandedMatrix{Float64}(undef, 1:N,1:N, (3,3), (3,3))
        @time C .= A .+ B
        @test C == A + B == A .+ B

        @test A + B isa typeof(A)
        @test A .+ B isa typeof(A)
        @test blockbandwidths(A+B) == blockbandwidths(A.+B) == (2,2)
        @test subblockbandwidths(A+B) == subblockbandwidths(A.+B) == (2,2)
        @time B .= A .+ B
        @test B == C

        C .= 2.0 .* A .+ B
        @test C == 2A+B == 2.0.*A .+ B
        @test 2A + B isa typeof(A)
        @test 2.0.*A .+ B isa typeof(A)
        @test blockbandwidths(2A+B) == blockbandwidths(2.0.*A .+ B) == (2,2)
        @test subblockbandwidths(2A+B) == subblockbandwidths(2.0.*A .+ B) == (2,2)
        B .= 2.0 .* A .+ B
        @test B == C
    end

    @testset "Degenerate bands" begin
        A = BandedBlockBandedMatrix{Float64}(undef, Fill(4,4), Fill(4,3), (2,0), (1,1)); A.data .= randn.();
        B = BandedBlockBandedMatrix{Float64}(undef, Fill(4,4), Fill(4,3), (1,-1), (1,1)); B.data .= randn.();
        @test A + B == B + A == Matrix(A) + Matrix(B)
        @test A - B == Matrix(A) - Matrix(B)
        @test B - A == Matrix(B) - Matrix(A)

        A = BandedBlockBandedMatrix{Float64}(undef, Fill(4,4), Fill(4,3), (0,2), (1,1)); A.data .= randn.();
        B = BandedBlockBandedMatrix{Float64}(undef, Fill(4,4), Fill(4,3), (-1,1), (1,1)); B.data .= randn.();
        @test A + B == B + A == Matrix(A) + Matrix(B)
        @test A - B == Matrix(A) - Matrix(B)
        @test B - A == Matrix(B) - Matrix(A)

        A = BlockBandedMatrix{Float64}(undef, Fill(4,4), Fill(4,3), (2,0)); A.data .= randn.();
        B = BandedBlockBandedMatrix{Float64}(undef, Fill(4,4), Fill(4,3), (1,-1), (1,1)); B.data .= randn.();
        @test A + B == B + A == Matrix(A) + Matrix(B)
        @test A - B == Matrix(A) - Matrix(B)
        @test B - A == Matrix(B) - Matrix(A)

        A = BlockBandedMatrix{Float64}(undef, Fill(4,4), Fill(4,3), (2,0)); A.data .= randn.();
        B = BlockBandedMatrix{Float64}(undef, Fill(4,4), Fill(4,3), (1,-1)); B.data .= randn.();
        @test A + B == B + A == Matrix(A) + Matrix(B)
        @test A - B == Matrix(A) - Matrix(B)
        @test B - A == Matrix(B) - Matrix(A)
    end
end


text:
@testset "gaussian" begin
    seed!(42)
    @testset "Likelihood" begin
        σ² = 1e-3
        l = GaussianLikelihood(σ²)
        f = 0.5
        y = 0.2
        @test AGP.noise(l) == σ²
        @test l(y, f) == pdf(Normal(f, sqrt(σ²)), y)
        @test loglikelihood(l, y, f) == logpdf(Normal(f, sqrt(σ²)), y)
        @test repr(l) == "Gaussian likelihood (σ² = $(AGP.noise(l)))"
        @test AGP.AugmentedKL(l, [], (;)) == 0.0
        @test AGP.n_latent(l) == 1
    end
    N, d = 20, 2
    k = SqExponentialKernel() ∘ ScaleTransform(10.0)
    σ = 0.1
    X, f = generate_f(N, d, k)
    y = f + σ * randn(N)
    floattypes = [Float64]
    @testset "GP" begin
        for floattype in floattypes
            model = GP(X, y, k; opt_noise=true, verbose=0)
            @test eltype(model) == floattype
            @test AGP.likelihood(model) isa GaussianLikelihood
            @test AGP.inference(model) isa Analytic
            @test AGP.getf(model) isa AGP.LatentGP
            @test AGP.n_latent(model) == 1
            L = AGP.objective(model, X, y)
            @test_nowarn train!(model, 10)
            @test L < AGP.objective(model, X, y)
            @test testconv(model, "Regression", X, f, y)
            @test all(proba_y(model, X)[2] .> 0)
        end
    end
    @testset "VGP" begin
        @test_throws ErrorException VGP(X, y, k, GaussianLikelihood(), AnalyticVI())
        @test_throws ErrorException VGP(X, y, k, GaussianLikelihood(), QuadratureVI())
        @test_throws ErrorException VGP(X, y, k, GaussianLikelihood(), MCIntegrationVI())
    end
    @testset "SVGP" begin
        Z = inducingpoints(KmeansAlg(M), X)
        for floattype in floattypes
            test_inference_SVGP(
                X,
                y,
                f,
                k,
                GaussianLikelihood(),
                GaussianLikelihood,
                floattype,
                Z,
                1,
                "Regression",
                AnalyticVI(),
                AnalyticSVI(10);
                valid=true,
            )
            @test_throws ErrorException SVGP(k, GaussianLikelihood(), QuadratureVI(), Z)
            @test_throws ErrorException SVGP(k, GaussianLikelihood(), MCIntegrationVI(), Z)
        end
    end
    @testset "MCGP" begin
        @test_throws ErrorException MCGP(X, y, k, GaussianLikelihood(), GibbsSampling())
        @test_throws ErrorException MCGP(X, y, k, GaussianLikelihood(), HMCSampling())
    end
end


text:
struct GenericBPE{T, IT, OT, G, C, Norm}
  sepsym::Union{T, Nothing}
  oldendsym::Union{T, Nothing}
  endsym::Union{T, Nothing}
  input_transform::IT
  output_transform::OT
  merging_rank::Dict{Tuple{T, T}, Int}
  cache::Dict{T, Vector{T}}
  glossaries::G
  codemap::C
  normalizer::Norm
end

function GenericBPE{T}(;sepsym=nothing, oldendsym = nothing, endsym = oldendsym,
                       input_transform = nothing, output_transform = nothing,
                       merging_rank = Dict{Tuple{T, T}, Int}(), cache = Dict{T, Vector{T}}(),
                       glossaries = nothing, codemap = nothing,
                       normalizer = nothing) where T
  return GenericBPE{T}(sepsym, oldendsym, endsym,
                       input_transform, output_transform,
                       merging_rank, cache,
                       glossaries, codemap, normalizer)
end

function GenericBPE(bpe::GenericBPE{T}; kws...) where T
  kws = kws.data
  sepsym = haskey(kws, :sepsym) ? kws.sepsym : bpe.sepsym
  oldendsym = haskey(kws, :oldendsym) ? kws.oldendsym : bpe.oldendsym
  endsym = haskey(kws, :endsym) ? kws.endsym : bpe.endsym
  input_transform = haskey(kws, :input_transform) ? kws.input_transform : bpe.input_transform
  output_transform = haskey(kws, :output_transform) ? kws.output_transform : bpe.output_transform
  merging_rank = haskey(kws, :merging_rank) ? kws.merging_rank : bpe.merging_rank
  cache = haskey(kws, :cache) ? kws.cache : bpe.cache
  glossaries = haskey(kws, :glossaries) ? kws.glossaries : bpe.glossaries
  codemap = haskey(kws, :codemap) ? kws.codemap : bpe.codemap
  normalizer = haskey(kws, :normalizer) ? kws.normalizer : bpe.normalizer
  return GenericBPE{T}(; sepsym, oldendsym, endsym,
                       input_transform, output_transform,
                       merging_rank, cache,
                       glossaries, codemap, normalizer)
end

function GenericBPE{T}(endsym, merging_rank = Dict{Tuple{T, T}, Int}(); input_transform = nothing, codemap = nothing,
                       glossaries = nothing, normalizer = nothing) where T
  return GenericBPE{T}(; oldendsym=endsym, endsym=endsym, input_transform = input_transform,
                       codemap = codemap, glossaries = glossaries, normalizer = normalizer)
end

function GenericBPE{T}(s, oe, e, i, o, m, c, g, cm, n) where T
  return GenericBPE{T, typeof(i), typeof(o), typeof(g), typeof(cm), typeof(n)}(s, oe, e, i, o, m, c, g, cm, n)
end

Base.eltype(::GenericBPE{T}) where T = T
Base.empty!(bpe::GenericBPE) = (empty!(bpe.cache); empty!(bpe.merging_rank))
Base.length(bpe::GenericBPE) = length(bpe.merging_rank)

read_bpefile(bfile::AbstractString; kws...) = open(io->read_bpefile(io; kws...), bfile)
function read_bpefile(io::IO; has_header::Bool=true, merge::Int = -1)
  rank = Dict{Tuple{String,String}, Int}()

  if has_header
    header = readline(io)
    esi = findlast(isequal("#endsym:"), header)
    oldsym = esi === nothing ? nothing : header[last(esi)+1:end]
  else
    oldsym = nothing
  end

  for (i, line) ∈ enumerate(eachline(io))
    merge < 0 || i <= merge || break
    pair = Tuple(intern.(split(line, ' ')))
    rank[pair] = i
  end
  return oldsym, rank
end

"""
  init_merges(bpe::GenericBPE, x)::Vector{Merge}

Generate the initial list of `Merge`s.
"""
function init_merges(bpe::GenericBPE{String}, x)
  ms = Vector{Merge{String}}(undef, length(x))
  for (i, idx) = enumerate(eachindex(x))
    ms[i] = Merge(SubString(x, idx, idx))
  end
  return ms
end

"""
  merges(bpe::GenericBPE, x)::Vector{Merge}

Return the buffer for merging subword.
"""
function merges(bpe::GenericBPE, x)
  buf = init_merges(bpe, x)
  if bpe.oldendsym !== nothing
    @inbounds buf[end] = Merge(buf[end], bpe.oldendsym)
  end
  return buf
end

"""
  bigram(ms::AbstractVector{Merge{T}}, i)

Get the bigram(`Tuple{T, T}`) of `ms` start at index `i`.
"""
bigram(ms::AbstractVector{Merge{T}}, i) where T = @inbounds ms[i], ms[i+1]

"""
  bigram(T, ms::AbstractVector{Merge}, i) where T

Get the bigram(`Tuple{T, T}`) of `ms` start at index `i` and convert `Merge` to type `T`.
`T(::Merge)` need to be defined.
"""
bigram(::Type{T}, ms::AbstractVector{Merge{T}}, i) where T =  T.(bigram(ms, i))

"""
  lowestrank(bpe::GenericBPE{T}, ms::AbstractVector{Merge{T}}) where T

Get index of the bigram in `ms` with lowest rank value and the value. return `(value, index)`.
"""
function lowestrank(bpe::GenericBPE{T}, ms::AbstractVector{Merge{T}}) where T
  m0 = bigram(T, ms, 1)
  lst = get(bpe.merging_rank, m0, typemax(Int))
  li = 1

  for i = 2:(length(ms) - 1)
    m = bigram(T, ms, i)
    r = get(bpe.merging_rank, m, typemax(Int))
    if r < lst
      m0 = m
      lst = r
      li = i
    end
  end
  return lst, li
end

"""
  merge!(ms::AbstractVector{M}, i) where M <: Merge

Merge the bigram to singal `Merge` in `ms`. All bigram equal to the one at index `i` will also be merged.
"""
function merge!(ms::AbstractVector{M}, i) where M <: Merge
  m1, m2 = bigram(ms, i)
  n = m1 * m2
  desidx = i
  len = length(ms) - 1
  i += 2
  @inbounds ms[desidx] = n
  @inbounds while i <= len
    desidx += 1
    m = ms[i]
    if m == m1 && ms[i+1] == m2
      ms[desidx] = n
      i += 2
    else
      ms[desidx] = m
      i += 1
    end
  end

  if i == length(ms)
    desidx += 1
    @inbounds ms[desidx] = ms[i]
  end

  return @inbounds @view(ms[1:desidx])
end

"""
  merge_loop!(bpe::GenericBPE, ms, x)

Iteratively merge subword that is mergeable.
"""
function merge_loop!(bpe::GenericBPE, ms, x)
  for _ in 1:length(ms)
    r, i = lowestrank(bpe, ms)
    r == typemax(Int) && return ms
    ms = merge!(ms, i)
    isone(length(ms)) && return ms
  end
  return ms
end

"""
  bpe_postprocess(bpe::GenericBPE{T}, ms::AbstractVector{Merge{T}}) where T

Convert list of `Merge` to `Vector` of `T` (the encoding result we want).
"""
bpe_postprocess(bpe::GenericBPE{T}, y::AbstractVector{Merge{T}}) where T = map(Base.Fix1(bpe_postprocess, bpe), y)

function bpe_postprocess(bpe::GenericBPE{T}, y::Merge{T}) where T
  transform(s) = T(change_extra(s, bpe.sepsym, bpe.oldendsym, bpe.endsym))
  return intern(transform(y))
end

function flatmap(T, f, xs)
  ys = Vector{T}()
  for x in xs
    y = f(x)
    append!(ys, y)
  end
  return ys
end

"""
  uncache_bpe(bpe::GenericBPE{T}, x::T) where {T}

Byte pair encode `x` without saving the result in `bpe.cache`.
"""
function uncache_bpe(bpe::GenericBPE{T, IT, OT, G}, x) where {T, IT, OT, G}
  !(G <: Nothing) && isgloss(bpe.glossaries, x) && return [bpe_postprocess(bpe, Merge(x, bpe.oldendsym))]
  ms = merges(bpe, x)
  length(ms) < 2 && return [bpe_postprocess(bpe, Merge(x, bpe.oldendsym))]
  y = merge_loop!(bpe, ms, x)
  return bpe_postprocess(bpe, y)
end

uncache_bpe(bpe::GenericBPE, xs::V) where {V <: AbstractVector} = flatmap(eltype(bpe), Base.Fix1(uncache_bpe, bpe), xs)

"""
  bytepairencode(bpe::GenericBPE{T}, x::T) where T

Byte pair encode `x` and save the result in `bpe.cache`.
"""
function bytepairencode(bpe::GenericBPE, x)
  haskey(bpe.cache, x) && return @inbounds bpe.cache[x]
  e = uncache_bpe(bpe, x)
  @inbounds bpe.cache[x] = e
  return e
end

bytepairencode(bpe::GenericBPE, xs::V) where {V <: AbstractVector} = flatmap(eltype(bpe), Base.Fix1(bytepairencode, bpe), xs)

function Base.show(io::IO, bpe::GenericBPE{T}) where T
  print(io, "GenericBPE{$T}(n_merge=$(length(bpe.merging_rank))")
  bpe.endsym === nothing || print(io, ", endsym=$(bpe.endsym)")
  bpe.sepsym === nothing || print(io, ", sepsym=$(bpe.sepsym)")
  bpe.oldendsym === nothing || print(io, ", oldendsym=$(bpe.oldendsym)")
  bpe.input_transform === nothing || print(io, ", input_transform=$(bpe.input_transform)")
  bpe.output_transform === nothing || print(io, ", output_transform=$(bpe.output_transform)")
  bpe.glossaries === nothing || print(io, ", glossaries=$(bpe.glossaries)")
  bpe.codemap === nothing || print(io, ", codemap=$(bpe.codemap)")
  bpe.normalizer === nothing || print(io, ", normalizer=$(bpe.normalizer)")
  print(io, ")")
end

function preprocess(bpe::GenericBPE{T, IT, OT, Nothing, C, Norm}, x) where {T, IT, OT, C, Norm}
  Norm <: Nothing || (x = bpe.normalizer(x))
  IT <: Nothing || (x = bpe.input_transform(x))
  C <: Nothing || (x = bpe.codemap(x))
  return x
end

function preprocess(bpe::GenericBPE{T, IT, OT, G, C, Norm}, x) where {T, IT, OT, G, C, Norm}
  Norm <: Nothing || (x = bpe.normalizer(x))
  x = bpe.glossaries(x)
  g = map(Base.Fix1(isgloss, bpe.glossaries), x)
  IT <: Nothing || begin
    ys = Vector{T}()
    for (xi, gi) in zip(x, g)
      if gi
        push!(ys, xi)
      else
        y = bpe.input_transform(xi)
        append!(ys, y)
      end
    end
    x = ys
  end
  C <: Nothing || (x = bpe.codemap(x))
  return x
end

function (bpe::GenericBPE{T, IT, OT})(x) where {T, IT, OT}
  x = preprocess(bpe, x)
  y  = bytepairencode(bpe, x)
  OT <: Nothing || (y = bpe.output_transform(y))
  return y
end


text:
# Includes explore_complexity()
using DataFrames, Statistics, Dates, CSV, Distributed, Printf, Random
export evolve_complexity, run_complexity_list, complexity_list
export run_explore_complexity, explore_complexity
export average_dataframes, extend_list_by_dups, goal_complexity_frequency_dataframe
export filter_goallist_by_complexity, complexity_freq_scatter_plot
export bin_value, bin_data, bin_counts
export kolmogorov_complexity, run_kolmogorov_complexity

# Distribution of complexities for circuits evolving into a given goal
# Results in data/10_31
function evolve_complexity( g::Goal, p::Parameters, num_circuits::Int64, maxsteps::Int64 )
  n_repeats = 20
  funcs = default_funcs(p.numinputs)
  total_steps = 0
  complexity_list = Float64[]
    c = random_chromosome(p,funcs)
    (c,step,worse,same,better,output,matched_goals,matched_goals_list) = mut_evolve( c, [g], funcs, maxsteps )
    total_steps += step
    count = 0
    while step == maxsteps && count < n_repeats
      (c,step,worse,same,better,output,matched_goals,matched_goals_list) = mut_evolve( c, [g], funcs, maxsteps )
      total_steps += step
      #println("total_steps: ",total_steps)
      count += 1
    end
    if step == maxsteps || count == n_repeats
      println("step == maxsteps in function evolve_to_random_goals()")
    end
    cmplx = complexity5(c)
  #println("(cmplx,total_steps): ",(cmplx,total_steps))
  (cmplx,total_steps)
end

function run_complexity_list( ngoals::Int64, p::Parameters, num_circuits::Int64, maxsteps::Int64; csvfile::String="" )
  goallist = randgoallist( ngoals, p.numinputs, p.numoutputs )
  run_complexity_list(goallist, p, num_circuits,maxsteps, csvfile=csvfile )
end

function run_complexity_list( goallist::GoalList, p::Parameters, num_circuits::Int64, maxsteps::Int64; csvfile::String="" )
  df = DataFrame()
  df.goal = Goal[]
  df.steps = Float64[]
  df.mean = Float64[]
  df.sddev = Float64[]
  df.max = Float64[]
  df.min = Float64[]
  df.q10 = Float64[]
  df.q90 = Float64[]
  df.q95 = Float64[]
  df.q90 = Float64[]
  df.q95 = Float64[]
  df.q99 = Float64[]
  extended_goallist = vcat([ fill(g,num_circuits) for g in goallist ]... )
  #println("length(extended_goallist): ",length(extended_goallist),"  extended_goallist: ",extended_goallist)
  complexity_lists = pmap( g->evolve_complexity( g, p, num_circuits, maxsteps ), extended_goallist )
  #println("complexity_lists: ",complexity_lists)
  ngoals = length(goallist)
  for i = 1:ngoals
    complexity_list = Float64[]
    steps_list = Int64[]
    for j = 1:num_circuits
      #println("(i,j): ",(i,j),"  ind: ",(i-1)*num_circuits+j)
      push!(complexity_list,complexity_lists[(i-1)*num_circuits+j][1])
      push!(steps_list,complexity_lists[(i-1)*num_circuits+j][2])
    end
    #println("steps_list: ",steps_list)
    row_tuple =
      (
        goallist[i],
        mean(steps_list),
        mean(complexity_list),
        std(complexity_list),
        maximum(complexity_list),
        minimum(complexity_list),
        quantile(complexity_list,0.10),
        quantile(complexity_list,0.90),
        quantile(complexity_list,0.95),
        quantile(complexity_list,0.99)
      )
    #println("row_tuple: ",row_tuple)
    push!(df,row_tuple)
  end
  hostname = chomp(open("/etc/hostname") do f read(f,String) end)
  if length(csvfile) > 0
    open( csvfile, "w" ) do f
      println(f,"# date and time: ",Dates.now())
      println(f,"# host: ",hostname," with ",nprocs()-1,"  processes: " )
      println(f,"# funcs: ", Main.CGP.default_funcs(p.numinputs))
      print_parameters(f,p,comment=true)
      println(f,"# maxsteps: ",maxsteps)
      println(f,"# num_circuits: ",num_circuits)
      CSV.write( f, df, append=true, writeheader=true )
    end
  else
    println("# date and time: ",Dates.now())
    println("# host: ",hostname," with ",nprocs()-1,"  processes: " )
    println("# funcs: ", Main.CGP.default_funcs(p.numinputs))
    print_parameters(p,comment=true)
    println("# maxsteps: ",maxsteps)
      println("# num_circuits: ",num_circuits)
  end
  df
end

# Run the next version of run_explore_complexity() nreps time parallelized by pmap().
# This version generates the goallist.
function run_explore_complexity( nreps::Int64, nruns::Int64, p::Parameters, ngoals::Int64, num_circuits::Int64, max_ev_steps::Int64;
      csvfile::String="", insert_gate_prob::Float64=0.0, delete_gate_prob::Float64=0.0 )
  goallist = randgoallist( ngoals, p.numinputs,p.numoutputs)
  run_explore_complexity( nreps, nruns, p, goallist, num_circuits, max_ev_steps, csvfile=csvfile )
end

# Run the next version of run_explore_complexity() nreps time parallelized by pmap().
# This version takes goallist as an argument.
function run_explore_complexity( nreps::Int64, nruns::Int64, p::Parameters, goallist::GoalList, num_circuits::Int64, max_ev_steps::Int64;
      csvfile::String="", insert_gate_prob::Float64=0.0, delete_gate_prob::Float64=0.0 )
  ngoals = length(goallist)
  result_dfs = DataFrame[]
  result_dfs = pmap(x->run_explore_complexity( nruns, p, goallist, num_circuits, max_ev_steps,
        insert_gate_prob=insert_gate_prob, delete_gate_prob=delete_gate_prob), collect(1:nreps))
  #result_dfs = map(x->run_explore_complexity( nruns, p, goallist, num_circuits, max_ev_steps,
  #      insert_gate_prob=insert_gate_prob, delete_gate_prob=delete_gate_prob), collect(1:nreps))
  println("length(result_dfs): ",length(result_dfs))
  df = average_dataframes( result_dfs )
  insertcols!(df, 1, :generation=>collect(1:size(df)[1]))
  cumm_unique_goals = zeros(Float64,size(df)[1])
  cumm_unique_goals[1] = df.unique_goals[1]
  for i = 2:size(df)[1]
    cumm_unique_goals[i] = df.unique_goals[i] + cumm_unique_goals[i-1]
  end
  insertcols!(df, size(df)[2]+1, :cumm_unique_goals=>cumm_unique_goals )
  hostname = chomp(open("/etc/hostname") do f read(f,String) end)
  if length(csvfile) > 0
    open( csvfile, "w" ) do f
      println(f,"# date and time: ",Dates.now())
      println(f,"# host: ",hostname," with ",nprocs()-1,"  processes: " )
      println(f,"# funcs: ", Main.CGP.default_funcs(p.numinputs))
      print_parameters(f,p,comment=true)
      println(f,"# nreps: ",nreps)
      println(f,"# nruns: ",nruns)
      println(f,"# num_circuits: ",num_circuits)
      println(f,"# max_ev_steps: ",max_ev_steps)
      println(f,"# ngoals: ",ngoals)
      CSV.write( f, df, append=true, writeheader=true )
    end
  end
  println("# date and time: ",Dates.now())
  println("# host: ",hostname," with ",nprocs()-1,"  processes: " )
  println("# funcs: ", Main.CGP.default_funcs(p.numinputs))
  print_parameters(p,comment=true)
  println("# nreps: ",nreps)
  println("# nruns: ",nruns)
  println("# num_circuits: ",num_circuits)
  println("# max_ev_steps: ",max_ev_steps)
  println("# ngoals: ",ngoals)
  df
end

# Run explore_complexity() nruns times with the discovered circuit_list from one run becoming the starting circuit_list for the next.
# Also, discoverd goals are removed from goallist on subsequent runs.
function run_explore_complexity( nruns::Int64, p::Parameters, goallist::GoalList, num_circuits::Int64, max_ev_steps::Int64;
      insert_gate_prob::Float64=0.0, delete_gate_prob::Float64=0.0 )
  df = DataFrame()
  done = false
  explore_complexity_tries = 0
  while !done
    explore_complexity_tries += 1
    #println("while !done loop.")
    df = DataFrame()
    df.numgates = Float64[]
    df.levelsback = Float64[]
    df.steps = Int64[]
    df.unique_goals = Int64[]
    df.mean_complexity = Float64[]
    df.max_complexity = Float64[]
    (circuit_list,goals_found,row_tuple) = explore_complexity( p, goallist, num_circuits, max_ev_steps,
        insert_gate_prob=insert_gate_prob, delete_gate_prob=delete_gate_prob )
    push!(df,row_tuple)
    reduced_goallist = goallist
    for i = 2:(nruns-1)
      extended_circuit_list = extend_list_by_dups( circuit_list, num_circuits )
      reduced_goallist = setdiff(reduced_goallist,goals_found)
      #println("run i: ",i,"  length(circuit_list): ",length(circuit_list),"  length(reduced_goallist): ",length(reduced_goallist))
      (circuit_list,goals_found,row_tuple) = explore_complexity( p, reduced_goallist, extended_circuit_list, max_ev_steps,
        insert_gate_prob=insert_gate_prob, delete_gate_prob=delete_gate_prob )
      if row_tuple[1] == 0
        println("i: ",i," break")
        break # Try again
      end
      push!(df,row_tuple)
    end
    extended_circuit_list = extend_list_by_dups( circuit_list, num_circuits )
    reduced_goallist = setdiff(reduced_goallist,goals_found)
    (circuit_list,goals_found,row_tuple) = explore_complexity( p, reduced_goallist, extended_circuit_list, max_ev_steps,
        insert_gate_prob=insert_gate_prob, delete_gate_prob=delete_gate_prob )
    #println("length(redued_goallist): ",length(reduced_goallist),"  row_tuple: ",row_tuple)
    push!(df,row_tuple)
    if row_tuple[3] > 0
      done = true
    end
  end
  println("explore complexity_tries: ",explore_complexity_tries)
  df
end

# This version generates circuit list by calling random_chromosome() to generate each circuit.
function explore_complexity( p::Parameters, goallist::GoalList, num_circuits::Int64, max_ev_steps::Int64;
      insert_gate_prob::Float64=0.0, delete_gate_prob::Float64=0.0 )
  funcs = default_funcs(p.numinputs)
  circuit_list = [ random_chromosome(p,funcs) for _=1:num_circuits ]
  explore_complexity( p, goallist, circuit_list, max_ev_steps,
      insert_gate_prob=insert_gate_prob, delete_gate_prob=delete_gate_prob )
end

# This is the explore_complexity() function that runs one generation starting with circuits in circuit_list.
# max_ev_steps  is the maximum number of steps in the evolution of a circuit to any of the goals in goalllist.
# Returns a triple (circuits_list, goals_found, row_tuple) where:
#   circuits_list is the list of evolved circuits that output one of the goals
#   goals_found is the list of goals that are found by evolution
#   row_tuple is the row of statistics that can be pushed onto a dataframe
function explore_complexity( p::Parameters, goallist::GoalList, circuit_list::Vector{Chromosome}, max_ev_steps::Int64;
      insert_gate_prob::Float64=0.0, delete_gate_prob::Float64=0.0 )
  #println("explore_complexity  insert_gate_prob: ",insert_gate_prob)
  funcs = default_funcs(p.numinputs)
  goals_found = Goal[]
  steps_list = Int64[]
  src_cmplx_list = Float64[]
  dest_cmplx_list = Float64[]
  circuits_list = Chromosome[]
  num_gates_sum = 0
  levelsback_sum = 0
  for c in circuit_list
    (new_c,step,worse,same,better,output,matched_goals,matched_goals_list) =
        mut_evolve( c, goallist, funcs, max_ev_steps, print_steps=false, insert_gate_prob=insert_gate_prob, delete_gate_prob=delete_gate_prob )
    if step < max_ev_steps
      #print("i: ",i,"  steps: ",step)
      #@printf("  goal: 0x%04x\n",matched_goals_list[1][1][1])
      num_gates_sum += new_c.params.numinteriors
      levelsback_sum += new_c.params.numlevelsback
      push!(goals_found, [matched_goals_list[1][1][1]])
      push!(steps_list,step)
      push!(src_cmplx_list, complexity5(c))
      push!(dest_cmplx_list, complexity5(new_c))
      push!(circuits_list, new_c )
    end
  end
  if length(steps_list) > 0
    num_gates_avg = num_gates_sum/length(steps_list)
    levelsback_avg = levelsback_sum/length(steps_list)
    row_tuple = (num_gates_avg, levelsback_avg, length(steps_list), length(unique(goals_found)),mean(dest_cmplx_list), maximum(dest_cmplx_list) )
  else
    row_tuple = (p.numinteriors, p.numlevelsback,0,0,0.0,0.0)
  end
  #println("length(circuits_list): ",length(circuits_list))
  return (circuits_list, goals_found, row_tuple)
end

# Averages a list of dataframes which should all be the same size with the same names.
# All columns should be numerical.  All columns of the result are Float64.
function average_dataframes( dataframes_list::Vector{DataFrame} )
  if length(dataframes_list) == 0
    return nothing
  end
  df = DataFrame()
  i = 1
  for nm in names(dataframes_list[1])
    insertcols!(df,i,nm=>zeros(Float64,size(dataframes_list[1])[1]))
    i += 1
  end
  #println("df: ",df)
  for d in dataframes_list
    @assert size(d) == size(dataframes_list[1])
    @assert names(d) == names(dataframes_list[1])
    #println("  d: ",d)
    for i = 1:size(dataframes_list[1])[2]
      df[:,i] = df[:,i] + d[:,i]
    end
  end
  for i = 1:size(dataframes_list[1])[2]
    df[:,i] = df[:,i]/length(dataframes_list)
  end
  df
end

# Extend vector v to new_length by concatenating v with itself.
function extend_list_by_dups( v::AbstractVector, new_length::Int64 )
  len = length(v)
  if len == 0
    return v
  end
  sum_lengths = len
  result = v
  while sum_lengths < new_length
    if len + sum_lengths <= new_length
      result = vcat(result,v)
      sum_lengths = sum_lengths+len
    else
      result = vcat(result,v[1:(new_length-sum_lengths)])
      sum_lengths = sum_lengths+new_length-len
    end
    #println("sum_lengths: ",sum_lengths,"  result: ",result)
  end
  result
end

# Works on windows if current working directory (by pwd()) is C:\\Users\\oldmtnbiker\\Dropbox\\evotech\\complexity\\data
function goal_complexity_frequency_dataframe( p::Parameters, gl_init::GoalList,
    cmplx_df_csvfile::String="../data/consolidate/geno_pheno_raman_df_all_9_13.csv",
    freq_df_csvfile::String="../data/counts/count_out_4x1_all_ints_10_10.csv" )
  gl_string = map(x->@sprintf("0x%x",x[1]),gl_init) # convert goals from Vector of MyInts to strings
  result_df = DataFrame()
  result_df.goal = gl_init
  cdf = read_dataframe(cmplx_df_csvfile)  # get datafrane with complex and counts
  result_df.complexity =  [cdf[cdf.goal.==gl_string[i],:complex][1] for i = 1:length(gl_string)] # add complexity
  fdf = read_dataframe(freq_df_csvfile)
  result_df.ints8_5 =  [fdf[fdf.goals.==gl_string[i],:ints8_5][1] for i = 1:length(gl_string)] # add complexity
  result_df.ints9_5 =  [fdf[fdf.goals.==gl_string[i],:ints9_5][1] for i = 1:length(gl_string)] # add complexity
  result_df.ints10_5 =  [fdf[fdf.goals.==gl_string[i],:ints10_5][1] for i = 1:length(gl_string)] # add complexity
  result_df.ints11_5 =  [fdf[fdf.goals.==gl_string[i],:ints11_5][1] for i = 1:length(gl_string)] # add complexity
  result_df.ints11_8 =  [fdf[fdf.goals.==gl_string[i],:ints11_8][1] for i = 1:length(gl_string)] # add complexity
  result_df
end

# Uses randgoallist to generate a random goal list of length gl_init_length.
# Reads a dataframe that contains complexities for all goals
# Creates a dataframe with goals and complexities, and sorts by complexities
# Returns the list of goals of length gl_final_length of greater complexity.
# Assumes goals have 1 component
# Works on windows if current working directory (by pwd()) is C:\\Users\\oldmtnbiker\\Dropbox\\evotech\\complexity\\data
function filter_goallist_by_complexity( p::Parameters, gl_init::GoalList, gl_final_length::Int64,
      cmplx_df_csvfile::String="../data/consolidate/geno_pheno_raman_df_all_9_13.csv")
  #gl_init = randgoallist(gl_init_length,p.numinputs, p.numoutputs )
  gl_string = map(x->@sprintf("0x%x",x[1]),gl_init) # convert goals from Vector of MyInts to strings
  result_df = DataFrame()
  result_df.goal = gl_init
  cdf = read_dataframe( cmplx_df_csvfile ) # get datafrane with :complex field
  result_df.complexity =  [cdf[cdf.goal.==gl_string[i],:complex][1] for i = 1:length(gl_string)] # add complexity
  sort!( result_df, order(:complexity, rev=true))
  result_df[1:gl_final_length,:goal]
end

# Possible figfile:  "11_7/complexities_$(length(gl))_random_goals.png"
function complexity_freq_scatter_plot( p::Parameters, gl::GoalList, inverse_power::Int64, count_field::Symbol=:ints11_8;
     denom::Int64=2, figfile::String="")
  inverse_root( x, power, denom ) = x^(1/power)/denom
  ngoals = length(gl)
  gcdf = goal_complexity_frequency_dataframe( p, gl )
  gcdf.crf = map(x->inverse_root(x,inverse_power,denom),gcdf[!,count_field])
  sort!(gcdf, order(:crf, rev=true))
  rp= Random.shuffle!(collect(1:size(gcdf)[1]))   # random permutation of indices
  root_types = ["$(inverse_power)th","square", "cube", "fourth"]
  root_type = inverse_power in [2,3,4] ? root_types[inverse_power] : root_types[1]
  xlabel = "dot diameter is $(root_type) root of goal frequency"
  title = "Complexities of $ngoals random goals"
  plt = Plots.scatter(rp, gcdf.complexity, markersize=gcdf.crf, label="", grid="none",xticks=[],xshowaxis=false,ylabel="complexity",xlabel=xlabel,markercolor=collect(1:ngoals),title=title)
  if length(figfile) > 0
    savefig(figfile)
  end
  plt
end

#  Bin values greater than or equal to v_min with boundaries [v_min, vmin+1.0/bin_fract_denom, v_min+2.0/bin_fract_denom, ... ]
#  Example:  to bin v into intervals  [2.8, 2.9), [2.9,3.0), [3.0,3.1) ... , using bin_value( v, 2.8, 10.0 )
#  values less than v_min will be mapped to 1.
function bin_value( v::Float64, v_min::Float64, bin_fract_denom::Float64 )
  result = Int(floor( bin_fract_denom*(v - v_min )))
  result >= 0.0 ? result+1 : 1    # Increase by 1 since Julia uses 1-based indexing
end

# Example:  let values = [0.3, 0.47, 0.57, 0.4, 0.76, 0.67, 0.04, 0.01, 0.13, 0.88]
# The bin_data( values, 10.0) bins values into intervals [0.0,0.1), [0.1,0.2), [0.2,0.3), . . .
# So bin_data( values, 10.0 ) = [2, 1, 0, 1, 2, 1, 1, 1, 1]  # transposed
function bin_data( values::Vector{Float64}, bin_fract_denom::Float64 )
  v_min = floor(bin_fract_denom*minimum(values))/bin_fract_denom
  ind_max = bin_value(maximum(values),v_min,bin_fract_denom)
  vdict = Dict{Int64,Int64}( i=>0 for i = 1:ind_max )
  for v in values
    vdict[bin_value(v, v_min, bin_fract_denom )] += 1
  end
  [ vdict[key] for key in sort(collect(keys(vdict))) ]
end

#function bin_counts( complexity_count_pairs::Vector{Tuple{Float64,Int64}}, bin_fract_denom::Float64 )
function bin_counts( values::Vector{Float64}, counts::Vector{Float64}, bin_fract_denom::Float64 )
  #values = map(x->x[1],complexity_count_pairs)
  @assert length(values) == length(counts)
  v_min = floor(bin_fract_denom*minimum(values))/bin_fract_denom
  ind_max = bin_value(maximum(values),v_min,bin_fract_denom)
  vdict = Dict{Int64,Int64}( i=>0 for i = 1:ind_max )
  for i in 1:length(values)
    vdict[bin_value(values[i], v_min, bin_fract_denom )] += counts[i]
  end
  [ vdict[key] for key in sort(collect(keys(vdict))) ]
end

# Run kolmogorov_complexity() for all goals in goal list gl (in parallel).
# Write dataframe to csvfile is that is given (keyword argument)
function run_kolmogorov_complexity( p::Parameters, gl::GoalList, max_goal_tries::Int64, max_ev_steps::Int64;
      csvfile::String="" )
  ngoals = length(gl)
  df = DataFrame()
  df.goal = Goal[]
  df.num_gates = Int64[]
  df.num_active_gates = Int64[]
  df.complexity = Float64[]
  df.tries = Int64[]
  df.avg_robustness = Float64[]
  df.avg_evolvability = Float64[]
  df.num_gates_exc = Int64[]
  result_list = pmap( g->kolmogorov_complexity( p, g, max_goal_tries, max_ev_steps ), gl )
  #result_list = map( g->kolmogorov_complexity( p, g, max_goal_tries, max_ev_steps ), gl )
  for r in result_list
    push!(df, r )
  end
  hostname = chomp(open("/etc/hostname") do f read(f,String) end)
  if length(csvfile) > 0
    open( csvfile, "w" ) do f
      println(f,"# date and time: ",Dates.now())
      println(f,"# host: ",hostname," with ",nprocs()-1,"  processes: " )
      println(f,"# funcs: ", Main.CGP.default_funcs(p.numinputs))
      print_parameters(f,p,comment=true)
      println(f,"# max_goal_tries: ",max_goal_tries)
      println(f,"# max_ev_steps: ",max_ev_steps)
      println(f,"# ngoals: ",ngoals)
      CSV.write( f, df, append=true, writeheader=true )
    end
  end
  println("# date and time: ",Dates.now())
  println("# host: ",hostname," with ",nprocs()-1,"  processes: " )
  println("# funcs: ", Main.CGP.default_funcs(p.numinputs))
  print_parameters(p,comment=true)
  println("# max_goal_tries: ",max_goal_tries)
  println("# max_ev_steps: ",max_ev_steps)
  println("# ngoals: ",ngoals)
  df
end

# Try to find the minimum number of gates to evolve a goal which I call the Kolmogorov complexity.
# Start with p.numinteriors and then deccrease the number of gates until the goal is found.
# Do max_goal_tries evolutions with each number of gates.
# Then once min number gates is found, do up to num_tries_multiplier*max_goal_tries further evolutions to get
#    max_goal_tries further approximations of Tononi complexity, robustness, and evolvability.
# If a chromosome with a smaller number of active gates is found in these further evolutions,
#    then num_gates is reset and the further evaluations of complexity, robustness, evolvability are restarted
# max_ev_steps is the maximum number of steps while doing mut_evolve()
function kolmogorov_complexity( p::Parameters, g::Goal, max_goal_tries::Int64, max_ev_steps::Int64 )
  num_tries_multiplier = 3   # used in second outer while loop
  println("goal: ",g)
  num_gates_exceptions = 0
  funcs = default_funcs(p.numinputs)
  complexities_list = Float64[]
  robust_evol_list = Tuple{Float64,Float64}[]
  num_gates = p.numinteriors+1   # decrement on the first iteration
  p_current = p
  found_c = Chromosome(p,[],[],[],0.0,0.0)  # dummy chromosome to establish scope
  goal_found = true
  while num_gates > 1 && goal_found  # terminates when no goal is found for this value of num_gates
    num_gates -= 1
    #println("num_gates: ",num_gates)
    p_current = Parameters( p.numinputs, p.numoutputs, num_gates, p.numlevelsback )
    goal_found = false
    tries = 0
    while !goal_found && tries < max_goal_tries
      tries += 1
      #println("inner while tries: ",tries)
      c = random_chromosome( p_current, funcs )
      (c,step,worse,same,better,output,matched_goals,matched_goals_list) =
          mut_evolve( c, [g], funcs, max_ev_steps, print_steps=false )
      if step < max_ev_steps
        outputs = output_values( c )
        if sort(outputs) != sort(g)
          println("g: ",g,"  outputs: ",outputs)
        end
        try
          @assert sort(outputs) == sort(g)
          catch(e)
          println("g: ",g,"  outputs: ",outputs)
        end
        goal_found = true
        found_c = deepcopy(c)
        num_gates = number_active_gates(found_c)
        p_current = Parameters( p.numinputs, p.numoutputs, num_gates, p.numlevelsback )
        println("circuit found for goal ",g," with num_gates = ",num_gates)
      end
    end
  end
  if num_gates >= 1 && !goal_found
    println("no goal found for goal ",g," for num_gates: ",num_gates)
    num_gates += 1  # now set to the minimum number of gates for successful circuit
  end
  push!(complexities_list, complexity5(found_c))
  push!(robust_evol_list, mutate_all( found_c, funcs,robustness_only=true))
  p_current = Parameters( p.numinputs, p.numoutputs, num_gates, p.numlevelsback )
  tries = 1
  iter = 0  # put a bound on iterations
  while iter < num_tries_multiplier*max_goal_tries && tries < max_goal_tries
    c = random_chromosome( p_current, funcs )
    #println("mut_evolve for goal ",g,"  with numints : ",c.params.numinteriors)
    (c,step,worse,same,better,output,matched_goals,matched_goals_list) =
        mut_evolve( c, [g], funcs, max_ev_steps, print_steps=true )
    if step < max_ev_steps
      outputs = output_values( c )
      #println("outputs: ",outputs,"  goal: ",g)
      @assert sort(outputs) == sort(g)
      if num_gates != number_active_gates(c)
        println("num_gates: ",num_gates,"  number_active_gates(c): ",number_active_gates(c))
        println("B  number gates not equal to number active gates for goal: ",g)
        print_build_chromosome(c)
        found_c = c =  remove_inactive( c )
        num_gates = number_active_gates(found_c)
        p_current = Parameters( p.numinputs, p.numoutputs, num_gates, p.numlevelsback )
        num_gates_exceptions += 1
        complexities_list = Float64[]  # restart complexities_list
        robust_evol_list = Tuple{Float64,Float64}[]  # restart robust_evol_list
        iter = 0
        tries = 1  # do more tries
      end
      push!(complexities_list, complexity5(c))
      push!(robust_evol_list, mutate_all( c, funcs,robustness_only=true))
      tries += 1
    end
    iter += 1
  end
  #println("tries: ",tries)
  avg_complexity = sum(complexities_list)/length(complexities_list)
  avg_robustness = sum( map(x->x[1],robust_evol_list))/length(robust_evol_list)
  avg_evolvability = sum( map(x->x[2],robust_evol_list))/length(robust_evol_list)
  (g,num_gates,number_active_gates(found_c),avg_complexity,tries,avg_robustness,avg_evolvability,num_gates_exceptions)
end


text:
module NetworkInference

using InformationMeasures
using Distributions
using Distributed
using Pkg, DelimitedFiles, SharedArrays, Statistics, StatsBase

export
    # Common types
    Node,
    Edge,
    InferredNetwork,
    # Network inference algorithms
    AbstractNetworkInference,
    CorrelationNetworkInference,
    MINetworkInference,
    CLRNetworkInference,
    PUCNetworkInference,
    PIDCNetworkInference,
    # Functions for inferring networks
    get_nodes,
    write_network_file,
    read_network_file,
    get_adjacency_matrix,
    infer_network

include("common.jl")
include("network_inference.jl")
include("infer_network.jl")
include("empirical_bayes_glue.jl")

# optional exports
if EB_EXISTS
    export # EmpiricalBayes glue functions
        to_index,
        make_priors,
        empirical_bayes
end

end # module


text:
module ChordalGraph

using LightGraphs

export CheckChordal, GreedyOrder!, chordal_cliques!, add_clique!, MinimalChordal!, max_cliques, fill!

include("chordal_extension.jl")

end


text:
using Test
using UnderwaterAcoustics
using UnderwaterAcoustics: amp2db, db2amp
using ForwardDiff

@testset "basic" begin

  @test soundspeed() ≈ 1539.0 atol=0.1
  @test soundspeed(; voidfrac=1e-5) ≈ 1402.1 atol=0.1
  @test soundspeed(; voidfrac=1.0) ≈ 340.0

  @test amp2db(absorption(10000, 1000.0, 35.0, 15.0)) ≈ -1.0 atol=0.5
  @test amp2db(absorption(50000)) ≈ -11.0 atol=0.5
  @test amp2db(absorption(100000)) ≈ -36.0 atol=0.5
  @test amp2db(absorption(100000, 1000.0, 38.5, 14.0, 0.0)) ≈ -40.0 atol=0.5
  @test amp2db(absorption(100000, 1000.0, 38.5, 14.0, 2000.0)) ≈ -30.0 atol=0.5
  @test amp2db(absorption(100000, 1000.0, 38.5, 14.0, 6000.0)) ≈ -16.0 atol=0.5

  @test waterdensity() ≈ 1022.7 atol=0.1

  @test reflectioncoef(0.0, 1200.0/1023.0, 1600.0/1540.0, 0.0) isa Complex
  @test reflectioncoef(0.0, 1200.0/1023.0, 1600.0/1540.0, 0.0) ≈ 0.0986 atol=0.0001
  @test amp2db(abs(reflectioncoef(0.0, 2.5, 2.5, 0.01374))) ≈ -2.8 atol=0.1
  @test amp2db(abs(reflectioncoef(0.0, 2.492, 1.3370, 0.01705))) ≈ -5 atol=0.5
  @test amp2db(abs(reflectioncoef(0.0, 1.195, 1.0179, 0.02158))) ≈ -20 atol=0.5
  @test amp2db(abs(reflectioncoef(1.22, 1.149, 0.9873, 0.00386))) ≈ -32 atol=0.5

  @test amp2db(surfaceloss(15.0, 20000.0, 80°)) ≈ -6.5 atol=0.1
  @test amp2db(surfaceloss(10.0, 20000.0, 80°)) ≈ -3.4 atol=0.1
  @test amp2db(surfaceloss(5.0, 20000.0, 80°)) ≈ -0.5 atol=0.1

  @test doppler(0.0, 50000.0) == 50000.0
  @test doppler(10.0, 50000.0) ≈ 50325 atol=0.5
  @test doppler(-10.0, 50000.0) ≈ 49675 atol=0.5

  @test bubbleresonance(100e-6) ≈ 32500.0 atol=0.1
  @test bubbleresonance(32e-6) ≈ 101562.5 atol=0.1
  @test bubbleresonance(100e-6, 10.0) ≈ 45962.0 atol=0.1

end

@testset "pm-core-basic" begin

  env = UnderwaterEnvironment()

  @test models() isa AbstractArray
  @test models(env) isa AbstractArray

  @test env isa UnderwaterAcoustics.BasicUnderwaterEnvironment
  @test altimetry(env) isa Altimetry
  @test bathymetry(env) isa Bathymetry
  @test ssp(env) isa SoundSpeedProfile
  @test salinity(env) isa Real
  @test seasurface(env) isa ReflectionModel
  @test seabed(env) isa ReflectionModel
  @test noise(env) isa NoiseModel

  @test altitude(altimetry(env), 0.0, 0.0) isa Real
  @test depth(bathymetry(env), 0.0, 0.0) isa Real
  @test soundspeed(ssp(env), 0.0, 0.0, 0.0) isa Real
  @test reflectioncoef(seasurface(env), 1000.0, 0.0) isa Complex
  @test reflectioncoef(seabed(env), 1000.0, 0.0) isa Complex
  sig = record(noise(env), 1.0, 44100.0)
  @test length(sig) == 44100
  @test sum(abs2.(sig)) > 0.0

  @test AcousticReceiver(0.0, 0.0) isa AcousticReceiver
  @test location(AcousticReceiver(100.0, 10.0, -50.0)) == (100, 10.0, -50.0)
  @test location(AcousticReceiver(100.0, -50.0)) == (100, 0.0, -50.0)

  src = AcousticSource(100.0, 10.0, -50.0, 4410.0; sourcelevel=1.0)
  sig = record(src, 1.0, 44100.0)
  @test src isa NarrowbandAcousticSource
  @test location(src) == (100.0, 10.0, -50.0)
  @test location(AcousticSource(100.0, -50.0, 1000.0)) == (100.0, 0.0, -50.0)
  @test nominalfrequency(src) == 4410.0
  @test phasor(src) == complex(1.0, 0.0)
  @test length(sig) == 44100
  @test eltype(sig) <: Complex
  @test sum(abs2.(sig))/44100.0 ≈ 1.0 atol=1e-4
  @test sig[1] != sig[2]
  @test sig[1] ≈ sig[11]
  @test sig[2] ≈ sig[12]
  src = Pinger(10.0, -10.0, 1000.0; sourcelevel=1.0, interval=0.5)
  sig = record(src, 1.0, 44100.0)
  @test src isa AcousticSource
  @test location(src) == (10.0, 0.0, -10.0)
  @test nominalfrequency(src) == 1000.0
  @test phasor(src) == complex(1.0, 0.0)
  @test length(sig) == 44100
  @test eltype(sig) <: Complex
  @test sum(abs2.(sig))/44100.0 ≈ 0.04 atol=1e-4
  @test maximum(abs2.(sig)) ≈ 1.0 atol=1e-4
  src = SampledAcousticSource(10.0, -10.0, ones(1000); fs=1000.0, frequency=100.0)
  @test src isa SampledAcousticSource
  @test location(src) == (10.0, 0.0, -10.0)
  @test nominalfrequency(src) == 100.0
  sig = record(src, 1.0, 1000.0)
  @test length(sig) == 1000
  @test eltype(sig) === Float64
  @test all(sig .== 1.0)
  sig = record(src, 1.0, 1000.0; start=-0.5)
  @test length(sig) == 1000
  @test eltype(sig) === Float64
  @test all(sig[1:500] .== 0.0)
  @test all(sig[501:1000] .== 1.0)
  sig = record(src, 1.0, 1000.0; start=0.5)
  @test length(sig) == 1000
  @test eltype(sig) === Float64
  @test all(sig[1:500] .== 1.0)
  @test all(sig[501:1000] .== 0.0)
  sig = record(src, 1.0, 1000.0; start=-2.0)
  @test length(sig) == 1000
  @test eltype(sig) === Float64
  @test all(sig .== 0.0)
  sig = record(src, 1.0, 1000.0; start=2.0)
  @test length(sig) == 1000
  @test eltype(sig) === Float64
  @test all(sig .== 0.0)
  @test_throws ArgumentError record(src, 1.0, 2000.0)
  src = SampledAcousticSource(10.0, 5.0, -10.0, cis.(2π * 1000 * (0:999) ./ 10000.0); fs=10000.0)
  @test src isa SampledAcousticSource
  @test location(src) == (10.0, 5.0, -10.0)
  @test nominalfrequency(src) == 1000.0
  sig = record(src, 0.1, 10000.0)
  @test length(sig) == 1000
  @test eltype(sig) === ComplexF64

  @test IsoSSP(1500.0) isa SoundSpeedProfile
  @test soundspeed(IsoSSP(1500.0), 100.0, 10.0, -50.0) == 1500.0
  @test MunkSSP() isa SoundSpeedProfile
  @test soundspeed(MunkSSP(), 0.0, 0.0, -2000.0) ≈ 1505.0 atol=1.0
  @test soundspeed(MunkSSP(), 0.0, 0.0, -3000.0) ≈ 1518.0 atol=1.0
  s = SampledSSP(0.0:500.0:1000.0, [1500.0, 1520.0, 1510.0])
  @test s isa SoundSpeedProfile
  @test soundspeed(s, 0.0, 0.0, 0.0) ≈ 1500.0
  @test soundspeed(s, 0.0, 0.0, -250.0) ≈ 1510.0
  @test soundspeed(s, 0.0, 0.0, -500.0) ≈ 1520.0
  @test soundspeed(s, 0.0, 0.0, -750.0) ≈ 1515.0
  @test soundspeed(s, 0.0, 0.0, -1000.0) ≈ 1510.0
  s = SampledSSP(0.0:500.0:1000.0, [1500.0, 1520.0, 1510.0], :smooth)
  @test s isa SoundSpeedProfile
  @test soundspeed(s, 0.0, 0.0, 0.0) ≈ 1500.0
  @test soundspeed(s, 0.0, 0.0, -250.0) > 1510.0
  @test soundspeed(s, 0.0, 0.0, -500.0) ≈ 1520.0
  @test soundspeed(s, 0.0, 0.0, -750.0) > 1515.0
  @test soundspeed(s, 0.0, 0.0, -1000.0) ≈ 1510.0

  @test ConstantDepth(20.0) isa Bathymetry
  @test depth(ConstantDepth(20.0), 0.0, 0.0) == 20.0
  @test maxdepth(ConstantDepth(20.0)) == 20.0
  b = SampledDepth(0.0:500.0:1000.0, [20.0, 25.0, 15.0])
  @test b isa Bathymetry
  @test depth(b, 0.0, 0.0) ≈ 20.0
  @test depth(b, 250.0, 0.0) ≈ 22.5
  @test depth(b, 500.0, 0.0) ≈ 25.0
  @test depth(b, 750.0, 0.0) ≈ 20.0
  @test depth(b, 1000.0, 0.0) ≈ 15.0
  @test maxdepth(b) ≈ 25.0
  b = SampledDepth(0.0:500.0:1000.0, [20.0, 25.0, 15.0], :smooth)
  @test b isa Bathymetry
  @test depth(b, 0.0, 0.0) ≈ 20.0
  @test depth(b, 250.0, 0.0) > 22.5
  @test depth(b, 500.0, 0.0) ≈ 25.0
  @test depth(b, 750.0, 0.0) > 20.0
  @test depth(b, 1000.0, 0.0) ≈ 15.0
  @test maxdepth(b) ≈ 25.0

  @test FlatSurface() isa Altimetry
  @test altitude(FlatSurface(), 0.0, 0.0) == 0.0
  a = SampledAltitude(0.0:500.0:1000.0, [0.0, 1.0, -1.0])
  @test a isa Altimetry
  @test altitude(a, 0.0, 0.0) ≈ 0.0
  @test altitude(a, 250.0, 0.0) ≈ 0.5
  @test altitude(a, 500.0, 0.0) ≈ 1.0
  @test altitude(a, 750.0, 0.0) ≈ 0.0
  @test altitude(a, 1000.0, 0.0) ≈ -1.0
  a = SampledAltitude(0.0:500.0:1000.0, [0.0, 1.0, -1.0], :smooth)
  @test a isa Altimetry
  @test altitude(a, 0.0, 0.0) ≈ 0.0 atol=1e-6
  @test altitude(a, 250.0, 0.0) > 0.5
  @test altitude(a, 500.0, 0.0) ≈ 1.0 atol=1e-6
  @test altitude(a, 1000.0, 0.0) ≈ -1.0 atol=1e-6

  @test ReflectionCoef(0.5 + 0.3im) isa ReflectionModel
  @test reflectioncoef(ReflectionCoef(0.5 + 0.3im), 1000.0, 0.0) == 0.5 + 0.3im
  @test RayleighReflectionCoef(1.0, 1.0) isa ReflectionModel
  @test RayleighReflectionCoef(1.0, 1.0, 0.0) isa ReflectionModel
  @test reflectioncoef(RayleighReflectionCoef(1.0, 1.0, 0.0), 1000.0, 0.0) ≈ 0.0
  @test reflectioncoef(RayleighReflectionCoef(0.0, 1.0, 0.0), 1000.0, 0.0) ≈ -1.0
  @test Rock isa RayleighReflectionCoef
  @test Pebbles isa RayleighReflectionCoef
  @test SandyGravel isa RayleighReflectionCoef
  @test CoarseSand isa RayleighReflectionCoef
  @test MediumSand isa RayleighReflectionCoef
  @test FineSand isa RayleighReflectionCoef
  @test VeryFineSand isa RayleighReflectionCoef
  @test ClayeySand isa RayleighReflectionCoef
  @test CoarseSilt isa RayleighReflectionCoef
  @test SandySilt isa RayleighReflectionCoef
  @test Silt isa RayleighReflectionCoef
  @test FineSilt isa RayleighReflectionCoef
  @test SandyClay isa RayleighReflectionCoef
  @test SiltyClay isa RayleighReflectionCoef
  @test Clay isa RayleighReflectionCoef
  @test Vacuum isa ReflectionModel
  @test reflectioncoef(Vacuum, 1000.0, 0.0) ≈ -1.0
  @test 0.0 < abs(reflectioncoef(Rock, 1000.0, 0.0)) < 1.0
  @test abs(reflectioncoef(Silt, 1000.0, 0.0)) < abs(reflectioncoef(SandyGravel, 1000.0, 0.0))
  @test abs(reflectioncoef(CoarseSilt, 1000.0, 0.0)) < abs(reflectioncoef(Rock, 1000.0, 0.0))
  @test SurfaceLoss(5.0) isa ReflectionModel
  @test SeaState0 isa SurfaceLoss
  @test SeaState1 isa SurfaceLoss
  @test SeaState2 isa SurfaceLoss
  @test SeaState3 isa SurfaceLoss
  @test SeaState4 isa SurfaceLoss
  @test SeaState5 isa SurfaceLoss
  @test SeaState6 isa SurfaceLoss
  @test SeaState7 isa SurfaceLoss
  @test SeaState8 isa SurfaceLoss
  @test SeaState9 isa SurfaceLoss
  @test 0.0 < abs(reflectioncoef(SeaState2, 1000.0, 0.0)) < 1.0
  @test abs(reflectioncoef(SeaState2, 1000.0, 0.0)) > abs(reflectioncoef(SeaState5, 1000.0, 0.0))
  @test abs(reflectioncoef(SeaState5, 1000.0, 0.0)) > abs(reflectioncoef(SeaState5, 2000.0, 0.0))

  rx = AcousticReceiverGrid2D(100.0, 2.0, 100, -100.0, 5.0, 10)
  @test rx isa AcousticReceiverGrid2D
  @test rx isa AbstractMatrix
  @test size(rx) == (100, 10)
  @test rx[1,1] == AcousticReceiver(100.0, -100.0)
  @test rx[end,end] == AcousticReceiver(298.0, -55.0)
  rx = AcousticReceiverGrid3D(100.0, 2.0, 100, 0.0, 1.0, 100, -100.0, 5.0, 10)
  @test rx isa AcousticReceiverGrid3D
  @test rx isa AbstractArray
  @test size(rx) == (100, 100, 10)
  @test rx[1,1,1] == AcousticReceiver(100.0, 0.0, -100.0)
  @test rx[end,end,end] == AcousticReceiver(298.0, 99.0, -55.0)

end

@testset "pm-pekeris" begin

  @test PekerisRayModel in models()

  env = UnderwaterEnvironment()
  pm = PekerisRayModel(env, 7)
  @test pm isa PekerisRayModel

  arr = arrivals(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiver(100.0, -10.0))
  @test arr isa AbstractArray{<:UnderwaterAcoustics.RayArrival}
  @test length(arr) == 7
  @test arr[1].time ≈ 0.0650 atol=0.0001
  @test arr[2].time ≈ 0.0657 atol=0.0001
  @test arr[3].time ≈ 0.0670 atol=0.0001
  @test all([arr[j].time > arr[j-1].time for j ∈ 2:7])
  @test abs(arr[1].phasor) ≈ 0.01 atol=0.001
  @test real(arr[2].phasor) < 0.0
  @test imag(arr[2].phasor) ≈ 0.0
  @test all([abs(arr[j].phasor) < abs(arr[j-1].phasor) for j ∈ 2:7])
  @test [(arr[j].surface, arr[j].bottom) for j ∈ 1:7] == [(0,0), (1,0), (0,1), (1,1), (1,1), (2,1), (1,2)]
  @test all([abs(arr[j].arrivalangle) == abs(arr[j].launchangle) for j ∈ 1:7])

  r = eigenrays(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiver(100.0, -10.0))
  @test r isa AbstractArray{<:UnderwaterAcoustics.RayArrival}
  @test length(r) == 7
  @test all([abs(r[j].arrivalangle) == abs(r[j].launchangle) for j ∈ 1:7])
  @test all([r[j].raypath[1] == (0.0, 0.0, -5.0) for j ∈ 1:7])
  @test all([r[j].raypath[end] == (100.0, 0.0, -10.0) for j ∈ 1:7])
  @test all([length(r[j].raypath) == r[j].surface + r[j].bottom + 2 for j ∈ 1:7])

  r = rays(pm, AcousticSource(0.0, -5.0, 1000.0), -60°:15°:60°, 100.0)
  @test r isa AbstractArray{<:UnderwaterAcoustics.RayArrival}
  @test length(r) == 9
  @test all([r[j].launchangle for j ∈ 1:9] .≈ -60°:15°:60°)
  @test all([abs(r[j].arrivalangle) == abs(r[j].launchangle) for j ∈ 1:9])
  @test all([r[j].raypath[1] == (0.0, 0.0, -5.0) for j ∈ 1:9])
  @test all([r[j].raypath[end][1] ≥ 100.0 for j ∈ 1:9])
  @test all([length(r[j].raypath) == r[j].surface + r[j].bottom + 2 for j ∈ 1:9])

  ir1 = impulseresponse(arr, 10000.0; reltime=true, approx=true)
  ir2 = impulseresponse(arr, 10000.0; reltime=false, approx=true)
  @test length(ir2) ≈ length(ir1) + round(Int, 10000.0 * arr[1].time) atol=1
  @test length(ir2) == round(Int, 10000.0 * arr[end].time) + 1
  @test sum(ir1 .!= 0.0) == 7
  @test sum(ir2 .!= 0.0) == 7
  ndx = findall(abs.(ir1) .> 0)
  @test (ndx .- 1) ./ 10000.0 ≈ [arr[j].time - arr[1].time for j ∈ 1:7] atol=1e-4
  ndx = findall(abs.(ir2) .> 0)
  @test (ndx .- 1) ./ 10000.0 ≈ [arr[j].time for j ∈ 1:7] atol=1e-4

  ir1a = impulseresponse(arr, 10000.0; reltime=true)
  ir2a = impulseresponse(arr, 10000.0; reltime=false)
  @test length(ir2a) ≈ length(ir1a) + round(Int, 10000.0 * arr[1].time) atol=1
  @test length(ir2a) ≥ length(ir2)
  @test sum(abs2.(ir1a))/sum(abs2.(ir1)) ≈ 1.0 atol=0.05
  @test sum(abs2.(ir2a))/sum(abs2.(ir2)) ≈ 1.0 atol=0.05

  @test length(impulseresponse(arr, 10000.0, 256; reltime=true, approx=true)) == 256
  @test length(impulseresponse(arr, 10000.0, 64; reltime=true, approx=true)) == 64
  @test length(impulseresponse(arr, 10000.0, 256; reltime=true, approx=false)) == 256
  @test length(impulseresponse(arr, 10000.0, 64; reltime=true, approx=false)) == 64
  @test length(impulseresponse(arr, 10000.0, 1024; reltime=false, approx=true)) == 1024
  @test length(impulseresponse(arr, 10000.0, 700; reltime=false, approx=true)) == 700
  @test length(impulseresponse(arr, 10000.0, 1024; reltime=false, approx=false)) == 1024
  @test length(impulseresponse(arr, 10000.0, 700; reltime=false, approx=false)) == 700

  env = UnderwaterEnvironment(ssp=IsoSSP(1500.0))
  pm = PekerisRayModel(env, 2)
  d = (√1209.0)/4.0
  x = transfercoef(pm, AcousticSource(0.0, -d, 1000.0), AcousticReceiver(100.0, -d))
  @test x isa Complex
  @test abs(x) ≈ 0.0 atol=0.0002
  x′ = transfercoef(pm, AcousticSource(0.0, -d, 1000.0), AcousticReceiver(100.0, -d); mode=:incoherent)
  @test x′ isa Complex
  @test imag(x′) == 0.0
  @test abs(x′) > 1/100.0
  d = (√2409.0)/8.0
  x = transfercoef(pm, AcousticSource(0.0, -d, 1000.0), AcousticReceiver(100.0, -d))
  @test abs(x) > abs(x′)
  y = transmissionloss(pm, AcousticSource(0.0, -d, 1000.0), AcousticReceiver(100.0, -d))
  @test -10 * log10(abs2(x)) ≈ y atol=0.1
  x = transfercoef(pm, AcousticSource(0.0, -d, 1000.0), AcousticReceiver(100.0, -d); mode=:incoherent)
  @test abs(x) ≈ abs(x′) atol=0.0001
  y = transmissionloss(pm, AcousticSource(0.0, -d, 1000.0), AcousticReceiver(100.0, -d); mode=:incoherent)
  @test -10 * log10(abs2(x)) ≈ y atol=0.1
  x1 = transfercoef(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiver(100.0, -5.0))
  x2 = transfercoef(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiver(100.0, -10.0))
  x3 = transfercoef(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiver(100.0, -15.0))
  x = transfercoef(pm, AcousticSource(0.0, -5.0, 1000.0), [AcousticReceiver(100.0, -d) for d ∈ 5.0:5.0:15.0])
  @test x isa AbstractVector
  @test [x1, x2, x3] == x
  x = transfercoef(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiverGrid2D(100.0, 0.0, 1, -5.0, -5.0, 3))
  @test x isa AbstractMatrix
  @test size(x) == (1, 3)
  @test [x1 x2 x3] == x
  x = transfercoef(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiverGrid2D(100.0, 10.0, 3, -5.0, -5.0, 3))
  @test x isa AbstractMatrix
  @test size(x) == (3, 3)
  @test [x1, x2, x3] == x[1,:]
  x1 = transmissionloss(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiver(100.0, -5.0))
  x2 = transmissionloss(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiver(100.0, -10.0))
  x3 = transmissionloss(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiver(100.0, -15.0))
  x = transmissionloss(pm, AcousticSource(0.0, -5.0, 1000.0), [AcousticReceiver(100.0, -d) for d ∈ 5.0:5.0:15.0])
  @test x isa AbstractVector
  @test [x1, x2, x3] == x
  x = transmissionloss(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiverGrid2D(100.0, 0.0, 1, -5.0, -5.0, 3))
  @test x isa AbstractMatrix
  @test size(x) == (1, 3)
  @test [x1 x2 x3] == x
  x = transmissionloss(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiverGrid2D(100.0, 10.0, 3, -5.0, -5.0, 3))
  @test x isa AbstractMatrix
  @test size(x) == (3, 3)
  @test [x1, x2, x3] == x[1,:]

  env = UnderwaterEnvironment()
  pm = PekerisRayModel(env, 7)
  tx = AcousticSource(0.0, -5.0, 1000.0)
  rx = AcousticReceiver(100.0, -10.0)
  sig = record(pm, tx, rx, 1.0, 44100.0)
  @test size(sig) == (44100,)
  sig = record(pm, tx, rx, 1.0, 44100.0; start=0.5)
  @test size(sig) == (44100,)
  sig = recorder(pm, tx, rx)(1.0, 44100.0)
  @test size(sig) == (44100,)
  sig = recorder(pm, tx, rx)(1.0, 44100.0; start=0.5)
  @test size(sig) == (44100,)
  tx = [AcousticSource(0.0, -5.0, 1000.0), AcousticSource(0.0, -10.0, 2000.0)]
  rx = AcousticReceiver(100.0, -10.0)
  sig = record(pm, tx, rx, 1.0, 44100.0)
  @test size(sig) == (44100,)
  sig = record(pm, tx, rx, 1.0, 44100.0; start=0.5)
  @test size(sig) == (44100,)
  sig = recorder(pm, tx, rx)(1.0, 44100.0)
  @test size(sig) == (44100,)
  sig = recorder(pm, tx, rx)(1.0, 44100.0; start=0.5)
  @test size(sig) == (44100,)
  tx = [AcousticSource(0.0, -5.0, 1000.0), AcousticSource(0.0, -10.0, 2000.0)]
  rx = [AcousticReceiver(100.0, -10.0), AcousticReceiver(100.0, -15.0)]
  sig = record(pm, tx, rx, 1.0, 44100.0)
  @test size(sig) == (44100,2)
  sig = record(pm, tx, rx, 1.0, 44100.0; start=0.5)
  @test size(sig) == (44100,2)
  sig = recorder(pm, tx, rx)(1.0, 44100.0)
  @test size(sig) == (44100,2)
  sig = recorder(pm, tx, rx)(1.0, 44100.0; start=0.5)
  @test size(sig) == (44100,2)
  tx = AcousticSource(0.0, -5.0, 1000.0)
  rx = [AcousticReceiver(100.0, -10.0), AcousticReceiver(100.0, -15.0)]
  sig = record(pm, tx, rx, 1.0, 44100.0)
  @test size(sig) == (44100,2)
  sig = record(pm, tx, rx, 1.0, 44100.0; start=0.5)
  @test size(sig) == (44100,2)
  sig = recorder(pm, tx, rx)(1.0, 44100.0)
  @test size(sig) == (44100,2)
  sig = recorder(pm, tx, rx)(1.0, 44100.0; start=0.5)
  @test size(sig) == (44100,2)

  env = UnderwaterEnvironment(noise=missing)
  pm = PekerisRayModel(env, 7)
  tx = Pinger(0.0, -5.0, 1000.0; interval=0.3)
  rx = AcousticReceiver(100.0, -10.0)
  sig1 = record(pm, tx, rx, 1.0, 44100.0)
  @test size(sig1) == (44100,)
  sig2 = record(pm, tx, rx, 1.0, 44100.0; start=0.5)
  @test size(sig2) == (44100,)
  @test sig1[22051:end] ≈ sig2[1:22050]
  rx = AcousticReceiver(100.0, -11.0)
  sig3 = record(pm, tx, rx, 1.0, 44100.0)
  @test size(sig3) == (44100,)
  @test !(sig1 ≈ sig3)
  rx = AcousticReceiver(100.0/√2, 100.0/√2, -10.0)
  sig3 = record(pm, tx, rx, 1.0, 44100.0)
  @test size(sig3) == (44100,)
  @test sig1 ≈ sig3
  tx = [Pinger(0.0, -5.0, 1000.0; interval=0.3), Pinger(1.0, -5.0, 2000.0; interval=0.5)]
  rx = AcousticReceiver(100.0, 0.0, -10.0)
  sig1 = record(pm, tx, rx, 1.0, 44100.0)
  rx = AcousticReceiver(100.0/√2, 100.0/√2, -10.0)
  rx = AcousticReceiver(-100.0, 0.0, -10.0)
  sig2 = record(pm, tx, rx, 1.0, 44100.0)
  @test !(sig1 ≈ sig2)

end

@testset "pm-bellhop" begin

  if Bellhop in models()

    env = UnderwaterEnvironment(seasurface=Vacuum)
    pm = Bellhop(env)
    @test pm isa Bellhop

    arr = arrivals(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiver(100.0, -10.0))
    @test arr isa AbstractArray{<:UnderwaterAcoustics.RayArrival}
    r = eigenrays(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiver(100.0, -10.0))
    @test r isa AbstractArray{<:UnderwaterAcoustics.RayArrival}
    r = rays(pm, AcousticSource(0.0, -5.0, 1000.0), -60°:15°:60°, 100.0)
    @test r isa AbstractArray{<:UnderwaterAcoustics.RayArrival}
    x = transfercoef(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiver(100.0, -5.0))
    @test x isa Complex
    y = transmissionloss(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiver(100.0, -5.0))
    @test -10 * log10(abs2(x)) ≈ y atol=0.1
    x = transfercoef(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiver(100.0, -5.0); mode=:incoherent)
    @test x isa Complex
    @test imag(x) == 0.0
    y = transmissionloss(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiver(100.0, -5.0); mode=:incoherent)
    @test -10 * log10(abs2(x)) ≈ y atol=0.1
    x1 = transfercoef(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiver(100.0, -5.0))
    x2 = transfercoef(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiver(100.0, -10.0))
    x3 = transfercoef(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiver(100.0, -15.0))
    x = transfercoef(pm, AcousticSource(0.0, -5.0, 1000.0), [AcousticReceiver(100.0, -d) for d ∈ 5.0:5.0:15.0])
    @test x isa AbstractVector
    @test [x1, x2, x3] == x
    x = transfercoef(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiverGrid2D(100.0, 0.0, 1, -5.0, -5.0, 3))
    @test x isa AbstractMatrix
    @test size(x) == (1, 3)
    @test [x1 x2 x3] == x
    x = transfercoef(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiverGrid2D(100.0, 10.0, 3, -5.0, -5.0, 3))
    @test x isa AbstractMatrix
    @test size(x) == (3, 3)
    @test [x1, x2, x3] == x[1,:]
    x1 = transmissionloss(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiver(100.0, -5.0))
    x2 = transmissionloss(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiver(100.0, -10.0))
    x3 = transmissionloss(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiver(100.0, -15.0))
    x = transmissionloss(pm, AcousticSource(0.0, -5.0, 1000.0), [AcousticReceiver(100.0, -d) for d ∈ 5.0:5.0:15.0])
    @test x isa AbstractVector
    @test [x1, x2, x3] == x
    x = transmissionloss(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiverGrid2D(100.0, 0.0, 1, -5.0, -5.0, 3))
    @test x isa AbstractMatrix
    @test size(x) == (1, 3)
    @test [x1 x2 x3] == x
    x = transmissionloss(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiverGrid2D(100.0, 10.0, 3, -5.0, -5.0, 3))
    @test x isa AbstractMatrix
    @test size(x) == (3, 3)
    @test [x1, x2, x3] == x[1,:]

    tx = AcousticSource(0.0, -5.0, 1000.0)
    rx = AcousticReceiver(100.0, -10.0)
    sig = record(pm, tx, rx, 1.0, 44100.0)
    @test size(sig) == (44100,)
    sig = record(pm, tx, rx, 1.0, 44100.0; start=0.5)
    @test size(sig) == (44100,)
    sig = recorder(pm, tx, rx)(1.0, 44100.0)
    @test size(sig) == (44100,)
    sig = recorder(pm, tx, rx)(1.0, 44100.0; start=0.5)
    @test size(sig) == (44100,)
    tx = [AcousticSource(0.0, -5.0, 1000.0), AcousticSource(0.0, -10.0, 2000.0)]
    rx = AcousticReceiver(100.0, -10.0)
    sig = record(pm, tx, rx, 1.0, 44100.0)
    @test size(sig) == (44100,)
    sig = record(pm, tx, rx, 1.0, 44100.0; start=0.5)
    @test size(sig) == (44100,)
    sig = recorder(pm, tx, rx)(1.0, 44100.0)
    @test size(sig) == (44100,)
    sig = recorder(pm, tx, rx)(1.0, 44100.0; start=0.5)
    @test size(sig) == (44100,)
    tx = [AcousticSource(0.0, -5.0, 1000.0), AcousticSource(0.0, -10.0, 2000.0)]
    rx = [AcousticReceiver(100.0, -10.0), AcousticReceiver(100.0, -15.0)]
    sig = record(pm, tx, rx, 1.0, 44100.0)
    @test size(sig) == (44100,2)
    sig = record(pm, tx, rx, 1.0, 44100.0; start=0.5)
    @test size(sig) == (44100,2)
    sig = recorder(pm, tx, rx)(1.0, 44100.0)
    @test size(sig) == (44100,2)
    sig = recorder(pm, tx, rx)(1.0, 44100.0; start=0.5)
    @test size(sig) == (44100,2)
    tx = AcousticSource(0.0, -5.0, 1000.0)
    rx = [AcousticReceiver(100.0, -10.0), AcousticReceiver(100.0, -15.0)]
    sig = record(pm, tx, rx, 1.0, 44100.0)
    @test size(sig) == (44100,2)
    sig = record(pm, tx, rx, 1.0, 44100.0; start=0.5)
    @test size(sig) == (44100,2)
    sig = recorder(pm, tx, rx)(1.0, 44100.0)
    @test size(sig) == (44100,2)
    sig = recorder(pm, tx, rx)(1.0, 44100.0; start=0.5)
    @test size(sig) == (44100,2)

    env = UnderwaterEnvironment(
      seasurface=Vacuum,
      ssp=SampledSSP(0.0:5.0:20.0, [1500.0, 1490.0, 1500.0, 1505.0, 1507.0]),
      altimetry=SampledAltitude(0.0:25.0:100.0, [0.0, -1.0, 0.0, -1.0, 0.0]),
      bathymetry=SampledDepth(0.0:25.0:100.0, [20.0, 17.0, 17.0, 19.0, 20.0])
    )
    pm = Bellhop(env)
    @test pm isa Bellhop
    r = eigenrays(pm, AcousticSource(0.0, -5.0, 5000.0), AcousticReceiver(100.0, -10.0))
    @test r isa AbstractArray{<:UnderwaterAcoustics.RayArrival}
    x = transmissionloss(pm, AcousticSource(0.0, -5.0, 5000.0), AcousticReceiverGrid2D(1.0, 1.0, 100, 0.0, -1.0, 20))
    @test x isa AbstractMatrix
    @test size(x) == (100, 20)

    struct TestAlt <: Altimetry end
    UnderwaterAcoustics.altitude(::TestAlt, x, y) = -1.0 + sin(2π*x/10.0)

    struct TestBathy <: Bathymetry end
    UnderwaterAcoustics.depth(::TestBathy, x, y) = 18.0 + 2*sin(2π*x/30.0)
    UnderwaterAcoustics.maxdepth(::TestBathy) = 20.0

    env = UnderwaterEnvironment(
      seasurface=Vacuum,
      ssp=MunkSSP(),
      altimetry=TestAlt(),
      bathymetry=TestBathy()
    )
    pm = Bellhop(env)
    @test pm isa Bellhop
    r = eigenrays(pm, AcousticSource(0.0, -5.0, 5000.0), AcousticReceiver(100.0, -10.0))
    @test r isa AbstractArray{<:UnderwaterAcoustics.RayArrival}
    x = transmissionloss(pm, AcousticSource(0.0, -5.0, 5000.0), AcousticReceiverGrid2D(1.0, 1.0, 100, 0.0, -1.0, 20))
    @test x isa AbstractMatrix
    @test size(x) == (100, 20)

  else
    @test_skip true
  end

end

@testset "pm-raysolver" begin

  @test RaySolver in models()

  env = UnderwaterEnvironment()
  pm = RaySolver(env)
  @test pm isa RaySolver

  arr = arrivals(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiver(100.0, -10.0))
  @test arr isa AbstractArray{<:UnderwaterAcoustics.RayArrival}
  @test length(arr) >= 7
  @test arr[1].time ≈ 0.0650 atol=0.0001
  @test arr[2].time ≈ 0.0657 atol=0.0001
  @test arr[3].time ≈ 0.0670 atol=0.0001
  @test all([arr[j].time > arr[j-1].time for j ∈ 2:7])
  @test abs(arr[1].phasor) ≈ 0.01 atol=0.001
  @test real(arr[2].phasor) < 0.0
  @test imag(arr[2].phasor) ≈ 0.0
  @test all([abs(arr[j].phasor) < abs(arr[j-1].phasor) for j ∈ 2:7])
  @test [(arr[j].surface, arr[j].bottom) for j ∈ 1:7] == [(0,0), (1,0), (0,1), (1,1), (1,1), (2,1), (1,2)]
  @test abs.([a.arrivalangle for a ∈ arr]) ≈ abs.([a.launchangle for a ∈ arr])

  r = eigenrays(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiver(100.0, -10.0))
  @test r isa AbstractArray{<:UnderwaterAcoustics.RayArrival}
  @test length(r) >= 7
  @test abs.([a.arrivalangle for a ∈ r]) ≈ abs.([a.launchangle for a ∈ r])
  @test all([r[j].raypath[1] == (0.0, 0.0, -5.0) for j ∈ 1:7])
  #@test all([r[j].raypath[end][k] .≈ (100.0, 0.0, -10.0)[k] for j ∈ 1:7, k ∈ 1:3])

  r = rays(pm, AcousticSource(0.0, -5.0, 1000.0), -60°:15°:60°, 100.0)
  @test r isa AbstractArray{<:UnderwaterAcoustics.RayArrival}
  @test length(r) == 9
  @test all([r[j].launchangle for j ∈ 1:9] .≈ -60°:15°:60°)
  @test abs.([a.arrivalangle for a ∈ r]) ≈ abs.([a.launchangle for a ∈ r])
  @test all([r[j].raypath[1] == (0.0, 0.0, -5.0) for j ∈ 1:9])
  @test r[4].raypath[end][1] ≥ 99.9
  @test r[5].raypath[end][1] ≥ 99.9
  @test r[6].raypath[end][1] ≥ 99.9
  @test r[7].raypath[end][1] ≥ 99.9

  env = UnderwaterEnvironment(ssp=IsoSSP(1500.0), seabed=RayleighReflectionCoef(1.0, 1.0))
  pm = RaySolver(env)
  d = (√1209.0)/4.0
  x = transfercoef(pm, AcousticSource(0.0, -d, 1000.0), AcousticReceiver(100.0, -d))
  @test x isa Complex
  @test abs(x) ≈ 0.0 atol=0.0002
  x′ = transfercoef(pm, AcousticSource(0.0, -d, 1000.0), AcousticReceiver(100.0, -d); mode=:incoherent)
  @test x′ isa Complex
  @test imag(x′) == 0.0
  @test abs(x′) > 1/100.0
  d = (√2409.0)/8.0
  x = transfercoef(pm, AcousticSource(0.0, -d, 1000.0), AcousticReceiver(100.0, -d))
  @test abs(x) > abs(x′)
  y = transmissionloss(pm, AcousticSource(0.0, -d, 1000.0), AcousticReceiver(100.0, -d))
  @test -10 * log10(abs2(x)) ≈ y atol=0.1
  x = transfercoef(pm, AcousticSource(0.0, -d, 1000.0), AcousticReceiver(100.0, -d); mode=:incoherent)
  @test abs(x) ≈ abs(x′) atol=0.0001
  y = transmissionloss(pm, AcousticSource(0.0, -d, 1000.0), AcousticReceiver(100.0, -d); mode=:incoherent)
  @test -10 * log10(abs2(x)) ≈ y atol=0.1
  x1 = transfercoef(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiver(100.0, -5.0))
  x2 = transfercoef(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiver(100.0, -10.0))
  x3 = transfercoef(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiver(100.0, -15.0))
  x = transfercoef(pm, AcousticSource(0.0, -5.0, 1000.0), [AcousticReceiver(100.0, -d) for d ∈ 5.0:5.0:15.0])
  @test x isa AbstractVector
  @test [x1, x2, x3] == x
  x = transfercoef(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiverGrid2D(100.0, 0.0, 1, -5.0, -5.0, 3))
  @test x isa AbstractMatrix
  @test size(x) == (1, 3)
  @test [x1 x2 x3] ≈ x atol=0.01
  y = transfercoef(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiverGrid2D(100.0, 10.0, 3, -5.0, -5.0, 3))
  @test y isa AbstractMatrix
  @test size(y) == (3, 3)
  @test x' ≈ y[1,:] atol=0.05
  x1 = transmissionloss(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiver(100.0, -5.0))
  x2 = transmissionloss(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiver(100.0, -10.0))
  x3 = transmissionloss(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiver(100.0, -15.0))
  x = transmissionloss(pm, AcousticSource(0.0, -5.0, 1000.0), [AcousticReceiver(100.0, -d) for d ∈ 5.0:5.0:15.0])
  @test x isa AbstractVector
  @test [x1, x2, x3] == x
  x = transmissionloss(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiverGrid2D(100.0, 0.0, 1, -5.0, -5.0, 3))
  @test x isa AbstractMatrix
  @test size(x) == (1, 3)
  @test [x1, x2] ≈ x[1:2] atol=1.5
  y = transmissionloss(pm, AcousticSource(0.0, -5.0, 1000.0), AcousticReceiverGrid2D(100.0, 10.0, 3, -5.0, -5.0, 3))
  @test y isa AbstractMatrix
  @test size(y) == (3, 3)
  @test x' ≈ y[1,:] atol=0.1

  tx = AcousticSource(0.0, -5.0, 1000.0)
  rx = AcousticReceiver(100.0, -10.0)
  sig = record(pm, tx, rx, 1.0, 44100.0)
  @test size(sig) == (44100,)
  sig = record(pm, tx, rx, 1.0, 44100.0; start=0.5)
  @test size(sig) == (44100,)
  sig = recorder(pm, tx, rx)(1.0, 44100.0)
  @test size(sig) == (44100,)
  sig = recorder(pm, tx, rx)(1.0, 44100.0; start=0.5)
  @test size(sig) == (44100,)
  tx = [AcousticSource(0.0, -5.0, 1000.0), AcousticSource(0.0, -10.0, 2000.0)]
  rx = AcousticReceiver(100.0, -10.0)
  sig = record(pm, tx, rx, 1.0, 44100.0)
  @test size(sig) == (44100,)
  sig = record(pm, tx, rx, 1.0, 44100.0; start=0.5)
  @test size(sig) == (44100,)
  sig = recorder(pm, tx, rx)(1.0, 44100.0)
  @test size(sig) == (44100,)
  sig = recorder(pm, tx, rx)(1.0, 44100.0; start=0.5)
  @test size(sig) == (44100,)
  tx = [AcousticSource(0.0, -5.0, 1000.0), AcousticSource(0.0, -10.0, 2000.0)]
  rx = [AcousticReceiver(100.0, -10.0), AcousticReceiver(100.0, -15.0)]
  sig = record(pm, tx, rx, 1.0, 44100.0)
  @test size(sig) == (44100,2)
  sig = record(pm, tx, rx, 1.0, 44100.0; start=0.5)
  @test size(sig) == (44100,2)
  sig = recorder(pm, tx, rx)(1.0, 44100.0)
  @test size(sig) == (44100,2)
  sig = recorder(pm, tx, rx)(1.0, 44100.0; start=0.5)
  @test size(sig) == (44100,2)
  tx = AcousticSource(0.0, -5.0, 1000.0)
  rx = [AcousticReceiver(100.0, -10.0), AcousticReceiver(100.0, -15.0)]
  sig = record(pm, tx, rx, 1.0, 44100.0)
  @test size(sig) == (44100,2)
  sig = record(pm, tx, rx, 1.0, 44100.0; start=0.5)
  @test size(sig) == (44100,2)
  sig = recorder(pm, tx, rx)(1.0, 44100.0)
  @test size(sig) == (44100,2)
  sig = recorder(pm, tx, rx)(1.0, 44100.0; start=0.5)
  @test size(sig) == (44100,2)

end

function ∂(f, x, i, ϵ)
  x1 = copy(x)
  x1[i] = x[i] + ϵ
  f1 = f(x1)
  x1[i] = x[i] - ϵ
  (f1 - f(x1)) / (2ϵ)
end

@testset "pm-∂pekeris" begin

  function ℳ(x)
    D, R, d1, d2, f, c = x
    env = UnderwaterEnvironment(ssp=IsoSSP(c), bathymetry=ConstantDepth(D))
    pm = PekerisRayModel(env, 7)
    transmissionloss(pm, AcousticSource(0.0, -d1, f), AcousticReceiver(R, -d2))
  end

  x = [20.0, 100.0, 5.0, 10.0, 5000.0, 1500.0]
  ∇ℳ = ForwardDiff.gradient(ℳ, x)
  for i ∈ 1:length(x)
    @test ∇ℳ[i] ≈ ∂(ℳ, x, i, 0.0001) atol=0.1
  end

  x = [25.0, 200.0, 10.0, 8.0, 1000.0, 1540.0]
  ∇ℳ = ForwardDiff.gradient(ℳ, x)
  for i ∈ 1:length(x)
    @test ∇ℳ[i] ≈ ∂(ℳ, x, i, 0.0001) atol=0.1
  end

end

@testset "pm-∂raysolver" begin

  function ℳ₁(x)
    D, R, d1, d2, f, c = x
    env = UnderwaterEnvironment(ssp=IsoSSP(c), bathymetry=ConstantDepth(D))
    pm = RaySolver(env)
    transmissionloss(pm, AcousticSource(0.0, -d1, f), AcousticReceiver(R, -d2))
  end

  function ℳ₂(x)
    D, R, d1, d2, f, c = x
    env = UnderwaterEnvironment(ssp=IsoSSP(c), bathymetry=ConstantDepth(D))
    pm = RaySolver(env)
    transmissionloss(pm, AcousticSource(0.0, -d1, f), AcousticReceiverGrid2D(R, 0.0, 1, -d2, 0.0, 1))[1,1]
  end

  x = [20.0, 100.0, 5.0, 10.0, 5000.0, 1500.0]
  ∇ℳ = ForwardDiff.gradient(ℳ₁, x)
  for i ∈ 1:length(x)
    # skip i = 2 because it is not yet supported
    i != 2 && @test ∇ℳ[i] ≈ ∂(ℳ₁, x, i, 0.0001) atol=0.1
  end

  x = [25.0, 200.0, 10.0, 8.0, 1000.0, 1540.0]
  ∇ℳ = ForwardDiff.gradient(ℳ₁, x)
  for i ∈ 1:length(x)
    # skip i = 2 because it is not yet supported
    i != 2 && @test ∇ℳ[i] ≈ ∂(ℳ₁, x, i, 0.0001) atol=0.1
  end

  x = [20.0, 100.0, 5.0, 10.0, 5000.0, 1500.0]
  ∇ℳ = ForwardDiff.gradient(ℳ₂, x)
  for i ∈ 1:length(x)
    @test ∇ℳ[i] ≈ ∂(ℳ₂, x, i, 0.0001) atol=0.1
  end

  x = [25.0, 200.0, 10.0, 8.0, 1000.0, 1540.0]
  ∇ℳ = ForwardDiff.gradient(ℳ₂, x)
  for i ∈ 1:length(x)
    @test ∇ℳ[i] ≈ ∂(ℳ₂, x, i, 0.0001) atol=0.1
  end

end


text:
module GeometricBase

    using Base: Callable, TwicePrecision

    using Test

    include("Config.jl")
    include("Utils.jl")

    export evaluate, evaluate!

    function evaluate end
    function evaluate! end

    export reset!, cut_periodic_solution!

    function reset! end
    function cut_periodic_solution! end

    export ntime, nsave, nsamples, nconstraints

    function ntime end
    function nsave end
    function nsamples end
    function nconstraints end

    export eachsample, eachtimestep, timesteps

    function eachsample end
    function eachtimestep end
    function timesteps end

    export invariants, parameters, periodicity

    function invariants end
    function parameters end
    function periodicity end

    include("types.jl")

    export NullInvariants, NullParameters, NullPeriodicity

    export OptionalAbstractArray, OptionalArray,
           OptionalCallable, OptionalFunction,
           OptionalNamedTuple, OptionalTuple

    export OptionalInvariants,
           OptionalParameters,
           OptionalPeriodicity

    export AbstractData

    export State, StateVector, SolutionVector

    include("solutions/dataseries.jl")

    export get_data!, set_data!
    export AbstractDataSeries, DataSeries,
           DataSeriesConstructor

    include("solutions/timeseries.jl")

    export TimeSeries, compute_timeseries!

    include("solutions/solution.jl")

    export Solution
    export counter, offset, lastentry

end


text:
using WENO4
using Test

@testset "WENO4" begin
    x = 0:0.2:20
    y = sin.(x)
    @test interpolate_weno4(x, x, y) ≈ y  atol=1e-10
    f(x) = 2.0 + 3.5 * x - 6.9 * x^2
    y = f.(x)
    @test interpolate_weno4(x, x, y) ≈ y  atol=1e-10
    # Constant edges
    @test interpolate_weno4([-1.0], x, y)[1] == y[1]
    @test interpolate_weno4([21.0], x, y)[1] == y[end]
    # Extrapolated edges, exact solution
    @test interpolate_weno4([-10.0], x, y, extrapolate=true)[1] ≈ f(-10) atol=1e-10
    @test interpolate_weno4([50.0], x, y, extrapolate=true)[1] ≈ f(50.) atol=1e-10
    # Zig-zag profile with enough points to fit polynomial, interpolation should
    # never be beyond extremes
    x = [0., 0.5, 0.75, 1., 1+1e-10, 1.25, 1.5, 2]
    y = [0, 0.2, 0.75, 1., -1, -0.75, -0.2, 0]
    new_x = 0.8:0.01:1.2
    new_y = interpolate_weno4(new_x, x, y)
    @test maximum(new_y) <= maximum(y)
    @test minimum(new_y) >= minimum(y)
    a = 1:0.1:2
    @test WENO4.binary_search(a, 1.25) == findfirst(x -> x>1.25, a) - 1
    @test WENO4.weno4_q(1, 1, 2, 3, 4, 1, 2.5, 3.5, 4.5) == (1.0, 1.5)
    @test WENO4.weno4_q(1, 1, 2, 3, 4, 0, 0, 0, 0) == (0.0, 0.0)
    @test WENO4.weno4_β(1, 2, 3, 4, 0, 0, 0, 0.5) == (0.0, 1.0)
    @test WENO4.weno4_β(1, 2, 3, 4, 0, 0, 0, 0) == (0.0, 0.0)
end