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

using Markdown
using InteractiveUtils

# ╔═╡ d09256dd-6c0d-4e28-9e54-3f7b3ca87ecb
begin
	using CSV, DataFrames
	using ScientificTypes
	using Tables
	using UrlDownload
	using SparseArrays
	using PlutoUI
end

# ╔═╡ 7e85780a-bd7e-4b5e-88ff-6053d64456f5
html"""

Tutorial - Part 1

Machine Learning in Julia

"""

# ╔═╡ d9536b16-bb2a-4f83-9666-9b9eeadc125b
PlutoUI.TableOfContents(title = "MLJ Tutorial - Part 1")

# ╔═╡ b67c62c2-15ab-4328-b795-033f6f2a0674
md"""
An introduction to the
[MLJ](https://alan-turing-institute.github.io/MLJ.jl/stable/)
toolbox.
"""

# ╔═╡ bc689638-fd19-4c9f-935f-ddf6a6bfbbdd
md"## Set-up"

# ╔═╡ 197fd00e-9068-46ea-af2a-25235e544a31
md"Inspect Julia version:"

# ╔═╡ f6d4f8c4-e441-45c4-8af5-148d95ea2900
VERSION

# ╔═╡ 45740c4d-b789-45dc-a6bf-47194d7e8e12
md"The following instantiates a package environment."

# ╔═╡ 42b0f1e1-16c9-4238-828a-4cc485149963
md"""
The package environment has been created using **Julia 1.6** and may not
instantiate properly for other Julia versions.
"""

# ╔═╡ 499cbc31-83ba-4583-ba1f-6363f43ec697
md"## General resources"

# ╔═╡ db5821e8-956a-4a46-95ea-2955abd45275
md"""
- [MLJ Cheatsheet](https://alan-turing-institute.github.io/MLJ.jl/dev/mlj_cheatsheet/)
- [Common MLJ Workflows](https://alan-turing-institute.github.io/MLJ.jl/dev/common_mlj_workflows/)
- [MLJ manual](https://alan-turing-institute.github.io/MLJ.jl/dev/)
- [Data Science Tutorials in Julia](https://juliaai.github.io/DataScienceTutorials.jl/)
"""

# ╔═╡ a9cee8a6-dcfb-4d97-b1b3-dd47e367ba54
md"# Part 1 - Data Representation"

# ╔═╡ da747378-24c9-4b9c-8d7e-8e449afd1c48
md"""
> **Goals:**
> 1. Learn how MLJ specifies it's data requirements using "scientific" types
> 2. Understand the options for representing tabular data
> 3. Learn how to inspect and fix the representation of data to meet MLJ requirements
"""

# ╔═╡ fe0ead17-c3f5-4f3c-a949-7f3bd292fe31
md"## Scientific types"

# ╔═╡ 4ff58cc5-ea99-4e9d-8514-99938a2932db
md"""
To help you focus on the intended *purpose* or *interpretation* of
data, MLJ models specify data requirements using *scientific types*,
instead of machine types. An example of a scientific type is
`OrderedFactor`. The other basic "scalar" scientific types are
illustrated below:
"""

# ╔═╡ 675753b1-786e-4344-a0de-1721c1bc2df2
html"""

	<img src="https://github.com/ablaom/MLJTutorial.jl/blob/dev/notebooks/01_data_representation/scitypes.png?raw=true"

"""

# ╔═╡ 74e167cd-bc94-486c-bca8-7eec7950fd82
md"""
A scientific type is an ordinary Julia type (so it can be used for
method dispatch, for example) but it usually has no instances. The
`scitype` function is used to articulate MLJ's convention about how
different machine types will be interpreted by MLJ models:
"""

# ╔═╡ 45c7e17a-63c8-4fb3-9873-1b1430e635cc
scitype(3.141)

# ╔═╡ 1a742c64-61e8-4fbf-b43e-d0ebe87da176
time = [2.3, 4.5, 4.2, 1.8, 7.1]

# ╔═╡ 44737448-20e7-464b-bb15-33279dbe6c3a
scitype(time)

# ╔═╡ f4ec6544-bdd3-4c58-9009-4cfb2012998f
md"""
To fix data which MLJ is interpreting incorrectly, we use the
`coerce` method:
"""

# ╔═╡ 5e3fc1a4-2e74-48c6-9864-152c603930ec
height0 = [185, 153, 163, 114, 180]

# ╔═╡ 6ec50b8f-9c57-4085-abd3-cb77d7a79683
scitype(height0)

# ╔═╡ 8e3fc777-327a-46e3-879e-dc128f3db7b3
height = coerce(height0, Continuous)

# ╔═╡ 91abff52-6801-423d-a367-9aa3c6cf34d0
md"""
Here's an example of data we would want interpreted as
`OrderedFactor` but isn't:
"""

# ╔═╡ 8d7b3ba5-a840-4377-bf32-94657e65284e
exam_mark0 = ["rotten", "great", "bla",  missing, "great"]

# ╔═╡ 7482fffc-7b40-4e84-9589-873d4bc96cca
scitype(exam_mark0)

# ╔═╡ 7424426a-cdab-4b1e-9afd-65a5b5facac9
exam_mark = coerce(exam_mark0, OrderedFactor)

# ╔═╡ dd8cef99-8f55-47a8-b6c8-54f5ed90a964
levels(exam_mark)

# ╔═╡ 41012c9e-8e71-4a9d-9292-fe822525fc77
md"Use `levels!` to put the classes in the right order:"

# ╔═╡ c1c65dab-f82a-400f-ae5e-4b83dcbc9675
levels!(exam_mark, ["rotten", "bla", "great"])

# ╔═╡ 0033787e-0796-41c6-a81d-fc52d08ea8e0
exam_mark[1] < exam_mark[2]

# ╔═╡ e33a8562-1689-43ad-8a29-32f11452654a
md"When sub-sampling, no levels are lost:"

# ╔═╡ c7135c01-4881-45f6-a5f1-151fcbe229a2
levels(exam_mark[1:2])

# ╔═╡ 0e1afb98-7b11-4eb8-81bc-2d5603785a09
md"""
**Note on binary data.** There is no separate scientific type for
binary data. Binary data is `OrderedFactor{2}` or
`Multiclass{2}`. If a binary measure like `truepositive` is a
applied to `OrderedFactor{2}` then the "positive" class is assumed
to appear *second* in the ordering. If such a measure is applied to
`Multiclass{2}` data, a warning is issued. A single `OrderedFactor`
can be coerced to a single `Continuous` variable, for models that
require this, while a `Multiclass` variable can only be one-hot
encoded.
"""

# ╔═╡ ffcd5bb0-b53e-4ac8-9d87-1d1fbb0e3997
md"## Two-dimensional data"

# ╔═╡ f4a155dc-f3f2-4c44-b950-fddef2a1fb10
md"""
Whenever it makes sense, MLJ Models generally expect two-dimensional
data to be *tabular*. All the tabular formats implementing the
[Tables.jl API](https://juliadata.github.io/Tables.jl/stable/) (see
this
[list](https://github.com/JuliaData/Tables.jl/blob/master/INTEGRATIONS.md))
have a scientific type of `Table` and can be used with such models.
"""

# ╔═╡ 308f453a-767e-4742-bab9-13a4b5e22567
md"""
### Column table
"""

# ╔═╡ 09cd28d1-eabf-488a-951b-64cb2a36c8e3
md"""
Probably the simplest example of a table is the julia native *column
table*, which is just a named tuple of equal-length vectors:
"""

# ╔═╡ 9428e4a9-8e51-4d13-b0e6-441461cc8770
column_table = (h=height, e=exam_mark, t=time)

# ╔═╡ b0d23945-5f6d-489e-8cb0-de601961bc02
scitype(column_table)

# ╔═╡ ed5001c5-e5b9-4a84-a87d-950c50fb2955
md"""
Notice the `Table{K}` type parameter `K` encodes the scientific
types of the columns. (This is useful when comparing table scitypes
with `<:`). To inspect the individual column scitypes, we use the
`schema` method instead:
"""

# ╔═╡ 91c12549-e518-47d3-8448-4bcb089096cd
sc = schema(column_table)

# ╔═╡ 62be5b63-a6a3-48b3-b23d-a60990b06fec
md"""
Specific information about a schema can be obtained by accessing the attributes `names`, `types` and `scitypes` directly:
"""

# ╔═╡ 4eaac509-d637-41bf-85b6-a8154be0f959
sc.names

# ╔═╡ fef11f2a-da56-430d-b56d-ac411bbca555
sc.types

# ╔═╡ c37d5a96-29a7-465e-9ed4-cd7542e9065a
sc.scitypes

# ╔═╡ 1257966c-8043-49de-a012-a84240254fb6
md"Here are five other examples of tables:"

# ╔═╡ 27415a96-42ec-4dc9-9428-713a24ceeda9
md"""
### Dictionary
"""

# ╔═╡ 3b776076-56c6-4215-bbdd-9799f7bb2e60
dict_table = Dict(:h => height, :e => exam_mark, :t => time)

# ╔═╡ 0cf291d4-4e7d-4b28-8562-56d54658ca00
schema(dict_table)

# ╔═╡ 919fcf78-227b-43c4-97a7-88e3af4f10ee
md"""
(To control column order here, instead use `LittleDict` from
OrderedCollections.jl.)
"""

# ╔═╡ e72b4214-b18f-40ba-a0c0-a78f67d659a0
md"""
### Row Table
"""

# ╔═╡ aa98759a-f43d-491e-b372-8c9866e51852
row_table = [(a=1, b=3.4),
             (a=2, b=4.5),
             (a=3, b=5.6)]

# ╔═╡ cfb975ac-d946-4204-8427-bccb0a0e1ae7
schema(row_table)

# ╔═╡ eba6d196-f627-49be-8f80-51fef8c5aa59
md"""
### DataFrame
"""

# ╔═╡ 74321d51-1831-4822-8f3d-b82f1e7b6f7a
df = DataFrames.DataFrame(column_table)

# ╔═╡ 6971b1f8-d607-42ae-ab04-579e5608ae53
schema(df) == schema(column_table)

# ╔═╡ ee1c9b0f-62f5-4dfb-89b8-6a303aca22b1
md"""
### CSV File
"""

# ╔═╡ fa257731-0521-4ffe-86cf-35420d9e6995
csv_file = urldownload("https://raw.githubusercontent.com/ablaom/"*
                     "MachineLearningInJulia2020/"*
                     "for-MLJ-version-0.16/data/horse.csv")

# ╔═╡ e6f680e3-3dd7-421c-bb76-0ab026336d20
 schema(csv_file)

# ╔═╡ 75f109af-7651-477a-a6f7-2e3eb635a58d
md"""
### Matrix
"""

# ╔═╡ 8134229e-c93c-482a-a29a-8d0445351914
md"""
Most MLJ models do not accept matrix in lieu of a table, but you can
wrap a matrix as a table:
"""

# ╔═╡ 6c7d2c63-9bec-4439-be65-28d0fcca50a8
matrix_table = Tables.table(rand(2,3))

# ╔═╡ ac973148-6168-473a-adb7-3880e5bd45b6
schema(matrix_table)

# ╔═╡ 5584fdb0-7a69-4f1f-9a2e-c468345d87d1
md"""
The matrix is *not* copied, only wrapped. Some models may perform
better if one wraps the adjoint of the transpose - see
[here](https://alan-turing-institute.github.io/MLJ.jl/dev/getting_started/#Observations-correspond-to-rows,-not-columns).
"""

# ╔═╡ 7cbaa6ba-e85a-4eba-b5fa-cb79ebf595ef
md"""
**Manipulating tabular data.** In this workshop we assume
familiarity with some kind of tabular data container (although it is
possible, in principle, to carry out the exercises without this.)
For a quick start introduction to `DataFrames`, see [this
tutorial](https://juliaai.github.io/DataScienceTutorials.jl/data/dataframe/).
"""

# ╔═╡ 9518f91d-655c-4e3c-91c5-524da38aa391
md"## Fixing scientific types in tabular data"

# ╔═╡ 03df98f9-b52b-4690-ad90-b4405b216771
md"""
To show how we can correct the scientific types of data in tables,
we introduce a cleaned up version of the UCI Horse Colic Data Set
(the cleaning work-flow is described
[here](https://juliaai.github.io/DataScienceTutorials.jl/end-to-end/horse/#dealing_with_missing_values)).
We already downloaded this data set immediately above.q
"""

