text:
# """
#     2-D vortex method.

#     # AUTHORSHIP
#         * Author    : Ryan Anderson
#         * Email     : rymanderson@gmail.com
#         * Created   : 2019
#         * License   : MIT License
# """
# module vm2d

# Import Modules
# import GeometricTools
# gt = GeometricTools
import LinearAlgebra
la = LinearAlgebra
import PyPlot
plt = PyPlot
include("vm2d_geometry.jl")
include("vm2d_particle.jl")
include("vm2d_plottools.jl")

# set freestream
function Uinf(x::Array{Float64,1})
    Uinf = [0.0,0.0,0.0]
    return Uinf
end

# # plot x-y view
# function plot_xy(particles::Array{Particle,1}; title::String="", save::Bool=false, name::String="default_sim", pause::Float64=-1.0)
#     # plot
#     plt.clf()
#     for particle in particles
#         plt.scatter(particle.x[1],particle.x[2],c=la.norm(particle.gamma))
#     end
#     plt.title(title)
#     plt.xlabel("x")
#     plt.ylabel("y")
#     plt.axis("equal")

#     # save frame
#     if save==true
#         plt.savefig("sim/"*name*".png")
#     end

#     # pause frame
#     if pause > 0.0
#         plt.pause(pause)
#     end
# end

# # plot y-z view
# function plot_yz(particles::Array{Particle,1}; title::String="", save::Bool=false, name::String="default_sim", pause::Float64=-1.0)
#     # plot
#     plt.clf()
#     for particle in particles
#         plt.scatter(particle.x[2],particle.x[3],c=la.norm(particle.gamma))
#     end
#     plt.title(title)
#     plt.xlabel("y")
#     plt.ylabel("z")
#     plt.axis("equal")

#     # save frame
#     if save==true
#         plt.savefig("sim/"*name*".png")
#     end

#     # pause frame
#     if pause > 0.0
#         plt.pause(pause)
#     end
# end

# Uinf(x) = [0.0,0.0,0.0] # default freestream velocity function
# # animate
# function animate(particles,timestart,timestep,numtimesteps; Uinf=Uinf)
#     currenttime = timestart
#     # set up plot
#     plt.plot()
#     plt.hot()
#     plt.axis("equal")
#     for particle in particles
#         plt.scatter(particle.x[1],particle.x[2],c=particle.gamma[3])
#     end
#     plt.title("time: "*string(round(currenttime; digits=3)))
#     # animate
#     for i in range(1,length=numtimesteps)
#         ## advance timestep
#         advance(particles,timestep,Uinf)
#         currenttime += timestep
#         ## update plot
#         plot_xy(particles; title="time: "*string(round(currenttime; digits=3)))
#     end
# end

# OPTION ONE: CIRCLE
# create circle geometry
# center = [0.0,0.0]
# radius = 0.1
# thetastart = 0.0
# thetaend = 2*pi
# refinement = 9
# points = circlexy(center,radius,refinement,thetastart,thetaend)

# # set vectorial vortex strengths
# gamma = [0,0,1.0]
# gammas = Array{Float64,1}[]
# for point in points
#     push!(gammas,gamma)
# end

# # set initial particle velocities
# vs = Array{Float64,1}[]
# for point in points
#     v = Uinf(point)
#     push!(vs,v)
# end

# particles = place(points,gammas,vs)

# # set simulation variables
# name = "vortexrings"
# global currenttime
# currenttime = 0.0
# timestep = 0.01
# numtimesteps = 150
# END OPTION 1

# OPTION 2: 2D VORTEX RINGS
# # create vortex ring geometry
# radius1 = 0.1
# radius2 = 0.1
# spacing = 0.05
# points = [[-spacing,radius1,0.0],[-spacing,-radius1,0.0],[0.0,radius2,0.0],[0.0,-radius2,0.0]]

# # set vectorial vortex strengths
# gamma = 0.1
# gamma_upper = [0,0,gamma]
# gamma_lower = [0,0,-gamma]
# gammas = [gamma_upper,gamma_lower,gamma_upper,gamma_lower]

# # set initial velocities
# v0 = 0.0
# vs = [[v0,0.0,0.0],[v0,0.0,0.0],[0.0,0.0,0.0],[0.0,0.0,0.0]]

# particles = place(points,gammas,vs)

# # set simulation variables
# name = "vortexrings"
# global currenttime
# currenttime = 0.0
# timestep = 0.01
# numtimesteps = 150
# END OPTION 2

# # OPTION 3: 3D VORTEX RINGS
# # create vortex ring geometry
# ## required arguments
# radius1 = 0.1
# radius2 = 0.1
# spacing = 0.05
# refinement = 9

# ## build geometry
# thetastart = 0.0
# thetaend = 2*pi

# center = [0.0,0.0,0.0]
# radius = radius1
# ring1 = circlexy(center,radius,refinement,thetastart,thetaend)

# center = [0.0,0.0,spacing]
# radius = radius2
# ring1 = circlexy(center,radius,refinement,thetastart,thetaend)

# # set vectorial vortex strengths
# gamma = 0.1
# gamma_upper = [0,0,gamma]
# gamma_lower = [0,0,-gamma]
# gammas = [gamma_upper,gamma_lower,gamma_upper,gamma_lower]

# # set initial velocities
# v0 = 0.0
# vs = [[v0,0.0,0.0],[v0,0.0,0.0],[0.0,0.0,0.0],[0.0,0.0,0.0]]

# particles = place(points,gammas,vs)

# # set simulation variables
# name = "vortexrings"
# global currenttime
# currenttime = 0.0
# timestep = 0.01
# numtimesteps = 150
# # END OPTION 3

# # plot
# ## OPTION 1: toggle 2D x-y plot
# plt.plot()
# plt.hot()
# plt.axis("equal")
# for particle in particles
#     plt.scatter(particle.x[1],particle.x[2],c=particle.gamma[3])
# end
# plt.title("time: "*string(round(currenttime; digits=3)))
# ## save frame
# plt.savefig("sim/"*name*"_"*string(currenttime)*".png")
# ## run simulation
# for i in range(1,length=numtimesteps)
#     ## advance timestep
#     advance(particles,timestep,Uinf)
#     currenttime += timestep
#     ## update plot
#     plt.clf()
#     for particle in particles
#         plt.scatter(particle.x[1],particle.x[2],c=particle.gamma[3])
#     end
#     plt.title("time: "*string(round(currenttime; digits=3)))
#     plt.axis("equal")
#     ## save frame
#     plt.savefig("sim/"*name*"_"*string(i)*".png")
#     println(particles)
# end


text:
# GFF3 Writer
# ===========

"""
    GFF3.Writer(output::IO)

Create a data writer of the GFF3 file format.

# Arguments:
* `output`: data sink
"""
mutable struct Writer{T<:IO} <: BioCore.IO.AbstractReader
    stream::T
end

function BioCore.IO.stream(writer::Writer)
    return writer.stream
end

function Base.write(writer::Writer, record::Record)
    return write(writer.stream, record, '\n')
end

ext: jl


text:
function FieldPoleCorrection(arg0::RealFieldElement, arg1::RealFieldElement)
    return FieldPoleCorrection((RealFieldElement, RealFieldElement), arg0, arg1)
end

function get_xp(obj::FieldPoleCorrection)
    return jcall(obj, "getXp", RealFieldElement, ())
end

function get_yp(obj::FieldPoleCorrection)
    return jcall(obj, "getYp", RealFieldElement, ())
end


text:
@testset "Matrix algebras" begin

  @testset "Construction" begin

    M1 = matrix(FlintQQ, [ 1 1; 0 1])
    M2 = matrix(FlintQQ, [ 1 0; 1 1])
    A = matrix_algebra(FlintQQ, [ M1, M2 ])
    @test dim(A) == 4
    @test basis_matrix(A) == identity_matrix(QQ, 4)

    QG = group_algebra(FlintQQ, small_group(2, 1))
    M1 = matrix(QG, [ 1 1; 0 1])
    M2 = matrix(QG, [ 1 0; 1 1])
    A = matrix_algebra(FlintQQ, QG, [ M1, M2 ])
    @test dim(A) == 4
  end

  @testset "Radical" begin

    mats = [[0 0 0 0 0 0 0 0 0 0;
             1 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0],
            [1 0 0 0 0 0 0 0 0 0;
             0 1 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0],
            [0 0 0 0 0 0 0 0 0 0;
             0 0 1 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0],
            [0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 1 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0],
            [0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 1 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0],
            [0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             1 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0],
            [0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 1 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0],
            [0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 1 0 0 0 0 0 0 0;
             0 0 0 1 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0],
            [0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 1 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0],
            [0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 1 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0],
            [0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             1 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0],
            [0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 1 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0],
            [0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 1 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0],
            [0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 1 0 0 0 0 0;
             0 0 0 0 0 1 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0],
            [0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 1 0 0 0 0 0;
             0 0 0 0 0 1 0 0 0 0;
             0 0 0 0 0 0 1 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0],
            [0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 1 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0],
            [0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             1 0 0 0 0 0 0 0 0 0],
            [0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 1 0 0 0 0 0 0 0],
            [0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 1 0 0 0 0 0],
            [0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 1 0 0],
            [0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 1 0 0;
             0 0 0 0 0 0 0 0 1 0],
            [0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 0 0 0;
             0 0 0 0 0 0 0 1 0 0;
             0 0 0 0 0 0 0 0 1 0;
             0 0 0 0 0 0 0 0 0 1]]

    mats2 = map(v -> matrix(FlintQQ, v), mats)
    A = matrix_algebra(FlintQQ, mats2)
    J = radical(A)
    B, = quo(A, J)
    @test dim(B) == 4
  end
end


text:
export sr1

