text:
# This file roughly corresponds to functions documented in the
# Assignments API: https://canvas.instructure.com/doc/api/assignments

"""
    Canvas.delete_assignment(c::Course, a::Assignment; kwargs...) -> Assignment

Delete the given assignment and return the former details. Return an [`Assignment`](@ref).

**Canvas API documentation**:
[Delete an assignment (`DELETE /api/v1/courses/:course_id/assignments/:id`)]
(https://canvas.instructure.com/doc/api/assignments#method.assignments.destroy)
"""
function delete_assignment(c::Course, a::Assignment; kwargs...)
    json = request("DELETE", "/api/v1$(Internals.pid(c))$(Internals.pid(a))"; kwargs...)
    return Assignment(json)
end

"""
    Canvas.assignments(c::Course; kwargs...) -> Vector{Assignment}, page_data

Return the paginated list of assignments for the course.
Return a vector of [`Assignment`](@ref)s and a dictionary with page data.

**Canvas API documentation**:
[List assignments (`GET /api/v1/courses/:course_id/assignments`)]
(https://canvas.instructure.com/doc/api/assignments#method.assignments_api.index)
"""
function assignments(c::Course; kwargs...)
    json, page_data = paged_request("GET", "/api/v1$(Internals.pid(c))/assignments"; kwargs...)
    return Assignment.(json), page_data
end

"""
    Canvas.assignments(c::Course, a::AssignmentGroup; kwargs...) -> Vector{Assignment}, page_data

Return the paginated list of assignments for the course and assignment group.
Return a vector of [`Assignment`](@ref)s and a dictionary with page data.

**Canvas API documentation**:
[List assignments (`GET /api/v1/courses/:course_id/assignment_groups/:assignment_group_id/assignments`)]
(https://canvas.instructure.com/doc/api/assignments#method.assignments_api.index)
"""
function assignments(c::Course, a::AssignmentGroup; kwargs...)
    json, page_data = paged_request("GET", "/api/v1$(Internals.pid(c))$(Internals.pid(a))/assignments"; kwargs...)
    return Assignment.(json), page_data
end

"""
    Canvas.assignments(u::User, c::Course; kwargs...) -> Vector{Assignment}, page_data

Return the paginated list of assignments for the specified user and course.
Return a vector of [`Assignment`](@ref)s and a dictionary with page data.

**Canvas API documentation**:
[List assignments for user(`GET /api/v1/users/:user_id/courses/:course_id/assignments`)]
(https://canvas.instructure.com/doc/api/assignments#method.assignments_api.user_index)
"""
function assignments(u::User, c::Course; kwargs...)
    json, page_data = paged_request("GET", "/api/v1$(Internals.pid(u))$(Internals.pid(c))/assignments"; kwargs...)
    return Assignment.(json), page_data
end

"""
    Canvas.assignment(c::Course, a::Assignment; kwargs...) -> Assignment

Return the assignment given by id. Return an [`Assignment`](@ref).

**Canvas API documentation**:
[Get a single assignment (`GET /api/v1/courses/:course_id/assignments/:id`)]
(https://canvas.instructure.com/doc/api/assignments#method.assignments_api.show)
"""
function assignment(c::Course, a::Assignment; kwargs...)
    json = request("GET", "/api/v1$(Internals.pid(c))$(Internals.pid(a))"; kwargs...)
    return Assignment(json)
end

"""
    Canvas.create_assignment(c::Course, a::Assignment; kwargs...) -> Assignment

Create a new assignment for the course and return the details. Return an [`Assignment`](@ref).

**Canvas API documentation**:
[Create an assignment (`POST /api/v1/courses/:course_id/assignments`)]
(https://canvas.instructure.com/doc/api/assignments#method.assignments_api.create)
"""
function create_assignment(c::Course; kwargs...)
    json = request("POST", "/api/v1$(Internals.pid(c))/assignments"; kwargs...)
    return Assignment(json)
end

"""
    Canvas.update_assignment(c::Course, a::Assignment; kwargs...) -> Assignment

Modify an existing assignment and return the new details. Return an [`Assignment`](@ref).

**Canvas API documentation**:
[Edit an assignment (`PUT /api/v1/courses/:course_id/assignments/:id`)]
(https://canvas.instructure.com/doc/api/assignments#method.assignments_api.update)
"""
function update_assignment(c::Course, a::Assignment; kwargs...)
    json = request("PUT", "/api/v1$(Internals.pid(c))$(Internals.pid(a))"; kwargs...)
    return Assignment(json)
end

"""
    Canvas.assignment_overrides(c::Course, a::Assignment; kwargs...) -> Vector{AssignmentOverride}, page_data

Return the paginated list of overrides for this course and assignment.
Return a vector of [`AssignmentOverride`](@ref)s and a dictionary with page data.

**Canvas API documentation**:
[List assignment overrides (`GET /api/v1/courses/:course_id/assignments/:assignment_id/overrides`)]
(https://canvas.instructure.com/doc/api/assignments#method.assignment_overrides.index)
"""
function assignment_overrides(c::Course, a::Assignment; kwargs...)
    json, page_data = paged_request("GET", "/api/v1$(Internals.pid(c))$(Internals.pid(a))/overrides"; kwargs...)
    return AssignmentOverride.(json), page_data
end

"""
    Canvas.assignment_override(c::Course, a::Assignment, o::AssignmentOverride; kwargs...) -> AssignmentOverride

Return the details of override specified by id. Return an [`AssignmentOverride`](@ref).

**Canvas API documentation**:
[Get a single assignment overrides (`GET /api/v1/courses/:course_id/assignments/:assignment_id/overrides/:id`)]
(https://canvas.instructure.com/doc/api/assignments#method.assignment_overrides.show)
"""
function assignment_override(c::Course, a::Assignment, o::AssignmentOverride; kwargs...)
    json = request("GET", "/api/v1$(Internals.pid(c))$(Internals.pid(a))$(Internals.pid(o))"; kwargs...)
    return AssignmentOverride(json)
end

"""
    Canvas.assignment_override(co::Union{Group,Section}, a::Assignment; kwargs...) -> AssignmentOverride

Return the override for the given group or section associated with the assignment (404 otherwise).
Return an [`AssignmentOverride`](@ref).

**Canvas API documentation**:
 - [Redirect to the assignment override for a group (`GET /api/v1/groups/:group_id/assignments/:assignment_id/override`)]
   (https://canvas.instructure.com/doc/api/assignments#method.assignment_overrides.group_alias)
 - [Redirect to the assignment override for a section (`GET /api/v1/sections/:course_section_id/assignments/:assignment_id/override`)]
   (https://canvas.instructure.com/doc/api/assignments#method.assignment_overrides.section_alias)
"""
function assignment_override(co::Union{Group,Section}, a::Assignment; kwargs...)
    json = request("GET", "/api/v1$(Internals.pid(co))$(Internals.pid(a))/override"; kwargs...)
    return AssignmentOverride(json)
end

"""
    Canvas.create_assignment_override(c::Course, a::Assignment; kwargs...) -> AssignmentOverride

Create an assignment override. One of `student_ids`, `group_id` or `course_section_id` must be present,
see API specification for details. Return the [`AssignmentOverride`](@ref).

**Canvas API documentation**:
[Create an assignment override (`POST /api/v1/courses/:course_id/assignments/:assignment_id/overrides`)]
(https://canvas.instructure.com/doc/api/assignments#method.assignment_overrides.create)
"""
function create_assignment_override(c::Course, a::Assignment; kwargs...)
    json = request("POST", "/api/v1$(Internals.pid(c))$(Internals.pid(a))/overrides"; kwargs...)
    return AssignmentOverride(json)
end

"""
    Canvas.update_assignment_override(c::Course, a::Assignment, o::AssignmentOverride; kwargs...) -> AssignmentOverride

Update the assignment override. All current overriden values must be supplied if the should
be retained, see API specification for details. Return the [`AssignmentOverride`](@ref).

**Canvas API documentation**:
[Update an assignment override (`PUT /api/v1/courses/:course_id/assignments/:assignment_id/overrides/:id`)]
(https://canvas.instructure.com/doc/api/assignments#method.assignment_overrides.update)
"""
function update_assignment_override(c::Course, a::Assignment, o::AssignmentOverride; kwargs...)
    json = request("PUT", "/api/v1$(Internals.pid(c))$(Internals.pid(a))$(Internals.pid(o))"; kwargs...)
    return AssignmentOverride(json)
end

"""
    Canvas.delete_assignment_override(c::Course, a::Assignment, o::AssignmentOverride; kwargs...) -> AssignmentOverride

Delete and override and return the former details. Return an [`AssignmentOverride`](@ref).

**Canvas API documentation**:
[Update an assignment override (`DELETE /api/v1/courses/:course_id/assignments/:assignment_id/overrides/:id`)]
(https://canvas.instructure.com/doc/api/assignments#method.assignment_overrides.update)
"""
function delete_assignment_override(c::Course, a::Assignment, o::AssignmentOverride; kwargs...)
    json = request("DELETE", "/api/v1$(Internals.pid(c))$(Internals.pid(a))$(Internals.pid(o))"; kwargs...)
    return AssignmentOverride(json)
end

# TODO:
# - Batch retrieve overrides in a course
# - Batch create overrides in a course
# - Batch update overrides in a course


text:
export GenomicPos,
    GenomicAxis

# ---- Genomic axis --------------------------------------------------------

"""
    GenomicAxis(;params...)

Create a genomic axis for the track plots.

# Parameters

- `:axis_line_width`
- `:axis_color`
- `:axis_font`
- `:axis_tick_length`

- `:width` width of track
- `:height` height of track

- `:margin_top`
- `:margin_bottom`
- `:margin_left`
- `:margin_right`

"""
struct GenomicAxis <: AbstractRefFrame
    # AbstractRefFrame
    transform::AbstractTransform
    children::Vector{AbstractRefFrame}
    params::Dict

    function GenomicAxis(;params...)
        new(Translation(0., 0.), AbstractRefFrame[], Dict(params))
    end
end

GenomicPos = GenomicAxis

function draw(cc::Cairo.CairoContext, tc::Vector{T}, gp::GenomicPos, xmin, xmax) where T <: AbstractRefFrame

    width = getparam(tc, :width)
    height = getparam(tc, :height, 1cm)

    units = Rect(xmin, 0, xmax - xmin, 1)
    bbox = Rect(getparam(tc, :margin_left),
                getparam(tc, :margin_top),
                width - getparam(tc, :margin_left) - getparam(tc, :margin_right),
                1)

    datactx = RefFrame(transform = Scaling(bbox, units))

    xticks = _ticks_default_linear([xmin, xmax])

    devx, devy = user2dev(vcat(tc, datactx), xticks, repmat([0], length(xticks)))

    axis_line_width = getparam(tc, :axis_line_width)
    Cairo.set_source_rgba(cc, rgba(getparam(tc, :axis_color))...)
    Cairo.set_font_face(cc, getparam(tc, :axis_font))
    Cairo.set_line_width(cc, axis_line_width)

    Cairo.move_to(cc, devx[1] - axis_line_width/2, devy[1])
    Cairo.line_to(cc, devx[end] + axis_line_width/2, devy[end])

    for i in 1:length(devx)
        Cairo.move_to(cc, devx[i], devy[i] - getparam(tc, :axis_tick_length))
        Cairo.line_to(cc, devx[i], devy[i])
        Cairo.text(cc, devx[i], devy[i] - 2*getparam(tc, :axis_tick_length), _format_ticklabel(xticks[i]),
                   halign="center", valign="bottom")
    end

    Cairo.stroke(cc)

    user2dev(gp, Rect(0.0, 0.0, width, height))
