text:
function initialize!(integrator,cache::AN5ConstantCache)
  integrator.kshortsize = 7
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  @inbounds for i in 2:integrator.kshortsize-1
    integrator.k[i] = zero(integrator.fsalfirst)
  end
  integrator.k[integrator.kshortsize] = integrator.fsallast
end

@muladd function perform_step!(integrator, cache::AN5ConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack z,l,m,c_LTE,dts,tsit5tab = cache
  # handle callbacks, rewind back to order one.
  if integrator.u_modified
    cache.order = 1
  end
  # Nordsieck form needs to build the history vector
  if cache.order == 1
    # Start the Nordsieck vector in one shot!
    perform_step!(integrator, tsit5tab, repeat_step)
    cache.order = 4
    z[1] = integrator.uprev
    z[2] = integrator.k[1]*dt
    z[3] = ode_interpolant(t,dt,nothing,nothing,integrator.k,tsit5tab,nothing,Val{2})*dt^2/2
    z[4] = ode_interpolant(t,dt,nothing,nothing,integrator.k,tsit5tab,nothing,Val{3})*dt^3/6
    z[5] = ode_interpolant(t,dt,nothing,nothing,integrator.k,tsit5tab,nothing,Val{4})*dt^4/24
    z[6] = zero(cache.z[6])
    fill!(dts, dt)
    perform_predict!(cache)
    cache.Δ = integrator.u - integrator.uprev
    update_nordsieck_vector!(cache)
    if integrator.opts.adaptive && integrator.EEst >= one(integrator.EEst)
      cache.order = 1
    end
  else
    # Reset time
    tmp = dts[6]
    for i in 5:-1:1
      dts[i+1] = dts[i]
    end
    dts[1] = dt
    dt != dts[2] && nordsieck_rescale!(cache)
    integrator.k[1] = z[2]/dt
    # Perform 5th order Adams method in Nordsieck form
    perform_predict!(cache)
    calc_coeff!(cache)
    isucceed = nlsolve_functional!(integrator, cache)
    if !isucceed
      # rewind Nordsieck vector
      integrator.force_stepfail = true
      nordsieck_rewind!(cache)
      return nothing
    end

    ################################### Error estimation

    if integrator.opts.adaptive
      atmp = calculate_residuals(cache.Δ, uprev, integrator.u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm)
      integrator.EEst = integrator.opts.internalnorm(atmp) * cache.c_LTE
      if integrator.EEst > one(integrator.EEst)
        for i in 1:5
          dts[i] = dts[i+1]
        end
        dts[6] = tmp
      end
    end

    # Correct Nordsieck vector
    cache.order = 5
    update_nordsieck_vector!(cache)

    ################################### Finalize

    integrator.k[2] = cache.z[2]/dt
  end
  return nothing
end

function initialize!(integrator, cache::AN5Cache)
  integrator.kshortsize = 7
  integrator.fsalfirst = cache.tsit5cache.k1; integrator.fsallast = cache.tsit5cache.k7 # setup pointers
  resize!(integrator.k, integrator.kshortsize)
  # Setup k pointers
  integrator.k[1] = cache.tsit5cache.k1
  integrator.k[2] = cache.tsit5cache.k2
  integrator.k[3] = cache.tsit5cache.k3
  integrator.k[4] = cache.tsit5cache.k4
  integrator.k[5] = cache.tsit5cache.k5
  integrator.k[6] = cache.tsit5cache.k6
  integrator.k[7] = cache.tsit5cache.k7
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
end

@muladd function perform_step!(integrator, cache::AN5Cache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p,uprev2 = integrator
  @unpack z,l,m,c_LTE,dts,tmp,ratetmp,atmp,tsit5cache = cache
  # handle callbacks, rewind back to order one.
  if integrator.u_modified
    cache.order = 1
  end
  # Nordsieck form needs to build the history vector
  if cache.order == 1
    ## Start the Nordsieck vector in two shots!
    perform_step!(integrator, tsit5cache, repeat_step)
    copyto!(tmp, integrator.u)
    cache.order = 4
    @. z[1] = integrator.uprev
    @. z[2] = integrator.k[1]*dt
    ode_interpolant!(z[3],t,dt,nothing,nothing,integrator.k,tsit5cache,nothing,Val{2})
    ode_interpolant!(z[4],t,dt,nothing,nothing,integrator.k,tsit5cache,nothing,Val{3})
    ode_interpolant!(z[5],t,dt,nothing,nothing,integrator.k,tsit5cache,nothing,Val{4})
    @. z[3] = z[3]*dt^2/2
    @. z[4] = z[4]*dt^3/6
    @. z[5] = z[5]*dt^4/24
    fill!(z[6], 0)
    fill!(dts, dt)
    perform_predict!(cache)
    @. cache.Δ = integrator.u - integrator.uprev
    update_nordsieck_vector!(cache)
    if integrator.opts.adaptive && integrator.EEst >= one(integrator.EEst)
      cache.order = 1
    end
  else
    # Reset time
    tmp = dts[6]
    for i in 5:-1:1
      dts[i+1] = dts[i]
    end
    dts[1] = dt
    # Rescale
    dt != dts[2] && nordsieck_rescale!(cache)
    @. integrator.k[1] = z[2]/dt
    # Perform 5th order Adams method in Nordsieck form
    perform_predict!(cache)
    calc_coeff!(cache)
    isucceed = nlsolve_functional!(integrator, cache)
    if !isucceed
      integrator.force_stepfail = true
      # rewind Nordsieck vector
      nordsieck_rewind!(cache)
      return nothing
    end

    ################################### Error estimation

    if integrator.opts.adaptive
      calculate_residuals!(atmp, cache.Δ, uprev, integrator.u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm)
      integrator.EEst = integrator.opts.internalnorm(atmp) * cache.c_LTE
      if integrator.EEst > one(integrator.EEst)
        for i in 1:5
          dts[i] = dts[i+1]
        end
        dts[6] = tmp
      end
    end

    # Correct Nordsieck vector
    cache.order = 5
    update_nordsieck_vector!(cache)

    ################################### Finalize

    @. integrator.k[2] = cache.z[2]/dt
  end
  return nothing
end

function initialize!(integrator,cache::JVODEConstantCache)
  integrator.kshortsize = 7
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal

  # Avoid undefined entries if k is an array of arrays
  integrator.fsallast = zero(integrator.fsalfirst)
  integrator.k[1] = integrator.fsalfirst
  @inbounds for i in 2:integrator.kshortsize-1
    integrator.k[i] = zero(integrator.fsalfirst)
  end
  integrator.k[integrator.kshortsize] = integrator.fsallast
end

@muladd function perform_step!(integrator, cache::JVODEConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  @unpack z,l,m,c_LTE,dts,tsit5tab = cache
  # handle callbacks, rewind back to order one.
  if integrator.u_modified || integrator.iter == 1
    cache.order = 1
    z[1] = integrator.uprev
    z[2] = f(uprev, p, t)*dt
    dts[1] = dt
  end
  # Reset time
  tmp = dts[13]
  for i = 12:-1:1
    dts[i+1] = dts[i]
  end
  dts[1] = dt
  dt != dts[2] && nordsieck_adjust!(integrator, cache)
  integrator.k[1] = z[2]/dt

  perform_predict!(cache)
  calc_coeff!(cache)
  isucceed = nlsolve_functional!(integrator, cache)
  if !isucceed
    # rewind Nordsieck vector
    integrator.force_stepfail = true
    nordsieck_rewind!(cache)
    return nothing
  end

  ################################### Error estimation
  if integrator.opts.adaptive
    atmp = calculate_residuals(cache.Δ, uprev, integrator.u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm)
    integrator.EEst = integrator.opts.internalnorm(atmp) * cache.c_LTE
    if integrator.EEst > one(integrator.EEst)
      for i in 1:12
        dts[i] = dts[i+1]
      end
      dts[13] = tmp
    end
  end

  ################################### Finalize
  nordsieck_finalize!(integrator, cache)
  nordsieck_prepare_next!(integrator, cache)
  integrator.k[2] = cache.z[2]/dt
  return nothing
end

function initialize!(integrator, cache::JVODECache)
  integrator.kshortsize = 7
  integrator.fsalfirst = cache.tsit5cache.k1; integrator.fsallast = cache.tsit5cache.k7 # setup pointers
  resize!(integrator.k, integrator.kshortsize)
  # Setup k pointers
  integrator.k[1] = cache.tsit5cache.k1
  integrator.k[2] = cache.tsit5cache.k2
  integrator.k[3] = cache.tsit5cache.k3
  integrator.k[4] = cache.tsit5cache.k4
  integrator.k[5] = cache.tsit5cache.k5
  integrator.k[6] = cache.tsit5cache.k6
  integrator.k[7] = cache.tsit5cache.k7
  integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
end

@muladd function perform_step!(integrator, cache::JVODECache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p,uprev2 = integrator
  @unpack z,l,m,c_LTE,dts,tmp,ratetmp,atmp,tsit5cache = cache
  # handle callbacks, rewind back to order one.
  if integrator.u_modified || integrator.iter == 1
    cache.order = 1
    @. z[1] = integrator.uprev
    f(z[2], uprev, p, t)
    @. z[2] = z[2]*dt
    dts[1] = dt
  end
  # Reset time
  tmp = dts[13]
  for i = 12:-1:1
    dts[i+1] = dts[i]
  end
  dts[1] = dt
  # Rescale
  dt != dts[2] && nordsieck_adjust!(integrator, cache)
  @. integrator.k[1] = z[2]/dt

  perform_predict!(cache)
  calc_coeff!(cache)
  isucceed = nlsolve_functional!(integrator, cache)
  # TODO: Handle NLsolve better
  if !isucceed
    integrator.force_stepfail = true
    # rewind Nordsieck vector
    nordsieck_rewind!(cache)
    return nothing
  end

  ################################### Error estimation

  if integrator.opts.adaptive
    calculate_residuals!(atmp, cache.Δ, uprev, integrator.u, integrator.opts.abstol, integrator.opts.reltol, integrator.opts.internalnorm)
    integrator.EEst = integrator.opts.internalnorm(atmp) * cache.c_LTE
    if integrator.EEst > one(integrator.EEst)
      for i in 1:12
        dts[i] = dts[i+1]
      end
      dts[13] = tmp
    end
  end

  ################################### Finalize

  nordsieck_finalize!(integrator, cache)
  nordsieck_prepare_next!(integrator, cache)
  @. integrator.k[2] = cache.z[2]/dt
  return nothing
end


text:
# Tested in Julia 1.5.1
using DelimitedFiles
data = readdlm("./data/adventInput06.txt", ',', skipblanks= false)
function CustomsYes(data)
    groups = Dict()
    yesAns = []
    c = 1
    lenpass = length(data)
    for i in 1:lenpass
        yans = data[i]
        if yans == "" || i == lenpass
            if i == lenpass
               push!(yesAns, collect(yans))
            end
            groups["group_$(c)"] = yesAns
            yesAns = []
            c += 1
        else
            push!(yesAns, collect(yans))
        end
    end
    groups
end
function YesCount(groups; setOp = union)
    yes_c = 0
    for (k,v) in groups
        yes_c += length(setOp(Set.(v)...))
    end
    yes_c
end
YesCount(CustomsYes(data); setOp = union)
YesCount(CustomsYes(data); setOp = intersect)

ext: jl


text:
# testing some of the functionalities of the fractional flow package

include("../FractionalFlow/FractionalFlow.jl")
FF = FractionalFlow

## 1- read the input file
core_props, fluids, rel_perms, core_flood = FF.read_frac_flow_input("input_file.json")

## 2- run a core flooding

max_stars_repo_name: simulkade/peteng

max_issues_repo_name: simulkade/peteng

max_forks_repo_name: simulkade/peteng


text:
@model binormal() = begin
  y ~ MvNormal(zeros(2), [1.0 0.1; 0.1 1.0])
end


text:
Pkg.add("RDatasets")
Pkg.add("Query")
Pkg.add("StatPlots")
Pkg.add("GR")

using RDatasets, Plots, Query, DataFrames

iris = dataset("datasets", "iris")

function SplitDFbySpecies(iris, species)
    df = @from i in iris begin
         @where i.Species == species
         @select {i.PetalLength, i.PetalWidth}
         @collect DataFrame
    end
    return df
end

function AttributeDescription(df, name)
    println(name, ":\n")
    colnames = names(df)
    ncols = size(colnames)[1]
    for i in 1:ncols
        println(colnames[i], ":")
        println(" mean: ", mean(df[i]))
        println(" var:  ", var(df[i]))
    end
    println()
end

function dfmean(df::DataFrame)
    return [mean(i) for i in df.columns]
end

function dfcov(df)
    cols = df.columns
    ncols = size(cols)[1]
    ndatas = size(cols[1])[1]
    array = zeros(ncols, ncols)
    mean = dfmean(df)
    for i in 1:ncols
        for j in 1:ncols
            array[i, j] = 1/ndatas * (cols[i] .- mean[i])' * (cols[j] .- mean[j])
        end
    end
    return array
end

function scale(df, iris)
    ncols = length(df.columns)
    means = dfmean(iris[[:PetalLength, :PetalWidth]])
    cov = dfcov(iris[[:PetalLength, :PetalWidth]])
    sigmas = [sqrt(cov[i, i]) for i in 1:ncols]

    scaler(x, mu, sigma) = (x .- mu) / sigma
    array = zeros(size(df))
    for i in 1:ncols
        array[:, i] = scaler(df[i], means[i], sigmas[i])
    end
    return array
end

setosa = SplitDFbySpecies(iris, "setosa")
versicolor = SplitDFbySpecies(iris, "versicolor")
virginica = SplitDFbySpecies(iris, "virginica")

gr()
scatter(setosa[1], setosa[2]; xlabel="Petal Length", ylabel="Petal Width",
        xlims=(1, 7), ylims=(-1, 4), m=(0.5, :s), label="setosa")
scatter!(versicolor[1], versicolor[2]; m=(0.5, :c), label="versicolor")
scatter!(virginica[1], virginica[2]; m=(0.5, :v), label="virginica")
savefig("figures/unscaled.png")

setosa_scaled = scale(setosa, iris)
versicolor_scaled = scale(versicolor, iris)
virginica_scaled = scale(virginica, iris)

scatter(setosa_scaled[:, 1], setosa_scaled[:, 2]; xlabel="Petal Length", ylabel="Petal Width",
        xlims=(-1.5, 1.5), ylims=(-1.5, 1.5), m=(0.5, :s), label="setosa")
scatter!(versicolor_scaled[:, 1], versicolor_scaled[:, 2]; m=(0.5, :c), label="versicolor")
scatter!(virginica_scaled[:, 1], virginica_scaled[:, 2]; m=(0.5, :v), label="virginica")
savefig("figures/scaled.png")


text:
using ChainRulesCore

function ChainRulesCore.rrule(
  ::typeof(split_network),
  tn::Matrix{ITensor};
  projector_center=default_projector_center(tn),
  rotation=false,
)
  function pullback(dtn_split::Matrix{ITensor})
    dtn = map(t -> replaceprime(t, 1 => 0), dtn_split)
    return (NoTangent(), dtn)
  end
  return split_network(tn; projector_center=projector_center, rotation=rotation), pullback
end


text:
# Arrays on GPU

typealias CUdeviceptr UInt64

immutable CuPtr
    p::CUdeviceptr

    CuPtr() = new(convert(CUdeviceptr, 0))
    CuPtr(p::CUdeviceptr) = new(p)
end

cubox(p::CuPtr) = cubox(p.p)

function cualloc(T::Type, len::Integer)
    a = CUdeviceptr[0]
    nbytes = Int(len) * sizeof(T)
    @cucall(cuMemAlloc, (Ptr{CUdeviceptr}, Csize_t), a, nbytes)
    return CuPtr(a[1])
end

function free(p::CuPtr)
    @cucall(cuMemFree, (CUdeviceptr,), p.p)
end

isnull(p::CuPtr) = (p.p == 0)

#################################################
#
#  CuArray: contiguous array on GPU
#
#################################################

type CuArray{T,N}
    ptr::CuPtr
    shape::NTuple{N,Int}
    len::Int
end

function CuArray(T::Type, len::Integer)
    n = Int(len)
    p = cualloc(T, n)
    CuArray{T,1}(p, (n,), n)
end

function CuArray{N}(T::Type, shape::NTuple{N,Int})
    n = prod(shape)
    p = cualloc(T, n)
    CuArray{T,N}(p, shape, n)
end

cubox(a::CuArray) = cubox(a.ptr)

length(g::CuArray) = g.len
size(g::CuArray) = g.shape
ndims{T,N}(g::CuArray{T,N}) = N
eltype{T,N}(g::CuArray{T,N}) = T

function size{T,N}(g::CuArray{T,N}, d::Integer)
    d >= 1 ? (d <= N ? g.shape[d] : 1) : error("Invalid index of dimension.")
end

function free(g::CuArray)
    if !isnull(g.ptr)
        free(g.ptr)
        g.ptr = CuPtr()
    end
end

function copy!{T}(dst::Array{T}, src::CuArray{T})
    if length(dst) != length(src)
        throw(ArgumentError("Inconsistent array length."))
    end
    nbytes = length(src) * sizeof(T)
    @cucall(cuMemcpyDtoH, (Ptr{Void}, CUdeviceptr, Csize_t), pointer(dst), src.ptr.p, nbytes)
    return dst
end

function copy!{T}(dst::CuArray{T}, src::Array{T})
    if length(dst) != length(src)
        throw(ArgumentError("Inconsistent array length."))
    end
    nbytes = length(src) * sizeof(T)
    @cucall(cuMemcpyHtoD, (CUdeviceptr, Ptr{Void}, Csize_t), dst.ptr.p, pointer(src), nbytes)
    return dst
end

CuArray{T,N}(a::Array{T,N}) = copy!(CuArray(T, size(a)), a)
to_host{T}(g::CuArray{T}) = copy!(Array(T, size(g)), g)


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

using Markdown
using InteractiveUtils

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

# ╔═╡ 42085492-ac8c-11eb-0620-adcb307077f1
html"""

Section 3.4

 Our first climate model

Lecture Video

"""

# ╔═╡ c7d387fa-cd19-458c-a45d-7893e8c21bbf
begin
    import Pkg
    Pkg.activate(mktempdir())

    Pkg.add([
        Pkg.PackageSpec(name="DifferentialEquations", version="6"),
        Pkg.PackageSpec(name="Plots", version="1"),
		Pkg.PackageSpec(name="PlutoUI", version="0.7"),
	    Pkg.PackageSpec(name="CSV", version="0.8"),
	    Pkg.PackageSpec(name="DataFrames", version="1")

         ])

    using DifferentialEquations, Plots, PlutoUI, LinearAlgebra, CSV, DataFrames
end

# ╔═╡ 14195fc4-40e1-4576-973a-69d649fddc02
TableOfContents(title="📚 Table of Contents", indent=true, depth=4, aside=true)

# ╔═╡ 92883ed9-5572-41fd-96c7-190279f90804
md"""
This lecture is based on the fall 2020 "Introduction to Climate Modeling" lecture by Henri Drake. [Henri's video](https://computationalthinking.mit.edu/Fall20/lecture20/),
[Henri's original notebook](https://github.com/mitmath/18S191/blob/Fall20/lecture_notebooks/week11/01_energy_balance_model.jl) where he describes the Energy Balance Model.
"""

# ╔═╡ 340fc7c4-128f-4476-b445-65005ffa9f5a
md"""
### Pre - Computational Thinking
(traditional ode "analytic" solutions)
"""

# ╔═╡ 7673c7b7-0921-4a7e-ae4d-e785c6391a0c
md"""
Notation:

``y=y(t)`` is a function of time $t$

Inititally $t=0$, $y(0)=y_0$

We write ``y'`` for ``dy/dt``.
"""

# ╔═╡ ebfeb0b0-d29e-4244-a0e3-a7185f3c8124
md"""
---
#### ``y'=`` constant:

``y' = a``
$(html"
") Solution: ``y(t) = at + y_0``

---
"""

# ╔═╡ 4bf60f27-8db5-42c1-83d8-750bec2060e5
md"""

#### ``y'=`` linear in y:

``y' = a - by``
$(html"
") Solution: ``y(t) =(- \frac{a}{b})(e^{-bt}-1) + y_0 e^{-bt}``

Comment:  Equilibrium obtained  (if $b\ne0$)  by solving y'=0 for y giving equilibrium $y=a/b$. Also
can be obtained by letting $t\rightarrow \infty$ in the solution killing the exponential terms.

---
"""

# ╔═╡ e2799d05-36a7-4547-8663-90708ce01ce8
md"""
#### ``y'=`` linear in y + "forcing term" in t:

``y' = a - by + f(t)``
$(html"
")
Solution: ``y(t) = e^{-bt} \left( y_0 + \int_0^t  e^{bu}(a+f(u)) du \right)``

---
"""

# ╔═╡ 4d58018b-752c-43f8-aafb-f11a95ca1874
md"""
### Computational Thinking
"""

# ╔═╡ 279b4748-829f-44a3-94da-da3abf32ff13
md"""
#### Review: Solving an ODE (y'=f(y)) in Julia:

0. Define f(y,p,t)
1. ODEProblem(  f, f(0), time_span, parameters)
2. solve
3. plot etc
"""

# ╔═╡ 578a80ad-fcf1-4d5c-81e5-205113bbf78d
sol = solve(  ODEProblem( f, y₀, (0,10.0), (a,b,forcing(c)) ) );

# ╔═╡ c1d87c42-e1e0-4e52-8023-6cc176266e86
f(y,(a,b,forcing(c)),t) = a - b*y + forcing(c)(t)

# ╔═╡ adb59adc-b30a-4dc4-bc3f-2804b7e02876
forcing(c)  =  t->c*t

# ╔═╡ 220420c1-e8cc-4ff6-8b79-ba2ec49c7695
y₀

# ╔═╡ 7bee45a5-fa62-455b-813a-3e5dcf430289
md"""
a = $(@bind a Slider(0:.01:10, show_value=true, default=0) )

b = $(@bind b Slider(0:.1:5, show_value=true, default=0) )

y₀ = $(@bind y₀ Slider(-5:.1:15, show_value=true, default=2.0) )

c = $(@bind c Slider(0:.1:5, show_value=true, default=0.0) )
"""

# ╔═╡ bda41881-75c7-4732-9a66-d7947607b1b6
begin
	p = plot( sol , legend=false, background_color_inside=:black , ylims=(0,10), lw=3, c=:red)
	# plot direction field:
		xs = Float64[]
		ys = Float64[]

	    lrx = LinRange( xlims(p)..., 30)
		for x in  lrx
			for y in LinRange( ylims(p)..., 30)
				v = [1, a - b*y + forcing(c)(x) ]
				v ./=  (20* (lrx[2]-lrx[1]))
				# plot!([x, x + v[1]], [y, y + v[2]], alpha=0.5, c=:gray)

				push!(xs, x - v[1], x + v[1], NaN)
				push!(ys, y - v[2], y + v[2], NaN)

			end
		end
		hline!( [b==0 ? 0 : a/b],c=:white,ls=:dash)
		plot!(xs, ys, alpha=0.7, c=:yellow)
	    ylabel!("y")
	    annotate!(-.5,y₀,text("y₀",color=:red))
	    title!("Solution to y'(t) = a - by + forcing(c)(t)")
end

# ╔═╡ aba7fc93-c0ac-4c9e-b975-18145f87707f
md"""#### 1) Background: climate physics

The simplest climate model can be conceptualized as:
\begin{align}
\text{\color{brown}{change in heat content}} = & + \text{\color{orange}{absorbed solar radiation (energy from the Sun's rays)}} \newline
& - \text{\color{blue}{outgoing thermal radiation (i.e. blackbody cooling to space)}}
\newline
& + \text{\color{grey}{human-caused greenhouse effect (trapped outgoing radiation)}}
\end{align}

where each of these is interpreted as an average over the entire globe (hence "zero-dimensional").
"""

# ╔═╡ b9eca29e-9028-4fd7-8c62-718a2dcf87d1
html""""""

# ╔═╡ dabca25c-e34b-4aed-8729-132d03c5bb45
md"""
### 1.1 Incoming 🌞: Absorbed solar radiation
(an example of $\mathrm{temp}$'=constant)
"""

# ╔═╡ d74936e9-b760-4add-b3e7-46d544064c16
md"""
In math we just write down a differential equation, but in the physical world there are physical variables to identify.

In our baking the earth example, we will identify the following quantities:

- Industrial Revolution Start: 1850
- Avg Temperature in 1850: 14.0 °C

- Solar Insolation $S=$1368 W/m^2:  energy from the sun

- Albedo or plentary reflectivity: α = 0.3

- atmosphere and upper-ocean heat capacity: C= 51  J/m^2/°C

Earth Baking Formula:
$(html"
")
 `` C\  \mathrm{temp}'(t) = S(1-α)/4 = `` $(S*(1-α)/4)
"""

# ╔═╡ 5123525a-3437-4b76-813c-8ad6b158f7f2
md"""
##### (Heating the earth nonstop)
"""

# ╔═╡ 2f19cbac-4f13-4c2b-9b11-fb92e8055527
begin
	plot(solve(p1),       legend = false,
		 background_color_inside = :black,
		                  xlabel = "years from $(1850)",
	                      ylabel = "Temperature °C")
	hline!( [temp₀,temp₀] ,c=:white,ls=:dash)
	annotate!( 80, 25+temp₀, text("Preindustrial Temperature = $(temp₀)°C",color=:white))
	title!("Absorbing Solar Radiation (only)")
end

# ╔═╡ 087f47b2-8283-4205-88f2-4d5883a340c2
md"""

At Earth's orbital distance from the Sun, the power of the Sun's rays that intercept the Earth is equal to
"""

# ╔═╡ 9c89c4e9-65ee-4424-bd74-17168b211797
S = 1368; # solar insolation [W/m^2]  (energy per unit time per unit area)

# ╔═╡ ca5b41d2-b511-486a-91fe-cceb8f7282c3
md"A small fraction"

# ╔═╡ 13984c61-4c34-40db-9043-fcff2721522e
α = 0.3; # albedo, or planetary reflectivity [unitless]

# ╔═╡ 868e19f4-d71e-4222-9bdd-470387991c67
md"""
of this incoming solar radiation is reflected back out to space (by reflective surfaces like white clouds, snow, and ice), with the remaining fraction $(1-\alpha)$ being absorbed.

Since the incoming solar rays are all approximately parallel this far from the Sun, the cross-sectional area of the Earth that intercepts them is just a disc of area $\pi R^{2}$. Since all of the other terms we will consider act on the entire surface area $4\pi R^{2}$ of the spherical Earth, the absorbed solar radiation *per unit surface area* (averaged over the entire globe) is reduced by a factor of 4.

![](https://www.open.edu/openlearn/ocw/pluginfile.php/101161/mod_oucontent/oucontent/890/639dcd57/ce3f1c3a/s250_3_002i.jpg)

The absorbed solar radiation per unit area is thus

$\textcolor{orange}{\text{absorbed solar radiation} \equiv \frac{S(1-\alpha)}{4}}$
"""

# ╔═╡ 0e4dedc5-0b1e-4e46-b943-284bfb2de57f
absorbed_solar_radiation = S*(1 - α)/4; # [W/m^2]

# ╔═╡ 91329865-f593-4388-a5d0-01af4be6e01d
begin
	C = 51.; # atmosphere and upper-ocean heat capacity [J/m^2/°C]
	temp₀ = 14.0 # preindustrial temperature [°C]
end

# ╔═╡ 7187ae25-239d-4752-898e-6674009b5de6
 p1 = ODEProblem( (temp,p,t)-> (1/C) * absorbed_solar_radiation, temp₀,  (0.0,170) )

# ╔═╡ ad1e294e-ad8a-48f9-b924-2474cce16aaf
md"""The heat content $C temp$ is determined by the temperature $temp$ (in Kelvin) and the heat capacity of the climate system. While we are interested in the temperature of the atmosphere, which has a very small heat capacity, its heat is closely coupled with that of the upper ocean, which has a much larger heat capacity of
"""

# ╔═╡ 2086e386-194b-4863-a9ee-178002925d42
md"""
The *change in heat content over time* is thus simply given by $\frac{d(CT)}{dt}$. Since the heat capacity of sea water hardly changes with temperature, we can rewrite this in terms of the change in temperature with time as:

$\color{brown}{\text{change in heat content } =\; C \frac{dtemp}{dt}}$
"""

# ╔═╡ 2b66ffaf-a9a3-4ca2-a6e9-732912c9d48e
md"""
### 1.2 Outgoing ♨ : thermal radiation
"""

# ╔═╡ b186d026-d02c-4f06-8097-aa54b363f4fb
md"""

The outgoing thermal radiation term $\mathcal{G}(T)$ (or "blackbody cooling to space") represents the combined effects of *negative feedbacks that dampen warming*, such as **blackbody radiation**, and *positive feedbacks that amplify warming*, such as the **water vapor feedback**.

Since these physics are too complicated to deal with here, we *linearize* the model
combining the incoming and the outgoing.

We assume that the preindustrial world was in energy balance, and thus
the equilibrium temperature is the preindustrial temperature.

Thus we assume

 temp'($t$) = $B$(temp($0$)-temp($t$))  for some value of $B$.
The minus sign in front of temp(t) indicating it restores equilibrium.
"""

# ╔═╡ edafca88-dd4a-4fe4-8ffc-bb50f7bd561d
md"""
The value that has been chosen is
"""

# ╔═╡ 9b97e41b-76ff-4cc2-90c1-446d59c4ece7
B = 1.3; # climate feedback parameter [W/m^2/°C],

# ╔═╡ af893e1e-3bd7-4d05-83b6-f5f7e751b5e3
md"""
start\_temp = $(@bind start_temp Slider(0:28; show_value=true, default=14))
"""

# ╔═╡ 333290b0-0198-4306-b19b-6d30df79a280
p2 = ODEProblem( (temp,p,t)-> (1/C) * B*(temp₀-temp), start_temp,  (0.0,170) )

# ╔═╡ eeaf3735-5139-4924-9fce-14df51a61042
begin
	plot(solve(p2),       legend = false,
		 background_color_inside = :black,
		                  xlabel = "years from start",
	                      ylabel = "Temperature °C",
	                      ylim = (0,30))
	hline!( [temp₀,temp₀] ,c=:white,ls=:dash)
	annotate!( 80, temp₀, text("Preindustrial Temperature = $(temp₀)°C",:bottom,color=:white))
	title!("Energy Balance Model (Healthy Earth)")
end

# ╔═╡ 20d6c513-2ca6-4dea-9092-156e2805d467
md"""
### 1.3 Greenhouse 🏭: Human-caused greenhouse effect

Empirically, the greenhouse effect is known to be a logarithmic function of gaseous carbon dioxide (CO₂) concentrations

$\color{grey}{\text{human-caused greenhouse effect}\; = {\mbox {(forcing\_coef)}} \ln \left( \frac{[\text{CO}₂]}{[\text{CO}₂]_{\text{PreIndust}}} \right),}$

How this depends on time into the future depends on human behavior!
Time is not modelled in the above equation.
"""

# ╔═╡ 692a6928-c7a2-4b61-a794-16bef5d4e919
md"where"

# ╔═╡ 6d1058bf-8a05-4b8e-835a-de9e95f567c4
forcing_coef = 5.0; # CO2 forcing coefficient [W/m^2]

# ╔═╡ 35184e6f-ac3e-480b-92e9-ed8baaaa833b
CO₂_PreIndust = 280.; # preindustrial CO2 concentration [parts per million; ppm];

# ╔═╡ 437faadd-0301-403a-bcd7-18ce279589d0
greenhouse_effect(CO₂) = forcing_coef * log(CO₂/CO₂_PreIndust)

# ╔═╡ 1a4d21bd-85ad-4935-913a-8992a8996db4
greenhouse_effect(CO₂(15))

# ╔═╡ ee6414b7-e92d-4055-af17-6b02f05c28cd
begin
	CO2_present = 420.
	CO2_range = 280*(2 .^(range(-1, stop=3,length=100)))
	plot(CO2_range, greenhouse_effect.(CO2_range), lw=2.5, label=nothing, color=:black)
	plot!([CO₂_PreIndust], [greenhouse_effect(CO₂_PreIndust)], marker=:., ms=6, linecolor=:white,
		color=:blue, lw=0, label="pre-industrial (PI)")
	plot!([CO2_present], [greenhouse_effect(CO2_present)], marker=:., ms=6, color=:red, linecolor=:white, lw=0, label="present day (2020)")
	plot!(xticks=[280, 280*2, 280*4, 280*8], legend=:bottomright, size=(400, 250))
	plot!(ylabel="Radiative forcing [W/m²]", xlabel="CO₂ concentration [ppm]")
end

# ╔═╡ fac5012a-960c-473c-bbf0-62c0be87f608
begin
	 #CO₂(t) = CO₂_PreIndust # no emissions
	 # CO₂(t) = CO₂_PreIndust * 1.01^t # test model
	 CO₂(t) = CO₂_PreIndust * (1+ (t/220)^3 )
end

# ╔═╡ 99629ec2-dc70-4253-b191-305bccc9f36b
p3 = ODEProblem( (temp,p,t)-> (1/C) * ( B*(temp₀-temp)  + greenhouse_effect(CO₂(t))    ) , start_temp,  (0.0,170) )

# ╔═╡ 6b2beeec-6383-42b3-b694-8d77b961c8a1
begin
	plot(solve(p3),       legend = false,
		 background_color_inside = :black,
		                  xlabel = "years from 1850",
	                      ylabel = "Temperature °C",
	                      ylim = (10,20))
	hline!( [temp₀,temp₀] ,c=:white,ls=:dash)
	annotate!( 80, temp₀, text("Preindustrial Temperature = $(temp₀)°C",:bottom,color=:white))
	title!("Model with CO₂")
end

# ╔═╡ 0b24f105-0166-4a41-97aa-156417d7203a
begin
	years = 1850:2020
	plot( years, CO₂.(years.-1850), lw=3, legend=false)
end

# ╔═╡ c6f8dcf6-950d-48ff-b040-55bdd347d74b
md"""
### Observations from Mauna Loa Volcano ![Mauna Loa Volcano](https://i.pinimg.com/originals/df/1a/e7/df1ae72cfd5e6d0d535c0ec99e708f6f.jpg)
"""

# ╔═╡ cd1dcbc4-2273-4eda-85d9-9af4fd71b3c1
md"""
information is available at
[https://www.ncei.noaa.gov/pub/data/paleo/icecore/antarctica/law/law2006.txt]
(https://www.ncei.noaa.gov/pub/data/paleo/icecore/antarctica/law/law2006.txt).
"""

# ╔═╡ c3b0b7fc-19be-4f04-8787-6349ab9bff7f
begin
	CO2_historical_data_url = "https://scrippsco2.ucsd.edu/assets/data/atmospheric/stations/in_situ_co2/monthly/monthly_in_situ_co2_mlo.csv"

	CO2_historical_data = CSV.read(download(CO2_historical_data_url), DataFrame, header=4, datarow=7, );

	first(CO2_historical_data, 11)
end

# ╔═╡ f5c2d0de-f5d2-43f2-bfca-41bd727b3ca9
md"""
Data is in the fifth column CO₂
"""

# ╔═╡ ef9d9512-cb40-4641-9d6a-5a36f5a4b33d
md"""
Oh no, missing data (-99.99)
"""

# ╔═╡ 2bc4c596-b313-49f6-84b5-5d53f2baf0e9
validrowsmask = CO2_historical_data[:, "     CO2"] .> 0

# ╔═╡ 648dfce3-5363-4539-8f41-f5d27bae4c6f
begin
	begin
		plot( CO2_historical_data[validrowsmask, "      Date"] , CO2_historical_data[validrowsmask, "     CO2"], label="Mauna Loa CO₂ data (Keeling curve)")
		plot!( years, CO₂.(years.-1850), lw=3 , label="Cubic Fit", legend=:topleft)
		xlabel!("year")
		ylabel!("CO₂ (ppm)")
		title!("CO₂ observations and fit")

	end
end

# ╔═╡ 288bfc6b-6848-40d5-916b-95a8a4248a4f
md"""
We will use this fit to compare against historical temperatures.
"""

# ╔═╡ a62ca8d1-321a-435e-a9df-f63d000376c7
begin
	plot(years,solp4.(years.-1850),lw=2,label="Predicted Temperature from model", legend=:topleft)
	xlabel!("year")
	ylabel!("Temp °C")

	plot!( parse.(Float64, T_df[:,1]), parse.(Float64, T_df[:,2]) .+ 14.15, color=:black, label="NASA Observations", legend=:topleft)
end

# ╔═╡ 3304174c-289d-47c5-b5ef-161b11e515eb
md"""
Climate feedback BB = $(@bind BB Slider(0:.1:4, show_value=true, default=B))

Ocean Heat Capacity CC =$(@bind CC Slider(10:.1:200, show_value=true, default=C))
"""

# ╔═╡ 13269c86-06a4-4354-b620-bdf3f6432294
p4 = ODEProblem( (temp,p,t)-> (1/CC) * ( BB*(temp₀-temp)  + greenhouse_effect(CO₂(t))    ) , start_temp,  (0.0,170) )

# ╔═╡ 28acb5a4-2a5f-49c5-9c78-deb40fdeed36
md""" #####  Best- and worst-case projections of future global warming

"""

# ╔═╡ 13b003d2-1fd4-4a4a-960c-4a1d9b673dc6
md"""Consider two divergent hypothetical futures:
1. a **low-emissions** world in which emissions decrease such that CO2 concentrations stay below 500 ppm by 2100 (known in climate circles as "RCP8.5") and
2. a **high-emissions** world in which emissions continue increasing and CO2 concentrations soar upwards of 1200 ppm ("RCP2.6").
"""

# ╔═╡ a2288816-3621-4871-9faf-3e9c78674969
md"""
![](https://raw.githubusercontent.com/mitmath/18S191/Spring21/notebooks/week12/predictthefuture.svg)
"""

# ╔═╡ 2b9c3427-c3ec-41ba-b3c5-268cf3dddccf
md"""
In the low-emissions scenario, the temperature increase stays below $ΔT = 2$ °C by 2100, while in the high-emissions scenario temperatures soar upwards of 3.5ºC above pre-industrial levels.
"""

# ╔═╡ fef0a76b-98dc-42a0-975b-0fba8a32e5cc
md"Although the greenhouse effect due to human-caused CO₂ emissions is the dominant forcing behind historical and future-projected warming, modern climate modelling considers a fairly exhaustive list of other forcing factors (aerosols, other greenhouse gases, ozone, land-use changes, etc.). The video below shows a breakdown of these forcing factors in a state-of-the-art climate model simulation of the historical period."

# ╔═╡ 2368f16c-9805-4dd6-a130-d831211f6155
html"""

"""

# ╔═╡ 6800c94d-e7f3-43a8-a823-1c550cb8fc95
solp4 = solve(p4)

# ╔═╡ b38727a9-72ff-499f-9751-ddfeee21959c
begin
	T_url = "https://data.giss.nasa.gov/gistemp/graphs/graph_data/Global_Mean_Estimates_based_on_Land_and_Ocean_Data/graph.txt";
	T_df = CSV.read(download(T_url),DataFrame, header=false, datarow=6,delim="     ");
    #T_df = T_df[:,[1,6]]

end

# ╔═╡ Cell order:
# ╟─42085492-ac8c-11eb-0620-adcb307077f1
# ╟─c7d387fa-cd19-458c-a45d-7893e8c21bbf
# ╟─14195fc4-40e1-4576-973a-69d649fddc02
# ╟─92883ed9-5572-41fd-96c7-190279f90804
# ╟─340fc7c4-128f-4476-b445-65005ffa9f5a
# ╟─7673c7b7-0921-4a7e-ae4d-e785c6391a0c
# ╟─ebfeb0b0-d29e-4244-a0e3-a7185f3c8124
# ╟─4bf60f27-8db5-42c1-83d8-750bec2060e5
# ╟─e2799d05-36a7-4547-8663-90708ce01ce8
# ╟─4d58018b-752c-43f8-aafb-f11a95ca1874
# ╟─279b4748-829f-44a3-94da-da3abf32ff13
# ╠═578a80ad-fcf1-4d5c-81e5-205113bbf78d
# ╠═c1d87c42-e1e0-4e52-8023-6cc176266e86
# ╠═adb59adc-b30a-4dc4-bc3f-2804b7e02876
# ╠═220420c1-e8cc-4ff6-8b79-ba2ec49c7695
# ╟─7bee45a5-fa62-455b-813a-3e5dcf430289
# ╟─bda41881-75c7-4732-9a66-d7947607b1b6
# ╟─aba7fc93-c0ac-4c9e-b975-18145f87707f
# ╟─b9eca29e-9028-4fd7-8c62-718a2dcf87d1
# ╟─dabca25c-e34b-4aed-8729-132d03c5bb45
# ╟─d74936e9-b760-4add-b3e7-46d544064c16
# ╟─5123525a-3437-4b76-813c-8ad6b158f7f2
# ╠═2f19cbac-4f13-4c2b-9b11-fb92e8055527
# ╟─087f47b2-8283-4205-88f2-4d5883a340c2
# ╠═9c89c4e9-65ee-4424-bd74-17168b211797
# ╟─ca5b41d2-b511-486a-91fe-cceb8f7282c3
# ╠═13984c61-4c34-40db-9043-fcff2721522e
# ╟─868e19f4-d71e-4222-9bdd-470387991c67
# ╠═0e4dedc5-0b1e-4e46-b943-284bfb2de57f
# ╠═91329865-f593-4388-a5d0-01af4be6e01d
# ╠═7187ae25-239d-4752-898e-6674009b5de6
# ╟─ad1e294e-ad8a-48f9-b924-2474cce16aaf
# ╟─2086e386-194b-4863-a9ee-178002925d42
# ╟─2b66ffaf-a9a3-4ca2-a6e9-732912c9d48e
# ╟─b186d026-d02c-4f06-8097-aa54b363f4fb
# ╟─edafca88-dd4a-4fe4-8ffc-bb50f7bd561d
# ╠═9b97e41b-76ff-4cc2-90c1-446d59c4ece7
# ╟─af893e1e-3bd7-4d05-83b6-f5f7e751b5e3
# ╠═333290b0-0198-4306-b19b-6d30df79a280
# ╠═eeaf3735-5139-4924-9fce-14df51a61042
# ╟─20d6c513-2ca6-4dea-9092-156e2805d467
# ╟─692a6928-c7a2-4b61-a794-16bef5d4e919
# ╠═6d1058bf-8a05-4b8e-835a-de9e95f567c4
# ╠═35184e6f-ac3e-480b-92e9-ed8baaaa833b
# ╠═437faadd-0301-403a-bcd7-18ce279589d0
# ╠═1a4d21bd-85ad-4935-913a-8992a8996db4
# ╟─ee6414b7-e92d-4055-af17-6b02f05c28cd
# ╠═fac5012a-960c-473c-bbf0-62c0be87f608
# ╠═99629ec2-dc70-4253-b191-305bccc9f36b
# ╠═6b2beeec-6383-42b3-b694-8d77b961c8a1
# ╠═0b24f105-0166-4a41-97aa-156417d7203a
# ╟─c6f8dcf6-950d-48ff-b040-55bdd347d74b
# ╟─cd1dcbc4-2273-4eda-85d9-9af4fd71b3c1
# ╠═c3b0b7fc-19be-4f04-8787-6349ab9bff7f
# ╟─f5c2d0de-f5d2-43f2-bfca-41bd727b3ca9
# ╟─ef9d9512-cb40-4641-9d6a-5a36f5a4b33d
# ╠═2bc4c596-b313-49f6-84b5-5d53f2baf0e9
# ╟─648dfce3-5363-4539-8f41-f5d27bae4c6f
# ╟─288bfc6b-6848-40d5-916b-95a8a4248a4f
# ╠═a62ca8d1-321a-435e-a9df-f63d000376c7
# ╟─3304174c-289d-47c5-b5ef-161b11e515eb
# ╠═13269c86-06a4-4354-b620-bdf3f6432294
# ╟─28acb5a4-2a5f-49c5-9c78-deb40fdeed36
# ╟─13b003d2-1fd4-4a4a-960c-4a1d9b673dc6
# ╟─a2288816-3621-4871-9faf-3e9c78674969
# ╟─2b9c3427-c3ec-41ba-b3c5-268cf3dddccf
# ╟─fef0a76b-98dc-42a0-975b-0fba8a32e5cc
# ╟─2368f16c-9805-4dd6-a130-d831211f6155
# ╠═6800c94d-e7f3-43a8-a823-1c550cb8fc95
# ╠═b38727a9-72ff-499f-9751-ddfeee21959c


text:
@doc doc"""
    SymmetricPositiveDefinite{N} <: Manifold

The manifold of symmetric positive definite matrices, i.e.

```math
\mathcal P(n) =
\bigl\{
x \in \mathbb R^{n\times n} :
\xi^\mathrm{T}x\xi > 0 \text{ for all } \xi \in \mathbb R^{n}\backslash\{0\}
\bigr\}
```

# Constructor

    SymmetricPositiveDefinite(n)

generates the manifold $\mathcal P(n) \subset \mathbb R^{n\times n}$
"""
struct SymmetricPositiveDefinite{N} <: Manifold end

SymmetricPositiveDefinite(n::Int) = SymmetricPositiveDefinite{n}()

@doc doc"""
    check_manifold_point(M::SymmetricPositiveDefinite, x; kwargs...)

checks, whether `x` is a valid point on the [`SymmetricPositiveDefinite`](@ref) `M`, i.e. is a matrix
of size `(N,N)`, symmetric and positive definite.
The tolerance for the second to last test can be set using the `kwargs...`.
"""
function check_manifold_point(M::SymmetricPositiveDefinite{N}, x; kwargs...) where {N}
    if size(x) != representation_size(M)
        return DomainError(
            size(x),
            "The point $(x) does not lie on $(M), since its size is not $(representation_size(M)).",
        )
    end
    if !isapprox(norm(x - transpose(x)), 0.0; kwargs...)
        return DomainError(
            norm(x),
            "The point $(x) does not lie on $(M) since its not a symmetric matrix:",
        )
    end
    if !all(eigvals(x) .> 0)
        return DomainError(
            norm(x),
            "The point $x does not lie on $(M) since its not a positive definite matrix.",
        )
    end
    return nothing
end

"""
    check_tangent_vector(M::SymmetricPositiveDefinite, x, v; kwargs... )

Check whether `v` is a tangent vector to `x` on the [`SymmetricPositiveDefinite`](@ref) `M`,
i.e. atfer [`check_manifold_point`](@ref)`(M,x)`, `v` has to be of same dimension as `x`
and a symmetric matrix, i.e. this stores tangent vetors as elements of the corresponding
Lie group. The tolerance for the last test can be set using the `kwargs...`.
"""
function check_tangent_vector(M::SymmetricPositiveDefinite{N}, x, v; kwargs...) where {N}
    mpe = check_manifold_point(M, x)
    mpe === nothing || return mpe
    if size(v) != representation_size(M)
        return DomainError(
            size(v),
            "The vector $(v) is not a tangent to a point on $(M) since its size does not match $(representation_size(M)).",
        )
    end
    if !isapprox(norm(v - transpose(v)), 0.0; kwargs...)
        return DomainError(
            size(v),
            "The vector $(v) is not a tangent to a point on $(M) (represented as an element of the Lie algebra) since its not symmetric.",
        )
    end
    return nothing
end

@doc doc"""
    injectivity_radius(M::SymmetricPositiveDefinite[, x])
    injectivity_radius(M::MetricManifold{SymmetricPositiveDefinite,LinearAffineMetric}[, x])
    injectivity_radius(M::MetricManifold{SymmetricPositiveDefinite,LogCholeskyMetric}[, x])

Return the injectivity radius of the [`SymmetricPositiveDefinite`](@ref).
Since `M` is a Hadamard manifold with respect to the [`LinearAffineMetric`](@ref) and the
[`LogCholeskyMetric`](@ref), the injectivity radius is globally $\infty$.
"""
injectivity_radius(M::SymmetricPositiveDefinite{N}, args...) where {N} = Inf

@doc doc"""
    manifold_dimension(M::SymmetricPositiveDefinite)

returns the dimension of
[`SymmetricPositiveDefinite`](@ref) `M`$=\mathcal P(n), n\in \mathbb N$, i.e.
````math
\dim \mathcal P(n) = \frac{n(n+1)}{2}
````
"""
@generated function manifold_dimension(M::SymmetricPositiveDefinite{N}) where {N}
    return div(N * (N + 1), 2)
end

"""
    mean(
        M::SymmetricPositiveDefinite,
        x::AbstractVector,
        [w::AbstractWeights,]
        method = GeodesicInterpolation();
        kwargs...,
    )

Compute the Riemannian [`mean`](@ref mean(M::Manifold, args...)) of `x` using
[`GeodesicInterpolation`](@ref).
"""
mean(::SymmetricPositiveDefinite, ::Any)

function mean!(
    M::SymmetricPositiveDefinite,
    y,
    x::AbstractVector,
    w::AbstractVector;
    kwargs...,
)
    return mean!(M, y, x, w, GeodesicInterpolation(); kwargs...)
end

@doc doc"""
    representation_size(M::SymmetricPositiveDefinite)

Return the size of an array representing an element on the
[`SymmetricPositiveDefinite`](@ref) manifold `M`, i.e. $n\times n$, the size of such a
symmetric positive definite matrix on $\mathcal M = \mathcal P(n)$.
"""
@generated representation_size(::SymmetricPositiveDefinite{N}) where {N} = (N, N)

@doc doc"""
    zero_tangent_vector(M::SymmetricPositiveDefinite,x)

returns the zero tangent vector in the tangent space of the symmetric positive
definite matrix `x` on the [`SymmetricPositiveDefinite`](@ref) manifold `M`.
"""
zero_tangent_vector(M::SymmetricPositiveDefinite, x) = zero(x)
function zero_tangent_vector(
    M::MetricManifold{SymmetricPositiveDefinite{N},T},
    x,
) where {N,T<:Metric}
    return zero(x)
end

zero_tangent_vector!(M::SymmetricPositiveDefinite{N}, v, x) where {N} = fill!(v, 0)


text:
# ------------------------------------------------------------------------------
# Space war game!
mutable struct Position
    x::Int
    y::Int
end

struct Size
    width::Int
    height::Int
end

# Widgets hold a single game piece
struct Widget
    name::String
    position::Position
    size::Size
end

# Define pretty print functions
# Base.show(io::IO, p::Position) = print(io, "(", p.x, ",", p.y, ")")
# Base.show(io::IO, s::Size) = print(io, s.width, " x ", s.height)
# Base.show(io::IO, w::Widget) = print(io, w.name, " at ", w.position, " size ", w.size)

# single-line functions
# Note the fact that we haven't annotated these functions with widget::Widget
# This is because that information is not strictly necessary. For method signatures
# types have no effect on performance. They are only useful to control dispatch.
move_up!(widget, v)    = widget.position.y -= v
move_down!(widget, v)  = widget.position.y += v
move_left!(widget, v)  = widget.position.x -= v
move_right!(widget, v) = widget.position.x += v

# long version
function move_up!(widget, v)
    widget.position.y -= v
end

# behavior without annotating type in arguments: error message is somewhat cryptic
move_up!(1, 2)
#= REPL
julia> move_up!(1, 2)
ERROR: type Int64 has no field position
=#

# Now, restart REPL.
# Re-evaluate the struct definitions.
# Then, define the following functions.
move_up!(widget::Widget, v::Int)    = widget.position.y -= v
move_down!(widget::Widget, v::Int)  = widget.position.y += v
move_left!(widget::Widget, v::Int)  = widget.position.x -= v
move_right!(widget::Widget, v::Int) = widget.position.x += v
move_up!(1, 2)

#= REPL
julia> move_up!(1, 2)
ERROR: MethodError: no method matching move_up!(::Int64, ::Int64)
Closest candidates are:
  move_up!(::Widget, ::Int64) at REPL[5]:1
=#

# Define pretty print functions
Base.show(io::IO, p::Position) = print(io, "(", p.x, ",", p.y, ")")
Base.show(io::IO, s::Size) = print(io, s.width, " x ", s.height)
Base.show(io::IO, w::Widget) = print(io, w.name, " at ", w.position, " size ", w.size)

#= REPL
julia> w = Widget("asteroid", Position(0, 0), Size(10, 20))
asteroid at (0,0) size 10 x 20

julia> move_up!(w, 10)
-10

julia> move_down!(w, 10)
0

julia> move_left!(w, 20)
-20

julia> move_right!(w, 20)
0

julia> print(w) # should be back to position (0,0)
asteroid at (0,0) size 10 x 20
=#

# ------------------------------------------------------------------------------
# Optional arguments

# Make a bunch of asteroids
function make_asteroids(N::Int, pos_range = 0:200, size_range = 10:30)
    pos_rand() = rand(pos_range)
    sz_rand() = rand(size_range)
    return [Widget("Asteroid #$i",
                Position(pos_rand(), pos_rand()),
                Size(sz_rand(), sz_rand()))
        for i in 1:N]
end

asteroids = make_asteroids(5)
#= REPL
julia> asteroids = make_asteroids(5)
5-element Array{Widget,1}:
Asteroid #1 at (129,172) size 20 x 17
Asteroid #2 at (93,143) size 29 x 22
Asteroid #3 at (157,152) size 26 x 14
Asteroid #4 at (120,3) size 26 x 15
Asteroid #5 at (5,178) size 20 x 26
=#

asteroids = make_asteroids(5, 1:10)
#= REPL
julia> asteroids = make_asteroids(5, 1:10)
5-element Array{Widget,1}:
Asteroid #1 at (4,10) size 10 x 20
Asteroid #2 at (4,5) size 27 x 14
Asteroid #3 at (3,1) size 12 x 16
Asteroid #4 at (5,9) size 28 x 26
Asteroid #5 at (4,3) size 11 x 28
=#

# --- how to find methods? ---

#= REPL
julia> make_asteroids
make_asteroids (generic function with 3 methods)

julia> make_asteroids(
make_asteroids(N::Int64) in Main at REPL[20]:2
make_asteroids(N::Int64, pos_range) in Main at REPL[20]:2
make_asteroids(N::Int64, pos_range, size_range) in Main at REPL[20]:2

julia> methods(make_asteroids)
# 3 methods for generic function "make_asteroids":
[1] make_asteroids(N::Int64) in Main at REPL[20]:2
[2] make_asteroids(N::Int64, pos_range) in Main at REPL[20]:2
[3] make_asteroids(N::Int64, pos_range, size_range) in Main at REPL[20]:2
=#

# --- fully specify arguments ---
#= REPL
julia> asteroids = make_asteroids(5, 100:5:200, 200:10:500)
5-element Array{Widget,1}:
 Asteroid #1 at (180,105) size 290 x 390
 Asteroid #2 at (145,115) size 410 x 200
 Asteroid #3 at (195,115) size 460 x 410
 Asteroid #4 at (200,185) size 450 x 370
 Asteroid #5 at (115,180) size 450 x 360
=#

# ------------------------------------------------------------------------------
# keyword arguments

function make_asteroids2(N::Int; pos_range = 0:200, size_range = 10:30)
    pos_rand() = rand(pos_range)
    sz_rand() = rand(size_range)
    return [Widget("Asteroid #$i",
                Position(pos_rand(), pos_rand()),
                Size(sz_rand(), sz_rand()))
        for i in 1:N]
end

asteroids = make_asteroids2(5, pos_range = 0:100:500)
#= REPL
julia> asteroids = make_asteroids2(5, pos_range = 0:100:500)
5-element Array{Widget,1}:
 Asteroid #1 at (200,500) size 29 x 25
 Asteroid #2 at (100,400) size 18 x 28
 Asteroid #3 at (400,500) size 13 x 10
 Asteroid #4 at (200,500) size 12 x 23
 Asteroid #5 at (200,500) size 16 x 12
=#

asteroids = make_asteroids2(5, size_range = 1:5, pos_range = 0:10:100)
#=
julia> asteroids = make_asteroids2(5, size_range = 1:5, pos_range = 0:10:100)
5-element Array{Widget,1}:
 Asteroid #1 at (10,90) size 4 x 1
 Asteroid #2 at (30,20) size 5 x 3
 Asteroid #3 at (10,70) size 3 x 5
 Asteroid #4 at (70,70) size 5 x 2
 Asteroid #5 at (0,60) size 3 x 3
 =#

function make_asteroids3(; N::Int, pos_range = 0:200, size_range = 10:30)
    pos_rand() = rand(pos_range)
    sz_rand() = rand(size_range)
    return [Widget("Asteroid #$i",
                Position(pos_rand(), pos_rand()),
                Size(sz_rand(), sz_rand()))
        for i in 1:N]
end

#= REPL
julia> make_asteroids3(N = 3)
3-element Array{Widget,1}:
 Asteroid #1 at (100,46) size 13 x 14
 Asteroid #2 at (168,132) size 26 x 20
 Asteroid #3 at (135,104) size 19 x 10
=#

# ------------------------------------------------------------------------------
# Variable arguments, aka "slurping"

# Shoot any number of targets
function shoot(from::Widget, targets::Widget...)
    println("Type of targets: ", typeof(targets))
    for target in targets
        println(from.name, " --> ", target.name)
    end
end

spaceship = Widget("Spaceship", Position(0, 0), Size(30,30))
target1 = asteroids[1]
target2 = asteroids[2]
target3 = asteroids[3]

shoot(spaceship, target1)
shoot(spaceship, target1, target2, target3)
# julia> shoot(spaceship, target1)
# Type of targets: Tuple{Widget}
# Spaceship --> Asteroid #1

# julia> shoot(spaceship, target1, target2, target3)
# Type of targets: Tuple{Widget,Widget,Widget}
# Spaceship --> Asteroid #1
# Spaceship --> Asteroid #2
# Spaceship --> Asteroid #3

# ------------------------------------------------------------------------------
# splatting

# Special arrangement before attacks
function triangular_formation!(s1::Widget, s2::Widget, s3::Widget)
    x_offset = 30
    y_offset = 50
    s2.position.x = s1.position.x - x_offset
    s3.position.x = s1.position.x + x_offset
    s2.position.y = s3.position.y = s1.position.y - y_offset
    (s1, s2, s3)
end

spaceships = [Widget("Spaceship $i", Position(0,0), Size(20, 50))
                for i in 1:3]
triangular_formation!(spaceships...);
spaceships

#= REPL
julia> spaceships = [Widget("Spaceship $i", Position(0,0), Size(20, 50))
                       for i in 1:3]
3-element Array{Widget,1}:
 Spaceship 1 at (0,0) size 20 x 50
 Spaceship 2 at (0,0) size 20 x 50
 Spaceship 3 at (0,0) size 20 x 50

julia> triangular_formation!(spaceships...);

julia> spaceships
3-element Array{Widget,1}:
 Spaceship 1 at (0,0) size 20 x 50
 Spaceship 2 at (-30,-50) size 20 x 50
 Spaceship 3 at (30,-50) size 20 x 50
=#

# ------------------------------------------------------------------------------
# first class functions

function random_move()
    return rand([move_up!, move_down!, move_left!, move_right!])
end

function random_leap!(w::Widget, move_func::Function, distance::Int)
    move_func(w, distance)
    return w
end

#= REPL
julia> spaceship = Widget("Spaceship", Position(0,0), Size(20,50))
Spaceship at (0,0) size 20 x 50

julia> random_leap!(spaceship, random_move(), rand(50:100))
Spaceship at (0,82) size 20 x 50

julia> random_leap!(spaceship, random_move(), rand(50:100))
Spaceship at (0,15) size 20 x 50

julia> random_leap!(spaceship, random_move(), rand(50:100))
Spaceship at (-63,15) size 20 x 50
=#

# ------------------------------------------------------------------------------
# anonymous functions

# using named function

function explode(x)
    println(x, " exploded!")
end

function clean_up_galaxy(asteroids)
    foreach(explode, asteroids)
end

#= REPL
julia> clean_up_galaxy(asteroids)
Asteroid #1 at (10,0) size 1 x 2 exploded!
Asteroid #2 at (30,50) size 1 x 1 exploded!
Asteroid #3 at (40,60) size 5 x 3 exploded!
Asteroid #4 at (0,20) size 3 x 1 exploded!
Asteroid #5 at (80,50) size 3 x 5 exploded!
=#

# using anonymous function
function clean_up_galaxy(asteroids)
    foreach(x -> println(x, " exploded!"), asteroids)
end

# bound to variable, and it can be used more than once
function clean_up_galaxy(asteroids, spaceships)
    ep = x -> println(x, " exploded!")
    foreach(ep, asteroids)
    foreach(ep, spaceships)
end

#= REPL
julia> clean_up_galaxy(asteroids, spaceships)
Asteroid #1 at (10,0) size 1 x 2 exploded!
Asteroid #2 at (30,50) size 1 x 1 exploded!
Asteroid #3 at (40,60) size 5 x 3 exploded!
Asteroid #4 at (0,20) size 3 x 1 exploded!
Asteroid #5 at (80,50) size 3 x 5 exploded!
Spaceship 1 at (0,0) size 20 x 50 exploded!
Spaceship 2 at (-30,-50) size 20 x 50 exploded!
Spaceship 3 at (30,-50) size 20 x 50 exploded!
=#

# ------------------------------------------------------------------------------
# do syntax

# Random healthiness function for testing
healthy(spaceship) = rand(Bool)

# make sure that the spaceship is healthy before any operation
function fire(f::Function, spaceship::Widget)
    if healthy(spaceship)
        f(spaceship)
    else
        println("Operation aborted as spaceship is not healthy")
    end
    return nothing
end

# use anonymous function

# julia> fire(s -> println(s, " launched missile!"), spaceship)
# Spaceship at (78,71) size 20 x 50 launched missile!

# julia> fire(s -> println(s, " launched missile!"), spaceship)
# Operation aborted as spaceship is not healthy

# what if we have a more complex operation?
fire(s -> begin
        move_up!(s, 100)
        println(s, " launched missile!")
        move_down!(s, 100)
    end, spaceship)

# cleaner syntax
fire(spaceship) do s
    move_up!(s, 100)
    println(s, " launched missile!")
    move_down!(s, 100)
end

# open/close house keeping
function process_file(func::Function, filename::AbstractString)
    ios = nothing
    try
        ios = open(filename)
        func(ios)
    finally
        close(ios)
    end
end

# ------------------------------------------------------------------------------
# Closure

# capture and hide variable
function create_attacker(name)
    spaceship = Widget(name, Position(0,0), Size(20,50))
    return attackee -> println(spaceship.name, " attacks ", attackee.name)
end
attack = create_attacker("Goose Fighter")
attack(asteroids[1])

# # motivation
# move_up!(spaceship, 100)
# move_down!(spaceship, 100)
# move_left!(spaceship, 100)
# move_right!(spaceship, 100)

# # create clousure by capturing the spaceship variable from outer scope
# up!    = v::Real -> move_up!(spaceship, v)
# down!  = v::Real -> move_down!(spaceship, v)
# left!  = v::Real -> move_left!(spaceship, v)
# right! = v::Real -> move_right!(spaceship, v)

# # cleaner code as there is no mention of spaceship!
# up!(100)
# down!(100)
# left!(100)
# right!(100)

# # using closure to implemnt callback
# function safe_move_up!(w::Widget, callback::Function)
#     try
#         # let's say this call may fail if it collides with an asteroid
#         move_up!(w, 10)
#     catch ex
#         callback(w, ex)
#     end
# end

# # how to use...
# crash_handler(widget, ex) = println(widget, " problem: ", ex)
# safe_move_up!(spaceship, crash_handler)


text:
# Basic types and constants for inspector-executor SpBLAS API

@enum(
    sparse_status_t,
    SPARSE_STATUS_SUCCESS           = 0,
    SPARSE_STATUS_NOT_INITIALIZED   = 1,
    SPARSE_STATUS_ALLOC_FAILED      = 2,
    SPARSE_STATUS_INVALID_VALUE     = 3,
    SPARSE_STATUS_EXECUTION_FAILED  = 4,
    SPARSE_STATUS_INTERNAL_ERROR    = 5,
    SPARSE_STATUS_NOT_SUPPORTED     = 6
)

@enum(
    sparse_operation_t,
    SPARSE_OPERATION_NON_TRANSPOSE       = 10,
    SPARSE_OPERATION_TRANSPOSE           = 11,
    SPARSE_OPERATION_CONJUGATE_TRANSPOSE = 12
)

@enum(
    sparse_matrix_type_t,
    SPARSE_MATRIX_TYPE_GENERAL            = 20,
    SPARSE_MATRIX_TYPE_SYMMETRIC          = 21,
    SPARSE_MATRIX_TYPE_HERMITIAN          = 22,
    SPARSE_MATRIX_TYPE_TRIANGULAR         = 23,
    SPARSE_MATRIX_TYPE_DIAGONAL           = 24,
    SPARSE_MATRIX_TYPE_BLOCK_TRIANGULAR   = 25,
    SPARSE_MATRIX_TYPE_BLOCK_DIAGONAL     = 26
)

@enum(
    sparse_index_base_t,
    SPARSE_INDEX_BASE_ZERO    = 0,
    SPARSE_INDEX_BASE_ONE     = 1
)

@enum(
    sparse_fill_mode_t,
    SPARSE_FILL_MODE_LOWER  = 40,
    SPARSE_FILL_MODE_UPPER  = 41,
    SPARSE_FILL_MODE_FULL   = 42
)

@enum(
    sparse_diag_type_t,
    SPARSE_DIAG_NON_UNIT    = 50,
    SPARSE_DIAG_UNIT        = 51
)

@enum(
    sparse_layout_t,
    SPARSE_LAYOUT_ROW_MAJOR    = 101,
    SPARSE_LAYOUT_COLUMN_MAJOR = 102
)

@enum(
    verbose_mode_t,
    SPARSE_VERBOSE_OFF      = 70,
    SPARSE_VERBOSE_BASIC    = 71,
    SPARSE_VERBOSE_EXTENDED = 72
)

@enum(
    sparse_memory_usage_t,
    SPARSE_MEMORY_NONE          = 80,
    SPARSE_MEMORY_AGGRESSIVE    = 81
)

@enum(
    sparse_request_t,
    SPARSE_STAGE_FULL_MULT            = 90,
    SPARSE_STAGE_NNZ_COUNT            = 91,
    SPARSE_STAGE_FINALIZE_MULT        = 92,
    SPARSE_STAGE_FULL_MULT_NO_VAL     = 93,
    SPARSE_STAGE_FINALIZE_MULT_NO_VAL = 94
)

struct matrix_descr
    type::sparse_matrix_type_t
    mode::sparse_fill_mode_t
    diag::sparse_diag_type_t
end


text:
"""
    AbstractKKTSystem{T, MT<:AbstractMatrix{T}}

Abstract type for KKT system.
"""
abstract type AbstractKKTSystem{T, MT<:AbstractMatrix{T}} end

"""
    AbstractUnreducedKKTSystem{T, MT} <: AbstractKKTSystem{T, MT}

Augmented KKT system associated to the linearization of the KKT
conditions at the current primal-dual iterate ``(x, s, y, z, ν, w)``.

The associated matrix is
```
[Wₓₓ  0  Aₑ'  Aᵢ'  -I   0 ]  [Δx]
[ 0   0   0   -I    0  -I ]  [Δs]
[Aₑ   0   0    0    0   0 ]  [Δy]
[Aᵢ  -I   0    0    0   0 ]  [Δz]
[V    0   0    0    X   0 ]  [Δν]
[0    W   0    0    0   S ]  [Δw]
```
with
* ``Wₓₓ``: Hessian of the Lagrangian.
* ``Aₑ``: Jacobian of the equality constraints
* ``Aᵢ``: Jacobian of the inequality constraints
* ``X = diag(x)``
* ``S = diag(s)``
* ``V = diag(ν)``
* ``W = diag(w)``
"""
abstract type AbstractUnreducedKKTSystem{T, MT} <: AbstractKKTSystem{T, MT} end

"""
    AbstractReducedKKTSystem{T, MT} <: AbstractKKTSystem{T, MT}

The reduced KKT system is a simplification of the original Augmented KKT system.
Comparing to [`AbstractUnreducedKKTSystem`](@ref)), `AbstractReducedKKTSystem` removes
the two last rows associated to the bounds' duals ``(ν, w)``.

At a primal-dual iterate ``(x, s, y, z)``, the matrix writes
```
[Wₓₓ + Σₓ   0    Aₑ'   Aᵢ']  [Δx]
[ 0         Σₛ    0    -I ]  [Δs]
[Aₑ         0     0     0 ]  [Δy]
[Aᵢ        -I     0     0 ]  [Δz]
```
with
* ``Wₓₓ``: Hessian of the Lagrangian.
* ``Aₑ``: Jacobian of the equality constraints
* ``Aᵢ``: Jacobian of the inequality constraints
* ``Σₓ = X⁻¹ V``
* ``Σₛ = S⁻¹ W``

"""
abstract type AbstractReducedKKTSystem{T, MT} <: AbstractKKTSystem{T, MT} end

"""
    AbstractCondensedKKTSystem{T, MT} <: AbstractKKTSystem{T, MT}

The condensed KKT system simplifies further the [`AbstractReducedKKTSystem`](@ref)
by removing the rows associated to the slack variables ``s`` and the inequalities.

At the primal-dual iterate ``(x, y)``, the matrix writes
```
[Wₓₓ + Σₓ + Aᵢ' Σₛ Aᵢ    Aₑ']  [Δx]
[         Aₑ              0 ]  [Δy]
```
with
* ``Wₓₓ``: Hessian of the Lagrangian.
* ``Aₑ``: Jacobian of the equality constraints
* ``Aᵢ``: Jacobian of the inequality constraints
* ``Σₓ = X⁻¹ V``
* ``Σₛ = S⁻¹ W``
"""
abstract type AbstractCondensedKKTSystem{T, MT} <: AbstractKKTSystem{T, MT} end

#=
    Templates
=#

"Number of primal variables associated to the KKT system."
function num_variables end

"""
    get_kkt(kkt::AbstractKKTSystem)::AbstractMatrix

Return a pointer to the KKT matrix implemented in `kkt`.
The pointer is passed afterward to a linear solver.
"""
function get_kkt end

"Get Jacobian matrix"
function get_jacobian end

"Get Hessian matrix"
function get_hessian end

"""
    initialize!(kkt::AbstractKKTSystem)

Initialize KKT system with default values.
Called when we initialize the `InteriorPointSolver` storing the current KKT system `kkt`.
"""
function initialize! end

"""
    build_kkt!(kkt::AbstractKKTSystem)

Assemble the KKT matrix before calling the factorization routine.

"""
function build_kkt! end

"""
    compress_hessian!(kkt::AbstractKKTSystem)

Compress the Hessian inside `kkt`'s internals.
This function is called every time a new Hessian is evaluated.

Default implementation do nothing.

"""
function compress_hessian! end

"""
    compress_jacobian!(kkt::AbstractKKTSystem)

Compress the Jacobian inside `kkt`'s internals.
This function is called every time a new Jacobian is evaluated.

By default, the function updates in the Jacobian the coefficients
associated to the slack variables.

"""
function compress_jacobian! end

"""
    jtprod!(y::AbstractVector, kkt::AbstractKKTSystem, x::AbstractVector)

Multiply with transpose of Jacobian and store the result
in `y`, such that ``y = A' x`` (with ``A`` current Jacobian).
"""
function jtprod! end

"""
    regularize_diagonal!(kkt::AbstractKKTSystem, primal_values::AbstractVector, dual_values::AbstractVector)

Regularize the values in the diagonal of the KKT system.
Called internally inside the interior-point routine.
"""
function regularize_diagonal! end

"""
    set_jacobian_scaling!(kkt::AbstractKKTSystem, scaling::AbstractVector)

Set the scaling of the Jacobian with the vector `scaling` storing
the scaling for all the constraints in the problem.

"""
function set_jacobian_scaling! end

"""
    is_inertia_correct(kkt::AbstractKKTSystem, n::Int, m::Int, p::Int)

Check if the inertia ``(n, m, p)`` returned by the linear solver is adapted
to the KKT system implemented in `kkt`.

"""
function is_inertia_correct end

# TODO: temporary
"Return true if KKT system is reduced."
function is_reduced end

"Nonzero in Jacobian"
function nnz_jacobian end

# TODO: we need these two templates as NLPModels does not implement
# a template for dense Jacobian and dense Hessian
"Dense Jacobian callback"
function jac_dense! end

"Dense Hessian callback"
function hess_dense! end

#=
    Generic functions
=#
function initialize!(kkt::AbstractKKTSystem)
    fill!(kkt.pr_diag, 1.0)
    fill!(kkt.du_diag, 0.0)
    fill!(kkt.hess, 0.0)
end

function regularize_diagonal!(kkt::AbstractKKTSystem, primal, dual)
    kkt.pr_diag .+= primal
    kkt.du_diag .= .-dual
end

Base.size(kkt::AbstractKKTSystem) = size(kkt.aug_com)
Base.size(kkt::AbstractKKTSystem, dim::Int) = size(kkt.aug_com, dim)

# Getters
get_kkt(kkt::AbstractKKTSystem) = kkt.aug_com
get_jacobian(kkt::AbstractKKTSystem) = kkt.jac
get_hessian(kkt::AbstractKKTSystem) = kkt.hess
get_raw_jacobian(kkt::AbstractKKTSystem) = kkt.jac_raw

# Fix variable treatment
function treat_fixed_variable!(kkt::AbstractKKTSystem{T, MT}) where {T, MT<:SparseMatrixCSC{T, Int32}}
    length(kkt.ind_fixed) == 0 && return
    aug = kkt.aug_com

    fixed_aug_diag = view(aug.nzval, aug.colptr[kkt.ind_fixed])
    fixed_aug_diag .= 1.0
    fixed_aug = view(aug.nzval, kkt.ind_aug_fixed)
    fixed_aug .= 0.0
    return
end
function treat_fixed_variable!(kkt::AbstractKKTSystem{T, MT}) where {T, MT<:Matrix{T}}
    length(kkt.ind_fixed) == 0 && return
    aug = kkt.aug_com
    @inbounds for i in kkt.ind_fixed
        aug[i, :] .= 0.0
        aug[:, i] .= 0.0
        aug[i, i]  = 1.0
    end
end

function is_inertia_correct(kkt::AbstractKKTSystem, num_pos, num_zero, num_neg)
    return (num_zero == 0) && (num_pos == num_variables(kkt))
end

function build_kkt!(kkt::AbstractKKTSystem{T, MT}) where {T, MT<:Matrix{T}}
    copyto!(kkt.aug_com, kkt.aug_raw)
    treat_fixed_variable!(kkt)
end

function build_kkt!(kkt::AbstractKKTSystem{T, MT}) where {T, MT<:SparseMatrixCSC{T, Int32}}
    transfer!(kkt.aug_com, kkt.aug_raw, kkt.aug_csc_map)
    treat_fixed_variable!(kkt)
end

compress_hessian!(kkt::AbstractKKTSystem) = nothing

include("rhs.jl")
include("sparse.jl")
include("dense.jl")


text:
using LinearAlgebra
using Test

include("test_solve_eye_plus_at_kron_b.jl")


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

using Markdown
using InteractiveUtils

# ╔═╡ 8af90e7a-0e5e-4a9d-b39e-f7603fc1e25e
begin
	ENV["LC_NUMERIC"]="C" # necessary for ensuring
	                      # proper parallel use of PyPlot and Triangulate
	using Triangulate, PyPlot, PlutoUI, Printf
	PyPlot.svg(true)
end;

# ╔═╡ 7b468c8a-ce16-11eb-3009-f9e4c79e6c17
md"""
# Triangulate examples
$(TableOfContents(title="",aside=false))
## Triangulations of point sets
"""

# ╔═╡ 83f07f94-7e35-42de-9214-715c92b6dddd
md"""
### Delaunay triangulation (DT)

The Delaunay triangulation of a set of points is a triangulation
where for each triangle, the interior of its circumcircle  does not contain any points
of the triangulation.

The Delaunay triangulation of a set of points in general position
(no 4 of them on a circle) is unique. At the same time, it is a
triangulation of the convex hull of these points.

Given an input list of points, without any further flags, Triangle creates
just this triangulation (the "Q" flag suppresses the text output of Triangle).
For this and the next examples, the input list of points is created randomly,
but on a raster, preventing the appearance of too close points.

"""

# ╔═╡ c56434f3-e0a4-462e-86d5-fcfbb6778c11
md"""
### DT with boundary
Same as the previous example, but in addition specify the "c" flag.

In this case, Triangle outputs an additional list of segments
describing the boundary of the convex hull. In fact this is a constrained
Delaunay triangulation (CDT) where the boundary segments
are the seen constraining edges which must appear in the output.

"""

# ╔═╡ f13ad00e-a3d3-462f-93d4-ebe1fe4a7b4d
md"""
### DT with Voronoi diagram

Same as the previous example, but instead of "c" specify the "v" flag.
In this case, Triangle outputs information about the Voronoi diagram
of the point set which is a structure dual to  the Delaunay triangulation.

The Voronoi cell around a point $p$ a point set $S$ is defined
as the set of points $x$ such that $|x-p|<|x-q|$ for all $q\in S$
such that $p\neq q$.
The Voronoi cells of boundary points
of the convex hull are of infinite size. The corners of the Voronoi
cells are the circumcenters of the triangles. They can be far
outside of the triangulated domain.

"""

# ╔═╡ 852e1ebf-3f29-472a-945e-53920d7653be
md"""
### Boundary conforming DT (BCDT)
Specify "c" flag for convex hull segments, "v" flag for Voronoi
and "D" flag for creating a boundary conforming Delaunay triangulation of
the point set.

In this case additional points ("Steiner points") are created which split
the boundary segments and ensure that all triangle circumcenters
lie within the convex hull.
Due to random input, there may be situations where Triangle fails with this task,
so we check for the corresponding exception.

"""

# ╔═╡ b95859af-ab21-4c41-b095-82b3ea39724a
md"""
### Constrained DT (CDT)
Constrained Delaunay triangulation (CDT) of a point set with
additional constraints given a priori. This is obtained when
specifying the "p" flag and an additional list of segments each described
by two points
which should become edges of the triangulation. Note that
the resulting triangulation is not Delaunay in the sense
given above.

"""

# ╔═╡ c5a87fb5-c264-498d-aa7c-91d9e6d7bfe2
md"""
## Triangulations of domains
"""

# ╔═╡ 7f2b4252-a827-4de4-abb6-42d1bbde44b9
md"""
### CDT of a domain

Specification is similar to that of the CDT of a point set.

The domain is given by a segment list specifying its boundary.

This is obtained by specifying the "p" flag.

"""

# ╔═╡ 6793711b-d476-45c6-8b2b-b1d0693d8e6e
md"""
### CDT with maximum area constraint

This constraint is specfied as a floating
point number given after the -a flag.
Be careful to not give it in the exponential format as Triangle would be unable to analyse it.
Therefore it is dangerous to provide it via string interpolation and it is better to convert it to a string before using `@sprintf`.

Specifying only the maximum area constraint does not prevent very thin
triangles from occuring at the boundary.

"""

# ╔═╡ 3f2f0b19-3c13-482d-abd5-e6325802c6e8
md"""
### BCDT with maximum area constraint

In addition to the area constraint specify the -D flag
in order to keep the triangle circumcenters  within the domain.

"""

# ╔═╡ 14145886-2a5f-4ae8-bf00-a5eaaa7e0568
md"""
### CDT with minimum angle condition

The "q" flag  allows to specify a minimum angle
constraint preventing skinny triangles.

This combination of flags, possibly with an additional "D" flag is recommended
when creating triangulations for finite element or finite volume methods.
It the mimimum angle is larger then 28.6 degrees, Triangle's algorithm may
run into an infinite loop.

"""

# ╔═╡ 61dc30ce-9256-4f20-9163-e8d473ac9e53
md"""
### Triangulation with refinement callback

A maximum area constraint is specified in the `unsuitable` callback
which is activated via the "u" flag if it has been passed before calling triangulate.

In addition, the "q" flag  allows to specify a minimum angle
constraint preventing skinny triangles.

"""

# ╔═╡ 30990760-92e1-4a26-98a5-b43be658dee7
md"""
### Triangulation of a heterogeneous domain

The segment list specifies its boundary and the inner boundary between subdomains.
An additional region list is specified which provides "region points" in `regionlist[1,:]`
and `regionlist[2,:]`.  These kind of mark the subdomains. `regionlist[3,:]` contains an attribute
which labels the subdomains. `regionlist[4,:]` contains a maximum area value. `size(regionlist,2)`
is the number of regions.

With the "A" flag, the subdomain labels are spread to all triangles in the corresponding
subdomains, becoming available in `triangleattributelist[1,:]`.
With the "a" flag, the area constraints are applied in the corresponding subdomains.

"""

# ╔═╡ 783975f5-cbab-4b23-9155-7a6276cd25df
md"""
### Triangulation of a domain with holes
The segment list specifies its boundary and the boundaries of the holes.
An additional hole list is specified which provides "hole points" in `holelist[1,:]`
and `holelist[2,:]`.

"""

# ╔═╡ 9a090bba-093b-4ca4-a186-1c43b52cd4ff
html"""<hr>"""

# ╔═╡ 9447e874-22ce-4b99-9037-e0d202430ee2
# Wrap "Pyplotting" into this function in order to shield calling code
# from all these peculiarities.
function pyplot(f;w=650,h=300)
	PyPlot.close()
	PyPlot.clf()
	fig=PyPlot.figure(1,dpi=100)
	fig.set_size_inches(w/100,h/100,forward=true)
	f()
	PyPlot.gcf()
end;

# ╔═╡ a5f7aca5-9e40-471a-bece-34498a804bd8
function example_convex_hull(;n=10,raster=10)
    triin=Triangulate.TriangulateIO()
    triin.pointlist=hcat(unique([ Cdouble[rand(1:raster)/raster, rand(1:raster)/raster] for i in 1:n])...)
    display(triin)
    (triout, vorout)=triangulate("Q", triin)
	pyplot() do
		plot_in_out(PyPlot,triin,triout,title="Convex hull",circumcircles=true)
	end
end;

# ╔═╡ d0e63ebd-9288-42d5-9735-3c94a2baa8e3
example_convex_hull(;n=10,raster=10)

# ╔═╡ 6f3f0014-42bc-4565-a03e-208fed8b8f48
function example_convex_hull_with_boundary(;n=10,raster=10)
    triin=Triangulate.TriangulateIO()
    triin.pointlist=hcat(unique([ Cdouble[rand(1:raster)/raster, rand(1:raster)/raster] for i in 1:n])...)
    display(triin)
    (triout, vorout)=triangulate("cQ", triin)
	pyplot() do
    	plot_in_out(PyPlot,triin,triout,title="Convex hull with boundary")
	end
end;

# ╔═╡ 6ec811ca-e6d3-43ea-8da9-02bb05060d8d
example_convex_hull_with_boundary(;n=10,raster=10)

# ╔═╡ 6e2d70c4-dc0f-4672-823a-930a50114811
function example_convex_hull_voronoi(;n=10,raster=10)
    triin=Triangulate.TriangulateIO()
    triin.pointlist=hcat(unique([ Cdouble[rand(1:raster)/raster, rand(1:raster)/raster] for i in 1:n])...)
    (triout, vorout)=triangulate("vQ", triin)
	pyplot() do
    plot_in_out(PyPlot,triin,triout,voronoi=vorout,title="Convex hull with Voronoi diagram")
	end
end;

# ╔═╡ 0b54833f-0458-4417-aa03-27c4fe2a873c
example_convex_hull_voronoi(;n=10,raster=10)

# ╔═╡ 1ee0af11-96e4-4929-926e-d7143ed1f791
function example_convex_hull_voronoi_delaunay(;n=10,raster=10)
    triin=Triangulate.TriangulateIO()
    triin.pointlist=hcat(unique([ Cdouble[rand(1:raster)/raster, rand(1:raster)/raster] for i in 1:n])...)
    try
        (triout, vorout)=triangulate("vcDQ", triin)
		pyplot() do
        plot_in_out(PyPlot,triin,triout,voronoi=vorout,title="Convex hull with Voronoi diagram")
		end
    catch err
        if typeof(err)==TriangulateError
            println("Triangle had some problem.")
            return true
        end
    end
end;

# ╔═╡ 4da9a598-bf59-4b3b-aea7-ee1c2721c206
 example_convex_hull_voronoi_delaunay(;n=10,raster=10)

# ╔═╡ 0d648d22-aa03-437b-b019-59bd693bc55e
function example_cdt(;n=10,raster=10)
    triin=Triangulate.TriangulateIO()
    triin.pointlist=hcat(unique([ Cdouble[rand(1:raster)/raster, rand(1:raster)/raster] for i in 1:n])...)
    npt=size(triin.pointlist,2)
    triin.segmentlist=Matrix{Cint}([1 2; npt-1 npt-2;  1 npt;]')
    triin.segmentmarkerlist=Vector{Cint}([2,3,4])
    (triout, vorout)=triangulate("pcQ", triin)
	pyplot() do
    plot_in_out(PyPlot,triin,triout,title="CDT")
	end
end;

# ╔═╡ 910ce428-4989-4b08-8037-887dae5c847e
example_cdt(;n=10,raster=10)

# ╔═╡ ccd8ed8d-7991-4a19-8a05-65761adc2fee
function example_domain_cdt()
    triin=Triangulate.TriangulateIO()
    triin.pointlist=Matrix{Cdouble}([0.0 0.0 ; 1.0 0.0 ; 1.0  1.0 ; 0.6 0.6; 0.0 1.0]')
    triin.segmentlist=Matrix{Cint}([1 2 ; 2 3 ; 3 4 ; 4 5 ; 5 1 ]')
    triin.segmentmarkerlist=Vector{Int32}([1, 2, 3, 4, 5])
    (triout, vorout)=triangulate("pQ", triin)
	pyplot() do
    plot_in_out(PyPlot,triin,triout,title="Domain CDT")
	end
end;

# ╔═╡ fc5264cf-b1f7-43c5-bcff-09e6274ca215
example_domain_cdt()

# ╔═╡ b71b1faf-ff3f-4404-8a2c-9de84fa498f7
function example_domain_cdt_area(;maxarea=0.05)
    triin=Triangulate.TriangulateIO()
    triin.pointlist=Matrix{Cdouble}([0.0 0.0 ; 1.0 0.0 ; 1.0  1.0 ; 0.6 0.6; 0.0 1.0]')
    triin.segmentlist=Matrix{Cint}([1 2 ; 2 3 ; 3 4 ; 4 5 ; 5 1 ]')
    triin.segmentmarkerlist=Vector{Int32}([1, 2, 3, 4, 5])
    area=@sprintf("%.15f",maxarea) # Don't use exponential format!
    (triout, vorout)=triangulate("pa$(area)Q", triin)
	pyplot() do
    plot_in_out(PyPlot,triin,triout,voronoi=vorout, title="Domain CDT with area constraint",circumcircles=true)
	end
end;

# ╔═╡ cd007961-7f0e-4c56-9c86-1a9827a71e3e
example_domain_cdt_area(;maxarea=0.05)

# ╔═╡ 8ac6cc0b-e4e7-4aa0-bc37-b22b15bf2a83
function example_domain_bcdt_area(;maxarea=0.05)
    triin=Triangulate.TriangulateIO()
    triin.pointlist=Matrix{Cdouble}([0.0 0.0 ; 1.0 0.0 ; 1.0  1.0 ; 0.6 0.6; 0.0 1.0]')
    triin.segmentlist=Matrix{Cint}([1 2 ; 2 3 ; 3 4 ; 4 5 ; 5 1 ]')
    triin.segmentmarkerlist=Vector{Int32}([1, 2, 3, 4, 5])
    area=@sprintf("%.15f",maxarea)
    (triout, vorout)=triangulate("pa$(area)DQ", triin)
	pyplot() do
    plot_in_out(PyPlot,triin,triout,voronoi=vorout, title="Boundary conforming Delaunay triangulation", circumcircles=true)

	end
end;

# ╔═╡ a0cb6060-8278-444a-a4e5-46055d98616c
example_domain_bcdt_area(;maxarea=0.05)

# ╔═╡ a7e97f1c-e091-4555-abc0-a71abd22dd8a
function example_domain_qcdt_area(;minangle=20, maxarea=0.05)
    triin=Triangulate.TriangulateIO()
    triin.pointlist=Matrix{Cdouble}([0.0 0.0 ; 1.0 0.0 ; 1.0  1.0 ; 0.6 0.6; 0.0 1.0]')
    triin.segmentlist=Matrix{Cint}([1 2 ; 2 3 ; 3 4 ; 4 5 ; 5 1 ]')
    triin.segmentmarkerlist=Vector{Int32}([1, 2, 3, 4, 5])
    area=@sprintf("%.15f",maxarea)
    angle=@sprintf("%.15f",minangle)
    (triout, vorout)=triangulate("pa$(area)q$(angle)", triin)
	pyplot() do
    plot_in_out(PyPlot,triin,triout,voronoi=vorout, title="Quality triangulation")
	end
end;

# ╔═╡ 4e7dad4d-7d3c-4201-a56f-28e3df51e885
example_domain_qcdt_area(;maxarea=0.05,minangle=20)

# ╔═╡ 083b8aef-67f1-489a-8eff-b70fb6dc9da4
function example_domain_localref(;minangle=20)
    center_x=0.6
    center_y=0.6
    localdist=0.1
    function unsuitable(x1,y1,x2,y2,x3,y3,area)
        bary_x=(x1+x2+x3)/3.0
        bary_y=(y1+y2+y3)/3.0
        dx=bary_x-center_x
        dy=bary_y-center_y
        qdist=dx^2+dy^2
        qdist>1.0e-5 && area>0.1*qdist
    end

    triunsuitable(unsuitable)
    triin=Triangulate.TriangulateIO()
    triin.pointlist=Matrix{Cdouble}([0.0 0.0 ; 1.0 0.0 ; 1.0  1.0 ; 0.6 0.6; 0.0 1.0]')
    triin.segmentlist=Matrix{Cint}([1 2 ; 2 3 ; 3 4 ; 4 5 ; 5 1 ]')
    triin.segmentmarkerlist=Vector{Int32}([1, 2, 3, 4, 5])
    angle=@sprintf("%.15f",minangle)
    (triout, vorout)=triangulate("pauq$(angle)Q", triin)
	pyplot() do
    plot_in_out(PyPlot,triin,triout,voronoi=vorout, title="Quality triangulation with local refinement")
	end
end;

# ╔═╡ 3c88e5cb-5e85-4146-a6d7-3682f4b9a892
example_domain_localref(;minangle=20)

# ╔═╡ f0f9922b-70ce-4866-bc41-349a26095ade
function example_domain_regions(;minangle=20)
    triin=Triangulate.TriangulateIO()
    triin.pointlist=Matrix{Cdouble}([0.0 0.0 ;0.5 0.0; 1.0 0.0 ; 1.0  1.0 ; 0.6 0.6; 0.0 1.0]')
    triin.segmentlist=Matrix{Cint}([1 2 ; 2 3 ;3 4 ; 4 5 ; 5 6 ; 6 1 ; 2 5]')
    triin.segmentmarkerlist=Vector{Int32}([1, 2, 3, 4, 5, 6, 7])
    triin.regionlist=Matrix{Cdouble}([0.2 0.8; 0.2 0.2; 1 2 ; 0.01 0.05])
    angle=@sprintf("%.15f",minangle)
    (triout, vorout)=triangulate("paAq$(angle)Q", triin)
	pyplot() do
    plot_in_out(PyPlot,triin,triout,voronoi=vorout, title="Hetero domain triangulation")
	end
end;

# ╔═╡ 58dbd0e3-b34f-44ec-905a-488c7bbd07ca
example_domain_regions(minangle=20)

# ╔═╡ 3dfe5acd-72f7-4510-b3a0-40364f679073
function example_domain_holes(;minangle=20,maxarea=0.001)
    triin=Triangulate.TriangulateIO()
    triin.pointlist=Matrix{Cdouble}([0.0 0.0;
                                     1.0 0.0;
                                     1.0 1.0;
                                     0.0 1.0;
                                     0.2 0.2;
                                     0.3 0.2;
                                     0.3 0.3;
                                     0.2 0.3;
                                     0.6 0.6;
                                     0.7 0.6;
                                     0.7 0.7;
                                     0.6 0.7;
                                     ]')
    triin.segmentlist=Matrix{Cint}([1 2; 2 3; 3 4; 4 1; 5 6; 6 7; 7 8; 8 5;  9 10; 10 11; 11 12; 12 9;]')
    triin.segmentmarkerlist=Vector{Int32}([1, 1,1,1, 2,2,2,2, 3,3,3,3])
    triin.holelist=[0.25 0.25; 0.65 0.65;]'
    area=@sprintf("%.15f",maxarea) # Don't use exponential format!
    angle=@sprintf("%.15f",minangle)
    (triout, vorout)=triangulate("pa$(area)q$(angle)Q", triin)

	pyplot() do
    plot_in_out(PyPlot,triin,triout,voronoi=vorout, title="Domain with holes")
	end
end;

# ╔═╡ 6545d149-d4dd-4f26-8f51-5db90f6b444d
example_domain_holes(minangle=20,maxarea=0.05)

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
Triangulate = "f7e6ffb2-c36d-4f8f-a77e-16e897189344"

[compat]
PlutoUI = "~0.7.34"
PyPlot = "~2.10.0"
Triangulate = "~2.1.1"
"""

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

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

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

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

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

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

[[Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

[[Conda]]
deps = ["Downloads", "JSON", "VersionParsing"]
git-tree-sha1 = "6cdc8832ba11c7695f494c9d9a1c31e90959ce0f"
uuid = "8f4d0f93-b110-5947-807f-2305c1781a2d"
version = "1.6.0"

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

[[DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

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

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

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

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

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

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

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.4.1"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.3"

[[LaTeXStrings]]
git-tree-sha1 = "f2355693d6778a178ade15952b7ac47a4ff97996"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.3.0"

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

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

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

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

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

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

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

[[MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.9"

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

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

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

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

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

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "13468f237353112a01b2d6b32f3d0f80219944aa"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.2"

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

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

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "2cf929d64681236a2e074ffafb8d568733d2e6af"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.3"

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

[[PyCall]]
deps = ["Conda", "Dates", "Libdl", "LinearAlgebra", "MacroTools", "Serialization", "VersionParsing"]
git-tree-sha1 = "71fd4022ecd0c6d20180e23ff1b3e05a143959c2"
uuid = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
version = "1.93.0"

[[PyPlot]]
deps = ["Colors", "LaTeXStrings", "PyCall", "Sockets", "Test", "VersionParsing"]
git-tree-sha1 = "14c1b795b9d764e1784713941e787e1384268103"
uuid = "d330b81b-6aea-500a-939a-2ce795aea3ee"
version = "2.10.0"

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

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

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

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

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

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

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

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

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

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

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

[[Triangle_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "bfdd9ef1004eb9d407af935a6f36a4e0af711369"
uuid = "5639c1d2-226c-5e70-8d55-b3095415a16a"
version = "1.6.1+0"

[[Triangulate]]
deps = ["DocStringExtensions", "Libdl", "Printf", "Test", "Triangle_jll"]
git-tree-sha1 = "0b011b75202d936d2f1af6215bf3b6cce26f2b7b"
uuid = "f7e6ffb2-c36d-4f8f-a77e-16e897189344"
version = "2.1.1"

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

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

[[VersionParsing]]
git-tree-sha1 = "58d6e80b4ee071f5efd07fda82cb9fbe17200868"
uuid = "81def892-9a0e-5fdd-b105-ffc91e053289"
version = "1.3.0"

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

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

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

# ╔═╡ Cell order:
# ╟─7b468c8a-ce16-11eb-3009-f9e4c79e6c17
# ╟─83f07f94-7e35-42de-9214-715c92b6dddd
# ╠═a5f7aca5-9e40-471a-bece-34498a804bd8
# ╠═d0e63ebd-9288-42d5-9735-3c94a2baa8e3
# ╟─c56434f3-e0a4-462e-86d5-fcfbb6778c11
# ╠═6f3f0014-42bc-4565-a03e-208fed8b8f48
# ╠═6ec811ca-e6d3-43ea-8da9-02bb05060d8d
# ╟─f13ad00e-a3d3-462f-93d4-ebe1fe4a7b4d
# ╠═6e2d70c4-dc0f-4672-823a-930a50114811
# ╠═0b54833f-0458-4417-aa03-27c4fe2a873c
# ╟─852e1ebf-3f29-472a-945e-53920d7653be
# ╠═1ee0af11-96e4-4929-926e-d7143ed1f791
# ╠═4da9a598-bf59-4b3b-aea7-ee1c2721c206
# ╟─b95859af-ab21-4c41-b095-82b3ea39724a
# ╠═0d648d22-aa03-437b-b019-59bd693bc55e
# ╠═910ce428-4989-4b08-8037-887dae5c847e
# ╟─c5a87fb5-c264-498d-aa7c-91d9e6d7bfe2
# ╟─7f2b4252-a827-4de4-abb6-42d1bbde44b9
# ╠═ccd8ed8d-7991-4a19-8a05-65761adc2fee
# ╠═fc5264cf-b1f7-43c5-bcff-09e6274ca215
# ╟─6793711b-d476-45c6-8b2b-b1d0693d8e6e
# ╠═b71b1faf-ff3f-4404-8a2c-9de84fa498f7
# ╠═cd007961-7f0e-4c56-9c86-1a9827a71e3e
# ╟─3f2f0b19-3c13-482d-abd5-e6325802c6e8
# ╠═8ac6cc0b-e4e7-4aa0-bc37-b22b15bf2a83
# ╠═a0cb6060-8278-444a-a4e5-46055d98616c
# ╟─14145886-2a5f-4ae8-bf00-a5eaaa7e0568
# ╠═a7e97f1c-e091-4555-abc0-a71abd22dd8a
# ╠═4e7dad4d-7d3c-4201-a56f-28e3df51e885
# ╟─61dc30ce-9256-4f20-9163-e8d473ac9e53
# ╠═083b8aef-67f1-489a-8eff-b70fb6dc9da4
# ╠═3c88e5cb-5e85-4146-a6d7-3682f4b9a892
# ╟─30990760-92e1-4a26-98a5-b43be658dee7
# ╠═f0f9922b-70ce-4866-bc41-349a26095ade
# ╠═58dbd0e3-b34f-44ec-905a-488c7bbd07ca
# ╟─783975f5-cbab-4b23-9155-7a6276cd25df
# ╠═3dfe5acd-72f7-4510-b3a0-40364f679073
# ╠═6545d149-d4dd-4f26-8f51-5db90f6b444d
# ╟─9a090bba-093b-4ca4-a186-1c43b52cd4ff
# ╠═8af90e7a-0e5e-4a9d-b39e-f7603fc1e25e
# ╠═9447e874-22ce-4b99-9037-e0d202430ee2
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
module PhaseOne

using MatrixOptInterface
import Simplex

const MatOI = MatrixOptInterface

include("Artificial.jl")
include("Cplex.jl")

@enum(
    Method,
    ARTIFICIAL,
    CPLEX,
)

function run(lp::MatOI.LPSolverForm, Tv::Type;
    method::Method = CPLEX, kwargs...)::Simplex.SpxData
    if method == ARTIFICIAL
        return Artificial.run(lp, Tv; kwargs...)
    elseif method == CPLEX
        return Cplex.run(lp, Tv; kwargs...)
    end
end

end


text:
module OceanRobots

using NetCDF, Dates, CFTime, DataFrames

include("read_data.jl")
export drifters_ElipotEtAl16

end # module


text:
# ====================
# Extend with_logger to work with  LoggingExtras.TeeLogger
# ====================

"""
    Logging.with_logger(@nospecialize(f::Function), demuxLogger::TeeLogger)

Extends the method [`Logging.with_logger`](@ref) to work for a `LoggingExtras.TeeLogger`
containing a `SeqLogger`.

### Note
This constitutes as type piracy and should be treated with caution.
But it is necessary such that also `TeeLogger`s that contain a `SeqLogger` do flush
the log events after exiting the `with_logger` function.
"""
function Logging.with_logger(@nospecialize(f::Function), demuxLogger::TeeLogger)
    result = Base.CoreLogging.with_logstate(f, Base.CoreLogging.LogState(demuxLogger))
    flush_events(demuxLogger)
    return result
end

"""
    flush_events(teeLogger::LoggingExtras.TeeLogger)

Extend `flush_events` to a work for a `LoggingExtras.TeeLogger` containing a `SeqLogger`.
"""
function flush_events(teeLogger::LoggingExtras.TeeLogger)
    loggers = teeLogger.loggers # LoggingExtras API
    for logger in loggers
        flush_events(logger)
    end
    return nothing
end
flush_events(::Logging.AbstractLogger) = nothing

"""
    event_property!(teelogger::LoggingExtras.TeeLogger; kwargs...)

Extend `event_property!` to a work for a `LoggingExtras.TeeLogger` containing a `SeqLogger`.
"""
function event_property!(teelogger::LoggingExtras.TeeLogger; kwargs...)
    loggers = teelogger.loggers
    for logger in loggers
        event_property!(logger; kwargs...)
    end
    return nothing
end
event_property!(logger::AbstractLogger; kwargs...) = nothing


text:
a, b, c = parse.(Int, split(readline()))
n = 0
for i in a:b
    if c % i == 0
        global n += 1
    end
end
println(n)


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id: DoublePipeIncr.mso								$
#*------------------------------------------------------------------
type Details_Main
	Details_Main()=begin
		N=outers.N
		new(
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of zones",
				:Default=>2
			)),
			area(Dict{Symbol,Any}(
				:Brief=>"Total Exchange Surface Area"
			)),
			fill(power(Dict{Symbol,Any}(
				:Brief=>"Incremental Duty",
				:Default=>7000,
				:Lower=>1e-6,
				:Upper=>1e10
			)),(N)),
			power(Dict{Symbol,Any}(
				:Brief=>"Total Duty",
				:Default=>7000,
				:Lower=>1e-6,
				:Upper=>1e10
			)),
			fill(heat_trans_coeff(Dict{Symbol,Any}(
				:Brief=>"Incremental Overall Heat Transfer Coefficient Clean",
				:Default=>1,
				:Lower=>1e-6,
				:Upper=>1e10
			)),(N)),
			fill(heat_trans_coeff(Dict{Symbol,Any}(
				:Brief=>"Incremental Overall Heat Transfer Coefficient Dirty",
				:Default=>1,
				:Lower=>1e-6,
				:Upper=>1e10
			)),(N)),
			[:N,],
			[:A,:Q,:Qtotal,:Uc,:Ud,]
		)
	end
	N::DanaInteger
	A::area
	Q::Array{power}
	Qtotal::power
	Uc::Array{heat_trans_coeff}
	Ud::Array{heat_trans_coeff}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Details_Main
function atributes(in::Details_Main,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="to be documented"
	fields[:Info]="to be documented"
	drive!(fields,_)
	return fields
end
Details_Main(_::Dict{Symbol,Any})=begin
	newModel=Details_Main()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
# save symbols
# compress_then_write(b::Vector{DateTime}, io) = begin
#     meta = compress_then_write(map(x->x.instant, b), io)
#     (type = DateTime, len = meta.len, metadata = meta)
# end
#
# # load a Symbol column
# column_loader(::Type{DateTime}, io, metadata) = begin
#     strs = column_loader(meta.type, io, metadata.metadata)
#     DateTime.(strs)
# end
#
# column_loader!(_, ::Type{DateTime}, io, metadata) = begin
#     column_loader(DateTime, io, metadata)
# end
#
some_elm(::Type{DateTime}) = DateTime(2000,1,1,1,1,1)

# turns out DateTime is BLosc Compressable; no need to do anything

max_stars_repo_name: xiaodaigh/JDF

max_issues_repo_name: xiaodaigh/JDF

max_forks_repo_name: xiaodaigh/JDF


text:
mitigation_curve_df = DataFrame(matched = Float64[], EOR_price = Float64[], Shipped_Mt = Float64[], Domestic_EOR_Mt = Float64[], Domestic_other_Mt = Float64[], Avg_Cost = Float64[], MAC = Float64[])
matched_df = DataFrame(Cost_sharing = Float64[], matched = Float64[], Oil_Price_per_ton = Float64[], Carbon_Price_per_ton = Float64[], Objective_value = Float64[])

#for  unit_utilization_price_dict["EOR"] = 22.5
        # 0 : 2.5 : 25
#for cost_sharing_factor = 0.1 : .1 : .6
for cost_sharing_factor = 0.3

     for Oil_price = 292 : 73/4 : 511

         for CO2_price = 20 : 2 : 40

             for matched = 1000 : 100 : 57800

#for mitigated = 1000 : 100 : 57800

Layout = Model(with_optimizer(Gurobi.Optimizer))

@variable(Layout,Matched_Domestic[i in sources, j in domestic_sites] >= 0)
@variable(Layout,Matched_Terminal[i in sources, k in CO2_shipping_terminals] >= 0)
@variable(Layout,Matched_Overseas[k in CO2_shipping_terminals, l in overseas_sites] >= 0)

@constraint(Layout,source_distribution[i in sources],
                   sum(Matched_Domestic[i,j] for j in domestic_sites) + sum(Matched_Terminal[i,k] for k in CO2_shipping_terminals) <= param_source_size[i])
@constraint(Layout,sink_accumulation_Inland[j in domestic_sites],
                   sum(Matched_Domestic[i,j] for i in sources) <= param_sink_limit[j])
@constraint(Layout,sink_accumulation_Offshore[l in overseas_sites],
                    sum(Matched_Overseas[k,l] for k in CO2_shipping_terminals) <= param_sink_limit[l])
@constraint(Layout,terminal_in_out_balance[k in CO2_shipping_terminals],
                    sum(Matched_Terminal[i,k] for i in sources) == sum(Matched_Overseas[k,l] for l in overseas_sites))
@constraint(Layout, matched_con,
                sum(Matched_Domestic[i,j] for i in sources, j in domestic_sites) + sum(Matched_Terminal[i,k] for i in sources, k in CO2_shipping_terminals) >= matched)

#= scope 1
@objective(Layout, Min, sum(param_capture_cost[i]*Matched_Domestic[i,j] for i in sources, j in domestic_sites)
                        + sum(param_capture_cost[i]*Matched_Terminal[i,k] for i in sources, k in CO2_shipping_terminals) # Capture cost
                        + sum(param_transport_cost[i]*Matched_Domestic[i,j]*param_source_sink_dist[i,j] for i in sources, j in domestic_sites)
                        + sum(param_transport_cost[i]*Matched_Terminal[i,k]*param_source_sink_dist[i,k] for i in sources, k in CO2_shipping_terminals) # Inland tranport cost from source to sink and terminal
                        + sum(param_shipping_cost[[k,l]]*Matched_Overseas[k,l] for k in CO2_shipping_terminals, l in overseas_sites) # Shipping cost
                        + storage_cost*sum(Matched_Domestic[i,j_2] for i in sources, j_2 in domestic_non_EOR_sites)  # Inland stroage cost
                        - unit_CO2_sales_revenue*(sum(Matched_Domestic[i,j_1] for i in sources,j_1 in domestic_EOR_sites)+
                                                 sum(Matched_Overseas[k,l] for k in CO2_shipping_terminals, l in overseas_sites)))# CO2 price for domestic storage both EOR and non EOR
=#
# scope 2
@objective(Layout, Min, sum(param_capture_cost[i]*Matched_Domestic[i,j] for i in sources, j in domestic_sites)
                        + sum(param_capture_cost[i]*Matched_Terminal[i,k] for i in sources, k in CO2_shipping_terminals) # Capture cost
                        + sum(param_transport_cost[i]*Matched_Domestic[i,j]*param_source_sink_dist[i,j] for i in sources, j in domestic_sites)
                        + sum(param_transport_cost[i]*Matched_Terminal[i,k]*param_source_sink_dist[i,k] for i in sources, k in CO2_shipping_terminals) # Inland tranport cost from source to sink and terminal
                        + sum(param_shipping_cost[[k,l]]*Matched_Overseas[k,l] for k in CO2_shipping_terminals, l in overseas_sites) # Shipping cost
                        + storage_cost*sum(Matched_Domestic[i,j] for i in sources, j in domestic_sites) # Domestic storage cost
                        - Oil_price*oil_per_ton_CO2*sum(Matched_Domestic[i,j_1] for i in sources, j_1 in domestic_EOR_sites) # Domestic oil revenue
                        - cost_sharing_factor*(Oil_price*oil_per_ton_CO2-storage_cost)*sum(Matched_Terminal[i,k] for i in sources, k in CO2_shipping_terminals) # Overseas oil revenue
                        - CO2_price*(sum(Matched_Domestic[i,j] for i in sources, j in domestic_sites)
                                    + sum(Matched_Terminal[i,k] for i in sources, k in CO2_shipping_terminals))) # Carbon revenue
#

optimize!(Layout)

Shipped = value(sum(Matched_Terminal[i,k] for i in sources, k in CO2_shipping_terminals))
EOR_domestic = value(sum(Matched_Domestic[i,j_1] for i in sources, j_1 in domestic_EOR_sites))
Inland_others = value(sum(Matched_Domestic[i,j_2] for i in sources, j_2 in domestic_non_EOR_sites))
marginal_abatement_cost = dual(matched_con)

push!(matched_df,(cost_sharing_factor,matched,Oil_price,CO2_price, objective_value(Layout)))

#Scope 1
push!(mitigation_curve_df,(matched,unit_CO2_sales_revenue,Shipped, EOR_domestic, Inland_others,
    objective_value(Layout)/matched, marginal_abatement_cost))
end
end
end
end

max_stars_repo_name: CUHK-IITR-EnergyTransition/CCUS

max_issues_repo_name: CUHK-IITR-EnergyTransition/CCUS

max_forks_repo_name: CUHK-IITR-EnergyTransition/CCUS


text:
############################################################################################
"""
$(TYPEDEF)

Stores summary information for diagnostics output that is included in all Baytes kernels.

# Fields
$(TYPEDFIELDS)
"""
struct BaseDiagnostics{P}
    "Log objective target result."
    ℓobjective::Float64
    "Temperature for target result"
	temperature::Float64
    "Prediction of future data."
	prediction::P
    "Current iteration of kernel."
    iter::Int64
    function BaseDiagnostics(
        ℓobjective::S,
        temperature::T,
        prediction::P,
        iter::Integer
    ) where {S<:AbstractFloat, T<:AbstractFloat, P}
        ArgCheck.@argcheck 0.0 <= temperature <= 1.0
        ArgCheck.@argcheck 0 <= iter
        return new{P}(ℓobjective, temperature, prediction, iter)
    end
end

############################################################################################
function get_prediction(diagnostics::BaseDiagnostics)
    return diagnostics.prediction
end

############################################################################################
"""
$(SIGNATURES)
Return `ℓobjective` scaled up to temperature == 1.

# Examples
```julia
```

"""
function untemper(diagnostics::BaseDiagnostics)
    return (1.0/diagnostics.temperature) * diagnostics.ℓobjective
end

############################################################################################
# Export
export BaseDiagnostics, untemper