text:
# distance(k::Stationary, X::AbstractMatrix, Y::AbstractMatrix) = pairwise(metric(k), X, Y)
distance(k::Stationary, x::AbstractVector, y::AbstractVector) = evaluate(metric(k), x, y)

distance(k::Stationary, X::AbstractMatrix, data::EmptyData) = distance(k, X)
distance(k::Isotropic, X::AbstractMatrix, data::IsotropicData) = data.R
function distance(k::Stationary, X::AbstractMatrix)
    nobs = size(X, 2)
    return distance!(Matrix{eltype(X)}(undef, nobs, nobs), metric(k), X)
end
distance!(dist::AbstractMatrix, k::Stationary, X::AbstractMatrix) = distance!(dist, metric(k), X)
function distance!(dist::AbstractMatrix, m::PreMetric, X::AbstractMatrix)
    dim, nobs = size(X)
    @inbounds @simd for i in 1:nobs
        dist[i,i] = 0.0
        for j in 1:i-1
            dist[i,j] = dist[j,i] = distij(m, X, i, j, dim)
        end
    end
    return dist
end
function distance(k::Stationary, X::AbstractMatrix, Y::AbstractMatrix)
    nobsx = size(X, 2)
    nobsy = size(Y, 2)
    return distance!(Matrix{eltype(Y)}(undef, nobsx, nobsy), metric(k), X, Y)
end
distance!(dist::AbstractMatrix, k::Stationary, X::AbstractMatrix, Y::AbstractMatrix) = distance!(dist, metric(k), X, Y)
function distance!(dist::AbstractMatrix, m::PreMetric, X::AbstractMatrix, Y::AbstractMatrix)
    dimx, nobsx = size(X)
    dimy, nobsy = size(Y)
    dimx == dimy || error("size(X, 1) != size(Y, 1)")
    @inbounds @simd for i in 1:nobsx
        for j in 1:nobsy
            dist[i,j] = distij(m, X, Y, i, j, dimx)
        end
    end
    return dist
end
function distance(k::StationaryARD, X::AbstractMatrix, data::StationaryARDData)
    distance(k, X)
end
@inline _SqEuclidean_ijk(X::AbstractMatrix,i::Int,j::Int,k::Int)=(X[k,i]-X[k,j])^2
@inline _SqEuclidean_ijk(X1::AbstractMatrix,X2::AbstractMatrix,i::Int,j::Int,k::Int)=(X1[k,i]-X2[k,j])^2
@inline function _SqEuclidean_ij(X::AbstractMatrix,i::Int,j::Int,dim::Int)
    s = 0.0
    @inbounds @simd for k in 1:dim
        s += _SqEuclidean_ijk(X,i,j,k)
    end
    return s
end
@inline function _SqEuclidean_ij(X1::AbstractMatrix,X2::AbstractMatrix,i::Int,j::Int,dim::Int)
    s = 0.0
    @inbounds @simd for k in 1:dim
        s += _SqEuclidean_ijk(X1,X2,i,j,k)
    end
    return s
end
# Squared Euclidean
# @inline distijk(dist::SqEuclidean,X::AbstractMatrix,i::Int,j::Int,k::Int)=_SqEuclidean_ijk(X,i,j,k)
# @inline distijk(dist::SqEuclidean,X1::AbstractMatrix,X2::AbstractMatrix,i::Int,j::Int,k::Int)=_SqEuclidean_ijk(X1,X2,i,j,k)
@inline distij(dist::SqEuclidean,X::AbstractMatrix,i::Int,j::Int,dim::Int)=_SqEuclidean_ij(X,i,j,dim)
@inline distij(dist::SqEuclidean,X1::AbstractMatrix,X2::AbstractMatrix,i::Int,j::Int,dim::Int)=_SqEuclidean_ij(X1,X2,i,j,dim)

# Euclidean
@inline distij(dist::Euclidean,X::AbstractMatrix,i::Int,j::Int,dim::Int)=√_SqEuclidean_ij(X,i,j,dim)
@inline function distij(dist::Euclidean,X1::AbstractMatrix,X2::AbstractMatrix,i::Int,j::Int,dim::Int)
    if X1 === X2 && i == j
        zero(eltype(X2))
    else
        √_SqEuclidean_ij(X1,X2,i,j,dim)
    end
end

@inline _WeightedSqEuclidean_ijk(weights::AbstractVector,X::AbstractMatrix,i::Int,j::Int,k::Int)=(X[k,i]-X[k,j])^2*weights[k]
@inline _WeightedSqEuclidean_ijk(weights::AbstractVector,X1::AbstractMatrix,X2::AbstractMatrix,i::Int,j::Int,k::Int)=(X1[k,i]-X2[k,j])^2*weights[k]
@inline function _WeightedSqEuclidean_ij(weights::AbstractVector,X::AbstractMatrix,i::Int,j::Int,dim::Int)
    s = 0.0
    @inbounds @simd for k in 1:dim
        s += _WeightedSqEuclidean_ijk(weights,X,i,j,k)
    end
    return s
end
@inline function _WeightedSqEuclidean_ij(weights::AbstractVector,X1::AbstractMatrix,X2::AbstractMatrix,i::Int,j::Int,dim::Int)
    s = 0.0
    @inbounds @simd for k in 1:dim
        s += _WeightedSqEuclidean_ijk(weights,X1,X2,i,j,k)
    end
    return s
end

# Weighted Squared Euclidean
@inline distijk(dist::WeightedSqEuclidean,X::AbstractMatrix,i::Int,j::Int,k::Int)=_WeightedSqEuclidean_ijk(dist.weights,X,i,j,k)
@inline distijk(dist::WeightedSqEuclidean,X1::AbstractMatrix,X2::AbstractMatrix,i::Int,j::Int,k::Int)=_WeightedSqEuclidean_ijk(dist.weights,X1,X2,i,j,k)
@inline distij(dist::WeightedSqEuclidean,X::AbstractMatrix,i::Int,j::Int,dim::Int)=_WeightedSqEuclidean_ij(dist.weights,X,i,j,dim)
@inline distij(dist::WeightedSqEuclidean,X1::AbstractMatrix,X2::AbstractMatrix,i::Int,j::Int,dim::Int)=_WeightedSqEuclidean_ij(dist.weights,X1,X2,i,j,dim)

# Weighted Euclidean
@inline dist2ijk(dist::WeightedEuclidean,X::AbstractMatrix,i::Int,j::Int,k::Int)=_WeightedSqEuclidean_ijk(dist.weights,X,i,j,k)
@inline dist2ijk(dist::WeightedEuclidean,X1::AbstractMatrix,X2::AbstractMatrix,i::Int,j::Int,k::Int)=_WeightedSqEuclidean_ijk(dist.weights,X1,X2,i,j,k)
@inline distij(dist::WeightedEuclidean,X::AbstractMatrix,i::Int,j::Int,dim::Int)=√_WeightedSqEuclidean_ij(dist.weights,X,i,j,dim)
@inline function distij(dist::WeightedEuclidean,X1::AbstractMatrix,X2::AbstractMatrix,i::Int,j::Int,dim::Int)
    if X1 === X2 && i == j
        zero(eltype(X2))
    else
        √_WeightedSqEuclidean_ij(dist.weights,X1,X2,i,j,dim)
    end
end


text:
function test_item_mean()
    println("-- Testing ItemMean recommender")

    data = DataAccessor(sparse([1 2 3; 4 5 6]))
    recommender = ItemMean(data)
    build!(recommender)
    actual = predict(recommender, 1, 1)

    @test actual == 2.5
end

test_item_mean()

ext: jl


text:
@test :S | :R == (:S,:R)
@test :S | [:R, :D] == (:S, [:R, :D])
@test [:S,:G] | [:R, :D] == ([:S,:G], [:R, :D])

@test P(:S|:R) == CPD(:S,:R)
@test P([:S,:G] | [:R, :D]) == CPD([:S,:G],[:R,:D])


text:
using Test, Onda, Random, Dates

@testset "$(repr(name)) serializer" for (name, options) in [(:lpcm, nothing),
                                         (Symbol("lpcm.zst"), Dict(:level => 2))]
    signal = Signal([:a, :b, :c], Nanosecond(0), Nanosecond(0), :unit, 0.25, -0.5, Int16,
                    50.5, name, options)
    samples = encode(Samples(signal, false, rand(MersenneTwister(1), 3, Int(50.5 * 10)))).data
    s = serializer(signal)
    bytes = serialize_lpcm(samples, s)
    name == :lpcm && @test bytes == reinterpret(UInt8, vec(samples))
    @test deserialize_lpcm(bytes, s) == samples
    @test deserialize_lpcm(IOBuffer(bytes), s) == samples
    io = IOBuffer()
    serialize_lpcm(io, samples, s)
    seekstart(io)
    @test take!(io) == bytes
end

@test_throws ArgumentError Onda.serializer_constructor_for_file_extension(Val(:extension))
@test_throws ErrorException Onda.register_file_extension_for_serializer(:extension, LPCM)


text:
using Test

function count(arr)
    if isempty(arr)
        return 0
    end
    return 1 + count(arr[2:end])
end

@test count([1,2,3,4]) == 4


text:
using ZMQ, Test

@info("Testing with ZMQ version $(ZMQ.version)")

@testset "ZMQ contexts" begin
	ctx=Context()
	@test ctx isa Context
	@test (ctx.io_threads = 2) == 2
	@test ctx.io_threads == 2
	ZMQ.close(ctx)

	#try to create socket with expired context
	@test_throws StateError Socket(ctx, PUB)
end

@testset "ZMQ sockets" begin
	s=Socket(PUB)
	@test s isa Socket
	ZMQ.close(s)

	s1=Socket(REP)
	s1.sndhwm = 1000
	s1.linger = 1
	s1.routing_id = "abcd"

	@test s1.routing_id == "abcd"
	@test s1.sndhwm === 1000
	@test s1.linger === 1
	@test s1.rcvmore === false

	s2=Socket(REQ)
	@test s1.type == REP
	@test s2.type == REQ

	ZMQ.bind(s1, "tcp://*:5555")
	ZMQ.connect(s2, "tcp://localhost:5555")

	msg = Message("test request")
	# Test similar() and copy() fixes in https://github.com/JuliaInterop/ZMQ.jl/pull/165
	# Note that we have to send this message to work around
	# https://github.com/JuliaInterop/ZMQ.jl/issues/166
	@test similar(msg, UInt8, 12) isa Vector{UInt8}
	@test copy(msg) == codeunits("test request")
	ZMQ.send(s2, msg)
	@test unsafe_string(ZMQ.recv(s1)) == "test request"
	ZMQ.send(s1, Message("test response"))
	@test unsafe_string(ZMQ.recv(s2)) == "test response"

	ZMQ.send(s2, "test request 2")
	@test ZMQ.recv(s1, String) == "test request 2"
	ZMQ.send(s1, Vector(codeunits("test response 2")))
	@test String(ZMQ.recv(s2, Vector{UInt8})) == "test response 2"
	ZMQ.send(s2, 3.14159)
	@test ZMQ.recv(s1, Float64) === 3.14159
	ZMQ.send(s1, [314159, 12345])
	@test ZMQ.recv(s2, Vector{Int}) == [314159, 12345]

	# Test task-blocking behavior
	c = Base.Condition()
	global msg_sent = false
	@async begin
		global msg_sent
		sleep(0.5)
		msg_sent = true
		ZMQ.send(s2, Message("test request"))
		@test (unsafe_string(ZMQ.recv(s2)) == "test response")
		notify(c)
	end

	# This will hang forver if ZMQ blocks the entire process since
	# we'll never switch to the other task
	@test unsafe_string(ZMQ.recv(s1)) == "test request"
	@test msg_sent == true
	ZMQ.send(s1, Message("test response"))
	wait(c)

	# Test _Message task-blocking behavior, similar to above
	c = Base.Condition()
	msg_sent = false
	@async begin
		global msg_sent
		sleep(0.5)
		msg_sent = true
		ZMQ.send(s2, "another test request")
		@test ZMQ.recv(s2, String) == "another test response"
		notify(c)
	end
	@test ZMQ.recv(s1, String) == "another test request"
	@test msg_sent == true
	ZMQ.send(s1, "another test response")
	wait(c)

	ZMQ.send(s2, Message("another test request"))
	msg = ZMQ.recv(s1)
	o=convert(IOStream, msg)
	seek(o, 0)
	@test String(take!(o)) == "another test request"
	ZMQ.send(s1) do io
		print(io, "buffer ")
		print(io, "this")
	end
	@test String(take!(ZMQ.recv(s2, IOBuffer))) == "buffer this"

	@testset "Message AbstractVector interface" begin
		m = Message("1")
		@test m[1]==0x31
		@test (m[1]=0x32) === 0x32
		@test unsafe_string(m)=="2"
		finalize(m)
	end

	# ZMQ.close(s1); ZMQ.close(s2) # should happen when context is closed
	ZMQ.close(ZMQ._context) # immediately close global context rather than waiting for exit
