text:
@testset "1260.shift-2-d-grid.jl" begin
    A = [1 2 3; 4 5 6; 7 8 9]
    @test shift_grid(A, 1) == [9 1 2; 3 4 5; 6 7 8]

    A = [
        3 8 1 9
        19 7 2 5
        4 6 11 10
        12 0 21 13
    ]

    A′ = [
        12 0 21 13
        3 8 1 9
        19 7 2 5
        4 6 11 10
    ]
    @test shift_grid(A, 4) == A′
end

ext: jl


text:
using POMDPModels
using Test
using POMDPTesting

let
    problem = LegacyGridWorld()

    policy = RandomPolicy(problem)

    sim = HistoryRecorder(rng=MersenneTwister(1), max_steps=1000)

    hist = simulate(sim, problem, policy, GridWorldState(1,1))

    for (s, a) in zip(state_hist(hist), action_hist(hist))
        td = transition(problem, s, a)
        @test sum(td.probs) ≈ 1.0 atol=0.01
        for p in td.probs
            @test p >= 0.0
        end
    end

    sv = convert_s(Array{Float64}, GridWorldState(1, 1, false), problem)
    @test sv == [1.0, 1.0, 0.0]
    sv = convert_s(Array{Float64}, GridWorldState(5, 3, false), problem)
    @test sv == [5.0, 3.0, 0.0]
    s = convert_s(GridWorldState, sv, problem)
    @test s == GridWorldState(5, 3, false)

    av = convert_a(Array{Float64}, :up, problem)
    @test av == [0.0]
    a = convert_a(Symbol, av, problem)
    @test a == :up

    @test GridWorldState(1,1,false) == GridWorldState(1,1,false)
    @test hash(GridWorldState(1,1,false)) == hash(GridWorldState(1,1,false))
    @test GridWorldState(1,2,false) != GridWorldState(1,1,false)
    @test GridWorldState(1,2,true) == GridWorldState(1,1,true)
    @test hash(GridWorldState(1,2,true)) == hash(GridWorldState(1,1,true))

    POMDPTesting.trans_prob_consistency_check(problem)

end


text:
module DijkstraIndexing

export @dsindex

rewrite(x::Number) = x;
rewrite(x::Symbol) = x;
rewrite(x::String) = x;
function rewrite(expr)
    if expr.head == :ref
        rewr_list = []
        for arg in expr.args[2:end]
            if arg == :(:)
                push!(rewr_list, :(:));
            else
                push!(rewr_list, rewrite(:(1 + $arg)));
            end
        end
        return Expr(expr.head, expr.args[1], rewr_list...);
    elseif expr.head == :(:)
        argl = [rewrite(:(0+$(expr.args[1])))];
        if length(expr.args) == 3
            # Step size.
            push!(argl, rewrite(:(0+$(expr.args[2]))));
        end
        push!(argl, rewrite(:($(expr.args[end])-1)));
        return Expr(expr.head, argl...);
    else
        return Expr(expr.head, [rewrite(i) for i in expr.args]...)
    end
end

macro dsindex(expr)
    return esc(rewrite(expr))
end

end # module


text:
using ClimaTimeSteppers, DiffEqBase, Test

using LinearAlgebra

include("ode_tests_common.jl")

const ArrayType = Array

a = 100
b = 1
c = 1 / 100
Δ = sqrt(4 * a * c - b^2)

α1, α2 = 1 / 4, 3 / 4
β1, β2, β3 = 1 / 3, 3 / 6, 1 / 6

function rhs!(du, u, param, t, α=true, β=false)
    @. du = $(α * cos(t)) * (a + b * u + c * u^2) + β * du
end

function rhs_linear!(dQ, Q, param, t; increment=false)
    if increment
        @. dQ += $cos(t) * b * Q
    else
        @. dQ = $cos(t) * b * Q
    end
end

#=
struct ODETestBasicLinBE <: AbstractBackwardEulerSolver end
ODESolvers.Δt_is_adjustable(::ODETestBasicLinBE) = true
(::ODETestBasicLinBE)(Q, Qhat, α, p, t) = @. Q = Qhat / (1 - α * $cos(t) * b)
=#

function rhs_nonlinear!(dQ, Q, param, t; increment=false)
    if increment
        @. dQ += $cos(t) * (a + c * Q^2)
    else
        @. dQ = $cos(t) * (a + c * Q^2)
    end
end

#=
function rhs_fast!(dQ, Q, ::Any, t; increment)
    if increment
        @. dQ += α1 * $cos(t) * (a + b * Q + c * Q^2)
    else
        @. dQ = α1 * $cos(t) * (a + b * Q + c * Q^2)
    end
end
function rhs_fast_linear!(dQ, Q, ::Nothing, t; increment)
    if increment
        @. dQ += α1 * $cos(t) * Q
    else
        @. dQ = α1 * $cos(t) * Q
    end
end
function rhs_slow!(dQ, Q, ::Nothing, t; increment)
    if increment
        @. dQ += α2 * $cos(t) * (a + b * Q + c * Q^2)
    else
        @. dQ = α2 * $cos(t) * (a + b * Q + c * Q^2)
    end
end
function rhs1!(dQ, Q, ::Any, t; increment)
    if increment
        @. dQ += β1 * $cos(t) * (a + b * Q + c * Q^2)
    else
        @. dQ = β1 * $cos(t) * (a + b * Q + c * Q^2)
    end
end
function rhs2!(dQ, Q, ::Any, t; increment)
    if increment
        @. dQ += β2 * $cos(t) * (a + b * Q + c * Q^2)
    else
        @. dQ = β2 * $cos(t) * (a + b * Q + c * Q^2)
    end
end
function rhs3!(dQ, Q, ::Nothing, t; increment)
    if increment
        @. dQ += β3 * $cos(t) * (a + b * Q + c * Q^2)
    else
        @. dQ = β3 * $cos(t) * (a + b * Q + c * Q^2)
    end
end
=#

function exactsolution(t, q0, t0)
    k = @. 2 * atan((2 * c * q0 + b) / Δ) / Δ - sin(t0)
    solution = @. (Δ * tan((k + sin(t)) * Δ / 2) - b) / (2 * c)
    return ArrayType(solution)
end

q0 = ArrayType === Array ? [1.0] : range(-1.0, 1.0, length = 303)
t0 = 0.1
finaltime = 1.1

Qinit = exactsolution(t0, q0, t0)
Q = similar(Qinit)
Qexact = exactsolution(finaltime, q0, t0)

@testset "Convergence/limited" begin
#=
    @testset "Explicit methods" begin
        dts = [2.0^(-k) for k in 3:4]
        errors = similar(dts)
        for (method, expected_order) in explicit_methods
            for (n, dt) in enumerate(dts)
                Q .= Qinit
                prob = IncrementingODEProblem(rhs!, Q, (t0, finaltime))
                solve(prob, method; dt=dt, adjustfinal=true)
                errors[n] = norm(Q - Qexact)
            end
            rates = log2.(errors[1:(end - 1)] ./ errors[2:end])
            @test isapprox(rates[end], expected_order; atol = 0.7)
        end
    end
    =#

    @testset "IMEX methods - split" begin
        dts = [2.0^(-k) for k in 5:10]
        errors = similar(dts)
        for (method, order) in imex_methods
            #for split_explicit_implicit in (false, true)
                #for variant in (LowStorageVariant(), NaiveVariant())
                    for (n, dt) in enumerate(dts)
                        Q .= Qinit
                        # rhs_arg! =
                        #   split_explicit_implicit ? rhs_nonlinear! : rhs!
                        prob = SplitODEProblem(rhs_linear!, rhs_nonlinear!, Q, (t0, finaltime))
                        solve(prob, method(DirectSolver); dt=dt, adjustfinal=true)
                        errors[n] = norm(Q - Qexact)
                        @show (log2(dt), norm(Q - Qexact))
                    end
                    rates = log2.(errors[1:(end - 1)] ./ errors[2:end])
                    @show rates
                    #=if variant isa LowStorageVariant && split_explicit_implicit
                        if method === ARK1ForwardBackwardEuler
                            expected_order = 1
                        else
                            expected_order = 2
                        end
                    else
                        expected_order = order
                    end=#
                    @test isapprox(rates[end], order; rtol = 0.3)
                #end
            #end
        end
    end
    @testset "IMEX methods - full" begin
        dts = [2.0^(-k) for k in 5:10]
        errors = similar(dts)
        for (method, order) in imex_methods
            #for split_explicit_implicit in (false, true)
                #for variant in (LowStorageVariant(), NaiveVariant())
                    for (n, dt) in enumerate(dts)
                        Q .= Qinit
                        # rhs_arg! =
                        #   split_explicit_implicit ? rhs_nonlinear! : rhs!
                        prob = ODEProblem(ODEFunction(rhs!, jvp=rhs_linear!), Q, (t0, finaltime))
                        solve(prob, method(DirectSolver); dt=dt, adjustfinal=true)
                        errors[n] = norm(Q - Qexact)
                        @show (log2(dt), norm(Q - Qexact))
                    end
                    rates = log2.(errors[1:(end - 1)] ./ errors[2:end])
                    @show rates
                    #=if variant isa LowStorageVariant && split_explicit_implicit
                        if method === ARK1ForwardBackwardEuler
                            expected_order = 1
                        else
                            expected_order = 2
                        end
                    else
                        expected_order = order
                    end=#
                    @test isapprox(rates[end], order; rtol = 0.3)
                #end
            #end
        end
    end

    # @testset "IMEX methods with direct solver" begin
    #     dts = [2.0^(-k) for k in 4:5]
    #     errors = similar(dts)
    #     for (method, order) in imex_methods
    #         for split_explicit_implicit in (false, true)
    #             for (n, dt) in enumerate(dts)
    #                 Q .= Qinit
    #                 rhs_arg! = split_explicit_implicit ? rhs_nonlinear! : rhs!
    #                 solver = method(
    #                     rhs_arg!,
    #                     rhs_linear!,
    #                     ODETestBasicLinBE(),
    #                     Q;
    #                     dt = dt,
    #                     t0 = t0,
    #                     split_explicit_implicit = split_explicit_implicit,
    #                     variant = NaiveVariant(),
    #                 )
    #                 solve!(Q, solver; timeend = finaltime)
    #                 errors[n] = norm(Q - Qexact)
    #             end
    #             rates = log2.(errors[1:(end - 1)] ./ errors[2:end])
    #             expected_order = order
    #             @test isapprox(rates[end], expected_order; rtol = 0.3)
    #         end
    #     end
    # end

    # @testset "MRRK methods with 2 rates" begin
    #     dts = [2.0^(-k) for k in 3:4]
    #     errors = similar(dts)
    #     for (slow_method, slow_expected_order) in slow_mrrk_methods
    #         for (fast_method, fast_expected_order) in fast_mrrk_methods
    #             for nsubsteps in (1, 3)
    #                 for (n, dt) in enumerate(dts)
    #                     Q .= Qinit
    #                     solver = Multirate(
    #                         (
    #                             slow_method(rhs_slow!, Q; dt = dt),
    #                             fast_method(rhs_fast!, Q; dt = dt / nsubsteps),
    #                         );
    #                         t0 = t0,
    #                     )
    #                     solve!(Q, solver; timeend = finaltime)
    #                     errors[n] = norm(Q - Qexact)
    #                 end

    #                 rates = log2.(errors[1:(end - 1)] ./ errors[2:end])
    #                 min_order = min(slow_expected_order, fast_expected_order)
    #                 max_order = max(slow_expected_order, fast_expected_order)
    #                 @test (
    #                     isapprox(rates[end], min_order; atol = 0.5) ||
    #                     isapprox(rates[end], max_order; atol = 0.5) ||
    #                     min_order <= rates[end] <= max_order
    #                 )
    #             end
    #         end
    #     end
    # end

    # @testset "MRRK methods with IMEX" begin
    #     dts = [2.0^(-k) for k in 3:4]
    #     errors = similar(dts)
    #     for (slow_method, slow_expected_order) in slow_mrrk_methods
    #         for (fast_method, fast_expected_order) in imex_methods
    #             for (n, dt) in enumerate(dts)
    #                 Q .= Qinit
    #                 solver = Multirate(
    #                     (
    #                         slow_method(rhs_slow!, Q; dt = dt),
    #                         fast_method(
    #                             rhs_fast!,
    #                             rhs_fast_linear!,
    #                             LinearBackwardEulerSolver(
    #                                 DivideLinearSolver();
    #                                 isadjustable = true,
    #                             ),
    #                             Q;
    #                             dt = dt,
    #                             split_explicit_implicit = false,
    #                         ),
    #                     ),
    #                     t0 = t0,
    #                 )
    #                 solve!(Q, solver; timeend = finaltime)
    #                 errors[n] = norm(Q - Qexact)
    #             end
    #             rates = log2.(errors[1:(end - 1)] ./ errors[2:end])
    #             min_order = min(slow_expected_order, fast_expected_order)
    #             max_order = max(slow_expected_order, fast_expected_order)
    #             @test (
    #                 isapprox(rates[end], min_order; atol = 0.5) ||
    #                 isapprox(rates[end], max_order; atol = 0.5) ||
    #                 min_order <= rates[end] <= max_order
    #             )
    #         end
    #     end
    # end

    # @testset "MRRK methods with 3 rates" begin
    #     dts = [2.0^(-k) for k in 3:4]
    #     errors = similar(dts)
    #     for (rate3_method, rate3_order) in slow_mrrk_methods
    #         for (rate2_method, rate2_order) in slow_mrrk_methods
    #             for (rate1_method, rate1_order) in fast_mrrk_methods
    #                 for nsubsteps in (1, 2)
    #                     for (n, dt) in enumerate(dts)
    #                         Q .= Qinit
    #                         solver = Multirate(
    #                             (
    #                                 rate3_method(rhs3!, Q, dt = dt),
    #                                 rate2_method(rhs2!, Q, dt = dt / nsubsteps),
    #                                 rate1_method(
    #                                     rhs1!,
    #                                     Q;
    #                                     dt = dt / nsubsteps^2,
    #                                 ),
    #                             );
    #                             dt = dt,
    #                             t0 = t0,
    #                         )
    #                         solve!(Q, solver; timeend = finaltime)
    #                         errors[n] = norm(Q - Qexact)
    #                     end
    #                     rates = log2.(errors[1:(end - 1)] ./ errors[2:end])
    #                     @test 3.8 <= rates[end]
    #                 end
    #             end
    #         end
    #     end
    # end

    # @testset "MIS methods" begin
    #     dts = [2.0^(-k) for k in 3:4]
    #     errors = similar(dts)
    #     for (method, expected_order) in mis_methods
    #         for fast_method in (LSRK54CarpenterKennedy,)
    #             for (n, dt) in enumerate(dts)
    #                 Q .= Qinit
    #                 solver = method(
    #                     rhs_slow!,
    #                     rhs_fast!,
    #                     fast_method,
    #                     4,
    #                     Q;
    #                     dt = dt,
    #                     t0 = t0,
    #                 )
    #                 solve!(Q, solver; timeend = finaltime)
    #                 errors[n] = norm(Q - Qexact)
    #             end
    #             rates = log2.(errors[1:(end - 1)] ./ errors[2:end])
    #             @test isapprox(rates[end], expected_order; atol = 0.6)
    #         end
    #     end
    # end

    # @testset "MRI GARK methods with 2 rates" begin
    #     dts = [2.0^(-k) for k in 3:4]
    #     errors = similar(dts)
    #     for (slow_method, expected_order) in mrigark_erk_methods
    #         for (fast_method, _) in fast_mrigark_methods
    #             for nsubsteps in (1, 3)
    #                 for (n, dt) in enumerate(dts)
    #                     dt /= 4 # Need a smaller dt to get convergence rate
    #                     Q .= Qinit
    #                     fastsolver =
    #                         fast_method(rhs_fast!, Q; dt = dt / nsubsteps)
    #                     solver = slow_method(
    #                         rhs_slow!,
    #                         fastsolver,
    #                         Q;
    #                         dt = dt,
    #                         t0 = t0,
    #                     )
    #                     solve!(Q, solver; timeend = finaltime)

    #                     errors[n] = norm(Q - Qexact)
    #                 end

    #                 rates = log2.(errors[1:(end - 1)] ./ errors[2:end])
    #                 @test isapprox(rates[end], expected_order; atol = 0.5)
    #             end
    #         end
    #     end
    # end

    # @testset "MRI GARK methods with 3 rates" begin
    #     dts = [2.0^(-k) for k in 3:4]
    #     errors = similar(dts)
    #     for (rate3_method, rate3_order) in mrigark_erk_methods
    #         for (rate2_method, rate2_order) in mrigark_erk_methods
    #             for (rate1_method, _) in fast_mrigark_methods
    #                 for nsubsteps in (1, 2)
    #                     for (n, dt) in enumerate(dts)
    #                         dt /= 4 # Need a smaller dt to get convergence rate
    #                         Q .= Qinit
    #                         solver1 =
    #                             rate1_method(rhs1!, Q; dt = dt / nsubsteps^2)
    #                         solver2 = rate2_method(
    #                             rhs2!,
    #                             solver1,
    #                             Q;
    #                             dt = dt / nsubsteps,
    #                         )
    #                         solver3 = rate3_method(
    #                             rhs3!,
    #                             solver2,
    #                             Q;
    #                             dt = dt,
    #                             t0 = t0,
    #                         )
    #                         solve!(Q, solver3; timeend = finaltime)
    #                         errors[n] = norm(Q - Qexact)
    #                     end
    #                     rates = log2.(errors[1:(end - 1)] ./ errors[2:end])
    #                     expected_order = min(rate3_order, rate2_order)
    #                     @test isapprox(rates[end], expected_order; atol = 0.5)
    #                 end
    #             end
    #         end
    #     end
    # end

    # @testset "MRI GARK implicit methods with 2 rates and linear solver" begin
    #     dts = [2.0^(-k) for k in 3:4]
    #     errors = similar(dts)
    #     for (slow_method, expected_order) in mrigark_irk_methods
    #         for (fast_method, _) in fast_mrigark_methods
    #             for nsubsteps in (1, 3)
    #                 for (n, dt) in enumerate(dts)
    #                     dt /= 4
    #                     Q .= Qinit
    #                     nsteps = ceil(Int, (finaltime - t0) / dt)
    #                     dt = (finaltime - t0) / nsteps
    #                     fastsolver =
    #                         fast_method(rhs_nonlinear!, Q; dt = dt / nsubsteps)
    #                     solver = slow_method(
    #                         rhs_linear!,
    #                         LinearBackwardEulerSolver(
    #                             DivideLinearSolver();
    #                             isadjustable = true,
    #                         ),
    #                         fastsolver,
    #                         Q;
    #                         dt = dt,
    #                         t0 = t0,
    #                     )
    #                     solve!(Q, solver; timeend = finaltime)

    #                     errors[n] = norm(Q - Qexact)
    #                 end

    #                 rates = log2.(errors[1:(end - 1)] ./ errors[2:end])
    #                 @test isapprox(rates[end], expected_order; atol = 0.5)
    #             end
    #         end
    #     end
    # end

    # @testset "MRI GARK implicit methods with 2 rates and custom solver" begin
    #     dts = [2.0^(-k) for k in 3:4]
    #     errors = similar(dts)
    #     for (slow_method, expected_order) in mrigark_irk_methods
    #         for (fast_method, _) in fast_mrigark_methods
    #             for nsubsteps in (1, 3)
    #                 for (n, dt) in enumerate(dts)
    #                     dt /= 4
    #                     Q .= Qinit
    #                     nsteps = ceil(Int, (finaltime - t0) / dt)
    #                     dt = (finaltime - t0) / nsteps
    #                     fastsolver =
    #                         fast_method(rhs_nonlinear!, Q; dt = dt / nsubsteps)
    #                     solver = slow_method(
    #                         rhs_linear!,
    #                         ODETestBasicLinBE(),
    #                         fastsolver,
    #                         Q;
    #                         dt = dt,
    #                         t0 = t0,
    #                     )
    #                     solve!(Q, solver; timeend = finaltime)

    #                     errors[n] = norm(Q - Qexact)
    #                 end

    #                 rates = log2.(errors[1:(end - 1)] ./ errors[2:end])
    #                 @test isapprox(rates[end], expected_order; atol = 0.5)
    #             end
    #         end
    #     end
    # end
