text:
refVals = []
refPrecs = []

#! format: off
# SC ========== Test number 1 reference values and precision match template. =======
# SC ========== /home/jmc/cliMa/cliMa_new_jmc/test/Ocean/SplitExplicit/simple_dbl_gyre.jl test reference values ======================================
# BEGIN SCPRINT
# varr - reference values (from reference run)
# parr - digits match precision (hand edit as needed)
#
# [
#  [ MPIStateArray Name, Field Name, Maximum, Minimum, Mean, Standard Deviation ],
#  [         :                :          :        :      :          :           ],
# ]
varr = [
 [  "oce Q_3D",   "u[1]", -1.41235234130601461366e-01,  2.44137911410471170059e-01, -1.17638792697003554538e-02,  4.76422973396204429974e-02 ],
 [  "oce Q_3D",   "u[2]", -3.71526973402696303328e-01,  2.63129391123315958811e-01, -3.84758873236049556144e-02,  4.34631279346227028526e-02 ],
 [  "oce Q_3D",      "η", -3.62315954238604787108e+00,  3.25374835544518203889e+00, -1.15530630772209001104e-04,  1.69499420692026547819e+00 ],
 [  "oce Q_3D",      "θ", -1.95989738813091111946e-02,  2.41980669630830433903e+01,  6.00291945213983790808e+00,  5.36502088796993170661e+00 ],
 [   "oce aux",      "w", -1.42436775039925483770e-03,  1.17093883909037276177e-03, -4.27143547438785408186e-07,  1.58700821328606480297e-04 ],
 [   "oce aux",   "pkin", -6.47863971232045656734e+01,  0.00000000000000000000e+00, -2.35470923626465094003e+01,  1.84884052655973327717e+01 ],
 [   "oce aux",    "wz0", -2.03056757204373681007e-04,  2.51217746755789830913e-04, -1.04288746719621225719e-08,  9.09212328801199519022e-05 ],
 [   "oce aux", "u_d[1]", -1.29380048170506856131e-01,  2.41847160349357936937e-01,  1.50957373152410783811e-04,  4.55063940872466113352e-02 ],
 [   "oce aux", "u_d[2]", -2.31151064740037659462e-01,  3.17638906635839823878e-01,  7.38946153332456550028e-05,  3.35035482074258969543e-02 ],
 [   "oce aux", "ΔGu[1]", -2.28025960163159889874e-05,  3.46723882659755934413e-06, -5.56600796147506152338e-07,  2.61685583178905381934e-06 ],
 [   "oce aux", "ΔGu[2]", -3.00923384704360015996e-06,  1.24960695834485056166e-05,  6.54681130814259427502e-06,  3.13708534178777629812e-06 ],
 [   "oce aux",      "y",  0.00000000000000000000e+00,  6.00000000000000093132e+06,  3.00000000000000000000e+06,  1.73273876310492726043e+06 ],
 [ "baro Q_2D",   "U[1]", -1.37905164694259781299e+02,  1.09296887114357033965e+02, -3.59834173192416955089e+01,  4.24653520384661575804e+01 ],
 [ "baro Q_2D",   "U[2]", -4.93706839547863808093e+02,  5.71157285218092027890e+01, -1.15735478409127139798e+02,  8.23847078816430951065e+01 ],
 [ "baro Q_2D",      "η", -3.62624136658312767878e+00,  3.25919860855558907176e+00, -1.15147341274860083972e-04,  1.69658109382101818241e+00 ],
 [  "baro aux",  "Gᵁ[1]", -1.04017164797926778969e-02,  6.84077880489479678294e-02,  1.66980238844251878058e-03,  7.85082570477715034618e-03 ],
 [  "baro aux",  "Gᵁ[2]", -3.74882087503455169175e-02,  9.02770154113080071367e-03, -1.96404339244277831300e-02,  9.41156556666298931002e-03 ],
 [  "baro aux", "U_c[1]", -1.36958883839043295438e+02,  1.09074110263089167461e+02, -3.57794851972332565992e+01,  4.23235946471143051895e+01 ],
 [  "baro aux", "U_c[2]", -4.92966277666792962009e+02,  5.72050394052077919582e+01, -1.15636183420627688179e+02,  8.22648167942200103653e+01 ],
 [  "baro aux",    "η_c", -3.62315954238604787108e+00,  3.25374835544518203889e+00, -1.15530630772203336148e-04,  1.69504995619627174541e+00 ],
 [  "baro aux", "U_s[1]", -1.37905164694259781299e+02,  1.09296887114357033965e+02, -3.59834173192416955089e+01,  4.24653520384661575804e+01 ],
 [  "baro aux", "U_s[2]", -4.93706839547863808093e+02,  5.71157285218092027890e+01, -1.15735478409127139798e+02,  8.23847078816430951065e+01 ],
 [  "baro aux",    "η_s", -3.62624136658312767878e+00,  3.25919860855558907176e+00, -1.15147341274860083972e-04,  1.69658109382101818241e+00 ],
 [  "baro aux",  "Δu[1]", -4.66680666769807577128e-04,  1.05054947340137826844e-02,  3.56177517668052924515e-03,  2.34394235403509801352e-03 ],
 [  "baro aux",  "Δu[2]", -2.98072019929583121103e-03,  4.75363215770287662193e-03,  1.22307915357422513150e-03,  1.31189027778619380499e-03 ],
 [  "baro aux", "η_diag", -3.62078008618114477457e+00,  3.25151939956774960194e+00, -1.15589924206263583179e-04,  1.69498993030730304987e+00 ],
 [  "baro aux",     "Δη", -8.83137470915729139165e-03,  5.97876216292370088468e-03,  5.92934342415881362251e-08,  2.08693123903101827171e-03 ],
 [  "baro aux",      "y",  0.00000000000000000000e+00,  6.00000000000000093132e+06,  3.00000000000000000000e+06,  1.73279575381979602389e+06 ],
]
parr = [
 [  "oce Q_3D",   "u[1]",    12,    12,    12,    12 ],
 [  "oce Q_3D",   "u[2]",    12,    12,    12,    12 ],
 [  "oce Q_3D",      "η",    12,    12,     8,    12 ],
 [  "oce Q_3D",      "θ",    12,    12,    12,    12 ],
 [   "oce aux",      "w",    12,    12,     8,    12 ],
 [   "oce aux",   "pkin",    12,    12,    12,    12 ],
 [   "oce aux",    "wz0",    12,    12,     8,    12 ],
 [   "oce aux", "u_d[1]",    12,    12,    12,    12 ],
 [   "oce aux", "u_d[2]",    12,    12,    12,    12 ],
 [   "oce aux", "ΔGu[1]",    12,    12,    12,    12 ],
 [   "oce aux", "ΔGu[2]",    12,    12,    12,    12 ],
 [   "oce aux",      "y",    12,    12,    12,    12 ],
 [ "baro Q_2D",   "U[1]",    12,    12,    12,    12 ],
 [ "baro Q_2D",   "U[2]",    12,    12,    12,    12 ],
 [ "baro Q_2D",      "η",    12,    12,     8,    12 ],
 [  "baro aux",  "Gᵁ[1]",    12,    12,    12,    12 ],
 [  "baro aux",  "Gᵁ[2]",    12,    12,    12,    12 ],
 [  "baro aux", "U_c[1]",    12,    12,    12,    12 ],
 [  "baro aux", "U_c[2]",    12,    12,    12,    12 ],
 [  "baro aux",    "η_c",    12,    12,     8,    12 ],
 [  "baro aux", "U_s[1]",    12,    12,    12,    12 ],
 [  "baro aux", "U_s[2]",    12,    12,    12,    12 ],
 [  "baro aux",    "η_s",    12,    12,     8,    12 ],
 [  "baro aux",  "Δu[1]",    12,    12,    12,    12 ],
 [  "baro aux",  "Δu[2]",    12,    12,    12,    12 ],
 [  "baro aux", "η_diag",    12,    12,     8,    12 ],
 [  "baro aux",     "Δη",     9,     9,     6,    10 ],
 [  "baro aux",      "y",    12,    12,    12,    12 ],
]
# END SCPRINT
# SC ====================================================================================

    append!(refVals ,[ varr ] )
    append!(refPrecs,[ parr ] )

#! format: on


text:
using ShapeFactory

# concentric arcs with random angles
# rad   :start radius
# delta :self explanatory
# num   :number of circles
#
#
function random_concentric(rad::Real, delta::Real, num::Integer, idx::Integer = 1)

    app = dispatch("CATIA.Application")
    prt = app.ActiveDocument.Part
    fac = prt.HybridShapeFactory

    rad_dim = prt.Parameters.CreateDimension("min_radius_idx$(idx)", "Length", 0)
    rad_dim.ValuateFromString("$(rad)mm")

    rdt_dim = prt.Parameters.CreateDimension("rad_delta_idx$(idx)", "Length", 0)
    rdt_dim.ValuateFromString("$(delta)mm")

    bdy = prt.HybridBodies.Add()
    bdy.Name = "concentric_circle_idx$(idx)"

    cen = fac.AddNewPointCoord(0, 0, 0)
    pln = prt.FindObjectByName("xy plane")

    rds = rad:delta:(rad + delta * num-1)
    a_a = rand(0:300, num)
    a_b = rand.((:).(a_a .+ 10, 360))

    con = fac.AddNewCircleCtrRadWithAngles.(cen, pln, false, rds, a_a, a_b)

    sub = prt.Relations.SubList(bdy, false)
    exp = x -> "min_radius_idx$(idx) + rad_delta_idx$(idx) * $(x)"
    sub.CreateFormula.("", "", getproperty.(con, :Radius), exp.(1:num))

    bdy.AppendHybridShape.(con)
    prt.Update()
end

ext: jl

max_stars_repo_name: zaika-generative/ShapeFactory

max_issues_repo_name: zaika-generative/ShapeFactory

max_forks_repo_name: zaika-generative/ShapeFactory


text:
using FastTransforms, Test

