text:
# A variant of SplitApplyCombine.groupfind using IdDict instead of Dictionaries.Dictionary
function _groupfind(container)
    T = keytype(container)
    inds = IdDict{eltype(container), Vector{T}}()
    @inbounds for i in keys(container)
        push!(get!(Vector{T}, inds, container[i]), i)
    end
    return inds
end

function _refs_pool(col::AbstractArray, reftype::Type{<:Integer}=UInt32)
    refs = refarray(col)
    pool = refpool(col)
    labeled = pool !== nothing
    if !labeled
        refs, invpool, pool = _label(col, eltype(col), reftype)
    end
    return refs, pool, labeled
end

# Obtain unique labels for row-wise pairs of values from a1 and a2 when mult is large enough
function _mult!(a1::AbstractArray, a2::AbstractArray, mult)
    z = zero(eltype(a1))
    @inbounds for i in eachindex(a1)
        x1 = a1[i]
        x2 = a2[i]
        # Handle missing values represented by zeros
        if iszero(x1) || iszero(x2)
            a1[i] = z
        else
            a1[i] += mult * (x2 - 1)
        end
	end
end

"""
    findcell(cols::VecColumnTable)
    findcell(names, data, esample=Colon())

Group the row indices of a collection of data columns
so that the combination of row values from these columns
are the same within each group.

Instead of directly providing the relevant portions of columns as
[`VecColumnTable`](@ref)``,
one may specify the `names` of columns from
`data` of any `Tables.jl`-compatible table type
over selected rows indicated by `esample`.
Note that unless `esample` covers all rows of `data`,
the row indices are those for the subsample selected based on `esample`
rather than those for the full `data`.

# Returns
- `IdDict{Tuple, Vector{Int}}`: a map from unique row values to row indices.
"""
function findcell(cols::VecColumnTable)
    ncol = size(cols, 2)
    isempty(cols) && throw(ArgumentError("no data column is found"))
    refs, pool, labeled = _refs_pool(cols[1])
    mult = length(pool)
    if ncol > 1
        # Make a copy to be used as cache
        labeled && (refs = copy(refs))
        @inbounds for n in 2:ncol
            refsn, pool, labeled = _refs_pool(cols[n])
            multn = length(pool)
            _mult!(refs, refsn, mult)
            mult = mult * multn
        end
    end
    return _groupfind(refs)
end

findcell(names, data, esample=Colon()) =
    findcell(subcolumns(data, names, esample))

"""
    cellrows(cols::VecColumnTable, refrows::IdDict)

A utility function for processing the object `refrows` returned by [`findcell`](@ref).
Unique row values from `cols` corresponding to
the keys in `refrows` are sorted lexicographically
and stored as rows in a new `VecColumnTable`.
Groups of row indices from the values of `refrows` are permuted to
match the order of row values and collected in a `Vector`.

# Returns
- `cells::VecColumnTable`: unique row values from columns in `cols`.
- `rows::Vector{Vector{Int}}`: row indices for each combination.
"""
function cellrows(cols::VecColumnTable, refrows::IdDict)
    isempty(refrows) && throw(ArgumentError("refrows is empty"))
    isempty(cols) && throw(ArgumentError("no data column is found"))
    ncol = length(cols)
    ncell = length(refrows)
    rows = Vector{Vector{Int}}(undef, ncell)
    columns = Vector{AbstractVector}(undef, ncol)
    for i in 1:ncol
        c = cols[i]
        if typeof(c) <: Union{ScaledArrOrSub, RotatingTimeArray}
            columns[i] = similar(c, ncell)
        else
            columns[i] = Vector{eltype(c)}(undef, ncell)
        end
    end
    refs = Vector{keytype(refrows)}(undef, ncell)
    r = 0
    @inbounds for (k, v) in refrows
        r += 1
        row1 = v[1]
        refs[r] = k
        for c in 1:ncol
            columns[c][r] = cols[c][row1]
        end
    end
    cells = VecColumnTable(columns, _names(cols), _lookup(cols))
    p = sortperm(cells)
    # Replace each column of cells with a new one in the sorted order
    @inbounds for i in 1:ncol
        columns[i] = cells[i][p]
    end
    # Collect rows in the same order as cells
    @inbounds for i in 1:ncell
        rows[i] = refrows[refs[p[i]]]
    end
    return cells, rows
end

"""
    settime(time::AbstractArray, step; start, stop, reftype, rotation)

Convert a column of time values to a [`ScaledArray`](@ref)
for representing discretized time periods of uniform length.
If `rotation` is specified (time values belong to multiple rotation groups),
a [`RotatingTimeArray`](@ref) is returned with the `time` field
being a [`ScaledArray`](@ref).
The returned array ensures well-defined time intervals for operations involving relative time
(such as [`lag`](@ref) and [`diff`](@ref)).
See also [`aligntime`](@ref).

# Arguments
- `time::AbstractArray`: the array containing time values.
- `step=nothing`: the length of each time interval; try `step=one(eltype(time))` if not specified.

# Keywords
- `start=nothing`: the first element of the `pool` of the returned [`ScaledArray`](@ref).
- `stop=nothing`: the last element of the `pool` of the returned [`ScaledArray`](@ref).
- `reftype::Type{<:Signed}=Int32`: the element type of the reference values for the [`ScaledArray`](@ref).
- `rotation=nothing`: rotation groups in a rotating sampling design.
"""
function settime(time::AbstractArray, step=nothing; start=nothing, stop=nothing,
        reftype::Type{<:Signed}=Int32, rotation=nothing)
    T = eltype(time)
    T <: ValidTimeType && !(T <: RotatingTimeValue) ||
        throw(ArgumentError("unaccepted element type $T from time column"))
    step === nothing && (step = one(T))
    time = ScaledArray(time, start, step, stop; reftype=reftype)
    rotation === nothing || (time = RotatingTimeArray(rotation, time))
    return time