end

# @testset "Explicit methods composition of solve!" begin
#     halftime = 1.0
#     finaltime = 2.0
#     dt = 0.075
#     for (method, _) in explicit_methods
#         Q .= Qinit
#         solver1 = method(rhs!, Q; dt = dt, t0 = t0)
#         solve!(Q, solver1; timeend = finaltime)

#         Q2 = similar(Q)
#         Q2 .= Qinit
#         solver2 = method(rhs!, Q2; dt = dt, t0 = t0)
#         solve!(Q2, solver2; timeend = halftime, adjustfinalstep = false)
#         solve!(Q2, solver2; timeend = finaltime)

#         @test Q2 == Q
#     end
# end


text:
module Symbolic

using ..Prim
using ..IDS: ID
using ..Space

export analyticalmean

include("symomega.jl")

end


text:
# single line

#=
Multi-
line
comment
=#

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
n = 2
m = 1

A = [1.0 1.0; 0.0 1.0]
B = [0.0; 1.0]

Q = [1.0 0.0; 0.0 1.0]
R = [1.0]

# x_init = [1.0; 0.0]

nz = n + m + n

function unpack(z)
	x2 = z[1:n]
	u1 = z[n .+ (1:m)]
	y2 = z[n + m .+ (1:n)]

	return x2, u1, y2
end

z0 = rand(nz)

function lagrangian(z, x1)
	x2, u1, y2 = unpack(z)

	L = 0.0

	# objective
	L += transpose(x2) * Q * x2
	L += (transpose(u1) * R * u1)[1]

	# constraints
	L += transpose(y2) * (A * x1 + B * u1[1] - x2)

	return L
end

@variables z_sym[1:nz]
@variables x1_sym[1:n]

L = lagrangian(z_sym, x1_sym)
L = simplify.(L)

dL = Symbolics.gradient(L, z_sym)
ddL = Symbolics.hessian(L, z_sym)
ddLx1 = Symbolics.jacobian(dL, x1_sym)

L_grad = eval(Symbolics.build_function(dL, z_sym, x1_sym)[1])
L_hess = eval(Symbolics.build_function(ddL, z_sym, x1_sym)[1])
L_hess_x1 = eval(Symbolics.build_function(ddLx1, z_sym, x1_sym)[1])

function solve(z0, x1)
	z = copy(z0)

	∇L = L_grad(z, x1)

	println()
	for i = 1:10
		r_norm = norm(∇L)
		println("iter ($i)")
		println("	grad norm: $(r_norm)")
		r_norm < 1.0e-8 && break

		∇²L = L_hess(z, x1)

		Δ = ∇²L \ ∇L

		α = 1.0
		ẑ = z - α * Δ

		∇L̂ = L_grad(ẑ, x1)

		iter = 1
		while norm(∇L̂) > norm(∇L)
			α *= 0.5
			ẑ = z - α * Δ
			∇L̂ = L_grad(ẑ, x1)
			iter += 1

			if iter > 100
				break
			end
		end

		z = ẑ
		∇L = ∇L̂
	end

	return z, ∇L, L_hess(z, x1)
end

z_sol, ∇L_sol, ∇²L_sol = solve(z0, x1)
K_direct = (∇²L_sol \ L_hess_x1(z_sol, x1))[n .+ (1:m), :]

# LQR solution
K, P = tvlqr([A], [B], [Q, Q], [R[1,1]])

norm(K_direct - K[1])

max_issues_repo_name: thowell/DirectMotionPlanning


text:
import Base.depwarn

for (fitfn, fittype) in ((:glm, :GeneralizedLinearModel), (:lm, :LinearModel), (:lmc, :(LinearModel{DensePredChol})))
    @eval begin
        function $fitfn(e::Expr, df, args...; kwargs...)
            depwarn($("$(string(fitfn))(e::Expr, df::AbstractDataFrame, ...) is deprecated, use fit($(string(fittype)), f::Formula, df::AbstractDataFrame, ...) instead"), $(Base.Meta.quot(fitfn)))
            eval(quote; import DataFrames; end)
            fit($fittype, DataFrames.Formula(e), df, args...; kwargs...)
        end

        function $fitfn(s::String, df, args...; kwargs...)
            depwarn($("$(string(fitfn))(e::String, df::AbstractDataFrame, ...) is deprecated, use fit($(string(fittype)), f::Formula, df::AbstractDataFrame, ...) instead"), $(Base.Meta.quot(fitfn)))
            eval(quote; import DataFrames; end)
            fit($fittype, DataFrames.Formula(parse(s)[1]), df, args...; kwargs...)
        end
    end
end

typealias LmMod LinearModel
typealias GlmMod GeneralizedLinearModel


text:
using Test, YaoArrayRegister, BitBasis, LinearAlgebra
using YaoBase
using Adapt

@testset "test ArrayReg constructors" begin
    @test ArrayReg(rand(4)) isa ArrayReg{2}
    @test_logs (:warn, "Input matrix element type is not `Complex`, got `Float64`") ArrayReg(
        rand(4, 4),
    )
    @test ArrayReg(rand(ComplexF64, 4)) isa ArrayReg{2}
    @test ArrayReg(arrayreg(bit"101")) == arrayreg(bit"101")

    st = rand(ComplexF64, 4, 8)
    @test state(adjoint(ArrayReg(st))) == adjoint(st)

    @test state(arrayreg(bit"101")) == reshape(onehot(bit"101"), :, 1)
    @test datatype(arrayreg(ComplexF32, bit"101")) == ComplexF32

    st = rand(ComplexF64, 4, 8)
    reg= ArrayReg{2}(st)
    @test similar(reg) isa ArrayReg
    @test nactive(similar(reg)) == 2
    m = randn(ComplexF64, 8, 8)
    @test similar(reg, m).state == m
    @test viewbatch(reg, 1) == reg
    @test_throws ErrorException viewbatch(reg, 2)
end

@testset "test BatchedArrayReg constructors" begin
    @test BatchedArrayReg(rand(4, 6), 3) isa BatchedArrayReg{2}
    @test_throws DimensionMismatch BatchedArrayReg(rand(4, 3), 2)
    @test_throws DimensionMismatch BatchedArrayReg(rand(5, 2), 2)
    @test BatchedArrayReg(rand(ComplexF64, 4, 3)) isa BatchedArrayReg{2}

    @test state(arrayreg(bit"101"; nbatch=2)) == repeat(reshape(onehot(bit"101"), :, 1), 1, 2)
    @test datatype(arrayreg(ComplexF32, bit"101"; nbatch=2)) == ComplexF32

    st = rand(ComplexF64, 4, 6)
    reg = BatchedArrayReg(st, 6)
    regt = transpose_storage(reg)
    @test regt.state isa Transpose
    @test regt == reg
    @test transpose_storage(regt).state isa Matrix
    @test transpose_storage(regt) == reg
    @test_throws ErrorException ArrayReg(reg)

    reg = zero_state(5; nbatch=1)
    @test reg isa BatchedArrayReg
    r2 = ArrayReg(reg)
    @test r2 isa ArrayReg
    @test BatchedArrayReg(r2) == reg

    st = rand(ComplexF64, 4, 8)
    reg= BatchedArrayReg{2}(st, 4)
    @test similar(reg) isa BatchedArrayReg
    @test nactive(similar(reg)) == 2
    @test nqubits(similar(reg)) == 3
    m = randn(ComplexF64, 8, 8)
    @test similar(reg, m).state == m
end

