text:
import StatsBase.predict

function predict(model::UnfoldLinearModel, newdata)
    # make a copy of it so we don't change it outside the function
     data = copy(newdata)
    #data = newdata
    # if only one formulas
    if typeof(model.X.formulas) <: FormulaTerm
        formulas = [model.X.formulas]
    else
        formulas = model.X.formulas
    end
    # add those as latencies [xxx change to basisfunction field] to newdata
    ##

    X = []
    fromTo = []
    times = []
    if !(typeof(formulas[1].rhs) <: unfold.TimeExpandedTerm)
        # mass univariate model. Not sure how I can multiple dispatch correctly :| Maybe I need 4 types after all
        X = designmatrix(UnfoldLinearModel,formulas[1],data)
        yhat = Array{Union{Missing,Float64}}(missing,size(model.beta,1),size(X.Xs,1),size(model.beta,2))
        #yhat = zeros(size(model.beta,1),size(X.Xs,1),size(model.beta,2))
        for ch = 1:size(model.beta,1)
            yhat[ch,:,:] = X.Xs * permutedims(model.beta[ch,:,:],(2,1))
        end
        # fake times because we never save it in the object
        times = repeat(range(1,length=size(yhat,3)),size(data,1))
        # reshape yhat

        yhat = reshape(permutedims(yhat,(1,3,2)),size(yhat,1),:)
        yhat = permutedims(yhat,(2,1))
    else

    for f = formulas

        # find out how long each designmatrix is
        n_range = length(f.rhs.basisfunction.times)
        # find out how much to shift so that X[1,:] is the the first "sample"
        n_negative = f.rhs.basisfunction.shiftOnset
        # generate the correct eventfields (default: latencies)
        data[:,f.rhs.eventfields[1]] = range(-n_negative + 1, step=n_range, length=size(data, 1))

        # get the model
        Xsingle = modelcols(f.rhs, data)
        # combine designmats
        append!(X, [Xsingle])

        # keep track of how long each event is
        append!(fromTo, [range(1, step=n_range, length=1+size(data, 1))])

        # keep track of the times
        append!(times, repeat(f.rhs.basisfunction.times[1:end], size(data, 1)))

    end

    # Concat them, but without introducing any overlap, we want the isolated responses
    Xconcat = blockdiag(X...)

    # calculate yhat
    # yhat x channels
    yhat =  Xconcat * model.beta'
    end

    # init the meta dataframe
    metaData = DataFrame([:times => vcat(times...),:basisname=>""])
    for c = names(newdata)
        metaData[:,c] .= newdata[1,c] # assign first element in order to have same column type
    end

    if !(typeof(formulas[1].rhs) <: unfold.TimeExpandedTerm)
        # for mass univariate we can make use of the knowledge that all events have the same length :)
        ntimes = size(model.beta,2)
        for c  = names(newdata)
            for row = 1:size(data,1)
                rowIx = (1. .+(row-1).*ntimes) .+ range(1.,length=ntimes) .-1
                metaData[ Int64.(rowIx)  ,c] .= data[row,c]
            end
        end

    else

    # shift variable to keep track of multiple basisfunctions
    shift = 0
    # for each basis function
    for (bIx,basisfun) = enumerate(fromTo)
        # go through all predictors
        for (i, fstart) in enumerate(basisfun[1:end-1])

            # the last one is skipped in the enumerate above
            fend = basisfun[i + 1]-1

            # couldn't figure out how to broadcast everything directly (i.e out[fstart:fend,names(newdata)] .= newdata[i,:])
            # copy the correct metadata
            for j = fstart:fend
                metaData[j+shift,names(newdata)] = newdata[i,:]
            end
            # add basisfunction name
            metaData[shift.+(fstart:fend),:basisname] .= formulas[bIx].rhs.basisfunction.name
        end
        # the next meta data has to be at the end
        shift += basisfun[end]-1

    end
    end

    out = DataFrame([:yhat => vec(reshape(Float64.(yhat), :, 1))])
    nchannel = size(yhat, 2)

    out.channel = repeat(1:nchannel, inner=size(yhat, 1))

    out = hcat(out, repeat(metaData, nchannel))
    return out
end


text:
# Workaround for GR warnings
ENV["GKSwstype"] = "100"

using Documenter, ImageQuilting

isCI = "CI" ∈ keys(ENV)

makedocs(
  format = Documenter.HTML(prettyurls=isCI),
  sitename = "ImageQuilting.jl",
  authors = "Júlio Hoffimann",
  pages = [
    "Home" => "index.md",
    "Examples" => "examples.md",
    "Voxel reuse" => "voxel-reuse.md",
    "GPU support" => "gpu-support.md",
    "About" => [
      "Author" => "about/author.md",
      "License" => "about/license.md",
      "Citation" => "about/citation.md"
    ]
  ]
)

deploydocs(repo="github.com/JuliaEarth/ImageQuilting.jl.git")

ext: jl


text:
"""Convenience iterable for state-dependent controllers"""
struct ControllerIterator{C,E,B}
    controller::C
    env::E
    T::Int
    trajectory::B
    plotiter::Int
    randstart::Bool
    alreadyran::Base.RefValue{Bool}

    function ControllerIterator(
        controller,
        env::AbstractEnvironment;
        T::Integer = 1000,
        plotiter::Integer = 1,
        randstart::Bool = false
    )
        trajectory = (
            states = Array(undef, statespace(env), T),
            observations = Array(undef, observationspace(env), T),
            actions = Array(undef, actionspace(env), T),
            rewards = Array(undef, rewardspace(env), T),
        )
        new{typeof(controller),typeof(env),typeof(trajectory)}(
            controller,
            env,
            T,
            trajectory,
            plotiter,
            randstart,
            Ref(false)
        )
    end
end

function Base.iterate(it::ControllerIterator, t::Int = 1)
    it.alreadyran[] && error("Cannot iterate on a ControllerIterator more than once")
    if t > it.T
        it.alreadyran[] = true
        return
    end

    rolloutstep!(it.controller, it.trajectory, it.env, t)
    if mod(t, it.plotiter) == 0
        rew = @views Tools.Line(it.trajectory.rewards[1:t], "Reward")
        eval = @views Tools.Line(it.trajectory.evaluations[1:t], "Eval")
        plt = Tools.termplot(
            rew,
            eval,
            title = "ControllerIterator Iteration=$t/$(it.T)",
            width = 40,
        )
        display(plt)
    end

    (t, it), t + 1
end

Base.length(it::ControllerIterator) = it.T

