text:
# Autogenerated wrapper script for libinput_jll for x86_64-linux-gnu
export libinput

using eudev_jll
using mtdev_jll
using libevdev_jll
JLLWrappers.@generate_wrapper_header("libinput")
JLLWrappers.@declare_library_product(libinput, "libinput.so.10")
function __init__()
    JLLWrappers.@generate_init_header(eudev_jll, mtdev_jll, libevdev_jll)
    JLLWrappers.@init_library_product(
        libinput,
        "lib/libinput.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


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

using Markdown
using InteractiveUtils

# ╔═╡ 628d5d5d-2753-47e3-a02b-21b4a89a159e
import Pkg

# ╔═╡ 9c58eb1a-313f-4f89-9958-7f33c7a072d3
Pkg.activate(".")

# ╔═╡ 301a110f-2fed-45d6-b7ae-b048a3b7209b
using Revise

# ╔═╡ 55bdcbfc-3e78-4db7-ab05-7a9abe9fd253
using DeltaRCWA

# ╔═╡ 87c60eb1-495f-424d-b31a-5ab47d0e1879
using LinearMaps

# ╔═╡ 65324c70-07b4-46b8-9d6f-3b7fc58d3fbf
using Plots

# ╔═╡ 0a627414-a5af-4fc3-852f-c98105e4d860
md"
## Using `DeltaRCWA`'s interface
These are the stages to using the solver, which needs this data
```julia
struct DeltaRCWAProblem{N, T<:SheetStack}
    stack::T
    pw::PlaneWaves{N}
    I₁::Array{ComplexF64, N}
    I₂::Array{ComplexF64, N}
end
```
as explained and demonstrated below:
### Defining the scattering modes
A struct called `PlaneWaves{N}` is used to specify the discretization and
periodicities of space in the unit cell where `DeltaRCWA` solves for the fields.
`N` refers to the number of dimensions along which the unit cell is periodic.
Along these periodic dimensions, the solutions can be expanded in the Fourier basis
due to Bloch's theorem (much like the [Kronig-Penny model](https://en.wikipedia.org/wiki/Particle_in_a_one-dimensional_lattice)).
The full structure is:
```julia
struct PlaneWaves{N}
    ω::Float64
    dims::NTuple{N, Tuple{Int64, Float64}}
    x⃗::NTuple{N, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}
    k⃗::NTuple{N, Frequencies{Float64}}
end
```
This data can be supplied with the following information (using convenience methods):
"

# ╔═╡ f0ca02cb-c7b0-42d8-b598-98499aa09f12
md"
#### Wave frequency
Choose an `ω::Float64` for the temporal frequency of the wave
"

# ╔═╡ 4026eb54-9b01-48ba-8bb2-3b2901a65558
ω₀ = 10.0

# ╔═╡ 6104f60d-d961-4963-b7b9-b8a6f5c891be
k₀ = ω₀ # ω/c (c=1)

# ╔═╡ cb2f3944-438f-4a42-81d7-00f9ddb19a67
ω = 1.1ω₀

# ╔═╡ bcaf5797-af13-443c-b670-84ea72d46f8d
md"
#### Unit cell dimensions and discretization
For each periodic dimension, give a number of grid points per unit cell and a size of
the unit cell: `dims::NTuple{N, Tuple{Int64,Float64}}`. This defines the dimensionality of the problem
"

# ╔═╡ 1e27a04a-0024-4dd4-8ee2-2cae0ef23d53
Nmodes = 50

# ╔═╡ 017beeeb-327c-4c6b-b60a-e0d5041149d0
λ = 2*pi/k₀

# ╔═╡ df75e240-a9f8-4543-be3f-bf034b152166
md"
#### Wave medium
Planewaves are the eigenmodes of uniform media. Define the medium properties with
either constructor
```julia
struct UniformMedium{ϵ, μ}
    function UniformMedium{ϵ, μ}() where {ϵ, μ}
        @assert ϵ isa Number
        @assert μ isa Number
        new{ϵ, μ}()
    end
end
const Vacuum = UniformMedium{1, 1}()
```

The default value is always vacuum
"

# ╔═╡ 22fd2b7b-4418-4cd6-9ff0-eecdf1775d9b
md"
#### Polarization
This is a trait of `PlaneWaves{N}`, and is relevant information only for 2D
photonic crystals with diagonal conductivity matrices, when the TE and TM
polarizations decouple. By default the solver computes the `TM()` modes when
`N=1`. `Coupled()` exists for `N=2`. Override the default with the line below
"

# ╔═╡ 324aec31-05f2-4fb2-8825-79b2402ec7ee
# DeltaRCWA.PolarizationStyle(::Type{PlaneWaves{1}}) = TE()

# ╔═╡ ce61314c-5ba5-4422-88ce-21f2b20e8110
md"
#### Incident mode amplitudes
Defining a distribution over the incident planewave modes allows for a solution of the
scattering problem. Two distributions must be specified: one at each port of the device.
"

# ╔═╡ d0d638f3-dc93-48fa-b95d-9fc8b20e22f7
### choose which mode to scatter
modeN = 1

# ╔═╡ 4a5de2d0-28e5-4963-b9b3-95fd13f2dd95
I₁ = [n == modeN ? 1 : 0 for n in 1:Nmodes]

# ╔═╡ 1bc23e37-fb20-4d00-95b4-aa7334e18578
I₂ = zeros(Nmodes)

# ╔═╡ 3e0209fe-f316-436d-b017-422844535579
md"
### Defining the scattering structure

#### Creating a `Sheet`
To create metasurface parameters for the unit cell of a sheet-like photonic crystal, create a struct that is a subtype
of `Sheet` and store all the geometric/metasurface parameters you need in
your struct to define the electric and magnetic impedances in the unit cell for it.
"

# ╔═╡ a4e6f529-b914-4445-a6ee-7a316e259b9e
struct ComplexExpSheet{T} <: Sheet
    θ::T # incidence angle
    θᵗ::T # transmission angle
	k::T # wavenumber
    d::T # designer parameter
end

# ╔═╡ 77281f06-fcc9-4f75-a415-217cf9e581c1
θ = -π/2.0

# ╔═╡ 1192e390-1f87-486b-b6b4-7447f07f407b
θᵗ = -π/3.0

# ╔═╡ 83eb9d5f-c836-4dd1-be15-1320d94cd743
d = cos(θᵗ)-cos(θ)

# ╔═╡ cef19064-86de-4e25-a5cc-4da77f4b6c2c
L = λ/abs(d) # unit-cell width/periodicity of ComplexExpSheet

# ╔═╡ 4125d1a4-2a57-431a-b7ea-ab8f44994143
dims = ((Nmodes, L), )

# ╔═╡ 6bfcdf25-1f4e-4724-a4bc-14a0222e2c2d
sheet = ComplexExpSheet{Float64}(θ, θᵗ, k₀, d)

# ╔═╡ a13facc7-06b6-4559-bf2f-a43afb3ffcfe
md"
#### Creating a `SheetStack`
In order to put 2 or more sheets in a sequence, separated by a uniform medium create an instance of a `SheetStack`.
Create a Tuple of the sheets you want to scatter off of and create a second Tuple
with the size of the Vacuum gap that separates each of the sheets and pass these to
the `SheetStack` constructor. Note that there is one gap fewer than the number of sheets.
```julia
struct SheetStack{L, T<:Tuple{Sheet, Vararg{Sheet, L}}, M<:Tuple{UniformMedium, UniformMedium, Vararg{UniformMedium, L}}}
    depths::Tuple{Vararg{Float64, L}}
    sheets::T
    media::M
end
```
"

# ╔═╡ 8297d683-f888-419c-812b-294397b16126
nsheets = 2

# ╔═╡ 989a6a21-71a8-4991-8a55-eb7f1d3d5f70
gap(x) = 2L

# ╔═╡ 73d6e43f-214b-47ea-b2fd-df60dc412ae1
stack = SheetStack(
	Tuple(gap(i) for i in 1:(nsheets-1)),
	Tuple(sheet for i in 1:nsheets),
)

# ╔═╡ af78faae-2d92-4035-81d5-9856d662e876
md"
#### Methods that a `RCWASheet` implements
Any `RCWASheet` may implement the following 16 methods
$(vec([Symbol(R, s, i, j) for j in (:ˣ, :ʸ), i in (:ˣ, :ʸ), s in (:ₑ, :ₘ), R in (:Z, :Y)]))
one for each component of the electric and magnetic admittance/impedance matrices.
These methods return the value of that admittance/impedance component at points in the
unit cell, thereby realizing the metasurface geometry from the parameters.
By default/fallback, these methods all return zero
(thus are either perfect insulators, `Y=0`, or are perfect conductors, `Z=0`)
unless you
[extend `DeltaRCWA`'s method](https://docs.julialang.org/en/v1/manual/modules/#using-and-import-with-specific-identifiers,-and-adding-methods)
to dispatch on your type.
In a 1D problem, the components used depend on the polarization of light.
Currently the implementation in 1D assumes the off-diagonal components are zero and
neglects them.

See the examples/documentation/help for the method signature
"

# ╔═╡ 030ef40d-043c-4e23-9fc4-7a5321cf67f9
function DeltaRCWA.Zₑˣˣ(sheet::ComplexExpSheet, x⃗::Tuple)
	-0.5sin(sheet.θ)*(1-exp(1im*sheet.k*sheet.d*x⃗[1]))
end

# ╔═╡ 6046fe03-6c7f-4497-a1e1-1368ab043489
function DeltaRCWA.Yₘʸʸ(sheet::ComplexExpSheet, x⃗::Tuple)
	-2sin(sheet.θ)*(1+exp(1im*sheet.k*sheet.d*x⃗[1]))
end

# ╔═╡ f2332c95-e784-41c6-a8b9-96df60bb58a6
md"
#### Response style
To choose between specifying material response parameters `Z` or `Y` (where one is the inverse of the other), you specify traits called `ElectricResponseStyle` and `MagneticResponseStyle` for your types which may return `Impedance()` or `Admittance()`.
The default values are
```julia
ElectricResponseStyle(::Type{<:RCWASheet}) = Impedance()
MagneticResponseStyle(::Type{<:RCWASheet}) = Admittance()
```
This lets you choose how to specify the material parameters (useful if one has a value of 0 (i.e. the other diverges)) exactly to prevent singular matrices.

Uncomment and run the code below to reproduce the same example with a different `ElectricResponseStyle`.
"

# ╔═╡ bf0fad44-d478-48c6-9c05-33cd6bf61a7b
# begin
# 	DeltaRCWA.ElectricResponseStyle(::Type{ComplexExpSheet}) = Admittance()
# 	function DeltaRCWA.Yₑˣˣ(sheet::ComplexExpSheet, x⃗)
# 		# this is 1 ./ Zₑˣˣ except with a small perturbation to be nonsingular
# 		2 / (-sin(sheet.θ)*(1e-10 + 1-exp(1im*sheet.k*sheet.d*x⃗[1])))
# 	end
# end

# ╔═╡ f89cf740-818a-45fc-8ec2-5c747fc65fa9
md"
### Problem and solution
Create `DeltaRCWAProblem`s, `solve` them, and receive a `DeltaRCWASolution`.
"

# ╔═╡ aad04ca7-b150-49c2-8a8a-4dc792959649
prob = DeltaRCWAProblem(sheet, dims, ω, I₁, I₂)

# ╔═╡ 786a6947-7082-4902-a625-8be4bd3e30e7
begin
	### Display the magnetic conductivity / electric resistivity  along sheet
	Z = [-0.5sin(sheet.θ)*(1-exp(1im*sheet.k*sheet.d*x)) for x in prob.pw.x⃗[1]]
	Y = [-2sin(sheet.θ)*(1+exp(1im*sheet.k*sheet.d*x)) for x in prob.pw.x⃗[1]]
	plot(prob.pw.x⃗[1],  real.(Z), label="Re(Zₑˣˣ)")
	plot!(prob.pw.x⃗[1], imag.(Z), label="Im(Zₑˣˣ)", ls=:dash)
	plot!(prob.pw.x⃗[1], real.(Y), label="Re(Yₘʸʸ)")
	plot!(prob.pw.x⃗[1], imag.(Y), label="Im(Yₘʸʸ)", ls=:dash)
end

# ╔═╡ 386f61fa-ab15-40bc-b9be-593622ad42da
sol = solve(prob)

# ╔═╡ fe8b8e5a-7de9-4655-ab19-076d58ce0143
stackprob = DeltaRCWAProblem(stack, dims, ω, I₁, I₂)

# ╔═╡ 20e74b68-23a4-47f1-acea-0b244efd25ae
stacksol = solve(stackprob, T=Matrix)

# ╔═╡ 2bc3991d-9626-4464-a532-5ab7705e454b
md"
### Analysis
The `DeltaRCWASolution` objects can have recipes to analyze and visualize them.
#### Plotting
`plot(::DeltaRCWASolution{1})` should just work!
"

# ╔═╡ 77400f50-4e25-4fe6-8a0a-16f6cf6cb150
plot(sol, combine=true, aspect_ratio=1, color=:RdBu)

# ╔═╡ d8d2123c-92ce-422e-bb68-bb09e689d44c
plot(stacksol, combine=true, color=:RdBu, aspect_ratio=1)

# ╔═╡ 10442fb5-d595-4ef6-a1e9-4d588609dbaf
md"""
## Example of 3D interface
"""

# ╔═╡ 4e412a25-1176-4db9-955b-130fd87e734b
Mmodes = 1 # additional points along y dimension

# ╔═╡ feb92d68-8034-44bf-aba1-bec640315d8f
NMmodes = Nmodes * Mmodes

# ╔═╡ 38dca2e1-2788-4bd9-ad9f-1cfb78bab06c
I₁ˣ = [n == m == modeN ? 0 : 0 for n in 1:Nmodes, m in 1:Mmodes]

# ╔═╡ b6fed619-11f8-40a2-80dd-123adaf9c77c
I₂ˣ = [n == m == modeN ? 0 : 0 for n in 1:Nmodes, m in 1:Mmodes]

# ╔═╡ 3aafea78-513f-4a5c-8a1d-5b689a358a4e
I₁ʸ = [n == m == modeN ? 1 : 0 for n in 1:Nmodes, m in 1:Mmodes]

# ╔═╡ 213628fd-77a1-4753-bc30-f78c36ddd801
I₂ʸ = [n == m == modeN ? 0 : 0 for n in 1:Nmodes, m in 1:Mmodes]

# ╔═╡ bf01b9d0-4437-422f-a769-05955a41e2b1
dims3d = ((Nmodes, L), (Mmodes, L/L))

# ╔═╡ c68fc7c6-5eae-4984-82bf-4368ad910610
prob3d = DeltaRCWAProblem(sheet, dims3d, ω, hcat(I₁ˣ, I₁ʸ), hcat(I₂ˣ, I₂ʸ))

# ╔═╡ b966638a-fca0-4431-ac6a-e094c03b138c
sol3d = solve(prob3d)

# ╔═╡ 7eb3f960-c89a-40c7-9c20-f157a2c394a1
z⃗ = range(-L, L, length=2Nmodes)

# ╔═╡ d487aace-25f6-4706-985a-c52bb5799612
z⃗₁ = z⃗[z⃗ .< 0]

# ╔═╡ 458e253d-f1f0-4e26-bd4a-831a44cd6722
z⃗₂ = z⃗[z⃗ .> 0]

# ╔═╡ d5d5ac2a-c94d-4695-8c19-63b9147c29fe
begin
	# quick plotting
	kz = DeltaRCWA._get_kz(sol.pw, sol3d.stack.media[1])
	CO₁ = rotr90(DeltaRCWA.bfft(exp.(-kz * transpose(im * z⃗₁)) .* sol3d.O₁[NMmodes.+(1:Nmodes)], 1))
	CI₂ = rotr90(DeltaRCWA.bfft(exp.(-kz * transpose(im * z⃗₂)) .* sol3d.I₂[NMmodes.+(1:Nmodes)], 1))
	CI₁ = rotr90(DeltaRCWA.bfft(exp.( kz * transpose(im * z⃗₁)) .* sol3d.I₁[NMmodes.+(1:Nmodes)], 1))
	CO₂ = rotr90(DeltaRCWA.bfft(exp.( kz * transpose(im * z⃗₂)) .* sol3d.O₂[NMmodes.+(1:Nmodes)], 1))
	heatmap(sol.pw.x⃗, z⃗, real.(cat(CI₁ + CO₁, CI₂ + CO₂; dims=1)), xguide="x", yguide="z", aspect_ratio=:equal,color=:RdBu)
end

# ╔═╡ 301ad780-83a1-46f3-a027-ef8a4ca198f6
prob3dstack = DeltaRCWAProblem(stack, dims3d, ω, hcat(I₁ˣ, I₁ʸ), hcat(I₂ˣ, I₂ʸ))

# ╔═╡ f155c239-2001-4fd1-84b2-558827c7cdd2
sol3dstack = solve(prob3dstack)

# ╔═╡ d7372f04-0ede-42ee-b9f5-0e96a2004d68
begin
	# quick plotting
	DO₁ = rotr90(DeltaRCWA.bfft(exp.(-kz * transpose(im * z⃗₁)) .* sol3dstack.O₁[NMmodes.+(1:Nmodes)], 1))
	DI₂ = rotr90(DeltaRCWA.bfft(exp.(-kz * transpose(im * z⃗₂)) .* sol3dstack.I₂[NMmodes.+(1:Nmodes)], 1))
	DI₁ = rotr90(DeltaRCWA.bfft(exp.( kz * transpose(im * z⃗₁)) .* sol3dstack.I₁[NMmodes.+(1:Nmodes)], 1))
	DO₂ = rotr90(DeltaRCWA.bfft(exp.( kz * transpose(im * z⃗₂)) .* sol3dstack.O₂[NMmodes.+(1:Nmodes)], 1))
	heatmap(sol.pw.x⃗, z⃗, real.(cat(DI₁ + DO₁, DI₂ + DO₂; dims=1)), xguide="x", yguide="z", aspect_ratio=:equal,color=:RdBu)
	end

# ╔═╡ Cell order:
# ╠═628d5d5d-2753-47e3-a02b-21b4a89a159e
# ╠═9c58eb1a-313f-4f89-9958-7f33c7a072d3
# ╠═301a110f-2fed-45d6-b7ae-b048a3b7209b
# ╠═55bdcbfc-3e78-4db7-ab05-7a9abe9fd253
# ╠═87c60eb1-495f-424d-b31a-5ab47d0e1879
# ╠═65324c70-07b4-46b8-9d6f-3b7fc58d3fbf
# ╟─0a627414-a5af-4fc3-852f-c98105e4d860
# ╟─f0ca02cb-c7b0-42d8-b598-98499aa09f12
# ╠═4026eb54-9b01-48ba-8bb2-3b2901a65558
# ╠═6104f60d-d961-4963-b7b9-b8a6f5c891be
# ╠═cb2f3944-438f-4a42-81d7-00f9ddb19a67
# ╟─bcaf5797-af13-443c-b670-84ea72d46f8d
# ╠═1e27a04a-0024-4dd4-8ee2-2cae0ef23d53
# ╠═017beeeb-327c-4c6b-b60a-e0d5041149d0
# ╠═cef19064-86de-4e25-a5cc-4da77f4b6c2c
# ╠═4125d1a4-2a57-431a-b7ea-ab8f44994143
# ╟─df75e240-a9f8-4543-be3f-bf034b152166
# ╟─22fd2b7b-4418-4cd6-9ff0-eecdf1775d9b
# ╠═324aec31-05f2-4fb2-8825-79b2402ec7ee
# ╟─ce61314c-5ba5-4422-88ce-21f2b20e8110
# ╠═d0d638f3-dc93-48fa-b95d-9fc8b20e22f7
# ╠═4a5de2d0-28e5-4963-b9b3-95fd13f2dd95
# ╠═1bc23e37-fb20-4d00-95b4-aa7334e18578
# ╟─3e0209fe-f316-436d-b017-422844535579
# ╠═a4e6f529-b914-4445-a6ee-7a316e259b9e
# ╠═77281f06-fcc9-4f75-a415-217cf9e581c1
# ╠═1192e390-1f87-486b-b6b4-7447f07f407b
# ╠═83eb9d5f-c836-4dd1-be15-1320d94cd743
# ╠═6bfcdf25-1f4e-4724-a4bc-14a0222e2c2d
# ╟─a13facc7-06b6-4559-bf2f-a43afb3ffcfe
# ╠═8297d683-f888-419c-812b-294397b16126
# ╠═989a6a21-71a8-4991-8a55-eb7f1d3d5f70
# ╠═73d6e43f-214b-47ea-b2fd-df60dc412ae1
# ╟─af78faae-2d92-4035-81d5-9856d662e876
# ╠═030ef40d-043c-4e23-9fc4-7a5321cf67f9
# ╠═6046fe03-6c7f-4497-a1e1-1368ab043489
# ╠═786a6947-7082-4902-a625-8be4bd3e30e7
# ╟─f2332c95-e784-41c6-a8b9-96df60bb58a6
# ╠═bf0fad44-d478-48c6-9c05-33cd6bf61a7b
# ╟─f89cf740-818a-45fc-8ec2-5c747fc65fa9
# ╠═aad04ca7-b150-49c2-8a8a-4dc792959649
# ╠═386f61fa-ab15-40bc-b9be-593622ad42da
# ╠═fe8b8e5a-7de9-4655-ab19-076d58ce0143
# ╠═20e74b68-23a4-47f1-acea-0b244efd25ae
# ╟─2bc3991d-9626-4464-a532-5ab7705e454b
# ╠═77400f50-4e25-4fe6-8a0a-16f6cf6cb150
# ╠═d8d2123c-92ce-422e-bb68-bb09e689d44c
# ╟─10442fb5-d595-4ef6-a1e9-4d588609dbaf
# ╠═4e412a25-1176-4db9-955b-130fd87e734b
# ╠═feb92d68-8034-44bf-aba1-bec640315d8f
# ╠═38dca2e1-2788-4bd9-ad9f-1cfb78bab06c
# ╠═b6fed619-11f8-40a2-80dd-123adaf9c77c
# ╠═3aafea78-513f-4a5c-8a1d-5b689a358a4e
# ╠═213628fd-77a1-4753-bc30-f78c36ddd801
# ╠═bf01b9d0-4437-422f-a769-05955a41e2b1
# ╠═c68fc7c6-5eae-4984-82bf-4368ad910610
# ╠═b966638a-fca0-4431-ac6a-e094c03b138c
# ╠═7eb3f960-c89a-40c7-9c20-f157a2c394a1
# ╠═d487aace-25f6-4706-985a-c52bb5799612
# ╠═458e253d-f1f0-4e26-bd4a-831a44cd6722
# ╠═d5d5ac2a-c94d-4695-8c19-63b9147c29fe
# ╠═301ad780-83a1-46f3-a027-ef8a4ca198f6
# ╠═f155c239-2001-4fd1-84b2-558827c7cdd2
# ╠═d7372f04-0ede-42ee-b9f5-0e96a2004d68

ext: jl


text:
#==========================================================
Code for comparing solvers

For comparing possibly flaky code, we need an ability to stop the code
if it runs for too long.
Code that does that is in the compare directory.

It uses threads for Julia code,
and has special code for the matlab solvers

===========================================================#

"""
    SolverTest(solver, name)

Encloses a solver with its name, so that we can compare it in tests
"""
struct SolverTest
    solver::Function
    name::String
end

"""
    initDictCol!(dic, name, typ)

For a dictionary in which each key indexes an array.
If dic does not contain an entry of `name`, create with set to `Array(typ,0)`.
"""
function initDictCol!(dic, name, typ)
    if ~haskey(dic,name)
        dic[name] = typ[]
    end
end

"""
`ret` is the answer returned by a speed test.
This pushed it into the dictionary on which we are storing the tests.
"""
function pushSpeedResult!(dic, name, ret)
    push!(dic["$(name)_solve"],ret[1])
    push!(dic["$(name)_build"],ret[2])
    push!(dic["$(name)_tot"],ret[1]+ret[2])
    push!(dic["$(name)_its"],ret[3])
    push!(dic["$(name)_err"],ret[4])
end

"""
    function speedTestLapSolvers{Tv,Ti}(solvers, dic, a::SparseMatrixCSC{Tv,Ti}, b::Array{Tv,1}; tol::Real=1e-2, maxits=Inf, maxtime=Inf, verbose=false)

Runs many Laplacians solvers.  Puts the build and solve time results into a dictionary dic.  It would be easiest to look at it via DataFrame(dic).  Returns the answer from the last solver.  `solvers` should be an array of `SolverTest`.
"""
function speedTestLapSolvers(solvers, dic, a::SparseMatrixCSC{Tv,Ti}, b::Array{Tv,1}; tol::Real=1e-2, maxits=10000, maxtime=1000, verbose=false, testName="") where {Tv,Ti}

    b = b .- mean(b)

    la = lap(a)

    it = Int[1]

    # init the dict, if it is the first time
    initDictCol!(dic, "nv", Int)
    initDictCol!(dic, "ne", Int)
    initDictCol!(dic, "hash_a", UInt64)
    initDictCol!(dic, "testName", String)

    solvecol(name) = "$(name)_solve"
    buildcol(name) = "$(name)_build"
    totcol(name) = "$(name)_tot"
    itscol(name) = "$(name)_its"
    errcol(name) = "$(name)_err"

    dic["names"] = [t.name for t in solvers]

    for solver in solvers
        name = solver.name
        initDictCol!(dic, solvecol(name), Float64)
        initDictCol!(dic, buildcol(name), Float64)
        initDictCol!(dic, totcol(name), Float64)
        initDictCol!(dic, itscol(name), Float64)
        initDictCol!(dic, errcol(name), Float64)
    end

    nv = size(a,1)
    ne = nnz(a)
    hash_a = hash(a)

    push!(dic["nv"],nv)
    push!(dic["ne"],ne)
    push!(dic["hash_a"],hash_a)
    push!(dic["testName"],testName)

    x = []

    for i in 1:length(solvers)
        solverTest = solvers[i]

        if verbose
            println()
            println(solverTest.name)
        end

        ret = testSolver(solverTest.solver, a, b, tol, maxits, verbose)

        if i == 1
            x = ret[5]
        end

        pushSpeedResult!(dic, solverTest.name, ret)
    end

    return x

end

function testSolver(solver, a, b, tol, maxits, verbose)

  try

    GC.gc()
    t0 = time()
    f = solver(a, tol=tol, maxits=maxits, verbose=verbose)
    build_time = time() - t0

    it = [0]
    GC.gc()

    t0 = time()
    x = f(b, pcgIts = it, tol=tol, maxits=maxits, verbose=verbose)
    solve_time = time() - t0

    err = norm(lap(a) * x .- b) / norm(b)

    ret = (solve_time, build_time, it[1], err, x)
    if verbose
      println((solve_time, build_time, it[1], err))
    end
    return ret
  catch
    println("Solver Error.")
    return (Inf, Inf, Inf, Inf)
  end

end


text:
using Test, BenchmarkTools, StaticArrays, LinearAlgebra, KiteUtils

using KiteModels, KitePodModels

const SEGMENTS = se().segments
if ! @isdefined kcu
    const kcu = KCU(se())
    const kps = KPS3(kcu)
end
res1 = zeros(SVector{SEGMENTS+1, KiteModels.KVec3})
res2 = deepcopy(res1)
if ! @isdefined res3
    const res3 = reduce(vcat, vcat(res1, res2))
end

@testset verbose = true "KPS3 tests...." begin

function set_defaults()
    KiteModels.clear!(kps)
    kps.set.l_tether = 150.0
    kps.set.elevation = 60.0
    kps.set.area = 20.0
    kps.set.rel_side_area = 50.0
    kps.set.v_wind = 8.0
    kps.set.mass = 11.4
    kps.set.damping =  2 * 473.0
    kps.set.alpha = 1.0/7
    kps.set.c_s = 0.6
end

function init_392()
    KiteModels.clear!(kps)
    kps.set.l_tether = 392.0
    kps.set.elevation = 70.0
    kps.set.area = 10.0
    kps.set.rel_side_area = 50.0
    kps.set.v_wind = 9.1
    kps.set.mass = 6.2
    kps.set.c_s = 0.6
end

set_defaults()

@testset "calc_rho             " begin
    @test isapprox(calc_rho(kps.am, 0.0), 1.225, atol=1e-5)
    @test isapprox(calc_rho(kps.am, 100.0), 1.210756, atol=1e-5)
end

@testset "calc_wind_factor     " begin
    set_defaults()
    @test isapprox(calc_wind_factor(kps.am, 6.0, Int(EXP)),   1.0, atol=1e-5)
    @test isapprox(calc_wind_factor(kps.am, 10.0, Int(EXP)),  1.0757037, atol=1e-5)
    @test isapprox(calc_wind_factor(kps.am, 100.0, Int(EXP)), 1.494685, atol=1e-5)
end

@testset "calc_cl              " begin
    @test isapprox(kps.calc_cl(-5.0), 0.150002588978, atol=1e-4)
    @test isapprox(kps.calc_cl( 0.0), 0.200085035326, atol=1e-4)
    @test isapprox(kps.calc_cl(10.0), 0.574103590856, atol=1e-4)
    @test isapprox(kps.calc_cl(20.0), 1.0, atol=1e-4)
end

@testset "test_calc_drag       " begin
    v_segment = KVec3(1.0, 2, 3)
    unit_vector = KVec3(2.0, 3.0, 4.0)
    rho = SimFloat(calc_rho(kps.am, 10.0))
    last_tether_drag = KVec3(0.0, 0.0, 0.0)
    v_app_perp = KVec3(0, -3.0, -4.0)
    area = 20.0
    kps.v_wind_tether .= [0.1, 0.2, 0.3]
    v_app_norm = KiteModels.calc_drag(kps, v_segment, unit_vector, rho, v_app_perp, area)
    @test v_app_norm ≈ 3.3674916481
    @test kps.last_tether_drag ≈ [-38506.7140169,  -57266.39520463, -76026.07639235]
    @test v_app_perp ≈ [ 35.1, 52.2, 69.3]
end

@testset "test_calc_aero_forces" begin
    set_defaults()
    kps.v_apparent .= KVec3(35.1, 52.2, 69.3)
    kps.v_wind .= kps.v_wind_gnd
    pos_kite = KVec3(30.0, 5.0, 100.0)
    v_kite = KVec3(3.0, 5.0, 2.0)
    rho = SimFloat(calc_rho(kps.am, 10.0))
    rel_steering = 0.1
    kps.beta = 0.1
    kps.psi = 0.2
    kps.param_cl = 0.2
    kps.param_cd = 1.0
    KiteModels.calc_aero_forces(kps, pos_kite, v_kite, rho, rel_steering)
    @test kps.v_apparent ≈ [5.0,  -5, -2]                                # same with Python
    @test kps.kite_y ≈ [ 0.64101597,  0.73258967, -0.22893427]           # same with Python
    @test kps.cor_steering ≈ 0.0250173783309
    @test kps.steering_force ≈ [-15.88482337, -18.15408385, 5.6731512 ]
    @test kps.last_force ≈ [-555.24319976, 544.82004621, 80.49946362]    # same with Python
end

@testset "test_calc_res        " begin
    set_defaults()
    i = 2
    pos1 = KVec3(30.0, 5.0, 100.0)
    pos2 = KVec3(30.0+10, 5.0+11, 100.0+20)
    vel1 = KVec3(3.0, 5.0, 2.0)
    vel2 = KVec3(3.0+0.1, 5.0+0.2, 2.0+0.3)
    mass = 9.0
    veld = KVec3(0.1, 0.3, 0.4)
    result = KVec3(0, 0, 0)
    kps.c_spring = 0.011
    kps.damping = 0.01
    kps.last_tether_drag = KVec3(5.0,6,7)
    kps.last_force = KVec3(-1.0, -2, -3)
    kps.v_app_perp = KVec3(0.1,0.22,0.33)
    kps.v_wind_tether .= [0.1, 0.2, 0.3]
    kps.segment_length = 10.0
    KiteModels.calc_res(kps, pos1, pos2, vel1, vel2, mass, veld, result, i)
    @test result ≈ [ -6.99178740e-03, 1.01297086e-01, 9.83843908e+00]
    i = SEGMENTS+1
    KiteModels.calc_res(kps, pos1, pos2, vel1, vel2, mass, veld, result, i)
    @test result ≈ [0.14263997, 0.41669208, 10.12449937]
end

@testset "test_calc_loop       " begin
    set_defaults()
    kps.last_tether_drag = KVec3(5.0,6,7)
    kps.last_force = KVec3(-1.0, -2, -3)
    kps.v_app_perp = KVec3(0.1,0.22,0.33)
    kps.v_wind_tether .= [0.1, 0.2, 0.3]
    kps.segment_length = 10.0
    kps.c_spring = kps.set.c_spring / kps.segment_length
    kps.damping  = kps.set.damping / kps.segment_length
    pos  = zeros(SVector{SEGMENTS+1, KVec3})
    for i in 1:SEGMENTS+1
        pos[i][3] = 5.0 * (i-1)
    end
    vel  = zeros(SVector{SEGMENTS+1, KVec3})
    posd = zeros(SVector{SEGMENTS+1, KVec3})
    veld = zeros(SVector{SEGMENTS+1, KVec3})
    res1 = zeros(SVector{SEGMENTS+1, KVec3})
    res2 = zeros(SVector{SEGMENTS+1, KVec3})
    @test kps.c_spring ≈ 61460.0
    @test kps.damping  ≈    94.6
    KiteModels.loop(kps, pos, vel, posd, veld, res1, res2)
    @test sum(res1) ≈ [0.0, 0.0, 0.0]
    @test isapprox(res2[7], [ -5.02874357e-02, -1.00574871e-01, -7.62063430e+02])
    @test isapprox(res2[6], [ -6.15356097e-03, -1.23071219e-02, 9.81000000e+00])
    @test isapprox(res2[5], [ -2.38000593e-03, -4.76001187e-03, 9.81000000e+00])
    @test isapprox(res2[2], [ -2.38139816e-03, -4.76279631e-03, 9.81000000e+00], rtol=1e-5)
    @test isapprox(res2[1], [0.0,0.0,0.0], rtol=1e-5)
end

@testset "test_calc_alpha      " begin
    v_app = KVec3(10,2,3)
    vec_z = normalize(KVec3(3,2,0))
    alpha = KiteModels.calc_alpha(v_app, vec_z)
    @test alpha ≈ -1.091003745821884
end

@testset "test_set_cl_cd       " begin
    alpha = deg2rad(10.0)
    KiteModels.set_cl_cd!(kps, alpha)
    @test kps.param_cl ≈ 0.5740976324353215
    @test kps.param_cd ≈ 0.12533689264639403
end

@testset "test_calc_set_cl_cd  " begin
    kps.kcu.depower = 0.236
    kps.kcu.set_depower = kps.kcu.depower
    KiteModels.set_depower_steering!(kps, get_depower(kps.kcu), get_steering(kps.kcu))
    v_app = KVec3(10,2,3)
    vec_c = KVec3(3,2,0)
    KiteModels.calc_set_cl_cd!(kps, vec_c, v_app)
    @test kps.param_cl ≈ -0.09238222805380593
    @test kps.param_cd ≈ 0.8117345278100984
end

@testset "test_clear           " begin
    kps.t_0 = 10.0
    KiteModels.clear!(kps)
    @test kps.t_0 == 0.0
    @test kps.v_reel_out == 0.0
    @test kps.last_v_reel_out == 0.0
end

# Inputs:
# State vector state_y   = pos1, pos2, ..., posn, vel1, vel2, ..., veln
# Derivative   der_yd    = vel1, vel2, ..., veln, acc1, acc2, ..., accn
# Output:
# Residual     res = res1, res2 = pos1,  ..., vel1, ...
@testset "test_residual!       " begin
    res1 = zeros(SVector{SEGMENTS, KVec3})
    res2 = deepcopy(res1)
    res = reduce(vcat, vcat(res1, res2))
    X = zeros(SimFloat, 2*kps.set.segments)
    y0, yd0 = KiteModels.init(kps, X)
    # println(y0)
    # println(yd0)
    p = kps
    t = 0.0
    clear!(kps)
    residual!(res, yd0, y0, p, t)
    res1 = res[1:3*SEGMENTS]
    res2 = res[3*SEGMENTS+1:end]
    @test res1 == zeros(3*(SEGMENTS))
    # TODO: add test for res2
    # println(res2)
end

@testset "test_set_v_reel_out  " begin
    v_reel_out = 1.1
    t_0 = 5.5
    KiteModels.set_v_reel_out!(kps, v_reel_out, t_0)
    @test kps.v_reel_out ≈ 1.1
    @test kps.t_0 ≈ 5.5
    clear!(kps)
end

@testset "test_set_depower_steering" begin
    depower  = 0.25
    steering = 0.1
    KiteModels.set_depower_steering!(kps, depower, steering)
    @test kps.depower == depower
    @test kps.steering ≈ 0.09034121603653548
end

@testset "test_init            " begin
    my_state = deepcopy(kps)
    y0, yd0 = KiteModels.init(my_state, zeros(SimFloat, 2*SEGMENTS), delta=1e-6)
    @test length(y0)  == (SEGMENTS) * 6
    @test length(yd0) == (SEGMENTS) * 6
    @test sum(y0)  ≈ 717.163369868302
    @test sum(yd0) ≈ 3.6e-5
    @test isapprox(my_state.param_cl, 0.574103590856, atol=1e-4)
    @test isapprox(my_state.param_cd, 0.125342896308, atol=1e-4)
end

function test_initial_condition(params::Vector)
    my_state = kps
    y0, yd0 = KiteModels.init(my_state, params)
    residual!(res3, yd0, y0, kps, 0.0)
    return norm(res3) # z component of force on all particles but the first
end

res = nothing
x= nothing
z= nothing
@testset "test_initial_residual" begin
    # global res, x, z
    init_392()
    initial_x =  [-1.52505,  -3.67761,  -5.51761,  -6.08916,  -4.41371,  0.902124,  0.366393,  0.909132,  1.27537,  1.1538,  0.300657,  -1.51768]
    res = test_initial_condition(initial_x)

    my_state = kps
    kps.set.l_tether = 392.0
    kps.set.elevation = 70.0
    kps.set.area = 10.0
    kps.set.v_wind = 9.1
    kps.set.mass = 6.2
    KiteModels.clear!(my_state)
    alpha = deg2rad(10.0)
    KiteModels.set_cl_cd!(kps, alpha)

    # println("state.param_cl: $(my_state.param_cl), state.param_cd: $(my_state.param_cd)")
    # println("res2: "); display(my_state.res2)
    # println("pos: "); display(my_state.pos)
    x = Float64[]
    z = Float64[]
    for i in 1:length(my_state.pos)
        push!(x, my_state.pos[i][1])
        push!(z, my_state.pos[i][3])
    end
    # println(norm(res))

    @test my_state.segment_length ≈ 65.33333333333333
    @test my_state.c_spring ≈ 9407.142857142859
    @test my_state.damping  ≈  14.479591836734695
    # @test isapprox(my_state.param_cl, 1.0641931441572074, atol=1e-4)
    # @test isapprox(my_state.param_cd, 0.22825898470541978, atol=1e-4)
    @test sum(my_state.res1) ≈ [0.0, 0.0, 0.0]
    @test my_state.res2[1]   ≈ [0.0, 0.0, 0.0]

    # @test isapprox(my_state.res2[2], [8.83559075e+00, -4.72588546e-07, -5.10109289e+00], rtol=3e-2)
    # @test isapprox(my_state.res2[3], [8.81318565e+00, -4.68864292e-07, -5.08829453e+00], rtol=1e-3)
    # @test isapprox(my_state.res2[7], [1.49735632e+01,  2.71870215e-06,  4.51115984e+01], rtol=1e-3)
    # println("res2: "); display(my_state.res2)
    # println("lift force: $(norm(my_state.lift_force)) N")
end

@testset "test_getters" begin
    x, y, z = kite_ref_frame(kps)
    @test all(x .≈ [-0.9070010101306292, 0.0, 0.4211284455151642])
    @test all(y .≈ [0.0, 1.0, 0.0])
    @test all(z .≈ [-0.4211284455151642, -0.0, -0.9070010101306292])
    @test all(orient_euler(kps) .≈ [1.5707963267948966, -0.4346891114736793, 1.5707963267948966])
    @test all(pos_kite(kps) .≈ [134.97402018366216, 0.0, 366.8418273480761])
    @test calc_elevation(kps) .≈ 1.2182337959242815 # 69.8 deg
    @test calc_azimuth(kps) ≈ 0
    @test calc_heading(kps) ≈ 0
    calc_course(kps) # the course for vel_kite=zero is undefined, so we cannot test it
end

@testset "test_find_steady_state" begin
   KiteModels.set_depower_steering!(kps, 0.25, 0.0)
   res1, res2 = find_steady_state!(kps; delta=1e-6, prn=false)
   @test norm(res2) < 1e-5                            # velocity and acceleration must be near zero
   pre_tension = KiteModels.calc_pre_tension(kps)
   @test pre_tension > 1.0001
   @test pre_tension < 1.01
   @test unstretched_length(kps) ≈ 392.0              # initial, unstreched tether lenght
   @test tether_length(kps) ≈ 392.1861381318156 # real, streched tether length
   @test winch_force(kps) ≈ 276.25751212763817        # initial force at the winch [N]
   lift, drag = lift_drag(kps)
   @test lift ≈ 443.63277537186394                    # initial lift force of the kite [N]
   @test drag ≈ 94.25218065939362                     # initial drag force of the kite [N]
   @test lift_over_drag(kps) ≈ 4.706870146326417      # initial lift-over-drag
   @test norm(v_wind_kite(kps)) ≈ 9.107670173739065   # initial wind speed at the height of the kite [m/s]
end

function run_benchmarks()
    println("\ncalc_rho:")
    show(@benchmark calc_rho(height) setup=(height=1.0 + rand() * 200.0))
    println("\ncalc_wind_factor:")
    show(@benchmark calc_wind_factor(state.am, height) setup=(height=rand() * 200.0))
    println("\ncalc_cl:")
    show(@benchmark calc_cl(α) setup=(α=(rand()-0.5) * 360.0))
    println("\ncalc_drag:")
    show(@benchmark KiteModels.calc_drag(state, v_segment, unit_vector, rho, v_app_perp,
                            area) setup=(v_segment = KVec3(1.0, 2, 3);
                            unit_vector = KVec3(2.0, 3.0, 4.0);
                            rho = calc_rho(10.0f0); state.last_tether_drag = KVec3(0.0, 0.0, 0.0);
                            v_app_perp =  KVec3(0, -3.0, -4.0); area=se().area))
    println("\ncalc_aero_forces:")
    show(@benchmark KPS3.calc_aero_forces(state, pos_kite, v_kite, rho, rel_steering) setup=(state.v_apparent .= KVec3(35.1,
                                        52.2, 69.3); pos_kite = KVec3(30.0, 5.0, 100.0);
                                        v_kite = KVec3(3.0, 5.0, 2.0);
                                        rho = SimFloat(calc_rho(10.0));  rel_steering = 0.1))
    println("\ncalc_res:")
    show(@benchmark KPS3.calc_res(state, pos1, pos2, vel1, vel2, mass, veld, result, i) setup=(i = 1;
                            pos1 = KVec3(30.0, 5.0, 100.0); pos2 = KVec3(30.0+10, 5.0+11, 100.0+20);
                            vel1 = KVec3(3.0, 5.0, 2.0); vel2 = KVec3(3.0+0.1, 5.0+0.2, 2.0+0.3);
                            mass = 9.0; veld = KVec3(0.1, 0.3, 0.4); result = KVec3(0, 0, 0)))
    println("\ncalc_loop")
    show(@benchmark KPS3.loop(state, pos, vel, posd, veld, res1, res2) setup=(pos = zeros(SVector{SEGMENTS+1, KVec3});
                            vel  = zeros(SVector{SEGMENTS+1, KVec3}); posd  = zeros(SVector{SEGMENTS+1, KVec3});
                            veld  = zeros(SVector{SEGMENTS+1, KVec3}); res1  = zeros(SVector{SEGMENTS+1, KVec3});
                            res2  = zeros(SVector{SEGMENTS+1, KVec3}) ))
    println("\nset_cl_cd")
    show(@benchmark KPS3.set_cl_cd!(state, alpha) setup= (alpha = 10.0))
    println("\ncalc_alpha")
    show(@benchmark KPS3.calc_alpha(v_app, vec_z) setup=(v_app = KVec3(10,2,3); vec_z = normalize(KVec3(3,2,0))))
    println("\ncalc_set_cl_cd")
    show(@benchmark KPS3.calc_set_cl_cd!(state, vec_c, v_app) setup=(v_app = KVec3(10,2,3); vec_c = KVec3(3,2,0)))
    println("\nresidual!")
    show(@benchmark residual!(res, yd, y, p, t) setup = (res1 = zeros(SVector{SEGMENTS, KVec3}); res2 = deepcopy(res1);
                                                            res = reduce(vcat, vcat(res1, res2)); pos = deepcopy(res1);
                                                            pos[1] .= [1.0,2,3]; vel = deepcopy(res1); y = reduce(vcat, vcat(pos, vel));
                                                            der_pos = deepcopy(res1); der_vel = deepcopy(res1); yd = reduce(vcat, vcat(der_pos, der_vel));
                                                            p = SciMLBase.NullParameters(); t = 0.0))
    println("\ntest_initial_condition")
    show(@benchmark res=test_initial_condition(initial_x) setup = (initial_x = (zeros(12))))
    println()
end

end
nothing


text:
#: Gradient operator via convolution

export gradient_operator, gradient_operator_batch

# Generic convolutional operator

struct RealConvolutionalOperator{T,N}<:AbstractLinearOperator{T,N,N}
    cdims::DenseConvDims
    stencil::AbstractArray{<:Real,N}
end

AbstractLinearOperators.domain_size(A::RealConvolutionalOperator) = NNlib.input_size(A.cdims)
AbstractLinearOperators.range_size(A::RealConvolutionalOperator) = NNlib.output_size(A.cdims)
AbstractLinearOperators.matvecprod(A::RealConvolutionalOperator{CT,N}, u::AbstractArray{T,N}) where {T<:Real,N,CT<:RealOrComplex{T}} = conv(u, A.stencil)
AbstractLinearOperators.matvecprod(A::RealConvolutionalOperator{T,N}, u::AbstractArray{T,N}) where {T<:Complex,N} = matvecprod(A, real(u))+1im*matvecprod(A, imag(u))
AbstractLinearOperators.matvecprod_adj(A::RealConvolutionalOperator{CT,N}, v::AbstractArray{T,N}) where {T<:Real,N,CT<:RealOrComplex{T}} = ∇conv_data(v, A.stencil, A.cdims)
AbstractLinearOperators.matvecprod_adj(A::RealConvolutionalOperator{T,N}, v::AbstractArray{T,N}) where {T<:Complex,N} = matvecprod_adj(A, real(v))+1im*matvecprod_adj(A, imag(v))

Flux.gpu(A::RealConvolutionalOperator{T,N}) where {T,N} = RealConvolutionalOperator{T,N}(A.cdims, gpu(A.stencil))
Flux.cpu(A::RealConvolutionalOperator{T,N}) where {T,N} = RealConvolutionalOperator{T,N}(A.cdims, cpu(A.stencil))

# Gradient operator

function gradient_operator(n::NTuple{dim,Int64}, h::NTuple{dim,S}; T::DataType=Float32) where {dim,S<:Real}

    # Gradient operator (w/out shaping)
    stencil = gradient_stencil(real(T).(h))
    cdims = DenseConvDims((n...,1,1), size(stencil))
    ∇_ = RealConvolutionalOperator{T,dim+2}(cdims, stencil)

    # Reshaping operator
    Rin  = reshape_operator(T, n, (n...,1,1))
    Rout = reshape_operator(T, ((n.-1)...,dim,1), ((n.-1)...,dim))

    return Rout*∇_*Rin

end

function gradient_operator_batch(n::NTuple{dim,Int64}, nc::Int64, nb::Int64, h::NTuple{dim,S}; T::DataType=Float32) where {dim,S<:Real}

    # Gradient operator (w/out shaping)
    stencil = gradient_stencil(real(T).(h))
    cdims = DenseConvDims((n...,1,nc*nb), size(stencil))
    ∇_ = RealConvolutionalOperator{T,dim+2}(cdims, stencil)

    # Reshaping operator
    Rin  = reshape_operator(T, (n...,nc,nb), (n...,1,nc*nb))
    Rout = reshape_operator(T, ((n.-1)...,dim,nc*nb), ((n.-1)...,dim*nc,nb))

    return Rout*∇_*Rin

end

# Gradient stencil utils

function gradient_stencil(h::NTuple{dim,T}) where {dim,T}
    k = tuple(2*ones(Int64, dim)...)
    stencil = zeros(T,k...,1,dim)
    for i = 1:dim
        idx = Vector{UnitRange{Int64}}(undef,dim); fill!(idx, 2:2)
        idx[i] = 1:2
        view(stencil, tuple(idx...)...,1,i)[1:2] .= [T(1)/h[i]; T(-1)/h[i]]
    end
    return stencil
end

max_stars_repo_name: grizzuti/FastSolversForWeightedTV

max_issues_repo_name: grizzuti/FastSolversForWeightedTV

max_forks_repo_name: grizzuti/FastSolversForWeightedTV


text:
import Base.Cartesian.@ntuple

const ParticleArray = AbstractArray{<:AbstractParticles}
const SomeKindOfParticles = Union{<:MonteCarloMeasurements.AbstractParticles, ParticleArray}

nparticles(p) = length(p)
nparticles(p::ParticleArray) = length(eltype(p))
nparticles(p::Type{<:ParticleArray}) = length(eltype(p))

vecindex(p,i) = getindex(p,i)
vecindex(p::ParticleArray,i) = getindex.(p,i)
vecindex(p::NamedTuple,i) = (; Pair.(keys(p), ntuple(j->arggetter(i,p[j]), fieldcount(typeof(p))))...)

function indexof_particles(args)
    inds = findall([a <: SomeKindOfParticles for a in args])
    inds === nothing && throw(ArgumentError("At least one argument should be <: AbstractParticles. If particles appear nested as fields inside an argument, see `with_workspace` and `Workspace`"))
    all(nparticles(a) == nparticles(args[inds[1]]) for a in args[inds]) || throw(ArgumentError("All p::Particles must have the same number of particles."))
    (inds...,)
    # TODO: test all same number of particles
end

@generated function Ngetter(args...)
    nargs = length(args)
    inds = indexof_particles(args)
    N = nparticles(args[inds[1]])
    :($N)
end

function arggetter(i,a::Union{SomeKindOfParticles, NamedTuple})
    vecindex(a,i)
end

function arggetter(i,a)
    a
end

macro bymap(ex)
    @capture(ex, f_(args__)) || error("expected a function call")
    fsym = string(f)
    nargs = length(args)
    quote
        N = Ngetter($(esc.(args)...))
        individuals = map(1:N) do i
            eargs = ($(esc.(args)...),)
            argsi = ntuple(j->arggetter(i,eargs[j]), $nargs)
            $(esc(f))(argsi...)
        end
        if ndims(individuals[1]) == 0
            Particles(individuals)
        elseif ndims(individuals[1]) == 1
            Particles(copy(reduce(hcat,individuals)'))
        elseif ndims(individuals[1]) == 2
            reshape(Particles(copy(reduce(hcat,vec.(individuals))')), size(individuals[1])...)
        else
            error("Output with dimension >2 is currently not supported by `@bymap`. Consider if `ℝⁿ2ℝⁿ_function($($fsym), $($args...))` works for your use case.")
        end
    end
end

macro bypmap(ex)
    @capture(ex, f_(args__)) || error("expected a function call")
    fsym = string(f)
    nargs = length(args)
    quote
        N = Ngetter($(esc.(args)...))
        individuals = pmap(1:N) do i
            eargs = ($(esc.(args)...),)
            argsi = ntuple(j->arggetter(i,eargs[j]), $nargs)
            $(esc(f))(argsi...)
        end
        if ndims(individuals[1]) == 0
            Particles(individuals)
        elseif ndims(individuals[1]) == 1
            Particles(copy(reduce(hcat,individuals)'))
        elseif ndims(individuals[1]) == 2
            reshape(Particles(copy(reduce(hcat,vec.(individuals))')), size(individuals[1])...)
        else
            error("Output with dimension >2 is currently not supported by `@bymap`. Consider if `ℝⁿ2ℝⁿ_function($($fsym), $($args...))` works for your use case.")
        end
    end
end


text:
#=

  Euler #3 in Julia.

  Problem 3
  """
  The prime factors of 13195 are 5, 7, 13 and 29.
  What is the largest prime factor of the number 600851475143 ?
  """

  This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

include("Euler.jl")

# 0.00354527s
function euler3a()
    n = 600851475143
    max = 0
    while (n > 1)
        m = round(Int,ceil(sqrt(n)))
        for i in 2:m
            if (n % i == 0 && isPrime(i))
                max = i;
                n /= i;
                continue;
            end
        end
    end
    return max;

end

# 0.00607204s
function euler3b()
    return maximum(prime_divisors(600851475143))
end

# 0.00169180s
function euler3c()
    return maximum(factors(600851475143))
end

# run_euler(euler3a)
# run_euler(euler3b)
run_euler(euler3c)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
using BenchmarkTools

"""target area: x=217..240, y=-126..-69"""
function day17(tx1=217, ty1=-126, tx2=240, ty2=-69)
    part = [0, 0]
    maxtuple, yfail = Tuple{Int,Int,Int}[], min(ty1, ty2)
    for i in isqrt(2tx1):tx2, j in ty1:(-ty1 - 1)
        x, y, xvel, yvel, maxy, success = 0, 0, i, j, 0, false
        for t in 1:1000
            x += xvel
            y += yvel
            y > maxy && (maxy = y)
            xvel -= sign(xvel)
            yvel -= 1
            if tx1 <= x <= tx2 && ty1 <= y <= ty2
                if !success
                    push!(maxtuple, (maxy, i, j))
                    success = true
                end
            end
            y < yfail && break
        end
    end
    part[1] = maximum([first(t) for t in maxtuple])
    part[2] = length(maxtuple)
    return part
end

part = day17()
println("Part 1: ", part[1])
println("Part 2: ", part[2])

@btime day17()


text:
using ITensorTTN
using Test

@testset "ITensorTTN.jl" begin
  include(joinpath("..", "examples", "basics.jl"))
end


text:
#=

  Killer Sudoku in Julia ConstraintSolver.jl

  http://en.wikipedia.org/wiki/Killer_Sudoku
  """
  Killer sudoku (also killer su doku, sumdoku, sum doku, addoku, or
  samunamupure) is a puzzle that combines elements of sudoku and kakuro.
  Despite the name, the simpler killer sudokus can be easier to solve
  than regular sudokus, depending on the solver's skill at mental arithmetic;
  the hardest ones, however, can take hours to crack.

  ...

  The objective is to fill the grid with numbers from 1 to 9 in a way that
  the following conditions are met:

    * Each row, column, and nonet contains each number exactly once.
    * The sum of all numbers in a cage must match the small number printed
      in its corner.
    * No number appears more than once in a cage. (This is the standard rule
      for killer sudokus, and implies that no cage can include more
      than 9 cells.)

  In 'Killer X', an additional rule is that each of the long diagonals
  contains each number once.
  """

  Here we solve the problem from the Wikipedia page, also shown here
  http://en.wikipedia.org/wiki/File:Killersudoku_color.svg

  The output is:
    2 1 5 6 4 7 3 9 8
    3 6 8 9 5 2 1 7 4
    7 9 4 3 8 1 6 5 2
    5 8 6 2 7 4 9 3 1
    1 4 2 5 9 3 8 6 7
    9 7 3 8 1 6 4 2 5
    8 2 1 7 3 9 5 4 6
    6 5 9 4 2 8 7 1 3
    4 3 7 1 6 5 2 8 9

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

#
# Basic Sudoku model
#
function sudoku(model,x,n)
    # Row and column constraints
    for rc = 1:n
        @constraint(model, x[rc,:] in CS.AllDifferent())
        @constraint(model, x[:,rc] in CS.AllDifferent())
    end

    # Cell constraints
    ns = round(Int, sqrt(n))
    n2 = ns-1
    for i in 1:ns:n, j in 1:ns:n
        @constraint(model, vec([x[i+k,j+l] for k in 0:n2, l in 0:n2]) in CS.AllDifferent())
    end

end

function killer_sudoku(problem,print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>true,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            "activity.decay" => 0.999, # default 0.999
                                                            "activity.max_probes" => 100, # default, 10
                                                            "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))
    n = 9
    # define the variables
    @variable(model, 1 <= x[1:n,1:n] <= n, Int)

    sudoku(model, x, n) # plain Sudoku model

    # The Killer Sudoku hints
    for hint in problem
        s,lst = hint
        xs = [x[t[1],t[2]] for t in lst]
        @constraint(model,s == sum(xs))
        @constraint(model, xs in CS.AllDifferent())
    end

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                display(x_val)

            end
        end
    else
        println("status:$status")
    end

    return status
end

#
# hints
#
killer_sudoku_problems = Dict(
   :1 => [
            #  [Sum, [list of indices in X]]
            [ 3, [[1,1], [1,2]]],
            [15, [[1,3], [1,4], [1,5]]],
            [22, [[1,6], [2,5], [2,6], [3,5]]],
            [ 4, [[1,7], [2,7]]],
            [16, [[1,8], [2,8]]],
            [15, [[1,9], [2,9], [3,9], [4,9]]],
            [25, [[2,1], [2,2], [3,1], [3,2]]],
            [17, [[2,3], [2,4]]],
            [ 9, [[3,3], [3,4], [4,4]]],
            [ 8, [[3,6], [4,6],[5,6]]],
            [20, [[3,7], [3,8],[4,7]]],
            [ 6, [[4,1], [5,1]]],
            [14, [[4,2], [4,3]]],
            [17, [[4,5], [5,5],[6,5]]],
            [17, [[4,8], [5,7],[5,8]]],
            [13, [[5,2], [5,3],[6,2]]],
            [20, [[5,4], [6,4],[7,4]]],
            [12, [[5,9], [6,9]]],
            [27, [[6,1], [7,1],[8,1],[9,1]]],
            [ 6, [[6,3], [7,2],[7,3]]],
            [20, [[6,6], [7,6], [7,7]]],
            [ 6, [[6,7], [6,8]]],
            [10, [[7,5], [8,4],[8,5],[9,4]]],
            [14, [[7,8], [7,9],[8,8],[8,9]]],
            [ 8, [[8,2], [9,2]]],
            [16, [[8,3], [9,3]]],
            [15, [[8,6], [8,7]]],
            [13, [[9,5], [9,6],[9,7]]],
            [17, [[9,8], [9,9]]]
            ]

)

@time killer_sudoku(killer_sudoku_problems[:1])

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
"""
    $(SIGNATURES)

Draw the upward belief from clique `cllb::Symbol` message on variable `lb::Symbol`.

Example:

```julia
plotUpMsgsAtCliq(tree, :x2, :x1)
```

Related

plotKDE, getUpMsgs
"""
function plotUpMsgsAtCliq(treel::AbstractBayesTree,
                          cllb::Symbol,
                          lb::Symbol;
                          show::Bool=true,
                          w=20cm, h=15cm,
                          levels::Int=1,
                          dims::Union{Vector{Int}, Nothing}=nothing )
  #

  cliq = getClique(treel, cllb)
  cliqoutmsg = getUpMsgs(cliq)
  bel = convert(BallTreeDensity, cliqoutmsg.belief[lb])
  plotKDE(bel, dims=dims)
end

function plotTreeUpwardMsgs(fgl::G,
                            bt::AbstractBayesTree;
                            N=300 ) where G <: AbstractDFG
  #
  len = length(bt.cliques)-1
  vv = Array{Gadfly.Compose.Context,1}(undef, len)
  #r = Array{RemoteRef,1}(len)
  i = 0
  for cliq in bt.cliques
      if cliq[1] == 1 println("No upward msgs from root."); continue; end
      @show cliq[2].attributes["label"]
      i+=1
      vv[i] = drawHorDens(fgl, cliq[2].attributes["debug"].outmsg.p, N)
  end
  vv
end

"""
    $SIGNATURES

Project (convolve) to and take product of variable in Bayes/Junction tree.

Notes
- assume cliq and var sym the same, unless both specified.
- `cliqsym` defines a frontal variable of a clique.
"""
function plotTreeProductUp(fgl::G,
                           treel::AbstractBayesTree,
                           cliqsym::Symbol,
                           varsym::Symbol=cliqsym;
                           levels::Int=1,
                           dims::Vector{Int}=Int[]  ) where G <: AbstractDFG
  #
  # build a subgraph copy of clique
  cliq = getClique(treel, cliqsym)
  syms = getCliqAllVarIds(cliq)
  subfg = buildSubgraph(fgl, syms)

  # add upward messages to subgraph
  msgs = fetchMsgsUpChildren(treel,cliq, TreeBelief)
  # @show typeof(msgs)
  addMsgFactors!.(subfg, msgs, IIF.UpwardPass)

  # predictBelief
  cllbl = cliq.attributes["label"]
  return plotLocalProduct(subfg, varsym, title="Tree Up $(cllbl) | ", levels=levels, dims=dims)
end

function plotTreeProductDown( fgl::AbstractDFG,
                              treel::AbstractBayesTree,
                              cliqsym::Symbol,
                              varsym::Symbol=cliqsym;
                              levels::Int=1  )
  #
  # build a subgraph copy of clique
  cliq = whichCliq(treel, cliqsym)
  syms = getCliqAllVarIds(cliq)
  subfg = buildSubgraphFromLabels!(fgl,syms)

  # add upward messages to subgraph
  msgs = getCliqParentMsgDown(treel,cliq)
  addMsgFactors!(subfg, msgs)

  # predictBelief
  cllbl = cliq.attributes["label"]
  return plotLocalProduct(subfg, varsym, title="Tree Dwn $(cllbl) | ", levels=levels)
end

"""
    $SIGNATURES

Overlay plot all upward messages from cliques.
"""
function plotCliqUpMsgs(fg::G,
                        tree::AbstractBayesTree,
                        sym::Symbol;
                        show::Bool=true,
                        dims::Vector{Int}=Int[],
                        levels::Int=1,
                        c=nothing,
                        title="up msgs on $(sym)"    ) where G <: AbstractDFG
  #
  # get all msgs
  allmsgs = getTreeCliqUpMsgsAll(tree)
  # stack messages by variable
  sckmsgs = stackCliqUpMsgsByVariable(tree, allmsgs)

  if !haskey(sckmsgs, sym)
    @warn "plotCliqUpMsgs -- tree does not have up messages for $sym."
    return nothing
  end

  # prepend current estimate too
  Xs = getBelief(fg, sym)
  # vectorize beliefs
  beliefs = BallTreeDensity[Xs;]
  lbls = String["curr,-1";]
  for msg in sckmsgs[sym]
    push!(beliefs, msg.belief)
    push!(lbls, "$(msg.cliqId),$(msg.depth)")
  end

  # ignoring legend and color information
  cc = c === nothing ? getColorsByLength(length(beliefs)) : c

  # plot and return
  plotKDE(beliefs, levels=levels, c=cc, title=title, legend=lbls)
end

"""
    $SIGNATURES

Plot the downward messages currently stored in a clique.
"""
function plotCliqDownMsgs(tree::AbstractBayesTree,
                          frnt::Symbol;
                          show::Bool=true,
                          levels::Int=2,
                          dims=nothing,
                          existing=nothing  )
  #
  cliq = getClique(tree,frnt)
  msgs = IIF.getMessageBuffer(cliq).downRx.belief

  PL = []

  for (key, beldim) in msgs
    npl = plotKDE(manikde!(beldim.varType, beldim.val), levels=levels, title="dwn msg $key", dims=dims)
    existing === nothing ? nothing : union!(npl.layers, existing.layers)
    push!(PL, npl)
  end

  existing === nothing ? nothing : push!(PL, existing)
  pl = vstack(PL...)

  folderpath = "/tmp/caesar/random/"
  filepath = folderpath*"downmsgs_cliq$(cliq.index).pdf"
  Base.mkpath(folderpath)
  pl |> PDF(filepath, 20cm, length(PL)*12cm)
  @async run(`evince $filepath`)
  pl
end

#


text:
# For each angular frequency ω there are a multitude of effective wavenumbers.
using EffectiveWaves

medium = Medium(ρ=1.0, c=1.0)

species = [
    Specie(ρ=0.2, r=1.0, c=0.2, volfrac=0.25)
];

ω = 0.6;
num_wavenumbers = 8; # calculate the 8 least attenuating effective wavenumbers

# often more wavenumbers a returned then asked for.
k_effs = wavenumbers(ω, medium, species; num_wavenumbers=num_wavenumbers)

using Plots
scatter(k_effs, lab="Effective wavenumbers")

# savefig("many_keffs.png")

# We can calculate the resulting effective field form each of these effective wavenumbers but just choosing an amplitude of 1 for each wave.

waves = [EffectiveWave(ω, k_eff, medium, species; extinction_rescale=false) for k_eff in k_effs]

# choose the mesh the size of the lowest attenuating effective wavenumber
x = range(0.; length=100, stop=2pi/abs(k_effs[1]))

plot(x,waves[1], xlabel = "x", hankel_indexes = 0:1)
plot!(x,waves[2], hankel_indexes = 0:1)
plot!(x,waves[3], hankel_indexes = 0:1)

savefig("many_fields.png")


text:
Array{T}(::UndefInitializer, ::Tuple{PosInfinity}) where T = throw(ArgumentError("Cannot create infinite Array"))
Array{T}(::UndefInitializer, ::Tuple{Integer, PosInfinity}) where {T,N} = throw(ArgumentError("Cannot create infinite Array"))
Array{T}(::UndefInitializer, ::Tuple{PosInfinity, Integer}) where {T,N} = throw(ArgumentError("Cannot create infinite Array"))
Matrix{T}(::UndefInitializer, ::Tuple{Integer, PosInfinity}) where T = throw(ArgumentError("Cannot create infinite Array"))
Matrix{T}(::UndefInitializer, ::Tuple{PosInfinity, Integer}) where T = throw(ArgumentError("Cannot create infinite Array"))
Array{T}(::UndefInitializer, ::Tuple{PosInfinity, PosInfinity}) where {T,N} = throw(ArgumentError("Cannot create infinite Array"))
Matrix{T}(::UndefInitializer, ::Tuple{PosInfinity, PosInfinity}) where T = throw(ArgumentError("Cannot create infinite Array"))

Array{T}(::UndefInitializer, ::PosInfinity) where T = throw(ArgumentError("Cannot create infinite Array"))

Array{T}(::UndefInitializer, ::PosInfinity, ::PosInfinity) where T = throw(ArgumentError("Cannot create infinite Array"))
Array{T}(::UndefInitializer, ::PosInfinity, ::Integer) where T = throw(ArgumentError("Cannot create infinite Array"))
Array{T}(::UndefInitializer, ::Integer, ::PosInfinity) where T = throw(ArgumentError("Cannot create infinite Array"))

Matrix{T}(::UndefInitializer, ::PosInfinity, ::PosInfinity) where T = throw(ArgumentError("Cannot create infinite Array"))
Matrix{T}(::UndefInitializer, ::PosInfinity, ::Integer) where T = throw(ArgumentError("Cannot create infinite Array"))
Matrix{T}(::UndefInitializer, ::Integer, ::PosInfinity) where T = throw(ArgumentError("Cannot create infinite Array"))

Vector{T}(::UndefInitializer, ::Tuple{PosInfinity}) where T = throw(ArgumentError("Cannot create infinite Array"))
Vector{T}(::UndefInitializer, ::PosInfinity) where T = throw(ArgumentError("Cannot create infinite Array"))

similar(A::AbstractArray, ::Type{T}, axes::Tuple{OneToInf{Int}}) where T = cache(Zeros{T}(axes))
similar(A::AbstractArray, ::Type{T}, axes::Tuple{OneToInf{Int},OneToInf{Int}}) where T = cache(Zeros{T}(axes))
similar(A::AbstractArray, ::Type{T}, dims::Tuple{PosInfinity}) where T = cache(Zeros{T}(dims))
similar(A::AbstractArray, ::Type{T}, dims::Tuple{PosInfinity,PosInfinity}) where T = cache(Zeros{T}(dims))
similar(A::AbstractArray, ::Type{T}, dims::Tuple{Integer,PosInfinity}) where T = cache(Zeros{T}(dims))
similar(A::AbstractArray, ::Type{T}, dims::Tuple{PosInfinity,Integer}) where T = cache(Zeros{T}(dims))

similar(::Type{<:AbstractArray{T}}, axes::Tuple{OneToInf{Int}}) where T = cache(Zeros{T}(axes))
similar(::Type{<:AbstractArray{T}}, axes::Tuple{OneToInf{Int},OneToInf{Int}}) where T = cache(Zeros{T}(axes))
similar(::Type{<:AbstractArray{T}}, axes::Tuple{OneToInf{Int},OneTo{Int}}) where T = cache(Zeros{T}(axes))
similar(::Type{<:AbstractArray{T}}, axes::Tuple{OneTo{Int},OneToInf{Int}}) where T = cache(Zeros{T}(axes))
similar(::Type{<:AbstractArray{T}}, dims::Tuple{PosInfinity}) where T = cache(Zeros{T}(dims))
similar(::Type{<:AbstractArray{T}}, dims::Tuple{PosInfinity,PosInfinity}) where T = cache(Zeros{T}(dims))
similar(::Type{<:AbstractArray{T}}, dims::Tuple{Integer,PosInfinity}) where T = cache(Zeros{T}(dims))
similar(::Type{<:AbstractArray{T}}, dims::Tuple{PosInfinity,Integer}) where T = cache(Zeros{T}(dims))

similar(arr::AbstractArray{T}, sz::Union{Infinity,Integer,AbstractUnitRange}...) where T = similar(arr, Base.to_shape(sz)...)
similar(arr::AbstractArray, ::Type{T}, sz::Union{Infinity,Integer,AbstractUnitRange}...) where T = similar(arr, T, Base.to_shape(sz)...)

for Typ in (:Zeros, :Ones)
    @eval begin
        @inline $Typ{T, N}(sz::Tuple{Vararg{Union{Infinity,Integer}, N}}) where {T, N} = $Typ{T,N}(oneto.(sz))
        @inline $Typ{T, N}(sz::Vararg{Union{Infinity,Integer}, N}) where {T, N} = $Typ{T,N}(sz)
        @inline $Typ{T}(sz::Vararg{Union{Infinity,Integer},N}) where {T, N} = $Typ{T, N}(sz)
        @inline $Typ{T}(sz::Infinity) where T = $Typ{T,1}(sz)
    end
end

Eye{T}(::Infinity) where T = Eye{T}(ℵ₀)
Eye{T}(::Infinity, ::Infinity) where T = Eye{T}(ℵ₀, ℵ₀)
Eye{T}(::Infinity, n::Integer) where T = Eye{T}(ℵ₀, n)
Eye{T}(m::Integer, ::Infinity) where T = Eye{T}(m, ℵ₀)
Eye(::Infinity) = Eye(ℵ₀)
Eye(::Infinity, ::Infinity) = Eye(ℵ₀, ℵ₀)
Eye(::Infinity, n::Integer) = Eye(ℵ₀, n)
Eye(m::Integer, ::Infinity) = Eye(m, ℵ₀)

@inline Fill(x, sz::Union{Infinity,Integer}...) = Fill(x, Base.to_shape(sz)...)
@inline Fill{T}(x, sz::Union{Infinity,Integer}...) where T = Fill{T}(x, Base.to_shape(sz)...)

zeros(sz::Union{Infinity, Integer, AbstractUnitRange}...) = zeros(Base.to_shape(sz)...)
zeros(::Type{T}, sz::Union{Infinity, Integer, AbstractUnitRange}...) where T = zeros(T, Base.to_shape(sz)...)
zeros(::Type{T}, ::Tuple{PosInfinity}) where T = cache(Zeros{T}(∞))
zeros(::Type{T}, nm::Tuple{Integer, PosInfinity}) where T = cache(Zeros{T}(nm...))
zeros(::Type{T}, nm::Tuple{PosInfinity, Integer}) where T = cache(Zeros{T}(nm...))
zeros(::Type{T}, nm::Tuple{PosInfinity, PosInfinity}) where T = cache(Zeros{T}(nm...))

ones(sz::Union{Infinity, Integer, AbstractUnitRange}...) = ones(Base.to_shape(sz)...)
ones(::Type{T}, ::Tuple{PosInfinity}) where T = cache(Ones{T}(∞))
ones(::Type{T}, nm::Tuple{Integer, PosInfinity}) where T = cache(Ones{T}(nm...))
ones(::Type{T}, nm::Tuple{PosInfinity, Integer}) where T = cache(Ones{T}(nm...))
ones(::Type{T}, nm::Tuple{PosInfinity, PosInfinity}) where T = cache(Ones{T}(nm...))

fill(x, sz::Union{Infinity, Integer, AbstractUnitRange}...) = fill(x, Base.to_shape(sz)...)
fill(x, ::Tuple{PosInfinity}) = cache(Fill(x,∞))
fill(x, nm::Tuple{Integer, PosInfinity}) = cache(Fill(x,nm...))
fill(x, nm::Tuple{PosInfinity, Integer}) = cache(Fill(x,nm...))
fill(x, nm::Tuple{PosInfinity, PosInfinity}) = cache(Fill(x,nm...))

# This gets called when infinit number of columns
axes_print_matrix_row(_, io, X, A, i, ::AbstractVector{<:PosInfinity}, sep) = nothing
print_matrix_row(io::IO, X::AbstractVecOrMat, A::Vector, i::Integer, cols::AbstractVector{<:PosInfinity}, sep::AbstractString) = nothing
Base.print_matrix_row(io::IO,
        X::Union{LayoutMatrix,
        LayoutVector,
        AbstractTriangular{<:Any,<:LayoutMatrix},
        AdjOrTrans{<:Any,<:LayoutMatrix},
        AdjOrTrans{<:Any,<:LayoutVector},
        HermOrSym{<:Any,<:LayoutMatrix},
        SubArray{<:Any,2,<:LayoutMatrix},
        Diagonal{<:Any,<:LayoutVector}}, A::Vector,
        i::Integer, cols::AbstractVector{<:PosInfinity}, sep::AbstractString) =
    axes_print_matrix_row(axes(X), io, X, A, i, cols, sep)
Base.print_matrix_row(io::IO,
        X::Union{AbstractFill{<:Any,1},
                 AbstractFill{<:Any,2},
                 Diagonal{<:Any,<:AbstractFill{<:Any,1}},
                 RectDiagonal,
                 AbstractTriangular{<:Any,<:AbstractFill{<:Any,2}}
                 }, A::Vector,
        i::Integer, cols::AbstractVector{<:PosInfinity}, sep::AbstractString) =
    axes_print_matrix_row(axes(X), io, X, A, i, cols, sep)

print_matrix_vdots(io::IO, vdots::AbstractString,
        A::Vector, sep::AbstractString, M::Integer, ::NotANumber, pad_right::Bool = true) = nothing

# Avoid infinite loops on maximum
Base.mapreduce_impl(f, op, A::AbstractArray, ifirst::Integer, ::PosInfinity) =
    throw(ArgumentError("Cannot call mapreduce on an infinite length $(typeof(A))"))

function show_delim_array(io::IO, itr::AbstractArray, op, delim, cl,
                          delim_one, i1, ::PosInfinity)
    print(io, op)
    l = 20
    if !show_circular(io, itr)
        recur_io = IOContext(io, :SHOWN_SET => itr)
        if !haskey(io, :compact)
          recur_io = IOContext(recur_io, :compact => true)
        end
        first = true
        i = i1
        if 20 >= i1
          while true
              if !isassigned(itr, i)
                  print(io, undef_ref_str)
              else
                  x = itr[i]
                  show(recur_io, x)
              end
              i += 1
              if i > l
                  print(io, delim)
                  print(io, ' ')
                  print(io, '…')
                  delim_one && first && print(io, delim)
                  break
              end
              first = false
              print(io, delim)
              print(io, ' ')
          end
        end
    end
    print(io, cl)
end

#####
# FillArrays
#####

# Lazy Broadacasting
for typ in (:Ones, :Zeros, :Fill)
    @eval begin
        BroadcastStyle(::Type{$typ{T,N,NTuple{N,<:OneToInf}}}) where {T,N} = LazyArrayStyle{N}()
        BroadcastStyle(::Type{$typ{T,2,<:Tuple{<:Any,<:OneToInf}}}) where {T} = LazyArrayStyle{2}()
        BroadcastStyle(::Type{$typ{T,2,<:Tuple{<:OneToInf,<:Any}}}) where {T} = LazyArrayStyle{2}()
    end
end

BroadcastStyle(::Type{<:Diagonal{T,<:AbstractFill{T,1,Tuple{OneToInf{I}}}}}) where {T,I} = LazyArrayStyle{2}()

## Support broadcast(*, ::AbstractFill, A)

function copy(bc::Broadcasted{<:BroadcastStyle,<:Any,typeof(*),<:Tuple{Ones{T,1,Tuple{OneToInf{Int}}},AbstractArray{V,N}}}) where {N,T,V}
    a,b = bc.args
    @assert bc.axes == axes(b)
    convert(AbstractArray{promote_type(T,V),N}, b)
end

function copy(bc::Broadcasted{<:BroadcastStyle,<:Any,typeof(*),<:Tuple{AbstractArray{T,N},Ones{V,1,Tuple{OneToInf{Int}}}}}) where {N,T,V}
    a,b = bc.args
    @assert bc.axes == axes(a)
    convert(AbstractArray{promote_type(T,V),N}, a)
end

function copy(bc::Broadcasted{<:BroadcastStyle,<:Any,typeof(*),<:Tuple{AbstractFill{T,1,Tuple{OneToInf{Int}}},AbstractArray{V,N}}}) where {N,T<:Number,V}
    a,b = bc.args
    @assert bc.axes == axes(b)
    getindex_value(a) * b
end

function copy(bc::Broadcasted{<:BroadcastStyle,<:Any,typeof(*),<:Tuple{AbstractArray{T,N},AbstractFill{V,1,Tuple{OneToInf{Int}}}}}) where {N,T,V<:Number}
    a,b = bc.args
    @assert bc.axes == axes(a)
    a * getindex_value(b)
end

function copy(bc::Broadcasted{<:BroadcastStyle,<:Any,typeof(*),<:Tuple{AbstractFill{T,1,Tuple{OneToInf{Int}}},AbstractArray{V,N}}}) where {N,T,V}
    a,b = bc.args
    @assert bc.axes == axes(b)
    Ref(getindex_value(a)) .* b # Use broadcast in-case a is array-valued
end

function copy(bc::Broadcasted{<:BroadcastStyle,<:Any,typeof(*),<:Tuple{AbstractArray{T,N},AbstractFill{V,1,Tuple{OneToInf{Int}}}}}) where {N,T,V}
    a,b = bc.args
    @assert bc.axes == axes(a)
    a .* Ref(getindex_value(b)) # Use broadcast in-case b is array-valued
end

# row Vector
function copy(bc::Broadcasted{<:BroadcastStyle,<:Any,typeof(*),<:Tuple{Adjoint{T,Ones{T,1,Tuple{OneToInf{Int}}}},AbstractMatrix{V}}}) where {T,V}
    a,b = bc.args
    @assert bc.axes == axes(b)
    convert(AbstractMatrix{promote_type(T,V)}, b)
end

function copy(bc::Broadcasted{<:BroadcastStyle,<:Any,typeof(*),<:Tuple{AbstractMatrix{T},Adjoint{V,Ones{V,1,Tuple{OneToInf{Int}}}}}}) where {T,V}
    a,b = bc.args
    @assert bc.axes == axes(a)
    convert(AbstractMatrix{promote_type(T,V)}, a)
end

function copy(bc::Broadcasted{<:BroadcastStyle,<:Any,typeof(*),<:Tuple{Adjoint{T,<:AbstractFill{T,1,Tuple{OneToInf{Int}}}},AbstractMatrix{V}}}) where {T,V}
    a,b = bc.args
    @assert bc.axes == axes(b)
    getindex_value(a') * b
end

function copy(bc::Broadcasted{<:BroadcastStyle,<:Any,typeof(*),<:Tuple{AbstractMatrix{T},Adjoint{V,<:AbstractFill{V,1,Tuple{OneToInf{Int}}}}}}) where {T,V}
    a,b = bc.args
    @assert bc.axes == axes(a)
    a * getindex_value(b')
end

#####
# Diagonal
#####

one(D::Diagonal{T,<:AbstractFill{T,1,Tuple{OneToInf{Int}}}}) where T = Eye{T}(size(D,1))
copy(D::Diagonal{T,<:AbstractFill{T,1,Tuple{OneToInf{Int}}}}) where T = D

BroadcastStyle(::Type{<:Diagonal{<:Any,<:AbstractInfUnitRange}}) = LazyArrayStyle{2}()
sub_materialize(::AbstractBandedLayout, V, ::Tuple{InfAxes,InfAxes}) = V
sub_materialize(::AbstractBandedLayout, V, ::Tuple{OneTo{Int},InfAxes}) = V
sub_materialize(::AbstractBandedLayout, V, ::Tuple{InfAxes,OneTo{Int}}) = V

##
# banded columns are padded
##

sublayout(::DiagonalLayout{L}, ::Type{<:Tuple{KR,Integer}}) where {L,KR<:InfAxes} =
    sublayout(PaddedLayout{UnknownLayout}(), Tuple{KR})
sublayout(::DiagonalLayout{L}, ::Type{<:Tuple{Integer,JR}}) where {L,JR<:InfAxes} =
    sublayout(PaddedLayout{UnknownLayout}(), Tuple{JR})
-

sublayout(::AbstractBandedLayout, ::Type{<:Tuple{KR,Integer}}) where {KR<:InfAxes} =
    sublayout(PaddedLayout{UnknownLayout}(), Tuple{KR})
sublayout(::AbstractBandedLayout, ::Type{<:Tuple{Integer,JR}}) where {JR<:InfAxes} =
    sublayout(PaddedLayout{UnknownLayout}(), Tuple{JR})

function sub_paddeddata(::AbstractBandedLayout, S::SubArray{T,1,<:AbstractMatrix,<:Tuple{InfAxes,Integer}}) where T
    P = parent(S)
    (kr,j) = parentindices(S)
    view(P,first(kr):last(colsupport(P,j)),j)
end

function sub_paddeddata(::AbstractBandedLayout, S::SubArray{T,1,<:AbstractMatrix,<:Tuple{Integer,InfAxes}}) where T
    P = parent(S)
    (k,jr) = parentindices(S)
    view(P,k,first(jr):last(rowsupport(P,k)))
end

#####
# Concat: Hcat/Vcat
#####

function getindex(f::Vcat{T,1}, k::PosInfinity) where T
    length(f) == ℵ₀ || throw(BoundsError(f,k))
    ∞
end

_gettail(k, a::Number, b...) = k ≤ 1 ? tuple(a, b...) : _gettail(k - length(a), b...)
_gettail(k, a, b...) = k ≤ length(a) ? tuple(a[k:end], b...) : _gettail(k - length(a), b...)
_vcat(a) = a
_vcat(a, b, c...) = Vcat(a, b, c...)
getindex(A::Vcat, r::InfUnitRange) = Base.invoke(getindex, Tuple{AbstractArray, Any}, A, r)
_unsafe_getindex(::IndexLinear, A::Vcat, r::InfUnitRange) = _vcat(_gettail(first(r), A.args...)...)

# some common cases not catched by LayoutArrays + ambiguities
for InfColMatrix in (:(SubArray{<:Any,2,<:Any,<:Tuple{Any,InfIndexRanges}}),
                     :(SubArray{<:Any,2,<:LayoutVecOrMat,<:Tuple{Any,InfIndexRanges}}),
                     :(AbstractFill{<:Any,2,Tuple{OneTo{Int},OneToInf{Int}}}))
    @eval begin
        Base.typed_vcat(::Type{T}, A::$InfColMatrix, B::AbstractVecOrMat...) where T = Vcat{T}(A, B...)
        Base.typed_vcat(::Type{T}, A::$InfColMatrix, B::LayoutVecOrMats, C::AbstractVecOrMat...) where T = Vcat{T}(A, B, C...)
        Base.hcat(A::Number, B::$InfColMatrix) = Hcat(A, B)
        Base.typed_hcat(::Type{T}, A::AbstractVecOrMat, B::$InfColMatrix) where T = Hcat{T}(A, B)
    end
end

for InfRowArray in (:(AbstractFill{<:Any,1,Tuple{OneToInf{Int}}}),
                    :(AbstractFill{<:Any,2,Tuple{OneToInf{Int},OneTo{Int}}}),
                    :(SubArray{<:Any,2,<:Any,<:Tuple{InfIndexRanges,Any}}),
                    :(SubArray{<:Any,2,<:LayoutVecOrMat,<:Tuple{InfIndexRanges,Any}}),
                    :(SubArray{<:Any,1,<:Any,<:Tuple{Any,InfIndexRanges}}),
                    :(SubArray{<:Any,1,<:LayoutMatrix,<:Tuple{Any,InfIndexRanges}}),
                    :(SubArray{<:Any,1,<:Any,<:Tuple{InfIndexRanges,Any}}),
                    :(SubArray{<:Any,1,<:LayoutMatrix,<:Tuple{InfIndexRanges,Any}}))
    @eval begin
        Base.typed_hcat(::Type{T}, A::$InfRowArray, B::AbstractVecOrMat...) where T = Hcat{T}(A, B...)
        Base.typed_hcat(::Type{T}, A::$InfRowArray, B::LayoutVecOrMats, C::AbstractVecOrMat...) where T = Hcat{T}(A, B, C...)
        Base.typed_hcat(::Type{T}, A::AbstractVecOrMat, B::$InfRowArray, C::AbstractVecOrMat...) where T = Hcat{T}(A, B, C...)
        Base.typed_hcat(::Type{T}, A::LayoutVecOrMat, B::$InfRowArray, C::AbstractVecOrMat...) where T = Hcat{T}(A, B, C...)
    end
end

ArrayLayouts.typed_hcat(::Type{T}, ::Tuple{InfiniteCardinal{0},Any}, A...) where T = Hcat{T}(A...)
ArrayLayouts.typed_hcat(::Type{T}, ::Tuple{InfiniteCardinal{0},InfiniteCardinal{0}}, A...) where T = Hcat{T}(A...)
ArrayLayouts.typed_hcat(::Type{T}, ::Tuple{Any,InfiniteCardinal{0}}, A...) where T = Hcat{T}(A...)

ArrayLayouts.typed_vcat(::Type{T}, ::Tuple{InfiniteCardinal{0},Any}, A...) where T = Vcat{T}(A...)
ArrayLayouts.typed_vcat(::Type{T}, ::Tuple{InfiniteCardinal{0},InfiniteCardinal{0}}, A...) where T = Vcat{T}(A...)
ArrayLayouts.typed_vcat(::Type{T}, ::Tuple{Any,InfiniteCardinal{0}}, A...) where T = Vcat{T}(A...)

##
# lazy sub_materialize
##

sub_materialize(_, V, ::Tuple{InfAxes}) = V
sub_materialize(_, V, ::Tuple{InfAxes,InfAxes}) = V
sub_materialize(_, V, ::Tuple{<:Any,InfAxes}) = V
sub_materialize(_, V, ::Tuple{InfAxes,Any}) = V

sub_materialize(::ApplyLayout{typeof(vcat)}, V::AbstractVector, ::Tuple{InfAxes}) = ApplyArray(V)
sub_materialize(::ApplyLayout{typeof(vcat)}, V::AbstractMatrix, ::Tuple{InfAxes,InfAxes}) = ApplyArray(V)
sub_materialize(::ApplyLayout{typeof(vcat)}, V::AbstractMatrix, ::Tuple{<:Any,InfAxes}) = ApplyArray(V)
sub_materialize(::ApplyLayout{typeof(vcat)}, V::AbstractMatrix, ::Tuple{InfAxes,Any}) = ApplyArray(V)

sub_materialize(::ApplyLayout{typeof(hcat)}, V, ::Tuple{InfAxes}) = V
sub_materialize(::ApplyLayout{typeof(hcat)}, V::AbstractMatrix, ::Tuple{InfAxes,InfAxes}) = ApplyArray(V)
sub_materialize(::ApplyLayout{typeof(hcat)}, V::AbstractMatrix, ::Tuple{<:Any,InfAxes}) = ApplyArray(V)
sub_materialize(::ApplyLayout{typeof(hcat)}, V::AbstractMatrix, ::Tuple{InfAxes,Any}) = ApplyArray(V)

sub_materialize(::PaddedLayout, v::AbstractVector{T}, ::Tuple{InfAxes}) where T =
    _padded_sub_materialize(v)

Base._unsafe_getindex(::IndexStyle, A::AbstractVector, r::InfAxes) = layout_getindex(A, r)
Base._unsafe_getindex(::IndexStyle, A::AbstractFill{<:Any,1}, r::InfAxes) = FillArrays._fill_getindex(A, r)
getindex(A::AbstractCachedVector, r::InfAxes) = layout_getindex(A, r)

Base._unsafe_getindex(::IndexStyle, A::AbstractMatrix, kr::InfAxes, jr::InfAxes) = layout_getindex(A, kr, jr)
Base._unsafe_getindex(::IndexStyle, A::AbstractMatrix, kr::Union{Real, AbstractArray}, jr::InfAxes) = layout_getindex(A, kr, jr)
Base._unsafe_getindex(::IndexStyle, A::AbstractMatrix, kr::InfAxes, jr::Union{Real, AbstractArray}) = layout_getindex(A, kr, jr)

Base._unsafe_getindex(::IndexStyle, A::AbstractFill{<:Any,2}, kr::InfAxes, jr::InfAxes) = FillArrays._fill_getindex(A, kr, jr)
Base._unsafe_getindex(::IndexStyle, A::AbstractFill{<:Any,2}, kr::Union{Real, AbstractArray}, jr::InfAxes) = FillArrays._fill_getindex(A, kr, jr)
Base._unsafe_getindex(::IndexStyle, A::AbstractFill{<:Any,2}, kr::InfAxes, jr::Union{Real, AbstractArray}) = FillArrays._fill_getindex(A, kr, jr)

@inline getindex(A::ApplyMatrix{<:Any,typeof(hcat)}, kr::InfAxes, j::Integer) = layout_getindex(A, kr, j)

Base.checkindex(::Type{Bool}, inds::AbstractUnitRange, I::AbstractFill) = Base.checkindex(Bool, inds, getindex_value(I))
LazyArrays.cache_getindex(::InfiniteCardinal{0}, A::AbstractVector, I, J...) = layout_getindex(A, I, J...)


text:
"""
        set_kit(proj_name, proj_path)

    Create a folder named `proj_name` in `proj_path`, as well as a Pluto notebook and an internal file structure inside `proj_path/proj_name`.

    File structure and use of the notebook are detailed [here](https://github.com/FellowsFreiesWissen/computational_notebooks).

    See also: [`Pluto`](https://juliapackages.com/p/pluto)

    # Example
    ```julia-repl
     julia> include("set_kit.jl")
     julia> set_kit("example_project", pwd())
    ```
    """
function set_kit(proj_name::String, proj_path::String)
    ## Check project path and whether it is empty
    println(string("Will create a project at and set your working directory to: ", joinpath(proj_path, proj_name)))
    println("Is that correct? Answer yes or no:")
    user_ans = readline()
    if user_ans == "no"
        println("Check help('set_kit') to set up your project's location.")
    else
        if isdir(joinpath(proj_path, proj_name))
            println("There already is a project with this name in this path. Move it or change its name.")
        else
            ## Create project in path
            cd(proj_path)
            ## Create file structure
            ### first level
            mkdir(proj_name)
            ### main folders of a project
            map(mkdir, map(x -> joinpath(proj_name, x),
                           ["results","text", "submission"]))
            ### main folders of results
            map(mkdir, map(x -> joinpath(joinpath(proj_name,"results"), x),
                           ["data", "scripts"]))
            map(mkdir, map(x -> joinpath(joinpath(proj_name,"results", "data"), x),
                           ["raw", "process", "metadata"]))
            ### main folders of text
            map(mkdir, map(x -> joinpath(joinpath(proj_name,"text"), x),
                           ["figures","tables", "supplementary", "references"]))

            ### main folder
            touch(joinpath(proj_path, proj_name, "README.txt"))
            open(joinpath(proj_path, proj_name, "README.txt"), "w") do io
                write(io, "This folder contains the set up for a reproducible workflow as described by https://github.com/FellowsFreiesWissen/computational_notebooks.git\n")
                write(io, "\n")
                write(io, "The file structure is organized as such:\n")
                write(io, "\n")
                write(io, "projet_name\n")
                write(io, "|-- README.txt\n")
                write(io, "|-- main.Rmd\n")
                write(io, "|-- results\n")
                write(io, "|   |-- README.txt\n")
                write(io, "|   |-- data\n")
                write(io, "|   |   |-- raw\n")
                write(io, "|   |   |-- process\n")
                write(io, "|   |   |-- metadata\n")
                write(io, "|   |-- scripts\n")
                write(io, "|-- text\n")
                write(io, "|   |-- README.txt\n")
                write(io, "|   |-- main.doc\n")
                write(io, "|   |-- figures\n")
                write(io, "|   |-- tables\n")
                write(io, "|   |-- supplementary\n")
                write(io, "|   |-- references\n")
                write(io, "|-- submission\n")
                write(io, "|   |-- README.txt\n")
                write(io, "|   |-- journal1\n")
                write(io, "|       |-- first\n")
                write(io, "|   |-- journal2\n")
                write(io, "|       |-- first\n")
                write(io, "|       |-- revisions\n")
            end

            ### results folder
            touch(joinpath(proj_path, proj_name, "results/README.txt"))
            open(joinpath(proj_path, proj_name, "results/README.txt"), "w") do io
                write(io, "This folder contains all files of results or their processing, organized in the following subfolders:")
                write(io, "\n")
                write(io, "`data/raw`: your raw data files. These should not be protected against any change after the first storage")
                write(io, "\n")
                write(io, "`data/process`: these are files generated by processing the raw data and these are the ones used in the analysis. It can be a copy of the raw data if that is already ready to use. Any processing of the raw data to generate the files here should be documented in the notebook.")
                write(io, "\n")
                write(io, "`data/metadata`: files containing information about the data that will be useful for future users, readers, and reviewers of your data (e.g. description of variables names, units, and values)")
                write(io, "`scripts`: all code used to process the data, and which, for some reason or another is not included in the notebook because they too cumbersome or not of upmost relevancy for comprehension.")
            end

            ### text folder
            touch(joinpath(proj_path, proj_name, "text/README.txt"))
            open(joinpath(proj_path, proj_name, "text/README.txt"), "w") do io
                write(io, "This folder contains the main text of the manuscript, folders containing the figures and tables (unformatted) to be included in it, as well as a folder with the supplementary material and one with references.")
            end

            ### submission folder
            touch(joinpath(proj_path, proj_name, "submission/README.txt"))
            open(joinpath(proj_path, proj_name, "submission/README.txt"), "w") do io
                write(io, "This folder contains the files specific to journal submissions, e.g. cover letters, submitted versions.")
            end

            ## Create the minimal notebook
            touch(joinpath(proj_path, proj_name, string(proj_name, ".jl")))
            open(joinpath(proj_path, proj_name, string(proj_name, ".jl")), "w") do io
                write(io, "### A Pluto.jl notebook ###")
                write(io, "\n\n")
                write(io, "using Markdown")
                write(io, "\n")
                write(io, "using InteractiveUtils")
                write(io, "\n\n")
                write(io, "# ╔═╡ 307335ba-cf46-11eb-28b0-f199cf048ae6")
                write(io, "\n")
                write(io, "md\" ## Your title here\"")
                write(io, "\n\n")
                write(io, "# ╔═╡ dfdfa33f-53c1-45e2-8727-86ac01a09398")
                write(io, "\n")
                write(io, "md\"\"\"")
                write(io, "\n\n")
                write(io, "## Brief intro to Pluto notebooks and Markdown syntax")
                write(io, "\n\n")
                write(io, "Feel free to skip this tutorial if you already know how Pluto noteboos work. Also, make sure to delete it once you submit/share your notebook.")
                write(io, "\n\n")
                write(io, "This is your notebook. When filling it, you should obey Pluto's own syntax to create cells containing code or text.")
                write(io, "\n\n")
                write(io, "To write one single line of narrative text, use simple quotes and include \"md\" before the quotes (as done for the title header above). To write several lines of text, use three pairs of quotes, as in this cell.")
                write(io, "\n\n")
                write(io, "Text written with no special markings will appear without any special formatting when you generate a `.html` or `.pdf` version of this file (click the \"Export\" button above to convert it).")
                write(io, "\n\n")
                write(io, "It is possible to include code inside Markdown cells (it is not executed, though):")
                write(io, "\n\n")
                write(io, "- either as `inline code`")
                write(io, "\n\n")
                write(io, "- or as a block of code:")
                write(io, "\n\n")
                write(io, "```{julia}")
                write(io, "\n")
                write(io, "print(\"This is an example of a code block...\")")
                write(io, "\n")
                write(io, "print(\"for several lines of code\")")
                write(io, "\n")
                write(io, "```")
                write(io, "\n\n")
                write(io, "It is also possible to include hyperlinks: for example, this [link to further details on the Markdown syntax used in Pluto notebooks](https://www.juliapackages.com/p/pluto).")
                write(io, "\n\n")
                write(io, "Cells for code do not need any special marking, unless the cell contains several lines of code. In that case, the code should be included inside a `begin ... end` block.")
                write(io, "\n\n")
                write(io, "```{julia}")
                write(io, "\n")
                write(io, "begin")
                write(io, "\n")
                write(io, "    # your lines")
                write(io, "\n")
                write(io, "    # of code")
                write(io, "\n")
                write(io, "end")
                write(io, "\n")
                write(io, "```")
                write(io, "\n\n")
                write(io, "To execute the code in code cell and have its results appear above it, click the \"play\" (\"Run cell\") button below the cells or place your cursor inside the chunk and press `Shift+Enter`.")
                write(io, "\n\n")
                write(io, "The outputs of a cell are always show, unless the cell is disabled (available in the `...` - \"Actions\" - button on top of the cell). To hide the contents of the cell itself, click the eye icon on top of the cell.")
                write(io, "\n")
                write(io, "\"\"\"")
                write(io, "\n\n")
                write(io, "# ╔═╡ 37c7bffc-cf46-11eb-19b8-4175e38e818b")
                write(io, "\n")
                write(io, "begin")
                write(io, "\n")
                write(io, "# We suggest having a chunk dedicated to variables containing the paths to the folders related to the project.")
                write(io, "\n")
                write(io, "# The code in this chunk is not relevant for the reader, and thus is not included in the knitted version (therefore, `include = FALSE`).")
                write(io, "\n")
                write(io, "# This is also useful if you are not using the folder structure suggested by this kit, and want to preserve your privacy.")
                write(io, "\n")
                write(io, "# Feel free to edit these paths to adapt them to your needs or not use this suggestion at all.")
                write(io, "\n\n")
                write(io, "\tdata_dir = joinpath(\"results\", \"data\", \"process\")  ## Do NOT play with stuff in data/raw. That is your back-up. Work only on `process`.")
                write(io, "\n")
                write(io, "\tscripts_dir = joinpath(\"results\", \"scripts\")")
                write(io, "\n")
                write(io, "\tsuppl_dir = joinpath(\"results\", \"supplementary\")")
                write(io, "\n")
                write(io, "\tfigs_dir = joinpath(\"text\", \"figures\")")
                write(io, "\n")
                write(io, "\ttabs_dir = joinpath(\"text\", \"tables\")")
                write(io, "\n")
                write(io, "end")
                write(io, "\n\n")
                write(io, "# ╔═╡ f5062ac5-18e0-42fd-8bb7-1602a8298779")
                write(io, "\n")
                write(io, "md\"## Load data\"")
                write(io, "\n\n")
                write(io, "# ╔═╡ d4465a06-7d32-11ec-3907-59e1f4d60e2f")
                write(io, "\n")
                write(io, "# Here, you write the code to read the data from data_dir")
                write(io, "\n\n")
                write(io, "# ╔═╡ 6d256dbe-edc9-4132-aac7-62a473601034")
                write(io, "\n")
                write(io, "md\"\"\"")
                write(io, "\n")
                write(io, "## Data analysis")
                write(io, "\n\n")
                write(io, "### Figure 1")
                write(io, "\n")
                write(io, "\"\"\"")
                write(io, "\n\n")
                write(io, "# ╔═╡ f12401f7-1c08-4800-a64d-88f91f6b59fe")
                write(io, "\n")
                write(io, "# Here, you write the code to create a graph (Figure 1) included in the main text.")
                write(io, "\n")
                write(io, "# To avoid the repetition of having the same figure here as in the text, include a line that saves the figure in the figures directory.")
                write(io, "\n\n")
                write(io, "# ╔═╡ 2af762fe-0173-47eb-98d7-8fa6d3b555d3")
                write(io, "\n")
                write(io, "md\"### Figure S1\"")
                write(io, "\n\n")
                write(io, "# ╔═╡ 93c7698d-6e5b-4803-b7c6-90d4da7cd9bf")
                write(io, "\n")
                write(io, "# your code to plot Figure S1")
                write(io, "\n\n")
                write(io, "# ╔═╡ Cell order:")
                write(io, "\n")
                write(io, "# ╠═307335ba-cf46-11eb-28b0-f199cf048ae6")
                write(io, "\n")
                write(io, "# ╠═dfdfa33f-53c1-45e2-8727-86ac01a09398")
                write(io, "\n")
                write(io, "# ╠═37c7bffc-cf46-11eb-19b8-4175e38e818b")
                write(io, "\n")
                write(io, "# ╠═f5062ac5-18e0-42fd-8bb7-1602a8298779")
                write(io, "\n")
                write(io, "# ╠═d4465a06-7d32-11ec-3907-59e1f4d60e2f")
                write(io, "\n")
                write(io, "# ╠═6d256dbe-edc9-4132-aac7-62a473601034")
                write(io, "\n")
                write(io, "# ╠═f12401f7-1c08-4800-a64d-88f91f6b59fe")
                write(io, "\n")
                write(io, "# ╠═2af762fe-0173-47eb-98d7-8fa6d3b555d3")
                write(io, "\n")
                write(io, "# ╠═93c7698d-6e5b-4803-b7c6-90d4da7cd9bf")

                println("Your project is ready to go")
            end
        end
    end
end


text:
# This file is part of InvertedFiles.jl

export SortedIntSet

struct SortedIntSet{VecInt}
    set::VecInt
end

@inline Base.eltype(::SortedIntSet{VecInt}) where VecInt = eltype(VecInt)
@inline Base.size(s::SortedIntSet) = size(s.set)
@inline Base.length(s::SortedIntSet) = length(s.set)
@inline Base.eachindex(s::SortedIntSet) = eachindex(s.set)
@inline Base.getindex(s::SortedIntSet, i) = @inbounds s.set[i]
@inline Intersections._get_key(s::SortedIntSet, i) = @inbounds s.set[i]
@inline Base.first(s::SortedIntSet) = @inbounds s[1]
@inline Base.last(s::SortedIntSet) = @inbounds s[end]


text:
export bissect_secA
function bissect_secA(h :: AbstractNLPModel;
                     t1 :: Float64 = h.meta.x0[1],
                     t2 :: Float64 = 100.0,
                     tol :: Float64 = 1e-7,
                     maxiter :: Int = 50,
                     verbose :: Bool = false)

    (length(h.meta.x0) > 1) && warn("Not a 1-D problem ")
      (t0, t1) = trouve_intervalle(h, t1, 2.5)
      γ = 0.8; t = t1; tp = t0; tqnp = t0
      hk = 0.0; hkm1 = 0.0; gkm1 = 0.0; hplus = 0.0
      iter = 0

      hk = obj(h, [t])[1]; gk = grad(h, [t])[1]
      hkm1 = obj(h, [tqnp])[1]; gkm1 = grad(h, [tqnp])[1]

      verbose &&
        @printf(" iter        tp        t         dN         gk          ")
      verbose &&
        @printf(" gplus        \n")
      verbose &&
        @printf(" %7.2e %7.2e %7.2e  %7.2e  %7.2e  %7.2e\n",
                  iter, tp, t, 0.0, gk, 0.0)

      while ((abs(gk) > tol) && (iter < maxiter)) || (iter == 0)
          s = t - tqnp
          y = gk - gkm1

          Γ = 3 * (gk + gkm1) * s - 6 * (hk - hkm1)
          if y * s + Γ < eps(Float64) * (s^2)
            yt = y
          else
            yt = y + Γ / s
          end
          dN = -gk * s / yt

          if ((tp - t) * dN > 0) && (dN / (tp - t) < γ)
            tplus = t + dN
            hplus = obj(h, [tplus])[1]
            gplus = grad(h, [tplus])[1]
            verbose && println("N")
          else
            tplus = (t + tp) / 2
            hplus = obj(h, [tplus])[1]
            gplus = grad(h, [tplus])[1]
            verbose && println("B")
          end

          if t > tp
            if gplus < 0.0
              tp = t; tqnp = t; t = tplus
            else
              tqnp = t; t = tplus
            end
          else
            if gplus > 0.0
              tp = t; tqnp = t; t = tplus
            else
              tqnp = t; t = tplus
            end
          end

          #mise à jour des valeurs
          hkm1 = hk; gkm1 = gk; hk = hplus; gk = gplus
          iter = iter + 1

          verbose &&
            @printf(" %7.2e %7.2e %7.2e  %7.2e  %7.2e  %7.2e\n",
                    iter, tp, t, dN, gk, gplus)
        end

      topt = t
      status = :NotSolved
      (abs(gk) < tol) && (status = :Optimal)
      (iter >= maxiter) && (status = :Tired)
      tired = iter > maxiter
      optimal = abs(gk) < tol
      return (topt, obj(h, [topt])[1], norm(gk, Inf), iter, optimal, tired, status, h.counters.neval_obj, h.counters.neval_grad, h.counters.neval_hess)
end