text:
"""
    Geom.band[(; orientation=:vertical)]

Draw bands across the plot canvas with a horizontal span specifed by `xmin` and `xmax` if `orientation` is `:vertical`, or a vertical span specified by `ymin` and `ymax` if the `orientation` is `:horizontal`.

This geometry is equivalent to [`Geom.rect`](@ref) with [`Stat.band`](@ref).
"""
band(; orientation=:vertical) = RectangularBinGeometry(Stat.band(orientation)) #TODO: use RectangularGeometry when it becomes available.

"""
    Geom.hband[()]

Draw horizontal bands across the plot canvas with a vertical span specified by `ymin` and `ymax` aesthetics.

This geometry is equivalent to [`Geom.band`](@ref) with `orientation` set to `:vertical`.
"""
hband() = band(orientation = :horizontal)

"""
    Geom.vband[()]

Draw vertical bands across the plot canvas with a horizontal span specified by `xmin` and `xmax` aesthetics.

This geometry is equivalent to [`Geom.band`](@ref).
"""
const vband = band


text:
using LinearAlgebra
using Makie

"""
λmap(p)
Sends a point on a 3-sphere to a point in the plane x₄=0 with the given point
in the form of 2 complex numbers representing a unit quaternion. This is the
stereographic projection of a 3-sphere. S³ ↦ R³
"""
function λmap(p)
    [real(p[1]), imag(p[1]), real(p[2])] ./ (1 - imag(p[2]))
end

"""
λ⁻¹map(p)
Sends a point on the plane back to a point on a unit sphere with the given
point. This is the inverse stereographic projection of a 3-sphere.
"""
function λ⁻¹map(p)
    x₁ = 2p[1] / (1 + p[1]^2 + p[2]^2 + p[3]^2)
    x₂ = 2p[2] / (1 + p[1]^2 + p[2]^2 + p[3]^2)
    x₃ = 2p[3] / (1 + p[1]^2 + p[2]^2 + p[3]^2)
    x₄ = (-1 + p[1]^2 + p[2]^2 + p[3]^2) / (1 + p[1]^2 + p[2]^2 + p[3]^2)
    [Complex(x₁, x₂), Complex(x₃, x₄)]
end

"""
S¹action(α, p)
Performs a group action corresponding to moving along the circumference of a
circle with the given angle and the point in the form of 2 complex numbers
representing a unit quaternion.
"""
function S¹action(α, p)
    [exp(im * α) * p[1], exp(im * α) * p[2]]
end

"""
convert_to_cartesian(p)
Converts a point in the geographic coordinate system to a point in the
cartesian one with the given point in radians.
"""
function convert_to_cartesian(p)
    [cos(p[2]) * cos(p[1]),
     cos(p[2]) * sin(p[1]),
     sin(p[2])]
end

"""
convert_to_geographic(p)
Converts a point in the cartesian coordinate system to a point in the
geographic one with the given point.
"""
function convert_to_geographic(p)
    r = sqrt(p[1]^2 + p[2]^2 + p[3]^2)
    if p[1] > 0
          ϕ = atan(p[2] / p[1])
    elseif p[2] > 0
          ϕ = atan(p[2] / p[1]) + pi
    else
          ϕ = atan(p[2] / p[1]) - pi
    end
    θ = asin(p[3] / r)
    [ϕ, θ]
end

"""
get_center(A, B, C)

Finds the center point of the fiber circle under stereographic projection
with the given 3 points on the circle circumference.
"""
function get_center(A, B, C)
    a = LinearAlgebra.norm(B - C)
    b = LinearAlgebra.norm(A - C)
    c = LinearAlgebra.norm(A - B)
    numerator = a^2 * (b^2 + c^2 - a^2) * A +
                b^2 * (a^2 + c^2 - b^2) * B +
                c^2 * (a^2 + b^2 - c^2) * C
    denominator = a^2 * (b^2 + c^2 - a^2) +
                  b^2 * (a^2 + c^2 - b^2) +
                  c^2 * (a^2 + b^2 - c^2)
    numerator / denominator
end

"""
get_flower(;N=4, A=.5, B=-pi/7, P=pi/2, Q=0, number=300)

Calculates the x, y and z points of a flower in the base space.
with the given number of petals N, the fattness of the petals A,
the height of the petals B, the latitude of the flower P,
the rotation of the flower Q, and the total number of points in the grid.
"""
function get_flower(;N=4, A=.5, B=-pi/7, P=pi/2, Q=0, number=300)
    N = 6
    A = .5
    B = -pi/7
    P = pi/3
    Q = 0
    t = range(0, stop = 2pi, length = number)
    az = 2pi .* t + A .* cos.(N .* 2pi .* t) .+ Q
    po = B .* sin.(N .* 2pi .* t) .+ P
    x = cos.(az).*sin.(po)
    y = sin.(az).*sin.(po)
    z = cos.(po)
    points = Array{Float64}(undef, number, 2)
    for i in 1:number
        points[i, :] = convert_to_geographic([x[i], y[i], z[i]])
    end
    points
end

"""
build_surface(scene, points, color; transparency, shading)

Builds a surface with the given scene, points, color, transparency and shading.
"""
function build_surface(scene,
                       points,
                       color;
                       transparency = false,
                       shading = true)
    surface!(scene,
             @lift($points[:, :, 1]),
             @lift($points[:, :, 2]),
             @lift($points[:, :, 3]),
             color = color,
             transparency = transparency,
             shading = shading)
end

"""
rotate3D_geographic(point, q)

Rotates a point in the 3D space with the given point and the unit quaternion.
"""
function rotate3D_geographic(point, q)
    c = convert_to_cartesian(point)
    p = Quaternion(c[1], c[2], c[3], 0.0)
    R = conj(q) * p * q
    convert_to_geographic([R[1], R[2], R[3]])
end

"""
rotate3D_cartesian(point, q)

Rotates a point in the 3D space with the given point and the unit quaternion.
"""
function rotate3D_cartesian(point, q)
    p = Quaternion(point..., 0.0)
    R = conj(q) * p * q
    [R[1], R[2], R[3]]
end

"""
get_fiber(point, segments, samples; r=0.025)

Calculates a torus of revolution for building a surface in a specific way with
the given point in the base space, the number of segments, the number of
samples and the radius of the smaller circle in the torus of revolution.
"""
function get_fiber(point, segments, samples; r=0.01)
    # Find 3 points on the circle
    b = λ⁻¹map(convert_to_cartesian(point))
    A = λmap(S¹action(pi / 6, b))
    B = λmap(S¹action(pi / 4, b))
    C = λmap(S¹action(pi / 3, b))
    # The fiber circle center
    Q = get_center(A, B, C)
    # The bigger radius
    R = norm(Q - A)
    # Get the normal to the plane containing the points
    n = LinearAlgebra.cross(A - Q, B - Q)
    n = n / LinearAlgebra.norm(n)
    # The initial normal to the circle
    i = [0.0, 0.0, 1.0]
    # The axis of rotation
    u = LinearAlgebra.cross(n, i)
    u = u / LinearAlgebra.norm(u)
    # The angle of rotation
    angle = acos(LinearAlgebra.dot(n, i)) / 2.0
    q = Quaternion(sin(angle)*u[1],
                   sin(angle)*u[2],
                   sin(angle)*u[3],
                   cos(angle))
    # Construct a torus of revolution grid
    manifold = Array{Float64}(undef, segments, samples, 3)
    for i in 1:segments
        for j in 1:samples
            longitude = i * 2pi / (segments - 1)
            latitude = j * 2pi / (samples - 1)
            x₁ = (Q[1] + (R + r * cos(longitude)) * cos(latitude)) / R
            x₂ = (Q[2] + (R + r * cos(longitude)) * sin(latitude)) / R
            x₃ = (Q[3] + r * sin(longitude)) / R
            manifold[i, j, :] = rotate3D_cartesian([x₁, x₂, x₃], q)
        end
    end
    manifold
end

# The scene object that contains other visual objects
universe = Scene(backgroundcolor = :white, show_axis=false, resolution = (360, 360))
#universe = Scene(backgroundcolor = :white, show_axis=false, resolution = (2560, 1440))

# Calculate a unit quaternion as the rotation axis
u = [sqrt(3)/3, sqrt(3)/3, sqrt(3)/3]
ϕ = Node(0.0)
q = @lift(Quaternion(sin($ϕ)*u[1],
                     sin($ϕ)*u[2],
                     sin($ϕ)*u[3],
                     cos($ϕ)))
segments = 36
samples = 72
number = 360
points = get_flower(number = number)
for i in 1:number
    rotated = @lift(rotate3D_geographic(points[i, :], $q))
    fiber = @lift(get_fiber($rotated, segments, samples))
    color = @lift begin
        x₁, x₂, x₃ = convert_to_cartesian($rotated)
        fill(RGBAf0(rand()/3+2x₁/3, rand()/3+2x₂/3, rand()/3+2x₃/3, 0.9),
             segments,
             samples)
    end
    build_surface(universe, fiber, color, shading = false)
end

# update eye position
eye_position, lookat, upvector = Vec3f0(0.01, 0, 5), Vec3f0(0), Vec3f0(0, 0, 1.0)
update_cam!(universe, eye_position, lookat)
universe.center = false # prevent scene from recentering on display
# Makie.save("gallery/porta.jpg", universe)
frames = 90
record(universe, "gallery/flower.gif") do io
    for i in 1:frames
        ϕ[] = i*2pi/frames # animate scene
        recordframe!(io) # record a new frame
    end
end

ext: jl

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
using Test
using ForestBiometrics
using DelimitedFiles

datapath = joinpath(@__DIR__, "data")
dat=readdlm(joinpath(datapath,"IEsubset_Data_CSV.csv"),',')

#select columns from example data
dbh=dat[2:35,10] #dbh; ignore header
spp=dat[2:35,9] #species; ignore header

# 2 parameter Wyckoff coefficients. Default values for IE variant of FVS
FVS_IE=Dict{String,Array{Float64}}(
"WP"=>[5.19988	-9.26718],
"WL"=>[4.97407	-6.78347],
"DF"=>[4.81519	-7.29306],
"GF"=>[5.00233	-8.19365],
"WH"=>[4.97331	-8.19730],
"RC"=>[4.89564	-8.39057],
"LP"=>[4.62171	-5.32481],
"ES"=>[4.9219	-8.30289],
"AF"=>[4.76537	-7.61062],
"PP"=>[4.9288	-9.32795],
"MH"=>[4.77951	-9.31743],
"WB"=>[4.97407	-6.78347],
"LM"=>[4.19200	-5.16510],
"LL"=>[4.76537	-7.61062],
"PI"=>[3.20000	-5.00000],
"RM"=>[3.20000	-5.00000],
"PY"=>[4.19200	-5.16510],
"AS"=>[4.44210	-6.54050],
"CO"=>[4.44210	-6.54050],
"MM"=>[4.44210	-6.54050],
"PB"=>[4.44210	-6.54050],
"OH"=>[4.44210	-6.54050],
"OS"=>[4.77951	-9.31743] )

wyckoff_test=[
100.3704728
61.35043145
11.31948442
11.99184897
71.56773647
75.48849752
4.662859967
38.1318738
86.57926895
89.09331657
83.29702964
82.02991682
59.36359424
92.24130804
80.79924609
95.99589413
89.52703552
73.04225223
95.88342114
59.38502782
4.662859967
4.662859967
97.69440543
96.76176795
106.4984871
75.82474076
100.3704728
105.4603241
92.79168745
78.84582366
79.6182448
77.22661356
92.79168745
86.208065]

#Wyckoff=(x,b)->4.5+exp(b[1]+(b[2]/(x+1)))
p=HeightDiameter(Wyckoff,FVS_IE)
wyckoff_out= calculate_height(p,dbh,spp)

@test isapprox(wyckoff_out, wyckoff_test)

#end test 1

user_eq_test=[
5.50000110437066
5.50000077727861
5.50000025616043
5.50000027222188
5.50000085504705
5.50000088557174
5.49997579696978
5.50000059480051
5.50000097615539
5.50000099792679
5.50000094851992
5.50002825335774
5.50002185801763
5.50003169790446
5.50000092800344
5.50003314012292
5.50000100174185
5.50000086646161
5.50003309510479
5.50000076242647
5.49997579696978
5.49997579696978
5.50003383536612
5.50003344996446
5.50000117082406
5.50000088821741
5.50000110437066
5.50000115899896
5.50003190202163
5.50002728005371
5.50002751279364
5.50002679851972
5.50003190202163
5.50002959423972]

p2=HeightDiameter((x,b)->4.5+x^b[1]^b[2],FVS_IE)

user_eq_out = calculate_height(p2,dbh,spp)

@test isapprox(user_eq_out, user_eq_test)
#end test 2


text:
import GridGames
using PrettyTables
using Dates

mutable struct Counter
    n::Int32
end
Base.:(==)(a::Counter, b::Counter) = a.n == b.n
set!(c::Counter, n) = c.n = n
increment!(c::Counter) = c.n += 1

mutable struct RunningAvg
    count::Int
    avg::Float32
    RunningAvg() = new(0, 0.)
end

function update_avg!(x::RunningAvg, amount)
    x.count += 1
    x.avg += (1 ./ x.count)*(amount - x.avg)
end

# @todo rename to MuZeroNetwork
abstract type Network end

abstract type AbstractAction end
abstract type AbstractGame end
abstract type AbstractState end
abstract type AbstractPlayer end

AbstractEnv = Union{AbstractGame,GridGames.AbstractGridGame}
AbstractEnvState = Union{AbstractState,GridGames.AbstractGridState}
AbstractEnvAction = Union{AbstractAction,GridGames.AbstractGridAction}
AbstractEnvPlayer = Union{AbstractPlayer,GridGames.AbstractGridPlayer}

abstract type AbstractBoardGame <: AbstractGame end
abstract type AbstractBoardState <: AbstractState end
abstract type AbstractBoardGameSVG <: AbstractBoardGame end
abstract type AbstractBoardStateSVG <: AbstractBoardState end
abstract type TrainingTest end
abstract type TrainingTests end
abstract type AbstractImage end
abstract type AbstractMuZeroGame end
abstract type AbstractLearner end

struct BoardGameSVG <: AbstractBoardGameSVG
    x::AbstractBoardGame
end

avgx(x,size) = mean.(collect(Iterators.partition(x, size)))

mutable struct MuZeroGame <: AbstractMuZeroGame
    game::AbstractEnv
    child_visits::Vector{Vector{Float16}}
    root_values::Vector{Float16}
    net_step::Integer
    function MuZeroGame(game::AbstractEnv, net_step=0)
        new(
            game,
            Vector{Int16}(),
            [],
            net_step,
        )
    end
end

# Sample position from game either uniformly or according to some priority.
function sample_position(game::MuZeroGame)
    to = length(history(game)) + 1
    rand(1:to)
end

ensuredirexists(dir) = mkpath(expanduser(dir))