function sr1(f::Function,J::Function,x::Vector;B=speye(length(x)),Delta=2.0, maxIter=20,atol=1e-8,out::Int=0,storeInterm::Bool=false,r=1e-8,
trt::Function=(f,fc,df,d2f,x,st,Delta)->trtest(f,fc,df,d2f,x,st,Delta))

    his = zeros(maxIter,5)
    X   = (storeInterm) ? zeros(length(x),maxIter) : []

    i = 1; flag = -1; alph=0.0
    while i<=maxIter
        fc = f(x)
        df  = J(x)

        his[i,1:2] = [fc norm(df)]
        if storeInterm; X[:,i] = x; end;
        if norm(df)<atol
            his  = his[1:i,:]
            flag = 0
            break
        end

        # get cauchy point
        curv = dot(df,B*df)
        if curv <= 0
            alph = Delta/norm(df)
        else
            alph = min(Delta/norm(df), norm(df)^2/curv)
        end
        st     = -alph*df
        xcp    = x + st
        flagnt = 0

        if norm(xcp-x) < Delta
            # get Newton point
            try
               sn = - B \ df
               xn = x + sn
               # check angle condition
               angl = dot(xn-xcp,xcp-x)
               if (angl>0) && (norm(xcp-x) < norm(sn)) && norm(sn) <= Delta
                  flagnt=1
                  st = sn # take newton point
               elseif (angl>0) # build dogleg
                  flagnt = 2
                  step = Delta/10; cnt = 1;

                  while (norm(xcp + cnt*step*(xn-xcp) - x) < Delta) && (cnt < 100)
                       cnt+=1
                  end
                  xt = xcp + (cnt-1)*step*(xn-xcp)
                  st = xt-x
               end
            catch err
                if ~isa(err,Base.LinAlg.SingularException)
                    throw(err)
                end
            end
        end

        # check for acceptance
        xn,Deltan,trflag,ared = trt(f,fc,df,B,x,st,Delta)
        his[i,3:5] = [trflag,flagnt,Delta]
        if out>0
            @printf "iter=%04d\t|f|=%1.2e\t|df|=%1.2e\tTR=%d\tnt?=%d\tDelta=%1.2f\n" i his[i,1] his[i,2] his[i,3] his[i,4] his[i,5]
        end
        # update Hessian approximation
        dfnew = J(x+st)
        yk    = dfnew - df
        if abs(dot(st,yk-B*st)) >= r*norm(st)*norm(yk-B*st)
	       B +=  ((yk-B*st)*(yk-B*st)')/dot(st,yk-B*st)
        end
		x     = xn
		Delta = Deltan
		i+=1
    end
    i = min(maxIter,i)

    if out>=0
        if flag==-1
            println(@sprintf("sr1 iterated maxIter (=%d) times but reached only atol of %1.2e instead of tol=%1.2e",i,his[i,2],atol))
        elseif out>1
            println(@sprintf("sr1 achieved desired atol of %1.2e at iteration %d.",atol,i))
        end
    end

    if storeInterm; X = X[:,1:i]; end
    return x,flag,his,X,B
end


text:
type DiscountingBondEngine{Y <: YieldTermStructure} <: PricingEngine
  yts::Y

  # function call(::Type{DiscountingBondEngine})
  #   new{YieldTermStructure}()
  # end
  #
  # function call{Y}(::Type{DiscountingBondEngine}, yts::Y)
  #   new{Y}(yts)
  # end
end

DiscountingBondEngine() = DiscountingBondEngine{NullYieldTermStructure}(NullYieldTermStructure())

function _calculate!{B <: Bond}(pe::DiscountingBondEngine, bond::B)
  yts = pe.yts
  valuation_date = yts.referenceDate
  value = npv(bond.cashflows, yts, valuation_date, valuation_date)
  bond.settlementValue = value
  # if the valuation_date is the same as the bond settlement date, then we don't have to recalculate

  return bond
end

clone(pe::DiscountingBondEngine, ts::TermStructure) = DiscountingBondEngine(ts)


text:
"""
Physical constants

The definition and default values are:

- `K₀ = -273.15`: absolute zero (°C)
- `R = 8.314`: universal gas constant (``J\\ mol^{-1}\\ K^{-1}``).
- `Rd = 287.0586`: gas constant of dry air (``J\\ kg^{-1}\\ K^{-1}``).
- `Dₕ₀ = 21.5e-6`: molecular diffusivity for heat at base temperature, applied in the integrated form of the
Fick’s Law of diffusion (``m^2\\ s^{-1}``). See eq. 3.10 from Monteith and Unsworth (2013).
- `Cₚ = 1013.0`: Specific heat of air at constant pressure (``J\\ K^{-1}\\ kg^{-1}``), also
known as efficiency of impaction of particles. See Allen et al. (1998), or Monteith and
Unsworth (2013). NB: bigleaf R package uses 1004.834 intead.
- `ε = 0.622`: ratio of molecular weights of water vapor and air. See Monteith and
Unsworth (2013).
- `λ₀ = 2.501`: latent heat of vaporization for water at 0 degree (``J\\ kg^{-1}``).
- `σ = 5.670373e-08` [Stefan-Boltzmann constant](https://en.wikipedia.org/wiki/Stefan%E2%80%93Boltzmann_law)
in (``W\\ m^{-2}\\ K^{-4}``).
- `Gbₕ_to_Gbₕ₂ₒ = 1.075`: conversion coefficient from conductance to heat to conductance to water
vapor.
- `Gsc_to_Gsw = 1.57`: conversion coefficient from stomatal conductance to CO₂ to conductance to water
vapor.
- `Gbc_to_Gbₕ = 1.32`: conversion coefficient from boundary layer conductance to CO₂ to heat.
- `Mₕ₂ₒ = 18.0e-3` (kg mol-1): Molar mass for water.

# References

Allen, Richard G., Luis S. Pereira, Dirk Raes, et Martin %J Fao Smith. 1998.
« Crop evapotranspiration-Guidelines for computing crop water requirements-FAO Irrigation
and drainage paper 56 » 300 (9): D05109.

Monteith, John, et Mike Unsworth. 2013. Principles of environmental physics: plants,
animals, and the atmosphere. Academic Press.
"""
Base.@kwdef struct Constants{T}
    K₀::T = -273.15
    R::T = 8.314
    Rd::T = 287.0586
    Dₕ₀::T = 21.2e-6
    Cₚ::T = 1013.0
    ε::T = 0.622
    λ₀::T = 2.501e6
    σ::T = 5.670373e-08
    Gbₕ_to_Gbₕ₂ₒ::T = 1.075
    Gsc_to_Gsw::T = 1.57
    Gbc_to_Gbₕ::T = 1.32
    Mₕ₂ₒ::T = 18.0e-3
    J_to_umol::T = 4.57
end


text:
"""
"""
function pair_arrays(a::AbstractArray,b::AbstractArray)
  lazy_map(tuple,a,b)
end

"""
"""
function unpair_arrays(pair::AbstractArray{<:Tuple{<:Any,<:Any}})
  a = lazy_map(first,pair)
  b = lazy_map(second,pair)
  a, b
end

second(x) = x[2]

function unpair_arrays(pair::LazyArray{<:Fill{typeof(tuple)}})
  pair.args
end


text:
# MathOptInterface constraints

##################################################
## Support constraints
MOI.supports_constraint(::Optimizer, ::Type{MOI.VariableIndex}, ::Type{MOI.ZeroOne}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.VariableIndex}, ::Type{MOI.Integer}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.VariableIndex}, ::Type{<:SS}) = true

function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.ScalarAffineFunction{Float64}},
    ::Type{MOI.LessThan{Float64}},
)
    return true
end
function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.ScalarAffineFunction{Float64}},
    ::Type{MOI.GreaterThan{Float64}},
)
    return true
end
function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.ScalarAffineFunction{Float64}},
    ::Type{MOI.EqualTo{Float64}},
)
    return true
end
function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.ScalarAffineFunction{Float64}},
    ::Type{MOI.Interval{Float64}},
)
    return true
end
function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.ScalarQuadraticFunction{Float64}},
    ::Type{MOI.LessThan{Float64}},
)
    return true
end
function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.ScalarQuadraticFunction{Float64}},
    ::Type{MOI.GreaterThan{Float64}},
)
    return true
end
function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.ScalarQuadraticFunction{Float64}},
    ::Type{MOI.EqualTo{Float64}},
)
    return true
end
function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.ScalarQuadraticFunction{Float64}},
    ::Type{MOI.Interval{Float64}},
)
    return true
end
MOI.supports_constraint(::Optimizer, ::Type{<:SF}, ::Type{<:SS}) = true
MOI.supports_constraint(::Optimizer, ::Type{VAF}, ::Type{<:VLS}) = true
MOI.supports_constraint(::Optimizer, ::Type{VOV}, ::Type{<:VLS}) = true
MOI.supports_constraint(::Optimizer, ::Type{VAF}, ::Type{MOI.SecondOrderCone}) = true
MOI.supports_constraint(::Optimizer, ::Type{VOV}, ::Type{MOI.SecondOrderCone}) = true
MOI.supports_constraint(::Optimizer, ::Type{VOV}, ::Type{MOI.Complements}) = true
MOI.supports_constraint(::Optimizer, ::Type{VAF}, ::Type{MOI.Complements}) = true

##################################################
# TODO: clean getters
## Getters
MOI.get(model::Optimizer, ::MOI.NumberOfConstraints) = KN_get_number_cons(model.inner)
function MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{MOI.VariableIndex,MOI.ZeroOne})
    return model.number_zeroone_constraints
end
function MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{MOI.VariableIndex,MOI.Integer})
    return model.number_integer_constraints
end
function MOI.get(
    model::Optimizer,
    ::MOI.NumberOfConstraints{MOI.VariableIndex,S},
) where {S<:SS}
    return sum(
        typeof.(collect(keys(model.constraint_mapping))) .==
        MOI.ConstraintIndex{MOI.VariableIndex,S},
    )
end
# TODO: a bit hacky, but that should work for MOI Test.
function MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{VAF,MOI.Nonnegatives})
    return sum(
        typeof.(collect(keys(model.constraint_mapping))) .==
        MOI.ConstraintIndex{VAF,MOI.Nonnegatives},
    )
end
function MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{VAF,MOI.Nonpositives})
    return sum(
        typeof.(collect(keys(model.constraint_mapping))) .==
        MOI.ConstraintIndex{VAF,MOI.Nonpositives},
    )
end
function MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{VAF,MOI.Zeros})
    return sum(
        typeof.(collect(keys(model.constraint_mapping))) .==
        MOI.ConstraintIndex{VAF,MOI.Zeros},
    )
end
function MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{VAF,MOI.SecondOrderCone})
    return sum(
        typeof.(collect(keys(model.constraint_mapping))) .==
        MOI.ConstraintIndex{VAF,MOI.SecondOrderCone},
    )
end
function MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{VOV,T}) where {T<:VLS}
    return sum(
        typeof.(collect(keys(model.constraint_mapping))) .== MOI.ConstraintIndex{VOV,T},
    )
end
function MOI.get(
    model::Optimizer,
    ::MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64},S},
) where {S<:LS}
    return sum(
        typeof.(collect(keys(model.constraint_mapping))) .==
        MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},S},
    )
end
function MOI.get(
    model::Optimizer,
    ::MOI.NumberOfConstraints{MOI.ScalarQuadraticFunction{Float64},S},
) where {S<:LS}
    return sum(
        typeof.(collect(keys(model.constraint_mapping))) .==
        MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{Float64},S},
    )
end

##################################################
# Generic constraint definition
#--------------------------------------------------
# Bound constraint on variables.
function MOI.add_constraint(
    model::Optimizer,
    vi::MOI.VariableIndex,
    lt::MOI.LessThan{Float64},
)
    (model.number_solved >= 1) && throw(AddConstraintError())
    check_inbounds(model, vi)
    if isnan(lt.upper)
        error("Invalid upper bound value $(lt.upper).")
    end
    if has_upper_bound(model, vi)
        error("Upper bound on variable $vi already exists.")
    end
    if is_fixed(model, vi)
        error("Variable $vi is fixed. Cannot also set upper bound.")
    end
    ub = check_value(lt.upper)
    model.variable_info[vi.value].has_upper_bound = true
    # By construction, MOI's indexing is the same as KNITRO's indexing.
    KN_set_var_upbnd(model.inner, vi.value - 1, ub)
    ci = MOI.ConstraintIndex{MOI.VariableIndex,MOI.LessThan{Float64}}(vi.value)
    model.constraint_mapping[ci] = convert(Cint, vi.value)
    return ci