end

ext: jl


text:
# Use this simple script to record an interaction happening on `figure`

framerate = 10
total_time = 30 # in seconds
framen = framerate*total_time

record(figure, "name.mp4"; framerate = framerate) do io
    for i = 1:framen
        sleep(1/framerate)
        recordframe!(io)
    end
end


text:
using SpineOpt
using JuMP
using Clp

input_db_url = ARGS[1]
output_db_url = ARGS[2]
m = run_spineopt(input_db_url, output_db_url; with_optimizer=optimizer_with_attributes(Clp.Optimizer))

# Show active variables and constraints
println("*** Active constraints: ***")
for key in keys(m.ext[:constraints])
    !isempty(m.ext[:constraints][key]) && println(key)
end
println("*** Active variables: ***")
for key in keys(m.ext[:variables])
    !isempty(m.ext[:variables][key]) && println(key)
end

max_stars_repo_name: iasonask/nordic-hydro-spineopt

max_issues_repo_name: iasonask/nordic-hydro-spineopt

max_forks_repo_name: iasonask/nordic-hydro-spineopt


text:
export RecMLP, reset!

immutable RecMLP{C} <: MLP
    layers::Vector{Int}
    data::MLPData{C}
end

function RecMLP{C}(ctx::C, layers::Vector{Int})
    numLayers = length(layers)
    dims = Vector{WeightDim}()
    for i in 1:numLayers - 1
        rows = (layers[i] + 1)
        for next in i + 1:numLayers
            rows += layers[next]
        end
        cols = layers[i + 1]
        push!(dims, WeightDim(rows, cols))
    end
    RecMLP{C}(layers, MLPData(RecMLP{C}, ctx, dims))
end

function makeInputVector(mlp::RecMLP, index)
    input = addBias!(mlp.data, mlp.data.signal[index])
    batchSize = dims(input, 0)
    outputs = Nullable{arrayType(mlp.data)}()
    for next in index + 1:mlp.data.numLayers
        output = mlp.data.signal[next]
        outputSize = mlp.layers[next]
        if isEmpty(output)
            output = constant(mlp.data.ctx, 0f0, batchSize, outputSize)
        elseif dims(output, 0) != batchSize
            error("Size of batches must be the same. Last batch size was $(dims(output, 0)) and current is $batchSize.")
        end
        outputs = isnull(outputs) ? output : joinArrays(1, get(outputs), output)
    end
    joinArrays(1, input, get(outputs))
end

reset!(mlp::RecMLP) = reset!(mlp.data)

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
export resnet12
using Flux

# follow the implementation of https://github.com/kjunelee/MetaOptNet.git
function resnet12(use_avgpool = false)
    layer1 = BasicBlock(3 => 64, drop_prob=0.1, block_size=3)
    layer2 = BasicBlock(64 => 160, drop_prob=0.1, block_size=3)
    layer3 = BasicBlock(160 => 320, drop_prob=0.1, block_size=3, use_dropblock=true)
    layer4 = BasicBlock(320 => 640, drop_prob=0.1, block_size=3, use_dropblock=true)
    if use_avgpool
        pool = MeanPool((5, 5), pad=1, stride=2)
    else
        pool = identity
    end
    Chain(layer1, layer2, layer3, layer4, pool, (x)->reshape(x, :, size(x)[end]))
end


text:
using LowRankModels
using Base.Test

###############################################################
# verify all examples run
include("../examples/runexamples.jl")

###############################################################
# verify all tests of specific loss functions run
include("prob_tests/runtests.jl")

################################################################################
# ScikitLearnBase test
import ScikitLearnBase

# Check that KMeans can correctly separate two non-overlapping Gaussians
srand(21)
gaussian1 = randn(100, 2) + 5
gaussian2 = randn(50, 2) - 10
A = vcat(gaussian1, gaussian2)

model = ScikitLearnBase.fit!(LowRankModels.KMeans(), A)
@test Set(sum(ScikitLearnBase.transform(model, A), 1)) == Set([100, 50])


text:
module CalendarPuzzleTests

using CalendarPuzzle
using Dates
using Test

@testset "board" begin
    board_7x7 = CalendarPuzzle.makeboard((7, 7), 0, 0)
    @test board_7x7 == [
        0  0  0  0  0  0 -1
        0  0  0  0  0  0 -1
        0  0  0  0  0  0  0
        0  0  0  0  0  0  0
        0  0  0  0  0  0  0
        0  0  0  0  0  0  0
        0  0  0 -1 -1 -1 -1
    ]
    board = CalendarPuzzle.makeboard()
    @test board == [
         0  0  0  0  0  0 -1 -1
         0  0  0  0  0  0 -1 -1
         0  0  0  0  0  0  0 -1
         0  0  0  0  0  0  0 -1
         0  0  0  0  0  0  0 -1
         0  0  0  0  0  0  0 -1
         0  0  0 -1 -1 -1 -1 -1
        -1 -1 -1 -1 -1 -1 -1 -1
    ]

    board_May09 = CalendarPuzzle.makeboard(Date(2021, 5, 9))
    @test board_May09 == [
         0  0  0  0 -1  0 -1 -1
         0  0  0  0  0  0 -1 -1
         0  0  0  0  0  0  0 -1
         0 -1  0  0  0  0  0 -1
         0  0  0  0  0  0  0 -1
         0  0  0  0  0  0  0 -1
         0  0  0 -1 -1 -1 -1 -1
        -1 -1 -1 -1 -1 -1 -1 -1
    ]
end

@testset "solve" begin
    @test all(CalendarPuzzle.hasanswer(dt) for dt in Date(2020,1,1):Day(1):Date(2020,12,31))
end

end


text:
n = parse(Int, readline())
a = parse.(Int, split(readline()))
for i in n:-1:1
    if i == n
        print(a[i])
    else
        print(" ", a[i])
    end
end
println()


