text:
#=
A BoundingBox implementation.

A BoundingBox is a volume
that surrounds all the elements of an array
where the elements have a certain property.

Axis aligned.  That is, the coordinates of the bounding box
are in the frame of the array.

Most often, for Bool arrays, where the property is isTrue()
But also for a trivial property "true"
i.e. the bounding box contains all the elements of the array.

A volume is defined by two points i.e. CartesianIndex.
See "Geometric Approximation Algorithms" by S. Har-Peled, 2008.
At least up to 3D it is computable in O(n3)
For 4D, two points define the <diameter>,
the length of the maximum projection onto a line.
Still not sure the below is correct for 4D.

Multi-dimensional.
=#

struct BoundingBox{DimensionCount}
    least::CartesianIndex{DimensionCount}
    most::CartesianIndex{DimensionCount}
end

#=
Outer constructor
From a simple array.
Property: trivial function that always returns true.
I.E. Returns the bounding box of the entire array.
=#
function BoundingBox(
        tensor) # ::AbstractArray{ValueType, DimensionCount})
    # where{ValueType, DimensionCount}
    #=

    =#
    # All the indices, sorted
    indices = CartesianIndices(tensor)
    # Assert they are in order !!!

    # More or less: upper left point and lower right point
    result = BoundingBox(indices[1], indices[end])

    @debug "BoundingBox" result
    return result
end

#=
Outer constructor
From a MaskedImage.
Returning the bounding box of trues.
=#
#tensor::MaskedImage{ValueType, DimensionCount}
#) where {ValueType, DimensionCount}

#=
Outer constructor
From a mask.
Returning the bound box of trues.
=#
function BoundingBox(
        mask::BitArray{DimensionCount}
    ) where {DimensionCount}
    #=

    =#
    # All the indices, sorted
    indices = CartesianIndices(mask)
    # Assert they are in order !!!

    # see logical masking on an image in Julia, SO
    # An array can be indexed by a BitArray
    indicesOfTrue = indices[mask]

    result = BoundingBox(indicesOfTrue[1], indicesOfTrue[end])
    @debug "Bounding box of mask" result
    return result
end

#=
Returns a new bounding box,
expanding given bounding box in all dimension by amount.
(Is shrink if amount is negative.)

Ensures least is not less than (1,1,...,1)
but does not ensure most is constrained.
If you need constrained most, subsequently use clamp()
=#
function expandOrShrink(
        boundingBox::BoundingBox{DimensionCount},   # mutated
        amount
    ) where{DimensionCount}

    #=
    Iteration not supported on CartesianIndex,
    so cast to Tuple and then back to CartesianIndex.
    =#

    #=
    <Subtract> from upper left.
    Across all dimensions, and clamp to <1> the least index.
    =#
    least = CartesianIndex(
        max.( (Tuple(boundingBox.least) .- amount),
            1)
        )

    #=
    <Add> to lower right.
    Across all dimensions.
    !!! Not clamp to <max of dimension>
    =#
    #=
    most = CartesianIndex(
        min.( (Tuple(boundingBox.most) .+ amount),
            1)
        )
    =#
    most = CartesianIndex(
            (Tuple(boundingBox.most) .+ amount)
        )

    result = BoundingBox(least, most)
    @debug "Expanded/shrunk boundingBox" result

    #=
    Ensure the bounding box least is > 1's
    but bounding box most may be out of bounds.
    =#
    return result
end

#=
Returns CartesianIndex of clamped most.
=#
function clampMostOfBoundingBox(boundingBox1, boundingBox2)
    most = min.(    Tuple(boundingBox1.most),
                    Tuple(boundingBox2.most)
                )

    result = CartesianIndex(most)
    return result
end

#=
Return new boundingBox which is contained in the image.

Assert given bounding box least is contained in the image.
Only need to clamp the most.
=#
function clampBoundingBoxToImage(
        boundingBox1::BoundingBox{DimensionCount},
        image   # ::Array{}
    ) where{DimensionCount}

    # In general, least is not (1,1,..1)
    least = boundingBox1.least

    #=
    Most is the smaller of the bounding box of the mask
    and the bounding box of the image.
    =#
    # OLD most = boundingBoxConstraint.most
    most = clampMostOfBoundingBox(boundingBox1, BoundingBox(image))

    result = BoundingBox(least, most)
    @debug "Clamped bounding box" result
    return result
end

#=
Return BitArray similar to image,
having trues in the rectangle defined by boundingBox.
=#
function bitMaskFromInShapeOf(
        boundingBox,
        image
    )
    # falses of similar shape
    mask = falses(size(image))

    # trues in the boundingBox
    cartesianIndices = cartesianIndicesFrom(boundingBox)
    mask[cartesianIndices] .= true

    @debug "Mask from boundingBox" mask
    @assert typeof(mask) <: BitArray
    return mask
end

#=
Return CartesianIndices
=#
function cartesianIndicesFrom(
        boundingBox
    )
    # Fail: result = CartesianIndices((boundingBox.least, boundingBox.most))
    # TODO generic on DimensionCount
    result = CartesianIndices(
            (boundingBox.least[1]:boundingBox.most[1],
             boundingBox.least[2]:boundingBox.most[2]
            )
        )
    @debug "CI of boundingBox" typeof(result) result
    return result

end

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
using SafeTestsets

const GROUP = get(ENV, "GROUP", "All")
const is_APPVEYOR = Sys.iswindows() && haskey(ENV,"APPVEYOR")
const is_TRAVIS = haskey(ENV,"TRAVIS")

#Start Test Script

@time begin
if GROUP == "All" || GROUP == "Interface"
    @time @safetestset "Utilities Tests" begin include("utils.jl") end
    @time @safetestset "Poisson example" begin include("../examples/poisson.jl") end
    @time @safetestset "Heat equation example" begin include("../examples/heat_equation.jl") end
    @time @safetestset "Robin Boundary Condition Operators" begin include("robin.jl") end
    @time @safetestset "JacVec Operators Interface" begin include("jacvec_operators.jl") end
    @time @safetestset "Composite Operators Interface" begin include("composite_operators_interface.jl") end
    @time @safetestset "BC and Coefficient Compositions" begin include("bc_coeff_compositions.jl") end
    @time @safetestset "Derivative Operators Interface" begin include("derivative_operators_interface.jl") end
    @time @safetestset "Ghost Derivative Operators Interface" begin include("ghost_derivative_operators_interface.jl") end
    @time @safetestset "Validate Regular Derivative Operators" begin include("regular_operator_validation.jl") end
    @time @safetestset "Validate and Compare Generic Operators" begin include("generic_operator_validation.jl") end
    @time @safetestset "Validate Boundary Padded Array Concretization" begin include("boundary_padded_array.jl") end
    @time @safetestset "Validate Higher Dimensional Boundary Extension" begin include("MultiDimBC_test.jl") end
    #@time @safetestset "2nd order check" begin include("2nd_order_check.jl") end
    #@time @safetestset "KdV" begin include("KdV.jl") end # KdV times out and all fails
    #@time @safetestset "Heat Equation" begin include("heat_eqn.jl") end
    @time @safetestset "Matrix-Free Operators" begin include("matrixfree.jl") end
    @time @safetestset "JacVec Operator Integration Test" begin include("jacvec_integration_test.jl") end
    @time @safetestset "Convolutions" begin include("convolutions.jl") end
    @time @safetestset "Differentiation Dimension" begin include("differentiation_dimension.jl") end
    @time @safetestset "Higher Dimensional Concretization" begin include("concretization.jl") end
    @time @safetestset "Coefficient Functions" begin include("coefficient_functions.jl") end
    @time @safetestset "Upwind Operator Interface" begin include("upwind_operators_interface.jl") end
    @time @safetestset "MOLFiniteDifference Interface" begin include("MOLtest.jl") end
    @time @safetestset "Basic SDO Examples" begin include("BasicSDOExamples.jl") end
    # @time @safetestset "Linear Complementarity Problem Examples" begin include("lcp.jl"); include("lcp_split.jl") end
end

if !is_APPVEYOR && (GROUP == "All" || GROUP == "Multithreading")
    @time @safetestset "2D and 3D fast multiplication" begin include("2D_3D_fast_multiplication.jl") end
end
end

ext: jl


text:
using Libm
using BenchmarkTools

xerf1 = linspace(-0.84375,0.84375,10000)
xerf2 = linspace(-2e-28,2e-28,10000)
xerf3 = linspace(-0.84375,0.84375,10000)
xerf4 = linspace(0.84375, 6,10000)
xerf5 = linspace(-6,-0.84375,10000);
xerf = union(xerf1,xerf2,xerf3,xerf4,xerf5);

xerfc1 = linspace(-0.84375,0.84375,10000)
xerfc2 = linspace(-2e-56,2e-56,10000)
xerfc3 = linspace(-0.25,0.25,10000)
xerfc4 = linspace(0.84375, 28,10000)
xerfc5 = linspace(-28,-0.84375,10000);
xercf = union(xerfc1,xerfc2,xerfc3,xerfc4,xerfc5)

t1 = @benchmark Libm.erf.($xerf)
t2 = @benchmark Base.erf.($xerf)
t3 = @benchmark Libm.erfc.($xercf)
t4 = @benchmark Base.erfc.($xercf)

println("***** Benchmark erf *****")
println()
println(ratio(mean(t1),mean(t2)))
println()
println("Details\n")
println(t1,"\n")
println(t2,"\n")

println("***** Benchmark erfc *****")
println()
println(ratio(mean(t3),mean(t4)))
println()
println("Details\n")
println(t3,"\n")
println(t4,"\n")


text:
using LinearAlgebra: diag, norm
using Quantica: Hamiltonian, ParametricHamiltonian

@testset "basic hamiltonians" begin
    presets = (LatticePresets.linear, LatticePresets.square, LatticePresets.triangular,
               LatticePresets.honeycomb, LatticePresets.cubic, LatticePresets.fcc,
               LatticePresets.bcc)
    types = (Float16, Float32, Float64, ComplexF16, ComplexF32, ComplexF64)
    ts = (1, 2.0I, @SMatrix[1 2; 3 4], 1.0f0*I)
    orbs = (Val(1), Val(1), Val(2), (Val(1), Val(2)))
    for preset in presets, lat in (preset(), unitcell(preset())), type in types
        E, L = dims(lat)
        dn0 = ntuple(_ -> 1, Val(L))
        for (t, o) in zip(ts, orbs)
            @test hamiltonian(lat, onsite(t) + hopping(t; range = 1), orbitals = o, type = type) isa Hamiltonian
            @test hamiltonian(lat, onsite(t) - hopping(t; dn = dn0), orbitals = o, type = type) isa Hamiltonian
        end
    end
    h = LatticePresets.honeycomb() |> hamiltonian(hopping(1, range = 1/√3))
    @test bloch(h) == h.harmonics[1].h
    # Inf range
    h = LatticePresets.square() |> unitcell(region = RegionPresets.square(5)) |>
        hamiltonian(hopping(1, range = Inf))
    @test Quantica.nhoppings(h) == 600

    h = LatticePresets.square() |> hamiltonian(hopping(1, dn = (10,0), range = Inf))
    @test Quantica.nhoppings(h) == 1
    @test isassigned(h, (10,0))
end

@testset "similarmatrix" begin
    types = (ComplexF16, ComplexF32, ComplexF64)
    lat = LatticePresets.honeycomb()
    for T in types
        h0 = hamiltonian(lat, onsite(I) + hopping(2I; range = 1), orbitals = (Val(1), Val(2)), type = T)
        hf = flatten(h0)
        hm = Matrix(h0)
        hs = (h0, hf, hm)
        As = (SparseMatrixCSC, SparseMatrixCSC, Matrix)
        Es = (SMatrix{2,2,T,4}, T, SMatrix{2,2,T,4})
        for (h, A, E) in zip(hs, As, Es)
            sh = similarmatrix(h)
            @test sh isa A{E}
            b1 = bloch!(similarmatrix(flatten(h)), flatten(h), (1,1))
            b2 = bloch!(similarmatrix(h, AbstractMatrix{T}), h, (1,1))
            @test isapprox(b1, b2)
            for T´ in types
                E´s = E <: SMatrix ? (SMatrix{2,2,T´,4}, T´) : (T´,)
                for E´ in E´s
                    s1 = similarmatrix(h, Matrix{E´})
                    s2 = similarmatrix(h, AbstractMatrix{E´})
                    @test s1 isa Matrix{E´}
                    @test s2 isa A{E´}
                end
            end
        end
    end
end

@testset "orbitals and sublats" begin
    orbs = (:a, (:a,), (:a, :b, 3), ((:a, :b), :c), ((:a, :b), (:c,)), (Val(2), Val(1)),
            (:A => (:a, :b), :D => :c), :D => Val(4))
    lat = LatticePresets.honeycomb()
    for orb in orbs
        @test hamiltonian(lat, onsite(I), orbitals = orb) isa Hamiltonian
    end
    @test hamiltonian(lat, onsite(I) + hopping(@SMatrix[1 2], sublats = :A =>:B),
                      orbitals = :A => Val(2)) isa Hamiltonian
    h1 = hamiltonian(lat, onsite(I) + hopping(@SMatrix[1 2], sublats = :A =>:B),
                      orbitals = :A => Val(2))
    h2 = hamiltonian(lat, onsite(I) + hopping(@SMatrix[1 2], sublats = (:A =>:B,)),
                      orbitals = :A => Val(2))
    @test bloch(h1, (1, 2)) == bloch(h2, (1, 2))
end

@testset "onsite dimensions" begin
    lat = LatticePresets.honeycomb()
    ts = (@SMatrix[1 2; 3 4], r -> @SMatrix[r[1] 2; 3 4], 2, r -> 2r[1],
          @SMatrix[1 2; 3 4], r -> @SMatrix[r[1] 2; 3 4])
    os = (Val(1), Val(1), Val(2), Val(2), Val(3), Val(3))
    for (t, o) in zip(ts, os)
        model = onsite(t)
        @test_throws DimensionMismatch hamiltonian(lat, model, orbitals = o)
    end

    ts = (@SMatrix[1 2; 3 4], r -> @SMatrix[r[1] 2; 3 4], 2, r -> 2r[1], 3I, r -> I, 3I)
    os = (Val(2), Val(2), Val(1), Val(1), Val(3), Val(3), (Val(1), Val(3)))
    for (t, o) in zip(ts, os)
        model = onsite(t)
        @test hamiltonian(lat, model, orbitals = o) isa Hamiltonian
    end
    @test bloch(hamiltonian(lat, onsite(3I), orbitals = (Val(1), Val(3))))[1,1] ==
        @SMatrix[3 0 0; 0 0 0; 0 0 0]
    @test bloch(hamiltonian(lat, onsite(3I), orbitals = (Val(1), Val(3))))[2,2] ==
        SMatrix{3,3}(3I)
