text:
struct ElectricField{T, N, S} <: AbstractArray{T, N}
    data::Array{<:StaticArray{Tuple{N}, T}, N}
    grid::Grid{T, N, S}
end

@inline size(ep::ElectricField{T, N, S}) where {T, N, S} = size(ep.data)
@inline length(ep::ElectricField{T, N, S}) where {T, N, S} = length(ep.data)
@inline getindex(ep::ElectricField{T, N, S}, I::Vararg{Int, N}) where {T, N, S} = getindex(ep.data, I...)
@inline getindex(ep::ElectricField{T, N, S}, i::Int) where {T, N, S} = getindex(ep.data, i)
@inline getindex(ep::ElectricField{T, N, S}, s::Symbol) where {T, N, S} = getindex(ep.grid, s)

function NamedTuple(ep::ElectricField{T, 3}) where {T}
    return (
        grid = NamedTuple(ep.grid),
        values = ep.data #* internal_efield_unit,
    )
end
Base.convert(T::Type{NamedTuple}, x::ElectricField) = T(x)

function ElectricField(nt::NamedTuple)
    grid = Grid(nt.grid)
    T = eltype(ustrip.(nt.values[1]))
    S = get_coordinate_system(grid)
    N = get_number_of_dimensions(grid)
    ef = Array{SVector{3, T}}(undef, size(grid)...)
    @inbounds for i in eachindex(ef)
        # ef[i] = ustrip.(uconvert.(internal_efield_unit, nt.values[i]))
        ef[i] = nt.values[i]
    end
    ElectricField{T, N, S}( ef, grid)
end
Base.convert(T::Type{ElectricField}, x::NamedTuple) = T(x)

function ElectricField(ep::ElectricPotential{T, 3, S}, pointtypes::PointTypes{T}) where {T, S}
    return ElectricField{T, 3, S}(get_electric_field_from_potential( ep, pointtypes ), ep.grid)
end

function get_magnitude_of_rφz_vector(vector::AbstractArray,cutoff=NaN)
    magn=0
    magn+=(vector[1])^2
    magn+=(vector[3])^2
    result = sqrt(magn)
    if isnan(cutoff)
        return result
    else
        if result >cutoff
            result=0
        end
        return result
    end
end

function get_electric_field_from_potential(ep::ElectricPotential{T, 3, :cylindrical}, pointtypes::PointTypes{T}, fieldvector_coordinates=:xyz)::ElectricField{T, 3, :cylindrical} where {T <: SSDFloat}
    p = ep.data
    axr::Vector{T} = collect(ep.grid.r)
    axφ::Vector{T} = collect(ep.grid.φ)
    axz::Vector{T} = collect(ep.grid.z)

    cyclic::T = ep.grid.φ.interval.right
    ef = Array{SVector{3, T}}(undef, size(p)...)
    for iz in 1:size(ef, 3)
        for iφ in 1:size(ef, 2)
            for ir in 1:size(ef, 1)
                ### r ###
                if ir == 1
                    Δp_r_1 = p[ir+1 ,iφ, iz] - p[ir ,iφ, iz]
                    d_r_1 = axr[ir+1] - axr[ir]
                    er = Δp_r_1 / d_r_1
                elseif ir == size(ef,1)
                    Δp_r_1 = p[ir ,iφ, iz] - p[ir-1 ,iφ, iz]
                    d_r_1 = axr[ir] - axr[ir-1]
                    er = Δp_r_1 / d_r_1
                else
                    Δp_r_1 = p[ir+1 ,iφ, iz]-p[ir ,iφ, iz]
                    Δp_r_2 = p[ir ,iφ, iz]-p[ir-1 ,iφ, iz]
                    d_r_1 = axr[ir+1]-axr[ir]
                    d_r_2 = axr[ir]-axr[ir-1]
                    er = ( Δp_r_1/d_r_1 + Δp_r_2/d_r_2) / 2
                end
                ### φ ###
                if iφ < 2
                    Δp_φ_1 = p[ir ,iφ+1, iz]-p[ir ,iφ, iz]
                    Δp_φ_2 = p[ir ,iφ, iz]-p[ir ,end, iz]
                    d_φ_1 = (axφ[iφ+1]-axφ[iφ]) * axr[ir]# to get the proper value in length units
                    d_φ_2 = (cyclic - axφ[end]) * axr[ir]
                    eφ = ( Δp_φ_1/d_φ_1 + Δp_φ_2/d_φ_2) / 2
                elseif iφ == size(ef,2)
                    Δp_φ_1 = p[ir ,1, iz]-p[ir ,iφ, iz]
                    Δp_φ_2 = p[ir ,iφ, iz]-p[ir ,iφ-1, iz]
                    d_φ_1 = (axφ[1]-axφ[iφ]) * axr[ir]# to get the proper value in length units
                    d_φ_2 = (axφ[iφ]-axφ[iφ-1]) * axr[ir]
                    eφ = ( Δp_φ_1/d_φ_1 + Δp_φ_2/d_φ_2) / 2
                else
                    Δp_φ_1 = p[ir ,iφ+1, iz]-p[ir ,iφ, iz]
                    Δp_φ_2 = p[ir ,iφ, iz]-p[ir ,iφ-1, iz]
                    d_φ_1 = (axφ[iφ+1]-axφ[iφ]) * axr[ir]# to get the proper value in length units
                    d_φ_2 = (axφ[iφ]-axφ[iφ-1]) * axr[ir]
                    eφ = ( Δp_φ_1/d_φ_1 + Δp_φ_2/d_φ_2) / 2
                end
                isinf(eφ) || isnan(eφ) ? eφ = 0.0 : nothing # for small radii and small distances(center of the grid) it would yield Infs or Nans
                if iz-1<1
                    Δp_z_1 = p[ir ,iφ, iz+1]-p[ir ,iφ, iz]
                    d_z_1 = axz[iz+1]-axz[iz]
                    ez = ( Δp_z_1/d_z_1 )
                elseif  iz+1 > size(ef,3)
                    Δp_z_1 = p[ir ,iφ, iz]-p[ir ,iφ, iz-1]
                    d_z_1 = axz[iz]-axz[iz-1]
                    ez = ( Δp_z_1/d_z_1 )
                else
                    Δp_z_1 = p[ir ,iφ, iz+1]-p[ir ,iφ, iz]
                    Δp_z_2 = p[ir ,iφ, iz]-p[ir ,iφ, iz-1]
                    d_z_1 = axz[iz+1]-axz[iz]
                    d_z_2 = axz[iz]-axz[iz-1]
                    ez = ( Δp_z_1/d_z_1 + Δp_z_2/d_z_2) / 2
                end
                if pointtypes[ir, iφ, iz] & update_bit == 0 # boundary points
                    if (1 < ir < size(pointtypes, 1))
                        if (pointtypes[ir - 1, iφ, iz] & update_bit > 0) && (pointtypes[ir + 1, iφ, iz] & update_bit > 0)
                            er = 0
                        elseif (pointtypes[ir - 1, iφ, iz] & update_bit > 0) || (pointtypes[ir + 1, iφ, iz] & update_bit > 0)
                            er *= 2
                        end
                    end
                    if (1 < iφ < size(pointtypes, 2))
                        if (pointtypes[ir, iφ - 1, iz] & update_bit > 0) && (pointtypes[ir, iφ + 1, iz] & update_bit > 0)
                            eφ = 0
                        elseif (pointtypes[ir, iφ - 1, iz] & update_bit > 0) || (pointtypes[ir, iφ + 1, iz] & update_bit > 0)
                            eφ *= 2
                        end
                    end
                    if (1 < iz < size(pointtypes, 3))
                        if (pointtypes[ir, iφ, iz - 1] & update_bit > 0) && (pointtypes[ir, iφ, iz + 1] & update_bit > 0)
                            ez = 0
                        elseif (pointtypes[ir, iφ, iz - 1] & update_bit > 0) || (pointtypes[ir, iφ, iz + 1] & update_bit > 0)
                            ez *= 2
                        end
                    end

                end
                ef[ir,iφ,iz] = [-er, -eφ, -ez]
            end
        end
    end
    if fieldvector_coordinates == :xyz
        ef = convert_field_vectors_to_xyz(ef, axφ)
    end
    return ElectricField(ef, pointtypes.grid)
end

function get_component_field(ef,component=:r,cutoff=NaN)
    components = [:r,:phi,:z]
    # component_index = findfirst(components,component)
    component_index = findfirst(x->x==component,components)
    ef_component = Array{Float32}(undef,size(ef,1),size(ef,2),size(ef,3))
    for iz in 1:size(ef, 3)
        for iφ in 1:size(ef, 2)
            for ir in 1:size(ef, 1)
                if !isnan(cutoff)
                    if abs(ef[ir,iφ ,iz][component_index]) >=cutoff
                        ef_component[ir,iφ,iz] = 0.0
                    else
                        ef_component[ir,iφ,iz] = ef[ir,iφ ,iz][component_index]
                    end
                else
                    ef_component[ir,iφ,iz] = ef[ir,iφ ,iz][component_index]
                end
            end
        end
    end
    return ef_component
end
function get_xyz_vector_from_rφz_vector(v::AbstractArray)::AbstractArray
    return [v[1]*cos(v[2]),v[1]*sin(v[2]),v[3]]
end

function convert_field_vectors_to_xyz(field::Array{SArray{Tuple{3},T,1,3},3}, φa::Array{T, 1})::Array{SVector{3, T},3} where {T}
    field_xyz = Array{SVector{3,T},3}(undef, size(field)...);
    for (iφ, φ) in enumerate(φa)
        Rα::SMatrix{3,3,T} = @SArray([cos(φ) -1*sin(φ) 0;sin(φ) cos(φ) 0;0 0 1])
        for iz in axes(field, 3)
            for ir in axes(field, 1)
                field_xyz[ir,iφ,iz] = Rα * field[ir,iφ,iz]
                # field_xyz[ir,iφ,iz] = get_xyz_vector_from_field_vector(field[ir,iφ,iz], φ)
            end
        end
    end
    return field_xyz
end
#
# function get_xyz_vector_from_field_vector(vector, α)
#         Rα = @SArray([cos(α) -1*sin(α) 0;sin(α) cos(α) 0;0 0 1])
#         result = Rα*vector
#         result
# end
function interpolated_scalarfield(ep::ScalarPotential{T, 3, :cylindrical}) where {T}
    knots = ep.grid.axes[1].ticks, cat(ep.grid.axes[2].ticks,T(2π),dims=1), ep.grid.axes[3].ticks#(grid.r, grid.φ, grid.z)
    ext_data = cat(ep.data, ep.data[:,1:1,:], dims=2)
    i = interpolate(knots, ext_data, Gridded(Linear()))
    vector_field_itp = extrapolate(i, (Interpolations.Line(), Periodic(), Interpolations.Line()))
    return vector_field_itp
end
function interpolated_scalarfield(ep::ScalarPotential{T, 3, :cartesian}) where {T}
    knots = ep.grid.axes#(grid.x, grid.y, grid.z)
    i = interpolate(knots, ep.data, Gridded(Linear()))
    vector_field_itp = extrapolate(i, (Interpolations.Line(), Interpolations.Line(), Interpolations.Line()))
    return vector_field_itp
end
#
# function interpolated_vectorfield(vectorfield::AbstractArray{<:SVector{3, T},3}, ep::ElectricPotential{T}) where {T}
#     knots = ep.grid.axes#(grid.r, grid.φ, grid.z)
#     i = interpolate(knots, vectorfield, Gridded(Linear()))
#     vectorfield_itp = extrapolate(i, Periodic())
#     return vectorfield_itp
# end
#
# function setup_interpolated_vectorfield(vectorfield, grid::CylindricalGrid{T}) where {T}
#     knots = grid.axes #(grid.r, grid.φ, grid.z)
#     i = interpolate(knots, vectorfield, Gridded(Linear()))
#     vector_field_itp = extrapolate(i, Periodic())
#     return vector_field_itp
# end

function get_interpolated_drift_field(velocityfield, grid::CylindricalGrid{T}) where {T}
    extended_velocityfield = cat(velocityfield, velocityfield[:,1:1,:], dims=2)
    knots = grid.axes[1].ticks, cat(grid.axes[2].ticks,T(2π),dims=1), grid.axes[3].ticks
    i = interpolate(knots, extended_velocityfield, Gridded(Linear()))
    velocity_field_itp = extrapolate(i, (Interpolations.Line(), Periodic(), Interpolations.Line()))
    return velocity_field_itp
end
function get_interpolated_drift_field(velocityfield, grid::CartesianGrid{T}) where {T}
    knots = grid.axes[:1].ticks, grid.axes[:2].ticks, grid.axes[:3].ticks
    i = interpolate(knots, velocityfield, Gridded(Linear()))
    velocity_field_itp = extrapolate(i, (Interpolations.Line(), Interpolations.Line(), Interpolations.Line()))
    return velocity_field_itp
end

include("plot_recipes.jl")

