text:
    using Revise, Parameters, Logging, Test
    using CUDA
    includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
    using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils, ..PrepareArrtoBool, ..BitWiseUtils
    using ..ResultListUtils, ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils, ..ScanForDuplicates, ..Housdorff

    mainArrDims= (512,512,500)

    arrGold = zeros(Int32,mainArrDims);
    arrSegm = zeros(Int32,mainArrDims);

        # so we should get all mins as 2 and all maxes as 5

        arrGold[9,2,5]= 2
        arrGold[3,6,4]= 2
        arrGold[3,3,4]= 2

        arrGold[4,4,4]= 2
        arrSegm[4,4,4]= 2

        arrSegm[3,18,4]= 2
        arrSegm[5,2,11]= 2

        gold3d = CuArray(arrGold)
        segm3d= CuArray(arrSegm)
        numberToLooFor =Int32(2)
        robustnessPercent = 0.9

        boolKernelArgs, mainKernelArgs,threadsBoolKern,blocksBoolKern ,threadsMainKern,blocksMainKern ,shmemSizeBool,shmemSizeMain=    preparehousedorfKernel(gold3d,segm3d,robustnessPercent,numberToLooFor)
        mainArrDims,dataBdim,metaData,metaDataDims,reducedGoldA,reducedSegmA,loopXinPlane,loopYinPlane,minxRes,maxxRes,minyRes,maxyRes,minzRes,maxzRes,fn,fp ,numberToLooFor,inBlockLoopXZIterWithPadding,shmemblockDataLoop,shmemblockDataLenght,loopAXFixed,loopBXfixed,loopAYFixed,loopBYfixed,loopAZFixed,loopBZfixed,loopdataDimMainX,loopdataDimMainY,loopdataDimMainZ,inBlockLoopX,inBlockLoopY,inBlockLoopZ,metaDataLength,loopMeta,loopWarpMeta,clearIterResShmemLoop,clearIterSourceShmemLoop,resShmemTotalLength,sourceShmemTotalLength=boolKernelArgs

        arrGold[dataBdim[1]*3+1,dataBdim[2]*4+1,dataBdim[3]*5+1]= 2
        arrSegm[dataBdim[1]*3+1,dataBdim[2]*4+1,dataBdim[3]*5+1]= 2

        golddd = CuArray(arrGold);
        segmmm= CuArray(arrSegm);

        # using CUDA
        # ff = CUDA.zeros(1)
        # function gggg(ff)
        # #   shmemblockData = @cuDynamicSharedMem(Float32,(dataBdim[1], dataBdim[2]))
        # #   shmemblockData = @cuDynamicSharedMem(Float32,(1,2))
        #    shmemblockData = CuDynamicSharedArray(Int, 6)
        #     return
        # end
        # @cuda threads=(32,2) blocks=2 gggg(ff)
        # ff[1]

        metaDataDims
        metaDataLength
        metaDataLength/blocksBoolKern
        loopMeta
        dataBdim
        Int64(blocksBoolKern)
        function locForKernel(goldGPU,segmGPU,mainArrDims,dataBdim,metaData,metaDataDims,reducedGoldA,reducedSegmA,loopXinPlane,loopYinPlane,minxRes,maxxRes,minyRes,maxyRes,minzRes,maxzRes,fn,fp ,numberToLooFor,inBlockLoopXZIterWithPadding,shmemblockDataLoop,shmemblockDataLenght,loopAXFixed,loopBXfixed,loopAYFixed,loopBYfixed,loopAZFixed,loopBZfixed,loopdataDimMainX,loopdataDimMainY,loopdataDimMainZ,inBlockLoopX,inBlockLoopY,inBlockLoopZ,metaDataLength,loopMeta,loopWarpMeta,clearIterResShmemLoop,clearIterSourceShmemLoop,resShmemTotalLength,sourceShmemTotalLength)

          @getBoolCubeKernel()

        return
    end
    minxRes[1]

        @cuda shmem=shmemSizeBool threads=threadsBoolKern blocks=blocksBoolKern locForKernel(golddd,segmmm,mainArrDims,dataBdim,metaData,metaDataDims,reducedGoldA,reducedSegmA,loopXinPlane,loopYinPlane,minxRes,maxxRes,minyRes,maxyRes,minzRes,maxzRes,fn,fp ,numberToLooFor,inBlockLoopXZIterWithPadding,shmemblockDataLoop,shmemblockDataLenght,loopAXFixed,loopBXfixed,loopAYFixed,loopBYfixed,loopAZFixed,loopBZfixed,loopdataDimMainX,loopdataDimMainY,loopdataDimMainZ,inBlockLoopX,inBlockLoopY,inBlockLoopZ,metaDataLength,loopMeta,loopWarpMeta,clearIterResShmemLoop,clearIterSourceShmemLoop,resShmemTotalLength,sourceShmemTotalLength)
          Int64( maxzRes[1])

        @test Int64(fn[])==3
        @test Int64(fp[])==2

       Int64(minxRes[1])
         Int64( maxxRes[])
         Int64(minyRes[])
          Int64(maxyRes[])
          Int64( minzRes[])
          Int64( maxzRes[])

          @test  Int64(minxRes[])==1
          @test  Int64( maxxRes[])==4
          @test  Int64(minyRes[])==1
          @test   Int64(maxyRes[])==5
          @test  Int64( minzRes[])==1
          @test  Int64( maxzRes[])==6

        #   mainArrDims
        # @test Int64(minxRes[])==2
        # @test Int64( maxxRes[])==5
        # @test Int64(minyRes[])==1
        # @test  Int64(maxyRes[])==9
        # @test  Int64( minzRes[])==2
        # @test  Int64( maxzRes[])==6
        # arrGold[32*3+1,32*4+1,32*5+1]
        # arrSegm[32*3+1,32*4+1,32*5+1]

Int64(sum(reducedGoldA))

Int64(reducedGoldA[4,4,1])
Int64(reducedSegmA[4,4,1])

Int64(reducedGoldA[3,18,1])

        @test  isBit1AtPos(reducedGoldA[9,2,1],5)==true

        @test  isBit1AtPos(reducedGoldA[3,3,1],4)==true
        @test  isBit1AtPos(reducedGoldA[3,3,1],4)==true
        @test  isBit1AtPos(reducedGoldA[4,4,1],4)==true
        @test  isBit1AtPos(reducedGoldA[5,5,1],5)==false
        @test  isBit1AtPos(reducedGoldA[9,2,1],5)==true
        @test  isBit1AtPos(reducedSegmA[3,18,1],4)==true
        @test  isBit1AtPos(reducedGoldA[3,18,1],4)==false

        @test  isBit1AtPos(reducedSegmA[dataBdim[1]*3+1,dataBdim[2]*4+1,6],1)==true

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

 using Revise, Parameters, Logging, Test
 using CUDA
 includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
 using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils, ..PrepareArrtoBool, ..BitWiseUtils
 using ..ResultListUtils, ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils, ..ScanForDuplicates, ..Housdorff

 mainArrDims= (60,60,60);

 mainArrCPU= zeros(UInt8,mainArrDims);
 refArrCPU = zeros(UInt8,mainArrDims);
 ##### we will create two planes 20 units apart from each
mainArrCPU[10:50,10:50,10].= 1;
refArrCPU[10:50,10:50,30].= 1;

     gold3d = CuArray(mainArrCPU)
     segm3d= CuArray(refArrCPU)
     numberToLooFor =Int32(1)
     robustnessPercent = 0.9

     boolKernelArgs, mainKernelArgs,threadsBoolKern,blocksBoolKern ,threadsMainKern,blocksMainKern ,shmemSizeBool,shmemSizeMain=    preparehousedorfKernel(gold3d,segm3d,robustnessPercent,numberToLooFor)
     mainArrDims,dataBdim,metaData,metaDataDims,reducedGoldA,reducedSegmA,loopXinPlane,loopYinPlane,minxRes,maxxRes,minyRes,maxyRes,minzRes,maxzRes,fn,fp ,numberToLooFor,inBlockLoopXZIterWithPadding,shmemblockDataLoop,shmemblockDataLenght,loopAXFixed,loopBXfixed,loopAYFixed,loopBYfixed,loopAZFixed,loopBZfixed,loopdataDimMainX,loopdataDimMainY,loopdataDimMainZ,inBlockLoopX,inBlockLoopY,inBlockLoopZ,metaDataLength,loopMeta,loopWarpMeta,clearIterResShmemLoop,clearIterSourceShmemLoop,resShmemTotalLength,sourceShmemTotalLength=boolKernelArgs

     golddd = CuArray(mainArrCPU);
     segmmm= CuArray(refArrCPU);

     @cuda shmem=shmemSizeBool threads=threadsBoolKern blocks=blocksBoolKern boolKernelLoad(golddd,segmmm,boolKernelArgs...)
     Int64(sum(mainArrCPU))
     Int64(fn[])
     @test Int64(fn[])==Int64(sum(mainArrCPU))


text:
using eFEM, Test

# parameter types

# problem types

# solution types

ext: jl


text:
import UUIDs: UUID, uuid1
import .ExpressionExplorer: SymbolsState, FunctionNameSignaturePair, FunctionName
import .Configuration

"The (information needed to create the) dependency graph of a notebook. Cells are linked by the names of globals that they define and reference. 🕸"
Base.@kwdef struct NotebookTopology
    nodes::Dict{Cell,ReactiveNode} = Dict{Cell,ReactiveNode}()
end

# `topology[cell]` is a shorthand for `get(topology, cell, ReactiveNode())`
# with the performance benefit of only generating ReactiveNode() when needed
function Base.getindex(topology::NotebookTopology, cell::Cell)::ReactiveNode
    get!(ReactiveNode, topology.nodes, cell)
end

"Like a [`Diary`](@ref) but more serious. 📓"
mutable struct Notebook
    "Cells are ordered in a `Notebook`, and this order can be changed by the user. Cells will always have a constant UUID."
    cells::Array{Cell,1}

    # i still don't really know what an AbstractString is but it makes this package look more professional
    path::AbstractString
    notebook_id::UUID
    topology::NotebookTopology

    # buffer will contain all unfetched updates - must be big enough
    pendingupdates::Channel

    executetoken::Token

    # per notebook compiler options
    # nothing means to use global session compiler options
    compiler_options::Union{Nothing,Configuration.CompilerOptions}
end
# We can keep 128 updates pending. After this, any put! calls (i.e. calls that push an update to the notebook) will simply block, which is fine.
# This does mean that the Notebook can't be used if nothing is clearing the update channel.
Notebook(cells::Array{Cell,1}, path::AbstractString, notebook_id::UUID) =
    Notebook(cells, path, notebook_id, NotebookTopology(), Channel(1024), Token(), nothing)

Notebook(cells::Array{Cell,1}, path::AbstractString=numbered_until_new(joinpath(new_notebooks_directory(), cutename()))) = Notebook(cells, path, uuid1())

function cell_index_from_id(notebook::Notebook, cell_id::UUID)::Union{Int,Nothing}
    findfirst(c -> c.cell_id == cell_id, notebook.cells)
end

const _notebook_header = "### A Pluto.jl notebook ###"
# We use a creative delimiter to avoid accidental use in code
# so don't get inspired to suddenly use these in your code!
const _cell_id_delimiter = "# ╔═╡ "
const _order_delimiter = "# ╠═"
const _order_delimiter_folded = "# ╟─"
const _cell_suffix = "\n\n"

emptynotebook(args...) = Notebook([Cell()], args...)

"""
Save the notebook to `io`, `file` or to `notebook.path`.

In the produced file, cells are not saved in the notebook order. If `notebook.topolgy` is up-to-date, I will save cells in _topological order_. This guarantees that you can run the notebook file outside of Pluto, with `julia my_notebook.jl`.

Have a look at our [JuliaCon 2020 presentation](https://youtu.be/IAF8DjrQSSk?t=1085) to learn more!
"""
function save_notebook(io, notebook::Notebook)
    println(io, _notebook_header)
    println(io, "# ", PLUTO_VERSION_STR)
    # Anything between the version string and the first UUID delimiter will be ignored by the notebook loader.
    println(io, "")
    println(io, "using Markdown")
    println(io, "using InteractiveUtils")
    # Super Advanced Code Analysis™ to add the @bind macro to the saved file if it's used somewhere.
    if any(occursin("@bind", c.code) for c in notebook.cells)
        println(io, "")
        println(io, "# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).")
        println(io, PlutoRunner.fake_bind)
    end
    println(io)

    # TODO: this can be optimised by caching the topological order:
    # maintain cache with ordered UUIDs
    # whenever a run_reactive! is done, move the found cells **down** until they are in one group, and order them topologically within that group. Errable cells go to the bottom.

    # the next call took 2ms for a small-medium sized notebook: (so not too bad)
    # 15 ms for a massive notebook - 120 cells, 800 lines
    notebook_topo_order = topological_order(notebook, notebook.topology, notebook.cells)

    cells_ordered = union(notebook_topo_order.runnable, keys(notebook_topo_order.errable))

    for c in cells_ordered
        println(io, _cell_id_delimiter, string(c.cell_id))
        print(io, c.code)
        print(io, _cell_suffix)
    end

    println(io, _cell_id_delimiter, "Cell order:")
    for c in notebook.cells
        delim = c.code_folded ? _order_delimiter_folded : _order_delimiter
        println(io, delim, string(c.cell_id))
    end
    notebook
end

function save_notebook(notebook::Notebook, path::String)
    open(path, "w") do io
        save_notebook(io, notebook)
    end
end

save_notebook(notebook::Notebook) = save_notebook(notebook, notebook.path)