end

"""
    aligntime(col::AbstractArray, time::ScaledArrOrSub)
    aligntime(col::AbstractArray, time::RotatingTimeArray)
    aligntime(data, colname::Union{Symbol,Integer}, timename::Union{Symbol,Integer})

Convert a column of time values `col` to a [`ScaledArray`](@ref) with a `pool`
that has the same first element and step size as the `pool` from
the [`ScaledArray`](@ref) `time`.
If `time` is a [`RotatingTimeArray`](@ref) with the `time` field being a [`ScaledArray`](@ref),
the returned array is also a [`RotatingTimeArray`](@ref)
with the `time` field being the converted [`ScaledArray`](@ref).
Alternative, the arrays may be specified with a `Tables.jl`-compatible `data` table
and column indices `colname` and `timename`.
See also [`settime`](@ref).

This is useful for representing all discretized time periods with the same scale
so that the underlying reference values returned by `DataAPI.refarray`
can be directly comparable across the columns.
"""
aligntime(col::AbstractArray, time::ScaledArrOrSub) = align(col, time)
aligntime(col::AbstractArray, time::RotatingTimeArray) =
    RotatingTimeArray(time.rotation, align(col, time.time))
aligntime(col::RotatingTimeArray, time::RotatingTimeArray) =
    RotatingTimeArray(time.rotation, align(col.time, time.time))

function aligntime(data, colname::Union{Symbol,Integer}, timename::Union{Symbol,Integer})
    checktable(data)
    return aligntime(getcolumn(data, colname), getcolumn(data, timename))
end

"""
    PanelStructure{R<:Signed, IP<:AbstractVector, TP<:AbstractVector}

Panel data structure defined by unique combinations of unit ids and time periods.
It contains the information required for certain operations such as
[`lag`](@ref) and [`diff`](@ref).
See also [`setpanel`](@ref).

# Fields
- `refs::Vector{R}`: reference values that allow obtaining time gaps by taking differences.
- `invrefs::Dict{R, Int}`: inverse map from `refs` to indices.
- `idpool::IP`: unique unit ids.
- `timepool::TP`: sorted unique time periods.
- `laginds::Dict{Int, Vector{Int}}`: a map from lag distances to vectors of indices of lagged values.
"""
struct PanelStructure{R<:Signed, IP<:AbstractVector, TP<:AbstractVector}
    refs::Vector{R}
    invrefs::Dict{R, Int}
    idpool::IP
    timepool::TP
    laginds::Dict{Int, Vector{Int}}
    function PanelStructure(refs::Vector{R}, idpool::IP, timepool::TP,
            laginds::Dict=Dict{Int, Vector{Int}}()) where {R,IP,TP}
        invrefs = Dict{R, Int}(ref=>i for (i, ref) in enumerate(refs))
        return new{R,IP,TP}(refs, invrefs, idpool, timepool, laginds)
    end
end

"""
    setpanel(data, idname, timename; step, reftype, rotation)
    setpanel(id::AbstractArray, time::AbstractArray; step, reftype, rotation)

Declare a [`PanelStructure`](@ref) which is required for certain operations
such as [`lag`](@ref) and [`diff`](@ref).
Unit IDs and time values can be provided either
as a table containing the relevant columns or as arrays.
`timestep` must be specified unless the `time` array is a [`ScaledArray`](@ref)
that is returned by [`settime`](@ref).

# Arguments
- `data`: a `Tables.jl`-compatible data table.
- `idname::Union{Symbol,Integer}`: the name of the column in `data` that contains unit IDs.
- `timename::Union{Symbol,Integer}`: the name of the column in `data` that contains time values.
- `id::AbstractArray`: the array containing unit IDs (only needed for the alternative method).
- `time::AbstractArray`: the array containing time values (only needed for the alternative method).

# Keywords
- `step=nothing`: the length of each time interval; try `step=one(eltype(time))` if not specified.
- `reftype::Type{<:Signed}=Int32`: the element type of the reference values for [`PanelStructure`](@ref).
- `rotation=nothing`: rotation groups in a rotating sampling design; use [`RotatingTimeValue`](@ref)s as reference values.

!!! note
    If the underlying data used to create the [`PanelStructure`](@ref) are modified.
    The changes will not be reflected in the existing instances of [`PanelStructure`](@ref).
    A new instance needs to be created with `setpanel`.
"""
function setpanel(id::AbstractArray, time::AbstractArray; step=nothing,
        reftype::Type{<:Signed}=Int32, rotation=nothing)
    eltype(time) <: ValidTimeType ||
        throw(ArgumentError("unaccepted element type $(eltype(time)) from time column"))
    length(id) == length(time) || throw(DimensionMismatch(
        "id has length $(length(id)) while time has length $(length(time))"))
    refs, idpool, labeled = _refs_pool(id)
    labeled && (refs = copy(refs); idpool = copy(idpool))
    time = settime(time, step; reftype=reftype, rotation=rotation)
    trefs = refarray(time)
    tpool = refpool(time)
    # Multiply 2 to create enough gaps between id groups for the largest possible lead/lag
    mult = 2 * length(tpool)
    _mult!(trefs, refs, mult)
    return PanelStructure(trefs, idpool, tpool)
end

function setpanel(data, idname::Union{Symbol,Integer}, timename::Union{Symbol,Integer};
        step=nothing, reftype::Type{<:Signed}=Int32)
    checktable(data)
    return setpanel(getcolumn(data, idname), getcolumn(data, timename);
        step=step, reftype=reftype)
