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

using Markdown
using InteractiveUtils

# ╔═╡ d65d8190-7aa2-11eb-0e64-837b0676b7ca
begin
	using Pkg
	Pkg.activate(".")
	using Plots, LaTeXStrings, PrettyTables
	include("tools.jl")
end;

# ╔═╡ 81a9741a-798f-11eb-22fc-7fd2f87480dc
md"""
# MATH 522 Numerical Analysis

Christoph Ortner, ortner@math.ubc.ca, University of British Columbia

A course on approximation theory and with focus on computational aspects and applications.
"""

# ╔═╡ d28341ea-798f-11eb-3819-a3bc46d12e77
md"""
# § 0. Introduction and Motivation

## § 0.1 Composite Trapezoidal Rule

Consider two functions defined on $[-\pi, \pi]$,
```math
	f_1(x) = \frac{1}{1 + x^2}, \qquad
	f_2(x) = \frac{1}{1 + \sin^2(x)}
```
"""

# ╔═╡ 4127884c-7991-11eb-17f4-c76f7a96a556
begin
	f1(x) = 1 / (1 + x^2)
	f2(x) = 1 / (1 + sin(x)^2)
	plot(f1, -π, π, lw=3, label = L"f_1(x) = 1/(1+x^2)")
	plot!(f2, -π, π, lw = 3, label = L"f_2(x) = 1/(1+\sin^2(x))",
		  size = (550, 250), legend = :outertopright )
end

# ╔═╡ 779c0428-7990-11eb-0437-0d930f606312
md"""
We approximate the integral $I[f_j] := \int_{-\pi}^\pi f_j(x) \,dx$ with a quadrature rule, the [composite trapezoidal rule](https://en.wikipedia.org/wiki/Trapezoidal_rule)
```math
\begin{aligned}
	I[f] \approx I_N[f] := \sum_{n = -N+1}^N \frac{2\pi}{2N}
		\cdot \frac{f(x_{n-1})) + f(x_n)}{2}
\end{aligned}
```
where $x_n = 2\pi n / (2N) = \pi n/N, n = -N, \dots, N$ are the quadrature nodes.

If $f \in C^1$ then it is not too difficult to show that on each sub-interval $(x_{n-1}, x_n)$ of length ``h \approx 1/N`` approximating ``f`` with a piecewise affine function yields an ``O(h^2)`` error and therefore the total error is expected to also scale like ``h^2 \approx N^{-2}``, i.e., we expect that
```math
  |I[f] - I_N[f]| \lesssim N^{-2}
```
"""

# ╔═╡ 8de1aef6-7abf-11eb-186a-f98c6e94f5d1
md"""
To test this numerically we implement the quadrature rule rewritten as follows:
```math
	I_N[f] = \frac{\pi}{2N} \big(f(-\pi) + f(\pi)\big)
		+ \frac{\pi}{N} \sum_{n = -N+1}^{N-1} f(\pi n/N)
```
"""

# ╔═╡ 07ad43fe-7993-11eb-28d0-a3a4b0efcf12
@doc raw"""
`trapezoidal_rule(f, N)` : composite trapezoidal rule on [-π, π],
```math
I_N[f] = \frac{\pi}{2N} \big(f(-\pi) + f(\pi)\big)
		+ \frac{\pi}{N} \sum_{n = -N+1}^{N-1} f(\pi n/N)
```

* `f` : function defining the integrand
* `N` : number of integration nodes is ``2N+1``
* Output: value for composite trapezoidal rule
"""
trapezoidal_rule(f, N) =  (
	    0.5*π/N * (f(-π) + f(π))
	+       π/N * sum( f(n*π/N) for n = -N+1:N-1 )  );

# ╔═╡ dc5cad34-7992-11eb-3fc4-a54869ad9b9f
begin
	NN = 3:3:30       # number of quadrature points N means 2N+1 points
	I1 = 2 * atan(π)  # exact value of ∫ f₁
	I2 = √2 * π       # exact value of ∫ f₂
	I1N = trapezoidal_rule.(f1, NN)   # trapezoidal rule approximations
	I2N = trapezoidal_rule.(f2, NN)
	E1N = abs.(I1N .- I1)   # errors
	E2N = abs.(I2N .- I2)
end;

# ╔═╡ 21134062-7ab5-11eb-1cdc-856a378f3c7d
md"""##### Convergence Trapezoidal rule
$(
  ata_table( (NN, "``N``", "%d"),
	       (I1N, "``I_N[f_1]``", "%1.5f"),
		   (E1N, "``E_N[f_1]``", "%1.1e"),
	       (I2N, "``I_N[f_2]``", "%1.5f"),
		   (E2N, "``E_N[f_2]``", "%1.1e"),
	        ))
* ``I_N[f]`` : trapezoidal rule approximation with ``2N+1`` quadrature points
* ``E_N[f] := |I_N[f] - \int_{-\pi}^\pi f(x) \,dx |`` : error
"""

# ╔═╡ edf219ca-7abf-11eb-3d52-83c82513ab0d
md"""Plotting the error will give us a more qualitative view ..."""

# ╔═╡ fc2e003a-7abf-11eb-38c3-2f5db67014e3
begin
	Px = plot(NN, E1N, lw=2, m=:o, ms=4, label = L"E_N[f_1]", yaxis = :log10)
	plot!(NN, E2N.+1e-16, lw=2, m=:o, ms=4, label = L"E_N[f_2]",
	      xlabel = L"N")
	P1 = plot!(deepcopy(Px), NN[3:end], 0.04*NN[3:end].^(-2), lw=2, ls=:dash,
				c=:black, label = L"N^{-2}", ylims = [1e-6, 1e-1], xaxis = :log10)
	P2 = plot!(Px, NN[2:6], 0.1*exp.(- 2 * log(1 + sqrt(2)) * NN[2:6]), lw=2,
			   c=:black, ls = :dash, label = L"e^{- 2 \alpha N}",
			   legend = :right, ylims = [1e-16, 1e-1])
	# alpha = log(sqrt(2)+1)
	plot(P1, P2, size = (500, 300))
end

# ╔═╡ 1f961e64-7abf-11eb-0429-5d79965db4ca
md"""
An unexpected outcome? By the end of the first part of this course we should
be able to explain this result.
"""

# ╔═╡ b1d433f0-7ac2-11eb-06ca-5ffe4fbf5bff
md"""
## §0.2 What is Approximation Theory

[Wikipedia:](https://en.wikipedia.org/wiki/Approximation_theory) In mathematics, approximation theory is concerned with how functions can best be approximated with simpler functions, and with quantitatively characterizing the errors introduced thereby. Note that what is meant by best and simpler will depend on the application.

For the purpose of computational mathematics we should start by asking what operations are available to us on a computer: +, -, *, /. Everything else must be built from those. This means that the only functions we can implement immediately are polynomials and rational functions:
```math
	p_N(x) = a_0 + a_1 x + \dots + a_N x^N, \qquad r_{NM}(x) = \frac{p_N(x)}{q_M(x)}
```
We could (maybe should?) build this entire course based on polynomials. Instead, I decided to use trigonometric polynomials; more on this in the next notebook.

In any programming language, including Julia which we are using in this course, when you call mathematical functions such as
```julia
exp, cos, sin, acos, log, ...
```
you are in fact evaluating a rational approximant that approximates this function to within machine precision (typically ``\epsilon \approx 10^{-16}`` for 64bit floating point accuracy).
"""

# ╔═╡ e5952ed8-88bc-11eb-08fe-27778602caba
md"""
Beyond implementing special functions, why should we approximate a general function ``f`` by a polynomial? There are many reasons:
* ``f`` might be expensive to evaluate, hence replacing it with a cheap but accurate "surrogate" ``p_N`` would give us computationally efficient access to ``f``
* ``f`` might be unknown, but we know it solves some equation, e.g. a PDE, ``L[f] = 0``. We may be able to construct an approximate equation for a polynomial ``L_N[p_N] = 0`` and prove that solving this implies ``p_N \approx f``.
* ``f`` might be unknown but we have some observations (data) about it. We might then "fit" a polynomial ``p_N`` to this data in order to infer further information about ``f``.

What implicit or explicit assumptions are we making in these tasks? How should we optimize our approximation parameters (the polynomial degree ``N`` in this case). What can we say about the accuracy of approximation, i.e. how close is ``p_N`` to ``f``? How can we optimally sample ``f`` to obtain good approximations? These are the kind of questions that approximation theory and numerical analysis are concerned with.
"""

# ╔═╡ 2caa010c-7aa8-11eb-1554-9b8864597364
md"""
## §0.4 Resources

### Julia

* https://julialang.org
* https://juliaacademy.com
* https://juliadocs.github.io/Julia-Cheat-Sheet/

Although you won't need it for this course, I recommend VS Code for serious work with Julia (unless you are already committed to another very good editor such as Emacs, Vim, Sublime etc. Atom is still a good choice but most development has now moved to VS Code.

### Pluto.jl

* https://github.com/fonsp/Pluto.jl
* https://www.wias-berlin.de/people/fuhrmann/SciComp-WS2021/assets/nb01-first-contact-pluto.html
* https://computationalthinking.mit.edu/Spring21/

### Course Material

* https://github.com/cortner/ApxThyApp.git

This contains Pluto notebooks, older Jupyter notebooks (I may choose to switch back to Jupyter if I decide I don't like Pluto), and my lecture notes which are still under development but at least the basic theory part is now maturing nicely. The lecture notes also contain further references and exercise.
"""

# ╔═╡ Cell order:
# ╟─d65d8190-7aa2-11eb-0e64-837b0676b7ca
# ╟─81a9741a-798f-11eb-22fc-7fd2f87480dc
# ╟─d28341ea-798f-11eb-3819-a3bc46d12e77
# ╟─4127884c-7991-11eb-17f4-c76f7a96a556
# ╟─779c0428-7990-11eb-0437-0d930f606312
# ╟─8de1aef6-7abf-11eb-186a-f98c6e94f5d1
# ╠═07ad43fe-7993-11eb-28d0-a3a4b0efcf12
# ╟─dc5cad34-7992-11eb-3fc4-a54869ad9b9f
# ╟─21134062-7ab5-11eb-1cdc-856a378f3c7d
# ╟─edf219ca-7abf-11eb-3d52-83c82513ab0d
# ╟─fc2e003a-7abf-11eb-38c3-2f5db67014e3
# ╟─1f961e64-7abf-11eb-0429-5d79965db4ca
# ╟─b1d433f0-7ac2-11eb-06ca-5ffe4fbf5bff
# ╟─e5952ed8-88bc-11eb-08fe-27778602caba
# ╟─2caa010c-7aa8-11eb-1554-9b8864597364

max_forks_repo_name: cortner/ApxThyApp


text:
A = zeros(3, 4, 5);
A[:] = 1:60
Ap = @view(PermutedDimsArray(A,(3,1,2))[:,1:2,1])';

ap_index = ArrayInterface.StrideIndex(Ap)
for x_i in axes(Ap, 1)
    for y_i in axes(Ap, 2)
        @test ap_index[x_i, y_i] == ap_index[x_i, y_i]
    end
end
@test @inferred(ArrayInterface.known_offsets(ap_index)) === ArrayInterface.known_offsets(Ap)
@test @inferred(ArrayInterface.known_offset1(ap_index)) === ArrayInterface.known_offset1(Ap)
@test @inferred(ArrayInterface.offsets(ap_index, 1)) === ArrayInterface.offset1(Ap)
@test @inferred(ArrayInterface.offsets(ap_index, static(1))) === ArrayInterface.offset1(Ap)
@test @inferred(ArrayInterface.known_strides(ap_index)) === ArrayInterface.known_strides(Ap)
@test @inferred(ArrayInterface.contiguous_axis(ap_index)) == 1
@test @inferred(ArrayInterface.contiguous_axis(ArrayInterface.StrideIndex{2,(1,2),nothing,NTuple{2,Int},NTuple{2,Int}})) == nothing
@test @inferred(ArrayInterface.stride_rank(ap_index)) == (1, 3)

ext: jl


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder
name = "ReadStat"
version = v"1.1.5"

sources = [
    GitSource("https://github.com/WizardMac/ReadStat.git",
              "69f55186ae615a14a3367ad5cd08b7829aa8f308"),
]

