text:
@testset "Production Problem" begin
    optimizer = POI.ParametricOptimizer(GLPK.Optimizer())

    c = [4.0, 3.0]
    A1 = [2.0, 1.0, 1.0]
    A2 = [1.0, 2.0, 1.0]
    b1 = 4.0
    b2 = 4.0

    x = MOI.add_variables(optimizer, length(c))

    @test typeof(x[1]) == MOI.VariableIndex

    w, cw = MOI.add_constrained_variable(optimizer, POI.Parameter(0))
    y, cy = MOI.add_constrained_variable(optimizer, POI.Parameter(0))
    z, cz = MOI.add_constrained_variable(optimizer, POI.Parameter(0))

    @test MOI.get(optimizer, MOI.VariablePrimal(), w) == 0

    for x_i in x
        MOI.add_constraint(optimizer, MOI.SingleVariable(x_i), MOI.GreaterThan(0.0))
    end

    cons1 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(A1, [x[1], x[2], y]), 0.0)
    MOI.add_constraint(optimizer, cons1, MOI.LessThan(b1))

    cons2 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(A2, [x[1], x[2], z]), 0.0)
    MOI.add_constraint(optimizer, cons2, MOI.LessThan(b2))

    @test cons1.terms[1].coefficient == 2
    @test POI.is_parameter_in_model(optimizer, cons2.terms[3].variable_index)

    obj_func = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([c[1], c[2], 3.0], [x[1], x[2], w]), 0.0)
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), obj_func)
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    MOI.optimize!(optimizer)

    MOI.get(optimizer, MOI.TerminationStatus())

    MOI.get(optimizer, MOI.PrimalStatus())

    @test isapprox.(MOI.get(optimizer, MOI.ObjectiveValue()), 28/3, atol = ATOL)
    @test isapprox.(MOI.get(optimizer, MOI.VariablePrimal(), x[1]), 4/3, atol = ATOL)
    @test isapprox.(MOI.get(optimizer, MOI.VariablePrimal(), x[2]), 4/3, atol = ATOL)

    @test isapprox(MOI.get(optimizer, MOI.ConstraintDual(), cy), -5/3, atol = ATOL)
    @test isapprox(MOI.get(optimizer, MOI.ConstraintDual(), cz), -2/3, atol = ATOL)
    @test isapprox(MOI.get(optimizer, MOI.ConstraintDual(), cw), 3.0, atol = ATOL)

    MOI.get(optimizer, MOI.VariablePrimal(), w)
    MOI.get(optimizer, MOI.VariablePrimal(), y)
    MOI.get(optimizer, MOI.VariablePrimal(), z)

    MOI.set(optimizer, MOI.ConstraintSet(), cw, POI.Parameter(2.0))
    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(1.0))
    MOI.set(optimizer, MOI.ConstraintSet(), cz, POI.Parameter(1.0))

    MOI.optimize!(optimizer)

    @test MOI.get(optimizer, MOI.VariablePrimal(), w) == 2.0
    @test MOI.get(optimizer, MOI.VariablePrimal(), y) == 1.0
    @test MOI.get(optimizer, MOI.VariablePrimal(), z) == 1.0

    @test isapprox.(MOI.get(optimizer, MOI.ObjectiveValue()), 13.0, atol = ATOL)
    @test MOI.get.(optimizer, MOI.VariablePrimal(), x) == [1.0, 1.0]

    @test isapprox(MOI.get(optimizer, MOI.ConstraintDual(), cy), -5/3, atol = ATOL)
    @test isapprox(MOI.get(optimizer, MOI.ConstraintDual(), cz), -2/3, atol = ATOL)
    @test isapprox(MOI.get(optimizer, MOI.ConstraintDual(), cw), 3.0, atol = ATOL)

    MOI.set(optimizer, MOI.ConstraintSet(), cw, POI.Parameter(0.0))

    MOI.optimize!(optimizer)

    @test MOI.get(optimizer, MOI.VariablePrimal(), w) == 0.0
    @test MOI.get(optimizer, MOI.VariablePrimal(), y) == 1.0
    @test MOI.get(optimizer, MOI.VariablePrimal(), z) == 1.0
    @test isapprox.(MOI.get(optimizer, MOI.ObjectiveValue()), 7, atol = ATOL)
    @test MOI.get.(optimizer, MOI.VariablePrimal(), x) == [1.0, 1.0]

    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(-5.0))
    MOI.optimize!(optimizer)
    @test isapprox.(MOI.get(optimizer, MOI.ObjectiveValue()), 12.0, atol = ATOL)
    @test MOI.get.(optimizer, MOI.VariablePrimal(), x) == [3.0, 0.0]

    @test isapprox(MOI.get(optimizer, MOI.ConstraintDual(), cy), 0, atol = ATOL)
    @test isapprox(MOI.get(optimizer, MOI.ConstraintDual(), cz), -4, atol = ATOL)

end

@testset "Production Problem variation on parameters for duals" begin
    optimizer = POI.ParametricOptimizer(GLPK.Optimizer())

    c = [4.0, 3.0]
    A1 = [2.0, 1.0, 3.0]
    A2 = [1.0, 2.0, 0.5]
    b1 = 4.0
    b2 = 4.0

    x = MOI.add_variables(optimizer, length(c))

    @test typeof(x[1]) == MOI.VariableIndex

    w, cw = MOI.add_constrained_variable(optimizer, POI.Parameter(0))
    y, cy = MOI.add_constrained_variable(optimizer, POI.Parameter(0))
    z, cz = MOI.add_constrained_variable(optimizer, POI.Parameter(0))

    @test MOI.get(optimizer, MOI.VariablePrimal(), w) == 0

    for x_i in x
        MOI.add_constraint(optimizer, MOI.SingleVariable(x_i), MOI.GreaterThan(0.0))
    end

    cons1 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(A1, [x[1], x[2], y]), 0.0)
    ci1 = MOI.add_constraint(optimizer, cons1, MOI.LessThan(b1))

    cons2 = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(A2, [x[1], x[2], z]), 0.0)
    ci2 = MOI.add_constraint(optimizer, cons2, MOI.LessThan(b2))

    @test cons1.terms[1].coefficient == 2
    @test POI.is_parameter_in_model(optimizer, cons2.terms[3].variable_index)

    obj_func = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([c[1], c[2], 2.0], [x[1], x[2], w]), 0.0)
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), obj_func)
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    MOI.optimize!(optimizer)

    MOI.get(optimizer, MOI.TerminationStatus())

    MOI.get(optimizer, MOI.PrimalStatus())

    @test isapprox(MOI.get(optimizer, MOI.ObjectiveValue()), 28/3, atol = ATOL)
    @test isapprox(MOI.get(optimizer, MOI.VariablePrimal(), x[1]), 4/3, atol = ATOL)
    @test isapprox(MOI.get(optimizer, MOI.VariablePrimal(), x[2]), 4/3, atol = ATOL)

    @test isapprox(MOI.get(optimizer, MOI.ConstraintDual(), cy), -5, atol = ATOL)
    @test isapprox(MOI.get(optimizer, MOI.ConstraintDual(), cz), -2/6, atol = ATOL)
    @test isapprox(MOI.get(optimizer, MOI.ConstraintDual(), cw), 2.0, atol = ATOL)

    @test isapprox(MOI.get(optimizer, MOI.ConstraintDual(), cy), 3*MOI.get(optimizer, MOI.ConstraintDual(), ci1), atol = 1e-4)
    @test isapprox(MOI.get(optimizer, MOI.ConstraintDual(), cz), 0.5*MOI.get(optimizer, MOI.ConstraintDual(), ci2), atol = 1e-4)

    MOI.set(optimizer, MOI.ConstraintSet(), cw, POI.Parameter(2.0))
    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(1.0))
    MOI.set(optimizer, MOI.ConstraintSet(), cz, POI.Parameter(1.0))

    MOI.optimize!(optimizer)

    @test isapprox.(MOI.get(optimizer, MOI.ObjectiveValue()), 7.0, atol = ATOL)
    @test MOI.get.(optimizer, MOI.VariablePrimal(), x) == [0.0, 1.0]

    @test isapprox(MOI.get(optimizer, MOI.ConstraintDual(), cy), -9.0, atol = ATOL)
    @test isapprox(MOI.get(optimizer, MOI.ConstraintDual(), cz), -0.0, atol = ATOL)
    @test isapprox(MOI.get(optimizer, MOI.ConstraintDual(), cw), 2.0, atol = ATOL)

    MOI.set(optimizer, MOI.ConstraintSet(), cw, POI.Parameter(0.0))

    MOI.optimize!(optimizer)

    @test isapprox.(MOI.get(optimizer, MOI.ObjectiveValue()), 3.0, atol = ATOL)
    @test MOI.get.(optimizer, MOI.VariablePrimal(), x) == [0.0, 1.0]

    @test isapprox(MOI.get(optimizer, MOI.ConstraintDual(), cy), -9.0, atol = ATOL)
    @test isapprox(MOI.get(optimizer, MOI.ConstraintDual(), cz), -0.0, atol = ATOL)
    @test isapprox(MOI.get(optimizer, MOI.ConstraintDual(), cw), 2.0, atol = ATOL)

    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(-5.0))
    MOI.optimize!(optimizer)
    @test isapprox.(MOI.get(optimizer, MOI.ObjectiveValue()), 14.0, atol = ATOL)
    @test isapprox(MOI.get.(optimizer, MOI.VariablePrimal(), x),[3.5, 0.0], atol = ATOL)

    @test isapprox(MOI.get(optimizer, MOI.ConstraintDual(), cy), 0, atol = ATOL)
    @test isapprox(MOI.get(optimizer, MOI.ConstraintDual(), cz), -2, atol = ATOL)

end


text:
using Optim, HTTP, GLM, LinearAlgebra, Random, Statistics, DataFrames, CSV, FreqTables
function q4()
    # Question 1
    url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2021/master/ProblemSets/PS3-gev/nlsw88w.csv"
    df = CSV.read(HTTP.get(url).body, DataFrame)
    X = [df.age df.white df.collgrad]
    Z = hcat(df.elnwage1, df.elnwage2, df.elnwage3, df.elnwage4,
    df.elnwage5, df.elnwage6, df.elnwage7, df.elnwage8)
    y = df.occupation

    function mlogit(Beta, X, Z,y)

        K = size(X,2)+1
        J = length(unique(y))
        N = length(y)
        bigY = zeros(N,J)
        for j=1:J
            bigY[:,j] = y.==j
        end
        bigBeta = [reshape(Beta,K,J-1) zeros(K)]

        bigZ = zeros(N,J)
        for j=1:J
            bigZ[:,j] = Z[:,j]-Z[:,J]
        end

        num = zeros(N,J)
        dem = zeros(N)
        for j=1:J
            XZ=cat(X,bigZ[:,j],dims=2)
            num[:,j] = exp.(XZ*bigBeta[:,j])
            dem .+= num[:,j]
        end

        P = num./repeat(dem,1,J)

        loglike = -sum( bigY.*log.(P) )

        return loglike
    end

    mlogit_hat_optim = optimize(b-> mlogit(b,X,Z, y), rand(7*(size(X,2)+1)), LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=50))
    println(mlogit_hat_optim.minimizer)

    # Question 2
    println

    # Question 3
    function nlogit(βwc,βbc,λwc,λbc, X, Z,y)

        K = size(X,2)+1
        J = length(unique(y))
        N = length(y)
        bigY = zeros(N,J)

        bigYwc = zeros(N,3)
        for j=1:3
            bigYwc[:,j] = y.==j
        end
        bigβwc = [reshape(βwc,K,3) zeros(K)]

        bigZwc = zeros(N,3)
        for j=1:3
            bigZwc[:,j] = Z[:,j]-Z[:,J]
        end

        numwc1 = zeros(N,3)
        demwc1 = zeros(N)
        numwc = zeros(N,3)
        for j=1:3
            XZwc=cat(X,bigZwc[:,j],dims=2)
            numwc1[:,j] = exp.((XZwc*bigβwc[:,j])./λwc)
            demwc1 .+= numwc1[:,j]
            demwc=(repeat(demwc1,1,3))^λwc
            numwc2=(repeat(demwc1,1,3))^(λwc-1)
            numwc[:,j]=numwc1[:,j].*numwc2
        end

        bigYbc = zeros(N,4)
        for j=1:4
            bigYbc[:,j] = y.==j+3
        end
        bigβbc = [reshape(βbc,K,4) zeros(K)]

        bigZbc = zeros(N,4)
        for j=1:4
            bigZbc[:,j] = Z[:,j+3]-Z[:,J]
        end

        numbc1 = zeros(N,4)
        dembc1 = zeros(N)
        numbc = zeros(N,4)
        for j=1:4
            XZbc=cat(X,bigZbc[:,j],dims=2)
            numbc1[:,j] = exp.((XZbc*bigβbc[:,j])./λbc)
            dembc1 .+= numbc1[:,j]
            dembc=(repeat(dembc1,1,4))^λbc
            numbc2=(repeat(dembc1,1,4))^(λbc-1)
            numbc[:,j]=numbc1[:,j].*numbc2
        end

        num=cat(numwc,numbc,dims=2)
        dem=1+demwc+dembc
        P = num./dem

        loglike = -sum( bigY.*log.(P) )

        return loglike
    end

    nlogit_hat_optim = optimize((βwc,βbc,λwc,λbc)->nlogit(βwc,βbc,λwc,λbc, X, Z,y), rand(7*(size(X,2)+1)+2), LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=50))
    println(nlogit_hat_optim.minimizer)

    end

    # Question 4
    q4()

