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

using Markdown
using InteractiveUtils

# ╔═╡ 9b7807e6-dc3e-11eb-0f5d-c9b81d4cf5b5
begin
	import Pkg
	Pkg.activate()

	using Revise
	using AstroTime, Measurements, EarthOrientation, Plots, BenchmarkTools
	using PlutoUI
	import Dates, ERFA, PyPlot
	PyPlot.svg(true)
	EarthOrientation.update()
end

# ╔═╡ 96f6bc01-e621-46aa-a3e5-57015b9c4e5b
html"<button onclick='present()'>present</button>"

# ╔═╡ d0fc98d1-d261-4a0b-b610-4819b4b56870
md"""
# A Short History of AstroTime.jl
Hello JuliaCon, my name is Helge Eichhorn. I am a software engineer at Telespazio Germany where I work on Space Mission Analysis and Design Tools and today I am going to tell you the short history of AstroTime.jl. First of all..."""

# ╔═╡ 64b85a57-1b94-4310-8c21-11e58bda1d9e
md"""
## What is AstroTime.jl?
In short, AstroTime.jl is an MIT-licensed library which gives you a high-accuracy `DateTime`-like data type that also supports different time scales. You might ask yourself...
"""

# ╔═╡ 4b00a62c-6d5f-413a-a8f7-8c87675ffc29
md"""
## What does that mean? Why would I even want different time scales?
Excellent question! But first, let's go back to the beginning.

[Stock footage - Big Bang CGI]

Okay, I meant the beginning of AstroTime.jl not the beginning of time.
But on the other hand, why not?

We will use the time at which I am recording this presentation as the starting point.
"""

# ╔═╡ 6bd1ccfa-f269-4f9d-8683-63bc10456d8a
today = now(TAIEpoch)

# ╔═╡ 97128ecc-2a60-4afe-b1e1-0cb76829129b
md"""
We will then define the age of the universe which is around 13.772 billion years.
"""

# ╔═╡ 0cc58f90-6284-4de4-92b5-3d56a85d2254
age_of_the_universe = 13.772e9years

# ╔═╡ 627bd6d6-ad82-44c4-811d-fdb095213844
md"""
So finally the date of the big bang is `today` minus the age of the universe.
"""

# ╔═╡ 4deefa9c-a50b-4bff-8c74-6c65a7864b76
big_bang = today - age_of_the_universe

# ╔═╡ f269b1c7-c924-4764-9f68-d0d5f2879705
md"""
This calendrical representation is of course not very meaningful on such a large time scale. I will talk more about calendar troubles later. Let's have a look at the Julian Date instead which is a continuous count of days since a specific epoch.
"""

# ╔═╡ d9724b53-f119-44fa-a4ba-3a7bd51f35c7
j2000(big_bang)

# ╔═╡ f377dd26-723d-4012-b1bf-26767e2194c1
md"""
This means the Big Bang occured roughly 5 trillion days before January 1st of the year 2000.

Claiming that the Big Bang occured at that exact moment might also be a bit of a stretch. Let's add some uncertainty with `Measurements.jl` then. More appropriately, the approximate age of the universe is 13.772 billion plus-minus 0.04 billion years.
"""

# ╔═╡ 25eea271-8589-497e-9279-7ed69beb6b2e
approximate_age_of_the_universe = (13.772e9 ± 0.040e9) * years

# ╔═╡ 06aa7a0d-fc27-4051-b3d5-923ef4f9e07b
md"""
With that we can calculate an approximate date for the Big Bang:
"""

# ╔═╡ 88d03e9b-560c-4b27-9a6b-2d6aeaf49946
approximate_big_bang = today - approximate_age_of_the_universe

# ╔═╡ cda813db-2975-491a-9cde-173ce82b2ad8
md"""
The calendar representation is still not helpful but the Julian Date reveals the uncertainty.
"""

# ╔═╡ e896ac12-0872-41c7-8073-da01f98ccb65
j2000(approximate_big_bang)

# ╔═╡ dc3ea266-258e-4f2e-8f85-eccfcd3ee4a2
md"""
Let's use the overly confident "exact" Big Bang date for another experiment. According to cosmologists, baryons such as protons and neutrons started forming at about $$10^{-6}$$ seconds after the Big Bang.
"""

# ╔═╡ 4b9b0c5d-bf85-4e4f-ae9e-ebc6abe68bdb
baryons = big_bang + 1e-6seconds

# ╔═╡ 87343448-d6bc-4256-bf7b-1a66945ceaaf
md"""
If we now take this timestamp and add the age of the universe and then substract today, we get back the same $$10^{-6}$$ seconds of elapsed time without any loss in accuracy.
"""

# ╔═╡ 4d8b385a-3b39-4e10-970b-223b96b2238d
baryons + age_of_the_universe - today

# ╔═╡ 2ebda195-3668-4da7-940c-ec9ab90537ad
md"""
Alright, this was a pretty silly example and I am quite sure that this is not how cosmological time works. So let's go forward several billion years to the beginning of AstroTime.jl.
"""

# ╔═╡ f3c30b90-a8df-434d-b522-f1f37ceb7444
md"""
## Part 1 – The Beginning

Whether you are tracking celestial phenomena or you are navigating a spacecraft to another world, you need accurate means of time keeping. Even if your needs are more down to Earth, in every field you will need the ability to time events and calculate the elapsed time between them.

At this point everybody realises that doing calculations with our western Gregorian calendar is absolutely insane. Months have different lengths by default and then there is the trouble with leap years.
The fundamental problem is that the Earth needs approximately 365.2422 days to orbit the Sun. By inserting leap years with 366 days the Gregorian calendar aims to make the average year 365.2425 days long to keep the calendar synchronized with Earth's actual revolution. Here is the rule:

> Every year that is exactly divisible by four is a leap year, except for years that
> are exactly divisible by 100, but these centurial years are leap years if they are
> exactly divisible by 400. For example, the years 1700, 1800, and 1900 are not leap
> years, but the years 1600 and 2000 are.[2]
>    — United States Naval Observatory

Messy.

It get's even more complicated if you need to consider historic events. If they occured before the introduction of the Gregorian calendar they are probably given in the preceding Julian calendar but they might also use a proleptic Gregorian calendar which is extended into the past.

For the record, AstroTime.jl uses the Julian calendar for dates preceding the introduction of the Gregorian calendar and a proleptic Julian calendar for dates before January 1, 1 AD. The `Dates` standard library on the other hand uses a proleptic Gregorian calendar.
"""

# ╔═╡ aee5aaa1-147d-4aa0-9704-d9894c48da28
AstroTime.AstroDates.calendar(AstroTime.Date(-753, 4, 21))

# ╔═╡ 77112dd1-df7e-4d22-83ab-3bc2de031383
AstroTime.AstroDates.calendar(AstroTime.Date(1415, 10, 25))

# ╔═╡ 7a0fe852-ee76-44f7-9f38-e80d40063120
AstroTime.AstroDates.calendar(AstroTime.Date(2012, 2, 14))

# ╔═╡ 339906f0-191f-4c3e-ab04-6fd6f724f4ce
md"""
A solution to this has been independently developed by multiple communities. Instead of using a calendar you just count specified time intervals since a start epoch.

One example is Unix time which counts seconds since Midnight on January 1st, 1970.

Another are Julian Dates which I have mentioned in the introductory example and are used by astronomers and astronautical engineers.

**Note:** Julian Dates should not be confused with the Julian calendar.
"""

# ╔═╡ cc006259-cac2-43ca-a77c-1b1a0ce6be63
Dates.datetime2unix(Dates.DateTime(today))

# ╔═╡ b3e8edf7-1b2f-4244-81b5-7ac9473c2e78
julian(today)

# ╔═╡ 55b4aa58-634e-4f79-853e-876a0f2156af
modified_julian(today)

# ╔═╡ f8353359-7638-420b-bf66-54a3df8565dc
j2000(today)

# ╔═╡ 9ed6cdae-902d-4233-95bf-d6b54bbd697f
md"""
They count days starting at noon on Monday, January 1, 4713 BC, in the proleptic Julian calendar. There are other variants that start at different epochs such as the aforementiond `J2000` which starts at January 1, 2000.

Why do they start at noon? Because the date shift should not be in the middle of the night when astronomers are out making their observations.

The advantages for doing math are clear. If you want to determine the time interval between two dates, you can directly subtract the Julian Dates from another and you do not need to care about whether or not there is a leap day in between.
"""

# ╔═╡ 29e00530-aa7a-4318-9b79-6b54f287d5ed
jd1 = julian(TAIEpoch(2020, 2, 28))

# ╔═╡ ce6168f3-d620-4e5e-8e46-51bea8cceb8c
jd2 = julian(TAIEpoch(2020, 3, 1))

# ╔═╡ 681586c1-0ec4-43b2-aebe-8d8004d92bb2
jd2 - jd1

# ╔═╡ 5efd6405-bf02-4dc1-87f5-42c4851b1b7c
md"""
The major disadvantage is also obvious. Julian Dates and also Unix timestamps are unparseable by most humans because they lack the calendarial structure that we are used to.

In the end, you need both.

Thus, my first experiments and incidentally also some of my first attempts with Julia while I was an intern at the ESA's European Space Operations Centre, tried to combine these two representations. This was back in 2013 before the introduction of the excellent `Dates` standard library. As you can see by looking at the code, I was not really certain about the data model and whether or not this was a good idea. *Spoiler alert:* It wasn't.
"""

# ╔═╡ bb827858-9f12-470e-86e6-0771c0524797
Show(MIME"image/png"(), read("figures/code1.png"))

# ╔═╡ 13d820b6-83f8-4d23-b894-d22d2413d00b
md"""
At this point in time, I was heavily invested in the Scientific Python ecosystem and still on the edge about Julia. I was certain that Julia would be an excellent fit for space science and engineering but humbled by the magnitude of the task of building such an ecosystem from scratch I started a PhD in a different field (*spoiler alert:* also not a great idea), worked on Python things, and procrastinated for a few years.
"""

# ╔═╡ 5a2ca7a1-13b6-45e7-8e8d-81c5efb92fdd
md"""
## Part 2 – Time Scales

Fast forward to early 2016. Two other guys and myself had just tried to merge our Python-based astrodynamics libraries into the one true library to rule them all. After a promising start we got hampered by clashing design sensibilities, bikeshedding about object hierarchies, and, what bothered me the most, massive performance problems everytime we tried to introduce an abstraction that went beyond a NumPy array. As a result we lost steam and the project was gradually abandoned.

Based on this experience, I decided to go all in with Julia. This time I did not try to write everything from scratch but tried to reuse as much existing code and know-how from other languages as possible. One of the nicest parts of working with Python was `Astropy`.

I was used to dealing with Julian Dates in the form of raw floating point values but Astropy provides a nice abstraction for dealing with all sorts of astronomical dates in form of the `Time` class.

```python
import numpy as np
from astropy.time import Time

times = ['1999-01-01T00:00:00.123456789', '2010-01-01T00:00:00']
t = Time(times, format='isot', scale='utc')
```

Building something as powerful in Julia would be the first step. But before we can talk about the next iteration of the data model, we need to discuss one of the initial questions: Why do we need different time scales at all?

According to Newtonian mechanics we only need to measure hard enought to find the one true time:

> Absolute, true and mathematical time, of itself, and from its own nature flows equably without regard to anything external, [...]
>	– Isaac Newton, Principia (1687)

This turned out to be a wicked problem. But 300 years later we are approaching this ideal with atomic clocks. An optical lattice strontium clock will gain or lose less than a second in 15 billion years which is longer than the age of the universe, as previously discussed.

International Atomic Time or TAI (from French "temps atomique international"), the weighted average of over 400 atomic clocks worldwide, gives us a high-precision time standard and our civil Coordinated Universal Time or UTC is also based on it. We also have global satellite navigation constellations such as GPS and Galileo which are essentially atomic clocks flying in space. This means that everybody with a modern smartphone has instant access those highly accurate clocks.

All is well? Not quite. Our daily lives are still governed by solar time and the passing of day and night which is defined by the rotation of the Earth. Unfortunately, Earth's rotation is much less stable than an atomic clock.

The main time scale in use based on Earth's rotation is Universal Time and more specifically UT1. The difference between TAI and UT1 is a measured quantity and one of the Earth Orientation Parameters published by the International Earth Rotation and Reference Systems Service (IERS). We can access these parameters with `EarthOrientation.jl`.
"""

# ╔═╡ c00d5a3e-54b9-483a-96b8-efc2bc4a8667
ep_rng = TAIEpoch(1975,1,1):1days:today

# ╔═╡ 49ae6271-8d46-4b2c-be77-7bc0c648091f
ΔUT1_TAI = getΔUT1_TAI.(value.(julian.(ep_rng)));

# ╔═╡ f8d73ace-06fe-4279-9e96-60ac72d6b71c
begin
	PyPlot.plot(Dates.DateTime.(ep_rng), ΔUT1_TAI, "tab:blue")
	PyPlot.gcf()
end

# ╔═╡ 5a609cc3-d9c3-45dd-b1cf-622da4f80146
md"""
We can see that TAI is ahead of UT1 because Earth's rotation is getting slower. This means that without countermeasures UTC which is based on atomic time would drift away from UT1 long term.

This is prevented by inserting leap seconds into the UTC time scale which keeps the difference between UTC and UT1 under 0.9 seconds.

If we plot this difference over time, the jumps become clearly visible.
"""

# ╔═╡ 83fdd3d8-aa97-4f26-bf55-44daf62af476
ΔUT1 = getΔUT1.(value.(julian.(ep_rng)));

# ╔═╡ e6b570fa-a15e-4510-a282-c71e4b590668
begin
	PyPlot.clf()
	PyPlot.plot(Dates.DateTime.(ep_rng), ΔUT1, "tab:blue")
	PyPlot.gcf()
end

