text:
module Formatters
    using Showoff

    function scientific(ticks::AbstractVector)
        Showoff.showoff(ticks, :scientific)
    end

    function plain(ticks::AbstractVector)
        try
            Showoff.showoff(ticks, :plain)
        catch e
            bt = Base.catch_backtrace()
            Base.showerror(stderr, e)
            Base.show_backtrace(stdout, bt)
            println("with ticks: ", ticks)
            String["-Inf", "Inf"]
        end
    end

end
using .Formatters

"""
    $(SIGNATURES)

Plots a 3-dimensional OldAxis.

## Attributes
$(ATTRIBUTES)
"""
@recipe(Axis3D) do scene

    q1 = qrotation(Vec3f(1, 0, 0), -0.5f0*pi)
    q2 = qrotation(Vec3f(0, 0, 1), 1f0*pi)
    tickrotations3d = (
        qrotation(Vec3f(0,0,1), -1.5pi),
        q2,
        qrotation(Vec3f(1, 0, 0), -0.5pi) * q2
    )
    axisnames_rotation3d = tickrotations3d
    tickalign3d = (
        (:left, :center), # x axis
        (:right, :center), # y axis
        (:right, :center), # z axis
    )
    axisnames_align3d = tickalign3d
    tick_color = RGBAf(0.5, 0.5, 0.5, 0.6)
    grid_color = RGBAf(0.5, 0.5, 0.5, 0.4)
    grid_thickness = 1
    axis_linewidth = 1.5
    gridthickness = ntuple(x-> 1f0, Val(3))
    axislinewidth = ntuple(x->1.5f0, Val(3))
    tsize = 5 # in percent
    Attributes(
        visible = true,
        showticks = (true, true, true),
        showaxis = (true, true, true),
        showgrid = (true, true, true),
        scale = Vec3f(1),
        padding = 0.1,
        inspectable = false,

        names = Attributes(
            axisnames = ("x", "y", "z"),
            textcolor = (:black, :black, :black),
            rotation = axisnames_rotation3d,
            textsize = (6.0, 6.0, 6.0),
            align = axisnames_align3d,
            font = lift(dim3, theme(scene, :font)),
            gap = 3
        ),

        ticks = Attributes(
            ranges_labels = (automatic, automatic),
            formatter = Formatters.plain,

            textcolor = (tick_color, tick_color, tick_color),

            rotation = tickrotations3d,
            textsize = (tsize, tsize, tsize),
            align = tickalign3d,
            gap = 3,
            font = lift(dim3, theme(scene, :font)),
        ),

        frame = Attributes(
            linecolor = (grid_color, grid_color, grid_color),
            linewidth = (grid_thickness, grid_thickness, grid_thickness),
            axislinewidth = (axis_linewidth, axis_linewidth, axis_linewidth),
            axiscolor = (:black, :black, :black),
        )
    )
end

isaxis(x) = false
isaxis(x::Axis3D) = true

const Limits{N} = NTuple{N, <:Tuple{<: Number, <: Number}}

function default_ticks(limits::Limits, ticks, scale_func::Function)
    default_ticks.(limits, (ticks,), scale_func)
end

default_ticks(limits::Tuple{Number, Number}, ticks, scale_func::Function) = default_ticks(limits..., ticks, scale_func)

function default_ticks(
        lmin::Number, lmax::Number,
        ticks::AbstractVector{<: Number}, scale_func::Function
    )
    scale_func.((filter(t -> lmin <= t <= lmax, ticks)))
end

function default_ticks(
        lmin::Number, lmax::Number, ::Automatic, scale_func::Function
    )
    # scale the limits
    scaled_ticks, mini, maxi = optimize_ticks(
        scale_func(lmin),
        scale_func(lmax);
        k_min = 4, # minimum number of ticks
        k_max = 8, # maximum number of ticks
    )
    length(scaled_ticks) == 1 && isnan(scaled_ticks[1]) && return [-Inf, Inf]
    scaled_ticks
end

function default_ticks(
        lmin::Number, lmax::Number, ticks::Integer, scale_func = identity
    )
    scaled_ticks, mini, maxi = optimize_ticks(
        scale_func(lmin),
        scale_func(lmax);
        k_min = ticks, # minimum number of ticks
        k_max = ticks, # maximum number of ticks
        k_ideal = ticks,
        # `strict_span = false` rewards cases where the span of the
        # chosen  ticks is not too much bigger than amin - amax:
        strict_span = false,
    )
    scaled_ticks
end

function default_ticks(x::Automatic, limits::Tuple, n)
    default_ticks(limits, n, identity)
end

function default_ticks(ticks::Tuple, limits::Tuple, n::Tuple)
    default_ticks.(ticks, (limits,), n)
end

default_ticks(ticks::Tuple, limits::Limits, n) = default_ticks.(ticks, limits, (n,))

default_ticks(ticks::Tuple, limits::Limits, n::Tuple) = default_ticks.(ticks, limits, n)

default_ticks(ticks::AbstractVector{<: Number}, limits, n) = ticks

function default_labels(x::NTuple{N, Any}, formatter::Function) where N
    default_labels.(x, formatter)
end

function default_labels(x::AbstractVector, y::AbstractVector, formatter::Function = Formatters.plain)
    default_labels.((x, y), formatter)
end

function default_labels(ticks::AbstractVector, formatter::Function = Formatters.plain)
    if applicable(formatter, ticks)
        return formatter(ticks) # takes the whole array
    elseif applicable(formatter, first(ticks))
        return formatter.(ticks)
    else
        error("Formatting function $(formatter) is neither applicable to $(typeof(ticks)) nor $(eltype(ticks)).")
    end
end

default_labels(x::Automatic, ranges, formatter) = default_labels(ranges, formatter)
default_labels(x::Tuple, ranges::Tuple, formatter) = default_labels.(x, ranges, (formatter,))
default_labels(x::Tuple, ranges, formatter) = default_labels.(x, (ranges,), (formatter,))
default_labels(x::AbstractVector{<: AbstractString}, ranges, formatter::Function) = x
default_labels(x::AbstractVector{<: AbstractString}, ranges::AbstractVector, formatter::Function) = x

function convert_arguments(::Type{<: Axis3D}, limits::Rect)
    e = (minimum(limits), maximum(limits))
    return (((e[1][1], e[2][1]), (e[1][2], e[2][2]), (e[1][3], e[2][3])),)
end

a_length(x::AbstractVector) = length(x)
a_length(x::Automatic) = x

function calculated_attributes!(::Type{<: Axis3D}, plot)
    ticks = plot.ticks
    args = (plot[1], ticks.ranges, ticks.labels, ticks.formatter)
    ticks[:ranges_labels] = lift(args...) do lims, ranges, labels, formatter
        num_ticks = labels === automatic ? automatic : a_length.(labels)
        ranges = default_ticks(ranges, lims, num_ticks)
        labels = default_labels(labels, ranges, formatter)
        (ranges, labels)
    end
    return
end

function labelposition(ranges, dim, dir, tgap, origin::StaticVector{N}) where N
    a, b = extrema(ranges[dim])
    whalf = Float32(((b - a) / 2))
    halfaxis = unit(Point{N, Float32}, dim) .* whalf

    origin .+ (halfaxis .+ (normalize(dir) * tgap))
end

_widths(x::Tuple{<: Number, <: Number}) = x[2] - x[1]
_widths(x) = Float32(maximum(x) - minimum(x))

to3tuple(x::Tuple{Any}) = (x[1], x[1], x[1])
to3tuple(x::Tuple{Any, Any}) = (x[1], x[2], x[2])
to3tuple(x::Tuple{Any, Any, Any}) = x
to3tuple(x) = ntuple(i-> x, Val(3))

function draw_axis3d(textbuffer, linebuffer, scale, limits, ranges_labels, args...)
    # make sure we extend all args to 3D
    ranges, ticklabels = ranges_labels
    args3d = to3tuple.(args)
    (
        showaxis, showticks, showgrid,
        axisnames, axisnames_color, axisnames_size, axisrotation, axisalign,
        axisnames_font, titlegap,
        gridcolors, gridthickness, axislinewidth, axiscolors,
        ttextcolor, trotation, ttextsize, talign, tfont, tgap,
        padding,
    ) = args3d # splat to names

    N = 3
    start!(textbuffer)
    start!(linebuffer)

    limit_widths = map(x-> x[2] - x[1], limits)
    # pad the drawn limits and use them as the ranges
    limits = map((lim, p)-> (lim[1] - p, lim[2] + p), limits, limit_widths .* padding)

    mini, maxi = first.(limits), last.(limits)

    origin = Point{N, Float32}(min.(mini, first.(ranges)))
    limit_widths = max.(last.(ranges), maxi) .- origin
    % = minimum(limit_widths) / 100 # percentage
    ttextsize = (%) .* ttextsize
    axisnames_size = (%) .* axisnames_size

    # index of the direction in which ticks and labels are drawn
    offset_indices = [ifelse(i != 2, mod1(i + 1, N), 1) for i in 1:N]
    # These need the real limits, not (%), to be scale-aware
    titlegap = 0.01limit_widths[offset_indices] .* titlegap
    tgap = 0.01limit_widths[offset_indices] .* tgap

    for i in 1:N
        axis_vec = unit(Point{N, Float32}, i)
        width = Float32(limit_widths[i])
        stop = origin .+ (width .* axis_vec)
        if showaxis[i]
            append!(linebuffer, [origin, stop], color = axiscolors[i], linewidth = axislinewidth[i])
        end
        if showticks[i]
            range = ranges[i]
            j = offset_indices[i]
            tickdir = unit(Vec{N, Float32}, j)
            offset2 = Float32(limit_widths[j] + tgap[i]) * tickdir
            for (j, tick) in enumerate(range)
                labels = ticklabels[i]
                if length(labels) >= j
                    str = labels[j]
                    if !isempty(str)
                        startpos = (origin .+ ((Float32(tick - origin[i]) * axis_vec)) .+ offset2)
                        push!(
                            textbuffer, str, startpos,
                            color = ttextcolor[i], rotation = trotation[i],
                            textsize = ttextsize[i], align = talign[i], font = tfont[i]
                        )
                    end
                end
            end
            if !isempty(axisnames[i])
                font = to_font(tfont[i])
                tick_widths = maximum(ticklabels[i]) do label
                    widths(text_bb(label, font, ttextsize[i]))[1]
                end / scale[j]
                pos = labelposition(ranges, i, tickdir, titlegap[i] + tick_widths, origin) .+ offset2
                push!(
                    textbuffer, to_latex(axisnames[i]), pos,
                    textsize = axisnames_size[i], color = axisnames_color[i],
                    rotation = axisrotation[i], align = axisalign[i], font = axisnames_font[i]
                )
            end
        end
        if showgrid[i]
            c = gridcolors[i]
            thickness = gridthickness[i]
            for _j = (i + 1):(i + N - 1)
                j = mod1(_j, N)
                dir = unit(Point{N, Float32}, j)
                range = ranges[j]
                for tick in range
                    offset = Float32(tick - origin[j]) * dir
                    append!(
                        linebuffer, [origin .+ offset, stop .+ offset],
                        color = c, linewidth = thickness
                    )
                end
            end
        end
        finish!(textbuffer); finish!(linebuffer)
    end
    return
end

function plot!(scene::SceneLike, ::Type{<: Axis3D}, attributes::Attributes, args...)
    axis = Axis3D(scene, attributes, args)
    # Disable any non linear transform for the axis plot!
    axis.transformation.transform_func[] = identity
    textbuffer = TextBuffer(axis, Point3, transparency = true, space = :data, inspectable = axis.inspectable)
    linebuffer = LinesegmentBuffer(axis, Point3, transparency = true, inspectable = axis.inspectable)

    tstyle, ticks, frame = to_value.(getindex.(axis, (:names, :ticks, :frame)))
    titlevals = getindex.(tstyle, (:axisnames, :textcolor, :textsize, :rotation, :align, :font, :gap))
    framevals = getindex.(frame, (:linecolor, :linewidth, :axislinewidth, :axiscolor))
    tvals = getindex.(ticks, (:textcolor, :rotation, :textsize, :align, :font, :gap))
    args = (
        getindex.(axis, (:showaxis, :showticks, :showgrid))...,
        titlevals..., framevals..., tvals..., axis.padding
    )
    map_once(
        draw_axis3d,
        Observable(textbuffer), Observable(linebuffer), scale(scene),
        axis[1], axis.ticks.ranges_labels, args...
    )
    push!(scene, axis)
    return axis
end

function axis3d!(scene::Scene, lims = data_limits(scene, isaxis); kw...)
    axis3d!(scene, Attributes(), lims; ticks = (ranges = automatic, labels = automatic), kw...)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ a2a0a0d7-b457-40bd-ae6b-da14cc557031
begin
	using PlutoUI
	bigbreak = HTML("<br>"^4)
	nothing
end

# ╔═╡ 6643612d-d1a6-4fdc-8398-e12682891eb0
using LinearAlgebra

# ╔═╡ 9426ed50-83ae-11ec-34f2-2f8e36fc8d05
md"""
# Composite types aka `structs`

If your project growths it gets more important to keep things modular.
"""

# ╔═╡ 869614f3-f243-4319-9f8f-64c3cc13eb97
md"## 1. Inmutable structs

The definition is rather easy. It needs name and a list of its variables (optionally with types)."

# ╔═╡ 70d2151b-0de1-41c2-a73a-7ee6f01a2d9f
struct Point
	x
	y
end