ext: jl


text:
function chambolle_pock(
    affine_sets::AffineSets,
    conic_sets::ConicSets,
    opt
)::Result

    # Initialize parameters
    p = Params()
    p.theta = opt.initial_theta
    p.adapt_level = opt.initial_adapt_level
    p.window = opt.convergence_window
    p.beta = opt.initial_beta
    p.time0 = time()
    p.norm_b = LinearAlgebra.norm(affine_sets.b, 2)
    p.norm_h = LinearAlgebra.norm(affine_sets.h, 2)
    p.norm_c = LinearAlgebra.norm(affine_sets.c, 2)
    p.rank_update, p.stop_reason, p.update_cont = 0, 0, 0
    p.stop_reason_string = "Not optimized"
    p.target_rank = 2 * ones(length(conic_sets.sdpcone))
    p.current_rank = 2 * ones(length(conic_sets.sdpcone))
    p.min_eig = zeros(length(conic_sets.sdpcone))
    p.dual_feasibility = -1.0
    p.dual_feasibility_check = false
    p.certificate_search = false
    p.certificate_search_min_iter = 0
    p.certificate_found = false
    sol = Array{Result}(undef, 0)
    arc_list = [
        EigSolverAlloc(Float64, sdp.sq_side, opt)
            for (idx, sdp) in enumerate(conic_sets.sdpcone)]
    ada_count = 0

    if opt.max_iter <= 0
        if length(conic_sets.socone) > 0 || length(conic_sets.sdpcone) > 0
            opt.max_iter_local = opt.max_iter_conic
        else
            opt.max_iter_local = opt.max_iter_lp
        end
    else
        opt.max_iter_local = opt.max_iter
    end

    # Print header
    if opt.log_verbose
        print_header_1()
        print_parameters(opt, conic_sets)
        print_constraints(affine_sets)
        if length(conic_sets.socone) + length(conic_sets.sdpcone) > 0
            print_prob_data(conic_sets)
        end
        print_header_2(opt)
    end

    @timeit "Init" begin

        # Scale objective function
        @timeit "normscale alloc" begin
            c_orig, var_ordering = preprocess!(affine_sets, conic_sets)
            A_orig, b_orig = copy(affine_sets.A), copy(affine_sets.b)
            G_orig, h_orig = copy(affine_sets.G), copy(affine_sets.h)
            rhs_orig = vcat(b_orig, h_orig)
        end

        # Diagonal preconditioning
        @timeit "equilibrate" begin
            if opt.equilibration
                M = vcat(affine_sets.A, affine_sets.G)
                UB = maximum(M)
                LB = minimum(M)
                if LB/UB <= opt.equilibration_limit
                    opt.equilibration = false
                end
            end
            if opt.equilibration_force
                opt.equilibration = true
            end
            if opt.equilibration
                @timeit "equilib inner" E, D = equilibrate!(M, affine_sets, opt)
                @timeit "equilib scaling" begin
                    M = E * M * D
                    affine_sets.A = M[1:affine_sets.p, :]
                    affine_sets.G = M[affine_sets.p + 1:end, :]
                    rhs = E * rhs_orig
                    affine_sets.b = rhs[1:affine_sets.p]
                    affine_sets.h = rhs[affine_sets.p + 1:end]
                    affine_sets.c = D * affine_sets.c
                end
            else
                E = LinearAlgebra.Diagonal(zeros(affine_sets.m + affine_sets.p))
                D = LinearAlgebra.Diagonal(zeros(affine_sets.n))
            end
        end

        # Scale the off-diagonal entries associated with p.s.d. matrices by √2
        @timeit "normscale" norm_scaling(affine_sets, conic_sets)

        # Initialization
        pair = PrimalDual(affine_sets)
        a = AuxiliaryData(affine_sets, conic_sets)
        map_socs!(pair.x, conic_sets, a)
        residuals = Residuals(p.window)

        # Diagonal scaling
        M = vcat(affine_sets.A, affine_sets.G)
        Mt = M'

        # Stepsize parameters and linesearch parameters
        if !opt.approx_norm
            @timeit "svd" if minimum(size(M)) >= 2
                try
                    spectral_norm = Arpack.svds(M, nsv=1)[1].S[1]
                catch
                    println("    WARNING: Failed to compute spectral norm of M, shifting to Frobenius norm")
                    spectral_norm = LinearAlgebra.norm(M)
                end
            else
                F = LinearAlgebra.svd!(Matrix(M))
                spectral_norm = maximum(F.S)
            end
        else
            spectral_norm = LinearAlgebra.norm(M)
        end
        if spectral_norm < 1e-10
            spectral_norm = 1.0
        end

        # Build struct for storing matrices
        mat = Matrices(M, Mt, affine_sets.c)

        # Initial primal and dual steps
        p.primal_step = 1. / spectral_norm
        p.primal_step_old = p.primal_step
        p.dual_step = p.primal_step

    end

    # Initialization
    if opt.advanced_initialization
        pair.x .= p.primal_step .* mat.c
        LinearAlgebra.mul!(a.Mx, mat.M, pair.x)
        LinearAlgebra.mul!(a.Mx_old, mat.M, pair.x_old)
    end

    # Fixed-point loop
    @timeit "CP loop" for k in 1:2*opt.max_iter_local

        # Update iterator
        p.iter = k

        # Primal step
        @timeit "primal" primal_step!(pair, a, conic_sets, mat, opt, p, arc_list, p.iter)

        # Linesearch (dual step)
        if opt.line_search_flag
            @timeit "linesearch" linesearch!(pair, a, affine_sets, mat, opt, p)
        else
            @timeit "dual step" dual_step!(pair, a, affine_sets, mat, opt, p)
        end

        # Compute residuals and update old iterates
        @timeit "residual" compute_residual!(residuals, pair, a, p, affine_sets)

        # Compute optimality gap and feasibility error
        @timeit "gap" compute_gap!(residuals, pair, a, affine_sets, p)

        if (opt.check_dual_feas && mod(k, opt.check_dual_feas_freq) == 0) ||
                (opt.log_verbose && mod(k, opt.log_freq) == 0 && opt.extended_log2)
            cc = ifelse(p.stop_reason == 6, 0.0, 1.0)*c_orig
            p.dual_feasibility = dual_feas(pair.y, conic_sets, affine_sets, cc, A_orig, G_orig, a)
            p.dual_feasibility_check = true
        else
            p.dual_feasibility_check = false
        end

        # Print progress
        if opt.log_verbose && mod(k, opt.log_freq) == 0
            print_progress(residuals, p, opt, p.dual_feasibility)
        end

        if p.iter < p.certificate_search_min_iter
            continue
        end

        if opt.certificate_search && p.certificate_search

            if p.stop_reason == 6 # Infeasible

                if residuals.dual_obj[k] > +opt.certificate_obj_tol

                    p.dual_feasibility = dual_feas(pair.y, conic_sets, affine_sets, 0*c_orig, A_orig, G_orig, a)
                    p.dual_feasibility_check = true

                    if p.dual_feasibility < opt.tol_feasibility_dual

                        if opt.log_verbose
                            println("---------------------------------------------------------------------------------------")
                            println("    Dual ray found")
                            println("---------------------------------------------------------------------------------------")
                        end

                        p.certificate_found = true
                        p.stop_reason_string *= " [Dual ray found]"

                        break
                    end
                end

            else p.stop_reason == 5 # Unbounded

                if residuals.prim_obj[k] < -opt.certificate_obj_tol

                    if residuals.feasibility[p.iter] < opt.tol_feasibility

                        if opt.log_verbose
                            println("---------------------------------------------------------------------------------------")
                            println("    Primal ray found")
                            println("---------------------------------------------------------------------------------------")
                        end

                        p.certificate_found = true
                        p.stop_reason_string *= " [Primal ray found]"

                        break
                    end
                end
            end

            if (
                residuals.prim_obj[k] < -opt.certificate_fail_tol &&
                residuals.dual_obj[k] < -opt.certificate_fail_tol &&
                residuals.feasibility[p.iter] < -opt.certificate_fail_tol
                ) || isnan(residuals.comb_residual[k])

                if opt.log_verbose
                    println("---------------------------------------------------------------------------------------")
                    println("    Failed to finds certificate")
                    println("---------------------------------------------------------------------------------------")
                end

                p.stop_reason_string *= " [Failed to find certificate]"

                break
            end
        end

        # Check convergence
        p.rank_update += 1
        if residuals.dual_gap[p.iter] <= opt.tol_gap && residuals.feasibility[p.iter] <= opt.tol_feasibility &&
                (!opt.check_dual_feas || p.dual_feasibility < opt.tol_feasibility_dual)
            if convergedrank(a, p, conic_sets, opt) && soc_convergence(a, conic_sets, pair, opt, p) && p.iter > opt.min_iter
                if !p.certificate_search
                    p.stop_reason = 1 # Optimal
                    p.stop_reason_string = "Optimal solution found"

                else
                    if opt.log_verbose
                        println("---------------------------------------------------------------------------------------")
                        println("    Failed to find certificate - type 2")
                        println("---------------------------------------------------------------------------------------")
                    end

                    p.stop_reason_string *= " [Failed to find certificate - type 2]"

                    break
                end
                break
            elseif p.rank_update > p.window
                p.update_cont += 1
                if p.update_cont > 0
                    for (idx, sdp) in enumerate(conic_sets.sdpcone)
                        if p.current_rank[idx] + opt.rank_slack >= p.target_rank[idx]
                            if min_eig(a, idx, p) > opt.tol_psd
                                if opt.rank_increment == 0
                                    p.target_rank[idx] = min(opt.rank_increment_factor * p.target_rank[idx], sdp.sq_side)
                                else
                                    p.target_rank[idx] = min(opt.rank_increment_factor + p.target_rank[idx], sdp.sq_side)
                                end
                            end
                        end
                    end
                    p.rank_update, p.update_cont = 0, 0
                end
            end
        # Check divergence
        elseif k > p.window && residuals.comb_residual[k - p.window] < residuals.comb_residual[k] && p.rank_update > p.window
            p.update_cont += 1
            if p.update_cont > opt.divergence_min_update
                full_rank_flag = true
                for (idx, sdp) in enumerate(conic_sets.sdpcone)
                    if p.target_rank[idx] < sdp.sq_side
                        full_rank_flag = false
                        p.rank_update, p.update_cont = 0, 0
                    end
                    if p.current_rank[idx] + opt.rank_slack >= p.target_rank[idx]
                        if min_eig(a, idx, p) > opt.tol_psd
                            if opt.rank_increment == 0
                                p.target_rank[idx] = min(opt.rank_increment_factor * p.target_rank[idx], sdp.sq_side)
                            else
                                p.target_rank[idx] = min(opt.rank_increment_factor + p.target_rank[idx], sdp.sq_side)
                            end
                        end
                    end
                end
            end
        # Adaptive stepsizes
        elseif residuals.primal_residual[k] > opt.tol_primal && residuals.dual_residual[k] < opt.tol_dual && k > p.window
            ada_count += 1
            if ada_count > opt.adapt_window
                ada_count = 0
                if opt.line_search_flag
                    p.beta *= (1. - p.adapt_level)
                    p.primal_step /= sqrt(1. - p.adapt_level)
                else
                    p.primal_step /= (1. - p.adapt_level)
                    p.dual_step *= (1. - p.adapt_level)
                end
                p.adapt_level *= opt.adapt_decay
            end
        elseif residuals.primal_residual[k] < opt.tol_primal && residuals.dual_residual[k] > opt.tol_dual && k > p.window
            ada_count += 1
            if ada_count > opt.adapt_window
                ada_count = 0
                if opt.line_search_flag
                    p.beta /= (1. - p.adapt_level)
                    p.primal_step *= sqrt(1. - p.adapt_level)
                else
                    p.primal_step *= (1. - p.adapt_level)
                    p.dual_step /= (1. - p.adapt_level)
                end
                p.adapt_level *= opt.adapt_decay
            end
        end

        # max_iter or time limit stop condition
        if p.iter >= opt.max_iter_local || time() - p.time0 >= opt.time_limit
            if p.iter > opt.min_iter_time_infeas &&
                max_abs_diff(residuals.dual_gap) < opt.infeas_stable_gap_tol &&
                residuals.dual_gap[k] > opt.infeas_limit_gap_tol # low gap but far from zero, say 10%
                if residuals.feasibility[p.iter] <= opt.tol_feasibility/100
                    p.stop_reason = 5 # Unbounded
                    p.stop_reason_string = "Problem declared unbounded due to lack of improvement"
                    if opt.certificate_search && !p.certificate_search
                        certificate_dual_infeasibility(affine_sets, p, opt)
                        push!(sol, cache_solution(pair, residuals, conic_sets, affine_sets, p, opt,
                            c_orig, A_orig, b_orig, G_orig, h_orig, D, E, var_ordering, a))
                        # println("6")
                    elseif opt.certificate_search && p.certificate_search
                        # error("6")
                    else
                        break
                    end
                elseif residuals.feasibility[p.iter] > opt.infeas_feasibility_tol
                    p.stop_reason = 6 # Infeasible
                    p.stop_reason_string = "Problem declared infeasible due to lack of improvement"
                    if opt.certificate_search && !p.certificate_search
                        certificate_infeasibility(affine_sets, p, opt)
                        push!(sol, cache_solution(pair, residuals, conic_sets, affine_sets, p, opt,
                            c_orig, A_orig, b_orig, G_orig, h_orig, D, E, var_ordering, a))
                        # println("7")
                    elseif opt.certificate_search && p.certificate_search
                        # error("7")
                    else
                        break
                    end
                end
            elseif p.iter >= opt.max_iter_local
                p.stop_reason = 3 # Iteration limit
                p.stop_reason_string = "Iteration limit of $(opt.max_iter_local) was hit"
                if opt.warn_on_limit
                    @warn("    WARNING: Iteration limit hit.")
                end
            else
                p.stop_reason = 2 # Time limit
                p.stop_reason_string = "Time limit hit, limit: $(opt.time_limit) time: $(time() - p.time0)"
                if opt.warn_on_limit
                    println("    WARNING: Time limit hit.")
                end
            end
            if p.iter >= opt.max_iter_local || time() - p.time0 >= opt.time_limit
                break
            end
        end

        # already looking of certificates
        if opt.certificate_search && p.certificate_search
            continue
        end

        # Dual obj growing too much
        if (p.iter > opt.min_iter_max_obj && residuals.dual_obj[k] >  opt.max_obj) || isnan(residuals.dual_obj[k])

            # Dual unbounded
            p.stop_reason = 6 # Infeasible
            p.stop_reason_string = "Infeasible: |Dual objective| = $(residuals.dual_obj[k]) > maximum allowed = $(opt.max_obj)"

            if opt.certificate_search && !p.certificate_search

                certificate_infeasibility(affine_sets, p, opt)
                push!(sol, cache_solution(pair, residuals, conic_sets, affine_sets, p, opt,
                    c_orig, A_orig, b_orig, G_orig, h_orig, D, E, var_ordering, a))

            else
                break
            end
        end

        # Primal obj growing too much
        if (p.iter > opt.min_iter_max_obj && residuals.prim_obj[k] < -opt.max_obj) || isnan(residuals.prim_obj[k])

            p.stop_reason = 5 # Unbounded
            p.stop_reason_string = "Unbounded: |Primal objective| = $(residuals.prim_obj[k]) > maximum allowed = $(opt.max_obj)"

            if opt.certificate_search && !p.certificate_search

                certificate_dual_infeasibility(affine_sets, p, opt)
                push!(sol, cache_solution(pair, residuals, conic_sets, affine_sets, p, opt,
                    c_orig, A_orig, b_orig, G_orig, h_orig, D, E, var_ordering, a))

            else
                break
            end
        end

        # Stalled feasibility with meaningful gap
        if (
            p.iter > opt.min_iter_max_obj &&
            residuals.dual_gap[k] > opt.infeas_limit_gap_tol && # Low gap but far from zero (~10%)
            residuals.feasibility[p.iter] > opt.infeas_feasibility_tol &&
            max_abs_diff(residuals.feasibility) < opt.infeas_stable_feasibility_tol
            )

            p.stop_reason = 6 # Infeasible
            p.stop_reason_string = "Infeasible: feasibility stalled at $(residuals.feasibility[p.iter])"

            if opt.certificate_search && !p.certificate_search

                certificate_infeasibility(affine_sets, p, opt)
                push!(sol, cache_solution(pair, residuals, conic_sets, affine_sets, p, opt,
                    c_orig, A_orig, b_orig, G_orig, h_orig, D, E, var_ordering, a))

            else
                break
            end
        end

        # Stalled gap at 100%
        if (
            p.iter > opt.min_iter_max_obj &&
            residuals.dual_gap[k] > 1-opt.infeas_gap_tol &&
            max_abs_diff(residuals.dual_gap) < opt.infeas_stable_gap_tol
            )

            if abs(residuals.dual_obj[k]) > abs(residuals.prim_obj[k]) && residuals.feasibility[p.iter] > opt.infeas_feasibility_tol

                # Dual unbounded
                p.stop_reason = 6 # Infeasible
                p.stop_reason_string = "Infeasible: duality gap stalled at 100 % with |Dual objective| >> |Primal objective|"

                if opt.certificate_search && !p.certificate_search

                    certificate_infeasibility(affine_sets, p, opt)
                    push!(sol, cache_solution(pair, residuals, conic_sets, affine_sets, p, opt,
                        c_orig, A_orig, b_orig, G_orig, h_orig, D, E, var_ordering, a))
                else
                    break
                end

            elseif abs(residuals.prim_obj[k]) > abs(residuals.dual_obj[k]) && residuals.feasibility[p.iter] <= opt.tol_feasibility

                p.stop_reason = 5 # Unbounded
                p.stop_reason_string = "Unbounded: duality gap stalled at 100 % with |Dual objective| << |Primal objective|"

                if opt.certificate_search && !p.certificate_search

                    certificate_dual_infeasibility(affine_sets, p, opt)
                    push!(sol, cache_solution(pair, residuals, conic_sets, affine_sets, p, opt,
                        c_orig, A_orig, b_orig, G_orig, h_orig, D, E, var_ordering, a))

                else
                    break
                end
            end
        end
    end

    # Compute results
    time_ = time() - p.time0

    # Print result
    if opt.log_verbose

        val = -1.0

        if opt.extended_log2
            cc = ifelse(p.stop_reason == 6, 0.0, 1.0)*c_orig
            val = dual_feas(pair.y, conic_sets, affine_sets, cc, A_orig, G_orig, a)
        end

        print_progress(residuals, p, opt, val)
        print_result(
            p.stop_reason,
            time_,
            residuals,
            length(p.current_rank) > 0 ? maximum(p.current_rank) : 0,
            p)
    end

    if opt.certificate_search && p.certificate_search

        @assert length(sol) == 1

        if p.certificate_found

            if p.stop_reason == 6
                c_orig .*= 0.0
            end

            pop!(sol)
            push!(sol, cache_solution(pair, residuals, conic_sets, affine_sets, p, opt,
                c_orig, A_orig, b_orig, G_orig, h_orig, D, E, var_ordering, a))
        end

    else
        @assert length(sol) == 0
        push!(sol, cache_solution(pair, residuals, conic_sets, affine_sets, p, opt,
        c_orig, A_orig, b_orig, G_orig, h_orig, D, E, var_ordering, a))
    end

    return sol[1]
