text:
tic()
@testset "Invariant distribution tests" begin
   dist = Normal()
   @testset "Time series length: $ts_length" for ts_length in 50:50:300
      @testset "Rep #$rep" for rep in 1:3

         is_invariant = false
         TS = Float64[]
         embedding = Array{Float64}(0, 0)
         while (!is_invariant)
            ts = rand(dist, ts_length)
            e = embed(ts, 3, 2)
            if invariantset(e)
               is_invariant = true
               TS = ts
               embedding = e
            end
         end

         println("length(ts) = ", ts_length, "\trep=", rep, "\tThe embedding forms an invariant set")
         isinvariant = true

         # Embed using all points except the last (to allow projection of all vertices
         # in the triangulation).
         points, simplex_inds = triangulate(embedding[1:end-1, :])
         image_points = embedding[2:end, :]

         @testset "Simplex vertices matches embedding[1:end-1, :]" begin
            @test points == embedding[1:end-1, :]
         end

         P = markovmatrix(points, image_points, simplex_inds)

         @testset "Markov matrix" begin
            @test all(sum(P, 2) .≈ 1.0)
         end

         invmeasure, inds_nonzero_simplices = invariantdist(P)

         @testset "Invariant measure" begin
            # There must be at least one simplex with nonzero measure; otherwise, the invariant
            # distribution cannot sum to 1.
            @test size(inds_nonzero_simplices, 1) >= 1

            # The invariant distribution must sum to 1 within reasonable tolerance.
            @test sum(invmeasure) ≈ 1.0 atol = 1/10^6
         end
            #
            # nonzero_simplices = simplex_inds[inds_nonzero_simplices, :]
            # vertices = points
            # centroids, radii = SimplexSplitting.centroids_radii2(vertices, nonzero_simplices)
            # #@show size(centroids), size(invmeasure)
            # n = 3
            # nonempty_bins, measure = get_nonempty_bins(centroids, invmeasure, [n, n, n])
            #
            # unique_nonempty_bins = unique(nonempty_bins, 1)
            #
            # @testset "Bin measure is 1" begin
            #    @test sum(measure) ≈ 1.0 atol = 1/10^6
            # end
            #
            # @testset "All centroids are assigned to bins" begin
            #    @test !any(nonempty_bins .== 0)
            # end
            #
            #
            # @testset "Joint distribution" begin
            #    joint = jointdist(nonempty_bins, measure)
            #    @show joint
            # end
            #
            # @testset "Marginal distribution" begin
            #    Py, Pxy, Pxz = marginaldists(unique_nonempty_bins, measure)
            #    @show Py
            #    @show Pxy
            #    @show Pxz
            # end
            #
            # joint = jointdist(nonempty_bins, measure)
            # Py, Pxy, Pxz = marginaldists(unique_nonempty_bins, measure)
            #
            # #@show size(joint), size(Py), size(Pxy), size(Pxz)
            #
            # @testset "Joint and marginals are the same size" begin
            #    sizes = [size(distribution, 1) for distribution in (joint, Py, Pxy, Pxz)]
            #    @test length(unique(sizes)) == 1
            # end
      end
   end
end
toc()


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ def88812-6b99-11eb-07fb-3f961d7ff6ae
begin
	using Plots
	plotly()
end

# ╔═╡ 67bd1b94-7f32-11eb-2869-5f075fca7fbc
begin
	using PlutoUI
	using StatsPlots
	using RDatasets
	using StatsBase
	using Distributions
	using LaTeXStrings
end

# ╔═╡ 253ad4f4-8253-11eb-03d8-d358d6075910
md"A mid-sem sample question:

Write a Julia snippet to sample 100 random variates denoted $\{x_1, x_2, \ldots, x_{100}\}$ from the standard normal distribution. Then compute $\dfrac{\sum_{i = 1}^{100} \sum_{j = 1}^{100} (x_i - x_j)^2}{100 \times 100}$."

# ╔═╡ 237d9290-8254-11eb-233a-4d362783d684
begin
	sample_ = rand(Normal(), 100)
	val_ = sum(sum([[(sample_[i] - sample_[j])^2 for i in 1:100] for j in 1:100]))/(100 * 100)
end

# ╔═╡ 8ce2cb16-8253-11eb-21a5-c59e5785950e
md"Run the above a few dozen times and build an intuition for what the value is for different samples in relation to $\sigma^2$ of the underlying distribution. Then deriva analytically $\mathbb{E}[(X - Y)^2]$, where $X$ and $Y$ are two independent random variables sampled from the same distribution with pdf $f$."

# ╔═╡ 22f18e82-6b7d-11eb-16eb-31758abdb1b0
md"# Notebook 6: _Measures of Two Variables_
"

# ╔═╡ 59607216-7b43-11eb-2cf7-3f74231a1839
md"
> 📕Reference material - Chapter 3: Using statistics to summarise data (Section 3.7) from Sheldon Ross.
"

# ╔═╡ 5f510e34-7f32-11eb-158a-afdded4ae4f7
iris = dataset("datasets", "iris");

# ╔═╡ 65c3f4a2-7f32-11eb-0eb4-97a2455d1ac7
diamonds = dataset("ggplot2", "diamonds");

# ╔═╡ b1716986-7e95-11eb-383d-39edf8743454
md"#### Statistics of two variables - Covariance"

# ╔═╡ 67e698c8-7e99-11eb-01c9-eb2703970ea4
md"So far we have been seeing the statistics of one variable. When there are two variables, we are interested in how they are related."

# ╔═╡ c6dc30e0-7e99-11eb-1781-7dba5912b724
md"Basic intuition is to check if increase in one of the variables also leads to an increase or decrase in another. _Do the two variables move together? Do they co-vary?_"

# ╔═╡ b79f4c20-7e99-11eb-3eba-492dc528735d
md"👉 Can you formulate a statistic that may help identify how two variables are related?"

# ╔═╡ 7b13d1a8-7e9a-11eb-39f2-5d3e9a95c811
iris_sample = select!(iris[sample(axes(iris, 1), 5; replace = false, ordered = true), :], Not(:Species))

# ╔═╡ a136b75a-7f04-11eb-1324-c9ea26353914
md"By just working with the petal length and width values, we cannot conclude if they are covarying. Notice how they are both of different scales as well."

# ╔═╡ ff01b9c2-7f02-11eb-2edd-5131a142f119
[mean(iris_sample.PetalLength) mean(iris_sample.PetalWidth)]

# ╔═╡ bbe2eec0-7f04-11eb-3f7f-594d57f134ec
md"Now that we know the mean values of the two columns are different, we can use it to compute offsets for individual values."

# ╔═╡ fc068802-7f03-11eb-36a0-f3764ac250ab
iris_offset = select(iris_sample, :PetalWidth, :PetalLength, :PetalWidth => (x -> x .- mean(x)) => :PW_offset, :PetalLength => (x -> x .- mean(x)) => :PL_offset)

# ╔═╡ eb742c80-7f04-11eb-0435-4f58142bcdc9
md"These offset values are now clearer: If two variables are to co-vary they must do so for these offset values."

# ╔═╡ 2f64101e-7f04-11eb-148b-d1c781173d4c
iris_covar = select(iris_offset, :PetalWidth, :PetalLength, :PW_offset, :PL_offset, [:PW_offset, :PL_offset] => ((x, y) -> x .* y) => :Covariance)

# ╔═╡ 98fcccfc-7f04-11eb-0b2d-6b1feaadfae4
md"All the values in the last column being positive shows that the two values actually co-vary (either both negative or both positive). A natural next step is to simply compute the average of these products."

# ╔═╡ 12e7d1a4-7f05-11eb-19d5-2b072cbafa05
mean(iris_covar.Covariance)

# ╔═╡ 2c8776c8-7f05-11eb-0226-bf1ef7542efd
md"Covariance of two sets of population values $\{X_1, X_2, \ldots, X_N\}$ and $\{Y_1, Y_2, \ldots, Y_N\}$ is given as
```math
\mathrm{COV}(X, Y) = \dfrac{\sum_{i = 1}^N (X_i - \mu_X)(Y_i - \mu_Y)}{N}
```
"

# ╔═╡ 9da0add4-7f05-11eb-1391-0f87db7d41c4
md"We can define this in terms of expectation as
```math
\mathrm{COV}(X, Y) = \mathbb{E}\left[(X - \mathbb{E}[X])(Y - \mathbb{E}[Y])\right]
```
"

# ╔═╡ a86b8dc4-7f05-11eb-3ca7-e5460fc15b96
md"Notice the parallels to the definition of $\sigma^2$ as $\mathbb{E}[(X - \mathbb{E}[X])^2]$"

# ╔═╡ 2e4af178-7f06-11eb-1247-3501f18174cc
md"👉Class exercise: Just like we rewrote $\sigma^2$ as $\mathbb{E}[X^2] - \mathbb{E}[X]^2$, can you rewrite $\mathrm{COV}(X, Y)$?"

# ╔═╡ e79a0a24-7f06-11eb-3c63-91d44fe78d49
md"Thus, we have $\mathrm{COV}(X, Y) = \mathbb{E}[XY] - \mathbb{E}[X]\mathbb{E}[Y]$"

# ╔═╡ 2be9abbe-824d-11eb-3b1d-2b1ddf16c195
md"👉 What is the covariance of two indpendent random variables $X$ and $Y$?"

# ╔═╡ 483778b4-824d-11eb-196c-2925d31dd603
md"👉 What is $\mathrm{COV}(aX, Y)$?"

# ╔═╡ b818699a-7f37-11eb-088c-a51bed3cbdc0
md"Here is a collection of some basic rules of covariance:

a) $\mathrm{COV}(X, X) = \sigma^2(X)$

b) If $X$ and $Y$ are independent then $\mathrm{COV}(X, Y) = 0$

c) $\mathrm{COV}(X, Y) = \mathrm{COV}(Y, X)$

d) $\mathrm{COV}(aX, Y) = a\mathrm{COV}(X, Y)$

e) $\mathrm{COV}(X + c, Y) = \mathrm{COV}(X, Y)$
"

# ╔═╡ 41664208-7f38-11eb-01f5-c32a33253f15
md"👉Homework exercise to prove these. All of these can be directly provded from the relation $\mathrm{COV}(X, Y) = \mathbb{E}[XY] - \mathbb{E}[X]\mathbb{E}[Y]$"

# ╔═╡ 6dfb60e6-7f38-11eb-070d-b99afff5f5a1
md"👉Class exercise: Show that $\mathrm{COV}(X + Y, Z) = \mathrm{COV}(X, Z) + \mathrm{COV}(Y, Z)$"

# ╔═╡ 2e32a1e2-7f39-11eb-0df0-812f689e1d1f
md"This rule can be generalised to the following
```math
\mathrm{COV}\left(\sum_i^m a_i X_i, \sum_j^n b_j Y_j\right) = \sum_i^m \sum_j^n a_i b_j \mathrm{COV}(X_i, Y_j)
```
"

# ╔═╡ c4d2c9d0-7f39-11eb-2ac9-373ed723a3f4
md"A neat application of the above is when we compute the variance of the sum of two random variables."

# ╔═╡ ce6afe4a-7f39-11eb-19f1-7157457700cd
md"
```math
\sigma^2(X + Y) = \mathrm{COV}(X+Y, X+Y)
```
```math
= \mathrm{COV}(X, X) + \mathrm{COV}(X, Y) + \mathrm{COV}(X, Y) + \mathrm{COV}(Y, Y)
```
```math
= \sigma^2(X) + 2\mathrm{COV}(X, Y) + \sigma^2(Y)
```
"

# ╔═╡ 175b1608-7f3a-11eb-0dfb-dd2ee73b9bac
md"Thus, we can interpret the covariance as half the value by which the variance of the sum increases from the sum of the variances for any two random variables."