end

function MOI.add_constraint(
    model::Optimizer,
    vi::MOI.VariableIndex,
    gt::MOI.GreaterThan{Float64},
)
    (model.number_solved >= 1) && throw(AddConstraintError())
    check_inbounds(model, vi)
    if isnan(gt.lower)
        error("Invalid lower bound value $(gt.lower).")
    end
    if has_lower_bound(model, vi)
        error("Lower bound on variable $vi already exists.")
    end
    if is_fixed(model, vi)
        error("Variable $vi is fixed. Cannot also set lower bound.")
    end
    lb = check_value(gt.lower)
    model.variable_info[vi.value].has_lower_bound = true
    # We assume that MOI's indexing is the same as KNITRO's indexing.
    KN_set_var_lobnd(model.inner, vi.value - 1, lb)
    ci = MOI.ConstraintIndex{MOI.VariableIndex,MOI.GreaterThan{Float64}}(vi.value)
    model.constraint_mapping[ci] = convert(Cint, vi.value)
    return ci
end

function MOI.add_constraint(
    model::Optimizer,
    vi::MOI.VariableIndex,
    set::MOI.Interval{Float64},
)
    (model.number_solved >= 1) && throw(AddConstraintError())
    check_inbounds(model, vi)
    if isnan(set.lower) || isnan(set.upper)
        error("Invalid lower bound value $(set.lower).")
    end
    if has_lower_bound(model, vi) || has_upper_bound(model, vi)
        error("Bounds on variable $vi already exists.")
    end
    if is_fixed(model, vi)
        error("Variable $vi is fixed. Cannot also set lower bound.")
    end
    lb = check_value(set.lower)
    ub = check_value(set.upper)
    model.variable_info[vi.value].has_lower_bound = true
    model.variable_info[vi.value].has_upper_bound = true
    # We assume that MOI's indexing is the same as KNITRO's indexing.
    KN_set_var_lobnd(model.inner, vi.value - 1, lb)
    KN_set_var_upbnd(model.inner, vi.value - 1, ub)
    ci = MOI.ConstraintIndex{MOI.VariableIndex,MOI.Interval{Float64}}(vi.value)
    model.constraint_mapping[ci] = convert(Cint, vi.value)
    return ci
end

function MOI.add_constraint(
    model::Optimizer,
    vi::MOI.VariableIndex,
    eq::MOI.EqualTo{Float64},
)
    (model.number_solved >= 1) && throw(AddConstraintError())
    check_inbounds(model, vi)
    if isnan(eq.value)
        error("Invalid fixed value $(eq.value).")
    end
    if has_lower_bound(model, vi)
        error("Variable $vi has a lower bound. Cannot be fixed.")
    end
    if has_upper_bound(model, vi)
        error("Variable $vi has an upper bound. Cannot be fixed.")
    end
    if is_fixed(model, vi)
        error("Variable $vi is already fixed.")
    end
    eqv = check_value(eq.value)
    model.variable_info[vi.value].is_fixed = true
    # We assume that MOI's indexing is the same as KNITRO's indexing.
    KN_set_var_fxbnd(model.inner, vi.value - 1, eqv)
    ci = MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{Float64}}(vi.value)
    model.constraint_mapping[ci] = convert(Cint, vi.value)
    return ci
end

function MOI.add_constraint(
    model::Optimizer,
    func::MOI.ScalarAffineFunction{Float64},
    set::LS,
)
    (model.number_solved >= 1) && throw(AddConstraintError())
    check_inbounds(model, func)
    # Add a single constraint in KNITRO.
    num_cons = KN_add_con(model.inner)
    # Add bound to constraint.
    if isa(set, MOI.LessThan{Float64})
        val = check_value(set.upper)
        KN_set_con_upbnd(model.inner, num_cons, val - func.constant)
    elseif isa(set, MOI.GreaterThan{Float64})
        val = check_value(set.lower)
        KN_set_con_lobnd(model.inner, num_cons, val - func.constant)
    elseif isa(set, MOI.EqualTo{Float64})
        val = check_value(set.value)
        KN_set_con_eqbnd(model.inner, num_cons, val - func.constant)
    end
    # Parse structure of constraint.
    indexvars, coefs = canonical_linear_reduction(func)
    KN_add_con_linear_struct(model.inner, num_cons, indexvars, coefs)
    # Add constraint to index.
    ci = MOI.ConstraintIndex{typeof(func),typeof(set)}(num_cons)
    model.constraint_mapping[ci] = num_cons
    return ci
end

function MOI.add_constraint(
    model::Optimizer,
    func::MOI.ScalarAffineFunction{Float64},
    set::MOI.Interval{Float64},
)
    (model.number_solved >= 1) && throw(AddConstraintError())
    check_inbounds(model, func)
    # Add a single constraint in KNITRO.
    num_cons = KN_add_con(model.inner)
    # Add upper bound.
    lb = check_value(set.lower)
    ub = check_value(set.upper)
    KN_set_con_lobnd(model.inner, num_cons, lb - func.constant)
    KN_set_con_upbnd(model.inner, num_cons, ub - func.constant)
    # Parse structure of constraint.
    indexvars, coefs = canonical_linear_reduction(func)
    KN_add_con_linear_struct(model.inner, num_cons, indexvars, coefs)
    # Add constraints to index.
    consi = MOI.ConstraintIndex{typeof(func),typeof(set)}(num_cons)
    model.constraint_mapping[consi] = num_cons
    return consi
end

function MOI.add_constraint(
    model::Optimizer,
    func::MOI.ScalarQuadraticFunction{Float64},
    set::SS,
)
    (model.number_solved >= 1) && throw(AddConstraintError())
    check_inbounds(model, func)
    # We add a constraint in KNITRO.
    num_cons = KN_add_con(model.inner)
    # Add upper bound.
    if isa(set, MOI.LessThan{Float64})
        val = check_value(set.upper)
        KN_set_con_upbnd(model.inner, num_cons, val - func.constant)
    elseif isa(set, MOI.GreaterThan{Float64})
        val = check_value(set.lower)
        KN_set_con_lobnd(model.inner, num_cons, val - func.constant)
    elseif isa(set, MOI.EqualTo{Float64})
        val = check_value(set.value)
        KN_set_con_eqbnd(model.inner, num_cons, val - func.constant)
    elseif isa(set, MOI.Interval{Float64})
        lb = check_value(set.lower)
        ub = check_value(set.upper)
        KN_set_con_lobnd(model.inner, num_cons, lb - func.constant)
        KN_set_con_upbnd(model.inner, num_cons, ub - func.constant)
    end
    # Parse linear structure of constraint.
    indexvars, coefs = canonical_linear_reduction(func)
    KN_add_con_linear_struct(model.inner, num_cons, indexvars, coefs)
    # Parse quadratic term.
    qvar1, qvar2, qcoefs = canonical_quadratic_reduction(func)
    # Take care that Knitro is 0-indexed!
    qvar1 .= qvar1 .- 1
    qvar2 .= qvar2 .- 1

    KN_add_con_quadratic_struct(model.inner, num_cons, qvar1, qvar2, qcoefs)
    # Add constraints to index.
    ci = MOI.ConstraintIndex{typeof(func),typeof(set)}(num_cons)
    model.constraint_mapping[ci] = num_cons
    return ci
end

function MOI.add_constraint(
    model::Optimizer,
    func::MOI.VectorAffineFunction,
    set::MOI.AbstractVectorSet,
)
    (model.number_solved >= 1) && throw(AddConstraintError())
    # TODO: add check inbounds for VectorAffineFunction.
    previous_col_number = number_constraints(model)
    num_cols = length(func.constants)
    # Add constraints inside KNITRO.
    index_cons = KN_add_cons(model.inner, num_cols)

    # Parse vector affine expression.
    indexcols, indexvars, coefs = canonical_vector_affine_reduction(func)
    # Reformate indexcols with current numbering.
    indexcols .+= previous_col_number

    # Load inside KNITRO.
    KN_add_con_linear_struct(model.inner, indexcols, indexvars, coefs)
    if isa(set, MOI.Nonnegatives)
        KN_set_con_lobnds(model.inner, num_cols, index_cons, -func.constants)
    elseif isa(set, MOI.Nonpositives)
        KN_set_con_upbnds(model.inner, num_cols, index_cons, -func.constants)
    elseif isa(set, MOI.Zeros)
        KN_set_con_eqbnds(model.inner, num_cols, index_cons, -func.constants)
    else
        error("Invalid set $set for VectorAffineFunction constraint")
    end
    # Add constraints to index.
    ci = MOI.ConstraintIndex{typeof(func),typeof(set)}(index_cons[1])
    model.constraint_mapping[ci] = index_cons
    return ci
end

# Add second order cone constraint.
function MOI.add_constraint(
    model::Optimizer,
    func::MOI.VectorAffineFunction,
    set::MOI.SecondOrderCone,
)
    (model.number_solved >= 1) && throw(AddConstraintError())
    # Add constraints inside KNITRO.
    index_con = KN_add_con(model.inner)

    # Parse vector affine expression.
    indexcoords, indexvars, coefs = canonical_vector_affine_reduction(func)
    constants = func.constants
    # Distinct two parts of secondordercone.
    # First row corresponds to linear part of SOC.
    indlinear = indexcoords .== 0
    indcone = indexcoords .!= 0
    ncoords = length(constants) - 1
    @assert ncoords == set.dimension - 1

    # Load Second Order Conic constraint.
    ## i) linear part
    KN_set_con_upbnd(model.inner, index_con, constants[1])
    KN_add_con_linear_struct(
        model.inner,
        index_con,
        indexvars[indlinear],
        -coefs[indlinear],
    )

    ## ii) soc part
    index_var_cone = indexvars[indcone]
    nnz = length(index_var_cone)
    index_coord_cone = convert.(Cint, indexcoords[indcone] .- 1)
    coefs_cone = coefs[indcone]
    const_cone = constants[2:end]

    KN_add_con_L2norm(
        model.inner,
        index_con,
        ncoords,
        nnz,
        index_coord_cone,
        index_var_cone,
        coefs_cone,
        const_cone,
    )

    # Add constraints to index.
    ci = MOI.ConstraintIndex{typeof(func),typeof(set)}(index_con)
    model.constraint_mapping[ci] = indexvars
    return ci
end

function MOI.add_constraint(model::Optimizer, func::VOV, set::T) where {T<:VLS}
    (model.number_solved >= 1) && throw(AddConstraintError())
    indv = convert.(Cint, [v.value for v in func.variables] .- 1)
    bnd = zeros(Float64, length(indv))

    if isa(set, MOI.Zeros)
        KN_set_var_fxbnds(model.inner, indv, bnd)
    elseif isa(set, MOI.Nonnegatives)
        KN_set_var_lobnds(model.inner, indv, bnd)
    elseif isa(set, MOI.Nonpositives)
        KN_set_var_upbnds(model.inner, indv, bnd)
    end

    ncons = MOI.get(model, MOI.NumberOfConstraints{VOV,T}())
    ci = MOI.ConstraintIndex{VOV,T}(ncons)
    model.constraint_mapping[ci] = indv
    return ci