@testset "Chebyshev transform"  begin
    @testset "Chebyshev points" begin
        @test @inferred(chebyshevpoints(10)) == @inferred(chebyshevpoints(Float64, 10))
        @test @inferred(chebyshevpoints(10, Val(2))) == @inferred(chebyshevpoints(Float64, 10, Val(2)))
        for T in (Float32, Float64, ComplexF32, ComplexF64)
            @test chebyshevpoints(T, 0) == T[]
            @test chebyshevpoints(T, 1) == T[0]

            n = 20
            @test @inferred(chebyshevpoints(T, n)) == [sinpi(convert(T,n-2k+1)/(2n)) for k=1:n]
            @test @inferred(chebyshevpoints(T, n, Val(2))) == [sinpi(convert(T,n-2k+1)/(2n-2)) for k=1:n]

            @test_throws MethodError chebyshevpoints(n, Val(-1))
            @test_throws ArgumentError chebyshevpoints(T, 0, Val(2))
            @test_throws ArgumentError chebyshevpoints(T, 1, Val(2))
        end
    end

    @testset "Chebyshev first kind points <-> first kind coefficients" begin
        for T in (Float32, Float64, ComplexF32, ComplexF64)
            n = 20
            p_1 = chebyshevpoints(T, n)
            f = exp.(p_1)
            f̌ = @inferred(chebyshevtransform(f))
            @test f̌ == chebyshevtransform!(copy(f))

            f̃ = x -> [cos(k*acos(x)) for k=0:n-1]' * f̌
            @test f̃(0.1) ≈ exp(T(0.1))
            @test @inferred(ichebyshevtransform(f̌)) ≈ ichebyshevtransform!(copy(f̌)) ≈ exp.(p_1)

            f̃ = copy(f)
            f̄ = copy(f̌)
            P = @inferred(plan_chebyshevtransform(f))
            @test @inferred(P*f) == f̌
            @test f == f̃
            @test_throws ArgumentError P * T[1,2]
            P = @inferred(plan_chebyshevtransform!(f))
            @test @inferred(P*f) == f̌
            @test f == f̌
            @test_throws ArgumentError P * T[1,2]
            Pi = @inferred(plan_ichebyshevtransform(f̌))
            @test @inferred(Pi*f̌) ≈ f̃
            @test f̌ == f̄
            @test_throws ArgumentError Pi * T[1,2]
            Pi = @inferred(plan_ichebyshevtransform!(f̌))
            @test @inferred(Pi*f̌) ≈ f̃
            @test f̌ ≈ f̃
            @test_throws ArgumentError Pi * T[1,2]

            @test chebyshevtransform(T[1]) == T[1]
            @test ichebyshevtransform(T[1]) == T[1]
            @test chebyshevtransform(T[]) == T[]
            @test ichebyshevtransform(T[]) == T[]
        end
    end
    @testset "Chebyshev second kind points <-> first kind coefficients" begin
        for T in (Float32, Float64, ComplexF32, ComplexF64)
            n = 20
            p_2 = chebyshevpoints(T, n, Val(2))
            f = exp.(p_2)
            f̌ = @inferred(chebyshevtransform(f, Val(2)))
            @test f̌ == chebyshevtransform!(copy(f), Val(2))

            f̃ = x -> [cos(k*acos(x)) for k=0:n-1]' * f̌
            @test f̃(0.1) ≈ exp(T(0.1))
            @test @inferred(ichebyshevtransform(f̌, Val(2))) ≈ ichebyshevtransform!(copy(f̌), Val(2)) ≈ exp.(p_2)

            P = @inferred(plan_chebyshevtransform!(f, Val(2)))
            Pi = @inferred(plan_ichebyshevtransform!(f, Val(2)))
            @test all(@inferred(P \ copy(f)) .=== Pi * copy(f))
            @test all(@inferred(Pi \ copy(f̌)) .=== P * copy(f̌))
            @test f ≈ P \ (P*copy(f)) ≈ P * (P\copy(f)) ≈ Pi \ (Pi*copy(f)) ≈ Pi * (Pi \ copy(f))

            f̃ = copy(f)
            f̄ = copy(f̌)
            P = @inferred(plan_chebyshevtransform(f, Val(2)))
            @test_throws ArgumentError P * T[1,2]
            @test P*f == f̌
            @test f == f̃
            P = @inferred(plan_chebyshevtransform!(f, Val(2)))
            @test P*f == f̌
            @test f == f̌
            @test_throws ArgumentError P * T[1,2]
            Pi = @inferred(plan_ichebyshevtransform(f̌, Val(2)))
            @test Pi*f̌ ≈ f̃
            @test f̌ == f̄
            @test_throws ArgumentError Pi * T[1,2]
            Pi = @inferred(plan_ichebyshevtransform!(f̌, Val(2)))
            @test Pi*f̌ ≈ f̃
            @test f̌ ≈ f̃
            @test_throws ArgumentError Pi * T[1,2]

            @test_throws ArgumentError chebyshevtransform(T[1], Val(2))
            @test_throws ArgumentError ichebyshevtransform(T[1], Val(2))
            @test_throws ArgumentError chebyshevtransform(T[], Val(2))
            @test_throws ArgumentError ichebyshevtransform(T[], Val(2))
        end
    end

    @testset "Chebyshev first kind points <-> second kind coefficients" begin
        for T in (Float32, Float64, ComplexF32, ComplexF64)
            n = 20
            p_1 = chebyshevpoints(T, n)
            f = exp.(p_1)
            f̌ = @inferred(chebyshevutransform(f))

            f̃ = x -> [sin((k+1)*acos(x))/sin(acos(x)) for k=0:n-1]' * f̌
            @test f̃(0.1) ≈ exp(T(0.1))
            @test ichebyshevutransform(f̌) ≈ exp.(p_1)

            f̃ = copy(f)
            f̄ = copy(f̌)
            P = @inferred(plan_chebyshevutransform(f))
            @test P*f == f̌
            @test f == f̃
            @test_throws ArgumentError P * T[1,2]
            P = @inferred(plan_chebyshevutransform!(f))
            @test P*f == f̌
            @test f == f̌
            @test_throws ArgumentError P * T[1,2]
            Pi = @inferred(plan_ichebyshevutransform(f̌))
            @test Pi*f̌ ≈ f̃
            @test f̌ == f̄
            @test_throws ArgumentError Pi * T[1,2]
            Pi = @inferred(plan_ichebyshevutransform!(f̌))
            @test Pi*f̌ ≈ f̃
            @test f̌ ≈ f̃
            @test_throws ArgumentError Pi * T[1,2]

            @test chebyshevutransform(T[1]) == T[1]
            @test ichebyshevutransform(T[1]) == T[1]
            @test chebyshevutransform(T[]) == T[]
            @test ichebyshevutransform(T[]) == T[]
        end
    end

    @testset "Chebyshev second kind points <-> second kind coefficients" begin
        for T in (Float32, Float64, ComplexF32, ComplexF64)
            n = 20
            p_2 = chebyshevpoints(T, n, Val(2))[2:end-1]
            f = exp.(p_2)
            f̌ = @inferred(chebyshevutransform(f, Val(2)))

            f̃ = x -> [sin((k+1)*acos(x))/sin(acos(x)) for k=0:n-3]' * f̌
            @test f̃(0.1) ≈ exp(T(0.1))
            @test @inferred(ichebyshevutransform(f̌, Val(2))) ≈ exp.(p_2)

            f̃ = copy(f)
            f̄ = copy(f̌)
            P = @inferred(plan_chebyshevutransform(f, Val(2)))
            @test @inferred(P*f) == f̌
            @test f == f̃
            @test_throws ArgumentError P * T[1,2]
            P = @inferred(plan_chebyshevutransform!(f, Val(2)))
            @test @inferred(P*f) == f̌
            @test f == f̌
            @test_throws ArgumentError P * T[1,2]
            Pi = @inferred(plan_ichebyshevutransform(f̌, Val(2)))
            @test @inferred(Pi*f̌) ≈ f̃
            @test f̌ == f̄
            @test_throws ArgumentError Pi * T[1,2]
            Pi = @inferred(plan_ichebyshevutransform!(f̌, Val(2)))
            @test @inferred(Pi*f̌) ≈ f̃
            @test f̌ ≈ f̃
            @test_throws ArgumentError Pi * T[1,2]

            @test_throws ArgumentError chebyshevutransform(T[1], Val(2))
            @test_throws ArgumentError ichebyshevutransform(T[1], Val(2))
            @test_throws ArgumentError chebyshevutransform(T[], Val(2))
            @test_throws ArgumentError ichebyshevutransform(T[], Val(2))
        end
    end

    @testset "matrix" begin
        X = randn(1,1)
        @test chebyshevtransform!(copy(X), Val(1)) == ichebyshevtransform!(copy(X), Val(1)) == X
        @test_throws ArgumentError chebyshevtransform!(copy(X), Val(2))
        @test_throws ArgumentError ichebyshevtransform!(copy(X), Val(2))

        X = randn(10,11)

        # manual 2D Chebyshev
        X̌ = copy(X)
        for j in axes(X̌,2)
            chebyshevtransform!(view(X̌,:,j))
        end
        for k in axes(X̌,1)
            chebyshevtransform!(view(X̌,k,:))
        end
        @test chebyshevtransform!(copy(X), Val(1)) ≈ X̌
        @test ichebyshevtransform!(copy(X̌), Val(1)) ≈ X

        # manual 2D Chebyshev
        X̌ = copy(X)
        for j in axes(X̌,2)
            chebyshevtransform!(view(X̌,:,j), Val(2))
        end
        for k in axes(X̌,1)
            chebyshevtransform!(view(X̌,k,:), Val(2))
        end
        @test chebyshevtransform!(copy(X), Val(2)) ≈ X̌
        @test ichebyshevtransform!(copy(X̌), Val(2)) ≈ X
    end

    @testset "Integer" begin
        @test chebyshevtransform([1,2,3]) == chebyshevtransform([1.,2,3])
        @test chebyshevtransform([1,2,3], Val(2)) == chebyshevtransform([1.,2,3], Val(2))
        @test ichebyshevtransform([1,2,3]) == ichebyshevtransform([1.,2,3])
        @test ichebyshevtransform([1,2,3], Val(2)) == ichebyshevtransform([1.,2,3], Val(2))

        @test chebyshevutransform([1,2,3]) == chebyshevutransform([1.,2,3])
        @test chebyshevutransform([1,2,3], Val(2)) == chebyshevutransform([1.,2,3], Val(2))
        @test ichebyshevutransform([1,2,3]) == ichebyshevutransform([1.,2,3])
        @test ichebyshevutransform([1,2,3], Val(2)) == ichebyshevutransform([1.,2,3], Val(2))
    end

    @testset "BigFloat" begin
        x = BigFloat[1,2,3]
        @test ichebyshevtransform(chebyshevtransform(x)) ≈ x
        @test plan_chebyshevtransform(x)x ≈ chebyshevtransform(x)
        @test plan_ichebyshevtransform(x)x ≈ ichebyshevtransform(x)
        @test plan_chebyshevtransform!(x)copy(x) ≈ chebyshevtransform(x)
        @test plan_ichebyshevtransform!(x)copy(x) ≈ ichebyshevtransform(x)
    end
end


text:
import Pkg

path_to_env = joinpath(@__DIR__, "../../")
Pkg.activate(path_to_env)

using CytofRepFAM
using BSON
using Distributions
using LaTeXStrings
using DataFrames
using CSV
using DelimitedFiles
import LinearAlgebra: Diagonal

include(joinpath(@__DIR__, "salso.jl"))
include(joinpath(@__DIR__, "dden_complete.jl"))

import PyCall, PyPlot, Seaborn
const plt = PyPlot.plt
PyPlot.matplotlib.use("Agg")
const sns = Seaborn

#= Interactive plot
PyPlot.matplotlib.use("TkAgg")
=#
#= Non-interactive plot
PyPlot.matplotlib.use("Agg")
=#

# General plot settings
rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
rcParams["font.size"] = 15
rcParams["xtick.labelsize"] = 15
rcParams["ytick.labelsize"] = 15
rcParams["figure.figsize"] = (6, 6)

# TODO
# Steal things from: Cytof5/sims/repfam_fs/test


text:
# parameters that govern counterfactual

remove_shock!(parameters, :global_sectoral_shock_njs)
remove_shock!(parameters, :idiosyncratic_shock_njs)


text:
# Users of this module should first read the document
# https://alan-turing-institute.github.io/MLJ.jl/dev/adding_models_for_general_use/
module MLJBase

export MLJType, Model, Supervised, Unsupervised
export Deterministic, Probabilistic, Interval
export DeterministicNetwork, ProbabilisticNetwork, UnsupervisedNetwork
export fit, update, clean!
export predict, predict_mean, predict_mode, fitted_params
export transform, inverse_transform, se, evaluate, best
export info, info_dict

export load_path, package_url, package_name, package_uuid  # model_traits.jl
export input_scitype, supports_weights                     # model_traits.jl
export target_scitype, output_scitype                      # model_traits.jl
export is_pure_julia, is_wrapper, prediction_type          # model_traits.jl
export params                                        # parameters.jl
export reconstruct, int, decoder, classes            # data.jl
export selectrows, selectcols, select, nrows         # data.jl
export table, levels_seen, matrix, container_type    # data.jl
export partition, unpack                             # data.jl
export @set_defaults                                 # utilities.jl
export @mlj_model                                    # mlj_model_macro.jl
export metadata_model, metadata_pkg                  # metadata_utilities
export HANDLE_GIVEN_ID, @more, @constant             # show.jl
export color_on, color_off                           # show.jl
export UnivariateFinite, average                     # distributions.jl
export SupervisedTask, UnsupervisedTask, MLJTask     # tasks.jl
export X_and_y, X_, y_, nrows, nfeatures             # tasks.jl
export info                                          # info.jl
export @load_boston, @load_ames, @load_iris          # datasets.jl
export @load_reduced_ames                            # datasets.jl
export @load_crabs                                   # datasets.jl
export orientation, reports_each_observation         # measures.jl
export is_feature_dependent                          # measures.jl
export default_measure, value                        # measures.jl
export mav, mae, rms, rmsl, rmslp1, rmsp, l1, l2     # measures.jl
export misclassification_rate, cross_entropy         # measures.jl
export BrierScore                                    # measures.jl

# methods from other packages to be rexported:
export pdf, mean, mode

# re-export of ScientificTypes (`Table` not exported):
export trait
export Scientific, Found, Unknown, Finite, Infinite
export OrderedFactor, Multiclass, Count, Continuous
export Binary, ColorImage, GrayImage, Image
export scitype, scitype_union, coerce, schema

# rexport from Random, Statistics, Distributions, CategoricalArrays:
export pdf, mode, median, mean, shuffle!, categorical, shuffle, levels, levels!
export std

import Base.==
import Base: @__doc__

using Tables
using OrderedCollections # already a dependency of StatsBase
import Distributions
import Distributions: pdf, mode
using CategoricalArrays
using OrderedCollections
import CategoricalArrays
using ScientificTypes
import ScientificTypes: trait

# to be extended:
import StatsBase: fit, predict, fit!
import Missings.levels

# from Standard Library:
using Statistics
using Random
using InteractiveUtils
using SparseArrays

## CONSTANTS

# the directory containing this file:
const srcdir = dirname(@__FILE__)
# horizontal space for field names in `MLJType` object display:
const COLUMN_WIDTH = 24
# how deep to display fields of `MLJType` objects:
const DEFAULT_SHOW_DEPTH = 0

include("utilities.jl")

## BASE TYPES

abstract type MLJType end
include("equality.jl") # equality for MLJType objects

## ABSTRACT MODEL TYPES

# for storing hyperparameters:
abstract type Model <: MLJType end

abstract type Supervised <: Model end
abstract type Unsupervised <: Model end

# supervised models that `predict` probability distributions are of:
abstract type Probabilistic <: Supervised end

# supervised models that `predict` point-values are of:
abstract type Deterministic <: Supervised end

# supervised models that `predict` intervals:
abstract type Interval <: Supervised end

# for models that are "exported" learning networks (return a Node as
# their fit-result; see MLJ docs:
abstract type ProbabilisticNetwork <: Probabilistic end
abstract type DeterministicNetwork <: Deterministic end
abstract type UnsupervisedNetwork <: Unsupervised end

## THE MODEL INTERFACE

# every model interface must implement a `fit` method of the form
# `fit(model, verbosity::Integer, training_args...) -> fitresult, cache, report`
# or, one the simplified versions
# `fit(model, training_args...) -> fitresult`
# `fit(model, X, ys...) -> fitresult`
fit(model::Model, verbosity::Integer, args...) = fit(model, args...), nothing, nothing

# each model interface may optionally overload the following refitting
# method:
update(model::Model, verbosity, fitresult, cache, args...) =
    fit(model, verbosity, args...)

# methods dispatched on a model and fit-result are called
# *operations*.  supervised models must implement a `predict`
# operation (extending the `predict` method of StatsBase).

# unsupervised methods must implement this operation:
function transform end

# unsupervised methods may implement this operation:
function inverse_transform end

# this operation can be optionally overloaded to provide access to
# fitted parameters (eg, coeficients of linear model):
fitted_params(::Model, fitresult) = (fitresult=fitresult,)

# probabilistic supervised models may also overload one or more of
# `predict_mode`, `predict_median` and `predict_mean` defined below.

# mode:
predict_mode(model::Probabilistic, fitresult, Xnew) =
    mode.(predict(model, fitresult, Xnew))

# mean:
predict_mean(model::Probabilistic, fitresult, Xnew) =
    mean.(predict(model, fitresult, Xnew))

# median:
predict_median(model::Probabilistic, fitresult, Xnew) =
    median.(predict(model, fitresult, Xnew))

# operations implemented by some meta-models:
function se end
function evaluate end
function best end

# a model wishing invalid hyperparameters to be corrected with a
# warning should overload this method (return value is the warning
# message):
clean!(model::Model) = ""

## TRAITS

"""

    info(object)

List the traits of an object, such as a model or a performance measure.

"""
info(object) = info(object, Val(ScientificTypes.trait(object)))

include("model_traits.jl")

# for unpacking the fields of MLJ objects:
include("parameters.jl")