# ╔═╡ fe62b640-8251-11eb-25fe-2b8aa69000ce
md"So we have seen three ways of computing covariance"

# ╔═╡ 27d06de0-7f07-11eb-23ce-4d6547ce31f6
mean(iris_covar.Covariance)

# ╔═╡ 0359bc70-7f07-11eb-2daa-d1f2a258d9e0
mean(iris_covar.PetalWidth .* iris_covar.PetalLength) - mean(iris_covar.PetalWidth) * mean(iris_covar.PetalLength)

# ╔═╡ a8c20bc6-8251-11eb-0d06-819a786e2b05
(var(iris_covar.PetalWidth .+ iris_covar.PetalLength, corrected=false) - var(iris_covar.PetalWidth, corrected=false) - var(iris_covar.PetalLength, corrected=false))/2

# ╔═╡ 9c42421a-7f0d-11eb-0172-f34b4adfa5f0
md"Notice that the covariance (like variance) has units that are product of the units of the individual variables. In the case of the iris dataset, the covariance we computed has the unit of length-squared."

# ╔═╡ 502bc23a-7f07-11eb-0b40-97e7e2a8a9ca
md"We can also directly use the `cov` function from `StatsBase`"

# ╔═╡ 2938c2ac-7f07-11eb-279c-ab02e04ac679
cov(iris_covar.PetalWidth, iris_covar.PetalLength, corrected=false)

# ╔═╡ 59828c88-7f07-11eb-1462-73e3850659fd
md"Notice we have used the `corrected = false` argument. This is because we are treating the dataset as a population. But in reality it is a sample and hence we need our correction."

# ╔═╡ a48951d0-7f07-11eb-0052-7d06cf26028c
md"Covariance of two sets of sample values $\{x_1, x_2, \ldots, x_n\}$ and $\{y_1, y_2, \ldots, y_n\}$ is given as
```math
\mathrm{COV}(x, y) = \dfrac{\sum_{i = 1}^n (x_i - \overline{x})(y_i - \overline{y})}{n - 1}
```
"

# ╔═╡ e0d54d10-7f07-11eb-214a-bd24b8accb1a
md"With this correction, the covariance is larger"

# ╔═╡ ec51f080-7f07-11eb-2559-5b4cf71a2d0e
cov(iris_covar.PetalWidth, iris_covar.PetalLength)

# ╔═╡ ef48b742-7f07-11eb-18e0-f1bcc96c86e0
md"We can then compute it for the entire dataframe instead of the small sample we were studying"

# ╔═╡ dc4d1258-7e99-11eb-25e2-ffcd05171460
cov(iris.PetalLength, iris.PetalWidth)

# ╔═╡ fa8d1044-7f07-11eb-247e-a18ecff0f4f3
md"Notice again that the large positive value of the covariance indicates that the two variable covary, This is easy to see in the following scatter plot."

# ╔═╡ 182c2ca2-7e9a-11eb-040a-21429d245b4f
scatter(iris.PetalLength, iris.PetalWidth)

# ╔═╡ 2348b81e-7f08-11eb-3070-25becdb70c5a
md"The shapes of the scatter plots reveals the covariance values."

# ╔═╡ 85e8d916-7f08-11eb-094c-17dbb38c5048
begin
	x1_ = rand(Normal(5, 2), 1000);
	y1_ = x1_ + rand(Normal(5, 1), 1000);
end

# ╔═╡ ed8c0606-7f08-11eb-13ec-d166c29a4d67
scatter(x1_, y1_, title="Covariance = " * string(cov(x1_, y1_)))

# ╔═╡ 175b39a2-7f09-11eb-02d7-c91ba9204c3c
begin
	x2_ = rand(Normal(5, 2), 1000);
	y2_ = rand(Normal(5, 2), 1000);
end

# ╔═╡ 1d4a9088-7f09-11eb-3754-07c4947ecad7
scatter(x2_, y2_, title="Covariance = " * string(cov(x2_, y2_)))

# ╔═╡ 27426ea6-7f09-11eb-3c2f-c9d5d5ed33b8
begin
	x3_ = rand(Normal(5, 2), 1000);
	y3_ = -x3_ + rand(Normal(5, 1), 1000);
end

# ╔═╡ 307bf5dc-7f09-11eb-3e23-8908558e17ca
scatter(x3_, y3_, title="Covariance = " * string(cov(x3_, y3_)))

# ╔═╡ 5a6fb146-7f09-11eb-02b3-0b2bc301eb95
md"So absolute value of covariance denotes the strength of the relationship while the sign denotes the direction (aligned is positive and oppositely aligned is negative). Two variables with 0 covariance are unrelated to each other."

# ╔═╡ 882610e4-7f09-11eb-142b-9d8ed5e88fcc
md"In real world data, we need to be careful about mixing data from different groups. For instance consider the covariance between petal length and sepal width acrsoss all flowers."

# ╔═╡ 0019040c-7e9a-11eb-3ab3-698eb0f98231
cov(iris.PetalLength, iris.SepalWidth)

# ╔═╡ 241c6342-7e9a-11eb-2053-c1984c19a023
scatter(iris.PetalLength, iris.SepalWidth)

# ╔═╡ 54ccb922-7f0a-11eb-3ba6-230b0489d928
md"The covariance value for all data together is negative. But the graph does show a positive slope."

# ╔═╡ 6273fb94-7f0a-11eb-3c1a-4dfa9d925400
md"We can gropuby the species column and compute covariance for each group separately"

# ╔═╡ b9d98800-7f09-11eb-3671-ed27dbde1c3a
combine(groupby(iris, :Species), [:PetalLength, :SepalLength] => ((x, y) -> cov(x, y))  => :Covariance)

# ╔═╡ 711d46a0-7f0a-11eb-15d7-b5b17afdeeaa
md"Each of the groups have a positive covariance (though not high) between the two variables."

# ╔═╡ 80fc7f8c-7f0a-11eb-3afd-4796a222b202
md"👉Repeat this exercise for sepal length and width variables"

# ╔═╡ 09b71b5a-7e9a-11eb-07dd-e7eabcdb50ef
cov(iris.SepalLength, iris.SepalWidth)

# ╔═╡ 2aaca4f6-7e9a-11eb-3071-fdd5483fb812
scatter(iris.SepalLength, iris.SepalWidth)

# ╔═╡ fb4be732-7f0a-11eb-2b25-b3b18d6cba8a
md"👉So far, the covariance seems like a good measure to relate two variables. Do you see a downside? Hint: When is covariance high?"

# ╔═╡ 031f2f2c-7f0c-11eb-25b7-eb7530aa62af
md"As an example consider the following two sets of covariance computations and then see the corresponding scatter plots."

# ╔═╡ 146d1aaa-7f0c-11eb-00ae-fb7ed270927b
begin
	x4_ = rand(Normal(5, 5), 1000);
	y4_ = 0.1 * x4_ + rand(Normal(5, 5), 1000);
	cov(x4_, y4_)
end

# ╔═╡ 331c74be-7f0c-11eb-1195-cf4e817f61d3
begin
	x5_ = rand(Normal(5, 0.5), 1000);
	y5_ = 1 * x5_ + rand(Normal(5, 0.5), 1000);
	cov(x5_, y5_)
end

# ╔═╡ 3cb5b896-7f0c-11eb-324d-6d745f9c3796
scatter(x4_, y4_, title="Covariance = " * string(cov(x4_, y4_)))

# ╔═╡ 4f001262-7f0c-11eb-2422-1ff23f6065a1
scatter(x5_, y5_, title="Covariance = " * string(cov(x5_, y5_)))

# ╔═╡ 6acc6cfc-7f0c-11eb-3522-174161cba17b
md"Clearly by comparing the covariance between these two sets we come to the wrong conclusion that the relationship is stronger for the first set. We thus need to **normalize** the covariance"

# ╔═╡ a8744e0a-7e95-11eb-1931-b1cdda813e0f
md"#### Statistics of two variables - Correlation"

# ╔═╡ 716cb26e-6b94-11eb-1cb0-b9c2f5608b5f
md"Correlation essentially fixes the problem we discovered - of normalizing covariance . 👉 Can you guess how?"

# ╔═╡ aff1b65e-7f0c-11eb-30de-45877373262c
md"The correlation between two sets of populations $X$ and $Y$ is given as
```math
\rho_{X, Y} = \dfrac{\mathrm{COV}(X, Y)}{\sigma_X \sigma_Y}
```
where $\sigma_X$ and $\sigma_Y$ are the standard deviations of $X$ and $Y$, respectively.
"

# ╔═╡ 3117ff70-7f0d-11eb-1c0a-95af44530f05
md"The corresponding equation for two samples $x$ and $y$ is given as
```math
\rho_{x, y} = \dfrac{\mathrm{COV}(x, y)}{s_x s_y}
```
where $s_x$ and $s_y$ are the standard deviations of $x$ and $y$, respectively.
"

# ╔═╡ a4933666-7f16-11eb-34fa-295373733b14
md"The above definition of correlation is called the _Pearson correlation coefficient_."

# ╔═╡ 5510b24e-7f0d-11eb-15a2-6714879d78f4
md"👉What is the unit of correlation?"

# ╔═╡ 93baa848-7f0e-11eb-2700-bfc3181e0f21
md"👉Class exercise: What are the minimum and maximum possible values of $\rho_{x, y}$?"

# ╔═╡ 9115a164-7f0f-11eb-3946-453ad1a450a0
md"Let us algebraically manipulate the formula of the correlation. Note that for simplicity I am dropping the limits of the summation. In all cases, they are going from $i = 1$ to $i = n$."

# ╔═╡ a731ef30-7f0e-11eb-24f0-5d102ea4a13a
md"
```math
\rho_{x, y} = \dfrac{\mathrm{COV}(x, y)}{s_x s_y} = \dfrac{\left(\sum((x_i - \overline{x})(y_i - \overline{y})\right)/(n - 1)}{\sqrt{\sum(x_i - \overline{x})^2} \sqrt{\sum(y_i - \overline{y})^2}/(n - 1)}
```
```math
\rho_{x, y}^2 = \dfrac{\left(\sum((x_i - \overline{x})(y_i - \overline{y})\right)^2}{\sum(x_i - \overline{x})^2 \sum(y_i - \overline{y})^2} = \dfrac{(\sum a_ib_i)^2}{\sum a_i^2 \sum b_i^2}
```
"