"Load a notebook without saving it or creating a backup; returns a `Notebook`. REMEMBER TO CHANGE THE NOTEBOOK PATH after loading it to prevent it from autosaving and overwriting the original file."
function load_notebook_nobackup(io, path)::Notebook
    firstline = String(readline(io))

    if firstline != _notebook_header
        error("File is not a Pluto.jl notebook")
    end

    file_VERSION_STR = readline(io)[3:end]
    if file_VERSION_STR != PLUTO_VERSION_STR
        # @info "Loading a notebook saved with Pluto $(file_VERSION_STR). This is Pluto $(PLUTO_VERSION_STR)."
    end

    collected_cells = Dict()

    # ignore first bits of file
    readuntil(io, _cell_id_delimiter)

    last_read = ""
    while !eof(io)
        cell_id_str = String(readline(io))
        if cell_id_str == "Cell order:"
            break
        else
            cell_id = UUID(cell_id_str)
            code_raw = String(readuntil(io, _cell_id_delimiter))
            # change Windows line endings to Linux
            code_normalised = replace(code_raw, "\r\n" => "\n")
            # remove the cell appendix
            code = code_normalised[1:prevind(code_normalised, end, length(_cell_suffix))]

            read_cell = Cell(cell_id, code)
            collected_cells[cell_id] = read_cell
        end
    end

    ordered_cells = Cell[]
    while !eof(io)
        cell_id_str = String(readline(io))
        o, c = startswith(cell_id_str, _order_delimiter),
        if length(cell_id_str) >= 36
            cell_id = let
                UUID(cell_id_str[end - 35:end])
            end
            next_cell = collected_cells[cell_id]
            next_cell.code_folded = startswith(cell_id_str, _order_delimiter_folded)
            push!(ordered_cells, next_cell)
        end
    end

    Notebook(ordered_cells, path)
end

function load_notebook_nobackup(path::String)::Notebook
    local loaded
    open(path, "r") do io
        loaded = load_notebook_nobackup(io, path)
    end
    loaded
end

"Create a backup of the given file, load the file as a .jl Pluto notebook, save the loaded notebook, compare the two files, and delete the backup of the newly saved file is equal to the backup."
function load_notebook(path::String, run_notebook_on_load::Bool=true)::Notebook
    backup_path = numbered_until_new(path; sep=".backup", suffix="", create_file=false)
    # local backup_num = 1
    # backup_path = path
    # while isfile(backup_path)
    #     backup_path = path * ".backup" * string(backup_num)
    #     backup_num += 1
    # end
    readwrite(path, backup_path)

    loaded = load_notebook_nobackup(path)
    # Analyze cells so that the initial save is in topological order
    update_caches!(loaded, loaded.cells)
    loaded.topology = updated_topology(loaded.topology, loaded, loaded.cells)
    save_notebook(loaded)
    # Clear symstates if autorun/autofun is disabled. Otherwise running a single cell for the first time will also run downstream cells.
    if run_notebook_on_load
        loaded.topology = NotebookTopology()
    end

    if only_versions_or_lineorder_differ(path, backup_path)
        rm(backup_path)
    else
        @warn "Old Pluto notebook might not have loaded correctly. Backup saved to: " backup_path
    end

    loaded
end

"""
Check if two savefiles are identical, up to their version numbers and a possible line shuffle.

If a notebook has not yet had all of its cells analysed, we can't deduce the topological cell order. (but can we ever??) (no)
"""
function only_versions_or_lineorder_differ(pathA::AbstractString, pathB::AbstractString)::Bool
    Set(readlines(pathA)[3:end]) == Set(readlines(pathB)[3:end])
end

function only_versions_differ(pathA::AbstractString, pathB::AbstractString)::Bool
    readlines(pathA)[3:end] == readlines(pathB)[3:end]
end

"Set `notebook.path` to the new value, save the notebook, verify file integrity, and if all OK, delete the old savefile. Normalizes the given path to make it absolute. Moving is always hard. 😢"
function move_notebook!(notebook::Notebook, newpath::String)
    # Will throw exception and return if anything goes wrong, so at least one file is guaranteed to exist.
    oldpath_tame = tamepath(notebook.path)
    newpath_tame = tamepath(newpath)
    save_notebook(notebook, oldpath_tame)
    save_notebook(notebook, newpath_tame)

    # @assert that the new file looks alright
    @assert only_versions_differ(oldpath_tame, newpath_tame)

    notebook.path = newpath_tame

    if oldpath_tame != newpath_tame
        rm(oldpath_tame)
    end
    notebook
end


text:
include("../rgi-setup.jl")
using BSON, DataStructures
function load_validation(region_str)
    #region_str ="03"
    region = parse(Int, region_str)

    # dir_output = "../output/validation/production-no-dhdt/"
    # dir_output = "../output/validation/production-dhdt/"
    #dir_output = "../output/validation-alps/production-dhdt/"

    if region in [4, 11]
        sub = "apr-20"
        dir_output = "../output/vali/validation-$(region_str)_fit-sigma/$sub"
    else
        sub = "apr-21"
        dir_output = "../output/vali/validation-$(region_str)_fit-sigma/$sub"
    end

    # calibration glaciers
    fl_fit = fit_target -> joinpath(dir_output, "reg$(region_str)_FIT_$(fit_target)")
    # validation glaciers with prior transfer
    fl_test1 = prior_choice -> joinpath(dir_output, "reg$(region_str)_PRIOR_$(prior_choice)-testgls")
    # test glaciers with prior transfer
    fl_test1_fitgls = prior_choice -> joinpath(dir_output, "reg$(region_str)_PRIOR_$(prior_choice)-fitgls")
    # validation glaciers with prior transfer and fit to IV
    fl_test2 = prior_choice -> joinpath(dir_output, "reg$(region_str)_PRIOR_$(prior_choice)_FIT_iv")

    fit_targets = [BM.FitTarget.h, BM.FitTarget.h_iv, BM.FitTarget.iv, BM.FitTarget.length]
    prior_choice = [:base, :h, :iv, :h_iv] # base means no prior transfer

    out_fit = OrderedDict()
    for ft in fit_targets
        tmp = BSON.load(fl_fit(ft)*".bson")[:out_fit]
        out_fit[Symbol(ft)] = tmp
    end

    out_test = OrderedDict()
    for p in prior_choice
        tmp = BSON.load(fl_test1(p)*".bson")[:out_test]
        out_test[(:length, p)] = tmp
    end
    for p in [:base, :h_iv]
        tmp = BSON.load(fl_test2(p)*".bson")[:out_test]
        out_test[(:iv, p)] = tmp
    end

    out_test_fitgls = OrderedDict()
    for p in prior_choice
        tmp = BSON.load(fl_test1_fitgls(p)*".bson")[:out_test]
        out_test_fitgls[(:length, p)] = tmp
    end

    return out_fit, out_test, out_test_fitgls
end


text:
module NearestNeighborDescent

using DataStructures
using Distances
using LightGraphs
using Reexport

include("knn_graph/KNNGraphs.jl")
@reexport using .KNNGraphs

include("utils.jl")
include("descent.jl")
include("search.jl")
export nndescent, search, local_join!

end # module


text:
T = 0.0

"""
    LanguageServerInstance(pipe_in, pipe_out, debug=false, env="", depot="")

Construct an instance of the language server.

Once the instance is `run`, it will read JSON-RPC from `pipe_out` and
write JSON-RPC from `pipe_in` according to the [language server
specification](https://microsoft.github.io/language-server-protocol/specifications/specification-3-14/).
For normal usage, the language server can be instantiated with
`LanguageServerInstance(stdin, stdout, false, "/path/to/environment")`.

# Arguments
- `pipe_in::IO`: Pipe to read JSON-RPC from.
- `pipe_out::IO`: Pipe to write JSON-RPC to.
- `debug::Bool`: Whether to log debugging information with `Base.CoreLogging`.
- `env::String`: Path to the
  [environment](https://docs.julialang.org/en/v1.2/manual/code-loading/#Environments-1)
  for which the language server is running. An empty string uses julia's
  default environment.
- `depot::String`: Sets the
  [`JULIA_DEPOT_PATH`](https://docs.julialang.org/en/v1.2/manual/environment-variables/#JULIA_DEPOT_PATH-1)
  where the language server looks for packages required in `env`.
"""
mutable struct LanguageServerInstance
    pipe_in
    pipe_out

    workspaceFolders::Set{String}
    documents::Dict{URI2,Document}

    debug_mode::Bool
    runlinter::Bool
    ignorelist::Set{String}
    isrunning::Bool

    packages::Dict
    env_path::String
    depot_path::String
    symbol_server::Union{Nothing,SymbolServer.SymbolServerProcess}
    format_options::DocumentFormat.FormatOptions

    function LanguageServerInstance(pipe_in, pipe_out, debug_mode::Bool = false, env_path = "", depot_path = "", packages = Dict())
        new(pipe_in, pipe_out, Set{String}(), Dict{URI2,Document}(), debug_mode, true, Set{String}(), false, packages, env_path, depot_path, nothing, DocumentFormat.FormatOptions())
    end
end

function send(message, server)
    message_json = JSON.json(message)

    write_transport_layer(server.pipe_out, message_json, server.debug_mode)
end

"""
    run(server::LanguageServerInstance)

Run the language `server`.
"""
function Base.run(server::LanguageServerInstance)
    server.symbol_server = SymbolServer.SymbolServerProcess(depot = server.depot_path, environment=server.env_path)

    @info "Started symbol server"
    server.packages = SymbolServer.getstore(server.symbol_server)
    @info "store set"
    kill(server.symbol_server)
    global T
    while true
        message = read_transport_layer(server.pipe_in, server.debug_mode)
        message_dict = JSON.parse(message)
        # For now just ignore response messages
        if haskey(message_dict, "method")
            server.debug_mode && (T = time())
            request = parse(JSONRPC.Request, message_dict)
            # server.isrunning && serverbusy(server)
            process(request, server)
            # server.isrunning && serverready(server)
        elseif get(message_dict, "id", 0)  == -100 && haskey(message_dict, "result")
            # set format options
            if length(message_dict["result"]) == length(fieldnames(DocumentFormat.FormatOptions))
                try
                    server.format_options = DocumentFormat.FormatOptions(message_dict["result"]...)
                catch
                end
            end
        end
    end
end

function serverbusy(server)
    write_transport_layer(server.pipe_out, JSON.json(Dict("jsonrpc" => "2.0", "method" => "window/setStatusBusy")), server.debug_mode)
end

function serverready(server)
    write_transport_layer(server.pipe_out, JSON.json(Dict("jsonrpc" => "2.0", "method" => "window/setStatusReady")), server.debug_mode)
end

function read_transport_layer(stream, debug_mode = false)
    header = String[]
    line = chomp(readline(stream))
    while length(line) > 0
        push!(header, line)
        line = chomp(readline(stream))
    end
    header_dict = Dict{String,String}()
    for h in header
        h_parts = split(h, ":")
        header_dict[chomp(h_parts[1])] = chomp(h_parts[2])
    end
    message_length = parse(Int, header_dict["Content-Length"])

    message = read(stream, message_length)
    message_str = String(message)
    debug_mode && @info "RECEIVED: $message_str"
    debug_mode && @info ""
    return message_str
end

function write_transport_layer(stream, response, debug_mode = false)
    global T
    response_utf8 = transcode(UInt8, response)
    n = length(response_utf8)
    write(stream, "Content-Length: $n\r\n\r\n")
    write(stream, response_utf8)
    debug_mode && @info "SENT: $response"
    debug_mode && @info string("TIME:", round(time()-T, sigdigits = 2))
end


text:
"""
Code for estimating the radial velocity based on projecting the flux onto the derivative of a template spectrum.

Author: Eric Ford
"""

""" Return mean flux (averaging over observations at different times, variance weighted) based on a common set of wavelengths.
   Inputs: flux & var (2d: pixel, time)
"""
function calc_mean_spectrum(flux::AbstractArray{T1,2}, var::AbstractArray{T2,2} ) where { T1<:Real, T2<:Real }
    flux_mean = vec(sum(flux./var,dims=2)./sum(1.0./var,dims=2))
end

""" Return mean numerical derivative (dflux/dlnλ) based on a common set of wavelengths.
Inputs:
- flux (2d)
- var (2d)
- λ (1d)
- chunk_map: Array of ranges specifying how each chunk maps into indexes of output derivative
Output:
- dflux/dlnλ: (1d)
 """
function calc_mean_dfluxdlnlambda(flux::AbstractArray{T1,2}, var::AbstractArray{T1,2}, λ::AbstractArray{T3,1},
        chunk_map::AbstractArray{URT,1}) where
    { T1<:Real, T2<:Real, T3<:Real, URT<:AbstractUnitRange} #, RT<:AbstractRange }
    flux_mean = calc_mean_spectrum(flux,var)
    deriv = Array{T1,1}(undef,length(flux_mean))
    #map(c->deriv[c] .= calc_dfluxdlnlambda(flux_mean[c],λ[c]),chunk_map )
    for c in chunk_map
         calc_dfluxdlnlambda!(deriv[c],flux_mean[c],λ[c])
    end
    return deriv
end

""" Return mean numerical derivative (d²flux/dlnλ²) based on a common set of wavelengths.
Inputs:
- flux (2d)
- var (2d)
- λ (1d)
- chunk_map: Array of ranges specifying how each chunk maps into indexes of output derivative
Output:
- d²flux/dlnλ²: (1d)
 """
 function calc_mean_d2fluxdlnlambda2(flux::AbstractArray{T1,2}, var::AbstractArray{T1,2}, λ::AbstractArray{T3,1},
        chunk_map::AbstractArray{URT,1}) where
    { T1<:Real, T2<:Real, T3<:Real, URT<:AbstractUnitRange} #, RT<:AbstractRange }
    flux_mean = calc_mean_spectrum(flux,var)
    deriv2 = Array{T1,1}(undef,length(flux_mean))
    #map(c->deriv2[c] .= calc_d2fluxdlnlambda2(flux_mean[c],λ[c]),chunk_map )
    for c in chunk_map
         calc_d2fluxdlnlambda2!(deriv2[c],flux_mean[c],λ[c])
    end
    return deriv2
end

