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

using Markdown
using InteractiveUtils

# ╔═╡ f30a37ec-6d8c-11ec-1c66-89962478f84e
begin

	# load required external packages -
	using PlutoUI
	using TOML
	using Dates
	using PQPolygonSDK

	# show -
	nothing
end

# ╔═╡ ee95e493-b7a6-43df-88c0-88ecf5028843
md"""
# PQPolygonSDK.jl
"""

# ╔═╡ b0aa81e5-adcd-401e-bce8-ee63960378ce
md"""
### Introduction

Data is the key to making great decisions, especially in finance. [Polygon.io](https://polygon.io) is a financial data warehouse with an application programming interface (API) that you can use to access all kinds of financial data and other information. The `PQPolygonSDK.jl` package is a software development kit (SDK) written in the [Julia programing language](https://julialang.org) for the [Polygon.io](https://polygon.io) application programming interface.

[Polygon.io](https://polygon.io) was founded by former Googlers in 2016. Since September of 2020, [Polygon.io](https://polygon.io) has provided both [free and paid access](https://polygon.io/pricing) to current and historical [stock](https://polygon.io/docs/stocks/getting-started), [option](https://polygon.io/docs/options/getting-started), [crypto](https://polygon.io/docs/crypto/getting-started) and [forex](https://polygon.io/docs/forex/getting-started) data.
Check out the [Polygon.io blog](https://polygon.io/blog/) for the latest updates and developments from [Polygon.io](https://polygon.io).

This notebook will discuss how to use the `PQPolygonSDK.jl` package to query the [Polygon.io REST API](https://polygon.io/docs/stocks). In particular, we'll:

* Download and analyze daily stock price data for an arbitrary ticker symbol and a specified date range
* Download and analyze crypto data for an arbitrary ticker and a specified date range
* Download information about a stock ticker, e.g., Apple (AAPL) and Microsoft (MSFT)
"""

# ╔═╡ fdc9fa10-4ac0-4055-ae8a-db3f16f2ae36

# ╔═╡ 66c5e100-ce9c-4387-a662-cd93cf76718b

# ╔═╡ 672868d2-c09a-4230-8d90-a09d7af7ae8d

# ╔═╡ dbd1eead-5e9f-4861-90a5-badcb5b2ad57

# ╔═╡ 2c357a88-e6da-420f-bb50-0a7bcd824de3
md"""
##### Installation and Requirements

`PQPolygonSDK.jl` can be installed, updated, or removed using the [Julia package management system](https://docs.julialang.org/en/v1/stdlib/Pkg/). To access the package management interface, open the [Julia REPL](https://docs.julialang.org/en/v1/stdlib/REPL/), and start the package mode by pressing `]`.
While in package mode, to install `PQPolygonSDK.jl`, issue the command:

```julia
(@v1.7.x) pkg> add PQPolygonSDK
```

Once installed, to use `PQPolygonSDK.jl` in your projects (like a Pluto notebook), issue the command:

```julia
julia> using PQPolygonSDK
```
"""

# ╔═╡ f87b1e0a-e6bb-4931-a590-15c2f2865572

# ╔═╡ 106d8bd6-beaf-4552-8a1c-e0536137ad43

# ╔═╡ 7fcd9ada-dd7c-4a02-9f12-b998ea9df3ad

# ╔═╡ 9c38031a-7a57-45ce-baa7-c13df489a3d8

# ╔═╡ 44891bab-deb2-4c55-8257-9d9ede10761c
md"""
### Overview of a Polygon API call using PQPolygonSDK

There are four steps required to make a Polygon.io API call using the `PQPolygonSDK.jl` package:

1. Create a `PQPolygonSDKUserModel` object
1. Create an instance of an `AbstractPolygonEndpointModel` object
1. Use the user and endpoint models to build the actual URL string
1. Call the API and get data back.

##### Build a user model
All [Polygon.io](https://polygon.io) application programming interface (API) calls in the `PQPolygonSDK.jl` package are started by creating a `PQPolygonSDKUserModel` object using the function:

```julia
model(userModelType::Type{PQPolygonSDKUserModel},
    options::Dict{String,Any}) -> PQPolygonSDKUserModel
```

where `options` dictionary holds `email` and `apikey` key-value pairs.

##### Build an endpoint model
The `PQPolygonSDKUserModel` object can then be passed as an argument to a `model` function:

```julia
model(apiModelType::Type{T}, userModel::PQPolygonSDKUserModel,
        options::Dict{String,Any}) -> AbstractPolygonEndpointModel where T<:AbstractPolygonEndpointModel
```

to build an endpoint-specific model where the information for a particular endpoint is passed into the `model` method via an `options` dictionary.

##### Encode the URL string
The endpoint specific model is then used to create the specific URL for the call using the `url` function:

```julia
# create the URL string for this endpoint -
my_stock_url_string = url(POLYGON_URL_STRING, stock_endpoint_model);
```

##### Make the API call
Finally, the URL and type of API model type are used to make the API call:

```julia
# make the API call -
(h_stock, df_stock) = api(PolygonAggregatesEndpointModel, my_stock_url_string);
```

The `api` method returns two pieces of data, a header dictionary holding technical information about the call (and potentially error information of the call failed). The data for the call is returned as a [DataFrame](https://dataframes.juliadata.org/stable/).

"""

# ╔═╡ 40db3a68-cd4e-45df-89a0-f9f1eb8d7227

# ╔═╡ 96fde2cf-c457-4f40-a8d3-1351a02abf49

# ╔═╡ f5398c1d-4f6f-4d10-8410-0ed9db82ae74

# ╔═╡ f650a3ae-56a3-4807-8e99-2a7d1693bff0

# ╔═╡ fc1cad78-b63e-445e-9cb1-5c4eb9f2fd89

# ╔═╡ 3c09ecd0-2a71-4f29-b8fc-69ab20bfaf93

# ╔═╡ 2a58e169-67fa-4222-b33b-17dc18eee230

# ╔═╡ bc73c133-eb6a-462c-90e9-497dfed184f9

# ╔═╡ 8ff5ac28-b6e6-42e4-b99d-2f26cda1e9be
md"""
### Aggregates Endpoint

[The Aggregates endpoint on Polygon.io](https://polygon.io/docs/stocks/get_v2_aggs_ticker__stocksticker__range__multiplier___timespan___from___to) can be used to access the Open/High/Low/Close (OHLC) price (along with other data like volume and volume-weighted average price or vwap) at different time resolutions ranging from minutes to yearly for stocks, options, cryptocurrencies and foreign exchange (forex). In this section, let's look at examples for stocks, cryptocurrencies, and forex. We'll save options for a later time.

A `PolygonAggregatesEndpointModel` is constructed using the `model` function:

```julia
model(apiModelType::Type{T}, userModel::PQPolygonSDKUserModel,
        options::Dict{String,Any}) -> AbstractPolygonEndpointModel
			where T<:AbstractPolygonEndpointModel
```

where the `options` dictionary holds the data for the API call as `key => value` pairs, see the [Aggregates endpoint documentation](https://polygon.io/docs/stocks/get_v2_aggs_ticker__stocksticker__range__multiplier___timespan___from___to) for the data that must be included.  Lastly, we use the `PolygonAggregatesEndpointModel` for stocks, crypto, forex, and options, where the ticker symbol convention for each type of asset is the only difference.

"""

# ╔═╡ 041c79b9-d930-4706-a020-8fc2183ef78c
md"""
##### User model
"""

# ╔═╡ fe7259dc-6d68-4561-b111-7b8eb62bef25
md"""
##### Download Equity/ETF data

Ticker symbols are unique alphabetical codes that represent a company, e.g., `AAPL` (Apple) or an
[exchange traded fund (ETF)](https://www.investopedia.com/articles/investing/122215/spy-spdr-sp-500-trust-etf.asp) such as `SPY`,  which tracks the [Standard & Poor's 500 Index](https://www.investopedia.com/terms/s/sp500.asp).

Ticker symbols for stocks, e.g., `AAPL` or `MSFT`, are relatively easy to find and understand. However, ticker symbols for crypto, forex, and options, are more complicated. For now, let's use the well-known ticker symbol `SPY` to demonstrate [the Aggregates endpoint on Polygon.io](https://polygon.io/docs/stocks/get_v2_aggs_ticker__stocksticker__range__multiplier___timespan___from___to).

However, [Polygon.io](https://polygon.io) has a collection of [reference data endpints](https://polygon.io/docs/stocks/get_v3_reference_tickers) for querying information about ticker symbols that we'll demonstrate below.
"""

# ╔═╡ a42d8dc9-489f-408a-b546-5437b425bfc4
#h_stock

# ╔═╡ 2dd6f756-9b3f-4d35-a6d8-d82b0c878b5d
# df_stock

# ╔═╡ dbc25db6-e592-4e7b-82aa-b7ea9dd7f1ee
md"""
##### Download Crypto data

For cryptocurrencies, the [Polygon.io](https://polygon.io) convention for ticker symbols is `X:<cyrpto>XXX`, where `crypto` denotes the particular cryptocurrency and `XXX` denotes the currency you want to see the price of the crypto measured in.

In the United States, we
would typically use United States Dollars `USD`, but in the Eurozone, the measurement currency would be in € or `EUR`. For example, for [bitcoin](https://bitcoin.org/) demarcated United States Dollars (USD), the [Polygon.io](https://polygon.io) ticker is given by: `X:BTCUSD`, while the [bitcoin](https://bitcoin.org/) ticker symbol demarcated in € would be: `X:BTCEUR`.

"""

# ╔═╡ 6617f059-c2ad-461e-8e83-73c85441ab5d
# h_crypto

# ╔═╡ 9368f7d5-48d2-445c-a55b-2fcc2bf21f8d
# df_crypto

# ╔═╡ eb846a5e-d78c-4208-a828-6591a2b96007

# ╔═╡ 236f21a7-aaa8-46ef-ae11-fb8c86d6872c

# ╔═╡ 374ab839-0ce6-47f4-97ca-cb345fd4493c

# ╔═╡ 590b11f9-97b3-4e6a-99eb-3cee1542a204

# ╔═╡ d1980d23-9c08-46c2-a49a-d761ff0e8866

# ╔═╡ 52bfdb59-4b1a-4314-be86-b1b1a4ba60b2

# ╔═╡ 1d13c045-2812-4bcd-8fda-5ac28f807497

# ╔═╡ 263dfd23-fa6a-41a5-a701-1910e7623453

# ╔═╡ 63b64f5b-793b-4ba1-9df0-ebab53fa1d31

# ╔═╡ 617ce76c-d5f2-4fd9-89c3-eef216285f81

# ╔═╡ 71316dc9-bf47-4e03-9330-b94139e765c1

# ╔═╡ 076b8e58-e54d-4bca-81db-bdb8bb34c41a

# ╔═╡ 96166eaa-561b-4ef5-b2d4-f649f0f2dae2
md"""
### Ticker Details Endpoint
The `PolygonTickerDetailsEndpointModel` model can be used to download details for a ticker symbol's company. This endpoint provides a general overview of the entity with information such as name, sector, exchange, logo, and similar companies.
"""

# ╔═╡ e86aa566-7c17-46c9-bb22-b135f6e5892f
md"""
##### Ticker details for AAPL
"""

# ╔═╡ bd00cfb0-8dae-4ae9-a8f4-2038ed1f4882
# h_ticker_data

# ╔═╡ 2de55011-527d-4a4a-836c-26e8a0803015
# df_ticker_data

# ╔═╡ 7a8df8f8-a7e1-43e3-a054-57ead5299835
# download logo -
# PlutoUI.Resource(df_ticker_data[!,:logo][1])

# ╔═╡ 08a5bf4e-e014-47ae-a947-cdd0abf44597

# ╔═╡ 14d185e3-162d-4570-98c6-454779005f4b

# ╔═╡ 282be644-c196-4fce-83fa-8483141e8671

# ╔═╡ 37e4cbfd-e4a4-4455-9756-7b086088ff28

# ╔═╡ 15f2fddb-01ab-45ce-a4d6-7c134fd6022e

# ╔═╡ eb14de30-9c5c-4887-9e72-e552a97eb232

# ╔═╡ 5444a0ea-1968-44de-a94e-50fca7f96bb9

# ╔═╡ 0c78d113-2b63-4bee-93fe-4061a0834898

# ╔═╡ 1e12c204-be9d-4222-b2b9-c38a6bd615cb

# ╔═╡ 52fd0a3b-d676-4b7d-bafe-bc4bdff4f3b1

# ╔═╡ 5edcd3be-6993-4085-8219-9b7d02ddc4d8

# ╔═╡ 09c54c83-57de-4df0-9212-823cdd6507e0
md"""
### Summary and Conclusions
In this notebook, we covered some of the `PQPolygonSDK.jl` package features. `PQPolygonSDK.jl` is a software development kit (SDK) written in the [Julia programing language](https://julialang.org) for the [Polygon.io](https://polygon.io) application programming interface.

* Design paradigm used to model [Polygon.io](https://polygon.io) api calls: `user model` $\rightarrow$ `endpoint model` $\rightarrow$ `api call`. The `PQPolygonSDK.jl` package
* Additional `endpoint models` are under development
* Data is returned from the methods as a `DataFrame`
* Metadata about an `api call` (including error messages) is returned as a `Dict`
"""

# ╔═╡ 2592fa9a-4194-4206-ac8b-e9727d63b550
TableOfContents(title="📚 Table of Contents", indent=true, depth=5, aside=true)

# ╔═╡ 95e08b94-ba07-47e5-8982-8719d1af8877
html"""
<style>
main {
    max-width: 1200px;
    width: 72%;
    margin: auto;
    font-family: "Roboto, monospace";
}

a {
    color: blue;
    text-decoration: none;
}

.H1 {
    padding: 0px 30px;
}
</style>"""