end

function linesearch!(
    pair::PrimalDual,
    a::AuxiliaryData,
    affine_sets::AffineSets,
    mat::Matrices,
    opt::Options,
    p::Params
)::Nothing

    p.primal_step = p.primal_step * sqrt(1. + p.theta)

    for i in 1:opt.max_linsearch_steps
        p.theta = p.primal_step / p.primal_step_old

        @timeit "linesearch 1" begin
            a.y_half .= pair.y .+ (p.beta * p.primal_step) .* ((1. + p.theta) .* a.Mx .- p.theta .* a.Mx_old)
        end

        @timeit "linesearch 2" begin
            copyto!(a.y_temp, a.y_half)
            box_projection!(a.y_half, affine_sets, p.beta * p.primal_step)
            a.y_temp .-= (p.beta * p.primal_step) .* a.y_half
        end

        @timeit "linesearch 3" LinearAlgebra.mul!(a.Mty, mat.Mt, a.y_temp)

        # In-place norm
        @timeit "linesearch 4" begin
            a.Mty .-= a.Mty_old
            a.y_temp .-= pair.y_old
            y_norm = LinearAlgebra.norm(a.y_temp)
            Mty_norm = LinearAlgebra.norm(a.Mty)
        end

        if sqrt(p.beta) * p.primal_step * Mty_norm <= opt.delta * y_norm
            break
        else
            p.primal_step *= opt.linsearch_decay
        end
    end

    # Reverte in-place norm
    a.Mty .+= a.Mty_old
    a.y_temp .+= pair.y_old

    copyto!(pair.y, a.y_temp)
    p.primal_step_old = p.primal_step
    p.dual_step = p.beta * p.primal_step

    return nothing
end

function dual_step!(
    pair::PrimalDual,
    a::AuxiliaryData,
    affine_sets::AffineSets,
    mat::Matrices,
    opt::Options,
    p::Params
)::Nothing

    @timeit "dual step 1" begin
        a.y_half .= pair.y .+ p.dual_step * (2. * a.Mx .- a.Mx_old)
    end

    @timeit "dual step 2" begin
        copyto!(a.y_temp, a.y_half)
        box_projection!(a.y_half, affine_sets, p.dual_step)
        a.y_temp .-= p.dual_step * a.y_half
    end

    @timeit "linesearch 3" LinearAlgebra.mul!(a.Mty, mat.Mt, a.y_temp)

    copyto!(pair.y, a.y_temp)
    p.primal_step_old = p.primal_step

    return nothing
