text:
#=
        ScalarQuadraticFunction -in- LessThan
        ScalarQuadraticFunction -in- GreaterThan
        ScalarQuadraticFunction -in- EqualTo

    TODO:
        - constraint sets
        - constraint functions
        - deleting constraints
=#
constrdict(model::LinQuadOptimizer, ::QCI{LE})  = cmap(model).q_less_than
constrdict(model::LinQuadOptimizer, ::QCI{GE})  = cmap(model).q_greater_than
constrdict(model::LinQuadOptimizer, ::QCI{EQ})  = cmap(model).q_equal_to

"""
    reduce_duplicates!(rows::Vector{T}, cols::Vector{T}, vals::Vector{S})

Given a matrix specified by row indices in `rows`, column indices in `cols` and
coefficients in `vals`, return new `rows`, `cols`, and `vals` vectors with
duplicate elements summed and any coefficients in the lower triangle moved to
the upper triangle.

This function swaps element `i` in `rows` and `cols` if `rows[i]>cols[i]`.

# Examples
```jldoctest
julia> reduce_duplicates!(
    [1,   2, 2, 2],  # rows
    [1,   1, 2, 2],  # cols
    [1, 0.5, 1, 1]   # vals
    )
([1, 1, 2], [1, 2, 2], [1.0, 0.5, 2.0])
```
"""
function reduce_duplicates!(rows::Vector{T}, cols::Vector{T}, vals::Vector{S}) where T where S
    @assert length(rows) == length(cols) == length(vals)
    for i in 1:length(rows)
        if rows[i] > cols[i]
            tmp = rows[i]
            rows[i] = cols[i]
            cols[i] = tmp
        end
    end
    findnz(sparse(rows, cols, vals))
end

"""
    canonical_reduction(model::LinQuadOptimizer, func::Quad)

Reduce a ScalarQuadraticFunction into five arrays, returned in the following
order:
 1. a vector of affine column indices
 2. a vector of affine coefficients
 3. a vector of quadratic row indices
 4. a vector of quadratic column indices
 5. a vector of quadratic coefficients
"""
function canonical_reduction(model::LinQuadOptimizer, func::Quad)
    quad_columns_1, quad_columns_2, quad_coefficients = reduce_duplicates!(
        [get_column(model, term.variable_index_1) for term in func.quadratic_terms],
        [get_column(model, term.variable_index_2) for term in func.quadratic_terms],
        [term.coefficient for term in func.quadratic_terms]
    )
    affine_columns = [get_column(model, term.variable_index) for term in func.affine_terms]
    affine_coefficients = [term.coefficient for term in func.affine_terms]
    return affine_columns, affine_coefficients, quad_columns_1, quad_columns_2, quad_coefficients
end

function MOI.addconstraint!(model::LinQuadOptimizer, func::Quad, set::S) where S <: Union{LE, GE, EQ}
    add_quadratic_constraint(model, func, set)
    model.last_constraint_reference += 1
    index = QCI{S}(model.last_constraint_reference)
    dict = constrdict(model, index)
    push!(model.qconstraint_primal_solution, NaN)
    push!(model.qconstraint_dual_solution, NaN)
    dict[index] = get_number_quadratic_constraints(model)
    return index
end

function add_quadratic_constraint(model::LinQuadOptimizer, func::Quad, set::S) where S<: Union{LE, GE, EQ}
    add_quadratic_constraint(model, func, backend_type(model, set),
                             MOIU.getconstant(set))
end

function add_quadratic_constraint(model::LinQuadOptimizer, func::Quad, sense, rhs::Float64)
    if abs(func.constant) > 0
        Compat.@warn("Constant in quadratic function. Moving into set")
    end
    (aff_cols, aff_coeffs, I, J, V) = canonical_reduction(model, func)
    add_quadratic_constraint!(model, aff_cols, aff_coeffs, rhs - func.constant,
                              sense, I, J, V)
end

function MOI.delete!(model::LinQuadOptimizer, index::QCI{<: LinSets})
    __assert_valid__(model, index)
    delete_constraint_name(model, index)
    dict = constrdict(model, index)
    row = dict[index]
    delete_quadratic_constraints!(model, row, row)
    deleteat!(model.qconstraint_primal_solution, row)
    deleteat!(model.qconstraint_dual_solution, row)
    # shift all the other references
    shift_references_after_delete_quadratic!(model, row)
    delete!(dict, index)
end

#=
    Constraint set of Linear function
=#

MOI.canget(::LinQuadOptimizer, ::MOI.ConstraintSet, ::Type{QCI{S}}) where S <: Union{LE, GE, EQ} = true
function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintSet, index::QCI{S}) where S <: Union{LE, GE, EQ}
    rhs = get_quadratic_rhs(model, model[index])
    S(rhs)
end
MOI.canget(::LinQuadOptimizer, ::MOI.ConstraintSet, ::Type{QCI{IV}}) = false

#=
    Constraint function of Linear function
=#

MOI.canget(::LinQuadOptimizer, ::MOI.ConstraintFunction, ::Type{<:QCI{<: LinSets}}) = true
function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintFunction, index::QCI{<: LinSets})
    aff_cols, aff_coeffs, Q = get_quadratic_constraint(model, model[index])
    affine_terms = map(
        (col, coef)->MOI.ScalarAffineTerm{Float64}(coef, model.variable_references[col]),
        aff_cols, aff_coeffs)
    rows = rowvals(Q)
    vals = nonzeros(Q)
    nrows, ncols = size(Q)
    quadratic_terms = MOI.ScalarQuadraticTerm{Float64}[]
    sizehint!(quadratic_terms, length(vals))
    for i = 1:ncols
        for j in nzrange(Q, i)
            row = rows[j]
            val = vals[j]
            push!(quadratic_terms, MOI.ScalarQuadraticTerm{Float64}(2*val,
                model.variable_references[row], model.variable_references[i]))
        end
    end
    Quad(affine_terms, quadratic_terms, 0.0) # constant was moved into set
end


text:
# Test linearity of sources
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#

using PyCall, PyPlot, JUDI.TimeModeling, Test, LinearAlgebra

## Set up model structure
n = (120,100)	# (x,y,z) or (x,z)
d = (10.,10.)
o = (0.,0.)

# Velocity [km/s]
v = ones(Float32,n) .+ 0.4f0
v[:,Int(round(end/2)):end] .= 4f0

# Slowness squared [s^2/km^2]
m = (1f0 ./ v).^2

# Setup info and model structure
nsrc = 1
ntComp = 250
info = Info(prod(n), nsrc, ntComp)	# number of gridpoints, number of experiments, number of computational time steps
model = Model(n,d,o,m)

## Set up receiver geometry
nxrec = 120
xrec = range(50f0,stop=1150f0,length=nxrec)
yrec = 0f0
zrec = range(50f0,stop=50f0,length=nxrec)

# receiver sampling and recording time
timeR = 1000f0	# receiver recording time [ms]
dtR = 4f0	# receiver sampling interval

# Set up receiver structure
recGeometry = Geometry(xrec,yrec,zrec;dt=dtR,t=timeR,nsrc=nsrc)

## Set up source geometries (cell array with source locations for each shot)
xsrc1 = 300f0
ysrc1 = 0f0
zsrc1 = 20f0

xsrc2 = 600f0
ysrc2 = 0f0
zsrc2 = 20f0

# source sampling and number of time steps
timeS = 1000f0
dtS = 4f0	# receiver sampling interval

# Set up source structure
srcGeometry1 = Geometry(xsrc1,ysrc1,zsrc1;dt=dtS,t=timeS)
srcGeometry2 = Geometry(xsrc2,ysrc2,zsrc2;dt=dtS,t=timeS)

# setup wavelet
f0 = 0.01f0
wavelet = ricker_wavelet(timeS,dtS,f0)

###################################################################################################

# Modeling operators
Pr = judiProjection(info,recGeometry)
Ps1 = judiProjection(info,srcGeometry1)
Ps2 = judiProjection(info,srcGeometry2)
F = judiModeling(info,model)
q1 = judiVector(srcGeometry1,wavelet)
q2 = judiVector(srcGeometry2,wavelet)

d1 = Pr*F*adjoint(Ps1)*q1
d2 = Pr*F*adjoint(Ps2)*q2
d3 = Pr*F*(adjoint(Ps1)*q1 + adjoint(Ps2)*q2)
d4 = Pr*F*(adjoint(Ps1)*q1 - adjoint(Ps2)*q2)

@test isapprox(norm(d3), norm(d1 + d2))
@test isapprox(norm(d4), norm(d1 - d2))

ext: jl


text:
@testset "Fitting" begin
  wl = geostatsimage("WalkerLake")
  TI = asarray(wl, :Z)[1:20,1:20]
  d = georef((z=TI,))
  g = EmpiricalVariogram(d, :z, maxlag=15.)

  # all fits lead to similar sill
  γ₁ = fit(GaussianVariogram, g)
  γ₂ = fit(SphericalVariogram, g)
  γ₃ = fit(ExponentialVariogram, g)
  γ₄ = fit(MaternVariogram, g)
  @test isapprox(sill(γ₁), 0.054, atol=1e-3)
  @test isapprox(sill(γ₂), 0.054, atol=1e-3)
  @test isapprox(sill(γ₃), 0.058, atol=1e-3)
  @test isapprox(sill(γ₄), 0.057, atol=1e-3)

  # best fit is a Gaussian variogram
  γ = fit(Variogram, g)
  @test γ isa GaussianVariogram
  @test isapprox(sill(γ), 0.054, atol=1e-3)

  if visualtests
    plts = map([γ₁,γ₂,γ₃,γ₄]) do γ
      plot(g, legend=false)
      plot!(γ, 0., 15.)
    end
    @test_reference "data/fitting.png" plot(plts...)
  end

  # make sure convenient methods work
  γ₁ = fit(GaussianVariogram, g, h -> 1/h)
  γ₂ = fit(Variogram, g, h -> 1/h)
  @test sill(γ₁) > 0
  @test sill(γ₂) > 0

  # unitful types
  wl = geostatsimage("WalkerLake")
  TI = asarray(wl, :Z)[1:20,1:20]
  d = georef((z=TI*u"K",))
  g = EmpiricalVariogram(d, :z, maxlag=15.)
  γ = fit(Variogram, g)
  @test unit(sill(γ)) == u"K^2"
  @test unit(nugget(γ)) == u"K^2"
end


text:
include("iterators.jl")
include("helper_functions.jl")
include("decays.jl")
include("callbacks.jl")


text:
using SparseArrays, Base64
import BlockArrays: _BlockArray

@testset "block constructors" begin
    ret = BlockArray{Float64}(undef, 1:3)
    fill!(ret, 0)
    @test Array(ret)  == zeros(6)

    ret = BlockArray{Float64,1}(undef, 1:3)
    fill!(ret, 0)
    @test Array(ret)  == zeros(6)

    ret = BlockArray{Float64,1,Vector{Float64}}(undef, 1:3)
    fill!(ret, 0)
    @test Array(ret)  == zeros(6)

    ret = BlockArray{Float64}(undef, BlockArrays.BlockSizes(1:3))
    fill!(ret, 0)
    @test Array(ret)  == zeros(6)

    ret = BlockArray{Float64,1}(undef, BlockArrays.BlockSizes(1:3))
    fill!(ret, 0)
    @test Array(ret)  == zeros(6)

    ret = BlockArray{Float64,1,Vector{Float64}}(undef, BlockArrays.BlockSizes(1:3))
    fill!(ret, 0)
    @test Array(ret)  == zeros(6)

    ret = BlockArrays._BlockArray([[0.0],[0.0,0.0],[0.0,0.0,0.0]], 1:3)
    @test Array(ret)  == zeros(6)

    ret = BlockArrays._BlockArray([[0.0],[0.0,0.0],[0.0,0.0,0.0]], BlockArrays.BlockSizes(1:3))
    @test Array(ret)  == zeros(6)

    ret = BlockArray{Float32}(undef_blocks, 1:3)
    @test eltype(ret.blocks) == Vector{Float32}
    @test_throws UndefRefError ret.blocks[1]

    ret = BlockArray{Float32,1}(undef_blocks, 1:3)
    @test eltype(ret.blocks) == Vector{Float32}
    @test_throws UndefRefError ret.blocks[1]

    ret = BlockArray{Float32,1,Vector{Float32}}(undef_blocks, 1:3)
    @test eltype(ret.blocks) == Vector{Float32}
    @test_throws UndefRefError ret.blocks[1]

    ret = BlockArray{Float32}(undef_blocks, 1:3, 1:3)
    @test eltype(ret.blocks) == Matrix{Float32}
    @test_throws UndefRefError ret.blocks[1]

    ret = BlockArray(undef_blocks, Vector{Float32}, 1:3)
    @test eltype(ret) == Float32
    @test eltype(ret.blocks) == Vector{Float32}
    @test_throws UndefRefError ret.blocks[1]

    ret = BlockArray{Float64}(undef, 1:3, 1:3)
    fill!(ret, 0)
    Matrix(ret) == zeros(6,6)

    ret = PseudoBlockArray{Float64}(undef, 1:3)
    fill!(ret, 0)
    @test Array(ret)  == zeros(6)

    ret = PseudoBlockArray{Float64,1}(undef, 1:3)
    fill!(ret, 0)
    @test Array(ret)  == zeros(6)

    ret = PseudoBlockArray{Float64,1,Vector{Float64}}(undef, 1:3)
    fill!(ret, 0)
    @test Array(ret)  == zeros(6)

    ret = PseudoBlockArray{Float64}(undef, 1:3, 1:3)
    fill!(ret, 0)
    Matrix(ret) == zeros(6,6)

    @test_throws DimensionMismatch BlockArray([1,2,3],[1,1])
end