function rolloutstep!(controller, traj, env, t)
    st = view(traj.states, .., t)
    ot = view(traj.observations, .., t)
    at = view(traj.actions, .., t)

    getstate!(st, env)
    getobservation!(ot, env)
    controller(at, st, ot)
    setaction!(env, at)

    step!(env)
    r = getreward(st, at, ot, env)
    done = isdone(st, ot, env)

    traj.rewards[t] = r

    traj
end


text:
using GHOST
setup()
licenses()


text:
##
# import Base: +, -, zero
import Base.Iterators
using Coverage, LinearAlgebra, StaticArrays, Zygote

# struct Point
#   x::Float64
#   y::Float64
# end
#
# width(p::Point) = p.x
# height(p::Point) = p.y
#
# a::Point + b::Point = Point(width(a) + width(b), height(a) + height(b))
# a::Point - b::Point = Point(width(a) - width(b), height(a) - height(b))
# dist(p::Point) = sqrt(width(p)^2 + height(p)^2)
#
# ##

@Zygote.adjoint (T :: Type{<:SVector})(xs :: Number ...) = T(xs...), dv -> (nothing, dv...)
@Zygote.adjoint (T :: Type{<:SVector})(x :: AbstractVector) = T(x), dv -> (nothing, dv)

@Zygote.adjoint enumerate(xs) = enumerate(xs), diys -> (map(last, diys),)

_ndims(::Base.HasShape{d}) where {d} = d
_ndims(x) = Base.IteratorSize(x) isa Base.HasShape ? _ndims(Base.IteratorSize(x)) : 1

@Zygote.adjoint function Iterators.product(xs...)
                    d = 1
                    Iterators.product(xs...), dy -> ntuple(length(xs)) do n
                        nd = _ndims(xs[n])
                        dims = ntuple(i -> i<d ? i : i+nd, ndims(dy)-nd)
                        d += nd
                        func = sum(y->y[n], dy; dims=dims)
                        ax = axes(xs[n])
                        reshape(func, ax)
                    end
                end

# @Zygote.adjoint width(p::Point) = p.x, x̄ -> (Point(x̄, 0),)
# @Zygote.adjoint height(p::Point) = p.y, ȳ -> (Point(0, ȳ),)
# @Zygote.adjoint Point(a, b) = Point(a, b), p̄ -> (p̄.x, p̄.y)
# zero(::Point) = Point(0, 0)

# ##
# xs = Point.(1:5, 5:9)
#
# function something(xs)
#     sum([ width(p1) + height(p2) for p1 in xs, p2 in xs ])
# end
#
# gradient(something, xs)

function sum2(op,arr)
    return sum(op,arr)
end

function sum2adj( Δ, op, arr )
    n = length(arr)
    g = x->Δ*Zygote.gradient(op,x)[1]
    return ( nothing, map(g,arr))
end

Zygote.@adjoint function sum2(op,arr)
    return sum2(op,arr),Δ->sum2adj(Δ,op,arr)
end

Zygote.refresh()

function KpG(kz::Real,gx::Real,gy::Real,gz::Real)
	# scale = ds.kpG[i,j,k].mag
	kpg = [-gx; -gy; kz-gz]
	mag = sqrt(sum2(abs2,kpg)) # norm(kpg)
	if mag==0
		n = [0.; 1.; 0.] # SVector(0.,1.,0.)
		m = [0.; 0.; 1.] # SVector(0.,0.,1.)
	else
		if kpg[1]==0. && kpg[2]==0.    # put n in the y direction if k+G is in z
			n = [0.; 1.; 0.] #SVector(0.,1.,0.)
		else                                # otherwise, let n = z x (k+G), normalized
			ntemp = [0.; 0.; 1.] × kpg  #SVector(0.,0.,1.) × kpg
			n = ntemp / sqrt(sum2(abs2,ntemp)) # norm(ntemp) #
		end
	end
	# m = n x (k+G), normalized
	mtemp = n × kpg
	m = mtemp / sqrt(sum2(abs2,mtemp)) #norm(mtemp) # sqrt( mtemp[1]^2 + mtemp[2]^2 + mtemp[3]^2 )
	return kpg, mag, m, n
end

function KpG2(kz::Float64,gx::Float64,gy::Float64,gz::Float64)
	# scale = ds.kpG[i,j,k].mag
	kpg = SVector{3,Float64}(-gx, -gy, kz-gz)
	mag = norm(kpg)
	if mag==0
		n = SVector{3,Float64}(0.,1.,0.)
		m = SVector{3,Float64}(0.,0.,1.)
	else
		if kpg[1]==0. && kpg[2]==0.    # put n in the y direction if k+G is in z
			n = SVector{3,Float64}(0.,1.,0.)
		else                                # otherwise, let n = z x (k+G), normalized
			ntemp = SVector{3,Float64}(0.,0.,1.) × kpg
			n = ntemp / norm(ntemp) #
		end
	end
	# m = n x (k+G), normalized
	mtemp = n × kpg
	m = mtemp / norm(mtemp) # sqrt( mtemp[1]^2 + mtemp[2]^2 + mtemp[3]^2 )
	return kpg, mag, m, n
end

# kpg, mag, m, n = KpG(1.5,3.5,6.5,0.0)
# KpG2(1.5,3.5,6.5,0.0)
# KpG(1.5,3.5,6.5,0.0)[2]
# KpG2(1.5,3.5,6.5,0.0)[2]
# Zygote.gradient((kz,gx,gy,gz)->KpG(kz,gx,gy,gz)[2],1.5,3.5,6.5,0.0)
# Zygote.gradient((kz,gx,gy,gz)->KpG2(kz,gx,gy,gz)[2],1.5,3.5,6.5,0.0)

function ∇KpG2(a,b,c,d)
	Zygote.gradient((kz,gx,gy,gz)->KpG(kz,gx,gy,gz)[2],a,b,c,d)
end

function ∇KpG22(a::Float64,b::Float64,c::Float64,d::Float64)
	Zygote.gradient((kz::Float64,gx::Float64,gy::Float64,gz::Float64)->KpG2(kz,gx,gy,gz)[2],a,b,c,d)
end

using ReverseDiff
function ∇KpG23(a::Float64,b::Float64,c::Float64,d::Float64)
	# ReverseDiff.gradient(input->KpG(input...)[2],[a,b,c,d])
	ReverseDiff.gradient!(results, compiled_f_tape, [a,b,c,d])
end

using ForwardDiff

ForwardDiff.gradient(input->KpG(input...)[2],[1.5,3.5,6.5,0.0])

function ∇KpG24(a::Real,b::Real,c::Real,d::Real)
	ForwardDiff.gradient(input->KpG(input...)[2],[a,b,c,d])
end

# function ∇KpG25(a::Real,b::Real,c::Real,d::Real)
# 	Zygote.gradient(a,b,c,d) do kz,gx,gy,gz
# 		Zygote.forwarddiff([kz,gx,gy,gz]) do (kz,gx,gy,gz)
# 			KpG(kz,gx,gy,gz)[2]
# 		end
# 	end
# end

