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

using Markdown
using InteractiveUtils

# ╔═╡ e2d82196-7945-11ec-2d31-ad911cabbfba
using CSV, DataFrames

# ╔═╡ 183ac353-4850-4cc6-be57-ba849fd50739
using DataFramesMeta

# ╔═╡ 403b0807-419c-4804-86b8-aca28751a1a6
using AstroAngles

# ╔═╡ 5337ca99-20c5-4418-8fff-76191d718522
begin
	using CondaPkg
	CondaPkg.add.(("astroquery", "astropy", "pandas")); CondaPkg.resolve()
	using PythonCall
end

# ╔═╡ 8293e156-e76e-41bc-a3cd-07877eb23663
using LinearAlgebra # Defines `norm` for us

# ╔═╡ b3df4b34-9287-4de0-bb26-e717ddad067d
using Latexify

# ╔═╡ 27c0f880-b775-458f-9092-a0d4d743d9ab
begin
	using PlutoUI # Handy for making widgets and such, like a table of contents
	import MarkdownLiteral: @mdx # Markdown + HTML = superpowers
end

# ╔═╡ 80664df1-a135-44f9-b83f-41029a503c77
@mdx """
# Observations at LCO

In this notebook, we are going to take an *ugly* SMF file and convert it into a gorgeous ``\\LaTeX`` table, complete with star coordinates pulled down from VizieR via [astroquery](https://astroquery.readthedocs.io/en/latest/) ✨
"""

# ╔═╡ e1ae732a-22de-4f5d-8216-064dd3e3a686
TableOfContents()

# ╔═╡ 46bee12a-81ca-43ea-932f-8cfa22ad1e9a
@mdx """
## Introduction 🤚

Las Campanas Observatory ([LCO](http://www.lco.cl/)) is home to several world-class telescopes that are used in a wide range of astronomy fields. In particular, The 6.5 meter twin [Magellan telescopes](https://obs.carnegiescience.edu/Magellan) provide essential facilities for observing exoplanet atmopsheres with their [IMACS and LDSS3C spectrographs](http://www.lco.cl/instruments/).

SMF files are essentially just text files with a list of target names and positions that are sent to the telescope to let the observatory know how to cut the mask for observations. They also include additonal information for the slits, but this will not be the focus for this notebook. Now this is all well and good, except for the fact that these SMF files tend to look something like this:
"""

# ╔═╡ 4a079964-5ac0-47a0-9ff0-c294df62772e
"NAME wasp50s
OBSERVER Cookie Monster
TITLE wasp50 science mask
MADE 2013-08-09  MaskGen  2.27.92
TELESCOPE 71161.1 0.80906 3 Magellan
INSTRUMENT IMACS_sc
DISPERSER  IMACS_grism_300  1  17.50000
POSITION 0 02:55:00.000 -10:57:00.00 2000.0  0.0 0.0 0.0
!Constellation (Eri) Eridanus
!# No rotator warnings above horizon.
!.OC wasp50s 02:55:00.000 -10:57:00.00 2000.0 0.0 0.0 -136.15 OFF  02:55:11.3 -10:41:29 2000.0 02:54:41.9 -11:12:17 2000.0
WAVELENGTH 6749.3
TEMPERATURE 12.0
EPOCH 2013.70211
WLIMIT 4000.1 7900.1
DREF  0
HANGLE    0
DATE  56548
SLIT WASP50    02:54:45.137 -10:53:53.03   75.570   64.529  3.450  3.795  3.795    0.00 # Vmag= 11.721
SLIT c06       02:54:49.392 -10:51:54.94   53.950  105.319  3.450  3.795  3.795    0.00 # Vmag= 11.585
SLIT c13       02:54:23.502 -11:02:28.19  185.851 -113.558  3.450  3.795  3.795    0.00 # Vmag= 13.075
SLIT c15       02:54:32.141 -10:54:36.94  141.733   49.389  3.450  3.795  3.795    0.00 # Vmag= 13.252
SLIT c18       02:54:56.565 -11:02:27.36   17.461 -113.018  3.450  3.795  3.795    0.00 # Vmag= 13.568
SLIT c20       02:54:14.107 -11:03:06.82  233.973 -127.107  3.450  3.795  3.795    0.00 # Vmag= 13.662
SLIT c21       02:54:52.547 -10:51:26.96   37.908  114.984  3.450  3.795  3.795    0.00 # Vmag= 13.668
SLIT c23       02:55:16.318 -10:53:50.69  -82.969   65.341  3.450  3.795  3.795    0.00 # Vmag= 13.708
SLIT c28       02:55:35.977 -10:56:20.97 -183.127   13.446  3.450  3.795  3.795    0.00 # Vmag= 13.885
HOLE ref015    02:54:53.542 -11:08:45.55   32.898 -244.204  1.725 1  0.862  0.862    0.00 # Vmag= 14.897
HOLE ref031    02:55:46.106 -10:51:22.87 -235.193  116.663  1.725 1  0.862  0.862    0.00 # Vmag= 15.430
HOLE ref036    02:54:05.105 -10:56:02.52  280.104   19.807  1.725 1  0.862  0.862    0.00 # Vmag= 15.490
HOLE ref042    02:54:40.376 -11:01:26.73   99.779  -92.124  1.725 1  0.862  0.862    0.00 # Vmag= 15.714
HOLE ref044    02:54:06.394 -10:50:53.83  273.783  126.843  1.725 1  0.862  0.862    0.00 # Vmag= 15.728
HOLE ref046    02:54:53.749 -10:55:28.22   31.764   31.665  1.725 1  0.862  0.862    0.00 # Vmag= 15.762
HOLE ref053    02:54:56.291 -10:45:54.96   18.908  230.088  1.725 1  0.862  0.862    0.00 # Vmag= 15.865
HOLE ref054    02:55:19.594 -10:49:31.69  -99.775  154.931  1.725 1  0.862  0.862    0.00 # Vmag= 15.894
HOLE ref058    02:55:44.049 -11:00:33.94 -224.424  -74.096  1.725 1  0.862  0.862    0.00 # Vmag= 15.997
HOLE ref068    02:55:08.761 -11:01:26.28  -44.527  -91.920  1.725 1  0.862  0.862    0.00 # Vmag= 16.266
HOLE ref071    02:54:25.846 -11:09:07.73  174.283 -252.376  1.725 1  0.862  0.862    0.00 # Vmag= 16.322
HOLE ref072    02:54:42.547 -10:52:41.78   88.772   89.150  1.725 1  0.862  0.862    0.00 # Vmag= 16.351
HOLE ref076    02:54:18.818 -10:46:32.58  210.387  217.464  1.725 1  0.862  0.862    0.00 # Vmag= 16.407
HOLE ref080    02:55:46.477 -10:59:38.07 -236.851  -54.783  1.725 1  0.862  0.862    0.00 # Vmag= 16.567
HOLE ref085    02:54:19.177 -10:58:22.60  207.882  -28.623  1.725 1  0.862  0.862    0.00 # Vmag= 17.101
<and then some random long string>
" |> Text

# ╔═╡ 3dcf92a9-21c2-4cfd-9089-e6fb501c7caf
@mdx """
!!! note
	More info about what each column means can be found [here](https://code.obs.carnegiescience.edu/maskgen/v214/smdfile.txt/view).
"""

# ╔═╡ 3941e471-6378-4f6b-a117-81aa08eff0c7
@mdx """
We need a way to ingest this highly unstructured file and extract some useful information from it. For this notebook, the final product will be this summary report that we can paste into a ``\\LaTeX`` document:
"""

# ╔═╡ af553292-ae3e-4966-9746-79de4b86a7e1
@mdx """
## Game plan 🏈
To generate the above, we are going to do the following:

1. Extract only the rows from the SMF file that refer to the slits (the rest just refer to other objects that we don't care about for this analysis)

2. Use the RA and Dec coordinates to find the corresponding star in each row using `astroquery`

3. Save each star's ID and a few other important quantities in a final report table

Alright, let's get started.
"""

# ╔═╡ 769c4b03-4ac0-4a25-9e7d-21825261c3dd
@mdx """
## 1. Read the SMF 🔎

We are first going to use [`CSV.jl`](https://github.com/JuliaData/CSV.jl) to read the SMF file into a `DataFrame` (provided by [`DataFrames.jl`](https://github.com/JuliaData/DataFrames.jl)):
"""

# ╔═╡ 89136b56-99f6-49bf-a823-6447fee208fc
df_SMF = CSV.read(download("https://github.com/icweaver/WASP-50b/raw/main/notebooks/data/raw_data/wasp50s.SMF"), DataFrame;
	header = false, # Wish there was one
	skipto = 19, # This is the row where the good stuff starts
	footerskip = 1, # Skip that random string at the bottom
	delim = ' ', # Spaces, spaces everywhere
	ignorerepeated = true, # Treat repeated delimiters as a single one
	comment = "HOLE", # Ignore all lines starting with this
	types = String, # Treat every column as a `String`
)

# ╔═╡ 16d610c6-6b7e-46a2-bd54-487392248302
@mdx """
Next let's just select the columns that we will need for the report. We use the `@select` macro from the [`DataFramesMeta.jl`](https://github.com/JuliaData/DataFramesMeta.jl) package to accomplish this in a convenient way:
"""

# ╔═╡ 8ceb1ae4-eecd-4e05-94cf-6e84749d6727
df = @select df_SMF begin
	:name = :Column2
	:_RAJ2000 = @. string(hms2deg(:Column3)) * "d" # astroquery needs columns
	:_DEJ2000 = @. string(dms2deg(:Column4)) * "d" # in this format
end

# ╔═╡ ed944c35-09a6-4b69-8ccc-20e238fa488d
@mdx """
!!! tip
	There is A LOT more that we can do with these kinds of macros. A nice overview from the creator of `DataFrames.jl` can be found [here](https://bkamins.github.io/julialang/2021/11/19/dfm.html).
"""

# ╔═╡ 07d51156-9b82-4b9b-b6a6-a8bf09966880
@mdx """
We also used `hms2deg` and `dms2deg` from [`AstroAngles.jl`](https://github.com/JuliaAstro/AstroAngles.jl) to convert the RA and Dec into a format that can automatically be read by `astroquery`, which we will do next.
"""

# ╔═╡ 21e927a2-ae8b-4623-ba17-e98ca3c04f60
@mdx """
## 2. Query the star catalog 📖

Here comes the fun part. We are going to feed this table into the Python package `astroquery`, which will query an all-sky star catalog for us and return the corresponding table of results. Let's start by setting up a minimal Python environment and download `astroquery` into it:
"""

# ╔═╡ 75a654c0-f15d-4cf0-b8bd-8f7a10de0e59
@mdx """
We used [CondaPkg.jl](https://github.com/cjdoris/CondaPkg.jl) to install a light-weight Python environment (via [MicroMamba.jl](https://github.com/cjdoris/MicroMamba.jl)) into a temporary folder, along with `astroquery` and a few other Python packages that will be needed for compatibility reasons later. We next use `PythonCall.jl` to allow us to seamlessly interact with Python objects from Julia.

Now that that's all set up, let's import some things!
"""

# ╔═╡ 5926658d-bcb1-477a-bebc-9414eaa3813c
@py begin
	import astroquery.vizier as vz
	import astropy.table as ap
end

# ╔═╡ 9a4f6266-f761-4eca-8eba-dd8f611a56bf
@mdx """
For SOME reason, the only flavor of table that `astroquery` accepts is an astropy `Table`, so we are going to do the following roundabout conversions to make the 🐍 happy:

```
DataFrame (julia) ---> DataFrame (pandas) ---> Table (astropy)
```

We accomplish the first step with `PythonCall.jl`'s `pytable` function.
"""

# ╔═╡ f9d21956-347b-48fc-abe2-04d733c92668
results_astropy = vz.Vizier.query_region(
	df |> pytable |> ap.Table.from_pandas,
	radius = "1 arcsec",
	catalog = "UCAC4",
	frame = "icrs", # I think this is default
)[0]

# ╔═╡ 9c710f1f-4029-4678-8ed8-654ffe107a59
@mdx """
!!! note
	We also could have done:

	```julia
	ap.Table.from_pandas(pytable(df))
	```

	for the conversion, but pipes `|>` can help make multi-step workflows more transparent.
"""

# ╔═╡ 905127e7-16ae-4d8d-8445-7197cc5dee2a
@mdx """
Alright, we have a table of results now! These are the list of objects (from the [UCAC4 catalog](https://vizier.u-strasbg.fr/viz-bin/VizieR?-source=I/322)) that, row-by-row from `df`, are within 1 arcsecond of the RA and Decs specified. We set this radius to be small enough so that only one target from the catalog (ideally the star we are searching for) matches the coords given.
"""

# ╔═╡ 1f2afbc8-5e04-4f18-b047-90775443d572
@mdx """
## 3. Generate the report ✏

Now that we have our list of targets, let's extract a few key columns and compute a color metric, for example. First let's convert the returned astropy `Table` back to a native Julia `DataFrame` to make things a bit easier to work with:
"""