# ╔═╡ 820238fb-a140-4f38-b58b-745df9b6cf5b
html"""
<script>
	// initialize -
	var section = 0;
	var subsection = 0;
	var headers = document.querySelectorAll('h3, h5');

	// main loop -
	for (var i=0; i < headers.length; i++) {

		var header = headers[i];
	    var text = header.innerText;
	    var original = header.getAttribute("text-original");
	    if (original === null) {

			// Save original header text
	        header.setAttribute("text-original", text);
	    } else {

			// Replace with original text before adding section number
	        text = header.getAttribute("text-original");
	    }

	    var numbering = "";
	    switch (header.tagName) {
	        case 'H3':
	            section += 1;
	            numbering = section + ".";
	            subsection = 0;
	            break;
	        case 'H5':
	            subsection += 1;
	            numbering = section + "." + subsection;
	            break;
	    }
		// update the header text
		header.innerText = numbering + " " + text;
	};
</script>"""

# ╔═╡ 5d57883d-c34a-4e5f-9106-b6cad4050a66
begin

	# setup paths -
	_PATH_TO_ROOT = pwd()
	_PATH_TO_CONFIG = joinpath(_PATH_TO_ROOT,"config")

	# load the configuration information -
	# alias the Polygon.io URL -
    POLYGON_URL_STRING = "https://api.polygon.io"

    # What we have here is the classic good news, bad news situation ...
    # Bad news: We don't check in our Polygon.io API key to GitHub (sorry).
    # Good news: Polygon.io API keys are free. check out: https://polygon.io -
    configuration_dictionary = TOML.parsefile(joinpath(_PATH_TO_CONFIG, "Configuration.toml"))

    # get some stuff from the configuration dictionary -
	# this needs to be YOUR api key, either hardcoded -or- from a config file
    POLYGON_API_KEY = configuration_dictionary["API"]["polygon_api_key"]
	POLYGON_API_EMAIL = configuration_dictionary["EMAIL"]["email"]

	# show -
	nothing
end

# ╔═╡ 08a4644d-46b8-48f3-a2ea-30a6474aa966
begin

	# create an options dictionary - holds data about the user (you)
	user_model_options_dict = Dict{String,Any}()
	user_model_options_dict["apikey"] = POLYGON_API_KEY 	# your API key goes here
	user_model_options_dict["email"] = POLYGON_API_EMAIL 	# your API email goes here

	# create my user model -
	my_user_model = model(PQPolygonSDKUserModel, user_model_options_dict);

	# show -
	nothing
end

# ╔═╡ e4d1c872-5b50-4f7d-aa2e-044c8aef2040
begin

	# what ticker do we want to download?
	stock_ticker_symbol = "SPY";

	# setup the parameters for the call to the endpoint -
	stock_endpoint_options = Dict{String,Any}()
	stock_endpoint_options["adjusted"] = true
	stock_endpoint_options["sortdirection"] = "asc"
	stock_endpoint_options["limit"] = 5000
	stock_endpoint_options["to"] = Date(2021, 12, 24)
	stock_endpoint_options["from"] = Date(2021, 12, 15)
	stock_endpoint_options["multiplier"] = 1
	stock_endpoint_options["timespan"] = "day"
	stock_endpoint_options["ticker"] = stock_ticker_symbol;
	stock_endpoint_model = model(PolygonAggregatesEndpointModel,
		my_user_model, stock_endpoint_options)

	# create the URL string for this endpoint -
	my_stock_url_string = url(POLYGON_URL_STRING, stock_endpoint_model);

	# make the API call -
	(h_stock, df_stock) = api(PolygonAggregatesEndpointModel, my_stock_url_string);

	# show -
	nothing;
end

# ╔═╡ 516589ec-3f40-4870-9a37-8e030cf1c951
begin

	# what crypto ticker do we want to download?
	crypto_ticker_symbol = "X:BTCUSD";

	# setup the parameters for the call to the endpoint -
	crypto_endpoint_options = Dict{String,Any}()
	crypto_endpoint_options["adjusted"] = true
	crypto_endpoint_options["sortdirection"] = "asc"
	crypto_endpoint_options["limit"] = 5000
	crypto_endpoint_options["to"] = Date(2021, 12, 24)
	crypto_endpoint_options["from"] = Date(2021, 12, 15)
	crypto_endpoint_options["multiplier"] = 1
	crypto_endpoint_options["timespan"] = "day"
	crypto_endpoint_options["ticker"] = crypto_ticker_symbol;
	crypto_endpoint_model = model(PolygonAggregatesEndpointModel,
		my_user_model, crypto_endpoint_options);

	# create the URL string -
	my_crypto_url_string = url(POLYGON_URL_STRING, crypto_endpoint_model);

	# make the API call -
	(h_crypto, df_crypto) = api(PolygonAggregatesEndpointModel, my_crypto_url_string);

	# show -
	nothing;
end

# ╔═╡ e0a7bb44-20f7-40b1-9597-a6f0ec1162f0
begin

	# now that we have the user_model, let's build an endpoint model -
	ticker_data_endpoint_options = Dict{String,Any}()
	ticker_data_endpoint_options["ticker"] = "AAPL"
	ticker_data_endpoint_model = model(PolygonTickerDetailsEndpointModel, my_user_model,
		ticker_data_endpoint_options);

	# build the ticker data url -
	my_ticker_data_url_string = url(POLYGON_URL_STRING, ticker_data_endpoint_model)

	# make the call -
	(h_ticker_data, df_ticker_data) = api(PolygonTickerDetailsEndpointModel,
		my_ticker_data_url_string)

	# show -
	nothing
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
PQPolygonSDK = "a9381516-e38f-4c81-935c-32707fb4df4c"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[compat]
PQPolygonSDK = "~0.1.3"
PlutoUI = "~0.7.30"
"""

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

julia_version = "1.7.1"
manifest_format = "2.0"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[deps.PQPolygonSDK]]
deps = ["CSV", "DataFrames", "Dates", "HTTP", "JSON", "TOML"]
git-tree-sha1 = "f5dfecd13178b9c9a00304f64591a375d1f3ce22"
uuid = "a9381516-e38f-4c81-935c-32707fb4df4c"
version = "0.1.3"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

# ╔═╡ Cell order:
# ╟─ee95e493-b7a6-43df-88c0-88ecf5028843
# ╟─b0aa81e5-adcd-401e-bce8-ee63960378ce
# ╟─fdc9fa10-4ac0-4055-ae8a-db3f16f2ae36
# ╟─66c5e100-ce9c-4387-a662-cd93cf76718b
# ╟─672868d2-c09a-4230-8d90-a09d7af7ae8d
# ╟─dbd1eead-5e9f-4861-90a5-badcb5b2ad57
# ╟─2c357a88-e6da-420f-bb50-0a7bcd824de3
# ╠═f30a37ec-6d8c-11ec-1c66-89962478f84e
# ╟─f87b1e0a-e6bb-4931-a590-15c2f2865572
# ╟─106d8bd6-beaf-4552-8a1c-e0536137ad43
# ╟─7fcd9ada-dd7c-4a02-9f12-b998ea9df3ad
# ╟─9c38031a-7a57-45ce-baa7-c13df489a3d8
# ╟─44891bab-deb2-4c55-8257-9d9ede10761c
# ╟─40db3a68-cd4e-45df-89a0-f9f1eb8d7227
# ╟─96fde2cf-c457-4f40-a8d3-1351a02abf49
# ╟─f5398c1d-4f6f-4d10-8410-0ed9db82ae74
# ╟─f650a3ae-56a3-4807-8e99-2a7d1693bff0
# ╟─fc1cad78-b63e-445e-9cb1-5c4eb9f2fd89
# ╟─3c09ecd0-2a71-4f29-b8fc-69ab20bfaf93
# ╟─2a58e169-67fa-4222-b33b-17dc18eee230
# ╟─bc73c133-eb6a-462c-90e9-497dfed184f9
# ╟─8ff5ac28-b6e6-42e4-b99d-2f26cda1e9be
# ╟─041c79b9-d930-4706-a020-8fc2183ef78c
# ╠═08a4644d-46b8-48f3-a2ea-30a6474aa966
# ╟─fe7259dc-6d68-4561-b111-7b8eb62bef25
# ╠═e4d1c872-5b50-4f7d-aa2e-044c8aef2040
# ╠═a42d8dc9-489f-408a-b546-5437b425bfc4
# ╠═2dd6f756-9b3f-4d35-a6d8-d82b0c878b5d
# ╟─dbc25db6-e592-4e7b-82aa-b7ea9dd7f1ee
# ╠═516589ec-3f40-4870-9a37-8e030cf1c951
# ╠═6617f059-c2ad-461e-8e83-73c85441ab5d
# ╠═9368f7d5-48d2-445c-a55b-2fcc2bf21f8d
# ╟─eb846a5e-d78c-4208-a828-6591a2b96007
# ╟─236f21a7-aaa8-46ef-ae11-fb8c86d6872c
# ╟─374ab839-0ce6-47f4-97ca-cb345fd4493c
# ╟─590b11f9-97b3-4e6a-99eb-3cee1542a204
# ╟─d1980d23-9c08-46c2-a49a-d761ff0e8866
# ╟─52bfdb59-4b1a-4314-be86-b1b1a4ba60b2
# ╟─1d13c045-2812-4bcd-8fda-5ac28f807497
# ╟─263dfd23-fa6a-41a5-a701-1910e7623453
# ╟─63b64f5b-793b-4ba1-9df0-ebab53fa1d31
# ╟─617ce76c-d5f2-4fd9-89c3-eef216285f81
# ╟─71316dc9-bf47-4e03-9330-b94139e765c1
# ╟─076b8e58-e54d-4bca-81db-bdb8bb34c41a
# ╟─96166eaa-561b-4ef5-b2d4-f649f0f2dae2
# ╟─e86aa566-7c17-46c9-bb22-b135f6e5892f
# ╠═e0a7bb44-20f7-40b1-9597-a6f0ec1162f0
# ╠═bd00cfb0-8dae-4ae9-a8f4-2038ed1f4882
# ╠═2de55011-527d-4a4a-836c-26e8a0803015
# ╠═7a8df8f8-a7e1-43e3-a054-57ead5299835
# ╟─08a5bf4e-e014-47ae-a947-cdd0abf44597
# ╟─14d185e3-162d-4570-98c6-454779005f4b
# ╟─282be644-c196-4fce-83fa-8483141e8671
# ╟─37e4cbfd-e4a4-4455-9756-7b086088ff28
# ╟─15f2fddb-01ab-45ce-a4d6-7c134fd6022e
# ╟─eb14de30-9c5c-4887-9e72-e552a97eb232
# ╟─5444a0ea-1968-44de-a94e-50fca7f96bb9
# ╟─0c78d113-2b63-4bee-93fe-4061a0834898
# ╟─1e12c204-be9d-4222-b2b9-c38a6bd615cb
# ╟─52fd0a3b-d676-4b7d-bafe-bc4bdff4f3b1
# ╟─5edcd3be-6993-4085-8219-9b7d02ddc4d8
# ╟─09c54c83-57de-4df0-9212-823cdd6507e0
# ╟─2592fa9a-4194-4206-ac8b-e9727d63b550
# ╟─95e08b94-ba07-47e5-8982-8719d1af8877
# ╟─820238fb-a140-4f38-b58b-745df9b6cf5b
# ╟─5d57883d-c34a-4e5f-9106-b6cad4050a66
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: Paliquant/PQPackageDocumentation

max_issues_repo_name: Paliquant/PQPackageDocumentation

max_forks_repo_name: Paliquant/PQPackageDocumentation


text:
#speed of light [m/s]
const 𝐜 = 299792458.0
#Planck constant [J*s]
const 𝐡 = 6.62607015e-34
#Boltzmann constant [J/Kelvin]
const 𝐤 = 1.38064852e-23
#Stefan-Boltzmann constant [W/m^2*Kelvin^4]
const 𝛔 = 5.67037442e-8
#gas constant [J/K*mole], equivalent to kB*Av
const 𝐑 = 8.31446262
#Pascals in 1 atm
const 𝐀 = 101325.0
#avogadros number [molecules/mole]
const 𝐍𝐚 = 6.02214076e23
#Dalton [kg] (mass of particle basically, also ≡ 1/𝐍𝐚/1000)
const 𝐃𝐚 = 1.66053907e-27
#gravitational constant [m^3/kg/s^2]
const 𝐆 = 6.6743e-11
#the Loschmidt number, but in molecules/cm^3 then squared [molecules^2/cm^6]
const 𝐋𝐨² = 7.21879268e38

#reference temperature of HITRAN database [K]
const 𝐓ᵣ = 296.0
#reference temperature [K] equivalent to 0 degrees Celsius
const 𝐓₀ = 273.15
#minimum pressure for certain temperature/pressure profiles
const 𝐏ₘᵢₙ = 1e-9

ext: jl


text:
function pmc(Cs::Vector{Config{D,S}}, M::Model, Obs::Vector{Thermo}, Grad::Vector{Gradient}, N::Int, Sto::Vector{StochasticMatrix}, num_chains::Int) where {D,S}

    Es = pmap(i -> mc(Cs[i], M, Obs[i], Grad[i], N, Sto[i]), 1:num_chains)

    E = mean([Es[i][1] for i = 1:size(Es,1)])
    Eσ = std([Es[i][1] for i = 1:size(Es,1)])

    mag = mean([Es[i][3] for i = 1:size(Es,1)])
    Mσ = std([Es[i][3] for i = 1:size(Es,1)])

    grad = mean([Es[i][2][:total] for i = 1:size(Es,1)])
    # sgrad = mean([Es[i][4] for i = 1:size(Es,1)])

    return E, Eσ, mag, Mσ, grad
end


text:
function naive_index_from_time(sample_rate, sample_time)
    # This stepping computation is prone to roundoff error, so we'll work in high precision
    sample_time_in_seconds = big(Dates.value(Nanosecond(sample_time))) //
                             big(TimeSpans.NS_IN_SEC)
    # At time 0, we are at index 1
    t = Rational{BigInt}(0 // 1)
    index = 1
    while true
        # Now step forward in time; one index, and time 1/sample_rate
        t += 1 // sample_rate
        index += 1
        if t > sample_time_in_seconds
            # we just passed it, so previous index is the last one before the time of interest
            return index - 1
        end
    end
end

# Modified from
# https://github.com/beacon-biosignals/TimeSpans.jl/blob/e3c999021336e51a08d118e6defb792e38ac1cc7/test/runtests.jl#L116-L126
@testset "index_and_error_from_time" begin
    for rate in (101 // 2, 1001 // 10, 200, 256, 1, 10)
        for sample_time in
            (Nanosecond(12345), Minute(5), Nanosecond(Minute(5)) + Nanosecond(1),
             Nanosecond(1), Nanosecond(10^6), Nanosecond(6970297031))
            # compute with a very simple algorithm
            index = naive_index_from_time(rate, sample_time)
            # Check against our `TimeSpans.index_from_time`:
            @test index ==
                  AlignedSpans.index_and_error_from_time(rate, sample_time, RoundDown)[1]
            # Works even if `rate` is in Float64 precision:
            @test index ==
                  AlignedSpans.index_and_error_from_time(Float64(rate), sample_time,
                                                         RoundDown)[1]
        end
    end
end

# https://github.com/beacon-biosignals/TimeSpans.jl/blob/e3c999021336e51a08d118e6defb792e38ac1cc7/test/runtests.jl#L93-L96
@testset "time_from_index" begin
    @test_throws ArgumentError AlignedSpans.time_from_index(200, 0)
    @test AlignedSpans.time_from_index(100, 1) == Nanosecond(0)
end

@testset "`n_samples`" begin
    times = [Nanosecond(12345), Minute(5), Nanosecond(Minute(5)) + Nanosecond(1),
             Nanosecond(1), Nanosecond(10^6), Nanosecond(6970297031)]
    append!(times, Nanosecond.(Second(1):Second(1):Second(20)))
    sort!(times)

    for rate in (101 // 2, 1001 // 10, 200, 256, 1, 10)
        for i in 1:length(times), j in (i + 1):length(times)
            span = TimeSpan(times[i], times[j])

            if duration(span) < Nanosecond(floor(Int, 10^9 * inv(rate)))
                @test n_samples(rate, TimeSpans.duration(span)) == 0
                continue
            end

            aligned_span = AlignedSpan(rate, span, RoundInward)
            n = n_samples(aligned_span)
            @test n == length(indices(aligned_span))

            # `n` is the *actual* number of samples, given the whole span
            # `n_samples(rate, TimeSpans.duration(span))` gives us the minimum
            # number of samples of any span of that duration.
            # Thus, we should have:
            @test n >= n_samples(rate, TimeSpans.duration(span)) >= (n - 1)
        end
    end
end


text:
using DataFrames
using EcoBase: AbstractAssemblage

"""
### Enumeration of levels that can exist / be calculated for a metacommunity.
"""
@enum DiversityLevel individualDiversity subcommunityDiversity communityDiversity typeDiversity typeCollectionDiversity metacommunityDiversity

"""
### Generates the function to calculate individual diversities