""" `calc_rvs_from_taylor_expansion( spectra )`
Inputs:
- spectra: SpectralTimeSeriesCommonWavelengths
Optional Arguments:
- mean: mean spectrum
- deriv: dmeanflux/dlnλ
- idx: range of pixel indcies to use for calculation
- equal_weight:  For now, spectra are equal weighted. (true)
Output named pair:
- rv: Vector of estimated radial velocities
- σ_rv: Vector of uncertainties in rv estimates
"""
function calc_rvs_from_taylor_expansion(spectra::STS; mean::MT = calc_mean_spectrum(spectra),
                deriv::DT = calc_mean_dfluxdlnlambda(spectra), idx::RT = 1:length(mean),
                equal_weight::Bool = true ) where {
                    STS<:AbstractSpectralTimeSeriesCommonWavelengths, T1<:Real, MT<:AbstractVector{T1},
                    T2<:Real, DT<:AbstractVector{T2}, RT<:AbstractUnitRange }
   @assert length(mean) == length(deriv)
   @assert size(spectra.flux,1) == length(mean)
   @assert minimum(idx) >=1
   @assert maximum(idx) <= length(mean)

   if equal_weight # Pixels equally-weighted
      norm = sum(abs2.(deriv[idx]))
      rv = sum((spectra.flux[idx,:].-mean[idx]).*deriv[idx],dims=1).*(speed_of_light_mps/norm)
      # TODO: WARN: Uncertinaties ignores correlations between pixels, particularly problematic when oversample pixels
      σ_rv = sqrt.(sum(spectra.var[idx,:].*abs2.(deriv[idx]),dims=1)) .*(speed_of_light_mps/norm)
  else # Pixels inverse variance weighted
      # TODO:  CHECK/FIX? math for inverse variance weighting.
      # Or maybe the math is right, but it's just a bad idea to have different weighting within one line/chunk
      @info "Warning: I think this is either wrong or a bad idea.  Need to check."
      norm = sum(abs2.(deriv[idx])./spectra.var[idx,:])
      rv = sum((spectra.flux[idx,:].-mean[idx]).*deriv[idx]./spectra.var[idx,:],dims=1).*(speed_of_light_mps/norm)
      σ_rv = sqrt.(sum(abs2.(deriv)./spectra.var[idx,:],dims=1)).*(speed_of_light_mps/norm)
   end
   return (rv=-vec(rv), σ_rv=vec(σ_rv))
end

""" `calc_rvs_from_taylor_expansion_alt( spectra )`
Experimental version of [calc_rvs_from_taylor_expansion](@ref).
"""
function calc_rvs_from_taylor_expansion_alt(spectra::STS; mean::MT = calc_mean_spectrum(spectra),
                deriv::DT = calc_mean_dfluxdlnlambda(spectra),
                deriv2::DT = calc_d2fluxdlnlambda2(spectra), idx::RT = 1:length(mean),
                equal_weight::Bool = true ) where {
                    STS<:AbstractSpectralTimeSeriesCommonWavelengths, T1<:Real, MT<:AbstractVector{T1},
                    T2<:Real, DT<:AbstractVector{T2}, RT<:AbstractUnitRange }
   @assert length(mean) == length(deriv)
   @assert size(spectra.flux,1) == length(mean)
   @assert minimum(idx) >=1
   @assert maximum(idx) <= length(mean)

   if equal_weight # Pixels equally-weighted
      #norm = sum(abs2.(deriv[idx]))
      norm = sum(deriv[idx].^2 .- (spectra.flux[idx,:].-mean[idx]).*deriv2[idx],dims=1)
      rv = sum(((spectra.flux[idx,:].-mean[idx]).*deriv[idx]),dims=1)./norm
      rv *= speed_of_light_mps
      # TODO: WARN: Uncertinaties ignores correlations between pixels, particularly problematic when oversample pixels
      σ_rv = sqrt.(sum(spectra.var[idx,:].*abs2.(deriv[idx]),dims=1)./norm)
      σ_rv *= speed_of_light_mps
  else # Pixels inverse variance weighted
      # TODO:  CHECK/FIX? math for inverse variance weighting.
      # Or maybe the math is right, but it's just a bad idea to have different weighting within one line/chunk
      @info "Warning: I think this is either wrong or a bad idea.  Need to check."
      norm = sum(abs2.(deriv[idx]./spectra.var[idx,:]))  # TODO comparable to norm_old, if new norm works, then check updating this also helps
      rv = sum((spectra.flux[idx,:].-mean[idx]).*deriv[idx]./spectra.var[idx,:],dims=1).*(speed_of_light_mps/norm)
      σ_rv = sqrt.(sum(abs2.(deriv)./spectra.var[idx,:],dims=1)).*(speed_of_light_mps/norm)
   end
   return (rv=-vec(rv), σ_rv=vec(σ_rv))
end

""" `calc_chunk_rvs_from_taylor_expansion( spectra )`
Inputs:
- spectra: SpectralTimeSeriesCommonWavelengths
Optional Arguments:
- mean: mean spectrum
- deriv: dmeanflux/dlnλ
- equal_weight:  For now, spectra are equal weighted. (true)
Output named pair:
- rv: Vector of vector of estimated radial velocities for each chunk
- σ_rv: Vector of vector of uncertainties in rv estimates for each chunk
"""
function calc_chunk_rvs_from_taylor_expansion(spectra::STS; mean::MT = calc_mean_spectrum(spectra),
                deriv::DT = calc_mean_dfluxdlnlambda(spectra),
                equal_weight::Bool = false ) where {
                    STS<:AbstractSpectralTimeSeriesCommonWavelengths, T1<:Real, MT<:AbstractVector{T1},
                    T2<:Real, DT<:AbstractVector{T2}, RT<:AbstractUnitRange }
   @assert length(mean) == length(deriv)
   @assert size(spectra.flux,1) == length(mean)

   map(idx->calc_rvs_from_taylor_expansion(spectra,mean=mean,deriv=deriv,idx=idx),spectra.chunk_map)
end


text:
using Plexifications
using Base.Test: @test, @test_throws

@test_throws ErrorException begin
    inv(HyperComplex(0, 0.0))
end

@test begin
    iszerodivisor(HyperComplex(0,0,1,2))
end

@test begin
    !iszerodivisor(HyperComplex(1,2,3,5))
end

@test begin
    a = Plexifications.randomBigFloat()
    isreal(HyperComplex(a))
end

@test begin
    a = 1
    !isreal(HyperComplex(a, a))
end

@test begin
    z = random(HyperComplex{Int})
    z == +(z)
end

@test begin
    a = Plexifications.randomBigFloat()
    real(HyperComplex(a)) == a
end

@test begin
    io = IOBuffer()
    show(io, HyperComplex(1, 2, 3, 4, 5, 6, 7, 8))
    l = String(take!(io))
    r = "[1: 1, H: 2, A: 3, HA: 4, B: 5, HB: 6, AB: 7, HAB: 8]"
    l == r
end

@test begin
    x = random(HyperComplex{BigFloat})
    y = random(HyperComplex{BigFloat})
    x != y
end

@test begin
    a = Plexifications.randomBigInt()
    z = random(HyperComplex{BigInt})
    l = a + z
    r = z + a
    l == r
end

@test begin
    a = Plexifications.randomBigInt()
    z = random(HyperComplex{BigInt})
    l = a - z
    r = -(z - a)
    l == r
end

@test begin
    x = random(HyperComplex{BigInt})
    y = random(HyperComplex{BigInt})
    l = x + y
    r = y + x
    l == r
end

@test begin
    x = random(DualComplex{BigInt})
    y = random(HyperComplex{BigInt})
    l = x + y
    r = y + x
    l == r
end

@test begin
    x = random(HyperComplex{BigInt})
    y = random(HyperComplex{BigInt})
    l = x * y
    r = y * x
    l == r
end

@test begin
    x = random(DualComplex{BigInt})
    y = random(HyperComplex{BigInt})
    l = x * y
    r = y * x
    l == r
end

@test begin
    x = random(HyperComplex{BigInt})
    y = random(HyperComplex{BigInt})
    l = x - y
    r = -(y - x)
    l == r
end

@test begin
    x = random(DualComplex{BigInt})
    y = random(HyperComplex{BigInt})
    l = x - y
    r = -(y - x)
    l == r
end

@test begin
    x = random(HyperComplex{BigInt})
    y = random(HyperComplex{BigInt})
    z = random(HyperComplex{BigInt})
    l = (x + y) + z
    r = x + (y + z)
    l == r
end

@test begin
    x = random(HyperComplex{BigInt})
    y = random(HyperComplex{BigInt})
    z = random(HyperComplex{BigInt})
    l = (x * y) * z
    r = x * (y * z)
    l == r
end

@test begin
    x = random(HyperComplex{BigInt})
    barA(barA(x)) == x
end

@test begin
    x = random(HyperComplex{BigInt})
    barB(barB(x)) == x
end

@test begin
    x = random(HyperComplex{BigInt})
    barH(barH(x)) == x
end

@test begin
    x = random(HyperComplex{BigInt})
    barB(barA(x)) == barA(barB(x))
end

@test begin
    x = random(HyperComplex{BigInt})
    barH(barA(x)) == barA(barH(x))
end

@test begin
    x = random(HyperComplex{BigInt})
    barB(barH(x)) == barH(barB(x))
end

@test begin
    x = random(HyperComplex{Rational{BigInt}})
    inv(inv(x)) == x
end

@test begin
    x = random(HyperComplex{BigInt})
    y = random(HyperComplex{BigInt})
    l = barA(x * y)
    r = barA(x) * barA(y)
    l == r
end

@test begin
    x = random(HyperComplex{BigInt})
    y = random(HyperComplex{BigInt})
    l = barA(x + y)
    r = barA(x) + barA(y)
    l == r
end

@test begin
    x = random(HyperComplex{BigInt})
    y = random(HyperComplex{BigInt})
    l = barA(x - y)
    r = barA(x) - barA(y)
    l == r
end

@test begin
    x = random(HyperComplex{BigInt})
    y = random(HyperComplex{BigInt})
    l = barB(x * y)
    r = barB(x) * barB(y)
    l == r
end

@test begin
    x = random(HyperComplex{BigInt})
    y = random(HyperComplex{BigInt})
    l = barB(x + y)
    r = barB(x) + barB(y)
    l == r
end

@test begin
    x = random(HyperComplex{BigInt})
    y = random(HyperComplex{BigInt})
    l = barB(x - y)
    r = barB(x) - barB(y)
    l == r
end

@test begin
    x = random(HyperComplex{BigInt})
    y = random(HyperComplex{BigInt})
    l = barH(x * y)
    r = barH(x) * barH(y)
    l == r
end

@test begin
    x = random(HyperComplex{BigInt})
    y = random(HyperComplex{BigInt})
    l = barH(x + y)
    r = barH(x) + barH(y)
    l == r
end

@test begin
    x = random(HyperComplex{BigInt})
    y = random(HyperComplex{BigInt})
    l = barH(x - y)
    r = barH(x) - barH(y)
    l == r
end

@test begin
    x = random(HyperComplex{Rational{BigInt}})
    y = random(HyperComplex{Rational{BigInt}})
    l = inv(x * y)
    r = inv(y) * inv(x)
    l == r
end

@test begin
    x = random(HyperComplex{Rational{BigInt}})
    y = random(HyperComplex{Rational{BigInt}})
    l = inv(x / y)
    r = y / x
    l == r
end

@test begin
    x = random(HyperComplex{Rational{BigInt}})
    y = random(HyperComplex{Rational{BigInt}})
    l = inv(x \ y)
    r = y \ x
    l == r
end

@test begin
    x = random(HyperComplex{BigInt})
    y = random(HyperComplex{BigInt})
    z = random(HyperComplex{BigInt})
    l = x * (y + z)
    r = (x * y) + (x * z)
    l == r
end

@test begin
    a = Plexifications.randomBigInt()
    x = random(HyperComplex{BigInt})
    y = random(HyperComplex{BigInt})
    l = a * (x + y)
    r = (a * x) + (a * y)
    l == r
end

@test begin
    a = Plexifications.randomBigInt()
    x = random(HyperComplex{BigInt})
    l = a * x
    r = x * a
    l == r
end

@test begin
    x = random(HyperComplex{BigInt})
    y = random(HyperComplex{BigInt})
    z = random(HyperComplex{BigInt})
    l = (x + y) * z
    r = (x * z) + (y * z)
    l == r
end

@test begin
    x = random(HyperComplex{BigInt})
    y = random(HyperComplex{BigInt})
    z = random(HyperComplex{BigInt})
    l = x * (y - z)
    r = (x * y) - (x * z)
    l == r
end

@test begin
    x = random(HyperComplex{BigInt})
    y = random(HyperComplex{BigInt})
    z = random(HyperComplex{BigInt})
    l = (x - y) * z
    r = (x * z) - (y * z)
    l == r
end

@test begin
    x = random(HyperComplex{BigInt})
    n = norm(x)
    n > big(0) || n == big(0)
end

@test begin
    x = random(HyperComplex{BigInt})
    y = random(HyperComplex{BigInt})
    l = norm(x * y)
    r = norm(x) * norm(y)
    l == r
end

@test_throws ErrorException begin
    inv(HyperComplex(0, 0.0))
end

@test_throws ErrorException begin
    random(HyperComplex{Int}) / 0
end

@test_throws ErrorException begin
    0 \ random(HyperComplex{Int})
end

@test begin
    a = Plexifications.randomBigRational()
    z = random(HyperComplex{Rational{BigInt}})
    l = z / a
    r = a \ z
    l == r
end

@test begin
    a = Plexifications.randomBigRational()
    z = random(HyperComplex{Rational{BigInt}})
    l = a / z
    r = z \ a
    l == r
end

@test begin
    z = random(HyperComplex{BigInt})
    z == z + zero(HyperComplex{BigInt})
end

@test begin
    z = random(HyperComplex{BigInt})
    z == z * one(HyperComplex{BigInt})
end

@test begin
    z = random(HyperComplex{BigInt})
    z == z + zero(z)
end

@test begin
    z = random(HyperComplex{BigInt})
    z == z * one(z)
end

@test begin
    x = random(HyperComplex{BigInt})
    y = random(HyperComplex{BigInt})
    l = commutator(x, y)
    r = zero(HyperComplex{BigInt})
    l == r
end

@test begin
    x = random(HyperComplex{BigInt})
    y = random(HyperComplex{BigInt})
    l = alternatorL(x, y)
    r = zero(HyperComplex{BigInt})
    l == r
end

@test begin
    x = random(HyperComplex{BigInt})
    y = random(HyperComplex{BigInt})
    l = alternatorR(x, y)
    r = zero(HyperComplex{BigInt})
    l == r
end

@test begin
    x = random(HyperComplex{BigInt})
    y = random(HyperComplex{BigInt})
    l = flexator(x, y)
    r = zero(HyperComplex{BigInt})
    l == r
end

@test begin
    x = random(HyperComplex{BigInt})
    y = random(HyperComplex{BigInt})
    z = random(HyperComplex{BigInt})
    l = associator(x, y, z)
    r = zero(HyperComplex{BigInt})
    l == r
end

@test begin
    x = random(HyperComplex{BigInt})
    y = random(HyperComplex{BigInt})
    z = random(HyperComplex{BigInt})
    l = jacobiator(x, y, z)
    r = 3 * x * y * z
    l == r
end