end

show(io::IO, ::PanelStructure) = print(io, "Panel Structure")

function show(io::IO, ::MIME"text/plain", panel::PanelStructure)
    println(io, "Panel Structure:")
    println(IOContext(io, :limit=>true, :displaysize=>(1, 80)), "  idpool:   ", panel.idpool)
    println(IOContext(io, :limit=>true, :displaysize=>(1, 80)), "  timepool:   ", panel.timepool)
    print(IOContext(io, :limit=>true, :displaysize=>(1, 80)), "  laginds:  ", panel.laginds)
end

"""
    findlag!(panel::PanelStructure, l::Integer=1)

Construct a vector of indices of the `l`th lagged values
for all id-time combinations of `panel`
and save the result in `panel.laginds`.
If a lagged value does not exist, its index is filled with 0.
See also [`ilag!`](@ref).
"""
function findlag!(panel::PanelStructure, l::Integer=1)
    abs(l) < length(panel.timepool) ||
        throw(ArgumentError("|l| must be smaller than $(length(panel.timepool)); got $l"))
    refs = panel.refs
    invrefs = panel.invrefs
    T = eltype(refs)
    inds = Vector{Int}(undef, size(refs))
    l = convert(T, l)
    z = zero(T)
    @inbounds for i in keys(refs)
        ref = refs[i]
        inds[i] = get(invrefs, ref-l, z)
    end
    panel.laginds[l] = inds
    return inds
end

"""
    findlead!(panel::PanelStructure, l::Integer=1)

Construct a vector of indices of the `l`th lead values
for all id-time combinations of `panel`
and save the result in `panel.laginds`.
If a lead value does not exist, its index is filled with 0.
See also [`ilead!`](@ref).
"""
findlead!(panel::PanelStructure, l::Integer=1) = findlag!(panel, -l)

"""
    ilag!(panel::PanelStructure, l::Integer=1)

Return a vector of indices of the `l`th lagged values
for all id-time combinations of `panel`.
The indices are retrieved from [`panel`](@ref) if they have been collected before.
Otherwise, they are created by calling [`findlag!`](@ref).
See also [`ilead!`](@ref).
"""
function ilag!(panel::PanelStructure, l::Integer=1)
    il = get(panel.laginds, l, nothing)
    return il === nothing ? findlag!(panel, l) : il
end

"""
    ilead!(panel::PanelStructure, l::Integer=1)

Return a vector of indices of the `l`th lead values
for all id-time combinations of `panel`.
The indices are retrieved from [`panel`](@ref) if they have been collected before.
Otherwise, they are created by calling [`findlead!`](@ref).
See also [`ilag!`](@ref).
"""
ilead!(panel::PanelStructure, l::Integer=1) = ilag!(panel, -l)

"""
    lag(panel::PanelStructure, v::AbstractArray, l::Integer=1; default=missing)

Return a vector of `l`th lagged values of `v` with missing values filled with `default`.
The `panel` structure is respected.
See also [`ilag!`](@ref) and [`lead`](@ref).
"""
function lag(panel::PanelStructure, v::AbstractArray, l::Integer=1; default=missing)
    length(v) == length(panel.refs) || throw(DimensionMismatch(
        "v has length $(length(v)) while expecting $(length(panel.refs))"))
    inds = ilag!(panel, l)
    out = default === missing ? similar(v, Union{eltype(v), Missing}) : similar(v)
    @inbounds for i in 1:length(v)
        out[i] = inds[i] == 0 ? default : v[inds[i]]
    end
    return out
end

"""
    lead(panel::PanelStructure, v::AbstractArray, l::Integer=1; default=missing)

Return a vector of `l`th lead values of `v` with missing values filled with `default`.
The `panel` structure is respected.
See also [`ilead!`](@ref) and [`lag`](@ref).
"""
lead(panel::PanelStructure, v::AbstractArray, l::Integer=1; default=missing) =
    lag(panel, v, -l, default=default)

function _diff!(dest::AbstractArray, v::AbstractArray, inds::AbstractArray, default)
    @inbounds for i in 1:length(v)
        dest[i] = inds[i] == 0 ? default : v[i] - v[inds[i]]
    end
end

"""
    diff!(dest::AbstractArray, panel::PanelStructure, v::AbstractArray; kwargs...)

Take the differences of `v` within observations for each unit in `panel`
and store the result in `dest`.
By default, it calculates the first differences.
See also [`diff`](@ref).

# Keywords
- `order::Integer=1`: the order of differences to be taken.
- `l::Integer=1`: the time interval between each pair of observations.
- `default=missing`: default values for indices where the differences do not exist.
"""
function diff!(dest::AbstractArray, panel::PanelStructure, v::AbstractArray;
        order::Integer=1, l::Integer=1, default=missing)
    length(dest) == length(v) || throw(DimensionMismatch(
        "dest has length $(length(dest)) while v has length $(length(v))"))
    0 < order < length(panel.timepool) || throw(ArgumentError(
        "order must be between 0 and $(length(panel.timepool)); got $order"))
    inds = get(panel.laginds, l, nothing)
    inds === nothing && (inds = findlag!(panel, l))
    _diff!(dest, v, inds, default)
    if order > 1
        cache = similar(dest)
        for i in 2:order
            copy!(cache, dest)
            _diff!(dest, cache, inds, default)
        end
    end
    return dest
end