end

@testset "hopping dimensions" begin
    lat = LatticePresets.honeycomb()
    ts = (@SMatrix[1 2; 2 3], (r,dr) -> @SMatrix[r[1] 2; 3 4], 2, (r,dr) -> 2r[1],
          @SMatrix[1 2], @SMatrix[1 ;2], @SMatrix[1 2], @SMatrix[1 2; 2 3])
    os = (Val(1), Val(1), Val(2), Val(2), (Val(2), Val(1)), (Val(2), Val(1)),
         (Val(2), Val(1)), (Val(2), Val(1)))
    ss = (missing, missing, missing, missing, :B => :A, :A => :B, missing, missing)
    for (t, o, s) in zip(ts, os, ss)
        model = hopping(t, sublats = s)
        @test_throws DimensionMismatch hamiltonian(lat, model, orbitals = o)
    end
    ts = (@SMatrix[1 2], @SMatrix[1 ;2])
    os = ((Val(2), Val(1)), (Val(2), Val(1)))
    ss = (:A => :B, :B => :A)
    for (t, o, s) in zip(ts, os, ss)
        model = hopping(t, sublats = s)
        @test hamiltonian(lat, model, orbitals = o) isa Hamiltonian
    end
    @test bloch(hamiltonian(lat, hopping(3I, range = 1/√3), orbitals = (Val(1), Val(2))))[2,1] ==
        @SMatrix[3 0; 0 0]
    @test bloch(hamiltonian(lat, hopping(3I, range = 1/√3), orbitals = (Val(1), Val(2))))[1,2] ==
        @SMatrix[3 0; 0 0]
end

@testset "hermiticity" begin
    lat = LatticePresets.honeycomb()
    @test !ishermitian(hamiltonian(lat, hopping(im, sublats = :A=>:B)))
    @test !ishermitian(hamiltonian(lat, hopping(1, sublats = :A=>:B)))
    @test !ishermitian(hamiltonian(lat, hopping(1, sublats = :A=>:B, dn = (-1,0))))
    @test ishermitian(hamiltonian(lat, hopping(1, sublats = :A=>:B, dn = (1,0))))
    @test !ishermitian(hamiltonian(lat, hopping(im)))
    @test ishermitian(hamiltonian(lat, hopping(1)))

    @test ishermitian(hamiltonian(lat, hopping(im, sublats = :A=>:B, plusadjoint = true)))
    @test ishermitian(hamiltonian(lat, hopping(1, sublats = :A=>:B, plusadjoint = true)))
    @test ishermitian(hamiltonian(lat, hopping(1, sublats = :A=>:B, dn = (1,0), plusadjoint = true)))
    @test ishermitian(hamiltonian(lat, hopping(im, plusadjoint = true)))
    @test ishermitian(hamiltonian(lat, hopping(1, plusadjoint = true)))
end

@testset "unitcell modifiers" begin
    h = LatticePresets.honeycomb() |> hamiltonian(hopping(1) + onsite(0)) |> unitcell(2, modifiers = (@onsite!((o, r) -> 1), @hopping!(h -> 1)))
    @test diag(bloch(h)) == ComplexF64[1, 1, 1, 1, 1, 1, 1, 1]
end

@testset "@onsite!" begin
    el = @SMatrix[1 2; 2 1]

    @test @onsite!(o -> 2o)(el) == 2el
    @test @onsite!(o -> 2o)(el, p = 2) == 2el
    @test @onsite!((o;) -> 2o)(el) == 2el
    @test @onsite!((o;) -> 2o)(el, p = 2) == 2el
    @test @onsite!((o; z) -> 2o)(el, z = 2) == 2el
    @test @onsite!((o; z) -> 2o)(el, z = 2, p = 2) == 2el
    @test @onsite!((o; z = 2) -> 2o)(el) == 2el
    @test @onsite!((o; z = 2) -> 2o)(el, p = 2) == 2el
    @test @onsite!((o; z = 2) -> 2o)(el, z = 1, p = 2) == 2el
    @test @onsite!((o; kw...) -> 2o)(el) == 2el
    @test @onsite!((o; kw...) -> 2o)(el, p = 2) == 2el
    @test @onsite!((o; z, kw...) -> 2o)(el, z = 2) == 2el
    @test @onsite!((o; z, kw...) -> 2o)(el, z = 2, p = 2) == 2el
    @test @onsite!((o; z, y = 2, kw...) -> 2o)(el, z = 2, p = 2) == 2el
    @test @onsite!((o; z, y = 2, kw...) -> 2o)(el, z = 2, y = 3, p = 2) == 2el

    r = SVector(0,0)

    @test @onsite!((o, r;) -> 2o)(el, r) == 2el
    @test @onsite!((o, r;) -> 2o)(el, r, p = 2) == 2el
    @test @onsite!((o, r; z) -> 2o)(el, r, z = 2) == 2el
    @test @onsite!((o, r; z) -> 2o)(el, r, z = 2, p = 2) == 2el
    @test @onsite!((o; z = 2) -> 2o)(el, r) == 2el
    @test @onsite!((o; z = 2) -> 2o)(el, r, p = 2) == 2el
    @test @onsite!((o; z = 2) -> 2o)(el, r, z = 1, p = 2) == 2el
    @test @onsite!((o, r; kw...) -> 2o)(el, r) == 2el
    @test @onsite!((o, r; kw...) -> 2o)(el, r, p = 2) == 2el
    @test @onsite!((o, r; z, kw...) -> 2o)(el, r, z = 2) == 2el
    @test @onsite!((o, r; z, kw...) -> 2o)(el, r, z = 2, p = 2) == 2el
    @test @onsite!((o, r; z, y = 2, kw...) -> 2o)(el, r, z = 2, p = 2) == 2el
    @test @onsite!((o, r; z, y = 2, kw...) -> 2o)(el, r, z = 2, y = 3, p = 2) == 2el

    @test @onsite!((o; z, y = 2, kw...) -> 2o) isa Quantica.UniformOnsiteModifier
    @test @onsite!((o, r; z, y = 2, kw...) -> 2o) isa Quantica.OnsiteModifier{2}

    @test parameters(@onsite!((o, r; z, y = 2, kw...) -> 2o)) == (:z, :y)
end

@testset "@hopping!" begin
    el = @SMatrix[1 2; 2 1]

    @test @hopping!(t -> 2t)(el) == 2el
    @test @hopping!(t -> 2t)(el, p = 2) == 2el
    @test @hopping!((t;) -> 2t)(el) == 2el
    @test @hopping!((t;) -> 2t)(el, p = 2) == 2el
    @test @hopping!((t; z) -> 2t)(el, z = 2) == 2el
    @test @hopping!((t; z) -> 2t)(el, z = 2, p = 2) == 2el
    @test @hopping!((t; z = 2) -> 2t)(el) == 2el
    @test @hopping!((t; z = 2) -> 2t)(el, p = 2) == 2el
    @test @hopping!((t; z = 2) -> 2t)(el, z = 1, p = 2) == 2el
    @test @hopping!((t; kw...) -> 2t)(el) == 2el
    @test @hopping!((t; kw...) -> 2t)(el, p = 2) == 2el
    @test @hopping!((t; z, kw...) -> 2t)(el, z = 2) == 2el
    @test @hopping!((t; z, kw...) -> 2t)(el, z = 2, p = 2) == 2el
    @test @hopping!((t; z, y = 2, kw...) -> 2t)(el, z = 2, p = 2) == 2el
    @test @hopping!((t; z, y = 2, kw...) -> 2t)(el, z = 2, y = 3, p = 2) == 2el

    r = dr = SVector(0,0)

    @test @hopping!((t, r, dr;) -> 2t)(el, r, dr) == 2el
    @test @hopping!((t, r, dr;) -> 2t)(el, r, dr, p = 2) == 2el
    @test @hopping!((t, r, dr; z) -> 2t)(el, r, dr, z = 2) == 2el
    @test @hopping!((t, r, dr; z) -> 2t)(el, r, dr, z = 2, p = 2) == 2el
    @test @hopping!((t; z = 2) -> 2t)(el, r, dr) == 2el
    @test @hopping!((t; z = 2) -> 2t)(el, r, dr, p = 2) == 2el
    @test @hopping!((t; z = 2) -> 2t)(el, r, dr, z = 1, p = 2) == 2el
    @test @hopping!((t, r, dr; kw...) -> 2t)(el, r, dr) == 2el
    @test @hopping!((t, r, dr; kw...) -> 2t)(el, r, dr, p = 2) == 2el
    @test @hopping!((t, r, dr; z, kw...) -> 2t)(el, r, dr, z = 2) == 2el
    @test @hopping!((t, r, dr; z, kw...) -> 2t)(el, r, dr, z = 2, p = 2) == 2el
    @test @hopping!((t, r, dr; z, y = 2, kw...) -> 2t)(el, r, dr, z = 2, p = 2) == 2el
    @test @hopping!((t, r, dr; z, y = 2, kw...) -> 2t)(el, r, dr, z = 2, y = 3, p = 2) == 2el

    @test @hopping!((t; z, y = 2, kw...) -> 2t) isa Quantica.UniformHoppingModifier
    @test @hopping!((t, r, dr; z, y = 2, kw...) -> 2t) isa Quantica.HoppingModifier{3}

    @test parameters(@hopping!((o, r, dr; z, y = 2, kw...) -> 2o)) == (:z, :y)
end

@testset "parametric" begin
    h = LatticePresets.honeycomb() |> hamiltonian(hopping(1) + onsite(2)) |> unitcell(10)
    T = typeof(h)
    @test parametric(h, @onsite!(o -> 2o))() isa T
    @test parametric(h, @onsite!((o, r) -> 2o))() isa T
    @test parametric(h, @onsite!((o, r; a = 2) -> a*o))() isa T
    @test parametric(h, @onsite!((o, r; a = 2) -> a*o))(a=1) isa T
    @test parametric(h, @onsite!((o, r; a = 2) -> a*o), @hopping!(t -> 2t))(a=1) isa T
    @test parametric(h, @onsite!((o, r) -> o), @hopping!((t, r, dr) -> r[1]*t))() isa T
    @test parametric(h, @onsite!((o, r) -> o), @hopping!((t, r, dr; a = 2) -> r[1]*t))() isa T
    @test parametric(h, @onsite!((o, r) -> o), @hopping!((t, r, dr; a = 2) -> r[1]*t))(a=1) isa T
    @test parametric(h, @onsite!((o, r; b) -> o), @hopping!((t, r, dr; a = 2) -> r[1]*t))(b=1) isa T
    @test parametric(h, @onsite!((o, r; b) -> o*b), @hopping!((t, r, dr; a = 2) -> r[1]*t))(a=1, b=2) isa T

    # Issue #35
    for orb in (Val(1), Val(2))
        h = LatticePresets.triangular() |> hamiltonian(hopping(I) + onsite(I), orbitals = orb) |> unitcell(10)
        ph = parametric(h, @onsite!((o, r; b) -> o+b*I), @hopping!((t, r, dr; a = 2) -> t+r[1]*I),
                       @onsite!((o, r; b) -> o-b*I), @hopping!((t, r, dr; a = 2) -> t-r[1]*I))
        @test isapprox(bloch(ph(a=1, b=2), (1, 2)), bloch(h, (1, 2)))
    end
    # Issue #37
    for orb in (Val(1), Val(2))
        h = LatticePresets.triangular() |> hamiltonian(hopping(I) + onsite(I), orbitals = orb) |> unitcell(10)
        ph = parametric(h, @onsite!(o -> o*cis(1)))
        @test ph()[1,1] ≈ h[1,1]*cis(1)
    end
    # Issue #54. Parametric Haldane model
    sK(dr::SVector) = sK(atan(dr[2],dr[1]))
    sK(ϕ) = 2*mod(round(Int, 6*ϕ/(2π)), 2) - 1
    ph = LatticePresets.honeycomb() |> hamiltonian(hopping(1)) |>
         parametric(@hopping!((t, r, dr; λ) ->  λ*im*sK(dr); sublats = :A=>:A),
                    @hopping!((t, r, dr; λ) -> -λ*im*sK(dr); sublats = :B=>:B))
    @test bloch(ph(λ=1), (π/2, -π/2)) == bloch(ph, (1, π/2, -π/2)) ≈ [4 1; 1 -4]
    # Non-numeric parameters
    ph = LatticePresets.honeycomb() |> hamiltonian(hopping(1)) |>
         parametric(@hopping!((t, r, dr; λ, k) ->  λ*im*sK(dr+k); sublats = :A=>:A),
                    @hopping!((t, r, dr; λ, k) -> -λ*im*sK(dr+k); sublats = :B=>:B))
    @test bloch(ph(λ=1, k=SA[1,0]), (π/2, -π/2)) == bloch(ph, (1, SA[1,0], π/2, -π/2)) ≈ [-4 1; 1 4]
    # Issue 61, type stability
    h = LatticePresets.honeycomb() |> hamiltonian(onsite(0))
    @inferred parametric(h, @onsite!((o;μ) -> o- μ))
    @inferred parametric(h, @onsite!(o->2o), @hopping!((t)->2t), @onsite!((o, r)->o+r[1]))
    @inferred parametric(h, @onsite!((o, r)->o*r[1]), @hopping!((t; p)->p*t), @onsite!((o; μ)->o-μ))
end

@testset "boolean masks" begin
    for b in ((), (1,1), 4)
        h1 = LatticePresets.honeycomb() |> hamiltonian(hopping(1) + onsite(2)) |>
             supercell(b, region = RegionPresets.circle(10))
        h2 = LatticePresets.honeycomb() |> hamiltonian(hopping(1) + onsite(2)) |>
             supercell(b, region = RegionPresets.circle(20))

        @test isequal(h1 & h2, h1)
        @test isequal(h1, h2) || !isequal(h1 & h2, h2)
        @test isequal(h1, h2) || !isequal(h1 | h2, h1)
        @test  isequal(h1 | h2, h2)

        @test isequal(unitcell(h1 & h2), unitcell(h1))
        @test isequal(h1, h2) || !isequal(unitcell(h1 & h2), unitcell(h2))
        @test isequal(h1, h2) || !isequal(unitcell(h1 | h2), unitcell(h1))
        @test isequal(unitcell(h1 | h2), unitcell(h2))

        h1 = h1.lattice
        h2 = h2.lattice

        @test isequal(h1 & h2, h1)
        @test isequal(h1, h2) || !isequal(h1 & h2, h2)
        @test isequal(h1, h2) || !isequal(h1 | h2, h1)
        @test  isequal(h1 | h2, h2)

        @test isequal(unitcell(h1 & h2), unitcell(h1))
        @test isequal(h1, h2) || !isequal(unitcell(h1 & h2), unitcell(h2))
        @test isequal(h1, h2) || !isequal(unitcell(h1 | h2), unitcell(h1))
        @test isequal(unitcell(h1 | h2), unitcell(h2))
    end