function ensure_pathto_exists(path)
    a = split(expanduser(path), "/")
    join(a[1:(end-1)], "/") |> mkpath
end

valornothing(d, k) = valornothing(d, (x)->x, k)
valornothing(d, fn, k) = haskey(d, k) ? fn(d[k]) : nothing
valormissing(d, k) = valormissing(d, (x)->x, k)
valormissing(d, fn, k) = haskey(d, k) ? fn(d[k]) : missing

lastbykeyornothing(d::Dict) = lastbykeyornothing(d,(x)->x)
function lastbykeyornothing(d::Dict, fn)
    if isempty(d)
        return nothing
    end
    last_key = keys(d) |> collect |> sort |> last
    fn(d[last_key])
end

isfileexpanduser(path) = isfile(expanduser(path))
isdirexpanduser(path) = isdir(expanduser(path))
readdirexpanduser(path) = readdir(expanduser(path), join=true)
serializeexpanduser(path, x) = serialize(expanduser(path), x)
deserializeexpanduser(path) = deserialize(expanduser(path))

getcheckpointstr(dt::DateTime) = getcheckpointstr(string(dt))

function getcheckpointstr(datetime_str::String)
    replace(
        datetime_str,
        r"[\:]"=>"-"
    )[1:(end-4)]
end


text:
#adapted from https://stackoverflow.com/questions/51038294/resize-image-using-nearest-neighborhood-with-cuda

function _CUDA_resize(pIn,pOut)
    _CUDA_resize(pIn,pOut,size(pIn,1),size(pIn,2),size(pOut,1),size(pOut,2))
end

function _CUDA_resize(pIn,pOut,w_in,h_in,w_out,h_out)

    index_i = blockDim().y * (blockIdx().y-1) + threadIdx().y
    index_j = blockDim().x * (blockIdx().x-1) + threadIdx().x

    stride_i = blockDim().y * gridDim().y
    stride_j = blockDim().x * gridDim().x

    for j=index_j:stride_j:h_out
        for i=index_i:stride_i:w_out
            jIn = div(j*h_in, h_out)
            iIn = div(i*w_in, w_out)
            @inbounds pOut[i,j] = pIn[iIn,jIn]
        end
    end
    return
end

function CUDA_resize(pIn,pOut)
    CUDA.@sync @cuda threads=(16,16) _CUDA_resize(pIn,pOut)
end

function _CUDA_resize4(pIn,pOut,w_in,h_in,w_out,h_out,n)

    c = blockIdx().y
    b = blockIdx().x
    a = threadIdx().x

    if ((a <= h_out)&&(b <= w_out))&&(c <=n)
        jIn = div(b*h_in, h_out)
        iIn = div(a*w_in, w_out)
        @inbounds pOut[a,b,1,c] = pIn[iIn,jIn,1,c]
    end

    return
end

function _CUDA_resize4(pIn,pOut)
    _CUDA_resize4(pIn,pOut,size(pIn,1),size(pIn,2),size(pOut,1),size(pOut,2),size(pIn,4))
end

function CUDA_resize4(pIn,pOut)
    numblocks_x = size(pOut,2)
    numblocks_y = size(pOut,4)

    CUDA.@sync @cuda threads=256 blocks=(numblocks_x,numblocks_y) _CUDA_resize4(pIn,pOut)
end

function _CUDA_normalize_images(pIn,meanImg,h_out,w_out,n)
    index_i = blockDim().y * (blockIdx().y-1) + threadIdx().y
    index_j = blockDim().x * (blockIdx().x-1) + threadIdx().x

    stride_i = blockDim().y * gridDim().y
    stride_j = blockDim().x * gridDim().x

    for k=1:n
        for j=index_j:stride_j:h_out
            for i=index_i:stride_i:w_out
                @inbounds pIn[i,j,1,k] = pIn[i,j,1,k] / 255
                @inbounds pIn[i,j,1,k] = pIn[i,j,1,k] - meanImg[i,j]
            end
        end
    end
    return
end

function _CUDA_normalize_images(pIn,meanImg)
    _CUDA_normalize_images(pIn,meanImg,size(pIn,1),size(pIn,2),size(pIn,4))
end

function CUDA_normalize_images(pIn,meanImg)
    CUDA.@sync @cuda threads=256 _CUDA_normalize_images(pIn,meanImg)
end

function CUDA_preprocess(pIn,pOut)
    n = size(pIn,4)

    w_in=size(pIn,1)
    h_in=size(pIn,2)
    w_out=size(pOut,1)
    h_out=size(pOut,2)

    CUDA.@sync @cuda threads=(16,16) _CUDA_preprocess(pIn,pOut,w_in,h_in,w_out,h_out,n)
end

function _CUDA_preprocess(pIn,pOut,w_in,h_in,w_out,h_out,n)
    index_i = blockDim().y * (blockIdx().y-1) + threadIdx().y
    index_j = blockDim().x * (blockIdx().x-1) + threadIdx().x

    stride_i = blockDim().y * gridDim().y
    stride_j = blockDim().x * gridDim().x

    for k=1:n
        for j=index_j:stride_j:h_out
            for i=index_i:stride_i:w_out
                jIn = div(j*h_in, h_out)
                iIn = div(i*w_in, w_out)
                @inbounds pOut[j,i,1,k] = Float32(pIn[iIn,jIn,1,k])
            end
        end
    end
    return
end

function _CUDA_blur_x(pIn,pOut,w_in,h_in,gauss,n)

    c = blockIdx().z
    b = blockIdx().y
    a = threadIdx().x

    while (a<=w_in)
        if ((a <= w_in)&&(b <= h_in))&&(c <=n)
            temp = 0.0
            for yy = 1:9

                new_y = yy-5 + a
                if new_y<1
                    new_y=abs(new_y)+1
                end
                if new_y>w_in
                    new_y = w_in - (new_y - w_in)
                end

                @inbounds temp = temp + Float32(pIn[new_y,b,1,c]) * gauss[yy]
            end
            @inbounds pOut[a,b,1,c] = temp
        end
        a += blockDim().x*blockIdx().x
    end
    return nothing
end

function _CUDA_blur_x(pIn,pOut,gauss)
    _CUDA_blur_x(pIn,pOut,size(pIn,1),size(pIn,2),gauss,size(pIn,4))
end

function CUDA_blur_x(pIn,pOut,gauss)
    numblocks_x = ceil(Int,size(pOut,1)/ 256)
    numblocks_y = size(pOut,2)
    numblocks_z = size(pOut,4)
    CUDA.@sync @cuda threads=256 blocks=(numblocks_x,numblocks_y,numblocks_z) _CUDA_blur_x(pIn,pOut,gauss)
end

function _CUDA_blur_y(pIn,pOut,gauss)
    _CUDA_blur_y(pIn,pOut,size(pIn,1),size(pIn,2),gauss,size(pIn,4))
end

function CUDA_blur_y(pIn,pOut,gauss)
    numblocks_x = ceil(Int,size(pOut,2)/ 256)
    numblocks_y = size(pOut,1)
    numblocks_z = size(pOut,4)
    CUDA.@sync @cuda threads=256 blocks=(numblocks_x,numblocks_y,numblocks_z) _CUDA_blur_y(pIn,pOut,gauss)
end

function _CUDA_blur_y(pIn,pOut,w_in,h_in,gauss,n)

    c = blockIdx().z
    b = blockIdx().y
    a = threadIdx().x

    while (a<=h_in)
        if ((a <= h_in)&&(b <= w_in))&&(c <=n)
            temp = 0.0
            for yy = 1:9

                new_y = yy-5 + a
                if new_y<1
                    new_y=abs(new_y)+1
                end
                if new_y>h_in
                    new_y = h_in - (new_y - h_in)
                end

                @inbounds temp = temp + Float32(pIn[b,new_y,1,c]) * gauss[yy]
            end
            @inbounds pOut[b,a,1,c] = temp
        end
        a += blockDim().x*blockIdx().x
    end
    return nothing
end

function _CUDA_flip_xy(pIn,pOut,w_in,h_in,n)

    c = blockIdx().y
    b = blockIdx().x
    a = threadIdx().x

    if ((a <= h_in)&&(b <= w_in))&&(c <=n)
        @inbounds pOut[a,b,1,c] = pIn[b,a,1,c]
    end

    return
end

function _CUDA_flip_xy(pIn,pOut)
    _CUDA_flip_xy(pIn,pOut,size(pIn,1),size(pIn,2),size(pIn,4))
end

function CUDA_flip_xy(pIn,pOut)
    numblocks_x = size(pOut,2)
    numblocks_y = size(pOut,4)

    CUDA.@sync @cuda threads=256 blocks=(numblocks_x,numblocks_y) _CUDA_flip_xy(pIn,pOut)
end


text:
using
  CUDA,
  FFTW,
  LinearAlgebra,
  Printf,
  JLD2,
  Documenter,
  Test

using
  FourierFlows,
  FourierFlows.Diffusion

using FourierFlows: parsevalsum2

using LinearAlgebra: mul!, ldiv!, norm

# the devices on which tests will run
devices = CUDA.functional() ? (CPU(), GPU()) : (CPU(),)

const rtol_fft = 1e-12
const rtol_output = 1e-12
const rtol_utils = 1e-13
const rtol_timesteppers = 1e-12

const steppers = [
  "ForwardEuler",
  "RK4",
  "ETDRK4",
  "AB3",
  "FilteredForwardEuler",
  "FilteredRK4",
  "FilteredETDRK4",
  "FilteredAB3"
]

# Run tests
include("createffttestfunctions.jl")

