text:
@testset "High-level wrapper" begin
    @testset "Generated structs" begin
        test(StructDefinition{true}, struct_by_name, :VkApplicationInfo, :(
            @auto_hash_equals struct ApplicationInfo <: HighLevelStruct
                next::Any
                application_name::String
                application_version::VersionNumber
                engine_name::String
                engine_version::VersionNumber
                api_version::VersionNumber
            end
        ))

        test(StructDefinition{true}, struct_by_name, :VkInstanceCreateInfo, :(
            @auto_hash_equals struct InstanceCreateInfo <: HighLevelStruct
                next::Any
                flags::UInt32
                application_info::OptionalPtr{ApplicationInfo}
                enabled_layer_names::Vector{String}
                enabled_extension_names::Vector{String}
            end
        ))

        test(StructDefinition{true}, struct_by_name, :VkXcbSurfaceCreateInfoKHR, :(
            @auto_hash_equals struct XcbSurfaceCreateInfoKHR <: HighLevelStruct
                next::Any
                flags::UInt32
                connection::Ptr{vk.xcb_connection_t}
                window::vk.xcb_window_t
            end
        ))
    end

    @testset "Conversion to low-level structs" begin
        test_ex(
            Constructor(
                StructDefinition{false}(struct_by_name(:VkApplicationInfo)),
                StructDefinition{true}(struct_by_name(:VkApplicationInfo)),
            ),
            :(
                _ApplicationInfo(x::ApplicationInfo) = _ApplicationInfo(x.application_version, x.engine_version, x.api_version; x.next, x.application_name, x.engine_name)
            ),
        )
        test_ex(
            Constructor(
                StructDefinition{false}(struct_by_name(:VkInstanceCreateInfo)),
                StructDefinition{true}(struct_by_name(:VkInstanceCreateInfo)),
            ),
            :(
                _InstanceCreateInfo(x::InstanceCreateInfo) = _InstanceCreateInfo(x.enabled_layer_names, x.enabled_extension_names; x.next, x.flags, application_info = convert_nonnull(_ApplicationInfo, x.application_info))
            ),
        )
    end

    @testset "Additional constructor" begin
        test(Constructor, StructDefinition{true}(StructDefinition{false}(struct_by_name(:VkApplicationInfo))), :(
            ApplicationInfo(application_version::VersionNumber, engine_version::VersionNumber, api_version::VersionNumber; next = C_NULL, application_name = "", engine_name = "") = ApplicationInfo(next, application_name, application_version, engine_name, engine_version, api_version)
        ))
    end
end


text:
using Statistics, OrderedCollections, Dates, Parameters

abstract type SData end

"""
InputData <: SData
Contains time series data (i.e., respecting a temporal order)
Attributes:
- region::String: optional information to specify the region data belongs to
- period_from::DateTime: date(year,month,day,hour) from which the data is from
- period_to::DateTime: date(year,month,day,hour) until which the data is
- resolution::Period: resolution of the original series (e.g., hours, days, weeks, etc.)
- nseries::Int: number of time series inputed (treated separately as if independent from each other)
- data::Dict{String,Array}: Dictionary with an entry for each attribute `[file name (attribute: e.g technology)]-[column name (node: e.g. location)]`, Each entry of the dictionary is a 2-dimensional `time-steps T x periods K`-Array holding the data
"""
mutable struct InputData
    region::String
    period_from::DateTime
    period_to::DateTime
    resolution::Period
    nseries::Int
    data::Dict{String,Array}
end

"""
SeriesInstance <: SData
Contains the series and respective attributes needed for performing the aggregation
Attributes:
- series::VecOrMat{T}: original time series
- block_size::Int: number of clusters to be aggregated at once each iteration
- stopping_k::Int: reduced number of clusters pre-established
- current_k::Int: current number of clusters
- dm::Symbol: discrepancy metric to be used to measure the difference between clusters to be aggregated
- rep_value::Symbol: function to assign a representative value for each cluster (by default, the mean of values inside the cluster)
- lseries::Int: length of series
- nseries::Int: number of series
"""
mutable struct SeriesInstance{T<:Float64,S<:Int}
    series::VecOrMat{T}
    block_size::S
    stopping_k::S
    current_k::S
    dm::Symbol
    rep_value::Symbol
    lseries::S
    nseries::S
end

"""
ClustInstance <: SData
Contains clusters information
Fields:
- k_cent::VecOrMat{T}:
- weights::Vector{<:Int}: 1-dimensional `periods K`-Array with the absolute weight for each period. E.g. for a yearly resolution time series of 365 days, sum(weights) = 365
- series_clust::Vector{<:Int}: clusters assignment for the original series (vector with the from-to relationship between the original series and the clusters)
- nclusters::Int: current number of clusters
- search_range::UnitRange: range in which the clusters can be merged, considering the amount of clusters and the block_size (generally equals to 1:(nclusters-block_size+1))
"""
mutable struct ClustInstance{T<:Float64,S<:Int}
    k_cent::VecOrMat{T}
    weights::Vector{S}
    series_clust::Vector{S}
    nclusters::S
    search_range::UnitRange
end

"""
AggregInstance <: SData
Contains the aggregation instance temporary values.
Attributes:
- merging_clust::UnitRange{<:Int}:
- series_comp::VecOrMat{T}:
- k_cent_comp::VecOrMat{T}:
"""
mutable struct AggregInstance{T<:Float64,S<:Int}
    merging_clust::UnitRange{S}
    series_comp::VecOrMat{T}
    k_cent_comp::VecOrMat{T}
end

"""
DistUpdate <: SData
Type to store the history of minimal distances.
"""

mutable struct DistUpdate{S<:Int}
    min_dist::S
    merging_clust::UnitRange{S}
end

"""
aggreg1D(series::VecOrMat{T}, marker::Vector{Bool}, rep_value::Symbol = :mean)
The function aggregates N one-dimensional vectors concomitantly following a marker gave as the input of the function. It assumes the collection of 1D vectors is vertically concatanated so forming a 2D matrix with the first dimension equal to the size of the vectors analysed and the second dimension equal to the number of vectors analysed.
Fields:
- series::VecOrMat{T}: collection of vectors
- marker::Vector{Bool}: vector with the flag on where to aggregate - so whenever marker[i] = 1 the aggregation happens between i and i+1
- rep_value::Symbol: value to be used for representation of the aggregated interval (we use mean by default)
"""
# (first method)
function aggreg1D(series::VecOrMat{T}, marker::Vector{Bool}, representation::Symbol = :mean) where {T<:Float64}
    #### Function to aggregate vectors under a representative value, normally used in hierarchical
    #### clustering. It merges a set of vectors (series) merging the state (i) with (i + 1) whenever there's a marker

    @assert length(size(series)) <= 2 "Series to be aggregated need to be represented as a Matrix (2D)"

    ## Initialization
    lseries = size(series,1)                      # Length of vectors
    nseries = size(series,2)                      # Number of vectors
    D = 1:nseries                                 # Range for the number of series
    class = collect(1:lseries)                    # Vector to mark the new classes
    new_v = copy(series)                          # New vector to be formed after aggregation
    aux = zeros(lseries)

    ## Error control
    @assert lseries == length(marker) "Different sizes"
    @assert any((marker .!= 0) .| (marker .!= 1)) "Inconsistent marker vector"

    ##Accounts the cumulative change in clusters order
    c_change = 0                        # Change in the number of clusters
    for i in 2:lseries
        c_change -= marker[i-1]
        aux[i] = c_change
    end

    ## Update values
    class .+= aux

    ## Calculate the representative value and replace it in the new_v
    if representation == :mean
        for i in 1:class[end], d in D
            new_v[class.==i,d] .= mean(series[class.==i,d])
        end
    elseif representation == :medoid
        for i in 1:class[end], d in D
            dist = pairwise(Euclidean(),series[class.==i,d])
            s_medoid = kmedoids(dist, 1)
            s_medoid = s_medoid.medoids[1]
            new_v[class.==i,d] .= series[class.==i,d][s_medoid]
        end
    else
        @assert false "Representation method not defined."
    end

    return new_v, class
end

"""
aggreg1D(series::VecOrMat{T}, rep_value::Symbol = :mean)
The function aggregates N one-dimensional vectors concomitantly following a marker gave as the input of the function. It assumes the collection of 1D vectors is vertically concatanated so forming a 2D matrix with the first dimension equal to the size of the vectors analysed and the second dimension equal to the number of vectors analysed.
Fields:
- series::VecOrMat{T}: collection of vectors
- rep_value::Symbol: value to be used for representation of the aggregated interval (we use mean by default)
"""
# (second method)
function aggreg1D(series::VecOrMat{T}, representation::Symbol = :mean) where {T<:Float64}
    #### Function to aggregate vectors under a representative value, normally used in hierarchical
    #### clustering. It merges a set of vectors (series) merging the state (i) with (i + 1) whenever there's a marker

    @assert length(size(series)) <= 2 "Series to be aggregated need to be represented as a Matrix (2D)"

    # Initialization
    lseries = size(series,1)                      # Length of vectors
    nseries = size(series,2)                      # Number of vectors
    D = 1:nseries                                 # Range for the number of series
    class = collect(1:lseries)                    # Vector to mark the new classes
    new_v = copy(series)                          # New vector to be formed after aggregation

    # Calculate the representative value and replace it in the new_v
    if representation == :mean
        for d in D
            new_v[:,d] .= mean(series[:,d])
        end
    elseif representation == :medoid
        for d in D
            dist = pairwise(Euclidean(),series[:,d])
            s_medoid = kmedoids(dist, 1)
            s_medoid = s_medoid.medoids[1]
            new_v[:,d] .= series[:,d][s_medoid]
        end
    else
        @assert false "Representation method not defined."
    end

    return new_v, class
end

"""
aggreg1D(series::VecOrMat{T}, marker::Vector{Bool}, rep_value::Symbol = :mean)
The function aggregates N one-dimensional vectors concomitantly following a marker gave as the input of the function. It assumes the collection of 1D vectors is vertically concatanated so forming a 2D matrix with the first dimension equal to the size of the vectors analysed and the second dimension equal to the number of vectors analysed.
Fields:
- series::VecOrMat{T}: collection of vectors
- marker::Vector{Bool}: vector indicating the assignment of the series
- rep_value::Symbol: value to be used for representation of the aggregated interval (we use mean by default)
"""
# (third method)
function aggreg1D(series::VecOrMat{T}, marker::Vector{Int}, representation::Symbol = :mean) where {T<:Float64}
    #### Function to aggregate vectors under a representative value, normally used in hierarchical
    #### clustering. It merges a set of vectors (series) merging the state (i) with (i + 1) whenever there's a marker

    @assert length(size(series)) <= 2 "Series to be aggregated need to be represented as a Matrix (2D)"

    # Initialization
    lseries = size(series,1)                      # Length of vectors
    nseries = size(series,2)                      # Number of vectors
    D = 1:nseries                                 # Range for the number of series
    class = collect(1:lseries)                    # Vector to mark the new classes
    new_lseries = maximum(marker)                 # New vectors length

    # Error control
    @assert lseries == length(marker) "Different sizes between series and assignment vector"

    # Finding the unique class vector
    class = unique(sort(marker))
    weights = [count(i->(i == j),marker) for j in class] .|> Int

    # Declaring the new VecOrMat object
    new_v = zeros(length(class))                          # New vector to be formed after aggregation

    # Calculate the representative value and replace it in the new_v
    if representation == :mean
        for i in 1:length(class), d in D
            new_v[i,d] .= mean(series[marker .== class[i],d])
        end
    elseif representation == :medoid
        for i in 1:class[end], d in D
            dist = pairwise(Euclidean(),series[marker .== class[i],d])
            s_medoid = kmedoids(dist, 1)
            s_medoid = s_medoid.medoids[1]
            new_v[i,d] .= series[marker .== class[i],d][s_medoid]
        end
    else
        @assert false "Representation method not defined."
    end

    return new_v, class, weights
end

"""
cdad(u::Vector{<:Float64},v::Vector{<:Float64})
Function to compute the cumulative distribution absolute difference from vector `u`
to vector `v` assuming both have an equal mass (i.e., both normalised) and
they are ordered.
"""
function cdad(u::Vector{T},v::Vector{T}) where {T<:Float64}
    # Error control
    @assert length(u) == length(v) "Vectors need to have an equal length."
    @assert sum(u) ≈ sum(v) "Vectors need to have an equal mass."

    # Compute the differences between pairs of successive values from u and v.
    deltas = u .- v

    # Calculate the CDFs of u and v using their weights, if specified.
    cdf = sum(abs.(cumsum(deltas)))

    return cdf
end

"""
load_series_instance(series::VecOrMat{T}, block_size::S, current_k::S, stopping_k::S = 1, dm::Symbol = :ward, rep_value::Symbol = :mean, lseries::S = size(series,1), nseries::S = size(series,2)
) where {T<:Float64, S<:Int}
Function to load the data provided and fit into a struct SeriesInstance().
"""
function load_series_instance(series::VecOrMat{T},
    block_size::S,
    current_k::S,
    stopping_k::S = 1,
    dm::Symbol = :ward,
    rep_value::Symbol = :mean,
    lseries::S = size(series,1),
    nseries::S = size(series,2)
) where {T<:Float64, S<:Int}

    @assert block_size <= current_k - stopping_k + 1 "Aggregation not possible: block_size $block_size and stopping_k $stopping_k need to be checked."

    # Create an object of type SeriesInstance
    _SeriesInstance = SeriesInstance(series,block_size,stopping_k,current_k,dm,rep_value,lseries,nseries)

    return _SeriesInstance
end