text:
@rule BIFM(:zprev, Marginalisation) (m_out::MultivariateNormalDistributionsFamily, m_in::MultivariateNormalDistributionsFamily, m_znext::MultivariateNormalDistributionsFamily, meta::BIFMMeta) = begin

    # fetch information from meta data
    A       = getA(meta)
    B       = getB(meta)
    C       = getC(meta)

    # fetch statistics from messages
    ξ_in, Λ_in          = weightedmean_precision(m_in)
    ξ_out, Λ_out        = weightedmean_precision(m_out)
    ξ_znext, Λ_znext    = weightedmean_precision(m_znext)

    # calculate intermediate quantities
    ξ_z = C' * ξ_out + ξ_znext
    Λ_z = C' * Λ_out * C + Λ_znext
    H = cholinv(Λ_in + B' * Λ_z * B)
    BHBt = B * H * B'
    ξ_ztilde = ξ_z + Λ_z * (B * (H * (-ξ_in - B' * ξ_z)))
    Λ_ztilde = Λ_z * (I - BHBt * Λ_z)

    # save required intermediate quantities
    setH!(meta, H)
    setξztilde!(meta, ξ_ztilde)
    setΛztilde!(meta, Λ_ztilde)
    setξz!(meta, ξ_z)
    setΛz!(meta, Λ_z)
    setBHBt!(meta, BHBt)

    # calculate outgoing message to zprev
    ξ_zprev = A' * ξ_ztilde
    Λ_zprev = A' * Λ_ztilde * A

    # Actual return type depends on meta object as well, so we explicitly cast the result here
    # Should be noop if type matches
    T = promote_type(eltype(m_out), eltype(m_in), eltype(m_znext))

    # return message
    return convert(MvNormalWeightedMeanPrecision{T}, MvNormalWeightedMeanPrecision(ξ_zprev, Λ_zprev))

end


text:
mutable struct OptimizationParameters
    d::Int # dimension of embedding
    # main optimization parameters
    max_iter::Int
    lr::Real
    beta_1::Real
    beta_2::Real
    weight_decay::Real
    eps::Real # adding to denominator of ADAM rescaling
    # early optimization parameters
    early_iter::Int
    early_lr::Real
    early_beta_1::Real
    early_beta_2::Real
    early_compression::Real
    early_exaggeration::Real
    # boolean flags
    dofast::Bool
    debias::Bool
    verbose::Bool
end

function OptimizationParameters(; d::Int = 2, max_iter::Int = 1024, lr::Real = 1,
            beta_1::Real = .9, beta_2::Real = .999, weight_decay::Real = 0,
            early_iter::Int = 64, early_lr::Real = lr, early_beta_1::Real = 0.5,
            early_beta_2::Real = 0.9, early_compression::Real = 1,
            early_exaggeration::Real = 12, eps::Real = 1e-8, verbose::Bool = false,
            debias::Bool = true, dofast::Bool = size(X, 2) ≥ fast_algorithm_min_size)
    early_exaggeration ≥ 1 || throw(DomainMismatch("early_exaggeration < 1"))
    OptimizationParameters(d, max_iter, lr, beta_1, beta_2, weight_decay, eps,
                early_iter, early_lr, early_beta_1, early_beta_2, early_compression,
                early_exaggeration, dofast, debias, verbose)
end


text:
struct BoxBound a::Float64; b::Float64; i::Int end # a * x[i] + b

resolve_lb(lb::Real) = iszero(lb) ? (-1.0, 0.0) : (-abs(inv(lb)),  1.0 * sign(lb)) # @code_warntype ✓
resolve_ub(ub::Real) = iszero(ub) ? ( 1.0, 0.0) : ( abs(inv(ub)), -1.0 * sign(ub)) # @code_warntype ✓

function boxbounds(lb::NTuple{ND}, ub::NTuple{ND}) where ND
    if @generated
        a = Vector{Expr}(undef, 2*ND)
        @inbounds for i in 1:ND
            a[i]    = :(BoxBound(resolve_lb(lb[$i])..., $i))
            a[i+ND] = :(BoxBound(resolve_ub(ub[$i])..., $i))
        end
        return quote
            $(Expr(:meta, :inline))
            @inbounds return $(Expr(:tuple, a...))
        end
    else
        return ntuple(i -> i > ND ? BoxBound(resolve_ub(ub[i - ND])..., i - ND) : BoxBound(resolve_lb(lb[i])..., i), 2*ND)
    end
end

eval_violation(x::VecI, bb::BoxBound) = max(0.0, bb.a * x[bb.i] + bb.b) # @code_warntype ✓
function eval_violation(x::VecI, cons::NTuple{NB,BoxBound}) where NB    # @code_warntype ✓
    if @generated
        a = Vector{Expr}(undef, NB)
        @inbounds for i in eachindex(a)
            a[i] = :(eval_violation(x, cons[$i]))
        end
        return quote
            $(Expr(:meta, :inline))
            @inbounds return $(Expr(:call, :+, a...))
        end
    else
        ret = 0.0
        @inbounds for i in eachindex(cons)
            ret += eval_violation(x, cons[i])
        end
        return ret
    end
end

# check feasibility of agent in throng, @code_warntype ✓
function check!(xnew::VecI, agents::VecI{Agent}, elites::VecI{Agent}, throng::VecI{Agent}, edx::Int, tdx::Int, fn::Function, cons::NTuple)
    violation = eval_violation(xnew, cons)
    violation > 0.0 && return @inbounds check!(xnew, violation, throng[tdx]) # x[new] is infeasible
    return check!(xnew, fcall(fn, xnew), agents, elites, throng, edx, tdx)   # x[new] is feasible
end

# Matchup for a feasible x[new] agent in throng, @code_warntype ✓
function check!(xnew::VecI, fnew::Real, agents::VecI{Agent}, elites::VecI{Agent}, throng::VecI{Agent}, edx::Int, tdx::Int)
    @inbounds xold = throng[tdx]
    # x[old] is infeasible
    if !xold.v
        xold.f = fnew
        xold.v = true
        xold.c = 0.0
        copy!(xold.x, xnew)
        return nothing
    end
    # x[old], x[new] are feasible, x[new] is better than/equals to x[best] (greedy strategy)
    if @inbounds !(elites[edx].f < fnew)
        xold.f = fnew
        copy!(xold.x, xnew)
        swap!(agents, edx, length(elites) + tdx)
        return nothing
    end
    # x[old], x[new] are feasible
    if !(xold.f < fnew)
        xold.f = fnew
        copy!(xold.x, xnew)
        return nothing
    end
end

# check feasibility of agent in elites, @code_warntype ✓
function check!(xnew::VecI, elites::VecI{Agent}, edx::Int, fn::Function, cons::NTuple)
    violation = eval_violation(xnew, cons)
    violation > 0.0 && return @inbounds check!(xnew, violation, elites[edx]) # x[new] is infeasible
    return check!(xnew, fcall(f, xnew), elites, edx)                         # x[new] is feasible
end

# Matchup for a feasible x[new] trial in elites, @code_warntype ✓
function check!(xnew::VecI, fnew::Real, elites::VecI{Agent}, edx::Int)
    @inbounds elite = elites[edx]
    # x[old] is infeasible
    if !elite.v
        elite.f = fnew
        elite.v = true
        elite.c = 0.0
        copy!(elite.x, xnew)
        return nothing
    end
    # x[old], x[new] are feasible, x[new] is better than/equals to x[best] (greedy strategy)
    if @inbounds !(elites[1].f < fnew)
        elite.f = fnew
        copy!(elite.x, xnew)
        swap!(elites, 1, edx)
        return nothing
    end
    # x[old], x[new] are feasible
    if !(elite.f < fnew)
        elite.f = fnew
        copy!(elite.x, xnew)
        return nothing
    end
end

# Matchup for an infeasible x[new] trial, here "fnew = violation", @code_warntype ✓
function check!(xnew::VecI, violation::Real, agent::Agent)
    # x[old], x[new] are infeasible, compare violation
    # There is no `else` condition, if x[old] is feasible, then a matchup is unnecessary.
    if !agent.v && !(agent.c < violation)
        agent.c = violation
        copy!(agent.x, xnew)
    end
    return nothing
end


text:
using MultiModalMuSig
using Test

K = 2
α = 0.1
η = 0.1
X = Matrix{Int}[
    [
        1 5;
        2 8
    ],
    [
        1 2;
        2 5
    ]
]

@testset "constructor" begin
    model = MultiModalMuSig.LDA(K, α, η, X)
    @test model.K == 2
    @test model.D == 2
    @test model.N == [13, 7]
    @test model.V == 2

    @test size(model.λ) == (2, 2)
    @test all(model.λ .> 0) == true

    @test size(model.γ) == (2, 2)
    @test all(model.γ .> 0) == true

    @test sum(model.ϕ[1], dims=1) ≈ ones(2)'

    model = MultiModalMuSig.LDA(K, α, η, 3, X)
    @test model.V == 3
    @test size(model.λ) == (3, 2)
end

@testset "update_ϕ!" begin
    Elnθ = [0.5  -1.1; 2.3 -0.7]
    Elnβ = [-0.2 -0.9; -1.1 0.3]

    ϕ = Array{Float64}(undef, 2, 2)
    # K1 W1
    ϕ[1, 1] = exp(Elnθ[1, 1] + Elnβ[1, 1])
    # K1 W2
    ϕ[1, 2] = exp(Elnθ[1, 1] + Elnβ[2, 1])
    # K2 W1
    ϕ[2, 1] = exp(Elnθ[2, 1] + Elnβ[1, 2])
    # K2 W2
    ϕ[2, 2] = exp(Elnθ[2, 1] + Elnβ[2, 2])

    ϕ[:, 1] ./= sum(ϕ[:, 1])
    ϕ[:, 2] ./= sum(ϕ[:, 2])

    model = MultiModalMuSig.LDA(K, α, η, X)
    model.Elnθ .= Elnθ
    model.Elnβ .= Elnβ
    MultiModalMuSig.update_ϕ!(model)

    @test model.ϕ[1] ≈ ϕ
end

@testset "update_γ!" begin
    ϕ = [0.4 0.2; 0.6 0.8]

    γ = Array{Float64}(undef, 2)
    γ[1] = α + (ϕ[1, 1] * X[1][1, 2]) + (ϕ[1, 2] * X[1][2, 2])
    γ[2] = α + (ϕ[2, 1] * X[1][1, 2]) + (ϕ[2, 2] * X[1][2, 2])

    Elnθ = Array{Float64}(undef, 2)
    Elnθ[1] = digamma(γ[1]) - digamma(γ[1] + γ[2])
    Elnθ[2] = digamma(γ[2]) - digamma(γ[1] + γ[2])

    model = MultiModalMuSig.LDA(K, α, η, X)
    model.ϕ[1] = ϕ
    MultiModalMuSig.update_γ!(model)

    @test model.γ[:, 1] ≈ γ
    @test model.Elnθ[:, 1] ≈ Elnθ
end

@testset "update_λ!" begin
    ϕ = [[0.4 0.2; 0.6 0.8], [0.1 0.6; 0.9 0.4]]

    λ = Array{Float64}(undef, 2, 2)
    λ[1, 1] = η + (ϕ[1][1, 1] * X[1][1, 2]) + (ϕ[2][1, 1] * X[2][1, 2])
    λ[2, 1] = η + (ϕ[1][1, 2] * X[1][2, 2]) + (ϕ[2][1, 2] * X[2][2, 2])
    λ[1, 2] = η + (ϕ[1][2, 1] * X[1][1, 2]) + (ϕ[2][2, 1] * X[2][1, 2])
    λ[2, 2] = η + (ϕ[1][2, 2] * X[1][2, 2]) + (ϕ[2][2, 2] * X[2][2, 2])

    Elnβ = Array{Float64}(undef, 2, 2)
    Elnβ[1, 1] = digamma(λ[1, 1]) - digamma(λ[1, 1] + λ[2, 1])
    Elnβ[2, 1] = digamma(λ[2, 1]) - digamma(λ[1, 1] + λ[2, 1])
    Elnβ[1, 2] = digamma(λ[1, 2]) - digamma(λ[1, 2] + λ[2, 2])
    Elnβ[2, 2] = digamma(λ[2, 2]) - digamma(λ[1, 2] + λ[2, 2])

    model = MultiModalMuSig.LDA(K, α, η, X)
    model.ϕ = ϕ
    MultiModalMuSig.update_λ!(model)

    @test model.λ ≈ λ
    @test model.Elnβ ≈ Elnβ
end

@testset "calculate_elbo" begin
    model = MultiModalMuSig.LDA(K, α, η, X)

    @test typeof(MultiModalMuSig.calculate_ElnPβ(model)) == Float64
    @test typeof(MultiModalMuSig.calculate_ElnPθ(model)) == Float64
    @test typeof(MultiModalMuSig.calculate_ElnPZ(model)) == Float64
    @test typeof(MultiModalMuSig.calculate_ElnPX(model)) == Float64

    @test typeof(MultiModalMuSig.calculate_ElnQβ(model)) == Float64
    @test typeof(MultiModalMuSig.calculate_ElnQθ(model)) == Float64
    @test typeof(MultiModalMuSig.calculate_ElnQZ(model)) == Float64

    @test MultiModalMuSig.calculate_elbo(model) < 0.0
end


text:
"""
    TelescopeIterator(v, n[, memory=typemax])

Iterator that is formally equivalent to
```julia
for i = 1:length(v), j = max(1,i-memory):i-1, k = max(1,j-memory):j-1, ...
    e = v[i,j,k,...]
end
```
for `n` loop variables.

# Example

```jldoctest
julia> ti = StrongFieldApproximation.TelescopeIterator(1:6, 2, 3)
StrongFieldApproximation.TelescopeIterator{UnitRange{Int64}, Int64}(1:6, 2, 3)

julia> for e in ti
       println(e)
       end
[2, 1]
[3, 1]
[4, 1]
[3, 2]
[4, 2]
[5, 2]
[4, 3]
[5, 3]
[6, 3]
[5, 4]
[6, 4]
[6, 5]
```
"""
struct TelescopeIterator{V<:AbstractVector, I<:Integer}
    v::V
    n::I
    memory::I
    function TelescopeIterator(v::V, n::I, memory::I=typemax(I)) where {V<:AbstractVector, I<:Integer}
        (n < 0 || memory < 0) && throw(ArgumentError("Invalid values for n = $n or memory = $memory"))
        new{V,I}(v, n, memory)
    end
end

function Base.iterate(ti::TelescopeIterator, state=[ti.n+1-i for i=1:ti.n])
    (isempty(ti.v) || iszero(ti.n) || ti.n > length(ti.v) || ti.n > 1 && iszero(ti.memory)) && return nothing

    nv = length(ti.v)
    any(>(nv), state) && return nothing

    val = ti.v[state]

    if ti.n == 1 && state[1] < nv
        state[1] += 1
        return val, state
    end

    for j = 1:ti.n
        state[j] > nv || j == ti.n && state[j] == nv-(ti.n-2) && return nothing
        nsj = state[j] + 1
        j < ti.n && nsj - state[j+1] ≤ ti.memory && nsj ≤ nv-(j-1) || j == ti.n || continue
        state[j] = nsj
        for i = j-1:-1:1
            state[i] = state[i+1]+1
        end
        break
    end

    val, state
end


text:
#!/usr/bin/env julia

using MIToS.Utils.Scripts

Args = parse_commandline(
# TO DO -----------------------------------------------------------------------
    description="""
    Creates a file in Stockholm format with the aligned columns from a Pfam Stockholm file.
    Insertions are deleted, as they are unaligned in a proﬁle HMM.
    The output file *.aligned.* contains UniProt residue numbers and original column numbers in its annotations.
    """,
    output=".aligned."
# -----------------------------------------------------------------------------
    )

set_parallel(Args["parallel"])

@everywhere begin

    const args = remotecall_fetch(1,()->Args)

    import MIToS.Utils.Scripts: script

    # TO DO -----------------------------------------------------------------------
    using MIToS.MSA

    function script(input::Union{Base.LibuvStream,  AbstractString},
                    args,
                    fh_out::Union{Base.LibuvStream, IO})
        try
            aln = readorparse(input, Stockholm, generatemapping=true, useidcoordinates=true, deletefullgaps=true)
            print(fh_out, aln, Stockholm)
        catch err
            warn(string("ERROR for ", input, ": ", err))
        end

    end
    # -----------------------------------------------------------------------------

end

runscript(args)


text:
module TestTransducers

import PerformanceTestTools
using Distributed: addprocs, nworkers
using Test

if get(ENV, "CI", "false") == "true"
    addprocs(3)
end
@info "Testing with:" nworkers()

@testset "$file" for file in sort([file for file in readdir(@__DIR__) if
                                   match(r"^test_.*\.jl$", file) !== nothing])

    # BlockArrays.jl 0.11 only support Julia 1.1 and above.
    VERSION < v"1.1" && file == "test_blockarrays.jl" && continue

    include(file)
end

PerformanceTestTools.@include_foreach(
    "threads/runtests.jl",
    [nothing, ["JULIA_NUM_THREADS" => Threads.nthreads() > 1 ? "1" : "2"]],
    parallel = true,
)

end  # module


text:
# --------------------------------------------------------------------------- #
# collect the un-hidden filenames available in a given folder

function getfname(pathtofolder)

    # recupere tous les fichiers se trouvant dans le repertoire cible
    allfiles = readdir(pathtofolder)

    # vecteur booleen qui marque les noms de fichiers valides
    flag = trues(size(allfiles))

    k=1
    for f in allfiles
        # traite chaque fichier du repertoire
        if f[1] != '.'
            # pas un fichier cache => conserver
            println("fname = ", f)
        else
            # fichier cache => supprimer
            flag[k] = false
        end
        k = k+1
    end

    # extrait les noms valides et retourne le vecteur correspondant
    finstances = allfiles[flag]
    return finstances
end

# --------------------------------------------------------------------------- #
# Loading an instance of SPP (format: OR-library)

function loadSPP(fname)

    f=open(fname)
    # lecture du nbre de contraintes (m) et de variables (n)
    m, n = parse.(Int, split(readline(f)) )
    # lecture des n coefficients de la fonction economique et cree le vecteur d'entiers c
    C = parse.(Int, split(readline(f)) )
    # lecture des m contraintes et reconstruction de la matrice binaire A
    A=zeros(Int, m, n)
    for i=1:m
        # lecture du nombre d'elements non nuls sur la contrainte i (non utilise)
        readline(f)
        # lecture des indices des elements non nuls sur la contrainte i
        for valeur in split(readline(f))
          j = parse(Int, valeur)
          A[i,j]=1
        end
    end
    close(f)
    return C, A
end

# --------------------------------------------------------------------------- #
# Deamon checking the feasibility of a solution given its vector x

function isAdmissible(C, A, x)

    vecSat = zeros(Int, size(A,1))
    vecUnit = ones(Int,size(A,1))
    z::Int64 = 0
    verbose = true
    var1 = findall(isequal(1), x[:])

    for j in var1
        vecSat = vecSat .+ A[:,j]
        z = z + C[j]
    end

    if findfirst(isequal(false), (vecSat .<= vecUnit)) != nothing
        println( "Feasible : no")
        @assert false "no-feasible solution detected"
    end
    println( "Feasible : yes | Σ(x_i) = ", length(var1), " ; z(x) = ", z)
    return true
end

max_stars_repo_name: xgandibleux/solveSPP

max_issues_repo_name: xgandibleux/solveSPP

max_forks_repo_name: xgandibleux/solveSPP


text:
function filter_cutest(problem)
    correct_size = 50 <= problem["variables"]["number"] + problem["constraints"]["number"] && problem["constraints"]["number"] >= 10 && problem["variables"]["number"] <= 600 && problem["constraints"]["number"] <= 1000
    regular = problem["derivative_order"] >= 2 && problem["regular"] == true
    if correct_size && regular
        return true
    else
      return false
    end
end

#problem_list = CUTEst.select(max_var=100, max_var=1000, min_con=100, max_con=3000)
problem_list = CUTEst.select(custom_filter=filter_cutest)
problem_list = convert(Array{String,1},problem_list)

for problem_name in problem_list
  cute = CUTEstModel(problem_name)

  J = jac(cute, cute.meta.x0)
  H = hess(cute, cute.meta.x0, y=ones(size(J,2)))
  @show sum(H .> 0.0) / sum(J .> 0.0)
  finalize(cute)
end


text:
function build_channel(ch::IonChannel; name = get_name(ch))
    @variables V(t) Ca(t)
    eqs, states, parameters, current, defaultmap = channel_dynamics(ch, V, Ca)
    return ODESystem(eqs, t, [V, Ca, states...], [parameters...];
                                            observed=current,
                                            name = name,
                                            defaults = defaultmap)
end

function build_channel(syn::Synapse; name = get_name(syn))
    @variables Vpre(t) Vpost(t)
    eqs, states, parameters, current, defaultmap = channel_dynamics(syn, Vpre, Vpost)
    return ODESystem(eqs, t, [Vpre, Vpost, states...], [parameters...];    observed = current,
                        name = name,
                        defaults = defaultmap)
end

# 10 nF/mm² for default specific membrane capacitance
function build_neuron(comp, channels; name = :unidentified_neuron)

    neuron_parameters = @parameters Cₘ τCa Ca∞
    neuron_states = @variables V(t) Ca(t)
    syns = [Num(Variable{Symbolics.FnType{Tuple{Any},Real}}(Symbol(:Isyn, i)))(t) for i in 1:comp.hooks]

    channel_systems = build_channel.(channels)

    summed_membrane_currents = sum(ionic_current(ch, ch_sys) for (ch, ch_sys) in zip(channels, channel_systems))

    summed_calcium_flux = sum(calcium_current(ch, ch_sys) for (ch, ch_sys) in zip(channels, channel_systems))

    connections = cat(
                    [voltage_hook(V, ch, chs) for (ch, chs) in zip(channels, channel_systems)],
                    [calcium_hook(Ca, ch, chs) for (ch, chs) in zip(channels, channel_systems)],
                    dims=1
                    )
    f = 0.094
    diffeqs =   cat(
                    [D(V) ~ (summed_membrane_currents + sum(syns))/Cₘ],
                    [D(Ca) ~ (1/τCa)*(-Ca + Ca∞ + f*summed_calcium_flux)],
                    dims=1
                    )

    eqs = cat(connections, diffeqs; dims=1)
    eqs = filter(x -> !(x == (Num(0) ~ Num(0))), eqs)

    channel_defs = Dict{SymbolicUtils.Symbolic{Real}, Float64}()
    for (ch, chs) in zip(channels, channel_systems)
        param_names = external_params(ch)
        if !isnothing(param_names)
            for el in param_names
                channel_defs[getproperty(chs, el)] = comp.parameters[el]
            end
        end
    end

    soma_defs = Dict{Num, Float64}()
    for el in neuron_parameters
        soma_defs[el] = comp.parameters[ModelingToolkit.tosymbol(el)]
    end
    state_defs = Dict(V => comp.initial_states[:V], Ca => comp.initial_states[:Ca])
    defaults = merge(soma_defs, channel_defs, state_defs)
    neur = ODESystem(eqs, t;
                    systems = channel_systems,
                    defaults = defaults,
                    name = name
                    )

    return neur
end

function build_group(arr_of_neurons; name = :group)
    eqs = Array{Equation,1}(undef,0)
    ODESystem(eqs, t, [], []; name=name, systems = arr_of_neurons)
end

"""
Issue: need to do all connections in one ḡChol

"""
function add_connection(group, pre, post, syn::Synapse; name = ModelingToolkit.getname(group), i=1)
    prename, postname = ModelingToolkit.getname.([pre, post])
    synapse = build_channel(syn; name = Symbol(prename, :to, postname, get_name(syn)))

    oldeqs = ModelingToolkit.get_eqs(group)
    neweqs = [
            get_pre(syn, pre) ~ my_pre(syn, synapse),
            get_post(syn, post) ~ my_post(syn, synapse),
            ionic_current(syn, synapse) ~ getproperty(post, Symbol(post_connector(syn), i))
    ]

    eqs = cat(oldeqs, neweqs, dims=1)
    systems = cat(ModelingToolkit.get_systems(group), synapse, dims=1)

    return connected = ODESystem(eqs, t, [], [];
                                name=name,
                                systems = systems)
end

max_stars_repo_name: wsphillips/NeuronBuilder

max_issues_repo_name: wsphillips/NeuronBuilder

max_forks_repo_name: wsphillips/NeuronBuilder


text:
# Policy for playing back a sequence of actions
mutable struct PlaybackPolicy <: Policy
    actions
    backup_policy::Policy
    logpdfs::Vector{Float64} # The log probabilities of each action
    i::Int64
end

# Constructor for the PlaybackPolicy
PlaybackPolicy(actions, backup_policy::Policy) = PlaybackPolicy(actions, backup_policy, [], 1)
PlaybackPolicy(actions, logpdfs::Vector{Float64}, backup_policy::Policy) = PlaybackPolicy(actions, backup_policy, logpdfs, 1)

# Action selection for the PlaybackPolicy
function POMDPs.action(p::PlaybackPolicy, s)
    a = p.i <= length(p.actions) ? p.actions[p.i] : action(p.backup_policy, s)
    p.i += 1
    a
end

# Get the logpdf of the history from the playback policy and the backup policy
function Distributions.logpdf(p::PlaybackPolicy, h)
    N = min(length(p.actions), length(h))
    @assert all(collect(action_hist(h))[1:N] .== p.actions[1:N])
    @assert length(p.actions) == length(p.logpdfs)
    if length(h) > N
        return sum(p.logpdfs) + sum(logpdf(p.backup_policy, view(h, N+1:length(h))))
    else
        return sum(p.logpdfs[1:N])
    end
end


text:
const BlockSparseMatrix{ElT,StoreT,IndsT} = BlockSparseTensor{ElT,2,StoreT,IndsT}
const DiagBlockSparseMatrix{ElT,StoreT,IndsT} = DiagBlockSparseTensor{ElT,2,StoreT,IndsT}
const DiagMatrix{ElT,StoreT,IndsT} = DiagTensor{ElT,2,StoreT,IndsT}

function _truncated_blockdim(S::DiagMatrix,
                             docut::Float64;
                             singular_values=false,
                             truncate=true)
  !truncate && return diaglength(S)
  newdim = 0
	val = singular_values ? getdiagindex(S,newdim+1)^2 : getdiagindex(S,newdim+1)
  while newdim+1 ≤ diaglength(S) && val > docut
    newdim += 1
    if newdim+1 ≤ diaglength(S)
      val = singular_values ? getdiagindex(S,newdim+1)^2 :
             getdiagindex(S,newdim+1)
    end
  end
  return newdim
end

"""
    svd(T::BlockSparseTensor{<:Number,2}; kwargs...)

svd of an order-2 BlockSparseTensor.

This function assumes that there is one block
per row/column, otherwise it fails.
This assumption makes it so the result can be
computed from the dense svds of seperate blocks.
"""
function LinearAlgebra.svd(T::BlockSparseMatrix{ElT};
                           kwargs...) where {ElT}
  alg::String = get(kwargs, :alg, "recursive")

  truncate = haskey(kwargs, :maxdim) || haskey(kwargs, :cutoff)

  Us = Vector{BlockSparseMatrix{ElT}}(undef, nnzblocks(T))
  Ss = Vector{DiagBlockSparseMatrix{real(ElT)}}(undef, nnzblocks(T))
  Vs = Vector{BlockSparseMatrix{ElT}}(undef, nnzblocks(T))

  # Sorted eigenvalues
  d = Vector{real(ElT)}()

  for n in 1:nnzblocks(T)
    b = nzblock(T, n)
    blockT = blockview(T, n)
    Ub, Sb, Vb = svd(blockT; alg = alg)
    Us[n] = Ub
    Ss[n] = Sb
    Vs[n] = Vb
    append!(d, vector(diag(Sb)))
  end

  # Square the singular values to get
  # the eigenvalues
  d .= d .^ 2
  sort!(d; rev=true)

  # Get the list of blocks of T
  # that are not dropped
  nzblocksT = nzblocks(T)

  dropblocks = Int[]
  if truncate
    truncerr,docut = truncate!(d; kwargs...)
    for n in 1:nnzblocks(T)
      blockdim = _truncated_blockdim(Ss[n],
                                     docut;
                                     singular_values=true,
                                     truncate=truncate)
      if blockdim == 0
        push!(dropblocks,n)
      else
        Strunc = tensor(Diag(store(Ss[n])[1:blockdim]),
                        (blockdim,blockdim))
        Us[n] = copy(Us[n][1:dim(Us[n],1),1:blockdim])
        Ss[n] = Strunc
        Vs[n] = copy(Vs[n][1:dim(Vs[n],1),1:blockdim])
      end
    end
    deleteat!(Us,dropblocks)
    deleteat!(Ss,dropblocks)
    deleteat!(Vs,dropblocks)
    deleteat!(nzblocksT,dropblocks)
    else
    truncerr,docut = 0.0,0.0
  end

  # The number of blocks of T remaining
  nnzblocksT = nnzblocks(T)-length(dropblocks)

  #
  # Put the blocks into U,S,V
  #

  nb1_lt_nb2 = (nblocks(T)[1] < nblocks(T)[2] ||
                (nblocks(T)[1] == nblocks(T)[2] &&
                 dim(T,1) < dim(T,2)))

  if nb1_lt_nb2
    uind = sim(ind(T,1))
  else
    uind = sim(ind(T,2))
  end

  deleteat!(uind,dropblocks)

  # uind may have too many blocks
  if nblocks(uind) > nnzblocksT
    resize!(uind,nnzblocksT)
  end

  for n in 1:nnzblocksT
    setblockdim!(uind, minimum(dims(Ss[n])), n)
  end

  if dir(uind) != dir(inds(T)[1])
    uind = dag(uind)
  end
  indsU = setindex(inds(T), dag(uind), 2)

  vind = sim(uind)
  if dir(vind) != dir(inds(T)[2])
    vind = dag(vind)
  end
  indsV = setindex(inds(T),dag(vind),1)
  indsV = permute(indsV,(2,1))

  indsS = setindex(inds(T),uind,1)
  indsS = setindex(indsS,vind,2)

  nzblocksU = Vector{Block{2}}(undef,nnzblocksT)
  nzblocksS = Vector{Block{2}}(undef,nnzblocksT)
  nzblocksV = Vector{Block{2}}(undef,nnzblocksT)

  for n in 1:nnzblocksT
    blockT = nzblocksT[n]

    blockU = (blockT[1],n)
    nzblocksU[n] = blockU

    blockS = (n,n)
    nzblocksS[n] = blockS

    blockV = (blockT[2],n)
    nzblocksV[n] = blockV
  end

  U = BlockSparseTensor(ElT, undef, nzblocksU, indsU)
  S = DiagBlockSparseTensor(real(ElT), undef, nzblocksS, indsS)
  V = BlockSparseTensor(ElT, undef, nzblocksV, indsV)

  for n in 1:nnzblocksT
    Ub, Sb, Vb = Us[n], Ss[n], Vs[n]

    blockU = nzblocksU[n]
    blockS = nzblocksS[n]
    blockV = nzblocksV[n]

    blockview(U, blockU) .= Ub
    blockviewS = blockview(S, blockS)
    for i in 1:diaglength(Sb)
      setdiagindex!(blockviewS, getdiagindex(Sb, i), i)
    end
    blockview(V, blockV) .= Vb
  end

  return U,S,V,Spectrum(d,truncerr)
end

_eigen_eltypes(T::Hermitian{ElT,<:BlockSparseMatrix{ElT}}) where {ElT} = real(ElT), ElT

_eigen_eltypes(T::BlockSparseMatrix{ElT}) where {ElT} = complex(ElT), complex(ElT)

function LinearAlgebra.eigen(T::Union{Hermitian{ElT,<:BlockSparseMatrix{ElT}},
                                      BlockSparseMatrix{ElT}};
                             kwargs...) where {ElT<:Union{Real,Complex}}
  truncate = haskey(kwargs,:maxdim) || haskey(kwargs,:cutoff)

  ElD, ElV = _eigen_eltypes(T)

  # Sorted eigenvalues
  d = Vector{real(ElT)}()

  b = nzblock(T, 1)
  all(==(b[1]), b) || error("Eigen currently only supports block diagonal matrices.")
  blockT = blockview(T, 1)
  Db, Vb = eigen(blockT)
  Ds = [Db]
  Vs = [Vb]
  append!(d, abs.(vector(diag(Db))))
  for n in 2:nnzblocks(T)
    b = nzblock(T, n)
    all(==(b[1]), b) || error("Eigen currently only supports block diagonal matrices.")
    blockT = blockview(T, n)
    Db, Vb = eigen(blockT)
    push!(Ds, Db)
    push!(Vs, Vb)
    append!(d, abs.(vector(diag(Db))))
  end

  dropblocks = Int[]
  sort!(d; rev=true, by=abs)

  if truncate
    truncerr,docut = truncate!(d; kwargs...)
    for n in 1:nnzblocks(T)
      blockdim = _truncated_blockdim(Ds[n],docut)
      if blockdim == 0
        push!(dropblocks,n)
      else
        Dtrunc = tensor(Diag(store(Ds[n])[1:blockdim]),
                        (blockdim,blockdim))
        Ds[n] = Dtrunc
        Vs[n] = copy(Vs[n][1:dim(Vs[n],1),1:blockdim])
      end
    end
    deleteat!(Ds,dropblocks)
    deleteat!(Vs,dropblocks)
  else
    truncerr = 0.0
  end

  # Get the list of blocks of T
  # that are not dropped
  nzblocksT = nzblocks(T)
  deleteat!(nzblocksT, dropblocks)

  # The number of blocks of T remaining
  nnzblocksT = nnzblocks(T) - length(dropblocks)

  #
  # Put the blocks into D, V
  #

  i1, i2 = inds(T)
  l = sim(i1)

  deleteat!(l, dropblocks)

  # l may have too many blocks
  if nblocks(l) > nnzblocksT
    resize!(l, nnzblocksT)
  end

  # Truncation may have changed
  # some block sizes
  for n in 1:nnzblocksT
    setblockdim!(l, minimum(dims(Ds[n])), n)
  end

  r = dag(sim(l))

  indsD = (l, r)
  indsV = (dag(i2), r)

  nzblocksD = Vector{Block{2}}(undef, nnzblocksT)
  nzblocksV = Vector{Block{2}}(undef, nnzblocksT)
  for n in 1:nnzblocksT
    blockT = nzblocksT[n]

    blockD = (n, n)
    nzblocksD[n] = blockD

    blockV = (blockT[1], n)
    nzblocksV[n] = blockV
  end

  D = DiagBlockSparseTensor(ElD, undef, nzblocksD, indsD)
  V = BlockSparseTensor(ElV, undef, nzblocksV, indsV)

  for n in 1:nnzblocksT
    Db, Vb = Ds[n], Vs[n]

    blockD = nzblocksD[n]
    blockviewD = blockview(D,blockD)
    for i in 1:diaglength(Db)
      setdiagindex!(blockviewD, getdiagindex(Db, i), i)
    end

    blockV = nzblocksV[n]
    blockview(V, blockV) .= Vb
  end

  return D, V, Spectrum(d, truncerr)
end

function LinearAlgebra.exp(T::Union{BlockSparseMatrix{ElT},
                                    Hermitian{ElT, <:BlockSparseMatrix{ElT}}}) where {ElT<: Union{Real,
                                                                                Complex}}
  expT = BlockSparseTensor(ElT, undef, nzblocks(T), inds(T))
  for n in 1:nnzblocks(T)
    b = nzblock(T,n)
    all(==(b[1]),b) || error("exp currently supports only block-diagonal matrices")
    blockT = blockview(T,n)
    blockview(expT,n) .= exp(blockT)
  end
  return expT
end


text:
struct AxisPlot
    axis
    plot::AbstractPlot
end

Base.show(io::IO, fap::FigureAxisPlot) = show(io, fap.figure)
Base.show(io::IO, ::MIME"text/plain", fap::FigureAxisPlot) = print(io, "FigureAxisPlot()")

Base.iterate(fap::FigureAxisPlot, args...) = iterate((fap.figure, fap.axis, fap.plot), args...)
Base.iterate(ap::AxisPlot, args...) = iterate((ap.axis, ap.plot), args...)

get_scene(ap::AxisPlot) = get_scene(ap.axis.scene)

function plot(P::PlotFunc, args...; axis = NamedTuple(), figure = NamedTuple(), kw_attributes...)
    # scene_attributes = extract_scene_attributes!(attributes)
    fig = Figure(; figure...)

    axis = Dict(pairs(axis))
    if haskey(axis, :type)
        axtype = axis[:type]
        pop!(axis, :type)
        ax = axtype(fig; axis...)
    else
        proxyscene = Scene()
        attrs = Attributes(kw_attributes)
        delete!(attrs, :show_axis)
        delete!(attrs, :limits)
        plot!(proxyscene, P, attrs, args...)
        if is2d(proxyscene)
            ax = Axis(fig; axis...)
        else
            ax = LScene(fig; scenekw = (camera=cam3d!, axis...))
        end
    end

    fig[1, 1] = ax
    p = plot!(ax, P, Attributes(kw_attributes), args...)
    FigureAxisPlot(fig, ax, p)
end

# without scenelike, use current axis of current figure

function plot!(P::PlotFunc, args...; kw_attributes...)
    ax = current_axis(current_figure())
    isnothing(ax) && error("There is no current axis to plot into.")
    plot!(P, ax, args...; kw_attributes...)
end

function plot(P::PlotFunc, gp::GridPosition, args...; axis = NamedTuple(), kwargs...)

    f = MakieLayout.get_top_parent(gp)

    c = contents(gp, exact = true)
    if !isempty(c)
        error("""
        You have used the non-mutating plotting syntax with a GridPosition, which requires an empty GridLayout slot to create an axis in, but there are already the following objects at this layout position:

        $(c)

        If you meant to plot into an axis at this position, use the plotting function with `!` (e.g. `func!` instead of `func`).
        If you really want to place an axis on top of other layoutables, make your intention clear and create it manually.
        """)
    end

    axis = Dict(pairs(axis))

    if haskey(axis, :type)
        axtype = axis[:type]
        pop!(axis, :type)
        ax = axtype(f; axis...)
    else
        proxyscene = Scene()
        plot!(proxyscene, P, Attributes(kwargs), args...)
        if is2d(proxyscene)
            ax = Axis(f; axis...)
        else
            ax = LScene(f; scenekw = (camera = cam3d!, axis...))
        end
    end

    gp[] = ax
    p = plot!(P, ax, args...; kwargs...)
    AxisPlot(ax, p)
end

function plot!(P::PlotFunc, gp::GridPosition, args...; kwargs...)

    c = contents(gp, exact = true)
    if !(length(c) == 1 && c[1] isa Union{Axis, LScene})
        error("There needs to be a single axis at $(gp.span), $(gp.side) to plot into.\nUse a non-mutating plotting command to create an axis implicitly.")
    end
    ax = first(c)
    plot!(P, ax, args...; kwargs...)
end

function plot(P::PlotFunc, gsp::GridSubposition, args...; axis = NamedTuple(), kwargs...)

    layout = GridLayoutBase.get_layout_at!(gsp.parent, createmissing = true)
    c = contents(gsp, exact = true)
    if !isempty(c)
        error("""
        You have used the non-mutating plotting syntax with a GridSubposition, which requires an empty GridLayout slot to create an axis in, but there are already the following objects at this layout position:

        $(c)

        If you meant to plot into an axis at this position, use the plotting function with `!` (e.g. `func!` instead of `func`).
        If you really want to place an axis on top of other layoutables, make your intention clear and create it manually.
        """)
    end

    fig = MakieLayout.get_top_parent(gsp)

    axis = Dict(pairs(axis))

    if haskey(axis, :type)
        axtype = axis[:type]
        pop!(axis, :type)
        ax = axtype(fig; axis...)
    else
        proxyscene = Scene()
        plot!(proxyscene, P, Attributes(kwargs), args...)

        if is2d(proxyscene)
            ax = Axis(fig; axis...)
        else
            ax = LScene(fig; scenekw = (camera = automatic, axis...))
        end
    end

    gsp.parent[gsp.rows, gsp.cols, gsp.side] = ax
    p = plot!(P, ax, args...; kwargs...)
    AxisPlot(ax, p)
end

function plot!(P::PlotFunc, gsp::GridSubposition, args...; kwargs...)

    layout = GridLayoutBase.get_layout_at!(gsp.parent, createmissing = false)

    gp = layout[gsp.rows, gsp.cols, gsp.side]

    c = contents(gp, exact = true)
    if !(length(c) == 1 && c[1] isa Union{Axis, LScene})
        error("There is not just one axis at $(gp).")
    end
    ax = first(c)
    plot!(P, ax, args...; kwargs...)
end


text:
mutable struct CompleteBinaryTree{T} <: AbstractTree{T}
    data::Vector{T}
    size::Int
    function CompleteBinaryTree{T}(size::Int) where T
        data = Vector{T}(undef, size)
        new(data, size)
    end
    function CompleteBinaryTree{T}(data) where T
        new(copy(data), length(data))
    end
end

CompleteBinaryTree(data::Vector) = CompleteBinaryTree{eltype(data)}(data)

# TODO add function such that for t(pos) = t.data[pos]

# Note: It is assumed that pos ∈ [1,tree.size]. Furthermore, if the result
# is out of bounds, zero is returned.

parent(t::CompleteBinaryTree, pos) = pos ÷ 2
leftchild(t::CompleteBinaryTree, pos) = 2*pos <= t.size ? 2*pos : 0
rightchild(t::CompleteBinaryTree, pos) = 2*pos + 1 <= t.size ? 2*pos + 1 : 0
leftsibling(t::CompleteBinaryTree, pos) = pos&1 == 1 ? pos - 1 : 0
rightsibling(t::CompleteBinaryTree, pos) = (pos&1 == 0) && (pos + 1 <= t.size) ? pos + 1 : 0
isleaf(t::CompleteBinaryTree, pos) = (pos > t.size ÷ 2) && (pos <= t.size)

swap!(t::CompleteBinaryTree, i, j) = swap!(t.data, i, j)

# TODO this breaks interface of array-based list, where append! is used
function append!(t::CompleteBinaryTree, item)
    if length(t.data) == t.size
        t.size += 1
        push!(t.data, item)
    else
        t.size += 1
        t.data[t.size] = item
    end
end


text:
# we write a function that solves Ux = y for x where U is upper triangular and y is a vector
function Usolve(U, y)

    # we get the sizes of U and y
    URows, UCols = size(U)
    yRows = size(y)[1]

    # we check that U and y ahve appropriate sizes
    if URows != UCols || yRows != URows
        error("The size of U or y is not appropriate")
    end

    # we initialize x to y
    x = copy(y);

    # we loop backwards over the rows of U
    for k = URows:-1:1

        # we check that we don't divide by zero
        if U[k,k] == 0
            error("There are zeros on the diagonal of U")
        end

        # we modify x accordingly
        for j = k+1:URows
            x[k] -= U[k,j] * x[j];
        end
        x[k] /= U[k,k];
    end

    # we return x
    return x
end

max_stars_repo_path: 1a/LU Decomposition/Usolve.jl

max_issues_repo_path: 1a/LU Decomposition/Usolve.jl

max_forks_repo_path: 1a/LU Decomposition/Usolve.jl


text:
include("polyagamma/je_polyagamma.jl")
include("mmle/je_mmle.jl")

function joint_estimate!(
    items::Vector{<:AbstractItem},
    examinees::Vector{<:AbstractExaminee},
    responses::Vector{Response};
    method = "mmle",
    kwargs...
    )
    if method == "pg"
        return joint_estimate_pg!(items, examinees, responses; kwargs...)
    elseif method =="mmle"
        return joint_estimate_mmle!(items, examinees, responses; kwargs...)
    end
end


text:
using Rimu.Hamiltonians: i_to_k, shift_lattice

"""
    MomPotentialFunction(v, potential)

Function that when combined with [`FullOneBodyTerm`](@ref) acts as a potential in momentum
space.

# Arguments

* `v`: a function or callable object that returns the strength of the potential for given
  component.
* `potential`: a `SVector` of potential values. The `potential` passed to the constructor
  must already be Fourier transformed!
"""
struct MomPotentialFunction{M,V}
    v::V
    potential::SVector{M,Float64}
end
function (f::MomPotentialFunction{M})(σ, p, q) where {M}
    return f.v(σ) * f.potential[mod(p - q, M) + 1]
end

"""
    RealPotentialFunction(v, potential)

Function that when combined with [`ParticleCountTerm`](@ref) acts as a potential in real
space.

# Arguments

* `v`: a function or callable object that returns the strength of the potential for given
  component.
* `potential`: a `SVector` of potential values.
"""
struct RealPotentialFunction{M,V}
    v::V
    potential::SVector{M,Float64}
end
function (f::RealPotentialFunction)(σ, i)
    return f.v(σ) * f.potential[i]
end

"""
    HarmonicPotential(v_ho)

Harmonic potential with strength `v_ho`. The argument `v_ho` can be a real number or a
vector of real numbers whose length matches the number of components in the Hamiltonian.

Can be used with both real and momentum space Hamiltonians.
"""
struct HarmonicPotential{V} <: ExtensionPrototype
    v_ho::V
end

function initialize(hp::HarmonicPotential, ham)
    address = starting_address(ham)
    M = num_modes(address)
    v_ho = ParameterColumnFunction(address, hp.v_ho)

    if basis(ham) ≡ MomentumSpace()
        fun = MomPotentialFunction(v_ho, momentum_space_harmonic_potential(M, 1))
        return FullOneBodyTerm(fun)
    elseif basis(ham) ≡ RealSpace()
        is = range(-fld(M, 2); length=M) # [-M÷2, M÷2) including left boundary
        js = shift_lattice(is) # shifted such that js[1] = 0
        potential = SVector{M}(float.(js).^2)
        fun = RealPotentialFunction(v_ho, potential)
        return ParticleCountTerm(fun)
    else
        throw(ArgumentError(
            "Basis `$(basis(ham))` is not compatible with `HarmonicPotential`"
        ))
    end
end

"""
    DeltaPotential(v)

Delta function potential with strength `v`. The argument `v` can be a real number or a
vector of real numbers whose length matches the number of components in the Hamiltonian.

Can be used with both real and momentum space Hamiltonians.
"""
struct DeltaPotential{V} <: ExtensionPrototype
    v::V
end

function initialize(dp::DeltaPotential, ham)
    address = starting_address(ham)
    M = num_modes(address)
    v = to_parameter_vector(address, dp.v)
    if basis(ham) ≡ MomentumSpace()
        term = FullOneBodyTerm(ParameterColumnFunction(address, v ./ M))
    else
        throw(ArgumentError(
            "Basis `$(basis(ham))` is not compatible with `DeltaPotential`"
        ))
    end
end

"""
    TCDeltaFunctionOneBody(M, cutoff, v, t)

The function used in the one-body term of the transcorrelated delta function potential.
"""
struct TCDeltaFunctionOneBody{M,V}
    corr::CorrelationFactor{M}
    w::WFunction{M}
    v::V
    t::Float64
end

function TCDeltaFunctionOneBody(M::Int, cutoff::Int, v, t)
    corr = CorrelationFactor(M, cutoff)
    w = WFunction(M, cutoff)
    return TCDeltaFunctionOneBody(corr, w, v, t)
end

function (d::TCDeltaFunctionOneBody{M})(σ, p, q) where {M}
    t = d.t
    v = d.v(σ)
    k = q - p
    k_pi = n_to_k(k, M)
    p_pi = i_to_k(p, M)
    cor_k = d.corr(k)

    # note: the factors of 2 and 4 come from the fact that ṽ does not have a factor of
    # 2 in the denominator.
    return v/M + 2v/M * (cor_k * k_pi + 2 * cor_k * p_pi) + 4*v^2/t * d.w(k)
end

"""
    TCDeltaFunctionTwoBody(M, cutoff, v, t)

The function used in the two-body term of the transcorrelated delta function potential.
"""
struct TCDeltaFunctionTwoBody{M,V}
    corr::CorrelationFactor{M}
    v::V
    u::Float64
    t::Float64
end
function TCDeltaFunctionTwoBody(M, cutoff, v::V, u, t) where {V}
    TCDeltaFunctionTwoBody(CorrelationFactor(M, cutoff; length=2M), v, u, t)
end

function (d::TCDeltaFunctionTwoBody{M})(σ, τ, s, r, q, p) where {M}
    u = d.u; v = d.v; t = d.t; corr = d.corr
    result_στ = (v(σ) * corr(p + q - r - s) * corr(r - q))
    result_τσ = (v(τ) * corr(p + q - r - s) * corr(s - p))
    return 16u/(t*M^2) * (result_στ + result_τσ)
end

"""
    TranscorrelatedDeltaPotential(v; cutoff, two_body_term=false)

Transcorrelated delta function potential with strength `v`. The argument `v` can be a real
number or a vector of real numbers whose length matches the number of components in the
Hamiltonian.

The `cutoff` defaults to the cutoff of the Hamiltonian is attached to, or 1.

Can only be used with momentum space Hamiltonians.
"""
struct TranscorrelatedDeltaPotential{V,C} <: ExtensionPrototype
    v::V
    cutoff::C
    two_body_term::Bool
end
function TranscorrelatedDeltaPotential(v; cutoff=nothing, two_body_term=false)
    return TranscorrelatedDeltaPotential(v, cutoff, two_body_term)
end

function initialize(tcd::TranscorrelatedDeltaPotential, ham)
    if basis(ham) ≡ MomentumSpace()
        address = starting_address(ham)
        M = num_modes(address)
        t = only(ham.t)
        u = only(ham.u)
        if ham isa Transcorrelated
            cutoff = isnothing(tcd.cutoff) ? ham.cutoff : tcd.cutoff
        else
            cutoff = isnothing(tcd.cutoff) ? 1 : tcd.cutoff
        end
        v = ParameterColumnFunction(address, tcd.v)

        term = FullOneBodyTerm(TCDeltaFunctionOneBody(M, cutoff, v, t))
        if tcd.two_body_term
            return term + FullTwoBodyTerm(TCDeltaFunctionTwoBody(M, cutoff, v, u, t))
        else
            return term
        end
    else
        throw(ArgumentError(
            "Basis `$(basis(ham))` is not compatible with `TranscorrelatedDeltaPotential`"
        ))
    end
end


text:
###########
# Futures #
###########

# Every future does have a location assigned to it: both a source and a
# location.

"""
    Future()
    Future(value::Any)
    Future(location::Location)
    Future(; kwargs...)

Constructs a new future, representing a value that has not yet been evaluated.
"""
function Future(;source::Location = None(), mutate_from::Union{<:AbstractFuture,Nothing}=nothing)
    # Generate new value id
    value_id = generate_value_id()

    # Create new Future and assign a location to it
    new_future = Future(nothing, value_id, false, true)
    sourced(new_future, source)
    destined(new_future, None())

    # TODO: Add Size location here if needed
    # Handle locations that have an associated value
    if source.src_name in ["None", "Client", "Value"]
        new_future.value = source.sample.value
        new_future.stale = false
    end

    if !isnothing(mutate_from)
        # Indicate that this future is the result of an in-place mutation of
        # some other value
        mutated(mutate_from, new_future)
    elseif source.src_name == "None"
        # For convenience, if a future is constructed with no location to
        # split from, we assume it will be mutated in the next code region
        # and mark it as mutated. This is pretty common since often when
        # we are creating new futures with None location it is as an
        # intermediate variable to store the result of some code region.
        #
        # Mutation can also be specified manually with mutate=true|false in
        # `partition` or implicitly through `Future` constructors
        mutated(new_future)
    end

    new_future
end

function Future(value::Any)
    location = if Base.summarysize(value) ≤ 4 * 1024
        Value(value)
    else
        Client(value)
    end

    # Create future, store value, and return
    Future(source=location)
end

"""
    Future(future::AbstractFuture)

Constructs a future from a future that was already created.

If the given future has not had its value mutated (meaning that the value
stored with it here on the client is the most up-to-date version of it), we use
its value to construct a new future from a copy of the value.

However, if the future has been mutated by some code region that has already
been recorded, we construct a new future with location `None` and mark it as
mutated. This is because presumably in the case that we _can't_ copy over the
given future, we would want to assign to it in the upcoming code region where
it's going to be used.
"""
function Future(fut::AbstractFuture; mutation::Function=identity)
    fut = convert(Future, fut)
    if !fut.stale
        # Copy over value
        new_future = Future(
            deepcopy(mutation(fut.value)),
            generate_value_id(),
            # If the future is not stale, it is not mutated in a way where
            # a further `compute` is needed. So we can just copy its value.
            false,
            false
        )

        # Copy over location
        located(new_future, deepcopy(get_location(fut)))

        new_future
    else
        Future()
    end
end

# convert(::Type{Future}, value::Any) = Future(value)
convert(::Type{Future}, fut::Future) = fut

get_location(value_id::ValueId) = get(get_job().locations, value_id, nothing)
get_location(fut::AbstractFuture) = get_location(convert(Future, fut).value_id)

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
"""
Storing the intermediate results  - like true positives
... and other constants that will be used to given metricks
"""
module  IntermediateData

using Parameters

tn += min(y1,y2)
tp += min(x1,x2)
fn += x1>x2 ? x1-x2 : 0
fp += x2>x1 ? x2-x1 : 0

end#IntermediateData


text:
module FrameFun

using RecipesBase, FillArrays, Reexport
using DomainSets, GridArrays
@reexport using BasisFunctions
# Submodules
include("FrameFunDomains/FrameFunDomains.jl")
@reexport using .FrameFunDomains
include("generic/Platforms/Platforms.jl")
@reexport using .Platforms
include("generic/Platforms/ParameterPaths/ParameterPaths.jl")
@reexport using .ParameterPaths

include("ExtensionFrames/ExtensionFrames.jl")
@reexport using .ExtensionFrames

include("generic/ApproximationProblems.jl")
using .ApproximationProblems

include("generic/FrameFunInterface/FrameFunInterface.jl")
@reexport using .FrameFunInterface

# Other platforms
include("platforms/AugmentationPlatforms.jl")
include("platforms/ExtensionFramePlatforms.jl")
include("platforms/WeightedSumPlatforms.jl")
include("platforms/SumPlatforms.jl")
include("platforms/FrameFunPlatforms.jl")
@reexport using .AugmentationPlatforms, .ExtensionFramePlatforms,
    .WeightedSumPlatforms, .SumPlatforms, .FrameFunPlatforms

include("FrameFunInterfaceExtension/FrameFunInterfaceExtension.jl")
using .FrameFunInterfaceExtension

export residual, abserror, maxerror, L2error
include("fun/error.jl")

include("applications/DiffEquations.jl")
@reexport using .DiffEquations
include("applications/WeightedApproximation.jl")
@reexport using .WeightedApproximation
include("applications/high_dimensional.jl")

include("generic/Adaptivity.jl")
@reexport using .Adaptivity

include("recipes.jl")

end # module


text:
using Test, ClassicalOrthogonalPolynomials, BandedMatrices, LinearAlgebra, LazyArrays, ContinuumArrays, LazyBandedMatrices
import ClassicalOrthogonalPolynomials: CholeskyJacobiBands
import LazyArrays: AbstractCachedMatrix
import LazyBandedMatrices: SymTridiagonal

@testset "Basic properties" begin
    # basis
    P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
    x = axes(P,1)
    J = jacobimatrix(P)
    # example weight
    w(x) = (1 - x^2)
    W = Symmetric(P \ (w.(x) .* P))
    # banded cholesky for symmetric-tagged W
    @test cholesky(W).U isa UpperTriangular
    # compute Jacobi matrix via cholesky
    Jchol = cholesky_jacobimatrix(w,P)
    @test Jchol isa SymTridiagonal
    # CholeskyJacobiBands object
    Cbands = CholeskyJacobiBands(W,P)
    @test Cbands isa CholeskyJacobiBands
    @test Cbands isa AbstractCachedMatrix
    @test getindex(Cbands,1,100) == getindex(Cbands,1,1:100)[100]
end

@testset "Comparison with Lanczos and Classical" begin
    @testset "Using Clenshaw for polynomial weights" begin
        @testset "w(x) = x^2*(1-x)" begin
            P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
            x = axes(P,1)
            J = jacobimatrix(P)
            wf(x) = x^2*(1-x)
            # compute Jacobi matrix via cholesky
            Jchol = cholesky_jacobimatrix(wf, P)
            # compute Jacobi matrix via classical recurrence
            Q = Normalized(Jacobi(1,2)[affine(0..1,Inclusion(-1..1)),:])
            Jclass = jacobimatrix(Q)
            # compute Jacobi matrix via Lanczos
            Jlanc = jacobimatrix(LanczosPolynomial(@.(wf.(x)),Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])))
            # Comparison with Lanczos
            @test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
            # Comparison with Classical
            @test Jchol[1:500,1:500] ≈ Jclass[1:500,1:500]
        end

        @testset "w(x) = (1-x^2)" begin
            P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
            x = axes(P,1)
            J = jacobimatrix(P)
            wf(x) = (1-x^2)
            # compute Jacobi matrix via cholesky
            Jchol = cholesky_jacobimatrix(wf, P)
            # compute Jacobi matrix via Lanczos
            Jlanc = jacobimatrix(LanczosPolynomial(@.(wf.(x)),Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])))
            # Comparison with Lanczos
            @test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
        end

        @testset "w(x) = (1-x^4)" begin
            P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
            x = axes(P,1)
            J = jacobimatrix(P)
            wf(x) = (1-x^4)
            # compute Jacobi matrix via cholesky
            Jchol = cholesky_jacobimatrix(wf, P)
            # compute Jacobi matrix via Lanczos
            Jlanc = jacobimatrix(LanczosPolynomial(@.(wf.(x)),Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])))
            # Comparison with Lanczos
            @test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
        end

        @testset "w(x) = (1.014-x^3)" begin
            P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
            x = axes(P,1)
            J = jacobimatrix(P)
            wf(x) = 1.014-x^4
            # compute Jacobi matrix via cholesky
            Jchol = cholesky_jacobimatrix(wf, P)
            # compute Jacobi matrix via Lanczos
            Jlanc = jacobimatrix(LanczosPolynomial(@.(wf.(x)),Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])))
            # Comparison with Lanczos
            @test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
        end
    end

    @testset "Using Clenshaw with exponential weights" begin
        @testset "w(x) = exp(x)" begin
            P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
            x = axes(P,1)
            J = jacobimatrix(P)
            wf(x) = exp(x)
            # compute Jacobi matrix via cholesky
            Jchol = cholesky_jacobimatrix(wf, P)
            # compute Jacobi matrix via Lanczos
            Jlanc = jacobimatrix(LanczosPolynomial(@.(wf.(x)),Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])))
            # Comparison with Lanczos
            @test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
        end

        @testset "w(x) = (1-x)*exp(x)" begin
            P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
            x = axes(P,1)
            J = jacobimatrix(P)
            wf(x) = (1-x)*exp(x)
            # compute Jacobi matrix via cholesky
            Jchol = cholesky_jacobimatrix(wf, P)
            # compute Jacobi matrix via Lanczos
            Jlanc = jacobimatrix(LanczosPolynomial(@.(wf.(x)),Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])))
            # Comparison with Lanczos
            @test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
        end

        @testset "w(x) = (1-x^2)*exp(x^2)" begin
            P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
            x = axes(P,1)
            J = jacobimatrix(P)
            wf(x) = (1-x^2)*exp(x^2)
            # compute Jacobi matrix via cholesky
            Jchol = cholesky_jacobimatrix(wf, P)
            # compute Jacobi matrix via Lanczos
            Jlanc = jacobimatrix(LanczosPolynomial(@.(wf.(x)),Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])))
            # Comparison with Lanczos
            @test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
        end

        @testset "w(x) = x*(1-x^2)*exp(-x^2)" begin
            P = Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])
            x = axes(P,1)
            J = jacobimatrix(P)
            wf(x) = x*(1-x^2)*exp(-x^2)
            # compute Jacobi matrix via cholesky
            Jchol = cholesky_jacobimatrix(wf, P)
            # compute Jacobi matrix via Lanczos
            Jlanc = jacobimatrix(LanczosPolynomial(@.(wf.(x)),Normalized(Legendre()[affine(0..1,Inclusion(-1..1)),:])))
            # Comparison with Lanczos
            @test Jchol[1:500,1:500] ≈ Jlanc[1:500,1:500]
        end
    end