end

function primal_step!(
    pair::PrimalDual,
    a::AuxiliaryData,
    cones::ConicSets,
    mat::Matrices,
    opt::Options,
    p::Params,
    arc_list,
    iter::Int64,
)::Nothing

    pair.x .-= p.primal_step .* (a.Mty .+ mat.c)

    # Projection onto the p.s.d. cone
    if length(cones.sdpcone) >= 1
        @timeit "sdp proj" psd_projection!(pair.x, a, cones, opt, p, arc_list, iter)
    end

    # Projection onto the second order cone
    if length(cones.socone) >= 1
        @timeit "soc proj" soc_projection!(pair.x, a, cones, opt, p)
    end

    @timeit "linesearch -1" LinearAlgebra.mul!(a.Mx, mat.M, pair.x)

    return nothing
end

function certificate_dual_infeasibility(affine_sets, p, opt)

    if opt.log_verbose
        println("-"^87)
        println("    Begin search for dual infeasibility certificate")
        println("-"^87)
    end

    fill!(affine_sets.b, 0.0)
    fill!(affine_sets.h, 0.0)

    certificate_parameters(p, opt)

    return nothing
end

function certificate_infeasibility(affine_sets, p, opt)

    if opt.log_verbose
        println("-"^87)
        println("    Begin search for infeasibility certificate")
        println("-"^87)
    end

    fill!(affine_sets.c, 0.0)

    certificate_parameters(p, opt)

    return nothing
end

function certificate_parameters(p, opt)
    p.certificate_search_min_iter = p.iter + 2 * opt.convergence_window + div(p.iter, 5) + 1000
    p.certificate_search = true
    opt.time_limit *= 1.1
    opt.max_iter_local = opt.max_iter_local + div(opt.max_iter_local, 10)
    return nothing
end

function cone_feas(v, cones, a, num = sqrt(2))
    sdp_viol = 0.0
    sdplen = psd_vec_to_square(v, a, cones, num) - 1
    for (idx, sdp) in enumerate(cones.sdpcone)
        if sdp.sq_side == 1
            sdp_viol = max(sdp_viol, -min(0.0, a.m[idx][1]))
        else
            fact = LinearAlgebra.eigen!(a.m[idx])
            sdp_viol = max(sdp_viol, -min(0.0, minimum(fact.values)))
        end
    end
    soc_viol = 0.0
    cont = sdplen
    for (idx, soc) in enumerate(cones.socone)
        len = soc.len
        push!(a.soc_s, view(v, cont + 1))
        s = v[cont+1]
        sdp_viol = max(sdp_viol, -min(0.0, s - LinearAlgebra.norm(view(v, cont + 2:cont + len))))
        cont += len
    end
    return max(sdp_viol, soc_viol), cont
end

function get_duals(y::Vector{T}, cones::ConicSets, affine::AffineSets, c, A, G) where T

    dual_eq = y[1:affine.p]
    dual_in = y[affine.p+1:end]

    dual_cone = + c + A' * dual_eq + G' * dual_in
    fix_diag_scaling(dual_cone, cones, 2.0)

    return dual_eq, dual_in, dual_cone
end

function dual_feas(y::Vector{T}, cones::ConicSets, affine::AffineSets, c, A, G, a) where T
    dual_eq, dual_in, dual_cone = get_duals(y, cones, affine, c, A, G)
    return dual_feas(dual_in, dual_cone, cones, a)
end
function dual_feas(dual_in::Vector{T}, dual_cone::Vector{T}, cones, a) where T

    ineq_viol = 0.0
    if length(dual_in) > 0
        ineq_viol = -min(0.0, minimum(dual_in))
    end

    cone_viol, cont = cone_feas(dual_cone, cones, a)

    zero_viol = 0.0
    dual_zr = dual_cone[(cont+1):end]
    if length(dual_zr) > 0
        zero_viol = maximum(abs.(dual_zr))
    end

    return max(cone_viol, ineq_viol, zero_viol)
end

function fix_diag_scaling(v, cones, num)
    # Remove diag scaling
    cont = 1
    @inbounds for sdp in cones.sdpcone, j in 1:sdp.sq_side, i in 1:j#j:sdp.sq_side
        if i != j
            v[cont] /= num
        end
        cont += 1
    end
end

function cache_solution(pair, residuals, conic_sets, affine_sets, p, opt,
    c, A, b, G, h, D, E, var_ordering, a)

    # Remove diag scaling
    fix_diag_scaling(pair.x, conic_sets, sqrt(2.0))

    # Remove equilibrating
    if opt.equilibration
        pair.x = D * pair.x
        pair.y = E * pair.y
    end

    slack_eq   = A * pair.x - b
    slack_ineq = G * pair.x - h

    dual_eq, dual_in, dual_cone =
        get_duals(pair.y, conic_sets, affine_sets, c, A, G)

    dual_feasibility = dual_feas(dual_in, dual_cone, conic_sets, a)

    return Result(
        p.stop_reason,
        p.stop_reason_string,
        pair.x[var_ordering],
        dual_cone[var_ordering],
        dual_eq,
        dual_in,
        slack_eq,
        slack_ineq,
        residuals.equa_feasibility,
        residuals.ineq_feasibility,
        residuals.prim_obj[p.iter],
        residuals.dual_obj[p.iter],
        residuals.dual_gap[p.iter],
        time() - p.time0,
        p.iter,
        sum(p.current_rank),
        residuals.feasibility[p.iter] <= opt.tol_feasibility,
        dual_feasibility <= opt.tol_feasibility_dual,
        p.certificate_found,
        1,
    )
end


max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
abstract SQLite3 <: DBI.DatabaseSystem

type SQLiteDatabaseHandle <: DBI.DatabaseHandle
    ptr::Ptr{Void}
    status::Cint
end

type SQLiteStatementHandle <: DBI.StatementHandle
    db::SQLiteDatabaseHandle
    ptr::Ptr{Void}
    executed::Int

    function SQLiteStatementHandle(
        db::SQLiteDatabaseHandle,
        ptr::Ptr{Void}
    )
        new(db, ptr, 0)
    end
end


text:
using Bolt
using Test
using DelimitedFiles
using LinearAlgebra

@testset "FFTLog" begin
    N = 64
    μ = 0
    q = 0.0
    r₀ = 1.0
    L = 8.0
    Nhalf = N ÷ 2
    n = range(-Nhalf,Nhalf,length=N)
    r = r₀ .* 10 .^ (n .* L ./ N )
    pl = Bolt.plan_fftlog(r, μ, q, 1.0; kropt=true)
    aₙ = r .^ (μ + 1) .* exp.(-r.^2 / 2)
    y = similar(r, ComplexF64)
    fftdata = readdlm("data/fftlog_example.txt", ' ', Float64, '\n')

    # test forward
    mul!(y, pl, aₙ)
    f_ref = fftdata[:,2]
    @test all(abs.(y .- f_ref) .< 1e-15)
    @test isapprox(y, f_ref)

    # test backward
    y2 = similar(r, ComplexF64)
    ldiv!(y2, pl, y)
    @test all(abs.(y2 .- aₙ) .< 1e-15)
end

##
@testset "RECFAST" begin
    recfastdata = readdlm("data/test_recfast_1.dat", ',', Float64, '\n', header=true)[1]
    z⃗, Xe_fort = recfastdata[:,1], recfastdata[:,2]
    𝕡 = CosmoParams(Σm_ν=0.0, N_ν=3.0)
    bg = Background(𝕡)
    𝕣 = Bolt.RECFAST(bg=bg, OmegaB=𝕡.Ω_b, Yp=𝕡.Y_p)
    xe_bespoke, Tmat_bespoke = Bolt.recfast_xe(𝕣; Nz=1000, zinitial=10000., zfinal=0.)
    @test all(abs.(Xe_fort .- xe_bespoke) .< 1e-5)
end


text:
"""
    AccelDesang
An action type with a longitudinal acceleration and a desired heading angle

# Fields
- `a::Float64` longitudinal acceleration [m/s^2]
- `ϕdes::Float64` desired heading angle
"""
struct AccelDesang
    a::Float64 # accel [m/s²]
    ϕdes::Float64 # Desired heading angle
end
Base.show(io::IO, a::AccelDesang) = @printf(io, "AccelDesang(%6.3f,%6.3f)", a.a, a.ϕdes)
Base.length(::Type{AccelDesang}) = 2

"""
    Base.convert(::Type{AccelDesang}, v::Vector{Float64})
Convert a vector containing the acceleration and desired heading angle
into AccelDesang type
"""
Base.convert(::Type{AccelDesang}, v::Vector{Float64}) = AccelDesang(v[1], v[2])

"""
    Base.copyto!(v::Vector{Float64}, a::AccelDesang)
Extract the AccelDesang components into a vector and return the vector
"""
function Base.copyto!(v::Vector{Float64}, a::AccelDesang)
    v[1] = a.a
    v[2] = a.ϕdes
    v
end

"""
    propagate(veh::Entity{VehicleState,D,I}, action::AccelDesang, roadway::Roadway, Δt::Float64; n_integration_steps::Int=4) where {D,I}
Propagate vehicle forward in time using a desired acceleration
and heading angle
"""
function propagate(veh::Entity{VehicleState,D,I}, action::AccelDesang, roadway::Roadway, Δt::Float64; n_integration_steps::Int=4) where {D,I}

    a = action.a # accel
    ϕdes = action.ϕdes # desired heading angle

    x = veh.state.posG.x
    y = veh.state.posG.y
    θ = veh.state.posG.θ
    v = veh.state.v

    δt = Δt/n_integration_steps

    for i in 1 : n_integration_steps

        posF = Frenet(VecSE2(x, y, θ), roadway)
        ω = ϕdes - posF.ϕ

        x += v*cos(θ)*δt
        y += v*sin(θ)*δt
        θ += ω*δt
        v += a*δt
    end

    posG = VecSE2(x, y, θ)
    VehicleState(posG, roadway, v)
end


text:
# Download links in `assets.csv`

dir = joinpath(@__DIR__, "Pluto", "frontend", "offline_assets")
rm(dir, force=true, recursive=true)
mkpath(dir)

for url in readlines(joinpath(@__DIR__, "assets.csv") )
    @info "Downloading: $url"
    file = touch(joinpath(dir, basename(url)))
    download(url, file)
end


text:
# types used by Ising module

abstract type AbstractSampler end
abstract type AbstractSamplerState end


text:
include("lib_for_robot.jl")
function task14(r)
    for i in 0:3
        x=movetoside(r,HorizonSide(i))
        rtn(r,reverse(HorizonSide(i)),x)
    end
end


text:
"""
	loadatmacs(locfile,globfile)
(Download &) read Atmcas local and global data

**Input**
* glofiles: vector of either URLs or files (name) with global Atmacs effect.
	If locfile is empty, then global model covering whole Earth is assumed to
	be used thus including local part + pressure). All files will be stacked
	into one output for global DataFrame.
* locfiles: vector of either URLs or files (name) with local Atmacs effect
	(+ pressure)

**Output**
* globaal atmacs, local atmacs DataFrames with DateTime, and all columns of
  input files

**Example**
```
# Download data automatically
lurl = ["http://atmacs.bkg.bund.de/data/results/iconeu/we_iconeu_70km.grav"];
gurl = ["http://atmacs.bkg.bund.de/data/results/icon/we_icon384_20deg.grav"];
glo,loc = loadatmacs(locfiles=lurl,glofiles=gurl);
# Or read data from already downloaded files
glo,loc = loadatmacs(glofiles=["/test/input/atmacs_glo.grav"],
					locfiles=["/test/input/atmacs_loc.grav"]);
# Or stack global files covering global+local part (`loc` will be empty)
glo,loc = loadatmacs(glofiles=["/test/input/atmacs_all_1.grav",
								"/test/input/atmacs_all_2.grav"]);
```
"""
function loadatmacs(;glofiles::Vector{String}=[""],locfiles::Vector{String}=[""])
	# Constant output channel names: http://atmacs.bkg.bund.de/data/read_me
	# Assumption: Global file contains local part only if local file not loaded
	atmacs_loc = DataFrame();
	atmacs_glo = DataFrame();
	if isempty(locfiles[1])
		channelsGlo = [:pressure,:local_newton,:global_newton,:total_loading];
	else
		channelsGlo = [:global_newton,:total_loading];
		channelsLoc = [:pressure,:local_newton,:regional_newton];
		atmacs_loc = loadatmacsparts(locfiles,channelsLoc);
	end
	atmacs_glo = loadatmacsparts(glofiles,channelsGlo);
	return atmacs_glo,atmacs_loc;