@testset "block indexing" begin
    BA_1 = BlockArray(undef_blocks, Vector{Float64}, [1,2,3])
    a_1 = rand(2)
    BA_1[Block(2)] = a_1
    @test BA_1[BlockIndex(2, 1)] == a_1[1]
    @test BA_1[BlockIndex(2, 2)] == a_1[2]
    @test BA_1[Block(2)] == a_1
    @test BA_1[2] == a_1[1]

    @test_throws DimensionMismatch (BA_1[Block(3)] = rand(4))
    @test_throws BlockBoundsError blockcheckbounds(BA_1, 4)
    @test_throws BlockBoundsError BA_1[Block(4)]

    BA_2 = BlockArray(undef_blocks, Matrix{Float64}, [1,2], [3,4])
    a_2 = rand(1,4)
    BA_2[Block(1,2)] = a_2
    @test BA_2[Block(1,2)] == a_2
    BA_2[Block(1,2)] = a_2

    @test BA_2[1,5] == a_2[2]
    @test_throws DimensionMismatch BA_2[Block(1,2)] = rand(1,5)
end

@testset "misc block tests" begin
    for BlockType in (BlockArray, PseudoBlockArray)
        a_1 = rand(6)
        BA_1 = BlockType(a_1, [1,2,3])
        @test Array(BA_1) == a_1
        @test nblocks(BA_1) == (3,)
        @test nblocks(BA_1,1) == 3
        @test eltype(similar(BA_1, Float32)) == Float32
        q = rand(1)
        BA_1[Block(1)] = q
        BA_1[BlockIndex(3, 2)] = a_1[5]
        @test BA_1[Block(1)] == q
        if BlockType == PseudoBlockArray
            q2 = zero(q)
            getblock!(q2, BA_1, 1)
            @test q2 == q
            @test_throws DimensionMismatch getblock!(zeros(2), BA_1, 1)
            fill!(q2, 0)
            getblock!(q2, BA_1, 1)
            @test q2 == q
        end
        fill!(BA_1, 1.0)
        @test BA_1 == ones(size(BA_1))
        ran = rand(size(BA_1)...)
        copyto!(BA_1, ran)
        @test BA_1 == ran

        a_1_sparse = sprand(6, 0.9)
        BA_1_sparse = BlockType(a_1_sparse, [1,2,3])
        @test Array(BA_1_sparse) == a_1_sparse
        BA_1_sparse[4] = 3.0
        @test BA_1_sparse[4] == 3.0

        a_2 = rand(3, 7)
        BA_2 = BlockType(a_2, [1,2], [3,4])
        @test Array(BA_2) == a_2
        @test nblocks(BA_2) == (2,2)
        @test nblocks(BA_2, 1) == 2
        @test nblocks(BA_2, 2, 1) == (2, 2)
        BA_2[BlockIndex((2,1), (2,2))] = a_2[3,2]
        @test eltype(similar(BA_2, Float32)) == Float32
        q = rand(1,4)
        BA_2[Block(1,2)] = q
        @test_throws DimensionMismatch BA_2[Block(1,2)] = rand(1,5)
        @test BA_2[Block(1,2)] == q
        if BlockType == PseudoBlockArray
            q2 = zero(q)
            getblock!(q2, BA_2, 1, 2)
            @test q2 == q
            @test_throws DimensionMismatch getblock!(zeros(1,5), BA_2, 1, 2)
        end
        fill!(BA_2, 1.0)
        @test BA_2 == ones(size(BA_2))
        ran = rand(size(BA_2)...)
        copyto!(BA_2, ran)
        @test BA_2 == ran

        a_2_sparse = sprand(3, 7, 0.9)
        BA_2_sparse = BlockType(a_2_sparse, [1,2], [3,4])
        @test Array(BA_2_sparse) == a_2_sparse
        BA_2_sparse[1,2] = 3.0
        @test BA_2_sparse[1,2] == 3.0

        a_3 = rand(3, 7,4)
        BA_3 = BlockType(a_3, [1,2], [3,4], [1,2,1])
        @test Array(BA_3) == a_3
        @test nblocks(BA_3) == (2,2,3)
        @test nblocks(BA_3, 1) == 2
        @test nblocks(BA_3, 3, 1) == (3, 2)
        @test nblocks(BA_3, 3) == 3
        BA_3[BlockIndex((1,1,1), (1,1,1))] = a_3[1,1,1]
        @test eltype(similar(BA_3, Float32)) == Float32
        q = rand(1,4,2)
        BA_3[Block(1,2,2)] = q
        @test BA_3[Block(1,2,2)] == q
        if BlockType == PseudoBlockArray
            q3 = zero(q)
            getblock!(q3, BA_3, 1, 2, 2)
            @test q3 == q
            @test_throws DimensionMismatch getblock!(zeros(1,3,2), BA_3, 1, 2,2)
        end
        fill!(BA_3, 1.0)
        @test BA_3 == ones(size(BA_3))
        ran = rand(size(BA_3)...)
        copyto!(BA_3, ran)
        @test BA_3 == ran
    end
end

@testset "convert" begin
    # Could probably be DRY'd.
    A = PseudoBlockArray(rand(2,3), [1,1], [2,1])
    C = convert(BlockArray, A)
    @test C == A == BlockArray(A)
    @test eltype(C) == eltype(A)

    C = convert(BlockArray{Float32}, A)
    @test C ≈ A ≈ BlockArray(A)
    @test eltype(C) == Float32

    C = convert(BlockArray{Float32, 2}, A)
    @test C ≈ A ≈ BlockArray(A)
    @test eltype(C) == Float32

    A = BlockArray(rand(2,3), [1,1], [2,1])
    C = convert(PseudoBlockArray, A)
    @test C == A == PseudoBlockArray(A)
    @test eltype(C) == eltype(A)

    C = convert(PseudoBlockArray{Float32}, A)
    @test C ≈ A ≈ PseudoBlockArray(A)
    @test eltype(C) == Float32

    C = convert(PseudoBlockArray{Float32, 2}, A)
    @test C ≈ A ≈ PseudoBlockArray(A)
    @test eltype(C) == Float32
end

@testset "string" begin
    A = BlockArray(rand(4, 5), [1,3], [2,3]);
    buf = IOBuffer()
    Base.showerror(buf, BlockBoundsError(A, (3,2)))
    @test String(take!(buf)) == "BlockBoundsError: attempt to access 2×2-blocked 4×5 BlockArray{Float64,2,Array{Float64,2}} at block index [3,2]"
end

if isdefined(Base, :flatten)
    flat = Base.flatten
else
    flat = Base.Iterators.flatten
end

replstrmime(x) = stringmime("text/plain", x)
@testset "replstring" begin
    @test replstrmime(BlockArray(collect(reshape(1:16, 4, 4)), [1,3], [2,2])) == "4×4 BlockArray{Int64,2,Array{Int64,2}}:\n 1  5  │   9  13\n ──────┼────────\n 2  6  │  10  14\n 3  7  │  11  15\n 4  8  │  12  16"
    design = zeros(Int16,6,9);
    A = BlockArray(design,[6],[4,5])
    @test replstrmime(A) == "6×9 BlockArray{Int16,2,Array{Int16,2}}:\n 0  0  0  0  │  0  0  0  0  0\n 0  0  0  0  │  0  0  0  0  0\n 0  0  0  0  │  0  0  0  0  0\n 0  0  0  0  │  0  0  0  0  0\n 0  0  0  0  │  0  0  0  0  0\n 0  0  0  0  │  0  0  0  0  0"
end

@testset "AbstractVector{Int} blocks" begin
    A = BlockArray(ones(6,6), 1:3, 1:3)
    @test A[1,1] == 1
    @test A[Block(2,3)] == ones(2,3)

    A = BlockArray(undef_blocks, Matrix{Float64}, 1:3, 1:3)
    A[Block(2,3)] = ones(2,3)
    @test A[Block(2,3)] == ones(2,3)
end

@testset "Block arithmetic" begin
    @test +(Block(1)) == Block(1)
    @test -(Block(1)) == Block(-1)
    @test Block(2) + Block(1) == Block(3)
    @test Block(2) + 1 == Block(3)
    @test 2 + Block(1) == Block(3)
    @test Block(2) - Block(1) == Block(1)
    @test Block(2) - 1 == Block(1)
    @test 2 - Block(1) == Block(1)
    @test 2*Block(1) == Block(2)
    @test Block(1)*2 == Block(2)

    @test isless(Block(1), Block(2))
    @test !isless(Block(1), Block(1))
    @test !isless(Block(2), Block(1))
    @test Block(1) < Block(2)
    @test Block(1) ≤ Block(1)
    @test Block(2) > Block(1)
    @test Block(1) ≥ Block(1)
    @test min(Block(1), Block(2)) == Block(1)
    @test max(Block(1), Block(2)) == Block(2)

    @test +(Block(1,2)) == Block(1,2)
    @test -(Block(1,2)) == Block(-1,-2)
    @test Block(1,2) + Block(2,3) == Block(3,5)
    @test Block(1,2) + 1 == Block(2,3)
    @test 1 + Block(1,2) == Block(2,3)
    @test Block(2,3) - Block(1,2) == Block(1,1)
    @test Block(1,2) - 1 == Block(0,1)
    @test 1 - Block(1,2) == Block(0,-1)
    @test 2*Block(1,2) == Block(2,4)
    @test Block(1,2)*2 == Block(2,4)

    @test isless(Block(1,1), Block(2,2))
    @test isless(Block(1,1), Block(2,1))
    @test !isless(Block(1,1), Block(1,1))
    @test !isless(Block(2,1), Block(1,1))
    @test Block(1,1) < Block(2,1)
    @test Block(1,1) ≤ Block(1,1)
    @test Block(2,1) > Block(1,1)
    @test Block(1,1) ≥ Block(1,1)
    @test min(Block(1,2), Block(2,2)) == Block(1,2)
    @test max(Block(1,2), Block(2,2)) == Block(2,2)

    @test convert(Int, Block(2)) == 2
    @test convert(Float64, Block(2)) == 2.0

    @test_throws MethodError convert(Int, Block(2,1))
    @test convert(Tuple{Int,Int}, Block(2,1)) == (2,1)
    @test convert(Tuple{Float64,Int}, Block(2,1)) == (2.0,1)
end

@testset "Strided array interface" begin
    A = PseudoBlockArray{Float64}(undef, 1:3, 1:3)
    fill!(A, 1)
    @test strides(A) == (1, size(A,1))
    x = randn(size(A,2))
    y = similar(x)
    @test BLAS.gemv!('N', 2.0, A, x, 0.0, y) ≈ 2A*x
end

@testset "lmul!/rmul!" begin
    A = PseudoBlockArray{Float64}(undef, 1:3)
    @test fill!(A, NaN) === A
    @test lmul!(0.0, A) === A
    @test Array(A) == zeros(6)
    @test fill!(A, NaN) === A
    @test rmul!(A, 0.0) === A
    @test Array(A) == zeros(6)

    A = BlockArray{Float64}(undef, 1:3)
    @test fill!(A, NaN) === A
    @test lmul!(0.0, A) === A
    @test Array(A) == zeros(6)
    @test fill!(A, NaN) === A
    @test rmul!(A, 0.0) === A
    @test Array(A) == zeros(6)
end

@testset "copy" begin
    A = PseudoBlockArray(randn(6), 1:3)
    B = copy(A)
    @test typeof(A) == typeof(B)
    @test blocksizes(A) == blocksizes(B)
    @test A == B
    B[1] = 2
    @test B[1] == 2
    @test A[1] ≠ 2

    A = BlockArray(randn(6), 1:3)
    B = copy(A)
    @test typeof(A) == typeof(B)
    @test blocksizes(A) == blocksizes(B)
    @test A == B
    B[1] = 2
    @test B[1] == 2
    @test A[1] ≠ 2
end


text:
using Test
using PandaModels; const _PdM = PandaModels
import PowerModels; const _PM = PowerModels

_PM.silence()

pdm_path = joinpath(dirname(pathof(PandaModels)), "..")
data_path = joinpath(pdm_path, "test", "data")

case_pm = joinpath(data_path, "test_pm.json")
case_pf_ac = joinpath(data_path, "test_pf_ac.json")
case_opf_ac = joinpath(data_path, "test_opf_ac.json")
case_opf_cl = joinpath(data_path, "test_opf_cl.json")
case_tnep_ac = joinpath(data_path, "test_tnep_ac.json")
case_ots_dc = joinpath(data_path, "test_ots_dc.json")
case_vd = joinpath(data_path, "test_vd.json")

@testset "PandaModels.jl" begin

        include("input.jl")
        include("call_powermodels.jl")
        include("call_pandamodels.jl")

end


text:
# Defines structus and functions relating to display.
# Defines `render` for all game objects.

function blendAlphaColors(x::SDL2.Color, y::SDL2.Color)
    xAlphaPercent = x.a / 255
    yAlphaPercent = y.a / 255
    r = round(x.r*xAlphaPercent + (y.r - Int32(x.r) * xAlphaPercent) * yAlphaPercent)
    g = round(x.g*xAlphaPercent + (y.g - Int32(x.g) * xAlphaPercent) * yAlphaPercent)
    b = round(x.b*xAlphaPercent + (y.b - Int32(x.b) * xAlphaPercent) * yAlphaPercent)
    a = round(x.a + (1-xAlphaPercent)*yAlphaPercent * 255)
    SDL2.Color(r,g,b,a)
end

""" Coordinate system to represent the screen, not the Game World.
Importantly, these are "upside-down" from WorldCoords: (0,0) is top-left, and
higher numbers go down and to the right.
"""
abstract type ScreenCoords end

""" ScreenPixelCoords are absolute space on screen, in actual pixels. """
struct ScreenPixelCoords <: ScreenCoords end

"""
    ScreenPixelPos(1200, 1150)
Absolute position on screen, in pixels. (0,0) is top-left of screen.
"""
struct ScreenPixelPos <: AbstractPos{ScreenPixelCoords}  # 0,0 == top-left
    x::Int
    y::Int
end
ScreenPixelPos(x::Number, y::Number) = ScreenPixelPos(convert.(Int, floor.((x,y)))...)

# TODO: CONSIDER SWITCHING TO DOWN-SCALING INSTEAD of up-scaling, so you can get
# more precise sizes (1.5 "pixels" -> 3 pixels). Or consider using Floats for UIPixelCoords.
""" UIPixelCoords are space on screen, in un-dpi-scaled "pixels". """
struct UIPixelCoords <: ScreenCoords end