∇KpG25(a,b,c,d) = Zygote.gradient(a,b,c,d) do kz,gx,gy,gz
	Zygote.forwarddiff([kz,gx,gy,gz]) do (kz,gx,gy,gz)
		KpG(kz,gx,gy,gz)[2]
	end
end

∇KpG2(1.5,3.5,6.5,0.0)
∇KpG22(1.5,3.5,6.5,0.0)
∇KpG23(1.5,3.5,6.5,0.0)
∇KpG24(1.5,3.5,6.5,0.0)
∇KpG25(1.5,3.5,6.5,0.0)

##
@btime ∇KpG2(1.5,3.5,6.5,0.0)
@btime ∇KpG22(1.5,3.5,6.5,0.0)
@btime ∇KpG23(1.5,3.5,6.5,0.0)
@btime ∇KpG24(1.5,3.5,6.5,0.0)
@btime ∇KpG25(1.5,3.5,6.5,0.0)
##

ReverseDiff.gradient(input->KpG(input...)[2],[1.5,3.5,6.5,0.0])

const f_tape = ReverseDiff.GradientTape(input->KpG(input...)[2], rand(4))
const compiled_f_tape = ReverseDiff.compile(f_tape)
results = similar(rand(4))
@btime ReverseDiff.gradient!(results, compiled_f_tape, [1.5,3.5,6.5,0.0])

##
using Tullio
Nz = 1
gx = collect(fftfreq(Nx,Nx/6.0))
gy = collect(fftfreq(Ny,Ny/4.0))
gz = collect(fftfreq(Nz,Nz/1.0))
kz = 1.476254662632502 #1.5
kpg_tuples = [KpG(kz,ggx,ggy,ggz) for ggx in gx, ggy in gy, ggz in gz]
kpg = [kpg_tuples[i,j,k][1][a] for a=1:3,i=1:Nx,j=1:Ny,k=1:Nz]
kpg_mag = [kpg_tuples[i,j,k][2] for i=1:Nx,j=1:Ny,k=1:Nz]
m = [kpg_tuples[i,j,k][3][a] for a=1:3,i=1:Nx,j=1:Ny,k=1:Nz]
n = [kpg_tuples[i,j,k][4][a] for a=1:3,i=1:Nx,j=1:Ny,k=1:Nz]
mn = [kpg_tuples[i,j,k][2+b][a] for a=1:3,b=1:2,i=1:Nx,j=1:Ny,k=1:Nz]
Ha = copy(reshape(H,(2,Nx,Ny,1)))
Hri = vcat(reim(reshape(Ha,(1,size(Ha)...)))...)
ε⁻¹

zxinds = [2; 1; 3]
zxscales = [-1.; 1.; 0.]
@tullio zxH[a,i,j,k] := zxscales[a] * Ha[b,i,j,k] * mn[zxinds[a],b,i,j,k] nograd=(zxscales,zxinds,m,n) verbose=true

using Zygote, ForwardDiff
function zx(Ha,mn)
	zxinds = [2; 1; 3]
	zxscales = [-1.; 1.; 0.]
	@tullio zxH[a,i,j,k] := zxscales[a] * Ha[b,i,j,k] * mn[zxinds[a],b,i,j,k] nograd=(zxscales,zxinds) #,m,n) #verbose=true # grad=Dual
end

zx(Ha,mn)
Zygote.gradient(abs2∘sum∘zx,Ha,mn)
@btime Zygote.gradient(abs2∘sum∘zx,$Ha,$mn)[1]

@btime fft(zx($Ha,$mn),(2:4))

@tullio e[a,i,j,k] :=  ε⁻¹[a,b,i,j,k] * fft(zx(Ha,mn),(2:4))[b,i,j,k] verbose=true

function ε⁻¹_op(ε⁻¹,Ha,mn)
	@tullio e[a,i,j,k] :=  ε⁻¹[a,b,i,j,k] * fft(zx(Ha,mn),(2:4))[b,i,j,k]
	return ifft(e,(2:4))
end

E = ε⁻¹_op(ε⁻¹,Ha,mn)
Zygote.gradient(abs2∘sum∘ε⁻¹_op,ε⁻¹,Ha,mn)
@btime Zygote.gradient(abs2∘sum∘ε⁻¹_op,$ε⁻¹,$Ha,$mn)

kxinds = [2; 1]
kxscales = [-1.; 1.]
@tullio kxE[b,i,j,k] := kxscales[b] * kpg_mag[i,j,k] * E[a,i,j,k] * mn[a,kxinds[b],i,j,k] nograd=(kxscales,kxinds) verbose=true
@tullio kxE[b,i,j,k] := kxscales[b] * kpg_mag[i,j,k] * ε⁻¹_op(ε⁻¹,Ha,mn)[a,i,j,k] * mn[a,kxinds[b],i,j,k] nograd=(kxscales,kxinds) verbose=true

function Mk(Ha,ε⁻¹,mn,kpg_mag)
	kxinds = [2; 1]
	kxscales = [-1.; 1.]
	@tullio kxE[b,i,j,k] := kxscales[b] * kpg_mag[i,j,k] * ε⁻¹_op(ε⁻¹,Ha,mn)[a,i,j,k] * mn[a,kxinds[b],i,j,k] nograd=(kxscales,kxinds)
end

@btime Mk($Ha,$ε⁻¹,$mn,$kpg_mag)
@tullio out[_] := conj.(Ha)[b,i,j,k] * kxscales[b] * kpg_mag[i,j,k] * ε⁻¹_op(ε⁻¹,Ha,mn)[a,i,j,k] * mn[a,kxinds[b],i,j,k] nograd=(kxscales,kxinds) verbose=true

function H_Mk_H(Ha,ε⁻¹,mn,kpg_mag)
	kxinds = [2; 1]
	kxscales = [-1.; 1.]
	@tullio out[_] := conj.(Ha)[b,i,j,k] * kxscales[b] * kpg_mag[i,j,k] * ε⁻¹_op(ε⁻¹,Ha,mn)[a,i,j,k] * mn[a,kxinds[b],i,j,k] nograd=(kxscales,kxinds)
	return abs(out[1])
end

H_Mk_H(Ha,ε⁻¹,mn,kpg_mag)
ω_mpb / H_Mk_H(Ha,ε⁻¹,mn,kpg_mag)
( ng_mpb - ( ω_mpb / H_Mk_H(Ha,ε⁻¹,mn,kpg_mag) ) ) / ng_mpb
Zygote.gradient(H_Mk_H,Ha,ε⁻¹,mn,kpg_mag)
@btime Zygote.gradient(H_Mk_H,$Ha,$ε⁻¹,$mn,$kpg_mag)
# 3.108 ms (37422 allocations: 7.25 MiB) !!!