end

@testset "unitcell seeds" begin
    p1 = SA[100,0]
    p2 = SA[0,20]
    lat = LatticePresets.honeycomb()
    model = hopping(1, range = 1/√3) + onsite(2)
    h1 = lat |> hamiltonian(model) |> supercell(region = r -> norm(r-p1)<3, seed = p1)
    h2 = lat |> hamiltonian(model) |> supercell(region = r -> norm(r-p2)<3, seed = p2)
    h = unitcell(h1 | h2)
    h3 = lat |> hamiltonian(model) |> unitcell(region = r -> norm(r-p1)<3 || norm(r-p2)<3, seed = p2)

    @test Quantica.nsites(h) == 130
    @test Quantica.nsites(h3) == 64
end


text:
abstract type RegressionModels end
#If only we could add methods to abstract types...
# (M::RegressionModel)(X) = RegressionOut(X, M)
# maybe in Julia 2.0? - Update: Julia 1.2+ allows this!!!

struct ClassicLeastSquares <: RegressionModels
    Coefficients::Array
    Bias::Bool
end

struct RidgeRegression <: RegressionModels
    CLS::ClassicLeastSquares
end

"""
    ClassicLeastSquares( X, Y; Bias = false )

Constructs a ClassicLeastSquares regression model of the form `Y` = A`X` with or without a `Bias` term. Returns a CLS object.
"""
function ClassicLeastSquares( X::Array, Y::Array; Bias = false )
    Z = (Bias) ? hcat( ones( size( X )[ 1 ] ), X ) : X
    return ClassicLeastSquares(Base.inv(Z' * Z) * Z' * Y, Bias)
end

"""
    OrdinaryLeastSquares( X, Y; Bias = false )

Makes a ClassicLeastSquares regression model of the form `Y` = A`X` with or without a `Bias` term. Returns a CLS object.
This is a wrapper function for CLS, because most other fields refer to this as OLS.
"""
function OrdinaryLeastSquares( X, Y; Bias = false )
    Z = (Bias) ? hcat( ones( size( X )[ 1 ] ), X ) : X
    return ClassicLeastSquares(Base.inv(Z' * Z) * Z' * Y, Bias)
end

function PredictFn(X, M::Union{ClassicLeastSquares, RidgeRegression})
    Z = ( M.Bias ) ? hcat( ones(size( X )[ 1 ] ), X ) : X
    return Z * M.Coefficients
end

"""
    (M::ClassicLeastSquares)(X)

Makes an inference from `X` using a ClassicLeastSquares object.
"""
(M::ClassicLeastSquares)(X) = PredictFn(X, M)

"""
    RidgeRegression( X, Y, Penalty; Bias = false )

Makes a RidgeRegression model of the form `Y` = A`X` with or without a `Bias` term and has an L2 `Penalty`. Returns a CLS object.
"""
function RidgeRegression( X, Y, Penalty; Bias = false )
    Y = forceMatrix(Y)
    Z = (Bias) ? hcat( ones( size( X )[ 1 ] ), X ) : X
    return RidgeRegression( ClassicLeastSquares( Base.inv( (Z' * Z) .+ (Penalty .* Diagonal( ones( size(Z)[2] ) ) ) ) * Z' * Y, Bias) )
end

"""
    (M::RidgeRegression)(X)

Makes an inference from `X` using a RidgeRegression object which wraps a ClassicLeastSquares object.
"""
(M::RidgeRegression)(X) = PredictFn(X, M.CLS)

struct KRR
    kernel::Kernel
    RR::RidgeRegression
end

"""
    KernelRidgeRegression( X, Y, Penalty; KernelParameter = 0.0, KernelType = "linear" )

Makes a KernelRidgeRegression model of the form `Y` = A`K` using a user specified Kernel("Linear", or "Guassian") and has an L2 `Penalty`.
Returns a KRR Wrapper for a CLS object.
"""
function KernelRidgeRegression( X, Y, Penalty; KernelParameter = 0.0, KernelType = "linear" )
    Kern = Kernel( KernelParameter, KernelType, X )
    return KRR(Kern, RidgeRegression( Kern(X), Y, Penalty ) )
end

"""
    (M::KRR)(X)

Makes an inference from `X` using a KRR object which wraps a ClassicLeastSquares object.
"""
(M::KRR)(X) = PredictFn(M.kernel(X), M.RR.CLS)

struct LSSVM
    kernel::Kernel
    RR::RidgeRegression
end

function formatlssvminput(X)
    Y = zeros(size(X) .+ 1)
    Y[ 1, 1 ] = 1.0
    Y[ 2:end, 2:end ] .= X
    return Y
end

"""
    LSSVM( X, Y, Penalty; KernelParameter = 0.0, KernelType = "linear" )

Makes a LSSVM model of the form `Y` = A`K` with a bias term using a user specified Kernel("linear", or "gaussian") and has an L2 `Penalty`.
Returns a LSSVM Wrapper for a CLS object.
"""
function LSSVM( X, Y, Penalty; KernelParameter = 0.0, KernelType = "linear" )
    Kern = Kernel( KernelParameter, KernelType, X )
    return LSSVM(   Kern,
                    RidgeRegression( formatlssvminput( Kern( X ) ),
                    vcat( 0.0, Y ),
                    Penalty ) )
end

"""
    (M::LSSVM)(X)

Makes an inference from `X` using a LSSVM object.
"""
(M::LSSVM)(X) = PredictFn( formatlssvminput( M.kernel(X)), M.RR.CLS )[2:end,:]

struct PrincipalComponentRegression <: RegressionModels
    PCA::PCA
    CLS::ClassicLeastSquares
end

"""
    PrincipalComponentRegression(PCAObject, Y )

Makes a PrincipalComponentRegression model object from a PCA Object and property value `Y`.
"""
function PrincipalComponentRegression(PCAObject::PCA, Y )
    return PrincipalComponentRegression(PCAObject, ClassicLeastSquares( PCAObject.Scores, forceMatrix(Y) ) )
end

PredictFn(X, M::PrincipalComponentRegression) = M.CLS( M.PCA( X ) )

"""
    (M::PrincipalComponentRegression)( X )

Makes an inference from `X` using a PrincipalComponentRegression object.
"""
(M::PrincipalComponentRegression)( X ) = PredictFn( X, M )

struct PartialLeastSquares <: RegressionModels
    XLoadings::Array
    XScores::Array
    YLoadings::Array
    YScores::Array
    XWeights::Array
    RWeights::Array
    Coefficients::Array
    Factors::Int
    XVariance::Array
    YVariance::Array
end

"""
    PartialLeastSquares( X, Y; Factors = minimum(size(X)) - 2, tolerance = 1e-8, maxiters = 200 )

Returns a PartialLeastSquares regression model object from arrays `X` and `Y`.

1. PARTIAL LEAST-SQUARES REGRESSION: A TUTORIAL PAUL GELADI and BRUCE R.KOWALSKI. Analytica Chimica Acta, 186, (1986) PARTIAL LEAST-SQUARES REGRESSION:
2. Martens H., NÊs T. Multivariate Calibration. Wiley: New York, 1989.
3. Re-interpretation of NIPALS results solves PLSR inconsistency problem. Rolf Ergon. Published in Journal of Chemometrics 2009; Vol. 23/1: 72-75
"""
function PartialLeastSquares( X, Y; Factors = minimum(size(X)) - 2, tolerance = 1e-8, maxiters = 200 )
    (Xrows, Xcols) = size(X)
    Y = forceMatrix(Y)
    (Yrows, Ycols) = size(Y)
    @assert Factors < (reduce(min, (Xrows, Xcols) ) - 1)
    Xd = copy(X); Yd = copy(Y)
    Xvar = Statistics.var(X, dims = 1)
    Yvar = Statistics.var(Y, dims = 1)
    Coefficients = []
    T = zeros(Xrows, Factors); t = zeros(Xrows); tprime = zeros(Xrows);
    U = zeros(Yrows, Factors); u = zeros(Yrows)
    P = zeros(Xcols, Factors); p = zeros(Xcols);
    Q = zeros(Ycols, Factors); q = zeros(Ycols)
    W = zeros(Xcols, Factors); w = zeros(Xrows);
    for factor in 1:Factors
        u = (Ycols == 1) ? Yd[:,1] : Yd[ :, argmax( [ (Yd[:,col]' * Yd[:,col])[1] for col in 1:Ycols] ) ]
        for iter in 1:maxiters
            w = X' * u
            w = w ./ sqrt.(w' * w)
            t = Xd * w
            q = Y' * (t ./ (t' * t) )
            diff = t .- tprime
            if (Ycols == 1) || ( sqrt.( diff' * diff )[1] < tolerance )
                break
            end
            u = Y * (q ./ (q' * q))
            tprime = t;
        end#End for ALS iterations
        tnorm = t ./ ( t' * t )
        p = X' * tnorm
        Xd = Xd .- ( t * p' )
        Yd = Yd .- ( t * q' )
        #Update Model Variables
        Q[:,factor] = q; T[:,factor] = t
        U[:,factor] = u; P[:,factor] = p
        W[:, factor] = w
    end#end for factors
    R = (Factors == 1) ? W : W * Base.inv( P' * W )
    #Use a more modern way to solve for the regression coefficients (2)
    Coefficients = R * Q'
    #An Equivalent way to obtain the regression coefficients.
    #Coefficients = W*Base.inv(W'*X'*X*W)*W'*X'*Y*Q
    return PartialLeastSquares(P, T, Q, U, W, R, Coefficients, Factors, Xvar, Yvar)
end

function PredictFn(X, M::PartialLeastSquares; Factors)
    Coeffs = []
    if Factors < M.Factors
        R = (Factors == 1) ? M.XWeights[:, 1:Factors] : M.XWeights[:, 1:Factors] * Base.inv( M.XLoadings[:, 1:Factors]' * M.XWeights[:, 1:Factors]  )
        Coeffs = R * M.YLoadings[:, 1:Factors]'
    else
        Coeffs = M.Coefficients
    end
    return X * Coeffs
end

"""
    (M::PartialLeastSquares)

Makes an inference from `X` using a PartialLeastSquares object.
"""
(M::PartialLeastSquares)(X; Factors = M.Factors) = PredictFn(X, M; Factors = Factors)

struct ELM <: RegressionModels
    Reservoir::Array
    Coefficients::Array
    Fn::Function
end

"""
    sigmoid(x)

Applies the sigmoid function to a scalar value X. Returns a scalar. Can be broad-casted over an Array.
"""
sigmoid(x) = 1.0 / (1.0 + exp(-1.0 * x))

"""
    ExtremeLearningMachine(X, Y, ReservoirSize = 10; ActivationFn = sigmoid)

Returns a ELM regression model object from arrays `X` and `Y`, with a user specified `ReservoirSize` and `ActivationFn`.

Extreme learning machine: a new learning scheme of feedforward neural networks. Guang-Bin Huang ; Qin-Yu Zhu ; Chee-Kheong Siew. 	2004 IEEE International Joint...
"""
function ExtremeLearningMachine(X, Y, ReservoirSize = 10; ActivationFn = sigmoid)
    W = randn( size( X )[2], ReservoirSize )
    return ELM( W, LinearAlgebra.pinv( ActivationFn.(X * W) ) * Y,
                ActivationFn)
end

PredictFn(X, M::ELM) = M.Fn.(X * M.Reservoir) * M.Coefficients;

"""
    (M::ELM)(X)

Makes an inference from `X` using a ELM object.
"""
(M::ELM)(X) = PredictFn(X, M)

mutable struct block
    value::Float64 #Avg value
    weight::Float64
    size::Int32
    previous::Int32
    next::Int32
end

"""
    MonotoneRegression(x, w = nothing)

Performs a monotone/isotonic regression on a vector x. This can be weighted
with a vector w.

Code was translated directly from:
Exceedingly Simple Monotone Regression. Jan de Leeuw. Version 02, March 30, 2017
"""
function MonotoneRegression(x::Array{Float64,1}, w = nothing)
    bins = length(x)
    retx = ones(bins)
    retx[:] = x[:]
    if isa(w, Nothing)
        w = ones(bins)
    end
    bins = bins + 1
    blocks = [ block(x[i], w[i], 1, i-1, i+1) for i in 1:(bins-1)]
    #I have a one off error somewhere an easy solution is just to pad
    push!(blocks, block( 0.0, w[end], 1, bins-1, bins+1) )
    active = 1
    continue_reg = true
    while( continue_reg )
        upsatisfied = false
        next = blocks[active].next
        if next == bins
            upsatisfied = true
        elseif blocks[next].value > blocks[active].value
            upsatisfied = true
        end
        if !upsatisfied
            ww = blocks[active].weight + blocks[next].weight
            nextnext = blocks[next].next
            wxactive = blocks[active].weight * blocks[active].value
            wxnext = blocks[next].weight * blocks[next].value
            blocks[active].value = (wxactive + wxnext) / ww
            blocks[active].weight = ww
            blocks[active].size += blocks[next].size
            blocks[active].next = nextnext
            if nextnext < bins
                blocks[nextnext].previous = active
            end
            blocks[next].size = 0
        end
        downsatisfied = false
        previous = blocks[active].previous

        if previous == 0
            downsatisfied = true
        elseif blocks[previous].value < blocks[active].value
            downsatisfied = true
        end
        if !downsatisfied
            ww = blocks[active].weight + blocks[previous].weight
            previousprevious = blocks[previous].previous
            wxactive = blocks[active].weight * blocks[active].value
            wxprevious = blocks[previous].weight * blocks[previous].value
            blocks[active].value = (wxactive + wxprevious) / ww
            blocks[active].weight = ww
            blocks[active].size += blocks[previous].size
            blocks[active].previous = previousprevious
            if previousprevious > 0
                blocks[previousprevious].next = active
            end
            blocks[previous].size = 0
        end
        if (blocks[active].next == bins) && (downsatisfied)
            continue_reg = false
        end
        if (upsatisfied && downsatisfied)
            active = next
        end
    end
    k = 1
    for i in 1:bins
        blksize = blocks[i].size;
        if (blksize > 0) && (i < bins )
            retx[k:(k+(blksize-1))] .= blocks[i].value;
            k += blksize
        end
    end
    return retx
end


text:
using BinDeps

#=
@BinDeps.setup
cmdstan = library_dependency("homebrew/science/cmdstan")

if Sys.isapple()
    using Homebrew

    provides(Homebrew.HB, "homebrew/science/cmdstan", cmdstan, os = :Darwin)
end
=#


text:
module TestExamples
using Test
dir = joinpath("..", "examples")
@testset "example $filename" for filename in readdir(dir)
    path = joinpath(dir, filename)
    @eval module $(Symbol("TestExample_$filename"))
    include($path)
    end
end
end#module


text:
module BenchmarkHelper

using BenchmarkTools

using Pkg
using GitHub

export PROJECT_DIR, BenchmarkSuite,
    benchmarks_info, reply_comment,
    benchmark_files, run_benchmarks, get_results

include("ns-utilities.jl")

const PROJECT_DIR = abspath(@__DIR__) |> dirname
const CURRENT_BRANCH = Ref("HEAD")
const CURRENT_JOB_ID = Ref("NEW-ID")
const CONFIG = Pkg.TOML.parsefile(joinpath(PROJECT_DIR, "benchmarks/config.toml"))

# Global Benchmark Suite
const BenchmarkSuite = BenchmarkGroup()

# Functions run in benchmark processes

function set_benchmark_info(branch, job_id)
    CURRENT_BRANCH[] = branch
    CURRENT_JOB_ID[] = job_id
end

function output_path()
    branch = replace(CURRENT_BRANCH[], r"[-/.]"=>"_")
    outdir= joinpath(PROJECT_DIR, "benchmarks", "output")
    mkpath(outdir)
    return joinpath(outdir, "BM-$(CURRENT_JOB_ID[])-$branch.txt")
end

function run_suite(tags_expr=:ALL)
    benchmarks = eval(Meta.parse("BenchmarkSuite[@tagged $tags_expr]"))
    t = run(benchmarks, samples=1)
    branch = CURRENT_BRANCH[]
    @info "benchmark result", t
    open(output_path(), "w") do io
        show(io, bmresults_to_dataframe(t))
    end
end

# Functions run in the dispather process

function gitcurrentbranch()
    branches = readlines(`git branch`)
    ind = findfirst(x -> x[1] == '*', branches)
    brname = branches[ind][3:end]
    if brname[1] == '('
        brname = readlines(`git rev-parse --short HEAD`)[1]
    end
    return brname
end

function command_info(event_data)
    comment_body = event_data["comment"]["body"]
    !occursin("@BayesBot", comment_body) && return (nothing, nothing, nothing)
    reg_bm_cmd = r"`\w+\(.*\)`"i
    bm_cmd = match(reg_bm_cmd, comment_body)
    bm_cmd == nothing && return (nothing, nothing, nothing)
    cmd, args, kwargs = parse_submission_string(bm_cmd.match)
    return (cmd, args, kwargs)
end

function benchmarks_info(event_data)
    cmd, args, kwargs = command_info(event_data)
    if cmd != "runbenchmarks"
        return false, nothing, nothing
    end

    branches =  haskey(kwargs, :vs) ? Meta.parse(kwargs[:vs]) : []
    if branches isa Expr && branches.head == :vect
        branches = branches.args
    elseif branches isa AbstractString
        branches = [branches]
    end
    target = "master"
    if haskey(event_data, "pull_request")
        target = event_data["pull_request"]["head"]["ref"]
    end
    if isempty(branches)
        if target != "master"
            branches = [target, "master"]
        end
    elseif length(branches) == 1
        if !in(target, branches)
            push!(branches, target)
        elseif !in("master", branches)
            push!(branches, "master")
        end
    end
    return true, Meta.parse(args[1]), branches
end

function target_benchmarks(bm_text)
    return filter(!isempty, map(strip, split(bm_text, ",")))
end

function reply_comment(event_data, body)
    repo = event_data["repository"]["full_name"]
    # auth = GitHub.authenticate(get(ENV, "GITHUB_TOKEN", ""))
    auth = GitHub.authenticate(get(ENV, "BENCHMARK_KEY", ""))
    params = Dict("body" => body)

    comment_kind = :issue
    comment_parent = nothing

    if haskey(event_data, "issue")
        comment_kind = :issue
        comment_parent = event_data["issue"]["number"]
    elseif haskey(event_data["comment"], "commit_id")
        comment_kind = :commit
        comment_parent = event_data["comment"]["commit_id"]
    end

    create_comment(repo, comment_parent, comment_kind;
                   params=params, auth=auth)
end

function benchmark_excluded(file)
    excludes = joinpath(PROJECT_DIR, "benchmarks", "excludes.txt")
    patterns = get(CONFIG, "EXCLUDES", [])
    for pat in patterns
        match(Regex(pat), file) != nothing && return true
    end
    return false
end

function benchmark_files()
    bm_files = []
    for (root, dirs, files) in walkdir(joinpath(PROJECT_DIR, "benchmarks"))
        for file in files
            push!(bm_files,joinpath(root, file))
        end
    end

    bm_files = map(bm_files) do path
        replace(path, joinpath(PROJECT_DIR, "benchmarks/") => "")
    end

    filter(bm_files) do file
        endswith(file, ".jl") && !benchmark_excluded(file)
    end
end

function run_benchmarks(tags, branches)
    job_id = CURRENT_JOB_ID[]
    for branch in branches
        branch = replace(branch, r"^.*?[:@#]"i => "")
        run(`git checkout .`)
        run(`git checkout $branch`)
        includes = ""
        for bm in benchmark_files()
            includes *= """include("$(joinpath("benchmarks", bm))")\n"""
        end

        code = """
        using Pkg;
        pkg"instantiate; add JSON GitHub DataFrames BenchmarkTools;"
        Pkg.build(verbose=true)
        push!(LOAD_PATH, joinpath("$(PROJECT_DIR)", "benchmarks"))
        using BenchmarkHelper
        BenchmarkHelper.set_benchmark_info("$(gitcurrentbranch())", "$(job_id)")
        $(includes)
        BenchmarkHelper.run_suite(:($tags))
        """
        job = `julia --project=$(PROJECT_DIR) -e $code`
        @info(job);
        run(job)
    end
end

function get_results(job_id, user)
    # TODO
    body = """
### Benchmark Results

Hi @$user, here are the results of your demanded benchmarks:
"""

    outdir= joinpath(PROJECT_DIR, "benchmarks", "output")
    for file in readdir(outdir)
        body *= "\nContent of file `$file`:\n"
        body *= "\n```\n"
        body *= read(joinpath(outdir, file), String)
        body *= "\n```\n"
    end
    return body
end

end


text:
using Test
using Latexify
using ModelingToolkit
using ReferenceTests

### Tips for generating latex tests:
### Latexify has an unexported macro:
###
### Latexify.@generate_test
###
### which generates a test using a given latexify function.
### For example:
###
### Latexify.@generate_test latexify([1, 2, 3], [4, 5, 6]; env=:mdtable)
###
### This puts a ready-made test in your clipboard which you can paste into the
### test file.
###
### Just be sure to remove all such macros before you commit a change since it
### will cause issues with Travis.

@parameters t σ ρ β
@variables x(t) y(t) z(t)
D = Differential(t)

eqs = [D(x) ~ σ*(y-x)*D(x-y)/D(z),
       0 ~ σ*x*(ρ-z)/10-y,
       D(z) ~ x*y^(2//3) - β*z]

# Latexify.@generate_test latexify(eqs)
@test_reference "latexify/10.tex" latexify(eqs)

@variables u[1:3](t)
@parameters p[1:3]
eqs = [D(u[1]) ~ p[3]*(u[2]-u[1]),
       0 ~ p[2]*p[3]*u[1]*(p[1]-u[1])/10-u[2],
       D(u[3]) ~ u[1]*u[2]^(2//3) - p[3]*u[3]]

@test_reference "latexify/20.tex" latexify(eqs)

eqs = [D(u[1]) ~ p[3]*(u[2]-u[1]),
       D(u[2]) ~ p[2]*p[3]*u[1]*(p[1]-u[1])/10-u[2],
       D(u[3]) ~ u[1]*u[2]^(2//3) - p[3]*u[3]]

@test_reference "latexify/30.tex" latexify(eqs)
@parameters t
@variables x(t)
D = Differential(t)
eqs = [D(x) ~ (1+cos(t))/(1+2*x)]

@test_reference "latexify/40.tex" latexify(eqs)

max_stars_repo_name: paulrschrater/ModelingToolsJulia

max_issues_repo_name: paulrschrater/ModelingToolsJulia

max_forks_repo_name: paulrschrater/ModelingToolsJulia


text:
"""
  `Boundary`

Structure holding sets of node and edge indices for one particular physical boundary.
"""
mutable struct Boundary
  Nodes::Set{Int64}
  Edges::Set{Int64}
end

"""
  `Mesh`

Structure for a triangular finite element mesh with volume and boundary markers.
"""
mutable struct Mesh
  "number of nodes"
  nnodes::Int64

  "number of triangles"
  nelems::Int64

  "number of physical (marked) edges"
  nedges::Int64

  "list of node coordinates"
  Nodes::Array{Array{Float64,1},1}

  "list of triangle node indices"
  Triangles::Array{Array{Int64,1},1}

  "list of physical edges"
  Edges::Array{Array{Int64,1},1}

  "dictionary of marked boundaries"
  Boundaries::Dict{Int64,Boundary}

  "dictionary of marked volume regions"
  Volumes::Dict{Int64,Set{Int64}}
end

"""
  `PDESystem`

Structure holding all information to describe simple PDEs with Dirichlet realized
resolved by multipliers.
"""
mutable struct PDESystem
  "stiffness matrix"
  A::SparseMatrixCSC{Float64,Int64}

  "load vector"
  b::Vector{Float64}

  "Dirichlet boundary values"
  bc::Vector{Float64}

  "Dirichlet nodes"
  DI::Set{Int64}

  "Dirichlet indices, in case of qdim!=1"
  vec_ind::Set{Int64}

  "dimension of vector-valued state"
  qdim::Int64

  "system matrix factorization"
  Factors

  "system of stiffness matrix and Dirichlet conditinons"
  SystemMatrix::SparseMatrixCSC{Float64,Int64}

  "Dirichlet projection matrix"
  B::SparseMatrixCSC{Float64,Int64}

  "solution vector"
  state::Vector{Float64}

  "rhs vector with source term and dirichlet conditions"
  rhs::Vector{Float64}

end


text:
abstract type AbstractDataset end
abstract type AbstractDataPoint <: Flux3D.AbstractDataPoint end
abstract type AbstractCustomDataset <: AbstractDataset end

Base.size(v::AbstractCustomDataset) = (v.length,)
Base.length(v::AbstractCustomDataset) = v.length

Base.size(v::AbstractDataset) = error("define Base.size of dataset type: $(typeof(v)).")
Base.length(v::AbstractDataset) = error("define Base.length of dataset type: $(typeof(v)).")

const default_root = normpath(@__DIR__, "..", "..", "datasets")


text:
module Pages

using HttpServer, WebSockets, URIParser, Mustache, JSON, DataFrames
using Compat; import Compat: String, @static

export Endpoint, Session, Callback, Request, URI, query_params

type Session
    id::String
    route::String
    client::WebSocket

    function Session(route::String)
        id = uppercase(randstring(16))
        while haskey(sessions,id)
            id = uppercase(randstring(16))
        end
        s = new(id,route)
        sessions[id] = s
        pages[route].sessions[id] = s
        s
    end

    function Session()
        id = uppercase(randstring(16))
        while haskey(sessions,id)
            id = uppercase(randstring(16))
        end
        s = new(id)
        sessions[id] = s
        s
    end
end
function Base.show(io::Base.IO,session::Session)
    print(io,"Session: ",
        "\n  ID: ",session.id,
        "\n  Route: ",isdefined(session,:route) ? session.route : "")
end
const sessions = Dict{String,Session}()

type Endpoint
    handler::Function
    route::String
    sessions::Dict{String,Session}

    function Endpoint(handler,route)
        p = new(handler,route,Dict{String,WebSocket}())
        !haskey(pages,route) || warn("Page $route already exists.")
        pages[route] = p
        finalizer(p, p -> delete!(pages, p.route))
        p
    end
end
function Base.show(io::Base.IO,endpoint::Endpoint)
    print(io,"Endpoint created at $(endpoint.route).")
end
const pages = Dict{String,Endpoint}() # url => page

include("callbacks.jl")
include("server.jl")
include("api.jl")
include("ijulia.jl")

end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
using QuadGK

# %%
a = 1.0
f = x -> sin(a*x)/(a*x)
@time quadgk(f, -1e8, 0, 1e8)
@time quadgk(f, -1e8, 0, 1e8)
@time quadgk(f, -1e8, 0, 1e8)

# %%
dump(f)

# %%
a = 1.0
g = let a = a; x -> sin(a*x)/(a*x) end
@time quadgk(g, -1e8, 0, 1e8)
@time quadgk(g, -1e8, 0, 1e8)
@time quadgk(g, -1e8, 0, 1e8)

# %%
dump(g)

# %%
function F(a = 1.0)
    quadgk(x -> sin(a*x)/(a*x), -1e8, 0, 1e8)
end

@time F()
@time F()
@time F()

# %%

max_stars_repo_path: 0020/function containing a global variable.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
# ============================
# Structs to define a zone
# ============================

"""
    Zone

Struct that defines a zone.

### Notes

The `zone_id` is a unique identifier denoting the geographical planning area
into which the stop falls. The numeral before the dash denotes a high-level
planning zone. The text after the dash denotes the subzone within the high-level zone.

### Examples

For example, "D-16.2C" corresponds to "D" = high-level planning zone
and "16.2C" a subzone.

```julia
julia> Z = Zone('D', "16.2C")
Zone('D', "16.2C")

julia> convert(String, Z)
"D16.2C"

julia> Zone("D16.2C")
Zone('D', "6.2C")
```

Other examples: "D-16.2C", "D-17.1C".
"""
struct Zone
    high::Char
    subzone::String
end

# assume that Z is of the correct form eg. "D-16.2C"
Zone(Z::String) = Zone(Z[1], Z[3:end])

Base.convert(::Type{String}, Z::Zone) = Z.high * Z.subzone;

# ============================
# Structs to define a package
# ============================

struct TravelTimes{MT}
    route_id::String
    stop_names::Vector{String}
    travel_times::MT
end

name(t::TravelTimes) = t.route_id
stop_names(t::TravelTimes) = t.stop_names
times_matrix(t::TravelTimes) = t.travel_times

# dimensions of a package
struct Dimensions
    depth::Float64
    height::Float64
    width::Float64
end

# expected time window for delivery of a package
struct TimeWindow
    start_date::Date
    start_time::Time
    end_date::Date
    end_time::Time
end

function is_defined(x::TimeWindow)
    got_start_date = x.start_date != Date("")
    got_end_date = x.end_date != Date("")
    got_start_time = x.start_time != Time("")
    got_end_time = x.end_time != Time("")

    return got_start_date && got_end_date && got_start_time && got_end_time
end

#=
abstract type AbstractScanStatus end

struct Delivered end <: AbstractScanStatus
struct Rejected end <: AbstractScanStatus
struct DeliveryAttempted end <: AbstractScanStatus
struct Unknown end <: AbstractScanStatus
=#

abstract type AbstractPackage end

struct Package <: AbstractPackage
    hex_hash::String
    time_window::TimeWindow
    planned_service::Float64 # UInt32 tiempo en entregar el paquete una vez llegado a la parada
    dimensions::Dimensions
    scan_status::Int64  # 0 = unknown, 1 = rejected, 2 = attempted, 3 = delivered
end

# by default the scan status is unknown (= 0)
Package(hex_hash::String, time_window::TimeWindow, planned_service::Float64, dimensions::Dimensions) = Package(hex_hash, time_window, planned_service, dimensions, 0)

name(P::Package) = P.hex_hash
time_window(P::Package) = P.time_window
has_time_window(P::Package) =  is_defined(P.time_window)

planned_service(P::Package) = P.planned_service
dimensions(P::Package) = P.dimensions

is_unknown(P::Package) = P.scan_status == 0
is_rejected(P::Package) = P.scan_status == 1
is_attempted(P::Package) = P.scan_status == 2
is_delivered(P::Package) = P.scan_status == 3

function volume(P::Package)
    d = P.dimensions
    return d.depth * d.height * d.width
end

# ============================
# Structs to define a stop
# ============================

abstract type AbstractStop end

# parada que solo contiene el nombre
struct BareStop <: AbstractStop
    id::String
end

name(x::BareStop) = x.id

"""
    Stop

Struct to store the information relative to a given stop.
"""
struct Stop <: AbstractStop
    id::String
    lat::Float64
    lng::Float64
    stop_type::Bool # 0 if Station, 1 if Dropoff
    zone_id::Zone
end

# constructor overloads
Stop(id::String, lat::Float64, lng::Float64, stop_type::Bool, zone_id::String) = Stop(id, lat, lng, stop_type, Zone(zone_id))
Stop(id::String, lat::Float64, lng::Float64, stop_type::Int, zone_id::String) = Stop(id, lat, lng, Bool(stop_type), Zone(zone_id))

# 0 if Station, 1 if Dropoff
is_station(x::Stop) = x.stop_type == 0
is_dropoff(x::Stop) = x.stop_type == 1;

latitude(x::Stop) = x.lat
longitude(x::Stop) = x.lng

zone(x::Stop) = convert(String, x.zone_id)
name(S::Stop) = S.id

function distance(X::AbstractStop, Y::AbstractStop)
    latX = latitude(X)
    latY = latitude(Y)
    lonX = longitude(X)
    lonY = longitude(Y)
    out = (latX - latY)^2 + (lonX - lonY)^2
    return sqrt(out)
end

# parada asociada a un vector de paquetes para entregar
struct Task{ST<:AbstractStop, PT<:AbstractPackage} <: AbstractStop
    stop::ST
    packages::Vector{PT}
end

packages(x::Task) = x.packages

# 0 if Station, 1 if Dropoff
is_station(x::Task) = is_station(x.stop)
is_dropoff(x::Task) = is_dropoff(x.stop)

latitude(x::Task) = latitude(x.stop)
longitude(x::Task) = longitude(x.stop)

zone(x::Task) = zone(x.stop)
name(x::Task) = name(x.stop)

# ================================
# Structs to define routes
# ================================

"""
    AbstractRoute

Abstract supertype for all route types.

### Notes

Requires to define the following fields:

`stops` -- vector of stops for the route
"""
abstract type AbstractRoute end

function station(R::AbstractRoute)
    S = R.stops
    idx = findall(is_station, S)
    if length(idx) == 0
        throw(ArgumentError("no station found"))
    elseif length(idx) > 1
        throw(ArgumentError("got more than one station"))
    end
    return S[idx[1]]
end

function stops(R::AbstractRoute; remove_station=false)
    S = R.stops
    if remove_station
        return filter(!is_station, S)
    else
        return S
    end
end

function zones(R::AbstractRoute; remove_station=false)
    Z = unique(zone.(R.stops))
    if remove_station
        return filter(x -> x != "0", Z)
    else
        return Z
    end
end

"""
    split_by_zone(R::AbstractRoute)

Return the set of stop indices associated to each zone.
"""
function split_by_zone(R::AbstractRoute; remove_station=true)

    S = stops(R, remove_station=remove_station)
    num_stops = length(S)

    zone_names = unique(zone.(S))
    num_zones = length(zone_names)

    # preallocate output vectors
    zone_idx = [Vector{Int}() for _ in 1:num_zones]
    stop_names = [Vector{String}() for _ in 1:num_zones]

    for i in 1:num_stops
        idx = findfirst(x -> x == zone(S[i]), zone_names)
        push!(zone_idx[idx], i)
        push!(stop_names[idx], name(S[i]))
    end
    return zone_idx, zone_names, stop_names
end

"""
    is_sorted_by_zone(R::AbstractRoute)

Return whether the stops in this route are sorted according to the zone.
"""
function is_sorted_by_zone(R::AbstractRoute)
    # TODO
end

"""
    Route

Struct to store the information relative to a given route.
"""
struct Route{ST<:AbstractStop} <: AbstractRoute
    hex_hash::String
    station_code::String
    date::Date
    departure_time::Time
    capacity::Float64 # UInt32
    stops::Vector{ST}
    score::Int64 # = 0 if the score has not been assigned, 1 = low, 2 = medium, 3 = high
end

# default constructor without score assigned
Route(hex_hash, station_code, date, departure_time, capacity, stops) = Route(hex_hash, station_code, date, departure_time, capacity, stops, 0)
name(R::Route) = R.hex_hash
capacity(R::Route) = R.capacity
score(R::Route) = R.score

is_score_unknown(R::Route) = R.score == 0
is_score_low(R::Route) = R.score == 1
is_score_medium(R::Route) = R.score == 2
is_score_high(R::Route) = R.score == 3

struct BareRoute{ST<:AbstractStop} <: AbstractRoute
    hex_hash::String
    stops::Vector{ST}
    score::Int64 # = 0 if the score has not been assigned, 1 = low, 2 = medium, 3 = high
end

BareRoute(hex_hash, stops) = BareRoute(hex_hash, stops, 0)

name(R::BareRoute) = R.hex_hash
score(R::BareRoute) = R.score


text:
function resample_by_1D(arr::AbstractArray{<:Complex, N}, new_size; normalize=true) where N
    return resample_by_1D_FT!(copy(arr), new_size, normalize=normalize)
end

function resample_by_1D(arr::AbstractArray{<:Real, N}, new_size; normalize=true) where N
    return real(resample_by_1D_FT!(Complex.(arr), new_size, normalize=normalize))
end

function resample_by_1D_FT!(arr::AbstractArray{<:Complex, N}, new_size; normalize=true) where N
    correction_factor = 1 # √2
    initial_length = length(arr)
    for d = 1:N
        ns = new_size[d]
        s = size(arr, d)
        # nothing to do in that dimension
        if ns == s
            continue
        end
        # go to fourier space
        arr = ffts!(arr, d)
        if ns > s
            out = zeros(eltype(arr), Base.setindex(size(arr), ns, d))
            center_set!(out, arr)
            # in the even case we need to fix hermitian property
            if iseven(s)
                l, r = get_indices_around_center(ns, s)
                inds_left = NDTools.slice_indices(axes(out), d, l)
                out[inds_left...] .*= 0.5 * correction_factor
                inds_right = NDTools.slice_indices(axes(out), d, r+1)
                out[inds_right...] .= out[inds_left...]
            end
            arr = iffts(out, d)
        else ns < s
            # extract a new view array with new size in that dimension
            arr_v = center_extract(arr, Base.setindex(size(arr), ns, d))
            # in the even case, we need to add the highest slice
            if iseven(ns)
                l, r = get_indices_around_center(s, ns)
                inds_left = NDTools.slice_indices(axes(arr_v), d, 1)
                inds_right = NDTools.slice_indices(axes(arr), d, r+1)
                arr_v[inds_left...] += arr[inds_right...]
                arr_v[inds_left...] ./= correction_factor
            end
            #overwrite old arr handle
            arr = iffts(arr_v, d)
        end
    end
    # normalize that values scale accordingly
    # this violates energy!
    if normalize
        arr .*= length(arr) ./ initial_length
    end
    return arr
end


text:
using JuMP

function conversion(model::Model)
    if num_variables(model) != 2
        error("Only models with two variables can be converted to GeoGebra")
    end

    x, y = all_variables(model)

    coefs = Vector{Tuple{Float64, Float64, Float64, Bool}}()
    for cst in all_constraints(model, AffExpr, MOI.GreaterThan{Float64})
        push!(coefs, (normalized_coefficient(cst, x), normalized_coefficient(cst, y), normalized_rhs(cst), true))
    end
    for cst in all_constraints(model, AffExpr, MOI.LessThan{Float64})
        push!(coefs, (normalized_coefficient(cst, x), normalized_coefficient(cst, y), normalized_rhs(cst), false))
    end

    #point_obj = (objective_function(model)[x], objective_function(model)[y])
    point_obj = (get(objective_function(model).terms, x, 0.0), get(objective_function(model).terms, y, 0.0))

    open("geogebra.xml", "w") do f
        write(f, """
<?xml version="1.0" encoding="utf-8"?>
<geogebra format="5.0" version="5.0.620.0" app="classic" platform="offline" id="F0B3856E-58E2-4C3F-B8ED-7C0AF991ACA8"  xsi:noNamespaceSchemaLocation="http://www.geogebra.org/apps/xsd/ggb.xsd" xmlns="" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" >
<gui>
	<window width="1848" height="1016" />
	<perspectives>
<perspective id="tmp">
	<panes>
		<pane location="" divider="0.18452380952380953" orientation="1" />
	</panes>
	<views>
		<view id="512" toolbar="0 | 1 501 5 19 , 67 | 2 15 45 18 , 7 37 | 514 3 9 , 13 44 , 47 | 16 51 | 551 550 11 ,  20 22 21 23 , 55 56 57 , 12 | 69 | 510 511 , 512 513 | 533 531 , 534 532 , 522 523 , 537 536 , 535 , 538 | 521 520 | 36 , 38 49 560 | 571 30 29 570 31 33 | 17 | 540 40 41 42 , 27 28 35 , 6 , 502" visible="false" inframe="false" stylebar="false" location="1,1,1" size="500" window="100,100,600,400" />
		<view id="4" toolbar="0 || 2020 , 2021 , 2022 || 2001 , 2003 , 2002 , 2004 , 2005 || 2040 , 2041 , 2042 , 2044 , 2043" visible="false" inframe="false" stylebar="false" location="1,1" size="300" window="100,100,600,400" />
		<view id="8" toolbar="1001 | 1002 | 1003  || 1005 | 1004 || 1006 | 1007 | 1010 || 1008 | 1009 || 6" visible="false" inframe="false" stylebar="false" location="1,3" size="300" window="100,100,600,400" />
		<view id="4097" visible="false" inframe="false" stylebar="true" location="1,1" size="344" window="100,100,700,550" />
		<view id="1" visible="true" inframe="false" stylebar="false" location="1" size="1499" window="100,100,600,400" />
		<view id="2" visible="true" inframe="false" stylebar="true" location="3" size="341" tab="ALGEBRA" window="100,100,600,400" />
		<view id="16" visible="false" inframe="false" stylebar="false" location="1" size="300" window="50,50,500,500" />
		<view id="32" visible="false" inframe="false" stylebar="true" location="1" size="300" window="50,50,500,500" />
		<view id="64" toolbar="0" visible="false" inframe="false" stylebar="false" location="1" size="480" window="50,50,500,500" />
		<view id="128" visible="false" inframe="false" stylebar="false" location="1" size="480" window="50,50,500,500" />
		<view id="70" toolbar="0 || 2020 || 2021 || 2022" visible="false" inframe="false" stylebar="true" location="1" size="900" window="50,50,500,500" />
	</views>
	<toolbar show="true" items="0 73 62 | 1 501 67 , 5 19 , 72 75 76 | 2 15 45 , 18 65 , 7 37 | 4 3 8 9 , 13 44 , 58 , 47 | 16 51 64 , 70 | 10 34 53 11 , 24  20 22 , 21 23 | 55 56 57 , 12 | 36 46 , 38 49  50 , 71  14  68 | 30 29 54 32 31 33 | 25 17 26 60 52 61 | 40 41 42 , 27 28 35 , 6" position="1" help="false" />
	<input show="true" cmd="true" top="algebra" />
	<dockBar show="false" east="false" />
</perspective>
	</perspectives>
	<labelingStyle  val="0"/>
	<font  size="16"/>
</gui>
<euclidianView>
	<viewNumber viewNo="1"/>
	<size  width="1499" height="963"/>
	<coordSystem xZero="302.2278178518016" yZero="553.9915207563241" scale="180.28883765485202" yscale="180.28883765485207"/>
	<evSettings axes="true" grid="true" gridIsBold="false" pointCapturing="3" rightAngleStyle="1" checkboxSize="26" gridType="3"/>
	<bgColor r="255" g="255" b="255"/>
	<axesColor r="0" g="0" b="0"/>
	<gridColor r="192" g="192" b="192"/>
	<lineStyle axes="1" grid="0"/>
	<axis id="0" show="true" label="" unitLabel="" tickStyle="1" showNumbers="true"/>
	<axis id="1" show="true" label="" unitLabel="" tickStyle="1" showNumbers="true"/>
</euclidianView>
<algebraView>
	<mode val="2"/>
	<auxiliary show="true"/>
	<collapsed val="0,3"/>
</algebraView>
<kernel>
	<continuous val="false"/>
	<usePathAndRegionParameters val="true"/>
	<decimals val="2"/>
	<angleUnit val="degree"/>
	<algebraStyle val="3" spreadsheet="0"/>
	<coordStyle val="0"/>
</kernel>
<tableview min="-2" max="2" step="1"/>
<scripting blocked="false" disabled="false"/>
<construction title="" author="" date="">
<element type="numeric" label="v">
	<value val="0"/>
	<show object="false" label="true"/>
	<objColor r="0" g="0" b="255" alpha="0.10000000149011612"/>
	<layer val="2"/>
	<labelMode val="1"/>
	<slider min="-5" max="5" absoluteScreenLocation="true" width="200" x="545" y="515" fixed="false" horizontal="true" showAlgebra="true"/>
	<lineStyle thickness="10" type="0" typeHidden="1"/>
	<animation speed="1" type="0" playing="false"/>
</element>
""")

    # Constraints
    allcstexp = ""
    for i in 1:length(coefs)
        const_expression = ""
        const_expression_eq = ""
        if !(coefs[i][1] in [0,1])
            const_expression *= "("
            if coefs[i][1] < 0
                const_expression *= "-"
            end
            if coefs[i][1] != -1
                const_expression *= "$(abs(coefs[i][1])) * x"
            else
                const_expression *= "x"
            end

            const_expression *= ") "
        else
            if coefs[i][1] == 1
                const_expression *= "x "
            end
        end

        if coefs[i][1] != 0 && coefs[i][2] != 0
            const_expression *= "+ "
        end

        if !(coefs[i][2] in [0,1])
            const_expression *= "("
            if coefs[i][2] < 0
                const_expression *= "-"
            end
            if coefs[i][2] != -1
                const_expression *= "$(abs(coefs[i][2])) * y"
            else
                const_expression *= "y"
            end

            const_expression *= ") "
        else
            if coefs[i][2] == 1
                const_expression *= "y "
            end
        end

        const_expression_eq = const_expression

        if coefs[i][4]
            const_expression *= "≥ "
        else
            const_expression *= "≤ "
        end

        const_expression *= "$(coefs[i][3])"
        const_expression_eq *= "= $(coefs[i][3])"

        allcstexp *= " ∧ " * const_expression

        write(f, """
<expression label="cst$i" exp="cst$i(x, y) = $const_expression"/>
<element type="functionnvar" label="cst$i">
	<show object="false" label="false" ev="4"/>
	<objColor r="0" g="0" b="0" alpha="0" fillType="4" hatchAngle="45" hatchDistance="20" inverseFill="true"/>
	<layer val="1"/>
	<labelMode val="0"/>
	<animation step="0.1" speed="1" type="0" playing="false"/>
	<lineStyle thickness="3" type="0" typeHidden="1"/>
</element>$(coefs[i][1] == 0 ? "\n<expression label=\"c$(i)\" exp=\"$(replace(replace(const_expression, "≤" => "="), "≥" => "="))\" type=\"line\"/>" : "")
<element type="line" label="c$(i)">
	<show object="false" label="true" ev="4"/>
	<objColor r="0" g="0" b="0" alpha="0"/>
	<layer val="1"/>
	<labelMode val="0"/>
	<fixed val="true"/>
	<coords x="$(coefs[i][1])" y="$(coefs[i][2])" z="$(-coefs[i][3])"/>
	<lineStyle thickness="3" type="0" typeHidden="1"/>
	<eqnStyle style="implicit"/>
</element>
""")
    end

    # Var bounds
    bounds = ""
    if has_lower_bound(x)
        bounds *= " ∧ x ≥ $(lower_bound(x))"
    end
    if has_upper_bound(x)
        bounds *= " ∧ x ≤ $(upper_bound(x))"
    end
    if has_lower_bound(y)
        bounds *= " ∧ y ≥ $(lower_bound(y))"
    end
    if has_upper_bound(y)
        bounds *= " ∧ y ≤ $(upper_bound(y))"
    end

    if bounds != ""
        write(f, """
<expression label="var" exp="var(x, y) = $(bounds[(sizeof(" ∧ ")+1):end])"/>
<element type="functionnvar" label="var">
	<show object="false" label="false" ev="4"/>
	<objColor r="0" g="0" b="0" alpha="0" fillType="4" hatchAngle="45" hatchDistance="20" inverseFill="true"/>
	<layer val="1"/>
	<labelMode val="0"/>
	<animation step="0.1" speed="1" type="0" playing="false"/>
	<lineStyle thickness="3" type="0" typeHidden="1"/>
</element>
""")
    end

    allcstexp *= bounds
    write(f, """
<expression label="cst" exp="cst(x, y) = $(allcstexp[(sizeof(" ∧ ")+1):end])"/>
<element type="functionnvar" label="cst">
	<show object="false" label="false" ev="4"/>
	<objColor r="204" g="153" b="255" alpha="0.28"/>
	<layer val="1"/>
	<labelMode val="0"/>
	<animation step="0.1" speed="1" type="0" playing="false"/>
	<lineStyle thickness="0" type="0" typeHidden="1"/>
</element>
""")

    # Objective
    write(f, """
<element type="point" label="OBJ">
	<show object="false" label="true" ev="4"/>
	<objColor r="125" g="125" b="255" alpha="0"/>
	<layer val="2"/>
	<labelMode val="0"/>
	<animation step="0.1" speed="1" type="1" playing="false"/>
	<coords x="$(point_obj[1])" y="$(point_obj[2])" z="1"/>
	<pointSize val="5"/>
	<pointStyle val="0"/>
</element>
<expression label="ZERO" exp="(0, 0)" type="point"/>
<element type="point" label="ZERO">
	<show object="false" label="true" ev="4"/>
	<objColor r="125" g="125" b="255" alpha="0"/>
	<layer val="0"/>
	<labelMode val="0"/>
	<animation step="0.1" speed="1" type="1" playing="false"/>
	<coords x="0" y="0" z="1"/>
	<pointSize val="5"/>
	<pointStyle val="0"/>
</element>
<command name="Ray">
	<input a0="ZERO" a1="OBJ"/>
	<output a0="a"/>
</command>
<element type="ray" label="a">
	<show object="false" label="false" ev="4"/>
	<objColor r="0" g="0" b="255" alpha="0.75"/>
	<layer val="2"/>
	<labelMode val="0"/>
	<coords x="$(-point_obj[2])" y="$(point_obj[1])" z="0"/>
	<lineStyle thickness="2" type="0" typeHidden="1"/>
	<outlyingIntersections val="false"/>
	<keepTypeOnTransform val="true"/>
</element>
<expression label="A" exp="(($(point_obj[1]) * v), ($(point_obj[2]) * v))" type="point" />
<element type="point" label="A">
	<show object="false" label="true" ev="4"/>
	<objColor r="125" g="125" b="255" alpha="0"/>
	<layer val="0"/>
	<labelMode val="0"/>
	<coords x="0" y="0" z="1"/>
	<pointSize val="5"/>
	<pointStyle val="0"/>
</element>
<command name="OrthogonalLine">
	<input a0="A" a1="a"/>
	<output a0="f"/>
</command>
<element type="line" label="f">
	<show object="false" label="false" ev="4"/>
	<objColor r="0" g="0" b="255" alpha="0"/>
	<layer val="2"/>
	<labelMode val="0"/>
	<coords x="$(-point_obj[1])" y="$(-point_obj[2])" z="0"/>
	<lineStyle thickness="5" type="0" typeHidden="1" opacity="178"/>
	<eqnStyle style="implicit"/>
</element>
""")

    write(f, """
</construction>
</geogebra>
""")
    end

    # Compress geogebra.xml
    run(`zip geogebra.zip geogebra.xml`)

    # Delete geogebra.xml
    rm("geogebra.xml")

    # Convert into geogebra.ggb
    mv("geogebra.zip", "geogebra.ggb", force=true)

    return nothing
end


text:
function  count_solution(p::Int)
	ans=0

	for a=1:p
		for b=a: floor(Int,(p-a)/2)
			c=p-a-b
			if a^2+b^2==c^2
				ans+=1
			end
		end
	end

	return ans

end

function compute()
	ans=-1
	index=0;
	for i=1:1000
		cnt=count_solution(i);
		if ans<cnt
			ans=cnt
			index=i
		end
	end

	return index
end

println(compute())


text:
using Random
#Faber Schauder Functions structure
include("faber.jl")

#structures for models and ZigZag sampling schemes
include("types.jl")

#functions for Faber Schauder expansion by finite elements representation
include("fs_expansion.jl")

#ZigZag sampler
include("zz_sampler.jl")

#plotting zigzag sampler
include("plotting.jl")

max_stars_repo_name: SebaGraz/ZZDiffusionBridge

max_issues_repo_name: SebaGraz/ZZDiffusionBridge

max_forks_repo_name: SebaGraz/ZZDiffusionBridge


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* Author: Rafael de Pelegrini Soares
#* $Id$
#*--------------------------------------------------------------------
type Capacitor
	Capacitor()=begin
		new(
			capacitance(),
			charge(),
			wire(Dict{Symbol,Any}(
				:Brief=>"Inlet",
				:PosX=>0.3978,
				:PosY=>0,
				:Symbol=>"_{in}"
			)),
			wire(Dict{Symbol,Any}(
				:Brief=>"Outlet",
				:PosX=>0.3965,
				:PosY=>1,
				:Symbol=>"_{out}"
			)),
			[
				:(diff(q) = inlet.i),
				:(inlet.V - outlet.V = (1/C) * q),
				:(outlet.i = inlet.i),
			],
			[
				"","","",
			],
			[:C,],
			[:q,:inlet,:outlet,]
		)
	end
	C::capacitance
	q::charge
	inlet::wire
	outlet::wire
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Capacitor
function setEquationFlow(in::Capacitor)
	addEquation(1)
	addEquation(2)
	addEquation(3)
end
function atributes(in::Capacitor,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/Capacitor"
	fields[:Brief]="Electrical Capacitor."
	drive!(fields,_)
	return fields
end
Capacitor(_::Dict{Symbol,Any})=begin
	newModel=Capacitor()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
using  Dates
using  Base.Threads
using  Printf
    function __GenFilename(typeStr::String,td::DateTime )
        curpath = abspath( pwd()* "/../") #"/../Jl-Aux/"
        minuts  =0;#Dates.minute(td)÷30 *30
        windStatuFileName = @sprintf("%s/log/%s-%s-%02d-%02d-%02d-%02d.lst",curpath,Main.machineName,typeStr,Dates.month(td), Dates.day(td),Dates.hour(td),Dates.minute(td) )
        #println(windStatuFileName);println();
        return windStatuFileName,Dates.hour(td)
    end

    mutable struct __logfile
        fileLock::SpinLock
        fileObj
        fileType::String
        switchTime::Int
        function __logfile( fileType1::String )
            # global filesBuffer
            fileLock = SpinLock()
            td1 = Dates.now()
            name,hour1 = __GenFilename(fileType1,td1)
            fileObj = open( name,"w")
            # fileType = fileType1
            switchTime = hour1
            new(fileLock,fileObj,fileType1,switchTime  )
        end
    end
    function __switchFile( logfileObj , td1 )
        close( logfileObj.fileObj)
        name,hour1 = __GenFilename(logfileObj.fileType,td1)
        logfileObj.fileObj = open( name,"w")
        # fileType = fileType1
        logfileObj.switchTime = Dates.hour(td1)
    end

    logfileObj = nothing

    function Init( fileType::String)
        global  logfileObj
        logfileObj= __logfile( fileType)
    end

    function LogEvent(info1)
        # return ;########### test only
        td1 = Dates.now()
        lock(logfileObj.fileLock)
            if Dates.hour( td1 ) != logfileObj.switchTime
                __switchFile( logfileObj , td1 )
            end
                # f1 = open( windStatuFileName,"a+")
            write( logfileObj.fileObj,info1 )
            flush(logfileObj.fileObj )
                # close( f1 )
            # @warn  info1
        unlock(logfileObj.fileLock)

    end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


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

"""
    KNearestSearch(domain, k; metric=Euclidean())

A method for searching `k` nearest neighbors in `domain`
according to `metric`.
"""
struct KNearestSearch{D,T} <: BoundedNeighborSearchMethod
  # input fields
  domain::D
  k::Int

  # state fields
  tree::T
end

function KNearestSearch(domain::D, k::Int; metric=Euclidean()) where {D}
  xs = [coordinates(centroid(domain, i)) for i in 1:nelements(domain)]
  tree = metric isa MinkowskiMetric ? KDTree(xs, metric) : BallTree(xs, metric)
  KNearestSearch{D,typeof(tree)}(domain, k, tree)
end

maxneighbors(method::KNearestSearch) = method.k

function search!(neighbors, pₒ::Point, method::KNearestSearch; mask=nothing)
  k       = method.k
  inds, _ = knn(method.tree, coordinates(pₒ), k, true)

  if mask ≠ nothing
    nneigh = 0
    @inbounds for i in 1:k
      if mask[inds[i]]
        nneigh += 1
        neighbors[nneigh] = inds[i]
      end
    end
  else
    nneigh = k
    @inbounds for i in 1:k
      neighbors[i] = inds[i]
    end
  end

  nneigh
end


text:
fig = Figure(fontsize = 28, backgroundcolor = :transparent);
# fig = Figure(font = assetpath("fonts", fontpath), fontsize=24);
ax = Axis(fig[1, 1]);

b = 0.05 # [0, 0.01, 0.05, 0.1, 0.15]
d_f = 2 # [2,10,20,30,40,50,75,100]

dft = filter([:driving_force, :b] => (x, y) -> x == d_f && y == b, df)

idxs = sortperm(dft.:S_ext)
xs = dft.:S_ext[idxs]

prot_b = dft.:b[idxs]
prot_tw = dft.:b[idxs] + dft.:tw[idxs]
prot_ts = dft.:b[idxs] + dft.:tw[idxs] + dft.:ts[idxs]
prot_l = dft.:b[idxs] + dft.:tw[idxs] + dft.:ts[idxs] + dft.:l[idxs]
prot_r = dft.:b[idxs] + dft.:tw[idxs] + dft.:ts[idxs] + dft.:l[idxs] + dft.:r[idxs]
prot_a = dft.:b[idxs] + dft.:tw[idxs] + dft.:ts[idxs] + dft.:l[idxs] + dft.:r + dft.:a[idxs]
prot_c =
    dft.:b[idxs] +
    dft.:tw[idxs] +
    dft.:ts[idxs] +
    dft.:l[idxs] +
    dft.:r +
    dft.:a[idxs] +
    dft.:c[idxs]

band!(
    ax,
    xs,
    zeros(length(prot_b)),
    prot_b ./ prot_c,
    color = colscheme[1],
    label = "Transporter",
)
band!(
    ax,
    xs,
    prot_b ./ prot_c,
    prot_tw ./ prot_c,
    color = colscheme[2],
    label = "Metabolic enzyme",
)
band!(
    ax,
    xs,
    prot_tw ./ prot_c,
    prot_ts ./ prot_c,
    color = colscheme[3],
    label = "Ribosomes",
)
band!(
    ax,
    xs,
    prot_ts ./ prot_c,
    prot_l ./ prot_c,
    color = colscheme[4],
    label = "Lipid synthesis enzyme",
)

#hidexdecorations!(ax, ticks = false, ticklabels = false, label=false)
#hideydecorations!(ax, ticks = false, ticklabels = false, label=false)
ax.xlabel = "External substrate concentration"
ax.ylabel = "Proteome fraction"
fig[1, 2] = Legend(fig, ax, "Protein", unique = true, framevisible = false)

fig


text:
function recursiveCycle(param::MGparam{VAL,IND},b::Array{VAL},x::Array{VAL},level::Int64) where {VAL,IND}

# println(string("Starting level ", level));

gmresTol = 1e-5;

n = size(b,1);
nrhs = size(b,2);
numCores = param.numCores;
As = param.As;
nlevels = length(As);

if level==nlevels # This actually does not need to happen unless one level only is used (i.e. exact solver).
	r = param.memCycle[level].r;
	r[:] = b;
	x = solveCoarsest(param,r,x);
	return x;
end

Ps = param.Ps;
Rs = param.Rs;
AT = As[level];

oneType = one(VAL);
zeroType = zero(VAL);
r = param.memCycle[level].r;

r[:] = b;
if norm(x)>0.0
   	SpMatMul(-oneType,AT,x,oneType,r,numCores)#  r -= A'*x;
end
D = param.relaxPrecs[level];
MM = identity;
Afun = identity;
if param.relaxType=="Jac-GMRES"
	y = param.memRelax[level].v_prec;
	MM(xx::Array{VAL}) = (return y .= D.*xx);
end

PT = Ps[level];
RT = Rs[level];

npresmth  = param.relaxPre(level);
npostsmth = param.relaxPost(level);

if param.relaxType=="Jac-GMRES"
    Afun = getAfun(AT,param.memRelax[level].Az,numCores)
	x = FGMRES_relaxation(Afun,r,x,npresmth,MM,gmresTol,false,numCores,param.memRelax[level])[1];
elseif param.relaxType == "VankaFaces"
	x = RelaxVankaFacesColor(AT,x,b,D,npresmth,numCores,param.Meshes[level],param.transferOperatorType=="SystemsFacesMixedLinear",FULL_VANKA);
elseif param.relaxType == "EconVankaFaces"
	x = RelaxVankaFacesColor(AT,x,b,D,npresmth,numCores,param.Meshes[level],param.transferOperatorType=="SystemsFacesMixedLinear",ECON_VANKA);
elseif param.relaxType=="hybridVankaFacesKaczmarz"
	x = RelaxHybridVanka(param.relaxPrecs[level], AT,x,b,npresmth,param.relaxPrecs[level].numCores,param.Meshes[level],
								param.transferOperatorType=="SystemsFacesMixedLinear",KACMARZ_VANKA);
else
	x = relax(AT,r,x,b,D,npresmth,numCores);
end

SpMatMul(-oneType,AT,x,zeroType,r,numCores); #  r = -A'*x;

addVectors(oneType,b,r); # r = r + b;

xc = param.memCycle[level+1].x;
xc[:] .= 0.0;
bc = param.memCycle[level+1].b;
bc = SpMatMul(RT,r,bc,numCores);
if level==nlevels-1
	# println("solving coarsest");
	xc = solveCoarsest(param,bc,xc);
else
	Ac = As[level+1];
    if param.cycleType == 'K'
		yzK = param.memKcycle[level].v_prec;
		AfunK = getAfun(Ac,param.memKcycle[level].Az,numCores);
		MMG(x) = (yzK .= 0.0; recursiveCycle(param,x,yzK,level+1)); # x does not change...
		xc = FGMRES_relaxation(AfunK,bc,xc,2,MMG,gmresTol,false,numCores,param.memKcycle[level])[1];
    else
		xc = recursiveCycle(param,bc,xc,level+1);# xc changes, bc does not change
		if param.cycleType=='W'
            xc = recursiveCycle(param,bc,xc,level+1);
		elseif param.cycleType=='F'
            param.cycleType='V';
            xc = recursiveCycle(param,bc,xc,level+1);
			param.cycleType='F';
        end
    end
end

x = SpMatMul(oneType,PT,xc,oneType,x,numCores); # x += PT'*xc;

r[:] = b;
SpMatMul(-oneType,AT,x,oneType,r,numCores); #  r -= A'*x;
if param.relaxType=="Jac-GMRES"
	Afun = getAfun(AT,param.memRelax[level].Az,numCores)
	x = FGMRES_relaxation(Afun,r,x,npostsmth,MM,gmresTol,false,numCores,param.memRelax[level])[1];
elseif param.relaxType == "VankaFaces"
	x = RelaxVankaFacesColor(AT,x,b,D,npostsmth,numCores,param.Meshes[level],param.transferOperatorType=="SystemsFacesMixedLinear",FULL_VANKA);
elseif param.relaxType == "EconVankaFaces"
	x = RelaxVankaFacesColor(AT,x,b,D,npostsmth,numCores,param.Meshes[level],param.transferOperatorType=="SystemsFacesMixedLinear", ECON_VANKA);
elseif param.relaxType=="hybridVankaFacesKaczmarz"
	x = RelaxHybridVanka(param.relaxPrecs[level], AT,x,b,npostsmth,param.relaxPrecs[level].numCores,param.Meshes[level],
								param.transferOperatorType=="SystemsFacesMixedLinear",KACMARZ_VANKA);
else
	x = relax(AT,r,x,b,D,npostsmth,numCores);
end

return x
end

function relax(AT::SparseMatrixCSC{VAL,IND},r::Array{VAL},x::Array{VAL},b::Array{VAL},d::Vector{VAL},numit::Int64,numCores::Int64) where {VAL,IND}
# x is actually the output... x and is being changed in the iterations.
# r does not end up accurate becasue we do not update it in the last iteration.
oneType = one(eltype(r));
zeroType = zero(eltype(r));
# nr0 = vecnorm(r);
for i=1:numit-1
	x .+= d.*r;
	SpMatMul(-oneType,AT,x,zeroType,r,numCores) # r = -A'*x
	addVectors(oneType,b,r); # r = r + b;
	# println("Reduced: ", vecnorm(r)/nr0)
end
x .+= d.*r;
return x
end

function solveCoarsest(param::MGparam{VAL,IND},b::Array{VAL},x::Array{VAL},doTranspose::Int64=0) where {VAL,IND}

if isa(param.LU,DomainDecompositionParam)
	AT = param.As[end];
	(x,param.LU) = solveDDSerial(AT,b,x,param.LU,1,doTranspose);
	return x;
elseif isa(param.LU,AbstractSolver)
	AT = param.As[end];
	(x,param.LU) = solveLinearSystem!(AT,b,x,param.LU);
	return x;
end

if param.coarseSolveType == "MUMPS"
	applyMUMPS(param.LU,b,xt,param.doTranspose);
elseif param.coarseSolveType == "GMRES"
	AT = param.As[end];
	maxIter = 1;
	tol = 0.01;
	out= -2;
	Afun = getAfun(AT,zeros(eltype(b),size(b)),param.numCores);
	d = param.LU;
	y = zeros(eltype(b),size(b));
	M2(xx::Array{VAL}) = (y .= d.*xx; return y);
	if size(b,2)==1
		b = vec(b);
		x[:] .= 0.0;
		x, flag,rnorm,iter = KrylovMethods.fgmres(Afun,b,10,tol = tol,maxIter = 1,M = M2,out=out,x=x);
	else
		x[:] .= 0.0;
		x, flag,rnorm,iter = KrylovMethods.blockFGMRES(Afun,b,10,tol = tol,maxIter = 1,M = M2,out=out,X=x);
	end
elseif param.coarseSolveType == "VankaFaces"
	AT = param.As[end];
	inner = 20;
	Afun = getAfun(AT,zeros(eltype(b),size(b)),param.numCores);
	out= -2;
	x = KrylovMethods.fgmres(Afun,b,inner;tol = 0.8,maxIter = 5,M = param.LU,x = x,flexible = true,out = out)[1];
	# println("Coarsest norm is: ",norm(Afun(y) - b)/norm(b));
else
	z = param.LU\b;
	x[:] = z;
end
return x;
end


text:
# The file test.mseed comes from an older IRIS libmseed, found by anowacki
# It has a more complicated structure than the test.mseed file in more recent
# versions of libmseed, which reads with no issues
printstyled("  mini-SEED\n", color=:light_green)

printstyled("    sample rate\n", color=:light_green)

# Tests from SEED manual v2.4 page 110
import SeisIO.SEED.update_dt!
r1 = [33, 330, 3306, -60, 1, -10, -1]
r2 = [10, 1, -10, 1, -10, 1, -10]
fs = [330.0, 330.0, 330.6, 1.0/60.0, 0.1, 0.1, 0.1]

for i = 1:length(r1)
  BUF.r1 = r1[i]
  BUF.r2 = r2[i]
  update_dt!(BUF)
  @test isapprox(1.0/BUF.dt, fs[i])
end

test_sac_file     = string(path, "/SampleFiles/SAC/test_le.sac")
test_mseed_file   = string(path, "/SampleFiles/SEED/test.mseed")
test_mseed_pat    = string(path, "/SampleFiles/SEED/t*.mseed")
mseed_vals_file   = string(path, "/SampleFiles/SEED/test_mseed_vals.txt")

printstyled("    file read\n", color=:light_green)

@test_throws ErrorException verified_read_data("mseed", test_sac_file)

S = verified_read_data("mseed", test_mseed_file, v=0, strict=false)
@test isequal(S.id[1], "NL.HGN.00.BHZ")
@test ≈(S.fs[1], 40.0)
@test ≈(S.gain[1], 1.0)
@test isequal(string(u2d(S.t[1][1,2]*μs)), "2003-05-29T02:13:22.043")
@test ≈(S.x[1][1:5], [ 2787, 2776, 2774, 2780, 2783 ])

# mseed with mmap
printstyled("    file read with mmap\n", color=:light_green)
Sm = read_data("mseed", test_mseed_file, v=0, memmap=true)
@test Sm == S

# Test breaks if memory-resident SeisIOBuf structure SEED is not reset
S1 = verified_read_data("mseed", test_mseed_file, v=0, strict=false)
if Sys.iswindows() == false
  S2 = verified_read_data("mseed", test_mseed_pat, v=0, strict=false)
  @test S2.src[1] == abspath(test_mseed_pat)
  S2.src = S1.src
  @test S == S1 == S2
end

mseed_vals = readdlm(mseed_vals_file, ',', comments=true, comment_char='#')
seedvals = Dict{String,Any}()
ntest = size(mseed_vals,1)
for i = 1:ntest
  seedvals[mseed_vals[i,1]] = Float32.(mseed_vals[i, 2:end])
end

if safe_isdir(path*"/SampleFiles/Restricted")
  printstyled("    mseed test files (time gaps, blockette types)\n", color=:light_green)
  redirect_stdout(out) do
    seed_support()
    files = ls(path*"/SampleFiles/Restricted/*mseed")
    for f in files
      println(stdout, "attempting to read ", f)
      S = SeisData()
      ae = any([occursin(i, f) for i in ("blkt2000", "detection.record", "text-encoded", "timing.500s")])
      if ae
        verified_read_data!(S, "mseed", f, v=3, allow_empty=true, strict=false)
      else
        verified_read_data!(S, "mseed", f, v=2, strict=false)
      end

      # Test that our encoders return the expected values
      (tmp, fname) = splitdir(f)
      if haskey(seedvals, fname)
        x = get(seedvals, :fname, Float32[])
        nx = lastindex(x)
        y = getindex(getfield(S, :x), 1)[1:nx]
        @test isapprox(x,y)
      end

      if occursin("text-encoded", f)
        @test haskey(S.misc[1], "seed_ascii") == true
        str = split(S.misc[1]["seed_ascii"][1], "\n", keepempty=false)
        @test occursin("Quanterra Packet Baler Model 14 Restart.", str[1])

      elseif occursin("detection.record",f )
        ev_rec = get(S.misc[1], "seed_event", "no record")[1]
        @test ev_rec == "2004,7,28,20,28,6,185,80.0,0.39999998,18.0,dilatation,1,3,2,1,4,0,2,0,Z_SPWWSS"

      elseif occursin("SHW.UW", f)
        println("testing read accuracy of SHW.UW with ", f)
        @test size(S.t[1],1) >= 158
        @test size(S.t[2],1) >= 8
        @test string(u2d(S.t[1][1,2]*μs)) == "1980-03-22T20:45:18.349"
        @test isequal(S.id, String[ "UW.SHW..EHZ", "UW.SHW..SHZ" ])
        @test ≈(S.fs, Float64[104.085000, 52.038997])
        @test ≈(S.x[1][1:5], Float64[-68.0, -57.0, -71.0, -61.0, -52.0])
        fnames = ls(path*"/SampleFiles/Restricted/1980*SHZ.D.SAC")

        C = S[2]
        @test w_time(t_win(C.t, C.fs), C.fs) == C.t
        t = t_win(C.t, C.fs)[:,1]
        W = Array{DateTime,1}(undef, 0)
        for i=1:length(t)
          # push!(W, round(u2d(t[i]*μs), Second))
          push!(W, u2d(t[i]*μs))
        end

        Y = Array{DateTime,1}(undef,0)
        for f in fnames
          seis = verified_read_data("sac", f)[1]
          push!(Y, u2d(seis.t[1,2]*μs))
        end
        #[round(u2d(i*μs), Second) for i in t]
        # println("W = ", string.(W))
        # println("Y = ", string.(Y))
        Δ = [abs(.001*(W[i]-Y[i]).value)*C.fs for i=1:length(Y)]
        @test maximum(Δ) < 1.0
      else
        @test isempty(S) == false
      end
    end
  end
else
  printstyled("  extended SEED tests skipped. (files not found; is this Appveyor?)\n", color=:green)
end

# Tests for unparseable blockette and data format
mseed_out = "test2.mseed"
io = open(test_mseed_file, "r")
buf = read(io)
close(io)

buf[40] = 0x03                # 3 blockettes follow
buf[53] = 0x13                # 19 = Steim-3
buf[67:68] .= [0x00, 0x54]    # byte 84
buf[85:86] .= [0x01, 0x90]    # [400]
buf[87:88] .= [0x00, 0x00]    # next blockette at "byte 0" means we're done
write(mseed_out, buf)

S1 = read_data("mseed", test_mseed_file, v=0)[1]
S2 = read_data("mseed", mseed_out, v=0)[1]

# Check that skipping an unparseable data type on a new channel doesn't
# affect channel start time or data read-in
printstyled("    unparseable data\n", color=:light_green)
δx = length(S1.x)-length(S2.x)
@test div(S2.t[1,2]-S1.t[1,2], round(Int64, sμ/S1.fs)) == length(S1.x)-length(S2.x)
@test S1.x[δx+1:end] == S2.x

# Check that bytes skipped are accurately logged
printstyled("    unparseable blockettes\n", color=:light_green)
@test any([occursin("3968 bytes skipped", i) for i in S2.notes])

# Done
rm(mseed_out)


text:
module QRPNN

using LinearAlgebra
using Quaternions
using Random
rng = MersenneTwister(1234);

function identity(X,a)
    return X
end

function potential(X,L)
    return 1 ./ ((1 + 1.e-8 .- X) .^ L)
end

function high_order(X,q)
    return (1 .+ X) .^ q;
end

function exponential(X,alpha)
    return exp.(alpha*X)
end

function train(f,f_params,U)
    N = size(U)[1]
    C = f(Array{Float64}(real(U'*U))/N,f_params)
    return U*inv(C)
end

function hopfield(U, Up, xinput, it_max = 1.e3, verbose=true)
    ### Quaternionic Hopfield Neural Network
    Name = "Quaternionic Hopfield Neural Network"

    N = size(U,1)
    tau = 1.e-6

    # Initialization
    x = copy(xinput)
    xold = copy(x)
    it = 0
    Error = 1+tau
    while (Error>tau)&&(it<it_max)
        it = it+1

        # Compute the next state
        a = U*(Up*x)
        x = a ./ abs.(a)

        Error = norm(x-xold)
        xold = copy(x)
    end

    if verbose == true
        if it_max<=it
            println(Name," failed to converge in ",it_max," iterations.")
        end
    end

    return x
end

function main(f,f_params, U, V, xinput, it_max = 1.e3, verbose = true)
    (N,K) = size(U)
    tau = 1.e-6

    # Initialization
    x = copy(xinput)
    xold = copy(x)
    it = 0
    Error = 1+tau
    while (Error>tau)&&(it<it_max)
        it = it+1

        # Compute the weights
        w = f(Array{Float64}(real((U'*x)))/N,f_params);

        # Compute the next state
        a = V*w
        x = a ./ abs.(a)

        Error = norm(x-xold)
        xold = copy(x)
    end
    if verbose == true
        if it_max<=it
            println(split(string(f),".")[end]," QRPNN failed to converge in ",it_max," iterations.")
        end
    end
    return x
end

function CIFAR2Quat(x)
    tau = 1.e-4
    phi_ang = (-pi+tau) .+ 2*(pi-tau)*x[:,:,1]
    psi_ang = (-pi/4+tau) .+ (pi/2-2*tau)*x[:,:,2]
    theta_ang = (-pi/2+tau) .+ (pi-2*tau)*x[:,:,3]
    return Quaternion.(cos.(phi_ang),sin.(phi_ang),zeros(32,32),zeros(32,32)).*
        Quaternion.(cos.(psi_ang),zeros(32,32),zeros(32,32),sin.(psi_ang)).*
        Quaternion.(cos.(theta_ang),zeros(32,32),sin.(theta_ang),zeros(32,32));
end

function imnoise(img,noise_std = 0.1)
    return clamp.(img + noise_std*randn(rng, Float64, size(img)),0,1)
end

function Quat2CIFAR(x)
    tau = 1.e-4
    N = size(x,1)
    img_x = zeros(N,3)
    for i=1:length(x)
        q = x[i]/abs(x[i])

        # Take the components of q
        a = q.s;
        b = q.v1;
        c = q.v2;
        d = q.v3;

        # Rodriguez entries
        R11 = a.^2+b.^2-c.^2-d.^2;
        R13 = 2*(b.*d+a.*c);
        R22 = a.^2-b.^2+c.^2-d.^2;
        R32 = 2*(c.*d+a.*b);
        R12 = 2*(b.*c-a.*d);
        R23 = 2*(c.*d-a.*b);
        R33 = a.^2-b.^2-c.^2+d.^2;

        # Psi, Psi, and theta angles
        psiq = real(asin(-R12)/2);
        phiq   = real(atan(R32,R22)/2);
        thetaq = real(atan(R13,R11)/2);
        if abs(abs(psiq)-pi/4)<=tau
            phiq = real(atan(-R23,R33)/2);
            thetaq = 0;
        end

        a1   = cos(phiq).*cos(psiq).*cos(thetaq)+sin(phiq).*sin(psiq).*sin(thetaq);
        if (a1*a<0)
            if phiq>=0
                phiq = phiq - pi;
            else
                phiq = phiq + pi;
            end
        end
        R = (phiq+pi-tau)/(2*(pi-tau));
        G = (psiq+pi/4-tau)/((pi/2-2*tau));
        B = (thetaq+pi/2-tau)/((pi-2*tau));
        img_x[i,:] = [R,G,B]
    end
    return img_x
end

end


text:
MatrixOrFunction = Union{AbstractMatrix,Function}
VectorOrFunction = Union{AbstractVector,Function}
ArrayOrFunction = Union{AbstractArray,Function}
RealOrFunction = Union{Real,Function}

#Delay behavoiur: stochastic or deterministic / Zero order holders or continuous delay

struct Delay{T}
    τ::T
end
(τ::Delay{<:Real})(t) = τ.τ
(τ::Delay{<:Function})(t) = τ.τ(t)
Delay(τ::T1, holdT::T1) where {T1 <: Real} = Delay(t->τ + t % holdT)
Base.convert(::Type{Delay}, τ::RealOrFunction) = Delay(τ)
Base.convert(::Type{Vector{<:Delay}}, τs::Vector{<:RealOrFunction}) = Delay.(τs)

abstract type Coefficients{d} end
abstract type CoefficientMatrix{d} <: Coefficients{d} end
abstract type AdditiveVector{d} <: Coefficients{d} end

Base.size(cs::T) where T<:CoefficientMatrix{d} where d = (d, d)
Base.size(cs::T) where T<:AdditiveVector{d} where d = (d,)
Base.size(cs::T, ::Integer) where T<:Coefficients{d} where d = d
Base.length(cs::T) where T<:AdditiveVector{d} where d = d
Base.length(cs::T) where T<:CoefficientMatrix{d} where d = d^2

struct ProportionalMX{d,T} <: CoefficientMatrix{d}
    MX::T # matrix
end
ProportionalMX(mx::mxT) where mxT <: Function = ProportionalMX{size(mx(0.), 1),mxT}(mx)
ProportionalMX(mx::mxT) where mxT <: SMatrix = ProportionalMX{size(mx, 1),mxT}(mx)
ProportionalMX(mx::mxT) where mxT <: AbstractMatrix{<:Real} = ProportionalMX(SMatrix{size(mx)...}(mx))
(PMX::ProportionalMX{d,<:Function})(t) where d = PMX.MX(t)
(PMX::ProportionalMX{d,<:AbstractMatrix{<:Real}})(t) where d = PMX.MX
Base.convert(::Type{ProportionalMX}, mx::mT) where mT <: AbstractMatrix{<:Real} = ProportionalMX(mx)

struct DelayMX{d,dT,bT} <: CoefficientMatrix{d}
    τ::Delay{dT} # delay
    MX::bT # matrix
end

DelayMX(τ::RealOrFunction, MX::MatrixOrFunction) = DelayMX(Delay(τ), MX)
DelayMX(τ::Delay{dT}, MX::mxT) where {dT,mxT<:Function} = DelayMX{size(MX(0.),1),dT,mxT}(τ, MX)
DelayMX(τ::Delay{dT}, MX::mxT) where {dT,mxT <: SMatrix} = DelayMX{size(MX, 1),dT,mxT}(τ, MX)

DelayMX(τ::Delay{dT}, MX::mxT) where {dT,mxT<:AbstractMatrix{<:Real}} = DelayMX(τ,SMatrix{size(MX)...}(MX))
(DMX::DelayMX{d,<:dT,<:AbstractMatrix{<:Real}})(t) where {d,dT} = DMX.MX
(DMX::DelayMX{d,<:dT,<:Function})(t) where {d,dT} = DMX.MX(t)

struct Additive{d,T} <: AdditiveVector{d}
    V::T
end
Additive(v::vT) where vT <: Function = Additive{length(v(0.)),vT}(v)
Additive(v::vT) where vT <: SVector{d} where d = Additive{d,vT}(v)
Additive(v::vT) where vT <: AbstractVector{<:Real} = Additive(SVector(v...))
Additive(d::Integer) = Additive{d,Vector{Nothing}}(Vector{Nothing}(undef, d))
(AV::Additive{d,<:AbstractVector{<:Real}})(t) where d = AV.V
(AV::Additive{d,<:Function})(t) where d = AV.V(t)
(AV::Additive{d,<:Array{<:Nothing}})(t) where d = zeros(d) # TODO check the effect on performance
Base.convert(::Type{Additive}, v::vT) where vT <: AbstractVector{<:Real} = Additive(v)

# Linear Delay Differential Equation Problem
abstract type AbstractLDDEProblem{d} end
struct LDDEProblem{d,AT <: ProportionalMX{d,<:MatrixOrFunction},BT <: DelayMX{d,<:Any,<:Any}, cT <: Additive{d,<:VectorOrFunction}} <: AbstractLDDEProblem{d}
    A::AT # Coefficient of the proportional term (present /non-discretised term)
    Bs::Vector{BT} # Coefficient of the delayed terms
    c::cT # Additive vector
    # d::Int64 # dimension of the state space
end
_problemsize(::LDDEProblem{d,AT,BT,cT}) where {d,AT,BT,cT} = d # size of the state space

LDDEProblem(A::AT, Bs::Vector{<:BT}, c::cT = Additive(size(A, 2))) where {AT <: ProportionalMX,BT <: DelayMX,cT <: Additive} = LDDEProblem{size(A, 2),AT,BT,cT}(A, Bs, c)
LDDEProblem(A::AT, B::BT, c::cT = Additive(size(A, 2))) where {AT <: ProportionalMX,BT <: DelayMX,cT <: Additive} = LDDEProblem{size(A, 2),AT,BT,cT}(A, [B], c)


text:
using Test
using EarthModels

@testset "Basic properties" begin
    let m = SteppedLayeredModel(1, 2, [0.5, 1], [1, 2], [3, 4], [5, 6], false,
            [], [], [], [], [], true, [10, 20], [100, 200])
        @test !isanisotropic(m)
        @test surface_radius(m) == 1.0
        @test radius(m, 0.25) == 0.75
        @test depth(m, 0.25) == 0.75
        @test hasattenuation(m)
    end
    let a = 1, n = 2, arr = rand(2), m = SteppedLayeredModel(a, n, [0.005 + 0.99rand(), 1],
            arr, arr, arr, true, arr, arr, arr, arr, arr, false, [], [])
        @test isanisotropic(m)
        @test !hasattenuation(m)
    end
end


text:
# # Fluspect

## load packages
using Land.CanopyLayers
using PlotPlants
FT = Float32;
#------------------------------------------------------------------------------

# ## Excitation wavelength
wls    = create_wave_length(FT);
can    = create_canopy_rt(FT);
rt_dim = create_rt_dims(can, wls);
leaf   = create_leaf_bios(FT, rt_dim);
fluspect!(leaf, wls);

_fig,_axes = create_canvas("Fluspect example"; ncol=3);
_ax1,_ax2,_ax3 = _axes;
_ax1.plot(wls.WL, leaf.ρ_SW, "k-", label="Reflectance");
_ax1.plot(wls.WL, leaf.τ_SW, "k:", label="Transmittance");
_ax2.contourf(wls.WLE, wls.WLF, leaf.Mb);
_ax3.contourf(wls.WLE, wls.WLF, leaf.Mf);
_ax1.legend(loc="upper right");
_ax1.set_ylim(0,0.65);
set_xlabels!(_axes, ["Wavelength (nm)", "Excitation wavelength (nm)",
                     "Excitation wavelength (nm)"], fontsize=12);
set_ylabels!(_axes, ["Reflectance or Transmittance", "SIF wavelength (nm)",
                     "SIF wavelength (nm)"], fontsize=12);
set_titles!(_axes; labels=["Spectrum", "Backward matrix", "Forward matrix"],
                   usetex=false);
_fig.set_tight_layout(true);
_fig
#------------------------------------------------------------------------------

# ## Change leaf properties
## here we change all the properties at the same time as an example
leaf.N   = 2.0;
leaf.Cab = 50.0;
leaf.Car = 15.0;
leaf.Ant = 0.1;
leaf.Cs  = 0.1;
leaf.Cw  = 0.02;
leaf.Cm  = 0.02;
leaf.Cx  = 0.1;
leaf.fqe = 0.8;
fluspect!(leaf, wls);

_fig,_axes = create_canvas("Change leaf properties"; ncol=3);
_ax1,_ax2,_ax3 = _axes;
_ax1.plot(wls.WL, leaf.ρ_SW, "k-", label="Reflectance");
_ax1.plot(wls.WL, leaf.τ_SW, "k:", label="Transmittance");
_ax2.contourf(wls.WLE, wls.WLF, leaf.Mb);
_ax3.contourf(wls.WLE, wls.WLF, leaf.Mf);
_ax1.legend(loc="upper right");
_ax1.set_ylim(0,0.65);
set_xlabels!(_axes, ["Wavelength (nm)", "Excitation wavelength (nm)",
                     "Excitation wavelength (nm)"], fontsize=12);
set_ylabels!(_axes, ["Reflectance or Transmittance", "SIF wavelength (nm)",
                     "SIF wavelength (nm)"], fontsize=12);
set_titles!(_axes; labels=["Spectrum", "Backward matrix", "Forward matrix"],
                   usetex=false);
_fig.set_tight_layout(true);
_fig
#------------------------------------------------------------------------------

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land