"""
load_clust_instance(k_cent::VecOrMat{T}, weights::Vector{S} = ones(size(k_cent,1)), series_clust::Vector{S}, search_range::UnitRange = collect(1:size(k_cent,1))
) where {T<:Float64,S<:Int}
Function to load the clusters attributes.
"""
function load_clust_instance(k_cent::VecOrMat{T},
    series_clust::Vector{S},
    weights::Vector{S} = ones(size(k_cent,1)),
    search_range::UnitRange = collect(1:size(k_cent,1))
) where {T<:Float64,S<:Int}

    nclusters = size(k_cent,1)

    # Create an object of type SeriesInstance
    _ClustInstance = ClustInstance(k_cent,weights,series_clust,nclusters,search_range)

    return _ClustInstance
end

"""
search_min_dist(_SeriesInstance, _ClustInstance)
Search minimal distance and return the marker, the min_dist position, and the clusters to be merged.
"""
function search_min_dist(_SeriesInstance, _ClustInstance)
    # Unpacking SeriesInstance
    series = _SeriesInstance.series
    block_size = _SeriesInstance.block_size
    dm = _SeriesInstance.dm
    rep_value = _SeriesInstance.rep_value
    nseries = _SeriesInstance.nseries

    # Unpacking ClustInstance
    series_clust = _ClustInstance.series_clust
    search_range = _ClustInstance.search_range

    # Vector with distances (to be updated as it goes)
    dist = Vector{Float64}(undef,length(search_range))

    # Sets
    N = 1:nseries
    K = copy(search_range)

    # Compute the distance for each aggregation (i.e., changing the merging_clust)
    @inbounds for k in K
        # Merging to be tested
        merging_clust = k:k+block_size-1
        # Create a temporary marker to merge the clusters tested
        marker_temp = [sc in merging_clust for sc in series_clust]
        # Part of series compared
        series_comp = series[marker_temp,:]
        # Centroids of the temporarily formed cluster (TODO: implement another method for aggreg1D receiving the clusters with respective weights)
        (k_cent_comp,) = aggreg1D(series_comp, rep_value)

        # Distance computation
        dist[k] = compute_dist(N, dm, series_comp, k_cent_comp)
    end

    # Find whenever the min_dist occurs first (i.e., using findmin()[2])
    ## TODO: implement the multiple merges (e.g., using findall())
    min_dist = findmin(dist)[2] |> Int
    merging_clust = min_dist:min_dist+block_size-1

    # Create a flag to the positions in the series that will be aggregated
    marker = update_marker(_SeriesInstance, _ClustInstance, min_dist)

    return marker, min_dist, merging_clust
end

"""
find_clusters!(_SeriesInstance, _ClustInstance, _DistUpdate)
Function calling the necessary functions for the data aggregation procedure (to be enveloped in a loop with the necessary conditions, e.g., while nclusters > stopping_k).
"""
function find_clusters!(_SeriesInstance, _ClustInstance, _DistUpdate)

    # Find the min distance and create a marker
    (marker, min_dist, merging_clust) = search_min_dist(_SeriesInstance, _ClustInstance)

    # Update _DistUpdate dictionary with the minimal distance found and the new marker
    _DistUpdate = merge(+, _DistUpdate, Dict(marker => DistUpdate(min_dist,merging_clust)))

    # Update clusters and series_clust
    update_clust!(_ClustInstance, _SeriesInstance, min_dist)

    # Update number of clusters k
    new_current_k = _ClustInstance.nclusters
    update_k!(_SeriesInstance, new_current_k)

    return new_current_k, _DistUpdate
end

"""
function compute_dist(N::UnitRange, dm::Symbol, series_comp::VecOrMat{T}, k_cent_comp::VecOrMat{T})  where {T <: Float64}
Compute the distances between two matrices (series_comp and k_cent_comp).
"""
function compute_dist(N::UnitRange, dm::Symbol, series_comp::VecOrMat{T}, k_cent_comp::VecOrMat{T})  where {T <: Float64}
    # Distances matrix designation (TODO: add more discrepancy metrics (e.g., DTW and others) and wrap their calculation in a function)
    if dm == :wd  # Wasserstein distance (or cumulative distribution absolute difference) between the agglomerated cluster and the original series
        dist = sum(cdad(series_comp[:,n], k_cent_comp[:,n]) for n in N)
    elseif dm == :ward  # Ward's cluster criterion (min variance proxy) for the original series when clustering `series_range`
        dist = sum(sum((series_comp[:,n] .- k_cent_comp[:,n]).^2 for n in N))
    else
        @assert false "Distance metric not defined."
    end

    return dist
end

"""
update_marker(_SeriesInstance, _ClustInstance, min_dist::T) where {T <: Int}
Update the clusters assignment (Bool marker indicating which positions are to be merged)
"""
## TODO: implement the case of simultaneous multiple merges
function update_marker(_SeriesInstance, _ClustInstance, min_dist::T) where {T <: Int}
    block_size = _SeriesInstance.block_size
    lseries = _SeriesInstance.lseries
    series_clust = _ClustInstance.series_clust

    # Create the new marker with 1 for the elements within the series to be aggregated
    marker = [series_clust[i] in min_dist:min_dist+block_size-1 for i in 1:lseries] |> Vector{Bool}

    return marker
end

## TODO: normalisation function (for example to use medoids representative values and the wasserstein distance as the discrepancy metric - via cdad)

"""
function replace_lines(mat::VecOrMat{T}, mergings::OrderedDict{UnitRange{Int64},Matrix{T}}) where {T <: Float64}
Replace lines of a matrix/vector by another matrix/vector with the same number of columns than the original matrix/vector.
"""
function replace_lines(mat::VecOrMat{T}, mergings::OrderedDict{UnitRange{Int64},Matrix{T}}) where {T <: Float64}
    nmergings = length(mergings)
    nrows = size(mat,1)
    ncols = size(mat,2)
    merg_keys = keys(mergings) .|> UnitRange

    # Error control for overlapping intervals and matrix scope
    lmerged = [] |> Vector{Int}
    for i in 1:nmergings
        lmerged = vcat(lmerged,collect(merg_keys[i]))
    end
    lmerged = sort(lmerged)
    @assert maximum(count(i -> i == lmerged[j], lmerged) for j in 1:length(lmerged)) == 1 "Overlapping merging intervals"
    @assert all([i in 1:nrows for i in lmerged]) "Lines to be merged out of matrix scope"

    # Flag to indicate if the merging was already transfered for the final matrix
    flag_merg = zeros(nmergings) |> Vector{Bool}

    # New matrix
    new_mat = Array{Float64}(undef,0,ncols)

    # Loop to merge the new lines
    for i in 1:nrows
        merg_pos = findmax([i in merg_keys[j] for j in 1:nmergings])[2]       # find the position in `mergings` where this merging happens
        if (i in lmerged) && (flag_merg[merg_pos] == 0)                       # case in which the line is in the mergings list and still not added to the new matrix
            new_mat = vcat(new_mat,mergings[merg_keys[merg_pos]])
            flag_merg[merg_pos] = 1
        elseif !(i in lmerged)
            new_mat = vcat(new_mat,reshape(mat[i,:],(1,ncols)))
        end
    end

    return new_mat
end

"""
function update_clust!(_ClustInstance, _SeriesInstance, min_dist::T) where {T <: Int}
Updates clusters attributes in the _ClustInstance object.
"""
function update_clust!(_ClustInstance, _SeriesInstance, min_dist::T) where {T <: Int}
    series = _SeriesInstance.series
    series_clust = _ClustInstance.series_clust
    nclusters = _ClustInstance.nclusters
    k_cent = _ClustInstance.k_cent

    block_size = _SeriesInstance.block_size
    rep_value = _SeriesInstance.rep_value
    nseries = _SeriesInstance.nseries
    lseries = _SeriesInstance.lseries

    # Ranges
    L = 1:lseries
    merging_clust = min_dist:min_dist + block_size - 1

    @assert length(merging_clust) >= 2 "unitary merging_clust not allowed"

    # Updating the centroids
    (k_change,) = aggreg1D(series[[sc in merging_clust for sc in series_clust],:], rep_value)
    k_change = reshape(k_change[1,:],(1,nseries))
    k_cent = replace_lines(k_cent,sort(Dict(merging_clust => k_change)))

    # Update series_clust
    change_marker = zeros(lseries) |> Vector{Int}                        # Mark whenever the change in the cluster assignment happens for the series_clust
    @inbounds for l in 2:L[end]
        if (series_clust[l] in merging_clust[2:end]) && (series_clust[l] != series_clust[l-1])
            change_marker[l] = 1
        end
    end

    # Transform change_marker to carry the cumulative change in the clusters order
    c_change = 0                        # Change in the number of clusters
    @inbounds for l in L
        c_change = c_change - change_marker[l]
        change_marker[l] = c_change
    end

    series_clust = (series_clust .+ change_marker) |> Vector{Int64}
    series_clust = series_clust |> Vector{Int64}

    # New number of clusters
    nclusters = series_clust[end]
    # Update weights
    weights = [sum(series_clust .== i) for i in 1:nclusters] .|> Int
    # Update search_range
    new_search_range = 1:(1+nclusters-block_size)

    # Update struct
    _ClustInstance.series_clust = series_clust
    _ClustInstance.weights = weights
    _ClustInstance.nclusters = nclusters
    _ClustInstance.k_cent = k_cent
    _ClustInstance.search_range = new_search_range

    return
end

"""
update_k!(_SeriesInstance, new_current_k::Int)
Updates number of clusters in the _SeriesInstance object.
"""
function update_k!(_SeriesInstance, new_current_k::Int)
    _SeriesInstance.current_k = new_current_k
    stopping_k = _SeriesInstance.stopping_k

    @assert new_current_k >= stopping_k "Number of clusters $new_current_k is now < then $stopping_k"
end

ext: jl


text:
using GeometryTypes, ColorTypes
using FactCheck
import Base.Test.@inferred

facts("GeometryTypes") do
    include("polygons.jl")
    include("hyperrectangles.jl")
    include("faces.jl")
    include("meshes.jl")
    include("distancefields.jl")
    include("primitives.jl")
    include("decompose.jl")
    include("simplerectangle.jl")
    include("hypersphere.jl")
    include("typeutils.jl")
    include("simplices.jl")
    include("convexhulls.jl")
    include("gjk.jl")
    include("lines.jl")
end

FactCheck.exitstatus()


text:
# Shaded line plots

import .LibCImPlot: PlotShaded

function PlotShaded(x::AbstractArray{T,1}, y1::AbstractArray{T,1}, y2::AbstractArray{T,1};
                               label_id::String = "", count::Integer = min(length(x), length(y2),
                               length(y2)), offset::Integer = 0, stride::Integer = 1) where {T <: ImPlotData}

    LibCImPlot.PlotShaded(label_id, x, y1, y2, count, offset, stride * sizeof(T))
end

function PlotShaded(x::AbstractArray{T1}, y1::AbstractArray{T2}, y2::AbstractArray{T3};
                    kwargs...) where {T1<:Real, T2<:Real, T3<:Real}
    LibCImPlot.PlotShaded(promote(x,y1,y2)...; kwargs...)
end

function PlotShaded(x::AbstractArray{T}, y::AbstractArray{T}, y_ref::T;
                    count::Integer = min(length(x),length(y)), offset::Integer = 0,
                    stride::Integer = 1, label_id::String = "") where {T <: ImPlotData}

    LibCImPlot.PlotShaded(label_id, x, y, count, y_ref, offset, stride * sizeof(T))
end

function PlotShaded(x::AbstractArray{T1}, y::AbstractArray{T2}, y_ref::T3;
                    kwargs...) where {T1<:Real, T2<:Real, T3<:Real}
    (x, y) = promote(x, y)
    y_ref = eltype(x)(y_ref)
    PlotShaded(x, y, y_ref; kwargs...)
end

function PlotShaded(x::UnitRange{<:Integer}, y1::AbstractArray{T}, y2::AbstractArray{T}; label_id::String="") where {T<:ImPlotData}
    count = length(x)
    offset = x.start >= 1 ? x.start - 1 : throw("Offset out of bounds")
    x = collect(T, x)
    LibCImPlot.PlotShaded(label_id, x, y1, y2, count, offset, sizeof(T))
end

function PlotShaded(x::UnitRange{<:Integer}, y::AbstractArray{T}, y_ref::T; label_id::String="") where {T <: ImPlotData}
    count = length(x)
    offset = x.start >= 1 ? x.start - 1 : throw("Offset out of bounds")
    x = collect(T, x)
    LibCImPlot.PlotShaded(label_id, x, y, count, y_ref, offset, sizeof(T))
end

function PlotShaded(x::StepRange, y1::AbstractArray{T}, y2::AbstractArray{T}; label_id::String="") where {T <: ImPlotData}
    x.stop < 1 && throw("Range out of bounds")
    count = length(x)
    offset = x.start >= 1 ? x.start -1 : throw("Range out of bounds")
    stride = x.step >= 1 ? x.step : throw("Stride must be greater than zero.")
    x = collect(T, x)
    LibCImPlot.PlotShaded(label_id, x, y1, y2, count, offset, stride * sizeof(T))
end

function PlotShaded(x::StepRange, y::AbstractArray{T}, y_ref::T; label_id::String = "") where {T <: ImPlotData}
    x.stop < 1 && throw("Range out of bounds")
    count = length(x)
    offset = x.start >= 1 ? x.start -1 : throw("Range out of bounds")
    stride = x.step >= 1 ? x.step : throw("Stride must be greater than zero.")

    x = collect(T, x)
    LibCImPlot.PlotShaded(label_id, x, y, count, y_ref, offset, stride * sizeof(T))
end