@test begin
    l = HyperComplex(Cmplex(1,2))
    r = HyperComplex(1,2)
    l == r
end

@test begin
    l = HyperComplex(DualComplex(1,2))
    r = HyperComplex(1,2,0)
    l == r
end

@test begin
    l = HyperComplex(1, 2, 3.0)
    r = HyperComplex(1.0, 2.0, 3.0)
    l == r
end

@test begin
    l = HyperComplex(1, 2, 3, 4.0)
    r = HyperComplex(1.0, 2.0, 3.0, 4.0)
    l == r
end

@test begin
    l = HyperComplex(1, 2, 3, 4, 5.0)
    r = HyperComplex(1.0, 2.0, 3.0, 4.0, 5.0)
    l == r
end

@test begin
    l = HyperComplex(1, 2, 3, 4, 5, 6.0)
    r = HyperComplex(1.0, 2.0, 3.0, 4.0, 5.0, 6.0)
    l == r
end

@test begin
    l = HyperComplex(1, 2, 3, 4, 5, 6, 7.0)
    r = HyperComplex(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0)
    l == r
end

@test begin
    l = HyperComplex(1, 2, 3, 4, 5, 6, 7, 8.0)
    r = HyperComplex(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0)
    l == r
end


text:
using Plots
using CSV, DataFrames
using Statistics
# %% Read in the fitting results of the Synchrotron Index (βs), Dust Spectral Index (βd) and Dust Temperature (Td)
milca_so = CSV.read("sosat/beta_results.csv", DataFrame)
println("SO-SAT only:")
println("- βs = ", mean(milca_so.beta_s), " ± ", std(milca_so.beta_s))
println("- βd = ", mean(milca_so.beta_d), " ± ", std(milca_so.beta_d))
println("- Td = ", mean(milca_so.T_d), " ± ", std(milca_so.T_d))
milca = CSV.read("sosatfyst/beta_results.csv", DataFrame)
println("SO-SAT+FYST:")
println("- βs = ", mean(milca.beta_s), " ± ", std(milca.beta_s))
println("- βd = ", mean(milca.beta_d), " ± ", std(milca.beta_d))
println("- Td = ", mean(milca.T_d), " ± ", std(milca.T_d))
# %% Scatter Plot of βd and Td. Save to "figure6.pdf"
p = scatter(
    milca_so.beta_d,
    milca_so.T_d,
    xlab = "Dust Spectral Index, βd",
    ylab = "Dust Temperature, Td [K]",
    lab = "SO-SAT only",
    legend = :best,
    xlims = [1.4, 1.9],
    markersize = 3,
    legendfontsize = 15,
    tickfontsize = 15,
    labelfontsize = 18,
)
p = scatter!(milca.beta_d, milca.T_d, lab = "SO-SAT+FYST")
savefig("figure6.pdf")
display(p)


text:
using SantasLittleHelpers
using Test
using Statistics, StatsBase, Dates, DataFrames

@testset "Simple Functions" begin
    x = rand(100)
    @test loggrowth(x)[end] ≈ 100*diff(log.(x))[end]
    @test autocorrelate(x) ≈ first(autocor(x,[1]))
    @test autocorrelate(abs,x) ≈ first(autocor(abs.(x),[1]))
    @test logdiff(x[1:4]) ≈ log(x[4])-log(x[1])

    @test monthtoquarter(12) == 4
    @test yrqtrfun(Date(2018,1,1)) == "2018Q1"

    @test simiterator(10,3)[end] == 4
    @test perioditerator(10,3)[end] == 1
end

@testset "Harder Functions" begin
    x = rand(100)

    varx = var(x)
    x2 = [sum(x[i-1:i]) for i in 2:100]
    testvr2 = var(x2) / (2*varx)

    @test varianceratio(x,1) == 1
    @test varianceratio(x,2) ≈ vr2(x)
    @test vr2(x) ≈ testvr2

    @test applyrolling(sum,x,-1:0,missing)[end] ≈ x2[end]
    @test applyrolling(makekernel(sum,-1:0),x,missing)[end] ≈ x2[end]
end

@testset "DataFrames" begin
    df = DataFrame(x=[1,2,NaN,3],y=[NaN,2,3,4])
    nantomissing!(df,:x)
    @test ismissing(df.x[3])
    nantomissing!(df)
    @test ismissing(df.y[1])
end


text:
function _convert_step_to_docker_commands(template::Template,
                                          step::PlaceholderDockerCopyFiles)::Vector{String}
    commands = String[]
    for i = 1:length(template.files.list)
        file = template.files.list[i]
        outside_src_full_path = "simplecontainergenerator_container_files/file_$(i)/$(file.dest_filename)"
        inside_dest_full_path = file.dest_full_path
        inside_dest_dir = file.dest_dir
        push!(commands, "RUN mkdir -p $(inside_dest_dir)")
        push!(commands, "COPY $(outside_src_full_path) $(inside_dest_full_path)")
    end
    return commands
end

function _convert_step_to_docker_commands(template::Template,
                                          step::DockerOnlyLine)::Vector{String}
    return String[step.full_contents]
end

function _convert_step_to_docker_commands(template::Template,
                                          step::RunStep)::Vector{String}
    return String["RUN $(step.command)"]
end

function _convert_step_to_docker_commands(template::Template,
                                          step::EnvStep)::Vector{String}
    return String["ENV $(step.name)=$(step.value)"]
end

# function _convert_step_to_docker_commands(template::Template,
#                                           step::AptStep)::Vector{String}
#     return String["RUN $(step.command)"]
# end

function _convert_step_to_docker_commands(template::Template,
                                          step::PlaceholderYumStep)::Vector{String}
    return String[]
end

function _docker_expand_template(template::Template)::String
    dockerfile_lines = String[]
    for step in template.steps.list
        lines_to_add::Vector{String} = _convert_step_to_docker_commands(template,
                                                                        step)::Vector{String}
        append!(dockerfile_lines, lines_to_add)
    end
    dockerfile_content = join(dockerfile_lines, "\n") * "\n"
    return dockerfile_content
end

function _write_all_docker_files(template::Template;
                                 output_directory::AbstractString)
    dockerfile_output_path = joinpath(output_directory, "Dockerfile")
    simplecontainergenerator_container_files_output_path = joinpath(output_directory, "simplecontainergenerator_container_files")

    dockerfile_content = _docker_expand_template(template)

    rm(dockerfile_output_path;
       force = true,
       recursive = true)
    mkpath(output_directory)
    open(dockerfile_output_path, "w") do io
        println(io, dockerfile_content)
    end
    @debug("Wrote file: \"$(dockerfile_output_path)\"")

    rm(simplecontainergenerator_container_files_output_path;
       force = true,
       recursive = true)
    mkpath(simplecontainergenerator_container_files_output_path)
    for i = 1:length(template.files.list)
        file = template.files.list[i]
        file_contents = file.contents
        file_output_path = joinpath(simplecontainergenerator_container_files_output_path,
                                    "file_$(i)",
                                    "$(file.dest_filename)")
        mkpath(dirname(file_output_path))
        open(file_output_path, "w") do io
            println(io, file_contents)
        end
        @debug("Wrote file: \"$(file_output_path)\"")
    end

    return output_directory
end


text:
const TEST_NAME = "Xorshift1024Star"

include("common.jl")

using RandomNumbers.Xorshifts

r = Xorshift1024Star(123)

test_all(r, 100_000_000)


text:
# Precompile
# ==========

# Parser
# ------

if VERSION < v"0.5-"
    precompile(Base.open, (ASCIIString, Type{Seq.FASTA},))
    precompile(Base.open, (ASCIIString, Type{Seq.FASTQ}, Type{Seq.QualityEncoding},))
    precompile(Base.open, (ASCIIString, Type{Intervals.BED},))
else
    precompile(Base.open, (String, Type{Seq.FASTA},))
    precompile(Base.open, (String, Type{Seq.FASTQ}, Type{Seq.QualityEncoding},))
    precompile(Base.open, (String, Type{Intervals.BED},))
end
precompile(Base.read, (Seq.FASTAParser{Seq.BioSequence},))
precompile(Base.read, (Seq.FASTAParser{Seq.DNASequence},))
precompile(Base.read, (Seq.FASTAParser{Seq.RNASequence},))
precompile(Base.read, (Seq.FASTAParser{Seq.AminoAcidSequence},))
precompile(Base.read, (Seq.FASTQParser{Seq.DNASequence},))
precompile(Base.read, (Intervals.BEDParser,))


text:
include("../../src/evaluation/diagnostics.jl")
using Test

data_actual = [1,1,0,0,1,0,0,0,1,1,1]
data_pred =   [1,1,1,0,1,1,1,0,0,0,1]
data_proba = [0.9, 0.8, 0.52, 0.2, 0.55, 0.52, 0.6, 0.1, 0.45, 0.3, 0.8]

@testset "diagnostics" begin
    @test true_positives(data_actual, data_pred) == 4.0
    @test round(true_positives(data_actual, data_pred, rate=true), digits=4) == 0.6667

    @test true_negatives(data_actual, data_pred) == 2.0
    @test round(true_negatives(data_actual, data_pred, rate=true), digits=4) == 0.3333

    @test false_positives(data_actual, data_pred) == 3.0
    @test false_positives(data_actual, data_pred, rate=true) == 0.6

    @test false_negatives(data_actual, data_pred) == 2.0
    @test round(false_negatives(data_actual, data_pred, rate=true), digits=4) == 0.3333

    @test round(precision_classification(data_actual, data_pred), digits=4) == 0.5714
    @test recall(data_actual, data_pred) == 0.5
    @test round(accuracy(data_actual, data_pred), digits=4) == 0.5455
    @test round(f1_score(data_actual, data_pred), digits=4) == 0.5333
    @test confusion_matrix(data_actual, data_pred) == [4.0 3.0; 2.0 2.0]
    @test round.(confusion_matrix(data_actual, data_pred, ratio=true), digits=4) == [0.3636 0.2727; 0.1818 0.1818]
    @test predict_from_threshold(data_proba, 0.5) == data_pred

    @test best_f1_search(data_actual, data_proba)[1:2] == (0.2, 0.8)
end


text:
### replace this with SQL later.
# function getexperimentpath(tag::String)
#     repo_path = pwd()
#     base_path = ""
#     experiment_name = ""
#
#     if tag[1:3] == "NRC"
#         base_path = "experiments/NRC"
#
#         if tag == "NRC_Glucose_2018"
#             experiment_name = "glucose/Sep-25-2018"
#
#         elseif tag == "NRC_DMEM_2012"
#             experiment_name = "dmem_medium/Oct-22-2012"
#
#         elseif tag ==
#     end
#
#     return joinpath(repo_path, base_path, experiment_name)
# end

# function getpath(base_path::String, experiment_name::String)::String
#     #repo_path = pwd()
#     repo_path = dirname(pathof(@__MODULE__))

#     return joinpath(repo_path, base_path, experiment_name)
# end

"""
Currently available options for folder_name are: BMRB-500-0.5mM, BMRB-500-2mM, and BMRB-700.
"""
function getBMRBexperimentpaths(folder_name::String; join_flag = true)

    src_path = dirname(pathof(NMRData))
    exp_folder = "../combination/$(folder_name)"

    return readdir(joinpath(src_path, exp_folder); join = join_flag)
end

"""
Currently available options for folder_name are: BMRB-500-0.5mM, BMRB-500-2mM, and BMRB-700.
Run getBMRBexperimentpaths(folder_name) for a list of available compound names.
"""
function getBMRBexperimentpath(folder_name::String, name::String; join_flag = true)

    paths = getBMRBexperimentpaths(folder_name::String; join_flag = join_flag)
    filter!(xx->occursin(name,xx), paths)

    return paths
end


text:
"""
    $(TYPEDEF)
An abstract type for all VLBI interfermetry data types. See [EHTVisibilityDatum](@ref) for an example.
"""
abstract type AbstractInterferometryDatum{T} end

abstract type AbstractVisibilityDatum{T} <: AbstractInterferometryDatum{T} end
abstract type AbstractLinearPolDatum{S<:AbstractVisibilityDatum, T} <: AbstractInterferometryDatum{T} end
abstract type AbstractCrossPolDatum{S,T} <: AbstractInterferometryDatum{T} end

abstract type ClosureProducts{T} <: AbstractInterferometryDatum{T} end

abstract type Observation{T} end

using DelimitedFiles
using AstroTime: modified_julian

export uvpositions, stations, getdata, arrayconfig,
       getuv, baselines, rescaleuv!, scantable

"""
    $(TYPEDEF)
This defined the abstract type for an array configuration. Namely, baseline
times, SEFD's, bandwidth, observation frequencies, etc.
"""
abstract type ArrayConfiguration end

"""
    $(TYPEDEF)
Stores all the non-visibility data products for an EHT array. This is useful when evaluating
model visibilities.

# $(FIELDS)
"""
struct EHTArrayConfiguration{F,T<:AbstractArray} <: ArrayConfiguration
    """
    Observing frequency (Hz)
    """
    frequency::F
    """
    Observing bandwith (Hz)
    """
    bandwidth::F
    """
    A struct array of `ArrayBaselineDatum` holding time, u, v, baselines.
    """
    data::T
end

"""
    $(TYPEDEF)
Array config file for closure quantities. This stores the design matrix `designmat`
that transforms from visibilties to closure products.
"""
struct ClosureConfig{A,D} <: ArrayConfiguration
    ac::A
    designmat::D
function ClosureConfig(ac, dmat)
    A = typeof(ac)
    sdmat = blockdiag(sparse.(dmat)...)
    D = typeof(sdmat)
    return new{A,D}(ac, sdmat)
end
end

"""
    $(SIGNATURES)
Get the u, v positions of the array.
"""
function getuv(ac::ArrayConfiguration)
    return ac.data.u, ac.data.v
end

function getuv(ac::ClosureConfig)
    return getuv(ac.ac)
end

"""
    $(SIGNATURES)
Get the u, v, time, freq of the array as a tuple.
"""
function uvtimefreq(ac::EHTArrayConfiguration)
    u,v = getuv(ac)
    t = ac.data.time
    ν = ac.frequency
    return u, v, t, fill(ν, length(u))
end

"""
    $(TYPEDEF)
A single datum of an observing array.
"""
struct ArrayBaselineDatum{T}
    time::T
    u::T
    v::T
    baseline::Tuple{Symbol, Symbol}
end