@testset "test $T initialization methods" for T in [ComplexF64, ComplexF32, ComplexF16]
    @testset "test product state" begin
        reg = product_state(T, bit"100")
        st = state(reg)
        @test nqudits(reg) == 3
        @test !(st isa Transpose)
        st2 = state(product_state(T, [0, 0, 1]))
        @test st2 ≈ st
        st = state(product_state(T, bit"100"; nbatch = 2, no_transpose_storage = true))
        @test !(st isa Transpose)
        st = state(product_state(T, bit"100"; nbatch = 2))
        @test st isa Transpose
        for k = 1:2
            @test st[:, k] ≈ onehot(T, bit"100")
        end

        st = state(product_state(T, 4, 0; nbatch = 3))
        for k = 1:3
            @test st[:, k] ≈ onehot(T, 4, 0)
        end
        @test eltype(product_state(Float64, 4, 0).state) == Float64
    end
    @testset "test zero state" begin
        st = state(zero_state(T, 3))
        @test !(st isa Transpose)
        st = state(zero_state(T, 3; nbatch = 2, no_transpose_storage = true))
        @test !(st isa Transpose)
        st = state(zero_state(T, 4; nbatch = 4))
        @test st isa Transpose
        for k = 1:4
            @test st[:, k] ≈ onehot(T, 4, 0)
        end
        @test eltype(zero_state(Float64, 4).state) == Float64
    end
    @testset "test rand state" begin
        st = state(rand_state(T, 3))
        @test !(st isa Transpose)
        st = state(rand_state(T, 3; nbatch = 2, no_transpose_storage = true))
        @test !(st isa Transpose)
        # NOTE: we only check if the state is normalized
        st = state(rand_state(T, 4, nbatch = 2))
        @test st isa Transpose
        for k = 1:2
            @test norm(st[:, k]) ≈ 1.0
        end
        @test eltype(rand_state(Float64, 4).state) == Float64
    end
    @testset "test uniform state" begin
        st = state(uniform_state(T, 3))
        @test !(st isa Transpose)
        st = state(uniform_state(T, 3; nbatch = 2, no_transpose_storage = true))
        @test !(st isa Transpose)
        st = state(uniform_state(T, 4; nbatch = 2))
        @test st isa Transpose
        for k = 1:2
            for each in st[:, k]
                @test each ≈ 1 / sqrt(16)
            end
        end
        @test eltype(uniform_state(Float64, 4).state) == Float64
    end
    @testset "test oneto" begin
        r1 = uniform_state(ComplexF64, 4)
        r2 = oneto(r1, 2)
        @test nactive(r2) == 2
        @test r1 |> oneto(2) |> nactive == 2
    end
    @testset "test repeat" begin
        r = repeat(arrayreg(T, bit"101"; nbatch=3), 4)
        @test nactive(r) == 3
        @test nbatch(r) == 12
    end
end

@testset "test YaoBase interface" begin
    @testset "test probs" begin
        r = rand_state(5; nbatch = 3)
        @test probs(r) ≈ abs2.(state(r))
        r = rand_state(5)
        @test probs(r) ≈ abs2.(state(r))
    end
    @testset "test batch iteration" begin
        r = arrayreg(bit"101"; nbatch= 3)
        for k = 1:3
            @test viewbatch(r, k) == arrayreg(bit"101")
        end
        # broadcast
        for each in r
            @test each == arrayreg(bit"101")
        end
    end
    @testset "test addbits!" begin
        @test addbits!(zero_state(3), 3) == zero_state(6)
        r = rand_state(3; nbatch = 2)
        @test addbits!(copy(r), 2) ≈ join(zero_state(2; nbatch = 2), r)
        r = rand_state(3; nbatch = 1)
        @test addbits!(copy(r), 2) ≈ join(zero_state(2; nbatch = 1), r)
    end
end

@testset "test statevec" begin
    r = arrayreg(bit"101"; nbatch=3)
    @test statevec(r) == r.state
    r = oneto(r, 2)
    @test statevec(r) == reshape(state(r), 4, 6)
end

@testset "test hypercubic" begin
    @test hypercubic(rand_state(3)) |> size == (2, 2, 2, 1)
end

# TODO: test concat multiple registers
@testset "test join" begin
    r1 = rand_state(6)
    r2 = rand_state(6)
    r3 = join(r2, r1)
    r4 = join(focus!(copy(r2), 1:2), focus!(copy(r1), 1:3))
    @test r4 |> relaxedvec ≈
          focus!(copy(r3), [1, 2, 3, 7, 8, 4, 5, 6, 9, 10, 11, 12]) |> relaxedvec
    reg5 = focus!(repeat(r1, 3), 1:3)
    reg6 = focus!(repeat(r2, 3), 1:2)
    @test (join(reg6, reg5)|>relaxedvec)[:, 1] ≈ r4 |> relaxedvec

    # manual trace
    r = join(arrayreg(bit"011"), zero_state(1))
    focus!(r, 2:4)
    @test sum(r.state, dims = 2) ≈ arrayreg(bit"011").state
end

@testset "YaoBlocks.jl/issues/21" begin
    st = rand(ComplexF64, 16, 2)
    r1 = ArrayReg(view(st, :, 1))
    r2 = ArrayReg(rand(ComplexF64, 16, 1))
    copyto!(r1, r2)
    @test r1 == r2
end

@testset "transpose copy" begin
    reg = rand_state(5; nbatch = 10)
    reg1 = copy(reg)
    reg2 = focus!(copy(reg), (3, 5))
    reg3 = relax!(copy(reg2), (3, 5))
    reg4 = oneto(reg1, 3)
    @test reg1.state isa Transpose
    @test reg2.state isa Transpose
    @test reg3.state isa Transpose
    @test reg4.state isa Transpose
    @test reg4.state ≈ oneto(reg, 3).state
    reg4.state[1] = 2.0
    @test reg.state[1] != 2.0
end

@testset "collapseto" begin
    reg = rand_state(4)
    reg2 = copy(reg)
    focus!(reg, (4, 2))
    collapseto!(reg, bit"01")
    relax!(reg, (4, 2), to_nactive = 4)
    instruct!(reg2, Const.P0, (2,))
    instruct!(reg2, Const.P1, (4,))
    normalize!(reg2)
    @test reg ≈ reg2
end

@testset "overlap between two focused reg" begin
    rega = rand_state(5)
    regb = rand_state(5)
    reg1 = focus!(copy(rega), (3, 2, 4))
    reg2 = focus!(copy(regb), (3, 2, 4))
    @test reg1' * reg2 ≈ rega' * regb
end

@testset "adapt" begin
    @test datatype(adapt(Array{ComplexF32}, zero_state(5))) == ComplexF32
end

@testset "von_neumann_entropy" begin
    reg = (product_state(bit"00000") + product_state(bit"11111")) / sqrt(2)
    @test von_neumann_entropy(reg, [2,3]) ≈ log(2)
    @test mutual_information(reg, [2,3], [5]) ≈ log(2)
end

@testset "qudit" begin
    reg = ArrayReg{3}(reshape(randn(ComplexF64, 9), :, 1))
    @test nlevel(reg) == 3
    @test nactive(reg) == 2
    @test_throws MethodError nqubits(reg)
    @test nqudits(reg) == 2

    # constructors
    reg = zero_state(2; nlevel=3)
    @test statevec(reg) == [1.0, 0.0im, 0, 0, 0, 0, 0, 0, 0]
    @test statevec(product_state(2, 1; nlevel=3)) == [0.0, 1.0+0.0im, 0, 0, 0, 0, 0, 0, 0]
    @test length(statevec(rand_state(2; nlevel=3))) == 9
    @test size(statevec(zero_state(2; nbatch=5, nlevel=3))) == (9,5)
end

@testset "more (push test coverage)" begin
    reg1 = transpose_storage(rand_state(3; nbatch=5))
    reg2 = transpose_storage(rand_state(3; nbatch=5))
    res = reg1' * reg2
    @test res ≈ (reg1') .* reg2

    reg = transpose_storage(focus!(rand_state(6), [2,3]))
    @test reg.state isa Transpose
    @test copy(reg).state isa Transpose

    copyto!(reg1', reg2')
    @test reg1 ≈ reg2
    @test reorder!(reg1, [2,3,1]).state ≈ reshape(permutedims(reshape(reg2.state, 2, 2, 2, 5), sortperm([2,3,1,4])), 8, 5)
end

@testset "most_probable" begin
    reg = arrayreg(ComplexF64[0.0, 0.6, 0.8, 0.0])
    reg2 = arrayreg(ComplexF64[0.0, 0.6, 0.0, 0.8])
    mp = most_probable(reg, 1)
    @test length(mp) == 1 && mp[1] === BitStr{2}(2)
    @test most_probable(reg, 2) == BitStr{2}.([2, 1])
    breg = BatchedArrayReg(reg, reg2)
    @test most_probable(breg, 2) == BitStr{2}.([2 3; 1 1])
end


text:
using RethinkDB
using Base.Test

# write your own tests here
@test 1 == 1

include("ast.jl")

max_stars_repo_name: mrcsparker/julia-rethinkdb

max_issues_repo_name: mrcsparker/julia-rethinkdb

max_forks_repo_name: mrcsparker/julia-rethinkdb


text:
using EchoJulia

filename = EK60_SAMPLE # or some EK60 RAW file name

ps = pings(filename) # Get the pings
ps38 = [p for p in ps if p.frequency == 38000] # Just 38 kHz pings
Sv38 = Sv(ps38) # Convert to a matrix of volume backscatter

# Show a quick echogram

eg(Sv38)

# Full echogram, preserving resolution

egshow(Sv38)

# More like Echoview?

eg(Sv38,cmap=EK500, vmin=-95,vmax=-50)

# Show a histogram

eghist(Sv38)

# Echogram as an RGB image for use with Images.jl and friends

img = imagesc(Sv38)

# Export to MATLAB format for further analysis

using MAT

matwrite("myfile.mat", Dict("Sv38" => Sv38))


text:
using LinearAlgebra

### SubArray
# `offset1` and `stride1` fields are calculated from parent indices.
# Setting them has no effect.
subarray_constructor(parent, indices, args...) = SubArray(parent, indices)

constructorof(::Type{<:SubArray}) = subarray_constructor

### ReinterpretArray
struct ReinterpretArrayConstructor{T} end
# `readable` and `writeable` are calculated from `T` and `parent`
# Setting them has no effect.
function (::ReinterpretArrayConstructor{T})(parent, args...) where T
    reinterpret(T, parent)
end

constructorof(::Type{<:Base.ReinterpretArray{T}}) where T = ReinterpretArrayConstructor{T}()

### PermutedDimsArray
struct PermutedDimsArrayConstructor{N,perm,iperm} end
# Parent must have the same N - it has to match length(perm)
function (::PermutedDimsArrayConstructor{N,perm,iperm})(parent::AA
) where {N,perm,iperm,AA<:AbstractArray{T,N}} where T
    PermutedDimsArray{T,N,perm,iperm,AA}(parent)
end

constructorof(::Type{<:PermutedDimsArray{<:Any,N,perm,iperm,<:Any}}) where {N,perm,iperm} =
    PermutedDimsArrayConstructor{N,perm,iperm}()

### Tridiagonal
function tridiagonal_constructor(dl::V, d::V, du::V) where {V<:AbstractVector{T}} where T
    Tridiagonal{T,V}(dl, d, du)
end
function tridiagonal_constructor(dl::V, d::V, du::V, du2::V) where {V<:AbstractVector{T}} where T
    Tridiagonal{T,V}(dl, d, du, du2)
end

# `du2` may be undefined, so we need a custom `getproperties` that checks `isdefined`
function getproperties(o::Tridiagonal)
    if isdefined(o, :du2)
        (dl=o.dl, d=o.d, du=o.du, du2=o.du2)
    else
        (dl=o.dl, d=o.d, du=o.du)
    end
end

constructorof(::Type{<:LinearAlgebra.Tridiagonal}) = tridiagonal_constructor

### LinRange
# `lendiv` is a calculated field
linrange_constructor(start, stop, len, lendiv) = LinRange(start, stop, len)

constructorof(::Type{<:LinRange}) = linrange_constructor


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

mutable struct PacketCaptureFilter <: SwaggerModel
    protocol::Any # spec type: Union{ Nothing, String } # spec name: protocol
    localIPAddress::Any # spec type: Union{ Nothing, String } # spec name: localIPAddress
    remoteIPAddress::Any # spec type: Union{ Nothing, String } # spec name: remoteIPAddress
    localPort::Any # spec type: Union{ Nothing, String } # spec name: localPort
    remotePort::Any # spec type: Union{ Nothing, String } # spec name: remotePort

    function PacketCaptureFilter(;protocol="Any", localIPAddress=nothing, remoteIPAddress=nothing, localPort=nothing, remotePort=nothing)
        o = new()
        validate_property(PacketCaptureFilter, Symbol("protocol"), protocol)
        setfield!(o, Symbol("protocol"), protocol)
        validate_property(PacketCaptureFilter, Symbol("localIPAddress"), localIPAddress)
        setfield!(o, Symbol("localIPAddress"), localIPAddress)
        validate_property(PacketCaptureFilter, Symbol("remoteIPAddress"), remoteIPAddress)
        setfield!(o, Symbol("remoteIPAddress"), remoteIPAddress)
        validate_property(PacketCaptureFilter, Symbol("localPort"), localPort)
        setfield!(o, Symbol("localPort"), localPort)
        validate_property(PacketCaptureFilter, Symbol("remotePort"), remotePort)
        setfield!(o, Symbol("remotePort"), remotePort)
        o
    end
end # type PacketCaptureFilter

const _property_map_PacketCaptureFilter = Dict{Symbol,Symbol}(Symbol("protocol")=>Symbol("protocol"), Symbol("localIPAddress")=>Symbol("localIPAddress"), Symbol("remoteIPAddress")=>Symbol("remoteIPAddress"), Symbol("localPort")=>Symbol("localPort"), Symbol("remotePort")=>Symbol("remotePort"))
const _property_types_PacketCaptureFilter = Dict{Symbol,String}(Symbol("protocol")=>"String", Symbol("localIPAddress")=>"String", Symbol("remoteIPAddress")=>"String", Symbol("localPort")=>"String", Symbol("remotePort")=>"String")
Base.propertynames(::Type{ PacketCaptureFilter }) = collect(keys(_property_map_PacketCaptureFilter))
Swagger.property_type(::Type{ PacketCaptureFilter }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_PacketCaptureFilter[name]))}
Swagger.field_name(::Type{ PacketCaptureFilter }, property_name::Symbol) =  _property_map_PacketCaptureFilter[property_name]

const _allowed_PacketCaptureFilter_protocol = ["TCP", "UDP", "Any"]

function check_required(o::PacketCaptureFilter)
    true
end

function validate_property(::Type{ PacketCaptureFilter }, name::Symbol, val)
    if name === Symbol("protocol")
        Swagger.validate_param(name, "PacketCaptureFilter", :enum, val, _allowed_PacketCaptureFilter_protocol)
    end
end


text:
module SystemImageBuilder

using FunctionalData

export buildimage, resetimage