# function kpg_mn(kpg,mag)
# 	if mag==0
# 		n = [0.; 1.; 0.]
# 		m = [0.; 0.; 1.]
# 	else
# 		if kpg[1]==0. && kpg[2]==0.    # put n in the y direction if k+G is in z
# 			n = [0.; 1.; 0.] #SVector(0.,1.,0.)
# 		else                                # otherwise, let n = z x (k+G), normalized
# 			ntemp = [0.; 0.; 1.] × kpg
# 			n = ntemp / sqrt(sum2(abs2,ntemp)) # norm(ntemp) #
# 		end
# 	end
# 	# m = n x (k+G), normalized
# 	mtemp = n × kpg
# 	m = mtemp / sqrt(sum2(abs2,mtemp))
# 	return hcat(m,n)
# end

Nx,Ny,Nz = size(ε⁻¹)[end-2:end]
Δz = 1.0
gx = Zygote.@ignore [[ggx;0.0;0.0] for ggx in collect(fftfreq(Nx,Nx/Δx))]
gy = Zygote.@ignore [[0.0;ggy;0.0] for ggy in collect(fftfreq(Ny,Ny/Δy))]
gz = Zygote.@ignore [[0.0;0.0;ggz] for ggz in collect(fftfreq(Nz,Nz/Δz))]
g⃗ = Zygote.@ignore [ [gx;gy;gz] for gx in fftfreq(Nx,Nx/Δx), gy in fftfreq(Ny,Ny/Δy), gz in fftfreq(Nz,Nz/Δz)]
g⃗ₜ_zero_mask = [ sum(abs2.(gg[1:2])) for gg in g⃗ ] .> 0.
g⃗ₜ_zero_mask! = .!(g⃗ₜ_zero_mask)
k⃗ = [0.;0.;kz]
@tullio kpg[a,i,j,k] := k⃗[a] - g⃗[i,j,k][a] nograd=g⃗ verbose=true
@tullio kpg_mag[i,j,k] := sqrt <| kpg[a,i,j,k]^2 verbose=true
zxinds = [2; 1; 3]
zxscales = [-1; 1. ;0.] #[[0. -1. 0.]; [-1. 0. 0.]; [0. 0. 0.]]
ŷ = [0.; 1. ;0.]
# @tullio kpg_nt[a,i,j,k] := zxscales[a] * kpg[zxinds[a],i,j,k] nograd=(zxscales,zxinds) verbose=true
@tullio kpg_nt[a,i,j,k] := zxscales[a] * kpg[zxinds[a],i,j,k] * g⃗ₜ_zero_mask[i,j,k] + ŷ[a] * g⃗ₜ_zero_mask![i,j,k]  nograd=(zxscales,zxinds,ŷ,g⃗ₜ_zero_mask,g⃗ₜ_zero_mask!) verbose=true
@tullio kpg_nmag[i,j,k] := sqrt <| kpg_nt[a,i,j,k]^2  verbose=true
@tullio kpg_n[a,i,j,k] := kpg_nt[a,i,j,k] / kpg_nmag[i,j,k]
xinds1 = [2; 3; 1]
xinds2 = [3; 1; 2]
@tullio kpg_mt[a,i,j,k] := kpg_n[xinds1[a],i,j,k] * kpg[xinds2[a],i,j,k] - kpg[xinds1[a],i,j,k] * kpg_n[xinds2[a],i,j,k] nograd=(xinds1,xinds2) verbose=true
@tullio kpg_mmag[i,j,k] := sqrt <| kpg_mt[a,i,j,k]^2  verbose=true
@tullio kpg_m[a,i,j,k] := kpg_mt[a,i,j,k] / kpg_mmag[i,j,k]
kpg_mn_basis = [[1. 0.] ; [0. 1.]]
@tullio kpg_mn[a,b,i,j,k] := kpg_mn_basis[b,1] * kpg_m[a,i,j,k] + kpg_mn_basis[b,2] * kpg_n[a,i,j,k] nograd=kpg_mn_basis verbose=true

function calc_kpg(kz,Δx,Δy,Δz,Nx,Ny,Nz)
	g⃗ = Zygote.@ignore [ [gx;gy;gz] for gx in fftfreq(Nx,Nx/Δx), gy in fftfreq(Ny,Ny/Δy), gz in fftfreq(Nz,Nz/Δz)]
	g⃗ₜ_zero_mask = Zygote.@ignore [ sum(abs2.(gg[1:2])) for gg in g⃗ ] .> 0.
	g⃗ₜ_zero_mask! = Zygote.@ignore .!(g⃗ₜ_zero_mask)
	ŷ = [0.; 1. ;0.]
	k⃗ = [0.;0.;kz]
	@tullio kpg[a,i,j,k] := k⃗[a] - g⃗[i,j,k][a] nograd=g⃗
	@tullio kpg_mag[i,j,k] := sqrt <| kpg[a,i,j,k]^2
	zxinds = [2; 1; 3]
	zxscales = [-1; 1. ;0.] #[[0. -1. 0.]; [-1. 0. 0.]; [0. 0. 0.]]
	@tullio kpg_nt[a,i,j,k] := zxscales[a] * kpg[zxinds[a],i,j,k] * g⃗ₜ_zero_mask[i,j,k] + ŷ[a] * g⃗ₜ_zero_mask![i,j,k]  nograd=(zxscales,zxinds,ŷ,g⃗ₜ_zero_mask,g⃗ₜ_zero_mask!)
	@tullio kpg_nmag[i,j,k] := sqrt <| kpg_nt[a,i,j,k]^2
	@tullio kpg_n[a,i,j,k] := kpg_nt[a,i,j,k] / kpg_nmag[i,j,k]
	xinds1 = [2; 3; 1]
	xinds2 = [3; 1; 2]
	@tullio kpg_mt[a,i,j,k] := kpg_n[xinds1[a],i,j,k] * kpg[xinds2[a],i,j,k] - kpg[xinds1[a],i,j,k] * kpg_n[xinds2[a],i,j,k] nograd=(xinds1,xinds2)
	@tullio kpg_mmag[i,j,k] := sqrt <| kpg_mt[a,i,j,k]^2
	@tullio kpg_m[a,i,j,k] := kpg_mt[a,i,j,k] / kpg_mmag[i,j,k]
	kpg_mn_basis = [[1. 0.] ; [0. 1.]]
	@tullio kpg_mn[a,b,i,j,k] := kpg_mn_basis[b,1] * kpg_m[a,i,j,k] + kpg_mn_basis[b,2] * kpg_n[a,i,j,k] nograd=kpg_mn_basis
	return kpg_mag, kpg_mn