end

function MOI.add_constraint(
    model::Optimizer,
    func::MOI.VectorOfVariables,
    set::MOI.SecondOrderCone,
)
    (model.number_solved >= 1) && throw(AddConstraintError())
    # Add constraints inside KNITRO.
    index_con = KN_add_con(model.inner)
    indv = [v.value - 1 for v in func.variables]

    KN_set_con_upbnd(model.inner, index_con, 0.0)
    KN_add_con_linear_struct(model.inner, index_con, indv[1], -1.0)

    indexVars = convert.(Cint, indv[2:end])
    nnz = length(indexVars)
    indexCoords = Cint[i for i in 0:(nnz-1)]
    coefs = ones(Float64, nnz)
    constants = zeros(Float64, nnz)

    KN_add_con_L2norm(
        model.inner,
        index_con,
        nnz,
        nnz,
        indexCoords,
        indexVars,
        coefs,
        constants,
    )

    # Add constraints to index.
    ci = MOI.ConstraintIndex{typeof(func),typeof(set)}(index_con)
    model.constraint_mapping[ci] = convert.(Cint, indv)
    return ci
end

##################################################
## Binary & Integer constraints.

# Define integer and boolean constraints.
function MOI.add_constraint(model::Optimizer, vi::MOI.VariableIndex, ::MOI.ZeroOne)
    indv = vi.value - 1
    check_inbounds(model, vi)
    model.number_zeroone_constraints += 1
    # Made the bounds explicit for KNITRO and take care of
    # preexisting MOI's bounds.
    lobnd = 0.0
    upbnd = 1.0
    KN_set_var_lobnd(model.inner, indv, lobnd)
    KN_set_var_upbnd(model.inner, indv, upbnd)
    KN_set_var_type(model.inner, vi.value - 1, KN_VARTYPE_BINARY)

    return MOI.ConstraintIndex{MOI.VariableIndex,MOI.ZeroOne}(vi.value)
end

function MOI.add_constraint(model::Optimizer, vi::MOI.VariableIndex, ::MOI.Integer)
    (model.number_solved >= 1) && throw(AddConstraintError())
    check_inbounds(model, vi)
    model.number_integer_constraints += 1
    KN_set_var_type(model.inner, vi.value - 1, KN_VARTYPE_INTEGER)
    return MOI.ConstraintIndex{MOI.VariableIndex,MOI.Integer}(vi.value)
end

##################################################
## Constraint DualStart
function MOI.supports(
    ::Optimizer,
    ::MOI.ConstraintDualStart,
    ::MOI.ConstraintIndex{<:SF,<:SS},
)
    return true
end
function MOI.set(
    model::Optimizer,
    ::MOI.ConstraintDualStart,
    ci::MOI.ConstraintIndex{<:SF,<:SS},
    value::Union{Real,Nothing},
)
    if isa(value, Real)
        KN_set_con_dual_init_values(model.inner, ci.value, Cdouble(value))
    else
        # By default, initial value is set to 0.
        KN_set_con_dual_init_values(model.inner, ci.value, Cdouble(0.0))
    end
    return
end

function MOI.supports(
    ::Optimizer,
    ::MOI.ConstraintDualStart,
    ::MOI.ConstraintIndex{MOI.VariableIndex,<:LS},
)
    return true
end
function MOI.set(
    model::Optimizer,
    ::MOI.ConstraintDualStart,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,<:LS},
    value::Union{Real,Nothing},
)
    if isa(value, Real)
        KN_set_var_dual_init_values(model.inner, ci.value, Cdouble(value))
    else
        # By default, initial value is set to 0.
        KN_set_var_dual_init_values(model.inner, ci.value, Cdouble(0.0))
    end
    return
end

##################################################
## Constraint naming
# TODO: dry supports with macros.
MOI.supports(::Optimizer, ::MOI.ConstraintName, ::Type{MOI.ConstraintIndex}) = true


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# SourceLocation struct for file position

"""
    SourceLocation

Represents a position in a file at a certain line and column.

# Fields

- `file_name::String`: the input file name
- `line_num::Int`: line position
- `col_num::Int`: column position
"""
Base.@kwdef mutable struct SourceLocation
    file_name::String = ""
    line_num::Int = 1
    col_num::Int = 0
end

@doc """
    SourceLocation(file_name::String, line_num::Int, col_num::Int)

Constructor for a [`SourceLocation`](@ref) instance.
""" SourceLocation(::String, ::Int, ::Int)

@doc """
    SourceLocation(; file_name::String = "", line_num::Int = 1, col_num::Int = 0)

Constructor for a [`SourceLocation`](@ref) instance.
""" SourceLocation(; ::String, ::Int, ::Int)

function Base.show(io::IO, ::MIME"text/plain", loc::SourceLocation)
    print(io, abspath(loc.file_name), ":", loc.line_num, ":", loc.col_num)
end

function Base.print(io::IO, loc::SourceLocation)
    print(io, abspath(loc.file_name), ":", loc.line_num, ":", loc.col_num)
end

function Base.copy(loc::SourceLocation)
    SourceLocation(loc.file_name, loc.line_num, loc.col_num)
end


text:
using MOCNeutronTransport
@testset "Tree" begin
    root = Tree("root" )
    @test root.data === "root"
    @test root.parent === nothing
    @test root.children === nothing
    l1_1 = Tree("L1_1", root)
    @test l1_1.data === "L1_1"
    @test l1_1.parent === root
    @test root.children == [l1_1]
    l1_2 = Tree("L1_2", root)
    @test l1_2.data === "L1_2"
    @test l1_2.parent === root
    @test root.children == [l1_1, l1_2]
end

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
# N is the number of dimensions of the prediction
abstract Loss{N}
typealias UnivariateLoss Loss{0}
typealias MultivariateLoss Loss{1}

## Abs loss (for regression)
#
#   loss(p, y) = |p - y|
#
immutable AbsLoss <: UnivariateLoss
end

value{T<:BlasReal}(::AbsLoss, p::T, y::T) = abs(p - y)
deriv{T<:BlasReal}(::AbsLoss, p::T, y::T) = sign(p - y)
value_and_deriv{T<:BlasReal}(::AbsLoss, p::T, y::T) = (r = p - y; (abs(r), sign(r)))

## Squared loss (for regression)
#
#   loss(p, y) := (1/2) * (p - y)^2
#
immutable SqrLoss <: UnivariateLoss
end

value{T<:BlasReal}(::SqrLoss, p::T, y::T) = half(abs2(p - y))
deriv{T<:BlasReal}(::SqrLoss, p::T, y::T) = p - y
value_and_deriv{T<:BlasReal}(::SqrLoss, p::T, y::T) = (r = p - y; v = half(abs2(r)); (v, r))

## Quantile loss (for quantile regression, asymmetric version of Abs loss)
#
#   loss(p, y) := t * (p - y)         ... (p >= y)
#               = (1 - t) * (y - p)   ... (p < y)
#
immutable QuantileLoss <: UnivariateLoss
    t::Float64

    function QuantileLoss(t::Real)
        zero(t) < t < one(t) || error("t must be a real value in (0, 1).")
        new(convert(Float64, t))
    end
end

function value{T<:BlasReal}(loss::QuantileLoss, p::T, y::T)
    t = convert(T, loss.t)
    p >= y ? t * (p - y) : (one(T) - t) * (y - p)
end

function deriv{T<:BlasReal}(loss::QuantileLoss, p::T, y::T)
    t = convert(T, loss.t)
    p > y ? t :
    p < y ? t - one(T) : zero(T)
end

function value_and_deriv{T<:BlasReal}(loss::QuantileLoss, p::T, y::T)
    t = convert(T, loss.t)
    if p > y
        (t * (p - y), t)
    elseif p < y
        c = t - one(T)
        (c * (p - y), c)
    else
        (zero(T), zero(T))
    end
end

## Epsilon Insensitive loss (for support vector regression)
#
#   loss(p, y) := 0                   ... abs(y - p) <= eps
#               = abs(y - p) - eps    ... otherwise
#
immutable EpsilonInsLoss <: UnivariateLoss
    epsilon::Float64

    function EpsilonInsLoss(epsilon::Real)
        new(convert(Float64, epsilon))
    end
end

function value{T<:BlasReal}(loss::EpsilonInsLoss, p::T, y::T)
    eps = convert(T, loss.epsilon)
    a = abs(p - y)
    a > eps ? a - eps : zero(T)
end

function deriv{T<:BlasReal}(loss::EpsilonInsLoss, p::T, y::T)
    eps = convert(T, loss.epsilon)
    r = p - y
    abs(r) > eps ? sign(r) : zero(T)
end

function value_and_deriv{T<:BlasReal}(loss::EpsilonInsLoss, p::T, y::T)
    eps = convert(T, loss.epsilon)
    r = p - y
    a = abs(r)
    a > eps ? (a - eps, sign(r)) : (zero(T), zero(T))
end

## Huber loss (for regression, smoothed version of Abs loss)
#
#   loss(p, y) := (1/2) * (p - y)^2      ... (|p - y| <= h)
#                 h * |p - y| - h^2/2    ... otherwise
#
immutable HuberLoss <: UnivariateLoss
    h::Float64

    function HuberLoss(h::Real)
        h > zero(h) || error("h must be a positive value.")
        new(convert(Float64, h))
    end
end

function value{T<:BlasReal}(loss::HuberLoss, p::T, y::T)
    h = convert(T, loss.h)
    a = abs(p - y)
    a <= h ? half(a * a) : h * a - half(h * h)
end

function deriv{T<:BlasReal}(loss::HuberLoss, p::T, y::T)
    h = convert(T, loss.h)
    r = p - y
    r > h ? h : r < -h ? -h : r
end

function value_and_deriv{T<:BlasReal}(loss::HuberLoss, p::T, y::T)
    h = convert(T, loss.h)
    r = p - y
    r > h ? (h * r - half(h * h), h) :
    r < -h ? (-h * r - half(h * h), -h) :
    (half(r * r), r)
end

## Hinge loss (for L1-SVM)
#
#   loss(p, y) := max(1 - y * p, 0)
#
immutable HingeLoss <: UnivariateLoss
end

value{T<:BlasReal}(::HingeLoss, p::T, y::T) = nonneg(one(T) - y * p)
deriv{T<:BlasReal}(::HingeLoss, p::T, y::T) = y * p < one(T) ? -y : zero(T)

function value_and_deriv{T<:BlasReal}(::HingeLoss, p::T, y::T)
    yp = y * p
    yp >= one(T) ? (zero(T), zero(T)) : (one(T) - yp, -y)
end

## Squared Hinge loss (for L2-SVM)
#
#   loss(p, y) := max(1 - y * p, 0)^2
#
immutable SqrHingeLoss <: UnivariateLoss
end