function requires(pkg)
    dir = Pkg.dir(pkg)
    filename = @p joinpath dir "REQUIRE"
    !isfile(filename) && return (pkg,[])
    reqs = @p Base.Pkg.Reqs.read filename | Base.Pkg.Reqs.parse | keys | collect | filter unequal "julia"
    pkg, reqs
end

function recursiverequirements(installed, included)
    req = @p map installed requires | Dict
    @p work included (x->req[x] = [])
    makerec(a,k) = try
        @p map a[k] (x->makerec(a,x)) | flatten | vcat a[k] | filter not*isempty | unique | sort
    catch e
        @show k; rethrow(e)
    end
    map(req, (k,v)->(k,makerec(req,k)))
end

try
    include(joinpath(JULIA_HOME, Base.DATAROOTDIR, "julia", "build_sysimg.jl"))
catch
    script = @p normpath JULIA_HOME ".." ".." "contrib" "build_sysimg.jl"
    include(script)
end

defaultexclude = ["Tk","PyPlot","PyCall","IJulia","SystemImageBuilder","WinRPM","RCall","LMDB","BinDeps",
    "MATLAB",
    "Lint" # makes julia startup really slow
    ]

sysimg = default_sysimg_path
if !isfile(sysimg*".ji")
    sysimg = joinpath(dirname(sysimg), "julia", "sys")
end
if !isfile(sysimg*".ji")
    error("$(sysimg).ji does not seem to be the correct path of sys.ji")
end

resetimage() = buildimage(reset = true)
function buildimage(;exclude = defaultexclude, include = [], targetpath = sysimg, reset = false)
    base_dir = dirname(Base.find_source_file("sysimg.jl"))
    userimg = @p joinpath base_dir "userimg.jl"
    try
        touch(userimg)
    catch
        println()
        println("Path $userimg is not writable. Did you install Julia using a package manager?")
        println("Can't proceed. Please either build Julia from source or download a binary from")
        println("  http://julialang.org/downloads/")
        println()
        return
    end

    if reset
        try rm(userimg) end
        println("##  SystemImageBuilder: building clean sys.ji ...\n")
        build_sysimg(sysimg, "native", force = true)
    else
        installed = sort([k for (k,v) in Pkg.installed()])
        req = recursiverequirements(installed, include)
        packages = filter(x->!in(x,exclude), installed)
        packages = filter(x->isempty(intersect(req[x],exclude)), packages)
        packages = [include; packages]

        println("##  SystemImageBuilder: the following $(length(packages)) packages will be included in the system image:\n")
        println(join(packages, ", "))

        skipped = setdiff(installed, packages)
        if !isempty(skipped)
            println("\n##  SystemImageBuilder: the following $(length(skipped)) packages will be skipped:\n")
            println(join(skipped, ", "))
        end
        println()

        @p map packages (x->"try Base.require(\"$x\") end") | unlines | write userimg

        println("##  SystemImageBuilder: invoking build_sysimg ...\n")
        mkpath(targetpath)
        build_sysimg(targetpath, "native", force = true)
    end
end

end # module


text:
# This file is a part of Julia, but is derived from
# https://github.com/google/double-conversion which has the following license
#
# Copyright 2006-2014, the V8 project authors. All rights reserved.
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are
# met:
#
#     * Redistributions of source code must retain the above copyright
#       notice, this list of conditions and the following disclaimer.
#     * Redistributions in binary form must reproduce the above
#       copyright notice, this list of conditions and the following
#       disclaimer in the documentation and/or other materials provided
#       with the distribution.
#     * Neither the name of Google Inc. nor the names of its
#       contributors may be used to endorse or promote products derived
#       from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
# A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
# OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
# LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
# THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

struct Float
    s::UInt64
    e::Int32
    de::Int32
end

Float() = Float(0,0,0)
Float(x,y) = Float(x,y,Int32(0))
Float(d::AbstractFloat) = Float(_significand(d), _exponent(d))

# Consts
const Float10MSBits = 0xFFC0000000000000 # used normalize(Float)
const FloatSignMask = 0x8000000000000000 # used in normalize(Float)
const FloatSignificandSize = Int32(64)

function normalize(v::Float)
    f = v.s
    e::Int32 = v.e
    while (f & Float10MSBits) == 0
        f <<= 10
        e -= 10
    end
    while (f & FloatSignMask) == 0
        f <<= 1
        e -= 1
    end
    return Float(f,e)
end
function normalize(v::Float64)
    s = _significand(v); e = _exponent(v)
    while (s & HiddenBit(Float64)) == 0
        s <<= UInt64(1)
        e -= Int32(1)
    end
    s <<= UInt64(FloatSignificandSize - SignificandSize(Float64))
    e -=  Int32( FloatSignificandSize - SignificandSize(Float64))
    return Float(s, e)
end

# Float128
#DenormalExponent(::Type{Float128}) = Int32(-ExponentBias(Float128) + 1)
#ExponentMask(::Type{Float128}) = 0x7fff0000000000000000000000000000
#PhysicalSignificandSize(::Type{Float128}) = Int32(112)
#SignificandSize(::Type{Float128}) = Int32(113)
#ExponentBias(::Type{Float128}) = Int32(0x00003fff + PhysicalSignificandSize(Float128))
#SignificandMask(::Type{Float128}) = 0x0000ffffffffffffffffffffffffffff
#HiddenBit(::Type{Float128}) = 0x00010000000000000000000000000000
#uint_t(d::Float128) = reinterpret(UInt128,d)
# Float64
DenormalExponent(::Type{Float64}) = Int32(-ExponentBias(Float64) + 1)
ExponentMask(::Type{Float64}) = 0x7FF0000000000000
PhysicalSignificandSize(::Type{Float64}) = Int32(52)
SignificandSize(::Type{Float64}) = Int32(53)
ExponentBias(::Type{Float64}) = Int32(0x3FF + PhysicalSignificandSize(Float64))
SignificandMask(::Type{Float64}) = 0x000FFFFFFFFFFFFF
HiddenBit(::Type{Float64}) = 0x0010000000000000
uint_t(d::Float64) = reinterpret(UInt64,d)
# Float32
DenormalExponent(::Type{Float32}) = Int32(-ExponentBias(Float32) + 1)
ExponentMask(::Type{Float32}) = 0x7F800000
PhysicalSignificandSize(::Type{Float32}) = Int32(23)
SignificandSize(::Type{Float32}) = Int32(24)
ExponentBias(::Type{Float32}) = Int32(0x7F + PhysicalSignificandSize(Float32))
SignificandMask(::Type{Float32}) = 0x007FFFFF
HiddenBit(::Type{Float32}) = 0x00800000
uint_t(d::Float32) = reinterpret(UInt32,d)
# Float16
DenormalExponent(::Type{Float16}) = Int32(-ExponentBias(Float16) + 1)
ExponentMask(::Type{Float16}) = 0x7c00
PhysicalSignificandSize(::Type{Float16}) = Int32(10)
SignificandSize(::Type{Float16}) = Int32(11)
ExponentBias(::Type{Float16}) = Int32(0x000f + PhysicalSignificandSize(Float16))
SignificandMask(::Type{Float16}) = 0x03ff
HiddenBit(::Type{Float16}) = 0x0400
uint_t(d::Float16) = reinterpret(UInt16,d)

function _exponent(d::T) where T<:AbstractFloat
  isdenormal(d) && return DenormalExponent(T)
  biased_e::Int32 = Int32((uint_t(d) & ExponentMask(T)) >> PhysicalSignificandSize(T))
  return Int32(biased_e - ExponentBias(T))
end
function _significand(d::T) where T<:AbstractFloat
  s = uint_t(d) & SignificandMask(T)
  return !isdenormal(d) ? s + HiddenBit(T) : s
end
isdenormal{T<:AbstractFloat}(d::T) = (uint_t(d) & ExponentMask(T)) == 0

function normalizedbound(f::AbstractFloat)
    v = Float(_significand(f),_exponent(f))
    m_plus = normalize(Float((v.s << 1) + 1, v.e - 1))
    if lowerboundaryiscloser(f)
        m_minus = Float((v.s << 2) - 1, v.e - 2)
    else
        m_minus = Float((v.s << 1) - 1, v.e - 1)
    end
    return Float(m_minus.s << (m_minus.e - m_plus.e), m_plus.e), m_plus
end
function lowerboundaryiscloser(f::T) where T<:AbstractFloat
    physical_significand_is_zero = (uint_t(f) & SignificandMask(T)) == 0
    return physical_significand_is_zero && (_exponent(f) != DenormalExponent(T))
end

(-)(a::Float,b::Float) = Float(a.s - b.s,a.e,a.de)

const FloatM32 = 0xFFFFFFFF

function (*)(this::Float,other::Float)
    a::UInt64 = this.s >> 32
    b::UInt64 = this.s & FloatM32
    c::UInt64 = other.s >> 32
    d::UInt64 = other.s & FloatM32
    ac::UInt64 = a * c
    bc::UInt64 = b * c
    ad::UInt64 = a * d
    bd::UInt64 = b * d
    tmp::UInt64 = (bd >> 32) + (ad & FloatM32) + (bc & FloatM32)
    # By adding 1U << 31 to tmp we round the final result.
    # Halfway cases will be round up.
    tmp += UInt64(1) << 31
    result_f::UInt64 = ac + (ad >> 32) + (bc >> 32) + (tmp >> 32)
    return Float(result_f,this.e + other.e + 64,this.de)
end

const CachedPowers = Float[
  Float(0xfa8fd5a0081c0288, -1220, -348),
  Float(0xbaaee17fa23ebf76, -1193, -340),
  Float(0x8b16fb203055ac76, -1166, -332),
  Float(0xcf42894a5dce35ea, -1140, -324),
  Float(0x9a6bb0aa55653b2d, -1113, -316),
  Float(0xe61acf033d1a45df, -1087, -308),
  Float(0xab70fe17c79ac6ca, -1060, -300),
  Float(0xff77b1fcbebcdc4f, -1034, -292),
  Float(0xbe5691ef416bd60c, -1007, -284),
  Float(0x8dd01fad907ffc3c, -980, -276),
  Float(0xd3515c2831559a83, -954, -268),
  Float(0x9d71ac8fada6c9b5, -927, -260),
  Float(0xea9c227723ee8bcb, -901, -252),
  Float(0xaecc49914078536d, -874, -244),
  Float(0x823c12795db6ce57, -847, -236),
  Float(0xc21094364dfb5637, -821, -228),
  Float(0x9096ea6f3848984f, -794, -220),
  Float(0xd77485cb25823ac7, -768, -212),
  Float(0xa086cfcd97bf97f4, -741, -204),
  Float(0xef340a98172aace5, -715, -196),
  Float(0xb23867fb2a35b28e, -688, -188),
  Float(0x84c8d4dfd2c63f3b, -661, -180),
  Float(0xc5dd44271ad3cdba, -635, -172),
  Float(0x936b9fcebb25c996, -608, -164),
  Float(0xdbac6c247d62a584, -582, -156),
  Float(0xa3ab66580d5fdaf6, -555, -148),
  Float(0xf3e2f893dec3f126, -529, -140),
  Float(0xb5b5ada8aaff80b8, -502, -132),
  Float(0x87625f056c7c4a8b, -475, -124),
  Float(0xc9bcff6034c13053, -449, -116),
  Float(0x964e858c91ba2655, -422, -108),
  Float(0xdff9772470297ebd, -396, -100),
  Float(0xa6dfbd9fb8e5b88f, -369, -92),
  Float(0xf8a95fcf88747d94, -343, -84),
  Float(0xb94470938fa89bcf, -316, -76),
  Float(0x8a08f0f8bf0f156b, -289, -68),
  Float(0xcdb02555653131b6, -263, -60),
  Float(0x993fe2c6d07b7fac, -236, -52),
  Float(0xe45c10c42a2b3b06, -210, -44),
  Float(0xaa242499697392d3, -183, -36),
  Float(0xfd87b5f28300ca0e, -157, -28),
  Float(0xbce5086492111aeb, -130, -20),
  Float(0x8cbccc096f5088cc, -103, -12),
  Float(0xd1b71758e219652c, -77, -4),
  Float(0x9c40000000000000, -50, 4),
  Float(0xe8d4a51000000000, -24, 12),
  Float(0xad78ebc5ac620000, 3, 20),
  Float(0x813f3978f8940984, 30, 28),
  Float(0xc097ce7bc90715b3, 56, 36),
  Float(0x8f7e32ce7bea5c70, 83, 44),
  Float(0xd5d238a4abe98068, 109, 52),
  Float(0x9f4f2726179a2245, 136, 60),
  Float(0xed63a231d4c4fb27, 162, 68),
  Float(0xb0de65388cc8ada8, 189, 76),
  Float(0x83c7088e1aab65db, 216, 84),
  Float(0xc45d1df942711d9a, 242, 92),
  Float(0x924d692ca61be758, 269, 100),
  Float(0xda01ee641a708dea, 295, 108),
  Float(0xa26da3999aef774a, 322, 116),
  Float(0xf209787bb47d6b85, 348, 124),
  Float(0xb454e4a179dd1877, 375, 132),
  Float(0x865b86925b9bc5c2, 402, 140),
  Float(0xc83553c5c8965d3d, 428, 148),
  Float(0x952ab45cfa97a0b3, 455, 156),
  Float(0xde469fbd99a05fe3, 481, 164),
  Float(0xa59bc234db398c25, 508, 172),
  Float(0xf6c69a72a3989f5c, 534, 180),
  Float(0xb7dcbf5354e9bece, 561, 188),
  Float(0x88fcf317f22241e2, 588, 196),
  Float(0xcc20ce9bd35c78a5, 614, 204),
  Float(0x98165af37b2153df, 641, 212),
  Float(0xe2a0b5dc971f303a, 667, 220),
  Float(0xa8d9d1535ce3b396, 694, 228),
  Float(0xfb9b7cd9a4a7443c, 720, 236),
  Float(0xbb764c4ca7a44410, 747, 244),
  Float(0x8bab8eefb6409c1a, 774, 252),
  Float(0xd01fef10a657842c, 800, 260),
  Float(0x9b10a4e5e9913129, 827, 268),
  Float(0xe7109bfba19c0c9d, 853, 276),
  Float(0xac2820d9623bf429, 880, 284),
  Float(0x80444b5e7aa7cf85, 907, 292),
  Float(0xbf21e44003acdd2d, 933, 300),
  Float(0x8e679c2f5e44ff8f, 960, 308),
  Float(0xd433179d9c8cb841, 986, 316),
  Float(0x9e19db92b4e31ba9, 1013, 324),
  Float(0xeb96bf6ebadf77d9, 1039, 332),
  Float(0xaf87023b9bf0ee6b, 1066, 340)]