for dev in devices

  @info "testing on " * string(typeof(dev))

  @time @testset "Grid tests" begin
    include("test_grid.jl")

    nx, Lx =  6, 2π
    ny, Ly =  8, 4π
    nz, Lz = 10, 3.0

    # Test 1D grid
    g₁ = OneDGrid(dev, nx, Lx)
    @test testnx(g₁, nx)
    @test testdx(dev, g₁)
    @test testdk(g₁)
    @test testx(g₁)
    @test testk(g₁)
    @test testkr(g₁)
    @test testdealias(g₁)
    @test testmakefilter(dev, g₁)

    # Test 2D rectangular grid
    g₂ = TwoDGrid(dev, nx, Lx, ny, Ly)
    @test testnx(g₂, nx)
    @test testny(g₂, ny)
    @test testdx(dev, g₂)
    @test testdy(dev, g₂)
    @test testdk(g₂)
    @test testdl(g₂)
    @test testx(g₂)
    @test testy(g₂)
    @test testk(g₂)
    @test testkr(g₂)
    @test testl(g₂)
    @test testgridpoints(dev, g₂)
    @test testdealias(g₂)
    @test testmakefilter(dev, g₂)

    # Test 3D parallelogram grid
    g₃ = ThreeDGrid(dev, nx, Lx, ny, Ly, nz, Lz)
    @test testnx(g₃, nx)
    @test testny(g₃, ny)
    @test testnz(g₃, nz)
    @test testdx(dev, g₃)
    @test testdy(dev, g₃)
    @test testdz(dev, g₃)
    @test testdk(g₃)
    @test testdl(g₃)
    @test testdm(g₃)
    @test testx(g₃)
    @test testy(g₃)
    @test testz(g₃)
    @test testk(g₃)
    @test testkr(g₃)
    @test testl(g₃)
    @test testm(g₃)
    @test testgridpoints(dev, g₃)
    @test testdealias(g₃)
    @test testmakefilter(dev, g₃)
    @test test_plan_flows_fftrfft(dev)

    # Test typed grids
    T = Float32
    @test test_plan_flows_fftrfft(dev, T=T)
    @test testtypedonedgrid(dev, nx, Lx; T=T)
    @test testtypedtwodgrid(dev, nx, Lx, ny, Ly; T=T)
    @test testtypedthreedgrid(dev, nx, Lx, ny, Ly, nz, Lz; T=T)

    # Test aliased fraction
    for aliased_fraction ∈ [0, 1/3, 1/2, 1/4]
      @test test_aliased_fraction(dev, aliased_fraction)
    end

    # Test show() methods
    @test repr(g₁) == "OneDimensionalGrid\n  ├─────────── Device: "*FourierFlows.griddevice(g₁)*"\n  ├──────── FloatType: Float64\n  ├────────── size Lx: 6.283185307179586\n  ├──── resolution nx: 6\n  ├── grid spacing dx: 1.0471975511965976\n  ├─────────── domain: x ∈ [-3.141592653589793, 2.094395102393195]\n  └─ aliased fraction: 0.3333333333333333"

    @test repr(g₂) == "TwoDimensionalGrid\n  ├───────────────── Device: "*FourierFlows.griddevice(g₂)*"\n  ├────────────── FloatType: Float64\n  ├────────── size (Lx, Ly): (6.283185307179586, 12.566370614359172)\n  ├──── resolution (nx, ny): (6, 8)\n  ├── grid spacing (dx, dy): (1.0471975511965976, 1.5707963267948966)\n  ├───────────────── domain: x ∈ [-3.141592653589793, 2.094395102393195]\n  |                          y ∈ [-6.283185307179586, 4.71238898038469]\n  └─ aliased fraction: 0.3333333333333333"

    @test repr(g₃) == "ThreeDimensionalGrid\n  ├───────────────────── Device: "*FourierFlows.griddevice(g₃)*"\n  ├────────────────── FloatType: Float64\n  ├────────── size (Lx, Ly, Lz): (6.283185307179586, 12.566370614359172, 12.566370614359172)\n  ├──── resolution (nx, ny, nz): (6, 8, 10)\n  ├── grid spacing (dx, dy, dz): (1.0471975511965976, 1.5707963267948966, 0.3)\n  ├────────────────────  domain: x ∈ [-3.141592653589793, 2.094395102393195]\n  |                              y ∈ [-6.283185307179586, 4.71238898038469]\n  |                              z ∈ [-1.5, 1.2]\n  └─ aliased fraction: 0.3333333333333333"

    # Test no dealiasing
    g₁ = OneDGrid(nx, Lx; aliased_fraction = 0)
    g₂ = TwoDGrid(nx, Lx, ny, Ly; aliased_fraction = 0)
    g₃ = ThreeDGrid(nx, Lx, ny, Ly, nz, Lz; aliased_fraction = 0)

    @test testnodealias(g₁)
    @test testnodealias(g₂)
    @test testnodealias(g₃)
  end

  @time @testset "FFT tests" begin
    include("test_fft.jl")

    # Test 1D grid
    nx = 32             # number of points
    Lx = 2π             # Domain width
    g1 = OneDGrid(dev, nx, Lx)

    @test test_fft_cosmx(g1)
    @test test_rfft_cosmx(g1)
    @test test_rfft_mul_cosmx(g1)

    # Test 2D rectangular grid
    nx, ny = 32, 64     # number of points
    Lx, Ly = 2π, 3π     # Domain width
    g2 = TwoDGrid(dev, nx, Lx, ny, Ly)

    @test test_fft_cosmxcosny(g2)
    @test test_rfft_cosmxcosny(g2)
    @test test_rfft_mul_cosmxcosny(g2)
    @test test_fft_sinmxny(g2)
    @test test_rfft_sinmxny(g2)
    @test test_rfft_mul_sinmxny(g2)

    # Test 3D parallelogram grid
    nx, ny, nz = 32, 30, 16     # number of points
    Lx, Ly, Lz = 2π, 3π, 4.0    # Domain width
    g3 = ThreeDGrid(dev, nx, Lx, ny, Ly, nz, Lz)

    @test test_fft_cosmxcosny(g3)
    @test test_rfft_cosmxcosny(g3)
    @test test_rfft_mul_cosmxcosny(g3)
    @test test_fft_sinmxny(g3)
    @test test_rfft_sinmxny(g3)
    @test test_rfft_mul_sinmxny(g3)
  end

  @time @testset "IFFT tests" begin
    include("test_ifft.jl")

    # Test 1D grid
    nx = 32             # number of points
    Lx = 2π             # Domain width
    g1 = OneDGrid(dev, nx, Lx)

    @test test_ifft_cosmx(g1)
    @test test_irfft_cosmx(g1)
    @test test_irfft_mul_cosmx(g1)

    # Test 2D rectangular grid
    nx, ny = 32, 64     # number of points
    Lx, Ly = 2π, 3π     # Domain width
    g2 = TwoDGrid(dev, nx, Lx, ny, Ly)

    @test test_ifft_cosmxcosny(g2)
    @test test_irfft_cosmxcosny(g2)
    @test test_irfft_mul_cosmxcosny(g2)
    @test test_ifft_sinmxny(g2)
    @test test_irfft_sinmxny(g2)
    @test test_irfft_mul_sinmxny(g2)

    # Test 3D parallelogram grid
    nx, ny, nz = 32, 30, 16     # number of points
    Lx, Ly, Lz = 2π, 3π, 4.0    # Domain width
    g3 = ThreeDGrid(dev, nx, Lx, ny, Ly, nz, Lz)

    @test test_ifft_cosmxcosny(g3)
    @test test_irfft_cosmxcosny(g3)
    @test test_irfft_mul_cosmxcosny(g3)
    @test test_ifft_sinmxny(g3)
    @test test_irfft_sinmxny(g3)
    @test test_irfft_mul_sinmxny(g3)
  end

  @time @testset "Timestepper tests" begin
    include("test_timesteppers.jl")
    for stepper in steppers
      @test constantdiffusiontest_stepforward(stepper, dev=dev)

      @test varyingdiffusiontest_stepforward(stepper, dev=dev)

      if FourierFlows.isexplicit(stepper)
        @test constantdiffusiontest_step_until(stepper, dev=dev)
      else
        @test_throws Exception constantdiffusiontest_step_until(stepper, dev=dev)
      end
    end
  end

  @time @testset "Utils tests" begin
    include("test_utils.jl")

    @test test_fltype()
    @test test_cxtype()
    @test test_innereltype()
    @test test_superzeros()
    @test test_supertuplezeros()
    @test test_supersize()
    @test test_arraytype(dev)
    @test test_arraytypeTdim(dev, Float32, 2)

    # Test on a rectangular grid
    nx, ny = 64, 128   # number of points
    Lx, Ly = 2π, 3π    # Domain width
    g = TwoDGrid(dev, nx, Lx, ny, Ly)
    x, y = gridpoints(g)
    k0, l0 = 2π/Lx, 2π/Ly

    # Real and complex-valued functions
    σ = 0.5
    f1 = @. exp(-(x^2 + y^2)/(2σ^2))
    f2 = @. (cos(2k0*x + 3l0*y^2) + im*sin(2k0*x + 3l0*y^2)) * (exp(-(x^2 + y^2)/(2σ^2)) + 2im*exp(-(x^2 + y^2)/(5σ^2)))

    # Sine/Exp waves
    k1, l1 = 2*k0,  6*l0
    k2, l2 = 3*k0, -3*l0
    sinkl1 = @. sin(k1*x + l1*y)
    sinkl2 = @. sin(k2*x + l2*y)
    expkl1 = @. cos(k1*x + l1*y) + im*sin(k1*x + l1*y)
    expkl2 = @. cos(k2*x + l2*y) + im*sin(k2*x + l2*y)

    # Analytical expression for the Jacobian of sin1 and sin2 and of exp1 and exp2
    Jsinkl1sinkl2 = @. (k1*l2-k2*l1)*cos(k1*x + l1*y)*cos(k2*x + l2*y)
    Jexpkl1expkl2 = @. (k2*l1-k1*l2)*(cos((k1+k2)*x + (l1+l2)*y)+im*sin((k1+k2)*x + (l1+l2)*y))

    @test test_parsevalsum(f1, g; realvalued=true)   # Real valued f with rfft
    @test test_parsevalsum(f1, g; realvalued=false)  # Real valued f with fft
    @test test_parsevalsum(f2, g; realvalued=false)  # Complex valued f with fft
    @test test_parsevalsum2(f1, g; realvalued=true)  # Real valued f with rfft
    @test test_parsevalsum2(f1, g; realvalued=false) # Real valued f with fft
    @test test_parsevalsum2(f2, g; realvalued=false) # Complex valued f with fft

    @test test_jacobian(sinkl1, sinkl1, 0*sinkl1, g)      # Test J(a, a) = 0
    @test test_jacobian(sinkl1, sinkl2, Jsinkl1sinkl2, g) # Test J(sin1, sin2) = Jsin1sin2
    @test test_jacobian(expkl1, expkl2, Jexpkl1expkl2, g) # Test J(exp1, exp2) = Jexp1exps2

    @test test_zeros()

    g = OneDGrid(dev, nx, Lx)
    σ = 0.5
    f1 = @. exp(-g.x^2 / 2σ^2)
    @test test_parsevalsum2(f1, g; realvalued=true)  # Real valued f with rfft
    @test test_parsevalsum2(f1, g; realvalued=false) # Real valued f with fft

    # Radial spectrum tests. Note that ahρ = ∫ ah ρ dθ.
    n = 128; δ = n/10                                       # Parameters
    ahkl_isotropic(k, l) = exp(-(k^2 + l^2) / 2δ^2)                     # a  = exp(-ρ²/2δ²)
    ahρ_isotropic(ρ) = 2π * ρ * exp(-ρ^2 / 2δ^2)                        # aᵣ = 2π ρ exp(-ρ²/2δ²)
    @test test_radialspectrum(dev, n, ahkl_isotropic, ahρ_isotropic)
    @test test_radialspectrum(dev, n, ahkl_isotropic, ahρ_isotropic; rfft=true)

    ahkl_anisotropic(k, l) = exp(-(k^2 + l^2) / 2δ^2) * k^2 / (k^2+l^2) # a  = exp(-ρ²/2δ²) cos²θ
    ahρ_anisotropic(ρ) = π * ρ * exp(-ρ^2/2δ^2)                         # aᵣ = π ρ exp(-ρ²/2δ²)
    @test test_radialspectrum(dev, n, ahkl_anisotropic, ahρ_anisotropic)
    @test test_radialspectrum(dev, n, ahkl_anisotropic, ahρ_anisotropic; rfft=true)
    @test test_ongrid(dev)
  end

  @time @testset "Problem instantiation tests" begin
    include("test_instantiate_problem.jl")

    for stepper in steppers
      @test instantiate_problem(dev, stepper)
      if occursin("Filtered", stepper)
        @test instantiate_problem_with_filter_kwargs(dev, stepper)
      end
    end
  end

  @time @testset "Diagnostics tests" begin
    include("test_diagnostics.jl")

    @test test_diagnosticsteps(dev, freq=1)
    @test test_diagnosticsteps(dev, freq=2)
    @test test_diagnosticsteps(dev, nsteps=100, freq=9, ndata=20)
    @test test_basicdiagnostics(dev)
    @test test_scalardiagnostics(dev, freq=1)
    @test test_scalardiagnostics(dev, freq=2)
    @test test_incrementdiagnostic(dev)
    @test test_extenddiagnostic(dev)
    @test test_getindex(dev)
  end

  @time @testset "Output tests" begin
    include("test_output.jl")

    @test test_withoutjld2()
    @test test_uniquepath()
    @test test_outputconstructor(dev)
    @test test_getindex(dev)
    @test test_saveproblem_saveoutput(dev)
    @test test_saveproblemTwoDGrid(dev)
    @test test_savediagnostic(dev)

    file = test_savefields(dev; parameter=2.2)
    @test file["params"]["parameter"] == 2.2
    @test_throws KeyError file["params"]["func"] == func
    @test_throws KeyError file["params"]["plan"] == plan
    close(file)
  end

  @time @testset "show() methods tests" begin
    prob = Diffusion.Problem(dev)

    struct Params1 <: AbstractParams
        κ1 :: Float64
        κ2 :: Float64
      func :: Function
    end

    struct Params2 <: AbstractParams
        κ1 :: Float64
      func :: Function
        κ2 :: Float64
    end

    func(x) = sin(x^2)

    params1 = Params1(1.0, 2.0, func)
    params2 = Params2(1.0, func, 2.0)

    prob1 = FourierFlows.Problem(prob.eqn, "RK4", prob.clock.dt, prob.grid, prob.vars, params1, dev)
    prob2 = FourierFlows.Problem(prob.eqn, "RK4", prob.clock.dt, prob.grid, prob.vars, params2, dev)
    prob = prob1 # prob2 is only useful for testing params show() method

    get_sol(prob) = prob.sol
    diag = Diagnostic(get_sol, prob)
    out = Output(prob, "output.jld2")

    @test repr(prob1.params) == "Parameters\n  ├───── parameter: κ1 -> Float64\n  ├───── parameter: κ2 -> Float64\n  └───── parameter: func -> Function\n"
    @test repr(prob2.params) == "Parameters\n  ├───── parameter: κ1 -> Float64\n  ├───── parameter: func -> Function\n  └───── parameter: κ2 -> Float64\n"

    if dev == CPU()
      @test repr(prob.vars) == "Variables\n  ├───── variable: c -> 128-element Vector{Float64}\n  ├───── variable: cx -> 128-element Vector{Float64}\n  ├───── variable: ch -> 65-element Vector{ComplexF64}\n  └───── variable: cxh -> 65-element Vector{ComplexF64}\n"
      @test repr(diag) == "Diagnostic\n  ├─── calc: get_sol\n  ├─── prob: FourierFlows.Problem{DataType, Vector{ComplexF64}, Float64, Vector{Int64}}\n  ├─── data: 101-element Vector{Vector{ComplexF64}}\n  ├────── t: 101-element Vector{Float64}\n  ├── steps: 101-element Vector{Int64}\n  ├─── freq: 1\n  └────── i: 1"
      @test repr(out) == "Output\n  ├──── prob: FourierFlows.Problem{DataType, Vector{ComplexF64}, Float64, Vector{Int64}}\n  ├──── path: output.jld2\n  └── fields: Dict{Symbol, Function}()"
    else
      @test repr(prob.vars) == "Variables\n  ├───── variable: c -> 128-element CuArray{Float64, 1, CUDA.Mem.DeviceBuffer}\n  ├───── variable: cx -> 128-element CuArray{Float64, 1, CUDA.Mem.DeviceBuffer}\n  ├───── variable: ch -> 65-element CuArray{ComplexF64, 1, CUDA.Mem.DeviceBuffer}\n  └───── variable: cxh -> 65-element CuArray{ComplexF64, 1, CUDA.Mem.DeviceBuffer}\n"
      @test repr(diag) == "Diagnostic\n  ├─── calc: get_sol\n  ├─── prob: FourierFlows.Problem{DataType, CuArray{ComplexF64, 1, CUDA.Mem.DeviceBuffer}, Float64, CuArray{Int64, 1, CUDA.Mem.DeviceBuffer}}\n  ├─── data: 101-element Vector{CuArray{ComplexF64, 1, CUDA.Mem.DeviceBuffer}}\n  ├────── t: 101-element Vector{Float64}\n  ├── steps: 101-element Vector{Int64}\n  ├─── freq: 1\n  └────── i: 1"
      @test repr(out) == "Output\n  ├──── prob: FourierFlows.Problem{DataType, CuArray{ComplexF64, 1, CUDA.Mem.DeviceBuffer}, Float64, CuArray{Int64, 1, CUDA.Mem.DeviceBuffer}}\n  ├──── path: output.jld2\n  └── fields: Dict{Symbol, Function}()"
    end

    @test repr(prob.eqn) == "Equation\n  ├──────── linear coefficients: L\n  │                              ├───type: Int64\n  │                              └───size: (65,)\n  ├───────────── nonlinear term: calcN!()\n  └─── type of state vector sol: ComplexF64"
    @test repr(prob.clock) == "Clock\n  ├─── timestep dt: 0.01\n  ├────────── step: 0\n  └──────── time t: 0.0"
    @test repr(prob) == "Problem\n  ├─────────── grid: grid (on " * FourierFlows.griddevice(prob.grid) * ")\n  ├───── parameters: params\n  ├────── variables: vars\n  ├─── state vector: sol\n  ├─────── equation: eqn\n  ├────────── clock: clock\n  └──── timestepper: RK4TimeStepper"
  end

  @time @testset "Doctests" begin
    doctest(FourierFlows)
  end

end # end loop over devices


text:
# Size of buffer
const DIG_BLOCK_SIZE = 1024 * 1024 # 1 MiB

# Clock selection
const CLOCK_INTERNAL = 1
const CLOCK_EXTERNAL = 0

# Channel selection for AD14 and AD16
const IN0 = 1
const IN1 = 2
const IN2 = 4
const IN3 = 8

# Channel selection for AD12
const AIN0     = 1 # Also called DESQ
const AIN1     = 2 # Also called DESI
const DESCLKIQ = 4 # Both inputs must be driven, I and Q are shorted together
const DESIQ    = 8 # Both inputs must be driven, I and Q are separate, increasing input BW

# Triggering options
const NO_TRIGGER               = 0
const WAVEFORM_TRIGGER         = 1
const SYNC_SELECTIVE_RECORDING = 2
const HETERODYNE               = 3
const TTL_TRIGGER_EDGE         = 4

const FALLING_EDGE = 0
const RISING_EDGE  = 1


text:
module DiffinDiffsBase

using Base: @propagate_inbounds
using CSV
using CodecZlib: GzipDecompressorStream
using DataAPI
using DataAPI: refarray, refpool, invrefpool
using Dates: Period, TimeType
using LinearAlgebra: Diagonal
using Missings: allowmissing, disallowmissing
using PooledArrays: _label
using PrettyTables: pretty_table
using Reexport
using StatsBase: CoefTable, Weights, stderror, uweights
using StatsFuns: normccdf, norminvcdf, tdistccdf, tdistinvcdf
@reexport using StatsModels
using StatsModels: Schema
@reexport using StatsProcedures
using StatsProcedures: ≊, _args_kwargs, _parse!
using StructArrays: StructArray
using Tables
using Tables: AbstractColumns, istable, columnnames, getcolumn