# ╔═╡ 6dd3ebb2-cb17-4d0a-895b-997a92b731e0
horse = DataFrames.DataFrame(csv_file) # convert to data frame

# ╔═╡ d11f7c90-acf2-421b-a525-e9b04a4bd246
md"""
From [the UCI
docs](http://archive.ics.uci.edu/ml/datasets/Horse+Colic) we can
surmise how each variable ought to be interpreted (a step in our
work-flow that cannot reliably be left to the computer):
"""

# ╔═╡ 604f6d68-7c31-468c-80f0-ba7781e173cf
md"""
variable                    | scientific type (interpretation)
----------------------------|-----------------------------------
`:surgery`                  | Multiclass
`:age`                      | Multiclass
`:rectal_temperature`       | Continuous
`:pulse`                    | Continuous
`:respiratory_rate`         | Continuous
`:temperature_extremities`  | OrderedFactor
`:mucous_membranes`         | Multiclass
`:capillary_refill_time`    | Multiclass
`:pain`                     | OrderedFactor
`:peristalsis`              | OrderedFactor
`:abdominal_distension`     | OrderedFactor
`:packed_cell_volume`       | Continuous
`:total_protein`            | Continuous
`:outcome`                  | Multiclass
`:surgical_lesion`          | OrderedFactor
`:cp_data`                  | Multiclass
"""

# ╔═╡ ea376156-57e3-43e1-9cbc-9a2c39793333
md"""
Let's see how MLJ will actually interpret the data, as it is
currently encoded:
"""

# ╔═╡ f878b48f-bc60-4e3c-b887-2b16710e5502
schema(horse)

# ╔═╡ 6d08247f-b9b8-456d-9452-5298a8a4a401
md"""
As a first correction step, we can get MLJ to "guess" the
appropriate fix, using the `autotype` method:
"""

# ╔═╡ 181aad5d-e272-49ec-b017-96f1602f2cad
autotype(horse)

# ╔═╡ f8e9f01f-74f6-4612-8be2-2a4017c45345
md"""
Okay, this is not perfect, but a step in the right direction, which
we implement like this:
"""

# ╔═╡ 43498690-c1cb-4dd4-a7ad-80d84f5b0070
begin
  coerce!(horse, autotype(horse));
  schema(horse)
end

# ╔═╡ 22c0ce1a-6a07-4b7a-8378-5aa686f0b407
md"All remaining `Count` data should be `Continuous`:"

# ╔═╡ 69a45801-0a1c-4f47-9f43-49763e8691a1
begin
  coerce!(horse, Count => Continuous);
  schema(horse)
end

# ╔═╡ f4734743-a383-45a2-bb0e-3a45761c733a
md"We'll correct the remaining truant entries manually:"

# ╔═╡ 2e0c0c2b-b011-48f4-9979-7044b2f2df33
begin
  coerce!(horse,
          :surgery               => Multiclass,
          :age                   => Multiclass,
          :mucous_membranes      => Multiclass,
          :capillary_refill_time => Multiclass,
          :outcome               => Multiclass,
          :cp_data               => Multiclass);
  schema(horse)
end

# ╔═╡ 1fc47a9c-b126-489d-b2a1-fbbde543f620
md"""
# Resources for Part 1

- From the MLJ manual:
   - [A preview of data type specification in
  MLJ](https://alan-turing-institute.github.io/MLJ.jl/dev/getting_started/#A-preview-of-data-type-specification-in-MLJ-1)
   - [Data containers and scientific types](https://alan-turing-institute.github.io/MLJ.jl/dev/getting_started/#Data-containers-and-scientific-types-1)
   - [Working with Categorical Data](https://alan-turing-institute.github.io/MLJ.jl/dev/working_with_categorical_data/)
- [Summary](https://juliaai.github.io/ScientificTypes.jl/dev/#Summary-of-the-default-convention) of the MLJ convention for representing scientific types
- [ScientificTypes.jl](https://juliaai.github.io/ScientificTypes.jl/dev/)
- From Data Science Tutorials:
    - [Data interpretation: Scientific Types](https://juliaai.github.io/DataScienceTutorials.jl/data/scitype/)
    - [Horse colic data](https://juliaai.github.io/DataScienceTutorials.jl/end-to-end/horse/)
- [UCI Horse Colic Data Set](http://archive.ics.uci.edu/ml/datasets/Horse+Colic)
"""

# ╔═╡ c9ac414c-ca43-460f-910b-f8a9a217eff2
md"# Exercises for Part 1"

# ╔═╡ 84a75325-e427-45c5-aa36-d2e8546da46a
md"## Exercise 1"

# ╔═╡ 1c47cc40-33f8-4552-88a1-a2bf91434413
md"""
Try to guess how each code snippet below will evaluate (uncomment to see the solutions):
"""

# ╔═╡ 4da068ed-3935-4728-a1cb-54e34396a855
# scitype(42)

# ╔═╡ 04c88aba-6c17-48eb-8037-0de5806e1a54
questions = ["who", "why", "what", "when"]

# ╔═╡ e5c60366-e03c-490f-9475-e53e5b3347d4
# scitype(questions)

# ╔═╡ dc8c4bd5-e753-4b7b-9961-a09632c33fb0
# elscitype(questions)

# ╔═╡ 34c94131-ef2a-4f86-b7cc-46f4ef988c68
t = (3.141, 42, "how")

# ╔═╡ a9238f1a-2751-46a8-9d69-62669915ed91
# scitype(t)

# ╔═╡ 8f5681ef-5432-4daf-90fb-7f2721f6fcc7
A = rand(2, 3)

# ╔═╡ 281afa2b-c539-4c61-8890-979691212d9b
#scitype(A)

# ╔═╡ 67054d46-550d-40e6-a6fa-eaf14df5d44b
#elscitype(A)

# ╔═╡ ec84d270-fcb9-4b4b-8025-5084804a9f6c
Asparse = sparse(A)

# ╔═╡ d3c65fbf-1a6c-4779-9e90-09023d201062
#scitype(Asparse)

# ╔═╡ 1c530c2e-f0b8-4e30-b7bb-2f33ef765cc0
C = coerce(A, Multiclass)

# ╔═╡ 09b804b7-30f7-4831-9625-01172c4a0081
#scitype(C)

# ╔═╡ 1e774630-8eac-4394-af4f-11c7de9e21dd
#elscitype(C)

# ╔═╡ 347d35af-60d3-4422-8dba-241d9b744683
v = [1, 2, missing, 4]

# ╔═╡ 3aac62ed-9e28-4db9-a172-d57784b35656
#scitype(v)

# ╔═╡ ed2aa555-53d0-4f61-a6e5-2f0b4dca5c59
#elscitype(v)

# ╔═╡ 28a31668-333e-434a-8550-20498aa03ed0
#scitype(v[1:2])

# ╔═╡ 2cc4182e-ffaa-47a6-9e7b-5b943cf6b560
md"""
Can you guess at the general behavior of
`scitype` with respect to tuples, abstract arrays and missing
values? The answers are
[here](https://github.com/juliaai/ScientificTypesBase.jl#2-the-scitype-and-scitype-methods)
(ignore "Property 1").
"""

# ╔═╡ 34a20115-36f5-47d5-bce6-0a1379cc1259
md"## Exercise 2"

# ╔═╡ b5866134-32b9-4ad9-9608-b5032c116833
md"""
Coerce the following vector to make MLJ recognize it as a vector of
ordered factors (with an appropriate ordering):
"""

# ╔═╡ 7f9003b5-95c7-418b-b473-4aa968e6937a
quality0 = ["good", "poor", "poor", "excellent", missing, "good", "excellent"]

# ╔═╡ 58ed39f5-3c80-460d-8d9e-644a1b3cc6b6
md"## Exercise 3 (fixing scitypes in a table)"

# ╔═╡ c8162baa-74b7-4046-ac0a-23ded81445da
md"""
Fix the scitypes for the [House Prices in King
County](https://mlr3gallery.mlr-org.com/posts/2020-01-30-house-prices-in-king-county/)
dataset:
"""

# ╔═╡ 2e797d14-c00b-4166-8535-1a088a6a3228
begin
  house_csv = urldownload("https://raw.githubusercontent.com/ablaom/"*
                          "MachineLearningInJulia2020/for-MLJ-version-0.16/"*
                          "data/house.csv");
  house = DataFrames.DataFrame(house_csv)
  first(house, 4)
end

# ╔═╡ 706005cd-a3e5-4b5b-a39f-7893c73eef39
md"""
(Two features in the original data set have been deemed uninformative
and dropped, namely `:id` and `:date`. The original feature
`:yr_renovated` has been replaced by the `Bool` feature `is_renovated`.)
"""

# ╔═╡ 7d15cc45-23a6-4a3d-bcca-51a27994a151
html"""

"""

# ╔═╡ 135dac9b-0bd9-4e1d-9b33-09b4b6661170
md"""
---

*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
ScientificTypes = "321657f4-b219-11e9-178b-2701a2544e81"
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
UrlDownload = "856ac37a-3032-4c1c-9122-f86d88358c8b"

[compat]
CSV = "~0.9.11"
DataFrames = "~1.3.1"
PlutoUI = "~0.7.27"
ScientificTypes = "~3.0.0"
Tables = "~1.6.1"
UrlDownload = "~1.0.0"
"""

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

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

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

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

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

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

[[CategoricalArrays]]
deps = ["DataAPI", "Future", "Missings", "Printf", "Requires", "Statistics", "Unicode"]
git-tree-sha1 = "c308f209870fdbd84cb20332b6dfaf14bf3387f8"
uuid = "324d7699-5711-5eae-9e2f-1d82baa6b597"
version = "0.10.2"

