text:
# Estimate derivatives numerically
# These use centered-differencing for higher accuracy
function derivative_numer{T<:Number}(func::Function, x::T, h::T)
    xp = x + h
    vp = func(xp)
    xm = x - h
    vm = func(xm)
    return (vp-vm)/(xp-xm)
end
function derivative2_numer{T<:Number}(func::Function, x::T, h::T)
    return (func(x + 2h) - 2*func(x) + func(x - 2h)) / 4h^2
end
derivative_numer{T<:Number}(func::Function, x::T, index::Int, h::T) = derivative_numer(func, x, h)  # compatibility with partial derivatives
derivative2_numer{T<:Number}(func::Function, x::T, index::Int, h::T) = derivative2_numer(func, x, h)
function derivative_numer{T<:Number}(func::Function, x::Array{T}, index::Int, h::T)
    xsave = x[index]
    xp = xsave + h
    xm = xsave - h
    x[index] = xp
    vp = func(x)
    x[index] = xm
    vm = func(x)
    x[index] = xsave
    return (vp-vm)/(xp-xm)
end
function derivative2_numer{T<:Number}(func::Function, x::Array{T}, index::(Int,Int,), h::T)
    if index[1] == index[2] # 2nd order derivative in one direction
        v0 = func(x...)
        xsave = x[index[1]]
        xp = xsave + 2h
        xm = xsave - 2h
        x[index[1]] = xp
        vp = func(x...)
        x[index[1]] = xm
        vm = func(x...)
        x[index[1]] = xsave
        deriv = (vp - 2v0 + vm) / 4h^2
    else # cross-derivative, 1st order in each direction
        xpp = [x[1]+h,x[2]+h]
        xpm = [x[1]+h,x[2]-h]
        xmp = [x[1]-h,x[2]+h]
        xmm = [x[1]-h,x[2]-h]
        deriv = (func(xpp...) - func(xpm...) - func(xmp...) + func(xmm...)) / 4h^2
    end
    return deriv
end
function derivative_numer{T<:Number}(func::Function, c::(T...,), index::Int, h::T)
    x = [c...]
    xsave = x[index]
    xp = xsave + h
    xm = xsave - h
    x[index] = xp
    vp = func(x...)
    x[index] = xm
    vm = func(x...)
    x[index] = xsave
    return (vp-vm)/(xp-xm)
end
function derivative2_numer{T<:Number}(func::Function, c::(T...,), index::(Int,Int,), h::T)
    x = [c...]
    return derivative2_numer(func, x, index, h)
end
function derivative_numer{T<:Number}(func::Function, x::T, h::Vector{T})
    d = zeros(T, length(h))
    for i = 1:length(h)
        d[i] = derivative_numer(func, x, h[i])
    end
    return d
end
function derivative2_numer{T<:Number}(func::Function, x, h::Vector{T})
    d = zeros(T, length(h))
    for i = 1:length(h)
        d[i] = derivative2_numer(func, x, h[i])
    end
    return d
end
function derivative_numer{T<:Number}(func::Function, x, index::Int, h::Vector{T})
    d = zeros(T, length(h))
    for i = 1:length(h)
        d[i] = derivative_numer(func, x, index, h[i])
    end
    return d
end
function derivative2_numer{T<:Number}(func::Function, x, index::(Int,Int,), h::Vector{T})
    d = zeros(T, length(h))
    for i = 1:length(h)
        d[i] = derivative2_numer(func, x, index, h[i])
    end
    return d
end
derivative_numer{T<:Number}(func::Function, x::T) = derivative_numer(func, x, (eps(max(abs(x),one(T))))^convert(T, 1/3))
derivative_numer{T<:Number}(func::Function, x::Array{T}, index::Int) = derivative_numer(func, x, index, (eps(max(abs(x[index]),one(T))))^convert(T, 1/3))
derivative_numer{T<:Number}(func::Function, c::(T...,), index::Int) = derivative_numer(func, c, index, (eps(max(abs(c[index]),one(T))))^convert(T, 1/3))
derivative2_numer{T<:Number}(func::Function, x::T) = derivative2_numer(func, x, (eps(max(abs(x),one(T))))^convert(T,1/3))
derivative2_numer{T<:Number}(func::Function, x::Array{T}, index::(Int,Int,)) = derivative2_numer(func, x, index, (eps(max(abs(x[index[1]]),one(T))))^convert(T, 1/3))
derivative2_numer{T<:Number}(func::Function, c::(T...,), index::(Int,Int,)) = derivative2_numer(func, c, index, (eps(max(abs(c[index[1]]),one(T))))^convert(T, 1/3))


text:
"""
    Agent(id, opinion, inclin_interact, check_regularity, desired_input_count)

Provide data structure for agents within agent-based simulation.

# Examples
```julia-repl
julia>using ABM4OSN

julia>ABM4OSN.Agent(0, 0, 0, 0, 0)
ABM4OSN.Agent(0, 0.0, 0.0, 0.0, 0.0, 0, true, 0, ABM4OSN.Post[], ABM4OSN.Post[], ABM4OSN.Post[], ABM4OSN.Post[])
```

# Arguments
- `id`: Agent identifier within simulation
- `opinion`: Opinion value between -1 and 1
- `inclin_interact`: Value for inclination to interact between 0 and 1
- `check_regularity`: Value for check regularity between 0 and 1
- `desired_input_count`: How many other agents an agent ideally wants to follow

See also: [Post](@ref), [generate_opinion](@ref), [generate_inclin_interact](@ref), [generate_check_regularity](@ref), [generate_desired_input_count](@ref)
"""
mutable struct Agent

    id::Int64
    opinion::Float64
    inclin_interact::Float64
    perceiv_publ_opinion::Float64
    check_regularity::Float64
    desired_input_count::Int64
    active::Bool
    inactive_ticks::Int16
    feed::Array{Post, 1}
    liked_posts::Array{Post, 1}
    disliked_posts::Array{Post, 1}
    shared_posts::Array{Post, 1}

    function Agent(
        id,
        opinion,
        inclin_interact,
        check_regularity,
        desired_input_count
    )

        if opinion < -1 || opinion > 1
            error("invalid opinion value")
        end
        if inclin_interact < 0
            error("invalid value for inclination to interact")
        end

        new(
            id,
            opinion,
            inclin_interact,
            opinion,
            check_regularity,
            desired_input_count,
            true,
            0,
            Array{Post, 1}(undef, 0),
            Array{Post, 1}(undef, 0),
            Array{Post, 1}(undef, 0),
            Array{Post, 1}(undef, 0)
        )

    end

end

# suppress output of include()
;

ext: jl


text:
module TestD
    @show "???!"
    d = 10

end


text:
module MT1DInversion
using ..AbstractOperator, ..CommonToAll
import ..AbstractOperator.get_misfit
import ..Model, ..Options
using ..MT1D
using Random, PyPlot

abstract type MT1DZ <: Operator1D
end

mutable struct MT1DZ_nodepthprior <: MT1DZ
    d_log10_ρ
    d_phase_deg
    σ_log10_ρ
    σ_phase_deg
    freqs
    zboundaries
    irxlayer
end

function MT1DZ_nodepthprior(;
                            d_log10_ρ = nothing,
                            d_phase_deg = nothing,
                            σ_log10_ρ = nothing,
                            σ_phase_deg = nothing,
                            freqs = nothing,
                            zboundaries = nothing,
                            irxlayer = nothing)
    @assert !isa(d_log10_ρ, Nothing)
    @assert !isa(d_phase_deg, Nothing)
    @assert !isa(σ_log10_ρ, Nothing)
    @assert !isa(σ_phase_deg, Nothing)
    @assert !isa(freqs, Nothing)
    @assert !isa(zboundaries, Nothing)
    @assert !isa(irxlayer, Nothing)
    @assert all(diff(zboundaries) .> 0)
    MT1DZ_nodepthprior(d_log10_ρ, d_phase_deg, σ_log10_ρ, σ_phase_deg, freqs, zboundaries, irxlayer)
end

mutable struct MT1DZ_depthprior <: MT1DZ
    d_log10_ρ
    d_phase_deg
    σ_log10_ρ
    σ_phase_deg
    freqs
    zboundaries
    irxlayer
    stretch
    stretchmodel
    low
    Δ
end

function makestretchop(F::MT1DZ_nodepthprior; ρlow=nothing, Δ=nothing)
    @assert !isa(ρlow, Nothing)
    @assert !isa(Δ, Nothing)
    @assert all(Δ .> 0)
    A = []
    for fn in fieldnames(typeof(F))
        push!(A, getfield(F, fn))
    end
    F = MT1DZ_depthprior(A..., true, copy(ρlow), ρlow, Δ)
end

function get_misfit(m::Model, opt::Options, F::MT1DZ)
    if stretchexists(F)
        F.stretchmodel[:] .= 10 .^(F.low + vec(m.fstar).*F.Δ)
        get_misfit(F.stretchmodel, opt, F)
    else
        get_misfit(10 .^m.fstar, opt, F)
    end
end
# above defined function and type signature MUST be defined

function get_misfit(ρ::AbstractArray, opt::Options, F::MT1DZ)
    if !opt.debug
        get_misfit(F.d_log10_ρ, F.d_phase_deg, F.σ_log10_ρ, F.σ_phase_deg, F.freqs, ρ, F.zboundaries, F.irxlayer)
    else
        0.0
    end
end

function get_misfit(d_log10_ρ, d_phase_deg, σ_log10_ρ, σ_phase_deg, freqs, ρ, z, irxlayer=1)
    Z = get_Z(freqs, ρ, z, irxlayer)
    ρₐ = MT1D.ρapp(freqs, Z)
    ϕ  = MT1D.phase(Z)
    r = [(d_log10_ρ - log10.(ρₐ))./σ_log10_ρ ; (d_phase_deg - ϕ)./σ_phase_deg]
    r'r/2
end

get_Z(freqs, ρ, z, irxlayer) = MT1D.Z_f(freqs, ρ, diff(z), irxlayer)

function add_noise(ρ, z, freqs; noisefrac=0.05, rseed=1, irxlayer=1)
    Random.seed!(rseed)
    Z = get_Z(freqs, ρ, z, irxlayer)
    σ_real = noisefrac*abs.(Z)/2
    noise = σ_real.*( randn(length(freqs)) + 1im*randn(length(freqs)) )
    # same as sqrt(2)*σ_real.*randn(ComplexF64, length(freqs))
    Z + noise
end

function create_synthetic(;ρ=nothing, zboundaries=nothing, freqs=nothing,
                        noisefrac=0.05, rseed=1, showplot=true,
                        logscaledepth=true, showfreq=false, gridalpha=0.5, irxlayer=1)
    @assert !isa(ρ, Nothing)
    @assert !isa(zboundaries, Nothing)
    @assert !isa(freqs, Nothing)
    Znoisy = add_noise(ρ, zboundaries, freqs, noisefrac=noisefrac, rseed=rseed, irxlayer=irxlayer)
    ρₐnoisy = MT1D.ρapp(freqs, Znoisy)
    d_log10_ρ = log10.(ρₐnoisy)
    d_phase_deg = MT1D.phase(Znoisy)
    σ_log10_ρ = noisefrac/log(10)
    σ_phase_deg = rad2deg(noisefrac/2)
    F = MT1DZ_nodepthprior(d_log10_ρ, d_phase_deg, σ_log10_ρ, σ_phase_deg, freqs, zboundaries, irxlayer)
    if showplot
        fig = MT1D.plotmodelcurve(1 ./freqs, ρ, zboundaries, logscaledepth=logscaledepth, showfreq=showfreq, irxlayer=irxlayer)
        plotdata(F, fig, iaxis=2, showfreq=showfreq, gridalpha=gridalpha)
    end
    F
end

function plotdata(F::MT1DZ; showfreq=false, gridalpha=0.5, figsize=(6,3))
    fig, ax = plt.subplots(1, 2, sharex=true, figsize=figsize)
    plotdata(F, fig, showfreq=showfreq, gridalpha=gridalpha)
    fig
end

function plotdata(F::MT1DZ, fig; iaxis=1, showfreq=false, gridalpha=0.5)
    xlabel, abcissa = MT1D.f_or_T(F.freqs, showfreq=showfreq)
    ax = fig.axes
    ax[iaxis].errorbar(abcissa, F.d_log10_ρ, F.σ_log10_ρ, linestyle="none", marker=".", elinewidth=0, capsize=3)
    ax[iaxis+1].errorbar(abcissa, F.d_phase_deg, F.σ_phase_deg, linestyle="none", marker=".", elinewidth=0, capsize=3)
    MT1D.labelaxis(xlabel, ax, iaxis, gridalpha=gridalpha)
    ax[iaxis].set_xscale("log")
    fig.tight_layout()
end

function plot_posterior(F::MT1DZ, M::AbstractArray; showfreq=false, gridalpha=0.5, logscaledepth=true,
                            figsize=(10,4), lcolor="nocolor", modelalpha=0.5)
    fig = figure(figsize=(figsize))
    s1 = subplot(131)
    s2 = subplot(132)
    s3 = subplot(133, sharex=s2)
    mnew = ones(length(M[1]))
    for m in M
        if stretchexists(F)
            mnew[:] = 10 .^(F.low + vec(m).*F.Δ)
        else
            mnew[:] .= 10 .^m[:]
        end
        MT1D.plotmodelcurve(1 ./F.freqs, mnew, F.zboundaries, fig, showfreq=showfreq, irxlayer=F.irxlayer,
                        gridalpha=gridalpha, logscaledepth=logscaledepth, lcolor=lcolor, modelalpha=modelalpha)
    end
    plotdata(F, fig, iaxis=2, showfreq=showfreq, gridalpha=gridalpha)
end