# for displaying objects of `MLJType`:
include("show.jl")

# convenience methods for manipulating categorical and tabular data
include("data.jl")

# probability distributions and methods not provided by
# Distributions.jl package:
include("distributions.jl")

include("info.jl")
include("datasets.jl") # importing CSV will also load datasets_requires.jl
include("tasks.jl")
include("measures.jl")

# mlj model macro to help define models
include("mlj_model_macro.jl")

# metadata utils
include("metadata_utilities.jl")

# __init__() function:
include("init.jl")

end # module


text:
function parse_do(ps::ParseState, ret)
    startbyte = ps.nt.startbyte
    start_col = ps.nt.startpos[2] - ret.span + 4

    # Parsing
    next(ps)
    kw = INSTANCE(ps)
    format_kw(ps)

    args = EXPR{TupleH}(EXPR[], - ps.nt.startbyte, Variable[], "")
    @default ps @closer ps comma @closer ps block while !closer(ps)
        @catcherror ps startbyte a = parse_expression(ps)
        push!(args.defs, Variable(Symbol(_arg_id(a).val), get_t(a), args))

        push!(args.args, a)
        if ps.nt.kind == Tokens.COMMA
            next(ps)
            format_comma(ps)
            push!(args.args, INSTANCE(ps))
        end
    end
    args.span += ps.nt.startbyte
    block = EXPR{Block}(EXPR[], 0, Variable[], "")
    @catcherror ps startbyte @default ps parse_block(ps, block, start_col)

    # Construction
    ret = EXPR{Do}(EXPR[ret, kw], ret.span - startbyte, Variable[], "")
    push!(ret.args, args)
    push!(ret.args, block)
    next(ps)
    push!(ret.args, INSTANCE(ps))
    ret.span += ps.nt.startbyte

    return ret
end


text:
	#
	#	Script for calculation of far-field pattern in the liquid case (sphere)
	#

	# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	# User configurable parameters
	# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	# Physical parameters
	rho10 =1/1000 ; # Density ratio
	k_0 = 5 ; # Wave number in media 0
	k_1 = 3 ; # Wave number in media 1
	a = 1 ; # Major semiaxis spheroid
	b = 0.25 ; # Minor semiaxis spheroid
	theta_inc = 2*pi/3 ; # Incidence angle in radians
	theta_incdeg= round(Int,180*theta_inc/pi);
	# Software parameters
	M = 8 ;
	method = 2 ;
	delta_eta = 0.00390625 ; # tuve que hacer lo de poner un paso como 1/2^8, de lo contrario da NaN

	# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	# Calculation
	# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	# Conversion theta to eta
	eta_inc = cos( theta_inc ) ;

	# Coefficient evaluation
	(Amn, Bmn) = Coeff_Obl( rho10, k_0, k_1, a, b, M, method, eta_inc ) ;

	# Far-Field pattern calculation
	Pattern = Pattern_LiquidObl(k_0, a, b, M, method, eta_inc, Amn, delta_eta ) ;

	# Saving to disk
	fileName = string("Obl_a_",a,"_b_",b,"_k0_",k_0,"_k1_",k_1,"_rho1rho0_",rho10,"_inc_",theta_incdeg,".dat");
	writedlm(fileName, Pattern , '\t' ) ;

	# Coefficient evaluation

	 AmnD=Coeff_OblDN( k_0, a, b, M, method, eta_inc, "D")

	# Far-Field pattern calculation
	Pattern = Pattern_LiquidObl( k_0, a, b, M, method, eta_inc, AmnD, delta_eta ) ;

	# Saving to disk
	fileName = string("Obl_a_",a,"_b_",b,"_k0_",k_0,"_k1_",k_1,"_rho1rho0_",rho10,"_inc_",theta_incdeg,"Dirichlet.dat");
	writedlm(fileName, Pattern , '\t' ) ;

	# caso comun
	rho10=3;

	# Coefficient evaluation
	(Amn, Bmn) = Coeff_Obl(rho10, k_0, k_1, a, b, M, method, eta_inc ) ;

	# Far-Field pattern calculation
	Pattern = Pattern_LiquidObl(k_0, a, b, M, method, eta_inc, Amn, delta_eta ) ;

	# Saving to disk
	fileName = string("Obl_a_",a,"_b_",b,"_k0_",k_0,"_k1_",k_1,"_rho1rho0_",rho10,"_inc_",theta_incdeg,".dat");
	writedlm(fileName, Pattern , '\t' ) ;


text:
module TestDualContainers

using Reagents.Internal: @trace
using Test
using ..Utils: run_concurrently

using ArgCheck: @check as @assert
include("../../../examples/dualcontainers.jl")

function test_fifo_data()
    @testset for f in [dual_queue, quack]
        test_fifo_data(f)
    end
end

function test_fifo_data(constructor)
    c = constructor(Int)
    put!(c, 1)
    put!(c, 2)
    put!(c, 3)
    @test [take!(c), take!(c), take!(c)] == [1, 2, 3]
end

function test_lifo_data()
    @testset for f in [dual_stack, steue]
        test_lifo_data(f)
    end
end

function test_lifo_data(constructor)
    c = constructor(Int)
    put!(c, 1)
    put!(c, 2)
    put!(c, 3)
    @test [take!(c), take!(c), take!(c)] == [3, 2, 1]
end

function put_take!(c, items)
    takeorder = Int[]
    tasks = map(1:3) do i
        t = @task begin
            y = take!(c)
            push!(takeorder, i)
            y
        end
        yield(t)
        t
    end
    for x in items
        put!(c, x)
    end
    taken = map(fetch, tasks)
    return (; taken, takeorder)
end

function test_fifo_reservation()
    @testset for f in [dual_queue, steue]
        test_fifo_reservation(f)
    end
end

function test_fifo_reservation(constructor)
    c = constructor(Int)
    (; taken, takeorder) = put_take!(c, [111, 222, 333])
    @test taken == [111, 222, 333]
    @test takeorder == [1, 2, 3]
end

function test_lifo_reservation()
    @testset for f in [dual_stack, quack]
        test_lifo_reservation(f)
    end
end

function test_lifo_reservation(constructor)
    c = constructor(Int)
    (; taken, takeorder) = put_take!(c, [111, 222, 333])
    @test taken == [333, 222, 111]
    @test takeorder == [3, 2, 1]
end

function test_many_items(; kwargs...)
    @testset for spawn in [false, true]
        test_many_items(spawn; kwargs...)
    end
end

function test_many_items(spawn::Bool; kwargs...)
    @testset for f in [dual_stack, dual_queue, quack, steue]
        test_many_items(f, spawn; kwargs...)
    end
end

function test_many_items(constructor, spawn::Bool, nrepeat = 1000)
    c = constructor(Int)

    sentinel = -1  # TODO: support `close` and get rid of this

    received = [Int[] for _ in 1:5]
    receivers = map(eachindex(received), received) do i, dest
        function receiver_task()
            while true
                y = take!(c)
                y == sentinel && break
                push!(dest, y)
            end
            @debug "`test_many_items`: Reciever $i done"
        end
    end

    nsenders = 5
    nfinished = Threads.Atomic{Int}(0)
    senders = map(1:nsenders) do i
        function sender_task()
            for j in 1:nrepeat
                put!(c, nsenders * (j - 1) + i)
            end
            if Threads.atomic_add!(nfinished, 1) + 1 == nsenders
                # Note: using `nfinished` instead of using a task that waits for
                # all senders, so that we can use a single `run_concurrently`
                # call (to make debugging easier).

                sleep(0.2)  # TODO: don't do this; support `close`
                # `sleep` here is just a "hint" but it's still bad...

                for _ in receivers
                    put!(c, sentinel)
                end
                @debug "`test_many_items`: All sent"
            end
        end
    end
    allsent = 1:nsenders*nrepeat

    run_concurrently([senders; receivers]; spawn)
    @debug "`test_many_items`: All sendres and receivers done"

    allreceived = reduce(vcat, received)
    nleft = 0
    ncleaned = 0
    while true
        ref = trytaking(c.items)()
        ref === nothing && break  # empty
        ncleaned += 1
        y = something(ref)[]
        y === nothing && continue  # taken
        push!(allreceived, something(y))
        nleft += 1
    end
    @debug "`test_many_items`: $nleft items not recieved; $ncleaned refs cleaned"

    @test length(allreceived) == length(allsent)
    @test sort!(allreceived) == allsent
end

end  # module


text:
__precompile__()

module WCS

export WCSTransform,
       pix_to_world, pix_to_world!,
       world_to_pix, world_to_pix!

import Base: convert, copy, deepcopy, getproperty, show, setproperty!, propertynames

using Base.Threads
const wcs_lock = SpinLock()

include(joinpath(@__DIR__, "..", "deps", "deps.jl"))

# -----------------------------------------------------------------------------
# constants defined in WCSLIB's wcshdr.h

# 'relax' flags for `from_header()`. See
# http://docs.astropy.org/en/stable/wcs/relax.html
# for descriptions.
const HDR_NONE     = 0x00000000
const HDR_ALL      = 0x000FFFFF
const HDR_REJECT   = 0x10000000
const HDR_STRICT   = 0x20000000

const HDR_CROTAia  = 0x00000001
const HDR_EPOCHa   = 0x00000002
const HDR_VELREFa  = 0x00000004
const HDR_CD00i00j = 0x00000008
const HDR_PC00i00j = 0x00000010
const HDR_PROJPn   = 0x00000020
const HDR_CD0i_0ja = 0x00000040
const HDR_PC0i_0ja = 0x00000080
const HDR_PV0i_0ma = 0x00000100
const HDR_PS0i_0ma = 0x00000200
const HDR_RADECSYS = 0x00000400
const HDR_VSOURCE  = 0x00000800
const HDR_DOBSn    = 0x00001000
const HDR_LONGKEY  = 0x00002000
const HDR_CNAMn    = 0x00004000
const HDR_AUXIMG   = 0x00008000
const HDR_ALLIMG   = 0x00010000

const HDR_IMGHEAD  = 0x00100000
const HDR_BIMGARR  = 0x00200000
const HDR_PIXLIST  = 0x00400000

# Constants defined in wcsmath.h
const UNDEFINED    = 987654321.0e99

# -----------------------------------------------------------------------------
# Utilities

function unsafe_store_vec!(p::Ptr{T}, v::AbstractVector{T}) where {T}
    for i = 1:length(v)
        unsafe_store!(p, v[i], i)
    end
end

# convert a string to a tuple of bytes
function convert_string(::Type{NTuple{N,UInt8}}, s::String) where {N}
    @assert length(s) < N
    @assert isascii(s)
    v = zeros(UInt8, N)  # intermediate array that we can fill
    copyto!(v, codeunits(s))
    (v...,)
end

# load an String from a tuple of bytes, truncating at first NULL
function convert_string(::Type{String}, v::NTuple{N,UInt8}) where {N}
    len = N

    # reduce length if we find a null
    i = 1
    for i = 1:N
        if v[i] == 0x00
            len = i - 1
            break
        end
    end
    s = Array{UInt8}(undef, len)
    copyto!(s, 1, v, 1, len)
    return String(s)  # wraps the array `s`
end

# load a String from a pointer, truncating at first NULL or maxlen
function convert_string(::Type{String}, ptr::Ptr{UInt8}, maxlen::Int)
    len = maxlen

    # reduce length if we find a null
    i = 1
    for i = 1:maxlen
        if unsafe_load(ptr, i) == 0x00
            len = i - 1
            break
        end
    end
    return unsafe_string(ptr, len)
end

macro check_type(k, v, t)
    :(typeof($(esc(v))) <: $t || error($(esc(k)), " must have type ", $t))
end

macro check_prop(k, prop, v, op, n)
    en = esc(n)
    :(($op)($prop($(esc(v))), $en) ||
      error($(esc(k)), " must have ", $prop, " ", $op, " ", $en))
end

macro same_size(a, b)
    :(size($(esc(a))) == size($(esc(b))) ||
      error($(string(a)), " must have the same dimensions as ", $(string(b))))
end

# -----------------------------------------------------------------------------
# error messaging

function get_error_message(i::Cint)
    msgptrs = cglobal((:wcs_errmsg, libwcs), Ptr{UInt8})
    msgptr = unsafe_load(msgptrs, i + 1)
    unsafe_string(msgptr)
end

function assert_ok(i::Cint)
    # Negative values aren't usually an error
    if i > 0
        error(get_error_message(i))
    end
end

function wcslib_version()
    vers = Array{Cint}(undef, 3)
    ccall((:wcslib_version, libwcs), Ptr{UInt8}, (Ptr{Cint},), vers)
    return VersionNumber(vers[1], vers[2], vers[3])
end

# -----------------------------------------------------------------------------
# WCSTransform

# mirror of `struct pvcard`
const PVCard = Tuple{Cint,Cint,Cdouble}  # i, m, value entries

# mirror of `struct pscard`
const PSCard = Tuple{Cint,Cint,NTuple{72,UInt8}}  # i, m, value entries

# mirror of `struct wcserr`
struct WCSErr
    status::Cint
    line_no::Cint
    _function::Ptr{UInt8}
    file::Ptr{UInt8}
    msg::NTuple{160,UInt8}
end

struct linprm
    flag::Cint
    naxis::Cint
    crpix::Ptr{Cdouble}
    pc::Ptr{Cdouble}
    cdelt::Ptr{Cdouble}
    dispre::Ptr{Cvoid}  # Ptr{disprm}
    disseq::Ptr{Cvoid}  # Ptr{disprm}
    piximg::Ptr{Cdouble}
    imgpix::Ptr{Cdouble}
    i_naxis::Cint
    unity::Cint
    affine::Cint
    simple::Cint
    err::Ptr{WCSErr}
    tmpcrd::Ptr{Cdouble}
    m_flag::Cint
    m_naxis::Cint
    m_crpix::Ptr{Cdouble}
    m_pc::Ptr{Cdouble}
    m_cdelt::Ptr{Cdouble}
    m_dispre::Ptr{Cvoid}  # Ptr{disprm}
    m_disseq::Ptr{Cvoid}  # Ptr{disprm}
end