function PlotShaded(x::OrdinalRange, y1::AbstractArray{T1}, y2::AbstractArray{T2};
                    kwargs...) where {T1<:Real, T2<:Real}
    LibCImPlot.PlotShaded(x,promote(y1,y2)...; kwargs...)
end

function PlotShaded(x::OrdinalRange, y::AbstractArray{T1}, y_ref::T2;
                    kwargs...) where {T1<:Real, T2<:Real}
    y_ref = eltype(y)(y_ref)
    LibCImPlot.PlotShaded(x, y, y_ref; kwargs...)
end


text:
"""
Add a variable to a Mathematical Programm. See JuMP `@variable` for examples
"""
macro variableMP(mp, args...)
    mp = esc(mp)
    # Pick out keyword arguments
    if Meta.isexpr(args[1], :parameters) # these come if using a semicolon
        kwargs = args[1]
        args = args[2:end]
    else
        kwargs = Expr(:parameters)
    end
    kwsymbol = VERSION < v"0.6.0-dev.1934" ? :kw : :(=) # changed by julia PR #19868
    append!(kwargs.args, collect(Iterators.filter(x -> Meta.isexpr(x, kwsymbol), collect(args)))) # comma separated
    args = collect(Iterators.filter(x -> !Meta.isexpr(x, kwsymbol), collect(args)))
    if length(kwargs.args) > 0
        error("Please do not use keyword argument in @variableMP")
    end

    # TODO(xhub) fix this with proper syntax
    # dummyconstr = Expr(:call, @variable, $(mp).emp.m, $(esc(args[1])))
    #then extend the arguments
    local len = length(args)
    if len == 0
        code = :( v = @eval @variable $(mp).emp.model_ds )
    elseif len == 1
        code = :( v = @variable $(mp).emp.model_ds $(esc(args[1])) )
    elseif len == 2
        code = :( v = @eval @variable $(mp).emp.model_ds $(esc(args[2])) $(esc(args[3])) )
    elseif len == 3
        code = :( v = @eval @variable $(mp).emp.model_ds $(esc(args[2])) $(esc(args[3])) $(esc(args[4])) )
    elseif len == 4
        code = :( v = @eval @variable $(mp).emp.model_ds $(esc(args[2])) $(esc(args[3])) $(esc(args[4])) $(esc(args[5])) )
    elseif len == 6
        code = :( v = @eval @variable $(mp).emp.model_ds $(esc(args[2])) $(esc(args[3])) $(esc(args[4])) $(esc(args[5])) $(esc(args[6])) )
    end
    quote
        if $len > 6
            error("unsupported syntax $(esc(args[2]))")
        end
        $code
        addvar!($(mp), v)
        v
    end
end

macro constraintFromExprMP(mp, expr)
    #dummyconstr = Expr(:call, esc(:(==)), esc(expr), 0)
    dummyconstr = Expr(:call, :(==), esc(expr), 0)
    code = :( cref = @constraint $(esc(mp)).emp.model_ds $dummyconstr )
    quote
        $code
        local gdix = addequ!($(esc(mp)), cref)
        (cref, gidx)
    end
end

macro constraintFromExprMP2(mp, constr)
    quote
        local cref = @constraint $(esc(mp)).emp.model_ds $(esc(constr))
        local gdix = addequ!($(esc(mp)), cref)
        (cref, gdix)
    end
end

macro NLconstraintFromExprMP(mp, expr)
    dummyconstr = Expr(:call, esc(:(==)), esc(expr), 0)
    #dummyconstr = Expr(:call, :(==), esc(expr), 0)
    code = :( cref = @NLconstraint $(esc(mp)).emp.model_ds $dummyconstr )
    quote
        $code
        local gidx = addequ!($(esc(mp)), cref)
        (cref, gidx)
    end
end
# TODO(xhub) add check warning for an already defined objective equation
"""
    @NLobjectiveMP(mp, sense, expr)

Add a nonlinear objective to a mathematical programm. Sense is either `:Min` or `:Max`.
"""
macro NLobjectiveMP(mp, sense, expr)
    dummyconstr = Expr(:call, esc(:(==)), esc(expr), 0)
    code = :( cref = @NLconstraint $(esc(mp)).emp.model_ds $dummyconstr )
    quote
        $code
        gidx = addequ!($(esc(mp)), cref)
        $(esc(mp)).objequ = gidx
        $(esc(mp)).sense = $(esc(sense))
        cref
    end
end

"""
    @objectiveMP(mp, sense, expr)

Add a linear objective to a mathematical programm. Sense is either `:Min` or `:Max`.
"""
macro objectiveMP(mp, sense, expr)
    dummyconstr = Expr(:call, esc(:(==)), esc(expr), 0)
    code = :( cref = @constraint $(esc(mp)).emp.model_ds $dummyconstr )
    quote
        $code
        gidx = addequ!($(esc(mp)), cref)
        $(esc(mp)).objequ = gidx
        $(esc(mp)).sense = $(esc(sense))
        cref
    end
end

"""
    @constraintMP(mp, expr)

Add a linear constraint to a mathematical programm
"""
macro constraintMP(mp, expr)
    quote
        cref = @constraint($(esc(mp)).emp.model_ds, $(esc(expr)))
        gdix = addequ!($(esc(mp)), cref)
        cref
    end
end

"""
    @constraintMP(mp, name, expr)

Add a linear constraint (with a identifier `name`) to a mathematical programm
"""
macro constraintMP(mp, name, expr)
    quote
        cref = @constraint($(esc(mp)).emp.model_ds, $(esc(name)), $(esc(expr)))
        gdix = addequ!($(esc(mp)), cref)
        cref
    end
end

"""
    @NLconstraintMP(mp, expr)

Add a nonlinear constraint to a mathematical programm
"""
macro NLconstraintMP(mp, expr)
    quote
        cref = @NLconstraint($(esc(mp)).emp.model_ds, $(esc(expr)))
        gdix = addequ!($(esc(mp)), cref)
        cref
    end
end

"""
    @NLconstraintMP(mp, name, expr)

Add a nonlinear constraint (with a identifier `name`) to a mathematical programm
"""
macro NLconstraintMP(mp, name, expr)
    quote
        cref = @NLconstraint($(esc(mp)).emp.model_ds, $(esc(name)), $(esc(expr)))
        gdix = addequ!($(esc(mp)), cref)
        cref
    end
end

"""
    vipair(mp, expr, var)

Add an affine variational inequality relationship between the variable `var` and the mapping `expr`

If there is no additional constraint of the feasible set of `var` besides lower and upper bounds,
then this is defines a Mixed Complementarity Problem
    expr ⟂ lb ≤ var ≤ ub
Or
    0 ∈ expr + N_[lb, ub] (var)
"""
function vipair(mp::MathPrgm, expr, var::JuMP.Variable)
#    cref, eidx = @constraintFromExprMP(mp, expr)
    cref, eidx = @constraintFromExprMP2(mp, expr == 0.)
    mp.matching[var.col] = eidx
    cref
end

function vipair(mp::MathPrgm, expr::Vector, var::Vector{JuMP.Variable})
    @assert length(expr) == length(var)
    for i in 1:length(var)
        vipair(mp, expr[i], var[i])
    end
end

macro vipair(mp, expr, var)
    quote
        return vipair($(esc(mp)), $(esc(expr)), $(esc(var)))
    end
end

"""
    @NLvipair(mp, expr, var)

Add a nonlinear variational inequality relationship between the variable `var` and the mapping `expr`

If there is no additional constraint of the feasible set of `var` besides lower and upper bounds,
then this is defines a Mixed Complementarity Problem
    expr ⟂ lb ≤ var ≤ ub
Or
    0 ∈ expr + N_[lb, ub] (var)
"""
macro NLvipair(mp, expr, var::JuMP.Variable)
    cref, eidx = @NLconstraintFromExprMP(mp, expr)
    mp.matching[var.col] = eidx
    cref
end

#function NLvipair(mp::MathPrgm, expr::Vector{Any}, var::Vector{JuMP.Variable})
#    @assert length(expr) == length(var)
#    for i in 1:length(var)
#        NLvipair(mp, expr[i], var[i])
#    end
#end


text:
using Test
using DomainIntegrals, DomainSets, BasisFunctions, GridArrays, FillArrays

using BasisFunctions.Test: generic_test_discrete_measure

@testset "Discrete measure" begin
    g = ScatteredGrid(rand(2))
    μ = discretemeasure(g)
    generic_test_discrete_measure(μ)
    @test weights(μ) isa FillArrays.Ones
    @test !isnormalized(μ)

    g2 = g×g
    μ = discretemeasure(g2)
    generic_test_discrete_measure(μ)
    @test Matrix(weights(μ))≈ones(2,2)
    @test !isnormalized(μ)

    μ = DiracWeight(rand())
    generic_test_discrete_measure(μ)
    @test weights(μ) isa FillArrays.Ones
    @test isnormalized(μ)

    μ = discretemeasure(PeriodicEquispacedGrid(3,0,1))
    @test isuniform(μ)
    generic_test_discrete_measure(μ)
    μ = discretemeasure(MidpointEquispacedGrid(3,0,1))
    @test isuniform(μ)
    generic_test_discrete_measure(μ)
    μ = discretemeasure(EquispacedGrid(3,0,1))
    @test isuniform(μ)
    generic_test_discrete_measure(μ)
    μ = BasisFunctions.NormalizedDiracComb(EquispacedGrid(3,0,1))
    @test isnormalized(μ)
    generic_test_discrete_measure(μ)
    μ = discretemeasure(MidpointEquispacedGrid(3,0,1),rand(3))
    generic_test_discrete_measure(μ)

    g = FourierGrid(10)
    mg = rescale(FourierGrid(10),-1, 1)
    m = forward_map(mg)
    μ = discretemeasure(mg)
    generic_test_discrete_measure(μ)
    @test forward_map(μ) == m
    @test supermeasure(μ) == discretemeasure(g)
    @test BasisFunctions.apply_map(discretemeasure(g),m) ≈ μ
    @test BasisFunctions.apply_map(μ,inverse(m))≈discretemeasure(g)

    g = PeriodicEquispacedGrid(3,-1,1)×PeriodicEquispacedGrid(4,-1,1)
    μ = discretemeasure(g)
    generic_test_discrete_measure(μ)
    @test μ isa BasisFunctions.DiscreteProductWeight
    @test components(μ) == (component(μ,1),component(μ,2))
end


text:
using LinearAlgebra, SparseArrays, SuiteSparse
mutable struct LinSolveFactorize{F}
  factorization::F
  A
end
LinSolveFactorize(factorization) = LinSolveFactorize(factorization,nothing)
function (p::LinSolveFactorize)(x,A,b,update_matrix=false)
  if update_matrix
    p.A = p.factorization(A)
  end
  if typeof(p.A) <: SuiteSparse.UMFPACK.UmfpackLU || typeof(p.factorization) <: typeof(lu)
    ldiv!(x,p.A,b) # No 2-arg form for SparseArrays!
  else
    x .= b
    ldiv!(p.A,x)
  end
end
DEFAULT_LINSOLVE = LinSolveFactorize(lu!)
function (p::LinSolveFactorize)(::Type{Val{:init}},f,u0_prototype)
  LinSolveFactorize(p.factorization,nothing)
end


text:
using Random, LinearAlgebra, TensorRefinement, Test

Random.seed!(1365)

@time @testset "Factorization orthogonalization" begin
	n = [2,1,3]; r = [5,4,3,7,8]
	d = length(n)
	L = length(r)-1
	@testset "Type: $S" for T ∈ (Float64, Float32), S ∈ (T, Complex{T})
		tol = 100*eps(T)
		U = decrand(n, r); Uf = factor(U); nrm = norm(Uf)
		@testset "pivot=$pivot, path=$path" for pivot ∈ (false, true), path ∈ ("forward", "backward")
			V = decqr!(deepcopy(U); pivot=pivot, path=path)
			Vf = factor(V)
			@test norm(Vf-Uf)/nrm ≈ 0 atol=tol
			(path == "backward") && decreverse!(V)
			p = decrank(V)
			for ℓ ∈ 1:L
				Q = reshape(V[ℓ], p[ℓ]*prod(n), p[ℓ+1])
				E = Q'*Q-I
				@test norm(E)/p[ℓ+1]^2 ≈ 0 atol=tol
			end
		end
    end
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

"""
    struct ForwardDiffAD

*Experimental feature, not part of stable public API.*

Constructors:

* ```ForwardDiffAD()```
"""
struct ForwardDiffAD <: DifferentiationAlgorithm end
export ForwardDiffAD

struct ForwardDiffGradient{
    F<:Function,
    VS<:AbstractValueShape,
    GS<:AbstractValueShape,
} <: GradientFunction
    _unshaped_f::F
    _input_shape::VS
    _grad_shape::GS
end

function ForwardDiffGradient(f::Function)
    input_shape = varshape(f)
    unshaped_f = unshaped(f)
    grad_shape = gradient_shape(input_shape)
    ForwardDiffGradient(unshaped_f, input_shape, grad_shape)
end

function _run_forwarddiff!(grad_f_x::AbstractVector{<:Real}, ::Type{T}, f::Function, x::AbstractVector{<:Real}) where {T<:Real}
    result = DiffResults.MutableDiffResult(zero(T), (grad_f_x,))

    # chunk = ForwardDiff.Chunk(x)
    # config = ForwardDiff.GradientConfig(f, x, chunk)
    # ForwardDiff.gradient!(result, f, x, config)
    ForwardDiff.gradient!(result, f, x)

    @assert DiffResults.gradient(result) === grad_f_x
    DiffResults.value(result)
end