# ╔═╡ c84dc04c-7f43-4f5b-b07a-1aaed19b7d10
md"""
I'll do a little bit of foreshadowing by saying that leap seconds have been the bane of my existence for the past few years and we will talk about those troubles later.

In summary, we need these three time scales to accurately deal with time on Earth. How does the picture change if we go to outer space?

In the previous discussion, I omitted an important detail about TAI. Namely that it is based on the passage of "proper time" on  Earth's geoid, an idealized representation of Earth's surface.

Why does that matter? Because thanks to Albert Einstein's theory of relativity and subsequent experiments proving it we know today that there is no absolute Newtonian time. Time is relative.

Due to gravitational time dilation a clock that is closer to a source of gravity will tick slower than a clock that is farther away. There is also time dilation due to relative velocity. A clock that is moving relative to a clock at rest will be measured to tick slower.

If you think about the implications of this when moving around the Solar System it gets complicated pretty fast.

To deal with these relativistic issues astronomers have defined additional time scales and the transformations between them:

- **Terrestrial Time (TT):** For observations from Earth's surface. Identical to TAI except for a fixed offset.
- **Geocentric Coordinate Time (TCG, from "Temps-Coordonnée Géocentrique"):** The time experienced by a clock at rest in a coordindate frame moving with the center of the Earth but unaffected by Earth's gravity and rotation.
- **Barycentric Coordinate Time (TCB, from "Temps-Coordonnée Barycentrique"):** The time experienced by a clock at rest in a coordindate frame moving with the center of mass of the Solar System but unaffected by the gravity of the Solar System.
- **Barycentric Dynamical Time (TDB, from "Temps Dynamique Barycentrique"):** Which is defined as a linear scaling from TCB and takes into account time dilation.

Thus, my goal for what would become AstroTime.jl was set, deal with the zoo of astronomical time scales in a human-friendly manner while taking inspiration from Astropy.

Internally, Astropy uses the ERFA C library which is an open-source derivative of the Standards of Fundamental Astronomy (SOFA) library published by the International Astronomical Union for most of its date and time functionality. Thanks to Mike Nolta and the JuliaAstro community, ERFA.jl, a Julia wrapper for ERFA, was available at this time and I decided to take the same approach.

ERFA models dates as two-part Julian dates. This means that a date consists of two arbitrarily split `Float64` values so that `jd = jd1 + jd2`. For example, by splitting between the integral day number and the fractional part of day.
"""

# ╔═╡ 127301a1-8d09-4881-a4cd-bdeca5cd3a8f
utc1, utc2 = ERFA.dtf2d("UTC", 2018, 2, 6, 20, 45, 0.0)

# ╔═╡ 4723271f-a96f-4fa9-aba3-ccb9207561db
md"""
This approach increases the accuracy because current Julian Dates are large numbers and thus more digits are available for the fractional part if the day number has been subtracted.

This was how my initial rewrite looked like. As you can see, I had just discovered metaprogramming. A dangerous time for a young padawan...
"""

# ╔═╡ dcd8c12a-3e80-41c7-8467-1a6538c95616
Show(MIME"image/png"(), read("figures/code2.png"))

# ╔═╡ 204dd602-470d-4952-9909-b115724c8620
md"""
While none of the code from this era has survived until the present day, an important idea that endured crystallized during this time.

**The computer should be in charge of keeping track of time scales.**

In the classical ERFA approach this is the user's job:
"""

# ╔═╡ 95b02b29-d2db-427e-b535-0c84dcfc493f
tai1, tai2 = ERFA.utctai(utc1, utc2)

# ╔═╡ 706e344d-4750-4536-8b2a-c82e749f581d
tt1, tt2 = ERFA.taitt(tai1, tai2)

# ╔═╡ e6333e26-0e8b-4191-9517-7bfb5df4c296
tcg1, tcg2 = ERFA.tttcg(tt1, tt2)

# ╔═╡ 50ff0c6d-4e59-48a7-8af4-2eac7cc250fe
md"""
This means that the user needs to be aware of the relationships between time scales.
"""

# ╔═╡ 3a814113-a528-4702-a3c0-728c9a5212ba
Show(MIME"image/png"(), read("figures/graph.png"))

# ╔═╡ 81d2aa16-a205-4f8c-a2ec-4ec0a35db927
md"""
Astropy's `Time` class already abstract this away but AstroTime.jl takes it one step further by leveraging Julia's expressive type system. All iterations of the `Epoch` type that followed have used the time scale as a type parameter:

```julia
abstract type TimeScale end

struct Epoch{T<:TimeScale}
	scale::T
	...
end
```

This enables the Julia compiler to figure out the necessary transformations between time scales automatically and we can directly convert to TCG without any intermediary steps.
"""

# ╔═╡ 2b32d793-3122-4a6c-a184-02ed3e3ff43e
tai = from_utc(2018, 2, 6, 20, 45, 0.0, scale=TAI)

# ╔═╡ c296abea-4511-4933-897d-acd597892221
typeof(tai)

# ╔═╡ 65dcf800-ede6-44db-8085-520549a8deb9
tcg = TCGEpoch(tai)

# ╔═╡ 93287052-6673-4245-b691-124d553e50f3
typeof(tcg)

# ╔═╡ e6d1c9ce-fab9-47b7-8dcc-844191d5489b
value(julian(tcg)) ≈ tcg1 + tcg2

# ╔═╡ 542b1bd3-39b4-4880-aef6-2aca13e3e420
md"""
This simple API has proven to be very powerful and has remained unchanged till today.

However, at this point the `Epoch` code was still part of `Astrodynamics.jl` which I presented at my first JuliaCon in Boston in 2016. Shortly thereafter I decided to put the time code into a separate `AstronomicalTime.jl` package which led to me being contacted by Kyle Barbary and Mosé Giordano from JuliaAstro and joining forces with them.

`AstronomicalTime.jl` was moved to `JuliaAstro` and renamed to `AstroTime.jl`. I also quit my job in academia and returned to the space industry.

With a solid foundation based on ERFA in place, we decided to start porting the relevant parts of ERFA from C to Julia. Our GSOC student Prakhar Srivastava took care of that and more in the summer of 2018.

All was well! Or was it?
"""

# ╔═╡ f419a06d-3dbc-41ab-847d-f7dc2a1f74ad
md"""
## Part 3 – Standing on the Shoulders of Giants

In mid 2018, my colleague Bernard Godard who is a Flight Dynamics Engineer at ESOC discovered the library and did some experiments. He used the calculation of transmission times for an interplanetary spacecraft such as Voyager 2 where there is a considerable amount of light time delay because of the large distances as an example.
"""

# ╔═╡ c77f4560-f501-45d8-b627-cb9f743125ec
reception_time = TDBEpoch("2021-07-01T00:00:00.00")

# ╔═╡ da3560e2-366e-4264-ada3-5fd100f0e78d
# round trip light time
rtlt_a = seconds(1.5days)

# ╔═╡ a128db96-ef58-474d-8815-b1015337dfb3
rtlt_b = rtlt_a + 1e-6seconds

# ╔═╡ 8a5f595c-8ac2-4f93-8e99-03bed6429387
velocity = 15.341 # km/s

# ╔═╡ b400d692-e232-4852-8e07-9a6b96bfd938
distance_travelled = value(seconds(rtlt_b - rtlt_a)) * velocity # km

# ╔═╡ d0abad31-b745-4fd9-9725-2ca2544c177d
transmission_time_a = reception_time + rtlt_a

# ╔═╡ b2230ea5-bf63-44f1-babe-0684207871fe
transmission_time_b = reception_time + rtlt_b

# ╔═╡ 15bfe97a-4ef7-4fb7-9463-574385c068c8
transmission_time_b - transmission_time_a

# ╔═╡ ee3cfe44-40f3-4cd6-942a-35cc0ac81381
md"""
The current version of AstroTime.jl preserves the time difference accurately. But back then the result was different...
"""

# ╔═╡ e80e6b32-c6f5-4ad1-bb13-0807686212ce
Show(MIME"image/png"(), read("figures/issue.png"))

# ╔═╡ b387b901-4ce1-4022-995d-999f5b000d4e
md"""
Since deep space exploration is one of my favorite topics and Julia shall become a multiplanetary programming language one day, this was not acceptable. Back to the drawing board...

After some discussions, we decided to adopt the same high-accuracy internal representation of dates that the open-source Java flight dynamics library Orekit uses. The result looked somewhat like this.

```julia
struct Epoch{T<:TimeScale}
	scale::T
	tai_second::Int64 # Integral seconds since 2000-01-01T00:00:00.00 TAI
	fraction::Float64 # Fraction of the current second
end
```

Which meant another rewrite of the library was in order. The first version was a straight port of the `AbsoluteDate` class and other auxiliary classes from Java to Julia.

The name of that class should give one pause though. The trick the Orekit authors apply is to convert all dates internally to the TAI time scale. This makes it easier to deal with the conversions between scales because you are dealing with 1-to-N conversions and not M-to-N conversions and it also avoids troubles with leap seconds because TAI is a nice and uniform time scale. It is also a reasonable thing to do if you are dealing with spacecraft in orbit around Earth.

Not so much for deep space missions, though, as Bernard informed me after I presented my rewrite to him. The problem boils down to gravitational time dilation again. If we consider a spacecraft far away from Earth, it will be in a completely different part of the Solar System's gravity well and experience a different gravitational potential than a spacecraft close Earth. Thus, there would be no straightforward relation between the proper relativistic time experienced by the spacecraft's clock and TAI on Earth.

Turns out an `AbsoluteDate` does not really work when everything is relative.

Another iteration later in early 2020, we had successfully moved back to relative transformations.

```julia
struct Epoch{T<:TimeScale}
	scale::T
	second::Int64 # Integral seconds since 2000-01-01T00:00:00.00 in `scale`
	fraction::Float64 # Fraction of the current second
end
```

And then leap seconds reared their ugly head again...

The problem with UTC is that there is a discontinuity in the time scale everytime a leap second is introduced. We learned the hard way that this created ambiguities in our data model during around leap seconds. Let me show you the problem:

| TAI | UTC | TAI second | Offset | UTC second |
|:----|:-----|:------|:------------|:--------------|
| 2017-01-01T00:00:35.000 | 2016-12-31T23:59:59.000 | 536500835 | 36 | 536500835 - 36 = 536500799 |
| 2017-01-01T00:00:36.000 | 2016-12-31T23:59:60.000 | 536500836 | **36/37** | **536500836 - 36 = 536500800 / 536500836 - 37 = 536500799** |
| 2017-01-01T00:00:37.000 | 2017-01-01T00:00:00.000 | 536500837 | 37 | 536500837 - 37 = 536500800 |
| 2017-01-01T00:00:38.000 | 2017-01-01T00:00:01.000 | 536500838 | 37 | 536500838 - 37 = 536500801 |

At this point, we had two choices: change the storage format once again which would incur another complete rewrite or relegate UTC to be an I/O format.

We chose the latter. Sorry UTC, but from now on this is a uniform-time-scales-only club...
"""

# ╔═╡ fac44891-7af3-447a-b930-fd1711a096f0
ep = from_utc(2018, 2, 6, 20, 45, 0.0) # This used to be `UTCEpoch(...)`

# ╔═╡ 5cfb21ed-dc03-4c37-94ee-bb8e3436dca3
typeof(ep)

# ╔═╡ a30edc9f-7f22-4bcd-9623-31e20e3ee08a
to_utc(String, ep)

# ╔═╡ 9badb9ee-9a46-4a62-b70d-616679f47af4
md"""
## Part 4 – Seeing Further

This is where we are today. We would not have been able to get this far without open source because we took inspiration and concrete implementation details from many different sources.

> If I have seen further it is by standing on the shoulders of Giants.
> – Isaac Newton

I want to use this part of the talk to highlight these giants and show where we have been able to make improvements because of the headstart they gave us.

I also want to give you ideas how you could use AstroTime.jl even if your goal is not to explore distant worlds.

### Astropy & Orekit

I have already talked at length about how Astropy and Orekit served as inspirations and starting points for AstroTime.jl's implementation. The founders of both projects are pioneers of the space science open-source community and doing it a great service.

Thanks to Julia, we have been able to surpass these projects in terms of performance. The following conversion does four separate conversion steps in a matter of a few hundred nanoseconds.
"""

# ╔═╡ bc10b3a4-f3ff-41d6-a4a4-1fc472015fd0
@benchmark TCBEpoch($today)

# ╔═╡ 1b1b43ed-893f-4e09-a289-698ab66d3107
md"""
As of this month, we have also been able to improve upon the accuracy of Orekit, which is another great milestone.

### AccurateArithmetic.jl

We achieved this by using error-free transforms wherever we can and also by tracking the associated floating point error.

```julia
struct Epoch{S<:TimeScale, T} <: Dates.AbstractDateTime
    scale::S
    second::Int64
    fraction::T
    error::T
    function Epoch{S}(second::Int64, fraction::T, error::T=zero(T)) where {S<:TimeScale, T<:AbstractFloat}
        return new{S, T}(S(), second, fraction, error)
    end
end
```

The algorithms were adapted from AccurateArithmetic.jl.

### Chrono.jl

You have seen me use the shorthand syntax for time periods.
"""

# ╔═╡ 83050f51-b7bf-48c0-8e04-21e4c1b0a354
today + 5seconds

# ╔═╡ 5d794ce9-6a63-404d-b588-d4795a8af199
md"""
I adopted this idea from the [Chrono.jl](https://github.com/FugroRoames/Chrono.jl) prototype by Andy Ferris and Chris Foster.

### Dates

AstroTime.jl is reusing all of the parsing and formatting machinery of the excellent `Dates` standard library and is generally highly compatible with it.

You can easily convert from `DateTime` to `Epoch` and vice-versa. Although, please note that `DateTime` only supports millisecond accuracy.
"""

# ╔═╡ 3ff48e3d-f4c1-4e4f-8d2e-f5b94f084366
ep_from_dt = TAIEpoch(Dates.DateTime(2018, 2, 6, 20, 45, 0, 0))

# ╔═╡ 7bd69f4b-fc08-42e9-a9b2-2177b55574fe
dt_from_ep = Dates.DateTime(ep_from_dt)

# ╔═╡ a50a4abf-9f2f-4937-a95d-70792f1698e0
md"""
There is currently an [issue](https://github.com/JuliaLang/julia/issues/37579) limiting the parser from `Dates` to milliseconds even though `Dates.Time` does support nanoseconds.

We have a workaround for that...
"""

# ╔═╡ fbb110b9-9651-4c66-ad5d-9bfddefaa2f2
Dates.Time(TAIEpoch("2018-02-06T20:45:00.123456789"))