const CachedPowersLength = length(CachedPowers)
const CachedPowersOffset = 348  # -1 * the first decimal_exponent.
const D_1_LOG2_10 = 0.30102999566398114  #  1 / lg(10)
# Difference between the decimal exponents in the table above.
const DecimalExponentDistance = 8
const MinDecimalExponent = -348
const MaxDecimalExponent = 340

function binexp_cache(min_exponent,max_exponent)
    k = ceil(Integer,(min_exponent+63)*D_1_LOG2_10)
    index = div(CachedPowersOffset+k-1,DecimalExponentDistance) + 1
    cp = CachedPowers[index+1]
    return cp
end

max_stars_repo_name: odow/julia

max_issues_repo_name: odow/julia

max_forks_repo_name: odow/julia


text:
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute(n::Int, limit::Int)::Int
    result = 0
    for i ∈ 1:n
        for j ∈ 1:i ÷ 2
            if binomial(i, j) > limit
                result += i - 2j + 1
                break
            end
        end
    end
    return result
end

compute(100, 1_000_000)

@benchmark compute(100, 1_000_000)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
using D91PVInversion
using IterativeSolvers
using Random
import PyPlot; const plt = PyPlot

# Define functional form of PV anomaly
function q′(x, y, z, A, L, H, p::Params)
    πc = p.π0 - 0.5
    Π = p.Π
    return A*exp(-(x^2 + y^2)/L^2)*exp(-(z-πc)^2/H^2)
end

# Set up inversion domain
params = Params(Float64)
domain = Domain(params, size = (3, 3, 3), x = (-1, 1)./2, y = (-1, 1)./2)
println(params)
println(domain)

# Create solvers
sψ = Solver(params)
sϕ = Solver(params)

# Allocate storage for fields and vectors
ψ, ϕ, q = allocate_fields(domain)
xψ, ∂ψ, bψ, xϕ, ∂ϕ, bϕ = allocate_rhs(domain)

# Set initial values
set_background_ψ!(ψ, domain, params)
set_background_ϕ!(ϕ, domain, params)

# Set PV field
# Set PV field
A = 0
L = 0.1
H = 0.1
set_q!(q, domain, params, q′, A, L, H, params)

