text:
import PyPlot
const plt = PyPlot

function eval_LF1d_p( α::Int64, L::Float64, grid_x::Array{Float64}, x)
    Nbasis = size(grid_x)[1]
    N = (Nbasis-1)/2
    f = 0.0
    for l = 1:Nbasis
        k = -N + (l-1)
        f = f + sqrt(1/L/Nbasis)*cos(2*pi*k*(x-grid_x[α])/L)
    end
    return f
end

function gen_T_matrix( L::Float64, grid_x::Array{Float64} )
    Nbasis = size(grid_x)[1]
    N = (Nbasis-1)/2
    # Transformation matrix
    T = Array{Complex128}(Nbasis,Nbasis)
    for α = 1:Nbasis
        for l = 1:Nbasis
            k = -N + (l-1)
            x = grid_x[α]
            T[α,l] = sqrt(1/Nbasis)*exp(2*pi*im*k*x/L)
        end
    end
    return T
end

# Example periodic function
function myfunc(L::Float64, x::Float64)
    ω = 2*pi/L
    f = cos(ω*x)*sin(2*ω*x)
    return f
end

function do_plot_myfunc( NptsPlot::Int64, L::Float64, lim )
    NptsPlot = 200
    x = Array{Float64}(linspace(lim[1], lim[2],NptsPlot));
    y = Array{Float64}(NptsPlot)
    for i = 1:NptsPlot
        y[i] = myfunc(L,x[i])
    end
    plt.clf()
    plt.grid()
    plt.plot(x,y)
    plt.savefig("myfunc.png", dpi=300)
end

function eval_from_ex_coefs( ex_coef::Array{Float64},
        L::Float64, grid_x::Array{Float64}, x::Float64 )
    Nbasis = size(ex_coef)[1]
    f = 0.0
    for i = 1:Nbasis
        f = f + ex_coef[i]*eval_LF1d_p( i, L, grid_x, x )
    end
    return f
end

function eval_FBR( l::Int64, L::Float64, N::Int64, x )
    k = -N + (l-1)
    f = 1/sqrt(L) * exp(im*2*pi*k*x/L)
    return f
end

function FBR_eval_from_ex_coefs( ex_coef, L, x )
    Nbasis = size(ex_coef)[1]
    N = Int64( (Nbasis-1)/2 )
    f = 0.0 + im*0.0
    for i = 1:Nbasis
        f = f + ex_coef[i]*eval_FBR( i, L, N, x )
    end
    return f
end

#function eval_DVR_from_FBR( α, T, x )
#    Nbasis = size(T)[1]
#    f = 0.0 + im*0.0
#    for l = 1:Nbasis
#        f = f + T'[l,α]*eval_FBR()
#    end
#end