end


text:
struct kNNRegression
    t
    y::Vector
    kernel::Function
    h::Float64
end

function kNNRegression(X::Matrix, y::Vector; kernel::Symbol = :epanechnikov, bandwidth::Real = 1.0, metric::Metric = Euclidean())
    return kNNRegression(KDTree(X, metric), y,SmoothingKernels.kernels[kernel], bandwidth)
end

function StatsBase.predict(m::kNNRegression, x::Matrix, kernel, bandwidth, k::Integer = 1)
    inds, dists = NearestNeighbors.knn(m.t, x, k)
    map(zip(inds,dists)) do z
        i, d = z
        w = bandwidth .* kernel.(bandwidth .* d)
        sum(w .* m.y[i])/sum(w)
    end
end

StatsBase.predict(m::kNNRegression, x::Matrix, k::Integer = 1) = predict(m, x, m.kernel, m.h, k)


text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni
#

"""
    at(r::Ray, t::Float64) :: Point

Compute the point along the ray's path at some distance from the origin.

Return a [`Point`](@ref) object representing the point in 3D space whose distance
from the ray's origin is equal to `t`, measured in units of the length of `Vec.dir`.

See also: [`Ray`](@ref), [`Point`](@ref), [`Vec`](@ref)
"""
at(r::Ray, t::Float64) = r.origin + r.dir * t

