text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# 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).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ ac87ced2-986f-4e6c-80b2-104b25c171c2
begin
	using CSV, JLD2, FITSIO, FileIO# , HDF5
	using DataFrames #, Query
	using PyCall
	using PlutoUI, PlutoTeachingTools
	eval(Meta.parse(code_for_check_type_funcs))
end

# ╔═╡ e0e60e09-3808-4d6b-a773-6ba59c02f517
md"""
# Astro 528, Lab 3, Exercise 2
# Benchmarking File I/O &
# (and calling Python packages)
"""

# ╔═╡ 8c1d81ab-d215-4972-afeb-7e00bf6063c2
md"""
For many applications, its important that we be able to read input data from a file and/or to write our outputes to files so they can be reused later.  Disk access is typically much slower than accessing system memory.  Therefore, disk access can easily become the limiting factor for a project.  In this set of exercises, you'll see examples of how to perform basic file I/O.  In the process, you'll also see how to import Python packages, call python functions and access data computed by python functions.

You'll be provided with most of the code you need, so that you can focus on comparing how much disk space and time is required by different file formats.  Near the end of the lab, you'll be asked to think about when each type of file format would be a good choice for you to use in your research projects.
"""

# ╔═╡ 1607eac9-e76f-4d1f-a9ce-981ce3be9bea
md"""
### Download some data
First, we're going to download some data from the web.  Julia has a built in `download` function that can be handy for this.  It relies on your system having some utilities already installed (e.g., `curl`, `wget` or `fetch`).  If you run this on a local system and run into trouble, then you can leave the cell below, and manually download the file to the data subdirectory.
"""

# ╔═╡ f27e1e8f-15eb-4754-a94c-7f37c54b871e
begin
	path = basename(pwd())=="test" ? "../data/" : "data/"
	if !isdir(path)  mkdir(path)  end  # make sure there's a data directory
	url = "https://exoplanetarchive.ipac.caltech.edu/data/KeplerData/Simulated/kplr_dr25_inj1_plti.txt"
	filename_ipac = joinpath(path,basename(url)) # extract the filename and prepend "data/"
	time_to_download = NaN
	if !isfile(filename_ipac)  # skip downloading if file already exists
	    time_to_download = @elapsed download(url,filename_ipac)
	end
end

# ╔═╡ 80f02c3a-6751-48df-92ec-13f5c7d8c71e
if !isnan(time_to_download)
	md"Downloading the file took $time_to_download seconds."
end

# ╔═╡ 624c9038-3008-4e78-a149-60796dacf9c0
md"""
Previously, everything you needed for an assignment was included in a GitHub repository.  So why did I make you download the file?

Notice the size of the file.  Git is great for tracking source code, but it wasn't really designed for working with large files (especially large _binary_ files).  Since  we're not going to be editing it, we'll simply download it once.  Besides, it's useful to know how to download a file from within a julia script and to compare the time required to download the file from the internet to the time required to read the file from disk.
"""

# ╔═╡ 2eb255d9-707d-4224-a0ce-fe90a1c69722
md"""
## Calling AstroPy to read data in unusual formats

Here I've picked a file containing the results of applying the pipeline for NASA's Kepler mission to [simulated data](https://exoplanetarchive.ipac.caltech.edu/docs/KeplerSimulated.html) in which the signals of simulated "planet's" have been injected into actual Kepler data.  This data set is the basis for computing the efficiency of the Kepler pipeline at detecting real planets.  This dataset has played an important role in enabling astronomers to estimate the occurrence rates of planets around other stars.
For documentation of its contents, you could read [its documentation](https://exoplanetarchive.ipac.caltech.edu/docs/KSCI-19110-001.pdf).  However, that's not necessary for this lab.  For now, we'll just do some basic manipulations of the file, and the details of its contents aren't important.  That said, it is important to know the _file format_.

This data file that we downloaded is in [IPAC format](https://irsa.ipac.caltech.edu/applications/DDGEN/Doc/ipac_tbl.html).
It would be tedious to learn the details of every file format that is used in astronomy, let alone to write our own code to read them.  Fortunately, there are packages that can read the most common file formats.  Sometimes astronomers use specialized/less common file formats which aren't always implemented in native Julia.  Reading such files is an example of something that [astropy](http://docs.astropy.org) is particularly good for.  It provides a function [`astropy.io.ascii.read`](http://docs.astropy.org/en/stable/io/ascii/) that will read a file in IPAC for us.
Since astropy is written in Python, we import the [`PyCall`](https://github.com/stevengj/PyCall.jl) package, so we can import python packages and call python functions from Julia.
Since reading a file from disk is typically limited by the rate of getting data from disk, rather than compute speed, it's usually not a problem that Python isn't particularly fast, when it's comes to reading files.
"""

# ╔═╡ 8571ccaf-5fbc-4593-82db-ead073a4074f
md"""
At the bottom of this notebook, I've run `using PyCall`.  Here, we'll import the required python module using `pyimport`.
"""

# ╔═╡ 442604d0-b490-404b-8b1c-7e1c01cebad7
astropy_io_ascii = pyimport("astropy.io.ascii")

# ╔═╡ 70f5fe63-fea1-4afc-b03d-f113b8bd621e
md"""
Now, we can review the documentation for [`astropy.io.ascii.read`](http://docs.astropy.org/en/stable/io/ascii/) and call that function load the data in our input file.  Let's use `@time`, so we can compare the time required to read various formats.
"""

# ╔═╡ 8cd5a772-a3aa-4166-92e6-39436d0d2278
with_terminal() do
	@time astropy_io_ascii.read(filename_ipac, format="ipac", fast_reader=false)
end

# ╔═╡ 074b00ce-dae6-4e4f-8673-49be73085327
md"""
On the plus side, the data was read in.  However, what is the type of the data that astropy read for us and stored into `data_from_astropy`?  Let's check.
"""

# ╔═╡ ced3ecc9-94ea-408f-bdef-d34802166663
begin
	time_to_read_with_astropy = @elapsed data_from_astropy = astropy_io_ascii.read(filename_ipac, format="ipac", fast_reader=false)
	typeof(data_from_astropy)
end

# ╔═╡ c62adfb6-6ef8-45fd-992f-07bca59f82cd
md"""
Since Python is a weakly-typed language, it's type is `PyObject`.  That can contain most anything!
That flexibility can be convenient, but it is also one of the reasons that Python is not a good language for high-performance computing.  To enable Julia to work efficiently with the data, we'll want Julia to know what type the data is and store a list of strictly-typed columns into a the data into a `DataFrame`.
The PyCall package provides an interface for For accessing data from PyObjects.  Because of the weak typing issue, the syntax and details can be a bit confusing.  If you're not intending to call Python for your class project, then there's no reason to worry about that these details.
"""

# ╔═╡ 03edae54-4368-4b81-8713-17f93bbb9ed9
protip(md"""
This list is just for students who are curious about accessing Python data and methods from Julia.

   - Given o::PyObject, o.attribute in Julia is equivalent to o.attribute in Python, with automatic type conversion. To get an attribute as a PyObject without type conversion, do o."attribute" instead. The keys(o::PyObject) function returns an array of the available attribute symbols.

   - Given o::PyObject, get(o, key) is equivalent to o[key] in Python, with automatic type conversion. To get as a PyObject without type conversion, do get(o, PyObject, key), or more generally get(o, SomeType, key). You can also supply a default value to use if the key is not found by get(o, key, default) or get(o, SomeType, key, default). Similarly, set!(o, key, val) is equivalent to o[key] = val in Python, and delete!(o, key) is equivalent to del o[key] in Python. For one or more integer indices, o[i] in Julia is equivalent to o[i-1] in Python.

   - You can call an o::PyObject via o(args...) just like in Python (assuming that the object is callable in Python). The explicit pycall form is still useful in Julia if you want to specify the return type.

   - pystr(o) and pyrepr(o) are analogous to str and repr in Python, respectively.

   - There's more information about accessing data in PyObjects (and other types to contain Python data) in the [PyCall documentation](https://github.com/JuliaPy/PyCall.jl#types).  (The above info is copied directly from there.)
""")

# ╔═╡ e0b82e6c-043c-4e66-a2a4-ae0f3e4e5404
md"""
For now, you can get a [`Dict` or "Dictionary"](https://docs.julialang.org/en/v1/base/collections/index.html#Dictionaries-1) by using `data_from_astropy.columns`.  The dictionary consists of a set of _keys_ (in this case strings), where each key is associated with a _value_ (in this case a Vector or 1-d array).
Often times, the data you want to work with can be represented as a table.  For efficiency's sake, it's usually best to represent these as a bunch of `Vector`'s, each containing one columns of data.  Using a `Dict` allows you to give the columns names (instead of just numbers, so you're less likely to access the wrong column) and allows each column to have a different type (again useful for Julia to optimize your code).
"""

# ╔═╡ aa62f1e5-6af5-43f6-be11-665ed25986c0
dict_from_astropy = data_from_astropy.columns

# ╔═╡ 540468dc-3f9f-457b-9b44-efb3ec036166
md"""
Most of the columns can be automatically converted to an array with a known Julia type.  For example, `data_from_astropy.columns["KIC_ID"]` returns data as an `Array{Int64,1}`.  However, value associated with `TCE_ID` is some `PyObject` that can't be automatically reinterpretted as a Julia array.
"""

# ╔═╡ 21e3adbb-0176-4963-bdfc-0ba2422af4bf
typeof(dict_from_astropy["KIC_ID"])

# ╔═╡ d27aef3b-bfe7-466a-9580-93fc53dcef18
typeof(dict_from_astropy["TCE_ID"])

# ╔═╡ 02f638d9-ec2f-4072-9fb3-e6fabac1b1e6
md"""Since many `KIC_ID`'s don't have an associated `TCE_ID`, there are many missing entries.  Python is trying to store a list with lots of empty entries efficiently, but PyCall doesn't (yet?) know how to deal with this `MaskedColumn` for us.  Working with the data when Julia can't know its type would be very inefficient.  Therefore, we want to create an array of Strings that allows Julia to represent this data more efficiently.  Technically, it will be an array where each element is either a `String` or a `missing`.  It took a little tinkering, but eventually, I figured out how to extract that data into an efficient Julia object.
"""

# ╔═╡ 481c5b7a-4660-4954-84f3-29d33bd73f3d
begin
	TCE_ID_list = map(x -> x != nothing ? x : missing, data_from_astropy.columns["TCE_ID"].data.tolist())
	TCE_ID_list_type = typeof(TCE_ID_list)
	TCE_ID_list
end

# ╔═╡ aba84742-a927-4ce5-9f7c-6e535db1ee47
protip(md"Note that the type of `TCE_ID_list` is a $TCE_ID_list_type.  Julia is able to store and compute on vectors with missing data very efficiently thanks to its rich type system.")

# ╔═╡ aa7c4dfe-a3d4-448b-a559-1bc7b338a1dc
md"""
Now let's replace the value associated with TCE_ID with this list.
"""

# ╔═╡ a97f9c8e-dbd3-4613-9c6e-c471340ea2d6
dict_from_astropy["TCE_ID"] = TCE_ID_list;

# ╔═╡ 7c7430a0-d16e-485b-bbd9-7f231fee853a
md"""
If we just wantted to access the data, then we could use the data stored in the columns as a dictionary.
However, a dictionary doesn't guarentee anything about the relationship of the value of different keys.
For example, in a table, each column should have the same number of rows.  Therefore, we'll switch from representing the data as a dictionary and start using a `DataFrame`, a type provided by the [DataFrames.jl package](https://github.com/JuliaData/DataFrames.jl).  A `DataFrame` can be thought of as a table, where the data for each column is stored as an array.   A `DataFrame` also provides some additional features to allow easy and efficient access and manipulation of the table that will come in useful later.
"""

# ╔═╡ cebe52b3-386c-4d31-8497-c19b6c742577
begin
	# First, we'll create and copy a small `DataFrame`, so the functions get compiled before we start timing things.
	small_dict = Dict(:a=>[1,2],:b=>[3.0,4.0],:c=>["hello",missing])
	small_df1 = DataFrame(small_dict)
	small_df2 = DataFrame(small_dict, copycols=false)
end;

# ╔═╡ 3d09d38a-ca72-43f2-bd85-17fc0e05a645
md"""
We'll use our existing dictionary to initialize a `DataFrame`.
"""