const ArrayQuadrangleDatum = NTuple{4, ArrayBaselineDatum{T}} where {T}
const ArrayTriangleDatum = NTuple{3, ArrayBaselineDatum{T}} where {T}

"""
    uvpositions(datum)
Get the uvp positions of an inferometric datum.
"""
uvpositions(D::AbstractVisibilityDatum) = D.u, D.v

"""
    $(SIGNATURES)
The main data product type in `Comrade` this stores the `data` which can be a StructArray
of any `AbstractInterferometryDatum` type.

# ($FIELDS)
"""
Base.@kwdef struct EHTObservation{F,T<:AbstractInterferometryDatum{F},S<:StructArray{T}, A, N} <: Observation{F}
    """
    StructArray of data productts
    """
    data::S
    """
    Array config holds ancillary information about array
    """
    config::A
    """
    modified julia date of the observation
    """
    mjd::N
    """
    RA of the observation in J2000 (deg)
    """
    ra::F
    """
    DEC of the observation in J2000 (deg)
    """
    dec::F
    """
    bandwidth of the observation (Hz)
    """
    bandwidth::F
    """
    frequency of the observation (Hz)
    """
    frequency::F
    """
    Common source name
    """
    source::Symbol
    """
    Time zone used.
    """
    timetype::Symbol = :UTC
end

getdata(obs::Observation, s::Symbol) = getproperty(getfield(obs, :data), s)

# Implement the tables interface
Tables.istable(::Type{<:EHTObservation}) = true
Tables.columnaccess(::Type{<:EHTObservation}) = true
Tables.columns(t::EHTObservation) = getfield(t, :data)

Tables.getcolumn(t::EHTObservation, ::Type{T}, col::Int, nm::Symbol) where {T} = getdata(t, nm)
Tables.getcolumn(t::EHTObservation, nm::Symbol) = getdata(t, nm)
Tables.getcolumn(t::EHTObservation, i::Int) = Tables.getcolumn(t, Tables.columnames(t)[i])
Tables.columnnames(t::EHTObservation) = propertynames(getfield(t, :data))

Base.getindex(data::EHTObservation, s::Symbol) = Tables.getcolumn(data, s)
Base.getindex(data::EHTObservation, i::Int) = data.data[i]
Base.getindex(data::EHTObservation, I...) = getindex(data.data, I...)
Base.length(data::EHTObservation) = length(data.data)

"""
    $(SIGNATURES)
Get all the stations in a observation. The result is a vector of symbols.
"""
function stations(d::EHTObservation{T,A}) where {T,A<:AbstractInterferometryDatum}
    bl = getdata(d, :baseline)
    s1 = first.(bl)
    s2 = last.(bl)
    return unique([s1..., s2...])
end

function Base.show(io::IO, d::EHTObservation{F,D}) where {F,D}
    println(io, "EHTObservation{$F,$D, ...}")
    println(io, "  source: ", d.source)
    println(io, "  mjd: ", d.mjd)
    println(io, "  frequency: ", d.frequency)
    println(io, "  bandwidth: ", d.bandwidth)
    println(io, "  stations: ", stations(d))
    println(io, "  nsamples: ", length(d))
end

"""
    getdata(obs::EHTObservation, s::Symbol)
Pass-through function that gets the array of `s` from the EHTObservation. For example
say you want the times of all measurement then

```julia
getdata(obs, :time)
```
"""

"""
    $(SIGNATURES)
A struct holding the information for a single measured visibility.

# $(FIELDS)

"""
Base.@kwdef struct EHTVisibilityDatum{T<:Number} <: AbstractVisibilityDatum{T}
    """
    real component of the visibility (Jy)
    """
    visr::T
    """
    imaginary component of the visibility (Jy)
    """
    visi::T
    """
    error of the visibility (Jy)
    """
    error::T
    """
    x-direction baseline length in λ
    """
    u::T
    """
    y-direction baseline length in λ
    """
    v::T
    """
    Time of the observation in hours
    """
    time::T
    """
    station baseline codes
    """
    baseline::NTuple{2,Symbol}
end

"""
    $(SIGNATURES)
Return the complex visibility of the visibility datum
"""
@inline function visibility(D::EHTVisibilityDatum{T}) where {T}
        return Complex{T}(D.visr, D.visi)
end

"""
    $(SIGNATURES)
Get the amplitude of a visibility datum
"""
function amplitude(D::EHTVisibilityDatum)
    amp = hypot(D.visr,D.visi)
    return EHTVisibilityAmplitudeDatum(amp, D.error,
                                       D.u, D.v,
                                       D.time,
                                       D.frequency,
                                       D.bandwidth,
                                       D.baseline
                                    )
end

"""
    $(SIGNATURES)
A struct holding the information for a single measured visibility amplitude.

# $(FIELDS)

"""
Base.@kwdef struct EHTVisibilityAmplitudeDatum{T<:Number} <: AbstractVisibilityDatum{T}
    """
    amplitude (Jy)
    """
    amp::T
    """
    error of the visibility amplitude (Jy)
    """
    error::T
    """
    x-direction baseline length in λ
    """
    u::T
    """
    y-direction baseline length in λ
    """
    v::T
    """
    Time of the observation in hours
    """
    time::T
    """
    station baseline codes
    """
    baseline::NTuple{2,Symbol}
end

function checktriangle(D1::EHTVisibilityDatum,
                       D2::EHTVisibilityDatum,
                       D3::EHTVisibilityDatum)
    b1 = D1.baseline
    b2 = D2.baseline
    b3 = D3.baseline
    l = length(unique([b1..., b2..., b3...]))
    @assert l == 3 "For a valid closure phase you need 3 unique stations not $l"
    @assert (D1.time == D2.time == D3.time) "For a valid closure phase the times need to match"

end

@inline function amplitude(D::EHTVisibilityAmplitudeDatum{T}) where {T}
    return D.amp
end

@inline function bispectrum(D1::EHTVisibilityDatum, D2::EHTVisibilityDatum, D3::EHTVisibilityDatum)
    visibility(D1)*visibility(D2)*visibility(D3)
end

"""
    rescaleuv!(data)
rescale the u-v lengths according to scale. This can be useful when you want the spatial scales
to be in 1/μas instead of 1/rad
"""
function rescaleuv!(data::EHTObservation{T,D}, scale) where {T,D<:Union{EHTVisibilityAmplitudeDatum,EHTVisibilityDatum}}
    data.data.u .= data.data.u.*scale
    data.data.v .= data.data.v.*scale
    return data
end

"""
    $(TYPEDEF)

A Datum for a single closure phase.

# $(FIELDS)

"""
Base.@kwdef struct EHTClosurePhaseDatum{T<:Number} <: ClosureProducts{T}
    """
    closure phase (rad)
    """
    phase::T
    """
    error of the closure phase assuming the high-snr limit
    """
    error::T
    """
    u (λ) of first station
    """
    u1::T
    """
    v (λ) of first station
    """
    v1::T
    """
    u (λ) of second station
    """
    u2::T
    """
    v (λ) of second station
    """
    v2::T
    """
    u (λ) of third station
    """
    u3::T
    """
    v (λ) of third station
    """
    v3::T
    """
    Measured time of closure phase in hours
    """
    time::T
    """
    station baselines used
    """
    triangle::NTuple{3,Symbol}
end

function rescaleuv!(data::EHTObservation{T,D}, scale) where {T,D<:EHTClosurePhaseDatum}
    data.data.u1 .= data.data.u1.*scale
    data.data.v1 .= data.data.v1.*scale
    data.data.u2 .= data.data.u2.*scale
    data.data.v2 .= data.data.v2.*scale
    data.data.u3 .= data.data.u3.*scale
    data.data.v3 .= data.data.v3.*scale
    return data
end

function stations(d::EHTObservation{T,A}) where {T,A<:EHTClosurePhaseDatum}
    bl = getdata(d, :triangle)
    return unique(vcat(collect.(bl)...))
end

"""
    $(SIGNATURES)
Computes the closure phase of the three visibility datums.

# Notes
We currently use the high SNR Gaussian error approximation for the closure phase.
In the future we may use the moment matching from Monte Carlo sampling.
"""
function closure_phase(D1::EHTVisibilityDatum,
                      D2::EHTVisibilityDatum,
                      D3::EHTVisibilityDatum)

    checktriangle(D1,D2,D3)

    amp1 = amplitude(D1).amp
    amp2 = amplitude(D2).amp
    amp3 = amplitude(D3).amp
    u1,v1 = uvpositions(D1)
    u2,v2 = uvpositions(D2)
    u3,v3 = uvpositions(D3)

    bis = bispectrum(D1, D2, D3)
    s12 = unique([D1.baseline..., D2.baseline...])
    s123 = unique([s12..., D3.baseline...])
    #Use the Gaussian approximation TODO hook this into Measurements.jl?
    err = sqrt((D1.error/amp1)^2 + (D2.error/amp2)^2 + (D3.error/amp3)^2)
    return EHTClosurePhaseDatum(angle(bis), err,
                                u1, v1, u2, v2, u3, v3,
                                time, s123)
end

function baselines(CP::EHTClosurePhaseDatum)
    tri = CP.triangle
    return ((tri[1],tri[2]), (tri[2], tri[3]), (tri[3], tri[1]))
end

uvpositions(datum::EHTClosurePhaseDatum) = (datum.u1, datum.v1, datum.u2, datum.v2, datum.u3, datum.v3)

"""
    $(TYPEDEF)

A Datum for a single log closure amplitude.

# $(FIELDS)

"""
Base.@kwdef struct EHTLogClosureAmplitudeDatum{T<:Number} <: ClosureProducts{T}
    """
    log-closure amplitude
    """
    amp::T
    """
    log-closure amplitude error in the high-snr limit
    """
    error::T
    """
    u (λ) of first station
    """
    u1::T
    """
    v (λ) of first station
    """
    v1::T
    """
    u (λ) of second station
    """
    u2::T
    """
    v (λ) of second station
    """
    v2::T
    """
    u (λ) of third station
    """
    u3::T
    """
    v (λ) of third station
    """
    v3::T
    """
    u (λ) of fourth station
    """
    u4::T
    """
    v (λ) of fourth station
    """
    v4::T
    """
    Observation time of the quadrangle
    """
    time::T
    """
    station codes for the quadrangle
    """
    quadrangle::NTuple{4,Symbol}
end

function rescaleuv!(data::EHTObservation{T,D}, scale) where {T,D<:EHTLogClosureAmplitudeDatum}
    data.data.u1 .= data.data.u1.*scale
    data.data.v1 .= data.data.v1.*scale
    data.data.u2 .= data.data.u2.*scale
    data.data.v2 .= data.data.v2.*scale
    data.data.u3 .= data.data.u3.*scale
    data.data.v3 .= data.data.v3.*scale
    data.data.u4 .= data.data.u4.*scale
    data.data.v4 .= data.data.v4.*scale
    return data
end

function baselines(CP::EHTLogClosureAmplitudeDatum)
    quad = CP.quadrangle
    return ((quad[1],quad[2]), (quad[3], quad[4]), (quad[1], quad[3]), (quad[2], quad[4]))
end

function stations(d::EHTObservation{T,A}) where {T,A<:EHTLogClosureAmplitudeDatum}
    bl = getdata(d, :quadrangle)
    return unique(vcat(collect.(bl)...))
end

uvpositions(datum::EHTLogClosureAmplitudeDatum) = (datum.u1, datum.v1, datum.u2, datum.v2, datum.u3, datum.v3, datum.u4, datum.v4)

# @doc raw"""
#     $(SIGNATURES)
# Compute the log-closure amplitude of 4 visibility datums.
# This uses the same ordering as eht-imaging i.e.

# ```math
# C_A = \log \frac{|V_{1}| |V_{2}|}{|V_{3}| |V_{4}|}
# ```

# # Notes

# Currently we use the high SNR Gaussian approximation to the thermal error.
# In the future we may use the moment matching approach using Monte Carlo sampling.
# """
# function logclosure_amplitude(D1::EHTVisibilityDatum,
#                               D2::EHTVisibilityDatum,
#                               D3::EHTVisibilityDatum,
#                               D4::EHTVisibilityDatum)
#     check_quadrangle(D1,D2,D3,D4)

#     amp1 = amplitude(D1)
#     amp2 = amplitude(D2)
#     amp3 = amplitude(D3)
#     amp4 = amplitude(D4)

#     u1,v1 = uvpositions(D1)
#     u2,v2 = uvpositions(D2)
#     u3,v3 = uvpositions(D3)
#     u4,v4 = uvpositions(D4)

#     #Construct the quadrangle
#     s12 = unique([D1.baseline..., D2.baseline...])
#     s123 = unique([s12..., D3.baseline...])
#     s1234 = unique([s123..., D4.baseline...])

#     # Now do the error propogation
#     err = hypot(D1.err/map1, D2.err/amp2, D3.err/amp3, D4.err/amp4)

#     #construct the lcamp
#     lcamp = log(amp1*amp2/(amp3*amp4))

#     return EHTLogClosureAmplitudeDatum(lcamp, err,
#                                        u1, v1,
#                                        u2, v2,
#                                        u3, v3,
#                                        u4, v4,
#                                        time, s1234)
# end

"""
    $(SIGNATURES)
Extract the array configuration from a visibility EHT observation.
"""
function arrayconfig(vis::EHTObservation)
    vis.config
end

function _arrayconfig(data, bandwidth, frequency)
    u = getproperty(data, :u)
    v = getproperty(data, :v)
    times = getproperty(data, :time)
    error = getproperty(data, :error)
    baseline = getproperty(data, :baseline)
    uvsamples = StructArray{ArrayBaselineDatum}(time=times,
                                        u=u,
                                        v=v,
                                        baseline=baseline,
                                        error_real=error,
                                        error_imag=error
                                    )
    return EHTArrayConfiguration(frequency, bandwidth, uvsamples)
end

const VisAmpDatum = Union{EHTVisibilityAmplitudeDatum, EHTVisibilityDatum}

"""
    `$(TYPEDEF)`
Wraps EHTObservation in a table that separates the observation into scans.
This implements the table interface and is used primarly for keeping track of
gains and closure information.
"""
struct ScanTable{O<:Union{Observation,ArrayConfiguration}, T, S}
    obs::O
    times::T
    scanind::S
end

Base.length(st::ScanTable) = length(st.times)
stations(st::ScanTable) = stations(st.obs)

struct Scan{T,I,S}
    time::T
    index::I
    scan::S
end

Base.length(s::Scan) = length(s.scan)