struct prjprm
    flag::Cint
    code::NTuple{4,UInt8}
    r0::Cdouble
    pv::NTuple{30,Cdouble}
    phi0::Cdouble
    theta0::Cdouble
    bounds::Cint
    name::NTuple{40,UInt8}
    category::Cint
    pvrange::Cint
    simplezen::Cint
    equiareal::Cint
    conformal::Cint
    _global::Cint
    divergent::Cint
    x0::Cdouble
    y0::Cdouble
    err::Ptr{WCSErr}
    padding::Ptr{Cvoid}
    w::NTuple{10,Cdouble}
    m::Cint
    n::Cint
    prjx2s::Ptr{Cvoid}
    prjs2x::Ptr{Cvoid}
end

struct celprm
    flag::Cint
    offset::Cint
    phi0::Cdouble
    theta0::Cdouble
    ref::NTuple{4,Cdouble}
    prj::prjprm
    euler::NTuple{5,Cdouble}
    latpreq::Cint
    isolat::Cint
    err::Ptr{WCSErr}
    padding::Ptr{Cvoid}
end

struct spcprm
    flag::Cint
    _type::NTuple{8,UInt8}
    code::NTuple{4,UInt8}
    crval::Cdouble
    restfrq::Cdouble
    restwav::Cdouble
    pv::NTuple{7,Cdouble}
    w::NTuple{6,Cdouble}
    isGrism::Cint
    padding1::Cint
    err::Ptr{WCSErr}
    padding2::Ptr{Cvoid}
    spxX2P::Ptr{Cvoid}
    spxP2S::Ptr{Cvoid}
    spxS2P::Ptr{Cvoid}
    spxP2X::Ptr{Cvoid}
end

# mirror of `wcsprm` struct in wcslib
"""
    WCSTransform(naxis; kwds...)

Construct a WCS transformation with the given number of axes `naxis`.
Keyword arguments can be passed to set various attributes of the
transform. Specifying keyword arguments is equivalent to setting them
after construction:

```
julia> wcs = WCSTransform(2; crpix=[1000., 1000.])
```

is equilvalent to:

```
julia> wcs = WCSTransform(2)

julia> wcs.crpix = [1000., 1000.]
```

# Properties
Below is the entire list of public properties for a `WCSTransform`

| Keyword | Type | Description |
|--------:|:-----------------------------------------|:--------------------------------------------------------------------------------|
| naxis | `Int` | Number of dimensions |
| crval | `Vector{Float}[naxis]` | coordinate value at reference point |
| crpix | `Vector{Float}[naxis]` | array location of the reference point in pixels |
| cdelt | `Vector{Float}[naxis]` | coordinate increment at reference point |
| crder | `Vector{Float}[naxis]` | random error in coordinate |
| csyer | `Vector{Float}[naxis]` | systematic error in coordinate |
| ctype | `Vector{String}[naxis]` | axis type (8 characters) |
| crota | `Vector{Float}[naxis]` | rotation from stated coordinate type |
| cunit | `Vector{String}[naxis]` | units of axes |
| cunit | `Vector{String}[naxis]` | names of axes |
| pc | `Matrix{Float}[naxis, naxis]` | linear transformation matrix |
| cd | `Matrix{Float}[naxis, naxis]` | linear transformation matrix (with scale) |
| equinox | `Float` | the equinox associated with dynamical equatorial or ecliptic coordinate systems |
| latpole | `Float` | the native latitude of the celestial pole |
| lonpole | `Float` | the native longitude of the celestial pole |
| mjdavg | `Float` | Modified Julian Date corresponding to `DATE-AVG` |
| mjdobs | `Float` | Modified Julian Date corresponding to `DATE-OBS` |
| restfrq | `Float` | rest frequency (Hz) |
| restwav | `Float` | rest wavelength (m) |
| velangl | `Float` | velocity angle |
| velosys | `Float` | relative radial velocity |
| zsource | `Float` | the redshift of the source |
| colnum | `Int` | column of FITS binary table associated with this WCS |
| dateavg | `String` | representative mid-point of the date of observation |
| dateobs | `String` | start of the date of observation |
| radesys | `String` | the equatorial or ecliptic coordinate system type |
| specsys | `String` | spectral reference frame (standard of rest) |
| ssysobs | `String` | spectral reference frame |
| ssyssrc | `String` | spectral reference frame for redshift |
| wcsname | `String` | name of this coordinate representation |
| obsgeo | `Vector{Float}[3]` or `Vector{Float}[6]` | location of the observer in a standard terrestrial reference frame |
| alt | `String` | character code for alternate coordinate descriptions |
"""
mutable struct WCSTransform
    flag::Cint
    naxis::Cint
    crpix::Ptr{Cdouble}
    pc::Ptr{Cdouble}
    cdelt::Ptr{Cdouble}
    crval::Ptr{Cdouble}
    cunit::Ptr{Cvoid}
    ctype::Ptr{Cvoid}
    lonpole::Cdouble
    latpole::Cdouble
    restfrq::Cdouble
    restwav::Cdouble
    npv::Cint
    npvmax::Cint
    pv::Ptr{PVCard}
    nps::Cint
    npsmax::Cint
    ps::Ptr{PSCard}
    cd::Ptr{Cdouble}
    crota::Ptr{Cdouble}
    altlin::Cint
    velref::Cint
    alt::NTuple{4,UInt8}
    colnum::Cint
    colax::Ptr{Cint}
    cname::Ptr{Cvoid}
    crder::Ptr{Cdouble}
    csyer::Ptr{Cdouble}
    czphs::Ptr{Cdouble}
    cperi::Ptr{Cdouble}
    wcsname::NTuple{72,UInt8}
    timesys::NTuple{72,UInt8}
    trefpos::NTuple{72,UInt8}
    trefdir::NTuple{72,UInt8}
    plephem::NTuple{72,UInt8}
    timeunit::NTuple{72,UInt8}
    dateref::NTuple{72,UInt8}
    mjdref::NTuple{2,Cdouble}
    timeoffs::Cdouble
    dateobs::NTuple{72,UInt8}
    datebeg::NTuple{72,UInt8}
    dateavg::NTuple{72,UInt8}
    dateend::NTuple{72,UInt8}
    mjdobs::Cdouble
    mjdbeg::Cdouble
    mjdavg::Cdouble
    mjdend::Cdouble
    jepoch::Cdouble
    bepoch::Cdouble
    tstart::Cdouble
    tstop::Cdouble
    xposure::Cdouble
    telapse::Cdouble
    timsyer::Cdouble
    timrder::Cdouble
    timedel::Cdouble
    timepixr::Cdouble
    obsgeo::NTuple{6,Cdouble}
    obsorbit::NTuple{72,UInt8}
    radesys::NTuple{72,UInt8}
    equinox::Cdouble
    specsys::NTuple{72,UInt8}
    ssysobs::NTuple{72,UInt8}
    velosys::Cdouble
    zsource::Cdouble
    ssyssrc::NTuple{72,UInt8}
    velangl::Cdouble
    ntab::Cint
    nwtb::Cint
    tab::Ptr{Cvoid}  # Ptr{tabprm}
    wtb::Ptr{Cvoid}  # Ptr{wtbarr}
    lngtyp::NTuple{8,UInt8}
    lattyp::NTuple{8,UInt8}
    lng::Cint
    lat::Cint
    spec::Cint
    cubeface::Cint
    types::Ptr{Cint}
    lin::linprm
    cel::celprm
    spc::spcprm
    err::Ptr{WCSErr}
    m_flag::Cint
    m_naxis::Cint
    m_crpix::Ptr{Cdouble}
    m_pc::Ptr{Cdouble}
    m_cdelt::Ptr{Cdouble}
    m_crval::Ptr{Cdouble}
    m_cunit::Ptr{Cvoid}
    m_ctype::Ptr{Cvoid}
    m_pv::Ptr{PVCard}
    m_ps::Ptr{PSCard}
    m_cd::Ptr{Cdouble}
    m_crota::Ptr{Cdouble}
    m_colax::Ptr{Cint}
    m_cname::Ptr{Cvoid}
    m_crder::Ptr{Cdouble}
    m_csyer::Ptr{Cdouble}
    m_czphs::Ptr{Cdouble}
    m_cperi::Ptr{Cdouble}
    m_tab::Ptr{Cvoid}  # Ptr{tabprm}
    m_wtb::Ptr{Cvoid}  # Ptr{wtbarr}

    function WCSTransform(naxis::Integer; kwds...)
        w = new(-1)
        status = ccall((:wcsini, libwcs), Cint,
                       (Cint, Cint, Ref{WCSTransform}),
                       1, naxis, w)
        assert_ok(status)
        finalizer(free!, w)
        for (k, v) in kwds
            setproperty!(w, k, v)
        end

        # wcsset() modifies the WCSTransform, so call it here so it doesn't
        # get called in wcss2p and wcsp2s. Otherwise, calling wcss2p / wcsp2s
        # on the same WCSTransform object from multiple threads would create
        # a race condition.
        status = ccall((:wcsset, libwcs), Cint, (Ref{WCSTransform},), w)
        assert_ok(status)

        return w
    end
end

function free!(w::WCSTransform)
    status = ccall((:wcsfree, libwcs), Cint, (Ref{WCSTransform},), w)
    assert_ok(status)
end

function propertynames(::WCSTransform, private::Bool = false)
    if private
        return fieldnames(WCSTransform)
    else
        return (:naxis, :cdelt, :crder, :crota, :crpix, :crval, :csyer,
                :cname, :ctype, :cunit, :cd, :pc, :equinox, :latpole, :lonpole,
                :mjdavg, :mjdobs, :restfrq, :restwav, :velangl, :velosys, :zsource, :colnum,
                :dateavg, :dateobs, :radesys, :specsys, :ssysobs, :ssyssrc, :wcsname,
                :obsgeo, :alt)
    end
end

"""
    obsfix(ctrl::Integer, wcs::WCSTransform)

Complete the `obsgeo` field `wcs` of observatory coordinates.  That is, if only the (x,y,z)
Cartesian coordinate triplet or the (l,b,h) geodetic coordinate triplet are set, then it
derives the other triplet from it. If both triplets are set, then it checks for consistency
at the level of 1 metre.

# Parameters

* `ctrl`: flag that controls behaviour if one triplet is defined and the other is only
  partially defined:
    * 0: Reset only the undefined elements of an incomplete coordinate triplet.
    * 1: Reset all elements of an incomplete triplet.
    * 2: Don't make any changes, check for consistency only. Returns an error if either of
      the two triplets is incomplete.
* `wcs`: Coordinate transformation parameters. Its `obsgeo` field may be changed.

# Returns

* -1: No change required (not an error).
*  0: Success.
*  1: Null wcsprm pointer passed.
*  5: Invalid parameter value.
"""
function obsfix(ctrl::Integer, wcs::WCSTransform)
    status = ccall((:obsfix, libwcs), Cint, (Cint, Ref{WCSTransform}),
                   ctrl, wcs)
    assert_ok(status)
end

# -----------------------------------------------------------------------------
# getting attributes of a WCSTransform
# These return newly-allocated memory (not views of the WCSTransform).

function getproperty(wcs::WCSTransform, k::Symbol)
    @assert getfield(wcs, :flag) != -1
    setfield!(wcs, :flag, Cint(0))
    naxis = getfield(wcs, :naxis)

    if k == :naxis
        v = getfield(wcs, k)

    # double[naxis]
    elseif k in (:cdelt, :crder, :crota, :crpix, :crval, :csyer)
        v = Array{Float64}(undef, naxis)
        unsafe_copyto!(pointer(v), getfield(wcs, k), naxis)

    # char[72,naxis]
    elseif k in (:cname, :ctype, :cunit)
        p = convert(Ptr{UInt8}, getfield(wcs, k))
        v = Array{String}(undef, naxis)
        for i = 1:naxis
            pi = p + 72 * (i - 1)  # Ptr{UInt8} to the i-th entry.
            v[i] = convert_string(String, pi, 72)
        end

    # PVCard[]
    elseif k === :pv
        error("pv getter not yet implemented")

    # PSCard[]
    elseif k === :ps
        error("ps getter not yet implemented")

    # double[naxis,naxis]
    elseif k in (:cd, :pc)
        v = Array{Cdouble}(undef, naxis, naxis)
        unsafe_copyto!(pointer(v), getfield(wcs, k), naxis * naxis)

    # double
    elseif k in (:equinox, :latpole, :lonpole, :mjdavg, :mjdobs,
                 :restfrq, :restwav, :velangl, :velosys, :zsource)
        v = getfield(wcs, k)

    # int
    elseif k === :colnum
        v = Int(getfield(wcs, k))

    # char[72]
    elseif k in (:dateavg, :dateobs, :radesys, :specsys, :ssysobs, :ssyssrc,
                 :wcsname)
        v = convert_string(String, getfield(wcs, k))

    # double[6]
    elseif k === :obsgeo
        v = getfield(wcs, k)  # Tuple{Cdouble, Cdouble, Cdouble}

    # char[4], but only uses first
    elseif k === :alt
        v = Char(getfield(wcs, :alt)[1])

    else
        error("no property \"$k\" available")
    end

    return v
end

# -----------------------------------------------------------------------------
# modifying a WCSTransform