# ╔═╡ 9257c879-d96f-4d2b-994e-d542617b0c65
with_terminal() do
	small_dict  # makes sure already compiled code for creating small DataFrame
	@time df_tmp = DataFrame(dict_from_astropy)
end

# ╔═╡ 307a7a29-afa6-4c4d-88ae-4f757aeba892
md"""
2a.  Look at how much memory was allocated during this line of code.
Do you think it involved allocating and copying lots of data from the dictionary?

"""

# ╔═╡ afdc442a-2931-4129-824e-98431e1d8be2
response_2a = missing

# ╔═╡ c42deea3-acf1-491a-97bd-fdc472917584
display_msg_if_fail(check_type_isa(:response_2a,response_2a,[AbstractString,Markdown.MD]))

# ╔═╡ 7c1f4a27-dfd7-428a-97ef-eac1784e8d6f
md"""By default `DataFrames` prioritizes _safe operations_ over efficient ones.  However, it also provides options for more efficient operations for times when users want to invest a little more time to get good performance.  For example, we can tell it not to copy the data, using `copycols=false`)."""

# ╔═╡ 97111ccc-8a16-4e7d-adf7-0dcd90524be3
with_terminal() do
		@time df = DataFrame(dict_from_astropy, copycols=false)
end

# ╔═╡ f45668b1-1cbd-48e6-9224-d469826e093b
md"""
2b.  Look at how much memory was allocated the second time.
Did Julia make a new copy of all of the data?  How did the time required compare to the first attempt?
"""

# ╔═╡ 52beb148-6214-4a03-a0cf-6c2a6da77d40
response_2b = missing

# ╔═╡ 1ecbd400-0157-45ce-b1b8-db764012ba5e
display_msg_if_fail(check_type_isa(:response_2b,response_2b,[AbstractString,Markdown.MD]))

# ╔═╡ 191ba96e-2573-4bc1-a352-46a66e0a5c4f
md"""
## Writing a CSV file

The IPAC format allows for significant metadata, but reading not-so-common file formats can be annoying.
Let's say that we'd like to write the data to a [CSV file](https://en.wikipedia.org/wiki/Comma-separated_values), so that it's easier for other programs to read in.  Using the CSV package, we can read a CSV file into a DataFrame (or several other tabular data structures) and write CSV files from a DataFrame with code like the following.
"""

# ╔═╡ 82a757ad-566d-4c1d-8b3d-366ffd980fb4
begin
	# Read & write a small test file so compilation time not included below
	n_small_df = 1024
	small_csv_filename = "random_numbers.csv"
	small_df = DataFrame(a=rand(n_small_df),b=rand(n_small_df) )
	CSV.write(joinpath(path,small_csv_filename),small_df)
	small_df_from_csv = CSV.read(joinpath(path,small_csv_filename),DataFrame)
end;

# ╔═╡ 3e550b71-4750-460b-be18-911a848a8f49
md"""
Now, let's write the DataFrame from our IPAC file to a CSV file.  We'll investigate how the filesize and time to read files compares.
"""

# ╔═╡ b0c3875b-ef07-4e92-a0a8-55f42b266c6b
begin
	filename_csv = replace(filename_ipac, ".txt" => ".csv")
	df = DataFrame(dict_from_astropy, copycols=false)
	with_terminal() do
		@time CSV.write(filename_csv,df)
	end
end

# ╔═╡ 58e021c6-500f-40d8-a388-a5732bc808b3
md"""Next, we'll compare the filesizes for the IPAC file and the CSV file containing the same data.

2c.  Which do you expect will be larger?  Why?
Once you've made a prediction, mouse over the hint box below to see the sizes of the two files.  If you were suprsised, try to explain what happened."""

# ╔═╡ 36bfd224-a13b-438c-9106-382c15a6a1d2
response_2c = missing

# ╔═╡ 7f4c3ab8-041e-4ef0-bb2e-320f5293cdde
display_msg_if_fail(check_type_isa(:response_2c,response_2c,[AbstractString,Markdown.MD]))

# ╔═╡ 195c2df4-7e1c-49f8-871f-f76b985dab6d
begin
	ipac_filesize = filesize(filename_ipac) /1024^2
	csv_filesize = filesize(filename_csv) /1024^2
	hint(md"The IPAC file size is $ipac_filesize MB versus $csv_filesize MB for the CSV.")
end

# ╔═╡ 32124f8a-f5cf-41d1-97a9-ce1d05145fde
md"""

## Reading CSV files

2d. Think back to how long it took to read in the file in IPAC format ($time_to_read_with_astropy seconds).  How long do you think it will take to read in the same data once it's been stored in CSV format?
"""

# ╔═╡ 507e1516-5433-49eb-831d-32426f30895e
response_2d = missing

# ╔═╡ eac67cc9-754b-4f7d-add8-93900a1b5b49
display_msg_if_fail(check_type_isa(:response_2d,response_2d,[AbstractString,Markdown.MD]))

# ╔═╡ 64224c6b-c5a0-44f2-b2a0-7f77759cb848
md"""
Now, try reading it in and see.
"""

# ╔═╡ 73c06bb1-be49-46bd-b7f1-c45cc56af7b4
md"""
I'm ready to benchmark reading a CSV file. $(@bind ready_read_csv CheckBox())
$(@bind go_read_csv Button("Rerun the benchmarks."))
"""

# ╔═╡ ebfaa677-829b-4bf9-bdbb-19f3c87dd3a4
if  ready_read_csv
	go_read_csv
	time_read_csv = @elapsed df_csv = CSV.read(filename_csv, DataFrame)
	md"It took $time_read_csv seconds to read the CSV file."
end

# ╔═╡ 945f5a55-3026-4497-9ece-8af878c87788
md"""
2e.  How did the time required to read the data in CSV format compare to the time to read the data in IPAC format?  If you were suprsised, try to explain what happened.
"""

# ╔═╡ 57397ee4-9efc-48b3-b640-d2b7a10da633
response_2e = missing

# ╔═╡ 8059a6a3-384a-4344-8a23-650ee0be10c2
display_msg_if_fail(check_type_isa(:response_2e,response_2e,[AbstractString,Markdown.MD]))

# ╔═╡ f5b93929-2c59-4360-8c41-97a1324ba455
md"2f. How do you think the sizes of the files in the two formats will compare?  Once you've made your prediction, mouse over the hint box.  If you were suprsised, try to explain what happened."

# ╔═╡ 122196fa-45ca-4031-85eb-4afd4782de9e
response_2f = missing

# ╔═╡ e9dc1456-616b-4e4b-b209-9f6ba4c48607
display_msg_if_fail(check_type_isa(:response_2f,response_2f,[AbstractString,Markdown.MD]))

# ╔═╡ f15d37a7-d962-4da0-977f-76729a3313be
hint(md"The IPAC file size is $ipac_filesize MB versus $csv_filesize MB for the CSV.")

# ╔═╡ 28f195c4-4f61-4873-85d6-b4e3aaa3660f
md"""
## Binary formats: HDF5/JLD2

There are numerous binary file formats that one could use.  Here, we'll try using JLD2 which is a subset of the [HDF5](https://www.hdfgroup.org/solutions/hdf5/) file format.  This means that when [Julia's JLD2 package](https://github.com/JuliaIO/JLD2.jl) writes jld2 files, they can be read by other programs that can read HDF5 files.  However, a generic HDF5 file is not a valid JLD2 file.  If you want to read a HDF5 file, then you can use Julia's [HDF5.jl package](https://github.com/JuliaIO/HDF5.jl).  The [FileIO.jl](https://github.com/JuliaIO/FileIO.jl) package provides a common interface for reading and writing from multiple file formats, including these and several others.

As before, we'll call each function once using a small DataFrame, just so they get compiled before we benchmark them.
"""

# ╔═╡ 3837e439-250b-4577-b0d7-93352ec19f6e
begin
	filename_jld2_small = replace(small_csv_filename, ".csv" => ".jld2")
	@save joinpath(path,filename_jld2_small) small_df
	small_df_from_jld2 = load(joinpath(path,filename_jld2_small), "small_df")
end;

# ╔═╡ df9b701f-d314-4512-b2ea-1f6ae015166c
response_2g = missing

# ╔═╡ 377f7527-a338-4526-bb24-9766c635e719
display_msg_if_fail(check_type_isa(:response_2g,response_2g,[AbstractString,Markdown.MD]))

# ╔═╡ 25f24754-16d4-4343-b54a-cf8ea1358ce9
begin
	small_jld2_filesize = filesize(joinpath(path,filename_jld2_small))/1024
	small_csv_filesize = filesize(joinpath(path,small_csv_filename))/1024
	hint(md"For the small random data set, the JLD2 file's size is $small_jld2_filesize KB versus $small_csv_filesize KB for the CSV file.")
end

# ╔═╡ 8255b5a9-cbe6-4604-bedd-0e366f311096
md"""
Let's check the filesizes for the JLD2 file to the CSV file.
The small CSV file we created is $small_jld2_filesize KB.

2g.  How large would you guess the JLD2 file will be?  Once you've made your prediction, mouse over the hint box.  If you were suprised, try to explain what happened.
"""

# ╔═╡ 570fd826-23fd-46ee-bdb4-58fb0c45719a
md"""
Now let's time how long it takes to save the data from IPAC into a JLD2 file.
"""

# ╔═╡ 691410bb-0472-4800-a90d-29ddf947de3e
begin
	filename_jld2 = replace(filename_ipac, ".txt" => ".jld2")
	with_terminal() do
		@time @save filename_jld2 df
	end
end

# ╔═╡ 8cbb1c90-bd94-44b5-80b6-81d38f3e6252
md"2h.  How long do you think it will take to load the data from the JLD2 file? "

# ╔═╡ c3065acf-6205-455f-ba74-ca51f3f6761b
response_2h = missing

# ╔═╡ 9c392be9-1505-40bc-a290-68085a1c2700
display_msg_if_fail(check_type_isa(:response_2h,response_2h,[AbstractString,Markdown.MD]))

# ╔═╡ d738bdcc-5f83-4dfa-a17f-e9ea23db2986
md"""
Now time how long it takes to load the data from the JLD2 file.
"""

# ╔═╡ d55ce157-099c-4c1b-94db-62918f04e5fe
md"""
I'm ready to benchmark reading a JLD2 file. $(@bind ready_read_jld2 CheckBox())
$(@bind go_read_jld2 Button("Rerun the benchmarks."))
"""

# ╔═╡ 206f464b-55fe-46aa-85b7-8f0246a0aaad
if ready_read_jld2
	go_read_jld2
	time_read_jld2 = @elapsed df_jld2_read_back_in = load(filename_jld2, "df")
	md"It took $time_read_jld2 seconds to read the JLD2 file."
end

# ╔═╡ 6f72d1b2-63f6-4301-8272-bb2d6d2d049e
md"""
2i.  How does the time required to read and write the JLD2 file compare to the time required to read the IPAC and CSV formatted files?  If you were suprised, try to explain what happend.
"""

# ╔═╡ 01201c37-0b79-46b1-a001-e716f5b3ba67
response_2i = missing

# ╔═╡ e3d37bc1-a119-4add-a111-899ee0caea05
display_msg_if_fail(check_type_isa(:response_2i,response_2i,[AbstractString,Markdown.MD]))

# ╔═╡ 6dff3f21-fac0-42e1-910a-f969a231374f
md"""
Next, we'll compare the filesizes for the JLD2 file to the CSV file.

2j.  How large would you guess the JLD2 file will be?  Once you've made your prediction, mouse over the hint box.  If you were suprised, try to explain what happened.
"""

# ╔═╡ b26b8253-e6cd-49f4-81c5-2a3c2963a37c
response_2j = missing

# ╔═╡ e5dc123f-1596-4311-9398-f0cfe80a5342
display_msg_if_fail(check_type_isa(:response_2j,response_2j,[AbstractString,Markdown.MD]))

# ╔═╡ 7b9d26c2-899e-45e9-b664-39d5f1adfe3f
begin
	jld2_filesize = filesize(filename_jld2) / 1024^2;
	hint(md"The JLD2 file size is $jld2_filesize MB versus $csv_filesize MB for the CSV.")
end