function (gf::ForwardDiffGradient)(v::Any)
    input_shape = gf._input_shape
    v_shaped = fixup_variate(input_shape, v)
    v_unshaped = unshaped(v_shaped)
    R = density_logval_type(v_unshaped, default_dlt())

    grad_f_unshaped = similar(v_unshaped)

    value = R(_run_forwarddiff!(grad_f_unshaped, R, gf._unshaped_f, v_unshaped))

    (value, gf._grad_shape(grad_f_unshaped))
end

function (gf::ForwardDiffGradient)(::typeof(!), grad_f::Any, v::Any)
    input_shape = gf._input_shape
    v_shaped = fixup_variate(input_shape, v)
    v_unshaped = unshaped(v_shaped)
    R = density_logval_type(v_unshaped, default_dlt())

    result = if isnothing(grad_f)
        R(gf._unshaped_f(v_unshaped))
    else
        grad_f_unshaped = unshaped(grad_f, gf._grad_shape)
        R(_run_forwarddiff!(grad_f_unshaped, R, gf._unshaped_f, v_unshaped))
    end
end

function bat_valgrad_impl(f::Function, algorithm::ForwardDiffAD)
    (result = ForwardDiffGradient(f),)
end


text:
function GitHubActions()
    result = GitHubActions(
        "github-actions[bot]",
        "41898282+github-actions[bot]@users.noreply.github.com",
    )
    return result
end

function TeamCity()
    result = TeamCity(
        "github-actions[bot]",
        "41898282+github-actions[bot]@users.noreply.github.com",
    )
    return result
end

function get_my_username(ci_cfg::Union{GitHubActions, TeamCity};
                         auth = nothing,
                         env::AbstractDict = ENV)
    return ci_cfg.username
end

function set_git_identity(ci_cfg::Union{GitHubActions, TeamCity})
    run(Cmd(String["git", "config", "user.name", strip(ci_cfg.username)]))
    run(Cmd(String["git", "config", "user.email", strip(ci_cfg.email)]))
    return nothing
end

function auto_detect_ci_service(; env::AbstractDict = ENV)
    if haskey(env, "GITHUB_REPOSITORY")
        return GitHubActions()
    else
         error("Could not detect CI service")
    end
end

function github_repository(ci_cfg::Union{GitHubActions, TeamCity};
                           env::AbstractDict = ENV)
    result = env["GITHUB_REPOSITORY"]::String
    return result
end

function github_token(ci_cfg::Union{GitHubActions, TeamCity};
                      env::AbstractDict = ENV)
    result = env["GITHUB_TOKEN"]::String
    return result
end


text:
# ---
# title: 106. Construct Binary Tree from Inorder and Postorder Traversal
# id: problem106
# author: zhwang
# date: 2022-03-02
# difficulty: Medium
# categories: Array, Tree, Depth-first Search
# link: <https://leetcode.com/problems/construct-binary-tree-from-inorder-and-postorder-traversal/description/>
# hidden: true
# ---
#
# Given inorder and postorder traversal of a tree, construct the binary tree.
#
# **Note:**
# You may assume that duplicates do not exist in the tree.
#
# For example, given
#
#
#
#     inorder = [9,3,15,20,7]
#     postorder = [9,15,7,20,3]
#
# Return the following binary tree:
#
#
#
#         3
#        / \
#       9  20
#         /  \
#        15   7
#
#
#
## @lc code=start
using LeetCode

function build_tree_inpost(inorder::AbstractArray, postorder::AbstractArray)::TreeNode
    root = TreeNode(last(postorder))
    pos = findfirst(==(root.val), inorder)
    pos != 1 && (
        root.left = build_tree_inpost(
            @view(inorder[1:(pos - 1)]), @view(postorder[1:(pos - 1)])
        )
    )
    pos != length(postorder) && (
        root.right = build_tree_inpost(
            @view(inorder[(pos + 1):end]), @view(postorder[pos:(end - 1)])
        )
    )
    return root
end
## @lc code=end


text:
struct PolynomialModP
    poly::Polynomial
    prime::Integer
    PolynomialModP() = new(Polynomial())
    PolynomialModP(prime::Integer) = new(Polynomial(),prime)
    PolynomialModP(poly::Polynomial,prime::Integer,dummy::Int) = new(poly,prime)
end

"""
Construct polynomial with a single term or multiple terms.
"""
function PolynomialModP(poly::Polynomial,prime::Integer)
    terms = Term[]
    for t in poly.terms
        if mod(t.coeff,prime) != 0
            push!(terms, mod(t,prime))
        end
    end
    p = Polynomial(terms,1)
    return PolynomialModP(p,prime,1)
end
"""
Push a new term into the polynomial.
"""
#Note that ideally this would throw and error if pushing another term of degree that is already in the polynomial
function push!(p::PolynomialModP, t::Term)
    iszero(t) && return #don't push a zero
    push!(p.poly.terms,t)
end

"""
Pop the leading term out of the polynomial.
"""
pop!(p::PolynomialModP)= popat!(p.poly.terms,1)
"""
Check if the polynomial is zero.
"""
iszero(p::PolynomialModP)::Bool = isempty(p.poly.terms)

"""
The negative of a polynomial.
"""

function -(p::PolynomialModP)::PolynomialModP
    pol = p.poly
    pol_out = Polynomial()
    for i in pol
        push!(pol_out,Term(-i.coeff,i.degree))
    end
    return PolynomialModP(pol_out,p.prime)
end
"""
Leading term of PolynomialModP
"""
function leading(p::PolynomialModP)::Term
    return p.poly.terms[1]
end
"""
Degree of PolynomialModP
"""
function degree(p::PolynomialModP)::Int
    if iszero(p)
        return 0
    else
        return leading(p.poly).degree
    end
end

"""
Add two polynomials.
"""
function +(p1::PolynomialModP, p2::PolynomialModP)::PolynomialModP
    prime = p1.prime
    p3 = PolynomialModP(prime)
    while !iszero(p1) && !iszero(p2)
        t1, t2 = leading(p1), leading(p2)
        if t1.degree == t2.degree
            degree = t1.degree
            lead_p1 = pop!(p1)
            lead_p2 = pop!(p2)
            if mod(t1.coeff+t2.coeff,prime) != 0
                push!(p3, Term(mod(lead_p1.coeff+lead_p2.coeff,prime),lead_p1.degree))
            end
        elseif t1.degree < t2.degree
            push!(p3,pop!(p2))
        else
            push!(p3,pop!(p1))
        end
    end
    while !iszero(p1)
        push!(p3,pop!(p1))
    end
    while !iszero(p2)
        push!(p3,pop!(p2))
    end
    return p3
end

"""
Add a polynomial and a term.
"""
+(p::PolynomialModP, t::Term) = PolynomialModP(p.polynomial + Polynomial(t),p.prime)
+(t::Term, p::PolynomialModP) = PolynomialModP(p.polynomial + t,p.prime)

"""
Add a polynomial and an integer.
"""
+(p::PolynomialModP, n::Int) = PolynomialModP(p.polynomial + Term(n,0),p.prime)
+(n::Int, p::PolynomialModP) = PolynomialModP(p.polynomial + Term(n,0),p.prime)
"""
Subtraction of two polynomials.
"""

function -(p1::PolynomialModP,p2::PolynomialModP)::PolynomialModP
    @assert p1.prime == p2.prime
    pol1 = p1.poly
    pol2 = p2.poly
    neg_pol2 = -pol2
    pol_out = pol1 + neg_pol2
    return PolynomialModP(pol_out,p1.prime)
end

"""
Multiply two polynomials.
"""
function *(p1::PolynomialModP, p2::PolynomialModP)::PolynomialModP
    @assert p1.prime == p2.prime
    p_out = Polynomial()
    for t in p1.poly
        p_out = p_out + (t * p2.poly)
    end
    return PolynomialModP(p_out,p1.prime)
end

"""
Division of 2 PolynomialModP types
"""
function divide(sp1::PolynomialModP,sp2::PolynomialModP)
    @assert sp1.prime == sp2.prime
    iszero(sp2.poly) && throw(DivideError())
    p1, p2 = deepcopy(sp1), deepcopy(sp2)
    q = PolynomialModP(sp1.prime)
    prev_degree = degree(p1)
    while degree(p1) ≥ degree(p2)
        h_sub = (leading(p1)÷leading(p2))(sp1.prime)
        h_sub_2 = Polynomial(h_sub)
        h = PolynomialModP(h_sub_2,sp1.prime )  #syzergy
        p1 = p1 - h*p2
        q = q + h
        prev_degree == degree(p1) && break
        prev_degree = degree(p1)
    end
    p1
    return q, p1
end

"""
Show function for PolynomialModP
"""
function show(io::IO, p::PolynomialModP)
    p = deepcopy(p.poly)
    if iszero(p)
        print(io,"0")
    else
        n = length(p.terms)
        for (i,t) in enumerate(p.terms)
            print(io, t, i != n ? " + " : "")
        end
    end
end
"""
Derivative function for PolynomialModP
"""
function derivative(p::PolynomialModP)::PolynomialModP
    der_p = Polynomial()
    for term in p.poly
        der_term = derivative(term)
        !iszero(der_term) && push!(der_p,der_term)
    end
    return PolynomialModP(der_p,p.prime)
end
"""
Creates the unit polynomial.
"""
#one(prime::Integer)::PolynomialModP = PolynomialModP(Polynomial(one(Term())),prime)
#one(p::PolynomialModP) = one(typeof(p))

"""
Power of a polynomial mod prime.
"""
function pow_mod(p::PolynomialModP, n::Int)
    n < 0 && error("No negative power")
    out = deepcopy(p)
    if n == 1
        return out
    else
        for _ in 2:n
            out *= p
        end
        return out
    end
end

"""
The extended euclid algorithm for polynomials modulo prime.
"""
#Trying to divide by 1 (r=1), need to fix and find out how to go from r being 1 to becoming 0.
#Maybe a remainder of 1 (0*prime + 1), so then you do euclid(1,0).

function extended_euclid_alg(a::PolynomialModP, b::PolynomialModP)
    @assert a.prime == b.prime
    if a.poly == Polynomial(Term(0,0))
        return b, 0, 1
    elseif b.poly == Polynomial(Term(0,0))
        return a, 1, 0
    else
        old_r, r = a, b
        old_s, s = PolynomialModP(Polynomial(Term(1,0)),a.prime), PolynomialModP(Polynomial(Term(0,0)),a.prime)
        old_t, t = PolynomialModP(Polynomial(Term(0,0)),a.prime), PolynomialModP(Polynomial(Term(1,0)),a.prime)
        while !iszero(r)
            q = divide(old_r, r) |> first
            old_r, r = r, old_r-q*r
            old_s, s = s, old_s-q*s
            old_t, t = t, old_t-q*t
        end
        g, s, t = old_r, old_s, old_t
        return g, s, t
    end
end

function test()
    a = PolynomialModP(Polynomial([Term(4,3),Term(2,2),Term(3,1),Term(2,0)]),7)
    b = PolynomialModP(Polynomial([Term(6,8),Term(3,4),Term(5,2),Term(2,1)]),7)
    ans = PolynomialModP(Polynomial([Term(3,11),Term(5,10),Term(4,9),Term(5,8),Term(5,7),Term(6,6),Term(1,5),Term(3,4),Term(5,3),Term(2,2),Term(4,1)]),7)
    @time mult_ans = a*b
    if ans.poly == mult_ans.poly
        return "Multiplication works"
    else
        return "Multiplication doesnt work"
    end
end
"""
Chinese remainder theorem on integers
"""

function chineseremainder(n::Array, a::Array)
    Π = prod(n)
    return mod(sum(ai * invmod(Π ÷ ni, ni) * (Π ÷ ni) for (ni, ai) in zip(n, a)), Π)
end

function CRT(a::Vector{PolynomialModP},n::Vector{Int64})::Polynomial
    c = Polynomial(Term(0,0))
    for k in max(degree(a[1]),degree(a[2])):-1:0
        if k == degree(a[1])
            if iszero(a[1])
                ak = 0
            else
                ak = pop!(a[1]).coeff
            end
        else
            ak = 0
        end
        if k == degree(a[2])
            if iszero(a[2])
                bk = 0
            else
                bk = pop!(a[2]).coeff
            end
        else
            bk = 0
        end
        ck = chineseremainder(n,[ak,bk])
        push!(c.terms,Term(ck,k))
    end
    return c
end

function *(a::Polynomial,b::Polynomial)
    a_height = 0
    for i in a.terms
        if abs(i.coeff) > a_height
            a_height = abs(i.coeff)
        end
    end
    b_height = 0
    for j in b.terms
        if abs(j.coeff) > b_height
            b_height = abs(j.coeff)
        end
    end
    B = 2*a_height*b_height*min(degree(a)+1,degree(b)+1)
    M = 3
    ap = PolynomialModP(a,M)
    bp = PolynomialModP(b,M)
    c = ap*bp
    while M < B
        p = nextprime(M+1)
        ap = PolynomialModP(a,p)
        bp = PolynomialModP(b,p)
        c2 = ap*bp
        c = PolynomialModP(CRT([c,c2],[M,p]),p)
        M = M*p
    end
    return smod(c)
end

"""
Symmetric mod to allow for negative terms
"""
function smod(a::PolynomialModP)::PolynomialModP
    prime = a.prime
    out = PolynomialModP(a.prime)
    for i in a.poly
        if mod(i.coeff,prime) > prime//2
            i = Term(i.coeff - prime, i.degree)
            push!(out.poly,i)
        else
            push!(out,i)
        end
    end
    return out
end

function pow_mod(a::PolynomialModP, n::Int)::PolynomialModP
    n_bit = bitstring(n)
    n_bit_split = split(n_bit,"")
    current = a
    out = PolynomialModP(Polynomial(Term(1,0)),a.prime)
    highest_one = 0
    for i in 1:64
        if n_bit_split[i] == "1"
            highest_one = i
            break
        end
    end
    for j in 64:-1:highest_one
        if n_bit_split[j] == "1"
            out *= current
        end
        current *= current
    end
    return out