"""
    diff(panel::PanelStructure, v::AbstractArray; kwargs...)

Return the differences of `v` within observations for each unit in `panel`.
By default, it calculates the first differences.
See also [`diff!`](@ref).

# Keywords
- `order::Integer=1`: the order of differences to be taken.
- `l::Integer=1`: the time interval between each pair of observations.
- `default=missing`: default values for indices where the differences do not exist.
"""
function diff(panel::PanelStructure, v::AbstractArray;
        order::Integer=1, l::Integer=1, default=missing)
    out = default === missing ? similar(v, Union{eltype(v), Missing}) : similar(v)
    diff!(out, panel, v, order=order, l=l, default=default)
    return out
end


text:
@testset "StraightLineHomotopy" begin
    PolyImpl.@polyvar x y a
    H = StraightLineHomotopy([y^2-x], [x^2+3.0y])
    @test typeof(H)<:StraightLineHomotopy{Float64, PolySystem{Float64}, PolySystem{Float64}}
    @test_throws ErrorException StraightLineHomotopy([x^2+3.1y], [a])
    @test_throws ErrorException StraightLineHomotopy([x^2+3y, y^2-x], [x^2+3y])

    @test evaluate(H, [2.0, 1.0], 1.0) == [-1.0]
    @test H([2.0, 1.0], 1.0) == [-1.0]
    @test H([2.0, 1.0], 0.0) == [7.0]

    K = substitute(H, :y=>1.0)
    @test K([1.0], 0.0) == [4.0]
    @test K([1.0], 1.0) == [0.0]

    J_H = differentiate(H)

    @test J_H([1.0, 1.0], 0.0) == [2 3]
    @test J_H([1.0, 1.0], 1.0) == [-1 2]

    ∂H∂t = dt(H)

    ## time derivate is independent of t
    @test ∂H∂t([1.0,1.0], 0.23) == [-4]

    @test degrees(H) == [2]
    @test startsystem(H) == PolySystem([y^2-x])
    @test targetsystem(H) == PolySystem([x^2+3.0y])

    H = StraightLineHomotopy([y^2-x, 3x + y], [x^2+3y, 2y-3x])

    @test homogenized(H) == false
    K = homogenize(H)
    @test ishomogenous(H) == false
    @test homogenized(K)
    @test ishomogenous(K)
    @test nvariables(K) == 3

    @test length(removepoly(H, 2)) == 1
end

ext: jl


text:
function sqlite3_sql(stmtptr::Ptr{Void})
    @windows_only return ccall(
        (:sqlite3_sql, sqlite3_lib),
        stdcall,
        Ptr{Uint8},
        (Ptr{Void}, ),
        stmtptr
    )
    @unix_only return ccall(
        (:sqlite3_sql, sqlite3_lib),
        Ptr{Uint8},
        (Ptr{Void}, ),
        stmtptr
    )
end

function sqlite3_last_insert_rowid(dbptr::Ptr{Void})
    @windows_only return ccall(
        (:sqlite3_last_insert_rowid, sqlite3_lib),
        stdcall,
        Int64,
        (Ptr{Void}, ),
        dbptr
    )
    @unix_only return ccall(
        (:sqlite3_last_insert_rowid, sqlite3_lib),
        Int64,
        (Ptr{Void}, ),
        dbptr
    )
end


text:
using GLVisualize, GeometryTypes

if !isdefined(:runtests)
    window = glscreen()
    timesignal = loop(linspace(0f0,1f0,360))
end
description = """
Simple animated contour plot. You need to press Ctrl to move the camera.
Most 2D camera's use the Ctrl modifier to make editing easier.
"""

# use the performance tips to speed this up
# (http://docs.julialang.org/en/release-0.4/manual/performance-tips/)
# the array is 512x512 after all
const N = 256
const range = linspace(-5f0, 5f0, N)
const data = zeros(Intensity{Float32}, N, N)

function contour_inner(i, x, y)
    Intensity{Float32}(sin(1.3*x*i)*cos(0.9*y)+cos(.8*x)*sin(1.9*y)+cos(y*.2*x))
end

function contourdata(t)
    for i=1:size(data, 1)
        for j=1:size(data, 2)
            @inbounds data[i,j] = contour_inner(t, range[i], range[j])
        end
    end
    data
end

renderable = visualize(map(contourdata, timesignal), color_norm=Vec2f0(-3, 3))

_view(renderable, window, camera=:orthographic_pixel)

if !isdefined(:runtests)
    renderloop(window)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ d08a411c-f2c2-11ea-21a6-ad13627dfed8
function sum(a::Int, b::Int)
	a + b
end

# ╔═╡ e20da1b8-f2c2-11ea-2add-59f09ccd472b
sum(1, 2)

# ╔═╡ Cell order:
# ╠═d08a411c-f2c2-11ea-21a6-ad13627dfed8
# ╠═e20da1b8-f2c2-11ea-2add-59f09ccd472b


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

using Markdown
using InteractiveUtils

# ╔═╡ 79d92320-2412-11eb-14a7-b1b5246bbc3c
begin
	using Revise, ImageView, Images, FFTW
	using Plots, PhysicalOptics
end

# ╔═╡ 645e2eac-2433-11eb-0fed-6982080e4951
md" ## Single slit diffraction

The position of the diffraction peaks is:

$$x_k = (2k+1)\frac{z \lambda}{2d}$$

where $z$ is the distance from the slit to the screen, $\lambda$ the wavelength and $d$ the slit width. $\lambda$ is the wavelength.

The analytical solution of the field behind a single slit is in 1D:

$$I \sim \left( \frac{\sin\left(\frac{\pi\,d\,x}{\lambda \, z}\right)}{\left(\frac{\pi\,d\,x}{\lambda \, z}\right)} \right)^2$$

"