# ╔═╡ fc01d57f-c90b-4231-96be-ddd48656d55e
md"""
## Flexible Format: FITS

Astronomers often use the [FITS file format](https://en.wikipedia.org/wiki/FITS).  Like [HDF
5](https://www.hdfgroup.org/solutions/hdf5/), it's a very flexible (e.g., it can store both text and binary data) and thus complicated file
 format.
Therefore, most languages call a common [FITSIO library written in C](https://heasarc.gsfc.nasa.gov/fitsio/), rather than implementing code themselves.  Indeed, that's what [Julia's FITSIO.jl package](https://github.com/JuliaAstro/FITSIO.jl) does.

Unfortunately, the FITSIO package isn't as polished as the others.  It expects a `Dict` rather than a `DataFrame`, and it can't handle missing values.  So I've provided some helper functions at the bottom of the notebook.  Also, FITS files have complicated headers, so I'll provide a function to read all the tabular data from a simple FITS file.  As usual, we'll use each function once, so that Julia compiles them before we start timing.
"""

# ╔═╡ 3b232365-f2fe-4edb-a39f-3e37c8cbb666
md"Now we can time how long it takes to write and read the data as FITS files."

# ╔═╡ 568862b3-6fce-426a-a9e2-e558adf3932a
md"""
I'm ready to benchmark reading a FITS file. $(@bind ready_read_fits CheckBox())
$(@bind go_read_fits Button("Rerun the benchmarks."))
"""

# ╔═╡ ae79ee89-d788-4612-8b1d-fc22d85c7744
md"2k.  How do the read/write times and file sizes for FITS compare to CSV and JLD2?"

# ╔═╡ 3745c237-ba48-4f2c-959e-441484244764
response_2k = missing

# ╔═╡ 6d061411-6f19-4119-aefb-cc380198b0ce
display_msg_if_fail(check_type_isa(:response_2k,response_2k,[AbstractString,Markdown.MD]))

# ╔═╡ 8def87d2-f10b-4a82-b353-a6477eeead9b
md"## Implications for your project"

# ╔═╡ c74b3105-f480-4688-b85f-3e7dff70da3b
md"""
2l.  How does the time required to read any of the above formats compare to the time required to download the files ($time_to_download seconds)?
Will your project need to transfer large files over the internet?  If so, very roughly how large do you expect they will be?
"""

# ╔═╡ 55438c09-1d94-4ff7-90c3-0cc6064a091e
response_2l = missing

# ╔═╡ bcc796c9-db11-4a09-a5f9-215127ac0938
display_msg_if_fail(check_type_isa(:response_2l,response_2l,[AbstractString,Markdown.MD]))

# ╔═╡ 92b2ecdd-9491-4f16-8cb9-bacbcf180280
md"""
2m.  Will your project need to read in any input files?  If so, what format are they in?

I assume everyone's project code will write at least some output files?  Very roughly, how large do you expect that they'll be?   What format(s) would be a good choice for your project?"""

# ╔═╡ 83216915-cdcf-4f0f-829f-a5ff4c4b8da0
response_2m = missing

# ╔═╡ 452872f9-2009-4733-b2d3-28f262ae19b7
display_msg_if_fail(check_type_isa(:response_2m,response_2m,[AbstractString,Markdown.MD]))

# ╔═╡ 29415ddc-e002-4f56-a169-95f7b1c36be9
md"# Helper Functions"

# ╔═╡ fb23d6c6-b812-4fe1-b224-0014bedbd43f
ChooseDisplayMode()

# ╔═╡ 1e53aa10-dff6-40d5-89e2-da194ffc2052
TableOfContents()

# ╔═╡ 14cca8ce-cc61-4fae-b871-21c3fd23d0ea
"Convert a DataFrame to a Dict, replacing missing values with 0 or an empty string."
function convert_dataframe_to_dict_remove_missing(df::DataFrame)
    d = Dict(map(n->"$n"=>               # create a dictionary
              ( any(ismissing.(df[!,n])) ? # if column contains a missing
                    map(x-> !ismissing(x) ? # search for missings
                        x :                 # leave non-missing values alone
                        ( (eltype(df[!,n]) <: Number) ? zero(eltype(n)) : "")
                        , df[!,n])            # but replace missing with 0 or ""
                : df[!,n] ), # if nothing is missing, just use column as is
            names(df) ))
end

# ╔═╡ 28fc8de4-749b-4093-b32f-c398f8d27d3d
"Write a DataFrame to a FITS file, replacing missing values with 0 or an empty string."
function write_dataframe_as_fits(filename::String, df::DataFrame)
    try
       dict = convert_dataframe_to_dict_remove_missing(df)
       fits_file = FITS(filename,"w")
       write(fits_file, dict )
       close(fits_file)
    catch
        @warn("There was a problem writing a dataframe to " * filename * ".")
    end
end

# ╔═╡ 57b422e5-0ad0-4674-bdd3-a8358bc7aaeb

"Read the columns of the first table from a FITS file into a Dict"
function read_fits_tables(filename::String)
    dict = Dict{String,Any}()
    fits_file = FITS(filename,"r")
    # fits_file[1] is image data, we're interested in the table
    @assert length(fits_file) >= 2
    header = read_header(fits_file[2])
    for i in 1:length(header)
        c = get_comment(header,i)
        if !occursin("label for field",c)
            continue
        end
        h = header[i]
        @assert typeof(h) == String
        try
            dict[h] = read(fits_file[2],h)
        catch
            @warn "# Problem reading table column " * h * "."
        end
    end
    close(fits_file)
    return dict
end

# ╔═╡ e05f16d6-eb50-49a9-bf14-95d63c9da7ff
begin
	filename_fits_small = replace(small_csv_filename, ".csv" => ".fits")
	write_dataframe_as_fits(joinpath(path,filename_fits_small),small_df)
	read_fits_tables(joinpath(path,filename_fits_small))
end;

# ╔═╡ 1992595f-9976-4b18-bf9c-df8a73d30dc8
begin
	filename_fits_small # make sure have already compiled functions
	filename_fits = replace(filename_ipac, ".txt" => ".fits")
	with_terminal() do
		@time write_dataframe_as_fits(filename_fits,df)
	end
end

# ╔═╡ 90d5244e-17be-4601-b922-8c254f1248bf
begin
	fits_filesize = filesize(filename_fits) /1024^2
	hint(md"The FITS file size is $fits_filesize MB versus $jld2_filesize MB for the JLD2 and $csv_filesize MB for the CSV.")
end

# ╔═╡ 52f9edd0-79b0-4a9a-9930-3a05d3aa2447
if ready_read_fits
	go_read_fits
	time_read_fits = @elapsed read_fits_tables(filename_fits)
	md"It took $time_read_fits seconds to read the FITS file."
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
FITSIO = "525bcba6-941b-5504-bd06-fd0dc1a4d2eb"
FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
PlutoTeachingTools = "661c6b06-c737-4d37-b85c-46df65de6f69"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"

[compat]
CSV = "~0.8.5"
DataFrames = "~1.2.2"
FITSIO = "~0.16.7"
FileIO = "~1.10.1"
JLD2 = "~0.4.13"
PlutoTeachingTools = "~0.1.3"
PlutoUI = "~0.7.9"
PyCall = "~1.92.3"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[CFITSIO]]
deps = ["CFITSIO_jll"]
git-tree-sha1 = "c860f5545064216f86aa3365ec186ce7ced6a935"
uuid = "3b1b4be9-1499-4b22-8d78-7db3344d1961"
version = "1.3.0"

[[CFITSIO_jll]]
deps = ["Artifacts", "JLLWrappers", "LibCURL_jll", "Libdl", "Pkg"]
git-tree-sha1 = "2fabb5fc48d185d104ca7ed7444b475705993447"
uuid = "b3e40c51-02ae-5482-8a39-3ace5868dcf4"
version = "3.49.1+0"

[[CSV]]
deps = ["Dates", "Mmap", "Parsers", "PooledArrays", "SentinelArrays", "Tables", "Unicode"]
git-tree-sha1 = "b83aa3f513be680454437a0eee21001607e5d983"
uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
version = "0.8.5"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "727e463cfebd0c7b999bbf3e9e7e16f254b94193"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.34.0"

[[Conda]]
deps = ["JSON", "VersionParsing"]
git-tree-sha1 = "299304989a5e6473d985212c28928899c74e9421"
uuid = "8f4d0f93-b110-5947-807f-2305c1781a2d"
version = "1.5.2"

[[Crayons]]
git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.0.4"

[[DataAPI]]
git-tree-sha1 = "ee400abb2298bd13bfc3df1c412ed228061a2385"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.7.0"

