text:
# ------------------------------------------------------------------------------------------
# # Introduction to Queryverse.jl
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# This notebook has all the code required to follow the tutorial by [David
# Anthoff](https://www.david-anthoff.com/) on 6/14/2018. This notebook _is not_ a self
# contained introduction to the
# [Queryverse.jl](https://github.com/davidanthoff/Queryverse.jl), instead it is meant to
# complement the video tutorial.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Logistics
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ![alt](slides/Slide3.PNG)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ![alt](slides/Slide4.PNG)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## End-to-end example
# ------------------------------------------------------------------------------------------

using Queryverse, VegaDatasets, IndexedTables

# ------------------------------------------------------------------------------------------
# The `dataset` function from the
# [VegaDatasets.jl](https://github.com/davidanthoff/VegaDatasets.jl) package returns a
# number of useful example datasets.
# ------------------------------------------------------------------------------------------

cars = dataset("cars")

# ------------------------------------------------------------------------------------------
# __WARNING__: The following cell will not run on juliabox, you need to run it on a local
# Jupyter instance.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Piping the `cars` dataset into
# [DataVoyager.jl](https://github.com/davidanthoff/DataVoyager.jl) opens a UI for data
# exploration and plotting.
# ------------------------------------------------------------------------------------------

cars |> Voyager()

# ------------------------------------------------------------------------------------------
# The default view shows a plot per column, picking an appropriate style based on the data
# type of the column.
# ![alt](screenshots/voyager1.png)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Dragging the quantitative wildcard fields on both the x and y channel shows as correlation
# plots.
# ![alt](screenshots/voyager2.png)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We can use [Query.jl](https://github.com/davidanthoff/Query.jl) operators to filter the
# `cars` dataset.
# ------------------------------------------------------------------------------------------

cars |>
@filter(_.Origin=="USA")

# ------------------------------------------------------------------------------------------
# The file IO features of the Queryverse allow us to save the filtered dataset as a CSV
# file.
# ------------------------------------------------------------------------------------------

cars |>
@filter(_.Origin=="USA") |>
save("us_cars.csv")

# ------------------------------------------------------------------------------------------
# We can also load CSV files from disc, and directly create a plot for that data with
# [VegaLite.jl](https://github.com/fredo-dedup/VegaLite.jl).
# ------------------------------------------------------------------------------------------

load("us_cars.csv") |>
@vlplot(:point, x=:Miles_per_Gallon, y=:Weight_in_lbs, color="Cylinders:n")

# ------------------------------------------------------------------------------------------
# ## File IO
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Queryverse provides a unified file IO API for a large number of tabular file formats.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ![alt](slides/Slide7.PNG)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Load
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ![alt](slides/Slide8.PNG)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Load a CSV file:
# ------------------------------------------------------------------------------------------

load("data/cars.csv")

# ------------------------------------------------------------------------------------------
# Load a feather file:
# ------------------------------------------------------------------------------------------

load("data/cars.feather")

# ------------------------------------------------------------------------------------------
# Load an Excel file (note how we have to specify the sheet we want to load):
# ------------------------------------------------------------------------------------------

load("data/cars.xlsx", "cars")

# ------------------------------------------------------------------------------------------
# We can also specify the exact range in the Excel file we want to load:
# ------------------------------------------------------------------------------------------

load("data/cars.xlsx", "cars!C1:G41")

# ------------------------------------------------------------------------------------------
# We can also load various stats files:
# ------------------------------------------------------------------------------------------

load("data/cars.dta")

load("data/cars.sas7bdat")

load("data/cars.sav")

# ------------------------------------------------------------------------------------------
# ### Save
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ![alt](slides/Slide9.PNG)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Remember, the `cars` variable is a table we previously loaded:
# ------------------------------------------------------------------------------------------

cars

# ------------------------------------------------------------------------------------------
# We can save it as a CSV file with this syntax:
# ------------------------------------------------------------------------------------------

save("cars1.csv", cars)

# ------------------------------------------------------------------------------------------
# We can also pipe the data into the `save` function:
# ------------------------------------------------------------------------------------------

cars |> save("cars2.csv")

# ------------------------------------------------------------------------------------------
# Depending on the file format, we can specify various options, for example the delimiter
# character for CSV files:
# ------------------------------------------------------------------------------------------

cars |> save("cars3.csv", delim=';')

# ------------------------------------------------------------------------------------------
# We can also save as a feather file:
# ------------------------------------------------------------------------------------------

cars |> save("cars4.feather")

# ------------------------------------------------------------------------------------------
# ### Combine load and save
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We can also load a file in one format and directly pipe it into the `save` function and
# save it in another format. Note that the data will never be materialized as say a
# `DataFrame` in that case.
# ------------------------------------------------------------------------------------------

load("data/cars.sas7bdat") |> save("cars5.csv")

# ------------------------------------------------------------------------------------------
# ## Tables
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ![alt](slides/Slide11.PNG)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ![alt](slides/Slide12.PNG)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ![alt](slides/Slide13.PNG)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Note how the load function returns a different type for each kind of file we can load. All
# of them implement the TableTraits.jl interface.
# ------------------------------------------------------------------------------------------

load("data/cars.csv") |> typeof

load("data/cars.feather") |> typeof

dataset("cars") |> typeof

# ------------------------------------------------------------------------------------------
# The constructors of most table types accept any table, and will automatically convert the
# table at construction. For example, the following code loads a feather file. The call to
# `load` returns a `FeatherFile` instance, which implements the TableTraits.jl interface.
# When we pass that to the `DataFrame` constructor, it automatically constructs a
# `DataFrame` from the data in the file.
# ------------------------------------------------------------------------------------------

DataFrame(load("data/cars.feather"))

# ------------------------------------------------------------------------------------------
# We can also pipe things into a `DataFrame`:
# ------------------------------------------------------------------------------------------

load("data/cars.feather") |> DataFrame

# ------------------------------------------------------------------------------------------
# We can also pipe things into an `IndexedTable`:
# ------------------------------------------------------------------------------------------

load("data/cars.feather") |> table

# ------------------------------------------------------------------------------------------
# Because both `DataFrame` and `IndexedTable` implement the TableTraits.jl interface, we can
# also do silly things like 1) load a feather file, 2) convert it into an `IndexedTable`,
# then 3) convert that `IndexedTable` into a `DataFrame`.
# ------------------------------------------------------------------------------------------

load("data/cars.feather") |> table |> DataFrame

# ------------------------------------------------------------------------------------------
# It works the other way around as well:
# ------------------------------------------------------------------------------------------

load("data/cars.feather") |> DataFrame |> table

# ------------------------------------------------------------------------------------------
# Because [VegaLite.jl](https://github.com/fredo-dedup/VegaLite.jl) accepts data as any
# TableTraits.jl table, we can plot `DataFrame`s:
# ------------------------------------------------------------------------------------------

load("data/cars.feather") |> DataFrame |> @vlplot(:point, x=:Horsepower, y=:Acceleration)

# ------------------------------------------------------------------------------------------
# Or we can plot `IndexedTable`s
# ------------------------------------------------------------------------------------------

load("data/cars.feather") |> table |> @vlplot(:point, x=:Horsepower, y=:Acceleration)

# ------------------------------------------------------------------------------------------
# Or we can skip the table types entirely and directly plot the file:
# ------------------------------------------------------------------------------------------

load("data/cars.feather") |> @vlplot(:point, x=:Horsepower, y=:Acceleration)

# ------------------------------------------------------------------------------------------
# We can also do lots of other combinations. Here we load a stats file into an
# `IndexedTable`:
# ------------------------------------------------------------------------------------------

it = load("data/cars.sas7bdat") |> table

# ------------------------------------------------------------------------------------------
# Then we convert it into a `DataFrame`:
# ------------------------------------------------------------------------------------------

df = it |> DataFrame

# ------------------------------------------------------------------------------------------
# We can pipe the `IndexedTable` into Voyager:
# ------------------------------------------------------------------------------------------

it |> Voyager()

# ------------------------------------------------------------------------------------------
# Or we can pipe the `DataFrame` into Voyager:
# ------------------------------------------------------------------------------------------

df |> Voyager()

# ------------------------------------------------------------------------------------------
# All the file IO functions work with anything that implements the TableTraits.jl interface,
# so we can e.g. save both an `IndexedTable` and a `DataFrame` as a feather file:
# ------------------------------------------------------------------------------------------

it |> save("cars8.feather")

df |> save("cars9.feather")

# ------------------------------------------------------------------------------------------
# ![alt](slides/Slide14.PNG)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## The pipe
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ![alt](slides/Slide16.PNG)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We can pipe in one line:
# ------------------------------------------------------------------------------------------

load("data/cars.csv") |> @filter(_.Origin=="Europe") |> DataFrame

# ------------------------------------------------------------------------------------------
# Or we can write each step on its own line:
# ------------------------------------------------------------------------------------------

load("data/cars.csv") |>
    @filter(_.Origin=="Europe") |>
    DataFrame

# ------------------------------------------------------------------------------------------
# The following example doesn't return anything at the end because the `save` function has
# no return value:
# ------------------------------------------------------------------------------------------

load("data/cars.csv") |>
    @filter(_.Origin=="Europe") |>
    save("cars10.csv")

# ------------------------------------------------------------------------------------------
# ![alt](slides/Slide17.PNG)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# But sometimes we want to continue the pipe, for example to display the table at the
# interactive prompt:
# ------------------------------------------------------------------------------------------

load("data/cars.csv") |>
    @filter(_.Origin=="Europe") |>
    @tee(save("cars10.csv"))

# ------------------------------------------------------------------------------------------
# Or to save the data as a CSV File and then create a plot:
# ------------------------------------------------------------------------------------------