Generates the function to calculate individual diversities for a
series of orders, represented as a vector of qs.

#### Arguments:

- `dm`: DiversityMeasure

#### Returns:

- Function which takes a single number or vector of values of
  parameter q, and returns the individual diversities for those
  values.
"""
individualDiversity

"""
### Generates the function to calculate subcommunity diversity

Generates the function to calculate subcommunity diversity for a
series of orders, represented as a vector of qs.

#### Arguments:

- `dm`: DiversityMeasure

#### Returns:

- Function which takes a single number or vector of values of
  parameter q, and returns the subcommunity diversities for those values.
"""
subcommunityDiversity

"""
### Generates the function to calculate metacommunity diversity

Generates the function to calculate metacommunity diversity for a
series of orders, represented as a vector of qs.

#### Arguments:

- `dm`: DiversityMeasure

#### Returns:

- Function which takes a single number or vector of values of
  parameter q, and returns the metacommunity diversities for those
  values.
"""
metacommunityDiversity

"""
    DiversityMeasure

This type is the abstract supertype of all diversity measure types.
DiversityMeasure subtypes allow you to calculate and cache any kind of
diversity of a metacommunity.
"""
abstract type DiversityMeasure{FP <: AbstractFloat,
                               AbMatrix <: AbstractMatrix,
                               DivArray <: AbstractArray,
                               MC <: AbstractAssemblage} end

"""
    getASCIIName(dm::DiversityMeasure)

Return the ASCII name of the DiversityMeasure

# Arguments:

- `dm`: DiversityMeasure

# Returns:

- String containing simple ASCII name of DiversityMeasure
"""
function getASCIIName(dm::DiversityMeasure)
    s = replace(string(typeof(dm)), "Diversity." => "")
    replace(s, r"{.*}$" => "")
end

"""
    getName(dm::DiversityMeasure)

Return the character corresponding to the DiversityMeasure.

# Arguments:

- `dm`: DiversityMeasure

# Returns:

- String containing unicode (greek) name of DiversityMeasure.
"""
function getName end

"""
    getFullName(dm::DiversityMeasure)

Return the full name of the DiversityMeasure.

# Arguments:

- `dm`: DiversityMeasure

# Returns:

- String containing full descriptive name of DiversityMeasure
"""
function getFullName end

"""
    _getmeta(dm::DiversityMeasure)

Return the metacommunity belonging to the DiversityMeasure.
"""
function _getmeta(dm::DiversityMeasure)
    return dm.meta
end

getsubcommunitynames(dm::DiversityMeasure) = getsubcommunitynames(_getmeta(dm))
gettypenames(dm::DiversityMeasure) = gettypenames(_getmeta(dm))
getdiversityname(dm::DiversityMeasure) = getdiversityname(_getmeta(dm))

(dl::DiversityLevel)(dm::DiversityMeasure) = getPartitionFunction(dm, dl)
(dl::DiversityLevel)(dm::DiversityMeasure, qs) = getPartitionFunction(dm, dl)(qs)

"""
    PowerMeanMeasure

This abstract DiversityMeasure subtype is the supertype of all
diversity measures which are straight power means. PowerMeanMeasure
subtypes allow you to calculate and cache any kind of diversity of a
metacommunity.
"""
abstract type PowerMeanMeasure{FP, AbMatrix, DivArray, MC} <:
    DiversityMeasure{FP, AbMatrix, DivArray, MC} end

"""
    RelativeEntropyMeasure

This abstract DiversityMeasure subtype is the supertype of all
diversity measures which are relative entropy-based diversity measures.
RelativeEntropyMeasure subtypes allow you to calculate and cache any
kind of diversity of a metacommunity.
"""
abstract type RelativeEntropyMeasure{FP, AbMatrix, DivArray, MC} <:
    DiversityMeasure{FP, AbMatrix, DivArray, MC} end

"""
    inddiv(measure::DiversityMeasure, q::Real)
    inddiv(measure::DiversityMeasure, qs::AbstractVector{Real})

Takes a diversity measure and single order or vector of orders, and
returns a DataFrame containing the individual diversities for those values.

# Arguments:

- `dm`: DiversityMeasure
- `q` / `qs`: a single order or a vector of orders

# Returns:

- Returns individual diversities of `dm` for a single order `q` or a
  vector of order `qs`.
"""
function inddiv end

@inline function inddiv(measure::DiversityMeasure, q::Real)
    raw = inddiv_raw(measure, q)
    types = gettypenames(measure)
    scn = getsubcommunitynames(measure)
    scs = reshape(scn, 1, length(scn))
    dfs = broadcast((div, tn, pn) ->
                    DataFrame(div_type=getdiversityname(measure),
                              measure=getASCIIName(measure),
                              q=q,
                              type_level="type", type_name=tn,
                              partition_level="subcommunity",
                              partition_name=pn,
                              diversity=div),
                    raw, types, scs)
    df = reduce(append!, dfs)
    cols = addedoutputcols(_getmeta(measure))
    if length(cols) > 0
        len = length(df)
        data = getaddedoutput(_getmeta(measure))
        for col in keys(cols)
            insert!(df, length(df) + 1, data[col], col)
        end
    end
    return df
end

@inline function inddiv(measure::DiversityMeasure, qs::AbstractVector)
    mapreduce(q -> inddiv(measure, q), append!, qs)
end

@inline function inddiv(meta::AbstractAssemblage, qs)
    mapreduce(dm -> inddiv(dm(meta), qs),
              append!,
              [RawAlpha, NormalisedAlpha,
               RawBeta, NormalisedBeta,
               RawRho, NormalisedRho, Gamma])
end

@inline function inddiv_raw(measure::DiversityMeasure, ::Real)
    measure.diversities
end

"""
    subdiv(measure::DiversityMeasure, q::Real)
    subdiv(measure::DiversityMeasure, qs::AbstractVector{Real})

Takes a diversity measure and single order or vector of orders, and
calculates and returns the subcommunity diversities for those values.

# Arguments:
- `dm`: DiversityMeasure
- `q` / `qs`: a single order or a vector of orders

# Returns:

- Returns subcommunity diversities of `dm` for a single order `q` or a
  vector of order `qs`.
"""
function subdiv end

@inline function subdiv(measure::DiversityMeasure, q::Real)
    raw = subdiv_raw(measure, q)
    scs = getsubcommunitynames(measure)
    dfs = broadcast((div, pn) -> DataFrame(div_type=getdiversityname(measure),
                                           measure=getASCIIName(measure), q=q,
                                           type_level="types", type_name="",
                                           partition_level="subcommunity",
                                           partition_name=pn,
                                           diversity=div),
                    raw, scs)
    df = reduce(append!, dfs)
    cols = addedoutputcols(_getmeta(measure))
    if length(cols) > 0
        len = length(df)
        data = getaddedoutput(_getmeta(measure))
        for col in keys(cols)
            insert!(df, length(df) + 1, data[col], col)
        end
    end
    return df
end

@inline function subdiv(measure::DiversityMeasure, qs::AbstractVector)
    mapreduce(q -> subdiv(measure, q), append!, qs)
end

@inline function subdiv(meta::AbstractAssemblage, qs)
    mapreduce(dm -> subdiv(dm(meta), qs),
              append!,
              [RawAlpha, NormalisedAlpha,
               RawBeta, NormalisedBeta,
               RawRho, NormalisedRho, Gamma])
end

@inline function subdiv_raw(measure::PowerMeanMeasure, q::Real)
    powermean(inddiv_raw(measure, q), one(q) - q, measure.abundances)
end

@inline function subdiv_raw(measure::RelativeEntropyMeasure, q::Real)
    powermean(inddiv_raw(measure, q), q - one(q), measure.abundances)
end

"""
    metadiv(measure::DiversityMeasure, q::Real)
    metadiv(measure::DiversityMeasure, qs::AbstractVector{Real})

Takes a diversity measure and single order or vector of orders, and
calculates and returns the metacommunity diversities for those values.

# Arguments:

- `dm`: DiversityMeasure
- `q` / `qs`: a single order or a vector of orders

# Returns:

- Returns metacommunity diversities of `dm` for a single order `q` or a
  vector of order `qs`.
"""
function metadiv end

@inline function metadiv(measure::DiversityMeasure, q::Real)
    raw = metadiv_raw(measure, q)
    df = DataFrame(div_type=getdiversityname(measure),
                   measure=getASCIIName(measure), q=q,
                   type_level="types", type_name="",
                   partition_level="metacommunity",
                   partition_name="",
                   diversity=raw)
   cols = addedoutputcols(_getmeta(measure))
   if length(cols) > 0
       len = length(df)
       data = getaddedoutput(_getmeta(measure))
       for col in keys(cols)
           insert!(df, length(df) + 1, data[col], col)
       end
   end
   return df
end

@inline function metadiv(measure::DiversityMeasure, qs::AbstractVector)
    mapreduce(q -> metadiv(measure, q), append!, qs)
end

@inline function metadiv(meta::AbstractAssemblage, qs)
    mapreduce(dm -> metadiv(dm(meta), qs),
              append!,
              [RawAlpha, NormalisedAlpha,
               RawBeta, NormalisedBeta,
               RawRho, NormalisedRho, Gamma])
end

@inline function metadiv_raw(measure::DiversityMeasure, q::Real)
    powermean(subdiv_raw(measure, q), one(q) - q, measure.weights)
end

function getPartitionFunction(measure::DiversityMeasure, level::DiversityLevel)
    if (level == individualDiversity)
        return function (qs)
            inddiv(measure, qs)
        end
    elseif (level == subcommunityDiversity)
        function (qs)
            subdiv(measure, qs)
        end
    elseif (level == metacommunityDiversity)
        function (qs)
            metadiv(measure, qs)
        end
    else
        error("Unrecognised diversity level")
    end
end

"""
    RawAlpha

Calculates raw alpha diversity (α) of all of the individuals in a
metacommunity, and caches them for subsequent analysis. This is a
subtype of PowerMeanMeasure, meaning that all composite diversity
measures are simple powermeans of the individual measures.

#### Constructor arguments:

- `meta`: a Metacommunity
"""
struct RawAlpha{FP, AbMatrix, DivArray, MC} <:
    PowerMeanMeasure{FP, AbMatrix, DivArray, MC}
    abundances::AbMatrix
    weights::Vector{FP}
    diversities::DivArray
    meta::MC
end

function RawAlpha(meta::M) where M <: AbstractAssemblage
    ab = getabundance(meta)
    ws = getweight(meta)
    value = getordinariness!(meta) .^ -1
    return RawAlpha{eltype(ab), typeof(ab),
                    typeof(value), M}(ab, ws, value, meta)
end

getName(::RawAlpha) = "α"
getFullName(::RawAlpha) = "raw alpha diversity"

"""
    NormalisedAlpha