end

function pow_mod(a::Polynomial, n::Int, p::Int)::Polynomial
    poly_mod = PolynomialModP(a,p)
    return pow_mod(poly_mod,n).poly
end

|


text:
#using Pkg
# pkg"add Oceananigans CairoMakie"

ENV["GKSwstype"] = "100"

pushfirst!(LOAD_PATH, @__DIR__)

using Printf
using Statistics

using Oceananigans
using Oceananigans.Units
using Oceananigans.OutputReaders: FieldTimeSeries
using Oceananigans.Grids: xnode, ynode, znode

const Lx = 1000kilometers # zonal domain length [m]
const Ly = 2000kilometers # meridional domain length [m]

# Architecture
architecture  = CPU()

# number of grid points
Nx = 100
Ny = 200
Nz = 35

# stretched grid
k_center = collect(1:Nz)
Δz_center = @. 10 * 1.104^(Nz - k_center)

const Lz = sum(Δz_center)

z_faces = vcat([-Lz], -Lz .+ cumsum(Δz_center))
z_faces[Nz+1] = 0

grid = RectilinearGrid(architecture = architecture,
                       topology = (Periodic, Bounded, Bounded),
                       size = (Nx, Ny, Nz),
                       halo = (3, 3, 3),
                       x = (0, Lx),
                       y = (0, Ly),
                       z = z_faces)

@info "Built a grid: $grid."

#####
##### Boundary conditions
#####

α  = 2e-4     # [K⁻¹] thermal expansion coefficient
g  = 9.8061   # [m/s²] gravitational constant
cᵖ = 3994.0   # [J/K]  heat capacity
ρ  = 999.8    # [kg/m³] reference density

parameters = (
    Ly = Ly,
    Lz = Lz,
    Qᵇ = 10/(ρ * cᵖ) * α * g,            # buoyancy flux magnitude [m² s⁻³]
    y_shutoff = 5/6 * Ly,                # shutoff location for buoyancy flux [m]
    τ = 0.2/ρ,                           # surface kinematic wind stress [m² s⁻²]
    μ = 1 / 30days,                      # bottom drag damping time-scale [s⁻¹]
    ΔB = 8 * α * g,                      # surface vertical buoyancy gradient [s⁻²]
    H = Lz,                              # domain depth [m]
    h = 1000.0,                          # exponential decay scale of stable stratification [m]
    y_sponge = 19/20 * Ly,               # southern boundary of sponge layer [m]
    λt = 7.0days                         # relaxation time scale [s]
)

@inline function buoyancy_flux(i, j, grid, clock, model_fields, p)
    y = ynode(Center(), j, grid)
    return ifelse(y < p.y_shutoff, p.Qᵇ * cos(3π * y / p.Ly), 0.0)
end

buoyancy_flux_bc = FluxBoundaryCondition(buoyancy_flux, discrete_form=true, parameters=parameters)

@inline function u_stress(i, j, grid, clock, model_fields, p)
    y = ynode(Center(), j, grid)
    return - p.τ * sin(π * y / p.Ly)
end

u_stress_bc = FluxBoundaryCondition(u_stress, discrete_form=true, parameters=parameters)

@inline u_drag(i, j, grid, clock, model_fields, p) = @inbounds - p.μ * p.Lz * model_fields.u[i, j, 1]
@inline v_drag(i, j, grid, clock, model_fields, p) = @inbounds - p.μ * p.Lz * model_fields.v[i, j, 1]

u_drag_bc = FluxBoundaryCondition(u_drag, discrete_form=true, parameters=parameters)
v_drag_bc = FluxBoundaryCondition(v_drag, discrete_form=true, parameters=parameters)

b_bcs = FieldBoundaryConditions(top = buoyancy_flux_bc)

u_bcs = FieldBoundaryConditions(top = u_stress_bc, bottom = u_drag_bc)
v_bcs = FieldBoundaryConditions(bottom = v_drag_bc)

#####
##### Coriolis
#####

const f = -1e-4
const β =  1e-11
coriolis = BetaPlane(f₀ = f, β = β)

#####
##### Forcing and initial condition
#####

@inline initial_buoyancy(z, p) = p.ΔB * (exp(z / p.h) - exp(-p.Lz / p.h)) / (1 - exp(-p.Lz / p.h))
@inline mask(y, p) = max(0.0, y - p.y_sponge) / (Ly - p.y_sponge)

@inline function buoyancy_relaxation(i, j, k, grid, clock, model_fields, p)
    timescale = p.λt
    y = ynode(Center(), j, grid)
    z = znode(Center(), k, grid)
    target_b = initial_buoyancy(z, p)
    b = @inbounds model_fields.b[i, j, k]

    return - 1 / timescale  * mask(y, p) * (b - target_b)
end

Fb = Forcing(buoyancy_relaxation, discrete_form = true, parameters = parameters)

# closure

κh = 0.5e-5 # [m²/s] horizontal diffusivity
νh = 30.0   # [m²/s] horizontal viscocity
κz = 0.5e-5 # [m²/s] vertical diffusivity
νz = 3e-4   # [m²/s] vertical viscocity

closure = AnisotropicDiffusivity(νh=νh, νz=νz, κh=κh, κz=κz)

convective_adjustment = ConvectiveAdjustmentVerticalDiffusivity(convective_κz = 1.0,
                                                                convective_νz = 0.0)

#####
##### Model building
#####

@info "Building a model..."

model = HydrostaticFreeSurfaceModel(architecture = architecture,
                                    grid = grid,
                                    free_surface = ImplicitFreeSurface(),
                                    momentum_advection = WENO5(grid = grid),
                                    tracer_advection = WENO5(grid = grid),
                                    buoyancy = BuoyancyTracer(),
                                    coriolis = coriolis,
                                    closure = (closure, convective_adjustment),
                                    tracers = :b,
                                    boundary_conditions = (b=b_bcs, u=u_bcs, v=v_bcs),
                                    forcing = (; b=Fb,)
                                    )

@info "Built $model."

#####
##### Initial conditions
#####

# resting initial condition
ε(σ) = σ * randn()
bᵢ(x, y, z) = parameters.ΔB * ( exp(z/parameters.h) - exp(-Lz/parameters.h) ) / (1 - exp(-Lz/parameters.h)) + ε(1e-8)

set!(model, b=bᵢ)

#####
##### Simulation building
#####

simulation = Simulation(model, Δt=Δt₀, stop_time=stop_time)

# add timestep wizard callback
wizard = TimeStepWizard(cfl=0.1, max_change=1.1, max_Δt=20minutes)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(20))

# add progress callback
wall_clock = [time_ns()]

function print_progress(sim)
    @printf("[%05.2f%%] i: %d, t: %s, wall time: %s, max(u): (%6.3e, %6.3e, %6.3e) m/s, next Δt: %s\n",
            100 * (sim.model.clock.time / sim.stop_time),
            sim.model.clock.iteration,
            prettytime(sim.model.clock.time),
            prettytime(1e-9 * (time_ns() - wall_clock[1])),
            maximum(abs, sim.model.velocities.u),
            maximum(abs, sim.model.velocities.v),
            maximum(abs, sim.model.velocities.w),
            prettytime(sim.Δt))

    wall_clock[1] = time_ns()

    return nothing
end

simulation.callbacks[:print_progress] = Callback(print_progress, IterationInterval(20))

#####
##### Diagnostics
#####

u, v, w = model.velocities
b = model.tracers.b

ζ = ComputedField(∂x(v) - ∂y(u))

B = AveragedField(b, dims=1)
V = AveragedField(v, dims=1)
W = AveragedField(w, dims=1)

b′ = b - B
v′ = v - V
w′ = w - W

v′b′ = AveragedField(v′ * b′, dims=1)
w′b′ = AveragedField(w′ * b′, dims=1)

outputs = (; b, ζ, w)

averaged_outputs = (; v′b′, w′b′, B)

#####
##### Build checkpointer and output writer
#####

simulation.output_writers[:checkpointer] = Checkpointer(model,
                                                        schedule = TimeInterval(100days),
                                                        prefix = "abernathey_channel",
                                                        force = true)

simulation.output_writers[:fields] = JLD2OutputWriter(model, outputs,
                                                      schedule = TimeInterval(5days),
                                                      prefix = "abernathey_channel",
                                                      field_slicer = nothing,
                                                      verbose = true,
                                                      force = true)

simulation.output_writers[:averages] = JLD2OutputWriter(model, averaged_outputs,
                                                        schedule = AveragedTimeInterval(1days, window=1days, stride=1),
                                                        prefix = "abernathey_channel_averages",
                                                        verbose = true,
                                                        force = true)

@info "Running the simulation..."

try
    run!(simulation, pickup=false)
catch err
    @info "run! threw an error! The error message is"
    showerror(stdout, err)
end

# #####
# ##### Visualization
# #####

using Plots

grid = RectilinearGrid(architecture = CPU(),
                       topology = (Periodic, Bounded, Bounded),
                       size = (grid.Nx, grid.Ny, grid.Nz),
                       halo = (3, 3, 3),
                       x = (0, grid.Lx),
                       y = (0, grid.Ly),
                       z = z_faces)

xζ, yζ, zζ = nodes((Face, Face, Center), grid)
xc, yc, zc = nodes((Center, Center, Center), grid)
xw, yw, zw = nodes((Center, Center, Face), grid)

j′ = round(Int, grid.Ny / 2)
y′ = yζ[j′]

b_timeseries = FieldTimeSeries("abernathey_channel.jld2", "b", grid=grid)
ζ_timeseries = FieldTimeSeries("abernathey_channel.jld2", "ζ", grid=grid)
w_timeseries = FieldTimeSeries("abernathey_channel.jld2", "w", grid=grid)

@show b_timeseries

anim = @animate for i in 1:length(b_timeseries.times)
    b = b_timeseries[i]
    ζ = ζ_timeseries[i]
    w = w_timeseries[i]

    b′ = interior(b) .- mean(b)
    b_xy = b′[:, :, grid.Nz]
    ζ_xy = interior(ζ)[:, :, grid.Nz]
    ζ_xz = interior(ζ)[:, j′, :]
    w_xz = interior(w)[:, j′, :]

    @show bmax = max(1e-9, maximum(abs, b_xy))
    @show ζmax = max(1e-9, maximum(abs, ζ_xy))
    @show wmax = max(1e-9, maximum(abs, w_xz))

    blims = (-bmax, bmax) .* 0.8
    ζlims = (-ζmax, ζmax) .* 0.8
    wlims = (-wmax, wmax) .* 0.8

    blevels = vcat([-bmax], range(blims[1], blims[2], length=31), [bmax])
    ζlevels = vcat([-ζmax], range(ζlims[1], ζlims[2], length=31), [ζmax])
    wlevels = vcat([-wmax], range(wlims[1], wlims[2], length=31), [wmax])

    xlims = (-grid.Lx/2, grid.Lx/2) .* 1e-3
    ylims = (0, grid.Ly) .* 1e-3
    zlims = (-grid.Lz, 0)

    w_xz_plot = contourf(xw * 1e-3, zw, w_xz',
                        xlabel = "x (km)",
                        ylabel = "z (m)",
                        aspectratio = 0.05,
                        linewidth = 0,
                        levels = wlevels,
                        clims = wlims,
                        xlims = xlims,
                        ylims = zlims,
                        color = :balance)

    ζ_xy_plot = contourf(xζ * 1e-3, yζ * 1e-3, ζ_xy',
                        xlabel = "x (km)",
                        ylabel = "y (km)",
                        aspectratio = :equal,
                        linewidth = 0,
                        levels = ζlevels,
                        clims = ζlims,
                        xlims = xlims,
                        ylims = ylims,
                        color = :balance)

    b_xy_plot = contourf(xc * 1e-3, yc * 1e-3, b_xy',
                        xlabel = "x (km)",
                        ylabel = "y (km)",
                        aspectratio = :equal,
                        linewidth = 0,
                        levels = blevels,
                        clims = blims,
                        xlims = xlims,
                        ylims = ylims,
                        color = :balance)

    w_xz_title = @sprintf("w(x, z) at t = %s", prettytime(ζ_timeseries.times[i]))
    ζ_xz_title = @sprintf("ζ(x, z) at t = %s", prettytime(ζ_timeseries.times[i]))
    ζ_xy_title = "ζ(x, y)"
    b_xy_title = "b(x, y)"

    layout = @layout [upper_slice_plot{0.2h}
                    Plots.grid(1, 2)]

    plot(w_xz_plot, ζ_xy_plot,  b_xy_plot, layout = layout, size = (1200, 1200), title = [w_xz_title ζ_xy_title b_xy_title])
end

mp4(anim, "abernathey_channel.mp4", fps = 8) # hide


text:
include("cheb.jl")

mutable struct Grid{I, F, V, M}
    H::F            # meters               Physical
   Ly::F            # meters
   Lϕ::F            # radians
 ϕmid::F            # radians
    a::F            # meters
    N::I            # Integer              Numerical
    y::V
    ϕ::V
    D::M
   D2::M
end

"""

Cartesian version of Grid function.

"""
function Grid(geometry::Cartesian;
     H = 500.0,
    Ly = 20e3,
  ϕmid = π/4,
     a = 6.3781e6,
     Lϕ = Ly/a,
     N = 100
 )

 D, y = cheb(N)
    y = Ly / 2 * y
    ϕ = ϕmid .+ y/a
    D = 2 / Ly * D
   D2 = D * D

 return Grid(H, Ly, Lϕ, ϕmid, a, N, y, ϕ, D, D2)
