text:
# ------------------------------------------------------------------------------------------
# # Julia for Data Science
# Prepared by [@nassarhuda](https://github.com/nassarhuda)! 😃
#
# `Last updated on 03/Jan/2020`
#
# In this tutorial, we will discuss why *Julia* is the tool you want to use for your data
# science applications.
#
# We will cover the following:
# * **Data**
# * Data processing
# * Visualization
#
# ### Data: Build a strong relationship with your data.
# Every data science task has one main ingredient, the _data_! Most likely, you want to use
# your data to learn something new. But before the _new_ part, what about the data you
# already have? Let's make sure you can **read** it, **store** it, and **understand** it
# before you start using it.
#
# Julia makes this step really easy with data structures and packages to process the data,
# as well as, existing functions that are readily usable on your data.
#
# The goal of this first part is get you acquainted with some Julia's tools to manage your
# data.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# First, let's download a csv file from github that we can work with.
#
# Note: `download` depends on external tools such as curl, wget or fetch. So you must have
# one of these.
# ------------------------------------------------------------------------------------------

P = download("https://raw.githubusercontent.com/nassarhuda/easy_data/master/programming_languages.csv","programminglanguages.csv")

# ------------------------------------------------------------------------------------------
# We can use shell commands like `ls` in Julia by preceding them with a semicolon.
# ------------------------------------------------------------------------------------------

;ls

# ------------------------------------------------------------------------------------------
# And there's the *.csv file we downloaded!
#
# Now let's load it into Julia
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We'll use the `DelimitedFiles` standard library package and its `readdlm()` function as
# shown
# below.
#
# (Today, the [CSV.jl](https://juliadata.github.io/CSV.jl/stable/) package is the
# recommended way to load CSVs in Julia. We can install it via `Pkg.add()`, and load .csv
# files using `CSV.read()`. This tutorial hasn't been updated to use CSV.jl yet.)
# ------------------------------------------------------------------------------------------

# using Pkg
# Pkg.add("CSV") # for CSV.read()
# Pkg.add("DelimitedFiles") # for readdlm

# using CSV
# P = CSV.read("programminglanguages.csv",header=true)
# or
using DelimitedFiles  # Standard library in Base

# ------------------------------------------------------------------------------------------
# By default, `readdlm` will fill an array with the data stored in the input .csv file. If
# we set the keyword argument `header` to `true`, we'll get a second output array for just
# the headers.
# ------------------------------------------------------------------------------------------

P,H = readdlm("programminglanguages.csv", ',', header=true)

P # stores the dataset

H # stores the header names

# ------------------------------------------------------------------------------------------
# Here we write our first small function.

# Now you can answer questions such as, "when was language X created?"
# ------------------------------------------------------------------------------------------

function language_created_year(P,language::String)
    loc = findfirst(P[:,2].==language)
    return P[loc,1]
end

language_created_year(P,"Julia")

language_created_year(P,"julia")

# ------------------------------------------------------------------------------------------
# As expected, this will not return what you want, but thankfully, string manipulation is
# really easy in Julia!
# ------------------------------------------------------------------------------------------

function language_created_year_v2(P,language::String)
    loc = findfirst(lowercase.(P[:,2]).==lowercase.(language))
    return P[loc,1]
end
language_created_year_v2(P,"julia")

# ------------------------------------------------------------------------------------------
# **Reading and writing to files is really easy in Julia.**

#
# You can use different delimiters with the function `readdlm`, from the `DelimitedFiles`
# package.

#
# To write to files, we can use `writedlm`.

#
# Let's write this same data to a file with a different delimiter.
# ------------------------------------------------------------------------------------------

writedlm("programming_languages_data.txt", P, '-')

# ------------------------------------------------------------------------------------------
# We can now check that this worked using a shell command to glance at the file,
# ------------------------------------------------------------------------------------------

;head -10 programming_languages_data.txt

# ------------------------------------------------------------------------------------------
# and also check that we can use `readdlm` to read our new text file correctly.
# ------------------------------------------------------------------------------------------

P_new_delim = readdlm("programming_languages_data.txt", '-');
P == P_new_delim

# ------------------------------------------------------------------------------------------
# ### Dictionaries
# Let's try to store the above data in a dictionary format!
#
# First, let's initialize an empty dictionary
# ------------------------------------------------------------------------------------------

dict = Dict{Integer,Vector{String}}()

# ------------------------------------------------------------------------------------------
# Here we told Julia that we want `dict` to only accept integers as keys and vectors of
# strings as values.
#
# However, we could have initialized an empty dictionary without providing this information
# (depending on our application).
# ------------------------------------------------------------------------------------------

dict2 = Dict()

# ------------------------------------------------------------------------------------------
# This dictionary takes keys and values of any type!
#
# Now, let's populate the dictionary with years as keys and vectors that hold all the
# programming languages created in each year as their values.
# ------------------------------------------------------------------------------------------

for i = 1:size(P,1)
    year,lang = P[i,:]

    if year in keys(dict)
        dict[year] = push!(dict[year],lang)
    else
        dict[year] = [lang]
    end
end

# ------------------------------------------------------------------------------------------
# Now you can pick whichever year you want and find what programming languages were invented
# in that year
# ------------------------------------------------------------------------------------------

dict[2003]

# ------------------------------------------------------------------------------------------
# ### DataFrames!
# *Shout out to R fans!*
# One other way to play around with data in Julia is to use a DataFrame.
#
# This requires loading the `DataFrames` package. Thankfully, this tutorial came with a
# Project.toml file that specifies exactly which version of DataFrames to install...
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Project.toml files
#
# For this tutorial (Julia for Data Science), you may have noticed that there are files in
# this folder called [`Project.toml`](/edit/introductory-tutorials/broader-topics-and-
# ecosystem/intro-to-julia-for-data-science/Project.toml) and
# [`Manifest.toml`](/edit/introductory-tutorials/broader-topics-and-ecosystem/intro-to-
# julia-for-data-science/Manifest.toml). These are files autogenerated by Julia's package
# manager, `Pkg`, that describe the _exact set of packages_ installed for a julia project,
# allowing you to share your work in a perfectly reproducible way.
#
# Jupyter was able to detect those `.toml` files, and so this notebook was automatically
# started with _this project activated!_ Note that this means any packages you add or remove
# inside this notebook will only affect this "Julia for Data Science" _project_.
#
# To install all of the package dependencies used in the rest of the tutorial, you only need
# to run this next cell (commands that start with `]` are package repl commands):
# ------------------------------------------------------------------------------------------

] instantiate

# ------------------------------------------------------------------------------------------
# You can read more about package manager commands, here:
# https://docs.julialang.org/en/v1/stdlib/Pkg/index.html
#
# **Now back to DataFrames!**
# ------------------------------------------------------------------------------------------

using DataFrames
df = DataFrame(year = P[:,1], language = P[:,2])

# ------------------------------------------------------------------------------------------
# You can access columns by header name, or column index.
#
# In this case, `df[1]` is equivalent to `df.year` or `df[!, :year]`.
#
# Note that if we want to index columns by header name, we precede the header name with a
# colon. In Julia, this means that the header names are treated as *symbols*.
# ------------------------------------------------------------------------------------------

df.year

# ------------------------------------------------------------------------------------------
# **`DataFrames` provides some handy features when dealing with data**
#
# First, it uses julia's "missing" type.
# ------------------------------------------------------------------------------------------

a = missing
typeof(a)

# ------------------------------------------------------------------------------------------
# Let's see what happens when we try to add a "missing" type to a number.
# ------------------------------------------------------------------------------------------

a + 1

# ------------------------------------------------------------------------------------------
# `DataFrames` provides the `describe` function, which can give you quick statistics about
# each column in your dataframe
# ------------------------------------------------------------------------------------------

describe(df)

# ------------------------------------------------------------------------------------------
# ### RDatasets
#
# We can use RDatasets to play around with pre-existing datasets
# ------------------------------------------------------------------------------------------

using RDatasets
iris = dataset("datasets", "iris")

# ------------------------------------------------------------------------------------------
# Note that data loaded with `dataset` is stored as a DataFrame. 😃
# ------------------------------------------------------------------------------------------

typeof(iris)

# ------------------------------------------------------------------------------------------
# The summary we get from `describe` on `iris` gives us a lot more information than the
# summary on `df`!
# ------------------------------------------------------------------------------------------

describe(iris)

# ------------------------------------------------------------------------------------------
# ### More on Missing Values
#
# Julia 1.0 and beyond has native support for `missing` values. (Before Julia 1.0, this was
# done via the DataArrays.jl package.)
# More information on using arrays with missing values can be found
# [in the Julia documentation](https://docs.julialang.org/en/v1/manual/missing/#Arrays-With-
# Missing-Values-1).
# ------------------------------------------------------------------------------------------

foods = ["apple", "cucumber", "tomato", "banana"]
calories = [missing,47,22,105]
typeof(calories)

using Statistics  # julia's standard library for stats
mean(calories)

# ------------------------------------------------------------------------------------------
# Missing values ruin everything! 😑
#
# Luckily we can ignore them with `skipmissing`!
# ------------------------------------------------------------------------------------------

mean(skipmissing(calories))

# ------------------------------------------------------------------------------------------
# Oh WAIT! Detour. How did I get the emoji there?
#
# Try this out:
#
# ```
# \:expressionless: + <TAB>
# ```
# ------------------------------------------------------------------------------------------

😑 = 0 # expressionless
😀 = 1
😞 = -1

# ------------------------------------------------------------------------------------------
# *Back to missing values*
# ------------------------------------------------------------------------------------------

prices = [0.85,1.6,0.8,0.6,]

dataframe_calories = DataFrame(item=foods,calories=calories)

dataframe_prices = DataFrame(item=foods,price=prices)