# ╔═╡ cfdc0104-2d25-4d62-9056-fa0a4ae20871
md"""
We have also added a few other small features that you might find useful.

For example, the day-of-year format for timestamps:
"""

# ╔═╡ ff63332f-13b6-4556-a1ee-0a3886344034
doy_ep = TAIEpoch("2018-37T20:45:00", Dates.dateformat"yyyy-DDDTHH:MM:SS.fff")

# ╔═╡ 345ec872-d1f0-4c58-87d9-22b3fd6f07a7
md"""
Or the aforementioned support for the Julian and proleptic Julian calendars.

We hope to be able to upstream a few of these improvements to `Dates` in the future.

To cut a long story short: AstroTime.jl should be able to serve as a drop-in replacement for Dates if you need higher resolution, deal with leap seconds or want to fly a spacecraft to Jupiter. And if not, I want to know...
"""

# ╔═╡ c3f5eeeb-2a6c-43be-b53b-ab435a45baa5
md"""
## Part 5 – The Undiscovered Country

Where do we go from here?

As of now, I consider AstroTime.jl to be fairly complete and fit for purpose. We have been using it in production for a while. Unfortunatly not for flying spacecraft to Jupiter but at least for designing missions to the Moon.

Even though there has been a lot of turmoil about the internal representations over the last few years, the public API has been remarkably stable. Thus, I am planning for the next major version to be 1.0.

Finally, let me say thank you for your attention and I hope that you found this presentation relatively interesting (pun intended). I encourage you to give AstroTime.jl a try by following the tutorial and you can find me on GitHub (@helgee), Discourse (@helgee), Twitter (@helge_e), and sometimes Slack and Zulip where I will be waiting for your feedback.

Thanks again and I hope to see you all in person next year.
"""

# ╔═╡ Cell order:
# ╠═9b7807e6-dc3e-11eb-0f5d-c9b81d4cf5b5
# ╟─96f6bc01-e621-46aa-a3e5-57015b9c4e5b
# ╟─d0fc98d1-d261-4a0b-b610-4819b4b56870
# ╟─64b85a57-1b94-4310-8c21-11e58bda1d9e
# ╟─4b00a62c-6d5f-413a-a8f7-8c87675ffc29
# ╠═6bd1ccfa-f269-4f9d-8683-63bc10456d8a
# ╟─97128ecc-2a60-4afe-b1e1-0cb76829129b
# ╠═0cc58f90-6284-4de4-92b5-3d56a85d2254
# ╟─627bd6d6-ad82-44c4-811d-fdb095213844
# ╠═4deefa9c-a50b-4bff-8c74-6c65a7864b76
# ╟─f269b1c7-c924-4764-9f68-d0d5f2879705
# ╠═d9724b53-f119-44fa-a4ba-3a7bd51f35c7
# ╟─f377dd26-723d-4012-b1bf-26767e2194c1
# ╠═25eea271-8589-497e-9279-7ed69beb6b2e
# ╟─06aa7a0d-fc27-4051-b3d5-923ef4f9e07b
# ╠═88d03e9b-560c-4b27-9a6b-2d6aeaf49946
# ╟─cda813db-2975-491a-9cde-173ce82b2ad8
# ╠═e896ac12-0872-41c7-8073-da01f98ccb65
# ╟─dc3ea266-258e-4f2e-8f85-eccfcd3ee4a2
# ╠═4b9b0c5d-bf85-4e4f-ae9e-ebc6abe68bdb
# ╟─87343448-d6bc-4256-bf7b-1a66945ceaaf
# ╠═4d8b385a-3b39-4e10-970b-223b96b2238d
# ╟─2ebda195-3668-4da7-940c-ec9ab90537ad
# ╟─f3c30b90-a8df-434d-b522-f1f37ceb7444
# ╠═aee5aaa1-147d-4aa0-9704-d9894c48da28
# ╠═77112dd1-df7e-4d22-83ab-3bc2de031383
# ╠═7a0fe852-ee76-44f7-9f38-e80d40063120
# ╟─339906f0-191f-4c3e-ab04-6fd6f724f4ce
# ╠═cc006259-cac2-43ca-a77c-1b1a0ce6be63
# ╠═b3e8edf7-1b2f-4244-81b5-7ac9473c2e78
# ╠═55b4aa58-634e-4f79-853e-876a0f2156af
# ╠═f8353359-7638-420b-bf66-54a3df8565dc
# ╟─9ed6cdae-902d-4233-95bf-d6b54bbd697f
# ╠═29e00530-aa7a-4318-9b79-6b54f287d5ed
# ╠═ce6168f3-d620-4e5e-8e46-51bea8cceb8c
# ╠═681586c1-0ec4-43b2-aebe-8d8004d92bb2
# ╟─5efd6405-bf02-4dc1-87f5-42c4851b1b7c
# ╟─bb827858-9f12-470e-86e6-0771c0524797
# ╟─13d820b6-83f8-4d23-b894-d22d2413d00b
# ╟─5a2ca7a1-13b6-45e7-8e8d-81c5efb92fdd
# ╠═c00d5a3e-54b9-483a-96b8-efc2bc4a8667
# ╠═49ae6271-8d46-4b2c-be77-7bc0c648091f
# ╠═f8d73ace-06fe-4279-9e96-60ac72d6b71c
# ╟─5a609cc3-d9c3-45dd-b1cf-622da4f80146
# ╠═83fdd3d8-aa97-4f26-bf55-44daf62af476
# ╠═e6b570fa-a15e-4510-a282-c71e4b590668
# ╟─c84dc04c-7f43-4f5b-b07a-1aaed19b7d10
# ╠═127301a1-8d09-4881-a4cd-bdeca5cd3a8f
# ╟─4723271f-a96f-4fa9-aba3-ccb9207561db
# ╟─dcd8c12a-3e80-41c7-8467-1a6538c95616
# ╟─204dd602-470d-4952-9909-b115724c8620
# ╠═95b02b29-d2db-427e-b535-0c84dcfc493f
# ╠═706e344d-4750-4536-8b2a-c82e749f581d
# ╠═e6333e26-0e8b-4191-9517-7bfb5df4c296
# ╟─50ff0c6d-4e59-48a7-8af4-2eac7cc250fe
# ╟─3a814113-a528-4702-a3c0-728c9a5212ba
# ╟─81d2aa16-a205-4f8c-a2ec-4ec0a35db927
# ╠═2b32d793-3122-4a6c-a184-02ed3e3ff43e
# ╠═c296abea-4511-4933-897d-acd597892221
# ╠═65dcf800-ede6-44db-8085-520549a8deb9
# ╠═93287052-6673-4245-b691-124d553e50f3
# ╠═e6d1c9ce-fab9-47b7-8dcc-844191d5489b
# ╟─542b1bd3-39b4-4880-aef6-2aca13e3e420
# ╟─f419a06d-3dbc-41ab-847d-f7dc2a1f74ad
# ╠═c77f4560-f501-45d8-b627-cb9f743125ec
# ╠═da3560e2-366e-4264-ada3-5fd100f0e78d
# ╠═a128db96-ef58-474d-8815-b1015337dfb3
# ╠═8a5f595c-8ac2-4f93-8e99-03bed6429387
# ╠═b400d692-e232-4852-8e07-9a6b96bfd938
# ╠═d0abad31-b745-4fd9-9725-2ca2544c177d
# ╠═b2230ea5-bf63-44f1-babe-0684207871fe
# ╠═15bfe97a-4ef7-4fb7-9463-574385c068c8
# ╟─ee3cfe44-40f3-4cd6-942a-35cc0ac81381
# ╟─e80e6b32-c6f5-4ad1-bb13-0807686212ce
# ╟─b387b901-4ce1-4022-995d-999f5b000d4e
# ╠═fac44891-7af3-447a-b930-fd1711a096f0
# ╠═5cfb21ed-dc03-4c37-94ee-bb8e3436dca3
# ╠═a30edc9f-7f22-4bcd-9623-31e20e3ee08a
# ╟─9badb9ee-9a46-4a62-b70d-616679f47af4
# ╠═bc10b3a4-f3ff-41d6-a4a4-1fc472015fd0
# ╟─1b1b43ed-893f-4e09-a289-698ab66d3107
# ╠═83050f51-b7bf-48c0-8e04-21e4c1b0a354
# ╟─5d794ce9-6a63-404d-b588-d4795a8af199
# ╠═3ff48e3d-f4c1-4e4f-8d2e-f5b94f084366
# ╠═7bd69f4b-fc08-42e9-a9b2-2177b55574fe
# ╟─a50a4abf-9f2f-4937-a95d-70792f1698e0
# ╠═fbb110b9-9651-4c66-ad5d-9bfddefaa2f2
# ╟─cfdc0104-2d25-4d62-9056-fa0a4ae20871
# ╠═ff63332f-13b6-4556-a1ee-0a3886344034
# ╟─345ec872-d1f0-4c58-87d9-22b3fd6f07a7
# ╟─c3f5eeeb-2a6c-43be-b53b-ab435a45baa5

max_stars_repo_name: helgee/juliacon

max_issues_repo_name: helgee/juliacon

max_forks_repo_name: helgee/juliacon


text:
@testset "Slice" begin
    a = [1,2,3,4]
    b = Slice(a, 2:4)
    @test b[1] == a[2]
    @test b[0] == a[1]
    b = b << -1
    @test length(b) == 4
    @test b[1] == a[1]
    b = b >> -1
    @test length(b) == 3
    @test b[3] == a[3]
    @test b[4] == a[4]
    b[4] = 10000
    @test a[4] == 10000

    @test_broken Slice(a, 3:2)
    @test_broken Slice(a, 0:4)
    @test_broken Slice(a, 1:5)
    b = Slice(a, 2:3)
    @test_broken b << -2
    @test_broken b << 2
    @test_broken b >> 2
    @test_broken b >> -2

    @test_broken b[4]
    @test_broken b[-1]
    @test_broken b[4] = 2
    @test_broken b[-1] = 2

    @test size(b) == (2,)
    @test interval(b) == 2:3

    c = Slice(b, 1:3)
    @test length(c) == 3
    @test c[3] == 10000
    c = Slice(b, 1:1)
    @test length(c) == 1
    @test b[1] == 2
    c = Slice(b, 0:3)
    @test sum(c .== a) == length(c)
    @test_broken Slice(b, -1:3)
    @test_broken Slice(b, 0:4)
    @test_broken Slice(b, -1:4)

end

ext: jl


text:
"""
Beta Uniform Mixture Model (BUM)

**Arguments**

- π0 : Contributing fraction of the uniform distribution to the full model
- α, β : Parameters of the Beta distribution, Float64, default: 0.5, 3.0

**Returns**

`MixtureModel`, as defined in the `Distributions` package, composed of

- a uniform distribution in the interval [0, 1], with weight/prior π0
- a Beta distribution with parameters α and β, with weight/prior 1-π0

**Examples**

```julia
bum = BetaUniformMixtureModel(0.2)
bum = BetaUniformMixtureModel(0.2, 0.5, 1.0)

using Distributions

rand(bum, 10)

pdf(bum, 0:0.05:1)

cdf(bum, 0:0.05:1)
```

"""
function BetaUniformMixtureModel end

function BetaUniformMixtureModel(π0::Float64, α::Float64 = 0.5, β::Float64 = 1.0)
    if !isin(π0, 0., 1.)
        throw(DomainError())
    end
    MixtureModel([Beta(α, β), Uniform()], [1-π0, π0])
end


text:
module test_91_CompareOneScalarSignalWithUnit

using ModiaResult
using Unitful
using DataFrames
ModiaResult.@usingModiaPlotPackage

t = range(0.0, stop=10.0, length=100)

result1 = DataFrame()
result2 = DataFrame()
result3 = DataFrame()

result1."time" = t*u"s"
result1."w"    = sin.(t)*u"rad/s"

result2."time" = t*u"s"
result2."w"    = 0.0005*u"rad/s" .+ sin.(t)*u"rad/s"

result3."time" = t
result3."w"    = sin.(t.+0.001)

(success2, diff2, diff_names2, max_error2, within_tolerance2) = ModiaResult.compareResults(result1, result2)
println("success2 = $success2, max_error2 = $max_error2, within_tolerance2 = $within_tolerance2")

(success3, diff3, diff_names3, max_error3, within_tolerance3) = ModiaResult.compareResults(result1, result3)
println("success3 = $success3, max_error3 = $max_error3, within_tolerance3 = $within_tolerance3")

plot(result1, "w", prefix="r1.")
plot(result2, "w", prefix="r2.", reuse=true)
plot(result3, "w", prefix="r3.", reuse=true)

plot(diff2, "w", prefix="diff2_", figure=2)

plot(diff3, "w", prefix="diff3_", figure=2, reuse=true)

end

max_stars_repo_name: ModiaSim/ModiaResult

max_issues_repo_name: ModiaSim/ModiaResult

max_forks_repo_name: ModiaSim/ModiaResult


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 1b737805-a411-4585-b215-d0f99eafac0c
using HypertextLiteral

# ╔═╡ 881b75d3-cebe-4a53-8cf5-beaeeddafd35
function skip_as_script(m::Module)
	if isdefined(m, :PlutoForceDisplay)
		return m.PlutoForceDisplay
	else
		isdefined(m, :PlutoRunner) && parentmodule(m) == Main
	end
end

# ╔═╡ dadf2f40-1764-47a4-b560-683b6479d77f
if skip_as_script(@__MODULE__)
	import Pkg
	Pkg.activate(Base.current_project(@__DIR__))
	Text("Project env active")
end

# ╔═╡ 38b7eeb9-80bb-4a3a-a2d2-809fc423625c
if skip_as_script(@__MODULE__)
	using PlutoUI
end

# ╔═╡ ba94c9b3-55bf-486a-8bf3-eb0b3e65e536
"""
	@skip_as_script expression

Marks a expression as Pluto-only, which means that it won't be executed when running outside Pluto. Do not use this for your own projects.
"""
macro skip_as_script(ex) skip_as_script(__module__) ? esc(ex) : nothing end

# ╔═╡ 9e56dbb9-2e93-4136-b10b-7b1ecb07e6f1
macro only_as_script(ex) skip_as_script(__module__) ? nothing : esc(ex) end

# ╔═╡ 5028fc01-4e14-4550-8ae6-48703b86dc21
import AbstractPlutoDingetjes.Bonds