##########################################################################################92

function fire_ray(Ocam::OrthogonalCamera, u::Float64, v::Float64)
    origin = Point(-1.0, (1.0 - 2 * u) * Ocam.a, 2 * v -1)
    direction = Vec(1.0, 0.0, 0.0)
    return Ocam.T*Ray(origin, direction, 1.0)
end # fire_ray

function fire_ray(Pcam::PerspectiveCamera, u::Float64, v::Float64)
    origin = Point(-Pcam.d, 0.0, 0.0)
    direction = Vec(Pcam.d, (1.0 - 2 * u) * Pcam.a,  2 * v - 1)
    return Pcam.T*Ray(origin, direction, 1.0)
end # fire_ray

function fire_ray(
                ImTr::ImageTracer,
                col::Int64, row::Int64,
                u_px::Float64=0.5, v_px::Float64=0.5
            )
    u = (col + u_px) / (ImTr.img.width)
    v = 1. - (row + v_px) / (ImTr.img.height)
    return fire_ray(ImTr.cam, u, v)
end

"""
    fire_ray(Ocam::OrthogonalCamera, u::Float64, v::Float64) :: Ray
    fire_ray(Pcam::PerspectiveCamera, u::Float64, v::Float64) :: Ray
    fire_ray(
            ImTr::ImageTracer,
            col::Int64, row::Int64,
            u_px::Float64=0.5, v_px::Float64=0.5
        ) :: Ray

Shoot one light `Ray` through the pixel (`col`, `row`) of `ImTr.img` image.
The parameters (`col`, `row`) are measured in a diffetent way compared to the
`HDRimage` struct: here, the bottom left corner is placed at `(0, 0)`.
The following diagram shows the convenction for an image with dimensions (`w`,`h`):
```ditaa
| (h,0)  (h,1)  (h,2)  ...  (h,w) |
|  ...    ...    ...   ...   ...  |
| (1,0)  (1,1)  (1,2)  ...  (1,w) |
| (0,0)  (0,1)  (0,2)  ...  (0,w) |
```

The optional values `u_px` and `v_px` specify where the ray should cross
the pixel; the convenction for their values are represented in the following
diagram as `(u_px, v_px)`:
```ditaa
(0, 1)                          (1, 1)
    +------------------------------+
    |                              |
    |                              |
    |                              |
    +------------------------------+
(0, 0)                          (1, 0)
```

See also: [`OrthogonalCamera`](@ref), [`PerspectiveCamera`](@ref),
[`Ray`](@ref), [`HDRimage`](@ref), [`ImageTracer`](@ref)
"""
fire_ray