end

"""

Spherical version of Grid function.

"""
function Grid(geometry::Spherical;
     H = 500.0,
    Lϕ = 3.135731e-3,
  ϕmid = π/4,
     a = 6.3781e6,
     Ly = Lϕ*a,
     N = 100
)

D, ϕ = cheb(N)
   ϕ = Lϕ / 2 * ϕ
   y = ϕ * a
   ϕ = ϕ .+ ϕmid
   D = 2 / Lϕ * D
  D2 = D * D

return Grid(H, Ly, Lϕ, ϕmid, a, N, y, ϕ, D, D2)
end


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

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

# ╔═╡ 1685eaa0-7c38-11eb-0dc2-8bfbea75d947
begin
	using PlutoUI
	PlutoUI.TableOfContents(aside=true)
end

# ╔═╡ c75546c0-47ae-4706-9a73-aad97700584d
# Traces of 5 random matrices. We need the package LinearAlgebra
begin
	using LinearAlgebra
	[tr(rand(n,n)) for n=1:5]
end

# ╔═╡ 1f472327-7696-4070-877b-36fa00b32d11
md"""
# Lightning Round - Basic Features and Commands

In this notebook, we go through basic constructs and commands.

__Competences__

The user should know to start `Julia`, how to exit, learn some features and be able to write simple programs.

__Credits__

This notebook is based on the [Julia Tutorial](https://github.com/JuliaLang/julia-tutorial).
"""

# ╔═╡ 5d1910ad-80e4-4ea0-972f-3cb47ae7f057
md"""
# Julia resources

Julia resources are accessible through the [Julia home page](http://julialang.org).

Please check [Packages](https://julialang.org/packages/), [Documentation](https://docs.julialang.org/en/v1/) and JuliaCon where you will also find links to videos from previous conferences.

Julia is open source and all routines are available on GitHub. You will learn how to make full use of this later in the course.

"""

# ╔═╡ 3275dbbd-8612-4395-adb9-1318cc85d136
md"""
# Execution
To execute cell use `Shift + Enter` or press __Run__ button (circled triangle at the bottom right of the cell).

Pluto notebook is __reactive__, so after any change all affected cells are recomputed.
"""

# ╔═╡ 8f8e9dba-72ba-41e1-b39a-aeef95402555
md"""
## Markdown cells

Possibility to write comments / code / formulas in `Markdown` cells, makes Pluto notebooks ideal for teaching and research. Text is written using _Julia Markdown_, which is _GitHub Markdown_ with additional understanding of basic `LaTeX`.

[Mastering (GitHub) Markdown](https://guides.github.com/features/mastering-markdown/) is a 3-minute read, another short and very good manual is at [link](http://daringfireball.net/projects/markdown/).
"""

# ╔═╡ d9196675-4ec2-4236-a7f1-bcc589e234ab
md"""
## Which version of `Julia` is running?

The output of this command is in the console.
"""

# ╔═╡ 4f2aee7c-6e08-404c-a561-18f1662e5837
versioninfo()

# ╔═╡ b46be41b-edc5-4321-850c-ecb82e13ac7a
md"""
## Quitting

Exiting from `julia>` kernel: `exit()`.

Killing Pluto server `2 × Ctrl C`.
"""

# ╔═╡ dbc23df0-c4f0-4efc-a452-3392e0a18bf1
# exit()

# ╔═╡ 9779cbcd-d500-4383-a7e3-ff909837aaa6
md"""
# Punctuation review

* `[...]` are for indexing, array constructors and __Comprehensions__
* `(...)` are __required__ for functions `quit()`, `tic()`, `toc()`, `help()`
* `{...}` are for arrays
* `#` is for comments
"""

# ╔═╡ 2b973dd0-1e52-44c9-bf22-be092080fa3e
md"""
# Basic indexing
"""

# ╔═╡ 334ea76b-a126-41b0-8945-ba166507b0a1
# Matrix with random entries between 0 and 1
A₀=rand(5,5)

# ╔═╡ 9ee1a335-8060-4c6c-b1c6-4e1699053086
A₀[1,1]

# ╔═╡ 9293e691-1b04-4989-a891-fbd1bd3306fe
# You can index into output directly
rand(5,5)[1:2,3:4]

# ╔═╡ 2d16ab49-9244-4ef6-b5e0-fa5304b367e4
md"""
## Indexing is elegant

If you want to compute the lower right $2\times 2$ block of $A^{10}$, in most languages you need to first compute $B=A^{10}$ and then index into $B$. In Julia, the command is simply
"""

# ╔═╡ 9e98e5ce-3ca2-4bf8-84f1-23dc3e793b9e
(A₀^10)[4:5,4:5] # Parenthesses around A^10 are necessary

# ╔═╡ 6acd3389-37b1-48a7-9a04-795372503bac
md"""
## Comprehensions - elegant array constructors
"""

# ╔═╡ 2078e00d-b3b2-4188-b5be-9835de2ab357
[i for i=1.0:5]

# ╔═╡ 0bc90ce5-9976-498d-930c-e34fea9c7193
x₀=1:10

# ╔═╡ 3d949c09-801b-434d-a972-f590e8de8704
[ x₀[i]+x₀[i+1] for i=1:9 ]

# ╔═╡ 963ea673-6469-4718-b348-a281578c7742
z = [Matrix{Int64}(I,n,n) for n=1:5]  # z is Array of Arrays

# ╔═╡ ae3fb606-c31d-4edd-9d4c-4c3982c86be4
# First element is a 1x1 Array
z[1]

# ╔═╡ e4790db9-99ee-4200-bc30-d80d0494f23f
# What is the fourth element?
z[4]

# ╔═╡ 9efef8f4-130b-4286-b109-a43e47a1d00b
# Another example of a comprehension
A=[ i+j for i=1:5, j=1:5 ]

# ╔═╡ b80eb1ee-0d74-47fe-8fbf-87140fad8e7c
# Notice the promotion
B=[ i+j for i=1:5, j=1.0:5 ]

# ╔═╡ de3245f1-8b93-4e9e-be8c-76f7f096907f
md"""
# Commands `ndims()` and `typeof()`
"""

# ╔═╡ 6cc8eb48-4dbb-4599-b126-b4b89f8376de
ndims(B)

# ╔═╡ cdff612c-ad32-4fc8-aa11-d4501f490f05
# z is a one-dimensional array
ndims(z)

# ╔═╡ 4725ae25-0800-442e-88b5-19243da2a3e9
# Array of Arrays
typeof(z)

# ╔═╡ 48da687f-935d-4262-ad0a-f5e8f334da9a
# z[5] is a two-dimensional array
typeof(z[5])

# ╔═╡ 0a918be5-070d-4f41-a09d-dc5d3238189b
typeof(A)

# ╔═╡ 84e88328-de83-49ed-89f8-f9ba13f1ea52
md"""
# Vectors are 1-dimensional arrays

See [Multi-dimensional arrays](https://docs.julialang.org/en/v1/manual/arrays/) for more.
"""

# ╔═╡ a232e8d5-1a28-4a85-9cd4-e58bb79b77a0
# This is 2-dimensional array
v=rand(5,1)

# ╔═╡ 851065e9-d701-46ae-9f22-d6ef0975694f
# This is an 1-dimensional array or vector
vv=vec(v)

# ╔═╡ 029ac818-a082-4454-98fc-b4a9b94b719a
# Notice that they are different
v==vv

# ╔═╡ af54e7ef-3512-4cdf-bea5-e7848299638d
# Again a promotion
v-vv

# ╔═╡ a5268ef8-425e-4b8a-b62a-44af4a4993b7
# This is again a vector
w=rand(5)

# ╔═╡ 24dcfabb-b9a1-4c60-afe2-2b29ab653b6d
# First column is a 5 x 1 matrix, second column is a vector of length 5
Mv=[v w]

# ╔═╡ 59bebc62-f680-4d08-9ed9-febe8d8e58d6
# # Matrix columns and rows are extracted as vectors
x=Mv[:,1]

# ╔═╡ 270b5108-2f91-497e-9f01-187c8c3f4d6a
y=Mv[:,2]

# ╔═╡ eb2954a3-b774-45bb-9776-8015daed12e1
x==v # The types differ

# ╔═╡ 7ad96bf2-08c4-43e7-8d79-5aecc90ff9d8
y==w

# ╔═╡ f56b5ed0-7c29-11eb-2b4f-3d73065ae542
# Row is also a vector
Mv[1,:]

# ╔═╡ 1d361187-3759-41a0-9760-a234bfc054f1
# Transpose of a matrix is adjoint
v'

# ╔═╡ 0dd59a04-f872-4c56-8a99-d586086b5802
# Transpose of a vector is adjoint
w'

# ╔═╡ 5cb26aec-503d-4b62-87aa-6d695a2e7c31
md"""
## 1D and 2D arrays
"""

# ╔═╡ c5552f19-1f48-4910-bf92-ef4d1816fcba
# Range
w₁=1.0:5

# ╔═╡ 29d32710-3f61-4dad-b1a2-41cca022de6c
A*w₁

# ╔═╡ 5a89699f-e40f-49a0-ae31-561466963ce3
# Vector
w₂=collect(1.0:5)

# ╔═╡ 19c4f76b-a9f2-4c78-b0a3-0f2bf3109b6f
A*w₂

# ╔═╡ d31925e5-5c02-4386-8181-d308b3f937f4
# This returns a 2-dimensional array - v is a 5 x 1 array
A*v

# ╔═╡ c919d9d4-079e-490a-b858-762bac13c648
# Usual matrix multiplication formula does not work. Why?
B₁=[A[i,:]*A[:,j] for i=1:5, j=1:5]

# ╔═╡ f4ba5280-7cd2-11eb-1cb5-e3bf3534f077
A

# ╔═╡ b5263af0-bd84-4086-9481-eaf9e8b6d402
# Rows and columns are both 1D vectors - must use dot product!
B₂=[A[i,:]⋅A[:,j] for i=1:5, j=1:5]

# ╔═╡ 19abf655-4f37-44aa-9557-e776fb9e5392
inv(B₂)

# ╔═╡ 30e36ce9-896c-4714-b82c-17736ed6b508
md"""
# `ones()`, `one()`, `Uniform Scaling()`, `zero()`, and `zeros()`

Notice that the output type depends on the argument. This is a general Julia feature called `Multiple dispatch`
and will be explained later in more detail.
"""

# ╔═╡ 91294cdb-1bd9-4e36-8244-19ed3f3d9101
# The output type depends on the argument. Float64 is the default.
ones(3,5), ones(Int,5)

# ╔═╡ 7b86d9a0-7c38-11eb-3ca5-df8e517f3a99
zeros(ComplexF64,3,5)

# ╔═╡ 6b28a250-7c38-11eb-200d-c7c0836084f0
zero(ComplexF64)

# ╔═╡ 99a7b5d0-7c38-11eb-291e-1deab770ea9a
one(Int)

# ╔═╡ d46bd8e2-0c27-4e76-98a8-85d73940ff92
zeros(3,5)

# ╔═╡ f3253ec4-3aaa-4b51-bded-4c37b17e6208
zeros(5)

# ╔═╡ a3ee2c40-7c38-11eb-0c7a-359758b3af8c
# I is automatically recognised
zeros(Int,3,3)+I

# ╔═╡ ae870870-7c38-11eb-1b1f-475e0119efa2
# Multiplicative identity
one(rand(3,3))

# ╔═╡ 1a6706f5-d464-4b4d-856a-c969aadd44e1
md"""
# Complex numbers

`i` is too valuable symbol for loops, so Julia uses `im` for the complex unit.
"""

# ╔═╡ 72faab81-f1ce-4d66-9b4d-348abaecd87d
im

# ╔═╡ 8a6ddee3-9f38-4074-83cf-64e919d2ba2e
2im

# ╔═╡ 8d62321c-985e-4399-9ea7-d89c04b430ac
typeof(2.0im)

# ╔═╡ 809b89d6-7868-4ca3-8012-ffd10aca0cd7
# Another way of defining complex numbers
complex(3,4)

# ╔═╡ c5b96be9-2b53-4cbc-aad4-2b0ef59b84df
# Notice promotion
complex(3,4.0)

# ╔═╡ e71d851f-89e9-4c4a-9e41-0c4f0a517013
# This produces an error (like in any other language)
sqrt(-1)

# ╔═╡ 13761e35-aeff-4463-8a3a-3be47262b889
# and this is fine.
sqrt(complex(-1))

# ╔═╡ 37c50436-f17c-4cbe-a8cd-41cc48555596
md"""
# Ternary operator
Let us define our version of the sign function
"""

# ╔═╡ 8c7c1862-4a10-498d-a4aa-bcd0ac00ad62
si(x) = (x>0) ? 1 : -1

# ╔═╡ 7f606319-7cbc-4338-9fcd-c00a5f6a32f0
si(-13)

# ╔═╡ d35fefc3-a529-4021-a91a-08bc5bf80ebb
md"""
This is equivalent to:
"""

# ╔═╡ d986ebe2-76ec-42bf-8c6e-d486841fbe36
function si₁(x)
    if x>0
        return 1
    else
        return -1
    end
end

# ╔═╡ 634626ed-739d-4900-b639-b1447ec9ee43
si(pi-8), si(0), si₁(0.0)

# ╔═╡ f9968b8a-6474-4f1c-b65c-e0bc2187a70e
md"""
The expressions can be nested:
"""

# ╔═╡ 3a1e5e6a-6235-498c-a3e4-c19e4536d63d
# now si(0) is 0
si₂(x) = (x>0) ? 1 : ((x<0) ? -1 : 0)

# ╔═╡ 1a842df8-274a-4ed9-afa7-c1c7430f9283
si₂(π-8), si₂(0)