# ╔═╡ 15eef4f3-53ac-45f0-953c-0e46e14f3e46
results = results_astropy.to_pandas() |> PyTable

# ╔═╡ ea756fff-853a-41b8-902b-a495cee1fd43
@mdx """
We generated a [`PyPandasDataFrame`](https://cjdoris.github.io/PythonCall.jl/stable/compat/#Tabular-data-and-Pandas) via `PythonCall.jl`'s `PyTable` function, which is basically a generic table type in Julia. We'd like to access the full feature set of `DataFrames.jl`, so we will construct one directly from `results`:
"""

# ╔═╡ 37d2bc2e-ee75-42a9-ba4f-df15c99caf48
_df_paper = DataFrame(
	Star = df.name,
	UCAC4_ID = results.UCAC4,
	# RAJ2000 = results.RAJ2000,
	# DecJ200 = results.DEJ2000,
	B = results.Bmag,
	V = results.Vmag,
	J = results.Jmag,
	K = results.Kmag,
)

# ╔═╡ d6f7f9eb-5776-4419-98e0-f18cd4cee12f
@mdx """
Now that we have our color magnitudes for each target, we compute a color-space metric ``D``:

```math
D \\equiv \\sqrt{
	  \\left[(B-V)_c-(B-V)_t\\right]^{2}

	+ \\left[(J-K)_c-(J-K)_t\\right]^{2}
}
```

relative to our main star, "WASP-50":
"""

# ╔═╡ 3a2891dc-36bf-4e1f-913a-4ba53b47eac6
"""
	compute_D(B, V, J, K, [ΔBVₜ=0.0, ΔJKₜ=0.0])

* `B, V, J, K`: comparison star magnitudes:

```math
	B_c, V_c, J_c, K_c
```

* `ΔBVₜ`, `ΔJKₜ`: target star delta mags (defaults to zero for each):

```math
(B-V)_{t}, (J-K)_{t}
```
"""
compute_D(B, V, J, K; ΔBVₜ=0.0, ΔJKₜ=0.0) = norm(((B-V)-ΔBVₜ, (J-K)-ΔJKₜ))

# ╔═╡ 96035871-6bcb-4703-9507-3365f92a93d9
@mdx """
and add it to our table:
"""

# ╔═╡ f3326ab8-21da-403b-b290-7694cefb63fe
df_paper = @chain _df_paper begin
	@aside ΔBVₜ, ΔJKₜ = _.B[1] - _.V[1], _.J[1] - _.K[1]
	@transform :D = compute_D.(:B, :V, :J, :K; ΔBVₜ, ΔJKₜ)
end

# ╔═╡ 8fbc74d9-c993-4757-865f-b7ec3d84495c
latexify(df_paper, latex=false, fmt="%.2f")

# ╔═╡ c9b37015-057e-4838-b510-5762819b2462
@mdx """
!!! note
	We used the `aside` macro from `DataFramesMeta.jl` to compute intermediate values needed for ``D`` before adding it to our table, which was all done in a chain of actions with `@chain` from the [`Chain.jl`](https://github.com/jkrumbiegel/Chain.jl) package that is automatically included (exported) with `DataFramesMeta.jl`.

	We could also have used bangs (!) to mutate tables in-place and save memory, e.g.:

	```julia
	@chain df_paper begin
		@aside ΔBVₜ, ΔJKₜ = _.B[1] - _.V[1], _.J[1] - _.K[1]
		@transform! :D = compute_D.(:B, :V, :J, :K; ΔBVₜ, ΔJKₜ)
	end
	```

	but for clarity, we are just creating a new table each time in this notebook.
"""

# ╔═╡ 06d8b523-4a05-43ec-973b-8fe55fa613a1
@mdx """
Ok, now we can display `df_paper` in a ``LaTeX`` context with [`Latexify.jl`](https://github.com/korsbo/Latexify.jl):
"""

# ╔═╡ 4f3afae4-7b3c-42ed-81a0-936e1348700a
t = latexify(df_paper, env=:tabular, latex=false, fmt="%.2f") |> Text

# ╔═╡ 875ca00d-0c0c-403e-bbf1-0d2b70a95107
t

# ╔═╡ 52100a8e-3adf-4337-9985-645a10d9e8cd
@mdx """
!!! note
	We used `latex=false` to remove the dollar signs around each number, because this will be used in a [deluxetable](https://journals.aas.org/aastexguide/#preamble_deluxetable) environment that can handle automatic math mode columns.

	To display the nicely rendered ``\\LaTeX`` table at the beginning of the document, we just used the default Markdown backend:

	```julia
	latexify(df_paper, latex=false, fmt="%.2f")
	```
"""

# ╔═╡ 34d1417b-2d56-4de0-800e-e65d5f9817fb
@mdx """
We now have a copy-and-pasteable table for a paper that we are totally not ignoring right now ✔
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
AstroAngles = "5c4adb95-c1fc-4c53-b4ea-2a94080c53d2"
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
CondaPkg = "992eb4ea-22a4-4c89-a5bb-47a3300528ab"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
DataFramesMeta = "1313f7d8-7da2-5740-9ea0-a2ca25f37964"
Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
MarkdownLiteral = "736d6165-7244-6769-4267-6b50796e6954"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
PythonCall = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"

[compat]
AstroAngles = "~0.1.3"
CSV = "~0.10.2"
CondaPkg = "~0.2.6"
DataFrames = "~1.3.2"
DataFramesMeta = "~0.10.0"
Latexify = "~0.15.9"
MarkdownLiteral = "~0.1.1"
PlutoUI = "~0.7.32"
PythonCall = "~0.5.1"
"""

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

julia_version = "1.7.1"
manifest_format = "2.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

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

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

[[deps.AstroAngles]]
git-tree-sha1 = "41621fa5ed5f7614b75eea8e0b3cfd967b284c87"
uuid = "5c4adb95-c1fc-4c53-b4ea-2a94080c53d2"
version = "0.1.3"

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

[[deps.CSV]]
deps = ["CodecZlib", "Dates", "FilePathsBase", "InlineStrings", "Mmap", "Parsers", "PooledArrays", "SentinelArrays", "Tables", "Unicode", "WeakRefStrings"]
git-tree-sha1 = "9519274b50500b8029973d241d32cfbf0b127d97"
uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
version = "0.10.2"

[[deps.Chain]]
git-tree-sha1 = "339237319ef4712e6e5df7758d0bccddf5c237d9"
uuid = "8be319e6-bccf-4806-a6f7-6fae938471bc"
version = "0.4.10"

[[deps.CodecZlib]]
deps = ["TranscodingStreams", "Zlib_jll"]
git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da"
uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
version = "0.7.0"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.CommonMark]]
deps = ["Crayons", "JSON", "URIs"]
git-tree-sha1 = "4aff51293dbdbd268df314827b7f409ea57f5b70"
uuid = "a80b9123-70ca-4bc0-993e-6e3bcb318db6"
version = "0.8.5"

[[deps.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 = "44c37b4636bc54afac5c574d2d02b625349d6582"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.41.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.CondaPkg]]
deps = ["JSON3", "Markdown", "MicroMamba", "Pkg", "TOML"]
git-tree-sha1 = "49589b0a9284fd8ea84c64a42e96739d1fe06751"
uuid = "992eb4ea-22a4-4c89-a5bb-47a3300528ab"
version = "0.2.6"

[[deps.Crayons]]
git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.1.1"

[[deps.DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

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

[[deps.DataFramesMeta]]
deps = ["Chain", "DataFrames", "MacroTools", "OrderedCollections", "Reexport"]
git-tree-sha1 = "ab4768d2cc6ab000cd0cec78e8e1ea6b03c7c3e2"
uuid = "1313f7d8-7da2-5740-9ea0-a2ca25f37964"
version = "0.10.0"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "3daef5523dd2e769dad2365274f760ff5f282c7d"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.11"

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

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

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

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

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

[[deps.FilePathsBase]]
deps = ["Compat", "Dates", "Mmap", "Printf", "Test", "UUIDs"]
git-tree-sha1 = "04d13bfa8ef11720c24e4d840c0033d145537df7"
uuid = "48062228-2e41-5def-b9a4-89aafe57970f"
version = "0.9.17"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

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

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

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.InlineStrings]]
deps = ["Parsers"]
git-tree-sha1 = "8d70835a3759cdd75881426fced1508bb7b7e1b6"
uuid = "842dd82b-1e85-43dc-bf29-5d0ee9dffc48"
version = "1.1.1"

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

[[deps.InvertedIndices]]
git-tree-sha1 = "bee5f1ef5bf65df56bdd2e40447590b272a5471f"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.1.0"

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

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

[[deps.JSON3]]
deps = ["Dates", "Mmap", "Parsers", "StructTypes", "UUIDs"]
git-tree-sha1 = "7d58534ffb62cd947950b3aa9b993e63307a6125"
uuid = "0f8b85d8-7281-11e9-16c2-39a750bddbf1"
version = "1.9.2"

[[deps.LaTeXStrings]]
git-tree-sha1 = "f2355693d6778a178ade15952b7ac47a4ff97996"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.3.0"

[[deps.Latexify]]
deps = ["Formatting", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "Printf", "Requires"]
git-tree-sha1 = "a8f4f279b6fa3c3c4f1adadd78a621b13a506bce"
uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
version = "0.15.9"

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

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

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

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

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

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

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