"""
    UIPixelPos(400, 450)
Position on screen, in un-dpi-scaled "pixels". (0,0) is top-left of screen.

These should be used whenever placing anything on the screen, since they are
indepdent of resolution-scaling.
"""
struct UIPixelPos <: AbstractPos{UIPixelCoords}  # 0,0 == top-left
    x::Int
    y::Int
end
UIPixelPos(x::Number, y::Number) = UIPixelPos(convert.(Int, floor.((x,y)))...)

+(a::UIPixelPos, b::UIPixelPos) = UIPixelPos(a.x+b.x, a.y+b.y)
# TODO: 1.0?
+(a::UIPixelPos, x::Number) = UIPixelPos(a.x+x, a.y+x)

""" Absolute size on screen, in pixels. Use with ScreenPixelPos. """
struct ScreenPixelDims <: AbstractDims{ScreenPixelCoords}
    w::Int32  # Int32 to match SDL
    h::Int32
end
""" Size on screen, in un-dpi-scaled "pixels". Use with UIPixelPos. """
struct UIPixelDims <: AbstractDims{UIPixelCoords}
    w::Int
    h::Int
end

mutable struct Camera
    pos::WorldPos
    # Note, these are in WorldPos size.
    w::Threads.Atomic{Float32}   # Note: These are Atomics, since they can be modified by the
    h::Threads.Atomic{Float32}   # windowEventWatcher callback, which can run in another thread!
end
Camera() = Camera(WorldPos(0,0),100,100)

screenCenter() = UIPixelPos(winWidth[]/2, winHeight[]/2)
screenCenterX() = winWidth[]/2
screenCenterY() = winHeight[]/2
screenOffsetFromCenter(x::Int,y::Int) = UIPixelPos(screenCenterX()+x,screenCenterY()+y)

dpiScale() = winWidth_highDPI[] / winWidth[];
worldScale(c::Camera) = dpiScale() * (winWidth[] / cam.w[]);
function toScreenPos(p::WorldPos, c::Camera)
    scale = worldScale(c)
    ScreenPixelPos(
        round(winWidth_highDPI[]/2.0 + scale*(p.x-c.pos.x)), round(winHeight_highDPI[]/2.0 - scale*(p.y-cam.pos.y)))
end
function toScreenPos(p::UIPixelPos, c::Camera)
    scale = dpiScale()
    ScreenPixelPos(round(scale*p.x), round(scale*p.y))
end
toScreenPos(p::ScreenPixelPos, c::Camera) = p
function toWorldPos(p::ScreenPixelPos, c::Camera)
    scale = worldScale(c)
    WorldPos((p.x - c.w[]/2.)/scale, -(p.y - c.h[]/2.)/scale)
end
function toWorldPos(p::UIPixelPos, c::Camera)
    toWorldPos(toScreenPos(p, cam), cam)
end
toWorldPos(p::WorldPos, c::Camera) = p
function toUIPixelPos(p::ScreenPixelPos, c::Camera)
    scale = dpiScale()
    UIPixelPos(round(p.x/scale), round(p.y/scale))
end
function toUIPixelPos(p::WorldPos, c::Camera)
    toUIPixelPos(toScreenPos(p, cam), cam)
end
toUIPixelPos(p::UIPixelPos, c::Camera) = p
function toScreenPixelDims(dims::UIPixelDims,c::Camera)
    scale = dpiScale()
    ScreenPixelDims(round(scale*dims.w), round(scale*dims.h))
end
function toScreenPixelDims(dims::WorldDims,c::Camera)
    scale = worldScale(c)
    ScreenPixelDims(round(scale*dims.w), round(scale*dims.h))
end
toScreenPixelDims(dims::ScreenPixelDims,c::Camera) = dims
function toUIPixelDims(dims::ScreenPixelDims,c::Camera)
    scale = dpiScale()
    UIPixelDims(round(dims.w/scale), round(dims.h/scale))
end
function toUIPixelDims(d::WorldDims, c::Camera)
    toUIPixelDims(toScreenDims(d, cam), cam)
end
toUIPixelDims(dims::UIPixelDims,c::Camera) = dims
function toWorldDims(dims::ScreenPixelDims,c::Camera)
    scale = worldScale(c)
    WorldDims(dims.w/scale, dims.h/scale)
end
function toWorldDims(d::UIPixelDims, c::Camera)
    toWorldDims(toScreenDims(d, cam), cam)
end
toWorldDims(dims::WorldDims,c::Camera) = dims

SetRenderDrawColor(renderer::Ptr{SDL2.Renderer}, c::SDL2.Color) = SDL2.SetRenderDrawColor(
    renderer, Int64(c.r), Int64(c.g), Int64(c.b), Int64(c.a))

# Convenience functions to allow `WorldPos(x,y)...` to become `x,y`
Base.iterate(p::AbstractPos, i=1) = if i==1 return (p.x,2) elseif i==2 return (p.y,3) else nothing end
Base.iterate(p::AbstractDims, i=1) = if i==1 return (p.w,2) elseif i==2 return (p.h,3) else nothing end
# To allow .+ (but doesn't work)
#Base.length(::Union{AbstractPos,AbstractDims}) = 2

topLeftPos(center::P, dims::D) where P<:AbstractPos{Coord} where D<:AbstractDims{Coord} where Coord<:ScreenCoords = P(center.x - dims.w/2.0, center.y - dims.h/2.0)  # positive is down
topLeftPos(center::P, dims::D) where P<:AbstractPos{Coord} where D<:AbstractDims{Coord} where Coord<:WorldCoords = P(center.x - dims.w/2.0, center.y + dims.h/2.)  # positive is up
rectOrigin(center::P, dims::D) where P<:AbstractPos{C} where D<:AbstractDims{C} where {C} = P(center.x - dims.w/2.0, center.y - dims.h/2.0)  # always minus

#makeSDL2RectFromOrigin(p::ScreenPixelPos, d::ScreenPixelDims) = SDL2.Rect(p..., d...)
#makeSDL2RectFromCenter(p::ScreenPixelPos, d::ScreenPixelDims) = SDL2.Rect(rectOrigin(p,d)..., d...)
#makeSDL2RectFromOrigin(p,d) = makeSDL2RectFromCenter(toScreenPos(p), toScreenPixelDims(d))
#makeSDL2RectFromCenter(p,d) =
function renderRectCentered(cam, renderer, center::AbstractPos{C}, dims::AbstractDims{C}, color; outlineColor=nothing) where C
    origin = topLeftPos(center, dims)
    renderRectFromOrigin(cam, renderer, origin, dims, color; outlineColor=outlineColor)
end
function renderRectFromOrigin(cam, renderer, origin::AbstractPos{C}, dims::AbstractDims{C}, color; outlineColor=nothing) where C
    screenPos = toScreenPos(origin, cam)
    rect = SDL2.Rect(screenPos.x, screenPos.y, toScreenPixelDims(dims, cam)...)
    if color != nothing
        SetRenderDrawColor(renderer, color)
        SDL2.RenderFillRect(renderer, Ref(rect) )
    end
    if outlineColor != nothing
        SetRenderDrawColor(renderer, outlineColor)
        SDL2.RenderDrawRect(renderer, Ref(rect) )
    end
end

function renderProgressBar(percent, cam::Camera, renderer, center::AbstractPos{C}, dims::D, color, bgColor, boxColor) where D<:AbstractDims{C} where C
    # bg
    renderRectCentered(cam, renderer, center, dims, bgColor)
    # progress
    origin = topLeftPos(center, dims)
    renderRectFromOrigin(cam, renderer, origin, D(round(percent * dims.w), dims.h), color)
    # outline
    renderRectCentered(cam, renderer, center, dims, nothing; outlineColor=boxColor)
end
function renderUnit(o::UnitTypes, playerColor, cam::Camera, renderer, dims::WorldDims, color)
    # First render the player color, then the unit color (for transparency)
    renderRectCentered(cam, renderer, o.pos, dims, playerColor)
    renderRectCentered(cam, renderer, o.pos, dims, color)

    # render health bar
    healthBarPos = WorldPos(o.pos.x, o.pos.y + dims.h/2 + healthBarRenderOffset)
    healthBarDims = WorldDims(healthBarRenderWidth, healthBarRenderHeight)
    renderProgressBar(health_percent(o), cam, renderer, healthBarPos,
          healthBarDims, healthBarColor, kBuildOpsBgColor,
          if toScreenPixelDims(healthBarDims, cam).h <= 2 # pixels
              nothing
          else
              healthBarOutlineColor
          end)
end
function render(o::Collector, playerColor, cam::Camera, renderer)
    unitColor = blendAlphaColors(playerColor, kCollectorColor)
    dims = WorldDims(collectorRenderWidth, collectorRenderWidth)
    tl = topLeftPos(o.pos, dims)
    for y in [tl.y, tl.y - dims.h]
        for x in tl.x+1 : 5 : tl.x-1 + dims.w
            nudge = Vector2D(rand(0.0:0.01:2.0,2)...)
            renderRectCentered(cam, renderer, WorldPos(x,y)+nudge, WorldDims(2,10), unitColor)
        end
    end
    for x in [tl.x, tl.x + dims.w]
        for y in tl.y-1 : -5 : tl.y+1 - dims.h
            nudge = Vector2D(rand(0.0:0.01:2.0,2)...)
            renderRectCentered(cam, renderer, WorldPos(x,y)+nudge, WorldDims(10,2), unitColor)
        end
    end
    renderUnit(o, playerColor, cam, renderer, dims, kCollectorColor)
end
function render(o::Fighter, playerColor, cam::Camera, renderer)
    dims = WorldDims(unitRenderWidth, unitRenderWidth)
    renderUnit(o, playerColor, cam, renderer, dims, kFighterColor)
end

abstract type AbstractButton end
mutable struct MenuButton <: AbstractButton
    enabled::Bool
    pos::UIPixelPos
    dims::UIPixelDims
    text::String
    callBack
end
mutable struct KeyButton <: AbstractButton
    enabled::Bool
    pos::UIPixelPos
    dims::UIPixelDims
    text::String
    callBack
end

mutable struct CheckboxButton
    toggled::Bool
    button::MenuButton
end

import Base.run
run(b::AbstractButton) = b.callBack()
function run(b::CheckboxButton)
    b.toggled = !b.toggled
    b.button.callBack(b.toggled)
end

# pointwise subtraction with bounds checking (floors to 0)
-(a::SDL2.Color, b::Int) = SDL2.Color(a.r-min(b,a.r), a.g-min(b,a.g), a.b-min(b,a.b), a.a-min(b,a.a))
SDL2.Color(1,5,1,1) - 2 == SDL2.Color(0,3,0,0)

function render(b::AbstractButton, cam::Camera, renderer, color, fontSize)
    if (!b.enabled)
         return
    end
    topLeft = topLeftPos(b.pos, b.dims)
    screenPos = toScreenPos(topLeft, cam)
    rect = SDL2.Rect(screenPos..., toScreenPixelDims(b.dims, cam)...)
    x,y = Int[0], Int[0]
    SDL2.GetMouseState(pointer(x), pointer(y))
    if clickedButton == b
        if mouseOnButton(UIPixelPos(x[],y[]),b,cam)
            color = color - 50
        else
            color = color - 30
        end
    else
        if mouseOnButton(UIPixelPos(x[],y[]),b,cam)
            color = color - 10
        end
    end
    SetRenderDrawColor(renderer, color)
    SDL2.RenderFillRect(renderer, Ref(rect) )
    renderText(renderer, cam, b.text, b.pos; fontSize = fontSize)
end

function render(b::MenuButton, cam::Camera, renderer)
    render(b, cam, renderer, kMenuButtonColor, kMenuButtonFontSize)
end
function render(b::KeyButton, cam::Camera, renderer)
    render(b, cam, renderer, kKeySettingButtonColor, kKeyButtonFontSize)
end

function render(b::CheckboxButton, cam::Camera, renderer)
    # Hack: move button text offcenter before rendering to accomodate checkbox
    offsetText = " "
    text_backup = b.button.text
    b.button.text = offsetText * b.button.text
    render(b.button, cam, renderer)
    b.button.text = text_backup

    # Render checkbox
    render_checkbox_square(b.button, 6, SDL2.Color(200,200,200, 255), cam, renderer)

    if b.toggled
        # Inside checkbox "fill"
        render_checkbox_square(b.button, 8, SDL2.Color(100,100,100, 255), cam, renderer)
    end
end

function render_checkbox_square(b::AbstractButton, border, color, cam, renderer)
    checkbox_radius = b.dims.h/2.0 - border  # (checkbox is a square)
    topLeft = topLeftPos(b.pos, b.dims)
    topLeft = topLeft + border
    screenPos = toScreenPos(topLeft, cam)
    screenDims = toScreenPixelDims(UIPixelDims(checkbox_radius*2, checkbox_radius*2), cam)
    rect = SDL2.Rect(screenPos..., screenDims...)
    SetRenderDrawColor(renderer, color)
    SDL2.RenderFillRect(renderer, Ref(rect) )
end

# ---- Text Rendering ----

fonts_cache = Dict()
txt_cache = Dict()
function sizeText(txt, fontName, fontSize)
    sizeText(txt, loadFont(dpiScale(), fontName, fontSize))
end
function sizeText(txt, font::Ptr{SDL2.TTF_Font})
   fw,fh = Cint[1], Cint[1]
   SDL2.TTF_SizeText(font, txt, pointer(fw), pointer(fh))
   return ScreenPixelDims(fw[1],fh[1])
end
function loadFont(scale, fontName, fontSize)
   fontSize = scale*fontSize
   fontKey = (fontName, Cint(round(fontSize)))
   if haskey(fonts_cache, fontKey)
       font = fonts_cache[fontKey]
   else
       font = SDL2.TTF_OpenFont(fontKey...)
       font == C_NULL && throw(ErrorException("Failed to load font '$fontKey'"))
       fonts_cache[fontKey] = font
   end
   return font
