text:
Getxyp(DM::AbstractDataModel) = (xdim(DM), ydim(DM), pdim(DM))
Getxyp(DS::AbstractDataSet, model::Function) = (xdim(DS),ydim(DS),pdim(DS,model))
Getxyp(DS::AbstractDataSet, M::ModelMap) = M.xyp

SymbolicArguments(args...) = SymbolicArguments(Getxyp(args...))
function SymbolicArguments(xyp::Tuple{Int,Int,Int})
    (xyp[1] == 1 ? (@variables x)[1] : (@variables x[1:xyp[1]])[1]),
    (xyp[2] == 1 ? (@variables y)[1] : (@variables y[1:xyp[2]])[1]),
    (@variables θ[1:xyp[3]])[1]
end
# function SymbolicArguments(xyp::Tuple{Int,Int,Int})
#     @variables X[1:xyp[1]] Y[1:xyp[2]] θ[1:xyp[3]] x y
#     X, Y, θ
#     (xyp[1] == 1 ? x : X), (xyp[2] == 1 ? y : Y), θ
# end

ToExpr(DM::AbstractDataModel; timeout::Real=5) = ToExpr(Predictor(DM), (xdim(DM), ydim(DM), pdim(DM)); timeout=timeout)
ToExpr(DS::AbstractDataSet, model::Function; timeout::Real=5) = ToExpr(model, (xdim(DS),ydim(DS),pdim(DS,model)); timeout=timeout)
ToExpr(DS::AbstractDataSet, M::ModelMap; timeout::Real=5) = ToExpr(M.Map, M.xyp; timeout=timeout)
ToExpr(M::ModelMap; timeout::Real=5) = ToExpr(M.Map, M.xyp; timeout=timeout)
ToExpr(M::ModelMap, xyp::Tuple{Int,Int,Int}; timeout::Real=5) = xyp == M.xyp ? ToExpr(M.Map, M.xyp; timeout=timeout) : throw("Inconsistent xyp information.")

function ToExpr(model::Function, xyp::Tuple{Int,Int,Int}; timeout::Real=5)
    X, Y, θ = SymbolicArguments(xyp)
    KillAfter(model, X, θ; timeout=timeout)
    # Add option for models which are already inplace
end

function SymbolicModel(DM::AbstractDataModel)
    expr = @suppress_err ToExpr(DM)
    isnothing(expr) ? "Unable to represent given model symbolically." : "y(x,θ) = $expr"
end

function SymbolicdModel(DM::AbstractDataModel)
    if !GeneratedFromSymbolic(dPredictor(DM))
        println("Given Model jacobian not symbolic. Trying to apply OptimizedDM() first.")
        odm = OptimizedDM(DM)
        if isnothing(ToExpr(odm))
            return "Unable to represent given jacobian symbolically."
        else
            X, Y, θ = SymbolicArguments(odm)
            return "(∂y/∂θ)(x,θ) = $(dPredictor(odm)(X, θ))"
        end
    else
        X, Y, θ = SymbolicArguments(DM)
        return "(∂y/∂θ)(x,θ) = $(dPredictor(DM)(X, θ))"
    end
end

function Optimize(DM::AbstractDataModel; inplace::Bool=false, timeout::Real=5, parallel::Bool=false, kwargs...)
    Optimize(Data(DM), Predictor(DM); inplace=inplace, timeout=timeout, parallel=parallel, kwargs...)
end
function Optimize(DS::AbstractDataSet, model::ModelOrFunction; inplace::Bool=false, timeout::Real=5, parallel::Bool=false, kwargs...)
    Optimize(model, Getxyp(DS, model); inplace=inplace, timeout=timeout, parallel=parallel, kwargs...)
end
function Optimize(M::ModelMap, xyp::Tuple{Int,Int,Int}; inplace::Bool=false, timeout::Real=5, parallel::Bool=false, kwargs...)
    xyp != M.xyp && throw("xyp inconsistent.")
    model, dmodel = Optimize(M.Map, xyp; inplace=inplace, timeout=timeout, parallel=parallel, kwargs...)
    ModelMap(model, M), ModelMap(dmodel, M)
end
function Optimize(model::Function, xyp::Tuple{Int,Int,Int}; inplace::Bool=false, timeout::Real=5, parallel::Bool=false, kwargs...)
    modelexpr = ToExpr(model, xyp; timeout=timeout) |> Symbolics.simplify
    isnothing(modelexpr) && return nothing, nothing

    X, Y, θ = SymbolicArguments(xyp)

    # Need to make sure that modelexpr is of type Vector{Num}, not just Num
    !(modelexpr isa AbstractVector{<:Num}) && (modelexpr = [modelexpr])
    derivative = Symbolics.jacobian(modelexpr, θ; simplify=true)

    ExprToModelMap(X, θ, modelexpr; inplace=inplace, parallel=parallel, IsJacobian=false, kwargs...), ExprToModelMap(X, θ, derivative; inplace=inplace, parallel=parallel, IsJacobian=true, kwargs...)
end

function ExprToModelMap(X::Union{Num,AbstractVector{<:Num}}, P::AbstractVector{Num}, modelexpr::Union{Num,AbstractArray{<:Num}};
                                                        inplace::Bool=false, parallel::Bool=false, IsJacobian::Bool=false, force_SA::Bool=IsJacobian, kwargs...)
    OptimizedModel = Builder(modelexpr, X, P; inplace=inplace, parallel=parallel, force_SA=force_SA, kwargs...)

    ### Pretty Function names
    if IsJacobian
        # THROWING AWAY KWARGS HERE!
        SymbolicModelJacobian(x::Union{Number,AbstractVector{<:Number}}, θ::AbstractVector{<:Number}; Kwargs...) = OptimizedModel(x, θ)
        function SymbolicModelJacobian!(y::Union{Number,AbstractMatrix{<:Number}}, x::Union{Number,AbstractVector{<:Number}}, θ::AbstractVector{<:Number}; Kwargs...)
            OptimizedModel(y, x, θ)
        end
        return inplace ? SymbolicModelJacobian! : SymbolicModelJacobian
    else
        # THROWING AWAY KWARGS HERE!
        SymbolicModel(x::Union{Number,AbstractVector{<:Number}}, θ::AbstractVector{<:Number}; Kwargs...) = OptimizedModel(x, θ)
        SymbolicModel!(y::Union{Number,AbstractVector{<:Number}}, x::Union{Number,AbstractVector{<:Number}}, θ::AbstractVector{<:Number}; Kwargs...) = OptimizedModel(y, x, θ)
        return inplace ? SymbolicModel! : SymbolicModel
    end
end

function OptimizedDM(DM::AbstractDataModel; kwargs...)
    model, dmodel = Optimize(DM; kwargs...)
    # Very simple models (ydim=1) typically slower after simplification using ModelingToolkit.jl / Symbolics.jl
    !isnothing(dmodel) ? DataModel(Data(DM), Predictor(DM), dmodel, MLE(DM), LogLikeMLE(DM)) : DM
end


text:
# # FFTW plans

# Here we test planning the FFTW and run some simple benchmarks.

# Here are the packages we need.

using FFTW
using Plots
using LinearAlgebra: mul!
using Test
using Random
using BenchmarkTools

@info "Threads: $(FFTW.nthreads())"

# ## Testing setup

# ### The spatial domain and discretization

L = 2π
κ₀ = 2π/L
N = 144 # 2^4 * 3^2
x = y = (L/N):(L/N):L
nothing

# ### Vorticity field for the tests

# We randomly excite a certain number of modes.
rng = Xoshiro(123)
num_modes = 8
vort = sum(
    [
        2κ₀^2 * (kx^2 + ky^2) * (
            ar * cos.(κ₀ * (kx * one.(y) * x' + ky * y * one.(x)'))
            - ai * sin.(κ₀ * (kx * one.(y) * x' + ky * y * one.(x)'))
        )
        for (kx, ky, ar, ai) in zip(
            rand(rng, 1:div(N,4), num_modes),
            rand(rng, 1:div(N,4), num_modes),
            10*rand(rng, num_modes),
            10*rand(rng, num_modes)
        )
    ]
)

vort_hat = rfft(vort)
nothing

# Visualizing the vorticity field

heatmap(x, y, vort, xlabel="x", ylabel="y", title="Vorticity field", titlefont=12)

# ### Planning the FFTWs

#
plan_estimate = plan_rfft(vort); # default is FFTW.ESTIMATE
plan_inv_estimate = plan_irfft(vort_hat, N);
plan_measure = plan_rfft(vort, flags=FFTW.MEASURE);
plan_inv_measure = plan_irfft(vort_hat, N, flags=FFTW.MEASURE);
plan_patient = plan_rfft(vort, flags=FFTW.PATIENT);
plan_inv_patient = plan_irfft(vort_hat, N, flags=FFTW.PATIENT);
plan_exhaustive = plan_rfft(vort, flags=FFTW.EXHAUSTIVE);
plan_inv_exhaustive = plan_irfft(vort_hat, N, flags=FFTW.EXHAUSTIVE);

# ## Sanity tests of the different plans

@testset "Planned FFTW with ESTIMATE" begin
    w_hat = plan_estimate * vort
    w_hat_org = copy(w_hat)
    @test w_hat ≈ rfft(vort)
    @test plan_inv_estimate * w_hat ≈ vort
    @test plan_inv_estimate * w_hat ≈ irfft(rfft(vort), N)
    w_hat_mul = similar(w_hat)
    mul!(w_hat_mul, plan_estimate, vort)
    @test w_hat ≈ w_hat_mul
    vort_back = similar(vort)
    mul!(vort_back, plan_inv_estimate, w_hat) # careful, inverse with mul! may mutate w_hat as well
    @test vort_back ≈ vort
end
nothing

#

@testset "Planned FFTW with MEASURE" begin
    w_hat = plan_measure * vort
    @test w_hat ≈ rfft(vort)
    @test plan_inv_measure * w_hat ≈ vort
    @test plan_inv_measure * w_hat ≈ irfft(rfft(vort), N)
    w_hat_back = similar(w_hat)
    mul!(w_hat_back, plan_measure, vort)
    @test w_hat ≈ w_hat_back
    vort_back = similar(vort)
    mul!(vort_back, plan_inv_measure, w_hat)  # careful, inverse with mul! may mutate w_hat as well
    @test vort_back ≈ vort
end
nothing

#

@testset "Planned FFTW with PATIENT" begin
    w_hat = plan_patient * vort
    @test w_hat ≈ rfft(vort)
    @test plan_inv_patient * w_hat ≈ vort
    @test plan_inv_patient * w_hat ≈ irfft(rfft(vort), N)
    w_hat_back = similar(w_hat)
    mul!(w_hat_back, plan_patient, vort)
    @test w_hat ≈ w_hat_back
    vort_back = similar(vort)
    mul!(vort_back, plan_inv_patient, w_hat)  # careful, inverse with mul! may mutate w_hat as well
    @test vort_back ≈ vort
end
nothing

#

@testset "Planned FFTW with EXHAUSTIVE" begin
    w_hat = plan_exhaustive * vort
    @test w_hat ≈ rfft(vort)
    @test plan_inv_exhaustive * w_hat ≈ vort
    @test plan_inv_exhaustive * w_hat ≈ irfft(rfft(vort), N)
    w_hat_back = similar(w_hat)
    mul!(w_hat_back, plan_exhaustive, vort)
    @test w_hat ≈ w_hat_back
    vort_back = similar(vort)
    mul!(vort_back, plan_inv_exhaustive, w_hat) # careful, inverse with mul! may mutate w_hat as well
    @test vort_back ≈ vort
end
nothing

# ## Timings and allocations

# Let us look at it with `@btime`. **Notice that `mul!` does not allocate and is slightly faster.**

@info "FFTW no plan"
@btime rfft(w) setup = (w = copy(vort));

@info "FFTW plan with ESTIMATE"
@btime p * w setup = (p = $plan_estimate; w = copy($vort));
@info "FFTW plan with ESTIMATE and mul!"
@btime mul!(w_hat, p, w) setup = (
    w_hat = copy($vort_hat);
    p = $plan_estimate;
    w = copy($vort)
);

@info "FFTW plan with MEASURE"
@btime p * w setup = (p = $plan_measure; w = copy($vort));
@info "FFTW plan with MEASURE and mul!"
@btime mul!(w_hat, p, w) setup = (
    w_hat = copy($vort_hat);
    p = $plan_measure;
    w = copy($vort)
);