end

@testset "ZMQ resource management" begin
    local leaked_req_socket, leaked_rep_socket
    ZMQ.Socket(ZMQ.REQ) do req_socket
	leaked_req_socket = req_socket

	ZMQ.Socket(ZMQ.REP) do rep_socket
	    leaked_rep_socket = rep_socket

	    ZMQ.bind(rep_socket, "inproc://tester")
	    ZMQ.connect(req_socket, "inproc://tester")

	    ZMQ.send(req_socket, "Mr. Watson, come here, I want to see you.")
	    @test unsafe_string(ZMQ.recv(rep_socket)) == "Mr. Watson, come here, I want to see you."
	    ZMQ.send(rep_socket, "Coming, Mr. Bell.")
	    @test unsafe_string(ZMQ.recv(req_socket)) == "Coming, Mr. Bell."
	end

	@test !ZMQ.isopen(leaked_rep_socket)
    end
    @test !ZMQ.isopen(leaked_req_socket)

    local leaked_ctx
    ZMQ.Context() do ctx
	leaked_ctx = ctx

	@test isopen(ctx)
    end
    @test !isopen(leaked_ctx)
end


text:
import .ForwardDiff

@generated function SLEEFPirates.tanh_fast(x::ForwardDiff.Dual{T,S,N}) where {T,S,N}
    quote
        $(Expr(:meta,:inline))
        t = tanh_fast(x.value)
        ∂t = vfnmadd_fast(t, t, one(S))
        p = x.partials
        ForwardDiff.Dual(t, ForwardDiff.Partials(Base.Cartesian.@ntuple $N n -> mul_fast(∂t, p[n])))
    end
end
@generated function SLEEFPirates.sigmoid_fast(x::ForwardDiff.Dual{T,S,N}) where {T,S,N}
    quote
        $(Expr(:meta,:inline))
        s = sigmoid_fast(x.value)
        ∂s = vfnmadd_fast(s,s,s)
        p = x.partials
        ForwardDiff.Dual(s, ForwardDiff.Partials(Base.Cartesian.@ntuple $N n -> mul_fast(∂s, p[n])))
    end
end
@generated function VectorizationBase.relu(x::ForwardDiff.Dual{T,S,N}) where {T,S,N}
    quote
        $(Expr(:meta,:inline))
        v = x.value
        z = zero(v)
        cmp = v < z
        r = ifelse(cmp, z, v)
        p = x.partials
        ForwardDiff.Dual(r, ForwardDiff.Partials(Base.Cartesian.@ntuple $N n -> ifelse(cmp, z, p[n])))
    end
end
@generated function init_dual(v::Tuple{Vararg{AbstractSIMD,A}}) where {A}
    res = Expr(:tuple)
    q = Expr(:block, Expr(:meta,:inline))
    for a ∈ 1:A
        v_a = Symbol(:v_,a)
        push!(q.args, Expr(:(=), v_a, Expr(:ref, :v, a)))
        partials = Expr(:tuple)
        for i ∈ 1:A
            push!(partials.args, Expr(:call, i == a ? :one : :zero, v_a))
        end
        push!(res.args, :(ForwardDiff.Dual($v_a, ForwardDiff.Partials($partials))))
    end
    push!(q.args, res)
    q
end
@generated function dual_store!(∂p::Tuple{Vararg{AbstractStridedPointer,A}}, p::AbstractStridedPointer, ∂v, im::Vararg{Any,N}) where {A,N}
    quote
        $(Expr(:meta,:inline))
        v = ∂v.value
        ∂ = ∂v.partials
        Base.Cartesian.@nextract $N im im
        Base.Cartesian.@ncall $N VectorizationBase.vnoaliasstore! p v im  # store
        Base.Cartesian.@nexprs $A a -> begin # for each of `A` partials
            ∂p_a = ∂p[a]
            ∂_a = ∂[a]
            Base.Cartesian.@ncall $N VectorizationBase.vnoaliasstore! ∂p_a ∂_a im # store
        end
        nothing
    end
end


text:
#!/usr/bin/env julia
using Plots, Random, Printf

function initCondGen() #get random initial conditions for mass/radius, position, and velocity
    m=rand(1:1500,3)./10 #3 random masses between 0.1 and 150 solar masses
    rad=m.^0.8 #3 radii based on masses in solar units
    m=m.*2e30 #convert to SI kg
    rad=rad.*7e8 #convert to SI m
    pos1=rand(-35:35,2) #random initial coordinates x & y for first body, AU
    function genPos2(pos1)
        accept2=false
        while accept2==false
            pos2=rand(-35:35,2) #random initial coordinates for second body, AU
            dist21=sqrt((pos1[1]-pos2[1])^2+(pos1[2]-pos2[2])^2)
            if (dist21*1.5e11)>(rad[1]+rad[2]) #they aren't touching
                accept2=true
                return pos2
            end
        end
    end
    pos2=genPos2(pos1)
    function genPos3(pos1,pos2)
        accept3=false
        while accept3==false
            pos3=rand(-35:35,2) #random initial coordinates for third body, AU
            dist31=sqrt((pos1[1]-pos3[1])^2+(pos1[2]-pos3[2])^2)
            dist32=sqrt((pos2[1]-pos3[1])^2+(pos2[2]-pos3[2])^2)
            if (dist31*1.5e11)>(rad[1]+rad[3]) && (dist32*1.5e11)>(rad[2]+rad[3]) #3rd isn't touching either
                accept3=true
                return pos3
            end
        end
    end
    pos3=genPos3(pos1,pos2)
    pos=[pos1[1],pos1[2],pos2[1],pos2[2],pos3[1],pos3[2]].*1.5e11 #convert accepted positions to SI, m
    v=rand(-7e3:7e3,6) #random x & y velocities with mag between -10 & 10 km/s, totally arbitrary...
    #r=[x1,y1,x2,y2,x3,y3,v1x,v1y,v2x,v2y,v3x,v3y]
    r=[pos[1],pos[2],pos[3],pos[4],pos[5],pos[6],v[1],v[2],v[3],v[4],v[5],v[6]]
    open("initCond.txt","w") do f #save initial conditions to file in folder where script is run
        write(f,"m1=$(@sprintf("%.1f",(m[1]/2e30))) m2=$(@sprintf("%.1f",(m[2]/2e30))) m3=$(@sprintf("%.1f",(m[3]/2e30))) (solar masses)\nv1x=$(v[1]/1e3) v1y=$(v[2]/1e3) v2x=$(v[3]/1e3) v2y=$(v[4]/1e3) v3x=$(v[5]/1e3) v3y=$(v[6]/1e3) (km/s)\nx1=$(pos1[1]) y1=$(pos1[2]) x2=$(pos2[1]) y2=$(pos2[2]) x3=$(pos3[1]) y3=$(pos3[2]) (AU from center)")
    end
    return r, rad, m
end

function dR(r,m) #function we will use RK4 on to approximate solution
    G=6.67408313131313e-11# Nm^2/kg^2
    M1,M2,M3=m[1],m[2],m[3] #kg
    x1,x2,x3=r[1],r[3],r[5] #m
    y1,y2,y3=r[2],r[4],r[6] #m

    c1,c2,c3=G*M1,G*M2,G*M3 #Nm^2/kg
    r1_2=sqrt((x1-x2)^2+(y1-y2)^2) #distance from 1->2, m
    r1_3=sqrt((x1-x3)^2+(y1-y3)^2) #distance from 1->3, m
    r2_3=sqrt((x2-x3)^2+(y2-y3)^2) #distance from 2->3, m

    v1X,v2X,v3X=r[7],r[9],r[11] #these are our change in position after dt (dr/dt*dt=dr)
    v1Y,v2Y,v3Y=r[8],r[10],r[12] #m after * dt

    #get change in velocity from accelerations (d^2r/dt^2*dt=dv/dt*dt=dv)
    dx1=-(c2*(x1-x2)/(r1_2^3))-(c3*(x1-x3)/(r1_3^3)) #d^2x/dt^2 for 1, m/s after * dt
    dx2=-(c1*(x2-x1)/(r1_2^3))-(c3*(x2-x3)/(r2_3^3)) #d^2x/dt^2 for 2, m/s
    dx3=-(c1*(x3-x1)/(r1_3^3))-(c2*(x3-x2)/(r2_3^3)) #d^2x/dt^2 for 3, m/s
    dy1=-(c2*(y1-y2)/(r1_2^3))-(c3*(y1-y3)/(r1_3^3)) #d^2y/dt^2 for 1, m/s
    dy2=-(c1*(y2-y1)/(r1_2^3))-(c3*(y2-y3)/(r2_3^3)) #d^2y/dt^2 for 2, m/s
    dy3=-(c1*(y3-y1)/(r1_3^3))-(c2*(y3-y2)/(r2_3^3)) #d^2y/dt^2 for 3, m/s

    return [v1X,v1Y,v2X,v2Y,v3X,v3Y,dx1,dy1,dx2,dy2,dx3,dy3]
end

function gen3Body(stopCond=[10,100],numSteps=10000) #default stop conditions of 10 yrs and 100 AU sep
    tStop=stopCond[1]*365*24*3600 #convert to SI s
    sepStop=stopCond[2]*1.5e11 #convert to SI m
    stop=false
    currentT=0
    t=range(0,stop=tStop,length=(numSteps+1)) #+1 because I don't want 0 to count
    stepSize=tStop/numSteps
    x1=zeros(length(t))
    y1=zeros(length(t))
    x2=zeros(length(t))
    y2=zeros(length(t))
    x3=zeros(length(t))
    y3=zeros(length(t))
    r,rad,m=initCondGen()
    min12=rad[1]+rad[2]
    min13=rad[1]+rad[3]
    min23=rad[2]+rad[3]
    i=1
    stopT=maximum(t)

    #implement RK4 to model solutions to differential equations
    while stop==false
        if currentT==stopT || currentT>stopT #in case of rounding error or something
            stop=true
        elseif i>numSteps #inf loop failsafe
            stop=true
            println("error: shouldn't have gotten here")
        else
            x1[i]=r[1]
            y1[i]=r[2]
            x2[i]=r[3]
            y2[i]=r[4]
            x3[i]=r[5]
            y3[i]=r[6]

            k1=stepSize*dR(r,m)
            k2=stepSize*dR(r.+0.5.*k1,m)
            k3=stepSize*dR(r.+0.5.*k2,m)
            k4=stepSize*dR(r.+k3,m)
            r+=(k1.+2.0*k2.+2.0.*k3.+k4)./6

            #check separation after each dt step
            sep12=sqrt((x1[i]-x2[i])^2+(y1[i]-y2[i])^2)
            sep13=sqrt((x1[i]-x3[i])^2+(y1[i]-y3[i])^2)
            sep23=sqrt((x3[i]-x2[i])^2+(y3[i]-y2[i])^2)

            if sep12<min12 || sep13<min13 || sep23<min23 || sep12>sepStop || sep13>sepStop || sep23>sepStop
                stop=true #stop if collision happens or body is ejected
                t=range(0,stop=currentT,length=i) #t should match pos vectors
                x1=x1[1:i] #don't want trailing zeros
                y1=y1[1:i]
                x2=x2[1:i]
                y2=y2[1:i]
                x3=x3[1:i]
                y3=y3[1:i]
            end
            i+=1
            currentT+=stepSize #next step
        end
    end
    return [x1,y1,x2,y2,x3,y3], t, m, rad
end

function getInteresting3Body(minTime=0) #in years, defaults to 0
    #sometimes random conditions result in a really short animation where things
    #just crash into each other/fly away, so this function throws away those
    yearSec=365*24*3600
    interesting=false
    i=1
    while interesting==false
        plotData,t,m,rad=gen3Body([50,150],50000)
        if (maximum(t)/yearSec)>minTime #only return if simulation runs for longer than minTime
            println(maximum(t)/yearSec) #tell me how many years we are simulating
            open("cron_log.txt","a") do f #for cron logging, a flag = append
                write(f,"$(maximum(t)/yearSec)\n")
            end
            return plotData,t,m,rad
            interesting=true
        elseif i>50 #computationally expensive so don't want to go forever
            interesting=true #render it anyways I guess because sometimes it's fun?
            println("did not find interesting solution in number of tries allotted, running anyways")
            println(maximum(t)/yearSec) #how many years simulation runs for
            open("cron_log.txt","a") do f #for cron logging
                write(f,"$(maximum(t)/yearSec)\n")
            end
            return plotData,t,m,rad
        end
        i+=1
    end
end

