text:
"""
    StaticStepRangeLen{T,B,S,E,L,F}
"""
abstract type StaticStepRangeLen{T,R,S,E,L,F} <: StaticOrdinalRange{T,R,S,E,L} end

@pure offset(::StaticStepRangeLen{T,R,S,E,L,F}) where {T,R,S,E,L,F<:SInteger} = F()
@pure offset(::Type{<:StaticStepRangeLen{T,R,S,E,L,F}}) where {T,R,S,E,L,F<:SInteger} = F()

@pure reference(::StaticStepRangeLen{T,R,S,E,L,F}) where {T,R<:SVal,S,E,L,F} = R()
@pure reference(::Type{<:StaticStepRangeLen{T,R,S,E,L,F}}) where {T,R<:SVal,S,E,L,F} = R()

@pure getrefhi(r::StaticStepRangeLen{T,TPVal{Hb,Lb}}) where {T,Hb,Lb} = Hb()::Hb
@pure getreflo(r::StaticStepRangeLen{T,TPVal{Hb,Lb}}) where {T,Hb,Lb} = Lb()::Lb

@pure getstephi(r::StaticStepRangeLen{T,R,TPVal{Hs,Ls}}) where {T,R,Hs,Ls} = Hs()::Hs
@pure getsteplo(r::StaticStepRangeLen{T,R,TPVal{Hs,Ls}}) where {T,R,Hs,Ls} = Ls()::Ls

@inline first(r::StaticStepRangeLen) = unsafe_getindex(r, SOne)

"StepSRangeLen - Static parametric type variant of StaticStepRangeLen"
struct StepSRangeLen{T,R,S,E,L,F} <: StaticStepRangeLen{T,R,S,E,L,F}
    function StepSRangeLen{T,R,S,E,L,F}() where {T,R,S,E,L,F}
        values(L) >= 0 || throw(ArgumentError("StepSRangeLen: length cannot be negative, got $(values(L))"))
        1 <= values(F) <= max(1, values(L)) || throw(ArgumentError("StaticStepRangeLen: offset must be in [1,$(values(L))], got $offset"))
        new{eltype(T),R,S,E,L,F}()  # ensure that T <: BaseType
    end
end

StepSRangeLen(ref::R, step::S, len::SInteger, offset::SInteger = SOne) where {R,S} =
    StepSRangeLen{eltype(ref+SZero*step)}(ref, step, len, offset)

StepSRangeLen{T}(ref::R, step::S, len::SInteger, offset::SInteger = SOne) where {T,R<:SReal,S<:SReal} =
    StepSRangeLen{T,R,S}(ref, step, ref + (len-offset) * step, len, offset)

function StepSRangeLen{T}(start::TPVal, step::TPVal, len::SInteger, offset::SInteger) where T
    u = len - offset
    shift_hi, shift_lo = u*gethi(step), u * getlo(step)
    x_hi, x_lo = add12(gethi(start), shift_hi)
    StepSRangeLen{T,typeof(start),typeof(step)}(start, step, ofeltype(T, x_hi + (x_lo + (shift_lo + getlo(start)))), len, offset)
end

StepSRangeLen{T,B,S}(start::B, step::S, stop::E, len::SInteger, offset::SInteger) where {T,B,S,E} =
    StepSRangeLen{T,B,S,E,typeof(len),typeof(offset)}()

"StepMRangeLen - Mutable variant of StaticStepRangeLen"
mutable struct StepMRangeLen{T,R,S,L,F} <: StaticStepRangeLen{T,R,S,Dynamic,L,F}
    ref::R       # reference value (might be smallest-magnitude value in the range)
    step::S      # step value
    len::L     # length of the range
    offset::F  # the index of ref

    function StepMRangeLen{T,R,S}(ref::R, step::S, len::Integer, offset::Integer = SOne) where {T,R,S}
        len >= 0 || throw(ArgumentError("length cannot be negative, got $len"))
        1 <= offset <= max(1,len) || throw(ArgumentError("StepMRangeLen: offset must be in [1,$len], got $offset"))
        new{T,R,S,typeof(len),typeof(offset)}(ref, step, len, offset)
    end
end

@inline offset(r::StepMRangeLen{T,R,S,L,F}) where {T,R,S,L,F<:BaseInteger} = getfield(r, :offset)
@inline offset!(r::StepMRangeLen{T,R,S,L,F}, f::F) where {T,R,S,L,F<:BaseInteger} = setfield!(r, :offset, f)

@inline reference(r::StepMRangeLen{T,R,S,L,F}) where {T,R<:BaseAny,S,L,F} = getfield(r, :ref)
@inline reference!(r::StepMRangeLen{T,R,S,L,F}, ref::R) where {T,R<:BaseAny,S,L,F} = setfield!(r, :ref, ref)

StepMRangeLen(ref::R, step::S, len::Integer, offset::Integer = SOne) where {R,S} =
    StepMRangeLen{eltype(ref+0*step)}(ref, step, len, offset)

StepMRangeLen{T}(ref::R, step::S, len::Integer, offset::Integer=SOne) where {T,R,S} =
    StepMRangeLen{T,R,S}(ref, step, len, offset)

@inline getrefhi(r::StepMRangeLen{T,TwicePrecision{Tb}}) where {T,Tb} = r.ref.hi::Tb
@inline getreflo(r::StepMRangeLen{T,TwicePrecision{Tb}}) where {T,Tb} = r.ref.lo::Tb

@inline getstephi(r::StepMRangeLen{T,B,TwicePrecision{Ts}}) where {T,B,Ts} = r.step.hi::Ts
@inline getsteplo(r::StepMRangeLen{T,B,TwicePrecision{Ts}}) where {T,B,Ts} = r.step.lo::Ts

@inline length!(r::StepMRangeLen{T,B,S,L}, len::L) where {T,B,S,L<:BaseInteger} = setfield!(r, :len, len)
@inline length(r::StepMRangeLen{T,B,S,L,F}) where {T,B,S,L,F} = getfield(r, :len)

step(r::StepMRangeLen{T,B,S}) where {T,B,S<:BaseAny} = getfield(r, :step)

StaticStepRangeLen(ref::SVal, step::SVal, len::SInteger, offset::SInteger = SOne) = StepSRangeLen(ref, step, len, offset)
StaticStepRangeLen{T}(ref::SVal, step::SVal, len::SInteger, offset::SInteger = SOne) where T = StepSRangeLen{T}(ref, step, len, offset)

StaticStepRangeLen(ref::R, step::S, len::Integer, offset::Integer = SOne) where {R,S} = StepMRangeLen(ref, step, len, offset)
StaticStepRangeLen{T}(ref::R, step::S, len::Integer, offset::Integer = SOne) where {T,R,S} = StepMRangeLen{T}(ref, step, len, offset)

showrange(io::IO, r::StaticStepRangeLen) = print(io, "$(first(r)):$(step(r)):$(last(r))")

isstatic(::StepSRangeLen) = true
isstatic(::Type{<:StepSRangeLen}) = true

isstatic(::StepMRangeLen) = false
isstatic(::Type{<:StepMRangeLen}) = false

#promote_rule(a::Type{StepRangeLen{T,R,S}}, ::Type{OR}) where {T,R,S,OR<:AbstractRange} =
#    promote_rule(a, StepRangeLen{eltype(OR), eltype(OR), eltype(OR)})
Base.StepRangeLen(r::StaticRange) = StepRangeLen{eltype(r)}(r)
Base.StepRangeLen{T}(r::StaticRange) where {T} =
    StepRangeLen(T(first(r)), T(step(r)), values(length(r)))
Base.StepRangeLen{T,R,S}(r::StaticRange) where {T,R,S} =
    StepRangeLen{T,R,S}(R(first(r)), S(step(r)), values(length(r)))

StepSRangeLen(r::AbstractRange) = StepSRangeLen{eltype(r)}(r)
StepSRangeLen{T}(r::AbstractRange) where {T} =
    StepSRangeLen(SVal(ofeltype(T, first(r))), SVal(ofeltype(T, step(r))), SVal(length(r)))
StepSRangeLen{T,R,S}(r::AbstractRange) where {T,R,S} =
    StepSRangeLen(SVal(R(first(r))), SVal(S(step(r))), SVal(length(r)))

StepMRangeLen(r::AbstractRange) = StepMRangeLen{eltype(r)}(r)
StepMRangeLen{T}(r::AbstractRange) where {T} =
    StepMRangeLen(ofeltype(T, first(r)), ofeltype(T, step(r)), length(r))
StepMRangeLen{T,R,S}(r::AbstractRange) where {T,R,S} =
    StepMRangeLen{T}(R(first(r)), S(step(r)), length(r))

@inline function _getindex_hiprec(r::StaticStepRangeLen{T,R}, i::Integer) where {T,R<:Union{<:TPVal,TwicePrecision}}
    u = i - offset(r)
    shift_hi, shift_lo = u*getstephi(r), u*getsteplo(r)
    x_hi, x_lo = add12(getrefhi(r), shift_hi)
    x_hi, x_lo = add12(x_hi, x_lo + (shift_lo + getreflo(r)))
    tpval(x_hi, x_lo)
end

@inline function _getindex_hiprec(r::StaticStepRangeLen{T,R}, i::Integer) where {T,R} # without rounding by T
    reference(r) + (i - offset(r)) * step(r)
end

@inline function Base.iterate(r::StepSRangeLen{T,R,S,E,L,F}) where {T,R,S,E,L,F}
    if values(L) < 1
        return nothing
    else
        return (unsafe_getindex(r, 1)::T, 1)::Tuple{T,Int}
    end
end

@inline function Base.iterate(r::StepSRangeLen{T,R,S,E,L,F}, i::Int) where {T,R,S,E,L,F}
    if i < values(L)
        return (unsafe_getindex(r, i+1), i + 1)::Tuple{T,Int}
    else
        return nothing
    end
end


text:
#**************************************************************************************
# test_Physics_Model.jl
# ====================== part of the GeoEfficiency.jl package.
#
#
#**************************************************************************************