# ╔═╡ b0c8db42-2439-11eb-3d9b-b7af3c007f4c
begin
	# define the parameters and define slit
	M = 1024
	slit = zeros(Float64, (1, M))
	d = 100e-6
	L_slit = 10e-3
	z = 0.2
	λ = 550e-9
	xs = LinRange(-L_slit / 2, L_slit / 2, M)
	slit[1, (abs.(xs) .< d / 2)] .= 1
	plot(xs, slit[1, :])
end

# ╔═╡ ff9e9bf6-2445-11eb-0ca7-0f84a22d3c02
begin
	# propagate the initial field and take abs2 for intensity
	slit_prop = abs2.(propagate(slit, L_slit, z, kernel=fresnel_kernel))
	# normalize for a better plot
	slit_prop ./= maximum(slit_prop)

	analytical = begin
		x_0 = π .* d .* xs ./ λ ./ z
		(sin.(x_0) ./ x_0).^2
	end

	# plots results
	plot(xs, slit_prop[1, :])
	plot!(xs, analytical)
end

# ╔═╡ 3425c3a0-2440-11eb-08af-570a0aa4a373
md"## Diffraction pattern of a square aperture"

# ╔═╡ 340e8258-2440-11eb-2390-ffb2a54e889c
begin
	N = 256
	rect = zeros(Float64, (N, N))
	L_rect = 0.5
	w_rect = 0.051*2
	x_rect = LinRange(-L_rect/2, L_rect/2, N)
	for (j, y) in enumerate(x_rect)
		for (i, x) in enumerate(x_rect)
			if abs(x) <= w_rect /2 && abs(y) <= w_rect / 2
				rect[i, j] = 1.0
			end
		end
	end
end

# ╔═╡ c798a162-2444-11eb-1a76-ffe0d355671a
begin
	z_rect = 2e3
	# fresnel kernel
	out_rec = propagate(rect, L_rect, z_rect, kernel=fresnel_kernel, λ=500e-9)
	# full Rayleigh Sommerfeld Kernel
	out_recrs = propagate(rect, L_rect, z_rect, kernel=rs_kernel, λ=500e-9)

	imshow([abs2.(out_rec) abs2.(out_recrs) abs2.(rect)])
end

# ╔═╡ e34619f8-2444-11eb-3ce6-41e9921d4b0a
begin
	plot(x_rect, abs2.(out_rec[129, :]))
	plot!(x_rect, abs2.(rect[129, :]))
	plot!(x_rect, abs2.(out_recrs[129, :]))
end

# ╔═╡ Cell order:
# ╠═79d92320-2412-11eb-14a7-b1b5246bbc3c
# ╟─645e2eac-2433-11eb-0fed-6982080e4951
# ╠═b0c8db42-2439-11eb-3d9b-b7af3c007f4c
# ╠═ff9e9bf6-2445-11eb-0ca7-0f84a22d3c02
# ╟─3425c3a0-2440-11eb-08af-570a0aa4a373
# ╠═340e8258-2440-11eb-2390-ffb2a54e889c
# ╠═c798a162-2444-11eb-1a76-ffe0d355671a
# ╠═e34619f8-2444-11eb-3ce6-41e9921d4b0a

max_stars_repo_name: JuliaPhysics/PhysicalOptics

max_issues_repo_name: JuliaPhysics/PhysicalOptics

max_forks_repo_name: JuliaPhysics/PhysicalOptics


text:
using MLJ
using PrettyPrinting
import DataFrames
import Statistics
using PyPlot
using StableRNGs

LGBMRegressor = @load LGBMRegressor

features, targets = @load_boston
features = DataFrames.DataFrame(features)
@show size(features)
@show targets[1:3]
first(features, 3) |> pretty

DataFrames.describe(features)

train, test = partition(collect(eachindex(targets)), 0.70, shuffle=true,
                        rng=StableRNG(52))

lgb = LGBMRegressor() #initialised a model with default params
lgbm = machine(lgb, features[train, :], targets[train, 1])
boostrange = range(lgb, :num_iterations, lower=2, upper=500)
curve = learning_curve!(lgbm, resampling=CV(nfolds=5),
                        range=boostrange, resolution=100,
                        measure=rms)

figure(figsize=(8,6))
plot(curve.parameter_values, curve.measurements)
xlabel("Number of rounds", fontsize=14)
ylabel("RMSE", fontsize=14)

lgb = LGBMRegressor() #initialised a model with default params
lgbm = machine(lgb, features[train, :], targets[train, 1])
learning_range = range(lgb, :learning_rate, lower=1e-3, upper=1, scale=:log)
curve = learning_curve!(lgbm, resampling=CV(nfolds=5),
                        range=learning_range, resolution=100,
                        measure=rms)

figure(figsize=(8,6))
plot(curve.parameter_values, curve.measurements)
xscale("log")
xlabel("Learning rate (log scale)", fontsize=14)
ylabel("RMSE", fontsize=14)

lgb = LGBMRegressor() #initialised a model with default params
lgbm = machine(lgb, features[train, :], targets[train, 1])

leaf_range = range(lgb, :min_data_in_leaf, lower=1, upper=50)

curve = learning_curve!(lgbm, resampling=CV(nfolds=5),
                        range=leaf_range, resolution=50,
                        measure=rms)

figure(figsize=(8,6))
plot(curve.parameter_values, curve.measurements)
xlabel("Min data in leaf", fontsize=14)
ylabel("RMSE", fontsize=14)

r1 = range(lgb, :num_iterations, lower=50, upper=100)
r2 = range(lgb, :min_data_in_leaf, lower=2, upper=10)
r3 = range(lgb, :learning_rate, lower=1e-1, upper=1e0)
tm = TunedModel(model=lgb, tuning=Grid(resolution=5),
                resampling=CV(rng=StableRNG(123)), ranges=[r1,r2,r3],
                measure=rms)