function baselines(scancp::Scan{A,B,C}) where {A,B,C<:StructArray{<:EHTClosurePhaseDatum}}
    tri = scancp.scan.triangle
    # organize the closure phase stations
    ant1 = getindex.(tri, 1)
    ant2 = getindex.(tri, 2)
    ant3 = getindex.(tri, 3)
    return ant1, ant2, ant3
end

function baselines(scancp::Scan{A,B,C}) where {A,B,C<:StructArray{<:EHTLogClosureAmplitudeDatum}}
    tri = scancp.scan.quadrangle
    # organize the closure phase stations
    ant1 = getindex.(tri, 1)
    ant2 = getindex.(tri, 2)
    ant3 = getindex.(tri, 3)
    ant4 = getindex.(tri, 4)
    return ant1, ant2, ant3, ant4
end

function baselines(scancp::Scan{A,B,C}) where {A,B,C<:StructArray{<:EHTVisibilityDatum}}
    bl = scancp.scan.baseline
    # organize the closure phase stations
    ant1 = first.(bl)
    ant2 = last.(bl)
    return ant1, ant2
end

function baselines(scancp::Scan{A,B,C}) where {A,B,C<:StructArray{<:EHTVisibilityAmplitudeDatum}}
    bl = scancp.scan.baseline
    # organize the closure phase stations
    ant1 = first.(bl)
    ant2 = last.(bl)
    return ant1, ant2
end

function stations(s::Scan)
    ants = baselines(s)
    stat = unique(vcat(ants...))
    return stat
end

function Base.show(io::IO, s::Scan)
    println(io, "VLBI Scan")
    println(io, "\tscan index: ", s.index)
    println(io, "\tscan time:  ", s.time)
    println(io, "\tstations: ", stations(s))
end

function Base.getindex(st::ScanTable, i::Int)
    istart = st.scanind[i]
    if i < length(st.scanind)
        iend = st.scanind[i+1]-1
    else
        iend = length(st.obs)
    end
    return Scan(st.times[i], (i, istart, iend), @view st.obs.data[istart:iend])
end

function Base.getindex(st::ScanTable, I)
    [getindex(st, i) for i in I]
end

function Base.getindex(scan::Scan, s::Symbol)
    getproperty(scan.scan, s)
end

function Base.getindex(scan::Scan, i::Int)
    scan.scan[i]
end

function Base.getindex(scan::Scan, i::AbstractVector{<:Union{Bool,Int}})
    Scan(scan.time, scan.index, scan.scan[i])
end

Base.first(st::ScanTable) = st[1]

"""
    `scantable(obs)`
Reorganizes the observation into a table of scans, where scan are defined by unique timestamps.
To access the data you can use scalar indexing

# Example

```julia
st = scantable(obs)
# Grab the first scan
scan1 = st[1]

# Acess the detections in the scan
scan1[1]

# grab e.g. the baselines
scan1[:baseline]
```

"""
function scantable(obs)
    times = obs[:time]
    scantimes = unique(times)
    scanind = Int[]
    for t in scantimes
        ind = findfirst(==(t), times)
        append!(scanind, ind)
    end
    return ScanTable(obs, scantimes, scanind)
end

include(joinpath(@__DIR__, "io.jl"))
include(joinpath(@__DIR__, "ehtim.jl"))


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 = "TVM"
version = v"0.5"

sources = [
    "https://github.com/dmlc/tvm.git" =>
    "f08015e7fde92c835907d4c9b7ad6d3f634e94a5",
    "https://github.com/dmlc/dlpack.git" =>
    "0acb731e0e43d15deee27b66f10e4c5b4e667913",
    "https://github.com/dmlc/dmlc-core.git" =>
    "808f485387f9a03f78fa9f1159f387d0d91b7a28",
    "https://github.com/agauniyal/rang.git" =>
    "cabe04d6d6b05356fa8f9741704924788f0dd762",
    "https://github.com/dmlc/HalideIR.git" =>
    "b257a9221ee1e5180d994b3488ddcc259b0ac157",
]