[[ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "d711603452231bad418bd5e0c91f1abd650cba71"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.11.3"

[[ChangesOfVariables]]
deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
git-tree-sha1 = "bf98fa45a0a4cee295de98d4c1462be26345b9a1"
uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
version = "0.1.2"

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

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

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

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

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

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

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

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

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

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

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

[[DensityInterface]]
deps = ["InverseFunctions", "Test"]
git-tree-sha1 = "80c3e8639e3353e5d2912fb3a1916b8455e2494b"
uuid = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
version = "0.4.0"

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

[[Distributions]]
deps = ["ChainRulesCore", "DensityInterface", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns", "Test"]
git-tree-sha1 = "6a8dc9f82e5ce28279b6e3e2cea9421154f5bd0d"
uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
version = "0.25.37"

[[DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

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

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

[[FillArrays]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
git-tree-sha1 = "8756f9935b7ccc9064c6eef0bff0ad643df733a3"
uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
version = "0.12.7"

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

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

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

[[HTTP]]
deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
git-tree-sha1 = "0fa77022fe4b511826b39c894c90daf5fce3334a"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.9.17"

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

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

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

[[IniFile]]
deps = ["Test"]
git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.0"

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

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

[[InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "a7254c0acd8e62f1ac75ad24d5db43f5f19f3c65"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.2"

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

[[IrrationalConstants]]
git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.1"

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

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

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

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

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

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

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

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

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

[[LogExpFunctions]]
deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "e5718a00af0ab9756305a0392832c8952c7426c1"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.6"

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

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

[[MbedTLS]]
deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
version = "1.0.3"

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

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

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

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

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

[[OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"

[[OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

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

[[PDMats]]
deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
git-tree-sha1 = "ee26b350276c51697c9c2d88a072b339f9f03d73"
uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
version = "0.11.5"

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

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

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

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

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

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

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

[[ProgressMeter]]
deps = ["Distributed", "Printf"]
git-tree-sha1 = "afadeba63d90ff223a6a48d2009434ecee2ec9e8"
uuid = "92933f4c-e287-5a05-a399-4b506db050ca"
version = "1.7.1"

[[QuadGK]]
deps = ["DataStructures", "LinearAlgebra"]
git-tree-sha1 = "78aadffb3efd2155af139781b8a8df1ef279ea39"
uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
version = "2.4.2"

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

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

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

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "8f82019e525f4d5c669692772a6f4b0a58b06a6a"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.2.0"

[[Rmath]]
deps = ["Random", "Rmath_jll"]
git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f"
uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
version = "0.7.0"

[[Rmath_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7"
uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
version = "0.3.0+0"

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

[[ScientificTypes]]
deps = ["CategoricalArrays", "ColorTypes", "Dates", "Distributions", "PrettyTables", "Reexport", "ScientificTypesBase", "StatisticalTraits", "Tables"]
git-tree-sha1 = "ba70c9a6e4c81cc3634e3e80bb8163ab5ef57eb8"
uuid = "321657f4-b219-11e9-178b-2701a2544e81"
version = "3.0.0"

[[ScientificTypesBase]]
git-tree-sha1 = "a8e18eb383b5ecf1b5e6fc237eb39255044fd92b"
uuid = "30f210dd-8aff-4c5f-94ba-8e64358c1161"
version = "3.0.0"

[[SentinelArrays]]
deps = ["Dates", "Random"]
git-tree-sha1 = "244586bc07462d22aed0113af9c731f2a518c93e"
uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
version = "1.3.10"

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

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

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

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

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

[[SpecialFunctions]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "e08890d19787ec25029113e88c34ec20cac1c91e"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "2.0.0"

[[StatisticalTraits]]
deps = ["ScientificTypesBase"]
git-tree-sha1 = "271a7fea12d319f23d55b785c51f6876aadb9ac0"
uuid = "64bff920-2084-43da-a3e6-9bb72801c0c9"
version = "3.0.0"

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

[[StatsAPI]]
git-tree-sha1 = "d88665adc9bcf45903013af0982e2fd05ae3d0a6"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.2.0"

[[StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "51383f2d367eb3b444c961d485c565e4c0cf4ba0"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.14"

[[StatsFuns]]
deps = ["ChainRulesCore", "InverseFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
git-tree-sha1 = "bedb3e17cc1d94ce0e6e66d3afa47157978ba404"
uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
version = "0.9.14"

[[SuiteSparse]]
deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"

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

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

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

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

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

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

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

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

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

[[UrlDownload]]
deps = ["HTTP", "ProgressMeter"]
git-tree-sha1 = "05f86730c7a53c9da603bd506a4fc9ad0851171c"
uuid = "856ac37a-3032-4c1c-9122-f86d88358c8b"
version = "1.0.0"

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

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

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

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

# ╔═╡ Cell order:
# ╟─7e85780a-bd7e-4b5e-88ff-6053d64456f5
# ╟─d9536b16-bb2a-4f83-9666-9b9eeadc125b
# ╟─b67c62c2-15ab-4328-b795-033f6f2a0674
# ╟─bc689638-fd19-4c9f-935f-ddf6a6bfbbdd
# ╟─197fd00e-9068-46ea-af2a-25235e544a31
# ╠═f6d4f8c4-e441-45c4-8af5-148d95ea2900
# ╟─45740c4d-b789-45dc-a6bf-47194d7e8e12
# ╟─42b0f1e1-16c9-4238-828a-4cc485149963
# ╠═d09256dd-6c0d-4e28-9e54-3f7b3ca87ecb
# ╟─499cbc31-83ba-4583-ba1f-6363f43ec697
# ╟─db5821e8-956a-4a46-95ea-2955abd45275
# ╟─a9cee8a6-dcfb-4d97-b1b3-dd47e367ba54
# ╟─da747378-24c9-4b9c-8d7e-8e449afd1c48
# ╟─fe0ead17-c3f5-4f3c-a949-7f3bd292fe31
# ╟─4ff58cc5-ea99-4e9d-8514-99938a2932db
# ╟─675753b1-786e-4344-a0de-1721c1bc2df2
# ╟─74e167cd-bc94-486c-bca8-7eec7950fd82
# ╠═45c7e17a-63c8-4fb3-9873-1b1430e635cc
# ╠═1a742c64-61e8-4fbf-b43e-d0ebe87da176
# ╠═44737448-20e7-464b-bb15-33279dbe6c3a
# ╟─f4ec6544-bdd3-4c58-9009-4cfb2012998f
# ╠═5e3fc1a4-2e74-48c6-9864-152c603930ec
# ╠═6ec50b8f-9c57-4085-abd3-cb77d7a79683
# ╠═8e3fc777-327a-46e3-879e-dc128f3db7b3
# ╟─91abff52-6801-423d-a367-9aa3c6cf34d0
# ╠═8d7b3ba5-a840-4377-bf32-94657e65284e
# ╠═7482fffc-7b40-4e84-9589-873d4bc96cca
# ╠═7424426a-cdab-4b1e-9afd-65a5b5facac9
# ╠═dd8cef99-8f55-47a8-b6c8-54f5ed90a964
# ╟─41012c9e-8e71-4a9d-9292-fe822525fc77
# ╠═c1c65dab-f82a-400f-ae5e-4b83dcbc9675
# ╠═0033787e-0796-41c6-a81d-fc52d08ea8e0
# ╟─e33a8562-1689-43ad-8a29-32f11452654a
# ╠═c7135c01-4881-45f6-a5f1-151fcbe229a2
# ╟─0e1afb98-7b11-4eb8-81bc-2d5603785a09
# ╟─ffcd5bb0-b53e-4ac8-9d87-1d1fbb0e3997
# ╟─f4a155dc-f3f2-4c44-b950-fddef2a1fb10
# ╟─308f453a-767e-4742-bab9-13a4b5e22567
# ╟─09cd28d1-eabf-488a-951b-64cb2a36c8e3
# ╠═9428e4a9-8e51-4d13-b0e6-441461cc8770
# ╠═b0d23945-5f6d-489e-8cb0-de601961bc02
# ╟─ed5001c5-e5b9-4a84-a87d-950c50fb2955
# ╠═91c12549-e518-47d3-8448-4bcb089096cd
# ╟─62be5b63-a6a3-48b3-b23d-a60990b06fec
# ╠═4eaac509-d637-41bf-85b6-a8154be0f959
# ╠═fef11f2a-da56-430d-b56d-ac411bbca555
# ╠═c37d5a96-29a7-465e-9ed4-cd7542e9065a
# ╟─1257966c-8043-49de-a012-a84240254fb6
# ╟─27415a96-42ec-4dc9-9428-713a24ceeda9
# ╠═3b776076-56c6-4215-bbdd-9799f7bb2e60
# ╠═0cf291d4-4e7d-4b28-8562-56d54658ca00
# ╟─919fcf78-227b-43c4-97a7-88e3af4f10ee
# ╟─e72b4214-b18f-40ba-a0c0-a78f67d659a0
# ╠═aa98759a-f43d-491e-b372-8c9866e51852
# ╠═cfb975ac-d946-4204-8427-bccb0a0e1ae7
# ╟─eba6d196-f627-49be-8f80-51fef8c5aa59
# ╠═74321d51-1831-4822-8f3d-b82f1e7b6f7a
# ╠═6971b1f8-d607-42ae-ab04-579e5608ae53
# ╟─ee1c9b0f-62f5-4dfb-89b8-6a303aca22b1
# ╠═fa257731-0521-4ffe-86cf-35420d9e6995
# ╠═e6f680e3-3dd7-421c-bb76-0ab026336d20
# ╟─75f109af-7651-477a-a6f7-2e3eb635a58d
# ╟─8134229e-c93c-482a-a29a-8d0445351914
# ╠═6c7d2c63-9bec-4439-be65-28d0fcca50a8
# ╠═ac973148-6168-473a-adb7-3880e5bd45b6
# ╟─5584fdb0-7a69-4f1f-9a2e-c468345d87d1
# ╟─7cbaa6ba-e85a-4eba-b5fa-cb79ebf595ef
# ╟─9518f91d-655c-4e3c-91c5-524da38aa391
# ╟─03df98f9-b52b-4690-ad90-b4405b216771
# ╠═6dd3ebb2-cb17-4d0a-895b-997a92b731e0
# ╟─d11f7c90-acf2-421b-a525-e9b04a4bd246
# ╟─604f6d68-7c31-468c-80f0-ba7781e173cf
# ╟─ea376156-57e3-43e1-9cbc-9a2c39793333
# ╠═f878b48f-bc60-4e3c-b887-2b16710e5502
# ╟─6d08247f-b9b8-456d-9452-5298a8a4a401
# ╠═181aad5d-e272-49ec-b017-96f1602f2cad
# ╟─f8e9f01f-74f6-4612-8be2-2a4017c45345
# ╠═43498690-c1cb-4dd4-a7ad-80d84f5b0070
# ╟─22c0ce1a-6a07-4b7a-8378-5aa686f0b407
# ╠═69a45801-0a1c-4f47-9f43-49763e8691a1
# ╟─f4734743-a383-45a2-bb0e-3a45761c733a
# ╠═2e0c0c2b-b011-48f4-9979-7044b2f2df33
# ╟─1fc47a9c-b126-489d-b2a1-fbbde543f620
# ╟─c9ac414c-ca43-460f-910b-f8a9a217eff2
# ╟─84a75325-e427-45c5-aa36-d2e8546da46a
# ╟─1c47cc40-33f8-4552-88a1-a2bf91434413
# ╠═4da068ed-3935-4728-a1cb-54e34396a855
# ╠═04c88aba-6c17-48eb-8037-0de5806e1a54
# ╠═e5c60366-e03c-490f-9475-e53e5b3347d4
# ╠═dc8c4bd5-e753-4b7b-9961-a09632c33fb0
# ╠═34c94131-ef2a-4f86-b7cc-46f4ef988c68
# ╠═a9238f1a-2751-46a8-9d69-62669915ed91
# ╠═8f5681ef-5432-4daf-90fb-7f2721f6fcc7
# ╠═281afa2b-c539-4c61-8890-979691212d9b
# ╠═67054d46-550d-40e6-a6fa-eaf14df5d44b
# ╠═ec84d270-fcb9-4b4b-8025-5084804a9f6c
# ╠═d3c65fbf-1a6c-4779-9e90-09023d201062
# ╠═1c530c2e-f0b8-4e30-b7bb-2f33ef765cc0
# ╠═09b804b7-30f7-4831-9625-01172c4a0081
# ╠═1e774630-8eac-4394-af4f-11c7de9e21dd
# ╠═347d35af-60d3-4422-8dba-241d9b744683
# ╠═3aac62ed-9e28-4db9-a172-d57784b35656
# ╠═ed2aa555-53d0-4f61-a6e5-2f0b4dca5c59
# ╠═28a31668-333e-434a-8550-20498aa03ed0
# ╟─2cc4182e-ffaa-47a6-9e7b-5b943cf6b560
# ╟─34a20115-36f5-47d5-bce6-0a1379cc1259
# ╟─b5866134-32b9-4ad9-9608-b5032c116833
# ╠═7f9003b5-95c7-418b-b473-4aa968e6937a
# ╟─58ed39f5-3c80-460d-8d9e-644a1b3cc6b6
# ╟─c8162baa-74b7-4046-ac0a-23ded81445da
# ╠═2e797d14-c00b-4166-8535-1a088a6a3228
# ╟─706005cd-a3e5-4b5b-a39f-7893c73eef39
# ╟─7d15cc45-23a6-4a3d-bcca-51a27994a151
# ╟─135dac9b-0bd9-4e1d-9b33-09b4b6661170
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
################################################################################
#
#  Field access
#
################################################################################

order(u::UnitGrpCtx) = u.order

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, U::UnitGrpCtx)
  print(io, "Unit group context of\n$(order(U))\n")
end
################################################################################
#
#  Initialization
#
################################################################################

function _unit_group_init(O::NfOrd)
  u = UnitGrpCtx{FacElem{nf_elem, AnticNumberField}}(O)
  return u
end

function _add_dependent_unit(U::UnitGrpCtx{S}, y::T; rel_only = false) where {S, T}
  K = nf(order(U))
  deg = degree(K)
  r1, r2 = signature(K)
  rr = r1 + r2
  r = rr - 1 # unit rank

  p = _rel_add_prec(U)

  #p = 64

  zz = Array{fmpz}(r + 1)

  @vprint :UnitGroup 1 "Adding dependent unit ... \n"

  @v_do :UnitGroup 1 pushindent()
  p, B = _conj_log_mat_cutoff_inv(U, p)
  @v_do :UnitGroup 1 popindent()
  @vprint :UnitGroup 2 "Precision is now $p\n"

  Ar = base_ring(B)

  b = ArbMatSpace(Ar, 1, r, false)()

  conlog = conjugates_arb_log(y, p)

  for i in 1:r
    b[1, i] = conlog[i]
  end

  # I should do this using lu decomposition and caching
  # The inversion could go wrong,
  # Then we again have to increase the precision

  inv_succesful = true

  @vprint :UnitGroup 3 "For $p element b: $b\n"
  v = b*B
  @vprint :UnitGroup 3 "For $p the vector v: $v\n"

  z = Array{fmpq}(r)

  rreg = arb()

  if isdefined(U, :tentative_regulator)
    rreg = U.tentative_regulator
  else
    @vprint :UnitGroup 1 "Computing regulator of independent units with 64 bits ... \n"
    rreg = regulator(U.units, 64)
    @vprint :UnitGroup 1 "done \n"
  end

  bound = _denominator_bound_in_relation(rreg, K)

  rel = Array{fmpz}(r + 1)
  for i in 1:r+1
    rel[i] = zero(FlintZZ)
  end

  @vprint :UnitGroup 2 "First iteration to find a rational relation ... \n"
  while !_find_rational_relation!(rel, v, bound)
    @vprint :UnitGroup 2 "Precision not high enough, increasing from $p to $(2*p)\n"
    p =  2*p
    p, B = _conj_log_mat_cutoff_inv(U, p)

    @assert p > 0

    conlog = conjugates_arb_log(y, p)

    # This is not totally fool proof.
    # It could still be the case that conlog has higher precision then the
    # elements in B.
    b = ArbMatSpace(base_ring(B), 1, r, false)()

    for i in 1:r
      b[1, i] = conlog[i]
    end

    @vprint :UnitGroup 3 "For $p element b: $b\n"

    v = b*B
    @vprint :UnitGroup 3 "For $p the vector v: $v\n"
  end

  @vprint :UnitGroup 3 "For $p rel: $rel\n"

  @vprint :UnitGroup 2 "Second iteration to check relation ... \n"

  while !_check_relation_mod_torsion(U.units, y, rel, p)
    @vprint :UnitGroup 2 "Precision not high enough, increasing from $p to $(2*p)\n"
    p = 2*p
    @assert p > 0
    p, B = _conj_log_mat_cutoff_inv(U, p)

    @assert p != 0

    conlog = conjugates_arb_log(y, p)

    b = ArbMatSpace(parent(conlog[1]), 1, r, false)()

    for i in 1:r
      b[1, i] = conlog[i]
    end

    @vprint :UnitGroup 3 "For $p element b: $b\n"
    v = b*B
    @vprint :UnitGroup 3 "For $p the vector v: $v\n"
    _find_rational_relation!(rel, v, bound)
    @vprint :UnitGroup 3 "For $p rel: $rel\n"
  end

  if rel_only
    return rel
  end

  if abs(rel[r + 1]) == 1 || rel[r + 1] == 0
    U.rel_add_prec = p
    return false
  end

  m = matrix(FlintZZ, reshape(rel, r + 1, 1))

  h, u = hnf_with_transform(m)

  @assert h[1,1] == 1

  u = inv(u)

  m = sub(u, 1:r+1, 2:r+1)

  U.units =  _transform(vcat(U.units, y), m)

  U.conj_log_mat_cutoff = Dict{Int, arb_mat}()
  U.conj_log_mat_cutoff_inv = Dict{Int, arb_mat}()
  U.tentative_regulator = regulator(U.units, 64)
  U.rel_add_prec = p
  @vprint :UnitGroup 1 "reduction of the new unit group...index improved by $(abs(rel[r+1]))\n"
  @vtime :UnitGroup 1 U.units = reduce(U.units, p)
  return true
end

function _conj_log_mat_cutoff(x::UnitGrpCtx, p::Int)
  #if haskey(x.conj_log_mat_cutoff,  p)
  #  @vprint :UnitGroup 2 "Conj matrix for $p cached\n"
  #  return x.conj_log_mat_cutoff[p]
  #else
    @vprint :UnitGroup 2 "Conj matrix for $p not cached\n"
    x.conj_log_mat_cutoff[p] = _conj_log_mat_cutoff(x.units, p)
    return x.conj_log_mat_cutoff[p]
  #end
end

function _conj_log_mat_cutoff_inv(x::UnitGrpCtx, p::Int)
  @vprint :UnitGroup 2 "Computing inverse of conjugates log matrix (starting with prec $p) ... \n"
  if haskey(x.conj_log_mat_cutoff_inv,  p)
    @vprint :UnitGroup 2 "Inverse matrix cached for $p\n"
    return p, x.conj_log_mat_cutoff_inv[p]
  else
    @vprint :UnitGroup 2 "Inverse matrix not cached for $p\n"
    @vprint :UnitGroup 2 "Trying to invert conj matrix with prec $p \n"
    local y
    try
      y = _conj_log_mat_cutoff(x, p)
    catch e
      #println("I failed to find the _conj_log_mat_cutoff(x, p) for $p")
      println("Error computing the cutoff matrix with precision $p")
      rethrow(e)
    end
    local yy
    try
      yy = inv(y)
      x.conj_log_mat_cutoff_inv[p] = yy
      return p, x.conj_log_mat_cutoff_inv[p]
    catch e
      if !(e isa ErrorException)
        rethrow(e)
      end
      @vprint :UnitGroup 2 "Could not invert with prec $p \n"
      @vprint :UnitGroup 2 "Increasing precision .. (error was $e)\n"
      @v_do :UnitGroup 2 pushindent()
      r = _conj_log_mat_cutoff_inv(x, 2 * p)
      @v_do :UnitGroup 2 popindent()
      return r
    end
  end
end

function _isindependent(u::UnitGrpCtx{T}, y::FacElem{T}) where T
  K = _base_ring(x[1])
  p = u.indep_prec

  deg = degree(K)
  r1, r2 = signature(K)
  rr = r1 + r2
  r = rr - 1 # unit rank

  # This can be made more memory friendly
  while true
    @assert p != 0

    A = _conj_log_mat(u.units, p)

    B = A*transpose(A)
    @vprint :UnitGroup 1 "Computing det of $(rows(B))x$(cols(B)) matrix with precision $(p) ... \n"
    d = det(B)

    y = (Ar(1)//Ar(r))^r * (Ar(21)//Ar(128) * log(Ar(deg))//(Ar(deg)^2))^(2*r)
    if isfinite(d) && ispositive(y - d)
      return false, p
    elseif isfinite(d) && ispositive(d)
      return true, p
    end
    p = 2*p
  end
end

function _rel_add_prec(U)
  return U.rel_add_prec
end

function  _conj_log_mat(x::Array{T, 1}, p::Int) where T
  conlog = conjugates_arb_log(x[1], p)

  r, s = signature(_base_ring(x[1]))
  rr = r + s

  A = zero_matrix(parent(conlog[1]), length(x), rr)

  for i in 1:rr
    A[1, i] = conlog[i]
  end

  Ar = base_ring(A)

  for k in 2:length(x)
    conlog = conjugates_arb_log(x[k], p)
    for i in 1:rr
      A[k, i] = conlog[i]
    end
  end
  return A
end

function _conj_log_mat_cutoff(x::Array{T, 1}, p::Int) where T
  r = length(x)

  conlog = Vector{Vector{arb}}(length(x))
  q = 2
  for i in 1:length(x)
    conlog[i] = conjugates_arb_log(x[i], p)
    for j in 1:length(conlog[i])
      q = max(q, bits(conlog[i][j]))
    end
  end

  A = zero_matrix(ArbField(q, false), r, r)

  for i in 1:r
    for j in 1:r
      A[i, j] = conlog[i][j]
    end
  end

  return A
end

# Powering function for fmpz exponents
function _pow(x::Array{T, 1}, y::Array{fmpz, 1}) where T
  K = _base_ring(x[1])

  zz = deepcopy(y)

  z = Array{fmpz}(length(x))

  for i in 1:length(x)
    z[i] = mod(zz[i], 2)
    zz[i] = zz[i] - z[i]
  end

  r = K(1)

  return zz
end

function _add_unit(u::UnitGrpCtx, x::FacElem{nf_elem, AnticNumberField})
  isindep, p = isindependent(vcat(u.units, [x]), u.indep_prec)
  u.indep_prec = max(p, u.indep_prec)
  if isindep
    push!(u.units, x)
  end
end

ext: jl


text:
export de440

register(DataDep("2003IERSConventions",
                 """
                 2003 IERS Conventions
                 Author: Dennis D. McCarthy and Gérard Petit.
                 Website: https://iers-conventions.obspm.fr/index.php

                 The IERS Conventions describes the reference systems realized by the IERS, in addition to
                 developing and maintaining the models and procedures used to support this endeavor.

                 Citation: IERS Conventions (2003). Dennis D. McCarthy and Gérard Petit.
                         (IERS Technical Note ; 32) Frankfurt am Main: Verlag des Bundesamts für Kartographie und Geodäsie,
                         2004. 127 pp., paperback, ISBN 3-89888-884-3 (print version)
                 """,
                 "https://iers-conventions.obspm.fr/archive/2003/2003IERSConventionsArchive.tar.gz",
                 "1c7ac5c46f49dbfc6e140a34228067b480796ac11e3348ab41905e5af119f708";
                 post_fetch_method=unpack))
# I would like to use a ManualDataDep here, but that doesn't work for unregistered packages for some reason
# Once we register, we will move those constants back here
register(DataDep("nu2000k",
                 """
                 NU2000K Nutation Model
                 Author: George H. Kaplan
                 Website: https://www.usno.navy.mil/USNO/astronomical-applications/software-products/novas/novas-c

                 The 488-term nu2000k nutation model utilizes the luni-solar fundamental argument, planetary longitude, and general
                 precession in longitude expressions from Simon et al. (1994) throughout. The NU2000K Nutation Model is fully documented in the
                 United States Naval Observatory's NOVAS-C 3.0 software package.
                 """, "https://github.com/kiranshila/nu2000k/archive/main.zip",
                 "48b070d61c464da2abcea6a2221831e8ddd78b8d611bb4d9bcc712a50be9eeff";
                 post_fetch_method=function (fn)
                     unpack(fn)
                     dir = "nu2000k-main"
                     innerfiles = readdir(dir)
                     mv.(joinpath.(dir, innerfiles), innerfiles)
                     rm(dir; recursive=true)
                     return rm("main.zip")
                 end))

"""
iers2003(file)

Retrives and gives the location of a `file` from the 2003 IERS Convention Chapter 5 datapack
"""
function iers2003(file::String)
    return resolve("2003IERSConventions/chapter5/$file", @__FILE__)
end

"""
nu2000k(file)

Retrives and gives the location of a `file` from the NU2000k nutation model
"""
function nu2000k(file::String)
    return resolve("nu2000k/$file", @__FILE__)
end

function __init__()
    return register(DataDep("de440",
                            """
                            The JPL Planetary and Lunar Ephemerides DE440 and DE441
                            Author: Ryan S. Park, William M. Folkner, James G. Williams, and Dale H. Boggs
                            Website: https://iopscience.iop.org/article/10.3847/1538-3881/abd414

                            The latest JPL ephemeris with fully consistent treatment of planetary and
                            lunar laser ranging data is DE440 (Park et al., 2021). The dynamical model
                            for DE440 includes a frictional damping between the fluid core and the
                            elastic mantle. This damping term is not suitable for extrapolation more
                            than several centuries into the past. In order to cover a longer time span,
                            the ephemeris DE441 was integrated without the lunar core/mantle damping term.
                            The positions of the planets for DE441 agree with the positions on DE440 to
                            within one meter over the time covered by DE440. For the Moon DE441 differs
                            from DE440 mainly in the estimated tidal damping term causing a difference
                            in along-track position of the Moon of ~10 meters 100 years from the present
                            and growing quadratically for times more than 100 years from present.

                            Citation: Ryan S. Park et al 2021 AJ 161 105
                            """,
                            "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de440.bsp",
                            "a4ce9bf9b3282becc9f4b2ac3cebe03a2ae7599981aabd7265fd8482fff7c4b5"))
end

"""
    de440()

Retrives and gives the location of the JPL DE440 Planetary Ephemeris
"""
function de440()
    return datadep"de440/de440.bsp"
end


text: const con = Contribution


text:
"""
    $(TYPEDEF)

`PointMass` is a simple environment useful for trying out and debugging new algorithms. The
task is simply to move a 2D point mass to a target position by applying x and y forces
to the mass.

# Spaces

* **State: (13, )**
* **Action: (2, )**
* **Observation: (6, )**
"""
struct PointMass{Sim, OSpace} <: AbstractMuJoCoEnvironment
    sim::Sim
    observationspace::OSpace
    function PointMass(sim::MJSim)
        observationspace = MultiShape(
            agent_xy_pos = VectorShape(mjtNum, 2),
            agent_xy_vel = VectorShape(mjtNum, 2),
            target_xy_pos = VectorShape(mjtNum, 2)
        )
        new{typeof(sim), typeof(observationspace)}(sim, observationspace)
    end
end

PointMass() = first(tconstruct(PointMass, 1))

function LyceumBase.tconstruct(::Type{PointMass}, N::Integer)
    modelpath = joinpath(@__DIR__, "pointmass.xml")
    Tuple(PointMass(s) for s in tconstruct(MJSim, N, modelpath, skip=1))
end

@inline getsim(env::PointMass) = env.sim

@inline LyceumBase.observationspace(env::PointMass) = env.observationspace

@propagate_inbounds function LyceumBase.getobservation!(obs, env::PointMass)
    @boundscheck checkaxes(observationspace(env), obs)
    dn = env.sim.dn
    shaped = observationspace(env)(obs)
    @uviews shaped @inbounds begin
        shaped.agent_xy_pos .= dn.xpos[:x, :agent], dn.xpos[:y, :agent]
        shaped.agent_xy_vel .= dn.qvel[:agent_x], dn.qvel[:agent_y]
        shaped.target_xy_pos .= dn.xpos[:x, :target], dn.xpos[:y, :target]
    end
    obs
end

@propagate_inbounds function LyceumBase.getreward(::Any, ::Any, obs, env::PointMass)
    @boundscheck checkaxes(observationspace(env), obs)
    shaped = observationspace(env)(obs)
    @uviews shaped @inbounds begin
        1.0 - euclidean(shaped.agent_xy_pos, shaped.target_xy_pos)
    end
end

@propagate_inbounds function LyceumBase.randreset!(rng::Random.AbstractRNG, env::PointMass)
    reset_nofwd!(env.sim)
    @inbounds begin
        env.sim.dn.qpos[:agent_x] = rand(rng) * 2.0 - 1.0
        env.sim.dn.qpos[:agent_y] = rand(rng) * 2.0 - 1.0
        env.sim.dn.qpos[:target_x] = rand(rng) * 2.0 - 1.0
        env.sim.dn.qpos[:target_y] = rand(rng) * 2.0 - 1.0
    end
    forward!(env.sim)
    env
end


text:
using SparseADRules:get_approximate_basis, low_rank_svd
@testset "low rank svd" begin
    for T in (Float64, ComplexF64)
        approx_rtol = 100*eps(real(T))
        for (n, m) in ((64,32), (32,32))
            for i = 1:5
                A = sprand(n, m, 0.1)
                r = rank(A) # sparse matrix not always low rank
                @testset "Q basis for A" begin
                    Q = get_approximate_basis(A, r)
                    @test norm(A - Q*Q'*A) < approx_rtol*norm(A)
                end

                @testset "svd test for A" begin
                    U, S, Vt = low_rank_svd(A, r)
                    A_approx = U * Diagonal(S) * Vt
                    @test norm(A - A_approx) < approx_rtol*norm(A)
                end

                M = rand(1, m)
                A_m = A .- M
                r_m = rank(A_m)
                @testset "Q basis for A - M" begin
                    Q = get_approximate_basis(A, r_m, 2, M)
                    @test norm(A_m - Q*Q'*A_m) < approx_rtol*norm(A_m)
                end

                @testset "svd test for A - M" begin
                    U, S, Vt = low_rank_svd(A, r_m, 2, M)
                    A_m_approx = U * Diagonal(S) * Vt
                    @test norm(A_m - A_m_approx) < approx_rtol*norm(A_m)
                end
            end
        end
    end
end


text:
function count_markers_temp(r::Robot)
    count = 0
    sum = 0
    side = Ost
    (n_sum, n_count) = line(r, side)
    sum += n_sum
    count += n_count
    while !isborder(r, Nord)
        move!(r, Nord)
        side = inverse(side)
        (n_sum, n_count) = line(r, side)
        sum += n_sum
        count += n_count
    end
    return (sum, count)
end

inverse(side::HorizonSide) = HorizonSide(mod(Int(side) + 2, 4))

# have a stride by the line shore and have all the markers sammeln and return tuple of temp and count
function line(r::Robot, side::HorizonSide)
    count = 0
    sum = 0
    while !isborder(r, side)
        if ismarker(r)
            count += 1
            sum += temperature(r)
        end
        move!(r, side)
    end

    if ismarker(r)
        count += 1
        sum += temperature(r)
    end

    return (sum, count) # tuple of digits
end
# вычисление
function middle(r::Robot)
    (sum, count) = count_markers_temp(r)
    if count == 0
        print(0)
    else
        print(sum / count)
    end
  end

max_stars_repo_name: NikitaDestroyer/Basics

max_issues_repo_name: NikitaDestroyer/Basics

max_forks_repo_name: NikitaDestroyer/Basics


text:
module CUDArt_gen  # the generated module

using Compat

const ext = joinpath(dirname(@__DIR__), "deps", "ext.jl")
if isfile(ext)
    include(ext)
else
    error("Unable to load dependency file $ext.\nPlease run Pkg.build(\"CUDAdrv\") and restart Julia.")
end
const libcudart = libcudart_path

# In the runtime API, these are all used only inside Ptrs,
# so these typealiases are safe (if you don't need access to
# struct elements)
include(joinpath(dirname(@__DIR__), "gen-6.5", "gen_libcudart_h.jl"))
include(joinpath(dirname(@__DIR__), "gen-6.5", "gen_libcudart.jl"))

function checkerror(code::Cuint)
    if code == cudaSuccess
        return nothing
    end
    # Because try/finally may disguise the source of the problem,
    # let's show a backtrace here
    warn("CUDA error triggered from:")
    Base.show_backtrace(STDOUT, backtrace())
    print("\n\n")
    throw(unsafe_string(cudaGetErrorString(code)))
end

# # Make cudaSetupArgument more convenient
# function cudaSetupArgument(arg, size, offset)
#   checkerror(ccall( (:cudaSetupArgument, libcudart), cudaError_t, (Ptr{Void}, Csize_t, Csize_t), &arg, size, offset))
# end

# Fix issues stemming from the inability to pass structs as args
const libwrapcuda = Libdl.find_library("libwrapcuda", [joinpath(dirname(@__DIR__), "deps")])
isempty(libwrapcuda) && error("Cannot find libwrapcuda")

function cudaMalloc3D(p::Array{cudaPitchedPtr,1}, ext::cudaExtent)
    checkerror(ccall((:wrapcudaMalloc3D, libwrapcuda), cudaError_t, (Ptr{cudaPitchedPtr}, Ptr{cudaExtent}), p, &ext))
end
function cudaMemset3D(pitchedDevPtr::cudaPitchedPtr, value::UInt8, extent::cudaExtent)
  checkerror(ccall( (:wrapcudaMemset3D, libwrapcuda), cudaError_t, (Ptr{cudaPitchedPtr}, Cint, Ptr{cudaExtent}), &pitchedDevPtr, value, &extent))
end

# # Issue warning, not error, on bad cudaFree
# function cudaFree(devPtr)
#     code = ccall( (:cudaFree, libcudart), cudaError_t, (Ptr{Void},), devPtr)
#     if code != cudaSuccess
#         warn(string(bytestring(cudaGetErrorString(code)), " on ", devPtr))
#         Base.show_backtrace(STDOUT, backtrace())
#     end
# end

# Convenience constructors
cudaPitchedPtr() = cudaPitchedPtr(C_NULL, 0, 0, 0)
cudaPitchedPtr(pp::cudaPitchedPtr, xsize::Integer) = cudaPitchedPtr(pp.ptr, pp.pitch, xsize, pp.ysize)

end


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

module CallGraph

using CompilerTools
import CompilerTools.DebugMsg
DebugMsg.init()

import Base.show
import CompilerTools.AstWalker
using CompilerTools.LambdaHandling
using CompilerTools.Helper
using CompilerTools.LivenessAnalysis
using Core: IntrinsicFunction
using DataStructures
using DataStructures.Queue
using ..ParallelAccelerator

type CallInfo
    func_sig                             # a tuple of (function, signature)
    could_be_global :: Set               # the set of indices of arguments to the function that are or could be global arrays
    array_parameters :: Dict{Int,Any}    # the set of indices of arguments to the function that are or could be array parameters to the caller
end

type FunctionInfo
    func :: GlobalRef
    ast  :: Expr
    sig  :: Tuple
    calls :: Array{CallInfo,1}
    array_params_set_or_aliased :: Set   # the indices of the parameters that are set or aliased
    can_parallelize :: Bool              # false if a global is written in this function
    recursive_parallelize :: Bool
    LambdaVarInfo :: LambdaVarInfo

    function FunctionInfo(f, a, s)
      new(f, a, s, CallInfo[], Set(), true, true, LambdaVarInfo())
    end
end

type extractStaticCallGraphState
    cur_func_sig
    mapNameFuncInfo :: Dict{Any,FunctionInfo}
    functionsToProcess :: Queue{FunctionInfo}

    calls :: Array{CallInfo,1}
    globalWrites :: Set
    LambdaVarInfo   :: Union{LambdaVarInfo, Void}
    array_params_set_or_aliased :: Set   # the indices of the parameters that are set or aliased
    cant_analyze :: Bool

    local_lambdas :: Dict{Symbol,LambdaInfo}

    function extractStaticCallGraphState(fs, mnfi, ftp)
        new(fs, mnfi, ftp, CallInfo[], Set(), nothing, Set(), false, Dict{Symbol,LambdaInfo}())
    end

    function extractStaticCallGraphState(fs, mnfi, ftp, calls, gw, ll)
        new(fs, mnfi, ftp, calls, gw, nothing, Set(), false, ll)
    end
end

function indexin(a, b::Array{Any,1})
    for i = 1:length(b)
        if b[i] == a
            return i
        end
    end
    throw(string("element not found in input array in indexin"))
end

function resolveFuncByName(cur_func_sig , func :: Symbol, call_sig_arg_tuple, local_lambdas)
    @dprintln(4,"resolveFuncByName ", cur_func_sig, " ", func, " ", call_sig_arg_tuple, " ", local_lambdas)
    ftyp = typeof(cur_func_sig[1])
    if ftyp == Function
        cur_module = Base.function_module(cur_func_sig[1])
    elseif ftyp == LambdaInfo
        cur_module = cur_func_sig[1].module
    else
        throw(string("Unsupported current function type in resolveFuncByName."))
    end
    @dprintln(4,"Starting module = ", Base.module_name(cur_module))

    if haskey(local_lambdas, func)
        return local_lambdas[func]
    end

    while true
        if isdefined(cur_module, func)
            return getfield(cur_module, func)
        else
            if cur_module == Main
                break
            else
                cur_module = Base.module_parent(cur_module)
            end
        end
    end

    if isdefined(Base, func)
        return getfield(Base, func)
    end
    if isdefined(Main, func)
        return getfield(Main, func)
    end
    throw(string("Failed to resolve symbol ", func, " in ", Base.function_name(cur_func_sig[1])))
end

function isGlobal(sym :: Symbol, state)
    return !haskey(state.LambdaVarInfo.var_defs, sym)
end

function getGlobalOrArrayParam(state, real_args)
    possibleGlobalArgs = Set()
    possibleArrayArgs  = Dict{Int,Any}()

    @dprintln(3,"getGlobalOrArrayParam ", real_args)
    for i = 1:length(real_args)
        atyp = typeof(real_args[i])
        @dprintln(3,"getGlobalOrArrayParam arg ", i, " ", real_args[i], " type = ", atyp)
        if atyp<:RHSvar
            aname = lookupVariableName(real_args[i], state.LambdaVarInfo)
            if isGlobal(aname, state)
                # definitely a global so add
                push!(possibleGlobalArgs, i)
            elseif CompilerTools.LambdaHandling.isInputParameter(aname, state.LambdaVarInfo)
                possibleArrayArgs[i] = aname
            end
        else
            # unknown so add
            push!(possibleGlobalArgs, i)
            possibleArrayArgs[i] = nothing
        end
    end

    [possibleGlobalArgs, possibleArrayArgs]
end

function processFuncCall(state, func_expr, call_sig_arg_tuple, possibleGlobals, possibleArrayParams)
    fetyp = typeof(func_expr)

    @dprintln(3,"processFuncCall ", func_expr, " ", call_sig_arg_tuple, " ", fetyp, " ", possibleGlobals)
    if fetyp<:RHSVar
        func = resolveFuncByName(state.cur_func_sig, lookupVariableName(func_expr, state.LambdaVarInfo), call_sig_arg_tuple, state.local_lambdas)
    elseif fetyp == TopNode
        return nothing
    elseif fetyp == DataType
        return nothing
    elseif fetyp == GlobalRef
        func = getfield(func_expr.mod, func_expr.name)
    elseif fetyp == Expr
        @dprintln(3,"head = ", func_expr.head)
        if isCall(func_expr) && isBaseFunc(getCallFunction(func_expr), :getfield)
            args = getCallArguments(func_expr)
            @dprintln(3,"args2 = ", args[1], " type = ", typeof(args[1]))
            @dprintln(3,"args3 = ", args[2], " type = ", typeof(args[2]))
            fsym = args[2]
            if typeof(fsym) == QuoteNode
                fsym = fsym.value
            end
            if typeof(args[1]) == GlobalRef
                func = getfield(eval(args[1]), fsym)
            else
                func = getfield(Main.(args[1]), fsym)
            end
        else
            func = eval(func_expr)
        end
        @dprintln(3,"eval of func_expr is ", func, " type = ", typeof(func))
    else
        throw(string("Unhandled func expression type ", fetyp," for ", func_expr, " in ", Base.function_name(state.cur_func_sig.func)))
    end

    ftyp = typeof(func)
    @dprintln(4,"After name resolution: func = ", func, " type = ", ftyp)
    if ftyp == DataType
        return nothing
    end
    assert(ftyp == Function || ftyp == IntrinsicFunction || ftyp == LambdaInfo)

    if ftyp == Function
        if !isgeneric(func)
            @dprintln(3,func, " is not generic.")
            return nothing
        end
        fgr = GlobalRef(Base.function_module(func, call_sig_arg_tuple), Base.function_name(func))
        fs  = (fgr, call_sig_arg_tuple)

        if !haskey(state.mapNameFuncInfo, fs)
            @dprintln(3,"Adding ", func, " to functionsToProcess")
            method = methods(func, call_sig_arg_tuple)
            if length(method) < 1
                @dprintln(3,"Skipping function not found. ", func, " ", Base.function_name(func), " ", call_sig_arg_tuple)
                return nothing
            else
                @dprintln(3,"method = ", method)
            end

            ct = ParallelAccelerator.Driver.code_typed(func, call_sig_arg_tuple)
            enqueue!(state.functionsToProcess, FunctionInfo(fgr, ct, call_sig_arg_tuple))
            @dprintln(3,state.functionsToProcess)
        end
        push!(state.calls, CallInfo(fs, possibleGlobals, possibleArrayParams))
    elseif ftyp == LambdaInfo
        fgr = GlobalRef(Base.function_module(func, call_sig_arg_tuple), Base.function_name(func))
        fs  = (fgr, call_sig_arg_tuple)

        if !haskey(state.mapNameFuncInfo, fs)
            @dprintln(3,"Adding ", func, " to functionsToProcess")

            method = methods(func, call_sig_arg_tuple)
            if length(method) < 1
                @dprintln(3,"Skipping function not found. ", func, " ", Base.function_name(func), " ", call_sig_arg_tuple)
                return nothing
            else
                @dprintln(3,"method = ", method)
            end

            ct = ParallelAccelerator.Driver.code_typed(func, call_sig_arg_tuple)
            enqueue!(state.functionsToProcess, FunctionInfo(fgr, ct, call_sig_arg_tuple))
            @dprintln(3,state.functionsToProcess)
        end
        push!(state.calls, CallInfo(fs, possibleGlobals, possibleArrayParams))
    end
end

function extractStaticCallGraphWalk(node::Expr,
                                    state::extractStaticCallGraphState,
                                    top_level_number::Int64,
                                    is_top_level::Bool,
                                    read :: Bool)
    @dprintln(4,"escgw: ", node, " type = Expr")

    head = node.head
    args = node.args
    if head == :lambda
        if state.LambdaVarInfo == nothing
            linfo, body = lambdaToLambdaVarInfo(node)
            state.LambdaVarInfo = linfo
            @dprintln(3, "LambdaVarInfo = ", state.LambdaVarInfo)
        else
            new_lambda_state = extractStaticCallGraphState(state.cur_func_sig, state.mapNameFuncInfo, state.functionsToProcess, state.calls, state.globalWrites, state.local_lambdas)
            AstWalker.AstWalk(node, extractStaticCallGraphWalk, new_lambda_state)
            return node
        end
    elseif head == :invoke || head == :call || head == :call1
        func_expr = getCallFunction(node)
        @dprintln(4,"func_expr = ", func_expr)
        if isBaseFunc(func_expr, :ccall)
            state.cant_analyze = true
            return node
        end
        call_args = getCallArguments(node)
        @dprintln(4,"call_args = ", call_args)
        call_sig = Expr(:tuple)
        call_sig.args = map(x -> CompilerTools.LambdaHandling.getType(x, state.LambdaVarInfo), call_args)
        call_sig_arg_tuple = eval(call_sig)
        @dprintln(4,"Arg tuple = ", call_sig_arg_tuple)
        if in(Any, call_sig.args)
            state.cant_analyze = true
            return node
        end

        pmap_name = Symbol("__pmap#39__")
        if isBaseFunc(func_expr, pmap_name)
            func_expr = call_args[3]
            assert(typeof(func_expr) == TypedVar)
            func_expr = toLHSVar(func_expr)
            call_sig = Expr(:tuple)
            push!(call_sig.args, ParallelIR.getArrayElemType(call_args[4]))
            call_sig_arg_tuple = eval(call_sig)
            @dprintln(3,"Found pmap with func = ", func_expr, " and arg = ", call_sig_arg_tuple)
            processFuncCall(state, func_expr, call_sig_arg_tuple, getGlobalOrArrayParam(state, [call_args[4]])...)
            return node
        elseif func_expr == :map
            func_expr = call_args[1]
            call_sig = Expr(:tuple)
            push!(call_sig.args, ParallelIR.getArrayElemType(call_args[2]))
            call_sig_arg_tuple = eval(call_sig)
            @dprintln(3,"Found map with func = ", func_expr, " and arg = ", call_sig_arg_tuple)
            processFuncCall(state, func_expr, call_sig_arg_tuple, getGlobalOrArrayParam(state, [call_args[2]])...)
            return node
        elseif func_expr == :broadcast!
            assert(length(call_args) == 4)
            func_expr = call_args[1]
            @dprintln(3,"ca types = ", typeof(call_args[2]), " ", typeof(call_args[3]), " ", typeof(call_args[4]))
            catype = call_args[2].typ
            assert(call_args[3].typ == catype)
            assert(call_args[4].typ == catype)
            call_sig = Expr(:tuple)
            push!(call_sig.args, ParallelIR.getArrayElemType(catype))
            push!(call_sig.args, ParallelIR.getArrayElemType(catype))
            call_sig_arg_tuple = eval(call_sig)
            @dprintln(3,"Found broadcast! with func = ", func_expr, " and arg = ", call_sig_arg_tuple)
            processFuncCall(state, func_expr, call_sig_arg_tuple, getGlobalOrArrayParam(state, call_args[2:4])...)
            return node
        elseif func_expr == :cartesianarray
            new_lambda_state = extractStaticCallGraphState(state.cur_func_sig, state.mapNameFuncInfo, state.functionsToProcess, state.calls, state.globalWrites, state.local_lambdas)
            AstWalker.AstWalk(call_args[1], extractStaticCallGraphWalk, new_lambda_state)
            for i = 2:length(call_args)
                AstWalker.AstWalk(call_args[i], extractStaticCallGraphWalk, state)
            end
            return node
#        elseif ParallelIR.isArrayset(func_expr) || func_expr == :setindex!
#            @dprintln(3,"arrayset or setindex! found")
#            array_name = call_args[1]
#            if in(toLHSVar(array_name), state.params)
#                @dprintln(3,"can't analyze due to write to array")
#                # Writing to an array parameter which could be a global passed to this function makes parallelization analysis unsafe.
#                #state.cant_analyze = true
#                push!(state.array_params_set_or_aliased, indexin(toLHSVar(array_name), state.params))
#            end
#            return node
        end

        processFuncCall(state, func_expr, call_sig_arg_tuple, getGlobalOrArrayParam(state, call_args)...)
    elseif head == :(=)
        lhs = args[1]
        rhs = args[2]
        rhstyp = typeof(rhs)
        @dprintln(3,"Found assignment: lhs = ", lhs, " rhs = ", rhs, " type = ", typeof(rhs))
        if rhstyp == LambdaInfo
            state.local_lambdas[lhs] = rhs
            return node
        elseif rhstyp<:RHSVar
            rhsname = lookupVariableName(rhs, state.LambdaVarInfo)
            if CompilerTools.LambdaHandling.isInputParameter(rhsname, state.LambdaVarInfo)
                rhsname_typ = state.types[rhsname][2]
                if isArrayType(rhsname_typ)
                    # if an array parameter is the right-hand side of an assignment then we give up on tracking whether some parallelization requirement is violated for this parameter
                    # state.cant_analyze = true
                    push!(state.array_params_set_or_aliased, indexin(rhsname, state.params))
                end
            end
        end
    end

# Don't make any changes to the AST, just record information in the state variable.
return CompilerTools.AstWalker.ASTWALK_RECURSE
end

function extractStaticCallGraphWalk(node::Symbol,
                                    state::extractStaticCallGraphState,
                                    top_level_number::Int64,
                                    is_top_level::Bool,
                                    read::Bool)
    @dprintln(4,"escgw: ", node, " type = Symbol")
    @dprintln(4,"Symbol ", node, " found in extractStaticCallGraphWalk. ", read, " ", isGlobal(node, state))
    if !read && isGlobal(node, state)
        @dprintln(3,"Detected a global write to ", node, " ", state)
        push!(state.globalWrites, node)
    end

    # Don't make any changes to the AST, just record information in the state variable.
    return CompilerTools.AstWalker.ASTWALK_RECURSE
end

function extractStaticCallGraphWalk(node::TypedVar,
                                    state::extractStaticCallGraphState,
                                    top_level_number::Int64,
                                    is_top_level::Bool,
                                    read::Bool)
    @dprintln(4,"escgw: ", node, " type = TypedVar")
    node_sym = lookupVariableName(node, state.LambdaVarInfo)
    @dprintln(4,"TypedVar ", node, " found in extractStaticCallGraphWalk.", read, " ", isGlobal(node_sym, state))
    if !read && isGlobal(node_sym, state)
        @dprintln(3,"Detected a global write to ", node_sym, " ", state.types)
        push!(state.globalWrites, node_sym)
    end
    if getType(node, state.LambdaVarInfo) == Function
        throw(string("Unhandled function in extractStaticCallGraphWalk."))
    end

    # Don't make any changes to the AST, just record information in the state variable.
    return CompilerTools.AstWalker.ASTWALK_RECURSE
end

function extractStaticCallGraphWalk(node::ANY,
                                    state::extractStaticCallGraphState,
                                    top_level_number::Int64,
                                    is_top_level::Bool,
                                    read::Bool)

    # Don't make any changes to the AST, just record information in the state variable.
    return CompilerTools.AstWalker.ASTWALK_RECURSE
end

function show(io::IO, cg :: Dict{Any,FunctionInfo})
    for i in cg
        println(io,"Function: ", i[2].func_sig[1], " ", i[2].func_sig[2], " local_par: ", i[2].can_parallelize, " recursive_par: ", i[2].recursive_parallelize, " array_params: ", i[2].array_params_set_or_aliased)
        for j in i[2].calls
            print(io,"    ", j)
            if haskey(cg, j.func_sig)
                jfi = cg[j.func_sig]
                print(io, " (", jfi.can_parallelize, ",", jfi.recursive_parallelize, ",", jfi.array_params_set_or_aliased, ")")
            end
            println(io,"")
        end
    end
end

function propagate(mapNameFuncInfo)
    # The can_parallelize field is purely local but the recursive_parallelize field includes all the children of this function.
    # Initialize the recursive state to the local state.
    for i in mapNameFuncInfo
        i[2].recursive_parallelize = i[2].can_parallelize
    end

    changes = true
    loop_iter = 1

    # Keep iterating so long as we find changes to a some recursive_parallelize field.
    while changes
        @dprintln(3,"propagate ", loop_iter)
        loop_iter = loop_iter + 1

        changes = false
        # For each function in the call graph
        for i in mapNameFuncInfo  # i of type FunctionInfo
            this_fi = i[2]

            @dprintln(3,"propagate for ", i)
            # If the function so far thinks it can parallelize.
            if this_fi.recursive_parallelize
                # Check everything it calls.
                for j in this_fi.calls # j of type CallInfo
                    @dprintln(3,"    call ", j)
                    if haskey(mapNameFuncInfo, j.func_sig)
                        @dprintln(3,"    call ", j, " found")
                        other_func = mapNameFuncInfo[j.func_sig]
                        # If one of its callees can't be parallelized then it can't either.
                        if !other_func.recursive_parallelize
                            # Record we found a change so the outer iteration must continue.
                            changes = true
                            # Record that the current function can't parallelize.
                            @dprintln(3,"Function ", i[1], " not parallelizable since callee ", j.func_sig," isn't.")
                            this_fi.recursive_parallelize = false
                            break
                        end
                        if !isempty(intersect(j.could_be_global, other_func.array_params_set_or_aliased))
                            # Record we found a change so the outer iteration must continue.
                            changes = true
                            # Record that the current function can't parallelize.
                            @dprintln(3,"Function ", i[1], " not parallelizable since callee ", j.func_sig," intersects at ", intersect(j.could_be_global, other_func.array_params_set_or_aliased))
                            this_fi.recursive_parallelize = false
                            break
                        end
                        for k in j.array_parameters
                            @dprintln(3,"    k = ", k)
                            arg_index = k[1]
                            arg_array_name = k[2]
                            if in(arg_index, other_func.array_params_set_or_aliased)
                                if typeof(arg_array_name) == Symbol
                                    push!(this_fi.array_params_set_or_aliased, indexin(arg_array_name, this_fi.params))
                                else
                                    for l = 1:length(this_fi.params)
                                        ptype = this_fi.types[this_fi.params[l]][2]
                                        if isArrayType(ptype)
                                            push!(this_fi.array_params_set_or_aliased, l)
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end

function extractStaticCallGraph(func :: GlobalRef, ast :: Expr, sig :: Tuple)
    functionsToProcess = Queue(FunctionInfo)
    enqueue!(functionsToProcess, FunctionInfo(func,ast,sig))
    mapNameFuncInfo = Dict{Any,FunctionInfo}()

    while !isempty(functionsToProcess)
        cur_func_sig = dequeue!(functionsToProcess)
        @dprintln(3,"functionsToProcess left = ", functionsToProcess)

        @dprintln(3,"Processing function ", cur_func_sig)

        state = extractStaticCallGraphState(cur_func_sig, mapNameFuncInfo, functionsToProcess)
        AstWalker.AstWalk(cur_func_sig.ast, extractStaticCallGraphWalk, state)
        @dprintln(4,state)
        cur_func_sig.calls = state.calls
        cur_func_sig.array_params_set_or_aliased = state.array_params_set_or_aliased
        cur_func_sig.can_parallelize = !state.cant_analyze && length(state.globalWrites) == 0
        cur_func_sig.recursive_parallelize = true
        cur_func_sig.LambdaVarInfo = state.LambdaVarInfo
        mapNameFuncInfo[(cur_func_sig.func, cur_func_sig.sig)] = cur_func_sig

#        if ftyp == Function
#            @dprintln(3,"Processing function ", the_func, " ", Base.function_name(the_func), " ", the_sig)
#
#            method = methods(the_func, the_sig)
#            if length(method) < 1
#                @dprintln(3,"Skipping function not found. ", the_func, " ", Base.function_name(the_func), " ", the_sig)
#                continue
#            end
#            ct = ParallelAccelerator.Driver.code_typed(the_func, the_sig)
#            @dprintln(4,ct)
#            state = extractStaticCallGraphState(cur_func_sig, mapNameFuncInfo, functionsToProcess)
#            AstWalker.AstWalk(ct, extractStaticCallGraphWalk, state)
#            @dprintln(4,state)
#            mapNameFuncInfo[cur_func_sig] = FunctionInfo(cur_func_sig, state.calls, state.array_params_set_or_aliased, !state.cant_analyze && length(state.globalWrites) == 0, true, state.LambdaVarInfo)
#        elseif ftyp == LambdaInfo
#            @dprintln(3,"Processing lambda static data ", the_func, " ", the_sig)
#            ast = Base.uncompressed_ast(the_func)
#            @dprintln(4,ast)
#            state = extractStaticCallGraphState(cur_func_sig, mapNameFuncInfo, functionsToProcess)
#            AstWalker.AstWalk(ast, extractStaticCallGraphWalk, state)
#            @dprintln(4,state)
#            mapNameFuncInfo[cur_func_sig] = FunctionInfo(cur_func_sig, state.calls, state.array_params_set_or_aliased, !state.cant_analyze && length(state.globalWrites) == 0, true, state.LambdaVarInfo)
#        end
    end

    propagate(mapNameFuncInfo)
    mapNameFuncInfo
end

use_extract_static_call_graph = 0
function use_static_call_graph(x)
    global use_extract_static_call_graph = x
end

end


text:
module sphere_dipole_examples
using FinEtools
using FinEtoolsAcoustics
using LinearAlgebra

function sphere_dipole_transient()
    println("The interior sphere accelerates in the alternately in the positive
    and negative x-direction, generating positive pressure ahead of it, negative
    pressure behind. Time-dependent simulation.
    ")
    rho = 1.21*phun("kg/m^3");# mass density
    c  = 343.0*phun("m/s");# sound speed
    bulk =  c^2*rho;
    a_amplitude=1.0*phun("mm/s^2");# amplitude of the  acceleration of the sphere
    # omega = 2000*phun("rev/s");      # frequency of the incident wave
    R = 50.0*phun("mm"); # radius of the interior sphere
    Ro = 8*R # radius of the external sphere
    P_amplitude = R*rho*a_amplitude; # pressure amplitude
    nref=2;
    nlayers=40;
    tolerance = Ro/(nlayers)/100
    frequency=1200.; # Hz
    omega=2*pi*frequency;
    dt=1.0/frequency/20;
    tfinal=90*dt;
    nsteps=round(tfinal/dt)+1;

    t0  =  time()

    # Hexahedral mesh
    fens,fes  =  H8sphere(R,nref);
    fens1,fes1  =  mirrormesh(fens, fes, [-1.0, 0.0, 0.0], [0.0, 0.0, 0.0],
    renumb =  r(c) = c[[1, 4, 3, 2, 5, 8, 7, 6]]);
    fens,newfes1,fes2 =  mergemeshes(fens1, fes1, fens, fes, tolerance)
    fes = cat(newfes1,fes2)

    # Derive the finite element sets for the boundary
    bfes  =  meshboundary(fes)
    # Outer spherical boundary
    function dout(xyz)
        return xyz/norm(xyz)
    end

    louter = selectelem(fens, bfes, facing = true, direction = dout)
    outer_fes = subset(bfes, louter);

    fens,fes  = H8extrudeQ4(fens, outer_fes, nlayers, (xyz, layer)->(R+layer/nlayers*(Ro-R))*xyz/norm(xyz));
    connected = findunconnnodes(fens, fes); fens, new_numbering = compactnodes(fens, connected);
    fess = renumberconn!(fes, new_numbering);

    geom  =  NodalField(fens.xyz)
    P  =  NodalField(zeros(FCplxFlt,size(fens.xyz,1),1))

    bfes  =  meshboundary(fes)

    # File  =   "Sphere.vtk"
    # vtkexportmesh(File, bfes.conn, geom.values, FinEtools.MeshExportModule.Q4)
    # @async run(`"paraview.exe" $File`)

    # File  =   "Sphere.vtk"
    # vtkexportmesh (File, fes.conn, fens.xyz, MeshExportModule.H8)
    # @async run(`"C:/Program Files (x86)/ParaView 4.2.0/bin/paraview.exe" $File`)

    linner = selectelem(fens, bfes, distance = R, from = [0.0 0.0 0.0], inflate = tolerance)
    louter = selectelem(fens, bfes, facing = true, direction = dout)

    println("Pre-processing time elapsed  =  ",time() - t0,"s")

    t1  =  time()

    numberdofs!(P)

    material = MatAcoustFluid(bulk,rho)
    femm  =  FEMMAcoust(IntegDomain(fes, GaussRule(3, 2)), material)

    S  =  acousticstiffness(femm, geom, P);
    C  =  acousticmass(femm, geom, P);

    abcfemm  =  FEMMAcoustSurf(IntegDomain(subset(bfes, louter), GaussRule(2, 2)), material)
    D  =  acousticABC(abcfemm, geom, P);

    # Inner sphere pressure loading
    function dipole(dpdn, xyz, J, label, t)
        n = cross(J[:,1],J[:,2]);
        n = vec(n/norm(n));
        dpdn[1] = -rho*a_amplitude*sin(omega*t)*n[1]
    end

    dipfemm  =  FEMMAcoustSurf(IntegDomain(subset(bfes, linner), GaussRule(2, 2)), material)

    # Solve
    P0 = deepcopy(P)
    P0.values[:] .= 0.0; # initially all pressure is zero
    vP0 = gathersysvec(P0);
    vP1 = zeros(eltype(vP0), size(vP0));
    vQ0 = zeros(eltype(vP0), size(vP0));
    vQ1 = zeros(eltype(vP0), size(vP0));
    t = 0.0;
    P1 = deepcopy(P0);

    fi  =  ForceIntensity(FCplxFlt, 1, (dpdn, xyz, J, label)->dipole(dpdn, xyz, J, label, t));
    La0 = distribloads(dipfemm, geom, P1, fi, 2);

    A = (2.0/dt)*S + D + (dt/2.)*C;

    step =0;
    while t <= tfinal
      step = step  + 1;
      println("Time $t ($(step)/$(round(tfinal/dt)+1))")
      t = t+dt;
      fi  = ForceIntensity(FCplxFlt, 1, (dpdn, xyz, J, label)->dipole(dpdn, xyz, J, label, t));
      La1 = distribloads(dipfemm, geom, P1, fi, 2);
      vQ1 = A\((2/dt)*(S*vQ0)-D*vQ0-C*(2*vP0+(dt/2)*vQ0)+La0+La1);
      vP1 = vP0 + (dt/2)*(vQ0+vQ1);
      vP0 = vP1;
      vQ0 = vQ1;
      P1 = scattersysvec!(P1, vec(vP1));
      P0 = P1;
      La0 = La1;
    end

    File  =   "sphere_dipole_1.vtk"
    vtkexportmesh(File, fes.conn, geom.values, FinEtools.MeshExportModule.VTK.H8;
     scalars = [( "realP", real(P1.values))])
    @async run(`"paraview.exe" $File`)

    true
end # sphere_dipole_transient

function allrun()
    println("#####################################################")
    println("# sphere_dipole_transient ")
    sphere_dipole_transient()
    return true
end # function allrun

end # module baffled_piston_examples


text:
a = RchA()
rchb(a)
b = Reachable2.B()
ModuleB.f(b)


text:
using .CUDA

asarray(x, arr::CuArray) where T = CuArray(fill(x, ()))
asarray(x::AbstractArray, y::CuArray) = x
asscalar(x::DenseCuArray) = Array(x)[]

Base.Array(x::Base.ReshapedArray{T,0,<:CuArray}) where T = Array(x.parent)

function get_output_array(xs::NTuple{N, DenseCuArray{<:Any,M} where M}, size; has_repeated_indices=true) where N
    CUDA.zeros(promote_type(map(eltype,xs)...), size)
end

CUDA.cudaconvert(A::EinArray{T}) where T = EinArray{T}(cudaconvert.(A.xs), A.x_indexers, A.y_indexer, A.size, A.ICIS, A.OCIS)
CUDA.cu(A::EinArray{T}) where T = EinArray{T}(cu.(A.xs), A.x_indexers, A.y_indexer, A.size, A.ICIS, A.OCIS)

for TP in [:Diag, :Repeat, :Duplicate, :DefaultRule]
    @eval function einsum(::$TP, ix, iy, x::DenseCuArray, size_dict::Dict{LT}) where LT
        loop_einsum(EinCode{((ix...,),),(iy...,)}(), (x,), size_dict)
    end
end

function einsum(::SimpleBinaryRule{('j',), ('j',), ()}, xs::NTuple{2, DenseCuArray})
    dropdims(reshape(xs[1],1,:) * xs[2]; dims=1)
end

function loop_einsum!(code::EinCode{ixs, iy},
                xs::NTuple{N, DenseCuArray{<:Any,M} where M},
                y::DenseCuArray{T,L}, size_dict::Dict{LT}) where {N,L,T, ixs, iy, LT}
    iy_ = _unique(LT,iy)
    NO = length(iy_)
    A = einarray(code, xs, size_dict)
    if NO == length(iy)
        y = reshape(y, fill(1, ndims(A)-NO)...,size(y)...)
        raw = Base.mapreducedim!(x->x, +, y, A)
        if ndims(A)-NO > 0  # fix 1.7 compatibility
            raw = dropdims(raw, dims=(1:ndims(A)-NO...,))
        end
        return raw
    else
        y_ = CUDA.zeros(T, size(A)[end-NO+1:end]...)
        y_ = reshape(y_, fill(1, ndims(A)-NO)...,size(y_)...)
        raw = Base.mapreducedim!(x->x, +, y_, A)
        if ndims(A)-NO > 0  # fix 1.7 compatibility
            raw = dropdims(raw, dims=(1:ndims(A)-NO...,))
        end
        return expanddims!(EinCode{((iy_...,),), iy}(), raw, y)
    end
end

@generated function expandind(::EinCode{ixs,iy}, ind) where {ixs, iy}
    ix = ixs[1]
    ids = map(ii->:(ind[$(findfirst(==(ii), ix))]), iy)
    Expr(:tuple, ids...)
end

function expanddims!(code::EinCode{ixs, iy}, x, y) where {LT,ixs, iy}
    nthreads = 256
    nblocks = cld(prod(size(x)), nthreads)
    CIS = CartesianIndices(x)
    @inline function kernel(code, y, x)
        i = (blockIdx().x-1) * blockDim().x + threadIdx().x
        i > length(x) && return nothing
        @inbounds yi = expandind(code, CIS[i].I)
        @inbounds y[CartesianIndex(yi)] = x[i]
        nothing
    end
    @cuda(blocks=nblocks, threads=nthreads, kernel(code, y, x))
    return y
end

function _batched_gemm(C1::Char, C2::Char, A::DenseCuArray{T1,3}, B::DenseCuArray{T2,3}) where {T1<:CuBlasFloat, T2<:CuBlasFloat}
    CUDA.CUBLAS.gemm_strided_batched(C1, C2, align_eltypes(A,B)...)
end

tensorpermute(A::DenseCuArray, perm) = permutedims(A, perm)
tensorpermute(A::DenseCuArray, perm::Tuple{}) = A

function einsum(::SimpleBinaryRule{(),(), ()}, xs::NTuple{2, DenseCuArray})
    asarray(Array(xs[1])[] * Array(xs[2])[], xs[1])
end

@info("OMEinsum loaded the CUDA module successfully")


text:
using Test, PolynomialBases

if !haskey(ENV, "JULIA_PKGEVAL") # sympy is not install on https://github.com/JuliaComputing/NewPkgEval.jl
  import SymPy
  x, α = SymPy.symbols("x, alpha")

  @test 0 == SymPy.simplify( laguerre(x, 0) - 1 )
  @test 0 == SymPy.simplify( laguerre(x, 1) - ( 1-x ) )
  @test 0 == SymPy.simplify( laguerre(x, 2) - ( x^2 - 4x + 2 ) / 2 )
  @test 0 == SymPy.simplify( laguerre(x, 3) - ( -x^3 + 9x^2 - 18x + 6 ) / 6 )
  @test 0 == SymPy.simplify( laguerre(x, 4) - ( x^4 - 16x^3 + 72x^2 - 96x + 24 ) / 24 )
  @test 0 == SymPy.simplify( laguerre(x, 5) - ( -x^5 + 25x^4 - 200x^3 + 600x^2 - 600x + 120 ) / 120 )
  @test 0 == SymPy.simplify( laguerre(x, 6) - ( x^6 - 36x^5 + 450x^4 - 2400x^3 + 5400x^2 - 4320x + 720 ) / 720 )

  @test 0 == SymPy.simplify( laguerre(x, 0, α) - 1 )
  # https://www.wolframalpha.com/input/?i=LaguerreL%5B1,+alpha,+x%5D
  @test 0 == SymPy.simplify( laguerre(x, 1, α) - ( α-x+1 ) )
  # https://www.wolframalpha.com/input/?i=LaguerreL%5B2,+alpha,+x%5D
  @test 0 == SymPy.simplify( laguerre(x, 2, α) - ( α^2 + 3α + x^2 - 2α*x - 4x + 2 ) / 2 )
  # https://www.wolframalpha.com/input/?i=LaguerreL%5B3,+alpha,+x%5D
  @test 0 == SymPy.simplify( laguerre(x, 3, α) - ( α^3+6α^2+11α-x^3+3α*x^2+9x^2-3α^2*x-15α*x-18x+6 ) / 6 )
  # https://www.wolframalpha.com/input/?i=LaguerreL%5B4,+alpha,+x%5D
  @test 0 == SymPy.simplify( laguerre(x, 4, α) - ( α^4+10α^3+35α^2+50α+x^4-4α*x^3-16x^3+6α^2*x^2+42α*x^2+72x^2-4α^3*x-36α^2*x-104α*x-96x+24 ) / 24 )

  @test_skip laguerre(x, 4)
  @test_skip laguerre(x, 4, α)
  @test_skip laguerre(x, 4, 2)
end
@inferred laguerre(10., 4)
@inferred laguerre(10., 4, 2.)
@inferred laguerre(10., 4, 2)


text:
using GaussianDistributions, Random
import Random.rand
import Distributions.logpdf
import GaussianDistributions: whiten, unwhiten

"""
    StartingPtPrior

Types inheriting from the abstract type `StartingPtPrior` indicate the prior
that is put on the starting point of the observed path following the dynamics
of some stochastic differential equation.
"""
abstract type StartingPtPrior{T} end

"""
    KnownStartingPt{T} <: StartingPtPrior

Indicates that the starting point is known and stores its value in `y`
"""
struct KnownStartingPt{T} <: StartingPtPrior{T}
    y::T
    KnownStartingPt(y::T) where T = new{T}(y)
end

"""
    GsnStartingPt{T,S} <: StartingPtPrior

Indicates that the starting point is equipped with a Gaussian prior with
mean `μ` and covariance matrix `Σ`. It also stores the most recently sampled
white noise `z` used to compute the starting point and a precision matrix
`Λ`:=`Σ`⁻¹. `μ₀` and `Σ₀` are the mean and covariance of the white noise

    GsnStartingPt(μ::T, Σ::S)

Base constructor. It initialises the mean `μ` and covariance `Σ` parameters and
`Λ` is set according to `Λ`:=`Σ`⁻¹
"""
struct GsnStartingPt{T,S} <: StartingPtPrior{T} where {S}
    μ::T
    Σ::S
    Λ::S
    μ₀::T
    Σ₀::UniformScaling

    GsnStartingPt(μ::T, Σ::S) where {T,S} = new{T,S}(μ, Σ, inv(Σ), 0*μ, I)
end

"""
    rand(G::GsnStartingPt, ρ)

Sample new white noise using Crank-Nicolson scheme with memory parameter `ρ` and
a previous value of the white noise stored inside object `G`
"""
function rand(G::GsnStartingPt, z, ρ=0.0)
    zᵒ = rand(Gaussian(G.μ₀, G.Σ₀))
    √(1-ρ)*zᵒ + √(ρ)*z # preconditioned Crank-Nicolson
end

"""
    rand(G::KnownStartingPt, ::Any)

If starting point is known then nothing can be sampled
"""
rand(G::KnownStartingPt, ::Any=nothing, ::Any=nothing) = nothing

"""
    start_pt(G::GsnStartingPt, P)

Compute a new starting point from the white noise for a given posterior
distribution obtained from combining prior `G` and the likelihood encoded by the
object `P`.
"""
function start_pt(z, G::GsnStartingPt, P::GuidPropBridge)
    μ_post = (P.H[1] + G.Λ) \ (P.Hν[1] + G.Λ * G.μ)
    Σ_post = inv(P.H[1] + G.Λ)
    Σ_post = 0.5 * (Σ_post + Σ_post') # remove numerical inaccuracies
    unwhiten(Σ_post, z) + μ_post
end

"""
    start_pt(G::GsnStartingPt, P)

Compute a new starting point from the white noise for a given prior
distribution `G`
"""
start_pt(z, G::GsnStartingPt) = unwhiten(G.Σ, z) + G.μ

"""
    start_pt(G::KnownStartingPt, P)

Return starting point
"""
start_pt(z, G::KnownStartingPt, P::GuidPropBridge) = G.y

"""
    start_pt(G::KnownStartingPt, P)

Return starting point
"""
start_pt(z, G::KnownStartingPt) = G.y

"""
    inv_start_pt(y, G::GsnStartingPt, P::GuidPropBridge)

Compute the driving noise that is needed to obtain starting point `y` under
prior `G` and likelihood `P`
"""
function inv_start_pt(y, G::GsnStartingPt, P::GuidPropBridge)
    μ_post = (P.H[1] + G.Λ) \ (P.Hν[1] + G.Λ * G.μ)
    Σ_post = inv(P.H[1] + G.Λ)
    Σ_post = 0.5 * (Σ_post + Σ_post')
    whiten(Σ_post, y-μ_post)
end

"""
    inv_start_pt(y, G::KnownStartingPt, P::GuidPropBridge)

Starting point known, no need for dealing with white noise
"""
inv_start_pt(y, G::KnownStartingPt, P::GuidPropBridge) = nothing

"""
    logpdf(::KnownStartingPt, y)

Nothing to do so long as `y` is equal to the known starting point inside `G`
"""
function logpdf(G::KnownStartingPt, y)
    (G.y == y) ? 0.0 : error("Starting point is known, but a different value ",
                             "was passed to logpdf.")
end

"""
    logpdf(::GsnStartingPt, y)

log-probability density function evaluated at `y` of a prior distribution `G`
"""
logpdf(G::GsnStartingPt, y) = logpdf(Gaussian(G.μ, G.Σ), y)


text:
module AstroChemistry

export solve_equilibrium_abundances, calc_abund_derived, init_abund, initialize_chemistry_network
export calc_coeff!
export Par, AbundTotal, Network
export fH2selfshield, fCOselfshield

using DelimitedFiles
using Printf
using StaticArrays
using DifferentialEquations
#using NLsolve
using Parameters
using Statistics

#using Sundials #for CVODE_BDF

include("OctTree.jl")
include("GadgetReader.jl")

include("shielding_functions.jl")
include("grain_recomb.jl")
include("constants.jl")
include("network.jl")
include("system_rate_eqs.jl")

end #module


text:
## Julia program for Generation of Gamma Random Variable
## author: weiya <szcfweiya@gmail.com>
## date: 2018-08-21

## sample from Ga(a, beta)
function rgamma_int(a::Int, beta)
    u = rand(a)
    return(-1.0 / beta * sum(log.(u)))
end

## density of Ga(alpha, beta)
include("func_gamma.jl")
function dgamma(x, alpha, beta)
    return(beta^alpha / lanczos_gamma(alpha) * x^(alpha-1) * exp(-1*beta*x))
end

## accept-reject algorithm
function rgamma(alpha = 1.5, beta = 2.1)
    a = Int(floor(alpha))
    b = a * beta / alpha
    #M = exp(a * (log(a) - 1) - alpha * (log(alpha) - 1))
    M = exp(alpha * (log(alpha) - 1) - a * (log(a) - 1))
    while true
        x = rgamma_int(a, b)
        u = rand()
        cutpoint = dgamma(x, alpha, beta)/(M*dgamma(x, a, b))
        if u <= cutpoint
            return x
        end
    end
end

# example
println(rgamma())

max_stars_repo_name: szcf-weiya/MCMC

max_issues_repo_name: szcf-weiya/MCMC

max_forks_repo_name: szcf-weiya/MCMC