##########################################################################################92

"""
    fire_all_rays!(
            ImTr::ImageTracer,
            func::Function,
            callback::Union{Nothing, Function} = nothing,
            callback_time_s::Float64 = 2.,
            callback_kwargs::String
            )

Shoot several light rays crossing each of the pixels in the `ImTr.img` image.

For each pixel in the `HDRimage` object fire one `Ray`, and pass it
to the function `func`, which must:
- accept a `Ray` as its only parameter
- return a `RGB{Float32}` color instance telling the color to
  assign to that pixel in the image.

If `callback` is not `nothing`, it must be a function accepting at least two
parameters named `col` and `row`.
This function is called periodically during the rendering, and the two mandatory
arguments are the row and column number of the last pixel that has been traced.

_**Pay Attention**_: Both the row and column are increased by one starting
from zero: first the row and then the column.

The time between two consecutive calls to the callback can be tuned using the
parameter `callback_time_s`. Any keyword argument passed to `fire_all_rays`
is passed to the callback.

See also: [`Ray`](@ref), [`HDRimage`](@ref), [`ImageTracer`](@ref)
"""
function fire_all_rays!(
            ImTr::ImageTracer,
            func::Function,
            callback::Union{Nothing, Function} = nothing,
            callback_time_s::Float64 = 2.,
            callback_kwargs::Union{Nothing, String} = nothing
            )
    last_call_time = time()  # use if @elapsed doesn't work propely for our pourpose

    #isnothing(callback) || callback(0, 0)

    # need a command from line, remember must be a squared number
    ImTr.samples_per_side <= 0 ? smp4side = 0 : smp4side = ImTr.samples_per_side
    pcg = ImTr.pcg
    if smp4side > 0
        for row in ImTr.img.height-1:-1:0, col in 0:ImTr.img.width-1
            cum_color = RGB{Float32}(0., 0., 0.)

            for inter_pixel_row in 0:smp4side-1, inter_pixel_col in 0:smp4side-1
                u_pixel = (inter_pixel_col + random(pcg)) / smp4side
                v_pixel = (inter_pixel_row + random(pcg)) / smp4side
                ray = fire_ray(ImTr, col, row, u_pixel, v_pixel)
                cum_color += func(ray)
            end

            set_pixel(ImTr.img,
                      col,
                      row,
                      cum_color * (1.0/smp4side^2)
            )

            # current_time = time() # use if @elapsed doesn't work propely for our pourpose
            t = time() - last_call_time# current_time - last_call_time # use if @elapsed doesn't work propely for our pourpose

            if (callback ≠ nothing) && (t > callback_time_s)
                callback(row, col)
                last_call_time = time() # use if @elapsed doesn't work propely for our pourpose
            end
        end
    else
        for row in ImTr.img.height-1:-1:0, col in 0:ImTr.img.width-1
            ray = fire_ray(ImTr, col, row)
            set_pixel(ImTr.img, col, row, func(ray))

            # current_time = time() # use if @elapsed doesn't work propely for our pourpose
            t = time() - last_call_time# current_time - last_call_time # use if @elapsed doesn't work propely for our pourpose

            if (callback ≠ nothing) && (t > callback_time_s)
                callback(row, col)
                last_call_time = time()    # use if @elapsed doesn't work propely for our pourpose
            end
        end
    end