# Plot fields
plt.rc("font", size = 8)
fig, axes = plt.subplots(
    figsize = (9.5, 2), nrows = 1, ncols = 5,
    sharey = true, sharex = true, constrained_layout = true
)
for i = 1:5
    im = axes[i].imshow(ψ[:,i-1,:]')
    plt.colorbar(im, ax = axes[i], location = "bottom", aspect = 5, shrink = 0.7)
    axes[i].set_title(string("ψ, j = ", i))
end
axes[1].invert_yaxis()
plt.show()

# Construct linear operators
Lψ = generate_Lψ(ϕ, domain)
Lϕ = generate_Lϕ(ψ, domain)

# Plot operators
plt.figure()
plt.imshow(Matrix(Lψ))
plt.show()
plt.figure()
plt.imshow(Matrix(Lϕ))
plt.show()

# Fill halos
fill_ψ_halos!(ψ, domain, params)
fill_ϕ_halos!(ϕ, domain, params)

# Re-plot fields
plt.rc("font", size = 8)
fig, axes = plt.subplots(
    figsize = (9.5, 2), nrows = 1, ncols = 5,
    sharey = true, sharex = true, constrained_layout = true
)
for i = 1:5
    im = axes[i].imshow(ψ[:,i-1,:]')
    plt.colorbar(im, ax = axes[i], location = "bottom", aspect = 5, shrink = 0.7)
    axes[i].set_title(string("ψ, j = ", i))
end
axes[1].invert_yaxis()
plt.show()

# Re-plot operators
plt.figure()
plt.imshow(Matrix(Lψ))
plt.show()
plt.figure()
plt.imshow(Matrix(Lϕ))
plt.show()

# BEGIN ITERATION HERE
# Set boundary contributions to ψ RHS
set_∂ψ!(∂ψ, ϕ, domain, params)

# Plot boundary contributions to ψ RHS
f∂ψ = field_from_rhs(∂ψ, domain)
fig, axes = plt.subplots(
    figsize = (9.5, 2), nrows = 1, ncols = 5,
    sharey = true, sharex = true, constrained_layout = true
)
for i = 1:5
    im = axes[i].imshow(f∂ψ[:,i-1,:]')
    plt.colorbar(im, ax = axes[i], location = "bottom", aspect = 5, shrink = 0.7)
    axes[i].set_title(string("∂ψ, j = ", i))
end
axes[1].invert_yaxis()
plt.show()

# Set ψ RHS
set_bψ!(bψ, ψ, ϕ, q, domain)

# Plot q and RHS
fbψ = field_from_rhs(bψ, domain)
fig, axes = plt.subplots(
    figsize = (9.5, 2), nrows = 1, ncols = 5,
    sharey = true, sharex = true, constrained_layout = true
)
for i = 1:5
    im = axes[i].imshow(q[:,i-1,:]')
    plt.colorbar(im, ax = axes[i], location = "bottom", aspect = 5, shrink = 0.7)
    axes[i].set_title(string("q, j = ", i))
end
axes[1].invert_yaxis()
plt.show()
fig, axes = plt.subplots(
    figsize = (9.5, 2), nrows = 1, ncols = 5,
    sharey = true, sharex = true, constrained_layout = true
)
for i = 1:5
    im = axes[i].imshow(fbψ[:,i-1,:]')
    plt.colorbar(im, ax = axes[i], location = "bottom", aspect = 5, shrink = 0.7)
    axes[i].set_title(string("bψ, j = ", i))
end
axes[1].invert_yaxis()
plt.show()

# Compute residual
rhs_from_field!(xψ, ψ, domain)
rψ = field_from_rhs(compute_residual(Lψ, xψ, ∂ψ, bψ, domain), domain)

# Plot residual
fig, axes = plt.subplots(
    figsize = (9.5, 2), nrows = 1, ncols = 5,
    sharey = true, sharex = true, constrained_layout = true
)
for i = 1:5
    im = axes[i].imshow(rψ[:,i-1,:]')
    plt.colorbar(im, ax = axes[i], location = "bottom", aspect = 5, shrink = 0.7)
    axes[i].set_title(string("ψ residual, j = ", i))
end
axes[1].invert_yaxis()
plt.show()

# Solve for ψ
Random.seed!(1234)
@. bψ = bψ - ∂ψ
bicgstabl!(xψ, Lψ, bψ; log = false, verbose = true)
relax!(ψ, xψ, sψ, domain; verbose = true)
fill_ψ_halos!(ψ, domain, params)
println("is_converged(sψ) = ", is_converged(sψ))

# Plot solution
fig, axes = plt.subplots(
    figsize = (9.5, 2), nrows = 1, ncols = 5,
    sharey = true, sharex = true, constrained_layout = true
)
for i = 1:5
    im = axes[i].imshow(ψ[:,i-1,:]')
    plt.colorbar(im, ax = axes[i], location = "bottom", aspect = 5, shrink = 0.7)
    axes[i].set_title(string("ψ, j = ", i))
end
axes[1].invert_yaxis()
plt.show()

# Set boundary contributions to ϕ RHS
set_∂ϕ!(∂ϕ, ψ, domain, params)

# Plot boundary contributions to ϕ RHS
f∂ϕ = field_from_rhs(∂ϕ, domain)
fig, axes = plt.subplots(
    figsize = (9.5, 2), nrows = 1, ncols = 5,
    sharey = true, sharex = true, constrained_layout = true
)
for i = 1:5
    im = axes[i].imshow(f∂ϕ[:,i-1,:]')
    plt.colorbar(im, ax = axes[i], location = "bottom", aspect = 5, shrink = 0.7)
    axes[i].set_title(string("∂ϕ, j = ", i))
end
axes[1].invert_yaxis()
plt.show()

# Set ψ RHS
set_bϕ!(bϕ, ψ, ϕ, q, domain)

# Plot q and RHS
fbϕ = field_from_rhs(bϕ, domain)
fig, axes = plt.subplots(
    figsize = (9.5, 2), nrows = 1, ncols = 5,
    sharey = true, sharex = true, constrained_layout = true
)
for i = 1:5
    im = axes[i].imshow(q[:,i-1,:]')
    plt.colorbar(im, ax = axes[i], location = "bottom", aspect = 5, shrink = 0.7)
    axes[i].set_title(string("q, j = ", i))
end
axes[1].invert_yaxis()
plt.show()
fig, axes = plt.subplots(
    figsize = (9.5, 2), nrows = 1, ncols = 5,
    sharey = true, sharex = true, constrained_layout = true
)
for i = 1:5
    im = axes[i].imshow(fbϕ[:,i-1,:]')
    plt.colorbar(im, ax = axes[i], location = "bottom", aspect = 5, shrink = 0.7)
    axes[i].set_title(string("bϕ, j = ", i))
end
axes[1].invert_yaxis()
plt.show()

# Compute residual
rhs_from_field!(xϕ, ϕ, domain)
rϕ = field_from_rhs(compute_residual(Lϕ, xϕ, ∂ϕ, bϕ, domain), domain)

# Plot residual
fig, axes = plt.subplots(
    figsize = (9.5, 2), nrows = 1, ncols = 5,
    sharey = true, sharex = true, constrained_layout = true
)
for i = 1:5
    im = axes[i].imshow(rϕ[:,i-1,:]')
    plt.colorbar(im, ax = axes[i], location = "bottom", aspect = 5, shrink = 0.7)
    axes[i].set_title(string("ϕ residual, j = ", i))
end
axes[1].invert_yaxis()
plt.show()

# Solve for ϕ
@. bϕ = bϕ - ∂ϕ
bicgstabl!(xϕ, Lϕ, bϕ; log = false, verbose = true)
relax!(ϕ, xϕ, sϕ, domain; verbose = true)
fill_ϕ_halos!(ϕ, domain, params)
println("is_converged(sϕ) = ", is_converged(sϕ))

# Plot solution
fig, axes = plt.subplots(
    figsize = (9.5, 2), nrows = 1, ncols = 5,
    sharey = true, sharex = true, constrained_layout = true
)
for i = 1:5
    im = axes[i].imshow(ϕ[:,i-1,:]')
    plt.colorbar(im, ax = axes[i], location = "bottom", aspect = 5, shrink = 0.7)
    axes[i].set_title(string("ϕ, j = ", i))
end
axes[1].invert_yaxis()
plt.show()


text:
"""
Base class for `RTree` node (`Branch` or `Leaf`).
"""
abstract type Node{T,N,V} end

dimtype(::Type{<:Node{T}}) where T = T
dimtype(node::Node) = dimtype(typeof(node))
Base.ndims(::Type{<:Node{<:Any, N}}) where N = N
Base.ndims(node::Node) = ndims(typeof(node))
mbrtype(::Type{<:Node{T,N}}) where {T,N} = Rect{T,N}
mbrtype(node::Node) = mbrtype(typeof(node))
mbr(node::Node) = node.mbr
parent(node::Node) = node.parent
hasparent(node::Node) = node.parent !== nothing
children(node::Node) = node.children
@inline Base.isempty(node::Node) = isempty(children(node))
@inline Base.length(node::Node) = length(children(node))
@inline Base.getindex(node::Node, i::Integer) = getindex(children(node), i)

# `node` position among its siblings in `children` vector of the parent or nothing
pos_in_parent(node::Node) =
    hasparent(node) ? findfirst(c -> c === node, children(parent(node))) : nothing

# update mbr to match the combined MBR of its children
syncmbr!(node::Node) =
    node.mbr = !isempty(node) ? mapreduce(mbr, combine, node.children) : empty(mbrtype(node))

# check whether type E could be used as the data element type in R-tree
# FIXME also need to check whether E could be converted to V of Node
check_eltype_rtree(::Type{E}, ::Type{N}) where {N<:Node, E} =
    check_hasmbr(mbrtype(N), E)

check_eltype_rtree(el::Any, node::Node) =
    check_eltype_rtree(typeof(el), typeof(node))

# equality check for R-tree elements (by MBR and, optionally, ID)
isequal_rtree(el::Any, reg::Region, id::Any = nothing) =
    isequal_rtree(idtrait(typeof(el)), el, reg, id)

isequal_rtree(::Type{<:HasID}, el::Any, reg::Region, key::Any) =
    isequal(id(el), convert(idtype(idtrait(typeof(el))), key)) && (mbr(el) == reg)

# if the tree has IDs but `nothing` supplied, skip ID check
isequal_rtree(::Type{<:HasID}, el::Any, reg::Region, ::Nothing = nothing) =
    (mbr(el) == reg)

isequal_rtree(::Type{HasNoID}, el::Any, reg::Region, ::Nothing = nothing) =
    (mbr(el) == reg)

"""
R-Tree leaf (level 1 node).
Its children are data elements of type `V`.
"""
mutable struct Leaf{T,N,V} <: Node{T,N,V}
    parent::Union{Node{T,N,V}, Nothing} # actually Branch{T,N,V}, but cannot do without forward decl
    mbr::Rect{T,N}
    children::Vector{V}

    function Leaf{T,N,V}(parent::Union{Nothing, Node{T,N,V}},
                         ::Type{V}) where {T,N,V}
        new{T,N,V}(parent, empty(Rect{T,N}), Vector{V}())
    end
end

level(node::Leaf) = 1 # always

"""
R-Tree node for levels above 1 (non-`Leaf`).
"""
mutable struct Branch{T,N,V} <: Node{T,N,V}
    parent::Union{Branch{T,N,V}, Nothing}
    level::Int  # level (thickness) in the R-tree (0=leaves)
    mbr::Rect{T,N}
    children::Union{Vector{Branch{T,N,V}},Vector{Leaf{T,N,V}}}

    function Branch{T,N,V}(parent::Union{Nothing, Branch{T,N,V}}, ::Type{C}) where {T,N,V, C<:Node{T,N,V}}
        new{T,N,V}(parent, 0, empty(Rect{T,N}), Vector{C}())
    end
end

level(node::Branch) = node.level

function Base.setindex!(node::Branch, child::Node, i::Integer)
    (level(node) == level(child) + 1) ||
        throw(SpatialIndexException("Parent ($(level(node)) and child ($(level(child))) levels don't match"))
    setindex!(children(node), child, i)
    child.parent = node
end


text:
# development testing for multithreaded convolution

using IncrementalInference
using Test
using TensorCast

##

@testset "Basic ContinuousScalar example to ensure multithreaded convolutions work..." begin

##

@show Threads.nthreads()

N = 100

# Start with an empty factor graph
fg = initfg()

# add the first node
addVariable!(fg, :x0, ContinuousScalar, N=N)

# this is unary (prior) factor and does not immediately trigger autoinit of :x0.
addFactor!(fg, [:x0], Prior(Normal(0,1)))

addVariable!(fg, :x1, ContinuousScalar, N=N)
# P(Z | :x1 - :x0 ) where Z ~ Normal(10,1)
addFactor!(fg, [:x0, :x1], LinearRelative(Normal(10.0,1)) , threadmodel=MultiThreaded)

pts_ = approxConv(fg, :x0x1f1, :x1, N=N)
@cast pts[i,j] := pts_[j][i]

@test 0.95*N <= sum(abs.(pts .- 10.0) .< 5.0)

##

end

#


text:
struct DayAheadIndices <: AbstractModelIndices
    hours::Vector{Int}
    plants::Vector{Plant}
    segments::Vector{Int}
    bids::Vector{Int}
    blockbids::Vector{Int}
    blocks::Vector{Int}
    hours_per_block::Vector{Vector{Int}}
end
plants(indices::DayAheadIndices) = indices.plants

function HydroModels.modelindices(data::DayAheadData, horizon::Horizon; areas::Vector{Area} = [0], rivers::Vector{River} = [:All])
    hours = collect(1:num_hours(horizon))
    plants = plants_in_areas_and_rivers(hydrodata(data), areas, rivers)
    if isempty(plants)
        error("No plants in given set of price areas and rivers")
    end
    segments = collect(1:2)
    bids = collect(1:length(data.bidlevels[1]))
    blockbids = collect(1:length(data.bidlevels[1])-2)
    hours_per_block = [collect(h:ending) for h in hours for ending in hours[h+data.regulations.blockminlength-1:end]]
    blocks = collect(1:length(hours_per_block))
    return DayAheadIndices(hours, plants, segments, bids, blockbids, blocks, hours_per_block)
end

function blockbidlevel(b, bidlevels, hours_per_block)
    return mean([bidlevels[i][2:end-1] for i in hours_per_block[b]])
end


text:
using Documenter
using FinancialDerivatives

makedocs(
    format = :html,
    sitename = "FinancialDerivatives.jl",
    pages = [
        "index.md",
        "getting_started.md",
        "api.md"
    ]
)

deploydocs(
    repo = "github.com/JuliaQuant/FinancialDerivatives.jl.git"
)


text:
module PosteriorAnalysis

using MacroTools
using Iterators
using AutoHashEquals
using Requires

import Base: getindex, keys, size, map, vcat, show, ndims

export
    # vars
    vartype,
    varsize,
    vectorview,
    # types
    PosteriorDraws,
    drawtype,
    vectorview,
    # operations
    addvars,
    dropvars,
    @pdmap

include("types.jl")
include("operations.jl")
@require Mamba include("mamba.jl")

end # module


text:
function solve(
    i::MinimumBudget{UniformMatroidInstance{Float64, Maximise}, Int},
    ::DynamicProgramming,
)
    # Consider items above δ.
    # Recurrence relation:
    #  V[m, t, δ] = max { V[m, t, δ + 1],  v_δ + V[m - 1, t - w_δ, δ + 1] }
    #                     \_ no item δ _/ \________ take item δ ________/
    V = Array{Float64, 3}(
        undef,
        i.instance.m,
        dimension(i) + 1,
        i.min_budget + 1,
    )
    S = Dict{Tuple{Int, Int, Int}, Vector{Int}}()

    # Initialise: µ == 1, just take the best element among [δ+1, d], as long as
    # the budget constraint is satisfied; δ == d, take nothing, whatever the
    # budget; β == 0, a simple m-set problem.
    for µ in 1:(i.instance.m) # β == 0
        for δ in (dimension(i) - 1):-1:0 # δ < d
            # Due to δ, there may be too few items (with respect to µ).
            # i.e., there are d - δ items to consider.
            m_ = min(dimension(i) - δ, µ)

            items =
                collect(
                    partialsortperm(
                        i.instance.rewards[(δ + 1):dimension(i)],
                        1:m_,
                        rev=true,
                    ),
                ) .+ δ
            V[µ, δ + 1, 0 + 1] = sum(i.instance.rewards[o] for o in items)
            S[µ, δ, 0] = items
        end
    end

    for β in 0:(i.min_budget)
        # δ == d
        for µ in 1:(i.instance.m)
            if β == 0
                # Optimum solution: take no object.
                V[µ, dimension(i) + 1, β + 1] = 0.0
                S[µ, dimension(i), β] = Int[]
            else
                # No solution.
                V[µ, dimension(i) + 1, β + 1] = -Inf
                S[µ, dimension(i), β] = [-1] # Nonsensical solution.
            end
        end

        # δ < d
        for δ in (dimension(i) - 1):-1:0
            # Filter the available objects (i.e. above δ+1 and satisfying the
            # budget constraint) and retrieve just their indices.
            all_objects = collect(enumerate(i.weights))
            all_objects = collect(filter(o -> o[1] >= δ + 1, all_objects))
            all_objects = collect(filter(o -> o[2] >= β, all_objects))
            all_objects = collect(t[1] for t in all_objects)

            # No feasible solution?
            if length(all_objects) == 0
                V[1, δ + 1, β + 1] = -Inf
                S[1, δ, β] = [-1] # Nonsensical solution.
                continue
            end

            # Take the best one available object.
            v, x = findmax(collect(i.instance.rewards[t] for t in all_objects))
            x = all_objects[x]

            V[1, δ + 1, β + 1] = v
            S[1, δ, β] = [x]
        end
    end

    # Dynamic part.
    for β in 1:(i.min_budget)
        for µ in 2:(i.instance.m)
            for δ in (dimension(i) - 1):-1:0
                remaining_budget_with_δ = max(0, β - i.weights[δ + 1])
                take_δ =
                    i.instance.rewards[δ + 1] +
                    V[µ - 1, δ + 1 + 1, remaining_budget_with_δ + 1]
                dont_take_δ = V[µ, δ + 1 + 1, β + 1]

                both_subproblems_infeasible =
                    -Inf == take_δ && -Inf == dont_take_δ

                if both_subproblems_infeasible
                    V[µ, δ + 1, β + 1] = -Inf
                    S[µ, δ, β] = [-1]
                elseif take_δ >= dont_take_δ
                    V[µ, δ + 1, β + 1] = take_δ
                    S[µ, δ, β] =
                        vcat(δ + 1, S[µ - 1, δ + 1, remaining_budget_with_δ])
                else
                    V[µ, δ + 1, β + 1] = dont_take_δ
                    S[µ, δ, β] = S[µ, δ + 1, β]
                end
            end
        end
    end

    S_refined = Dict(β => S[i.instance.m, 0, β] for β in 0:(i.min_budget))
    return MinBudgetedUniformMatroidSolution(
        i,
        S[i.instance.m, 0, i.min_budget],
        S_refined,
    )
end


text:
#=  Here I try to simulate a realistic college-admissions market, where schools
    start with arbitrary cutoffs and adjust them over time in response to demand.
    At the same time, schools also try to improve their programs, so the equilibrium
    cutoff vector varies dynamically and we get an interesting plot. =#

using DeferredAcceptance
using Plots

m = 6

γ = rand(m)
γ ./= sum(γ)
capacities = 0.8 * ones(m) / m
cutoffs_container = zeros(m, 1)
avg_improvement = .95 .+ 0.1 .* rand(m)
sort!(avg_improvement, rev=true)

"""
Example of a custom demand function. Same as that in `../MNLuniform/` but the
qualities get randomly perturbed each iteration.
"""
function MNL_single_score_demand(cutoffs)
    global cutoffs_container, γ
    cutoffs_container = hcat(cutoffs_container, cutoffs)

    sort_order = sortperm(cutoffs)
    cutoffs[sort_order]

    demands = zeros(m)

    prob_of_th = diff([cutoffs[sort_order]; 1])

    for c in 1:m, d in c:m     # For each score threshold
        demands[sort_order[c]] += prob_of_th[d] *
                                  γ[sort_order[c]] / sum(γ[sort_order[1:d]])
    end

    # println("Qualities:  ", γ)
    γ .*= avg_improvement .+ .1 .* randn(m)
    γ ./= sum(γ)

    return demands
end

# Don't worry about the maxit exceeded warning; equilibrium is in the distant future.
nonatomictatonnement(MNL_single_score_demand, capacities, maxit=300, β = 0., verbose=false)

p = plot(cutoffs_container'[2:end, :],
         label = reshape(["p$i, avg. imp. $(round(avg_improvement[i], digits=3))" for i in 1:m], 1, :),
         ls = [:dashdot :dot :solid],
         c = [:crimson :dodgerblue :rebeccapurple :olivedrab :goldenrod :teal],
         title = "Tâtonnement process with $m schools, one score,\nequal capacities, MNL choice, random quality improvements",
         titlefontsize = 12,
         xlabel = "iteration",
         ylabel = "score cutoff",
         ylim = (0, 1),
         legend = :bottom)

# savefig(p, string("examples/tatonnementovertime/plot.pdf"))
# savefig(p, string("examples/tatonnementovertime/plot.png"))

max_stars_repo_name: maxkapur/DeferredAcceptance

max_issues_repo_name: maxkapur/DeferredAcceptance

max_forks_repo_name: maxkapur/DeferredAcceptance


text:
using TimerOutputs

export TimedLinearOperator

mutable struct TimedLinearOperator{T, F, Ft, Fct} <: AbstractLinearOperator{T}
  timer::TimerOutput
  op::AbstractLinearOperator{T}
  prod!::F
  tprod!::Ft
  ctprod!::Fct
end

TimedLinearOperator{T}(timer::TimerOutput, op::AbstractLinearOperator{T},
                       prod!::F, tprod!::Ft, ctprod!::Fct
                       ) where {T,F,Ft,Fct} = TimedLinearOperator{T,F,Ft,Fct}(timer, op, prod!, tprod!, ctprod!)

"""
    TimedLinearOperator(op)
Creates a linear operator instrumented with timers from TimerOutputs.
"""
function TimedLinearOperator(op::AbstractLinearOperator{T}) where {T}
  timer = TimerOutput()
  prod!(res, x, α, β) = @timeit timer "prod" op.prod!(res, x, α, β)
  tprod!(res, x, α, β) = @timeit timer "tprod" op.tprod!(res, x, α, β)
  ctprod!(res, x, α, β) = @timeit timer "ctprod" op.ctprod!(res, x, α, β)
  TimedLinearOperator{T}(timer, op, prod!, tprod!, ctprod!)
end

TimedLinearOperator(op::AdjointLinearOperator) = adjoint(TimedLinearOperator(op.parent))
TimedLinearOperator(op::TransposeLinearOperator) = transpose(TimedLinearOperator(op.parent))
TimedLinearOperator(op::ConjugateLinearOperator) = conj(TimedLinearOperator(op.parent))

for fn ∈ (
  :size,
  :shape,
  :symmetric,
  :issymmetric,
  :hermitian,
  :ishermitian,
  :nprod,
  :ntprod,
  :nctprod,
  :increase_nprod,
  :increase_ntprod,
  :increase_nctprod,
  :reset!,
)
  @eval begin
    $fn(A::TimedLinearOperator) = $fn(A.op)
  end
end

function show(io::IO, op::TimedLinearOperator)
  show(io, op.op)
  show(io, op.timer)
end


text:
@testset "StabChain unit tests" begin
    a,b = [perm"(1,2,4,3)(5)", perm"(1,2,5,4)"]
    sc = StabilizerChain([a,b])
    @test length(sc) == 1
    pt, S, Δ = sc[1]
    @test pt == 1
    @test S isa Vector{<:Perm}
    @test S == [a,b]
    @test Δ isa Schreier
    @test collect(Δ) == [1,2,4,5,3]

    # first point on the orbit:
    u = Δ[1]
    # the first Schreier generator is trivial
    @test isone(u*a*getinv(Δ, 1^a))
    # but the second is not
    @test u*b*getinv(Δ, 1^b) == perm"(2,5)(3,4)"
    c = perm"(2,5)(3,4)"
    h, new_depth = sift(c, sc)
    @test h == c
    @test new_depth == 2

    # so we extend the base
    @test firstmoved(c) == 2
    push!(sc, firstmoved(c))

    # pushing extends base...
    @test length(sc) == 2
    @test length(sc.base) == 2
    @test sc.base[2] == firstmoved(c)
    @test length(sc.sgs) == 2

    # and allocates generating set and the transversal
    @test sc.sgs[2] == Perm{Int}[] # we've just extended basis here!
    @test length(sc.transversals) == 2 # the sc.transversals[2] is right now a garbage

    # add c to generators at depth new_depth=2
    push!(sc, c, new_depth, recompute=true)
    # pushing the generator adds to the basis ath the given depth
    @test sc.sgs[2] == [c]
    # and (re)computes the orbit at the given depth:
    @test sc.transversals[2] == Schreier([c], sc.base[2])

    # c should sift now to the identity at depth 3
    x, d = sift(c, sc)
    @test isone(x)
    @test d == 3

    u = Δ[2]
    @test u == a
    # Schreier generators are trivial
    @test u*a == Δ[2^a]
    @test u*b == Δ[2^b]
    # so we move to the next point on the orbit
    u = Δ[4] # β = 4
    @test u*a == Δ[4^a]
    @test u*b*getinv(Δ, 4^b) == a^2*b
    d = a^2*b
    @test sift(d, sc) == (d, 2)

    push!(sc, d, 2, recompute=true)

    @test collect(sc.transversals[2]) == [2,5,3,4]

    x, d = sift(a^2*b, sc)
    @test isone(x)
    @test d == 3

    @test order(sc) == 20

    m = match(r"size (\d+)", string(sc))
    @test parse(Int, m.captures[1]) == 20

    m = match(r"Orbit:\s+(\[.*\])", sprint(show, MIME"text/plain"(), sc))
    @test Meta.parse(m.captures[1]) |> eval == collect(sc.transversals[1])
end

@testset "PermGroups" begin
    @test PermGroup([perm"(1)"]) isa PermGroup
	@test PermGroup(perm"(1)") isa PermGroup
    G = PermGroup(perm"(1)")
    @test isdefined(G, :stabchain) == false
    @test order(G) == 1
	@test order(G) isa BigInt
	@test order(Int, G) isa Int
    @test isdefined(G, :stabchain) == true

	G = PermGroup(perm"(1,2,3,4)", perm"(1,2)(4)")
    @test StabilizerChain(G) isa StabilizerChain
    sc = StabilizerChain(G)
    @test length(sc) == 3
    @test order(sc) == factorial(4)
	@test order(sc) isa BigInt
	@test order(Int, sc) isa Int

    @test length(base(G)) == 3
    @test order(G) == factorial(4)

    SN(n) = [Perm(circshift(collect(1:n), -1)), Perm([[2,1];3:n])]

    for N in 6:20
        S = SN(N)
        S = [S; [Perm(randperm(N)) for _ in 1:3]]
        G = PermGroup(S)
        @test order(G) == factorial(N)
    end

    G = PermGroup(perm"(1,2,3,4)", perm"(1,2)")

    m = match(r"order (\d+)", string(G))
    @test m === nothing

    @test perm"(1,3)" in G
    @test perm"(1,5)" ∉ G

    m = match(r"order (\d+)", string(G))
    @test parse(Int, m.captures[1]) == 24

    A = PermGroup(perm"(1,2,3)", perm"(2,3,4)")
    @test order(A) == 12
    m = match(r"order (\d+)", string(A))
    @test parse(Int, m.captures[1]) == 12

    @test perm"(1,2)" ∉ A
    @test perm"(1,3,4)" in A
end

@testset "Iterate over PermGroup" begin
	K1 = PermGroup([perm"(5,6)", perm"(1,2,3,4,5,6)"]) # Symmetric group on 6 symbols
	elements = [g for g in K1]
    @test elements isa Vector{Perm{Int64}}
	uniq_elements = unique(elements)
    @test order(K1) == length(uniq_elements) == 720
    @test uniq_elements == elements

	K2 = PermGroup(Perm{Int16}[perm"(3,4,5)", perm"(1,2,3,4,5)"]) # Alternating group on 5 symbols
	elements = [g for g in K2]
    @test elements isa Vector{Perm{Int16}}
	uniq_elements = unique(elements)
    @test order(K2) == length(uniq_elements) == 60
	@test uniq_elements == elements
end


text:
export ShamanskiiDir

function ShamanskiiDir(nlp            :: AbstractNLPModel,
                       nlp_at_x       :: NLPAtX,
                       x              :: AbstractVector,
                       Nwtdirection   :: Function)

    dₙ = Nwtdirection(nlp_at_x.Hx, nlp_at_x.gx, verbose = false)

    xₚ = x + dₙ

    dS = Nwtdirection(nlp_at_x.Hx, grad(nlp, xₚ), verbose = false)

    return xₚ, dS
end


text:
export NoRegularizer, TikonovRegularizer, FrobeniusNormRegularizer, Type1, Type2, RollingMemory, RestartedMemory, AndersonAccelerator
export QRDecomp, NormalEquations

"""
    AbstractRegularizer

Abstract supertype for Anderson Acceleration Type-II regularisation schemes.
"""
abstract type AbstractRegularizer end
struct TikonovRegularizer <: AbstractRegularizer end
struct FrobeniusNormRegularizer <: AbstractRegularizer end
struct NoRegularizer <: AbstractRegularizer end

"""
    AbstractLeastSquaresMethod

Abstract supertype for least squares solution method when Type-II acceleration is used
"""
abstract type AbstractLeastSquaresMethod end
struct NormalEquations <: AbstractLeastSquaresMethod end
struct QRDecomp <: AbstractLeastSquaresMethod end

"""
    AbstractBroydenType

Abstract supertype for the Broyden type of the accelerator.
"""
abstract type AbstractBroydenType end
struct Type1 <: AbstractBroydenType end
struct Type2{LSM} <: AbstractBroydenType
  function Type2{LSM}() where {LSM <: AbstractLeastSquaresMethod}
    return new{LSM}()
  end
end
Type2() = Type2{QRDecomp}()

"""
    AbstractMemory

Abstract supertype for the memory management of the accelerator.
"""
abstract type AbstractMemory end

"""
    RollingMemory

The accelerator will append new iterates to the history and discard the oldest iterate.
"""
struct RollingMemory <: AbstractMemory end

"""
    RestartedMemory

The accelerator will delete the history once the memory buffers are full.
"""
struct RestartedMemory <: AbstractMemory end

"""
    AndersonAccelerator{T, BT, MT, RE} <: AbstractAccelerator

Accelerator object implementing Anderson Acceleration. Parameterized by:

 - T: AbstractFloat, floating-point type
 - BT: Broyden-type, i.e. Type1 or Type2
 - MT: AbstractMemory, how full memory buffers are handled
 - RE: AbstractRegularizer

"""
mutable struct AndersonAccelerator{T, BT, MT, RE}  <: AbstractAccelerator
  init_phase::Bool
  mem::Int
  min_mem::Int
  dim::Int
  iter::Int
  num_accelerated_steps::Int
  x_last::Vector{T}
  g_last::Vector{T}
  f::Vector{T}
  f_last::Vector{T}
  eta::Vector{T}
  F::Matrix{T}
  X::Matrix{T}
  G::Matrix{T}
  M::Matrix{T}
  Q::Matrix{T}
  R::Matrix{T}
  λ::T # regularisation parameter
  success::Bool # a flag to indicate whether the last attempted acceleration was successful
  acceleration_status::Vector{Tuple{Int, Symbol}}
  num_safe_accepted::Int
  num_safe_declined::Int
  activate_logging::Bool

  function AndersonAccelerator{T, BT, MT, RE}() where {T <: AbstractFloat, RE <: AbstractRegularizer, BT <: AbstractBroydenType, MT <: AbstractMemory}
    new(true, 0, 3, 0, 0, 0, zeros(T, 1), zeros(T, 1), zeros(T, 1),  zeros(T, 1), zeros(T, 1), zeros(T, 1), zeros(T, 1, 1), zeros(T, 1, 1), zeros(T, 1, 1), zeros(T, 1, 1), zeros(T, 1, 1), zeros(T, 1, 1), zero(T), false, Vector{Tuple{Int, Symbol}}(undef, 0), 0, 0, false)
  end

  function AndersonAccelerator{T, BT, MT, RE}(dim::Int; mem::Int = 10, min_mem::Int = 3, λ::T = T(1e-8), activate_logging::Bool = false) where {T <: AbstractFloat, RE <: AbstractRegularizer, BT <: AbstractBroydenType, MT <: AbstractMemory}
    mem, min_mem = arguments_check(mem, min_mem, dim, λ)

    F = zeros(T, dim, mem)
    X = zeros(T, dim, mem)
    G = zeros(T, dim, mem)
    M = zeros(T, mem, mem)
    Q = zeros(T, 0, 0)
    R = zeros(T, 0, 0)
    x_last = zeros(T, dim)
    new(true, mem, min_mem, dim, 0, 0, x_last, zeros(T, dim), zeros(T, dim),  zeros(T, dim), zeros(T, mem), F, X, G, M, Q, R, λ, false, Vector{Tuple{Int, Symbol}}(undef, 0), 0, 0, activate_logging)
  end

  # specific constructor for Type 2 with QR decomposition
  function AndersonAccelerator{T, Type2{QRDecomp}, RestartedMemory, RE}(dim::Int; mem::Int = 10, min_mem::Int = 3, λ::T = T(1e-8), activate_logging::Bool = false) where {T <: AbstractFloat, RE <: AbstractRegularizer}

    mem, min_mem = arguments_check(mem, min_mem, dim, λ)

    # for QR decomposition we don't need the caches F, X, M
    x_last = zeros(T, 0)
    F = zeros(T, 0, 0)
    X = zeros(T, 0, 0)
    G = zeros(T, dim, mem)
    M = zeros(T, 0, 0)
    Q = zeros(T, dim, mem)
    R = zeros(T, mem, mem)

    new(true, mem, min_mem, dim, 0, 0, x_last, zeros(T,dim), zeros(T, dim), zeros(T, dim),  zeros(T, dim), F, X, G, M, Q, R, λ, false, Vector{Tuple{Int, Symbol}}(undef, 0), 0, 0, activate_logging)
  end

end

# define some default constructors for parameters
AndersonAccelerator(args...; kwargs...)  = AndersonAccelerator{Float64, Type2{QRDecomp}, RestartedMemory, NoRegularizer}(args...; kwargs...)
AndersonAccelerator{T}(args...; kwargs...) where {T <: AbstractFloat} = AndersonAccelerator{T,  Type2{QRDecomp}, RestartedMemory, NoRegularizer}(args...; kwargs...)
AndersonAccelerator{BT}(args...; kwargs...) where {BT <: AbstractBroydenType} = AndersonAccelerator{Float64, BT, RestartedMemory, NoRegularizer}(args...; kwargs...)
AndersonAccelerator{M}(args...; kwargs...) where {M <: AbstractMemory} = AndersonAccelerator{Float64, Type2{NormalEquations}, M, NoRegularizer}(args...; kwargs...)
AndersonAccelerator{RE}(args...; kwargs...) where {RE <: AbstractRegularizer} = AndersonAccelerator{Float64,  Type2{NormalEquations}, RestartedMemory, RE}(args...; kwargs...)

# restrict certain combination of parameters
AndersonAccelerator{T, Type2{QRDecomp}, RollingMemory, RE}(dim::Int; mem::Int = 10, min_mem::Int = 3, λ::T = T(1e-8), activated::Bool = true) where {T <: AbstractFloat, RE <: AbstractRegularizer}= throw(ArgumentError("Anderson - Type 2 with QR decomposition and RollingMemory not yet implemented."))
AndersonAccelerator{T, Type2{QRDecomp}, RestartedMemory, RE}(dim::Int; mem::Int = 10, min_mem::Int = 3, λ::T = T(1e-8), activated::Bool = true) where {T <: AbstractFloat, RE <: Union{FrobeniusNormRegularizer, TikonovRegularizer}} = throw(ArgumentError("Anderson - Type 2 with QR decomposition and regularisation not yet implemented."))

get_type(::AndersonAccelerator{T, BT, M, RE}) where {T, RE, BT, M} = BT
get_memory(::AndersonAccelerator{T, BT, M, RE}) where {T, RE, BT, M} = M

get_regularizer(::AndersonAccelerator{T, BT, M, RE}) where {T, RE, BT, M} = RE
was_successful(aa::AndersonAccelerator) = aa.success

function arguments_check(mem::Int, min_mem::Int, dim::Int, λ::T) where {T <: AbstractFloat}
  mem <= 1 && throw(ArgumentError("Memory has to be bigger than one."))
  min_mem <= 1 && throw(ArgumentError("Minimum memory has to be bigger than one."))
  dim <= 0 && throw(ArgumentError("Dimension has to be a positive integer."))
  (λ <= 0 || λ >= 1) && throw(ArgumentError("Regularisation parameter λ has to be between 0 and 1."))

  # mem shouldn't be bigger than the dimension
  mem = min(mem, dim)
  min_mem = min(min_mem, mem)

  return mem, min_mem
end

function log!(aa::AndersonAccelerator, iter::Int, status::Symbol)
  aa.activate_logging && push!(aa.acceleration_status, (iter, status))
  if status == :acc_guarded_declined
    aa.num_safe_declined += 1
  elseif status == :acc_guarded_accepted
    aa.num_accelerated_steps += 1
    aa.num_safe_accepted += 1
  elseif status == :acc_unguarded
    aa.num_accelerated_steps += 1
  end
  return nothing
end

"Put accelerator into a fresh state, i.e. forget about past history."
function restart!(aa::AndersonAccelerator{T, BT, MT, RE}) where {T <: AbstractFloat, RE <: AbstractRegularizer, BT <: AbstractBroydenType, MT <: AbstractMemory}
  # For performance reasons we avoid the deletion operation here, we just reset the aa.iter pointer which
  # tracks the number of filled columns
  # aa.F .= 0;
  # aa.X .= 0;
  # aa.G .= 0;
  # aa.Q .= 0;
  # aa.R .= 0;

  # aa.f .= 0; we keep it if a solver needs to query f for safeguarding
  aa.f_last .= 0;
  aa.g_last .= 0;
  aa.x_last .= 0;
  aa.eta .= 0;

  aa.iter = 0
  aa.init_phase = true
  aa.success = false

end

"Reset the pointer that tracks the number of valid cols in the cached history of past vectors."
function empty_caches!(aa::AndersonAccelerator)
  #to save time we can leave the old data in here, this is a potential source of bugs if the indexing of valid data columns becomes invalid
  # aa.F .= 0;
  # aa.X .= 0;
  # aa.G .= 0;
  aa.iter = 0 #this is important as for the RestartedMemory holds information on how many rows are full with valid data
end

"""
  update!(aa, g, x, iter)

- Update history of accelerator `aa` with iterates g = g(xi)
- Computes residuals f = x - g
- The iteration `iter` is passed in for logging
"""
function update!(aa::AndersonAccelerator{T, BT, MT, RE}, g::AbstractVector{T}, x::AbstractVector{T}, iter::Int) where {T <: AbstractFloat, RE <: AbstractRegularizer, BT <: AbstractBroydenType, MT <: AbstractMemory}
    # compute residual
    @. aa.f = x - g

    if aa.init_phase
      set_prev_vectors!(aa.x_last, aa.g_last, aa.f_last, x, g, aa.f)
      aa.init_phase = false
      return nothing
    end

    j = (aa.iter % aa.mem) + 1 # (aa.iter % aa.mem) number of cols filled, j is the next col where data should be entered

    if j == 1 && aa.iter != 0
      # for the RestartedMemory approach we want to flush the data cache matrices and start from scratch
      # for the RollingMemory approach we just start overwriting the oldest data
      apply_memory_approach!(aa, iter)
    end

    # store Δx, Δg, Δf in X, G, F
    fill_caches!(j, aa.X, aa.G, aa.F, x, g, aa.f, aa.x_last, aa.g_last, aa.f_last)

    # set previous values for next iteration
    set_prev_vectors!(aa.x_last, aa.g_last, aa.f_last, x, g, aa.f)

    aa.iter += 1
  return nothing
end

# This method is used for TypeII with QR decomposition
function update!(aa::AndersonAccelerator{T, Type2{QRDecomp}, RestartedMemory, RE}, g::AbstractVector{T}, x::AbstractVector{T}, iter::Int) where {T <: AbstractFloat, RE <: AbstractRegularizer}

    # compute residual
    @. aa.f = x - g

    if aa.init_phase
      set_prev_vectors!(aa.g_last, aa.f_last, g, aa.f)
      aa.init_phase = false
      return nothing
    end

    j = (aa.iter % aa.mem) + 1 # (aa.iter % aa.mem) number of cols filled, j is the next col where data should be entered

    if j == 1 && aa.iter != 0
      # for a RestartedMemory approach we want to flush the data cache matrices and start from scratch
      apply_memory_approach!(aa, iter)
    end

    #G[:, j] = Δg
    fill_delta!(j, aa.G, g, aa.g_last)

    # use f_last memory to store Δf = f - f_last
    compute_Δf!(aa.f_last, aa.f)

    # QR decomposition step
    qr!(aa.Q, aa.R, aa.f_last, j)

    # set previous values for next iteration
    set_prev_vectors!(aa.g_last, aa.f_last, g, aa.f)

    aa.iter += 1
  return nothing
end

# adapted from https://github.com/JuliaNLSolvers/NLsolve.jl
# add new column to QR-Factorisation: F = QR
function qr!(Q::AbstractMatrix{T}, R::AbstractMatrix{T}, Δf::AbstractVector{T}, j::Int) where {T <: AbstractFloat}

  n, m = size(Q)
  n == length(Δf) || throw(DimensionMismatch())
  m == LinearAlgebra.checksquare(R) || throw(DimensionMismatch())
  1 ≤ j ≤ m || throw(ArgumentError())

  @inbounds for k in 1:(j - 1)
    qk = view(Q, :, k)
    ip = dot(qk, Δf)
    R[k, j] = ip

    axpy!(-ip, qk, Δf)
  end
  @inbounds begin
    nrm_f = norm(Δf, 2)
    R[j, j] = nrm_f
    @. Q[:, j] = Δf / nrm_f
  end
  return nothing
end

" Update the history matrices X = [Δxi, Δxi+1, ...], G = [Δgi, Δgi+1, ...] and F = [Δfi, Δfi+1, ...]."
function fill_caches!(j::Int, X::AbstractMatrix{T}, G::AbstractMatrix{T}, F::AbstractMatrix{T}, x::AbstractVector{T}, g::AbstractVector{T}, f::AbstractVector{T}, x_last::AbstractVector{T}, g_last::AbstractVector{T}, f_last::AbstractVector{T}) where {T <: AbstractFloat}
  # fill memory with deltas
  @. X[:, j] = x - x_last # Δx
  @. G[:, j] = g - g_last # Δg
  @. F[:, j] = f - f_last # Δf
  return nothing
end

" Update a single history matrices V = [vgi, vgi+1, ...] ."
function fill_delta!(j::Int, V::AbstractMatrix{T}, v::AbstractVector{T}, v_last::AbstractVector{T}) where {T <: AbstractFloat}
  @inbounds for i in eachindex(v)
    V[i, j] = v[i] - v_last[i]
  end
end

"Compute Δf = f - f_last and store result in f_last."
function compute_Δf!(f_last::AbstractVector{T}, f::AbstractVector{T}) where {T <: AbstractFloat}

  @inbounds for i in eachindex(f_last)
    f_last[i] *= -one(T)
    f_last[i] += f[i]
  end
end

"Store a copy of the last x, g, f to be able to compute Δx, Δg, Δf at the next step."
function set_prev_vectors!(x_last::AbstractVector{T}, g_last::AbstractVector{T}, f_last::AbstractVector{T}, x::AbstractVector{T}, g::AbstractVector{T}, f::AbstractVector{T}) where {T <: AbstractFloat}
  @. x_last = x
  @. g_last = g
  @. f_last = f
  return nothing
end

function set_prev_vectors!(g_last::AbstractVector{T}, f_last::AbstractVector{T}, g::AbstractVector{T}, f::AbstractVector{T}) where {T <: AbstractFloat}
  @. g_last = g
  @. f_last = f
  return nothing
end

"Depending on the AndersonAccelerator parameter AbstractMemory, dispatch on the correct method that handles the case when the memory buffers are full."
function apply_memory_approach!(aa::AndersonAccelerator{T, BT, RestartedMemory, RE}, iter) where {T, RE, BT}
    empty_caches!(aa)
    aa.activate_logging && log!(aa, iter, :memory_full)
    return true
end
apply_memory_approach!(aa::AndersonAccelerator{T, BT, RollingMemory, RE}, iter) where {T, RE, BT} = nothing

"Anderson type-1: eta = M^{-1} X' where M = (X'F)."
function assemble_inv_matrix!(M::AbstractMatrix{T}, X::AbstractMatrix{T}, F::AbstractMatrix{T}, aa::AndersonAccelerator{T, Type1, ME, RE}) where {T <: AbstractFloat, RE <: AbstractRegularizer, ME <: AbstractMemory}
  mul!(M, X', F)
end
"Anderson type-2: eta = M^{-1} X' where M = (F'F)."
function assemble_inv_matrix!(M::AbstractMatrix{T}, X::AbstractMatrix{T}, F::AbstractMatrix{T}, aa::AndersonAccelerator{T, Type2{NormalEquations}, ME, RE}) where {T <: AbstractFloat, RE <: AbstractRegularizer, ME <: AbstractMemory}
  mul!(M, F', F)
end

"Recombine past iterates to compute an accelerated point. Overwrite g with accelerated point."
function accelerate!(g::AbstractVector{T}, x::AbstractVector{T}, aa::AndersonAccelerator{T, BT, ME, RE}, iter::Int) where {T <: AbstractFloat, BT <: AbstractBroydenType, ME <: AbstractMemory, RE <: AbstractRegularizer}
  aa.success = false
  l = min(aa.iter, aa.mem) #number of columns filled with data
  if l < aa.min_mem
     aa.activate_logging && push!(aa.acceleration_status, (iter, :not_enough_cols))
     return nothing
  end

  eta = view(aa.eta, 1:l)
  X = view(aa.X, :, 1:l)
  M = view(aa.M, 1:l, 1:l)
  G = view(aa.G, :, 1:l)
  F = view(aa.F, :, 1:l)

  # compute eta depending on method type
  # Type-1: eta = (X'F)^{-1}X'f
  # Type-2: eta = (F'F)^{-1}F'f
  assemble_inv_matrix!(M, X, F, aa)
  initialise_eta!(eta, aa, X, F)
  regularize!(M, X, F, aa)

  # aa.eta = aa.M   eta
  info = solve_linear_sys!(M, eta, aa)

  if (info < 0 || norm(eta, 2) > 1e4)
    if info < 0
      aa.activate_logging && push!(aa.acceleration_status, (iter, :fail_singular))
    elseif norm(eta, 2) > 1e4
      aa.activate_logging && push!(aa.acceleration_status, (iter, :fail_eta_norm))
    end
    return nothing
  else
    # calculate the accelerated candidate point and overwrite g
    # g = g - G * eta
    # BLAS.gemv!('N', -one(T), G, eta, one(T), g)
    mul!(g, G, eta)
    @. g = aa.g_last - g

    aa.success = true
    return nothing
  end
end

"Recombine past iterates to compute an accelerated point. Overwrite g with accelerated point. Uses QR-decomposition."
function accelerate!(g::AbstractVector{T}, x::AbstractVector{T}, aa::AndersonAccelerator{T, Type2{QRDecomp}, RestartedMemory, RE}, iter::Int) where {T <: AbstractFloat, RE <: AbstractRegularizer}
  aa.success = false
  l = min(aa.iter, aa.mem) #number of columns filled with data
  if l < aa.min_mem
     aa.activate_logging && push!(aa.acceleration_status, (iter, :not_enough_cols))
     return nothing
  end

  eta = view(aa.eta, 1:l)
  G = view(aa.G, :, 1:l)
  Q = view(aa.Q, :, 1:l)
  R = UpperTriangular(view(aa.R, 1:l, 1:l))
  # solve least squares problem ||f_k - η Fk ||_2 where Fk = QR
  # initialise_eta!(eta, aa, X, F)
  mul!(eta, Q', aa.f)

  info = solve_linear_sys!(R, eta, aa)
  if info < 0
    aa.activate_logging && push!(aa.acceleration_status, (iter, :fail_singular))
    return nothing
  end

  # TODO: maybe replace this with a check of the condition number of R
  nrm_eta = norm(eta, 2)
  if isnan(nrm_eta) || nrm_eta > 1e4
    aa.activate_logging && push!(aa.acceleration_status, (iter, :fail_eta_norm))
    return nothing
  else
    # calculate the accelerated candidate point
    #  g .= g - G * eta
    mul!(g, G, eta)
    @. g = aa.g_last - g

    aa.success = true
    return nothing
  end
end

"Anderson Type 1: Initialise eta = X'f."
function initialise_eta!(eta::AbstractVector{T}, aa::AndersonAccelerator{T, Type1}, X::AbstractMatrix{T}, F::AbstractMatrix{T}) where {T <: AbstractFloat}
  mul!(eta, X', aa.f)
end

"Anderson Type 2: Initialise eta = F'f."
function initialise_eta!(eta::AbstractVector{T}, aa::AndersonAccelerator{T, Type2{NormalEquations}}, X::AbstractMatrix{T}, F::AbstractMatrix{T}) where {T <: AbstractFloat}
  mul!(eta, F', aa.f)
end

"Use Tikonov - Regularizer for the Least squares matrix M."
function regularize!(M::AbstractMatrix{T}, X::AbstractMatrix{T}, F::AbstractMatrix{T}, aa::AndersonAccelerator{T, BT, ME, TikonovRegularizer}) where {T <: AbstractFloat, BT <: AbstractBroydenType, ME <: AbstractMemory}
  # add regularisation term
  for i = 1:size(M, 1)
    M[i, i] += aa.λ
  end
  return nothing
end

"Choose regularisation parameter based on the frobenius norms of the matrices X, F (Fu, Zhang, Boyd, 2019)."
function regularize!(M::AbstractMatrix{T}, X::AbstractMatrix{T}, F::AbstractMatrix{T}, aa::AndersonAccelerator{T, BT, ME, FrobeniusNormRegularizer}) where {T <: AbstractFloat, BT <: AbstractBroydenType, ME <: AbstractMemory}
  β = aa.λ * (norm(X)^2 + norm(F)^2)
  # add regularisation term
  for i = 1:size(M, 1)
    M[i, i] += β
  end
  return nothing
end

regularize!(M::AbstractMatrix{T}, X::AbstractMatrix{T}, F::AbstractMatrix{T}, aa::AndersonAccelerator{T, BT, ME, NoRegularizer}) where {T <: AbstractFloat, BT <: AbstractBroydenType, ME <: AbstractMemory} = nothing

function solve_linear_sys!(M::AbstractMatrix{T}, eta::AbstractVector{T}, aa::AndersonAccelerator{T, BT, ME, RE}) where {T <: AbstractFloat, BT <: AbstractBroydenType, ME <: AbstractMemory, RE <: AbstractRegularizer}
  # BLAS gesv! wrapper with error handling
  # solve A X = B and for X and save result in B
  try
    LinearAlgebra.LAPACK.gesv!(M, eta)
    return 1
  catch
      return -1
  end
end

function solve_linear_sys!(R::AbstractMatrix{T}, eta::AbstractVector{T}, aa::AndersonAccelerator{T, Type2{QRDecomp}, ME, NoRegularizer}) where {T <: AbstractFloat, ME <: AbstractMemory}
  try
     LinearAlgebra.ldiv!(R, eta)
      return 1
  catch
      return -1
    end
end

max_issues_repo_name: oxfordcontrol/COSMOAccelerators

max_forks_repo_name: oxfordcontrol/COSMOAccelerators


text:
# MIDI constants
const MIN_MIDI_PITCH = 0
const MAX_MIDI_PITCH = 127
const MIN_MIDI_VELOCITY = 1
const MAX_MIDI_VELOCITY = 127

# Default quantization constants
const QUANTIZE_CUTOFF = 0.5
const DEFAULT_STEPS_PER_SECOND = 100

# MIDI rhythm related constants
const DEFAULT_QPM = 120.0 # Quarter notes per minute
const DEFAULT_TPQ = 220   # Ticks per quarter note
const DEFAULT_STEPS_PER_BAR = 16
const DEFAULT_STEPS_PER_QUARTER = 4


text:
export error_if_not

"""
    error_if_not(checkFn, errorFn)

Creates an `error_if_not` operator, which performs a check for every emission on the source Observable with `checkFn`.
If `checkFn` returns `false`, the operator sends an `error` event and unsubscribes from the observable.

Note: `error_if_not` is an alias for `error_if` operator with inverted `checkFn`.

# Arguments
- `checkFn`: check function with `(data) -> Bool` signature
- `errorFn`: error object generating function with `(data) -> Any` signature, optional

# Producing

Stream of type `<: Subscribable{L}` where `L` refers to type of source stream

# Examples
```jldoctest
using Rocket

source = from([1, 2, 3]) |> error_if_not((data) -> data < 2, (data) -> "CustomError")

subscription = subscribe!(source, lambda(
    on_next  = (d) -> println("Next: ", d),
    on_error = (e) -> println("Error: ", e),
    on_complete = () -> println("Completed")
))
;

# output
Next: 1
Error: CustomError
```

See also: [`error_if`](@ref), [`error_if_empty`](@ref), [`default_if_empty`](@ref), [`logger`](@ref)
"""
error_if_not(checkFn::F, errorFn::E = nothing) where { F, E } = error_if((d) -> !checkFn(d), errorFn)