text:
	println("	GRDINFO")
	G=gmt("grdmath", "-R0/10/0/10 -I1 5");
	r=gmt("grdinfo -C", G);
	@assert(r[1].data[1:1,1:10] == [0.0  10.0  0.0  10.0  5.0  5.0  1.0  1.0  11.0  11.0])
	r2=grdinfo(G, C=true, V=:q);
	@assert(r[1].data == r2[1].data)
	grdinfo(mat2grid(rand(4,4)));				# Test the doubles branch in grid_init
	grdinfo(mat2grid(rand(Float32,4,4)));		# Test the float branch in grid_init
	nx = 3; ny = 2; rang = [1, nx, 1, ny, 0.0,1]; x = collect(1.0:nx); y = collect(1.0:ny); inc = [1.0, 1];
	zz = rand(1:nx,ny,nx); z = zeros(Float32, ny+4, nx+4); z[3:end-2, 3:end-2] = zz;
	G = GMT.GMTgrid("","",0,rang, inc,0,NaN,"","","",x,y,[0.],deepcopy(collect(z')),"","","","","TRB",2)
	grdinfo(G);

	println("	GRD2CPT")
	G=gmt("grdmath", "-R0/10/0/10 -I2 X");
	C=grd2cpt(G);
	grd2cpt(G, cptname="lixo.cpt")

	# GRD2XYZ (It's tested near the end)
	#D=grd2xyz(G); # Use G of previous test
	gmtwrite("lixo.grd", G)
	D1=grd2xyz(G);
	D2=grd2xyz("lixo.grd");
	@assert(sum(D1[1].data) == sum(D2[1].data))

	println("	GRD2KML")
	G=gmt("grdmath", "-R0/10/0/10 -I1 X -fg");
	grd2kml(G, I="+", N="NUL", V="q", Vd=dbg2)

	println("	GRDBLEND")
	G3=gmt("grdmath", "-R5/15/0/10 -I1 X Y -Vq");
	G2=grdblend(G,G3);

	println("	GRDCLIP")
	G2=grdclip(G,above="5/6", low=[2 2], between=[3 4.5 4]);	 # Use G of previous test
	@test_throws ErrorException("Wrong number of elements in S option") G2=grdclip(G,above="5/6", low=[2], between=[3 4 4.5]);
	@test_throws ErrorException("OPT_S: argument must be a string or a two elements array.") G2=grdclip(G,above=5, low=[2 2]);

	println("	GRDCONTOUR")
	G = gmt("grdmath -R-15/15/-15/15 -I0.3 X Y HYPOT DUP 2 MUL PI MUL 8 DIV COS EXCH NEG 10 DIV EXP MUL =");
	C = grdcontour(G, C="+0.7", D=[]);
	@assert((size(C[1].data,1) == 21) && abs(-0.6 - C[1].data[1,1]) < 1e-8)
	# Do the same but write the file on disk first
	gmt("write lixo.grd", G)
	GG = gmt("read -Tg lixo.grd");
	C = grdcontour("lixo.grd", C="+0.7", D=[]);
	@assert((size(C[1].data,1) == 21) && abs(-0.6 - C[1].data[1,1]) < 1e-8)
	r = grdcontour("lixo.grd", cont=10, A=(int=50,labels=(font=7,)), G=(dist="4i",), L=(-1000,-1), W=((contour=1,pen="thinnest,-"), (annot=1, pen="thin,-")), T=(gap=("0.1i","0.02i"),), Vd=dbg2);
	@test startswith(r, "grdcontour lixo.grd  -JX" * split(GMT.def_fig_size, '/')[1] * "/0" * " -Baf -BWSen -L-1000/-1 -A50+f7 -Gd4i -T+d0.1i/0.02i -Wcthinnest,- -Wathin,- -R-15/15/-15/15 -C10")
	r = grdcontour("lixo.grd", A="50+f7p", G="d4i", W=((contour=1,pen="thinnest,-"), (annot=1, pen="thin,-")), Vd=dbg2);
	@test startswith(r, "grdcontour lixo.grd  -JX" * split(GMT.def_fig_size, '/')[1] * "/0" * " -Baf -BWSen -A50+f7p -Gd4i -Wcthinnest,- -Wathin,-")
	G = GMT.peaks()
	cpt = makecpt(T="-6/8/1");
	grdcontour(G, axis="a", fmt="png", color=cpt, pen="+c", X=1, Y=1, N=true, U=[])
	grdcontour!(G, axis="a", color=cpt, pen="+c", X=1, Y=1, N=true, Vd=dbg2)
	grdcontour!("", G, axis="a", color=cpt, pen="+c", X=1, Y=1, N=cpt, Vd=dbg2)

	println("	GRDCUT")
	G=gmt("grdmath", "-R0/10/0/10 -I1 X Y MUL");
	G2=grdcut(G, limits=[3 9 2 8]);
	G2=grdcut("lixo.grd", limits=[3 9 2 8], V=:q);	# lixo.grd was written above in the gmtwrite test
	G2=grdcut(data="lixo.grd", limits=[3 9 2 8], V=:q);
	G2=grdcut(data=G, limits=[3 9 2 8]);

	# GRDEDIT
	grdedit(G, C=true);

	println("	GRDFFT")
	G2=grdfft(G, upward=800); 	# Use G of previous test
	G2=grdfft(G, G, E=[]);

	println("	GRDFIL")
	G1=grdmask([3 3], R="0/6/0/6", I=1, N="10/NaN/NaN", S=0);
	G2=grdfill(G, algo=:n);

	println("	GRDFILTER")
	G2=grdfilter(G, filter="m600", distflag=4, inc=0.5); # Use G of previous test

	println("	GRDGRADIENT")
	G2=grdgradient(G, azim="0/270", normalize="e0.6");
	G2=grdgradient(G, azim="0/270", normalize="e0.6", Q=:save, Vd=dbg2);

	println("	GRDHISTEQ")
	G2 = grdhisteq(G, gaussian=[]);	# Use G of previous test

	if (GMTver > v"6.1.1")
		println("	GRDINTERPOLATE")
		G = grdinterpolate("cube.nc", T=4)
		C = grdinterpolate("cube.nc", T="3/4/0.25");
		grdinfo(C)
	end

	println("	GRDLANDMASK")
	G2 = grdlandmask(R="-10/4/37/45", res=:c, inc=0.1);
	G2 = grdlandmask("-R-10/4/37/45 -Dc -I0.1");			# Monolithitc

	println("	GRDMASK")
	G2 = grdmask([10 20; 40 40; 70 20; 10 20], R="0/100/0/100", out_edge_in=[100 0 0], I=2);

	println("	GRDPASTE")
	G3 = grdmath("-R10/20/0/10 -I1 X");
	G2 = grdpaste(G,G3);

	println("	GRDPROJECT")
	# GRDPROJECT	-- Works but does not save projection info in header
	G2 = grdproject(G, proj="u29/1:1", F=[], C=[]); 		# Use G of previous test
	G2 = grdproject("-Ju29/1:1 -F -C", G);					# Monolithic

	println("	GRDSAMPLE")
	G2 = grdsample(G, inc=0.5);		# Use G of previous test

	println("	GRDTREND")
	G  = gmt("grdmath", "-R0/10/0/10 -I1 X Y MUL");
	G2 = grdtrend(G, model=3);
	W = mat2grid(ones(Float32, size(G.z,1), size(G.z,2)));
	W = mat2grid(rand(16,16), x=11:26, y=1:16);
	mat2grid([3 4 5; 1 2 5; 5 5 5], reg=:pixel, x=[1 3], y=[1 3]);
	G2 = grdtrend(G, model=3, diff=[], trend=true);
	#G2 = grdtrend(G, model="3+r", W=W);
	G2 = grdtrend(G, model="3+r", W=(W,0), Vd=dbg2);

	println("	GRDTRACK")
	#G = gmt("grdmath -R-15/15/-15/15 -I0.3 X Y HYPOT DUP 2 MUL PI MUL 8 DIV COS EXCH NEG 10 DIV EXP MUL =");
	G = gmt("grdmath -R-2/2/-2/2 -I1 1");
	gmtwrite("lixo.grd", G)
	D = grdtrack([0 0], G);
	@assert(D[1].data == [0.0 0 1])
	D = grdtrack([0 0], G="lixo.grd");
	@assert(D[1].data == [0.0 0 1])
	D = grdtrack("lixo.grd", [0 0]);
	D = grdtrack(G, [0 0]);
	D = grdtrack([0 0], G=G);
	D = grdtrack([0 0], G=(G,G));
	@assert(D[1].data == [0.0 0 1 1])

	println("	GRDVECTOR")
	G = gmt("grdmath -R-2/2/-2/2 -I0.1 X Y R2 NEG EXP X MUL");
	dzdy = gmt("grdmath ? DDY", G);
	dzdx = gmt("grdmath ? DDX", G);
	grdvector(dzdx, dzdy, I=0.2, vector=(len=0.25, stop=1, norm=0.65, shape=0.5), G=:black, W="1p", S=12)
	grdvector!(dzdx, dzdy, I=0.2, vector=(len=0.25, stop=1, norm=0.65, shape=0.5), W="1p", S=12, Vd=dbg2)
	r = grdvector!("",dzdx, dzdy, I=0.2, vector=(len=0.25, stop=1, norm=0.65), W="1p", S=12, Vd=dbg2);
	@test startswith(r, "grdvector  -R -J -I0.2 -S12 -Q0.25+e+n0.65 -W1p")
	r = grdvector!("", 1, 2, I=0.2, vec="0.25+e+n0.66", W=1, S=12, Vd=dbg2);
	@test startswith(r, "grdvector  -R -J -I0.2 -S12 -Q0.25+e+n0.66 -W1")

	println("	GRDVOLUME")
	grdvolume(G);


text:
using BlockDiagonalMatrices
using BlockDiagonalMatrices: isequal_blocksizes
using Random
using Test

@testset "blockdiagonal.jl" begin
    rng = MersenneTwister(123456)
    N1, N2, N3 = 3, 4, 5
    N = N1 + N2 + N3
    b1 = BlockDiagonal([rand(rng, N1, N1), rand(rng, N2, N2), rand(rng, N3, N3)])
    b2 = BlockDiagonal([rand(rng, N1, N1), rand(rng, N3, N3), rand(rng, N2, N2)])
    b3 = BlockDiagonal([rand(rng, N1, N1), rand(rng, N2, N2), rand(rng, N2, N2)])
    A = rand(rng, N, N + N1)
    B = rand(rng, N + N1, N + N2)
    A′, B′ = A', B'
    a = rand(rng, N)
    b = rand(rng, N + N1)

    @testset "AbstractArray" begin
        X = rand(2, 2); Y = rand(3, 3)

        @test size(b1) == (N, N)
        @test size(b1, 1) == N && size(b1, 2) == N

        eqs = []
        for i in 1:size(b1, 1)
            for j in 1:size(b1, 2)
                push!(eqs, b1[i, j] ≈ Matrix(b1)[i, j])
            end
        end
        @test all(eqs)

        @testset "BlockDiagonal does not copy" begin
            Bxy = BlockDiagonal([X, Y])
            X[1] = 1.1
            @test Bxy[1] == 1.1
            Bxy2 = BlockDiagonal(Bxy)
            Bxy[2] = 2.2
            @test X[2] == Bxy[2] == Bxy2[2] == 2.2
        end

        @testset "parent" begin
            @test parent(b1) isa Vector{<:AbstractMatrix}
            @test parent(BlockDiagonal([X, Y])) == [X, Y]
        end

        @testset "similar" begin
            @test similar(b1) isa BlockDiagonal
            @test size(similar(b1)) == size(b1)
            @test size.(blocks(similar(b1))) == size.(blocks(b1))
        end

        @testset "zero" begin
            @test typeof(zero(b1)) == typeof(b1)
            @test size.(blocks(zero(b1))) == size.(blocks(b1))
            @test zero(b1) == 0 * b1
        end

        @testset "setindex!" begin
            X = BlockDiagonal([rand(Float32, 5, 5), rand(Float32, 3, 3)])
            X[10] = Int(10)
            @test X[10] === Float32(10.0)
            X[3, 3] = Int(9)
            @test X[3, 3] === Float32(9.0)
            # Should not allow setting value outside on-diagonal blocks to non-zero
            @test_throws ArgumentError X[1, 7] = 1
        end
    end  # AbstractArray

    @testset "isequal_blocksizes" begin
        @test isequal_blocksizes(b1, b1) == true
        @test isequal_blocksizes(b1, similar(b1)) == true
        @test isequal_blocksizes(b1, b2) == false
    end

    @testset "blocks size" begin
        B = BlockDiagonal([rand(3, 3), rand(4, 4)])
        @test nblocks(B) == 2
        @test blocksizes(B) == [(3, 3), (4, 4)]
        @test blocksize(B, 2) == blocksizes(B)[2] == blocksize(B, 2, 2)
    end

    @testset "Equality" begin
        # Equality
        @test b1 == b1
        @test b1 == Matrix(b1)
        @test Matrix(b1) == b1

        # Inequality
        @test b1 != b2
        @test b1 != Matrix(b2)
        @test Matrix(b1) != b2

        # Approximate equality
        @test b1 ≈ b1
        @test Matrix(b1) ≈ b1
        @test b1 ≈ Matrix(b1)
    end  # Equality

    @testset "Non-Square Matrix" begin
	A1 = ones(2, 4)
	A2 = 2 * ones(3, 2)
	B1 = BlockDiagonal([A1, A2])
        B2 = [A1 zeros(2, 2); zeros(3, 4) A2]

	@test B1 == B2
	# Dimension check
	@test sum(size.(B1.blocks, 1)) == size(B2, 1)
	@test sum(size.(B1.blocks, 2)) == size(B2, 2)
    end  # Non-Square Matrix

    @testset "copy" begin
        bc = similar(b1)

        copy!(bc, b1)
        @test bc == b1

        c = copy(b1)
        @test c == b1

        @test_throws DimensionMismatch copy!(b2, b1)
    end
end

ext: jl


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    AnisotropicBall(radii, angles; convention=TaitBryanExtr)

An anisotropic ball with `radii` and `angles` according to the rotation `convention`.

- For 2D balls (i.e. ellipses), there are two radii and one rotation angle.
- For 3D balls (i.e. ellipsoids), there are three radii and three rotation angles.

The list of available conventions can be found with `subtypes(RotationConvention)`.
"""
struct AnisotropicBall{V,A,C,M} <: MetricBall
  # input fields
  radii::V
  angles::A
  convention::C

  # state fields
  metric::M
end

function AnisotropicBall(radii::V, angles::A; convention::C=TaitBryanExtr) where {V,A,C}
  Dim, nangles = length(radii), length(angles)
  valid = (Dim == 3 && nangles == 3) || (Dim == 2 && nangles == 1)
  @assert valid "invalid number of radii/angles"

  # invert radii if necessary
  invert = mainaxis(convention) == :Y
  ranges = invert ? [radii[i] for i in reverse(1:Dim,1,2)] : radii

  # scaling matrix
  Λ = Diagonal(SVector{Dim}(one(eltype(ranges))./ranges.^2))

  # rotation matrix
  R = rotmat(angles, convention)

  # ellipsoid metric
  metric = Mahalanobis(Symmetric(R'*Λ*R))

  AnisotropicBall{V,A,C,typeof(metric)}(radii, angles, convention, metric)
end

metric(ball::AnisotropicBall) = ball.metric

radii(ball::AnisotropicBall) = ball.radii

function Base.show(io::IO, ball::AnisotropicBall)
  r = ball.radii
  a = ball.angles
  c = ball.convention
  T = eltype(r)
  print(io, "AnisotropicBall{$T}($r, $a, $c)")
end


text:
immutable Categorical <: DiscreteUnivariateDistribution
    prob::Vector{Float64}
    drawtable::DiscreteDistributionTable
    function Categorical{T <: Real}(p::Vector{T})
        if length(p) <= 1
            error("Categorical: there must be at least two categories")
        end
        sump = 0.0
        for i in 1:length(p)
            if p[i] < 0.0
                error("Categorical: probabilities must be non-negative")
            end
            sump += p[i]
        end
        for i in 1:length(p)
            p[i] /= sump
        end
        new(p, DiscreteDistributionTable(p))
    end
end

function Categorical(d::Integer)
    if d <= 1
        error("d must be greater than 1")
    end
    prob = Array(Float64, d)
    fill!(prob, 1.0 / d)
    Categorical(prob)
end

function cdf(d::Categorical, x::Integer)
    if !insupport(d, x)
        error("$x is not in the support")
    else
        p = 0.0
        for i in 1:x
            p += d.prob[i]
        end
        return p
    end
end

entropy(d::Categorical) = pventropy(d.prob)

function insupport(d::Categorical, x::Real)
    return isinteger(x) && 1 <= x <= length(d.prob) && d.prob[x] != 0.0
end

function mean(d::Categorical)
    s = 0.0
    for i in 1:length(d.prob)
        s += i * d.prob[i]
    end
    return s
end

function median(d::Categorical)
    p, n = 0.0, length(d.prob)
    i = 0
    while p < 0.5 && i <= n
        i += 1
        p += d.prob[i]
    end
    return i
end

function mgf(d::Categorical, t::AbstractVector)
    s = 0.0
    for i in 1:length(d.prob)
        s += d.prob[i] * exp(t[i])
    end
    return s
end

function cf(d::Categorical, t::AbstractVector)
    s = 0.0 + 0.0im
    for i in 1:length(d.prob)
        s += d.prob[i] * exp(im * t[i])
    end
    return s
end

modes(d::Categorical) = [indmax(d.prob)]

pdf(d::Categorical, x::Real) = !insupport(d, x) ? 0.0 : d.prob[x]

rand(d::Categorical) = draw(d.drawtable)

var(d::Categorical) = var(d, mean(d))

function var(d::Categorical, m::Number)
    s = 0.0
    for i in 1:length(d.prob)
        s += (i - m)^2 * d.prob[i]
    end
    return s
end


text:
export minimize!, solve!

# Eventually, this can take in arguments to specify solving method.
# For now, backslash is the only supported method.
function solve!(problem_type::Symbol, objective::SumSquaresExpr, constraints::Array{EqConstraint, 1})
  p = Problem(problem_type, objective, constraints)
  # return backslash_solve!(p)
  backslash_solve!(p)
  if p.status == "KKT singular"
    error("KKT system is singular. No unique solution")
  elseif p.status == "solved"
    return p.optval
  else
    error("Unrecognized problem status")
  end
end

minimize!(objective::SumSquaresExpr, constraints::Array{EqConstraint, 1}) = solve!(:minimize, objective, constraints)
minimize!(objective::SumSquaresExpr, constraint::EqConstraint) = minimize!(objective, [constraint])
minimize!(objective::SumSquaresExpr, constraints::EqConstraint...) = minimize!(objective, [constraints...])
minimize!(objective::SumSquaresExpr) = minimize!(objective, EqConstraint[])

solve!(constraints::Array{EqConstraint, 1}) = solve!(:solve, sum_squares(), constraints)
solve!(constraint::EqConstraint) = solve!([constraint])
solve!(constraints::EqConstraint...) = solve!([constraints...])


text:
module Day13

using AdventOfCode2019

function day13(input::String = readInput(joinpath(@__DIR__, "..", "data", "day13.txt")))
    program = parse.(Int, split(input, ","))
    output = AdventOfCode2019.IntCode.run_program!(copy(program), nothing, nothing, nothing)
    part1 = sum(output[3:3:length(output)] .== 2)

    part2 = _play_game(program)
    return [part1, part2]
end

function _play_game(program::Array{Int,1})
    program[1] = 2
    input = Channel{Int}(10)
    output = Channel{Int}(10)
    t = @task AdventOfCode2019.IntCode.run_program!(copy(program), input, output, nothing)
    schedule(t)

    currentScore = 0
    paddlePos = (-1, 0)
    ballPos = (-1, 0)
    while !istaskdone(t)
        x = take!(output)
        y = take!(output)
        id = take!(output)
        if id == 4
            ballPos = (x, y)
            if ballPos[1] < paddlePos[1]
                put!(input, -1)
            elseif ballPos[1] > paddlePos[1]
                put!(input, 1)
            else
                put!(input, 0)
            end
        elseif id == 3
            paddlePos = (x, y)
        end

        if x == -1 && y == 0
            currentScore = id
        end
    end
    close(input)
    close(output)
    return currentScore
end

end # module


text:
"""
**FUNÇÃO ÁREA LATERAL**
Sintaxe:
arealateral(raio, altura_cilindro)
"""

function arealateral(raio, altura)
    return  2*pi*raio*altura
end


text:
@testset "Views" begin
  @testset "Domain" begin
    d = RegularGrid(10, 10)
    X = coordinates(d)

    v1 = view(d, 1:10)
    v2 = view(d, 1:10)
    v3 = view(d, 1:20)
    @test v1 == v2
    @test v1 != v3
    @test v2 != v3

    v = view(d, 1:10)
    @test nelms(v) == 10
    @test coordinates(v) == X[:,1:10]
    @test collect(v) isa PointSet

    @test sprint(show, v) == "10 DomainView{Float64,2}"
    @test sprint(show, MIME"text/plain"(), v) == "10 DomainView{Float64,2}\n 0.0  1.0  2.0  3.0  4.0  5.0  6.0  7.0  8.0  9.0\n 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0"
  end

  @testset "Data" begin
    d = georef((z=rand(100), w=rand(100)))
    T = values(d)
    X = coordinates(d)

    v1 = view(d, 1:10)
    v2 = view(d, 1:10)
    v3 = view(d, 1:20)
    @test v1 == v2
    @test v1 != v3
    @test v2 != v3

    v = view(d, 1:10)
    @test nelms(v) == 10
    @test coordinates(v) == X[:,1:10]
    @test collect(v) isa SpatialData

    @test sprint(show, v) == "10 DataView{Float64,1}"
    @test sprint(show, MIME"text/plain"(), v) == "10 DomainView{Float64,1}\n  variables\n    └─w (Float64)\n    └─z (Float64)"

    d = georef((z=rand(100), w=rand(100)))
    v1 = view(d, shuffle(1:100))
    c1 = collect(v1)
    v2 = view(v1, [:z])
    c2 = view(c1, [:z])
    @test values(v2) == values(c2)
  end
end


text:
using Xpress, MathOptInterface, Test

const MOI = MathOptInterface
const MOIT = MathOptInterface.Test

const MOIT = MOI.Test
const MOIU = MOI.Utilities

const OPTIMIZER = Xpress.Optimizer()
const BRIDGED_OPTIMIZER = MOI.Bridges.full_bridge_optimizer(
    Xpress.Optimizer(), Float64)

const CONFIG = MOIT.TestConfig()

@testset "SolverName" begin
    @test MOI.get(OPTIMIZER, MOI.SolverName()) == "Xpress"
end

@testset "supports_default_copy_to" begin
    @test MOIU.supports_default_copy_to(OPTIMIZER, true)
end

@testset "Unit Tests" begin
    MOIT.basic_constraint_tests(OPTIMIZER, CONFIG;
            exclude = [
                (MOI.SingleVariable, MOI.Integer),
                (MOI.SingleVariable, MOI.EqualTo{Float64}),
                (MOI.SingleVariable, MOI.Interval{Float64}),
                (MOI.SingleVariable, MOI.GreaterThan{Float64}),
                (MOI.SingleVariable, MOI.LessThan{Float64}),
            ]
        )
        MOIT.unittest(OPTIMIZER, CONFIG, [
            "solve_qp_edge_cases",    # tested below
            "solve_qcp_edge_cases"    # tested below
        ])
    #    MOIT.modificationtest(BRIDGED_OPTIMIZER, CONFIG)
end
#=
@testset "Linear tests" begin
    @testset "Default Solver"  begin
        MOIT.contlineartest(OPTIMIZER, MOIT.TestConfig(basis = true), [
            # This requires an infeasiblity certificate for a variable bound.
            "linear12"
        ])
    end
    @testset "No certificate" begin
        MOIT.linear12test(OPTIMIZER, MOIT.TestConfig(infeas_certificates=false))
    end
end

@testset "Quadratic tests" begin
    MOIT.contquadratictest(OPTIMIZER, MOIT.TestConfig(atol=1e-3, rtol=1e-3), [
        "ncqcp"  # Gurobi doesn't support non-convex problems.
    ])
end

@testset "Conic tests" begin
    MOIT.lintest(OPTIMIZER, CONFIG)
    MOIT.soctest(OPTIMIZER, MOIT.TestConfig(duals = false, atol=1e-3), ["soc3"])
    MOIT.soc3test(
        OPTIMIZER,
        MOIT.TestConfig(duals = false, infeas_certificates = false, atol = 1e-3)
    )
    MOIT.rsoctest(OPTIMIZER, MOIT.TestConfig(duals = false, atol=1e-3))
    MOIT.geomeantest(OPTIMIZER, MOIT.TestConfig(duals = false, atol=1e-3))
end

@testset "Integer Linear tests" begin
    MOIT.intlineartest(OPTIMIZER, CONFIG, [
        # Indicator sets not supported.
        "indicator1", "indicator2", "indicator3"
    ])
end

@testset "ModelLike tests" begin
    @testset "default_objective_test" begin
        MOIT.default_objective_test(OPTIMIZER)
    end

    @testset "default_status_test" begin
        MOIT.default_status_test(OPTIMIZER)
    end

    @testset "nametest" begin
        MOIT.nametest(OPTIMIZER)
    end

    @testset "validtest" begin
        MOIT.validtest(OPTIMIZER)
    end

    @testset "emptytest" begin
        MOIT.emptytest(OPTIMIZER)
    end

    @testset "orderedindicestest" begin
        MOIT.orderedindicestest(OPTIMIZER)
    end
    @testset "copytest" begin
        OPTIMIZER_2 = Xpress.Optimizer()
        MOIT.copytest(OPTIMIZER, OPTIMIZER_2 )
    end
end

=#


text:
#!/usr/bin/env julia

# Given some model parameters, plot the channel maps and
# integrated spectrum

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "--config"
    help = "a YAML configuration file"
    default = "config.yaml"
    "--linear"
    help = "Plot the regular linear channel maps."
    action = :store_true
    "--log"
    help = "Plot the channel maps with log stretch."
    action = :store_true
    "--blur"
    help = "Plot the channel maps, convolved with a fake beam."
    action = :store_true
    "--spectrum"
    help = "Plot the spatially-integrated spectrum."
    action = :store_true
end

parsed_args = parse_args(ARGS, s)

import YAML
config = YAML.load(open(parsed_args["config"]))

using DiskJockey.constants
using DiskJockey.image
using DiskJockey.model
using HDF5

import PyPlot.plt
using LaTeXStrings
import Images

species = config["species"]
transition = config["transition"]
lam0 = lam0s[species*transition]
model = config["model"]

# Read the incination, position angle, and (grid) radius of the disk. Plot ellipses.
pars = config["parameters"]
r_c = pars["r_c"]
incl = pars["incl"]
PA = pars["PA"]
dpc = pars["dpc"]

# from matplotlib.patches import Ellipse
#     ax[:add_artist](PyPlot.matplotlib[:patches][:Ellipse](xy=xy, width=BMIN, height=BMAJ, angle=BPA, facecolor="0.8", linewidth=0.2))
width = r_c/dpc # [arcsec]
height = width * cosd(incl) # [arcsec]

# cmap = plt[:get_cmap]("viridis")
# cmap = plt[:get_cmap]("inferno")
cmap = plt[:get_cmap]("plasma")

# First contour is at 3 sigma, and then contours go up (or down) in multiples of spacing
function get_levels(rms::Float64, vmax::Float64, spacing=3)
    levels = Float64[]

    val = 3 * rms
    while (val < vmax)
        append!(levels, [val])
        val += rms * spacing
    end

    return levels
end

# function plot_beam(ax, BMAJ, BMIN, xy=(1,-1))
#     BMAJ = 3600. * header["BMAJ"] # [arcsec]
#     BMIN = 3600. * header["BMIN"] # [arcsec]
#     BPA =  header["BPA"] # degrees East of North
#     # from matplotlib.patches import Ellipse
#     ax[:add_artist](PyPlot.matplotlib[:patches][:Ellipse](xy=xy, width=BMIN, height=BMAJ, angle=BPA, facecolor="0.8", linewidth=0.2))
# end

# Plot the channel maps using sky convention
"""Plot the channel maps using the sky convention. If log is true, plot intensity using
a log scale."""
function plot_chmaps(img::image.SkyImage; log=false, contours=true, fname="channel_maps_sky.png")

    if log
        ldata = log10(img.data + 1e-99)
        vvmax = maximum(ldata)
        norm = PyPlot.matplotlib[:colors][:Normalize](vvmax - 6, vvmax)
    else
        vvmax = maxabs(img.data)
        # println(vmin, " ", vmax, " ", vvmax)
        norm = PyPlot.matplotlib[:colors][:Normalize](0, vvmax)

        if contours
            levels = get_levels(rms, vvmax)
        end
    end

    (im_ny, im_nx) = size(img.data)[1:2] # y and x dimensions of the image

    # Image needs to be flipped along RA dimension
    ext = (img.ra[end], img.ra[1], img.dec[1], img.dec[end])

    # Figure out how many plots we'll have.
    ncols = 8
    nrows = ceil(Int, nlam/ncols)

    xx = 1.5 * 9
    dx = 1.5
    yy = (nrows + 1) * 1.5
    dy = 1.5

    fig, ax = plt[:subplots](nrows=nrows, ncols=ncols, figsize=(xx, yy))

    for row=1:nrows
        for col=1:ncols
            iframe = col + (row - 1) * ncols

            if col != 1 || row != nrows
                ax[row, col][:xaxis][:set_ticklabels]([])
                ax[row, col][:yaxis][:set_ticklabels]([])
            else
                ax[row, col][:set_xlabel](L"$\Delta \alpha$ ('')")
                ax[row, col][:set_ylabel](L"$\Delta \delta$ ('')")
            end

            if iframe > nlam
                # Plot a blank square if we run out of channels
                ax[row, col][:imshow](zeros((im_ny, im_nx)), cmap=cmap, vmin=0, vmax=20, extent=ext, origin="lower")

            else
                #Flip the frame for Sky convention
                frame = flipdim(img.data[:,:,iframe], 2)

                if log
                    frame += 1e-15 #Add a tiny bit so that we don't have log10(0)
                    lframe = log10(frame)
                    im = ax[row, col][:imshow](lframe, extent=ext, interpolation="none", origin="lower", cmap=cmap, norm=norm)

                else
                    im = ax[row, col][:imshow](frame, extent=ext, interpolation="none", origin="lower", cmap=cmap, norm=norm)

                    if contours
                        ax[row, col][:contour](frame, origin="lower", colors="k", levels=levels, extent=ext, linestyles="solid", linewidths=0.2)
                    end

                end

                ax[row, col][:add_artist](PyPlot.matplotlib[:patches][:Ellipse]((0,0), width, height, PA, linewidth=0.15, facecolor="none", edgecolor="w"))

                if iframe==1
                    # Plot the colorbar
                    cax = fig[:add_axes]([(xx - 0.35 * dx)/xx, (yy - 1.5 * dy)/yy, (0.1 * dx)/xx, dy/yy])
                    cbar = fig[:colorbar](mappable=im, cax=cax)

                    cbar[:ax][:tick_params](labelsize=6)
                    fig[:text](0.99, (yy - 1.7 * dy)/yy, "Jy/beam", size=8, ha="right")
                end

                ax[row, col][:annotate](@sprintf("%.1f", vels[iframe]), (0.1, 0.8), xycoords="axes fraction", size=8)
            end

        end
    end

    fig[:subplots_adjust](hspace=0.00, wspace=0.00, top=(yy - 0.5 * dy)/yy, bottom=(0.5 * dy)/yy, left=(0.5 * dx)/xx, right=(xx - 0.5 * dy)/xx)

    plt[:savefig](fname, dpi=600)

end

# Plot the spatially-integrated spectrum
function plot_spectrum(img::image.SkyImage; fname="spectrum.png")

    fig = plt[:figure]()
    ax = fig[:add_subplot](111)

    spec = imToSpec(img)

    ax[:plot](vels, spec[:,2], ls="steps-mid")

    ax[:set_ylabel](L"$f_\nu$ [Jy]")
    ax[:set_xlabel](L"$v$ [km/s]")

    fig[:subplots_adjust](left=0.15, bottom=0.15, right=0.85)

    plt[:savefig](fname)

    # Calculate integrated line intensity.
    tot = image.integrateSpec(spec, lam0)
    println("Total line flux ", tot, " Jy / km / s")
end

pars = convert_dict(config["parameters"], config["model"])

# Assume that image has already been generated by synthesize_model.jl
im = imread()
skim = imToSky(im, pars.dpc)

# Do the velocity conversion here for plot labels
global nlam = length(skim.lams)
# convert wavelengths to velocities
global vels = c_kms * (skim.lams .- lam0)/lam0

# Make sure we are always plotting channels from blueshift to redshift.
if vels[2] < vels[1]
  println("Flipping channel maps to plot in order of increasing velocity.")
  skim.data = skim.data[:,:,end:-1:1]
  vels = vels[end:-1:1]
end

if parsed_args["linear"]
    plot_chmaps(skim, fname="chmaps_linear.png", contours=false)
end

if parsed_args["log"]
    plot_chmaps(skim, fname="chmaps_log.png", log=true, contours=false)
end

if parsed_args["blur"]
    beam = config["beam"]
    rms = beam["rms"] # Jy/beam
    BMAJ = beam["BMAJ"]/2 # semi-major axis [arcsec]
    BMIN = beam["BMIN"]/2 # semi-minor axis [arcsec]
    BAVG = (BMAJ + BMIN)/2
    BPA = beam["BPA"] # position angle East of North [degrees]

    println("Beam sigma ", BAVG, " [arcsec]")

    arcsec_ster = (4.25e10)
    # Convert beam from arcsec^2 to Steradians
    global rms = rms/(pi * BMAJ * BMIN) * arcsec_ster

    println("bluring maps")
    sk_blur = blur(skim, [BAVG, BAVG])

    println("Plotting blured maps")
    plot_chmaps(sk_blur, fname="chmaps_blur.png", log=false, contours=true)
end

if parsed_args["spectrum"]
    plot_spectrum(skim)
end

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
module Mathes

import Statistics: mean
import LinearAlgebra: svd, diagm, diag

"""
    pca(data)
perform PCA using SVD
inputs:
    - data: M x N matrix of input data. (M dimensions, N trials)
outputs:
    - PC: each column is a principle component
    - V: M x 1 matrix of variances
"""
function pca(data::Array{T,2}) where T
    X = data .- mean(data, dims=2)
    Y = X' ./ sqrt(T(size(X,2)-1))
    U,S,PC = svd(Y)
    S = diagm(0=>S)
    V = S .* S

    # find the least variance vector
    indexList = sortperm(diag(V); rev=true)

    PCs = map(x->PC[:,x], indexList)
    return PCs, diag(V)[indexList]
end

end # end of module


text:
using Pkg;
Pkg.activate("");
using CUTEst, Dates, NLPModels, SolverBenchmark

using ARCTR, JSOSolvers, NLPModelsIpopt

nmax = 10000
problems = readlines("list_problems_$nmax.dat")
cutest_problems = (CUTEstModel(p) for p in problems)

max_time = 300.0 # 3600.0
max_ev = typemax(Int)
max_iter = typemax(Int)
atol = 1e-5
rtol = 1e-6

solvers = Dict(
    :ARCqKOp10201z025 =>
        nlp -> ARCqKOp(
            nlp,
            verbose = false,
            atol = atol,
            rtol = rtol,
            max_time = max_time,
            max_iter = max_iter,
            shifts = 10.0 .^ (collect(-10.0:1.0:20.0)),
            ζ = 0.25,
        ),
    :ARCqKOp10201z05 =>
        nlp -> ARCqKOp(
            nlp,
            verbose = false,
            atol = atol,
            rtol = rtol,
            max_time = max_time,
            max_iter = max_iter,
            shifts = 10.0 .^ (collect(-10.0:1.0:20.0)),
            ζ = 0.5,
        ),
    :ARCqKOp10201z075 =>
        nlp -> ARCqKOp(
            nlp,
            verbose = false,
            atol = atol,
            rtol = rtol,
            max_time = max_time,
            max_iter = max_iter,
            shifts = 10.0 .^ (collect(-10.0:1.0:20.0)),
            ζ = 0.75,
        ),
    :ARCqKOp10201z1 =>
        nlp -> ARCqKOp(
            nlp,
            verbose = false,
            atol = atol,
            rtol = rtol,
            max_time = max_time,
            max_iter = max_iter,
            shifts = 10.0 .^ (collect(-10.0:1.0:20.0)),
            ζ = 1.0,
        ),
)
rmv_problems = [
    "SBRYBND",
    "GENHUMPS",
    "NELSONLS",
    "NONMSQRT",
    "FMINSRF2",
    "JIMACK",
    "EIGENCLS",
    "INDEF",
    "FLETCBV3",
    "FMINSURF",
    "FLETCHBV",
    "PARKCH",
    "EIGENBLS",
]
stats = bmark_solvers(
    solvers,
    cutest_problems,
    skipif = prob -> (prob.meta.name in rmv_problems),
)

using JLD2
@save "$(today())_$(prod(String.(keys(solvers)) .* "_"))cutest_$(string(length(problems)))_$(nmax).jld2" stats


text:
using Query
using DataFrames
using Test

@testset "Pipe Syntax" begin

    df = DataFrame(a=[1,2,3], b=[3.,2.,1.], c=["a", "b", "c"])

    df2 = df |> @query(i, begin
        @where i.a>2
        @select {i.c, i.b}
    end) |> DataFrame

    @test df2 isa DataFrame
    @test size(df2) == (1,2)
    @test df2[1,:c] == "c"
    @test df2[1,:b] == 1.
end


text:
"""
    Fermi.HartreeFock

Module for running Hartree--Fock computations in Fermi.

# Methods

    > Fermi.HartreeFock.RHF
    > Fermi.HartreeFock.UHF
"""
module HartreeFock
# Import Fermi basics
using LinearAlgebra: hermitian, vcat
using Fermi
using Fermi.Options
using Fermi.Integrals
using Fermi.Orbitals

function hf_header()
    output(repeat("=",80))
    output("|{:33}{:^12}{:33}|", "", "Hartree-Fock", "")
    output("|{:34}{:^9}{:34}|", "", "Module  by","")
    output("|{:25}{:^28}{:25}|", "", "G.J.R Aroeira and M.M. Davis", "")
    output(repeat("=",80))
end

function uhf_header()
    output(repeat("=",80))
    output("|{:33}{:^12}{:33}|", "", "Hartree-Fock", "")
    output("|{:34}{:^9}{:34}|", "", "Module  by","")
    output("|{:17}{:^44}{:17}|", "", "G.J.R Aroeira, M.M. Davis, and S.M. Goodlett", "")
    output(repeat("=",80))
end

"""
    Fermi.HartreeFock.AbstractHFWavefunction

Abstract type common to all Hartree-Fock wave functions.

Struct tree

**AbstractHFWavefunction** <: AbstractWavefunction
"""
abstract type AbstractHFWavefunction <: Fermi.AbstractWavefunction end

# Different Hartree-Fock methods are included here:
# Restricted Hartree--Fock
include("RHF/RHF.jl")

# UHF
include("UHF/UHF.jl")

end #module

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
using Test
using SafeTestsets

@safetestset "SpinWaveTheory" begin
    include("SpinWaveTheory.jl")
end


text:
include("viscous_penalty.jl")

@doc """

Calculate fluxes at edge cubature points using face-based form.
eqn.flux_face, eqn.xflux, eqn.yflux will be updated.

Input:
  mesh :
  sbp  :
  eqn  :
  opts :
Output:

"""->
function calcViscousFlux_interior(mesh::AbstractDGMesh{Tmsh},
                                  sbp::AbstractOperator,
                                  eqn::NSData{Tsol, Tres, Tdim},
                                  opts) where {Tmsh, Tsol, Tres, Tdim}

  Ma      = eqn.params.euler_params.Ma
  Re      = eqn.params.Re
  gamma_1 = eqn.params.euler_params.gamma_1
  Pr      = 0.72
  coef_nondim = Ma/Re
  interfaces  = sview(mesh.interfaces, :)
  nfaces      = length(mesh.interfaces)
  p    = opts["order"]
  dq   = Array{Tsol}(mesh.numDofPerNode, mesh.numNodesPerFace)
  dqn  = Array{Tsol}(Tdim, mesh.numDofPerNode, mesh.numNodesPerFace)
  GtL  = zeros(Tsol, mesh.numDofPerNode, mesh.numDofPerNode, Tdim, Tdim, mesh.numNodesPerFace)
  GtR  = zeros(Tsol, mesh.numDofPerNode, mesh.numDofPerNode, Tdim, Tdim, mesh.numNodesPerFace)
  pMat  = zeros(Tsol, mesh.numDofPerNode, mesh.numDofPerNode, mesh.numNodesPerFace)
  Fv_faceL = zeros(Tsol, Tdim, mesh.numDofPerNode, mesh.numNodesPerFace)
  Fv_faceR = zeros(Tsol, Tdim, mesh.numDofPerNode, mesh.numNodesPerFace)
  Fv_avg   = zeros(Tsol, Tdim, mesh.numDofPerNode, mesh.numNodesPerFace)
  vecfluxL = zeros(Tsol, Tdim, mesh.numDofPerNode, mesh.numNodesPerFace)
  vecfluxR = zeros(Tsol, Tdim, mesh.numDofPerNode, mesh.numNodesPerFace)

  sbpface = mesh.sbpface
  sat_type = opts["SAT_type"]
  # this one is Harmann's definition
  const_tii = (p + 1.0)*(p + Tdim)/(2.0*Tdim)
  # const_tii = calcTraceInverseInequalityConst(sbp, sbpface)
  area_sum = sview(eqn.area_sum, :)

  params = eqn.params
  # sigma = calcTraceInverseInequalityConst(sbp, sbpface)
  # println("rho_max = ", sigma)

  for f = 1:nfaces    # loop over faces

    flux  = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerFace)
    face = interfaces[f]
    elemL = face.elementL
    elemR = face.elementR
    faceL = face.faceL
    faceR = face.faceR
    permL = sview(sbpface.perm, :, faceL)
    permR = sview(sbpface.perm, :, faceR)

    nrm_xy = ro_sview(mesh.nrm_face, :, :, f)

    # We need viscous flux and diffusion tensor on interfaces, and there
    # are different ways to compute them. For viscous flux:
    # 1) Fv = Fv(q, ∇q). Since we already have Q on face nodes, if 𝛻Q is also available
    # on interface, then we are done. This way is probably more consistent with computation
    # of other terms like Fv(q_b, ∇ϕ)
    # 2) we can compute viscous flux on volume nodes and then interpolate to interface nodes.
    # It's logically simple but computationally expensive.

    # compute viscous flux and diffusion tensor
    # Arrayviews returns a n x 1 x m array (ie. 3D, not 2D), so don't use it
    q_faceL = sview(eqn.q_face, :, 1, :, f)
    q_faceR = sview(eqn.q_face, :, 2, :, f)       # TODO: use sview instead of slice
    q_elemL = sview(eqn.q, :, :, elemL)
    q_elemR = sview(eqn.q, :, :, elemR)
    calcDiffusionTensor(eqn.params, q_faceL, GtL)
    calcDiffusionTensor(eqn.params, q_faceR, GtR)

    # one way to compute Fv_face
    # calcFvis_interiorFace(mesh, sbp, f, q_elemL, q_elemR, Fv_face)

    # compute the face derivatives first, i.e., we first compute
    # the derivatives at element nodes, and then do the interpolation.
    dqdx_face  = Array{Tsol}(Tdim, mesh.numDofPerNode, 2, mesh.numNodesPerFace)
    dqdx_elemL = Array{Tsol}(Tdim, mesh.numDofPerNode, mesh.numNodesPerElement)
    dqdx_elemR = Array{Tsol}(Tdim, mesh.numDofPerNode, mesh.numNodesPerElement)
    calcGradient(mesh, sbp, elemL, q_elemL, dqdx_elemL)
    calcGradient(mesh, sbp, elemR, q_elemR, dqdx_elemR)

    #
    # TODO: we can consider the first 2 dimension as a single dimension,
    # then we will not need slice here any more.
    #
    # for d = 1 : Tdim
      # dqdxL = sview(dqdx_elemL, d, :, :)
      # dqdxR = sview(dqdx_elemR, d, :, :)
      # dqdx_f = sview(dqdx_face, d, :, :, :)
      # interiorfaceinterpolate(sbpface, face, dqdxL, dqdxR, dqdx_f)
    # end

    # # Now both G and dqdx are avaiable at face nodes
    # dqdx_faceL = sview(dqdx_face, :, :, 1, :)
    # dqdx_faceR = sview(dqdx_face, :, :, 2, :)
    # calcFvis(params, GtL, dqdx_faceL, Fv_faceL)
    # calcFvis(params, GtR, dqdx_faceR, Fv_faceR)

    Fv_face = zeros(Tsol, Tdim, mesh.numDofPerNode, 2, mesh.numNodesPerFace)
    jacL = ro_sview(mesh.jac, :, elemL)
    jacR = ro_sview(mesh.jac, :, elemR)
    dxidxL = ro_sview(mesh.dxidx, :,:,:,elemL)
    dxidxR = ro_sview(mesh.dxidx, :,:,:,elemR)
    calcFaceFvis(params, sbp, sbpface, q_elemL, q_elemR, dxidxL, jacL, dxidxR, jacR, face, Fv_face)
    Fv_faceL = sview(Fv_face, :,:,1,:)
    Fv_faceR = sview(Fv_face, :,:,2,:)
    # diffL = maximum(abs.(real(view(Fv_face, :, :, 1, :) - Fv_faceL)))
    # diffR = maximum(abs.(real(view(Fv_face, :, :, 2, :) - Fv_faceR)))
    # if (diffL > 1.e-8)
        # println(diffL)
    # end

    cmptIPMat(mesh, sbp, eqn, opts, f, GtL, GtR, pMat)

    # Start to compute fluxes. We have 3 terms on interfaces:
    # 1) {Fv}⋅[ϕ]
    # 2) {G^T ∇ϕ}⋅[q]
    # 3) δ{G}[q]:[ϕ]

    # q jump
    for n = 1 : mesh.numNodesPerFace
      for iDof = 1 : mesh.numDofPerNode
        dq[iDof, n] = q_faceL[iDof, n] - q_faceR[iDof, n]
      end
    end

    # average viscous flux on face
    for n = 1 : mesh.numNodesPerFace
      for iDof = 2 : mesh.numDofPerNode
        for d = 1 : Tdim
          Fv_avg[d, iDof, n] = 0.5 * (Fv_faceL[d, iDof, n] + Fv_faceR[d, iDof, n] )
        end
      end
    end

    # Finally, everything is ready, let's compute fluxes, or penalties

    # This part computes the contribution of
    #     ∫ {G^T∇ϕ}:[q] dΓ = ∫ ∇ϕ⋅F dΓ ,
    # where
    #     [q] = (q+ - q-) ⊗ n = Δq⊗n ,
    # Then we can consider Δq⊗n as ∇q and F as viscous flux.
    fill!(vecfluxL, 0.0)
    fill!(vecfluxR, 0.0)
    for n = 1 : mesh.numNodesPerFace
      for iDof = 1 : mesh.numDofPerNode
        #
        # sum up columns of each row
        #
        for iDim = 1 : Tdim
          # vecfluxL[iDim, iDof, n] = 0.0
          # vecfluxR[iDim, iDof, n] = 0.0
          for jDim = 1 : Tdim
            tmpL = 0.0
            tmpR = 0.0
            for jDof = 1 : mesh.numDofPerNode
              tmpL += GtL[iDof, jDof, iDim, jDim, n]
              tmpR += GtR[iDof, jDof, iDim, jDim, n]
            end
            vecfluxL[iDim, iDof, n] += tmpL * nrm_xy[jDim, n]
            vecfluxR[iDim, iDof, n] += tmpR * nrm_xy[jDim, n]
          end
          vecfluxL[iDim,iDof,n] *=  dq[iDof,n]
          vecfluxR[iDim,iDof,n] *=  dq[iDof,n]
        end
      end
    end

    # δ{G}[q]:n, contributing to  δ{G}[q]:[ϕ]
    for n = 1:mesh.numNodesPerFace
      for iDof = 2 : mesh.numDofPerNode
        for jDof = 1 : mesh.numDofPerNode
          flux[iDof, n] +=  pMat[iDof, jDof, n]*dq[jDof, n]
        end
      end
    end

    # {Fv}⋅n, contributing to {Fv}⋅[ϕ]
    for n = 1:mesh.numNodesPerFace
      for iDof = 2 : mesh.numDofPerNode
        for iDim = 1 : Tdim
          flux[iDof, n] -= Fv_avg[iDim, iDof, n]*nrm_xy[iDim,n]
        end
      end
    end
    # accumulate fluxes
    for n = 1:mesh.numNodesPerFace
      for iDof = 2 : Tdim+2
        for iDim = 1 : Tdim
          eqn.vecflux_faceL[iDim, iDof, n, f] -=  vecfluxL[iDim, iDof, n]*coef_nondim
          eqn.vecflux_faceR[iDim, iDof, n, f] -=  vecfluxR[iDim, iDof, n]*coef_nondim
        end
        eqn.flux_face[iDof, n, f] += flux[iDof, n]*coef_nondim
      end
    end
  end # end of loop over all interfaces

  return nothing
end # end of function calcViscousFlux_interior

@doc """
Now actually we are integrating
∫ G∇ϕ:[q] dΓ
where G is the diffusion tensor and q is the solution variable. It is not
immediately in the 2nd form. A alternative (or better) way to do this
integral is as follows
∫ ∇ϕ⋅(G^t (qn))

Input:
Output:

"""->
function evalFaceIntegrals_vector(mesh::AbstractDGMesh{Tmsh},
                                  sbp::AbstractOperator,
                                  eqn::NSData{Tsol, Tres, Tdim},
                                  opts) where {Tmsh, Tsol, Tres, Tdim}
  # This part computes ∫ ∇ϕ⋅F  dΓ,
  sbpface = mesh.sbpface
  DxL = Array{Tmsh}(mesh.numNodesPerElement, mesh.numNodesPerElement, Tdim)
  DxR = Array{Tmsh}(mesh.numNodesPerElement, mesh.numNodesPerElement, Tdim)

  GtL = zeros(Tsol, mesh.numDofPerNode, mesh.numDofPerNode, Tdim, Tdim, mesh.numNodesPerFace)
  GtR = zeros(Tsol, mesh.numDofPerNode, mesh.numDofPerNode, Tdim, Tdim, mesh.numNodesPerFace)

  R = sview(sbpface.interp[:,:])
  w = sview(sbpface.wface, :)
  res = sview(eqn.res, :,:,:)

  numNodes_elem = mesh.numNodesPerElement    # number of Nodes per elemet
  numNodes_face = mesh.numNodesPerFace       # number of nodes on interfaces
  stencilsize = sbpface.stencilsize        # size of stencil for interpolation

  RDxL = Array{Tmsh}(mesh.numNodesPerFace, mesh.numNodesPerElement, Tdim)
  RDxR = Array{Tmsh}(mesh.numNodesPerFace, mesh.numNodesPerElement, Tdim)
  FvL  = zeros(Tsol, Tdim, mesh.numDofPerNode, mesh.numDofPerNode, mesh.numNodesPerFace, mesh.numNodesPerElement)
  FvR  = zeros(Tsol, Tdim, mesh.numDofPerNode, mesh.numDofPerNode, mesh.numNodesPerFace, mesh.numNodesPerElement)
  GtRDxL = zeros(Tsol, mesh.numDofPerNode, mesh.numDofPerNode, mesh.numNodesPerFace, mesh.numNodesPerElement)
  GtRDxR = zeros(Tsol, mesh.numDofPerNode, mesh.numDofPerNode, mesh.numNodesPerFace, mesh.numNodesPerElement)
  nrm    = Array{Tmsh}(Tdim, mesh.numNodesPerFace)
  dq     =  Array{Tsol}(mesh.numDofPerNode, mesh.numNodesPerFace)

  nfaces = length(mesh.interfaces)

  for f = 1 : nfaces
    face = mesh.interfaces[f]

    elemL = face.elementL
    elemR = face.elementR
    faceL = face.faceL
    faceR = face.faceR
    pL = sview(sbpface.perm, :, faceL)
    pR = sview(sbpface.perm, :, faceR)

    # compute RDx
    calcDx(mesh, sbp, elemL, DxL)
    calcDx(mesh, sbp, elemR, DxR)

    for i = 1 : length(RDxL)
      RDxL[i] = 0.0
      RDxR[i] = 0.0
    end

    for d =  1 : Tdim
      for row = 1 : numNodes_face
        rowR = sbpface.nbrperm[row, face.orient]
        for col = 1 : numNodes_elem
          for s = 1 : stencilsize
            RDxL[row, col, d] +=  R[s, row]  * DxL[pL[s], col, d]
            RDxR[row, col, d] +=  R[s, rowR] * DxR[pR[s], col, d]
          end
        end
      end
    end

    vecfluxL = sview(eqn.vecflux_faceL,:,:,:,f)
    vecfluxR = sview(eqn.vecflux_faceR,:,:,:,f)
    for i = 1 : numNodes_elem
      for j = 1 : numNodes_face
        for iDof = 2 : mesh.numDofPerNode
          tmpL = 0.0
          tmpR = 0.0
          for iDim = 1 : Tdim
            tmpL += RDxL[j, i, iDim] * vecfluxL[iDim, iDof, j]
            tmpR += RDxR[j, i, iDim] * vecfluxR[iDim, iDof, j]
          end
          res[iDof, i, elemL] += tmpL * w[j]
          res[iDof, i, elemR] += tmpR * w[j]
        end
      end
    end
  end

  return nothing
end

"""

Integrate ∫ ∇ϕ⋅F dΩ
TODO: consider combine it together with `weakdifferentiate`

Input:
  mesh   :
  sbp    :
  eqn    :
  res    :
Output:

"""
function weakdifferentiate2!(mesh::AbstractMesh{Tmsh},
                             sbp::AbstractOperator{Tsbp},
                             eqn::NSData{Tsol, Tres, Tdim},
                             res::AbstractArray{Tres,3}) where {Tmsh, Tsbp, Tsol, Tres, Tdim}
  @assert (sbp.numnodes ==  size(res,2))

  dim             = Tdim
  numElems        = mesh.numEl
  numNodesPerElem = mesh.numNodesPerElement
  numDofsPerNode  = mesh.numDofPerNode

  gamma_1 = eqn.params.euler_params.gamma_1
  Pr = 0.72
  Ma = eqn.params.euler_params.Ma
  Re = eqn.params.Re
  coef_nondim = Ma/Re

  Qx = Array{Tsbp}(numNodesPerElem, numNodesPerElem, dim)
  Fv = zeros(Tres, Tdim, numDofsPerNode, numNodesPerElem)
  w = sview(sbp.w, :)

  for elem = 1 : numElems
    # compute viscous flux
    q      = sview(eqn.q, :, :, elem)
    dxidx = sview(mesh.dxidx, :,:,:,elem)
    jac      = sview(mesh.jac, :, elem)

    calcFvis_elem(eqn.params, sbp, q, dxidx, jac, Fv)

    calcQx(mesh, sbp, elem, Qx)

    for d = 1 : dim
      for i = 1 : sbp.numnodes
        for j = 1 : sbp.numnodes
          for iDof = 2 : numDofsPerNode
            res[iDof, i, elem] -=  coef_nondim * Qx[j,i,d] * Fv[d, iDof, j] / jac[j]
          end
        end
      end
    end
  end
end

"""
  This flux function throws an error. Useful for defaults.
"""
mutable struct ErrorFlux <: FluxType
end

function (obj::ErrorFlux)(params::ParamType,
              uL::AbstractArray{Tsol,1},
              uR::AbstractArray{Tsol,1},
              aux_vars::AbstractVector{Tres},
              nrm::AbstractArray,
              F::AbstractArray{Tres}) where {Tsol, Tres}

  error("ErrorFlux called.")
  return nothing
end

"""
  Calls the [`SIPG`](@ref) (viscous) flux
"""
mutable struct SIPGViscousFlux <: FluxType
end

function (obj::SIPGViscousFlux)(params::ParamType,
              sbp::AbstractOperator,
              sbpface,    # TODO: type
              uL::AbstractArray{Tsol,1},
              uR::AbstractArray{Tsol,1},
              dxidxL,     # TODO: type
              jacL,       # TODO: type
              dxidxR,     # TODO: type
              jacR,       # TODO: type
              face,       # TODO: type
              F::AbstractVector{Tres}) where {Tsol, Tres}

  # calcViscousFlux_SIPG(params, uL, uR, aux_vars, nrm, F)    # this is the inviscid flux signature, needs to be changed
  calcViscousFlux_SIPG(params, sbp, sbpface, uL, uR, dxidxL, jacL, dxidxR, jacR, face, F)
  return nothing
end

"""
### NavierStokesMod.FluxDict

  This dictonary maps the names of the fluxes (Strings) to the
  functor object itself.  All flux functors should be added to the dictionary.

  TODO: document signature of the functors here

"""
global const FluxDict = Dict{String, FluxType}(
"SIPGViscousFlux" => SIPGViscousFlux(),
"ErrorFlux" => ErrorFlux(),
)

"""
### NavierStokesMod.getFluxFunctors

  This function retrieves the flux functors from the dictonary and
  stores them to eqn.viscous_flux_func.

  Inputs:
    mesh: an AbstractDGMesh
    sbp
    eqn
    opts
"""
function getFluxFunctors(mesh::AbstractDGMesh, sbp, eqn, opts)

  name = opts["Viscous_flux_name"]
  eqn.viscous_flux_func = FluxDict[name]

  return nothing
end


text:
module PortAudio

using alsa_plugins_jll: alsa_plugins_jll
import Base:
    close,
    eltype,
    getindex,
    getproperty,
    IteratorSize,
    isopen,
    iterate,
    length,
    read,
    read!,
    show,
    showerror,
    write
using Base.Iterators: flatten, repeated
using Base.Threads: @spawn
using libportaudio_jll: libportaudio
using LinearAlgebra: transpose!
import SampledSignals: nchannels, samplerate, unsafe_read!, unsafe_write
using SampledSignals: SampleSink, SampleSource
using Suppressor: @capture_err

export devices, PortAudioStream

include("libportaudio.jl")

using .LibPortAudio:
    paBadStreamPtr,
    Pa_CloseStream,
    PaDeviceIndex,
    PaDeviceInfo,
    PaError,
    PaErrorCode,
    Pa_GetDefaultInputDevice,
    Pa_GetDefaultOutputDevice,
    Pa_GetDeviceCount,
    Pa_GetDeviceInfo,
    Pa_GetErrorText,
    Pa_GetHostApiInfo,
    Pa_GetStreamReadAvailable,
    Pa_GetStreamWriteAvailable,
    Pa_GetVersion,
    Pa_GetVersionText,
    PaHostApiTypeId,
    Pa_Initialize,
    paInputOverflowed,
    Pa_IsStreamStopped,
    paNoFlag,
    Pa_OpenStream,
    paOutputUnderflowed,
    Pa_ReadStream,
    PaSampleFormat,
    Pa_StartStream,
    Pa_StopStream,
    PaStream,
    PaStreamParameters,
    Pa_Terminate,
    Pa_WriteStream

# for structs and strings
# PortAudio will return C_NULL instead of erroring
# so we need to handle these errors
function safe_load(result, an_error)
    if result == C_NULL
        throw(an_error)
    end
    unsafe_load(result)
end

# for functions that retrieve an index, throw a key error if it doesn't exist
function safe_key(a_function, an_index)
    safe_load(a_function(an_index), KeyError(an_index))
end

# error numbers are integers, while error codes are an @enum
function get_error_text(error_number)
    unsafe_string(Pa_GetErrorText(error_number))
end

struct PortAudioException <: Exception
    code::PaErrorCode
end
function showerror(io::IO, exception::PortAudioException)
    print(io, "PortAudioException: ")
    print(io, get_error_text(PaError(exception.code)))
end

# for integer results, PortAudio will return a negative number instead of erroring
# so we need to handle these errors
function handle_status(error_number; warn_xruns = true)
    if error_number < 0
        error_code = PaErrorCode(error_number)
        if error_code == paOutputUnderflowed || error_code == paInputOverflowed
            # warn instead of error after an xrun
            # allow users to disable these warnings
            if warn_xruns
                @warn("libportaudio: " * get_error_text(error_number))
            end
        else
            throw(PortAudioException(error_code))
        end
    end
    error_number
end

function write_debug(debug_file)
    debug_log = open(io -> read(io, String), debug_file, "r")
    if !isempty(debug_log)
        @debug debug_log
    end
end

function initialize()
    # ALSA will throw extraneous warnings on start-up
    # send them to debug instead
    mktemp() do debug_file, debug_io
        redirect_stderr(debug_io) do
            handle_status(Pa_Initialize())
        end
        write_debug(debug_file)
    end
end

function terminate()
    handle_status(Pa_Terminate())
end

# alsa needs to know where the configure file is
function seek_alsa_conf(folders)
    for folder in folders
        if isfile(joinpath(folder, "alsa.conf"))
            return folder
        end
    end
    throw(ArgumentError("Could not find alsa.conf in $folders"))
end

function __init__()
    if Sys.islinux()
        config_folder = "ALSA_CONFIG_DIR"
        if config_folder ∉ keys(ENV)
            ENV[config_folder] =
                seek_alsa_conf(("/usr/share/alsa", "/usr/local/share/alsa", "/etc/alsa"))
        end
        # the plugin folder will contain plugins for, critically, PulseAudio
        plugin_folder = "ALSA_PLUGIN_DIR"
        if plugin_folder ∉ keys(ENV) && alsa_plugins_jll.is_available()
            ENV[plugin_folder] = joinpath(alsa_plugins_jll.artifact_dir, "lib", "alsa-lib")
        end
    end
    initialize()
    atexit(() -> terminate())
end

function versioninfo(io::IO = stdout)
    println(io, unsafe_string(Pa_GetVersionText()))
    println(io, "Version: ", Pa_GetVersion())
end

# bounds for when a device is used as an input or output
struct Bounds
    max_channels::Int
    low_latency::Float64
    high_latency::Float64
end

struct PortAudioDevice
    name::String
    host_api::String
    default_sample_rate::Float64
    index::PaDeviceIndex
    input_bounds::Bounds
    output_bounds::Bounds
end

function PortAudioDevice(info::PaDeviceInfo, index)
    PortAudioDevice(
        unsafe_string(info.name),
        # replace host api code with its name
        unsafe_string(safe_key(Pa_GetHostApiInfo, info.hostApi).name),
        info.defaultSampleRate,
        index,
        Bounds(
            info.maxInputChannels,
            info.defaultLowInputLatency,
            info.defaultHighInputLatency,
        ),
        Bounds(
            info.maxOutputChannels,
            info.defaultLowOutputLatency,
            info.defaultHighInputLatency,
        ),
    )
end

name(device::PortAudioDevice) = device.name
function show(io::IO, device::PortAudioDevice)
    print(io, repr(name(device)))
    print(io, ' ')
    print(io, device.input_bounds.max_channels)
    print(io, '→')
    print(io, device.output_bounds.max_channels)
end

function get_default_input_index()
    handle_status(Pa_GetDefaultInputDevice())
end

function get_default_output_index()
    handle_status(Pa_GetDefaultOutputDevice())
end

# we can look up devices by index or name
function get_device(index::Integer)
    PortAudioDevice(safe_key(Pa_GetDeviceInfo, index), index)
end

function get_device(device_name::AbstractString)
    for device in devices()
        potential_match = name(device)
        if potential_match == device_name
            return device
        end
    end
    throw(KeyError(device_name))
end

"""
    devices()

List the devices available on your system.
Devices will be shown with their internal name, and maximum input and output channels.
"""
function devices()
    # need to use 0 indexing for C
    map(get_device, 0:(handle_status(Pa_GetDeviceCount()) - 1))
end

# we can handle reading and writing from buffers in a similar way
function read_or_write(a_function, buffer, use_frames = buffer.frames_per_buffer)
    handle_status(
        # because we're calling Pa_ReadStream and Pa_WriteStream from separate threads,
        # we put a lock around these calls
        lock(
            let a_function = a_function,
                pointer_to = buffer.pointer_to,
                data = buffer.data,
                use_frames = use_frames
                () -> a_function(pointer_to, data, use_frames)
            end,
            buffer.stream_lock,
        ),
        warn_xruns = buffer.warn_xruns,
    )
end

"""
    abstract type PortAudio.Scribe end

A scribe must implement the following:

  - A method for [`PortAudio.get_input_type`](@ref)
  - A method for [`PortAudio.get_output_type`](@ref)
  - A method to call itself on two arguments: a [`PortAudio.Buffer`](@ref) and an input of the input type.
    This method must return an output of the output type.
    This method should make use of [`PortAudio.read_buffer!`](@ref) and [`PortAudio.write_buffer`](@ref).
"""
abstract type Scribe end

abstract type SampledSignalsScribe <: Scribe end

"""
    struct PortAudio.SampledSignalsReader

A [`PortAudio.Scribe`](@ref) that will use the `SampledSignals` package to manage reading data from PortAudio.
"""
struct SampledSignalsReader <: SampledSignalsScribe end

"""
    struct PortAudio.SampledSignalsReader

A [`PortAudio.Scribe`](@ref) that will use the `SampledSignals` package to manage writing data to PortAudio.
"""
struct SampledSignalsWriter <: SampledSignalsScribe end

"""
    PortAudio.get_input_type(scribe::PortAudio.Scribe, Sample)

Get the input type of a [`PortAudio.Scribe`](@ref) for samples of type `Sample`.
"""
function get_input_type(::SampledSignalsScribe, Sample)
    # SampledSignals input_channel will be a triple of the last 3 arguments to unsafe_read/write
    # we will already have access to the stream itself
    Tuple{Array{Sample, 2}, Int, Int}
end

"""
    PortAudio.get_input_type(scribe::PortAudio.Scribe, Sample)

Get the output type of a [`PortAudio.Scribe`](@ref) for samples of type `Sample`.
"""
function get_output_type(::SampledSignalsScribe, Sample)
    # output is the number of frames read/written
    Int
end

# the julia buffer is bigger than the port audio buffer
# so we need to split it up into chunks
# we do this the same way for both reading and writing
function split_up(
    buffer,
    julia_buffer,
    already,
    frame_count,
    whole_function,
    partial_function,
)
    frames_per_buffer = buffer.frames_per_buffer
    # when we're done, we'll have written this many frames
    goal = already + frame_count
    # this is what we'll have left after doing all complete chunks
    left = frame_count % frames_per_buffer
    # this is how many we'll have written after doing all complete chunks
    even = goal - left
    foreach(
        let whole_function = whole_function,
            buffer = buffer,
            julia_buffer = julia_buffer,
            frames_per_buffer = frames_per_buffer
            already -> whole_function(
                buffer,
                frames_per_buffer,
                julia_buffer,
                (already + 1):(already + frames_per_buffer)
            )
        end,
        # start at the already, keep going until there is less than a chunk left
        already:frames_per_buffer:(even - frames_per_buffer),
        # each time we loop, add chunk frames to already
        # after the last loop, we'll reach "even"
    )
    # now we just have to read/write what's left
    if left > 0
        partial_function(buffer, left, julia_buffer, (even + 1):goal)
    end
    frame_count
end

# the full version doesn't have to make a view, but the partial version does
function full_write!(buffer, count, julia_buffer, julia_range)
    @inbounds transpose!(buffer.data, view(julia_buffer, julia_range, :))
    write_buffer(buffer, count)
end

function partial_write!(buffer, count, julia_buffer, julia_range)
    @inbounds transpose!(view(buffer.data, :, 1:count), view(julia_buffer, julia_range, :))
    write_buffer(buffer, count)
end

function (writer::SampledSignalsWriter)(buffer, arguments)
    split_up(buffer, arguments..., full_write!, partial_write!)
end

# similar to above
function full_read!(buffer, count, julia_buffer, julia_range)
    read_buffer!(buffer, count)
    @inbounds transpose!(view(julia_buffer, julia_range, :), buffer.data)
end

function partial_read!(buffer, count, julia_buffer, julia_range)
    read_buffer!(buffer, count)
    @inbounds transpose!(view(julia_buffer, julia_range, :), view(buffer.data, :, 1:count))
end

function (reader::SampledSignalsReader)(buffer, arguments)
    split_up(buffer, arguments..., full_read!, partial_read!)
end

"""
    struct PortAudio.Buffer{Sample}

A `PortAudio.Buffer` contains everything you might need to read or write data from or to PortAudio.
The `data` field contains the raw data in the buffer.
Use [`PortAudio.write_buffer`](@ref) to write data to PortAudio, and [`PortAudio.read_buffer!`](@ref) to read data from PortAudio.
"""
struct Buffer{Sample}
    stream_lock::ReentrantLock
    pointer_to::Ptr{PaStream}
    data::Array{Sample}
    number_of_channels::Int
    frames_per_buffer::Int
    warn_xruns::Bool
end

function Buffer(
    stream_lock,
    pointer_to,
    number_of_channels;
    Sample = Float32,
    frames_per_buffer = 128,
    warn_xruns = true,
)
    Buffer{Sample}(
        stream_lock,
        pointer_to,
        zeros(Sample, number_of_channels, frames_per_buffer),
        number_of_channels,
        frames_per_buffer,
        warn_xruns,
    )
end

eltype(::Type{Buffer{Sample}}) where {Sample} = Sample
nchannels(buffer::Buffer) = buffer.number_of_channels

"""
    PortAudio.write_buffer(buffer, use_frames = buffer.frames_per_buffer)

Write a number of frames (`use_frames`) from a [`PortAudio.Buffer`](@ref) to PortAudio.
"""
function write_buffer(buffer::Buffer, use_frames = buffer.frames_per_buffer)
    read_or_write(Pa_WriteStream, buffer, use_frames)
end

"""
    PortAudio.read_buffer!(buffer::Buffer, use_frames = buffer.frames_per_buffer)

Read a number of frames (`use_frames`) from PortAudio to a [`PortAudio.Buffer`](@ref).
"""
function read_buffer!(buffer, use_frames = buffer.frames_per_buffer)
    read_or_write(Pa_ReadStream, buffer, use_frames)
end

# the messanger will send tasks to the scribe
# the scribe will read/write from the buffer
struct Messanger{Sample, Scribe, Input, Output}
    device_name::String
    buffer::Buffer{Sample}
    scribe::Scribe
    input_channel::Channel{Input}
    output_channel::Channel{Output}
end

eltype(::Type{Messanger{Sample}}) where {Sample} = Sample
name(messanger::Messanger) = messanger.device_name
nchannels(messanger::Messanger) = nchannels(messanger.buffer)

# the scribe will be running on a separate thread in the background
# alternating transposing and
# waiting to pass inputs and outputs back and forth to PortAudio
function send(messanger)
    buffer = messanger.buffer
    scribe = messanger.scribe
    input_channel = messanger.input_channel
    output_channel = messanger.output_channel
    while true
        input = try
            take!(input_channel)
        catch an_error
            # if the input channel is closed, the scribe knows its done
            if an_error isa InvalidStateException && an_error.state === :closed
                break
            else
                rethrow(an_error)
            end
        end
        put!(output_channel, scribe(buffer, input))
    end
end

# convenience method
has_channels(something) = nchannels(something) > 0

# create the messanger, and start the scribe on a separate task
function messanger_task(
    device_name,
    buffer::Buffer{Sample},
    scribe::Scribe,
    debug_io,
) where {Sample, Scribe}
    Input = get_input_type(scribe, Sample)
    Output = get_output_type(scribe, Sample)
    input_channel = Channel{Input}(0)
    output_channel = Channel{Output}(0)
    # unbuffered channels so putting and taking will block till everyone's ready
    messanger = Messanger{Sample, Scribe, Input, Output}(
        device_name,
        buffer,
        scribe,
        input_channel,
        output_channel,
    )
    # we will spawn new threads to read from and write to port audio
    # while the reading thread is talking to PortAudio, the writing thread can be setting up, and vice versa
    # start the scribe thread when its created
    # if there's channels at all
    # we can't make the task a field of the buffer, because the task uses the buffer
    task = Task(let messanger = messanger, debug_io = debug_io
        # xruns will return an error code and send a duplicate warning to stderr
        # since we handle the error codes, we don't need the duplicate warnings
        # so we send them to a debug log
        () -> redirect_stderr(let messanger = messanger
            () -> send(messanger)
        end, debug_io)
    end)
    # makes it able to run on a separate thread
    task.sticky = false
    if has_channels(buffer)
        schedule(task)
        # output channel will close when the task ends
        bind(output_channel, task)
    else
        close(input_channel)
        close(output_channel)
    end
    messanger, task
end

function close_messanger_task(messanger, task)
    if has_channels(messanger)
        # this will shut down the channels, which will shut down the thread
        close(messanger.input_channel)
        # wait for tasks to finish to make sure any errors get caught
        wait(task)
        # output channel will close because it is bound to the task
    end
end

#
# PortAudioStream
#

struct PortAudioStream{SinkMessanger, SourceMessanger}
    sample_rate::Float64
    # pointer to the c object
    pointer_to::Ptr{PaStream}
    sink_messanger::SinkMessanger
    sink_task::Task
    source_messanger::SourceMessanger
    source_task::Task
    debug_file::String
    debug_io::IOStream
end

# portaudio uses codes instead of types for the sample format
const TYPE_TO_FORMAT = Dict{Type, PaSampleFormat}(
    Float32 => 1,
    Int32 => 2,
    # Int24 => 4,
    Int16 => 8,
    Int8 => 16,
    UInt8 => 3,
)

function make_parameters(
    device,
    channels,
    latency;
    Sample = Float32,
    host_api_specific_stream_info = C_NULL,
)
    if channels == 0
        # if we don't need any channels, we don't need the source/sink at all
        C_NULL
    else
        Ref(
            PaStreamParameters(
                device.index,
                channels,
                TYPE_TO_FORMAT[Sample],
                latency,
                host_api_specific_stream_info,
            ),
        )
    end
end

function fill_max_channels(kind, device, bounds, channels; adjust_channels = false)
    max_channels = bounds.max_channels
    if channels === maximum
        max_channels
    elseif channels > max_channels
        if adjust_channels
            max_channels
        else
            throw(
                DomainError(
                    channels,
                    "$channels exceeds maximum $kind channels for $(name(device))",
                ),
            )
        end
    else
        channels
    end
end

function combine_default_sample_rates(
    input_device,
    input_sample_rate,
    output_device,
    output_sample_rate,
)
    if input_sample_rate != output_sample_rate
        throw(
            ArgumentError(
                """
Default sample rate $input_sample_rate for input $(name(input_device)) disagrees with
default sample rate $output_sample_rate for output $(name(output_device)).
Please specify a sample rate.
""",
            ),
        )
    end
    input_sample_rate
end

# this is the top-level outer constructor that all the other outer constructors end up calling
"""
    PortAudioStream(input_channels = 2, output_channels = 2; options...)
    PortAudioStream(duplex_device, input_channels = 2, output_channels = 2; options...)
    PortAudioStream(input_device, output_device, input_channels = 2, output_channels = 2; options...)

Audio devices can either be `PortAudioDevice` instances as returned by [`devices`](@ref), or strings with the device name as reported by the operating system.
Set `input_channels` to `0` for an output only stream; set `output_channels` to `0` for an input only steam.
If you pass the function `maximum` instead of a number of channels, use the maximum channels allowed by the corresponding device.
If a single `duplex_device` is given, it will be used for both input and output.
If no devices are given, the system default devices will be used.

The `PortAudioStream` type supports all the stream and buffer features defined [SampledSignals.jl](https://github.com/JuliaAudio/SampledSignals.jl) by default.
For example, if you load SampledSignals with `using SampledSignals` you can read 5 seconds to a buffer with `buf = read(stream, 5s)`, regardless of the sample rate of the device.
`write(stream, stream)` will set up a loopback that will read from the input and play it back on the output.

Options:

  - `adjust_channels = false`: If set to `true`, if either `input_channels` or `output_channels` exceeds the corresponding device maximum, adjust down to the maximum.
  - `call_back = C_NULL`: The PortAudio call-back function.
    Currently, passing anything except `C_NULL` is unsupported.
  - `eltype = Float32`: Sample type of the audio stream
  - `flags = PortAudio.paNoFlag`: PortAudio flags
  - `frames_per_buffer = 128`: the number of frames per buffer
  - `input_info = C_NULL`: host API specific stream info for the input device.
    Currently, passing anything except `C_NULL` is unsupported.
  - `latency = nothing`: Requested latency. Stream could underrun when too low, consider using the defaults.  If left as `nothing`, use the defaults below:
    - For input/output only streams, use the default high latency.
    - For duplex streams, use the max of the default high latency for the input and output devices.
  - `output_info = C_NULL`: host API specific stream info for the output device.
    Currently, passing anything except `C_NULL` is unsupported.
  - `reader = PortAudio.SampledSignalsReader()`: the scribe that will read input.
    Defaults to a [`PortAudio.SampledSignalsReader`](@ref).
    Users can pass custom scribes; see [`PortAudio.Scribe`](@ref).
  - `samplerate = nothing`: Sample rate. If left as `nothing`, use the defaults below:
    - For input/output only streams, use the corresponding input/output device default sample rate.
    - For duplex streams, use the default sample rate if the default sample rates for the input and output devices match, otherwise throw an error.
  - `warn_xruns = true`: Display a warning if there is a stream overrun or underrun, which often happens when Julia is compiling, or with a particularly large GC run.
    Only affects duplex streams.
  - `writer = PortAudio.SampledSignalsWriter()`: the scribe that will write output.
    Defaults to a [`PortAudio.SampledSignalsReader`](@ref).
    Users can pass custom scribes; see [`PortAudio.Scribe`](@ref).

## Examples:

Set up an audio pass-through from microphone to speaker

```julia
julia> using PortAudio, SampledSignals

julia> stream = PortAudioStream(2, 2; warn_xruns = false);

julia> try
            # cancel with Ctrl-C
            write(stream, stream, 2s)
        finally
            close(stream)
        end
```

Use `do` syntax to auto-close the stream

```julia
julia> using PortAudio, SampledSignals

julia> PortAudioStream(2, 2; warn_xruns = false) do stream
            write(stream, stream, 2s)
        end
```

Open devices by name

```julia
using PortAudio, SampledSignals
PortAudioStream("Built-in Microph", "Built-in Output"; warn_xruns = false) do stream
    write(stream, stream, 2s)
end
2 s
```

Record 10 seconds of audio and save to an ogg file

```julia
julia> using PortAudio, SampledSignals, LibSndFile

julia> PortAudioStream(2, 0; warn_xruns = false) do stream
            buf = read(stream, 10s)
            save(joinpath(tempname(), ".ogg"), buf)
        end
2 s
```
"""
function PortAudioStream(
    input_device::PortAudioDevice,
    output_device::PortAudioDevice,
    input_channels = 2,
    output_channels = 2;
    eltype = Float32,
    adjust_channels = false,
    call_back = C_NULL,
    flags = paNoFlag,
    frames_per_buffer = 128,
    input_info = C_NULL,
    latency = nothing,
    output_info = C_NULL,
    reader = SampledSignalsReader(),
    samplerate = nothing,
    stream_lock = ReentrantLock(),
    user_data = C_NULL,
    warn_xruns = true,
    writer = SampledSignalsWriter(),
)
    debug_file, debug_io = mktemp()
    input_channels_filled = fill_max_channels(
        "input",
        input_device,
        input_device.input_bounds,
        input_channels;
        adjust_channels = adjust_channels,
    )
    output_channels_filled = fill_max_channels(
        "output",
        output_device,
        output_device.output_bounds,
        output_channels;
        adjust_channels = adjust_channels,
    )
    if input_channels_filled > 0
        if output_channels_filled > 0
            if latency === nothing
                latency = max(
                    input_device.input_bounds.high_latency,
                    output_device.output_bounds.high_latency,
                )
            end
            if samplerate === nothing
                samplerate = combine_default_sample_rates(
                    input_device,
                    input_device.default_sample_rate,
                    output_device,
                    output_device.default_sample_rate,
                )
            end
        else
            if latency === nothing
                latency = input_device.input_bounds.high_latency
            end
            if samplerate === nothing
                samplerate = input_device.default_sample_rate
            end
        end
    else
        if output_channels_filled > 0
            if latency === nothing
                latency = output_device.output_bounds.high_latency
            end
            if samplerate === nothing
                samplerate = output_device.default_sample_rate
            end
        else
            throw(ArgumentError("Input or output must have at least 1 channel"))
        end
    end
    # we need a mutable pointer so portaudio can set it for us
    mutable_pointer = Ref{Ptr{PaStream}}(0)
    handle_status(
        Pa_OpenStream(
            mutable_pointer,
            make_parameters(
                input_device,
                input_channels_filled,
                latency;
                host_api_specific_stream_info = input_info,
            ),
            make_parameters(
                output_device,
                output_channels_filled,
                latency;
                Sample = eltype,
                host_api_specific_stream_info = output_info,
            ),
            samplerate,
            frames_per_buffer,
            flags,
            call_back,
            user_data,
        ),
    )
    pointer_to = mutable_pointer[]
    handle_status(Pa_StartStream(pointer_to))
    PortAudioStream(
        samplerate,
        pointer_to,
        # we need to keep track of the tasks so we can wait for them to finish and catch errors
        messanger_task(
            output_device.name,
            Buffer(
                stream_lock,
                pointer_to,
                output_channels_filled;
                Sample = eltype,
                frames_per_buffer = frames_per_buffer,
                warn_xruns = warn_xruns,
            ),
            writer,
            debug_io,
        )...,
        messanger_task(
            input_device.name,
            Buffer(
                stream_lock,
                pointer_to,
                input_channels_filled;
                Sample = eltype,
                frames_per_buffer = frames_per_buffer,
                warn_xruns = warn_xruns,
            ),
            reader,
            debug_io,
        )...,
        debug_file,
        debug_io,
    )
end

# handle device names given as streams
function PortAudioStream(
    in_device_name::AbstractString,
    out_device_name::AbstractString,
    arguments...;
    keywords...,
)
    PortAudioStream(
        get_device(in_device_name),
        get_device(out_device_name),
        arguments...;
        keywords...,
    )
end

# if one device is given, use it for input and output
function PortAudioStream(
    device::Union{PortAudioDevice, AbstractString},
    input_channels = 2,
    output_channels = 2;
    keywords...,
)
    PortAudioStream(device, device, input_channels, output_channels; keywords...)
end

# use the default input and output devices
function PortAudioStream(input_channels = 2, output_channels = 2; keywords...)
    in_index = get_default_input_index()
    out_index = get_default_output_index()
    PortAudioStream(
        get_device(in_index),
        get_device(out_index),
        input_channels,
        output_channels;
        keywords...,
    )
end

# handle do-syntax
function PortAudioStream(do_function::Function, arguments...; keywords...)
    stream = PortAudioStream(arguments...; keywords...)
    try
        do_function(stream)
    finally
        close(stream)
    end
end

function close(stream::PortAudioStream)
    # closing is tricky, because we want to make sure we've read exactly as much as we've written
    # but we have don't know exactly what the tasks are doing
    # for now, just close one and then the other
    close_messanger_task(stream.source_messanger, stream.source_task)
    close_messanger_task(stream.sink_messanger, stream.sink_task)
    pointer_to = stream.pointer_to
    # only stop if it's not already stopped
    if !Bool(handle_status(Pa_IsStreamStopped(pointer_to)))
        handle_status(Pa_StopStream(pointer_to))
    end
    handle_status(Pa_CloseStream(pointer_to))
    # close the debug log and then read the file
    # this will contain duplicate xrun warnings mentioned above
    close(stream.debug_io)
    write_debug(stream.debug_file)
end

function isopen(pointer_to::Ptr{PaStream})
    # we aren't actually interested if the stream is stopped or not
    # instead, we are looking for the error which comes from checking on a closed stream
    error_number = Pa_IsStreamStopped(pointer_to)
    if error_number >= 0
        true
    else
        PaErrorCode(error_number) != paBadStreamPtr
    end
end
isopen(stream::PortAudioStream) = isopen(stream.pointer_to)

samplerate(stream::PortAudioStream) = stream.sample_rate
function eltype(
    ::Type{<:PortAudioStream{<:Messanger{Sample}, <:Messanger{Sample}}},
) where {Sample}
    Sample
end

# these defaults will error for non-sampledsignal scribes
# which is probably ok; we want these users to define new methods
read(stream::PortAudioStream, arguments...) = read(stream.source, arguments...)
read!(stream::PortAudioStream, arguments...) = read!(stream.source, arguments...)
write(stream::PortAudioStream, arguments...) = write(stream.sink, arguments...)
function write(sink::PortAudioStream, source::PortAudioStream, arguments...)
    write(sink.sink, source.source, arguments...)
end

function show(io::IO, stream::PortAudioStream)
    # just show the first type parameter (eltype)
    print(io, "PortAudioStream{")
    print(io, eltype(stream))
    println(io, "}")
    print(io, "  Samplerate: ", samplerate(stream), "Hz")
    # show source or sink if there's any channels
    sink = stream.sink
    if has_channels(sink)
        print(io, "\n  ")
        show(io, sink)
    end
    source = stream.source
    if has_channels(source)
        print(io, "\n  ")
        show(io, source)
    end
end

#
# PortAudioSink & PortAudioSource
#

# Define our source and sink types
# If we had multiple inheritance, then PortAudioStreams could be both a sink and source
# Since we don't, we have to make wrappers instead
for (TypeName, Super) in ((:PortAudioSink, :SampleSink), (:PortAudioSource, :SampleSource))
    @eval struct $TypeName{InputMessanger, OutputMessanger} <: $Super
        stream::PortAudioStream{InputMessanger, OutputMessanger}
    end
end

# provided for backwards compatibility
# only defined for SampledSignals scribes
function getproperty(
    stream::PortAudioStream{
        <:Messanger{<:Any, <:SampledSignalsWriter},
        <:Messanger{<:Any, <:SampledSignalsReader},
    },
    property::Symbol,
)
    if property === :sink
        PortAudioSink(stream)
    elseif property === :source
        PortAudioSource(stream)
    else
        getfield(stream, property)
    end
end

function nchannels(source_or_sink::PortAudioSource)
    nchannels(source_or_sink.stream.source_messanger)
end
function nchannels(source_or_sink::PortAudioSink)
    nchannels(source_or_sink.stream.sink_messanger)
end
function samplerate(source_or_sink::Union{PortAudioSink, PortAudioSource})
    samplerate(source_or_sink.stream)
end
function eltype(
    ::Type{
        <:Union{
            <:PortAudioSink{<:Messanger{Sample}, <:Messanger{Sample}},
            <:PortAudioSource{<:Messanger{Sample}, <:Messanger{Sample}},
        },
    },
) where {Sample}
    Sample
end
function isopen(source_or_sink::Union{PortAudioSink, PortAudioSource})
    isopen(source_or_sink.stream)
end
name(source_or_sink::PortAudioSink) = name(source_or_sink.stream.sink_messanger)
name(source_or_sink::PortAudioSource) = name(source_or_sink.stream.source_messanger)

# could show full type name, but the PortAudio part is probably redundant
# because these will usually only get printed as part of show for PortAudioStream
kind(::PortAudioSink) = "sink"
kind(::PortAudioSource) = "source"
function show(io::IO, sink_or_source::Union{PortAudioSink, PortAudioSource})
    print(
        io,
        nchannels(sink_or_source),
        " channel ",
        kind(sink_or_source),
        ": ",
        # put in quotes
        repr(name(sink_or_source)),
    )
end

# both reading and writing will outsource to the readers or writers
# so we just need to pass inputs in and take outputs out
# SampledSignals can take care of this feeding for us
function exchange(messanger, arguments...)
    put!(messanger.input_channel, arguments)
    take!(messanger.output_channel)
end

# these will only work with sampledsignals scribes
function unsafe_write(
    sink::PortAudioSink{<:Messanger{<:Any, <:SampledSignalsWriter}},
    julia_buffer::Array,
    already,
    frame_count,
)
    exchange(sink.stream.sink_messanger, julia_buffer, already, frame_count)
end

function unsafe_read!(
    source::PortAudioSource{<:Any, <:Messanger{<:Any, <:SampledSignalsReader}},
    julia_buffer::Array,
    already,
    frame_count,
)
    exchange(source.stream.source_messanger, julia_buffer, already, frame_count)
end

end # module PortAudio


text:
"""
    function elephant()

An elephant out and about

# Example
```jldoctest
julia> cowsay("Toot!", cow=Cowsay.elephant)
 _______
< Toot! >
 -------
 \\     /\\  ___  /\\
  \\   // \\/   \\/ \\\\
     ((    O O    ))
      \\\\ /     \\ //
       \\/  | |  \\/
        |  | |  |
        |  | |  |
        |   o   |
        | |   | |
        |m|   |m|
```
"""
function elephant(;eyes="oo", tongue="  ", thoughts="\\")
the_cow =
"""
 $thoughts     /\\  ___  /\\
  $thoughts   // \\/   \\/ \\\\
     ((    O O    ))
      \\\\ /     \\ //
       \\/  | |  \\/
        |  | |  |
        |  | |  |
        |   o   |
        | |   | |
        |m|   |m|
"""
return the_cow
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module TypesTest

using Test
using Dates

# Date internal algorithms
@testset "totaldays" begin
    @test DatesPlus.totaldays(0, 2, 28) == -307
    @test DatesPlus.totaldays(0, 2, 29) == -306
    @test DatesPlus.totaldays(0, 3, 1) == -305
    @test DatesPlus.totaldays(0, 12, 31) == 0
    # Rata Die Days # start from 0001-01-01
    @test DatesPlus.totaldays(1, 1, 1) == 1
    @test DatesPlus.totaldays(1, 1, 2) == 2
    @test DatesPlus.totaldays(2013, 1, 1) == 734869
end
@testset "daysinmonth" begin
    @test DatesPlus.daysinmonth(2000, 1) == 31
    @test DatesPlus.daysinmonth(2000, 2) == 29
    @test DatesPlus.daysinmonth(2000, 3) == 31
    @test DatesPlus.daysinmonth(2000, 4) == 30
    @test DatesPlus.daysinmonth(2000, 5) == 31
    @test DatesPlus.daysinmonth(2000, 6) == 30
    @test DatesPlus.daysinmonth(2000, 7) == 31
    @test DatesPlus.daysinmonth(2000, 8) == 31
    @test DatesPlus.daysinmonth(2000, 9) == 30
    @test DatesPlus.daysinmonth(2000, 10) == 31
    @test DatesPlus.daysinmonth(2000, 11) == 30
    @test DatesPlus.daysinmonth(2000, 12) == 31
    @test DatesPlus.daysinmonth(2001, 2) == 28
end
@testset "isleapyear" begin
    @test DatesPlus.isleapyear(1900) == false
    @test DatesPlus.isleapyear(2000) == true
    @test DatesPlus.isleapyear(2004) == true
    @test DatesPlus.isleapyear(2008) == true
    @test DatesPlus.isleapyear(0) == true
    @test DatesPlus.isleapyear(1) == false
    @test DatesPlus.isleapyear(-1) == false
    @test DatesPlus.isleapyear(4) == true
    @test DatesPlus.isleapyear(-4) == true
end
# Create "test" check manually
y = DatesPlus.Year(1)
m = DatesPlus.Month(1)
w = DatesPlus.Week(1)
d = DatesPlus.Day(1)
h = DatesPlus.Hour(1)
mi = DatesPlus.Minute(1)
s = DatesPlus.Second(1)
ms = DatesPlus.Millisecond(1)
@testset "DateTime construction by parts" begin
    test = DatesPlus.DateTime(DatesPlus.UTM(63492681600000))
    @test DatesPlus.DateTime(2013) == test
    @test DatesPlus.DateTime(2013, 1) == test
    @test DatesPlus.DateTime(2013, 1, 1) == test
    @test DatesPlus.DateTime(2013, 1, 1, 0) == test
    @test DatesPlus.DateTime(2013, 1, 1, 0, 0) == test
    @test DatesPlus.DateTime(2013, 1, 1, 0, 0, 0) == test
    @test DatesPlus.DateTime(2013, 1, 1, 0, 0, 0, 0) == test

    @test DatesPlus.DateTime(y) == DatesPlus.DateTime(1)
    @test DatesPlus.DateTime(y, m) == DatesPlus.DateTime(1, 1)
    @test DatesPlus.DateTime(y, m, d) == DatesPlus.DateTime(1, 1, 1)
    @test DatesPlus.DateTime(y, m, d, h) == DatesPlus.DateTime(1, 1, 1, 1)
    @test DatesPlus.DateTime(y, m, d, h, mi) == DatesPlus.DateTime(1, 1, 1, 1, 1)
    @test DatesPlus.DateTime(y, m, d, h, mi, s) == DatesPlus.DateTime(1, 1, 1, 1, 1, 1)
    @test DatesPlus.DateTime(y, m, d, h, mi, s, ms) == DatesPlus.DateTime(1, 1, 1, 1, 1, 1, 1)
    @test DatesPlus.DateTime(DatesPlus.Day(10), DatesPlus.Month(2), y) == DatesPlus.DateTime(1, 2, 10)
    @test DatesPlus.DateTime(DatesPlus.Second(10), DatesPlus.Month(2), y, DatesPlus.Hour(4)) == DatesPlus.DateTime(1, 2, 1, 4, 0, 10)
    @test DatesPlus.DateTime(DatesPlus.Year(1), DatesPlus.Month(2), DatesPlus.Day(1),
                         DatesPlus.Hour(4), DatesPlus.Second(10)) == DatesPlus.DateTime(1, 2, 1, 4, 0, 10)
end

@testset "Date construction by parts" begin
    test = DatesPlus.Date(DatesPlus.UTD(734869))
    @test DatesPlus.Date(2013) == test
    @test DatesPlus.Date(2013, 1) == test
    @test DatesPlus.Date(2013, 1, 1) == test
    @test DatesPlus.Date(y) == DatesPlus.Date(1)
    @test DatesPlus.Date(y, m) == DatesPlus.Date(1, 1)
    @test DatesPlus.Date(y, m, d) == DatesPlus.Date(1, 1, 1)
    @test DatesPlus.Date(m) == DatesPlus.Date(1, 1, 1)
    @test DatesPlus.Date(d, y) == DatesPlus.Date(1, 1, 1)
    @test DatesPlus.Date(d, m) == DatesPlus.Date(1, 1, 1)
    @test DatesPlus.Date(m, y) == DatesPlus.Date(1, 1, 1)
    @test DatesPlus.Date(DatesPlus.Day(10), DatesPlus.Month(2), y) == DatesPlus.Date(1, 2, 10)
end

@testset "Time construction by parts" begin
    t = DatesPlus.Time(DatesPlus.Nanosecond(82800000000000))
    @test DatesPlus.Time(23) == t
    @test DatesPlus.Time(23, 0) == t
    @test DatesPlus.Time(23, 0, 0) == t
    @test DatesPlus.Time(23, 0, 0, 0) == t
    @test DatesPlus.Time(23, 0, 0, 0, 0) == t
    @test DatesPlus.Time(23, 0, 0, 0, 0, 0) == t
    us = DatesPlus.Microsecond(1)
    ns = DatesPlus.Nanosecond(1)
    @test DatesPlus.Time(h) == DatesPlus.Time(1)
    @test DatesPlus.Time(h, mi) == DatesPlus.Time(1, 1)
    @test DatesPlus.Time(h, mi, s) == DatesPlus.Time(1, 1, 1)
    @test DatesPlus.Time(h, mi, s, ms) == DatesPlus.Time(1, 1, 1, 1)
    @test DatesPlus.Time(h, mi, s, ms, us) == DatesPlus.Time(1, 1, 1, 1, 1)
    @test DatesPlus.Time(h, mi, s, ms, us, ns) == DatesPlus.Time(1, 1, 1, 1, 1, 1)
    @test DatesPlus.Time(us, h, s, ns, mi, ms) == DatesPlus.Time(1, 1, 1, 1, 1, 1)
    @test DatesPlus.Time(DatesPlus.Second(10), DatesPlus.Minute(2), us, DatesPlus.Hour(4)) == DatesPlus.Time(4, 2, 10, 0, 1)
    @test DatesPlus.Time(DatesPlus.Hour(4), DatesPlus.Second(10), DatesPlus.Millisecond(15),
                     DatesPlus.Microsecond(20), DatesPlus.Nanosecond(25)) == DatesPlus.Time(4, 0, 10, 15, 20, 25)
end

@testset "various input types for Date/DateTime" begin
    test = DatesPlus.Date(1, 1, 1)
    @test DatesPlus.Date(Int8(1), Int8(1), Int8(1)) == test
    @test DatesPlus.Date(UInt8(1), UInt8(1), UInt8(1)) == test
    @test DatesPlus.Date(Int16(1), Int16(1), Int16(1)) == test
    @test DatesPlus.Date(UInt8(1), UInt8(1), UInt8(1)) == test
    @test DatesPlus.Date(Int32(1), Int32(1), Int32(1)) == test
    @test DatesPlus.Date(UInt32(1), UInt32(1), UInt32(1)) == test
    @test DatesPlus.Date(Int64(1), Int64(1), Int64(1)) == test
    @test DatesPlus.Date('\x01', '\x01', '\x01') == test
    @test DatesPlus.Date(true, true, true) == test
    @test_throws ArgumentError DatesPlus.Date(false, true, false)
    @test DatesPlus.Date(false, true, true) == test - DatesPlus.Year(1)
    @test_throws ArgumentError DatesPlus.Date(true, true, false)
    @test DatesPlus.Date(UInt64(1), UInt64(1), UInt64(1)) == test
    @test DatesPlus.Date(-1, UInt64(1), UInt64(1)) == test - DatesPlus.Year(2)
    @test DatesPlus.Date(Int128(1), Int128(1), Int128(1)) == test
    @test_throws InexactError DatesPlus.Date(170141183460469231731687303715884105727, Int128(1), Int128(1))
    @test DatesPlus.Date(UInt128(1), UInt128(1), UInt128(1)) == test
    @test DatesPlus.Date(big(1), big(1), big(1)) == test
    @test DatesPlus.Date(big(1), big(1), big(1)) == test
    # Potentially won't work if can't losslessly convert to Int64
    @test DatesPlus.Date(BigFloat(1), BigFloat(1), BigFloat(1)) == test
    @test DatesPlus.Date(complex(1), complex(1), complex(1)) == test
    @test DatesPlus.Date(Float64(1), Float64(1), Float64(1)) == test
    @test DatesPlus.Date(Float32(1), Float32(1), Float32(1)) == test
    @test DatesPlus.Date(Float16(1), Float16(1), Float16(1)) == test
    @test DatesPlus.Date(Rational(1), Rational(1), Rational(1)) == test
    @test_throws InexactError DatesPlus.Date(BigFloat(1.2), BigFloat(1), BigFloat(1))
    @test_throws InexactError DatesPlus.Date(1 + im, complex(1), complex(1))
    @test_throws InexactError DatesPlus.Date(1.2, 1.0, 1.0)
    @test_throws InexactError DatesPlus.Date(1.2f0, 1.f0, 1.f0)
    @test_throws InexactError DatesPlus.Date(3//4, Rational(1), Rational(1)) == test

    # Months, days, hours, minutes, seconds, and milliseconds must be in range
    @test_throws ArgumentError DatesPlus.Date(2013, 0, 1)
    @test_throws ArgumentError DatesPlus.Date(2013, 13, 1)
    @test_throws ArgumentError DatesPlus.Date(2013, 1, 0)
    @test_throws ArgumentError DatesPlus.Date(2013, 1, 32)
    @test_throws ArgumentError DatesPlus.DateTime(2013, 0, 1)
    @test_throws ArgumentError DatesPlus.DateTime(2013, 13, 1)
    @test_throws ArgumentError DatesPlus.DateTime(2013, 1, 0)
    @test_throws ArgumentError DatesPlus.DateTime(2013, 1, 32)
    @test_throws ArgumentError DatesPlus.DateTime(2013, 1, 1, 25)
    @test_throws ArgumentError DatesPlus.DateTime(2013, 1, 1, -1)
    @test_throws ArgumentError DatesPlus.DateTime(2013, 1, 1, 0, -1)
    @test_throws ArgumentError DatesPlus.DateTime(2013, 1, 1, 0, 60)
    @test_throws ArgumentError DatesPlus.DateTime(2013, 1, 1, 0, 0, -1)
    @test_throws ArgumentError DatesPlus.DateTime(2013, 1, 1, 0, 0, 60)
    @test_throws ArgumentError DatesPlus.DateTime(2013, 1, 1, 0, 0, 0, -1)
    @test_throws ArgumentError DatesPlus.DateTime(2013, 1, 1, 0, 0, 0, 1000)
    @test_throws ArgumentError DatesPlus.Time(24)
    @test_throws ArgumentError DatesPlus.Time(-1)
    @test_throws ArgumentError DatesPlus.Time(0, -1)
    @test_throws ArgumentError DatesPlus.Time(0, 60)
    @test_throws ArgumentError DatesPlus.Time(0, 0, -1)
    @test_throws ArgumentError DatesPlus.Time(0, 0, 60)
    @test_throws ArgumentError DatesPlus.Time(0, 0, 0, -1)
    @test_throws ArgumentError DatesPlus.Time(0, 0, 0, 1000)
    @test_throws ArgumentError DatesPlus.Time(0, 0, 0, 0, -1)
    @test_throws ArgumentError DatesPlus.Time(0, 0, 0, 0, 1000)
    @test_throws ArgumentError DatesPlus.Time(0, 0, 0, 0, 0, -1)
    @test_throws ArgumentError DatesPlus.Time(0, 0, 0, 0, 0, 1000)
end
a = DatesPlus.DateTime(2000)
b = DatesPlus.Date(2000)
c = DatesPlus.Time(0)
@testset "DateTime traits" begin
    @test DatesPlus.calendar(a) == DatesPlus.ISOCalendar
    @test DatesPlus.calendar(b) == DatesPlus.ISOCalendar
    @test eps(DateTime) == DatesPlus.Millisecond(1)
    @test eps(Date) == DatesPlus.Day(1)
    @test eps(Time) == DatesPlus.Nanosecond(1)
    @test eps(a) == DatesPlus.Millisecond(1)
    @test eps(b) == DatesPlus.Day(1)
    @test eps(c) == DatesPlus.Nanosecond(1)
    @test zero(DateTime) == DatesPlus.Millisecond(0)
    @test zero(Date) == DatesPlus.Day(0)
    @test zero(Time) == DatesPlus.Nanosecond(0)
    @test zero(a) == DatesPlus.Millisecond(0)
    @test zero(b) == DatesPlus.Day(0)
    @test zero(c) == DatesPlus.Nanosecond(0)
    @test string(typemax(DatesPlus.DateTime)) == "146138512-12-31T23:59:59"
    @test string(typemin(DatesPlus.DateTime)) == "-146138511-01-01T00:00:00"
    @test typemax(DatesPlus.DateTime) - typemin(DatesPlus.DateTime) == DatesPlus.Millisecond(9223372017043199000)
    @test string(typemax(DatesPlus.Date)) == "252522163911149-12-31"
    @test string(typemin(DatesPlus.Date)) == "-252522163911150-01-01"
    @test string(typemax(DatesPlus.Time)) == "23:59:59.999999999"
    @test string(typemin(DatesPlus.Time)) == "00:00:00"
    @test isfinite(DatesPlus.Date)
    @test isfinite(DatesPlus.DateTime)
    @test isfinite(DatesPlus.Time)
    @test c == c
    @test c == (c + DatesPlus.Hour(24))
    @test hash(c) == hash(c + DatesPlus.Hour(24))
    @test hash(c + DatesPlus.Nanosecond(10)) == hash(c + DatesPlus.Hour(24) + DatesPlus.Nanosecond(10))
end
@testset "Date-DateTime conversion/promotion" begin
    global a, b, c, d
    @test DatesPlus.DateTime(a) == a
    @test DatesPlus.Date(a) == b
    @test DatesPlus.DateTime(b) == a
    @test DatesPlus.Date(b) == b
    @test a == b
    @test a == a
    @test b == a
    @test b == b
    @test !(a < b)
    @test !(b < a)
    c = DatesPlus.DateTime(2000)
    d = DatesPlus.Date(2000)
    @test ==(a, c)
    @test ==(c, a)
    @test ==(d, b)
    @test ==(b, d)
    @test ==(a, d)
    @test ==(d, a)
    @test ==(b, c)
    @test ==(c, b)
    b = DatesPlus.Date(2001)
    @test b > a
    @test a < b
    @test a != b
    @test DatesPlus.Date(DatesPlus.DateTime(DatesPlus.Date(2012, 7, 1))) == DatesPlus.Date(2012, 7, 1)
end

@testset "min and max" begin
    for (a, b) in [(DatesPlus.Date(2000), DatesPlus.Date(2001)),
                    (DatesPlus.Time(10), DatesPlus.Time(11)),
                    (DatesPlus.DateTime(3000), DatesPlus.DateTime(3001)),
                    (DatesPlus.Week(42), DatesPlus.Week(1972)),
                    (DatesPlus.Quarter(3), DatesPlus.Quarter(52))]
        @test min(a, b) == a
        @test min(b, a) == a
        @test min(a) == a
        @test max(a, b) == b
        @test max(b, a) == b
        @test max(b) == b
        @test minmax(a, b) == (a, b)
        @test minmax(b, a) == (a, b)
        @test minmax(a) == (a, a)
    end
end

@testset "issue #31524" begin
    dt1 = Libc.strptime("%Y-%M-%dT%H:%M:%SZ", "2018-11-16T10:26:14Z")
    dt2 = Libc.TmStruct(14, 30, 5, 10, 1, 99, 3, 40, 0)

    time = Time(dt1)
    @test typeof(time) == Time
    @test time == DatesPlus.Time(10, 26, 14)

    date = Date(dt2)
    @test typeof(date) == Date
    @test date == DatesPlus.Date(1999, 2, 10)

    datetime = DateTime(dt2)
    @test typeof(datetime) == DateTime
    @test datetime == DatesPlus.DateTime(1999, 2, 10, 5, 30, 14)

end

@testset "timedwait" begin
    @test timedwait(() -> false, Second(0); pollint=Millisecond(1)) === :timed_out
end

end


text:
using ProjSplx
using Base.Test

srand(0)

# KKT residuals for projection onto simplex.
function kkt_splx_pass(x, b; c = ones(b), τ=1.0)
    ϵ = eps(1.0)^(2/3)*maximum(abs,b)
    r = x - b
    γ = dot(x,r)/τ
    z = r - γ*c
    t1 = τ - ϵ ≤ dot(c,x) ≤ τ + ϵ
    t2 = minimum(x) ≥ 0
    t3 = abs(dot(x, z)) ≤ ϵ
    return t1 & t2 & t3
end

# Opt gap for projection onto 1-norm ball.
function kkt_nrm1_pass(x, b; c = ones(b), τ=1.0)
    ϵ = eps(1.0)^(2/3)*maximum(abs, b)
    r = b - x
    prObj = 0.5*dot(r,r)
    duObj = dot(b,r) - 0.5*dot(r,r) - τ*norm(r./c,Inf)
    fea = norm(c.*x,1) ≤ τ+ϵ
    gap = (prObj - duObj) ≤ ϵ
    return gap & fea
end

@testset "Projection tests" begin

    @testset "Unit simplex" begin
        n = 110
        τ = rand()
        b = randn(n)

        # in place version
        x = copy(b)
        projsplx!(x, τ)
        @test kkt_splx_pass(x, b, τ=τ) == true

        # return version
        x = projsplx(b, τ)
        @test kkt_splx_pass(x, b, τ=τ) == true
    end

    @testset "Weighted unit simplex" begin
        n = 101
        τ = rand()
        b = randn(n)
        c = rand(n)

        # in place version
        x = copy(b)
        projsplx!(x, c, τ)
        @test kkt_splx_pass(x, b, c=c, τ=τ) == true

        # return version
        x = projsplx(b, c, τ)
        @test kkt_splx_pass(x, b, c=c, τ=τ) == true
    end

    n = 99
    b = randn(n)
    c = rand(n)

    @testset "1-norm ball (from outside)" begin
        τ = 0.1*norm(b,1)

        # in-place version
        x = copy(b)
        projnorm1!(x, τ)
        @test kkt_nrm1_pass(x, b, τ=τ) == true

        # return version
        x = projnorm1(b, τ)
        @test kkt_nrm1_pass(x, b, τ=τ) == true
    end

    @testset "1-norm ball (from inside)" begin
        τ = 1.1*norm(b,1)

        # in-place version
        x = copy(b)
        projnorm1!(x, τ)
        @test kkt_nrm1_pass(x, b, τ=τ) == true

        # return version
        x = projnorm1(b, τ)
        @test kkt_nrm1_pass(x, b, τ=τ) == true
    end

    @testset "weighted 1-norm ball" begin
        τ = 0.1*norm(b,1)

        # in-place version
        x = copy(b)
        projnorm1!(x, c, τ)
        @test kkt_nrm1_pass(x, b, c=c, τ=τ) == true

        # return version
        x = projnorm1(b, c, τ)
        @test kkt_nrm1_pass(x, b, c=c, τ=τ) == true

    end

end