load("data/cars.csv") |>
    @filter(_.Origin=="Europe") |>
    @tee(save("cars10.csv")) |>
    @vlplot(:point, x=:Acceleration, y=:Horsepower)

# ------------------------------------------------------------------------------------------
# We can also pipe the plot to the `save` function, i.e. it is not just for tabular data!
# ------------------------------------------------------------------------------------------

load("data/cars.csv") |>
    @filter(_.Origin=="Europe") |>
    @tee(save("cars10.csv")) |>
    @vlplot(:point, x=:Acceleration, y=:Horsepower) |>
    save("fig1.png")

# ------------------------------------------------------------------------------------------
# If we want to save the plot in multiple formats and see it at the end, we can use repeated
# `@tee` calls:
# ------------------------------------------------------------------------------------------

load("data/cars.csv") |>
    @filter(_.Origin=="Europe") |>
    @tee(save("cars10.csv")) |>
    @vlplot(:point, x=:Acceleration, y=:Horsepower) |>
    @tee(save("fig1.png")) |>
    @tee(save("fig1.pdf"))

# ------------------------------------------------------------------------------------------
# ## Queries
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Basics
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ![alt](slides/Slide19.PNG)
# ------------------------------------------------------------------------------------------

load("data/cars.feather") |> @filter(_.Origin=="Europe") |> @orderby(_.Horsepower)

# ------------------------------------------------------------------------------------------
# Queries are not restricted to tabular data.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ![alt](slides/Slide20.PNG)
# ------------------------------------------------------------------------------------------

1:8 |> @filter(_%2==0) |> @orderby_descending(_)

# ------------------------------------------------------------------------------------------
# If we pipe things into `collect` at the end, we generate an array:
# ------------------------------------------------------------------------------------------

1:8 |> @filter(_%2==0) |> @orderby_descending(_) |> collect

# ------------------------------------------------------------------------------------------
# We can also start from a `Dict`:
# ------------------------------------------------------------------------------------------

Dict(:a=>4, :b=>2, :c=>8)

# ------------------------------------------------------------------------------------------
# And filter the `Dict`:
# ------------------------------------------------------------------------------------------

Dict(:a=>4, :b=>2, :c=>8) |> @filter(_[2]>3)

# ------------------------------------------------------------------------------------------
# And collect the pairs into an array:
# ------------------------------------------------------------------------------------------

Dict(:a=>4, :b=>2, :c=>8) |> @filter(_[2]>3) |> collect

# ------------------------------------------------------------------------------------------
# Or we can collec them into another `Dict`. Not that this will only work if the pipe
# generates a stream of `Pair`s:
# ------------------------------------------------------------------------------------------

Dict(:a=>4, :b=>2, :c=>8) |> @filter(_[2]>3) |> Dict

# ------------------------------------------------------------------------------------------
# We can also start from a generator:
# ------------------------------------------------------------------------------------------

(i=>i^2 for i in 1:10) |> @filter(_[2]>3) |> @orderby(_[1]) |> collect

# ------------------------------------------------------------------------------------------
# ### @filter
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# `@filter` removes elements from the input stream based on some condition.
# ------------------------------------------------------------------------------------------

1:10 |> @filter(_>5)

# ------------------------------------------------------------------------------------------
# ![alt](slides/Slide21.PNG)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# You don't have to use the anonymous function shorthand, you can also use normal anonymous
# functions (but why would you, more typing!):
# ------------------------------------------------------------------------------------------

1:10 |> @filter(i -> i>5)

# ------------------------------------------------------------------------------------------
# ### @map
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# `@map` applies a transformation function to each element from the input stream. Here we
# square every element:
# ------------------------------------------------------------------------------------------

1:10 |> @map(_^2)

# ------------------------------------------------------------------------------------------
# Here we create a `Pair` for each input element, so that we can pipe things into a `Dict`:
# ------------------------------------------------------------------------------------------

1:10 |> @map(_=>_^2) |> Dict

# ------------------------------------------------------------------------------------------
# Here we extract the `Origin` column for the `cars` table:
# ------------------------------------------------------------------------------------------

cars |> @map(_.Origin)

# ------------------------------------------------------------------------------------------
# #### Brief detour about named tuples
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# First, lets store the stream of named tuples in an array:
# ------------------------------------------------------------------------------------------

r = cars |> collect

# ------------------------------------------------------------------------------------------
# We then access the first named tuple from that array:
# ------------------------------------------------------------------------------------------

nt = r[1]

# ------------------------------------------------------------------------------------------
# And look at the type:
# ------------------------------------------------------------------------------------------

typeof(nt)

# ------------------------------------------------------------------------------------------
# We can use the `fieldnames` function to inquire about the names of the fields of the named
# tuple:
# ------------------------------------------------------------------------------------------

fieldnames(nt)

# ------------------------------------------------------------------------------------------
# We can access a field with the `.` notation:
# ------------------------------------------------------------------------------------------

nt.Origin

# ------------------------------------------------------------------------------------------
# Or by indexing, the way we would index into a normal tuple:
# ------------------------------------------------------------------------------------------

nt[3]

# ------------------------------------------------------------------------------------------
# So the example from above extracts the `Origin` field for each row it gets, at which point
# we have a stream of `String`s.
# ------------------------------------------------------------------------------------------

cars |> @map(_.Origin)

# ------------------------------------------------------------------------------------------
# ![alt](slides/Slide22.PNG)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We can construct named tuples with the `{}` syntax. Note that this only works within
# Query.jl macros. Whenever we create named tuples in a `@map`, we are creating a table,
# because a table is a stream of named tuples:
# ------------------------------------------------------------------------------------------

1:10 |> @map({foo=_, bar=_^2})

# ------------------------------------------------------------------------------------------
# Lets remember the columns of the `cars` dataset:
# ------------------------------------------------------------------------------------------

cars

# ------------------------------------------------------------------------------------------
# We can extract only two columns:
# ------------------------------------------------------------------------------------------

cars |> @map({_.Name, _.Year})

# ------------------------------------------------------------------------------------------
# We can extract two columns and rename one:
# ------------------------------------------------------------------------------------------

cars |> @map({_.Name, Foo=_.Year})

# ------------------------------------------------------------------------------------------
# Because we are producing a table with the `@map` command, we can save things as a CSV
# file:
# ------------------------------------------------------------------------------------------

cars |> @map({_.Name, Foo=_.Year}) |> save("cars11.csv")

# ------------------------------------------------------------------------------------------
# We can also construct `Pair`s from a table and then store things as a `Dict`:
# ------------------------------------------------------------------------------------------

cars |> @map(_.Name => Date(_.Year)) |> Dict

# ------------------------------------------------------------------------------------------
# ### @take and @drop
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# `@drop(n)` drops the first `n` elements from the input stream:
# ------------------------------------------------------------------------------------------

1:10 |> @drop(3)

# ------------------------------------------------------------------------------------------
# `@take(n)` takes the first `n` elements from the input stream:
# ------------------------------------------------------------------------------------------

1:10 |> @drop(3) |> @take(4)

# ------------------------------------------------------------------------------------------
# That of course also works with tables. To view only row 4-7 we can do:
# ------------------------------------------------------------------------------------------

cars |> @drop(3) |> @take(4)

# ------------------------------------------------------------------------------------------
# ### Sorting
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# The `@orderby` operators accepts an anonymous function that extracts a key for each input
# element. The elements are then sorted by that key. For example, to sort the `cars` dataset
# by the `Origin` column, we can write this code:
# ------------------------------------------------------------------------------------------

cars |> @orderby(_.Origin)

# ------------------------------------------------------------------------------------------
# We can also sort in descending order:
# ------------------------------------------------------------------------------------------

cars |> @orderby_descending(_.Origin)

# ------------------------------------------------------------------------------------------
# There are lots of rows that have the same value for `Origin`. We can use the `@thenby`
# command to specify a secondary sort key:
# ------------------------------------------------------------------------------------------

cars |> @orderby(_.Origin) |> @thenby(_.Year)

# ------------------------------------------------------------------------------------------
# You can chain as many `@thenby` operators as you want:
# ------------------------------------------------------------------------------------------

cars |> @orderby(_.Origin) |> @thenby(_.Year) |> @thenby_descending(_.Horsepower)

# ------------------------------------------------------------------------------------------
# The key extraction function is not restricted to just returning a value, it can also call
# arbitrary julia functions. For example, here we sort the `cars` dataset by the character
# length of the `Origin` column:
# ------------------------------------------------------------------------------------------

cars |> @orderby(length(_.Origin))

# ------------------------------------------------------------------------------------------
# ### @groupby
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ![alt](slides/Slide23.PNG)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Here we group by the `Origin` column:
# ------------------------------------------------------------------------------------------

cars |> @groupby(_.Origin)

# ------------------------------------------------------------------------------------------
# We can look at the keys of each group:
# ------------------------------------------------------------------------------------------

cars |> @groupby(_.Origin) |> @map(_.key)

# ------------------------------------------------------------------------------------------
# We can also look at how many elements there are in each group:
# ------------------------------------------------------------------------------------------

cars |> @groupby(_.Origin) |> @map(length(_))

# ------------------------------------------------------------------------------------------
# We can also extract e.g. the third element for each group:
# ------------------------------------------------------------------------------------------

cars |> @groupby(_.Origin) |> @map(_[3])

# ------------------------------------------------------------------------------------------
# In the next example we are creating a table that shows us how many elements we have in
# each group as a table:
# ------------------------------------------------------------------------------------------

cars |> @groupby(_.Origin) |> @map({Origin=_.key, Count=length(_)})

# ------------------------------------------------------------------------------------------
# ![alt](slides/Slide24.PNG)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# The anonymous function that returns the group key can call any julia function. For
# example, we can group things by the length of the `Origin` column:
# ------------------------------------------------------------------------------------------