# ------------------------------------------------------------------------------------------
# We can also `join` two dataframes together
# ------------------------------------------------------------------------------------------

DF = join(dataframe_calories,dataframe_prices,on=:item)

# ------------------------------------------------------------------------------------------
# ### FileIO
# ------------------------------------------------------------------------------------------

using FileIO
julialogo = download("https://avatars0.githubusercontent.com/u/743164?s=200&v=4","julialogo.png")

# ------------------------------------------------------------------------------------------
# Again, let's check that this download worked!
# ------------------------------------------------------------------------------------------

;ls

# ------------------------------------------------------------------------------------------
#
# Next, let's load the Julia logo, stored as a .png file
#
# **NOTE: You may see errors below, that certain Image packages could not be found. If so:**
#  - This is because these packages are specific to your OS, so aren't installed by default.
#  - Simply run the suggested commands to install the packages, then try again.
# ------------------------------------------------------------------------------------------

# These commands may vary depending on your operating system.
#import Pkg; Pkg.add("QuartzImageIO")
#import Pkg; Pkg.add("ImageMagick")

X1 = load("julialogo.png")

# ------------------------------------------------------------------------------------------
# We see here that Julia stores this logo as an array of colors.
# ------------------------------------------------------------------------------------------

@show typeof(X1);
@show size(X1);

# ------------------------------------------------------------------------------------------
# And if we load the Images package, it will display in Jupyter as an image:
# ------------------------------------------------------------------------------------------

using Images

X1

# ------------------------------------------------------------------------------------------
# ### File types
# In Julia, many file types are supported so you do not have to transfer a file you from
# another language to a text file before you read it.
#
# *Some packages that achieve this:*
# MAT CSV NPZ JLD FASTAIO
#
#
# Let's try using MAT to write a file that stores a matrix.
# ------------------------------------------------------------------------------------------

using MAT

A = rand(5,5)
matfile = matopen("densematrix.mat", "w")
write(matfile, "A", A)
close(matfile)

# ------------------------------------------------------------------------------------------
# Now try opening densematrix.mat with MATLAB!
# ------------------------------------------------------------------------------------------

newfile = matopen("densematrix.mat")
read(newfile,"A")

names(newfile)

close(newfile)

max_stars_repo_path: .nbexports/introductory-tutorials/broader-topics-and-ecosystem/intro-to-julia-for-data-science/1. Julia for Data Science - Data.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
module SNNs
import SpikingNeuralNets: AbstractSNN, memory, bias
using LightGraphs

export SNN, memory, bias

"""
A simple data-type for generalized spiking neural networks.

# Arguments
- `graph::G`: the graph structure of neurons (vertices) and weights (edges)
- `bias::Vector{Real}`: the activation (potential) bias of each neuron
- `I::Vector{Bool}`: specifies whether each neuron is an input neuron
- `O::Vector{Bool}`: specifies whether each neuron is an output neuron
- `m::Integer`: the length of configuration memory
- `S::Vector{Dict{ST, ST}}`: a memory S[neuron][time] of spike counts
- `V::Vector{Vector{VT}}`: a length m voltage memory V[neuron][time],
    where V[1] is the least recent voltage and V[end] is the current recent voltage
"""
struct SNN{VT, ST, G} <: AbstractSNN{VT, ST, G}
    graph::G
    bias::Vector{VT}
    I::Vector{Bool}
    O::Vector{Bool}
    m::ST
    S::Vector{Dict{ST, ST}}
    V::Vector{Vector{VT}}

    """
        SNN{VT,ST,G}(graph::G; bias::Vector{<:VT}=zeros(nv(graph)), m::Integer=1,
                     S::Vector{<:Dict{<:ST, <:ST}}=[Dict{ST, ST}() for _ in vertices(graph)],
                     V::Vector{Vector{<:VT}}=[zeros(VT, nv(graph)) for _ in 1:m])

    Create a new SNN{VT,ST,G} with the given graph and starting configuration.
    """
    function SNN{VT,ST,G}(graph::G; bias::Vector{<:VT}=zeros(VT, nv(graph)), m::Integer=1,
                          S::Vector{<:Dict{<:ST, <:ST}}=[Dict{ST, ST}() for _ in vertices(graph)],
                          V::Vector{<:Vector{<:VT}}=[zeros(VT, nv(graph)) for _ in 1:m]) where {VT<:Real, ST<:Integer, G<:AbstractGraph{ST}}
        if m < 1
            error("Configuration memory length (m) must be >= 1.")
        end

        # Input neurons are those without incoming connections
        I = map(n -> length(inneighbors(graph, n)) == 0, vertices(graph))
        # Output neurons are those without outgoing connections
        O = map(n -> length(outneighbors(graph, n)) == 0, vertices(graph))
        return new(graph, bias, I, O, m, S, V)
    end
end

"""
    SNN{VT,ST}(graph::G, bias::Vector{Real}=zeros(nv(graph)), m::Integer=1,
               S::Vector{Dict{<:ST, <:ST}}=[Dict{ST, ST}() for _ in vertices(graph)],
               VT::Vector{<:Vector{<:VT}}=[zeros(VT, nv(graph)) for _ in 1:m])

Create a new SNN{VT,ST,G} with the given graph and starting configuration.
Shorthand for the more explicit call SNN{VT,ST,G}(graph, m, S, V)
"""
function SNN{VT, ST}(graph::G; bias::Vector{VT}=zeros(VT, nv(graph)), m::Integer=1,
                     S::Vector{<:Dict{<:ST, <:ST}}=[Dict{ST, ST}() for _ in vertices(graph)],
                     V::Vector{<:Vector{<:VT}}=[zeros(VT, nv(graph)) for _ in 1:m]) where {VT<:Real, ST<:Integer, G<:AbstractGraph}
    return SNN{VT, ST, G}(graph; bias=bias, m=m, S=S, V=V)
end

"""
    memory(snn::AbstractSNN)::Integer

Return the size of the network `snn`'s voltage memory.
"""
function memory(snn::SNN{VT,ST})::ST where {VT<:Real, ST<:Integer}
    return snn.m
end

"""
    bias(snn::AbstractSNN)::Vector{<:Real}

Return the bias for all neurons in the network `snn`.
"""
function bias(snn::SNN{VT})::Vector{VT} where {VT<:Real}
    return snn.bias
end

end

ext: jl


text:
using VLKeggSDK

# test reaction -
reaction_phrase = "C00533 + C00014 + 3 C00126 <=> C05361 + C00001 + 3 C00125"

# go -
std = extract_stoichiometric_dictionary(reaction_phrase)


text:
shape(od::AbstractOddFrame) = [length(od.labels), length(od.columns[1])]
show(od::OddFrame) = od.head(5)


text:
"""
    Alias

Contains an alias for a GraphQL field or query. `Alias`es can be used in
the `output_fields` keyword argument as well as directly instead of query,
mutation and subscription names.

# Examples

Using an `Alias` instead of the query name:

```jldoctest alias_label; setup=:(using GraphQLClient)
julia> client = Client("https://countries.trevorblades.com");

julia> alias = Alias("country_alias", "country");

julia> query(client, alias, query_args=Dict("code"=>"BR"), output_fields=["name"]).data
Dict{String, Any} with 1 entry:
  "country_alias" => Dict{String, Any}("name"=>"Brazil")

```

Using an `Alias` in `output_fields`:

```jldoctest alias_label; setup=:(using GraphQLClient)
julia> field_alias = Alias("country_name_alias", "name");

julia> query(client, "country", query_args=Dict("code"=>"BR"), output_fields=[field_alias]).data
Dict{String, Any} with 1 entry:
  "country" => Dict{String, Any}("country_name_alias"=>"Brazil")
```
"""
struct Alias
    alias::String
    name::String
end
Base.show(io::IO, ::MIME"text/plain", a::Alias) = print(io, "Alias: ", a.alias, ":", a.name)
Base.show(io::IOBuffer, a::Alias) = print(io, a.alias * ":" * a.name) # for string interpolation

"""
    get_name(name::AbstractString)
    get_name(alias::Alias)

To be used when the GraphQL field name is required, which will
either be inputted as a string or in the `name` field of an
`Alias`.
"""
get_name(name::AbstractString) = name
get_name(alias::Alias) = alias.name

"""
    GQLEnum

When using `direct_write=true` in queries and mutations, ENUMs must be wrapped
in this type to ensure that they are not wrapped in quotes in the query string.

See [`directly_write_query_args`](@ref) for more information and examples.
"""
struct GQLEnum
    value::String
end


text:
function tfidf(X::SparseMatrixCSC)
    # item = document
    # user = term
    U, I, R = findnz(X)
    n_users, n_items = size(X)

    bincount = zeros(n_users)
    for u in U
        bincount[u] += 1
    end

    idf = log.(n_items ./ (bincount .+ 1))

    for j = 1:length(U)
        R[j] = R[j] * idf[U[j]]
    end

    return sparse(U, I, R)
end

function bm25(X::SparseMatrixCSC, k1 = 1.2, b = 0.75)
    U, I, R = findnz(X)
    n_users, n_items = size(X)

    bincount = zeros(n_users)
    for u in U
        bincount[u] += 1
    end

    idf = log.((n_items .- bincount .+ 0.5) ./ (bincount .+ 0.5))

    dl = zeros(n_items)
    for i in I
        dl[i] += 1
    end
    avgdl = sum(dl) / n_items

    for j = 1:length(U)
        R[j] = R[j] * idf[U[j]] * (k1 + 1) / (1 + k1 * (1 - b + b * dl[I[j]] / avgdl))
    end

    return sparse(U, I, R)
end