end

"""
Auxiliary function for loading either global or local part of Atmacs
"""
function loadatmacsparts(files::Vector{String},channels::Vector{Symbol})
	atmacs = DataFrame();
	for i in files
		downto = FileTools.downfile(i);
		if !isempty(downto)
			dataall = readdlm(downto)
			temp = DataFrame(datetime=pattern2time.(dataall[:,1],"hour"));
			for j in 1:length(channels)
				temp[channels[j]] = dataall[:,j+1];
			end
			# delete downloaded file
			FileTools.delfile(i);
			# Stack data
			atmacs = stackframes(atmacs,temp,corroffset=true,maxval=NaN,maxtime=Dates.Hour(7));
		end
	end
	return atmacs;
end

"""
	stackframes(frameall,frameapp;maxtime,maxval)
Stack dataframas assuming identical number of columns but with unknown
overlapping or gap between DataFrame datetime vectors

**Input**
* frameall: main DataFrame that should be extended (must contain datetime column)
* frameapp: DataFrame to be appended at the end of `frameall`
* maxtime: maximum time difference allowed for stacking (if > than not stacked)
* corroffset: correct offset between individual columns of input DataFrame (see `maxval`)
* maxval: maximum allowed offset between values in columns for stacking (if > || `corroffset=false` than not stacked)
> can be either one absolute value (Float64) applied to all columns or vector with same number
of columns (-datetime) as input DataFrames. Will be only applied if overlapping
exists. Set to NaN for stacking without offset check.

**Output**
* stacked DataFrame with identical number of columns as `frameall`

**Example**
```
f1 = DataFrame(datetime=[DateTime(2010,1,1,0),DateTime(2010,1,1,1),
           DateTime(2010,1,1,2),DateTime(2010,1,1,3)],
		   grav=[10.,11.,12.,13.]);
f2 = DataFrame(datetime=[DateTime(2010,1,1,2),DateTime(2010,1,1,3),
           DateTime(2010,1,1,4),DateTime(2010,1,1,5)],
		   grav=[22.,23.,24.,25.]);
fall = stackframes(f1,f2,maxtime=Dates.Hour(1),maxval=NaN,corroffset=true);
# Following call will not be stacked as maximum offset is exceeded
fnot = stackframes(f1,f2,maxval=1.0)
```
"""
function stackframes(frameall::DataFrame,frameapp::DataFrame;
						maxtime=Dates.Year(9999),maxval=NaN,corroffset::Bool=true)
	# return current file if the total file is still empty (first run)
	if isempty(frameall)
		return frameapp;
	else
		# Check date (for overlapping or missing data)
    	r = findall(frameall[:datetime][end] .== frameapp[:datetime]);
	    if isempty(r) # No such time exist (no overlapping) => check how big is the gap
			# get temporal resolution (res) and gap between files (diff)
        	time_diff = frameapp[:datetime][1] - frameall[:datetime][end];
        	time_res = frameall[:datetime][end] - frameall[:datetime][end-1];
			# Start with spetial case where files such as Atmacs are stacked
			if (time_diff<=maxtime) && isnan(maxval) && corroffset
				return vcat(frameall,frameapp);
			elseif time_diff <= maxtime
				temp = DataFrame(datetime = frameall[:datetime][end]+time_res)
				for i in names(frameapp)
					if i != :datetime
						temp[i] = NaN;
					end
				end
	        	return vcat(frameall,temp,frameapp);
			else
				return frameall;
			end
    	else # Overlapping exists
			# First, prepare temporary variable with `frameapp` values
			temp = DataFrame(datetime = frameapp[:datetime][r[1]+1:end])
			stackornot = true;
			# Prepare maximum allowed offset: convert to vector
			if typeof(maxval) == Float64
				maxoffset = zeros(size(frameall,2)-1) .+ maxval;
			end
			# and set to infinity if NaN on input
			maxoffset[isnan.(maxoffset)] .= Inf;
			# Apply offset (if set so)
			for (j,i) in enumerate(names(frameapp))
				if i != :datetime
					offset = frameall[i][end] - frameapp[i][r[1]];
					if abs(offset) > maxoffset[j-1]
						stackornot = false;
					elseif corroffset
						temp[i] = frameapp[i][r[1]+1:end] .+ offset;
					else
						temp[i] = frameapp[i][r[1]+1:end];
					end
				end
			end
			if stackornot
				return vcat(frameall,temp);
			else
				return frameall;
			end
		end
    end
end

"""
Auxiliary Function to convert time pattern (e.g., yyyymmdd) to datetime value
"""
function pattern2time(time::Float64,resolution::String)
	# Declare output
	year,month,day,hour,minute,second = (0,0,0,0,0,0)
	if resolution == "month"
		year = floor(Int,time/1e+2);
		month = floor(Int,time-year*1e+2);
		day = 1;
	elseif resolution == "day"
	    year = floor(Int,time/1e+4);
	    month= floor(Int,(time - year*1e+4)/1e+2);
	    day  = floor(Int,(time - year*1e+4 - month*1e+2));
	elseif resolution == "hour"
	    year = floor(Int,time/1e+6);
	    month= floor(Int,(time - year*1e+6)/1e+4);
	    day  = floor(Int,(time - year*1e+6 - month*1e+4)/1e+2);
	    hour = floor(Int,(time - year*1e+6 - month*1e+4 - day*1e+2));
	elseif resolution == "minute"
	    year = floor(Int,time/1e+8);
	    month= floor(Int,(time - year*1e+8)/1e+6);
	    day  = floor(Int,(time - year*1e+8 - month*1e+6)/1e+4);
	    hour = floor(Int,(time - year*1e+8 - month*1e+6 - day*1e+4)/1e+2);
	    minute=floor(Int,(time - year*1e+8 - month*1e+6 - day*1e+4 - hour*1e+2));
	elseif resolution == "second"
	    year = floor(Int,time/1e+10);
	    month= floor(Int,(time - year*1e+10)/1e+8);
	    day  = floor(Int,(time - year*1e+10 - month*1e+8)/1e+6);
	    hour = floor(Int,(time - year*1e+10 - month*1e+8 - day*1e+6)/1e+4);
	    minute=floor(Int,(time - year*1e+10 - month*1e+8 - day*1e+6 - hour*1e+4)/1e+2);
	    second=round(Int,time - year*1e+10 - month*1e+8 - day*1e+6 - hour*1e+4 - minute*1e+2);
	end
	return DateTime(year,month,day,hour,minute,second);
end
function pattern2time(time::Int,resolution::String)
	pattern2time(float(time),resolution);
end


text:
@testset "normalize_par_angles" begin
    @test normalize_par_angles([-30, -20, -10]) == [330, 340, 350]
    @test normalize_par_angles([-10, 50, 51]) == [350, 410, 411]
    @test normalize_par_angles([20, 40, 60]) == [20, 40, 60]

    x = [-10, 50, 51]
    normalize_par_angles!(x)
    @test x == [350, 410, 411]
end


text:
##### Beginning of file

# Parts of this file are based on:
# 1. https://github.com/JuliaPackaging/Git.jl/blob/master/deps/build.jl

import Conda

function _default_git_cmd()::String
    result::String = lowercase(strip("git"))
    return result
end

function _found_default_git()::Bool
    default_git_cmd::String = _default_git_cmd()
    found_default_git::Bool = try
        success(`$(default_git_cmd) --version`)
    catch
        false
    end
    return found_default_git
end

function _install_git()::String
    result::String = _install_git_conda()
    return result
end

function _install_git_conda()::String
    @info("Attempting to install Git using Conda.jl...")
    environment::Symbol = :MirrorUpdater
    Conda.add("git", environment)
    @info("Successfully installed Git using Conda.jl.")
    git_cmd::String = strip(
        joinpath(
            Conda.bin_dir(environment),
            "git",
            )
        )
    run(`$(git_cmd) --version`)
    return git_cmd
end

function _build_git()::String
    install_git::Bool = lowercase(strip(get(ENV, "INSTALL_GIT", "false"))) ==
        lowercase(strip("true"))
    found_default_git::Bool = _found_default_git()
    if install_git
        @info("INSTALL_GIT is true, so I will now install git.")
        git_cmd = _install_git()
    elseif found_default_git
        @info("I found git on your system, so I will use that git.")
        git_cmd = _default_git_cmd()
    else
        @info("I did not find git on your system, so I will now install git.")
        git_cmd = _install_git()
    end
    return git_cmd
end

function _build_mirrorupdater()::Nothing
    git_cmd = _build_git()
    build_jl_file_path = strip(
        abspath(
            strip(
                @__FILE__
                )
            )
        )
    @debug(
        "deps/build.jl: ",
        build_jl_file_path,
        )
    deps_directory = strip(
        abspath(
            strip(
                dirname(
                    strip(
                        build_jl_file_path
                        )
                    )
                )
            )
        )
    @debug(
        "deps:",
        deps_directory,
        )
    deps_jl_file_path = strip(
        abspath(
            joinpath(
                strip(deps_directory),
                strip("deps.jl"),
                )
            )
        )
    @debug(
        "deps/deps.jl:",
        deps_jl_file_path,
        )
    open(deps_jl_file_path, "w") do f
        line_1::String = "git_cmd = \"$(strip(string(git_cmd)))\""
        @info("Writing line 1 to deps.jl: ", line_1,)
        println(f, line_1)
    end
    return nothing
end

_build_mirrorupdater()

##### End of file


text:
export LayerNorm
using Statistics
using Knet.Train20: param

# TODO: this is slow, need a kernel, maybe https://github.com/tensorflow/tensorflow/pull/6205/files
# TODO: support other axes

"""
    LayerNorm(size::Integer; eps=1e-5)
    LayerNorm(γ, β, ϵ)

References:
* [Ba, Kiros and Hinton 2016](https://arxiv.org/abs/1607.06450) Layer Normalization
* torch.nn.LayerNorm
* tf.keras.layers.LayerNormalization
"""
struct LayerNorm; γ; β; ϵ; end

function LayerNorm(dmodel; eps=1e-5)
    γ = param(dmodel; init=ones)
    β = param(dmodel; init=zeros)
    LayerNorm(γ, β, eps)
end

function (l::LayerNorm)(x, o...)
    μ = mean(x,dims=1)
    σ = std(x,mean=μ,dims=1)
    ϵ = eltype(x)(l.ϵ)
    l.γ .* (x .- μ) ./ (σ .+ ϵ) .+ l.β # TODO: doing x .- μ twice?
end

# function (l::LayerNorm)(x::MaskedArray, o...)
#     MaskedArray(l(x.array), x.mask) # TODO: shouldn't normalization ignore masked values?
# end


text:
using ModelingToolkit

# Bounds
@variables u [bounds = (-1, 1)]
@test getbounds(u) == (-1, 1)
@test hasbounds(u)

@variables y
@test !hasbounds(y)

# Disturbance
@variables u [disturbance = true]
@test isdisturbance(u)

@variables y
@test !isdisturbance(y)

# Tunable
@parameters u [tunable = true]
@test istunable(u)

@parameters y
@test !istunable(y)

# Distributions
struct FakeNormal end
d = FakeNormal()
@parameters u [dist = d]
@test hasdist(u)
@test getdist(u) == d