# Bash recipe for building across all platforms
script = raw"""
# Start by moving things into the proper places
for submodule in dlpack dmlc-core rang HalideIR; do
    rm -rf "${WORKSPACE}/srcdir/tvm/3rdparty/${submodule}"
    mv "${WORKSPACE}/srcdir/${submodule}" "${WORKSPACE}/srcdir/tvm/3rdparty/${submodule}"
done

# Next, enter TVM and start building it
cd $WORKSPACE/srcdir/tvm

mkdir build && cd build
cmake .. -DCMAKE_INSTALL_PREFIX="${prefix}" \
         -DCMAKE_TOOLCHAIN_FILE="${CMAKE_TARGET_TOOLCHAIN}" \
         -DCMAKE_
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64; libc=:glibc),

    # Fails with can't find `execinfo.h`
    #Linux(:x86_64; libc=:musl),
    MacOS(:x86_64),

    # dll import woes
    #Windows(:x86_64),
]

# The products that we will ensure are always built
products = Product[
    LibraryProduct("libtvm", :libtvm),
    LibraryProduct("libtvm_topi", :libtvm_topi),
    LibraryProduct("libtvm_runtime", :libtvm_runtime),
    LibraryProduct("libnnvm_compiler", :libnnvm_compiler),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "Zlib_jll",
    "LLVM_jll",
]

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


text:
#
# BaseFields are just wrappers around arrays tagged with additional
# information in the form of:
#  * the `B` parameter, for information which is needed at
#    compile-time (right now, just the basis)
#  * the `metadata` field, for information which is only needed at
#    run-time
#

abstract type Proj end
broadcastable(proj::Proj) = Ref(proj)

struct BaseField{B, M<:Proj, T, A<:AbstractArray{T}} <: Field{B, T}
    arr :: A
    metadata :: M
    function (::Type{F})(arr::A, metadata::M) where {B,M<:Proj,T,A<:AbstractArray{T},F<:BaseField{B}}
        check_field_consistency(B(), arr, metadata)
        new{B,M,T,A}(arr, metadata) :: F
    end
end

# can be overrriden for specific types to check things like the
# dimensions of arr being consistent with metadata, etc...
check_field_consistency(::Any, ::Any, ::Any) = ()

typealias_def(::Type{F}) where {B,M,T,A,F<:BaseField{B,M,T,A}} = "BaseField{$(typealias(B)),$(typealias(A)),$(typealias(M))}"

## array interface
# even though the underlying data in BaseFields is multi-dimensional,
# they act as a 1D vector with all their entries splayed out in a
# column
size(f::BaseField) = (length(f.arr),)
lastindex(f::BaseField, i::Int) = lastindex(f.arr, i)
@propagate_inbounds getindex(f::BaseField, I::Union{Int,Colon,AbstractArray}...) = getindex(f.arr, I...)
@propagate_inbounds setindex!(f::BaseField, X, I::Union{Int,Colon,AbstractArray}...) = (setindex!(f.arr, X, I...); f)
similar(f::BaseField{B}, ::Type{T}) where {B,T} = BaseField{B}(similar(f.arr, T), f.metadata)
copy(f::BaseField{B}) where {B} = BaseField{B}(copy(f.arr), f.metadata)

## promotion
function promote(f₁::BaseField{B₁}, f₂::BaseField{B₂}) where {B₁,B₂}
    get_metadata_strict(f₁.metadata, f₂.metadata) # for now just check compatible
    B = typeof(promote_basis_generic(B₁(), B₂()))
    B(f₁), B(f₂)
end

## broadcasting

# any broadcast expression that contains at least one BaseField will
# have a broadcast style of BaseFieldStyle{S,B}. the S is the
# broadcast style for the underlying arrays and B is the B parameter
# of the result BaseField. S and B are promoted recursively across all
# the arguments according to the set of rules below.
struct BaseFieldStyle{S,B} <: AbstractArrayStyle{1} end
BroadcastStyle(::Type{F}) where {B,M,T,A,F<:BaseField{B,M,T,A}} =
    BaseFieldStyle{typeof(BroadcastStyle(A)),B}()
BroadcastStyle(::BaseFieldStyle{S₁,B₁}, ::BaseFieldStyle{S₂,B₂}) where {S₁,B₁,S₂,B₂} =
    BaseFieldStyle{typeof(result_style(S₁(), S₂())), typeof(promote_basis_strict(B₁(),B₂()))}()
BroadcastStyle(S::BaseFieldStyle, ::DefaultArrayStyle{0}) = S

# with the Broadcasted object created, we now compute the answer
function materialize(bc::Broadcasted{BaseFieldStyle{S,B}}) where {S,B}

    # first, recursively go through the broadcast arguments and figure
    # out the metadata of the result, using the
    # promote_metadata_strict rules
    metadata = get_metadata_strict(bc)

    # then "preprocess" all the arguments. this unwraps all of the
    # BaseFields in the expression into just the underlying arrays,
    # and turns things which were ImplicitFields into actual arrays
    # (which are computed by dispatching on the now-known S, B, and
    # metadata of the result)
    bc′ = preprocess((BaseFieldStyle{S,B}(), metadata), bc)

    # the arguments of bc′ are now all normal arrays, so convert it to
    # the broadcast style S that goes along with them
    bc″ = convert(Broadcasted{S}, bc′)

    # run the normal array broadcast, and wrap in the right
    # result type
    BaseField{B}(materialize(bc″), metadata)

end

function materialize!(dst::BaseField{B}, bc::Broadcasted{BaseFieldStyle{S,B′}}) where {B,B′,S}

    (B == B′) || error("Can't broadcast a $(typealias(B′)) into a $(typealias(B))")

    # for inplace broadcasting, we don't need to compute B or the
    # metadata from the broadcasted object, we just take it from the
    # destination BaseField. otherwise its the same as materialize above
    bc′ = preprocess((BaseFieldStyle{S,B}(), dst.metadata), bc)
    bc″ = convert(Broadcasted{S}, bc′)
    materialize!(dst.arr, bc″)
    dst

end

# the default preprocessing, which just unwraps the underlying array.
# this doesn't dispatch on the first argument, but custom BaseFields
# are free to override this and dispatch on it if they need
preprocess(::Any, f::BaseField) = f.arr

# we re-wrap each Broadcasted object as we go through preprocessing
# because some array types do special things here (e.g. CUDA wraps
# bc.f in a CUDA.cufunc)
preprocess(dest::Tuple{BaseFieldStyle{S,B},M}, bc::Broadcasted) where {S,B,M} =
    broadcasted(S(), bc.f, preprocess_args(dest, bc.args)...)

# recursively go through a Broadcasted object's arguments and compute
# the final metadata according to the promote_metadata_strict rules.
# we use the select_known_rule machinery (see util.jl) to make it so
# promote_metadata_strict_rule only needs one argument order defined
# (similar to Base.promote_rule)
get_metadata_strict(x, rest...)      = promote_metadata_strict(get_metadata_strict(x), get_metadata_strict(rest...))
get_metadata_strict(bc::Broadcasted) = get_metadata_strict(bc.args...)
get_metadata_strict(f ::BaseField)   = f.metadata
get_metadata_strict(  ::Any)         = nothing
get_metadata_strict()                = nothing

promote_metadata_strict(x) = x
promote_metadata_strict(x, y) = select_known_rule(promote_metadata_strict_rule, x, y)
promote_metadata_strict_rule(metadata,   ::Nothing) = metadata
promote_metadata_strict_rule(::Nothing,  ::Nothing) = nothing
promote_metadata_strict_rule(::Any,      ::Any) = Unknown()

## properties
getproperty(f::BaseField, s::Symbol)           = getproperty(f,Val(s))
getproperty(f::BaseField,  ::Val{:arr})        = getfield(f,:arr)
getproperty(f::BaseField,  ::Val{:metadata})   = getfield(f,:metadata)
getproperty(f::BaseField,  ::Val{s}) where {s} = getfield(getfield(f,:metadata),s)
propertynames(f::BaseField) = (fieldnames(typeof(f))..., fieldnames(typeof(f.metadata))...)

## CMBLensing-specific stuff
global_rng_for(::Type{BaseField{B,M,T,A}}) where {B,M,T,A} = global_rng_for(A)
fieldinfo(f::BaseField) = f # for backwards compatibility
get_storage(f::BaseField) = typeof(f.arr)
adapt_structure(to, f::BaseField{B}) where {B} = BaseField{B}(adapt(to, f.arr), adapt(to, f.metadata))
hash(f::BaseField, h::UInt64) = foldr(hash, (typeof(f), cpu(f.arr), f.metadata), init=h)

#
default_proj(::Type{F}) where {F<:BaseField{<:Any,<:Proj}} = Base.unwrap_unionall(F).parameters[2].ub
make_field_aliases("Base", Proj)

# useful for enumerating some cases below and in plotting
_sub_components = [
    (:BaseMap,        (:Ix=>:, :I=>:)),
    (:BaseFourier,    (:Il=>:, :I=>:)),
    (:BaseQUMap,      (:Qx=>1, :Ux=>2, :Q =>1, :U=>2, :P=>:)),
    (:BaseQUFourier,  (:Ql=>1, :Ul=>2, :Q =>1, :U=>2, :P=>:)),
    (:BaseEBMap,      (:Ex=>1, :Bx=>2, :E =>1, :B=>2, :P=>:)),
    (:BaseEBFourier,  (:El=>1, :Bl=>2, :E =>1, :B=>2, :P=>:)),
    (:BaseIQUMap,     (:Ix=>1, :Qx=>2, :Ux=>3, :I=>1, :Q=>2, :U=>3, :P=>2:3, :IP=>:)),
    (:BaseIQUFourier, (:Il=>1, :Ql=>2, :Ul=>3, :I=>1, :Q=>2, :U=>3, :P=>2:3, :IP=>:)),
    (:BaseIEBMap,     (:Ix=>1, :Ex=>2, :Bx=>3, :I=>1, :E=>2, :B=>3, :P=>2:3, :IP=>:)),
    (:BaseIEBFourier, (:Il=>1, :El=>2, :Bl=>3, :I=>1, :E=>2, :B=>3, :P=>2:3, :IP=>:)),
]

# sub-components like f.Ix, f.Q, f.P, etc...
for (F, props) in _sub_components
    for (k, I) in props
        if String(k)[end] in "xl"
            if I == (:)
                @eval getproperty(f::$F, ::Val{$(QuoteNode(k))}) = getfield(f,:arr)
            else
                @eval getproperty(f::$F, ::Val{$(QuoteNode(k))}) = view(getfield(f,:arr), pol_slice(f,$I)...)
            end
        else
            if I == (:)
                @eval getproperty(f::$F, ::Val{$(QuoteNode(k))}) = f
            else
                B = (k == :P) ? Basis2Prod{basis(@eval($F)).parameters[end-1:end]...} : basis(@eval($F)).parameters[end]
                @eval getproperty(f::$F, ::Val{$(QuoteNode(k))}) = BaseField{$B}(_reshape_batch(view(getfield(f,:arr), pol_slice(f,$I)...)), f.metadata)
            end
        end
    end
end


text:
using MPI
using Test

MPI.Init()
MPI.finalize_atexit()

excluded = ["runtests.jl"]

@testset "Trilinos tests" begin
  @testset "$f" for f in filter(fname -> fname ∉ excluded, readdir())
    println("======== Starting tests from $f ========")
    include(f)
  end
  sleep(1)
end


text:
module StokesParticles

using figtree_jll

include("struct_defs.jl")
include("inside_checking_functions.jl")
include("cell_lists.jl")
include("fgt.jl")
include("forces.jl")
include("update_particles.jl")
include("params.jl")
include("io.jl")

export generate_cell_list,
       update_cell_list!

export Point2D,
       Particle2D,
       LineWall,
       CircleWall,
       ArcWall

export update_particles_nofluid!,
       update_particles_nofluid_nocl!

export scratch_data,
       sp_params

export write_circ, read_circ,
       write_circ2, read_circ2,
       write_wall, read_wall,
       combine_circ

end # module


text:
function compute_effective_coefficents!(ν_eff, T, mol::MolecularConstants{FT}) where {FT}#molecules::Array{MolecularConstants{FT}}) where {FT}
    #@unpack Y = mol
    @unpack α̅,  γ̅, α_prime, γ_prime, ϵ, ϵ_prime = mol.effCoeff
    @unpack α̅₀₀, γ̅₀₀, ω₀, α_b, α_c, α₀₀_prime, γ₀₀_prime = mol.PolTensor
    #Computing α̅
    α̅ = α̅₀₀*(1 + α_b*T + α_c*T^2)/(1-(2π*c*ν_eff/ω₀)^2)
    γ̅ = γ̅₀₀
    ϵ = α̅/γ̅
    α_prime = α₀₀_prime * sqrt(mol.Y[1,2]/mol.Y[2,1]) # see Eqs (36a-39b) of Buldakov et al. 1996 (Spectrochimica Acta Part A)
    γ_prime = γ₀₀_prime * sqrt(mol.Y[1,2]/mol.Y[2,1]) # see Eqs (36a-39b) of Buldakov et al. 1996 (Spectrochimica Acta Part A)
    γ_C_Rayl = 3/(45(ϵ)^2+4)
    γ_C_RotRaman = 3/4
    ϵ_prime = α_prime/γ_prime
    γ_C_VibRaman =  3/(45(ϵ_prime)^2+4)
    γ_C_RoVibRaman = 3/4
    rho_depol_Rayl = 2γ_C_Rayl/(1 + γ_C_Rayl) # depolarization factor, as defined in Table 1 (Spurr 2006)
    rho_depol_RotRaman = 2γ_C_RotRaman/(1 + γ_C_RotRaman) # depolarization factor, as defined in Table 1 (Spurr 2006)
    rho_depol_VibRaman = 2γ_C_VibRaman/(1 + γ_C_VibRaman) # depolarization factor, as defined in Table 1 (Spurr 2006)
    rho_depol_RoVibRaman = 2γ_C_RoVibRaman/(1 + γ_C_RoVibRaman) # depolarization factor, as defined in Table 1 (Spurr 2006)

    @pack! mol.effCoeff = rho_depol_Rayl, rho_depol_RotRaman, rho_depol_VibRaman, rho_depol_RoVibRaman,
        α̅, γ̅, α_prime, γ_prime, ϵ, ϵ_prime, γ_C_Rayl, γ_C_RotRaman, γ_C_VibRaman, γ_C_RoVibRaman
end

#Compute elastic scattering cross-section (Cabannes line)
function compute_σ_Rayl_coeff!(mol::MolecularConstants{FT}) where {FT}#ν, molecules::Array{MolecularConstants{FT}}) where {FT}
    @unpack α̅, γ_C_Rayl, σ_Rayl_coeff = mol.effCoeff
    σ_Rayl_coeff = 128π^5 * α̅^2 * (1+2*γ_C_Rayl)/(3-4*γ_C_Rayl)# * ν^4
    @show σ_Rayl_coeff,  α̅^2
    @pack! mol.effCoeff = σ_Rayl_coeff,α̅,γ_C_Rayl
end

#Compute elastic scattering cross-section (Cabannes line)
function compute_σ_Rayl_VibRaman_coeff_hires!(T, mol::MolecularConstants{FT}; Jmax=30) where {FT}#ν, molecules::Array{MolecularConstants{FT}}) where {FT}
    @unpack α̅, γ̅, α_prime, γ_prime, E_vJ, σ_Rayl_coeff_hires, σ_VibRaman_coeff_0to1_hires, σ_VibRaman_coeff_1to0_hires, Δν̃_Rayl_coeff_hires, Δν̃_VibRaman_coeff_0to1_hires, Δν̃_VibRaman_coeff_1to0_hires = mol.effCoeff
    @unpack gₛ = mol
    σ_Rayl_coeff_hires = OffsetArray(zeros(FT, Jmax+1), 0:Jmax);
    Δν̃_Rayl_coeff_hires = OffsetArray(zeros(FT, Jmax+1), 0:Jmax);
    σ_VibRaman_coeff_0to1_hires = OffsetArray(zeros(FT, Jmax+1), 0:Jmax);
    σ_VibRaman_coeff_1to0_hires = OffsetArray(zeros(FT, Jmax+1), 0:Jmax);
    Δν̃_VibRaman_coeff_0to1_hires = OffsetArray(zeros(FT, Jmax+1), 0:Jmax);
    Δν̃_VibRaman_coeff_1to0_hires = OffsetArray(zeros(FT, Jmax+1), 0:Jmax);

    kᵥ = 128*π^5
    #σ_Rayl_coeff_hires = 128π^5 * α̅^2 * (1+2*γ_C_Rayl)/(3-4*γ_C_Rayl)# * ν^4
    Z_pf = 0
    for Ji in 0:Jmax
        b_JJ   = Ji*(Ji+1)/((2Ji-1)*(2Ji+3))
        g_N = gₛ[2-Ji%2]

        # Rayleigh
        vi = 0
        vf = 0
        Ni_by_N = exp(-h*c*E_vJ[vi,Ji]/(k_B*T))
        γ_C = 3/(4+45(α̅/(b_JJ*γ̅))^2)
        Z_pf += g_N * (2Ji+1) * (exp(-h*c*E_vJ[0,Ji]/(k_B*T)) + exp(-h*c*E_vJ[1,Ji]/(k_B*T))) #accounting for both the ground- and the first excited vibrational state - depending on the ambient temperature (e.g. on hot exoplanets), more states may have to be accounted for
        # Molecular energy level
        Δν =  E_vJ[vi, Ji]
        # Wavenumber change of scattered light
        Δνₛ = Δν

        σ_Rayl_coeff_hires[Ji] = kᵥ * g_N * (2Ji+1) * Ni_by_N * α̅^2 * (1+2*γ_C)/(3-4*γ_C)# * νₛ[vi, Ji, vf, Jf]^4
        Δν̃_Rayl_coeff_hires[Ji] = Δνₛ

        # Vib. Raman 0->1
        vi = 0
        vf = 1
        Ni_by_N = exp(-h*c*E_vJ[vi,Ji]/(k_B*T))
        γ_C = 3/(4+45(α_prime/(b_JJ*γ_prime))^2)
        # Molecular energy level
        Δν =  E_vJ[vf, Ji] - E_vJ[vi, Ji]
        # Wavenumber change of scattered light
        Δνₛ = -Δν

        σ_VibRaman_coeff_0to1_hires[Ji] = kᵥ * g_N * (2Ji+1) * Ni_by_N * α_prime^2 * (1+2*γ_C)/(3-4*γ_C)# * νₛ[vi, Ji, vf, Jf]^4
        Δν̃_VibRaman_coeff_0to1_hires[Ji] = Δνₛ
        @show Ji, α_prime^2, (1+2*γ_C)/(3-4*γ_C), α_prime^2 * (1+2*γ_C)/(3-4*γ_C)
        # Vib. Raman 1->0
        vi = 1
        vf = 0
        Ni_by_N = exp(-h*c*E_vJ[vi,Ji]/(k_B*T))

        # Molecular energy level
        Δν =  E_vJ[vf, Ji] - E_vJ[vi, Ji]
        # Wavenumber change of scattered light
        Δνₛ = -Δν

        σ_VibRaman_coeff_1to0_hires[Ji] = kᵥ * g_N * (2Ji+1) * Ni_by_N * α_prime^2 * (1+2*γ_C)/(3-4*γ_C)# * νₛ[vi, Ji, vf, Jf]^4
        Δν̃_VibRaman_coeff_1to0_hires[Ji] = Δνₛ
    end
    σ_Rayl_coeff_hires /= Z_pf
    σ_VibRaman_coeff_0to1_hires /= Z_pf
    σ_VibRaman_coeff_1to0_hires /= Z_pf
    #@show Z_pf
    @pack! mol.effCoeff = σ_Rayl_coeff_hires, σ_VibRaman_coeff_0to1_hires, σ_VibRaman_coeff_1to0_hires, Δν̃_Rayl_coeff_hires, Δν̃_VibRaman_coeff_0to1_hires, Δν̃_VibRaman_coeff_1to0_hires
end

#Compute energy levels [in wavenumbers [cm^{-1}]] corresponding to v={0, 1, 2} and J={0, 1, 2,..., 10}
function compute_energy_levels!(mol::MolecularConstants{FT}; vmax=2, Jmax=30) where {FT}#molecules::Array{MolecularConstants{FT}}) where {FT}
    #for mol in molecules
        @unpack E_vJ = mol.effCoeff
        @unpack Y = mol
        E_vJ = OffsetArray(zeros(FT, vmax+1,Jmax+1), 0:vmax, 0:Jmax);
        for v in 0:vmax, J in 0:Jmax
            E₁ = J*(J+1)
            for k in 1:5
                E₂ = (v+0.5)^(k-1)
                for l in 1:5
                    E₃ = E₁^(l-1) * E₂ * Y[k,l]
                    E_vJ[v,J] += E₃
                end
            end
        end
        @pack! mol.effCoeff = E_vJ
    #end
end

#Compute vibrational Raman scattering coefficient (for Δν=±1)
function compute_σ_VibRaman_coeff!(T, mol::MolecularConstants{FT}; vmax=2, Jmax=30) where {FT}#ν, molecules::Array{MolecularConstants{FT}}) where {FT}
    @unpack α_prime, γ_C_VibRaman, E_vJ,σ_VibRaman_coeff_0to1,σ_VibRaman_coeff_1to0, Δν̃_VibRaman_coeff_0to1, Δν̃_VibRaman_coeff_1to0 = mol.effCoeff

    #σ_Rayl_coeff = (128/3)π^5 * α̅^2 * F_King# * ν^4
    #Stokes
    vi = 0
    Δν̃ = E_vJ[1,0]-E_vJ[0,0]
    Nvib = (1 - exp(-h*c*Δν̃/(k_B*T)))^(-1)
    #sq_bₖ = h/(8π²c*E_νJ[vi,0])#bₖ = √(h/8π²cνₖ)
    σ_VibRaman_coeff_0to1 = 128π^5 * α_prime^2 * Nvib * (1+2*γ_C_VibRaman)/(3-4*γ_C_VibRaman)
    Δν̃_VibRaman_coeff_0to1 = -Δν̃
    #Anti-Stokes
    vi = 1

    Nvib = (exp(h*c*Δν̃/(k_B*T))-1)^(-1)
    #sq_bₖ = h/(8π²c*E_νJ[vi,0])#bₖ = √(h/8π²cνₖ)
    σ_VibRaman_coeff_1to0 = 128π^5 * α_prime^2 * Nvib * (1+2*γ_C_VibRaman)/(3-4*γ_C_VibRaman)
    Δν̃_VibRaman_coeff_1to0 = Δν̃

    @show σ_VibRaman_coeff_0to1,σ_VibRaman_coeff_1to0
    @pack! mol.effCoeff = σ_VibRaman_coeff_0to1, σ_VibRaman_coeff_1to0, Δν̃_VibRaman_coeff_0to1, Δν̃_VibRaman_coeff_1to0
end

function compute_σ_RoVibRaman_coeff!(T, mol::MolecularConstants{FT}; vmax=2, Jmax=30) where {FT}#molecules::Array{MolecularConstants{FT}}) where {FT}
    kᵥ = (256/27)*π^5
    #for mol in molecules
    @unpack γ̅, γ_prime,
        E_vJ, σ_RoRaman_coeff_JtoJm2,σ_RoRaman_coeff_JtoJp2,
        σ_RoVibRaman_coeff_0to1_JtoJm2,σ_RoVibRaman_coeff_0to1_JtoJp2,
        σ_RoVibRaman_coeff_1to0_JtoJm2,σ_RoVibRaman_coeff_1to0_JtoJp2,
        Δν̃_RoRaman_coeff_JtoJm2, Δν̃_RoRaman_coeff_JtoJp2,
        Δν̃_RoVibRaman_coeff_0to1_JtoJm2, Δν̃_RoVibRaman_coeff_0to1_JtoJp2,
        Δν̃_RoVibRaman_coeff_1to0_JtoJm2, Δν̃_RoVibRaman_coeff_1to0_JtoJp2 = mol.effCoeff
    @unpack gₛ = mol

    #σ_VibRaman_coeff = OffsetArray(zeros(FT, vmax+1,vmax+1), 0:vmax, 0:vmax);
    σ_RoVibRaman_coeff_0to1_JtoJm2 = OffsetArray(zeros(FT, Jmax+1), 0:Jmax); #zeros(3, 11, 3, 11) #(vi, ji, vf, jf)
    σ_RoVibRaman_coeff_0to1_JtoJp2 = OffsetArray(zeros(FT, Jmax+1), 0:Jmax);
    σ_RoVibRaman_coeff_1to0_JtoJm2 = OffsetArray(zeros(FT, Jmax+1), 0:Jmax); #zeros(3, 11, 3, 11) #(vi, ji, vf, jf)
    σ_RoVibRaman_coeff_1to0_JtoJp2 = OffsetArray(zeros(FT, Jmax+1), 0:Jmax);
    σ_RoRaman_coeff_JtoJm2 = OffsetArray(zeros(FT, Jmax+1), 0:Jmax);
    σ_RoRaman_coeff_JtoJp2 = OffsetArray(zeros(FT, Jmax+1), 0:Jmax);

    Δν̃_RoVibRaman_coeff_0to1_JtoJm2 = OffsetArray(zeros(FT, Jmax+1), 0:Jmax); #zeros(3, 11, 3, 11) #(vi, ji, vf, jf)
    Δν̃_RoVibRaman_coeff_0to1_JtoJp2 = OffsetArray(zeros(FT, Jmax+1), 0:Jmax);
    Δν̃_RoVibRaman_coeff_1to0_JtoJm2 = OffsetArray(zeros(FT, Jmax+1), 0:Jmax); #zeros(3, 11, 3, 11) #(vi, ji, vf, jf)
    Δν̃_RoVibRaman_coeff_1to0_JtoJp2 = OffsetArray(zeros(FT, Jmax+1), 0:Jmax);
    Δν̃_RoRaman_coeff_JtoJm2 = OffsetArray(zeros(FT, Jmax+1), 0:Jmax);
    Δν̃_RoRaman_coeff_JtoJp2 = OffsetArray(zeros(FT, Jmax+1), 0:Jmax);
    #Δνₛ = OffsetArray(zeros(FT, vmax+1,Jmax+1,vmax+1,Jmax+1), 0:vmax, 0:Jmax, 0:vmax, 0:Jmax); #zeros(3, 11, 3, 11)

    Z_pf = 0
    for Ji in 0:Jmax
        b_JJm2 = 3Ji*(Ji-1)/(2*(2Ji+1)*(2Ji-1))
        #b_JJ   = Ji*(Ji+1)/((2Ji-1)*(2Ji+3))
        b_JJp2 = 3*(Ji+1)*(Ji+2)/(2*(2Ji+1)*(2Ji+3))
        g_N = gₛ[2-Ji%2]

        #Rotational Raman scattering
        vi = 0
        vf = 0
        Ni_by_N = exp(-h*c*E_vJ[vi,Ji]/(k_B*T))
        @show vi, Ji, Ni_by_N
        Z_pf += g_N * (2Ji+1) * (exp(-h*c*E_vJ[0,Ji]/(k_B*T)) + exp(-h*c*E_vJ[1,Ji]/(k_B*T))) #accounting for both the ground- and the first excited vibrational state - depending on the ambient temperature (e.g. on hot exoplanets), more states may have to be accounted for

        Jf=Ji-2
        if (Ji-2)>=0 #anti-Stokes
            #Change in molecular energy level
            Δν = E_vJ[vf, Jf] - E_vJ[vi, Ji]
            #Wavenumber change of scattered light
            Δνₛ = - Δν

            σ_RoRaman_coeff_JtoJm2[Ji] = kᵥ * g_N * (2Ji+1) * b_JJm2 * Ni_by_N * γ̅^2# * νₛ[vi, Ji, vf, Jf]^4
            Δν̃_RoRaman_coeff_JtoJm2[Ji] = Δνₛ
        else
            σ_RoRaman_coeff_JtoJm2[Ji] = 0
            Δν̃_RoRaman_coeff_JtoJm2[Ji] = 0
        end
        Jf=Ji+2
        if (Ji+2)<=Jmax #Stokes
            #Change in molecular energy level
            Δν = E_vJ[vf, Jf] - E_vJ[vi, Ji]
            #Wavenumber change of scattered light
            Δνₛ = - Δν

            σ_RoRaman_coeff_JtoJp2[Ji] = kᵥ * g_N * (2Ji+1) * b_JJp2 * Ni_by_N * γ̅^2# * νₛ[vi, Ji, vf, Jf]^4
            Δν̃_RoRaman_coeff_JtoJp2[Ji] = Δνₛ
        else
            σ_RoRaman_coeff_JtoJp2[Ji] = 0
            Δν̃_RoRaman_coeff_JtoJp2[Ji] = 0
        end

        #Rovibrational Raman scattering: Stokes
        vi = 0
        vf = 1

        Ni_by_N = exp(-h*c*E_vJ[vi,Ji]/(k_B*T))

        Jf=Ji-2
        #b_v = (1/2π)*(h/(2c*E_vJ[vf,Jf]))^0.5
        if (Ji-2)>=0 #anti-Stokes
            #Change in molecular energy level
            Δν = E_vJ[vf, Jf] - E_vJ[vi, Ji]
            #Wavenumber change of scattered light
            Δνₛ = - Δν

            σ_RoVibRaman_coeff_0to1_JtoJm2[Ji] = kᵥ * g_N * (2Ji+1) * b_JJm2 * Ni_by_N * (γ_prime)^2# * νₛ[vi, Ji, vf, Jf]^4
            Δν̃_RoVibRaman_coeff_0to1_JtoJm2[Ji] = Δνₛ
        else
            σ_RoVibRaman_coeff_0to1_JtoJm2[Ji] = 0
            Δν̃_RoVibRaman_coeff_0to1_JtoJm2[Ji] = 0
        end
        Jf=Ji+2
        #b_v = (1/2π)*(h/(2c*E_vJ[vf,Jf]))^0.5
        if (Ji+2)<=Jmax #Stokes
            #Change in molecular energy level
            Δν = E_vJ[vf, Jf] - E_vJ[vi, Ji]
            #Wavenumber change of scattered light
            Δνₛ = - Δν

            σ_RoVibRaman_coeff_0to1_JtoJp2[Ji] = kᵥ * g_N * (2Ji+1) * b_JJp2 * Ni_by_N * (γ_prime)^2# * νₛ[vi, Ji, vf, Jf]^4
            Δν̃_RoVibRaman_coeff_0to1_JtoJp2[Ji] = Δνₛ
        else
            σ_RoVibRaman_coeff_0to1_JtoJp2[Ji] = 0
            Δν̃_RoVibRaman_coeff_0to1_JtoJp2[Ji] = 0
        end

        #Rovibrational Raman scattering: anti-Stokes
        vi = 1
        vf = 0

        Ni_by_N = exp(-h*c*E_vJ[vi,Ji]/(k_B*T))

        Jf=Ji-2
        if (Ji-2)>=0 #anti-Stokes
            #Change in molecular energy level
            Δν = E_vJ[vf, Jf] - E_vJ[vi, Ji]
            #Wavenumber change of scattered light
            Δνₛ = - Δν

            σ_RoVibRaman_coeff_1to0_JtoJm2[Ji] = kᵥ * g_N * (2Ji+1) * b_JJm2 * Ni_by_N * (γ_prime)^2# * νₛ[vi, Ji, vf, Jf]^4
            Δν̃_RoVibRaman_coeff_1to0_JtoJm2[Ji] = Δνₛ
        else
            σ_RoVibRaman_coeff_1to0_JtoJm2[Ji] = 0
            Δν̃_RoVibRaman_coeff_1to0_JtoJm2[Ji] = 0
        end
        Jf=Ji+2
        if (Ji+2)<=Jmax #Stokes
            #Change in molecular energy level
            Δν = E_vJ[vf, Jf] - E_vJ[vi, Ji]
            #Wavenumber change of scattered light
            Δνₛ = - Δν

            σ_RoVibRaman_coeff_1to0_JtoJp2[Ji] = kᵥ * g_N * (2Ji+1) * b_JJp2 * Ni_by_N * (γ_prime)^2# * νₛ[vi, Ji, vf, Jf]^4
            Δν̃_RoVibRaman_coeff_1to0_JtoJp2[Ji] = Δνₛ
        else
            σ_RoVibRaman_coeff_1to0_JtoJp2[Ji] = 0
            Δν̃_RoVibRaman_coeff_1to0_JtoJp2[Ji] = 0
        end
    end
    σ_RoRaman_coeff_JtoJm2 /= Z_pf
    σ_RoRaman_coeff_JtoJp2 /= Z_pf
    σ_RoVibRaman_coeff_0to1_JtoJm2 /= Z_pf
    σ_RoVibRaman_coeff_0to1_JtoJp2 /= Z_pf
    σ_RoVibRaman_coeff_1to0_JtoJm2 /= Z_pf
    σ_RoVibRaman_coeff_1to0_JtoJp2 /= Z_pf

    @show Z_pf
    @pack! mol.effCoeff = σ_RoRaman_coeff_JtoJm2, σ_RoRaman_coeff_JtoJp2,
        σ_RoVibRaman_coeff_0to1_JtoJm2, σ_RoVibRaman_coeff_0to1_JtoJp2,
        σ_RoVibRaman_coeff_1to0_JtoJm2, σ_RoVibRaman_coeff_1to0_JtoJp2,
        Δν̃_RoRaman_coeff_JtoJm2, Δν̃_RoRaman_coeff_JtoJp2,
        Δν̃_RoVibRaman_coeff_0to1_JtoJm2, Δν̃_RoVibRaman_coeff_0to1_JtoJp2,
        Δν̃_RoVibRaman_coeff_1to0_JtoJm2, Δν̃_RoVibRaman_coeff_1to0_JtoJp2
end

max_stars_repo_name: RadiativeTransfer/InelasticScattering

max_issues_repo_name: RadiativeTransfer/InelasticScattering

max_forks_repo_name: RadiativeTransfer/InelasticScattering


text:
# This file was generated, do not modify it. # hide
ridge_model = Ridge(lambda=0.1)
ridge = machine(ridge_model, W, z)
ẑ = predict(ridge, W)

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
require("test.jl")

using TimeSeries

my_tests = ["test/plot.jl"]

print_with_color(:cyan, "Running tests: ")
println("")

for my_test in my_tests
    print_with_color(:magenta, "**   ")
    print_with_color(:blue, "$my_test")
    println("")
    include(my_test)
end


text:
using ObjectTracker
using Base.Test

@time @testset "Test blob creation" begin include("test_blobs.jl") end
@time @testset "Test object tracking" begin include("test_tracker.jl") end


text:
using SIMD

include("sorting-network-parameters.jl")

"""
    sort_net(input::Vararg{T, L}) where {L,T}