mtm = machine(tm, features, targets)
fit!(mtm, rows=train);

best_model = fitted_params(mtm).best_model
@show best_model.learning_rate
@show best_model.min_data_in_leaf
@show best_model.num_iterations

predictions = MLJ.predict(mtm, rows=test)
rms_score = round(rms(predictions, targets[test, 1]), sigdigits=4)

@show rms_score

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
module TestModuleContext

using Test
using ..Utils: @getbool

const IN_CONTEXT = Ref(false)

function module_context(f)
    old = IN_CONTEXT[]
    IN_CONTEXT[] = @getbool("use_module_context")
    try
        f()
    finally
        IN_CONTEXT[] = old
    end
end

function test()
    @test IN_CONTEXT[]
end

end  # module


text:
struct GDALError <: Exception
    class::CPLErr
    code::Cint
    msg::String
end

function GDALError()
    class = getlasterrortype()
    code = getlasterrorno()
    msg = getlasterrormsg()
    GDALError(class, code, msg)
end

const throw_class = (CE_Failure, CE_Fatal)

function Base.showerror(io::IO, err::GDALError)
    err = string("GDALError (", err.class, ", code ", err.code, "):\n\t", err.msg)
    println(io, err)
end

function gdaljl_errorhandler(class::CPLErr, errno::Cint, errmsg::Cstring)
    # function signature needs to match the one in __init__, and the signature
    # of the callback for a custom error handler in the GDAL docs
    if class in throw_class
        throw(GDALError(class, errno, unsafe_string(errmsg)))
    end
    return C_NULL
end

function checknull(ptr)
    if ptr == C_NULL && getlasterrortype() in throw_class
        throw(GDALError())
    end
    ptr
end


text:
"""
    relative_energy_regain(ground_truth, rec)

Calculates the relative energy regain between the `ground_truth`
and the reconstruction.
Assumes that both arrays are 2 dimensional

# Reference
* [Rainer Heintzmann, \"Estimating missing information by maximum likelihood deconvolution\"](https://www.sciencedirect.com/science/article/abs/pii/S0968432806001272)
"""
function relative_energy_regain(ground_truth, rec)
    T = eltype(ground_truth)
    # go to fourier space
    ground_truth_fft = fft(ground_truth)
    rec_fft = fft(rec)

    # a dict to store the values for certain frequencies
    # we store a list since some (rounded) frequencies occur more than once
    ΔE_R_dict = Dict{T, Vector{T}}()
    E_R_dict = Dict{T, Vector{T}}()

    # round the frequencies to 4 digits, alternative would be to bin
    round4(x) = T(round(x, digits=3))

    # iterate over the frequencies and calculate the relative energy regain
    for (i₂, f₂) in enumerate(fftfreq(size(rec_fft, 2)))
        for (i₁, f₁) in enumerate(fftfreq(size(rec_fft, 1)))
            f_res = round4(√(f₁^2 + f₂^2))
            Δ_E_R = abs2(ground_truth_fft[i₁, i₂] - rec_fft[i₁, i₂])
            E_R = abs2(ground_truth_fft[i₁, i₂])

            update_dict_list!(ΔE_R_dict, f_res, Δ_E_R)
            update_dict_list!(E_R_dict, f_res, E_R)
        end
    end

    # finally transform everything into a list of frequencies and
    # a list of relative energy regains
    freqs = T[]
    G_R_list = T[]
    for f in sort(T.(keys(ΔE_R_dict)))
        push!(freqs, f)
        mean_ΔE_r = mean(ΔE_R_dict[f])
        mean_E_r = mean(E_R_dict[f])
        push!(G_R_list, (mean_E_r - mean_ΔE_r) / mean_E_r)
    end

    return freqs, G_R_list
end

"""
    update_dict_list!(d, k, v)

Updates the dict `d` which stores a list.
If `k` is in the keys of `d` we simply push `v` to the list
otherwise create a new list `[v]`
"""
function update_dict_list!(d, k, v)
    if haskey(d, k)
        push!(d[k], v)
    else
        d[k] = [v]
    end
    return d
end

"""
    normalized_cross_correlation(ground_truth, measured)

Calculates the normalized cross correlation.

External links:
* [Wikipedia](https://en.wikipedia.org/wiki/Sombrero_function)
* [StatsBase.jl](https://juliastats.org/StatsBase.jl/stable/signalcorr/#StatsBase.crosscor)
"""
function normalized_cross_correlation(ground_truth, measured)
    fl(x) = collect(Iterators.flatten(x))
    ground_truth = fl(ground_truth)
    measured = fl(measured)

    ncc = crosscor(ground_truth, measured, [0], demean=true)[begin]
    return ncc
end


text:
# uniform interface for model estimation

export Estimator, MLEstimator
export nsamples, estimate

abstract type Estimator{D<:Distribution} end

nsamples(e::Estimator{D}, x::Array) where {D<:UnivariateDistribution} = length(x)
nsamples(e::Estimator{D}, x::Matrix) where {D<:MultivariateDistribution} = size(x, 2)

mutable struct MLEstimator{D<:Distribution} <: Estimator{D} end
MLEstimator(::Type{D}) where {D<:Distribution} = MLEstimator{D}()

estimate(e::MLEstimator{D}, x) where {D<:Distribution} = fit_mle(D, x)
estimate(e::MLEstimator{D}, x, w) where {D<:Distribution} = fit_mle(D, x, w)


text:
using FDM, Test, Random, Printf, LinearAlgebra

@testset "FDM" begin
    include("methods.jl")
    include("numerics.jl")
    include("grad.jl")
end