# ╔═╡ e9d5e84d-b6b1-434d-9580-c08b0c45b9dd
import AbstractPlutoDingetjes

# ╔═╡ ccbe4423-801b-4c20-ad2d-ee89d5cfa859
md"""
# Holding intermediate updates
"""

# ╔═╡ 85918609-5e1f-4040-99be-61c2dd8ff654
md"""
## The magic
"""

# ╔═╡ 957858f8-4ace-4cae-bb16-49569daa9869
const set_input_value_compat = HypertextLiteral.JavaScript("""
(() => {
	let result = null
	try {
	result = setBoundElementValueLikePluto
} catch (e) {
	result = ((input, new_value) => {
	// fallback in case https://github.com/fonsp/Pluto.jl/pull/1755 is not available
    if (new_value == null) {
        //@ts-ignore
        input.value = new_value
        return
    }
    if (input instanceof HTMLInputElement) {
        switch (input.type) {
            case "range":
            case "number": {
                if (input.valueAsNumber !== new_value) {
                    input.valueAsNumber = new_value
                }
                return
            }
            case "date": {
                if (input.valueAsDate == null || Number(input.valueAsDate) !== Number(new_value)) {
                    input.valueAsDate = new_value
                }
                return
            }
            case "checkbox": {
                if (input.checked !== new_value) {
                    input.checked = new_value
                }
                return
            }
            case "file": {
                // Can't set files :(
                return
            }
        }
    } else if (input instanceof HTMLSelectElement && input.multiple) {
        for (let option of Array.from(input.options)) {
            option.selected = new_value.includes(option.value)
        }
        return
    }
    //@ts-ignore
    if (input.value !== new_value) {
        //@ts-ignore
        input.value = new_value
    }
})
}
return result
})()
""")

# ╔═╡ 19613f3f-5825-45a4-8951-8ff1043d0867
begin
	Base.@kwdef struct ConfirmBond
		element::Any
		secret_key::String=String(rand('a':'z', 10))
	end

	function Base.show(io::IO, m::MIME"text/html", cb::ConfirmBond)
		if !AbstractPlutoDingetjes.is_supported_by_display(io, Bonds.transform_value)
			return Base.show(io, m, HTML("❌ You need to update Pluto to use this PlutoUI element."))
		end
		output = @htl(
			"""$(
				cb.element
			)<button id=$(cb.secret_key)>Confirm</button
			>""")
		Base.show(io, m, output)
	end

	function Bonds.initial_value(cb::ConfirmBond)
		Bonds.initial_value(cb.element)
	end
	function Bonds.validate_value(cb::ConfirmBond, from_js)
		Bonds.validate_value(cb.element, from_js)
	end
	function Bonds.transform_value(cb::ConfirmBond, from_js)
		Bonds.transform_value(cb.element, from_js)
	end
	function Bonds.possible_values(cb::ConfirmBond)
		Bonds.possible_values(cb.element)
	end

end

# ╔═╡ a20da18f-7a74-43ca-9b66-1f3b82efa0c3
"""
```julia
confirm(element::Any)
```

Normally, when you move a [`Slider`](@ref) or type in a [`TextField`](@ref), all intermediate values are sent back to `@bind`.

By wrapping an input element in `confirm`, you get a button to manually control when the value is sent, intermediate updates are hidden from Pluto.

One case where this is useful is a notebook that does a long computation based on a `@bind` input.

# Examples

```julia
@bind x confirm(Slider(1:100))

x == 1 # (initially)
```

Now, when you move the slider, nothing happens, until you click the `"Confirm"` button, and `x` is set to the new slider value.

> The result looks like:
>
> ![screenshot of running the code above in pluto](https://user-images.githubusercontent.com/6933510/145615211-2731f1f5-5c7d-4aac-9aa5-8afe89e24a6b.png)

# See also

You can combine this with [`PlutoUI.combine`](@ref)!

```julia
@bind speeds confirm(
	combine() do Child
		@htl(""\"
		<h3>Wind speeds</h3>
		<ul>
		\$([
			@htl("<li>\$(name): \$(Child(name, Slider(1:100)))")
			for name in ["North", "East", "South", "West"]
		])
		</ul>
		""\")
	end
)
```

> The result looks like:
>
> ![screenshot of running the code above in pluto](https://user-images.githubusercontent.com/6933510/145614965-7a1e8630-4766-4589-8a84-b022bdfb09fc.gif)

"""
function confirm(element::Any)
	ConfirmBond(;
		element=element,
	)
end

# ╔═╡ ad5cffa5-313c-4de9-9360-005365b40780
export confirm

# ╔═╡ 6c8a03e4-7d8e-4aa4-a750-7b815622147d
md"""
## Examples
"""

# ╔═╡ f08680b2-ed28-4fac-838c-2eca7e75c6dc
@skip_as_script xb = @bind x confirm(Slider([sin, cos, tan, sqrt, floor]))

# ╔═╡ 41d64a95-6c03-4a5d-9fd7-496d0e6c0346
@skip_as_script xb

# ╔═╡ 8dcb2498-00cf-49a8-8074-301fe88b76ea
@skip_as_script x(234)

# ╔═╡ fc92b3fc-6143-477c-a413-84dcd1b4cfc0
@bind t confirm(html"<input>")

# ╔═╡ 363a65ec-218c-43a2-b740-8061fac25011
t

# ╔═╡ 801fb021-73a0-4114-a36a-328e84f00b51
@skip_as_script @bind speeds identity(
	combine() do Child
		@htl("""

Wind speeds

		$([
			@htl("
$(name): $(Child(name, Slider(1:100)))")
			for name in ["North", "East", "South", "West"]
		])

		""")
	end
)

# ╔═╡ d7985844-5944-42b9-ad41-599cd72eea82
@skip_as_script speeds

# ╔═╡ 5a0196d0-e19f-4202-b36d-18ab9be839b3
md"""
## Tests
"""

# ╔═╡ 3d3a6abb-bea7-41e2-862d-9536a9687ea5
md"""
### Initial value & transform
"""

# ╔═╡ 451c3c11-c3dc-4d6e-b9e5-c3904e405812
@skip_as_script begin
	itvs = []

	itb = @bind itv confirm(Slider([sin, cos]))
end

# ╔═╡ bb9eddc3-208c-4113-8988-311eb1dcfcd4
@skip_as_script push!(itvs, itv)

# ╔═╡ 7de50435-e962-4405-b06d-83fbe9436cff
md"""
### Combine inside combine
"""