# ╔═╡ 6d579f98-7f0f-11eb-26ea-c381d36153f3
md"The Cauchy Schwarz inequality provides a bound on the RHS. It states:
```math
\left(\sum a_i b_i\right)^2 \leq \left(\sum a_i^2\right) \left(\sum b_i^2\right)
```
For a proof see [this article](https://mathworld.wolfram.com/CauchysInequality.html).
"

# ╔═╡ cdbe10ce-7f0f-11eb-303e-8f9f3527ce25
md"Applying this back, we have $\rho^2_{x, y} \leq 1$ and thus $\rho_{x, y} \in [-1, 1]$"

# ╔═╡ e33ba8e4-7f0f-11eb-36b2-039693234b96
md"👉When are the limits obtained?"

# ╔═╡ ebaeb804-7f0f-11eb-1621-216682cc1e62
md"The equality holds when $a_i = cb_i + \delta$ for all $i$ for some constant $c$. So the extreme values of correlation occur when $\rho_{x, y} = \mathrm{Sign}(c)$ if $\dfrac{(x_i - \overline{x})}{(y_i - \overline{y})} = c$ for all values in the sample"

# ╔═╡ 5adeedd8-7f11-11eb-185e-eb4b7521074b
cslider = @bind c_slider html"<input type=range min=-50 max=45 step=1>"

# ╔═╡ f1b4a228-7f10-11eb-3149-8f86ab1fd053
begin
	x6_ = rand(Normal(), 1000)
	y6_ = x6_ .* c_slider
	cor(x6_, y6_)
end

# ╔═╡ b633765a-7f10-11eb-2c80-f3cd7302ab5a
md"So, $\rho$ is unitless and has maximum and minimum values which allow us to compare across data fields or sets. We can check this with our two datasets from before."

# ╔═╡ e8bdb338-7f10-11eb-0b2d-47d7f358334f
scatter(x4_, y4_, title="Covariance = " * string(round(cov(x4_, y4_) * 1000)/1000) * ". Correlation = " * string(round(cor(x4_, y4_) * 1000)/1000), label=false)

# ╔═╡ e3cf7068-7f11-11eb-0645-1d988c17ea0f
scatter(x5_, y5_, title="Covariance = " * string(round(cov(x5_, y5_) * 1000)/1000) * ". Correlation = " * string(round(cor(x5_, y5_) * 1000)/1000), label=false)

# ╔═╡ 0e4169dc-7f12-11eb-1f09-695a9cbff347
md"The bottom case has now much higher correlation because we have scaled for the low standard deviation of the two individual variables."

# ╔═╡ 6d5a3c1e-7f12-11eb-244f-2787db2ff86b
md"We can use the `cor` function for dataframe columns also."

# ╔═╡ 4f405416-7f12-11eb-072f-b5ebf6ecdb89
cor(iris.PetalLength, iris.PetalWidth)

# ╔═╡ f678bd9a-7f12-11eb-3b82-0748786b4b75
md"For dataframes, it is usually helpful to compute pair-wise covariance and correlation between a set of quantitative variables. There is no built-in method for this, so we will use this one from [this forum](https://discourse.julialang.org/t/covariance-from-dataframe-or-timearray/48375/2)."

# ╔═╡ fd483aba-8584-11eb-1af1-fd7ba22276ba
md"👉What do you expect in the diagonals of such a matrix/"

# ╔═╡ 6c2b36ea-7f12-11eb-3b91-5f6186dd5ce3
function covmat(df)
   nc = ncol(df)
   t = zeros(nc, nc)
   for (i, c1) in enumerate(eachcol(df))
	   for (j, c2) in enumerate(eachcol(df))
		   sx, sy = skipmissings(c1, c2)
		   t[i, j] = cov(collect(sx), collect(sy))
	   end
   end
   return t
   end;

# ╔═╡ 3ae44970-7f13-11eb-3222-25121d19dfa0
covmat(select(iris, Not(:Species)))

# ╔═╡ 1387ed38-7f14-11eb-0f7b-c3deea1d8d76
md"Notice that the matrix is symmetric and the values on the diagonal are all positive"

# ╔═╡ 67b909b0-7f13-11eb-36ca-f74db0715854
propertynames(iris)

# ╔═╡ 11ef7960-8585-11eb-2dc3-7fa091595dbb
md"Now let us do the same with correlation. 👉What do you expect in the diagonals of such a matrix/"

# ╔═╡ 0376056a-7f14-11eb-0af4-45c83b91ff94
function cormat(df)
   nc = ncol(df)
   t = zeros(nc, nc)
   for (i, c1) in enumerate(eachcol(df))
	   for (j, c2) in enumerate(eachcol(df))
		   sx, sy = skipmissings(c1, c2)
		   t[i, j] = cor(collect(sx), collect(sy))
	   end
   end
   return t
   end;

# ╔═╡ 07f0fc76-7f14-11eb-14ef-ed4a723c9f39
cormat(select(iris, Not(:Species)))

# ╔═╡ 1cf0720a-7f14-11eb-28af-c73c1f94fa63
md"Notice that the values on the diagonal are now all 1. The sign of covariance and correlation is the same, but they have all been scaled up. It is easier to conclude on relatedness of variables with the correlation."

# ╔═╡ 681cbce8-7f14-11eb-30d6-17d3f7f5f116
cormat(select(diamonds, [:X, :Y, :Z, :Price]))

# ╔═╡ 5ce49422-7f14-11eb-3959-4fb5d18d1f38
md"The dimensions of the diamonds along different axes are highly correlated. Naturally the price also is correlated with the dimensions."

# ╔═╡ 76ff24b4-7f15-11eb-33c1-6118de89b12f
md"Rule of thumb: If $|\rho| > 0.5$, then we say the two samples are highly correlated. For $0.3 < |\rho| < 0.5$, we say there is moderate correlation. And for $|\rho| < 0.3$, the variables are uncorrelated or weakly correlated."

# ╔═╡ 62023c38-7f15-11eb-039c-e5020554a66f
md"👉What does a high absolute value of the correlation coefficient actually mean?"

# ╔═╡ 71394f5c-7f15-11eb-156d-27290791641f
md"If $|\rho|$ is high then we are confident of a **linear** relationship between the two samples."

# ╔═╡ 1fc41aca-7f16-11eb-331d-072d0e72cf16
md"This means that even if there is an exact _non-linear_ relationship between two variables, the correlation coefficient can be close to 0."

# ╔═╡ 35f0324a-7f16-11eb-159b-1f6056dad3f9
md"Consider the following case where there is a quadratic relationship between the two variables, but the correlation close to 0."

# ╔═╡ 577ceeb8-7f15-11eb-3700-cd396e2622b1
scatter(x5_, (x5_ .- mean(x5_)).^2, title="Correlation = " * string(cor(x5_, (x5_ .- mean(x5_)).^2)))

# ╔═╡ 59545e4e-7f16-11eb-03ae-a796f324110d
md"Thus the correlation coefficient is only a way to check a linear relatedness between variables. Often this itself is very useful in practice. For instance, does adminstering a vaccine reduce the spread of a disease? Or does adding a new modification to the engine increase range of an electric car?"

# ╔═╡ 9f311f58-7f16-11eb-2909-cb836b098441
md"But we must recognise that in most of machine learning today, we are interested in complex non-linear relationships. In such cases, the correlation coefficient may not be of much significance."

# ╔═╡ 3e702d8e-7f33-11eb-271d-5956d16a7e7d
md"👉 Can you think of a linear relationship that you would study in most ML problems?"

# ╔═╡ ed9d08b4-7f32-11eb-04b1-ef5a467dabea
md"👉Do you see a major limitation with this measure?"

# ╔═╡ 09a0093e-7f34-11eb-2210-ff1d72e9637a
md"We know that the mean (and hence offset from mean) and the standard deviation of any sample is adversely affected by outliers. So, it is likely that the Pearson correlation coefficient is also sensitive to outliers.

Notice that the argument that we are dividing the covariance by the two standard deviations is to only allow comparison across variables, not counter any potential effect of outliers."

# ╔═╡ 3fdc0494-7f34-11eb-139d-ad7fa8c098e2
begin
	x7_ = rand(Normal(5, 1), 20);
	y7_ = 1 * x7_ + rand(Normal(5, 0.3), 20);

	x8_ = append!(copy(x7_), 10)
	y8_ = append!(copy(y7_), 9)
end

# ╔═╡ 1eb842ea-7f35-11eb-3e5b-753535a8bfea
md"Click on the checkbox to add outliers"

# ╔═╡ 8608e458-7f34-11eb-3b92-290fea88ac18
scatter(x7_, y7_, title="Correlation = " * string(cor(x7_, y7_)))

# ╔═╡ bc17e118-7f34-11eb-1189-95315fab5ce6
scatter(x8_, y8_, title="Correlation = " * string(cor(x8_, y8_)))

# ╔═╡ 176b14f8-7f36-11eb-1dde-2f1ca7c9af38
md"Correlation drops significantly due to a single outlier"

# ╔═╡ 79d21678-7f36-11eb-3e14-9568a7781ebe
md"👉Can you think of a way in which correlation coefficient increases due to outliers?"

# ╔═╡ 1e02e264-7f36-11eb-12f1-81411831a1b3
begin
	x9_ = rand(Normal(5, 1), 20);
	y9_ = 0.5 * x7_ + rand(Normal(5, 2), 20);

	x10_ = append!(copy(x9_), [20])
	y10_ = append!(copy(y9_), [20])
end

# ╔═╡ a109437e-7f36-11eb-0bdf-ab2bf0b8130a
scatter(x9_, y9_, title="Correlation = " * string(cor(x9_, y9_)))

# ╔═╡ a9195216-7f36-11eb-25cf-ab7eee4ca256
scatter(x10_, y10_, title="Correlation = " * string(cor(x10_, y10_)))

# ╔═╡ ee876554-7f36-11eb-1e02-b98451ccf26f
md"👉OK so we have a problem with outliers. What would your defence be?"

# ╔═╡ 769f6de6-8585-11eb-3a9f-450df3ab75a0
md"💡Hint: What did we do in the case of mean?"

# ╔═╡ 6cd8b81a-7f17-11eb-1346-270ab26a1df4
md"#### Spearman Correlation"

# ╔═╡ e3b09fdc-7f32-11eb-3b6d-2bca422fd7f3
md"Our solution with the outlier sensitivity of mean was to use median which uses the ranking of variables rather than their values."

# ╔═╡ 14ee007e-7f37-11eb-20b4-1df87a90f48d
md"👉 Can you think of correlating two variables with ranks?"

# ╔═╡ 515ed8d2-7f3d-11eb-1f0e-f920f7799263
md"The basic idea is to compute the ranks of two set of samples and then compute the Pearson's correlation coefficient of the ranks."

# ╔═╡ 60057c6a-7f3d-11eb-0f32-595b689315a3
x12_ = rand(Normal(5, 1), 5)

# ╔═╡ afff7064-7f50-11eb-05b1-d392d177ace5
y12_ = 0.5 * x12_ + rand(Normal(5, 2), 5)

# ╔═╡ bc51fde0-7f3d-11eb-02a5-9331000ad2ff
cor(x12_, y12_)

# ╔═╡ ced6f668-7f50-11eb-2fa2-755dc4b48e74
sort(x12_)

# ╔═╡ 5aac63e2-7f50-11eb-1888-a5b1980fea98
[indexin(x12_[i], sort(x12_))[1] for i in 1:length(x12_)]

# ╔═╡ 10f3be2a-7f51-11eb-1d4f-ff3943118e1f
md"The following function replaces each value in an input array by its rank in the sorted version of the array"

# ╔═╡ d8afc632-7f50-11eb-08a5-83914cf962d0
findpos(arr) = [indexin(arr[i], sort(arr))[1] for i in 1:length(arr)]

# ╔═╡ e6cad106-7f50-11eb-1069-8520b81a05fd
x12_rank = findpos(x12_)

# ╔═╡ ed174580-7f50-11eb-29ac-e925387f1779
y12_rank = findpos(y12_)

# ╔═╡ 2222b07a-7f51-11eb-11f9-3f8c191c620e
md"Once we have the ranks for both the vectors, we can compute the Spearman correlation coefficient as the Pearson correlation coefficient on the ranks."

# ╔═╡ fcf549a2-7f50-11eb-05d3-11fbdf7f3984
cor(x12_rank, y12_rank)

# ╔═╡ 4fa95fbc-7f51-11eb-3f79-05e2c0408808
md"We can put all this together in a single function"

# ╔═╡ 3fc8a3a0-7f51-11eb-2a6c-6385c55d2d2d
cor_s(x, y) = cor(findpos(x), findpos(y))

# ╔═╡ 4ec8c7f4-7f51-11eb-133e-bbc12393f533
cor_s(x12_, y12_)

# ╔═╡ 8fbc0f54-7f52-11eb-0ab8-fb60067b43e0
md"Let us now revisit the examples where outliers sensitively affected the Pearson correlation coefficient"

# ╔═╡ a97b5ff8-7f52-11eb-220e-55f4aceb21de
md"Case 1: ρ-Pearson sharply reduced by a single outlier."

# ╔═╡ b2052776-7f52-11eb-3790-ff5c917d9da7
md"Case 2: ρ-Pearson sharply increased by an outlier"

# ╔═╡ ddad5ac6-7f52-11eb-148b-91e5b7a5df8f
md"In both cases the Spearman correlation coefficient remains more robust in the face of the outliers. This is clear why - the outliers have only limited effect on the rank of the values which are then used to compute correlation."

# ╔═╡ c17e66ac-8585-11eb-2cc8-3fd3501f33ce
md"👉Do you see one other advantage of Spearman correlation coefficient over the Pearson correlation coefficient?"

# ╔═╡ 0959a256-7f53-11eb-2afe-c90157d032f3
md"Since the Spearman correlation depends only on the rank ordering it can also detect correlations with non-linear _monotonic_ functions."

# ╔═╡ 73f12006-7f53-11eb-23e9-47e725e79336
begin
	x16_ = rand(Uniform(-10, 10), 100)
	y16_ = sinh.(x16_)
end

# ╔═╡ 5684ace4-7f54-11eb-197a-bd770705eced
md"👉Where is this particularly important?"

# ╔═╡ 0868668c-7f53-11eb-0025-eb7e12f385dc
hint(text, title) = Markdown.MD(Markdown.Admonition("hint", title, [text]))

# ╔═╡ c6efb902-8255-11eb-16f0-6ffc51cf9224
hint(md"
```julia
begin
	sample_ = rand(Normal(), 100)
	val_ = sum(sum([[(sample_[i] - sample_[j])^2 for i in 1:100] for j in 1:100]))/(100 * 100)
end
```", "Solution")

# ╔═╡ 1dc8ee70-8255-11eb-3d0e-53dda549c76b
hint(md"$\mathbb{E}[(X - Y)^2] = \mathbb{E}[X^2] + \mathbb{E}[Y^2] - 2\mathbb{E}[XY]$ But we know that $X$ and $Y$ are independent and so $\mathbb{E}[XY] = \mathbb{E}[X]\mathbb{E}[Y] = \mu^2$. We also know that $\mathbb{E}[X^2] = {E}[Y^2] = \mu^2 + \sigma^2$. Hence, the above expectation is equal to $2\sigma^2$. Hence, it makes sense to get values around 2 for the previous code.", "Solution")

# ╔═╡ 5acfe6ae-7f06-11eb-182b-9b83107e0c4c
hint(md"
```math
\mathbb{E}[(X - \mathbb{E}[X])(Y - \mathbb{E}[Y])] = \mathbb{E}[XY - X\mathbb{E}[Y] - Y\mathbb{E}[X] + \mathbb{E}[X]\mathbb{E}[Y]]
```
```math
= \mathbb{E}[XY] - \mathbb{E}[X]\mathbb{E}[Y] - \mathbb{E}[Y]\mathbb{E}[X] + \mathbb{E}[X]\mathbb{E}[Y] =  \mathbb{E}[XY] - \mathbb{E}[X]\mathbb{E}[Y]
```", "Rewriting covariance")

# ╔═╡ ac68800c-7f38-11eb-0c1b-79e8d2ab8ecf
hint(md"
```math
\mathrm{COV}(X + Y, Z) = \mathbb{E}[(X + Y)Z] - \mathbb{E}[X+Y]\mathbb{E}[Z]
```
```math
= \mathbb{E}[XZ] + \mathbb{E}[YZ] - \mathbb{E}[X]\mathbb{E}[Z] - \mathbb{E}[Y]\mathbb{E}[Z]
```
```math
= (\mathbb{E}[XZ] - \mathbb{E}[X]\mathbb{E}[Z]) + (\mathbb{E}[YZ]  - \mathbb{E}[Y]\mathbb{E}[Z])
```
```math
= \mathrm{COV}(X, Z) + \mathrm{COV}(Y, Z)
```", "Covariance sum")

# ╔═╡ baa84258-7f0b-11eb-03d4-c5e0d9547508
hint(md"Covariance is high when the signs of the offsets of the two variable match often (good part), but also when the two variables have large variances individually (bad part)", "Issue with covariance")

# ╔═╡ c6bd6268-7f0d-11eb-2094-0fdc2a8fbd88
hint(md"Correlation is unit-less as we have scaled the product of offsets from mean by the product of two standard deviations.", "Unit of covariance")

# ╔═╡ 48065f1c-7f33-11eb-1888-13ba63eea532
hint(md"Usually whenever we do regression, we have the actual value and the estimate. And we do expect an identity relation between them. In such cases we use the coefficient of determination $R^2$ as a goodness of fit. But in some cases we are fine with a high correlation between the actual and estimated value.", "Use of correlation in ML")

# ╔═╡ 5ff7afb0-7f54-11eb-0f8c-bb7af9a0de90
hint(md"Let us say we are building a ML estimator for efficacy of different manufacturing parameters for yield. Even if this function is non-linear but monotonic in the design parameters, then we can use the Spearman correlation coefficient to set the parameters for increasing the yield..", "Spearman correlation - usage")

# ╔═╡ 14eca90a-7f52-11eb-1257-c3eb50212780
disp(number) = string(round(number*1000)/1000)

# ╔═╡ 58917ca4-7f51-11eb-2399-53cc4269ba1b
scatter(x7_, y7_, title="ρ-Pearson = " * disp(cor(x7_, y7_)) * " ρ-Spearman = " * disp(cor_s(x7_, y7_)))

# ╔═╡ 37d45bea-7f52-11eb-2302-2da751a91c33
scatter(x8_, y8_, title="ρ-Pearson = " * disp(cor(x8_, y8_)) * " ρ-Spearman = " * disp(cor_s(x8_, y8_)))

# ╔═╡ c1ab36ca-7f52-11eb-2965-8914dacd089b
scatter(x9_, y9_, title="ρ-Pearson = " * disp(cor(x9_, y9_)) * " ρ-Spearman = " * disp(cor_s(x9_, y9_)))

# ╔═╡ cc35a918-7f52-11eb-1bc5-db77f47abf9f
scatter(x10_, y10_, title="ρ-Pearson = " * disp(cor(x10_, y10_)) * " ρ-Spearman = " * disp(cor_s(x10_, y10_)))

# ╔═╡ 8494b7a6-7f53-11eb-0554-e3f7afb68fc5
scatter(x16_, y16_, title="ρ-Pearson = " * disp(cor(x16_, y16_)) * " ρ-Spearman = " * disp(cor_s(x16_, y16_)), legend=false)

# ╔═╡ Cell order:
# ╟─253ad4f4-8253-11eb-03d8-d358d6075910
# ╟─c6efb902-8255-11eb-16f0-6ffc51cf9224
# ╠═237d9290-8254-11eb-233a-4d362783d684
# ╟─8ce2cb16-8253-11eb-21a5-c59e5785950e
# ╟─1dc8ee70-8255-11eb-3d0e-53dda549c76b
# ╟─22f18e82-6b7d-11eb-16eb-31758abdb1b0
# ╟─59607216-7b43-11eb-2cf7-3f74231a1839
# ╠═def88812-6b99-11eb-07fb-3f961d7ff6ae
# ╠═67bd1b94-7f32-11eb-2869-5f075fca7fbc
# ╠═5f510e34-7f32-11eb-158a-afdded4ae4f7
# ╠═65c3f4a2-7f32-11eb-0eb4-97a2455d1ac7
# ╟─b1716986-7e95-11eb-383d-39edf8743454
# ╟─67e698c8-7e99-11eb-01c9-eb2703970ea4
# ╟─c6dc30e0-7e99-11eb-1781-7dba5912b724
# ╟─b79f4c20-7e99-11eb-3eba-492dc528735d
# ╠═7b13d1a8-7e9a-11eb-39f2-5d3e9a95c811
# ╟─a136b75a-7f04-11eb-1324-c9ea26353914
# ╠═ff01b9c2-7f02-11eb-2edd-5131a142f119
# ╟─bbe2eec0-7f04-11eb-3f7f-594d57f134ec
# ╠═fc068802-7f03-11eb-36a0-f3764ac250ab
# ╟─eb742c80-7f04-11eb-0435-4f58142bcdc9
# ╠═2f64101e-7f04-11eb-148b-d1c781173d4c
# ╟─98fcccfc-7f04-11eb-0b2d-6b1feaadfae4
# ╠═12e7d1a4-7f05-11eb-19d5-2b072cbafa05
# ╟─2c8776c8-7f05-11eb-0226-bf1ef7542efd
# ╟─9da0add4-7f05-11eb-1391-0f87db7d41c4
# ╟─a86b8dc4-7f05-11eb-3ca7-e5460fc15b96
# ╟─2e4af178-7f06-11eb-1247-3501f18174cc
# ╟─5acfe6ae-7f06-11eb-182b-9b83107e0c4c
# ╟─e79a0a24-7f06-11eb-3c63-91d44fe78d49
# ╟─2be9abbe-824d-11eb-3b1d-2b1ddf16c195
# ╟─483778b4-824d-11eb-196c-2925d31dd603
# ╟─b818699a-7f37-11eb-088c-a51bed3cbdc0
# ╟─41664208-7f38-11eb-01f5-c32a33253f15
# ╟─6dfb60e6-7f38-11eb-070d-b99afff5f5a1
# ╟─ac68800c-7f38-11eb-0c1b-79e8d2ab8ecf
# ╟─2e32a1e2-7f39-11eb-0df0-812f689e1d1f
# ╟─c4d2c9d0-7f39-11eb-2ac9-373ed723a3f4
# ╟─ce6afe4a-7f39-11eb-19f1-7157457700cd
# ╟─175b1608-7f3a-11eb-0dfb-dd2ee73b9bac
# ╟─fe62b640-8251-11eb-25fe-2b8aa69000ce
# ╠═27d06de0-7f07-11eb-23ce-4d6547ce31f6
# ╠═0359bc70-7f07-11eb-2daa-d1f2a258d9e0
# ╠═a8c20bc6-8251-11eb-0d06-819a786e2b05
# ╟─9c42421a-7f0d-11eb-0172-f34b4adfa5f0
# ╟─502bc23a-7f07-11eb-0b40-97e7e2a8a9ca
# ╠═2938c2ac-7f07-11eb-279c-ab02e04ac679
# ╟─59828c88-7f07-11eb-1462-73e3850659fd
# ╟─a48951d0-7f07-11eb-0052-7d06cf26028c
# ╟─e0d54d10-7f07-11eb-214a-bd24b8accb1a
# ╠═ec51f080-7f07-11eb-2559-5b4cf71a2d0e
# ╟─ef48b742-7f07-11eb-18e0-f1bcc96c86e0
# ╠═dc4d1258-7e99-11eb-25e2-ffcd05171460
# ╟─fa8d1044-7f07-11eb-247e-a18ecff0f4f3
# ╠═182c2ca2-7e9a-11eb-040a-21429d245b4f
# ╟─2348b81e-7f08-11eb-3070-25becdb70c5a
# ╠═85e8d916-7f08-11eb-094c-17dbb38c5048
# ╠═ed8c0606-7f08-11eb-13ec-d166c29a4d67
# ╠═175b39a2-7f09-11eb-02d7-c91ba9204c3c
# ╠═1d4a9088-7f09-11eb-3754-07c4947ecad7
# ╠═27426ea6-7f09-11eb-3c2f-c9d5d5ed33b8
# ╠═307bf5dc-7f09-11eb-3e23-8908558e17ca
# ╟─5a6fb146-7f09-11eb-02b3-0b2bc301eb95
# ╟─882610e4-7f09-11eb-142b-9d8ed5e88fcc
# ╠═0019040c-7e9a-11eb-3ab3-698eb0f98231
# ╠═241c6342-7e9a-11eb-2053-c1984c19a023
# ╟─54ccb922-7f0a-11eb-3ba6-230b0489d928
# ╟─6273fb94-7f0a-11eb-3c1a-4dfa9d925400
# ╠═b9d98800-7f09-11eb-3671-ed27dbde1c3a
# ╟─711d46a0-7f0a-11eb-15d7-b5b17afdeeaa
# ╟─80fc7f8c-7f0a-11eb-3afd-4796a222b202
# ╠═09b71b5a-7e9a-11eb-07dd-e7eabcdb50ef
# ╠═2aaca4f6-7e9a-11eb-3071-fdd5483fb812
# ╟─fb4be732-7f0a-11eb-2b25-b3b18d6cba8a
# ╟─baa84258-7f0b-11eb-03d4-c5e0d9547508
# ╟─031f2f2c-7f0c-11eb-25b7-eb7530aa62af
# ╠═146d1aaa-7f0c-11eb-00ae-fb7ed270927b
# ╠═331c74be-7f0c-11eb-1195-cf4e817f61d3
# ╠═3cb5b896-7f0c-11eb-324d-6d745f9c3796
# ╠═4f001262-7f0c-11eb-2422-1ff23f6065a1
# ╟─6acc6cfc-7f0c-11eb-3522-174161cba17b
# ╟─a8744e0a-7e95-11eb-1931-b1cdda813e0f
# ╟─716cb26e-6b94-11eb-1cb0-b9c2f5608b5f
# ╟─aff1b65e-7f0c-11eb-30de-45877373262c
# ╟─3117ff70-7f0d-11eb-1c0a-95af44530f05
# ╟─a4933666-7f16-11eb-34fa-295373733b14
# ╟─5510b24e-7f0d-11eb-15a2-6714879d78f4
# ╟─c6bd6268-7f0d-11eb-2094-0fdc2a8fbd88
# ╟─93baa848-7f0e-11eb-2700-bfc3181e0f21
# ╟─9115a164-7f0f-11eb-3946-453ad1a450a0
# ╟─a731ef30-7f0e-11eb-24f0-5d102ea4a13a
# ╟─6d579f98-7f0f-11eb-26ea-c381d36153f3
# ╟─cdbe10ce-7f0f-11eb-303e-8f9f3527ce25
# ╟─e33ba8e4-7f0f-11eb-36b2-039693234b96
# ╟─ebaeb804-7f0f-11eb-1621-216682cc1e62
# ╠═5adeedd8-7f11-11eb-185e-eb4b7521074b
# ╠═f1b4a228-7f10-11eb-3149-8f86ab1fd053
# ╟─b633765a-7f10-11eb-2c80-f3cd7302ab5a
# ╟─e8bdb338-7f10-11eb-0b2d-47d7f358334f
# ╟─e3cf7068-7f11-11eb-0645-1d988c17ea0f
# ╟─0e4169dc-7f12-11eb-1f09-695a9cbff347
# ╟─6d5a3c1e-7f12-11eb-244f-2787db2ff86b
# ╠═4f405416-7f12-11eb-072f-b5ebf6ecdb89
# ╟─f678bd9a-7f12-11eb-3b82-0748786b4b75
# ╟─fd483aba-8584-11eb-1af1-fd7ba22276ba
# ╠═6c2b36ea-7f12-11eb-3b91-5f6186dd5ce3
# ╠═3ae44970-7f13-11eb-3222-25121d19dfa0
# ╟─1387ed38-7f14-11eb-0f7b-c3deea1d8d76
# ╠═67b909b0-7f13-11eb-36ca-f74db0715854
# ╟─11ef7960-8585-11eb-2dc3-7fa091595dbb
# ╠═0376056a-7f14-11eb-0af4-45c83b91ff94
# ╠═07f0fc76-7f14-11eb-14ef-ed4a723c9f39
# ╟─1cf0720a-7f14-11eb-28af-c73c1f94fa63
# ╠═681cbce8-7f14-11eb-30d6-17d3f7f5f116
# ╟─5ce49422-7f14-11eb-3959-4fb5d18d1f38
# ╟─76ff24b4-7f15-11eb-33c1-6118de89b12f
# ╟─62023c38-7f15-11eb-039c-e5020554a66f
# ╟─71394f5c-7f15-11eb-156d-27290791641f
# ╟─1fc41aca-7f16-11eb-331d-072d0e72cf16
# ╟─35f0324a-7f16-11eb-159b-1f6056dad3f9
# ╠═577ceeb8-7f15-11eb-3700-cd396e2622b1
# ╟─59545e4e-7f16-11eb-03ae-a796f324110d
# ╟─9f311f58-7f16-11eb-2909-cb836b098441
# ╟─3e702d8e-7f33-11eb-271d-5956d16a7e7d
# ╟─48065f1c-7f33-11eb-1888-13ba63eea532
# ╟─ed9d08b4-7f32-11eb-04b1-ef5a467dabea
# ╟─09a0093e-7f34-11eb-2210-ff1d72e9637a
# ╠═3fdc0494-7f34-11eb-139d-ad7fa8c098e2
# ╟─1eb842ea-7f35-11eb-3e5b-753535a8bfea
# ╠═8608e458-7f34-11eb-3b92-290fea88ac18
# ╠═bc17e118-7f34-11eb-1189-95315fab5ce6
# ╟─176b14f8-7f36-11eb-1dde-2f1ca7c9af38
# ╟─79d21678-7f36-11eb-3e14-9568a7781ebe
# ╠═1e02e264-7f36-11eb-12f1-81411831a1b3
# ╠═a109437e-7f36-11eb-0bdf-ab2bf0b8130a
# ╠═a9195216-7f36-11eb-25cf-ab7eee4ca256
# ╟─ee876554-7f36-11eb-1e02-b98451ccf26f
# ╟─769f6de6-8585-11eb-3a9f-450df3ab75a0
# ╟─6cd8b81a-7f17-11eb-1346-270ab26a1df4
# ╟─e3b09fdc-7f32-11eb-3b6d-2bca422fd7f3
# ╟─14ee007e-7f37-11eb-20b4-1df87a90f48d
# ╟─515ed8d2-7f3d-11eb-1f0e-f920f7799263
# ╠═60057c6a-7f3d-11eb-0f32-595b689315a3
# ╠═afff7064-7f50-11eb-05b1-d392d177ace5
# ╠═bc51fde0-7f3d-11eb-02a5-9331000ad2ff
# ╠═ced6f668-7f50-11eb-2fa2-755dc4b48e74
# ╠═5aac63e2-7f50-11eb-1888-a5b1980fea98
# ╟─10f3be2a-7f51-11eb-1d4f-ff3943118e1f
# ╠═d8afc632-7f50-11eb-08a5-83914cf962d0
# ╠═e6cad106-7f50-11eb-1069-8520b81a05fd
# ╠═ed174580-7f50-11eb-29ac-e925387f1779
# ╟─2222b07a-7f51-11eb-11f9-3f8c191c620e
# ╠═fcf549a2-7f50-11eb-05d3-11fbdf7f3984
# ╟─4fa95fbc-7f51-11eb-3f79-05e2c0408808
# ╠═3fc8a3a0-7f51-11eb-2a6c-6385c55d2d2d
# ╠═4ec8c7f4-7f51-11eb-133e-bbc12393f533
# ╟─8fbc0f54-7f52-11eb-0ab8-fb60067b43e0
# ╟─a97b5ff8-7f52-11eb-220e-55f4aceb21de
# ╟─58917ca4-7f51-11eb-2399-53cc4269ba1b
# ╟─37d45bea-7f52-11eb-2302-2da751a91c33
# ╟─b2052776-7f52-11eb-3790-ff5c917d9da7
# ╟─c1ab36ca-7f52-11eb-2965-8914dacd089b
# ╟─cc35a918-7f52-11eb-1bc5-db77f47abf9f
# ╟─ddad5ac6-7f52-11eb-148b-91e5b7a5df8f
# ╟─c17e66ac-8585-11eb-2cc8-3fd3501f33ce
# ╟─0959a256-7f53-11eb-2afe-c90157d032f3
# ╠═73f12006-7f53-11eb-23e9-47e725e79336
# ╠═8494b7a6-7f53-11eb-0554-e3f7afb68fc5
# ╟─5684ace4-7f54-11eb-197a-bd770705eced
# ╟─5ff7afb0-7f54-11eb-0f8c-bb7af9a0de90
# ╠═0868668c-7f53-11eb-0025-eb7e12f385dc
# ╠═14eca90a-7f52-11eb-1257-c3eb50212780

ext: jl


text:
struct F84abs <: F84
  κ::Float64
  β::Float64
  πA::Float64
  πC::Float64
  πG::Float64
  πT::Float64
  function F84abs(κ::Float64, β::Float64,
                  πA::Float64, πC::Float64, πG::Float64, πT::Float64)
    if κ <= 0.
      error("F84 parameter α must be positive")
    elseif β <= 0.
      error("F84 parameter β must be positive")
    elseif sum([πA,πC,πG,πT]) != 1.0
      error("F84 frequencies must sum to 1.0")
    elseif any([πA,πC,πG,πT] .<= 0.0)
      error("F84 frequencies must be positive")
    end
    new(κ, β, πA, πC, πG, πT)
  end
end

function show(io::IO, object::F84abs)
  print(io, "\r\e[0m\e[1mF\e[0melsenstein 19\e[1m84\e[0m substitution model (absolute rate form)
κ = $(object.κ), β = $(object.β), π = [$(object.πA), $(object.πC), $(object.πG), $(object.πT)]")
end

F84(κ, β, πA, πC, πG, πT) = F84abs(κ, β, πA, πC, πG, πT)

@inline function Q(mod::F84abs)
  β = mod.β; κ = mod.κ
  πA = _πA(mod); πC = _πC(mod); πG = _πG(mod); πT = _πT(mod)
  πR = _πR(mod); πY = _πY(mod)

  α₁ = (1.0 + κ / πY) * β
  α₂ = (1.0 + κ / πR) * β

  Q₁  = β * πA
  Q₂  = α₂ * πA
  Q₃  = β * πC
  Q₄  = α₁ * πC
  Q₅  = α₂ * πG
  Q₆  = β * πG
  Q₇  = β * πT
  Q₈  = α₁ * πT
  Q₉  = -(Q₃ + Q₅ + Q₇)
  Q₁₀ = -(Q₁ + Q₆ + Q₈)
  Q₁₁ = -(Q₂ + Q₃ + Q₇)
  Q₁₂ = -(Q₁ + Q₄ + Q₆)

  return Qmatrix(Q₉,  Q₁,  Q₂,  Q₁,
                 Q₃,  Q₁₀, Q₃,  Q₄,
                 Q₅,  Q₆,  Q₁₁, Q₆,
                 Q₇,  Q₈,  Q₇,  Q₁₂)
end

@inline function P(mod::F84abs, t::Float64)
  if t < 0.0
    error("t must be positive")
  end
  β = mod.β; κ = mod.κ
  πA = _πA(mod); πC = _πC(mod); πG = _πG(mod); πT = _πT(mod)
  πR = _πR(mod); πY = _πY(mod)

  α₁ = (1.0 + κ / πY) * β
  α₂ = (1.0 + κ / πR) * β

  e₁ = exp(-β * t)
  e₂ = exp(-(πR * α₂ + πY * β) * t)
  e₃ = exp(-(πY * α₁ + πR * β) * t)

  P₁  = πA + (πA * πY / πR) * e₁ + (πG / πR) * e₂
  P₂  = πC + (πC * πR / πY) * e₁ + (πT / πY) * e₃
  P₃  = πG + (πG * πY / πR) * e₁ + (πA / πR) * e₂
  P₄  = πT + (πT * πR / πY) * e₁ + (πC / πY) * e₃
  P₅  = πA * (1 - e₁)
  P₆  = πA + (πA * πY / πR) * e₁ - (πA / πR) * e₂
  P₇  = πC * (1 - e₁)
  P₈  = πC + (πC * πR / πY) * e₁ - (πC / πY) * e₃
  P₉  = πG + (πG * πY / πR) * e₁ - (πG / πR) * e₂
  P₁₀ = πG * (1 - e₁)
  P₁₁ = πT * (1 - e₁)
  P₁₂ = πT + (πT * πR / πY) * e₁ - (πT / πY) * e₃

  return Pmatrix(P₁,  P₅,  P₆,  P₅,
                 P₇,  P₂,  P₇,  P₈,
                 P₉,  P₁₀, P₃,  P₁₀,
                 P₁₁, P₁₂, P₁₁, P₄)
end


text:
module Activations

export Activation, forward, backward, parms, ∇parms, Sigmoid, Tanh, ReLU

include("./tensor_dt.jl")
using ..TensorDT: Tensor

include("./abstract_layers.jl")
import ..AbstractLayers: AbstractLayer, AL, forward, backward, parms, ∇parms

## ======================================================================
## Activation Layer
## ======================================================================
struct Activation <: AbstractLayer
  name::Symbol
  fn::Function
  der_fn::Function
  #
  _type::Symbol
  store::Dict

  function Activation(name, fn, der_fn)
    new(name, fn, der_fn, :activation, Dict())
  end
end

function forward(self::Activation, input::Tensor)::Tensor
  self.store[self.name] = self.fn.(input)
  self.store[self.name]
end

function backward(self::Activation, ∇p::Tensor)::Tensor
  self.der_fn(self.store[self.name]) .* ∇p
end

parms(self::Activation) = []
∇parms(self::Activation) = []

### Sigmoid (Activation) function
function Sigmoid()
  σ = z -> 1. ./ (1. .+ exp.(-z))
  der_σ = z -> σ(z) .* (1 .- σ(z))
  Activation(:sigmoid, σ, der_σ)
end

### Tanh (Activation) function
function Tanh()
  tanₕ = z -> (x = exp.(z); y = exp.(-z); (x .- y) ./ (x .+ y))
  der_tanₕ = z -> 1 .- tanₕ.(z) .^ 2
  Activation(:tanh, tanₕ, der_tanₕ)
end

### ReLU (Activation) function
function ReLU()
  relu_fn = z -> max(zero(eltype(z)), z)
  der_relu = z -> (z₀ =	zero(eltype(z)); z .< z₀ ? z₀ : one(eltype(z)))
  Activation(:relu, relu_fn, der_relu)
end

end  ## Module

max_stars_repo_path: Data Science From Scratch/src/activations.jl

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks


text:
MOI.copy!(dest::LinQuadOptimizer, src::MOI.ModelLike; copynames=false)    = MOIU.defaultcopy!(dest, src, copynames)
MOI.copy!(dest::MOI.ModelLike,    src::LinQuadOptimizer; copynames=false) = MOIU.defaultcopy!(dest, src, copynames)
MOI.copy!(dest::LinQuadOptimizer, src::LinQuadOptimizer; copynames=false) = MOIU.defaultcopy!(dest, src, copynames)

MOI.canget(::LinQuadOptimizer, ::MOI.ListOfVariableAttributesSet) = true
MOI.get(::LinQuadOptimizer, ::MOI.ListOfVariableAttributesSet) = MOI.AbstractVariableAttribute[]
MOI.canget(::LinQuadOptimizer, ::MOI.ListOfModelAttributesSet) = true
MOI.get(m::LinQuadOptimizer, ::MOI.ListOfModelAttributesSet) = [
    MOI.ObjectiveSense(),
    MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
    (MOI.get(m, MOI.Name()) != "") ? MOI.Name() : nothing
]
MOI.canget(::LinQuadOptimizer, ::MOI.ListOfConstraintAttributesSet{F,S}) where F where S = true
MOI.get(::LinQuadOptimizer, ::MOI.ListOfConstraintAttributesSet{F,S}) where F where S = MOI.AbstractConstraintAttribute[]


text:
module MicrobiomeMixedModels

using Mamba
using Distributions
using LinearAlgebra
using YAML
using RCall

using StatsBase: Histogram, Weights, sample
using DataFrames: DataFrame, rename!, names!

import Base: convert, maximum
import Distributions: params, partype, mean, var, mode, pdf, logpdf, rand
import StatsBase: fit

export
    GeneralizedInverseGaussian,
    params,
    partype,
    mean,
    var,
    mode,
    pdf,
    logpdf,
    rand,
    MIMIX,
    MIMIXGaussian,
    MIMIXNoFactors,
    permanova,
    fit,
    get_model,
    get_inits,
    get_post,
    clr,
    iclr,
    parse_config,
    load_config

abstract type MicrobiomeModel end

struct MIMIX <: MicrobiomeModel
    factors::Int
    function MIMIX(factors::Int)
        @assert factors > zero(factors)
        new(factors)
    end
end

struct MIMIXGaussian <: MicrobiomeModel
    factors::Int
    function MIMIXGaussian(factors::Int)
        @assert factors > zero(factors)
        new(factors)
    end
end

struct MIMIXNoFactors <: MicrobiomeModel end

include("generalized-inverse-gaussian.jl")  # distribution used in MIMIX
include(joinpath("models", "mimix.jl"))
include(joinpath("models", "mimix-gaussian.jl"))
include(joinpath("models", "mimix-no-factors.jl"))
include(joinpath("models", "permanova.jl"))
include("utils.jl")

end  # module

max_stars_repo_name: nsgrantham/mimix

max_issues_repo_name: nsgrantham/mimix

max_forks_repo_name: nsgrantham/mimix


text:
# maxabstest.jl

using BenchmarkTools

N = isempty(ARGS) ? 500 : parse(Int, ARGS[1])

function max_abs(a, b)
    m = zero(promote_type(eltype(a), eltype(b)))
    for i in eachindex(a, b)
        tmp = abs(a[i] - b[i])
        tmp > m && (m = tmp)
    end
    m
end

a = randn(N, N)
b = randn(N, N)

@show VERSION
@show Threads.nthreads()
println()
@show N
print("simple for loop:     ")
@btime max_abs($a, $b)
print("mapreduce abs∘- max: ")
@btime mapreduce(abs∘-, max, $a, $b)
print("maximum(abs, a - b): ")
@btime maximum(abs, $a - $b)
print("maximum(generator):  ")
@btime maximum(abs(i - j) for (i, j) in zip($a, $b))
print("maximum splat(abs∘-):")
@btime maximum(Base.splat(abs∘-), zip($a, $b))
using Tullio
print("Tullio:              ")
@btime @tullio (max) _ := abs($a[i] - $b[i])
print("Tullio (LoopVect.):  ")
using LoopVectorization
@btime @tullio (max) _ := abs($a[i] - $b[i])

function max_abs_turbo(a, b)
    m = zero(promote_type(eltype(a), eltype(b)))
    @turbo for i in eachindex(a, b)
        m = max(m, abs(a[i] - b[i]))
    end
    m
end

function max_abs_tturbo(a, b)
    m = zero(promote_type(eltype(a), eltype(b)))
    @tturbo for i in eachindex(a, b)
        m = max(m, abs(a[i] - b[i]))
    end
    m
end

print("LoopVect. @turbo:    ")
@btime max_abs_turbo($a, $b)
print("LoopVect. @tturbo:   ")
@btime max_abs_tturbo($a, $b)

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
module FractionalDelayFilter

using Distributions
using LinearAlgebra
using DSP

export filtord, getfdfcoef, fdfilter, fdfilter!

function filtord(D; maxorder=-1)
    # Input arguments:
    #   - D: delay in sample
    # Output: Filter order that satisfies causal filter
    if maxorder < 0
        return convert(Int, floor(2D+1))
    else
        return min(convert(Int, floor(2D+1)), maxorder)
    end
end

function getfdfcoef(N, D)
    # Input arguments:
    #   - N: Filter Order
    #   - D: delay in sample
    # Output:
    #   - fdfcoef: Filter coefficient
    #   - M: Non-zero sample index

    U = vandermonde(N+1)
    Q = invcramer(U)  # inverse matrix by using Cramer's rule
    T = transmat(N)
    Qf = T*Q

    # Fractional delay: d
    # The first non-zero sample index: M
    d = D%1
    if mod(N,2)==0  # even
        M = convert(Int, round(D+0.5,RoundNearestTiesAway) - N/2)
    else  # odd
        M = convert(Int, floor(D) - (N-1)/2)
    end
    fdfcoef = transpose(transpose(d.^(0:N)) * Qf)
    return (fdfcoef, M)
end

function fdfilter(in, fdfcoef)
    # By using the modified Farrow structure, a value of inputted delay is updated as follow: D'=D+1
    # This means that the first value of the output of the filter is useless.
    return conv(in, fdfcoef)[2:end]
end

function fdfilter!(delayedsignal, in, fdfcoef, nonzeroindex)
    delayedsignal[1+nonzeroindex:length(in)+nonzeroindex+length(fdfcoef)-2] = conv(in, fdfcoef)[2:end]
end

function fdfilter(in::Array{T}, fdfcoef, nonzeroindex) where T
    out = zeros(T, length(in)+length(fdfcoef)-2+nonzeroindex)
    fdfilter!(out, in, fdfcoef, nonzeroindex)
    return out
end

function vandermonde(N::Int)
    out = repeat(convert.(Float64, collect(0:N-1)), 1, N)
    for n = 0:N-1
        @. out[:,n+1] ^= n
    end
    return out
end

function invcramer(in::Array{T}) where T
    N = size(in,1)
    out = zeros(T, N, N)
    for m = 1:N
        for n = 1:N
            out[n,m] = (-1)^(m+n) * det(in[[1:m-1; m+1:N],[1:n-1; n+1:N]])
        end
    end
    out /= det(in)
    return out
end

function transmat(N::Int)
    T = zeros(Float64, N+1, N+1)
    for n = 0:N
        for m = 0:N
            if n >= m
                T[m+1,n+1] = round(N/2,RoundNearestTiesAway)^(n-m) * binomial(n,m)
            end
        end
    end
    return T
end

end


text:
struct ODEBacksolveSensitivityFunction{C<:AdjointDiffCache,Alg<:BacksolveAdjoint,uType,pType,fType<:DiffEqBase.AbstractDiffEqFunction} <: SensitivityFunction
  diffcache::C
  sensealg::Alg
  discrete::Bool
  y::uType
  prob::pType
  f::fType
  noiseterm::Bool
end

function ODEBacksolveSensitivityFunction(g,sensealg,discrete,sol,dg,f;noiseterm=false)
  diffcache, y = adjointdiffcache(g,sensealg,discrete,sol,dg,f;quad=false,noiseterm=noiseterm)

  return ODEBacksolveSensitivityFunction(diffcache,sensealg,discrete,
                                         y,sol.prob,f,noiseterm)
end

function (S::ODEBacksolveSensitivityFunction)(du,u,p,t)
  @unpack y, prob, discrete = S

  λ,grad,_y,dλ,dgrad,dy = split_states(du,u,t,S)

  copyto!(vec(y), _y)

  if S.noiseterm
    if length(u) == length(du)
      vecjacobian!(dλ, y, λ, p, t, S, dgrad=dgrad,dy=dy)
    elseif length(u) != length(du) &&  StochasticDiffEq.is_diagonal_noise(prob) && !isnoisemixing(S.sensealg)
      vecjacobian!(dλ, y, λ, p, t, S, dy=dy)
      jacNoise!(λ, y, p, t, S, dgrad=dgrad)
    else
      jacNoise!(λ, y, p, t, S, dgrad=dgrad, dλ=dλ, dy=dy)
    end
  else
    vecjacobian!(dλ, y, λ, p, t, S, dgrad=dgrad, dy=dy)
  end
  dλ .*= -1

  discrete || accumulate_cost!(dλ, y, p, t, S, dgrad)
  return nothing
end

# u = λ' # for the RODE case
function (S::ODEBacksolveSensitivityFunction)(du,u,p,t,W)
  @unpack y, prob, discrete = S

  λ,grad,_y,dλ,dgrad,dy = split_states(du,u,t,S)
  copyto!(vec(y), _y)

  vecjacobian!(dλ, y, λ, p, t, S, dgrad=dgrad, dy=dy,W=W)
  dλ .*= -one(eltype(λ))

  discrete || accumulate_cost!(dλ, y, p, t, S, dgrad)
  return nothing
end

function split_states(du,u,t,S::ODEBacksolveSensitivityFunction;update=true)
  @unpack y, prob = S
  idx = length(y)

  λ     = @view u[1:idx]
  grad  = @view u[idx+1:end-idx]
  _y    = @view u[end-idx+1:end]

  if length(u) == length(du)
    # ODE/Drift term and scalar noise
    dλ    = @view du[1:idx]
    dgrad = @view du[idx+1:end-idx]
    dy    = @view du[end-idx+1:end]

  elseif length(u) != length(du) &&  StochasticDiffEq.is_diagonal_noise(prob) && !isnoisemixing(S.sensealg)
    # Diffusion term, diagonal noise, length(du) =  u*m
    idx1 = [length(u)*(i-1)+i for i in 1:idx] # for diagonal indices of [1:idx,1:idx]
    idx2 = [(length(u)+1)*i-idx for i in 1:idx] # for diagonal indices of [end-idx+1:end,1:idx]

    dλ    = @view du[idx1]
    dgrad = @view du[idx+1:end-idx,1:idx]
    dy    = @view du[idx2]

  elseif length(u) != length(du) &&  StochasticDiffEq.is_diagonal_noise(prob) && isnoisemixing(S.sensealg)
    # Diffusion term, diagonal noise, (as above but can handle mixing noise terms)
    idx2 = [(length(u)+1)*i-idx for i in 1:idx] # for diagonal indices of [end-idx+1:end,1:idx]

    dλ    = @view du[1:idx,1:idx]
    dgrad = @view du[idx+1:end-idx,1:idx]
    dy    = @view du[idx2]

  elseif typeof(du) <: AbstractMatrix
    # non-diagonal noise
    dλ    = @view du[1:idx, 1:idx]
    dgrad = @view du[idx+1:end-idx,1:idx]
    dy    = @view du[end-idx+1:end, 1:idx]
  end
  λ,grad,_y,dλ,dgrad,dy
end

# g is either g(t,u,p) or discrete g(t,u,i)
@noinline function ODEAdjointProblem(sol,sensealg::BacksolveAdjoint,
                                     g,t=nothing,dg=nothing;
                                     checkpoints=sol.t,
                                     callback=CallbackSet(),
                                     z0=nothing,
                                     M=nothing,
                                     nilss=nothing,
                                     tspan=sol.prob.tspan,
                                     kwargs...)
  # add homogenous adjoint for NILSAS by explicitly passing a z0 and nilss::NILSSSensitivityFunction
  @unpack f, p, u0 = sol.prob

  # check if solution was terminated, then use reduced time span
  terminated = false
  if hasfield(typeof(sol),:retcode)
    if sol.retcode == :Terminated
      tspan = (tspan[1], sol.t[end])
      terminated = true
    end
  end

  tspan = reverse(tspan)

  discrete = t !== nothing

  numstates = length(u0)
  numparams = p === nothing || p === DiffEqBase.NullParameters() ? 0 : length(p)

  len = length(u0)+numparams

  if z0===nothing
    λ = p === nothing || p === DiffEqBase.NullParameters() ? similar(u0) : one(eltype(u0)) .* similar(p, len)
    λ .= false
  else
    λ = nothing
  end

  sense = ODEBacksolveSensitivityFunction(g,sensealg,discrete,sol,dg,f)

  if z0!==nothing
    sense = NILSASSensitivityFunction{isinplace(f),typeof(nilss),typeof(sense),typeof(M)}(nilss,sense,M,discrete)
  end

  init_cb = t !== nothing && tspan[1] == t[end]
  cb, duplicate_iterator_times = generate_callbacks(sense, g, λ, t, callback, init_cb,terminated)
  checkpoints = ischeckpointing(sensealg, sol) ? checkpoints : nothing
  if checkpoints !== nothing
    cb = backsolve_checkpoint_callbacks(sense, sol, checkpoints, cb, duplicate_iterator_times)
  end

  if z0===nothing
    z0 = [vec(zero(λ)); vec(sense.y)]
  end
  original_mm = sol.prob.f.mass_matrix
  zzz(A, m, n) = fill!(similar(A, m, n), zero(eltype(original_mm)))
  if original_mm === I || original_mm === (I,I)
    mm = I
  else
    sense.diffcache.issemiexplicitdae && @warn "`BacksolveAdjoint` is likely to fail on semi-explicit DAEs, if memory is a concern, please consider using InterpolatingAdjoint(checkpoint=true) instead."
    II = Diagonal(I, numparams)
    Z1 = zzz(original_mm, numstates, numstates+numparams)
    Z2 = zzz(original_mm, numparams, numstates)
    mm = [copy(original_mm')   Z1
          Z2                   II  Z2
          Z1                       original_mm]
  end
  jac_prototype = sol.prob.f.jac_prototype
  if !sense.discrete || jac_prototype === nothing
    adjoint_jac_prototype = nothing
  else
    J = jac_prototype
    Ja = copy(J')
    II = Diagonal(I, numparams)
    Z1 = zzz(J, numstates, numstates+numparams)
    Z2 = zzz(J, numparams, numstates)
    adjoint_jac_prototype = [Ja       Z1
                             Z2       II     Z2
                             Z1              J]
  end
  odefun = ODEFunction(sense, mass_matrix=mm, jac_prototype=adjoint_jac_prototype)
  return ODEProblem(odefun,z0,tspan,p,callback=cb)
end

@noinline function SDEAdjointProblem(sol,sensealg::BacksolveAdjoint,
                                     g,t=nothing,dg=nothing;
                                     checkpoints=sol.t,
                                     callback=CallbackSet(),
                                     corfunc_analytical=nothing,diffusion_jac=nothing, diffusion_paramjac=nothing,
                                     kwargs...)
  @unpack f, p, u0, tspan = sol.prob
  # check if solution was terminated, then use reduced time span
  terminated = false
  if hasfield(typeof(sol),:retcode)
    if sol.retcode == :Terminated
      tspan = (tspan[1], sol.t[end])
      terminated = true
    end
  end

  tspan = reverse(tspan)
  discrete = t !== nothing

  p === DiffEqBase.NullParameters() && error("Your model does not have parameters, and thus it is impossible to calculate the derivative of the solution with respect to the parameters. Your model must have parameters to use parameter sensitivity calculations!")

  numstates = length(u0)
  numparams = length(p)

  len = length(u0)+numparams
  λ = one(eltype(u0)) .* similar(p, len)

  if StochasticDiffEq.alg_interpretation(sol.alg) == :Stratonovich
    sense_drift = ODEBacksolveSensitivityFunction(g,sensealg,discrete,sol,dg,sol.prob.f)
  else
    transformed_function = StochasticTransformedFunction(sol,sol.prob.f,sol.prob.g,corfunc_analytical)
    drift_function = ODEFunction(transformed_function)
    sense_drift = ODEBacksolveSensitivityFunction(g,sensealg,discrete,sol,dg,drift_function)
  end

  diffusion_function = ODEFunction(sol.prob.g, jac=diffusion_jac, paramjac=diffusion_paramjac)
  sense_diffusion = ODEBacksolveSensitivityFunction(g,sensealg,discrete,sol,dg,diffusion_function;noiseterm=true)

  init_cb = t !== nothing && tspan[1] == t[end]
  cb, duplicate_iterator_times = generate_callbacks(sense_drift, g, λ, t, callback, init_cb,terminated)
  checkpoints = ischeckpointing(sensealg, sol) ? checkpoints : nothing
  if checkpoints !== nothing
    cb = backsolve_checkpoint_callbacks(sense_drift, sol, checkpoints, cb, duplicate_iterator_times)
  end

  z0 = [vec(zero(λ)); vec(sense_drift.y)]

  original_mm = sol.prob.f.mass_matrix
  if original_mm === I
    mm = I
  else
    sense_drift.diffcache.issemiexplicitdae && @warn "`BacksolveAdjoint` is likely to fail on semi-explicit DAEs, if memory is a concern, please consider using InterpolatingAdjoint(checkpoint=true) instead."
    len2 = length(z0)
    mm = zeros(len2, len2)
    idx = 1:numstates
    copyto!(@view(mm[idx, idx]), sol.prob.f.mass_matrix')
    idx = numstates+1:numstates+1+numparams
    copyto!(@view(mm[idx, idx]), I)
    idx = len+1:len2
    copyto!(@view(mm[idx, idx]), sol.prob.f.mass_matrix)
  end

  sdefun = SDEFunction(sense_drift,sense_diffusion,mass_matrix=mm)

  # replicated noise
  _sol = deepcopy(sol)
  backwardnoise = reverse(_sol.W)

  if StochasticDiffEq.is_diagonal_noise(sol.prob) && typeof(sol.W[end])<:Number
    # scalar noise case
    noise_matrix = nothing
  else
    noise_matrix = similar(z0,length(z0),numstates)
    noise_matrix .= false
  end

  return SDEProblem(sdefun,sense_diffusion,z0,tspan,p,
    callback=cb,
    noise=backwardnoise,
    noise_rate_prototype = noise_matrix
    )
end

@noinline function RODEAdjointProblem(sol,sensealg::BacksolveAdjoint,
                                     g,t=nothing,dg=nothing;
                                     checkpoints=sol.t,
                                     callback=CallbackSet(),
                                     kwargs...)
  @unpack f, p, u0, tspan = sol.prob
  # check if solution was terminated, then use reduced time span
  terminated = false
  if hasfield(typeof(sol),:retcode)
    if sol.retcode == :Terminated
      tspan = (tspan[1], sol.t[end])
      terminated = true
    end
  end
  tspan = reverse(tspan)
  discrete = t !== nothing

  p === DiffEqBase.NullParameters() && error("Your model does not have parameters, and thus it is impossible to calculate the derivative of the solution with respect to the parameters. Your model must have parameters to use parameter sensitivity calculations!")

  numstates = length(u0)
  numparams = length(p)

  len = length(u0)+numparams
  λ = one(eltype(u0)) .* similar(p, len)

  sense = ODEBacksolveSensitivityFunction(g,sensealg,discrete,sol,dg,f;noiseterm=false)

  init_cb = t !== nothing && tspan[1] == t[end]
  cb, duplicate_iterator_times = generate_callbacks(sense, g, λ, t, callback, init_cb,terminated)
  checkpoints = ischeckpointing(sensealg, sol) ? checkpoints : nothing
  if checkpoints !== nothing
    cb = backsolve_checkpoint_callbacks(sense, sol, checkpoints, cb, duplicate_iterator_times)
  end

  z0 = [vec(zero(λ)); vec(sense.y)]

  original_mm = sol.prob.f.mass_matrix
  if original_mm === I
    mm = I
  else
    sense.diffcache.issemiexplicitdae && @warn "`BacksolveAdjoint` is likely to fail on semi-explicit DAEs, if memory is a concern, please consider using InterpolatingAdjoint(checkpoint=true) instead."
    len2 = length(z0)
    mm = zeros(len2, len2)
    idx = 1:numstates
    copyto!(@view(mm[idx, idx]), sol.prob.f.mass_matrix')
    idx = numstates+1:numstates+1+numparams
    copyto!(@view(mm[idx, idx]), I)
    idx = len+1:len2
    copyto!(@view(mm[idx, idx]), sol.prob.f.mass_matrix)
  end

  rodefun = RODEFunction(sense,mass_matrix=mm)

  # replicated noise
  _sol = deepcopy(sol)
  backwardnoise = reverse(_sol.W)

  if StochasticDiffEq.is_diagonal_noise(sol.prob) && typeof(sol.W[end])<:Number
    # scalar noise case
    noise_matrix = nothing
  else
    noise_matrix = similar(z0,length(z0),numstates)
    noise_matrix .= false
  end

  return RODEProblem(rodefun,z0,tspan,p,
    callback=cb,
    noise=backwardnoise,
    noise_rate_prototype = noise_matrix
    )
end

function backsolve_checkpoint_callbacks(sensefun, sol, checkpoints, callback, duplicate_iterator_times=nothing)
  prob = sol.prob
  if duplicate_iterator_times !== nothing
    _checkpoints = filter(x->x ∉ duplicate_iterator_times[1], checkpoints)
  else
    _checkpoints = checkpoints
  end
  cur_time = Ref(length(_checkpoints))
  affect! = let sol=sol, cur_time=cur_time, idx=length(prob.u0)
    function (integrator)
      _y = reshape(@view(integrator.u[end-idx+1:end]), axes(prob.u0))
      sol(_y, integrator.t)
      u_modified!(integrator,true)
      cur_time[] -= 1
      return nothing
    end
  end

  cb = PresetTimeCallback(_checkpoints,affect!)
  return CallbackSet(cb,callback)
end

function backsolve_checkpoint_callbacks(sensefun::NILSASSensitivityFunction, sol, checkpoints, callback, duplicate_iterator_times=nothing)
  prob = sol.prob
  if duplicate_iterator_times !== nothing
    _checkpoints = filter(x->x ∉ duplicate_iterator_times[1], checkpoints)
  else
    _checkpoints = checkpoints
  end
  cur_time = Ref(length(_checkpoints))
  affect! = let sol=sol, cur_time=cur_time
    function (integrator)
      _y = integrator.u.x[3]
      sol(_y, integrator.t)
      u_modified!(integrator,true)
      cur_time[] -= 1
      return nothing
    end
  end

  cb = PresetTimeCallback(_checkpoints,affect!)
  return CallbackSet(cb,callback)
end


text:
is_diagonal_noise(prob::DEProblem) = false
is_diagonal_noise(prob::AbstractRODEProblem{uType,tType,isinplace,ND}) where {uType,tType,isinplace,ND} = ND <: Nothing

isinplace(prob::AbstractODEProblem{uType,tType,iip}) where {uType,tType,iip} = iip
isinplace(prob::AbstractSteadyStateProblem{uType,iip}) where {uType,iip} = iip
isinplace(prob::AbstractRODEProblem{uType,tType,iip,ND}) where {uType,tType,iip,ND} = iip
isinplace(prob::AbstractDDEProblem{uType,tType,lType,iip}) where {uType,tType,lType,iip} = iip
isinplace(prob::AbstractDAEProblem{uType,duType,tType,iip}) where {uType,duType,tType,iip} = iip
isinplace(prob::AbstractNoiseProblem) = isinplace(prob.noise)
isinplace(::SplitFunction{iip}) where iip = iip


text:
"""
    transform!(node::Node, args..., <keyword arguments>)
    transform(node::Node, args..., <keyword arguments>)

Transform (mutate) an MTG (`node`) in place (`transform!`) or on a copy (`transform`) to add
attributes specified by `args...`.

# Arguments

- `node::Node`: An MTG node (*e.g.* the whole mtg returned by `read_mtg()`).
- `args::Any`: the transformations (see details)
- <keyword arguments>:

    - `scale = nothing`: The scale to filter-in (i.e. to keep). Usually a Tuple-alike of integers.
    - `symbol = nothing`: The symbol to filter-in. Usually a Tuple-alike of Strings.
    - `link = nothing`: The link with the previous node to filter-in. Usually a Tuple-alike of Char.
    - `filter_fun = nothing`: Any filtering function taking a node as input, e.g. [`isleaf`](@ref).
    - `ignore_nothing = false`: filter-out the nodes with `nothing` values for the given
    attributes used as inputs (apply only to the form :var_name => ...)

# Returns

`transform!`: Nothing, mutates the (sub-)tree in-place.
`transform`: A mutated copy of `node`.

# Notes

Carefull, `transform` is much slower than `transform!` because it makes a copy of the whole
MTG each time.

# Details

The interface of the function is inspired from the one used in
[`DataFrames.jl`](https://dataframes.juliadata.org/stable/), but adapted to an MTG.

The `args...` provided can be of the following forms:

1. a `:var_name => :new_var_name` pair. This form is used to rename an attribute name
2. a `:var_name => function => :new_var_name` or
`[:var_name1, :var_name2...] => function => :new_var_name` pair. The variables are declared
as a Symbol or a String (or a vector of), and they are passed as positional
arguments to the function. The new attribute name is optional and is automatically generated
if not provided by concatenating the source column name(s) and the function name if any.
3. a `function => :new_var_name` form that applies a function to a node and puts the results
in a new attribute. This form is usually applied when searching ancestors or descendants values.
4. a `function` form that applies a mutating function to a node, without expecting any output.
This form is adapted when using a function that already mutates the node, without the need to
return anything, *e.g.* [`branching_order!`](@ref).

Carefull to the form you use! Form 2 expect a function that takes one or more node
attributes (== variables) as inputs, while form 3 and 4 expect a function that takes a node.

# Examples

```julia
file = joinpath(dirname(dirname(pathof(MultiScaleTreeGraph))),"test","files","simple_plant.mtg")
mtg = read_mtg(file)

# We can use transform to apply a function over all nodes (same as using [`traverse!`](@ref))
transform!(mtg,  x -> isleaf(x) ? println(x.name," is a leaf") : nothing)
node_5 is a leaf
node_7 is a leaf

# We can compute a new variable based on another. For example to know if the value of the
# `:Length` attribute is provided or not, we can do:
transform!(mtg, :Length => isnothing)
# To check the values we first call [`get_attributes`](@ref) to know the new variable name:
get_attributes(mtg)
# And then we get the values using [`descendants`](@ref)
descendants(mtg, :Length_isnothing, self = true)
# Or DataFrame:
DataFrame(mtg, :Length_isnothing)

# We can also set the attribute name ourselves like so:
transform!(mtg, :Length => isnothing => :no_length)
descendants(mtg, :no_length, self = true)

# We can provide anonymous functions if we want to:
transform!(mtg, :Length => (x -> isnothing(x)) => :no_length)
descendants(mtg, :no_length, self = true)

# When a node does not have an attribute, it returns `nothing`. Most basic functions do not
# handle those very well, e.g.:
transform!(mtg, :Length => log)
# It does not work because some nodes have no value for `:Length`.
# To remove automatically the nodes with `nothing` values, use `ignore_nothing`:
transform!(mtg, :Length => log => :log_length, ignore_nothing = true)
descendants(mtg, :log_length, self = true)

# Or you could handle these manually in your function if you prefer:
transform!(mtg, :Length => (x -> x === nothing ? nothing : log(x)) => :log_length2)
descendants(mtg, :log_length2, self = true)

# Another way is to give a filtering function as an argument:
transform!(mtg, :Length => log => :log_length, filter_fun = x -> x[:Length] !== nothing)

# We can use more than one attribute as input to our function like so:
transform!(
    mtg,
    [:Width, :Length] => ((x, y) -> (x/2)^2 * π * y) => :volume,
    filter_fun = x -> x[:Length] !== nothing && x[:Width] !== nothing
)
descendants(mtg, :volume, self = true)

# Note that `filter_fun` filter the node, so we use the node[:attribute] notation here.

# We can also chain operations, and they will be executed sequentially so we can use variables
# computed on the instruction just before:
density = 0.6
transform!(
    mtg,
    [:Width, :Length] => ((x, y) -> (x/2)^2 * π * y) => :vol,
    :vol => (x -> x * density) => :biomass,
    filter_fun = x -> x[:Length] !== nothing && x[:Width] !== nothing
)
DataFrame(mtg, [:vol, :biomass])

# We can also rename a variable like so:
transform!(
    mtg,
    :biomass => :mass,
    filter_fun = x -> x[:Length] !== nothing && x[:Width] !== nothing
)
DataFrame(mtg, [:vol, :mass])

# Finally, we can use variables from ancestors/descendants using the `function => :new_var` form:
function get_mass_descendants(x)
    masses = descendants(x, :mass, ignore_nothing = true)
    if length(masses) == 0
        nothing
    else
        sum(masses)
    end
end

transform!(
    mtg,
    get_mass_descendants => :mass_beared
)
DataFrame(mtg, [:mass, :mass_beared])
```
"""
transform!, transform

function transform!(
    mtg::Node,
    args...;
    scale = nothing,
    symbol = nothing,
    link = nothing,
    filter_fun = nothing,
    ignore_nothing = false
)

    check_filters(mtg, scale = scale, symbol = symbol, link = link)
    filter_fun_ = filter_fun

    for nc in args
        if nc isa Symbol || nc isa String
            # The expression is just a variable name, we ignore it and pas sto the next expression
            continue
        elseif nc isa Base.Callable
            # If we provide just a function, it is applied to the node directly, so it must handle
            # a node as the first argument.
            # ?NOTE: Here the function takes a node as input
            fun_ = nc
        elseif first(nc) isa Symbol && last(nc) isa Symbol
            # `Name => new name` form, i.e. :x => :y.
            # ?NOTE: Here we just rename an attribute
            fun_ = x -> rename!(x, nc)
        elseif first(nc) isa Base.Callable && last(nc) isa Symbol
            # `function => new name` form, i.e. `node -> sum(descendants(node, :var)) => :newvar`.
            # ?NOTE: Here the function takes a node as input
            fun, newname = nc
            fun_ = x -> x[newname] = fun(x)
        elseif last(nc) isa Pair
            # `Name => function => new name` form, i.e. :x => sqrt => :x_sq
            # ?NOTE: Here the function takes one or more attributes as input
            col_idx, (fun, newname) = nc

            if !isa(col_idx, Vector) && !isa(col_idx, Tuple)
                col_idx = [col_idx]
            end

            fun_ = x -> x[newname] = fun([x[i] for i in col_idx]...)

            # Add a filter to the filtering function: checks if the node has the attributes
            # and if not, filter-out the node (in case ignore_nothing == true)
            filter_fun_ = filter_fun_nothing(filter_fun, ignore_nothing, col_idx)
        else
            # `Name => function` form, i.e. :x => sqrt
            # ?NOTE: Here the function takes one or more attributes as input
            col_idx, fun = nc

            newname, col_idx = col_name_from_call(col_idx, fun)

            fun_ = x -> x[newname] = fun([x[i] for i in col_idx]...)

            # Add a filter to the filtering function: checks if the node has the attributes
            # and if not, filter-out the node (in case ignore_nothing == true)
            filter_fun_ = filter_fun_nothing(filter_fun, ignore_nothing, col_idx)
        end

        traverse!(
            mtg,
            fun_;
            scale = scale,
            symbol = symbol,
            link = link,
            filter_fun = filter_fun_
        )
    end
end

function transform(
    mtg::Node,
    args...;
    scale = nothing,
    symbol = nothing,
    link = nothing,
    filter_fun = nothing,
    ignore_nothing = false
)

    new_mtg = deepcopy(mtg)

    transform!(
        new_mtg,
        args...;
        scale = scale,
        symbol = symbol,
        link = link,
        filter_fun = filter_fun,
        ignore_nothing = ignore_nothing
    )

    return new_mtg
end

function col_name_from_call(col_idx, fun)
    if !isa(col_idx, Vector)
        col_idx = [col_idx]
    end
    fnname = Symbol(fun)
    fnname_string = String(fnname)
    col_idx_name = join([String(i) for i in col_idx], "_")

    if startswith(fnname_string, '#')
        newname = Symbol(string(col_idx_name, "_function"))
    else
        newname = Symbol(join([col_idx_name, String(fnname)], "_"))
    end

    return newname, col_idx
end