@info "FFTW plan with PATIENT"
@btime p * w setup = (p = $plan_patient; w = copy($vort));
@info "FFTW plan with PATIENT and mul!"
@btime mul!(w_hat, p, w) setup = (
    w_hat = copy($vort_hat);
    p = $plan_patient;
    w = copy($vort)
);

@info "FFTW plan with EXHAUSTIVE"
@btime p * w setup = (p = $plan_exhaustive; w = copy($vort));
@info "FFTW plan with EXHAUSTIVE and mul!"
@btime mul!(w_hat, p, w) setup = (
    w_hat = copy($vort_hat);
    p = $plan_exhaustive;
    w = copy($vort)
);

# ## Direct and inverse transforms

# Here we compare the timings between direct and inverse transforms. We restrict this to using `mul!`, since this is what we will use at the end. **Notice the direct transform is slightly faster than the inverse transform.**

# Only one note of **warning:** when using `mul!` with an inverse transform plan, the last argument may mutate. Type `@edit mul!(vort, plan_inv_estimate.p, vort_hat)` to see where this happens.

@info "Direct and inverse FFTW plan with ESTIMATE and mul!"
@btime mul!(w_hat, p, w) setup = (
    w_hat = similar($vort_hat);
    p = $plan_estimate;
    w = copy($vort)
);
@btime mul!(w, p, w_hat) setup = (
    w = similar($vort);
    p = $plan_inv_estimate;
    w_hat = copy($vort_hat)
);

@info "Direct and inverse FFTW plan with MEASURE and mul!"
@btime mul!(w_hat, p, w) setup = (
    w_hat = similar($vort_hat);
    p = $plan_measure;
    w = copy($vort)
);
@btime mul!(w, p, w_hat) setup = (
    w = similar($vort);
    p = $plan_inv_measure;
    w_hat = copy($vort_hat)
);

@info "Direct and inverse FFTW plan with PATIENT and mul!"
@btime mul!(w_hat, p, w) setup = (
    w_hat = similar($vort_hat);
    p = $plan_patient;
    w = copy($vort)
);
@btime mul!(w, p, w_hat) setup = (
    w = similar($vort);
    p = $plan_inv_patient;
    w_hat = copy($vort_hat)
);

@info "Direct and inverse FFTW plan with EXHAUSTIVE and mul!"
@btime mul!(w_hat, p, w) setup = (
    w_hat = similar($vort_hat);
    p = $plan_exhaustive;
    w = copy($vort)
);
@btime mul!(w, p, w_hat) setup = (
    w = similar($vort);
    p = $plan_inv_exhaustive;
    w_hat = copy($vort_hat)
);

ext: jl


text:
struct GridWorldMDP <: MDP{GWPos, Symbol}
    g::SimpleGridWorld
end

GridWorldMDP(;size::Tuple{Int, Int}           = (10,10),
              rewards::Dict{GWPos, Float64}   = Dict(GWPos(4,3)=>-10.0, GWPos(4,6)=>-5.0, GWPos(9,3)=>10.0, GWPos(8,8)=>3.0),
              terminate_from::Set{GWPos}      = Set(keys(rewards)),
              tprob::Float64                  = 0.7,
              discount::Float64               = 0.95) = GridWorldMDP(SimpleGridWorld(size = size, rewards = rewards, terminate_from = terminate_from, tprob = tprob, discount = discount))

POMDPs.gen(mdp::GridWorldMDP, s, a, rng = Random.GLOBAL_RNG) = (sp = rand(rng, transition(mdp.g, s, a )), r = POMDPs.reward(mdp.g, s, a))

function POMDPs.initialstate(mdp::GridWorldMDP)
    function istate(rng::AbstractRNG)
        while true
            x, y = rand(rng, 1:mdp.g.size[1]), rand(rng, 1:mdp.g.size[2])
            !(GWPos(x,y) in mdp.g.terminate_from) && return GWPos(x,y)
        end
    end
    return ImplicitDistribution(istate)
end

POMDPs.convert_s(::Type{AbstractArray}, s::GWPos, mdp::GridWorldMDP) = Float32.([s...])

POMDPs.actions(mdp::GridWorldMDP) = POMDPs.actions(mdp.g)

POMDPs.isterminal(mdp::GridWorldMDP, s) = isterminal(mdp.g, s)
POMDPs.discount(mdp::GridWorldMDP) = discount(mdp.g)

function render(mdp::GridWorldMDP, s=GWPos(7,5), a=nothing; color = s->10.0*POMDPs.reward(mdp.g, s), policy= nothing)
    img = POMDPModelTools.render(mdp.g, (s = s,), color = color, policy = isnothing(policy) ? nothing : FunctionPolicy((s) ->  action(policy, convert_s(AbstractArray, s, mdp))))
    tmpfilename = "/tmp/out.png"
    img |> PNG(tmpfilename, 1cm .* mdp.g.size...)
    load(tmpfilename)
end

max_forks_repo_name: ancorso/POMDPGym


text:
# Author: Ashwani Rathee
# Program to add numbers

result = 0
println("Enter 5 numbers line by line")

for number in 1:5
    num = readline()
    num = parse(Int64, num)
    global result+= num

end

println("The sum is :", result)

max_stars_repo_name: ashwani-rathee/Learning--Julia

max_issues_repo_name: ashwani-rathee/Learning--Julia

max_forks_repo_name: ashwani-rathee/Learning--Julia


text:
include("../src/HMat2d.jl");

function laplacian2d(nx,ny)
    M = speye(nx*ny);
    for x=1:nx
        for y=1:ny
            s = x+(y-1)*nx;
            M[s,s] = 4;
            if x > 1
                M[s,s-1] = -1;
            end
            if x < nx
                M[s,s+1] = -1;
            end
            if y > 1
                M[s,s-nx] = -1;
            end
            if y < ny
                M[s,s+nx] = -1;
            end
        end
    end
    return M;
end

EPS = 1e-6;
MaxRank = 2;
minn = 4;

n=16;
A = full(laplacian2d(n,n));
Z2C = Z2Cmapper([n,n],minn,reshape(1:n^2,n,n));
A = inv(full(A[Z2C,Z2C]));
AH16 = HMat2dd2h(A, [n,n], [n,n], EDGE, [0,0], [0,0], 1, EPS, MaxRank, minn);

n=32;
A = full(laplacian2d(n,n));
Z2C = Z2Cmapper([n,n],minn,reshape(1:n^2,n,n));
A = inv(full(A[Z2C,Z2C]));
AH32 = HMat2dd2h(A, [n,n], [n,n], EDGE, [0,0], [0,0], 1, EPS, MaxRank, minn);

n=64;
A = full(laplacian2d(n,n));
Z2C = Z2Cmapper([n,n],minn,reshape(1:n^2,n,n));
A = inv(full(A[Z2C,Z2C]));
AH64 = HMat2dd2h(A, [n,n], [n,n], EDGE, [0,0], [0,0], 1, EPS, MaxRank, minn);