function compute_similarity(
    X::SparseMatrixCSC,
    topK::Int,
    shrink::Float64,
    normalize::Bool,
    normalize_similarity::Bool,
)
    # (user, item)^T * (user, item) -> (item, item)
    # Return S[i, j] where j is full items, and i is related items at topK
    if shrink < 0
        throw(ArgumentError("shrink must be 0 or positive."))
    end

    n_users, n_items = size(X)
    topK = min(topK, n_items - 1)

    simJ = Vector{Int64}(undef, topK * n_items)
    simI = Vector{Int64}(undef, topK * n_items)
    simS = Vector{Float64}(undef, topK * n_items)

    norms = sqrt.(sum(X .^ 2, dims = 1))
    norms = dropdims(norms, dims = 1)

    # to speed up, cache non zero indices
    nonzero_I_R = [findnz(X[u, :]) for u = 1:n_users]

    Threads.@threads for j = 1:n_items
        Uj, Rj = findnz(X[:, j])
        simj = zeros(n_items)
        for (u, ruj) in zip(Uj, Rj)
            Iu, Ri = nonzero_I_R[u]
            for (i, rui) in zip(Iu, Ri)
                s = rui * ruj
                if normalize
                    s /= norms[j] * norms[i] + shrink + 1e-6
                end
                simj[i] += s
            end
        end
        arg_sort_i = sortperm(simj, rev = true)[2:topK+1]
        simI[(1+(j-1)*topK):j*topK] = arg_sort_i
        simS[(1+(j-1)*topK):j*topK] = simj[arg_sort_i]
        simJ[(1+(j-1)*topK):j*topK] = fill(j, length(arg_sort_i))
    end

    similarity = sparse(simI, simJ, simS)

    if normalize_similarity
        # see M. Deshpande and G. Karypis (2004)
        # https://doi.org/10.1145/963770.963776
        for i = 1:size(similarity)[2]
            similarity[:, i] /= sum(similarity[:, i])
        end
    end
    return similarity
end

function predict_u2i(
    similarity::SparseMatrixCSC,
    user_history::Vector,
    n::Int64;
    drop_history::Bool = false,
)
    user_history = sparse(user_history)
    return predict_u2i(similarity, user_history, n; drop_history = drop_history)
end

function predict_u2i(
    similarity::SparseMatrixCSC,
    user_history::SparseVector,
    n::Int64;
    drop_history::Bool = false,
)
    pred_iidx, pred_score = findnz(similarity * user_history)
    permsocre = sortperm(pred_score, rev = true)
    pred = pred_iidx[permsocre]

    viewed_item, _ = findnz(user_history)
    # this is very slow
    if drop_history
        pred = filter(p -> !(p in viewed_item), pred)
    end
    n = min(n, length(pred))
    return pred[1:n]
end

function predict_i2i(similarity::SparseMatrixCSC, iidx::Int, n::Int)
    pred_iidx, pred_score = findnz(similarity[:, iidx])
    n = min(n, length(pred_iidx))
    return pred_iidx[sortperm(pred_score, rev = true)][1:n]
end


text:
module AuditoryFilters

using DSP
import Base: convert, filt
export erb_space, make_erb_filterbank, erb_filterbank, compute_modulation_cfs, make_modulation_filter, modulation_filterbank, gammatonegram,
	   ERBFilterbank, ModulationFilterbank

abstract Filterbank

nchannels(f::Filterbank) = length(f.filters)

function filt(fb::Filterbank, x)
	output = zeros(length(x), length(fb.filters))
	@inbounds begin
		for k = 1:length(fb.filters)
			output[:, k] = filt(fb.filters[k], x)
		end
	end
	return output
end

include("ERBFilterbank.jl")
include("ModulationFilterbank.jl")
include("Gammatonegram.jl")

end #module


text:
include("../common.jl")

name = "CXSparse"