@parameters y
@test !hasdist(y)

## System interface
@parameters t
Dₜ = Differential(t)
@variables x(t)=0 u(t)=0 [input = true] y(t)=0 [output = true]
@parameters T [tunable = true, bounds = (0, Inf)]
@parameters k [tunable = true, bounds = (0, Inf)]
@parameters k2
eqs = [Dₜ(x) ~ (-k2 * x + k * u) / T
       y ~ x]
sys = ODESystem(eqs, t, name = :tunable_first_order)

p = tunable_parameters(sys)
sp = Set(p)
@test k ∈ sp
@test T ∈ sp
@test k2 ∉ sp
@test length(p) == 2

lb, ub = getbounds(p)
@test lb == [0, 0]
@test ub == [Inf, Inf]

b = getbounds(sys)
@test b[T] == (0, Inf)

p = tunable_parameters(sys, default = true)
sp = Set(p)
@test k ∈ sp
@test T ∈ sp
@test k2 ∈ sp
@test length(p) == 3


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

module Display

using UUIDs
import LibGit2

using ..Types

const PackageEntry = Types.PackageEntry

const colors = Dict(
    ' ' => :white,
    '+' => :light_green,
    '-' => :light_red,
    '↑' => :light_yellow,
    '~' => :light_yellow,
    '↓' => :light_magenta,
    '?' => :red,
)
const color_dark = :light_black

function git_file_stream(repo::LibGit2.GitRepo, spec::String; fakeit::Bool=false)::IO
    blob = try LibGit2.GitBlob(repo, spec)
    catch err
        err isa LibGit2.GitError && err.code == LibGit2.Error.ENOTFOUND || rethrow()
        fakeit && return devnull
    end
    iob = IOBuffer(LibGit2.content(blob))
    close(blob)
    return iob
end