cars |> @groupby(length(_.Origin)) |> @map({OriginLength=_.key, Count=length(_)})

# ------------------------------------------------------------------------------------------
# ![alt](slides/Slide25.PNG)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We can also supply a second anonymous function that gets applied to each element from the
# source before that element gets placed into a group. This can be helpful if we want to
# take e.g. the mean for each group:
# ------------------------------------------------------------------------------------------

cars |> @groupby(_.Origin, _.Acceleration) |> @map({Origin=_.key, MeanAcceleration=mean(_)})

# ------------------------------------------------------------------------------------------
# ![alt](slides/Slide26.PNG)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ![alt](slides/Slide27.PNG)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# The following example shows an alternative way to computing the mean of a given column for
# each group:
# ------------------------------------------------------------------------------------------

cars |> @groupby(_.Origin) |> @map({Origin=_.key, MeanAcceleration=mean(_..Acceleration)})

# ------------------------------------------------------------------------------------------
# This alternative style is especially useful if we want to compute aggregates for multiple
# columns:
# ------------------------------------------------------------------------------------------

region_stats = cars |>
    @groupby(_.Origin) |>
    @map({
        Region=_.key,
        MeanAcceleration=mean(_..Acceleration),
        MinCylinders=minimum(_..Cylinders)
    })

# ------------------------------------------------------------------------------------------
# ### @join
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Joins combine elements from two sources into one stream of elements.
# ------------------------------------------------------------------------------------------

cars

region_stats

# ------------------------------------------------------------------------------------------
# Here we pick all row combinations where the `Origin` and `Region` column match, and then
# do some computations:
# ------------------------------------------------------------------------------------------

cars |>
    @join(region_stats, _.Origin, _.Region, {_.Name, reg_ex_acc=_.Acceleration - __.MeanAcceleration}) |>
    @orderby_descending(_.reg_ex_acc)

# ------------------------------------------------------------------------------------------
# ### When is processing happening?
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# The following line does zero data processing because we don't display the results (note
# the `;` at the end of the line):
# ------------------------------------------------------------------------------------------

q = load("data/cars.feather") |>
    @filter(_.Origin=="USA") |>
    @map({_.Name, _.Cylinders}) |>
    @take(5);

# ------------------------------------------------------------------------------------------
# Only when we start to iterate the results is any data read from disc and processed:
# ------------------------------------------------------------------------------------------

for row in q
    println(row)
end

# ------------------------------------------------------------------------------------------
# ![alt](slides/Slide28.PNG)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# The following example is completely streaming. Each row is read from disc, filtered,
# projected and saved into the CSV file before the next row is read. In some sense this is
# resembling loop fusion.
# ------------------------------------------------------------------------------------------

load("data/cars.feather") |>
    @filter(_.Origin=="USA") |>
    @map({_.Name, _.Cylinders}) |>
    save("cars12.csv")

# ------------------------------------------------------------------------------------------
# ## Plotting & visual exploration
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Lets look at the `cars` dataset again:
# ------------------------------------------------------------------------------------------

cars

# ------------------------------------------------------------------------------------------
# We can pipe it into a plot. We must first specify the kind of `mark` we want to use. In
# this example we are using a point mark. Vega-Lite will draw a point for each row in our
# dataset. Because we haven't specified any encodings, all points are drawn on top of each
# other:
# ------------------------------------------------------------------------------------------

cars |> @vlplot(:point)

# ------------------------------------------------------------------------------------------
# We can start to make the plot more useful by encoding the `x` position of the points to
# the `Miles_per_Gallon` column:
# ------------------------------------------------------------------------------------------

cars |> @vlplot(:point, x=:Miles_per_Gallon)

# ------------------------------------------------------------------------------------------
# If we also encode the `y` channel of the mark, we start to get a useful plot:
# ------------------------------------------------------------------------------------------

cars |> @vlplot(:point, x=:Miles_per_Gallon, y=:Acceleration)

# ------------------------------------------------------------------------------------------
# We can further split things up by coloring each point based on the content of the `Origin`
# column:
# ------------------------------------------------------------------------------------------

cars |> @vlplot(:point, x=:Miles_per_Gallon, y=:Acceleration, color=:Origin)

# ------------------------------------------------------------------------------------------
# We can also facet the plot so that we generate one plot per `Origin` value. We now use the
# `Cylinders` column for the color of the points:
# ------------------------------------------------------------------------------------------

cars |> @vlplot(:point, x=:Miles_per_Gallon, y=:Acceleration, column=:Origin, color=:Cylinders)

# ------------------------------------------------------------------------------------------
# Note how Vega-Lite automatically picked a continous color scale for `Cylinders` because
# that column is a numeric column. If we instead want things drawn as a nominal legend, we
# can specify the type of the encoding with a shorthand syntax:
# ------------------------------------------------------------------------------------------

cars |> @vlplot(:point, x=:Miles_per_Gallon, y=:Acceleration, column=:Origin, color="Cylinders:n")

# ------------------------------------------------------------------------------------------
# Another type of encoding is an `ordinal` encoding:
# ------------------------------------------------------------------------------------------

cars |> @vlplot(:point, x=:Miles_per_Gallon, y=:Acceleration, column=:Origin, color="Cylinders:o")

# ------------------------------------------------------------------------------------------
# We could also use different shapes of points for each `Origin`:
# ------------------------------------------------------------------------------------------

cars |> @vlplot(:point, x=:Miles_per_Gallon, y=:Acceleration, shape=:Origin, color="Cylinders:o")

# ------------------------------------------------------------------------------------------
# The plots are a little small, so lets increase the size:
# ------------------------------------------------------------------------------------------

cars |> @vlplot(:point, x=:Miles_per_Gallon, y=:Acceleration, shape=:Origin, color="Cylinders:o", height=300, width=300)

# ------------------------------------------------------------------------------------------
# Another type of plot is a histogram:
# ------------------------------------------------------------------------------------------

cars |> @vlplot(:bar, x={:Miles_per_Gallon, bin=true}, y="count()")

# ------------------------------------------------------------------------------------------
# We can easily split this up per `Origin` with the `column` channel:
# ------------------------------------------------------------------------------------------

cars |> @vlplot(:bar, x={:Miles_per_Gallon, bin=true}, y="count()", column=:Origin)

# ------------------------------------------------------------------------------------------
# Using Query.jl, we can manipulate all sorts of information into a table and then plot it:
# ------------------------------------------------------------------------------------------

readdir() |> @groupby(length(_)) |> @map({length=_.key, count=length(_)}) |> @vlplot(:bar, x=:count, y="length:o")

# ------------------------------------------------------------------------------------------
# Lets assign a plot to the variable `p`
# ------------------------------------------------------------------------------------------

p = cars |> @vlplot(:bar, x={:Miles_per_Gallon, bin=true}, y="count()", column=:Origin)

# ------------------------------------------------------------------------------------------
# And look at the type of `p`:
# ------------------------------------------------------------------------------------------

typeof(p)

# ------------------------------------------------------------------------------------------
# Plots can be saved in lots of different image formats:
# ------------------------------------------------------------------------------------------

p |> save("foo.png")

p |> save("foo.pdf")

p |> save("foo.svg")

p |> save("foo.eps")

# ------------------------------------------------------------------------------------------
# __WARNING__: The following line only works with `master` of FileIO.jl right now. Hopefully
# a new version will be tagged in a few days.
#
# We can also save a plot as a `.vegalite` file, the native Vega-Lite JSON format:
# ------------------------------------------------------------------------------------------

p |> save("foo.vegalite")

# ------------------------------------------------------------------------------------------
# Voyager uses Vega-Lite under the hood, so the two packages are tightly integrated. Lets
# pipe the `cars` dataset into Voyager again and create a plot in Voyager:
# ------------------------------------------------------------------------------------------

v = cars |> Voyager

# ------------------------------------------------------------------------------------------
# ![alt](screenshots/voyager3.png)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# We can now access the current plot in the Voyager window with the `[]` syntax:
# ------------------------------------------------------------------------------------------

v[]

# ------------------------------------------------------------------------------------------
# Note that this just returns a plot instance, as if we had created it with the `@vlplot`
# macro:
# ------------------------------------------------------------------------------------------

typeof(v[])

# ------------------------------------------------------------------------------------------
# We can therefore save it as files:
# ------------------------------------------------------------------------------------------

v[] |> save("foo2.pdf")

# ------------------------------------------------------------------------------------------
# __WARNING__: The following cells only works with the `master` branch of FileIO.jl.
#
# Here we can save the plot as a Vega-Lite file.
# ------------------------------------------------------------------------------------------

v[] |> save("foo.vegalite")

# ------------------------------------------------------------------------------------------
# We can then later load that file again, which returns a standard plot instance, as if we
# had created it with the `@vlplot` macro:
# ------------------------------------------------------------------------------------------

load("foo.vegalite")

# ------------------------------------------------------------------------------------------
# We can use this for interesting manipulations. For example, we can load a saved plot, but
# then pipe different data into it to create an updated plot:
# ------------------------------------------------------------------------------------------

cars |>
    @filter(_.Origin in ("USA", "Japan")) |>
    load("foo.vegalite")

# ------------------------------------------------------------------------------------------
# ## Missing values
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# The Queryverse uses [DataValues.jl](https://github.com/davidanthoff/DataValues.jl) for its
# missing value representation.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Conclusion
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ![alt](slides/Slide32.PNG)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ![alt](slides/Slide33.PNG)
# ------------------------------------------------------------------------------------------

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
function ω(xyz)
  θϕ = xyz2θϕ(xyz)
  θ,ϕ = θϕ
  cos(θ)*cos(ϕ)^2
end

function ωθϕ(θϕ)
  θ,ϕ = θϕ
  cos(θ)*cos(ϕ)^2
end