end

function ng1(ω,ε⁻¹,Δx,Δy,Δz)
	Nx,Ny,Nz = size(ε⁻¹)[end-2:end]
	H,kz = solve_k(ω, ε⁻¹)
	kpg_mag, kpg_mn = calc_kpg(kz,Δx,Δy,Δz,Nx,Ny,Nz)
	Ha = copy(reshape(H,(2,Nx,Ny,Nz)))
	n_g = ω / H_Mk_H(Ha,ε⁻¹,kpg_mn,kpg_mag)
end

ng1(0.65,ε⁻¹,6.0,4.0,1.0)
@btime ng1(0.65,$ε⁻¹,6.0,4.0,1.0)
# 910.910 ms (152467 allocations: 237.08 MiB)
Zygote.gradient(ng1,0.65,ε⁻¹,6.0,4.0,1.0)
@btime Zygote.gradient(ng1,0.65,ε⁻¹,6.0,4.0,1.0)
# 4.695 s (49137386 allocations: 6.82 GiB)

function ng2(H,ω,ε⁻¹,kpg_mn,kpg_mag)
	Ha = copy(reshape(H,(2,Nx,Ny,Nz)))
	n_g = ω / H_Mk_H(Ha,ε⁻¹,kpg_mn,kpg_mag)
end

ng2(H,ω,ε⁻¹,kpg_mn,kpg_mag)
@btime ng2($H,0.65,$ε⁻¹,$kpg_mn,$kpg_mag)
# 1.526 ms (36971 allocations: 2.63 MiB)
Zygote.gradient(ng2,H,0.65,ε⁻¹,kpg_mn,kpg_mag)
@btime Zygote.gradient(ng2,$H,0.65,$ε⁻¹,$kpg_mn,$kpg_mag)
# 3.181 ms (37379 allocations: 7.44 MiB)

function ng3(Ha,ω,ε⁻¹,kpg_mn,kpg_mag)
	n_g = ω / H_Mk_H(Ha,ε⁻¹,kpg_mn,kpg_mag)
end

ng3(Ha,ω,ε⁻¹,kpg_mn,kpg_mag)
@btime ng3($Ha,0.65,$ε⁻¹,$kpg_mn,$kpg_mag)
# 1.596 ms (36962 allocations: 2.45 MiB)
Zygote.gradient(ng3,Ha,0.65,ε⁻¹,kpg_mn,kpg_mag)
@btime Zygote.gradient(ng3,$Ha,0.65,$ε⁻¹,$kpg_mn,$kpg_mag)
# 3.164 ms (37356 allocations: 7.25 MiB)

function ng4(ω,ε⁻¹,kpg_mn,kpg_mag)
	H,kz = solve_k(ω, ε⁻¹)
	Ha = copy(reshape(H,(2,Nx,Ny,Nz)))
	n_g = ω / H_Mk_H(Ha,ε⁻¹,kpg_mn,kpg_mag)
end

ng4(ω,ε⁻¹,kpg_mn,kpg_mag)
@btime ng4(0.65,$ε⁻¹,$kpg_mn,$kpg_mag)
# 971.820 ms (80838 allocations: 245.70 MiB)
Zygote.gradient(ng4,0.65,ε⁻¹,kpg_mn,kpg_mag)
@btime Zygote.gradient(ng4,0.65,$ε⁻¹,$kpg_mn,$kpg_mag)
# 5.081 s (48676432 allocations: 6.76 GiB)

function ng5(kz,H,ω,ε⁻¹,Δx,Δy,Δz)
	Nx,Ny,Nz = size(ε⁻¹)[end-2:end]
	kpg_mag, kpg_mn = calc_kpg(kz,Δx,Δy,Δz,Nx,Ny,Nz)
	Ha = copy(reshape(H,(2,Nx,Ny,Nz)))
	n_g = ω / H_Mk_H(Ha,ε⁻¹,kpg_mn,kpg_mag)
end

ng5(kz,H,ω,ε⁻¹,6.0,4.0,1.0)
@btime ng5($kz,$H,0.65,$ε⁻¹,6.0,4.0,1.0)
# 2.655 ms (55445 allocations: 5.64 MiB)
Zygote.gradient(ng5,kz,H,0.65,ε⁻¹,6.0,4.0,1.0)
@btime Zygote.gradient(ng5,$kz,$H,0.65,$ε⁻¹,6.0,4.0,1.0)
# 6.204 ms (56614 allocations: 14.89 MiB)

@tullio ε_approx[i,j,k] := 3 / ε⁻¹[a,a,i,j,k]
@tullio d_approx[b,i,j,k] := zx(Ha,kpg_mn)[b,i,j,k] * 3 / ε⁻¹[a,a,i,j,k]

function kx_t2c(H,mn,kpg_mag)
	kxscales = [-1.; 1.]
    @tullio d[a,i,j,k] := kxscales[b] * H[b,i,j,k] * mn[a,b,i,j,k] * kpg_mag[i,j,k] nograd=kxscales
end

d⃗ = kx_t2c(Ha,kpg_mn,kpg_mag)
@btime @tullio e⃗[a,i,j,k] := d⃗[b,i,j,k] * ε⁻¹[a,b,i,j,k] # verbose=true
@btime @tullio e⃗[a,i,j,k] := d⃗[b,i,j,k] * ε⁻¹_mpb[i,j,k][a,b] # verbose=true

function ftest(d⃗)
	# @tullio e⃗[a,i,j,k] := d⃗[b,i,j,k] * ε⁻¹_mpb[i,j,k][a,b]
	@tullio e⃗[a,i,j,k] := d⃗[b,i,j,k] * ε⁻¹[a,b,i,j,k]
end

@btime ftest($d⃗)
# 344.180 μs (6154 allocations: 672.33 KiB)  # with StaticArrays ε⁻¹
# 94.073 μs (3 allocations: 288.13 KiB)  # with normal Array ε⁻¹