function status(ctx::Context, pkgs::Vector{PackageSpec}=PackageSpec[];
                diff::Bool=false, mode::PackageMode=PKGMODE_PROJECT, use_as_api=false)
    env = ctx.env
    project₀ = project₁ = env.project
    manifest₀ = manifest₁ = env.manifest
    mdiff = nothing

    pkgfilter = (entry) -> begin
        for pkg in pkgs
            if pkg.uuid == entry.uuid
                return true
            end
            if mode == PKGMODE_MANIFEST || mode == PKGMODE_COMBINED
                # also look for pkg's deps
                info = manifest_info(ctx.env, pkg.uuid)
                for uuid in values(info.deps)
                    if uuid == entry.uuid
                        return true
                    end
                end
            end
        end
        return false
    end
    filter_pkgs = length(pkgs) > 0

    if !use_as_api
        pkg = ctx.env.pkg
        if pkg !== nothing
           printstyled("Project "; color=Base.info_color(), bold=true)
           println(pkg.name, " v", pkg.version)
        end
    end
    if diff
        if env.git === nothing
            @warn "diff option only available for environments in git repositories, ignoring."
        else # env.git !== nothing
            try
                LibGit2.with(LibGit2.GitRepo(env.git)) do repo
                    git_path = LibGit2.path(repo)
                    project_path = relpath(env.project_file, git_path)
                    manifest_path = relpath(env.manifest_file, git_path)
                    project₀ = read_project(git_file_stream(repo, "HEAD:$project_path", fakeit=true))
                    manifest₀ = read_manifest(git_file_stream(repo, "HEAD:$manifest_path", fakeit=true))
                end
            catch
                @warn "could not read project from HEAD, displaying absolute status instead."
            end
        end
    end
    if mode == PKGMODE_PROJECT || mode == PKGMODE_COMBINED
        # TODO: handle project deps missing from manifest
        m₀ = Dict(uuid => entry for (uuid, entry) in manifest₀ if (uuid in values(project₀.deps)))
        m₁ = Dict(uuid => entry for (uuid, entry) in manifest₁ if (uuid in values(project₁.deps)))
        mdiff = manifest_diff(ctx, m₀, m₁)
        filter_pkgs && filter!(pkgfilter, mdiff)
        if !use_as_api
            printpkgstyle(ctx, :Status, pathrepr(env.project_file), #=ignore_indent=# true)
            print_diff(ctx, mdiff, #=status=# true)
        end
    end
    if mode == PKGMODE_MANIFEST
        mdiff = manifest_diff(ctx, manifest₀, manifest₁)
        filter_pkgs && filter!(pkgfilter, mdiff)
        if !use_as_api
            printpkgstyle(ctx, :Status, pathrepr(env.manifest_file), #=ignore_indent=# true)
            print_diff(ctx, mdiff, #=status=# true)
        end
    elseif mode == PKGMODE_COMBINED
        combined = values(merge(project₀.deps, project₁.deps))
        m₀ = Dict(uuid => entry for (uuid, entry) in manifest₀ if (uuid in combined))
        m₁ = Dict(uuid => entry for (uuid, entry) in manifest₁ if (uuid in combined))
        c_diff = filter!(x->x.old != x.new, manifest_diff(ctx, m₀, m₁))
        filter_pkgs && filter!(pkgfilter, c_diff)
        if !isempty(c_diff)
            if !use_as_api
                printpkgstyle(ctx, :Status, pathrepr(env.manifest_file), #=ignore_indent=# true)
                print_diff(ctx, c_diff, #=status=# true)
            end
            mdiff = Base.vcat(c_diff, mdiff)
        end
    end
    return mdiff
end

function print_project_diff(ctx::Context, env0::EnvCache, env1::EnvCache)
    pm0 = Dict(uuid => entry for (uuid, entry) in env0.manifest if (uuid in values(env0.project.deps)))
    pm1 = Dict(uuid => entry for (uuid, entry) in env1.manifest if (uuid in values(env1.project.deps)))
    diff = filter!(x->x.old != x.new, manifest_diff(ctx, pm0, pm1))
    if isempty(diff)
        printstyled(color = color_dark, " [no changes]\n")
    else
        print_diff(ctx, diff)
    end
end

function print_manifest_diff(ctx::Context, env₀::EnvCache, env₁::EnvCache)
    diff = manifest_diff(ctx, env₀.manifest, env₁.manifest)
    diff = filter!(x->x.old != x.new, diff)
    if isempty(diff)
        printstyled(color = color_dark, " [no changes]\n")
    else
        print_diff(ctx, diff)
    end
end

struct VerInfo
    hash::Union{SHA1,Nothing}
    path::Union{String,Nothing}
    ver::Union{VersionNumber,Nothing}
    pinned::Bool
    repo::Union{Types.GitRepo, Nothing}
end

revstring(str::String) = occursin(r"\b([a-f0-9]{40})\b", str) ? str[1:7] : str

vstring(ctx::Context, a::VerInfo) =
    string((a.ver == nothing && a.hash != nothing) ? "[$(string(a.hash)[1:16])]" : "",
           a.ver != nothing ? "v$(a.ver)" : "",
           a.path != nothing ? " [$(pathrepr(a.path))]" : "",
           a.repo != nothing ? " #$(revstring(a.repo.rev)) ($(a.repo.url))" : "",
           a.pinned == true ? " ⚲" : "",
           )

Base.:(==)(a::VerInfo, b::VerInfo) =
    a.hash == b.hash && a.ver == b.ver && a.pinned == b.pinned && a.repo == b.repo &&
    a.path == b.path

≈(a::VerInfo, b::VerInfo) = a.hash == b.hash &&
    (a.ver == nothing || b.ver == nothing || a.ver == b.ver) &&
    (a.pinned == b.pinned) &&
    (a.repo == nothing || b.repo == nothing || a.repo == b.repo) &&
    (a.path == b.path)

struct DiffEntry
    uuid::UUID
    name::String
    old::Union{VerInfo,Nothing}
    new::Union{VerInfo,Nothing}
end

function print_diff(io::IO, ctx::Context, diff::Vector{DiffEntry}, status=false)
    same = all(x.old == x.new for x in diff)
    some_packages_not_downloaded = false
    for x in diff
        pkgid = Base.PkgId(x.uuid, x.name)
        package_downloaded = pkgid in keys(Base.loaded_modules) ||
                             Base.locate_package(pkgid) !== nothing
        # Package download detection doesnt work properly when runn running targets
        ctx.currently_running_target && (package_downloaded = true)
        if x.old !== nothing && x.new !== nothing
            if x.old ≈ x.new
                verb = ' '
                vstr = vstring(ctx, x.new)
            else
                if (x.old.hash === nothing || x.new.hash === nothing || x.old.hash != x.new.hash) &&
                    x.old.ver != x.new.ver
                    verb = x.old.ver == nothing || x.new.ver == nothing ||
                           x.old.ver == x.new.ver ? '~' :
                           x.old.ver < x.new.ver  ? '↑' : '↓'
                elseif x.old.ver == x.new.ver && x.old.pinned != x.new.pinned ||
                    x.old.path != x.new.path ||
                    x.old.repo != nothing || x.new.repo != nothing
                    verb = '~'
                else
                    verb = '?'
                end
                vstr = (x.old.ver == x.new.ver && x.old.pinned == x.new.pinned &&
                        x.old.repo == x.new.repo && x.old.path == x.new.path) ?
                      vstring(ctx, x.new) :
                      vstring(ctx, x.old) * " ⇒ " * vstring(ctx, x.new)
            end
        elseif x.new != nothing
            verb = '+'
            vstr = vstring(ctx, x.new)
        elseif x.old != nothing
            verb = '-'
            vstr = vstring(ctx, x.old)
        else
            verb = '?'
            vstr = "[unknown]"
        end
        v = same ? "" : " $verb"
        if verb != '-' && status && !package_downloaded
            printstyled(io, "→", color=:red)
        else
            print(io, " ")
        end
        if verb != '-'
            package_downloaded || (some_packages_not_downloaded = true)
        end
        printstyled(io, " [$(string(x.uuid)[1:8])]"; color = color_dark)
        printstyled(io, "$v $(x.name) $vstr\n"; color = colors[verb])
    end
    if isempty(diff)
        str = ctx.env.git === nothing ? "  (empty environment)\n" : "  (no changes since last commit)\n"
        printstyled(str, color = color_dark)
    end
    if status && some_packages_not_downloaded
        @warn "Some packages (indicated with a red arrow) are not downloaded, use `instantiate` to instantiate the current environment"
    end
end
# TODO: Use the Context stream
print_diff(ctx::Context, diff::Vector{DiffEntry}, status=false) = print_diff(stdout, ctx, diff, status)

function name_ver_info(entry::PackageEntry)
    entry.name, VerInfo(
        entry.tree_hash,
        entry.path,
        entry.version,
        entry.pinned,
        entry.repo.url === nothing ? nothing : entry.repo, # TODO
        )
end

function manifest_diff(ctx::Context, manifest0::Dict, manifest1::Dict)
    diff = DiffEntry[]
    for uuid in union(keys(manifest0), keys(manifest1))
        name₀ = name₁ = v₀ = v₁ = nothing
        haskey(manifest0, uuid) && ((name₀, v₀) = name_ver_info(manifest0[uuid]))
        haskey(manifest1, uuid) && ((name₁, v₁) = name_ver_info(manifest1[uuid]))
        name₀ === nothing && (name₀ = name₁)
        name₁ === nothing && (name₁ = name₀)
        if name₀ == name₁
            push!(diff, DiffEntry(uuid, name₀, v₀, v₁))
        else
            push!(diff, DiffEntry(uuid, name₀, v₀, nothing))
            push!(diff, DiffEntry(uuid, name₁, nothing, v₁))
        end
    end
    sort!(diff, by=x->(x.uuid in keys(ctx.stdlibs), x.name, x.uuid))
end

end # module


text:
"""Projection of Ω to subspace

Conceptually if ω::Ω represents the unit hypercube id -> [0, 1], an ωπ::ΩProj
represets the element `ω(id)`.  However:
- The application is lazy to actually get element in `[0, 1]`, we use `memrand`
- The application is lazy and can be 'undone', ω can be recovered with parenωπ(ωπ)
- ΩProj supports rand which (i) resolves ωπ at its index, (ii) increments the id
"""
mutable struct ΩProj{O, I} <: Ω{I}
  ω::O
  id::I
end

increment!(ωπ::ΩProj) = ωπ.id = increment(ωπ.id)
parentω(ωπ::ΩProj) = ωπ.ω

# "memrand ωπ::ΩProj to a concrete value"
# function memrand end

@spec rand(ωπ::ΩProj, args...) _res == memrand(_pre(ωπ), args...) "result is resolution of ωπ"
@spec rand(ωπ::ΩProj, args...) _res == ωπ.id == increment(_pre(ωπ.id)) "ωπ id is incremented"

@inline randinc!(ωπ, res) = (increment!(ωπ); res)

# Try Again #
@inline Base.rand(ωπ::ΩProj; rng = Random.GLOBAL_RNG) =
  randinc!(ωπ, memrand(ωπ.ω, ωπ.id, Float64; rng = rng))
@inline Base.rand(ωπ::ΩProj, dims::Dims; rng = Random.GLOBAL_RNG) =
  randinc!(ωπ, memrand(ωπ.ω, ωπ.id, Float64, dims; rng = rng))
@inline Base.rand(ωπ::ΩProj, dim::Integer, dims::Integer...; rng = Random.GLOBAL_RNG) =
  randinc!(ωπ, memrand(ωπ.ω, ωπ.id, Float64, Dims((dim, dims...)); rng = rng))

@inline Base.rand(ωπ::ΩProj, ::Type{T}; rng = Random.GLOBAL_RNG) where T =
  randinc!(ωπ, memrand(ωπ.ω, ωπ.id, T; rng = rng))
@inline Base.rand(ωπ::ΩProj, ::Type{T}, dims::Dims; rng = Random.GLOBAL_RNG) where T =
  randinc!(ωπ, memrand(ωπ.ω, ωπ.id, T, dims; rng = rng))
@inline Base.rand(ωπ::ΩProj, ::Type{T}, dim::Integer, dims::Integer...; rng = Random.GLOBAL_RNG) where T=
  randinc!(ωπ, memrand(ωπ.ω, ωπ.id, T, Dims((dim, dims...)); rng = rng))

@inline Base.rand(ωπ::ΩProj, arr::Array; rng = Random.GLOBAL_RNG) =
  randinc!(ωπ, memrand(ωπ.ω, ωπ.id, arr; rng = rng))
@inline Base.rand(ωπ::ΩProj, arr::Array, dims::Dims; rng = Random.GLOBAL_RNG) =
  randinc!(ωπ, memrand(ωπ.ω, ωπ.id, arr, dims; rng = rng))
@inline Base.rand(ωπ::ΩProj, arr::Array, dim::Integer, dims::Integer...; rng = Random.GLOBAL_RNG) =
  randinc!(ωπ, memrand(ωπ.ω, ωπ.id, arr, Dims((dim, dims...)); rng = rng))

@inline Base.rand(ωπ::ΩProj, ur::UnitRange; rng = Random.GLOBAL_RNG) =
  randinc!(ωπ, memrand(ωπ.ω, ωπ.id, ur; rng = rng))
@inline Base.rand(ωπ::ΩProj, ur::UnitRange, dims::Dims; rng = Random.GLOBAL_RNG) =
  randinc!(ωπ, memrand(ωπ.ω, ωπ.id, ur, dims; rng = rng))
@inline Base.rand(ωπ::ΩProj, ur::UnitRange, dim::Integer, dims::Integer...; rng = Random.GLOBAL_RNG) =
  randinc!(ωπ, memrand(ωπ.ω, ωπ.id, ur, Dims((dim, dims...)); rng = rng))

# Fail Again

# memrand(ωπ::ΩProj, args...) = memrand(ωπ.ω, ωπ.id, args...)

# Ideally
# Base.rand(ωπ::ΩProj, args...; rng = Random.GLOBAL_RNG) =
#   (res = memrand(ωπ.ω, ωπ.id, args...; rng = rng); increment!(ωπ); res)

# e.g. rand((1,2,3))
# Base.rand(ωπ::ΩProj, dims::Dims; rng = Random.GLOBAL_RNG) =
#   (res = memrand(ωπ.ω, ωπ.id, dims; rng = rng); increment!(ωπ); res)

# # e.g. rand(1,2,3)
# Base.rand(ωπ::ΩProj, dims::Integer...; rng = Random.GLOBAL_RNG) =
#   (res = memrand(ωπ.ω, ωπ.id, dims...; rng = rng); increment!(ωπ); res)

# # Base.rand(ωπ::ΩProj; rng = rng) =
# #   (res = memrand(ωπ.ω, dims...); increment!(ωπ); res)

# # e.g: rand(Float64)
# Base.rand(ωπ::ΩProj, ::Type{X} = Float64; rng = Random.GLOBAL_RNG) where {X} =
#   (res = memrand(ωπ.ω, ωπ.id, X; rng = rng); increment!(ωπ); res)

# Projection #

Base.getindex(ωπ::ΩProj{O, I}, i::I) where {O, I} =
  ΩProj{O, I}(ωπ.ω, combine(ωπ.id, i))

Base.getindex(ωπ::ΩProj{O, I}, i::SI) where {O, I, SI} =
  ΩProj{O, I}(ωπ.ω, append(ωπ.id, i))

Base.getindex(ωπ::ΩProj{O, Paired}, i::Int) where O = ΩProj{O, Paired}(ωπ.ω, pair(ωπ.id, i))

function Base.getindex(sω::O, i::Int) where {I, O <: ΩBase{I}}
  ΩProj{O, I}(sω, base(I, i))
end


text:
"Z-normalize the input array"
function znorm(x::AbstractVector)
    x = x .- mean(x)
    x ./= std(x, mean=0, corrected=false)
end

function znorm!(x::AbstractVector)
    x .= x .- mean(x)
    x ./= std(x, mean=0, corrected=false)
end

"""
    meanstd(x::AbstractArray)

Return mean and std of `x`. Is faster but less accurate than built-in functions.
"""
function meanstd(x::AbstractArray{T}) where T
    s = ss = zero(float(T))
    n = length(x)
    @avx for i in eachindex(x)
        s  += x[i]
        ss += x[i]^2
    end
    m   = s/n
    sig = √(ss/n - m^2)
    m, sig
end

import Base: @propagate_inbounds, @boundscheck, getindex
import LinearAlgebra: normalize

abstract type AbstractNormalizer{T,N} <: AbstractArray{T,N} end

advance!(x) = 0 # Catches everything that's not a normalizer

setup_normalizer(n::Type{Nothing}, q, y) = q, y

normalize(::Type{Nothing}, q) = q

(z::AbstractNormalizer)(q) = normalize(z, q)
(::Type{N})(q) where N <: AbstractNormalizer = normalize(N, q)

abstract type AbstractZNormalizer{T,N} <: AbstractNormalizer{T,N} end

"""
    ZNormalizer{T} <: AbstractNormalizer{T, 1}

Utility object that normalizes an input array on the fly.
Works like a vector, index into window with `z[i]`, index into normalized window with `z[!,i]`.

# Arguments:
- `x::Vector{T}`: The vecetor to operate on
- `n::Int`: The lenght of a window
- `μ::T`: mean over window
- `σ::T`: std over window
- `s::T`: sum over window
- `ss::T`: sum^2 over window
- `i::Int`: the current index
- `buffer::Vector{T}`: stores calculated normalized values, do not access this, it might not be complete
- `bufi::Int`: index of last normalized value
"""
mutable struct ZNormalizer{T,N} <: AbstractZNormalizer{T,N}
    x::Array{T,N}
    n::Int
    μ::T
    σ::T
    s::T
    ss::T
    i::Int
    buffer::Array{T,N}
    bufi::Int
end

function ZNormalizer(x::AbstractVector{T}, n) where T
    @assert length(x) >= n
    s = ss = zero(T)
    @avx for i in 1:n
        s += x[i]
        ss += x[i]^2
    end
    μ = s/n
    σ = sqrt(max(0, ss/n - μ^2))
    buffer = similar(x, n)
    ZNormalizer(x, n, μ, σ, s, ss, 0, buffer, 0)
end

function ZNormalizer(x::AbstractMatrix{T}, n) where T
    @assert length(x) >= n
    m = size(x,1)*n
    s = ss = zero(T)
    for i in 1:n
        s += sum(x[!,i])
        ss += sum(abs2, x[!,i])
    end
    μ = s/m
    σ = sqrt(max(0, ss/m - μ^2))
    buffer = similar(x, size(x,1), n)
    ZNormalizer(x, n, μ, σ, s, ss, 0, buffer, 0)
end

function normalize(::Type{ZNormalizer}, q::AbstractArray)
    q = q .- mean(q)
    q ./= std(q, corrected=false, mean=0)
end

setup_normalizer(z::Type{ZNormalizer}, q::AbstractVecOrMat, y::AbstractVecOrMat) = normalize(z, q), ZNormalizer(y, lastlength(q))

@propagate_inbounds function advance!(z::ZNormalizer{T,1}) where T
    if z.i == 0
        return z.i = 1
    end
    @boundscheck if z.i + z.n > length(z.x)
        return z.i += 1
    end
    z.bufi = 0

    # Remove old point
    x = z.x[z.i]
    z.s -= x
    z.ss -= x^2

    # Add new point
    x = z.x[z.i+z.n]
    z.s += x
    z.ss += x^2
    z.μ = z.s/z.n
    z.σ = sqrt(max(z.ss/z.n - z.μ^2, 0))
    z.i += 1
end

@propagate_inbounds function advance!(z::ZNormalizer{T,2}) where T
    if z.i == 0
        return z.i = 1
    end
    @boundscheck if z.i + z.n > lastlength(z.x)
        return z.i += 1
    end
    z.bufi = 0

    # Remove old point
    x = z.x[!, z.i]
    z.s -= sum(x)
    z.ss -= sum(abs2, x)

    # Add new point
    x = z.x[!, z.i+z.n]
    z.s += sum(x)
    z.ss += sum(abs2, x)
    m = length(z)
    z.μ = z.s/m
    z.σ = sqrt(max(z.ss/m - z.μ^2, 0))
    z.i += 1
end

@inline @propagate_inbounds function getindex(z::ZNormalizer{<:Any, 1}, i)
    @boundscheck 1 <= i <= z.n || throw(BoundsError(z,i))
    xi = i+z.i-1
    @boundscheck xi <= length(z.x) || throw(BoundsError(z,i))
    z.x[xi]
end

@inline @propagate_inbounds function getindex(z::ZNormalizer{<:Any, 1}, ::typeof(!), i::Int, inorder = i == z.bufi + 1)
    y = (z[i]-z.μ) / z.σ
    if inorder
        z.bufi = i
        z.buffer[i] = y
    end
    y
end

@inline @propagate_inbounds function getindex(z::ZNormalizer{<:Any, 1}, ::typeof(!), i::AbstractRange)
    @boundscheck (i[1] == z.i && length(i) == z.n) || throw(ArgumentError("ZNormalizers can only be indexed by ranges corresponding to their current state. Got range $i but state was $(z.i) corresponding to range $(z.i):$(z.i+z.n-1)"))
    z
end

@inline @propagate_inbounds function getindex(z::ZNormalizer{T,2}, ::typeof(!), i::Int, inorder = i == z.bufi + 1) where T
    j = inorder ? i : z.n
    xj = z.i + i - 1
    μ,σ = z.μ, z.σ + eps(T)
    @avx for k = 1:size(z.x, 1)
        z.buffer[k, j] = (z.x[k, xj] - μ) / σ
    end
    if inorder
        z.bufi = i
    end
    z.buffer[!, j]
end

Statistics.mean(z::AbstractZNormalizer) = z.μ
Statistics.std(z::AbstractZNormalizer) = z.σ

SlidingDistancesBase.lastlength(z::AbstractZNormalizer) = z.n
Base.length(z::AbstractZNormalizer) = length(z.buffer)
Base.size(z::ZNormalizer, args...) = size(z.buffer, args...)
actuallastlength(x::AbstractZNormalizer) = lastlength(x.x)

# Multi dim ==================================================================================
mutable struct DiagonalZNormalizer{T} <: AbstractZNormalizer{T,2}
    x::Array{T,2}
    n::Int
    μ::Array{T,1}
    σ::Array{T,1}
    s::Array{T,1}
    ss::Array{T,1}
    i::Int
    buffer::Array{T,2}
    bufi::Int
end

function DiagonalZNormalizer(x::AbstractArray{T,2}, n) where T
    @assert length(x) >= n
    s  = zeros(T, size(x,1))
    ss = zeros(T, size(x,1))
    @inbounds @simd for i in 1:n
        s  .+= x[!, i]
        ss .+= x[!, i].^2
    end
    μ = s./n
    σ = sqrt.(max.(ss./n .- μ.^2, 0))
    buffer = similar(x, size(x,1), n)
    DiagonalZNormalizer(x, n, μ, σ, s, ss, 0, buffer, 0)
end

function normalize(::Type{DiagonalZNormalizer}, q::AbstractMatrix)
    q = q .- mean(q, dims=2) # TODO: this will cause a ton of allocations
    q ./= (std(q, dims=2, corrected=false) .+ eps(eltype(q)))
end

setup_normalizer(z::Type{DiagonalZNormalizer}, q, y) = normalize(z, q), DiagonalZNormalizer(y, lastlength(q))

@propagate_inbounds function advance!(z::DiagonalZNormalizer{T}) where T

    if z.i == 0
        return z.i = 1
    end
    @boundscheck if z.i + z.n > length(z.x)
        return z.i += 1
    end
    z.bufi = 0

    # Remove old point
    x = z.x[!, z.i]
    @avx z.s .-= x
    @avx z.ss .-= x.^2

    # Add new point
    x = z.x[!, z.i+z.n]
    @avx z.s .+= x
    @avx z.ss .+= x.^2
    @avx z.μ .= z.s./z.n
    @avx z.σ .= sqrt.(z.ss./z.n .- z.μ.^2)
    z.i += 1
end

@inline @propagate_inbounds function getindex(z::AbstractNormalizer{<:Any, 2}, i::Union{Number, AbstractRange}, j)
    @boundscheck 1 <= j <= z.n || throw(BoundsError(z,j))
    @boundscheck 1 <= 1 <= size(z.x, 1) || throw(BoundsError(z,i))
    xj = j+z.i-1
    @boundscheck xj <= lastlength(z.x) || throw(BoundsError(z,j))
    z.x[i, xj]
end

@inline @propagate_inbounds function getindex(z::DiagonalZNormalizer{T}, ::typeof(!), i::Int, inorder = i == z.bufi + 1) where T
    j = inorder ? i : z.n
    xj = z.i + i - 1
    @avx for k = 1:size(z.x, 1)
        z.buffer[k, j] = (z.x[k, xj] - z.μ[k]) / (z.σ[k] + eps(T))
    end
    if inorder
        z.bufi = i
    end
    z.buffer[!, j]
end

@inline @propagate_inbounds function getindex(z::AbstractNormalizer{<:Any, 2}, ::typeof(!), i::AbstractRange)
    @boundscheck (i[1] == z.i && length(i) == z.n) || throw(ArgumentError("Normalizers can only be indexed by ranges corresponding to their current state. Got range $i but state was $(z.i) corresponding to range $(z.i):$(z.i+z.n-1)"))
    z
end

Base.Matrix(z::AbstractNormalizer{<:Any, 2}) = z.x[:,z.i:z.i+z.n-1]

Base.length(z::DiagonalZNormalizer) = size(z.x,1) * z.n
Base.size(z::AbstractNormalizer{<:Any, 2}) = (size(z.x,1), z.n)

# Norm normalizer ==================================================================================

abstract type AbstractNormNormalizer{T,N} <: AbstractNormalizer{T,N} end

mutable struct NormNormalizer{T,N} <: AbstractNormNormalizer{T,N}
    x::Array{T,N}
    n::Int
    σ::T
    ss::T
    i::Int
    buffer::Array{T,N}
    bufi::Int
end

mutable struct SqNormNormalizer{T,N} <: AbstractNormNormalizer{T,N}
    x::Array{T,N}
    n::Int
    σ::T
    ss::T
    i::Int
    buffer::Array{T,N}
    bufi::Int
end

function NormNormalizer(x::AbstractArray{T,N}, n) where {T,N}
    @assert length(x) >= n
    ss = zero(T)
    @inbounds for i in 1:n
        ss += sum(abs2, x[!, i])
    end
    σ = sqrt(ss)
    if N == 1
        buffer = similar(x, n)
    else
        buffer = similar(x, size(x,1), n)
    end
    NormNormalizer(x, n, σ, ss, 0, buffer, 0)
end

function SqNormNormalizer(x::AbstractArray{T,N}, n) where {T,N}
    @assert length(x) >= n
    ss = zero(T)
    @inbounds for i in 1:n
        ss += sum(abs2, x[!, i])
    end
    σ = ss
    if N == 1
        buffer = similar(x, n)
    else
        buffer = similar(x, size(x,1), n)
    end
    SqNormNormalizer(x, n, σ, ss, 0, buffer, 0)
end

function normalize(::Type{N}, q::AbstractArray) where N <: AbstractNormNormalizer
    power = N == NormNormalizer ? 1 : 2
    q ./ (norm(q)^power .+ eps(eltype(q)))
end

setup_normalizer(z::Type{N}, q, y) where N <: AbstractNormNormalizer = normalize(N, q), N(y, lastlength(q))

@propagate_inbounds function advance!(z::AbstractNormNormalizer{T}) where T

    if z.i == 0
        return z.i = 1
    end
    @boundscheck if z.i + z.n > length(z.x)
        return z.i += 1
    end
    z.bufi = 0

    # Remove old point
    x = z.x[!, z.i]
    z.ss -= sum(abs2, x)

    # Add new point
    x = z.x[!, z.i+z.n]
    z.ss += sum(abs2, x)
    z.σ = z isa SqNormNormalizer ? z.ss : sqrt(z.ss)
    z.i += 1
end

@inline @propagate_inbounds function getindex(z::AbstractNormNormalizer{T}, ::typeof(!), i::Int, inorder = i == z.bufi + 1) where T
    j = inorder ? i : z.n
    xj = z.i + i - 1
    σ = z.σ + eps(T)
    @avx for k = 1:size(z.x, 1)
        z.buffer[k, j] = z.x[k, xj] / σ
    end
    if inorder
        z.bufi = i
    end
    z.buffer[!, j]
end

LinearAlgebra.norm(z::NormNormalizer) = z.σ
LinearAlgebra.norm(z::SqNormNormalizer) = sqrt(z.σ)

Base.length(z::AbstractNormNormalizer{<:Any,2}) = size(z.x,1) * z.n
Base.length(z::AbstractNormNormalizer{<:Any,1}) = z.n
SlidingDistancesBase.lastlength(z::AbstractNormNormalizer) = z.n
actuallastlength(x::AbstractNormNormalizer) = lastlength(x.x)

# ============================================================================================

for T in [ZNormalizer, DiagonalZNormalizer, NormNormalizer, SqNormNormalizer]
    @eval @inline @propagate_inbounds function normalize(::Type{$T}, z::$T)
        if z.bufi == z.n
            return z.buffer
        end
        for i = z.bufi+1:z.n
            z[!, i, true] # This populates the buffer
        end
        return z.buffer
    end
end


text:
# traj_folder = joinpath(dirname(pathof(TrajectoryOptimization)),"..")
# urdf_folder = joinpath(traj_folder, "dynamics/urdf")
# urdf_doublependulum = joinpath(urdf_folder, "doublependulum.urdf")
#
# acrobot_model = Model(urdf_doublependulum,[0.;1.]) # underactuated, only control for second "elbow" joint

@with_kw struct Acrobot{T} <: AbstractModel
    l::SVector{2,T} = @SVector [1.0, 1.0]
    m::SVector{2,T} = @SVector [1.0, 1.0]
    J::SVector{2,SVector{3,T}} = @SVector [(@SVector ones(3)), (@SVector ones(3))]
end

function dynamics(model::Acrobot, x, u)
    g = 9.81
    m1,m2 = model.m
    l1,l2 = model.l
    J1,J2 = model.J
    θ1,    θ2    = x[1], x[2]
    θ1dot, θ2dot = x[3], x[4]
    s1,c1 = sincos(θ1)
    s2,c2 = sincos(θ2)
    c12 = cos(θ1 + θ2)

    # mass matrix
    m11 = m1*l1^2 + J1[3] + m2*(l1^2 + l2^2 + 2*l1*l2*c2) + J2[3]
    m12 = m2*(l2^2 + l1*l2*c2 + J2[3])
    m22 = l2^2*m2 + J2[3]
    M = @SMatrix [m11 m12; m12 m22]

    # bias term
    tmp = l1*l2*m2*s2
    b1 = -(2 * θ1dot * θ2dot + θ2dot^2)*tmp
    b2 = tmp * θ1dot^2
    B = @SVector [b1, b2]

    # friction
    c = 1.0
    C = @SVector [c*θ1dot, c*θ2dot]

    # gravity term
    g1 = ((m1 + m2)*l2*c1 + m2*l2*c12) * g
    g2 = m2*l2*c12*g
    G = @SVector [g1, g2]

    # equations of motion
    τ = @SVector [0, u[1]]
    θddot = M\(τ - B - G - C)
    return @SVector [θ1dot, θ2dot, θddot[1], θddot[2]]
end

Base.size(::Acrobot) = 4,1


text:
"""
    paretoshapeplot(values; kwargs...)
    paretoshapeplot!(values; kwargs...)

Plot shape parameters of fitted Pareto tail distributions for diagnosing convergence.

# Arguments

  - `values`: may be either a vector of Pareto shape parameters or a [`PSISResult`](@ref).

# Keywords

  - `showlines=false`: if `true`, plot horizontal lines indicating relevant Pareto shape
    thresholds are drawn. See [`PSISResult`](@ref) for explanation of thresholds.
  - `backend::Symbol`: backend to use for plotting, defaulting to `:Plots`, unless `:Makie`
    is available.

All remaining keywords are passed to the plotting backend.

See [`psis`](@ref), [`PSISResult`](@ref).

!!! note

    Plots.jl or a Makie.jl backend must be loaded to use these functions.

# Examples

```julia
using PSIS, Distributions, Plots
proposal = Normal()
target = TDist(7)
x = rand(proposal, 100, 1_000)
log_ratios = logpdf.(target, x) .- logpdf.(proposal, x)
result = psis(log_ratios)
```

Plot with Plots.jl.

```julia
using Plots
plot(result; showlines=true)
```

Plot with GLMakie.jl.

```julia
using GLMakie
plot(result; showlines=true)
```
"""
paretoshapeplot, paretoshapeplot!

"""
    plotting_backend!(backend::Symbol)

Set default plotting backend. Valid values are `:Plots` and `:Makie`.
"""
plotting_backend!(backend::Symbol) = Base.setindex!(PLOTTING_BACKEND, backend)

function paretoshapeplot(args...; backend::Symbol=PLOTTING_BACKEND[], kw...)
    return _paretoshapeplot(Val(backend), args...; kw...)
end
function paretoshapeplot!(args...; backend::Symbol=PLOTTING_BACKEND[], kw...)
    return _paretoshapeplot!(Val(backend), args...; kw...)
end