# ╔═╡ 829f90ac-c4d9-4af9-984a-b4e52a075460
@skip_as_script cb2b = @bind wowz confirm(confirm(Slider(1:1//3:10)))

# ╔═╡ a1a9b22f-9df7-4c5f-a73d-38e184da7c35
@skip_as_script let
	sleep(1)
	wowz
end

# ╔═╡ c7a32ab2-2a5b-485b-9795-4a743604de82
@skip_as_script cb2b

# ╔═╡ Cell order:
# ╟─881b75d3-cebe-4a53-8cf5-beaeeddafd35
# ╟─ba94c9b3-55bf-486a-8bf3-eb0b3e65e536
# ╟─9e56dbb9-2e93-4136-b10b-7b1ecb07e6f1
# ╟─dadf2f40-1764-47a4-b560-683b6479d77f
# ╠═1b737805-a411-4585-b215-d0f99eafac0c
# ╠═5028fc01-4e14-4550-8ae6-48703b86dc21
# ╠═e9d5e84d-b6b1-434d-9580-c08b0c45b9dd
# ╟─ccbe4423-801b-4c20-ad2d-ee89d5cfa859
# ╠═ad5cffa5-313c-4de9-9360-005365b40780
# ╠═a20da18f-7a74-43ca-9b66-1f3b82efa0c3
# ╟─85918609-5e1f-4040-99be-61c2dd8ff654
# ╠═19613f3f-5825-45a4-8951-8ff1043d0867
# ╟─957858f8-4ace-4cae-bb16-49569daa9869
# ╟─6c8a03e4-7d8e-4aa4-a750-7b815622147d
# ╠═38b7eeb9-80bb-4a3a-a2d2-809fc423625c
# ╠═f08680b2-ed28-4fac-838c-2eca7e75c6dc
# ╠═41d64a95-6c03-4a5d-9fd7-496d0e6c0346
# ╠═8dcb2498-00cf-49a8-8074-301fe88b76ea
# ╠═fc92b3fc-6143-477c-a413-84dcd1b4cfc0
# ╠═363a65ec-218c-43a2-b740-8061fac25011
# ╠═801fb021-73a0-4114-a36a-328e84f00b51
# ╠═d7985844-5944-42b9-ad41-599cd72eea82
# ╟─5a0196d0-e19f-4202-b36d-18ab9be839b3
# ╟─3d3a6abb-bea7-41e2-862d-9536a9687ea5
# ╠═451c3c11-c3dc-4d6e-b9e5-c3904e405812
# ╠═bb9eddc3-208c-4113-8988-311eb1dcfcd4
# ╟─7de50435-e962-4405-b06d-83fbe9436cff
# ╠═829f90ac-c4d9-4af9-984a-b4e52a075460
# ╠═a1a9b22f-9df7-4c5f-a73d-38e184da7c35
# ╠═c7a32ab2-2a5b-485b-9795-4a743604de82


text:
version = v"12.0.0"
rocm_version = v"4.2.0"

include("../../L/LLVM/common.jl")

build_tarballs(ARGS, configure_build(ARGS, version; experimental_platforms=true,
                                     git_path="https://github.com/RadeonOpenCompute/llvm-project",
                                     git_ver="b204d7f0cae65b6cd4446eec50fc1fb675d582af",
                                     custom_name="ROCmLLVM",
                                     custom_version=rocm_version, static=true,
                                     platform_filter=Sys.islinux)...;
               preferred_gcc_version=v"7", preferred_llvm_version=v"9", julia_compat="1.7")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using LightXML
using JLD
using HDF5
"""
Usage:
    julia ims_cp.jl imaris_root_dir exp_folder_name imaris_creation_name destination

Version  Commit
0.1      first version
TODO: combine scene and image file
"""

"copy and rename imaris batch result from imaris database with Folder and Creation name"
function ims_cp(ims_root::String, ims_array::String, ims_creation::String, output::String)
    id2name = Dict()
    id2path= Dict()
    # load ims TreeData.xml
    TREEDATE = parse_file("$ims_root/TreeData.xml")
    # enter root node
    xroot = root(TREEDATE)

    # local -> mRNA -> ims_array
    for folder in child_elements(xroot)
        if (name(folder) =="Folder") && (attribute(folder, "name") == "local")
            for folder2 in child_elements(folder)
                if (name(folder2) == "Folder") && (attribute(folder2, "name") == "mRNA")
                    for folder3 in child_elements(folder2)
                        # enter single experiment folder
                        #println(name(folder3)," ", attribute(folder3, "name"))
                        if (name(folder3) == "Folder") && (attribute(folder3, "name") == ims_array)
                            # enter "Image" node for each folder
                            for e in child_elements( folder3 )
                                # travel all image
                                if name(e) == "Image"
                                    image_name = attribute(e, "name")
                                    image_id = attribute(e, "id")
                                    image_path = replace( attribute(e, "path"),
                                                         "file:///C%3A/Imaris_Data_qblab/"=>"" )
                                    push!(id2name, image_id=>image_name)
                                    push!(id2path, image_id=>image_path)
                                elseif name(e) == "CreationParameter"
                                    # get CreationParameter id
                                    creation_parameters = attribute(e, "CreationParameters")
                                    open("$output/$ims_creation.txt", "w") do io
                                        for line in split(creation_parameters, "> <")
                                            println(io, line)
                                        end
                                    end
                                end
                            end
                            # cp batch result
                            for e in child_elements( folder3 )
                                if (name(e) == "BatchRun") && (attribute(e, "name") == "$ims_array $ims_creation")
                                    #println(e["batchrunresults"]["result"])
                                    batchrunresults = e["batchrunresults"][1]["result"]
                                    for result in batchrunresults
                                        if attribute(result, "status") == "FINISHED" && !isempty(attribute(result, "output"))
                                            id = attribute(result, "input")
                                            outputpath = attribute(result, "outputpath")
                                            src_scene = replace(outputpath, "file:/"=>"$ims_root")
                                            src_ims = "$ims_root/$(id2path[id])"
                                            out_ims = "$output/$(id2name[id]).ims";
                                            print(id2name[id], " ")
                                            cp(src_ims, out_ims, force=true)
                                            ims_f = h5open(out_ims, "r+")
                                            scene_f = h5open(src_scene, "r")
                                            if haskey(ims_f, "Scene")
                                                HDF5.delete_object(ims_f, "Scene")
                                            end
                                            if haskey(ims_f, "Scene8")
                                                HDF5.delete_object(ims_f, "Scene8")
                                            end
                                            HDF5.copy_object(scene_f, "Scene", ims_f, "Scene")
                                            HDF5.copy_object(scene_f, "Scene8", ims_f, "Scene8")
                                            close(ims_f)
                                            close(scene_f)
                                        end
                                    end
                                else
                                    #println("$ims_creation no found!")
                                end
                            end
                        else
                            #println("$ims_array no found!")
                        end
                    end
                end
            end
        end
    end
end

if length(ARGS) != 4
    println("Usage: julia ims_cp.jl imaris_root_dir exp_folder_name imaris_creation_name destination")
else
    ims_cp(ARGS[1], ARGS[2], ARGS[3], ARGS[4])
end


text:
α_data = [-10.0 -5.0 0.0 5.0 10.0 15.0 20.0 25.0 30.0 35.0 40.0 45.0]  # deg
β_data = [0.0 5.0 10.0 15.0 20.0 25.0 30.0]  # deg

de_data = [-24.0 -12.0 0.0 12.0 24.0]  # deg

# UTILS for interpolations
function get_interp_idx(val, coeff, min_, max_, off_ = 3)
    s = val * coeff
    k = floor(Int, s)
    k = max(min_, k)
    k = min(k, max_)
    da = s - float(k)
    k += off_
    l = k + Int(sign(da))
    return k, l, da
end

function get_interp_idx2(val, coeff, min_, max_)
    s = abs(val) * coeff
    k = floor(Int, s)
    k = max(min_, k)
    k = min(k, max_)
    da = s - float(k)
    k += 1
    l = k + Int(sign(da))
    return k, l, da
end

function interp1d(val, coeff, min, max, data)
    k, l, da = get_interp_idx(val, coeff, min, max)
    res = data[k] + abs(da) * (data[l] - data[k])
    return res
end

function interp2d(
    val1,
    val2,
    coeff1,
    coeff2,
    min1,
    min2,
    max1,
    max2,
    data,
    off1 = 3,
    off2 = 3,
)
    k, l, da = get_interp_idx(val1, coeff1, min1, max1, off1)
    m, n, de = get_interp_idx(val2, coeff2, min2, max2, off2)

    t = data[k, m]
    u = data[k, n]
    v = t + abs(da) * (data[l, m] - t)
    w = u + abs(da) * (data[l, n] - u)
    res = v + (w - v) * abs(de)

    return res
end

function interp2d2(val1, val2, coeff1, coeff2, min1, min2, max1, max2, data)
    k, l, da = get_interp_idx(val1, coeff1, min1, max1)
    m, n, de = get_interp_idx2(val2, coeff2, min2, max2)

    t = data[k, m]
    u = data[k, n]
    v = t + abs(da) * (data[l, m] - t)
    w = u + abs(da) * (data[l, n] - u)
    res = v + (w - v) * abs(de)
    res *= sign(val2)

    return res
end

# DAMPING COEFFICIENTS
CXq_data = [-0.267 -0.110 0.308 1.340 2.080 2.910 2.760 2.050 1.500 1.490 1.830 1.210]
CYr_data = [0.882 0.852 0.876 0.958 0.962 0.974 0.819 0.483 0.590 1.210 -0.493 -1.040]
CYp_data = [-0.108 -0.108 -0.188 0.110 0.258 0.226 0.344 0.362 0.611 0.529 0.298 -0.227]  # <<-- Morelli
# CYp_data = [-0.108 -0.108 -0.188 0.110 0.258 0.226 0.344 0.362 0.611 0.529 0.298 -2.27]  # <<-- Stevens
CZq_data =
    [-8.800 -25.800 -28.900 -31.400 -31.200 -30.700 -27.700 -28.200 -29.000 -29.800 -38.300 -35.300]
Clr_data = [-0.126 -0.026 0.063 0.113 0.208 0.230 0.319 0.437 0.680 0.100 0.447 -0.330]
Clp_data =
    [-0.360 -0.359 -0.443 -0.420 -0.383 -0.375 -0.329 -0.294 -0.230 -0.210 -0.120 -0.100]
Cmq_data =
    [-7.210 -5.40 -5.230 -5.260 -6.110 -6.640 -5.690 -6.000 -6.200 -6.400 -6.600 -6.000]  # <<-- Morelli
# Cmq_data =
#     [-7.210 -0.540 -5.230 -5.260 -6.110 -6.640 -5.690 -6.000 -6.200 -6.400 -6.600 -6.000]  # <<-- Stevens
Cnr_data =
    [-0.380 -0.363 -0.378 -0.386 -0.370 -0.453 -0.550 -0.582 -0.595 -0.637 -1.020 -0.840]
Cnp_data = [0.061 0.052 0.052 -0.012 -0.013 -0.024 0.050 0.150 0.130 0.158 0.240 0.150]

damp_data =
    [
        CXq_data
        CYr_data
        CYp_data
        CZq_data
        Clr_data
        Clp_data
        Cmq_data
        Cnr_data
        Cnp_data
    ]'
damp_data = damp_data'

# CX
# rows: de, columns: alpha
CX_data = [
    -0.099 -0.081 -0.081 -0.063 -0.025 0.044 0.097 0.113 0.145 0.167 0.174 0.166
    -0.048 -0.038 -0.040 -0.021 0.016 0.083 0.127 0.137 0.162 0.177 0.179 0.167
    -0.022 -0.020 -0.021 -0.004 0.032 0.094 0.128 0.130 0.154 0.161 0.155 0.138
    -0.040 -0.038 -0.039 -0.025 0.006 0.062 0.087 0.085 0.100 0.110 0.104 0.091
    -0.083 -0.073 -0.076 -0.072 -0.046 0.012 0.024 0.025 0.043 0.053 0.047 0.040
]
CX_data = CX_data'

# CZ
# columns: alpha
CZ_data =
    [0.770 0.241 -0.100 -0.416 -0.731 -1.053 -1.366 -1.646 -1.917 -2.120 -2.248 -2.229]

# Cm
# rows: de columns: alpha
Cm_data = [
    0.205 0.168 0.186 0.196 0.213 0.251 0.245 0.238 0.252 0.231 0.198 0.192
    0.081 0.077 0.107 0.110 0.110 0.141 0.127 0.119 0.133 0.108 0.081 0.093
    -0.046 -0.020 -0.009 -0.005 -0.006 0.010 0.006 -0.001 0.014 0.000 -0.013 0.032
    -0.174 -0.145 -0.121 -0.127 -0.129 -0.102 -0.097 -0.113 -0.087 -0.084 -0.069 -0.006
    -0.259 -0.202 -0.184 -0.193 -0.199 -0.150 -0.160 -0.167 -0.104 -0.076 -0.041 -0.005
]
Cm_data = Cm_data'

# Cl
# Stevens
Cl_data = [
    0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
    -0.001 -0.004 -0.008 -0.012 -0.016 -0.019 -0.020 -0.020 -0.015 -0.008 -0.013 -0.015
    -0.003 -0.009 -0.017 -0.024 -0.030 -0.034 -0.040 -0.037 -0.016 -0.002 -0.010 -0.019
    -0.001 -0.010 -0.020 -0.030 -0.039 -0.044 -0.050 -0.049 -0.023 -0.006 -0.014 -0.027
    0.000 -0.010 -0.022 -0.034 -0.047 -0.046 -0.059 -0.061 -0.033 -0.036 -0.035 -0.035
    0.007 -0.010 -0.023 -0.034 -0.049 -0.046 -0.068 -0.071 -0.060 -0.058 -0.062 -0.059
    0.009 -0.011 -0.023 -0.037 -0.050 -0.047 -0.074 -0.079 -0.091 -0.076 -0.077 -0.076
]

# Morelli version
# found differences in Morelli Cl_data for columns 6-9
# Cl_data = [ 0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000  0.000;
#            -0.001 -0.004 -0.008 -0.012 -0.016 -0.022 -0.022 -0.021 -0.015 -0.008 -0.013 -0.015;
#            -0.003 -0.009 -0.017 -0.024 -0.030 -0.041 -0.045 -0.040 -0.016 -0.002 -0.010 -0.019;
#            -0.001 -0.010 -0.020 -0.030 -0.039 -0.054 -0.057 -0.054 -0.023 -0.006 -0.014 -0.027;
#             0.000 -0.010 -0.022 -0.034 -0.047 -0.060 -0.069 -0.067 -0.033 -0.036 -0.035 -0.035;
#             0.007 -0.010 -0.023 -0.034 -0.049 -0.063 -0.081 -0.079 -0.060 -0.058 -0.062 -0.059;
#             0.009 -0.011 -0.023 -0.037 -0.050 -0.068 -0.089 -0.088 -0.091 -0.076 -0.077 -0.076]

Cl_data = Cl_data'

# rows: beta columns=alpha
clda_data = [
    -0.041 -0.052 -0.053 -0.056 -0.050 -0.056 -0.082 -0.059 -0.042 -0.038 -0.027 -0.017
    -0.041 -0.053 -0.053 -0.053 -0.050 -0.051 -0.066 -0.043 -0.038 -0.027 -0.023 -0.016
    -0.042 -0.053 -0.052 -0.051 -0.049 -0.049 -0.043 -0.035 -0.026 -0.016 -0.018 -0.014
    -0.040 -0.052 -0.051 -0.052 -0.048 -0.048 -0.042 -0.037 -0.031 -0.026 -0.017 -0.012
    -0.043 -0.049 -0.048 -0.049 -0.043 -0.042 -0.042 -0.036 -0.025 -0.021 -0.016 -0.011
    -0.044 -0.048 -0.048 -0.047 -0.042 -0.041 -0.020 -0.028 -0.013 -0.014 -0.011 -0.010
    -0.043 -0.049 -0.047 -0.045 -0.042 -0.037 -0.003 -0.013 -0.010 -0.003 -0.007 -0.008
]
clda_data = clda_data'

# rows: beta columns=alpha
cldr_data = [
    0.005 0.017 0.014 0.010 -0.005 0.009 0.019 0.005 -0.000 -0.005 -0.011 0.008
    0.007 0.016 0.014 0.014 0.013 0.009 0.012 0.005 0.000 0.004 0.009 0.007
    0.013 0.013 0.011 0.012 0.011 0.009 0.008 0.005 -0.002 0.005 0.003 0.005  # <<-- Stevens
  # 0.013 0.013 0.011 0.012 0.011 0.009 0.008 0.005  0.000 0.005 0.003 0.005  # <<-- Morelli
    0.018 0.015 0.015 0.014 0.014 0.014 0.014 0.015 0.013 0.011 0.006 0.001
    0.015 0.014 0.013 0.013 0.012 0.011 0.011 0.010 0.008 0.008 0.007 0.003
    0.021 0.011 0.010 0.011 0.010 0.009 0.008 0.010 0.006 0.005 0.000 0.001
    0.023 0.010 0.011 0.011 0.011 0.010 0.008 0.010 0.006 0.014 0.020 0.000
]
cldr_data = cldr_data'

# Cn
# rows: beta columns: alpha
Cn_data = [
    0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
    0.018 0.019 0.018 0.019 0.019 0.018 0.013 0.007 0.004 -0.014 -0.017 -0.033
    0.038 0.042 0.042 0.042 0.043 0.039 0.030 0.017 0.004 -0.035 -0.047 -0.057
    0.056 0.057 0.059 0.058 0.058 0.053 0.032 0.012 0.002 -0.046 -0.071 -0.073
    0.064 0.077 0.076 0.074 0.073 0.057 0.029 0.007 0.012 -0.034 -0.065 -0.041
    0.074 0.086 0.093 0.089 0.080 0.062 0.049 0.022 0.028 -0.012 -0.002 -0.013
    0.079 0.090 0.106 0.106 0.096 0.080 0.068 0.030 0.064 0.015 0.011 -0.001
]
Cn_data = Cn_data'

# rows: beta columns=alpha
cnda_data = [
    0.001 -0.027 -0.017 -0.013 -0.012 -0.016 0.001 0.017 0.011 0.017 0.008 0.016
    0.002 -0.014 -0.016 -0.016 -0.014 -0.019 -0.021 0.002 0.012 0.015 0.015 0.011  # <<-- Stevens
  # 0.002 -0.014 -0.016 -0.016 -0.014 -0.019 -0.021 0.002 0.012 0.015 0.016 0.011  # <<-- Morelli
    -0.006 -0.008 -0.006 -0.006 -0.005 -0.008 -0.005 0.007 0.004 0.007 0.006 0.006
    -0.011 -0.011 -0.010 -0.009 -0.008 -0.006 0.000 0.004 0.007 0.010 0.004 0.010
    -0.015 -0.015 -0.014 -0.012 -0.011 -0.008 -0.002 0.002 0.006 0.012 0.011 0.011
    -0.024 -0.010 -0.004 -0.002 -0.001 0.003 0.014 0.006 -0.001 0.004 0.004 0.006
    -0.022 0.002 -0.003 -0.005 -0.003 -0.001 -0.009 -0.009 -0.001 0.003 -0.002 0.001
]
cnda_data = cnda_data'

# rows: beta columns=alpha
cndr_data = [
    -0.018 -0.052 -0.052 -0.052 -0.054 -0.049 -0.059 -0.051 -0.030 -0.037 -0.026 -0.013
    -0.028 -0.051 -0.043 -0.046 -0.045 -0.049 -0.057 -0.052 -0.030 -0.033 -0.030 -0.008
    -0.037 -0.041 -0.038 -0.040 -0.040 -0.038 -0.037 -0.030 -0.027 -0.024 -0.019 -0.013
    -0.048 -0.045 -0.045 -0.045 -0.044 -0.045 -0.047 -0.048 -0.049 -0.045 -0.033 -0.016
    -0.043 -0.044 -0.041 -0.041 -0.040 -0.038 -0.034 -0.035 -0.035 -0.029 -0.022 -0.009
    -0.052 -0.034 -0.036 -0.036 -0.035 -0.028 -0.024 -0.023 -0.020 -0.016 -0.010 -0.014
    -0.062 -0.034 -0.027 -0.028 -0.027 -0.027 -0.023 -0.023 -0.019 -0.009 -0.025 -0.010
]
cndr_data = cndr_data'

"""
   dampings(α)

Calculates CX_q, CY_r, CY_p, CZ_q, Cl_r, Cl_p, Cm_q, Cn_r, Cn_p for the given
angle of attack, α (deg).
"""
function damp(α)
    dampings = Array{Float64}(undef, 9)
    for ii = 1:9
        dampings[ii] = interp1d(α, 0.2, -1, 8, damp_data[ii, :])
    end
    return dampings
end

CYdr(dr) = 0.086 / DR_MAX * dr
CYda(da) = 0.021 / DA_MAX * da
CYβ(β) = -0.02 * β

CZα(α) = interp1d(α, 0.2, -1, 8, CZ_data)
CZαβ(α, β) = CZα(α) * (1.0 - (β / 57.3)^2)
CZde(de) = -0.19 / DE_MAX * de

CX(α, de) = interp2d(α, de, 0.2, 1.0 / 12, -1, -1, 8, 1, CX_data)
CY(β, da, dr) = CYβ(β) + CYda(da) + CYdr(dr)
CZ(α, β, de) = CZαβ(α, β) + CZde(de)

CL(α, β) = interp2d2(α, β, 0.2, 0.2, -1, 1, 8, 5, Cl_data)
DLDA(α, β) = interp2d(α, β, 0.2, 0.1, -1, -2, 8, 2, clda_data, 3, 4)
DLDR(α, β) = interp2d(α, β, 0.2, 0.1, -1, -2, 8, 2, cldr_data, 3, 4)

CM(α, de) = interp2d(α, de, 0.2, 1.0 / 12, -1, -1, 8, 1, Cm_data)

CN(α, β) = interp2d2(α, β, 0.2, 0.2, -1, 1, 8, 5, Cn_data)
DNDA(α, β) = interp2d(α, β, 0.2, 0.1, -1, -2, 8, 2, cnda_data, 3, 4)
DNDR(α, β) = interp2d(α, β, 0.2, 0.1, -1, -2, 8, 2, cndr_data, 3, 4)

function calculate_aero_forces_moments(ac::F16, dss::DSState, controls, xcg, qbar, S, b, c)

    # Unpack controls
    THTL = controls[1]  # 0-1
    EL = controls[2]  # deg
    AIL = controls[3]  # deg
    RDR = controls[4]  # deg

    # Assign state
    VT = get_tas(dss)
    ALPHA = get_α(dss) * RAD2DEG
    BETA = get_β(dss) * RAD2DEG
    PSI, THETA, PHI = get_euler_angles(dss)
    P, Q, R = get_ang_vel_body(dss)

    # Look-up tables and component buildup
    CXT = CX(ALPHA, EL)
    CYT = CY(BETA, AIL, RDR)
    CZT = CZ(ALPHA, BETA, EL)

    DAIL = AIL / DA_MAX
    DRDR = RDR / DR_MAX

    CLT = CL(ALPHA, BETA) + DLDA(ALPHA, BETA) * DAIL + DLDR(ALPHA, BETA) * DRDR
    CMT = CM(ALPHA, EL)
    CNT = CN(ALPHA, BETA) + DNDA(ALPHA, BETA) * DAIL + DNDR(ALPHA, BETA) * DRDR

    # Add damping derivatives
    CBTA = cos(get_β(dss))
    U = VT * cos(get_α(dss)) * CBTA
    V = VT * sin(get_β(dss))
    W = VT * sin(get_α(dss)) * CBTA
    TVT = 0.5 / VT
    B2V = B * TVT
    CQ = CBAR * Q * TVT

    D = damp(ALPHA)
    # Non-dimensional coefficients
    CXT = CXT + CQ * D[1]
    CYT = CYT + B2V * (D[2] * R + D[3] * P)
    CZT = CZT + CQ * D[4]
    CLT = CLT + B2V * (D[5] * R + D[6] * P)
    CMT = CMT + CQ * D[7] + CZT * (XCGR - xcg)
    CNT = CNT + B2V * (D[8] * R + D[9] * P) - CYT * (XCGR - xcg) * CBAR / B

    # Correct differences in gravity acceleration between Stevens and constants.jl
    CXT, CYT, CZT = [CXT, CYT, CZT] .* (GD * FT2M) / gD

    # Dimensional forces and moments
    qbarS = qbar * S

    Fx = qbarS * CXT
    Fy = qbarS * CYT
    Fz = qbarS * CZT
    L = qbarS * b * CLT
    M = qbarS * c * CMT
    N = qbarS * b * CNT

    return [Fx, Fy, Fz, L, M, N]
end


text:
const HeaderBinary = Tuple{Vector{UInt8}, Vector{UInt8}}
const Header = Tuple{String, String}

mutable struct DynamicTable
    table::Vector{HeaderBinary}
    size::Int
    max_size::Int

    # 4096 is the initial value of SETTINGS_HEADER_TABLE_SIZE (rfc7540)
    DynamicTable(max_size::Int=4096) = new(HeaderBinary[], 0, max_size)
end

function consolidate_table!(table::DynamicTable)
    while table.size > table.max_size
        last_header = pop!(table.table)
        table.size -= length(last_header[1]) + length(last_header[2]) + 32
    end
end

function add_header!(table::DynamicTable, name::Vector{UInt8}, value::Vector{UInt8})
    table.size += length(name) + length(value) + 32
    pushfirst!(table.table, (name, value))
    consolidate_table!(table)
end
add_header!(table::DynamicTable, header::HeaderBinary) =
    add_header!(table, header[1], header[2])
add_header!(table::DynamicTable, header::Header) =
    add_header!(table, bytearr(header[1]), bytearr(header[2]))

function set_max_table_size!(table::DynamicTable, size::Int)
    table.max_size = size
    consolidate_table!(table)
end

const STATIC_TABLE = [
     (bytearr(":authority"                    ), bytearr(""             ))
     (bytearr(":method"                       ), bytearr("GET"          ))
     (bytearr(":method"                       ), bytearr("POST"         ))
     (bytearr(":path"                         ), bytearr("/"            ))
     (bytearr(":path"                         ), bytearr("/index.html"  ))
     (bytearr(":scheme"                       ), bytearr("http"         ))
     (bytearr(":scheme"                       ), bytearr("https"        ))
     (bytearr(":status"                       ), bytearr("200"          ))
     (bytearr(":status"                       ), bytearr("204"          ))
     (bytearr(":status"                       ), bytearr("206"          ))
     (bytearr(":status"                       ), bytearr("304"          ))
     (bytearr(":status"                       ), bytearr("400"          ))
     (bytearr(":status"                       ), bytearr("404"          ))
     (bytearr(":status"                       ), bytearr("500"          ))
     (bytearr("accept-charset"                ), bytearr(""             ))
     (bytearr("accept-encoding"               ), bytearr("gzip, deflate"))
     (bytearr("accept-language"               ), bytearr(""             ))
     (bytearr("accept-ranges"                 ), bytearr(""             ))
     (bytearr("accept"                        ), bytearr(""             ))
     (bytearr("access-control-allow-origin"   ), bytearr(""             ))
     (bytearr("age"                           ), bytearr(""             ))
     (bytearr("allow"                         ), bytearr(""             ))
     (bytearr("authorization"                 ), bytearr(""             ))
     (bytearr("cache-control"                 ), bytearr(""             ))
     (bytearr("content-disposition"           ), bytearr(""             ))
     (bytearr("content-encoding"              ), bytearr(""             ))
     (bytearr("content-language"              ), bytearr(""             ))
     (bytearr("content-length"                ), bytearr(""             ))
     (bytearr("content-location"              ), bytearr(""             ))
     (bytearr("content-range"                 ), bytearr(""             ))
     (bytearr("content-type"                  ), bytearr(""             ))
     (bytearr("cookie"                        ), bytearr(""             ))
     (bytearr("date"                          ), bytearr(""             ))
     (bytearr("etag"                          ), bytearr(""             ))
     (bytearr("expect"                        ), bytearr(""             ))
     (bytearr("expires"                       ), bytearr(""             ))
     (bytearr("from"                          ), bytearr(""             ))
     (bytearr("host"                          ), bytearr(""             ))
     (bytearr("if-match"                      ), bytearr(""             ))
     (bytearr("if-modified-since"             ), bytearr(""             ))
     (bytearr("if-none-match"                 ), bytearr(""             ))
     (bytearr("if-range"                      ), bytearr(""             ))
     (bytearr("if-unmodified-since"           ), bytearr(""             ))
     (bytearr("last-modified"                 ), bytearr(""             ))
     (bytearr("link"                          ), bytearr(""             ))
     (bytearr("location"                      ), bytearr(""             ))
     (bytearr("max-forwards"                  ), bytearr(""             ))
     (bytearr("proxy-authenticate"            ), bytearr(""             ))
     (bytearr("proxy-authorization"           ), bytearr(""             ))
     (bytearr("range"                         ), bytearr(""             ))
     (bytearr("referer"                       ), bytearr(""             ))
     (bytearr("refresh"                       ), bytearr(""             ))
     (bytearr("retry-after"                   ), bytearr(""             ))
     (bytearr("server"                        ), bytearr(""             ))
     (bytearr("set-cookie"                    ), bytearr(""             ))
     (bytearr("strict-transport-security"     ), bytearr(""             ))
     (bytearr("transfer-encoding"             ), bytearr(""             ))
     (bytearr("user-agent"                    ), bytearr(""             ))
     (bytearr("vary"                          ), bytearr(""             ))
     (bytearr("via"                           ), bytearr(""             ))
     (bytearr("www-authenticate"              ), bytearr(""             ))
 ]

@noinline bounderr() =
    throw(DecodeError("Index out of bound."))
@noinline indexerr() =
    throw(DecodeError("Index greater than sum of both static and dynamic tables."))

function get_header(table::DynamicTable, index)
    # IETF's table indexing is 1-based.
    index <= length(STATIC_TABLE) && return STATIC_TABLE[index]
    (index > length(STATIC_TABLE) + table.max_size) && indexerr()
    (index > length(STATIC_TABLE) + length(table.table)) && bounderr()
    return table.table[index - length(STATIC_TABLE)]
end


text:
using spc
using Test

@testset "spc.jl" begin
    # Write your tests here.
end

max_stars_repo_name: DaanVanHauwermeiren/spc

max_issues_repo_name: DaanVanHauwermeiren/spc

max_forks_repo_name: DaanVanHauwermeiren/spc


text:
module FastGMM

include("GMMTypes.jl")
include("FastGMMCuda.jl")

end


text:
module Libtask

using Libdl
using Libtask_jll

export CTask, consume, produce, TArray, tzeros, tfill, TRef

function __init__()
    @static if VERSION < v"1.6.0"
        push!(Libdl.DL_LOAD_PATH,
              Libtask_jll.get_libtask_julia_path() |> dirname)
        Libdl.dlopen(Libtask_jll.libtask_julia_path)
    end
end

include("memlayout/main.jl")
include("ctask.jl")
include("tarray.jl")
include("tref.jl")

end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test, Mmap

file = tempname()
write(file, "Hello World\n")
t = b"Hello World"
@test Mmap.mmap(file, Array{UInt8,3}, (11,1,1)) == reshape(t,(11,1,1))
gc(); gc()
@test Mmap.mmap(file, Array{UInt8,3}, (1,11,1)) == reshape(t,(1,11,1))
gc(); gc()
@test Mmap.mmap(file, Array{UInt8,3}, (1,1,11)) == reshape(t,(1,1,11))
gc(); gc()
@test Mmap.mmap(file, Array{UInt8,3}, (11,0,1)) == Array{UInt8}((0,0,0))
@test Mmap.mmap(file, Vector{UInt8}, (11,)) == t
gc(); gc()
@test Mmap.mmap(file, Array{UInt8,2}, (1,11)) == t'
gc(); gc()
@test Mmap.mmap(file, Array{UInt8,2}, (0,12)) == Array{UInt8}((0,0))
m = Mmap.mmap(file, Array{UInt8,3}, (1,2,1))
@test m == reshape(b"He",(1,2,1))
finalize(m); m=nothing; gc()

# constructors
@test length(@inferred Mmap.mmap(file)) == 12
@test length(@inferred Mmap.mmap(file, Vector{Int8})) == 12
@test length(@inferred Mmap.mmap(file, Matrix{Int8}, (12,1))) == 12
@test length(@inferred Mmap.mmap(file, Matrix{Int8}, (12,1), 0)) == 12
@test length(@inferred Mmap.mmap(file, Matrix{Int8}, (12,1), 0; grow=false)) == 12
@test length(@inferred Mmap.mmap(file, Matrix{Int8}, (12,1), 0; shared=false)) == 12
@test length(@inferred Mmap.mmap(file, Vector{Int8}, 12)) == 12
@test length(@inferred Mmap.mmap(file, Vector{Int8}, 12, 0)) == 12
@test length(@inferred Mmap.mmap(file, Vector{Int8}, 12, 0; grow=false)) == 12
@test length(@inferred Mmap.mmap(file, Vector{Int8}, 12, 0; shared=false)) == 12
s = open(file)
@test length(@inferred Mmap.mmap(s)) == 12
@test length(@inferred Mmap.mmap(s, Vector{Int8})) == 12
@test length(@inferred Mmap.mmap(s, Matrix{Int8}, (12,1))) == 12
@test length(@inferred Mmap.mmap(s, Matrix{Int8}, (12,1), 0)) == 12
@test length(@inferred Mmap.mmap(s, Matrix{Int8}, (12,1), 0; grow=false)) == 12
@test length(@inferred Mmap.mmap(s, Matrix{Int8}, (12,1), 0; shared=false)) == 12
@test length(@inferred Mmap.mmap(s, Vector{Int8}, 12)) == 12
@test length(@inferred Mmap.mmap(s, Vector{Int8}, 12, 0)) == 12
@test length(@inferred Mmap.mmap(s, Vector{Int8}, 12, 0; grow=false)) == 12
@test length(@inferred Mmap.mmap(s, Vector{Int8}, 12, 0; shared=false)) == 12
close(s)
@test_throws ErrorException Mmap.mmap(file, Vector{Ref}) # must be bit-type
gc(); gc()

s = open(f->f,file,"w")
@test Mmap.mmap(file) == Array{UInt8}(0) # requested len=0 on empty file
@test Mmap.mmap(file,Vector{UInt8},0) == Array{UInt8}(0)
m = Mmap.mmap(file,Vector{UInt8},12)
m[:] = b"Hello World\n"
Mmap.sync!(m)
finalize(m); m=nothing; gc()
@test open(x->read(x, String),file) == "Hello World\n"

s = open(file, "r")
close(s)
@test_throws Base.UVError Mmap.mmap(s) # closed IOStream
@test_throws ArgumentError Mmap.mmap(s,Vector{UInt8},12,0) # closed IOStream
@test_throws SystemError Mmap.mmap("")

# negative length
@test_throws ArgumentError Mmap.mmap(file, Vector{UInt8}, -1)
# negative offset
@test_throws ArgumentError Mmap.mmap(file, Vector{UInt8}, 1, -1)

for i = 0x01:0x0c
    @test length(Mmap.mmap(file, Vector{UInt8}, i)) == Int(i)
end
gc(); gc()

sz = filesize(file)
m = Mmap.mmap(file, Vector{UInt8}, sz+1)
@test length(m) == sz+1 # test growing
@test m[end] == 0x00
finalize(m); m=nothing; gc()
sz = filesize(file)
m = Mmap.mmap(file, Vector{UInt8}, 1, sz)
@test length(m) == 1
@test m[1] == 0x00
finalize(m); m=nothing; gc()
sz = filesize(file)
# test where offset is actually > than size of file; file is grown with zeroed bytes
m = Mmap.mmap(file, Vector{UInt8}, 1, sz+1)
@test length(m) == 1
@test m[1] == 0x00
finalize(m); m=nothing; gc()

s = open(file, "r")
m = Mmap.mmap(s)
@test_throws ReadOnlyMemoryError m[5] = UInt8('x') # tries to setindex! on read-only array
finalize(m); m=nothing; gc()

write(file, "Hello World\n")

s = open(file, "r")
m = Mmap.mmap(s)
close(s)
finalize(m); m=nothing; gc()
m = Mmap.mmap(file)
s = open(file, "r+")
c = Mmap.mmap(s)
d = Mmap.mmap(s)
c[1] = UInt8('J')
Mmap.sync!(c)
close(s)
@test m[1] == UInt8('J')
@test d[1] == UInt8('J')
finalize(m); finalize(c); finalize(d)
m=nothing; c=nothing; d=nothing; gc()

write(file, "Hello World\n")

s = open(file, "r")
@test isreadonly(s) == true
c = Mmap.mmap(s, Vector{UInt8}, (11,))
@test c == b"Hello World"
finalize(c); c=nothing; gc()
c = Mmap.mmap(s, Vector{UInt8}, (UInt16(11),))
@test c == b"Hello World"
finalize(c); c=nothing; gc()
@test_throws ArgumentError Mmap.mmap(s, Vector{UInt8}, (Int16(-11),))
@test_throws ArgumentError Mmap.mmap(s, Vector{UInt8}, (typemax(UInt),))
close(s)
s = open(file, "r+")
@test isreadonly(s) == false
c = Mmap.mmap(s, Vector{UInt8}, (11,))
c[5] = UInt8('x')
Mmap.sync!(c)
close(s)
s = open(file, "r")
str = readline(s)
close(s)
@test startswith(str, "Hellx World")
finalize(c); c=nothing; gc()

c = Mmap.mmap(file)
@test c == b"Hellx World\n"
finalize(c); c=nothing; gc()
c = Mmap.mmap(file, Vector{UInt8}, 3)
@test c == b"Hel"
finalize(c); c=nothing; gc()
s = open(file, "r")
c = Mmap.mmap(s, Vector{UInt8}, 6)
@test c == b"Hellx "
close(s)
finalize(c); c=nothing; gc()
c = Mmap.mmap(file, Vector{UInt8}, 5, 6)
@test c == b"World"
finalize(c); c=nothing; gc()

s = open(file, "w")
write(s, "Hello World\n")
close(s)

# test Mmap.mmap
m = Mmap.mmap(file)
tdata = b"Hello World\n"
for i = 1:12
    @test m[i] == tdata[i]
end
@test_throws BoundsError m[13]
finalize(m); m=nothing; gc()

m = Mmap.mmap(file,Vector{UInt8},6)
@test m[1] == b"H"[1]
@test m[2] == b"e"[1]
@test m[3] == b"l"[1]
@test m[4] == b"l"[1]
@test m[5] == b"o"[1]
@test m[6] == b" "[1]
@test_throws BoundsError m[7]
finalize(m); m=nothing; gc()

m = Mmap.mmap(file,Vector{UInt8},2,6)
@test m[1] == b"W"[1]
@test m[2] == b"o"[1]
@test_throws BoundsError m[3]
finalize(m); m = nothing; gc()

s = open(file, "w")
write(s, [0xffffffffffffffff,
          0xffffffffffffffff,
          0xffffffffffffffff,
          0x000000001fffffff])
close(s)
s = open(file, "r")
@test isreadonly(s)
b = @inferred Mmap.mmap(s, BitArray, (17,13))
@test Base._check_bitarray_consistency(b)
@test b == trues(17,13)
@test_throws ArgumentError Mmap.mmap(s, BitArray, (7,3))
close(s)
s = open(file, "r+")
b = Mmap.mmap(s, BitArray, (17,19))
@test Base._check_bitarray_consistency(b)
rand!(b)
Mmap.sync!(b)
b0 = copy(b)
@test Base._check_bitarray_consistency(b0)
close(s)
s = open(file, "r")
@test isreadonly(s)
b = Mmap.mmap(s, BitArray, (17,19))
@test Base._check_bitarray_consistency(b)
@test b == b0
close(s)
finalize(b); finalize(b0)
b = nothing; b0 = nothing
gc()

open(file,"w") do f
    write(f,UInt64(1))
    write(f,UInt8(1))
end
@test filesize(file) == 9
m = Mmap.mmap(file, BitArray, (72,))
@test Base._check_bitarray_consistency(m)
@test length(m) == 72
finalize(m); m = nothing; gc()
rm(file)

# Mmap.mmap with an offset
A = rand(1:20, 500, 300)
fname = tempname()
s = open(fname, "w+")
write(s, size(A,1))
write(s, size(A,2))
write(s, A)
close(s)
s = open(fname)
m = read(s, Int)
n = read(s, Int)
A2 = Mmap.mmap(s, Matrix{Int}, (m,n))
@test A == A2
seek(s, 0)
A3 = Mmap.mmap(s, Matrix{Int}, (m,n), convert(Int64, 2*sizeof(Int)))
@test A == A3
A4 = Mmap.mmap(s, Matrix{Int}, (m,150), convert(Int64, (2+150*m)*sizeof(Int)))
@test A[:, 151:end] == A4
close(s)
finalize(A2); finalize(A3); finalize(A4)
A2 = A3 = A4 = nothing
gc()
rm(fname)

# Mmap.Anonymous
m = Mmap.Anonymous()
@test m.name == ""
@test !m.readonly
@test m.create
@test isopen(m)
@test isreadable(m)
@test iswritable(m)

m = Mmap.mmap(Vector{UInt8}, 12)
@test length(m) == 12
@test all(m .== 0x00)
@test m[1] === 0x00
@test m[end] === 0x00
m[1] = 0x0a
Mmap.sync!(m)
@test m[1] === 0x0a
m = Mmap.mmap(Vector{UInt8}, 12; shared=false)
m = Mmap.mmap(Vector{Int}, 12)
@test length(m) == 12
@test all(m .== 0)
@test m[1] === 0
@test m[end] === 0
m = Mmap.mmap(Vector{Float64}, 12)
@test length(m) == 12
@test all(m .== 0.0)
m = Mmap.mmap(Matrix{Int8}, (12,12))
@test size(m) == (12,12)
@test all(m == zeros(Int8, (12,12)))
@test sizeof(m) == prod((12,12))
n = similar(m)
@test size(n) == (12,12)
n = similar(m, (2,2))
@test size(n) == (2,2)
n = similar(m, 12)
@test length(n) == 12
@test size(n) == (12,)
finalize(m); m = nothing; gc()

# test #14885
file = tempname()
touch(file)
open(file, "r+") do s
    A = Mmap.mmap(s, Vector{UInt8}, (10,), 0)
    Mmap.sync!(A)
    finalize(A); A = nothing; gc()
    A = Mmap.mmap(s, Vector{UInt8}, (10,), 1)
    Mmap.sync!(A)
    finalize(A); A = nothing; gc()
end
rm(file)

max_stars_repo_name: anowacki/julia

max_issues_repo_name: anowacki/julia

max_forks_repo_name: anowacki/julia


text:
using Test, YaoArrayRegister

@testset "test fidelity" begin
    reg = rand_state(3)
    reg_ = rand_state(3)
    reg2 = repeat(reg, 3)
    @test fidelity(reg, reg) ≈ [1]
    @test fidelity(reg, reg_)[] < 1
    @test fidelity(reg2, reg2) ≈ [1, 1, 1]

    # mix
    reg4 = join(reg, reg)
    reg5 = join(reg_, reg_)
    focus!(reg4, 1:3)
    focus!(reg5, 1:3)
    @test isapprox(fidelity(reg, reg_)[], fidelity(reg4, reg5)[], atol = 1e-5)
    @test isapprox.(fidelity(reg, reg_)[], fidelity(repeat(reg4, 3), repeat(reg5, 3)), atol = 1e-5) |>
          all
end

@testset "test trace distance" begin
    reg = rand_state(3)
    reg_ = rand_state(3)
    reg2 = repeat(reg, 3)
    dm = ρ(reg)
    dm_ = ρ(reg_)
    dm2 = ρ(reg2)
    @test reg |> probs ≈ dm |> probs
    @test isapprox(tracedist(dm, dm), tracedist(reg, reg), atol = 1e-5)
    @test isapprox(tracedist(dm, dm_), tracedist(reg, reg_), atol = 1e-5)
    @test isapprox(tracedist(dm2, dm2), tracedist(reg2, reg2), atol = 1e-5)

    # mix
    reg4 = join(reg, reg)
    reg5 = join(reg_, reg_)
    focus!(reg4, 1:3)
    focus!(reg5, 1:3)
    dm4 = reg4 |> density_matrix
    dm5 = reg5 |> density_matrix
    @test isapprox(tracedist(dm, dm_)[], tracedist(dm4, dm5)[], atol = 1e-5)
    @test isapprox.(
        tracedist(dm, dm_)[],
        tracedist(repeat(reg4, 3) |> density_matrix, repeat(reg5, 3) |> density_matrix),
        atol = 1e-5,
    ) |> all
end

@testset "purify" begin
    reg = rand_state(6)
    reg_p = purify(reg |> ρ)
    @test reg_p |> isnormalized
    @test reg_p |> exchange_sysenv |> probs |> maximum ≈ 1
    reg_p = purify(reg |> ρ; nbit_env = 0)
    @test Yao.fidelity(reg, reg_p) ≈ [1]

    reg = rand_state(6; nbatch = 10)
    reg_p = purify(reg |> ρ)
    @test reg_p |> isnormalized
    @test reg_p |> exchange_sysenv |> probs |> maximum ≈ 1
    reg_p = purify(reg |> ρ; nbit_env = 0)
    @test Yao.fidelity(reg, reg_p) ≈ ones(10)
    reg_p = purify(reg |> ρ; nbit_env = 2)
    @test reg_p |> nqubits == 8
end


text:
function coef_G(::Type{T}, rn::Integer) where T <: AbstractFloat
    G = ones(T, rn, rn)
    for j = 1:rn
        for i = j+1:2:rn
            G[i,j] = -G[i,j]
        end
    end
    return G
end

coef_G(rn) = coef_G(Float64, rn)

function coef_K(::Type{T}, rn::Integer, rnm1::Integer) where T <:AbstractFloat
    K = ones(T, rn, rnm1)
    for j = 1:rnm1
        for i = 2:2:rn
            K[i,j] = -K[i,j]
        end
    end
    return K
end

coef_K(rn, rnm1) = coef_K(Float64, rn, rnm1)

function coef_Hn!(r::Integer, n::Integer, ℓbar_hi::Integer, t::OffsetVector{T},
                 M::Integer, store::Store{T},
                 version=2) where T <: AbstractFloat
    Hn = OffsetVector{Matrix{T}}(undef, 0:ℓbar_hi)
    for ℓbar = 0:ℓbar_hi
        Hn[ℓbar] = Array{T}(undef, r, r)
    end
    coef_Hn!(Hn, n, ℓbar_hi, t, M, store, version)
    return Hn
end

function coef_Hn!(Hn::OffsetVector{Matrix{T}}, n::Integer, ℓbar_hi::Integer,
                 t::OffsetVector{T}, M::Integer,
                 store::Store{T}, version=2) where T <: AbstractFloat
    N = length(t) - 1
    @argcheck 1 ≤ n ≤ N
    @argcheck 0 ≤ ℓbar_hi ≤ n-1
    @argcheck length(Hn) ≥ ℓbar_hi + 1
    rn = size(Hn[0], 1)
    @argcheck size(Hn[0], 2) == rn
    @argcheck rn ≤ store.rmax
    coef_H0_uniform!(Hn[0], store)
    kn = t[n] - t[n-1]
    α = store.α
    for j = 1:rn, i = 1:rn
        Hn[0][i,j] *= kn^α
    end
    if ℓbar_hi ≥ 1
        coef_Hn1!(Hn[1], n, t, M, store)
    end
    if ℓbar_hi ≥ 2
        if version == 1
            coef_Hn_ver1!(Hn, n, 2:ℓbar_hi, t, M, store)
        elseif version == 2
            coef_Hn_ver2!(Hn, n, 2:ℓbar_hi, t, M, store)
        else
            throw(ArgumentError("version must be 1 or 2"))
        end
    end
end

function coef_Hn1!(rn::Integer, rnm1::Integer, n::Integer, t::OffsetVector{T},
                   M::Integer, store::Store{T}) where T <: AbstractFloat
    Hn1 = Array{T}(undef, rn, rnm1)
    coef_Hn1!(Hn1, n, t, M, store)
    return Hn1
end

function coef_Hn1!(Hn1::Matrix{T}, n::Integer, t::OffsetVector{T}, M::Integer,
                   store::Store{T}) where T <: AbstractFloat
    rn, rnm1 = size(Hn1)
    @argcheck rnm1 ≤ store.rmax
    A = view(store.A, 1:rnm1)
    B = view(store.B, 1:rnm1)
    C = view(store.C, 1:rn, 1:rnm1)
    α = store.α
    kn = t[n] - t[n-1]
    knm1 = t[n-1] - t[n-2]
    ρn = kn / knm1

    σ, wσ = rule(store.legendre[M])
    Ψ = view(store.Ψ, 1:rnm1, 1:M)
    legendre_polys!(Ψ, σ)
    c = (1+ρn)^(1-α)
    for j = 1:rnm1
        s = zero(T)
        for m = 1:M
            s += wσ[m] * ( 2ρn + 1 - σ[m] )^(α-1) * Ψ[j,m]
        end
        A[j] = c * s
    end

    Mσ = ceil(Integer, rnm1/2)
    σ, wσ = rule(store.jacobi1[Mσ])
    Ψ = view(store.Ψ, 1:rnm1, 1:Mσ)
    legendre_polys!(Ψ, σ)
    for j = 1:rnm1
        s = zero(T)
        for m = 1:Mσ
            s+= wσ[m] * Ψ[j,m]
        end
        B[j] = c * s
    end

    Mτ = ceil(Integer, rn-1)
    Mσ = Mτ
    Mz = M
    z, wz = rule(store.unitlegendre[Mz])
    σ, wσ = rule(store.jacobi3[Mσ])
    τ, wτ = rule(store.jacobi4[Mτ])
    dΨ = view(store.dΨ, 1:rn, 1:Mτ)
    deriv_legendre_polys!(dΨ, τ)
    for j = 1:rnm1
        Ψ = view(store.Ψ, 1:j)
        for i = 1:rn
            outer = zero(T)
            for mτ = 1:Mτ
                τ_ = τ[mτ]
                inner = zero(T)
                for mz = 1:Mz
                    z_ = z[mz]
                    legendre_polys!(Ψ, 1 - z_*(1+τ_))
                    inner += wz[mz] * ( ρn + z_ )^(α-1) * Ψ[j]
                end
                outer += wτ[mτ] * dΨ[i,mτ] * inner
            end
            C[i,j] = c * outer
        end
    end
    Ψ = view(store.Ψ, 1:rnm1, 1:Mσ)
    legendre_polys!(Ψ, σ)
    for j = 1:rnm1
        for i = 1:rn
            dΨ = view(store.dΨ, 1:i)
            outer = zero(T)
            for mσ = 1:Mσ
                σ_ = σ[mσ]
                inner = zero(T)
                for mz = 1:Mz
                    z_ = z[mz]
                    deriv_legendre_polys!(dΨ, z_*(1-σ_)-1)
                    inner += wz[mz] * ( ρn*z_ + 1 )^(α-1) * dΨ[i]
                end
                outer += wσ[mσ] * Ψ[j,mσ] * inner
            end
            C[i,j] += c * outer
        end
    end

    D = ( kn + knm1 ) / 2
    c = D^(α-1) * knm1 / ( 2Γ(α) )
    for j = 1:rnm1
        pow = one(T)
        for i = 1:rn
            pow = -pow
            Hn1[i,j] = c * ( A[j] + pow * B[j] - C[i,j] )
        end
    end
end

function coef_Hn_ver1!(Hn::OffsetVector{Matrix{T}}, n::Integer, rng::UnitRange,
                       t::OffsetVector{T}, M::Integer,
                       store::Store{T}) where T <: AbstractFloat
    @argcheck rng.start ≥ 2
    rn = size(Hn[rng.start], 1)
    kn = t[n] - t[n-1]

    σ, wσ = rule(store.legendre[M])
    τ, wτ = σ, wσ
    Ψ  = view(store.Ψ,  1:store.rmax, 1:M)
    dΨ = view(store.dΨ, 1:store.rmax, 1:M)
    legendre_polys!(Ψ, σ)
    deriv_legendre_polys!(dΨ, τ)
    α = store.α
    c = 1 / ( 2Γ(α) )
    for ℓbar in rng
        rℓ = size(Hn[ℓbar], 2)
        @argcheck rℓ ≤ store.rmax
        ℓ = n - ℓbar
        kℓ = t[ℓ] - t[ℓ-1]
        two_D = (t[n]+t[n-1]) - (t[ℓ]+t[ℓ-1])

        A = view(store.A, 1:rℓ)
        B = view(store.B, 1:rℓ)
        C = view(store.C, 1:rn, 1:rℓ)

        for j = 1:rℓ
            s = zero(T)
            for m = 1:M
                Δ = ( kn - σ[m]*kℓ ) / two_D
                s += wσ[m] * ( 1 + Δ )^(α-1) * Ψ[j,m]
            end
            A[j] = s
        end

        for j = 1:rℓ
            s = zero(T)
            for m = 1:M
                Δ = ( -kn - σ[m]*kℓ ) / two_D
                s += wσ[m] * ( 1 + Δ )^(α-1) * Ψ[j,m]
            end
            B[j] = s
        end

        for j = 1:rℓ, i = 1:rn
            outer = zero(T)
            for mτ = 1:M
                inner = zero(T)
                for mσ = 1:M
                    Δ = ( τ[mτ]*kn - σ[mσ]*kℓ ) / two_D
                    inner += wσ[mσ] * (1+Δ)^(α-1) * Ψ[j,mσ]
                end
                outer += wτ[mτ] * dΨ[i,mτ] * inner
            end
            C[i,j] = outer
        end

        D = two_D / 2
        for j = 1:rℓ
            pow = one(T)
            for i = 1:rn
                pow = -pow
                Hn[ℓbar][i,j] = c * D^(α-1) * kℓ * (
                                A[j] + pow * B[j] - C[i,j] )
            end
        end
    end
end

function coef_Hn_ver2!(Hn::OffsetVector{Matrix{T}}, n::Integer, rng::UnitRange,
                       t::OffsetVector{T}, M::Integer,
                       store::Store{T}) where T <: AbstractFloat
    @argcheck rng.start ≥ 2
    rn = size(Hn[rng.start], 1)
    kn = t[n] - t[n-1]

    σ, wσ = rule(store.legendre[M])
    τ, wτ = σ, wσ
    Ψ  = view(store.Ψ,  1:store.rmax, 1:M)
    legendre_polys!(Ψ, σ)
    α = store.α
    c = - (1-α) / ( 4Γ(α) )
    for ℓbar in rng
        rℓ = size(Hn[ℓbar], 2)
        @argcheck rℓ ≤ store.rmax
        ℓ = n - ℓbar
        kℓ = t[ℓ] - t[ℓ-1]
        two_D = (t[n]+t[n-1]) - (t[ℓ]+t[ℓ-1])
        D = two_D / 2

        for j = 1:rℓ, i = 1:rn
            outer = zero(T)
            for mτ = 1:M
                inner = zero(T)
                for mσ = 1:M
                    Δ = ( τ[mτ]*kn - σ[mσ]*kℓ ) / two_D
                    inner += wσ[mσ] * (1+Δ)^(α-2) * Ψ[j,mσ]
                end
                outer += wτ[mτ] * Ψ[i,mτ] * inner
            end
            Hn[ℓbar][i,j] = c * kn * kℓ * D^(α-2) * outer
        end
    end
end


text:
# Set up the environment to run this example. Make sure you're within the folder that this
# file lives in.
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

using Stheno, Plots, Random
gr()

# Define model.
σ², ω, T = 1e0, 1.0, 25.0
gpc = GPC()
f₁ = periodic(GP(EQ(), gpc), ω)
f₂ = GP(0.1 * EQ(), gpc)
f₃ = f₁ + f₂

# Sample from marginal process to generate toy data.
rng = MersenneTwister(123456)
S = 25
x = range(0.0, T; length=100)
xp = range(-2.5, T + 2.5; length=500)
fx = f₃(x, σ²)
y = rand(rng, fx)

# Plots for stuff.
f₁′_plot, f₂′_plot, f₃′_plot = plot(), plot(), plot();

posterior_plot = plot();
approx_in_marginal_posterior_plot = plot();
approx_in_latents_posterior_plot = plot();

##################################### Exact Inference ######################################

# Compute the posterior processes, sample from them, and compute marginals.
@show logpdf(fx, y)
f₁′, f₂′, f₃′ = (f₁, f₂, f₃) | (fx ← y)
f₁′xp, f₂′xp, f₃′xp = rand(rng, [f₁′(xp, 1e-6), f₂′(xp, 1e-6), f₃′(xp, 1e-6)], S)

ms₁ = marginals(f₁′(xp))
ms₂ = marginals(f₂′(xp))
ms₃ = marginals(f₃′(xp))
μ₁′, σ₁′ = mean.(ms₁), std.(ms₁)
μ₂′, σ₂′ = mean.(ms₂), std.(ms₂)
μ₃′, σ₃′ = mean.(ms₃), std.(ms₃)

items = [
    (μ₁′, σ₁′, f₁′xp, :red, "exact", f₁′_plot),
    (μ₂′, σ₂′, f₂′xp, :red, "", f₂′_plot),
    (μ₃′, σ₃′, f₃′xp, :red, "", f₃′_plot),
]

# Posterior marginal variance.
for (μ, σ, _, colour, name, plt) in items
    plot!(plt, xp, [μ, μ];
        linewidth=0.0,
        fillrange=[μ .- 3 .* σ, μ .+ 3 * σ],
        fillalpha=0.3,
        fillcolor=colour,
        label="",
    );
end

# Posterior samples.
for (μ, σ, f, colour, name, plt) in items
    plot!(plt, xp, f,
        linecolor=colour,
        linealpha=0.2,
        label="",
    );
end

# Posterior mean.
for (μ, σ, f, colour, name, plt) in items
    plot!(plt, xp, μ;
        linecolor=colour,
        linewidth=2.0,
        label=name,
    );
end

# Plot observations.
scatter!(posterior_plot, x, y;
    markercolor=:blue,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.7,
    label="",
);

##################### Approximate inference with pseudo-data in f₃ #########################

# Compute approximate posterior process.
M₃ = 15;
Z₃ = range(0.0, T; length=M₃);
Z₃ = 0:T
Z₃ = 0:0.5:T
# Z₃ = 0:0.25:
u₃ = f₃(Z₃);
# μ′u₃, Σ′u₃ = Stheno.optimal_q(f₃(X), ŷ, u₃, σ_noise);
# q_u₃ = Stheno.Titsias(u₃, μ′u₃, Σ′u₃, gpc);
pseudo_obs = Stheno.PseudoObs(f₃(x) ← y, u₃)
f₁′u₃ = f₁ | pseudo_obs
f₂′u₃ = f₂ | pseudo_obs
f₃′u₃ = f₃ | pseudo_obs
f₁′u₃xp, f₂′u₃xp, f₃′u₃xp = rand(rng, [f₁′u₃(xp), f₂′u₃(xp), f₃′u₃(xp)], S);

μ₁′u₃, σ₁′u₃ = mean.(marginals(f₁′u₃(Xp))), std.(marginals(f₁′u₃(Xp)))
μ₂′u₃, σ₂′u₃ = mean.(marginals(f₂′u₃(Xp))), std.(marginals(f₂′u₃(Xp)))
μ₃′u₃, σ₃′u₃ = mean.(marginals(f₃′u₃(Xp))), std.(marginals(f₃′u₃(Xp)))

@show elbo(fx, y, u₃);

items = [
    (μ₁′u₃, σ₁′u₃, f₁′u₃xp, :green, "Z in f3", f₁′_plot),
    (μ₂′u₃, σ₂′u₃, f₂′u₃xp, :green, "", f₂′_plot),
    (μ₃′u₃, σ₃′u₃, f₃′u₃xp, :green, "", f₃′_plot),
];

# Posterior marginal variance.
for (μ, σ, _, colour, name, plt) in items
    plot!(plt, Xp, [μ, μ];
        linewidth=0.0,
        fillrange=[μ .- 3 .* σ, μ .+ 3 * σ],
        fillalpha=0.3,
        fillcolor=colour,
        label="");
end

# # Posterior samples.
# for (μ, σ, f, colour, name, plt) in items
#     plot!(plt, Xp, f,
#         linecolor=colour,
#         linealpha=0.2,
#         label="");
# end

# Posterior mean.
for (μ, σ, f, colour, name, plt) in items
    plot!(plt, Xp, μ;
        linecolor=colour,
        linewidth=2.0,
        label=name,
    );
end

# Plot observations and pseudo-input locations.
scatter!(approx_in_marginal_posterior_plot, X, ŷ;
    markercolor=:blue,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.7,
    label="",
);
scatter!(approx_in_marginal_posterior_plot, Z₃, zeros(M₃);
    markercolor=:black,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.8,
    label="Z",
);

############# Perform approximate inference by placing pseudo-data in f₁ and f₂ ############

# Compute approximate posterior process.
M₁, M₂ = 10, 5
Z₁, Z₂ = linspace(0.0, 1 / ω, M₁), linspace(0.0, T, M₂)
u₁₂ = BlockGP([f₁(Z₁), f₂(Z₂)])
μ′u, Σ′u = Stheno.optimal_q(f₃(X), ŷ, u₁₂, σ_noise)
conditioner = Stheno.Titsias(u₁₂, μ′u, Σ′u, gpc)
f₁′u₁₂, f₂′u₁₂, f₃′u₁₂ = (f₁, f₂, f₃) | conditioner

f₁′u₁₂Xp, f₂′u₁₂Xp, f₃′u₁₂Xp = rand(rng, [f₁′u₁₂(Xp), f₂′u₁₂(Xp), f₃′u₁₂(Xp)], S)
μ₁′u₁₂, σ₁′u₁₂ = marginals(f₁′u₁₂(Xp))
μ₂′u₁₂, σ₂′u₁₂ = marginals(f₂′u₁₂(Xp))
μ₃′u₁₂, σ₃′u₁₂ = marginals(f₃′u₁₂(Xp))

@show elbo(f₃(X), ŷ, u₁₂, σ_noise)

items = [
    (μ₁′u₁₂, σ₁′u₁₂, f₁′u₁₂Xp, :blue, "Z in f1 and f2", f₁′_plot),
    (μ₂′u₁₂, σ₂′u₁₂, f₂′u₁₂Xp, :blue, "", f₂′_plot),
    (μ₃′u₁₂, σ₃′u₁₂, f₃′u₁₂Xp, :blue, "", f₃′_plot),
]

# Posterior marginal variance.
for (μ, σ, _, colour, name, plt) in items
    plot!(plt, Xp, [μ, μ];
        linewidth=0.0,
        fillrange=[μ .- 3 .* σ, μ .+ 3 * σ],
        fillalpha=0.3,
        fillcolor=colour,
        label="",
    )
end

# # Posterior samples.
# for (μ, σ, f, colour, name, plt) in items
#     plot!(plt, Xp, f,
#         linecolor=colour,
#         linealpha=0.2,
#         label="");
# end

# Posterior mean.
for (μ, σ, f, colour, name, plt) in items
    plot!(plt, Xp, μ;
        linecolor=colour,
        linewidth=2.0,
        label=name,
    )
end

# Plot observations and pseudo-input locations.
scatter!(approx_in_latents_posterior_plot, X, ŷ;
    markercolor=:blue,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.7,
    label="",
);
scatter!(approx_in_latents_posterior_plot, Z₁, zeros(M₁);
    markercolor=:black,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.8,
    label="Z₁",
);
scatter!(approx_in_latents_posterior_plot, Z₂, zeros(M₂);
    markercolor=:magenta,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.8,
    label="Z₂",
);

# display(posterior_plot);
# display(approx_in_marginal_posterior_plot);
# display(approx_in_latents_posterior_plot);

plot(f₁′_plot, f₂′_plot, f₃′_plot; layout=(3, 1))