# function kx_t2c!(ds::MaxwellData,m::Array{Float64,4},n::Array{Float64,4},kpg_mag::Array{Float64,3})::AbstractArray{ComplexF64,4}
function kx_t2c!(ds::MaxwellData,mn::Array{Float64,5},kpg_mag::Array{Float64,3})::AbstractArray{ComplexF64,4}
	# kxscales = [-1.; 1.]
	# @tullio ds.d[a,i,j,k] = kxscales[b] * ds.H[b,i,j,k] * mn[a,b,i,j,k] * kpg_mag[i,j,k] nograd=kxscales
	# @tullio ds.d[a,i,j,k] = ( ds.H[2,i,j,k] * n[a,i,j,k] - ds.H[1,i,j,k] * m[a,i,j,k] ) * kpg_mag[i,j,k]
	@tullio ds.d[a,i,j,k] = ( ds.H[2,i,j,k] * mn[a,2,i,j,k] - ds.H[1,i,j,k] * mn[a,1,i,j,k] ) * kpg_mag[i,j,k]
	# @tullio ds.d[a,i,j,k] = ( ds.H[2,i,j,k] * n[a,i,j,k] - ds.H[1,i,j,k] * m[a,i,j,k] )
	# @tullio ds.d[a,i,j,k] *=  kpg_mag[i,j,k]

    # @inbounds for i=1:ds.Nx,j=1:ds.Ny,k=1:ds.Nz
    #     scale = -ds.kpG[i,j,k].mag #-ds.kpG[i,j,k].mag
    #     ds.d[1,i,j,k] = ( ds.H[1,i,j,k] * ds.kpG[i,j,k].n[1] - ds.H[2,i,j,k] * ds.kpG[i,j,k].m[1] ) * scale
    #     ds.d[2,i,j,k] = ( ds.H[1,i,j,k] * ds.kpG[i,j,k].n[2] - ds.H[2,i,j,k] * ds.kpG[i,j,k].m[2] ) * scale
    #     ds.d[3,i,j,k] = ( ds.H[1,i,j,k] * ds.kpG[i,j,k].n[3] - ds.H[2,i,j,k] * ds.kpG[i,j,k].m[3] ) * scale
    # end
    return ds.d
end

kx_t2c!(ds,mn,kpg_mag)
@btime kx_t2c!($ds,$mn,$kpg_mag)
kx_t2c!(ds,m,n,kpg_mag)
@btime kx_t2c!($ds,$m,$n,$kpg_mag)

##
eigind = 1
P = LinearMap(x -> H[:,eigind] * dot(H[:,eigind],x),length(H[:,eigind]),ishermitian=true)
A = M̂(ε⁻¹,mn,kpg_mag,𝓕) - ω^2 * I
b = ( I  -  P ) * H̄[:,eigind]
λ⃗₀ = IterativeSolvers.bicgstabl(A,b,l)
# ωₖ = real( ( H[:,eigind]' * Mₖ(H[:,eigind], ε⁻¹,kz,gx,gy,gz) )[1]) / ω # ds.ω²ₖ / ( 2 * ω )
Ha = reshape(H,(2,Nx,Ny,Nz))
ωₖ = H_Mₖ_H(Ha,ε⁻¹,mn,kpg_mag) / ω
# Hₖ =  ( I  -  P ) * * ( Mₖ(H[:,eigind], ε⁻¹,ds) / ω )
ω̄  =  ωₖ * k̄
λ⃗₀ -= P*λ⃗₀ - ω̄  * H
Ha_F = 𝓕 * kx_t2c(H,mn,kpg_mag) #fft(kcross_t2c(Ha,kz,gx,gy,gz),(2:4))
λ₀ = reshape(λ⃗₀,(2,Nx,Ny,Nz))
λ₀_F  = 𝓕 * kx_t2c(λ₀,mn,kpg_mag) #fft(kcross_t2c(λ₀,kz,gx,gy,gz),(2:4))
# ε̄ ⁻¹ = ( 𝓕 * kcross_t2c(λ₀,ds) ) .* ( 𝓕 * kcross_t2c(Ha,ds) )
# ε⁻¹_bar = [ Diagonal( real.(λ₀_F[:,i,j,kk] .* Ha_F[:,i,j,kk]) ) for i=1:Nx,j=1:Ny,kk=1:Nz]
ε⁻¹_bar = [ Diagonal( real.(λ₀_F[:,i,j,kk] .* Ha_F[:,i,j,kk]) )[a,b] for a=1:3,b=1:3,i=1:Nx,j=1:Ny,kk=1:Nz]

##

function zyg_zcross_t2c(Hin,kz,gx,gy,gz)
    # Hout = zeros(ComplexF64,(3,Nx,Ny,Nz))
    Nx,Ny,Nz = size(Hin)[2:4]
    Hout = Zygote.Buffer(Hin,3,Nx,Ny,Nz)
	@inbounds for i=1:Nx,j=1:Ny,k=1:Nz
		kpg, mag, m, n = Zygote.forwarddiff(kz) do kz
                            zyg_KplusG(kz,gx[i],gy[j],gz[k])
                        end
        Hout[1,i,j,k] = -Hin[1,i,j,k] * m[2] - Hin[2,i,j,k] * n[2]
        Hout[2,i,j,k] =  Hin[1,i,j,k] * m[1] + Hin[2,i,j,k] * n[1]
        Hout[3,i,j,k] = 0.0
    end
    return copy(Hout)
end

function zyg_kcross_t2c(Hin,kz,gx,gy,gz)
    # Hout = Array{ComplexF64}(undef,(3,Nx,Ny,Nz))
    Nx,Ny,Nz = size(Hin)[2:4]
    Hout = Zygote.Buffer(Hin,3,Nx,Ny,Nz)
	@inbounds for i=1:Nx,j=1:Ny,k=1:Nz
		kpg, mag, m, n = Zygote.forwarddiff(kz) do kz
                            zyg_KplusG(kz,gx[i],gy[j],gz[k])
                        end
        Hout[1,i,j,k] = ( Hin[1,i,j,k] * n[1] - Hin[2,i,j,k] * m[1] ) * -mag
        Hout[2,i,j,k] = ( Hin[1,i,j,k] * n[2] - Hin[2,i,j,k] * m[2] ) * -mag
        Hout[3,i,j,k] = ( Hin[1,i,j,k] * n[3] - Hin[2,i,j,k] * m[3] ) * -mag
    end
    return copy(Hout)
end

function zyg_kcross_c2t(Hin,kz,gx,gy,gz)
    # Hout = Array{ComplexF64}(undef,(2,Nx,Ny,Nz))
    Nx,Ny,Nz = size(Hin)[2:4]
    Hout = Zygote.Buffer(Hin,2,Nx,Ny,Nz)
	@inbounds for i=1:Nx,j=1:Ny,k=1:Nz
		kpg, mag, m, n = Zygote.forwarddiff(kz) do kz
                            zyg_KplusG(kz,gx[i],gy[j],gz[k])
                        end
        at1 = Hin[1,i,j,k] * m[1] + Hin[2,i,j,k] * m[2] + Hin[3,i,j,k] * m[3]
        at2 = Hin[1,i,j,k] * n[1] + Hin[2,i,j,k] * n[2] + Hin[3,i,j,k] * n[3]
        Hout[1,i,j,k] =  -at2 * mag
        Hout[2,i,j,k] =  at1 * mag
    end
    return copy(Hout)