function getLims(pos,padding) #determines plot limits at each frame, padding in units of pos
    x=[pos[1],pos[3],pos[5]]
    xMin=minimum(x)
    xMax=maximum(x)
    dx=xMax-xMin
    y=[pos[2],pos[4],pos[6]]
    yMin=minimum(y)
    yMax=maximum(y)
    dy=yMax-yMin
    if dx>dy
        #use x for square
        xlims=[xMin-padding,xMax+padding]
        ylims=[yMin-padding,yMin+dx+padding]
    else
        #use y for square
        xlims=[xMin-padding,xMin+dy+padding]
        ylims=[yMin-padding,yMax+padding]
    end
    return xlims,ylims
end

function getColors(m,c) #places colors of objects according to mass/size
    #c=[:biggest,:medium,:smallest] (order of input colors)
    maxM=maximum(m)
    minM=minimum(m)
    colors=[:blue,:blue,:blue] #testing
    if m[1]==maxM
        colors[1]=c[1]
        if m[2]==minM
            colors[2]=c[3]
            colors[3]=c[2]
        else
            colors[3]=c[3]
            colors[2]=c[2]
        end
    elseif m[2]==maxM
        colors[2]=c[1]
        if m[1]==minM
            colors[1]=c[3]
            colors[3]=c[2]
        else
            colors[3]=c[3]
            colors[1]=c[2]
        end
    else
        colors[3]=c[1]
        if m[1]==minM
            colors[1]=c[3]
            colors[2]=c[2]
        else
            colors[2]=c[3]
            colors[1]=c[2]
        end
    end
    return colors
end

function makeCircleVals(r,center=[0,0])
    xOffset=center[1]
    yOffset=center[2]
    xVals=[r*cos(i)+xOffset for i=0:(pi/64):(2*pi)]
    yVals=[r*sin(i)+yOffset for i=0:(pi/64):(2*pi)]
    return xVals,yVals
end

plotData,t,m,rad=getInteresting3Body(15)
c=[:DodgerBlue,:Gold,:Tomato] #most massive to least massive, also roughly corresponds to temp
colors=getColors(m,c)
#adding fake stars
numStars=2500
starsX=zeros(numStars)
starsY=zeros(numStars)
for i=1:numStars
    num=rand(-200:200,2) #box size is 70 AU but we need some extra padding for movement
    starsX[i]=num[1]
    starsY[i]=num[2]
end

#this new way runs significantly faster (~2x improvement over @anim)
#Downside is it spams folder with png images of every frame and must manually compile with ffmpeg
#Comment out and use older way (after this below) if performance/specific formatting is not an issue

plotLoadPath="/home/kirk/Documents/3Body/tmpPlots/"
threeBodyAnim=Animation(plotLoadPath,String[])
for i=1:33:length(t) #this makes animation scale ~1 sec/year with other conditions
    GR.inline("png") #added to eneable cron/jobber compatibility, also this makes frames generate WAY faster? Prior to adding this when run from cron/jobber frames would stop generating at 408 for some reason.
    gr(legendfontcolor = plot_color(:white)) #legendfontcolor=:white plot arg broken right now (at least in this backend)
    print("$(@sprintf("%.2f",i/length(t)*100)) % complete\r") #output percent tracker
    pos=[plotData[1][i],plotData[2][i],plotData[3][i],plotData[4][i],plotData[5][i],plotData[6][i]] #current pos
    limx,limy=getLims(pos./1.5e11,10) #convert to AU, 10 AU padding
    p=plot(plotData[1][1:i]./1.5e11,plotData[2][1:i]./1.5e11,label="",linecolor=colors[1],linealpha=max.((1:i) .+ 10000 .- i,2500)/10000) #plot orbits up to i
    p=plot!(plotData[3][1:i]./1.5e11,plotData[4][1:i]./1.5e11,label="",linecolor=colors[2],linealpha=max.((1:i) .+ 10000 .- i,2500)/10000) #linealpha argument causes lines to decay
    p=plot!(plotData[5][1:i]./1.5e11,plotData[6][1:i]./1.5e11,label="",linecolor=colors[3],linealpha=max.((1:i) .+ 10000 .- i,2500)/10000) #example: alpha=max.((1:i) .+ 100 .- i,0) causes only last 100 to be visible
    p=scatter!(starsX,starsY,markercolor=:white,markersize=:1,label="") #fake background stars
    star1=makeCircleVals(rad[1],[plotData[1][i],plotData[2][i]]) #generate circles with appropriate sizes for each star
    star2=makeCircleVals(rad[2],[plotData[3][i],plotData[4][i]]) #at current positions
    star3=makeCircleVals(rad[3],[plotData[5][i],plotData[6][i]])
    p=plot!(star1[1]./1.5e11,star1[2]./1.5e11,label="$(@sprintf("%.1f", m[1]./2e30))",color=colors[1],fill=true) #plot star circles with labels
    p=plot!(star2[1]./1.5e11,star2[2]./1.5e11,label="$(@sprintf("%.1f", m[2]./2e30))",color=colors[2],fill=true)
    p=plot!(star3[1]./1.5e11,star3[2]./1.5e11,label="$(@sprintf("%.1f", m[3]./2e30))",color=colors[3],fill=true)
    p=plot!(background_color=:black,background_color_legend=:transparent,foreground_color_legend=:transparent,
        background_color_outside=:white,aspect_ratio=:equal,legendtitlefontcolor=:white) #formatting for plot frame
    p=plot!(xlabel="x: AU",ylabel="y: AU",title="Random Three Body Problem\nt: $(@sprintf("%0.2f",t[i]/365/24/3600)) yrs after start",
        legend=:best,xaxis=("x: AU",(limx[1],limx[2]),font(9,"Courier")),yaxis=("y: AU",(limy[1],limy[2]),font(9,"Courier")),
        grid=false,titlefont=font(14,"Courier"),size=(720,721),legendfontsize=8,legendtitle="Mass (in solar masses)",legendtitlefontsize=8) #add in axes/title/legend with formatting
    frame(threeBodyAnim,p) #generate the frame
    closeall() #close plots
end

#threeBodyFile="3Body_fps30.mp4"

#crf is compression value (17 or 18 "visually lossless"), pix_fmt is for twitter specific vid req, -b:v specifies target bitrate, -vcodec specifies codec (h264 in this case) -y says overwrite existing file
#run( `ffmpeg -framerate 30 -i $plotLoadPath"%06d.png" -vcodec libx264 -pix_fmt yuv420p -profile:v high -b:v 2048K -y -vf "scale=720:720,setdar=1/1" $threeBodyFile` ) #-vf scale=720:72 -crf 25
#run( `ffmpeg -framerate 30 -i $plotLoadPath"%06d.png" -c:v libx264 -preset slow -coder 1 -movflags +faststart -g 15 -crf 18 -pix_fmt yuv420p -profile:v high -y -bf 2 -vf "scale=720:720,setdar=1/1" $threeBodyFile` ) #all this bullshit to hopefully satisfy twitter requirements

#NOTE: moved ffmpeg commands to shell script

#old (simpler) way of generating animation
#uncomment and use this way if you just want a simple animation saved and don't
#care about performance/specific formatting of video.

# threeBodyAnim=@animate for i=1:length(t)
#     gr(legendfontcolor = plot_color(:white)) #plot arg broken right now in Julia
#     print("$(@sprintf("%.2f",i/length(t)*100)) % complete\r") #output percent tracker
#     pos=[plotData[1][i],plotData[2][i],plotData[3][i],plotData[4][i],plotData[5][i],plotData[6][i]] #current pos
#     limx,limy=getLims(pos./1.5e11,5) #convert to AU, 5 AU padding
#     plot(plotData[1][1:i]./1.5e11,plotData[2][1:i]./1.5e11,label="",linecolor=colors[1])
#     plot!(plotData[3][1:i]./1.5e11,plotData[4][1:i]./1.5e11,label="",linecolor=colors[2])
#     plot!(plotData[5][1:i]./1.5e11,plotData[6][1:i]./1.5e11,label="",linecolor=colors[3])
#     scatter!(starsX,starsY,markercolor=:white,markersize=:1,label="") #fake background stars
#     star1=makeCircleVals(rad[1],[plotData[1][i],plotData[2][i]])
#     star2=makeCircleVals(rad[2],[plotData[3][i],plotData[4][i]])
#     star3=makeCircleVals(rad[3],[plotData[5][i],plotData[6][i]])
#     plot!(star1[1]./1.5e11,star1[2]./1.5e11,label="$(@sprintf("%.1f", m[1]./2e30))",color=colors[1],fill=true)
#     plot!(star2[1]./1.5e11,star2[2]./1.5e11,label="$(@sprintf("%.1f", m[2]./2e30))",color=colors[2],fill=true)
#     plot!(star3[1]./1.5e11,star3[2]./1.5e11,label="$(@sprintf("%.1f", m[3]./2e30))",color=colors[3],fill=true)
#     plot!(background_color=:black,background_color_legend=:transparent,background_color_outside=:white,aspect_ratio=:equal,legendtitlefontcolor=:white) #legendfontcolor=:white
#     plot!(xlabel="x: AU",ylabel="y: AU",title="Random Three Body Problem\nt: $(@sprintf("%0.2f",t[i]/365/24/3600)) yrs after start",
#         legend=:best,xaxis=("x: AU",(limx[1],limx[2]),font(12,"Courier")),yaxis=("y: AU",(limy[1],limy[2]),font(12,"Courier")),
#         grid=false,titlefont=font(24,"Courier"),size=(720,720),legendfontsize=12,legendtitle="Mass (in solar masses)",legendtitlefontsize=14)
#     end every 25

#mp4(threeBodyAnim,"3Body_fps30.mp4",fps=30)
#OR
#gif(threeBodyAnim,"3Body_fps30.gif",fps=30)


text:
using EyeFix
using Base.Test

dummy = [
1 1 2
2 2 3
3 3 4
4 9 10
5 10 11
6 11 12
7 20 20
]

res1 = IVT(dummy; velthresh = 3)
res2 = IVT(dummy[:, 1], dummy[:, 2:3]; velthresh = 3)
res3 = IVT(dummy[:, 1], dummy[:, 2], dummy[:, 3]; velthresh = 3)
shouldbe = [1.0, 1.0, 1.0, 2.0, 2.0, 2.0, 3.0]

@assert res1[:, 4] == shouldbe
@assert res2[:, 4] == shouldbe
@assert res3[:, 4] == shouldbe

res1 = IDT(dummy)
res2 = IDT(dummy[:, 1], dummy[:, 2:3])
res3 = IDT(dummy[:, 1], dummy[:, 2], dummy[:, 3])
shouldbe = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0]

@assert res1[:, 4] == shouldbe
@assert res2[:, 4] == shouldbe
@assert res3[:, 4] == shouldbe
res = IDT(dummy; durthresh = 10000, disthresh = 100000)

dd = dummy[:, 2:3]
res1 = kmeans(dd; k = 2)
res2 = kmeans(dd[:, 1], dd[:, 2]; k = 2)

res = GMM(dummy[:, 1]; k = 2)


text:
# CWT (continuous wavelet transform)
# cwt(Y::AbstractVector, ::ContWave)

@doc """
     wave = cwt(Y::AbstractArray{T,N}, c::CWT{W, S, WT}, daughters, rfftPlan =
             plan_rfft([1]), fftPlan = plan_fft([1])) where {N, T<:Real,
                                                             S<:Real,
                                                             U<:Number,
                                                             W<:WaveletBoundary,
                                                             WT<:ContWaveClass}

  return the continuous wavelet transform along the first axis with averaging.
  `wave`, is (signalLength)×(nscales+1)×(previous dimensions), of type `T` of
  `Y`. `averagingLength` defines the number of octaves (powers of 2) that are
  replaced by an averaging function. This has form `averagingType`, which can be
  one of `Father()` or `Dirac()`- in the `Father()` case, it uses the same form
  as for the wavelets, while the `Dirac` uses a constant window. If you have
  sampling information, you will need to scale wave by δt^(1/p). The default
  assumption is that the sampling rate is 2kHz.

"""
function cwt(Y::AbstractArray{T,N}, cWav::CWT, daughters, fftPlans = 1) where {N,T}
    @assert typeof(N) <: Integer
    # vectors behave a bit strangely, so we reshape them
    if N == 1
        Y = reshape(Y, (length(Y), 1))
    end
    n1 = size(Y, 1)

    _, nScales, _ = getNWavelets(n1, cWav)
    #construct time series to analyze, pad if necessary
    x = reflect(Y, boundaryType(cWav)()) #this function is defined below

    # check if the plans we were given are dummies or not
    x̂, fromPlan = prepSignalAndPlans(x, cWav, fftPlans)
    # If the vector isn't long enough to actually have any other scales, just
    # return the averaging
    if nScales <= 0 || size(daughters, 2) == 1
        daughters = daughters[:, 1:1]
        nScales = 1
    end

    if isAnalytic(cWav.waveType)
        OutType = ensureComplex(T)
    else
        OutType = T
    end

    wave = zeros(OutType, size(x)..., nScales)  # result array
    # faster if we put the example index on the outside loop through all scales
    # and compute transform
    if isAnalytic(cWav.waveType)
        if eltype(x) <: Real
            analyticTransformReal!(wave, daughters, x̂, fromPlan, cWav.averagingType)
        else
            analyticTransformComplex!(wave, daughters, x̂, fromPlan, cWav.averagingType)
        end
    else
        otherwiseTransform!(wave, daughters, x̂, fromPlan, cWav.averagingType)
    end
    wave = permutedims(wave, [1, ndims(wave), (2:(ndims(wave)-1))...])
    ax = axes(wave)
    wave = wave[1:n1, ax[2:end]...]
    if N == 1
        wave = dropdims(wave, dims = 3)
    end

    return wave