function main(L::Float64, N::Int64)
    @printf("L = %f\n", L)
    @printf("N = %d\n", N)
    Nbasis = 2*N + 1
    @printf("Nbasis = %d\n", Nbasis)
    # Create grid
    grid_x = Array{Float64}(Nbasis)
    for α = 1:Nbasis
        grid_x[α] = L/Nbasis*(α-N-1)
    end
    Δ = L/Nbasis

    f1 = eval_LF1d_p( 1, L, grid_x, grid_x[2] )
    f2 = eval_LF1d_p( 2, L, grid_x, grid_x[2] )
    println("\nShould be close to zero: ", f1)
    println("\nShould be close to one: ", f2*sqrt(Δ))

    T = gen_T_matrix(L, grid_x)
    print("\nShould be close to zero: ")
    println(sum(abs.(T'-inv(T)))/Nbasis^2)

    ex_coefs = Array{Float64}(Nbasis)
    for i = 1:Nbasis
        ex_coefs[i] = myfunc(L, grid_x[i])
    end

    x = 2.0
    println("\nShould be close to each other:")
    println( myfunc(L, x) )
    println( eval_from_ex_coefs(ex_coefs*sqrt(Δ), L, grid_x, x) )

    FBR_ex_coefs = inv(T)*ex_coefs*sqrt(Δ)
    println( FBR_eval_from_ex_coefs(FBR_ex_coefs, L, x) )

    # FBR_ex_coefs = ifft(ex_coefs)
    # println( FBR_eval_from_ex_coefs(FBR_ex_coefs, L, x) )

end

main(10.0, 30)


text:
using Test
import ArchGDAL;
const AG = ArchGDAL;

@testset "test_dataset.jl" begin
    @testset "Test methods for raster dataset" begin
        AG.read("data/utmsmall.tif") do dataset
            @testset "Method 1" begin
                AG.copy(
                    dataset,
                    filename = "/vsimem/utmcopy.tif",
                ) do copydataset
                    @test AG.ngcp(copydataset) == 0
                    AG.getband(copydataset, 1) do band
                        @test AG.noverview(band) == 0
                        AG.buildoverviews!(copydataset, Cint[2, 4, 8])
                        @test AG.noverview(band) == 3
                        AG.copywholeraster!(dataset, copydataset)
                        return nothing
                    end
                end
            end
            @testset "Method 2" begin
                copydataset = AG.copy(dataset, filename = "/vsimem/utmcopy.tif")
                @test AG.ngcp(copydataset) == 0
                @test AG.noverview(AG.getband(copydataset, 1)) == 0
                AG.buildoverviews!(copydataset, Cint[2, 4, 8])
                @test AG.noverview(AG.getband(copydataset, 1)) == 3
                AG.copywholeraster!(dataset, copydataset)
            end
            AG.copyfiles("GTiff", "/vsimem/utmcopy2.tif", "/vsimem/utmcopy.tif")
            AG.update("/vsimem/utmcopy2.tif") do copydataset
                @test AG.ngcp(copydataset) == 0
                @test AG.noverview(AG.getband(copydataset, 1)) == 3
                AG.copywholeraster!(
                    dataset,
                    copydataset,
                    options = ["COMPRESS=LZW"],
                )
                return nothing
            end
        end
    end

    @testset "Test methods for vector dataset" begin
        dataset1 = AG.read("data/point.geojson")
        @test AG.nlayer(dataset1) == 1
        layer1 = AG.getlayer(dataset1, 0)
        @test AG.nfeature(layer1) == 4
        AG.getlayer(dataset1, 0) do layer1
            @test AG.nfeature(layer1) == 4
        end
        @test AG.getgeotransform(dataset1) ≈ [0, 1, 0, 0, 0, 1]
        # the following test covers an example of `macro cplerr(code, message)`
        @test_throws ErrorException AG.setproj!(dataset1, "nonsensestring")

        dataset2 = AG.create(AG.getdriver("Memory"))
        @test AG.nlayer(dataset2) == 0
        layer2 = AG.copy(layer1, dataset = dataset2, name = "copy")
        @test AG.nlayer(dataset2) == 1
        @test AG.nfeature(layer2) == 4
        @test AG.getname(layer2) == "copy"

        layer3a = AG.getlayer(dataset2, "copy")
        @test AG.nlayer(dataset2) == 1
        @test AG.nfeature(layer3a) == 4
        @test AG.getname(layer3a) == "copy"
        AG.getlayer(dataset2, "copy") do layer3a
            @test AG.nfeature(layer3a) == 4
            @test AG.getname(layer3a) == "copy"
        end

        layer3b = AG.copy(layer3a)
        @test AG.nlayer(dataset2) == 1 # layer3b is not associated with dataset2
        @test AG.getname(layer3b) == "copy(copy)"

        AG.copy(layer3b, dataset = dataset2)
        @test AG.nlayer(dataset2) == 2
        AG.deletelayer!(dataset2, 1)
        @test AG.nlayer(dataset2) == 1

        # the following test is to cover an example of `macro ogrerr(code, message)`
        @test_throws ErrorException AG.deletelayer!(dataset2, -1)

        dataset4 = AG.create(tempname(), driver = AG.getdriver("KML"))
        @test AG.nlayer(dataset4) == 0
        layer4 = AG.createlayer(
            name = "layer4",
            dataset = dataset4,
            geom = AG.wkbLineString,
        )
        @test AG.nlayer(dataset4) == 1

        AG.create(tempname(), driver = AG.getdriver("KML")) do dataset5
            @test AG.nlayer(dataset5) == 0
            AG.createlayer(
                name = "layer5",
                dataset = dataset5,
                geom = AG.wkbLineString,
            ) do layer5
                @test AG.nfeature(layer5) == 0
            end
            @test AG.nlayer(dataset5) == 1
        end

        layer5 = AG.createlayer()
        @test AG.getname(layer5) == ""
        @test AG.nfeature(layer5) == 0
        @test AG.nfield(layer5) == 0
        @test AG.ngeom(layer5) == 1

        AG.create(AG.getdriver("Memory")) do dataset6
            for i in 1:20
                AG.createlayer(name = "layer$(i - 1)", dataset = dataset6)
            end
            @test AG.nlayer(dataset6) == 20
            @test sprint(print, dataset6) == """
            GDAL Dataset (Driver: Memory/Memory)
            File(s):

            Number of feature layers: 20
              Layer 0: layer0 (wkbUnknown)
              Layer 1: layer1 (wkbUnknown)
              Layer 2: layer2 (wkbUnknown)
              Layer 3: layer3 (wkbUnknown)
              Layer 4: layer4 (wkbUnknown)
              Remaining layers:
                layer5, layer6, layer7, layer8, layer9,
                layer10, layer11, layer12, layer13, layer14,
                layer15, layer16, layer17, layer18, layer19, """
        end
    end
end

ext: jl


text:
#=
Challenge #261 [Easy] verifying 3x3 magic squares

Description

A 3x3 magic square is a 3x3 grid of the numbers 1-9 such that each row, column,
and major diagonal adds up to 15. Here's an example:

8 1 6
3 5 7
4 9 2

The major diagonals in this example are 8 + 5 + 2 and 6 + 5 + 4. (Magic squares
have appeared here on r/dailyprogrammer before, in #65 [Difficult] in 2012.)
Write a function that, given a grid containing the numbers 1-9, determines whether
it's a magic square. Use whatever format you want for the grid, such as a
2-dimensional array, or a 1-dimensional array of length 9, or a function that
takes 9 arguments. You do not need to parse the grid from the program's input,
but you can if you want to. You don't need to check that each of the 9 numbers
appears in the grid: assume this to be true.

Example inputs/outputs

[8, 1, 6, 3, 5, 7, 4, 9, 2] => true
[2, 7, 6, 9, 5, 1, 4, 3, 8] => true
[3, 5, 7, 8, 1, 6, 4, 9, 2] => false
[8, 1, 6, 7, 5, 3, 4, 9, 2] => false

Optional bonus 1

Verify magic squares of any size, not just 3x3.

Optional bonus 2

Write another function that takes a grid whose bottom row is missing, so it only
has the first 2 rows (6 values). This function should return true if it's possible
to fill in the bottom row to make a magic square. You may assume that the numbers
given are all within the range 1-9 and no number is repeated. Examples:

[8, 1, 6, 3, 5, 7] => true
[3, 5, 7, 8, 1, 6] => false

Hint: it's okay for this function to call your function from the main challenge.
This bonus can also be combined with optional bonus 1. (i.e. verify larger magic
squares that are missing their bottom row.)

( source: https://redd.it/4dccix )
=#

read_matrix( data::DataType=Any, input::IO=STDIN, spacing::Char=' ' ) = readdlm(
input, spacing, data, header=false, skipblanks=true, ignore_invalid_chars=true,
quotes=false, comments=false )

function pseudo_magic_square( grid::Matrix{ Int } )
    local value::Bool = loop::Bool = false
    local missing::Int = tmp::Int = 0
    local last::Array{ Int, 1 } = Int[]

    for i in 1:size( grid, 2 )
        column = sum( grid[ :, i ] )
        missing = 15 - column
        push!( last, missing )

        if 9 >= missing
            for j in 1:size( grid, 1 ), k in 1:size( grid, 2 )
                if grid[ j, k ] == missing
                    loop = true
                    break
                end
            end
        else
            loop = true
        end

        if true == loop
            break
        end
    end

    if false == loop
        for i in 1:2
            tmp += grid[ i,i ]
        end
        tmp += last[ 3 ]
        if 15 == tmp
            value = true
        end
    end

    return value
end

function verify_magic_square( grid::Matrix{ Int } )
    local value::Bool = false
    local line::Integer = column::Integer = 0

    if 15 == trace( grid )  #   sums the value of the principal diagonal
        for i in 1:size( grid, 1 )
            line = sum( grid[ i, : ] )
            column = sum( grid[ :, i ] )
            if line != 15 || column != 15
                break
            elseif size( grid, 1 ) == i   #   if is the last verification and it's valid
                value = true
            end
        end
    end

    return value
end

#=
    I've modified the input format. So, that only be possibles input like:

        8 1 6                       8 1 6
        3 5 7           or          3 5 7
        4 9 2
=#
function main( )
    local input::Matrix{ Int } = read_matrix( Int )

    println( @time ( size( input, 1 ) == size( input, 2 ) ) ?
                     verify_magic_square( input ):
                     pseudo_magic_square( input ) )
end

main( )

max_stars_repo_name: farmfirst/dailyprogrammer

max_issues_repo_name: farmfirst/dailyprogrammer

max_forks_repo_name: farmfirst/dailyprogrammer


text:
__precompile__(true)
module LevelSetEvaluation

using Reexport
@reexport using AutoRisk
# manually include an AutoRisk config script
include(joinpath(Pkg.dir("AutoRisk"), "scripts", "collection",
    "heuristic_dataset_config.jl"))

@reexport using Coord

include("level_set_evaluator.jl")
include("auto_evaluator.jl")
include("uav_evaluator.jl")

end # module


text:
struct MagneticPartitioner end

function partition_stripe(A::SparseMatrixCSC, k, ::MagneticPartitioner, Π; kwargs...)
    @inbounds begin
        (m, n) = size(A)
        Π = convert(MapPartition, Π)
        asg = undefs(Int, n)
        for j = 1:n
            if A.colptr[j] == A.colptr[j + 1]
                asg[j] = rand(1:k)
                #asg[j] = k
            else
                i = A.rowval[rand(A.colptr[j]:A.colptr[j + 1] - 1)]
                #i = A.rowval[A.colptr[j]]
                asg[j] = Π.asg[i]
            end
        end
        return MapPartition(k, asg)
    end
end

struct GreedyBottleneckPartitioner{Mdl}
    mdl::Mdl
end

function partition_stripe(A::SparseMatrixCSC{Tv, Ti}, K, method::GreedyBottleneckPartitioner{Mdl}, Π; adj_A = nothing, kwargs...) where {Tv, Ti, Mdl <: AbstractSecondaryConnectivityModel}
    @inbounds begin
        (m, n) = size(A)
        if adj_A === nothing
            adj_A = adjointpattern(A)
        end
        n_vertices = undefs(Ti,K)
        n_pins = undefs(Ti,K)
        n_local_nets = undefs(Ti,K)
        n_remote_nets = undefs(Ti,K)
        cst = undefs(cost_type(method.mdl),K)

        Π_dmn = convert(DomainPartition, Π)
        Π_map = convert(MapPartition, Π)

        m, n = size(A)
        hst = zeros(Ti, n)
        for k = 1:K
            s = Π_dmn.spl[k]
            s′ = Π_dmn.spl[k + 1]
            n_vertices_k = s′ - s
            n_pins_k = 0
            n_local_nets_k = 0
            n_remote_nets_k = 0
            for _s = s:(s′ - 1)
                _i = Π_dmn.prm[_s]
                q = adj_A.colptr[_i]
                q′ = adj_A.colptr[_i + 1]
                n_pins_k += q′ - q
                for _q = q : q′ - 1
                    j = adj_A.rowval[_q]
                    if hst[j] < s
                        n_remote_nets_k += 1
                    end
                    hst[j] = s
                end
            end
            n_vertices[k] = n_vertices_k
            n_pins[k] = n_pins_k
            n_local_nets[k] = n_local_nets_k
            n_remote_nets[k] = n_remote_nets_k
            cst[k] = method.mdl(n_vertices_k, n_pins_k, n_local_nets_k, n_remote_nets_k, k)
        end

        asg = ones(Ti, n)

        for j = randperm(n)
            best_c = -Inf
            best_k = 0
            for q = A.colptr[j] : A.colptr[j + 1] - 1
                i = A.rowval[q]
                k = Π_map.asg[i]
                c = cst[k]
                if c > best_c
                    best_c = c
                    best_k = k
                end
            end
            if best_k != 0
                asg[j] = best_k
                n_local_nets[best_k] += 1
                n_remote_nets[best_k] -= 1
                cst[best_k] = method.mdl(n_vertices[best_k], n_pins[best_k], n_local_nets[best_k], n_remote_nets[best_k], best_k)
            end
        end

        return MapPartition(K, asg)
    end
end

function partition_stripe(A::SparseMatrixCSC{Tv, Ti}, K, method::GreedyBottleneckPartitioner{Mdl}, Π::SplitPartition; adj_A = nothing, adj_net=nothing, kwargs...) where {Tv, Ti, Mdl <: AbstractSecondaryConnectivityModel}
    @inbounds begin
        (m, n) = size(A)
        if adj_A === nothing
            adj_A = adjointpattern(A)
        end
        n_vertices = undefs(Ti,K)
        n_pins = undefs(Ti,K)
        n_local_nets = undefs(Ti,K)
        n_remote_nets = undefs(Ti,K)
        cst = undefs(cost_type(method.mdl),K)

        m, n = size(A)

        if adj_net === nothing
            hst = zeros(Ti, n)
            for k = 1:K
                i = Π.spl[k]
                i′ = Π.spl[k + 1]
                n_vertices_k = i′ - i
                n_pins_k = adj_A.colptr[i′] - adj_A.colptr[i]
                n_local_nets_k = 0
                n_remote_nets_k = 0
                for _i = i:(i′ - 1)
                    q = adj_A.colptr[_i]
                    q′ = adj_A.colptr[_i + 1]
                    for _q = q : q′ - 1
                        j = adj_A.rowval[_q]
                        if hst[j] < i
                            n_remote_nets_k += 1
                        end
                        hst[j] = i
                    end
                end
                n_vertices[k] = n_vertices_k
                n_pins[k] = n_pins_k
                n_local_nets[k] = n_local_nets_k
                n_remote_nets[k] = n_remote_nets_k
                cst[k] = method.mdl(n_vertices_k, n_pins_k, n_local_nets_k, n_remote_nets_k, k)
            end
        else
            for k = 1:K
                i = Π.spl[k]
                i′ = Π.spl[k + 1]
                n_vertices_k = i′ - i
                n_pins_k = adj_A.colptr[i′] - adj_A.colptr[i]
                n_local_nets_k = 0
                n_remote_nets_k = adj_net[i, i′]
                n_vertices[k] = n_vertices_k
                n_pins[k] = n_pins_k
                n_local_nets[k] = n_local_nets_k
                n_remote_nets[k] = n_remote_nets_k
                cst[k] = method.mdl(n_vertices_k, n_pins_k, n_local_nets_k, n_remote_nets_k, k)
            end
        end

        Π = convert(MapPartition, Π)
        asg = ones(Ti, n)

        for j = randperm(n)
            best_c = typemin(cost_type(method.mdl))
            best_k = 0
            for q = A.colptr[j] : A.colptr[j + 1] - 1
                i = A.rowval[q]
                k = Π.asg[i]
                c = cst[k]
                if c > best_c
                    best_c = c
                    best_k = k
                end
            end
            if best_k != 0
                asg[j] = best_k
                n_local_nets[best_k] += 1
                n_remote_nets[best_k] -= 1
                cst[best_k] = method.mdl(n_vertices[best_k], n_pins[best_k], n_local_nets[best_k], n_remote_nets[best_k], best_k)
            end
        end

        return MapPartition(K, asg)
    end
end


text:
mutable struct Instruction
    op::String
    arg::Int
end

function instructionParser(line)
    parts = split(line, " ")
    op = parts[1]
    arg = parse(Int, parts[2])
    return Instruction(op, arg)
end

function execute(program)
    # print("executing ",program)
    accumulator = 0
    loc = length(program)
    executedInstructions::Vector = zeros(Int,loc)
    pointer = 1
    while pointer <= length(program)
        if executedInstructions[pointer] == 0
            instruction = program[pointer]
            if instruction.op == "nop"
                executedInstructions[pointer] += 1
                pointer += 1
            elseif instruction.op == "acc"
                accumulator += instruction.arg
                executedInstructions[pointer] += 1
                pointer += 1
            elseif instruction.op == "jmp"
                executedInstructions[pointer] += 1
                pointer += instruction.arg
            end
        else
            println(" result=", 0)
            return 0
        end
    end
    println(" result=", accumulator)
    return accumulator
end

function changeInstruction(program, pointer)
    instruction = program[pointer]
    currOp = instruction.op
    if currOp == "nop"
        instruction.op = "jmp"
    else
        instruction.op = "nop"
    end
    program[pointer] = instruction
    #println("Changed instruction to ", instruction.op)
    return program
end

function fix(program)
    pointer = 1
    result = 0
    while result == 0
        # println("pointer=", pointer)
        result = execute(program)
        if result == 0
            changeInstruction(program, pointer)
            result = execute(program)
            if (result == 0)
                changeInstruction(program, pointer)
                pointer += 1
                while program[pointer].op == "acc"
                    if pointer < length(program)
                        pointer += 1
                        # println("skipping acc")
                    end
                end
            end
        end
    end
    println("success",result)
end

program = []
open("data.txt") do file
    for line in eachline(file)
        instruction = instructionParser(line)
        push!(program, instruction)
    end
end

fix(program)


text:
#  Copyright 2022 California Institute of Technology
#
#  Licensed under the Apache License, Version 2.0 (the "License");
#  you may not use this file except in compliance with the License.
#  You may obtain a copy of the License at
#
#      http://www.apache.org/licenses/LICENSE-2.0
#
#  Unless required by applicable law or agreed to in writing, software
#  distributed under the License is distributed on an "AS IS" BASIS,
#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#  See the License for the specific language governing permissions and
#  limitations under the License.
#
# Author: Philip G. Brodrick, philip.g.brodrick@jpl.nasa.gov

using DataFrames
using CSV
using Interpolations
using Logging
using Statistics
using Plots
using ModularIndices
using NMF
using MultivariateStats
using Clustering

function nanargmax(input::Array)
    x = copy(input)
    x[isnan.(x)] .= -Inf;
    return argmax(x);
end

function nanargmin(input::Array)
    x = copy(input)
    x[isnan.(x)] .= Inf;
    return argmin(x);
end

function read_envi_wavelengths(filename::String)
    header_name = splitext(filename)[1] * ".hdr"
    header = readlines(header_name)
    found = false
    for line in header
        if occursin("wavelength = {", line) || occursin("wavelength= {", line) ||  occursin("wavelength={", line)
            header = line
            found = true
            break
        end
    end

    if !found
        @error "No wavelength found in " * header_name
        return nothing
    end

    wavelengths = [parse(Float64, strip(x)) for x in split(split(split(header, "{")[2], "}")[1],",")]
    return wavelengths
end

function get_good_bands_mask(wavelengths::Array{Float64}, wavelength_pairs)
    good_bands = ones(Bool, length(wavelengths))

    for wvp in wavelength_pairs
        wavelength_diff = wavelengths .- wvp[1]
        wavelength_diff[wavelength_diff .< 0] .= maximum(filter(!isnan, wavelength_diff))
        lower_index = nanargmin(wavelength_diff)

        wavelength_diff = wvp[2] .- wavelengths
        wavelength_diff[wavelength_diff .< 0] .= maximum(filter(!isnan, wavelength_diff))
        upper_index = nanargmin(wavelength_diff)
        good_bands[lower_index:upper_index] .= false
    end

    return good_bands
end

mutable struct SpectralLibrary
    file_name::String
    class_header_name::String
    spectral_starting_column::Int64
    truncate_end_columns::Int64
    class_valid_keys
    scale_factor::Float64
    wavelength_regions_ignore
    SpectralLibrary(file_name::String, class_header_name::String, spectral_starting_column::Int64 = 2, truncate_end_columns::Int64 = 0, class_valid_keys = nothing,
                    scale_factor = 1.0, wavelength_regions_ignore= [[0,440],[1310,1490],[1770,2050],[2440,2880]]) =
                    new(file_name, class_header_name, spectral_starting_column, truncate_end_columns, class_valid_keys, scale_factor, wavelength_regions_ignore)

    spectra
    classes
    good_bands
    wavelengths
end

function load_data!(library::SpectralLibrary)

    df = DataFrame(CSV.File(library.file_name))

    try
        library.spectra = Matrix(df[:,library.spectral_starting_column:end-library.truncate_end_columns])
    catch e
        throw(ArgumentError("Could not read spectra from endmember library.  Try adjusting spectral_starting_column or truncate_end_columns, or reworking the library."))
    end

    try
        library.classes = convert(Array{AbstractString}, df[!,library.class_header_name])
    catch e
        throw(ArgumentError("Could not read classes from endmember library using class key "*library.class_header_name *
                            ".  Try adjusting class_header_name, or reworking the library."))
    end

    try
        library.wavelengths = parse.(Float64, names(df[:,library.spectral_starting_column:end-library.truncate_end_columns]))
    catch e
        throw(ArgumentError("Could not read wavelengths from endmember library.  Try adjusting class_header_name, or reworking the library."))
    end

    wl_order = sortperm(library.wavelengths)
    library.spectra = library.spectra[:,wl_order]
    library.wavelengths = library.wavelengths[wl_order]

    if isnothing(library.class_valid_keys)
        library.class_valid_keys = unique(library.classes)
    end

    good_bands = get_good_bands_mask(library.wavelengths, library.wavelength_regions_ignore)
    library.good_bands = good_bands

    return library
end

function save_data(library::SpectralLibrary, output_filename::String; class_label_name::String = "Label")

    df = DataFrame()
    insertcols!(df, 1, class_label_name => library.classes)
    for (_wv, wv) in enumerate(library.wavelengths)
        insertcols!(df, 1 + _wv, string(wv) => library.spectra[:,_wv])
    end
    CSV.write(output_filename, df)

end

function filter_by_class!(library::SpectralLibrary)
    if isnothing(library.class_valid_keys)
        @info "No class valid keys provided, no filtering occuring"
        return
    end

    valid_classes = zeros(Bool, size(library.spectra)[1])
    for cla in library.class_valid_keys
        valid_classes[library.classes .== cla] .= true
    end

    library.spectra = library.spectra[valid_classes,:]
    library.classes = library.classes[valid_classes]
end

function remove_wavelength_region_inplace!(library::SpectralLibrary, set_as_nans::Bool=false)
    if set_as_nans
        library.spectra[:,.!library.good_bands] .= NaN
        library.wavelengths[.!library.good_bands] .= NaN
    else
        library.spectra = library.spectra[:, library.good_bands]
        library.wavelengths = library.wavelengths[library.good_bands]
        library.good_bands = library.good_bands[library.good_bands]
    end
end

function interpolate_library_to_new_wavelengths!(library::SpectralLibrary, new_wavelengths::Array{Float64})
    old_spectra = copy(library.spectra)

    library.spectra = zeros((size(library.spectra)[1], length(new_wavelengths)))
    for _s in 1:size(old_spectra)[1]
        fit = LinearInterpolation(library.wavelengths, old_spectra[_s,:], extrapolation_bc=Flat());
        library.spectra[_s,:] = fit(new_wavelengths)
    end
    library.wavelengths = new_wavelengths

    good_bands = get_good_bands_mask(library.wavelengths, library.wavelength_regions_ignore)
    library.good_bands = good_bands
end

function scale_library!(library::SpectralLibrary, scaling_factor=nothing)
    if isnothing(scaling_factor)
        library.spectra =  library.spectra ./ library.scale_factor
    else
        library.spectra /= scaling_factor
    end
end

function brightness_normalize!(library::SpectralLibrary)
    library.spectra = library.spectra ./ sqrt.(mean(library.spectra[:,library.good_bands].^2, dims=2))
end

function split_library(library::SpectralLibrary, split_fraction::Float64)
    perm = randperm((library.spectra)[1])

    split_1 = perm[1:round(split_fraction * length(perm))]
    split_2 = perm[round(split_fraction * length(perm)):end]

    output_library_1 = deepcopy(library)
    output_library_2 = deepcopy(library)

    output_library_1.spectra = output_library_1.spectra[split_1,:]
    output_library_2.spectra = output_library_2.spectra[split_2,:]

    output_library_1.classes = output_library_1.classes[split_1]
    output_library_2.classes = output_library_2.classes[split_2]

    return output_library_1, output_library_2

end

function reduce_endmembers_nmf!(library::SpectralLibrary, max_endmembers_per_class::Int64)

    reduced_library = []
    new_classes = copy(library.classes)[1:max_endmembers_per_class*size(library.class_valid_keys)[1]]
    for (_c, cla) in enumerate(library.class_valid_keys)
        library_subset = library.spectra[library.classes .== cla, library.good_bands]

        out_spectra = zeros(max_endmembers_per_class, size(library.spectra)[2])
        out_spectra[:,:] .= NaN

        r = nnmf(library_subset, max_endmembers_per_class; maxiter=500, tol=1.0e-2)
        out_spectra[:,library.good_bands] = r.H

        push!(reduced_library,out_spectra)
        new_classes[(_c-1) * max_endmembers_per_class + 1:_c * max_endmembers_per_class] .= cla
    end
    library.spectra = cat(reduced_library...,dims=1)
    library.classes = new_classes
end

function reduce_endmembers_pca!(library::SpectralLibrary, max_endmembers_per_class::Int64)

    reduced_library = []
    new_classes = copy(library.classes)[1:max_endmembers_per_class*size(library.class_valid_keys)[1]]
    for (_c, cla) in enumerate(library.class_valid_keys)
        library_subset = library.spectra[library.classes .== cla, library.good_bands]

        out_spectra = zeros(max_endmembers_per_class, size(library.spectra)[2])
        out_spectra[:,:] .= NaN

        model = fit(PCA, library_subset', maxoutdim=max_endmembers_per_class, pratio=1)
        out_spectra[:,library.good_bands] = projection(model)'

        push!(reduced_library,out_spectra)
        new_classes[(_c-1) * max_endmembers_per_class + 1:_c * max_endmembers_per_class] .= cla
    end
    library.spectra = cat(reduced_library...,dims=1)
    library.classes = new_classes
end

function reduce_endmembers_kmeans!(library::SpectralLibrary, max_endmembers_per_class::Int64)

    reduced_library = []
    new_classes = copy(library.classes)[1:max_endmembers_per_class*size(library.class_valid_keys)[1]]
    for (_c, cla) in enumerate(library.class_valid_keys)
        library_subset = library.spectra[library.classes .== cla, library.good_bands]

        out_spectra = zeros(max_endmembers_per_class, size(library.spectra)[2])
        out_spectra[:,:] .= NaN

        model = kmeans(library_subset', max_endmembers_per_class, maxiter=1000)
        out_spectra[:,library.good_bands] = model.centers'

        push!(reduced_library,out_spectra)
        new_classes[(_c-1) * max_endmembers_per_class + 1:_c * max_endmembers_per_class] .= cla
    end
    library.spectra = cat(reduced_library...,dims=1)
    library.classes = new_classes
end

max_stars_repo_name: fotxoa-geo/SpectralUnmixing

max_issues_repo_name: fotxoa-geo/SpectralUnmixing

max_forks_repo_name: fotxoa-geo/SpectralUnmixing


text:
# elastic model 3d
function write_mor(fid::IOStream, model::elmod3d)
    write(fid,model.wf.vx)
    write(fid,model.wf.vy)
    write(fid,model.wf.vz)
    write(fid,model.wf.txx)
    write(fid,model.wf.tyy)
    write(fid,model.wf.tzz)
    write(fid,model.wf.txy)
    write(fid,model.wf.txz)
    write(fid,model.wf.tyz)
    write(fid,model.pml.PVxBTxx)
    write(fid,model.pml.PVxBTxy)
    write(fid,model.pml.PVxBTxz)
    write(fid,model.pml.PVyBTxy)
    write(fid,model.pml.PVyBTyy)
    write(fid,model.pml.PVyBTyz)
    write(fid,model.pml.PVzBTxz)
    write(fid,model.pml.PVzBTyz)
    write(fid,model.pml.PVzBTzz)
    write(fid,model.pml.PTxxBVx)
    write(fid,model.pml.PTxxBVy)
    write(fid,model.pml.PTxxBVz)
    write(fid,model.pml.PTyyBVx)
    write(fid,model.pml.PTyyBVy)
    write(fid,model.pml.PTyyBVz)
    write(fid,model.pml.PTzzBVx)
    write(fid,model.pml.PTzzBVy)
    write(fid,model.pml.PTzzBVz)
    write(fid,model.pml.PTxyBVx)
    write(fid,model.pml.PTxyBVy)
    write(fid,model.pml.PTxzBVx)
    write(fid,model.pml.PTxzBVz)
    write(fid,model.pml.PTyzBVy)
    write(fid,model.pml.PTyzBVz)
end

# acoustic model 3d
function write_mor(fid::IOStream, model::acmod3d)
    write(fid,model.wf.vx)
    write(fid,model.wf.vy)
    write(fid,model.wf.vz)
    write(fid,model.wf.tpp)
    write(fid,model.pml.PVxBTpp)
    write(fid,model.pml.PVyBTpp)
    write(fid,model.pml.PVzBTpp)
    write(fid,model.pml.PTppBVx)
    write(fid,model.pml.PTppBVy)
    write(fid,model.pml.PTppBVz)
end

# elastic model 2d
function write_mor(fid::IOStream, model::elmod2d)
    write(fid,model.wf.vx)
    write(fid,model.wf.vz)
    write(fid,model.wf.txx)
    write(fid,model.wf.tzz)
    write(fid,model.wf.txz)
    write(fid,model.pml.PVxBTxx)
    write(fid,model.pml.PVxBTxz)
    write(fid,model.pml.PVzBTxz)
    write(fid,model.pml.PVzBTzz)
    write(fid,model.pml.PTxxBVx)
    write(fid,model.pml.PTxxBVz)
    write(fid,model.pml.PTzzBVx)
    write(fid,model.pml.PTzzBVz)
    write(fid,model.pml.PTxzBVx)
    write(fid,model.pml.PTxzBVz)
end

# acoustic model 2d
function write_mor(fid::IOStream, model::acmod2d)
    write(fid,model.wf.vx)
    write(fid,model.wf.vz)
    write(fid,model.wf.tpp)
    write(fid,model.pml.PVxBTpp)
    write(fid,model.pml.PVzBTpp)
    write(fid,model.pml.PTppBVx)
    write(fid,model.pml.PTppBVz)
end


text:
include("prop_logic_theory.jl")
include("prover.jl")

ex = rewrite(:(((p => q) ∧ (r => s) ∧ (p ∨ r)) => (q ∨ s)), impl)
@test prove(t, ex, 3, 7)

@test @areequal t true ((¬p == p) == false)
@test @areequal t true ((¬p == ¬p) == true)
@test @areequal t true ((¬p ∨ ¬p) == ¬p) (¬p ∨ p) ¬(¬p ∧ p)
@test @areequal t true ((p => (p ∨ p)))
@test @areequal t true ((p => (p ∨ p)) == ((¬(p) ∧ q) => q)) == true

# Frege's theorem
@test @areequal t true (p => (q => r)) => ((p => q) => (p => r))

# Demorgan's
@test @areequal t true (¬(p ∨ q) == (¬p ∧ ¬q))

# Consensus theorem
@test @areequal t ((x ∧ y) ∨ (¬x ∧ z) ∨ (y ∧ z))   ((x ∧ y) ∨ (¬x ∧ z))

# @timev areequal(t, :((x ∧ y) ∨ (¬x ∧ z) ∨ (y ∧ z)), :((x ∧ y) ∨ (¬x ∧ z)))
# @timev areequal(t, :((x ∧ y) ∨ (¬x ∧ z) ∨ (y ∧ z)), :((x ∧ y) ∨ (¬babo ∧ z)))
# @timev areequalmagic(t, :((x ∧ y) ∨ (¬x ∧ z) ∨ (y ∧ z)),   :((x ∧ y) ∨ (¬x ∧ z)))
# @timev areequalmagic(t, :((x ∧ y) ∨ (¬x ∧ z) ∨ (y ∧ z)),   :((babo ∧ y) ∨ (¬x ∧ z)))

# TODO proof strategies?
# Constructive Dilemma

# @test @areequal (t ∪ [@rule :p => true]) true (((p => q) ∧ (r => s)) ∧ (p ∨ r)) => (q ∨ s)

# @test areequal(t, true, :(¬(((¬p ∨ q) ∧ (¬r ∨ s)) ∧ (p ∨ r)) ∨ (q ∨ s)))

# using Metatheory.EGraphs.Schedulers
# Metatheory.options.verbose = true
# Metatheory.options.printiter = true
# ex = rewrite(:(((p => q) ∧ (r => s) ∧ (p ∨ r)) => (q ∨ s)), impl)
# g = EGraph(ex)
# params = SaturationParams(timeout=10, eclasslimit=5000, scheduler=SimpleScheduler)#, scheduler=ScoredScheduler)
# @profview saturate!(g, t, params)

# ex = :((p ∧ q ∨ (¬q ∧ (p => r) ∧ q)) => (s ∧ q))
# prove(t, ex, 3)

# extran = addanalysis!(g, ExtractionAnalysis, astsize)
# ex = extract!(g, extran)
# println(ex)
# #
# #
# ex = rewrite(:(((p => q) ∧ (r => s) ∧ (p ∨ r)) => (q ∨ s)), impl)
# g = EGraph(ex)
# @timev areequal(g, t, ex, true; timeout=10, eclasslimit=5000, scheduler=ScoredScheduler, schedulerparams=(32, 1))
#

# @profiler saturate!(g, t; timeout=8, eclasslimit=5000)
# exit(0)

# ex = rewrite(:(((p => p) ∧ (r => z) ∧ (p ∨ r)) => (q ∨ s)), impl)
# @test false == prove(t, ex, 4)

# g = EGraph(:(((p => q) ∧ (r => s) ∧ (p ∨ r)) => (q ∨ s)))
# @time saturate!(g, t; timeout=30, eclasslimit=Inf)
#
# in_same_set(g.uf, g.root, addexpr!(g, true).id) |> println
#
# struct LogicAnalysis <: AbstractAnalysis
#     egraph::EGraph
#     logic_egraph::EGraph
# end


text:
using .ReverseDiff

(ev::Eval)(A::ReverseDiff.TrackedArray, args...) = ReverseDiff.track(ev, A, args...)
(ev::Eval)(A, B::ReverseDiff.TrackedArray, args...) = ReverseDiff.track(ev, A, B, args...)
(ev::Eval)(A::ReverseDiff.TrackedArray, B::ReverseDiff.TrackedArray, args...) = ReverseDiff.track(ev, A, B, args...)

ReverseDiff.@grad function (ev::Eval)(args...)
    Z = ev.fwd(ReverseDiff.value.(args)...)
    Z, Δ -> begin
        isnothing(ev.rev) && error("no gradient definition here!")
        ev.rev(Δ, ReverseDiff.value.(args)...)
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "AMCL"
version = v"2.0.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/apache/incubator-milagro-crypto-c/archive/$version.tar.gz", "48fc7e0d7ad5edbfe8a779dc56f28813c25798ddfd60f1d94ccf1203eacd5645"),
]

# Bash recipe for building across all platforms
script = raw"""
env | sort
cd "$WORKSPACE/srcdir"/incubator-milagro-crypto-c-*
mkdir build
cd build
cmake -G "Ninja" -DCMAKE_INSTALL_PREFIX="$prefix" -DCMAKE_TOOLCHAIN_FILE="${CMAKE_TARGET_TOOLCHAIN}" -DBUILD_BENCHMARKS=off -DBUILD_DOCS=off -DBUILD_EXAMPLES=off -DBUILD_PYTHON=off -DBUILD_TESTING=off -DWORD_SIZE=$nbits ..
ninja
ninja install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libamcl_bls_BLS24", :libamcl_bls_BLS24),
    LibraryProduct("libamcl_bls_BLS381", :libamcl_bls_BLS381),
    LibraryProduct("libamcl_bls_BLS383", :libamcl_bls_BLS383),
    LibraryProduct("libamcl_bls_BLS461", :libamcl_bls_BLS461),
    LibraryProduct("libamcl_bls_BLS48", :libamcl_bls_BLS48),
    LibraryProduct("libamcl_bls_BN254", :libamcl_bls_BN254),
    LibraryProduct("libamcl_bls_BN254CX", :libamcl_bls_BN254CX),
    LibraryProduct("libamcl_bls_FP256BN", :libamcl_bls_FP256BN),
    LibraryProduct("libamcl_bls_FP512BN", :libamcl_bls_FP512BN),
    LibraryProduct("libamcl_core", :libamcl_core),
    LibraryProduct("libamcl_curve_ANSSI", :libamcl_curve_ANSSI),
    LibraryProduct("libamcl_curve_BLS24", :libamcl_curve_BLS24),
    LibraryProduct("libamcl_curve_BLS381", :libamcl_curve_BLS381),
    LibraryProduct("libamcl_curve_BLS383", :libamcl_curve_BLS383),
    LibraryProduct("libamcl_curve_BLS461", :libamcl_curve_BLS461),
    LibraryProduct("libamcl_curve_BLS48", :libamcl_curve_BLS48),
    LibraryProduct("libamcl_curve_BN254", :libamcl_curve_BN254),
    LibraryProduct("libamcl_curve_BN254CX", :libamcl_curve_BN254CX),
    LibraryProduct("libamcl_curve_BRAINPOOL", :libamcl_curve_BRAINPOOL),
    LibraryProduct("libamcl_curve_C25519", :libamcl_curve_C25519),
    LibraryProduct("libamcl_curve_C41417", :libamcl_curve_C41417),
    LibraryProduct("libamcl_curve_ED25519", :libamcl_curve_ED25519),
    LibraryProduct("libamcl_curve_FP256BN", :libamcl_curve_FP256BN),
    LibraryProduct("libamcl_curve_FP512BN", :libamcl_curve_FP512BN),
    LibraryProduct("libamcl_curve_GOLDILOCKS", :libamcl_curve_GOLDILOCKS),
    LibraryProduct("libamcl_curve_HIFIVE", :libamcl_curve_HIFIVE),
    LibraryProduct("libamcl_curve_NIST256", :libamcl_curve_NIST256),
    LibraryProduct("libamcl_curve_NIST384", :libamcl_curve_NIST384),
    LibraryProduct("libamcl_curve_NIST521", :libamcl_curve_NIST521),
    LibraryProduct("libamcl_curve_NUMS256E", :libamcl_curve_NUMS256E),
    LibraryProduct("libamcl_curve_NUMS256W", :libamcl_curve_NUMS256W),
    LibraryProduct("libamcl_curve_NUMS384E", :libamcl_curve_NUMS384E),
    LibraryProduct("libamcl_curve_NUMS384W", :libamcl_curve_NUMS384W),
    LibraryProduct("libamcl_curve_NUMS512E", :libamcl_curve_NUMS512E),
    LibraryProduct("libamcl_curve_NUMS512W", :libamcl_curve_NUMS512W),
    LibraryProduct("libamcl_curve_SECP256K1", :libamcl_curve_SECP256K1),
    LibraryProduct("libamcl_mpin_BLS24", :libamcl_mpin_BLS24),
    LibraryProduct("libamcl_mpin_BLS381", :libamcl_mpin_BLS381),
    LibraryProduct("libamcl_mpin_BLS383", :libamcl_mpin_BLS383),
    LibraryProduct("libamcl_mpin_BLS461", :libamcl_mpin_BLS461),
    LibraryProduct("libamcl_mpin_BLS48", :libamcl_mpin_BLS48),
    LibraryProduct("libamcl_mpin_BN254", :libamcl_mpin_BN254),
    LibraryProduct("libamcl_mpin_BN254CX", :libamcl_mpin_BN254CX),
    LibraryProduct("libamcl_mpin_FP256BN", :libamcl_mpin_FP256BN),
    LibraryProduct("libamcl_mpin_FP512BN", :libamcl_mpin_FP512BN),
    LibraryProduct("libamcl_pairing_BLS24", :libamcl_pairing_BLS24),
    LibraryProduct("libamcl_pairing_BLS381", :libamcl_pairing_BLS381),
    LibraryProduct("libamcl_pairing_BLS383", :libamcl_pairing_BLS383),
    LibraryProduct("libamcl_pairing_BLS461", :libamcl_pairing_BLS461),
    LibraryProduct("libamcl_pairing_BLS48", :libamcl_pairing_BLS48),
    LibraryProduct("libamcl_pairing_BN254", :libamcl_pairing_BN254),
    LibraryProduct("libamcl_pairing_BN254CX", :libamcl_pairing_BN254CX),
    LibraryProduct("libamcl_pairing_FP256BN", :libamcl_pairing_FP256BN),
    LibraryProduct("libamcl_pairing_FP512BN", :libamcl_pairing_FP512BN),
    LibraryProduct("libamcl_rsa_2048", :libamcl_rsa_2048),
    LibraryProduct("libamcl_rsa_3072", :libamcl_rsa_3072),
    LibraryProduct("libamcl_rsa_4096", :libamcl_rsa_4096),
    LibraryProduct("libamcl_wcc_BLS24", :libamcl_wcc_BLS24),
    LibraryProduct("libamcl_wcc_BLS381", :libamcl_wcc_BLS381),
    LibraryProduct("libamcl_wcc_BLS383", :libamcl_wcc_BLS383),
    LibraryProduct("libamcl_wcc_BLS461", :libamcl_wcc_BLS461),
    LibraryProduct("libamcl_wcc_BLS48", :libamcl_wcc_BLS48),
    LibraryProduct("libamcl_wcc_BN254", :libamcl_wcc_BN254),
    LibraryProduct("libamcl_wcc_BN254CX", :libamcl_wcc_BN254CX),
    LibraryProduct("libamcl_wcc_FP256BN", :libamcl_wcc_FP256BN),
    LibraryProduct("libamcl_wcc_FP512BN", :libamcl_wcc_FP512BN),
    LibraryProduct("libamcl_x509", :libamcl_x509),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
#= NOTE:

- general rule: use latex command for either block or inline, avoid mixing (1, 2, 3)
- if a command defines a block, use empty lines to separate otherwise it will be plugged in a paragraph (3)
- if mixing inline and block, you may get unexpected stuff
    - ok (4.a)
    - nok (4.b)
- nesting (5)
=#

@testset "1/noblock nop" begin
    h = raw"""
        \newcommand{\foo}{A **B** C}
        1 \foo 2
        """ |> fd2html
    @test h // "
1 A B C 2
"
end

@testset "2/block nop" begin
    h = raw"""
        \newcommand{\foo}{```julia
        x = 5
        @show x
        ```}
        1 \foo 2
        """ |> fd2html
    @test h // """
1
x = 5
                  @show x
 2
"""
end

@testset "3/block wp" begin
    h = raw"""
        \newcommand{\foo}{```julia
        x = 5
        @show x
        ```}
        1

        \foo

        2""" |> fd2html
    @test h // """
1

x = 5
                  @show x

2
"""
end

@testset "4/mixing" begin
    # OK (4.a)
    h = raw"""
        \newcommand{\foo}{A $$ x = 5 $$ B}
        1 \foo 2
        """ |> fd2html
    @test h // raw"""

1 A

               \[ x = 5 \]

B 2

               """
    # XXX NOK (4.b)
    h = raw"""
       \newcommand{\foo}{A $$ x = 5 $$ B}
       1

       \foo

       2
       """ |> fd2html
    @test h // raw"""

1

               A

               \[ x = 5 \]

B

2
"""
end

@testset "5/nesting" begin
    h = raw"""
        \newcommand{\foo}{A `B` C}
        \newcommand{\bar}{1 \foo 2}
        aa \bar bb
        """ |> fd2html
    @test h // "
aa 1 A B C 2 bb
"
    h = raw"""
        \newcommand{\foo}[1]{@@b #1 @@}
        \newcommand{\bar}[1]{A `!#1` \foo{g} C}
        aa \bar{hh} bb
        """ |> fd2html
    @test h // """

aa A hh
g
 C bb

               """
    # nesting with block // mixing
    h = raw"""
        \newcommand{\foo}[1]{```python
        !#1
        ```}
        \newcommand{\bar}[1]{ABC \foo{!#1} DEF}
        aa \bar{x=1} bb
        """ |> fd2html
    @test h // raw"""

aa ABC
x=1
 DEF bb

               """
end

@testset "6/div-mix" begin
    s = raw"""
        \newcommand{\note}[1]{@@note #1 @@}

        \note{A}

        \note{A `B` C}

        \note{A @@cc B @@ D}

        \note{A @@cc B `D` E @@ F}
        """ |> fd2html_td
    @test isapproxstr(s, """

A

A B C

A

B

D

A

B D E

F

        """)
end


text:
using SpecialFunctions: erfc, erf

"""
@article{riveros1993review,
  title={Review of ϕ(ρz) curves in electron probe microanalysis},
  author={Riveros, Jose and Castellano, Gustavo},
  journal={X-Ray Spectrometry},
  volume={22},
  number={1},
  pages={3--10},
  year={1993},
  publisher={Wiley Online Library}
}

The instruction in Packwood1991 in Electron Probe Quantitation is:
\"For compounds weight averaging is used for all appropriate variables: Z, Z/A, η, and (Z/A)log(1.166(E0-Ec)/2J)\"
"""
struct Riveros1993 <: MatrixCorrection
    material::Material
    E0::Float64  # Beam energy
    subshell::Union{Nothing,AtomicSubShell}
    Ea::Float64 # Sub-shell energy
    γ::Float64
    ϕ0::Float64
    α::Float64
    β::Float64

    function Riveros1993(mat::Material, ashell::Union{AtomicSubShell,Nothing}, ea::Float64, e0::Float64)
        e0k, eck, u0 = 0.001 * e0, 0.001 * ea, e0 / ea
        @assert u0 >= 1.0
        αz(elm) = (2.14e5*z(elm)^1.16/(a(elm)*e0k^1.25))*sqrt(log(1.166*e0/J(Berger1982,elm))/(e0k-eck))
        βz(elm) = (1.1e5*z(elm)^1.5)/((e0k-eck)*a(elm))
        weightavg(f, mat) = sum(f(elm)*mat[elm] for elm in keys(mat))
        ηm = η(mat, e0) # Use Donovan's averaging (not weight averaging...)
        ϕ0 = 1.0 + (ηm*u0*log(u0))/(u0-1.0)
        γ = (1.0 + ηm)*(u0*log(u0))/(u0-1.0)
        α = weightavg(αz, mat)
        β = weightavg(βz, mat)
        @assert α > 0.0
        @assert β > 0.0
        return new(mat, e0, ashell, ea, γ, ϕ0, α, β)
    end
end

ϕ(rv::Riveros1993, ρz) = exp(-(rv.α*ρz)^2)*(rv.γ - (rv.γ-rv.ϕ0)*exp(-rv.β*ρz))

function F(rv::Riveros1993)
    gg = rv.β/(2.0*rv.α)
    return gg < 22.3 ? (sqrt(π)*(rv.γ - exp(gg^2)*(rv.γ - rv.ϕ0)*erfc(gg)))/(2.0*rv.α) : 1.0
end

function Fχ(rv::Riveros1993, θtoa::Float64)
    @assert isnothing(rv.subshell)  "Use only for continuum correction"
    χm = χ(material(rv), rv.Ea, θtoa)
    ff, gg = χm/(2.0*rv.α), (rv.β + χm)/(2.0*rv.α)
    return gg < 22.3 ? (sqrt(π)*(exp(ff^2)*rv.γ*rv.α*(1.0 - erf(ff)) -
            exp(gg^2)*(rv.γ - rv.ϕ0)*rv.α*(1.0 - erf(gg))))/(2.0*rv.α^2) :
            0.0
end

function Fχ(rv::Riveros1993, xray::CharXRay, θtoa::Float64)
    @assert !isnothing(rv.subshell) "Use only for characteristic correction"
    @assert inner(xray) == rv.subshell
    χm = χ(material(rv), xray, θtoa)
    ff, gg = χm/(2.0*rv.α), (rv.β + χm)/(2.0*rv.α)
    return gg < 22.3 ? (sqrt(π)*(exp(ff^2)*rv.γ*rv.α*(1.0 - erf(ff)) -
           exp(gg^2)*(rv.γ - rv.ϕ0)*rv.α*(1.0 - erf(gg))))/(2.0*rv.α^2) :
           0.0
end

function Fχp(rv::Riveros1993, xray::CharXRay, θtoa::Float64, τ::Float64)
    @assert isnothing(rv.subshell) || (inner(xray) == rv.subshell)
    χm = χ(material(rv), xray, θtoa)
    @assert rv.α > 0.0 && χm > 0.0
    ff, gg = χm/(2.0*rv.α), (rv.β + χm)/(2.0*rv.α)
    return gg < 22.3 ? (sqrt(π)*rv.α*(exp(ff^2)*rv.γ*
            (-((χm*erf(ff^2))/χm) +
            ((2*rv.α^2*τ + χm)*erf((2.0*rv.α^2*τ + χm)/(2.0*rv.α)))/ (2*rv.α^2*τ + χm)) -
           exp(gg^2)*rv.γ*
            (-(((rv.β + χm)*erf(gg))/abs(rv.β + χm)) +
              ((rv.β + 2.0*rv.α^2*τ + χm)*erf(abs(rv.β + 2*rv.α^2*τ + χm))/(2.0*rv.α))/
               abs(rv.β + 2.0*rv.α^2*τ + χm))))/(2.0*rv.α^2) : 0.0
end

matrixcorrection(::Type{Riveros1993}, mat::Material, ashell::AtomicSubShell, e0::Float64) = Riveros1993(mat, ashell, energy(ashell), e0)

continuumcorrection(::Type{Riveros1993}, mat::Material, ea::Float64, e0::Float64) = Riveros1993(mat, nothing, ea, e0)

Base.range(::Type{Riveros1993}, mat::Material, e0::Float64, inclDensity=true) = range(XPP, mat, e0, inclDensity)

NeXLCore.minproperties(::Type{Riveros1993}) = (:BeamEnergy, :TakeOffAngle)

max_stars_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_issues_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_forks_repo_name: NicholasWMRitchie/NeXLMatrixCorrection


text:
@testset "BLAS" begin
    @testset "gemm" begin
        dims = 3:5
        for m in dims, n in dims, p in dims, tA in ('N', 'T'), tB in ('N', 'T')
            α = randn()
            A = randn(tA === 'N' ? (m, n) : (n, m))
            B = randn(tB === 'N' ? (n, p) : (p, n))
            C = gemm(tA, tB, α, A, B)
            ȳ = randn(size(C)...)
            rrule_test(
                gemm,
                ȳ,
                (tA, nothing),
                (tB, nothing),
                (α, randn()),
                (A, randn(size(A))),
                (B, randn(size(B))),
            )
        end
    end
    @testset "gemv" begin
        for n in 3:5, m in 3:5, t in ('N', 'T')
            α = randn()
            A = randn(m, n)
            x = randn(t === 'N' ? n : m)
            y = α * (t === 'N' ? A : A') * x
            ȳ = randn(size(y)...)
            rrule_test(
                gemv,
                ȳ,
                (t, nothing),
                (α, randn()),
                (A, randn(size(A))),
                (x, randn(size(x))),
            )
        end
    end
end


text:
# Run dynamic programming only within a band.
# Formally, a[i] and b[j] can be matched only if -L ≤ j - i ≤ U.
function affinegap_banded_global_align{T}(a, b, L::Int, U::Int, submat::AbstractSubstitutionMatrix{T}, gap_open_penalty::T, gap_extend_penalty::T)
    m = length(a)
    n = length(b)
    go = gap_open_penalty
    ge = gap_extend_penalty
    goe = go + ge
    L = min(L, m)
    U = min(U, n)
    # band width
    W = L + U + 1
    trace = Matrix{Trace}(W, n + 1)
    H = Vector{T}(W)
    E = Vector{T}(W)
    # In order to save the working space, the matrices are vertically sheared.
    # Namely, the coordinate is transformed as (i, j) → (i-j+U, j), where
    # (i, j) points to the subproblem of the optimal alignment between two
    # prefixes, a[1:i] and b[1:j], within the band.
    # The position (i, j) is accessed with:
    #   * trace[i-j+U+1,j+1]
    #   * H[i-j+U+1]
    #   * E[i-j+U+1]  (j ≥ 1)
    @inbounds begin
        # (i, j) = (0, 0)
        H[0-0+U+1] = T(0)
        trace[0-0+U+1,0+1] = TRACE_NONE
        for i in 1:L
            # j = 0
            H[i-0+U+1] = affinegap_score(i, go, ge)
            trace[i-0+U+1,0+1] = i == 1 ? TRACE_GAPOPEN_B : TRACE_GAPEXTD_B
        end
        # add gap_extend_penalty to avoid overflow for integers
        minimum = typemin(T) + ge
        for j in 1:n
            if j ≤ U
                # i = 0
                H[0-j+U+1] = affinegap_score(j, go, ge)
                trace[0-j+U+1,j+1] = j == 1 ? TRACE_GAPOPEN_A : TRACE_GAPEXTD_A
            end
            F = T(0)
            # vertical bounds along the j-th column
            lo = max(1, j - U)
            hi = min(m, j + L)
            for i in lo:hi
                # gap in the sequence A
                if j == 1 || i == j + L
                    e = minimum
                    gap_a = TRACE_NONE
                else
                    e = H[i-(j-1)+U+1] - goe
                    gap_a = TRACE_GAPOPEN_A
                    e′ = E[i-(j-1)+U+1] - ge
                    if e′ == e
                        gap_a |= TRACE_GAPEXTD_A
                    elseif e′ > e
                        gap_a  = TRACE_GAPEXTD_A
                        e = e′
                    end
                end
                # gap in the sequence B
                if i == 1 || j == i + U
                    f = minimum
                    gap_b = TRACE_NONE
                else
                    f = H[(i-1)-j+U+1] - goe
                    gap_b = TRACE_GAPOPEN_B
                    f′ = F - ge
                    if f′ == f
                        gap_b |= TRACE_GAPEXTD_B
                    elseif f′ > f
                        gap_b  = TRACE_GAPEXTD_B
                        f = f′
                    end
                end
                # match
                h = H[i-j+U+1] + submat[a[i],b[j]]
                # find the best score and its trace
                best = max(e, f, h)
                t = TRACE_NONE
                e == best && (t |= gap_a)
                f == best && (t |= gap_b)
                h == best && (t |= TRACE_MATCH)
                # update
                E[i-j+U+1] = e
                F = f
                H[i-j+U+1] = best
                trace[i-j+U+1,j+1] = t
            end
        end
    end
    return H[m-n+U+1], trace
end

function affinegap_banded_global_traceback(a, b, L, U, trace, endpos)
    L = min(L, length(a))
    U = min(U, length(b))
    # gap/character counts (reversed order)
    counts_a = [0, 0]
    counts_b = [0, 0]
    # if gap extension is selected in the previous traceback step, either
    # gap extension or gap open in that direction should be selected.
    gapext_a = false
    gapext_b = false
    i, j = endpos
    while i ≥ 1 && j ≥ 1
        @assert !(gapext_a && gapext_b)
        t = trace[i-j+U+1,j+1]
        if gapext_a
            if j ≥ 2 && t & TRACE_GAPEXTD_A > 0
                @gapext a
            elseif t & TRACE_GAPOPEN_A > 0
                @gapopen a
            end
        elseif gapext_b
            if i ≥ 2 && t & TRACE_GAPEXTD_B > 0
                @gapext b
            elseif t & TRACE_GAPOPEN_B > 0
                @gapopen b
            end
        elseif t & TRACE_MATCH > 0
            @match
        elseif j ≥ 2 && t & TRACE_GAPEXTD_A > 0
            @gapext a
        elseif t & TRACE_GAPOPEN_A > 0
            @gapopen a
        elseif i ≥ 2 && t & TRACE_GAPEXTD_B > 0
            @gapext b
        elseif t & TRACE_GAPOPEN_B > 0
            @gapopen b
        end
        # do not come here
        @assert false
    end
    while j ≥ 1 @gap a end
    while i ≥ 1 @gap b end
    reverse!(counts_a)
    reverse!(counts_b)
    return GappedSequence(a, counts_a), GappedSequence(b, counts_b)
end

function isinband(i, j, L, U, a, b)
    return 0 ≤ i ≤ length(a) && 0 ≤ j ≤ length(b) && i - L ≤ j ≤ i + U
end


text:
export process_raw, process_raw!, process_fft, rfft
export onebit!, onebit, remove_response!, remove_response, remove_amp!, remove_amp
export clip, clip!, clamp, clamp!, mute!, mute, phase
import Base: clamp, clamp!
import FFTW: rfft
"""
    process_raw!(S,fs)

Pre-process raw seismic data.

- Removes mean from each channel in `S`.
- Detrends each channel in `S`.
- Downsamples data to sampling rate `fs`
- Phase-shifts data to begin at 00:00:00.0

# Arguments
- `S::SeisData`: SeisData structure.
- `fs::Real`: Sampling rate to downsample `S`.
"""
function process_raw!(S::SeisData, fs::Real; ϕshift::Bool=true)
    merge!(S)
    ungap!(S)
	detrend!(S)         # remove mean & trend from channel
	taper!(S)                      # taper channel ends
	if fs ∉ S.fs
		filtfilt!(S,fh=Float64(fs/2),rt="Lowpass")    # lowpass filter before downsampling
	end
	resample!(S,fs=Float64(fs)) # downsample to lower fs
	taper!(S)
    phase_shift!(S, ϕshift=ϕshift) # timing offset from sampling period
    return nothing
end
process_raw(S::SeisData, fs::Real;
           ϕshift::Bool=true) = (U = deepcopy(S);
           process_raw!(U,fs, ϕshift=ϕshift); return U)

"""
   process_raw!(C,fs)

Pre-process raw SeisChannel.

- Removes mean from data in `C`.
- Detrends data in `C`.
- Downsamples data to sampling rate `fs`
- Phase-shifts data to begin at 00:00:00.0

# Arguments
- `C::SeisChannel`: SeisChannel structure.
- `fs::Real`: Sampling rate to downsample `C`.
"""
function process_raw!(C::SeisChannel, fs::Real; ϕshift::Bool=true)
	ungap!(C)
    detrend!(C)         # remove mean & trend from channel
    taper!(C)         # taper channel ends
    if fs != C.fs
    	filtfilt!(C,fh=Float64(fs/2),rt="Lowpass")    # lowpass filter before downsampling
	end
	resample!(C,Float64(fs)) # downsample to lower fs
	taper!(C)
	phase_shift!(C, ϕshift=ϕshift) # timing offset from sampling period
	return nothing
end
process_raw(C::SeisChannel, fs::Real;
        ϕshift::Bool=true) = (U = deepcopy(C);
        process_raw!(U,fs, ϕshift=ϕshift); return U)

"""

 rfft(R)

Computes windowed rfft of ambient noise data. Returns FFTData structure.

Overloads the rfft function from FFTW.

# Arguments
- `R::RawData`: RawData structure
"""
function rfft(R::RawData,dims::Int=1)
    FFT = rfft(R.x,dims)
    return FFTData(R.name, R.id,R.loc, R.fs, R.gain, R.freqmin, R.freqmax,
                 R.cc_len, R.cc_step, R.whitened, R.time_norm, R.resp,
                 R.misc, R.notes, R.t, FFT)
end

"""
    phase(A::AbstractArray)

Extract instantaneous phase from signal A.

For time series `A`, its analytic representation ``S = A + H(A)``, where
``H(A)`` is the Hilbert transform of `A`. The instantaneous phase ``e^{iθ}``
of `A` is given by dividing ``S`` by its modulus: ``e^{iθ} = \\frac{S}{|S|}``
For more information on Phase Cross-Correlation, see:
[Ventosa et al., 2019](https://pubs.geoscienceworld.org/ssa/srl/article-standard/570273/towards-the-processing-of-large-data-volumes-with).
"""
function phase(A::AbstractArray)
	Nrows = size(A,1)
    T = eltype(A)
    f = similar(A,fftouttype(T))
    fill!(f,fftouttype(T)(0))
    f[1:Nrows÷2 + 1 + isodd(Nrows),:] .= rfft(A,1)
    f[2:Nrows÷2  + isodd(Nrows),:] .*= T(2.)
	f[1:Nrows÷2 + 1 + isodd(Nrows),:] ./= abs.(f[1:Nrows÷2 + 1 + isodd(Nrows),:])
	ind = findall(isweird.(f[1:Nrows÷2 + isodd(Nrows) + 1,:]))
	if length(ind) > 0
		fill!(@view(f[ind]),fftouttype(T)(0))
	end
    return f
end

"""

 phase(R)

Computes windowed analytic signal of ambient noise data. Returns FFTData structure.

# Arguments
- `R::RawData`: RawData structure
"""
function phase(R::RawData)
    FFT = phase(R.x)
    return FFTData(R.name, R.id,R.loc, R.fs, R.gain, R.freqmin, R.freqmax,
                 R.cc_len, R.cc_step, R.whitened, R.time_norm, R.resp,
                 R.misc, R.notes, R.t, FFT)
end

isweird(x) =  isnan(x) .| isinf(x)

"""

  mute(A,factor)

Set high amplitudes in array `A` to zero.
Uses median of envelope of `A` to find outliers.
"""
function mute!(A::AbstractArray,factor::Real=3)
    T = eltype(A)
    envelope = abs.(hilbert(A))
    levels = mean(envelope,dims=1)
    level = factor .* median(levels)
    A[envelope .> level] .= T(0)
    return nothing
end
mute(A::AbstractArray,factor::Real=3) = (U = deepcopy(A); mute!(U,factor);
     return U)
mute!(R::RawData,factor::Real=3) = mute!(R.x,factor)
mute(R::RawData,factor::Real=3) = (U = deepcopy(R); mute!(U,factor);
     return U)

"""
  clip(A,factor)

Truncate array A at `factor` times the root mean square of each column.

#Arguments
- `A::AbstractArray`: N-d time series array
- `factor::Real`:
- `f::Function`: Input statistical function (e.g. rms, var, std, mad)
- `dims`: Dimension of `A` to apply clipping (defaults to 1)
"""
function clip!(A::AbstractArray{T,N}, factor::Real; f::Function=std,dims=1) where {T,N}
    if N == 1
        high = f(A) .* factor
        clamp!(@view(A[:]),-high,high)
    else
        high = f(@view(A[:,:]),dims=dims) .* factor
        for ii = 1:size(A,2)
            clamp!(@view(A[:,ii]),-high[ii],high[ii])
        end
    end
    return nothing
end
clip(A::AbstractArray, factor::Real; f::Function=std, dims=1) = (U = deepcopy(A);
     clip!(U,factor,f=f,dims=dims);return U)
clip!(R::RawData,factor::Real;f::Function=std,dims=1) = clip!(R.x,factor,f=f,dims=dims)
clip(R::RawData,factor::Real;f::Function=std,dims=1) = (U = deepcopy(R);
     clip!(U.x,factor,f=f,dims=dims); return U)
clamp!(R::RawData,val::Real) = clamp!(R.x,-abs(val),abs(val))
clamp!(R::RawData,lo::Real,hi::Real) = lo < hi ? clamp!(R.x,lo,hi) : throw(ArgumentError("lo value $lo must be less than hi value $hi."))
clamp(R::RawData,val::Real) = (U = deepcopy(R); clamp!(U,val);return U)
clamp(R::RawData,lo::Real,hi::Real) = (U = deepcopy(R); clamp!(U,lo,hi);return U)

"""
  remove_amp!(R)

Filter raw data based on amplitude.
"""
function remove_amp!(R::RawData; max_std::Real=10.)
    # remove nonzero columns
    zeroind = nonzero(R.x)
    if length(zeroind) == 0
        throw(ErrorException("All values in Rawdata == 0"))
    elseif size(R.x,2) != length(zeroind)
        R.x = R.x[:,zeroind]
        R.t = R.t[zeroind]
    end

    # amplitude threshold indices
    stdind = std_threshold(R.x,max_std)
    if length(stdind) == 0
        throw(ErrorException("All columns in Rawdata contain values > max_std"))
    elseif size(R.x,2) != length(stdind)
        R.x = R.x[:,stdind]
        R.t = R.t[stdind]
    end
    return nothing
end
remove_amp(R::RawData; max_std::Real=10.) = (U = deepcopy(R);remove_amp!(U,
          max_std=max_std); return U)

"""

  onebit!(R)

One-bit amplitude modification of RawData `R`.
"""
function onebit!(R::RawData)
    R.x .= sign.(R.x)
    return nothing
end
onebit(R::RawData) = (U = deepcopy(R); onebit!(U);return U)

"""
nonzero(A)

Find indices of all nonzero columns in array `A`.
"""
function nonzero(A::AbstractArray)
    Nrows, Ncols = size(A)
    ind = Int64[]
    sizehint!(ind,Ncols)
    for ii = 1:Ncols
        for jj = 1:Nrows
            if !iszero(A[jj,ii])
                append!(ind,ii)
                break
            end
        end
    end
    return ind
end


text:
##### Mass tendencies

#####
##### First order fluxes
#####

function flux(::Mass, ::Advect, lm::AtmosLinearModel, args)
    @unpack state = args
    return state.ρu
end

##### Momentum tendencies

#####
##### First order fluxes
#####

struct LinearPressureGradient <: TendencyDef{Flux{FirstOrder}} end

function flux(::Momentum, ::LinearPressureGradient, lm::AtmosLinearModel, args)
    @unpack state, aux = args
    s = state.ρu * state.ρu'
    pad = SArray{Tuple{size(s)...}}(ntuple(i -> 0, length(s)))
    pL = linearized_pressure(lm.atmos, state, aux)
    return pad + pL * I
end

#####
##### Sources (Momentum)
#####
function source(
    ::Momentum,
    s::Gravity,
    lm::AtmosAcousticGravityLinearModel,
    args,
)
    @unpack direction, state, aux = args
    if direction === VerticalDirection || direction === EveryDirection
        ∇Φ = ∇gravitational_potential(lm.atmos.orientation, aux)
        return -state.ρ * ∇Φ
    end
    FT = eltype(state)
    return SVector{3, FT}(0, 0, 0)
end

##### Energy tendencies

#####
##### First order fluxes
#####

struct LinearEnergyFlux <: TendencyDef{Flux{FirstOrder}} end

function flux(::Energy, ::LinearEnergyFlux, lm::AtmosAcousticLinearModel, args)
    @unpack state, aux = args
    ref = aux.ref_state
    e_pot = gravitational_potential(lm.atmos.orientation, aux)
    return ((ref.ρe + ref.p) / ref.ρ - e_pot) * state.ρu
end

function flux(
    ::Energy,
    ::LinearEnergyFlux,
    lm::AtmosAcousticGravityLinearModel,
    args,
)
    @unpack state, aux = args
    ref = aux.ref_state
    return ((ref.ρe + ref.p) / ref.ρ) * state.ρu
end


text:
"""
    scatterlines(xs, ys, [zs]; kwargs...)

Plots `scatter` markers and `lines` between them.

## Attributes
$(ATTRIBUTES)
"""
@recipe(ScatterLines) do scene
    s_theme = default_theme(scene, Scatter)
    l_theme = default_theme(scene, Lines)
    Attributes(
        color = l_theme.color,
        colormap = l_theme.colormap,
        colorrange = get(l_theme.attributes, :colorrange, automatic),
        linestyle = l_theme.linestyle,
        linewidth = l_theme.linewidth,
        markercolor = s_theme.color,
        markercolormap = s_theme.colormap,
        markercolorrange = get(s_theme.attributes, :colorrange, automatic),
        markersize = s_theme.markersize,
        strokecolor = s_theme.strokecolor,
        strokewidth = s_theme.strokewidth,
        marker = s_theme.marker,
    )
end

function plot!(p::Combined{scatterlines, <:NTuple{N, Any}}) where N
    lines!(p, p[1:N]...;
        color = p.color,
        linestyle = p.linestyle,
        linewidth = p.linewidth,
        colormap = p.colormap,
        colorrange = p.colorrange,
    )
    scatter!(p, p[1:N]...;
        color = p.markercolor,
        strokecolor = p.strokecolor,
        strokewidth = p.strokewidth,
        marker = p.marker,
        markersize = p.markersize,
        colormap = p.markercolormap,
        colorrange = p.markercolorrange,
    )
end


text:
using Sniffles
using Test
using JuMP
using GLPK

@testset "Sniffles" begin

end


text:
using Distributions, SpecialFunctions, Random
"""
GParticleFilter: pf for the Gaussian task (section 2.3.1)
"""
struct GParticleFilter <: GLearner
    nparticles::Int
    Neffthrs::Float64
    changeprobability::Float64
    m::Float64
    χ_0::Float64
    ν_0::Float64
    sigma::Float64 # std of observations
    particlesswitch::Array{Bool, 1} # wannabe nparticles
    weights::Array{Float64, 1} # wannabe nparticles
    Sgmparticles::Array{Float64, 1} # wannabe nparticles, Bayes Factor Surprise
    Sgm_n::Array{Float64, 1} # Bayes Factor Surprise
    χ_n::Array{Float64, 1}
    ν_n::Array{Float64, 1}
    mu_n::Array{Float64, 1} # parameter we are interested in
    var_n::Array{Float64, 1}
    γ_n::Array{Float64, 1}
    Ssh_n::Array{Float64, 1} # Shannon surprise
    seed::Any
    rng::MersenneTwister
    nresamplings::Array{Float64, 1} # track the number of times we resample
    Py_0_n::Array{Float64, 1}
end
function GParticleFilter(;nparticles = 6, changeprobability = .001,
                        mu_0 = 0., sigma_0 = 1., sigma = 0.1, seed = 3,
                        Neffthrs = nparticles/2.)
    particlesswitch = Array{Bool, 1}(undef, nparticles)
    weights = Array{Float64, 1}(undef, nparticles)
    Sgmparticles = Array{Float64, 1}(undef, nparticles)
    χ_n = Array{Float64, 1}(undef, nparticles)
    ν_n = Array{Float64, 1}(undef, nparticles)
    rng = MersenneTwister(seed)
    ν_0 = calcν(sigma, sigma_0)
    χ_0 = calcχ(mu_0, ν_0, sigma)
    particlesswitch .= false; weights .= 1. /nparticles
    Sgmparticles .= 0.; Sgm_n = [1.]
    χ_n .= copy(χ_0); ν_n .= copy(ν_0)
    m = changeprobability/ (1. - changeprobability)
    var_n_0 = sigma^2 / ν_0
    mu_n_0 = sigma * χ_0 / ν_0
    GParticleFilter(nparticles, Neffthrs, changeprobability, m, χ_0, ν_0,
                    sigma, particlesswitch, weights, Sgmparticles, Sgm_n,
                    χ_n, ν_n, [mu_n_0], [var_n_0], [m/(1. + m)], [1.], seed, rng,
                    [0.], [1.])
end
export GParticleFilter
function update!(learnerG::GParticleFilter, y)
    learnerG.particlesswitch .= false
    calcSgmparticles!(learnerG, y)
    calcSgm!(learnerG)
    calcγ!(learnerG)
    calcShannon!(learnerG, y)
    getweights!(learnerG, y)
    sampleparticles!(learnerG, y)
    # Evaluate Neff = approx nr of particles that have a weight which meaningfully contributes to the probability distribution.
    Neff = 1. /sum(learnerG.weights .^2)
    if Neff <= learnerG.Neffthrs; resample!(learnerG); end
    update_χ_ν!(learnerG, y)
    computemu_n!(learnerG)
    computevar_n!(learnerG)
end
export update!
function calcSgmparticles!(learnerG::GParticleFilter, y)
    for i in 1:learnerG.nparticles
        learnerG.Sgmparticles[i] = calcdoublefratio(learnerG.χ_0,
                                    learnerG.χ_n[i] + calcφ(learnerG, y),
                                    learnerG.χ_0 + calcφ(learnerG, y),
                                    learnerG.χ_n[i],
                                    learnerG.ν_0, learnerG.ν_n[i] + 1.,
                                    learnerG.ν_0 + 1., learnerG.ν_n[i])
        learnerG.Sgmparticles[i] = checkInf(learnerG.Sgmparticles[i])
    end
end
# Calculate Bayes Factor Surprise
function calcSgm!(learnerG::GParticleFilter)
    if any(isequal.(learnerG.Sgmparticles, 0.))
        learnerG.Sgm_n[1] = 0.
    else
        SgmInvWeighted = learnerG.weights ./ learnerG.Sgmparticles
        learnerG.Sgm_n[1] = 1. / sum(SgmInvWeighted)
    end
end
function getweights!(learnerG::GParticleFilter, y)
    weightscopy = deepcopy(learnerG.weights)
    for i in 1:learnerG.nparticles
        wBterm = calcwBterm(learnerG, weightscopy, y, i)
        learnerG.weights[i] *= (1. - learnerG.γ_n[1]) * wBterm + learnerG.γ_n[1]
    end
    weigthsnorm = learnerG.weights ./ sum(learnerG.weights)
    if ~all(isapprox(learnerG.weights, weigthsnorm))#, atol = 0.0001))
        println("---- Here! ------------ ")
        @show learnerG.sigma
        @show learnerG.changeprobability
        @show learnerG.seed
        @show learnerG.weights
        @show weigthsnorm
    end
end
export getweights!
function calcwBterm(learnerG::GParticleFilter, weightscopy, y, i)
    wBdenom = 0.
    for j in 1:learnerG.nparticles
        if j == i
            wBdenom += weightscopy[i]
        else
            wBdoublefratio = calcwBdoublefratio(learnerG, y, j, i)
            wBdoublefratio = checkInf(wBdoublefratio)
            wBdenom += weightscopy[j] * wBdoublefratio
        end
    end
    wBterm = 1. / wBdenom
end
function calcwBdoublefratio(learnerG::Union{GParticleFilter, GBayesFilter, GMPN, GSOR},
                            y, j, i)
    calcdoublefratio(learnerG.χ_n[j], learnerG.χ_n[i] + calcφ(learnerG, y),
                    learnerG.χ_n[i], learnerG.χ_n[j] + calcφ(learnerG, y),
                    learnerG.ν_n[j], learnerG.ν_n[i] + 1.,
                    learnerG.ν_n[i], learnerG.ν_n[j] + 1.)
end
function sampleparticles!(learnerG::GParticleFilter, y)
    for i in 1:learnerG.nparticles
        particlestayprobability = computeproposaldistribution(learnerG, i)
        # Draw and possibly update
        r = rand(learnerG.rng)
        if r > particlestayprobability
            learnerG.particlesswitch[i] = true
        end
    end
end
export sampleparticles!
function computeproposaldistribution(learnerG::GParticleFilter, i)
    particlestayprobability = 1. /(1. + (learnerG.m * learnerG.Sgmparticles[i]))
end
export computeproposaldistribution
function update_χ_ν!(learnerG::GParticleFilter, y)
    for i in 1:learnerG.nparticles
        if learnerG.particlesswitch[i] # if new hidden state
            learnerG.χ_n[i] = copy(learnerG.χ_0)
            learnerG.ν_n[i] = copy(learnerG.ν_0)
        end
        learnerG.χ_n[i] += calcφ(learnerG, y)
        learnerG.ν_n[i] += 1
    end
end
export update_χ_ν!
function resample!(learnerG::GParticleFilter)
    learnerG.nresamplings[1] += 1.
    tempcopyparticleweights = deepcopy(learnerG.weights)
    d = Categorical(tempcopyparticleweights)
    tempcopyparticlesswitch = deepcopy(learnerG.particlesswitch)
    tempcopyχ_n = deepcopy(learnerG.χ_n)
    tempcopyν_n = deepcopy(learnerG.ν_n)
    for i in 1:learnerG.nparticles
        sampledindex = rand(learnerG.rng, d)
        learnerG.particlesswitch[i] = deepcopy(tempcopyparticlesswitch[sampledindex])
        learnerG.weights[i] = 1. /learnerG.nparticles
        learnerG.χ_n[i] = deepcopy(tempcopyχ_n[sampledindex])
        learnerG.ν_n[i] = deepcopy(tempcopyν_n[sampledindex])
    end
end
export resample!
function computemu_n!(learnerG::Union{GParticleFilter, GBayesFilter, GMPN, GSOR})
    mweighted = learnerG.weights .* (learnerG.sigma * learnerG.χ_n ./ learnerG.ν_n)
    learnerG.mu_n[1] = sum(mweighted)
end
function computevar_n!(learnerG::Union{GParticleFilter, GBayesFilter, GMPN, GSOR})
    μ_i = learnerG.sigma * learnerG.χ_n ./ learnerG.ν_n
    var_i = learnerG.sigma ./ learnerG.ν_n
    learnerG.var_n[1] = sum(learnerG.weights .* (μ_i .^2 + var_i)) - learnerG.mu_n[1]^2
end


text:
# For now dummy runtests.jl
# (As StanSample is available in test environment)


text:
##
## Join / merge
##

# Like similar, but returns a array that can have missings and is initialized with missings
similar_missing(dv::AbstractArray{T}, dims::Union{Int, Tuple{Vararg{Int}}}) where {T} =
    fill!(similar(dv, Union{T, Missing}, dims), missing)

const OnType = Union{Symbol, NTuple{2,Symbol}, Pair{Symbol,Symbol}}

# helper structure for DataFrames joining
struct DataFrameJoiner{DF1<:AbstractDataFrame, DF2<:AbstractDataFrame}
    dfl::DF1
    dfr::DF2
    dfl_on::DF1
    dfr_on::DF2
    left_on::Vector{Symbol}
    right_on::Vector{Symbol}

    function DataFrameJoiner{DF1, DF2}(dfl::DF1, dfr::DF2,
                                       on::Union{<:OnType, AbstractVector{<:OnType}}) where {DF1, DF2}

        on_cols = isa(on, Vector) ? on : [on]
        if eltype(on_cols) == Symbol
            left_on = on_cols
            right_on = on_cols
        else
            left_on = [first(x) for x in on_cols]
            right_on = [last(x) for x in on_cols]
        end
        new(dfl, dfr, dfl[left_on], dfr[right_on], left_on, right_on)
    end
end

DataFrameJoiner(dfl::DF1, dfr::DF2, on::Union{<:OnType, AbstractVector{<:OnType}}) where
    {DF1<:AbstractDataFrame, DF2<:AbstractDataFrame} =
    DataFrameJoiner{DF1,DF2}(dfl, dfr, on)

# helper map between the row indices in original and joined table
struct RowIndexMap
    "row indices in the original table"
    orig::Vector{Int}
    "row indices in the resulting joined table"
    join::Vector{Int}
end

Base.length(x::RowIndexMap) = length(x.orig)

# composes the joined data table using the maps between the left and right
# table rows and the indices of rows in the result

function compose_joined_table(joiner::DataFrameJoiner, kind::Symbol,
                              left_ixs::RowIndexMap, leftonly_ixs::RowIndexMap,
                              right_ixs::RowIndexMap, rightonly_ixs::RowIndexMap;
                              makeunique::Bool=false)
    @assert length(left_ixs) == length(right_ixs)
    # compose left half of the result taking all left columns
    all_orig_left_ixs = vcat(left_ixs.orig, leftonly_ixs.orig)

    ril = length(right_ixs)
    lil = length(left_ixs)
    loil = length(leftonly_ixs)
    roil = length(rightonly_ixs)

    if loil > 0
        # combine the matched (left_ixs.orig) and non-matched (leftonly_ixs.orig) indices of the left table rows
        # preserving the original rows order
        all_orig_left_ixs = similar(left_ixs.orig, lil + loil)
        @inbounds all_orig_left_ixs[left_ixs.join] = left_ixs.orig
        @inbounds all_orig_left_ixs[leftonly_ixs.join] = leftonly_ixs.orig
    else
        # the result contains only the left rows that are matched to right rows (left_ixs)
        all_orig_left_ixs = left_ixs.orig # no need to copy left_ixs.orig as it's not used elsewhere
    end
    # permutation to swap rightonly and leftonly rows
    right_perm = vcat(1:ril, ril+roil+1:ril+roil+loil, ril+1:ril+roil)
    if length(leftonly_ixs) > 0
        # compose right_perm with the permutation that restores left rows order
        right_perm[vcat(right_ixs.join, leftonly_ixs.join)] = right_perm[1:ril+loil]
    end
    all_orig_right_ixs = vcat(right_ixs.orig, rightonly_ixs.orig)

    # compose right half of the result taking all right columns excluding on
    dfr_noon = without(joiner.dfr, joiner.right_on)

    nrow = length(all_orig_left_ixs) + roil
    @assert nrow == length(all_orig_right_ixs) + loil
    ncleft = ncol(joiner.dfl)
    cols = Vector{AbstractVector}(undef, ncleft + ncol(dfr_noon))
    # inner and left joins preserve non-missingness of the left frame
    _similar_left = kind == :inner || kind == :left ? similar : similar_missing
    for (i, col) in enumerate(eachcol(joiner.dfl))
        cols[i] = _similar_left(col, nrow)
        copyto!(cols[i], view(col, all_orig_left_ixs))
    end
    # inner and right joins preserve non-missingness of the right frame
    _similar_right = kind == :inner || kind == :right ? similar : similar_missing
    for (i, col) in enumerate(eachcol(dfr_noon))
        cols[i+ncleft] = _similar_right(col, nrow)
        copyto!(cols[i+ncleft], view(col, all_orig_right_ixs))
        permute!(cols[i+ncleft], right_perm)
    end
    res = DataFrame(cols, vcat(names(joiner.dfl), names(dfr_noon)),
                    makeunique=makeunique, copycols=false)

    if length(rightonly_ixs.join) > 0
        # some left rows are missing, so the values of the "on" columns
        # need to be taken from the right
        for (on_col_ix, on_col) in enumerate(joiner.left_on)
            # fix the result of the rightjoin by taking the nonmissing values from the right table
            offset = nrow - length(rightonly_ixs.orig) + 1
            copyto!(res[on_col], offset, view(joiner.dfr_on[on_col_ix], rightonly_ixs.orig))
        end
    end
    if kind ∈ (:right, :outer) && !isempty(rightonly_ixs.join)
        # At this point on-columns of the result allow missing values, because
        # right-only rows were filled with missing values when processing joiner.dfl
        # However, when the right on-column (plus the left one for the outer join)
        # does not allow missing values, the result should also disallow them.
        for (on_col_ix, on_col) in enumerate(joiner.left_on)
            LT = eltype(joiner.dfl_on[on_col_ix])
            RT = eltype(joiner.dfr_on[on_col_ix])
            if !(RT >: Missing) && (kind == :right || !(LT >: Missing))
                res[on_col] = disallowmissing(res[on_col])
            end
        end
    end
    return res
end

# map the indices of the left and right joined tables
# to the indices of the rows in the resulting table
# if `nothing` is given, the corresponding map is not built
function update_row_maps!(left_table::AbstractDataFrame,
                          right_table::AbstractDataFrame,
                          right_dict::RowGroupDict,
                          left_ixs::Union{Nothing, RowIndexMap},
                          leftonly_ixs::Union{Nothing, RowIndexMap},
                          right_ixs::Union{Nothing, RowIndexMap},
                          rightonly_mask::Union{Nothing, Vector{Bool}})
    # helper functions
    @inline update!(ixs::Nothing, orig_ix::Int, join_ix::Int, count::Int = 1) = nothing
    @inline function update!(ixs::RowIndexMap, orig_ix::Int, join_ix::Int, count::Int = 1)
        n = length(ixs.orig)
        resize!(ixs.orig, n+count)
        ixs.orig[n+1:end] .= orig_ix
        append!(ixs.join, join_ix:(join_ix+count-1))
        ixs
    end
    @inline update!(ixs::Nothing, orig_ixs::AbstractArray, join_ix::Int) = nothing
    @inline function update!(ixs::RowIndexMap, orig_ixs::AbstractArray, join_ix::Int)
        append!(ixs.orig, orig_ixs)
        append!(ixs.join, join_ix:(join_ix+length(orig_ixs)-1))
        ixs
    end
    @inline update!(ixs::Nothing, orig_ixs::AbstractArray) = nothing
    @inline update!(mask::Vector{Bool}, orig_ixs::AbstractArray) = (mask[orig_ixs] .= false)

    # iterate over left rows and compose the left<->right index map
    right_dict_cols = ntuple(i -> right_dict.df[i], ncol(right_dict.df))
    left_table_cols = ntuple(i -> left_table[i], ncol(left_table))
    next_join_ix = 1
    for l_ix in 1:nrow(left_table)
        r_ixs = findrows(right_dict, left_table, right_dict_cols, left_table_cols, l_ix)
        if isempty(r_ixs)
            update!(leftonly_ixs, l_ix, next_join_ix)
            next_join_ix += 1
        else
            update!(left_ixs, l_ix, next_join_ix, length(r_ixs))
            update!(right_ixs, r_ixs, next_join_ix)
            update!(rightonly_mask, r_ixs)
            next_join_ix += length(r_ixs)
        end
    end
end

# map the row indices of the left and right joined tables
# to the indices of rows in the resulting table
# returns the 4-tuple of row indices maps for
# - matching left rows
# - non-matching left rows
# - matching right rows
# - non-matching right rows
# if false is provided, the corresponding map is not built and the
# tuple element is empty RowIndexMap
function update_row_maps!(left_table::AbstractDataFrame,
                          right_table::AbstractDataFrame,
                          right_dict::RowGroupDict,
                          map_left::Bool, map_leftonly::Bool,
                          map_right::Bool, map_rightonly::Bool)
    init_map(df::AbstractDataFrame, init::Bool) = init ?
        RowIndexMap(sizehint!(Vector{Int}(), nrow(df)),
                    sizehint!(Vector{Int}(), nrow(df))) : nothing
    to_bimap(x::RowIndexMap) = x
    to_bimap(::Nothing) = RowIndexMap(Vector{Int}(), Vector{Int}())

    # init maps as requested
    left_ixs = init_map(left_table, map_left)
    leftonly_ixs = init_map(left_table, map_leftonly)
    right_ixs = init_map(right_table, map_right)
    rightonly_mask = map_rightonly ? fill(true, nrow(right_table)) : nothing
    update_row_maps!(left_table, right_table, right_dict, left_ixs, leftonly_ixs, right_ixs, rightonly_mask)
    if map_rightonly
        rightonly_orig_ixs = findall(rightonly_mask)
        rightonly_ixs = RowIndexMap(rightonly_orig_ixs,
                                    collect(length(right_ixs.orig) +
                                            (leftonly_ixs === nothing ? 0 : length(leftonly_ixs)) .+
                                            (1:length(rightonly_orig_ixs))))
    else
        rightonly_ixs = nothing
    end

    return to_bimap(left_ixs), to_bimap(leftonly_ixs), to_bimap(right_ixs), to_bimap(rightonly_ixs)
end

"""
    join(df1, df2; on = Symbol[], kind = :inner, makeunique = false,
         indicator = nothing, validate = (false, false))

Join two `DataFrame` objects

### Arguments

* `df1`, `df2` : the two AbstractDataFrames to be joined

### Keyword Arguments

* `on` : A column, or vector of columns to join df1 and df2 on. If the column(s)
    that df1 and df2 will be joined on have different names, then the columns
    should be `(left, right)` tuples or `left => right` pairs, or a vector of
    such tuples or pairs. `on` is a required argument for all joins except for
    `kind = :cross`

* `kind` : the type of join, options include:

  - `:inner` : only include rows with keys that match in both `df1`
    and `df2`, the default
  - `:outer` : include all rows from `df1` and `df2`
  - `:left` : include all rows from `df1`
  - `:right` : include all rows from `df2`
  - `:semi` : return rows of `df1` that match with the keys in `df2`
  - `:anti` : return rows of `df1` that do not match with the keys in `df2`
  - `:cross` : a full Cartesian product of the key combinations; every
    row of `df1` is matched with every row of `df2`

* `makeunique` : if `false` (the default), an error will be raised
  if duplicate names are found in columns not joined on;
  if `true`, duplicate names will be suffixed with `_i`
  (`i` starting at 1 for the first duplicate).

* `indicator` : Default: `nothing`. If a `Symbol`, adds categorical indicator
   column named `Symbol` for whether a row appeared in only `df1` (`"left_only"`),
   only `df2` (`"right_only"`) or in both (`"both"`). If `Symbol` is already in use,
   the column name will be modified if `makeunique=true`.

* `validate` : whether to check that columns passed as the `on` argument
   define unique keys in each input data frame (according to [`isequal`](@ref)).
   Can be a tuple or a pair, with the first element indicating whether to
   run check for `df1` and the second element for `df2`.
   By default no check is performed.

For the three join operations that may introduce missing values (`:outer`, `:left`,
and `:right`), all columns of the returned data table will support missing values.

When merging `on` categorical columns that differ in the ordering of their levels, the
ordering of the left `DataFrame` takes precedence over the ordering of the right `DataFrame`

### Result

* `::DataFrame` : the joined DataFrame

### Examples

```julia
name = DataFrame(ID = [1, 2, 3], Name = ["John Doe", "Jane Doe", "Joe Blogs"])
job = DataFrame(ID = [1, 2, 4], Job = ["Lawyer", "Doctor", "Farmer"])

join(name, job, on = :ID)
join(name, job, on = :ID, kind = :outer)
join(name, job, on = :ID, kind = :left)
join(name, job, on = :ID, kind = :right)
join(name, job, on = :ID, kind = :semi)
join(name, job, on = :ID, kind = :anti)
join(name, job, kind = :cross)

job2 = DataFrame(identifier = [1, 2, 4], Job = ["Lawyer", "Doctor", "Farmer"])
join(name, job2, on = (:ID, :identifier))
join(name, job2, on = :ID => :identifier)
```

"""
function Base.join(df1::AbstractDataFrame,
                   df2::AbstractDataFrame;
                   on::Union{<:OnType, AbstractVector{<:OnType}} = Symbol[],
                   kind::Symbol = :inner, makeunique::Bool=false,
                   indicator::Union{Nothing, Symbol} = nothing,
                   validate::Union{Pair{Bool, Bool}, Tuple{Bool, Bool}}=(false, false))
    if indicator !== nothing
        indicator_cols = ["_left", "_right"]
        for i in 1:2
            while (haskey(index(df1), Symbol(indicator_cols[i])) ||
                   haskey(index(df2), Symbol(indicator_cols[i])) ||
                   Symbol(indicator_cols[i]) == indicator)
                 indicator_cols[i] *= 'X'
            end
        end
        df1 = copy(df1, copycols=false)
        df1[Symbol(indicator_cols[1])] = trues(nrow(df1))
        df2 = copy(df2, copycols=false)
        df2[Symbol(indicator_cols[2])] = trues(nrow(df2))
    end

    if kind == :cross
        (on == Symbol[]) || throw(ArgumentError("Cross joins don't use argument 'on'."))
        return crossjoin(df1, df2, makeunique=makeunique)
    elseif on == Symbol[]
        throw(ArgumentError("Missing join argument 'on'."))
    end

    joiner = DataFrameJoiner(df1, df2, on)

    # Check merge key validity
    left_invalid = validate[1] ? any(nonunique(joiner.dfl, joiner.left_on)) : false
    right_invalid = validate[2] ? any(nonunique(joiner.dfr, joiner.right_on)) : false

    if left_invalid && right_invalid
        first_error_df1 = findfirst(nonunique(joiner.dfl, joiner.left_on))
        first_error_df2 = findfirst(nonunique(joiner.dfr, joiner.right_on))
        throw(ArgumentError("Merge key(s) are not unique in both df1 and df2. " *
                            "First duplicate in df1 at $first_error_df1. " *
                            "First duplicate in df2 at $first_error_df2"))
    elseif left_invalid
        first_error = findfirst(nonunique(joiner.dfl, joiner.left_on))
        throw(ArgumentError("Merge key(s) in df1 are not unique. " *
                            "First duplicate at row $first_error"))
    elseif right_invalid
        first_error = findfirst(nonunique(joiner.dfr, joiner.right_on))
        throw(ArgumentError("Merge key(s) in df2 are not unique. " *
                            "First duplicate at row $first_error"))
    end

    if kind == :inner
        joined = compose_joined_table(joiner, kind, update_row_maps!(joiner.dfl_on, joiner.dfr_on,
                                                                     group_rows(joiner.dfr_on),
                                                                     true, false, true, false)...,
                                      makeunique=makeunique)
    elseif kind == :left
        joined = compose_joined_table(joiner, kind, update_row_maps!(joiner.dfl_on, joiner.dfr_on,
                                                            group_rows(joiner.dfr_on),
                                                            true, true, true, false)...,
                                      makeunique=makeunique)
    elseif kind == :right
        joined = compose_joined_table(joiner, kind, update_row_maps!(joiner.dfr_on, joiner.dfl_on,
                                                            group_rows(joiner.dfl_on),
                                                            true, true, true, false)[[3, 4, 1, 2]]...,
                                      makeunique=makeunique)
    elseif kind == :outer
        joined = compose_joined_table(joiner, kind, update_row_maps!(joiner.dfl_on, joiner.dfr_on,
                                                            group_rows(joiner.dfr_on),
                                                            true, true, true, true)...,
                                      makeunique=makeunique)
    elseif kind == :semi
        # hash the right rows
        dfr_on_grp = group_rows(joiner.dfr_on)
        # iterate over left rows and leave those found in right
        left_ixs = Vector{Int}()
        sizehint!(left_ixs, nrow(joiner.dfl))
        dfr_on_grp_cols = ntuple(i -> dfr_on_grp.df[i], ncol(dfr_on_grp.df))
        dfl_on_cols = ntuple(i -> joiner.dfl_on[i], ncol(joiner.dfl_on))
        @inbounds for l_ix in 1:nrow(joiner.dfl_on)
            if findrow(dfr_on_grp, joiner.dfl_on, dfr_on_grp_cols, dfl_on_cols, l_ix) != 0
                push!(left_ixs, l_ix)
            end
        end
        joined = joiner.dfl[left_ixs, :]
    elseif kind == :anti
        # hash the right rows
        dfr_on_grp = group_rows(joiner.dfr_on)
        # iterate over left rows and leave those not found in right
        leftonly_ixs = Vector{Int}()
        sizehint!(leftonly_ixs, nrow(joiner.dfl))
        dfr_on_grp_cols = ntuple(i -> dfr_on_grp.df[i], ncol(dfr_on_grp.df))
        dfl_on_cols = ntuple(i -> joiner.dfl_on[i], ncol(joiner.dfl_on))
        @inbounds for l_ix in 1:nrow(joiner.dfl_on)
            if findrow(dfr_on_grp, joiner.dfl_on, dfr_on_grp_cols, dfl_on_cols, l_ix) == 0
                push!(leftonly_ixs, l_ix)
            end
        end
        joined = joiner.dfl[leftonly_ixs, :]
    else
        throw(ArgumentError("Unknown kind of join requested: $kind"))
    end

    if indicator !== nothing
        left = joined[Symbol(indicator_cols[1])]
        right = joined[Symbol(indicator_cols[2])]

        refs = UInt8[coalesce(l, false) + 2 * coalesce(r, false) for (l, r) in zip(left, right)]
        indicatorcol = CategoricalArray{String,1}(refs, CategoricalPool{String,UInt8}(["left_only", "right_only", "both"]))
        unique_indicator = indicator
        try_idx = 0
        while haskey(joined, unique_indicator)
            try_idx += 1
            unique_indicator = Symbol(string(indicator, "_", try_idx))
        end
        joined[unique_indicator] = indicatorcol

        deletecols!(joined, Symbol(indicator_cols[1]))
        deletecols!(joined, Symbol(indicator_cols[2]))
    end

    return joined
end

function crossjoin(df1::AbstractDataFrame, df2::AbstractDataFrame; makeunique::Bool=false)
    r1, r2 = size(df1, 1), size(df2, 1)
    colindex = merge(index(df1), index(df2), makeunique=makeunique)
    cols = Any[[repeat(c, inner=r2) for c in eachcol(df1)];
               [repeat(c, outer=r1) for c in eachcol(df2)]]
    DataFrame(cols, colindex, copycols=false)
end


text:
module MM1_Example

using EventSimulation
using Random
using Test

# Objectives of the example:
# * show basic use of monitor function

mutable struct MMQueue <: AbstractState
    ar::Float64 # arrival rate
    sr::Float64 # service rate
    len::Int    # queue length
    busy::Bool  # is server busy?
    load_data::Dict{Int, Float64} # store for monitor data
end

# here we accumulate information about queue lengths and periods
function monitor(s, Δ)
    ld = s.state.load_data
    load = s.state.len + s.state.busy
    if haskey(ld, load)
        ld[load] += Δ
    else
        ld[load] = Δ
    end
end

# customer arrival event
function arrival(s)
    if s.state.busy # server is working?
        s.state.len += 1
    else
        s.state.busy = true
        register!(s, leave, randexp()/s.state.sr)
    end
end

# customer leave the system event
function leave(s)
    if s.state.len > 0 # any customers waiting?
        s.state.len -= 1
        register!(s, leave, randexp()/s.state.sr)
    else
        s.state.busy = false
    end
end

function run(ar, sr)
    q = MMQueue(ar, sr, 0, false, Dict{Int, Float64}())
    s = Scheduler(q, Float64, monitor)
    repeat_register!(s, arrival, x -> randexp()/ar)
    go!(s, 1_000_000)
    ks = sort(collect(keys(q.load_data)))
    s = sum(values(q.load_data))
    ρ = ar / sr
    maxdiff = 0.0 # calculated for accuracy testing purposes
    for k in ks
        em = round(q.load_data[k]/s, digits=4)
        th = round((1-ρ)*ρ^k, digits=4)
        maxdiff = max(maxdiff, abs(em-th))
        println("$k\t=> ", rpad(em,6), "\t", th)
    end
    return(maxdiff)
end

@testset "Test of M/M/1 queue with monitor" begin
    @test run(1.5, 3.5) < 0.005
end

end # module