end

function zyg_ε⁻¹_dot(Hin,ε⁻¹)
    # Hout = similar(Hin)
    Nx,Ny,Nz = size(Hin)[2:4]
    Hout = Zygote.Buffer(Hin,3,Nx,Ny,Nz)
    @inbounds for i=1:Nx,j=1:Ny,k=1:Nz
        Hout[1,i,j,k] =  ε⁻¹[1,1,i,j,k]*Hin[1,i,j,k] + ε⁻¹[2,1,i,j,k]*Hin[2,i,j,k] + ε⁻¹[3,1,i,j,k]*Hin[3,i,j,k]
        Hout[2,i,j,k] =  ε⁻¹[1,2,i,j,k]*Hin[1,i,j,k] + ε⁻¹[2,2,i,j,k]*Hin[2,i,j,k] + ε⁻¹[3,2,i,j,k]*Hin[3,i,j,k]
        Hout[3,i,j,k] =  ε⁻¹[1,3,i,j,k]*Hin[1,i,j,k] + ε⁻¹[2,3,i,j,k]*Hin[2,i,j,k] + ε⁻¹[3,3,i,j,k]*Hin[3,i,j,k]
        # Hout[1,i,j,k] =  ε⁻¹[i,j,k][1,1]*Hin[1,i,j,k] + ε⁻¹[i,j,k][1,2]*Hin[2,i,j,k] + ε⁻¹[i,j,k][1,3]*Hin[3,i,j,k]
        # Hout[2,i,j,k] =  ε⁻¹[i,j,k][2,1]*Hin[1,i,j,k] + ε⁻¹[i,j,k][2,2]*Hin[2,i,j,k] + ε⁻¹[i,j,k][2,3]*Hin[3,i,j,k]
        # Hout[3,i,j,k] =  ε⁻¹[i,j,k][3,1]*Hin[1,i,j,k] + ε⁻¹[i,j,k][3,2]*Hin[2,i,j,k] + ε⁻¹[i,j,k][3,3]*Hin[3,i,j,k]
    end
    return copy(Hout)
end


text:
using Documenter

using STMOZOO
using STMOZOO.SudokuSolver

makedocs(sitename="STMO ZOO",
    format = Documenter.HTML(),
    modules=[SudokuSolver], # add your module
    pages=Any[
        "Sudoku"=> "man/sudoku.md",  # add the page to your documentation
    ])

#=
deploydocs(
            repo = "github.com/michielstock/STMOZOO.jl.git",
        )
=#


text:
function tcm_estimateHALOMNL(N, assortments, transaction_counts)
	M = size(assortments, 1)

	total_trans_by_assortment = sum(transaction_counts, 2)

	m_HALOMNL = Model(solver = IpoptSolver(print_level=1))

	@variable(m_HALOMNL, -20 <= u[1:N-1] <= 20)
	@variable(m_HALOMNL, -20 <= alph[1:N-1, 1:N-1] <= 20)

	@NLobjective(m_HALOMNL, Max, sum(sum( transaction_counts[m,p] * (u[p] + sum(alph[q,p]*(1 - assortments[m,q]) for q in 1:N-1 ) ) for p in 1:N-1) for m in 1:M)
							- sum(total_trans_by_assortment[m]*log(1+sum(assortments[m,p]* exp(u[p] + sum(alph[q,p]*(1 - assortments[m,q]) for q in 1:N-1 )) for p in 1:N-1)) for m in 1:M) )

	start_time = time()
	solve(m_HALOMNL)
	end_time = time()
	elapsed_time = end_time - start_time

	u_val = getvalue(u)
	alpha_val = getvalue(alph)

	obj = getobjectivevalue(m_HALOMNL)

	return u_val, alpha_val, obj, elapsed_time
end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
using Test
using FakePkg

@testset "FakePkg" begin
    @test FakePkg.CASTED_COMMANDS["main"].root.name == "pkg"

    @test FakePkg.command_main(["-h"]) == 0
    @test FakePkg.command_main(["--help"]) == 0
    @test FakePkg.command_main(["-V"]) == 0
    @test FakePkg.command_main(["--version"]) == 0

    @test FakePkg.command_main(["add", "-h"]) == 0
    @test FakePkg.command_main(["rm", "-h"]) == 0

    @test FakePkg.command_main(["rm", "CCC", "-h"]) == 0
    @test FakePkg.command_main(["add", "CCC", "-h"]) == 0

    @test FakePkg.command_main(["add", "ABC"]) == 0
    @test FakePkg.command_main(["rm", "ABC"]) == 0
    @test FakePkg.command_main(["activate", "fake", "-s"]) == 0
    @test FakePkg.command_main(["activate", "fake", "--shared"]) == 0

    @test FakePkg.command_main(["registry", "add", "abc"]) == 0
    @test FakePkg.command_main(["registry", "rm", "abc"]) == 0
end


text:
using Documenter
using GeneralAstrodynamics

makedocs(
    format=Documenter.HTML(),
    sitename="GeneralAstrodynamics.jl",
    authors = "Joey Carpinelli",
    pages=[
        "Quick Start" => [
            "Getting Started" => "index.md",
            "Docstrings" => "docstrings.md"
        ]
    ]
)

deploydocs(
    target = "build",
    repo="github.com/cadojo/GeneralAstrodynamics.jl.git",
    branch = "gh-pages",
    deps   = nothing,
    make   = nothing,
    devbranch = "main",
    versions = ["stable" => "v^", "manual", "v#.#", "v#.#.#"],
)


text:
# Copyright (c) 2015 Aki Vehtari, Tuomas Sivula
# Original Matlab version by Aki Vehtari. Translation to Python
# by Tuomas Sivula.

# This software is distributed under the GNU General Public
# License (version 3 or later); please refer to the file
# License.txt, included with the software, for details.

"""
    psisloo(log_lik, [wcpp, wtrunc])

PSIS leave-one-out log predictive densities.

Computes the log predictive densities given posterior samples of the log
likelihood terms p(y_i|\theta^s) in input parameter `log_lik`. Returns a
sum of the leave-one-out log predictive densities `loo`, individual
leave-one-out log predictive density terms `loos` and an estimate of Pareto
tail indeces `ks`. If tail index k>0.5, variance of the raw estimate does
not exist and if tail index k>1 the mean of the raw estimate does not exist
and the PSIS estimate is likely to have large variation and some bias.

# Arguments
* `log_lik::AbstractArray`: Array of size n x m containing n posterior samples of the log likelihood terms p(y_i|\theta^s).
* `wcpp::Real`: Percentage of samples used for GPD fit estimate (default is 20).
* `wtrunc::Float64`: Positive parameter for truncating very large weights to n^wtrunc. Providing False or 0 disables truncation. Default values is 3/4.

# Returns
* `loo::Real`: sum of the leave-one-out log predictive densities.
* `loos::AbstractArray`: Individual leave-one-out log predictive density terms.* `ks::AbstractArray`: Estimated Pareto tail indeces.
"""
# Compute LOO and standard error
function psisloo(log_lik::AbstractArray, wcpp::Int64=20, wtrunc::Float64=3/4)

    # log raw weights from log_lik
    lw = copy(log_lik)

    # compute Pareto smoothed log weights given raw log weights
    lwp, ks = psislw(-lw, wcpp, wtrunc)

    lwp += lw
    loos = reshape(logsumexp(lwp; dims=1), size(lwp, 2))
    loo = sum(loos)

    return loo, loos, ks