# spherical surface gradient computed analytically
function gradω(xyz)
  θϕr = xyz2θϕr(xyz)
  θ,ϕ,r = θϕr
  dθ = -sin(θ)*cos(ϕ)
  dϕ = -2*cos(θ)*cos(ϕ)*sin(ϕ)
  spherical_to_cartesian_matrix(θϕr)⋅VectorValue(dθ,dϕ,0)
end

function compute_error_weak_grad(model,order,degree,ls=BackslashSolver())
  # Setup geometry
  Ω=Triangulation(model)
  dΩ=Measure(Ω,degree)
  dω=Measure(Ω,degree,ReferenceDomain())

  # Setup H(div) spaces
  reffe_rt = ReferenceFE(raviart_thomas, Float64, order)
  V = FESpace(model, reffe_rt ; conformity=:HDiv)
  U = TrialFESpace(V)

  # L2 spaces
  reffe_lgn = ReferenceFE(lagrangian, Float64, order)
  S = FESpace(model, reffe_lgn; conformity=:L2)
  R = TrialFESpace(S)

  # Project the analytical field ω onto the L2 space
  a1(u,v) = ∫(v*u)dΩ
  b1(v)   = ∫(v*ω)dΩ
  op      = AffineFEOperator(a1,b1,S,R)
  fels    = LinearFESolver(ls)
  wh      = solve(fels,op)

  # Compute weak grad of wh
  a2(u,v) = ∫(v⋅u)dΩ
  b2(v)   = ∫((-wh)*DIV(v))dω
  op      = AffineFEOperator(a2,b2,U,V)
  fels    = LinearFESolver(ls)
  gradwh  = solve(fels,op)

  e    = gradwh-gradω
  dive = divergence(gradwh)-laplacian_unit_sphere(ωθϕ)∘xyz2θϕ
  sqrt(sum(∫(e⋅e+dive*dive)dΩ))
end

ext: jl

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
# ------------------------------------------------------------------------------------------
# ## Neural networks
#
# Now that we know what neurons are, we are ready for the final step: the neural network!. A
# neural network is literally made out of a network of neurons that are connected together.
#
# So far, we have just looked at single neurons, that only have a single output.
# What if we want multiple outputs?
#
#
# ### Multiple output models
#
# What if we wanted to distinguish between apples, bananas, *and* grapes? We could use
# *vectors* of `0` or `1` values to symbolize each output.
#
#
#
# The idea of using vectors is that different directions in the space of outputs encode
# information about different types of inputs.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Now we extend our previous model to give multiple outputs by repeating it with different
# weights. For the first element of the array we'd use:
#
# $$\sigma(x;w^{(1)},b^{(1)}) := \frac{1}{1 + \exp(-w^{(1)} \cdot x + b^{(1)})};$$
#
# then for the second we'd use
#
# $$\sigma(x;w^{(2)},b^{(2)}) := \frac{1}{1 + \exp(-w^{(2)} \cdot x + b^{(2)})};$$
#
# and if you wanted $n$ outputs, you'd have for each one
#
# $$\sigma(x;w^{(i)},b^{(i)}) := \frac{1}{1 + \exp(-w^{(i)} \cdot x + b^{(i)})}.$$
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Notice that these equations are all the same, except for the parameters, so we can write
# this model more succinctly, as follows. Let's write $b$ in an array:
#
# $$b=\left[\begin{array}{c}
# b_{1}\\
# b_{2}\\
# \vdots\\
# b_{n}
# \end{array}\right]$$
#
# and put our array of weights as a matrix:
#
# $$ \mathsf{W}=\left[\begin{array}{c}
# \\
# \\
# \\
# \\
# \end{array}\begin{array}{cccc}
# w_{1}^{(1)} & w_{2}^{(1)} & \ldots & w_{n}^{(1)}\\
# w_{1}^{(2)} & w_{2}^{(2)} & \ldots & w_{n}^{(2)}\\
# \vdots & \vdots &  & \vdots\\
# w_{1}^{(n)} & w_{2}^{(n)} & \ldots & w_{n}^{(n)}
# \end{array}\right]
# $$
#
# We can write this all in one line as:
#
# $$\sigma(x;w,b)= \left[\begin{array}{c}
# \sigma^{(1)}\\
# \sigma^{(2)}\\
# \vdots\\
# \sigma^{(n)}
# \end{array}\right] = \frac{1}{1 + \exp(-\mathsf{W} x + b)}$$
#
# $\mathsf{W} x$ is the operation called "matrix multiplication"
#
# [Show small matrix multiplication]
#
# It takes each column of weights and does the dot product against $x$ (remember, that's how
# $\sigma^{(i)}$ was defined) and spits out a vector from doing that with each column. The
# result is a vector, which makes this version of the function give a vector of outputs
# which we can use to encode larger set of choices.
#
# Matrix multiplication is also interesting since **GPUs (Graphics Processing Units, i.e.
# graphics cards) are basically just matrix multiplication machines**, which means that by
# writing the equation this way, the result can be calculated really fast.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# This "multiple input and multiple output" version of the sigmoid function is known as a
# *layer of neurons*.
#
# Previously we worked with a single neuron, which we visualized as
#
#
#
# where we have two pieces of data (green) coming into a single neuron (pink) that returned
# a single output. We could use this single output to do binary classification - to identify
# an image of a fruit as `1`, meaning banana or as `0`, meaning not a banana (or an apple).
#
# To do non-binary classification, we can use a layer of neurons, which we can visualize as
#
#
#
# We now have stacked a bunch of neurons on top of each other to hopefully work together and
# train to output results of more complicated features.
#
# We still have two input pieces of data, but now have several neurons, each of which
# produces an output for a given binary classification:
# * neuron 1: "is it an apple?"
# * neuron 2: "is it a banana?"
# * neuron 3: "is it a grape?"
# ------------------------------------------------------------------------------------------

max_stars_repo_path: .nbexports/introductory-tutorials/broader-topics-and-ecosystem/intro-to-ml/15. ML - Intro to neural networks.jl


text:
########################################################
# Functions related to fitting case series to SIRX model
########################################################

import Base.summary
using DataFrames

struct CaseModelFitSIRX <: CaseModelFitResult
    κ::Float64
    κ0::Float64
    IXRatio::Float64
    inputs::NamedTuple
    fit::LsqFit.LsqFitResult
end

"""
fitCaseModel(nt::Int64,cases::Vector{Float64},
    N::Float64,R0Free::Float64,TInfected::Float64,
    p0::Vector{Float64})

- nt: Integer with length of series
- cases: Vector of case number series
- N: Population size
- R0Free: User supplied R0 in freely spreading population
- TInfected: Average time of infection
- p0: Floating point vector with initial values of estimates
"""
function fitCaseModel(nt::Int64,cases::Vector{Float64},
    N::Float64,R0Free::Float64,TInfected::Float64,
    p0::Vector{Float64})
    inputs = (nt=nt,C0=cases[1],N=N,R0Free=R0Free,TInfected=TInfected,C=cases)
    model(t,p) = logCaseModelUnknown(nt,exp.(p),N,cases[1],R0Free,TInfected)
    fit = curve_fit(model,(1:nt)*1.0,log.(cases),log.(p0))
    return CaseModelFitSIRX( exp(fit.param[1]), exp(fit.param[2]),
    exp(fit.param[3]), inputs, fit )
end

#######################################################
function fitCaseModel_ps(nt::Int64,cases::Vector{Float64},
    N::Float64,R0Free::Float64,TInfected::Float64,p0::Vector{Float64})
    inputs = (nt=nt,C0=cases[1],N=N,R0Free=R0Free,TInfected=TInfected,C=cases)
    logcases = log.(cases)
    f(p) = mse(logcases,logCaseModelUnknown(nt,exp.(p),N,cases[1],
               R0Free,TInfected),nt-3.0)
    fit = optimize(f,log.(p0),ParticleSwarm(lower=-3.0.*[1.0,1.0,0.0],
                   upper=3.0.*[0.0,0.0,1.0],n_particles=10),
                   Optim.Options(iterations=10000))
    # fit = optimize(f,log.(p0),NelderMead())
#    fit = optimize(f,log.(p0),-# 3.0.*[1.0,1.0,1.0],3.0.*[1.0,1.0,1.0],SAMIN())
    param = Optim.minimizer(fit)

    return ( κ=exp(param[1]), κ0= exp(param[2]),
             IXRatio=exp(param[3]), inputs=inputs, fit=fit )
end

function fitCaseModel_nm(nt::Int64,cases::Vector{Float64},
    N::Float64,R0Free::Float64,TInfected::Float64,p0::Vector{Float64})
    inputs = (nt=nt,C0=cases[1],N=N,R0Free=R0Free,TInfected=TInfected,C=cases)
    logcases = log.(cases)
    f(p) = mse(logcases,logCaseModelUnknown(nt,exp.(p),N,cases[1],
               R0Free,TInfected),nt-3.0)
    fit = optimize(f,log.(p0),NelderMead())

    param = Optim.minimizer(fit)

    return ( κ=exp(param[1]), κ0= exp(param[2]),
             IXRatio=exp(param[3]), inputs=inputs, fit=fit )
end

function mse(y::Vector{Float64},yhat::Vector{Float64},df::Float64)
    return sum( (y-yhat).^2  )/df
end

#######################################################

function summary(fit::CaseModelFitSIRX)
    r0eff = R0Eff(getParams(fit.κ,fit.κ0,fit.inputs.R0Free,
            fit.inputs.TInfected,SIRX()))
    df = exp.(DataFrame(LsqFit.confidence_interval(fit.fit)))
    rename!(df,[:lo,:hi])
    param = [ "κ", "κ0", "IXRatio" ]
    est = [ fit.κ, fit.κ0, fit.IXRatio ]
    insertcols!(df,1,parameter=param)
    insertcols!(df,3,est=est)
    return (coef=df,R0Eff=r0eff,mse=LsqFit.mse(fit.fit))