end
function createText(renderer, cam, txt, fontName, fontSize)
   font = loadFont(dpiScale(), fontName, fontSize)
   txtKey = (font, txt)
   if haskey(txt_cache, txtKey)
       tex = txt_cache[txtKey]
   else
       text = SDL2.TTF_RenderText_Blended(font, txt, SDL2.Color(20,20,20,255))
       tex = SDL2.CreateTextureFromSurface(renderer,text)
       #SDL2.FreeSurface(text)
       txt_cache[txtKey] = tex
   end

   fw,fh = Cint[1], Cint[1]
   SDL2.TTF_SizeText(font, txt, pointer(fw), pointer(fh))
   fw,fh = fw[1],fh[1]

   return tex, ScreenPixelDims(fw, fh)
end
@enum TextAlign centered leftJustified rightJustified
function renderText(renderer, cam::Camera, txt::String, pos::UIPixelPos
                    ; fontName = defaultFontName,
                     fontSize=defaultFontSize, align::TextAlign = centered)
   tex, fDims = createText(renderer, cam, txt, fontName, fontSize)
   renderTextSurface(renderer, cam, pos, tex, toUIPixelDims(fDims,cam), align)
end

function renderTextSurface(renderer, cam::Camera, pos::AbstractPos{C},
                           tex::Ptr{SDL2.Texture}, dims::AbstractDims{C}, align::TextAlign) where C
   screenPos = toScreenPos(pos, cam)
   screenDims = toScreenPixelDims(dims, cam)
   x,y, fw, fh = screenPos.x, screenPos.y, screenDims.w, screenDims.h
   renderPos = SDL2.Rect(0,0,0,0)
   if align == centered
       renderPos = SDL2.Rect(Int(floor(x-fw/2.0)), Int(floor(y-fh/2.0)), fw,fh)
   elseif align == leftJustified
       renderPos = SDL2.Rect(Int(floor(x)), Int(floor(y-fh/2.0)), fw,fh)
   else # align == rightJustified
       renderPos = SDL2.Rect(Int(floor(x-fw)), Int(floor(y-fh/2.0)), fw,fh)
   end
   SDL2.RenderCopy(renderer, tex, C_NULL, pointer_from_objref(renderPos))
   #SDL2.DestroyTexture(tex)
end

function hcat_render_text(lines, renderer, cam, gap, pos::UIPixelPos;
         fixedWidth=nothing, fontName=defaultFontName, fontSize=defaultFontSize)
    numLines = size(lines)[1]
    if fixedWidth != nothing
        widths = fill(fixedWidth,size(lines))
    else
        widths = [toUIPixelDims(sizeText(line, fontName, fontSize), cam).w for line in lines]
    end
    totalWidth = sum(widths) + gap*(numLines-1)
    runningWidth = 0
    leftMostPos = pos.x - totalWidth/2.0
    text_centers = []
    for i in 1:numLines
        linePos = leftMostPos + runningWidth
        leftPos = UIPixelPos(linePos, pos.y)
        renderText(renderer, cam, lines[i], leftPos; fontName=fontName, fontSize=fontSize, align=leftJustified)
        runningWidth += widths[i] + gap
        push!(text_centers, UIPixelPos(leftPos.x + widths[i]÷2, leftPos.y))
    end
    return text_centers
end

#  ------- Image rendering ---------

function render(t::Ptr{SDL2.Texture}, pos::AbstractPos{C}, cam::Camera, renderer; size::Union{Cvoid, AbstractDims{C}} = nothing) where C
    if (t == C_NULL) return end
    pos = toScreenPos(pos, cam)
    if size != nothing
        size = toScreenPixelDims(size, cam)
        w = size.w
        h = size.h
    else
        w,h,access = Cint[1], Cint[1], Cint[1]
        format = Cuint[1]
        SDL2.QueryTexture( t, format, access, w, h );
        w,h = w[], h[]
    end
    rect = SDL2.Rect(pos.x - w÷2,pos.y - h÷2,w,h)
    SDL2.RenderCopy(renderer, t, C_NULL, pointer_from_objref(rect))
end


text:
using KeplerTools
using ProfileView

include(joinpath(@__DIR__, "..", "data", "kerbol_system.jl"))
rand_orbits = Orbit[]
rand_start_times = Float64[]
rand_end_times = Float64[]
for i=1:1000000
    orb = Orbit(moho.orbit.a + 3*duna.orbit.a*rand(),
                rand(),
                π*rand(),
                2π*rand(),
                2π*rand(),
                2π*rand(),
                1e6*rand(),
                sun
    )
    starttime = orb.epoch + 1e6*rand()
    endtime = starttime + orb.period*rand()
    push!(rand_orbits, orb)
    push!(rand_start_times, starttime)
    push!(rand_end_times, endtime)
end

function prof_p_lambert(rand_orbits::Vector{Orbit}, rand_start_times, rand_end_times, N=1000000)
    for i=1:N
        torb = p_lambert(rand_orbits[i], rand_orbits[i], rand_start_times[i], rand_end_times[i])
    end
end

ProfileView.@profview prof_p_lambert(rand_orbits, rand_start_times, rand_end_times)


text:
function _clear_wd(gl_repo)
    for path in _readdir(gl_repo; join = true)
        endswith(path, ".git") && continue
        _rm(path)
    end
end

function hard_pull(gl::GitLink; verbose = false, clearwd = true, tries = 1)

    for t in 1:max(tries, 1)

        ignorestatus = true
        gl_repo = repo_dir(gl)
        url = remote_url(gl)
        mkpath(gl_repo)

        # check remote connection
        rhash = _remote_HEAD_hash(url)
        isempty(rhash) && continue

        # enforce repo format (lazy method)
        _format_repo!(gl_repo, url; verbose)

        # clear wd
        clearwd && _clear_wd(gl_repo)

        if _check_gitdir(gl_repo)
            # if all ok
            _run("git -C $(gl_repo) fetch 2>&1"; verbose, ignorestatus)
            _run("git -C $(gl_repo) reset --hard FETCH_HEAD 2>&1"; verbose, ignorestatus)
        else
            # clone
            _rm(gl_repo)
            mkpath(gl_repo)
            _run("git -C $(gl_repo) clone --depth=1 $(url) $(gl_repo) 2>&1"; verbose, ignorestatus)
            _run("git -C $(gl_repo) fetch 2>&1"; verbose, ignorestatus)
            _run("git -C $(gl_repo) reset --hard FETCH_HEAD 2>&1"; verbose, ignorestatus)
        end

        # check success
        chash = _HEAD_hash(gl_repo)
        if rhash != chash
            _rm(gl_repo) # something fail
            continue
        end
        return true
    end

    return false
end


text:
import ArchGDAL, GDAL;
const AG = ArchGDAL
using Test

@testset "test_gdalutilities_errors.jl" begin
    AG.read("data/utmsmall.tif") do ds_small
        @testset "GDAL Error" begin
            @test_throws GDAL.GDALError AG.gdalinfo(
                ds_small,
                ["-novalidoption"],
            )
            @test_throws GDAL.GDALError AG.unsafe_gdaltranslate(
                ds_small,
                ["-novalidoption"],
            )
            @test_throws GDAL.GDALError AG.unsafe_gdalbuildvrt(
                [ds_small],
                ["-novalidoption"],
            )
            @test_throws GDAL.GDALError AG.unsafe_gdaldem(
                ds_small,
                "hillshade",
                ["-novalidoption"],
            )
            @test_throws GDAL.GDALError AG.unsafe_gdalnearblack(
                ds_small,
                ["-novalidoption"],
            )
            @test_throws GDAL.GDALError AG.unsafe_gdalwarp(
                [ds_small],
                ["-novalidoption"],
            )
        end
    end

    @testset "Interactive data/utmsmall.tif" begin
        ds_small = AG.read("data/utmsmall.tif")
        @test_throws GDAL.GDALError AG.gdalinfo(ds_small, ["-novalidoption"])
        @test_throws GDAL.GDALError AG.unsafe_gdaltranslate(
            ds_small,
            ["-novalidoption"],
        )
        @test_throws GDAL.GDALError AG.unsafe_gdalbuildvrt(
            [ds_small],
            ["-novalidoption"],
        )
        @test_throws GDAL.GDALError AG.unsafe_gdaldem(
            ds_small,
            "hillshade",
            ["-novalidoption"],
        )
        @test_throws GDAL.GDALError AG.unsafe_gdalnearblack(
            ds_small,
            ["-novalidoption"],
        )
        @test_throws GDAL.GDALError AG.unsafe_gdalwarp(
            [ds_small],
            ["-novalidoption"],
        )
    end
end


text:
using ClangCompiler.JLLEnvs
using Test
using Base.BinaryPlatforms

@testset "JLLEnv" begin
    p = HostPlatform()
    @info p
    @info os(p), arch(p), libc(p)
    @test_nowarn JLLEnvs.get_default_env(; is_cxx=true, version=v"7.0.1")
end