function setproperty!(wcs::WCSTransform, k::Symbol, v)
    @assert getfield(wcs, :flag) != -1
    setfield!(wcs, :flag, Cint(0))
    naxis = getfield(wcs, :naxis)

    # double[naxis]
    if k in (:cdelt, :crder, :crota, :crpix, :crval, :csyer)
        @check_type k v Vector{Float64}
        @check_prop k length v (==) naxis
        unsafe_copyto!(getfield(wcs, k), pointer(v), naxis)

    # char[72,naxis]
    elseif k in (:cname, :ctype, :cunit)
        @check_type k v Vector{String}
        @check_prop k length v (==) naxis
        for el in v; @assert isascii(el); end

        p = convert(Ptr{UInt8}, getfield(wcs, k))
        for i in 1:naxis
            pi = p + 72 * (i - 1)  # Ptr{UInt8} to i-th entry.
            n = length(v[i])
            @assert n < 72
            unsafe_copyto!(pi, pointer(v[i]), n)
            unsafe_store!(pi, 0x00, n + 1)  # trailing null
        end

    # PVCard[]
    elseif k === :pv
        @check_prop k length v (<=) getfield(wcs, :npvmax)
        npv = length(v)
        for i in 1:npv  # TODO: what if a convert fails halfway through?
            unsafe_store!(getfield(wcs, :pv), convert(PVCard, v[i]), i)
        end
        setfield!(wcs, :npv, Cint(npv))

    # PSCard[]
    elseif k === :ps
        @check_prop k length v (<=) wcs.npsmax
        nps = length(v)
        for i in 1:nps
            unsafe_store!(wcs.ps, convert(PSCard, v[i]), i)
        end
        wcs.nps = nps

    # double[naxis,naxis]
    elseif k in (:cd, :pc)
        @check_type k v Matrix{Float64}
        @check_prop k size v (==) (naxis, naxis)
        unsafe_store_vec!(getfield(wcs, k), vec(v'))

    # double
    elseif k in (:equinox, :latpole, :lonpole, :mjdavg, :mjdobs,
                 :restfrq, :restwav, :velangl, :velosys, :zsource)
        setfield!(wcs, k, convert(Float64, v))

    # int
    elseif k === :colnum
        setfield!(wcs, k, convert(Cint, v))

    # char[72]
    elseif k in (:dateavg, :dateobs, :radesys, :specsys, :ssysobs, :ssyssrc,
                 :wcsname)
        setfield!(wcs, k, convert_string(NTuple{72,UInt8}, v))

    # double[6]
    elseif k === :obsgeo
        @check_type k v Vector{Float64}
        if length(v) == 3
            # For compatibility with previous versions of WCSLIB which required a 3-element
            # vector.
            setfield!(wcs, :obsgeo, (v[1], v[2], v[3], UNDEFINED, UNDEFINED, UNDEFINED))
            obsfix(0, wcs)
        else
            @check_prop k length v (==) 6
            setfield!(wcs, :obsgeo, (v[1], v[2], v[3], v[4], v[5], v[6]))
            obsfix(2, wcs)
        end

    # char[4], but only uses first
    elseif k === :alt
        @check_type k v Char
        (('A' <= v <= 'Z') || v == ' ') || error("alt must be 'A'-'Z' or ' '")
        setfield!(wcs, :alt, (UInt8(v), 0x00, 0x00, 0x00))

    else
        error("no property \"$k\" available")
    end
end

function deepcopy(wcs::WCSTransform)
    dst = WCSTransform(wcs.naxis)
    p = convert(Ptr{Cint}, C_NULL)
    ccall((:wcssub, libwcs), Cint,
          (Cint, Ref{WCSTransform}, Ptr{Cint}, Ptr{Cint}, Ref{WCSTransform}),
          0, wcs, p, p, dst)
    return dst
end

# The default Julia `copy` would be a gotcha, as it would leave two
# WCSTransforms linked by their underlying C-allocated arrays;
# only allow deepcopy.
copy(w::WCSTransform) = deepcopy(w)

# TODO: more info here.
function show(io::IO, wcs::WCSTransform)
    print(io, "WCSTransform(naxis=$(wcs.naxis))")
end

# -----------------------------------------------------------------------------
# transforms

"""
    pix_to_world(wcs, pixcoords)

Convert the array of pixel coordinates `pixcoords` to world coordinates
according to the WCSTransform `wcs`. `pixcoords` should be a 2-d array
where \"pixcoords[:, i]\" is the i-th set of coordinates, or a 1-d array
representing a single set of coordinates.

The return value is the same shape as `pixcoords`.
"""
pix_to_world(wcs::WCSTransform, pixcoords::VecOrMat{Float64}) =
    pix_to_world!(wcs, pixcoords, similar(pixcoords))

"""
    pix_to_world!(wcs, pixcoords, worldcoords[; stat=, imcoords=, phi=, theta=])

Convert the array of pixel coordinates `pixcoords` to world coordinates
according to the WCSTransform `wcs`, storing the result in the
`worldcoords` and `stat` arrays. `pixcoords` should be a 2-d array where
\"pixcoords[:, i]\" is the i-th set of coordinates, or a 1-d array
representing a single set of coordinates. `worldcoords` must be
the same size and type as `pixcoords`.

If given, the arrays `stat`, `imcoords`, `phi`, `theta` will be used
to store intermediate results. Their sizes and types must all match
`pixcoords`, except for `stat` which should be the same size but of type
Cint (typically Int32).
"""
function pix_to_world!(wcs::WCSTransform, pixcoords::VecOrMat{Float64},
                       worldcoords::VecOrMat{Float64};
                       stat = similar(pixcoords, Cint),
                       imcoords = similar(pixcoords),
                       phi = similar(pixcoords),
                       theta = similar(pixcoords))
    nelem = size(pixcoords, 1)
    ncoords = size(pixcoords, 2)
    if nelem < wcs.naxis
        error("size(pixcoords, 1) must be greater than or equal to naxis")
    end
    @same_size worldcoords pixcoords
    @same_size imcoords pixcoords
    @same_size phi pixcoords
    @same_size theta pixcoords
    @same_size stat pixcoords
    ccall((:wcsp2s, libwcs), Cint,
          (Ref{WCSTransform}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble},
           Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}),
          wcs, ncoords, nelem, pixcoords, imcoords, phi, theta, worldcoords,
          stat)
    return worldcoords
end

"""
    world_to_pix(wcs, worldcoords)

Convert the array of world coordinates `worldcoords` to pixel coordinates
according to the WCSTransform `wcs`. `worldcoords` is a 2-d array
where \"worldcoords[:, i]\" is the i-th set of coordinates, or a 1-d array
representing a single set of coordinates.

The return value is the same size as `worldcoords`.
"""
world_to_pix(wcs::WCSTransform, worldcoords::VecOrMat{Float64}) =
    world_to_pix!(wcs, worldcoords, similar(worldcoords))

"""
    world_to_pix!(wcs, worldcoords, pixcoords[; stat=, phi=, theta=, imcoords=])

Convert the array of pixel coordinates `worldcoords` to pixel coordinates
according to the WCSTransform `wcs`, storing the result in the
`pixcoords` array. `worldcoords` should be a 2-d array where
\"worldcoords[:, i]\" is the i-th set of coordinates, or a 1-d array
representing a single set of coordinates. `pixcoords` must be
the same size and type as `worldcoords`.

If given, the arrays `stat`, `imcoords`, `phi`, `theta` will be used
to store intermediate results. Their sizes and types must all match
`worldcoords`, except for `stat` which should be the same size but of type
Cint (typically Int32).
"""
function world_to_pix!(wcs::WCSTransform, worldcoords::VecOrMat{Float64},
                       pixcoords::VecOrMat{Float64};
                       stat = similar(pixcoords, Cint),
                       phi = similar(pixcoords),
                       theta = similar(pixcoords),
                       imcoords = similar(pixcoords))
    nelem = size(worldcoords, 1)
    ncoords = size(worldcoords, 2)
    if nelem < wcs.naxis
        error("size(worldcoords, 1) must be greater than or equal to naxis")
    end
    @same_size pixcoords worldcoords
    @same_size imcoords worldcoords
    @same_size phi worldcoords
    @same_size theta worldcoords
    @same_size stat worldcoords
    ccall((:wcss2p, libwcs), Cint,
          (Ref{WCSTransform}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble},
           Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint}),
          wcs, ncoords, nelem, worldcoords, phi, theta, imcoords, pixcoords,
          stat)
    return pixcoords
end

# -----------------------------------------------------------------------------
# WCSTransform <--> header

"""
    from_header(header[; relax=WCS.HDR_ALL, ctrl=0, ignore_rejected=false, table=false])

Parse the FITS image header in the String `header`, returning a
`Vector{WCSTransform}` giving all the transforms defined in the header.
The `relax` determines the treatment of non-standard keywords. The default is
to accept all known non-standard keywords. Use `relax=WCS.HDR_NONE` to
ignore all non-standard keywords. Use, e.g.,
`relax=(WCS.HDR_RADECSYS & WCS.HDR_CROTAia)` to only accept selected
non-standard keywords.
"""
function from_header(header::String; relax::Integer = HDR_ALL, ctrl::Integer = 0,
                     ignore_rejected::Bool = false, table::Bool = false)
    @assert ctrl >= 0  # < 0 modifies the header
    nkeyrec::Integer = div(length(header), 80)
    nreject = Ref{Cint}(0)
    nwcs = Ref{Cint}(0)
    wcsptr = Ref{Ptr{WCSTransform}}(0)
    keysel = 0
    status = Cint(0)

    # wcsbth & wcspih are not thread-safe; see
    # http://www.atnf.csiro.au/people/mcalabre/WCS/wcslib/threads.html
    lock(wcs_lock)
    if table
        colsel = convert(Ptr{Cint}, C_NULL)
        status = ccall((:wcsbth, libwcs), Cint,
                       (Ptr{UInt8}, Cint, Cint, Cint, Cint, Ptr{Cint},
                        Ref{Cint}, Ref{Cint}, Ref{Ptr{WCSTransform}}),
                       header, nkeyrec, relax, ctrl, keysel, colsel,
                       nreject, nwcs, wcsptr)
    else
        status = ccall((:wcspih, libwcs), Cint,
                       (Ptr{UInt8}, Cint, Cint, Cint, Ref{Cint}, Ref{Cint},
                        Ref{Ptr{WCSTransform}}),
                       header, nkeyrec, relax, ctrl, nreject, nwcs, wcsptr)
    end
    unlock(wcs_lock)
    assert_ok(status)
    p = wcsptr[]
    result = WCSTransform[unsafe_load(p, i) for i = 1:nwcs[]]

    # Free the array of pointers allocated by wcslib. (But not the
    # structs they point to, which we've `unsafe_load`ed into a
    # Julia-allocated array!)
    # TODO: this could call the wrong `free` function if wcslib is linked
    # against a different libc than Julia.
    Libc.free(p)

    # For each of the WCSTransforms, register a finalizer and finish
    # initialization of the struct by calling wcsset. This avoids race
    # conditions between threads using the same WCSTransform.
    lock(wcs_lock)
    for w in result
        finalizer(free!, w)
        status = ccall((:wcsset, libwcs), Cint, (Ref{WCSTransform},), w)
        assert_ok(status)
    end
    unlock(wcs_lock)

    if !ignore_rejected && nreject[] != 0
        error("$(nreject[]) WCS transformations were rejected; " *
              "use ignore_rejected = true keyword to ignore")
    end
    return result
end

"""
    to_header(wcs[; relax=WCS.HDR_NONE])

Encode the WCSTransform `wcs` as a FITS header string. The `relax` keyword
controls how non-standard extensions to the WCS standard are handled.
"""
function to_header(wcs::WCSTransform; relax::Integer = HDR_NONE)
    nkeyrec = Ref{Cint}(0)
    hdrptr = Ref{Ptr{UInt8}}(C_NULL)
    status = ccall((:wcshdo, libwcs), Cint,
                   (Cint, Ref{WCSTransform}, Ref{Cint}, Ref{Ptr{UInt8}}),
                   relax, wcs, nkeyrec, hdrptr)
    assert_ok(status)
    header = unsafe_string(hdrptr[])
    Libc.free(hdrptr[])
    return header
end

function Base.getindex(wcs::WCSTransform, k::Symbol)
    Base.depwarn("`wcs[:$k]` is deprecated, use `wcs.$k` instead.", :getindex)
    getproperty(wcs, k)
end
function Base.setindex!(wcs::WCSTransform, v, k::Symbol)
    Base.depwarn("`wcs[:$k] = $v` is deprecated, use `wcs.$k = $v` instead.", :getindex)
    setproperty!(wcs, k, v)
end

end  # module


text:
# This file was generated, do not modify it. # hide
using MLJ
import Statistics
using PrettyPrinting
using StableRNGs

MLJ.color_off() # hide
X, y = @load_iris;


text:
export constructtriad
export constructtorus
export constructsphere
export constructcylinder
export constructbox
export constructwhirl
export constructframe
export getplane

"""
    constructtriad(q [; length])

Construct a triad with the given configuration `q` and `length`.
"""
function constructtriad(q::Biquaternion; length::Float64 = 1.0)
    map(x -> x + gettranslation(q), rotate(ℝ³.([[0; 0; 0], [length; 0; 0],
                                                [0; 0; 0], [0; length; 0],
                                                [0; 0; 0], [0; 0; length]]),
                                           getrotation(q)))
end

"""
    constructtorus(q, r, R)

Construct a torus of revolution with the given configuration `q`, the smaller radius `r`
and the bigger radius `R`.
"""
function constructtorus(q::Biquaternion,
                        r::Real,
                        R::Real;
                        segments::Int = 36)
    array = Array{ℝ³,2}(undef, segments, segments)
    for i in 1:segments
        for j in 1:segments
            ϕ = i * 2pi / (segments - 1)
            θ = j * 2pi / (segments - 1)
            x₁ = (R + r * cos(ϕ)) * cos(θ)
            x₂ = (R + r * cos(ϕ)) * sin(θ)
            x₃ = r * sin(ϕ)
            array[i, j] = ℝ³(x₁, x₂, x₃)
        end
    end
    map(x -> x + gettranslation(q), rotate(array, getrotation(q)))
end

"""
    constructsphere(q, radius)

Construct a sphere with the given configuration `q` and `radius`.
"""
function constructsphere(q::Biquaternion,
                         radius::Real;
                         segments::Int = 36)
    array = Array{ℝ³,2}(undef, segments, segments)
    lspace = collect(range(float(-pi), stop = float(pi), length = segments))
    for i in 1:segments
        for j in 1:segments
            ϕ = lspace[i]
            θ = lspace[j] / 2
            array[i, j] = ℝ³(Cartesian(Geographic(radius, ϕ, θ)))
        end
    end
    map(x -> x + gettranslation(q), rotate(array, getrotation(q)))
end