end

function fitted(fit::CaseModelFitSIRX,nt::Int64=0)
    if (nt==0)
        nt = fit.inputs.nt
    end
    return exp.(logCaseModelUnknown(nt,
                exp.(fit.fit.param), fit.inputs.N, fit.inputs.C0,
                fit.inputs.R0Free,fit.inputs.TInfected))
end

function estimatedStates(fit::CaseModelFitSIRX,nt::Int64=0)
    if (nt==0)
        nt = fit.inputs.nt
    end
    d = getParams(exp(fit.fit.param[1]),
                  exp(fit.fit.param[2]), fit.inputs.R0Free,
                  fit.inputs.TInfected, SIRX() )
    return estimatedStates( nt, fit.inputs.N, fit.inputs.C0,
           exp(fit.fit.param[3]), d )
end


text:
export
    solve,
    low_energy_spectrum

function low_energy_spectrum(
    ig::IsingGraph,
    Dcut::Int,
    var_ϵ::T,
    max_sweeps::Int,
    dβ::T,
    β::T,
    schedule::Symbol,
    num_states::Int
) where T <: Number

    igp = prune(ig)
    ψ = MPS(igp, Dcut, var_ϵ, max_sweeps, dβ, β, schedule)
    states, probs, ldp = solve(ψ, num_states)

    en = energy.(states, Ref(igp))
    idx = sortperm(en)
    Solution(en[idx], states[idx], probs[idx], [0], ldp)
end

function solve(ψ::AbstractMPS, keep::Int)
    @assert keep > 0 "Number of states has to be > 0"
    T = eltype(ψ)

    keep_extra = keep
    lpCut = -Inf
    k = 1

    if keep < prod(rank(ψ)) keep_extra += 1 end

    lprob = zeros(T, k)
    states = fill([], 1, k)
    left_env = ones(eltype(ψ), 1, k)

    for (i, M) ∈ enumerate(ψ)
        _, d, b = size(M)

        pdo = ones(T, k, d)
        lpdo = zeros(T, k, d)
        LL = zeros(T, b, k, d)
        config = zeros(Int, i, k, d)

        for j ∈ 1:k
            L = left_env[:, j]

            for σ ∈ local_basis(d)
                m = idx(σ)
                LL[:, j, m] = L' * M[:, m, :]
                pdo[j, m] = dot(LL[:, j, m], LL[:, j, m])
                config[:, j, m] = vcat(states[:, j]..., σ)
                LL[:, j, m] = LL[:, j, m] / sqrt(pdo[j, m])
            end
            pdo[j, :] = pdo[j, :] / sum(pdo[j, :])
            lpdo[j, :] = log.(pdo[j, :]) .+ lprob[j]
        end

        perm = collect(1 : k * d)
        k = k * d

        if k > keep_extra
            k = keep_extra
            partialsortperm!(perm, vec(lpdo), 1:k, rev=true)
            lprob = vec(lpdo)[perm]
            lpCut < last(lprob) ? lpCut = last(lprob) : ()
        end

        lprob = vec(lpdo)[perm]
        @cast A[α, (l, d)] |= LL[α, l, d]
        left_env = A[:, perm]
        @cast B[β, (l, d)] |= config[β, l, d]
        states = B[:, perm]
    end
    Vector.(eachcol(states[:, 1:keep])), lprob[1:keep], lpCut
end

function _apply_bias!(
    ψ::AbstractMPS,
    ig::LabelledGraph,
    dβ::Number,
    i::Int
)
    M = ψ[i]
    h = get_prop(ig, i, :h)
    v = exp.(-0.5 * dβ * h * local_basis(ψ, i))
    @cast M[x, σ, y] = M[x, σ, y] * v[σ]
    ψ[i] = M
end