function value{T<:BlasReal}(::SqrHingeLoss, p::T, y::T)
    yp = y * p
    yp >= one(T) ? zero(T) : abs2(nonneg(one(T) - yp))
end

deriv{T<:BlasReal}(::SqrHingeLoss, p::T, y::T) = y * p < one(T) ? 2(p - y) : zero(T)

function value_and_deriv{T<:BlasReal}(::SqrHingeLoss, p::T, y::T)
    yp = y * p
    yp >= one(T) ? (zero(T), zero(T)) : (abs2(one(T) - yp), 2(p - y))
end

## Smoothed HingeLoss
#
#   loss(p, y) := 0              ... y * p > 1 + h
#                 1 - y * p      ... y * p < 1 - h
#                 (1 + h - y * p)^2 / 4h   ... otherwise
#
#  Reference
#
#   O. Chapelle, "Training a Support Vector Machine in the Primal", Neural Computation.
#
immutable SmoothedHingeLoss <: UnivariateLoss
    h::Float64

    function SmoothedHingeLoss(h::Real)
        h > zero(h) || error("h must be a positive value.")
        new(convert(Float64, h))
    end
end

function value{T<:BlasReal}(loss::SmoothedHingeLoss, p::T, y::T)
    h = convert(T, loss.h)
    yp = y * p
    yp >= one(T) + h ? zero(T) :
    yp <= one(T) - h ? one(T) - yp :
    abs2(one(T) + h - yp) / 4h
end

function deriv{T<:BlasReal}(loss::SmoothedHingeLoss, p::T, y::T)
    h = convert(T, loss.h)
    yp = y * p
    yp >= one(T) + h ? zero(T) :
    yp <= one(T) - h ? -y :
    y * (yp - one(T) - h) / 2h
end

function value_and_deriv{T<:BlasReal}(loss::SmoothedHingeLoss, p::T, y::T)
    h = convert(T, loss.h)
    yp = y * p
    if yp >= one(T) + h
        (zero(T), zero(T))
    elseif yp <= one(T) - h
        (one(T) - yp, -y)
    else
        z = yp - one(T) - h
        (abs2(z) / 4h, y * z / 2h)
    end
end

## Logistic loss (for logistic regression)
#
#   loss(p, y) := log(1 + exp(-y * p))
#
immutable LogisticLoss <: UnivariateLoss
end

value{T<:BlasReal}(::LogisticLoss, p::T, y::T) =
    (yp = y * p; yp >= zero(T) ? log1p(exp(-yp)) : log1p(exp(yp)) - yp)

function deriv{T<:BlasReal}(::LogisticLoss, p::T, y::T)
    yp = y * p
    if yp >= zero(T)
        e = exp(-yp)
        -y * e / (one(T) + e)
    else
        -y / (one(T) + exp(yp))
    end
end

function value_and_deriv{T<:BlasReal}(::LogisticLoss, p::T, y::T)
    yp = y * p
    if yp >= zero(T)
        e = exp(-yp)
        (log1p(e), -y * e / (one(T) + e))
    else
        e = exp(yp)
        (log1p(e) - yp, -y / (one(T) + e))
    end
end

## general functions for multivariate loss

grad{T<:BlasReal}(loss::MultivariateLoss, p::StridedVector{T}, y) =
    grad!(loss, zeros(T, length(p)), p, y)

value_and_grad{T<:BlasReal}(loss::MultivariateLoss, p::StridedVector{T}, y) =
    value_and_grad!(loss, zeros(T, length(p)), p, y)

## SumLoss
#
#   loss(p, y) := sum_k loss.intern(p[k], y[k])
#
immutable SumLoss{L<:UnivariateLoss} <: MultivariateLoss
    intern::L
end

SumLoss{L<:UnivariateLoss}(loss::L) = SumLoss{L}(loss)

typealias SumSqrLoss SumLoss{SqrLoss}
SumSqrLoss() = SumLoss{SqrLoss}(SqrLoss())

function value{T<:BlasReal}(s::SumLoss, p::StridedVector{T}, y::StridedVector{T})
    loss = s.intern
    k = length(p)
    @_checkdims k == length(y)
    s = value(loss, p[1], y[1])
    @inbounds for i = 2:k
        s += value(loss, p[i], y[i])
    end
    s
end

function grad!{T<:BlasReal}(s::SumLoss, g::StridedVector{T}, p::StridedVector{T}, y::StridedVector{T})
    loss = s.intern
    k = length(p)
    @_checkdims k == length(g) == length(y)
    @inbounds for i = 1:k
        g[i] = deriv(loss, p[i], y[i])
    end
    g
end

function value_and_grad!{T<:BlasReal}(s::SumLoss, g::StridedVector{T}, p::StridedVector{T}, y::StridedVector{T})
    loss = s.intern
    k = length(p)
    @_checkdims k == length(g) == length(y)
    s, dv = value_and_deriv(loss, p[1], y[1])
    g[1] = dv
    @inbounds for i = 2:k
        v, dv = value_and_deriv(loss, p[i], y[i])
        s += v
        g[i] = dv
    end
    (s, g)
end

## Multinomial logistic loss (for Multinomial logistic regression)
#
#   loss(p, y) := log(sum_k exp(p[k])) - p[y]
#
immutable MultiLogisticLoss <: MultivariateLoss
end

function value{T<:BlasReal}(::MultiLogisticLoss, p::StridedVector{T}, y::Integer)
    k = length(p)
    pmax = maximum(p)
    s = zero(T)
    @inbounds for i = 1:k
        s += exp(p[i] - pmax)
    end
    pmax + log(s) - p[y]
end

function grad!{T<:BlasReal}(::MultiLogisticLoss, g::StridedVector{T}, p::StridedVector{T}, y::Integer)
    k = length(p)
    @_checkdims length(g) == k
    pmax = maximum(p)
    s = zero(T)
    @inbounds for i = 1:k
        pi = exp(p[i] - pmax)
        g[i] = pi
        s += pi
    end
    @inbounds for i = 1:k
        g[i] /= s
    end
    g[y] -= one(T)
    g
end

function value_and_grad!{T<:BlasReal}(::MultiLogisticLoss, g::StridedVector{T}, p::StridedVector{T}, y::Integer)
    k = length(p)
    @_checkdims length(g) == k

    pmax = maximum(p)
    p_y = p[y]  # g and p can be the same array, so have to cache p[y] before p is overwritten

    s = zero(T)
    @inbounds for i = 1:k
        pi = exp(p[i] - pmax)
        g[i] = pi
        s += pi
    end

    @inbounds for i = 1:k
        g[i] /= s
    end
    g[y] -= one(T)
    return (pmax - p_y + log(s), g)
end


text:
module fastPHASE

using Pkg.Artifacts
using SnpArrays
using DelimitedFiles
using CSV
using DataFrames

export fastPHASE_EXE,
    fastphase_estim_param,
    process_fastphase_output,
    fastphase,
    merge_knockoffs_with_original

if Sys.isapple()
    fastPHASE_EXE = joinpath(artifact"fastPHASE", "fastPHASE")
elseif Sys.islinux()
    fastPHASE_EXE = joinpath(artifact"fastPHASE", "fastPHASE")
else
    error("fastPHASE only supports Linux and MacOS")
end
run(`chmod +x $fastPHASE_EXE`) # fastPHASE executable by default is read-only

include("utilities.jl")

end # module


text:
@testset "779.k-th-symbol-in-grammar.jl" begin
    @test kth_grammar(30, 235236) == true
    @test kth_grammar(30, 755232) == true
    @test kth_grammar(1, 1) == 0
    @test kth_grammar(2, 1) == 0
end


text:
using WormlikeChain
using Documenter

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

makedocs(;
    modules=[WormlikeChain],
    authors="Nathan Zimmerberg <nhz2@cornell.edu> and contributors",
    repo="https://github.com/nhz2/WormlikeChain.jl/blob/{commit}{path}#{line}",
    sitename="WormlikeChain.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://nhz2.github.io/WormlikeChain.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/nhz2/WormlikeChain.jl",
    devbranch = "main"
)


text:
module SortMark

export make_df, compute!, reproduce, stat!

using DataFrames: DataFrame, DataFrameRow
using Random: shuffle, shuffle!
using Base.Order
using Base.Sort: Algorithm
using Statistics
using HypothesisTests
using ProgressMeter
include("print.jl")

## Types
const Ints = [Int64, Int8, Int32, Int16, Int128]
const UInts = unsigned.(Ints)
const Floats = [Float64, Float32, Float16]
const BitTypes = vcat(Ints, Floats, UInts, Char, Bool)

matched_UInt(::Type{Float64}) = UInt64
matched_UInt(::Type{Float32}) = UInt32
matched_UInt(::Type{Float16}) = UInt16
matched_UInt(T::Type{<:Signed}) = unsigned(T)
matched_UInt(T::Type{<:Unsigned}) = T
matched_UInt(::Type{Char}) = UInt32

## Special Values
special_values(T) = T
function special_values(T::Type{<:AbstractFloat})
    U = matched_UInt(T)
    T[-π, -1.0, -1/π, 1/π, 1.0, π, -0.0, 0.0, Inf, -Inf, NaN, -NaN,
        prevfloat(T(0)), nextfloat(T(0)), prevfloat(T(Inf)), nextfloat(T(-Inf)),
        reinterpret(T, reinterpret(U, T(NaN)) | 0x73), reinterpret(T, reinterpret(U, -T(NaN)) | 0x37)]
end
function special_values(T::Type{<:Union{Signed, Unsigned}})
    T[17, -T(17), 0, -one(T), 1, typemax(T), typemin(T), typemax(T)-1, typemin(T)+1]
end
special_values(T::Type{Bool}) = [true, false]

bit_random(T::Type{Bool}, len::Integer) = rand(T, len)
bit_random(T::Type, len::Integer) = reinterpret.(T, rand(matched_UInt(T), len))

## Sources (perhaps chagne this to a module-level function with dispatch by symbol)
const sources = Dict(
    :simple => (T, len) -> rand(T, len),
    :special => (T, len) -> rand(special_values(T), len),
    :tricky => (T, len) -> vcat(
        rand(special_values(T), len-2*(len÷3)),
        rand(T, len÷3),
        bit_random(T, len÷3)),
    :small_positive => (T::Type{<:Real}, len) -> rand(T.(1:min(1_000, typemax(T))), len)
)

## Lengths
function lengths(min=1, max=100_000, spacing=3)
    length = 1+round(Integer, (log(max)-log(min))/log(spacing))
    source = exp.(range(log(min), log(max), length=length))
    sort!(collect(Set(round.(Integer, source))))
end

## Orders
const orders = [Forward, Reverse]