sources = [
    GitSource("https://github.com/DrTimothyAldenDavis/SuiteSparse.git",
              "538273cfd53720a10e34a3d80d3779b607e1ac26"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/SuiteSparse
FLAGS+=(INSTALL="${prefix}" INSTALL_LIB="${libdir}" INSTALL_INCLUDE="${prefix}/include" CFOPENMP=)

if [[ ${target} == *mingw32* ]]; then
    FLAGS+=(UNAME=Windows)
    FLAGS+=(LDFLAGS="${LDFLAGS} -L${libdir} -shared")
else
    FLAGS+=(UNAME="$(uname)")
    FLAGS+=(LDFLAGS="${LDFLAGS} -L${libdir}")
fi

make -j${nproc} -C CXSparse "${FLAGS[@]}" library CFOPENMP="$CFOPENMP"
make -j${nproc} -C CXSparse "${FLAGS[@]}" install CFOPENMP="$CFOPENMP"

install_license LICENSE.txt
"""

platforms = supported_platforms(;experimental=true)

dependencies = [
    Dependency("SuiteSparse_jll")
]

products = [
    LibraryProduct("libcxsparse", :libcxsparse)
]
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
using Suppressor

@testset "Integration tests" begin
    #cubic lattice antiferromagnet in a magnetic field
    a1 = (1.0, 0.0, 0.0)
    a2 = (0.0, 1.0, 0.0)
    a3 = (0.0, 0.0, 1.0)
    uc = UnitCell(a1,a2,a3)
    b1 = addBasisSite!(uc, (0.0, 0.0, 0.0))
    addInteraction!(uc, b1, b1, [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0], (1, 0, 0))
    addInteraction!(uc, b1, b1, [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0], (0, 1, 0))
    addInteraction!(uc, b1, b1, [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0], (0, 0, 1))
    setField!(uc, b1, [0.0, 0.0, 0.2])
    cubiclattice = Lattice(uc, (4,4,4))

    m = MonteCarlo(cubiclattice, 2.0, 10000, 1000000, seed=UInt(0))
    @suppress run!(m)
    e,e2 = means(m.observables.energy)
    @test isapprox(e, -2.4850, rtol=1e-3)
    @test isapprox(e2, 6.1794, rtol=1e-3)

    #triangular lattice ferromagnet in a field
    a1 = (3/2, sqrt(3)/2)
    a2 = (3/2, -sqrt(3)/2)
    uc = UnitCell(a1,a2)
    b1 = addBasisSite!(uc, (0.0, 0.0))
    addInteraction!(uc, b1, b1, [-1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0], (1, 0))
    addInteraction!(uc, b1, b1, [-1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0], (0, 1))
    addInteraction!(uc, b1, b1, [-1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0], (1, -1))
    setField!(uc, b1, [0.0, 0.0, 0.5])
    triangularlattice = Lattice(uc, (8,8))

    m = MonteCarlo(triangularlattice, 2.0, 10000, 1000000, seed=UInt(0))
    @suppress run!(m)
    e,e2 = means(m.observables.energy)
    @test isapprox(e, -2.9806, rtol=1e-3)
    @test isapprox(e2, 8.8882, rtol=1e-3)

    #honeycomb lattice
    J1 = -1.1 #NN Heisenberg
    J2 = -0.3 #2NN Heisenberg
    K = -3.0 #NN Kitaev
    G = -0.5 #NN Gamma
    D = -0.6 #2NN Dzyaloshinskii-Moriya
    A = -0.8 #Onsite anisotropy

    M1x = [J1+K 0.0 0.0; 0.0 J1 G; 0.0 G J1]
    M1y = [J1 0.0 G; 0.0 J1+K 0.0; G 0.0 J1]
    M1z = [J1 G 0.0; G J1 0.0; 0.0 0.0 J1+K]
    M2 = [J2 D/sqrt(3.0) -D/sqrt(3.0); -D/sqrt(3.0) J2 D/sqrt(3.0); D/sqrt(3.0) -D/sqrt(3.0) J2]
    a1 = (3/2, sqrt(3)/2)
    a2 = (3/2, -sqrt(3)/2)
    uc = UnitCell(a1,a2)
    b1 = addBasisSite!(uc, (0.0, 0.0))
    b2 = addBasisSite!(uc, (1.0, 0.0))
    addInteraction!(uc, b1, b2, M1x, (0, 0))
    addInteraction!(uc, b1, b2, M1y, (0, -1))
    addInteraction!(uc, b1, b2, M1z, (-1, 0))
    addInteraction!(uc, b1, b1, M2, (1, 0))
    addInteraction!(uc, b1, b1, collect(transpose(M2)), (1, -1))
    addInteraction!(uc, b1, b1, M2, (0, -1))
    addInteraction!(uc, b2, b2, collect(transpose(M2)), (1, 0))
    addInteraction!(uc, b2, b2, M2, (1, -1))
    addInteraction!(uc, b2, b2, collect(transpose(M2)), (0, -1))
    setInteractionOnsite!(uc, b1, A*ones(3,3))
    setInteractionOnsite!(uc, b2, A*ones(3,3))
    honeycomblattice = Lattice(uc, (8,8))

    m = MonteCarlo(honeycomblattice, 0.5, 100000, 1000000, seed=UInt(0))
    @suppress run!(m)
    e,e2 = means(m.observables.energy)
    @test isapprox(e, -4.3855, rtol=1e-3)
    @test isapprox(e2, 19.2992, rtol=1e-3)
end


text:
module Output
import HDF5
using H5Zblosc
import Logging
import Base: getindex, show
using EllipsisNotation
import EllipsisNotation: Ellipsis
import Printf: @sprintf
import Luna: Scans, Utils, @hlock
import Pidfile: mkpidlock

abstract type AbstractOutput end

"Output handler for writing only to memory"
mutable struct MemoryOutput{sT, S} <: AbstractOutput
    save_cond::sT
    yname::AbstractString  # Name for solution (e.g. "Eω")
    tname::AbstractString  # Name for propagation direction (e.g. "z")
    saved::Integer  # How many points have been saved so far
    data::Dict{String, Any}  # The actual data
    statsfun::S  # Callable, returns dictionary of statistics
end

function MemoryOutput(tmin, tmax, saveN::Integer, statsfun=nostats;
                      yname="Eω", tname="z", script=nothing)
    save_cond = GridCondition(tmin, tmax, saveN)
    MemoryOutput(save_cond, yname, tname, statsfun, script)
end

function MemoryOutput(save_cond, yname, tname, statsfun=nostats, script=nothing)
    data = Dict{String, Any}()
    data["stats"] = Dict{String, Any}()
    data["meta"] = Dict{String, Any}()
    data["meta"]["sourcecode"] = Utils.sourcecode()
    data["meta"]["git_commit"] = Utils.git_commit()
    if !isnothing(script)
        data["meta"]["script_code"] = script
    end
    MemoryOutput(save_cond, yname, tname, 0, data, statsfun)
end

function initialise(o::MemoryOutput, y)
    dims = init_dims(size(y), o.save_cond)
    o.data[o.yname] = Array{ComplexF64}(undef, dims)
    o.data[o.tname] = Array{Float64}(undef, (dims[end],))
end

"getindex works interchangeably so when switching from one Output to
another, subsequent code can stay the same"
getindex(o::MemoryOutput, ds::AbstractString) = o.data[ds]
getindex(o::MemoryOutput, ds::AbstractString, I...) = o.data[ds][I...]

show(io::IO, o::MemoryOutput) = print(io, "MemoryOutput$(collect(keys(o.data)))")

"""Calling the output handler saves data in the arrays
    Arguments:
        y: current function value
        t: current propagation point
        dt: current stepsize
        yfun: callable which returns interpolated function value at different t
    Note that from RK45.jl, this will be called with yn and tn as arguments.
"""
function (o::MemoryOutput)(y, t, dt, yfun)
    save, ts = o.save_cond(y, t, dt, o.saved)
    append_stats!(o, o.statsfun(y, t, dt))
    !haskey(o.data, o.yname) && initialise(o, y)
    while save
        s = size(o.data[o.yname])
        if s[end] < o.saved+1
            o.data[o.yname] = fastcat(o.data[o.yname], yfun(ts))
            push!(o.data[o.tname], ts)
        else
            idcs = fill(:, ndims(y))
            o.data[o.yname][idcs..., o.saved+1] = yfun(ts)
            o.data[o.tname][o.saved+1] = ts
        end
        o.saved += 1
        save, ts = o.save_cond(y, t, dt, o.saved)
    end
end

function append_stats!(o::MemoryOutput, d)
    for (k, v) in pairs(d)
        append_stat!(o, k, v)
    end
end

function append_stat!(o::MemoryOutput, name, value::Number)
    if ~haskey(o.data["stats"], name)
        o.data["stats"][name] = [value]
    else
        push!(o.data["stats"][name], value)
    end
end

function append_stat!(o::MemoryOutput, name, value::AbstractArray)
    if ~haskey(o.data["stats"], name)
        dims = size(value)
        o.data["stats"][name] = reshape(value, (size(value)..., 1))
    else
        o.data["stats"][name] = fastcat(o.data["stats"][name], value)
    end
end

"Calling the output on a dictionary writes the items to the array"
function (o::MemoryOutput)(d::Dict; force=false, meta=false, group=nothing)
    for (k, v) in pairs(d)
        o(k, v; force=force, meta=meta, group=group)
    end
end

"Calling the output with a key, value pair writes the value to the array."
function (o::MemoryOutput)(key::AbstractString, val; force=false, meta=false, group=nothing)
    parent = meta ? o.data["meta"] : o.data
    if haskey(parent, key)
            if force
                Logging.@warn("Key $key already exists and will be overwritten.")
            else
                error("Key $key already present in dataset.")
        end
    end
    if !isnothing(group)
        if !haskey(parent, group)
            parent[group] = Dict{String, Any}()
        end
        parent[group][key] = val
    else
        parent[key] = val
    end
end

function tofile(fpath, o::MemoryOutput)
    Utils.save_dict_h5(fpath, o.data)
end

function fastcat(A, v)
    Av = vec(A)
    append!(Av, vec(v))
    dims = size(A)
    return reshape(Av, (dims[1:end-1]..., dims[end]+1))
end

"Output handler for writing to an HDF5 file"
mutable struct HDF5Output{sT, S} <: AbstractOutput
    fpath::AbstractString  # Path to output file
    save_cond::sT  # callable, determines when data is saved and where it is interpolated
    yname::AbstractString  # Name for solution (e.g. "Eω")
    tname::AbstractString  # Name for propagation direction (e.g. "z")
    saved::Integer  # How many points have been saved so far
    statsfun::S  # Callable, returns dictionary of statistics
    stats_tmp::Vector{Dict{String, Any}}  # Temporary storage for statistics between saves
    compression::Bool # whether to use compression
    cache::Bool # whether to cache latest solution point (for continuing after interrupt)
    cachehash::UInt64 # safety hash to prevent cache-continuing for different propagations
    readonly::Bool
end

"Simple constructor"
function HDF5Output(fpath, tmin, tmax, saveN::Integer, statsfun=nostats;
                    yname="Eω", tname="z", compression=false, script=nothing, cache=true,
                    readonly=false)
    save_cond = GridCondition(tmin, tmax, saveN)
    HDF5Output(fpath, save_cond, yname, tname, statsfun, compression, script, cache, readonly)
end

"Internal constructor - creates the file"
function HDF5Output(fpath, save_cond, yname, tname, statsfun, compression,
                    script=nothing, cache=true, readonly=false)
    if isfile(fpath) && cache
        @hlock HDF5.h5open(fpath, "cw") do file
            if HDF5.haskey(file["meta"], "cache")
                saved = read(file["meta"]["cache"]["saved"])
                chash = hash((sort(keys(file["stats"])), size(file[yname])[1:end-1]))
            else
                error("cached HDF5Output created, file exists, but has no cache")
            end
        end
    elseif !readonly
        if isfile(fpath)
            Logging.@warn("output file $(fpath) already exists and will be overwritten!")
            rm(fpath)
        end
        fdir, fname = splitdir(fpath)
        isdir(fdir) || mkpath(fdir)
        @hlock HDF5.h5open(fpath, "cw") do file
            HDF5.create_group(file, "stats")
            HDF5.create_group(file, "meta")
            file["meta"]["sourcecode"] = Utils.sourcecode()
            file["meta"]["git_commit"] = Utils.git_commit()
            if !isnothing(script)
                file["meta"]["script_code"] = script
            end
            if cache
                HDF5.create_group(file["meta"], "cache")
            end
        end
        chash = UInt64(0)
        saved = 0
    else
        chash = UInt64(0)
        saved = 0
    end
    stats0 = Vector{Dict{String, Any}}()
    HDF5Output(fpath, save_cond, yname, tname, saved, statsfun, stats0,
               compression, cache, chash, readonly)
end

function HDF5Output(fpath::AbstractString)
    isfile(fpath) || error("Cannot open read-only HDF5Output: file not found")
    HDF5Output(fpath, 0, 0, 1; readonly=true)
end

function initialise(o::HDF5Output, y)
    ydims = size(y)
    idims = init_dims(ydims, o.save_cond)
    cdims = collect(idims)
    dims = Tuple(cdims)
    chdims = (dims[1:end-1]..., 1) # Chunk size is that of one z-point
    mdims = copy(cdims)
    mdims[end] = -1
    maxdims = Tuple(mdims)
    @hlock HDF5.h5open(o.fpath, "r+") do file
        if o.compression
            HDF5.create_dataset(file, o.yname, HDF5.datatype(ComplexF64), (dims, maxdims),
                          chunk=chdims, blosc=3)
        else
            HDF5.create_dataset(file, o.yname, HDF5.datatype(ComplexF64), (dims, maxdims),
                          chunk=chdims)
        end
        HDF5.create_dataset(file, o.tname, HDF5.datatype(Float64), ((dims[end],), (-1,)),
                      chunk=(1,))
        statsnames = sort(collect(keys(o.stats_tmp[end])))
        o.cachehash = hash((statsnames, size(y)))
        file["meta"]["cachehash"] = o.cachehash
        if o.cache
            file["meta"]["cache"]["t"] = typemin(0.0)
            file["meta"]["cache"]["dt"] = typemin(0.0)
            file["meta"]["cache"]["y"] = y
            file["meta"]["cache"]["saved"] = 0
        end
    end
end

# for single String index, read whole data set
function getindex(o::HDF5Output, idx::AbstractString)
    @hlock HDF5.h5open(o.fpath, "r") do file
        read(file[idx])
    end
end

# more indices -> read slice of data
function getindex(o::HDF5Output, ds::AbstractString,
                  I::Union{AbstractRange, Int, Colon, Ellipsis}...)
    @hlock HDF5.h5open(o.fpath, "r") do file
        file[ds][to_indices(file[ds], I)...]
    end
end

# indexing with an array, e.g. o["Eω", :, [1, 2, 3]] has to be handled separately
function getindex(o::HDF5Output, ds::AbstractString,
                  I::Union{AbstractRange, Int, Colon, Array, Ellipsis}...)
    if count(isa.(I, Array)) > 1
        error("Only one dimension can be index with an array.")
    end
    @hlock HDF5.h5open(o.fpath, "r") do file
        dset = file[ds]
        idcs = to_indices(dset, I)
        adim = findfirst(isa.(idcs, Array)) # which of the indices is the array
        arr = idcs[adim] # the array itself
        T = eltype(dset)
        ret = Array{T}(undef, map(length, idcs))
        Ilo = idcs[1:adim-1]
        Ihi = idcs[adim+1:end]
        for ii in eachindex(arr)
            ret[Ilo..., ii, Ihi...] .= dset[Ilo..., arr[ii], Ihi...]
        end
        ret
    end
end

function show(io::IO, o::HDF5Output)
    if isfile(o.fpath)
        fields = @hlock HDF5.h5open(o.fpath) do file
            keys(file)
        end
        print(io, "HDF5Output$(fields)")
    else
        print(io, "HDF5Output[FILE DELETED]")
    end
end

"""Calling the output handler writes data to the file
    Arguments:
        y: current function value
        t: current propagation point
        dt: current stepsize
        yfun: callable which returns interpolated function value at different t
    Note that from RK45.jl, this will be called with yn and tn as arguments.
"""
function (o::HDF5Output)(y, t, dt, yfun)
    o.readonly && error("Cannot add data to read-only output!")
    save, ts = o.save_cond(y, t, dt, o.saved)
    push!(o.stats_tmp, o.statsfun(y, t, dt))
    if save
        @hlock HDF5.h5open(o.fpath, "r+") do file
            !HDF5.haskey(file, o.yname) && initialise(o, y)
            statsnames = sort(collect(keys(o.stats_tmp[end])))
            cachehash = hash((statsnames, size(y)))
            cachehash == o.cachehash || error(
                "the hash for this propagation does not agree with cache in file")
            while save
                s = collect(size(file[o.yname]))
                idcs = fill(:, length(s)-1)
                if s[end] < o.saved+1
                    s[end] += 1
                    HDF5.set_extent_dims(file[o.yname], Tuple(s))
                end
                file[o.yname][idcs..., o.saved+1] = yfun(ts)
                s = collect(size(file[o.tname]))
                if s[end] < o.saved+1
                    s[end] += 1
                    HDF5.set_extent_dims(file[o.tname], Tuple(s))
                end
                file[o.tname][o.saved+1] = ts
                o.saved += 1
                save, ts = o.save_cond(y, t, dt, o.saved)
            end
            append_stats!(file["stats"], o.stats_tmp)
            o.stats_tmp = Vector{Dict{String, Any}}()
            if o.cache
                write(file["meta"]["cache"]["t"], t)
                write(file["meta"]["cache"]["dt"], dt)
                write(file["meta"]["cache"]["y"], y)
                write(file["meta"]["cache"]["saved"], o.saved)
            end
        end
    end
end

function append_stats!(parent, a::Array{Dict{String,Any},1})
    N = length(a)
    names = HDF5.keys(parent)
    for (k, v) in pairs(a[1])
        if ~(k in names)
            create_dataset(parent, k, v)
        end
        s = collect(size(parent[k]))
        curN = s[end]
        if ~(k in names)
            curN -= 1 # new dataset - overwrite initial value
        end
        s[end] += N
        if ~(k in names)
            s[end] -= 1 # new dataset - overwrite initial value
        end
        HDF5.set_extent_dims(parent[k], Tuple(s))
        for ii = 1:N
            parent[k][fill(:, ndims(a[ii][k]))..., curN+ii] = a[ii][k]
        end
    end
end

function create_dataset(parent, name, x::Number)
    @hlock HDF5.create_dataset(parent, name, HDF5.datatype(typeof(x)), ((1,), (-1,)),
                  chunk=(1,))
end

function create_dataset(parent, name, x::AbstractArray)
    dims = (size(x)..., 1)
    maxdims = (size(x)..., -1)
    @hlock HDF5.create_dataset(parent, name, HDF5.datatype(eltype(x)), (dims, maxdims),
                  chunk=dims)
end

"Calling the output on a dictionary writes the items to the file"
function (o::HDF5Output)(d::AbstractDict; force=false, meta=false, group=nothing)
    o.readonly && error("Cannot add data to read-only output!")
    @hlock HDF5.h5open(o.fpath, "r+") do file
        parent = meta ? file["meta"] : file
        for (k, v) in pairs(d)
            if HDF5.haskey(parent, k)
                if force
                    Logging.@warn("Dataset $k already present in file $(o.fpath)"*
                                  " and will be overwritten")
                    HDF5.delete_object(parent, k)
                else
                    Logging.@warn("File $(o.fpath) already has dataset $k. Pass force=true"*
                              " to overwrite")
                end
            end
            isa(v, BitArray) && (v = Array{Bool, 1}(v))
            if !isnothing(group)
                if !HDF5.haskey(parent, group)
                    HDF5.create_group(parent, group)
                end
                if HDF5.haskey(parent[group], k)
                    write(parent[group][k], v)
                else
                    parent[group][k] = v
                end
            else
                if HDF5.haskey(parent, k)
                    write(parent[k], v)
                else
                    parent[k] = v
                end
            end
        end
    end
end

"Calling the output on a key, value pair writes the value to the file"
function (o::HDF5Output)(key::AbstractString, val; force=false, meta=false, group=nothing)
    o.readonly && error("Cannot add data to read-only output!")
    @hlock HDF5.h5open(o.fpath, "r+") do file
        parent = meta ? file["meta"] : file
        if HDF5.haskey(parent, key)
            if force
                Logging.@warn("Dataset $key already present in file $(o.fpath)"*
                              " and will be overwritten")
                HDF5.delete_object(parent, key)
            else
                Logging.@warn("File $(o.fpath) already has dataset $(key). Pass force=true"*
                              " to overwrite")
            end
        end
        isa(val, BitArray) && (val = Array{Bool, 1}(val))
        if !isnothing(group)
            if !HDF5.haskey(parent, group)
                HDF5.create_group(parent, group)
            end
            if HDF5.haskey(parent[group], key)
                write(parent[group][key], val)
            else
                parent[group][key] = val
            end
        else
            if HDF5.haskey(parent, key)
                write(parent[key], val)
            else
                parent[key] = val
            end
        end
    end
end

"""
    check_cache(o::HDF5Output, y, t, dt)

Check for an existing cached propagation in the output `o` and return this cache if present.
"""
function check_cache(o::HDF5Output, y, t, dt)
    if !o.cache || !haskey(o["meta"]["cache"], "t")
        return y, t, dt
    end
    tc = o["meta"]["cache"]["t"]
    if tc < t
        return y, t, dt
    end
    yc = o["meta"]["cache"]["y"]
    dtc = o["meta"]["cache"]["dt"]
    return yc, tc, dtc
end

# For other outputs (e.g. MemoryOutput or another function), checking the cache does nothing.
check_cache(o, y, t, dt) = y, t, dt

"Condition callable that distributes save points evenly on a grid"
struct GridCondition
    grid::Vector{Float64}
    saveN::Integer
end

function GridCondition(tmin, tmax, saveN)
    GridCondition(range(tmin, stop=tmax, length=saveN), saveN)
end

function (cond::GridCondition)(y, t, dt, saved)
    save = (saved < cond.saveN) && cond.grid[saved+1] <= t
    return save, save ? cond.grid[saved+1] : 0
end

"Condition which saves every native point of the propagation"
function always(y, t, dt, saved)
    return true, t
end

"Condition which saves every nth native point"
function every_nth(n)
    i = 0
    cond = let i = i, n = n
        function condition(y, t, dt, saved)
            save = i % n == 0
            i += 1
            return save, t
        end
    end
    return cond
end

"""Making initial array dimensions.
For a GridCondition, we know in advance how many points there will be.
"""
function init_dims(ydims, save_cond::GridCondition)
    return (ydims..., save_cond.saveN)
end

"For other conditions, we do not know in advance."
function init_dims(ydims, save_cond)
    return (ydims..., 1)
end

function nostats(args...)
    return Dict{String, Any}()
end

"""
    ScanHDF5Output(scan, scanidx, args...; fname=nothing, fdir=nothing, kwargs...)

Create an [`HDF5Output`](@ref) for the given `scan` at the current `scanidx` and automatically
save the scan arrays and current values of the scan variables in the file. If given,
`fdir` is used as a directory in which to store the scan output. `fname` can be used to
manually name files. The running scan index will be appended to `fname` for each file.
"""
function ScanHDF5Output(scan, scanidx, args...; fname=nothing, fdir=nothing, kwargs...)
    fpath = isnothing(fname) ?
        Scans.makefilename(scan, scanidx) :
        Scans.makefilename(fname, scanidx)
    if !isnothing(fdir)
        fpath = joinpath(fdir, fpath)
        if !isdir(fdir)
            mkpath(fdir)
        end
    end
    savescan(HDF5Output(fpath, args...; kwargs...), scan, scanidx)
end

function ScanMemoryOutput(scan, scanidx, args...; kwargs...)
    savescan(MemoryOutput(args...; kwargs...), scan, scanidx)
end

function savescan(out, scan, scanidx)
    out("scanidx", scanidx, meta=true)
    vars = Dict{String, Any}()
    arrays = Dict{String, Any}()
    order = String[]
    shape = Int[] # scan shape
    # create grid of scan points
    for (var, arr) in zip(scan.variables, scan.arrays)
        push!(order, string(var))
        push!(shape, length(arr))
        arrays[string(var)] = arr
        vars[string(var)] = Scans.getvalue(scan, var, scanidx)
    end
    out(vars; meta=true, group="scanvars")
    out(arrays; meta=true, group="scanarrays")
    out("scanorder", order; meta=true)
    out("scanshape", shape; meta=true)
    out
end

"""
    @ScanHDF5Output(scan, scanidx, args...)

Create an [`HDF5Output`](@ref) for the given `scan` at the current `scanidx` and automatically
save the running script, scan arrays, and current values of the scan variables in the file.
All arguments, including keyword arguments, after `scanidx` are identical to `HDF5Output`.
"""
macro ScanHDF5Output(scan, scanidx, args...)
    code = ""
    try
        script = string(__source__.file)
        code = open(script, "r") do file
            read(file, String)
        end
    catch
    end
    for arg in args
        if isa(arg, Expr) && arg.head == :(=)
            arg.head = :kw
        end
    end
    ex = Expr(:call, :ScanHDF5Output, esc(scan), esc(scanidx))
    for arg in args
        push!(ex.args, esc(arg))
    end
    push!(ex.args, Expr(:kw, :script, code))
    ex
end

# Auto-generate @MemoryOutput and @HDF5Output macros
for op in (:MemoryOutput, :HDF5Output)
    eval(
        quote
            macro $op(args...)
                script_code = ""
                try
                    script = string(__source__.file)
                    code = open(script, "r") do file
                        read(file, String)
                    end
                    script_code = script*"\n"*code
                catch
                end
                for arg in args
                    if isa(arg, Expr) && arg.head == :(=)
                        arg.head = :kw
                    end
                end
                exp = Expr(:call, $op)
                for arg in args
                    push!(exp.args, esc(arg))
                end
                push!(exp.args, Expr(:kw, :script, script_code))
                exp
            end
        end
    )
end

"""
    scansave(scan, scanidx; grid, stats, fpath, script, kwargs...)

While running the given `scan`, save the variables given as keyword arguments into the scan grid as determined from the variables of the `scan`. Special keyword arguments are:

- `grid::AbstractGrid`: Save the simulation grid in a dictionary (but only once)
- `stats`: The statistics dictionary from a simulation is saved in scan grid and `NaN`-padded to account for variable lengths in the output arrays
- `fpath`: Path to the file. Defaults to the scan name plus "_collected"
- `script`: Path to the Julia scrpt file running the scan. Can be grabbed automatically using the macro [`@scansave`](@ref).
"""
function scansave(scan, scanidx; stats=nothing, fpath=nothing,
                                 grid=nothing, script=nothing, kwargs...)
    fpath = isnothing(fpath) ? "$(scan.name)_collected.h5" : fpath
    lockpath = joinpath(Utils.cachedir(), "$(basename(fpath)).scanlock")
    isdir(Utils.cachedir()) || mkpath(Utils.cachedir())
    pidlock = mkpidlock(lockpath)
    if !isfile(fpath)
        # First save - set up file structure
        @hlock HDF5.h5open(fpath, "cw") do file
            group = HDF5.create_group(file, "scanvariables")
            order = String[]
            shape = Int[] # scan shape
            # create grid of scan points
            for (k, var) in zip(scan.variables, scan.arrays)
                group[string(k)] = var
                push!(order, string(k))
                push!(shape, length(var))
            end
            file["scanorder"] = order
            if !isnothing(stats)
                group = HDF5.create_group(file, "stats")
                for (k, v) in pairs(stats)
                    dims = (size(v)..., shape...)
                    #= last dimension of a statistics array is number of steps,
                        so save in chunks of 100 steps =#
                    fixeddims_v = size(v)[1:end-1]
                    mdims = (fixeddims_v..., -1, shape...)
                    chdims = (fixeddims_v..., 100, fill(1, length(shape))...)
                    HDF5.create_dataset(group, k, HDF5.datatype(eltype(v)), (dims, mdims),
                                  chunk=chdims)
                end
                group["valid_length"] = zeros(Int, shape...)
            end
            if !isnothing(grid)
                group = HDF5.create_group(file, "grid")
                for (k, v) in pairs(grid)
                    isa(v, BitArray) && (v = Array{Bool, 1}(v))
                    group[k] = v
                end
            end
            if !isnothing(script)
                script_code = ""
                try
                    code = open(script, "r") do file
                        read(file, String)
                    end
                    script_code = script*"\n"*code
                catch
                end
                file["script"] = script_code
            end
            # deal with other keyword arguments (additional quantities to be saved)
            for (k, v) in kwargs
                # dimensions of the array
                dims = (size(v)..., shape...)
                # chunk size is dimension of one array
                chdims = (size(v)..., fill(1, length(shape))...)
                HDF5.create_dataset(file, string(k), HDF5.datatype(eltype(v)), (dims, dims),
                              chunk=chdims)
            end
        end
    end
    @hlock HDF5.h5open(fpath, "r+") do file
        scanshape = Tuple([length(ai) for ai in scan.arrays])
        cidcs = CartesianIndices(scanshape)
        scanidcs = Tuple(cidcs[scanidx])
        if !isnothing(stats)
            for (k, v) in pairs(stats)
                if size(v)[end] > size(file["stats"][k])[ndims(v)]
                    #= new point has more stats points than before - extend dataset
                        stats arrays are of shape (N1, N2,... Ns) where N1 etc are fixed and
                        Ns depends on the number of steps
                        stats *datasets" have shape (N1, N2,... Ns, Nx, Ny,...) where Nx, Ny...
                        are the lengths of the scan arrays (see scanshape above)=#
                    oldlength = size(file["stats"][k])[ndims(v)] # current Ns
                    newlength = size(v)[end] # new Ns
                    newdims = (size(v)..., scanshape...) # (N1, N2,..., new Ns, Nx, Ny,...)
                    HDF5.set_extent_dims(file["stats"][k], newdims) # set new dimensions
                    # For existing shorter arrays, fill everything above their length with NaN
                    nanidcs = (fill(:, (ndims(v)-1))..., oldlength+1:newlength)
                    allscan = fill(:, length(scanshape)) # = (1:Nx, 1:Ny,...)
                    file["stats"][k][nanidcs..., allscan...] = NaN
                end
                # Stats array has shape (N1, N2,... Ns) - fill everything up to Ns with data
                sidcs = (fill(:, (ndims(v)-1))..., 1:size(v)[end])
                file["stats"][k][sidcs..., scanidcs...] = v
                # save number of valid points we just saved
                file["stats"]["valid_length"][scanidcs...] = size(v)[end]
                # fill everything after Ns with NaN
                nanidcs = (fill(:, (ndims(v)-1))..., size(v)[end]+1:size(file["stats"][k])[ndims(v)])
                file["stats"][k][nanidcs..., scanidcs...] = NaN
            end
        end
        for (k, v) in pairs(kwargs)
            sidcs = fill(:, ndims(v))
            file[string(k)][sidcs..., scanidcs...] = v
        end
    end
    close(pidlock)
end

"""
    @scansave(scan, scanidx; kwargs...)

Like [`scansave`](@ref) but also saves the script being run automatically.
"""
macro scansave(scan, scanidx, kwargs...)
    global script = string(__source__.file)
    ex = :(scansave($(esc(scan)), $(esc(scanidx)),
                    script=script))
    for arg in kwargs
        if isa(arg, Expr) && arg.head == :(=)
            arg.head = :kw
            push!(ex.args, esc(arg))
        else
            # To a macro, arguments and keyword arguments look the same, so check manually
            error("@scansave only accepts keyword arguments")
        end
    end
    ex
end
end


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

module DefineLog

export addObservers

using AdaptiveStressTesting
using SISLES
using SISLES: Encounter, EncounterDBN, WorldModel, Sensor, CollisionAvoidanceSystem,
    DynamicModel, PilotResponse, GenerativeModel
using RLESUtils, Loggers

using ..DefineSave

const ENABLE_ROUNDING = true
const ROUND_NDECIMALS = 9

function addObservers(sim::ACASX_GM)
    clearObservers!(sim)
    log = TaggedDFLogger()
    #separate folders for each aircraft.  logs can be different.  e.g. 2 different CAS
    for i = 1:sim.params.num_aircraft
        #data
        add_folder!(log, "CAS_info_$i", cas_info_types(sim.cas[i]),  cas_info_names(sim.cas[i]))
        add_folder!(log, "Initial_$i", initial_types(sim.em),  initial_names(sim.em))
        add_folder!(log, "Command_$i", command_types(sim.em.output_commands[i]),
            command_names(sim.em.output_commands[i]))
        add_folder!(log, "Sensor_$i", sensor_types(sim.sr[i]),  sensor_names(sim.sr[i]))
        add_folder!(log, "CAS_$i", cas_types(sim.cas[i]),  cas_names(sim.cas[i]))
        add_folder!(log, "Response_$i", response_types(sim.pr[i]),  response_names(sim.pr[i]))
        add_folder!(log, "Dynamics_$i", dynamics_types(sim.dm[i]),  dynamics_names(sim.dm[i]))
        add_folder!(log, "WorldModel_$i", worldmodel_types(sim.wm),  worldmodel_names(sim.wm))

        #units
        names = cas_info_names(sim.cas[i])
        units = cas_info_units(sim.cas[i])
        add_folder!(log, "CAS_info_units_$i", fill(String, length(units)), names)
        push!(log, "CAS_info_units_$i", units)
        names = initial_names(sim.em)
        units = initial_units(sim.em)
        add_folder!(log, "Initial_units_$i", fill(String, length(units)), names)
        push!(log, "Initial_units_$i", units)
        names = command_names(sim.em.output_commands[i])
        units = command_units(sim.em.output_commands[i])
        add_folder!(log, "Command_units_$i", fill(String, length(units)), names)
        push!(log, "Command_units_$i", units)
        names = sensor_names(sim.sr[i])
        units = sensor_units(sim.sr[i])
        add_folder!(log, "Sensor_units_$i", fill(String, length(units)), names)
        push!(log, "Sensor_units_$i", units)
        names = cas_names(sim.cas[i])
        units = cas_units(sim.cas[i])
        add_folder!(log, "CAS_units_$i", fill(String, length(units)), names)
        push!(log, "CAS_units_$i", units)
        names = response_names(sim.pr[i])
        units = response_units(sim.pr[i])
        add_folder!(log, "Response_units_$i", fill(String, length(units)), names)
        push!(log, "Response_units_$i", units)
        names = dynamics_names(sim.dm[i])
        units = dynamics_units(sim.dm[i])
        add_folder!(log, "Dynamics_units_$i", fill(String, length(units)), names)
        push!(log, "Dynamics_units_$i", units)
        names = worldmodel_names(sim.wm)
        units = worldmodel_units(sim.wm)
        add_folder!(log, "WorldModel_units_$i", fill(String, length(units)), names)
        push!(log, "WorldModel_units_$i", units)
    end
    #data
    add_folder!(log, "logProb", logprob_types(), logprob_names())
    add_folder!(log, "run_info", run_info_types(), run_info_names())

    #units
    names = logprob_names()
    units = logprob_units()
    add_folder!(log, "logProb_units", fill(String, length(units)), names)
    push!(log, "logProb_units", units)
    names = run_info_names()
    units = run_info_units()
    add_folder!(log, "run_info_units", fill(String, length(units)), names)
    push!(log, "run_info_units", units)

    addObserver(sim, "CAS_info",   x->log_cas_info!(log, x))
    addObserver(sim, "Initial",    x->log_initial!(log, x))
    addObserver(sim, "Command",    x->log_command!(log, x))
    addObserver(sim, "Sensor",     x->log_sensor!(log, x))
    addObserver(sim, "CAS",        x->log_cas!(log, x))
    addObserver(sim, "Response",   x->log_response!(log, x))
    addObserver(sim, "Dynamics",   x->log_adm!(log, x))
    addObserver(sim, "WorldModel", x->log_wm!(log, x))
    addObserver(sim, "logProb",   x->log_logProb!(log, x))
    addObserver(sim, "run_info",   x->log_run_info!(log, x))

    TrajLog(log)
end

function log_cas_info!(log::TaggedDFLogger, args)
    #[CAS version]
    i, cas = args
    push!(log, "CAS_info_$i", cas_info_data(cas))
end

cas_info_names(cas::Union{ACASX_CCAS,ACASX_ADD}) = String["version"]
cas_info_types(cas::Union{ACASX_CCAS,ACASX_ADD}) = [String]
cas_info_units(cas::Union{ACASX_CCAS,ACASX_ADD}) = String["n/a"]
cas_info_data(cas::Union{ACASX_CCAS,ACASX_ADD}) = Any[cas.version]

function log_initial!(log::TaggedDFLogger, args)
    #[aircraft_number, time_index, initial]
    i, t, aem = args
    @assert t == 0
    push!(log, "Initial_$i", initial_data(i, aem))
end

initial_names(aem::CorrAEMDBN) = String["v", "x", "y", "z", "psi", "theta", "phi", "v_d"]
initial_types(aem::CorrAEMDBN) = [Float64,Float64,Float64,Float64, Float64,Float64,
    Float64,Float64]
initial_units(aem::CorrAEMDBN) = String["ft/s", "ft", "ft", "ft", "deg", "deg", "deg", "ft/s^2"]

function initial_data(i, aem::CorrAEMDBN)
    # initial
    # airspeed, ft/s, double
    # north position, ft, double
    # east position, ft, double
    # altitude, ft, double
    # heading angle, degrees, double
    # flight path angle, degrees, double
    # roll angle, degrees, double
    # airspeed acceleration, ft/s^2, double

    t, x, y, h, v, psi = aem.aem.dynamic_states[i, 1, :]
    t_n, x_n, y_n, h_n, v_n, psi_n, = aem.aem.dynamic_states[i, 2, :]

    t, v_d, h_d, psi_d = aem.aem.states[i, 1, :]

    theta = atand((h_n - h) / norm(x_n - x, y_n - y)) |> to_plusminus_180
    phi = 0.0

    return round_floats(Any[v, x, y, h, psi, theta, phi, v_d], ROUND_NDECIMALS,
                      enable = ENABLE_ROUNDING)
end

initial_names(aem::Union{StarDBN,SideOnDBN}) = String["v", "x", "y", "z", "psi", "theta", "phi", "v_d"]
initial_types(aem::Union{StarDBN,SideOnDBN}) = [Float64,Float64,Float64,Float64,Float64,Float64,
    Float64,Float64]
initial_units(aem::Union{StarDBN,SideOnDBN}) = String["ft/s", "ft", "ft", "ft", "deg", "deg", "deg", "ft/s^2"]

function initial(i, aem::Union{StarDBN,SideOnDBN})
    # initial
    # airspeed, ft/s, double
    # north position, ft, double
    # east position, ft, double
    # altitude, ft, double
    # heading angle, degrees, double
    # flight path angle, degrees, double
    # roll angle, degrees, double
    # airspeed acceleration, ft/s^2, double

    is = aem.initial_states[i]
    t, x, y, h, v, psi = is.t, is.x, is.y, is.h, is.v, is.psi

    L, v_d, h_d, psi_d = aem.initial_commands[i]

    theta = atand(h_d / v)
    phi = 0.0

    return round_floats(Any[v, x, y, h, psi, theta, phi, v_d],
        ROUND_NDECIMALS, enable=ENABLE_ROUNDING)
end

function log_command!(log::TaggedDFLogger, args)
    #[aircraft_number, time_index, command]
    i, t, command = args
    push!(log, "Command_$i", command_data(t, command))
end

command_names(command::CorrAEMCommand) = String["t", "h_d", "v_d", "psi_d"]
command_types(command::CorrAEMCommand) = [Int64, Float64, Float64, Float64]
command_units(command::CorrAEMCommand) = String["index", "ft/s", "ft/s^2", "deg/s"]

function command_data(t::Int64, command::CorrAEMCommand)
    round_floats(Any[t, command.h_d, command.v_d, command.psi_d], ROUND_NDECIMALS,
        enable=ENABLE_ROUNDING)
end

function log_adm!(log::TaggedDFLogger, args)
    #[aircraft_number, time_index, dynamics model]
    i, t, adm = args
    push!(log, "Dynamics_$i", dynamics_data(t, adm))
end

dynamics_names(adm::SimpleADM) = String["t", "x", "y", "h", "v", "psi"]
dynamics_types(adm::SimpleADM) = [Int, Float64, Float64, Float64, Float64, Float64]
dynamics_units(adm::SimpleADM) = String["s", "ft", "ft", "ft", "ft/s", "deg"]

function dynamics_data(t::Int64, adm::SimpleADM)
    s = adm.state
    round_floats(Any[t, s.x, s.y, s.h, s.v, s.psi], ROUND_NDECIMALS, enable=ENABLE_ROUNDING)
end

dynamics_names(adm::LLADM) = String["t", "v", "N", "E", "h", "psi", "theta", "phi", "a"]
dynamics_types(adm::LLADM) = [Int, Float64, Float64, Float64, Float64, Float64, Float64, Float64,
    Float64]
dynamics_units(adm::LLADM) = String["s", "ft/s", "ft", "ft", "ft", "rad", "rad", "rad", "ft/s^2"]

function dynamics_data(t::Int64, adm::LLADM)
    s = adm.state
    round_floats(Any[t, s.v, s.N, s.E, s.h, s.psi, s.theta, s.phi, s.a], ROUND_NDECIMALS,
        enable=ENABLE_ROUNDING)
end

function log_wm!(log::TaggedDFLogger, args)
    #[time_index, world model]
    t, wm = args
    for i = 1:wm.number_of_aircraft
        push!(log, "WorldModel_$i", worldmodel_data(t, wm, i)) #[t, x,y,z,vx,vy,vz]
    end
end

worldmodel_names(wm::AirSpace) = String["t", "x", "y", "z", "vx", "vy", "vz"]
worldmodel_types(wm::AirSpace) = [Int, Float64, Float64, Float64, Float64, Float64, Float64]
worldmodel_units(wm::AirSpace) = String["s", "ft", "ft", "ft", "ft/s", "ft/s", "ft/s"]

function worldmodel_data(t::Int64, wm::AirSpace, aircraft_number::Int)
    s = wm.states[aircraft_number]
    round_floats(Any[t, s.x, s.y, s.h, s.vx, s.vy, s.vh], ROUND_NDECIMALS, enable=ENABLE_ROUNDING)
end

function log_sensor!(log::TaggedDFLogger, args)
    #[aircraft_number, time_index, cas]
    i, t, sr = args
    push!(log, "Sensor_$i", sensor_data(t, sr))
end

sensor_names(sr::Void) = String["t"]
sensor_types(sr::Void) = Type[Int64]
sensor_units(sr::Void) = String["index"]
sensor_data(t::Int64, sr::Void) = Any[] #input=[empty]

function sensor_names(sr::SimpleTCASSensor)
    v = String["t", "r", "r_d", "a", "a_d", "num_intruders"]
    for i = 1:length(sr.output.h)
        push!(v,"h_$i", "h_d_$i")
    end
    v
end
function sensor_types(sr::SimpleTCASSensor)
    v = [Int, Float64, Float64, Float64, Float64, Int]
    for i = 1:length(sr.output.h)
        push!(v, Float64, Float64)
    end
    v
end
function sensor_units(sr::SimpleTCASSensor)
    v = String["index", "ft", "ft/s", "ft", "ft/s", "num"]
    for i = 1:length(sr.output.h)
        push!(v, "ft", "ft/s")
    end
    v
end
function sensor_data(t::Int64, sr::SimpleTCASSensor)
    #input=[r,r_d, a, a_d],[num_intruders],[h,h_d for each intruder]
    out = sr.output
    num_intruders = length(out.h)
    hhd = Float64[]

    for (h, h_d) in zip(out.h, out.h_d)
        push!(hhd, h, h_d)
    end
    round_floats(vcat(Any[t, out.r, out.r_d, out.a, out.a_d],
                           num_intruders,
                           hhd), ROUND_NDECIMALS, enable = ENABLE_ROUNDING)
end

function sensor_names(sr::ACASXSensor)
    v = String["t", "dz", "z", "psi", "h", "modes", "num_intruders"]
    for i = 1:length(sr.output.intruders)
        v = vcat(v, String["valid_$i", "id_$i", "modes_$i", "sr_$i", "chi_$i", "z_$i",
            "cvc_$i", "vrc_$i", "vsb_$i"])
    end
    v
end
function sensor_types(sr::ACASXSensor)
    v = [Int64, Float64, Float64, Float64, Float64, Int64, Int64]
    for i = 1:length(sr.output.intruders)
        v = vcat(v, [Bool, Int64, Int64, Float64, Float64, Float64, Int64, Int64, Int64])
    end
    v
end
function sensor_units(sr::ACASXSensor)
    v = String["index", "ft/s", "ft", "rad", "ft", "n/a", "num"]
    for i = 1:length(sr.output.intruders)
        push!(v, "bool", "num", "n/a", "ft", "rad", "ft", "n/a", "n/a", "n/a")
    end
    v
end

function sensor_data(t::Int64, sr::ACASXSensor)
    #input=[dz,z,psi,h,modes],[num_intruders],[valid,id,modes,
        #sr,chi,z,cvc,vrc,vsb for each intruder]
    out = sr.output
    own = out.ownInput
    num_intruders = length(out.intruders)
    intr_vec = Any[]
    for intr in out.intruders
        push!(intr_vec, intr.valid, intr.id, intr.modes, intr.sr, intr.chi, intr.z,
            intr.cvc, intr.vrc, intr.vsb)
    end
    round_floats(vcat(Any[t, own.dz, own.z, own.psi, own.h, own.modes], num_intruders,
        intr_vec), ROUND_NDECIMALS, enable=ENABLE_ROUNDING)
end

function log_cas!(log::TaggedDFLogger, args)
    #[aircraft_number, time_index, cas]
    i, t, cas = args
    #state=[alarm,target_rate,dh_min,dh_max]
    push!(log, "CAS_$i", cas_data(t, cas))
end

function cas_names(cas::Union{ACASX_CCAS,ACASX_ADD})
    vcat(["t", "ra_active"], ["ownInput.dz", "ownInput.z", "ownInput.psi", "ownInput.h",
        "ownInput.modes"],
         [["intruderInput[$i].valid", "intruderInput[$i].id", "intruderInput[$i].modes",
         "intruderInput[$i].sr", "intruderInput[$i].chi", "intruderInput[$i].z",
         "intruderInput[$i].cvc", "intruderInput[$i].vrc", "intruderInput[$i].vsb",
         "intruderInput[$i].equipage", "intruderInput[$i].quant",
         "intruderInput[$i].sensitivity_index", "intruderInput[$i].protection_mode"]
         for i = 1:cas.max_intruders]...,
         ["ownOutput.cc", "ownOutput.vc", "ownOutput.ua", "ownOutput.da", "ownOutput.target_rate",
          "ownOutput.turn_off_aurals", "ownOutput.crossing", "ownOutput.alarm", "ownOutput.alert",
          "ownOutput.dh_min", "ownOutput.dh_max", "ownOutput.sensitivity_index", "ownOutput.ddh"],
         [["intruderOutput[$i].id", "intruderOutput[$i].cvc", "intruderOutput[$i].vrc",
          "intruderOutput[$i].vsb", "intruderOutput[$i].tds", "intruderOutput[$i].code"]
          for i = 1:cas.max_intruders]...)
end
function cas_types(cas::Union{ACASX_CCAS,ACASX_ADD})
    vcat([Int64, Bool], [Float64, Float64, Float64, Float64, Int64],
         [[Bool, Int64, Int64, Float64, Float64, Float64, Int64, Int64, Int64,
         Int64, Int64,  Int64, Int64] for i = 1:cas.max_intruders]...,
         [Int64, Int64, Int64, Int64, Float64, Bool, Bool, Bool, Bool,
         Float64, Float64, Int64, Float64], [[Int64, Int64, Int64,
          Int64, Float64, Int64] for i = 1:cas.max_intruders]...)
end
function cas_units(cas::Union{ACASX_CCAS,ACASX_ADD})
    vcat(["index", "bool"], ["ft/s", "ft", "rad", "ft/s", "num"],
        [["bool", "num", "num", "ft", "rad", "ft", "n/a", "n/a", "n/a", "enum",
           "ft", "num", "num"] for i = 1:cas.max_intruders]...,
         ["n/a", "n/a", "n/a", "n/a", "ft/s", "bool", "bool", "bool", "boobool",
          "ft/s", "ft/s", "num", "ft/s^2"],
         [["num", "n/a", "n/a", "n/a", "n/a", "n/a"]
          for i = 1:cas.max_intruders]...)
end
function cas_data(t::Int64, cas::Union{ACASX_CCAS,ACASX_ADD}) #log everything
    ra_active = (cas.output.dh_min > -9999.0 || cas.output.dh_max < 9999.0)::Bool
    ownInput = Any[cas.input.ownInput.dz,
             cas.input.ownInput.z,
             cas.input.ownInput.psi,
             cas.input.ownInput.h,
             cas.input.ownInput.modes]

    intruderInputs = vcat([Any[cas.input.intruders[i].valid,
                             cas.input.intruders[i].id,
                             cas.input.intruders[i].modes,
                             cas.input.intruders[i].sr,
                             cas.input.intruders[i].chi,
                             cas.input.intruders[i].z,
                             cas.input.intruders[i].cvc,
                             cas.input.intruders[i].vrc,
                             cas.input.intruders[i].vsb,
                             cas.input.intruders[i].equipage.val,
                             cas.input.intruders[i].quant,
                             cas.input.intruders[i].sensitivity_index,
                             cas.input.intruders[i].protection_mode]
                          for i=1:length(cas.input.intruders)]...)

    ownOutput = Any[cas.output.cc,
                   cas.output.vc,
                   cas.output.ua,
                   cas.output.da,
                   cas.output.target_rate,
                   cas.output.turn_off_aurals,
                   cas.output.crossing,
                   cas.output.alarm,
                   cas.output.alert,
                   cas.output.dh_min,
                   cas.output.dh_max,
                   cas.output.sensitivity_index,
                   cas.output.ddh]

    intruderOutputs = vcat([Any[cas.output.intruders[i].id,
                              cas.output.intruders[i].cvc,
                              cas.output.intruders[i].vrc,
                              cas.output.intruders[i].vsb,
                              cas.output.intruders[i].tds,
                              cas.output.intruders[i].code]
                          for i=1:length(cas.output.intruders)]...)

    round_floats(vcat(t, ra_active,
                           ownInput,
                           intruderInputs,
                           ownOutput,
                           intruderOutputs), ROUND_NDECIMALS, enable = ENABLE_ROUNDING)
end

function log_response!(log::TaggedDFLogger, args)
    #[aircraft_number, time_index, cas]
    i, t, pr = args
    push!(log, "Response_$i", response_data(t, pr))
end

response_names(pr::StochasticLinearPR) = String["t", "state", "displayRA", "response",
    "v_d", "h_d", "psi_d", "logProb"]
response_types(pr::StochasticLinearPR) = [Int64, String, String, String, Float64, Float64,
    Float64, Float64]
response_units(pr::StochasticLinearPR) = String["enum", "enum", "enum", "ft/s^2", "ft/s",
    "deg/s", "float"]

function response_data(t::Int64, pr::StochasticLinearPR)
    round_floats(Any[t, pr.state,
                          pr.displayRA,
                          pr.response,
                          pr.output.t,
                          pr.output.v_d,
                          pr.output.h_d,
                          pr.output.psi_d,
                          pr.output.logProb], ROUND_NDECIMALS, enable = ENABLE_ROUNDING)
end

response_names(pr::LLDetPR) = String["t", "state", "timer", "t_s", "v_d", "h_d", "psi_d", "logProb"]
response_types(pr::LLDetPR) = [Int64, Symbol, Int64, Float64, Float64, Float64, Float64, Float64]
response_units(pr::LLDetPR) = String["int", "enum", "s", "s", "ft/s^2", "ft/s", "deg/s", "float"]

function response_data(t::Int64, pr::LLDetPR)
    round_floats(Any[t, pr.state,
                          pr.timer,
                          pr.output.t,
                          pr.output.v_d,
                          pr.output.h_d,
                          pr.output.psi_d,
                          pr.output.logProb], ROUND_NDECIMALS, enable = ENABLE_ROUNDING)
end

function log_run_info!(log::TaggedDFLogger, args)
  # called only once when isEndState is true
  reward, md_time, hmd, vmd, nmac = args
  push!(log, "run_info", Any[reward, md_time, hmd, vmd, nmac])
end

run_info_names() = String["reward", "md_time", "hmd", "vmd", "nmac"]
run_info_types() = [Float64, Float64, Float64, Float64, Bool]
run_info_units() = String["float", "s", "ft", "ft", "bool"]

function log_logProb!(log::TaggedDFLogger, args)
    #[time_index, logProb]
    t, logProb = args
    push!(log, "logProb", round_floats(Any[t, logProb], ROUND_NDECIMALS,
        enable=ENABLE_ROUNDING))
end

logprob_names() = String["t", "logprob"]
logprob_types() = [Int64, Float64]
logprob_units() = String["int", "n/a"]

function round_floats(v::Vector{Any}, ndigits::Int64; enable::Bool=true)
  out = v
  if enable
    out = map(v) do x
      return isa(x, AbstractFloat) ? round(x, ndigits) : x
    end
  end
  return out
end

#mods x to the range [-b, b]
function to_plusminus_b(x::AbstractFloat, b::AbstractFloat)
  z = mod(x, 2 * b)
  return (z > b) ? (z - 2 * b) : z
end

to_plusminus_180(x::AbstractFloat) = to_plusminus_b(x, 180.0)

end #module

max_stars_repo_name: sisl/RLESCAS

max_issues_repo_name: sisl/RLESCAS

max_forks_repo_name: sisl/RLESCAS