text:
N = 51
model, obj = Dynamics.dubinscar
obj_con = ConstrainedObjective(obj, u_min=-10, u_max=10)
solver = Solver(model,obj_con,N=N)
solver.opts.verbose = true
n,m = get_sizes(solver)

U0 = rand(m,N)
res, stats = solve(solver, U0)
res.μ[1]

solver_warm = Solver(model,obj_con,N=N)
solver_warm.opts.verbose = true
solver_warm.opts.penalty_initial = res.μ[1][1]
λ = deepcopy(res.λ)
push!(λ, res.λN)
res_warm, stats_warm = _solve(solver_warm, to_array(res.U))
res_warm, stats_warm = _solve(solver_warm, to_array(res.U),λ=λ)

to_array(res.X) == to_array(res_warm.X)
to_array(res.U) == to_array(res_warm.U)
to_array(res.C) == to_array(res_warm.C)
to_array(res.λ) == to_array(res_warm.λ)
(res.CN) == (res_warm.CN)
(res.λN) == (res_warm.λN)

max_violation(res_warm)
max_violation(res)

to_array(res.λ)
p = plot()
plot_trajectory!(res)
cost(solver, res)

U0_warm = to_array(res.U)
res_warm = init_results(solver, X0, U0_warm, λ=λ)
rollout!(res_warm,solver)
cost(solver,res_warm)

to_array(res.μ)
to_array(res_warm.μ)

U = to_array(res.U)
solve(solver, U)
res.λ


text:
##
# cauchymoment implements the moments of the cauchy transform of a space
#  note that it starts from k==1
##

cauchymoment(S...) = stieltjesmoment(S...)/(-2π*im)

stieltjesmoment(S::WeightedJacobi{T,D},n::Int,z) where {T,D} = stieltjesjacobimoment(S.space.a,S.space.b,n,z)
stieltjesmoment(S::WeightedJacobi{T,D},z) where {T,D} = stieltjesjacobimoment(S.space.a,S.space.b,z)

stieltjesmoment(S::WeightedJacobiQ{T,D},n::Int,z) where {T,D} = stieltjesjacobimoment(S.space.a,S.space.b,n,z)
stieltjesmoment(S::WeightedJacobiQ{T,D},z) where {T,D} = stieltjesjacobimoment(S.space.a,S.space.b,z)

stieltjesmoment(S::Jacobi,n::Int,z) = stieltjesjacobimoment(S.a,S.b,n,z)
stieltjesmoment(S::Jacobi,z) = stieltjesjacobimoment(S.a,S.b,z)
stieltjesmoment(S::JacobiQ,n::Int,z) = stieltjesjacobimoment(S.a,S.b,n,z)
stieltjesmoment(S::JacobiQ,z) = stieltjesjacobimoment(S.a,S.b,z)

normalization(n::Int,α::Real,β::Real) = 2^(α+β)*gamma(n+α+1)*gamma(n+β+1)/gamma(2n+α+β+2)
stieltjesjacobimoment(α::Real,β::Real,n::Int,z) =
    (x = 2/(1-z);normalization(n,α,β)*HypergeometricFunctions.mxa_₂F₁(n+1,n+α+1,2n+α+β+2,x))
stieltjesjacobimoment(α::Real,β::Real,z) = stieltjesjacobimoment(α,β,0,z)

function logjacobimoment(α::Real,β::Real,n::Int,z)
    x = 2/(1-z)
    if n == 0
        2normalization(0,α,β)*(log(z-1)-dualpart(_₂F₁(dual(zero(α)+eps(α+β),one(β)),α+1,α+β+2,x)))
        # For testing purposes only, should be equivalent to above within radius of convergence
        # 2normalization(0,α,β)*(log(z-1)-(α+1)/(α+β+2)*x*_₃F₂(α+2,α+β+3,x))
    else
        -2normalization(n,α,β)/n*(-x)^n*_₂F₁(n,n+α+1,2n+α+β+2,x)
    end
end

logjacobimoment(α::Real,β::Real,z) = logjacobimoment(α,β,0,z)

function logabsjacobimoment(α::Real,β::Real,n::Int,z)
    x = 2/(1-z)
    if n == 0
        2normalization(0,α,β)*(logabs(z-1)-real(dualpart(_₂F₁(dual(zero(α)+eps(α+β),one(β)),α+1,α+β+2,x))))
        # For testing purposes only, should be equivalent to above within radius of convergence
        # 2normalization(0,α,β)*(logabs(z-1)-real((α+1)/(α+β+2)*x*_₃F₂(α+2,α+β+3,x)))
    else
        -2normalization(n,α,β)/n*real((-x)^n*_₂F₁(n,n+α+1,2n+α+β+2,x))
    end
end

logabsjacobimoment(α::Real,β::Real,z) = logabsjacobimoment(α,β,0,z)

stieltjeslegendremoment(n::Int,z) = stieltjesjacobimoment(zero(real(eltype(z))),zero(real(eltype(z))),n,z)
stieltjeslegendremoment(z) = stieltjeslegendremoment(0,z)

logabslegendremoment(z) = real(z)*logabs((z+1)/(z-1))-imag(z)*angle((z+1)/(z-1))+logabs(z^2-1)-2
                        # real(z*log((z+1)/(z-1)))+logabs(z^2-1)-2

#=
#TODO: this is for x^k but we probably want P_k