function plotpriorenv(F::MT1DZ_depthprior; ax=nothing, lw=2, lc="r", plotlinear=true)
    if isa(ax, Nothing)
        fig = figure()
        ax = gca()
        plotlinear ? ax.set_xlabel("ohm-m") : ax.set_xlabel("Log₁₀ ohm-m")
        ax.set_ylabel("Depth m")
        ax.invert_yaxis()
        fig.suptitle("Prior bounds")
    end
    zlast = F.zboundaries[end] + diff(F.zboundaries)[end]
    low, high = copy(F.low), F.low + F.Δ
    plotlinear && (low = 10 .^low; high = 10 .^high)
    ax.step([low; low[end]], [F.zboundaries; zlast], linewidth=lw, color=lc)
    ax.step([high; high[end]], [F.zboundaries; zlast], linewidth=lw, color=lc)
    plt.tight_layout()
end

end

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text: println("Wow, such empty")


text:
function test_linalg(AT)
    @testset "linear algebra" begin
        @testset "adjoint and transpose" begin
            @test compare(adjoint, AT, rand(Float32, 32, 32))
            @test compare(adjoint!, AT, rand(Float32, 32, 32), rand(Float32, 32, 32))
            @test compare(transpose, AT, rand(Float32, 32, 32))
            @test compare(transpose!, AT, rand(Float32, 32, 32), rand(Float32, 32, 32))
            @test compare((x,y)->copyto!(x, adjoint(y)), AT, rand(Float32, 32, 32), rand(Float32, 32, 32))
            @test compare((x,y)->copyto!(x, transpose(y)), AT, rand(Float32, 32, 32), rand(Float32, 32, 32))
            @test compare(transpose!, AT, Array{Float32}(undef, 32, 32), rand(Float32, 32, 32))
            @test compare(transpose!, AT, Array{Float32}(undef, 128, 32), rand(Float32, 32, 128))
        end

        @testset "copyto! for triangular" begin
            ga = Array{Float32}(undef, 128, 128)
            gb = AT{Float32}(undef, 128, 128)
            rand!(gb)
            copyto!(ga, UpperTriangular(gb))
            @test ga == Array(collect(UpperTriangular(gb)))
            ga = Array{Float32}(undef, 128, 128)
            gb = AT{Float32}(undef, 128, 128)
            rand!(gb)
            copyto!(ga, LowerTriangular(gb))
            @test ga == Array(collect(LowerTriangular(gb)))
        end

        @testset "permutedims" begin
            @test compare(x -> permutedims(x, (2, 1)), AT, rand(Float32, 2, 3))
            @test compare(x -> permutedims(x, (2, 1, 3)), AT, rand(Float32, 4, 5, 6))
            @test compare(x -> permutedims(x, (3, 1, 2)), AT, rand(Float32, 4, 5, 6))
            @test compare(x -> permutedims(x, [2,1,4,3]), AT, randn(ComplexF64,3,4,5,1))
        end

        @testset "issymmetric/ishermitian" begin
            n = 128
            areal = randn(n,n)/2
            aimg  = randn(n,n)/2

            @testset for eltya in (Float32, Float64, ComplexF32, ComplexF64)
                a = convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal)
                asym = transpose(a) + a        # symmetric indefinite
                aherm = a' + a                 # Hermitian indefinite
                @test issymmetric(asym)
                @test ishermitian(aherm)
            end
        end
        @testset "Array + Diagonal" begin
            n = 128
            A = AT{Float32}(undef, n, n)
            d = AT{Float32}(undef, n)
            rand!(A)
            rand!(d)
            D = Diagonal(d)
            B = A + D
            @test collect(B) ≈ collect(A) + collect(D)
        end

        @testset "$f! with diagonal $d" for (f, f!) in ((triu, triu!), (tril, tril!)),
                                            d in -2:2
            A = randn(10, 10)
            @test f(A, d) == Array(f!(AT(A), d))
        end

        @testset "matrix multiplication" begin
            for (a,b) in [((3,4),(4,3)), ((3,), (1,3)), ((1,3), (3))], T in supported_eltypes()
                @test compare(*, AT, rand(T, a), rand(T, b))

                if length(a) > 1
                    @test compare(*, AT, transpose(rand(T, reverse(a))), rand(T, b))
                    @test compare(*, AT, adjoint(rand(T, reverse(a))), rand(T, b))
                end

                if length(b) > 1
                    @test compare(*, AT, rand(T, a), transpose(rand(T, reverse(b))))
                    @test compare(*, AT, rand(T, a), adjoint(rand(T, reverse(b))))
                end

                if length(a) > 1 && length(b) > 1
                    @test compare(*, AT, transpose(rand(T, reverse(a))), transpose(rand(T, reverse(b))))
                    @test compare(*, AT, adjoint(rand(T, reverse(a))), adjoint(rand(T, reverse(b))))
                end

                @test compare(rmul!, AT, rand(T, a), Ref(rand(T)))
                @test compare(lmul!, AT, Ref(rand(T)), rand(T, b))
            end
        end
    end
end


text:
using EMPCA
using Test