text:
using MOCNeutronTransport
@testset "Triangle6_3D" begin
    for T in [Float32, Float64]
        @testset "Constructors" begin
            p₁ = Point_3D(T, 0)
            p₂ = Point_3D(T, 1)
            p₃ = Point_3D(T, 1, 1)
            p₄ = Point_3D(T, 1//2)
            p₅ = Point_3D(T, 1, 1//2)
            p₆ = Point_3D(T, 1//2, 1//2)
            tri6 = Triangle6_3D((p₁, p₂, p₃, p₄, p₅, p₆))
            @test tri6.points == (p₁, p₂, p₃, p₄, p₅, p₆)

            # single constructor
            tri6 = Triangle6_3D(p₁, p₂, p₃, p₄, p₅, p₆)
            @test tri6.points == (p₁, p₂, p₃, p₄, p₅, p₆)
        end

        @testset "Methods" begin
            p₁ = Point_3D(T, 0)
            p₂ = Point_3D(T, 1)
            p₃ = Point_3D(T, 1, 1)
            p₄ = Point_3D(T, 1//2)
            p₅ = Point_3D(T, 1, 1//2)
            p₆ = Point_3D(T, 1//2, 1//2)
            tri6 = Triangle6_3D((p₁, p₂, p₃, p₄, p₅, p₆))

            # interpolation
            tri6(0, 0) ≈ p₁
            tri6(1, 0) ≈ p₂
            tri6(0, 1) ≈ p₃
            tri6(1//2, 0) ≈ p₄
            tri6(1//2, 1//2) ≈ p₅
            tri6(0, 1//2) ≈ p₆
            tri6(1//2, 1//2) ≈ Point_3D(T, 1//2, 1//2)

            # area
            p₁ = Point_3D(T, 0)
            p₂ = Point_3D(T, 2)
            p₃ = Point_3D(T, 2, 2)
            p₄ = Point_3D(T, 1, 1//4)
            p₅ = Point_3D(T, 3, 1)
            p₆ = Point_3D(T, 1, 1)
            tri6 = Triangle6_3D((p₁, p₂, p₃, p₄, p₅, p₆))

            # 3D default
            p₂ = Point_3D(T, 2, 0, 3)
            tri6 = Triangle6_3D((p₁, p₂, p₃, p₄, p₅, p₆))
            @test isapprox(area(tri6; N = 79), 6.328781460309, atol=1.0e-6)

            # real_to_parametric
            p₁ = Point_3D(T, 0)
            p₂ = Point_3D(T, 2)
            p₃ = Point_3D(T, 2, 2)
            p₄ = Point_3D(T, 1, 1//4)
            p₅ = Point_3D(T, 3, 1)
            p₆ = Point_3D(T, 1, 1)
            tri6 = Triangle6_3D((p₁, p₂, p₃, p₄, p₅, p₆))
            p = real_to_parametric(p₁, tri6)
            @test tri6(p[1], p[2]) ≈ p₁
            p = real_to_parametric(p₂, tri6)
            @test tri6(p[1], p[2]) ≈ p₂
            p = real_to_parametric(p₃, tri6)
            @test tri6(p[1], p[2]) ≈ p₃
            p = real_to_parametric(p₄, tri6)
            @test tri6(p[1], p[2]) ≈ p₄
            p = real_to_parametric(p₅, tri6)
            @test tri6(p[1], p[2]) ≈ p₅
            p = real_to_parametric(p₆, tri6)
            @test tri6(p[1], p[2]) ≈ p₆

            # intersect
            p₁ = Point_3D(T, 0)
            p₂ = Point_3D(T, 2)
            p₃ = Point_3D(T, 2, 2)
            p₄ = Point_3D(T, 1, 1//4)
            p₅ = Point_3D(T, 3, 1)
            p₆ = Point_3D(T, 1, 1)
            tri6 = Triangle6_3D((p₁, p₂, p₃, p₄, p₅, p₆))
            l = LineSegment_3D(Point_3D(T, 1, 1//2, -2),
                               Point_3D(T, 1, 1//2,  2))
            intersection = l ∩ tri6

            # 1 intersection
            @test intersection[1]
            @test intersection[2] == 1
            @test intersection[3] ≈ Point_3D(T, 1, 1//2, 0)

            # 0 intersection
            l = LineSegment_3D(Point_3D(T, 1, 0, -2),
                               Point_3D(T, 1, 0,  2))
            intersection = l ∩ tri6
            @test !intersection[1]
            @test intersection[2] == 0

            # 2 intersections
            p₂ = Point_3D(T, 2, 0, 3)
            tri6 = Triangle6_3D((p₁, p₂, p₃, p₄, p₅, p₆))

            l = LineSegment_3D(Point_3D(T, 1, -2, -1//5),
                               Point_3D(T, 1,  2, -1//5))
            intersection = l ∩ tri6
            @test intersection[1]
            @test intersection[2] == 2
            @test norm(intersection[3] - Point_3D(T, 1, 0.3116568318033659, -1//5)) < 1e-3
            @test norm(intersection[4] - Point_3D(T, 1, 0.7706348348626366, -1//5)) < 1e-3
        end
    end
end

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
function factor_KKT!(ws::COSMO.Workspace{T}) where {T <: AbstractFloat}
	settings = ws.settings
	p = ws.p
	if length(ws.M) == 0
		# KKT matrix
		ws.M = [p.P + settings.sigma*I SparseMatrixCSC(p.A'); p.A -I]
	end
	n = p.model_size[2]
	m = p.model_size[1]
	@inbounds @simd for i = (n + 1):(n + m)
		ws.M[i, i] = -one(T) / ws.ρvec[i - n]
	end

	# Do LDLT Factorization: A = LDL^T
	if settings.verbose_timing
		ws.times.factor_time += @elapsed ws.F = ldlt(ws.M)
	else
		ws.F = ldlt(ws.M)
	end
	return nothing
end


text:
"""
Calculate distances and distributions of logical operator and stabilizer weights for a
[`TensorNetworkCode`](@ref), see [`tn_distance`](@ref) and [`tn_operator_weights`](@ref).
"""
module TNDistance

#imports
using ..TensorNetworkCodes: TensorNetworkCode
using ..TensorNetworkCodes: all_cosets, code_to_Itensor, identity_coset
using ..TensorNetworkCodes: nodes, node_indices, node_types
using ..TensorNetworkCodes: num_qubits
using ..TensorNetworkCodes: pauli_weight
using ITensors: Index, ITensor, array, contract, dim, hastags, inds
using StatsPlots: groupedbar

#exports
export OperatorWeights
export plot_operator_weights
export tn_distance, tn_operator_weights

"""
    OperatorWeights(stabilizer_weights,all_operator_weights)

Stores operator weight distribution for stabilizers and all code operators
(stabilizers + all logical representatives).  Also stores distance.  Fields:

    stabilizer_weights::Array{Int64,1}
    all_operator_weights::Array{Int64,1}
    distance::Int64

Note that since julia is 1-indexed `stabilizer_weights[j]` gives the number of
stabilizers of weight j-1.
"""
struct OperatorWeights
    stabilizer_weights::Array{Int64,1}
    all_operator_weights::Array{Int64,1}
    distance::Int64

    function OperatorWeights(stabilizer_weights,all_operator_weights)
        logical_weights = all_operator_weights - stabilizer_weights
        index = findfirst(x->x>0,logical_weights)

        if index == nothing
            distance = 0
        else
            distance =  index -1
        end
        new(stabilizer_weights,all_operator_weights,distance)
    end
end
OperatorWeights() = OperatorWeights(Int64[0],Int64[1])

"""
    simple_weight_tensor(physical_indices,out_index) -> ITensor

Returns a tensor that, when contracted with a code tensor, gives the weights
of all Pauli operators described by that code tensor.
"""
function _simple_weight_tensor(
        physical_indices::Array{Index{Int64},1},
        out_index::Index{Int64})

    all_indices = vcat(physical_indices,[out_index])
    dims = dim.(all_indices)
    array = fill(0,dims...)

    n = length(physical_indices)

    for γ in 0:4^n-1
        index = digits!(zeros(Int64,n),γ,base=4)
        w = pauli_weight(index)
        index = index .+ 1

        full_index = vcat(index,w+1)
        array[full_index...] = 1
    end

    return ITensor(array,all_indices...)
end

"""
    _addition_tensor(in_index1,in_index2,out_index)->ITensor

Gives a tensor that is one if `out_index` equals  the sum of `in_index1`
and `in_index2`.  Warning: since Julia is one indexed, zero is replaced by
1, so that the tensor is nonzero if `in_index1 + in_index2 - 1 = out_index`.
"""
function _addition_tensor(in_index1,in_index2,out_index)

    in_dim1 = dim(in_index1)
    in_dim2 = dim(in_index2)
    out_dim = dim(out_index)
    dims = [in_dim1,in_dim2,out_dim]

    array = fill(0,dims...)

    for α in 1:in_dim1, β in 1:in_dim2, γ in 1:out_dim
        if α + β - 2 == γ - 1
            array[α,β,γ] = 1
        end
    end

    return ITensor(array,in_index1,in_index2,out_index)
end

function _addition_tensor(in_index1,in_index2)

    in_dim1 = dim(in_index1)
    in_dim2 = dim(in_index2)
    out_dim = in_dim1 + in_dim2 -2 + 1  #(indices start from zero in Julia)
    out_index = Index(out_dim,"out_index")

    return _addition_tensor(in_index1,in_index2,out_index)
end

"""
    _tn_weights(code::TensorNetworkCode;cosets=all_cosets,
    truncate_to=num_qubits(code)+1) -> ITensor

Returns a tensor describing the weights of all operators in the `cosets` (either
`all_cosets` or `identity_coset`) of
the code.  (`cosets` is applied to each `seed_code` separately.)
"""
function _tn_weights(
        code::TensorNetworkCode;
        cosets=all_cosets,
        truncate_to=num_qubits(code)+1)

    code_tensors = ITensor[]

    for node in nodes(code)
        if node < 0
            code_type = node_types(code,node)
            seed_code = code.seed_codes[code_type]

            k = length(seed_code.logicals)/2
            logical_indices = [Index(4,"logical") for _ in 1:k]
            physical_indices = node_indices(code,node)
            code_tensor = code_to_Itensor(seed_code,logical_indices,physical_indices)

            code_tensor = cosets(code_tensor)

            push!(code_tensors,code_tensor)
        end
    end

    for α in 1:length(code_tensors)
        indices = inds(code_tensors[α])  #tuple
        physical_indices = filter(x->hastags(x,"physical"),indices)
        num_phys = length(physical_indices)
        output_index = Index(num_phys + 1,"weight")

        local_weight_tensor = _simple_weight_tensor([physical_indices...],output_index)
        code_tensors[α] = code_tensors[α]*local_weight_tensor
    end

    indices = inds(code_tensors[1])
    β = findfirst(x->hastags(x,"weight"),indices)
    in_index1 = indices[β]

    for α in 2:length(code_tensors)
        indices = inds(code_tensors[α])
        β = findfirst(x->hastags(x,"weight"),indices)
        in_index2 = indices[β]

        in_dim1 = dim(in_index1)
        in_dim2 = dim(in_index2)
        total_out_dim = in_dim1 + in_dim2 -2 + 1

        if total_out_dim <= truncate_to
            out_dim = total_out_dim
        else
            out_dim = truncate_to
        end

        out_index = Index(out_dim,"outind")

        add_tensor = _addition_tensor(in_index1,in_index2,out_index)
        code_tensors[α] = code_tensors[α] * add_tensor

        in_index1 = out_index
    end

    return contract(code_tensors)

end

"""
    tn_operator_weights(code::TensorNetworkCode;truncate_to=num_qubits(code)+1)
    -> OperatorWeights

Returns [`OperatorWeights`](@ref) which includes the number of stabilizers
of each weight, as well as the number of logical representatives
(excluding identity) of each weight, as well as the code distance.

# Examples
```jldoctest
julia> using TensorNetworkCodes.TNDistance

julia> code = TensorNetworkCode(five_qubit_code());

julia> code = contract(code,code,[[1,1],[3,3]]);

julia> tn_operator_weights(code)
OperatorWeights([1, 0, 0, 0, 9, 0, 6], [1, 0, 9, 24, 99, 72, 51], 2)
```
"""
function tn_operator_weights(
        code::TensorNetworkCode;
        truncate_to=num_qubits(code)+1)

    stabilizer_weights = array(_tn_weights(
            code,
            cosets=identity_coset,
        truncate_to=truncate_to))

    all_operator_weights = array(_tn_weights(
            code,
            cosets=all_cosets,
            truncate_to=truncate_to))

    return OperatorWeights(stabilizer_weights,all_operator_weights)
end

"""
    tn_distance(code::TensorNetworkCode;truncate_to=num_qubits(code)+1)
    -> Int64

Returns the code distance calculated by contracting a tensor network.  Returns
zero for stabilizer states (codes with no logical qubits).

# Examples
```jldoctest
julia> using TensorNetworkCodes.TNDistance

julia> code = TensorNetworkCode(steane_code());

julia> code = contract(code,code,[[1,1],[3,3]]); # contract two copies

julia> tn_distance(code) # this code has poor distance!
2
```
"""
function tn_distance(
        code::TensorNetworkCode;
        truncate_to=num_qubits(code)+1)

    op_weights = tn_operator_weights(
        code::TensorNetworkCode;
        truncate_to=truncate_to)

    return op_weights.distance
end

"""
    plot_operator_weights(
        weights::OperatorWeights; truncate_to=length(weights.stabilizer_weights)
    )

Plots a bar plot (on a log scale) of the number of operators of each weight.
This includes all code operators (logicals plus stabilizers) and all stabilizers plotted
separately.
"""
function plot_operator_weights(
    weights::OperatorWeights;
    truncate_to=length(weights.stabilizer_weights)
)

    all_ops = weights.all_operator_weights[1:truncate_to]
    stabs = weights.stabilizer_weights[1:truncate_to]

    max_y_val = maximum(all_ops)

    logical_weights = all_ops - stabs

    # for log plot, so zero doesn't give -infinity
    logical_weights = logical_weights .+ 0.001
    all_ops = all_ops .+ 0.001
    stabs = stabs .+ 0.001

    bar_data = hcat(
        stabs,
        logical_weights)

    L = length(logical_weights)

    labels = repeat(
        ["Stabilizers", "Logicals"],
        inner = L)

    return groupedbar(
        repeat(collect(0:L-1),outer = 2),
        bar_data,
        bar_position = :dodge,
        group = labels,
        xlabel = "Weight",
        ylabel = "Number",
        yaxis=:log,
        title = "Number of operators of given weight",
        bar_width = 0.67,
        lw = 0,
        framestyle = :box,
        legend = :topleft,
        ylims = (0.1,10*max_y_val),
        xticks = (0:L-1))
end

end


text:
## Effect differences ##

const EFFECT_FUNCS = Dict{Symbol,Function}()

"Convert effect formula to a state difference (additions, deletions, etc.)"
function effect_diff(domain::Domain, state::State, effect::Term)
    effect_fn! = get(EFFECT_FUNCS, effect.name, add_effect!)
    return effect_fn!(GenericDiff(), domain, state, effect)
end

function add_effect!(diff::GenericDiff,
                     domain::Domain, state::State, effect::Term)
    if effect isa Compound # Evaluated all nested functions
        args = Term[a isa Compound ? Const(evaluate(domain, state, a)) : a
                    for a in effect.args]
        effect = Compound(effect.name, args)
    end
    push!(diff.add, effect)
    return diff
end

function del_effect!(diff::GenericDiff,
                     domain::Domain, state::State, effect::Term)
    effect = effect.args[1]
    if effect isa Compound # Evaluated all nested functions
        args = Term[a isa Compound ? Const(evaluate(domain, state, a)) : a
                    for a in effect.args]
        effect = Compound(effect.name, args)
    end
    push!(diff.del, effect)
    return diff
end
EFFECT_FUNCS[:not] = del_effect!

function and_effect!(diff::GenericDiff,
                     domain::Domain, state::State, effect::Term)
    for eff in get_args(effect)
        combine!(diff, effect_diff(domain, state, eff))
    end
    return diff
end
EFFECT_FUNCS[:and] = and_effect!

function when_effect!(diff::GenericDiff,
                      domain::Domain, state::State, effect::Term)
    cond, eff = effect.args[1], effect.args[2]
    if !satisfy(domain, state, cond) return diff end
    return combine!(diff, effect_diff(domain, state, eff))
end
EFFECT_FUNCS[:when] = when_effect!

function forall_effect!(diff::GenericDiff,
                        domain::Domain, state::State, effect::Term)
    cond, eff = effect.args[1], effect.args[2]
    for s in satisfiers(domain, state, cond)
        combine!(diff, effect_diff(domain, state, substitute(eff, s)))
    end
    return diff
end
EFFECT_FUNCS[:forall] = forall_effect!

function assign_effect!(diff::GenericDiff,
                        domain::Domain, state::State, effect::Term)
    term, val = effect.args[1], effect.args[2]
    diff.ops[term] = val
    return diff
end
EFFECT_FUNCS[:assign] = assign_effect!

function modify_effect!(diff::GenericDiff,
                        domain::Domain, state::State, effect::Term)
    term, val = effect.args[1], effect.args[2]
    op = GLOBAL_MODIFIERS[effect.name]
    diff.ops[term] = Compound(op, Term[term, val])
    return diff
end
for name in keys(GLOBAL_MODIFIERS)
    EFFECT_FUNCS[name] = modify_effect!
end

## Precondition differences ##

const PRECOND_FUNCS = Dict{Symbol,Function}()

"Convert precondition formula to a state difference."
function precond_diff(domain::Domain, state::State, precond::Term)
    precond_fn! = get(PRECOND_FUNCS, precond.name, add_precond!)
    return precond_fn!(GenericDiff(), domain, state, precond)
end

function add_precond!(diff::GenericDiff,
                      domain::Domain, state::State, precond::Term)
    push!(diff.add, precond)
    return diff
end

function del_precond!(diff::GenericDiff,
                      domain::Domain, state::State, precond::Term)
    push!(diff.del, precond.args[1])
    return diff
end
PRECOND_FUNCS[:not] = del_precond!

function and_precond!(diff::GenericDiff,
                      domain::Domain, state::State, precond::Term)
    for pre in get_args(precond)
        combine!(diff, precond_diff(domain, state, pre))
    end
    return diff
end
PRECOND_FUNCS[:and] = and_precond!

function forall_precond!(diff::GenericDiff,
                         domain::Domain, state::State, precond::Term)
    cond, pre = precond.args[1], precond.args[2]
    for s in satisfiers(domain, state, cond)
        combine!(diff, precond_diff(domain, state, substitute(pre, s)))
    end
    return diff
end
PRECOND_FUNCS[:forall] = forall_precond!


text:
module Finance
    export finance

    finance(n) = n*(n+1)*(n+2)÷2

end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
module BioinformaticsBISC195

export normalizeDNA,
       composition,
       gc_content,
       complement,
       reverse_complement,
       parse_fasta,
       isDNA,
       kmercount,
       kmerdistance,
       kmercollecting,
       kmercombining,
       monthlycomparison,
       shortendates,
       fasta_header,
       nwscore,
       nwsetupmatrix,
       nwscorematrix,
       swsetupmatrix,
       swscorematrix,
       nwalign,
       swalign
#=
include("assignment3code.jl")
include("assignment4code.jl")
include("assignment5code.jl")
include("assignment6code.jl")
=#

include("AlignmentAlgs.jl")

# # uncomment the following line if you intend to use BioSequences types
# using BioSequences

"""
    normalizeDNA(::AbstractString)

Ensures that a sequence only contains valid bases
(or `'N'` for unknown bases).
Returns a String.
"""
function normalizeDNA(seq)
    seq = uppercase(string(seq))
    seq = replace(seq, "Y" => "N")
    seq = replace(seq, "W" => "N")
    seq = replace(seq, "K" => "N")
    seq = replace(seq, "M" => "N")
    seq = replace(seq, "B" => "N")
    seq = replace(seq, "D" => "N")
    seq = replace(seq, "H" => "N")
    seq = replace(seq, "V" => "N")
    seq = replace(seq, "R" => "N")
    seq = replace(seq, "S" => "N")
    for base in seq
        # note: `N` indicates an unknown base
        occursin(base, "AGCTN") || error("invalid base $base")
    end

    return seq # change to `return LongDNASeq(seq)` if you want to try to use BioSequences types
end

"""
    composition(sequence)

    Counts the number of each type of base
    in a DNA sequence and returns a Dict of those counts for A, C, G, T, and N

    Ex.
    julia> composition("ACCGGGTTTTN")
   Dict{Char,Int64} with 5 entries:
     'A' => 1
     'G' => 3
     'T' => 4
     'N' => 1
     'C' => 2

   julia> composition("AAX")
   ERROR: Invalid base, X
"""
function composition(sequence)
    sequence = uppercase(sequence)
    bases = Dict('A' => 0, 'G' => 0, 'T' => 0, 'N' => 0, 'C' => 0)

    for base in sequence
        if base == 'A'
            bases['A'] = bases['A'] + 1
        elseif base == 'C'
            bases['C'] = bases['C'] + 1
        elseif base == 'G'
            bases['G'] = bases['G'] + 1
        elseif base == 'T'
            bases['T'] = bases['T'] + 1
        elseif base == 'N'
            bases['N'] = bases['N'] + 1
        else
            throw(ErrorException("Invalid base $base"))
        end
    end

    return bases
end

"""
    gc_content(seq=)

    Calculates the GC ratio of a DNA sequence.
    The GC ratio is the total number of G and C bases divided by the total length of the sequence.
    For more info about GC content, see here: https://en.wikipedia.org/wiki/GC-content
    Ambiguous bases (N) not included in calculations.

    Ex.

   julia> gc_content("ATNG")
   0.25

   julia> gc_content("ccccggggn")
   0.8888888888888888

"""
function gc_content(seq)
    seq = normalizeDNA(seq)
    if !isDNA(seq)
        throw(ErrorException("Invalid sequence $seq"))
    end
    seq = normalizeDNA(seq)

    #seqlength = length(seq)

    ## count the number of G's
    gs = count(==('G'), seq)
    ## count the number of C's
    cs = count(==('C'), seq)

    ts = count(==('T'), seq)
    as = count(==('A'), seq)

    return (gs + cs) / (gs + cs + ts + as)
end

"""
    complement(base)

Get the DNA complement of the provided base:

    A <-> T
    G <-> C

Accepts uppercase or lowercase `String` or `Char`,
but always returns an uppercase `Char`.
If a valid base is not provided, the function throws an error.
"""
function complement(base::Char)
    complements = Dict("A" => 'T',
                       "T" => 'A',
                       "G" => 'C',
                       "C" => 'G')

    base = uppercase(string(base))

    !(base in keys(complements)) && error("Invalid base $base")
    return complements[base]
end

"""
    complement(base)

    Get the DNA complement of the provided sequence

    Ex.
    julia> complement("ATTN")
    "TAAN"

    julia> complement("ATTAGC")
    "TAATCG"

    Accepts Strings only, returns string
    If a valid base is not provided, the function throws an error.
"""
function complement(seq::String)
    seq = uppercase(seq)
    returnSeq = ""
    for base in seq
        if base == 'N'
            returnSeq = returnSeq * base
        else
            returnSeq = returnSeq * complement(base)
        end
    end

    return returnSeq

end

"""
    reverse_complement(sequence)

Takes a DNA sequence and returns the reverse complement
of that sequence.

Takes lowercase or uppercase sequences,
but always returns uppercase.

Examples
≡≡≡≡≡≡≡≡≡≡
    julia> reverse_complement("AAATTT")
    "AAATTT"

    julia> reverse_complement("GCAT")
    "ATGC"

    julia> rc = reverse_complement("TTGGG");

    julia> println(rc)
    CCCAA

    julia> reverse_complement("ATTAGC")
    "GCTAAT"

    julia> reverse_complement("ATN")
    "NAT"
"""
function reverse_complement(sequence)
    reverseseq = complement(sequence)
    reverseseq = reverse(reverseseq)
    return reverseseq
end

"""
 function parse_fasta(path)

Reads a fasta-formated file and returns 2 vectors,
one containing headers,
the other containing the entire sequence as a `String`.

Note: function validates DNA sequences for correctness.

Example
≡≡≡≡≡≡≡≡≡
julia> ex1 = parse_fasta("data/ex1.fasta");

julia> ex1[1]
2-element Array{String,1}:
   "ex1.1 | easy"
   "ex1.2 | multiline"

julia> ex1[2]
2-element Array{String,1}:
   "AATTATAGC"
   "CGCCCCCCAGTCGGATT"

julia> ex2 = parse_fasta("data/ex2.fasta");
ERROR: invalid base H
"""
function parse_fasta(path)
 headersVect = Vector{String}()
 seqVect = Vector{String}()
 toBeCombined = Vector{String}()

 for line in eachline(path)
     #@info "line is " line
     #@info "to be combined is" toBeCombined
     if startswith(line, '>')
         push!(headersVect, line[2:end])
         if size(toBeCombined,1) > 0
             push!(seqVect, join(toBeCombined))
             toBeCombined = Vector()
         end
     elseif !isDNA(normalizeDNA(line))
        throw(ErrorException("Invalid base $base"))
     elseif isDNA(normalizeDNA(line))
         push!(toBeCombined, line)

     end
     #@info "to be combined is" toBeCombined
 end
 push!(seqVect, join(toBeCombined))

 return (headersVect, seqVect)
end

"""
    isDNA(sequence)

Examples
========

    isDNA("atGG")
    true

    isDNA("ATYYG")
    false

Checks the given sequence, and if one of its letters is not a base as stored in the bases array, then returns false.
otherwise, returns true.
NOTE: does not make use of normalizeDNA to correct the DNA into a uniform format with accepted bases; that will need to be done before isDNA is used
"""
function isDNA(sequence)
    bases = ['A', 'C', 'G', 'T', 'N']
    sequence = uppercase(sequence)
    for letter in sequence
        if !(letter ∈ bases)
            return false
        end
    end
    return true
end

"""
   kmercount(sequence, k)

Finds all kmers in a sequence,
returning a dictionary of those kmers
and the number of times they appear in the sequence.

Examples
≡≡≡≡≡≡≡≡≡≡

    julia> kmercount("ggg", 3)
    Dict{Any,Any} with 1 entry:
    "GGG" => 1

    julia> kmercount("ATATATATA", 4)
    Dict{Any,Any} with 2 entries:
    "TATA" => 3
    "ATAT" => 3

    julia> kmercount("ATATATATAx", 4)
    ERROR: Invalid base X encountered

    julia> kmercount("A", 2)
    ERROR: k must be a positive integer less than the length of the sequence
"""
function  kmercount(sequence, k)
    1 <= k <= length(sequence) || error("k must be a positive integer less than the length of the sequence")
    kmers = Dict() # initialize dictionary

    stopindex = length(sequence)-k+1

    for i in 1:stopindex
        kmer = sequence[i:i+k-1] # Change to index the sequence from i to i+k-1
        kmer = normalizeDNA(kmer)
        if !occursin("N", kmer)
            if haskey(kmers, kmer)
            #       add 1 to the value referenced by that kmer
                kmers[kmer] = kmers[kmer] + 1
            #   otherwise
            else
            #       make a new entry in the dictionary with a value of 1
                kmers[kmer] = 1
                #setindex!(kmers, 1, kmer)
            end
        end
        #println(kmers)
    end
    return kmers
end

"""
    kmerdistance(set1, set2)

Takes two kmer sets and calculates the distance between them, using distance metric that is 1 - (length of intersection / length of union).
Or alternatively, ((length of set1 - set2) + (length of set2 - set1)) / length of union.
Note: does not check for validity of kmers

Examples
========

kmerdistance(["AAA", "AGT"], ["AAA", "GGG"])
0.666667 (2/3)

kmerdistance(["AAA", "AGT"], ["TTT", "GGG"])
1

kmerdistance(["AAA", "ATT"], ["AGT", "AAA", "GTGG"])
0.75

"""
function kmerdistance(set1, set2)
    intersectLength = length(intersect(set1, set2))
    unionLength = length(union(set1, set2))
    return 1 - (intersectLength / unionLength)
end

"""
    kmercollecting(set, k)

    Takes in a set of sequences and k, then returns a collection of vectors containing all kmers of said length k in the entire set.
    For example, it can take all sequences that were parsed from a certain file, and it returns all kmers of length k that exist in that file.

    Examples
    ========
    ex1 = parse_fasta("C:/Users/Jen/Documents/my_repo/BioinformaticsBISC195/data/ex1.fasta")
    (Any["ex1.1 | easy", "ex1.2 | multiline"], Any["AATTATAGC", "CGCCCCCCAGTCGGATT"])

    kmercollecting(ex1[2], 3)
    2-element Vector{Any}:
    Any["AAT", "ATT", "TAT", "ATA", "TAG", "TTA", "AGC"]
    Any["AGT", "CCC", "CGG", "ATT", "GAT", "GCC", "CAG", "GTC", "CGC", "GGA", "CCA", "TCG"]

    set1 = ("ACCGGTT", "AAA")
    kmercollecting(set1, 3)
    2-element Vector{Any}
    Any["ACC", "CCG", "CGG", "GGT", "GTT"]
    Any["AAA"]

"""
function kmercollecting(set, k)
    map(set) do seq
        Set(seq[i:(i+k-1)] for i in 1:(length(seq)-k+1))
    end
end

#check that the sequences are organized by date
"""
    monthlycomparison(original, allSeq, totalMonths)

    Takes the string of original sequence that all others will be compared against as well as a collection of all the entire parsed fasta file and the total months to be as input.
    Note: at the moment only works if the original string of data is from Dec 2019
    With more time I would modify this to be much more flexible, ie you would be able to put in you own start and end date rather than hard coding the original dates like I did here.
    Note 1: Also depends on the fact that there is at least one datapoint for each month, which there is in the big file of for all sequences in Asia, but is not necessarily always true
    Note 2: the swscorematrix format is very time consuming so for my purposes the totalMonths must be fairly low or the length of the sequences very short in order to reduce comparisons

    Examples
    ========
    covgen_asiashort.fasta is a version of the covgen_asia file with only 4 data points listed, from month 2019-21 thru 2020-03
    headers1 = ["Should Work | 2019-12", "Should work|2020-01-01| extra", "should not work | 2020", "should work|2020-02", "should work|2020-03-2"]
    seq1 = ["ACCTTGGA", "ACCCTGGA", "ACCCTGG", "ACCCTGC", "AACCTGC"]
    data1 = [headers1, seq1]

    monthlycomparison(data1[2][1], data1, 3)
    3-element Vector{Int64}:
    6
    5
    4

    monthlycomparison(data1[2][1], data1, 5)
    ERROR: Date 2020-04 does not exist in file; cannot make complete graph.

    monthlycomparison(data1[2][1], data1, 5)
    ERROR: invalid base I

"""
function monthlycomparison(original, allSeq, totalMonths)
	# stores all of dates for each sequences in its own file for easier traversal later
    if isDNA(normalizeDNA(original)) == false
        throw(ErrorException("Original sequences is not valid DNA"))
    end
    dates = shortendates(allSeq[1])
    # I originally tried to put all of the shortendates code inside monthly comparions, but I ended up making them two separate functions for easier testing

	# filling out dateIndices to hold ints that represent the indices of the next correct sequences to grab
    # will be used to grab the correct sequences from allSeq to that will be used to score the differences per month
	dateIndices = Vector{Int}()
	# targetMonth and Year are hard coded for the original sequence having a date of Dec 2019
    targetMonth = 1
	targetYear = 2020
	for i in 1:totalMonths
		targetString = string(targetYear, "-", lpad(targetMonth,2,"0"))
		# @info targetString
		# @info findfirst(targetString .== dates)
        if findfirst(targetString .== dates) === nothing
            throw(ErrorException("Date $targetString does not exist in file; cannot make complete graph."))
        end
		push!(dateIndices, findfirst(targetString .== dates))
		if targetMonth == 12
			targetMonth = 1
			targetYear = targetYear+1
		else
			targetMonth = targetMonth+1
		end
	end

	score = Vector{Int}()
	for index in dateIndices
		push!(score, maximum(swscorematrix(original, allSeq[2][index])))
	end
	return score
end

"""
    shortendates(dates)

    Takes in a collection of parsed header lines and returns a collection of just the dates of these headers, in the same order, using only year-month format.
    Does not save or return any dates that only have a year listed.
    Designed as a helper function for monthlycomparison
    Note: strongly depends on the format of the headers in the fasta file, namely that the date is the second item the fasta header and the date is in the order year-month-day.

    Examples
    ========
    dates1 = ["Should Work | 2031-04", "Should work|2000-12-01| extra", "should not work | 2020", "should work|2001-01"]
    shortendates(dates1)
    3-element Vector{String}:
    ["2031-04", "2000-12", "2001-01"]

    dates2 = ["should not work| XX-XX", "should work|XXXX-XX", "should work   |  XXXXXXXextra"]
    shortendates(dates2)
    2-element Vector{String}:
    "XXXX-XX"
    "XXXXXXX"
"""
function shortendates(parsed_headers)
    shortDates = Vector{String}()
	for headers in parsed_headers
		header_components = fasta_header(headers)
		if length(header_components[2]) >= 7 #if its shorter than this then it does not have a month listed, so I can't use it and don't need it in the vector
			push!(shortDates, header_components[2][1:7]) #1-7 gives only the year and month, which is all that is needed to search for the right sequence
		end
	end

    return shortDates
end

"""
 fasta_header(header)

Divides a fasta header into its component parts,
removing any leading or trailing spaces.
NOTE: does not check for or depend on having a '>' to mark as a header

Example
≡≡≡≡≡≡≡≡≡
 julia> fasta_header("M0002 |China|Homo sapiens")
 ("M0002", "China", "Homo sapiens")

 julia> fasta_header("AAATTC")
 Error: Invalid header (headers must start with '>')

 julia> fasta_header("Another sequence")
 ("Another sequence",)

 julia> fasta_header("headers| can | have| any number | of | info blocks")
 ("headers", "can", "have", "any number", "of", "info blocks")
"""
function fasta_header(header)
   # I ended up not using  returnTuple = (strip(component) for component in splitVect) bc it returned a Base.generator that was hard to deal with,
   # but I don't think this way constantly creates and destroys a tuple over and over again either
   splitVect = split(header, "|")
   return Tuple(map(strip, splitVect))
end

end # module BioinformaticsBISC195


text:
# These tests are closely based on the server tests in HTTP.jl:
# https://github.com/JuliaWeb/HTTP.jl/blob/master/test/server.jl

using HTTP, Joseki, JSON, Sockets
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

@testset "Joseki" begin

    # Simple server
    function pow(req::HTTP.Request)
        j = HTTP.queryparams(HTTP.URI(req.target))
        if !(haskey(j, "x")&haskey(j, "y"))
            return error_responder(req, "You need to specify values for x and y!")
        end
        x = parse(Float32, j["x"])
        y = parse(Float32, j["y"])
        json_responder(req, x^y)
    end

    endpoints = [
        (pow, "GET", "/pow")
    ]
    server = Joseki.server(endpoints) # Default middleware

    @info "Starting test server"
    server_task = @async HTTP.serve(server, ip"127.0.0.1", 8000; verbose=false)
    sleep(1.0)

    res = HTTP.get("http://localhost:8000/pow/?x=2&y=3")
    @test res.status == 200
    @test HTTP.hasheader(res, "Access-Control-Allow-Methods", "GET,PUT,POST,DELETE,PATCH,OPTIONS")

    @info "Shutting down test server"

    put!(server.in, HTTP.Servers.KILL)
    sleep(2)
    @test istaskdone(server_task)
end


text:
export
    instant_runtime_from_instances,
    retrieve_values_from_instant_runtime!,
    initial_instant_runtime,
    instant_name,
    dynamic_name_and_time,
    instant_node

"""
Create a name in an instant network corresponding to the given dynamic name and time.
"""
function instant_name(dynamic_name::Symbol, time::Number)::Symbol
    return Symbol(dynamic_name, "_", time)
end

"""
Create a dynamic name and time from an instant node. T is the time type.
"""
function dynamic_name_and_time(instant_node::Node, T = Int)::Tuple{Symbol, T}
    instant_name = collect(repr(instant_node.name))
    # Need to handle two different representations of symbols
    if length(instant_name) >= 6 && instant_name[1:6] == ['S', 'y', 'm', 'b', 'o', 'l']
        instant_name = instant_name[9:length(instant_name)-2] # strip parens and quotes
    else
        instant_name = instant_name[2:length(instant_name)] # strip leading colon
    end
    len = length(instant_name)
    i = findlast(x -> x == '_', instant_name)
    dynamic_name = Symbol(String(instant_name[1:i-1]))
    time = parse(T, String(instant_name[i+1:len]))
    return (dynamic_name, time)
end

"""
Create an instant node from a dynamic variable instance.
"""
function instant_node(dyninst::VariableInstance)
    name = instant_name(get_name(dyninst), get_time(dyninst))
    model = SimpleModel(get_sfunc(dyninst))
    var = Variable(name, model)
    return var
end

"""
Create an instant node from a dynamic placeholder instance.
"""
function instant_node(dyninst::PlaceholderInstance{O}) where O
    name = instant_name(get_name(dyninst), get_time(dyninst))
    ph = Placeholder{O}(name)
    return ph
end

"""
    instant_runtime_from_instances(runtime::DynamicRuntime, instances::Vector{Instance})

Create an instant runtime from the given instances in the given dynamic runtime.

This runtime has an instant network that contains a variable for each instance in `insts`,
tagged with the time of the instance.
The network also contains a placeholder for each instance in `placeholder_insts`.
The function also instantiates the variables in the instant runtime and stores any runtime
values from the dynamic runtime with the corresponding instances in the instant runtime.
This function is useful for running instant algorithms on a time window
for dynamic reasoning.
"""
function instant_runtime_from_instances(dynrun::DynamicRuntime, dyninsts::Vector{Instance})
    dynnet = get_network(dynrun)
    forward_index = Dict{Instance, Node}()
    back_index = Dict{Node, Instance}()
    placeholders = Placeholder[]
    variables = Variable[]
    nodes = Node[]

    for dyninst in dyninsts
        node = instant_node(dyninst)
        if get_node(dyninst) isa Variable
            push!(variables, node)
        else
            push!(placeholders, node)
        end
        forward_index[dyninst] = node
        back_index[node] = dyninst
        push!(nodes, node)
    end

    instgraph = VariableGraph()
    for node in nodes
        nodepars = Node[]
        dyninst = back_index[node]
        insttime = get_time(dyninst)
        dynnode = get_node(dyninst)
        dynpars = get_transition_parents(dynnet, dynnode)
        for dynpar in dynpars
            # Find the most recent parent instance equal or before this variable's instance
            parinst = latest_instance_before(dynrun, dynpar, insttime, dynpar != dynnode)
            if isnothing(parinst)
                error("Variable does not have parent in instances")
            elseif !(parinst in dyninsts)
                # should be a placeholder
                ph = Placeholder{output_type(dynpar)}(get_name(dynpar))
                parinst = PlaceholderInstance(ph, get_time(parinst))
            end
            push!(nodepars, forward_index[parinst])
        end
        instgraph[node] = nodepars
    end

    instnet = InstantNetwork(variables, instgraph, placeholders)
    instrun = Runtime(instnet)
    ensure_all!(instrun)

    for ((dyninst, valuename), value) in dynrun.values
        if dyninst in dyninsts
            instinst = current_instance(instrun, forward_index[dyninst])
            set_value!(instrun, instinst, valuename, value)
        end
    end

    for (k,v) in get_state(dynrun)
        set_state!(instrun, k, v)
    end

    return instrun
end

"""
    Creates an instant runtime for the first time step.
"""
function initial_instant_runtime(dynrun::DynamicRuntime)
    variables = Variable[]
    placeholders = Placeholder[]
    dynnet = get_network(dynrun)
    node_index = Dict{Symbol, Node}() # Maps names of dynamic nodes to instant instances
    instgraph = VariableGraph()
    for dynnode in topsort(get_initial_graph(dynnet))
        nodename = get_name(dynnode)
        dyninst = current_instance(dynrun, dynnode)
        instvar = instant_node(dyninst)
        if instvar isa Variable
            push!(variables, instvar)
        else
            push!(placeholders, instvar)
        end
        node_index[nodename] = instvar
        dynpars = get_initial_parents(dynnet, dynnode)
        instpars = Node[]
        for dynpar in dynpars
            instpar = node_index[get_name(dynpar)]
            push!(instpars, instpar)
        end
        instgraph[instvar] = instpars
    end
    instnet = InstantNetwork(variables, instgraph, placeholders)
    instrun = Runtime(instnet)
    for node in values(node_index)
        instantiate!(instrun, node, 0)
    end
    return instrun
end

"""
    retrieve_values_from_instant_runtime!(dynrun::DynamicRuntime, instrun::InstantRuntime)

    Retrieve values in a dynamic runtime from an instant runtime constructed
    using `instant_runtime_from_instances`.
"""
function retrieve_values_from_instant_runtime!(dynrun::DynamicRuntime{T},
        instrun::InstantRuntime) where T
    index = Dict{Instance, Instance}()
    for node in get_nodes(get_network(instrun))
        # If instrun has been constructed correctly from dynrun,
        # there should be no runtime errors in this code.
        instinst = current_instance(instrun, node)
        (dynname, dyntime) = dynamic_name_and_time(node, T)
        dynnode = get_node(get_network(dynrun), dynname)
        dyninst = get_instance(dynrun, dynnode, dyntime)
        index[instinst] = dyninst
    end

    # TODO: This assumes that the values in the dynamic runtime take the same value as in the instant runtime.
    # This is not necessarily the case, e.g. for :particles, which contains samples that are dictionaries from node name to
    # node value. The node names in the instant runtime and dynamic runtime are different.
    # For this, importance stores beliefs with nodes in the instant runtime that can get translated.
    # A better general method is needed.
    for ((instinst, valuename), value) in instrun.values
        set_value!(dynrun, index[instinst], valuename, value)
    end
    for (k, v) in get_state(instrun)
        set_state!(dynrun, k, v)
    end
end


text:
module trackme
    expression = :(:x + :y -1)
end

max_stars_repo_name: hustf/EquationEditor

max_issues_repo_name: hustf/EquationEditor

max_forks_repo_name: hustf/EquationEditor


text:
module Pipe
#Reflow piped things replacing the _ in the next section

export @pipe

function rewrite(ff::Expr,target)
    function replace(arg::Any)
        arg #Normally do nothing
    end
    function replace(arg::Symbol)
        if arg==:_
            target
        else
            arg
        end
    end
    function replace(arg::Expr)
        rep = copy(arg)
        rep.args = map(replace,rep.args)
        rep
    end

    if (ff.head==:call)
        rep_args = map(replace,ff.args)
        if ff.args != rep_args
            #_ subsitution
            ff.args=rep_args
            return ff
        end
    end
    #No subsitution was done (either cos not a call, or cost no _ found)
    #Apply to a function that is being returned by ff, (ff could be a function call or something more complex)
    rewrite_apply(ff,target)
end

function rewrite_apply(ff, target)
    #function application
    :($ff($target))
end

function rewrite(ff::Symbol, target)
    rewrite_apply(ff,target)
end

function funnel(ee::Any) #Could be a Symbol could be a literal
    #first (left most) input
    ee
end

function funnel(ee::Expr)
    if (ee.args[1]==:|>)
        ff = ee.args[3]
        target = funnel(ee.args[2]) #Recurse

        rewrite(ff,target)
    else
        #Not in a piping situtation
        ee #make no change
    end
end

macro pipe(ee)
    esc(funnel(ee))
end

end


text:
using Base58
using Base.Test

@testset "Base58" begin

    test_data = hcat(
        [b"",                b""],
        [[0x00],             b"1"],
        [[0x00, 0x00],       b"11"],
        [b"hello world",     b"StV1DL6CwTryKyV"],
        [b"\0\0hello world", b"11StV1DL6CwTryKyV"],
        [nothing,            b"3vQOB7B6uFg4oH"],
        [b""" !"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\\]^_`abcdefghijklmnopqrstuvwxyz{|}~""",
                             b"3WSNuyEGf19K7EdeCmokbtTAXJwJUdvg8QXxAacYC7kR1bQoYeWVr5iMUHvxvv4FCFY48dVUrX6BrFLod6zsEhHU9NciUXFM17h1qtveYD7ocbnXQyuY84An9nAvEjdt6H"]
    )

    # The weird eval magic gives the value of i if the test fails
    for i in 1:size(test_data, 2)
        if(test_data[1, i] != nothing)
            @eval @test Base58.base58encode($(test_data)[1, $i]) == $(test_data)[2, $i]
        end
    end

    for i in 1:size(test_data, 2)
        if(test_data[1, i] == nothing)
            @eval @test_throws ArgumentError Base58.base58decode($(test_data)[2, $i])
        else
            @eval @test Base58.base58decode($(test_data)[2, $i]) == $(test_data)[1, $i]
        end
    end

    for i in 1:size(test_data, 2)
        if(test_data[1, i] == nothing)
            @eval @test_throws ArgumentError Base58.base58decode([b"   "..., $(test_data)[2, $i]...])
        else
            @eval @test Base58.base58decode([b"    "..., $(test_data)[2, $i]...]) == $(test_data)[1, $i]
        end
    end

    @test base58decode(b"     ") == UInt8[]
    @test base58decode(b"11111") == [0x00 for i in 1:5]
end

@testset "Base58Check" begin

    # https://en.bitcoin.it/wiki/Technical_background_of_version_1_Bitcoin_addresses
    address = [0x00, 0x01, 0x09, 0x66, 0x77, 0x60, 0x06, 0x95, 0x3D, 0x55, 0x67, 0x43,
               0x9E, 0x5E, 0x39, 0xF8, 0x6A, 0x0D, 0x27, 0x3B, 0xEE]
    address_encoded = b"16UwLL9Risc3QfPqBUvKofHmBQ7wMtjvM"
    wrong_address_encoded = b"16UwLL9Risc3cfPqBUvKofHmBQ7wMtjvM"

    @test Base58.base58checkencode(address) == address_encoded
    @test Base58.base58checkdecode(address_encoded, true) == address
    @test Base58.base58checkdecode(address_encoded, false) == address
    @test_throws ArgumentError Base58.base58checkdecode(wrong_address_encoded, true)

end


text:
using CUDNN: cudnnConvolutionBackwardBias, cudnnConvolutionBackwardFilter, cudnnConvolutionBackwardData, cudnnConvolutionForward
using CUDNN: ConvolutionDescriptor, CUDNN_CROSS_CORRELATION, CUDNN_CONVOLUTION
using CUDArt

corrd = ConvolutionDescriptor(padding=(1,1), mode=CUDNN_CROSS_CORRELATION)
convd = ConvolutionDescriptor(mode=CUDNN_CONVOLUTION)

x = reshape(Float64[1:9],3,3,1,1,1); tx = CudaArray(x)
w = reshape(Float64[1:4],2,2,1,1,1); tw = CudaArray(w)
@show size(x) #3x3
@show size(w) #2x2
println("x:")
println(x)
println("w:")
println(w)
println("corr y:")
println(to_host(cudnnConvolutionForward(tx, tw;convDesc=corrd)))
println("conv y:")
println(to_host(cudnnConvolutionForward(tx, tw;convDesc=convd)))
println()

#dy = rand(size(ty)); tdy = CudaArray(dy) #2x2
dy = reshape(Float64[5 6 7 8],2,2,1,1); tdy = CudaArray(dy) #2x2
@show size(dy)
println("dy:")
println(dy)
dx = zeros(x); tdx = CudaArray(dx) #3x3
dw = zeros(w); tdw = CudaArray(dw) #2x2

println("conv dw:")
cudnnConvolutionBackwardFilter(tx, tdy, tdw; convDesc=convd)
println(to_host(tdw))
println("corr dw:")
cudnnConvolutionBackwardFilter(tx, tdy, tdw; convDesc=corrd)
println(to_host(tdw))
println("conv dx:")
cudnnConvolutionBackwardData(tw, tdy, tdx; convDesc=convd)
println(to_host(tdx))
println("corr dx:")
cudnnConvolutionBackwardData(tw, tdy, tdx; convDesc=corrd)
println(to_host(tdx))

println("testing pooling..")
using CUDNN: cudnnPoolingForward, cudnnPoolingBackward, cudnnGetPoolingNdForwardOutputDim, PoolingDescriptor, CUDNN_POOLING_MAX
x = reshape(Float64[5 2 3 4 5 6 7 7 3 1 4 5], 3, 4, 1, 1); tx = CudaArray(x)
println("x:"); println(x)
dims = (2,2)
@show dims
pd1 = PoolingDescriptor(dims; padding=map(x->0,dims), stride=dims, mode=CUDNN_POOLING_MAX)
ydims = cudnnGetPoolingNdForwardOutputDim(pd1, tx)
y = zeros(ydims); ty = CudaArray(y)
cudnnPoolingForward(pd1, tx, ty)
println(to_host(ty))
dy = reshape(Float64[5 3 2 6],2,2,1,1); tdy = CudaArray(dy)
dx = zeros(x); tdx = CudaArray(dx)
cudnnPoolingBackward(pd1,ty,tdy,tx,tdx)
println(to_host(tdx))
println(to_host(ty))


text:
abstract type Helmholtz3DOp <: MaxwellOperator3D end
abstract type Helmholtz3DOpReg <: MaxwellOperator3DReg end
"""
```
∫_Γ dx ∫_Γ dy \\left(α G g(x) n_x ⋅ n_y f(y) + β G \\mbox{curl} g(x) ⋅ \\mbox{curl} f(y) \\right)
```
with ``G(x,y) = \\frac{e^{-γ |x-y|}}{4 π |x-y|}``
"""
struct HH3DHyperSingularFDBIO{T,K} <: Helmholtz3DOp
    "coefficient of the weakly singular term"
    alpha::T
    "coefficient of the hyper singular term"
    beta::T
    "`im*κ` with `κ` the wave number"
    gamma::K
end

HH3DHyperSingularFDBIO(gamma) = HH3DHyperSingularFDBIO(gamma^2, one(gamma), gamma)
scalartype(op::HH3DHyperSingularFDBIO) = promote_type(typeof(op.alpha), typeof(op.beta), typeof(op.gamma))

"""
```math
a(u,v) = α ∬_{Γ×Γ} u(x) G_{γ}(|x-y|) v(y)
```
with ``G_{γ}(r) = \\frac{e^{-γr}}{4πr}``.
"""
struct HH3DSingleLayerFDBIO{T,K} <: Helmholtz3DOp
    alpha::T
    gamma::K
end

struct HH3DSingleLayerReg{T,K} <: Helmholtz3DOpReg
    alpha::T
    gamma::K
end

struct HH3DSingleLayerSng{T,K} <: Helmholtz3DOp
    alpha::T
    gamma::K
end

struct HH3DDoubleLayer{T,K} <: Helmholtz3DOp
    alpha::T
    gamma::K
end

struct HH3DDoubleLayerTransposed{T,K} <: Helmholtz3DOp
    alpha::T
    gamma::K
end

function quaddata(op::Helmholtz3DOp, test_refspace::LagrangeRefSpace,
        trial_refspace::LagrangeRefSpace, test_elements, trial_elements)

    test_eval(x)  = test_refspace(x,  Val{:withcurl})
    trial_eval(x) = trial_refspace(x, Val{:withcurl})

    # The combinations of rules (6,7) and (5,7 are) BAAAADDDD
    # they result in many near singularity evaluations with any
    # resemblence of accuracy going down the drain! Simply don't!
    # (same for (5,7) btw...).
    # test_qp = quadpoints(test_eval,  test_elements,  (6,))
    # bssi_qp = quadpoints(trial_eval, trial_elements, (7,))

    test_qp = quadpoints(test_eval,  test_elements,  (4,))
    bsis_qp = quadpoints(trial_eval, trial_elements, (7,))

    return test_qp, bsis_qp
end

function quaddata(op::Helmholtz3DOp, test_refspace::subReferenceSpace,
        trial_refspace::subReferenceSpace, test_elements, trial_elements)

    test_qp = quadpoints(test_refspace,  test_elements,  (4,))
    bsis_qp = quadpoints(trial_refspace, trial_elements, (7,))

    return test_qp, bsis_qp
end

function quadrule(op::HH3DSingleLayerFDBIO,
        test_refspace::LagrangeRefSpace{T,0} where T,
        trial_refspace::LagrangeRefSpace{T,0} where T,
        i, test_element, j, trial_element, quadrature_data)

    tol, hits = sqrt(eps(eltype(eltype(test_element.vertices)))), 0
    for t in test_element.vertices
        for s in trial_element.vertices
            norm(t-s) < tol && (hits +=1; break)
    end end

    test_quadpoints  = quadrature_data[1]
    trial_quadpoints = quadrature_data[2]

    hits != 0 && return WiltonSEStrategy(
        test_quadpoints[1,i],
        DoubleQuadStrategy(
            test_quadpoints[1,i],
            trial_quadpoints[1,j]))

    return DoubleQuadStrategy(
        quadrature_data[1][1,i],
        quadrature_data[2][1,j])
end

function quadrule(op::HH3DSingleLayerFDBIO, test_refspace::subReferenceSpace,
    trial_refspace::subReferenceSpace, i, test_element, j, trial_element, quadrature_data)

    # tol, hits = 1e-10, 0
    # for t in getelementVertices(test_element)
    #     for s in getelementVertices(trial_element)
    #         norm(t-s) < tol && (hits +=1; break)
    # end end
    #
    # test_quadpoints  = quadrature_data[1]
    # trial_quadpoints = quadrature_data[2]

    # hits != 0 && return WiltonSEStrategy(
    #     test_quadpoints[1,i],
    #     DoubleQuadStrategy(
    #         test_quadpoints[1,i],
    #         trial_quadpoints[1,j]))

    return DoubleQuadStrategy(
        quadrature_data[1][1,i],
        quadrature_data[2][1,j])
    # return SauterSchwabStrategy(hits)
     # test_qd = qd[1]
     # trial_qd = qd[2]
     #
     # DoubleQuadStrategy(
     #    test_qd[1,i], # rule 1 on test element j
     #    trial_qd[1,j] # rule 1 on trial element i
     # )
end

function quadrule(op::Helmholtz3DOp,
    test_refspace::RefSpace, trial_refspace::RefSpace,
    i, test_element, j, trial_element, quadrature_data)

    test_quadpoints  = quadrature_data[1]
    trial_quadpoints = quadrature_data[2]

    return DoubleQuadStrategy(
        quadrature_data[1][1,i],
        quadrature_data[2][1,j])
end

function quadrule(op::Helmholtz3DOp,
    test_refspace::RTRefSpace, trial_refspace::RTRefSpace,
    i, test_element, j, trial_element, quadrature_data)

    test_quadpoints  = quadrature_data[1]
    trial_quadpoints = quadrature_data[2]

    return DoubleQuadStrategy(
        quadrature_data[1][1,i],
        quadrature_data[2][1,j])
end

function integrand(op::HH3DHyperSingularFDBIO,
        kernel, test_values, test_element, trial_values, trial_element)

    α = op.alpha
    β = op.beta

    G = kernel.green

    g, curlg = test_values
    f, curlf = trial_values

    nx = normal(test_element)
    ny = normal(trial_element)

    α*dot(nx,ny)*g*f*G + β*dot(curlg,curlf)*G
end

HH3DSingleLayerFDBIO(gamma) = HH3DSingleLayerFDBIO(one(gamma), gamma)

regularpart(op::HH3DSingleLayerFDBIO) = HH3DSingleLayerReg(op.alpha, op.gamma)
singularpart(op::HH3DSingleLayerFDBIO) = HH3DSingleLayerSng(op.alpha, op.gamma)

function integrand(op::Union{HH3DSingleLayerFDBIO,HH3DSingleLayerReg},
        kernel, test_values, test_element, trial_values, trial_element)

    α = op.alpha
    G = kernel.green

    g = test_values.value
    f = trial_values.value

    α*dot(g, G*f)
end

function innerintegrals!(op::HH3DSingleLayerSng, test_neighborhood,
        test_refspace::LagrangeRefSpace{T,0} where {T},
        trial_refspace::LagrangeRefSpace{T,0} where {T},
        test_elements, trial_element, zlocal, quadrature_rule::WiltonSEStrategy, dx)

    γ = op.gamma
    α = op.alpha

    s1, s2, s3 = trial_element.vertices

    x = cartesian(test_neighborhood)
    n = normalize((s1-s3)×(s2-s3))
    ρ = x - dot(x - s1, n) * n

    scal, vec = WiltonInts84.wiltonints(s1, s2, s3, x, Val{1})
    ∫G = (scal[2] - γ*scal[3] + 0.5*γ^2*scal[4]) / (4π)

    zlocal[1,1] += α * ∫G * dx
    return nothing
end

function integrand(biop::HH3DDoubleLayer,
        kernel, fp, mp, fq, mq)

    nq = normal(mq)
    fp[1] * dot(nq, -kernel.gradgreen) * fq[1]
end

function integrand(biop::HH3DDoubleLayerTransposed,
        kernel, fp, mp, fq, mq)

    np = normal(mp)
    fp[1] * dot(np, kernel.gradgreen) * fq[1]
end

module Helmholtz3D

    using ..BEAST
    const Mod = BEAST

    function singlelayer(;
        gamma=nothing,
        wavenumber=nothing,
        alpha=nothing)

        if (gamma == nothing) && (wavenumber == nothing)
            error("Supply one of (not both) gamma or wavenumber")
        end

        if (gamma != nothing) && (wavenumber != nothing)
            error("Supply one of (not both) gamma or wavenumber")
        end

        if gamma == nothing
            if iszero(real(wavenumber))
                gamma = -imag(wavenumber)
            else
                gamma = im*wavenumber
            end
        end

        @assert gamma != nothing
        alpha == nothing && (alpha = one(real(typeof(gamma))))

        Mod.HH3DSingleLayerFDBIO(alpha,gamma)
    end

    hypersingular(;
            gamma=error("propagation constant is a required argument"),
            alpha=gamma^2,
            beta=one(gamma)) =
        Mod.HH3DHyperSingularFDBIO(alpha, beta, gamma)

    planewave(;
            direction=error("direction is a required argument"),
            wavenumber=error("wavenumber is a required arguement"),
            amplitude=one(eltype(direction))) =
        Mod.HH3DPlaneWave(direction, wavenumber)

    doublelayer(;gamma=error("gamma missing"), alpha=one(gamma)) =
        Mod.HH3DDoubleLayer(alpha, gamma)

    doublelayer_transposed(;gamma=error("gamma missing"), alpha=one(gamma)) =
        Mod.HH3DDoubleLayerTransposed(alpha, gamma)
end

export Helmholtz3D


text:
mutable struct LaserFeatures
  t::Float64
  feats::Array{Float64,2}
end


text:
import java.util.ArrayList;

import accrue.cryptoerase.runtime.Condition;

class C {
    static final Condition c = new Condition();
    void m() {
	ArrayList foo = new ArrayList();
	Object{L /c H} o = new Object();
	foo.add(o);
	Object{L} out = foo.get(0);
    }
    public static void main(String[] args) {
	new C().m();
    }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
### Rules
using Test
using Cassette
import Cassette: tag, untag, Tagged, metadata, hasmetadata, istagged
using SparseArrays
using SparsityDetection
import SparsityDetection: abstract_run, HessianSparsityContext,
                          JacobianSparsityContext, TermCombination, HessInput
using Test

Term(i...) = TermCombination(Set([Dict(j=>1 for j in i)]))

function jactester(f, Y, X, args...)
    ctx, val = jacobian_sparsity(f, Y, X, args...; raw=true)
end

jactestmeta(args...) = jactester(args...)[1].metadata
jactestval(args...) = jactester(args...) |> ((ctx,val),) -> untag(val, ctx)
jactesttag(args...) = jactester(args...) |> ((ctx,val),) -> metadata(val, ctx)

function hesstester(f, X, args...)
    ctx, terms, val = hessian_sparsity(f, X, args...; raw=true)
end

hesstestmeta(args...) = hesstester(args...)[1].metadata
hesstestval(args...)  = hesstester(args...) |> ((ctx,terms,val),) -> untag(val, ctx)
hesstesttag(args...)  = hesstester(args...) |> ((ctx,terms,val),) -> metadata(val, ctx)
hesstestterms(args...)  = hesstester(args...) |> ((ctx,terms,val),) -> terms

Base.show(io::IO, ::Type{<:Cassette.Context}) = print(io, "ctx")
Base.show(io::IO, ::Type{<:Tagged{<:Any, A}}) where {A} = print(io, "tagged{", A, "}")


text:
using SafeTestsets

@safetestset "All Test Sets" begin
    include("test_sets.jl")
end # @safetestset "All Test Sets"


text:
include("graph_utils.jl")
include("uniform_diforest.jl")
include("tree_preference_diforest.jl")
include("prufer_codes.jl")