# Bash recipe for building across all platforms
script = raw"""
# Add `gettext` for `autogen.sh`
apk add gettext-dev

# GCC builds complain about string truncation, but we don't care
if [[ ${target} != *apple* ]] && [[ ${target} != *freebsd* ]]; then
    export CFLAGS="${CFLAGS} -Wno-stringop-truncation"
fi

# Windows doesn't search ${prefix}/include?
export CPPFLAGS="${CPPFLAGS} -I${prefix}/include"

cd $WORKSPACE/srcdir/ReadStat/
./autogen.sh
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    ExecutableProduct("readstat", :readstat),
    LibraryProduct("libreadstat", :libreadstat),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Libiconv_jll"),
    Dependency("Zlib_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"8")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# SCILa - FakeNews Project
# 2021-05-25

# SETUP - always needed
using RCall, DataFrames, MixedModels, Serialization, LinearAlgebra, StatsModels

dat = rcopy(R"readRDS('./data/SCILa_FakeNews_1.rds')");
describe(dat)
levels(dat.UNL)
dat[:, [:unl1, :unl2, :unl3, :unl4]]

contr = merge(
	   Dict(nm => Grouping() for nm in (:Subj, :Item)),
	   Dict(nm => EffectsCoding() for nm in (:Frame, :Truth, :Source, :See, :Che)),
           Dict(:Gen => EffectsCoding(base="F")),
           Dict(:Edu => EffectsCoding(base="6")),
           Dict(:UNL => HypothesisCoding( [ -3 -3 +2 +2 +2   # unl1: Conflict vs. no conflict [x frame x truth]
                                            +2  0  0 -1 -1   # unl2: Conflict (UA vs KZ users) [x frame]
	         			    -2 +3 +3 -2 -2   # unl3: RU vs. (KZ+UA) [x  frame x loyalty]
	           			     0  0  0 +1 -1], # unl4: RU vs. UA language [x frame] | UA users
		                   levels=["kz ru ru", "ru kz ru", "ru ua ru", "ua ru ru", "ua ru ua"],
			           labels=["unl1", "unl2", "unl3", "unl4"])));
# Excursion: contr[:UNL]

## hypothesis matrix
using StatsModels, LinearAlgebra
cmat=
 [-3 -3 +2 +2 +2
  +2  0  0 -1 -1
  -2 +3 +3 -2 -2
   0  0  0 +1 -1]

## pseudoinverse matrix = contrast matrix
cmat_pi = pinv(cmat)

cmat_pi_v2 =
StatsModels.ContrastsMatrix(HypothesisCoding(cmat),
  ["kz_ru_ru", "ru_kz_ru", "ru_ua_ru", "ua_ru_ru", "ua_ru_ua"]).matrix

cmat_pi == cmat_pi_v2

## extract hypothesis matrix for constructed `ContrastMatrix`
cmat_v2 =
StatsModels.hypothesis_matrix(StatsModels.ContrastsMatrix(HypothesisCoding(cmat),
  ["kz_ru_ru", "ru_kz_ru", "ru_ua_ru", "ua_ru_ru", "ua_ru_ua"]))

# MODEL SELECTION - may be skipped
# only-varying-intercepts LMM:  m_ovi - SELECTED (selection with lme4)
f_ovi = @formula(answer ~ 1 + UNL*Frame*Truth + UNL*Frame*Source + UNL*Truth*Source + Frame*Truth*Source +
                              p_loy_c*UNL*Frame +  age_c+Gen+edu1+edu2+See+Che+
                              (UNL+Frame+Truth) & (age_c+Gen+edu1+edu2+See) +
                              (1  | Subj) +  (1  | Item));
m_ovi = fit(MixedModel, f_ovi, dat, contrasts=contr);

# full zero-correlation parameter LMM: m_zcp1
f_zcp1 = @formula(answer ~ 1 + UNL*Frame*Truth + UNL*Frame*Source + UNL*Truth*Source + Frame*Truth*Source +
                               p_loy_c*UNL*Frame + age_c+Gen+edu1+edu2+See+Che+
                              (UNL+Frame+Truth) & (age_c+Gen+edu1+edu2+See) +
                              zerocorr(1 + Frame + Truth + Source + See + Che | Subj) +
                              zerocorr(1 + unl4 + p_loy_c + Source + Gen + edu1 + edu2 + See + Che | Item));
m_zcp1 = fit(MixedModel, f_zcp1, dat, contrasts=contr);
issingular(m_zcp1)
VarCorr(m_zcp1)

# reduced zero-correlation parameter LMM: m_zcp - SELECTED
f_zcp =  @formula(answer ~ 1 + UNL*Frame*Truth + UNL*Frame*Source + UNL*Truth*Source + Frame*Truth*Source +
                               p_loy_c*UNL*Frame + age_c+Gen+edu1+edu2+See+Che+
                              (UNL+Frame+Truth) & (age_c+Gen+edu1+edu2+See) +
                              zerocorr(1 + Frame + Truth + See + Che | Subj) +
                              zerocorr(1 + unl4 + p_loy_c + Gen + edu1 + edu2 + See + Che | Item));
m_zcp = fit(MixedModel, f_zcp, dat, contrasts=contr);
issingular(m_zcp)
VarCorr(m_zcp)

MixedModels.likelihoodratiotest(m_zcp, m_zcp1)

# complex LMM  m_cpx_subj_all does not converge correctly
f_cpx_subj_all =  @formula(answer ~ 1 + UNL*Frame*Truth + UNL*Frame*Source + UNL*Truth*Source + Frame*Truth*Source +
                                    p_loy_c*UNL*Frame + age_c+Gen+edu1+edu2+See+Che+
                                   (UNL+Frame+Truth) & (age_c+Gen+edu1+edu2+See) +
                               (1 + Frame + Truth + See + Che | Subj) +
                       zerocorr(1 + unl4 + p_loy_c + Gen + edu1 + edu2 + See + Che | Item));
m_cpx_subj_all = fit(MixedModel, f_cpx_subj_all, dat, contrasts=contr);
issingular(m_cpx_subj_all) # true
VarCorr(m_cpx_subj_all)

# complex LMM with only a subset of Subj CPs -- SELECTED
f_cpx_subj =  @formula(answer ~ 1 + UNL*Frame*Truth + UNL*Frame*Source + UNL*Truth*Source + Frame*Truth*Source +
                                    p_loy_c*UNL*Frame + age_c+Gen+edu1+edu2+See+Che+
                                   (UNL+Frame+Truth) & (age_c+Gen+edu1+edu2+See) +
                               (1 | Subj) + (0 + Frame + Truth + See + Che | Subj) +
                       zerocorr(1 + unl4 + p_loy_c + Gen + edu1 + edu2 + See + Che | Item));
m_cpx_subj = fit(MixedModel, f_cpx_subj, dat, contrasts=contr);
issingular(m_cpx_subj) # false
VarCorr(m_cpx_subj)

# complex LMM with Item CPs added
f_cpx_subj_item =  @formula(answer ~ 1 + UNL*Frame*Truth + UNL*Frame*Source + UNL*Truth*Source + Frame*Truth*Source +
                                    p_loy_c*UNL*Frame + age_c+Gen+edu1+edu2+See+Che+
                                   (UNL+Frame+Truth) & (age_c+Gen+edu1+edu2+See) +
                               (1 | Subj) + (0 + Frame + Truth + See + Che | Subj) +
                               (1 + unl4 + p_loy_c + Gen + edu1 + edu2 + See + Che | Item));
m_cpx_subj_item = fit(MixedModel, f_cpx_subj_item, dat, contrasts=contr);
issingular(m_cpx_subj_item)  # true
VarCorr(m_cpx_subj_item)

# Compare LMMs: note m_cpx_subj_all did not converge correctly
MixedModels.likelihoodratiotest(m_ovi, m_zcp, m_cpx_subj, m_cpx_subj_all)
MixedModels.likelihoodratiotest(m_ovi, m_zcp, m_cpx_subj, m_cpx_subj_item)
#=
─────────────────────────────────────────────────────
     model-dof    -2 logLik        χ²  χ²-dof  P(>χ²)
─────────────────────────────────────────────────────
[1]         85  236452.3222
[2]         95  236136.5412  315.7810      10  <1e-61
[3]        102  236108.1477   28.3935       7  0.0002
[4]        130  236063.1876   44.9601      28  0.0223
─────────────────────────────────────────────────────
=#

lmms=[ "m_ovi", "m_zcp", "m_cpx_subj", "m_cpx_subj_item"];
mods = [m_ovi, m_zcp, m_cpx_subj, m_cpx_subj_item];
gof_summary = DataFrame(model=lmms,
                        dof=dof.(mods), deviance=deviance.(mods),
                        AIC = aic.(mods), BIC = bic.(mods))

#=
 Row │ model            dof    deviance   AIC        BIC
     │ String           Int64  Float64    Float64    Float64
─────┼─────────────────────────────────────────────────────────
   1 │ m_ovi               85  2.36452e5  2.36622e5  2.37383e5
   2 │ m_zcp               95  2.36137e5  2.36327e5  2.37177e5
   3 │ m_cpx_subj         102  2.36108e5  2.36312e5  2.37225e5
   4 │ m_cpx_subj_item    130  2.36063e5  2.36323e5  2.37487e5
=#

# Final model selection
m_cpx = m_cpx_subj;

#show(stdout, m1)
show(MIME("text/markdown"), m_cpx)  # show(MIME("text/xelatex"), m1)
VarCorr(m_cpx)

# write to file
open("1_SCILa_FakeNewsParams_MM.md", "w") do io
    show(io, MIME("text/markdown"), m_ovi)
    show(io, MIME("text/markdown"), m_zcp)
    show(io, MIME("text/markdown"), m_cpx)
end

# Store lmms
serialize("./fits/m_ovi.jls", m_ovi)
serialize("./fits/m_zcp.jls", m_zcp)
serialize("./fits/m_cpx.jls", m_cpx)

#= FOLLOW-UP ANALYSES
2021-05-28: 4-covariate interaction "unl3 &  frm  & trth & p_loy_c" not sign.
2021-05-29:
=#

m_ovi = deserialize("./fits/m_ovi.jls");
m_zcp = deserialize("./fits/m_zcp.jls");
m_cpx = deserialize("./fits/m_cpx.jls"); # m_cpx = m_cpx_subj_2 !

MixedModels.likelihoodratiotest(m_ovi, m_zcp, m_cpx, m_cpx_subj_item)

lmms=[ "m_ovi", "m_zcp", "m_cpx", "m_cpx_subj_item"];
mods = [m_ovi, m_zcp, m_cpx, m_cpx_subj_item];
gof_summary = DataFrame(model=lmms, dof=dof.(mods), deviance=deviance.(mods),
                        AIC = aic.(mods), BIC = bic.(mods))

#= Following up 2 three-variable interactions
  + unl1 x frame x poly_c
  + unl1 x frame x truth
 =#

## Establish equivalence with f_cpx_subj_2 (renamed above to f_cpx)
f_cpx_2 =  @formula(answer ~ 1 + (unl1+unl2+unl3+unl4)*frm*trth +
                                       (unl1+unl2+unl3+unl4)*frm*src +
                                       (unl1+unl2+unl3+unl4)*trth*src + frm*trth*src +
                               p_loy_c*(unl1+unl2+unl3+unl4)*frm + age_c+gen+edu1+edu2+see+che+
                                       (unl1+unl2+unl3+unl4+frm+trth) & (age_c+gen+edu1+edu2+see) +
                               (1 | Subj) + (0 + frm + trth + see + che | Subj) +
                       zerocorr(1 + unl4 + p_loy_c + gen + edu1 + edu2 + see + che | Item));
m_cpx_2 = fit(MixedModel, f_cpx_2, dat, contrasts=contr);
m_cpx_2
issingular(m_cpx_2) # false
VarCorr(m_cpx_2)
VarCorr(m_cpx)
deviance.([m_cpx, m_cpx_2])

## Add four-variable interaction
f_cpx_3 =  @formula(answer ~ 1 +      (unl1+unl2+unl3+unl4)*frm*trth +
                                      (unl1+unl2+unl3+unl4)*frm*src +
                                      (unl1+unl2+unl3+unl4)*trth*src + frm*trth*src +
                              p_loy_c*(unl1+unl2+unl3+unl4)*frm + age_c+gen+edu1+edu2+see+che+
                                      (unl1+unl2+unl3+unl4+frm+trth) & (age_c+gen+edu1+edu2+see) +
                                       unl3 & frm & trth & p_loy_c  +
                            (1 | Subj) + (0 + frm + trth + see + che | Subj)  +
                    zerocorr(1 + unl4 + p_loy_c + gen + edu1 + edu2 + see + che | Item));

m_cpx_3 = fit(MixedModel, f_cpx_3, dat, contrasts=contr);
m_cpx_3
issingular(m_cpx_3)
VarCorr(m_cpx_3)

MixedModels.likelihoodratiotest( m_cpx_2, m_cpx_3)

#=
───────────────────────────────────────────────────
     model-dof    -2 logLik      χ²  χ²-dof  P(>χ²)
───────────────────────────────────────────────────
[1]         98  236112.9616
[2]         99  236112.9613  0.0003       1  0.9869
───────────────────────────────────────────────────
=#

max_stars_repo_name: hse-scila/fakenews-project

max_issues_repo_name: hse-scila/fakenews-project

max_forks_repo_name: hse-scila/fakenews-project


text:
# Model
include_model("hopper3D")

# Horizon
T = 31

# Time step
tf = 3.0
h = tf / (T - 1)

# Bounds
_uu = Inf * ones(model.m)
_ul = zeros(model.m)
_ul[model.idx_u] .= -Inf
ul, uu = control_bounds(model, T, _ul, _uu)

# Initial and final states
mrp_init = MRP(UnitQuaternion(RotZ(0.0) * RotY(0.0) * RotX(0.0)))

z_h = 0.0
q1 = [0.0, 0.0, 0.5, mrp_init.x, mrp_init.y, mrp_init.z, 0.5]
x1 = [q1; q1]
qT = [1.0, 1.0, 0.5 + z_h, mrp_init.x, mrp_init.y, mrp_init.z, 0.5]
xT = [qT; qT]

xl, xu = state_bounds(model, T,
    [model.qL; model.qL], [model.qU; model.qU],
    x1 = x1, xT = xT)

# Objective
Qq = Diagonal([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0])
Q = cat(0.5 * Qq, 0.5 * Qq, dims = (1, 2))
QT = cat(0.5 * Qq, 100.0 * Diagonal(ones(model.nq)), dims = (1, 2))
R = Diagonal([1.0e-1, 1.0e-1, 1.0e-3, zeros(model.m - model.nu)...])

obj_tracking = quadratic_tracking_objective(
    [t < T ? Q : QT for t = 1:T],
    [R for t = 1:T-1],
    [xT for t = 1:T],
    [zeros(model.m) for t = 1:T])
obj_penalty = PenaltyObjective(1.0e5, model.m)
obj = MultiObjective([obj_tracking, obj_penalty])

# Constraints
include_constraints("contact")
con_contact = contact_constraints(model, T)

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con_contact)

# Trajectory initialization
x0 = [x1 for t = 1:T] #linear_interpolation(x1, x1, T) # linear interpolation on state
u0 = [1.0e-5 * rand(model.m) for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem

@time z̄ , info = solve(prob, copy(z0), tol = 1.0e-3, c_tol = 1.0e-3, mapl = 5)

check_slack(z̄, prob)
x̄, ū = unpack(z̄, prob)

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
open(vis)
visualize!(vis, model, state_to_configuration(x̄), Δt = h)


text:
function s(n)
    sum([parse(Int, x) for x in n])
end

function solve()
    l = split(readline())
    maximum([s(x) for x in l])
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
include("../mesh/mesh_clear.jl")

@testset "Simple values always same result" begin
    expected_result = 0.2
    @test fh(2, 3) == expected_result
    @test fh(0.2, 3) == expected_result
    @test fh(0.3, -0.1) == expected_result
end

@testset "Exception on wrong type of argument" begin
    @test_throws MethodError fh(1 + 1im, 3)
    @test_throws MethodError fh(1 + 2im, 3 + 3im)
    @test_throws MethodError fh(3 + 3im, 1 + 1im)
end

max_stars_repo_name: jstarczewski/julia-num

max_issues_repo_name: jstarczewski/julia-num

max_forks_repo_name: jstarczewski/julia-num


text: Dict{String,Any}("a" => Dict{String,Any}("b" => Dict{String,Any}()))

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
@testset "915.partition-array-into-disjoint-intervals.jl" begin
    @test partition_disjoint([5, 0, 3, 8, 6]) == 3
    @test partition_disjoint([1, 1, 1, 0, 6, 12]) == 4
    @test partition_disjoint([20, 34, 34, 17, 11, 31, 20, 10, 28, 31, 10,
        19, 9, 16, 23, 10, 38, 6, 40, 4, 38, 33, 40, 27, 48, 41, 48, 46, 44, 50]) == 24
end


text:
using CUDA

function runSimulation!(simulation)
    runSimulation!(simulation, TimerOutput())
end

function runSimulation!(simulation, to::TimerOutput)

    @synctimeit to "run simulation" begin

        # Unpack time data
        timeData = simulation.timeData
        maxsteps = timeData.maxsteps

        # Now start the time loop
        println("Starting time loop ...")

        kiter = 0
        for n = 1:maxsteps

            simulationConverged, Δh_max, max_index = doExplicitIter!(n, simulation, to::TimerOutput)
            kiter += 1

            if simulationConverged
                println()
                println("-- simulation converged --")
                println("Δh_max = ", Δh_max, " at ", max_index)
                println("nr. of iters = ", kiter)
                println()

                return true, kiter
            end
        end

        return false, kiter
    end # end timer
end

function doTimeStep!(simulation)
    doTimeStep!(simulation, TimerOutput())
end

function doTimeStep!(simulation, to)
    # Unpack data
    state       = simulation.state
    timeData    = simulation.timeData

    Δt          = timeData.Δt
    time        = timeData.time

    prepareAndSolve!(simulation, to)

    # Update the old to the new solution
    @synctimeit to "update state" begin
        updateState!(state)
        simulation.timeData.time += Δt
    end

    nothing
end

function doExplicitIter!(n, simulation)
    doExplicitIter!(n, simulation, TimerOutput())
end

function doExplicitIter!(n, simulation, to)

    # Unpack data
    state       = simulation.state
    timeData    = simulation.timeData
    solverData  = simulation.solverData

    Δt          = timeData.Δt
    time        = timeData.time

    prepareAndSolve!(simulation, to)

    # check convergence
    @synctimeit to "convergence check" begin
        simulationConverged, Δh_max, max_index = checkConvergence!(state, solverData)
    end

    @synctimeit to "print iterations" begin
        if mod(n, 100) == 0
            println("iter = ", n, ", Δh_max = ", Δh_max, " at ", max_index)
        end
    end

    # Update the old to the new solution
    @synctimeit to "update state" begin
        updateState!(state)
        time += Δt
    end

    return simulationConverged, Δh_max, max_index
end

function prepareAndSolve!(simulation, to)
        # Unpack data
        grid        = simulation.grid
        model       = simulation.model
        state       = simulation.state
        parameters  = simulation.parameters
        timeData    = simulation.timeData
        time        = timeData.time

        # Add the sources
        @synctimeit to "set rhs" begin
            setSources!(time, model.source)
            setRecharge!(time, model.recharge)
        end

        @synctimeit to "set_boundaries" begin
        bc = model.boundaryConditions
            setBoundaries!(grid, state, bc)
        end

            @synctimeit to "solve pressure" begin
            pressureEquation!(grid, model, state, parameters, timeData)
        end

        @synctimeit to "solve darcy" begin
            darcyEquation!(grid, model, state, parameters, timeData)
        end
end

function updateState!(state)
    state.h.parent .= state.hⁿ⁺¹.parent
    state.u.parent .= state.uⁿ⁺¹.parent
    state.v.parent .= state.vⁿ⁺¹.parent
    state.w.parent .= state.wⁿ⁺¹.parent
end

function checkConvergence!(state, solverData)
    solverData.Δh.parent .= abs.(state.hⁿ⁺¹.parent - state.h.parent)
    Δh_max, max_index = findMaximum(solverData.Δh.parent)
    if Δh_max > solverData.ΔhConv
        return false, Δh_max, max_index
    end

    # convergence
    return true, Δh_max, max_index
end

function findMaximum(h)
    return findmax(h)
end

function findMaximum(h::CuArray)
    m = reduce(max, h)
    return m, -1
end


text:
abstract type SpinSystem end

# Spin system

struct XYZSystem{T} <: SpinSystem
    J::Matrix{T}
    h::VecOrMat{T}
    anisotropy::Vector{T}
end

get_dt(system::XYZSystem) = 1/maximum(maximum(system.J, dims=1))/4

function build_hamiltonian!(system::XYZSystem)
    J, h, anisotropy = system.J, system.h, system.anisotropy
    (ds, s, p, t) -> hamiltonian!(ds, s, J, h, anisotropy, t)
end

# Spin system with finite temperature

struct XYZSystem_t <: SpinSystem
    J_t
    h_t
    anisotropy
end

function build_hamiltonian!(system::XYZSystem_t)
    J_t, h_t, anisotropy = system.J_t, system.h_t, system.anisotropy
    (ds, s, p, t) -> hamiltonian!(ds, s, J_t(t), h_t(t), anisotropy, t)
end

get_dt(system::XYZSystem_t) = 1/maximum(maximum(system.J_t(0.0), dims=1))/4

function reduction(u, data, i)
    (u .+ sum(data), false)
end

function dtwa_prob_func(prob, Ψ)
    prob.u0 .= permutedims(dtwa_state(Ψ))
    return prob
end

function conserve_spin_norm(resid,u,p,t)
    for i in 1:size(u, 1)
        u_i = @view u[i, :]
        resid[i, 1] = norm(u_i) - sqrt(3.)
    end
    resid[:, 2:end] .= 0
end

conserve_spin_callback = ManifoldProjection(conserve_spin_norm)

function DTWAProblem(system, Ψ, saveat; kwargs...)
    dt = get(kwargs, :dt, get_dt(system))
    hamiltonian! = build_hamiltonian!(system)
    Ψ_ini = permutedims(dtwa_state(Ψ))
    prob = ODEProblem(hamiltonian!, Ψ_ini, (saveat[1], saveat[end]),
        dense=false,
        dt=dt,
        force_dtmin=true,
        adaptive=true,
        saveat=saveat
    )

    EnsembleProblem(
        prob,
        prob_func = (prob, i, repeat) -> dtwa_prob_func(prob, Ψ),
        output_func = (sol, i) -> (vcat(mean.(sol.u, dims=1)...), false),
        reduction = reduction,
        u_init = zeros(Float64, (length(saveat), 3)),
    )
end

function DTWAProblem_noreduction(system::SpinSystem, Ψ, saveat; kwargs...)
    dt = get(kwargs, :dt, get_dt(system))
    hamiltonian! = build_hamiltonian!(system)
    Ψ_ini = permutedims(dtwa_state(Ψ))
    prob = ODEProblem(hamiltonian!, Ψ_ini, (saveat[1], saveat[end]),
        dense=false,
        dt=dt,
        force_dtmin=true,
        adaptive=true,
        saveat=saveat
    )

    EnsembleProblem(
        prob,
        prob_func = (prob, i, repeat) -> dtwa_prob_func(prob, Ψ),
        output_func = (sol,i) -> (sol.u, false),
        reduction = (u, data, I)->(append!(u, data), false),
        u_init = []
    )
end

function solve_dtwa(ensemble_prob, trajectories; kwargs...)
    tspan = ensemble_prob.prob.tspan
    nsteps = size(ensemble_prob.u_init, 1)
    sol = solve(
        ensemble_prob,
        Tsit5(),
        # AutoTsit5(Rosenbrock23()),
        # EnsembleSerial(),
        # EnsembleThreads(),
        # EnsembleDistributed(),
        trajectories=trajectories;
        #saveat=collect(LinRange(tspan[1], tspan[2], nsteps));
        kwargs...
    )
    sol.u ./= trajectories
end

function solve_dtwa(ensemble_prob, tracectories, s::Sequence; kwargs...)
    cb = IterativeCallback(s)
    solve_dtwa(ensemble_prob, tracectories; callback=cb, kwargs...)
end

max_stars_repo_name: titusfranz/DTWA

max_issues_repo_name: titusfranz/DTWA

max_forks_repo_name: titusfranz/DTWA


text:
using Distributions, JLD, CSV, DataFrames,Pkg, DataStructures

austin_stations_pre = CSV.File("raw/austin-hospitals.csv") |> DataFrame

# Yeesian has his coords in 3-stations.csv in this string format:
# "GeoInterface.Point([-77.00514742400335,38.83089499410552])"
formatted_coordinates = [string("GeoInterface.Point([", austin_stations_pre["longitude"][i], ", ", austin_stations_pre["latitude"][i], "])") for i in 1:size(austin_stations_pre["longitude"], 1)]

austin_stations_post = DataFrame(geometry = formatted_coordinates,
                                 name = austin_stations_pre["hospital_names"]
                                 )

CSV.write("formatted/austin_hospitals_formatted.csv", austin_stations_post)

max_stars_repo_name: michaelhilborn/AmbulanceDeployement

max_issues_repo_name: joshua-ong/AmbulanceDeployment

max_forks_repo_name: joshua-ong/AmbulanceDeployment


text:
getmark(x::T) where {T<:Base.IEEEFloat} = Marker0

getmark(x::T) where {T<:MarkableFloat{0}} = Marker0
getmark(x::T) where {T<:MarkableFloat{1}} = Marker1
getmark(x::T) where {T<:MarkableFloat{2}} = Marker2
getmark(x::T) where {T<:MarkableFloat{3}} = Marker3

getmarkvalue(x::T) where {T<:Base.IEEEFloat} = 0

getmarkvalue(x::T) where {T<:MarkableFloat{0}} = 0
getmarkvalue(x::T) where {T<:MarkableFloat{1}} = 1
getmarkvalue(x::T) where {T<:MarkableFloat{2}} = 2
getmarkvalue(x::T) where {T<:MarkableFloat{3}} = 3

@inline setmarkvalue(x::T, value::Int) where {T<:Union{Base.IEEEFloat, MarkableFloat}} =
    value in 0:3 ? setmark(x, Marker{value}) : throw(ErrorException("marker value out of range ($value)"))

setmark(x::Float64, ::Type{Marker0}) = x
setmark(x::Float32, ::Type{Marker0}) = x
setmark(x::Float16, ::Type{Marker0}) = x

setmark(x::Float64, ::Type{Marker1}) = MarkableF641(x)
setmark(x::Float32, ::Type{Marker1}) = MarkableF321(x)
setmark(x::Float16, ::Type{Marker1}) = MarkableF161(x)

setmark(x::Float64, ::Type{Marker2}) = MarkableF642(x)
setmark(x::Float32, ::Type{Marker2}) = MarkableF322(x)
setmark(x::Float16, ::Type{Marker2}) = MarkableF162(x)

setmark(x::Float64, ::Type{Marker3}) = MarkableF643(x)
setmark(x::Float32, ::Type{Marker3}) = MarkableF323(x)
setmark(x::Float16, ::Type{Marker3}) = MarkableF163(x)

setmark(x::MarkableF640, ::Type{Marker0}) = x
setmark(x::MarkableF320, ::Type{Marker0}) = x
setmark(x::MarkableF160, ::Type{Marker0}) = x

setmark(x::MarkableF640, ::Type{Marker1}) = MarkableF641(x)
setmark(x::MarkableF320, ::Type{Marker1}) = MarkableF321(x)
setmark(x::MarkableF160, ::Type{Marker1}) = MarkableF161(x)

setmark(x::MarkableF640, ::Type{Marker2}) = MarkableF642(x)
setmark(x::MarkableF320, ::Type{Marker2}) = MarkableF322(x)
setmark(x::MarkableF160, ::Type{Marker2}) = MarkableF162(x)

setmark(x::MarkableF640, ::Type{Marker3}) = MarkableF643(x)
setmark(x::MarkableF320, ::Type{Marker3}) = MarkableF323(x)
setmark(x::MarkableF160, ::Type{Marker3}) = MarkableF163(x)

setmark(x::MarkableF641, ::Type{Marker0}) = MarkableF640(x)
setmark(x::MarkableF321, ::Type{Marker0}) = MarkableF320(x)
setmark(x::MarkableF161, ::Type{Marker0}) = MarkableF160(x)

setmark(x::MarkableF641, ::Type{Marker1}) = x
setmark(x::MarkableF321, ::Type{Marker1}) = x
setmark(x::MarkableF161, ::Type{Marker1}) = x

setmark(x::MarkableF641, ::Type{Marker2}) = MarkableF642(x)
setmark(x::MarkableF321, ::Type{Marker2}) = MarkableF322(x)
setmark(x::MarkableF161, ::Type{Marker2}) = MarkableF162(x)

setmark(x::MarkableF641, ::Type{Marker3}) = MarkableF643(x)
setmark(x::MarkableF321, ::Type{Marker3}) = MarkableF323(x)
setmark(x::MarkableF161, ::Type{Marker3}) = MarkableF163(x)

setmark(x::MarkableF642, ::Type{Marker0}) = MarkableF640(x)
setmark(x::MarkableF322, ::Type{Marker0}) = MarkableF320(x)
setmark(x::MarkableF162, ::Type{Marker0}) = MarkableF160(x)

setmark(x::MarkableF642, ::Type{Marker1}) = MarkableF641(x)
setmark(x::MarkableF322, ::Type{Marker1}) = MarkableF321(x)
setmark(x::MarkableF162, ::Type{Marker1}) = MarkableF161(x)

setmark(x::MarkableF642, ::Type{Marker2}) = x
setmark(x::MarkableF322, ::Type{Marker2}) = x
setmark(x::MarkableF162, ::Type{Marker2}) = x

setmark(x::MarkableF642, ::Type{Marker3}) = MarkableF643(x)
setmark(x::MarkableF322, ::Type{Marker3}) = MarkableF323(x)
setmark(x::MarkableF162, ::Type{Marker3}) = MarkableF163(x)

setmark(x::MarkableF643, ::Type{Marker0}) = MarkableF640(x)
setmark(x::MarkableF323, ::Type{Marker0}) = MarkableF320(x)
setmark(x::MarkableF163, ::Type{Marker0}) = MarkableF160(x)

setmark(x::MarkableF643, ::Type{Marker1}) = MarkableF641(x)
setmark(x::MarkableF323, ::Type{Marker1}) = MarkableF321(x)
setmark(x::MarkableF163, ::Type{Marker1}) = MarkableF161(x)

setmark(x::MarkableF643, ::Type{Marker2}) = MarkableF642(x)
setmark(x::MarkableF323, ::Type{Marker2}) = MarkableF322(x)
setmark(x::MarkableF163, ::Type{Marker2}) = MarkableF162(x)

setmark(x::MarkableF643, ::Type{Marker3}) = x
setmark(x::MarkableF323, ::Type{Marker3}) = x
setmark(x::MarkableF163, ::Type{Marker3}) = x

unmark(x::Float64) = x
unmark(x::Float32) = x
unmark(x::Float16) = x

unmark(x::MarkableF640) = x
unmark(x::MarkableF320) = x
unmark(x::MarkableF160) = x

unmark(x::MarkableF641) = Float64(x)
unmark(x::MarkableF321) = Float32(x)
unmark(x::MarkableF161) = Float16(x)

unmark(x::MarkableF642) = Float64(x)
unmark(x::MarkableF322) = Float32(x)
unmark(x::MarkableF162) = Float16(x)

unmark(x::MarkableF643) = Float64(x)
unmark(x::MarkableF323) = Float32(x)
unmark(x::MarkableF163) = Float16(x)


text:
using Test

tests = [
    "format_checksum",
    "citation_generation",
    "UCI",
    "GitHub",
    "DataDryad",
    "DataOneV2/KNB",
    "DataOneV2/TERN",
    "CKAN",
    "DataCite",
    "Figshare",
    "JSONLD",
    "combine",
]

@testset "DataDepGenerators" begin
    for filename in tests
        @testset "$filename" begin
            include(filename * ".jl")
        end
    end
end


text:
function precondition(obj::Preconditioner, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(obj, "precondition", Vector{jdouble}, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end


text:
@testset "Ranges" begin
    @testset "UnitRange" begin
        for T in (halfinttypes..., halfuinttypes..., :BigHalfInt)
            @eval @test $T(1/2):$T(5) isa UnitRange{$T}
            @eval @test $T(1/2):big(5) isa UnitRange{BigHalfInt}
            @eval @test first($T(1/2):$T(5)) == $T(1/2)
            @eval @test last($T(1/2):$T(5)) == $T(9/2)
            @eval @test length($T(1/2):$T(5)) == 5
            @eval @test length($T(1/2):$T(0)) == 0
            @eval @test length($T(1/2):$T(1/2)) == 1
            @eval @test length($T(1/2):$T(1)) == 1
            @eval @test $T(5/2) ∉ $T(1):$T(5)
            @eval @test $T(5/2) ∈ $T(1/2):$T(5)
            @eval @test $T(2) ∈ $T(1):$T(5)
            @eval @test $T(2) ∉ $T(1/2):$T(5)
            @eval @test 2.5 ∉ $T(1):$T(5)
            @eval @test 2.5 ∈ $T(1/2):$T(5)
            @eval @test 5//2 ∉ $T(1):$T(5)
            @eval @test 5//2 ∈ $T(1/2):$T(5)
            @eval @test 2 ∈ $T(1):$T(5)
            @eval @test 2 ∉ $T(1/2):$T(5)

            @eval @test any(isinteger, $T(1):$T(5))
            @eval @test !any(isinteger, $T(1):$T(0))
            @eval @test !any(isinteger, $T(1/2):$T(5))
            @eval @test !any(isinteger, $T(3/2):$T(0))
            @eval @test !any(!isinteger, $T(1):$T(5))
            @eval @test !any(!isinteger, $T(1):$T(0))
            @eval @test any(!isinteger, $T(1/2):$T(5))
            @eval @test !any(!isinteger, $T(3/2):$T(0))

            @eval @test all(isinteger, $T(1):$T(5))
            @eval @test all(isinteger, $T(1):$T(0))
            @eval @test !all(isinteger, $T(1/2):$T(5))
            @eval @test all(isinteger, $T(3/2):$T(0))
            @eval @test !all(!isinteger, $T(1):$T(5))
            @eval @test all(!isinteger, $T(1):$T(0))
            @eval @test all(!isinteger, $T(1/2):$T(5))
            @eval @test all(!isinteger, $T(3/2):$T(0))
        end
        for T in (halfinttypes..., :BigHalfInt)
            @eval @test first($T(-1/2):$T(1)) == $T(-1/2)
            @eval @test last($T(-1/2):$T(1)) == $T(1/2)
            @eval @test length($T(-1/2):$T(1)) == 2
        end
    end

    @testset "StepRange" begin
        for T in (halfinttypes..., halfuinttypes..., :BigHalfInt)
            @eval @test $T(2):$T(1/2):$T(5) isa StepRange{$T}
            @eval @test $T(2):$T(2):$T(9/2) isa StepRange{$T}
            @eval @test $T(2):big(2):$T(9/2) isa StepRange{BigHalfInt}
            @eval @test big(2):2:$T(9/2) isa StepRange{BigHalfInt}

            @eval @test $T(2):$T(2):$T(9/2) == $T(2):$T(2):$T(4)
            @eval @test first($T(2):$T(2):$T(9/2)) == 2
            @eval @test step($T(2):$T(2):$T(9/2)) == 2
            @eval @test last($T(2):$T(2):$T(9/2)) == 4
            @eval @test length($T(2):$T(2):$T(9/2)) == 2

            @eval @test first($T(2):$T(1/2):$T(5)) == 2
            @eval @test step($T(2):$T(1/2):$T(5)) == $T(1/2)
            @eval @test last($T(2):$T(1/2):$T(5)) == 5
            @eval @test length($T(2):$T(1/2):$T(5)) == 7

            @eval @test $T(5/2) ∈ $T(1/2):$T(2):$T(7)
            @eval @test $T(5/2) ∉ $T(1/2):$T(5/2):$T(7)
            @eval @test $T(2) ∈ $T(0):$T(1):$T(5/2)
            @eval @test $T(2) ∉ $T(1/2):$T(1):$T(5/2)
            @eval @test 2.5 ∈ $T(1/2):$T(2):$T(7)
            @eval @test 2.5 ∉ $T(1/2):$T(5/2):$T(7)
            @eval @test 5//2 ∈ $T(1/2):$T(2):$T(7)
            @eval @test 5//2 ∈ $T(1/2):2:$T(7)
            @eval @test 5//2 ∈ $T(1/2):2:7
            @eval @test 5//2 ∉ $T(1/2):$T(5/2):$T(7)
            @eval @test 2 ∈ $T(0):$T(1):$T(5/2)
            @eval @test 2 ∉ $T(1/2):$T(1):$T(5/2)
        end
        for T in (halfinttypes..., :BigHalfInt)
            @eval @test $T(2):$T(-1):$T(1/2) == $T(2):$T(-1):$T(1)
            @eval @test first($T(2):$T(-1):$T(1/2)) == 2
            @eval @test step($T(2):$T(-1):$T(1/2)) == -1
            @eval @test last($T(2):$T(-1):$T(1/2)) == 1
            @eval @test length($T(2):$T(-1):$T(1/2)) == 2

            @eval @test $T(5/2) ∉ $T(1/2):$T(-1):$T(-1/2)
            @eval @test 2.5 ∉ $T(1/2):$T(-1):$T(-1/2)
            @eval @test 5//2 ∉ $T(1/2):$T(-1):$T(-1/2)
            @eval @test 5//2 ∉ $T(1/2):-1:$T(-1/2)
        end
    end

    @testset "reverse" begin
        for T in (halfinttypes..., :BigHalfInt)
            # There is currently a discussion on how to treat Unsigned types:
            # * https://github.com/JuliaLang/julia/issues/29576
            # * https://github.com/JuliaLang/julia/pull/29842
            @eval @test @inferred(reverse($T(1/2):$T(9/2))) isa StepRange{$T}
            @eval @test @inferred(reverse($T(9/2):$T(-1):$T(1/2))) isa StepRange{$T}
            @eval @test reverse($T(1/2):$T(9/2)) == $T(9/2):$T(-1):$T(1/2)
            @eval @test reverse($T(9/2):$T(-1):$T(1/2)) == $T(1/2):$T(1):$T(9/2)
            @eval @test reverse($T(1/2):$T(1/2):$T(3)) == $T(3):$T(-1/2):$T(1/2)
            @eval @test reverse($T(-1/2):$T(2):$T(3/2)) == $T(3/2):$T(-2):$T(-1/2)
            @eval @test isempty(reverse($T(3/2):$T(1/2)))
            @eval @test isempty(reverse($T(3/2):$T(1):$T(1/2)))
            @eval @test isempty(reverse($T(3/2):$T(1/2):$T(1)))
            @eval @test isempty(reverse($T(7/2):$T(-2):$T(11/2)))
        end
    end

    @testset "intersect" begin
        for T in (halfinttypes..., halfuinttypes..., :BigHalfInt)
            @eval @test @inferred(($T(1):$T(3)) ∩ ($T(1/2):$T(3))) isa UnitRange{$T}
            @eval @test @inferred((big(1):big(3)) ∩ ($T(1/2):$T(3))) isa UnitRange{BigHalfInt}
            @eval @test @inferred(($T(1/2):$T(3)) ∩ (big(1):big(3))) isa UnitRange{BigHalfInt}
            @eval @test @inferred(($T(3/2):$T(1):$T(5)) ∩ ($T(1):$T(3))) isa StepRange{$T}
            @eval @test @inferred(($T(3/2):$T(1):$T(5)) ∩ (big(1):big(3))) isa StepRange{BigHalfInt}
            @eval @test @inferred((big(1):big(3)) ∩ ($T(3/2):$T(1):$T(5))) isa StepRange{BigHalfInt}
            @eval @test @inferred((big(1):big(3):big(7)) ∩ ($T(3/2):$T(5))) isa StepRange{BigHalfInt}
            @eval @test @inferred(($T(3/2):$T(5)) ∩ (big(1):big(3):big(7))) isa StepRange{BigHalfInt}
            @eval @test @inferred(($T(3/2):$T(1/2):$T(5)) ∩ ($T(1):$T(3):$T(5))) isa StepRange{$T}
            @eval @test @inferred((big(-1):big(3):big(5)) ∩ ($T(3/2):$T(1/2):$T(5))) isa StepRange{BigHalfInt}
            @eval @test @inferred(($T(3/2):$T(1/2):$T(5)) ∩ (big(-1):big(3):big(5))) isa StepRange{BigHalfInt}
            @eval @test isempty(($T(1):$T(3)) ∩ ($T(1/2):$T(3)))
            @eval @test ($T(1/2):$T(2)) ∩ ($T(3/2):$T(3)) == $T(3/2):$T(3/2)
            @eval @test isempty(($T(1/2):$T(2)) ∩ (1:2))
            @eval @test isempty((1:2) ∩ ($T(1/2):$T(2)))
            @eval @test ($T(2):$T(5)) ∩ (1:3) == 2:3
            @eval @test (1:3) ∩ ($T(2):$T(5)) == 2:3
            @eval @test ($T(2):$T(1):$T(5)) ∩ (1:3) == 2:1:3
            @eval @test (1:3) ∩ ($T(2):$T(1):$T(5)) == 2:1:3
            @eval @test isempty(($T(3/2):$T(1):$T(5)) ∩ (1:3))
            @eval @test isempty((1:3) ∩ ($T(3/2):$T(1):$T(5)))
            @eval @test ($T(2):$T(1/2):$T(5)) ∩ (1:3) == 2:1:3
            @eval @test (1:3) ∩ ($T(2):$T(1/2):$T(5)) == 2:1:3
            @eval @test ($T(3/2):$T(1/2):$T(5)) ∩ (1:3) == 2:1:3
            @eval @test (1:3) ∩ ($T(3/2):$T(1/2):$T(5)) == 2:1:3
            @eval @test ($T(2):$T(1):$T(5)) ∩ ($T(1):$T(3)) == $T(2):$T(1):$T(3)
            @eval @test ($T(1):$T(3)) ∩ ($T(2):$T(1):$T(5)) == $T(2):$T(1):$T(3)
            @eval @test isempty(($T(3/2):$T(1):$T(5)) ∩ ($T(1):$T(3)))
            @eval @test isempty(($T(1):$T(3)) ∩ ($T(3/2):$T(1):$T(5)))
            @eval @test ($T(2):$T(1/2):$T(5)) ∩ ($T(1):$T(3)) == $T(2):$T(1):$T(3)
            @eval @test ($T(1):$T(3)) ∩ ($T(2):$T(1/2):$T(5)) == $T(2):$T(1):$T(3)
            @eval @test ($T(3/2):$T(1/2):$T(5)) ∩ ($T(1):$T(3)) == $T(2):$T(1):$T(3)
            @eval @test ($T(1):$T(3)) ∩ ($T(3/2):$T(1/2):$T(5)) == $T(2):$T(1):$T(3)
            @eval @test ($T(3/2):$T(1):$T(6)) ∩ ($T(1):$T(3/2):$T(6)) == $T(5/2):$T(3):$T(11/2)
        end
        for T in (halfinttypes..., :BigHalfInt)
            @eval @test ($T(3/2):$T(1/2):$T(5)) ∩ (-1:3:5) == 2:3:5
            @eval @test (-1:3:5) ∩ ($T(3/2):$T(1/2):$T(5)) == 2:3:5
            @eval @test ($T(2):$T(1):$T(5)) ∩ (-1:3:5) == 2:3:5
            @eval @test (-1:3:5) ∩ ($T(2):$T(1):$T(5)) == 2:3:5
            @eval @test isempty(($T(3/2):$T(1):$T(5)) ∩ (-1:3:5))
            @eval @test isempty((-1:3:5) ∩ ($T(3/2):$T(1):$T(5)))
            @eval @test isempty(($T(3/2):$T(1):$T(5)) ∩ ($T(-1):$T(3):$T(5)))
            @eval @test isempty(($T(-1):$T(3):$T(5)) ∩ ($T(3/2):$T(1):$T(5)))
            @eval @test ($T(3/2):$T(1):$T(5)) ∩ ($T(-1):$T(3/2):$T(5)) == $T(7/2):$T(3):$T(7/2)
            # These tests should in theory also work for Unsigned types, but due to some
            # bugs with ranges of Unsigned integers, they currently don’t (see below)
            @eval @test ($T(5/2):$T(3/2):$T(7)) ∩ ($T(1):$T(8)) == $T(4):$T(3):$T(7)
            @eval @test ($T(1):$T(8)) ∩ ($T(5/2):$T(3/2):$T(7)) == $T(4):$T(3):$T(7)
            @eval @test isempty((1:3:7) ∩ ($T(3/2):$T(5)))
            @eval @test isempty(($T(3/2):$T(5)) ∩ (1:3:7))
            @eval @test (1:3:7) ∩ ($T(3):$T(5)) == 4:3:4
            @eval @test ($T(3):$T(5)) ∩ (1:3:7) == 4:3:4
        end
        # Ranges of UInts are currently broken, cf. e.g.
        # * https://github.com/JuliaLang/julia/issues/29576
        # * https://github.com/JuliaLang/julia/issues/29801
        # * https://github.com/JuliaLang/julia/issues/29810
        # These (or other) bugs related to ranges of Unsigned integers prevent some of the
        # following tests from passing
        for T in halfuinttypes
            @eval @test_broken ($T(5/2):$T(3/2):$T(7)) ∩ ($T(1):$T(8)) == $T(4):$T(3):$T(7)
            @eval @test_broken ($T(1):$T(8)) ∩ ($T(5/2):$T(3/2):$T(7)) == $T(4):$T(3):$T(7)
        end
        for T in (:HalfUInt8, :HalfUInt16)
            @eval @test isempty((1:3:7) ∩ ($T(3/2):$T(5)))
            @eval @test isempty(($T(3/2):$T(5)) ∩ (1:3:7))
            @eval @test (1:3:7) ∩ ($T(3):$T(5)) == 4:3:4
            @eval @test ($T(3):$T(5)) ∩ (1:3:7) == 4:3:4
        end
        @static if Sys.WORD_SIZE == 64
            @test isempty((1:3:7) ∩ (HalfUInt32(3/2):HalfUInt32(5)))
            @test isempty((HalfUInt32(3/2):HalfUInt32(5)) ∩ (1:3:7))
            @test (1:3:7) ∩ (HalfUInt32(3):HalfUInt32(5)) == 4:3:4
            @test (HalfUInt32(3):HalfUInt32(5)) ∩ (1:3:7) == 4:3:4
        else
            @test_broken isempty((1:3:7) ∩ (HalfUInt32(3/2):HalfUInt32(5)))
            @test_broken isempty((HalfUInt32(3/2):HalfUInt32(5)) ∩ (1:3:7))
            @test_broken (1:3:7) ∩ (HalfUInt32(3):HalfUInt32(5)) == 4:3:4
            @test_broken (HalfUInt32(3):HalfUInt32(5)) ∩ (1:3:7) == 4:3:4
        end
        for T in (:HalfUInt64, :HalfUInt128)
            @eval @test_broken isempty((1:3:7) ∩ ($T(3/2):$T(5)))
            @eval @test_broken isempty(($T(3/2):$T(5)) ∩ (1:3:7))
            @eval @test_broken (1:3:7) ∩ ($T(3):$T(5)) == 4:3:4
            @eval @test_broken ($T(3):$T(5)) ∩ (1:3:7) == 4:3:4
        end
    end
end


text:
"""
    c