##Data frame (this could be much faster if it constructed by column instead of row)
"""
    make_df(algs = [MergeSort, QuickSort]; ...)

Make a dataframe where each row is a sorting task specified by the cartesian product of
keyword arguments.

# Arguments

- `algs::AbstractVector{<:Base.Sort.Algorithm} = [MergeSort, QuickSort]`:
    sorting algorithms to test
- `unstable::AbstractVector{<:Base.Sort.Algorithm} = [QuickSort]`:
    which of the algorithms are allowed to be unstable

Axes of the cartesian product
- `Types::AbstractVector{<:Type} = SortMark.BitTypes`:
    element type
- `lens::AbstractVector{<:Integer} = SortMark.lengths()`:
    number of elements to be sorted
- `orders::AbstractVector{<:Ordering} = SortMark.orders`:
    orders to sort by
- `sources::AbstractDict{<:Any, <:Function} = SortMark.sources`:
    generation procedures to create input data

Benchmarking time
- `seconds::Union{Real, Nothing} = .005`:
    maximum benchmarking time for each row. Compute sum(df.seconds) for an estimated
    benchmarking runtime
- `samples::Union{Nothing, Integer} = nothing`:
    maximum number of samples for each row. Compute sum(df.seconds) for an estimated
    benchmarking runtime

Setting `seconds` or `samples` to `nothing` removes that limit.
"""
function make_df(algs::AbstractVector{<:Algorithm} = [MergeSort, QuickSort];
    Types::AbstractVector{<:Type} = BitTypes,
    lens::AbstractVector{<:Integer} = lengths(),
    orders::AbstractVector{<:Ordering} = orders,
    sources::AbstractDict{<:Any, <:Function} = sources, unstable = [QuickSort],
    seconds::Union{Real, Nothing} = .005,
    samples::Union{Integer, Nothing} = seconds==nothing ? 1 : nothing)

    df = DataFrame(ContainerType=Type{<:AbstractVector}[], Type=Type[], len=Int[],
        order=Ordering[], source_key=Symbol[], source=Function[], algs=Vector{<:Algorithm}[], unstable=Vector{<:Algorithm}[],
        evals=Int[], samples=Union{Nothing, Int}[], seconds=Union{Nothing, Float64}[], data=DataFrame[],
        errors=Dict[])

    for Type in Types, len in lens, order in orders, (source_key, source) in pairs(sources)
        if source_key == :small_positive && !(Type <: Real); continue; end
        push!(df, (Vector, Type, len, order, source_key, source, algs, unstable,
                   max(1, 1_000 ÷ max(len, 1)), samples, seconds, DataFrame([Float64[] for _ in algs], Symbol.(algs)),
                   Dict()))
    end
    df
end

##Compute data
function target!(vec, alg, order)
    sort!(vec, alg=alg, order=order)
end

const PASSED_TESTS = Ref(0)

"""
    compute!(df::DataFrame; verbose=true, fail_fast=true)

Test and run benchmarks for every row in `df`.

Benchmark results are saved in df.data.
"""
function compute!(df::DataFrame; verbose=true, fail_fast=true)
    start_time = time()
    benchmark_start_time = sum(sum(sum.(eachcol(d))) for d in df.data)
    start_passed_tests = PASSED_TESTS[]
    empty!.(df.errors)

    rows = shuffle(eachrow(df))#TODO make this appear immediately
    nominal_runtime = sum([s==nothing ? .001 : s for s in df.seconds])

    dt = max(.1, min(1, nominal_runtime/100))
    message = "$(round(Integer, nominal_runtime))s, Progress: "

    try
        @showprogress dt message for row in rows
            compute!(row, fail_fast)
        end
    catch
        printstyled("\nFailed test set stored as SortMark.fail. Reproduce it's error with reproduce().\n", color=Base.debug_color())
        rethrow()
    end

    np, nf, ne = PASSED_TESTS[]-start_passed_tests, 0, 0
    for dict in df.errors, vect in values(dict), err in vect
        if err isa AssertionError
            nf += 1
        else
            ne += 1
        end
    end

    end_time = time()
    benchmark_end_time = sum(sum(sum.(eachcol(d))) for d in df.data)

    if !all(isempty.(df.errors))
        global fail = deepcopy(df[findfirst((!).(isempty.(df.errors))), :])

        color = Base.error_color()
        printstyled("ERROR: Some tests did not pass: $np passed, $nf failed, $ne errored.\n", color=color)
        verbose && Print._print(df, np, end_time-start_time, benchmark_end_time-benchmark_start_time, color)
        printstyled("SortMark.fail is an example of a failed set of tests. Reproduce it's error with reproduce().\n", color=Base.debug_color())
    elseif verbose
        Print._print(df, np, end_time-start_time, benchmark_end_time-benchmark_start_time, :normal)
    end

    df
end