Applies a sorting network of size L to the input elements, returning a sorted tuple.

The elements must support the `min` and `max` functions. In the case of `SIMD.Vec` objects each "lane" across the
vectors will be sorted. Therefore with L vectors of size N this function will produce N sorted sequences of size L,
after the data is transposed (see transpose_vecs).
"""

@generated function sort_net(input::Vararg{T, L})::NTuple{T, L} where {L,T}

    ex = Expr[Expr(:meta, :inline)]

    net_params = sorting_network_parameters[L]

    nsteps = length(net_params)

    for t in 1:L
        a1 = Symbol("input_", 0, "_", t)
        push!(ex, :($a1 = input[$t]))
    end

    for st in 1:nsteps
        touched = [x for t in net_params[st] for x in t]
        untouched = setdiff(1:L, touched)

        for t in untouched
            a1 = Symbol("input_",st-1, "_", t)
            b1 = Symbol("input_",st, "_", t)
            push!(ex, :($b1 = $a1))
        end

        for t in net_params[st]
            a1 = Symbol("input_",st-1, "_", t[1])
            a2 = Symbol("input_", st-1, "_", t[2])
            b1 = Symbol("input_", st, "_", t[1])
            b2 = Symbol("input_", st, "_", t[2])
            push!(ex, :($b1 = min($a1, $a2)))
            push!(ex, :($b2 = max($a1, $a2)))
        end
    end

    push!(ex,
          Expr(:tuple, ntuple(t->Symbol("input_", nsteps, "_", t), Val(L))...))

    quote $(ex...) end
end


text:
struct GUIManager{WM<:AbstractWindowManager}
    wm::WM
    callbacks::Dictionary{Widget,WidgetCallbacks}
end

GUIManager(wm::AbstractWindowManager, callbacks=Dictionary{Widget,WidgetCallbacks}()) = GUIManager{typeof(wm)}(wm, widgets, callbacks)

widgets(gm::GUIManager) = keys(gm.callbacks)

WindowAbstractions.run(gm::GUIManager, mode::ExecutionMode; kwargs...) = run(gm.wm, mode, (ed) -> execute_callback(gm, ed); kwargs...)

function _execute_callback(gm::GUIManager, ed::EventDetails{T}) where {T}
    wm = gm.wm
    act = action(ed)()

    # handle drag/drop operations
    if !isempty(wm.active_buttons) && act in (PointerMoves(), ButtonReleased())
        buttons = collect(values(wm.active_buttons))
        src_idx = argmin(map(x -> x.time, buttons))
        src_ed = buttons[src_idx]
        src_w = find_target(gm, src_ed)
        if !isnothing(src_w) && hypot((Point(src_ed.location) - Point(ed.location))...) > 15.
            w = find_target(gm, ed)
            if act == PointerMoves()
                execute_callback(callbacks(gm, src_w).on_drag, (src_w, src_ed, w, ed))
            elseif act == ButtonReleased()
                @debug "Drop operation"
                execute_callback(callbacks(gm, src_w).on_drop, (src_w, src_ed, w, ed))
            end
        end
        return
    end

    w = find_target(gm, ed)
    !isnothing(w) || return
    cbs = callbacks(gm, w)
    if act == PointerMoves()
        execute_callback(cbs.on_pointer_move, (w, ed))
    elseif act == ButtonPressed() && ed.data.button ∉ [ButtonScrollUp(), ButtonScrollDown()]
        @debug "Button pressed ($(ed.data.button))"
        execute_callback(cbs.on_mouse_button_pressed, (w, ed))
        if !isempty(wm.history)
            pressed = ed.data.button
            last = get(wm.history, act, nothing)
            if !isnothing(last) && last.data.button == pressed
                Δt = ed.time - last.time
                if Δt ≤ 0.5
                    @debug "Double click (Δt=$Δt, $pressed)"
                    execute_callback(cbs.on_double_click, (w, ed))
                end
            end
        end
    else
        cb = getproperty(cbs, callback_symbol(typeof(act)))
        execute_callback(cb, (w, ed))
    end
end

function WindowAbstractions.execute_callback(gm::GUIManager, ed::EventDetails{T}) where {T}
    _execute_callback(gm, ed)
    execute_callback(gm.wm, ed)
end