Calculates normalised alpha diversity (ᾱ) of all of the individuals in
a metacommunity, and caches them for subsequent analysis. This is a
subtype of PowerMeanMeasure, meaning that all composite diversity
measures are simple powermeans of the individual measures.

#### Constructor arguments:

- `meta`: a Metacommunity
"""
struct NormalisedAlpha{FP, AbMatrix, DivArray, MC} <:
    PowerMeanMeasure{FP, AbMatrix, DivArray, MC}
    abundances::AbMatrix
    weights::Vector{FP}
    diversities::DivArray
    meta::MC
end

function NormalisedAlpha(meta::M) where M <: AbstractAssemblage
    ab = getabundance(meta)
    ws = getweight(meta)
    value = ws' ./ getordinariness!(meta)
    return NormalisedAlpha{eltype(ab), typeof(ab),
                           typeof(value), M}(ab, ws, value, meta)
end

getName(::NormalisedAlpha) = "ᾱ"
getFullName(::NormalisedAlpha) = "normalised alpha diversity"

"""
    RawBeta

Calculates distinctiveness (β, raw beta diversity) of all of the individuals in a
metacommunity, and caches them for subsequent analysis. This is a
subtype of RelativeEntropyMeasure, meaning that subcommunity and type
composite diversity measures are relative entropies, and their
composite types are powermeans of those measures.

#### Constructor arguments:

- `meta`: a Metacommunity
"""
struct RawBeta{FP, AbMatrix, DivArray, MC} <:
    RelativeEntropyMeasure{FP, AbMatrix, DivArray, MC}
    abundances::AbMatrix
    weights::Vector{FP}
    diversities::DivArray
    meta::MC
end

function RawBeta(meta::M) where M <: AbstractAssemblage
    ab = getabundance(meta)
    ws = getweight(meta)
    value = getordinariness!(meta) ./ getmetaordinariness!(meta)
    return RawBeta{eltype(ab), typeof(ab),
                   typeof(value), M}(ab, ws, value, meta)
end

const Distinctiveness = RawBeta

getName(::RawBeta) = "β"
getFullName(::RawBeta) = "distinctiveness"

"""
    NormalisedBeta

Calculates normalised beta diversity (β̄) of all of the individuals in
a metacommunity, and caches them for subsequent analysis. This is a
subtype of RelativeEntropyMeasure, meaning that subcommunity and type
composite diversity measures are relative entropies, and their
composite types are powermeans of those measures.

#### Constructor arguments:

- `meta`: a Metacommunity
"""
struct NormalisedBeta{FP, AbMatrix, DivArray, MC} <:
    RelativeEntropyMeasure{FP, AbMatrix, DivArray, MC}
    abundances::AbMatrix
    weights::Vector{FP}
    diversities::DivArray
    meta::MC
end

function NormalisedBeta(meta::M) where M <: AbstractAssemblage
    ab = getabundance(meta)
    ws = getweight(meta)
    value = getordinariness!(meta) ./ (getmetaordinariness!(meta) .* ws')
    return NormalisedBeta{eltype(ab), typeof(ab),
                          typeof(value), M}(ab, ws, value, meta)
end

getName(::NormalisedBeta) = "β̄"
getFullName(::NormalisedBeta) = "effective number of subcommunities"

"""
    RawRho

Calculates redundancy (ρ, raw beta diversity) of all of the
individuals in a metacommunity, and caches them for subsequent
analysis. This is a subtype of PowerMeanMeasure, meaning that all
composite diversity measures are simple powermeans of the individual
measures.

#### Constructor arguments:

- `meta`: a Metacommunity
"""
struct RawRho{FP, AbMatrix, DivArray, MC} <:
    PowerMeanMeasure{FP, AbMatrix, DivArray, MC}
    abundances::AbMatrix
    weights::Vector{FP}
    diversities::DivArray
    meta::MC
end

function RawRho(meta::M) where M <: AbstractAssemblage
    ab = getabundance(meta)
    ws = getweight(meta)
    value = getmetaordinariness!(meta) ./ getordinariness!(meta)
    return RawRho{eltype(ab), typeof(ab),
                  typeof(value), M}(ab, ws, value, meta)
end

const Redundancy = RawRho

getName(::RawRho) = "ρ"
getFullName(::RawRho) = "redundancy"

"""
    NormalisedRho

Calculates redundancy (ρ̄, normalised beta diversity) of all of the
individuals in a metacommunity, and caches them for subsequent
analysis. This is a subtype of PowerMeanMeasure, meaning that all
composite diversity measures are simple powermeans of the individual
measures.

#### Constructor arguments:

- `meta`: a Metacommunity
"""
struct NormalisedRho{FP, AbMatrix, DivArray, MC} <:
    PowerMeanMeasure{FP, AbMatrix, DivArray, MC}
    abundances::AbMatrix
    weights::Vector{FP}
    diversities::DivArray
    meta::MC
end

function NormalisedRho(meta::M) where M <: AbstractAssemblage
    ab = getabundance(meta)
    ws = getweight(meta)
    value = (getmetaordinariness!(meta) .* ws') ./ getordinariness!(meta)
    return NormalisedRho{eltype(ab), typeof(ab),
                         typeof(value), M}(ab, ws, value, meta)
end

const Representativeness = NormalisedRho

getName(::NormalisedRho) = "ρ̄"
getFullName(::NormalisedRho) = "representativeness"

"""
    Gamma

Calculates gamma diversity (γ) of all of the individuals in a
metacommunity, and caches them for subsequent analysis. This is a
subtype of PowerMeanMeasure, meaning that all composite diversity
measures are simple powermeans of the individual measures.

#### Constructor arguments:

- `meta`: a Metacommunity
"""
struct Gamma{FP, AbMatrix, DivArray, MC} <:
    PowerMeanMeasure{FP, AbMatrix, DivArray, MC}
    abundances::AbMatrix
    weights::Vector{FP}
    diversities::DivArray
    meta::MC
end

function Gamma(meta::M) where M <: AbstractAssemblage
    ab = getabundance(meta)
    ws = getweight(meta)
    value = fill!(similar(ws), 1)' ./ getmetaordinariness!(meta)
    return Gamma{eltype(ab), typeof(ab), typeof(value), M}(ab, ws, value, meta)
end

getName(::Gamma) = "γ"
getFullName(::Gamma) = "gamma diversity"

RecipesBase.@recipe function f(var::Tuple{<: DiversityMeasure,
                                          <: Real})
    title := getFullName(var[1]) * " (q = $(var[2])) - " * getdiversityname(var[1]) *
        " diversity"
    colorbar_title := getASCIIName(var[1])
    subdiv(var...)[:diversity], getcoords(places(_getmeta(var[1])))
end


text:
# Before running this, please make sure to activate and instantiate the environment
# corresponding to [this `Project.toml`](https://raw.githubusercontent.com/alan-turing-institute/DataScienceTutorials.jl/master/Project.toml) and [this `Manifest.toml`](https://raw.githubusercontent.com/alan-turing-institute/DataScienceTutorials.jl/master/Manifest.toml)
# so that you get an environment which matches the one used to generate the tutorials:
#
# ```julia
# cd("DataScienceTutorials") # cd to folder with the *.toml
# using Pkg; Pkg.activate("."); Pkg.instantiate()
# ```

# Before running this, please make sure to activate and instantiate the environment# corresponding to [this `Project.toml`](https://raw.githubusercontent.com/alan-turing-institute/MLJTutorials/master/Project.toml) and [this `Manifest.toml`](https://raw.githubusercontent.com/alan-turing-institute/MLJTutorials/master/Manifest.toml)# so that you get an environment which matches the one used to generate the tutorials:## ```julia# cd("MLJTutorials") # cd to folder with the *.toml# using Pkg; Pkg.activate("."); Pkg.instantiate()# ```
# Classification of fraudulent/not credit card transactions (imbalanced data)# By Kristian Bjarnason. The original script can be found [here](https://github.com/kbjarnason/credit-card-fraud-classification)# This version implements train, test, val split rather than just train, test.
# ***
# using Pkg; Pkg.add("Revise")# Revise wasn't installed - so installed it, but it triggers update to instantiatd environment (Project.toml,Manifest.toml)
#Relevant packages
using Revise, DataFrames, CSV, LinearAlgebra, Dates, Statistics, MLJ, MLJBase, MLJModels, MLJLinearModels, Plots, Flux
using UrlDownload, MLBase, StatsBase# , ROC, EvalCurves, ==> problem with these last two packages
using Flux:outdims, activations, @epochs, throttle
using Flux.Data

# using AUC # add git@github.com:paulstey/AUC.jl.git
# ***
# ## *Data Preparation*
# Divide the sample into two equal sub-samples. Keep the proportion of frauds the same in each sub-sample (246 frauds in each).# Use one sub-sample to estimate (train) your models and the second one to evaluate the out-of-sample performance of each model.
# ***
# Import data
data = urldownload("https://storage.googleapis.com/download.tensorflow.org/data/creditcard.csv")
data = DataFrame(data)

# Let's inspect the (types of) variables contained in the DataFrame
schema(data)

# The Time column is not relevant to our analysis, we drop it.
select!(data, Not(:Time))

# Let's get a summary of the remaining data.
describe(data)

# Note that the Amount variable spans a wide range of values.# To reduce variation in the data, we take the natural logarithm. Note that some values are 0 and that log(0) will raise an error. We add 1e-6 so no values are 0 prior to being log transformed.# We transform in place using '!'
data[!,:Amount] = log.(data[!,:Amount] .+ 1e-6)

# ***
# Let's unpack the dataframe and create separate frames for our target variable and features.
y, X = unpack(data, ==(:Class), col -> true)

# And partition between training and test observations
train, test = partition(eachindex(y), 0.8, shuffle=true, rng=111)

# Important: train and test are one dimensional arrays and respectively contain the row indices of the train and test sets' observations
# Setup train and test arrays/vectors
Xtrain = selectrows(X, train)
Xtest = selectrows(X, test)
ytrain = selectrows(y, train)
ytest = selectrows(y, test)

# Let's also create categorical versions of the target variable, for use in the LogisticClassifier
ytrain_cat = categorical(selectrows(y, train))
ytest_cat = categorical(selectrows(y, test))

# Create standardised features for SVM and NN
stand_model = Standardizer()
std_m_fit = fit!(machine(stand_model, Xtrain))
Xtrain_std = MLJModels.transform(std_m_fit, Xtrain)
Xtest_std = MLJModels.transform(std_m_fit, Xtest)

# ***
# ##Estimation of models##
# Estimation of three different models:# 1. logit# 2. support vector machines# 3. neural network.
# ### Logit
# **Initial logit classification with lambda = 1.0**
@load LogisticClassifier pkg=MLJLinearModels

model_logit = LogisticClassifier(lambda=1.0)
m_logit = machine(model_logit, Xtrain, ytrain_cat)
fit!(m_logit)

# **Predictions**
# The prediction from the LogisticClassifier is a probabilistic output, i.e. for each observation in the sample it attaches a probability to each of the possible values of the target.# To recover a deterministic output, we take the mode of the distribution.
yhat_logit_p = MLJBase.predict(m_logit, Xtest)
yhat_logit = categorical(mode.(yhat_logit_p))

# How does this model perform?
@show confusion_matrix(yhat_logit, ytest_cat)
@show misclassification_rate(yhat_logit, ytest_cat)

# (Looks like it's not too bad)# Let' see if we can do even better by doing a little tuning.
#CSV.write("yhat_logit.csv", yhat_logit)
#CSV.write("yhat_logit_p.csv", yhat_logit_p)

# ***
# ### Tuned logit
# Still LogisticClassifier but implementing hyperparameter tuning.
model_logit = @load LogisticClassifier pkg=MLJLinearModels
r = range(model_logit, :lambda, lower=1e-6, upper=100, scale=:log)

# The tuning strategy is specified below; note the use of a probabilistic measure (cross_entropy).
self_tuning_logit_model = TunedModel(model=model_logit,
                                                  resampling = CV(nfolds=3),
                                                  tuning = Grid(resolution=10),
                                                  range = r,
                                                  measure = cross_entropy)

self_tuning_logit = machine(self_tuning_logit_model, Xtrain, ytrain_cat)
fit!(self_tuning_logit)

# Predictions
yhat_logit_tuned_p = MLJBase.predict(self_tuning_logit, Xtest)
yhat_logit_tuned = categorical(mode.(yhat_logit_tuned_p))

# Let's take a look at the misclassification_rate. It is, surprisingly, slightly higher than the one calculated for the non tuned model.
@show misclassification_rate(yhat_logit_tuned, ytest_cat)

# ### Support Vector Machine
# ***
# Initial SVM classification with cost = 1.0
@load SVC

# To fit the SVM, we declare a pipeline which comprises both a standardizer and the model
model_svm = @pipeline Std_SVC(std_model = Standardizer(),
                                      svc = SVC())

svc = machine(model_svm, Xtrain, ytrain_cat)

fit!(svc)

yhat_svm = MLJBase.predict(svc, Xtest_std)
#0.00163
@show confusion_matrix(yhat_svm, ytest_cat)
@show misclassification_rate(yhat_svm, ytest_cat)

# CSV.write("yhat_svm.csv", yhat_svm)
# ***
# Tuned SVM
#@load SVC #does this command need to be repeated?

model_svm = @pipeline Std_SVC(std_model = Standardizer(),
                              svc = SVC())

r = range(model_svm, :(svc.cost), lower=0.0, upper=2.5, scale=:linear)
iterator(r,6)

self_tuning_svm_model = TunedModel(model=model_svm,
                                   resampling = CV(nfolds=3),
                                   tuning = Grid(resolution=6),
                                   range = r,
                                   measure = MLJ.precision)

self_tuning_svm = machine(self_tuning_svm_model, Xtrain, ytrain_cat)

fit!(self_tuning_svm)

@show fitted_params(self_tuning_svm).best_model
@show report(self_tuning_svm)

yhat_svm_tuned = MLJBase.predict(self_tuning_svm, Xtest_std)

@show misclassification_rate(yhat_svm_tuned, ytest_cat)
@show confusion_matrix(yhat_svm_tuned, ytest_cat)

#CSV.write("yhat_svm_tuned.csv", yhat_svm_tuned)

# ### Neural Network
#oversample fraudulent cases (since data so imbalanced)
X_train_oversampled = vcat(X_train,repeat(data_fraud[1:Int(nrow(data_fraud)/2),1:29], 100))
y_train_oversampled = vcat(y_train,repeat(data_fraud[1:Int(nrow(data_fraud)/2),30], 100))

stand_model = Standardizer()
X_train_oversampled_std = MLJModels.transform(fit!(machine(stand_model, X_train_oversampled)), X_train_oversampled)

data1 = DataLoader(Array(X_train_oversampled_std)', y_train_oversampled, batchsize=2048, shuffle=true)

n_inputs = ncol(X_train_oversampled)
n_outputs = 1
n_hidden1 = 16

m = Chain(
          Dense(n_inputs, n_hidden1, relu),
          Dropout(0.1),
          Dense(n_hidden1, n_outputs, σ)
          )

loss(x, y) = Flux.tversky_loss(m(x), y, β=0.9) #tversky loss uses precision and recall, slower calc than crossentropy

# loss(x, y) = Flux.crossentropy(m(x), y)
ps = Flux.params(m)

# opt = ADAM()
opt = Descent()

@epochs 50 Flux.train!(loss, ps, data1, opt)

yhat_nn_p = vec(m(Array(X_test_std)'))
yhat_nn = categorical(Int.(yhat_nn_p .<= 0.5))

cm_nn = confusion_matrix(yhat_nn, y_test)
misclassification_rate(yhat_nn, y_test)

CSV.write("yhat_nn.csv", yhat_nn)

#%%md
*OOS results*
#%%
#if needed can reload from here
yhat_logit_tuned = CSV.read("yhat_logit_tuned.csv")
yhat_svm_tuned = CSV.read("yhat_svm_tuned.csv")
yhat_nn = CSV.read("yhat_nn.csv")

#%%md
Misclassification rate
#%%
misclassification_rate(yhat_logit_tuned, y_test)
misclassification_rate(yhat_svm_tuned, y_test)
misclassification_rate(yhat_nn, y_test)

#%%md
Confusion matrix
#%%
cm_logit = confusion_matrix(yhat_logit_tuned, y_test)
cm_svm = confusion_matrix(yhat_svm_tuned , y_test)
cm_nn = confusion_matrix(yhat_nn, y_test)

#%%md
ROC curves
#%%
#Due to different data output structure had to use different packages for ROC curves
plot(roc_curve(yhat_logit_tuned_p, y_test))

# plot(ROC.roc(pdf.(yhat_logit_tuned_p), y_test, 1)))
plot(ROC.roc(yhat_nn_p, y_test, 1))

# don't have score vectors for SVM# plot(roc_curve(yhat_svm_p,y_test))# plot(roc_curve(yhat_nn_p,y_test))
#how to plot this??
MLBase.roc(y_test_int, yhat_nn_p)

#%%md
Precision-Recall curve
#%%
plot(prcurve(pdf.(yhat_logit_tuned_p,1), y_test_int))
plot(prcurve(yhat_nn_p, y_test_int))

# don't have score vectors for SVM# plot(prcurve(pdf.(yhat_svm_p,1), y_test_int))
#END

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
using DiffEqBiological
using Random
using Distributions

prokaryote_model = @reaction_network ProkaryoteRN begin
    # order of appearnce: DNA, P₂, DNAxP₂, RNA, P
    θ1, DNA + P₂ --> DNAxP₂
    θ2, DNAxP₂ --> DNA + P₂
    θ3, DNA --> DNA + RNA
    θ4, RNA --> RNA + P
    θ5, 2P --> P₂
    θ6, P₂ --> 2*P
    θ7, RNA --> 0*RNA
    θ8, P --> 0*P
end θ1 θ2 θ3 θ4 θ5 θ6 θ7 θ8

function simulate_data(u0, tt, rates)
    prob = DiscreteProblem(u0, tt, rates)
    jump_prob = JumpProblem(prob, Direct(), prokaryote_model)
    sol = solve(jump_prob, SSAStepper())
    sol
end

u0 = [5,8,5,8,8]
θ = [0.1, 0.7, 0.35, 0.2, 0.1, 0.9, 0.3, 0.1]
tt = (0.0, 100.0)

Random.seed!(4)
sol = simulate_data(u0, tt, θ)
obs_time = collect(tt[1]:1:tt[end])
# permute to obtain (RNA, P, P₂, DNA, DNAxP₂) ordering
XX = map(x->x[[4,5,2,1,3]], sol.(obs_time))

# Let's add noise
#obs = [o[2] + 2*o[3] for o in XX[2:end]] .+ rand(Normal(0.0, 2.0), length(obs)-1)
# Alternative
obs = [[o[1] + rand(Normal(0.0, 1.0)), o[2] + 2*o[3]+rand(Normal(0.0, 2.0))] for o in XX[2:end]]

# save the data
OUT_DIR = joinpath(Base.source_dir(), "..", "..", "output9")

# save the process
open(joinpath(OUT_DIR, "prokaryote_custom.dat"), "w") do f
    write(f, "time RND P P2 DNA DNAxP2\n")
    for i in 1:length(obs_time)
        data = join(obs[i], " ")
        line = join([obs_time[i], " ", data, "\n"])
        write(f, line)
    end
end

# save just partial and noisy Observations
open(joinpath(OUT_DIR, "prokaryote_custom.dat"), "w") do f
    write(f, "time P+2*P2\n")
    for i in 2:length(obs_time)
        data = join(obs[i-1], " ")
        line = join([obs_time[i], " ", data, "\n"])
        write(f, line)
    end
end

# alternative save
open(joinpath(OUT_DIR, "prokaryote_custom.dat"), "w") do f
    write(f, "time RNA P+2*P2\n")
    for i in 2:length(obs_time)
        data = join(obs[i-1], " ")
        line = join([obs_time[i], " ", data, "\n"])
        write(f, line)
    end
end

# Ignore below....

Random.seed!(4)
θˣ = [0.1, 0.7, 0.35, 0.2, 0.1, 0.9, 0.3, 0.1]
K = 10.0
P˟ = Prokaryote(θ_init..., K)
x0, dt, T = ℝ{4}([8.0,8.0,8.0,5.0]), 1/5000, 100.0

XX, _ = simulate_segment(ℝ{8}(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
                        x0, P˟, 0.0:dt:T)

skip = 5000

obs_time, obs_vals = XX.tt[1:skip:end], [o[2] + 2*o[3]+rand(Normal(0.0, 2.0)) for o in XX.yy[1:skip:end]]

obs_time
obs_vals

# save just partial and noisy Observations
open(joinpath(OUT_DIR, "prokaryote_custom.dat"), "w") do f
    write(f, "time P+2*P2\n")
    for i in 2:length(obs_time)
        data = join(obs_vals[i-1], " ")
        line = join([obs_time[i], " ", data, "\n"])
        write(f, line)
    end
end


text:
function check_samptypes(coms::Vector{ImagineSignal}, rig::AbstractString)
    #check that all analog and digital entries are acceptable for this rig
    tm = rigtemplate(rig)
    chans0 = map(daq_channel, tm)
    chans = map(daq_channel, coms)
    for chan in chans
        if intervals(getdaqchan(coms, chan)) != intervals(getdaqchan(tm, chan))
            error("The sample value intervals for channel '$chan' are incompatible with the $rig rig.")
        end
    end
    #check that all have equal sample rate
    rs = map(samprate,coms)
    if !all(rs.==rs[1])
        error("All commands must use equal sample rates.  This can be set per-channel with `set_samprate!`")
    end
    #check that all have equal number of samples
    nsamps = map(length, getoutputs(coms))
    if !all(nsamps.==nsamps[1])
        error("All output commands must have an equal number of samples.  Check this with `length(com::ImagineSignal)`")
    end
    return true
end

#In order to be "sufficient" a command must include a positioner trace, at least one camera trace, and at least one laser trace
function check_sufficiency(coms::Vector{ImagineSignal})
    if length(findcameras(coms)) == 0
        error("No camera commands were found")
    end
    if length(findpositioners(coms)) == 0
        error("No positioner commands were found")
    end
    if length(findlasers(coms)) == 0
        error("No laser commands were found")
    end
end

#Check whether user has modified fixed names
function check_fixed_names(coms::Vector{ImagineSignal}, rig::String)
    fixed_names = FIXED_NAMES[rig]
    name_lookup = DEFAULT_DAQCHANS_TO_NAMES[rig]
    for c in coms
        nm = name(c)
        dc = daq_channel(c)
        if nm != name_lookup[dc] && in(name_lookup[dc], fixed_names)
            error("DAQ Channel $(dc) must have the name $(name_lookup[dc]) but instead it is named $nm")
        end
    end
    return true
end

#checks that all commands in the vector have the same rig, and the rig is recognized
function check_rig_names(coms::Vector{ImagineSignal})
    rig = rig_name(coms[1])
    if !all(map(rig_name, coms) .== rig)
        error("The set of commands to be written must all be targeted to the same rig")
    end
    if !in(rig, RIGS)
        error("Unrecognized rig: $rig")
    end
    return true
end

#Check for invalid DAQ channel names
function check_valid_channels(coms::Vector{ImagineSignal}, rig::String)
    name_lookup = DEFAULT_DAQCHANS_TO_NAMES[rig]
    for c in map(daq_channel, coms)
        if !haskey(name_lookup, c)
            error("The channel $c is not accessible on this rig")
        end
    end
    nms = map(name, coms)
    chan_nms = map(daq_channel, coms)
    if length(unique(nms)) != length(nms)
        error("Found one or more duplicate channel names")
    end
    if length(unique(chan_nms)) != length(chan_nms)
        error("Found one or more duplicate DAQ channel identifiers")
    end
    return true
end

function validate_group(sigs::Vector{ImagineSignal}; check_is_sufficient = true)
    check_rig_names(sigs)
    rig = rig_name(first(sigs))
    check_valid_channels(sigs, rig)
    check_fixed_names(sigs, rig)
    if check_is_sufficient
        check_sufficiency(sigs)
    end
    check_samptypes(sigs, rig)
    return true
end

max_forks_repo_name: HolyLab/ImagineInterface


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

using Markdown
using InteractiveUtils

# ╔═╡ f80fd13d-207a-458a-b326-77256fe946e0
using Plots

# ╔═╡ bda0f8d4-5ab0-4c65-9d79-71ccbfe8dd3a
using SpecialFunctions

# ╔═╡ b608f792-dd91-4bd2-97ab-baac7b7d2f9c
using Dierckx

# ╔═╡ 240e07ba-1508-485c-8fe1-d33e8fa52804
using SpecialFunctions

# ╔═╡ 4e16dbbe-b0fc-11eb-2349-ad7ccd0ff56b
md"""
# Simple questions and simple answers about the Boys function
Rick Muller

A large part of electronic structure theory lies in computing the repulsive
interaction between electrons, and a large part of electron repulsion
is captured in the Boys function that computes the repulsion between
partially overlapping electronic wave functions.

The Boys function is given by [^boys50]

$Fm(T)=\int_0^1 u^{2m}\exp(-Tu^2)du$

and was introduced by Boys in 1950. Since that time, a surprising amount
of work has gone into efficiently computing [^gjp91] this function. The
best of these methods [^gjp91] divides the integration over $u$ into multiple slices
and approximates the integrand in this region using Taylor or Chebyshev
polynomials. Efficient programs [^libint] to evaluate these terms have detailed
tables with nodes and weights for the interpolation and quadrature schemes.

One says "surprising" because we very much know what electron repulsion
looks like in crude terms: at long distances, the interaction should look like Coulombic $1/r$ repulsion, and at shorter distances the interaction rolls
over to a constant value with zero slope. This behavior was approximated by
Mataga and Nishimoto [^mn57] as

$F_{MN}(r) \sim \frac{1}{a+r}$

which we plot below.

"""

# ╔═╡ 8a736133-17f2-446d-9ab5-4bc46bdd556b
fmn(r,a=1) = 1/(a+r);

# ╔═╡ bf014265-210c-412f-9af4-5a654f557de9
begin
	rs = 0:100
	mnvals =
	plot(rs,[fmn(r,1) for r in rs],label="fmn(r,a=1)",
		title="Behavior of the Mataga-Nishimoto approximation",
		xlabel="R (au)")
end

# ╔═╡ d8061039-367c-42ee-8a65-a32dca37fc67
md"""
We can use the following simple definition for the Boys function taken from
[^libint]:
"""

# ╔═╡ bc53663d-1f9f-43a6-8140-e8eeac1618bb
"Reference Fm from libint"
function fboys(m,T,eps = 1e-10)
    denom = (m + 0.5)
    term = exp(-T) /2denom
    old_term = 0.0
    sum = term
    while true
        denom += 1
        old_term = term
        term = old_term * T / denom
        sum += term
        term > sum*eps || old_term < term || break
    end
    return sum
end;

# ╔═╡ 1a0bffb2-8656-44db-b3b3-6cf6b2494b7c
md"This is the version that uses the gamma function:"

# ╔═╡ b6b757ad-7fca-46c5-b42f-20f94ae4cbd7
"gammainc returns the lower incomplete gamma function"
gammainc(a,x) = gamma(a)*gamma_inc(a,x)[1]

# ╔═╡ 145a0007-b8d7-46a1-9cb8-ee063057527c
"Boys Fgamma function, using the lower incomplete gamma function."
function Fgamma(m,x,SMALL=1e-18,Tcrit=20.0)
    # Note, most programs use a much larger value for Tcrit (117)
    mhalf = m+0.5
    x = max(x,SMALL) # Evidently needs underflow protection
    if x>Tcrit
        retval = sqrt(pi/2)*factorial2(2m-1)/(2x)^mhalf
    else
        retval = 0.5*x^-mhalf*gammainc(mhalf,x)
    end
    return retval
end

# ╔═╡ f6fa9869-b454-458d-95e3-818a9112cc52
md"""
The asymptotic behavior of the Boys function for large $T$ is known
to be