function _apply_exponent!(
    ψ::AbstractMPS,
    ig::LabelledGraph,
    dβ::Number,
    i::Int,
    j::Int,
    last::Int
)
    M = ψ[j]
    D = typeof(M).name.wrapper(I(physical_dim(ψ, i)))

    J = get_prop(ig, i, j, :J)
    C = exp.(-0.5 * dβ * local_basis(ψ, i) * J * local_basis(ψ, j)')

    if j == last
        @cast M̃[(x, a), σ, b] := C[x, σ] * M[a, σ, b]
    else
        @cast M̃[(x, a), σ, (y, b)] := C[x, σ] * D[x, y] * M[a, σ, b]
    end
    ψ[j] = M̃
end

function _apply_projector!(ψ::AbstractMPS, i::Int)
    M = ψ[i]
    D = typeof(M).name.wrapper(I(physical_dim(ψ, i)))
    @cast M̃[a, σ, (y, b)] := D[σ, y] * M[a, σ, b]
    ψ[i] = M̃
end

function _apply_nothing!(ψ::AbstractMPS, l::Int, i::Int)
    M = ψ[l]
    D = typeof(M).name.wrapper(I(physical_dim(ψ, i)))
    @cast M̃[(x, a), σ, (y, b)] := D[x, y] * M[a, σ, b]
    ψ[l] = M̃
end

function purifications(χ::T, ϕ::T) where {T <: AbstractMPS}
    S = promote_type(eltype(χ), eltype(ϕ))
    ψ = MPS(S, length(ϕ))
    for (i, (A, B)) ∈ enumerate(zip(χ, ϕ))
        @cast C[(l, l̃), σ, (r, r̃)] := A[l, σ, r] * B[l̃, σ, r̃]
        ψ[i] = C
    end
    ψ
end
purifications(χ) = purifications(χ, χ)

_holes(l::Int, nbrs::Vector) = setdiff(l+1:last(nbrs), nbrs)

function _apply_gates(
    ρ::AbstractMPS,
    ig::IsingGraph,
    Dcut::Int,
    var_ϵ::Number,
    sweeps::Int,
    dβ::Number
)
    is_right = true
    for i ∈ 1:nv(ig)
        _apply_bias!(ρ, ig, dβ, i)
        is_right = false

        nbrs = unique_neighbors(ig, i)
        if !isempty(nbrs)
            _apply_projector!(ρ, i)
            for j ∈ nbrs
                _apply_exponent!(ρ, ig, dβ, i, j, last(nbrs))
            end
            for l ∈ _holes(i, nbrs)
                _apply_nothing!(ρ, l, i)
            end
        end

        if bond_dimension(ρ) > Dcut
            ρ = SpinGlassTensors.compress(ρ, Dcut, var_ϵ, sweeps)
            is_right = true
        end
    end

    if !is_right SpinGlassTensors.canonise!(ρ, :right) end
    ρ
end

function SpinGlassTensors.MPS(
    ig::IsingGraph,
    Dcut::Int,
    var_ϵ::Number,
    sweeps::Int,
    schedule::Vector{<:Number}
)
    ρ = HadamardMPS(ig)
    for dβ ∈ schedule
        ρ =_apply_gates(ρ, ig, Dcut, var_ϵ, sweeps, dβ)
    end
    ρ
end

#TODO: WTF?
function SpinGlassTensors.MPS(
    ig::IsingGraph,
    Dcut::Int,
    var_ϵ::Number,
    sweeps::Int,
    dβ::Number,
    β::Number,
    schedule::Symbol
)
    ρ = HadamardMPS(ig)
    is_right = true

    if schedule == :log
        k = Int(ceil(log2(β/dβ)))
        ρ = _apply_gates(ρ, ig, Dcut, var_ϵ, sweeps, β/(2^k))
        for _ ∈ 1:k
            ρ = purifications(ρ)
            if bond_dimension(ρ) > Dcut
                ρ = compress(ρ, Dcut, var_ϵ, sweeps)
                is_right = true
            end
        end
    elseif schedule == :lin
        k = ceil(β/dβ)
        ρ = _apply_gates(ρ, ig, Dcut, var_ϵ, sweeps, β/k)
        ρ0 = copy(ρ)
        for _ ∈ 1:Int(k)
            ρ = purifications(ρ, ρ0)
            if bond_dimension(ρ) > Dcut
                ρ = compress(ρ, Dcut, var_ϵ, sweeps)
                is_right = true
            end
        end
    end
    if !is_right SpinGlassTensors.canonise!(ρ, :right) end
    ρ
end

function HadamardMPS(::Type{T}, ig::IsingGraph) where T <: Number
    MPS([
        fill(one(T), r) ./ sqrt(T(r))
        for r ∈ values(get_prop(ig, :rank))
    ])
end
HadamardMPS(ig) = HadamardMPS(Float64, ig)


text:
@params mutable struct BESOResult{T}
    topology::AbstractVector{T}
    objval::T
    change::T
    converged::Bool
    fevals::Int
end

"""
The BESO algorithm, see [HuangXie2010](@cite).
"""
@params struct BESO{T} <: TopOptAlgorithm
    comp::Compliance
    vol::Volume
    vol_limit
    filter
    vars::Vector{T}
    topology::Vector{T}
    er::T
    maxiter::Int
    p::T
    sens::AbstractVector{T}
    old_sens::AbstractVector{T}
    obj_trace::MVector{<:Any, T}
    tol::T
    sens_tol::T
    result::BESOResult{T}
end
Base.show(::IO, ::MIME{Symbol("text/plain")}, ::BESO) = println("TopOpt BESO algorithm")

function BESO(comp::Compliance, vol::Volume, vol_limit, filter; maxiter = 200, tol = 0.0001, p = 3.0, er = 0.02, sens_tol = tol/100, k = 10)
    solver = comp.solver
    T = eltype(solver.vars)
    solver = comp.solver
    topology = zeros(T, getncells(solver.problem.ch.dh.grid))
    result = BESOResult(topology, T(NaN), T(NaN), false, 0)
    black = solver.problem.black
    white = solver.problem.white
    nvars = length(topology) - sum(black) - sum(white)
    vars = zeros(T, nvars)
    sens = zeros(T, nvars)
    old_sens = zeros(T, nvars)
    obj_trace = zeros(MVector{k, T})

    return BESO(comp, vol, vol_limit, filter, vars, topology, er, maxiter, p, sens, old_sens, obj_trace, tol, sens_tol, result)
end

update_penalty!(b::BESO, p::Number) = (b.p = p)

function (b::BESO)(x0 = copy(b.obj.solver.vars))
    T = eltype(x0)
    @unpack sens, old_sens, er, tol, maxiter = b
    @unpack obj_trace, topology, sens_tol, vars = b
    @unpack solver = b.comp
    @unpack varind, black, white = solver.problem
    @unpack total_volume, cellvolumes = b.vol
    V = b.vol_limit
    k = length(obj_trace)

    # Initialize the topology
    for i in 1:length(topology)
        if black[i]
            topology[i] = 1
        elseif white[i]
            topology[i] = 0
        else
            topology[i] = round(x0[varind[i]])
            vars[varind[i]] = topology[i]
        end
    end

    # Calculate the current volume fraction
    true_vol = vol = dot(topology, cellvolumes) / total_volume
    # Main loop
    change = T(1)
    iter = 0
    setpenalty!(solver, b.p)
    while (change > tol || true_vol > V) && iter < maxiter
        iter += 1
        if iter > 1
            old_sens .= sens
        end
        vol = max(vol*(1-er), V)
        for j in max(2, k-iter+2):k
            obj_trace[j-1] = obj_trace[j]
        end
        obj_trace[k], pb = Zygote.pullback(x -> b.comp(b.filter(x)), vars)
        sens = pb(1.0)[1]
        rmul!(sens, -1)
        if iter > 1
            @. sens = (sens + old_sens) / 2
        end
        l1, l2 = minimum(sens), maximum(sens)
        while (l2 - l1) / l2 > sens_tol
            th = (l1 + l2) / 2
            for i in 1:length(topology)
                if !black[i] && !white[i]
                    topology[i] = T(sign(sens[varind[i]] - th) > 0)
                    vars[varind[i]] = topology[i]
                end
            end
            if dot(topology, cellvolumes) - vol * total_volume > 0
                l1 = th
            else
                l2 = th
            end
        end
        true_vol = dot(topology, cellvolumes) / total_volume
        if iter >= k
            l = sum(@view obj_trace[1:k÷2])
            h = sum(@view obj_trace[k÷2+1:k])
            change = abs(l-h)/h
        end
    end

    objval = obj_trace[k]
    converged = change <= tol
    fevals = iter
    @pack! b.result = change, objval, converged, fevals

    return b.result
end


text:
using LinearAlgebra
using StaticArrays
using ConstrainedDynamics
using ConstrainedDynamics: State

function simplependulum2D(steps::Int; Δt::Real = 0.01, θstart::Real = 0., ωstart::Real = 0., Δm::Real = 1.0, ΔJ::Real = 1.0, friction::Real = 0, threadlock = nothing)
    joint_axis = [1.0; 0.0; 0.0]
    m = 1.0Δm
    g = -9.81
    l = 1.0
    r = 0.01
    p2 = [0.0;0.0;l / 2] # joint connection point
    ΔT = Δt * steps
    if threadlock === nothing
        mech, origin, link1 = _simplependulum2Dmech(r, l, m, ΔJ, joint_axis, p2, g, Δt)
    else
        lock(threadlock)  # ConstrainedDynamics uses globals with possible data races to initialize mechanisms -> lock
        try
            mech, origin, link1 = _simplependulum2Dmech(r, l, m, ΔJ, joint_axis, p2, g, Δt)
        finally
            unlock(threadlock)
        end
    end
    q1 = UnitQuaternion(RotX(θstart))
    setPosition!(origin, link1; p2=p2, Δq=q1)
    setVelocity!(origin, link1; p2=p2, Δω=SA[ωstart, 0., 0.])
    initialstates = [deepcopy(body.state) for body in mech.bodies]
    if friction != 0
        function control!(mech, _)
            setForce!(mech.bodies[1], τ=-SA[1.,0,0]friction*mech.bodies[1].state.ωc[1])
        end
        storage = simulate!(mech, ΔT, control!, record = true)
    else
        storage = simulate!(mech, ΔT, record = true)
    end
    return storage, mech, initialstates
end

function _simplependulum2Dmech(r, l, m, ΔJ, joint_axis, p2, g, Δt)
    origin = Origin{Float64}()
    link1 = Cylinder(r, l, m)
    link1.J = I*1/12*m*l^2*ΔJ + zeros(3,3)  # Zeros to set broadcast shape
    joint_between_origin_and_link1 = EqualityConstraint(Revolute(origin, link1, joint_axis; p2=p2))
    links = [link1]
    constraints = [joint_between_origin_and_link1]
    return Mechanism(origin, links, constraints, g=g, Δt=Δt), origin, link1
end

function doublependulum2D(steps::Int; Δt::Real = 0.01, θstart::Vector{<:Real} = [0., 0.], ωstart::Vector{<:Real} = [0., 0.], Δm::Vector{<:Real} = ones(2),
                          ΔJ::Vector{<:Real} = ones(2), friction::Vector{<:Real} = zeros(2), threadlock = nothing)
    # Parameters
    l1 = 1.0
    l2 = 1.0
    m = ones(2) .* Δm
    x, y = .1, .1
    vert11 = [0.;0.;l1 / 2]
    vert12 = -vert11
    vert21 = [0.;0.;l2 / 2]
    joint_axis = [1.0; 0.0; 0.0]
    ΔT = Δt * steps
    # Initial orientation
    q1 = UnitQuaternion(RotX(θstart[1]))
    q2 = UnitQuaternion(RotX(θstart[2]))
    if threadlock === nothing
        mech, origin, link1, link2 = _doublependulum2Dmech(x, y, l1, l2, m, ΔJ, joint_axis, vert11, vert12, vert21, Δt)
    else
        lock(threadlock)
        try
            mech, origin, link1, link2 = _doublependulum2Dmech(x, y, l1, l2, m, ΔJ, joint_axis, vert11, vert12, vert21, Δt)
        finally
            unlock(threadlock)
        end
    end
    setPosition!(origin,link1, p2=vert11, Δq=q1)
    setPosition!(link1, link2, p1=vert12, p2=vert21, Δq=q2)
    setVelocity!(origin, link1; p2=vert11, Δω=SA[ωstart[1], 0., 0.])
    setVelocity!(link1, link2; p1=vert12, p2=vert21, Δω=SA[ωstart[2], 0., 0.])
    initialstates = [deepcopy(body.state) for body in mech.bodies]
    if (any(friction .!= 0))
        function control!(mech, _)
            setForce!(mech.bodies[1], τ=-SA[1.,0,0]friction[1]*mech.bodies[1].state.ωc[1])
            setForce!(mech.bodies[2], τ=-SA[1.,0,0]friction[2]*mech.bodies[2].state.ωc[1])
        end
        storage = simulate!(mech, ΔT, control!, record = true)
    else
        storage = simulate!(mech, ΔT, record = true)
    end
    return storage, mech, initialstates
end

function _doublependulum2Dmech(x, y, l1, l2, m, ΔJ, joint_axis, vert11, vert12, vert21, Δt)
    origin = Origin{Float64}()
    link1 = Box(x, y, l1, m[1], color = RGBA(1., 1., 0.))
    link2 = Box(x, y, l2, m[2], color = RGBA(1., 1., 0.))
    link1.J = I*1/12*m[1]*l1^2*ΔJ[1] + zeros(3,3)  # Zeros to set broadcast shape
    link2.J = I*1/12*m[2]*l2^2*ΔJ[2] + zeros(3,3)
    socket0to1 = EqualityConstraint(Revolute(origin, link1, joint_axis; p2=vert11))
    socket1to2 = EqualityConstraint(Revolute(link1, link2, joint_axis; p1=vert12, p2=vert21))
    links = [link1;link2]
    constraints = [socket0to1;socket1to2]
    mech = Mechanism(origin, links, constraints, Δt=Δt)
    return mech, origin, link1, link2
end

function cartpole(steps::Int; Δt::Real = 0.01, xstart::Real = 0., θstart::Real = 0., vstart::Real = 0., ωstart::Real = 0.,
                  Δm::Vector{<:Real} = ones(2), ΔJ::Real = 1., friction::Vector{<:Real} = zeros(2), threadlock = nothing)
    xaxis = [1.0; 0.0; 0.0]
    yaxis = [0.0; 1.0; 0.0]
    g = -9.81
    l = 0.5
    r = 0.01
    m = ones(2) .* Δm
    p01 = [0.0; 0.0; 0.0] # joint connection point
    p12 = [0.0; 0.0; 0.0]
    p21 = [0.0; 0.0; l/2]
    ΔT = steps*Δt
    x, y, z = 0.2, 0.3, 0.1
    if threadlock === nothing
        mech, origin, link1, link2 = _cartpolemech(x, y, z, r, l, m, ΔJ, xaxis, yaxis, p01, p12, p21, Δt)
    else
        lock(threadlock)
        try
            mech, origin, link1, link2 = _cartpolemech(x, y, z, r, l, m, ΔJ, xaxis, yaxis, p01, p12, p21, Δt)
        finally
            unlock(threadlock)
        end
    end
    setPosition!(origin, link1; p2=p01, Δx=SA[0., xstart, 0.])
    setPosition!(link1, link2; p1=p12, p2=p21, Δq=UnitQuaternion(RotX(θstart)))
    setVelocity!(origin, link1; p2=p01, Δv=SA[0., vstart, 0.])
    setVelocity!(link1, link2; p1=p12, p2=p21, Δω=SA[ωstart, 0., 0.])
    initialstates = [deepcopy(body.state) for body in mech.bodies]
    if (any(friction .!= 0))
        function control!(mech, _)
            setForce!(mech.bodies[1], F=-SA[0,1.,0]friction[1]*mech.bodies[1].state.vc[2])
            setForce!(mech.bodies[2], τ=-SA[1.,0,0]friction[2]*mech.bodies[2].state.ωc[1])
        end
        storage = simulate!(mech, ΔT, control!, record = true)
    else
        storage = simulate!(mech, ΔT, record = true)
    end
    return storage, mech, initialstates
end

function _cartpolemech(x, y, z, r, l, m, ΔJ, xaxis, yaxis, p01, p12, p21, Δt)
    origin = Origin{Float64}()
    link1 = Box(x, y, z, m[1])
    link2 = Cylinder(r, l, m[2])
    # link1 inertia doesn't matter, no rotation
    link2.J = I*1/12*m[2]*l^2*ΔJ + zeros(3,3)  # Only xx matters
    joint_origin_link1 = EqualityConstraint(Prismatic(origin, link1, yaxis; p2=p01))
    joint_link1_link2 = EqualityConstraint(Revolute(link1, link2, xaxis; p1=p12, p2=p21))
    links = [link1, link2]
    constraints = [joint_origin_link1, joint_link1_link2]
    mech = Mechanism(origin, links, constraints, Δt=Δt)
    return mech, origin, link1, link2
end

function fourbar(steps::Int; Δt::Real = 0.01, θstart::Vector{<:Real} = [0., 0.], Δm::Vector{<:Real} = ones(4), ΔJ::Vector{<:Real} = ones(4),
                 friction::Vector{<:Real} = zeros(4), threadlock = nothing)
    # Parameters
    ex = [1.;0.;0.]
    l = 1.0
    m = ones(4) .* Δm
    x, y = .1, .1
    vert11 = [0.;0.;l / 2]
    vert12 = -vert11
    ΔT = steps*Δt
    # Initial orientation
    q = UnitQuaternion(RotX(θstart[2]))
    qoff = UnitQuaternion(RotX(θstart[1]))
    Δq1 = q * qoff
    Δq2 = q

    if threadlock === nothing
        mech, origin, links = _fourbarmech(x, y, l, m, ΔJ, ex, vert11, vert12, Δt)
    else
        lock(threadlock)
        try
            mech, origin, links = _fourbarmech(x, y, l, m, ΔJ, ex, vert11, vert12, Δt)
        finally
            unlock(threadlock)
        end
    end
    setPosition!(origin, links[1], p2 = vert11, Δq = Δq1)
    setPosition!(links[1], links[2], p1 = vert12, p2 = vert11, Δq = inv(Δq2) * inv(Δq2))
    setPosition!(links[1], links[3], p1 = vert11, p2 = vert11, Δq = inv(Δq2) * inv(Δq2))
    setPosition!(links[3], links[4], p1 = vert12, p2 = vert11, Δq = Δq2 * Δq2)
    initialstates = [deepcopy(body.state) for body in mech.bodies]
    if (any(friction .!= 0))
        control!(mech, _) = for i in 1:4 setForce!(mech.bodies[i], τ=-SA[1.,0,0]friction[i]*mech.bodies[i].state.ωc[1]) end
        storage = simulate!(mech, ΔT, control!, record = true)
    else
        storage = simulate!(mech, ΔT, record = true)
    end
    return storage, mech, initialstates
end

function _fourbarmech(x, y, l, m, ΔJ, ex, vert11, vert12, Δt)
    origin = Origin{Float64}()
    links = [Box(x, y, l, m[i], color = RGBA(1., 1., 0.)) for i in 1:4]
    for (i, link) in enumerate(links)
        link.J = I*1/12*m[i]*l^2*ΔJ[i] + zeros(3,3)  # Zeros to set broadcast shape
    end
    j1 = EqualityConstraint(Revolute(origin, links[1], ex; p2=vert11))
    j2 = EqualityConstraint(Revolute(links[1], links[2], ex; p1=vert12, p2=vert11), Cylindrical(links[1], links[3], ex; p1=vert11, p2=vert11))
    j3 = EqualityConstraint(Revolute(links[3], links[4], ex; p1=vert12, p2=vert11))
    j4 = EqualityConstraint(Revolute(links[2], links[4], ex; p1=vert12, p2=vert12))
    constraints = [j1, j2, j3, j4]
    mech = Mechanism(origin, links, constraints, Δt=Δt)
    return mech, origin, links
end


text:
# Scaling matrix

immutable ScalMat{T<:Real} <: AbstractPDMat{T}
  dim::Int
  value::T
  inv_value::T
end

ScalMat(d::Int,v::Real) = ScalMat{typeof(one(v)/v)}(d, v, one(v) / v)

### Conversion
# can rewrite convert(T, a.value) as T(a.value) when we drop v0.3 support
convert{T<:Real}(::Type{ScalMat{T}}, a::ScalMat) = ScalMat(a.dim, convert(T, a.value))

### Basics

dim(a::ScalMat) = a.dim
full(a::ScalMat) = diagm(fill(a.value, a.dim))
diag(a::ScalMat) = fill(a.value, a.dim)

### Arithmetics

function pdadd!(r::Matrix, a::Matrix, b::ScalMat, c)
    @check_argdims size(r) == size(a) == size(b)
    if is(r, a)
        _adddiag!(r, b.value * c)
    else
        _adddiag!(copy!(r, a), b.value * c)
    end
    return r
end

*{T<:Real}(a::ScalMat, c::T) = ScalMat(a.dim, a.value * c)
/{T<:Real}(a::ScalMat{T}, c::T) = ScalMat(a.dim, a.value / c)
*(a::ScalMat, x::StridedVecOrMat) = a.value * x
\(a::ScalMat, x::StridedVecOrMat) = a.inv_value * x

### Algebra

inv(a::ScalMat) = ScalMat(a.dim, a.inv_value, a.value)
logdet(a::ScalMat) = a.dim * log(a.value)
eigmax(a::ScalMat) = a.value
eigmin(a::ScalMat) = a.value

### whiten and unwhiten

function whiten!(r::StridedVecOrMat, a::ScalMat, x::StridedVecOrMat)
    @check_argdims dim(a) == size(x, 1)
    c = sqrt(a.inv_value)
    for i = 1:length(x)
        @inbounds r[i] = x[i] * c
    end
    return r
end

function unwhiten!(r::StridedVecOrMat, a::ScalMat, x::StridedVecOrMat)
    @check_argdims dim(a) == size(x, 1)
    c = sqrt(a.value)
    for i = 1:length(x)
        @inbounds r[i] = x[i] * c
    end
    return r
end

### quadratic forms

quad(a::ScalMat, x::Vector) = sumabs2(x) * a.value
invquad(a::ScalMat, x::Vector) = sumabs2(x) * a.inv_value

quad!(r::AbstractArray, a::ScalMat, x::Matrix) = colwise_sumsq!(r, x, a.value)
invquad!(r::AbstractArray, a::ScalMat, x::Matrix) = colwise_sumsq!(r, x, a.inv_value)

### tri products

function X_A_Xt(a::ScalMat, x::StridedMatrix)
    @check_argdims dim(a) == size(x, 2)
    gemm('N', 'T', a.value, x, x)
end

function Xt_A_X(a::ScalMat, x::StridedMatrix)
    @check_argdims dim(a) == size(x, 1)
    gemm('T', 'N', a.value, x, x)
end

function X_invA_Xt(a::ScalMat, x::StridedMatrix)
    @check_argdims dim(a) == size(x, 2)
    gemm('N', 'T', a.inv_value, x, x)
end

function Xt_invA_X(a::ScalMat, x::StridedMatrix)
    @check_argdims dim(a) == size(x, 1)
    gemm('T', 'N', a.inv_value, x, x)
end


text:
```julia
using CatmullRom, Plots

xt(t) = 16 * sin(t)^3
yt(t) = 13 * cos(t) - 5 * cos(2*t) - 2 * cos(3*t) - cos(4*t)

ts = collect(range(0.0, length=15, stop=2.0*pi));
n_between = 28;

xs=xt.(ts); ys=yt.(ts);
points = collect(zip(xs,ys));

crpts = catmullrom(points, n_between);
cxs,cys=crpts;

plot(xs, ys, linecolor=:lightgreen, linewidth=7, size=(500,500), xaxis=nothing, yaxis=nothing, legend=nothing)
plot!(cxs, cys, linecolor=:black, linewidth=2, size=(500,500), xaxis=nothing, yaxis=nothing, legend=nothing)

arcpoints = (8, 64)
crpts = catmullrom_by_arclength(points, arcpoints);
cxs,cys=crpts;

plot(xs, ys, linecolor=:lightgreen, linewidth=7, size=(500,500), xaxis=nothing, yaxis=nothing, legend=nothing)
plot!(cxs, cys, linecolor=:black, linewidth=2, size=(500,500), xaxis=nothing, yaxis=nothing, legend=nothing)

```


text:
# This file was generated, do not modify it. # hide
fp = fitted_params(mach)
@show round.(fp.coefs[1:3], sigdigits=3)
@show round(fp.intercept, sigdigits=3)

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


text:
using SafeTestsets
using Test

@testset "JuliaProject.jl" begin
    @time @safetestset "Doubling" begin include("doubling.jl") end
end


text:
using LinearAlgebra

"""
    kiops(tstops, A, u; kwargs...) -> (w, stats)

Evaluate a linear combinaton of the ``φ`` functions evaluated at ``tA`` acting on
vectors from ``u``, that is

```math
  w(i) = φ_0(t[i] A) u[:, 1] + φ_1(t[i] A) u[:, 2] + φ_2(t[i] A) u[:, 3] + ...
```

The size of the Krylov subspace is changed dynamically during the integration.
The Krylov subspace is computed using the incomplete orthogonalization method.

Arguments:
  - `tstops`     - Array of `tstop`
  - `A`          - the matrix argument of the ``φ`` functions
  - `u`          - the matrix with columns representing the vectors to be multiplied by the ``φ`` functions

Keyword arguments:
  - `tol`        - the convergence tolerance required (default: 1e-7)
  - `mmin`, `mmax` - let the Krylov size vary between mmin and mmax (default: 10, 128)
  - `m`          - an estimate of the appropriate Krylov size (default: mmin)
  - `iop`        - length of incomplete orthogonalization procedure (default: 2)
  - `ishermitian` -  whether ``A`` is Hermitian (default: ishermitian(A))
  - `opnorm`     -  the operator norm of ``A`` (default: opnorm(A, Inf))
  - `task1`      - if true, divide the result by 1/T^p

Returns:
  - `w`        - the linear combination of the ``φ`` functions evaluated at ``tA`` acting on the vectors from ``u``
  - `stats[1]` - number of substeps
  - `stats[2]` - number of rejected steps
  - `stats[3]` - number of Krylov steps
  - `stats[4]` - number of matrix exponentials
  - `stats[5]` - the Krylov size of the last substep

`n` is the size of the original problem
`p` is the highest index of the ``φ`` functions

References:
* Gaudreault, S., Rainwater, G. and Tokman, M., 2018. KIOPS: A fast adaptive Krylov subspace solver for exponential integrators. Journal of Computational Physics. Based on the PHIPM and EXPMVP codes (http://www1.maths.leeds.ac.uk/~jitse/software.html). https://gitlab.com/stephane.gaudreault/kiops.
* Niesen, J. and Wright, W.M., 2011. A Krylov subspace method for option pricing. SSRN 1799124
* Niesen, J. and Wright, W.M., 2012. Algorithm 919: A Krylov subspace algorithm for evaluating the ``φ``-functions appearing in exponential integrators. ACM Transactions on Mathematical Software (TOMS), 38(3), p.22
"""
function kiops(tau_out, A, u; mmin::Int = 10, mmax::Int = 128, m::Int=min(mmin, mmax),
               tol::Real=1e-7, opnorm=LinearAlgebra.opnorm(A, Inf), iop::Int=2,
               ishermitian::Bool=LinearAlgebra.ishermitian(A), task1::Bool = false)
    n, ppo = size(u, 1), size(u, 2)
    p = ppo - 1

    if p == 0
        p = 1
        # Add extra column of zeros
        u = [u zero(u)]
    end

    # Preallocate matrix
    TA, Tb = eltype(A), eltype(u)
    T = promote_type(TA, Tb)
    Ks = KrylovSubspace{T, ishermitian ? real(T) : T}(n, m, p)

    step    = 0
    krystep = 0
    ireject = 0
    reject  = 0
    exps    = 0
    sgn     = sign(tau_out[end])
    tau_now = 0
    tau_end = abs(tau_out[end])
    j       = 0

    numSteps = size(tau_out, 2)

    # Initial condition
    w     = zeros(n, numSteps)
    w_aug = zeros(p)
    copyto!(@view(w[:, 1]), @view(u[:, 1]))

    # Normalization factors
    normU = norm(@view(u[:, 2:end]), 1)
    if ppo > 1 && normU > 0
        ex = ceil(log2(normU))
        nu = exp2(-ex)
        mu = exp2(ex)
    else
        nu = 1
        mu = 1
    end

    # Flip the rest of the u matrix
    u_flip = reverse(@view(u[:, 2:end]), dims = 2)
    rmul!(u_flip, nu)

    # Compute and initial starting approximation for the step size
    tau = tau_end

    # Setting the safety factors and tolerance requirements
    if tau_end > 1
        gamma = 0.2
        gamma_mmax = 0.1
    else
        gamma = 0.9
        gamma_mmax = 0.6
    end
    delta = 1.4

    # Used in the adaptive selection
    oldm = -1; oldtau = NaN; omega = NaN
    orderold = true; kestold = true

    l=1

    local beta, kest
    while tau_now < tau_end
        oldj = Ks.m
        arnoldi!(Ks, (A, u_flip), (w, w_aug); opnorm=opnorm, ishermitian=ishermitian, iop=iop, init=j, t=tau_now, mu=mu, l=l, m=m)
        V = getfield(Ks, :V)
        H = getfield(Ks, :H)
        j = Ks.m
        happy = j < oldj
        beta = Ks.beta

        # To obtain the phi_1 function which is needed for error estimate
        H[1, j + 1] = 1

        # Save h_j+1,j and remove it temporarily to compute the exponential of H
        nrm = H[j + 1, j]
        H[j + 1, j] = 0

        # Compute the exponential of the augmented matrix
        F = exp(sgn * tau * H[1:j + 1, 1:j + 1])
        exps = exps + 1

        # Restore the value of H_{m+1,m}
        H[j + 1, j] = nrm

        if happy
            # Happy breakdown wrap up
            omega   = 0
            tau_new = min(tau_end - (tau_now + tau), tau)
            m_new   = m
            happy   = false
        else
            # Local truncation error estimation
            err = abs(beta * nrm * F[j, j + 1])

            # Error for this step
            oldomega = omega
            omega = tau_end * err / (tau * tol)

            # Estimate order
            if m == oldm && tau != oldtau && ireject >= 1
                order = max(1, log(omega/oldomega) / log(tau/oldtau))
                orderold = false
            elseif orderold || ireject == 0
                orderold = true
                order = j/4
            else
                orderold = true
            end
            # Estimate k
            if m != oldm && tau == oldtau && ireject >= 1
                kest = max(1.1, (omega/oldomega)^(1/(oldm-m)))
                kestold = false
            elseif kestold || ireject == 0
                kestold = true
                kest = 2
            else
                kestold = true
            end

            if omega > delta
                remaining_time = tau_end - tau_now
            else
                remaining_time = tau_end - (tau_now + tau)
            end

            # Krylov adaptivity

            same_tau = min(remaining_time, tau)
            tau_opt  = tau * (gamma / omega)^(1 / order)
            tau_opt  = min(remaining_time, max(tau/5, min(5*tau, tau_opt)))

            m_opt = ceil(Int, j + log(omega / gamma) / log(kest))
            m_opt = max(mmin, min(mmax, max(3÷4*m, min(m_opt, cld(4, 3)*m))))

            if j == mmax
                if omega > delta
                    m_new = j
                    tau_new = tau * (gamma_mmax / omega)^(1 / order)
                    tau_new = min(tau_end - tau_now, max(tau/5, tau_new))
                else
                    tau_new = tau_opt
                    m_new = m
                end
            else
                m_new = m_opt
                tau_new = same_tau
            end
        end

        # Check error against target
        if omega <= delta
            tau_now, l, j, reject, ireject, step = kiops_update_solution!(tau_now, tau, tau_out, w, l, V, F, H, beta, j, n, step, numSteps, reject, ireject)
        else
            # Nope, try again
            ireject = ireject + 1

            # Restore the original matrix
            H[1, j + 1] = 0
        end

        oldtau = tau
        tau    = tau_new

        oldm = m
        m    = m_new
    end

    # FIXME: phiHan doesn't seem right
    if tau_out[1]!=1 && task1
        wl = @view w[:, l]
        if length(tau_out)==1
            @. wl = wl*(1/tau_out[l])^p
        else
            tmp = maximum(abs, u, dims=2)
            phiHan = map(first, findall(!iszero, tmp))

            if isempty(phiHan)
                phiHan = length(u)
            else
                phiHan = phiHan .- 1
            end

            for l = 1:numSteps
                @. wl = wl*(1/tau_out[l]^phiHan)
            end
        end
    end

    m_ret=m

    stats = (step, reject, krystep, exps, m_ret)

    return w, stats
end

Base.@propagate_inbounds function kiops_update_solution!(tau_now, tau, tau_out, w, l, V, F, H, beta, j, n, step, numSteps, reject, ireject)
    # Yep, got the required tolerance update
    reject = reject + ireject
    step = step + 1

    # Udate for tau_out in the interval (tau_now, tau_now + tau)
    blownTs = 0
    nextT = tau_now + tau
    for k = l:numSteps
        if abs(tau_out[k]) < abs(nextT)
            blownTs = blownTs + 1
        end
    end

    if blownTs != 0
        # Copy current w to w we continue with.
        copyto!(@view(w[:, l+blownTs]), @view(w[:, l]))

        for k = 0:blownTs - 1
            tauPhantom = tau_out[l+k] - tau_now
            F2 = exp(sign(tau_out[end]) * tauPhantom * @view(H[1:j, 1:j]))
            mul!(@view(w[:, l+k]), @view(V[1:n, 1:j]), @view(F2[1:j, 1]))
            rmul!(@view(w[:, l+k]), beta)
        end

        # Advance l.
        l = l + blownTs
    end

    # Using the standard scheme
    mul!(@view(w[:, l]), @view(V[1:n, 1:j]), @view(F[1:j, 1]))
    rmul!(@view(w[:, l]), beta)

    # Update tau_out
    tau_now = tau_now + tau

    j = 0
    ireject = 0
    return tau_now, l, j, reject, ireject, step
end

max_stars_repo_name: JuliaDiffEq/ExponentialUtilities

max_issues_repo_name: JuliaDiffEq/ExponentialUtilities

max_forks_repo_name: JuliaDiffEq/ExponentialUtilities


text:
using JuliaLSP

push!(ARGS, "arg")
JuliaLSP.julia_main()