end

export
    psisloo


text:
using DelimitedFiles

x = readdlm("data_2.dat", '\t')

global horizontal = 0
global depth = 0
for val in x
    checkval = val[1]
    if checkval == 'f'
        global horizontal += parse(Int, val[9])
    elseif checkval == 'd'
        global depth += parse(Int, val[6])
    elseif checkval == 'u'
        global depth -= parse(Int, val[4])
    end
end

println(horizontal * depth)


text:
println("Started julia processes")
using ParallelTest

# do something in your system to start processes:
# salloc -N 3 --ntasks-per-node=5
wrkers = ParallelTest.machines()

@everywhere using ParallelTest

println("make everybody say hello")

@everywhere ParallelTest.sayhello()

println("make everybody do some math")

pmap( i->ParallelTest.domath(i), [100 for j in 1:length(workers())] )

# println("make everybody pass a memory test")

pmap( i->ParallelTest.doBIGmath(10000,10000), 1:length(workers()) )
#pmap( i->doBIGmath(15000,10000), 1:length(workers()) )

#println("trying parallel for loop with $(nprocs()) processes")
#println("numworkers: $(length(workers()))")
#println("workers: $(workers())")

# println("svd call on master")
t0 = @elapsed ParallelTest.serial();
t0 = @elapsed ParallelTest.serial();
println("serial call takes $t0")
t1 = @elapsed ParallelTest.parallel();
t1 = @elapsed ParallelTest.parallel();
println("parallel call takes $t1")

println(" quitting ")

quit()


text:
#SignalProcessing: Vector operations
#-------------------------------------------------------------------------------

#==Helper functions
===============================================================================#
evenlength(v::AbstractArray) = (0 == (length(v) % 2))

#==Obtain a data vector (TODO: move to MDDatatsets?)
===============================================================================#

datavec(s::Symbol, args...; kwargs...) = datavec(DS{s}(), args...; kwargs...)
datavec(::DS{T}, args...; kwargs...) where {T<:Symbol} = throw(ArgumentError("Unknown vector type: $T"))

#==Range generators
===============================================================================#

#==timespace: Generates time vector (range) from sampling period & fundamental.
Automatically computes # of time steps, N.
Why? Computation of N can easily be off-by-1 - if specified directly...
NOTE:
  -First argument is targeted argument (will minimize amount of rounding).
  -Cannot represent tfund exactly (because maximum(timespace(...)) = tfund-ts).
   You would therefore need to define a new "range" type to represent tfund
   exactly. TODO??==#

timespace(primary::Symbol, v1, secondary::Symbol, v2; kwargs...) =
	timespace(DS{primary}(), DS{secondary}(), v1, v2; kwargs...)
timespace(::DS{T1}, ::DS{T2}, args...; kwargs...) where {T1, T2} =
	throw(ArgumentError("timespace does not support combination (:$T1, :$T2)."))

#Sampling period (timestep), ts is most important parameter:
#Will throw error if "suggested" fundamental period does not land on the timestep:
function timespace(::DS{:ts}, ::DS{:tfund}, ts, tfund; tstart=0)
	ABSTOL = .05 #WANTCONST
	ns = tfund/ts
	ins = round(Int, ns)
	if abs(ns-ins) > ABSTOL
		throw("tfund must be approx. an integer multiple of ts.")
	end
	return tstart .+ StepRangeLen{DataFloat}(0:ts:((ins-1)*ts))
end

#Fundamental period, tfund, is most important parameter.
#Will throw error if fundamental period does not land on the "suggested" timestep:
function timespace(::DS{:tfund}, ::DS{:ts}, tfund, ts; tstart=0)
	ABSTOL = .05 #WANTCONST
	ns = tfund/ts
	ins = round(Int, ns)
	if abs(ns-ins) > ABSTOL
		throw("tfund must be approx. an integer multiple of ts.")
	end
	ts = tfund/ins #Re-compute most accurate possible version of timestep
	return tstart .+ StepRangeLen{DataFloat}(0:ts:((ins-1)*ts))
end

function timetofreq(t::StepRangeLen{DataFloat})
	n = length(t)
	Δt = step(t)
	tfund = Δt*n
	Δf = 1/tfund
	fmax = (floor(n/2)/n)/Δt
	return 0:Δf:fmax
end

function freqtotime(f::StepRangeLen{DataFloat}, teven::Bool)
	tfund = 1/step(f)
	n = 2*length(f) - 1
	if teven
		n -= 1
	end
	Δt = tfund/n
	return timespace(:ts, Δt, :tfund, tfund)
end

#Last line


text:
get_x(vec::SA.SVector) = vec[1]
get_y(vec::SA.SVector) = vec[2]

rotate_plus_90(vec::SA.SVector{2}) = typeof(vec)(-vec[2], vec[1])
rotate_minus_90(vec::SA.SVector{2}) = typeof(vec)(vec[2], -vec[1])
rotate_180(vec::SA.SVector{2}) = -vec

rotate(x::T, y::T, c::T, s::T) where {T} = SA.SVector(c * x - s * y, s * x + c * y)
rotate(vec::SA.SVector{2, T}, theta::T) where {T} = rotate(vec[1], vec[2], cos(theta), sin(theta))
rotate(vec::SA.SVector{2, T}, dir::SA.SVector{2, T}) where {T} = rotate(vec[1], vec[2], dir[1], dir[2])

get_relative_direction(dir1::SA.SVector{2, T}, dir2::SA.SVector{2, T}) where {T} = typeof(dir1)(dir2[1] * dir1[1] + dir2[2] * dir1[2], dir2[2] * dir1[1] - dir2[1] * dir1[2])
invert_relative_direction(dir::SA.SVector{2}) = typeof(dir)(dir[1], -dir[2])

function invert(pos::SA.SVector{2, T}, dir::SA.SVector{2, T}) where {T}
    inv_dir = invert_relative_direction(dir)
    inv_pos = -rotate(pos, inv_dir)
    return inv_pos, inv_dir
end