$F_\mathrm{asymp}(T,m) = \sqrt{\frac{\pi}{2}}(2m-1)!!/(2T)^{m+0.5}$

which leads to a slightly more complicated approximation form of

$F_\mathrm{approx}(T)=\sqrt{\frac{\pi}{2}}\frac{1}{a\sqrt{\pi/2}+\sqrt{2T}}$
"""

# ╔═╡ e7454913-e406-4b4a-856f-47103a7bf23d
fasymp(m,T) = sqrt(0.5*pi)*prod(2m-1:-2:1)*(2T)^(-m-0.5);

# ╔═╡ cdc2c6c9-e567-4ee6-a002-8fea660e35fa
function fapprox(T,a=1)
	pre = sqrt(0.5*pi)
	return pre/(a*pre + (2T)^(0.5))
end;

# ╔═╡ b5f42a47-eb6c-4b59-9827-a07b81618018
md"""
This approximation leads to good agreement with the Boys function.
"""

# ╔═╡ 105717a4-f481-479d-bfbb-dc35cc3c0997
begin
	plot(rs,[fboys(0,r) for r in rs],label="fboys(0,r)")
	plot!(rs,[fapprox(r,1) for r in rs],label="fapprox(r,a=1)")
	#title!("Parameter-free approximation of Boys function")
end

# ╔═╡ 0b513de7-cca5-45ca-a4fa-2f2604eb7da9
md"""
We note that although the behavior of `fapprox` plotted above is not exact,
it is essentially parameter-free when `a=1`, and thus represents a very efficient starting point for further analysis.
"""

# ╔═╡ 4dd0b1e6-61f0-4daa-80f1-625dc793b33a
md"""
As an aside, we certainly don't mean to denigrate the excellent work that has
gone into the development of the approximations methods [^gjp91] or their
efficient implementation [^libint]. These are remarkably accurate and fast
methods that have been derived with great thought and applied repeatedly to
chemical calculations with excellent results. Our intent here is to merely
understand whether a more intuitive approach is possible to this admittedly
solved problem.
"""

# ╔═╡ fdd2691f-f26f-4c26-8cc6-0d8edaa5278a
md"""
Instead of

$F_\mathrm{approx}(T)=\sqrt{\frac{\pi}{2}}\frac{1}{a\sqrt{\pi/2}+\sqrt{2T}}$

we can expand the denominator in a term that will go to zero at large r:

$F_\mathrm{approx2}(T)=\sqrt{\frac{\pi}{2}}\frac{1}{a\exp(-bT)\sqrt{\pi/2}+\sqrt{2T}}$
"""

# ╔═╡ a0a09c9f-c6d2-422f-b32a-8e3f8131b3c3
function fapprox2(T,a=1,b=1)
	pre = sqrt(0.5*pi)
	return pre/(a*pre*exp(-b*T) + (2T)^(0.5))
end;

# ╔═╡ ae9658a3-9241-423d-9daa-d6644e8a206a
begin
	rs2=0:0.1:15
	plot(rs2,[fboys(0,r) for r in rs2],label="fboys(0,r)",ylim=(0,1))
	plot!(rs2,[fasymp(0,r) for r in rs2],label="fasymp(0,r)")
	plot!(rs2,[fapprox(r,1) for r in rs2],label="fapprox(r,1)")
	plot!(rs2,[fapprox2(r,1,1.5) for r in rs2],label="fapprox2(r,a=1,b=1.5)")
	#title!("Parameter-free approximation of Boys function")
end

# ╔═╡ 8c8adf9b-5eff-4942-9a49-95a4aa1e2002
md"""
With very few parameters (essentially just estimating the `b=1.5` fit) it's possible to get a very good reproduction of the Boys function.

If rather than $exp(-bT)$ we assume there is some function of $T$ that multiplies
the decay of the constant term, what would that function look like?
"""

# ╔═╡ 9428670f-5cf7-4c9e-bb93-dd0747a34182
boysdenom(T) = 1/fboys(0,T) - sqrt(2T)/sqrt(pi/2);

# ╔═╡ 771d4eeb-72e6-4863-a114-ff203a42455b
begin
	rs3 = 0:0.1:3
	plot(rs3,boysdenom.(rs3),label="Boys denom decay")
	plot!(rs3,[exp(-T-0.56) for T in rs3],label="-T-0.56")
	plot!(rs3,[exp(-4T) for T in rs3],label="-4T")
	#plot!(rs3,[exp(-4T)/2+exp(-T-0.56)/2 for T in rs3],label="avg")
end

# ╔═╡ 997e6018-bb0e-4933-8b13-b904e96fc819
md"Plot this as `fapprox3`"

# ╔═╡ 894135da-1ef9-4170-aa49-1b6307e398eb
function fapprox3(T,a=1)
	pre = sqrt(0.5*pi)
	return pre/(a*pre*exp(-T-0.56) + (2T)^(0.5))
end;

# ╔═╡ 9dfb09bb-22d8-44f2-9c4b-98777e43de64
begin
	plot(rs3,[fboys(0,r) for r in rs3],label="fboys(0,r)")
	plot!(rs3,[fapprox2(r) for r in rs3],label="fapprox2(r)")
	plot!(rs3,[fapprox3(r) for r in rs3],label="fapprox3(r)")
end

# ╔═╡ d526caa1-ff7e-4ec5-ab59-2a10e9ad9c5e
md"""
## Plot of multiple Fm(T) values

Plot the family of m values on a semilog-y plot. Very regular behavior:
"""

# ╔═╡ 3923126f-93be-48f2-bd37-3ef50588ff47
begin
	rs4 = 0:0.1:10
	plot(rs4,fboys.(0,rs4),yaxis=:log,legend=false)
	plot!(rs4,Fgamma.(0,rs4))
	plot!(rs4,fboys.(1,rs4))
	plot!(rs4,fboys.(2,rs4))
	plot!(rs4,fboys.(3,rs4))
	plot!(rs4,fboys.(4,rs4))
	plot!(rs4,fboys.(5,rs4))
	plot!(rs4,fboys.(6,rs4))
end

# ╔═╡ eb05a798-8eb2-4ff1-b88f-b58babe64891
denom(T,m) = 1/fboys(m,T);

# ╔═╡ 84df8774-a1f6-49f7-b84f-c1482b51dd9f
denomg(T,m) = 1/Fgamma(m,T);

# ╔═╡ 08c5dc9b-0cb1-4bea-bd97-4872f46e63bf
begin
	let rs4 = 0:0.1:5
	plot(rs4,denom.(0,rs4),yaxis=:log,legend=false)
	plot!(rs4,denom.(1,rs4))
	plot!(rs4,denom.(2,rs4))
	plot!(rs4,denom.(3,rs4))
	plot!(rs4,denom.(4,rs4))
	plot!(rs4,denom.(5,rs4))
	plot!(rs4,denom.(6,rs4))
	end
end

# ╔═╡ 347ca2e4-b932-4451-a12e-15ca7e23d052
md"""
## Interpolation
The other option that is appealing is to simply tabulate the values of the function
for the m-values that matter [0,12] on a relatively fine mesh and then interpolate
them either linearly or using Interpolate.jl.

Tasks:
- Express Fm(T) on a mesh
	- Decide on linear or log
	- Might have to be linear, since we definitely want values at T=0.