end

function ensureComplex(T)
    if T <: Real
        return Complex{T}
    else
        return T
    end
end

# there are 4 cases to deal with
#       input Type | Real | Complex
#       analytic?  |----------------
#             yes  | both | fft
#              no  | rfft | fft
#  Analytic on Real input
function prepSignalAndPlans(x::AbstractArray{T}, cWav::CWT{W,S,WaTy,N,true}, fftPlans) where {T<:Real,W,S,WaTy,N}
    # analytic wavelets that are being applied on real inputs
    if fftPlans isa Tuple{<:AbstractFFTs.Plan{<:Real},<:AbstractFFTs.Plan{<:Complex}}
        # they handed us the right kind of thing, so no need to make new ones
        x̂ = fftPlans[1] * x
        fromPlan = fftPlans[2]
    else
        toPlan = plan_rfft(x, 1)
        x̂ = toPlan * x
        fromPlan = plan_fft(x, 1)
    end
    return x̂, fromPlan
end

#  Non-analytic on Real input
function prepSignalAndPlans(x::AbstractArray{T}, cWav::CWT{W,S,WaTy,N,false}, fftPlans) where {T<:Real,W,S,WaTy,N}
    # real wavelets that are being applied on real inputs
    if fftPlans isa AbstractFFTs.Plan{<:Real}
        # they handed us the right kind of thing, so no need to make new ones
        x̂ = fftPlans * x
        fromPlan = fftPlans
    else
        fromPlan = plan_rfft(x, 1)
        x̂ = fromPlan * x
    end
    return x̂, fromPlan
end
# complex input
function prepSignalAndPlans(x::AbstractArray{T}, cWav, fftPlans) where {T<:Complex}
    # real wavelets that are being applied on real inputs
    if fftPlans isa AbstractFFTs.Plan{<:Complex}
        # they handed us the right kind of thing, so no need to make new ones
        x̂ = fftPlans * x
        fromPlan = fftPlans
    else
        fromPlan = plan_fft(x, 1)
        x̂ = fromPlan * x
    end
    return x̂, fromPlan
end

# analytic on real data with an averaging function
function analyticTransformReal!(wave, daughters, x̂, fftPlan, ::Union{Father,Dirac})
    outer = axes(x̂)[2:end]
    n1 = size(x̂, 1)
    isSourceEven = mod(size(wave, 1) + 1, 2)
    # the averaging function isn't analytic, so we need to do both positive and
    # negative frequencies
    @views tmpWave = x̂ .* daughters[:, 1]
    @views wave[(n1+1):end, outer..., 1] = reverse(conj.(tmpWave[2:end-isSourceEven, outer...]), dims = 1)
    @views wave[1:n1, outer..., 1] = tmpWave
    @views wave[:, outer..., 1] = fftPlan \ (wave[:, outer..., 1])  # averaging
    for j = 2:size(daughters, 2)
        @views wave[1:n1, outer..., j] = x̂ .* daughters[:, j]
        wave[:, outer..., j] = fftPlan \ (wave[:, outer..., j])  # wavelet transform
    end
end

# analytic on complex data with an averaging function
function analyticTransformComplex!(wave, daughters, x̂, fftPlan, ::Union{Father,Dirac})
    outer = axes(x̂)[2:end]
    n1 = size(daughters, 1)
    isSourceEven = mod(size(wave, 1) + 1, 2)
    # the averaging function isn't analytic, so we need to do both positive and
    # negative frequencies
    @views positiveFreqs = x̂[1:n1, outer...] .* daughters[:, 1]
    @views negativeFreqs = x̂[(n1-isSourceEven+1):end, outer...] .* reverse(conj.(daughters[2:end, 1]))
    @views wave[(n1-isSourceEven+1):end, outer..., 1] = negativeFreqs
    @views wave[1:n1, outer..., 1] = positiveFreqs
    @views wave[:, outer..., 1] = fftPlan \ (wave[:, outer..., 1])  # averaging
    for j = 2:size(daughters, 2)
        @views wave[1:n1, outer..., j] = x̂[1:n1, outer...] .* daughters[:, j]
        @views wave[:, outer..., j] = fftPlan \ (wave[:, outer..., j])  # wavelet transform
    end
end

function analyticTransformComplex!(wave, daughters, x̂, fftPlan, averagingType)
    outer = axes(x̂)[2:end]
    n1 = size(x̂, 1)
    for j = 1:size(daughters, 2)
        @views wave[1:n1, outer..., j] = x̂[1:n1, outer...] .* daughters[:, j]
        @views wave[:, outer..., j] = fftPlan \ (wave[:, outer..., j])  # wavelet transform
    end
end

# analytic on real data without an averaging function
function analyticTransformReal!(wave, daughters, x̂, fftPlan, ::NoAve)
    outer = axes(x̂)[2:end]
    n1 = size(x̂, 1)
    # the no averaging version
    for j = 1:size(daughters, 2)
        wave[1:n1, outer..., j] = x̂ .* daughters[:, j]
        wave[:, outer..., j] = fftPlan \ (wave[:, outer..., j])  # wavelet transform
    end
end

function otherwiseTransform!(wave::AbstractArray{<:Real}, daughters, x̂, fromPlan, averagingType)
    # real wavelets on real data: that just makes sense
    outer = axes(x̂)[2:end]
    n1 = size(x̂, 1)
    for j = 1:size(daughters, 2)
        @views tmp = x̂ .* daughters[:, j]
        @views wave[:, outer..., j] = fromPlan \ tmp  # wavelet transform
    end
end

# if it isn't analytic, the output is complex only if the input is complex
function otherwiseTransform!(wave::AbstractArray{<:Complex}, daughters, x̂, fromPlan, averagingType)
    # applying a real transform to complex data is maybe a bit odd, but you do you
    outer = axes(x̂)[2:end]
    n1 = size(daughters, 1)
    isSourceEven = mod(size(fromPlan, 1) + 1, 2)
    for j = 1:size(daughters, 2)
        @views wave[1:n1, outer..., j] = @views x̂[1:n1, outer...] .* daughters[:, j]
        @views wave[n1-isSourceEven+1:end, outer..., j] = x̂[n1-isSourceEven+1:end, outer...] .* reverse(conj.(daughters[2:end, j]))
        @views wave[:, outer..., j] = fromPlan \ (wave[:, outer..., j])  # wavelet transform
    end
end

function reflect(Y, bt)
    n1 = size(Y, 1)
    if typeof(bt) <: ZPBoundary
        base2 = ceil(Int, log2(n1))   # power of 2 nearest to N
        x = cat(Y, zeros(2^(base2) - n1, size(Y)[2:end]...), dims = 1)
    elseif typeof(bt) <: SymBoundary
        x = cat(Y, reverse(Y, dims = 1), dims = 1)
    else
        x = Y
    end
    return x
end

function cwt(Y::AbstractArray{T}, c::CWT{W}; varArgs...) where {T<:Number,W<:WaveletBoundary}
    daughters, ω = computeWavelets(size(Y, 1), c; varArgs...)
    return cwt(Y, c, daughters)
end

cwt(Y::AbstractArray{T}, w::ContWaveClass; varargs...) where {T<:Number} = cwt(Y, CWT(w); varargs...)
cwt(Y::AbstractArray{T}) where {T<:Real} = cwt(Y, Morlet())

abstract type InverseType end
struct DualFrames <: InverseType end
struct NaiveDelta <: InverseType end
struct PenroseDelta <: InverseType end

"""
    icwt(res::AbstractArray, cWav::CWT, inverseStyle::InverseType=PenroseDelta())
Compute the inverse wavelet transform using one of three dual frames. The default uses delta functions with weights chosen via a least squares method, the `PenroseDelta()` below. This is chosen as a default because the Morlet wavelets tend to fail catastrophically using the canonical dual frame (the `dualFrames()` type).

    icwt(res::AbstractArray, cWav::CWT, inverseStyle::PenroseDelta)
Return the inverse continuous wavelet transform, computed using the simple dual frame ``β_jδ_{ji}``, where ``β_j`` is chosen to solve the least squares problem ``\\|Ŵβ-1\\|_2^2``, where ``Ŵ`` is the Fourier domain representation of the `cWav` wavelets. In both this case and `NaiveDelta()`, the fourier transform of ``δ`` is the constant function, thus this least squares problem.

    icwt(res::AbstractArray, cWav::CWT, inverseStyle::NaiveDelta)
Return the inverse continuous wavelet transform, computed using the simple dual frame ``β_jδ_{ji}``, where ``β_j`` is chosen to negate the scale factor ``(^1/_s)^{^1/_p}``. Generally less accurate than choosing the weights using `PenroseDelta`. This is the method discussed in Torrence and Compo.

    icwt(res::AbstractArray, cWav::CWT, inverseStyle::dualFrames)
Return the inverse continuous wavelet transform, computed using the canonical dual frame ``\\tilde{\\widehat{ψ}} = \\frac{ψ̂_n(ω)}{∑_n\\|ψ̂_n(ω)\\|^2}``. The algorithm is to compute the cwt again, but using the canonical dual frame; consiquentially, it is the most computationally intensive of the three algorithms, and typically the best behaved. Will be numerically unstable if the high frequencies of all of the wavelets are too small however, and tends to fail spectacularly in this case.

"""
function icwt(res::AbstractArray, cWav::CWT, ::PenroseDelta)
    Ŵ = computeWavelets(size(res, 1), cWav)[1]
    β = computeDualWeights(Ŵ, cWav)
    testDualCoverage(β, Ŵ)
    compXRecon = sum(res .* β, dims = 2)
    imagXRecon = irfft(im * rfft(imag.(compXRecon), 1), size(compXRecon, 1)) # turns out the dual frame for the imaginary part is rather gross in the time domain
    return imagXRecon + real.(compXRecon)
end

function icwt(res::AbstractArray, cWav::CWT, ::NaiveDelta)
    Ŵ = computeWavelets(size(res, 1), cWav)[1]
    β = computeNaiveDualWeights(Ŵ, cWav, size(res, 1))
    testDualCoverage(β, Ŵ)
    compXRecon = sum(res .* β, dims = 2)
    imagXRecon = irfft(im * rfft(imag.(compXRecon), 1), size(compXRecon, 1)) # turns out the dual frame for the imaginary part is rather gross in the time domain
    return imagXRecon + real.(compXRecon)
end

function icwt(res::AbstractArray, cWav::CWT, ::DualFrames)
    Ŵ = computeWavelets(size(res, 1), cWav)[1]
    canonDualFrames = explicitConstruction(Ŵ)
    testDualCoverage(canonDualFrames)
    convolved = cwt(res, cWav, canonDualFrames)
    ax = axes(convolved)
    @views xRecon = sum(convolved[:, i, i, ax[4:end]...] for i = 1:size(Ŵ, 2))
    return xRecon
end

icwt(Y::AbstractArray, w::ContWaveClass; varargs...) = icwt(Y, CWT(w), PenroseDelta(); varargs...)
icwt(Y::AbstractArray; varargs...) = icwt(Y, Morlet(), PenroseDelta(); varargs...)

# CWT (continuous wavelet transform directly) TODO: direct if sufficiently small


text:
"""
    function Scale(ratio)

A simple wrapper that will multiply inputs by the ratio.

```jldoctest
julia> using AudioSchedules

julia> Scale(3)(2)
6
```
"""
struct Scale{Ratio}
    ratio::Ratio
end

(scale::Scale)(thing) = thing * scale.ratio

export Scale

const DIGITS = of(:maybe, raw("-")), of(:some, short(:digit))
const QUOTIENT = pattern(
    CONSTANTS.start,
    of(:maybe, capture(DIGITS..., name = "numerator")),
    of(:maybe, raw("/"), capture(DIGITS..., name = "denominator")),
    of(:maybe, "o", capture(DIGITS..., name = "octave")),
    CONSTANTS.stop,
)

get_parse(something, _) = parse(Int, something)
get_parse(::Nothing, default) = default