Speed of light = `299,792,458m/s`
"""
const c = 299792458  # m/s

"""
    k

Boltzman constant = `1.38×10⁻²³ J/K`
"""
const k = 1.380649e-23  # J/K

"""
    Rₑ

Rₑ (Volumetric Mean Radius of the Earth) = `6371.000e3m`

From [NASA Earth Fact Sheet](https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html)
"""
const Rₑ = 6371.000e3  # meters

"""
    L1_freq

L1 carrier frequency = `1575.42e6Hz`
"""
const L1_freq = 1575.42e6  # Hz

"""
    L2_freq

L2 carrier frequency = `1227.60e6Hz`
"""
const L2_freq = 1227.60e6  # Hz

"""
    L5_freq

L5 Frequency `1176.45e6Hz`
"""
const L5_freq = 1176.45e6  # Hz

"""
    h_parms_tcxo

TCXXO typical oscillator h-parameters. Taken from Table 1 in
*Digital GNSS PLL DesignConditioned on Thermal andOscillator Phase Noise* by
Curran & Murphy located [https://ieeexplore.ieee.org/document/6129629](https://ieeexplore.ieee.org/document/6129629)

Order: [h₋₂, h₋₁, h₀, h₁, h₂]
"""
const h_parms_tcxo = [[2.0e-20, 7.0e-21, 2.0e-19,       0,       0],
                      [2.0E-20, 1.0E-20, 1.0e-21,       0,       0],
                      [3.8e-21, 1.8e-19, 9.4e-20,       0,       0],
                      [2.4e-22, 2.4e-21, 3.9e-22, 2.5e-23, 2.0e-26],
                      [8.5e-22, 4.5e-20, 3.5e-20, 6.0e-23,       0],
                      [2.5e-23, 3.5e-21, 1.9e-21, 6.0e-23, 2.8e-24]]

"""
    h_parms_ocxo

OCXO typical oscillator h-parameters. Taken from Table 1 in
*Digital GNSS PLL DesignConditioned on Thermal andOscillator Phase Noise* by
Curran & Murphy located [https://ieeexplore.ieee.org/document/6129629](https://ieeexplore.ieee.org/document/6129629)

Order: [h₋₂, h₋₁, h₀, h₁, h₂]
"""
const h_parms_ocxo = [[4.0e-26, 3.6e-25, 2.6e-22, 2.0e-23, 2.2e-26],
                      [4.0e-23, 2.0e-21, 8.0e-20,       0,       0]]

"""
    h_parms_cesium

Cesium typical oscillator h-parameters. Taken from Table 1 in
*Digital GNSS PLL DesignConditioned on Thermal andOscillator Phase Noise* by
Curran & Murphy located [https://ieeexplore.ieee.org/document/6129629](https://ieeexplore.ieee.org/document/6129629)

Order: [h₋₂, h₋₁, h₀, h₁, h₂]
"""
const h_parms_cesium = [4.0e-29, 7.0e-23, 2.0e-20,      0,       0]

"""
    g1b1()

Data whose filenames contain the `g1b1` string are collected at 25 Msps,
with the center frequency set to 1.57GHz which may contain the following
signals:

- GPS L1
- Galileo E1bc
- Beidou B1i

Arguments: `None`

Returns:

- `(f_s, f_if, f_center, L1_freq, sigtype)::Tuple`: where
    * `f_s::Float64`: sampling rate of data `(25MHz)`
    * `f_if::Float64`: absolute value of the difference between the signal
                       carrier frequency and receiver center frequency `(5.42MHz)`
    * `f_center::Float64`: receiver center frequency `(1.57GHz)`
    * `L1_freq::Float64`: GPS L1 carrier frequency `(1.57542GHz)`
    * `sigtype::Val{:l1ca}`: `[DEPRICATED]` used in previous implementation of
                             GNSSTools and will be removed soon
"""
function g1b1()
    f_s = 25e6  # Hz
    f_center = 1.57e9  # Hz
    f_if = abs(L1_freq-f_center)  # Hz
    sigtype = Val(:l1ca)
    return (f_s, f_if, f_center, L1_freq, sigtype)
end

"""
    g5()

Data whose filenames contain the `g5` string are collected at 25 Msps,
with the center frequency set to 1.17645GHz which may contain the following
signals:

- GPS L5
- Galileo E5a

Arguments: `None`

Returns:

- `(f_s, f_if, f_center, L5_freq, sigtype)::Tuple`: where
    * `f_s::Float64`: sampling rate of data `(25MHz)`
    * `f_if::Float64`: absolute value of the difference between the signal
                       carrier frequency and receiver center frequency `(0Hz)`
    * `f_center::Float64`: receiver center frequency `(1.17645GHz)`
    * `L1_freq::Float64`: GPS L1 carrier frequency `(1.17645GHz)`
    * `sigtype::Val{:l5q}`: `[DEPRICATED]` used in previous implementation of
                            GNSSTools and will be removed soon
"""
function g5()
    f_s = 25e6  # Hz
    f_center = 1.17645e9  # Hz
    f_if = abs(L5_freq-f_center)  # Hz
    sigtype = Val(:l5q)
    return (f_s, f_if, f_center, L5_freq, sigtype)
end

"""
    WGS84

Struct type that hols the WGS84 constants.

Fields:

- `a::Float64`: semi-major axis of Earth `(6378137.0m)`
- `f::Float64`: flattening factor `(1/298.257223563)`
- `ω::Float64`: Earth's angular velocity `(7.2921151467e-5rad/s)`
- `μ::Float64`: also known as `GM`, is Earth's gravitational constant
                `(3.986005e14m³/s²)`
"""
struct WGS84
    a::Float64
    f::Float64
    ω::Float64
    μ::Float64
end

"""
    initWGS84()

Initializes the WGS84 constants and stores them in a `WGS84` struct.

Arguments: `None`

Returns:

- `WGS84` struct
"""
function initWGS84()
    return WGS84(6378137.0,  # meters
                 1/298.257223563,
                 7.2921151467e-5,  # rad/s
                 3.986005e14)  # m³/s²
end

"""
    wgs84

Initialized constant struct that contains the WGS84 constants.
"""
const wgs84 = initWGS84()

"""
    parseGPSData(data_file)

Parse the GPS data file from AGI. See function `getCurrentGPSNORADIDS` below.

Required Arguments:

- `data_file::String`: file name of data file

Returns:

- `GPSData::Dict`: multi-layer `Dict` containing the history of SVs for each PRN
                   see `getCurrentGPSNORADIDS` for more detailed information
"""
function parseGPSData(data_file)
    file = open(data_file, "r")
    lines = readlines(file)
    close(file)
    N = size(lines)[1]
    GPSData = Dict{Int,Dict{Int,Dict{String,Any}}}()
    for i in 23:N
        line = lines[i]
        if length(line) > 1
            prn = parse(Int, line[3:4])
            sv = parse(Int, line[7:8])
            satnum = parse(Int, line[11:15])
            block = line[18:20]
            start_week = parse(Int, line[27:30])
            start_toa = parse(Int, line[32:37])
            start_year = parse(Int, line[39:42])
            start_month = parse(Int, line[44:45])
            start_day = parse(Int, line[47:48])
            start_date = (start_year, start_month, start_day)
            start_julian_date = DatetoJD(start_date..., 0, 0, 0)
            active = occursin("Active", line)
            notes = line[77:end]
            if active
                end_week = missing
                end_toa = missing
                end_date = missing
                end_julian_date = missing
            else
                end_week = parse(Int, line[53:56])
                end_toa = parse(Int, line[58:63])
                end_year = parse(Int, line[65:68])
                end_month = parse(Int, line[70:71])
                end_day = parse(Int, line[73:74])
                end_date = (end_year, end_month, end_day)
                end_julian_date = DatetoJD(end_date..., 0, 0, 0)
            end
            sv_info = Dict("satnum"=>satnum,
                           "block"=>block,
                           "start_week"=>start_week,
                           "start_toa"=>start_toa,
                           "start_date"=>start_date,
                           "start_julian_date"=>start_julian_date,
                           "end_week"=>end_week,
                           "end_toa"=>end_toa,
                           "end_date"=>end_date,
                           "end_julian_date"=>end_julian_date,
                           "active"=>active,
                           "status"=>notes)
            try
                GPSData[prn][sv] = sv_info
            catch
                GPSData[prn] = Dict(sv=>sv_info)
            end
        end
    end
    return GPSData
end

"""
    getCurrentGPSNORADIDS(file_url)

Download GPS SV data from AGI and parse. Returns a dictionary with historical
SVs for each PRN number.

Required Arguments:

- `file_url::String`: the file URL `(Default = GPS_data_file_url)
    * See `const GPS_data_file_url` below

Returns:

- `GPSData::Dict`: multi-layer `Dict` containing the history of SVs for each PRN
    * top level is sorted by PRN number
    * `GPSData[prn]` yields another set of dictionaries, sorted by SV
    * `GPSData[prn][sv]` yields a `Dict` containing attributes for that specific
      SV such as
        + `"active"`: `Bool` value if set to `true`, indicates that SV is still
                      in service and the fields `"end_week"`, `"end_toa"`,
                      `"end_date"`, and `"end_julian_date"` are all set to
                      `misssing`
        + `"status"`: `String` that either says `"active"`, if `"active"` is
                      `true` or `"Removed from almanac"` if it is `false`
        + `"satnum"`: the SV catalog number for Space-Track
        + `"block"`: the satellite version number
        + `"start_week"`: GPS week when SV goes into service
        + `"start_toa"`: GPS seconds in the week when SV goes into service
        + `"start_date"`: `Tuple` in the format (year, month, day) representing
                          the date when the SV goes into service
        + `"start_julian_date"`: `"start_date"` above converted to the Julian
                                 Date when the SV goes into service
        + `"end_week"`: GPS week when SV goes out of service
        + `"end_toa"`: GPS seconds in the week when SV goes out of service
        + `"end_date"`: `Tuple` in the format (year, month, day) representing
                        the date when the SV goes out of service
        + `"end_julian_date"`: `"start_date"` above converted to the Julian
                               Date when the SV goes out of service
"""
function getCurrentGPSNORADIDS(file_url)
    # Search for existing `.GNSSTools` folder in user home directory and create
    # one if it does not already exist
    directory = string(homedir(), "/.GNSSTools")
    if ~isdir(directory)
        mkpath(directory)
    end
    # Check for `GPSData.txt` file to parse and if it does not exist, download
    # it from the url, `file_url`
    ids_file = string(directory, "/GPSData.txt")
    if ~isfile(ids_file)
        run(`curl -o $(ids_file) $(file_url)`)
    end
    # Parse GPSData file and return the `Dict` containing the history of
    # Sv usage of each PRN number
    GPSData = parseGPSData(ids_file)
    return GPSData
end

"""
    GPS_data_file_url

URL location of the historical and current GPS PRN data file which contains info
on each SV that used each PRN number throughout the U.S. GPS constellation's
history.

url = `ftp://ftp.agi.com/pub/Catalog/Almanacs/SEM/GPSData.txt`
"""
const GPS_data_file_url = "ftp://ftp.agi.com/pub/Catalog/Almanacs/SEM/GPSData.txt"

"""
    GPSData

Dictionary containing information about each SV that used each PRN numbers in
the U.S. GPS constellation throughout its history. The format of the dictionary
is `GPSData[prn][sv]["[parameter]"]`. See `getCurrentGPSNORADIDS` for more
detail.
"""
const GPSData = getCurrentGPSNORADIDS(GPS_data_file_url)

"""
    updateCurrentGPSNORADIDS()

Redownload the historical GPS PRN data file. It will be reparsed when
`GNSSTools` is imported in a new Julia session.
"""
function updateCurrentGPSNORADIDS()
    directory = string(homedir(), "/.GNSSTools")
    ids_file = string(directory, "/GPSData.txt")
    run(`curl -o $(ids_file) $(file_url)`)
end


text:
"""
        function inferparameters!(model,samples,trans,stepsize::Float64,D::Diagonal,q₀;serowaningrate = 1/180,num_chains = 1)

Infer the unknown transmission and observation parameters for the `k`th county in by `country_model`. Inference is performed by drawing `samples` number of replicates from the posterior distribution with log-density
    given in the model with sero-waning rate fixed. The parameters and their transformation are encoded in `trans`. Initial parameter search starts at `q₀` in the transformed domain defined by `trans`.
    whilst the initial kinetic energy search begins at `D` and the HMC step size is fixed to be `stepsize`. Inference output is stored in-place.

    Serowaning rate can be added, as well as additional chains for cross-comparison and validation.

"""
function inferparameters!(model,samples,trans,stepsize::Float64,D::Diagonal,q₀;serowaningrate = 1/365,num_chains = 1,σ = 0.16,ω = 1/180)
        println("Starting MCMC parameter inference for $(model.areaname)")
        n = length(q₀)
        l_area(x) = model.log_likelihood(x,model,serowaningrate;σ = σ,ω = ω) + model.log_priors(x)
        ℓ = TransformedLogDensity(trans, l_area)#transformed log-likelihood
        ∇ℓ = LogDensityProblems.ADgradient(:ForwardDiff, ℓ)#transformed log-likelihood gradient wrt the parameters

        results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇ℓ, samples,
                                initialization = (q = q₀,κ=GaussianKineticEnergy(D),ϵ = stepsize),
                                warmup_stages = fixed_stepsize_warmup_stages(M=Symmetric),
                                reporter = NoProgressReport())
        # results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇ℓ, samples,
        #                         initialization = (q = q₀,κ=GaussianKineticEnergy(D)))

        transformed_results = TransformVariables.transform.(trans,results.chain)
        val = zeros(length(transformed_results),length(transformed_results[1]),1)

        for i = 1:size(val,1),j = 1:size(val,2)
                val[i,j,1] = transformed_results[i][j]
        end

        for chain_reps = 2:num_chains
                println("Running additional chain $(chain_reps) for $(model.areaname)")
                results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇ℓ, samples,
                                        initialization = (q = q₀,κ=GaussianKineticEnergy(D),ϵ = stepsize),
                                        warmup_stages = fixed_stepsize_warmup_stages(M=Symmetric),
                                        reporter = NoProgressReport())
                transformed_results = TransformVariables.transform.(trans,results.chain)
                for i = 1:size(val,1),j = 1:size(val,2)
                        val[i,j,chain_reps] = transformed_results[i][j]
                end
        end

        chn = Chains(val,[String(k) for k in keys(transformed_results[1])])

        model.MCMC_results = MCMCResults(chn,
                                        [l_area(transformed_results[i]) - model.log_priors(transformed_results[i]) for i = 1:length(transformed_results)],
                                        results.tree_statistics)

        return nothing
end

"""
        function inferparameters!(model,samples,trans,stepsize::Float64,D::Diagonal,searchrange;serowaningrate = 1/180,num_chains = 1)

Infer the unknown transmission and observation parameters for the `k`th county in by `country_model`. Inference is performed by drawing `samples` number of replicates from the posterior distribution with log-density
    given in the model with sero-waning rate fixed. The parameters and their transformation are encoded in `trans`. Initial parameter search has a pre-optimisation using the default adaptive differential evolution algorithm implemented by the `BlackBoxOptim` package,
    whilst the initial kinetic energy search begins at `D` and the HMC step size is fixed to be `stepsize`. Inference output is stored in-place.

    serowaning rate can be added, as well as additional chains for cross-comparison and validation.
"""
function inferparameters!(model,samples,trans,stepsize::Float64,D::Diagonal;serowaningrate = 1/365,num_chains = 1,σ = 0.16,ω = 1/180)
        println("Searching for a good initial condition for $(model.areaname)")
        l_area(x) = model.log_likelihood(x,model,serowaningrate;σ = σ,ω = ω) + model.log_priors(x)
        f(x) = -transform_logdensity(trans, l_area,x)
        searchrange = fill((-3.,3.),TransformVariables.dimension(trans))
        res = bboptimize(f; SearchRange = searchrange,PopulationSize=1000,MaxSteps=30000,TraceMode = :silent)
        q₀ = best_candidate(res)
        inferparameters!(model,samples,trans,stepsize,D,q₀;serowaningrate=serowaningrate,num_chains=num_chains,σ = σ,ω = ω)
        return nothing
end


text:
"""

    addtodiag!(A, v)

Efficiently add the vector `v` in-place to the main diagonal of the matrix `A` (overwriting
the latter). `length(v)` must match one of the dimensions of `A`, with the other dimension
being at least as large, and the element types of `v` must be promotable to those of `A`.
If so, yield the same result as `A .= diagm(v) .+ A`.

# Examples
```jldoctest
julia> A = [1.0 2.0 3.0; 2.0 4.0 6.0; 3.0 6.0 9.0];

julia> v = [1.0; 2.0; 3.0];

julia> addtodiag!(A, v)
3×3 Array{Float64,2}:
 2.0  2.0   3.0
 2.0  6.0   6.0
 3.0  6.0  12.0
```

See also: [`addtodiag`](@ref), [`subtractfromdiag!`](@ref), [`subtractfromdiag`](@ref)

"""
function addtodiag!(A::AbstractArray{<:Number, 2}, v::AbstractArray{<:Number, 1})

    if (length(v) == size(A, 1) && length(v) <= size(A, 2))
        @inbounds for jj=1:size(A, 1)
            A[jj, jj] += v[jj];
        end
    elseif (length(v) == size(A, 2) && length(v) <= size(A, 1))
        @inbounds for jj=1:size(A, 2)
            A[jj, jj] += v[jj];
        end
    else
        error(string("addtodiag!(): Length of v (currently ", length(v),
                     ") must match at least one dimension of A (currently", size(A), ")"));
    end

end

"""

    addtodiag(A, v)

Efficiently add the vector `v` not in-place to the main diagonal of the matrix `A`,
returning the result. `length(v)` must match one of the dimensions of `A`, with the other
dimension being at least as large, and the element types of `v` must be promotable to those
of `A`. If so, yield the same result as `diagm(v) .+ A`.

# Examples
```jldoctest
julia> A = [1.0 2.0 3.0; 2.0 4.0 6.0; 3.0 6.0 9.0];

julia> v = [1.0; 2.0; 3.0];

julia> addtodiag(A, v)
3×3 Array{Float64,2}:
 2.0  2.0   3.0
 2.0  6.0   6.0
 3.0  6.0  12.0
```

See also: [`addtodiag!`](@ref), [`subtractfromdiag!`](@ref), [`subtractfromdiag`](@ref)

"""
function addtodiag(A::AbstractArray{<:Number, 2}, v::AbstractArray{<:Number, 1})
    B = deepcopy(A);
    addtodiag!(B, v);
    return B;
end

"""

    subtractfromdiag!(A, v)

Efficiently subtract the main diagonal of the matrix `A` from the vector `v` in-place
(overwriting the latter). `length(v)` must match one of the dimensions of `A`, with the
other dimension being at least as large, and the element types of `v` must be promotable to
those of `A`. If so, yield the same result as `A .= diagm(v) .- A`.

# Examples
```jldoctest
julia> A = [1.0 2.0 3.0; 2.0 4.0 6.0; 3.0 6.0 9.0];

julia> v = [1.0; 2.0; 3.0];

julia> subtractfromdiag!(A, v)
3×3 Array{Float64,2}:
  0.0  -2.0  -3.0
 -2.0  -2.0  -6.0
 -3.0  -6.0  -6.0
```

See also: [`addtodiag!`](@ref), [`addtodiag`](@ref), [`subtractfromdiag`](@ref)

"""
function subtractfromdiag!(A::AbstractArray{<:Number, 2}, v::AbstractArray{<:Number, 1})
    lmul!(-1.0, A);
    addtodiag!(A, v);
end

"""

    subtractfromdiag(A, v)

Efficiently subtract the main diagonal of the matrix `A` from the vector `v` not in-place,
returning the result. `length(v)` must match one of the dimensions of `A`, with the other
dimension being at least as large, and the element types of `v` must be promotable to those
of `A`. If so, yield the same result as `diagm(v) .- A`.

# Examples
```jldoctest
julia> A = [1.0 2.0 3.0; 2.0 4.0 6.0; 3.0 6.0 9.0];

julia> v = [1.0; 2.0; 3.0];

julia> subtractfromdiag(A, v)
3×3 Array{Float64,2}:
 1.0  2.0  3.0
 2.0  4.0  6.0
 3.0  6.0  9.0
```

See also: [`addtodiag!`](@ref), [`addtodiag`](@ref), [`subtractfromdiag!`](@ref)

"""
function subtractfromdiag(A::AbstractArray{<:Number, 2}, v::AbstractArray{<:Number, 1})
    B = deepcopy(A);
    subtractfromdiag!(B, v);
    return B;
end

"""

    mulMatDiagleft!(A, v)

Efficiently multiply each row of the matrix `A` by the corresponding element of the vector
`v` in-place (overwriting the former). The element types of `v` must be promotable to those
of `A`. If so, yield the same result as `A .= diagm(v) * A`.

# Examples
```jldoctest
julia> A = [1.0 2.0 3.0; 2.0 4.0 6.0; 3.0 6.0 9.0];

julia> v = [1.0; 2.0; 3.0];

julia> mulMatDiagleft!(A, v)
3×3 Array{Float64,2}:
 1.0   2.0   3.0
 4.0   8.0  12.0
 9.0  18.0  27.0
```

See also: [`mulMatDiagleft`](@ref), [`mulMatDiagright!`](@ref), [`mulMatDiagright`](@ref)

"""
mulMatDiagleft!(A::AbstractArray{<:Number, 2}, v::AbstractArray{<:Number, 1}) = broadcast!(*, A, A, v);

"""

    mulMatDiagleft(A, v)

Efficiently multiply each row of the matrix `A` by the corresponding element of the vector
`v` not in-place, returning the result. The element types of `v` must be promotable to those
of `A`. If so, yield the same result as `diagm(v) * A`.

# Examples
```jldoctest
julia> A = [1.0 2.0 3.0; 2.0 4.0 6.0; 3.0 6.0 9.0];

julia> v = [1.0; 2.0; 3.0];

julia> mulMatDiagleft(A, v)
3×3 Array{Float64,2}:
 1.0   2.0   3.0
 4.0   8.0  12.0
 9.0  18.0  27.0
```

See also: [`mulMatDiagleft!`](@ref), [`mulMatDiagright!`](@ref), [`mulMatDiagright`](@ref)

"""
mulMatDiagleft(A::AbstractArray{<:Number, 2}, v::AbstractArray{<:Number, 1}) = broadcast(*, A, v);

"""

    mulMatDiagright!(A, v)

Efficiently multiply each column of the matrix `A` by the corresponding element of the
vector `v` in-place (overwriting the former). The element types of `v` must be promotable to
those of `A`, and `v` must have the same length as the number of columns of `A`. If so,
yield the same result as `A .= A * diagm(v)`.

# Examples
```jldoctest

julia> A = [1.0 2.0 3.0; 2.0 4.0 6.0; 3.0 6.0 9.0];

julia> v = [1.0; 2.0; 3.0];

julia> mulMatDiagright!(A, v)
3×3 Array{Float64,2}:
 1.0   4.0   9.0
 2.0   8.0  18.0
 3.0  12.0  27.0
```

See also: [`mulMatDiagleft!`](@ref), [`mulMatDiagleft`](@ref), [`mulMatDiagright`](@ref)

"""
function mulMatDiagright!(A::AbstractArray{<:Number, 2}, v::AbstractArray{<:Number, 1})

    if (length(v) != size(A, 2))
        error("mulMatDiagright!(): v (currently of length ", length(v),
              ") must have the same length as the number of columns in A (currently ",
              size(A, 2), ")");
    end

    @inbounds for jj=1:size(A, 2)
        @inbounds for ii=1:size(A, 1)
            A[ii, jj] *= v[jj];
        end
    end

end

"""

    mulMatDiagright(A, v)

Efficiently multiply each column of the matrix `A` by the corresponding element of the
vector `v` not in-place, returning the result. The element types of `v` must be promotable
to those of `A`. If so, yield the same result as `A * diagm(v)`.

# Examples
```jldoctest
julia> A = [1.0 2.0 3.0; 2.0 4.0 6.0; 3.0 6.0 9.0];

julia> v = [1.0; 2.0; 3.0];

julia> mulMatDiagright(A, v)
3×3 Array{Float64,2}:
 1.0   4.0   9.0
 2.0   8.0  18.0
 3.0  12.0  27.0
```

See also: [`mulMatDiagleft!`](@ref), [`mulMatDiagleft`](@ref), [`mulMatDiagright!`](@ref)

"""
function mulMatDiagright(A::AbstractArray{<:Number, 2}, v::AbstractArray{<:Number, 1})
    B = deepcopy(A);
    mulMatDiagright!(B, v);
    return B;
end

"""

    vecNto3N(v)

Return a vector `w` that is 3 times the length of `v`, such that each 3-element block
`w[3*ii-2:3*ii]` has all elements identical to `v[ii]`. Mathematically, the vector
```math
v = \\begin{bmatrix}
v_1 &
v_2 &
\\ldots &
v_N
\\end{bmatrix}^{\\top}
```
is mapped to the vector
```math
w = \\begin{bmatrix}
v_1 &
v_1 &
v_1 &
v_2 &
v_2 &
v_2 &
\\ldots &
v_N &
v_N &
v_N
\\end{bmatrix}^{\\top}
```
# Examples
```jldoctest
julia> v = [1.0; 2.0; 3.0];

julia> vecNto3N(v)
9-element Array{Float64,1}:
 1.0
 1.0
 1.0
 2.0
 2.0
 2.0
 3.0
 3.0
 3.0
```

"""
vecNto3N(v::AbstractArray{<:Number, 1}) = vec(repeat(transpose(v), 3, 1));

"""

    sympart!(A)

Convert a square matrix `A` into its symmetric part in-place (overwriting `A`), yielding the
same result as `A .= (A .+ transpose(A))./2`.

# Examples
```jldoctest
julia> A = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0];

julia> B = 1im .* A;

julia> sympart!(A)
3×3 Array{Float64,2}:
 1.0  3.0  5.0
 3.0  5.0  7.0
 5.0  7.0  9.0

julia> sympart!(B)
3×3 Array{Complex{Float64},2}:
 0.0+1.0im  0.0+3.0im  0.0+5.0im
 0.0+3.0im  0.0+5.0im  0.0+7.0im
 0.0+5.0im  0.0+7.0im  0.0+9.0im
```

See also: [`sympart`](@ref), [`hermpart!`](@ref), [`hermpart`](@ref), [`ahermpart!`](@ref),
[`ahermpart`](@ref)

"""
function sympart!(A::AbstractArray{<:Number, 2})

    if (size(A, 1) != size(A, 2))
        error("sympart!(): A must be a square matrix");
    end

    @inbounds for jj=1:size(A, 2)
        @inbounds for ii=(jj + 1):size(A, 1)
            A[ii, jj] = (A[ii, jj] + A[jj, ii])/2;
            A[jj, ii] = A[ii, jj];
        end
    end

end

"""

    sympart(A)

Convert a square matrix `A` into its symmetric part not in-place, returning the result, and
yielding the same result as `(A .+ transpose(A))./2`.

# Examples
```jldoctest
julia> A = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0];

julia> B = 1im .* A;

julia> sympart(A)
3×3 Array{Float64,2}:
 1.0  3.0  5.0
 3.0  5.0  7.0
 5.0  7.0  9.0

julia> sympart(B)
3×3 Array{Complex{Float64},2}:
 0.0+1.0im  0.0+3.0im  0.0+5.0im
 0.0+3.0im  0.0+5.0im  0.0+7.0im
 0.0+5.0im  0.0+7.0im  0.0+9.0im
```

See also: [`sympart!`](@ref), [`hermpart!`](@ref), [`hermpart`](@ref), [`ahermpart!`](@ref),
[`ahermpart`](@ref)

"""
function sympart(A::AbstractArray{<:Number, 2})
    B = deepcopy(A);
    sympart!(B);
    return B;
end

"""

    hermpart!(A)

Convert a square matrix `A` into its Hermitian part in-place (overwriting `A`), yielding the
same result as `A .= (A .+ adjoint(A))./2`. (This automatically uses `sympart!()` if `A` is
real to avoid conflict with promotion to complex types.)

# Examples
```jldoctest
julia> A = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0];

julia> B = 1im .* A;

julia> hermpart!(A)
3×3 Array{Float64,2}:
 1.0  3.0  5.0
 3.0  5.0  7.0
 5.0  7.0  9.0

julia> hermpart!(B)
3×3 Array{Complex{Float64},2}:
 0.0+0.0im  0.0-1.0im  0.0-2.0im
 0.0+1.0im  0.0+0.0im  0.0-1.0im
 0.0+2.0im  0.0+1.0im  0.0+0.0im
```

See also: [`sympart!`](@ref), [`sympart`](@ref), [`hermpart`](@ref), [`ahermpart!`](@ref),
[`ahermpart`](@ref)

"""
function hermpart!(A::AbstractArray{<:Complex{<:Real}, 2})

    if (size(A, 1) != size(A, 2))
        error("hermpart!(): A must be a square matrix");
    end

    @inbounds for jj=1:size(A, 2)
        A[jj, jj] = real(A[jj, jj]);
        @inbounds for ii=(jj + 1):size(A, 1)
            A[ii, jj] = (A[ii, jj] + conj(A[jj, ii]))/2;
            A[jj, ii] = conj(A[ii, jj]);
        end
    end

end
hermpart!(A::AbstractArray{<:Real, 2}) = sympart!(A);

"""

    hermpart(A)

Convert a square matrix `A` into its Hermitian part not in-place, returning the result, and
yielding the same result as `(A .+ adjoint(A))./2`. (Automatically use `sympart()` if `A` is
real.)

# Examples
```jldoctest
julia> A = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0];

julia> B = 1im .* A;

julia> hermpart(A)
3×3 Array{Float64,2}:
 1.0  3.0  5.0
 3.0  5.0  7.0
 5.0  7.0  9.0

julia> hermpart(B)
3×3 Array{Complex{Float64},2}:
 0.0+0.0im  0.0-1.0im  0.0-2.0im
 0.0+1.0im  0.0+0.0im  0.0-1.0im
 0.0+2.0im  0.0+1.0im  0.0+0.0im
```

See also: [`sympart!`](@ref), [`sympart`](@ref), [`hermpart!`](@ref), [`ahermpart!`](@ref),
[`ahermpart`](@ref)

"""
function hermpart(A::AbstractArray{<:Complex{<:Real}, 2})
    B = deepcopy(A);
    hermpart!(B);
    return B;
end
hermpart(A::AbstractArray{<:Real, 2}) = sympart(A);

"""

    ahermpart!(A)

Convert a square matrix `A` into its antiHermitian part in-place (overwriting `A`), yielding
the same result as `A .= (A .- adjoint(A))./2im`. (Note that the result is Hermitian.)

# Examples
```jldoctest
julia> B = 1im .* [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0];

julia> ahermpart!(B)
3×3 Array{Complex{Float64},2}:
 1.0+0.0im  3.0-0.0im  5.0-0.0im
 3.0+0.0im  5.0+0.0im  7.0-0.0im
 5.0+0.0im  7.0+0.0im  9.0+0.0im
```

See also: [`sympart!`](@ref), [`sympart`](@ref), [`hermpart!`](@ref), [`hermpart`](@ref),
[`ahermpart`](@ref)

"""
function ahermpart!(A::AbstractArray{<:Complex{<:Real}, 2})
    lmul!(-1im, A);
    hermpart!(A);
end

"""

    ahermpart(A)

Convert a square matrix `A` into its antiHermitian part not in-place, returning the result,
and yielding the same result as `(A .- adjoint(A))./2im`. (Note that the result is
Hermitian.)

# Examples
```jldoctest
julia> B = 1im .* [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0];

julia> ahermpart(B)
3×3 Array{Complex{Float64},2}:
 1.0+0.0im  3.0-0.0im  5.0-0.0im
 3.0+0.0im  5.0+0.0im  7.0-0.0im
 5.0+0.0im  7.0+0.0im  9.0+0.0im
```

See also: [`sympart!`](@ref), [`sympart`](@ref), [`hermpart!`](@ref), [`hermpart`](@ref),
[`ahermpart!`](@ref)

"""
function ahermpart(A::AbstractArray{<:Complex{<:Real}, 2})
    B = deepcopy(A);
    ahermpart!(B);
    return B;
end


text:
using Sockets
server = listen(2001)

@async begin
    while true
        sock = accept(server)
        @async while isopen(sock)
            line = readline(sock, keep=true)
            write(sock, line)
            write(stdout, line)
        end
    end
end

conn = connect(2001)
call(buff::IOBuffer) = @async begin
    write(buff, "\r\n")
    write(conn, take!(buff))
    write(stdout, readline(conn, keep=true))
    println(".........")
end

max_stars_repo_name: rsd-consultorias/scripts-engenharia

max_issues_repo_name: rsd-consultorias/scripts-engenharia

max_forks_repo_name: rsd-consultorias/scripts-engenharia


text:
# Invertible Hyperbolic network from Keegan et al (2018)
# Author: Gabrio Rizzuti, rizzuti.gabrio@gatech.edu
# Date: August 2020

export NetworkHyperbolic1D, clear_grad!, get_params

"""
    H = NetworkHyperbolic1D(nx, ny, n_in, batchsize, n_hidden, depth; k1=3, k2=3, p1=1, p2=1, s1=1, s2=1)
 Create an hyperbolic network for data-driven generative modeling based
 on the change of variables formula.
 *Input*:
 - `nx`, `ny`, `n_in`, `batchsize`: spatial dimensions, number of channels and batchsize of input tensors `X`
 - `n_hidden`: number of hidden units in residual blocks
 - `depth`: number network layers
 - `k1`, `k2`: kernel size for first and third residual layer (`k1`) and second layer (`k2`)
 - `p1`, `p2`: respective padding sizes for residual block layers
 - `s1`, `s2`: respective strides for residual block layers
 *Output*:
 - `H`: HINT network
 *Usage:*
 - Forward mode: `Z, logdet = H.forward(X)`
 - Inverse mode: `X = H.inverse(Z)`
 - Backward mode: `ΔX, X = H.backward(ΔZ, Z)`
 *Trainable parameters:*
 - None in `H` itself
 - Trainable parameters in activation normalizations `H.AN[i]`, and in coupling layers `H.CL[i]`, where `i` ranges from `1` to `depth`.
 See also: [`ActNorm`](@ref), [`CouplingLayerHINT`](@ref), [`get_params`](@ref), [`clear_grad!`](@ref)
"""
mutable struct NetworkHyperbolic1D <: InvertibleNetwork
    AN::AbstractArray{ActNorm, 1}
    CL::AbstractArray{HyperbolicLayer, 2}
    logdet::Bool
    is_reversed::Bool
end

@Flux.functor NetworkHyperbolic1D

# Constructor
function NetworkHyperbolic1D(n_in::Int64, batchsize::Int64, n_hidden::Int64, depth::Int64; logdet=true, depthHL=1)

    AN = Array{ActNorm}(undef, depth)
    CL = Array{HyperbolicLayer, 2}(undef, depth, depthHL)

    # Create layers
    for j = 1:depth
        AN[j] = ActNorm(n_in; logdet=logdet)
        for k = 1:depthHL
            CL[j, k] = HyperbolicLayer(1, 1, Int(n_in/2), batchsize, 1, 1, 0; action="same", α=1f0, hidden_factor=Int(n_hidden/n_in))
        end
    end

    return NetworkHyperbolic1D(AN, CL, logdet, false)
end

# Forward pass and compute logdet
function forward(X, H::NetworkHyperbolic1D; logdet=nothing)
    isnothing(logdet) ? logdet = (H.logdet && ~H.is_reversed) : logdet = logdet

    depth, depthHL = size(H.CL)
    logdet_ = 0f0
    for j=1:depth
        logdet ? (X, logdet1) = H.AN[j].forward(X) : X = H.AN[j].forward(X)
        X_prev, X_curr = tensor_split(X)
        for k = 1:depthHL
            X_prev, X_curr = H.CL[j, k].forward(X_prev, X_curr)
        end
        X = tensor_cat(X_prev, X_curr)
        logdet && (logdet_ += logdet1)
    end
    logdet ? (return X, logdet_) : (return X)
end

# Inverse pass and compute gradients
function inverse(Z, H::NetworkHyperbolic1D; logdet=nothing)
    isnothing(logdet) ? logdet = (H.logdet && H.is_reversed) : logdet = logdet

    depth, depthHL = size(H.CL)
    logdet_ = 0f0
    for j = depth:-1:1
        Z_curr, Z_new = tensor_split(Z)
        for k = depthHL:-1:1
            Z_curr, Z_new = H.CL[j, k].inverse(Z_curr, Z_new)
        end
        Z = tensor_cat(Z_curr, Z_new)
        logdet ? (Z, logdet2) = H.AN[j].inverse(Z; logdet=true) : Z = H.AN[j].inverse(Z; logdet=false)
        logdet && (logdet_ += (logdet1 + logdet2))
    end
    logdet ? (return Z, logdet_) : (return Z)
end

# Backward pass and compute gradients
function backward(ΔZ, Z, H::NetworkHyperbolic1D)
    depth, depthHL = size(H.CL)
    for j = depth:-1:1
        ΔZ_curr, ΔZ_new = tensor_split(ΔZ)
        Z_curr, Z_new = tensor_split(Z)
        for k = depthHL:-1:1
            ΔZ_curr, ΔZ_new, Z_curr, Z_new = H.CL[j, k].backward(ΔZ_curr, ΔZ_new, Z_curr, Z_new)
        end
        ΔZ = tensor_cat(ΔZ_curr, ΔZ_new)
        Z = tensor_cat(Z_curr, Z_new)
        ΔZ, Z = H.AN[j].backward(ΔZ, Z)
    end
    return ΔZ, Z
end

# Clear gradients
function clear_grad!(H::NetworkHyperbolic1D)
    depth, depthHL = size(H.CL)
    for j=1:depth
        clear_grad!(H.AN[j])
        for k = 1:depthHL
            clear_grad!(H.CL[j, k])
        end
    end
end

# Get parameters
function get_params(H::NetworkHyperbolic1D)
    depth, depthHL = size(H.CL)
    p = []
    for j=1:depth
        p = cat(p, get_params(H.AN[j]); dims=1)
        for k = 1:depthHL
            p = cat(p, get_params(H.CL[j, k]); dims=1)
        end
    end
    return p
end

# Set is_reversed flag in full network tree
function tag_as_reversed!(H::NetworkHyperbolic1D, tag::Bool)
    depth, depthHL = size(H.CL)
    H.is_reversed = tag
    for j=1:depth
        tag_as_reversed!(H.AN[j], tag)
        for k = 1:depthHL
            tag_as_reversed!(H.CL[j, k], tag)
        end
    end
    return H
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions to compute the nutation according to IAU-76/FK5.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
#   [2] Vallado, D. A (06-Feb-2018). Consolidated Errata of Fundamentals of
#       Astrodynamics and Applications 4th Ed.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export nutation_fk5

################################################################################
#                                  Constants
################################################################################

################################################################################
#
# 1980 IAU Theory of Nutation Coefficients
# ========================================
#
# Those coefficients can be found on:
#
#   Seidelmann, P. K. 1980 IAU theory of nutation - The final report of the IAU
#   Working Group on Nutation. Celestial Mechanics, vol. 27, p. 79-106.
#
# However, the .dat file that was used to create this matrix was obtained from
# this website:
#
#   http://hpiers.obspm.fr/eop-pc/models/nutations/nut.html
#
# Notice that the order of them is not equal to that presented in the original
# paper, but this have no impact when computing the nutation. Moreover, this is
# the same order that is presented in [1](p. 1043).
#
################################################################################

const nut_coefs_1980 = [
  # Notation used in [1, p. 226].
  #
  # an1  an2  an3  an4  an5        Ai            Bi         Ci            Di
  #                  Units:  [0.0001"]  [0.0001"/JC]  [0.0001"]  [0.0001"/JC]
      0    0    0    0    1  -171996.0       -174.2    92025.0           8.9;
      0    0    2   -2    2   -13187.0         -1.6     5736.0          -3.1;
      0    0    2    0    2    -2274.0         -0.2      977.0          -0.5;
      0    0    0    0    2     2062.0          0.2     -895.0           0.5;
      0   -1    0    0    0    -1426.0          3.4       54.0          -0.1;
      1    0    0    0    0      712.0          0.1       -7.0           0.0;
      0    1    2   -2    2     -517.0          1.2      224.0          -0.6;
      0    0    2    0    1     -386.0         -0.4      200.0           0.0;
      1    0    2    0    2     -301.0          0.0      129.0          -0.1;
      0   -1    2   -2    2      217.0         -0.5      -95.0           0.3;
     -1    0    0    2    0      158.0          0.0       -1.0           0.0;
      0    0    2   -2    1      129.0          0.1      -70.0           0.0;
     -1    0    2    0    2      123.0          0.0      -53.0           0.0;
      1    0    0    0    1       63.0          0.1      -33.0           0.0;
      0    0    0    2    0       63.0          0.0       -2.0           0.0;
     -1    0    2    2    2      -59.0          0.0       26.0           0.0;
     -1    0    0    0    1      -58.0         -0.1       32.0           0.0;
      1    0    2    0    1      -51.0          0.0       27.0           0.0;
     -2    0    0    2    0      -48.0          0.0        1.0           0.0;
     -2    0    2    0    1       46.0          0.0      -24.0           0.0;
      0    0    2    2    2      -38.0          0.0       16.0           0.0;
      2    0    2    0    2      -31.0          0.0       13.0           0.0;
      2    0    0    0    0       29.0          0.0       -1.0           0.0;
      1    0    2   -2    2       29.0          0.0      -12.0           0.0;
      0    0    2    0    0       26.0          0.0       -1.0           0.0;
      0    0    2   -2    0      -22.0          0.0        0.0           0.0;
     -1    0    2    0    1       21.0          0.0      -10.0           0.0;
      0    2    0    0    0       17.0         -0.1        0.0           0.0;
      0    2    2   -2    2      -16.0          0.1        7.0           0.0;
     -1    0    0    2    1       16.0          0.0       -8.0           0.0;
      0    1    0    0    1      -15.0          0.0        9.0           0.0;
      1    0    0   -2    1      -13.0          0.0        7.0           0.0;
      0   -1    0    0    1      -12.0          0.0        6.0           0.0;
      2    0   -2    0    0       11.0          0.0        0.0           0.0;
     -1    0    2    2    1      -10.0          0.0        5.0           0.0;
      1    0    2    2    2       -8.0          0.0        3.0           0.0;
      0   -1    2    0    2       -7.0          0.0        3.0           0.0;
      0    0    2    2    1       -7.0          0.0        3.0           0.0;
      1    1    0   -2    0       -7.0          0.0        0.0           0.0;
      0    1    2    0    2        7.0          0.0       -3.0           0.0;
     -2    0    0    2    1       -6.0          0.0        3.0           0.0;
      0    0    0    2    1       -6.0          0.0        3.0           0.0;
      2    0    2   -2    2        6.0          0.0       -3.0           0.0;
      1    0    0    2    0        6.0          0.0        0.0           0.0;
      1    0    2   -2    1        6.0          0.0       -3.0           0.0;
      0    0    0   -2    1       -5.0          0.0        3.0           0.0;
      0   -1    2   -2    1       -5.0          0.0        3.0           0.0;
      2    0    2    0    1       -5.0          0.0        3.0           0.0;
      1   -1    0    0    0        5.0          0.0        0.0           0.0;
      1    0    0   -1    0       -4.0          0.0        0.0           0.0;
      0    0    0    1    0       -4.0          0.0        0.0           0.0;
      0    1    0   -2    0       -4.0          0.0        0.0           0.0;
      1    0   -2    0    0        4.0          0.0        0.0           0.0;
      2    0    0   -2    1        4.0          0.0       -2.0           0.0;
      0    1    2   -2    1        4.0          0.0       -2.0           0.0;
      1    1    0    0    0       -3.0          0.0        0.0           0.0;
      1   -1    0   -1    0       -3.0          0.0        0.0           0.0;
     -1   -1    2    2    2       -3.0          0.0        1.0           0.0;
      0   -1    2    2    2       -3.0          0.0        1.0           0.0;
      1   -1    2    0    2       -3.0          0.0        1.0           0.0;
      3    0    2    0    2       -3.0          0.0        1.0           0.0;
     -2    0    2    0    2       -3.0          0.0        1.0           0.0;
      1    0    2    0    0        3.0          0.0        0.0           0.0;
     -1    0    2    4    2       -2.0          0.0        1.0           0.0;
      1    0    0    0    2       -2.0          0.0        1.0           0.0;
     -1    0    2   -2    1       -2.0          0.0        1.0           0.0;
      0   -2    2   -2    1       -2.0          0.0        1.0           0.0;
     -2    0    0    0    1       -2.0          0.0        1.0           0.0;
      2    0    0    0    1        2.0          0.0       -1.0           0.0;
      3    0    0    0    0        2.0          0.0        0.0           0.0;
      1    1    2    0    2        2.0          0.0       -1.0           0.0;
      0    0    2    1    2        2.0          0.0       -1.0           0.0;
      1    0    0    2    1       -1.0          0.0        0.0           0.0;
      1    0    2    2    1       -1.0          0.0        1.0           0.0;
      1    1    0   -2    1       -1.0          0.0        0.0           0.0;
      0    1    0    2    0       -1.0          0.0        0.0           0.0;
      0    1    2   -2    0       -1.0          0.0        0.0           0.0;
      0    1   -2    2    0       -1.0          0.0        0.0           0.0;
      1    0   -2    2    0       -1.0          0.0        0.0           0.0;
      1    0   -2   -2    0       -1.0          0.0        0.0           0.0;
      1    0    2   -2    0       -1.0          0.0        0.0           0.0;
      1    0    0   -4    0       -1.0          0.0        0.0           0.0;
      2    0    0   -4    0       -1.0          0.0        0.0           0.0;
      0    0    2    4    2       -1.0          0.0        0.0           0.0;
      0    0    2   -1    2       -1.0          0.0        0.0           0.0;
     -2    0    2    4    2       -1.0          0.0        1.0           0.0;
      2    0    2    2    2       -1.0          0.0        0.0           0.0;
      0   -1    2    0    1       -1.0          0.0        0.0           0.0;
      0    0   -2    0    1       -1.0          0.0        0.0           0.0;
      0    0    4   -2    2        1.0          0.0        0.0           0.0;
      0    1    0    0    2        1.0          0.0        0.0           0.0;
      1    1    2   -2    2        1.0          0.0       -1.0           0.0;
      3    0    2   -2    2        1.0          0.0        0.0           0.0;
     -2    0    2    2    2        1.0          0.0       -1.0           0.0;
     -1    0    0    0    2        1.0          0.0       -1.0           0.0;
      0    0   -2    2    1        1.0          0.0        0.0           0.0;
      0    1    2    0    1        1.0          0.0        0.0           0.0;
     -1    0    4    0    2        1.0          0.0        0.0           0.0;
      2    1    0   -2    0        1.0          0.0        0.0           0.0;
      2    0    0    2    0        1.0          0.0        0.0           0.0;
      2    0    2   -2    1        1.0          0.0       -1.0           0.0;
      2    0   -2    0    1        1.0          0.0        0.0           0.0;
      1   -1    0   -2    0        1.0          0.0        0.0           0.0;
     -1    0    0    1    1        1.0          0.0        0.0           0.0;
     -1   -1    0    2    1        1.0          0.0        0.0           0.0;
      0    1    0    1    0        1.0          0.0        0.0           0.0;]

################################################################################
#                                  Functions
################################################################################

"""
    nutation_fk5(JD_TT::Number, n_max::Number = 106, nut_coefs_1980::Matrix = nut_coefs_1980)

Compute the nutation parameters at the Julian Day `JD_TT` [Terrestrial Time]
using the 1980 IAU Theory of Nutation. The coefficients are `nut_coefs_1980`
that must be a matrix in which each line has the following syntax **[1]**(p.
1043):

    an1  an2  an3  an4  an5  Ai  Bi  Ci  Di

where the units of `Ai` and `Ci` are [0.0001"] and the units of `Bi` and `Di`
are [0.0001"/JC]. The user can also specify the number of coefficients `n_max`
that will be used when computing the nutation. If `n_max` is omitted, the it
defaults to 106.

# Returns

* The mean obliquity of the ecliptic [rad].
* The nutation in obliquity of the ecliptic [rad].
* The nutation in longitude [rad].

# References

- **[1]**: Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
    Microcosm Press, Hawthorn, CA, USA.
"""
function nutation_fk5(JD_TT::Number, n_max::Number = 106, nut_coefs_1980::Matrix = nut_coefs_1980)
    # Check inputs.
    if n_max > 106
        @warn("The maximum number of coefficients to compute nutation using IAU-76/FK5 theory is 106.")
        n_max = 106
    elseif n_max <= 0
        @warn("n_max must greater than 0. The default value will be used (106).")
        n_max = 106
    end

    # Compute the Julian Centuries from `JD_TT`.
    T_TT = (JD_TT - JD_J2000) / 36525

    # Auxiliary variables
    # ===================

    d2r = π / 180

    # Mean obliquity of the ecliptic
    # ==============================

    # Compute the mean obliquity of the ecliptic [°].
    mϵ_1980 = @evalpoly(T_TT, 23.439291, -0.0130042, -1.64e-7, +5.04e-7)

    # Reduce to the interval [0, 2π]°.
    mϵ_1980 = mod(mϵ_1980, 360) * d2r

    # Delaunay parameters of the Sun and Moon
    # =======================================

    # Evaluate the Delaunay parameters associated with the Moon and the Sun
    # in the interval [0,2π]°.
    #
    # The parameters here were updated as stated in the errata [2].
    r = 360

    M_m = @evalpoly(
        T_TT,
        +134.96298139,
        +(1325r + 198.8673981),
        +0.0086972,
        +1.78e-5
    )
    M_m = mod(M_m, 360) * d2r

    M_s = @evalpoly(
        T_TT,
        +357.52772333,
        +(99r + 359.0503400),
        -0.0001603,
        -3.3e-6
    )
    M_s = mod(M_s, 360) * d2r

    u_Mm = @evalpoly(
        T_TT,
        +93.27191028,
        +(1342r + 82.0175381),
        -0.0036825,
        +3.1e-6
    )
    u_Mm = mod(u_Mm, 360) * d2r

    D_s = @evalpoly(
        T_TT,
        +297.85036306,
        +(1236r + 307.1114800),
        -0.0019142,
        +5.3e-6
    )
    D_s = mod(D_s, 360) * d2r

    Ω_m = @evalpoly(
        T_TT,
        +125.04452222,
        -(5r + 134.1362608),
        +0.0020708,
        +2.2e-6
    )
    Ω_m = mod(Ω_m, 360) * d2r

    # Nutation in longitude and obliquity
    # ===================================

    # Compute the nutation in the longitude and in obliquity.
    ΔΨ_1980 = 0.0
    Δϵ_1980 = 0.0

    @inbounds for i = 1:n_max
        # Unpack values.
        an1 = nut_coefs_1980[i,1]
        an2 = nut_coefs_1980[i,2]
        an3 = nut_coefs_1980[i,3]
        an4 = nut_coefs_1980[i,4]
        an5 = nut_coefs_1980[i,5]
        Ai  = nut_coefs_1980[i,6]
        Bi  = nut_coefs_1980[i,7]
        Ci  = nut_coefs_1980[i,8]
        Di  = nut_coefs_1980[i,9]

        a_pi = an1 * M_m + an2 * M_s + an3 * u_Mm + an4 * D_s + an5 * Ω_m

        sin_a_pi, cos_a_pi = sincos(a_pi)

        ΔΨ_1980 += (Ai + Bi * T_TT) * sin_a_pi
        Δϵ_1980 += (Ci + Di * T_TT) * cos_a_pi
    end

    # The nutation coefficients in `nut_coefs_1980` lead to angles with unit
    # 0.0001". Hence, we must convert to [rad].
    ΔΨ_1980 *= 0.0001 / 3600 * d2r
    Δϵ_1980 *= 0.0001 / 3600 * d2r

    # Return the values.
    return mϵ_1980, Δϵ_1980, ΔΨ_1980
end