"""
    constructcylinder(q, height, radius)

Construct a cylinder with the given configuration `q`, `height` and `radius`.
"""
function constructcylinder(q::Biquaternion,
                           height::Float64,
                           radius::Float64;
                           segments::Int = 36)
    array = Array{ℝ³,2}(undef, segments, segments)
    lspace = collect(range(float(-pi), stop = float(pi), length = segments))
    for i in 1:segments
        for j in 1:segments
            α = lspace[j]
            x₁ = radius * cos(α)
            x₂ = radius * sin(α)
            x₃ = (i - 1) * height / (segments - 1)
            array[i, j] = ℝ³(x₁, x₂, x₃)
        end
    end
    map(x -> x + gettranslation(q), rotate(array, getrotation(q)))
end

"""
    constructbox(p, q, s)

Construct a box with the given translation `p`, rotation `q` and scale `s`.
"""
function constructbox(p::ℝ³, q::Quaternion, s::ℝ³)
    x, y, z = vec(s) ./ 2
    points = [ℝ³(-x, -y, -z) ℝ³(x, -y, -z) ℝ³(-x, -y, -z) ℝ³(x, -y, -z);
              ℝ³(-x, -y, +z) ℝ³(x, -y, +z) ℝ³(-x, -y, +z) ℝ³(x, -y, +z);
              ℝ³(-x, +y, +z) ℝ³(x, +y, +z) ℝ³(-x, +y, +z) ℝ³(x, +y, +z);
              ℝ³(-x, +y, -z) ℝ³(x, +y, -z) ℝ³(-x, +y, -z) ℝ³(x, +y, -z)]
    map(x -> x + p, rotate(points, q))
end

"""
    constructwhirl(points,
                   s2tos3map,
                   s2tos2map;
                   [top, [bottom, [s3rotation, [config, [segments [, scale]]]]]])

Construct a whirl with the given `points` in the base space, , a map from the base space to
the total space `s2tos3map`, a conformal map that takes the base space to itself `s2tos2map`
, `top` and `bottom` U(1) actions, S³ rotation `s3rotation`, configuration `config`, the
number of `segments` and `scale`.
"""
function constructwhirl(points::Array{<:S²,1},
                        s2tos3map,
                        s2tos2map;
                        top::S¹ = U1(-pi),
                        bottom::S¹ = U1(pi),
                        s3rotation::S³ = Quaternion(1, 0, 0, 0),
                        config::Biquaternion = Biquaternion(ℝ³(0, 0, 0)),
                        segments::Int = 36,
                        scale::Real = 1.0)
    array = Array{ℝ³,2}(undef, segments, length(points))
    lspace = range(angle(top), stop = angle(bottom), length = segments)
    for (i, α) in enumerate(lspace)
        for (j, p) in enumerate(points)
            array[i, j] = compressedλmap(rotate(S¹action(s2tos3map(s2tos2map(p)), U1(α)),
                                                s3rotation)) * scale
        end
    end
    applyconfig(array, config)
end

"""
    constructframe(circle, s2tos3map, s2tos2map; [s3rotation, [config, [segments]]])

Construct a frame that is like an S³ cross-section under stereographic projection with the
given `circle` that determines the phase angle in the fiber space, a map from the base space
to the total space `s2tos3map`, a conformal map that takes the base space to itself
`s2tos2map`, S³ rotation `s3rotation`, configuration `config` and the number of `segments`.
"""
function constructframe(circle::S¹,
                        s2tos3map,
                        s2tos2map;
                        s3rotation::S³ = Quaternion(1, 0, 0, 0),
                        config::Biquaternion = Biquaternion(ℝ³(0, 0, 0)),
                        segments::Int = 36)
    array = Array{ℝ³,2}(undef, segments, segments)
    lspaceϕ = collect(range(float(-pi), stop = float(pi), length = segments))
    lspaceθ = collect(range(float(pi / 2), stop = float(-pi / 2), length = segments))
    for (i, θ) in enumerate(lspaceθ)
        for (j, ϕ) in enumerate(lspaceϕ)
            p = Geographic(1, ϕ, θ)
            array[i, j] = compressedλmap(rotate(S¹action(s2tos3map(s2tos2map(p)), circle),
                                                s3rotation))
        end
    end
    applyconfig(array, config)
end

getplane(p::Array{Float64}, q::Quaternion, s::Array{Float64}) = begin
    plane = Array{Float64,3}(undef, 2, 2, 3)
    plane[1, 1, :] = vec(rotate(ℝ³([0.0; -1.0; 1.0]), q)) .* s + p
    plane[1, 2, :] = vec(rotate(ℝ³([0.0; 1.0; 1.0]), q)) .* s + p
    plane[2, 1, :] = vec(rotate(ℝ³([0.0; -1.0; -1.0]), q)) .* s + p
    plane[2, 2, :] = vec(rotate(ℝ³([0.0; 1.0; -1.0]), q)) .* s + p
    plane
end


text:
"""
    global_move(mc::MC, m::Model, conf) -> accepted::Bool

Propose a global move for configuration `conf`.
Returns wether the global move has been accepted or not.
"""
global_move(mc::MC, m::Model, conf) = false


text:
# Import Turing and Distributions.
using Turing, Distributions

# Import RDatasets.
using RDatasets

# Import MCMCChains, Plots, and StatPlots for visualizations and diagnostics.
using MCMCChains, Plots, StatsPlots

# Functionality for splitting and normalizing the data.
using MLDataUtils: shuffleobs, splitobs, rescale!

# Functionality for evaluating the model predictions.
using Distances

# Set a seed for reproducibility.
using Random
Random.seed!(0)

# Hide the progress prompt while sampling.
Turing.setprogress!(false);

# Import the "Default" dataset.
data = RDatasets.dataset("datasets", "mtcars");

# Show the first six rows of the dataset.
first(data, 6)

size(data)

# Remove the model column.
select!(data, Not(:Model))

# Split our dataset 70%/30% into training/test sets.
trainset, testset = splitobs(shuffleobs(data), 0.7)

# Turing requires data in matrix form.
target = :MPG
train = Matrix(select(trainset, Not(target)))
test = Matrix(select(testset, Not(target)))
train_target = trainset[:, target]
test_target = testset[:, target]

# Standardize the features.
μ, σ = rescale!(train; obsdim=1)
rescale!(test, μ, σ; obsdim=1)

# Standardize the targets.
μtarget, σtarget = rescale!(train_target; obsdim=1)
rescale!(test_target, μtarget, σtarget; obsdim=1);

# Bayesian linear regression.
@model function linear_regression(x, y)
    # Set variance prior.
    σ₂ ~ truncated(Normal(0, 100), 0, Inf)

    # Set intercept prior.
    intercept ~ Normal(0, sqrt(3))

    # Set the priors on our coefficients.
    nfeatures = size(x, 2)
    coefficients ~ MvNormal(nfeatures, sqrt(10))

    # Calculate all the mu terms.
    mu = intercept .+ x * coefficients
    return y ~ MvNormal(mu, sqrt(σ₂))
end

model = linear_regression(train, train_target)
chain = sample(model, NUTS(0.65), 3_000);

plot(chain)

describe(chain)

# Import the GLM package.
using GLM

# Perform multiple regression OLS.
train_with_intercept = hcat(ones(size(train, 1)), train)
ols = lm(train_with_intercept, train_target)

# Compute predictions on the training data set
# and unstandardize them.
p = GLM.predict(ols)
train_prediction_ols = μtarget .+ σtarget .* p

# Compute predictions on the test data set
# and unstandardize them.
test_with_intercept = hcat(ones(size(test, 1)), test)
p = GLM.predict(ols, test_with_intercept)
test_prediction_ols = μtarget .+ σtarget .* p;

# Make a prediction given an input vector.
function prediction(chain, x)
    p = get_params(chain[200:end, :, :])
    targets = p.intercept' .+ x * reduce(hcat, p.coefficients)'
    return vec(mean(targets; dims=2))
end

# Calculate the predictions for the training and testing sets
# and unstandardize them.
p = prediction(chain, train)
train_prediction_bayes = μtarget .+ σtarget .* p
p = prediction(chain, test)
test_prediction_bayes = μtarget .+ σtarget .* p

# Show the predictions on the test data set.
DataFrame(; MPG=testset[!, target], Bayes=test_prediction_bayes, OLS=test_prediction_ols)

println(
    "Training set:",
    "\n\tBayes loss: ",
    msd(train_prediction_bayes, trainset[!, target]),
    "\n\tOLS loss: ",
    msd(train_prediction_ols, trainset[!, target]),
)

println(
    "Test set:",
    "\n\tBayes loss: ",
    msd(test_prediction_bayes, testset[!, target]),
    "\n\tOLS loss: ",
    msd(test_prediction_ols, testset[!, target]),
)

let
    bayes_train_loss = msd(train_prediction_bayes, trainset[!, target])
    bayes_test_loss = msd(test_prediction_bayes, testset[!, target])
    ols_train_loss = msd(train_prediction_ols, trainset[!, target])
    @assert bayes_train_loss < 5.5
    @assert bayes_test_loss < 15.5
    @assert isapprox(bayes_train_loss, ols_train_loss; atol=1)
end

if isdefined(Main, :TuringTutorials)
    Main.TuringTutorials.tutorial_footer(WEAVE_ARGS[:folder], WEAVE_ARGS[:file])
end

max_stars_repo_name: TuringLang/TuringTutorialsOutput

max_issues_repo_name: TuringLang/TuringTutorialsOutput

max_forks_repo_name: TuringLang/TuringTutorialsOutput


text:
include("physical_constants.jl")

"Calculate delta0 from Tc and a custom value of the BCS constant."
function get_delta0(tc, bcs=1.76)
    return bcs*kb_eV*tc
end

"Calculate Tc from delta0 and a custom value of the BCS constant."
function get_tc(delta0, bcs=1.76)
    return delta0/(bcs*kb_eV)
end

"Calculate the BCS constant given a Tc and a delta0."
function get_bcs(tc, delta0)
    return delta0/(kb_eV*tc)
end

"Calculate the BCS coherence length from the Fermi velocity and delta0."
function get_ksi0(vf, delta0)
    return vf*hbar_eV/(π*delta0)
end

"Calculate the Fermi velocity from the BCS coherence length and delta0."
function get_vf(ksi0, delta0)
    return ksi0*π*delta0/hbar_eV
end

"""Calculate some unitless quantities needed for the computation from
physical data.

Parameters
----------
tc : float
    Critical temperature in Kelvin.

vf : float
    Fermi velocity in m/s.

mfp : float
    Mean free path in m.

london0 : float
    London penetration depth at zero temperature in m.

bcs : float (optional)
    The constant that gives the zero-temperature superconducting energy gap
    delta0 according to the equation delta0 = bcs*kB*Tc, where kB is
    Boltzmann's constant and Tc is the superconducting critical temperature.
    functionault value is the Bardeen-Cooper-Schrieffer value of 1.76.

ksi0 : float (optional)
    The BCS coherence length in m. ksi0 = hbar*vf/(pi*delta0), where hbar is
    Planck's constant divided by 2*pi, vf is the Fermi velocity, and delta0
    is the zero-temprature superconducting energy gap.

delta0 : float (optional)
    The zero-temperature superconducting energy gap in electron volts.
    delta0 = bcs*kB*Tc, where kB is Boltzmann's constant and Tc is the
    superconducting critical temperature. functionault value is the
    Bardeen-Cooper-Schrieffer value of 1.76.

Keyword Arguments
-----------------
verbose : int
    Whether to print out some debugging information. Very much a work in
    progress. functionault is 0.

    * 0: Don't print out any debugging information
    * 1: Print out minimal debugging information
    * 2: Print out full output from quad routines

Returns
-------
output_dict : dict
    output_dict contains the following keys: bcs, vf, x0, x1. These are all
    necessary inputs to impedance.cmplx_impedance."""
function init_from_physical_data(tc, vf, london0, mfp, bcs=1.76, ksi0=nothing, delta0=nothing, verbose=false)

    if delta0 == nothing
        delta0 = bcs*kb_eV*tc
    end

    if ksi0 == nothing
        ksi0 = vf*hbar_eV/(π*delta0)
        ksi0_calc = nothing
    else
        ksi0_calc = vf*hbar_eV/(π*delta0)
    end

    #Reduced lengths
    xk = ksi0/london0
    xm = mfp/london0

    # if verbose > 1
    #     if ksi0_calc != nothing
    #         print "calculated ksi0 = " + str(ksi0_calc*1e9)+" nm"
    #         print "supplied ksi0 = " + str(ksi0*1e9)+" nm"
    #     else
    #         print "calculated ksi0 = " + str(ksi0*1e9)+" nm"
    #     end
    #
    #     print "calculated delta0 = " + str(delta0*1e6) + " ueV"
    #     print "xk = ksi0/london0 = " + str(xk)
    #     print "xm = mfp/london0 = " + str(xm)
    #     print "xk/xm = ksi0/mfp = " + str(xk/xm)
    # end

    output_dict = Dict(:bc=>bcs,
                    :vf=>vf,
                    :xk=>xk,
                    :xm=>xm)

    return output_dict
end


text:
"""
SACPlot.jl provides routines for plotting SAC traces.
"""
module SACPlot

using Dates
using Printf
import DSP

import Plots

using SAC
import SAC: SACtr

const SACArray = AbstractArray{SACtr}

export
    plot1,
    p1,
    plot2,
    p2,
    plotpm,
    ppm,
    plotrs,
    plotrs!,
    plotspec,
    prs,
    prs!

const TIME_PICKS = [:a, :t0, :t1, :t2, :t3, :t4, :t5, :t6, :t7, :t8, :t9, :f]
const NAME_PICKS = Symbol.("k" .* String.(SACPlot.TIME_PICKS))

"Maximum number of samples to show by default"
const sacplot_qdp_thresh = 100_000
"Default plotting parameters for traces"
const sacplot_p1_defaults = Dict{Symbol,Any}()
"Default plotting parameters for record section"
const sacplot_prs_defaults = Dict{Symbol,Any}()