n=16;
v = rand(n^2,5);
u = hvecmat(v',AH16);
tic();
for i=1:100
u = hvecmat(v',AH16);
end
toc();

n=32;
v = rand(n^2,5);
u = hvecmat(v',AH32);
tic();
for i=1:100
u = hvecmat(v',AH32);
end
toc();

n=64;
v = rand(n^2,5);
u = hvecmat(v',AH64);
tic();
for i=1:100
u = hvecmat(v',AH64);
end
toc();

max_stars_repo_name: YingzhouLi/HMat

max_issues_repo_name: YingzhouLi/HMat


text:
# from https://github.com/JuliaLang/julia/pull/36425
function detectwsl()
    Sys.islinux() &&
    isfile("/proc/sys/kernel/osrelease") &&
    occursin(r"Microsoft|WSL"i, read("/proc/sys/kernel/osrelease", String))
end

function browser(url::AbstractString)
    if Sys.isapple()
        Base.run(`open $url`)

    elseif Sys.iswindows() || detectwsl()
        Base.run(`powershell.exe Start $url`)

    elseif Sys.islinux()
        Base.run(`xdg-open $url`)

    else
        printstyled("[FlowAtlas] ", color = :blue)
        println("Could not open default browser")
    end
end

max_stars_repo_name: gszep/flow-explorer

max_issues_repo_name: gszep/flow-explorer

max_forks_repo_name: gszep/flow-explorer


text:
set_ui!(ui) = (setupfolders(); ui[] = hbox(loadfrommemory(ui), hskip(20px), loadbutton(ui)))

set_ui!(ui, t::AbstractString) = set_ui!(ui, loadtable(t))

function set_ui!(ui, t::IndexedTable)
    setupfolders()
    d_obs = Observable{Any}(loadfrommemory(ui))
    s = loadbutton(ui)
    plot_options = dropdownrow(t)
    updater = Observable(0)
    selectors = TableWidgets.selectors(t, updater)
    filtered_data() = (updater[] += 1; selectors[])
    categoricalstyle, continuousstyle = stylechoosers(t)
    style = vcat(categoricalstyle, continuousstyle)
    plot_command = button("Plot")
    spreadsheet_command = button("Spreadsheet")
    save_plot_button = PlotSaver()
    save_table_button = TableSaver()
    plot_buttons = hbox(plot_command, hskip(20px), save_plot_button.button, save_plot_button.name)
    table_buttons = hbox(spreadsheet_command,
                         hskip(20px),
                         save_table_button.button,
                         save_table_button.name,
                         hskip(20px),
                         save_table_button.checkbox)
    plt = Observable{Any}(default_plot())
    spreadsheet =  Observable{Any}(dom"div"(""))
    plt_kwargs = textbox("Insert plot attributes")
    smoother = slider(1:100, label = "smoothing")
    on(x -> spreadsheet[] = build_spreadsheet(filtered_data(), style), observe(spreadsheet_command))
    on(x -> savefig(plt[], filename(save_plot_button)), observe(save_plot_button))
    onany((x, y) -> plt[] = build_plot(filtered_data(), plot_options, style, plt_kwargs, y), observe(plot_command), observe(smoother))
    on(x -> (_save(filtered_data(), style, filename(save_table_button), isselected(save_table_button)); d_obs[] = loadfrommemory(ui)),
        observe(save_table_button))

    left_menu = tabulator(OrderedDict(zip(
        ["Load", "Filter", "Style"],
        [
            dom"div"(pad(1em, d_obs), pad(1em, s)),
            pad(1em, selectors),
            pad(1em, dom"div.columns"(layout.((categoricalstyle, continuousstyle))...))
        ])), index=1)

    plot_view = dom"div"(
        layout(plot_options),
        vskip(2em),
        plot_buttons,
        plt,
        plt_kwargs,
        smoother
    )
    spreadsheet_view = dom"div"(
        table_buttons,
        vskip(20px),
        spreadsheet,
    )

    right_menu = tabulator(OrderedDict("Table" => spreadsheet_view, "Graph" => plot_view), index = 1)

    ui[] = dom"div.columns"(
        dom"div.column.is-two-fifths.has-background-light[style=height:100vh;overflow-y:scroll;overflow-x:hidden]"(
            pad(1em, left_menu)
        ),
        dom"div.column[style=height:100vh;overflow-y:scroll;overflow-x:hidden]"(
            pad(1em, right_menu)
        )
    )
end

function get_ui(args...; kwargs...)
    ui = Observable{Any}("")
    set_ui!(ui, args...; kwargs...)
    ui
end

# To do: make it work not only in Juno but in every other environment

#_display(t) = display(t)
function launch(args...; env = :blink, page = "/", kwargs...)
    ui = get_ui(args...; kwargs...)
    _display(dom"div"(ui), env; page = page)
end

function _display(t, env; page = "/")
    if env == :blink
        w = Window()
        body!(w, t)
        w
    elseif env == :juno
        display(t)
    elseif env == :mux
        webio_serve(WebIO.page(page, req -> t))
    else
        error("Only blink, juno and mux supported so far")
    end
end


text:
include("../ds.jl")
using LinearAlgebra
using Distributions
# using PDMats

struct compact_mnm_dist <: distibution_sample
    α::AbstractArray{Float64,1}
end

#topic_modeling_dist
function log_likelihood!(r::AbstractArray,x::AbstractArray, distibution_sample::compact_mnm_dist , group::Int64 = -1)
    if length(distibution_sample.α) == 0
        return
    end
    @inbounds for i in eachindex(r)
        r[i] = 0.0
        @inbounds for j=1:size(x,1)
            if x[j,i] > 0
                r[i] += distibution_sample.α[Int64.(x[j,i])]
            end
        end
    end
end


text:
export RangeSlider

struct RangeSlider
    range::AbstractRange
    left::Number
    right::Number
    show_value::Bool
end

RangeSlider(range::AbstractRange; left=first(range), right=last(range), show_value=false) =
left > right ? error("Left value cannot be bigger than right") : RangeSlider(range, left, right, show_value)

"""
A `RangeSlider` is a two thumb slider which returns the range between the two thumbs.

If you set `show_value` to `false`, the slider won't display its value.

By default `show_value` is set to `true`

Use `left` and `right` to set default values for the slider thumbs.

## Examples
`@bind range RangeSlider(1:100)`

`@bind range RangeSlider(1:0.1:100, show_value=false)`

`@bind range RangeSlider(1:0.1:100, left=25, right=75)`
"""

function show(io::IO, ::MIME"text/html", slider::RangeSlider)
    js = read(joinpath(PKG_ROOT_DIR, "assets", "rangeslider.js"), String)
    css = read(joinpath(PKG_ROOT_DIR, "assets", "rangeslider.css"), String)

    result = """
    <multi-range-slider>

            <input type="range" id="input-left" step="$(step(slider.range))"
            min="$(first(slider.range))" max="$(last(slider.range))"
            value="$(slider.left)">
            <input type="range" id="input-right" step="$(step(slider.range))"
            min="$(first(slider.range))" max="$(last(slider.range))"
            value="$(slider.right)">

        $(slider.show_value ? """<output id=\"slider-output\"></output>""" : "")

    </multi-range-slider>
    """

    print(io, result)
end

get(slider::RangeSlider) = collect(slider.left:step(slider.range):slider.right)


text:
module TestLGMM

using Test
import TreeParzen: LogGMM

# LGMM1_lpdf dimenstion mismatch
@test_throws DimensionMismatch LogGMM.LGMM1_lpdf(
    [[1.0 2.0] [3.0 4.0]], [1.0, 0.0], [0.0, 1.0, 2.0], [10.0]
)

end
true


text:
# This file is for functions that explictly mess with the dimensions of a NameDimsArray

function Base.dropdims(nda::NamedDimsArray; dims)
    numerical_dims = dim(nda, dims)
    data = dropdims(parent(nda); dims=numerical_dims)
    L = remaining_dimnames_after_dropping(dimnames(nda), numerical_dims)
    return NamedDimsArray{L}(data)
end

function Base.permutedims(nda::NamedDimsArray{L}, perm) where {L}
    numerical_perm = dim(nda, perm)
    new_names = permute_dimnames(L, numerical_perm)

    return NamedDimsArray{new_names}(permutedims(parent(nda), numerical_perm))
end

function Base.selectdim(nda::NamedDimsArray, s::Symbol, i)
    return selectdim(nda, dim(nda, s), i)
end

for f in (
    :(Base.transpose),
    :(Base.adjoint),
    :(Base.permutedims),
    :(LinearAlgebra.pinv)
)
    @eval function $f(nda::NamedDimsArray{L}) where (L)
        data = $f(parent(nda))
        new_names = if ndims(nda) == 1 # vector input
            (:_, first(L))
        elseif ndims(data) == 1 # vector output
            (last(L),)
        else
            (last(L), first(L))
        end
        return NamedDimsArray{new_names}(data)
    end
end

# reshape
# For now we only implement the version that drops dimension names
Base.reshape(x::NamedDimsArray, d::Vararg{Union{Colon, Int}}) = reshape(parent(x), d)


text:
using Compose

compose(context(),
    ngon(0.15, 0.15, 0.08, 5),
    star(0.35, 0.15, 0.08, 5, 0.3),
    xgon(0.55, 0.15, 0.08, 5, 0.3)
) |> SVG("polygon_forms.svg", 3inch, 3inch)


text:
export submodules, minimal_submodules, maximal_submodules, composition_series,
       composition_factors, meataxe

################################################################################
#
#  Tools for MeatAxe
#
################################################################################

function pivot(M::MatElem, i::Int)
  ind = 1
  while ind <= ncols(M) && iszero(M[i, ind])
    ind += 1
  end
  if ind > ncols(M)
    error("Zero row! Can't find a pivot")
  end
  return ind
end
#
# Given a matrix $M$ in echelon form and a vector, it returns
# the vector reduced with respect to $M$
#
function cleanvect(M::T, v::T) where {T}
  @assert nrows(v)==1
  w = deepcopy(v)
  if iszero(v)
    return w
  end
  nr = nrows(M)
  nc = ncols(M)
  for i=1:nrows(M)
    ind = 1
    while ind <= nc && iszero(M[i,ind])
      ind += 1
    end
    if ind > nc
      continue
    end
    if iszero(w[1,ind])
      continue
    end
    mult = divexact(w[1, ind], M[i, ind])
    w[1, ind] = base_ring(M)(0)
    for k = ind+1:ncols(M)
      w[1, k]-= mult*M[i, k]
    end
  end
  return w

end

function cleanvect!(M::T, v::T) where T
  if v == M
    zero!(v)
    return v
  end
  @assert nrows(v) == 1
  w = v
  if iszero(v)
    return w
  end
  R = base_ring(M)
  for i=1:nrows(M)
    if iszero_row(M,i)
      continue
    end
    ind=1
    Miind = M[i, ind]
    while iszero(Miind)
      ind+=1
      Miind = M[i, ind]
    end
    w1ind = w[1, ind]
    if iszero(w1ind)
      continue
    end
    mult=divexact(w1ind, Miind)
    w[1, ind] = R(0)
    for k=ind+1:ncols(M)
      c = M[i, k]
      if !iszero(c)
        w[1,k] -= mult * c
      end
    end
  end
  return v
end

function reduce_mod_rref(M::T, v::T) where T
  @assert nrows(v)==1
  v1 = deepcopy(v)
  reduce_mod_rref!(M, v1)
  return v1
end

function reduce_mod_rref!(M::T, w::T) where {T}
  @assert nrows(w)==1
  if iszero(w)
    return w
  end
  ind =1
  R = base_ring(M)
  for i=1:nrows(M)
    if iszero_row(M, i)
      break
    end
    while iszero(M[i, ind])
      ind += 1
    end
    if iszero(w[1, ind])
      continue
    end
    mult = divexact(w[1,ind], M[i,ind])
    w[1,ind] = zero(R)
    for k = ind+1:ncols(M)
      w[1,k] -= mult*M[i,k]
    end
  end
  return nothing
end

@doc Markdown.doc"""
    closure(C::T, G::Array{T,1}) where T <: MatElem

Given a matrix C representing a subspace of K^n and a list of matrices G representing endomorphisms of K^n,
the function returns a matrix representing the closure of the subspace under the action, i.e. the smallest
subspace of K^n invariant under the endomorphisms.
"""
function closure(C::T, G::Array{T,1}) where {T}
  if nrows(C) != 1
    rref!(C)
  else
    # Do the rref by hand
    for k in 1:ncols(C)
      c = C[1, k]
      if !iszero(c)
        if isone(c)
          break
        end
        C[1, k] = one(base_ring(C))
        cinv = inv(c)
        for j in (k + 1):ncols(C)
          C[1, j] = cinv * C[1, j]
        end
        break
      end
    end
  end
  i=1
  nc = ncols(C)
  while i <= nrows(C)
    w=view(C, i:i, 1:ncols(C))
    for j=1:length(G)
      res = cleanvect!(C, w*G[j]) # cleanvect(C, w*G[j]) but we
                                  # can do it inplace since w*G[j] is fresh
      if !iszero(res)
        C = vcat(C, res)
        if nrows(C) == nc
          i = ncols(C)+1
          break
        end
      end
    end
    i+=1
  end

  if nrows(C) == 1
    r = 0
    for k in 1:ncols(C)
      c = C[1, k]
      if !iszero(c)
        r = 1
        if isone(c)
          break
        end
        C[1, k] = one(base_ring(C))
        cinv = inv(c)
        for j in (k + 1):ncols(C)
          C[1, j] = cinv * C[1, j]
        end
        break
      end
    end
  else
    r = rref!(C)
  end

  if r != nrows(C)
    C = sub(C, 1:r, 1:ncols(C))
  end
  return C
end

#
#  Function to obtain the action of G on the quotient and on the submodule
#
function clean_and_quotient(M::T, N::T, pivotindex::Set{Int}) where {T}
  coeff = zero_matrix(base_ring(M), nrows(N), nrows(M))
  for i=1:nrows(N)
    for j=1:nrows(M)
      if iszero_row(M,j)
        continue
      end
      ind=1
      while iszero(M[j,ind])
        ind+=1
      end
      coeff[i,j]=divexact(N[i,ind], M[j,ind])
      for s=1:ncols(N)
        N[i,s]-=coeff[i,j]*M[j,s]
      end
    end
  end
  vec= zero_matrix(base_ring(M), nrows(N), ncols(M)-length(pivotindex))
  for i=1:nrows(N)
    pos=0
    for s=1:ncols(M)
      if !(s in pivotindex)
        pos+=1
        vec[i,pos]=N[i,s]
      end
    end
  end
  return coeff, vec
end

function clean_and_quotient_quo(M::T, N::T, pivotindex::Set{Int}) where {T}
  for i=1:nrows(N)
    for j=1:nrows(M)
      if iszero_row(M, j)
        continue
      end
      ind=1
      while iszero(M[j,ind])
        ind+=1
      end
      r = divexact(N[i,ind], M[j,ind])
      N[i, ind] = 0
      for s=ind+1:ncols(N)
        N[i,s] -= r*M[j,s]
      end
    end
  end
  vec = zero_matrix(base_ring(M), nrows(N), ncols(M)-length(pivotindex))
  for i=1:nrows(N)
    pos=0
    for s=1:ncols(M)
      if !(s in pivotindex)
        pos+=1
        vec[i,pos]=N[i,s]
      end
    end
  end
  return vec
end

#  Restriction of the action to the submodule generated by C and the quotient
function __split(C::T, G::Vector{T}) where T <: MatElem
  action_sub = Vector{T}(undef, length(G))
  action_quo = Vector{T}(undef, length(G))
  pivotindex = Set{Int}()
  non_zero_rows = Int[]
  for i = 1:nrows(C)
    if iszero_row(C, i)
      continue
    end
    push!(non_zero_rows, i)
    ind = 1
    while iszero(C[i, ind])
      ind += 1
    end
    push!(pivotindex, ind)
  end
  N = zero_matrix(base_ring(C), ncols(C), ncols(C))
  ind = 1
  for i in non_zero_rows
    _copy_matrix_into_matrix(N, ind, 1, view(C, i:i, 1:ncols(C)))
    ind += 1
  end
  for j = 1:ncols(C)
    if j in pivotindex
      continue
    end
    N[ind, j] = 1
    ind += 1
  end
  Ninv = inv(N)
  for i = 1:length(G)
    Mat = N*G[i]*Ninv
    action_sub[i] = view(Mat, 1:length(pivotindex), 1:length(pivotindex))
    action_quo[i] = view(Mat, length(pivotindex)+1:ncols(C), length(pivotindex)+1:ncols(C))
  end
  return ModAlgAss(action_sub), ModAlgAss(action_quo), pivotindex
end

#  Restriction of the action to the submodule generated by C
function _actsub(C::T, G::Vector{T}) where {T <: MatElem}
  esub = Vector{T}(undef, length(G))
  pivotindex = Set{Int}()
  for i=1:nrows(C)
    ind = 1
    while iszero(C[i, ind])
      ind += 1
    end
    push!(pivotindex, ind)
  end
  for a=1:length(G)
    subm, vec = clean_and_quotient(C, C*G[a], pivotindex)
    esub[a] = subm
  end
  return ModAlgAss(esub)
end

#  Restriction of the action to the quotient by the submodule generated by C
function _actquo(C::T, G::Vector{T}) where {T <: MatElem}
  equot = Vector{T}(undef, length(G))
  pivotindex = Set{Int}()
  for i=1:nrows(C)
    ind = 1
    while iszero(C[i,ind])
      ind += 1
    end
    push!(pivotindex,ind)
  end
  for a=1:length(G)
    s = zero_matrix(base_ring(C), ncols(G[1]) - length(pivotindex), ncols(G[1]) - length(pivotindex))
    pos = 0
    for i=1:nrows(G[1])
      if !(i in pivotindex)
        vec = clean_and_quotient_quo(C, sub(G[a], i:i, 1:nrows(G[1])), pivotindex)
        for j=1:ncols(vec)
          s[i - pos, j] = vec[1, j]
        end
      else
        pos += 1
      end
    end
    equot[a] = s
  end
  return ModAlgAss(equot), pivotindex
end

#  Function that determine if two modules are isomorphic, provided that the first is irreducible
function isisomorphic(M::ModAlgAss{S, T, V}, N::ModAlgAss{S, T, V}) where {S, T, V}
  @assert M.isirreducible == 1
  @assert base_ring(M) == base_ring(N)
  @assert length(M.action) == length(N.action)
  if dimension(M) != dimension(N)
    return false
  end

  if M.dimension==1
    return M.action==N.action
  end

  K = base_ring(M)
  Kx, x = PolynomialRing(K, "x", cached=false)

  if length(M.action) == 1
    f = charpoly(Kx, M.action[1])
    g = charpoly(Kx, N.action[1])
    if f==g
      return true
    else
      return false
    end
  end
  rel = _relations(M,N)
  return iszero(rel[N.dimension, N.dimension])

end

function dual_space(M::ModAlgAss{S, T, V}) where {S, T, V}
  G = T[transpose(g) for g in M.action]
  return ModAlgAss(G)
end

function _subst(f::Nemo.PolyElem{T}, a::MatElem{T}) where {T <: Nemo.RingElement}
   n = degree(f)
   if n < 0
      return similar(a)
   elseif n == 0
      return coeff(f, 0) * identity_matrix(base_ring(a), nrows(a))
   elseif n == 1
      return coeff(f, 0) * identity_matrix(base_ring(a), nrows(a)) + coeff(f, 1)*a
   end
   d1 = isqrt(n)
   d = div(n, d1)
   A = powers(a, d)
   s = coeff(f, d1*d)*A[1]
   for j = 1:min(n - d1*d, d - 1)
      c = coeff(f, d1*d + j)
      if !iszero(c)
         s += c*A[j + 1]
      end
   end
   for i = 1:d1
      s *= A[d + 1]
      s += coeff(f, (d1 - i)*d)*A[1]
      for j = 1:min(n - (d1 - i)*d, d - 1)
         c = coeff(f, (d1 - i)*d + j)
         if !iszero(c)
            s += c*A[j + 1]
         end
      end
   end
   return s
end

#################################################################
#
#  MeatAxe
#
#################################################################

@doc Markdown.doc"""
    meataxe(M::ModAlgAss) -> Bool, MatElem

Given module M, returns true if the module is irreducible (and the identity matrix) and false if the space is reducible, togheter with a basis of a submodule
"""
function meataxe(M::ModAlgAss{S, T, V}) where {S, T, V}

  K = base_ring(M)
  Kx, x = PolynomialRing(K, "x", cached = false)
  n = dimension(M)
  @assert n > 0
  H = M.action
  if n == 1
    M.isirreducible = 1
    return true, identity_matrix(K, n)
  end

  G = T[x for x in H if !iszero(x)]
  @assert typeof(G) == typeof(H)

  if isempty(G)
    return false, matrix(K, 1, n, V[one(K) for i = 1:n])
  end

  if isone(length(G))
    A = G[1]
    poly = minpoly(Kx, A)
    sq = factor_squarefree(poly)
    lf = factor(first(keys(sq.fac)))
    t = first(keys(lf.fac))
    if degree(t)==n
      M.isirreducible = 1
      return true, identity_matrix(K, n)
    else
      N = _subst(t, A)
      null, kern = kernel(N, side = :left)
      B = closure(sub(kern, 1:1, 1:n), G)
      return false, B
    end
  end

  #  Adding generators to obtain randomness
  Gt = T[transpose(x) for x in G]
  cnt = 0
  while true
    cnt += 1
    if cnt > 1000
      error("Too many attempts")
    end
    # At every step, we add a generator to the group.
    new_gen = G[rand(1:length(G))]*G[rand(1:length(G))]
    while iszero(new_gen)
      mul!(new_gen, G[rand(1:length(G))], G[rand(1:length(G))])
    end
    push!(G, new_gen)

    # Choose a random combination of the generators of G
    A = zero_matrix(K, n, n)
    for i = 1:length(G)
      add!(A, A, rand(K)*G[i])
    end

    # Compute the characteristic polynomial and, for irreducible factor f, try the Norton test
    poly = minpoly(Kx, A)
    sqfpart = keys(factor_squarefree(poly).fac)
    for el in sqfpart
      sq = el
      i = 1
      while !isone(sq)
        f = gcd(powmod(x, order(K)^i, sq)-x,sq)
        sq = divexact(sq, f)
        lf = factor(f)
        for t in keys(lf.fac)
          N = _subst(t, A)
          a, kern = kernel(N, side = :left)
          @assert a > 0
          #  Norton test
          B = closure(sub(kern, 1:1, 1:n), M.action)
          if nrows(B) != n
            M.isirreducible= 2
            return false, B
          end
          aa, kernt = kernel(transpose(N), side = :left)
          @assert aa == a
          Bt = closure(sub(kernt, 1:1, 1:n), Gt)
          if nrows(Bt) != n
            aa, Btnu = kernel(Bt)
            subst = transpose(Btnu)
            #@assert nrows(subst)==nrows(closure(subst,G))
            M.isirreducible = 2
            return false, subst
          end
          if degree(t) == a
            # f is a good factor, irreducibility!
            M.isirreducible = 1
            return true, identity_matrix(K, n)
          end
        end
        i += 1
      end
    end
  end
end

################################################################################
#
#  Composition series
#
################################################################################

@doc Markdown.doc"""
    composition_series(M::ModAlgAss) -> Array{MatElem,1}

Given a Fq[G]-module M, it returns a composition series for M, i.e. a
sequence of submodules such that the quotient of two consecutive element is irreducible.
"""
function composition_series(M::ModAlgAss{S, T, V}) where {S, T, V}

  if M.isirreducible == 1 || M.dimension == 1
    return [identity_matrix(base_ring(M.action[1]), dimension(M))]
  end

  bool, C = meataxe(M)
  #
  #  If the module is irreducible, we return a basis of the space
  #
  if bool
    return [identity_matrix(base_ring(M.action[1]), dimension(M))]
  end
  #
  #  The module is reducible, so we call the algorithm on the quotient and on the subgroup
  #
  G = M.action
  K = M.base_ring

  rref!(C)

  esub,equot,pivotindex=__split(C,G)
  sub_list = composition_series(esub)
  quot_list = composition_series(equot)
  #
  #  Now, we have to write the submodules of the quotient and of the submodule in terms of our basis
  #
  list=Vector{T}(undef, length(sub_list)+length(quot_list))
  for i=1:length(sub_list)
    list[i]=sub_list[i]*C
  end
  for z=1:length(quot_list)
    s=zero_matrix(K,nrows(quot_list[z]), ncols(C))
    for i=1:nrows(quot_list[z])
      pos=0
      for j=1:ncols(C)
        if j in pivotindex
          pos+=1
        else
          s[i,j]=quot_list[z][i,j-pos]
        end
      end
    end
    list[length(sub_list)+z]=vcat(C,s)
  end
  return list
end

################################################################################
#
#  Composition factors
#
################################################################################

function _composition_factors_cyclic(M::ModAlgAss{S, T, V}; dimension::Int = -1) where {S, T, V}
  @assert length(M.action) == 1
  mat = M.action[1]
  pols, basis_transf, gens = _rational_canonical_form_setup(mat)
  factors, gens_polys_mults = refine_for_jordan(pols, gens, mat)
  mats = Tuple{T, Int}[]
  for i = 1:length(factors)
    for j = 1:length(gens_polys_mults)
      el = gens_polys_mults[j]
      if el[2] != factors[i]
        continue
      end
      if dimension != -1 && degree(factors[i]) != dimension
        continue
      end
      push!(mats, (jordan_block(factors[i], 1), el[3]))
    end
  end
  res = Tuple{typeof(M), Int}[]
  done = falses(length(mats))
  for i = 1:length(mats)
    if done[i]
      continue
    end
    done[i] = false
    m = mats[i][2]
    for j = i+1:length(mats)
      if mats[i][1] == mats[j][1]
        done[j] = true
        m += mats[j][2]
      end
    end
    AA = ModAlgAss(T[mats[i][1]])
    AA.isirreducible = 1
    push!(res, (AA, m))
  end
  return res
end

@doc Markdown.doc"""
    composition_factors(M::ModAlgAss)

Given a Fq[G]-module M, it returns, up to isomorphism, the composition factors of M with their multiplicity,
i.e. the isomorphism classes of modules appearing in a composition series of M
"""
function composition_factors(M::ModAlgAss{S, T, V}; dimension::Int=-1) where {S, T, V}
  if M.isirreducible == 1 || M.dimension == 1
    if dimension != -1
      if M.dimension == dimension
        return Tuple{typeof(M), Int}[(M,1)]
      else
        return Tuple{typeof(M), Int}[]
      end
    else
      return Tuple{typeof(M), Int}[(M,1)]
    end
  end

  if isone(length(M.action))
    return _composition_factors_cyclic(M, dimension = dimension)
  end

  K = base_ring(M)
  bool, C = meataxe(M)
  #  If the module is irreducible, we just return a basis of the space
  if bool
    if dimension != -1
      if M.dimension==dimension
        return Tuple{typeof(M), Int}[(M,1)]
      else
        return Tuple{typeof(M), Int}[]
      end
    else
      return Tuple{typeof(M), Int}[(M,1)]
    end
  end
  G = M.action

  #  The module is reducible, so we call the algorithm on the quotient and on the subgroup
  rref!(C)
  sub, quot, pivotindex = __split(C, G)
  sub_list = composition_factors(sub, dimension = dimension)
  quot_list = composition_factors(quot, dimension = dimension)
  #  Now, we check if the factors are isomorphic
  done = falses(length(quot_list))
  for i = 1:length(sub_list)
    for j = 1:length(quot_list)
      if !done[j] && isisomorphic(sub_list[i][1], quot_list[j][1])
        sub_list[i]=(sub_list[i][1], sub_list[i][2]+quot_list[j][2])
        done[j] = true
        break
      end
    end
  end
  for j = 1:length(done)
    if !done[j]
      push!(sub_list, quot_list[j])
    end
  end
  return sub_list

end

function eigenspace_as_matrix(M::MatElem{T}, lambda::T) where T <: FieldElem
  N = sub(M, 1:nrows(M), 1:ncols(M))
  for i = 1:ncols(N)
    N[i, i] -= lambda
  end
  d, res = Hecke.left_kernel(N)
  return view(res, 1:d, 1:ncols(res))
end

function _relations_dim_1(M::ModAlgAss{S, T, V}, N::ModAlgAss{S, T, V}) where {S, T, V}
  @assert M.isirreducible == 1
  @assert dimension(M) == 1

  K = base_ring(M)
  G = M.action
  H = N.action
  subs = eigenspace_as_matrix(H[1], G[1][1, 1])
  if iszero(nrows(subs))
    return zero_matrix(K, 0, dimension(N))
  end
  for i = 2:length(H)
    newsubs = eigenspace_as_matrix(H[i], G[i][1, 1])
    if iszero(nrows(newsubs))
      return zero_matrix(K, 0, dimension(N))
    end
    subs = _intersect(subs, newsubs)
    if iszero(nrows(subs))
      return zero_matrix(K, 0, dimension(N))
    end
  end
  return subs

end

function _relations(M::ModAlgAss{S, T, V}, N::ModAlgAss{S, T, V}) where {S, T, V}
  @assert M.isirreducible == 1
  G = M.action
  H = N.action
  K = base_ring(M)
  n = dimension(M)

  B = zero_matrix(K, n, n)
  B[1,1] = K(1)
  sys = zero_matrix(K,2*dimension(N),dimension(N))
  matrices = T[identity_matrix(K, dimension(N))]

  X = zero_matrix(K, n, n)
  X[1,1] = K(1)
  i = 1
  k = 1
  while i <= k
    w = sub(B, i:i, 1:n)
    for j = 1:length(G)
      v = w*G[j]
      res = cleanvect(X, v)#reduce_mod_rref(X, v)
      if !iszero(res)
        k += 1
        for s = 1:n
          X[k, s] = res[1, s]
          B[k, s] = v[1, s]
        end
        rref!(X)
        push!(matrices, matrices[i]*H[j])
      else
        fl, x = can_solve(B, v, side = :left)
        @assert fl
        A = sum(T[x[1, q]*matrices[q] for q = 1:k])
        A = sub!(A, A, matrices[i]*H[j])
        for s = 1:N.dimension
          for t = 1:N.dimension
            sys[N.dimension+s, t] = A[t,s]
          end
        end
        rref!(sys)
      end
    end
    if !iszero(sys[N.dimension, N.dimension])
      break
    end
    i += 1
  end
  return view(sys, 1:N.dimension, 1:N.dimension)
end

#Finds the irreducible submodules of N isomorphic to M
#M must be irreducible
function irreducible_submodules(N::ModAlgAss{S, T, V}, M::ModAlgAss{S, T, V}) where {S, T, V}
  @assert M.isirreducible == 1
  K = M.base_ring
  if dimension(M) == 1
    kern = _relations_dim_1(M, N)
    a = rref!(kern)
    #@assert closure(kern, N.action) == kern
    if iszero(a)
      return T[]
    elseif isone(a)
      kern = view(kern, 1:a, 1:ncols(kern))
      return T[kern]
    end
  else
    rel = _relations(M,N)
    if !iszero(rel[N.dimension, N.dimension])
      return T[]
    end
    a, kern = nullspace(rel)
    kern = transpose(kern)
    if a == 1
      return T[closure(kern, N.action)]
    end
  end
  vects = T[sub(kern, i:i, 1:N.dimension) for i=1:a]
  if dimension(M) == 1
    return _submodules_direct_sum_dim_1(vects, N)
  end
  #First sieve: The vectors may be dependent when considering the structure of G-module.
  #So I reduce the number of vectors by looking at the subspace they generate.
  #If there is an inclusion, then I can remove them.
  to_reduce = closure(vects[1], N.action)
  final_vect_list = T[vects[1]]
  for i = 2:length(vects)
    w = cleanvect(to_reduce, vects[i])
    if !iszero(w)
      push!(final_vect_list, vects[i])
      #to_reduce = vcat(vects[i], to_reduce)
      to_reduce = closure(vcat(vects[i], to_reduce), N.action)
    end
  end
  if isone(length(final_vect_list))
    return T[closure(final_vect_list[1], N.action)]
  end
  # Now, I have a list of generators as G-modules.
  # Therefore I need to consider combinations as G-modules
  return _submodules_direct_sum(final_vect_list, N)

end

function _all_combinations(M::MatElem{T}) where T
  K = base_ring(M)
  els = collect(x for x in K)
  @assert fits(Int, fmpz(length(els))^nrows(M))
  res = Vector{typeof(M)}(undef, length(els)^nrows(M))
  ind = 1
  m = zero_matrix(K, 1, nrows(M))
  for i in CartesianIndices(Tuple(1:length(els) for i in 1:nrows(M)))
    for j = 1:nrows(M)
      m[1, j] = els[i[j]]
    end
    res[ind] = m * M
    ind += 1
  end
  return res
end

function _isclosed_dim1(C, G)
  # Do the rref by hand
  for k in 1:ncols(C)
    c = C[1, k]
    if !iszero(c)
      if isone(c)
        break
      end
      C[1, k] = one(base_ring(C))
      cinv = inv(c)
      for j in (k + 1):ncols(C)
        C[1, j] = cinv * C[1, j]
      end
      break
    end
  end
  n = similar(C)
  for i = 1:length(G)
    n = mul!(n, C, G[i])
    cleanvect!(C, n)
    if !iszero(n)
      return false
    end
  end
  return true
end

function _submodules_direct_sum_dim_1(gens::Vector{T}, N::ModAlgAss{S, T, V}) where {S, T, V}
  K = base_ring(N)
  #The module I am working with is the direct sum of the submodules of N generated by the vectors in gens
  #I list all the vectors of the subspaces
  all_combinations = Vector{T}[_all_combinations(x) for x in gens]
  res = Vector{T}()
  for i = 1:length(gens)-1
    #I have to list all the elements that have 1 in the first component and all the possible elements in the other.
    it = CartesianIndices(Tuple(1:length(all_combinations[j]) for j = i+1:length(gens)))
    non_zero_el = all_combinations[i][1]
    ind = 2
    while iszero(non_zero_el)
      non_zero_el = all_combinations[i][ind]
      ind += 1
    end
    for I in it
      m = non_zero_el
      for s = 1:length(I)
        m += all_combinations[i+s][I[s]]
      end
      push!(res, m)
    end
  end
  non_zero_el = all_combinations[length(gens)][1]
  ind = 2
  while iszero(non_zero_el)
    non_zero_el = all_combinations[length(gens)][ind]
    ind += 1
  end
  push!(res, non_zero_el)
  return res
end

function _submodules_direct_sum(gens::Vector{T}, N::ModAlgAss{S, T, V}) where {S, T, V}
  K = base_ring(N)
  #The module I am working with is the direct sum of the submodules of N generated by the vectors in gens
  closures = T[closure(x, N.action) for x in gens]
  #I list all the vectors of the subspaces
  all_combinations = Vector{T}[_all_combinations(x) for x in closures]
  res = Vector{T}()
  for i = 1:length(gens)-1
    #I have to list all the elements that have 1 in the first component and all the possible elements in the other.
    it = CartesianIndices(Tuple(1:length(all_combinations[j]) for j = i+1:length(gens)))
    non_zero_el = all_combinations[i][1]
    ind = 2
    while iszero(non_zero_el)
      non_zero_el = all_combinations[i][ind]
      ind += 1
    end
    for I in it
      m = non_zero_el
      for s = 1:length(I)
        m += all_combinations[i+s][I[s]]
      end
      candidate = closure(m, N.action)
      if nrows(candidate) == nrows(closures[1])
        push!(res, candidate)
      end
    end
  end
  non_zero_el = all_combinations[length(gens)][1]
  ind = 2
  while iszero(non_zero_el)
    non_zero_el = all_combinations[length(gens)][ind]
    ind += 1
  end
  push!(res, closure(non_zero_el, N.action))
  return res
end

function _minimal_submodules(M::ModAlgAss{S, T, V}, dim::Int=M.dimension+1, lf = Tuple{ModAlgAss{S, T, V}, Int}[]) where {S, T, V}
  K = base_ring(M)
  n = dimension(M)
  if isone(M.isirreducible)
    if dim >= n
      return Tuple{T, ModAlgAss{S, T, V}}[(identity_matrix(K, n), M)]
    else
      return Tuple{T, ModAlgAss{S, T, V}}[]
    end
  end
  if isempty(lf)
    lf = composition_factors(M)
  end
  if meataxe(M)[1]
    if dim >= n
      return Tuple{T, ModAlgAss{S, T, V}}[(identity_matrix(K, n), M)]
    else
      return Tuple{T, ModAlgAss{S, T, V}}[]
    end
  end
  if dim != n+1
    lf = Tuple{ModAlgAss{S, T, V}, Int}[x for x in lf if x[1].dimension == dim]
  end
  list = Tuple{T, ModAlgAss{S, T, V}}[]
  for x in lf
    irr_subs = irreducible_submodules(M, x[1])
    for y in irr_subs
      push!(list, (y, x[1]))
    end
  end
  return list
end

@doc Markdown.doc"""
    minimal_submodules(M::ModAlgAss)

Given a Fq[G]-module M, it returns all the minimal submodules of M
"""
function minimal_submodules(M::ModAlgAss{S, T, V}, dim::Int=M.dimension+1, lf = Tuple{ModAlgAss{S, T, V}, Int}[]) where {S, T, V}
  return T[x[1] for x in _minimal_submodules(M, dim, lf)]
end

@doc Markdown.doc"""
    maximal_submodules(M::ModAlgAss)

Given a $G$-module $M$, it returns all the maximal submodules of M
"""
function maximal_submodules(M::ModAlgAss{S, T, V}, index::Int=M.dimension, lf = Tuple{ModAlgAss{S, T, V}, Int}[]) where {S, T, V}

  M_dual = dual_space(M)
  minlist = minimal_submodules(M_dual, index+1, lf)
  maxlist = Array{T, 1}(undef, length(minlist))
  for j=1:length(minlist)
    maxlist[j]=transpose(nullspace(minlist[j])[2])
  end
  return maxlist

end

@doc Markdown.doc"""
    submodules(M::ModAlgAss)

Given a $G$-module $M$, it returns all the submodules of M

"""
function submodules(M::ModAlgAss{S, T, V}) where {S, T, V}

  K = base_ring(M)
  list = T[]
  if iszero(dimension(M))
    return list
  end
  if M.dimension == 1
    return T[zero_matrix(K, 1, 1), identity_matrix(K, 1)]
  end
  lf = composition_factors(M)
  minlist = minimal_submodules(M, M.dimension+1, lf)
  for x in minlist
    rref!(x)
    N, pivotindex = _actquo(x, M.action)
    ls = submodules(N)
    for a in ls
      s=zero_matrix(K,nrows(a), M.dimension)
      for t=1:nrows(a)
        pos=0
        for j=1:M.dimension
          if j in pivotindex
            pos+=1
          else
            s[t,j]=a[t,j-pos]
          end
        end
      end
      push!(list,vcat(x,s))
    end
  end
  for i = 1:length(list)
    x = list[i]
    rk = rref!(x)
    if rk < nrows(x)
      nx = sub(x, 1:rk, 1:ncols(x))
      list[i] = nx
    end
  end
  push!(list, zero_matrix(K, 0, M.dimension))
  push!(list, identity_matrix(K, M.dimension))
  append!(list, minlist)
  res = T[list[1]]
  for i = 2:length(list)
    x = list[i]
    found = false
    for j = 1:length(res)
      if x == res[j]
        found = true
        break
      end
    end
    if !found
      push!(res, x)
    end
  end
  return res

end

@doc Markdown.doc"""
    submodules(M::ModAlgAss, index::Int)

Given a $G$-module $M$, it returns all the submodules of M of index q^index, where q is the order of the field
"""
function submodules(M::ModAlgAss{S, T, V}, index::Int; comp_factors = Tuple{ModAlgAss{S, T, V}, Int}[]) where {S, T, V}
  K = base_ring(M)
  if index == M.dimension
    return T[zero_matrix(K, 1, M.dimension)]
  end
  list = T[]
  list_hash = UInt64[]
  if index >= div(M.dimension, 2)
    if index == M.dimension-1
      if isempty(comp_factors)
        lf = composition_factors(M, dimension = 1)
      else
        lf = comp_factors
      end
      list = minimal_submodules(M, 1, lf)
      return list
    end
    if isempty(comp_factors)
      lf=composition_factors(M)
    else
      lf=comp_factors
    end
    if length(lf) == 1 && dimension(lf[1][1]) == 1
      return _submodules_primary(M, M.dimension-index, lf[1][1])
    end
    for i=1:M.dimension-index-1
      minlist = _minimal_submodules(M, i, lf)
      for X in minlist
        x, Sub = X
        N, pivotindex = _actquo(x, M.action)
        #  Recover the composition factors of the quotient
        lf1 = Tuple{ModAlgAss{S, T, V}, Int}[]
        found = false
        for j = 1:length(lf)
          if !found && isisomorphic(lf[j][1], Sub)
            if !isone(lf[j][2])
              push!(lf1, (lf[j][1], lf[j][2]-1))
            end
            found = true
          else
            push!(lf1, lf[j])
          end
        end
        #  Recursively ask for submodules and write their bases in terms of the given set of generators
        ls = submodules(N, index, comp_factors = lf1)
        length_list = length(list)
        for a in ls
          s = _lift(a, x, pivotindex)
          rref!(s)
          k = hash(s)
          found = false
          for j = 1:length_list
            if k == list_hash[j] && list[j] == s
              found = true
              break
            end
          end
          if !found
            push!(list, s)
            push!(list_hash, k)
          end
        end
      end
    end
    append!(list, minimal_submodules(M, M.dimension-index, lf))
  else
    #  Duality
    M_dual=dual_space(M)
    dlist=submodules(M_dual, M.dimension-index)
    list=T[transpose(nullspace(x)[2]) for x in dlist]
  end
  return list

end

function _lift(a::T, x::T, pivotindex::Set{Int}) where T
  s = zero_matrix(base_ring(a), nrows(a)+nrows(x), ncols(x))
  for t = 1:nrows(a)
    pos = 0
    #Using that a in in rref
    for j = 1:ncols(x)
      if j in pivotindex
        pos += 1
      elseif j >= t
        s[t, j] = a[t, j-pos]
      end
    end
  end
  for t = nrows(a)+1:nrows(s)
    for j = 1:ncols(s)
      s[t, j] = x[t-nrows(a), j]
    end
  end
  return s
end

function _submodules_primary(M::ModAlgAss{S, T, V}, dimension::Int, composition_factor::ModAlgAss{S, T, V}) where {S, T, V}
  list = irreducible_submodules(M, composition_factor)
  K = base_ring(M)
  for i = 2:dimension
    list_new = Set{T}()
    for x in list
      N, pivotindex = _actquo(x, M.action)
      subs = irreducible_submodules(N, composition_factor)
      length_list = length(list_new)
      for a in subs
        s = _lift(a, x, pivotindex)
        rref!(s)
        push!(list_new, s)
      end
    end
    list = collect(list_new)
  end
  return list
end

function powmod(f::Zmodn_poly, e::fmpz, g::Zmodn_poly)
  if fits(Int, e)
    return powmod(f, Int(e), g)
  else
    _e = BigInt()
    z = parent(f)()
    ccall((:fmpz_get_mpz, libflint), Nothing, (Ref{BigInt}, Ref{fmpz}), _e, e)
    ccall((:nmod_poly_powmod_mpz_binexp, libflint), Nothing,
          (Ref{Zmodn_poly}, Ref{Zmodn_poly}, Ref{BigInt}, Ref{Zmodn_poly}),
           z, f, e, g)
    return z
  end
end


text:
#===============================================================================
This file contains an independent implementation of the excellent  line search
algorithm described in:

	"Line Search Algorithms with Guaranteed Sufficient Decrease"
	Jorge J. Moré and David J. Thuente
    ACM Transactions on Mathematical Software (TOMS) (1994)
	doi: 10.1145/192115.192132

Features:
  - Finite termination
  - Wolfe conditions
  - α stays in minstep/maxstep
Rwquires:
  - c1 ≤ c2
  - ϕ is C1
===============================================================================#

struct MoreThuente{T1, T2} <: LineSearch
    ratio::T1
	c::T1
	maxiter::T2
	verbose::Bool
end

#===============================================================================
                            Trial value selection
	Notation:
		αc is the minimizer of the cubic that interpolates fl, ft, gl, gt
		αq is the minimizer of the quadratic that interpolates fl, ft, gl
		αs is the minimizer of the quadratic that interpolates fl, gl, gt
===============================================================================#
#===============================================================================
	CASE 1: a) ft > fl

	Compute αc and αq and set αt⁺ = ifelse(|αc-αl|<|αq-αl, αc, (αq+αc)/2).
	This choice places us closer to αl - the step length that gives us the
	lowest function value so far.
===============================================================================#
αt⁺ = abs(αc - αl) < abs(αq - αl) ? αc : (αq + αc)/2
#===============================================================================
	CASE 2: a) ft ≤ fl
	        b) gt×gl < 0
===============================================================================#
αt⁺ = abs(αc - αt) < abs(αs - αt) ? αc : αs
#===============================================================================
	CASE 3: a) ft ≤ fl
            b) gt×gl ≥ 0
            c) |gt| ≤ |gl|

	αt⁺ :  cubic that interpolates fl, ft, gl, gt may not have a minimizer
===============================================================================#
if cubic tend to inf for alpha inc and is αc > αt?
	αt⁺ = ifelse(abs(αc - αt) < abs(αs - αt)), αc, αs
else
	αt⁺ = αs
end
#===============================================================================
	CASE 4: a) ft ≤ fl
	        b) gt×gl ≥ 0
	        c) |gt| < |gl|

	αt⁺ : minimizer of cubic that interpolates fu, ft, gu, gt
===============================================================================#

function find_steplength(ls::Union{Backtracking, TwoPointQuadratic}, φ::T, d, x,
	                     λ, #
						 ϕ_0::Tf,
						 ∇f_x,
						 dϕ_0=dot(d, ∇f_x);
						 minlength,
						 maxlength) where {T, Tf}
	ratio, c, maxiter, verbose = Tf(ls.ratio), Tf(ls.c), ls.maxiter, ls.verbose

    if verbose
        println("Entering line search with step size: ", λ)
        println("Initial value: ", ϕ_0)
        println("Value at first step: ", φ(λ))
    end

    t = -c*dϕ_0

    iter, α, β = 0, λ, λ # iteration variables
	f_α = φ(α)   # initial function value

	is_solved = isfinite(f_α) && f_α <= ϕ_0 + c*α*t

    while !is_solved && iter <= maxiter
        iter += 1
        β, α, f_α = _solve_pol(ls, φ, ϕ_0, dϕ_0, α, f_α, ratio)
		is_solved = isfinite(f_α) && f_α <= ϕ_0 + c*α*t
    end

	ls_success = iter >= maxiter ? false : true

    if verbose
		!ls_success && println("maxiter exceeded in backtracking")
        println("Exiting line search with step size: ", α)
        println("Exiting line search with value: ", f_α)
    end
    return α, f_α, ls_success
end


text:
using Revise
using Plots
using Basins
using DynamicalSystems

# SHRIMALI, M. D., PRASAD, A., RAMASWAMY, R., & FEUDEL, U. (2008). THE NATURE OF ATTRACTOR BASINS IN MULTISTABLE SYSTEMS. International Journal of Bifurcation and Chaos, 18(06), 1675–1688. doi:10.1142/s0218127408021269
function coupled_logistic_maps(dz,z, p, n)
    θ = z[3]; xn = z[1]; yn = z[2]
    α=3.5; β=0.1; β₂=β*(4/α-1); β₁=0.01; ω = (√5-1)/2
    dz[1]=α*(1 +β₂*cos(2π*θ))*xn*(1-xn)+β₁*(yn-xn)
    dz[2]=α*(1 +β₂*cos(2π*θ))*yn*(1-yn)+β₁*(xn-yn)
    dz[3]=mod(θ+ ω,1)
    return
end

# dummy function to keep the initializator happy
function coupled_logistic_maps_J(J,z0, p, n)
   return
end
ds = DiscreteDynamicalSystem(coupled_logistic_maps,[1., -1., 0.], [] , coupled_logistic_maps_J)
integ  = integrator(ds)

xg=range(0.,0.999999,length=100)
yg=range(0.,0.999999,length=100)

@time bsn=Basins.basins_general(xg, yg, integ; dt=1)
plot(xg,yg,bsn.basin',seriestype=:heatmap)

# Estimation using box counting
ind  = findall(iseven.(bsn.basin) .== true)
basin_test = deepcopy(bsn.basin)
[basin_test[k] =basin_test[k]+1 for k in ind ]
bd = box_counting_dim(xg,yg, basin_test)
α = 2 - bd


text:
function run!{T}(method::SuperlinearDimer, E, dE, x0::Vector{T}, v0::Vector{T})
   # read all the parameters
   @unpack maximum_translation, max_num_rot, trial_angle, trial_trans_step, extrapolate, translation_method,
           tol_trans, tol_rot, maxnumdE, len, precon_prep!, verbose = method
  # @show translation_method
  P=method.precon
   # initialise variables
   x, v = copy(x0), copy(v0)
   nit = 0
   numdE, numE = 0, 0
   log = DimerLog()
   if translation_method == "CG"
      cg_init = true
      direction_old = nothing
      cg_direction = nothing
   elseif translation_method == "LBFGS"
      reset_hessian = true
      s = []
      y = []
      rho = []
      x0, v0 = x, v
      vec0 = x
   else
      error("unknown translation_method")
   end

   # and just start looping
   if verbose >= 2
      @printf(" nit |  n∇E  |   |∇E|_∞    |∇R|_∞    λ        \n")
      @printf("-----|-------|--------------------------------\n")
   end
   P = precon_prep!(P, x)
   v /= sqrt(dot(v, P, v))
   normalize(a) = a / norm(a)
   parallel_vector(a, b) = dot(a, b) * b    # assume b is normalised
   perpendicular_vector(a, b) = a - parallel_vector(a, b)
   rotate_vectors(a, b, phi) =    # do a, b need to be normalised??? CHECK!
         normalize(a * cos(phi) + b * sin(phi)) * norm(a),
         normalize(b * cos(phi) - a * sin(phi)) * norm(b)
   curv = 0.0
   for nit = 0:2*maxnumdE
      # normalise v
      P = precon_prep!(P, x)
      # evaluate gradients, and more stuff
      dE0 = dE(x)
      numdE += 1

      # rotation step

      stoprot = false
      nrot = 0
      dE1e = nothing      # extrapolated guess for dE1
      res_rot = 0.0
      while !stoprot     # ROTATION STEP

          nrot += 1

          if dE1e == nothing
             dE1 = dE(x + len * v)
             numdE += 1
          else
             dE1 = dE1e
          end
          dE2 = dE(x - len * v)
          numdE += 1

          dE1a = copy(dE1)     # dE1a store for later use (does it need to be a copy?)

          Hxv = (dE1 - dE2) / 2.0 / len
          curv = dot(Hxv, v)
          f_rot_A = perpendicular_vector(Hxv, v)   # rotation residual!
          res_rot = norm(f_rot_A, Inf)

          if res_rot <= tol_rot
             stoprot = true
          else                  # single rotation step
             n_A = v
             rot_unit_A = normalize(f_rot_A)

             # get the curvature and its derivative
             c0 = curv
             c0d = dot(dE1 - dE2, rot_unit_A) / len

             # trial rotation
             n_B, rot_unit_B = rotate_vectors(n_A, rot_unit_A, trial_angle)
             v = n_B

             # update gradients
             dE1 = dE(x + len * v)
             dE2 = dE(x - len * v)
             numdE += 2

             dE1b = copy(dE1)

             # get the curvature's derivative
             c1d = dot(dE1 - dE2, rot_unit_B) / len

             # calculate Fourier coefficients
             a1 = c0d * cos(2.0 * trial_angle) - c1d / (2.0 * sin(2.0 * trial_angle))
             b1 = 0.5 * c0d
             a0 = 2 * (c0 - a1)

             # estimate rotation angle
             rotangle = atan(b1 / a1) / 2.0

             # make sure its not the maximum
             cmin = a0 / 2.0 + a1 * cos(2 * rotangle) + b1 * sin(2 * rotangle)
             if c0 < cmin
                rotangle += pi / 2.0
             end

             # rotate into the estimated lowest eigenmode
             n_C, rot_unit_C = rotate_vectors(n_A, rot_unit_A, rotangle)
             v = n_C

             curv = cmin

             # force extrapolation to reduce number of calls
             if extrapolate
                dE1e = sin(trial_angle - rotangle) / sin(trial_angle) * dE1a + sin(rotangle) / sin(trial_angle) * dE1b + (1 - cos(rotangle) - sin(rotangle) * tan(trial_angle / 2.0)) * dE0
             else
                dE1e = nothing
             end

             if !stoprot
                if nrot >= max_num_rot
                   stoprot = true
                end
             end
          end
      end

      # TRANSLATION STEP

      # choose e-vec following or standard dimer direction
      if curv > 0
         f0p = -parallel_vector(-dE0, v)
      else
         f0p = -dE0 - 2.0 * parallel_vector(-dE0, v)
      end
      direction = f0p

      # mix up directions according to CG or LBFGS (or nothing)
      if translation_method == "CG"
         if cg_init
            cg_init = false
            direction_old = direction
            cg_direction = direction
         end
         old_norm = dot(direction_old, direction_old)
         # Polak-Ribiere Conjugate Gradient
         if old_norm != 0.0
            betaPR = dot(direction, (direction - direction_old)) / old_norm
         else
            betaPR = 0.0
         end
         if betaPR < 0.0
            betaPR = 0.0
         end
         cg_direction = direction + cg_direction * betaPR
         direction_old = direction

         direction = cg_direction
      elseif translation_method == "LBFGS"
         error("LBFGS not implemented")
      end

      direction = normalize(direction)

      if curv > 0.0
         xstep = direction * maximum_translation
      else
         dE0t = dE(x + direction * trial_trans_step)
         numdE += 1
         f0tp = -dE0t - 2.0 * parallel_vector(-dE0t, v)
         # f0p: current translation residual (not necessarily search direction!)
         # f0tp: translation residual at trial step
         # F: some kind of slope in search direction???
         # C: curvature in search direction?
         F = dot(f0tp + f0p, direction) / 2.0
         C = dot(f0tp - f0p, direction) / trial_trans_step
         xstep = ( -F / C + trial_trans_step / 2.0 ) * direction
         if norm(xstep) > maximum_translation
            xstep = direction * maximum_translation
         end
      end

      x += xstep

      # translation residual, store history
      res_trans = vecnorm(dE0, Inf)
      push!(log, numE, numdE, res_trans, res_rot)
      if verbose >= 2
         @printf("%4d | %4d  |  %1.2e  %1.2e  %1.2e \n", nit, numdE, res_trans, res_rot, curv)
      end
      if res_trans <= tol_trans
         if verbose >= 1
            println("SuperlinearDimer terminates succesfully after $(nit) iterations")
         end
         return x, v, log
      end
      if numdE > maxnumdE
         if verbose >= 1
            println("SuperlinearDimer terminates unsuccesfully due to numdE >= $(maxnumdE)")
         end
         return x, v, log
      end

   end
   error("why am I here?")
   return x, v, log
end


text:
using Profile
using ProfileView
function myfunc()
           A = rand(200, 200, 400)
           maximum(A)
       end

myfunc()
@profile myfunc()

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
## TODO: revamp these older functions to newer API in IIF/RoME/Caesar (Nov2018)
# standardize a common front-end design with feature tracking and real-time integration
# part of DFG discussion

function triggerPose(x, xprev, Tnow, Tprev,
                    distrule, timerule, yawrule)

  if norm(x[1:2]-xprev[1:2]) >= distrule
    @show Tnow, round(x,digits=2), round(xprev,digits=2)
    @show norm(x[1:2]-xprev[1:2])
    return 1
  end
  if abs(x[3]-xprev[3]) >= yawrule
    @show Tnow, round(x,digits=2), round(xprev,digits=2)
    @show abs(x[3]-xprev[3])
    return 2
  end
  if Tnow-Tprev > timerule
    return 3
  end
  return 0
end

mutable struct GenericInSituSystem{T}
  xprev::Array{Float64,1}
  x::Array{Float64,1}
  dOdo::Dict{Int,Array{Float64,1}}
  FeatAssc::Dict{Int, Dict{Int,Array{Float64,1}}}
  Tprev::Float64
  T0::Float64
  poseid::Int
  wTbk1::Array{Float64,2}
  bk1Tbk::Array{Float64,2}
  lstlaseridx::Int
  trackers::Dict{Int,T}
end

const InSituSystem = GenericInSituSystem{Feature}

function makeInSituSys(x::Array{Float64,1}, bfts0::Array{Float64,2})
  dOdo = Dict{Int,Array{Float64,1}}()
  FeatAssc = Dict{Int, Dict{Int,Array{Float64,1}}}()
  Tprev = 0.0
  T0 = 0.0
  wTbk1 = SE2(x)
  bk1Tbk = SE2(zeros(3))
  poseid = 1
  dOdo[poseid] = [x[1];x[2];x[3];T0;0]
  lstlaseridx = 1
  trackers = initTrackersFrom(bfts0)
  return InSituSystem(
  x,
  x,
  dOdo,
  FeatAssc,
  Tprev,
  T0,
  poseid,
  wTbk1,
  bk1Tbk,
  lstlaseridx,
  trackers
  )
end

function makeGenericInSituSys(x::Array{Float64,1})
  dOdo = Dict{Int,Array{Float64,1}}()
  FeatAssc = Dict{Int, Dict{Int,Array{Float64,1}}}()
  Tprev = 0.0
  T0 = 0.0
  wTbk1 = SE2(x)
  bk1Tbk = SE2(zeros(3))
  poseid = 1
  dOdo[poseid] = [x[1];x[2];x[3];T0;0]
  lstlaseridx = 1
  trackers = Dict{Int,Any}()
  return GenericInSituSystem{Any}(
  x,
  x,
  dOdo,
  FeatAssc,
  Tprev,
  T0,
  poseid,
  wTbk1,
  bk1Tbk,
  lstlaseridx,
  trackers
  )
end

# doesn't work in full general case yet, still requires reset to zero at each new pose
function poseTrigAndAdd!(instSys::GenericInSituSystem, Ts::Float64,
                        distrule::Float64, timerule::Float64, yawrule::Float64;
                        xprev=zeros(3), auxtrig::Bool=false)
  rule = triggerPose(instSys.x, xprev, Ts, instSys.Tprev, distrule, timerule, yawrule)
  if rule != 0 || auxtrig
    instSys.bk1Tbk = SE2(instSys.x)
    n = [instSys.x;[Ts;rule]]
    instSys.poseid +=1
    instSys.dOdo[instSys.poseid] = n
    instSys.wTbk1 = instSys.wTbk1*instSys.bk1Tbk
    instSys.Tprev = Ts
    instSys.x[1] = 0.0; instSys.x[2]=0.0; instSys.x[3]=0.0
    return true
  end
  return false
end

function poseTrigAndAdd!(instSys::InSituSystem, Ts::Float64,
                        distrule::Float64, timerule::Float64, yawrule::Float64;
                        xprev=zeros(3), auxtrig::Bool=false)

  error(" Yeah dehann, do it properly")

end

function processTreeTrackersUpdates!(instSys::InSituSystem, lsrFeats::Dict{Int,LaserFeatures},
                                    Ts::Float64, b1Dxb::Array{Float64,1}, DBG=Union{}; dbgflag=true)
  propAllTrackers!(instSys.trackers, b1Dxb, [0.05;0.05;0.004])
  newlsridx, Ta = getFeatsAtT(lsrFeats, Ts, prev=instSys.lstlaseridx)
  if newlsridx != instSys.lstlaseridx
    # we have a new laser scan available and should update the trackers
    instSys.lstlaseridx = newlsridx
    bfts = lsrFeats[instSys.lstlaseridx].feats
    hardassc = assocMeasWFeats!(instSys.trackers, bfts)
    measUpdateTrackers!(instSys.trackers, hardassc, [0.5;0.05]) # TODO, noise here is not in polar coords!! something wrong
    if dbgflag
      DBG[instSys.lstlaseridx] = deepcopy(instSys.trackers)
    end
  end

  nothing
end

function advOdoByRules(DRS::Array{Float64,2}, lsrFeats::Dict{Int,LaserFeatures};
                        distrule=20.0, timerule=30.0, yawrule=pi/3.0, trkfeats=true)
  # DBG = Dict{Int,Dict{Int,Array{BallTreeDensity,1}}}()
  # DBG = Dict{Int,Dict{Int,Feature}}()
  DBG = Union{}

  bfts0 = lsrFeats[1].feats
  instSys = makeInSituSys(zeros(3), bfts0)
  fdict = Dict{Int,Array{Float64,1}}()
  for f in instSys.trackers
    # if length(f[2].lastz) < 3
    #   error("HEY!! $(f[1])")
    # end
    fdict[f[2].id] = f[2].lastz
  end
  instSys.FeatAssc[instSys.poseid] = fdict
  for i in 1:size(DRS,1)
    # do odo
    dt = DRS[i,1]-instSys.T0
    whlspd, strang = compensateRawDRS(vec(DRS[i,:]))
    bTbm = SE2(instSys.x)
    instSys.x=uteOdomEasy(instSys.x, whlspd, strang, dt)
    bTbp = SE2(instSys.x)

    # propagate, and maybe measure update, feature trackers
    if trkfeats
      #@show keys(instSys.trackers)
      bmTbp = vec(se2vee(inv(bTbm)*bTbp))
      @time processTreeTrackersUpdates!(instSys, lsrFeats, DRS[i,1], bmTbp, DBG,dbgflag=false)
    end
    # subsample motion tracks for constructing the factor graph SLAM solution
    if poseTrigAndAdd!(instSys, DRS[i,1], distrule, timerule, yawrule)
      fdict = Dict{Int,Array{Float64,1}}()
      for f in instSys.trackers
        mpt = vec(Base.mean(getPoints(f[2].bel),2))
        r,b = c2p(mpt)
        fdict[f[2].id] = [r;b;f[2].lastz[3]]
      end
      instSys.FeatAssc[instSys.poseid] = fdict
    end
    # advance previous time point
    instSys.T0 = DRS[i,1]
  end
  return instSys.dOdo, instSys.FeatAssc, DBG
end

# function loadVicPrkDataset(filename::AbstractString="datasets/VicPrk.jld")
#   DRS,GPS,LsrFeats,d,f = jldopen(filename, "r") do file
#     read(file, "DRS")
#     read(file, "GPS")
#     read(file, "LsrFeats")
#     read(file, "d")
#     read(file, "f")
#   end
#   return DRS,GPS,LsrFeats,d,f
# end


text:
using Base.Test
using QuantumOptics

@testset "schroedinger" begin

N = 3
Ncutoff = 2
T = [0.:0.1:1.;]

Ω = [0 2 3;
     2 0 1;
     3 1 0]

ω = [1., 1.2, 1.5]

basis_fock = FockBasis(Ncutoff)
basis = tensor([basis_fock for i=1:N]...)

a = destroy(basis_fock)
at = create(basis_fock)
I = identityoperator(basis_fock)

psi0 = tensor([coherentstate(basis_fock, i%Ncutoff) for i=1:N]...)

# Interaction picture
Hrot = SparseOperator[]
for i=1:N, j=1:N
    if i==j
        continue
    end
    h = embed(basis, [i,j], [a, Ω[i,j]*at])
    push!(Hrot, h)
end
Hrot = sum(Hrot)

# Schroedinger picture
function f(t, psi)
    H = SparseOperator[embed(basis, i, ω[i]*at*a) for i=1:N]
    for i=1:N, j=1:N
        if i==j
            continue
        end
        h = embed(basis, [i,j], [a, exp(1im*(ω[i]-ω[j])*t)*Ω[i,j]*at])
        push!(H, h)
    end
    sum(H)
end

tout, psi_rot_t = timeevolution.schroedinger(T, psi0, Hrot)
tout, psi_t = timeevolution.schroedinger_dynamic(T, psi0, f)

for (i, t) in enumerate(tout)
    R = prod([embed(basis, i, expm(1im*ω[i]*t*full(at*a))) for i=1:N])
    psi_rot = psi_rot_t[i]
    psi = psi_t[i]
    # @test abs(dagger(psi_rot)*R*psi) < 1e-5
    rho = dm(psi)
    rho_rot = dm(psi_rot)
    @test tracedistance(rho_rot, full(R)*rho*dagger(full(R))) < 1e-5
end

function fout(t, psi)
 deepcopy(psi)
end
t_fout, psi_fout = timeevolution.schroedinger(T, psi0, Hrot; fout=fout)
@test t_fout == tout && psi_fout == psi_rot_t

end # testset


text:
function b2c!(wc::AbstractVector, c::AbstractVector, α,
              prev=Array(eltype(wc),length(wc)))
    T = eltype(wc)

    fill!(wc, zero(T))
    desired_order = length(wc) - 1

    m1 = length(c)-1
    for i=-m1:0
        copy!(prev, wc)
        wc[1] = c[-i+1]

        if desired_order >= 1
            wc[2] = (1.0-α*α)*prev[1] + α*prev[2]
        end
        for m=3:desired_order+1
            @inbounds wc[m] = prev[m-1] + α*(prev[m] - wc[m-1])
        end
    end

    wc
end

b2c(c::AbstractVector, order, α) = b2c!(Array(eltype(c), order+1), c, α)


text:
## Skating system for final placement

"""
    skating_combined(dances, results_single_dances)

Calculates the final placement of competitors in a dancing competition based on the skating
rules. `dances` is a vector of strings, containing the names of the dances.
`results_single_dances` is a vector of DataFrames, containing the results for the individual
dances. See [`skating_single_dance`](@ref) for the correct input format.

Outputs a DataFrame with the accumulated skating tableau as strings, a vector containing
tableaus containing the application of Rule 10 and Rule 11, the skating results for the
individual dances as well as a DataFrame with the starter numbers and final places.
"""
function skating_combined(dances, results_single_dances, places, reports)
    rule10_counts, rule10_sums = prepare_sums(places, size(places, 1))

    places[!, :Sum] = sum(eachcol(places[!, Not(:Number)]))
    places[!, :Place] .= 0.0
    places_text = string.(copy(places))
    rule10_text = string.(copy(rule10_counts))
    rule10_text[!, Not(:Number)] .= "-"
    rule10_text[!, :Place] .= "-"

    rule11_table = get_rule11_table(results_single_dances)
    rule11_text = copy(rule10_text)

    current_place = 1
    while current_place <= size(places, 1)
        idx = findall(==(minimum(places[!, :Sum])), places[!, :Sum])
        if length(idx) == 1
            # @info "Minimum sum"
            write_places!(places, places_text, idx, current_place)
            current_place += 1
        else
            # Rule 10
            id = findall(==(maximum(rule10_counts[idx, current_place+1])), rule10_counts[idx, current_place+1])
            if length(id) == 1
                # @info "Minimum better places"
                index = idx[id]
                write_places!(places, places_text, index, current_place)
                rule10_text[idx, current_place+1] .= string.(rule10_counts[idx, current_place+1])
                rule10_text[index, :Place] .= string(current_place)
                current_place += 1
            else
                # @info "Minimum better summed places"
                i = findall(==(minimum(rule10_sums[idx[id], current_place+1])), rule10_sums[idx[id], current_place+1])
                if length(i) == 1
                    index = idx[id[i]]
                    write_places!(places, places_text, index, current_place)
                    rule10_text[idx[id], current_place+1] .= string.(rule10_counts[idx[id], current_place+1]) .* "*" .* "(" .* string.(rule10_sums[idx[id], (current_place+1)]) .* ")"
                    rule10_text[index, :Place] .= string(current_place)
                    current_place += 1
                else
                    # Rule 11
                    # @info "Rule 11"
                    skating_text, skating_result = skating_single_dance(rule11_table[idx[id[i]], :];
                        initial_place = current_place, initial_column = current_place, depth = size(places, 1))
                    if length(i) == 2
                        index = idx[id[i]]
                        write_places!(places, places_text, index, skating_result[!, :Place], skating_text[!, :Place], current_place)
                        rule11_text[index, (current_place+1):end] .= skating_text[!, (current_place+1):end]
                        current_place += 2
                    else
                        j = findall(==(minimum(skating_result[!, :Place])), skating_result[!, :Place])
                        if length(j) == 1
                            index = idx[id[i[j]]]
                            write_places!(places, places_text, index, skating_result[j, :Place], skating_text[j, :Place], current_place)
                            rule11_text[idx[id[i]], (current_place+1):(end-1)] .= skating_text[!, (current_place+1):(end-1)]
                            rule11_text[index, :Place] .= skating_text[j, :Place]
                            current_place += 1
                        else
                            # @info "Rule 11 Tie"
                            # shared place
                            index = idx[id[i[j]]]
                            place_list = range(current_place, length = length(j))
                            place = mean(place_list)
                            write_places!(places, places_text, index, place, current_place)
                            rule11_text[index, (current_place+1):(end-1)] .= skating_text[!, (current_place+1):(end-1)]
                            rule11_text[index, :Place] .= string(place)
                            current_place = maximum(place_list) + 1
                        end
                    end
                end
            end
        end
    end
    return places_text, (; Rule10 = rule10_text, Rule11 = rule11_text), (; zip(Symbol.(dances), reports)...), places[!, [:Number, :Place]]
end

function skating_combined(dances, results_single_dances)
    return skating_combined(dances, results_single_dances, get_places(dances, results_single_dances)...)
end

"""
    get_places(dances, results_single_dances)

Collects the skating results for a vector of dances and result DataFrames.
"""
function get_places(dances, results_single_dances)
    places = DataFrame(Number = results_single_dances[1][!, :Number])
    reports = []
    for (i, df) in enumerate(results_single_dances)
        report, pl = skating_single_dance(df)
        push!(reports, report)
        insertcols!(places, Symbol(dances[i]) => pl[!, :Place])
    end
    return places, reports
end

"""
    get_rule11_table(results_single_dances)

Prepares a DataFrame with all judgements concatenated into a single table for use when
applying Rule 11.
"""
function get_rule11_table(results_single_dances)
    numbers = DataFrame(Number = results_single_dances[1][!, :Number])
    df_cat = hcat([df[!, Not(:Number)] for df in results_single_dances]..., makeunique = true)
    return hcat(numbers, df_cat)
end

"""
    write_places!(places, places_text, index, current_place)
    write_places!(places, places_text, index, place_to_write, current_place)
    write_places!(places, places_text, index, place_to_write, place_to_write_text, current_place)

Write the result of the skating procedure to the DataFrames containing the final calculation.
"""
function write_places!(places, places_text, index, place_to_write, place_to_write_text, current_place)
    places[index, :Place] .= place_to_write
    places[index, :Sum] .= 1000 + current_place
    places_text[index, :Place] .= place_to_write_text
end
function write_places!(places, places_text, index, place_to_write, current_place)
    write_places!(places, places_text, index, place_to_write, string.(place_to_write), current_place)
end
function write_places!(places, places_text, index, current_place)
    write_places!(places, places_text, index, current_place, current_place)
end


text:
using Random
using HDF5 # ProgressMeter
# τ >= 2  period and n sites
periodicSignal(n, τ) = [i % τ ==0 ? 1.0 : 0.0 for i in 0:n-1]

function signal_plus_noise(signal, σ)
    noise = σ .* abs.(randn(length(signal)))
    signalx =  signal .+ noise
    signalx./maximum(signalx)
end
function signal_plus_noise(signal, σ, doble)
    noise = σ .* randn(length(signal))
    signalx =  signal .+ noise
    minsig = minimum(signalx)
    signalx = signalx .- minsig
    signalx ./maximum(signalx)
end

function gen_samples_signal(signal, σ, num_samples)
    set_samples = zeros(length(signal), num_samples)
    for i in 1:num_samples
        set_samples[:, i] = signal_plus_noise(signal, σ)
    end
    set_samples
end

function gen_samples_signal(signal, σ, num_samples, doble)
    set_samples = zeros(length(signal), num_samples)
    for i in 1:num_samples
        set_samples[:, i] = signal_plus_noise(signal, σ, doble)
    end
    set_samples
end

function gen_data_set(all_signals, σ, num_samples)
    set_samples = gen_samples_signal(all_signals[1, :], σ, num_samples)
    labels = fill(1, num_samples)
    for i in 2:size(all_signals)[1]
        set_samplesx = gen_samples_signal(all_signals[i, :], σ, num_samples)
        set_samples = hcat(set_samples, set_samplesx)
        labelsx = fill(i, num_samples)
        labels = vcat(labels,labelsx)
    end
    set_samples, labels
end
function gen_data_set(all_signals, σ, num_samples, doble)
    set_samples = gen_samples_signal(all_signals[1, :], σ, num_samples, doble)
    labels = fill(1, num_samples)
    for i in 2:size(all_signals)[1]
        set_samplesx = gen_samples_signal(all_signals[i, :], σ, num_samples, doble)
        set_samples = hcat(set_samples, set_samplesx)
        labelsx = fill(i, num_samples)
        labels = vcat(labels,labelsx)
    end
    set_samples, labels
end

function gen_signals(n_sites, periods)
    all_signals = zeros(length(periods), n_sites)
    i = 1
    for τ in periods
        all_signals[i, :] = periodicSignal(n_sites, τ)
        i += 1
    end
    all_signals
end

max_stars_repo_name: lazarusA/noisySignals

max_issues_repo_name: lazarusA/noisySignals

max_forks_repo_name: lazarusA/noisySignals


text:
"""
    @fastkw function f(x, y, ...; a=A, b=B, ...)
        ...
    end

Construct a fast, type-stable set of dispatches for `f` such that
"""
macro fastkw(expr)
    (fname, kw_args, std_args, body) = function_def(expr)

    kw_names = map(kw -> esc(kw.args[1]), kw_args)
    kw_defaults = map(kw -> esc(kw.args[2]), kw_args)

    std_names = map(esc, std_args)
    std_defaults = nothing # TODO

    out_exprs = Vector{Any}()

    # First make a series of definitions for any standard arguments with default values
    # (TODO extract and support default values)

    # Next dispatch the final "standard" function to a keyword version.
    push!(out_exprs, quote
        function $(esc(fname))($(std_names)...)
            Base.@_propagate_inbounds_meta
            $(esc(fname))($(vcat((KWSentinal,), std_args)...))
        end
    end)

    # Now populate the keyword arguments one at a time
    kw_types = (KWSentinal, )
    for i = 1:length(kw_args)
        push!(out_exprs, quote
            function $(esc(fname))(kw::Tuple{$(kw_types...)}, $(std_names...))
                Base.@_propagate_inbounds_meta
                $(esc(fname))((kw..., KW{$(kw_names[i])}($(esc(kw_defaults[i])))), $(std_args...))
            end
        end)
        kw_types = (kw_types..., KW{kw_names[i].args[1]})
    end

    # Finally, when all keywords are populated, fill in local variables and given body
    kw_bindings = map(kw_names) do kw_name
        :(KW($(kw_name)) = kw[Val{$(kw_name)}])
    end

    metas = nothing # TODO: extract meta's from beginning of body and insert in correct place.
    main_body = esc(body)

    push!(out_exprs, quote
        function $(esc(fname))(kw::Tuple{$(kw_types...)}, $(std_args...))
            $metas
            $(Expr(:block, kw_bindings...))
            $main_body
        end
    end)

    # Currently don't yet support additional unnamed keyword arguments (TODO)
    # TODO: support :stagedfunction
    # TODO: support where clauses

    return Expr(:block, out_exprs...)
end

function_def(expr) = error("Expected a function definition")

function function_def(expr::Expr)
    if expr.head == :function || expr.head == :(=) || expr.head == :kw
        @assert length(expr.args) == 2
        body = expr.args[2]

        if (expr.args[1] isa Expr) && (expr.args[1].head == :call)
            fname = expr.args[1].args[1]
            args = expr.args[1].args[2:end]
            if (length(args) > 0) && (args[1] isa Expr) && (args[1].head == :parameters)
                kw_args = args[1].args
                std_args = args[2:end]
            else
                kwargs = Any[]
                std_args = args
            end
        else
            error("Expected a function definition")
        end

    else
        error("Expected a function definition")
    end
    return (fname, kw_args, std_args, body)
end


text:
# Special-Ordered-Set type 1

MOI.supports_constraint(o::Optimizer, ::Type{VECTOR}, ::Type{SOS1}) = true

function MOI.add_constraint(o::Optimizer, func::VECTOR, set::SOS1)
    allow_modification(o)

    varrefs = [VarRef(vi.value) for vi in func.variables]
    cr = add_special_ordered_set_type1(o.mscip, varrefs, set.weights)
    ci = CI{VECTOR, SOS1}(cr.val)
    register!(o, ci)
    register!(o, cons(o, ci), cr)
    return ci
end

function MOI.delete(o::Optimizer, ci::CI{VECTOR, SOS1})
    allow_modification(o)
    delete!(o.constypes[VECTOR, SOS1], ConsRef(ci.value))
    delete!(o.reference, cons(o, ci))
    delete(o.mscip, ConsRef(ci.value))
    return nothing
end

function MOI.get(o::Optimizer, ::MOI.ConstraintFunction, ci::CI{VECTOR, SOS1})
    c = cons(o, ci)
    nvars::Int = SCIPgetNVarsSOS1(o, c)
    vars = unsafe_wrap(Array{Ptr{SCIP_VAR}}, SCIPgetVarsSOS1(o, c), nvars)
    return VECTOR([VI(ref(o, v).val) for v in vars])
end

function MOI.get(o::Optimizer, ::MOI.ConstraintSet, ci::CI{VECTOR, SOS1})
    c = cons(o, ci)
    nvars::Int = SCIPgetNVarsSOS1(o, c)
    weights = unsafe_wrap(Array{Float64}, SCIPgetWeightsSOS1(o, c), nvars)
    return SOS1(weights)
end

# Special-Ordered-Set type 2

MOI.supports_constraint(o::Optimizer, ::Type{VECTOR}, ::Type{SOS2}) = true

function MOI.add_constraint(o::Optimizer, func::VECTOR, set::SOS2)
    allow_modification(o)

    varrefs = [VarRef(vi.value) for vi in func.variables]
    cr = add_special_ordered_set_type2(o.mscip, varrefs, set.weights)
    ci = CI{VECTOR, SOS2}(cr.val)
    register!(o, ci)
    register!(o, cons(o, ci), cr)
    return ci
end

function MOI.delete(o::Optimizer, ci::CI{VECTOR, SOS2})
    allow_modification(o)
    delete!(o.constypes[VECTOR, SOS2], ConsRef(ci.value))
    delete!(o.reference, cons(o, ci))
    delete(o.mscip, ConsRef(ci.value))
    return nothing
end

function MOI.get(o::Optimizer, ::MOI.ConstraintFunction, ci::CI{VECTOR, SOS2})
    c = cons(o, ci)
    nvars::Int = SCIPgetNVarsSOS2(o, c)
    vars = unsafe_wrap(Array{Ptr{SCIP_VAR}}, SCIPgetVarsSOS2(o, c), nvars)
    return VECTOR([VI(ref(o, v).val) for v in vars])
end

function MOI.get(o::Optimizer, ::MOI.ConstraintSet, ci::CI{VECTOR, SOS2})
    c = cons(o, ci)
    nvars::Int = SCIPgetNVarsSOS2(o, c)
    weights = unsafe_wrap(Array{Float64}, SCIPgetWeightsSOS2(o, c), nvars)
    return SOS2(weights)
end