function compute!(row::DataFrameRow, fail_fast=true)
    Base.require_one_based_indexing(row.algs)

    gen() = row.ContainerType(row.source(row.Type, row.len))

    if row.samples == nothing && row.seconds == nothing
        println("WARNING: row=$row has unbounded samples=$samples and seconds=$seconds.
Because of this absence of an end condition it will hang.")
    end

    sample = 1

    start_time = time()
    while (row.samples == nothing || sample <= row.samples) &&
          (row.seconds == nothing || time() < start_time+row.seconds)

        perm = shuffle!(collect(axes(row.algs, 1)))

        times = trial(row, row.algs[perm], [gen() for _ in 1:row.evals], fail_fast)

        push!(row.data, times[invperm(perm)])

        sample += 1
    end

    row
end
function trial(row, algs, inputs, fail_fast)
    Base.require_one_based_indexing(algs)

    stable_output = nothing
    test_index = rand(eachindex(inputs))

    times = []
    order = row.order
    for alg in algs
        cinputs = copy.(inputs)
        outputs = similar(inputs)
        t = @elapsed for i in eachindex(inputs, cinputs, outputs)
            try
                outputs[i] = target!(cinputs[i], alg, order)
            catch x
                log_error!(row, x, alg, copy(inputs[i]), fail_fast)
            end
        end
        push!(times, t/length(inputs))

        try
            test_sorted(outputs[test_index], cinputs[test_index], inputs[test_index], row.order)
            if alg ∉ row.unstable
                if stable_output == nothing
                    stable_output = outputs[test_index]
                else
                    test_matched(stable_output, outputs[test_index])
                end
            end
        catch x
            log_error!(row, x, alg, copy(inputs[test_index]), fail_fast)
        end
    end

    times
end

function test_sorted(output, cinput, input, order)
    @assert cinput === output
    PASSED_TESTS[] += 1
    @assert issorted(output, order=order)
    PASSED_TESTS[] += 1
    @assert typeof(input) == typeof(output)
    PASSED_TESTS[] += 1
    @assert axes(input) == axes(output)
    PASSED_TESTS[] += 1
    nothing
end

function test_matched(a, b)
    @assert all(a .=== b)
    PASSED_TESTS[] += 1
    nothing
end

function log_error!(row, exception, alg, input, fail_fast)
    entry = (exception, input)
    if alg ∉ keys(row.errors)
        row.errors[alg] = []
    end
    push!(row.errors[alg], entry)
    if fail_fast || exception isa InterruptException
        global fail = row
        throw(exception)
    end
end

fail = nothing
function reproduce(row=fail, alg=nothing)
    if fail == nothing || isempty(row.errors)
        return "no errors"
    end
    if alg == nothing
        alg = first(keys(row.errors))
    end

    x, input = first(row.errors[alg])
    cinput = copy(input)
    println("sort!($cinput, "* (row.order != Forward ? "order=$(row.order), " : "")*"alg=$alg)")
    output = target!(cinput, alg, row.order)
    test_sorted(output, cinput, input, row.order)
    test_matched(output, sort(input, alg=MergeSort))

    "no errors"
end

##Statistics
"""
    function stat!(df, a=1, b=2)

Compute comparative stats for the `a`th and `b`th algorithm tested in `df`.

Returns the 95% confidence interval for the ratio of runtimes a/b for each row.
"""
function stat!(df, a=1, b=2)
    df.log_test = [length(eachrow(d)) <= 1 ? missing : OneSampleTTest(d[!,a]-d[!,b]) for d in [log.(frame) for frame in df.data]]
    df.pvalue = [ismissing(t) ? missing : pvalue(t) for t in df.log_test]
    df.point_estimate = [ismissing(t) ? missing : exp(mean(confint(t))) for t in df.log_test]
    df.confint = [ismissing(t) ? missing : exp.(confint(t)) for t in df.log_test]
end

end #module

max_stars_repo_name: LilithHafner/SortMark

max_issues_repo_name: LilithHafner/SortMark


text:
using PiecewiseDeterministicMarkovProcesses, Random, DifferentialEquations
	const PDMP = PiecewiseDeterministicMarkovProcesses

function F_tcp!(ẋ, xc, xd, parms, t)
	if mod(xd[1],2)==0
		 ẋ[1] =  1.0
		 ẋ[2] = -1.0 * xc[2]
	else
		 ẋ[1] = -1.0 * xc[1]
		 ẋ[2] =  1.0
	end
	nothing
end

R(x) = x

function R_tcp!(rate, xc, xd, parms, t, issum::Bool)
		rate[1] = R(xc[1]) +  R(xc[2])
		rate[2] = parms[1] * xd[1] * xc[1]
	if issum == false
		return 0.
	else
		return sum(rate)
	end
end

xc0 = [0.05, 0.075]
xd0 = [0, 1]

nu_tcp = [1 0;0 -1]
parms = [0.1]
tf = 10000.0
nj = 1000

Random.seed!(1234)
	problem = PDMP.PDMPProblem(F_tcp!, R_tcp!, nu_tcp, xc0, xd0, parms, (0.0, tf))
	result1 = @time PDMP.solve(problem, CHV(Tsit5()); n_jumps = nj, save_positions = (false, false))

Random.seed!(1234)
	result2 = PDMP.solve(problem, CHV(:cvode); n_jumps = nj, save_positions = (false, false))

Random.seed!(1234)
	result3 = PDMP.solve(problem, CHV(:lsoda); n_jumps = 2, save_positions = (false, false))

#test auto-differentiation
# result4 = @time PDMP.solve(problem, CHV(Rodas5()))

####################################################################################################
# DEBUG DEBUG
#
# algo = CHV(Tsit5())
# xd1 = zeros(Float64, length(xc0)+1)
# xdd1 = similar(xd1)
#
# J = zeros(3,3)
# algo(xdd1,xd1,problem.caract,0.)
#
# vf = (dx,x) -> algo(dx,x,problem.caract,0.)
# #it works!
# vf(xdd1,xd1)
#
# ForwardDiff.jacobian!(J, vf, xdd1, xd1)


text:
# extracted from Meshes.jl

"""
Evaluation method used to obtain a point along
a Bézier curve from a parametric expression.
"""
abstract type BezierEvalMethod end

"""
Approximate evaluation using Horner's method.
Recommended for a large number of control points,
if you can afford a precision loss.
See https://en.wikipedia.org/wiki/Horner%27s_method.
"""
struct Horner <: BezierEvalMethod end

struct BezierCurve{M<:BezierEvalMethod} <: Curve
  method::M
end

BezierCurve() = BezierCurve(Horner())

degree(p) = length(p) - 1

"""
Apply Horner's method on the monomial representation of the
Bézier curve B = ∑ᵢ aᵢtⁱ with i ∈ [0, n], n the degree of the
curve, aᵢ = binomial(n, i) * pᵢ * t̄ⁿ⁻ⁱ and t̄ = (1 - t).
Horner's rule recursively reconstructs B from a sequence bᵢ
with bₙ = aₙ and bᵢ₋₁ = aᵢ₋₁ + bᵢ * t until b₀ = B.
"""
function (curve::BezierCurve{Horner})(t, points)
  if t < 0 || t > 1
    throw(DomainError(t, "b(t) is not defined for t outside [0, 1]."))
  end
  d = length(points)
  T = eltype(eltype(points))
  t̄ = one(T) - t
  n = degree(points)
  pₙ = last(points)
  aₙ = pₙ

  # initialization with i = n + 1, so bᵢ₋₁ = bₙ = aₙ
  bᵢ₋₁ = aₙ
  cᵢ₋₁ = one(T)
  t̄ⁿ⁻ⁱ = one(T)
  for i in n:-1:1
    cᵢ₋₁ *= i / (n - i + one(T))
    pᵢ₋₁ = points[i]
    t̄ⁿ⁻ⁱ *= t̄
    aᵢ₋₁ = cᵢ₋₁ * pᵢ₋₁ * t̄ⁿ⁻ⁱ
    bᵢ₋₁ = aᵢ₋₁ + bᵢ₋₁ * t
  end

  b₀ = bᵢ₋₁
  b₀
end


text:
using JSServe
using WGLMakie
using GeometryBasics
using FileIO
using WGLMakie: px
using JSServe: @js_str, onjs, App, Slider
using JSServe.DOM

JSServe.browser_display()

set_theme!(resolution=(1200, 800))

app = App() do
    return hbox(
        vbox(
            scatter(1:4, color=1:4),
            scatter(1:4, color=rand(RGBAf, 4)),
            scatter(1:4, color=rand(RGBf, 4)),
            scatter(1:4, color=:red),
            scatter(1:4)
        ),
        vbox(
            scatter(1:4, marker='☼'),
            scatter(1:4, marker=['☼', '◒', '◑', '◐']),
            scatter(1:4, marker="☼◒◑◐"),
            # scatter(1:4, marker=rand(RGBf, 10, 10), markersize=20px),
            scatter(1:4, markersize=20px),
            scatter(1:4, markersize=20, markerspace=Pixel),
            scatter(1:4, markersize=LinRange(20, 60, 4), markerspace=Pixel),
            scatter(1:4, marker='▲', markersize=0.3, rotations=LinRange(0, pi, 4)),
        )
    )
end

display(app)

app = App() do
    return DOM.div(
        meshscatter(1:4, color=1:4),
        meshscatter(1:4, color=rand(RGBAf, 4)),
        meshscatter(1:4, color=rand(RGBf, 4)),
        meshscatter(1:4, color=:red),
        meshscatter(rand(Point3f, 10), color=rand(RGBf, 10)),
        meshscatter(rand(Point3f, 10), marker=Pyramid(Point3f(0), 1f0, 1f0)),
    )
end

display(app)

app = App() do
    x = Point2f[(1, 1), (2, 2), (3, 2), (4, 4)]
    points = connect(x, LineFace{Int}[(1, 2), (2, 3), (3, 4)])
    return DOM.div(
        linesegments(1:4),
        linesegments(1:4, linestyle=:dot),
        linesegments(1:4, linestyle=[0.0, 1.0, 2.0, 3.0, 4.0]),
        linesegments(1:4, color=1:4),
        linesegments(1:4, color=rand(RGBf, 4), linewidth=10),
        linesegments(points)
    )
end

display(app)

app = App() do
    x = Point2f[(1, 1), (2, 2), (3, 2), (4, 4)]
    points = connect(x, LineFace{Int}[(1, 2), (2, 3), (3, 4)])
    return DOM.div(
        lines(1:4),
        lines(1:4, linestyle=:dot),
        lines(1:4, linestyle=[0.0, 1.0, 2.0, 3.0, 4.0]),
        lines(1:4, color=1:4),
        lines(1:4, color=rand(RGBf, 4), linewidth=10),
        lines(points)
    )
end

display(app)

app = App() do
    data = AbstractPlotting.peaks()
    return hbox(vbox(
        surface(-10..10, -10..10, data, show_axis=false),
        surface(-10..10, -10..10, data, color=rand(size(data)...))),
        vbox(surface(-10..10, -10..10, data, color=rand(RGBf, size(data)...)),
        surface(-10..10, -10..10, data, colormap=:magma, colorrange=(0.0, 2.0)),
    ))
end

display(app)

app = App() do
    return vbox(
        image(rand(10, 10)),
        heatmap(rand(10, 10)),
        image(rand(RGBAf, 10, 10)),
        heatmap(rand(RGBAf, 10, 10)),
        image(rand(RGBf, 10, 10)),
        heatmap(rand(RGBf, 10, 10)),
    )
end

display(app)

app = App() do
    return hbox(vbox(
        volume(rand(4, 4, 4), isovalue=0.5, isorange=0.01, algorithm=:iso),
        volume(rand(4, 4, 4), algorithm=:mip),
        volume(1..2, -1..1, -3..(-2), rand(4, 4, 4), algorithm=:absorption)),
        vbox(
        volume(rand(4, 4, 4), algorithm=Int32(5)),
        volume(rand(RGBAf, 4, 4, 4), algorithm=:absorptionrgba),
        contour(rand(4, 4, 4)),
    ))
end

display(app)

app = App() do
    cat = FileIO.load(MakieGallery.assetpath("cat.obj"))
    tex = FileIO.load(MakieGallery.assetpath("diffusemap.png"))
    return hbox(vbox(
        AbstractPlotting.mesh(Circle(Point2f(0), 1f0)),
        AbstractPlotting.poly(decompose(Point2f, Circle(Point2f(0), 1f0)))), vbox(
        AbstractPlotting.mesh(cat, color=tex),
        AbstractPlotting.mesh([(0.0, 0.0), (0.5, 1.0), (1.0, 0.0)]; color=[:red, :green, :blue], shading=false)
    ))
end

display(app)

function n_times(f, n=10, interval=0.5)
    obs = Observable(f(1))
    @async for i in 2:n
        try
            obs[] = f(i)
            sleep(interval)
        catch e
            @warn "Error!" exception=CapturedException(e, Base.catch_backtrace())
        end
    end
    return obs
end

function n_times(f, n=10, interval=1)
    obs = Observable(f(1))
    @async for i in 2:n
        try
            obs[] = f(i)
            sleep(interval)
        catch e
            @warn "Error!" exception=CapturedException(e, Base.catch_backtrace())
        end
    end
    return obs
end

app = App() do
    s1 = annotations(n_times(i-> map(j-> ("$j", Point2f(j*30, 0)), 1:i)), textsize=20,
                      limits=Rect2f(30, 0, 320, 50))
    s2 = scatter(n_times(i-> Point2f.((1:i).*30, 0)), markersize=20px,
                  limits=Rect2f(30, 0, 320, 50))
    s3 = linesegments(n_times(i-> Point2f.((2:2:2i).*30, 0)), limits=Rect2f(30, 0, 620, 50))
    s4 = lines(n_times(i-> Point2f.((2:2:2i).*30, 0)), limits=Rect2f(30, 0, 620, 50))
    return hbox(s1, s2, s3, s4)
end

display(app)

app = App() do
    outer_padding = 30
    scene, layout = layoutscene(
        outer_padding, resolution=(1200, 700),
        backgroundcolor=RGBf(0.98, 0.98, 0.98)
    )
    ax1 = layout[1, 1] = LAxis(scene, title="Sine")
    xx = 0:0.2:4pi
    line1 = lines!(ax1, sin.(xx), xx, color=:red)
    scat1 = scatter!(ax1, sin.(xx) .+ 0.2 .* randn.(), xx,
        color=(:red, 0.5), markersize=15px, marker='■')
    return scene
end

display(app)

app = App() do
    scene = Scene(resolution=(1000, 1000));
    campixel!(scene);

    maingl = GridLayout(scene, alignmode=Outside(30))

    las = Array{LAxis,2}(undef, 4, 4)

    for i in 1:4, j in 1:4
        las[i, j] = maingl[i, j] = LAxis(scene)
    end

    las[4, 1].attributes.aspect = AxisAspect(1)
    las[4, 2].attributes.aspect = AxisAspect(2)
    las[4, 3].attributes.aspect = AxisAspect(0.5)
    las[4, 4].attributes.aspect = nothing
    las[1, 1].attributes.maxsize = (Inf, Inf)
    las[1, 2].attributes.aspect = nothing
    las[1, 3].attributes.aspect = nothing

    subgl = gridnest!(maingl, 1, 1)
    cb1 = subgl[:, 2] = LColorbar(scene, width=30, height=Relative(0.66))

    subgl2 = gridnest!(maingl, 1:2, 1:2)
    cb2 = subgl2[:, 3] = LColorbar(scene, width=30, height=Relative(0.66))

    subgl3 = gridnest!(maingl, 1:3, 1:3)
    cb3 = subgl3[:, 4] = LColorbar(scene, width=30, height=Relative(0.66))

    subgl4 = gridnest!(maingl, 1:4, 1:4)
    cb4 = subgl4[:, 5] = LColorbar(scene, width=30, height=Relative(0.66))
    return scene

end

display(app)

app = App() do
    sl = JSServe.Slider(1:10)
    rect = Rect2f(0, -5, 1025, 10)
    chars = [collect('a':'z'); 0:9;]
    char2 = [collect('A':'Z'); 0:9;]
    tex1 = map(x->map(j-> ("$(chars[rand(1:length(chars))])", Point2f(j*30, 0)), 1:36), sl)
    tex2 = map(x->map(j-> ("$(char2[rand(1:length(char2))])", Point2f(j*30, 1)), 1:36), sl)
    scene = annotations(tex1, textsize=20, limits=rect, show_axis=false)
    annotations!(scene, tex2, textsize=20, limits=rect, show_axis=false)
    return DOM.div(sl, scene)
end

display(app)


text:
@params mutable struct GESOResult{T}
    topology::AbstractVector{T}
    objval::T
    change::T
    converged::Bool
    fevals::Int
end

"""
The GESO algorithm, see [LiuYiLiShen2008](@cite).
"""
struct GESO <: TopOptAlgorithm
    comp::Compliance
    vol::Volume
    vol_limit::Any
    filter::Any
    vars::AbstractVector
    topology::AbstractVector
    Pcmin::Any
    Pcmax::Any
    Pmmin::Any
    Pmmax::Any
    Pen::Any
    string_length::Int
    var_volumes::AbstractVector
    cum_var_volumes::AbstractVector
    order::AbstractVector{Int}
    genotypes::BitArray{2}
    children::BitArray{2}
    var_black::BitVector
    maxiter::Int
    penalty::Any
    sens::AbstractVector
    old_sens::AbstractVector
    obj_trace::MVector{10}
    tol::Any
    sens_tol::Any
    result::GESOResult
end
Base.show(::IO, ::MIME{Symbol("text/plain")}, ::GESO) = println("TopOpt GESO algorithm")

function GESO(
    comp::Compliance,
    vol::Volume,
    vol_limit,
    filter;
    maxiter=1000,
    tol=0.001,
    p=3.0,
    Pcmin=0.6,
    Pcmax=1.0,
    Pmmin=0.5,
    Pmmax=1.0,
    Pen=3.0,
    sens_tol=tol / 100,
    string_length=4,
    k=10,
)
    penalty = comp.solver.penalty
    setpenalty!(penalty, p)
    solver = comp.solver
    T = eltype(solver.vars)
    nel = getncells(solver.problem.ch.dh.grid)
    @unpack white, black = solver.problem
    nvars = nel - sum(black) - sum(white)
    vars = zeros(T, nvars)

    topology = zeros(T, nel)
    result = GESOResult(topology, T(NaN), T(NaN), false, 0)
    sens = zeros(T, nvars)
    old_sens = zeros(T, nvars)
    obj_trace = zeros(MVector{k,T})
    var_volumes = vol.cellvolumes[.!black .& .!white]
    cum_var_volumes = zeros(T, nvars)
    order = zeros(Int, nvars)
    genotypes = trues(string_length, nvars)
    children = trues(string_length, nvars)
    var_black = trues(nvars)

    return GESO(
        comp,
        vol,
        vol_limit,
        filter,
        vars,
        topology,
        Pcmin,
        Pcmax,
        Pmmin,
        Pmmax,
        Pen,
        string_length,
        var_volumes,
        cum_var_volumes,
        order,
        genotypes,
        children,
        var_black,
        maxiter,
        penalty,
        sens,
        old_sens,
        obj_trace,
        tol,
        sens_tol,
        result,
    )
end

function Utilities.setpenalty!(b::GESO, p::Number)
    b.penalty.p = p
    return b
end

function get_progress(current_volume, total_volume, design_volume)
    return clamp(
        min(
            (total_volume - current_volume) / (total_volume - design_volume),
            current_volume / design_volume,
        ),
        0,
        1,
    )
end

function get_probs(b::GESO, Prg)
    return (
        b.Pcmin + (b.Pcmax - b.Pcmin) * Prg^b.Pen, b.Pmmin + (b.Pmmax - b.Pmmin) * Prg^b.Pen
    )
end

function crossover!(children, genotypes, i, j)
    for k in 1:size(genotypes, 1)
        r = rand()
        if r < 0.5
            children[k, i] = genotypes[k, i]
        else
            children[k, i] = genotypes[k, j]
        end
    end
    return nothing
end

function update!(var_black, children, genotypes, Pc, Pm, high_class, mid_class, low_class)
    topology_changed = false
    while !topology_changed
        for i in high_class
            r = rand()
            j = i
            if length(high_class) > 1
                if r < Pc
                    while i == j
                        j = rand(high_class)
                    end
                elseif r < 0.5 + 0.5 * Pc
                    j = rand(mid_class)
                else
                    j = rand(low_class)
                end
            else
                if r < 0.5
                    j = rand(mid_class)
                else
                    j = rand(low_class)
                end
            end
            crossover!(children, genotypes, i, j)
        end
        for i in mid_class
            r = rand()
            j = i
            if length(mid_class) > 1
                if r < Pc
                    while i == j
                        j = rand(mid_class)
                    end
                elseif r < 0.5 + 0.5 * Pc
                    j = rand(high_class)
                else
                    j = rand(low_class)
                end
            else
                if r < 0.5 + 0.5 * Pc
                    j = rand(high_class)
                else
                    j = rand(low_class)
                end
            end
            crossover!(children, genotypes, i, j)
        end
        for i in low_class
            r = rand()
            j = i
            if length(low_class) > 1
                if r < Pc
                    while i == j
                        j = rand(low_class)
                    end
                elseif r < 0.5 + 0.5 * Pc
                    j = rand(mid_class)
                else
                    j = rand(high_class)
                end
            else
                if r < 0.5
                    j = rand(mid_class)
                else
                    j = rand(high_class)
                end
            end
            crossover!(children, genotypes, i, j)
        end
        genotypes .= children

        for i in high_class
            for j in 1:size(genotypes, 1)
                r = rand()
                if r < Pm && !genotypes[j, i]
                    genotypes[j, i] = !genotypes[j, i]
                end
            end
            if any(@view genotypes[:, i]) != var_black[i]
                var_black[i] = !var_black[i]
                topology_changed = true
            end
        end
        for i in mid_class
            for j in 1:size(genotypes, 1)
                r = rand()
                if r < Pm && genotypes[j, i]
                    genotypes[j, i] = !genotypes[j, i]
                end
            end
            if any(@view genotypes[:, i]) != var_black[i]
                var_black[i] = !var_black[i]
                topology_changed = true
            end
        end
        for i in low_class
            for j in 1:size(genotypes, 1)
                r = rand()
                if r < Pm && genotypes[j, i]
                    genotypes[j, i] = !genotypes[j, i]
                end
            end
            if any(@view genotypes[:, i]) != var_black[i]
                var_black[i] = !var_black[i]
                topology_changed = true
            end
        end
    end

    return var_black
end

function (b::GESO)(x0=copy(b.comp.solver.vars); seed=NaN)
    @unpack sens, old_sens, tol, maxiter = b
    @unpack obj_trace, topology, sens_tol, vars = b
    @unpack Pcmin, Pcmax, Pmmin, Pmmax, Pen = b
    @unpack string_length, genotypes, children, var_black = b
    @unpack cum_var_volumes, var_volumes, order = b
    @unpack varind, black, white = b.comp.solver.problem
    @unpack total_volume, cellvolumes, fixed_volume = b.vol
    T = eltype(x0)
    V = b.vol_limit
    design_volume = V * total_volume

    nel = length(x0)
    nvars = length(vars)

    # Set seed
    isnan(seed) || Random.seed!(seed)

    # Initialize the topology
    for i in 1:length(topology)
        if black[i]
            topology[i] = 1
        elseif white[i]
            topology[i] = 0
        else
            topology[i] = round(x0[varind[i]])
            vars[varind[i]] = topology[i]
        end
    end

    check(x) = x > design_volume - fixed_volume
    #rrmax = clamp(1 - design_volume/current_volume, 0, 1)
    current_volume = dot(vars, var_volumes) + fixed_volume
    vol = current_volume / total_volume
    # Main loop
    change = T(1)
    iter = 0
    while (change > tol || vol > V) && iter < maxiter
        iter += 1
        if iter > 1
            old_sens .= sens
        end
        for j in max(2, 10 - iter + 2):10
            obj_trace[j - 1] = obj_trace[j]
        end
        obj_trace[10], pb = Zygote.pullback(x -> b.comp(b.filter(x)), vars)
        sens = pb(1.0)[1]
        rmul!(sens, -1)
        if iter > 1
            @. sens = (sens + old_sens) / 2
        end

        # Classify the cells by their sensitivities
        sortperm!(order, sens; rev=true)
        accumulate!(+, cum_var_volumes, view(var_volumes, order))
        N1 = findfirst(check, cum_var_volumes) - 1
        N2 = (nel - N1) ÷ 2
        N3 = nvars - N1 - N2
        high_class = @view order[1:N1]
        mid_class = @view order[(N1 + 1):(N1 + N2)]
        low_class = @view order[(N1 + N2 + 1):end]

        # Crossover and mutation
        Prg = get_progress(current_volume, total_volume, design_volume)
        Pc, Pm = get_probs(b, Prg)
        vars .= update!(
            var_black, children, genotypes, Pc, Pm, high_class, mid_class, low_class
        )

        # Update crossover and mutation probabilities
        current_volume = dot(vars, var_volumes) + fixed_volume
        vol = current_volume / total_volume

        if iter >= 10
            l = sum(@view obj_trace[1:5])
            h = sum(@view obj_trace[6:10])
            change = abs(l - h) / h
        end
    end

    for i in 1:length(topology)
        if black[i]
            topology[i] = 1
        elseif white[i]
            topology[i] = 0
        else
            topology[i] = vars[varind[i]]
        end
    end

    b.result.objval = obj_trace[10]
    b.result.change = change
    b.result.converged = change <= tol
    b.result.fevals = iter

    return b.result
end


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# =================================================================================================
# Test writing a SPIX image.
@testset "pix_write_spix/pix_read_spix" begin
    local pix = pix_with()

    # ---------------------------------------------------------------------------------------------
    # To/From a file.
    local filename = safe_tmp_file()
    @test pix_write_spix(filename, pix) == true
    @test pix_read_spix(filename) != nothing
    rm(filename)

    # ---------------------------------------------------------------------------------------------
    # To/From an IOBuffer.
    local buffer = IOBuffer()
    @test pix_write_spix(buffer, pix) == true
    @test position(buffer) != 0
    seekstart(buffer)
    @test pix_read_spix(buffer) != nothing

    # ---------------------------------------------------------------------------------------------
    # To/From an IOStream
    local filename = safe_tmp_file()
    local file     = open(filename, read=true, write=true)
    @test pix_write_spix(file, pix) == true
    @test position(file) != 0
    seekstart(file)
    @test pix_read_spix(file) != nothing
    close(file)
    rm(filename)

    # ---------------------------------------------------------------------------------------------
    # To/From a byte array.
    local data  = spix_with()
    local bytes = pix_write_spix(pix)
    @test length(bytes) > 0
    @test pix_read_spix(bytes) != nothing
end

# =================================================================================================
# Test reading empty data.
@testset "pix_read_spix empty data" begin
    @suppress begin
        # -----------------------------------------------------------------------------------------
        # From a file.
        local filename = safe_tmp_file()
        @test pix_read_spix(filename) == nothing
        @test pix_read_spix("xyzzy.spix") == nothing
        rm(filename)

        # -----------------------------------------------------------------------------------------
        # From an IOBuffer.
        local buffer = IOBuffer()
        @test pix_read_spix(buffer) == nothing

        # -----------------------------------------------------------------------------------------
        # From an IOStream
        local filename = safe_tmp_file()
        local file     = open(filename, read=true, write=true)
        @test pix_read_spix(file) == nothing
        close(file)
        rm(filename)

        # -----------------------------------------------------------------------------------------
        # To/From a byte array.
        local bytes = Vector{UInt8}()
        @test pix_read_spix(bytes) == nothing
    end
end

# =================================================================================================
# Test writing a Pix that has been freed.
@testset "pix_write_spix /w freed Pix" begin
    local err = "Pix has been freed."
    local pix = pix_with()
    pix_delete!(pix)

    # ---------------------------------------------------------------------------------------------
    # To a file.
    local filename = safe_tmp_file()
    @test (@test_logs (:error, err) pix_write_spix(filename, pix)) == false
    rm(filename)

    # ---------------------------------------------------------------------------------------------
    # To an IOBuffer.
    local buffer = IOBuffer()
    @test (@test_logs (:error, err) pix_write_spix(buffer, pix)) == false

    # ---------------------------------------------------------------------------------------------
    # To a byte array.
    @test (@test_logs (:error, err) pix_write_spix(pix)) == nothing
end