#=
    for row in ImTr.img.height-1:-1:0, col in 0:ImTr.img.width-1
        cum_color = RGB{Float32}(0., 0., 0.)

        # need a command from line, remember must be a squared number
        ImTr.samples_per_side <= 0 ? smp4side = 0 : smp4side = ImTr.samples_per_side

        #t = @elapsed
        if smp4side > 0
            for inter_pixel_row in 0:smp4side-1, inter_pixel_col in 0:smp4side-1
                u_pixel = (inter_pixel_col + random(ImTr.pcg)) / smp4side
                v_pixel = (inter_pixel_row + random(ImTr.pcg)) / smp4side
                ray = fire_ray(ImTr, col, row, u_pixel, v_pixel)
                cum_color += func(ray)
            end
            set_pixel(ImTr.img, col, row, cum_color * (1.0/smp4side^2) )
        else
            ray = fire_ray(ImTr, col, row)
            set_pixel(ImTr.img, col, row, func(ray))
        end

        current_time = time() # use if @elapsed doesn't work propely for our pourpose
        t = current_time - last_call_time # use if @elapsed doesn't work propely for our pourpose

        if (callback ≠ nothing) && (t > callback_time_s)
            callback(row, col)
            last_call_time = current_time    # use if @elapsed doesn't work propely for our pourpose
        end

    end