[[deps.MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.9"

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

[[deps.MarkdownLiteral]]
deps = ["CommonMark", "HypertextLiteral"]
git-tree-sha1 = "0d3fa2dd374934b62ee16a4721fe68c418b92899"
uuid = "736d6165-7244-6769-4267-6b50796e6954"
version = "0.1.1"

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

[[deps.MicroMamba]]
deps = ["Pkg", "Scratch"]
git-tree-sha1 = "e61632c723f4d0fafa5419a506f08baafe973391"
uuid = "0b3b1443-0f03-428d-bdfb-f27f9c1191ea"
version = "0.1.4"

[[deps.Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

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

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

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

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

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

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "92f91ba9e5941fc781fecf5494ac1da87bdac775"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.0"

[[deps.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"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "ae6145ca68947569058866e443df69587acc1806"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.32"

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

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

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

[[deps.PythonCall]]
deps = ["CondaPkg", "Dates", "Libdl", "MacroTools", "Markdown", "Pkg", "Requires", "Serialization", "Tables", "UnsafePointers"]
git-tree-sha1 = "402a4d004a9e7c9656758f479f456205f7f6a3ce"
uuid = "6099a3de-0909-46bc-b1f4-468b9a2dfc0d"
version = "0.5.1"

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

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

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.3.0"

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

[[deps.Scratch]]
deps = ["Dates"]
git-tree-sha1 = "0b4b7f1393cff97c33891da2a0bf69c6ed241fda"
uuid = "6c6a2e73-6563-6170-7368-637461726353"
version = "1.1.0"

[[deps.SentinelArrays]]
deps = ["Dates", "Random"]
git-tree-sha1 = "15dfe6b103c2a993be24404124b8791a09460983"
uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
version = "1.3.11"

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

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

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

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

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

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

[[deps.StructTypes]]
deps = ["Dates", "UUIDs"]
git-tree-sha1 = "d24a825a95a6d98c385001212dc9020d609f2d4f"
uuid = "856f2bd8-1eba-4b0a-8007-ebc267875bd4"
version = "1.8.1"

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

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

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

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

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

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

[[deps.URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

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

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

[[deps.UnsafePointers]]
git-tree-sha1 = "c81331b3b2e60a982be57c046ec91f599ede674a"
uuid = "e17b2a0c-0bdf-430a-bd0c-3a23cae4ff39"
version = "1.0.0"

[[deps.WeakRefStrings]]
deps = ["DataAPI", "InlineStrings", "Parsers"]
git-tree-sha1 = "c69f9da3ff2f4f02e811c3323c22e5dfcb584cfa"
uuid = "ea10d353-3f73-51f8-a26c-33c1cb351aa5"
version = "1.4.1"

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

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

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

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

# ╔═╡ Cell order:
# ╟─80664df1-a135-44f9-b83f-41029a503c77
# ╠═e1ae732a-22de-4f5d-8216-064dd3e3a686
# ╟─46bee12a-81ca-43ea-932f-8cfa22ad1e9a
# ╟─4a079964-5ac0-47a0-9ff0-c294df62772e
# ╟─3dcf92a9-21c2-4cfd-9089-e6fb501c7caf
# ╟─3941e471-6378-4f6b-a117-81aa08eff0c7
# ╟─8fbc74d9-c993-4757-865f-b7ec3d84495c
# ╟─875ca00d-0c0c-403e-bbf1-0d2b70a95107
# ╟─af553292-ae3e-4966-9746-79de4b86a7e1
# ╟─769c4b03-4ac0-4a25-9e7d-21825261c3dd
# ╠═e2d82196-7945-11ec-2d31-ad911cabbfba
# ╠═89136b56-99f6-49bf-a823-6447fee208fc
# ╟─16d610c6-6b7e-46a2-bd54-487392248302
# ╠═183ac353-4850-4cc6-be57-ba849fd50739
# ╠═403b0807-419c-4804-86b8-aca28751a1a6
# ╠═8ceb1ae4-eecd-4e05-94cf-6e84749d6727
# ╟─ed944c35-09a6-4b69-8ccc-20e238fa488d
# ╟─07d51156-9b82-4b9b-b6a6-a8bf09966880
# ╟─21e927a2-ae8b-4623-ba17-e98ca3c04f60
# ╠═5337ca99-20c5-4418-8fff-76191d718522
# ╟─75a654c0-f15d-4cf0-b8bd-8f7a10de0e59
# ╠═5926658d-bcb1-477a-bebc-9414eaa3813c
# ╟─9a4f6266-f761-4eca-8eba-dd8f611a56bf
# ╠═f9d21956-347b-48fc-abe2-04d733c92668
# ╟─9c710f1f-4029-4678-8ed8-654ffe107a59
# ╟─905127e7-16ae-4d8d-8445-7197cc5dee2a
# ╟─1f2afbc8-5e04-4f18-b047-90775443d572
# ╠═15eef4f3-53ac-45f0-953c-0e46e14f3e46
# ╟─ea756fff-853a-41b8-902b-a495cee1fd43
# ╠═37d2bc2e-ee75-42a9-ba4f-df15c99caf48
# ╟─d6f7f9eb-5776-4419-98e0-f18cd4cee12f
# ╠═8293e156-e76e-41bc-a3cd-07877eb23663
# ╠═3a2891dc-36bf-4e1f-913a-4ba53b47eac6
# ╟─96035871-6bcb-4703-9507-3365f92a93d9
# ╠═f3326ab8-21da-403b-b290-7694cefb63fe
# ╟─c9b37015-057e-4838-b510-5762819b2462
# ╟─06d8b523-4a05-43ec-973b-8fe55fa613a1
# ╠═b3df4b34-9287-4de0-bb26-e717ddad067d
# ╠═4f3afae4-7b3c-42ed-81a0-936e1348700a
# ╟─52100a8e-3adf-4337-9985-645a10d9e8cd
# ╟─34d1417b-2d56-4de0-800e-e65d5f9817fb
# ╠═27c0f880-b775-458f-9092-a0d4d743d9ab
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using MultivariateStats
using LinearAlgebra
using Test
import Statistics: mean, cov
import Random

@testset "Whitening" begin

    Random.seed!(34568)

    ## data

    d = 3
    n = 5

    X = rand(d, n)
    mval = vec(mean(X, dims=2))
    C = cov(X, dims=2)
    C0 = copy(C)

    emax = maximum(eigvals(C))
    rc = 0.01
    Cr = C + (emax * rc) * Matrix(I, d, d)

    ## cov_whitening

    cf = cholesky(Hermitian(C, :U))
    W = cov_whitening(cf)
    @test istriu(W)
    @test W'C * W ≈ Matrix(I, d, d)

    cf = cholesky(Hermitian(C, :L))
    W = cov_whitening(cf)
    @test istriu(W)
    @test W'C * W ≈ Matrix(I, d, d)

    W = cov_whitening(C)
    @test C == C0
    @test istriu(W)
    @test W'C * W ≈ Matrix(I, d, d)

    W = cov_whitening(C, 0.0)
    @test C == C0
    @test istriu(W)
    @test W'C * W ≈ Matrix(I, d, d)

    W = cov_whitening(C, rc)
    @test istriu(W)
    @test W'Cr * W ≈ Matrix(I, d, d)

    # Whitening from data

    f = fit(Whitening, X)
    W = f.W
    @test isa(f, Whitening{Float64})
    @test mean(f) === f.mean
    @test istriu(W)
    @test W'C * W ≈ Matrix(I, d, d)
    @test transform(f, X) ≈ W' * (X .- f.mean)

    f = fit(Whitening, X; regcoef=rc)
    W = f.W
    @test W'Cr * W ≈ Matrix(I, d, d)

    f = fit(Whitening, X; mean=mval)
    W = f.W
    @test W'C * W ≈ Matrix(I, d, d)

    f = fit(Whitening, X; mean=0)
    Cx = (X * X') / (n - 1)
    W = f.W
    @test W'Cx * W ≈ Matrix(I, d, d)

    # invsqrtm

    R = invsqrtm(C)
    @test C == C0
    @test R ≈ inv(sqrt(C))

end

ext: jl


text:
using Documenter
using MarketData

makedocs(
    format = :html,
    sitename = "MarketData.jl",
    pages = [
        "index.md",
        "getting_started.md",
        "test_data.md",
        "company_financial_series.md",
        "downloads.md",
    ]
)

deploydocs(
    repo = "github.com/JuliaQuant/MarketData.jl.git",
    julia  = "1.0",
    latest = "master",
    target = "build",
    deps = nothing,  # we use the `format = :html`, without `mkdocs`
    make = nothing,  # we use the `format = :html`, without `mkdocs`
)


text:
# GRAPE Amplitude Limiting via penalty
# given a set of limit control amplitudes, calculate a penalty (order ~1) when the limit is reached by the current control amplitudes. As the limit is approached, given the option for a gradually increasing penalty
#
# Questions:
# 1. given u and ulim, how do we calculate the penalty?
# 2. user function?
# 3. parameters?
# 4. do we average the penalty over all control amplitudes, or just sum them up?
#
# Discussion:
# the fidelity error used as objective function has range between 0 and 1, therefore, it makes sense for the penalty to be of order 1. (right?) We do not want the penatly to completely overtake the fidelity error, especially if we do have an option for smaller penalties near the limit.

# After much thinking, the best penalty function is:
#penalty(u,ulim=1,p=1) = p*u^10/(1000000*(ulim^10-u^10))
#penaltyprime(u,ulim=1,p=1) = 10p*abs(u)^9*ulim^10/(1000000*(u^10-ulim^10)^2)
# The parameter p is chosend such that the penalty ~ 1E-9 for u/ulim = 0.5
#penalty(u,ulim=1,p=1) = 1/(1+exp(-1000p*(u-ulim)/ulim))
#penaltyprime(u,ulim=1,p=1) = (x=penalty(u,ulim,p); 1000x*(1-x))

penalty(u,ulim=1,p=1) = (abs(u)/ulim)^100p
penaltyprime(u,ulim=1,p=1) = 100p*(abs(u)/ulim)^(100p-1)

struct AmplitudeLimiter{T<:Union{Array,Number}} <: PenaltyFunction
    limit::T
    p::Float64
end

function objective(L::AmplitudeLimiter{T},u) where {T<:Array}
    N = length(u)
    N == length(L.limit) || throw(DimensionMismatch("Control and limit amplitudes do not have the same size!"))
    res = 0.0
    for i = 1:N
        @inbounds res += penalty(u[i],L.limit[i],L.p)
    end
    return res/N
end

function objective(L::AmplitudeLimiter{T},u) where {T<:Number}
    return sum(u->penalty(u,L.limit,L.p),u)/length(u)
end

function gradient!(L::AmplitudeLimiter{T},fp,u) where {T<:Array}
    N = length(u)
    N == length(L.limit) || throw(DimensionMismatch("Control and limit amplitudes do not have the same size!"))
    for i = 1:N
        @inbounds fp[i] += penaltyprime(u[i],L.limit[i],L.p)/N
    end
    return fp
end

function gradient!(L::AmplitudeLimiter{T},fp,u) where {T<:Number}
    N = length(u)
    for i = 1:N
        @inbounds fp[i] += penaltyprime(u[i],L.limit,L.p)/N
    end
    return fp
end


text:
module TuringMachine
export load_input, load_program, set_up, simulate

using DataStructures, DelimitedFiles

"""
    load_input(file_name)

Takes a text file containing a continuous string of `n` integers (no spaces) and returns an `n`-element
array, where each element is an integer. The array is ordered as the input string.

Input
- `file_name`::string : path of text file containing an input string of `n` integers (no spaces)

Output
- ::`n`-element Array{Char, 1} : output array of integers, ordered as the input string
"""
function load_input(file_name)
    file = readdlm(file_name, String)
    return collect(file[])
end

"""
    load_program(file_name)

Takes a text file containing a Turing program. Please see docs for how to write such a program.

Input
- `file_name`::string : path of text file containing the Turing program to simulate

Output
- `program`::Dict{Any,Any} with `m` entries: dictionary with keys of [states,read_cells] and values of [next_state,write_cell,movement]
"""
function load_program(file_name)
    program_raw = readdlm(file_name)
    rows, cols = size(program_raw)
    program = Dict{Any,Any}()
    init_state = []
    halt_state = []
    for i = 1:rows
        if program_raw[i,1] == "//"
            continue
        else
            if program_raw[i,1] == "init:"
                push!(init_state, program_raw[i,2])
            elseif program_raw[i,1] == "halt:"
                push!(halt_state, program_raw[i,2])
            else
                card = split(program_raw[i,1], r",")
                state_read = (String(card[1]),String(card[2]))
                instruction = (String(card[3]),String(card[4]),String(card[5]))
                program[state_read] = instruction
            end
        end
    end
    return program, init_state[], halt_state[]
end

"""
    set_up(init_state, input)

Performs initial set-up for the simulation of single tape. The `length(input)` must be equal of greater
than 1. The input will be placed on tape to the right of the head.

Input
- `input`::`n`-element Array{Char, 1} : initial data to insert to the right of the head
- `init_state`::String : string label for initial/starting state of the Turing machine

Outputs
- `tape_left`::Deque{Int} : initial values for the stack to the left of the head
- `tape_right`::Deque{Int} : initial values for the stack to the right of the head
"""
function set_up(init_state, input)
    # set up two stacks to simulate a single tape
    tape_left = Deque{Int}() # simulates tape lying to the left of the head
    tape_right = Deque{Int}() # simulates tape lying to the right of the head
    for i = 1:length(input)
        push!(tape_right, parse(Int, input[i]))
    end
    state = init_state
    return state, tape_left, tape_right
end

"""
    read_tape(tape_right)

Reads a tape and returns the current value.

Input
- `tape_right`::Deque{Int}() : values for the stack to the right of the head

Output
- `read_cell`::String : value for the current cell being read by the head
"""
function read_tape(tape_right)
    read_cell = front(tape_right)
    return read_cell
end

"""
    get_instruction(state, read_cell, program)

Inputs
- `state`::String : curent state of the Turing machine
- `read_cell`::String : value for the current cell being read by the head
- `program`::Dict{Any,Any} with `m` entries: dictionary with keys of [states,read_cells] and values of [next_state,write_cell,movement]

Outputs
- `next_state`::String : next state for the Turing machine to move to
- `write_cell`::String : value for the head to write on the current cell
- `movement`::String : movement instruction for the head (left="<", right=">", any other symbol doesn't move)
"""
function get_instruction(state, read_cell, program)
    instruction = program[(string(state), string(read_cell))]
    next_state = instruction[1]
    write_cell = instruction[2]
    movement = instruction[3]
    return next_state, write_cell, movement
end

"""
    write_move!(movement, write_cell, tape_left, tape_right)

Performs the operations of writing at the current position then moving the head left or right.

Inputs
- `movement`::String : movement instruction for the head (left="<", right=">", any other symbol doesn't move)
- `write_cell`::String : value for the head to write on the current cell
- `tape_left`::Deque{Int} : values for the stack to the left of the head
- `tape_right`::Deque{Int} : values for the stack to the right of the head

Outputs
- `tape_left`::Deque{Int} : updated values for the stack to the left of the head (after writing/movement)
- `tape_right`::Deque{Int} : updated values for the stack to the right of the head (after writing/movement)
"""
function write_move!(movement, write_cell, tape_left, tape_right)
    if isempty(tape_right)
        return tape_left, tape_right
    end
    popfirst!(tape_right)
    pushfirst!(tape_right, parse(Int, write_cell))
    if movement == "<"
        if isempty(tape_left)
            pushfirst!(tape_right, "_")
        else
            carry = pop!(tape_left)
            pushfirst!(tape_right, carry)
        end
    elseif movement == ">"
        if isempty(tape_right)
            push!(tape_left, "_")
        else
            carry = popfirst!(tape_right)
            push!(tape_left, carry)
        end
    end
    return tape_left, tape_right
end

"""
    simulate(state, program, tape_left, tape_right)

Simulates one step of the Turing machine: reading, writing, and moving the head.

Inputs
- `state`::String : curent state of the Turing machine
- `program`::Dict{Any,Any} with `m` entries: dictionary with keys of [states,read_cells] and values of [next_state,write_cell,movement]
- `tape_left`::Deque{Int} : values for the stack to the left of the head
- `tape_right`::Deque{Int} : values for the stack to the right of the head

Outputs
- `state`::String : updated state of the Turing machine after one simulation step
- `tape_left`::Deque{Int} : updated values for the stack to the left of the head after one simulation step
- `tape_right`::Deque{Int} : updated values for the stack to the right of the head after one simulation step
"""
function simulate(state, program, tape_left, tape_right)
    if isempty(tape_right)
        read_cell = "_"
    else
        read_cell = read_tape(tape_right)
    end
    state, write_cell, movement = get_instruction(state, read_cell, program)
    tape_left, tape_right = write_move!(movement, write_cell, tape_left, tape_right)
    return state, tape_left, tape_right
end

end # module


text:
module Identifiers

export Identifier, Sedol, Cusip, Isin, Ticker, Figi, Index

abstract type Identifier <: AbstractString end

"""
    Sedol(x::String)

Create a Sedol Identifier.

# Example
```jldoctest; setup = :(using FinancialSymbology)
julia> sedol_id = Sedol("B0YQ5W0")
"B0YQ5W0"
```
"""
struct Sedol <: Identifier
    s::AbstractString
end

"""
    Cusip(x::String)

Create a Cusip Identifier.

# Example
```jldoctest; setup = :(using FinancialSymbology)
julia> cusip_id = Cusip("037833100")
"037833100"
```
"""
struct Cusip <: Identifier
    s::AbstractString
end

"""
    Isin(x::String)

Create a ISIN Identifier.

# Example
```jldoctest; setup = :(using FinancialSymbology)
julia> isin_id = Isin("US0378331005")
"US0378331005"
```
"""
struct Isin <: Identifier
    s::AbstractString
end

"""
    Figi(x::String)

Create a FIGI Identifier.

# Example
```jldoctest; setup = :(using FinancialSymbology)
julia> figi_id = Figi("BBG001S5N8V8")
"BBG001S5N8V8"
```
"""
struct Figi <: Identifier
    s::AbstractString
end

"""
    Ticker(x::String)

Create a Ticker Identifier.

# Example
```jldoctest; setup = :(using FinancialSymbology)
julia> ticker_id = Ticker("AAPL US Equity")
"AAPL US Equity"
```
"""
struct Ticker <: Identifier
    s::AbstractString
    Ticker(s::AbstractString) = new(join(split(s), ' '))
end

"""
    Index(x::String)

Create a Index Identifier.

# Example
```jldoctest; setup = :(using FinancialSymbology)
julia> index_id = Index("990100")
"990100"
```
"""
struct Index <: Identifier
    s::AbstractString
end

Sedol(s::Identifier)= Sedol(s.s)
Sedol(s::Sedol) = s
Sedol(s::Missing) = missing

Cusip(s::Identifier) = Cusip(s.s)
Cusip(s::Cusip) = s
Cusip(s::Missing) = missing

Isin(s::Identifier)= Isin(s.s)
Isin(s::Isin) = s
Isin(s::Missing) = missing

Figi(s::Identifier) = Figi(s.s)
Figi(s::Figi) = s
Figi(s::Missing) = missing

Ticker(s::Identifier) = Ticker(s.s)
Ticker(s::Ticker) = s
Ticker(s::Missing) = missing

Index(s::Identifier) = Index(s.s)
Index(s::Ticker) = s
Index(s::Missing) = missing

end # module

max_stars_repo_name: robfs/figi

max_issues_repo_name: robfs/figi

max_forks_repo_name: robfs/figi


text:
# Autogenerated wrapper script for DuckDB_jll for aarch64-apple-darwin
export libduckdb

JLLWrappers.@generate_wrapper_header("DuckDB")
JLLWrappers.@declare_library_product(libduckdb, "@rpath/libduckdb.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libduckdb,
        "lib/libduckdb.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@testset "controller" begin
    @testset "options" begin
        @system SControllerOptions(Controller) begin
            a ~ preserve(extern)
            b ~ ::int(override)
        end
        o = (a=1, b=2)
        s = instance(SControllerOptions; options=o)
        @test s.a' == 1
        @test s.b == 2
    end

    @testset "config placeholder" begin
        @system SControllerConfigPlaceholder(Controller) begin
            a => 1 ~ preserve(parameter)
        end
        o = 0 => :a => 2
        s = instance(SControllerConfigPlaceholder; config=o)
        @test s.a' == 2
    end

    @testset "seed" begin
        @system SControllerSeed(Controller) begin
            a => 0 ± 1 ~ preserve
        end
        s1 = instance(SControllerSeed; seed=0)
        @test s1.a' == 0.6791074260357777
        s2 = instance(SControllerSeed; seed=0)
        @test s1.a' == s2.a'
    end

    @testset "no seed" begin
        @system SControllerNoSeed(Controller) begin
            a => 0 ± 1 ~ preserve
        end
        s1 = instance(SControllerNoSeed; seed=nothing)
        s2 = instance(SControllerNoSeed; seed=nothing)
        @test s1.a' != s2.a'
    end
end


text:
module CLagrangianFESpacesTests

using Test
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.TensorValues
using Gridap.Arrays
using Gridap.FESpaces
using Gridap.CellData

domain = (0,1,0,1)
partition = (2,2)
grid = CartesianGrid(domain,partition)

x = collect1d(get_node_coordinates(grid))

V = CLagrangianFESpace(Float64,grid)
matvecdata = ([],[],[])
matdata = ([],[],[])
vecdata = ([],[])
test_single_field_fe_space(V,matvecdata,matdata,vecdata)

@test get_cell_dofs(V) === get_cell_nodes(grid)
@test V.dof_to_node == collect(1:num_nodes(grid))
@test V.dof_to_comp == ones(Int,num_nodes(grid))
@test V.node_and_comp_to_dof == V.dof_to_node

u(x) = x[1]+x[2]
uh = interpolate(u,V)
uhx = get_free_values(uh)
ux = u.(x)
@test uhx ≈ ux

V = CLagrangianFESpace(VectorValue{2,Float64},grid)
test_single_field_fe_space(V,matvecdata,matdata,vecdata)

@test get_cell_dofs(V) == [[1,3,7,9,2,4,8,10], [3,5,9,11,4,6,10,12], [7,9,13,15,8,10,14,16], [9,11,15,17,10,12,16,18]]
@test V.dof_to_node == [1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9]
@test V.dof_to_comp == [1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2]
@test V.node_and_comp_to_dof == VectorValue{2,Int}[(1,2),(3,4),(5,6),(7,8),(9,10),(11,12),(13,14),(15,16),(17,18)]

u(x) = VectorValue(x[1]+x[2],x[2])
uh = interpolate(u,V)
uhx = get_free_values(uh)
ux = u.(x)
@test uhx ≈ collect1d(reinterpret(ux))

end # module


text:
export another_a, another_b,another_c

function another_a()
	1
end

another_b= 5

another_c(x)=x

type another_d <: AbstractString
	s
end

max_stars_repo_name: oxinabox/julia-vim-omnicomplete

max_issues_repo_name: oxinabox/julia-vim-omnicomplete

max_forks_repo_name: oxinabox/julia-vim-omnicomplete


text:
export AbstractVariable, degree
export is_clamped, is_marginalisation, is_moment_matching
export FoldLeftProdStrategy, FoldRightProdStrategy, CustomProdStrategy
export getmarginal, getmarginals, setmarginal!, setmarginals!, name, as_variable
export setmessage!, setmessages!

using Rocket

abstract type AbstractVariable end

## Variable collection type

abstract type AbstractVariableCollectionType end

struct VariableIndividual <: AbstractVariableCollectionType end

struct VariableVector <: AbstractVariableCollectionType
    index :: Int
end

struct VariableArray <: AbstractVariableCollectionType
    index :: CartesianIndex
end

indexed_name(::VariableIndividual, name::Symbol) = string(name)
indexed_name(seq::VariableVector, name::Symbol)  = string(name, "_", seq.index)
indexed_name(array::VariableArray, name::Symbol) = string(name, "_", join(array.index.I, "_"))

indexed_name(randomvar::AbstractVariable) = indexed_name(collection_type(randomvar), name(randomvar))

## Messages to Marginal product strategies

struct FoldLeftProdStrategy end
struct FoldRightProdStrategy end

struct CustomProdStrategy{F}
    prod_callback_generator :: F
end

"""
    messages_prod_fn(strategy, prod_constraint, form_constraint, form_check_strategy)

Returns a suitable prod computation function for a given strategy and constraints

See also: [`FoldLeftProdStrategy`](@ref), [`FoldRightProdStrategy`](@ref), [`CustomProdStrategy`](@ref)
"""
function messages_prod_fn end

messages_prod_fn(::FoldLeftProdStrategy, prod_constraint, form_constraint, form_check_strategy)       = prod_foldl_reduce(prod_constraint, form_constraint, form_check_strategy)
messages_prod_fn(::FoldRightProdStrategy, prod_constraint, form_constraint, form_check_strategy)      = prod_foldr_reduce(prod_constraint, form_constraint, form_check_strategy)
messages_prod_fn(strategy::CustomProdStrategy, prod_constraint, form_constraint, form_check_strategy) = strategy.prod_callback_generator(prod_constraint, form_constraint, form_check_strategy)

function marginal_prod_fn(strategy, prod_constraint, form_constraint, form_check_strategy)
    return let prod_fn = messages_prod_fn(strategy, prod_constraint, form_constraint, form_check_strategy)
        return (messages) -> as_marginal(prod_fn(messages))
    end
end

## Common functions

function degree end

add_pipeline_stage!(variable::AbstractVariable, stage) = error("Its not possible to add a new pipeline stage for $(variable)")

# Helper functions
# Getters

getmarginal(variable::AbstractVariable)                                      = getmarginal(variable, SkipInitial())
getmarginal(variable::AbstractVariable, skip_strategy::MarginalSkipStrategy) = apply_skip_filter(_getmarginal(variable), skip_strategy)

getmarginals(variables::AbstractArray{ <: AbstractVariable })                                      = getmarginals(variables, SkipInitial())
getmarginals(variables::AbstractArray{ <: AbstractVariable }, skip_strategy::MarginalSkipStrategy) = collectLatest(map(v -> getmarginal(v, skip_strategy), variables))

## Setters

### Marginals

setmarginal!(variable::AbstractVariable, marginal) = setmarginal!(getmarginal(variable, IncludeAll()), marginal)

setmarginals!(variables::AbstractArray{ <: AbstractVariable }, marginal::Distribution)    = _setmarginals!(Base.HasLength(), variables, Iterators.repeated(marginal, length(variables)))
setmarginals!(variables::AbstractArray{ <: AbstractVariable }, marginals)                 = _setmarginals!(Base.IteratorSize(marginals), variables, marginals)

function _setmarginals!(::Base.IteratorSize, variables::AbstractArray{ <: AbstractVariable }, marginals)
    foreach(zip(variables, marginals)) do (variable, marginal)
        setmarginal!(variable, marginal)
    end
end

function _setmarginals!(::Any, variables::AbstractArray{ <: AbstractVariable }, marginals)
    error("setmarginals!() failed. Default value is neither an iterable object nor a distribution.")
end

### Messages

setmessage!(variable::AbstractVariable, index::Int, message) = setmessage!(messageout(variable, index), message)
setmessage!(variable::AbstractVariable, message)             = foreach(i -> setmessage!(variable, i, message), 1:degree(variable))

setmessages!(variables::AbstractArray{ <: AbstractVariable }, message::Distribution)    = _setmessages!(Base.HasLength(), variables, Iterators.repeated(message, length(variables)))
setmessages!(variables::AbstractArray{ <: AbstractVariable }, messages)                 = _setmessages!(Base.IteratorSize(messages), variables, messages)

function _setmessages!(::Base.IteratorSize, variables::AbstractArray{ <: AbstractVariable }, messages)
    foreach(zip(variables, messages)) do (variable, message)
        setmessage!(variable, message)
    end
end

function _setmessages!(::Any, variables::AbstractArray{ <: AbstractVariable }, marginals)
    error("setmessages!() failed. Default value is neither an iterable object nor a distribution.")
end

##

function name(variable::AbstractVariable)
    return variable.name
end

## Helper functions

as_variable(x)                   = constvar(gensym(:as_var), x)
as_variable(v::AbstractVariable) = v
as_variable(t::Tuple)            = map(as_variable, t)


text:
function SecularAndHarmonic(arg0::jint, arg1::Vector{jdouble})
    return SecularAndHarmonic((jint, Vector{jdouble}), arg0, arg1)
end

function add_point(obj::SecularAndHarmonic, arg0::AbsoluteDate, arg1::jdouble)
    return jcall(obj, "addPoint", void, (AbsoluteDate, jdouble), arg0, arg1)
end

function approximate_as_polynomial_only(obj::SecularAndHarmonic, arg0::jint, arg1::AbsoluteDate, arg2::jint, arg3::jint, arg4::AbsoluteDate, arg5::AbsoluteDate, arg6::jdouble)
    return jcall(obj, "approximateAsPolynomialOnly", Vector{jdouble}, (jint, AbsoluteDate, jint, jint, AbsoluteDate, AbsoluteDate, jdouble), arg0, arg1, arg2, arg3, arg4, arg5, arg6)
end

function fit(obj::SecularAndHarmonic)
    return jcall(obj, "fit", void, ())
end

function get_fitted_parameters(obj::SecularAndHarmonic)
    return jcall(obj, "getFittedParameters", Vector{jdouble}, ())
end

function get_harmonic_amplitude(obj::SecularAndHarmonic)
    return jcall(obj, "getHarmonicAmplitude", jdouble, ())
end

function get_reference_date(obj::SecularAndHarmonic)
    return jcall(obj, "getReferenceDate", AbsoluteDate, ())
end

function mean_derivative(obj::SecularAndHarmonic, arg0::AbsoluteDate, arg1::jint, arg2::jint)
    return jcall(obj, "meanDerivative", jdouble, (AbsoluteDate, jint, jint), arg0, arg1, arg2)
end

function mean_second_derivative(obj::SecularAndHarmonic, arg0::AbsoluteDate, arg1::jint, arg2::jint)
    return jcall(obj, "meanSecondDerivative", jdouble, (AbsoluteDate, jint, jint), arg0, arg1, arg2)
end

function mean_value(obj::SecularAndHarmonic, arg0::AbsoluteDate, arg1::jint, arg2::jint)
    return jcall(obj, "meanValue", jdouble, (AbsoluteDate, jint, jint), arg0, arg1, arg2)
end

function osculating_derivative(obj::SecularAndHarmonic, arg0::AbsoluteDate)
    return jcall(obj, "osculatingDerivative", jdouble, (AbsoluteDate,), arg0)
end

function osculating_second_derivative(obj::SecularAndHarmonic, arg0::AbsoluteDate)
    return jcall(obj, "osculatingSecondDerivative", jdouble, (AbsoluteDate,), arg0)
end

function osculating_value(obj::SecularAndHarmonic, arg0::AbsoluteDate)
    return jcall(obj, "osculatingValue", jdouble, (AbsoluteDate,), arg0)
end

function reset_fitting(obj::SecularAndHarmonic, arg0::AbsoluteDate, arg1::Vector{jdouble})
    return jcall(obj, "resetFitting", void, (AbsoluteDate, Vector{jdouble}), arg0, arg1)
end


text:
module RelPerm

#TODO add capillary pressure here -- see equations 8.11 and 8.12 in "introduction to reservoir simulation using matlab/gnu octave" for the brooks-corey and van genuchten forms of capillary pressure

module Corey
	function S_wn(S_w, S_wir, S_nir)
		return (S_w - S_wir) / (1 - S_wir - S_nir)
	end
	function kr_w(S_w, S_wir, S_nir, N_w, k0_w)
		if S_w < S_wir
			return 0.0
		elseif S_w > 1 - S_nir
			return k0_w
		else
			return k0_w * (S_wn(S_w, S_wir, S_nir)) ^ N_w
		end
		#return k0_w * (S_wn(S_w, S_wir, S_nir)) ^ N_w
	end
	function kr_n(S_n, S_wir, S_nir, N_n, k0_n)
		if S_n < S_nir
			return 0.0
		elseif S_n > 1 - S_wir
			return k0_n
		else
			return k0_n * (1 - S_wn(1 - S_n, S_wir, S_nir)) ^ N_n
		end
		#return k0_n * (1 - S_wn(1 - S_n, S_wir, S_nir)) ^ N_n
	end
end

end


text:
include("handle_event_tests/test_unhandled_event.jl")
include("handle_event_tests/test_parent_child_event_handling.jl")


text:
using Plots;
using Gadfly

W_fejer(omega,N) = 1/N * (sin(omega*N/2)/sin(omega/2))^2

N= 30
omega=range(-pi+1e-3,pi-1e-3,length=10000)
W0 = W_fejer.(1e-13,N)
plot(omega,W_fejer.(omega,N)./W0,yaxis=:log,label="",
xlabel="omega",ylabel="W_B(omega)/W(0)")
savefig("/home/thomas/atom-amd64/usr/share/atom/example/FejerKernel")


text:
module TextParse

using CodecZlib, WeakRefStrings, Dates, Nullables, DoubleFloats

include("VectorBackedStrings.jl")
include("lib/compat.jl")
include("util.jl")
include("field.jl")
include("record.jl")

include("utf8optimizations.jl")

include("guesstype.jl")
include("csv.jl")

end # module


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

struct __enum_EventType <: ProtoEnum
    JudgeResult::Int32
    MatchStart::Int32
    MatchStop::Int32
    FirstHalfStart::Int32
    SecondHalfStart::Int32
    OvertimeStart::Int32
    PenaltyShootoutStart::Int32
    __enum_EventType() = new(0,1,2,3,4,5,6)
end #struct __enum_EventType
const EventType = __enum_EventType()

struct __enum_JudgeResultEvent_ResultType <: ProtoEnum
    PlaceKick::Int32
    GoalKick::Int32
    PenaltyKick::Int32
    FreeKickRightTop::Int32
    FreeKickRightBot::Int32
    FreeKickLeftTop::Int32
    FreeKickLeftBot::Int32
    __enum_JudgeResultEvent_ResultType() = new(0,1,2,3,4,5,6)
end #struct __enum_JudgeResultEvent_ResultType
const JudgeResultEvent_ResultType = __enum_JudgeResultEvent_ResultType()

mutable struct JudgeResultEvent <: ProtoType
    _type::Int32
    offensive_team::Int32
    reason::AbstractString
    JudgeResultEvent(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct JudgeResultEvent

mutable struct EventArguments <: ProtoType
    judge_result::JudgeResultEvent
    field::Field
    EventArguments(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct EventArguments
const __oneofs_EventArguments = Int[1,1]
const __oneof_names_EventArguments = [Symbol("argument")]
meta(t::Type{EventArguments}) = meta(t, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, __oneofs_EventArguments, __oneof_names_EventArguments, ProtoBuf.DEF_FIELD_TYPES)

export EventType, EventArguments, JudgeResultEvent_ResultType, JudgeResultEvent


text:
include("../examples/rijke.jl")
#include("../src/clvs.jl")
using PyPlot
using OrdinaryDiffEq
using LinearAlgebra
using JLD
using SharedArrays
using Distributed
function compute_rayleigh_criterion(n_samples)
    d_u = 3
    Jray = SharedArray{Float64}(n_samples)
    Jray .= 0.
    beta = LinRange(6.0,9.0,n_samples)
    nSteps = 20000
    nRunup = 400000
    answer = @distributed for i = 1:n_samples
    	println("beta = ", beta[i])
    	s = [beta[i], 0.2]
    	u = Rijke_ODE(rand(N), s, nRunup)
    	for j = 1:nSteps
        	u = Rijke_ODE(u, s, 1)
			Jray[i] += 0.5*sum(u[Ng+1:2*Ng].*zetaj.*u[Ng+1:2*Ng])/nSteps
    	end
    	@show Jray[i]
    end
    wait(answer)
    save("../data/attractor/Jray.jld", "beta", beta,
     "Jray", Jray)
end
function plot_rayleigh()
	filename = string("../data/attractor/",
					  "Jray_actual.jld")
	X = load(filename)
	Jray = X["Jray"]
	beta = X["beta"]
	fig, ax = subplots(1,1)
	ax.plot(beta, Jray, "b.", ms=4.0)
	ax_in = fig.add_axes([0.2,0.6,0.2,0.2])
	ax_in.plot(beta[6.7 .<= beta .<= 7.3],
			Jray[6.7 .<= beta .<= 7.3], "b.",ms=4.0)
	ax.xaxis.set_tick_params(labelsize=25)
	ax.yaxis.set_tick_params(labelsize=25)
	ax_in.xaxis.set_tick_params(labelsize=20)
	ax_in.yaxis.set_tick_params(labelsize=20)

	ax.set_xlabel(L"$\beta$", fontsize=25)
	ax.set_ylabel(L"$\langle J_{\rm ray}\rangle$", fontsize=25)
	ax.grid(true)
	ax_in.grid(true)
end

max_stars_repo_name: nishaChandramoorthy/energies

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies


text:
"""
    Monin_Obukhov_length(Tair, pressure, ustar, H; constants)
    Monin_Obukhov_length!(df;constants=BigleafConstants())
    Monin_Obukhov_length(df;constants=BigleafConstants())

calculates the Monin-Obukhov length.

# Arguments
- `Tair`      : Air temperature (degC)
- `pressure`  : Atmospheric pressure (kPa)
- `ustar`     : Friction velocity (m s-1)
- `H`         : Sensible heat flux (W m-2)
- `df`        : DataFrame containing the above variables
optional
- `constants=`[`BigleafConstants`](@ref)`()`: Dictionary with entries
  - `Kelvin` - conversion degree Celsius to Kelvin
  - `cp` - specific heat of air for constant pressure (J K-1 1)
  - `k` - von Karman constant (-)
  - `g` - gravitational acceleration (m s-2)

# Details
The Monin-Obukhov length (L) is given by:

``L = - (\\rho * cp * ustar^3 * Tair) / (k * g * H)``

where ``\\rho`` is air density (kg m-3).

# Note
Note that L gets very small for very low ustar values with implications
for subsequent functions using L as input. It is recommended to filter
data and exclude low ustar values (ustar < ~0.2) beforehand.

# Value
Monin-Obukhov length L (m). The non-mutating DataFrame variant returns a vector,
the mutating variant add or modifies column `:MOL`.

# References
Foken, T, 2008: Micrometeorology. Springer, Berlin, Germany.

# See also
[`stability_parameter`](@ref)

```@example; output = false
Monin_Obukhov_length(
  Tair=25,pressure=100,
  ustar=seq(0.2,1,0.1),H=seq(40,200,20))
```
"""
function Monin_Obukhov_length(Tair::FT, pressure, ustar, H; constants=BigleafConstants()
  ) where FT
  rho  = air_density(Tair, pressure; constants)
  TairK = Tair + oftype(Tair, constants.Kelvin)
  #MOL  = (-rho*constants.cp*ustar^3*TairK) / (constants.k*constants.g*H)
  MOL  = (-rho*FT(constants.cp)*ustar^3*TairK) / (FT(constants.k)*FT(constants.g)*H)
end
function Monin_Obukhov_length!(df;constants=BigleafConstants())
  fr = (args...) -> Monin_Obukhov_length(args...; constants)
  transform!(df, SA[:Tair, :pressure, :ustar, :H] => ByRow(fr) => :MOL)
end
# function Monin_Obukhov_length(df::DFTable; kwargs...)
#     tmp = Monin_Obukhov_length.(df.Tair, df.pressure, df.ustar, df.H; kwargs...)
# end

"""
    stability_parameter(z,d,MOL)
    stability_parameter(z,d,Tair, pressure, ustar, H; constants)
    stability_parameter!(df::AbstractDataFrame; z,d, MOL=nothing, constants)

calculates stability parameter "zeta", a parameter characterizing stratification in the
lower atmosphere.

# Arguments
- `z`   : height (m)
- `d`   : Zero-plane displacement height (m)
- `MOL` : Monin-Obukhov-length L (m)
- `df`  : DataFrame containting the variables required by [`Monin_Obukhov_length`](@ref)
optional
- `constants=`[`BigleafConstants`](@ref)`()`

In the second variant and if `MOL=nothing` in the DataFrame variants, MOL is computed by
[`Monin_Obukhov_length`](@ref).

# Details
The stability parameter ``\\zeta`` is given by:

``\\zeta = (z - d) / L``

where L is the Monin-Obukhov length (m), calculated by
[`Monin_Obukhov_length`](@ref). The displacement height can
be estimated from the function [`roughness_parameters`](@ref).

# Value
``\\zeta``: stability parameter (-). The nonmutainting DataFrame variant returns a vector.
The mutating variant modifies or adds column [`:zeta`].

```jldoctest; output = false
using DataFrames
df = DataFrame(Tair=25.0, pressure=100.0, ustar=0.2:0.1:1.0, H=40:20.0:200)
z=40;d=15
zeta = stability_parameter.(z,d, df.Tair, df.pressure, df.ustar, df.H)
all(zeta .< 0)
# output
true
```
"""
function stability_parameter(z,d,MOL)
  zeta = (z - d) / MOL
end
function stability_parameter(z,d,Tair, pressure, ustar::FT, H;
  constants=BigleafConstants()) where FT
  MOL = Monin_Obukhov_length(Tair, pressure, ustar, H; constants);
  stability_parameter(FT(z),FT(d),MOL)
end
function stability_parameter!(df::AbstractDataFrame; z,d, MOL=nothing,
  constants=BigleafConstants())
  if isnothing(MOL)
    ft = (args...) -> stability_parameter.(z,d,args...; constants)
    transform!(df, SA[:Tair, :pressure, :ustar, :H] => ft => :zeta)
  else
    ft = () -> stability_parameter.(z,d,MOL)
    transform!(df, [] => ft => :zeta)
  end
end
# function stability_parameter(df::DFTable; z,d, MOL=nothing,
#   constants=BigleafConstants())
#   if isnothing(MOL) MOL = Monin_Obukhov_length(df; constants); end
#   stability_parameter.(z,d,MOL)
# end

abstract type StabilityCorrectionMethod end
struct Dyer1970 <: StabilityCorrectionMethod end
struct Businger1971 <: StabilityCorrectionMethod end
struct NoStabilityCorrection <: StabilityCorrectionMethod end

"""
    stability_correction(zeta;
      stab_formulation=Dyer1970())
    stability_correction(z,d, Tair,pressure,ustar,H; constants,
      stab_formulation=Dyer1970())
    stability_correction!(df; zeta, z, d;
      stab_formulation=Dyer1970(), constants =BigleafConstants())

Integrated Stability Correction Functions for Heat and Momentum

# Arguments
- `zeta`             : Stability parameter zeta (-)
- `Tair`,`pressure`,`ustar`,`H` : see [`Monin_Obukhov_length`](@ref)
- `z`,`d`            : see [`stability_parameter`](@ref)
- `df`  : DataFrame containting the variables required by [`Monin_Obukhov_length`](@ref)
- `stab_formulation` : Formulation for the stability function. Either
            `Dyer1970()`, or `Businger1971()` or `NoStabilityCorrection()`

In the second and third form computes `zeta` by [`stability_parameter`](@ref) and
[`Monin_Obukhov_length`](@ref) and requires respective arguments.

# Details
These dimensionless values are needed to correct deviations
from the exponential wind profile under non-neutral conditions.
The functions give the integrated form of the universal functions. They
depend on the value of the stability parameter ``\\zeta``,
a function of heigh `z`,
which can be calculated from the function [`stability_parameter`](@ref).
The integration of the universal functions is:

``\\psi = -x * \\zeta``

for stable atmospheric conditions (``\\zeta`` >= 0), and

``\\psi = 2 * log( (1 + y(\\zeta)) / 2)``

for unstable atmospheric conditions (``\\zeta`` < 0).

The different formulations differ in their value of x and y.

# Value
a NamedTuple with the following columns:
- `psi_h`: the value of the stability function for heat and water vapor (-)
- `psi_m`: the value of the stability function for momentum (-)

# References
- Dyer, A_J., 1974: A review of flux-profile relationships.
  Boundary-Layer Meteorology 7, 363-372.
- Dyer, A. J., Hicks, B_B., 1970: Flux-Gradient relationships in the
  constant flux layer. Quart. J. R. Meteorol. Soc. 96, 715-721.
- Businger, J_A., Wyngaard, J. C., Izumi, I., Bradley, E. F., 1971:
  Flux-Profile relationships in the atmospheric surface layer.
  J. Atmospheric Sci. 28, 181-189.
- Paulson, C_A., 1970: The mathematical representation of wind speed
  and temperature profiles in the unstable atmospheric surface layer.
  Journal of Applied Meteorology 9, 857-861.
  Foken, T, 2008: Micrometeorology. Springer, Berlin, Germany.

# Examples
```jldoctest; output = false
using DataFrames
zeta = -2:0.5:0.5
df2 = DataFrame(stability_correction.(zeta; stab_formulation=Businger1971()))
propertynames(df2) == [:psi_h, :psi_m]
# output
true
```
"""
function stability_correction(zeta::FT; stab_formulation=Dyer1970()) where FT
  # integration of universal functions (after Paulson1970 and Foken 2008)
  stab_formulation isa NoStabilityCorrection && return(
    (psi_h = zero(zeta), psi_m = zero(zeta)))
  ismissing(zeta) && return((psi_h = missing, psi_m = missing))
  is_stable = zeta >= 0
  if is_stable
    x_h::FT, x_m::FT = get_stability_coefs_stable(stab_formulation)
    psi_h = x_h * zeta
    psi_m = x_m * zeta
  else
    y_h::FT, y_m::FT = get_stability_coefs_unstable(stab_formulation, zeta)
    psi_h = 2 * log( (1 + y_h ) / 2)
    psi_m = 2 * log( (1 + y_m ) / 2) +
                      log( ( 1 + y_m^2 ) / 2)
                      -2 * atan(y_m) + pi/2
  end
  (;psi_h, psi_m)
end

get_stability_coefs_stable(::Businger1971) = (x_h = -7.8, x_m = -6)
function get_stability_coefs_unstable(::Businger1971, zeta)
  y_h = 0.95 * ( 1 - 11.6 * zeta)^0.5
  y_m = (1 - 19.3*zeta)^0.25
  (;y_h, y_m)
end

get_stability_coefs_stable(::Dyer1970) = (x_h = -5, x_m = -5)
function get_stability_coefs_unstable(::Dyer1970, zeta)
  y_h       = (1 - 16 * zeta)^0.5
  y_m       = (1 - 16 * zeta)^0.25
  (;y_h, y_m)
end

function stability_correction(z,d, Tair::Union{Missing,Number},pressure,ustar::FT,H;
  stab_formulation=Dyer1970(), constants=BigleafConstants()) where FT
  stab_formulation isa NoStabilityCorrection && return(
    (psi_h = zero(FT), psi_m = zero(FT)))
  FT <: Missing && return(missing)
  MOL = Monin_Obukhov_length(Tair,pressure,ustar,H; constants)
  zeta  = stability_parameter(FT(z),FT(d),MOL)
  psis = stability_correction(zeta; stab_formulation)
  psis
end
# function stability_correction(Tair,pressure,ustar,H, z,d; kwargs...)
#   Tables.columns(
#     stability_correction.(Tair,pressure,ustar,H, z,d; kwargs...))
# end
# function stability_correction(df::DFTable, z,d; kwargs...)
#   tmp1 = stability_correction.(df.Tair, df.pressure, df.ustar, df.H, z, d; kwargs...)
#   tmp2 = Tables.columns(tmp1)
# end
function stability_correction!(df; zeta=nothing, z=nothing, d=nothing,
  stab_formulation=Dyer1970(), constants =BigleafConstants())
  # cannot dispatch on keyword argument, hence need if-clause
  if stab_formulation isa NoStabilityCorrection
    zeroT = :ustar in propertynames(df) ? zero(eltype(skipmissing(df.ustar))) : zero(z)
    df[!,:psi_h] .= zeroT
    df[!,:psi_m] .= zeroT
    return(df)
  end
  if !isnothing(zeta)
    ft = () -> stability_correction.(zeta; stab_formulation)
    transform!(df, [] => ft => AsTable)
  else
    ft = (args...) -> stability_correction.(z, d, args...; stab_formulation, constants)
    transform!(df, SA[:Tair,:pressure,:ustar,:H] => ft => AsTable)
  end
end

function stability_correction(df::DFTable; z, d,
  stab_formulation=Dyer1970(), constants =BigleafConstants()
  )
  # do not provide zeta, because can simply invoke stability_correction.(zeta)
  if stab_formulation isa NoStabilityCorrection
    zeroT = zero(eltype(skipmissing(df.ustar)))
    rows = map(Tables.rows(df)) do row
      (psi_h = zeroT, psi_m = zeroT)
    end
    return(Tables.columns(rows))
  end
  Tables.columns(stability_correction.(
    z, d, df.Tair, df.pressure, df.ustar, df.H; stab_formulation, constants))
end


text:
using FileIO,MAT,NeuroAnalysis,Images,Statistics,StatsBase,Plots,ProgressMeter,Interact,MsgPack,YAML

function imagepatch(img,ppd,sizedeg;topleftdeg=nothing)
    sizepx = round.(Int,sizedeg.*ppd)
    isize = size(img)
    isizepx = isize[1:2]
    any(sizepx.>isizepx) && error("Patch$sizepx is larger than image$isizepx.")

    maxtopleftpx = isizepx .- sizepx .+ 1
    if isnothing(topleftdeg)
        topleftpx = rand.(map(i->1:i,maxtopleftpx))
    end
    ppx = map((i,j)->i:i+j-1,topleftpx,sizepx)
    if length(isize)==2
        patch = img[ppx...]
    else
        patch = img[ppx...,map(i->1:i,isize[3:end])...]
    end
    return patch
end

function newimageset(imgdbroot,imgfile,ppd;n=100,sizedeg=(3,3),sizepx=(32,32))
    imgfiles=[]
    for (root,dirs,files) in walkdir(imgdbroot)
        append!(imgfiles,joinpath.(root,filter(f->occursin(imgfile,f),files)))
    end
    isempty(imgfiles) && error("No valid image file found in the image database.")

    imgpatchs=[]
    p = Progress(n,desc="Sampling Images ... ")
    for i in 1:n
        file = rand(imgfiles)
        ext = splitext(file)[2]
        if ext == ".mat"
            img = first(matread(file)).second
        end
        push!(imgpatchs,imresize_antialiasing(imagepatch(img,ppd,sizedeg),sizepx))
        next!(p)
    end
    return imgpatchs
end

function combineimageset(imgsets...;ps=fill(1.0/length(imgsets),length(imgsets)),n=100,c=nothing,normfun=x->clampscale(x,3))
    ns = round.(Int,ps.*n)
    imgs = []
    p = Progress(sum(ns),desc="Picking Images ... ")
    for j in eachindex(ns)
        is = sample(1:length(imgsets[j]),ns[j],replace=false)
        if isnothing(c)
            for i in is
                t = imgsets[j][i][:,:,rand(1:3)]
                any(isnan.(t)) && continue
                push!(imgs,normfun(t))
                next!(p)
            end
        end
    end
    finish!(p)
    return imgs
end

function checkimageset!(imgset;alpha=0.8,n=Inf,sdrange=[0.23,0.33],unique=true)
    di = []
    for i in eachindex(imgset)
        any(isnan.(imgset[i])) && push!(di,i)
    end
    deleteat!(imgset,di);di=[]
    for i in eachindex(imgset)
        any(isinf.(imgset[i])) && push!(di,i)
    end
    deleteat!(imgset,di);di=[]
    for i in eachindex(imgset)
        sdrange[1] < std(imgset[i]) < sdrange[2] || push!(di,i)
    end
    deleteat!(imgset,di)
    if unique
        l=length(imgset);di=falses(l)
        p = Progress(l-1,desc="Checking Images ... ")
        for i in 1:(l-1)
            di[i] && continue
            Threads.@threads for j in (i+1):l
                di[j] && continue
                cor(vec(imgset[i]),vec(imgset[j])) > alpha && (di[j]=true)
            end
            next!(p)
        end
        deleteat!(imgset,findall(di));finish!(p)
    end

    l = length(imgset)
    l > n && deleteat!(imgset,(n+1):l)
    imgset
end

function saveunityrawtexture(fn,imgset)
    ds = size(imgset[1]);imgs=[]
    if length(ds) < 3
        pbit = 8
        for i in imgset
            # Unity Texture2D RawData unfold from left->right, then bottom->up
            ps = vec(reverse(permutedims(i),dims=2))
            # Pack Gray in UInt8
            ps = [reinterpret(N0f8(p)) for p in ps]
            push!(imgs,ps)
        end
    else
        pbit = 32
        for i in imgset
            # Unity Texture2D RawData unfold from left->right, then bottom->up
            ps = vec(reverse(permutedims(i),dims=2))
            # Pack RGBA32 in UInt32
            ps = [parse(UInt32,"0x"*hex(RGBA{N0f8}(p,p,p,1),:rrggbbaa)) for p in ps]
            push!(imgs,ps)
        end
    end
    imgdata = Dict("ImageSize"=>[map(i->UInt16(i),ds)...],"Images"=>imgs)
    # MessagePack Binary File
    write("$fn.mpis$pbit", pack(imgdata))
    return imgdata
end

stimuliroot = "../NaturalStimuli"
## Imageset from McGillCCIDB
idb = "McGillCCIDB"
itype = "LMS"
ppd = 90 # should be close to UPennNIDB
n = 60000
sizedeg = (3,3)
sizepx = (64,64)
imgset_MG = newimageset(joinpath(stimuliroot,idb),Regex("\\w*$itype.mat"),ppd;n,sizedeg,sizepx)
save(joinpath(stimuliroot,"$(idb)_$(itype)_n$(n)_sizedeg$(sizedeg)_sizepx$(sizepx).jld2"),"imgset",imgset_MG)
imgset_MG = load(joinpath(stimuliroot,"$(idb)_$(itype)_n$(n)_sizedeg$(sizedeg)_sizepx$(sizepx).jld2"),"imgset")

## Imageset from UPennNIDB
idb = "UPennNIDB"
itype = "LMS"
ppd = 92
n = 60000
sizedeg = (3,3)
sizepx = (64,64)
imgset_UP = newimageset(joinpath(stimuliroot,idb),Regex("\\w*$itype.mat"),ppd;n,sizedeg,sizepx)
save(joinpath(stimuliroot,"$(idb)_$(itype)_n$(n)_sizedeg$(sizedeg)_sizepx$(sizepx).jld2"),"imgset",imgset_UP)
imgset_UP = load(joinpath(stimuliroot,"$(idb)_$(itype)_n$(n)_sizedeg$(sizedeg)_sizepx$(sizepx).jld2"),"imgset")

timgset = imgset_UP
@manipulate for ii in 1:length(timgset), ci in 1:size(timgset[1],3)
    t = timgset[ii][:,:,ci]
    f = scaleminmax(extrema(t)...)
    Gray.(f.(t))
end

## Combine Imageset
normfun = (x;c=0.5,r=0.5,low=0.1,high=0.9)->begin
    xx = clamp.(x,quantile(vec(x),low),quantile(vec(x),high))
    xx .-= mean(xx)
    r*xx/maximum(abs.(xx)) .+ c
end
c=0.5 # mean luminance
r=0.5 # max extrema relative to c
imgset = combineimageset(imgset_MG,imgset_UP;ps=[0.7,0.3],n=40000,normfun=x->normfun(x;c,r))

timgset=imgset
@manipulate for ii in 1:length(timgset), ci in 1:size(timgset[1],3)
    t = timgset[ii][:,:,ci]
    ta = adjust_histogram(t, Equalization(nbins = 1024))
    ps,f1,f2 = powerspectrum2(t,sizepx[1]/sizedeg[1],freqrange=[-6,6])
    psa,f1,f2 = powerspectrum2(ta,sizepx[1]/sizedeg[1],freqrange=[-6,6])
    ps = log10.(ps);psa=log10.(psa);plims=(minimum([ps;psa]),maximum([ps;psa]))

    p = plot(layout=(3,2),leg=false,size=(2*300,3*300))
    heatmap!(p[1,1],t,color=:grays,aspect_ratio=1,frame=:none,yflip=true,clims=(0,1))
    histogram!(p[2,1],vec(t),xlims=[0,1],xticks=0:0.1:1);vline!(p[2,1],[mean(t),std(t)],lw=3,color=[:hotpink,:seagreen])

    heatmap!(p[1,2],ta,color=:grays,aspect_ratio=1,frame=:none,yflip=true,clims=(0,1))
    histogram!(p[2,2],vec(ta),xlims=[0,1],xticks=0:0.1:1);vline!(p[2,2],[mean(ta),std(ta)],lw=3,color=[:hotpink,:seagreen])

    heatmap!(p[3,1],f2,f1,ps,aspect_ratio=:equal,frame=:semi,color=:plasma,clims=plims)
    heatmap!(p[3,2],f2,f1,psa,aspect_ratio=:equal,frame=:semi,color=:plasma,clims=plims)
end

## Check and Save Imageset
sdrange=[0.23,0.33]
checkimageset!(imgset;sdrange)

imgsetname = "n$(length(imgset))_sizedeg$(sizedeg)_c$(c)_r$(r)_sd$(sdrange)"
save(joinpath(stimuliroot,"$imgsetname.jld2"),"imgset",imgset)
saveunityrawtexture(joinpath(stimuliroot,imgsetname),imgset)
matwrite(joinpath(stimuliroot,"$imgsetname.mat"),Dict("imgset"=>imgset))
imgset = load(joinpath(stimuliroot,"$imgsetname.jld2"),"imgset")

## Gaussian Noise Imageset
gnormfun = x->begin
    xx = clamp.(x,-2.5,2.5)
    xx .-= mean(xx)
    0.5*xx/maximum(abs.(xx)) .+ 0.5
end
gimgset = [gnormfun(randn(sizepx)) for _ in 1:10000]

## Hartley Subspace Imageset
hs = hartleysubspace(kbegin=0.2,kend=5.0,dk=0.2,addhalfcycle=true,shape=:circle)
himgset = map(i -> begin
          ss = cas2sin(i...)
          grating(θ = ss.θ, sf = ss.f, phase = ss.phase, size = sizedeg, ppd = 30)
          end, hs)

## Imageset Mean Spectrum
pss,f1,f2 = powerspectrums2(imgset,sizepx[1]/sizedeg[1],freqrange=[-6,6])
gpss,gf1,gf2 = powerspectrums2(gimgset,sizepx[1]/sizedeg[1],freqrange=[-6,6])
hpss,hf1,hf2 = powerspectrums2(himgset,30,freqrange=[-6,6])

mps = log10.(reduce((i,j)->i.+j,pss))/length(pss)
gmps = log10.(reduce((i,j)->i.+j,gpss))/length(gpss)
hmps = log10.(reduce((i,j)->i.+j,hpss))/length(hpss)

mps[f1.==0,f2.==0].=minimum(mps)
gmps[gf1.==0,gf2.==0].=minimum(gmps)
hmps[hf1.==0,hf2.==0].=minimum(hmps)

plotmeanpowerspectrum = () -> begin
    p = plot(layout=(3,1),leg=false,size=(300,3*300))
    heatmap!(p[1],f2,f1,mps,aspect_ratio=:equal,frame=:none,color=:plasma,title="Natural Image")
    heatmap!(p[2],gf2,gf1,gmps,aspect_ratio=:equal,frame=:none,color=:plasma,title="Gaussian Noise")
    heatmap!(p[3],hf2,hf1,hmps,aspect_ratio=:equal,frame=:none,color=:plasma,title="Hartley Subspace")
    p
end

plotmeanpowerspectrum()
foreach(ext->savefig(joinpath(@__DIR__,"ImageSets_Mean_PowerSpectrum.$ext")),[".png",".svg"])

max_stars_repo_name: VLAB-ExperimentSystem/VLab

max_issues_repo_name: VLAB-ExperimentSystem/VLab

max_forks_repo_name: VLAB-ExperimentSystem/VLab


text:
using YAML
using ArgParse

function add_args(s::ArgParseSettings)
    @add_arg_table(
        s,
        "--cfg", arg_type = String, default = "cfg/default.yaml",
        "--vinput", arg_type = Int64, default = 16,
        "--noise", arg_type = Float64, default = 5.0,
        "--theta_plus", arg_type = Float64, default = 0.4,
        "--refrac_e", arg_type = Float64, default = 15.0,
        "--refrac_i", arg_type = Float64, default = 4.0,
        "--ttheta", arg_type = Float64, default = 4.7,
        "--tdop", arg_type = Float64, default = 6.2,
        "--lr", arg_type = Float64, default = 6.0,
        "--target", arg_type = Float64, default = 0.4,
        "--wrandom", arg_type = Bool, default = true,
        "--wstart", arg_type = Float64, default = 0.9,
        "--winh", arg_type = Float64, default = 0.4,
        "--wmax", arg_type = Float64, default = 3.5,
        "--mu", arg_type = Float64, default = 0.1
    )
    s
end

function get_config(argtable::ArgParseSettings)
    argtable = add_args(argtable)
    args = parse_args(argtable)
    cfg = YAML.load_file(args["cfg"])

    for k in ["noise", "vinput", "theta_plus", "refrac_e", "refrac_i", "lr",
              "target", "wrandom", "wstart", "winh", "wmax", "mu"]
        cfg[k] = args[k]
    end

    for k in ["ttheta", "tdop"]
        cfg[k] = args[k] * 100
    end
    args, cfg
end


text:
import .Metaheuristics
const MH = Metaheuristics

Base.@kwdef struct MetaH{T<:MH.Algorithm} <: ComradeSoss.AbstractOptimizer
    alg::T=MH.ECA(options=MH.Options(f_calls_limit=10^5))
end

function optimize(opt::MetaH, lj::Soss.ConditionalModel)
    tc = ascube(lj)
    lower = transform(tc, zeros(dimension(tc)))
    upper = transform(tc, 0.999*ones(dimension(tc)))
    upflat, _ = ParameterHandling.flatten(upper)
    lowflat, unflatten = ParameterHandling.flatten(lower)

    bounds = [lowflat upflat]'
    f(x) = -logdensity(lj, unflatten(x))
    res = MH.optimize(f, bounds, opt.alg)
    show(res)

    return unflatten(MH.minimizer(res)), -MH.minimum(res)
end


text:
# License for this file: MIT (expat)
# Copyright 2021, DLR Institute of System Dynamics and Control
#
# This file is part of module ModiaPlot

###################################################################################
#
# Generic definitions that might be specialized from special result data structures
#
###################################################################################

"""
    hasSignal(result, name)

Returns `true` if signal `name` is available in `result`.
"""
hasSignal(result, name) = haskey(result, name)

"""
    (isConstant, signal) = getRawSignal(result, name)

Returns result time series `signal` of `name` (an error is raised, if `name` is not known).

If `isConstant=false`, then `signal[i]` is the `value` of the signal at time instant `i`.

If `isConstant=true`, then `signal` is the `value` of the signal and this value
holds for all time instants.

`typeof(value)` must be either `<:Number` or `<:AbstractArray` with `eltype(value) <: Number`.
"""
getRawSignal(result, name) = (false, result[name])

"""
    getNames(result)

Return a vector of the names that are present in result.
"""
getNames(result) = collect( keys(result) )

"""
    getDefaultHeading(result)

Return default heading.
"""
getDefaultHeading(result) = ""

"""
    (isConstant, sigSize, sigElType, sigUnit) = getSignalInfo(result, name)

Return information about a signal, given the `name` of the signal:

- isConstant: = true, if signal is constant.
- sigSize: size(signal)
- sigElType: Element type of signal without unit
- sigUnit: Unit of signal
"""
function getSignalInfo(result, name)
    (isConstant, signal) = getRawSignal(result,name)
    if isnothing(signal)
        return (isConstant, -1, nothing, nothing, nothing)
    elseif ismissing(signal)
        return (isConstant, -1, missing, missing, missing)
    elseif length(signal) == 0
        return (isConstant, -1, nothing, nothing, nothing)
    end
    value       = isConstant ? signal : signal[1]
    valueSize   = size(value)
    valueUnit   = unit(value[1])

    if typeof(value) <: MonteCarloMeasurements.Particles
        elTypeAsString = string(typeof(ustrip.(value[1])))
        nparticles     = length(value)
        valueElType    = "MonteCarloMeasurements.Particles{" * elTypeAsString * ",$nparticles}"
    elseif typeof(value) <: MonteCarloMeasurements.StaticParticles
        elTypeAsString = string(typeof(ustrip.(value[1])))
        nparticles     = length(value)
        valueElType    = "MonteCarloMeasurements.StaticParticles{" * elTypeAsString * ",$nparticles}"
    else
        valueElType = typeof( ustrip.(value) )
    end
    nTime = length(signal)
    return (isConstant, nTime, valueSize, valueElType, valueUnit)
end

###################################################################################
#
# Utility functions to return the desired part of a signal and the heading
#
###################################################################################

"""
    (signal, isConstant, arrayName, arrayIndices) = getSignal(result, name)

Return the signal defined by `name::AbstractString`.
`name` may include an array range, such as "a.b.c[2:3,5]".
In this case `arrayName` is the name without the array indices,
such as `"a.b.c"` and arrayIndices is a tuple with the array indices,
such as `(2:3, 5)`. Otherwise `arrayName = name, arrayIndices=()`.

In case the signal is not known or `name` cannot be interpreted,
`(nothing, false, name, ())` is returned.

If `isConstant=false`, then `signal[i]` is the `value` of the signal at time instant `i`.

If `isConstant=true`, then `signal` is the `value` of the signal and this value
holds for all time instants.

`typeof(value)` is either `<:Number` or `<:AbstractArray` with `eltype(value) <: Number`.

The following `Number` types are currently supported for signal elements sigElement:

1. `convert(Float64, sigElement)` is supported
   (for example Float32, Float64, DoubleFloat, Rational, Int32, Int64, Bool).

2. Measurements.Measurement{<Type of (1)>}.

3. MonteCarloMeasurements.StaticParticles{<Type of (1)>}.

4. MonteCarloMeasurements.Particles{<Type of (1)>}.
"""
function getSignal(result, name::AbstractString)
    sigPresent = false
    if hasSignal(result, name)
        (isConstant, sig) = getRawSignal(result, name)
        if !( isnothing(sig) || ismissing(sig) || length(sig) == 0 )
            sigPresent   = true
            value        = isConstant ? sig : sig[1]
            arrayName    = name
            arrayIndices = ndims(value) == 0 ? () : Tuple(1:Int(ni) for ni in size(value))
        end

    else
        # Handle signal arrays, such as a.b.c[3] or a.b.c[2:3, 1:5, 3]
        if name[end] == ']'
            i = findlast('[', name)
            if i >= 2
                arrayName = name[1:i-1]
                indices   = name[i+1:end-1]
                if hasSignal(result, arrayName)
                    (isConstant, sig2) = getRawSignal(result, arrayName)
                    if !( isnothing(sig2) || ismissing(sig2) || length(sig2) == 0 )
                        # Determine indices as tuple
                        arrayIndices = eval( Meta.parse( "(" * indices * ",)" ) )

                        # Extract sub-matrix
                        if isConstant
                            sig = getindex(sig2, arrayIndices...)
                        else
                            sig = [getindex(sig2[i], arrayIndices...) for i in eachindex(sig2)]
                        end
                        sigPresent = true
                    end
                end
            end
        end
    end

    if sigPresent
        return (sig, isConstant, arrayName, arrayIndices)
    else
        return (nothing, false, name, ())
    end
end

"""
    (signal, isConstant, arrayName, arrayIndices) = getSignalWithWarning(result, name)

Calls getSignal(result,name) and prints a warning message if `signal == nothing`
"""
function getSignalWithWarning(result,name)
    (signal, isConstant, arrayName, arrayIndices) = getSignal(result,name)
    if isnothing(signal)
        @warn "ModiaPlot: \"$name\" is not correct or is not defined or has no values."
    end
    return (signal, isConstant, arrayName, arrayIndices)
end

appendUnit2(name, unit) = unit == "" ? name : string(name, " [", unit, "]")

function appendUnit(name, value)
    if typeof(value) <: MonteCarloMeasurements.StaticParticles ||
       typeof(value) <: MonteCarloMeasurements.Particles
        appendUnit2(name, string(unit(value.particles[1])))
    else
        appendUnit2(name, string(unit(value)))
    end
end

"""
    (xsig, xsigLegend, ysig, ysigLegend, yIsConstant) = getPlotSignal(result, xsigName, ysigName)

Given the `AbstractString` key of the x-signal (`xsigName`) and of the y-signal (`ysigName`), returns

- `xsig`: the x-axis vector,
- `ysig`: the y-axis vector(s),
- `xsigLegend`: the legend of the x-axis vector as String,
- `ysigLegend`: the legend of the y-axis vector(s) as String vector,
- `yIsConstant`: whether ysig is a constant signal (= true) or not (=false).

If an element of `ysigName` is a scalar, then `ysig` is a vector
and `ysig[i]` is the value at time instant `i`.

If an element of `ysigName` is an array, then `ysig` is a matrix, where
ysig[i,:] are the array values at time instant `i`.
"""
function getPlotSignal(result, xsigName::AbstractString, ysigName::AbstractString)
    # Get x-axis signal
    (xsig, xIsConstant, dummy1, dummy2) = getSignalWithWarning(result, xsigName)
    if isnothing(xsig)
        @goto ERROR
    elseif xIsConstant
        @warn "ModiaPlot: \"$xsigName\" is a constant and this is not allowed for the x-axis variable."
        @goto ERROR
    elseif ndims(xsig[1]) != 0
        @warn "ModiaPlot: \"$xsigName\" does not characterize a scalar variable as needed for the x-axis."
        @goto ERROR
    elseif !( typeof(xsig[1]) <: Number )
        @warn "ModiaPlot: \"$xsigName\" has no Number type values, but values of type " * string(typeof(xsig[1])) * "."
        @goto ERROR
    elseif typeof(xsig[1]) <: Measurements.Measurement
        @warn "ModiaPlot: \"$xsigName\" is a Measurements.Measurement type and this is not supported for the x-axis."
        @goto ERROR
    elseif typeof(xsig[1]) <: MonteCarloMeasurements.StaticParticles
        @warn "ModiaPlot: \"$xsigName\" is a MonteCarloMeasurements.StaticParticles type and this is not supported for the x-axis."
        @goto ERROR
    elseif typeof(xsig[1]) <: MonteCarloMeasurements.Particles
        @warn "ModiaPlot: \"$xsigName\" is a MonteCarloMeasurements.Particles type and this is not supported for the x-axis."
        @goto ERROR
    end
    xsigLegend = appendUnit(xsigName, xsig[1])

    # Get y-axis signals
    (ysig, yIsConstant, yArrayName, yArrayIndices) = getSignalWithWarning(result, ysigName)
    if isnothing(ysig)
        @goto ERROR
    elseif !yIsConstant && length(ysig) != length(xsig)
        println("before warning message")
        @warn "ModiaPlot: \"$xsigName\" (= x-axis) has length " * string(length(xsig)) *
              " but \"$ysigName\" (= y-axis) has length " * string(length(ysig)) * "."
        @goto ERROR
    end

    value = yIsConstant ? ysig : ysig[1]

    if !(typeof(value) <: Number || typeof(value) <: AbstractArray)
        @warn "ModiaBase: \"$ysigName\" has no Number or AbstractArray values, but values of type " * string(typeof(value))
        @goto ERROR
    end

    if ndims(value) == 0
        # ysigName is a scalar variable
        ysigLegend = [appendUnit(ysigName, value)]
        if yIsConstant
            # Construct a constant time series with two points at the first and the last value of the time vector
            xsig = [xsig[1], xsig[end]]
            ysig = [ysig   , ysig     ]
        end

    else
        # ysigName is an array variable
        if yIsConstant
            xsig = [xsig[1], xsig[end]]
            ysig = reshape(ysig, :)
            ysig = collect([ysig ysig]')
        else
            nvalue = length(value)
            ysig2  = ysig
            ysig   = zeros(eltype(value), length(xsig), nvalue)
            for (i, value_i) in enumerate(ysig2)
                for j in 1:nvalue
                    ysig[i,j] = ysig2[i][j]
                end
            end

        end

        nLegend    = size(ysig,2)
        ysigLegend = [yArrayName * "[" for i = 1:nLegend]
        i = 1
        ySizeLength = Int[]
        for j1 in eachindex(yArrayIndices)
            push!(ySizeLength, length(yArrayIndices[j1]))
            i = 1
            if j1 == 1
                for j2 in 1:div(nLegend, ySizeLength[1])
                    for j3 in yArrayIndices[1]
                        ysigLegend[i] *= string(j3)
                        i += 1
                    end
                end
            else
                ncum = prod( ySizeLength[1:j1-1] )
                for j2 in yArrayIndices[j1]
                    for j3 = 1:ncum
                        ysigLegend[i] *= "," * string(j2)
                        i += 1
                    end
                end
            end
        end
        value = ysig[1]
        for i = 1:nLegend
            ysigLegend[i] *= appendUnit("]", value)
        end
    end

    return (ustrip.(xsig), xsigLegend, ustrip.(ysig), ysigLegend, yIsConstant)

    @label ERROR
    return (nothing, nothing, nothing, nothing, nothing)
end

getHeading(result, heading) = heading != "" ? heading : getDefaultHeading(result)

"""
    prepend!(prefix, ysigLegend)

Add `prefix` string in front of every element of the `ysigLegend` string-Vector.
"""
function prepend!(prefix::AbstractString, ysigLegend::Vector{AbstractString})
   for i in eachindex(ysigLegend)
      ysigLegend[i] = prefix*ysigLegend[i]
   end
   return ysigLegend
end

###################################################################################
#
# Exported functions
#
###################################################################################

"""
    printResultInfo(result)

Print names, sizes, units and element types of the signals that are stored in result.
"""
function printResultInfo(result)
    if isnothing(result)
        @info "The call of ModiaPlot.printResultInfo(result) is ignored, since the first argument is nothing."
        return
    end

    resultInfoTable = DataFrames.DataFrame(name=String[], unit=String[], constant=Bool[], varSize=String[], eltype=String[])

    nTime = -1
    first = true
    firstTimeSignal = ""
    for name in getNames(result)
        (isConstant, nTime2, varSize, elType, unit) = getSignalInfo(result, name)
        if !isConstant && !isnothing(varSize) && !ismissing(varSize)
            if first
                first = false
                nTime = nTime2
                firstTimeSignal = name
            elseif nTime != nTime2
                @warn "\"$name\" has $nTime2 time points but \"$firstTimeSignal\" has $nTime time points."
            end
        end
        push!(resultInfoTable, [name, string(unit), isConstant, string(varSize), string(elType)] )
    end

    # eltypes=false, truncate=50: Requires and older version of DataFrames
    show(stdout, resultInfoTable, summary=false, rowlabel=Symbol("#"), allcols=true)
    if nTime != -1
        println("\n... Number of time points = ", nTime)
    else
        println()
    end
end

"""
    plot(result, names;
         heading = "", grid = true, xAxis = "time",
         figure = 1, prefix = "", reuse = false, maxLegend = 10,
         minXaxisTickLabels = false,
         MonteCarloAsArea = true)

Generate **line plots** of the signals of the `result` data structure that are
identified with the `names` keys.

By default `result` is any type of dictionary with key type `::AbstractString` and value type `::Any`.
By providing own implementations of some access-functions, other `result` data structures can be used
(for details see xxx).

Argument `names` defines the diagrams to be drawn and the result data to be included in the respective diagram:

- If `names` is a **String**, generate one diagram with one time series of the variable with key `names`.

- If `names` is a **Tuple** of Strings, generate one diagram with the time series of the variables
  with the keys given in the tuple.

- If names is a **Vector** or a **Matrix** of **Strings** and/or **Tuples**,
  generate a vector or matrix of diagrams.

Note, the names (and their units, if available in the result) are automatically used as legends in the
respective diagram.

A signal variable identified by a `String` key can be a scalar of type `<:Number`
or an array of element type `<:Number`. A signal can be either a **constant** or a **time series**
where for every time instant one value of the signal variable
is stored in `result` (time instants are along the first dimension of an array).
In case of a constant, a constant line is drawn for the first to the last time
instant.

Note, before passing data to the plot package,
it is converted to Float64. This allows to, for example, also plot rational numbers,
even if not supported by the plot package. `Measurements.Measurement{xxx}` is specially handled.

# Optional Arguments
- `heading::AbstractString`: Optional heading above the diagram.

- `grid::Bool`: = true, to display a grid.

- `xAxis::AbstractString`: Name of x-axis.

- `figure::Int`: Integer identifier of the window in which the diagrams shall be drawn.

- `prefix::AbstractString`: String that is appended in front of every legend label
  (useful especially if `reuse=true`).

- `reuse::Bool`: If figure already exists and reuse=false, clear the figure before adding the plot.
   Otherwise, include the plot in the existing figure without removing the curves present in the figure.

- `maxLegend::Int`: If the number of legend entries in one plot command `> maxLegend`,
  the legend is suppressed.
  All curves have still their names as labels. In PyPlot, the curves can be inspected by
  their names by clicking in the toolbar of the plot on button `Edit axis, curve ..`
  and then on `Curves`.

- `minXaxisTickLabels::Bool`: = true, if xaxis tick labels shall be
  removed in a vector or array of plots, if not the last row
  (useful when including plots in a document).
  = false, x axis tick labels are always shown (useful when interactively zooming into a plot).

- `MonteCarloAsArea::Bool`: = true, if MonteCarloMeasurements values are shown with the mean value
  and the area between the minimum and the maximum value of all particles.
  = false, if all particles of MonteCarloMeasurements values are shown (e.g. if a value has 2000 particles,
  then 2000 curves are shown in the diagram).

# Examples
```julia
using ModiaPlot
using Unitful

# Construct result data
t = range(0.0, stop=10.0, length=100);
result = Dict{AbstractString,Any}();
result["time"] = t*u"s";
result["phi"]  = sin.(t)*u"rad";
result["w"]    = cos.(t)*u"rad/s";
result["a"]    = 1.2*sin.(t)*u"rad/s^2";
result["r"]    = hcat(0.4 * cos.(t), 0.5 * sin.(t), 0.3*cos.(t))*u"m";

# 1 signal in one diagram (legend = "phi [rad]")
plot(result, "phi")

# 3 signals in one diagram
plot(result, ("phi", "w", "a"), figure=2)

# 3 diagrams in form of a vector (every diagram has one signal)
plot(result, ["phi", "w", "r"], figure=3)

# 4 diagrams in form of a matrix (every diagram has one signal)
plot(result, ["phi" "w";
              "a"   "r[2]" ], figure=4)

# 2 diagrams in form of a vector
plot(result, [ ("phi", "w"), ("a") ], figure=5)

# 4 diagrams in form of a matrix
plot(result, [ ("phi",)           ("phi", "w");
               ("phi", "w", "a")  ("r[2:3]",)     ],figure=6)

# Plot w=f(phi) in one diagram
plot(result, "w", xAxis="phi", figure=7)

# Append signal of the next simulation run to figure=1
# (legend = "Sim 2: phi [rad]")
result["phi"] = 0.5*result["phi"];
plot(result, "phi", prefix="Sim 2: ", reuse=true)
```

Example of a matrix of plots:

![Matrix of plots](../resources/images/matrix-of-plots.png)
"""
plot(result, names::AbstractString; kwargs...) = plot(result, [names]        ; kwargs...)
plot(result, names::Symbol        ; kwargs...) = plot(result, [string(names)]; kwargs...)
plot(result, names::Tuple         ; kwargs...) = plot(result, [names]        ; kwargs...)
plot(result, names::AbstractVector; kwargs...) = plot(result, reshape(names, length(names), 1); kwargs...)


text:
using AddTwoDigit

try
    @assert add_two_digit(1.0,2.0) == 3.0
catch e
    throw(string(e.msg))
end