function get_electric_field_from_potential(ep::ElectricPotential{T, 3, :cartesian}, pointtypes::PointTypes{T})::ElectricField{T, 3, :cartesian} where {T <: SSDFloat}
    axx::Vector{T} = collect(ep.grid[:x])
    axy::Vector{T} = collect(ep.grid[:y])
    axz::Vector{T} = collect(ep.grid.z)
    axx_ext::Vector{T} = get_extended_ticks(ep.grid[:x])
    axy_ext::Vector{T} = get_extended_ticks(ep.grid[:y])
    axz_ext::Vector{T} = get_extended_ticks(ep.grid.z)

    ef::Array{SVector{3, T}} = Array{SVector{3, T}}(undef, size(ep.data))

    for ix in eachindex(axx)
        for iy in eachindex(axy)
            for iz in eachindex(axz)
                if ix - 1 < 1
                    Δp_x_1::T = ep.data[ix + 1, iy, iz] - ep.data[ix, iy, iz]
                    d_x_1::T = axx[ix + 1] - axx[ix]
                    ex::T =  Δp_x_1 / d_x_1
                elseif ix + 1 > size(ef, 1)
                    Δp_x_1 = ep.data[ix, iy, iz] - ep.data[ix - 1, iy, iz]
                    d_x_1 = axx[ix] - axx[ix - 1]
                    ex = Δp_x_1 / d_x_1
                else
                    Δp_x_1 = ep.data[ix + 1, iy, iz] - ep.data[ix ,iy, iz]
                    Δp_x_2::T = ep.data[ix, iy, iz] - ep.data[ix - 1, iy, iz]
                    d_x_1 = axx[ix + 1] - axx[ix]
                    d_x_2::T = axx[ix] - axx[ix - 1]
                    ex = (Δp_x_1 / d_x_1 + Δp_x_2 / d_x_2) / 2
                end

                if iy - 1 < 1
                    Δp_y_1::T = ep.data[ix, iy + 1, iz] - ep.data[ix ,iy, iz]
                    d_y_1::T = axy[iy + 1] - axy[iy]
                    ey::T =  Δp_y_1 / d_y_1
                elseif iy + 1 > size(ef, 2)
                    Δp_y_1 = ep.data[ix, iy, iz] - ep.data[ix, iy - 1, iz]
                    d_y_1 = axy[iy] - axy[iy - 1]
                    ey = Δp_y_1 / d_y_1
                else
                    Δp_y_1 = ep.data[ix, iy + 1, iz] - ep.data[ix ,iy, iz]
                    Δp_y_2::T = ep.data[ix, iy, iz] - ep.data[ix, iy - 1, iz]
                    d_y_1 = axy[iy + 1] - axy[iy]
                    d_y_2::T = axy[iy] - axy[iy - 1]
                    ey = (Δp_y_1 / d_y_1 + Δp_y_2 / d_y_2) / 2
                end

                if iz - 1 < 1
                    Δp_z_1::T = ep.data[ix, iy, iz + 1] - ep.data[ix, iy, iz]
                    d_z_1::T = axz[iz + 1] - axz[iz]
                    ez::T =  Δp_z_1 / d_z_1
                elseif iz + 1 > size(ef, 3)
                    Δp_z_1 = ep.data[ix, iy, iz] - ep.data[ix, iy, iz - 1]
                    d_z_1 = axz[iz] - axz[iz - 1]
                    ez = Δp_z_1 / d_z_1
                else
                    Δp_z_1 = ep.data[ix, iy, iz + 1] - ep.data[ix ,iy, iz]
                    Δp_z_2::T = ep.data[ix, iy, iz] - ep.data[ix, iy, iz - 1]
                    d_z_1 = axz[iz + 1] - axz[iz]
                    d_z_2::T = axz[iz] - axz[iz - 1]
                    ez = (Δp_z_1 / d_z_1 + Δp_z_2 / d_z_2) / 2
                end

                if pointtypes[ix, iy, iz] & update_bit == 0 # boundary points
                    if (1 < ix < size(pointtypes, 1))
                        if (pointtypes[ix - 1, iy, iz] & update_bit > 0) && (pointtypes[ix + 1, iy, iz] & update_bit > 0)
                            ex = 0
                        elseif (pointtypes[ix - 1, iy, iz] & update_bit > 0) || (pointtypes[ix + 1, iy, iz] & update_bit > 0)
                            ex *= 2
                        end
                    end
                    if (1 < iy < size(pointtypes, 2))
                        if (pointtypes[ix, iy - 1, iz] & update_bit > 0) && (pointtypes[ix, iy + 1, iz] & update_bit > 0)
                            ey = 0
                        elseif (pointtypes[ix, iy - 1, iz] & update_bit > 0) || (pointtypes[ix, iy + 1, iz] & update_bit > 0)
                            ey *= 2
                        end
                    end
                    if (1 < iz < size(pointtypes, 3))
                        if (pointtypes[ix, iy, iz - 1] & update_bit > 0) && (pointtypes[ix, iy, iz + 1] & update_bit > 0)
                            ez = 0
                        elseif (pointtypes[ix, iy, iz - 1] & update_bit > 0) || (pointtypes[ix, iy, iz + 1] & update_bit > 0)
                            ez *= 2
                        end
                    end
                end
                ef[ix, iy, iz] = @SVector [-ex, -ey, -ez]
            end
        end
    end
    return ElectricField(ef, pointtypes.grid)
end

function get_electric_field_strength(ef::ElectricField{T}) where {T <: SSDFloat}
    efs::Array{T, 3} = Array{T, 3}(undef, size(ef.data))
    @inbounds for i in eachindex(ef.data)
        efs[i] = norm(ef.data[i])
    end
    return efs
end


text:
import CSV

candidates = CSV.read("data/candidates.csv")
total = CSV.read("data/totalvotes.csv")

ext: jl

max_stars_repo_name: TnTo/Electoral-Flows

max_issues_repo_name: TnTo/Electoral-Flows

max_forks_repo_name: TnTo/Electoral-Flows


text:
#=
image_geom.jl
Methods related to an image geometry for image reconstruction
2017-10-02 Samuel Rohrer, University of Michigan
2019-03-05 Jeff Fessler, Julia 1.1 + tests
2019-06-23 Jeff Fessler, overhaul
2020-07-21 Jeff Fessler, D-dimensional
=#

export ImageGeom, image_geom, cbct
export image_geom_circle
export image_geom_ellipse

#using MIRT: downsample2, downsample3
using FillArrays: Trues
using ImageTransformations: imresize

"""
    ImageGeom{D}

- `dims::Dims{D}` image dimensions
- `deltas::NTuple{D,Float32}` pixel sizes
- `offsets::NTuple{D,Float32}` unitless
- `mask::AbstractArray{Bool,D}` logical mask

Image geometry "struct" with essential parameters
"""
struct ImageGeom{D}
	dims::Dims{D} # image dimensions
	deltas::NTuple{D,Float32} # pixel sizes
	offsets::NTuple{D,Float32} # unitless
	mask::AbstractArray{Bool,D} # logical mask

	function ImageGeom{D}(
		dims::Dims{D},
		deltas::NTuple{D,Real},
		offsets::NTuple{D,Real},
		mask::AbstractArray{Bool,D},
	) where {D}
		any(dims .<= 0) && throw("dims must be positive")
		any(deltas .== 0) && throw("deltas must be nonzero")
		size(mask) == dims ||
			throw(DimensionMismatch("mask size $(size(mask)) vs dims $dims"))
		new{D}(dims, Float32.(deltas), Float32.(offsets), mask)
	end
end

"""
    ImageGeom{D}(dims, deltas, offsets, [, mask])
Constructor for `ImageGeom`
Default `mask` is `Trues(dims)` which is akin to `trues(dims)`.
"""
ImageGeom{D}(
	dims::Dims{D},
	deltas::NTuple{D,Real},
	offsets::NTuple{D,Real},
) where {D} = ImageGeom{D}(dims, deltas, offsets, Trues(dims))