# ╔═╡ 24eabc0b-496b-4159-8da6-29910b85ea4f
md"There will always be a default constructor:"

# ╔═╡ 996a3423-bd55-49d7-9401-d270b6e109c1
p = Point(1.0, 2.0)

# ╔═╡ 909cc091-b321-4752-9f7a-62d1aa3a3a3d
p.x

# ╔═╡ be82d829-d6c3-4822-8f57-3f0727776817
p.y

# ╔═╡ 3528ea22-89da-4330-9b66-49fdfef87687
# p.x = 10   # does not work!

# ╔═╡ 260566f0-c909-41a2-9ebf-ddb8103635f4

# ╔═╡ 8f62d5bd-764b-49e5-9b00-05d9f90332f9
md"Two typical use cases of inmutable structs:

- Parameters for a solver which do not change during the simulation.
- Implementation of useful syntax, by methods functions for your type.
- Realisation of an interface for your data. (Construction of a mutable struct is often almost zero overhead.)
"

# ╔═╡ 1515c0ab-f7ed-4273-9b64-c810abbaf941
md"### Example (`FastDiagonal`)"

# ╔═╡ 5790a006-640d-4751-a690-b139654049c6
md"Let's implement a matrix like type, which only stores diagonal elements."

# ╔═╡ 62ee2d5b-d757-47db-9019-37622330e9bc
begin
	struct FastDiag{T}
	    data::Vector{T}
		z
	end

	import Base: *, +
	function *(A::FastDiag{T}, x::Vector{T}) where T
	    return A.data .* x
	end

	function +(A::FastDiag{T}, B::FastDiag{T}) where T
	    return FastDiag{T}( A.data + B.data )
	end

	nothing # just to surpress Pluto output
end

# ╔═╡ 8fabcd94-1c71-420c-b2b9-d10a6eca76dc
begin
	Base.@kwdef struct Foo
		x
		y::Float64 = 1.0
	end
end

# ╔═╡ bbe30b74-8e29-4b3f-bf04-e05e67c82dbf
FastDiag([1], 1)

# ╔═╡ dbd7131b-05a9-4d70-8864-d5ca9b2fe9c5
md"It is always good, if your type is `concrete`. That means Julia can predict each variables type in advance."

# ╔═╡ 3c7c5961-014d-47c6-8ffa-cd06198682ce
isconcretetype(FastDiag{Float64})

# ╔═╡ b9a4a4ac-c07f-4f73-aa48-dbfa123a267a

# ╔═╡ cb9e96e1-597a-4b0e-912a-1972c84edd5f
begin
    N = 1000
    A = [i == j ? rand() : 0.0 for i in 1:N, j in 1:N]
	A_diag = collect(diag(A))
    A_fast = FastDiag(diag(A))

    x = rand(N)
end

# ╔═╡ 4bdf13bd-a27c-4fda-9807-eb12b1132ea0

# ╔═╡ 0fbd5d9a-30dc-4e73-9bff-c088768de5ab
md"Let's first make sure that our methods work correctly:"

# ╔═╡ 4a9409d1-d04d-459c-93a2-fd47315a1c27
(A+A) * x == (A_fast+A_fast) * x

# ╔═╡ 1dc516fd-b6de-466c-8cc2-ec56301cebba
md"Now, let's measure performance:"

# ╔═╡ 49044d10-11a2-44e7-9a86-481ddfd595c4
with_terminal() do
	printstyled("Full matrix             :", color = :red)
	@time A * x

	printstyled("Fast digonal type       :", color = :green)
	@time A_fast * x

	printstyled("On the fly construction :", color = :green)
	@time FastDiag(A_diag) * x
end

# ╔═╡ 35345440-550c-4723-9aec-e08142ff0802
md"### Success 🚀

By the way, there is a build in Diagonal type `Diagonal` in the `LinearAlgebra` standard package.

It has basically the same definition as our example!"

# ╔═╡ 7c12b757-8ca5-4d2b-92b0-d204f94e6d0b

# ╔═╡ 292ee470-e286-47ea-8065-92ad2fc824d7
Diagonal(A_diag) * x

# ╔═╡ 55485242-77c4-43ec-996d-c1177ab34ecf
bigbreak

# ╔═╡ ee6fa527-4d43-4b51-baa3-2cb2fa57d81b
md"### Are inmutable structs constant values? NO!"

# ╔═╡ f491e000-1f10-443d-bf18-18a413f8cf94
begin
	data = [1,2,3]
    D = FastDiag(data)

	D.data[2] = 100

	(;data, D)
end

# ╔═╡ f8afb728-5080-41f0-a0c4-9937416493b8
begin
	struct InM
		x
	end

	mutable struct Mut
		y::InM
	end
end

# ╔═╡ c13de59d-508e-4472-a602-17d334337ca4
md"
The variables (here `data`) are inmutable.

But if the value is of an mutable type, than we can mutate the underlying data indeed!"

# ╔═╡ 2dc06114-7f29-4021-943d-e2b14f787e3d
bigbreak

# ╔═╡ 95fbd51f-62a0-4167-9478-452420d19369
md"## 2. Mutable structs"

# ╔═╡ 91185b65-573b-42de-8238-3cd046aa5b98
md"We just have to add the keyword `mutable`. That's it!"

# ╔═╡ fa237ef5-f1df-4826-8062-81e061be64a2
mutable struct SimulationState
	X::Vector{Float64}
	t::Float64
end

# ╔═╡ f01b0192-2a01-4e04-905b-3500c84c52c1
s = SimulationState( [10.0, 20.0], 0.0)

# ╔═╡ 5e55db0d-6245-47e6-ac65-6a37d66758ae
s.t = 10.0

# ╔═╡ 1ab4ba54-b86c-46ec-8852-326223fce26b
s

# ╔═╡ 72fec38d-4c24-4356-97de-71a2914b821f
bigbreak

# ╔═╡ 82209778-7d13-49fd-9812-1d534ad79a49
md"## 3. Parametric types and functions"

# ╔═╡ 94b65eff-9df8-4621-9a92-74ae0a0c8052
md"Warning: This is more an advanced feature. If you want to use it, your will have to study it a bit more!"

# ╔═╡ 5fa76a28-b1d4-4c94-a7b8-5057115caa7a
struct Unpredictable
	x
end

# ╔═╡ aa7189d4-53a5-4c54-bfdb-c8b242be2613
struct Predictable{T}
	x::T
end

# ╔═╡ 7fc21ab5-52d9-4b12-a0c9-780cf2c62c75

# ╔═╡ 2c21a7fe-3223-46d4-af30-296045c905ec
a = Unpredictable(rand([1.0, 2, "3"]))

# ╔═╡ d7a61c32-0ed9-4d04-91e1-d6849057df2d
b = Predictable(rand([1.0, 2, "3"]))

# ╔═╡ fc842dae-373f-4ec5-815e-e25193fce8fe
with_terminal() do
	println("a = ", a.x, "\t     typeof(a) = ", typeof(a))
	println("b = ", b.x, "\t     typeof(b) = ", typeof(b))
end

# ╔═╡ 9c4c826b-f1a8-43f9-a50c-455826223402
bigbreak

# ╔═╡ 8ad53b1a-55c3-4a85-9f06-0ecd01305120
square(s) = s.x^2

# ╔═╡ 4f962f52-99e4-4d2c-a16c-7b187b03d9d6
md"#### Problem: Type instabilities:"

# ╔═╡ e02cf43c-3389-4020-b2cc-57f201cee0aa
md"The output of `@code_warntype` shows whenever some types are unpredictable at compile time:"

# ╔═╡ 0bfee97a-6b2c-42d7-b095-b356f5858bc0
with_terminal() do
	@code_warntype square(a)
end

# ╔═╡ 4120837d-749b-4813-b53f-0ab22bceacee
with_terminal() do
	@code_warntype square(b)
end

# ╔═╡ 08facfc2-4b2b-4055-891b-b8ea6b0b61d0
bigbreak

# ╔═╡ 27352174-d77c-4d61-bcf4-b897e451c5f6
md"### Parametric functions

You can define always a domain of types which you would like to allow.

**Syntax:**
```julia
function fnc(x::T, y::S) where {T, S <: AbstractVector}
   #
end
```
"

# ╔═╡ e355c3b9-470d-447a-8208-cf99dfbceadc
function pseudo_dot(x::Vector{T}, y::Vector{S}) where {T<:Number, S<:Number}
	res = zero(T)
	for i in 1:min(length(x), length(y))
		res += x[i] * y[i]
	end
	return res
end

# ╔═╡ 99f36b6d-ade6-4f37-8dfb-76c967c445d9
pseudo_dot( [1, 2, 3], [1.0, 1im])

# ╔═╡ cebfa1b4-6876-41ac-b219-032f061f44ac
bigbreak

# ╔═╡ 665b754d-be3e-4982-8e03-02d1b5fbbf39
md"## 4. Abstract types

You might want to define a type hiearachy for your types!

We cannot define relations directly between our custom structs. Instead, we use abstract types to carry the hierachy information!