function __init__()
    props = Dict(:legend => false,
                 :grid => false,
                 :display => false,
                 :framestyle => :box)
    for (k, v) in props
        sacplot_p1_defaults[k] = v
    end
    props = Dict(:legend => false,
                 :grid => false,
                 :display => false,
                 :framestyle => :box)
    for (k, v) in props
        sacplot_prs_defaults[k] = v
    end
end

"""
    plot1(s::Array{SACtr}; kwargs...) -> ::Plots.Plot
    plot1(s::SACtr; kwargs...) -> ::Plots.Plot

Create a plot of the SAC trace(s) `s` and return a `Plots.Plot` object.

## Keyword arguments

| Keyword   | Example values | Description |
|:----------|:---------------|:------------|
|`label`    |`[:kcmpnm, :user0]`| A single or array of symbols listing the header values to show.|
|`line`     |`(:black, 2)`   | An argument to Plots defining the line property to use.  Arrays will apply in turn to each trace.|
|`over`     |`true`          | If `true`, plot on top of the active Plots plot object|
|`picks`    |`(["SKS"], [1610])`| Add arrival time picks to plot from a tuple of two arrays, one giving the phase name and the second the arrival time.|
|`qdp`      |`false`         | Plot every single sample if `false`.  (Default produces 'quick-and-dirty plot'.)|
|`relative` |`true`, `:user0`| Plot all times relative to beginning of trace if `true` or a header value if a `Symbol`|
|`xlim`     |`10:12`         | A range, tuple or array giving the limits in time to plot.  `NaN` as a limit uses the limits of the data.|
|`ylim`     |`(NaN, 2e-9)`   | As for `xlim`, but controlling the y-axis.  `:all` will force all traces to have the same scale.|
|`xlabel`   |`"Time / s"`    | Set the independent axis label|
|`ylabel`   |`"Amplitude / nm"`| Set the dependent axis label|
"""
function plot1(a::SACArray;
               label=:default, line=(:black,), over=false, picks=([],[]), qdp=true,
               relative=false, title="", xlabel="Time / s", xlim=[NaN, NaN],
               ylabel=nothing, ylim=[NaN, NaN], plots_kwargs...)
    # Check arguments
    check_lims(xlim, "SACPlot.plot1")
    typeof(ylim) <: Union{String, Symbol} || check_lims(ylim, "SACPlot.plot1")
    # Make sure any single label is in an array
    typeof(label) == Symbol && (label = [label])
    for l in label
        if l != :default
            l in [:kzdate, :kztime] && continue # Special labels
            l in fieldnames(a[1]) || error("Label header '$l' is not a valid header variable")
        end
    end
    # Make plots
    n = length(a)
    b, e, depmin, depmax = lims(a, relative)
    if ! isnan(xlim[1]); b = xlim[1]; end
    if ! isnan(xlim[end]); e = xlim[end]; end
    if ! (typeof(ylim) <: Union{String,Symbol})
        if ! isnan(ylim[1]); depmin = ylim[1]; end
        if ! isnan(ylim[end]); depmax = ylim[end]; end
    end
    # Downsample plot
    iskip = qdp_skip(a, qdp)
    # Turn off x labels for all but bottom trace
    p = over ? Plots.current() : Plots.plot(layout=(n,1), grid=false, legend=false,
        frame_style=:box; plots_kwargs...)
    for i = 1:n
        # Which samples to plot
        inds = qdp_inds(a[i], iskip)
        # Traces
        offset = lims_offset(a[i], relative)
        Plots.plot!(p[i], SAC.time(a[i])[inds] .- offset, a[i].t[inds],
            xlim=(b, e), xticks=(i==n ? true : false), l=line)
        if typeof(ylim) <: Union{String, Symbol}
            if Symbol(lowercase(string(ylim))) == :all
                Plots.ylims!(p[i], (depmin, depmax))
            end
        elseif ! all([isnan(ylim[1]), isnan(ylim[end])])
            Plots.ylims!(p[i], (depmin, depmax))
        end
        y1, y2 = Plots.ylims(p[i])
        # Add picks
        for (tp, kp) in zip(TIME_PICKS, NAME_PICKS)
            t, k = getfield(a[i], tp), strip(getfield(a[i], kp))
            if t != SAC.sac_rnull
                t -= offset
                Plots.plot!(p[i], [t, t], [y1, y2], l=:blue)
                if tp in [:a, :f] && k == SAC.sac_cnull
                    k = uppercase("$tp")
                end
                k != SAC.sac_cnull &&
                    Plots.plot!(p[i],
                        ann=(t, y1+0.05*(y2-y1), Plots.text(k, 8, :bottom, :left, :blue)))
            end
        end
        # Add custom picks
        phase_names, times = picks
        length(times) == length(phase_names) ||
            error("`picks` must be a tuple of phase name and time, with both arrays the same length")
        if length(times) > 0
            ilabel = 0
            for (k, t) in zip(phase_names, times)
                if b <= t <= e
                    Plots.plot!(p[i], [t, t], [y1, y2], l=:red)
                    Plots.plot!(p[i], ann=(t, y1+(1-0.05*ilabel)*(y2-y1), Plots.text(k, 8, :top, :left, :red)))
                    ilabel = (ilabel + 1)%4
                end
            end
        end
        # Add text annotation
        if label != :nothing
            if :default in label
                if all(Bool[getfield(a[i], f) != SAC.sac_inull for f in [:nzyear, :nzhour, :nzmin, :nzsec, :nzmsec]])
                    date = Date(Date(0) + Dates.Year(a[i].nzyear) + Dates.Day(a[i].nzjday - 1))
                    day = string(Dates.dayofmonth(date))
                    month = Dates.monthname(date)[1:3]
                else
                    day = month = "?"
                end
                label_string = strip(a[i].kevnm)*"\n"*strip(a[i].kstnm)*"."*
                               strip(a[i].knetwk)*"."*strip(a[i].kcmpnm) * "\n" *
                         @sprintf("%04d.%03dT%02d:%02d:%02d:%03d (%s %s)", a[i].nzyear,
                                  a[i].nzjday, a[i].nzhour, a[i].nzmin, a[i].nzsec,
                                  a[i].nzmsec, day, month)
            else
                label_string = ""
                for l in label
                    label_string *= string(l) * ": " * strip(string(getfield(a[i], l))) * "\n"
                end
            end
            Plots.plot!(p[i], ann=(b+0.97*(e-b), y1+0.95*(y2-y1),
                Plots.text(label_string, 8, :top, :right, :black)), annotationguide=:below)
        end
    end
    title != nothing && Plots.title!(p[1], title)
    xlabel != nothing && Plots.xlabel!(p[end], xlabel)
    ylabel != nothing && Plots.ylabel!(p[end÷2], ylabel)
    p
end

# Single-trace version
plot1(s::SACtr; kwargs...) = plot1([s]; kwargs...)
# Iterable which contains SAC traces
plot1(args...; kwargs...) = plot1(SACtr[args...]; kwargs...)

p1 = plot1

"""
    plot2(::Array{SACtr}; relative=false)

Plot all traces in array of SAC traces `a` on the same plot
"""
function plot2(a::SACArray; relative=false, legend=false)
    p = Plots.plot(legend=legend, show=false)
    offset = relative ? a[:b] : zeros(length(a))
    for i = 1:length(a)
        Plots.plot!(p, SAC.time(a[i]).-offset[i], a[i].t)
    end
    Plots.plot!(p, xlim=lims(a, relative)[1:2], show=true)
    p
end

# Single-trace version
plot2(s::SACtr) = plot2([s])

p2 = plot2

"""
    plotpm(::Array{SACtr}; xlim=[NaN, NaN], geog=false, kwargs...) -> ::Plots.Plot

Plot the particle motion for a pair of orthogonal components, within
the time window `xlim[1]` to `xlim[2]` if provided.

If `geog` is `true`, then plot the particle motion as relative to north,
with north up the page, and mark on the component directions.

Pass any further plotting commands to `Plots` with `kwargs...`.
"""
function plotpm(a::SACArray; xlim=[NaN, NaN], geog=false, kwargs...)
    angle_tol = 0.1
    length(a) == 2 || error("plotpm: Can only plot two components")
    angle = SAC.angle_difference(a[1][:cmpaz], a[2][:cmpaz])
    mod(angle - 90., 180.) <= angle_tol ||
        error("plotpm: Components must be orthogonal")
    # Set so that t1 is clockwise of t2 (e.g., t1 is E and t2 is N)
    if angle - 90. <= angle_tol
        t1, t2 = a[1], a[2]
    else
        t1, t2 = a[2], a[1]
    end
    amax = sqrt(maximum(a[1].t.^2 + a[2].t.^2))
    p = Plots.plot(xlim=(-amax,amax), ylim=(-amax,amax), aspect_ratio=:equal,
        legend=false, kwargs...)
    if geog
        az1, az2 = deg2rad(t1[:cmpaz]), deg2rad(t2[:cmpaz])
        Plots.plot!(p, amax.*[sin(az1), 0, sin(az2)], amax.*[cos(az1), 0, cos(az2)],
            line=(:darkgray, 2, :dash))
        for t in (t1, t2)
            rotation = mod(180-t[:cmpaz], 180) - 90
            font = Plots.Font("sans-serif", 8, :hcenter, :bottom,
                rotation, Plots.RGB(0.0,0.0,0.0))
            Plots.annotate!(p, 0.5amax*sind(t[:cmpaz]), 0.5amax*cosd(t[:cmpaz]),
                Plots.PlotText(t[:kcmpnm], font))
        end
        amp = sqrt.(t1.t.^2 + t2.t.^2)
        phi = atan2.(t1.t, t2.t) # Angle from t2->t1 (e.g., N->E)
        Plots.plot!(p, amp.*sin.(phi + az2), amp.*cos.(phi + az2),
            xlabel="East", ylabel="North")
    else
        Plots.plot!(p, t1.t, t2.t,
            xlabel=strip(t1.kcmpnm) * " (" * string(t1.cmpaz) * ")",
            ylabel=strip(t2.kcmpnm) * " (" * string(t2.cmpaz) * ")")
    end
    p
end

plotpm(s1::SACtr, s2::SACtr; kwargs...) = plotpm([s1, s2]; kwargs...)
ppm = plotpm

const _plotrs_first_run = Ref(true)

#TODO: Make plotting limits and automatic scaling more sensible when some traces
#      have very different amplitudes
"""
    plotrs!(p::Plots.Plot, s, align=:o; plots_kwargs...) -> ::Plots.Plot
    plotrs(s, align=:o; prs_kwargs..., plots_kwargs...) -> ::Plots.Plot

Plot a record section of the traces in `s`, aligned on the time given by
`align`.  This may be a header name, given as a Symbol, or an array of times.

By default, the y-axis is distance (`:gcarc`), but this can by any header name
passed as a symbol to the `y` keywords argument, or an arbitrary array of values
(for instance, to plot equally spaced in some order).

Other plotting options are controlled by keywords arguments `prs_kwargs`:

## Keyword arguments

|Name   |Type          |Description|
|:------|:-------------|:----------|
|`dw`   |Range or array|Set *distance window* (or other y-axis variable) for plotting|
|`fill` |`Tuple`       |Set colour fill.  Pass a 2-tuple of (positive_colour, negative_colour) or (pos, neg, level) to fill from/to a certain level between -1 and 1|
|`label`|`Symbol`      |Label traces with header value or array of values|
|`qdp`  |`Bool`        |If `true`, reduce the number of points plotted for speed. (Default is automatic)|
|`reverse`|`Bool`      |If `true`, reverse the direction of the y-axis.  (Default for `y=:gcarc`)|
|`tw`   |Range or array|Set *time window* for plotting|
|`scale`|Real          |Scaling factor for traces (can be negative to reverse polarity)|
|`line` |`Any`         |Argument passed to Plots specifying style for lines.  Use `false` to not plot lines.|
|`y`    |`Symbol` or array|Header value or array of values to use for y-axis|

Keyword arguments to be passed to Plots are those appended after any of the above;
e.g., to plot a record section aligned on the `A` header, and set the size of the
plot window using the Plots `size` argument, you can do:

```
julia> prs(sample(:array) |> cut(:a, -30, :a, 30), :a, scale=0.5, size=(500,1000))
```
"""
function plotrs!(p::Union{Plots.Plot,Plots.Subplot}, s::SACArray, align=0.;
        tw=[nothing, nothing], dw=[nothing, nothing], y=:gcarc, line=(:black,),
        scale::Real=1.0, reverse=nothing, fill=(nothing, nothing),
        qdp=sacplot_qdp_thresh, label=nothing,
        kwargs...)
    maxamp = maximum(abs, (s[:depmax]; s[:depmin]))
    y_shift = _y_shifts(s, y)
    d = _x_shifts(s, align)
    reverse = if reverse == nothing
        if y == :gcarc true else false end
    else
        reverse
    end
    scale = scale*abs(minimum(y_shift) - maximum(y_shift))/10
    reverse && (scale *= -1) # Set positive values to still plot 'up' the page if reversed axis
    # Determine a limit on total number of points to plot
    iskip = qdp_skip(s, qdp)
    # Traces
    for i in 1:length(s)
        inds = qdp_inds(s[i], iskip)
        t = (SAC.time(s[i]) .+ d[i])[inds]
        trace = y_shift[i] .+ s[i].t[inds].*scale/maxamp
        # FIXME: Filled wiggles not working with Plots yet
        if _plotrs_first_run[] && any([fill[1], fill[2]] .!= nothing)
            _plotrs_first_run[] = false
            @warn("Filled wiggles do not display perfectly yet")
        end

        fill_level = length(fill) >= 3 ? y_shift[i]+fill[3]*scale : y_shift[i]
        if fill[1] != nothing
            Plots.plot!(p, t, [tt  < fill_level ? tt : NaN for tt in trace],
                fillrange=fill_level, c=fill[1], line=Plots.stroke(0))
        end
        if fill[2] != nothing
            Plots.plot!(p, t, [tt >= fill_level ? tt : NaN for tt in trace],
                fillrange=fill_level, c=fill[2], line=Plots.stroke(0))
        end

        line != false && Plots.plot!(p, t, trace, l=line)
    end
    # x limits
    t1 = minimum(s[:b] + d)
    t2 = maximum(s[:e] + d)
    if !(tw[1] == tw[end] == nothing)
        tw[1] != nothing && (t1 = tw[1])
        tw[end] != nothing && (t2 = tw[end])
    end
    Plots.xlims!(p, t1, t2)
    # y limits
    if !(dw[1] == dw[end] == nothing)
        d1, d2 = Plots.ylims(p)
        dw[1] != nothing && (d1 = dw[1])
        dw[end] != nothing && (d2 = dw[end])
        Plots.ylims!(p, d1, d2)
    end
    reverse && Plots.yflip!(p)
    # Labels
    if label != nothing
        label_text = if label isa AbstractArray
            length(label) == length(s) ||
                throw(ArgumentError("`label` vector must be same length as number of traces"))
            string.(label)
        elseif label isa Symbol
            strip.(s[label])
        else
            throw(ArgumentError("Don't know how to plot labels for input: $label"))
        end
        t1, t2 = Plots.xlims(p)
        for i in 1:length(s)
            Plots.annotate!(p, t2, y_shift[i], Plots.text(label_text[i], 8, :left, :black),
                right_margin=10Plots.mm)
        end
    end
    Plots.plot!(p; kwargs...)
    p