#These formulae are from mathematica
function stieltjesjacobimoment(α,β,k::Integer,z)
    if α == β == 0
        return stieltjeslegendremoment(k,z)
    elseif isapprox(α,0) && isapprox(β,0.5)
        if k==1
            return -2*sqrt(2)*(sqrtatansqrt(2/(z-1))-1)
        elseif k==2
            return 2*sqrt(2)*(1/3*(3z-2)-z*sqrtatansqrt(2/(z-1)))
        end
    elseif isapprox(α,0) && isapprox(β,-0.5)
        if k==1
            return 2sqrt(1/(z-1))*atan(sqrt(2/(z-1)))
        elseif k==2
            return z/sqrt(1-z)*log((1+1/sqrt(1-z))*(1-sqrt(2)/sqrt(1-z))/((1-1/sqrt(1-z))*(1+sqrt(2)/sqrt(1-z))))-
                2*sqrt(2)-2z/sqrt(1-z)*asinh(sqrt(-1/z))
        end
    elseif !isapprox(α,0) && isapprox(β,0.)
        return (-1)^k*stieltjesjacobimoment(β,α,k,-z)
    end
    error("stieltjesmoment not implemented for JacobiWeight "*string(α)*string(β))
end

stieltjesmoment{DD}(S::JacobiWeight{Chebyshev{DD},DD},k::Integer,z)=stieltjesjacobimoment(S.β,S.α,k,mobius(S,z))

function stieltjesmoment{T,DD}(S::JacobiWeight{Jacobi{T,DD},DD},k::Integer,z)
    z=mobius(S,z)
    J=S.space

    if k==1
        return stieltjesjacobimoment(S.β,S.α,k,z)
    elseif k==2
        if J.a==J.b
            return (1+J.a)*stieltjesjacobimoment(S.α,S.β,k,z)
        else
            return (J.a-J.b)/2*stieltjesjacobimoment(S.α,S.β,1,z) + (2+J.a+J.b)/2*stieltjesjacobimoment(S.α,S.β,2,z)
        end
    end
    error("stieltjesmoment not implemented for JacobiWeight "*string(S.α)*string(S.β))
end
=#

hilbertmoment(S::Space,k::Integer,x) = -real(stieltjesmoment(S,k,(x)⁺)+stieltjesmoment(S,k,(x)⁻))/(2π)


text:
@testset "Observables tests" begin
    a1 = (3/2, sqrt(3)/2)
    a2 = (3/2, -sqrt(3)/2)
    basis = [(0.0, 0.0), (1.0, 0.0)]
    uc = UnitCell(a1, a2)
    b1 = addBasisSite!(uc, basis[1])
    b2 = addBasisSite!(uc, basis[2])
    addInteraction!(uc, b1, b2, [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0], (0, 0))
    addInteraction!(uc, b1, b2, [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0], (0, -1))
    addInteraction!(uc, b1, b2, [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0], (-1, 0))
    setField!(uc, b1, [0.0, 0.0, -1.0])
    setField!(uc, b2, [0.0, 0.0, -1.0])
    lattice = Lattice(uc, (4,4))
    for i in 1:length(lattice)
        setSpin!(lattice, i, (0.0,0.0,1.0))
    end

    obs = Observables(lattice)
    energy = -42.0

    SpinMC.performMeasurements!(obs, lattice, energy)

    @test means(obs.energy)[1] ≈ energy / length(lattice)
    @test means(obs.energy)[2] ≈ (energy / length(lattice))^2
    @test mean(obs.magnetization) ≈ 1.0
    @test mean(obs.magnetizationVector) ≈ [0.0, 0.0, 1.0]
    @test all(mean(obs.correlation) .≈ ones(length(lattice)))
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""ResourceQuotaStatus defines the enforced hard limits and observed use.

    IoK8sApiCoreV1ResourceQuotaStatus(;
        hard=nothing,
        used=nothing,
    )

    - hard::Dict{String, IoK8sApimachineryPkgApiResourceQuantity} : Hard is the set of enforced hard limits for each named resource. More info: https://kubernetes.io/docs/concepts/policy/resource-quotas/
    - used::Dict{String, IoK8sApimachineryPkgApiResourceQuantity} : Used is the current observed total usage of the resource in the namespace.
"""
mutable struct IoK8sApiCoreV1ResourceQuotaStatus <: SwaggerModel
    hard::Any # spec type: Union{ Nothing, Dict{String, IoK8sApimachineryPkgApiResourceQuantity} } # spec name: hard
    used::Any # spec type: Union{ Nothing, Dict{String, IoK8sApimachineryPkgApiResourceQuantity} } # spec name: used

    function IoK8sApiCoreV1ResourceQuotaStatus(;hard=nothing, used=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1ResourceQuotaStatus, Symbol("hard"), hard)
        setfield!(o, Symbol("hard"), hard)
        validate_property(IoK8sApiCoreV1ResourceQuotaStatus, Symbol("used"), used)
        setfield!(o, Symbol("used"), used)
        o
    end
end # type IoK8sApiCoreV1ResourceQuotaStatus

const _property_map_IoK8sApiCoreV1ResourceQuotaStatus = Dict{Symbol,Symbol}(Symbol("hard")=>Symbol("hard"), Symbol("used")=>Symbol("used"))
const _property_types_IoK8sApiCoreV1ResourceQuotaStatus = Dict{Symbol,String}(Symbol("hard")=>"Dict{String, IoK8sApimachineryPkgApiResourceQuantity}", Symbol("used")=>"Dict{String, IoK8sApimachineryPkgApiResourceQuantity}")
Base.propertynames(::Type{ IoK8sApiCoreV1ResourceQuotaStatus }) = collect(keys(_property_map_IoK8sApiCoreV1ResourceQuotaStatus))
Swagger.property_type(::Type{ IoK8sApiCoreV1ResourceQuotaStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1ResourceQuotaStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ResourceQuotaStatus }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ResourceQuotaStatus[property_name]

function check_required(o::IoK8sApiCoreV1ResourceQuotaStatus)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1ResourceQuotaStatus }, name::Symbol, val)
end