[[DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "d785f42445b63fc86caa08bb9a9351008be9b765"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.2.2"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.10"

[[DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[FITSIO]]
deps = ["CFITSIO", "Printf", "Reexport", "Tables"]
git-tree-sha1 = "85b66005c9d16d3d27c9d06fd3be5e25074128da"
uuid = "525bcba6-941b-5504-bd06-fd0dc1a4d2eb"
version = "0.16.7"

[[FileIO]]
deps = ["Pkg", "Requires", "UUIDs"]
git-tree-sha1 = "256d8e6188f3f1ebfa1a5d17e072a0efafa8c5bf"
uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
version = "1.10.1"

[[Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[InvertedIndices]]
deps = ["Test"]
git-tree-sha1 = "15732c475062348b0165684ffe28e85ea8396afc"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.0.0"

[[IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[JLD2]]
deps = ["DataStructures", "FileIO", "MacroTools", "Mmap", "Pkg", "Printf", "Reexport", "TranscodingStreams", "UUIDs"]
git-tree-sha1 = "59ee430ac5dc87bc3eec833cc2a37853425750b4"
uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
version = "0.4.13"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[LaTeXStrings]]
git-tree-sha1 = "c7f1c695e06c01b95a67f0cd1d34994f3e7db104"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.2.1"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "0fb723cd8c45858c22169b2e42269e53271a6df7"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.7"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "2ca267b08821e86c5ef4376cffed98a46c2cb205"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.1"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "bfd7d8c7fd87f04543810d9cbd3995972236ba1b"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "1.1.2"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PlutoTeachingTools]]
deps = ["LaTeXStrings", "Markdown", "PlutoUI", "Random"]
git-tree-sha1 = "e2b63ee022e0b20f43fcd15cda3a9047f449e3b4"
uuid = "661c6b06-c737-4d37-b85c-46df65de6f69"
version = "0.1.4"

[[PlutoUI]]
deps = ["Base64", "Dates", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "Suppressor"]
git-tree-sha1 = "44e225d5837e2a2345e69a1d1e01ac2443ff9fcb"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.9"

[[PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "cde4ce9d6f33219465b55162811d8de8139c0414"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.2.1"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "0d1245a357cc61c8cd61934c07447aa569ff22e6"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.1.0"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[PyCall]]
deps = ["Conda", "Dates", "Libdl", "LinearAlgebra", "MacroTools", "Serialization", "VersionParsing"]
git-tree-sha1 = "169bb8ea6b1b143c5cf57df6d34d022a7b60c6db"
uuid = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
version = "1.92.3"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "5f6c21241f0f655da3952fd60aa18477cf96c220"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.1.0"

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[SentinelArrays]]
deps = ["Dates", "Random"]
git-tree-sha1 = "54f37736d8934a12a200edea2f9206b03bdf3159"
uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
version = "1.3.7"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[Suppressor]]
git-tree-sha1 = "a819d77f31f83e5792a76081eee1ea6342ab8787"
uuid = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
version = "0.2.0"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "d0c690d37c73aeb5ca063056283fde5585a41710"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.5.0"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[TranscodingStreams]]
deps = ["Random", "Test"]
git-tree-sha1 = "216b95ea110b5972db65aa90f88d8d89dcb8851c"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.6"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[VersionParsing]]
git-tree-sha1 = "80229be1f670524750d905f8fc8148e5a8c4537f"
uuid = "81def892-9a0e-5fdd-b105-ffc91e053289"
version = "1.2.0"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─e0e60e09-3808-4d6b-a773-6ba59c02f517
# ╟─8c1d81ab-d215-4972-afeb-7e00bf6063c2
# ╟─1607eac9-e76f-4d1f-a9ce-981ce3be9bea
# ╠═f27e1e8f-15eb-4754-a94c-7f37c54b871e
# ╠═80f02c3a-6751-48df-92ec-13f5c7d8c71e
# ╟─624c9038-3008-4e78-a149-60796dacf9c0
# ╟─2eb255d9-707d-4224-a0ce-fe90a1c69722
# ╟─8571ccaf-5fbc-4593-82db-ead073a4074f
# ╠═442604d0-b490-404b-8b1c-7e1c01cebad7
# ╟─70f5fe63-fea1-4afc-b03d-f113b8bd621e
# ╠═8cd5a772-a3aa-4166-92e6-39436d0d2278
# ╟─074b00ce-dae6-4e4f-8673-49be73085327
# ╠═ced3ecc9-94ea-408f-bdef-d34802166663
# ╟─c62adfb6-6ef8-45fd-992f-07bca59f82cd
# ╟─03edae54-4368-4b81-8713-17f93bbb9ed9
# ╟─e0b82e6c-043c-4e66-a2a4-ae0f3e4e5404
# ╠═aa62f1e5-6af5-43f6-be11-665ed25986c0
# ╟─540468dc-3f9f-457b-9b44-efb3ec036166
# ╠═21e3adbb-0176-4963-bdfc-0ba2422af4bf
# ╠═d27aef3b-bfe7-466a-9580-93fc53dcef18
# ╟─02f638d9-ec2f-4072-9fb3-e6fabac1b1e6
# ╠═481c5b7a-4660-4954-84f3-29d33bd73f3d
# ╟─aba84742-a927-4ce5-9f7c-6e535db1ee47
# ╟─aa7c4dfe-a3d4-448b-a559-1bc7b338a1dc
# ╠═a97f9c8e-dbd3-4613-9c6e-c471340ea2d6
# ╟─7c7430a0-d16e-485b-bbd9-7f231fee853a
# ╟─cebe52b3-386c-4d31-8497-c19b6c742577
# ╟─3d09d38a-ca72-43f2-bd85-17fc0e05a645
# ╠═9257c879-d96f-4d2b-994e-d542617b0c65
# ╟─307a7a29-afa6-4c4d-88ae-4f757aeba892
# ╠═afdc442a-2931-4129-824e-98431e1d8be2
# ╟─c42deea3-acf1-491a-97bd-fdc472917584
# ╟─7c1f4a27-dfd7-428a-97ef-eac1784e8d6f
# ╠═97111ccc-8a16-4e7d-adf7-0dcd90524be3
# ╟─f45668b1-1cbd-48e6-9224-d469826e093b
# ╠═52beb148-6214-4a03-a0cf-6c2a6da77d40
# ╟─1ecbd400-0157-45ce-b1b8-db764012ba5e
# ╟─191ba96e-2573-4bc1-a352-46a66e0a5c4f
# ╠═82a757ad-566d-4c1d-8b3d-366ffd980fb4
# ╟─3e550b71-4750-460b-be18-911a848a8f49
# ╠═b0c3875b-ef07-4e92-a0a8-55f42b266c6b
# ╟─58e021c6-500f-40d8-a388-a5732bc808b3
# ╠═36bfd224-a13b-438c-9106-382c15a6a1d2
# ╟─7f4c3ab8-041e-4ef0-bb2e-320f5293cdde
# ╟─195c2df4-7e1c-49f8-871f-f76b985dab6d
# ╟─32124f8a-f5cf-41d1-97a9-ce1d05145fde
# ╠═507e1516-5433-49eb-831d-32426f30895e
# ╟─eac67cc9-754b-4f7d-add8-93900a1b5b49
# ╟─64224c6b-c5a0-44f2-b2a0-7f77759cb848
# ╟─73c06bb1-be49-46bd-b7f1-c45cc56af7b4
# ╟─ebfaa677-829b-4bf9-bdbb-19f3c87dd3a4
# ╟─945f5a55-3026-4497-9ece-8af878c87788
# ╠═57397ee4-9efc-48b3-b640-d2b7a10da633
# ╟─8059a6a3-384a-4344-8a23-650ee0be10c2
# ╟─f5b93929-2c59-4360-8c41-97a1324ba455
# ╠═122196fa-45ca-4031-85eb-4afd4782de9e
# ╟─e9dc1456-616b-4e4b-b209-9f6ba4c48607
# ╟─f15d37a7-d962-4da0-977f-76729a3313be
# ╟─28f195c4-4f61-4873-85d6-b4e3aaa3660f
# ╠═3837e439-250b-4577-b0d7-93352ec19f6e
# ╟─8255b5a9-cbe6-4604-bedd-0e366f311096
# ╠═df9b701f-d314-4512-b2ea-1f6ae015166c
# ╟─377f7527-a338-4526-bb24-9766c635e719
# ╟─25f24754-16d4-4343-b54a-cf8ea1358ce9
# ╟─570fd826-23fd-46ee-bdb4-58fb0c45719a
# ╠═691410bb-0472-4800-a90d-29ddf947de3e
# ╟─8cbb1c90-bd94-44b5-80b6-81d38f3e6252
# ╠═c3065acf-6205-455f-ba74-ca51f3f6761b
# ╟─9c392be9-1505-40bc-a290-68085a1c2700
# ╟─d738bdcc-5f83-4dfa-a17f-e9ea23db2986
# ╟─d55ce157-099c-4c1b-94db-62918f04e5fe
# ╟─206f464b-55fe-46aa-85b7-8f0246a0aaad
# ╟─6f72d1b2-63f6-4301-8272-bb2d6d2d049e
# ╠═01201c37-0b79-46b1-a001-e716f5b3ba67
# ╟─e3d37bc1-a119-4add-a111-899ee0caea05
# ╟─6dff3f21-fac0-42e1-910a-f969a231374f
# ╠═b26b8253-e6cd-49f4-81c5-2a3c2963a37c
# ╟─e5dc123f-1596-4311-9398-f0cfe80a5342
# ╟─7b9d26c2-899e-45e9-b664-39d5f1adfe3f
# ╟─fc01d57f-c90b-4231-96be-ddd48656d55e
# ╠═e05f16d6-eb50-49a9-bf14-95d63c9da7ff
# ╟─3b232365-f2fe-4edb-a39f-3e37c8cbb666
# ╠═1992595f-9976-4b18-bf9c-df8a73d30dc8
# ╟─568862b3-6fce-426a-a9e2-e558adf3932a
# ╟─52f9edd0-79b0-4a9a-9930-3a05d3aa2447
# ╟─90d5244e-17be-4601-b922-8c254f1248bf
# ╟─ae79ee89-d788-4612-8b1d-fc22d85c7744
# ╠═3745c237-ba48-4f2c-959e-441484244764
# ╟─6d061411-6f19-4119-aefb-cc380198b0ce
# ╟─8def87d2-f10b-4a82-b353-a6477eeead9b
# ╟─c74b3105-f480-4688-b85f-3e7dff70da3b
# ╠═55438c09-1d94-4ff7-90c3-0cc6064a091e
# ╟─bcc796c9-db11-4a09-a5f9-215127ac0938
# ╟─92b2ecdd-9491-4f16-8cb9-bacbcf180280
# ╠═83216915-cdcf-4f0f-829f-a5ff4c4b8da0
# ╟─452872f9-2009-4733-b2d3-28f262ae19b7
# ╟─29415ddc-e002-4f56-a169-95f7b1c36be9
# ╟─fb23d6c6-b812-4fe1-b224-0014bedbd43f
# ╠═1e53aa10-dff6-40d5-89e2-da194ffc2052
# ╠═ac87ced2-986f-4e6c-80b2-104b25c171c2
# ╠═14cca8ce-cc61-4fae-b871-21c3fd23d0ea
# ╠═28fc8de4-749b-4093-b32f-c398f8d27d3d
# ╠═57b422e5-0ad0-4674-bdd3-a8358bc7aaeb
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
# sMSC test
dc = hcat(2* sin.(0:pi/24:8*pi), 3* sin.(0:pi/24:8*pi));
sMSC_dc = sMSC(dc, 48);
@test all(abs.(sMSC_dc) .< 10e-10)

# TDE test
x = [0,0,0,5,0,0,5]
@test all(TDE(x, 3, 2) .== hcat([5,0,0,5],[0,0,0,5]))
# Dim 2
@test all(TDE(reshape(x, 7, 1), 3, 2)        .== reshape(hcat([5,0,0,5],[0,0,0,5]), (4,2)))
# Dim 3
@test all(TDE(reshape(x, 7, 1, 1), 3, 2)     .== reshape(hcat([5,0,0,5],[0,0,0,5]), (4,1,2)))
# Dim 4
@test all(TDE(reshape(x, 7, 1, 1, 1), 3, 2)  .== reshape(hcat([5,0,0,5],[0,0,0,5]), (4,1,1,2)))

# mw_AVG
x = [4.0, 1.0, 1.0, 1.0, 4.0]
@test all(round.(Int, mw_AVG(x, 3)[2:4]) .== [2,1,2])

# mw_VAR
x = [3.0, -3.0, 4.0, -4.0, 5.0]
@test all(round.(Int, mw_VAR(x, 2)[2:5] .* 2) .== [20,25,34, 41])

# EWMA
x = [10.0, 0.0, 0.0, 0.0, 0.0, 10.0, 0.0]
@test all(round.(Int, EWMA(x, 0.5) * 4) .== [20, 10, 5, 2, 1, 21, 10])

# map moving window
x = randn(144, 1)
@test all(reshape(EWMA(reshape(x, 48, 3)')', 144, 1) .== mapMovingWindow(EWMA, x,ObsPerYear = 48, windowsize = 1, numoutvars = 1))
@test all(reshape(EWMA(reshape(x, 48, 3)', 0.5)', 144, 1) .== mapMovingWindow(EWMA, x, 0.5,ObsPerYear = 48, windowsize = 1, numoutvars = 1))

# get median cycle
 x = sind.(0:8:719)
 @test all(get_MedianCycle(x, 45) .== x[1:45])

ext: jl


text:
module TypedTables

using Tables
using SplitApplyCombine

using Base: @propagate_inbounds, @pure, OneTo, Fix2
import Tables.columns, Tables.rows

export Table, FlexTable, columns, rows, columnnames, ncolumns, showtable

# GetProperty
struct GetProperty{name}
end
@inline GetProperty(name::Symbol) = GetProperty{name}()

@inline function Base.getproperty(sym::Symbol)
	return GetProperty(sym)
end

@inline (::GetProperty{name})(x) where {name} = getproperty(x, name)

# Resultant element type of given column arrays
@generated function _eltypes(a::NamedTuple{names, T}) where {names, T <: Tuple{Vararg{AbstractArray}}}
    Ts = []
    for V in T.parameters
        push!(Ts, eltype(V))
    end
    return NamedTuple{names, Tuple{Ts...}}
end

_ndims(a::NamedTuple{<:Any, T}) where {T} = _ndims(T)
_ndims(::Type{<:Tuple{Vararg{AbstractArray{<:Any, n}}}}) where {n} = n

include("Table.jl")
include("FlexTable.jl")
include("columnops.jl")
include("show.jl")

end # module


text:
#include("../src/likelihood_v2.jl")
include("../src/likelihood_and_diff.jl")

using Distributions, TransformVariables

gamma_mv(mu,std) = Gamma(mu^2/std^2,std^2/mu)
gammainv_mv(mu,std) = Gamma(mu^2/std^2,(std^2/mu)^(-1))

prior_bet = Beta(100,1.01010101)
prior_epsilon = Normal(4,1.5)
prior_theta = Beta(3,1.5)
prior_sig = gamma_mv(2,0.5)
prior_phi = gamma_mv(2,0.5)
prior_phi_pi = gamma_mv(1,0.4)
prior_phi_y = gamma_mv(0.5/4,0.2)
prior_rho_v = Beta(0.88,0.88)
prior_s2 = InverseGamma(3.76,2.76)

#gammainv(alpha,mu) = InverseGamma(alpha,(alpha-1)*mu)

#Plots.plot(3.75:0.05:3.8, std.(gammainv.(3.75:0.05:3.8,1)))

# function posterior(par,data)
#     aa,bb,llh = log_like_dsge(par,data)
#     return llh+logpdf(prior_bet,par[2])+logpdf(prior_epsilon,par[3])+logpdf(prior_theta,par[4])+logpdf(prior_sig,par[5])+logpdf(prior_s2,par[6])+logpdf(prior_phi,par[7])+logpdf(prior_phi_pi,par[8])+logpdf(prior_phi_y,par[9])+logpdf(prior_rho_v,par[10])
# end

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
"""
    AbstractCompartimentalModel

An abstract type representing an arbitrary compartimental model.
"""
abstract type  AbstractCompartimentalModel end

num_states(::AbstractCompartimentalModel) = _NI("num_states")

# Metapopulation

rate(i, state, em::Metapopulation{AbstractCompartimentalModel}) = _NI("rate")

function init_rates!(a, state, mp::Metapopulation{T}) where T <: AbstractCompartimentalModel
    N = nv(mp.h)
    for k in 1:N
        a[k] = rate(k, state, mp)
    end
end

output_elem_size(mp::Metapopulation{T}) where T <: AbstractCompartimentalModel = (nv(mp.h), num_states(mp.dynamics))

function init_output(mp::Metapopulation{T}, state, nmax) where T <: AbstractCompartimentalModel
    M = nv(mp.h)
    m = num_states(mp.dynamics)
    output = [Array{Int,2}(undef, M, m) for i in 1:nmax]
    output[1] .= state
    return output
end

function update_output!(output, state, n, k, mp::Metapopulation{T}) where T <: AbstractCompartimentalModel
    output[n] .= state
end

function finalize_output(mp::Metapopulation{T}, output) where T <: AbstractCompartimentalModel
    M = nv(mp.h)
    n = length(output)
    m = num_states(mp.dynamics)
    out = [Array{Int,2}(undef, n, M) for i in 1:m]
    for i in 1:n, j in 1:m
        out[j][i,:] .= output[i][:,j]
    end
    return out
end

function finalize_meanfield(mp::Metapopulation{T}, sol) where T <: AbstractCompartimentalModel
    m = num_states(mp.dynamics)
    n = length(sol.u)
    M = nv(mp.h)
    output = [Array{Float64,2}(undef, n, M) for i in 1:m]
    for i in 1:n, j in 1:m
        output[j][i,:] .= sol.u[i][:,j]
    end
    return sol.t, output
end

# Contact Process

function init_rates!(a, state, cp::ContactProcess{T}) where T <: AbstractCompartimentalModel
    N = nv(cp.g)
    for k in 1:N
        a[k] = rate(k, state, cp)
    end
end

function init_output(cp::ContactProcess{T}, state, nmax) where T <: AbstractCompartimentalModel
    m = num_states(cp.dynamics)
    output = [Array{Int,1}(undef, m) for i in 1:nmax]
    for k in 1:m
        output[1][k] = length(filter(x->x==k, state))
    end
    return output
end

function update_output!(output, state, n, k, cp::ContactProcess{T}) where T <: AbstractCompartimentalModel
    m = num_states(cp.dynamics)
    for k in 1:m
        output[n][k] = length(filter(x->x==k, state))
    end
end

function finalize_output(cp::ContactProcess{T}, output) where T <: AbstractCompartimentalModel
    n = length(output)
    m = num_states(cp.dynamics)
    out = [Array{Int,1}(undef, n) for i in 1:m]
    for i in 1:n, j in 1:m
        out[j][i] = output[i][j]
    end
    return out
end

output_elem_size(cp::ContactProcess{<:AbstractCompartimentalModel}) = num_states(cp.dynamics)
output_type(::ContactProcess{<:AbstractCompartimentalModel}) = Int

function init_state_mf(cp::ContactProcess{T}, x0) where T <: AbstractCompartimentalModel
    N = nv(cp.g)
    m = num_states(cp.dynamics)
    state = zeros(N,m)
    for i in 1:N
        state[i, x0[i]] = 1.0
    end
    return state
end

function finalize_meanfield(cp::ContactProcess{T}, sol) where T <: AbstractCompartimentalModel
    N = nv(cp.g)
    m = num_states(cp.dynamics)
    n = length(sol.u)
    output = [Array{Float64,2}(undef, n, N) for i in 1:m]
    for i in 1:n, j in 1:m
        output[j][i,:] .= sol.u[i][:,j]
    end
    return sol.t, output
end

# Metaplex

function init_state(mpx::Metaplex{T}, state_0) where T <: AbstractCompartimentalModel
    N = nv(mpx.g)
    M = nv(mpx.h)
    m = num_states(mpx.dynamics)
    x_i = copy(state_0[1]) # epidemic state
    x_μ = copy(state_0[2]) # location
    popcounts = zeros(Int, M, m)
    for k in 1:N
        popcounts[x_μ[k], x_i[k]] += 1
    end
    return Dict(:state_i =>x_i, :state_μ => x_μ, :popcounts => popcounts)
end

function init_rates!(a, state, mpx::Metaplex{T}) where T <: AbstractCompartimentalModel
    N = nv(mpx.g)
    for k in 1:N
        a[k] = rate(k, state, mpx)
    end
end

output_elem_size(mp::Metaplex{T}) where T <: AbstractCompartimentalModel = (nv(mp.h), num_states(mp.dynamics))

function init_output(mpx::Metaplex{T}, state, nmax) where T <: AbstractCompartimentalModel
    M = nv(mpx.h)
    m = num_states(mpx.dynamics)
    output = [Array{Int,2}(undef, M, m) for i in 1:nmax]
    output[1] .= state[:popcounts]
    return output
end

function update_output!(output, state, n, k, mpx::Metaplex{T}) where T <: AbstractCompartimentalModel
    output[n] .= state[:popcounts]
end

finalize_output(mpx::Metaplex{T}, output) where T <: AbstractCompartimentalModel = finalize_output(Metapopulation(mpx), output)

function init_state_mf(mpx::Metaplex{T}, x0) where T <: AbstractCompartimentalModel
    N = nv(mpx.g)
    M = nv(mpx.h)
    m = num_states(mpx.dynamics)
    x_i = x0[1]
    x_μ = x0[2]
    state = zeros(m,N,M)
    for i in 1:N
        state[x_i[i],i,x_μ[i]] = 1.0
    end
    return state
end

function finalize_meanfield(mpx::Metaplex{T}, sol) where T <: AbstractCompartimentalModel
    N = nv(mpx.g)
    M = nv(mpx.h)
    m = num_states(mpx.dynamics)
    n = length(sol.u)
    output = [Array{Float64,2}(undef, n, M) for i in 1:m]
    for i in 1:n, j in 1:m, k in 1:M
        output[j][i,k] = sum(sol.u[i][j,:,k])
    end
    return sol.t, output
end

# Heterogeneous Metapopulation

function init_rates!(a, state, mp::HeterogeneousMetapopulation{T}) where T <: AbstractCompartimentalModel
    N = nv(mp.h)
    for k in 1:N
        a[k] = rate(k, state, mp)
    end
end

output_elem_size(mp::HeterogeneousMetapopulation{T}) where T <: AbstractCompartimentalModel = (nv(mp.h), num_states(mp.dynamics))

function init_output(mp::HeterogeneousMetapopulation{T}, state, nmax) where T <: AbstractCompartimentalModel
    M = nv(mp.h)
    m = num_states(mp.dynamics)
    output = [Array{Int,2}(undef, M, m) for i in 1:nmax]
    output[1] .= state
    return output
end

function update_output!(output, state, n, k, mp::HeterogeneousMetapopulation{T}) where T <: AbstractCompartimentalModel
    output[n] .= state
end

function finalize_output(mp::HeterogeneousMetapopulation{T}, output) where T <: AbstractCompartimentalModel
    M = nv(mp.h)
    n = length(output)
    m = num_states(mp.dynamics)
    out = [Array{Int,2}(undef, n, M) for i in 1:m]
    for i in 1:n, j in 1:m
        out[j][i,:] .= output[i][:,j]
    end
    return out
end

function finalize_meanfield(mp::HeterogeneousMetapopulation{T}, sol) where T <: AbstractCompartimentalModel
    m = num_states(mp.dynamics)
    n = length(sol.u)
    M = nv(mp.h)
    output = [Array{Float64,2}(undef, n, M) for i in 1:m]
    for i in 1:n, j in 1:m
        output[j][i,:] .= sol.u[i][:,j]
    end
    return sol.t, output
end


text:
#=
   Locomotive problem.

   From Think Bayes, page 22ff
   """
   I found the locomotive problem in Frederick Mosteller's "Fifty Challenging
   Problems in Probability with Solutions" (Dover, 1987):
      'A railroad numbers its locomotives in order 1..N. One day you
       see a locomotive with the number 60. Estimate how many loco-
       motives the railroad has.'

   ...

   The mean of the posterior is 333, so that might be a good guess if you wanted to
   minimize error.
   """

   As the book (Think Bayes) mentions, the model is very sensitive to
   the value of maxInt. Here are some results for different maxInt:

   Here are the values from this model for different values of maxInt
   (using MH())
     maxInt:100  mean(n): 78.23232
     maxInt:200  mean(n): 115.87732
     maxInt:500  mean(n): 206.42599
     maxInt:1000  mean(n): 335.66271
     maxInt:2000  mean(n): 552.68673

   These are the exact values from WebPPL model using the enumerate sampler.
   maxInt  n
   -------------------------
    100    77.8166306057923
    200   115.84577808279282
    500   207.2393675826458
   1000   334.04414386751915
   2000   553.5237331955558

   Mosteller's solution: 2*(60-1)+1: 119.

  cf german_tank.jl
     german_tank_int.jl
     ~/blog/locomotive_problem.blog
     ~/webppl/locomotive_problem.wppl
=#

using Turing
include("jl_utils.jl")

@model function locomotive_problem(ys,maxInt=1000)
    len = length(ys)
    n ~ DiscreteUniform(1,maxInt)
    ys ~ filldist(DiscreteUniform(1,n),len)
end

# Run the model with different values of maxInt
function run_locomotive_problem(ys,maxInt)
   model = locomotive_problem(ys,maxInt)

    # chns = sample(model, Prior(), 100_000)
    chns = sample(model, MH(), 100_000)
    # chns = sample(model, PG(15), 10_000)
    # chns = sample(model, SMC(), 10_000)
    # chns = sample(model, IS(), 100_000)
    # display(chns)

    println("maxInt:$maxInt  mean(n): ", mean_val(chns,:n))
end

ys = [60]
for maxInt in [100,200,500,1000,2000]
    run_locomotive_problem(ys,maxInt)
end

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
plotstuff_linear(args...) = println("Install package Plots.jl to plot results in the end of demo_linear")

@require Plots begin
function plotstuff_linear(x,u,cost,totalcost)
    p = Plots.subplot(n=4,nc=2)
    Plots.plot!(p[1,1],x', title="State Trajectories",legend=false)
    Plots.plot!(p[1,2],cost',c=:black,linewidth=3, title="Cost", xlabel="Time step")
    Plots.plot!(p[2,1],u',title="Control signals")
    Plots.plot!(p[2,2],totalcost,title="Total cost")
end
end

function demo_linear()
    println("Running linear demo function for DifferentialDynamicProgramming.jl")

    # make stable linear dynamics
    h = .01         # time step
    n = 10          # state dimension
    m = 2           # control dimension
    A = randn(n,n)
    A = A-A'        # skew-symmetric = pure imaginary eigenvalues
    A = expm(h*A)   # discrete time
    B = h*randn(n,m)

    # quadratic costs
    Q = h*eye(n)
    R = .1*h*eye(m)

    # control limits
    lims = [] #ones(m,1)*[-1 1]*.6

    T        = 1000              # horizon
    x0       = ones(n,1)        # initial state
    u0       = .1*randn(m,T)     # initial controls

    # optimization problem
    N   = T+1
    fx  = A
    fu  = B
    cxx = Q
    cxu = zeros(size(B))
    cuu = R

    function lin_dyn_df(x,u,Q,R)
        u[isnan(u)] = 0
        cx  = Q*x
        cu  = R*u
        fxx=fxu=fuu = []
        return fx,fu,fxx,fxu,fuu,cx,cu,cxx,cxu,cuu
    end
    function lin_dyn_f(x,u,A,B,Q,R)
        u[isnan(u)] = 0
        f = A*x + B*u
        c = 0.5*sum(x.*(Q*x)) + 0.5*sum(u.*(R*u))
        return f,c
    end

    function lin_dyn_fT(x,Q)
        c = 0.5*sum(x.*(Q*x))
        return c
    end

    f(x,u,i)   = lin_dyn_f(x,u,A,B,Q,R)
    fT(x)      = lin_dyn_fT(x,Q)
    df(x,u,i)  = lin_dyn_df(x,u,Q,R)
    # plotFn(x)  = plot(squeeze(x,2)')

    # run the optimization
    @time x, u, L, Vx, Vxx, cost, otrace = iLQG(f,fT,df, x0, u0, lims=lims, plotFn= x -> 0 );

    totalcost = [ t.cost for t in otrace]
    iters = sum(totalcost .> 0)
    totalcost = [ otrace[i].cost for i=1:iters]

    plotstuff_linear(x,u,cost,totalcost)

end


text:
const TIME_SERIES_STORAGE_FILE = "time_series_storage.h5"
const VALIDATION_DESCRIPTOR_FILE = "validation_descriptors.json"

"""
    mutable struct SystemData <: InfrastructureSystemsType
        components::Components
        "Masked components are attached to the system for overall management purposes but
        are not exposed in the standard library calls like [`get_components`](@ref).
        Examples are components in a subsystem."
        masked_components::Components
        time_series_params::TimeSeriesParameters
        validation_descriptors::Vector
        time_series_storage::TimeSeriesStorage
        time_series_storage_file::Union{Nothing, String}
        internal::InfrastructureSystemsInternal
    end

Container for system components and time series data
"""
mutable struct SystemData <: InfrastructureSystemsType
    components::Components
    masked_components::Components
    time_series_params::TimeSeriesParameters
    time_series_storage::TimeSeriesStorage
    validation_descriptors::Vector
    internal::InfrastructureSystemsInternal
end

"""
Construct SystemData to store components and time series data.

# Arguments
- `validation_descriptor_file = nothing`: Optionally, a file defining component validation
  descriptors.
- `time_series_in_memory = false`: Controls whether time series data is stored in memory or
  in a file.
- `time_series_directory = nothing`: Controls what directory time series data is stored in.
  Default is tempdir().
- `compression = CompressionSettings()`: Controls compression of time series data.
"""
function SystemData(;
    validation_descriptor_file = nothing,
    time_series_in_memory = false,
    time_series_directory = nothing,
    compression = CompressionSettings(),
)
    if isnothing(validation_descriptor_file)
        validation_descriptors = Vector()
    else
        validation_descriptors = read_validation_descriptor(validation_descriptor_file)
    end

    ts_storage = make_time_series_storage(;
        in_memory = time_series_in_memory,
        directory = time_series_directory,
        compression = compression,
    )
    components = Components(ts_storage, validation_descriptors)
    masked_components = Components(ts_storage, validation_descriptors)
    return SystemData(
        components,
        masked_components,
        TimeSeriesParameters(),
        ts_storage,
        validation_descriptors,
        InfrastructureSystemsInternal(),
    )
end

function SystemData(
    time_series_params,
    validation_descriptors,
    time_series_storage,
    internal,
)
    components = Components(time_series_storage, validation_descriptors)
    masked_components = Components(time_series_storage, validation_descriptors)
    return SystemData(
        components,
        masked_components,
        time_series_params,
        time_series_storage,
        validation_descriptors,
        internal,
    )
end

"""
Adds time_series from a metadata file or metadata descriptors.

# Arguments
- `data::SystemData`: system
- `::Type{T}`: type of the component associated with time series data; may be abstract
- `metadata_file::AbstractString`: metadata file for time series
  that includes an array of TimeSeriesFileMetadata instances or a vector.
- `resolution::DateTime.Period=nothing`: skip time_series that don't match this resolution.
"""
function add_time_series_from_file_metadata!(
    data::SystemData,
    ::Type{T},
    metadata_file::AbstractString;
    resolution = nothing,
) where {T <: InfrastructureSystemsComponent}
    metadata = read_time_series_file_metadata(metadata_file)
    return add_time_series_from_file_metadata!(data, T, metadata; resolution = resolution)
end

"""
Adds time series data from a metadata file or metadata descriptors.

# Arguments
- `data::SystemData`: system
- `file_metadata::Vector{TimeSeriesFileMetadata}`: metadata for time series
- `resolution::DateTime.Period=nothing`: skip time_series that don't match this resolution.
"""
function add_time_series_from_file_metadata!(
    data::SystemData,
    ::Type{T},
    file_metadata::Vector{TimeSeriesFileMetadata};
    resolution = nothing,
) where {T <: InfrastructureSystemsComponent}
    cache = TimeSeriesParsingCache()
    for metadata in file_metadata
        if resolution === nothing || metadata.resolution == resolution
            add_time_series_from_file_metadata_internal!(data, T, cache, metadata)
        end
    end
    return
end

"""
Add time series data to a component.

# Arguments
- `data::SystemData`: SystemData
- `component::InfrastructureSystemsComponent`: will store the time series reference
- `time_series::TimeSeriesData`: Any object of subtype TimeSeriesData

Throws ArgumentError if the component is not stored in the system.

"""
function add_time_series!(
    data::SystemData,
    component::InfrastructureSystemsComponent,
    time_series::TimeSeriesData;
    skip_if_present = false,
)
    metadata_type = time_series_data_to_metadata(typeof(time_series))
    ts_metadata = metadata_type(time_series)
    _attach_time_series_and_serialize!(
        data,
        component,
        ts_metadata,
        time_series;
        skip_if_present = skip_if_present,
    )
    return
end

"""
Add the same time series data to multiple components.

# Arguments
- `data::SystemData`: SystemData
- `components`: iterable of components that will store the same time series reference
- `time_series::TimeSeriesData`: Any object of subtype TimeSeriesData

This is significantly more efficent than calling `add_time_series!` for each component
individually with the same data because in this case, only one time series array is stored.

Throws ArgumentError if a component is not stored in the system.
"""
function add_time_series!(data::SystemData, components, time_series::TimeSeriesData)
    metadata_type = time_series_data_to_metadata(typeof(time_series))
    ts_metadata = metadata_type(time_series)
    for component in components
        _attach_time_series_and_serialize!(data, component, ts_metadata, time_series)
    end
end

function _attach_time_series_and_serialize!(
    data::SystemData,
    component::InfrastructureSystemsComponent,
    ts_metadata::T,
    ts::TimeSeriesData;
    skip_if_present = false,
) where {T <: TimeSeriesMetadata}
    _validate_component(data, component)
    check_add_time_series(data.time_series_params, ts)
    check_read_only(data.time_series_storage)
    if has_time_series(component, T, get_name(ts))
        skip_if_present && return
        throw(ArgumentError("time_series $(typeof(ts)) $(get_name(ts)) is already stored"))
    end

    serialize_time_series!(
        data.time_series_storage,
        get_uuid(component),
        get_name(ts_metadata),
        ts,
    )
    add_time_series!(component, ts_metadata, skip_if_present = skip_if_present)
    # Order is important. Set this last in case exceptions are thrown at previous steps.
    set_parameters!(data.time_series_params, ts)
    return
end

function add_time_series_from_file_metadata_internal!(
    data::SystemData,
    ::Type{T},
    cache::TimeSeriesParsingCache,
    file_metadata::TimeSeriesFileMetadata,
) where {T <: InfrastructureSystemsComponent}
    set_component!(file_metadata, data, InfrastructureSystems)
    component = file_metadata.component
    time_series = make_time_series!(cache, file_metadata)
    add_time_series!(data, component, time_series)
    return
end

"""
Remove the time series data for a component.
"""
function remove_time_series!(
    data::SystemData,
    ::Type{T},
    component::InfrastructureSystemsComponent,
    name::String,
) where {T <: TimeSeriesData}
    type = time_series_data_to_metadata(T)
    time_series = get_time_series(type, component, name)
    uuid = get_time_series_uuid(time_series)
    if remove_time_series_metadata!(component, type, name)
        remove_time_series!(data.time_series_storage, uuid, get_uuid(component), name)
    end

    return
end

"""
Return a time series from TimeSeriesFileMetadata.

# Arguments
- `cache::TimeSeriesParsingCache`: cached data
- `ts_file_metadata::TimeSeriesFileMetadata`: metadata
- `resolution::{Nothing, Dates.Period}`: skip any time_series that don't match this resolution
"""
function make_time_series!(
    cache::TimeSeriesParsingCache,
    ts_file_metadata::TimeSeriesFileMetadata,
)
    info = add_time_series_info!(cache, ts_file_metadata)
    return ts_file_metadata.time_series_type(info)
end

function add_time_series_info!(
    cache::TimeSeriesParsingCache,
    metadata::TimeSeriesFileMetadata,
)
    time_series = _add_time_series_info!(cache, metadata)
    info = TimeSeriesParsedInfo(metadata, time_series)
    @debug "Added TimeSeriesParsedInfo" _group = LOG_GROUP_TIME_SERIES metadata
    return info
end

"""
Checks that the component exists in data and the UUID's match.
"""
function _validate_component(
    data::SystemData,
    component::T,
) where {T <: InfrastructureSystemsComponent}
    name = get_name(component)
    comp = get_component(T, data.components, name)
    if isnothing(comp)
        comp = get_masked_component(T, data, name)
        if comp === nothing
            throw(ArgumentError("no $T with name=$name is stored"))
        end
    end

    user_uuid = get_uuid(component)
    ps_uuid = get_uuid(comp)
    if user_uuid != ps_uuid
        throw(
            ArgumentError(
                "comp UUID doesn't match, perhaps it was copied?; " *
                "$T name=$(get_name(component)) user=$user_uuid system=$ps_uuid",
            ),
        )
    end
end

function compare_values(x::SystemData, y::SystemData)::Bool
    match = true
    for name in fieldnames(SystemData)
        if name == :components || name == :masked_components
            # Not deserialized in IS.
            continue
        end
        val_x = getfield(x, name)
        val_y = getfield(y, name)
        if name == :time_series_storage && typeof(val_x) != typeof(val_y)
            # TODO 1.0: workaround for not being able to convert Hdf5TimeSeriesStorage to
            # InMemoryTimeSeriesStorage
            continue
        end
        if !compare_values(val_x, val_y)
            @error "SystemData field=$name does not match" getfield(x, name) getfield(
                y,
                name,
            )
            match = false
        end
    end

    return match
end

function remove_component!(::Type{T}, data::SystemData, name) where {T}
    return remove_component!(T, data.components, name)
end

function remove_component!(data::SystemData, component)
    remove_component!(data.components, component)
end

function remove_components!(::Type{T}, data::SystemData) where {T}
    return remove_components!(T, data.components)
end

"""
Removes the component from the main container and adds it to the masked container.
"""
function mask_component!(data::SystemData, component::InfrastructureSystemsComponent)
    remove_component!(data.components, component, remove_time_series = false)
    set_time_series_storage!(component, nothing)
    add_masked_component!(
        data,
        component,
        skip_validation = true,  # validation has already occurred
        allow_existing_time_series = true,
    )
end

function clear_time_series!(data::SystemData)
    clear_time_series!(data.time_series_storage)
    clear_time_series!(data.components)
    reset_info!(data.time_series_params)
end

function iterate_components_with_time_series(data::SystemData)
    return Iterators.Flatten((
        iterate_components_with_time_series(data.components),
        iterate_components_with_time_series(data.masked_components),
    ))
end

"""
Removes all time series of a particular type from a System.

# Arguments
- `data::SystemData`: system
- `type::Type{<:TimeSeriesData}`: Type of time series objects to remove.
"""
function remove_time_series!(data::SystemData, ::Type{T}) where {T <: TimeSeriesData}
    for component in iterate_components_with_time_series(data)
        for ts in get_time_series_multiple(component, type = T)
            remove_time_series!(data, typeof(ts), component, get_name(ts))
        end
    end
    counts = get_time_series_counts(data)
    if counts[3] == 0 # no more forecast objects
        if counts[2] == 0 # no more static time series objects
            reset_info!(data.time_series_params)
        else
            reset_info!(data.time_series_params.forecast_params)
        end
    end
end

"""
Returns an iterator of TimeSeriesData instances attached to the system.

Note that passing a filter function can be much slower than the other filtering parameters
because it reads time series data from media.

Call `collect` on the result to get an array.

# Arguments
- `data::SystemData`: system
- `filter_func = nothing`: Only return time_series for which this returns true.
- `type = nothing`: Only return time_series with this type.
- `name = nothing`: Only return time_series matching this value.
"""
function get_time_series_multiple(
    data::SystemData,
    filter_func = nothing;
    type = nothing,
    name = nothing,
)
    Channel() do channel
        for component in iterate_components_with_time_series(data)
            for time_series in
                get_time_series_multiple(component, filter_func; type = type, name = name)
                put!(channel, time_series)
            end
        end
    end
end

"""
Return a tuple of counts of components with time series and total time series and forecasts.
"""
function get_time_series_counts(data::SystemData)
    component_count = 0
    static_time_series_count = 0
    forecast_count = 0
    for component in iterate_components_with_time_series(data)
        component_count += 1
        _ts_count, _forecast_count = get_num_time_series(component)
        static_time_series_count += _ts_count
        forecast_count += _forecast_count
    end

    return (component_count, static_time_series_count, forecast_count)
end

"""
Transform all instances of SingleTimeSeries to DeterministicSingleTimeSeries.

Any existing DeterministicSingleTimeSeries forecasts will be deleted even if the inputs are
invalid.
"""
function transform_single_time_series!(
    data::SystemData,
    ::Type{T},
    horizon::Int,
    interval::Dates.Period,
) where {T <: DeterministicSingleTimeSeries}
    remove_time_series!(data, DeterministicSingleTimeSeries)
    params = nothing
    set_params = false
    for component in iterate_components_with_time_series(data)
        if params === nothing
            params = get_single_time_series_transformed_parameters(
                component,
                T,
                horizon,
                interval,
            )
            # This will throw if there is another forecast type with conflicting parameters.
            check_params_compatibility(data.time_series_params, params)
        end

        if transform_single_time_series_internal!(component, T, params) && !set_params
            set_parameters!(data.time_series_params, params)
            set_params = true
        end
    end
end

"""
Set the component value in metadata by looking up the category in module.
This requires that category be a string version of a component's abstract type.
Modules can override for custom behavior.
"""
function set_component!(metadata::TimeSeriesFileMetadata, data::SystemData, mod::Module)
    category = getfield(mod, Symbol(metadata.category))
    if isconcretetype(category)
        metadata.component =
            get_component(category, data.components, metadata.component_name)
        if isnothing(metadata.component)
            @warn "no component category=$category name=$(metadata.component_name)"
        end
    else
        # Note: this could dispatch to higher-level modules that reimplement it.
        components = get_components_by_name(category, data, metadata.component_name)
        if length(components) == 0
            @warn "no component category=$category name=$(metadata.component_name)"
            metadata.component = nothing
        elseif length(components) == 1
            metadata.component = components[1]
        else
            throw(
                DataFormatError(
                    "duplicate names type=$(category) name=$(metadata.component_name)",
                ),
            )
        end
    end
end

"""
Parent object should call this prior to serialization so that SystemData can store the
appropriate path information for the time series data.
"""
function prepare_for_serialization!(
    data::SystemData,
    filename::AbstractString;
    force = false,
)
    directory = dirname(filename)
    if !isdir(directory)
        mkpath(directory)
    end

    files = [
        filename,
        joinpath(directory, TIME_SERIES_STORAGE_FILE),
        joinpath(directory, VALIDATION_DESCRIPTOR_FILE),
    ]
    for file in files
        if !force && isfile(file)
            error("$file already exists. Set force=true to overwrite.")
        end
    end

    ext = get_ext(data.internal)
    ext["serialization_directory"] = directory
    ext["basename"] = splitext(basename(filename))[1]
end

function serialize(data::SystemData)
    @debug "serialize SystemData" _group = LOG_GROUP_SERIALIZATION
    json_data = Dict()
    for field in (:components, :masked_components, :time_series_params, :internal)
        json_data[string(field)] = serialize(getfield(data, field))
    end

    ext = get_ext(data.internal)
    if !haskey(ext, "serialization_directory")
        error("prepare_for_serialization! was not called")
    end
    directory = pop!(ext, "serialization_directory")
    base = pop!(ext, "basename")
    isempty(ext) && clear_ext!(data.internal)

    time_series_base_name = base * "_" * TIME_SERIES_STORAGE_FILE
    time_series_storage_file = joinpath(directory, time_series_base_name)
    serialize(data.time_series_storage, time_series_storage_file)
    json_data["time_series_storage_file"] = time_series_base_name
    json_data["time_series_storage_type"] = string(typeof(data.time_series_storage))

    descriptor_base_name = base * "_" * VALIDATION_DESCRIPTOR_FILE
    descriptor_file = joinpath(directory, descriptor_base_name)
    descriptors = Dict("struct_validation_descriptors" => data.validation_descriptors)
    text = JSON3.write(descriptors)
    open(descriptor_file, "w") do io
        write(io, text)
    end
    json_data["validation_descriptor_file"] = descriptor_base_name
    json_data["version_info"] = serialize_julia_info()
    return json_data
end

function deserialize(
    ::Type{SystemData},
    raw::Dict;
    time_series_read_only = false,
    time_series_directory = nothing,
)
    @debug "deserialize" raw _group = LOG_GROUP_SERIALIZATION
    time_series_params = deserialize(TimeSeriesParameters, raw["time_series_params"])
    # The code calling this function must have changed to this directory.
    if !isfile(raw["time_series_storage_file"])
        error("time series file $(raw["time_series_storage_file"]) does not exist")
    end
    if !isfile(raw["validation_descriptor_file"])
        error(
            "validation descriptor file $(raw["validation_descriptor_file"]) does not exist",
        )
    end
    validation_descriptors = read_validation_descriptor(raw["validation_descriptor_file"])

    # TODO 1.0: need to address this limitation
    if strip_module_name(raw["time_series_storage_type"]) == "InMemoryTimeSeriesStorage"
        @info "Deserializing with InMemoryTimeSeriesStorage is currently not supported. Using HDF"
        #hdf5_storage = Hdf5TimeSeriesStorage(raw["time_series_storage_file"], true)
        #time_series_storage = InMemoryTimeSeriesStorage(hdf5_storage)
    end
    time_series_storage = from_file(
        Hdf5TimeSeriesStorage,
        raw["time_series_storage_file"];
        read_only = time_series_read_only,
        directory = time_series_directory,
    )

    internal = deserialize(InfrastructureSystemsInternal, raw["internal"])
    @debug "deserialize" _group = LOG_GROUP_SERIALIZATION validation_descriptors time_series_storage internal
    sys = SystemData(
        time_series_params,
        validation_descriptors,
        time_series_storage,
        internal,
    )
    # Note: components need to be deserialized by the parent so that they can go through
    # the proper checks.
    return sys
end

# Redirect functions to Components and TimeSeriesContainer

add_component!(data::SystemData, component; kwargs...) =
    add_component!(data.components, component; kwargs...)

add_masked_component!(data::SystemData, component; kwargs...) = add_component!(
    data.masked_components,
    component;
    allow_existing_time_series = true,
    kwargs...,
)

iterate_components(data::SystemData) = iterate_components(data.components)

get_component(::Type{T}, data::SystemData, args...) where {T} =
    get_component(T, data.components, args...)

function get_component(data::SystemData, uuid::Base.UUID)
    for component in get_components(InfrastructureSystemsComponent, data)
        if get_uuid(component) == uuid
            return component
        end
    end

    @error "no component with UUID $uuid is stored"
    return nothing
end

function get_components(
    ::Type{T},
    data::SystemData,
    filter_func::Union{Nothing, Function} = nothing,
) where {T}
    return get_components(T, data.components, filter_func)
end

get_components_by_name(::Type{T}, data::SystemData, args...) where {T} =
    get_components_by_name(T, data.components, args...)

function get_masked_components(
    ::Type{T},
    data::SystemData,
    filter_func::Union{Nothing, Function} = nothing,
) where {T}
    return get_components(T, data.masked_components, filter_func)
end

get_masked_components_by_name(::Type{T}, data::SystemData, args...) where {T} =
    get_components_by_name(T, data.masked_components, args...)

get_masked_component(::Type{T}, data::SystemData, name) where {T} =
    get_component(T, data.masked_components, name)

function get_masked_component(data::SystemData, uuid::Base.UUID)
    for component in get_masked_components(InfrastructureSystemsComponent, data)
        if get_uuid(component) == uuid
            return component
        end
    end

    @error "no component with UUID $uuid is stored"
    return nothing
end

get_forecast_initial_times(data::SystemData) =
    get_forecast_initial_times(data.time_series_params)
get_forecast_total_period(data::SystemData) =
    get_forecast_total_period(data.time_series_params)
get_forecast_window_count(data::SystemData) =
    get_forecast_window_count(data.time_series_params)
get_forecast_horizon(data::SystemData) = get_forecast_horizon(data.time_series_params)
get_forecast_initial_timestamp(data::SystemData) =
    get_forecast_initial_timestamp(data.time_series_params)
get_forecast_interval(data::SystemData) = get_forecast_interval(data.time_series_params)
get_time_series_resolution(data::SystemData) =
    get_time_series_resolution(data.time_series_params)

clear_components!(data::SystemData) = clear_components!(data.components)

check_components(data::SystemData) = check_components(data.components)
check_component(data::SystemData, component) = check_component(data.components, component)

get_compression_settings(data::SystemData) =
    get_compression_settings(data.time_series_storage)

set_name!(data::SystemData, component, name) = set_name!(data.components, component, name)


text:
using WAV

function load_sound(file, n0)

    x, fs = wavread(file)
    n = length(x)

    if contains(file,"bird.wav")
        sel = [1:6000; 12500:15000; 22500:24000; 32500:34000]
        keep = filter(x -> !(x in sel), 1:length(x))
        x=x[keep]
    end

    if (n0 != 0) && (n0 < n)
        x = x[1:n0]
    end

    return x/maximum(x)
end


text:
"""
    FBMKernel(; h::Real=0.5)

Fractional Brownian motion kernel with Hurst index `h`.

# Definition

For inputs ``x, x' \\in \\mathbb{R}^d``, the fractional Brownian motion kernel with
[Hurst index](https://en.wikipedia.org/wiki/Hurst_exponent#Generalized_exponent)
``h \\in [0,1]`` is defined as
```math
k(x, x'; h) =  \\frac{\\|x\\|_2^{2h} + \\|x'\\|_2^{2h} - \\|x - x'\\|^{2h}}{2}.
```
"""
struct FBMKernel{T<:Real} <: Kernel
    h::Vector{T}
    function FBMKernel(h::Real)
        @check_args(FBMKernel, h, zero(h) ≤ h ≤ one(h), "h ∈ [0, 1]")
        return new{typeof(h)}([h])
    end
end

FBMKernel(; h::Real=0.5) = FBMKernel(h)

@functor FBMKernel

function (κ::FBMKernel)(x::AbstractVector{<:Real}, y::AbstractVector{<:Real})
    modX = sum(abs2, x)
    modY = sum(abs2, y)
    modXY = sqeuclidean(x, y)
    h = first(κ.h)
    return (modX^h + modY^h - modXY^h) / 2
end

function (κ::FBMKernel)(x::Real, y::Real)
    return (abs2(x)^first(κ.h) + abs2(y)^first(κ.h) - abs2(x - y)^first(κ.h)) / 2
end

function Base.show(io::IO, κ::FBMKernel)
    return print(io, "Fractional Brownian Motion Kernel (h = ", first(κ.h), ")")
end

_fbm(modX, modY, modXY, h) = (modX^h + modY^h - modXY^h) / 2

_mod(x::AbstractVector{<:Real}) = abs2.(x)
_mod(x::ColVecs) = vec(sum(abs2, x.X; dims=1))
_mod(x::RowVecs) = vec(sum(abs2, x.X; dims=2))

function kernelmatrix(κ::FBMKernel, x::AbstractVector)
    modx = _mod(x)
    modxx = pairwise(SqEuclidean(), x)
    return _fbm.(modx, modx', modxx, κ.h)
end

function kernelmatrix!(K::AbstractMatrix, κ::FBMKernel, x::AbstractVector)
    modx = _mod(x)
    pairwise!(K, SqEuclidean(), x)
    K .= _fbm.(modx, modx', K, κ.h)
    return K
end

function kernelmatrix(κ::FBMKernel, x::AbstractVector, y::AbstractVector)
    modxy = pairwise(SqEuclidean(), x, y)
    return _fbm.(_mod(x), _mod(y)', modxy, κ.h)
end

function kernelmatrix!(
    K::AbstractMatrix, κ::FBMKernel, x::AbstractVector, y::AbstractVector
)
    pairwise!(K, SqEuclidean(), x, y)
    K .= _fbm.(_mod(x), _mod(y)', K, κ.h)
    return K
end

function kernelmatrix_diag(κ::FBMKernel, x::AbstractVector)
    modx = _mod(x)
    modxx = colwise(SqEuclidean(), x)
    return _fbm.(modx, modx, modxx, κ.h)
end

function kernelmatrix_diag(κ::FBMKernel, x::AbstractVector, y::AbstractVector)
    modxy = colwise(SqEuclidean(), x, y)
    return _fbm.(_mod(x), _mod(y), modxy, κ.h)
end


text:
"""
Gallery of `SeisForwExpt`.
# Arguments
* `attrib::Symbol` : name
  * `=:acou_homo1`
"""
function SeisForwExpt(attrib::Symbol)
	@assert attrib in [:acou_homo1]

	model=Seismic(:acou_homo1);
	acqgeom =Geom(model.mgrid,:xwell);
	tgrid=range(0.0,stop=2.0,length=1000)
	wav=GeoPhyInv.Utils.Wavelets.ricker(10.0, tgrid, tpeak=0.25, );
	# source wavelet for modelling
	acqsrc=GeoPhyInv.Acquisition.Src_fixed(acqgeom.nss,1,[:P],wav,tgrid);

	vp0=mean(GeoPhyInv.Models.χ(model.χvp,model.ref.vp,-1))
	ρ0=mean(GeoPhyInv.Models.χ(model.χρ,model.ref.ρ,-1))

	pa=GeoPhyInv.Fdtd.Param(npw=1,model=model,
	    acqgeom=[acqgeom], acqsrc=[acqsrc],
		sflags=[2], rflags=[1],
		    tgridmod=tgrid, verbose=true);
	return pa

end


text:
module FEFMM_Tests
    using Test
    using FEFMM
    using UnicodePlots
    using LinearAlgebra: norm

    include("testmediums.jl")

    @testset "FEFMM Tests" begin
        @testset "Solve Piecewise Quadratic" begin
            @test isapprox(FEFMM.solve_piecewise_quadratic(1.0,1.0,1.0), 2.0)
            @test isapprox(FEFMM.solve_piecewise_quadratic(1.0,1.0,1.0,1.0,1.0), 1.0+1/sqrt(2))
            @test isapprox(FEFMM.solve_piecewise_quadratic(1.0,1.0,100.0,1.0,1.0), 2.0)
            @test isapprox(FEFMM.solve_piecewise_quadratic(2.0,1.0,1.0,100.0,1.0), 1.0)
            @test isapprox(FEFMM.solve_piecewise_quadratic(1.0,1.0,1.0,1.0,1.0,1.0,1.0), 1.0+1/sqrt(3))
            @test isapprox(FEFMM.solve_piecewise_quadratic(1.0,2.0,3.0,100.0,10.0,1.0,1.0), 2.0/3.0)
            @test isapprox(FEFMM.solve_piecewise_quadratic(1.0,2.0,3.0,10.0,100.0,1.0,1.0), 2.0/3.0)
            @test isapprox(FEFMM.solve_piecewise_quadratic(1.0,2.0,3.0,10.0,1.0,100.0,1.0), 1.0)
            @test isapprox(FEFMM.solve_piecewise_quadratic(1.0,2.0,3.0,1.0,10.0,100.0,1.0), 2.0)
        end

        @testset "Basic FEFMM Checks" begin
            κ21d = ones(101)
            dx1d = [0.1]
            xs1d = CartesianIndex(1)
            τa1d = 0:0.1:10
            @test all(isapprox.(τa1d, FEFMM.fefmm(κ21d, dx1d, xs1d)[1]))
            κ22d = ones(101,101)
            dx2d = [0.1,0.1]
            xs2d = CartesianIndex(1,1)
            τaxis = 0:0.1:10
            τa2d = sqrt.(τaxis'.^2 .+ τaxis.^2)
            @test all(isapprox.(τa2d, FEFMM.fefmm(κ22d, dx2d, xs2d)[1]))
            κ23d = ones(101,101,101)
            dx3d = [0.1,0.1,0.1]
            xs3d = CartesianIndex(1,1,1)
            τa3d = zeros(101,101,101)
            for k in 1:101
                for j in 1:101
                    for i in 1:101
                        τa3d[i,j,k] = sqrt(τaxis[i]^2+τaxis[j]^2+τaxis[k]^2)
                    end
                end
            end
            @test all(isapprox.(τa3d, FEFMM.fefmm(κ23d, dx3d, xs3d)[1]))
        end

        @testset "Constant Gradient of Squared Slowness 2D" begin
            hlist = [1/20, 1/40, 1/80, 1/160, 1/320]#, 1/640, 1/1280]
            linorm = zeros(length(hlist))
            l2norm = zeros(length(hlist))
            for i = 1:length(hlist)
                (x0, k2, t_exact) = const_k2_2D(hlist[i])
                (t_pred, ordering) = FEFMM.fefmm(k2, [hlist[i], hlist[i]], x0)
                err = t_exact.-t_pred
                linorm[i] = norm(err, Inf)
                l2norm[i] = norm(err, 2)/sqrt(length(err))
            end

            plt = lineplot(log10.(hlist),
                           log10.(linorm),
                           title = "Constant Gradient of Squared Slowness 2D",
                           name = "L_inf norm",
                           xlabel = "log_10(h)",
                           ylabel = "log_10(Error)",
                           xlim = [floor(minimum(log10.(hlist))), ceil(maximum(log10.(hlist)))],
                           ylim = [floor(minimum(log10.([linorm; l2norm]))), ceil(maximum(log10.([linorm; l2norm])))],
                           color = :blue)

            lineplot!(plt, log10.(hlist), log10.(l2norm), name = "mean L2 norm", color = :red)
            println(plt)
            @test reverse(linorm) == sort(linorm) #test for decreasing error with decreasing h
            @test reverse(l2norm) == sort(l2norm)
        end

        @testset "Constant Gradient of Squared Velocity 2D" begin
            hlist = [1/20, 1/40, 1/80, 1/160, 1/320]#, 1/640, 1/1280]
            linorm = zeros(length(hlist))
            l2norm = zeros(length(hlist))
            for i = 1:length(hlist)
                (x0, k2, t_exact) = const_v2_2D(hlist[i])
                (t_pred, ordering) = FEFMM.fefmm(k2, [hlist[i], hlist[i]], x0)
                err = t_exact.-t_pred
                linorm[i] = norm(err, Inf)
                l2norm[i] = norm(err, 2)/sqrt(length(err))
            end

                plt = lineplot(log10.(hlist),
                            log10.(linorm),
                            title = "Constant Gradient of Squared Velocity 2D",
                            name = "L_inf norm",
                            xlabel = "log_10(h)",
                            ylabel = "log_10(Error)",
                            xlim = [floor(minimum(log10.(hlist))), ceil(maximum(log10.(hlist)))],
                            ylim = [floor(minimum(log10.([linorm; l2norm]))), ceil(maximum(log10.([linorm; l2norm])))],
                            color = :blue)

            lineplot!(plt, log10.(hlist), log10.(l2norm), name = "mean L2 norm", color = :red)
            println(plt)
            @test reverse(linorm) == sort(linorm) #test for decreasing error with decreasing h
            @test reverse(l2norm) == sort(l2norm)
        end

        @testset "Gaussian Factor 2D" begin
            hlist = [1/20, 1/40, 1/80, 1/160, 1/320]#, 1/640, 1/1280]
            linorm = zeros(length(hlist))
            l2norm = zeros(length(hlist))
            for i = 1:length(hlist)
                (x0, k2, t_exact) = gaussian_factor_2D(hlist[i])
                (t_pred, ordering) = FEFMM.fefmm(k2, [hlist[i], hlist[i]], x0)
                err = t_exact.-t_pred
                linorm[i] = norm(err, Inf)
                l2norm[i] = norm(err, 2)/sqrt(length(err))
            end

                plt = lineplot(log10.(hlist),
                            log10.(linorm),
                            title = "Gaussian Factor 2D",
                            name = "L_inf norm",
                            xlabel = "log_10(h)",
                            ylabel = "log_10(Error)",
                            xlim = [floor(minimum(log10.(hlist))), ceil(maximum(log10.(hlist)))],
                            ylim = [floor(minimum(log10.([linorm; l2norm]))), ceil(maximum(log10.([linorm; l2norm])))],
                            color = :blue)

            lineplot!(plt, log10.(hlist), log10.(l2norm), name = "mean L2 norm", color = :red)
            println(plt)
            @test reverse(linorm) == sort(linorm) #test for decreasing error with decreasing h
            @test reverse(l2norm) == sort(l2norm)
        end

        @testset "Constant Gradient of Squared Slowness 3D" begin
            hlist = [1/10, 1/20, 1/40, 1/80]#, 1/160, 1/320]
            linorm = zeros(length(hlist))
            l2norm = zeros(length(hlist))
            for i = 1:length(hlist)
                (x0, k2, t_exact) = const_k2_3D(hlist[i])
                (t_pred, ordering) = FEFMM.fefmm(k2, [hlist[i], hlist[i], hlist[i]], x0)
                err = t_exact.-t_pred
                linorm[i] = norm(err, Inf)
                l2norm[i] = norm(err, 2)/sqrt(length(err))
            end

            plt = lineplot(log10.(hlist),
                           log10.(linorm),
                           title = "Constant Gradient of Squared Slowness 3D",
                           name = "L_inf norm",
                           xlabel = "log_10(h)",
                           ylabel = "log_10(Error)",
                           xlim = [floor(minimum(log10.(hlist))), ceil(maximum(log10.(hlist)))],
                           ylim = [floor(minimum(log10.([linorm; l2norm]))), ceil(maximum(log10.([linorm; l2norm])))],
                           color = :blue)

            lineplot!(plt, log10.(hlist), log10.(l2norm), name = "mean L2 norm", color = :red)
            println(plt)
            @test reverse(linorm) == sort(linorm) #test for decreasing error with decreasing h
            @test reverse(l2norm) == sort(l2norm)
        end

        @testset "Constant Gradient of Squared Velocity 3D" begin
            hlist = [1/10, 1/20, 1/40, 1/80]#, 1/160, 1/320]
            linorm = zeros(length(hlist))
            l2norm = zeros(length(hlist))
            for i = 1:length(hlist)
                (x0, k2, t_exact) = const_v2_3D(hlist[i])
                (t_pred, ordering) = FEFMM.fefmm(k2, [hlist[i], hlist[i], hlist[i]], x0)
                err = t_exact.-t_pred
                linorm[i] = norm(err, Inf)
                l2norm[i] = norm(err, 2)/sqrt(length(err))
            end

            plt = lineplot(log10.(hlist),
                           log10.(linorm),
                           title = "Constant Gradient of Squared Velocity 3D",
                           name = "L_inf norm",
                           xlabel = "log_10(h)",
                           ylabel = "log_10(Error)",
                           xlim = [floor(minimum(log10.(hlist))), ceil(maximum(log10.(hlist)))],
                           ylim = [floor(minimum(log10.([linorm; l2norm]))), ceil(maximum(log10.([linorm; l2norm])))],
                           color = :blue)

            lineplot!(plt, log10.(hlist), log10.(l2norm), name = "mean L2 norm", color = :red)
            println(plt)
            @test reverse(linorm) == sort(linorm) #test for decreasing error with decreasing h
            @test reverse(l2norm) == sort(l2norm)
        end

        @testset "Gaussian Factor 3D" begin
            hlist = [1/10, 1/20, 1/40, 1/80]#, 1/160, 1/320]
            linorm = zeros(length(hlist))
            l2norm = zeros(length(hlist))
            for i = 1:length(hlist)
                (x0, k2, t_exact) = gaussian_factor_3D(hlist[i])
                (t_pred, ordering) = FEFMM.fefmm(k2, [hlist[i], hlist[i], hlist[i]], x0)
                err = t_exact.-t_pred
                linorm[i] = norm(err, Inf)
                l2norm[i] = norm(err, 2)/sqrt(length(err))
            end

            plt = lineplot(log10.(hlist),
                           log10.(linorm),
                           title = "Gaussian Factor 3D",
                           name = "L_inf norm",
                           xlabel = "log_10(h)",
                           ylabel = "log_10(Error)",
                           xlim = [floor(minimum(log10.(hlist))), ceil(maximum(log10.(hlist)))],
                           ylim = [floor(minimum(log10.([linorm; l2norm]))), ceil(maximum(log10.([linorm; l2norm])))],
                           color = :blue)

            lineplot!(plt, log10.(hlist), log10.(l2norm), name = "mean L2 norm", color = :red)
            println(plt)
            @test reverse(linorm) == sort(linorm) #test for decreasing error with decreasing h
            @test reverse(l2norm) == sort(l2norm)
        end

    end
end