end
plotrs(s::SACArray, args...; kwargs...) =
    plotrs!(Plots.plot(; sacplot_prs_defaults...), s, args...; kwargs...)
prs! = plotrs!
prs = plotrs
@doc (@doc plotrs!) plotrs

# Routines to calculate offsets in x and y for record section plots
_y_shifts(s::SACArray, y::Symbol) = s[y]
_y_shifts(s::SACArray, y::AbstractArray) =
    length(y) == length(s) ? y : error("Length of `y` not the same as number of traces")
_x_shifts(s::SACArray, t::Symbol) = -s[t]
_x_shifts(s::SACArray, t::AbstractArray) =
    length(s) == length(t) ? -t : error("Length of `align` not the same as number of traces")
_x_shifts(s::SACArray, t::Real) = -t*ones(length(s))

"""
    plotspec(s::SACtr; len=trace_length/20, overlap=4len/5, trace=true, trace_frac, f=identity, plots_kwargs...) -> ::Plots.Plot

Plot a spectrogram of the SAC trace `s` with windows of length `len` s which overlap by
`overlap` s.

By default, show the trace on a separate panel above the spectrogram; set `trace`
to `false` to only return the plot of the spectrogram.  `trace_frac` sets the proportion
of the plot to be taken up by the trace.  The power is trasformed by the function `f`
elementwise.  (For example, to show the natural log, pass `f=log`.)

`plots_kwargs...` are passed to the `Plots.heatmap` function which shows the spectrogram.
"""
function plotspec(s::SACtr;
        len=(s[:e] - s[:b])/20, overlap=0.8len,
        trace=true, trace_frac=0.3, f=identity, plot_kwargs...)
    n = round(Int, len/s[:delta])
    noverlap = clamp(round(Int, overlap/s[:delta]), 0, n-1)
    spec = DSP.spectrogram(s.t, n, noverlap, fs=1/s.delta)
    dspec = step(spec.time)/2
    colours = Plots.cgrad([:white,:blue,:green,:yellow,:orange,:red])
    pspec = Plots.heatmap(spec.time .+ s[:b], spec.freq, f.(spec.power),
        xlabel="Time / s", ylabel="Frequency / Hz", cbar=false, c=colours,
        framestyle=:box;
        plot_kwargs...)
    if trace
        p = Plots.plot(
            p1(cut(s, s[:b]+spec.time[1]-dspec, s[:b]+spec.time[end]+dspec), xlabel=""),
            pspec,
            layout=Plots.grid(2, 1, heights=(trace_frac, 1 - trace_frac)),
            xticks=nothing # Remove time from all x-axes
        )
        Plots.plot!(p[end], xticks=:auto) # Add time only to bottom subplot x-axis
    else
        pspec
    end
end

"""
    lims(::Array{SACtr}, relative=false) -> b, e, depmin, depmax

Get the minimum and maximum times, `b` and `e`, for an array of SAC traces.
If `relative` is true, these times are all relative to the beginning of the trace,
not the zero time.
"""
function lims(a::SACArray, relative=false)
    n = length(a)
    offset = lims_offset(a[1], relative)
    b = a[1].b - offset
    e = a[1].e - offset
    depmin = a[1].depmin
    depmax = a[1].depmax
    for i = 2:n
        offset = lims_offset(a[i], relative)
        if a[i].b - offset < b; b = a[i].b - offset; end
        if a[i].e - offset > e; e = a[i].e - offset; end
        if a[i].depmin < depmin; depmin = a[i].depmin; end
        if a[i].depmax > depmax; depmax = a[i].depmax; end
    end
    return b, e, depmin, depmax
end
lims(s::SACtr) = s.b, s.e, s.depmin, s.depmax
lims_offset(s::SACtr, relative::Bool) = relative ? s.b : zero(SAC.SACFloat)
lims_offset(s::SACtr, relative::Symbol) = getfield(s, relative)

"""
    check_lims(a, routine="SACPlot.check_lims")

Throw an error if `a` is not a length-2 array with the first value lower
or equal to the second, if both are not NaN.
"""
function check_lims(a, routine="SACPlot.check_lims")
    a[1] > a[2] && # NB: (NaN {>,<,==} NaN) == false
        error(routine * ": Lower plot limit value must be less than upper")
end

"""
    qdp_inds(s::SACtr, iskip) -> inds

Return `inds`, a range containing which indices to use for 'quick-and-dirty
plotting' when only each `iskip`th point is needed.
"""
qdp_inds(s::SACtr, iskip::Integer) = 1:min(iskip, s.npts-1):s.npts

"""
    qdp_skip(a::Union{SACtr,Array{SACtr}}, thresh::Integer=$(sacplot_qdp_thresh))

Return `iskip`, such that only the `iskip`th sample need be plotted when using
'quick-and-dirty plotting'.
"""
qdp_skip(a::Union{SACtr,SACArray}, thresh::Integer=sacplot_qdp_thresh) =
    max(1, sum(a[:npts])÷thresh)
qdp_skip(a::Union{SACtr,SACArray}, thresh::Nothing) = qdp_skip(a)
qdp_skip(a::Union{SACtr,SACArray}, thresh::Bool) = thresh ? qdp_skip(a) : 1

end # module SACPlot


text:
using ..Common
using Rectangle

import ..Common: CDDate
import Base: convert, promote_rule, Vector, String
import Rectangle: Rect

function convert(::Type{CDTextString}, xstr::CosXString)
    feff = [LATIN_F, LATIN_E, LATIN_F, LATIN_F]
    FEFF = [LATIN_UPPER_F, LATIN_UPPER_E, LATIN_UPPER_F, LATIN_UPPER_F]
    prefix = xstr.val[1:4]
    hasPrefix = (prefix == feff || prefix == FEFF)
    isUTF16   = hasPrefix || prefix[1:2] == UInt8[0x30, 0x30]
    data = get(xstr)
    buffer = data |> String |> hex2bytes
    if isUTF16
        len2 = div(length(buffer), 2)
        utf_16_arr = zeros(UInt16, hasPrefix ? len2-1 : len2)
        utf_16_data = reinterpret(UInt8, utf_16_arr)
        if (0x04030201 == ENDIAN_BOM)
            for i=1:len2
                (buffer[2i-1], buffer[2i]) = (buffer[2i], buffer[2i-1])
            end
        end
        hasPrefix ? copyto!(utf_16_data, 1, buffer, 3, 2len2-2) :
            copyto!(utf_16_data, 1, buffer, 1, 2len2)
        str = transcode(CDTextString, utf_16_arr)
    else
        str = CDTextString(PDFEncodingToUnicode(buffer))
    end
    return str
end

String(xstr::CosXString) = convert(String, xstr)

convert(::Type{Vector{UInt8}}, xstr::CosXString) =
    xstr |> get |> String |> hex2bytes

convert(::Type{Vector{UInt8}}, str::CosLiteralString) = str |> get

Vector{T}(str::CosString) where {T <: UInt8} = convert(Vector{UInt8}, str)

Vector{UInt8}(::CosNullType) = Vector{UInt8}(undef, 0)

convert(::Type{CDTextString}, lstr::IDD{CosLiteralString}) =
    CDTextString(PDFEncodingToUnicode(get(lstr)))

String(lstr::IDD{CosLiteralString}) = convert(String, lstr)

convert(::Type{T}, i::CosInt) where {T <: Number} = convert(T, get(i))

convert(::Type{T}, f::CosFloat) where {T <: Number} = convert(T, get(f))

convert(::Type{CosFloat}, i::CosInt) = i.val |> Float32 |> CosFloat

promote_rule(::Type{CosFloat}, ::Type{CosInt}) = CosFloat

function Rect(a::CosArray)
    v = promote(get(a, true)...)
    @assert length(v) == 4
    return Rect{typeof(v[1])}(v...)
end

convert(::Type{CDDate}, ls::T) where {T <: CosString} = CDDate(CDTextString(ls))

CDDate(str::CosString) = convert(CDDate, str)

convert(::Type{CDTextString}, name::CosName) =
    CDTextString(split(String(name.val), '_'; limit=2)[2])

String(name::CosName) = convert(String, name)


text:
method(r::OptimizationResults) = r.method
minimizer(r::OptimizationResults) = r.minimizer
minimum(r::OptimizationResults) = r.minimum
iterations(r::OptimizationResults) = r.iterations
iteration_limit_reached(r::OptimizationResults) = r.iteration_converged
trace(r::OptimizationResults) = length(r.trace) > 0 ? r.trace : error("No trace in optimization results. To get a trace, run optimize() with store_trace = true.")

function x_trace(r::UnivariateOptimizationResults)
	tr = trace(r)
	!haskey(tr[1].metadata, "minimizer") && error("Trace does not contain x. To get a trace of x, run optimize() with extended_trace = true")
    [ state.metadata["minimizer"] for state in tr ]
end
function x_lower_trace(r::UnivariateOptimizationResults)
	tr = trace(r)
	!haskey(tr[1].metadata, "x_lower") && error("Trace does not contain x. To get a trace of x, run optimize() with extended_trace = true")
    [ state.metadata["x_lower"] for state in tr ]
end
x_lower_trace(r::MultivariateOptimizationResults) = error("x_lower_trace is not implemented for $(method(r)).")
function x_upper_trace(r::UnivariateOptimizationResults)
	tr = trace(r)
	!haskey(tr[1].metadata, "x_upper") && error("Trace does not contain x. To get a trace of x, run optimize() with extended_trace = true")
    [ state.metadata["x_upper"] for state in tr ]
end
x_upper_trace(r::MultivariateOptimizationResults) = error("x_upper_trace is not implemented for $(method(r)).")

function x_trace(r::MultivariateOptimizationResults)
	tr = trace(r)
	!haskey(tr[1].metadata, "x") && error("Trace does not contain x. To get a trace of x, run optimize() with extended_trace = true")
    [ state.metadata["x"] for state in tr ]
end

f_trace(r::OptimizationResults) = [ state.value for state in trace(r) ]
g_norm_trace(r::OptimizationResults) = error("g_norm_trace is not implemented for $(method(r)).")
g_norm_trace(r::MultivariateOptimizationResults) = [ state.g_norm for state in trace(r) ]

f_calls(r::OptimizationResults) = r.f_calls

g_calls(r::OptimizationResults) = error("g_calls is not implemented for $(method(r)).")
g_calls(r::MultivariateOptimizationResults) = r.g_calls

converged(r::UnivariateOptimizationResults) = r.converged
converged(r::MultivariateOptimizationResults) = r.x_converged || r.f_converged || r.g_converged
x_converged(r::OptimizationResults) = error("x_converged is not implemented for $(method(r)).")
x_converged(r::MultivariateOptimizationResults) = r.x_converged
f_converged(r::OptimizationResults) = error("f_converged is not implemented for $(method(r)).")
f_converged(r::MultivariateOptimizationResults) = r.f_converged
g_converged(r::OptimizationResults) = error("g_converged is not implemented for $(method(r)).")
g_converged(r::MultivariateOptimizationResults) = r.g_converged

x_tol(r::OptimizationResults) = error("x_tol is not implemented for $(method(r)).")
x_tol(r::MultivariateOptimizationResults) = r.x_tol
f_tol(r::OptimizationResults) = error("f_tol is not implemented for $(method(r)).")
f_tol(r::MultivariateOptimizationResults) = r.f_tol
g_tol(r::OptimizationResults) = error("g_tol is not implemented for $(method(r)).")
g_tol(r::MultivariateOptimizationResults) = r.g_tol

initial_state(r::OptimizationResults) = error("initial_state is not implemented for $(method(r)).")
initial_state(r::MultivariateOptimizationResults) = r.initial_x

lower_bound(r::OptimizationResults) = error("lower_bound is not implemented for $(method(r)).")
lower_bound(r::UnivariateOptimizationResults) = r.initial_lower
upper_bound(r::OptimizationResults) = error("upper_bound is not implemented for $(method(r)).")
upper_bound(r::UnivariateOptimizationResults) = r.initial_upper

rel_tol(r::OptimizationResults) = error("rel_tol is not implemented for $(method(r)).")
rel_tol(r::UnivariateOptimizationResults) = r.rel_tol
abs_tol(r::OptimizationResults) = error("abs_tol is not implemented for $(method(r)).")
abs_tol(r::UnivariateOptimizationResults) = r.abs_tol


text:
deletesortedfirst!(a::Vector{T}, i::T) where T = deleteat!(a, searchsortedfirst(a, i))
funcmap(funcs::Vector{Function}, x::Vector) = map((f,xi) -> f(xi), funcs, x)
NamedTuple(d::Dict) = (;zip(keys(d), values(d))...)

function bisection(fn; a=0., b=1., ϵ=1e-05)
    μ = 0.
    while ((b - a) > ϵ)
        μ = 0.5*(a + b)
        if fn(μ) > 0
            b = μ
        else
            a = μ
        end
    end
    return 0.5*(a + b)
end