@debug("Point")
@testset "Point" begin
	@test G.id(Point(5,3)) == "Point[Height=5.0, Rho=3.0]"
	@test show(Point(5,3)) === nothing

	local pnt11 = Point(1, 1)
	local pnt5 = Point(5.0, 0.0)
	local pnt51 = Point(5.0, 1.0)

	@test Point(5).Height === 5.0
	@test Point(5).Rho    === 0.0
	@test Point(-5).Height === -5.0
	@test Point(-5).Rho    === 0.0
	@test Point(-5,-1).Height === -5.0
	@test Point(-5,-1).Rho    === -1.0
	@test_throws MethodError     Point(1+0im,4)
	@test_throws MethodError     Point(4, 1+0im)
	@test_throws MethodError     Point(1+0im, 1+0im)
	@test_throws ErrorException  Point(5.0, 0.0).Height = 1.0
	@test_throws ErrorException  Point(5.0, 0.0).Rho    = 1.0

	@test Point(5)      === pnt5
	@test Point(5//1)   === pnt5
	@test Point(5.0)    === pnt5

	@test Point(5, 0)   === pnt5
	@test Point(5, 0//1)=== pnt5
	@test Point(5, 0.0) === pnt5

	@test Point(5//1, 0)    === pnt5
	@test Point(5//1, 0//1) === pnt5
	@test Point(5//1, 0.0)  === pnt5

	@test Point(5.0, 0)    === pnt5
	@test Point(5.0, 0//1) === pnt5
	@test Point(5.0, 0.0)  === pnt5

	@test Point(5, 1)   === pnt51
	@test Point(5, 1//1)=== pnt51
	@test Point(5, 1.0) === pnt51

	@test Point(5//1, 1)    === pnt51
	@test Point(5//1, 1//1) === pnt51
	@test Point(5//1, 1.0)  === pnt51

	@test Point(5.0, 1)    === pnt51
	@test Point(5.0, 1//1) === pnt51
	@test Point(5.0, 1.0)  === pnt51

	local pnt6 = Point(pnt5, 5)
	@test pnt6.Height == pnt5.Height
	@test pnt6.Rho == 5.0

	pnt6 = Point(5, pnt5)
	@test pnt6.Height == 5.0
	@test pnt6.Rho == pnt5.Rho

	@test H.exec_console_unattended(Point, [5, 1]) == pnt51
	@test H.exec_console_unattended(Point, [5, 0]) == pnt5
end #testset_Point

@debug("CylDetector")
@testset "CylDetector" begin

	@test G.id(CylDetector(5, 3)) == "CylDetector[CryRadius=5.0, CryLength=3.0]"
	@test show(CylDetector(5, 3)) === nothing

	let cyl0 = CylDetector(5), cyl1 = Detector(5)

		@test_throws ErrorException  cyl0.CryRadius = 1
		@test_throws ErrorException  cyl0.CryLength = 1
		@test_throws MethodError     CylDetector(1+0im)
		@test_throws MethodError     CylDetector(1+0im,4)
		@test_throws MethodError     CylDetector(4, 1+0im)
		@test_throws MethodError     CylDetector(1+0im, 1+0im)
		@test_throws G.InValidDetectorDim  CylDetector(0)
		@test_throws G.InValidDetectorDim  CylDetector(-5)
		@test_throws G.InValidDetectorDim  CylDetector(5, -1)
		@test isa(cyl1, Detector)
		@test isa(cyl1, CylDetector)
		@test cyl1.CryRadius === 5.0
		@test cyl1.CryLength === 0.0
	end #let

	let cyl0 = CylDetector(5),
		cyl1 = Detector(5),
		cyl2 = Detector(5, 0),
		cyl3 = Detector(5, 0.0),
		cyl4 = Detector(5.0, 0),
		cyl5 = Detector(5.0, 0.0),
		cyl6 = Detector(5, 0, 0),
		cyl7 = Detector(5, 0, 0, 0),
		cyl8 = Detector(5//1, 0, 0, 0),
		cyl9 = Detector(5//1)

		@test cyl0 === cyl1
		@test cyl1 === cyl2
		@test cyl2 === cyl3
		@test cyl3 === cyl4
		@test cyl4 === cyl5
		@test cyl5 === cyl6
		@test cyl6 === cyl7
		@test cyl7 === cyl8
		@test cyl8 === cyl9
		@test G.volume(CylDetector(5.0,1))  <=   G.volume(CylDetector(15.0,1))
		@test G.volume(CylDetector(10.0,1)) <=   G.volume(CylDetector(15.0,1))
		@test G.volume(CylDetector(15.0))   <=   G.volume(CylDetector(10.0))
		@test CylDetector(5.0,1) < CylDetector(15.0,1)
	end #let

	let detectors = [Detector(6,2),Detector(5,1), Detector(7,10)]

		@test eltype(Vector{Detector}(detectors)) === Detector
		@test Vector{Detector}(detectors) == detectors
	end #let

	@debug("CylDetector()")
	@test H.exec_console_unattended(CylDetector, [5, 1]) == CylDetector(5, 1)
	@test H.exec_console_unattended(CylDetector, [5, 0])  == CylDetector(5)
end #testset_CylDetector

@debug("BoreDetector")
@testset "BoreDetector" begin

		@test G.id(BoreDetector(5,3,2)) == "BoreDetector[CryRadius=5.0, CryLength=3.0, HoleRadius=2.0]"
  		@test show(BoreDetector(5,3,2)) === nothing

	let bore0 = BoreDetector(5,4,3),
		bore1 = Detector(5,4,3)

		@test_throws ErrorException  bore0.CryRadius = 1
		@test_throws ErrorException  bore0.CryLength = 1
		@test_throws ErrorException  bore0.HoleRadius= 1
		@test_throws MethodError     BoreDetector(1+1im,4,3)
		@test_throws MethodError     BoreDetector(5+1im,4,3)
		@test isa(bore1, Detector)
		@test isa(bore1, BoreDetector)
		@test 5.0 == bore1.CryRadius
		@test 4.0 == bore1.CryLength
		@test 3.0 == bore1.HoleRadius
	end #let

	let bore0 = BoreDetector(5,4,3),
		bore1 = Detector(5,4,3),
		bore2 = Detector(5.0,4,3),
		bore3 = Detector(5,4,3,0),
		bore4 = Detector(5,4.0,3),
		bore5 = Detector(5.0,4,3.0),
		bore6 = Detector(5.0,4.0,3.0),
		bore7 = Detector(5.0,4.0,3.0, 0),
		bore8 = Detector(5//1,4.0,3.0, 0)

		@test bore0 === bore1
		@test bore1 === bore2
		@test bore2 === bore3
		@test bore3 === bore4
		@test bore4 === bore5
		@test bore5 === bore6
		@test bore6 === bore7
		@test bore7 === bore8
		@test G.volume(Detector(5.0,1,.1))  <=   G.volume(Detector(15.0,1,0.1))
		@test G.volume(Detector(10.0,1,.0)) <=   G.volume(Detector(15.0,1,0.1))
		@test Detector(5.0,1,0.1) < Detector(15.0,1,0.1)
	end #let

	let	detectors = [Detector(6,2,1), Detector(5,1,.2), Detector(7,10,5)]

		@test eltype(Vector{Detector}(detectors)) === Detector
		@test Vector{Detector}(detectors) == detectors
	end #let

	@debug("BoreDetector()")
	let bore0 = BoreDetector(5,4,3)

		@test H.exec_console_unattended(BoreDetector, [5, 4, 3]) === bore0
		@test H.exec_console_unattended(BoreDetector, [5, 4, 6, 3]) === bore0
	end #let

end #testset_BoreDetector

@debug("WellDetector")
@testset "WellDetector" begin
	@test G.id(WellDetector(5, 3, 2, 1)) == "WellDetector[CryRadius=5.0, CryLength=3.0, HoleRadius=2.0, HoleDepth=1.0]"
	@test show(WellDetector(5, 3, 2, 1)) === nothing

	let	Well0 = WellDetector(5,4,3,2),
		Well1 = Detector(5,4,3,2),
		Well2 = Detector(5.0,4,3,2),
		Well3 = Detector(5,4.0,3,2),
		Well4 = Detector(5,4,3.0,2),
		Well5 = Detector(5,4,3,2.0),
		Well6 = Detector(5.0,4,3.0,2),
		Well7 = Detector(5,4.0,3,2.0),
		Well8 = Detector(5//1,4,3,2)

		@test_throws ErrorException  Well0.CryRadius = 1
		@test_throws ErrorException  Well0.CryLength = 1
		@test_throws ErrorException  Well0.HoleRadius= 1
		@test_throws ErrorException  Well0.HoleDepth = 1
		@test_throws MethodError     WellDetector(1+1im,4,3,2)
		@test_throws MethodError     WellDetector(5+1im,4,3,2)
		@test isa(Well1, Detector)
		@test isa(Well1, WellDetector)
		@test 5.0 == Well1.CryRadius
		@test 4.0 == Well1.CryLength
		@test 3.0 == Well1.HoleRadius
		@test 2.0 == Well1.HoleDepth

		@test Well0 === Well1
		@test Well1 === Well2
		@test Well2 === Well3
		@test Well3 === Well4
		@test Well4 === Well5
		@test Well5 === Well6
		@test Well6 === Well7
		@test Well7 === Well8
	end #let

	@test G.volume(Detector(5.0, 2, 0.1, 1))  <=   G.volume(Detector(15.0, 2, 0.1, 1))
	@test G.volume(Detector(10.0, 2, 0.1, 1)) <=   G.volume(Detector(15.0, 2, 0.1, 1))
	@test Detector(5.0, 2, 0.1, 1) < Detector(15.0, 2, 0.1, 1)

	let	detectors = [Detector(6, 2, 1, 0.1), Detector(5, 1, 0.2, 0.1), Detector(7, 10, 5, 0.1)]

		@test eltype(Vector{Detector}(detectors)) === Detector
		@test Vector{Detector}(detectors) == detectors
	end #let

	@debug("WellDetector()")
	let Well0 = WellDetector(5, 4, 3, 2)

		@test H.exec_console_unattended(WellDetector, [5, 4, 3, 2])    == Well0
		@test H.exec_console_unattended(WellDetector, [5, 4, 3, 6, 2]) == Well0
	end #let
end #testset_WellDetector

@debug("Detector")
@testset "Detector" begin

	@test_throws MethodError  Detector(1+1im)
	@test_throws MethodError  Detector(5+1im, 0)
	@test_throws MethodError  Detector(1+1im, 1)
	@test_throws MethodError  Detector(5+1im, 4, 3)
	@test_throws MethodError  Detector(1+1im, 4, 3, 2)
	@test_throws MethodError  Detector(5+1im, 4, 3, 2)
	@test_throws MethodError  Detector(1+1im, 4, 3, 2)
	@test_throws MethodError  Detector(4, 1+1im, 3, 2)
	@test_throws MethodError  Detector(4, 3, 2+1im, 1)
	@test_throws MethodError  Detector(4, 3, 2, 1+1im)

	@test Detector(5)       === CylDetector(5)
	@test Detector(5,0)     === CylDetector(5)
	@test Detector(5,0,0,0) === CylDetector(5)
	@test Detector(5,1)     === CylDetector(5,1)
	@test Detector(5,1,0)   === CylDetector(5,1)
	@test Detector(5,1,0,0) === CylDetector(5,1)
	@test Detector(5,2,1)   === BoreDetector(5,2, 1)
	@test Detector(5,2,1,0) === BoreDetector(5,2, 1)
	@test Detector(5,4,3,2) === WellDetector(5,4,3,2)

	let cyl0 = CylDetector(5),
		bore0 = BoreDetector(5, 4, 3),
		Well0 = WellDetector(5, 4, 3, 2)

    	@test Detector(cyl0)   === cyl0
		@test Detector(bore0)  === bore0
		@test Detector(Well0)  === Well0
	end #let
	for i=1:1000
		let det1 = Detector(rand()),
			det2 = Detector(rand(2)...),
			det3 = Detector(rand(15:20), rand(10:14), rand(1:13)),
			det4 = Detector(rand(15:20), rand(10:14), rand(9:13), rand(1:9))

			@test Detector(det1)   === det1
			@test Detector(det2)   === det2
			@test Detector(det3)   === det3
			@test Detector(det4)   === det4
		end #let
	end #for

	@test G.id(Detector(5,3,2,1)) == "WellDetector[CryRadius=5.0, CryLength=3.0, HoleRadius=2.0, HoleDepth=1.0]"
	@test G.id(Detector(5,3,2)) == "BoreDetector[CryRadius=5.0, CryLength=3.0, HoleRadius=2.0]"
	@test G.id(Detector(5,3)) == "CylDetector[CryRadius=5.0, CryLength=3.0]"
	@test G.id(Detector(5)) == "CylDetector[CryRadius=5.0, CryLength=0.0]"

	@debug("Detector(detectors::Vector{<:Detector})")
	@testset "Detector(detectors::Vector{<:Detector})" begin
	mix_detectors = [Detector(6,2), Detector(5,1,.2), Detector(7,10,5,.1)]
	cyl_detectors = [CylDetector(2,6), Detector(6,2), CylDetector(2), Detector(6)]

		@test eltype(Vector{Detector}(mix_detectors)) 	=== Detector
		@test eltype(Vector{Detector}(cyl_detectors)) 	=== Detector

		@test Vector{Detector}(mix_detectors) == mix_detectors
		@test Vector{Detector}(cyl_detectors) == cyl_detectors

		@test eltype(Detector(mix_detectors)) 	=== Detector
		@test eltype(Detector(cyl_detectors)) 	=== Detector

		@test Detector(mix_detectors) == mix_detectors
		@test Detector(cyl_detectors) == cyl_detectors
	end #Detector(detectors::Vector{<:Detector})

	@debug("Detector()")
	@test H.exec_console_unattended(Detector, [5, 4, 3, 6, 2]) 		== Detector(5, 4, 3, 2)
	@test H.exec_console_unattended(Detector, [5, 4, 0])				== Detector(5, 4)
	@test H.exec_console_unattended(Detector, "5\n" * "4\n" * "0\n") 	== Detector(5, 4)
	@test H.exec_console_unattended(Detector, [5, 4, 33, 3, 22, 2]) 		== Detector(5, 4, 3, 2)
end #testset_Detector

@debug("source")
@testset "source" begin
	setSrcToPoint(true)
	let	pnt1 = Point(5)

		@test source(pnt1) == (pnt1, 0.0, 0.0)
	end #let

	setSrcToPoint(false)
	write(stdin.buffer,
		"1\n" * "0\n" * #=axial point=#
		"2\n" * "3\n" #=SrcRadius=2 SrcHeight=3=#)
	@test  source() == (Point(1.0, 0.0), 2.0, 3.0)

	setSrcToPoint(false)
	write(stdin.buffer,
		"1\n" * "5\n" * #=non-axial point=#
		"2\n" * "3\n" #=SrcRadius=2 SrcHeight=3=#)
	@test  source() == (Point(1.0, 0.0), 2.0, 3.0) # axial anchor point

	setSrcToPoint(false)
	write(stdin.buffer,
		"1\n" * "5\n" * #=non-axial point=#
		"0\n" 			#=SrcRadius=2=# )
	@test  source() == (Point(1.0, 5.0), 0.0, 0.0) # non-axial  point source

	setSrcToPoint(true)
	write(stdin.buffer,	"1\n" * "5\n") # non-axial point
	@test  source() == (Point(1.0, 5.0), 0.0, 0.0)
end #testset_source

@debug("Invalid Detector Dimensions")
@testset "Invalid Detector Dimensions $dim"  for dim =
Number[0, -1, 0//1, -1//1, -e, 0.0, -1.0, -Inf, Inf,]

	@test_throws G.InValidDetectorDim	CylDetector(dim)
	@test_throws G.InValidDetectorDim	CylDetector(dim, 0)
	@test_throws G.InValidDetectorDim	CylDetector(dim, 1)
	@test_throws G.InValidDetectorDim	BoreDetector(dim, 1, 0.2)
	@test_throws G.InValidDetectorDim	BoreDetector(dim, 1, 0)
	@test_throws G.InValidDetectorDim	WellDetector(dim, 2, 1, 0)

	@test_throws G.InValidDetectorDim	Detector(dim)
	@test_throws G.InValidDetectorDim	Detector(dim, 0)
	@test_throws G.InValidDetectorDim	Detector(dim, 1)
	@test_throws G.InValidDetectorDim	Detector(dim, 0, 0)

	dim == 0.0 && break
	@test_throws G.InValidDetectorDim	CylDetector(5, dim)
	@test_throws G.InValidDetectorDim  	BoreDetector(5, dim, 0)
	@test_throws G.InValidDetectorDim	CylDetector(5, dim)
	@test_throws G.InValidDetectorDim	BoreDetector(5, dim, 0)
	@test_throws G.InValidDetectorDim	BoreDetector(5,dim, 1)
	@test_throws G.InValidDetectorDim	WellDetector(5,dim, 1, 0.1)
end #testset_Invalid_Detector_Dimensions

ext: jl

max_stars_repo_name: DrKrar/GeoEfficiency

max_issues_repo_name: DrKrar/GeoEfficiency

max_forks_repo_name: DrKrar/GeoEfficiency


text:
type SDF_VD <: SDFDynamicSystem
    omg::Float64
    zta::Float64
    xi::Float64
    alpha::Float64
    u0::Float64
    v0::Float64

    SDF_VD(omg::Float64, zta::Float64, xi::Float64, alpha::Float64) = new(omg, zta, xi, alpha, 0.0, 0.0)
    SDF_VD(omg::Float64, zta::Float64, xi::Float64) = new(omg, zta, xi, 1.0, 0.0, 0.0)
end

function SDF_VD(m::Float64, c::Float64, k::Float64, c_d::Float64, alpha::Float64)
    omg = sqrt(k/m)
    zta = c/(2.0*m*omg)
    xi = c_d/(2.0*m*omg)
    SDF_VD(omg, zta, xi, alpha)
end

function systemEquation(s::SDF_VD, p::Function)
    function dydt(y::AbstractArray{Float64,1}, t::Float64)
        if s.alpha == 1.0
            dvdt = p(t)-2.0*(s.zta + s.xi)*s.omg*y[2]-s.omg*s.omg*y[1]
        else
            dvdt = p(t)-2.0*s.zta*s.omg*y[2]-2.0*s.xi*s.omg*sign(y[2])*abs(y[2])^s.alpha-s.omg*s.omg*y[1]
        end
        [y[2]; dvdt]
    end
    dydt
end

function response(s::SDF_VD, p::Function, t::AbstractArray{Float64,1})
    y0 = [s.u0; s.v0]

    f = systemEquation(s, p)

    y = ode4(f, y0, t)

    u = y[:,1]
    v = y[:,2]
    fdz = 2.0*s.zta*s.omg*v
    fdx = 2.0*s.xi*s.omg*sign(v).*abs(v).^s.alpha
    a = map(p,t)-fdz-fdx-s.omg*s.omg*u

    [u v a fdz fdx]
end

max_stars_repo_name: Panchatantra/DSR

max_issues_repo_name: Panchatantra/DSR

max_forks_repo_name: Panchatantra/DSR


text:
using Test
using JustSayIt
using JustSayIt.API
using PyCall
import JustSayIt: DEFAULT_MODEL_NAME, TYPE_MODEL_NAME, MODELDIR_PREFIX, DEFAULT_NOISES, DIGITS_ENGLISH
import JustSayIt: init_jsi, finalize_jsi, voicearg_f_names, voiceargs, recognizer

# Test setup
@voiceargs space=>(valid_input=["world", "universe"]) hello(space::String) = println("hello $space")

@enum Name julia python
@voiceargs (n1=>(valid_input_auto=true), n2=>(valid_input_auto=true)) hi(n1::Name, n2::Name) = println("hi $n1 and $n2")

@voiceargs (
    nr=>(valid_input=[keys(DIGITS_ENGLISH)...], interpret_function=Keyboard.interpret_digits, use_max_speed=true),
    name=>(valid_input_auto=true),
    question=>(model=TYPE_MODEL_NAME, ignore_unknown=true, vararg_end="end", vararg_max=10, vararg_timeout=5.0)
) function ask(nr::Integer, name::Name, question::String...)
    println("[Q$nr] Hi $name, could you please $(join(question," "))?")
end

commands = Dict("help"  => Help.help,
                "hello" => hello)
modeldirs = Dict(DEFAULT_MODEL_NAME => joinpath(MODELDIR_PREFIX, "vosk-model-small-en-us-0.15"),
                 TYPE_MODEL_NAME    => joinpath(MODELDIR_PREFIX, "vosk-model-small-en-us-0.15"))
init_jsi(commands, modeldirs, DEFAULT_NOISES)

@testset "$(basename(@__FILE__))" begin
    @testset "1. voicearg dictionary" begin
        @testset "functions" begin
            @test :hello in voicearg_f_names()
            @test :hi in voicearg_f_names()
            @test :ask in voicearg_f_names()
        end;
        @testset "voiceargs" begin
            @test issetequal(keys(voiceargs(:hello)), [:space])
            @test issetequal(keys(voiceargs(:hi)), [:n1, :n2])
            @test issetequal(keys(voiceargs(:ask)), [:nr, :name, :question])
        end;
        @testset "kwargs" begin
            @test issetequal(keys(voiceargs(:hello)[:space]), [:recognizer, :valid_input])
            @test issetequal(keys(voiceargs(:hi)[:n1]), [:recognizer, :valid_input, :valid_input_auto])
            @test issetequal(keys(voiceargs(:hi)[:n2]), [:recognizer, :valid_input, :valid_input_auto])
            @test issetequal(keys(voiceargs(:ask)[:nr]), [:recognizer, :valid_input, :interpret_function, :use_max_speed])
            @test issetequal(keys(voiceargs(:ask)[:name]), [:recognizer, :valid_input, :valid_input_auto])
            @test issetequal(keys(voiceargs(:ask)[:question]), [:model, :ignore_unknown, :vararg_end, :vararg_max, :vararg_timeout])
        end;
        @testset "kwarg content" begin
            @testset "recognizers" begin
                @test isa(voiceargs(:hello)[:space][:recognizer], PyObject)
                @test isa(voiceargs(:hi)[:n1][:recognizer], PyObject)
                @test isa(voiceargs(:hi)[:n2][:recognizer], PyObject)
                @test isa(voiceargs(:ask)[:nr][:recognizer], PyObject)
                @test isa(voiceargs(:ask)[:name][:recognizer], PyObject)
            end;
            @testset "valid_input" begin
                @test voiceargs(:hello)[:space][:valid_input] == ["world", "universe"]
                @test voiceargs(:hi)[:n1][:valid_input] == ["julia", "python"]
                @test voiceargs(:hi)[:n2][:valid_input] == ["julia", "python"]
                @test voiceargs(:ask)[:nr][:valid_input] == [keys(DIGITS_ENGLISH)...]
                @test voiceargs(:ask)[:name][:valid_input] == ["julia", "python"]
            end;
            @testset "valid_input_auto" begin
                @test voiceargs(:hi)[:n1][:valid_input_auto] == true
                @test voiceargs(:hi)[:n2][:valid_input_auto] == true
                @test voiceargs(:ask)[:name][:valid_input_auto] == true
            end;
            @testset "interpret_function" begin
                @test voiceargs(:ask)[:nr][:interpret_function] == Keyboard.interpret_digits
            end;
            @testset "use_max_speed" begin
                @test voiceargs(:ask)[:nr][:use_max_speed] == true
            end;
            @testset "ignore_unknown" begin
                @test voiceargs(:ask)[:question][:ignore_unknown] == true
            end;
            @testset "model" begin
                @test voiceargs(:ask)[:question][:model] == TYPE_MODEL_NAME
            end;
            @testset "vararg_end" begin
                @test voiceargs(:ask)[:question][:vararg_end] == "end"
            end;
            @testset "vararg_max" begin
                @test voiceargs(:ask)[:question][:vararg_max] == 10
            end;
            @testset "vararg_timeout" begin
                @test voiceargs(:ask)[:question][:vararg_timeout] == 5.0
            end;
        end;
    end;
    @testset "2. recognizers" begin
        @test isa(recognizer(:hello, :space), PyObject)
        @test isa(recognizer(:hi, :n1), PyObject)
        @test isa(recognizer(:hi, :n2), PyObject)
        @test isa(recognizer(:ask, :nr), PyObject)
        @test isa(recognizer(:ask, :name), PyObject)
    end;
    @testset "3. voicearg wrapper" begin
        @testset "defined" begin
            @test length(methods(hello)) == 2
            @test length(methods(hi)) == 2
            @test length(methods(ask)) == 2
        end;
    end;
end;

# Test tear down
finalize_jsi()


text:
"""
This module provides `findvolume` methods to find the volume at a given
pressure, energy, or bulk modulus with(out) units.
"""
module Find

using InteractiveUtils: subtypes
using Unitful: AbstractQuantity, ustrip

using Roots: find_zero,
             AbstractBracketing,
             AbstractNonBracketing,
             AbstractHalleyLikeMethod,
             AbstractNewtonLikeMethod,
             AbstractAlefeldPotraShi,
             AbstractBisection,
             AbstractSecant,
             Brent,
             Newton,
             ConvergenceFailed

import ..EquationForm
using ..Collections: EquationOfState, apply

export findvolume

"""
    findvolume(form, eos, y, x0, method)
    findvolume(form, eos, y, x0::Union{AbstractVector,Tuple})

Find a volume which leads to the given pressure, energy, or bulk modulus based on an `eos`.

# Arguments
- `form::EquationForm`: an `EquationForm` instance.
- `eos::EquationOfState`: an equation of state. If it has units, `y` and `x0` must also have.
- `y`: a pressure, energy, or bulk modulus.
- `x0`: can be either a range of volumes (`Vector`, `Tuple`, etc.) or just a single volume.
    Units can be provided if necessary.
- `method::Roots.AbstractUnivariateZeroMethod`: a method used to find the root of an equation.
    If it is omitted, the algorithm will traverse all possible methods of
    [Roots.jl](https://github.com/JuliaMath/Roots.jl). And the `x0` parameter must be
    an array or a tuple, of which only the maximum and minimum values will be used in the
    root-finding process.
"""
function findvolume(form::EquationForm, eos::EquationOfState, y, x0, method)
    f = v -> apply(form, eos, v) - y
    return find_zero(f, x0, method)
end # function findvolume
function findvolume(
    form::EquationForm,
    eos::EquationOfState,
    y,
    x0::Union{AbstractVector,Tuple},
)
    for T in [subtypes(AbstractBisection); subtypes(AbstractAlefeldPotraShi)]
        @info("Using method \"$T\"...")
        try
            # `maximum` and `minimum` also works with `AbstractQuantity`s.
            return findvolume(form, eos, y, (minimum(x0), maximum(x0)), T())
        catch e
            @info("Method \"$T\" failed because of $e.")
            continue
        end
    end
    for T in [
        Brent
        subtypes(AbstractHalleyLikeMethod)
        Newton
        subtypes(AbstractSecant)
    ]
        @info("Using method \"$T\"...")
        try
            return findvolume(form, eos, y, (minimum(x0) + maximum(x0)) / 2, T())
        catch e
            @info("Method \"$T\" failed because of $e.")
            continue
        end
    end
end # function findvolume

end


text:
using BenchmarkTools
using ITensors
using LinearAlgebra

function main(; d=20, order=4)
  BLAS.set_num_threads(1)
  ITensors.Strided.disable_threads()

  println("#################################################")
  println("# order = ", order)
  println("# d = ", d)
  println("#################################################")
  println()

  @show Threads.nthreads()
  @show Sys.CPU_THREADS
  @show ITensors.blas_get_num_threads()
  @show ITensors.Strided.get_num_threads()
  println()

  i(n) = Index(QN(0) => d, QN(1) => d; tags="i$n")
  is = IndexSet(i, order ÷ 2)
  A = randomITensor(is'..., dag(is)...)
  B = randomITensor(is'..., dag(is)...)

  ITensors.disable_threaded_blocksparse()

  println("Serial contract:")
  @disable_warn_order begin
    C_contract = @btime $A' * $B samples = 5
  end
  println()

  println("Threaded contract:")
  @disable_warn_order begin
    ITensors.enable_threaded_blocksparse()
    C_threaded_contract = @btime $A' * $B samples = 5
    ITensors.disable_threaded_blocksparse()
  end
  println()
  @show C_contract ≈ C_threaded_contract
  return nothing
end


text:
# Prog   : solana_airdrop.jl
# Auth   : MrChaos
# E-Mail : mrchaos@naver.com
# 백그라운드 실행 : nohup julia solana_airdrop.jl > log &

using Printf:@printf
using Dates

# true : devnet, false : testnet
isdev = true

"Run a Cmd object, returning the stdout & stderr contents plus the exit code"
function execute(cmd::Cmd)
  out = Pipe()
  err = Pipe()

  process = run(pipeline(ignorestatus(cmd), stdout=out, stderr=err))
  close(out.in)
  close(err.in)

  (
    stdout = String(read(out)),
    stderr = String(read(err)),
    code = process.exitcode
  )
end
function execute(str::String)
  cmd = Cmd(Vector{String}(split(str)))
  return execute(cmd)
end
function get_balance(sbal::String)
    regex = r"([0-9\.]+[\s]).*"
    m = match(regex,sbal)
    return parse(Float64,strip(m[1]))
end

if isdev
  url = "https://api.devnet.solana.com"
  airdrop_account = "9B5XszUGdMaxCZ7uSQhPzdks5ZQSmWxrmzCSvtJ6Ns6g"
  amount = 2
  max_amount = 1000000
else
  url = "https://api.testnet.solana.com"
  airdrop_account = "4ETf86tK7b4W72f27kNLJLgRWi9UfJjgH4koHGUXMFtn"
  amount = 1
  max_amount = 2000000
end
amount_sum = 0

# Airdrop을 모으는 지갑주소
my_account = "92GHwaKBP4PnkTZur2XrW8atw5GGth7YVv9g3z9Y8akc"

# Airdrop을 받는 임시 지갑
tmp_key = "tmp-account.json"
gas_fee = 0.0001

wait_airdrop = 10   # airdrop 성공후 대기 시간
wait_check = 60*5  # 잔고 부족 후 대기 시간
wait_next_max =  60*60*24  # max amount 달성 후 대기 시간

cmd_wallet_create = "solana-keygen new --no-passphrase -o $(tmp_key)"

if isfile(tmp_key)
  code_wallet = 0
else
  # 지갑생성
  out,error,code_wallet = execute(cmd_wallet_create)
end

while true
  # 지갑생성 성공
  println(string("-------------- ",now()," --------------"))
  flush(stdout)
  if code_wallet == 0
    # airdrop address에서 잔고 조회
    out,err,code = execute("solana balance $(airdrop_account)  --url $(url)")
    # 잔고조회 성공
    if code == 0
      # 잔고
      b = get_balance(out)
      @printf("Airdrop balance : %.10f\n", b)
      flush(stdout)
      # 잔고가 1보다 큰경우 airdrop요청
      if b > 1
        # 지갑주소 가져오기
        out,err,code=execute("solana-keygen pubkey $(tmp_key)")
        if code == 0
          recv_account = string(strip(out))
          out,err,code = execute("solana airdrop $(amount) $(recv_account)  --url $(url)")
          # 에어드롭 성공
          if code == 0
            println("에어드롭성공(SUM:$(amount_sum)) ($recv_account) : $(out) SOL")
            flush(stdout)
          # 에어드롭 실패
          # 잔고가 있는데 에러가 나면 최대치에 해당 했다고 가정 하고 새로운 주소를 만들고 다시 한다.
          else
            # 현재 지갑의 잔고조회
            out,err,code = execute("solana balance $(recv_account)  --url $(url)")
            if code == 0
              b = get_balance(out)
              b = b - gas_fee
              if b <= 0
                println("전송잔고 부족 : ",b)
                flush(stdout)
                sleep(wait_airdrop)
                continue
              end
              # 잔고전송
              pay = "solana transfer --allow-unfunded-recipient --keypair $(tmp_key) --url $(url)  $(my_account) $(b)"
              out,err,code = execute(pay)
              if code==0
                println("성공 : $(my_account)로 $(b) SOL 전송")
                flush(stdout)
                global amount_sum
                amount_sum = amount_sum + b
                if amount_sum >= max_amount
                  amount_sum = 0
                  println("장기수면")
                  flush(stdout)
                  sleep(wait_next_max)
                end
              else
                println("실패 : $(my_account)로 $(b) SOL 전송")
                flush(stdout)
              end
              # 기존지갑파일 삭제
              if isfile(tmp_key)
                mv(tmp_key,string(tmp_key,".OLD");force=true)
              end
              # 지갑생성
              global code_wallet
              out,err,code_wallet = execute(cmd_wallet_create)
            end
          end
        #지갑주소 가져오기 실패
        else
          println("Airdrop받을 지갑주소 가져오기 실패 : $(out) : $(err)")
          flush(stdout)
        end
      # 잔고가 1보다 작은 경우
      else
        @printf("Airdrop잔고부족 : %.10f\n", b)
        flush(stdout)
        sleep(wait_check)
      end
    else
      println("Airdrop 잔고조회에러")
      flush(stdout)
    end
  else
    global code_wallet
    if isfile(tmp_key)
      code_wallet = 0
    else
      # 지갑생성
      out,error,code_wallet = execute(cmd_wallet_create)
      println("지갑생성에러:",err)
      flush(stdout)
    end
  end
  flush(stdout)
  sleep(wait_airdrop)
end


text:
#
# Example of a set of reference values that StateCheck will use to test current vlaues against.
# These are generated by the function scprintref() and parsed by the scdocheck()
#

#! format: off
varr = [
 [ "My gradients", "ν∇u[1]",  1.34348869323730468750e-04,  9.84732866287231445313e-01,  5.23545503616333007813e-01,  3.08209930764271777814e-01 ],
 [ "My gradients", "ν∇u[2]",  1.16317868232727050781e-01,  9.92088317871093750000e-01,  4.83800649642944335938e-01,  2.83350456014221541157e-01 ],
 [ "My gradients", "ν∇u[3]",  1.05845928192138671875e-03,  9.51775908470153808594e-01,  4.65474426746368408203e-01,  2.73615551085745090099e-01 ],
 [ "My gradients", "ν∇u[4]",  5.97668886184692382813e-02,  9.68048095703125000000e-01,  5.42618036270141601563e-01,  2.81570862027933854765e-01 ],
 [ "My gradients", "ν∇u[5]",  8.31030607223510742188e-02,  9.35931921005249023438e-01,  5.05405902862548828125e-01,  2.46073509972619536290e-01 ],
 [ "My gradients", "ν∇u[6]",  3.09681892395019531250e-02,  9.98341441154479980469e-01,  4.54375565052032470703e-01,  3.09461067853178561915e-01 ],
 [ "My gradients", "κ∇θ[1]",  8.47448110580444335938e-02,  9.94180679321289062500e-01,  5.27157366275787353516e-01,  2.92455951648181833313e-01 ],
 [ "My gradients", "κ∇θ[2]",  1.20514631271362304688e-02,  9.93527650833129882813e-01,  4.71063584089279174805e-01,  2.96449027197666359346e-01 ],
 [ "My gradients", "κ∇θ[3]",  8.14980268478393554688e-02,  9.55443382263183593750e-01,  5.05038917064666748047e-01,  2.77201022741208891187e-01 ],
 [    "My fields",   "u[1]",  4.31410233294131639781e-02,  9.97140933049696531754e-01,  4.62139750850942054861e-01,  3.23076684924287371725e-01 ],
 [    "My fields",   "u[2]",  1.01416659908237782872e-02,  9.14712023896926407218e-01,  4.76160523012988778913e-01,  2.71443440757963339038e-01 ],
 [    "My fields",   "θ[1]",  6.58965491052394547467e-02,  9.73216404386510802738e-01,  4.60007166313864512830e-01,  2.87310472114545079059e-01 ],
]
parr = [
 [ "My gradients", "ν∇u[1]",    16,     7,    16,     0 ],
 [ "My gradients", "ν∇u[2]",    16,     7,    16,     0 ],
 [ "My gradients", "ν∇u[3]",    16,     7,    16,     0 ],
 [ "My gradients", "ν∇u[4]",    16,     7,    16,     0 ],
 [ "My gradients", "ν∇u[5]",    16,     7,    16,     0 ],
 [ "My gradients", "ν∇u[6]",    16,     7,    16,     0 ],
 [ "My gradients", "κ∇θ[1]",    16,    16,    16,     0 ],
 [ "My gradients", "κ∇θ[2]",    16,    16,    16,     0 ],
 [ "My gradients", "κ∇θ[3]",    16,    16,    16,     0 ],
 [    "My fields",   "u[1]",    16,    16,    16,     0 ],
 [    "My fields",   "u[2]",    16,    16,    16,     0 ],
 [    "My fields",   "θ[1]",    16,    16,    16,     0 ],
]
#! format: on

max_stars_repo_name: leios/CLIMA

max_issues_repo_name: leios/CLIMA

max_forks_repo_name: leios/CLIMA


text:
# knn
module kNNClassifier
	import NearestNeighbors, Distances
	export kNNModel, kNN_train, kNN_exec

	using j4pr: countapp

	# kNNModel struct {T - type of neighbour, S - type of labels, U - type of distance}
	struct kNNModel{T<:Real, S<:Vector{<:Real}, U<:NearestNeighbors.NNTree, V<:Distances.Metric}
		k::T 						# k: T <: Real (if k is Int, do classic kNN, if k is Float, do range search)
		targets::S					# targets: Vector{T}, if T <: Float64, regression , classification otherwise
		treedata::U					# training data: a tree of sorts
		metric::V					# distance: Distances.Metric
		fsmooth::Function 				# smoother function::Function
		priors::Dict{Int,Float64}
	end

	# Aliases
	const kNNModelClassificationNN{T<:Int, S<:Vector{Int}} = kNNModel{T,S}
	const kNNModelClassificationR{T<:Float64, S<:Vector{Int}} = kNNModel{T,S}
	const kNNModelRegressionNN{T<:Int, S<:Vector{Float64}} = kNNModel{T,S}
	const kNNModelRegressionR{T<:Float64, S<:Vector{Float64}} = kNNModel{T,S}

	# Printers
	Base.show(io::IO, m::kNNModelClassificationNN) = print(io, "KNN model, classification, NN, k=$(m.k)")
	Base.show(io::IO, m::kNNModelClassificationR) = print(io, "KNN model, classification, range, r=$(m.k)")
	Base.show(io::IO, m::kNNModelRegressionNN) = print(io, "KNN model, regression, NN, k=$(m.k)")
	Base.show(io::IO, m::kNNModelRegressionR) = print(io, "KNN model, regression, range, r=$(m.k)")

	# Train method (classification)
	_parse_smooth_clf_(smooth::Symbol)::Function = begin
		# Check for smoothing option
		if smooth == :none fsmooth = no_smoother!
		elseif smooth == :ml fsmooth = ml_smoother!
		elseif smooth == :laplace fsmooth = laplace_smoother!
		elseif smooth == :mest fsmooth = mest_smoother!
		elseif smooth == :dist fsmooth = dist_smoother!
		else
			warn("Unrecognized posterior smoothing function, defaulting to :none")
			fsmooth = no_smoother!
		end
		return fsmooth
	end

	_parse_smooth_reg_(smooth::Symbol)::Function = begin
		# Check for smoothing option
		if smooth == :ml fsmooth = ml_smoother_r!
		elseif smooth == :dist fsmooth = dist_smoother_r!
		else
			warn("Unrecognized regression smoothing function, defaulting to `:ml`")
			fsmooth = ml_smoother_r!
		end
		return fsmooth
	end

	# Train methods (classification)
	kNN_train(x::Matrix{Float64}, y::Vector{Int}, k::Int; metric::Distances.Metric=Distances.Euclidean(), leafsize::Int=10, smooth::Symbol=:none) =
		return kNNModel(k, y, NearestNeighbors.KDTree(x, metric, leafsize=leafsize), metric, _parse_smooth_clf_(smooth), Dict(yi=>sum(yi.==y)/length(y) for yi in unique(y)))
	kNN_train(x::Matrix{Float64}, y::Vector{Int}, k::Float64; metric::Distances.Metric=Distances.Euclidean(), leafsize::Int=10, smooth::Symbol=:none) =
		return kNNModel(k, y, NearestNeighbors.BallTree(x, metric, leafsize=leafsize), metric, _parse_smooth_clf_(smooth), Dict(yi=>sum(yi.==y)/length(y) for yi in unique(y)))

	# Train methods (regression)
	kNN_train(x::Matrix{Float64}, y::Vector{Float64}, k::Int; metric::Distances.Metric=Distances.Euclidean(), leafsize::Int=10, smooth::Symbol=:ml) =
		return kNNModel(k, y, NearestNeighbors.KDTree(x, metric, leafsize=leafsize), metric, _parse_smooth_reg_(smooth), Dict{Int, Float64}())
	kNN_train(x::Matrix{Float64}, y::Vector{Float64}, k::Float64; metric::Distances.Metric=Distances.Euclidean(), leafsize::Int=10, smooth::Symbol=:ml) =
		return kNNModel(k, y, NearestNeighbors.BallTree(x, metric, leafsize=leafsize), metric, _parse_smooth_reg_(smooth), Dict{Int, Float64}())

	# Execution methods (classification)
	kNN_exec(m::kNNModelClassificationNN, x::Matrix{Float64}, classes::Vector{Int})=begin
		idx, d = NearestNeighbors.knn(m.treedata, x, m.k, true)
		out = zeros(Float64, length(classes), size(x,2))
		@inbounds for j = 1:size(x,2)
			m.fsmooth(view(out,:,j), m.targets[idx[j]], classes, d[j], m.priors)
		end
		return out
	end

	kNN_exec(m::kNNModelClassificationR, x::Matrix{Float64}, classes::Vector{Int})=begin
		idx = NearestNeighbors.inrange(m.treedata, x, m.k, false)
		out = zeros(Float64, length(classes), size(x,2))
		@inbounds for j = 1:size(x,2)
			d = zeros(Float64, size(idx[j]))
			@inbounds for i = 1: length(d)
				d[i] = Distances.evaluate(m.metric, m.treedata.data[idx[j][i]], x[:,j])
			end
			m.fsmooth(view(out,:,j), m.targets[idx[j]], classes, d, m.priors)
		end
		return out
	end

	# Execution methods (regression)
	kNN_exec(m::kNNModelRegressionNN, x::Matrix{Float64})=begin
		idx, d = NearestNeighbors.knn(m.treedata, x, m.k, true)
		out = zeros(Float64, 1, size(x,2))
		@inbounds for j = 1:size(x,2)
			m.fsmooth(view(out,:,j), m.targets[idx[j]], d[j])
		end
		return out
	end

	kNN_exec(m::kNNModelRegressionR, x::Matrix{Float64})=begin
		idx = NearestNeighbors.inrange(m.treedata, x, m.k, false)
		out = zeros(Float64, 1, size(x,2))
		@inbounds for j = 1:size(x,2)
			d = zeros(Float64, length(idx[j]))
			@inbounds for i = 1: length(d)
				d[i] = Distances.evaluate(m.metric, m.treedata.data[idx[j][i]], x[:,j])
			end
			m.fsmooth(view(out,:,j), m.targets[idx[j]], d)
		end
		return out
	end

	# Posterior smoothers (classification); parameters: labels of neighbors, unique labels in dataset and distances
	# All smoothers return a vector yu in length
	no_smoother!(out::T where T<:AbstractVector{Float64}, y::Vector{Int}, yu::Vector{Int}, d=zeros(length(y)), priors=Dict{Int,Float64}()) = begin
		n = countapp(y,yu)
		@inbounds out[findmax(n)[2]] = 1.0
		return out
	end
	ml_smoother!(out::T where T<:AbstractVector{Float64}, y::Vector{Int}, yu::Vector{Int}, d=zeros(length(y)), priors=Dict{Int,Float64}()) = begin
		n = countapp(y,yu)
		@fastmath out[:] = n/length(y)
		return out
	end
	laplace_smoother!(out::T where T<:AbstractVector{Float64}, y::Vector{Int}, yu::Vector{Int}, d=zeros(length(y)), priors=Dict{Int,Float64}()) = begin
		n = countapp(y,yu)
		@fastmath out[:] = (n+1.0)/(length(y)+length(yu))
		return out
	end

	mest_smoother!(out::T where T<:AbstractVector{Float64}, y::Vector{Int}, yu::Vector{Int}, d=zeros(length(y)), priors=Dict{Int,Float64}()) = begin
		n = countapp(y,yu)
		p = [priors[i] for i in yu]
		m = 10.0
		@fastmath out[:] = (n+m*p)/(length(y)+m)
		return out
	end

	dist_smoother!(out::T where T<:AbstractVector{Float64}, y::Vector{Int}, yu::Vector{Int}, d::Vector{Float64}, priors=Dict{Int,Float64}()) = begin
		dc = zeros(Float64, size(yu))
		m = length(yu)
		n = length(y)
		for j in 1:m
			for i in 1:n
				@inbounds if y[i] == yu[j]
					@fastmath dc[j] += 1/(eps()+d[i])
				end
			end
		end
		out[:] = dc / sum(dc)
		return out
	end

	# Posterior smoothers (regression)
	ml_smoother_r!(out::T where T<:AbstractVector{Float64}, y::Vector{Float64}, d=zeros(length(y))) = begin
		@fastmath out[1] = mean(y)
		return out
	end
	dist_smoother_r!(out::T where T<:AbstractVector{Float64}, y::Vector{Float64}, d::Vector{Float64}) = begin
		t1 = @. y*(1/(eps()+d))
		t2 = @. 1/(eps()+d)
		out[1] = sum(t1)/sum(t2)
		return out
	end
end

##########################
# FunctionCell Interface #
##########################
"""
	knn(k=1 [;kwargs])

Constructs an untrained cell that when piped data inside, returns a knn classifier trained
function cell based on the input data and labels.

# Arguments
  * `k::Int` represents the number of neighbours to be considered
  * `k::Float64=1.0` represents a range in which neighboring samples are to be searched (if no neighbours are found,
  equal class posteriors are generated)

# Keyword arguments
  * `metric::Distances.Metric` metric from `Distances.jl` which specifies how distances are calculated between samples (default `Distances.Euclidean()`)
  * `leafsize::Int` Size of the leafs for the tree structures holding the training data (default `10`)
  * `smooth::Symbol` method for smoothing class posteriors; available: `:none`, `:ml`, `:laplace`,`mest (m=10)` and `:dist` (default `:none`)

The posterior smoothing is as such: `:none` return probability `1.0` for the closest class, `:ml` returns class-wise mean probabilities or mean of
neighbouring values for regression, `:laplace` is similar to `:ml` but adds some constants, `:mest` regularizes towards class priors and `:dist` used
the distances from sample to neighbours to obtain the posterior in the case of classification while for regression it weighs neighbor values in inverse
proportion to their corresponding distance.

For more information on posterior smoothing:
	[1] L. Kuncheva "Combining Pattern Classifiers 2'nd Ed." 2014, ISBN 978-1-118-31523-1

Read the `NearestNeighbours.jl` and `Distances.jl` documentation for more information.
"""
knn(k::Real=1; kwargs...) = FunctionCell(knn, (k,), ModelProperties(), kwtitle(k isa Int ? "$k-NN classifier" : "$k-range-NN classifier", kwargs); kwargs...)

############################
# DataCell/Array Interface #
############################
"""
	knn(x, k=1 [;kwargs])

Trains a knn classifcation model that using the data `x` and `k` neighbors.
"""
# Training
knn(x::T where T<:CellDataL, k::Real=1; kwargs...) = knn((getx!(x), gety(x)), k; kwargs...)
knn(x::Tuple{T,S} where T<:AbstractVector where S<:AbstractVector, k::Real=1; kwargs...) = knn((mat(x[1], LearnBase.ObsDim.Constant{2}()), x[2]), k; kwargs...)
knn(x::Tuple{T,S} where T<:AbstractMatrix where S<:AbstractVector, k::Real=1; kwargs...) = begin

	@assert nobs(x[1]) == nobs(x[2]) "[knn] Expected $(nobs(x[1])) labels/values, got $(nobs(x[2]))."

	# Transform labels first
	enc = labelencn(x[2])
	yenc = label2ind.(x[2],enc)

	# Train model
	knndata = kNNClassifier.kNN_train(getobs(x[1]), yenc, k; kwargs...)

	# Build model properties
	modelprops = ModelProperties(nvars(x[1]), length(enc.label), enc)

	FunctionCell(knn, Model(knndata, modelprops), kwtitle(k isa Int ? "$k-NN classifier" : "$k-range-NN classifier", kwargs))

end

# Execution
knn(x::T where T<:CellData, model::Model{<:kNNClassifier.kNNModel}) =
	datacell(knn(getx!(x), model), gety(x))

knn(x::T where T<:AbstractVector, model::Model{<:kNNClassifier.kNNModel}) =
	knn(mat(x, LearnBase.ObsDim.Constant{2}()), model)

knn(x::T where T<:AbstractMatrix, model::Model{<:kNNClassifier.kNNModel}) =
	kNNClassifier.kNN_exec(model.data, getobs(x), collect(1:model.properties.odim))

##########################
# FunctionCell Interface #
##########################
"""
	knnr(k=1 [;kwargs])

Constructs an untrained cell that when piped data inside, returns a knn trained regressor
function cell based on the input data and labels.

# Arguments
  * `k::Int` represents the number of neighbours to be considered
  * `k::Float64=1.0` represents a range in which neighboring samples are to be searched (if no neighbours are found,
  equal class posteriors are generated)

# Keyword arguments
  * `metric::Distances.Metric` metric from `Distances.jl` which specifies how distances are calculated between samples (default `Distances.Euclidean()`)
  * `leafsize::Int` Size of the leafs for the tree structures holding the training data (default `10`)
  * `smooth::Symbol` method for smoothing class posteriors; available: `:ml`, `:dist` (default `:ml`)

Read the `NearestNeighbours.jl` and `Distances.jl` documentation for more information.
"""
knnr(k::Real=1; kwargs...) = FunctionCell(knnr, (k,), ModelProperties(), kwtitle(k isa Int ? "$k-NN regressor" : "$k-range-NN regressor", kwargs); kwargs...)

############################
# DataCell/Array Interface #
############################
"""
	knnr(x, k=1 [;kwargs])

Trains a knn regression model that using the data `x` and `k` neighbors.
"""
# Training
knnr(x::T where T<:CellDataL, k::Real=1; kwargs...) = knnr((getx!(x), gety(x)), k; kwargs...)
knnr(x::Tuple{T,S} where T<:AbstractVector where S<:AbstractVector, k::Real=1; kwargs...) =
	knnr((mat(x[1], LearnBase.ObsDim.Constant{2}()), x[2]), k; kwargs...)
knnr(x::Tuple{T,S} where T<:AbstractMatrix where S<:AbstractVector, k::Real=1; kwargs...) = begin

	@assert nobs(x[1]) == nobs(x[2]) "[knnr] Expected $(nobs(x[1])) labels/values, got $(nobs(x[2]))."

	# Train model
	knndata = kNNClassifier.kNN_train(getobs(x[1]), getobs(x[2]), k; kwargs...)

	FunctionCell(knnr, Model(knndata, ModelProperties(nvars(x[1]),1)), kwtitle(k isa Int ? "$k-NN regressor" : "$k-range-NN regressor", kwargs))

end

# Execution
knnr(x::T where T<:CellData, model::Model{<:kNNClassifier.kNNModel}) = datacell(knnr(getx!(x), model), gety(x))
knnr(x::T where T<:AbstractVector, model::Model{<:kNNClassifier.kNNModel}) = knnr(mat(x, LearnBase.ObsDim.Constant{2}()), model)
knnr(x::T where T<:AbstractMatrix, model::Model{<:kNNClassifier.kNNModel}) = kNNClassifier.kNN_exec(model.data, getobs(x))


text:
module Runtests

# Run all tests with SilentNoPlot (so not plots)

import ModiaResult
using Test

@testset "Test ModiaResult/test" begin
    ModiaResult.usePlotPackage("SilentNoPlot")
    include("include_all.jl")
    ModiaResult.usePreviousPlotPackage()
end

end

max_stars_repo_name: ModiaSim/ModiaResult

max_issues_repo_name: ModiaSim/ModiaResult

max_forks_repo_name: ModiaSim/ModiaResult


text:
function build_time_graph()
    g = time_graph(Int, 1)
    add_edge!(g, 0, 1)
    add_edge!(g, 1, 2)
    add_edge!(g, 2, 3)
    add_edge!(g, 2, 4)
    add_edge!(g, 3, 4)
    add_edge!(g, 3, 5)
    add_edge!(g, 0, 2)
    add_edge!(g, 1, 0)
    add_edge!(g, 3, 1)
    add_edge!(g, 4, 0)
    g
end

g = build_time_graph()

p = shortest_path(g, 0, 5, verbose = true)

p2 = Path([make_node(g,0), make_node(g,2), make_node(g,3), make_node(g,5)])

@test p == p2

@test shortest_distance(g,  0,  5) == 3


text:
export foo, bar

greet() = print("I love Anna!")

f1() = print("it works!")

foo(x::T, y::T) where T <: Real = x + y - 5
bar(z::Float64) = foo(sqrt(z), z)


text:
function check(numbers, idx)
    for i in idx - 25:idx - 2
        for j in i:idx - 1
            if numbers[i] + numbers[j] == numbers[idx]
                return true
            end
        end
    end
    return false
end

function run(s)
    numbers = readlines(IOBuffer(s))
    numbers = [parse(Int, n) for n in numbers]
    # find invalid number
    invalid_number = 0
    for idx in 26:length(numbers)
        s = numbers[idx]
        if !check(numbers, idx)
            invalid_number = numbers[idx]
            break
        end
    end

    # find matching sequence
    i = 1
    j = 2
    s = sum(numbers[i:j])
    while s != invalid_number
        if s < invalid_number
            j += 1
            s += numbers[j]
        else
            s -= numbers[i]
            i += 1
        end
    end
    return minimum(numbers[i:j]) + maximum(numbers[i:j])
end

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

function main()
    res, time, memory = @timed run(ARGS[1])
    println("_duration:$(time * 1000)")
    println(res)
end

main()


text:
using JLD3
using Test

@testset "JLD3.jl" begin
    # Write your tests here.
end


text:
export OceanBC,
    VelocityBC,
    VelocityDragBC,
    TemperatureBC,
    Impenetrable,
    Penetrable,
    NoSlip,
    FreeSlip,
    KinematicStress,
    Insulating,
    TemperatureFlux

using StaticArrays

using ..BalanceLaws
using ..DGMethods.NumericalFluxes

"""
    OceanBC(velocity    = Impenetrable(NoSlip())
            temperature = Insulating())

The standard boundary condition for OceanModel. The default options imply a "no flux" boundary condition.
"""
Base.@kwdef struct OceanBC{M, T}
    velocity::M = Impenetrable(NoSlip())
    temperature::T = Insulating()
end

abstract type VelocityBC end
abstract type VelocityDragBC end
abstract type TemperatureBC end

"""
    Impenetrable(drag::VelocityDragBC) :: VelocityBC

Defines an impenetrable wall model for velocity. This implies:
  - no flow in the direction normal to the boundary, and
  - flow parallel to the boundary is subject to the `drag` condition.
"""
struct Impenetrable{D <: VelocityDragBC} <: VelocityBC
    drag::D
end

"""
    Penetrable(drag::VelocityDragBC) :: VelocityBC

Defines an penetrable wall model for velocity. This implies:
  - no constraint on flow in the direction normal to the boundary, and
  - flow parallel to the boundary is subject to the `drag` condition.
"""
struct Penetrable{D <: VelocityDragBC} <: VelocityBC
    drag::D
end

"""
    NoSlip() :: VelocityDragBC

Zero velocity at the boundary.
"""
struct NoSlip <: VelocityDragBC end

"""
    FreeSlip() :: VelocityDragBC

No surface drag on velocity parallel to the boundary.
"""
struct FreeSlip <: VelocityDragBC end

"""
    KinematicStress(stress) :: VelocityDragBC

Applies the specified kinematic stress on velocity normal to the boundary.
Prescribe the net inward kinematic stress across the boundary by `stress`,
a function with signature `stress(problem, state, aux, t)`, returning the flux (in m²/s²).
"""
struct KinematicStress{S} <: VelocityDragBC
    stress::S

    function KinematicStress(stress::S = nothing) where {S}
        new{S}(stress)
    end
end

kinematic_stress(problem, y, ρ₀) = @SVector [0, 0] # fallback for generic problems
kinematic_stress(problem, y, ρ₀, ::Nothing) = kinematic_stress(problem, y, ρ₀)
kinematic_stress(problem, y, ρ₀, stress) = stress(y)

"""
    Insulating() :: TemperatureBC

No temperature flux across the boundary
"""
struct Insulating <: TemperatureBC end

"""
    TemperatureFlux(flux) :: TemperatureBC

Prescribe the net inward temperature flux across the boundary by `flux`,
a function with signature `flux(problem, state, aux, t)`, returning the flux (in m⋅K/s).
"""
struct TemperatureFlux <: TemperatureBC end

# these functions just trim off the extra arguments
function _ocean_boundary_state!(
    nf::Union{NumericalFluxFirstOrder, NumericalFluxGradient},
    bc,
    ocean,
    Q⁺,
    A⁺,
    n,
    Q⁻,
    A⁻,
    t,
    _...,
)
    return ocean_boundary_state!(nf, bc, ocean, Q⁺, A⁺, n, Q⁻, A⁻, t)
end

function _ocean_boundary_state!(
    nf::NumericalFluxSecondOrder,
    bc,
    ocean,
    Q⁺,
    D⁺,
    A⁺,
    n,
    Q⁻,
    D⁻,
    A⁻,
    t,
    _...,
)
    return ocean_boundary_state!(nf, bc, ocean, Q⁺, D⁺, A⁺, n, Q⁻, D⁻, A⁻, t)
end

max_forks_repo_name: jm-c/CLIMA


text:
using JSON, Compat, Mmap, HDF5, Dates, DataStructures, YAML

"""
    load_data(filename)

Load data file
"""
# add functionality to read new Auspex data files (format on develop branch)

function load_data(filename::AbstractString; load_var = false)
    # maintain original data format when variance is not avilable for compatibility with plotters and other data analysis. Use variance for tomography
    if load_var
        datasets    = Dict{String, Dict{String, Array{Any}}}()
    else
        datasets    = Dict{String, Array{Any}}()
    end
    descriptors = Dict{String,Any}()
    # per qubit
    for group in readdir(filename)
        if load_var
            datasets[group] = Dict{String, Any}()
        end
        datafile = filter(x-> occursin("data.dat", x), readdir(joinpath(filename, group)))
        varfile = filter(x-> occursin("variance.dat", x), readdir(joinpath(filename, group)))
        metafile = filter(x-> occursin("meta", x), readdir(joinpath(filename, group)))
        metadata = JSON.parsefile(joinpath(filename, group, metafile[1]), dicttype=DataStructures.OrderedDict)
        open(joinpath(filename, group, datafile[1]), "r") do op
            dims = length(metadata["shape"])>1 ? tuple(metadata["shape"]...) : tuple(metadata["shape"][1],1)
            data = Mmap.mmap(op, Array{Complex{Float64},1}, prod(dims))
            data = reshape(data, reverse(dims))
            if load_var
                datasets[group]["data"] = data
            else
                datasets[group] = data
            end
        end
        if load_var
            if ~isempty(varfile)
                open(joinpath(filename, group, varfile[1]), "r") do op
                    dims = length(metadata["shape"])>1 ? tuple(metadata["shape"]...) : tuple(metadata["shape"][1],1)
                    data = Mmap.mmap(op, Array{Complex{Float64},1}, prod(dims))
                    data = reshape(data, reverse(dims))
                    datasets[group]["variance"] = data
                end
            else
              @warn("Variance not available for dataset: $(group)")
            end
        end
        descriptors[group] = metadata
    end
    return datasets, descriptors
end

"""
    load_data(datapath, filenum; subdir, auspex)

Search file number filenum in folder datapath/subdir. Subdir default is current date in "yymmdd"
auspex: boolean for file format (default = true)
"""
function load_data(datapath::AbstractString, filenum::Int, subdir=Dates.format(Dates.today(),"yymmdd"), h5=false; load_var=false)
  # optionally, search for filenum instead of filename
  # search in a subdirectory with today's date, if not specified
  #regexpr = auspex? r"(\d{4})(.h5)" : r"(\d+)(_\w+)(.h5)"
  ext = h5 ? ".h5" : ".auspex"
  regexpr = Regex("(.+)($(lpad(filenum, 4, "0")))(\\$ext)")
  files = filter(x -> occursin(regexpr, x), readdir(joinpath(datapath, subdir)))
  if length(files) == 0
    println("No matching files found!")
    return
  elseif length(files) > 1
    println("More than one matching file found!")
    return
  end
  filename = joinpath(datapath, subdir, files[1])
  h5 ? load_h5_data(filename) : load_data(filename, load_var=load_var)
end

## TO BE UPDATED
# """
# 	load_latest_data()
#     load_latest_data(ind_from_last)
#
# Load latest data file, or ind_from_last files before the latest.
# """
# function load_latest_data(logpath::AbstractString, ind_from_last::Int = 0)
# 	df = CSV.read(logpath, delim = "\t"; use_mmap = false)
# 	@assert ind_from_last < length(df.columns[1]) "Not enough data files in the history."
# 	filename = df.columns[1][end-ind_from_last]
# 	load_auspex_data(filename)
# end

"""
    load_h5_data(filename)

Load hdf5 data and descriptor in legacy auspex format

FILE STRUCTURE:
file.hdf5:
   header/ (with yml settings, if saved)
   group1/
      data/
         axis_1 (contains all tuples visited on this axis)
         axis_2 (contains all tuples visited on this axis)
         quantity_1 (the measured values)
         quantity_2 (the measured values)
      descriptor (contains references to the axes datasets below)
      axis_1 (contains the values of the axis  see note #1)
      axis_2 (contains the values of the axis - see note #1)
   group2/
      ... same structure as above ...

Note #1: Any axes referenced by the descriptor can themselves contain
         references to multiple subaxes. This indicates that the axis
         is "unstructured", i.e. that multiple coordinates where changed
         simultaneously.
"""
function load_h5_data(filename::AbstractString)

    datasets    = Dict{String, Dict{String, Vector{Any}}}()
    descriptors = Dict{String, Vector{Dict{String,Any}}}()
	header =   Dict{String, Any}()
	#info on file, such as name, ...?

	header["filename"] = filename
    h5open(filename, "r") do f
        # Find all of the group names, which will correspond to qubits when
        # using exp_factory inside of auspex. The default group name is simply
        # "main".
        group_names = names(f)
		#load measurement/instrument settings
		if "header" in group_names
			header["settings"] = YAML.load(read(attrs(f["header"])["settings"]))
		end

        for group_name in filter!(name -> name != "header", group_names)
            g = f[group_name]

            # Read in the descriptor
            desc_refs = read(g["descriptor"])

            desc = Vector{Dict{String,Any}}()
            for desc_ref in desc_refs
                axis = Dict{String,Any}()
                if reinterpret(Bool, read(attrs(g[desc_ref])["unstructured"]))
                    axis_refs      = read(g[desc_ref])
                    axis["points"] = [read(g[axis_ref]) for axis_ref in axis_refs]
                    axis["name"]   = [read(attrs(g[axis_ref])["name"]) for axis_ref in axis_refs]
                    axis["unit"]   = [read(attrs(g[axis_ref])["unit"]) for axis_ref in axis_refs]
                else
                    axis["points"] = read(g[desc_ref])
                    axis["name"]   = read(attrs(g[desc_ref])["name"])
                    axis["unit"]   = read(attrs(g[desc_ref])["unit"])
                end
                # See if we need to get metadata:
                if "metadata" in names(attrs(g[desc_ref]))
                    meta_ref = read(attrs(g[desc_ref])["metadata"])
                    meta_enum_ref = read(attrs(g[desc_ref])["metadata_enum"])
                    axis["metadata"] = [read(g[meta_enum_ref])[m+1] for m in read(g[meta_ref])]
                end

                push!(desc, axis)
            end

            descriptors[group_name] = desc

            # Read in the actual data
            datacol_names = names(g["data"])
            datasets[group_name] = Dict{String, Vector{Any}}()
            for datacol_name in datacol_names
              dataset = read(g["data"][datacol_name])
              # heuristic to read complex numbers stored as HDF5Compounds
              datasets[group_name][datacol_name] = convert_dataset(dataset)
            end
        end
    end
    datasets, descriptors, header
end

convert_dataset(ds) = ds
# heuristic to read complex numbers stored as HDF5Compounds
function convert_dataset(ds::Vector{NamedTuple{(:r, :i), Tuple{Float32, Float32}}})
  return [complex(d.data[1], d.data[2]) for d in ds]
end


text:
@testset "Distribution types" begin
    # Test the types and the sizes of the distributions
    @testset "Distributions D1Q3" begin
        N = 10
        typelist = [Float16 Float32 Float64 Real]
        for type in typelist
            test_distribution = JuSwalbe.DistributionD1Q3(zeros(type,N),zeros(type,N),zeros(type,N))
            test_dict = struct2dict(test_distribution)
            @test isa(test_distribution, JuSwalbe.DistributionD1Q3)
            for key in keys(test_dict)
                @test isa(test_dict[key], Array{type,1})
                @test length(test_dict[key]) == N
            end
        end
    end
    @testset "Distributions D2Q9" begin
        N, M = 10, 14
        typelist = [Float16 Float32 Float64 Real]
        for type in typelist
            test_distribution = JuSwalbe.DistributionD2Q9(zeros(type,(N,M)),
                                                          zeros(type,(N,M)),
                                                          zeros(type,(N,M)),
                                                          zeros(type,(N,M)),
                                                          zeros(type,(N,M)),
                                                          zeros(type,(N,M)),
                                                          zeros(type,(N,M)),
                                                          zeros(type,(N,M)),
                                                          zeros(type,(N,M)))
            test_dict = struct2dict(test_distribution)
            @test isa(test_distribution, JuSwalbe.DistributionD2Q9)
            for key in keys(test_dict)
                @test isa(test_dict[key], Array{type,2})
                @test size(test_dict[key]) == (N,M)
            end
        end
    end
end

max_stars_repo_name: Zitzeronion/JuSwalbe

max_issues_repo_name: Zitzeronion/JuSwalbe


text:
markable(x) = false
markable(x::T) where {T<:Base.IEEEFloat} = true
markable(x::T) where {N, T<:MarkableFloat{N}} = true

unmarked(x::T) where {T<:Base.IEEEFloat} = true
unmarked(x::T) where {N, T<:MarkableFloat{N}} = false
unmarked(x::T) where {T<:MarkableFloat{0}} = true

marked(x::T) where {T<:Base.IEEEFloat} = false
marked(x::T) where {N, T<:MarkableFloat{N}} = true
marked(x::T) where {T<:MarkableFloat{0}} = false

has0marks(x::T) where {T<:Base.IEEEFloat} = true
has0or1marks(x::T) where {T<:Base.IEEEFloat} = true
has1mark(x::T) where {T<:Base.IEEEFloat} = false
has2marks(x::T) where {T<:Base.IEEEFloat} = false
has3marks(x::T) where {T<:Base.IEEEFloat} = false
has1or2marks(x::T) where {T<:Base.IEEEFloat} = false
has1or3marks(x::T) where {T<:Base.IEEEFloat} = false
has2or3marks(x::T) where {T<:Base.IEEEFloat} = false

has0marks(x::T) where {T<:MarkableFloat{0}} = true
has0marks(x::T) where {N, T<:MarkableFloat{N}} = false
has1mark(x::T) where {T<:MarkableFloat{1}} = true
has1mark(x::T) where {N, T<:MarkableFloat{N}} = false
has2marks(x::T) where {T<:MarkableFloat{2}} = true
has2marks(x::T) where {N, T<:MarkableFloat{N}} = false
has3marks(x::T) where {T<:MarkableFloat{3}} = true
has3marks(x::T) where {N, T<:MarkableFloat{N}} = false

has0or1marks(x::T) where {T<:MarkableFloat{0}} = true
has0or1marks(x::T) where {T<:MarkableFloat{1}} = true
has0or1marks(x::T) where {N, T<:MarkableFloat{N}} = false

has1or2marks(x::T) where {T<:MarkableFloat{1}} = true
has1or2marks(x::T) where {T<:MarkableFloat{2}} = true
has1or2marks(x::T) where {N, T<:MarkableFloat{N}} = false

has0or2marks(x::T) where {T<:MarkableFloat{0}} = true
has0or2marks(x::T) where {T<:MarkableFloat{2}} = true
has0or2marks(x::T) where {N, T<:MarkableFloat{N}} = false

has1or3marks(x::T) where {T<:MarkableFloat{1}} = true
has1or3marks(x::T) where {T<:MarkableFloat{3}} = true
has1or3marks(x::T) where {N, T<:MarkableFloat{N}} = false

has2or3marks(x::T) where {T<:MarkableFloat{2}} = true
has2or3marks(x::T) where {T<:MarkableFloat{3}} = true
has2or3marks(x::T) where {N, T<:MarkableFloat{N}} = false


text:
module TMvNormals
using Distributions, Statistics, LinearAlgebra
import Distributions: cdf, pdf, mean, cov, rand
using Parameters

include("mvnormal.jl")

"""
Truncated Multi-variate Normal distribution with lower limits [a] and upper limits [b].

To avoid limit issues, single- and un-truncated dimensions are handled as being truncated
at μ±10σ.  This shouldn't affect numerical results, as the code is also designed to error
out if the cdf routine is too inaccurate.
"""
struct TMvNormal <: ContinuousMultivariateDistribution
    μ::AbstractVector
    Σ::AbstractMatrix
    a::AbstractVector
    b::AbstractVector
    𝒩::AbstractMvNormal
    α::Float64
    K::Int
end

Base.show(io::IO, d::TMvNormal) =
    Distributions.show_multline(io, d, [(:𝒩, d.𝒩), (:lower, round.(d.a; digits=3)), (:upper, round.(d.b; digits=3))])

Base.broadcastable(d::TMvNormal) = Ref(d)
# function Base.length(::TMvNormal)
#     return 1
# end

# function Base.iterate(x::TMvNormal, state=1)
#     if state == 1
#         return (x, 2)
#     end
# end

"""
More convenient way to initialize the distribution with μ, Σ, a, and b
"""
function _TMvNormal_ab(μ::AbstractVector, Σ::AbstractMatrix, a::AbstractVector, b::AbstractVector)
    Σ = (Σ + Σ') / 2
    𝒩 = MvNormal(μ, Σ)
    a = ifelse.(
        isfinite.(a),
        a,
        μ - 10 * sqrt.(diag(Σ))
    )
    b = ifelse.(
        isfinite.(b),
        b,
        μ + 10 * sqrt.(diag(Σ))
    )
    return TMvNormal(
        μ,Σ,a,b,𝒩,
        cdf(𝒩, a, b),
        length(μ)
)
end

function TMvNormal(μ::AbstractVector{Float64}, Σ::AbstractMatrix, vecs::Vararg{AbstractVector})
    if length(μ) == length(vecs) > 2
        # Case 1: Sequence of lower, upper bounds
        return _TMvNormal_ab(μ, Σ, collect(map(first, vecs)), collect(map(last, vecs)))
    elseif length(vecs) == 2
        if all(map(issorted, vecs)) && all(map(==(2) ∘ length, vecs))
            # Case 2: 2d distribution, so it's ambiguous whether we have a lower and upper bound
            return _TMvNormal_ab(μ, Σ, collect(map(first, vecs)), collect(map(last, vecs)))
        elseif all(vecs[1] .<= vecs[2])
            # Case 3: 2d distribution and it's unambiguous
            return _TMvNormal_ab(μ, Σ, vecs[1], vecs[2])
        end
    else
        error("TMvNormal: invalid arguments")
    end
end

"""
More convenient way to initialize the distribution with 𝒩, a, and b
"""
function TMvNormal(𝒩::MvNormal, a::AbstractVector, b::AbstractVector)
    a = ifelse.(
        isfinite.(a),
        a,
        𝒩.μ - 10 * sqrt.(diag(𝒩.Σ))
    )
    b = ifelse.(
        isfinite.(b),
        b,
        𝒩.μ + 10 * sqrt.(diag(𝒩.Σ))
    )
    return TMvNormal(
        𝒩.μ, 𝒩.Σ, a, b, 𝒩,
        cdf(𝒩, a, b),
        length(𝒩.μ)
    )
end

"""
Worst possible way to draw samples from the distribution
"""
function rand(d::TMvNormal, N::Int)
    @unpack a, b, K, 𝒩 = d
    x = Array{Float64}(undef, (N, K))
    for i = 1:N
        _trial = 0
        while true
            _trial += 1
            _x = rand(𝒩)
            if all(a .< _x .< b)
                x[i, :] = _x
                break
            elseif _trial > 100
                error("Failed to generate random sample")
        end
    end
    end
    return x
end

function cov2cor(Σ::AbstractMatrix)
    Σ = (Σ + Σ') / 2
    D = Diagonal((diag(Σ).^-0.5))
    return D * Σ * D
end

function cdf(𝒩::MvNormal, a::AbstractVector, b::AbstractVector)
    K = length(𝒩.μ)
    if K == 1
        𝒩 = Normal(𝒩.μ[1], √𝒩.Σ[1,1])
        return ifelse(b[1] == Inf, 1, cdf(𝒩, b[1])) - ifelse(a[1] == -Inf, 0, cdf(𝒩, a[1]))
    end

    if !all(a .<= b)
        error("a must be less than or equal to b")
    end

    val, err = qsimvnv(
        𝒩.Σ,
        a - 𝒩.μ,
        b - 𝒩.μ
    )
    if val == 0 || err < val / 100
        return val
    else
        error("Relative error exceeds 1%")
    end
end

function cdf(𝒩::MvNormal, b::AbstractVector)
    a = fill(-Inf, size(b))
    return cdf(𝒩, a, b)
end

function pdf(d::TMvNormal, x::AbstractVector)
    @unpack a, b, α, 𝒩 = d
    if all(a .<= x .<= b)
        return pdf(𝒩, x) / α
    else
        return 0
    end
end

"""
Calculate the marginal univariate density of a truncated multivariate normal.
Formula taken from Cartinhour (1990)
https://doi.org/10.1080/03610929008830197
"""
function pdf(d::TMvNormal, x::Number, dim::Number=1)
    @unpack a, b, K, α, μ, 𝒩, Σ = d
    idx = filter(!=(dim), 1:K)

    if !(a[dim] <= x <= b[dim])
        return 0
    end

    # get A̲₁, added a Symmetric call to handle small numerical errors
    Σ̲ = Symmetric(inv(inv(Matrix(Σ))[idx, idx]))
    a̲ = a[idx]
    b̲ = b[idx]
    σ̲ = Σ[idx, dim]
    # display(a)
    S(_x) = α^-1 * cdf(MvNormal(μ[idx] + (_x - μ[dim]) * σ̲ / Σ[dim,dim], Σ̲), a̲, b̲)
    return S(x) * pdf(Normal(μ[dim], √Σ[dim,dim]), x)
end

"""
Calculate the marginal bivariate density of a truncated multivariate normal.
Formula taken from Manjunath Wilhelm (2021)
https://doi.org/10.35566/JBDS%2FV1N1%2FP2
"""
function pdf(d::TMvNormal, x::AbstractVector, margin::AbstractVector)
    @unpack α, K, μ, Σ, a, b = d
    margin = unique(margin)

    if length(margin) == 1
        return pdf(d, x[1], margin[1])
    elseif length(margin) != 2
        error("margin must be a vector of length 1 or 2, i got $(length(margin))")
    end

    q, r = margin

    if !all(a[[q, r]] .<= x .<= b[[q, r]])
        return 0
    end

    idx = filter(!∈(margin), 1:K)

    D = sqrt.(diag(Σ))
    R = cov2cor(Σ)

    # Manjunath Wilhelm (2021) says to use a z-transformed normal distribution for the marginal, but their code in `tmvtnorm.R` does not
    # Simulations agree with code, not paper
    ϕ = pdf(MvNormal(μ[[q, r]], Σ[[q, r], [q, r]]), x)

    # Multivariate regression coefficients
    β(s, q, r) = (R[s, q] - R[q, r] * R[s,r]) / (1 - R[q, r]^2)

    # Partial correlation coefficients
    function ρ(i, j, control::AbstractVector=[])
        if control == []
            return R[i,j]
        end
        (ρ(i, j, control[2:end]) - ρ(i, control[1], control[2:end]) * ρ(j, control[1], control[2:end])) / √((1 - ρ(j, control[1], control[2:end])^2) * (1 - ρ(i, control[1], control[2:end])^2))
    end

    function ρ(i, j, control::Int)
        ρ(i, j, [control])
    end

    # ρ(s, q, r) = β(s, q, r) * √(1 - R[q, r]^2) / √(1 - R[s, r]^2)

    a = (a - μ) ./ D
    b = (b - μ) ./ D

    c = (x - μ[[q, r]]) ./ D[[q, r]]

    A(_q, _r, _s) = (a[_s] - β(_s, _q, _r) * c[1] - β(_s, _r, _q) * c[2]) / √((1 - R[_s, _q]^2) * (1 - ρ(_s, _r, _q)^2))
    B(_q, _r, _s) = (b[_s] - β(_s, _q, _r) * c[1] - β(_s, _r, _q) * c[2]) / √((1 - R[_s, _q]^2) * (1 - ρ(_s, _r, _q)^2))

    R₂ = Array{Float64}(undef, K, K)
    for i in idx,j in idx
        R₂[i,j] = ρ(i, j, [q,r])
    end
    R₂ = R₂[idx,idx]

    if K - 2 > 0
        Φᵈ⁻² = cdf(MvNormal(zeros(K - 2), R₂), A.(q, r, idx), B.(q, r, idx))
#         display(A.(q, r, idx))
#         display(B.(q, r, idx))
#         display(R₂)
    else
        Φᵈ⁻² = 1
    end
    return α^-1 * ϕ * Φᵈ⁻²
end

function mean(d::TMvNormal, dim::Int)
    @unpack μ, K, Σ, a, b = d
    μ[dim] + reduce(
        +,
        Σ[1:K, dim] .* (pdf.(d, a, 1:K) - pdf.(d, b, 1:K))
    )
end

function mean(d::TMvNormal)
    @unpack K = d
    return mean.(d, 1:K)
end

function cov(d::TMvNormal, dims::AbstractVector{Int})

    if length(dims) != 2
        error("covariance matrix must be two dimensional")
    end

    @unpack K, μ, Σ, a, b = d

    D = sqrt.(diag(Σ))
    Σ
    a = (a - μ) # ./ D
    b = (b - μ) # ./ D
    d̂ = TMvNormal(zeros(K), Σ, a, b)
    μ̂ = mean(d̂)

    i, j = dims

    first_sum = reduce(
        +,
        map(
            k -> Σ[i,k] / Σ[k,k] * Σ[j,k] *
            (
                ifelse(isfinite(a[k]), a[k] * pdf(d̂, a[k], k), 0) -
                ifelse(isfinite(b[k]), b[k] * pdf(d̂, b[k], k), 0)
            ),
            1:K
        )
    )

    function inner_summand(q, k)
        pdf_term1 = pdf(d̂, [a[k], a[q]], [k, q]) - pdf(d̂, [a[k], b[q]], [k, q])
        pdf_term2 = pdf(d̂, [b[k], a[q]], [k, q]) - pdf(d̂, [b[k], b[q]], [k, q])
        (Σ[j, q] - Σ[k,q] * Σ[j,k] / Σ[k,k]) * (pdf_term1 - pdf_term2)
    end

    second_sum = reduce(
        +,
        map(
            k -> reduce(
                +,
                map(
                    q -> Σ[i, k] * inner_summand(q, k),
                    filter(!=(k), 1:K)
                )
            ),
            1:K
        )
    )

    return (Σ[i,j] + first_sum + second_sum - μ̂[i] * μ̂[j])
end

function cov(d::TMvNormal)
    @unpack K = d
    C = Array{Float64}(undef, K, K)
    for i = 1:K
        for j = i:K
            C[i,j] = cov(d, [i, j])
            C[j,i] = C[i,j]
    end
    end
    return Symmetric(C, :U)
end

"""
Test formulas using the 2 examples given in Manjunath Wilhelm (2021)
"""
function test_moments()
    println("2-d Example")
    let d = TMvNormal([0.5, 0.5], [1 1.2; 1.2 2], [-1, -Inf], [0.5, 1])
        ref_mean = [-0.152, -0.388]
        ref_cov = [0.163 0.161; 0.161 0.606]
        println("Computed mean:")
        display(round.(mean(d); digits=3))
        println("Reference mean:")
        display(round.(ref_mean; digits=3))
        println("Computed covariance:")
        display(round.(cov(d); digits=3))
        println("Reference covariance:")
        display(round.(ref_cov; digits=3))
    end

    println("3-d Example")
    let d = TMvNormal([0,0,0], [1.1 1.2 0; 1.2 2 -0.8; 0 -0.8 3], [-1, -Inf, -Inf], [0.5, Inf, Inf])
        ref_mean = [-0.210, -0.229, -0.0]
        ref_cov = [0.174 0.190 0.0; 0.190 0.898 -0.8; 0 -0.8 3.0]
        println("Computed mean:")
        display(round.(mean(d); digits=3))
        println("Reference mean:")
        display(round.(ref_mean; digits=3))
        println("Computed covariance:")
        display(round.(cov(d); digits=3))
        println("Reference covariance:")
        display(round.(ref_cov; digits=3))
    end
end

export TMvNormal, cov2cor, mean, cov, pdf, cdf, rand
end


text:
include("building.jl")

# ================
# BLDF01 - BDS01
# ================

# dense-time options
δ_dense = 0.003
algo_dense = GLGM06(δ=δ_dense)

# discrete-time options
approx_model = "nobloating"
δ_discrete = 0.01
algo_discrete = GLGM06(δ=δ_discrete, approx_model=approx_model)

# single run
sol_BLDF01_dense = solve(build_TV, alg=algo_dense, T=time_horizon)
# sol_BLDF01_discrete = solve(build_TV, alg=algo_discrete, T=time_horizon)
#
# # benchmark
# SUITE["Build"]["BLDF01-BDS01", "dense"] =
#     @benchmarkable solve($build_TV, $𝑂_BLDF01, op=GLGM06($𝑂_dense_BLDF01))
# SUITE["Build"]["BLDF01-BDS01", "discrete"] =
#     @benchmarkable solve($build_TV, $𝑂_BLDF01, op=GLGM06($𝑂_discrete_BLDF01))
#
# # ================
# # BLDC01 - BDS01
# # ================
#
# # general options
# 𝑂_BLDC01 = Options(:T => time_horizon, :mode => "check", :property => pBLDC01)
#
# # algorithm-specific options
# 𝑂_dense_BLDC01 = Options(:vars => [25], :partition => [1:24, [25], 26:48, [49]],
#                          :δ => 0.005, :block_options_init => LazySets.LinearMap)
# 𝑂_discrete_BLDC01 =
#     merge(𝑂_dense_BLDC01, Options(:discretization => "nobloating", :δ => 0.005))
#
# # single run to verify that specification holds
# sol_BLDC01_dense = solve(build_CONST, 𝑂_BLDC01, op=GLGM06(𝑂_dense_BLDC01))
# @assert sol_BLDC01_dense.satisfied
# sol_BLDC01_discrete = solve(build_CONST, 𝑂_BLDC01, op=GLGM06(𝑂_discrete_BLDC01))
# @assert sol_BLDC01_discrete.satisfied
#
# # register in benchmark suite
# SUITE["Build"]["BLDC01-BDS01", "dense"] =
#     @benchmarkable solve($build_CONST, $𝑂_BLDC01, op=GLGM06($𝑂_dense_BLDC01))
# SUITE["Build"]["BLDC01-BDS01", "discrete"] =
#     @benchmarkable solve($build_CONST, $𝑂_BLDC01, op=GLGM06($𝑂_discrete_BLDC01))

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
module BrokenRecord

export playback

using Base.Threads: nthreads, threadid

using BSON: BSON
using HTTP: HTTP, Header, Layer, Response, URI, body_was_streamed, header, insert_default!,
    mkheaders, nobody, remove_default!, request, request_uri, stack, top_layer
using JLD2: JLD2
using JLSO: JLSO
using JSON: JSON
using YAML: YAML

const FORMAT = v"1"
const DEFAULTS = Dict(
    :path => "",
    :ignore_headers => [],
    :ignore_query => [],
    :extension => "yml",
)
const STATE = map(1:nthreads()) do i
    (; responses=Response[], ignore_headers=String[], ignore_query=String[])
end

drop_keys(keys) = p -> !(p.first in keys)

get_state() = STATE[threadid()]

function reset_state()
    state = get_state()
    empty!(state.responses)
    empty!(state.ignore_headers)
    empty!(state.ignore_query)
end

"""
    configure!(;
        path=nothing,
        extension=nothing,
        ignore_headers=nothing,
        ignore_query=nothing,
    )

Set options globally so that you needn't pass keywords to every [`playback`](@ref) call.

## Keywords

- `path`: Path to the directory that contains data files.
  Any path you pass to [`playback`](@ref) will be relative to this path.
- `extension`: File extension for data files, which determines the storage backend used.
  The default is `"yml"`, which produces YAML files.
- `ignore_headers`: Names of headers to remove from requests.
- `ignore_query`: Names of query string parameters to remove from requests.
"""
function configure!(; path=nothing, extension=nothing, ignore_headers=nothing, ignore_query=nothing)
    path === nothing || (DEFAULTS[:path] = path)
    extension === nothing || (DEFAULTS[:extension] = extension)
    ignore_headers === nothing || (DEFAULTS[:ignore_headers] = ignore_headers)
    ignore_query === nothing || (DEFAULTS[:ignore_query] = ignore_query)
    return nothing
end

"""
    playback(f, path; ignore_headers=[], ignore_query=[])

Run `f` while either recording to or playing back from the file at `path`.
See [`configure!`](@ref) for more information on the available keywords.
"""
function playback(
    f, path;
    ignore_headers=DEFAULTS[:ignore_headers], ignore_query=DEFAULTS[:ignore_query],
)
    reset_state()
    state = get_state()
    append!(state.ignore_headers, ignore_headers)
    append!(state.ignore_query, ignore_query)

    path = joinpath(DEFAULTS[:path], replace(path, isspace => "_"))
    storage, path = get_storage(path, DEFAULTS[:extension])

    before_layer, custom_layer = if isfile(path)
        top_layer(stack()), PlaybackLayer
    else
        Union{}, RecordingLayer
    end

    before(custom_layer, storage, path)
    insert_default!(before_layer, custom_layer)
    return try
        f()
    finally
        remove_default!(before_layer, custom_layer)
        after(custom_layer, storage, path)
    end
end

include("recording.jl")
include("playback.jl")
include("storage.jl")

end


text:
module FHMMs

# Factorial Hidden Markov Models (FHMM)

# Reference:
# [Ghahramani1997] Ghahramani, Zoubin, and Michael I. Jordan.
# "Factorial hidden Markov models."
# Machine learning 29.2-3 (1997): 245-273.

export FHMM, fit!, fit_sv!

# Factorial HMM (FHMM)
type FHMM
    K::Int # number of states
    M::Int # number of factors
    # contributions to the means of Gaussian, shape (D, K, M)
    W::Array{Float64, 3}
    # initial state probability, shape (K, M)
    π::Array{Float64, 2}
    # transition maticies, shape (K, K, M)
    P::Array{Float64, 3}
    # covariance matrix, shape (D, D)
    C::Array{Float64, 2}

    function FHMM(D::Int, K::Int, M::Int)
        new(K, M,
            rand(D, K, M),
            ones(K, M) ./ (K*M),
            ones(K, K, M) ./ K,
            rand(D, D))
    end
end

for fname in ["completely_factorized_vfhmm",
              "structured_vfhmm"]
    include(string(fname, ".jl"))
end

end # module


text:
### related to signal analysis/synthesis, geometric transforms.

function rad2freq(r,fs)
    fn = fs/2
    #r/π = x/fn, x is in Hz.
    return fn*r/π
end

function freq2rad(x,fs)
    fn = fs/2

    return x*π/fn
end

# evaluates H(e^{j*ω}).
function computeDTFTviaformula(h,ω)
    N = length(h)
    return sum( h[n+1]*exp(-im*ω*n) for n = 0:N-1 )
end
# # test code
# ω = 99
# h = randn(5)
# Hn = h[1] + h[2]*exp(-im*ω) + h[3]*exp(-im*2*ω) + h[4]*exp(-im*3*ω) + h[5]*exp(-im*4*ω)
# Ht = computeDTFTviaformula(h,ω)
#
# println("discrepancy ", abs(Hn-Ht))

# length = 2N+1
function getdiracdeltaseq(N)
    out = zeros(Float64,2*N+1)
    out[N+1] = 1.0

    return out
end

# get highpass filter from a lowpass one.
function spectralinversion(h::Vector{T}) where T
    @assert isodd(length(h))

    out = collect( -h[i] for i = 1:length(h))
    out[div(length(h),2)+1] += one(T)

    return out
end

# assume x starts at time 0.
function upsample(x::Vector{T}, M::Int)::Vector{T} where T

    out = Vector{T}(undef, length(x)*M)
    #out = Vector{T}(undef, (length(x)-1)*M +1)
    fill!(out,zero(T))

    for n = 0:length(x)-1
        out[n*M+1] = x[n+1]
    end

    return out
end

function downsample(x::Vector{T}, M::Int)::Tuple{Vector{T},Vector{T}} where T

    x_padded = copy(x)
    Q = rem(length(x),M)
    if Q != 0
        #push!(out,zero(T))
        push!(x_padded,zeros(T,M-Q)...)
    end

    N = div(length(x_padded),M)

    out = Vector{T}(undef, N)

    for n = 0:N-1
        out[n+1] = x_padded[n*M+1]
    end

    return out, x_padded
end

function linearconv(  x::Vector{T},
                h::Vector{T}) where T

    out = Vector{T}(undef, length(x)+length(h)-1)
    for i = 1:length(x)+length(h)-1

        out[i] = sum( h[k]*getfinitesignalatindex(x,i-k+1) for k = 1:length(h) )
        # i-k+1 instead of i-k because Julia is 1-indexing.
    end

    return out
end

# remove the padding due to linear convolution. There will still be a bit of
#   shift probably due to phase-shift of system. This reasoning is untested.
# g is the filter.
# y is the signal.
function forcesamesizeafterlinearconv(g::Vector{T},y::Vector{T})::Vector{T} where T

    offset1 = div(length(g)-1,2)
    offset2 = length(g)-1-offset1

    return y[offset1:end-offset2]
end

# in progress. optimized for speed. eventually package pynogram without allocating N_chs worth of output.
# ignore entries of output that require out-of-bound h.
function fastlinearconv(a::Vector{T}, b::Vector{T})::Vector{T} where T
    m = length(a)
    n = length(b)

    out = zeros(T, m+n-1)
    @inbounds for j = 1:m
        for k = 1:n
            out[j+k-1] += a[j]*b[k]
        end
    end
    return out
end

function circconv(  x::Vector{T},
                h::Vector{T}) where T

    out = Vector{T}(undef, length(x)+length(h))
    for i = 1:length(x)+length(h)

        out[i] = sum( h[k]*x[mod(i-k,length(x))+1] for k = 1:length(h) )
    end

    return out
end

max_stars_repo_name: RoyCCWang/RWJuliaUtilities

max_issues_repo_name: RoyCCWang/RWJuliaUtilities

max_forks_repo_name: RoyCCWang/RWJuliaUtilities


text:
using LinearAlgebra
using NDTensors
using TBLIS

let
  d = 25

  nthreads = 4

  A = randomTensor(d, d, d, d)
  B = randomTensor(d, d, d, d)
  C_blas = randomTensor(d, d, d, d)

  labelsA = (1, -1, 2, -2)
  labelsB = (-2, 4, 3, -1)
  labelsC = (1, 2, 3, 4)

  println("Contracting (d x d x d x d) tensors A * B -> C, d = ", d)
  @show labelsA
  @show labelsB
  @show labelsC

  #
  # Use BLAS
  #

  NDTensors.disable_tblis()
  BLAS.set_num_threads(nthreads)

  println()
  println("Using BLAS with $nthreads threads")

  time_blas = @belapsed NDTensors.contract!($C_blas, $labelsC, $A, $labelsA, $B, $labelsB) samples =
    100

  println()
  println("Time (BLAS) = ", time_blas, " seconds")

  #
  # Use TBLIS
  #

  NDTensors.enable_tblis()
  TBLIS.set_num_threads(nthreads)

  println()
  println("Using TBLIS with $(TBLIS.get_num_threads()) threads")

  C_tblis = randomTensor(d, d, d, d)

  time_tblis = @belapsed NDTensors.contract!($C_tblis, $labelsC, $A, $labelsA, $B, $labelsB) samples =
    100

  println()
  println("Time (TBLIS) = ", time_tblis, " seconds")

  println()
  @show C_blas ≈ C_tblis

  println()
  println("Time (TBLIS) / Time (BLAS) = ", time_tblis / time_blas)
end


text:
using MappedArrays
using Test

#ambs = detect_ambiguities(ImageCore,AxisArrays,Base,Core)
#ambs = setdiff(detect_ambiguities(ImageAxes,ImageCore,AxisArrays,Base,Core), ambs)
#=
if !isempty(ambs)
    println("Ambiguities:")
    for a in ambs
        println(a)
    end
end
@test isempty(ambs)
=#

#@traitfn has_time_axis(::AA) where {AA<:AxisArray;  HasTimeAxis{AA}} = true
#@traitfn has_time_axis(::AA) where {AA<:AxisArray; !HasTimeAxis{AA}} = false

@testset "no units, no time" begin
    A = AxisArray(reshape(1:12, 3, 4), x = (1:3), y = (1:4))
    @test @inferred(timeaxis(A)) === nothing
    @test !has_time_axis(A)
    @test timedim(A) == 0
    @test @inferred(pixelspacing(A)) === (1,1)
    @test @inferred(spacedirections(A)) === ((1,0),(0,1))
    @test @inferred(coords_spatial(A)) === (1,2)
    @test spatialorder(A) === (:x, :y)  # TODO: make this inferrable
    @test @inferred(size_spatial(A)) === (3,4)
    @test @inferred(indices_spatial(A)) === (Base.OneTo(3), Base.OneTo(4))
    assert_timedim_last(A)
    @test map(istimeaxis, AxisArrays.axes(A)) == (false,false)

    @test @inferred(timeaxis(rand(3,5))) == nothing

    # deprecate
    @test nimages(A) == 1
end

@testset "units, no time" begin
    mm = u"mm"     # in real use these should be global consts
    m = u"m"
    A = AxisArray(reshape(1:12, 3, 4), x = (1mm:1mm:3mm), y = (1m:2m:7m))
    @test @inferred(timeaxis(A)) === nothing
    @test !has_time_axis(A)
    @test timedim(A) == 0
    @test nimages(A) == 1
    @test @inferred(pixelspacing(A)) === (1mm,2m)
    @test spacedirections(A) === ((1mm,0m),(0mm,2m))   # TODO: make this inferrable
    @test @inferred(coords_spatial(A)) === (1,2)
    @test spatialorder(A) === (:x,:y)
    @test @inferred(size_spatial(A)) === (3,4)
    @test @inferred(indices_spatial(A)) === (Base.OneTo(3),Base.OneTo(4))
    assert_timedim_last(A)
    @test map(istimeaxis, AxisArrays.axes(A)) == (false,false)
end

@testset "units, time" begin
    s = u"s" # again, global const
    axt = Axis(1s:1s:4s)
    A = AxisArray(reshape(1:12, 3, 4), Axis(1:3), axt)
    @test @inferred(timeaxis(A)) === axt
    @test has_time_axis(A)
    @test timedim(A) == 2
    @test nimages(A) == 4
    @test @inferred(pixelspacing(A)) === (1,)
    @test @inferred(spacedirections(A)) === ((1,),)
    @test @inferred(coords_spatial(A)) === (1,)
    @test spatialorder(A) === (:x,)
    @test @inferred(size_spatial(A)) === (3,)
    @test @inferred(indices_spatial(A)) === (Base.OneTo(3),)
    assert_timedim_last(A)
    @test map(istimeaxis, AxisArrays.axes(A)) == (false,true)
end

@testset "units, time first" begin
    s = u"s" # global const
    axt = Axis{:time}(1s:1s:4s)
    A = AxisArray(reshape(1:12, 4, 3), axt, Axis{:x}(1:3))
    @test @inferred(timeaxis(A)) === axt
    @test has_time_axis(A)
    @test timedim(A) == 1
    @test nimages(A) == 4
    @test @inferred(pixelspacing(A)) === (1,)
    @test @inferred(spacedirections(A)) === ((1,),)
    @test @inferred(coords_spatial(A)) === (2,)
    @test spatialorder(A) === (:x,)
    @test @inferred(size_spatial(A)) === (3,)
    @test @inferred(indices_spatial(A)) === (Base.OneTo(3),)
    @test_throws ErrorException assert_timedim_last(A)
    @test map(istimeaxis, AxisArrays.axes(A)) == (true,false)
end

# Possibly-ambiguous functions
@testset "ambig" begin
    A = AxisArray(rand(RGB{N0f8},3,5), :x, :y)
    @test isa(convert(Array{RGB{N0f8},2}, A), Array{RGB{N0f8},2})
    @test isa(convert(Array{Gray{N0f8},2}, A), Array{Gray{N0f8},2})
end

@testset "internal" begin
    A = AxisArray(rand(RGB{N0f8},3,5), :x, :y)
    @test ImageAxes.axtype(A) == Tuple{Axis{:x,Base.OneTo{Int}}, Axis{:y,Base.OneTo{Int}}}
end

# For testing streaming with a non-AxisArray parent
module TestStreaming
using AxisArrays, ImageAxes

struct AVIStream
    dims::NTuple{3,Int}
end
Base.ndims(::AVIStream) = 3
Base.size(A::AVIStream) = A.dims
AxisArrays.axisnames(::Type{AS}) where {AS<:AVIStream} = (:y, :x, :time)
AxisArrays.axes(A::AVIStream) = (Axis{:y}(Base.OneTo(A.dims[1])),
                                 Axis{:x}(Base.OneTo(A.dims[2])),
                                 Axis{:time}(Base.OneTo(A.dims[3])))
ImageAxes.StreamIndexStyle(::Type{AVIStream}, ::Type{typeof(read!)}) =
    IndexIncremental()

end

@testset "streaming" begin
    P = AxisArray([0 0 0 0;
                   1 2 3 4;
                   0 0 0 0], :x, :time)
    f!(dest, a) = (dest[1] = dest[3] = -0.2*a[2]; dest[2] = 0.6*a[2]; dest)
    # Next inference was special-cased for v0.6
    S = @inferred(StreamingContainer{Float64}(f!, P, Axis{:time}()))
    @test @inferred(axes(S)) === (Base.OneTo(3), Base.OneTo(4))
    @test @inferred(size(S)) == (3,4)
    @test @inferred(axes(S, 2)) === Base.OneTo(4)
    @test @inferred(size(S, 1)) === 3
    @test @inferred(length(S)) == 12
    @test @inferred(axisnames(S)) == (:x, :time)
    @test @inferred(axisvalues(S)) === (Base.OneTo(3), Base.OneTo(4))
    @test axisdim(S, Axis{:x}) == axisdim(S, Axis{:x}(1:2)) == axisdim(S, Axis{:x,UnitRange{Int}}) == 1
    @test axisdim(S, Axis{:time}) == 2
    @test_throws ErrorException axisdim(S, Axis{:y})
    @test axisdim(S, Axis{2}) == 2
    @test_throws ErrorException axisdim(S, Axis{3})
    @test @inferred(timeaxis(S)) === Axis{:time}(Base.OneTo(4))
    @test nimages(S) == 4
    @test @inferred(coords_spatial(S)) == (1,)
    @test @inferred(indices_spatial(S)) == (Base.OneTo(3),)
    @test @inferred(size_spatial(S)) == (3,)
    @test @inferred(spatialorder(S)) == (:x,)
    assert_timedim_last(S)
    for i = 1:4
        @test @inferred(S[:,i]) == [-0.2,0.6,-0.2]*i
        @test @inferred(S[2,i]) === 0.6*i
        @test @inferred(S[Axis{:time}(i)]) == [-0.2,0.6,-0.2]*i
        @test @inferred(S[Axis{:time}(i),Axis{:x}(2)]) === 0.6*i
        @test @inferred(S[Axis{:x}(2),Axis{:time}(i)]) === 0.6*i
    end
    buf = zeros(3)
    @test @inferred(getindex!(buf, S, :, 2)) == [-0.2,0.6,-0.2]*2
    @test StreamIndexStyle(S) === IndexAny()
    @test StreamIndexStyle(zeros(2,2)) === IndexAny()
    # Non-AbstractArray parent
    @test_throws DimensionMismatch StreamingContainer{UInt8}(read!, TestStreaming.AVIStream((1080,1920,10000)), Axis{:foo}())
    V = StreamingContainer{UInt8}(read!, TestStreaming.AVIStream((1080,1920,10000)), Axis{:time}())
    @test size(V) == (1080,1920,10000)
    @test axisnames(V) == (:y, :x, :time)
    @test StreamIndexStyle(V) === IndexIncremental()
    # internal
    @test ImageAxes.streamingaxisnames(S) == (:time,)
    @test ImageAxes.filter_streamed((1,2), S) == (2,)
end

nothing


text:
using Documenter, uCSV

makedocs(
    root = Pkg.dir("uCSV", "docs"),
    source = "src",
    build = "build",
    doctest = true,
    modules = [uCSV],
    format   = :html,
    sitename = "uCSV.jl",
    checkdocs = :all,
    debug = true,
    pages = Any["Home" => "index.md",
                "Manual" => Any["man/defaults.md",
                                "man/headers.md",
                                "man/dataframes.md",
                                "man/delimiters.md",
                                "man/missingdata.md",
                                "man/declaring-column-element-types.md",
                                "man/declaring-column-vector-types.md",
                                "man/international.md",
                                "man/customparsers.md",
                                "man/quotes-escapes.md",
                                "man/comments-skiplines.md",
                                "man/malformed.md",
                                "man/url.md",
                                "man/compressed.md",
                                "man/unsupported.md",
                                "man/write.md",
                                "man/benchmarks.md"]])

deploydocs(
    repo = "github.com/cjprybol/uCSV.jl.git",
    julia = "1.0",
    target = "build",
    deps = nothing,
    make = nothing)


text:
# ------------------------------------------------------------------------------------------
# # Linear algebra in Julia
# Based on work by Andreas Noack Jensen
#
# ## Basic linalg ops
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# First let's define a random matrix
# ------------------------------------------------------------------------------------------

A = rand(1:4,3,3)

# ------------------------------------------------------------------------------------------
# Define a vector of ones
# ------------------------------------------------------------------------------------------

x = fill(1.0, (3))

# ------------------------------------------------------------------------------------------
# Notice that $A$ has type Array{Int64,2} but $x$ has type Array{Float64,1}. Julia defines
# the aliases Vector{Type}=Array{Type,1} and Matrix{Type}=Array{Type,2}.
#
# Many of the basic operations are the same as in other languages
# #### Multiplication
# ------------------------------------------------------------------------------------------

b = A*x

# ------------------------------------------------------------------------------------------
# #### Transposition
# As in other languages `A'` is the conjugate transpose
# ------------------------------------------------------------------------------------------

A'

# ------------------------------------------------------------------------------------------
# and we can get the transpose with
# ------------------------------------------------------------------------------------------

transpose(A)

# ------------------------------------------------------------------------------------------
# #### Transposed multiplication
# Julia allows us to write this without *
# ------------------------------------------------------------------------------------------

A'A

# ------------------------------------------------------------------------------------------
# #### Solving linear systems
# The problem $Ax=b$ for ***square*** $A$ is solved by the \ function.
# ------------------------------------------------------------------------------------------

A\b

# ------------------------------------------------------------------------------------------
# ## Special Matrix Structures
#
# Matrix structure is very important in linear algebra. To see *how* important it is, let's
# work with a larger linear system. Use the LinearAlgebra standard package to get access to
# structured matrices:
# ------------------------------------------------------------------------------------------

using LinearAlgebra

n = 1000
A = randn(n,n);

# ------------------------------------------------------------------------------------------
# Julia can often infer special matrix structure
# ------------------------------------------------------------------------------------------

Asym = A + A'
issymmetric(Asym)

# ------------------------------------------------------------------------------------------
# but sometimes floating point error might get in the way.
# ------------------------------------------------------------------------------------------

Asym_noisy = copy(Asym)
Asym_noisy[1,2] += 5eps()

issymmetric(Asym_noisy)

# ------------------------------------------------------------------------------------------
# Luckily we can declare structure explicitly with, for example, `Diagonal`, `Triangular`,
# `Symmetric`, `Hermitian`, `Tridiagonal` and `SymTridiagonal`.
# ------------------------------------------------------------------------------------------

Asym_explicit = Symmetric(Asym_noisy);

# ------------------------------------------------------------------------------------------
# Let's compare how long it takes Julia to compute the eigenvalues of `Asym`, `Asym_noisy`,
# and `Asym_explicit`
# ------------------------------------------------------------------------------------------

@time eigvals(Asym);

@time eigvals(Asym_noisy);

@time eigvals(Asym_explicit);

# ------------------------------------------------------------------------------------------
# In this example, using `Symmetric()` on `Asym_noisy` made our calculations about `5x` more
# efficient :)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### A big problem
# Using the `Tridiagonal` and `SymTridiagonal` types to store tridiagonal matrices makes it
# possible to work with potentially very large tridiagonal problems. The following problem
# would not be possible to solve on a laptop if the matrix had to be stored as a (dense)
# `Matrix` type.
# ------------------------------------------------------------------------------------------

n = 1_000_000;
A = SymTridiagonal(randn(n), randn(n-1));
@time eigmax(A)

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
module BasicNT

using Primes

export primesieve, primes2,
         nextprime2, prevprime2,
         primefactors, twinprimes,
         coprime, linmod,
         ordermod, primroot,
         eulerphi

  include("primes.jl")

end # module