=#

    nothing
end # fire_all_rays!


text:
# This is the Kuncir (1961) adaptive simpson algorithm.
function simpson(a::T, b::T, f::Function, I_of_f::T, i::T, eps::T, N::Int64) where {T <: Real}

# value a, b, eps, N;
# integer N;
# real a, b, I_of_f, i, eps;
# real procedure f;
# comment This procedure integrates the function f(x) using a modified
# Simpson's Rule quadrature formula. The quadrature is performed over j
# subintervals of [a,b] forming the total area I_of_f. Convergence in each
# subinterval of length (b-a)/2^n is indicated when the relative difference
# between successive three-point and five-point area approximations
#  A_{3,j} = (b-a)(g_o + 4g_2 + g_4)/(3*2^{n+1})
#  A_{5,j} = (b-a)(g_o + 4g_1 + 2g_2 + 4g_3 + g_4)/(3*2{n+2})
# is less than or equal to an appropriate portion of the over-all tolerance
# eps (i.e., |(A_{5,j} - A_{3,j})/A_{5,j}| \le  \eps/2^n with n \le N).
# SIMPSON will reduce the size of each interval until this condition is satisfied.
# Complete integration over [a,b] is indicated by i = b. A value
# a =< i < b indicates that the integration was terminated, leaving I_of_f the true
# area under f in [a,i]. Further integration over [i,b] will necessitate either
# the assignment of a larger N, a larger eps, or an integral substitution reducing
# the slope of the integrand in that interval. It is recommended that this
# procedure be used between known integrand maxima and minima.

g=Array{Float64}(5)
A=Array{Float64}(3)
S=Array{Float64}(N,3)
I_of_f=zero(T)
i=zero(T)
m=0; n=0
g[1] = f(a)
g[3] = f(0.5*(a+b))
g[5] = f(b)
# println("g[1]: ",g[1]," g[3]: ",g[3]," g[5]: ",g[5])
bma = b-a
A[1] = 0.5*bma*(g[1]+4*g[3]+g[5])
@label AA
d = 2^n
h = 0.25*bma/d
g[2] = f(a+h*(4*m+1))
g[4] = f(a+h*(4*m+3))
A[2] = h*(g[1]+4*g[2]+g[3])
A[3] = h*(g[3]+4*g[4]+g[5])
#if abs(((A[2]+A[3])-A[1])/(A[2]+A[3])) > eps/d
#if abs((A[2]+A[3])-A[1]) > eps/d
if abs((A[2]+A[3])-A[1]) > 3*eps
#  println("A[2]+A[3])-A[1]: ",A[2]+A[3]-A[1]," 3*eps: ",3*eps," d: ",d)
  m *=2
  n +=1
  if n > N
    @goto CC
  end
  A[1] = A[2]
  S[n,1] = A[3]
  S[n,2] = g[4]
  S[n,3] = g[5]
  g[5] = g[3]
  g[3] = g[2]
  @goto AA
else
  I_of_f += (A[2]+A[3])/3
#  println("I(f): ",I_of_f)
  m += 1
  i = a+m*bma/d
  @label BB
  if iseven(m)
    m = div(m,2)
    n -= 1
    @goto BB
  end
  if (m != 1) || (n != 0)
    A[1] = S[n,1]
    g[1] = g[5]
    g[3] = S[n,2]
    g[5] = S[n,3]
    @goto AA
  end
end
@label CC
#println("I(f): ",I_of_f)
return I_of_f
end