function q_str(interval_string)
    a_match = match(QUOTIENT, interval_string)
    if a_match === nothing
        throw(ParseError("Can't parse interval $interval_string"))
    end
    get_parse(a_match["numerator"], 1) // get_parse(a_match["denominator"], 1) *
    (2 // 1)^get_parse(a_match["octave"], 0)
end

"""
    q"interval"

Create a musical interval. You can specify a numerator (which defaults to 1), denominator (which defaults to 1), and an octave shift (which defaults to 0).

```jldoctest
julia> using AudioSchedules

julia> q"1"
1//1

julia> q"3/2"
3//2

julia> q"2/3o1"
4//3

julia> q"2/3o-1"
1//3

julia> q"o2"
4//1

julia> q"1 + 1"
ERROR: LoadError: Base.Meta.ParseError("Can't parse interval 1 + 1")
[...]
```
"""
macro q_str(interval_string::AbstractString)
    esc(q_str(interval_string))
end
export @q_str

# TODO: use FFT here?
struct SawTooth{overtones} end

"""
    SawTooth(overtones)

Build a saw-tooth wave from its partials, starting with the fundamental (1), up to
`overtones`.

To increase richness but also buziness, increase `overtones`.

```jldoctest
julia> using AudioSchedules

julia> SawTooth(3)(π / 4)
0.9185207636218614
```
"""
SawTooth(overtones) = SawTooth{overtones}()

export SawTooth

const ADJUST = 2 / pi

function (saw::SawTooth{overtones})(an_angle) where {overtones}
    ADJUST * sum(ntuple(let an_angle = an_angle
        function (overtone)
            sin_fast(overtone * an_angle) / overtone
        end
    end, Val{overtones}()))
end

max_stars_repo_name: bramtayl/AudioSchedulers

max_forks_repo_name: bramtayl/AudioSchedulers


text:
"
Author: Daniel O'Connor
Email: dt.oconnor@outlook.com

File re-creating previous SVD work and new work.
"
include("samplers/tools.jl")
include("samplers/SA.jl")
include("samplers/SVMC.jl")
include("samplers/SVMC_TF.jl")
include("samplers/pimc.jl")
include("dynamic_models/spin_vector_O2_dynamics.jl")
include("dynamic_models/spin_vector_O3_dynamics.jl")

using Random
using ..SA
using ..SVMC
using ..SVMCTF
using ..PIMC
using ..SamplingTools
using ..SpinVectorO2Dynamics
using ..SpinVectorO3Dynamics
using LinearAlgebra
using DifferentialEquations
using Profile
using BenchmarkTools
using NPZ
using Interpolations
using Plots
using LsqFit

function closed_O2()
    """
    simulating code from here
    https://arxiv.org/pdf/1305.4904.pdf
    """
    h = Dict()
    j = Dict()
    for i =1:8
        if i <= 4
            h[i] = -1.0
            j[(i, i+4)] = -1.0
            j[(i, (i%4) + 1)] = -1.0
        else
            h[i] = 1.0
        end
    end

    s = collect(0:0.001:1)
    max_field = 1.0
    a_sch= max_field .* (1.0 .- s)
    b_sch= max_field .* s
    a_fn = SamplingTools.monotonic_interpolation(s, a_sch)
    b_fn = SamplingTools.monotonic_interpolation(s, b_sch)

    # schedule_data = npzread("data\\DWAVE_schedule.npz")
    # a_fn = SamplingTools.monotonic_interpolation(schedule_data[:, 1], 0.5.*schedule_data[:, 2])
    # b_fn = SamplingTools.monotonic_interpolation(schedule_data[:, 1], 0.5.*schedule_data[:, 3])
    # a_sch = a_fn.(s)
    # b_sch = b_fn.(s)

    run_time = 1000.0
    anneal_time = 1000.0    # for T = 1000 in paper
    result, sol = SpinVectorO2Dynamics.solve(h, j, anneal_time=anneal_time,
                                         run_time=run_time, a_fn=a_fn,
                                         b_fn=b_fn, saveat=s)

    run_time = 1000.0
    anneal_time = 200.0     # for T = 200 in paper
    result2, sol2 = SpinVectorO2Dynamics.solve(h, j, anneal_time=anneal_time,
                                      run_time=run_time, a_fn=a_fn,
                                      b_fn=b_fn, saveat=s)

    final_state, energy = result
    final_state2, energy2 = result2
    # ts is the times at which are saved on the ODE solver
    ts = run_time .* s
    thetas_in = zeros(Float64, length(ts))
    thetas_out = zeros(Float64, length(ts))
    thetas_in2 = zeros(Float64, length(ts))
    thetas_out2 = zeros(Float64, length(ts))
    for i = eachindex(ts)
        thetas_in[i] = mean(sol.u[i].x[1][1:4]) / pi
        thetas_out[i] = mean(sol.u[i].x[1][5:8]) / pi
        thetas_in2[i] = mean(sol2.u[i].x[1][1:4]) / pi
        thetas_out2[i] = mean(sol2.u[i].x[1][5:8]) / pi
    end

    print("T = 200 ns, solution = $(final_state2), energy = $(energy2)\n",
          "T = 1000 ns, solution = $(final_state), energy = $(energy) ")

    plot(ts, thetas_in, xaxis="t / ns", yaxis="θ / π", legend=false, ylims=[-0.1, 1.1], linecolor=:red)
    plot!(ts, thetas_out, linecolor=:red, legend=false)
    plot!(ts, thetas_in2, linecolor=:black, legend=false)
    plot!(ts, thetas_out2, linecolor=:black, legend=false)
end

function simulate_O3()
    """
    simulating code from here (plots are inverted due to the different
    signs in the Hamiltonian)
    https://arxiv.org/pdf/1403.4228.pdf
    """
    h = Dict()
    j = Dict()
    nspins = 8
    for i =1:8
        if i <= 4
            h[i] = -1.0
            j[(i, i+4)] = -1.0
            j[(i, (i%4) + 1)] = -1.0
        else
            h[i] = 1.0
        end
    end

    temp = 2.226 / (2*pi)

    s = collect(0:0.0001:1)
    num_steps = length(s)
    # approximate the schedule
    #  0.1099 * b(s=1) = temp = 17 mK
    #  0.2834 * b(s≈0.541) = temp
    a0 = 34.25 / (2*pi)
    b0 = 1 / (2*pi)
    bt = temp / 0.1099
    at = 1e-9   # assumption
    # form of b = quadratic
    b_fn = (x) -> (2.28 * x^2) + (0.78 * x) +0.01 # + (1/(2*pi))
    sa = [0.0, 0.38, 0.541, 0.8, 0.9, 1]
    a = [a0, 5/(2*pi), temp, 3e-4, 2e-6, 1e-9]
    a_fn = SamplingTools.monotonic_interpolation(sa, a)
    # a_fn = (x) -> 0.024074802432616997 * exp(5.427221233043444*(1.0 - x).^2) -
    #         0.024074802432616997

    # schedule_data = npzread("data\\DWAVE_schedule.npz")
    # a_fn = SamplingTools.monotonic_interpolation(schedule_data[:, 1], 0.5.*schedule_data[:, 2])
    # b_fn = SamplingTools.monotonic_interpolation(schedule_data[:, 1], 0.5.*schedule_data[:, 3])
    # a_sch = a_fn.(s)
    # b_sch = b_fn.(s)

    # # DQA schedule
    # a_sch, b_sch = SamplingTools.DQA_schedule(1, 2*N, a_sch, b_sch,
    #                                           sx=0.2, cx=0.0, c1=0.0)
    # # convert to functions
    # a_fn = Dict([(i,  SamplingTools.monotonic_interpolation(s, a_sch[i])) for i=1:nspins])
    # b_fn = Dict([(i,  SamplingTools.monotonic_interpolation(s, b_sch[i])) for i=1:nspins])

    # max_field = a0
    # a_sch= max_field .* (1.0 .- s)
    # b_sch= max_field .* s
    # a_fn = SamplingTools.monotonic_interpolation(s, a_sch)
    # b_fn = SamplingTools.monotonic_interpolation(s, b_sch)
    #
    run_time = 20000.0
    anneal_time = 20000.0
    # open system enables multi-threading for parallel computation
    num_samples = 10
    open_sys = false
    friction_constant = 1e-6
    bloch = true
    if open_sys
        atol, reltol = 1e-2, 1e-2
    else
        atol, reltol = 1e-9, 1e-6
    end

    # Long anneal times require very low tolerances
    results, sol = SpinVectorO3Dynamics.solve(h, j, anneal_time=anneal_time,
                                             run_time=run_time, a_fn=a_fn,
                                             b_fn=b_fn, temp=temp,
                                             friction_constant=friction_constant,
                                             trajectories=num_samples,
                                             saveat=s,
                                             open_system=open_sys, Bloch=bloch,
                                             abstol=atol, reltol=reltol,
                                             initial_theta_state=fill(0.5 * pi, 8),
                                             initial_phi_state=zeros(Float64, 8))

    # ts is the times at which are saved on the ODE solver
    ts = LinRange(0.0, run_time, num_steps)
    # plot(sol)
    if open_sys
        final_states, energies, occurences = results
        print("Solutions = $(final_states) \n energies = $(energies) \n occurences = $(occurences)")
        med = DifferentialEquations.EnsembleAnalysis.timeseries_point_median(sol, ts)
        # low = DifferentialEquations.EnsembleAnalysis.timeseries_point_quantile(sol, 0.025, ts)
        # high = DifferentialEquations.EnsembleAnalysis.timeseries_point_quantile(sol, 0.975, ts)
        print(size(med), size(s))
        magnetization = med[16:24, :]'
    else
        final_state, energy = results
        magnetization = zeros(Float64, num_steps, length(h))
        if !bloch
            print("Solution = $(final_state) \n energy = $(energy) \n final θs = $(sol.u[end].x[1])")
            for i = 1:num_steps
                magnetization[i, :] = cos.(sol.u[i].x[1])
            end
        else
            print("Solution = $(final_state) \n energy = $(energy) \n final ⟨σ⟩ = $(sol.u[end][(2*nspins + 1):(3*nspins)])")
            for i = 1:num_steps
                magnetization[i, :] = sol.u[i][(2*nspins + 1):(3*nspins)]
            end
        end
    end
    plot(s, magnetization, xaxis="\$ t  / t_f \$", yaxis="\$ \\langle \\sigma ^z \\rangle  \$", legend=false)
end

function main_dynamics()
    d = 0.1
    N = 2
    h, j = SamplingTools.PFC_dict(d, N=N)
    num_samples = 1000
    temp = 0.01226 * 1.380649e-23 / (1e9 * 6.62607015e-34)

    s = collect(0:0.001:1)
    num_steps = length(s)
    max_field = 3.0
    a_sch= max_field .* (1.0 .- s)
    b_sch= max_field .* s
    a_fn = SamplingTools.monotonic_interpolation(s, a_sch)
    b_fn = SamplingTools.monotonic_interpolation(s, b_sch)

    schedule_data = npzread("data\\DWAVE_schedule.npz")
    a_fn = SamplingTools.monotonic_interpolation(schedule_data[:, 1], 0.5.*schedule_data[:, 2])
    b_fn = SamplingTools.monotonic_interpolation(schedule_data[:, 1], 0.5.*schedule_data[:, 3])
    a_sch = a_fn.(s)
    b_sch = b_fn.(s)

    # # DQA schedule
    # a_sch, b_sch = SamplingTools.DQA_schedule(1, 2*N, a_sch, b_sch,
    #                                           sx=0.2, cx=0.0, c1=0.0)
    # # convert to functions
    # a_fn = Dict([(i,  SamplingTools.monotonic_interpolation(s, a_sch[i])) for i=1:nspins])
    # b_fn = Dict([(i,  SamplingTools.monotonic_interpolation(s, b_sch[i])) for i=1:nspins])

    run_time = 20000.0
    anneal_time = 20000.0
    # open system enables multi-threading for parallel computation
    open_sys = false
    friction_constant = 2e-2

    # anneal_time = 1000.0

    results, sol = SpinVectorO2Dynamics.solve(h, j, anneal_time=anneal_time,
                                             run_time=run_time, a_fn=a_fn,
                                             b_fn=b_fn, temp=temp, saveat=s,
                                             friction_constant=friction_constant,
                                             trajectories=num_samples,
                                             open_system=open_sys,
                                             abstol=1e-7, reltol=1e-6)

    # ts is the times at which are saved on the ODE solver
    ts = LinRange(0.0, run_time, num_steps)
    if open_sys
        final_states, energies, occurences = results
        print("Solutions = $(final_states) \n energies = $(energies) \n occurences = $(occurences)")
        med = DifferentialEquations.EnsembleAnalysis.timeseries_point_median(sol, ts)
        # low = DifferentialEquations.EnsembleAnalysis.timeseries_point_quantile(sol, 0.025, ts)
        # high = DifferentialEquations.EnsembleAnalysis.timeseries_point_quantile(sol, 0.975, ts)
        magnetization = cos.(med[1:(2*N), :]')

    else
        final_state, energy = results
        print("Solution = $(final_state) \n energy = $(energy) \n final θs = $(sol.u[end].x[1])")
        magnetization = zeros(Float64, num_steps, length(h))
        for i = 1:num_steps
            magnetization[i, :] = cos.(sol.u[i].x[1])
        end
    end
    plot(ts, magnetization, xaxis="\$ t \\, / \\,  ns \$", yaxis="\$ \\langle \\sigma ^z \\rangle  \$", legend=false, ylims=[-1, 1])
end

max_stars_repo_name: dtoconnor/Annealing-Samplers

max_issues_repo_name: dtoconnor/Annealing-Samplers

max_forks_repo_name: dtoconnor/Annealing-Samplers


text:
export same_sign,
       rectify

"""
    same_sign(A::AbstractArray{N}; [optimistic]::Bool=false) where {N}

Check whether all elements of the given array have the same sign.

### Input

- `A`          -- array
- `optimistic` -- (optional; default: `false`) flag for expressing that the
                  expected result is `true`

### Output

`true` if and only if all elements in `M` have the same sign.

### Algorithm

If `optimistic` is `false`, we check the sign of the first element and compare
to the sign of all elements.

If `optimistic` is `true`, we compare the absolute element sum with the sum of
the absolute of the elements; this is faster if the result is `true` because
there is no branching.

```math
    |\\sum_i A_i| = \\sum_i |A_i|
```
"""
function same_sign(A::AbstractArray{N}; optimistic::Bool=false) where {N}
    if optimistic
        return sum(abs, A) == abs(sum(A))
    else
        if isempty(A)
            return true
        end
        @inbounds if first(A) >= zero(N)
            return all(e -> e >= zero(N), A)
        else
            return all(e -> e <= zero(N), A)
        end
    end
end

"""
    _rationalize(::Type{T}, x::AbstractVecOrMat{N}, tol::Real) where {T<:Integer, N<:AbstractFloat}

Approximate an array of floating point numbers as a rational vector with entries
of the given integer type.

### Input

- `T`   -- (optional, default: `Int`) integer type to represent the rationals
- `x`   -- vector of floating point numbers
- `tol` -- (optional, default: `eps(N)`) tolerance the result at entry `i`
           will differ from `x[i]` by no more than `tol`

### Output

An array of type `Rational{T}` where the `i`-th entry is the rationalization
of the `i`-th component of `x`.

### Notes

See also [`Base.rationalize`](@ref).
"""
function _rationalize(::Type{T}, x::AbstractVecOrMat{N}, tol::Real) where {T<:Integer, N<:AbstractFloat}
    return rationalize.(Ref(T), x, Ref(tol))
end

# method extensions
_rationalize(::Type{T}, x::AbstractVecOrMat{N}; tol::Real=eps(N)) where {T<:Integer, N<:AbstractFloat} =  _rationalize(T, x, tol)
_rationalize(x::AbstractVecOrMat{N}; kwargs...) where {N<:AbstractFloat} = _rationalize(Int, x; kwargs...)

# fallback for scalars
_rationalize(::Type{T}, x::AbstractFloat, tol::Real) where {T<:Integer} = rationalize(T, x, tol)

# nested vectors
function _rationalize(::Type{T}, x::AbstractVecOrMat{<:AbstractVecOrMat{N}}, tol::Real) where {T<:Integer, N<:AbstractFloat}
    return _rationalize.(Ref(T), x, Ref(tol))
end

"""
    rectify(x::AbstractArray{N}) where {N<:Real}

Rectify an array, i.e., take the element-wise maximum with zero.

### Input

- `x` -- array

### Output

A copy of the array where each negative entry is replaced by zero.
"""
function rectify(x::AbstractArray{N}) where {N<:Real}
    return map(xi -> max(xi, zero(N)), x)
end


text:
"""
    detectPWPeaks(signal::Vector{<:Real}, fs::Real; windowLength::Int = 300, tuning::Real = 0.7)

Find Peaks in a pulsewave signal with an algorithm proposed by Nenova, B., & Iliev, I. (2010). An automated algorithm for fast pulse wave detection. International Journal Bioautomation, 14(3), 203.

# Args:

* signal: Pulswave data
* fs: Sampling Frequency [Hz]
* windowLenght: Optional length of the window in ms. (standard: 300ms)
* tuning: fine tuning to filter Peaks below a certain height

# Return:

Returns a Vector containing the Peak indices. Use signal[indices] to access the Peak values.
"""
function detectPWPeaks(signal::Vector{<:Real}, fs::Real; windowLength::Int = 300, tuning::Real = 0.7)
    # bandpass filter
    lowPass = Lowpass(12, fs = fs)
    highPass = Highpass(0.5, fs = fs)
    butterWorth = Butterworth(1)
    filterLowPass = digitalfilter(lowPass, butterWorth)
    filterHighPass = digitalfilter(highPass, butterWorth)
    lowFiltered = filt(filterLowPass, signal)
    filtered = filt(filterHighPass, lowFiltered)

    # getting right window length

    window = Int(ceil(windowLength/1000 * fs))
    winLen = 1:window:length(signal)
    vals = Vector{Float64}(undef, length(winLen) - 1)
    indices = Vector{Int}(undef, length(winLen) - 1)
    for ii = 1:length(winLen)-1
        (val, ind) = findmax(filtered[winLen[ii]:winLen[ii+1]-1])
        vals[ii] = val
        indices[ii] = ind + (ii - 1) * window
    end
    #remove all below mean
    #avg = mean(filtered)

    avg2 = tuning * mean(vals)
    deleteat!(indices, findall(x -> x < avg2, vals))
    deleteat!(vals, findall(x -> x < avg2, vals))

    #remove all with a dist of less than 200ms?
    indDiff = diff(indices)
    for ii = 1:length(indDiff)
        if indDiff[ii] <= window
            if vals[ii] > vals[ii+1]
                vals[ii+1] = 0
            else
                vals[ii] = 0
            end
        end
    end
    deleteat!(indices, findall(x -> x == 0, vals))

    return (indices)
end


text:
using .TimeZones

# We're just adding one timezone specific method for constructing and formatting a ZonedDateTime
function _format_datetime(dt::DateTime, fmt::AbstractString, tz::Dates.TimeZone)
    return Dates.format(astimezone(ZonedDateTime(dt, localzone()), tz), fmt)
end

# Since we already default to local zone in the base case explicitly convert a
# ZonedDateTime to be consistent.
function _format_datetime(zdt::ZonedDateTime, fmt::AbstractString, tz::Nothing)
    return Dates.format(astimezone(zdt, localzone()), fmt)
end

# Convert ZonedDateTimes to a specified output tz
function _format_datetime(zdt::ZonedDateTime, fmt::AbstractString, tz::Dates.TimeZone)
    return Dates.format(astimezone(zdt, tz), fmt)
end


text:
function fetch_project(repo::LibGit2.GitRepo)
   for remote in LibGit2.remotes(repo)
      LibGit2.fetch(repo; remote=remote)
   end
end

function add_remote(repo::LibGit2.GitRepo; fork=nothing)
   pkg = basename(LibGit2.path(repo))
   remote = isnothing(fork) ? "origin" : fork
   isnothing(LibGit2.lookup_remote(repo, remote)) || return

   @info "$pkg: adding remote '$remote'"
   url = pkg_url(pkg, fork=fork)
   if !github_repo_exists(pkg,fork)
      @warn "$pkg: no such repository at org '$fork'"
      return
   end
   LibGit2.GitRemote(repo, remote, url,
                              "+refs/heads/*:refs/remotes/$remote/*")
   LibGit2.set_remote_push_url(repo, remote, pkg_giturl(pkg, fork=fork))
   LibGit2.fetch(repo; remote=remote)
end

add_remote(pkgdir::AbstractString; kwargs...) =
   add_remote(LibGit2.GitRepo(pkgdir); kwargs...)

function create_tracking_branch(repo::LibGit2.GitRepo, branch::AbstractString; remote::AbstractString="origin")
   pkg = basename(LibGit2.path(repo))
   @info "$pkg: creating branch '$branch' tracking from '$remote'"
   hash = LibGit2.GitHash(repo, "refs/remotes/$remote/$branch")
   LibGit2.branch!(repo, branch, string(hash); track="$remote")
   # for some reason LibGit2.jl does not set the correct upstream branch...
   LibGit2.with(LibGit2.GitConfig, repo) do cfg
      LibGit2.set!(cfg, "branch.$branch.remote", remote)
   end
end

function create_branch(repo::LibGit2.GitRepo, branch::AbstractString; start="origin/master")
   pkg = basename(LibGit2.path(repo))
   @info "$pkg: creating branch '$branch' (based on '$start')"
   fetch_project(repo)
   commit = start == "HEAD" ?
      LibGit2.head(repo) :
      LibGit2.lookup_branch(repo,"$start")
   isnothing(commit) && (commit = LibGit2.lookup_branch(repo,"$start", true))
   if !isnothing(commit)
      LibGit2.branch!(repo, branch, string(LibGit2.GitHash(commit)))
   else
      @warn "$pkg: start point '$start' not found"
   end
end

create_branch(pkgdir::AbstractString, branch; kwargs...) = create_branch(LibGit2.GitRepo(pkgdir), branch; kwargs...)

function update_project(repo::LibGit2.GitRepo)
   pkg = basename(LibGit2.path(repo))
   head = LibGit2.head(repo)
   upstream = LibGit2.upstream(head)
   if !isnothing(upstream)
      @info "$pkg: updating from '$(LibGit2.shortname(upstream))'"
      upstream = LibGit2.GitAnnotated(repo, upstream)
      fetch_project(repo)
      LibGit2.ffmerge!(repo, upstream)
   else
      @warn "$pkg: branch '$(LibGit2.shortname(head))' has no upstream"
   end
end

update_project(pkgdir::AbstractString) = update_project(LibGit2.GitRepo(pkgdir))

function checkout_branch(pkg::AbstractString, devdir::AbstractString; branch::AbstractString, fork=nothing, url=nothing)
   @info "$pkg: checkout branch '$branch'"
   # libgit2 is quite confusing...
   repo = LibGit2.GitRepo(devdir)
   LibGit2.fetch(repo; remote="origin")
   remote = LibGit2.lookup_remote(repo, isnothing(fork) ? "origin" : fork)
   # add new remote for fork if necessary
   if !isnothing(fork)
      add_remote(repo; fork=fork)
   end
   remote = LibGit2.name(remote)
   # check for branch
   commit = LibGit2.lookup_branch(repo,"$branch")
   if commit != nothing
      # switch branch and update
      LibGit2.branch!(repo, branch)
      update_project(devdir)
   else
      # create new tracking branch
      create_tracking_branch(repo, branch; remote=remote)
   end
end

function clone_project(pkg::AbstractString, devdir::AbstractString; branch::AbstractString, fork=nothing, url=nothing)
   @info "$pkg: cloning to '$devdir' (branch '$branch'):"
   if !isnothing(fork)
      # we create the main repo as well as the fork as upstream
      repo = LibGit2.clone(pkg_url(pkg; full=true), devdir)
      add_remote(repo; fork=fork)
      create_tracking_branch(repo, branch; remote=fork)
   else
      repo = LibGit2.clone(url, devdir, branch=branch)
   end
   LibGit2.set_remote_push_url(repo,"origin", pkg_giturl(pkg))
end

function checkout_project(pkg::AbstractString, dir::AbstractString; branch::AbstractString="master", fork=nothing)
   (url, branch, fork) = find_branch(pkg, branch; fork=fork)
   devdir = abspath(joinpath(dir, pkg))
   if isdir(devdir)
      checkout_branch(pkg, devdir; branch=branch, fork=fork, url=url)
   else
      clone_project(pkg, devdir; branch=branch, fork=fork, url=url)
   end
   return devdir
end

function pkg_array_to_dict(pkgs::Array{String})
   pkgdict = Dict{String,Any}()
   for pkg in pkgs
      pkgre = match(r"(\w+)(?:#(.*))?", pkg)
      push!(pkgdict, pkgre[1] => pkgre[2])
   end
   return pkgdict
end

function oscar_develop(pkgs::Dict{String,Any}; dir=default_dev_dir, branch::AbstractString="master", fork=nothing, active_repo=nothing)
   mkpath(dir)
   active_pkg = pkg_from_repo(active_repo)
   withenv("JULIA_PKG_DEVDIR"=>"$dir") do
      Pkg.activate(joinpath(dir,"project")) do
         @info "populating development directory '$dir':"
         if !isnothing(active_pkg)
            # during CI we always need to dev the active checkout
            @info "  reusing current dir for $active_pkg"
            Pkg.develop(path=".")
         end
         for (pkg, pkgbranch) in pkgs
            if pkg === active_pkg
               continue
            else
               if pkgbranch == "release"
                  # make sure we have that package added explicitly
                  Pkg.add(pkg)
               else
                  isnothing(pkgbranch) && (pkgbranch=branch)
                  devdir = checkout_project(pkg, dir; branch=pkgbranch, fork=fork)
                  Pkg.develop(path=devdir)
               end
            end
         end
      end
   end
   @info "Please start julia with:\njulia --project=$(abspath(joinpath(dir,"project")))"
end

oscar_develop(pkgs::Array{String}; kwargs...) =
   oscar_develop(pkg_array_to_dict(pkgs); kwargs...)

function oscar_update(; dir=default_dev_dir)
   for pkgdir in joinpath.(dir, pkg_names(dir))
      if isdir(joinpath(pkgdir,".git"))
         update_project(pkgdir)
      end
   end
   # update project after switching branches
   Pkg.activate(joinpath(dir,"project")) do
      Pkg.update()
   end
end

function oscar_add_remotes(pkgs::Array{String}, fork::AbstractString; dir=default_dev_dir)
   for pkgdir in joinpath.(dir, pkgs)
      if isdir(joinpath(pkgdir, ".git"))
         add_remote(pkgdir; fork=fork)
      end
   end
end

oscar_add_remotes(fork::AbstractString; dir=default_dev_dir) =
   oscar_add_remotes(pkg_names(dir), fork; dir=dir)

function oscar_branch(pkgs::Array{String}, branch::AbstractString; dir=default_dev_dir, start="origin/master")
   # if path is not there, develop it first
   if !ispath(dir)
     oscar_develop(pkgs; dir=dir)
   end
   for pkgdir in joinpath.(dir, pkgs)
      if isdir(joinpath(pkgdir, ".git"))
         create_branch(pkgdir, branch; start=start)
      end
   end
   # update project after switching branches
   Pkg.activate(joinpath(dir,"project")) do
      Pkg.update()
   end
end

oscar_branch(branch::AbstractString; dir=default_dev_dir, start="origin/master") =
   oscar_branch(pkg_names(dir), branch; dir=dir, start=start)


text:
using WingBiology

using JSON

results = Dict()

for g in (4e-3, 7e-3, 1e-2)
    c = 1.8288
    b = 2 * 6.096

    e = 0.33 - 0.25
    xCG = 0.43 - 0.25

    m = 35.71
    Ixx = 8.64

    GJ = 0.99e6
    EIxx = 9.77e6

    ρ = 1.020
    α = 0.05

    acft = Aircraft()

    ls = WingSectInfo(
        [0.0, - b / 2, 0.0],
        c;
        e = e,
        xCG = xCG,
        m = m,
        Ixx = Ixx,
        EIy = EIxx,
        GJ = GJ,
        g = g
    )
    cs = WingSectInfo(
        [0.0, 0.0, 0.0],
        c;
        e = e,
        xCG = xCG,
        m = m,
        Ixx = Ixx,
        EIy = EIxx,
        GJ = GJ,
        g = g
    )
    rs = WingSectInfo(
        [0.0, b / 2, 0.0],
        c;
        e = e,
        xCG = xCG,
        m = m,
        Ixx = Ixx,
        EIy = EIxx,
        GJ = GJ,
        g = g
    )

    lwing_pts, lwing_sects = interpolate!(acft, 30, ls, cs)
    rwing_pts, rwing_sects = interpolate!(acft, 30, cs, rs)

    add_link!(acft, rwing_pts[1])

    gres = Dict()
    results[
        "$g"
    ] = gres

    for U∞ = 100.0:5.0:220.0
        @show g, U∞

        vres = Dict()
        gres[
            "$U∞"
        ] = vres

        q = get_state(acft, U∞; α = α)

        #=
        Arnoldi_solve!(acft, q, U∞; α = α, ρ = ρ)

        plot_aircraft(acft, q)
        =#

        A, info = assumed_modes_solve(acft, q, U∞; ρ = ρ, α = α, n_eig = 4)

        λs = info.λ

        #=
        for (i, λ) in enumerate(λs)
            @show i, λ
        end
        =#

        vres["real_parts"] = real.(λs)
        vres["imaginary_parts"] = imag.(λs)
    end
end

let fobj = open("goland_results.json", "w")
    JSON.print(fobj, results)

    close(fobj)
end

max_stars_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_issues_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_forks_repo_name: Equipe-AeroDesign-ITA/WingBiology


text:
# Parser functionality for folders adhering to the
# ETSF Nanoquanta file format, for details see http://www.etsf.eu/fileformats/

export EtsfFolder

import JSON

struct EtsfFolder
    folder
    gsr
    den
    eig
    pspmap
end

"""
Initialize a EtsfFolder from the path to the folder which contains
the data in the ETSF Nanoquanta format.
"""
function EtsfFolder(folder::AbstractString)
    if !isfile(joinpath(folder, "out_GSR.nc"))
        error("Did not find file $folder/out_GSR.nc")
    end
    gsr = Dataset(joinpath(folder, "out_GSR.nc"))

    den = nothing
    eig = nothing
    pspmap = Dict{Int, Any}()

    if isfile(joinpath(folder, "out_DEN.nc"))
        den = Dataset(joinpath(folder, "out_DEN.nc"))
    end
    if isfile(joinpath(folder, "out_EIG.nc"))
        eig = Dataset(joinpath(folder, "out_EIG.nc"))
    end
    if isfile(joinpath(folder, "pspmap.json"))
        dict = JSON.parsefile(joinpath(folder, "pspmap.json"))
        pspmap = Dict{Int, Any}(parse(Int, k) => dict[k] for k in keys(dict))
    end

    EtsfFolder(folder, gsr, den, eig, pspmap)
end

load_lattice(T, folder::EtsfFolder) = Mat3{T}(folder.gsr["primitive_vectors"][:])

"""
Load a DFTK-compatible atoms object from the ETSF folder.
Use the scalar type `T` to represent the data.
"""
function load_atoms(T, folder::EtsfFolder)
    atoms = Dict{ElementPsp, Vector{Vec3{T}}}()
    n_species = length(folder.gsr["atomic_numbers"])
    for ispec in 1:n_species
        atnum = Int(folder.gsr["atomic_numbers"][ispec])
        spec = ElementPsp(atnum, psp=load_psp(folder.pspmap[atnum]))

        mask_species = findall(isequal(ispec), folder.gsr["atom_species"][:])
        positions = folder.gsr["reduced_atom_positions"][:, mask_species]
        atoms[spec] = [Vec3{T}(positions[:, m]) for m in mask_species]
    end
    collect(pairs(atoms))
end

"""
Load a DFTK-compatible model object from the ETSF folder.
Use the scalar type `T` to represent the data.
"""
function load_model(T, folder::EtsfFolder; magnetic_moments=[])
    # Parse functional ...
    functional = Vector{Symbol}()
    ixc = Int(folder.gsr["ixc"][:])
    if ixc == 0
        # No XC functional
    elseif ixc == 1 || ixc == -20
        push!(functional, :lda_xc_teter93)
    elseif ixc == -1007
        push!(functional, :lda_x, :lda_c_vwn)
    elseif ixc == -101130 || ixc == 11
        push!(functional, :gga_x_pbe, :gga_c_pbe)
    else
        error("Unknown Functional ixc==$ixc")
    end

    # Parse smearing-scheme ...
    smearing = strip(join(folder.gsr["smearing_scheme"][:]))
    Tsmear = 0
    smearing_function = nothing
    if smearing == "none"
        # No smearing
    elseif smearing == "Fermi-Dirac"
        smearing_function = Smearing.FermiDirac()
    elseif smearing == "Methfessel and Paxton"
        smearing_function = Smearing.MethfesselPaxton1()
    elseif smearing == "gaussian"
        smearing_function = Smearing.Gaussian()
    else
        error("Unknown Smearing scheme: $smearing")
    end
    smearing_function !== nothing && (Tsmear = folder.gsr["smearing_width"][:])

    n_spin = size(folder.gsr["eigenvalues"], 3)
    spin_polarization = n_spin == 2 ? :collinear : :none
    if spin_polarization == :collinear && isempty(magnetic_moments)
        @warn("load_basis will most likely fail for collinear spin if magnetic_moments " *
              "not explicitly specified.")
    end

    # Build model
    lattice = load_lattice(T, folder)
    atoms = load_atoms(T, folder)
    model_DFT(Array{T}(lattice), atoms, functional;
              smearing=smearing_function, temperature=Tsmear,
              spin_polarization=spin_polarization, magnetic_moments=magnetic_moments)
end

"""
Load a DFTK-compatible basis object from the ETSF folder.
Use the scalar type `T` to represent the data.
"""
function load_basis(T, folder::EtsfFolder; magnetic_moments=[])
    model = load_model(T, folder, magnetic_moments=magnetic_moments)
    atoms = load_atoms(T, folder)

    Ecut = folder.gsr["kinetic_energy_cutoff"][:]
    kcoords = Vec3{T}.(eachcol(folder.gsr["reduced_coordinates_of_kpoints"]))

    # Try to determine whether this is a shifted kpoint mesh or not
    if length(kcoords) > 1
        ksmallest = sort(filter(k -> all(k .≥ 0), kcoords), by=norm)[1]
    else
        ksmallest = kcoords[1]
    end
    kshift = [(abs(k) > 10eps(T)) ? 0.5 : 0.0 for k in ksmallest]

    kgrid = Vector{Int}(folder.gsr["monkhorst_pack_folding"])
    kcoords_new, ksymops, _ = bzmesh_ir_wedge(kgrid, model.symmetries, kshift=kshift)
    @assert kcoords_new ≈ normalize_kpoint_coordinate.(kcoords)

    fft_size = size(folder.den["density"])[2:4]
    PlaneWaveBasis(model, Ecut, kcoords, ksymops, fft_size=fft_size,
                   kgrid=kgrid, kshift=kshift)
end

"""
Load a DFTK-compatible density object from the ETSF folder.
Use the scalar type `T` to represent the data.
"""
function load_density(T, folder::EtsfFolder; magnetic_moments=[])
    basis = load_basis(T, folder, magnetic_moments=magnetic_moments)

    n_real_cpx = size(folder.den["density"], 1)
    @assert n_real_cpx == 1

    n_spin = size(folder.den["density"], 5)
    @assert n_spin in (1, 2)

    if n_spin == 1
        @assert basis.model.spin_polarization == :none
        ρtot_real = Array{T}(folder.den["density"][1, :, :, :, 1])
        @assert basis.fft_size == size(ρ_real)

        return ρ_from_total_and_spin(ρtot_real)
    else
        @assert basis.model.spin_polarization == :collinear
        ρtot_real = Array{T}(folder.den["density"][1, :, :, :, 1])
        ρα_real = Array{T}(folder.den["density"][1, :, :, :, 2])
        @assert basis.fft_size == size(ρtot_real)
        @assert basis.fft_size == size(ρα_real)

        return ρ_from_total_and_spin(ρtot_real, 2ρα_real - ρtot_real)
    end
end


text:
using Revise
using RowTables
using DataFrames
using BenchmarkTools

# using DataFramesUtils

# Same example as in runtests
# ds = [
#     Dict(:a=>2,:b=>1,:c=>7)
#     Dict(:a=>10,:b=>4,:c=>7)
#     Dict(:a=>1,:b=>4,:c=>8)
#     Dict(:a=>10,:b=>4,:c=>1)
#     Dict(:a=>9,:b=>6,:c=>4)
#     Dict(:a=>10,:b=>4,:c=>3)]

# rt = RowTable(ds,[:a,:b,:c])
# df = DataFrame(rt)

# df2 = DataFramesUtils.gentestdf()
# rt2 = RowTable(df2)

rt3 = RowTable([[rand(),rand(),rand()] for i in 1:10^3], [:a,:b,:c]);

function mkbigdf(nr,nc)
    cols = Any[]
    for i in 1:nc
        push!(cols,rand(nr))
    end
    syms = [Symbol("x",i) for i in 1:nc]
    DataFrame(cols,syms)
end

function randsort!(obj::RowTable)
    sort!(obj,cols=[rand(1:size(obj,2))])
    nothing
end

function randsort!(obj::DataFrame)
    sort!(obj, [rand(1:size(obj,2))])
    nothing
end

function mkobjs()
    df = DataFrame(A = [45], B = ["F"])
    rt = RowTable(df)
    return (df,rt)
end

function appendrows(obj,n)
    nobj = copy(obj)
    for i in 1:n
        push!(nobj,[32, "F"])
    end
    nobj
end

nothing


text:
Base.@kwdef struct ProgressBar
    block::Block
    ratio::Float64
    function ProgressBar(block, ratio)
        ( ratio < 0 || ratio > 1 ) && error("Got $ratio. ProgressBar ratio must be in [0, 1].")
        new(block, ratio)
    end
end

function draw(pg::ProgressBar, rect::Rect, buf::Buffer)
    draw(pg.block, rect, buf)

    inner_area = inner(pg.block, rect)

    center = height(inner_area) ÷ 2 + top(inner_area)

    crayon = Crayon(foreground = (255, 255, 255), background = (0, 0, 0))

    for y in top(inner_area):bottom(inner_area)
        for x in left(inner_area):right(inner_area)
            if x <= pg.ratio * (right(inner_area) - left(inner_area) + 1)
                set(buf, x, y, crayon)
            else
                set(buf, x, y, inv(crayon))
            end
        end
    end

    label = "$(round(pg.ratio * 100))%"
    middle = Int(round( (width(inner_area) - length(label)) / 2 + left(inner_area) ))
    set(buf, middle, center, label)

end


text:
using Dates, NanoDates
using NanoDates
using Dates: value, CompoundPeriod, toms, tons, UTM, UTD
using NanoDates: nanosecs, NanoDate0, DateTime0, Date0, Time0
using Aqua, Test

Aqua.test_ambiguities([NanoDates, Base, Core])
Aqua.test_unbound_args(NanoDates)
Aqua.test_undefined_exports(NanoDates)
Aqua.test_project_extras(NanoDates)
Aqua.test_stale_deps(NanoDates; ignore=[:Aqua])
Aqua.test_deps_compat(NanoDates)

include("constants.jl")

@testset "nanodate" begin
    include("nanodate.jl")
end

@testset "accessors" begin
    include("accessors.jl")
end

@testset "convert" begin
    include("convert.jl")
end

@testset "compare" begin
    include("compare.jl")
end

@testset "interop" begin
    include("interop.jl")
end

@testset "arith" begin
    include("arith.jl")
end

@testset "strings" begin
    include("strings.jl")
end


text:
using DifferentiableStateSpaceModels: vech, inv_vech
using Random
const UTILITIES = BenchmarkGroup()
Random.seed!(1234)
const A = rand(10, 10)
const Asym = A' * A
const vech_A = vech(Asym)
UTILITIES["vech"] = @benchmarkable vech($Asym)
UTILITIES["inv_vech"] = @benchmarkable inv_vech($vech_A)

# return for the suite
UTILITIES

# Uncomment to run on its own, otherwise nests into full benchmarking
# results = run(UTILITIES; verbose = true)
# results["vech"] # etc.


text:
export @comment, @%, @import, @do, @until, @def

macro comment(expr) end
macro %(expr) end

@eval begin
    "@import Foo as f"
    macro $(:import)(mod::Symbol, as::Symbol, alias::Symbol)
        as ≠ :as && error("@import expected `as` got `$as`")
        Expr(:toplevel, Expr(:import, mod), :(const $alias = $mod)) |> esc
    end

    "do-when (do-while) loop"
    macro $(:do)(block, when::Symbol, condition)
        when ≠ :when && error("@do expected `when` got `$s`")
        quote
            let
                $block
                while $condition
                    $block
                end
            end
        end |> esc
    end
end

"until loop"
macro until(condition, block)
    quote
        while !$condition
            $block
        end
    end |> esc
end

macro def(name, definition)
    return quote
        macro $(esc(name))()
            esc($(quot(definition)))
        end
    end
end


text:
using AdversarialGridworld

## Generate training and test sets
N_train = 30
seed = 0
tprob_test = 0.7
folder = "data"

pwd()
cd(folder)

# Generate training tasks
train_tasks = generate_task_set(seed, N_train, tprob_test = tprob_test, folder = "diffmaps", save = true, keep_map = false, verbose = true)


text:
using VegaStreams
using Test

@testset "VegaStreams.jl" begin
    vls = vegastream()
    for row in enumerate(randn(100))
        sleep(0.01)
        push!(vls, row)
    end

    vls = vegastream(:point)
    for (x, y) in enumerate(randn(100))
        sleep(0.01)
        push!(vls, (x=x, y=y))
    end
end


text:
using ModelingToolkit, MethodOfLines, LinearAlgebra, OrdinaryDiffEq
using DomainSets

using Plots

local sol
@time begin #@testset "Test 01: Brusselator equation 2D" begin
       @parameters x y t
       @variables u(..) v(..)
       Dt = Differential(t)
       Dx = Differential(x)
       Dy = Differential(y)
       Dxx = Differential(x)^2
       Dyy = Differential(y)^2

       brusselator_f(x, y, t) = (((x-0.3)^2 + (y-0.6)^2) <= 0.1^2) * (t >= 1.1) * 5.

       x_min = y_min = t_min = 0.0
       x_max = y_max = 1.0
       t_max = 11.5

       α = 10.

       u0(x,y,t) = 22(y*(1-y))^(3/2)
       v0(x,y,t) = 27(x*(1-x))^(3/2)

       eq = [Dt(u(x,y,t)) ~ 1. + v(x,y,t)*u(x,y,t)^2 - 4.4*u(x,y,t) + α*(Dxx(u(x,y,t)) + Dyy(u(x,y,t))) + brusselator_f(x, y, t),
             Dt(v(x,y,t)) ~ 3.4*u(x,y,t) - v(x,y,t)*u(x,y,t)^2 + α*(Dxx(v(x,y,t)) + Dyy(v(x,y,t)))]

       domains = [x ∈ Interval(x_min, x_max),
                  y ∈ Interval(y_min, y_max),
                  t ∈ Interval(t_min, t_max)]

       bcs = [u(x,y,0) ~ u0(x,y,0),
              u(0,y,t) ~ u(1,y,t),
              u(x,0,t) ~ u(x,1,t),

              v(x,y,0) ~ v0(x,y,0),
              v(0,y,t) ~ v(1,y,t),
              v(x,0,t) ~ v(x,1,t)]

       @named pdesys = PDESystem(eq,bcs,domains,[x,y,t],[u(x,y,t),v(x,y,t)])

       # Method of lines discretization
       dx = 1/4
       dy = 1/4

       order = 2

       discretization = MOLFiniteDifference([x=>dx, y=>dy], t, approx_order=order)

       # Convert the PDE problem into an ODE problem
       generate_code(pdesys,discretization, "brusselator_2D_ode.jl")
end


text:
"""
Parsers that do not consume any input can inherit `Assertion{S,T}`.
!!! note
    TODO: allow to keep state and return wrapped get
"""
abstract type Assertion{S,T} <: CombinedParser{S,T} end
@inline _leftof(str,i,parser::Assertion,x...) = i
@inline _rightof(str,i,parser::Assertion,x...) = i

"""
    Base.get(parser::Assertion{MatchState, <:Assertion}, sequence, till, after, i, state)

Most assertions return the assertion parser as a result
([`AtStart`](@ref), [`AtEnd`](@ref),
[`Always`](@ref), [`Never`](@ref),
[`NegativeLookahead`](@ref), [`NegativeLookbehind`](@ref)).
"""
Base.get(parser::Assertion{MatchState, <:Assertion}, sequence, till, after, i, state) =
    parser

export AtStart, AtEnd
"""
    AtStart()

Parser succeding if and only if at index 1 with `result_type` `AtStart`.

```jldoctest
julia> AtStart()
re"^"

```
"""
struct AtStart <: Assertion{MatchState,AtStart} end
regex_inner(x::AtStart) = "^"
iterate_state(parser::AtStart, sequence, till, posi, next_i, state::Nothing) =
    next_i == 1 ? (next_i, MatchState()) : nothing

print_constructor(io::IO, x::AtStart) = print(io,"AtStart")

"""
    AtEnd()

Parser succeding if and only if at last index with `result_type` `AtEnd`.

```jldoctest
julia> AtEnd()
re"\$"

```
"""
struct AtEnd <: Assertion{MatchState,AtEnd} end
regex_inner(x::AtEnd) = "\$"
iterate_state(parser::AtEnd, sequence, till, posi, next_i, state::Nothing) =
    next_i > till ? (next_i, MatchState()) : nothing
print_constructor(io::IO, x::AtEnd) = print(io,"AtEnd")

export Never
"""
    Never()

Assertion parser matching never.

```jldoctest
julia> Never()
re"(*FAIL)"

```
"""
struct Never <: Assertion{MatchState,Never} end
regex_prefix(x::Never) = "(*"
regex_inner(x::Never) = "FAIL"
regex_suffix(x::Never) = ")"
iterate_state(x::Never,str,posi, next_i,till,state::Nothing) =
    nothing

export Always
"""
    Always()

Assertion parser matching always and not consuming any input.
Returns `Always()`.

```jldoctest
julia> Always()
re""

```
"""
struct Always <: Assertion{MatchState,Always}
end
Base.show(io::IO,x::Always) = print(io,"re\"\"")
children(x::Union{Never,Always}) = tuple()
regex_prefix(x::Always) = ""
regex_inner(x::Always) = ""
regex_suffix(x::Always) = ""
iterate_state(parser::Always, str, till, posi, next_i, s::Nothing) =
    next_i, MatchState()

Base.show(io::IO, x::Union{AtStart,AtEnd,Never,Always}) =
    print(io,"re\"",regex_string(x),"\"")

@inline iterate_state(t::Union{AtStart,AtEnd,Never,Always}, str, till, posi, next_i, state::MatchState) = nothing

"""
An assertion with an inner parser, like WrappedParser interface.
"""
abstract type WrappedAssertion{S,T} <: Assertion{S,T} end
children(x::WrappedAssertion) = children(x.parser)
regex_suffix(x::WrappedAssertion) = regex_suffix(x.parser)*")"
regex_inner(x::WrappedAssertion) = regex_inner(x.parser)

function print_constructor(io::IO,x::WrappedAssertion)
    if !hasregex(x.parser)
        print_constructor(io, x.parser)
        print(io, " |> ")
    end
    print(io, constructor_name(x))
end

export PositiveLookahead
"""
    PositiveLookahead(parser)

Parser that succeeds if and only if `parser` succeeds, but consumes no input.
The match is returned.
Useful for checks like "must be followed by `parser`, but don't consume its match".

```jldoctest
julia> la=PositiveLookahead("peek")
re"(?=peek)"

julia> parse(la*AnyChar(),"peek")
("peek", 'p')

```
"""
@auto_hash_equals struct PositiveLookahead{S,T,P} <: WrappedAssertion{S,T}
    parser::P
    PositiveLookahead(p_,reversed=true) =
        let p = parser(p_)
            new{Tuple{Int,state_type(p)},result_type(p),typeof(p)}(p)
        end
end
regex_prefix(x::PositiveLookahead) = "(?="*regex_prefix(x.parser)

function iterate_state(t::PositiveLookahead, str, till, posi, next_i, state)
    r = iterate_state(t.parser, str, till, posi, tuple_pos(state,posi), tuple_state(state))
    if r === nothing
        nothing
    else
        next_i, r
    end
end
function Base.get(parser::PositiveLookahead, sequence, till, after, i, state)
    after_ = tuple_pos(state)
    get(parser.parser, sequence, till, after_, i, tuple_state(state))
end

export NegativeLookahead
"""
    NegativeLookahead(parser)

Parser that succeeds if and only if `parser` does not succeed, but consumes no input.
`parser` is returned as match.
Useful for checks like "must not be followed by `parser`, don't consume its match".

```jldoctest
julia> la = NegativeLookahead("peek")
re"(?!peek)"

julia> parse(la*AnyChar(),"seek")
(re"(?!peek)", 's')

```
"""
@auto_hash_equals struct NegativeLookahead{P} <: WrappedAssertion{MatchState,NegativeLookahead{P}}
    parser::P
    NegativeLookahead(p_,reversed=true) =
        let p = parser(p_)
            new{typeof(p)}(p)
        end
end
regex_prefix(x::NegativeLookahead) = "(?!"*regex_prefix(x.parser)
function iterate_state(t::NegativeLookahead, str, till, posi, next_i, state::Nothing)
    r = iterate_state(t.parser, str, till, posi, next_i, nothing)
    if r === nothing
        next_i,MatchState()
    else
        nothing
    end
end

@inline iterate_state(t::NegativeLookahead, str, till, posi, next_i, state::MatchState) = nothing

export Lookahead

"""
    Lookahead(does_match::Bool, p)

[`PositiveLookahead`](@ref) if `does_match==true`,
[`NegativeLookahead`](@ref) otherwise.
"""
function Lookahead(does_match::Bool, p_)
    p = parser(p_)
    if does_match
        PositiveLookahead(p)
    else
        NegativeLookahead(p)
    end
end

@deprecate look_ahead(does_match,p) Lookahead(does_match, p)