import Base: ==, +, -, *, isless, show, parent, view, diff
import Base: eltype, firstindex, lastindex, getindex, iterate, length
import Missings: allowmissing, disallowmissing
import StatsBase: coef, vcov, confint, nobs, dof_residual, responsename, coefnames, weights,
    coeftable
import StatsModels: concrete_term, schema, termvars, lag, lead
import StatsProcedures: required, default, transformed, combinedargs, copyargs, _show_args

# Reexport objects from StatsBase
export coef, vcov, stderror, confint, nobs, dof_residual, responsename, coefnames, weights,
    coeftable

export cb,
       ≊,
       exampledata,

       RotatingTimeValue,
       rotatingtime,
       RotatingTimeArray,

       VecColumnTable,
       VecColsRow,
       subcolumns,
       apply,
       apply_and!,
       apply_and,
       TableIndexedMatrix,

       ScaledArray,
       ScaledVector,
       ScaledMatrix,
       scale,
       align,

       TreatmentSharpness,
       SharpDesign,
       sharp,
       AbstractTreatment,
       DynamicTreatment,
       dynamic,

       ParallelCondition,
       Unconditional,
       unconditional,
       CovariateConditional,
       ParallelStrength,
       Exact,
       exact,
       Approximate,
       AbstractParallel,
       UnspecifiedParallel,
       unspecifiedpr,
       TrendParallel,
       TrendOrUnspecifiedPR,
       NeverTreatedParallel,
       nevertreated,
       NotYetTreatedParallel,
       notyettreated,
       istreated,
       istreated!,

       TermSet,
       termset,
       eachterm,
       TreatmentTerm,
       treat,

       findcell,
       cellrows,
       settime,
       aligntime,
       PanelStructure,
       setpanel,
       findlag!,
       findlead!,
       ilag!,
       ilead!,
       diff!,
       diff,

       CheckData,
       GroupTreatintterms,
       GroupXterms,
       GroupContrasts,
       CheckVars,
       GroupSample,
       MakeWeights,

       DiffinDiffsEstimator,
       DefaultDID,
       did,
       didspec,
       @did,
       AbstractDIDResult,
       DIDResult,
       AggregatedDIDResult,
       vce,
       treatment,
       outcomename,
       treatnames,
       treatcells,
       ntreatcoef,
       treatcoef,
       treatvcov,
       coefinds,
       ncovariate,
       agg,
       SubDIDResult,
       TransformedDIDResult,
       TransSubDIDResult,
       lincom,
       rescale,
       ExportFormat,
       StataPostHDF,
       getexportformat,
       setexportformat!,
       post!

include("tables.jl")
include("utils.jl")
include("time.jl")
include("ScaledArrays.jl")
include("treatments.jl")
include("parallels.jl")
include("terms.jl")
include("operations.jl")
include("procedures.jl")
include("did.jl")

end


text:
# Autogenerated wrapper script for opusfile_jll for i686-w64-mingw32
export libopusfile, libopusurl