@testset "will EMPCA.jl work at all?" begin
    nobs = 100
    nvar = 200
    nvec = 3
    data = zeros(nobs, nvar)

    #- Generate data
    x = (0:(nvar - 1)) ./ (nvar - 1) .* (2 * π)
    for i in 1:nobs
        for k in 1:nvec
            data[i, :] += 5.0*nvec//(k)^2 * randn() * sin.(x*(k))
        end
    end

    #- Add noise
    sigma = ones(size(data))
    for i in Int.(1:(nobs//10))
        sigma[i, :] *= 5
        sigma[i, 1:(Int(floor(nvar/4)) + 1)] *= 5
    end

    weights = 1.0 ./ sigma.^2
    noisy_data = data + (randn(size(sigma)) .* sigma)
    m0 = empca.empca(noisy_data, weights, niter=20)
    @test isreal(m0.R2())
    m1 = empca.lower_rank(noisy_data, weights, niter=20)
    @test isreal(m1.R2())
    m2 = empca.classic_pca(noisy_data)
    @test isreal(m2.R2())
end


text:
#
# Co-occurence data structures
#
CooccurenceDict{T<:AbstractFloat} = DataStructures.DefaultDict{NTuple{2, Int}, T, T}
CooccurenceVector{T<:AbstractFloat} = Vector{Tuple{Int, Int, T}}

# TODO: better name
# TODO: not sure how the coccurences thing works here yet
"""
weightedsums counts the co-occurences as specified
in the GloVe paper.

Returns a `CooccurenceDict` with the weighted co-occurence sum.
"""
function weightedsums(::Type{T}, table::LookupTable, tokens::Vector;
  window::Int=10) where T<:AbstractFloat

    codict = CooccurenceDict{T}(T(0.0))

    for i = 1:length(tokens)
      @inbounds mtok = tokens[i]
      # skip tokens not in the table
      if !haskey(table, mtok)
        continue
      end
      mtok_id = table[mtok]

      for j = max(1, i - window):i-1
        @inbounds ctok = tokens[j]
        if !haskey(table, ctok)
          continue
        end
        ctok_id = table[ctok]

        # The farther away the context token is from the
        # main token the less it contributes
        dist = i - j
        incr = T(1.0 / dist)

        p1 = tuple(mtok_id, ctok_id)
        p2 = tuple(ctok_id, mtok_id)

        # Symmetric context
        codict[p1] += incr
        codict[p2] += incr
      end
    end
    codict
end

# Reading from files is kinda of slow still. Recommending to just use readall
# TODO: make this faster
# function make_cooccur(table::LookupTable, filepath::AbstractString; window::Int=10)
#     codict = DataStructures.DefaultDict(NTuple{2,Int}, Float64, 0.0)
#     open(filepath) do f
#         for line = eachline(f)
#             tokens = split(line)
#             for i = 1:length(tokens)
#                 @inbounds mtok = tokens[i]
#                 mtok_id = table[mtok]
#
#                 for j = max(1, i - window):i-1
#                     @inbounds ctok = tokens[j]
#                     ctok_id = table[ctok]
#
#                     # The farther away the context token is from the
#                     # main token the less it contributes
#                     dist = i - j
#                     incr = 1.0 / dist
#
#                     p1 = tuple(mtok_id, ctok_id)
#                     p2 = tuple(ctok_id, mtok_id)
#
#                     # Symmetric context
#                     codict[p1] += incr
#                     codict[p2] += incr
#                 end
#             end
#         end
#     end
#     codict
# end


text:
# simple print statement
println("Hello World")

# julia REPL can be used as a calculator
4+11

# simple function to print hello world
function helloworld()
    println("Hello world from chapter-1.jl")
end
helloworld()

# include hello.jl
include("hello.jl")

# helloworld from hello.jl
helloworld()

max_stars_repo_name: PacktPublishing/Learning-Julia

max_issues_repo_name: PacktPublishing/Learning-Julia

max_forks_repo_name: PacktPublishing/Learning-Julia


text:
#-----------------------------------------------#
#					dirdiff.jl					#
#-----------------------------------------------#
# Command line execution script for analysing	#
# the differences in directory contents, used	#
# mainly to compare SD card versions which		#
# contain `.bmp` and config files.				#
#						~						#
#				Pierre Steen 2021				#
#===============================================#

## Packages

using CSV
using FileIO
using DataFrames
using ImageMagick
using Images
using ImageIO
using ArgParse
using OutputCollectors

## Custom Type Definitions

"""
	Differences{T<:AbstractDict}

Type which contains
"""
struct Differences{T<:AbstractDict}
	unchanged::T
	reassigned::T
	renamed::T
	removed::T
	added::T
end

## `Differences` Type Constructor

"""
	calcdifferences(collection1, collection2)

Calculates and constructs a `::Differences` datatype instance for the input
"""
function calcdifferences(collection1, collection2)::Differences
	return Differences(
		unchangedpairs(collection1, collection2),
		reassignedpairs(collection1, collection2),
		renamedpairs(collection1, collection2),
		diffpairs(collection1, collection2),	# removed (v51 to v54)
		diffpairs(collection2, collection1)		# added
	)
end

## ArgParse Setup

"""
	parse_commandline()

Parses and returns a dictionary of parsed command line arguments and flags.
"""
function parsecommandline()
	s = ArgParseSettings()

	@add_arg_table! s begin
		# "--bmp", "-b"
			# help = "an option with an argument"
		# "--flag1"
			# help = "an option without argument, i.e. a flag"
			# action = :store_true
		"dir1"
		# dir1 represents the old directory version passed in as a string
			help = "references directory - old version"
			arg_type = String
			required = true
		"dir2"
		# dir2 represents the new directory version passed in as a string
			help = "target directory - new version"
			arg_type = String
			required = true
	end

	return parse_args(s)
end

## Analysis

"""
	encodeimages(collection::Dict{String, Matrix{ColorTypes.RGB{FixedPointNumbers.N0f8}}})

Encode a dictionary of `filename::String => image::Matrix{ColorTypes.RGB{FixedPointNumbers.N0f8}}` to a dictionary of `id::Int => image::Matrix{ColorTypes.RGB{FixedPointNumbers.N0f8}}` pairs.
Ensure image entries are unique before adding to dictionary.
"""
function encodeimages(collection::T) where {T<:AbstractDict}
	# encoded = Dict{Int, Matrix{ColorTypes.RGB{FixedPointNumbers.N0f8}}}() # initialise empty dictionary
	encoded = Dict()
	count = 1
	for pair in collection
		# ensure image not already encoded for (such that allunique(values(encoded)) = true)
		if !in(pair[2], values(encoded))
			push!(encoded, count => pair[2]) # add unique image to encode dictionary
			count += 1 # only increment count if image has been added
		end
	end

	return encoded
end

"""
	unchangedpairs(collection1::T, collection2::T) where {T<:AbstractDict}

Calculates which pairs are unchanged when moving comparing `collection2` to `collection1`, where both collections are dictionaries of `(file_name => file_content)` pairs.
"""
function unchangedpairs(collection1::T, collection2::T) where {T<:AbstractDict}
	# unchanged = Dict{String, Matrix{ColorTypes.RGB{FixedPointNumbers.N0f8}}}()
	unchanged = Dict()
	for pair in collection2
		if in(pair, collection1)
			push!(unchanged, pair)
		end
	end

	return unchanged
end

"""
	reassignedpairs(collection1::T, collection2::T) where {T<:AbstractDict}
"""
function reassignedpairs(collection1::T, collection2::T) where {T<:AbstractDict}
	# reassigned = Dict{String, Matrix{ColorTypes.RGB{FixedPointNumbers.N0f8}}}()
	reassigned = Dict()
	for pair in collection2
		if in(pair[1], keys(collection1)) && !in(pair[2], values(collection1))
			push!(reassigned, pair)
		end
	end

	return reassigned
end

"""
	renamedpairs(collection1::T, collection2::T) where {T<:AbstractDict}
"""
function renamedpairs(collection1::T, collection2::T) where {T<:AbstractDict}
	# renamed = Dict{String, Matrix{ColorTypes.RGB{FixedPointNumbers.N0f8}}}()
	renamed = Dict()
	for pair in collection2
		if in(pair[2], values(collection1)) && !in(pair[1], keys(collection1))
			push!(renamed, pair)
		end
	end

	return renamed
end

"""
	diffpairs(collection1::T, collection2::T) where {T<:AbstractDict}

Works for both removed and added pairs, switch argument order to change functionality.
"""
function diffpairs(collection1::T, collection2::T) where {T<:AbstractDict}
	# removed = Dict{String, Matrix{ColorTypes.RGB{FixedPointNumbers.N0f8}}}()
	removed = Dict()
	for pair in collection1
		if !in(pair[2], values(collection2)) && !in(pair, collection2) && !in(pair[1], keys(collection2))
			push!(removed, pair)
		end
	end

	return removed
end

## Shell diff and `.csv` file build

"""
	argstodirs(parsed_args)

Assembles the output of `OutputCollector -> parsedcommandline()` to an array of directory strings.
"""
function argstodirs(parsed_args)
	return (parsed_args["dir1"], parsed_args["dir2"])
end

"""
	streamfilenames(directories)

Navigate to ref and target directory then run shell command:
```
find . -type f
```
which will create and output stream of all the file names in the current directory.
Do this for each directory we want to analyse, use OutputCollector to collect stdout
and stderr IO streams.
Resulting output string will be formatted somewhat like
```
output = "<filename>\n<filename>\n<filename>\n..."
```
we can therefore use `x -> split(x, "\n")` to split the large string into sub-string filenames, removing the newline chars.
Finally we return an array containing two sub-string vectors.
"""
function streamfilenames(directories)
	returnto = @__DIR__
	filenames = []
	script = """
	#!/bin/bash
	find . -type f
	"""

	for dir in directories
		cd(dir)
		oc = OutputCollector(`bash -c $script`; verbose=false)
		merge(oc)
		if isempty(collect_stderr(oc))
			push!(filenames, String.(split(collect_stdout(oc), "\n")[1:(end-1)]))
		else
			break
		end
	end
	cd(returnto)

	return (filenames[1], filenames[2])
end

# refactored file name shell command stream collector: ------------!!!
# function streamfilenames(directory)
	# locald = @__DIR__
	# script = """
	# #!/bin/bash
	# find . -type f
	# """
#
	# cd(directory)
	# oc = OutputCollector(`bash -c $script`; verbose=false)
	# @assert isempty(collect_stderr(oc))
#
	# names = String.(split(collect_stdout(oc)[1:(end-1)], "\n"))
	# cd(locald)
#
	# return names
# end

"""
REFACTOR TO: filterformat(mixed_fileformats; format="bmp") --------!!!

	filterformat(mixed_filenames)

Split `.txt` and `.bmp` files from array of mixed string file names.
"""
function filterformat(mixed_filenames)
	bmp_filenames = String[]
	txt_filenames = String[]
	for filename in mixed_filenames
		if split(filename, ".")[3] != "txt"
			push!(bmp_filenames, filename)
		else
			push!(txt_filenames, filename)
		end
	end

	return bmp_filenames, txt_filenames
end

"""
	formatdifferences(analysis::Differences)

Formats the directory file content difference analysis results into a form which can later be parsed and used for programmatically determining the necessary file downloads, name changes and deletions to perform to convert the group of files from the reference directory to the target directory.
"""
# function formatdifferences(images::Tuple{T,T}, text::Tuple{T,T}, analysis::Differences) where {T<:AbstractDict}
function formatdifferences(image_analysis::Differences, text_analysis::Differences)
	# write a file (csv?, txt? JSON?)
	# data that needs to be conveyed:
	# reference file:	<filename>.<extension>
	# diff analysis:	<unchanged/reassigned/renamed/removed/added>
	#
end

## Tests

function testdiff(diff)
	for pair in diff.unchanged
		@assert !in(pair, diff.reassigned)
		@assert !in(pair, diff.renamed)
		@assert !in(pair, diff.removed)
		@assert !in(pair, diff.added)
	end
	for pair in diff.reassigned
		@assert !in(pair, diff.unchanged)
		@assert !in(pair, diff.renamed)
		@assert !in(pair, diff.removed)
		@assert !in(pair, diff.added) # fails for v51, v54 diff
	end
	for pair in diff.renamed
		@assert !in(pair, diff.unchanged)
		@assert !in(pair, diff.reassigned)
		@assert !in(pair, diff.removed)
		@assert !in(pair, diff.added)
	end
	for pair in diff.removed
		@assert !in(pair, diff.unchanged)
		@assert !in(pair, diff.reassigned)
		@assert !in(pair, diff.renamed)
		@assert !in(pair, diff.added)
	end
	for pair in diff.added
		@assert !in(pair, diff.unchanged)
		@assert !in(pair, diff.reassigned) # fails for v51, v54 diff -- overlap between added & reassigned sets
		@assert !in(pair, diff.renamed)
		@assert !in(pair, diff.removed)
	end
end

## Main Function

# Next steps:
# - [x] split .bmp txt file
# - [x] construct Pair(filename => image) and Pair(filename => content) dictionaries
# - [x] calculate changes and construct Differences instance
# - []
# - []
function main()
	# read args then capture filenames from bash IO stream using OutputCollector
	# parsed_args = parsecommandline()

	# use when testing from julia> REPL
	parsed_args = Dict(
		"dir1" => "/home/pb/Documents/work/dev-ops/sdcards/CHGCTV51/",
		"dir2" => "/home/pb/Documents/work/dev-ops/sdcards/CHGCTV54/"
	)

	# generate directory strings, then make filename string vector tuple
	dir_paths = argstodirs(parsed_args)
	ref_names, tgt_names = streamfilenames(dir_paths)

	ref_bmp, ref_txt = filterformat(ref_names) # filter reference and target filenames lists into `.bmp` and `.txt`
	tgt_bmp, tgt_txt = filterformat(tgt_names)

	# load images and text files into dictionaries, keeping track of (filename => contents) pairs

	# images
	ref_imgs = Dict(ref_bmp .=> load.(dir_paths[1] .* ref_bmp))
	tgt_imgs = Dict(tgt_bmp .=> load.(dir_paths[2] .* tgt_bmp))

	# txt
	ref_text = Dict(ref_txt .=> read.(dir_paths[1] .* ref_txt, String))
	tgt_text = Dict(tgt_txt .=> read.(dir_paths[2] .* tgt_txt, String))

	# calculate differences between reference and target files, construct Differences data type
	imgs_diff = calcdifferences(ref_imgs, tgt_imgs) # ::Differences
	text_diff = calcdifferences(ref_text, tgt_text)
end

# @show main() # program execution
@show out = main()

max_stars_repo_name: zceepst/dirdiff

max_issues_repo_name: zceepst/dirdiff

max_forks_repo_name: zceepst/dirdiff


text:
# shignal_dft_shifted.jl
##########################
fs=Int(1000);
Ts=1/fs;
N=Int(12000);
T=N/fs;
#t=(0:Int(N-1))*Ts;

x=zeros(Float64, N)
for m=1:N
    x[m]=sin(2*Ts*(m-1)*(Ts*(m-1)-3)*(Ts*(m-1)-6)*(Ts*(m-1)-9)*(Ts*(m-1)-12))
end

using FFTW, DSP

hatx=fft(x);
k = (0:length(hatx)-1)*50/length(hatx);
G=broadcast(abs, hatx);
kshift = (-length(x)/2:length(x)/2-1)*(50/length(x));
Gshift = fftshift(G);

using Plots

plot(kshift, Gshift,
    xaxis="k = -N/2, ..., 0, ..., N/2",
    yaxis="|F[x](k)|",
    title="DFT of x(t)=sin(2t(t-3)(t-6)(t-9)(t-12))",
    grid=false,
    legend=false)
plot!(xticks = ([-25 0 25;], [-6000 0 6000]))
plot!(yticks = ([0 200 400;], [0 200 400]))
savefig("signal_dft_shifted.png")


text:
# this script removes duplicate sections from log files, in favor of the lower
#      one (in the file)
# ARGS[1] = file name
# ARGS[2] = column of file to inspect
# ARGS[3] = output file name (optional), if not specified, "_cleaned" is
#           appended to the file name (before the extension)

using ODLCommonTools

function clean_log(fname::AbstractString, col::Int, output_fname::AbstractString)

  print("loading data...")
  data = readdlm(fname)
  println("done")

  @assert eltype(data) <: Number  # no arrays of strings allowed

  if size(data, 1) == 1  # handle special case
    writedlm(output_fname, data)
    return nothing
  end

  rnges = find_src_ranges(data, col)
  dest_rnges = find_dest_ranges(data, rnges, col)

  println("checking overlap")
  check_overlap(data, rnges, dest_rnges, col)

  print("copying data...")
  data2 = copy_blocks(data, rnges, dest_rnges, col)
  println("done")

  print("writing data...")
  writedlm(output_fname, data2)
  println("done")

  return nothing
end

"""
  Returns an array of Range objects that describe the monotonic blocks.
  Monotonicity is measured using the data in column col.

  **Inputs**

   * data: the array of data
   * col: the index of the column

  **Outputs**

   * rnges: array of Range objects containing the indices of each monotonic
            block
"""
function find_src_ranges(data::AbstractArray, col::Integer)

  rnges = Array{UnitRange{Int}}(0)

  start_idx = 1
  for i=2:size(data, 1)
    if data[i, col] <= data[i-1, col]
      rng = start_idx:(i-1)
      push!(rnges, rng)
      start_idx = i
    end
  end

  # add the last range
  rng = start_idx:size(data, 1)
  push!(rnges, rng)

  return rnges
end

"""
  This function figures out the ranges on the destination array corresponding
  to each block in the source array.  Non-used blocks have the range -1:-1.
  Not all unused blocks are guaranteed to have this sentinel value

  **Inputs**

   * data: the data array
   * rnges: source ranges, each must have montonically increasing indices in the
            data arrange
  * col: column of data array to look at

  **Outputs*

   * dest_rnges: ranges in the destination array
"""
function find_dest_ranges(data::AbstractArray, rnges::Array{UnitRange{Int}}, col::Integer)

  # figure out the destination ranges for each block
  # this assumes the blocks will be copies first to last, overwriting previous
  # data

  nblocks = length(rnges)
  dest_rnges = Array{UnitRange{Int}}(nblocks)
  dest_rng_idx = 1

  # first range starts at 1
  dest_rng = 1:length(rnges[1])
  dest_rnges[dest_rng_idx] = dest_rng
  dest_rng_idx += 1

  # find where range i should start overwriting range i - 1
  block_offset = 0  # global offset where block - 1 starts
  block_offsets = zeros(Int, nblocks)  # store all of them, so we can
                                       # search older blocks
  blocks_used = Array{Bool}(nblocks)
  fill!(blocks_used, true)
  for i=2:length(rnges)
    rng = rnges[i]
    block_offsets[i-1] = block_offset

    # search for first index in next block in the current block
    first_val = data[ first(rng), col]

    # search older blocks until match found
    found_block = 0
    found_idx = 0
    for block=(i-1):-1:1
      vals = sub(data, rnges[block], col)

      # TODO: this could be faster with a log(N) search
      first_idx = findfirst(vals, first_val)
      # if the index is not found, that means the input is not monotonic.
      # We would have to search the earlier blocks for the matching point
      if first_idx != 0
        found_block = block
        found_idx = first_idx
        break
      end
    end

    # construct the offset in the frame of reference of older block
    dest_rng = found_idx:(found_idx + length(rng) - 1)

    # shift it to the global frame of reference
    dest_rng += block_offsets[found_block]
    dest_rnges[dest_rng_idx] = dest_rng
    dest_rng_idx += 1

    # update the offset for the next iteration
    block_offset = block_offsets[found_block] + found_idx - 1
  end

  # TODO: check that the highest value is in the last block
  # I think that might break the algorithm

  vals = sub(data, :, col)
  max_val = maximum(vals)

  vals2 = sub(data, rnges[end], col)
  max_val2 = maximum(vals)

  if max_val != max_val2
    throw(ErrorException("Last interval does not contain the maximum index"))
  end

  # remove any blocks that exceed the length of the destination array
  nrows = last(dest_rnges[end])
  for i=1:nblocks
    rng = dest_rnges[i]
    if last(rng) > nrows
      dest_rnges[i] = -1:-1
    end
  end

  return dest_rnges
end

"""
  Copies the blocks from the source array into the destination array, which
  is returned

  **Inputs**

   * data: the data
   * src_rnges: monotonic ranges of data
   * dest_rnges: image of each src range in the destination array
   * col: column of data that is monotonic

  **Outputs**

   * data2: data extracted from data such that it is monotonic in the column
            col
"""
function copy_blocks(data::AbstractArray,
                     src_rnges::Array{UnitRange{Int}},
                     dest_rnges::Array{UnitRange{Int}}, col::Integer)

  # copy the blocks to the new array
  nrows = last(dest_rnges[end])
  data2 = zeros(eltype(data), nrows, size(data, 2))

  nblocks = length(src_rnges)

  for i=1:nblocks
    if dest_rnges[i] == -1:-1
      continue
    end
    src_block = sub(data, src_rnges[i], :)
    dest_block = sub(data2, dest_rnges[i], :)
    copy!(dest_block, src_block)
  end

  return data2
end

"""
  Checks that sections of blocks that are mapped to the same rows in the
  destination array are the same, throws an exception otherwise

  **Inputs**

   * data: the data
   * src_rnges: monotonic ranges of data
   * dest_rnges: image of each src range in the destination array
   * col: column of data that is monotonic
"""
function check_overlap(data::AbstractArray,
                       src_rnges::Array{UnitRange{Int}},
                       dest_rnges::Array{UnitRange{Int}}, col::Integer)

  nblocks = length(src_rnges)
  for i=1:nblocks
#    println("i = ", i)
    src_rng_i = src_rnges[i]
    dest_rng_i = dest_rnges[i]

    if dest_rng_i == -1:-1
      continue
    end
    for j=1:nblocks
      src_rng_j = src_rnges[j]
      dest_rng_j = dest_rnges[j]

      if dest_rng_j == -1:-1
        continue
      end

      overlap = intersect(dest_rng_i, dest_rng_j)
      overlap_start = first(overlap)
      overlap_last = last(overlap)

      if length(overlap) == 0
        continue
      end

      # figure out image of overlap in the two ranges
      begin_offset_i = overlap_start - first(dest_rng_i)
      end_offset_i = last(dest_rng_i) - overlap_last
      image_i = (first(src_rng_i) + begin_offset_i):(last(src_rng_i) - end_offset_i)
      @assert length(image_i) == length(overlap)

      begin_offset_j = overlap_start - first(dest_rng_j)
      end_offset_j = last(dest_rng_j) - overlap_last
      image_j = (first(src_rng_j) + begin_offset_j):(last(src_rng_j) - end_offset_j)

      @assert length(image_j) == length(overlap)
      compare_blocks(data, image_i, image_j, col, 1e-13)

    end
  end

  return nothing
end

"""
  Verifies that the specified rows in data are identical.

  **Inputs**

   * data: the data
   * image_i: first range of rows
   * image_j: second set of rows
   * col: column used to short
   * tol: tolerance for equality test

  **Outputs**

   * none

  This function throws an exception if the rows are not identical.
"""
function compare_blocks(data::AbstractArray, image_i::UnitRange{Int},
                        image_j::UnitRange{Int}, col::Integer,
                        tol::Number)

  @assert length(image_i) == length(image_j)
  # compare col, then the rest of the block
#=
  failflag1 = false
  for i=1:length(image_i)
    diff = abs(data[image_i[i], col] - data[image_j[i], col])
    failflag1 =  (diff > tol) || failflag1
#    @assert diff < tol
  end

  if failflag1
    throw(ErrorException("overlap does not refer to same range"))
  end
=#
  failflag = false
  for j=1:size(data, 2)
    for i=1:length(image_i)
      val = abs(data[ image_i[i], j] - data[ image_j[i], j] ) > tol
      failflag = val || failflag
    end
  end

  if failflag
    println(STDERR, "lines ", image_i, " and ", image_j, " should be identical but are not")
    throw(ErrorException(""))
  end

  return nothing
end

#------------------------------------------------------------------------------
# Run the code

if length(ARGS) < 2 || length(ARGS) > 3
  println(STDERR, "Usage: julia clean_log.jl fname column [output_fname]")
end

if length(ARGS) == 2
  fstem, fext = split_fname(ARGS[1])
  output_fname = string(fstem, "_cleaned", fext)
else
  output_fname = ARGS[3]
end

clean_log(ARGS[1], parse(Int, ARGS[2]), output_fname)


text:
# This is part of Quiz.jl. License is MIT.

using Quiz
using Test

println("Quiz type testing...")

#----------------------------------------------------------------------------------------------#
#                                     Type hyeranchy tests                                     #
#----------------------------------------------------------------------------------------------#

@testset "Quiz.type.abstract:                                                     " begin
    @test qEngine <: QUIZ
    @test qProblem <: QUIZ
    @test qElem <: qProblem
    @test qT <: qElem
    @test qC <: qElem
    @test qV <: qElem
end

#----------------------------------------------------------------------------------------------#
#                                       Type usage tests                                       #
#----------------------------------------------------------------------------------------------#

@testset "Quiz.type.utils:                                                        " begin
    the_f = [Float16(0.1), 0.1f0, 0.1, π, 2//3]
    the_i = [0x08, 0x0016, 0x00000032, Int8(8), Int16(16), Int32(32), 64]
    the_u = [u"m", u"kJ/kg", u"°C", u"K", unit(1)]
    for U in the_u
        for I in the_i
            v, u = Quiz.val_unit(I*U)
            @test v isa N64
            @test v isa Int64
            @test v == Int64(I)
            @test u == U
            @test v*u isa Q64 || v*u isa N64
            @test v*u isa A64
        end
        for F in the_f
            v, u = Quiz.val_unit(F*U)
            @test v isa N64
            @test v isa Float64
            @test v ≈ float(F)
            @test u == U
            @test v*u isa Q64 || v*u isa N64
            @test v*u isa A64
        end
    end
end

@testset "Quiz.type.qT instantiation:                                             " begin
    t = "This is a test"
    a = qT(t)
    @test t == a()
end

@testset "Quiz.type.qC instantiation:                                             " begin
    the_f = [Float16(0.1), 0.1f0, 0.1, π, 2//3]
    the_i = [0x08, 0x0016, 0x00000032, Int8(8), Int16(16), Int32(32), 64]
    the_u = [u"m", u"kJ/kg", u"°C", u"K", unit(1)]
    for U in the_u
        for V in vcat(the_i, the_f)
            # Two-argument constructor
            c = qC(:test, V*U)
            @test c isa qC
            @test c.nam isa Symbol
            @test c.val isa N64
            @test c.uni isa Unitful.Units
            @test c.val*c.uni ≈ V*U
            # Three-argument constructor
            c = qC(:test, V, U)
            @test c isa qC
            @test c.nam isa Symbol
            @test c.val isa N64
            @test c.uni isa Unitful.Units
            @test c.val*c.uni ≈ V*U
            # Three-argument constructor, unit misplaced
            c = qC(:test, V*U, unit(1))
            @test c isa qC
            @test c.nam isa Symbol
            @test c.val isa N64
            @test c.uni isa Unitful.Units
            @test c.val*c.uni ≈ V*U
            # Equivalence tests
            @test qC(:test, V*U) == qC(:test, V, U)
            @test qC(:test, V*U) == qC(:test, V*U, unit(1))
        end
    end
end

@testset "Quiz.type.qC usage:                                                     " begin
    the_f = [Float16(0.1), 0.1f0, 0.1, π, 2//3]
    the_i = [0x08, 0x0016, 0x00000032, Int8(8), Int16(16), Int32(32), 64]
    the_u = [u"m", u"kJ/kg", u"°C", u"K", unit(1)]
    for U in the_u
        for V in vcat(the_i, the_f)
            # Functor, the, and iE:
            c = qC(:test, V*U)
            @test c() isa Pair{Symbol, <:Q64}
            @test the(c) isa Q64
            @test the(c) == Quantity{typeof(c.val), dimension(c.uni), typeof(c.uni)}(V*U)
            @test iE(c) isa Expr
            @test iE(c).head == :(=)
            @test iE(c).args[1] == c.nam == :test
            @test iE(c).args[2] == the(c)
            eval(iE(c))
            @test test isa Q64
            @test test == the(c)
        end
    end
end

@testset "Quiz.type.qV instantiation:                                             " begin
    the_f = [Float16(0.1), 0.1f0, 0.1, π, 2//3]
    the_i = [0x08, 0x0016, 0x00000032, Int8(8), Int16(16), Int32(32), 64]
    the_u = [u"m", u"kJ/kg", u"°C", u"K", unit(1)]
    for U in the_u
        for V in vcat(the_i, the_f)
            VEC = V * (10:2:20)
            # Two-argument constructor
            v = qV(:test, collect(VEC)*U)
            @test v isa qV
            @test v.nam isa Symbol
            @test v.val isa Vector{<:N64}
            @test v.uni isa Unitful.Units
            @test v.val ≈ collect(VEC)
            # Three-argument constructor
            v = qV(:test, VEC, U)
            @test v isa qV
            @test v.nam isa Symbol
            @test v.val isa Vector{<:N64}
            @test v.uni isa Unitful.Units
            @test v.val ≈ collect(VEC)
            # Three-argument constructor, unit misplaced
            v = qV(:test, collect(VEC)*U, unit(1))
            @test v isa qV
            @test v.nam isa Symbol
            @test v.val isa Vector{<:N64}
            @test v.uni isa Unitful.Units
            @test v.val ≈ collect(VEC)
            # Equivalence tests
            @test qV(:test, collect(VEC)*U).val ≈ qV(:test, VEC, U).val
            @test qV(:test, collect(VEC)*U).val ≈ qV(:test, collect(VEC)*U, unit(1)).val
        end
    end
end

@testset "Quiz.type.qV usage:                                                     " begin
    the_f = [Float16(0.1), 0.1f0, 0.1, π, 2//3]
    the_i = [0x08, 0x0016, 0x00000032, Int8(8), Int16(16), Int32(32), 64]
    the_u = [u"m", u"kJ/kg", u"°C", u"K", unit(1)]
    for U in the_u
        for V in vcat(the_i, the_f)
            # Functor, the, and iE:
            VEC = collect(V * (10:2:20))
            v = qV(:test, VEC*U)
            @test v() isa Pair{Symbol, Vector{Quantity{T,D,U}}} where {T,D,U}
            @test the(v) isa Vector{<:Q64}
            @test [i.val for i in the(v)] ≈ VEC
            @test iE(v) isa Expr
            @test iE(v).head == :(=)
            @test iE(v).args[1] == v.nam == :test
            eval(iE(v))
            @test test isa Q64
            @test test ∈ the(v)
        end
    end
end


text:
const SETTINGS_FILE_NAME = "TexAssistant.toml"

is_setting_file(file) = isfile(file) && endswith(basename(file), SETTINGS_FILE_NAME)

function current_proj(path=pwd())
    !isdir(path) && return path
    if !(ishomedir(path) || isrootdir(path))
        setfile = joinpath(path, SETTINGS_FILE_NAME)
        isfile(setfile) && return setfile
        return current_proj(dirname(path))
    end
    return path
end

function load_settings(;path = pwd())
    setfile = current_proj(path)
    is_setting_file(setfile) ? TOML.parsefile(setfile) : Dict{String, Any}()
end

function load_setting(k, dflt; path = pwd())
    sets = load_settings(;path)
    get(sets, k, dflt)
end

load_setting(k; path=pwd()) = load_setting(k, ""; path)


text:
"""
Compatibility layer for transferring from FixedSizeArrays. This provides
alternative definitions of `Vec`, `Mat`, `Point`, `FixedVectorNoTuple`, `@fsa`,
etc, using StaticArrays as a backend.
The type definitions are not "perfect" matches because the type parameters are
different. However, it should cover common method signatures and constructors.
"""
module FixedSizeArrays

using StaticArrays

export FixedArray
export FixedVector
export FixedMatrix
export Mat
export Vec
export Point
export @fsa
export FixedVectorNoTuple

const FixedArray = StaticArray
const FixedVector = StaticVector
const FixedMatrix = StaticMatrix
const Mat = SMatrix
const FixedVectorNoTuple = FieldVector

function fsa_ast(ex)
    @assert isa(ex, Expr)
    if ex.head == :vect # Vector
        return Expr(:call, SVector{length(ex.args)}, Expr(:tuple, ex.args...))
    elseif ex.head == :hcat # 1 x n
        s1 = 1
        s2 = length(ex.args)
        return Expr(:call, SMatrix{s1, s2}, Expr(:tuple, ex.args...))
    elseif ex.head == :vcat
        if isa(ex.args[1], Expr) && ex.args[1].head == :row # n x m
            # Validate
            s1 = length(ex.args)
            s2s = map(i -> ((isa(ex.args[i], Expr) && ex.args[i].head == :row) ? length(ex.args[i].args) : 1), 1:s1)
            s2 = minimum(s2s)
            if maximum(s2s) != s2
                error("Rows must be of matching lengths")
            end

            exprs = [ex.args[i].args[j] for i = 1:s1, j = 1:s2]
            return Expr(:call, SMatrix{s1, s2}, Expr(:tuple, exprs...))
        else # n x 1
            return Expr(:call, SMatrix{length(ex.args), 1}, Expr(:tuple, ex.args...))
        end
    end
end
macro fsa(ex)
    expr = fsa_ast(ex)
    esc(expr)
end

macro fixed_vector(name, parent)
    esc(quote
        struct $(name){S, T} <: $(parent){S, T}
            data::NTuple{S, T}

            function $(name){S, T}(x::NTuple{S,T}) where {S, T}
                new{S, T}(x)
            end

            function $(name){S, T}(x::NTuple{S,Any}) where {S, T}
                new{S, T}(StaticArrays.convert_ntuple(T, x))
            end
        end
        size_or(::Type{$(name)}, or) = or
        eltype_or(::Type{$(name)}, or) = or
        eltype_or(::Type{$(name){S, T} where S}, or) where {T} = T
        eltype_or(::Type{$(name){S, T} where T}, or) where {S} = or
        eltype_or(::Type{$(name){S, T}}, or) where {S, T} = T

        size_or(::Type{$(name){S, T} where S}, or) where {T} = or
        size_or(::Type{$(name){S, T} where T}, or) where {S} = Size{(S,)}()
        size_or(::Type{$(name){S, T}}, or) where {S, T} = (S,)
        # Array constructor
        @inline function $(name){S}(x::AbstractVector{T}) where {S, T}
            @assert S <= length(x)
            $(name){S, T}(ntuple(i-> x[i], Val(S)))
        end
        @inline function $(name){S, T1}(x::AbstractVector{T2}) where {S, T1, T2}
            @assert S <= length(x)
            $(name){S, T1}(ntuple(i-> T1(x[i]), Val(S)))
        end

        @inline function $(name){S, T}(x) where {S, T}
            $(name){S, T}(ntuple(i-> T(x), Val(S)))
        end

        @inline function $(name){S}(x::T) where {S, T}
            $(name){S, T}(ntuple(i-> x, Val(S)))
        end
        @inline function $(name){1, T}(x::T) where T
            $(name){1, T}((x,))
        end
        @inline $(name)(x::NTuple{S}) where {S} = $(name){S}(x)
        @inline $(name)(x::T) where {S, T <: Tuple{Vararg{Any, S}}} = $(name){S, StaticArrays.promote_tuple_eltype(T)}(x)
        @inline function $(name){S}(x::T) where {S, T <: Tuple}
            $(name){S, StaticArrays.promote_tuple_eltype(T)}(x)
        end
        $(name){S, T}(x::FixedVector) where {S, T} = $(name){S, T}(Tuple(x))
        @generated function (::Type{$(name){S, T}})(x::$(name)) where {S, T}
            idx = [:(x[$i]) for i = 1:S]
            quote
                $($(name)){S, T}($(idx...))
            end
        end
        @generated function convert(::Type{$(name){S, T}}, x::$(name)) where {S, T}
            idx = [:(x[$i]) for i = 1:S]
            quote
                $($(name)){S, T}($(idx...))
            end
        end
        @generated function (::Type{SV})(x::FixedVector) where SV <: $(name)
            len = size_or(SV, size(x))[1]
            if length(x) == len
                :(SV(Tuple(x)))
            elseif length(x) > len
                elems = [:(x[$i]) for i = 1:len]
                :(SV($(Expr(:tuple, elems...))))
            else
                error("Static Vector too short: $x, target type: $SV")
            end
        end

        Base.@pure StaticArrays.Size(::Type{$(name){S, Any}}) where {S} = Size(S)
        Base.@pure StaticArrays.Size(::Type{$(name){S, T}}) where {S,T} = Size(S)

        Base.@propagate_inbounds function Base.getindex(v::$(name){S, T}, i::Int) where {S, T}
            v.data[i]
        end
        @inline Base.Tuple(v::$(name)) = v.data
        @inline Base.convert(::Type{$(name){S, T}}, x::NTuple{S, T}) where {S, T} = $(name){S, T}(x)
        @inline function Base.convert(::Type{$(name){S, T}}, x::Tuple) where {S, T}
            $(name){S, T}(convert(NTuple{S, T}, x))
        end

        @generated function StaticArrays.similar_type(::Type{SV}, ::Type{T}, s::Size{S}) where {SV <: $(name), T, S}
            if length(S) === 1
                $(name){S[1], T}
            else
                StaticArrays.default_similar_type(T,s(),Val{length(S)})
            end
        end

    end)
end

function Base.isnan(x::FixedArray)
    for elem in x
        isnan(elem) && return true
    end
    false
end

function unit(::Type{T}, i::Integer) where T <: FixedVector
    T(ntuple(Val(length(T))) do j
        ifelse(i == j, 1, 0)
    end)
end

export unit

function Base.extrema(a::AbstractVector{T}) where T <: FixedVector
    ET = eltype(T)
    reduce((x, v)-> (min.(x[1], v), max.(x[2], v)), a; init = (T(typemax(ET)), T(typemin(ET))))
end
function Base.minimum(a::AbstractVector{T}) where T <: FixedVector
    reduce((x, v)-> min.(x[1], v), a; init=T(typemax(eltype(T))))
end
function Base.maximum(a::AbstractVector{T}) where T <: FixedVector
    reduce((x, v)-> max.(x[1], v), a; init=T(typemin(eltype(T))))
end

@fixed_vector Vec FixedVector
@fixed_vector Point FixedVector

end


text:
"""
    get_boards(cred::TrelloCred)

List all Trello boards on your account.
(On the account assoicated with the credentials.)
"""
function get_boards(cred::TrelloCred)
    boards = get_request(cred, "/1/members/me/boards")
    # Technically board name is not garenteed to be unique, but the edge case where it
    # is not is so rare that I'ld rather have this cleaner API.
    return indexed_collection(boards)
end

"""
    create_board(cred::TrelloCred, name; desc="", id_organization="")

Create a new board on the `cred` account.
Board is created completely empty (no lists, nor labels).
If `id_organization` is provided then it will be created within that organization.
"""
function create_board(cred::TrelloCred, name; desc="", id_organization="")
    post_request(cred, "/1/boards/",
        name=name,
        desc=desc,
        idOrganization=id_organization,
        defaultLists=false,
        defaultLabels=false,
    )
end

"""
    delete_board(cred::TrelloCred, board_id|board)

Deletes the board.
"""
delete_board(cred::TrelloCred, board_id::AbstractString) = delete_request(cred, "/1/boards/$board_id")
delete_board(cred::TrelloCred, board) = delete_board(cred, board.id)

"""
    create_list(credit::TrelloCred, board_id|board, name, pos="bottom")

Create a list within the given board, with the given `name`,
at position given by `pos`
"""
function create_list(cred::TrelloCred, board_id::AbstractString, name, pos="bottom")
    resp = post_request(cred, "/1/boards/$(board_id)/lists"; name=name, pos=pos)
    return resp
end
create_list(cred::TrelloCred, board, args...) = create_list(cred, board.id, args...)

"""
    get_lists(cred::TrelloCred, board_id|board)

Get all the lists within a given Trello board.
"""
function get_lists(cred::TrelloCred, board_id::AbstractString)
    raw = get_request(cred, "/1/boards/$(board_id)/lists")
    return indexed_collection(sortbypos(filteroutclosed(raw)))
end
get_lists(cred::TrelloCred, board) = get_lists(cred, board.id)

"""
    get_labels(cred::TrelloCred, board_id|board)

Get all the labels within a given Trello board.
"""
function get_labels(cred::TrelloCred, board_id::AbstractString)
    raw = get_request(cred, "/1/boards/$(board_id)/labels")
    return indexed_collection(raw)
end
get_labels(cred::TrelloCred, board) = get_labels(cred, board.id)

"""
    create_card(cred::TrelloCred, list_id|list, name; desc="", label_ids::Vector=[])

Create a card, on the given list with the name, description and labels as specified.

See Trello docs:
 - https://developers.trello.com/reference#cards-2
 - https://api.trello.com/1/cards?idList=idList&keepFromSource=all
"""
function create_card(cred::TrelloCred, list_id::AbstractString, name; desc="", label_ids::Vector=[])
    return post_request(cred, "/1/cards";
        idList=list_id,
        name =name,
        desc=desc,
        idLabels = join(label_ids, ",")
    )
end
function create_card(cred::TrelloCred, list, args...; kwargs...)
    return create_card(cred, list.id, args...; kwargs...)
end

"""
    get_cards(cred::TrelloCred, list_id|list)

Returns all the cards on a given list.
"""
function get_cards(cred::TrelloCred, list_id::AbstractString)
    raw = get_request(cred, "/1/lists/$(list_id)/cards")
    cards = sortbypos(filteroutclosed(raw))
    return indexed_collection(cards)
end
get_cards(cred::TrelloCred, list) = get_cards(cred, list.id)


text:
module TerminalSpinners

using Printf

export Spinner, autospin, spin!, @spin

using Crayons

include("AnsiCodes.jl")
#=
     (minutes, seconds) = fldmod(elapsed, 60)
        (hours, minutes) = fldmod(minutes, 60)
        if hours == 0
            printstyled(s.stream, @sprintf("[%02d:%02d]", minutes, seconds); color=Base.info_color())
        else
            printstyled(s.stream, @sprintf("[%02d:%02d:%02d]", hours, minutes, seconds); color=Base.info_color())
        end
=#

ESC = "\e"
CSI = "\e["
DEC_RST  = "l"
DEC_SET  = "h"
DEC_TCEM = "?25"

function remove_ansi_characters(str::String)
    r = r"\x1B(?:[@-Z\\-_]|\[[0-?]*[ -/]*[@-~])"
    return replace(str, r => "")
end

const TTY_HIDE_CURSOR = CSI * DEC_TCEM * DEC_RST
const TTY_SHOW_CURSOR = CSI * DEC_TCEM * DEC_SET
const TTY_CLEAR_LINE = CSI * "2K"

include("spinners.jl")

tostring(msg) = msg
tostring(f::Function) = f()

Base.@kwdef mutable struct Spinner{IO_t <: IO}
    frames::Vector{String} = FORMATS[:dots][:frames]
    freq::Float64 = 10.0 # [1/s]
    msg::Any = ""
    stream::IO_t = stderr
    timer::Union{Nothing, Timer} = nothing
    hidecursor::Bool = true
    silent::Bool=false
    enabled::Bool= silent ? false : stream isa Base.TTY && !haskey(ENV, "CI")
    frame_idx::Int=1
    indent::Int=0
    color::Union{Crayon, Symbol} = Base.info_color()
end

getframe(s::Spinner) = s.frames[s.frame_idx]

function advance_frame!(s::Spinner)
    frame = s.frames[s.frame_idx]
    s.frame_idx = s.frame_idx == length(s.frames) ? 1 : s.frame_idx + 1
    return frame
end

function erase_and_reset(io::IO)
    AnsiCodes.erase_line!(io)
    AnsiCodes.cursor_horizontal_absolute!(io, 1)
end

function print_spinner(s::Spinner, io::IO)
    color = s.color
    if color isa Symbol
        printstyled(io, getframe(io); color=color)
    else
        print(io, color(getframe(io)))
    end
end

function next_line!(s::Spinner)
    return sprint(; context=s.stream) do io
        erase_and_reset(io)
        print(io, " "^s.indent)
        printstyled(io, getframe(s))
        msg = tostring(s.msg)
        print(io, " ", msg)
    end
end

function render!(s::Spinner)
    s.silent && return
    print(s.stream, next_line!(s))
    advance_frame!(s)
end

function start!(s::Spinner)
    s.silent && return

    if !s.enabled
        println(s.stream, s.msg)
        return
    end

    s.hidecursor && print(s.stream, TTY_HIDE_CURSOR)

    t = Timer(0.0; interval=1/s.freq) do timer
        try
            render!(s)
        catch e
            close(timer)
            @show e
        end
    end
    s.timer = t
    return s
end

function stop!(s::Spinner)
    if s.timer !== nothing
        close(s.timer)
    end
    s.frame_idx = 1
    if !s.enabled || s.silent
        return
    end
    s.hidecursor && print(s.stream, TTY_SHOW_CURSOR)
end

function success!(s)
    stop!(s)
    return sprint(; context=s.stream) do io
        erase_and_reset(io)
        print(io, " "^s.indent)
        print(tostring(s.success))
        msg = tostring(s.msg)
        print(io, " ", msg)
    end
    stop!(s)
    print(s.stream, tostring(msg))
end

macro spin(s, work)
    return quote
        spin(() -> $(esc(work)), $(esc(s)))
    end
end

function spin(f, s::Spinner)
    start!(s)
    try
        v = f()
        success!(s)
        v
    catch
        fail!(s)
        rethrow()
    end
end

end # module


text:
using Test, HTTP, JSON, Sockets
using bjoernfalconserver: @js

function test_server(host, port)

make_request(endpoint, body) = body === nothing ?
    HTTP.request("GET", "http://$host:$port/$endpoint") :
    HTTP.request("POST", "http://$host:$port/$endpoint", ["Content-Type" => "application/json"], body)

# test on response body
function test_response(f::Function, endpoint, body = nothing)
    r = make_request(endpoint, body)
    @test r.status == 200
    @test !isempty(r.body)
    body = String(r.body)
    f(body)
end

test_response("ping", nothing) do body
    @test body == "pong"
end

let data = @js { foo: "bar" }
    test_response("ping", data) do body
        @test body == @js { data }
    end
end

end # end test_server

function tryconnect(host, port)
    try
        return (Sockets.connect(host, port); true)
    catch err
        err isa Base.IOError && occursin(r"ECONNREFUSED", err.msg) && return false
        rethrow(err)
    end
end

@testset "server" begin
    host = get(ENV, "PYCALLSERVER_HOST", "0.0.0.0")
    port = parse(Int, get(ENV, "PYCALLSERVER_PORT", "8000"))
    while !tryconnect(host, port) end
    test_server(host, port)
end


text:
const SNV_SUBSTITUTION = 0.98
const SNV_INSERTION = 0.01
const SNV_DELETION = 0.01
const per_base_qual_adjust = rand(5, 1000)

# simulate the sequencing process based on the config
# return a fastq pair
function pair_end_seq(dna_template, config)
    const fixed_index = "ATCGATCG"
    index1 = fixed_index
    index2 = fixed_index
    if config["random_index1"]
        index1 = random_index()
    end
    if config["random_index2"]
        index2 = random_index()
    end
    # simulate the duplication rate
    read_num = Int(round((rand() * config["duplication_rate"])/0.5))
    reads = []
    for i in 1:read_num
        r1seq, r1qual = sequence_simulation(dna_template, config)
        r2seq, r2qual = sequence_simulation(~dna_template, config)
        r1name, r2name = name_simulation(true, index1, index2)
        read1 = FastqRead(r1name, r1seq, "+", r1qual)
        read2 = FastqRead(r2name, r2seq, "+", r2qual)
        pair = FastqPair(read1, read2)
        push!(reads, pair)
    end
    return reads
end

function random_index(index_len = 8)
    atcg = ['A', 'T', 'C', 'G']
    rnd = round(rand(index_len) * 3) + 1
    arr = [atcg[Int(rnd[i])] for i in 1:index_len]
    return join(arr)
end

function name_simulation(pairend, index1, index2)
    const device = "SEQMAKER"
    const run = 1
    const chip = SEQMAKER_VERSION
    lane = Int(round(rand() * 3)) + 1
    tile = Int(round(rand() * 99)) + 1
    x = Int(round(rand() * 9999)) + 1
    y = Int(round(rand() * 9999)) + 1

    firstpart = "@$device:$run:$chip:$lane:$tile:$x:$y"
    if pairend
        return "$firstpart 1:N:0:$index1", "$firstpart 2:N:0:$index2"
    else
        return "$firstpart 1:N:0:$index1"
    end
end

function sequence_simulation(dna_template, config)
    seq_arr = Char[ nt for nt in dna_template.seq]
    readlen = min(length(dna_template), config["readlen"])

    # simulate quality raise (by coordination) and quality drop
    const qual_curve = [0.95, 1.06, 1.14, 1.16, 1.168,
        1.176, 1.173, 1.172, 1.172, 1.170,
        1.168, 1.165, 1.163,
        1.159, 1.158, 1.153,
        1.148, 1.150, 1.143,
        1.137, 1.133,
        1.127, 1.130, 1.120,
        1.116, 1.112,
        1.105,
        1.096, 1.097, 1.093,
        1.087, 1.083, 1.081,
        1.079, 1.075, 1.073,
        1.066, 1.063,
        1.054, 1.051,
        1.046,
        1.038, 1.03, 1.021, 1.01, 1.02, 1.0]
    const base_index = Dict('A'=>1, 'T'=>2, 'C'=>3, 'G'=>4, 'N'=>5)
    # simulate last base with low quality
    const qual_last_base = 0.9
    qual_adjust_ratio = [0.0 for i in 1:readlen]
    # interpolation with qual_curve
    for i in 1:readlen
        if i == readlen
            qual_adjust_ratio[i] = qual_last_base
            continue
        end
        curve_len = length(qual_curve)
        pos = curve_len*i/readlen
        p = pos - floor(pos)
        left = min(curve_len, Int(floor(pos)) + 1)
        right = min(curve_len, left + 1)
        # linear interpolation
        qual_adjust_ratio[i] = qual_curve[left] * (1-p) + qual_curve[right] * p
    end

    rand_err = rand(readlen)
    has_error = false
    for i in 1:readlen
        if rand_err[i] < config["seq_error_rate"] / qual_adjust_ratio[i]
            has_error = true
            break
        end
    end

    rand_qual = rand(readlen)
    min_qual = config["normal_base_qual"]["min"]
    max_qual = config["normal_base_qual"]["max"]
    quals = rand_qual * min_qual + (1.0 - rand_qual) * max_qual
    # adjust the quality and make phred qual string
    qual_arr = [Char(Int( round(quals[i] * qual_adjust_ratio[i] * (per_base_qual_adjust[base_index[seq_arr[i]], i]*0.02+0.99) ) )+33) for i in 1:readlen]

    if !has_error
        sequence = Sequence(ASCIIString(seq_arr[1:readlen]))
        quality = Quality(ASCIIString(qual_arr))
        return sequence, quality
    else
        seqlen = 0
        dnapos = 0
        qualpos = 0
        for seq_err in rand_err
            # no error for this base
            if seq_err >= config["seq_error_rate"]
                seqlen += 1
                dnapos += 1
                qualpos += 1
                seq_arr[seqlen] = dna_template.seq[dnapos]
            else
                err = rand()
                # substitution
                if err < SNV_SUBSTITUTION
                    seqlen += 1
                    dnapos += 1
                    qualpos += 1
                    seq_arr[seqlen] = get_other_base(dna_template.seq[dnapos])
                    qual_arr[seqlen] = simulate_error_qual(config["seq_error_qual"]["min"], config["seq_error_qual"]["max"])
                # insertion
                elseif err < SNV_SUBSTITUTION + SNV_INSERTION
                    seqlen += 1
                    seq_arr[seqlen] = get_rand_base()
                    qual_arr[seqlen] = simulate_error_qual(config["seq_error_qual"]["min"], config["seq_error_qual"]["max"])
                # deletion
                elseif err < SNV_SUBSTITUTION + SNV_INSERTION + SNV_DELETION
                    dnapos += 1
                end
            end
            if seqlen>=readlen || dnapos>=length(dna_template)
                break
            end
        end
        sequence = Sequence(ASCIIString(seq_arr[1:seqlen]))
        quality = Quality(ASCIIString(qual_arr[1:seqlen]))
        return sequence, quality
    end
end

function simulate_error_qual(minqual, maxqual)
    rnd = rand()
    q = Int(round(minqual*rnd + maxqual * (1.0 - rnd)))
    return Char(q+33)
end

function get_other_base(base)
    const others = Dict(
        'A'=>['T', 'C', 'G'],
        'T'=>['A', 'C', 'G'],
        'C'=>['T', 'A', 'G'],
        'G'=>['T', 'C', 'A']
        )
    if !haskey(others, uppercase(base))
        return 'N'
    end
    rnd = Int(round(rand() * 2)) + 1
    return others[uppercase(base)][rnd]
end

function get_rand_base()
    const bases = ['A', 'T', 'C', 'G']
    rnd = Int(round(rand() * 3)) + 1
    return bases[rnd]
end


text:
using Millboard
using Test

Millboard.set_table_mode(:grid)
Base.eval(:(have_color = false))

board = ["Lorem ipsum dolor sit amet"]
@test """
+---+----------------------------+
|   |                          1 |
+===+============================+
| 1 | Lorem ipsum dolor sit amet |
+---+----------------------------+""" == table(board) |> string

board = ["Lorem ipsum\ndolor sit amet"]
@test """
+---+----------------+
|   |              1 |
+===+================+
|   |    Lorem ipsum |
| 1 | dolor sit amet |
+---+----------------+""" == table(board) |> string

@test """
+-------+----------------+
|       |              1 |
+=======+================+
|       |    Lorem ipsum |
| first | dolor sit amet |
+-------+----------------+""" == table(board, rownames=["first"]) |> string

@test """
+-----------+----------------+
|           |          first |
|           |         column |
+===========+================+
|           |    Lorem ipsum |
| first row | dolor sit amet |
+-----------+----------------+""" == table(board, colnames=["first\ncolumn"], rownames=["first row"]) |> string

board = ["Lorem ipsum\ndolor sit amet" 42]
@test """
+-----------+----------------+----+
|           |          first |  2 |
|           |         column |    |
+===========+================+====+
|           |    Lorem ipsum | 42 |
| first row | dolor sit amet |    |
+-----------+----------------+----+""" == table(board, colnames=["first\ncolumn"], rownames=["first row"]) |> string

Base.eval(:(have_color = true))


text:
using NumericalDebugging, Test, ForwardDiff

@testset "nonfinite" begin
    @test nonfinite(Inf)
    @test nonfinite(NaN)
    @test !nonfinite(-9)
    @test !nonfinite("a fish")
    @test nonfinite([1.0, Inf])
    @test !nonfinite([1.0, 2.0])
    @test nonfinite(ForwardDiff.Dual(1.0, Inf))
    @test nonfinite(ForwardDiff.Dual(NaN, 1.0))
    @test !nonfinite(ForwardDiff.Dual(0.0, 1.0))
end

let a = Inf
    @test_logs (:warn, "nonfinite: a") @check_if nonfinite a b = 9
    NumericalDebugging.ERROR[] = true
    @test_throws Exception @check_if nonfinite a
end


text:
function solve_010()
    vec = prime_sieve(2000000)
    return sum(vec)
end


text:
T=3
tsname="joDFT"
@testset "$tsname" begin
for t=1:T # start test loop
    m=4^t
    v1=rand(ComplexF64,m)
    mv1=rand(ComplexF64,m,m)
    v2=rand(ComplexF64,m,m)
    vv2=vec(v2)
    mv2=rand(ComplexF64,m,m,m)
    vmv2=reshape(mv2,m*m,m)

    A1=joDFT(m;DDT=ComplexF64)
    A2=joDFT(m,m;DDT=ComplexF64)
    verbose && println("$tsname ($m[,$m]) - planned")
    @testset "$m [x $m]" begin
        @test isadjoint(joDFT(m))[1]
        @test isadjoint(joDFT(m;centered=true))[1]
        @test islinear(joDFT(m))[1]
        @test islinear(joDFT(m;centered=true))[1]
        @test isadjoint(joDFT(m,m))[1]
        @test isadjoint(joDFT(m,m;centered=true))[1]
        @test islinear(joDFT(m,m))[1]
        @test islinear(joDFT(m,m;centered=true))[1]

        @test norm(A1*v1-fft(v1)/sqrt(m))<joTol
        @test norm(A1\v1-ifft(v1)*sqrt(m))<joTol
        @test norm(adjoint(A1)*v1-ifft(v1)*sqrt(m))<joTol
        @test norm((adjoint(A1)*A1)*v1-v1)<joTol
        @test norm(A1*mv1-fft(mv1,1)/sqrt(m))<joTol
        @test norm(A1\mv1-ifft(mv1,1)*sqrt(m))<joTol
        @test norm(adjoint(A1)*mv1-ifft(mv1,1)*sqrt(m))<joTol
        @test norm((adjoint(A1)*A1)*mv1-mv1)<joTol

        @test norm(A2*vv2-vec(fft(v2))/sqrt(m^2))<joTol
        @test norm(A2\vv2-vec(ifft(v2))*sqrt(m^2))<joTol
        @test norm(adjoint(A2)*vv2-vec(ifft(v2))*sqrt(m^2))<joTol
        @test norm((adjoint(A2)*A2)*vv2-vv2)<joTol
        @test norm(vec(A2*vmv2)-vec(fft(mv2,1:2))/sqrt(m^2))<joTol
        @test norm(vec(A2\vmv2)-vec(ifft(mv2,1:2))*sqrt(m^2))<joTol
        @test norm(vec(adjoint(A2)*vmv2)-vec(ifft(mv2,1:2))*sqrt(m^2))<joTol
        @test norm((adjoint(A2)*A2)*vmv2-vmv2)<joTol
    end

    A1=joDFT(m;planned=false,DDT=ComplexF64)
    A2=joDFT(m,m;planned=false,DDT=ComplexF64)
    verbose && println("$tsname ($m[,$m]) - not planned")
    @testset "$m [x $m]" begin
        @test isadjoint(joDFT(m;planned=false))[1]
        @test isadjoint(joDFT(m;planned=false,centered=true))[1]
        @test islinear(joDFT(m;planned=false))[1]
        @test islinear(joDFT(m;planned=false,centered=true))[1]
        @test isadjoint(joDFT(m,m;planned=false))[1]
        @test isadjoint(joDFT(m,m;planned=false,centered=true))[1]
        @test islinear(joDFT(m,m;planned=false))[1]
        @test islinear(joDFT(m,m;planned=false,centered=true))[1]

        @test norm(A1*v1-fft(v1)/sqrt(m))<joTol
        @test norm(A1\v1-ifft(v1)*sqrt(m))<joTol
        @test norm(adjoint(A1)*v1-ifft(v1)*sqrt(m))<joTol
        @test norm((adjoint(A1)*A1)*v1-v1)<joTol
        @test norm(A1*mv1-fft(mv1,1)/sqrt(m))<joTol
        @test norm(A1\mv1-ifft(mv1,1)*sqrt(m))<joTol
        @test norm(adjoint(A1)*mv1-ifft(mv1,1)*sqrt(m))<joTol
        @test norm((adjoint(A1)*A1)*mv1-mv1)<joTol

        @test norm(A2*vv2-vec(fft(v2))/sqrt(m^2))<joTol
        @test norm(A2\vv2-vec(ifft(v2))*sqrt(m^2))<joTol
        @test norm(adjoint(A2)*vv2-vec(ifft(v2))*sqrt(m^2))<joTol
        @test norm((adjoint(A2)*A2)*vv2-vv2)<joTol
        @test norm(vec(A2*vmv2)-vec(fft(mv2,1:2))/sqrt(m^2))<joTol
        @test norm(vec(A2\vmv2)-vec(ifft(mv2,1:2))*sqrt(m^2))<joTol
        @test norm(vec(adjoint(A2)*vmv2)-vec(ifft(mv2,1:2))*sqrt(m^2))<joTol
        @test norm((adjoint(A2)*A2)*vmv2-vmv2)<joTol
    end

end # end test loop
end


text:
export Dense

contract_t = 0.0

struct Dense{T} <: TensorStorage
  data::Vector{T}
  Dense{T}(data::Vector) where {T} = new{T}(convert(Vector{T},data))
  Dense{T}(size::Integer) where {T} = new{T}(Vector{T}(undef,size))
  Dense{T}(x::Number,size::Integer) where {T} = new{T}(fill(convert(T,x),size))
  Dense{T}() where {T} = new{T}(Vector{T}())
end

data(D::Dense) = D.data
length(D::Dense) = length(data(D))
eltype(D::Dense) = eltype(data(D))
Base.getindex(D::Dense,i::Integer) = getindex(data(D),i)
Base.setindex!(D::Dense,x,i::Integer) = setindex!(data(D),x,i)
#TODO: this should do proper promotions of the storage data
#e.g. ComplexF64*Dense{Float64} -> Dense{ComplexF64}
*(D::Dense{T},x::S) where {T,S<:Number} = Dense{promote_type(T,S)}(x*data(D))
*(x::Number,D::Dense) = D*x

Base.promote_type(::Type{Dense{T1}},::Type{Dense{T2}}) where {T1,T2} = Dense{promote_type(T1,T2)}

convert(::Type{Dense{T}},D::Dense) where {T} = Dense{T}(data(D))

storage_convert(::Type{Dense{T}},D::Dense,is::IndexSet) where {T} = convert(Dense{T},D)

# convert to complex
storage_complex(D::Dense{T}) where {T} = Dense{complex(T)}(complex(data(D)))

copy(D::Dense{T}) where {T} = Dense{T}(copy(data(D)))

function storage_outer(D1::Dense{T},is1::IndexSet,D2::Dense{S},is2::IndexSet) where {T, S <:Number}
  return Dense{promote_type(T,S)}(vec(data(D1)*transpose(data(D2)))),IndexSet(is1,is2)
end

storage_convert(::Type{Array},D::Dense,is::IndexSet) = reshape(data(D),dims(is))

function storage_convert(::Type{Array},
                         D::Dense,
                         ois::IndexSet,
                         nis::IndexSet)
  P = calculate_permutation(nis,ois)
  A = reshape(data(D),dims(ois))
  return permutedims(A,P)
end

storage_fill!(D::Dense,x::Number) = fill!(data(D),x)

function storage_getindex(Tstore::Dense{T},
                          Tis::IndexSet,
                          vals::Union{Int, AbstractVector{Int}}...) where {T}
  return getindex(reshape(data(Tstore),dims(Tis)),vals...)
end

function storage_setindex!(Tstore::Dense,
                           Tis::IndexSet,
                           x::Union{<:Number, AbstractArray{<:Number}},
                           vals::Union{Int, AbstractVector{Int}}...)
  return setindex!(reshape(data(Tstore),dims(Tis)),x,vals...)
end

# For use in _add!
using Base.Cartesian: @nexprs,
                      @ntuple,
                      @nloops

#
# A generalized permutedims!(P,B,perm) that also allows
# a function to be applied elementwise (defaults to addition of
# the elements of P and the permuted elements of B)
#
# Based off of the permutedims! implementation in Julia's base:
# https://github.com/JuliaLang/julia/blob/91151ab871c7e7d6689d1cfa793c12062d37d6b6/base/multidimensional.jl#L1355
#
@generated function _add!(P::Array{T,N},
                          B::Array{S,N},
                          perm,
                          f = (x,y)->x+y) where {T,S,N}
  quote
    Base.checkdims_perm(P, B, perm)

    #calculates all the strides
    native_strides = Base.size_to_strides(1, size(B)...)
    strides_1 = 0
    @nexprs $N d->(strides_{d+1} = native_strides[perm[d]])

    #Creates offset, because indexing starts at 1
    offset = 1 - sum(@ntuple $N d->strides_{d+1})

    ind = 1
    @nexprs 1 d->(counts_{$N+1} = strides_{$N+1}) # a trick to set counts_($N+1)
    @nloops($N, i, P,
            d->(counts_d = strides_d), # PRE
            d->(counts_{d+1} += strides_{d+1}), # POST
            begin # BODY
                sumc = sum(@ntuple $N d->counts_{d+1})
                @inbounds P[ind] = f(P[ind],B[sumc+offset])
                ind += 1
            end)

    return P
  end
end

function _add!(Bstore::Dense,
               Bis::IndexSet,
               Astore::Dense,
               Ais::IndexSet,
               x::Number = 1)
  p = calculate_permutation(Bis,Ais)
  Adata = data(Astore)
  Bdata = data(Bstore)
  if is_trivial_permutation(p)
    if x == 1
      Bdata .+= Adata
    else
      Bdata .= Bdata .+ x .* Adata
    end
  else
    reshapeBdata = reshape(Bdata,dims(Bis))
    reshapeAdata = reshape(Adata,dims(Ais))
    if x == 1
      _add!(reshapeBdata,reshapeAdata,p)
    else
      _add!(reshapeBdata,reshapeAdata,p,(a,b)->a+x*b)
    end
  end
end

function storage_copyto!(Bstore::Dense,
                         Bis::IndexSet,
                         Astore::Dense,
                         Ais::IndexSet,
                         x::Number = 1)
  p = calculate_permutation(Bis,Ais)
  Adata = data(Astore)
  Bdata = data(Bstore)
  if is_trivial_permutation(p)
    if x == 1
      Bdata .= Adata
    else
      Bdata .= x .* Adata
    end
  else
    reshapeBdata = reshape(Bdata,dims(Bis))
    reshapeAdata = reshape(Adata,dims(Ais))
    if x == 1
      _add!(reshapeBdata,reshapeAdata,p,(a,b)->b)
    else
      _add!(reshapeBdata,reshapeAdata,p,(a,b)->x*b)
    end
  end
end

function storage_mult!(Astore::Dense,
                       x::Number)
  Adata = data(Astore)
  rmul!(Adata, x)
end

function storage_mult(Astore::Dense,
                      x::Number)
  Bstore = copy(Astore)
  storage_mult!(Bstore, x)
  return Bstore
end

# For Real storage and complex scalar, promotion
# of the storage is needed
function storage_mult(Astore::Dense{T},
                      x::S) where {T<:Real,S<:Complex}
  Bstore = convert(Dense{promote_type(S,T)},Astore)
  storage_mult!(Bstore, x)
  return Bstore
end

# TODO: make this a special version of storage_add!()
# Make sure the permutation is optimized
function storage_permute!(Bstore::Dense, Bis::IndexSet,
                          Astore::Dense, Ais::IndexSet)
  p = calculate_permutation(Bis,Ais)
  Adata = data(Astore)
  Bdata = data(Bstore)
  if is_trivial_permutation(p)
    Bdata .= Adata
  else
    reshapeBdata = reshape(Bdata,dims(Bis))
    permutedims!(reshapeBdata,reshape(Adata,dims(Ais)),p)
  end
end

function storage_dag(Astore::Dense,Ais::IndexSet)
  return dag(Ais),storage_conj(Astore)
end

function storage_scalar(D::Dense)
  length(D)==1 && return D[1]
  throw(ErrorException("Cannot convert Dense -> Number for length of data greater than 1"))
end

function _contract(Cinds::IndexSet, Clabels::Vector{Int},
                   Astore::Dense{SA}, Ainds::IndexSet, Alabels::Vector{Int},
                   Bstore::Dense{SB}, Binds::IndexSet, Blabels::Vector{Int}) where {SA<:Number,SB<:Number}
  SC = promote_type(SA,SB)

  # Convert the arrays to a common type
  # since we will call BLAS
  Astore = convert(Dense{SC},Astore)
  Bstore = convert(Dense{SC},Bstore)

  Adims = dims(Ainds)
  Bdims = dims(Binds)
  Cdims = dims(Cinds)

  # Create storage for output tensor
  Cstore = Dense{SC}(prod(Cdims))

  Adata = reshape(data(Astore),Adims)
  Bdata = reshape(data(Bstore),Bdims)
  Cdata = reshape(data(Cstore),Cdims)

  if(length(Alabels)==0)
    contract_scalar!(Cdata,Clabels,Bdata,Blabels,Adata[1])
  elseif(length(Blabels)==0)
    contract_scalar!(Cdata,Clabels,Adata,Alabels,Bdata[1])
  else
    props = CProps(Alabels,Blabels,Clabels)
    compute!(props,Adata,Bdata,Cdata)
    _contract_dense_dense!(Cdata,props,Adata,Bdata)
  end

  return Cstore
end

function storage_svd(Astore::Dense{T},
                     Lis::IndexSet,
                     Ris::IndexSet;
                     kwargs...) where {T}
  maxdim::Int = get(kwargs,:maxdim,min(dim(Lis),dim(Ris)))
  mindim::Int = get(kwargs,:mindim,1)
  cutoff::Float64 = get(kwargs,:cutoff,0.0)
  absoluteCutoff::Bool = get(kwargs,:absoluteCutoff,false)
  doRelCutoff::Bool = get(kwargs,:doRelCutoff,true)
  utags::String = get(kwargs,:utags,"Link,u")
  vtags::String = get(kwargs,:vtags,"Link,v")
  fastSVD::Bool = get(kwargs,:fastSVD,false)

  if fastSVD
    MU,MS,MV = svd(reshape(data(Astore),dim(Lis),dim(Ris)))
  else
    MU,MS,MV = recursiveSVD(reshape(data(Astore),dim(Lis),dim(Ris)))
  end
  MV = conj!(MV)

  P = MS.^2
  #@printf "  Truncating with maxdim=%d cutoff=%.3E\n" maxdim cutoff
  truncate!(P;mindim=mindim,
              maxdim=maxdim,
              cutoff=cutoff,
              absoluteCutoff=absoluteCutoff,
              doRelCutoff=doRelCutoff)
  dS = length(P)
  if dS < length(MS)
    MU = MU[:,1:dS]
    resize!(MS,dS)
    MV = MV[:,1:dS]
  end

  u = Index(dS,utags)
  v = settags(u,vtags)
  Uis,Ustore = IndexSet(Lis...,u),Dense{T}(vec(MU))
  #TODO: make a diag storage
  Sis,Sstore = IndexSet(u,v),Diag{Vector{Float64}}(MS)
  Vis,Vstore = IndexSet(Ris...,v),Dense{T}(Vector{T}(vec(MV)))

  return (Uis,Ustore,Sis,Sstore,Vis,Vstore)
end

function storage_eigen(Astore::Dense{T},
                       Lis::IndexSet,
                       Ris::IndexSet;
                       kwargs...) where {T}
  maxdim::Int = get(kwargs,:maxdim,min(dim(Lis),dim(Ris)))
  mindim::Int = get(kwargs,:mindim,1)
  cutoff::Float64 = get(kwargs,:cutoff,0.0)
  absoluteCutoff::Bool = get(kwargs,:absoluteCutoff,false)
  doRelCutoff::Bool = get(kwargs,:doRelCutoff,true)
  tags::TagSet = get(kwargs,:lefttags,"Link,u")
  lefttags::TagSet = get(kwargs,:lefttags,tags)
  righttags::TagSet = get(kwargs,:righttags,prime(lefttags))

  dim_left = dim(Lis)
  dim_right = dim(Ris)
  MD,MU = eigen(Hermitian(reshape(data(Astore),dim_left,dim_right)))

  # Sort by largest to smallest eigenvalues
  p = sortperm(MD; rev = true)
  MD = MD[p]
  MU = MU[:,p]

  #@printf "  Truncating with maxdim=%d cutoff=%.3E\n" maxdim cutoff
  truncate!(MD;maxdim=maxdim,
              cutoff=cutoff,
              absoluteCutoff=absoluteCutoff,
              doRelCutoff=doRelCutoff)
  dD = length(MD)
  if dD < size(MU,2)
    MU = MU[:,1:dD]
  end

  #TODO: include truncation parameters as keyword arguments
  u = Index(dD,lefttags)
  v = settags(u,righttags)
  Uis,Ustore = IndexSet(Lis...,u),Dense{T}(vec(MU))
  Dis,Dstore = IndexSet(u,v),Diag{Vector{Float64}}(MD)
  return (Uis,Ustore,Dis,Dstore)
end

# TODO: move this to a general "linear_algebra.jl" file?
# i.e. a common place for custom linear algebra functionality
# for matrices
function polar(A::Matrix)
  U,S,V = svd(A) # calls LinearAlgebra.svd()
  return U*V',V*Diagonal(S)*V'
end

function storage_qr(Astore::Dense{T},
                    Lis::IndexSet,
                    Ris::IndexSet;
                    kwargs...) where {T}
  tags::TagSet = get(kwargs,:tags,"Link,u")
  dim_left = dim(Lis)
  dim_right = dim(Ris)
  MQ,MP = qr(reshape(data(Astore),dim_left,dim_right))
  dim_middle = min(dim_left,dim_right)
  u = Index(dim_middle,tags)
  #Must call Matrix() on MQ since the QR decomposition outputs a sparse
  #form of the decomposition
  Qis,Qstore = IndexSet(Lis...,u),Dense{T}(vec(Matrix(MQ)))
  Pis,Pstore = IndexSet(u,Ris...),Dense{T}(vec(Matrix(MP)))
  return (Qis,Qstore,Pis,Pstore)
end

function storage_polar(Astore::Dense{T},
                       Lis::IndexSet,
                       Ris::IndexSet) where {T}
  dim_left = dim(Lis)
  dim_right = dim(Ris)
  MQ,MP = polar(reshape(data(Astore),dim_left,dim_right))
  dim_middle = min(dim_left,dim_right)
  Uis = prime(Ris)
  Qis,Qstore = IndexSet(Lis...,Uis...),Dense{T}(vec(MQ))
  Pis,Pstore = IndexSet(Uis...,Ris...),Dense{T}(vec(MP))
  return (Qis,Qstore,Pis,Pstore)
end


text:
import Pkg; Pkg.activate("../../")
using Turing

# Simulate data.
N = 1000
mu_true = [5.0 10.0 2.0]
sigma_true = 0.6
y = randn(N, length(mu_true)) * sigma_true .+ mu_true  # size(y) = (1000, 3)

# Model definition.
@model function demo_model(y)
    K = size(y, 2)
    mu ~ filldist(Normal(0, 10), K)
    sigma ~ LogNormal(0, 2)
    for k in 1:K
        y[:, k] .~ Normal(mu[k], sigma)
    end
end

# Sample from posterior via NUTS.
nadapt, target_accept_ratio = (1000, 0.8)
iterations = 1000 + nadapt
chain = sample(demo_model(y), NUTS(nadapt, target_accept_ratio), iterations)

# Get sigma
chain[:sigma]  # works for univariate parameters
# 2-dimensional AxisArray{Float64,2,...} with axes:
#     :iter, 1:1:1000
#     :chain, 1:1
# And data, a 1000×1 Array{Float64,2}:
#  0.5891087703064164
#  0.5982078191819813
#  0.6014525731322462
#  0.6049135394153814
#  0.6135096911950264
#  ⋮
#  0.6045248767595608
#  0.6043856761562517
#  0.6197436168542185
#  0.6036546321663407

# FIXME: Get mu
# chain[:mu]  # broken

# Workarouund for getting array parameters
# function extract_array_param(chain, sym)
#     rgx = Regex("$(sym)\\[\\d+(,\\d+)*\\]")
#     syms = filter(x -> x != nothing, match.(rgx, String.(chain.name_map[1])))
#     syms = map(x -> Symbol(x.match), syms)
#     return chain[syms]
# end
#
# extract_array_param(chain, :mu)

# This has been fixed.
group(chain, :mu)

# FIXME: Get log-likelihood.
loglike = logprob"y=y | model=demo_model, chain=chain"  # broken

max_stars_repo_name: luiarthur/TuringBnpBenchmarks

max_issues_repo_name: luiarthur/TuringBnpBenchmarks

max_forks_repo_name: luiarthur/TuringBnpBenchmarks


text:
export AbstractSGStkeModel, NoSGStke, SGStkeModel

abstract type AbstractSGStkeModel end

struct NoSGStke <: AbstractSGStkeModel end
struct SGStkeModel <: AbstractSGStkeModel end

###########################
## default option: NoSGStke
###########################

vars_state(::AbstractSGStkeModel, ::AbstractStateType, FT) = @vars()

precompute(::AbstractSGStkeModel, atmos::AtmosModel, args, ts, ::Source) = NamedTuple()

function atmos_init_aux!(
    ::AbstractSGStkeModel,
    ::AtmosModel,
    aux::Vars,
    geom::LocalGeometry,
)
    nothing
end

function atmos_nodal_update_auxiliary_state!(
    ::AbstractSGStkeModel,
    m::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
)
    nothing
end

function compute_gradient_flux!(
    ::AbstractSGStkeModel,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    nothing
end

function compute_gradient_argument!(
    ::AbstractSGStkeModel,
    m::AtmosModel,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    nothing
end

function wavespeed_sgstke!(
    ::AbstractSGStkeModel,
    wavespeed::Vars,
    nM,
    state::Vars,
    aux::Vars,
    t::Real,
)
    nothing
end

###########################
## SGStkeModel
###########################

vars_state(sgstke::SGStkeModel, ::Prognostic, FT) = @vars(ρe_SGS::FT)
vars_state(sgstke::SGStkeModel, ::Gradient, FT) = @vars(e_SGS::FT, θ_liq_ice::FT)
vars_state(sgstke::SGStkeModel, ::GradientFlux, FT) =
    @vars(∇e_SGS::SVector{3, FT}, ∇θ_liq_ice::SVector{3, FT}, ∇u::SMatrix{3, 3, FT, 9})
vars_state(sgstke::SGStkeModel, ::Auxiliary, FT) = @vars(Δ::FT)

function precompute(::SGStkeModel, atmos::AtmosModel, args, ts, ::Source)
    @unpack state, aux, diffusive, t = args
    ν, D_t, τ = turbulence_tensors(atmos, state, diffusive, aux, t)
    turbulence = (ν = ν, D_t = D_t, τ = τ)
    return turbulence
end

@inline function atmos_init_aux!(
    ::SGStkeModel,
    ::AtmosModel,
    aux::Vars,
    geom::LocalGeometry,
)
    aux.sgstke.Δ = lengthscale(geom)
    nothing
end

"""
Set the variable to take the gradient of.
"""
function compute_gradient_argument!(
    ::SGStkeModel,
    atmos::AtmosModel,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    transform.sgstke.e_SGS = state.sgstke.ρe_SGS * (1 / state.ρ)
    ts = recover_thermo_state(atmos, state, aux)
    transform.sgstke.θ_liq_ice = liquid_ice_pottemp(ts)
end

compute_gradient_flux!(::AbstractSGStkeModel, _...) = nothing
"""
Compute the gradient of the variable.
"""
function compute_gradient_flux!(
    ::SGStkeModel,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    diffusive.sgstke.∇e_SGS = ∇transform.sgstke.e_SGS
    diffusive.sgstke.∇θ_liq_ice = ∇transform.sgstke.θ_liq_ice
    diffusive.sgstke.∇u = ∇transform.u
end

# AtmosLESConfigType's default numerical_flux_first_order is RusanovNumericalFlux(),
# and it requires a `flux_first_order!` and `wavespeed` method for the balance law.

function wavespeed_sgstke!(
    ::SGStkeModel,
    wavespeed::Vars,
    nM,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ρinv = 1 / state.ρ
    u = ρinv * state.ρu
    uN = abs(dot(nM, u))

    wavespeed.sgstke.ρe_SGS = uN

    return nothing
end

#####
##### Tendency specification
#####
eq_tends(pv::SGSTKE, m::SGStkeModel, tt::Flux{SecondOrder}) =
    (Diffusion(),)

eq_tends(pv::SGSTKE, m::SGStkeModel, tt::Source) =
    (ShearProduction(), BuoyancyProduction(), Dissipation(),)
#eq_tends(pv::SGSTKE, m::SGStkeModel, tt::Flux{SecondOrder}) = nothing


text:
## wiener filter

@doc doc"""
    argmaxf_lnP(ϕ,                ds::DataSet; kwargs...)
    argmaxf_lnP(ϕ, θ::NamedTuple, ds::DataSet; kwargs...)
    argmaxf_lnP(Lϕ,               ds::DataSet; kwargs...)

Computes either the Wiener filter at fixed $\phi$, or a sample from this slice
along the posterior.

Keyword arguments:

* `which` — `:wf`, `:sample`, or `fluctuation` to compute 1) the Wiener filter,
  i.e. the best-fit of $\mathcal{P}(f\,|\,\phi,d)$, 2) a sample from
  $\mathcal{P}(f\,|\,\phi,d)$, or 3) a sample minus the Wiener filter, i.e. the
  fluctuation on top of the mean.
* `guess` — starting guess for `f` for the conjugate gradient solver
* `conjgrad_kwargs` — Passed to the inner call to [`conjugate_gradient`][@ref]

"""
argmaxf_lnP(ϕ::Field,                ds::DataSet; kwargs...) = argmaxf_lnP(cache(ds.L(ϕ),ds.d), NamedTuple(), ds; kwargs...)
argmaxf_lnP(ϕ::Field, θ::NamedTuple, ds::DataSet; kwargs...) = argmaxf_lnP(cache(ds.L(ϕ),ds.d), θ,            ds; kwargs...)

function argmaxf_lnP(
    Lϕ,
    θ::NamedTuple,
    ds::DataSet;
    which = :wf,
    guess = nothing,
    preconditioner = :diag,
    conjgrad_kwargs = (tol=1e-1,nsteps=500)
)

    @unpack d, Cn, Cn̂, Cf, M, M̂, B, B̂, P = ds(θ)
    D = batchsize(d)

    Δ = d - nonCMB_data_components(θ,ds)
    b = 0
    if (which in (:wf, :sample))
        b += Lϕ'*B'*P'*M'*(Cn\Δ)
    end
    if (which in (:fluctuation, :sample))
        b += Cf\simulate(batch(Cf,D)) + Lϕ'*B'*P'*M'*(Cn\simulate(batch(Cn,D)))
    end

    A_diag  = pinv(Cf) +     B̂' *  M̂'*pinv(Cn̂)*M̂ * B̂
    A_zeroϕ = pinv(Cf) +     B'*P'*M'*pinv(Cn̂)*M*P*B
    A       = pinv(Cf) + Lϕ'*B'*P'*M'*pinv(Cn)*M*P*B*Lϕ

    A_preconditioner = @match preconditioner begin
        :diag  => A_diag
        :zeroϕ => FuncOp(op⁻¹ = (b -> conjugate_gradient(A_diag, A_zeroϕ, b, 0*b, tol=1e-1)))
        _      => error("Unrecognized preconditioner='$preconditioner'")
    end

    conjugate_gradient(A_preconditioner, A, b, (guess==nothing ? 0*b : guess); conjgrad_kwargs...)

end

@doc doc"""
    Σ(ϕ::Field,  ds; [conjgrad_kwargs])
    Σ(Lϕ,        ds; [conjgrad_kwargs])

An operator for the data covariance, Cn + P*M*B*L*Cf*L'*B'*M'*P', which can
applied and inverted. `conjgrad_kwargs` are passed to the underlying call to
`conjugate_gradient`.
"""
Σ(ϕ::Field, ds; kwargs...) = Σ(ds.L(ϕ), ds; kwargs...)
function Σ(Lϕ, ds; conjgrad_kwargs=(tol=1e-1,nsteps=500))
    @unpack d,P,M,B,Cn,Cf,Cn̂,B̂,M̂ = ds
    SymmetricFuncOp(
        op   = x -> (Cn + P*M*B*Lϕ*Cf*Lϕ'*B'*M'*P')*x,
        op⁻¹ = x -> conjugate_gradient((Cn̂ .+ M̂*B̂*Cf*B̂'*M̂'), Σ(Lϕ, ds), x; conjgrad_kwargs...)
    )
end

@doc doc"""

    MAP_joint(ds::DataSet; kwargs...)

Compute the maximum a posteriori (i.e. "MAP") estimate of the joint posterior,
$\mathcal{P}(f,\phi,\theta\,|\,d)$, or compute a quasi-sample.

Keyword arguments:

* `ϕstart` — Starting point of the maximizer *(default:* $\phi=0$*)*
* `nsteps` — The number of iterations for the maximizer.
* `lbfgs_rank` — The maximum rank of the LBFGS approximation to the
   Hessian *(default: 5)*
* `conjgrad_kwargs` — Passed to the inner call to [`conjugate_gradient`][@ref]
* `progress` — whether to show progress bar
* `Nϕ` — Noise to use in the initial approximation to the Hessian. Can also
   give `Nϕ=:qe` to use the quadratic estimate noise *(default:* `:qe`*)*
* `quasi_sample` — `false` *(default)* to compute the MAP, `true` to iterate
   quasi-samples, or an integer to compute a fixed-seed quasi-sample.

Returns a tuple `(f, ϕ, tr)` where `f` is the best-fit (or quasi-sample) field,
`ϕ` is the lensing potential, and `tr` contains info about the run.

"""
MAP_joint(ds::DataSet; kwargs...) = MAP_joint(NamedTuple(), ds; kwargs...)
function MAP_joint(
    θ :: NamedTuple,
    ds :: DataSet;
    nsteps = 20,
    lbfgs_rank = 5,
    Nϕ = :qe,
    ϕstart = nothing,
    progress::Bool = true,
    verbosity = (0,0),
    conjgrad_kwargs = (tol=1e-1,nsteps=500),
    quasi_sample = false,
    preconditioner = :diag,
    aggressive_gc = fieldinfo(ds.d).Nx>=1024 & fieldinfo(ds.d).Ny>=1024,
    αtol = nothing,
    αmax = nothing,
)

    if αtol!=nothing || αmax!=nothing
        @warn "αtol and αmax are deprecated and will be removed in the future. This is now handled automatically."
    end

    dsθ = copy(ds(θ))
    dsθ.G = 1 # MAP estimate is invariant to G so avoid wasted computation

    ϕ = Map(ϕstart==nothing ? zero(diag(ds.Cϕ)) : ϕstart)

    # compute approximate inverse ϕ Hessian used in gradient descent, possibly
    # from quadratic estimate
    if (Nϕ == :qe)
        Nϕ = quadratic_estimate(dsθ).Nϕ/2
    end
    Hϕ⁻¹ = (Nϕ == nothing) ? dsθ.Cϕ : pinv(pinv(dsθ.Cϕ) + pinv(Nϕ))

    tr = []
    argmaxf_lnP_kwargs = (
        which = (quasi_sample==false) ? :wf : :sample,
        preconditioner = preconditioner,
        conjgrad_kwargs = (hist=(:i,:res), progress=false, conjgrad_kwargs...),
    )
    pbar = Progress(nsteps, (progress ? 0 : Inf), "MAP_joint: ")

    f, = argmaxf_lnP((ϕstart==nothing ? 1 : ϕ), θ, dsθ; argmaxf_lnP_kwargs...)
    f°, = mix(f, ϕ, dsθ)
    ϕ, = optimize(
        ϕ -> (
            sum(unbatch(-2lnP(:mix,f°,ϕ,dsθ))),
            gradient(ϕ->-2lnP(:mix,f°,ϕ,dsθ), ϕ)[1]
        ),
        Map(ϕ),
        OptimKit.LBFGS(
            lbfgs_rank;
            maxiter = nsteps,
            verbosity = verbosity[1],
            linesearch = OptimKit.HagerZhangLineSearch(verbosity=verbosity[2])
        ),
        finalize! = function (ϕ,χ²,g,i)
            if isa(quasi_sample,Int)
                seed!(global_rng_for(f),quasi_sample)
            end
            (f, hist) = argmaxf_lnP(
                ϕ, θ, dsθ;
                guess = f,
                argmaxf_lnP_kwargs...
            )
            aggressive_gc && gc()
            f°, = mix(f, ϕ, dsθ)
            g .= gradient(ϕ->-2lnP(:mix,f°,ϕ,dsθ), ϕ)[1]
            χ²s = -2lnP(:mix,f°,ϕ,dsθ)
            χ² = sum(unbatch(χ²s))
            next!(pbar, showvalues=[("step",i), ("χ²",χ²s), ("Ncg",length(hist))])
            push!(tr, @namedtuple(f,ϕ,lnPcur=-χ²/2))
            ϕ, χ², g
        end;
        inner = (_,ξ1,ξ2)->sum(unbatch(dot(ξ1,ξ2))),
        precondition = (_,η)->Map(Hϕ⁻¹*η)
    )

    f, = argmaxf_lnP(ϕ, θ, dsθ; guess=f, argmaxf_lnP_kwargs...)

    f,ϕ,tr

end

OptimKit._scale!(η::Field, β) = η .*= β
OptimKit._add!(η::Field, ξ::Field, β) = η .+= β .* ξ

@doc doc"""

    MAP_marg(ds; kwargs...)

Compute the maximum a posteriori (i.e. "MAP") estimate of the marginal posterior,
$\mathcal{P}(\phi,\theta\,|\,d)$.

"""
MAP_marg(ds::DataSet; kwargs...) = MAP_marg(NamedTuple(), ds; kwargs...)
function MAP_marg(
    θ,
    ds :: DataSet;
    ϕstart = nothing,
    Nϕ = :qe,
    nsteps = 10,
    nsteps_with_meanfield_update = 4,
    conjgrad_kwargs = (tol=1e-1,nsteps=500),
    α = 0.2,
    weights = :unlensed,
    Nsims = 50,
    Nbatch = 1,
    progress::Bool = true,
    aggressive_gc = fieldinfo(ds.d).Nx >=512 & fieldinfo(ds.d).Ny >=512
)

    ds = (@set ds.G = 1)
    @unpack Cf, Cϕ, Cf̃, Cn̂ = ds
    T = eltype(Cf)

    # compute approximate inverse ϕ Hessian used in gradient descent, possibly
    # from quadratic estimate
    if (Nϕ == :qe); Nϕ = quadratic_estimate(ds).Nϕ/2; end
    Hϕ⁻¹ = (Nϕ == nothing) ? Cϕ : pinv(pinv(Cϕ) + pinv(Nϕ))

    ϕ = (ϕstart != nothing) ? ϕstart : ϕ = zero(diag(Cϕ))
    tr = []
    state = nothing
    pbar = Progress(nsteps, (progress ? 0 : Inf), "MAP_marg: ")

    for i=1:nsteps
        aggressive_gc && gc()
        g, state = δlnP_δϕ(
            ϕ, θ, ds,
            use_previous_MF = i>nsteps_with_meanfield_update,
            Nsims=Nsims, Nbatch=Nbatch, weights=weights,
            progress=false, return_state=true, previous_state=state,
            conjgrad_kwargs=conjgrad_kwargs, aggressive_gc=aggressive_gc
        )
        ϕ += T(α) * Hϕ⁻¹ * g
        push!(tr, @dict(i,g,ϕ))
        next!(pbar, showvalues=[
            ("step",i),
            ("Ncg", length(state.gQD.hist)),
            ("Ncg_sims", i<=nsteps_with_meanfield_update ? length(state.gQD_sims[1].hist) : "(MF not updated)"),
            ("α",α)
        ])
    end

    return ϕ, tr

end