"""
image_geom_help( ; io)
"""
function image_geom_help( ; io::IO = isinteractive() ? stdout : devnull)
	print(io, "propertynames:\n")
	print(io, propertynames(image_geom(nx=1, dx=1)))

	print(io,
	"
	Derived values for 2D (and 3D) images:

	is3	is it 3d (nz > 0?)
	fovs (|dx|*nx, |dy|*ny, ...)
	np	sum(mask) = # pixels to be estimated

	dim	(nx, ny, [nz])
	x	1D x coordinates of each pixel
	y	1D y coordiantes of each pixel
	wx	(nx - 1/2) * dx + offset_x
	wy	(ny - 1/2) * dy + offset_y
	wz	(nz - 1/2) * dz + offset_z

	xg	x coordinates of each pixel as a grid (2D or 3D)
	yg	y coordinates of each pixel as a grid (2D or 3D)

	u	1D frequency domain coordinates of each pixel (kx)
	v	1D frequency domain coordinates of each pixel (ky)
	w	1D frequency domain coordinates of each pixel (kz)

	ug	2D or 3D grid of frequency domain coordinates
	vg	2D or 3D grid of frequency domain coordinates
	fg	NamedTuple of 2D or 3D frequency coordinates
	ones	ones(dim)
	zeros	zeros(dim)

	mask_outline	binary image showing 2D outline of mask

	Derived values for 3D images:

	z	z coordinates of each pixel
	zg	z coordinates of each pixel as a grid (3D only)
	wg	3D grid of frequency domain coordinates

	mask_or	[nx ny] logical 'or' of mask

	Methods:

	embed(x)	turn short column(s) into array(s)
	maskit(x)	opposite of embed
	shape(x)	reshape long column x to [nx ny [nz]] array
	unitv(...)	j=j | i=(ix,iy) | c=(cx cy)
				j: single index from 1 to length(z)
				i: (ix,iy[,iz]) index from 1 to nx,ny
				c: (cx,cy[,cz]) index from +/- n/2 center at floor(n/2)+1
	circ(rx=,ry=,cx=,cy=)	circle of given radius and center (cylinder in 3D)
	plot(jim)	plot the image geometry using the `jim` function

	Methods that return a new `ImageGeom:`

	down(down::Int)		down-sample geometry by given factor
	over(over::Int)		over-sample geometry by given factor
	expand_nz(nz_pad)	expand image geometry in z by nz_pad on both ends
	\n")
end

"""
    MIRT_cbct_ig

Structure suitable for passing to C routines `cbct_*`
based on the struct `cbct_ig` found in `cbct,def.h`
"""
struct MIRT_cbct_ig
	nx::Cint
	ny::Cint
	nz::Cint
	dx::Cfloat
	dy::Cfloat
	dz::Cfloat
	offset_x::Cfloat
	offset_y::Cfloat
	offset_z::Cfloat
	mask2::Ptr{Cuchar}	# [nx,ny] 2D support mask: 0 or 1 .. nthread
	iy_start::Ptr{Cint} # [nthread] for mask2
	iy_end::Ptr{Cint}	# [nthread] for mask2
end

"""
cbct(ig::ImageGeom; nthread=1)
constructor for `MIRT_cbct_ig`
"""
function cbct(ig::ImageGeom{3} ; nthread::Int=1)
	iy_start = [0]
	iy_end = [ig.ny]
	nthread != 1 && throw("only nthread=1 for now due to iy_start")
	return MIRT_cbct_ig(
		Cint(ig.nx), Cint(ig.ny), Cint(ig.nz), Cfloat(ig.dx),
		Cfloat(ig.dy), Cfloat(ig.dz),
		Cfloat(ig.offset_x), Cfloat(ig.offset_y), Cfloat(ig.offset_z),
		pointer(UInt8.(ig.mask_or)),
		pointer(Cint.(iy_start)), pointer(Cint.(iy_end)),
	)
end

# see https://docs.julialang.org/en/stable/manual/constructors/
"""
    ig = image_geom(...)

Constructor for `ImageGeom`

# Arguments
- `nx::Int = 128`
- `ny::Int = nx`
- `dx::Real = ?` (must specify one of `dx` or `fov`)
- `dy::Real = -dx`
- `offset_x::Real = 0` (unitless)
- `offset_y::Real = 0` (unitless)
- `fov::Real = ?` (if specified, then `nx*dx=ny*dy`)
- `nz::Int = 0`
- `dz::Real = ?` (need one of `dz` or `zfov` if `nz > 0`)
- `zfov::Real = ?` (if specified, then `nz*dz`)
- `offset_z::Real = 0` (unitless)
- `offsets::Symbol = :none` or :dsp
- `mask::Union{Symbol,AbstractArray{Bool}} = :all` | `:circ` | `:all_but_edge_xy`
"""
function image_geom( ;
	nx::Int = 128,
	ny::Int = nx,
	dx::Real = NaN,
	dy::Real = NaN,
	offset_x::Real = 0,
	offset_y::Real = 0,
	fov::Real = NaN,
	nz::Int = 0,
	dz::Real = NaN,
	zfov::Real = NaN,
	offset_z::Real = 0,
	offsets::Symbol = :none,
	mask::Union{Symbol,AbstractArray{Bool}}	= :all,
)

	# handle optional arguments

	is3 = !isnan(dz) || !isnan(zfov) || (nz > 0)

	# offsets
	if offsets === :dsp
		# check that no nonzero offset was specified
		(offset_x != 0 || offset_y != 0) && throw("offsets usage incorrect")
		offset_x = 0.5
		offset_y = 0.5
		if is3
			offset_z != 0 && throw("offset_z usage incorrect")
			offset_z = 0.5
		end
	elseif offsets != :none
		throw("offsets $offsets")
	end

	# transverse plane (2d) distances
	if true
		if isnan(fov)
			isnan(dx) && throw("dx or fov required")
		#	fov = nx * dx
		else
			!isnan(dx) && throw("dx and fov?")
			dx = fov / nx
			dy = fov / ny
		end

		dy = isnan(dy) ? -dx : dy # default dy
	end

	# 3D geometry
	if is3
		if isnan(zfov)
			isnan(dz) && throw("dz or zfov required")
		#	zfov = nz * dz
		else
			!isnan(dz) && throw("dz and zfov?")
			dz = zfov / nz
		end
	end

	# mask
	if mask === :all
		mask = is3 ? Trues(nx, ny, nz) : Trues(nx, ny)
	elseif mask === :circ
		mask = image_geom_circle(nx, ny, dx, dy)
		if nz > 0
			mask = repeat(mask, 1, 1, nz)
		end
	elseif mask === :all_but_edge_xy
		mask = trues(nx,ny,max(nz,1))
		mask[  1,   :, :] .= false
		mask[end,   :, :] .= false
		mask[  :,   1, :] .= false
		mask[  :, end, :] .= false
		if !is3
			mask = mask[:,:,1]
		end
	elseif isa(mask, Symbol)
		throw("mask symbol $mask")
	elseif size(mask,1) != nx || size(mask,2) != ny || (is3 && size(mask,3) != nz)
		throw("mask size $(size(mask)), nx=$nx ny=$ny nz=$nz")
	end

	# return the object
	return is3 ?
		ImageGeom{3}((nx, ny, nz), (dx, dy, dz),
			(offset_x, offset_y, offset_z), mask) :
		ImageGeom{2}((nx, ny), (dx, dy),
			(offset_x, offset_y), mask)
end

"""
(n,d) = image_geom_down_round()
helper function needed to downsample `image_geom`
"""
function _down_round(
	val::NTuple{D,Real},
	dd::NTuple{D,Real},
	down::NTuple{D,Real},
) where {D}
	out = val ./ down
	# for non-divisors make dim a multiple of 2
	fun = out -> out == round(out) ? out : 2 * round(out / 2)
	out = fun.(out)
	dd = dd .* down
	return Int.(out), dd
end

"""
    ig_down = ig_downsample(ig, down::Tuple{Int})
down sample an image geometry by the factor `down`
cf `image_geom_downsample`
"""
function ig_downsample(ig::ImageGeom{D}, down::NTuple{D,Int}) where {D}

	# call the down round function
	down_dim, deltas = _down_round(ig.dims, ig.deltas, down)
	# adjust to "pixel" units
	down_offsets = ig.offsets ./ down

	# carefully down-sample the mask
	mdim = size(ig.mask)
	if ig.is3
		if down_dim .* down == mdim
			down_mask = downsample3(ig.mask, down) .> 0
		else
			throw("bug: bad mask size. need to address mask downsampling")
		end
	else
		if down_dim .* down == mdim
			down_mask = downsample2(ig.mask, down) .> 0
		else
			down_mask = imresize(ig.mask, down_dim...) .> 0
		end
	end

	return ImageGeom{D}(down_dim, deltas, down_offsets, down_mask)
end

ig_downsample(ig::ImageGeom{D}, down::Int) where {D} =
	ig_downsample(ig, ntuple(i->down, D))

"""
ig_new = image_geom_expand_nz(ig::ImageGeom{3}, nz_pad::Int)
pad both ends
"""
function image_geom_expand_nz(ig::ImageGeom{3}, nz_pad::Int)
	out_nz = ig.nz + 2*nz_pad
	out_mask = cat(dims=3, repeat(ig.mask[:,:,1], 1, 1, nz_pad),
		ig.mask, repeat(ig.mask[:,:,end], 1, 1, nz_pad),
	)
	return ImageGeom{3}((ig.dims[1], ig.dims[2], out_nz),
		ig.deltas, ig.offsets, out_mask,
	)
end

"""
ig_over = image_geom_over(ig::ImageGeom, over::Int)
over-sample an image geometry by the factor `over`
"""
function image_geom_over(ig::ImageGeom{D}, over::Int) where {D}
	if ig.mask == Trues(ig.dims) || all(ig.mask .== true)
		mask_over = Trues(ig.dims .* over)
	else
		mask_over = imresize(ig.mask, ig.dims .* over) .> 0
	end
	return ImageGeom{D}(ig.dims .* over, ig.deltas ./ over,
		ig.offsets .* over, mask_over,
	)
end

# ellipse that just inscribes the rectangle
# but keeping a 1 pixel border due to ASPIRE regularization restriction
function image_geom_ellipse(
	nx::Int, ny::Int, dx::Real, dy::Real ;
	rx::Real = min(abs((nx/2-1)*dx), abs((ny/2-1)*dy)),
	ry::Real = min(abs((nx/2-1)*dx), abs((ny/2-1)*dy)),
	cx::Real = 0, cy::Real = 0, over::Int=2,
)
	ig = image_geom(nx=nx, ny=ny, dx=dx, dy=dy)
	circ = ellipse_im(ig, [cx cy rx ry 0 1], oversample=over) .> 0
	return circ
end

# default is a circle that just inscribes the square
# but keeping a 1 pixel border due to ASPIRE regularization restriction
function image_geom_circle(
	nx::Int, ny::Int, dx::Real, dy::Real ;
	rx::Real = min(abs((nx/2-1)*dx), abs((ny/2-1)*dy)),
	ry::Real = rx, cx::Real = 0, cy::Real = 0, nz::Int=0, over::Int=2,
)
	ig = image_geom(nx=nx, ny=ny, dx=dx, dy=dy)
	circ = ellipse_im(ig, [cx cy rx ry 0 1], oversample=over) .> 0
	if nz > 0
		circ = repeat(circ, 1, 1, nz)
	end
	return circ
end

"""
out = image_geom_add_unitv(z::AbstractArray ; j=?, i=?, c=?)

add a unit vector to an initial array `z` (typically of zeros)

# options (use at one of these):
- `j` single index from 1 to length(z)
- `i` (ix, iy [,iz]) index from 1 to nx,ny
- `c` (cx, cy [,cz]) index from +/- n/2 center at floor(n/2)+1

default with no arguments gives unit vector at center `c=(0, 0 [,0])`
"""
function image_geom_add_unitv(
	z::AbstractArray{T, D} ; # starts with zeros()
	j::Int = 0,
	i::NTuple{D,Int} = ntuple(i -> 0, D),
	c::NTuple{D,Int} = ntuple(i -> 0, D),
) where {T <: Number, D}

	out = copy(z)

	(j < 0 || j > length(z)) && throw("bad j $j")
	if 1 <= j <= length(z)
		any(i .!= 0) && throw("i $i")
		any(c .!= 0) && throw("c $c")
		out[j] += one(T)
	elseif all(1 .<= i .<= size(z))
		any(c .!= 0) && throw("c $c")
		out[i...] += one(T)
	else
		any(i .!= 0) && throw("i $i")
		tmp = c .+ Int.(floor.(size(out) ./ 2)) .+ 1
		out[tmp...] += one(T)
	end

	return out
end

"""
image_geom_plot(ig, how ; kwargs...)
The `how` argument should be `jim` to be useful.
"""
image_geom_plot(ig::ImageGeom{2}, how::Function ; kwargs...) =
	how(ig.x, ig.y, ig.mask, "(nx,ny)=$(ig.nx),$(ig.ny)" ; kwargs...)
image_geom_plot(ig::ImageGeom{3}, how::Function ; kwargs...) =
	how(ig.x, ig.y, ig.mask_or,
		"(dx,dy,dz)=$(ig.dx),$(ig.dy),$(ig.dz)" ; kwargs...)

"""
    show(io::IO, ig::ImageGeom)
    show(io::IO, ::MIME"text/plain", ig::ImageGeom)
"""
Base.show(io::IO, ig::ImageGeom) =
	print(io, "ImageGeom: $(ig.dim)")

function Base.show(io::IO, ::MIME"text/plain", ig::ImageGeom)
	ir_dump(io, ig)
end

# Extended properties

image_geom_fun0 = Dict([
	(:help, ig -> image_geom_help()),

	(:ndim, ig -> length(ig.dims)),
	(:nx, ig -> ig.dims[1]),
	(:ny, ig -> ig.ndim >= 2 ? ig.dims[2] : 0),
	(:nz, ig -> ig.ndim >= 3 ? ig.dims[3] : 0),
	(:is3, ig -> ig.nz > 0),
	(:dim, ig -> ig.dims),
	(:fovs, ig -> abs.(ig.deltas) .* ig.dims),

	(:zeros, ig -> zeros(Float32, ig.dim)),
	(:ones, ig -> ones(Float32, ig.dim)),

	(:dx, ig -> ig.deltas[1]),
	(:dy, ig -> ig.ndim >= 2 ? ig.deltas[2] : Float32(0)),
	(:dz, ig -> ig.ndim >= 3 ? ig.deltas[3] : Float32(0)),

	(:offset_x, ig -> ig.offsets[1]),
	(:offset_y, ig -> ig.ndim >= 2 ? ig.offsets[2] : Float32(0)),
	(:offset_z, ig -> ig.ndim >= 3 ? ig.offsets[3] : Float32(0)),

	(:x, ig -> ((0:ig.nx-1) .- ig.wx) * ig.dx),
	(:y, ig -> ((0:ig.ny-1) .- ig.wy) * ig.dy),
	(:z, ig -> ((0:ig.nz-1) .- ig.wz) * ig.dz),

	(:wx, ig -> (ig.nx - 1)/2 + ig.offset_x),
	(:wy, ig -> (ig.ny - 1)/2 + ig.offset_y),
	(:wz, ig -> (ig.nz - 1)/2 + ig.offset_z),

	(:xg, ig -> ig.is3 ? repeat(ig.x, 1, ig.ny, ig.nz) :
						repeat(ig.x, 1, ig.ny)),
	(:yg, ig -> ig.is3 ? repeat(ig.y', ig.nx, 1, ig.nz) :
						repeat(ig.y', ig.nx, 1)),
	(:zg, ig -> ig.is3 ? repeat(reshape(ig.z, 1, 1, ig.nz), ig.nx, ig.ny, 1) :
						zeros(ig.nx, ig.ny)),

	# DFT frequency sample grid
	(:u, ig -> (-ig.nx/2:(ig.nx/2-1)) / (ig.nx*ig.dx)),
	(:v, ig -> (-ig.ny/2:(ig.ny/2-1)) / (ig.ny*ig.dy)),
	(:w, ig -> (-ig.nz/2:(ig.nz/2-1)) / (ig.nz*ig.dz)),

	(:ug, ig -> ig.is3 ? repeat(ig.u, 1, ig.ny, ig.nz) :
						repeat(ig.u, 1, ig.ny)),
	(:vg, ig -> ig.is3 ? repeat(ig.v', ig.nx, 1, ig.nz) :
						repeat(ig.v', ig.nx, 1)),
	(:wg, ig -> ig.is3 ? repeat(reshape(ig.w, 1, 1, ig.nz), ig.nx, ig.ny, 1) :
						zeros(ig.nx, ig.ny)),
	(:fg, ig -> ig.is3 ? (ig.ug, ig.vg, ig.wg) : (ig.ug, ig.vg)),

	(:np, ig -> sum(ig.mask)),
	(:mask_or, ig -> mask_or(ig.mask)),
	(:mask_outline, ig -> mask_outline(ig.mask)),

	# simple functions

	(:plot, ig -> ((how::Function ; kwargs...) ->
		image_geom_plot(ig, how ; kwargs...))),
	(:embed, ig -> (x::AbstractArray{<:Number} -> embed(x, ig.mask))),
	(:maskit, ig -> (x::AbstractArray{<:Number} -> maskit(x, ig.mask))),
	(:shape, ig -> (x::AbstractArray{<:Number} -> reshape(x, ig.dim))),
	(:unitv, ig -> (( ; kwargs...) -> image_geom_add_unitv(ig.zeros ; kwargs...))),
	(:circ, ig -> (( ; kwargs...) ->
		image_geom_circle(ig.nx, ig.ny, ig.dx, ig.dy, nz=ig.nz ; kwargs...))),

	# functions that return new geometry:

	(:down, ig -> (down::Int -> ig_downsample(ig, down))),
	(:over, ig -> (over::Int -> image_geom_over(ig, over))),
	(:expand_nz, ig -> (nz_pad::Int -> image_geom_expand_nz(ig, nz_pad))),

	])

# Tricky overloading here!

Base.getproperty(ig::ImageGeom, s::Symbol) =
		haskey(image_geom_fun0, s) ? image_geom_fun0[s](ig) :
		getfield(ig, s)

Base.propertynames(ig::ImageGeom) =
	(fieldnames(typeof(ig))..., keys(image_geom_fun0)...)


text:
module ViewHelper

using Genie, Genie.Helpers, SearchLight, Genie.Router
using Dates
using Tasks

export output_flash, task_form, task_tr

function output_flash(params::Dict{Symbol,Any}) :: String
  ! isempty( flash(params) ) ? """
$(flash(params))
""" : ""
end

function task_form(action; task::Tasks.Task = Tasks.Task(content = "", deadline = Date(now())))
    id_field = isnull(task.id) ? "" : """<input type="hidden" name="task_id" value="$(get(task.id))"/>"""
    """
    <form action="$(action)" method="POST" id="task-form">

        $(id_field)
        <input id="task-deadline" class="form-control" type="date" name="task_deadline" value="$(task.deadline)"/>
        <input id="task-content" class="form-control" type="text" name="task_content" placeholder="Write down your TODO" value="$(task.content)" style="width: 50%" />

          <input type="submit" value="$(isnull(task.id) ? "Create" : "Update")" class="btn btn-sm btn-primary"/>

    </form>
    """
end

function task_tr(task)
    today = Date(now())
    """

            <input type="checkbox" id="$(get(task.id))" data-id="$(get(task.id))"
                   class="checkbox done-checkbox" $(task.done ? "checked" : "") />

            <span id="deadline_$(get(task.id))" class="$(task.done ? "done" : task.deadline <= today ? "deadline" : "doing")">
                  $(task.deadline)

            <span id="content_$(get(task.id))" class="$(task.done ? "done" : task.deadline <= today ? "deadline" : "doing")">
                  $(task.content)

              Edit
              <button data-id="$(get(task.id))" class="btn btn-sm btn-danger delete-button">Delete</button>

    """
end

end

max_stars_repo_name: ky-s/MyGenieApp

max_issues_repo_name: ky-s/MyGenieApp

max_forks_repo_name: ky-s/MyGenieApp


text:
mutable struct Scenario{T}
    p::Int
    model::AbstractGameModel
    roadway::Roadway{T}
    player::Vector{Player{T}}
end

function Scenario(model::AbstractGameModel, roadway::Roadway{T}, players::Vector{<:Player}) where {T}
    p = model.p
    @assert p == length(players)
    sort!(players, by = x -> x.id)
    if [players[i].id for i=1:p] != (1:p)
        # @show "Reindexing the players"
        for i = 1:p
            players[i].id = i
        end
    end
    sce = Scenario{T}(p,model,roadway,players)
    return sce
end

################################################################################
# Getters
################################################################################

function get_state(sce::Scenario{T}) where {T}
    model = sce.model
    n = model.n
    p = model.p
    pz = model.pz
    x0 = zeros(n)
    for i = 1:p
        x0[pz[i]] = sce.player[i].x0
    end
    return SVector{n,T}(x0)
end

function get_control_bounds(sce::Scenario{T}) where {T}
    model = sce.model
    m = model.m
    p = model.p
    pu = model.pu
    u_min = zeros(m)
    u_max = zeros(m)
    for i = 1:p
        u_min[pu[i]] = sce.player[i].u_min
        u_max[pu[i]] = sce.player[i].u_max
    end
    return SVector{m,T}(u_min), SVector{m,T}(u_max)
end

function Algames.GameObjective(N::Int, sce::Scenario{T}) where {T}
    model = sce.model
    n = model.n
    p = model.p

    Q = [sce.player[i].Q for i=1:p]
    R = [sce.player[i].R for i=1:p]
    xf = [sce.player[i].xf for i=1:p]
    uf = [sce.player[i].uf for i=1:p]
    game_obj = Algames.GameObjective(Q, R, xf, uf, N, model)

    radius = [sce.player[i].r_cost for i=1:p]
    μ = [sce.player[i].μ for i=1:p]
    add_collision_cost!(game_obj, radius, μ)
    return game_obj
end

function Algames.GameConstraintValues(N::Int, sce::Scenario{T}) where {T}
    model = sce.model
    roadway = sce.roadway
    p = model.p
    probsize = ProblemSize(N,model)
    game_con = Algames.GameConstraintValues(probsize)

    radius = [sce.player[i].r_col for i=1:p]
    add_collision_avoidance!(game_con, radius)

    u_min, u_max = get_control_bounds(sce)
    if u_min != -Inf*ones(model.m) || u_max != Inf*ones(model.m) # Avoid adding void constraints
        add_control_bound!(game_con, u_max, u_min)
    end

    for i = 1:p
        player = sce.player[i]
        v_min = player.v_min
        v_max = player.v_max
        if v_min != -Inf || v_max != Inf # Avoid adding void constraints
            add_velocity_bound!(model, game_con, i, v_max, v_min)
        end
    end

    for i = 1:p
        player = sce.player[i]
        lane = roadway.lane[player.lane_id]
        if length(lane.circle) > 0 # Avoid adding void constraints
            ri = radius[i]
            xc = [c.x for c in lane.circle]
            yc = [c.y for c in lane.circle]
            rc = [c.r for c in lane.circle]
            add_circle_constraint!(game_con, i, xc, yc, rc .+ ri)
        end
    end

    for i = 1:p
        player = sce.player[i]
        lane = roadway.lane[player.lane_id]
        walls = deepcopy(lane.wall)
        ri = radius[i]
        for j = 1:length(lane.wall)
            walls[j].p1 -= ri*walls[j].v
            walls[j].p2 -= ri*walls[j].v
        end
        add_wall_constraint!(game_con, i, walls)
    end

    return game_con
end

function Algames.GameProblem(N::Int, dt::T, sce::Scenario{T}, opts::Options{T}) where {T}
    x0 = get_state(sce)
    game_obj = GameObjective(N, sce)
    game_con = GameConstraintValues(N, sce)
    prob = Algames.GameProblem(N, dt, x0, sce.model, opts, game_obj, game_con)
    return prob
end


text:
rm("Manifest.toml";force=true)

using Pkg

if ARGS[1] == "full"
    pkgs = ["MadNLPHSL","MadNLPPardiso","MadNLPMumps","MadNLPGPU","MadNLPGraph","MadNLPKrylov"]
elseif ARGS[1] == "basic"
    pkgs = ["MadNLPMumps","MadNLPGraph","MadNLPKrylov"]
else
    error("proper argument should be given - full or basic")
end

Pkg.develop(PackageSpec(path=joinpath(@__DIR__,"..")))
Pkg.develop(PackageSpec(path=joinpath(@__DIR__,"..","lib","MadNLPTests")))
Pkg.develop.([PackageSpec(path=joinpath(@__DIR__,"..","lib",pkg)) for pkg in pkgs])
Pkg.build()

Pkg.test.("MadNLP", coverage=true)
Pkg.test.(pkgs, coverage=true)


text:
using GridWorld

rows = [
    [7:22],[4:22],[2:22],
    [2:22],[3:22],[4:22],
    [5:20],[5:17],[5:16],
    [5:14],[5:13],[4:13],
    [3:13],[1:13],[1:13],
]

function getgridmat(rows)
    m = falses(15, 22)
    for (row, ranges) in enumerate(rows)
        for range in ranges
            m[row, range] .= true
        end
    end
    m
end

mat = getgridmat(rows)

function cellfn(row,col)
    if (row,col) == (7,7)
        return Label("●"; color="red", size=22)
    end
    if !mat[row,col]
        return nothing
    end
    bg_color = nothing
    if row == size(mat, 1)
        bg_color = "#98b8eb"
    elseif col == size(mat, 2)
        bg_color = "#92BD8E"
    end
    Cell(
        nothing,
        line_thickness=1,
        bg_color=bg_color
    )
end

g4fn() = Grid(
    size(mat, 1),
    size(mat, 2),
    cellfn,
    line_thickness=0,
    title_right=Label("Finish Line", color="#333", size=22),
    header=Label("Starting Line", color="#333", size=22),
    margin_top=30,
    margin_bottom=30,
    margin_left=20,
    margin_right=20,
    line_color="red",
    cell_size=20,
)
g4 = g4fn()


text:
import GreyDecision.Utility: makeminmax
import GreyDecision.GreyNumbers: GreyNumber
import GreyDecision.Topsis: topsis

@testset "Topsis sample" begin
    atol = 0.001
    decmat = [
        GreyNumber(1.0, 2.0) GreyNumber(2.0, 3.0) GreyNumber(3.0, 4.0);
        GreyNumber(2.0, 3.0) GreyNumber(1.0, 2.0) GreyNumber(3.0, 4.0);
        GreyNumber(3.0, 4.0) GreyNumber(2.0, 3.0) GreyNumber(1.0, 2.0);
    ]
    w = [0.5, 0.4, 0.1]
    fns = makeminmax([maximum, maximum, minimum])

    scores = topsis(decmat, w, fns)
    @test isapprox(scores[1], GreyNumber{Float64}(0.23506, 0.83613), atol = atol)
    @test isapprox(scores[2], GreyNumber{Float64}(0.24317, 1.00294), atol = atol)
    @test isapprox(scores[3], GreyNumber{Float64}(0.10851, 1.23913), atol = atol)
end

@testset "Grey Topsis with real numbers" begin
    atol = 0.001
    num6 = GreyNumber{Float64}(6.0, 6.0)
    num7 = GreyNumber{Float64}(7.0, 7.0)
    num8 = GreyNumber{Float64}(8.0, 8.0)
    num9 = GreyNumber{Float64}(9.0, 9.0)
    decmat = [
        num9 num7 num6 num7;
        num8 num7 num9 num6;
        num7 num8 num6 num6
    ]
    w = [4.0, 2, 6, 8]
    fns = makeminmax([maximum, maximum, maximum, maximum])
    result = topsis(decmat, w, fns)
    @test isapprox(result[1], GreyNumber{Float64}(0.38768, 0.38768), atol = atol)
    @test isapprox(result[2], GreyNumber{Float64}(0.65032, 0.65032), atol = atol)
    @test isapprox(result[3], GreyNumber{Float64}(0.08347, 0.08347), atol = atol)
end


text:
function parse_kw(ps)
    k = ps.t.kind
    if k == Tokens.IF
        return @default ps @closer ps block parse_if(ps)
    elseif k == Tokens.LET
        return @default ps @closer ps block parse_let(ps)
    elseif k == Tokens.TRY
        return @default ps @closer ps block parse_try(ps)
    elseif k == Tokens.FUNCTION
        return @default ps @closer ps block parse_function(ps)
    elseif k == Tokens.MACRO
        return @default ps @closer ps block parse_macro(ps)
    elseif k == Tokens.BEGIN
        return @default ps @closer ps block parse_begin(ps)
    elseif k == Tokens.QUOTE
        return @default ps @closer ps block parse_quote(ps)
    elseif k == Tokens.FOR
        return @default ps @closer ps block parse_for(ps)
    elseif k == Tokens.WHILE
        return @default ps @closer ps block parse_while(ps)
    elseif k == Tokens.BREAK
        return INSTANCE(ps)
    elseif k == Tokens.CONTINUE
        return INSTANCE(ps)
    elseif k == Tokens.IMPORT
        return parse_imports(ps)
    elseif k == Tokens.IMPORTALL
        return parse_imports(ps)
    elseif k == Tokens.USING
        return parse_imports(ps)
    elseif k == Tokens.EXPORT
        return parse_export(ps)
    elseif k == Tokens.MODULE ||  k == Tokens.BAREMODULE
        return @default ps @closer ps block parse_module(ps)
    elseif k == Tokens.CONST
        return @default ps parse_const(ps)
    elseif k == Tokens.GLOBAL
        return @default ps parse_global(ps)
    elseif k == Tokens.LOCAL
        return @default ps parse_local(ps)
    elseif k == Tokens.RETURN
        return @default ps parse_return(ps)
    elseif k == Tokens.END
        return parse_end(ps)
    elseif k == Tokens.ELSE || k == Tokens.ELSEIF || k == Tokens.CATCH || k == Tokens.FINALLY
        push!(ps.errors, Error((ps.t.startbyte:ps.t.endbyte) .+ 1 , "Unexpected end."))
        return ErrorToken(IDENTIFIER(ps))
    elseif k == Tokens.ABSTRACT
        return @default ps parse_abstract(ps)
    elseif k == Tokens.PRIMITIVE
        return @default ps parse_primitive(ps)
    elseif k == Tokens.TYPE
        return IDENTIFIER(ps)
    elseif k == Tokens.STRUCT
        return @default ps @closer ps block parse_struct(ps, false)
    elseif k == Tokens.MUTABLE
        return @default ps @closer ps block parse_mutable(ps)
    elseif k == Tokens.OUTER
        return IDENTIFIER(ps)
    end
end
# Prefix

function parse_const(ps::ParseState)
    kw = KEYWORD(ps)
    arg = parse_expression(ps)

    return EXPR{Const}(Any[kw, arg])
end

function parse_global(ps::ParseState)
    kw = KEYWORD(ps)
    arg = parse_expression(ps)

    return EXPR{Global}(Any[kw, arg])
end

function parse_local(ps::ParseState)
    kw = KEYWORD(ps)
    arg = parse_expression(ps)

    return EXPR{Local}(Any[kw, arg])
end

function parse_return(ps::ParseState)
    kw = KEYWORD(ps)
    args = closer(ps) ? NOTHING : parse_expression(ps)

    return EXPR{Return}(Any[kw, args])
end

# One line

@addctx :abstract function parse_abstract(ps::ParseState)
    # Switch for v0.6 compatability
    if ps.nt.kind == Tokens.TYPE
        kw1 = KEYWORD(ps)
        kw2 = KEYWORD(next(ps))
        sig = @closer ps block parse_expression(ps)
        ret = EXPR{Abstract}(Any[kw1, kw2, sig, accept_end(ps)])
    else
        kw = KEYWORD(ps)
        sig = parse_expression(ps)
        ret = EXPR{Abstract}(Any[kw, sig])
    end
    return ret
end

@addctx :primitive function parse_primitive(ps::ParseState)
    if ps.nt.kind == Tokens.TYPE
        kw1 = KEYWORD(ps)
        kw2 = KEYWORD(next(ps))
        sig = @closer ps ws @closer ps wsop parse_expression(ps)
        arg = @closer ps block parse_expression(ps)

        ret = EXPR{Primitive}(Any[kw1, kw2, sig, arg, accept_end(ps)])
    else
        ret = IDENTIFIER(ps)
    end
    return ret
end

function parse_imports(ps::ParseState)
    kw = KEYWORD(ps)
    kwt = is_import(kw) ? Import :
          is_importall(kw) ? ImportAll :
          Using
    tk = ps.t.kind

    arg = parse_dot_mod(ps)

    if ps.nt.kind != Tokens.COMMA && ps.nt.kind != Tokens.COLON
        ret = EXPR{kwt}(vcat(kw, arg))
    elseif ps.nt.kind == Tokens.COLON
        ret = EXPR{kwt}(vcat(kw, arg))
        push!(ret, OPERATOR(next(ps)))

        arg = parse_dot_mod(ps, true)
        append!(ret, arg)
        while ps.nt.kind == Tokens.COMMA
            accept_comma(ps, ret)
            arg = parse_dot_mod(ps, true)
            append!(ret, arg)
        end
    else
        ret = EXPR{kwt}(vcat(kw, arg))
        while ps.nt.kind == Tokens.COMMA
            accept_comma(ps, ret)
            arg = parse_dot_mod(ps)
            append!(ret, arg)
        end
    end

    return ret
end

function parse_export(ps::ParseState)
    args = Any[KEYWORD(ps)]
    append!(args, parse_dot_mod(ps))

    while ps.nt.kind == Tokens.COMMA
        push!(args, PUNCTUATION(next(ps)))
        arg = parse_dot_mod(ps)[1]
        push!(args, arg)
    end

    return EXPR{Export}(args)
end

# Block

@addctx :begin function parse_begin(ps::ParseState)
    sb = ps.t.startbyte
    kw = KEYWORD(ps)
    blockargs = parse_block(ps, Any[], (Tokens.END,), true)
    if isempty(blockargs)
        block = EXPR{Block}(blockargs, 0, 0)
    else
        fullspan = ps.nt.startbyte - sb - kw.fullspan
        block = EXPR{Block}(blockargs, fullspan, fullspan - last(blockargs).fullspan + last(blockargs).span)
    end
    # return EXPR{Begin}(Any[kw, EXPR{Block}(blockargs), accept_end(ps)])
    ender = accept_end(ps)
    fullspan1 = ps.nt.startbyte - sb
    return EXPR{Begin}(Any[kw, block, ender], fullspan1, fullspan1 - ender.fullspan + ender.span)
end

@addctx :quote function parse_quote(ps::ParseState)
    kw = KEYWORD(ps)
    blockargs = parse_block(ps)
    return EXPR{Quote}(Any[kw, EXPR{Block}(blockargs), accept_end(ps)])
end

@addctx :function function parse_function(ps::ParseState)
    kw = KEYWORD(ps)
    sig = @closer ps inwhere @closer ps ws parse_expression(ps)

    if sig isa EXPR{InvisBrackets} && !(sig.args[2] isa EXPR{TupleH})
        istuple = true
        sig = EXPR{TupleH}(sig.args)
    elseif sig isa EXPR{TupleH}
        istuple = true
    else
        istuple = false
    end

    while ps.nt.kind == Tokens.WHERE && ps.ws.kind != Tokens.NEWLINE_WS
        sig = @closer ps inwhere @closer ps ws parse_compound(ps, sig)
    end

    blockargs = parse_block(ps)

    if isempty(blockargs)
        if sig isa EXPR{Call} || sig isa WhereOpCall || (sig isa BinarySyntaxOpCall && !(is_exor(sig.arg1))) || istuple
            args = Any[sig, EXPR{Block}(blockargs)]
        else
            args = Any[sig]
        end
    else
        args = Any[sig, EXPR{Block}(blockargs)]
    end

    ret = EXPR{FunctionDef}(Any[kw])
    for a in args
        push!(ret, a)
    end
    accept_end(ps, ret)
    return ret
end

@addctx :macro function parse_macro(ps::ParseState)
    sb  = ps.t.startbyte
    kw = KEYWORD(ps)
    sig = @closer ps ws parse_expression(ps)
    sb1  = ps.nt.startbyte
    blockargs = parse_block(ps)

    if sig isa IDENTIFIER
        ender = accept_end(ps)
        fullspan1 = ps.nt.startbyte - sb
        ret = EXPR{Macro}(Any[kw, sig, ender], fullspan1, fullspan1 - ender.fullspan + ender.span)
    elseif isempty(blockargs)
        ender = accept_end(ps)
        fullspan1 = ps.nt.startbyte - sb
        ret = EXPR{Macro}(Any[kw, sig, EXPR{Block}([]), ender], fullspan1, fullspan1 - ender.fullspan + ender.span)
    else
        fullspan = ps.nt.startbyte - sb1
        block = EXPR{Block}(blockargs, fullspan, fullspan - last(blockargs).fullspan + last(blockargs).span)
        ender = accept_end(ps)
        fullspan1 = ps.nt.startbyte - sb
        ret = EXPR{Macro}(Any[kw, sig, block, ender], fullspan1, fullspan1 - ender.fullspan + ender.span)
    end
    return ret
end

# loops
@addctx :for function parse_for(ps::ParseState)
    sb  = ps.t.startbyte
    kw = KEYWORD(ps)
    ranges = parse_ranges(ps)
    sb1  = ps.nt.startbyte
    blockargs = parse_block(ps)

    # return EXPR{For}(Any[kw, ranges, EXPR{Block}(blockargs), accept_end(ps)])
    if isempty(blockargs)
        block = EXPR{Block}(blockargs, 0, 0)
    else
        fullspan = ps.nt.startbyte - sb1
        block = EXPR{Block}(blockargs, fullspan, fullspan - last(blockargs).fullspan + last(blockargs).span)
    end
    ender = accept_end(ps)
    fullspan1 = ps.nt.startbyte - sb
    return EXPR{For}(Any[kw, ranges, block, ender], fullspan1, fullspan1 - ender.fullspan + ender.span)
end

@addctx :while function parse_while(ps::ParseState)
    sb = ps.t.startbyte
    kw = KEYWORD(ps)
    cond = @closer ps ws parse_expression(ps)
    sb1 = ps.nt.startbyte
    blockargs = parse_block(ps)

    # return EXPR{While}(Any[kw, cond, EXPR{Block}(blockargs), accept_end(ps)])
    if isempty(blockargs)
        block = EXPR{Block}(blockargs, 0, 0)
    else
        fullspan = ps.nt.startbyte - sb1
        block = EXPR{Block}(blockargs, fullspan, fullspan - last(blockargs).fullspan + last(blockargs).span)
    end
    ender = accept_end(ps)
    fullspan1 = ps.nt.startbyte - sb
    return EXPR{While}(Any[kw, cond, block, ender], fullspan1, fullspan1 - ender.fullspan + ender.span)
end

# control flow

"""
    parse_if(ps, ret, nested=false, puncs=[])

Parse an `if` block.
"""
@addctx :if function parse_if(ps::ParseState, nested = false)
    # Parsing
    kw = KEYWORD(ps)
    if ps.ws.kind == NewLineWS || ps.ws.kind == SemiColonWS
        push!(ps.errors, Error((ps.ws.startbyte:ps.ws.endbyte) .+ 1 , "Missing conditional in if statement."))
        cond = ErrorToken()
    else
        cond = @closer ps ws parse_expression(ps)
    end
    ifblockargs = parse_block(ps, Any[], (Tokens.END, Tokens.ELSE, Tokens.ELSEIF))

    if nested
        ret = EXPR{If}(Any[cond, EXPR{Block}(ifblockargs)])
    else
        ret = EXPR{If}(Any[kw, cond, EXPR{Block}(ifblockargs)])
    end

    elseblockargs = Any[]
    if ps.nt.kind == Tokens.ELSEIF
        push!(ret, KEYWORD(next(ps)))
        push!(elseblockargs, parse_if(ps, true))
    end
    elsekw = ps.nt.kind == Tokens.ELSE
    if ps.nt.kind == Tokens.ELSE
        push!(ret, KEYWORD(next(ps)))
        parse_block(ps, elseblockargs)
    end

    # Construction
    if !(isempty(elseblockargs) && !elsekw)
        push!(ret, EXPR{Block}(elseblockargs))
    end
    !nested && accept_end(ps, ret)

    return ret
end

@addctx :let function parse_let(ps::ParseState)
    args = Any[KEYWORD(ps)]
    if !(ps.ws.kind == NewLineWS || ps.ws.kind == SemiColonWS)
        arg = @closer ps range @closer ps ws  parse_expression(ps)
        if ps.nt.kind == Tokens.COMMA
            arg = EXPR{Block}(Any[arg])
            while ps.nt.kind == Tokens.COMMA
                accept_comma(ps, arg)
                startbyte = ps.nt.startbyte
                nextarg = @closer ps comma @closer ps ws parse_expression(ps)
                push!(arg, nextarg)
            end
        end
        push!(args, arg)
    end

    blockargs = parse_block(ps)
    push!(args, EXPR{Block}(blockargs))
    accept_end(ps, args)

    return EXPR{Let}(args)
end

@addctx :try function parse_try(ps::ParseState)
    kw = KEYWORD(ps)
    ret = EXPR{Try}(Any[kw])

    tryblockargs = parse_block(ps, Any[], (Tokens.END, Tokens.CATCH, Tokens.FINALLY))
    push!(ret, EXPR{Block}(tryblockargs))

    #  catch block
    if ps.nt.kind == Tokens.CATCH
        next(ps)
        push!(ret, KEYWORD(ps))
        # catch closing early
        if ps.nt.kind == Tokens.FINALLY || ps.nt.kind == Tokens.END
            caught = FALSE
            catchblock = EXPR{Block}(Any[])
        else
            if ps.ws.kind == SemiColonWS || ps.ws.kind == NewLineWS
                caught = FALSE
            else
                caught = @closer ps ws parse_expression(ps)
            end

            catchblockargs = parse_block(ps, Any[], (Tokens.END, Tokens.FINALLY))
            if !(caught isa IDENTIFIER || caught == FALSE || (caught isa UnarySyntaxOpCall && caught.arg1 isa OPERATOR && caught.arg1.kind == Tokens.EX_OR))
                pushfirst!(catchblockargs, caught)
                caught = FALSE
            end
            catchblock = EXPR{Block}(catchblockargs)
        end
    else
        caught = FALSE
        catchblock = EXPR{Block}(Any[])
    end
    push!(ret, caught)
    push!(ret, catchblock)

    # finally block
    if ps.nt.kind == Tokens.FINALLY
        if isempty(catchblock.args)
            ret.args[4] = FALSE
        end
        push!(ret, KEYWORD(next(ps)))
        finallyblockargs = parse_block(ps)
        push!(ret, EXPR{Block}(finallyblockargs))
    end

    push!(ret, accept_end(ps))
    return ret
end

@addctx :do function parse_do(ps::ParseState, @nospecialize(ret))
    kw = KEYWORD(next(ps))

    args = EXPR{TupleH}(Any[])
    @closer ps comma @closer ps block while !closer(ps)
        a = parse_expression(ps)
        push!(args, a)
        if ps.nt.kind == Tokens.COMMA
            accept_comma(ps, args)
        end
    end

    blockargs = parse_block(ps)

    return EXPR{Do}(Any[ret, kw, args, EXPR{Block}(blockargs), accept_end(ps)])
end

# modules

@addctx :module function parse_module(ps::ParseState)
    sb = ps.t.startbyte
    kw = KEYWORD(ps)
    @assert kw.kind == Tokens.MODULE || kw.kind == Tokens.BAREMODULE # work around julia issue #23766
    if ps.nt.kind == Tokens.IDENTIFIER
        arg = IDENTIFIER(next(ps))
    else
        arg = @precedence ps 15 @closer ps ws parse_expression(ps)
    end
    sb1 = ps.nt.startbyte

    blockargs = parse_block(ps, Any[], (Tokens.END,), true)
    # block = EXPR{Block}(blockargs)

    # return EXPR{(is_module(kw) ? ModuleH : BareModule)}(Any[kw, arg, block, accept_end(ps)])
    if isempty(blockargs)
        block = EXPR{Block}(blockargs, 0, 0)
    else
        fullspan = ps.nt.startbyte - sb1
        block = EXPR{Block}(blockargs, fullspan, fullspan - last(blockargs).fullspan + last(blockargs).span)
    end
    ender = accept_end(ps)
    fullspan1 = ps.nt.startbyte - sb
    return EXPR{(is_module(kw) ? ModuleH : BareModule)}(Any[kw, arg, block, ender], fullspan1, fullspan1 - ender.fullspan + ender.span)
end

function parse_mutable(ps::ParseState)
    if ps.nt.kind == Tokens.STRUCT
        kw = KEYWORD(ps)
        next(ps)
        ret = parse_struct(ps, true)
        pushfirst!(ret, kw)
        update_span!(ret)
    else
        ret = IDENTIFIER(ps)
    end
    return ret
end

@addctx :struct function parse_struct(ps::ParseState, mutable)
    sb = ps.t.startbyte
    kw = KEYWORD(ps)
    sig = @closer ps ws parse_expression(ps)
    sb1 = ps.nt.startbyte
    blockargs = parse_block(ps)

    # return EXPR{mutable ? Mutable : Struct}(Any[kw, sig, EXPR{Block}(blockargs), accept_end(ps)])
    if isempty(blockargs)
        block = EXPR{Block}(blockargs, 0, 0)
    else
        fullspan = ps.nt.startbyte - sb1
        block = EXPR{Block}(blockargs, fullspan, fullspan - last(blockargs).fullspan + last(blockargs).span)
    end
    ender = accept_end(ps)
    fullspan1 = ps.nt.startbyte - sb
    return EXPR{mutable ? Mutable : Struct}(Any[kw, sig, block, ender], fullspan1, fullspan1 - ender.fullspan + ender.span)
end


text:
# phylo/manipulation.jl
# =====================
#
# Types and methods for phylogenetic trees.
#
# Part of the Bio.Phylo module.
#
# This file contains methods for manipulating the structure of phylogenies.
#
# This file is a part of BioJulia. License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

# Internal functions

"""
    unconnected_clades(tree::Phylogeny)

**Unexported function for internal use.**

Gets the internal clade vertices which are not connected to either parents or
children.
"""
function unconnected_clades(tree::Phylogeny)
    cv = clades(tree)
    idx = find(degree(tree.graph, cv) .== 0)
    return cv[idx]
end

"""
    subtree_roots(tree::Phylogeny)

**Unexported function for internal use.**

Gets clade vertices from the phylogeny which are roots of detached subphylogenies.
Such vertices have children, but no parent, and are not THE root of the
phylogeny (n + 1). Such subtrees are often created after pruning of trees.
"""
function subtree_roots(tree::Phylogeny)
    cv = clades(tree)
    idx = find(haschildren(tree, cv) & !hasparent(tree, cv))
    return cv[idx]
end

"""
    disconnect_root(tree::Phylogeny)

**Unexported function for internal use.**

Unconnects the root vertex from all of its children.
"""
function disconnect_root!(tree::Phylogeny)
    r = root(tree)
    for i in children(tree, r)
        destroy_branch!(tree, Edge(r, i))
    end
    return tree
end

"""
    delete!(tree::Phylogeny, vertex::Int, preserve_bl::Bool = false)

Delete a node from a phylogenetic tree.
"""
function Base.delete!(tree::Phylogeny, vertex::Int, preserve_bl::Bool = false)
    p = Phylogenies.parent(tree, vertex)
    # Delete the connection to parent but remember the branchlength of the
    # deleted branch.
    parentedge = Edge(p, vertex)
    lentoparent = preserve_bl ? branchlength(tree, parentedge) : 0.0
    destroy_branch!(tree, parentedge)
    for branch in child_branches(tree, vertex)
        newbl = branchlength(tree, branch) + lentoparent
        add_branch!(tree, Edge(parent, dst(branch)), newbl)
        destroy_branch!(tree, Edge(vertex, dst(branch)))
    end
    return tree
end


text:
# Defazio, Bach, Lacoste-Julien: SAGA: A fast incremental gradient method
# with support for non-strongly convex composite objectives.
# In: Advances in neural information processing systems, pp. 1646–1654 (2014).
#
# Reddi, Sra, Poczos, and Smola, "Proximal stochastic methods for nonsmooth
# nonconvex finite-sum optimization." In Advances in Neural Information
# Processing Systems, pp. 1145–1153 (2016).
#
# Schmidt, Le Roux, Bach, "Minimizing finite sums with the stochastic average gradient"
# Mathematical Programming, 162(1-2), 83-112 (2017).
#

using LinearAlgebra
using ProximalOperators
using ProximalAlgorithms.IterationTools
using Printf
using Base.Iterators
using Random

export solution

include("SAGA_basic.jl")

struct SAGA{R<:Real}
    γ::Maybe{R}
    maxit::Int
    verbose::Bool
    freq::Int
    SAG_flag::Bool
    function SAGA{R}(;
        γ::Maybe{R} = nothing,
        maxit::Int = 10000,
        verbose::Bool = false,
        freq::Int = 1000,
        SAG_flag::Bool = false,
    ) where {R}
        @assert γ === nothing || γ > 0
        @assert maxit > 0
        @assert freq > 0
        new(γ, maxit, verbose, freq, SAG_flag)
    end
end

function (solver::SAGA{R})(
    x0::AbstractArray{C};
    F = nothing,
    g = ProximalOperators.Zero(),
    L = nothing,
    N = N,
) where {R,C<:RealOrComplex{R}}

    stop(state::SAGA_basic_state) = false
    disp(it, state) = @printf "%5d | %.3e  \n" it state.γ

    F === nothing && (F = fill(ProximalOperators.Zero(), (N,)))

    # dispatching the structure
    iter = SAGA_basic_iterable(F, g, x0, N, L, solver.γ, solver.SAG_flag)
    iter = take(halt(iter, stop), solver.maxit)
    iter = enumerate(iter)
    num_iters, state_final = nothing, nothing
    for (it_, state_) in iter  # unrolling the iterator
        # see https://docs.julialang.org/en/v1/manual/interfaces/index.html
        if solver.verbose && mod(it_, solver.freq) == 0
            disp(it_, state_)
        end
        num_iters, state_final = it_, state_
    end
    if solver.verbose && mod(num_iters, solver.freq) !== 0
        disp(num_iters, state_final)
    end # for the final iteration
    return solution(state_final), num_iters
end

"""
    SAGA([γ, maxit, verbose, freq])

Instantiate the SAGA algorithm  for solving convex optimization problems of the form

    minimize 1/N sum_{i=1}^N f_i(x) + g(x)

If `solver = SAGA(args...)`, then the above problem is solved with

	solver(x0, [F, g, N, L])

where F is an array containing f_i's, x0 is the initial point, and L is an array of
smoothness moduli of f_i's; it is optional when γ is provided.

Optional keyword arguments are:
* `γ`: stepsize
* `L`: an array of smoothness moduli of f_i's
* `maxit::Integer` (default: `10000`), maximum number of iterations to perform.
* `verbose::Bool` (default: `true`), whether or not to print information during the iterations.
* `freq::Integer` (default: `100`), frequency of verbosity.

References:

[1] Defazio, Bach, Lacoste-Julien, "SAGA: A fast incremental gradient method
with support for non-strongly convex composite objectives"
In: Advances in neural information processing systems, pp. 1646–1654 (2014).

[2] Reddi, Sra, Poczos, and Smola, "Proximal stochastic methods for nonsmooth
nonconvex finite-sum optimization" In Advances in Neural Information
Processing Systems, pp. 1145–1153 (2016).
"""

SAGA(::Type{R}; kwargs...) where {R} = SAGA{R}(; kwargs...)
SAGA(; kwargs...) = SAGA(Float64; kwargs...)

"""
If `solver = SAGA(args...)`, then

    itr = iterator(solver, x0, [F, g, N, L])

is an iterable object. Note that [maxit, verbose, freq] fields of the solver are ignored here.

The solution at any given state can be obtained using solution(state), e.g.,
for state in Iterators.take(itr, maxit)
    # do something using solution(state)
end

See https://docs.julialang.org/en/v1/manual/interfaces/index.html
and https://docs.julialang.org/en/v1/base/iterators/ for a list of iteration utilities
"""

function iterator(
    solver::SAGA{R},
    x0::AbstractArray{C};
    F = nothing,
    g = ProximalOperators.Zero(),
    L = nothing,
    N = N,
) where {R,C<:RealOrComplex{R}}

    F === nothing && (F = fill(ProximalOperators.Zero(), (N,)))
    # dispatching the iterator
    iter = SAGA_basic_iterable(F, g, x0, N, L, solver.γ, solver.SAG_flag)

    return iter
end

"""
    SAG([γ, maxit, verbose, freq])

Instantiate the SAG algorithm  for solving convex optimization problems of the form

    minimize 1/N sum_{i=1}^N f_i(x) + g(x)

If `solver = SAG(args...)`, then the above problem is solved with

    solver(x0, [F, g, N, L])

where F is an array containing f_i's, x0 is the initial point, and L is an array of
smoothness moduli of f_i's; it is optional when γ is provided.

Optional keyword arguments are:
* `γ`: stepsize
* `L`: an array of smoothness moduli of f_i's
* `maxit::Integer` (default: `10000`), maximum number of iterations to perform.
* `verbose::Bool` (default: `true`), whether or not to print information during the iterations.
* `freq::Integer` (default: `100`), frequency of verbosity.

References:

[1] Schmidt, Le Roux, Bach, "Minimizing finite sums with the stochastic average gradient"
Mathematical Programming, 162(1-2), 83-112 (2017).

"""

"""
If `solver = SAG(args...)`, then

    itr = iterator(solver, x0, [F, g, N, L])

is an iterable object. Note that [maxit, verbose, freq] fields of the solver are ignored here.

The solution at any given state can be obtained using solution(state), e.g.,
for state in Iterators.take(itr, maxit)
    # do something using solution(state)
end

See https://docs.julialang.org/en/v1/manual/interfaces/index.html
and https://docs.julialang.org/en/v1/base/iterators/ for a list of iteration utilities
"""

SAG(::Type{R}; kwargs...) where {R} = SAGA{R}(; kwargs..., SAG_flag = true)
SAG(; kwargs...) = SAG(Float64; kwargs...)

max_issues_repo_name: kul-forbes/CIAOAlgorithms

max_forks_repo_name: kul-forbes/CIAOAlgorithms


text:
function transform(input::AbstractDict)
    output = Dict{Char,Int}()
    for (k, v) in input
        for vi in v
            output[lowercase(vi)] = k
        end
    end
    return output
end

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
include("EllCrv/EllCrv.jl")
include("EllCrv/Misc.jl")
include("EllCrv/QQ.jl")
include("EllCrv/Finite.jl")
include("EllCrv/Heights.jl")
include("EllCrv/MordellWeilRank.jl")
include("EllCrv/Isogeny.jl")


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# SubString and RevString types

## substrings reference original strings ##

"""
    SubString(s::AbstractString, i::Integer, j::Integer=endof(s))
    SubString(s::AbstractString, r::UnitRange{<:Integer})

Like [`getindex`](@ref), but returns a view into the parent string `s`
within range `i:j` or `r` respectively instead of making a copy.

# Examples
```jldoctest
julia> SubString("abc", 1, 2)
"ab"

julia> SubString("abc", 1:2)
"ab"

julia> SubString("abc", 2)
"bc"
```
"""
struct SubString{T<:AbstractString} <: AbstractString
    string::T
    offset::Int
    endof::Int

    function SubString{T}(s::T, i::Int, j::Int) where T<:AbstractString
        i > j && return new(s, i - 1, 0) # always allow i > j as it is consistent with getindex
        isvalid(s, i) || throw(BoundsError(s, i))
        isvalid(s, j) || throw(BoundsError(s, j))
        new(s, i-1, j-i+1)
    end
end

SubString(s::T, i::Int, j::Int) where {T<:AbstractString} = SubString{T}(s, i, j)
SubString(s::AbstractString, i::Integer, j::Integer=endof(s)) = SubString(s, Int(i), Int(j))
SubString(s::AbstractString, r::UnitRange{<:Integer}) = SubString(s, first(r), last(r))

function SubString(s::SubString, i::Int, j::Int)
    # always allow i > j as it is consistent with getindex
    i > j && return SubString(s.string, s.offset + i, s.offset + j)
    i >= 1 || throw(BoundsError(s, i))
    j <= endof(s) || throw(BoundsError(s, j))
    SubString(s.string, s.offset + i, s.offset + j)
end

SubString(s::AbstractString) = SubString(s, 1, endof(s))
SubString{T}(s::T) where {T<:AbstractString} = SubString{T}(s, 1, endof(s))

String(p::SubString{String}) =
    unsafe_string(pointer(p.string, p.offset+1), nextind(p, p.endof)-1)

sizeof(s::SubString{String}) = s.endof == 0 ? 0 : nextind(s, s.endof) - 1

# TODO: length(s::SubString) = ??
# default implementation will work but it's slow
# can this be delegated efficiently somehow?
# that may require additional string interfaces
function length(s::SubString{String})
    return s.endof==0 ? 0 : Int(ccall(:u8_charnum, Csize_t, (Ptr{UInt8}, Csize_t),
                                      pointer(s), nextind(s, s.endof) - 1))
end

function next(s::SubString, i::Int)
    if i < 1 || i > s.endof
        throw(BoundsError(s, i))
    end
    c, i = next(s.string, i+s.offset)
    c, i-s.offset
end

function getindex(s::SubString, i::Int)
    if i < 1 || i > s.endof
        throw(BoundsError(s, i))
    end
    getindex(s.string, i+s.offset)
end

endof(s::SubString) = s.endof

function isvalid(s::SubString, i::Integer)
    return (start(s) <= i <= endof(s)) && isvalid(s.string, s.offset+i)
end

nextind(s::SubString, i::Integer) = nextind(s.string, i+s.offset)-s.offset
prevind(s::SubString, i::Integer) = prevind(s.string, i+s.offset)-s.offset

function getindex(s::AbstractString, r::UnitRange{Int})
    checkbounds(s, r) || throw(BoundsError(s, r))
    SubString(s, first(r), last(r))
end

function cmp(a::SubString{String}, b::SubString{String})
    na = sizeof(a)
    nb = sizeof(b)
    c = ccall(:memcmp, Int32, (Ptr{UInt8}, Ptr{UInt8}, UInt),
              pointer(a), pointer(b), min(na,nb))
    c < 0 ? -1 : c > 0 ? +1 : cmp(na,nb)
end

# don't make unnecessary copies when passing substrings to C functions
cconvert(::Type{Ptr{UInt8}}, s::SubString{String}) = s
cconvert(::Type{Ptr{Int8}}, s::SubString{String}) = s
function unsafe_convert(::Type{Ptr{R}}, s::SubString{String}) where R<:Union{Int8, UInt8}
    convert(Ptr{R}, pointer(s.string)) + s.offset
end

## reversed strings without data movement ##

struct RevString{T<:AbstractString} <: AbstractString
    string::T
end

endof(s::RevString) = endof(s.string)
length(s::RevString) = length(s.string)
sizeof(s::RevString) = sizeof(s.string)

function next(s::RevString, i::Int)
    n = endof(s); j = n-i+1
    (s.string[j], n-prevind(s.string,j)+1)
end

"""
    reverse(s::AbstractString) -> AbstractString

Reverses a string.

Technically, this function reverses the codepoints in a string, and its
main utility is for reversed-order string processing, especially for reversed
regular-expression searches.  See also [`reverseind`](@ref) to convert indices
in `s` to indices in `reverse(s)` and vice-versa, and [`graphemes`](@ref)
to operate on user-visible "characters" (graphemes) rather than codepoints.

# Examples
```jldoctest
julia> reverse("JuliaLang")
"gnaLailuJ"

julia> reverse("ax̂e") # combining characters can lead to surprising results
"êxa"

julia> join(reverse(collect(graphemes("ax̂e")))) # reverses graphemes
"ex̂a"
```
"""
reverse(s::AbstractString) = RevString(s)
reverse(s::RevString) = s.string

## reverse an index i so that reverse(s)[i] == s[reverseind(s,i)]

"""
    reverseind(v, i)

Given an index `i` in [`reverse(v)`](@ref), return the corresponding index in `v` so that
`v[reverseind(v,i)] == reverse(v)[i]`. (This can be nontrivial in cases where `v` contains
non-ASCII characters.)

# Examples
```jldoctest
julia> r = reverse("Julia")
"ailuJ"

julia> for i in 1:length(r)
           print(r[reverseind("Julia", i)])
       end
Julia
```
"""
reverseind(s::AbstractString, i) = chr2ind(s, length(s) + 1 - ind2chr(reverse(s), i))
reverseind(s::RevString, i::Integer) = endof(s) - i + 1
reverseind(s::SubString{String}, i::Integer) =
    reverseind(s.string, nextind(s.string, endof(s.string))-s.offset-s.endof+i-1) - s.offset

"""
    repeat(s::AbstractString, r::Integer)

Repeat a string `r` times. This can equivalently be accomplished by calling [`s^r`](@ref ^).

# Examples
```jldoctest
julia> repeat("ha", 3)
"hahaha"
```
"""
repeat(s::AbstractString, r::Integer) = repeat(convert(String, s), r)

"""
    ^(s::Union{AbstractString,Char}, n::Integer)

Repeat a string or character `n` times.
The [`repeat`](@ref) function is an alias to this operator.

# Examples
```jldoctest
julia> "Test "^3
"Test Test Test "
```
"""
(^)(s::Union{AbstractString,Char}, r::Integer) = repeat(s,r)

pointer(x::SubString{String}) = pointer(x.string) + x.offset
pointer(x::SubString{String}, i::Integer) = pointer(x.string) + x.offset + (i-1)

max_stars_repo_name: xorJane/julia

max_issues_repo_name: xorJane/julia

max_forks_repo_name: xorJane/julia


text:
using Test
include(joinpath(dirname(@__FILE__), "../../src/includeall.jl"))

d = Dict()
default_slm_kwargs!(d)
d1 = deepcopy(d)
d1[:increasing] = true

@testset "Default SLM keyword arguments" begin
    for (i, di) in enumerate([d, d1])
        @test di[:degree] == 3
        @test di[:scaling] == true
        @test di[:knots] == 6
        @test di[:C2] == true
        @test di[:λ] == 1e-4
        if i == 1
            @test di[:increasing] == false
        else
            @test di[:increasing] == true
        end
        @test di[:decreasing] == false
        @test isempty(di[:increasing_intervals])
        @test isempty(di[:decreasing_intervals])
        @test isa(di[:increasing_intervals], Matrix{Float64})
        @test isa(di[:decreasing_intervals], Matrix{Float64})
        @test di[:concave_up] == false
        @test di[:concave_down] == false
        @test isempty(di[:concave_up_intervals])
        @test isempty(di[:concave_down_intervals])
        @test isa(di[:concave_up_intervals], Matrix{Float64})
        @test isa(di[:concave_down_intervals], Matrix{Float64})
        @test isnan(di[:left_value])
        @test isnan(di[:right_value])
        @test isnan(di[:min_value])
        @test isnan(di[:max_value])
        @test di[:min_max_sample_points] == [.017037, .066987, .1465, .25, .37059,
                                             .5, .62941, .75, .85355, .93301, .98296]
        @test isempty(di[:init])
        @test di[:use_sparse] == false
        @test di[:use_lls]
    end
end


text:
using RecipesBase, StatsBase, Measures

# # Plot Recipes
@recipe function f(output::SisMcmcOutput)
    model = output.model
    sample = output.sample
    x = map(x -> model.dist(model.mode, x), sample)
    xguide --> "Sample"
    yguide --> "Distance from Mode"
    legend --> false
    size --> (800, 300)
    margin --> 5mm
    x
end

@recipe function f(output::SimMcmcOutput)
    model = output.model
    sample = output.sample
    x = map(x -> model.dist(model.mode, x), sample)
    xguide --> "Sample"
    yguide --> "Distance from Mode"
    legend --> false
    size --> (800, 300)
    margin --> 5mm
    x
end

@recipe function f(output::SpfMcmcOutput)
    model = output.model
    sample = output.sample
    x = map(x -> model.dist(model.mode, x), sample)
    xguide --> "Sample"
    yguide --> "Distance from Mode"
    legend --> false
    size --> (800, 300)
    margin --> 5mm
    x
end

@userplot LengthPlot

@recipe function f(h::LengthPlot)
    input = h.args
    @assert length(input) == 1 "Accepts only a single argument"
    @assert typeof(input[1]) == SpfMcmcOutput "Input must be of type SpfMcmcOutput"
    size --> (800, 400)
    legend --> false
    data = length.(input[1].sample)
    @series begin
        seriestype := :bar
        yguide --> "Proportion"
        xguide --> "Path Length"
        xticks --> 1:maximum(data)
        xlims --> [0, maximum(data)+0.5]
        proportionmap(data)
    end

end

@userplot SummaryPlot
@recipe function f(h::SummaryPlot)
    output = h.args  # A mcmc_output type
    model = output[1].model
    sample = output[1].sample
    layout := (3,1)
    legend --> false
    xguide := "Sample"
    size --> (800, 600)
    margin --> 5mm

    @series begin
        seriestype := :line
        yguide --> "Distance from Mode"
        map(x -> model.dist(model.mode, x), sample)
    end
    @series begin
        seriestype := :line
        yguide --> "Num. Interactions"
        map(length, sample)
    end
    @series begin
        seriestype := :line
        yguide --> "Mean Interaction Length"
        map(x -> mean(length.(x)), sample)
    end
end


text:
using Compat.SparseArrays

const ZeroCones = Union{MOI.EqualTo, MOI.Zeros}
const LPCones = Union{MOI.GreaterThan, MOI.LessThan, MOI.Nonnegatives, MOI.Nonpositives}

_dim(s::MOI.AbstractScalarSet) = 1
_dim(s::MOI.AbstractVectorSet) = MOI.dimension(s)

# Computes cone dimensions
constroffset(cone::Cone, ci::CI{<:MOI.AbstractFunction, <:ZeroCones}) = ci.value
function _allocateconstraint!(cone::Cone, f, s::ZeroCones)
    ci = cone.f
    cone.f += _dim(s)
    ci
end
constroffset(cone::Cone, ci::CI{<:MOI.AbstractFunction, <:LPCones}) = cone.f + ci.value
function _allocateconstraint!(cone::Cone, f, s::LPCones)
    ci = cone.l
    cone.l += _dim(s)
    ci
end
constroffset(cone::Cone, ci::CI{<:MOI.AbstractFunction, <:MOI.SecondOrderCone}) = cone.f + cone.l + ci.value
function _allocateconstraint!(cone::Cone, f, s::MOI.SecondOrderCone)
    push!(cone.qa, s.dimension)
    ci = cone.q
    cone.q += _dim(s)
    ci
end
constroffset(cone::Cone, ci::CI{<:MOI.AbstractFunction, <:MOI.PositiveSemidefiniteConeTriangle}) = cone.f + cone.l + cone.q + ci.value
function _allocateconstraint!(cone::Cone, f, s::MOI.PositiveSemidefiniteConeTriangle)
    push!(cone.sa, s.dimension)
    ci = cone.s
    cone.s += _dim(s)
    ci
end
constroffset(cone::Cone, ci::CI{<:MOI.AbstractFunction, <:MOI.ExponentialCone}) = cone.f + cone.l + cone.q + cone.s + ci.value
function _allocateconstraint!(cone::Cone, f, s::MOI.ExponentialCone)
    ci = 3cone.ep
    cone.ep += 1
    ci
end
constroffset(optimizer::SCSOptimizer, ci::CI) = constroffset(optimizer.cone, ci::CI)
MOIU.canallocateconstraint(::SCSOptimizer, ::Type{<:SF}, ::Type{<:SS}) = true
function MOIU.allocateconstraint!(optimizer::SCSOptimizer, f::F, s::S) where {F <: MOI.AbstractFunction, S <: MOI.AbstractSet}
    CI{F, S}(_allocateconstraint!(optimizer.cone, f, s))
end

# Vectorized length for matrix dimension n
sympackedlen(n) = (n*(n+1)) >> 1
# Matrix dimension for vectorized length n
sympackeddim(n) = div(isqrt(1+8n) - 1, 2)
trimap(i::Integer, j::Integer) = i < j ? trimap(j, i) : div((i-1)*i, 2) + j
trimapL(i::Integer, j::Integer, n::Integer) = i < j ? trimapL(j, i, n) : i + div((2n-j) * (j-1), 2)
function _sympackedto(x, n, mapfrom, mapto)
    @assert length(x) == sympackedlen(n)
    y = similar(x)
    for i in 1:n, j in 1:i
        y[mapto(i, j)] = x[mapfrom(i, j)]
    end
    y
end
sympackedLtoU(x, n=sympackeddim(length(x))) = _sympackedto(x, n, (i, j) -> trimapL(i, j, n), trimap)
sympackedUtoL(x, n=sympackeddim(length(x))) = _sympackedto(x, n, trimap, (i, j) -> trimapL(i, j, n))

function sympackedUtoLidx(x::AbstractVector{<:Integer}, n)
    y = similar(x)
    map = sympackedLtoU(1:sympackedlen(n), n)
    for i in eachindex(y)
        y[i] = map[x[i]]
    end
    y
end

# Build constraint matrix
scalecoef(rows, coef, minus, s) = minus ? -coef : coef
scalecoef(rows, coef, minus, s::Union{MOI.LessThan, Type{<:MOI.LessThan}, MOI.Nonpositives, Type{MOI.Nonpositives}}) = minus ? coef : -coef
function _scalecoef(rows, coef, minus, d, rev)
    scaling = minus ? -1 : 1
    scaling2 = rev ? scaling / sqrt(2) : scaling * sqrt(2)
    output = copy(coef)
    diagidx = IntSet()
    for i in 1:d
        push!(diagidx, trimap(i, i))
    end
    for i in 1:length(output)
        if rows[i] in diagidx
            output[i] *= scaling
        else
            output[i] *= scaling2
        end
    end
    output
end
function scalecoef(rows, coef, minus, s::MOI.PositiveSemidefiniteConeTriangle)
    _scalecoef(rows, coef, minus, s.dimension, false)
end
function scalecoef(rows, coef, minus, ::Type{MOI.PositiveSemidefiniteConeTriangle})
    _scalecoef(rows, coef, minus, sympackeddim(length(rows)), true)
end
_varmap(f) = map(vi -> vi.value, f.variables)
_constant(s::MOI.EqualTo) = s.value
_constant(s::MOI.GreaterThan) = s.lower
_constant(s::MOI.LessThan) = s.upper
constrrows(::MOI.AbstractScalarSet) = 1
constrrows(s::MOI.AbstractVectorSet) = 1:MOI.dimension(s)
constrrows(optimizer::SCSOptimizer, ci::CI{<:MOI.AbstractScalarFunction, <:MOI.AbstractScalarSet}) = 1
constrrows(optimizer::SCSOptimizer, ci::CI{<:MOI.AbstractVectorFunction, <:MOI.AbstractVectorSet}) = 1:optimizer.cone.nrows[constroffset(optimizer, ci)]
MOIU.canloadconstraint(::SCSOptimizer, ::Type{<:SF}, ::Type{<:SS}) = true
MOIU.loadconstraint!(optimizer::SCSOptimizer, ci, f::MOI.SingleVariable, s) = MOIU.loadconstraint!(optimizer, ci, MOI.ScalarAffineFunction{Float64}(f), s)
function MOIU.loadconstraint!(optimizer::SCSOptimizer, ci, f::MOI.ScalarAffineFunction, s::MOI.AbstractScalarSet)
    a = sparsevec(_varmap(f), f.coefficients)
    # sparsevec combines duplicates with + but does not remove zeros created so we call dropzeros!
    dropzeros!(a)
    offset = constroffset(optimizer, ci)
    row = constrrows(s)
    i = offset + row
    # The SCS format is b - Ax ∈ cone
    # so minus=false for b and minus=true for A
    setconstant = _constant(s)
    optimizer.cone.setconstant[offset] = setconstant
    constant = f.constant - setconstant
    optimizer.data.b[i] = scalecoef(row, constant, false, s)
    append!(optimizer.data.I, fill(i, length(a.nzind)))
    append!(optimizer.data.J, a.nzind)
    append!(optimizer.data.V, scalecoef(row, a.nzval, true, s))
end
MOIU.loadconstraint!(optimizer::SCSOptimizer, ci, f::MOI.VectorOfVariables, s) = MOIU.loadconstraint!(optimizer, ci, MOI.VectorAffineFunction{Float64}(f), s)
orderval(val, s) = val
function orderval(val, s::MOI.PositiveSemidefiniteConeTriangle)
    sympackedUtoL(val, s.dimension)
end
orderidx(idx, s) = idx
function orderidx(idx, s::MOI.PositiveSemidefiniteConeTriangle)
    sympackedUtoLidx(idx, s.dimension)
end
function MOIU.loadconstraint!(optimizer::SCSOptimizer, ci, f::MOI.VectorAffineFunction, s::MOI.AbstractVectorSet)
    A = sparse(f.outputindex, _varmap(f), f.coefficients)
    # sparse combines duplicates with + but does not remove zeros created so we call dropzeros!
    dropzeros!(A)
    colval = zeros(Int, length(A.rowval))
    for col in 1:A.n
        colval[A.colptr[col]:(A.colptr[col+1]-1)] = col
    end
    @assert !any(iszero.(colval))
    offset = constroffset(optimizer, ci)
    rows = constrrows(s)
    optimizer.cone.nrows[offset] = length(rows)
    i = offset + rows
    # The SCS format is b - Ax ∈ cone
    # so minus=false for b and minus=true for A
    optimizer.data.b[i] = scalecoef(rows, orderval(f.constant, s), false, s)
    append!(optimizer.data.I, offset + orderidx(A.rowval, s))
    append!(optimizer.data.J, colval)
    append!(optimizer.data.V, scalecoef(A.rowval, A.nzval, true, s))
end

function MOIU.allocatevariables!(optimizer::SCSOptimizer, nvars::Integer)
    optimizer.cone = Cone()
    VI.(1:nvars)
end

function MOIU.loadvariables!(optimizer::SCSOptimizer, nvars::Integer)
    cone = optimizer.cone
    m = cone.f + cone.l + cone.q + cone.s + 3cone.ep + cone.ed
    I = Int[]
    J = Int[]
    V = Float64[]
    b = zeros(m)
    c = zeros(nvars)
    optimizer.data = Data(m, nvars, I, J, V, b, 0., c)
end

MOIU.canallocate(::SCSOptimizer, ::MOI.ObjectiveSense) = true
function MOIU.allocate!(optimizer::SCSOptimizer, ::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    optimizer.maxsense = sense == MOI.MaxSense
end
MOIU.canallocate(::SCSOptimizer, ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}) = true
function MOIU.allocate!(::SCSOptimizer, ::MOI.ObjectiveFunction, ::MOI.ScalarAffineFunction) end

MOIU.canload(::SCSOptimizer, ::MOI.ObjectiveSense) = true
function MOIU.load!(::SCSOptimizer, ::MOI.ObjectiveSense, ::MOI.OptimizationSense) end
MOIU.canload(::SCSOptimizer, ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}) = true
function MOIU.load!(optimizer::SCSOptimizer, ::MOI.ObjectiveFunction, f::MOI.ScalarAffineFunction)
    c0 = full(sparsevec(_varmap(f), f.coefficients, optimizer.data.n))
    optimizer.data.objconstant = f.constant
    optimizer.data.c = optimizer.maxsense ? -c0 : c0
end

function MOI.optimize!(optimizer::SCSOptimizer)
    cone = optimizer.cone
    m = optimizer.data.m
    n = optimizer.data.n
    A = sparse(optimizer.data.I, optimizer.data.J, optimizer.data.V)
    b = optimizer.data.b
    objconstant = optimizer.data.objconstant
    c = optimizer.data.c
    optimizer.data = nothing # Allows GC to free optimizer.data before A is loaded to SCS
    sol = SCS_solve(SCS.Indirect, m, n, A, b, c, cone.f, cone.l, cone.qa, cone.sa, cone.ep, cone.ed, cone.p)
    ret_val = sol.ret_val
    primal = sol.x
    dual = sol.y
    slack = sol.s
    objval = (optimizer.maxsense ? -1 : 1) * dot(c, primal) + objconstant
    optimizer.sol = Solution(ret_val, primal, dual, slack, objval)
end


text:
"""
# module Doolittle

# Examples

```jldoctest
julia>
```
"""
module Doolittle

using LinearAlgebra: I, kron, UpperTriangular

using APMAE4001.Matrix: issquare
using APMAE4001.LUFactorization.Substitution

export doolittle_lu,
    polysolve

function row_vector(n::Int, k::Int)
    v = zeros(n)
    v[k] = one(eltype(v))
    return v
end

"""
    elimination_factors(n::Int, k::Int, A::AbstractMatrix)

Generate a vector of factors to eliminate the `k`th column.
Combine this function with `row_vector` to to generate an general elementary elimination matrix.

# Arguments
- `n::Integer`: the size of the matrix `A`.
- `k::Integer`: the step count in the elimination process. It must satisfies `1 ≤ k ≤ n - 1`.
- `A::AbstractMatrix`: the matrix from `(k - 1)`th step of the elimination process.
"""
function elimination_factors(n::Int, k::Int, A::AbstractMatrix)
    v = zeros(promote_type(eltype(A), Float64), n)
    for j in (k + 1):n
        v[j] = A[j, k] / A[k, k]  # elimination factor
    end
    v
end

"""
    general_elementary_elimination_matrix(lk::AbstractVector, ek::AbstractVector)

Elimination matrix ``E_{k}`` do row operations to eliminate the `k`th column.
"""
general_elementary_elimination_matrix(lk::AbstractVector, ek::AbstractVector) = I - kron(lk, transpose(ek))  # `I` is the identity matrix.

"""
    doolittle_lu(A::AbstractMatrix; include_l::Bool=true, include_e::Bool=false)

Referenced from [Wikipedia "LU decomposition" page](https://en.wikipedia.org/wiki/LU_decomposition#Doolittle_algorithm).
"""
function doolittle_lu(A::AbstractMatrix; include_l::Bool = true, include_e::Bool = false)
    issquare(A) || throw(DimensionMismatch("The matrix `A` is not a square matrix!"))
    n = size(A, 1)

    intermediate_elementary_elimination_matrices = []
    for k in 1:(n - 1)
        Ek = general_elementary_elimination_matrix(elimination_factors(n, k, A), row_vector(n, k))
        pushfirst!(intermediate_elementary_elimination_matrices, Ek)  # Keep this order!
        A = Ek * A
    end
    # Now `A` is the U matrix.
    E = prod(*, intermediate_elementary_elimination_matrices)
    include_e && return E, A
    include_l && return inv(E), A  # L matrix is `inv(E)`.
    return A
end

function polysolve(A::AbstractMatrix, b::AbstractVector)
    E, U = doolittle_lu(A, include_e = true)
    substitution(UpperTriangular(U), E * b)
end

end


text:
#Template for a perfect model example:

# Import external modules

# Import CES modules
using EnsembleKalmanProcesses.EnsembleKalmanProcessModule
using EnsembleKalmanProcesses.Observations
using EnsembleKalmanProcesses.ParameterDistributionStorage
using EnsembleKalmanProcesses.DataStorage

# include model or data
include("GModel.jl")

##
## Create your perfect model "truth"
##

true_params = 0.0

GModel.set_properties(initial_conditions, boundary_conditions, fixed_params)

# Obtain samples of truth
truth_data = GModel.run_G_ensemble(true_params, GModel_settings, ensemble_size = 1)

# Calculate the observational noise on the data.
noise_covariance = prescribed_observational_noise_covariance

##
## Store the observation and covariance
##

truth = Observations.Obs(truth_data, observational_noise_covariance, data_names)
truth_sample = truth.samples[1]

##
## Now define the priors:
##
prior_dict = Dict("distribution" => Parameterized(Normal(0, 1)), "constraint" => bounded(1, 2), "name" => "")
priors = ParameterDistribution(prior_dict)

##
## Perform Ensemble_Kalman_Inversion
##

# initial ensemble of parameters
params = construct_initial_ensemble(priors, N_ens; rng_seed = rng_seed)

# constructor for the EKI object
ekiobj =
    EnsembleKalmanProcessModule.EnsembleKalmanProcess(initial_params, truth_sample, truth.obs_noise_cov, Inversion())

#outer EKI loop
for i in 1:N_iterations
    # theta^n_i -> G(theta^n_i)
    physical_params_i = transform_unconstrained_to_constrained(priors, params)
    g_ens = GModel.run_G_ensemble(physical_params_i, lorenz_settings_G)

    # G(theta^n_i) -> theta^{n+1}_i
    EnsembleKalmanProcessModule.update_ensemble!(ekiobj, g_ens)

    err[i] = get_error(ekiobj)[end] #mean((params_true - mean(params_i,dims=2)).^2)

end

##
## store results
##

u_stored = get_u(ekiobj, return_array = false)
g_stored = get_g(ekiobj, return_array = false)

#solution:
final_params = get_u_final(ekiobj)
solution = mean(final_params, dims = 2)

physical_final_params = transform_unconstrained_to_constrained(priors, final_ensemble)
physical_solution = mean(physical_final_soluation)


text:
# Add PATH (adapt!)
push!(LOAD_PATH, "/Volumes/cfranken/code/gitHub/Land/src/");
#push!(LOAD_PATH, "/Volumes/cfranken/code/gitHub/Land/src/Utils/");
#push!(LOAD_PATH, "/Volumes/cfranken/code/gitHub/LSM-SPAM/src/tools/");

# Use Plots:
using Plots
#using PhotoStructs
pyplot()

using Revise
using Land
using BenchmarkTools
using Statistics
using Parameters

@unpack CanopyRTMod = Land
@unpack FT,leafbio, canopy, angles, canOpt, canRad,sunRad,soil = Land.CanopyRTMod

wl = CanopyRTMod.wl;
wle = CanopyRTMod.wle;
wlf = CanopyRTMod.wlf;

#@time CanopyRTMod.optis.nr[3:20] .= 2

leaf = leafbio{FT, length(wl), length(wle), length(wlf),length(wle)*length(wlf)}();

using Parameters
@unpack N, Cab, Car, Ant, Cs, Cw, Cm, ρ_SW, τ_SW = CanopyRTMod.leaf;

@btime CanopyRTMod.fluspect!(CanopyRTMod.leaf, CanopyRTMod.optis);

contourf(wle, wlf, CanopyRTMod.leaf.Mb)

contourf(wle, wlf, CanopyRTMod.leaf.Mf)

arrayOfLeaves = Array{leafbio{FT,length(wl), length(wle), length(wlf),length(wle)*length(wlf)}, 1}(undef, CanopyRTMod.canopy.nlayers)
for i = 1:CanopyRTMod.canopy.nlayers
    #@show i
    arrayOfLeaves[i] = leafbio{FT, length(wl), length(wle), length(wlf),length(wle)*length(wlf)}()
    CanopyRTMod.fluspect!(arrayOfLeaves[i], CanopyRTMod.optis)
end

CanopyRTMod.angles.tto=-88
CanopyRTMod.angles.psi=50
@time CanopyRTMod.computeCanopyGeomProps!(canopy, angles,canOpt);
@time CanopyRTMod.computeCanopyMatrices!(arrayOfLeaves,canOpt);
@time CanopyRTMod.RTM_SW!(canopy, canOpt, canRad,sunRad, CanopyRTMod.soil);
CanOpt2 = deepcopy(canOpt)
@time CanopyRTMod.computeSIF_Fluxes!(arrayOfLeaves, canOpt, canRad, canopy, CanopyRTMod.soil);

# Just a test
@time CanopyRTMod.compCanopyOptsExact!(arrayOfLeaves,canOpt, canopy.lidf);

plot(CanOpt2.w[:,10])
plot!(CanOpt2.vb[:,1])
plot!(canOpt.w[:,1]/36)
plot!(canOpt.vb[:,1]/36)

CanopyRTMod.soil.albedo_SW[:] .=0.2;

reflVIS = []
reflNIR = []
SIF_FR = []
SIF_FR1 = []
SIF_FR2 = []
SIF_FR3 = []
SIF_FR4 = []
SIF_R1 = []
SIF_R2 = []
SIF_R3 = []
SIF_R4 = []
SIF_R  = []
Pso = []
Po = []
Ps = []
ko = []
CanopyRTMod.angles.tts=30
CanopyRTMod.angles.psi=0
CanopyRTMod.canopy.LAI = 3
for VZA=-89.5:0.5:89.5
    CanopyRTMod.angles.tto=VZA
    CanopyRTMod.computeCanopyGeomProps!(canopy, angles,canOpt);
    CanopyRTMod.computeCanopyMatrices!(arrayOfLeaves,canOpt);
    CanopyRTMod.RTM_SW!(canopy, canOpt, canRad,sunRad, CanopyRTMod.soil);
    CanopyRTMod.computeSIF_Fluxes!(arrayOfLeaves, canOpt, canRad, canopy, CanopyRTMod.soil);
    push!(reflVIS, canRad.alb_obs[28])
    push!(reflNIR, canRad.alb_obs[52])
    push!(SIF_R , canRad.SIF_obs[8])
    push!(SIF_FR, canRad.SIF_obs[20])
    push!(SIF_FR1, canRad.SIF_obs_sunlit[20])
    push!(SIF_FR2, canRad.SIF_obs_scattered[20])
    push!(SIF_FR3, canRad.SIF_obs_soil[20])
    push!(SIF_FR4, canRad.SIF_obs_shaded[20])
    push!(SIF_R1, canRad.SIF_obs_sunlit[8])
    push!(SIF_R2, canRad.SIF_obs_scattered[8])
    push!(SIF_R3, canRad.SIF_obs_soil[8])
    push!(SIF_R4, canRad.SIF_obs_shaded[8])
    push!(Pso, canOpt.Pso[1])
    push!(Ps, canOpt.Ps[1])
    push!(Po, canOpt.Po[1])
    push!(ko, canOpt.ko[1])
end

CanopyRTMod.angles.tto=89
CanopyRTMod.computeCanopyGeomProps!(canopy, angles,canOpt);
CanopyRTMod.computeCanopyMatrices!(arrayOfLeaves,canOpt);
CanopyRTMod.RTM_SW!(canopy, canOpt, canRad,sunRad, CanopyRTMod.soil);
CanopyRTMod.computeSIF_Fluxes!(arrayOfLeaves, canOpt, canRad, canopy, CanopyRTMod.soil);
plot(canRad.SIF_hemi)
plot!(canRad.SIF_sum)

VZA=-89.5:0.5:89.5
plot(VZA, ko, label="ko")

@show wl[28]
@show wl[52]
CanopyRTMod.angles.tts = 48

plot(VZA, Pso, label="Pso")
plot!(VZA, Po, label="Po")
plot!(VZA, Ps, label="Ps")
#@show wlf

plot(VZA, reflVIS, label="Red Reflectance")
plot!(VZA, SIF_FR/100, label="Red SIF")

plot(VZA, reflNIR, label="NIR Reflectance")
plot!(VZA, SIF_FR, label="Far Red SIF")

plot(VZA, SIF_FR1, label="FR SIF Sunlit")
plot!(VZA, SIF_FR2, label="FR SIF scattered")
plot!(VZA, SIF_FR3, label="FR SIF soil")
plot!(VZA, SIF_FR4, label="FR SIF shaded")
plot!(VZA, SIF_FR, label="FR SIF total")
ylims!(0,5)

plot(VZA, SIF_R1, label="R SIF Sunlit")
plot!(VZA, SIF_R2, label="R SIF scattered")
plot!(VZA, SIF_R3, label="R SIF soil")
plot!(VZA, SIF_R4, label="R SIF shaded")
plot!(VZA, SIF_R, label="R SIF ")
ylims!(0,2)

plot(VZA, SIF_FR./reflNIR./mean(SIF_FR./reflNIR), label="FR/refl")
plot!(VZA, SIF_R./reflVIS./mean(SIF_R./reflVIS), label="R/refl")
#plot!(VZA, , label="Far Red SIF")

CanopyRTMod.soil.albedo_SW[:] .=0.2;
reflRed_SZA = []
reflNIR_SZA = []
CanopyRTMod.canopy.Ω = 1.0
CanopyRTMod.angles.tto=0.2
CanopyRTMod.canopy.LAI=2.
for SZA=0.0:1:75
    CanopyRTMod.angles.tts=SZA
    CanopyRTMod.computeCanopyGeomProps!(canopy, angles,canOpt);
    CanopyRTMod.computeCanopyMatrices!(arrayOfLeaves,canOpt);
    CanopyRTMod.RTM_SW!(canopy, canOpt, canRad,sunRad, CanopyRTMod.soil);
    push!(reflRed_SZA, canRad.alb_direct[28])
    push!(reflNIR_SZA, canRad.alb_direct[52])
end

SZA=0:1:75
plot(SZA, reflNIR_SZA)

SZA=0:1:75
plot(SZA, reflRed_SZA)

# Test plots from Christiaan's papers

reflVIS = Float32[]
reflNIR = Float32[]
SIF_FR = Float32[]
SIF_R  = Float32[]
CanopyRTMod.angles.tts=48
CanopyRTMod.angles.psi=0
CanopyRTMod.canopy.LAI=3.22
for psi=0:360
    CanopyRTMod.angles.psi=psi
    for VZA=0:1:85
        CanopyRTMod.angles.tto=VZA

        CanopyRTMod.computeCanopyGeomProps!(canopy, angles,canOpt);
        CanopyRTMod.computeCanopyMatrices!(arrayOfLeaves,canOpt);
        CanopyRTMod.RTM_SW!(canopy, canOpt, canRad,sunRad, CanopyRTMod.soil);
        CanopyRTMod.computeSIF_Fluxes!(arrayOfLeaves, canOpt, canRad, canopy, CanopyRTMod.soil);
        push!(reflVIS, canRad.alb_obs[28])
        push!(reflNIR, canRad.alb_obs[52])
        push!(SIF_R , canRad.SIF_obs[8])
        push!(SIF_FR, canRad.SIF_obs[20])
    end
end

A = reshape(reflNIR, ( 86,361));
B = reshape(reflVIS, ( 86,361));
SIFFER = reshape(SIF_R, ( 86,361));
SIFFER_FR = reshape(SIF_FR, ( 86,361));

using Plots; pyplot()

#heatmap(A, cmap=)
hm = contourf(deg2rad.(collect((0:360))),collect(0:1:85),  A,  proj=:polar, color=:viridis, alpha=0.5)

#heatmap(A, cmap=)
hm = contourf(deg2rad.(collect((0:360))),collect(0:1:85),  B,  proj=:polar, color=:viridis)

hm = contourf(deg2rad.(collect((0:360))),collect(0:1:85),  SIFFER, proj=:polar, color=:viridis)

hm = contourf(deg2rad.(collect((0:360))),collect(0:1:85),  SIFFER_FR, proj=:polar, color=:viridis)

plot(CanopyRTMod.litab, cumsum(canopy.lidf))

litab_bnd  = FT[[0.,10.,20.,30.,40.,50.,60.,70.,80.,82.,84.,86.,88.] [10.,20.,30.,40.,50.,60.,70.,80.,82.,84.,86.,88., 90.] ];

cumsum(canopy.lidf)

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


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

using Markdown
using InteractiveUtils

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

# ╔═╡ c5eda0dc-215e-11eb-22d1-d76f1b9e047c
using Revise, Compose, Viznet, PlutoUI, TropicalTensors

# ╔═╡ 08f66f3c-2183-11eb-2260-cb73d013a3a9
begin
	# helper functions
	Compose.set_default_graphic_size(14cm, 14cm)
	function svgstring(c::Context)
		io = IOBuffer()
		SVG(io, Compose.default_graphic_width, Compose.default_graphic_height, true, :none)(c)
		String(take!(io))
	end
end;

# ╔═╡ 2824ef46-2165-11eb-23c9-a3335e6a6dba
md"## Define a lattice"

# ╔═╡ 47cbb81c-2169-11eb-0b07-9ff0d4a4c6b9
@bind lattice_type Select(["square", "chimera"]; default="square")

# ╔═╡ 9725485e-2161-11eb-2cf2-99d30b6a6dd6
md"""
``L_x`` = $(@bind Lx Slider(2:(lattice_type == "square" ? 20 : 6); default=(lattice_type == "square" ? 10 : 4), show_value=true))
"""

# ╔═╡ b320dc3a-2161-11eb-3edf-a1521d67bae0
md"""
``L_y`` = $(@bind Ly Slider(2:(lattice_type == "square" ? 20 : 6); default=(lattice_type == "square" ? 10 : 4), show_value=true))
"""

# ╔═╡ 126df144-215f-11eb-116e-a179ec34b675
lt = lattice_type == "square" ? SquareLattice(Lx, Ly) : ChimeraLattice(Lx, Ly)

# ╔═╡ 93f2d69e-21f5-11eb-13da-df4f1cc701fe
DownloadButton(svgstring(showlattice(lt)), "lattice_$(lattice_type)_Lx$(Lx)_Ly$Ly.svg")

# ╔═╡ 370cda48-2165-11eb-1514-d31e6fb4eb55
md"## Define the coupling"

# ╔═╡ b7cdbb3a-2162-11eb-0f41-6d89bcb699a6
md"""
coupling type = $(@bind coupling_type Select(["±1", "normal"]; default="±1"))
"""

# ╔═╡ 2cd5c642-2164-11eb-3e69-fd9f6da760fb
Js = if coupling_type == "±1"
	rand([-1.0,1], length(sgbonds(lt)))
elseif coupling_type == "normal"
	randn(length(sgbonds(lt)))
end;

# ╔═╡ 64f4a1ca-2163-11eb-20e8-cb6cfef00ec5
md"""
magnetic field = $(@bind magnetic_field Select(["abscent", "±1", "normal"]; default="abscent"))
"""

# ╔═╡ 1a3456ac-215f-11eb-3c67-39ef0ab19296
hs = if magnetic_field == "±1"
	rand([-1.0,1], length(sgvertices(lt)))
elseif magnetic_field == "normal"
	randn(length(sgvertices(lt)))
elseif magnetic_field == "abscent"
	zeros(Float64, length(sgvertices(lt)));
end;

# ╔═╡ 46644ab6-215f-11eb-08f5-378294e5fe24
# yellow -> negative coupling
# blue -> positive coupling
sg = Spinglass(lt, Js, hs)

# ╔═╡ b1036188-21f8-11eb-39cc-a5c83f86483b
DownloadButton(svgstring(viz(sg)), "spinglass_$(lattice_type)_Lx$(Lx)_Ly$Ly.svg")

# ╔═╡ 4c38a7f6-2165-11eb-15fd-7559fe9bbe86
md"## Find the optimal configuration"

# ╔═╡ 246f7f8c-2170-11eb-2659-a9ae4bc4022f
res = solve(CountingTropical{Float64}, sg; usecuda=false)

# ╔═╡ 6ff289d6-2166-11eb-3a5c-4f1e686b2810
md"The minimum energy is $(res.n), degenercy is $(res.c)"

# ╔═╡ 3f5c4dba-2162-11eb-14b9-632cdafa348d
optimal_config = Reversible.opt_config(sg)

# ╔═╡ e0469c76-21f8-11eb-39fa-d399a9f6654d
DownloadButton(svgstring(viz(optimal_config)), "optimalconfig_$(lattice_type)_Lx$(Lx)_Ly$Ly.svg")

# ╔═╡ Cell order:
# ╠═c5eda0dc-215e-11eb-22d1-d76f1b9e047c
# ╠═08f66f3c-2183-11eb-2260-cb73d013a3a9
# ╟─2824ef46-2165-11eb-23c9-a3335e6a6dba
# ╟─47cbb81c-2169-11eb-0b07-9ff0d4a4c6b9
# ╟─9725485e-2161-11eb-2cf2-99d30b6a6dd6
# ╟─b320dc3a-2161-11eb-3edf-a1521d67bae0
# ╠═126df144-215f-11eb-116e-a179ec34b675
# ╟─93f2d69e-21f5-11eb-13da-df4f1cc701fe
# ╟─370cda48-2165-11eb-1514-d31e6fb4eb55
# ╟─b7cdbb3a-2162-11eb-0f41-6d89bcb699a6
# ╠═2cd5c642-2164-11eb-3e69-fd9f6da760fb
# ╟─64f4a1ca-2163-11eb-20e8-cb6cfef00ec5
# ╠═1a3456ac-215f-11eb-3c67-39ef0ab19296
# ╠═46644ab6-215f-11eb-08f5-378294e5fe24
# ╟─b1036188-21f8-11eb-39cc-a5c83f86483b
# ╟─4c38a7f6-2165-11eb-15fd-7559fe9bbe86
# ╠═246f7f8c-2170-11eb-2659-a9ae4bc4022f
# ╟─6ff289d6-2166-11eb-3a5c-4f1e686b2810
# ╠═3f5c4dba-2162-11eb-14b9-632cdafa348d
# ╟─e0469c76-21f8-11eb-39fa-d399a9f6654d