using Ogg_jll
using Opus_jll
using OpenSSL_jll
JLLWrappers.@generate_wrapper_header("opusfile")
JLLWrappers.@declare_library_product(libopusfile, "libopusfile-0.dll")
JLLWrappers.@declare_library_product(libopusurl, "libopusurl-0.dll")
function __init__()
    JLLWrappers.@generate_init_header(Ogg_jll, Opus_jll, OpenSSL_jll)
    JLLWrappers.@init_library_product(
        libopusfile,
        "bin\\libopusfile-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libopusurl,
        "bin\\libopusurl-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
const crc3_rotor_R = 0.10204 # m
const crc3_rotor_radii = [10.20, 20.41, 30.61, 40.82, 51.02, 61.22, 66.33, 71.43, 76.53, 81.63, 86.73, 91.84, 96.94, 100.00] .* crc3_rotor_R / 100.0
const crc3_rotor_chord = [0.1195, 0.1340, 0.1607, 0.1920, 0.2142, 0.2200, 0.2190, 0.2100, 0.2039, 0.1911, 0.1758, 0.1581, 0.1272, 0.0982] .* crc3_rotor_R
const crc3_rotor_twist = [30.35612, 34.12497, 28.91561, 21.20027, 15.47515, 12.6852, 11.21356, 10.46311, 10.0537, 9.066763, 7.5527, 7.540607, 8.273884, 8.540983] .* pi/180

const crc3_rotor_x = 65.1925e-3 - 0.28717 * 0.3048/4 # m
const crc3_rotor_y = 242.55e-3/2 # m
const crc3_rotor_rhub = crc3_rotor_radii[1] # m
const crc3_rotor_rtip = crc3_rotor_radii[end] # m
const crc3_rotor_D = 2 * crc3_rotor_R
const crc3_rotor_cant = 10*pi/180 # inboard inclination of the rotor axes

const crc3_rotor_contours = [
    "e212.dat",
    "e212.dat",
    "e212.dat",
    "e212.dat",
    "e212.dat",
    "e212.dat",
    "e212.dat",
    "e212.dat",
    "e212.dat",
    "e212.dat",
    "e212.dat",
    "e212.dat",
    "e212.dat",
    "e212.dat"
]
# const crc3_rotor_contours = [
#     "contour10.20.csv",
#     "contour20.41.csv",
#     "contour30.61.csv",
#     "contour40.82.csv",
#     "contour51.02.csv",
#     "contour61.22.csv",
#     "contour66.33.csv",
#     "contour71.43.csv",
#     "contour76.53.csv",
#     "contour81.63.csv",
#     "contour86.73.csv",
#     "contour91.84.csv",
#     "contour96.94.csv",
#     "contour100.00.csv"
# ]
const crc3_rotor_polars = [
    "contour10.20_polar_5000.0_0.03.csv",
    "contour20.41_polar_5000.0_0.14.csv",
    "contour51.02_polar_5000.0_0.46.csv",
    "contour51.02_polar_5000.0_0.46.csv",
    "contour51.02_polar_5000.0_0.46.csv",
    "contour61.22_polar_5000.0_0.57.csv",
    "contour66.33_polar_5000.0_0.62.csv",
    "contour66.33_polar_5000.0_0.62.csv",
    "contour66.33_polar_5000.0_0.62.csv",
    "contour66.33_polar_5000.0_0.62.csv",
    "contour66.33_polar_5000.0_0.62.csv",
    "contour66.33_polar_5000.0_0.62.csv",
    "contour96.94_polar_5000.0_0.95.csv",
    "contour100.00_polar_5000.0_1.0.csv"
]

const crc3_wing_spacing = 0.254 # distance from lower wing to upper wing
const crc3_wing_b = 20 * 0.0254 # m
const crc3_wing_c = 0.28717 * 0.3048 # m

const crc3_mass = 3 / 2.2 # kg
const crc3_JInertia = 0.0111111 # kg-m^2

max_stars_repo_name: byuflowlab/MultifidelityTrajectories

max_issues_repo_name: byuflowlab/MultifidelityTrajectories

max_forks_repo_name: byuflowlab/MultifidelityTrajectories


text:
using Animated_CMA_ES, Plots, Revise
function banana(a,b)
  x->(a-x[1])^2+b*(x[2]-x[1]^2)^2
end
f = banana(1.0,1.0)
#f(x)=x[1]^2+x[2]^2

res=cma_es_animated(f,[0.0,6.0],xrange=[-2., 3.], yrange=[-3.,8.],tol=1e-3)
gif(res[3], "cma_es_banana.gif", fps=2)

res=cma_es_animated(f,[0.0,6.0],xrange=[-2., 3.], yrange=[-3.,8.],tol=1e-3)

g(x)=x[1]^2+x[2]^2
res=cma_es_animated(g,[-3.0,-3.0],xrange=[-5., 3.], yrange=[-5.,3.],tol=1e-3,maxiter=20)

gif(res[3], "cma_es_square.gif", fps=2)


text:
module color

import Term: NAMED_COLORS, nospaces, COLORS_16b, remove_brackets, ANSICode, int, CODES, CODES_16BIT_COLORS

export NamedColor, BitColor, RGBColor, get_color

# ----------------------------- types definition ----------------------------- #
"""
    AbstractColor

Abstract color type.
"""
abstract type AbstractColor end

Base.show(io::IO, color::AbstractColor) = print(io, "$(typeof(color))")

struct NamedColor <: AbstractColor
    color::String
end

struct BitColor <: AbstractColor
    color::String
end

struct RGBColor <: AbstractColor
    r::Int
    g::Int
    b::Int
end

function RGBColor(s)
    to_number(x) = '.' ∈ x ? parse(Float64, x) : parse(Int, x)
    r, g, b = to_number.(match(RGB_REGEX, s).captures)
    if r < 1 || g < 1 || b < 1
        r *= 255
        g *= 255
        b *= 255
    end
    return RGBColor(int(r), int(g), int(b))
end

"""
    ANSICode(color; bg::Bool=false)

Create ANSI tags for colors.
"""
function ANSICode(color::NamedColor; bg::Bool = false)
    Δ = bg ? 40 : 30
    v = CODES[color.color]
    return ANSICode("\e[$(Δ + v)m", "\e[$(Δ+9)m")
end

function ANSICode(color::BitColor; bg::Bool = false)
    Δ = bg ? 48 : 38
    v = CODES_16BIT_COLORS[color.color]
    return ANSICode("\e[$Δ;5;$(v)m", "\e[$(Δ+1)m")
end

function ANSICode(color::RGBColor; bg::Bool = false)
        Δ = bg ? 48 : 38
        rgb = "$(color.r);$(color.g);$(color.b)"
        return ANSICode("\e[$Δ;2;$(rgb)m", "\e[$(Δ+1)m")
end

# --------------------------------- is color? -------------------------------- #

RGB_REGEX = r"\(\s*([\d\.]{1,3})\s*,\s*([\d\.]{1,3})\s*,\s*([\d\.]{1,3})\s*\)"
HEX_REGEX = r"#(?:[0-9a-fA-F]{3}){1,2}$"

"""
    is_named_color(string::String)::Bool

Check if a string represents a named color.
"""
is_named_color(string)::Bool = string ∈ NAMED_COLORS || string ∈ COLORS_16b

"""
    is_rgb_color(string::String)::Bool

Check if a string represents a RGB color.
"""
is_rgb_color(string)::Bool = !occursin("on_", string) && occursin(RGB_REGEX, string)

"""
    is_hex_color(string::String)::Bool

Check if a string represents a hex color.
"""
is_hex_color(string)::Bool = !occursin("on_", string) && occursin(HEX_REGEX, string)

"""
    is_color(string::String)::Bool

Check if a string represents color information, of any type.
"""
is_color(string)::Bool = is_named_color(string) || is_rgb_color(string) || is_hex_color(string)

"""
    is_background(string::String)::Bool

Check if a string represents background color information, of any type.
"""
function is_background(string)::Bool
    stripped = nospaces(string)
    length(stripped) < 3 && return false
    return stripped[1:3] == "on_" && is_color(stripped[4:end])
end

# --------------------------------- get color -------------------------------- #
"""
    hex2rgb(hex::String)

Converts a string hex color code to a RGB color
"""
function hex2rgb(hex)::RGBColor
    to_int(h) = parse(Int, h; base = 16)
    r, g, b = [to_int(hex[i:(i + 1)]) for i in (2, 4, 6)]
    return RGBColor(r, g, b)
end

"""
    get_color(string::String; bg=false)::AbstractColor

Extract a color type from a string with color information.
"""
function get_color(string; bg = false)::AbstractColor
    if bg
        string = nospaces(string)[4:end]
    end

    if is_named_color(string)
        if string ∈ COLORS_16b
            idx = findfirst((c) -> c == string, COLORS_16b)[1]
            return BitColor(COLORS_16b[idx])
        else
            idx = findfirst((c) -> c == string, NAMED_COLORS)[1]
            return NamedColor(NAMED_COLORS[idx])
        end
    elseif is_rgb_color(string)
        return RGBColor(string)
    else
        # convert hex to rgb
        return hex2rgb(string)
    end
end

"""
Convert HSL values to RGB values, return a markup color string.
"""
function hsl2rgb(h, s, l)
    C = (1 - abs(2*l - 1)) * s
    X = C * (1 - abs(mod(h/60, 2) - 1))
    M = l - C/2

    if 0 ≤ h < 60
        r, g, b = C, X, 0
    elseif 60 ≤ h < 120
        r, g, b = X, C, 0
    elseif 120 ≤ h < 180
        r, g, b = 0, C, X
    elseif 180 ≤ h < 240
        r, g, b = 0, X, C
    elseif 240 ≤ h < 300
        r, g, b = X, 0, C
    elseif 300 ≤ h ≤ 360
        r, g, b = C, 0, X
    end

    r = (Int ∘ round)((r+M) * 255)
    g = (Int ∘ round)((g+M) * 255)
    b = (Int ∘ round)((b+M) * 255)

    return "($r, $g, $b)"
end
end


text:
"""
    writelines(path::String, rows; nworkers = 1, mode = "w")

Write `rows` to JSONLines file `path`

* `path`: Path to output file
* `rows`: `Tables.jl` compatible data
* Keyword Arguments:
    * `nworkers=1`: Number of threads to use for parsing to JSONLines
    * `mode="w"`: Mode the file is opened in. [See I/O and Network](https://docs.julialang.org/en/v1/base/io-network/)
"""
function writelines(path::String, rows; nworkers = 1, mode = "w")
    if mode ∉ ["a", "a+", "w", "w+", "r+"]
        throw(ArgumentError("Cannot open file in mode: $mode"))
    end
    if mode == "r+" && !isfile(path)
        throw(ArgumentError("File $path does not exist"))
    end
    if !ispath(dirname(abspath(path)))
        throw(ArgumentError("""opening file "$path": No such file or directory"""))
    end
    if nworkers > 1
        rows = Tables.rowtable(rows)
        _twriterows(path, rows, nworkers, mode)
    else
        rows = Tables.namedtupleiterator(rows)
        _writerows(path, rows, mode)
    end
end

function _twriterows(path::String, rows, nworkers, mode)
    len = length(rows)
    parts = partitionrows(len, nworkers)
    bufs = [IOBuffer() for _ in 1:nworkers]
    @sync for i in 1:nworkers
        @spawn begin
            crows = @inbounds(rows[parts[i]])
            for row in crows
                JSON3.write(bufs[i], row)
                write(bufs[i], '\n')
            end
        end
    end
    open(path, mode) do io
        for buf in bufs
            write(io, take!(buf))
        end
    end
end

function _writerows(path, rows, mode)
    open(path, mode) do io
        for row in rows
            JSON3.write(io, row)
            write(io, '\n')
        end
    end
end


text:
#=
An implementation of the Laplacians and SDDM solvers of Koutis, Miller and Peng
=#

#=
global KMP_SAVEMATS=false
global KMP_MATS=[]
global KMP_FS=[]
=#

mutable struct IJVS
    i::Array{Int64,1}
    j::Array{Int64,1}
    v::Array{Float64,1}  # wt of edge
    s::Array{Float64,1}  # stretch of edge
end

"""Parameters for the KMP solver"""
struct KMPParams
    frac::Float64  # fraction to decrease at each level
    iters::Int64   # iters of PCG to apply between levels
    treeScale::Float64 # scale tree by treeScale*log_2 (n) * aveStretch
    n0::Int64 # the number of edges at which to go direct
    treeAlg # :akpw or :rand
end

defaultKMPParams = KMPParams(1/36, 6, 0.125, 600, :akpw)

# this is just for Laplacians, not general SDD
struct elimLeafNode
    nodeid::Int64
    parent::Int64
    wtDeg::Float64
end

# this is just for Laplacians, not general SDD
struct elimDeg2Node
    nodeid::Int64
    nbr1::Int64
    nbr2::Int64
    wt1::Float64
    wt2::Float64
end

# The tree must be in DFS order
# marked is 1 if flagged for possible elimination,
# and set to 2 if we do eliminate it
# is just for Laplacians matrices, not general SDD
function elimDeg12(t, marked)

    # make sure root is not marked
    marked[1] = 0

    n = size(t,1)

    deg = Vector{Int64}(undef, n)
    for v in 1:n
        deg[v] = t.colptr[v+1] - t.colptr[v]
    end

    elims1 = elimLeafNode[]

    for v in n:-1:2

        if (deg[v] == 1 && marked[v] == 1)
            parent = t.rowval[t.colptr[v]];
            wt = t.nzval[t.colptr[v]];
            push!(elims1,elimLeafNode(v,parent,wt))

            deg[parent] = deg[parent] - 1
            marked[v] = 2
            deg[v] = 0
        end
    end

    elims2 = elimDeg2Node[]

    subt = triu(t)

    for v in n:-1:2

        if (deg[v] == 2 && marked[v] == 1)

            parent = t.rowval[t.colptr[v]];

            # to ident the child, enumerate to find one uneliminated, which we check by marked
            kidind = t.colptr[v]+1
            kid = t.rowval[kidind]
            while deg[kid] == 0
                kidind = kidind+1
                kid = t.rowval[kidind]

                if kidind >= t.colptr[v+1]
                    error("went of the end of the kid list without finding node to elim from")
                end
            end

            wt1 = t.nzval[t.colptr[v]];
            wt2 = t.nzval[kidind]

            push!(elims2,elimDeg2Node(v,parent,kid,wt1,wt2))
            marked[v] = 2

            newwt = 1/(1/wt1 + 1/wt2)

            # now that we've found the kid, go up the chain until done
            while (deg[parent] == 2 && marked[parent] == 1)
                v = parent
                parent = t.rowval[t.colptr[v]];
                wt1 = t.nzval[t.colptr[v]];
                wt2 = newwt

                push!(elims2,elimDeg2Node(v,parent,kid,wt1,wt2))
                marked[v] = 2

                newwt = 1/(1/wt1 + 1/wt2)
            end

            # now, hack the tree to adjust parent and wt of node kid
            subt.rowval[subt.colptr[kid]] = parent
            subt.nzval[subt.colptr[kid]] = newwt

        end
    end

    subt = subt + subt'

    ind = findall(marked.<2)
    subt = subt[ind,ind]

    return elims1, elims2, ind, subt
end

function forwardSolve(b, elims1, elims2)

    y = copy(b)

    for i in 1:length(elims1)
        y[elims1[i].parent] += y[elims1[i].nodeid]
    end

    for i in 1:length(elims2)
        wtsum = elims2[i].wt1 + elims2[i].wt2
        y[elims2[i].nbr1] += y[elims2[i].nodeid]*elims2[i].wt1 / wtsum
        y[elims2[i].nbr2] += y[elims2[i].nodeid]*elims2[i].wt2 / wtsum
    end

    return y

end

function backSolve(x, y, elims1, elims2)

    for i in length(elims2):-1:1
        node = elims2[i].nodeid
        wtsum = elims2[i].wt1 + elims2[i].wt2

        x[node] = (elims2[i].wt1*x[elims2[i].nbr1] + elims2[i].wt2*x[elims2[i].nbr2] + y[node])/wtsum
    end

    for i in length(elims1):-1:1
        node = elims1[i].nodeid
        x[node] = x[elims1[i].parent] + y[node]/elims1[i].wtDeg
    end

end

# subtract off the mean from a vector, in place
function subMean!(x::Array{Float64,1})
    n = size(x,1)
    mn = mean(x)
    for i in 1:n,
        x[i] = x[i] - mn
    end
end

"""
    sddmSolver = KMPSDDMSolver(mat; verbose, tol, maxits, maxtime, pcgIts, params::KMPParams)

Solves linear equations in symmetric, diagonally dominant matrices with non-positive off-diagonals.  Based on the paper "Approaching optimality for solving SDD systems" by Koutis, Miller, and Peng, SIAM Journal on Computing, 2014.
"""
function KMPSDDMSolver(mat; verbose=false,
                      tol::Real=1e-6, maxits::Integer=1000, maxtime=Inf, pcgIts=Int[], params::KMPParams=defaultKMPParams)

    n = size(mat,1)
    s = mat*ones(n)

    dmat = diag(mat)
    s = sparse(max.(s,0.0) .* (s .> (dmat*1e-12)))

    if (s == 0)
        error("Matrix was not diagonally dominant.")
    end

    # Force symmetric and diagonal zero
    a = triu(abs.(mat),1)
    a = a + a'

    a1 = [sparse([0 s']); [s a]]

    f1 = KMPLapSolver(a1, verbose=verbose, tol=tol, maxits=maxits, maxtime=maxtime, pcgIts=pcgIts, params=params)

    tol_=tol
    maxits_=maxits
    maxtime_=maxtime
    verbose_=verbose
    pcgIts_=pcgIts

    f = function(b; tol=tol_, maxits=maxits_, maxtime=maxtime_, verbose=verbose_, pcgIts=pcgIts_)

        b1 = [-sum(b);b]
        x1 = f1(b1; verbose=verbose, tol=tol, maxits=maxits, maxtime=maxtime, pcgIts=pcgIts)
        x = x1[2:end] .- x1[1]

        return x

    end

    return f
end

"""
    lapSolver = KMPLapSolver(A; verbose, tol, maxits, maxtime, pcgIts, params::KMPParams)

Solves linear equations in the Laplacian of graph with adjacency matrix `A`.

Based on the paper "Approaching optimality for solving SDD systems" by Koutis, Miller, and Peng, SIAM Journal on Computing, 2014.
"""
function KMPLapSolver(a; verbose=false,
                      tol::Real=1e-6, maxits::Integer=1000, maxtime=Inf, pcgIts=Int[], params::KMPParams=defaultKMPParams)

        return lapWrapComponents(KMPLapSolver1, a, verbose=verbose, tol=tol, maxits=maxits, maxtime=maxtime, pcgIts=pcgIts, params=params)

end

# KMPLapSolver drops right in to this after doing some checks and splitting on components
function KMPLapSolver1(a; verbose=false,
                      tol::Real=1e-6, maxits::Integer=1000, maxtime=Inf, pcgIts=pcgIts, params::KMPParams=defaultKMPParams)

    if (a.n <= params.n0)
        if verbose
            println("The graph is small.  Solve directly")
        end

        return chol_lap(a)
    end

    if (nnz(a) == 2*(a.n - 1))
        if verbose
            println("The graph is a tree.  Solve directly")
        end

        return chol_lap(a)
    end

    if params.treeAlg == :rand
        tree = randishPrim(a)
    else
        tree = akpw(a)
        if verbose
            println("akpw stretch : ", sum(comp_stretches(tree,a))/nnz(a))
        end
    end

    n = size(a,1);

    # if for some reason the graph is a tree, this will fail
    # because the stretches will sum to zero
    # so, default to a direct method

    ord::Array{Int64,1} = Laplacians.dfsOrder(tree)

    # these lines could be MUCH faster
    aord = symPermuteCSC(a,ord)
    tord = symPermuteCSC(tree,ord)

    la = lap(aord)

    #=
    if KMP_SAVEMATS
        KMP_MATS = []
        push!(KMP_MATS,la)
    end
    =#

    fsub = KMPLapPrecon(aord, tord, params, verbose=verbose)

    tol_=tol
    maxits_=maxits
    maxtime_=maxtime
    verbose_=verbose
    pcgIts_=pcgIts

    f = function(b; tol=tol_, maxits=maxits_, maxtime=maxtime_, verbose=verbose_, pcgIts=pcgIts_)

        bord = b[ord] .- mean(b)

        xord = pcg(la, bord, fsub, tol=tol, maxits=maxits, maxtime=maxtime, verbose=verbose, pcgIts=pcgIts)

        x = zeros(Float64,n)
        x[ord] = xord
        subMean!(x) # x = x - mean(x)
        return x
    end

    #=
    if KMP_SAVEMATS
        KMP_FS = []
        push!(KMP_FS,f)
    end
    =#

    return f

end

function KMPLapPrecon(a, t, params; verbose=false)
    n = size(a,1);

    rest = a-t;

    st = comp_stretches(t,rest);
    aveStretch = sum(st)/nnz(rest)

    if params.treeScale > 0

        targetStretch = 1/(params.treeScale*log(n)/log(2))

        fac = aveStretch/targetStretch
        tree = fac*t;

    else

        targetStretch = 1.0

        fac = aveStretch/targetStretch
        tree = t;

    end

    if verbose
        println("aveStretch : ", aveStretch, " fac : ", fac)
    end

    (ai,aj,av) = findnz(triu(rest))
    (si,sj,sv) = findnz(triu(st))
    sv = sv ./ fac

    ijvs = IJVS(ai,aj,av,sv)

    f = KMPLapPreconSub(tree, ijvs, targetStretch, 0, params, verbose=verbose)

    return f
end

function KMPLapPreconSub(tree, ijvs::IJVS, targetStretch::Real, level::Int, params::KMPParams; verbose=false)

    # problem: are forming la before sampling.  should be other way around, at least for top level!
    # that is, we are constructing Heavy, and I don't want to!

    m = size(ijvs.i,1)
    n = size(tree,1)

    if verbose
        println("level ", level, ". Dimension ", n, " off-tree edges : ", m)
    end

    # if is nothing in ijvs
    if m == 0
        la = lap(tree)
        return chol_lap(a)
    end

    ijvs1 = stretchSample(ijvs,targetStretch,params.frac)

    if (length(ijvs1.i) <= params.n0)

        # solve directly

        rest = sparse(ijvs1.i,ijvs1.j,ijvs1.v,n,n)

        return chol_lap(rest + rest' + tree)

    else

        marked = ones(Int64,n)
        marked[ijvs1.i] .= 0
        marked[ijvs1.j] .= 0

        elims1, elims2, ind::Array{Int64,1}, subtree = elimDeg12(tree, marked)

        map = zeros(Int64,n)

        n1 = length(ind)

        map[ind] = collect(1:n1)
        ijvs1.i = map[ijvs1.i]
        ijvs1.j = map[ijvs1.j]

        rest = sparse(ijvs1.i,ijvs1.j,ijvs1.v,n1,n1)
        la1 = lap(rest + rest' + subtree)

        fsub = KMPLapPreconSub(subtree, ijvs1, targetStretch, level+1, params, verbose=verbose)

        f = function(b::Array{Float64,1})
            subMean!(b) # b = b - mean(b)

            y = forwardSolve(b, elims1, elims2)
            ys = y[ind]

            xs = pcg(la1, ys, fsub, tol=0, maxits=params.iters)

            x = zeros(Float64,n)
            x[ind] = xs

            backSolve(x, y, elims1, elims2)
            subMean!(x) # x = x - mean(x)

            return x
        end

    end

    #=
    if KMP_SAVEMATS
        push!(KMP_FS,f)
    end
    =#

    return f

end

#=
goal is to downsample edges by frac.
However, those whose stretches are larger than stretchTarget
might just get their weight reduced instead

This will mostly be the "uniform sampling" envisioned by KMP.
=#
function stretchSample(ijvs::IJVS,stretchTarget::Float64,frac::Float64)

    sampi = Int64[]
    sampj = Int64[]
    sampv = Float64[]
    samps = Float64[]

    m = size(ijvs.i,1)

    stot = sum(ijvs.s)

    # fac = m * frac / stot
    # p = ijvs.s[ind] * fac

    for ind in 1:m
        if ijvs.s[ind] <= stretchTarget

            if rand() < frac
                push!(sampi,ijvs.i[ind])
                push!(sampj,ijvs.j[ind])
                push!(sampv,ijvs.v[ind])
                push!(samps,ijvs.s[ind])
            end

        elseif ijvs.s[ind] <= stretchTarget/frac

            if rand() < frac *  ijvs.s[ind]/ stretchTarget

                push!(sampi,ijvs.i[ind])
                push!(sampj,ijvs.j[ind])
                push!(sampv,ijvs.v[ind] * stretchTarget / ijvs.s[ind])
                push!(samps,stretchTarget)
            end

        else

            push!(sampi,ijvs.i[ind])
            push!(sampj,ijvs.j[ind])
            push!(sampv,ijvs.v[ind]*frac)
            push!(samps,ijvs.s[ind]*frac)

        end
    end
    return IJVS(sampi,sampj,sampv, samps)

end


text:
# # Chapter 2. Interval Arithmetic
# ## 2.4. Floating Point Interval Arithmetic
# ### 2.4.1. A ~~MATLAB~~ Julia Implementation of Interval Arithmetic

#-
using FastRounding
using Printf

const SysFloat = Union{Float32, Float64} # same as FastRounding

#-
struct Interval{T<:SysFloat}
    lo::T
    hi::T
    function Interval(lo::T, hi::T) where T<:SysFloat
        lo <= hi ? new{T}(lo, hi) : throw(ArgumentError("The endpoints do not define an interval."))
    end
end

Interval(lo::Real, hi::Real) = Interval(Float64(lo, RoundDown), Float64(hi, RoundUp))
Interval(x::Real) = Interval(x, x)

function Base.show(io::IO, a::Interval{T}) where T<:SysFloat
    print(io, "Interval{$(T)}(", @sprintf("%17.17f", a.lo), ", ", @sprintf("%17.17f", a.hi), ")")
end

#-
a = Interval(3, 4)

#-
b = Interval(2, 5)

#-
c = Interval(1)

#-
Interval(3.0f0, 4.0f0)

#-
Interval(1//10)

#-
Interval(π)

#-
## NG
## Interval(2.0, 1.0)

# Define sum of two intervals.

#-
function Base.:+(a::Interval, b::Interval)
    lo = add_round(a.lo, b.lo, RoundDown)
    hi = add_round(a.hi, b.hi, RoundUp)
    Interval(lo, hi)
end

Base.:+(a::Real, b::Interval) = Interval(a) + b
Base.:+(a::Interval, b::Real) = a + Interval(b)

Base.:+(a::Interval) = a

#-
Interval(1, 2) + Interval(3, 4)

#-
Interval(1, 2) + Interval(3, 3)

#-
Interval(1, 2) + 3

#-
3 + Interval(1, 2)

# Define other unary and binary functions

#-
function Base.:-(a::Interval, b::Interval)
    lo = sub_round(a.lo, b.hi, RoundDown)
    hi = sub_round(a.hi, b.lo, RoundUp)
    Interval(lo, hi)
end

Base.:-(a::Real, b::Interval) = Interval(a) - b
Base.:-(a::Interval, b::Real) = a - Interval(b)

Base.:-(a::Interval) = 0 - a

function Base.:*(a::Interval, b::Interval)
    lo = min(mul_round(a.lo, b.lo, RoundDown),
             mul_round(a.lo, b.hi, RoundDown),
             mul_round(a.hi, b.lo, RoundDown),
             mul_round(a.hi, b.hi, RoundDown))
    hi = max(mul_round(a.lo, b.lo, RoundUp),
             mul_round(a.lo, b.hi, RoundUp),
             mul_round(a.hi, b.lo, RoundUp),
             mul_round(a.hi, b.hi, RoundUp))
    Interval(lo, hi)
end

Base.:*(a::Real, b::Interval) = Interval(a) * b
Base.:*(a::Interval, b::Real) = a * Interval(b)

function Base.:/(a::Interval, b::Interval)
    if (b.lo <= 0.0) & (0.0 <= b.hi)
        throw(ArgumentError("Denominator straddles zero."))
    end

    lo = min(div_round(a.lo, b.lo, RoundDown),
             div_round(a.lo, b.hi, RoundDown),
             div_round(a.hi, b.lo, RoundDown),
             div_round(a.hi, b.hi, RoundDown))
    hi = max(div_round(a.lo, b.lo, RoundUp),
             div_round(a.lo, b.hi, RoundUp),
             div_round(a.hi, b.lo, RoundUp),
             div_round(a.hi, b.hi, RoundUp))
    Interval(lo, hi)
end

Base.:/(a::Real, b::Interval) = Interval(a) / b
Base.:/(a::Interval, b::Real) = a / Interval(b)

# Unary

#-
+a

#-
-a

# Binary

#-
a + b

#-
a - b

#-
a * b

#-
a / b

#-
Interval(1/10)

#-
Interval(1)/10

# Sub-distributivity

#-
Interval(1//10)

#-
a = Interval(-1, 1)
b = Interval(-1, 0)
c = Interval(3, 4)

#-
a * (b + c)

#-
a * b + a * c

# Defial equal

#-
import Base.==

==(a::Interval, b::Interval) = a.lo == b.lo && a.hi == b.hi

#-
a = Interval(1, 10)
b = Interval(-2, 3)
c = Interval(3, 5)

a == a

#-
a == b

#-
a == c

#-
a != a

#-
a != b

#-
a != c

# Inclusion

#-
Base.issubset(a::Interval, b::Interval) = b.lo <= a.lo && a.hi <= b.hi
Base.:⊊(a::Interval, b::Interval) = a ⊆ b && a != b

#-
a ⊆ a

#-
a ⊆ b

#-
b ⊆ a

#-
a ⊊ c

#-
c ⊊ a

#-
a ⊊ a

#-
a ⊇ a

#-
a ⊇ c

#-
a ⊋ a

#-
a ⊋ c

# Intersection

#-
function Base.:∩(a::Interval, b::Interval)
    if a.hi < b.lo || b.hi < a.lo
        nothing
    else
        Interval(max(a.lo, b.lo), min(a.hi, b.hi))
    end
end

#-
a = Interval(1, 3)
b = Interval(4, 5)
c = Interval(2, 5)

a ∩ b

#-
isnothing(a ∩ c)

#-
isnothing(a ∩ b)

# Union

#-
Base.:∪(a::Interval, b::Interval) = Interval(min(a.lo, b.lo), max(a.hi, b.hi))

#-
a ∪ c

#-
a ∪ b

max_stars_repo_name: matsueushi/study-validated-numerics

max_issues_repo_name: matsueushi/study-validated-numerics

max_forks_repo_name: matsueushi/study-validated-numerics


text:
module Package1

require("modules_ext.jl")
using LibA
using LibB.varB
using LibC: varC, funcC
import LibD.funcD
importall LibE

export Type1, perc

include("file1.jl")
include("file2.jl")

# code
type Type1
    total
end

perc(a::Type) = a.total * 0.01

end

max_stars_repo_path: Module 1/Chapter06/Old/modules.jl

max_issues_repo_path: Module 1/Chapter06/Old/modules.jl

max_forks_repo_path: Module 1/Chapter06/Old/modules.jl


text:
# Configuration Parameters
"""
    sdad_setup!(;db_user::AbstractString = "",
                 db_pwd::AbstractString = "",
                 token::AbstractString = "")
This function writes the configurations to `confs/config.simple`
# Arguments
- `db_user`: Your username for the `postgis_2` [database](http://sdad.policy-analytics.net:8080/?pgsql=postgis_2&db=sdad_data&ns=onet) (i.e., your UVA computing ID).
- `db_pwd`: The password for the `postgis_2` [database](http://sdad.policy-analytics.net:8080/?pgsql=postgis_2&db=sdad_data&ns=onet).
- `onet_auth`: Basic authentication [token](https://services.onetcenter.org/developer/) for the O*Net REST API.
# Examples
## Good!
```
julia> sdad_setup!(db_user = "jbs3hp",
                   db_pwd = "MyVerySafePwd",
                   onet_auth = "ASZ1QqATfYXdOLFnr6EVfxyHymx4kTGMLQgXCs==")
😃
```
## Incomplete!
```
julia> sdad_setup!(db_user = "jbs3hp",
                   db_pwd = "MyVerySafePwd")
Warning: github_login has not been defined.
😞
```
# Updating the configuration file.
```
julia> using ONet
julia> sdad_setup!(db_user = "jbs3hp",
                   db_pwd = "MyVerySafePwd",
                   onet_auth = "ASZ1QqATfYXdOLFnr6EVfxyHymx4kTGMLQgXCs==")
julia> exit()
> julia
julia> using ONet # config up-to-date
```
!!! note
    Updating the configuration file requires a restart of the session.
"""
function sdad_setup!(;db_user::AbstractString = "",
                      db_pwd::AbstractString = "",
                      onet_auth::AbstractString = "")
    args = ["db_user", "db_pwd", "onet_auth"]
    isdir(joinpath(dirname(@__FILE__), "..", "confs")) || mkdir(joinpath(dirname(@__FILE__), "..", "confs"))
    isfile(joinpath(dirname(@__FILE__), "..", "confs", "config.simple")) ||
        touch(joinpath(dirname(@__FILE__), "..", "confs", "config.simple"))
    conf = ConfParse(joinpath(dirname(@__FILE__), "..", "confs", "config.simple"),
                     "simple")
    parse_conf!(conf)
    for (key, val) ∈ zip(args, [db_user, db_pwd, onet_auth])
        isempty(val) || commit!(conf, key, val)
    end
    save!(conf)
    notdefined = filter(key -> !haskey(conf, key),
                        args)
    foreach(key -> @warn("$key has not been defined."),
            notdefined)
    if isinteractive()
        if isempty(notdefined)
            println("😃")
        else
            println("😞")
        end
    end
end
isfile(joinpath(dirname(@__FILE__), "..", "confs", "config.simple")) ||
    sdad_setup!(db_user = get(ENV, "db_user", ""),
                db_pwd = get(ENV, "db_pwd", ""),
                onet_auth = get(ENV, "onet_auth", ""))
const conf = ConfParse(joinpath(dirname(@__FILE__), "..", "confs", "config.simple"),
                       "simple");
parse_conf!(conf);
"""
    db_user::String
The username for the [database](http://sdad.policy-analytics.net:8080/?pgsql=postgis_2&db=sdad_data&ns=onet) (i.e., your UVA computing ID).
"""
const db_user = haskey(conf, "db_user") ? retrieve(conf, "db_user") : "";
"""
    db_pwd::String
The password for the [database](http://sdad.policy-analytics.net:8080/?pgsql=postgis_2&db=sdad_data&ns=onet).
"""
const db_pwd = haskey(conf, "db_pwd") ? retrieve(conf, "db_pwd") : "";
"""
    onet_auth::String
The application O*Net basic authentication token.
"""
const onet_auth = haskey(conf, "onet_auth") ? retrieve(conf, "onet_auth") : "";
"""
    headers = ["Authorization" => "Basic \$onet_auth",
               "Accept" => "application/json"]
Headers to be passed to the ONet REST API.
"""
const headers = ["Authorization" => "Basic $onet_auth",
                 "Accept" => "application/json"]
"""
    db_host = "postgis_2"
Host for the [database](http://sdad.policy-analytics.net:8080/?pgsql=postgis_2&db=sdad_data&schema=onet).
"""
const db_host = "postgis_2";
"""
    db_port = 5432
Port for the `postgis_2` in the [database](http://sdad.policy-analytics.net:8080/?pgsql=postgis_2&db=sdad_data&schema=onet).
"""
const db_port = 5432
"""
    dbname = "sdad_data"
Database for the SDAD data.
"""
const dbname = "sdad_data";


text:
function loaddep(dir)
    fn = first(readdir(dir))
    return CSV.read(joinpath(dir, fn), missingstring="NA")
end

pokes = loaddep(PokesDatasetFolder)
streaks = loaddep(StreaksDatasetFolder)
#
# println(names(pokes))
# t = pokes[1,:Session]
# pokes.Reward
# sort(pokes, :Poke)
#
# function minibout(sides)
#     # check if animal changed side and count consecutive omissions
#     s = findfirst(t -> t != sides[1], sides)
#     return if isnothing(s)
#         (omissionsplusone = length(sides), leaves = false)
#     else
#         # we discard the rest of the data which are the beginning of the next streak (if doesn't start with a reward)
#         (omissionsplusone = s-1, leaves = true)
#     end
# end
#
# function only(v)
#     s = first(v)
#     @assert all(==(s), v)
#     return s
# end
#
# function bouts(session, carry = [:MouseID]; IPIthreshold = 1)
#     sorted = sort(session, :Poke)
#     sorted.RewardCount = cumsum(sorted.Reward)
#     sorted.IPI = sorted.PokeIn .- lag(sorted.PokeOut)
#     fs = filter(row -> row.RewardCount >= 1, sorted) #discard omissions before first reward
#     bts = by(fs, :RewardCount, :Side => minibout)
#     extra_vars = Dict{Symbol, Any}(c => c => only for c in carry)
#     extra_vars[:MaxIPI] = (:Pre_Interpoke, :Side) => function (nt)
#         ipis, sides = nt.Pre_Interpoke, nt.Side
#         acc = missing
#         t = first(sides)
#         for (ipi, side) in zip(ipis, sides)
#             side != t && break
#             if ismissing(acc) && !isi || !ismis
#         findfirst(t -> !ismissing(t) && t > IPIthreshold, ipis)
#     extra_info = by(fs, :RewardCount; extra_vars...) # add info from columns in carry
#     data = join(bts, extra_info, on = :RewardCount)
#     filter(row -> isnothing(row.LongIPI) || row.LongIPI)
# end
#
# by(pokes, :Session) do dd
#     bouts(dd, [:MouseID, :Protocol, :Wall, :Streak, :Date])
# end
#
# bouts(filter(row -> row.Session == t, pokes), [:MouseID])
# extra_info = Dict(c => c => first for c in [:MouseID])
# by(pokes, :Session; extra_info...)
#
# pokes.Pre_Interpoke


text:
function calc_cost(reactor, R_0, B_0)
  this_term = 1 + reactor.b / ( reactor.epsilon * R_0 * reactor.kappa )
  that_term = 1 + reactor.b / ( reactor.epsilon * R_0 )

  K_TF = 10 * this_term * that_term * reactor.epsilon^2 * reactor.kappa
  K_C = K_TF / reactor.P_F

  K_C * (R_0^3) * (B_0^2)
end


text:
using LeafAreaIndex
LAI = LeafAreaIndex
using Test
import FixedPointNumbers

## CALIBRATION ##

Rmax1 = 3
N = 2Rmax1 + 1
testimg1 = reshape(collect(1:N^2), N, N)

test1cal = CameraLens((N, N), (Rmax1+1, Rmax1+1), Rmax1, [1/(pi/2)])

@test test1cal.fθρ(0) ≈ 0
@test test1cal.fρθ(0) ≈ 0
@test test1cal.fθρ(pi/2) ≈ Rmax1
@test test1cal.fρθ(Rmax1) ≈ π/2
@test test1cal.params.imgsize == size(testimg1)

@test test1cal.sort_ind == testimg1[test1cal.sort_ind]
# @test testimg1[test1cal.spiral_ind] == [13,9,14,19,18,17,12,7,8,15,23,11,3]

# @test test1cal.ρ²sort == UInt32[0,1,1,1,1,2,2,2,2,4,4,4,4]
#@test test1cal.ϕsort ≈ collect([0, -π/2:π/2:π, -3π/4:π/2:3π/4, -π/2:π/2:π])

## THRESHOLDING ##

testimg2 = 0.5 * ones(FixedPointNumbers.Normed{UInt16, 16}, 5, 5)
testimg2[2:4, 2:4] .= 0.7
testimg2[3, 3] = 1
@test threshold(testimg2, LAI.MinimumThreshold()) > 0.5

# circqueue
let
	a = [1; 2; 3]
	c = LeafAreaIndex.circqueue(a)
	@test LeafAreaIndex.pushshift!(c, 5) == 1
	@test LeafAreaIndex.pushshift!(c, 6) == 2
	@test LeafAreaIndex.pushshift!(c, 7) == 3
	@test LeafAreaIndex.pushshift!(c, 8) == 5
end


text:
#!/Users/yoyoyoju/.julia/v0.5/Contents/Resources/julia/bin/julia

workspace()
using Gallium
using HoTRG_lemon.SimulatorModule
using HoTRG_lemon.SpinModule
using HoTRG_lemon.LatticeModule
using TensorOperations
using TensorMatrices_lemon

#= debug
Gallium.breakpoint(renormalizeSpace!)
Gallium.breakpoint(renormalizeTrotter!)
=#

dimM = 5
trotterparameter = 0.01
trotteriteration = 10

isinginfo = SpinInfo("quantum_ising_2",1.0e-13,0.0)
trotterinfo = TrotterInfo(trotterparameter, trotteriteration)
ising = QuantumIsingModel(isinginfo, trotterinfo)
lattice = buildLattice("quantum_2d_fractal", ising)
simulator = buildSimulator(lattice,dimM)

# sqLattice = buildLattice("quantum_2d_square", ising)
# sqSimulator = buildSimulator(sqLattice,dimM)

# println(simulator())
# println(sqSimulator())

# fieldrange = linspace(0.1, 3.5, 10)
fieldrange = 0.1
simulatorQuantum(fieldrange, simulator; filename = "log_mag.txt", printlog="mag")
 # simulatorQuantum(fieldrange, sqSimulator; filename = "q2s_m6_t30_ori_m7.txt")
#tensorT = getTensorT(simulator)[1]
#tenmatMMd = getTenmatMMd(tensorT)
#Ul, lambdaVector, Uld = svd(tenmatMMd.matrix)
#trunUl = truncMatrixU(Ul,dimM)
#tensorU = matU2tenU(trunUl)


text:
import Base: eltype

# definition of the `Status` struct.
# TODO: this bakes in a lot of assumptions about ordering and padding
if !@isdefined(Status)
    let
        struct_fields = Any[]
        offset = 0
        i = 0

        while offset < MPI_Status_Source_offset
            push!(struct_fields, :($(Symbol(:_pad, i+=1))::Cint))
            offset += 4
        end
        @assert offset == MPI_Status_Source_offset
        push!(struct_fields, :(source::Cint))
        offset += 4

        while offset < MPI_Status_Tag_offset
            push!(struct_fields, :($(Symbol(:_pad, i+=1))::Cint))
            offset += 4
        end
        @assert offset == MPI_Status_Tag_offset
        push!(struct_fields, :(tag::Cint))
        offset += 4

        while offset < MPI_Status_Error_offset
            push!(struct_fields, :($(Symbol(:_pad, i+=1))::Cint))
            offset += 4
        end
        @assert offset == MPI_Status_Error_offset
        push!(struct_fields, :(error::Cint))
        offset += 4

        while offset < MPI_Status_size
            push!(struct_fields, :($(Symbol(:_pad, i+=1))::Cint))
            offset += 4
        end

        @eval struct Status
            $(struct_fields...)
        end
    end

    @assert fieldoffset(Status, findfirst(isequal(:source), fieldnames(Status))) == MPI_Status_Source_offset
    @assert fieldoffset(Status, findfirst(isequal(:tag), fieldnames(Status))) == MPI_Status_Tag_offset
    @assert fieldoffset(Status, findfirst(isequal(:error), fieldnames(Status))) == MPI_Status_Error_offset
    @assert sizeof(Status) == MPI_Status_size
end

"""
    MPI.Status

The status of an MPI receive communication. It has 3 accessible fields

- `source`: source of the received message
- `tag`: tag of the received message
- `error`: error code. This is only set if a function returns multiple statuses.

Additionally, the accessor function [`MPI.Get_count`](@ref) can be used to determine the
number of entries received.
"""
Status

Get_source(status::Status) = Int(status.source)
Get_tag(status::Status) = Int(status.tag)
Get_error(status::Status) = Int(status.error)

"""
    MPI.Request

An MPI Request object, representing a non-blocking communication. This also contains a
reference to the buffer used in the communication to ensure it isn't garbage-collected
during communication.
"""
@mpi_handle Request buffer

const REQUEST_NULL = _Request(MPI_REQUEST_NULL, nothing)
Request() = Request(REQUEST_NULL.val, nothing)

"""
    ismessage, (status|nothing) = Iprobe(src::Integer, tag::Integer, comm::Comm)

Blocks until there is a message that can be received matching `src`, `tag` and
`comm`. Returns the corresponding [`Status`](@ref) object.

# External links
$(_doc_external("MPI_Probe"))
"""
function Probe(src::Integer, tag::Integer, comm::Comm)
    stat_ref = Ref{Status}()
    @mpichk ccall((:MPI_Probe, libmpi), Cint,
          (Cint, Cint, MPI_Comm, Ptr{Status}),
          src, tag, comm, stat_ref)
    stat_ref[]
end

"""
    ismessage, (status|nothing) = Iprobe(src::Integer, tag::Integer, comm::Comm)

Checks if there is a message that can be received matching `src`, `tag` and `comm`. If so,
returns a tuple `true` and a [`Status`](@ref) object, otherwise returns a tuple `false, nothing`.

# External links
$(_doc_external("MPI_Iprobe"))
"""
function Iprobe(src::Integer, tag::Integer, comm::Comm)
    flag = Ref{Cint}()
    stat_ref = Ref{Status}()
    @mpichk ccall((:MPI_Iprobe, libmpi), Cint,
          (Cint, Cint, MPI_Comm, Ptr{Cint}, Ptr{Status}),
          src, tag, comm, flag, stat_ref)
    if flag[] == 0
        return false, nothing
    end
    true, stat_ref[]
end

"""
    MPI.Get_count(status::Status, T)

The number of entries received. `T` should match the argument provided by the receive call that set the status variable.

If the number of entries received exceeds the limits of the count parameter, then it returns `MPI_UNDEFINED`.

# External links
$(_doc_external("MPI_Get_count"))
"""
function Get_count(stat::Status, datatype::Union{MPI_Datatype, Datatype})
    count = Ref{Cint}()
    @mpichk ccall((:MPI_Get_count, libmpi), Cint,
                  (Ptr{Status}, MPI_Datatype, Ptr{Cint}),
                  Ref(stat), datatype, count)
    Int(count[])
end
Get_count(stat::Status, ::Type{T}) where {T} = Get_count(stat, mpitype(T))

"""
    Send(buf, [count::Integer, [datatype::Datatype,]]
         dest::Integer, tag::Integer, comm::Comm) where T

Perform a blocking send of `count` elements of type `datatype` from `buf` to MPI
rank `dest` of communicator `comm` using the message tag `tag`

If not provided, `datatype` and `count` are derived from the element type and length of
`buf`, respectively.

# External links
$(_doc_external("MPI_Send"))
"""
function Send(buf, count::Integer, datatype::Union{Datatype, MPI_Datatype},
              dest::Integer, tag::Integer, comm::Comm)
    # int MPI_Send(const void* buf, int count, MPI_Datatype datatype, int dest,
    #              int tag, MPI_Comm comm)
    @mpichk ccall((:MPI_Send, libmpi), Cint,
          (MPIPtr, Cint, MPI_Datatype, Cint, Cint, MPI_Comm),
          buf, count, datatype, dest, tag, comm)
    return nothing
end

Send(buf, count::Integer, dest::Integer, tag::Integer, comm::Comm) =
    Send(buf, count, mpitype(eltype(buf)), dest, tag, comm)
Send(buf::AbstractArray, dest::Integer, tag::Integer, comm::Comm) =
    Send(buf, length(buf), dest, tag, comm)

"""
    Send(obj::T, dest::Integer, tag::Integer, comm::Comm) where T

Complete a blocking send of `obj` to MPI rank `dest` of communicator `comm`
using with the message tag `tag`.
"""
function Send(obj::T, dest::Integer, tag::Integer, comm::Comm) where T
    buf = [obj]
    Send(buf, dest, tag, comm)
end

"""
    send(obj, dest::Integer, tag::Integer, comm::Comm)

Complete a blocking send of using a serialized version of `obj` to MPI rank
`dest` of communicator `comm` using with the message tag `tag`.
"""
function send(obj, dest::Integer, tag::Integer, comm::Comm)
    buf = MPI.serialize(obj)
    Send(buf, dest, tag, comm)
end

"""
    Isend(buf, [count::Integer, [datatype::Datatype,]]
          dest::Integer, tag::Integer, comm::Comm)

Starts a nonblocking send of `count` elements of type `datatype` from `buf` to
MPI rank `dest` of communicator `comm` using with the message tag `tag`

If not provided, `datatype` and `count` are derived from the element type and length of
`buf`, respectively.

Returns the [`Request`](@ref) object for the nonblocking send.

# External links
$(_doc_external("MPI_Isend"))
"""
function Isend(buf, count::Integer, datatype::Union{Datatype, MPI_Datatype},
               dest::Integer, tag::Integer, comm::Comm)
    req = Request()
    # int MPI_Isend(const void* buf, int count, MPI_Datatype datatype, int dest,
    #               int tag, MPI_Comm comm, MPI_Request *request)
    @mpichk ccall((:MPI_Isend, libmpi), Cint,
          (MPIPtr, Cint, MPI_Datatype, Cint, Cint, MPI_Comm, Ptr{MPI_Request}),
                  buf, count, datatype, dest, tag, comm, req)
    req.buffer = buf
    return req
end

Isend(buf, count::Integer, dest::Integer, tag::Integer, comm::Comm) =
    Isend(buf, count, mpitype(eltype(buf)), dest, tag, comm)
Isend(buf::AbstractArray, dest::Integer, tag::Integer, comm::Comm) =
    Isend(buf, length(buf), dest, tag, comm)

"""
    Isend(obj::T, dest::Integer, tag::Integer, comm::Comm) where T

Starts a nonblocking send of `obj` to MPI rank `dest` of communicator `comm`
using with the message tag `tag`.

Returns the commication `Request` for the nonblocking send.
"""
function Isend(obj::T, dest::Integer, tag::Integer, comm::Comm) where T
    buf = [obj]
    Isend(buf, dest, tag, comm)
end

"""
    isend(obj, dest::Integer, tag::Integer, comm::Comm)

Starts a nonblocking send of using a serialized version of `obj` to MPI rank
`dest` of communicator `comm` using with the message tag `tag`.

Returns the commication `Request` for the nonblocking send.
"""
function isend(obj, dest::Integer, tag::Integer, comm::Comm)
    buf = MPI.serialize(obj)
    Isend(buf, dest, tag, comm)
end

"""
    Recv!(buf, [count::Integer, [datatype::Datatype,]]
          src::Integer, tag::Integer, comm::Comm)

Completes a blocking receive of up to `count` elements of type `datatype` into `buf`
from MPI rank `src` of communicator `comm` using with the message tag `tag`

If not provided, `datatype` and `count` are derived from the element type and length of
`buf`, respectively.

Returns the [`Status`](@ref) of the receive.

# External links
$(_doc_external("MPI_Recv"))
"""
function Recv!(buf, count::Integer, datatype::Union{Datatype,MPI_Datatype}, src::Integer,
               tag::Integer, comm::Comm)
    stat_ref = Ref{Status}()
    # int MPI_Recv(void* buf, int count, MPI_Datatype datatype, int source,
    #              int tag, MPI_Comm comm, MPI_Status *status)
    @mpichk ccall((:MPI_Recv, libmpi), Cint,
                  (MPIPtr, Cint, MPI_Datatype, Cint, Cint, MPI_Comm, Ptr{Status}),
                  buf, count, datatype, src, tag, comm, stat_ref)
    return stat_ref[]
end

Recv!(buf, count::Integer, src::Integer, tag::Integer, comm::Comm) =
    Recv!(buf, count, mpitype(eltype(buf)), src, tag, comm)
Recv!(buf::AbstractArray, src::Integer, tag::Integer, comm::Comm) =
    Recv!(buf, length(buf), src, tag, comm)

function Recv(::Type{T}, src::Integer, tag::Integer, comm::Comm) where T
    buf = Ref{T}()
    stat = Recv!(buf, 1, src, tag, comm)
    (buf[], stat)
end

function recv(src::Integer, tag::Integer, comm::Comm)
    stat = Probe(src, tag, comm)
    count = Get_count(stat, UInt8)
    buf = Array{UInt8}(undef, count)
    stat = Recv!(buf, Get_source(stat), Get_tag(stat), comm)
    (MPI.deserialize(buf), stat)
end

"""
    Irecv!(buf, [count::Integer, [datatype::Datatype,]]
           src::Integer, tag::Integer, comm::Comm) where T

Starts a nonblocking receive of up to `count` elements of type `datatype` into `buf`
from MPI rank `src` of communicator `comm` using with the message tag `tag`

Returns the [`Request`](@ref) for the nonblocking receive.

# External links
$(_doc_external("MPI_Irecv"))
"""
function Irecv!(buf, count::Integer, datatype::Union{Datatype, MPI_Datatype},
                    src::Integer, tag::Integer, comm::Comm) where T
    req = Request()
    # int MPI_Irecv(void* buf, int count, MPI_Datatype datatype, int source,
    #               int tag, MPI_Comm comm, MPI_Request *request)
    @mpichk ccall((:MPI_Irecv, libmpi), Cint,
                  (MPIPtr, Cint, MPI_Datatype, Cint, Cint, MPI_Comm, Ptr{MPI_Request}),
                  buf, count, datatype, src, tag, comm, req)
    req.buffer = buf
    return req
end

Irecv!(buf, count::Integer, src::Integer, tag::Integer, comm::Comm) =
    Irecv!(buf, count, mpitype(eltype(buf)), src, tag, comm)
Irecv!(buf::AbstractArray, src::Integer, tag::Integer, comm::Comm) =
    Irecv!(buf, length(buf), src, tag, comm)

function irecv(src::Integer, tag::Integer, comm::Comm)
    (flag, stat) = Iprobe(src, tag, comm)
    if !flag
        return (false, nothing, nothing)
    end
    count = Get_count(stat, UInt8)
    buf = Array{UInt8}(undef, count)
    stat = Recv!(buf, Get_source(stat), Get_tag(stat), comm)
    (true, MPI.deserialize(buf), stat)
end

"""
    Sendrecv(sendbuf, [sendcount::Integer, [sendtype::Union{Datatype, MPI_Datatype}]],
             dest::Integer, sendtag::Integer,
             recvbuf, [recvcount::Integer, [recvtype::Union{Datatype, MPI_Datatype}]],
             source::Integer, recvtag::Integer,
             comm::Comm)

Complete a blocking send-receive operation over the MPI communicator `comm`. Send
`sendcount` elements of type `sendtype` from `sendbuf` to the MPI rank `dest` using message
tag `tag`, and receive `recvcount` elements of type `recvtype` from MPI rank `source` into
the buffer `recvbuf` using message tag `tag`. Return a [`Status`](@ref) object.

If not provided, `sendtype`/`recvtype` and `sendcount`/`recvcount` are derived from the
element type and length of `sendbuf`/`recvbuf`, respectively.

# External links
$(_doc_external("MPI_Sendrecv"))
"""
function Sendrecv(sendbuf, sendcount::Integer, sendtype::Union{Datatype, MPI_Datatype},   dest::Integer, sendtag::Integer,
                  recvbuf, recvcount::Integer, recvtype::Union{Datatype, MPI_Datatype}, source::Integer, recvtag::Integer,
                  comm::Comm)
    # int MPI_Sendrecv(const void *sendbuf, int sendcount, MPI_Datatype sendtype, int dest,   int sendtag,
    #                        void *recvbuf, int recvcount, MPI_Datatype recvtype, int source, int recvtag,
    #                    MPI_Comm comm, MPI_Status *status)
    stat_ref = Ref{Status}()
    @mpichk ccall((:MPI_Sendrecv, libmpi), Cint,
                  (MPIPtr, Cint, MPI_Datatype, Cint, Cint,
                   MPIPtr, Cint, MPI_Datatype, Cint, Cint,
                 MPI_Comm, Ptr{Status}),
                   sendbuf, sendcount, sendtype, dest,   sendtag,
                   recvbuf, recvcount, recvtype, source, recvtag, comm, stat_ref)
    return stat_ref[]
end

function Sendrecv(sendbuf, sendcount::Integer,   dest::Integer, sendtag::Integer,
                  recvbuf, recvcount::Integer, source::Integer, recvtag::Integer,
                  comm::Comm)
    return Sendrecv(sendbuf, sendcount, mpitype(eltype(sendbuf)), dest,   sendtag,
                    recvbuf, recvcount, mpitype(eltype(recvbuf)), source, recvtag, comm)
end
function Sendrecv(sendbuf::AbstractArray,   dest::Integer, sendtag::Integer,
                  recvbuf::AbstractArray, source::Integer, recvtag::Integer,
                  comm::Comm)
    return Sendrecv(sendbuf, length(sendbuf), dest,   sendtag,
                    recvbuf, length(recvbuf), source, recvtag, comm)
end

"""
    status = Wait!(req::Request)

Block until the request `req` is complete and deallocated. Returns the [`Status`](@ref) of the request.

# External links
$(_doc_external("MPI_Wait"))
"""
function Wait!(req::Request)
    stat_ref = Ref{Status}()
    # int MPI_Wait(MPI_Request *request, MPI_Status *status)
    @mpichk ccall((:MPI_Wait, libmpi), Cint,
                  (Ptr{MPI_Request}, Ptr{Status}),
                  req, stat_ref)
    req.buffer = nothing
    stat
end

"""
    iscomplete, (status|nothing) = Test!(req::Request)

Check if the request `req` is complete. If so, the request is deallocated and a tuple of
`true` and the [`Status`](@ref) of the request is returned. Otherwise a tuple of
`false, nothing` is returned.

# External links
$(_doc_external("MPI_Test"))
"""
function Test!(req::Request)
    flag = Ref{Cint}()
    stat_ref = Ref{Status}()
    # int MPI_Test(MPI_Request *request, int *flag, MPI_Status *status)
    @mpichk ccall((:MPI_Test, libmpi), Cint,
                  (Ptr{MPI_Request}, Ptr{Cint}, Ptr{Status}),
                  req, flag, stat_ref)
    if flag[] == 0
        return (false, nothing)
    end
    req.buffer = nothing
    (true, stat_ref[])
end

"""
    statuses = Waitall!(reqs::Vector{Request})

Block until all the requests in the array `reqs` are complete. Returns an array
of the [`Status`](@ref) objects corresponding to each request.

# External links
$(_doc_external("MPI_Waitall"))
"""
function Waitall!(reqs::Vector{Request})
    count = length(reqs)
    reqvals = [reqs[i].val for i in 1:count]
    stats = Array{Status}(undef, count)
    # int MPI_Waitall(int count, MPI_Request array_of_requests[],
    #                 MPI_Status array_of_statuses[])
    @mpichk ccall((:MPI_Waitall, libmpi), Cint,
                  (Cint, Ptr{MPI_Request}, Ptr{Status}),
                  count, reqvals, stats)
    for i in 1:count
        reqs[i].val = reqvals[i]
        reqs[i].buffer = nothing
    end
    return stats
end

"""
    iscomplete, (statuses|nothing) = Testall!(reqs::Vector{Request})

Check if all the requests in the array `reqs` are complete. If so, the requests are
deallocated and a tuple of `true` and an array of the [`Status`](@ref) objects
corresponding to each request is returned. Otherwise no requests are modified a tuple of
`false, nothing` is returned.

# External links
$(_doc_external("MPI_Testall"))
"""
function Testall!(reqs::Vector{Request})
    count = length(reqs)
    reqvals = [reqs[i].val for i in 1:count]
    flag = Ref{Cint}()
    stats = Array{Status}(undef, count)
    # int MPI_Testall(int count, MPI_Request array_of_requests[], int *flag,
    #                 MPI_Status array_of_statuses[])
    @mpichk ccall((:MPI_Testall, libmpi), Cint,
                  (Cint, Ptr{MPI_Request}, Ptr{Cint}, Ptr{Status}),
                  count, reqvals, flag, stats)
    if flag[] == 0
        return (false, nothing)
    end
    for i in 1:count
        reqs[i].val = reqvals[i]
        reqs[i].buffer = nothing
    end
    (true, stats)
end

"""
    (index, status) = Waitany!(reqs::Vector{Request})

Blocks until one of the requests in the array `reqs` is complete: if more than one is
complete, one is chosen arbitrarily. The request is deallocated and a tuple of the index
of the completed request and its [`Status`](@ref) is returned.

# External links
$(_doc_external("MPI_Waitany"))
"""
function Waitany!(reqs::Vector{Request})
    count = length(reqs)
    reqvals = [reqs[i].val for i in 1:count]
    ind = Ref{Cint}()
    stat_ref = Ref{Status}()
    # int MPI_Waitany(int count, MPI_Request array_of_requests[], int *index,
    #                 MPI_Status *status)
    @mpichk ccall((:MPI_Waitany, libmpi), Cint,
                  (Cint, Ptr{MPI_Request}, Ptr{Cint}, Ptr{Status}),
                  count, reqvals, ind, stat_ref)
    index = Int(ind[]) + 1
    reqs[index].val = reqvals[index]
    reqs[index].buffer = nothing
    (index, stat_ref[])
end

"""
    iscomplete, (index, status | 0, nothing) = Testany!(reqs::Vector{Request})

Check if any one of the requests in the array `reqs` is complete: if more than one is
complete, one is chosen arbitrarily. If so, the request is deallocated an a tuple of
`true`, its index and its [`Status`](@ref) is returned. Otherwise a tuple of `false, 0,
nothing` is returned.

# External links
$(_doc_external("MPI_Testany"))
"""
function Testany!(reqs::Vector{Request})
    count = length(reqs)
    reqvals = [req.val for req in reqs]
    ind = Ref{Cint}()
    flag = Ref{Cint}()
    stat_ref = Ref{Status}()
    # int MPI_Testany(int count, MPI_Request array_of_requests[], int *index,
    #                 int *flag, MPI_Status *status)
    @mpichk ccall((:MPI_Testany, libmpi), Cint,
                  (Cint, Ptr{MPI_Request}, Ptr{Cint}, Ptr{Cint}, Ptr{Status}),
                  count, reqvals, ind, flag, stat_ref)
    if flag[] == 0 || ind[] == MPI_UNDEFINED
        return (false, 0, nothing)
    end
    index = Int(ind[]) + 1
    reqs[index].val = reqvals[index]
    reqs[index].buffer = nothing
    (true, index, stat_ref[])
end

"""
    (indices, statuses) = Waitsome!(reqs::Vector{Request})

Block until at least one of the requests in the array `reqs` is complete. The completed
requests are deallocated, and a tuple of their indices in `reqs` and their corresponding
[`Status`](@ref) objects are returned.

# External links
$(_doc_external("MPI_Waitsome"))
"""
function Waitsome!(reqs::Vector{Request})
    count = length(reqs)
    reqvals = [reqs[i].val for i in 1:count]
    outcnt = Ref{Cint}()
    inds = Array{Cint}(undef, count)
    stats = Array{Status}(undef,count)
    # int MPI_Waitsome(int incount, MPI_Request array_of_requests[],
    #                  int *outcount, int array_of_indices[],
    #                  MPI_Status array_of_statuses[])
    @mpichk ccall((:MPI_Waitsome, libmpi), Cint,
          (Cint, Ptr{MPI_Request}, Ptr{Cint}, Ptr{Cint}, Ptr{Status}),
          count, reqvals, outcnt, inds, stats)
    outcount = Int(outcnt[])
    # This can happen if there were no valid requests
    if outcount == MPI_UNDEFINED
        outcount = 0
    end
    indices = Array{Int}(undef, outcount)
    for i in 1:outcount
        ind = Int(inds[i])+1
        reqs[ind].val = reqvals[ind]
        reqs[ind].buffer = nothing
        indices[i] = ind
    end
    (indices, stats)
end

"""
    (indices, statuses) = Testsome!(reqs::Vector{Request})

Similar to [`Waitsome!`](@ref) except that it returns immediately: if no operations have
completed then `indices` and `statuses` will be empty.

# External links
$(_doc_external("MPI_Testsome"))
"""
function Testsome!(reqs::Vector{Request})
    count = length(reqs)
    reqvals = [req.val for req in reqs]
    outcnt = Ref{Cint}()
    inds = Array{Cint}(undef, count)
    stats = Array{Status}(undef, count)
    # int MPI_Testsome(int incount, MPI_Request array_of_requests[],
    #                  int *outcount, int array_of_indices[],
    #                  MPI_Status array_of_statuses[])
    @mpichk ccall((:MPI_Testsome, libmpi), Cint,
                  (Cint, Ptr{MPI_Request}, Ptr{Cint}, Ptr{Cint}, Ptr{Status}),
                  count, reqvals, outcnt, inds, stats)
    outcount = outcnt[]
    # This can happen if there were no valid requests
    if outcount == MPI_UNDEFINED
        outcount = 0
    end
    indices = Array{Int}(undef, outcount)
    for i in 1:outcount
        ind = Int(inds[i])+1
        reqs[ind].val = reqvals[ind]
        reqs[ind].buffer = nothing
        indices[i] = ind
    end
    (indices, stats)
end

"""
    Cancel!(req::Request)

Marks a pending [`Irecv!`](@ref) operation for cancellation (cancelling a [`Isend`](@ref)
is deprecated). Note that the request is not deallocated, and can still be queried using
the test or wait functions.

# External links
$(_doc_external("MPI_Cancel"))
"""
function Cancel!(req::Request)
    # int MPI_Cancel(MPI_Request *request)
    @mpichk ccall((:MPI_Cancel, libmpi), Cint, (Ptr{MPI_Request},), req)
    req.buffer = nothing
    nothing
end


text:
""
function parse_file(file::String)
    if endswith(file, ".m")
        pm_data = PowerModels.parse_matpower(file)
    elseif endswith(lowercase(file), ".raw")
        warn(LOGGER, "This feature is incomplete, and will currently only return RAW data from a PTI file")
        pm_data = PowerModels.parse_pti(file)
    else
        pm_data = parse_json(file)
    end

    check_network_data(pm_data)

    return pm_data
end

""
function check_network_data(data::Dict{String,Any})
    make_per_unit(data)
    check_transformer_parameters(data)
    check_voltage_angle_differences(data)
    check_thermal_limits(data)
    check_branch_directions(data)
    check_bus_types(data)
    check_dcline_limits(data)
    check_voltage_setpoints(data)
    check_cost_functions(data)
end


text:
abstract type AbstractCivecm end

eigvals(obj::AbstractCivecm) = eigvals(convert(VAR, obj))

function loglikelihood(obj::AbstractCivecm)
    nobs = size(obj.endogenous, 1) - obj.lags
    -nobs*(logdet(cholesky!(residualvariance(obj))) - size(obj.endogenous, 2)*(log(nobs) - 1 - log(2π)))/2
end
# loglikelihood(A::StridedMatrix) = -0.5*size(A, 1)*(2sum(log(abs(diag(qr(A).factors)/sqrt(size(A,1))))) - size(A, 2)*(log(size(A, 1)) - 1 - log(2π)))
loglikelihood(A::StridedMatrix) =
    -size(A, 1)*(logdet(Statistics.covzm(A, 1, corrected=false)) - size(A, 2)*(log(size(A, 1)) - 1 - log(2π)))/2

# aic(obj::Civecm) = 2*(npars(obj) - loglikelihood(obj))

function residualvariance(obj::AbstractCivecm)
    mresiduals = residuals(obj)
    mOmega = mresiduals'mresiduals / size(mresiduals, 1)
    return mOmega
end

normalitytest(obj::AbstractCivecm) = normalitytest(residuals(obj))

## LR test
mutable struct LRTest{T<:AbstractCivecm}
    H0::T
    HA::T
    value::Float64
    df::Int64
end

lrtest(obj0::AbstractCivecm, objA::AbstractCivecm, df::Integer) = LRTest(obj0, objA, 2*(loglikelihood(objA) - loglikelihood(obj0)), df)

function bootstrap(obj::LRTest, reps::Integer, simH0 = true)
    lrvals = Vector{Float64}(undef, reps)
    bootH0 = copy(obj.H0)
    bootHA = copy(obj.HA)
    bootHA.endogenous .= bootH0.endogenous
    bootVAR = simH0 ? convert(VAR, bootH0) : convert(VAR, bootHA)
    iT, p = size(obj.H0.Z0)
    # simResiduals = residuals(obj.H0)
    simResiduals = residuals(obj.HA)
    mbr = mean(simResiduals, dims=1)
    bi = Vector{Float64}(undef, iT)
    bootResiduals = similar(simResiduals)
    for i = 1:reps
        bi[:] = randn(iT)
        bootResiduals[:] = copy(simResiduals)
        for k = 1:p
            for j = 1:iT
                 bootResiduals[j,k] -= mbr[k]
                 bootResiduals[j,k] *= bi[j]
                 bootResiduals[j,k] += mbr[k]
            end
        end
        bootH0.endogenous[:] = simulate(bootVAR, bootResiduals)
        auxilliaryMatrices!(bootH0)
        auxilliaryMatrices!(bootHA)
        estimate!(bootH0)
        estimate!(bootHA)
        lrvals[i] = -2*(loglikelihood(bootH0) - loglikelihood(bootHA))
    end
    return lrvals
end


text:
using GeometricIntegrators
using Weave

GeometricIntegrators.set_config(:verbosity, 0)

Weave.set_chunk_defaults!(:echo => false, :results => "raw")

include("../src/lotka-volterra-2d.jl")

weave("../weave/lotka-volterra-2d-spark-pinternal.jmd",
         out_path = "src/lotka-volterra-2d",
         doctype = "github",
         mod = LotkaVolterra2dSingularSPARK)

weave("../weave/lotka-volterra-2d-spark-pmidpoint.jmd",
         out_path = "src/lotka-volterra-2d",
         doctype = "github",
         mod = LotkaVolterra2dSingularSPARK)

weave("../weave/lotka-volterra-2d-spark-psymmetric.jmd",
         out_path = "src/lotka-volterra-2d",
         doctype = "github",
         mod = LotkaVolterra2dSingularSPARK)

weave("../weave/lotka-volterra-2d-spark-plobatto.jmd",
         out_path = "src/lotka-volterra-2d",
         doctype = "github",
         mod = LotkaVolterra2dSingularSPARK)

weave("../weave/lotka-volterra-2d-spark-glvprk.jmd",
         out_path = "src/lotka-volterra-2d",
         doctype = "github",
         mod = LotkaVolterra2dSingularSPARK)

weave("../weave/lotka-volterra-2d-spark-glspark.jmd",
         out_path = "src/lotka-volterra-2d",
         doctype = "github",
         mod = LotkaVolterra2dSingularSPARK)

weave("../weave/lotka-volterra-2d-spark-lobspark.jmd",
         out_path = "src/lotka-volterra-2d",
         doctype = "github",
         mod = LotkaVolterra2dSingularSPARK)