# ╔═╡ 28455bbd-5bfc-4f98-be4f-c26311f7e327
md"""
# Typing
Special mathematical (LaTeX) symbols can be used (like $\alpha$, $\Xi$, $\pi$, $\oplus$, $\cdot$, etc.). The symbol in both, the notebook and command line version, is produced by writing
LaTeX command followed by `<Tab>`.

Subscripts and superscripts are written as e.g., x\\_m`<TAB>`, x\^3`<TAB>`.

Julia uses full UNICODE set as variable names.
"""

# ╔═╡ 312bff26-67c8-4b5c-b61f-b8d66ece718a
Ξ = 8; Ψ  = 6; Γ = Ξ ⋅ Ψ; 🐮=√2

# ╔═╡ f0c39851-7ac9-4959-a48e-2612211de173
typeof(Γ)

# ╔═╡ a1b56d48-45c2-4483-81bb-f5829ff10e7e
begin
	ω₁=7; xᵏ=23
	ω₁*xᵏ
end

# ╔═╡ 6394216f-fcae-4dc8-8420-432d0b273621
md"""
# Writing a program and running a file

In Julia, the output of command is not displayed, unless explicitely required.
Only the output of the last command is given.

To display results you can use commands `@show` or `println()`.

Consider the file `deploy.jl` with the following code
```julia
n=map(Int,ARGS[1])          # take one integer argument
println(rand(1:n,n,n))  # generate and print n x n matrix of random integers between 1 and n
@show b=3               # set b to 3 and show the result
c=4                     # set c to 4
```
Running the program in the unix shell gives
```julia
$ julia deploy.jl 5
[1 3 2 4 1
 5 3 1 1 4
 5 4 2 2 5
 3 1 2 3 4
 4 4 5 4 4]
b = 3 => 3
```
Notice that the result of the last command (_c_) is not displayed.

Similarly, the program can be converted to executable and run directly,
without referencing `julia` in the command line.
The refernece to `julia` must be added in the first line, as in the file `deploy1.jl`:
```julia
#!/usr/bin/julia
n=map(Int,ARGS[1])
println(rand(1:n,n,n))
@show b=3
c=4
```
In the shell do:
```
$ chmod +x deploy1.jl
$ ./deploy1.jl 5
[4 5 3 2 5
 4 2 1 5 1
 3 2 4 5 1
 2 4 4 3 1
 3 4 5 3 3]
b = 3 => 3
```

Finally, to run the same program in `julia` shell or `IJulia`, the input has to be changed, as in the file `deploy2.jl`:
```julia
n=parse(Int,readline())
println(rand(1:n,n,n))
@show b=3
c=4
```
__Notice that now the result of the last line is displayed by default__  - in this case it is `4`, the values of `c`. The output of the random matrix and of `b` is forced.

Input in Pluto notebook must be handled differently:
"""

# ╔═╡ 505f8a60-7c3d-11eb-1ce1-c9b7cf6e6c6a
@bind s TextField()

# ╔═╡ 34d0a3ae-7c3d-11eb-0d94-35472ded534e
n = (s=="" ? "" : parse(Int,s))

# ╔═╡ 97de847a-2457-47f1-96a2-6f1214e81a5e
md"""
# Running external programs and unix pipe

## `run()` - calling external program
"""

# ╔═╡ 67d1137f-146b-49bd-94cb-cc666102cb15
md"""
Notice, that this is not a gret help, Julia has much better commands for this.
"""

# ╔═╡ 6f06c3c6-53e2-4374-aeaa-fbcad7db7204
# This calls the unix Calendar program
# run(`cal`)

# ╔═╡ fe3af7f1-bde9-4a27-99b6-665386b38af9
# The pipe is '|>' instead of usual '|'
# run(pipeline(`cal`,`grep Sa`))

# ╔═╡ 081ccb66-0b04-41b1-bae5-9d3e1f9b5c69
md"""
## `ccall()` - calling `C` program
"""

# ╔═╡ 37e4a4d1-99a7-49ff-8daf-6b204aa4a3b5
# Simple version
ccall(:clock,Int,())

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

[compat]
PlutoUI = "~0.7.10"
"""

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

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

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

[[HypertextLiteral]]
git-tree-sha1 = "72053798e1be56026b81d4e2682dbe58922e5ec9"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.0"

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

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

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

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

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

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

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

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "438d35d2d95ae2c5e8780b330592b6de8494e779"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.3"

[[PlutoUI]]
deps = ["Base64", "Dates", "HypertextLiteral", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "Suppressor"]
git-tree-sha1 = "26b4d16873562469a0a1e6ae41d90dec9e51286d"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.10"

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

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

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

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

[[Suppressor]]
git-tree-sha1 = "a819d77f31f83e5792a76081eee1ea6342ab8787"
uuid = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
version = "0.2.0"

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

# ╔═╡ Cell order:
# ╟─1f472327-7696-4070-877b-36fa00b32d11
# ╠═1685eaa0-7c38-11eb-0dc2-8bfbea75d947
# ╟─5d1910ad-80e4-4ea0-972f-3cb47ae7f057
# ╟─3275dbbd-8612-4395-adb9-1318cc85d136
# ╟─8f8e9dba-72ba-41e1-b39a-aeef95402555
# ╟─d9196675-4ec2-4236-a7f1-bcc589e234ab
# ╠═4f2aee7c-6e08-404c-a561-18f1662e5837
# ╟─b46be41b-edc5-4321-850c-ecb82e13ac7a
# ╠═dbc23df0-c4f0-4efc-a452-3392e0a18bf1
# ╟─9779cbcd-d500-4383-a7e3-ff909837aaa6
# ╟─2b973dd0-1e52-44c9-bf22-be092080fa3e
# ╠═334ea76b-a126-41b0-8945-ba166507b0a1
# ╠═9ee1a335-8060-4c6c-b1c6-4e1699053086
# ╠═9293e691-1b04-4989-a891-fbd1bd3306fe
# ╟─2d16ab49-9244-4ef6-b5e0-fa5304b367e4
# ╠═9e98e5ce-3ca2-4bf8-84f1-23dc3e793b9e
# ╟─6acd3389-37b1-48a7-9a04-795372503bac
# ╠═2078e00d-b3b2-4188-b5be-9835de2ab357
# ╠═c75546c0-47ae-4706-9a73-aad97700584d
# ╠═0bc90ce5-9976-498d-930c-e34fea9c7193
# ╠═3d949c09-801b-434d-a972-f590e8de8704
# ╠═963ea673-6469-4718-b348-a281578c7742
# ╠═ae3fb606-c31d-4edd-9d4c-4c3982c86be4
# ╠═e4790db9-99ee-4200-bc30-d80d0494f23f
# ╠═9efef8f4-130b-4286-b109-a43e47a1d00b
# ╠═b80eb1ee-0d74-47fe-8fbf-87140fad8e7c
# ╠═de3245f1-8b93-4e9e-be8c-76f7f096907f
# ╠═6cc8eb48-4dbb-4599-b126-b4b89f8376de
# ╠═cdff612c-ad32-4fc8-aa11-d4501f490f05
# ╠═4725ae25-0800-442e-88b5-19243da2a3e9
# ╠═48da687f-935d-4262-ad0a-f5e8f334da9a
# ╠═0a918be5-070d-4f41-a09d-dc5d3238189b
# ╟─84e88328-de83-49ed-89f8-f9ba13f1ea52
# ╠═a232e8d5-1a28-4a85-9cd4-e58bb79b77a0
# ╠═851065e9-d701-46ae-9f22-d6ef0975694f
# ╠═029ac818-a082-4454-98fc-b4a9b94b719a
# ╠═af54e7ef-3512-4cdf-bea5-e7848299638d
# ╠═a5268ef8-425e-4b8a-b62a-44af4a4993b7
# ╠═24dcfabb-b9a1-4c60-afe2-2b29ab653b6d
# ╠═59bebc62-f680-4d08-9ed9-febe8d8e58d6
# ╠═270b5108-2f91-497e-9f01-187c8c3f4d6a
# ╠═eb2954a3-b774-45bb-9776-8015daed12e1
# ╠═7ad96bf2-08c4-43e7-8d79-5aecc90ff9d8
# ╠═f56b5ed0-7c29-11eb-2b4f-3d73065ae542
# ╠═1d361187-3759-41a0-9760-a234bfc054f1
# ╠═0dd59a04-f872-4c56-8a99-d586086b5802
# ╠═5cb26aec-503d-4b62-87aa-6d695a2e7c31
# ╠═c5552f19-1f48-4910-bf92-ef4d1816fcba
# ╠═29d32710-3f61-4dad-b1a2-41cca022de6c
# ╠═5a89699f-e40f-49a0-ae31-561466963ce3
# ╠═19c4f76b-a9f2-4c78-b0a3-0f2bf3109b6f
# ╠═d31925e5-5c02-4386-8181-d308b3f937f4
# ╠═c919d9d4-079e-490a-b858-762bac13c648
# ╠═f4ba5280-7cd2-11eb-1cb5-e3bf3534f077
# ╠═b5263af0-bd84-4086-9481-eaf9e8b6d402
# ╠═19abf655-4f37-44aa-9557-e776fb9e5392
# ╟─30e36ce9-896c-4714-b82c-17736ed6b508
# ╠═91294cdb-1bd9-4e36-8244-19ed3f3d9101
# ╠═7b86d9a0-7c38-11eb-3ca5-df8e517f3a99
# ╠═6b28a250-7c38-11eb-200d-c7c0836084f0
# ╠═99a7b5d0-7c38-11eb-291e-1deab770ea9a
# ╠═d46bd8e2-0c27-4e76-98a8-85d73940ff92
# ╠═f3253ec4-3aaa-4b51-bded-4c37b17e6208
# ╠═a3ee2c40-7c38-11eb-0c7a-359758b3af8c
# ╠═ae870870-7c38-11eb-1b1f-475e0119efa2
# ╟─1a6706f5-d464-4b4d-856a-c969aadd44e1
# ╠═72faab81-f1ce-4d66-9b4d-348abaecd87d
# ╠═8a6ddee3-9f38-4074-83cf-64e919d2ba2e
# ╠═8d62321c-985e-4399-9ea7-d89c04b430ac
# ╠═809b89d6-7868-4ca3-8012-ffd10aca0cd7
# ╠═c5b96be9-2b53-4cbc-aad4-2b0ef59b84df
# ╠═e71d851f-89e9-4c4a-9e41-0c4f0a517013
# ╠═13761e35-aeff-4463-8a3a-3be47262b889
# ╟─37c50436-f17c-4cbe-a8cd-41cc48555596
# ╠═8c7c1862-4a10-498d-a4aa-bcd0ac00ad62
# ╠═7f606319-7cbc-4338-9fcd-c00a5f6a32f0
# ╟─d35fefc3-a529-4021-a91a-08bc5bf80ebb
# ╠═d986ebe2-76ec-42bf-8c6e-d486841fbe36
# ╠═634626ed-739d-4900-b639-b1447ec9ee43
# ╟─f9968b8a-6474-4f1c-b65c-e0bc2187a70e
# ╠═3a1e5e6a-6235-498c-a3e4-c19e4536d63d
# ╠═1a842df8-274a-4ed9-afa7-c1c7430f9283
# ╟─28455bbd-5bfc-4f98-be4f-c26311f7e327
# ╠═312bff26-67c8-4b5c-b61f-b8d66ece718a
# ╠═f0c39851-7ac9-4959-a48e-2612211de173
# ╠═a1b56d48-45c2-4483-81bb-f5829ff10e7e
# ╟─6394216f-fcae-4dc8-8420-432d0b273621
# ╠═505f8a60-7c3d-11eb-1ce1-c9b7cf6e6c6a
# ╠═34d0a3ae-7c3d-11eb-0d94-35472ded534e
# ╠═97de847a-2457-47f1-96a2-6f1214e81a5e
# ╟─67d1137f-146b-49bd-94cb-cc666102cb15
# ╠═6f06c3c6-53e2-4374-aeaa-fbcad7db7204
# ╠═fe3af7f1-bde9-4a27-99b6-665386b38af9
# ╟─081ccb66-0b04-41b1-bae5-9d3e1f9b5c69
# ╠═37e4a4d1-99a7-49ff-8daf-6b204aa4a3b5
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_path: src/Module A - Short Julia Course/L01a Lightning Round.jl

max_stars_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_issues_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_forks_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course


text:
struct Entry
    title::String
    source_location::LineNumberNode
    code::Expr
    func::Function
    used_functions::Set{Symbol}
end

function Entry(title, source_location, code, func)
    used_functions = Set{Symbol}()
    MacroTools.postwalk(code) do x
        if @capture(x, f_(xs__))
            push!(used_functions, Symbol(string(f)))
        end
        if @capture(x, f_(xs__) do; body__; end)
            push!(used_functions, Symbol(string(f)))
        end
        return x
    end
    return Entry(title, source_location, code, func, used_functions)
end

function unique_name(entry::Entry)
    loc = entry.source_location
    filename = splitext(basename(string(loc.file)))[1]
    sep = isempty(entry.title) ? "" : "_"
    return replace(string(filename, "_", loc.line, sep, entry.title), " " => "")
end

function nice_title(entry::Entry)
    isempty(entry.title) || return entry.title
    return unique_name(entry)
end

const DATABASE = Dict{String, Entry}()

function cell_expr(name, code, source)
    key = string(source.file, ":", source.line)
    return quote
        closure = () -> $(esc(code))
        entry = ReferenceTests.Entry(
            $(string(name)),
            $(QuoteNode(source)),
            $(QuoteNode(code)),
            closure
        )
        ReferenceTests.DATABASE[$key] = entry
        # display(closure())
    end
end

macro cell(name, code)
    return cell_expr(name, code, __source__)
end

macro cell(code)
    return cell_expr("", code, __source__)
end

"""
    save_result(path, object)