It is up to convention/documentation what an abstract type actually represents.
Famous example `AbstractArray`, see also [Documentation - Interfaces](https://docs.julialang.org/en/v1/manual/interfaces/#Indexing)
"

# ╔═╡ 7bf7c3c4-85cd-4d86-bbb1-933a7d5c8328

# ╔═╡ 40d2b89f-802d-441d-8053-ffd95feb22a0
md"Let's consider this example"

# ╔═╡ 7f67b94e-495d-4605-b0f4-7978a323a798
abstract type AbstractAlgorithm end

# ╔═╡ 9b214cfc-a24e-4e4d-86e3-71451f647268

# ╔═╡ 4f55481e-7884-4dd8-92ab-14a9fedd2840
md"#### Our `AbstractAlgorithm` Interface:

We assume that everything which is more specific than `AbstractAlgorithm`
has at least a function
```julia
minimise(::Function, ::AbstractAlgorithm)
```

Now, how could someone implement an algorithm? Just define the needed functions!
"

# ╔═╡ 685f46ae-e3dd-4514-bba2-393602831d3a

# ╔═╡ b03b9de7-022a-4ba3-924a-92fdbf599453
f(x::AbstractAlgorithm) = x.n_steps

# ╔═╡ 7f40cabc-8ce8-4ecc-b22e-1dc62bd8b648
struct RandomSolver <: AbstractAlgorithm

	domain::Tuple{Float64, Float64}
	n_steps::Int64

end

# ╔═╡ 7e32559f-7148-472a-bff0-2777bd0df755
random(d::Tuple{Float64,Float64}) = d[1] + (d[2]-d[1]) * rand()

# ╔═╡ f6c9269f-56a1-4e23-a7c7-93783c149e84
function minimise(f, r::RandomSolver)
	minimum( f( random(r.domain) ) for i in 1:r.n_steps )
end

# ╔═╡ 75c1c92f-1e03-4ec3-86dc-65cfb2bfe35e

# ╔═╡ d2e013da-c752-4fe8-a020-c67a5e89c397
bigbreak

# ╔═╡ 9705e4c1-1fdf-40aa-be08-b2df34715d11
md"## Exercises"

# ╔═╡ bfb72e16-8e42-430b-921a-cec2818eaf24

# ╔═╡ 5d2e637e-3f40-49d2-9f4b-59b7c1ee6afe
md"""
7.1

Define your own type, which represents a point on a circle.
Implement a add function, so that two points on a circle can be added!

It should allow operations such as
```julia
CirclePoint( π ) + CirclePoint( π ) == CirclePoint( 0.0 )
```

"""

# ╔═╡ ef533b3f-b705-461f-a590-b0a11fa05f52

# ╔═╡ 38599bbf-6e45-49d8-a117-2d1822688bed
bigbreak

# ╔═╡ 16dc8591-a6b3-4256-bbd9-093548352449
md"""
7.2

Implement your own an `AbstractAlgorithm`, call it `GridAlgorithm` which is compatible with the code below!

The algorithm should try to find the minimum of a function by evaluation of equidistant grid points in it's domain.

"""

# ╔═╡ ac3f9a7f-ecdf-486a-a3df-a7bdd3d91012
struct GridAlgorithm

end

# ╔═╡ c577d146-2b43-4820-956e-c708861e8881
# minimise( x -> x^2,   GridAlgorithm( (-10, 10),  10000) )

# ╔═╡ d638bc5e-44b5-401b-a4ce-e47552c28960

# ╔═╡ f2d11973-64b6-47b5-b34a-c4ae9405be0b
md"Imagine someone wrote code to allow solution of multidimensional functions.

*Does this `ProcutAlgorithm` still work with your implementation?*"

# ╔═╡ 97c78e41-21e5-46b0-b06f-c754143cae0c
struct ProductAlgorithm <: AbstractAlgorithm
	algs::Vector{AbstractAlgorithm}
end

# ╔═╡ 4e2e930c-8e2b-4bd3-98f1-ae35055d7b81
function minimise(f::Function, r::ProductAlgorithm)
	N = length(r.algs)
	x0 = zeros(N)

	for i in 1:N
		f_i = x -> f( [( j == i ) ? x : x0[j] for j in 1:N] )

		x_i = minimise( f_i, r.algs[i] )
		x0[i] = x_i
	end

	return x0
end

# ╔═╡ d69d58e1-a9d2-4933-a78e-a5cecc0c0c82
minimise( x -> x^2,   RandomSolver( (-10, 10),  10000) )

# ╔═╡ ced56d45-1a1d-4389-bb9d-5067e8500628
minimise( x -> x[1]^2 + x[2]^2,
			ProductAlgorithm( [RandomSolver( (-10.0, 10),  1000),
								RandomSolver( (-10.0, 10),  1000)] ))

# ╔═╡ 7383fa13-825b-4e65-af4a-3e90533ca342
bigbreak

# ╔═╡ c09a103b-d086-4718-90d5-9cd59d972286
md"""
7.3 Optional: (Object oriented programming)

Let's imagine three datatypes
"""

# ╔═╡ af24c0dd-b37f-420a-83bd-59e2ccddf17c
abstract type AbstractScientist end

# ╔═╡ fe617e26-b9cb-4252-bb14-45c1514e6c0f
struct Scientist <: AbstractScientist
	name
	field
end

# ╔═╡ e4a7b64f-8c2a-4077-aefb-8e8b9f22f65b
struct Mathematician <: AbstractScientist
	name
	field
	specialisation
end

# ╔═╡ a2fdbb49-dc6c-4ca7-affa-1639b62b7e10
md"""
Task: Complete the code below, such that scientists say "Hi, my name is Anna, I work in the chemisty department".

But mathematicians just say, "Hi, I'm Axel. Do you like optimisation."
"""

# ╔═╡ 7fd35958-c367-4002-9cea-482e0d1f495a
begin
	# your solution
	sayhi(x) = "Hello" # define the function appropriately

end

# ╔═╡ c106b584-471d-4ef7-b14d-ea2c82b9385b
with_terminal() do
	s = Scientist("Anna", "chemistry")
	m = Mathematician("Axel", "mathematics", "optimisation")

	println( sayhi(s) )
	println( sayhi(m) )
	nothing
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
PlutoUI = "~0.7.32"
"""

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

julia_version = "1.7.1"
manifest_format = "2.0"

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

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

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

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

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

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

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

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

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

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

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

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

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

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

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

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

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

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

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

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

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

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

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

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

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

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

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "0b5cfbb704034b5b4c1869e36634438a047df065"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.1"

[[deps.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"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

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

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

# ╔═╡ Cell order:
# ╟─a2a0a0d7-b457-40bd-ae6b-da14cc557031
# ╟─9426ed50-83ae-11ec-34f2-2f8e36fc8d05
# ╟─869614f3-f243-4319-9f8f-64c3cc13eb97
# ╠═70d2151b-0de1-41c2-a73a-7ee6f01a2d9f
# ╟─24eabc0b-496b-4159-8da6-29910b85ea4f
# ╠═996a3423-bd55-49d7-9401-d270b6e109c1
# ╠═909cc091-b321-4752-9f7a-62d1aa3a3a3d
# ╠═be82d829-d6c3-4822-8f57-3f0727776817
# ╠═3528ea22-89da-4330-9b66-49fdfef87687
# ╟─260566f0-c909-41a2-9ebf-ddb8103635f4
# ╟─8f62d5bd-764b-49e5-9b00-05d9f90332f9
# ╟─1515c0ab-f7ed-4273-9b64-c810abbaf941
# ╟─5790a006-640d-4751-a690-b139654049c6
# ╠═62ee2d5b-d757-47db-9019-37622330e9bc
# ╠═8fabcd94-1c71-420c-b2b9-d10a6eca76dc
# ╠═bbe30b74-8e29-4b3f-bf04-e05e67c82dbf
# ╟─dbd7131b-05a9-4d70-8864-d5ca9b2fe9c5
# ╠═3c7c5961-014d-47c6-8ffa-cd06198682ce
# ╟─b9a4a4ac-c07f-4f73-aa48-dbfa123a267a
# ╠═6643612d-d1a6-4fdc-8398-e12682891eb0
# ╠═cb9e96e1-597a-4b0e-912a-1972c84edd5f
# ╠═4bdf13bd-a27c-4fda-9807-eb12b1132ea0
# ╟─0fbd5d9a-30dc-4e73-9bff-c088768de5ab
# ╠═4a9409d1-d04d-459c-93a2-fd47315a1c27
# ╟─1dc516fd-b6de-466c-8cc2-ec56301cebba
# ╠═49044d10-11a2-44e7-9a86-481ddfd595c4
# ╟─35345440-550c-4723-9aec-e08142ff0802
# ╟─7c12b757-8ca5-4d2b-92b0-d204f94e6d0b
# ╠═292ee470-e286-47ea-8065-92ad2fc824d7
# ╟─55485242-77c4-43ec-996d-c1177ab34ecf
# ╟─ee6fa527-4d43-4b51-baa3-2cb2fa57d81b
# ╠═f491e000-1f10-443d-bf18-18a413f8cf94
# ╠═f8afb728-5080-41f0-a0c4-9937416493b8
# ╟─c13de59d-508e-4472-a602-17d334337ca4
# ╟─2dc06114-7f29-4021-943d-e2b14f787e3d
# ╟─95fbd51f-62a0-4167-9478-452420d19369
# ╟─91185b65-573b-42de-8238-3cd046aa5b98
# ╠═fa237ef5-f1df-4826-8062-81e061be64a2
# ╠═f01b0192-2a01-4e04-905b-3500c84c52c1
# ╠═5e55db0d-6245-47e6-ac65-6a37d66758ae
# ╠═1ab4ba54-b86c-46ec-8852-326223fce26b
# ╟─72fec38d-4c24-4356-97de-71a2914b821f
# ╟─82209778-7d13-49fd-9812-1d534ad79a49
# ╟─94b65eff-9df8-4621-9a92-74ae0a0c8052
# ╠═5fa76a28-b1d4-4c94-a7b8-5057115caa7a
# ╠═aa7189d4-53a5-4c54-bfdb-c8b242be2613
# ╟─7fc21ab5-52d9-4b12-a0c9-780cf2c62c75
# ╠═2c21a7fe-3223-46d4-af30-296045c905ec
# ╠═d7a61c32-0ed9-4d04-91e1-d6849057df2d
# ╠═fc842dae-373f-4ec5-815e-e25193fce8fe
# ╟─9c4c826b-f1a8-43f9-a50c-455826223402
# ╠═8ad53b1a-55c3-4a85-9f06-0ecd01305120
# ╟─4f962f52-99e4-4d2c-a16c-7b187b03d9d6
# ╟─e02cf43c-3389-4020-b2cc-57f201cee0aa
# ╠═0bfee97a-6b2c-42d7-b095-b356f5858bc0
# ╠═4120837d-749b-4813-b53f-0ab22bceacee
# ╟─08facfc2-4b2b-4055-891b-b8ea6b0b61d0
# ╟─27352174-d77c-4d61-bcf4-b897e451c5f6
# ╠═e355c3b9-470d-447a-8208-cf99dfbceadc
# ╠═99f36b6d-ade6-4f37-8dfb-76c967c445d9
# ╟─cebfa1b4-6876-41ac-b219-032f061f44ac
# ╟─665b754d-be3e-4982-8e03-02d1b5fbbf39
# ╟─7bf7c3c4-85cd-4d86-bbb1-933a7d5c8328
# ╟─40d2b89f-802d-441d-8053-ffd95feb22a0
# ╠═7f67b94e-495d-4605-b0f4-7978a323a798
# ╟─9b214cfc-a24e-4e4d-86e3-71451f647268
# ╟─4f55481e-7884-4dd8-92ab-14a9fedd2840
# ╟─685f46ae-e3dd-4514-bba2-393602831d3a
# ╠═b03b9de7-022a-4ba3-924a-92fdbf599453
# ╠═7f40cabc-8ce8-4ecc-b22e-1dc62bd8b648
# ╠═7e32559f-7148-472a-bff0-2777bd0df755
# ╠═f6c9269f-56a1-4e23-a7c7-93783c149e84
# ╟─75c1c92f-1e03-4ec3-86dc-65cfb2bfe35e
# ╠═d69d58e1-a9d2-4933-a78e-a5cecc0c0c82
# ╟─d2e013da-c752-4fe8-a020-c67a5e89c397
# ╟─9705e4c1-1fdf-40aa-be08-b2df34715d11
# ╟─bfb72e16-8e42-430b-921a-cec2818eaf24
# ╟─5d2e637e-3f40-49d2-9f4b-59b7c1ee6afe
# ╠═ef533b3f-b705-461f-a590-b0a11fa05f52
# ╟─38599bbf-6e45-49d8-a117-2d1822688bed
# ╟─16dc8591-a6b3-4256-bbd9-093548352449
# ╠═ac3f9a7f-ecdf-486a-a3df-a7bdd3d91012
# ╠═c577d146-2b43-4820-956e-c708861e8881
# ╟─d638bc5e-44b5-401b-a4ce-e47552c28960
# ╟─f2d11973-64b6-47b5-b34a-c4ae9405be0b
# ╠═97c78e41-21e5-46b0-b06f-c754143cae0c
# ╠═4e2e930c-8e2b-4bd3-98f1-ae35055d7b81
# ╠═ced56d45-1a1d-4389-bb9d-5067e8500628
# ╟─7383fa13-825b-4e65-af4a-3e90533ca342
# ╟─c09a103b-d086-4718-90d5-9cd59d972286
# ╠═af24c0dd-b37f-420a-83bd-59e2ccddf17c
# ╠═fe617e26-b9cb-4252-bb14-45c1514e6c0f
# ╠═e4a7b64f-8c2a-4077-aefb-8e8b9f22f65b
# ╟─a2fdbb49-dc6c-4ca7-affa-1639b62b7e10
# ╟─7fd35958-c367-4002-9cea-482e0d1f495a
# ╠═c106b584-471d-4ef7-b14d-ea2c82b9385b
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

ext: jl

max_stars_repo_path: slides/8 - Composite types.jl

max_stars_repo_name: SteffenPL/Julia-for-mathematicians

max_issues_repo_path: slides/8 - Composite types.jl

max_issues_repo_name: SteffenPL/Julia-for-mathematicians

max_forks_repo_path: slides/8 - Composite types.jl

max_forks_repo_name: SteffenPL/Julia-for-mathematicians


text:
@testset "execution" begin

function tester(A)
    idx = thread_position_in_grid_1d()
    A[idx] = Int(5)
    return nothing
end

bufferSize = 8
bufferA = MtlArray{Int,1}(undef, tuple(bufferSize), storage=Shared)
vecA = unsafe_wrap(Vector{Int}, bufferA.buffer, tuple(bufferSize))

@testset "basic execution and synchronization" begin
    @metal threads=(bufferSize) tester(bufferA.buffer)
    synchronize()
    @test all(vecA .== Int(5))
end

@testset "device synchronization" begin
    t = @async begin
        @metal threads=(bufferSize) tester(bufferA.buffer)
    end
    wait(t)
    device_synchronize()
    @test all(vecA .== Int(5))
end

@testset "launch params" begin
    vecA .= 0
    @metal threads=(2) tester(bufferA.buffer)
    synchronize()
    @test all(vecA == Int.([5, 5, 0, 0, 0, 0, 0, 0]))
    vecA .= 0

    @metal grid=(3) threads=(2) tester(bufferA.buffer)
    synchronize()
    @test all(vecA == Int.([5, 5, 5, 5, 5, 5, 0, 0]))
    vecA .= 0

    @test_throws InexactError @metal threads=(-2) tester(bufferA.buffer)
    @test_throws InexactError @metal grid=(-2) tester(bufferA.buffer)
    @test_throws ArgumentError @metal threads=(1025) tester(bufferA.buffer)
    @test_throws ArgumentError @metal threads=(1000,2) tester(bufferA.buffer)
end

@testset "argument passing" begin
    @testset "buffer argument" begin
        function kernel(ptr)
            unsafe_store!(ptr, 42)
            return
        end

        a = MtlArray([1])
        @metal kernel(pointer(a))
        @test Array(a)[] == 42
    end

    @testset "scalar argument" begin
        function kernel(ptr, val)
            unsafe_store!(ptr, val)
            return
        end

        a = MtlArray([1])
        @metal kernel(pointer(a), 42)
        @test Array(a)[] == 42
    end

    @testset "array argument" begin
        function kernel(ptr, vals)
            unsafe_store!(ptr, vals[1])
            return
        end

        a = MtlArray([1])
        @metal kernel(pointer(a), (42,))
        @test Array(a)[] == 42
    end

    @testset "struct argument" begin
        function kernel(ptr, vals)
            unsafe_store!(ptr, vals[1] + vals[2])
            return
        end

        a = MtlArray([1])
        @metal kernel(pointer(a), (20, Int32(22)))
        @test Array(a)[] == 42
    end

    @testset "indirect struct argument" begin
        function kernel(obj)
            unsafe_store!(obj[1], obj[2])
            return
        end

        a = MtlArray([1])
        @metal kernel((pointer(a), 42))
        @test Array(a)[] == 42
    end

    @testset "nested indirect struct argument" begin
        function kernel(obj)
            unsafe_store!(obj[1][1], obj[2])
            return
        end

        a = MtlArray([1])
        @metal kernel(((pointer(a), 0), 42))
        @test Array(a)[] == 42
    end

    @testset "array in struct argument" begin
        function kernel(obj)
            unsafe_store!(obj[1], obj[2][1]+obj[2][2])
            return
        end

        a = MtlArray([1])
        @metal kernel((pointer(a), (20,22)))
        @test Array(a)[] == 42
    end
end

end


text:
module ScientificTest

f(x) = 4x
g(x,y) = 2y + x

f(x::Int64) = 2

f(x::Number) = _f(x) + 2
_f(x) = 2+x

include("utils.jl")

export f

end # module


text:
function packtablegen(lengths, rarity_pct)
    weights = rarity_pct ./ lengths
    ans = Vector{ Tuple{Int,Int} }(undef, sum(lengths) )
    wvec = zeros(Float64, sum(lengths) )
    i = 0
    for (rarity,len) in enumerate( lengths)
        for j in 1:len
            i += 1
            ans[i] = (rarity, j )
            wvec[i] = weights[rarity]
        end
    end
    ans, StatsBase.weights(wvec)
end

max_stars_repo_name: asieev/HearthstoneSim

max_issues_repo_name: asieev/HearthstoneSim

max_forks_repo_name: asieev/HearthstoneSim


text:
struct WeekDayEffectProblem{Y,N,I,W,L,IL,P} <: SurveyProblem
    y::Y
    n::N
    indices::I
    weekday::W
    lnk::L
    invlnk::IL
    prior::P
end

function WeekDayEffectProblem(y,n,indices, weekday, lnk, invlnk)
    prior = PriorWeekDayEffectProblem(n, invlnk)
    WeekDayEffectProblem(y,n,indices, weekday, lnk, invlnk, prior)
end

struct PriorWeekDayEffectProblem{P,L,F,S,W,N,I}
    prior_p::P
    prior_μ0::L
    prior_ϕ::F
    prior_σ::S
    prior_weekday_effect::W
    n::N
    invlnk::I
end

function PriorWeekDayEffectProblem(n, invlnk)
    prior_p  = Beta()
    prior_μ0 = truncated(Normal(5,10),0,Inf)
    prior_ϕ  = GammaMeanCv(20,.5)
    prior_σ  = truncated(Normal(0.1, .05), 0., Inf)
    prior_weekday_effect = Dirichlet(7, 1)
    PriorWeekDayEffectProblem( prior_p, prior_μ0, prior_ϕ, prior_σ, prior_weekday_effect, n, invlnk)
end

function Base.rand(priors::PriorWeekDayEffectProblem)
    @unpack prior_p, prior_μ0, prior_ϕ, prior_σ, n, invlnk, prior_weekday_effect = priors
    p        = rand(prior_p)
    μ0       = rand(prior_μ0)
    ϕ        = rand(prior_ϕ)
    σ        = rand(prior_σ)
    latent_μ = rand( CovidSurvey.RandomWalk(n-1, σ, invlnk(μ0)) )
    weekday_effect = rand(prior_weekday_effect)
    return (; p, μ0, ϕ, σ, latent_μ, weekday_effect)
end

function Distributions.logpdf(prior::PriorWeekDayEffectProblem, θ)
    @unpack p, μ0, ϕ, σ, latent_μ, weekday_effect = θ
    @unpack prior_p, prior_μ0, prior_ϕ, prior_σ, n, invlnk, prior_weekday_effect = prior
    ℓ  = logpdf(prior_p, p)
    ℓ += logpdf(prior_μ0, μ0)
    ℓ += logpdf(prior_ϕ, ϕ)
    ℓ += logpdf(prior_σ, σ)
    ℓ += logpdf(CovidSurvey.RandomWalk(n-1, σ, invlnk(μ0)), latent_μ)
    ℓ += logpdf(prior_weekday_effect, weekday_effect)
    return ℓ
end

function lnk!(μs, problem::WeekDayEffectProblem, θ)
    @unpack lnk, weekday = problem
    @unpack μ0, latent_μ, weekday_effect = θ

    μs[1] = μ0 * weekday_effect[first(weekday)]
    for (i, x) in enumerate(latent_μ)
        μs[i+1] = lnk(x) * weekday_effect[weekday[i+1]] * 7
    end
    return nothing
end

function (problem::WeekDayEffectProblem)(θ)
    @unpack p, ϕ = θ               # extract the parameters
    @unpack y, n, indices, prior = problem
    !(ϕ > 0) && (@error "not ϕ>0"; return -Inf)
    !(0 <= p <= 1) && (@error "not 0 <= p <= 1"; return -Inf)

    ℓ  = logpdf(prior, θ)

    μ = Vector{typeof(ϕ)}(undef, n)
    lnk!(μ, problem, θ)
    !all(0 .<= μ ) && (@error "not all(0 .<= μ )"; return -Inf)

    d = likelihood_dist(p, μ, ϕ)
    ℓ += logpdf(arraydist(d[indices]),y)
    return ℓ
end

function StatsBase.predict(problem::WeekDayEffectProblem, θ)
    @unpack p, ϕ = θ               # extract the parameters
    @unpack y, n, indices = problem

    μ = Vector{typeof(ϕ)}(undef, n)
    lnk!(μ, problem, θ)

    d = likelihood_dist(p, μ, ϕ)
    y = rand(arraydist(d[indices]))
    return y
end

function problem_transformation(p::WeekDayEffectProblem)
    as((
        p = as𝕀,
        μ0 = as(Real, 0, 20),
        ϕ = asℝ₊,
        σ = asℝ₊,
        latent_μ = as(Array, p.n-1),
        weekday_effect = TransformVariables.UnitSimplex(7) #QUESTION
    ))
end

function generate_quantities(problem::WeekDayEffectProblem, θ::NamedTuple)
    @unpack n = problem

    μ = zeros(n)
    lnk!(μ, problem, θ)
    return μ
end

function WeekDayEffectProblem(data, col)
    y = data[:,col]
    indices = date2index( data.date )
    weekday = dayofweek.(data.date)
    n = last(indices)
    lnk = KLogistic(50.)
    invlnk = KLogit(50.)
    return WeekDayEffectProblem(y, n, indices, weekday, lnk, invlnk)
end

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
#=
Create caches for terms in the equations here.
=#

abstract type AbstractPDECache end

struct ConvectiveDerivativeCache{VT,DVT} <: AbstractPDECache
  V :: VT
  Vtf :: DVT
  DVf :: DVT
  VDVf :: DVT
end

struct RHSCache{VT} <: AbstractPDECache
  dv :: VT
end

struct PressureCache{VT,PT} <: AbstractPDECache
  dvp :: VT
  dvf :: VT
  lp :: PT
end

struct VelocityCache{WT,VT,FT,VBT,SBT} <: AbstractPDECache
  Wn :: WT
  Vf :: VT
  Sc :: FT
  Δus :: VBT
  Sb :: SBT
end

struct StreamfunctionCache{ST} <: AbstractPDECache
  Sn :: ST
end

struct ConstraintOperatorCache{VT,VBT} <: AbstractPDECache
  Vv :: VT
  Δus :: VBT
end

struct DoubleLayerCache{VT,VBT} <: AbstractPDECache
  Vv :: VT
  Δus :: VBT
end

struct SurfaceVelocityCache{VBT} <: AbstractPDECache
  Vb :: VBT
end

struct SurfaceTractionCache{VBT,SBT} <: AbstractPDECache
  Vb :: VBT
  Δus :: VBT
  Sb :: SBT
end


text:
using InstantiateFromURL
github_project("QuantEcon/quantecon-notebooks-julia", version = "0.5.0")
# github_project("QuantEcon/quantecon-notebooks-julia", version = "0.5.0", instantiate = true) # uncomment to force package installation

max_stars_repo_name: Shunsuke-Hori/lecture-source-jl

max_issues_repo_name: Shunsuke-Hori/lecture-source-jl

max_forks_repo_name: Shunsuke-Hori/lecture-source-jl


text:
using PSCOPF

using Test
using JuMP
using Dates
using DataStructures

@testset verbose=true "test_tso_unit_priority" begin

    #=
    TSOOutFO tries to deviate the least from the reference market schedule.
    Then, while respecting this deviation, chooses the cheapest units.

    TS: [11h]
    S: [S1]
                        bus 1
                        |
    (limitable) wind_1_1|load_2
    Pmin=0, Pmax=100    | S1: 55
    Csta=0, Cprop=150   |
      S1: 20            |
                        |
    (limitable) wind_1_2|
    Pmin=0, Pmax=100    |
    Csta=0, Cprop=1     |
      S1: 25            |
                        |
    (pilotable) prod_1_1|
    Pmin=0, Pmax=100    |
    Csta=0, Cprop=10    |
                        |
    (pilotable) prod_1_2|
     Pmin=0, Pmax=100   |
     Csta=0, Cprop=15   |
                        |
    =#
    @testset "tso_deviates_the_least_from_market" begin
        TS = [DateTime("2015-01-01T11:00:00")]
        ech = DateTime("2015-01-01T07:00:00")
        network = PSCOPF.Networks.Network()
        # Buses
        PSCOPF.Networks.add_new_bus!(network, "bus_1")
        # Limitables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "wind_1_1", PSCOPF.Networks.LIMITABLE,
                                                0., 100.,
                                                0., 150.,
                                                Dates.Second(0), Dates.Second(0))
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "wind_1_2", PSCOPF.Networks.LIMITABLE,
                                                0., 100.,
                                                0., 1.,
                                                Dates.Second(0), Dates.Second(0))
        # Pilotables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "prod_1_1", PSCOPF.Networks.PILOTABLE,
                                                0., 100.,
                                                0., 10.,
                                                Dates.Second(0), Dates.Second(0))
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "prod_1_2", PSCOPF.Networks.PILOTABLE,
                                                0., 100.,
                                                0., 15.,
                                                Dates.Second(0), Dates.Second(0))
        # Uncertainties
        uncertainties = PSCOPF.Uncertainties()
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S1", 20.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_2", DateTime("2015-01-01T11:00:00"), "S1", 25.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S1", 55.)
        # firmness
        firmness = PSCOPF.Firmness(
                    SortedDict{String, SortedDict{Dates.DateTime, PSCOPF.DecisionFirmness} }(),
                    SortedDict("wind_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE),
                                "wind_1_2" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE),
                                "prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE),
                                "prod_1_2" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE), )
                    )
        # initial generators state : No need because all pmin=0 => ON by default
        generators_init_state = SortedDict{String, PSCOPF.GeneratorState}()

        @testset "tso_follows_market_even_if_its_expensive" begin
            #=
            Market uses the more expensive unit prod_1_2
            =#
            context = PSCOPF.PSCOPFContext(network, TS, PSCOPF.PSCOPF_MODE_1,
                                            generators_init_state,
                                            uncertainties, nothing)
            context.market_schedule = PSCOPF.Schedule(PSCOPF.Market(), Dates.DateTime("2015-01-01T10:45:00"), SortedDict(
                    "wind_1_1" => PSCOPF.GeneratorSchedule("wind_1_1",
                        SortedDict{Dates.DateTime, PSCOPF.UncertainValue{PSCOPF.GeneratorState}}(),
                        SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.UncertainValue{Float64}(missing,
                                                                                SortedDict("S1"=>20.)))
                        ),
                    "wind_1_2" => PSCOPF.GeneratorSchedule("wind_1_2",
                        SortedDict{Dates.DateTime, PSCOPF.UncertainValue{PSCOPF.GeneratorState}}(),
                        SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.UncertainValue{Float64}(missing,
                                                                                SortedDict("S1"=>25.)))
                    ),
                    "prod_1_1" => PSCOPF.GeneratorSchedule("prod_1_1",
                        SortedDict{Dates.DateTime, PSCOPF.UncertainValue{PSCOPF.GeneratorState}}(),
                        SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.UncertainValue{Float64}(missing,
                                                                                SortedDict("S1"=>0.)))
                    ),
                    "prod_1_2" => PSCOPF.GeneratorSchedule("prod_1_2",
                        SortedDict{Dates.DateTime, PSCOPF.UncertainValue{PSCOPF.GeneratorState}}(),
                        SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.UncertainValue{Float64}(missing,
                                                                                SortedDict("S1"=>10.)))
                    ),
                    )
            )

            tso = PSCOPF.TSOOutFO()
            result = PSCOPF.run(tso, ech, firmness,
                        PSCOPF.get_target_timepoints(context),
                        context)
            PSCOPF.update_tso_schedule!(context, ech, result, firmness, tso)

            # Solution is optimal
            @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL

            @test value(result.objective_model.deltas) < 1e-09 # followed the market
            @test value(result.lol_model.p_loss_of_load["bus_1", TS[1], "S1"]) < 1e-09

            @test 20. ≈ PSCOPF.get_prod_value(context.tso_schedule, "wind_1_1", TS[1], "S1")
            @test 25. ≈ PSCOPF.get_prod_value(context.tso_schedule, "wind_1_2", TS[1], "S1")
            @test PSCOPF.get_prod_value(context.tso_schedule, "prod_1_1", TS[1], "S1") < 1e-09
            @test 10. ≈ PSCOPF.get_prod_value(context.tso_schedule, "prod_1_2", TS[1], "S1")

            # pmin = 0 for prod_1_1 & prod_2_1
            @test ismissing( PSCOPF.get_commitment_value(context.tso_schedule, "prod_1_1", TS[1], "S1") )
            @test ismissing( PSCOPF.get_commitment_value(context.tso_schedule, "prod_1_2", TS[1], "S1") )
        end

        @testset "tso_follows_market_while_choosing_the_cheapest" begin
            #=
            Market uses the more expensive unit prod_1_2
            5MW are missing => To respect EOD, TSO will increase prod by 5
            => deviation from market
            TSO will use prod_1_1 & prod_1_2 cause that is cheaper
            =#
            context = PSCOPF.PSCOPFContext(network, TS, PSCOPF.PSCOPF_MODE_1,
                                            generators_init_state,
                                            uncertainties, nothing)
            context.market_schedule = PSCOPF.Schedule(PSCOPF.Market(), Dates.DateTime("2015-01-01T10:45:00"), SortedDict(
                    "wind_1_1" => PSCOPF.GeneratorSchedule("wind_1_1",
                        SortedDict{Dates.DateTime, PSCOPF.UncertainValue{PSCOPF.GeneratorState}}(),
                        SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.UncertainValue{Float64}(missing,
                                                                                SortedDict("S1"=>20.)))
                        ),
                    "wind_1_2" => PSCOPF.GeneratorSchedule("wind_1_2",
                        SortedDict{Dates.DateTime, PSCOPF.UncertainValue{PSCOPF.GeneratorState}}(),
                        SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.UncertainValue{Float64}(missing,
                                                                                SortedDict("S1"=>25.)))
                    ),
                    "prod_1_1" => PSCOPF.GeneratorSchedule("prod_1_1",
                        SortedDict{Dates.DateTime, PSCOPF.UncertainValue{PSCOPF.GeneratorState}}(),
                        SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.UncertainValue{Float64}(missing,
                                                                                SortedDict("S1"=>0.)))
                    ),
                    "prod_1_2" => PSCOPF.GeneratorSchedule("prod_1_2",
                        SortedDict{Dates.DateTime, PSCOPF.UncertainValue{PSCOPF.GeneratorState}}(),
                        SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.UncertainValue{Float64}(missing,
                                                                                SortedDict("S1"=>5.)))
                    ),
                    )
            )

            tso = PSCOPF.TSOOutFO()
            result = PSCOPF.run(tso, ech, firmness,
                        PSCOPF.get_target_timepoints(context),
                        context)
            PSCOPF.update_tso_schedule!(context, ech, result, firmness, tso)

            # Solution is optimal
            @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL

            @test 5. ≈ value(result.objective_model.deltas)
            @test value(result.lol_model.p_loss_of_load["bus_1", TS[1], "S1"]) < 1e-09

            @test 20. ≈ PSCOPF.get_prod_value(context.tso_schedule, "wind_1_1", TS[1], "S1")
            @test 25. ≈ PSCOPF.get_prod_value(context.tso_schedule, "wind_1_2", TS[1], "S1")
            @test 5. ≈ PSCOPF.get_prod_value(context.tso_schedule, "prod_1_1", TS[1], "S1")
            @test 5. ≈ PSCOPF.get_prod_value(context.tso_schedule, "prod_1_2", TS[1], "S1")

            # pmin = 0 for prod_1_1 & prod_2_1
            @test ismissing( PSCOPF.get_commitment_value(context.tso_schedule, "prod_1_1", TS[1], "S1") )
            @test ismissing( PSCOPF.get_commitment_value(context.tso_schedule, "prod_1_2", TS[1], "S1") )
        end

        @testset "tso_empty_market_schedule" begin
            #=
            If market_schedule is empty
            consider that reference is 0.
            Deviate the least from market (not important in this case since cutting conso is highly penalized)
            Choose cheapest alternative
            NOTE : unlike market, limitables with high Cprop will be chosen first (cause we pay the non-usage ie capping)
            =#
            context = PSCOPF.PSCOPFContext(network, TS, PSCOPF.PSCOPF_MODE_1,
                                            generators_init_state,
                                            uncertainties, nothing)

            tso = PSCOPF.TSOOutFO()
            result = PSCOPF.run(tso, ech, firmness,
                        PSCOPF.get_target_timepoints(context),
                        context)
            PSCOPF.update_tso_schedule!(context, ech, result, firmness, tso)

            # Solution is optimal
            @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL

            @test 55. ≈ value(result.objective_model.deltas)
            @test value(result.lol_model.p_loss_of_load["bus_1", TS[1], "S1"]) < 1e-09

            @test 20. ≈ PSCOPF.get_prod_value(context.tso_schedule, "wind_1_1", TS[1], "S1") #Cprop=150 : Highlu penalized if not used
            @test 25. ≈ PSCOPF.get_prod_value(context.tso_schedule, "wind_1_2", TS[1], "S1") #Cprop=1
            @test 10. ≈ PSCOPF.get_prod_value(context.tso_schedule, "prod_1_1", TS[1], "S1") #Cprop=10
            @test PSCOPF.get_prod_value(context.tso_schedule, "prod_1_2", TS[1], "S1") < 1e-09 #Cprop=15

            # pmin = 0 for prod_1_1 & prod_2_1
            @test ismissing( PSCOPF.get_commitment_value(context.tso_schedule, "prod_1_1", TS[1], "S1") )
            @test ismissing( PSCOPF.get_commitment_value(context.tso_schedule, "prod_1_2", TS[1], "S1") )
        end

    end

    #=
    TS: [11h]
    S: [S1]
                        bus 1
                        |
    (limitable) wind_1_1|load_2
    Pmin=0, Pmax=100    | S1: 55
    Csta=0, Cprop=1     |
      S1: 45            |
                        |
    (pilotable) prod_1_1|
    Pmin=20, Pmax=100   |
    Csta=0, Cprop=10    |
                        |
    =#
    @testset "tso_respects_pmin" begin
        TS = [DateTime("2015-01-01T11:00:00")]
        ech = DateTime("2015-01-01T07:00:00")
        network = PSCOPF.Networks.Network()
        # Buses
        PSCOPF.Networks.add_new_bus!(network, "bus_1")
        # Limitables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "wind_1_1", PSCOPF.Networks.LIMITABLE,
                                                0., 100.,
                                                0., 1.,
                                                Dates.Second(0), Dates.Second(0))
        # Pilotables : have a Pmin but no start cost
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "prod_1_1", PSCOPF.Networks.PILOTABLE,
                                                20., 100.,
                                                0., 10.,
                                                Dates.Second(0), Dates.Second(0))

        # Uncertainties
        uncertainties = PSCOPF.Uncertainties()
        PSCOPF.add_uncertainty!(uncertainties, ech, "wind_1_1", DateTime("2015-01-01T11:00:00"), "S1", 45.)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", DateTime("2015-01-01T11:00:00"), "S1", 55.)
        # firmness
        firmness = PSCOPF.Firmness(
                    SortedDict("prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE),
                            ),
                    SortedDict("wind_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE),
                                "prod_1_1" => SortedDict(Dates.DateTime("2015-01-01T11:00:00") => PSCOPF.FREE), )
                    )
        # initial generators state
        generators_init_state = SortedDict(
            "prod_1_1" => PSCOPF.OFF,
        )

        context = PSCOPF.PSCOPFContext(network, TS, PSCOPF.PSCOPF_MODE_1,
                                        generators_init_state,
                                        uncertainties, nothing)
        context.market_schedule = PSCOPF.Schedule(PSCOPF.Market(), ech, SortedDict(
            "wind_1_1" => PSCOPF.GeneratorSchedule("wind_1_1",
                SortedDict(),
                SortedDict( TS[1] => PSCOPF.UncertainValue{Float64}(missing,
                                                                    SortedDict("S1"=>45.)))
            ),
            "prod_1_1" => PSCOPF.GeneratorSchedule("prod_1_1",
                SortedDict(TS[1] => PSCOPF.UncertainValue{PSCOPF.GeneratorState}(PSCOPF.ON,
                                                                                SortedDict("S1"=>PSCOPF.ON))),
                SortedDict(TS[1] => PSCOPF.UncertainValue{Float64}(missing,
                                                                    SortedDict("S1"=>10.))) # does not respect Pmin=20
                ),
            )
        )
        tso = PSCOPF.TSOOutFO()
        result = PSCOPF.run(tso, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)
        PSCOPF.update_tso_schedule!(context, ech, result, firmness, tso)

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL

        @test 20. ≈ value(result.objective_model.deltas)
        @test value(result.lol_model.p_loss_of_load["bus_1", TS[1], "S1"]) < 1e-09

        @test 35. ≈ PSCOPF.get_prod_value(context.tso_schedule, "wind_1_1", TS[1], "S1")
        @test PSCOPF.ON == PSCOPF.get_commitment_value(context.tso_schedule, "prod_1_1", TS[1], "S1")
        @test 20. ≈ PSCOPF.get_prod_value(context.tso_schedule, "prod_1_1", TS[1], "S1")

        @test value(result.objective_model.start_cost) < 1e-09 # prod_1_1 started by the market
        @test value(result.objective_model.prop_cost) ≈ (
              ((45. - 35.) * 1 + 20. * 10.)
        )
    end

end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
using LinearAlgebra

function wait_for_key(prompt)
    print(stdout, prompt)
    read(stdin, 1)
    return nothing
end

function dsyevd!(jobz, uplo, n, A, lda, w, work, lwork, iwork, liwork, info)
    ccall(
        (:dsyevd_64_, "libopenblas64_"),
        Cvoid,
        (
            Ref{UInt8},
            Ref{UInt8},
            Ref{Int},
            Ref{Float64},
            Ref{Int},
            Ref{Float64},
            Ref{Float64},
            Ref{Int},
            Ref{Int},
            Ref{Int},
            Ref{Int},
        ),
        jobz,
        uplo,
        n,
        A,
        lda,
        w,
        work,
        lwork,
        iwork,
        liwork,
        info,
    )
    if info[] != 0
        @warn("dsyevd failed: info = $(info[])")
    end
    return nothing
end

#=
 DSYRK  performs one of the symmetric rank k operations

    C := alpha*A*A**T + beta*C,

 or

    C := alpha*A**T*A + beta*C,

 where  alpha and beta  are scalars, C is an  n by n
 symmetric matrix and  A  is an  n by k  matrix in the
 first case and a  k by n  matrix in the second case.
=#
function dsyrk!(uplo, trans, n, k, alpha, A, lda, beta, C, ldc)
    ccall(
        (:dsyrk_64_, "libopenblas64_"),
        Cvoid,
        (
            Ref{UInt8},
            Ref{UInt8},
            Ref{Int},
            Ref{Int},
            Ref{Float64},
            Ref{Float64},
            Ref{Int},
            Ref{Float64},
            Ref{Float64},
            Ref{Int},
        ),
        uplo,
        trans,
        n,
        k,
        alpha,
        A,
        lda,
        beta,
        C,
        ldc,
    )
    return nothing
end

struct ProjPSD
    nmax::Int
    jobz::Char
    w::Vector{Float64}
    work::Vector{Float64}
    lwork::Int
    iwork::Vector{Int}
    liwork::Int
    info::Base.RefValue{Int}

    function ProjPSD(nmax)
        n = nmax
        A = Symmetric(zeros(1, 1))

        jobz = 'V'
        lda = n
        w = zeros(n)
        info = Ref{Int}(0)

        # Perform an optimal workspace query
        lwork = -1
        liwork = -1
        work = zeros(1)
        iwork = zeros(Int, 1)
        dsyevd!(
            jobz,
            A.uplo,
            n,
            A.data,
            lda,
            w,
            work,
            lwork,
            iwork,
            liwork,
            info,
        )
        lwork = Int(real(work[1]))
        liwork = iwork[1]
        resize!(work, lwork)
        resize!(iwork, liwork)

        new(nmax, jobz, w, work, lwork, iwork, liwork, info)
    end
end

function (obj::ProjPSD)(A::Symmetric, Aplus::Symmetric, Aminus::Symmetric)
    lda, n = size(A)

    @assert n <= obj.nmax

    dsyevd!(
        obj.jobz,
        A.uplo,
        n,
        A.data,
        lda,
        obj.w,
        obj.work,
        obj.lwork,
        obj.iwork,
        obj.liwork,
        obj.info,
    )

    λ, V = obj.w, A.data

    # V = V*Diagonal(sqrt.(abs.(λ)))
    k = 0
    @inbounds for j = 1:n
        tmp = λ[j]
        if tmp < 0.0
            k = j
            tmp *= -1.0
        end
        tmp = sqrt(tmp)
        for i = 1:n
            V[i, j] *= tmp
        end
    end

    # A = V*V'
    # A := alpha*V*V**T + beta*A, alpha = 1.0, beta = 0.0
    if k > 0
        dsyrk!(Aminus.uplo, 'N', n, k, 1.0, V, n, 0.0, Aminus.data, n)
    else
        Aminus.data .= 0.0
    end

    if k < n
        dsyrk!(
            Aplus.uplo,
            'N',
            n,
            n - k,
            1.0,
            pointer(V, n * k + 1),
            n,
            0.0,
            Aplus.data,
            n,
        )
    else
        Aplus.data .= 0.0
    end

    return nothing
end

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

max_stars_repo_name: krislock/IR-FISTA

max_issues_repo_name: krislock/IR-FISTA

max_forks_repo_name: krislock/IR-FISTA


text:
# file for misc loopset transforms

function hoist_constant_memory_accesses!(ls::LoopSet)
  hoist_stores = false
  for op ∈ operations(ls)
    if isload(op)
      length(getindicesonly(op)) == 0 && host_constant_load!(ls, op)
    elseif isstore(op) && iszero(length(getindicesonly(op)))
      hoist_stores = true
    end
  end
  hoist_stores && return hoist_constant_memory_accesses_nocheck!(ls)
  ls.preamble
end

function hoist_constant_memory_accesses_nocheck!(ls::LoopSet)
  post = Expr(:block)
  for op ∈ operations(ls)
    if isstore(op) && length(getindicesonly(op)) == 0
      hoist_constant_store!(post, ls, op)
    end
  end
  post
end
function hoist_constant_vload!(ls::LoopSet, op::Operation)
  op.instr = LOOPCONSTANT
  op.node_type = constant
  add_constant_vload!(ls, op, ArrayReferenceMetaPosition(op.ref, parents(op), loopdependencies(op), reduceddependencies(op), name(op)), elementbytes)
end

function return_empty_reductinit(op::Operation, var::Symbol)
  for (i,opp) ∈ enumerate(parents(op))
    if (name(opp) === var) && (length(reduceddependencies(opp)) == 0) && (length(loopdependencies(opp)) == 0) && (length(children(opp)) == 1)
      return opp
    end
    opcheck = return_empty_reductinit(opp, var)
    opcheck === opp || return opcheck
  end
  return op
end

function constant_symbol!(ls::LoopSet, op::Operation)
  # hack
  # relowers, but should make it work
  # TODO: DRY with `lower_licm_constants!` from `src/codegen/lower_constants.jl`
  skip_constant(instruction(op)) || return instruction(op).instr
  idcheck = identifier(op)
  symname = constantopname(op)
  for (id, sym) ∈ ls.preamble_symsym
    (idcheck ≢ nothing) && ((idcheck == id) && continue)
    pushpreamble!(ls, Expr(:(=), symname, sym))
    return symname
    # setconstantop!(ls, op,  sym)
    # setconstantop!(ls, op, Expr(:call, lv(:maybeconvert), ls.T, sym))
  end
  for (id,(intval,intsz,signed)) ∈ ls.preamble_symint
    (idcheck ≢ nothing) && ((idcheck == id) && continue)
    if intsz == 1
      pushpreamble!(ls, Expr(:(=), symname, intval % Bool))
    else
      pushpreamble!(ls, Expr(:(=), symname, sizeequivalent_symint_expr(intval, signed)))
    end
    return symname
  end
  for (id,floatval) ∈ ls.preamble_symfloat
    (idcheck ≢ nothing) && ((idcheck == id) && continue)
    pushpreamble!(ls, Expr(:(=), symname, Expr(:call, lv(:sizeequivalentfloat), ELTYPESYMBOL, floatval)))
    return symname
  end
  for (id,typ) ∈ ls.preamble_zeros
    (idcheck ≢ nothing) && ((idcheck == id) && continue)
    instruction(op) === LOOPCONSTANT || continue
    if typ == IntOrFloat
      pushpreamble!(ls, Expr(:(=), symname, Expr(:call, :zero, ELTYPESYMBOL)))
    elseif typ == HardInt
      pushpreamble!(ls, Expr(:(=), symname, Expr(:call, lv(:zerointeger), ELTYPESYMBOL)))
    else#if typ == HardFloat
      pushpreamble!(ls, Expr(:(=), symname, Expr(:call, lv(:zerofloat), ELTYPESYMBOL)))
    end
    return symname
  end
  for (id,f) ∈ ls.preamble_funcofeltypes
    (idcheck ≢ nothing) && ((idcheck == id) && continue)
    pushpreamble!(ls, Expr(:(=), symname, Expr(:call, reduction_zero(f), ELTYPESYMBOL)))
    return symname
  end
  throw("Constant operation symbol not found.")
end

function hoist_constant_store!(q::Expr, ls::LoopSet, op::Operation)
  op.instruction = DROPPEDCONSTANT
  op.node_type = constant

  opr = only(parents(op))
  while opr.instruction.instr === :identity
    opr.instruction = DROPPEDCONSTANT
    opr.node_type = constant
    opr = only(parents(opr))
  end
  push!(ls.outer_reductions, identifier(opr))

  initop = return_empty_reductinit(opr, name(opr))
  # @show last(ls.preamble.args)
  init = constant_symbol!(ls, initop)
  # @show last(ls.preamble.args)
  pushpreamble!(ls, Expr(:(=), outer_reduct_init_typename(opr), Expr(:call, lv(:typeof), init)))
  qpre = Expr(:block)
  push!(q.args, Expr(:call, lv(:unsafe_store!), Expr(:call, lv(:pointer), op.ref.ptr), outer_reduction_to_scalar_reduceq!(qpre, opr, init)))
  length(qpre.args) == 0 || pushpreamble!(ls, qpre) # creating `Expr` and pushing because `outer_reduction_to_scalar_reduceq!` uses `pushfirst!(q.args`, and we don't want it at the start of the preamble
  return nothing
end


text:
using Dilemma
using Test

@testset "Dilemma" begin
    # ...
end


text:
# Arnoldi/Lanczos iteration algorithms (with custom IOP)

#######################################
# Output type/cache
"""
    KrylovSubspace{T}(n,[maxiter=30]) -> Ks

Constructs an uninitialized Krylov subspace, which can be filled by `arnoldi!`.

The dimension of the subspace, `Ks.m`, can be dynamically altered but should
be smaller than `maxiter`, the maximum allowed arnoldi iterations.

    getV(Ks) -> V
    getH(Ks) -> H

Access methods for the (extended) orthonormal basis `V` and the (extended)
Gram-Schmidt coefficients `H`. Both methods return a view into the storage
arrays and has the correct dimensions as indicated by `Ks.m`.

    resize!(Ks, maxiter) -> Ks

Resize `Ks` to a different `maxiter`, destroying its contents.

This is an expensive operation and should be used scarcely.
"""
mutable struct KrylovSubspace{B, T, U}
    m::Int        # subspace dimension
    maxiter::Int  # maximum allowed subspace size
    augmented::Int# length of the augmented part
    beta::B       # norm(b,2)
    V::Matrix{T}  # orthonormal bases
    H::Matrix{U}  # Gram-Schmidt coefficients (real for Hermitian matrices)
    KrylovSubspace{T,U}(n::Integer, maxiter::Integer=30, augmented::Integer=false) where {T,U} =
        new{real(T), T, U}(maxiter, maxiter, augmented, zero(real(T)), Matrix{T}(undef, n + augmented, maxiter + 1),
                           fill(zero(U), maxiter + 1, maxiter + !iszero(augmented)))
    KrylovSubspace{T}(args...) where {T} = KrylovSubspace{T,T}(args...)
end
getV(Ks::KrylovSubspace) = @view(Ks.V[:, 1:Ks.m + 1])
getH(Ks::KrylovSubspace) = @view(Ks.H[1:Ks.m + 1, 1:Ks.m+!iszero(Ks.augmented)])
function Base.resize!(Ks::KrylovSubspace{B,T,U}, maxiter::Integer) where {B,T,U}
    isaugmented = !iszero(Ks.augmented)
    V = Matrix{T}(undef, size(Ks.V, 1), maxiter + 1)
    H = fill(zero(U), maxiter + 1, maxiter + isaugmented)
    if isaugmented
      copyto!(@view(V[axes(Ks.V)...]), Ks.V)
      copyto!(@view(H[axes(Ks.H)...]), Ks.H)
    end
    Ks.V = V; Ks.H = H
    Ks.m = Ks.maxiter = maxiter
    return Ks
end
function Base.show(io::IO, ::MIME"text/plain", Ks::KrylovSubspace)
    io′ = IOContext(io, :limit => true, :displaysize => 3 .* displaysize(io) .÷ 7)
    println(io, "$(Ks.m)-dimensional Krylov subspace with fields")
    println(io, "beta: $(Ks.beta)")
    println(io, "V:")
    show(io′, "text/plain", getV(Ks))
    println(io)
    println(io, "H:")
    show(io′, "text/plain", getH(Ks))
end

#######################################
# Arnoldi/Lanczos with custom IOP
## High-level interface
"""
    arnoldi(A,b[;m,tol,opnorm,iop]) -> Ks

Performs `m` anoldi iterations to obtain the Krylov subspace K_m(A,b).

The n x (m + 1) basis vectors `getV(Ks)` and the (m + 1) x m upper Hessenberg
matrix `getH(Ks)` are related by the recurrence formula

```
v_1=b,\\quad Av_j = \\sum_{i=1}^{j+1}h_{ij}v_i\\quad(j = 1,2,\\ldots,m)
```

`iop` determines the length of the incomplete orthogonalization procedure [^1].
The default value of 0 indicates full Arnoldi. For symmetric/Hermitian `A`,
`iop` will be ignored and the Lanczos algorithm will be used instead.

Refer to `KrylovSubspace` for more information regarding the output.

Happy-breakdown occurs whenver `norm(v_j) < tol * opnorm`, in this case
the dimension of `Ks` is smaller than `m`.

[^1]: Koskela, A. (2015). Approximating the matrix exponential of an
advection-diffusion operator using the incomplete orthogonalization method. In
Numerical Mathematics and Advanced Applications-ENUMATH 2013 (pp. 345-353).
Springer, Cham.
"""
function arnoldi(A, b; m=min(30, size(A, 1)), ishermitian=LinearAlgebra.ishermitian(A), kwargs...)
    TA, Tb = eltype(A), eltype(b)
    T = promote_type(TA, Tb)
    Ks = KrylovSubspace{T, ishermitian ? real(T) : T}(length(b), m)
    arnoldi!(Ks, A, b; m=m, ishermitian=ishermitian, kwargs...)
end

## Low-level interface

@inline function applyA!(y, A::AT, x, V, j, n, p) where AT
    # We cannot add `@inbounds` to `mul!`, because it is provided by the user.
    AT <: Tuple #= augmented =# || (mul!(y, A, x); return)
    A, B = A
    @inbounds begin
        # V[1:n, j + 1] = A * @view(V[1:n, j]) + B * @view(V[n+1:n+p, j])
        mul!(@view(V[1:n, j + 1]), A, @view(V[1:n, j]))
        BLAS.gemm!('N', 'N', 1.0, B, @view(V[n+1:n+p, j]), 1.0, @view(V[1:n, j + 1]))
        copyto!(@view(V[n+1:n+p-1, j + 1]), @view(V[n+2:n+p, j]))
        V[end, j + 1] = 0
    end
    return
end

function checkdims(A, b, V)
    isaugmented = b isa Tuple # isaugmented
    if isaugmented
        b′, b_aug = b
        n, p = length(b′), length(b_aug)
        _A = first(A)
    else
        n, p = size(V, 1), 0
        _A, b′, b_aug = A, b, nothing
    end
    length(b′) == size(_A,1) == size(_A,2) == size(V, 1)-p || throw(DimensionMismatch("length(b′) [$(length(b′))] == size(_A,1) [$(size(_A,1))] == size(_A,2) [$(size(_A,2))] == size(V, 1)-p [$(size(V, 1)-p)] doesn't hold"))
    return b′, b_aug, n, p
end

##############################
# Utilities
##############################
"""
    firststep!(Ks, V, H, b) -> nothing

Compute the first step of Arnoldi or Lanczos iteration.
"""
function firststep!(Ks::KrylovSubspace, V, H, b)
    @inbounds begin
        fill!(H, zero(eltype(H)))
        Ks.beta = norm(b)
        @. V[:, 1] = b / Ks.beta
    end
    return
end

"""
    firststep!(Ks, V, H, b, b_aug, t, mu, l) -> nothing

Compute the first step of Arnoldi or Lanczos iteration of augmented system.
"""
function firststep!(Ks::KrylovSubspace, V, H, b, b_aug, t, mu, l)
    @inbounds begin
        n, p = length(b), length(b_aug)
        for k=1:p-1
            i = p - k
            b_aug[k] = t^i/factorial(i) * mu
        end
        b_aug[p] = mu

        # Initialize the matrices V and H
        fill!(H, 0)

        # Normalize initial vector (this norm is nonzero)
        bl = @view b[:, l]
        Ks.beta = beta = sqrt(bl'bl + b_aug'b_aug)

        # The first Krylov basis vector
        @. V[1:n, 1]     = bl / beta
        @. V[n+1:n+p, 1] = b_aug / beta
    end
    return
end

##############################
# Arnoldi
##############################
"""
    arnoldi_step!(j, iop, n, A, V, H)

Take the `j`:th step of the Lanczos iteration.
"""
function arnoldi_step!(j::Integer, iop::Integer, A::AT,
                       V::AbstractMatrix{T}, H::AbstractMatrix{U},
                       n::Int=-1, p::Int=-1) where {AT,T,U}
    x, y = @view(V[:, j]), @view(V[:, j+1])
    applyA!(y, A, x, V, j, n, p)
    @inbounds for i = max(1, j - iop + 1):j
        α = H[i, j] = coeff(U, dot(@view(V[:, i]), y))
        axpy!(-α, @view(V[:, i]), y)
    end
    β = H[j+1, j] = norm(y)
    @. y /= β
    return β
end

"""
    arnoldi!(Ks,A,b[;tol,m,opnorm,iop,init]) -> Ks

Non-allocating version of `arnoldi`.
"""
function arnoldi!(Ks::KrylovSubspace{B, T1, U}, A::AT, b;
                  tol::Real=1e-7, m::Int=min(Ks.maxiter, size(A, 1)),
                  ishermitian::Bool=LinearAlgebra.ishermitian(A isa Tuple ? first(A) : A),
                  opnorm=nothing, iop::Int=0,
                  init::Int=0, t::Number=NaN, mu::Number=NaN, l::Int=-1) where {B, T1 <: Number, U <: Number, AT}
    ishermitian && return lanczos!(Ks, A, b; tol=tol, m=m, init=init, t=t, mu=mu, l=l)
    m > Ks.maxiter ? resize!(Ks, m) : Ks.m = m # might change if happy-breakdown occurs
    @inbounds V, H = getV(Ks), getH(Ks)
    b′, b_aug, n, p = checkdims(A, b, V)
    if iszero(init)
        isaugmented = AT <: Tuple
        isaugmented ? firststep!(Ks::KrylovSubspace, V, H, b′, b_aug, t, mu, l) : firststep!(Ks::KrylovSubspace, V, H, b)
        init = 1
    end
    iszero(iop) && (iop = m)
    for j = init:m
        beta = arnoldi_step!(j, iop, A, V, H, n, p)
        if beta < tol # happy-breakdown
            Ks.m = j
            break
        end
    end
    return Ks
end

##############################
# Lanczos
##############################

"""
    lanczos_step!(j, m, n, A, V, H)

Take the `j`:th step of the Lanczos iteration.
"""
function lanczos_step!(j::Integer, A,
                       V::AbstractMatrix{T},
                       u::AbstractVector{U},
                       v::AbstractVector{B},
                       n::Int=-1, p::Int=-1) where {B,T,U}
    x, y = @view(V[:, j]),@view(V[:, j+1])
    applyA!(y, A, x, V, j, n, p)
    α = u[j] = coeff(U, dot(x, y))
    axpy!(-α, x, y)
    j > 1 && axpy!(-v[j-1], @view(V[:, j-1]), y)
    β = v[j] = norm(y)
    @. y /= β
    return β
end

"""
    coeff(::Type,α)

Helper functions that returns the real part if that is what is
required (for Hermitian matrices), otherwise returns the value
as-is.
"""
coeff(::Type{T},α::T) where {T} = α
coeff(::Type{U},α::T) where {U<:Real,T<:Complex} = real(α)

"""
    realview(::Type, V) -> real view of `V`
"""
realview(::Type{R}, V::AbstractVector{C}) where {R,C<:Complex} =
    @view(reinterpret(R, V)[1:2:end])
realview(::Type{R}, V::AbstractVector{R}) where {R} = V

"""
    lanczos!(Ks,A,b[;tol,m,opnorm]) -> Ks

A variation of `arnoldi!` that uses the Lanczos algorithm for
Hermitian matrices.
"""
function lanczos!(Ks::KrylovSubspace{B, T1, U}, A::AT, b;
                  tol=1e-7, m=min(Ks.maxiter, size(A, 1)),
                  opnorm=nothing,
                  init::Int=0, t::Number=NaN, mu::Number=NaN, l::Int=-1) where {B, T1 <: Number, U <: Number, AT}
    m > Ks.maxiter ? resize!(Ks, m) : Ks.m = m # might change if happy-breakdown occurs
    @inbounds V, H = getV(Ks), getH(Ks)
    b′, b_aug, n, p = checkdims(A, b, V)
    if iszero(init)
        isaugmented = AT <: Tuple
        isaugmented ? firststep!(Ks::KrylovSubspace, V, H, b′, b_aug, t, mu, l) : firststep!(Ks::KrylovSubspace, V, H, b)
        init = 1
    end
    @inbounds begin
        u = @diagview(H)
        # `v` is always real, even though `u` may (in general) be complex.
        v = realview(B, @diagview(H,-1))
    end
    for j = 1:m
        if tol > lanczos_step!(j, A, V, u, v, n, p)
            # happy-breakdown
            Ks.m = j
            break
        end
    end
    @inbounds copyto!(@diagview(H,1), v[1:end-1])
    return Ks
end


text:
module PaddedViews
using Base: OneTo, tail
using OffsetArrays

export PaddedView, paddedviews, sym_paddedviews

"""
    datapadded = PaddedView(fillvalue, data, padded_axes)
    datapadded = PaddedView(fillvalue, data, padded_axes, data_axes)
    datapadded = PaddedView(fillvalue, data, sz)
    datapadded = PaddedView(fillvalue, data, sz, first_datum)
    datapadded = PaddedView{T}(args...)

Create a padded version of the array `data`, where any elements within
the span of `padded_axes` not assigned in `data` will have value
`fillvalue`.

Supply `data_axes` to specify an alterate set of axes for `data`, effectively
relocating `data` to a different set of indices.
This is shorthand for

    offsetdata = OffsetArray(data, data_axes)
    datapadded = PaddedView(fillvalue, offsetdata, padded_axes)

using the [OffsetArrays](https://github.com/JuliaArrays/OffsetArrays.jl) package.

Alternately, the padded array size `sz` can be specified, in which case `datapadded`
starts indexing at 1.
One may optionally specify the location of the `[1, 1, ...]` element of `data` with
`first_datum`.
Specifically, `datapadded[first_datum...]` corresponds to `data[1, 1, ...]`.
`first_datum` defaults to all-1s.

The view eltype `T` is optional. If not specified, then in most cases, `T` is inferred to be
`eltype(data)`. In cases when `fillvalue` can't be converted to `eltype(data)`, `T` will be
promoted the one that does. For example, when `fillvalue == nothing` and `eltype(data) == Float32`,
the inferred eltype `T` will be `Union{Nothing, Float32}`.

# Example

```jldoctest
julia> using PaddedViews

julia> a = collect(reshape(1:9, 3, 3))
3×3 $(Array{Int,2}):
 1  4  7
 2  5  8
 3  6  9

julia> PaddedView(-1, a, (4, 5))
4×5 PaddedView(-1, ::$(Array{Int,2}), (Base.OneTo(4), Base.OneTo(5))) with eltype $Int:
  1   4   7  -1  -1
  2   5   8  -1  -1
  3   6   9  -1  -1
 -1  -1  -1  -1  -1

julia> PaddedView(-1, a, (1:5,1:5), (2:4,2:4))
5×5 PaddedView(-1, OffsetArray(::$(Array{Int,2}), 2:4, 2:4), (1:5, 1:5)) with eltype $Int with indices 1:5×1:5:
 -1  -1  -1  -1  -1
 -1   1   4   7  -1
 -1   2   5   8  -1
 -1   3   6   9  -1
 -1  -1  -1  -1  -1

julia> PaddedView(-1, a, (0:4, 0:4))
5×5 PaddedView(-1, ::$(Array{Int,2}), (0:4, 0:4)) with eltype $Int with indices 0:4×0:4:
 -1  -1  -1  -1  -1
 -1   1   4   7  -1
 -1   2   5   8  -1
 -1   3   6   9  -1
 -1  -1  -1  -1  -1

julia> PaddedView(-1, a, (5,5), (2,2))
5×5 PaddedView(-1, OffsetArray(::$(Array{Int,2}), 2:4, 2:4), (Base.OneTo(5), Base.OneTo(5))) with eltype $Int:
 -1  -1  -1  -1  -1
 -1   1   4   7  -1
 -1   2   5   8  -1
 -1   3   6   9  -1
 -1  -1  -1  -1  -1
```
"""
struct PaddedView{T,N,I,A} <: AbstractArray{T,N}
    fillvalue::T
    data::A
    indices::I

    function PaddedView{T,N,I,A}(fillvalue,
                                 data,
                                 indices::NTuple{N,AbstractUnitRange}) where {T,N,I,A}
        ndims(data) == N || throw(DimensionMismatch("data and indices should have the same dimension, instead they're $(ndims(data)) and $N."))
        new{T,N,I,A}(fillvalue, data, indices)
    end
end

function PaddedView(fillvalue::FT, data::AbstractArray{T}, args...) where {FT, T}
    PaddedView{filltype(FT, T)}(fillvalue, data, args...)
end

function PaddedView(fillvalue::FT,
                    data::AbstractArray{T,N},
                    padded_inds::NTuple{N,AbstractUnitRange},
                    data_inds::NTuple{N,AbstractUnitRange}) where {FT,T,N}
    PaddedView{filltype(FT, T)}(fillvalue, data, padded_inds, data_inds)
end

function PaddedView{FT}(fillvalue,
                        data::AbstractArray{T,N},
                        indices) where {FT,T,N}
    PaddedView{FT,N,typeof(indices),typeof(data)}(convert(FT, fillvalue), data, indices)
end

function PaddedView{FT}(fillvalue,
                        data::AbstractArray{T,N},
                        sz::Tuple{Integer,Vararg{Integer}}) where {FT,T,N}
    inds = map(OneTo, sz)
    PaddedView{FT,N,typeof(inds),typeof(data)}(convert(FT, fillvalue), data, inds)
end

# No need to export this
#
# **Regardless of accuracy**, when `data` cannot represent the exact meaning of fillvalue
# object, the type of `fillvalue` should be lifted to a common type.
# Examples of this are: `Nothing`, `Missing`, `ColorTypes.Colorant`
#
# Since eltype of `data` will be lazily promoted to the filltype, it's likely to hit a
# performance issue if we abuse `filltype` by also considering storage type (e.g,
# `Float32`/`Float64`). Thus we shouldn't add lines like:
#     `filltype(::Type{FT}, ::Type{T}) where {FT<:Real, T<:Real} = promote_type(FT, T)`
# Ref: https://github.com/JuliaArrays/PaddedViews.jl/pull/25#issuecomment-610039569
filltype(::Type, ::Type{T}) where T = T
filltype(::Type{FT}, ::Type{T}) where {FT<:Union{Nothing, Missing}, T} = Union{FT, T}
filltype(::Type{FT}, ::Type{T}) where {FT, T<:Union{Nothing, Missing}} = Union{FT, T}
# ambiguity patch
filltype(::Type{FT}, ::Type{T}) where {FT<:Union{Nothing, Missing}, T<:Union{Nothing, Missing}} = Union{FT, T}

# This method eliminates an ambiguity between the two below it
function PaddedView{FT}(fillvalue,
                        data::AbstractArray{T,0},
                        ::Tuple{},
                        ::Tuple{}) where {FT,T}
    return PaddedView{FT}(fillvalue, data, ())
end

function PaddedView{FT}(fillvalue,
                        data::AbstractArray{T,N},
                        padded_inds::NTuple{N,AbstractUnitRange},
                        data_inds::NTuple{N,AbstractUnitRange}) where {FT,T,N}
    off_data = OffsetArray(data, data_inds...)
    return PaddedView{FT}(fillvalue, off_data, padded_inds)
end

function PaddedView{FT}(fillvalue,
                        data::AbstractArray{T,N},
                        sz::NTuple{N,Integer},
                        first_datum::NTuple{N,Integer}) where {FT,T,N}
    padded_inds = map(OneTo, sz)
    data_inds   = map((ax, o)->ax.+o, axes(data), first_datum .- 1)
    return PaddedView{FT}(fillvalue, data, padded_inds, data_inds)
end

Base.axes(A::PaddedView) = A.indices
@inline Base.axes(A::PaddedView, d::Integer) = d <= ndims(A) ? A.indices[d] : default_axes(A.indices)
default_axes(::NTuple{N,I}) where {N,I<:AbstractUnitRange} = convert(I, OneTo(1))
default_axes(::Tuple{}) = OneTo(1)
default_axes(::Any) = OneTo(1)

Base.size(A::PaddedView) = map(length, axes(A))

Base.parent(A::PaddedView) = A.data

Base.@propagate_inbounds function Base.getindex(A::PaddedView{T,N}, i::Vararg{Int,N}) where {T,N}
    @boundscheck checkbounds(A, i...)
    if Base.checkbounds(Bool, A.data, i...)
        return convert(T, A.data[i...])
    end
    return A.fillvalue
end

"""
    Aspad = paddedviews(fillvalue, A1, A2, ...; [dims])

Pad the arrays `A1`, `A2`, ..., to a common size or set of axes,
chosen as the span of axes enclosing all of the input arrays.

The padding is applied to one direction in dimensions `dims`. For example, values are filled to
bottom-right part of the new array in two-dimensional case. Use [`sym_paddedviews`](@ref) if
_both_ directions need to be padded.

The axes of original array `A` will be preserved in the padded result `Ap`, hence it's true
that `Ap[CartesianIndices(A)] == A`.

# Example:
```jldoctest
julia> using PaddedViews

julia> a1 = reshape([1, 2, 3], 3, 1)
3×1 $(Array{Int,2}):
 1
 2
 3

julia> a2 = [4 5 6]
1×3 $(Array{Int,2}):
 4  5  6

julia> a1p, a2p = paddedviews(-1, a1, a2);

julia> a1p
3×3 PaddedView(-1, ::$(Array{Int,2}), (Base.OneTo(3), Base.OneTo(3))) with eltype $Int:
 1  -1  -1
 2  -1  -1
 3  -1  -1

julia> a2p
3×3 PaddedView(-1, ::$(Array{Int,2}), (Base.OneTo(3), Base.OneTo(3))) with eltype $Int:
  4   5   6
 -1  -1  -1
 -1  -1  -1

julia> a1p[CartesianIndices(a1)]
3×1 $(Array{Int,2}):
 1
 2
 3
```

`dims` keyword allows padding only for specified dimensions.

```jldoctest; setup=:(using PaddedViews)
julia> a1 = reshape(collect(1:9), 3, 3)
3×3 $(Matrix{Int}):
 1  4  7
 2  5  8
 3  6  9

julia> a2 = [4 5;6 7]
2×2 $(Matrix{Int}):
 4  5
 6  7

julia> a1f, a2f = paddedviews(-1, a1, a2; dims=1);

julia> a2f
3×2 PaddedView(-1, ::$(Matrix{Int}), (Base.OneTo(3), Base.OneTo(2))) with eltype $(Int):
  4   5
  6   7
 -1  -1

julia> a1f, a2f = paddedviews(-1, a1, a2; dims=(1,2));

julia> a2f
3×3 PaddedView(-1, ::$(Matrix{Int}), (Base.OneTo(3), Base.OneTo(3))) with eltype $(Int):
  4   5  -1
  6   7  -1
 -1  -1  -1
```
"""
function paddedviews(fillvalue,  As::AbstractArray...; dims=1:ndims(first(As)))
    inds = outerinds(As...)
    map(A->PaddedView(fillvalue, A, _extended_axes(A, inds, dims)), As)
end
_extended_axes(A, inds, ::Nothing) = inds
function _extended_axes(A, inds, dims::Int)
    map((r1, r2, d)->d ? r2 : r1, axes(A), inds, ntuple(i->i==dims, ndims(A)))
end
function _extended_axes(A, inds, dims::Tuple)
    map((r1, r2, d)->d ? r2 : r1, axes(A), inds, ntuple(i->i in collect(dims), ndims(A)))
end
function _extended_axes(A, inds, dims::UnitRange{Int})
    map((r1, r2, d)->d ? r2 : r1, axes(A), inds, ntuple(i->i in collect(dims), ndims(A)))
end

# Zero, one, and two arrays are common, improve inferrability
paddedviews(fillvalue) = ()
paddedviews(fillvalue, A1::AbstractArray) = (PaddedView(fillvalue, A1, outerinds(A1)),)
function paddedviews(fillvalue, A1::AbstractArray, A2::AbstractArray)
    inds = outerinds(A1, A2)
    PaddedView(fillvalue, A1, inds), PaddedView(fillvalue, A2, inds)
end

# This is an (unexported) optimization if you're supplying the arrays from a vector.
# MosaicViews uses this.
function paddedviews_itr(fillvalue, itr)
    inds = outerinds(itr...)
    [PaddedView(fillvalue, A, inds) for A in itr]
end

@inline outerinds(A::AbstractArray, Bs...) = _outerinds(axes(A), Bs...)
@inline _outerinds(inds, A::AbstractArray, Bs...) =
    _outerinds(__outerinds(inds, axes(A)), Bs...)
_outerinds(inds) = inds
@inline function __outerinds(inds1::NTuple{N,I}, inds2) where {N,I<:AbstractUnitRange}
    map((i1, i2) -> padrange(i1, i2), inds1, inds2)
end

# This shouldn't be reached due to the `paddedviews(fillvalue)` method above,
# but it's here in case anyone extends `paddedviews` for types other than AbstractArrays.
# See https://github.com/JuliaImages/ImageCore.jl/pull/32#discussion_r111545756
outerinds() = error("must supply at least one array with concrete axes")

padrange(i1::T, i2::T) where T<: OneTo = OneTo(max(last(i1), last(i2)))
padrange(i1::T, i2::T) where T<: AbstractUnitRange = convert(T, min(first(i1),first(i2)):max(last(i1), last(i2)))
padrange(i1::AbstractRange, i2::AbstractRange) = padrange(convert(UnitRange{Int}, i1), convert(UnitRange{Int}, i2))

"""
    Aspad = sym_paddedviews(fillvalue, A1, A2, ...; [dims])

Pad the arrays `A1`, `A2`, ..., to a common size or set of axes, chosen as the span of axes
enclosing all of the input arrays.

The padding is applied to both directions in dimensions `dims`, which means original array
located at the center the padded result. Use [`paddedviews`](@ref) if only one direction
need to be padded.

The axes of original array `A` will be preserved in the padded result `Ap`, hence it's true
that `Ap[CartesianIndices(A)] == A`.

```jldoctest
julia> using PaddedViews

julia> a1 = reshape([1, 2, 3], 3, 1)
3×1 $(Array{Int,2}):
 1
 2
 3

julia> a2 = [4 5 6]
1×3 $(Array{Int,2}):
 4  5  6

julia> a1p, a2p = sym_paddedviews(-1, a1, a2);

julia> a1p
3×3 PaddedView(-1, ::$(Array{Int,2}), (1:3, 0:2)) with eltype $Int with indices 1:3×0:2:
 -1  1  -1
 -1  2  -1
 -1  3  -1

julia> a2p
3×3 PaddedView(-1, ::$(Array{Int,2}), (0:2, 1:3)) with eltype $Int with indices 0:2×1:3:
 -1  -1  -1
  4   5   6
 -1  -1  -1

julia> a1p[CartesianIndices(a1)]
3×1 $(Array{Int,2}):
 1
 2
 3
```

`dims` keyword allows padding only for specified dimensions.

```jldoctest; setup=:(using PaddedViews)
julia> a1 = reshape(collect(1:9), 3, 3)
 3×3 $(Matrix{Int}):
  1  4  7
  2  5  8
  3  6  9

julia> a2 = reshape([5, 6], 2, 1)
 2×1 $(Matrix{Int}):
  5
  6

julia> a1f, a2f = sym_paddedviews(-1, a1, a2; dims=1);

julia> a2f
 3×1 PaddedView(-1, ::$(Matrix{Int}), (1:3, 1:1)) with eltype $(Int) with indices 1:3×1:1:
   5
   6
  -1

julia> a1f, a2f = sym_paddedviews(-1, a1, a2; dims=(1,2));

julia> a2f
 3×3 PaddedView(-1, ::$(Matrix{Int}), (1:3, 0:2)) with eltype $(Int) with indices 1:3×0:2:
  -1   5  -1
  -1   6  -1
  -1  -1  -1
 ```
"""
function sym_paddedviews(fillvalue, As::AbstractArray...; dims=1:ndims(first(As)))
    inds = outerinds(As...)
    map(As) do A
        PaddedView(fillvalue, A, _sym_pad_inds(A, inds, dims))
    end
end
sym_paddedviews(fillvalue) = ()
sym_paddedviews(fillvalue, A::AbstractArray) = (PaddedView(fillvalue, A, _sym_pad_inds(A, outerinds(A))),)
function sym_paddedviews(fillvalue, A1::AbstractArray, A2::AbstractArray)
    inds = outerinds(A1, A2)
    PaddedView(fillvalue, A1, _sym_pad_inds(A1, inds)), PaddedView(fillvalue, A2, _sym_pad_inds(A2, inds))
end

function sym_paddedviews_itr(fillvalue, itr)
    inds = outerinds(itr...)
    [PaddedView(fillvalue, A, _sym_pad_inds(A, inds)) for A in itr]
end

function _sym_pad_inds(A, inds, dims=1:ndims(A))
    inds = _extended_axes(A, inds, dims)
    A_axes = axes(A)
    map(A_axes, inds) do ax, i
        pad_sz = length(i) - length(ax)
        offset = pad_sz ÷ 2
        first(ax)-offset:last(ax)+pad_sz-offset
    end
end

function Base.showarg(io::IO, A::PaddedView, toplevel)
    print(io, "PaddedView(", A.fillvalue, ", ")
    Base.showarg(io, parent(A), false)
    print(io, ", (", join(A.indices, ", "))
    print(io, ndims(A) == 1 ? ",))" : "))")
    toplevel && print(io, " with eltype ", eltype(A))
end

end # module


text:
#=
Compute all mixed Nash equilibria of a 2-player (non-degenerate) normal
form game by support enumeration.

Julia version of QuantEcon.py/support_enumeration.py

References
----------
B. von Stengel, "Equilibrium Computation for Two-Player Games in
Strategic and Extensive Form," Chapter 3, N. Nisan, T. Roughgarden, E.
Tardos, and V. Vazirani eds., Algorithmic Game Theory, 2007.
=#

using LinearAlgebra: LAPACKException, SingularException
using QuantEcon: next_k_array!

"""
    support_enumeration(g)

Compute mixed-action Nash equilibria with equal support size
for a 2-player normal form game by support enumeration. For a
non-degenerate game input, these are all the Nash equilibria.

The algorithm checks all the equal-size support pairs; if the
players have the same number n of actions, there are 2n choose n
minus 1 such pairs. This should thus be used only for small games.

# Arguments

- `g::NormalFormGame{2,T}`: 2-player NormalFormGame instance.

# Returns

- `::Vector{NTuple{2,Vector{S}}}`: Mixed-action Nash equilibria that are found,
  where `S` is Float if `T` is Int or Float, and Rational if `T` is Rational.
"""
function support_enumeration(g::NormalFormGame{2,T}) where T
    S = typeof(zero(T)/one(T))
    c = Channel{Tuple{Vector{S},Vector{S}}}(0)
    task = support_enumeration_task(c, g)
    bind(c, task)
    schedule(task)
    NEs = collect(c)

    return NEs

end

"""
    support_enumeration_task(c, g)

Task version of `support_enumeration`.

# Arguments

- `c::Channel`: Channel to be binded with the support enumeration task.
- `g::NormalFormGame{2}`: 2-player NormalFormGame instance.

# Returns

- `::Task`: Runnable task for generating Nash equilibria.
"""
function support_enumeration_task(c::Channel,
                                  g::NormalFormGame{2})

    task = Task(
        () -> _support_enumeration_producer(c,
                                            (g.players[1].payoff_array,
                                             g.players[2].payoff_array))
    )

    return task
end

"""
    _support_enumeration_producer(c, payoff_matrices)

Main body of `support_enumeration_task`.

# Arguments

- `c::Channel`: Channel to be binded with the support enumeration task.
- `payoff_matrices::NTuple{2,Matrix{T}}`: Payoff matrices of player 1 and
  player 2, where `T<:Real`.

# Puts

- `NTuple{2,Vector{S}}`: Tuple of Nash equilibrium mixed actions, where `S` is
  Float if `T` is Int or Float, and Rational if `T` is Rational.
"""
function _support_enumeration_producer(c::Channel,
                                       payoff_matrices
                                       ::NTuple{2,Matrix{T}}) where T<:Real

    nums_actions = size(payoff_matrices[1], 1), size(payoff_matrices[2], 1)
    n_min = min(nums_actions...)
    S = typeof(zero(T)/one(T))

    for k = 1:n_min
        supps = (collect(1:k), Vector{Int}(undef, k))
        actions = (Vector{S}(undef, k), Vector{S}(undef, k))
        A = Matrix{S}(undef, k+1, k+1)
        b = Vector{S}(undef, k+1)
        while supps[1][end] <= nums_actions[1]
            supps[2][:] = collect(1:k)
            while supps[2][end] <= nums_actions[2]
                if _indiff_mixed_action!(A, b, actions[2],
                                         payoff_matrices[1],
                                         supps[1], supps[2])
                    if _indiff_mixed_action!(A, b, actions[1],
                                             payoff_matrices[2],
                                             supps[2], supps[1])
                        out = (zeros(S, nums_actions[1]),
                               zeros(S, nums_actions[2]))
                        for (p, (supp, action)) in enumerate(zip(supps,
                                                                 actions))
                            out[p][supp] = action
                        end
                        put!(c, out)
                    end
                end
                next_k_array!(supps[2])
            end
            next_k_array!(supps[1])
        end
    end

end

function _solve!(A::Matrix{T}, b::Vector{T}) where T <: Union{Float64,Float32}
    r = 0
    try
        LAPACK.gesv!(A, b)
    catch LAPACKException
        r = 1
    end
    return r
end

@inline function _solve!(A::Matrix{Rational{T}},
                         b::Vector{Rational{T}}) where T <: Integer
    r = 0
    try
        b[:] = ldiv!(lu!(A), b)
    catch SingularException
        r = 1
    end
    return r
end

"""
    _indiff_mixed_action!(A, b, out, payoff_matrix, own_supp, opp_supp)

Given a player's payoff matrix `payoff_matrix`, an array `own_supp`
of this player's actions, and an array `opp_supp` of the opponent's
actions, each of length k, compute the opponent's mixed action whose
support equals `opp_supp` and for which the player is indifferent
among the actions in `own_supp`, if any such exists. Return `true`
if such a mixed action exists and actions in `own_supp` are indeed
best responses to it, in which case the outcome is stored in `out`;
`false` otherwise. Arrays `A` and `b` are used in intermediate
steps.

# Arguments

- `A::Matrix{T}`: Matrix of shape (k+1, k+1) used in intermediate steps, where
  `T<:Real`.
- `b::Vector{T}`: Vector of length k+1 used in intermediate steps, where
  `T<:Real`.
- `out::Vector{T}`: Vector of length k to store the nonzero values of the
  desired mixed action, where `T<:Real`.
- `payoff_matrix::Matrix`: The player's payoff matrix, of shape (m, n).
- `own_supp::Vector{Int}`: Vector containing the player's action indices, of
  length k.
- `opp_supp::Vector{Int}`: Vector containing the opponent's action indices, of
  length k.

# Returns

- `::Bool`: `true` if a desired mixed action exists and `false` otherwise.
"""
function _indiff_mixed_action!(A::Matrix{T}, b::Vector{T},
                               out::Vector{T},
                               payoff_matrix::Matrix,
                               own_supp::Vector{Int},
                               opp_supp::Vector{Int}) where T<:Real

    m = size(payoff_matrix, 1)
    k = length(own_supp)

    for j in 1:k, i in 1:k
        A[i, j] = payoff_matrix[own_supp[i], opp_supp[j]]
    end
    A[1:end-1, end] .= -one(T)
    A[end, 1:end-1] .= one(T)
    A[end, end] = zero(T)
    b[1:end-1] .= zero(T)
    b[end] = one(T)

    r = _solve!(A, b)
    r == 0 || return false  # A: singular

    for i in 1:k
        b[i] <= zero(T) && return false
    end

    out[:] = b[1:end-1]
    val = b[end]

    if k == m
        return true
    end

    own_supp_flags = falses(m)
    own_supp_flags[own_supp] .= true

    for i = 1:m
        if !own_supp_flags[i]
            payoff = zero(T)
            for j = 1:k
                payoff += payoff_matrix[i, opp_supp[j]] * out[j]
            end
            if payoff > val
                return false
            end
        end
    end

    return true
end