- Write a function to compute the maximum error with that discretization.
"""

# ╔═╡ b15df2d5-f407-4d2d-959f-18b21c8a3b60
ts = 0:0.01:100 # Intentially too coarse a mesh

# ╔═╡ ae817c36-d9b8-4d14-8a0e-91f83185d0ad
table_0 = fboys.(0,ts)

# ╔═╡ c615c015-34f6-441f-a285-0c7b8926660b
f0 = Spline1D(ts,fboys.(0,ts))

# ╔═╡ 85999318-1639-4677-bd9c-81723b4c809f
f0(0.5)

# ╔═╡ 61e753cb-f6c5-4301-8248-da8c21ba9186
fboys(0,0.5)

# ╔═╡ 24273947-c334-4508-a8e3-eab88100a745
function compute_error(steps)
	errs = []
	for step in steps
		mesh = 0:step:100
		f = Spline1D(mesh,fboys.(0,mesh))
		finemesh = 0:0.01:100
		error(x) = abs(fboys(0,x)-f(x))
		push!(errs,maximum(error(x) for x in finemesh))
	end
	return errs
end

# ╔═╡ 31c80276-35e5-4dc8-9a26-9600805a32da
steps = [0.05,0.1,0.2,0.5,1.0]

# ╔═╡ e3d4d6c7-c54b-408e-8575-331c7db3fec3
errs = compute_error(steps)

# ╔═╡ ca6f1d12-53c9-4df1-9194-aedb9480ec7c
md"Plot the dependence of the maximum error as a function of the 1d mesh spacing"

# ╔═╡ d73c0af0-559b-450a-b429-6d18991c839a
plot(steps,errs,yaxis=:log)

# ╔═╡ e0954971-c7a3-4847-a45a-ff22adc71490
md"Preliminary timing suggests that this is slow"

# ╔═╡ 88d466c8-1f3f-4c96-9507-18dc8cc123d1
@time f0(0.123)

# ╔═╡ 53659193-15dd-4bdc-820b-43422887b580
@time fboys(0,0.123)

# ╔═╡ 4161d94c-d545-4105-84e0-c92fd1344b05
md"## Implementing computation of a series of m-values for Fm(T)"

# ╔═╡ 44107e8f-116c-475f-b6d4-2b5aecc2f1e8
md"""
## References
[^boys50]: S.F. Boys. "Electronic Wave Functions. I. A General Method of Calculation for the Stationary States of Any Molecular System." *Proc. R. Soc. Lond.*, A200, 542-554 (1950).
[^gjp91]: P.M.W. Gill, B.G. Johnson, "Two-Electron Integrals Over Gaussian *s* Functions". J.A. Pople. *IJQC*. 40, 745-752, 1991.
[^libint]: E.F. Valeev. "Libint: A library for the evaluation of molecular integrals of many-body operators over Gaussian functions." [http://libint.valeyev.net/](http://libint.valeyev.net/)
[^mn57]: K. Nishimoto, N. Mataga. "Electronic Structure and Spectra of Some Nitrogen Heterocycles." *Z. Phys.*, 12, 335-338 (1957).
"""

# ╔═╡ Cell order:
# ╠═f80fd13d-207a-458a-b326-77256fe946e0
# ╠═bda0f8d4-5ab0-4c65-9d79-71ccbfe8dd3a
# ╟─4e16dbbe-b0fc-11eb-2349-ad7ccd0ff56b
# ╠═8a736133-17f2-446d-9ab5-4bc46bdd556b
# ╟─bf014265-210c-412f-9af4-5a654f557de9
# ╟─d8061039-367c-42ee-8a65-a32dca37fc67
# ╠═bc53663d-1f9f-43a6-8140-e8eeac1618bb
# ╟─1a0bffb2-8656-44db-b3b3-6cf6b2494b7c
# ╠═145a0007-b8d7-46a1-9cb8-ee063057527c
# ╠═b6b757ad-7fca-46c5-b42f-20f94ae4cbd7
# ╟─f6fa9869-b454-458d-95e3-818a9112cc52
# ╠═e7454913-e406-4b4a-856f-47103a7bf23d
# ╠═cdc2c6c9-e567-4ee6-a002-8fea660e35fa
# ╟─b5f42a47-eb6c-4b59-9827-a07b81618018
# ╟─105717a4-f481-479d-bfbb-dc35cc3c0997
# ╟─0b513de7-cca5-45ca-a4fa-2f2604eb7da9
# ╟─4dd0b1e6-61f0-4daa-80f1-625dc793b33a
# ╟─fdd2691f-f26f-4c26-8cc6-0d8edaa5278a
# ╠═a0a09c9f-c6d2-422f-b32a-8e3f8131b3c3
# ╠═ae9658a3-9241-423d-9daa-d6644e8a206a
# ╟─8c8adf9b-5eff-4942-9a49-95a4aa1e2002
# ╠═9428670f-5cf7-4c9e-bb93-dd0747a34182
# ╠═771d4eeb-72e6-4863-a114-ff203a42455b
# ╟─997e6018-bb0e-4933-8b13-b904e96fc819
# ╠═894135da-1ef9-4170-aa49-1b6307e398eb
# ╠═9dfb09bb-22d8-44f2-9c4b-98777e43de64
# ╠═d526caa1-ff7e-4ec5-ab59-2a10e9ad9c5e
# ╠═3923126f-93be-48f2-bd37-3ef50588ff47
# ╠═eb05a798-8eb2-4ff1-b88f-b58babe64891
# ╠═84df8774-a1f6-49f7-b84f-c1482b51dd9f
# ╠═08c5dc9b-0cb1-4bea-bd97-4872f46e63bf
# ╟─347ca2e4-b932-4451-a12e-15ca7e23d052
# ╠═b608f792-dd91-4bd2-97ab-baac7b7d2f9c
# ╠═b15df2d5-f407-4d2d-959f-18b21c8a3b60
# ╠═ae817c36-d9b8-4d14-8a0e-91f83185d0ad
# ╠═c615c015-34f6-441f-a285-0c7b8926660b
# ╠═85999318-1639-4677-bd9c-81723b4c809f
# ╠═61e753cb-f6c5-4301-8248-da8c21ba9186
# ╠═24273947-c334-4508-a8e3-eab88100a745
# ╠═31c80276-35e5-4dc8-9a26-9600805a32da
# ╠═e3d4d6c7-c54b-408e-8575-331c7db3fec3
# ╟─ca6f1d12-53c9-4df1-9194-aedb9480ec7c
# ╠═d73c0af0-559b-450a-b429-6d18991c839a
# ╟─e0954971-c7a3-4847-a45a-ff22adc71490
# ╠═88d466c8-1f3f-4c96-9507-18dc8cc123d1
# ╠═53659193-15dd-4bdc-820b-43422887b580
# ╠═4161d94c-d545-4105-84e0-c92fd1344b05
# ╠═44107e8f-116c-475f-b6d4-2b5aecc2f1e8

max_stars_repo_path: tools/simpleboys backup 1.jl

max_issues_repo_path: tools/simpleboys backup 1.jl

max_forks_repo_path: tools/simpleboys backup 1.jl


text:
include("boundary_mps.jl")
include("projectors.jl")
include("chain_rules.jl")


text:
########################################################################
# Atoms of the 1-norm ball.
########################################################################

"""
    OneBallAtom(n, z)

Atom in the OneBall of dimension n.
"""
mutable struct OneBallAtom{T1<:Int64, T2<:SparseVector{Float64, Int64}} <: AbstractAtom
    n::T1
    z::T2
    function OneBallAtom(z::SparseVector{Float64, Int64})
        n = z.n
        new{Int64, SparseVector{Float64, Int64}}(n, z)
    end
end

"""
Multiply a one-norm atom by a linear map.
"""
Base.:(*)(M::AbstractLinearOp, a::OneBallAtom{Int64, SparseVector{Float64, Int64}}) = M*a.z
LinearAlgebra.mul!(Ma::AbstractVector{Float64}, M::LinearOp, a::OneBallAtom{Int64, SparseVector{Float64, Int64}}) = mul!(Ma, M, a.z)
Base.vec(a::OneBallAtom{Int64, SparseVector{Float64, Int64}}) = a.z
Base.length(a::OneBallAtom{Int64, SparseVector{Float64, Int64}}) = a.n

########################################################################
# Atomic set for the 1-norm.
########################################################################

"""
    OneBall(n, maxrank=n)

Atomic set defined by 1-norm ball in `n` variables.
The atomic set takes two optional parameters:

`maxrank` is the maximum dimension of the face exposed by a vector.
"""
struct OneBall{T<:Int64} <: AbstractAtomicSet
    n::T
    maxrank::T
    function OneBall(n::Int64, maxrank::Int64)
        n ≥ 1 || throw(DomainError(n,"n must be a positive integer"))
        n ≥ maxrank ≥ 1 || throw(DomainError(maxrank,"maxrank must be ≥ 1"))
        new{Int64}(n, maxrank)
    end
end
OneBall(n::Int64; maxrank=n) = OneBall(n, maxrank)

"""
    gauge(A::OneBall, x::Vector)

Gives the 1-norm of `x`.
"""
gauge(::OneBall{Int64}, x::Vector{Float64}) = norm(x,1)
gauge(A::OneBall{Int64}, x::Matrix{Float64}) = gauge(A, vec(x))

"""
    support(A::OneBall, z::Vector)

Gives the inf-norm of `z`.
"""
support(::OneBall{Int64}, z::Vector{Float64}) = norm(z,Inf)

"""
    expose!(A::OneBall, z::Vector, a::OneBallAtom; tol=1e-2)

Obtain an atom in the face exposed by the vector `z`.
The atom `a` is overwritten.
"""
function expose!(A::OneBall{Int64}, z::Vector{Float64}, a::OneBallAtom{Int64, SparseVector{Float64, Int64}}; tol=1e-2)
    n = A.n
    @inbounds a.z .= spzeros(n)
    zmax = support(A, z)
    if zmax ≥ 1e-12
        nnz = 0
        @inbounds for i in 1:n
            val = z[i]
            if zmax - abs(val) < tol*zmax
                nnz += 1
                push!(a.z.nzind, i)
                push!(a.z.nzval, sign(val))
            end
        end
        a.z.nzval ./= nnz
    end
    return nothing
end

function expose(A::OneBall{Int64}, z::Vector{Float64})
    a = OneBallAtom(spzeros(A.n))
    expose!(A, z, a)
    return a
end

Base.length(A::OneBall{Int64}) = A.n
rank(A::OneBall{Int64}) = A.maxrank
atom_name(A::OneBall{Int64}) = "one-norm ball"
atom_description(A::OneBall{Int64}) = "{ x ∈ ℝⁿ | ||x||₁ ≤ 1 }"
atom_parameters(A::OneBall{Int64}) = "n = $(length(A)); rank = $(A.maxrank)"

########################################################################
# Face of the one-norm ball.
########################################################################

mutable struct OneBallFace{T1<:Int64, T2<:SparseMatrixCSC{Float64, Int64}} <: AbstractFace
    n::T1
    k::T1
    S::T2
    function OneBallFace(S::SparseMatrixCSC{Float64, Int64})
        n, k = size(S)
        new{Int64, SparseMatrixCSC{Float64, Int64}}(n, k, S)
    end
end

"""
    face!(A, z, F)

(Inplace) Return a face of the atomic set `A` exposed by the vector `z`. The dimension of
the exposed face is equal to `k=rank(A)`.
"""
function face!(A::OneBall{Int64}, z::Vector{Float64}, F::OneBallFace{Int64, SparseMatrixCSC{Float64, Int64}})
    k = A.maxrank
    idx = partialsortperm(abs.(z), 1:k, rev=true)
    val = [sign(z[i]) for i in idx]
    @inbounds F.S.rowval .= idx
    @inbounds F.S.nzval .= val
    return nothing
end

function face(A::OneBall{Int64}, z::Vector{Float64})
    n, k = A.n, A.maxrank
    F = OneBallFace( sparse(collect(1:k), collect(1:k), ones(k), n, k) )
    face!(A, z, F)
    return F
end

"""
Given a face and a set of weights, reveal the corresponding
point on the face.
"""
Base.:(*)(F::OneBallFace{Int64, SparseMatrixCSC{Float64, Int64}}, c::Vector{Float64}) = F.S*c
Base.:(*)(M::AbstractLinearOp, F::OneBallFace{Int64, SparseMatrixCSC{Float64, Int64}}) = M*F.S
Base.:(*)(λ::Float64, F::OneBallFace{Int64, SparseMatrixCSC{Float64, Int64}}) = λ*F.S
Base.length(F::OneBallFace{Int64, SparseMatrixCSC{Float64, Int64}}) = F.n
rank(F::OneBallFace{Int64, SparseMatrixCSC{Float64, Int64}}) = F.k
vec(F::OneBallFace{Int64, SparseMatrixCSC{Float64, Int64}}) = F.S
face_name(F::OneBallFace{Int64, SparseMatrixCSC{Float64, Int64}}) = "one-norm ball"
face_parameters(F::OneBallFace{Int64, SparseMatrixCSC{Float64, Int64}}) = "rank = $(rank(F)); n = $(length(F))"


text:
# input messages (to julia)
const MSG_INPUT_NULL              = 0
const MSG_INPUT_EVAL              = 1
const MSG_INPUT_REPLAY_HISTORY    = 2
const MSG_INPUT_GET_USER          = 3
const MSG_INPUT_AUTHENTICATE      = 4

# output messages (to the browser)
const MSG_OUTPUT_WELCOME          = 1
const MSG_OUTPUT_READY            = 2
const MSG_OUTPUT_MESSAGE          = 3
const MSG_OUTPUT_OTHER            = 4
const MSG_OUTPUT_EVAL_INPUT       = 5
const MSG_OUTPUT_FATAL_ERROR      = 6
const MSG_OUTPUT_EVAL_INCOMPLETE  = 7
const MSG_OUTPUT_EVAL_RESULT      = 8
const MSG_OUTPUT_EVAL_ERROR       = 9
const MSG_OUTPUT_PLOT             = 10
const MSG_OUTPUT_GET_USER         = 11
const MSG_OUTPUT_HTML             = 12


text:
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute(path::String)::String
    keylog = ∪(split(read(path, String), "\n")[1:end - 1])
    password = ""
    while length(keylog) ≠ 1
        candidates = ∪(i[1] for i ∈ keylog)
        for candidate ∈ collect(candidates)
            for key ∈ keylog
                if candidate ∈ key[2:end]
                    setdiff!(candidates, candidate)
                    break
                end
            end
        end
        password = string(password, candidates[1])
        for key ∈ collect(keylog)
            if candidates[1] == key[1]
                if length(key) == 3 union!(keylog, [key[2:end]]) end
                setdiff!(keylog, [key])
            end
        end
    end
    return string(password, keylog[1])
end

compute("problems/0079/p079_keylog.txt")

@benchmark compute("problems/0079/p079_keylog.txt")

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
export @plasticity_str
macro plasticity_str(plasticity)
  str="function plasticity!(net::Network$timestr)\n"
  plasticity=replace(plasticity, r"(?<=^)(.+)(?=\[.*])"m, s"net.\1")
  plasticity=replace(plasticity, r"(?<==)(.+)(?=\[.*])"m, s"net.\1")
  plasticity=replace(plasticity,"j",":")
  plasticity=replace(plasticity,"i","net.i")
  str*=plasticity*"\n"
  str*="end\n"
  print(str)
  eval(parse(str))
end
# general case
# wmax = 0.01
# plasticity"""
# apre[j,i] += Apre[j,i];
# w[j,i] = clip(w[j,i]+apost[j,i],0,wmax)
# apost[i,j]+=Apost[i,j];
# w[i,j]=clip(w[i,j]+apost[i,j],0,wmax)"""


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

using Markdown
using InteractiveUtils

# ╔═╡ 0c45e9aa-e0ba-11eb-0f56-af8192849f70
begin
  using Pkg
  Pkg.add([
    "PlutoUI",
    "Plots",
    "PlotlyJS",
    "WebIO",
  ])
  using PlutoUI
  using Plots
  #plotly()
  gr()
  using PlotlyJS
  using WebIO
end

# ╔═╡ 87226ee7-7adc-4fab-a86a-c242a85f14d2
let
  u = [1. 0.]'
  x = reshape(range(-1, 1; length=11), 1, :)
  u * x
end

# ╔═╡ 86e93269-6f70-4ebb-8390-74c4f8111cce
let
  u = [1. 0.]'
  v = [0. 1.]'
  Tu = [3. 4.]'
  Tv = [-1. 0.]'
  n_pts = 100
  #x = reshape(range(-1, 1; length=n_pts÷2), 1, :)
  x = range(-1, 1; length=n_pts÷2)
  y = .√(1 .- x.^2)
  y = [y; -y[2:end-1]]
  x = [x; x[2:end-1]]
  Plots.scatter(x, y,
    aspect_ratio=:equal,
    background_color=:black
  )
end

# ╔═╡ c7f38ef0-263a-45a7-b9ac-8b9c165419b2
function T(coeff, basis)
  return basis * coeff
end

# ╔═╡ ade94a6d-a68d-4c82-9dba-89fd70fd025e
let
  coeff = [1 0 -1  0
           0 1  0 -1]
  Tu = [3. 4.]'
  Tv = [-1. 0.]'
  basis = [Tu Tv]
  T(coeff, basis)
end

# ╔═╡ 7e348501-80aa-43a0-8b82-4cd0d8081df1
let
  n_pts = 100
  x = reshape(range(-1, 1; length=n_pts÷2), 1, :)
  y = .√(1 .- x.^2)
  size(y), size(-y[2:end-1]), size(-y[1:end, 2:end-1])
end

# ╔═╡ 814749da-b5de-4247-9eba-58b2055dd4fd
let
  n_pts = 100
  x = reshape(range(-1, 1; length=n_pts÷2), 1, :)
  y = .√(1 .- x.^2)
  hcat(y, -y[1:end, 2:end-1])
end

# ╔═╡ 6e4da27d-4861-4d3a-a3af-0c966685d698
md"""
The following recreates the same circle as above but with more flexibility.
"""

# ╔═╡ ac39c912-7480-461f-97d6-376ce1f01cc1
let
  u = [1. 0.]'
  v = [0. 1.]'
  n_pts = 100
  x = reshape(range(-1, 1; length=n_pts÷2), 1, :)
  y = .√(1 .- x.^2)
  #x = [x, x[2:end-1]]
  coeff = vcat(
    hcat(x, x[1:end, 2:end-1]),
    hcat(y, -y[1:end, 2:end-1]),
  )
  pts = [u v] * coeff
  Plots.scatter(
    pts[1, 1:end],
    pts[2, 1:end],
    aspect_ratio=:equal,
    background_color=:black
  )
end

# ╔═╡ f943106b-1d37-4240-95f4-087b958d62e3
let
  Tu = [3. 4.]'
  Tv = [-1. 0.]'
  basis = [Tu Tv]
  n_pts = 100
  x = reshape(range(-1, 1; length=n_pts÷2), 1, :)
  y = .√(1 .- x.^2)
  coeff = vcat(
    hcat(x, x[1:end, 2:end-1]),
    hcat(y, -y[1:end, 2:end-1]),
  )
  pts = T(coeff, basis)
  Plots.scatter(
    pts[1, 1:end],
    pts[2, 1:end],
    aspect_ratio=:equal,
    background_color=:black
  )
end

# ╔═╡ 56348137-1190-4a96-9c2b-260b5c5b574a
size([1 2 3; 4 5 6], 2)

# ╔═╡ cb58d7fd-0aa1-4371-bcc1-89dd476e4cb1
[1 2 3; 4 5 6][1:end, [2]]

# ╔═╡ b8a8190b-2ffd-4ba0-866c-ed12477b672e
let
  Tu = [3. 4.]'
  Tv = [-1. 0.]'
  basis = [Tu Tv]
  n_pts = 100
  x = reshape(range(-1, 1; length=n_pts÷2), 1, :)
  y = .√(1 .- x.^2)
  coeff = vcat(
    hcat(x, x[1:end, 2:end-1]),
    hcat(y, -y[1:end, 2:end-1]),
  )
  plt = Plots.plot(
    1,
    xlim = (-30, 30),
    ylim = (-30, 30),
    aspect_ratio=:equal,
    background_color=:black,
    title = "Image of unit sphere in 2D",
    marker = 2,
  )
  # build an animated gif by pushing new points to the plot, saving every 3 frame
  @gif for col = 1:size(coeff, 2)
    pt = T(coeff[1:end, [col]], basis)
    push!(plt, pt[1, 1], pt[2, 1])
  end every 2
end

# ╔═╡ ccdca454-c175-4f67-bd4a-8663e4bd9756
md"""
Quite neat now, but it would become better if we make the motion
- direction-consistent
- zoomed-in
- no label
"""

# ╔═╡ c0ecedae-c66e-46a3-989e-7fc9642577fc
let
  Tu = [3. 4.]'
  Tv = [-1. 0.]'
  basis = [Tu Tv]
  n_pts = 100
  x = reshape(range(-1, 1; length=n_pts÷2), 1, :)
  y = .√(1 .- x.^2)
  coeff = vcat(
    hcat(x, x[1:end, end-1:-1:2]),
    hcat(y, -y[1:end, end-1:-1:2]),
  )
  plt = Plots.plot(
    1,
    xlim=(-7, 7),
    ylim=(-7, 7),
    aspect_ratio=:equal,
    background_color=:black,
    title="Image of unit sphere in 2D",
    marker=2,
    label=false,
  )
  # build an animated gif by pushing new points to the plot, saving every 3 frame
  @gif for col = 1:size(coeff, 2)
    pt = T(coeff[1:end, [col]], basis)
    push!(plt, pt[1, 1], pt[2, 1])
  end every 1
end

# ╔═╡ 2e27720a-1ea1-48ba-85fd-ba7c12355a9d
let
  Tu = [3. 4.]'
  Tv = [0. 0.]'
  basis = [Tu Tv]
  n_pts = 100
  x = reshape(range(-1, 1; length=n_pts÷2), 1, :)
  y = .√(1 .- x.^2)
  coeff = vcat(
    hcat(x, x[1:end, end-1:-1:2]),
    hcat(y, -y[1:end, end-1:-1:2]),
  )
  plt = Plots.plot(
    1,
    xlim=(-7, 7),
    ylim=(-7, 7),
    aspect_ratio=:equal,
    background_color=:black,
    title="Image of unit sphere in 2D",
    marker=2,
    label=false,
  )
  # build an animated gif by pushing new points to the plot, saving every 3 frame
  @gif for col = 1:size(coeff, 2)
    pt = T(coeff[1:end, [col]], basis)
    push!(plt, pt[1, 1], pt[2, 1])
  end every 1
end

# ╔═╡ ee76cf5b-ada9-416d-be94-5451eafce775
md"""
Let's try with 3D.

It's a little hard to do the GIF with 3D, but let's at least start with a scatter plot.
"""

# ╔═╡ 7d0d395b-8c9f-4bd6-bc2d-5ff9af0e326d
[x for x in rand(range(-1,1;length=1_000), (2, 100))]

# ╔═╡ f9729911-9fb6-499e-b1e3-4e60775a8f4f
let
  square_xy = rand(range(-1,1;length=1_000), (2, 100))
  circle_xy = [square_xy[1:end, col] for col in 1:size(square_xy, 2) if sum(square_xy[1:end, col].^2) <= 1]
  #size(circle_xy)
  #square_xy, [col for col in 1:size(square_xy, 2) if sum(square_xy[1:end, col].^2) <= 1]
  #square_xy, [sum(square_xy[1:end, col].^2) for col in 1:size(square_xy, 2)]
end

# ╔═╡ b73ed7bc-4d9c-4f83-813e-9cf76ec602a8
md"""
### Stupidity: Polar Coordinate
I was so stupid as to forget the useful **_polar coordinate_**. Let's re-code the above using polar coordinate! I believe it's going to be simpler.
"""

# ╔═╡ e7640a11-24c0-49e8-b31e-58e8b2cff272
float(rand(-10:10, (2,2)))

# ╔═╡ d82027bb-dfe0-4fa7-8406-5212b62ab757
let
  Tu = [3. 4.]'
  Tv = [1. 0.]'
  basis = [Tu Tv]
  n_angles = 100
  ϵ = 1e-9
  #θ = range(0, 2π-ϵ; length=n_angles)
  coeff = zeros((2, n_angles))
  for (col, θ) in enumerate(range(0, 2π-ϵ; length=n_angles))
    coeff[1:end, col] = [cos(θ), sin(θ)]
  end
  pts = T(coeff, basis)
  plt = Plots.scatter(
    1,
    xlim=(-7, 7),
    ylim=(-7, 7),
    aspect_ratio=:equal,
    background_color=:black,
    title="Image of unit sphere in 2D",
    marker=2,
    label=false,
  )
  @gif for col = 1:size(coeff, 2)
    pt = pts[1:end, col]
    push!(plt, pt[1], pt[2])
  end every 1
end

# ╔═╡ 212c258a-3542-48ea-8259-ce9b282e41c6
md"""
#### ``T: \mathbb{R}^3 \to \mathbb{R}^2``
"""

# ╔═╡ 50a2def3-affa-4396-9041-a93131ffc7bd
let
  Tu = [3. 4.]'
  Tv = [-1. 0.]'
  Tw = [1. -2.]'
  basis = [Tu Tv Tw]
  n_angles = 50
  ϵ = 1e-9
  coeff = zeros((3, n_angles^2))
  counter = 1
  for ϕ in range(0, π-ϵ; length=n_angles), θ in range(0, 2π-ϵ; length=n_angles)
    coeff[1:end, counter] = [sin(ϕ)*cos(θ), sin(ϕ)*sin(θ), cos(ϕ)]
    counter += 1
  end
  pts = T(coeff, basis)
  size(pts)
end

# ╔═╡ b72db565-9c93-442f-8284-852c01af6f7b
let
  Tu = [3. 4.]'
  Tv = [-1. 0.]'
  Tw = [1. -2.]'
  basis = [Tu Tv Tw]
  n_angles = 50
  ϵ = 1e-9
  coeff = zeros((3, n_angles^2))
  counter = 1
  for ϕ in range(0, π-ϵ; length=n_angles), θ in range(0, 2π-ϵ; length=n_angles)
    coeff[1:end, counter] = [sin(ϕ)*cos(θ), sin(ϕ)*sin(θ), cos(ϕ)]
    counter += 1
  end
  pts = T(coeff, basis)
  plt = Plots.scatter(
    1,
    xlim=(-7, 7),
    ylim=(-7, 7),
    aspect_ratio=:equal,
    background_color=:black,
    title="Image of unit sphere in 2D",
    #marker=2,
    label=false,
  )
  @gif for col = 1:size(coeff, 2)
    pt = pts[1:end, col]
    push!(plt, pt[1], pt[2])
  end every 1
end

# ╔═╡ 07539502-d5ef-48a0-9f2f-ae740f0a6c01
let
  #ρ = 1
  n_angles = 100
  ϵ = 1e-9
  coeff = zeros((3, n_angles^2))
  counter = 1
  #for (i, ϕ) in enumerate(range(0, π-ϵ; length=n_angles)), (j, θ) in enumerate(range(0, 2π-ϵ; length=n_angles))
  for ϕ in range(0, π-ϵ; length=n_angles), θ in range(0, 2π-ϵ; length=n_angles)
    coeff[1:end, counter] = [sin(ϕ)*cos(θ), sin(ϕ)*sin(θ), cos(ϕ)]
    counter += 1
  end
  basis = float(rand(-10:10, (3,3)))
  pts = T(coeff, basis)
  Plots.scatter(
    pts[1, 1:end],
    pts[2, 1:end],
    pts[3, 1:end],
    aspect_ratio=:equal,
    background_color=:black,
    label=false,
  )
end

# ╔═╡ f7b1407b-3c9a-4a66-ae2c-259fc6a2e1aa
md"""
Plotting like this is not very helpful.

This link <https://plotly.com/julia/3d-scatter-plots/> provides better scatter plots.

> It seems that extra effort is needed in order to use `PlotlyJS`, so I disabled it.
"""

# ╔═╡ a31dff2c-766c-4d52-8fb3-f9d058f310c3
PlotlyJS.plot, PlotlyJS.scatter

# ╔═╡ 0cedae98-fe16-47c9-a735-7c86af0c803a
let
  #ρ = 1
  n_angles = 100
  ϵ = 1e-9
  coeff = zeros((3, n_angles^2))
  counter = 1
  #for (i, ϕ) in enumerate(range(0, π-ϵ; length=n_angles)), (j, θ) in enumerate(range(0, 2π-ϵ; length=n_angles))
  for ϕ in range(0, π-ϵ; length=n_angles), θ in range(0, 2π-ϵ; length=n_angles)
    coeff[1:end, counter] = [sin(ϕ)*cos(θ), sin(ϕ)*sin(θ), cos(ϕ)]
    counter += 1
  end
  basis = float(rand(-10:10, (3,3)))
  pts = T(coeff, basis)
  #plotlyjs()
  PlotlyJS.Plot(PlotlyJS.scatter3d(
    x=pts[1, 1:end],
    y=pts[2, 1:end],
    z=pts[3, 1:end],
    height=400,
    #width=900,
  ))
end

# ╔═╡ 259744e6-4b87-43fb-ae9e-88e02c2c0f35
md"""
`PlotlyJS`'s 3D plot in this case is too small. Let's plot using `@gif` again.
"""

# ╔═╡ 03774574-6ebe-4146-875d-f654a88457c9
let
  #ρ = 1
  n_angles = 50
  ϵ = 1e-9
  coeff = zeros((3, n_angles^2))
  counter = 1
  #for (i, ϕ) in enumerate(range(0, π-ϵ; length=n_angles)), (j, θ) in enumerate(range(0, 2π-ϵ; length=n_angles))
  for ϕ in range(0, π-ϵ; length=n_angles), θ in range(0, 2π-ϵ; length=n_angles)
    coeff[1:end, counter] = [sin(ϕ)*cos(θ), sin(ϕ)*sin(θ), cos(ϕ)]
    counter += 1
  end
  basis = float(rand(-10:10, (3,3)))
  pts = T(coeff, basis)
  sct = Plots.plot3d(
    1,
    xlim=(-15, 15),
    ylim=(-20, 20),
    zlim=(-7, 7),
    aspect_ratio=:equal,
    background_color=:black,
    label=false,
  )
  @gif for col in 1:size(pts, 2)
    push!(sct, pts[1, col], pts[2, col], pts[3, col])
  end every 5
end

# ╔═╡ 6c790f78-464d-49c2-9621-d3677cd70d9e

# ╔═╡ Cell order:
# ╠═0c45e9aa-e0ba-11eb-0f56-af8192849f70
# ╠═87226ee7-7adc-4fab-a86a-c242a85f14d2
# ╠═86e93269-6f70-4ebb-8390-74c4f8111cce
# ╠═c7f38ef0-263a-45a7-b9ac-8b9c165419b2
# ╠═ade94a6d-a68d-4c82-9dba-89fd70fd025e
# ╠═7e348501-80aa-43a0-8b82-4cd0d8081df1
# ╠═814749da-b5de-4247-9eba-58b2055dd4fd
# ╟─6e4da27d-4861-4d3a-a3af-0c966685d698
# ╠═ac39c912-7480-461f-97d6-376ce1f01cc1
# ╠═f943106b-1d37-4240-95f4-087b958d62e3
# ╠═56348137-1190-4a96-9c2b-260b5c5b574a
# ╠═cb58d7fd-0aa1-4371-bcc1-89dd476e4cb1
# ╠═b8a8190b-2ffd-4ba0-866c-ed12477b672e
# ╟─ccdca454-c175-4f67-bd4a-8663e4bd9756
# ╠═c0ecedae-c66e-46a3-989e-7fc9642577fc
# ╠═2e27720a-1ea1-48ba-85fd-ba7c12355a9d
# ╟─ee76cf5b-ada9-416d-be94-5451eafce775
# ╠═7d0d395b-8c9f-4bd6-bc2d-5ff9af0e326d
# ╠═f9729911-9fb6-499e-b1e3-4e60775a8f4f
# ╟─b73ed7bc-4d9c-4f83-813e-9cf76ec602a8
# ╠═e7640a11-24c0-49e8-b31e-58e8b2cff272
# ╠═d82027bb-dfe0-4fa7-8406-5212b62ab757
# ╟─212c258a-3542-48ea-8259-ce9b282e41c6
# ╠═50a2def3-affa-4396-9041-a93131ffc7bd
# ╠═b72db565-9c93-442f-8284-852c01af6f7b
# ╠═07539502-d5ef-48a0-9f2f-ae740f0a6c01
# ╟─f7b1407b-3c9a-4a66-ae2c-259fc6a2e1aa
# ╠═a31dff2c-766c-4d52-8fb3-f9d058f310c3
# ╠═0cedae98-fe16-47c9-a735-7c86af0c803a
# ╟─259744e6-4b87-43fb-ae9e-88e02c2c0f35
# ╠═03774574-6ebe-4146-875d-f654a88457c9
# ╠═6c790f78-464d-49c2-9621-d3677cd70d9e