Helper, to more easily save all kind of results from the test database
"""
function save_result(path::String, scene::Makie.FigureLike)
    FileIO.save(path * ".png", scene)
end

function save_result(path::String, stream::VideoStream)
    FileIO.save(path * ".mp4", stream)
end

function save_result(path::String, object)
    FileIO.save(path, object)
end

function load_database()
    empty!(DATABASE)
    include(joinpath(@__DIR__, "tests/text.jl"))
    include(joinpath(@__DIR__, "tests/attributes.jl"))
    include(joinpath(@__DIR__, "tests/examples2d.jl"))
    include(joinpath(@__DIR__, "tests/examples3d.jl"))
    include(joinpath(@__DIR__, "tests/short_tests.jl"))
    return DATABASE
end

function database_filtered(title_excludes = [], nice_title_excludes = []; functions=[])
    database = ReferenceTests.load_database()
    return filter(database) do (name, entry)
        !(entry.title in title_excludes) &&
        !(nice_title(entry) in nice_title_excludes) &&
        !any(x-> x in entry.used_functions, functions)
    end
end


text:
"""
    fit_LRMoE(Y, X, α_init, model; ...)

Fit an LRMoE model.

# Arguments
- `Y`: A matrix of response.
- `X`: A matrix of covariates.
- `α`: A matrix of logit regression coefficients.
- `model`: A matrix specifying the expert functions.

# Optional Arguments
- `expusure`: an array of numerics, indicating the time invertal over which the count data (if applicable) are collected.
    If `nothing` is provided, it is set to 1.0 for all observations. It is assumed that all continuous expert functions are
    not affected by `exposure`.
- `exact_Y`: `true` or `false` (default), indicating if `Y` is observed exactly or with censoring and truncation.
- `penalty`: `true` (default) or `false`, indicating whether penalty is imposed on the magnitude of parameters.
- `pen_α`: a numeric penalty on the magnitude of logit regression coefficients. Default is 1.0.
- `pen_params`: an array of penalty term on the magnitude of parameters of component distributions/expert functions.
- `ϵ`: Stopping criterion on loglikelihood (stop when the increment is less than `ϵ`). Default is 0.001.
- `α_iter_max`: Maximum number of iterations when updating `α`. Default is 5.
- `ecm_iter_max`: Maximum number of iterations of the ECM algorithm. Default is 200.
- `grad_jump`: **IN DEVELOPMENT**
- `grad_seq`: **IN DEVELOPMENT**
- `print_steps`: `true` (default) or `false`, indicating whether intermediate updates of parameters should be logged.

# Return Values
- `model_result.α_fit`: Fitted values of logit regression coefficients `α`.
- `model_result.comp_dist`: Fitted parameters of expert functions.
- `converge`: `true` or `false`, indicating whether the fitting procedure has converged.
- `iter`: Number of iterations passed in the fitting function.
- `ll`: Loglikelihood of the fitted model (with penalty on the magnitude of parameters).
- `ll_np`: Loglikelihood of the fitted model (without penalty on the magnitude of parameters).
- `AIC`: Akaike Information Criterion (AIC) of the fitted model.
- `BIC`: Bayesian Information Criterion (BIC) of the fitted model.
"""
function fit_LRMoE(Y, X, α_init, model;
                    exposure = nothing,
                    exact_Y = false,
                    penalty = true, pen_α = 1.0, pen_params = nothing,
                    ϵ = 1e-03, α_iter_max = 5, ecm_iter_max = 200,
                    grad_jump = true, grad_seq = nothing,
                    print_steps = true)

    # Convert possible dataframes to arrays
    # Y = Array(Y)
    # X = Array(X)

    if penalty == false
        pen_α = Inf
        pen_params = [LRMoE.no_penalty_init.(model[k,:]) for k in 1:size(model)[1]]
    elseif isnothing(pen_params)
        pen_params = [LRMoE.penalty_init.(model[k,:]) for k in 1:size(model)[1]]
    end

    if isnothing(exposure)
        exposure = fill(1.0, size(X)[1])
    end

    if exact_Y == true
        tmp = fit_exact(Array(Y), Array(X), Array(α_init), model;
                        exposure = exposure,
                        penalty = penalty, pen_α = pen_α, pen_params = pen_params,
                        ϵ = ϵ, α_iter_max = α_iter_max, ecm_iter_max = ecm_iter_max,
                        grad_jump = grad_jump, grad_seq = grad_seq,
                        print_steps = print_steps)
    else
        tmp = fit_main(Array(Y), Array(X), Array(α_init), model;
                        exposure = exposure,
                        penalty = penalty, pen_α = pen_α, pen_params = pen_params,
                        ϵ = ϵ, α_iter_max = α_iter_max, ecm_iter_max = ecm_iter_max,
                        grad_jump = grad_jump, grad_seq = grad_seq,
                        print_steps = print_steps)
    end

    model_result = LRMoESTD(tmp.α_fit, tmp.model_fit)

    fit_result = LRMoESTDFit(model_result, tmp.converge, tmp.iter, tmp.ll, tmp.ll_np, tmp.AIC, tmp.BIC)
end

max_stars_repo_name: sparktseung/JuliaLRMoE

max_issues_repo_name: sparktseung/JuliaLRMoE

max_forks_repo_name: sparktseung/JuliaLRMoE


text:
export legendre_coeff, chebyshev_coeff, jacobi_coeff
export laguerre_coeff, hermite_coeff

function depwarn(old, new)
    msg = """
    $old will be removed in GaussQuadrature 0.4.
    Use $new instead."""
    warn(msg)
end

function legendre_coeff{T<:AbstractFloat}(::Type{T},
                       n::Integer, endpt::EndPt)
    depwarn("legendre_coeff", "legendre_coefs")
    μ0 = convert(T, 2.0)
    a = zeros(T, n)
    b = zeros(T, n)
    for i = 1:n
        b[i] = i / sqrt(convert(T, 4*i^2-1))
    end
    return a, b, μ0
end

function chebyshev_coeff{T<:AbstractFloat}(::Type{T},
                        n::Integer, kind::Integer, endpt::EndPt)
    depwarn("chebyshev_coeff", "chebyshev_coefs")
    muzero = convert(T, pi)
    half = convert(T, 0.5)
    a = zeros(T, n)
    b = fill(half, n)
    if kind == 1
        b[1] = sqrt(half)
    elseif kind == 2
        muzero /= 2
    else
        error("Unsupported value for kind")
    end
    return a, b, muzero
end

function jacobi_coeff{T<:AbstractFloat}(n::Integer, alpha::T,
                                        beta::T, endpt::EndPt)
    depwarn("jacobi_coeff", "jacobi_coefs")
    ab = alpha + beta
    i = 2
    abi = ab + 2
    μ0 = 2^(ab+1) * exp(
             lgamma(alpha+1) + lgamma(beta+1) - lgamma(abi) )
    a = zeros(T, n)
    b = zeros(T, n)
    a[1] = ( beta - alpha ) / abi
    b[1] = sqrt( 4*(alpha+1)*(beta+1) / ( (ab+3)*abi*abi ) )
    a2b2 = beta*beta - alpha*alpha
    for i = 2:n
        abi = ab + 2i
        a[i] = a2b2 / ( (abi-2)*abi )
        b[i] = sqrt( 4i*(alpha+i)*(beta+i)*(ab+i) /
                     ( (abi*abi-1)*abi*abi ) )
    end
    return a, b, μ0
end

function laguerre_coeff{T<:AbstractFloat}(n::Integer, alpha::T,
                                          endpt::EndPt)
    depwarn("laguerre_coeff", "laguerre_coefs")
    @assert endpt in [neither, left]
    μ0 = gamma(alpha+1)
    a = zeros(T, n)
    b = zeros(T, n)
    for i = 1:n
        a[i] = 2i - 1 + alpha
        b[i] = sqrt( i*(alpha+i) )
    end
    return a, b, μ0
end

function hermite_coeff{T<:AbstractFloat}(::Type{T}, n::Integer)
    depwarn("hermite_coeff", "hermite_coefs")
    μ0 = sqrt(convert(T, pi))
    a = zeros(T, n)
    b = zeros(T, n)
    for i = 1:n
        iT = convert(T, i)
        b[i] = sqrt(iT/2)
    end
    return a, b, μ0
end

function custom_gauss_rule{T<:AbstractFloat}(lo::T, hi::T,
         a::Array{T,1}, b::Array{T,1}, μ0::T, endpt::EndPt,
         maxits::Integer=maxiterations[T])
    #
    # On entry:
    #
    # a, b hold the coefficients (as given, for instance, by
    # legendre_coeff) in the three-term recurrence relation
    # for the orthonormal polynomials p_0, p_1, p_2, ... , that is,
    #
    #    b[j] p (x) = (x-a[j]) p   (x) - b[j-1] p   (x).
    #          j                j-1              j-2
    #
    # μ0 holds the zeroth moment of the weight function, that is
    #
    #          / hi
    #         |
    #    μ0 = | w(x) dx.
    #         |
    #         / lo
    #
    # On return:
    #
    # x, w hold the points and weights.
    #
    warn("This method will be removed in GaussQuadrature 0.4")
    n = length(a)
    @assert length(b) == n
    if endpt == left
        if n == 1
            a[1] = lo
        else
            a[n] = solve(n, lo, a, b) * b[n-1]^2 + lo
        end
    elseif endpt == right
        if n == 1
            a[1] = hi
        else
            a[n] = solve(n, hi, a, b) * b[n-1]^2 + hi
        end
    elseif endpt == both
        if n == 1
            error("Must have at least two points for both ends.")
        end
        g = solve(n, lo, a, b)
        t1 = ( hi - lo ) / ( g - solve(n, hi, a, b) )
        b[n-1] = sqrt(t1)
        a[n] = lo + g * t1
    end
    w = zero(a)
    steig!(a, b, w, maxits)
    for i = 1:n
        w[i] = μ0 * w[i]^2
    end
    idx = sortperm(a)
    return a[idx], w[idx]
end

function steig!{T<:AbstractFloat}(d::Array{T,1}, e::Array{T,1},
                                  z::Array{T,1}, maxits::Integer)
    #
    # Finds the eigenvalues and first components of the normalised
    # eigenvectors of a symmetric tridiagonal matrix by the implicit
    # QL method.
    #
    # d[i]   On entry, holds the ith diagonal entry of the matrix.
    #        On exit, holds the ith eigenvalue.
    #
    # e[i]   On entry, holds the [i+1,i] entry of the matrix for
    #        i = 1, 2, ..., n-1.  (The value of e[n] is not used.)
    #        On exit, e is overwritten.
    #
    # z[i]   On exit, holds the first component of the ith normalised
    #        eigenvector associated with d[i].
    #
    # maxits The maximum number of QL iterations.
    #
    # Martin and Wilkinson, Numer. Math. 12: 377-383 (1968).
    # Dubrulle, Numer. Math. 15: 450 (1970).
    # Handbook for Automatic Computation, Vol ii, Linear Algebra,
    #        pp. 241-248, 1971.
    #
    # This is a modified version of the Eispack routine imtql2.
    #
    n = length(z)
    z[1] = 1
    z[2:n] = 0
    e[n] = 0

    if n == 1 # Nothing to do for a 1x1 matrix.
        return
    end
    for l = 1:n
        for j = 1:maxits
            # Look for small off-diagonal elements.
            m = n
            for i = l:n-1
                if abs(e[i]) <= eps(T) * ( abs(d[i]) + abs(d[i+1]) )
                    m = i
                    break
                end
            end
            p = d[l]
            if m == l
                continue
            end
            if j == maxits
                msg = @sprintf("No convergence after %d iterations", j)
                msg *= " (try increasing maxits)"
                error(msg)
            end
            # Form shift
            g = ( d[l+1] - p ) / ( 2 * e[l] )
            r = hypot(g, one(T))
            g = d[m] - p + e[l] / ( g + copysign(r, g) )
            s = one(T)
            c = one(T)
            p = zero(T)
            for i = m-1:-1:l
                f = s * e[i]
                b = c * e[i]
                if abs(f) <  abs(g)
                    s = f / g
                    r = hypot(s, one(T))
                    e[i+1] = g * r
                    c = one(T) / r
                    s *= c
                else
                    c = g / f
                    r = hypot(c, one(T))
                    e[i+1] = f * r
                    s = one(T) / r
                    c *= s
                end
                g = d[i+1] - p
                r = ( d[i] - g ) * s + 2 * c * b
                p = s * r
                d[i+1] = g + p
                g = c * r - b
                # Form first component of vector.
                f = z[i+1]
                z[i+1] = s * z[i] + c * f
                z[i]   = c * z[i] - s * f
            end # loop over i
            d[l] -= p
            e[l] = g
            e[m] = zero(T)
        end # loop over j
    end # loop over l
end

function orthonormal_poly{T<:AbstractFloat}(x::Array{T,1},
                         a::Array{T,1}, b::Array{T,1}, μ0::T)
    # p[i,j] = value at x[i] of orthonormal polynomial of degree j-1.
    m = length(x)
    n = length(a)
    p = zeros(T, m, n+1)
    c = one(T) / sqrt(μ0)
    rb = one(T) / b[1]
    for i = 1:m
        p[i,1] = c
        p[i,2] = rb * ( x[i] - a[1] ) * c
    end
    for j = 2:n
       rb = one(T) / b[j]
       for i = 1:m
           p[i,j+1] = rb * ( (x[i]-a[j]) * p[i,j]
                                - b[j-1] * p[i,j-1] )
       end
    end
    return p
end