text:
using JSON
using DataFrames
using Query
using CSV

# extract json from 7z
function extract_json(data_path, filepaths_to_extract)
    cd(data_path)
    try
        if typeof(filepaths_to_extract) <: AbstractArray
            for fp in filepaths_to_extract
                Base.run(`7z x $fp -aos`)
            end
        else
            Base.run(`7z x $filepaths_to_extract -aos`)
        end
    catch e
    end
    cd("..")
    return true
end

# remove json files (use after processing to save disk space)
function remove_json(data_path)
    to_remove = joinpath.(data_path, [i for i in readdir(data_path) if occursin(".json", i)])
    for i in to_remove
        rm(i)
    end
    return true
end

# read json file from given path
function read_json_sample(path)
    data_sample = Dict()
    open(path, "r") do file
        data_sample = JSON.parse(file)
    end
    return data_sample
end

# test if all elements in a list are equal [FIND REFERENCE]
allequal(x) = all(y -> y == x[1], x)

function read_and_split_raw_data(jsonpath)
    metadata = read_json_sample(jsonpath)
    results = pop!(metadata)
    return metadata, results
end

function metadata_to_df(metadata)
    metadata_df_list = DataFrame[]
    for i in 1:length(metadata)
        tmp_df_list = DataFrame[]
        for j in 1:length(metadata[i])
            push!(tmp_df_list, DataFrame(metadata[i][j]))
        end
        tmp_df = vcat(tmp_df_list..., cols=:union)
        push!(metadata_df_list, tmp_df)
    end
    metadata_df = vcat(metadata_df_list..., cols=:union)
    return metadata_df
end

function strip_protocol(s)
    if startswith(s, "http://")
        return s[8:end]
    elseif startswith(s, "https://")
        return s[9:end]
    else
        return s
    end
end

function strip_network(s)
    if startswith(s, "www.")
        return s[5:end]
    elseif startswith(s, "de-de.")
        return s[7:end]
    elseif startswith(s, "de.")
        return s[4:end]
    else
        return s
    end
end

function results_to_df(results)
    results_df_list = DataFrame[]
    for result_list in results
        tmp_url_list = [strip_network(strip_protocol(i["sourceUrl"])) for i in result_list["result"]]
        tmp_medium_list = [try i["medium"] catch e "NA" end for i in result_list["result"]]
        tmp_domain_list = [split(i, "/")[1] for i in tmp_url_list]
        # skip faulty result lists
        if allequal(tmp_url_list)
            continue
        end
        tmp_df = DataFrame(sourceUrl = tmp_url_list, medium = tmp_medium_list, domain = tmp_domain_list)
        tmp_df[:, :result_hash] .= result_list["result_hash"]
        tmp_df[:, :rank] = 1:length(tmp_url_list)
        push!(results_df_list, tmp_df)
    end
    results_df = vcat(results_df_list..., cols=:union)
    return results_df
end

function reduce_dataset(df)
    df = df |> @filter(_.search_type == "search") |> DataFrame
    cols_to_delete = intersect(names(df), ["search_type", "plugin_id", "plugin_version"])
    select!(df, Not(cols_to_delete))
    return df
end

##### ----------------------------- #####
##### RAW DATA PROCESSING PROCEDURE #####
##### ----------------------------- #####

if !("processed" in readdir())
    mkdir("processed")
end
pattern = ".7z"
file_endings = Dict(
    "7z" => ".7z",
    "json" => ".json"
)
data_path = "data"
archive_list = [i for i in readdir(data_path) if occursin(file_endings["7z"], i)]
filenames_map = Dict(
    "Alice Weidel" => "aliceweidel",
    "CSU" => "csu",
    "Christian Lindner" => "christianlindner",
    "Cem Özdemir" => "cemoezdemir",
    "Dietmar Bartsch" => "dietmarbartsch",
    "SPD" => "spd",
    "FDP" => "fdp",
    "CDU" => "cdu",
    "Bündnis90/Die Grünen" => "gruene",
    "Alexander Gauland" => "alexandergauland",
    "Angela Merkel" => "angelamerkel",
    "Martin Schulz" => "martinschulz",
    "AfD" => "afd",
    "Katrin Göring-Eckardt" => "katringoeringeckardt",
    "Die Linke" => "linke",
    "Sahra Wagenknecht" => "sahrawagenknecht"
)
keywords = collect(keys(filenames_map))
for archive in archive_list
    # read and seperate metadata and search results
    extract_json(data_path, archive)
    jsonpath = joinpath.(data_path, [i for i in readdir(data_path) if occursin(file_endings["json"], i)])[1]
    metadata, results = read_and_split_raw_data(jsonpath);
    remove_json(data_path)

    # put everything into one dataframe
    metadata = metadata_to_df(metadata)
    results = results_to_df(results)
    day_df = innerjoin(results, metadata, on = :result_hash)
    day_df = reduce_dataset(day_df)

    # create datasets by keyword
    for k in keywords
        filename = filenames_map[k] * ".csv"
        keyword_df = day_df |> @filter(_.keyword == k) |> DataFrame
        if filename in readdir("processed")
            CSV.write(joinpath("processed", filename), keyword_df, append = true)
        else
            CSV.write(joinpath("processed", filename), keyword_df)
        end
    end
end

##### ----------- #####
##### OTHER STUFF #####
##### ----------- #####

# yet another data issue: &sa...&ved -> some urls with google specific paths seem broken
# fix broken urls (or don't break them in the first place)
# remove "sa=.*"r


text:
@testset "Wirtinger" begin
    w = Wirtinger(1+1im, 2+2im)
    @test wirtinger_primal(w) == 1+1im
    @test wirtinger_conjugate(w) == 2+2im
    @test w + w == Wirtinger(2+2im, 4+4im)

    @test w + One() == w + 1 == w + Thunk(()->1) == Wirtinger(2+1im, 2+2im)
    @test w * One() == One() * w == w
    @test w * 2 == 2 * w == Wirtinger(2 + 2im, 4 + 4im)

    # TODO: other + methods stack overflow
    @test_throws ErrorException w*w
    @test_throws ArgumentError extern(w)
    for x in w
        @test x === w
    end
    @test broadcastable(w) == w
    @test_throws MethodError conj(w)
end

ext: jl


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/bound.jl
# Computes interval bounds of various functions.
#############################################################################

###
### AFFINE FUNCTIONS
###
function lower_interval_bound(m::GlobalOptimizer, f::AffineFunctionIneq)
    fL = f.constant
    for (c, j) in f.terms
        xL = _lower_bound(FullVar(), m, j)
        xU = _upper_bound(FullVar(), m, j)
        fL += (c > 0.0) ? c*xL : c*xU
    end
    return fL
end

function interval_bound(m::GlobalOptimizer, f::Union{AffineFunctionEq,AffineFunctionIneq})
    fL = fU = f.constant
    for (c, j) in f.terms
        xL = _lower_bound(FullVar(), m, j)
        xU = _upper_bound(FullVar(), m, j)
        if c > 0.0
            fL += c*xL
            fU += c*xU
        else
            fL += c*xU
            fU += c*xL
        end
    end
    return fL, fU
end

###
### QUADRATIC FUNCTIONS
###

function lower_interval_bound(m::GlobalOptimizer, f::BufferedQuadraticIneq)
    fval = Interval{Float64}(f.func.constant)
    for t in f.func.affine_terms
        c = t.coefficient
        j = t.variable.value
        xL = _lower_bound(FullVar(), m, j)
        xU = _upper_bound(FullVar(), m, j)
        fval += c > 0.0 ? c*xL : c*xU
    end
    for t in f.func.quadratic_terms
        c = t.coefficient
        i = t.variable_1.value
        j = t.variable_2.value
        xL = _lower_bound(FullVar(), m, i)
        xU = _upper_bound(FullVar(), m, i)
        if i == j
            if c > 0.0
                fval += (0.0 < xL) ? 0.5*c*xL*xL : ((xU <= 0.0) ? 0.5*c*xU*xU : 0.0)
            else
                fval += (xL < xU) ? 0.5*c*xU*xU : 0.5*c*xL*xL
            end
        else
            yL = _lower_bound(FullVar(), m, j)
            yU = _upper_bound(FullVar(), m, j)
            fval += c*Interval{Float64}(xL, xU)*Interval{Float64}(yL, yU)
        end
    end
    return fval.lo
end
function interval_bound(m::GlobalOptimizer, f::BufferedQuadraticIneq)

    fval = Interval{Float64}(f.func.constant)
    for t in f.func.affine_terms
        c = t.coefficient
        j = t.variable.value
        xL = _lower_bound(FullVar(), m, j)
        xU = _upper_bound(FullVar(), m, j)
        fval += c*Interval(xL, xU)
    end
    for t in f.func.quadratic_terms
        c = t.coefficient
        i = t.variable_1.value
        j = t.variable_2.value
        xL = _lower_bound(FullVar(), m, i)
        xU = _upper_bound(FullVar(), m, i)
        if i == j
            fval += 0.5*c*pow(Interval(xL, xU), 2)
        else
            yL = _lower_bound(FullVar(), m, j)
            yU = _upper_bound(FullVar(), m, j)
            fval += c*Interval(xL, xU)*Interval(yL, yU)
        end
    end
    return fval.lo, fval.hi
end

function interval_bound(m::GlobalOptimizer, f::BufferedQuadraticEq)

    fval = Interval{Float64}(f.func.constant)
    for t in f.func.affine_terms
        c = t.coefficient
        j = t.variable.value
        xL = _lower_bound(FullVar(), m, j)
        xU = _upper_bound(FullVar(), m, j)
        fval += c*Interval(xL, xU)
    end
    for t in f.func.quadratic_terms
        c = t.coefficient
        i = t.variable_1.value
        j = t.variable_2.value
        xL = _lower_bound(FullVar(), m, i)
        xU = _upper_bound(FullVar(), m, i)
        if i == j
            fval += 0.5*c*pow(Interval(xL, xU), 2)
        else
            yL = _lower_bound(FullVar(), m, j)
            yU = _upper_bound(FullVar(), m, j)
            fval += c*Interval(xL, xU)*Interval(yL, yU)
        end
    end
    return fval.lo, fval.hi
end

###
### SECOND-ORDER CONE
###
function lower_interval_bound(m::GlobalOptimizer, d::BufferedSOC)

    sol_branch_map = m._sol_to_branch_map
    lo_bnds = n.lower_variable_bounds
    up_bnds = n.upper_variable_bounds
    vec_of_vi = d.variables.variables

    norm_bound = Interval(0.0)
    for i = 2:length(vec_of_vi)
        mapped_vi = @inbounds sol_branch_map[vec_of_vi[i].value]
        x = Interval{Float64}(lo_bnds[mapped_vi], up_bnds[mapped_vi])
        norm_bound += pow(x, 2)
    end
    norm_bound = sqrt(norm_bound)

    mapped_vi = @inbounds sol_branch_map[vec_of_vi[1].value]
    lower_bound = norm_bound.lo -(@inbounds up_bnds[mapped_vi])

    return lower_bound
end

###
### NONLINEAR FUNCTIONS
###
function lower_interval_bound(m::GlobalOptimizer, d::BufferedNonlinearFunction{V,N,T}) where {V,N,T}
    !has_value(d) && forward_pass!(m._working_problem._relaxed_evaluator, d)
    is_num(d) ? num(d) : interval(d).lo
end
function interval_bound(m::GlobalOptimizer, d::BufferedNonlinearFunction{V,N,T}) where {V,N,T}
    !has_value(d) && forward_pass!(m._working_problem._relaxed_evaluator, d)
    v = is_num(d) ? Interval{Float64}(num(d)) : interval(d)
    return v.lo, v.hi
end

is_feasible(m::GlobalOptimizer, f::Union{AFI,BQI}) = lower_interval_bound(m, f) <= 0.0
function is_feasible(m::GlobalOptimizer, f::Union{AFE,BQE})
    l, u = interval_bound(m, f)
    l <= 0.0 <= u
end
function is_feasible(m::GlobalOptimizer, f::BufferedNonlinearFunction{V,N,T}) where {V,N,T}
    l, u = interval_bound(m, f)
    feasible_flag = (u >= lower_bound(f))
    feasible_flag && (l <= upper_bound(f))
end

bound_objective(m::GlobalOptimizer, f::BufferedNonlinearFunction) = interval_bound(m, f)
bound_objective(m::GlobalOptimizer, f::AffineFunctionIneq) = interval_bound(m, f)
bound_objective(m::GlobalOptimizer, f::BufferedQuadraticIneq) = interval_bound(m, f)
function bound_objective(m::GlobalOptimizer, f::VI)
    vval = f.value
    l = lower_bound(FullVar(), m, vval)
    u = upper_bound(FullVar(), m, vval)
    return l, u
end

bound_objective(t::ExtensionType, m::GlobalOptimizer) = bound_objective(m, m._working_problem._objective)
bound_objective(m::GlobalOptimizer{R,Q,S}) where {R,Q,S<:ExtensionType} = bound_objective(_ext(m), m)

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
using Documenter, ChainedFixes

makedocs(;
    modules=[ChainedFixes],
    format=Documenter.HTML(),
    pages=[
        "ChainedFixes" => "index.md",
    ],
    sitename="ChainedFixes.jl",
    authors="Zachary P. Christensen",
)

deploydocs(
    repo = "github.com/Tokazama/ChainedFixes.jl.git",
)


text:
struct CYBERTRUCK{T} <: Model
    n::Int
    m::Int
    d::Int

    mass::T
    J::T
    μ::T
    g::T

    nq
    nu
    nc
    nf
    nb
	ns

    idx_u
    idx_λ
    idx_b
    idx_ψ
    idx_η
    idx_s
end

# Dimensions
nq = 4 # configuration dim
nu = 2 # control dim
nc = 1 # number of contact points
nf = 4 # number of faces for friction cone pyramid
nb = nc * nf
ns = 1

# Parameters
g = 9.81   # gravity
μ = 0.5    # coefficient of friction
mass = 1.0 # body mass
J = 1.0    # body inertia

n = 2 * nq
m = nu + nc + nb + nc + nb + ns
d = 0

idx_u = (1:nu)
idx_λ = nu .+ (1:nc)
idx_b = nu + nc .+ (1:nb)
idx_ψ = nu + nc + nb .+ (1:nc)
idx_η = nu + nc + nb + nc .+ (1:nb)
idx_s = nu + nc + nb + nc + nb .+ (1:ns)

# Methods
M_func(model::CYBERTRUCK, q) = Diagonal(@SVector [model.mass, model.mass, model.mass, model.J])
G_func(model::CYBERTRUCK, q) = @SVector [0.0, 0.0, model.mass * model.g, 0.0]

function ϕ_func(::CYBERTRUCK, q)
    return @SVector [q[3]]
end

N_func(::CYBERTRUCK, q) = @SMatrix [0.0 0.0 1.0 0.0]

P_func(::CYBERTRUCK, q) = @SMatrix [1.0 0.0 0.0 0.0;
                                     -1.0 0.0 0.0 0.0;
                                     0.0 1.0 0.0 0.0;
                                     0.0 -1.0 0.0 0.0]

B_func(::CYBERTRUCK, q) = @SMatrix [cos(q[4]) sin(q[4]) 0.0 0.0;
                                    0.0 0.0 0.0 1.0]

function friction_cone(model::CYBERTRUCK, u)
	λ = u[model.idx_λ]
	b = u[model.idx_b]
	return @SVector [model.μ * λ[1] - sum(b)]
end

function maximum_dissipation(model::CYBERTRUCK, x⁺, u, h)
	q3 = x⁺[model.nq .+ (1:model.nq)]
	q2 = x⁺[1:model.nq]
	ψ = u[model.idx_ψ]
	ψ_stack = ψ[1] * ones(4)
	η = u[model.idx_η]
	return P_func(model, q3) * (q3 - q2) / h + ψ_stack - η
end

function fd(model::CYBERTRUCK, x⁺, x, u, w, h, t)
	q3 = view(x⁺, model.nq .+ (1:model.nq))
	q2⁺ = view(x⁺, 1:model.nq)
	q2⁻ = view(x, model.nq .+ (1:model.nq))
	q1 = view(x, 1:model.nq)
	u_ctrl = view(u, model.idx_u)
	λ = view(u, model.idx_λ)
	b = view(u, model.idx_b)

    [q2⁺ - q2⁻;
    ((1.0 / h) * (M_func(model, q1) * (SVector{4}(q2⁺) - SVector{4}(q1))
    - M_func(model, q2⁺) * (SVector{4}(q3) - SVector{4}(q2⁺)))
    + transpose(B_func(model, q3)) * SVector{2}(u_ctrl)
    + transpose(N_func(model, q3)) * SVector{1}(λ)
    + transpose(P_func(model, q3)) * SVector{4}(b)
    - h * G_func(model, q2⁺))]
end

model = CYBERTRUCK(n, m, d,
				   mass, J, μ, g,
				   nq,
  				   nu,
  				   nc,
  				   nf,
  				   nb,
  				   ns,
  				   idx_u,
  				   idx_λ,
  				   idx_b,
  				   idx_ψ,
  				   idx_η,
  				   idx_s)

function visualize!(vis, model::CYBERTRUCK, q;
		Δt = 0.1)

    obj_path = joinpath(pwd(), "src/models/cybertruck/cybertruck.obj")
    mtl_path = joinpath(pwd(), "src/models/cybertruck/cybertruck.mtl")

    ctm = ModifiedMeshFileObject(obj_path, mtl_path, scale=0.1)
    setobject!(vis["cybertruck"], ctm)
    settransform!(vis["cybertruck"], LinearMap(RotZ(pi) * RotX(pi / 2.0)))

    anim = MeshCat.Animation(convert(Int,floor(1.0 / Δt)))

    for t = 1:length(q)
        MeshCat.atframe(anim, t) do
            settransform!(vis["cybertruck"], compose(Translation(q[t][1:3]...),
				LinearMap(RotZ(q[t][4] + pi) * RotY(0.0) * RotX(pi / 2.0))))
        end
    end
    # settransform!(vis["/Cameras/default"], compose(Translation(-1, -1, 0),LinearMap(RotZ(pi/2))))
    MeshCat.setanimation!(vis, anim)
end

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
using NLSolversObjectives
using Test

@testset "NLSolversObjectives.jl" begin
    # Write your own tests here.
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Print function of the html backend.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# Low-level function to print the table using the text backend.
function _pt_html(io::IO, pinfo::PrintInfo;
                  tf::HTMLTableFormat = tf_html_default,
                  highlighters::Union{HTMLHighlighter,Tuple} = (),
                  linebreaks::Bool = false,
                  noheader::Bool = false,
                  nosubheader::Bool = false,
                  sortkeys::Bool = false,
                  standalone::Bool = true)

    # Unpack fields of `pinfo`.
    data                    = pinfo.data
    header                  = pinfo.header
    id_cols                 = pinfo.id_cols
    id_rows                 = pinfo.id_rows
    num_rows                = pinfo.num_rows
    num_cols                = pinfo.num_cols
    num_printed_cols        = pinfo.num_printed_cols
    num_printed_rows        = pinfo.num_printed_rows
    header_num_rows         = pinfo.header_num_rows
    header_num_cols         = pinfo.header_num_cols
    show_row_number         = pinfo.show_row_number
    row_number_column_title = pinfo.row_number_column_title
    show_row_names          = pinfo.show_row_names
    row_names               = pinfo.row_names
    row_name_alignment      = pinfo.row_name_alignment
    row_name_column_title   = pinfo.row_name_column_title
    alignment               = pinfo.alignment
    cell_alignment          = pinfo.cell_alignment
    formatters              = pinfo.formatters
    compact_printing        = pinfo.compact_printing
    title                   = pinfo.title
    title_alignment         = pinfo.title_alignment
    header_alignment        = pinfo.header_alignment
    header_cell_alignment   = pinfo.header_cell_alignment
    cell_first_line_only    = pinfo.cell_first_line_only
    renderer                = pinfo.renderer

    # Unpack fields of `tf`.
    css         = tf.css
    table_width = tf.table_width

    # Let's create a `IOBuffer` to write everything and then transfer to `io`.
    buf_io = IOBuffer()
    buf    = IOContext(buf_io)

    !noheader && num_cols != header_num_cols &&
    error("The header length must be equal to the number of columns.")

    # Additional processing necessary if the user wants to print the header.
    if !noheader
        # If the user do not want to print the sub-header but wants to print the
        # header, then just force the number of rows in header to be 1.
        if nosubheader
            # Now, `header` will be a view of the first line of the matrix that
            # has the header.
            header = @view header[1:header_num_rows:end]
            header_num_rows = 1
        end
    end

    # Make sure that `highlighters` is always a tuple.
    !(highlighters isa Tuple) && (highlighters = (highlighters,))

    # Get the string which is printed when `print` is called in each element of
    # the matrix. Notice that we must create only the matrix with the printed
    # rows and columns.
    header_str = Matrix{String}(undef, header_num_rows, num_printed_cols)
    data_str   = Matrix{String}(undef, num_printed_rows, num_printed_cols)

    @inbounds for i = 1:num_printed_cols
        # Index of the i-th printed column in `data`.
        ic = id_cols[i]

        if !noheader
            for j = 1:header_num_rows
                header_str[j,i] =
                    _parse_cell_html(header[(ic-1)*header_num_rows + j],
                                     compact_printing = compact_printing,
                                     renderer = Val(:print))
            end
        end

        for j = 1:num_printed_rows
            # Index of the j-th printed row in `data`.
            jr = id_rows[j]

            # Apply the formatters.
            data_ij = isassigned(data,jr,ic) ? data[jr,ic] : undef

            for f in formatters
                data_ij = f(data_ij, jr, ic)
            end

            data_str[j,i] = _parse_cell_html(data_ij;
                                             cell_first_line_only = cell_first_line_only,
                                             compact_printing = compact_printing,
                                             linebreaks = linebreaks,
                                             renderer = renderer)
        end
    end

    # Variables to store information about indentation
    # ==========================================================================

    il = 0 # ......................................... Current indentation level
    ns = 2 # ........................ Number of spaces in each indentation level

    # Print HTML header
    # ==========================================================================

    if standalone
        _aprintln(buf, """
                  <!DOCTYPE html>

                  <meta charset=\"UTF-8\">

                  """)
    end

    _aprintln(buf, "
", il, ns)
    il += 1

    # Table title
    # ==========================================================================

    if length(title) > 0
        style = Dict{String,String}("text-align" => _html_alignment[title_alignment])
        _aprintln(buf, _styled_html("caption", title, style), il, ns)
    end

    # If there is no column or row to be printed, then just exit.
    if (num_printed_cols == 0) || (num_printed_rows == 0)
        @goto print_to_output
    end

    # Data header
    # ==========================================================================

    # Header and sub-header texts
    # --------------------------------------------------------------------------

    if !noheader
        @inbounds @views for i = 1:header_num_rows
            if (i == 1) && (header_num_rows == 1)
                _aprintln(buf, "
", il, ns)
            elseif i == 1
                _aprintln(buf, "
", il, ns)
            elseif i == header_num_rows
                _aprintln(buf, "
", il, ns)
            else
                _aprintln(buf, "
", il, ns)
            end
            il += 1

            # The text "Row" must appear only on the first line.
            if show_row_number
                if i == 1
                    _aprintln(buf, "" *
                              row_number_column_title *
                              "	", il, ns)
                else
                    _aprintln(buf, "	", il, ns)
                end
            end

            # The row name column title must appear only on the first line.
            if show_row_names
                if i == 1
                    style = Dict{String,String}("text-align" => _html_alignment[row_name_alignment])
                    row_name_title_html = _styled_html("th",
                                                       row_name_column_title,
                                                       style;
                                                       class = "rowName")
                    _aprintln(buf, row_name_title_html, il, ns)
                else
                    _aprintln(buf, "	", il, ns)
                end
            end

            for j = 1:num_printed_cols
                # Index of the j-th printed column in `data`.
                jc = id_cols[j]

                # Check the alignment of this cell.
                alignment_ij = header_alignment[jc]

                for f in header_cell_alignment
                    aux = f(header, i, jc)

                    if aux ∈ (:l, :c, :r, :L, :C, :R, :s, :S)
                        alignment_ij = aux
                        break
                    end
                end

                # If alignment is `:s`, then we must use the column alignment.
                alignment_ij ∈ (:s,:S) && (alignment_ij = alignment[jc])

                # Alignment of this cell.
                style = Dict{String,String}("text-align" => _html_alignment[alignment_ij])

                _aprintln(buf, _styled_html("th", header_str[i,j], style), il, ns)
            end
            il -= 1
            _aprintln(buf, "
", il, ns)
        end
    end

    # Data
    # ==========================================================================

    @inbounds @views for i = 1:num_printed_rows
        ir = id_rows[i]

        _aprintln(buf, "
", il, ns)
        il += 1

        if show_row_number
            _aprintln(buf, "" * string(ir) * "	", il, ns)
        end

        if show_row_names
            row_name_i_str = _parse_cell_html(row_names[i];
                                              cell_first_line_only = false,
                                              compact_printing = compact_printing,
                                              linebreaks = false,
                                              renderer = renderer)
            style = Dict{String,String}("text-align" => _html_alignment[row_name_alignment])
            row_name_i_html = _styled_html("td", row_name_i_str, style;
                                           class = "rowName")
            _aprintln(buf, row_name_i_html, il, ns)
        end

        for j = 1:num_printed_cols
            jc = id_cols[j]

            # Check the alignment of this cell.
            alignment_ij = alignment[jc]

            for f in cell_alignment
                aux = f(_getdata(data), ir, jc)

                if aux ∈ [:l, :c, :r, :L, :C, :R]
                    alignment_ij = aux
                    break
                end
            end

            # Alignment of this cell.
            style = Dict{String,String}("text-align" => _html_alignment[alignment_ij])

            # If we have highlighters defined, then we need to verify if this
            # data should be highlight.
            for h in highlighters
                if h.f(_getdata(data), ir, jc)
                    merge!(style, Dict(h.fd(h,_getdata(data),i,j)))
                    break
                end
            end

            _aprintln(buf, _styled_html("td", data_str[i,j], style), il, ns)
        end

        il -= 1
        _aprintln(buf, "
", il, ns)
    end

    @label print_to_output

    # Print HTML footer
    # ==========================================================================

    il -= 1
    _aprintln(buf, "
", il, ns)
    if standalone
        _aprintln(buf, """

                  """, il, ns)
    end

    # Print the buffer into the io.
    # ==========================================================================

    # If we are printing to `stdout`, then wrap the output in a `HTML` object.
    if (io == stdout) || ( (io isa IOContext) && (io.io == stdout) )
        display("text/html", HTML(String(take!(buf_io))))
    else
        print(io, String(take!(buf_io)))
    end

    return nothing
end


text:
# early exit
    # if (max(primfeas, dualfeas) < stoptol) | (iter == maxiter)
    #   breakyes = 1
    # end
    # if (max(primfeas, dualfeas) < sqrt(stoptol)) & (dualfeas > 1.5 * minimum(runhist.dualfeas[max(iter - 49, 1):iter])) & (iter > 150)
    #   breakyes = 2
    # end

# tune ρ
    # if (ρReset > 0) & (rem(iter, 10)==0)
    #   if (primfeas < 0.5 * dualfeas)
    #     ρ = 0.7 * ρ
    #   elseif (primfeas > 2 * dualfeas)
    #     ρ = 1.3 * ρ
    #   end
    # end
    # if (breakyes > 0)
    #   @printf("\n break = %1.0f\n", breakyes)
    #   break;
# end

# function sdpProjection(mat::Array{Float64, 2})
#   λ, X    = eigs(KrylovMethods(), mat, nev = 20)
#   # @show(λ)
#   # @show(X)
#   # @show(size(λ))
#   # @show(size(X))
#   return project(λ, X)
#   # posEigenValuesIndex = findall(x -> x > 0, λ);
#   # posEigenValues      = λ[posEigenValuesIndex];
#   # posEigenVectors     = X[:,posEigenValuesIndex];
#   # projectedMatrix     = posEigenVectors * diagm(0 => posEigenValues) *posEigenVectors';
#   # return projectedMatrix;
# end

# function sdpProjection0(data)
#   eigDecomposition    = eigen(data);
#   posEigenValuesIndex = findall(x -> x>0,eigDecomposition.values);
#   posEigenValues      = eigDecomposition.values[posEigenValuesIndex];
#   posEigenVectors     = eigDecomposition.vectors[:,posEigenValuesIndex];
#   projectedMatrix     = posEigenVectors * diagm(0 => posEigenValues) *posEigenVectors';
#   return projectedMatrix;
# end

# # function logisticLoss(x,y)
# #   f_x = Losses.σ.(x);
# #   return -sum(y .* log.(f_x) + (1 .- y) .* log.(1 .- f_x));
# # end

  # return X * diagm(0 => Λ) * X'

  # TODO: further extract positive eigen values


text:
using Game2048
using BenchmarkTools

@benchmark move(board, LEFT) setup=(board=rand(UInt64))

@benchmark move(board, left) setup=(board=rand(UInt64))
@benchmark move(board, right) setup=(board=rand(UInt64))
@benchmark move(board, up) setup=(board=rand(UInt64))
@benchmark move(board, down) setup=(board=rand(UInt64))


text:
using FunctionOperators, Test

@testset "Constructors - StructDefs.jl" begin
    @testset "Proper constructors" begin
        @testset "Keyword constructors" begin
            @test FunctionOperator{Float64}(name = "Op₁",
                forw = x -> x, backw = x -> x,
                inDims = (1,), outDims = (1,)) isa FunOp
            @test FunctionOperator{Float64}(
                forw = x -> x, backw = x -> x,
                inDims = (1,), outDims = (1,)) isa FunOp
            @test FunctionOperator{Float64}(name = "Op₁",
                forw = x -> x,
                inDims = (1,), outDims = (1,)) isa FunOp
            @test FunctionOperator{Float64}(forw = x -> x,
                inDims = (1,), outDims = (1,)) isa FunOp
            @test FunctionOperator{Float64}(forw = (buffer,x) -> x,
                inDims = (1,), outDims = (1,)) isa FunOp
            op = FunctionOperator{Float64}(forw = x -> x, inDims = (1,), outDims = (1,))
            @test op.getScale() == nothing
        end
        @testset "Positional constructors" begin
            @test FunctionOperator{Float64}("Op₁", x -> x, x -> x,
                (1,), (1,)) isa FunOp
            @test FunctionOperator{Float64}(x -> x, x -> x,
                (1,), (1,)) isa FunOp
            @test FunctionOperator{Float64}("Op₁", x -> x,
                (1,), (1,)) isa FunOp
            @test FunctionOperator{Float64}(x -> x,
                (1,), (1,)) isa FunOp
        end
    end
    @testset "Missing value" begin
        @test_throws ErrorException FunctionOperator{Float64}(name = "Op₁",
            backw = x -> x,
            inDims = (1,), outDims = (1,))
        @test_throws ErrorException FunctionOperator{Float64}(name = "Op₁",
            forw = x -> x, backw = x -> x,
            outDims = (1,))
        @test_throws ErrorException FunctionOperator{Float64}(name = "Op₁",
            forw = x -> x, backw = x -> x,
            inDims = (1,))
    end
    @testset "No arguments for forw" begin
        @testset "Keyword constructors" begin
            @test_throws AssertionError FunctionOperator{Float64}(name = "Op₁",
                forw = () -> x, backw = x -> x,
                inDims = (1,), outDims = (1,))
            @test_throws AssertionError FunctionOperator{Float64}(
                forw = () -> x, backw = x -> x,
                inDims = (1,), outDims = (1,))
            @test_throws AssertionError FunctionOperator{Float64}(name = "Op₁",
                forw = () -> x,
                inDims = (1,), outDims = (1,))
            @test_throws AssertionError FunctionOperator{Float64}(forw = () -> x,
                inDims = (1,), outDims = (1,))
        end
        @testset "Positional constructors" begin
            @test_throws AssertionError FunctionOperator{Float64}("Op₁", () -> x, x -> x,
                (1,), (1,))
            @test_throws AssertionError FunctionOperator{Float64}(() -> x, x -> x,
                (1,), (1,))
            @test_throws AssertionError FunctionOperator{Float64}("Op₁", () -> x,
                (1,), (1,))
            @test_throws AssertionError FunctionOperator{Float64}(() -> x,
                (1,), (1,))
        end
    end
    @testset "Too many arguments for forw" begin
        @testset "Keyword constructors" begin
            @test_throws AssertionError FunctionOperator{Float64}(name = "Op₁",
                forw = (x,y,z) -> x, backw = x -> x,
                inDims = (1,), outDims = (1,))
            @test_throws AssertionError FunctionOperator{Float64}(
                forw = (x,y,z) -> x, backw = x -> x,
                inDims = (1,), outDims = (1,))
            @test_throws AssertionError FunctionOperator{Float64}(name = "Op₁",
                forw = (x,y,z) -> x,
                inDims = (1,), outDims = (1,))
            @test_throws AssertionError FunctionOperator{Float64}(forw = (x,y,z) -> x,
                inDims = (1,), outDims = (1,))
        end
        @testset "Positional constructors" begin
            @test_throws AssertionError FunctionOperator{Float64}("Op₁", (x,y,z) -> x,
                x -> x, (1,), (1,))
            @test_throws AssertionError FunctionOperator{Float64}((x,y,z) -> x, x -> x,
                (1,), (1,))
            @test_throws AssertionError FunctionOperator{Float64}("Op₁", (x,y,z) -> x,
                (1,), (1,))
            @test_throws AssertionError FunctionOperator{Float64}((x,y,z) -> x,
                (1,), (1,))
        end
    end
    @testset "Forw has more arguments as backw" begin
        @testset "Keyword constructors" begin
            @test_throws AssertionError FunctionOperator{Float64}(name = "Op₁",
                forw = (x,y) -> x.^3, backw = x -> x,
                inDims = (1,), outDims = (1,))
            @test_throws AssertionError FunctionOperator{Float64}(
                forw = (x,y) -> x.^3, backw = x -> x,
                inDims = (1,), outDims = (1,))
        end
        @testset "Positional constructors" begin
            @test_throws AssertionError FunctionOperator{Float64}("Op₁", (x,y) -> x.^3,
                x -> x, (1,), (1,))
            @test_throws AssertionError FunctionOperator{Float64}((x,y) -> x.^3, x -> x,
                (1,), (1,))
        end
    end
    @testset "Unique default name" begin
        @test FunctionOperator{Float64}(x->x,(1,),(1,)) ≠
            FunctionOperator{Float64}(x->x,(1,),(1,))
    end
    @testset "Undefined backw" begin
        Op₁ = FunctionOperator{Float64}(x->x,(1,),(1,))
        @test_throws ErrorException Op₁' * [1.]
    end
end;


text:
# partial constraint development
using Test
using IncrementalInference
# going to introduce two new constraint mutable structs
import IncrementalInference: getSample
using Statistics

mutable struct DevelopPartial <: IncrementalInference.FunctorSingleton
  x::Distribution
  partial::Tuple # should rather be static types or templates for performance Tuple{Int, Int} etc.
end
getSample(dpl::DevelopPartial, N::Int=1) = (rand(dpl.x, N)', )

mutable struct DevelopDim2 <: IncrementalInference.FunctorSingleton
  x::Distribution
end
getSample(dpl::DevelopDim2, N::Int=1) = (rand(dpl.x, N), )

N=100 # 50
fg = initfg()

v1 = addVariable!(fg,:x1,ContinuousMultivariate(2,manifolds=(:Euclid, :Euclid)),N=N)

pr = DevelopDim2(MvNormal([0.0;0.0], 0.01*Matrix{Float64}(LinearAlgebra.I,2,2)))
f1  = addFactor!(fg,[:x1],pr)

dp = DevelopPartial(Normal(2.0, 1.0),(1,))
f2  = addFactor!(fg,[:x1],dp)

@testset "test evaluation of full constraint prior" begin
  pts = evalFactor2(fg, f1, v1.label, N=N)
  @test size(pts,1) == 2
  @test size(pts,2) == N
  @test norm(Statistics.mean(pts,dims=2)[1] .- [0.0]) < 0.3
end

memcheck = getVal(v1)

@testset "test evaluation of partial constraint prior" begin
  X1pts = getVal(v1)
  pts = evalFactor2(fg, f2, v1.label, N=N)

  @test size(pts, 1) == 2
  @test size(pts, 2) == N
  @test norm(Statistics.mean(pts,dims=2)[1] .- [2.0]) < 0.75
  # ensure the correct response from
  @test norm(X1pts[1,:] - pts[1,:]) > 2.0
  @test norm(X1pts[2,:] - pts[2,:]) < 1e-10
  @test norm(X1pts - memcheck) < 1e-10
end

@testset "test solving of factor graph" begin
  tree = wipeBuildNewTree!(fg)

  inferOverTree!(fg,tree, N=N)
  pts = getVal(fg, :x1)
  @test norm(Statistics.mean(pts,dims=2)[1] .- [0.0]) < 0.25
  @test norm(Statistics.mean(pts,dims=2)[2] .- [0.0]) < 0.25
end
# plotKDE(getVertKDE(fg, :x1),levels=3)

mutable struct DevelopPartialPairwise <: IncrementalInference.FunctorPairwise
  x::Distribution
  partial::Tuple
  DevelopPartialPairwise(x::Distribution) = new(x, (2,))
end
getSample(dpl::DevelopPartialPairwise, N::Int=1) = (rand(dpl.x, N)', )

function (dp::DevelopPartialPairwise)(res::Vector{Float64},
                              userdata::FactorMetadata,
                              idx::Int,
                              meas::Tuple, #{RowVector{Float64,Array{Float64,1}}}, #Tuple{Array{Float64,2}},
                              x1::Array{Float64},
                              x2::Array{Float64}  )
  #
  res[1] = meas[1][1,idx] - (x2[2,idx]-x1[2,idx])
  nothing
end

v2 = addVariable!(fg,:x2,ContinuousMultivariate(2),N=N)

dpp = DevelopPartialPairwise(Normal(10.0, 1.0))
f3  = addFactor!(fg,[:x1;:x2],dpp)

dp2 = DevelopPartial( Normal(-20.0, 1.0), (1,) )
f4  = addFactor!(fg,[:x2;], dp2)

@testset "test evaluation of multiple simultaneous partial constraints" begin
global fg

ensureAllInitialized!(fg)
valx2 = getVal(fg, :x2)
pts = approxConv(fg, :x1x2f1, :x2, N=N) # evalFactor2(fg, f3, v2.index, N=N)
@test size(pts,1) == 2
@test norm(Statistics.mean(pts,dims=2)[2] .- [10.0]) < 3.0
@test norm(valx2[1,:] - pts[1,:]) < 1e-5

pts = approxConv(fg, :x2f1, :x2, N=N) # evalFactor2(fg, f4, v2.index, N=N)
@test size(pts,1) == 2
@test norm(Statistics.mean(pts,dims=2)[1] .- [-20.0]) < 0.75
@test (Statistics.std(pts,dims=2)[1] .- 1.0) < 0.4

end

# keep previous values to ensure funciton evaluation is modifying correct data fields

@warn "restore findRelatedFromPotential as testset!"
# @testset "test findRelatedFromPotential..." begin
# global v2, fg, f3, f4, N

thefac = getFactor(fg, :x1x2f1)

X2lpts = getVal(getVariable(fg, :x2))
keepaside, = findRelatedFromPotential(fg, thefac, :x2, N)
@test Ndim(keepaside) == 2
lpts = KernelDensityEstimate.getPoints(keepaside)
@test size(lpts,2) == N

@show X2lpts[2,95:100]
@show lpts[2,95:100]
@show getPoints(keepaside)

# DevelopPartialPairwise must only modify the second dimension of proposal distribution on X2
@test norm(X2lpts[1,:] - lpts[1,:]) < 1e-10
# @test norm(X2lpts[2,:] - lpts[2,:]) > 1e-10 # 10*N*0.5 # why so big?
memcheck = getVal(v2)
@test norm(X2lpts - memcheck) < 1e-10

X2lpts = getVal(v2)
p4, = findRelatedFromPotential(fg, f4, v2.label, N)
@test Ndim(p4) == 2
lpts = KernelDensityEstimate.getPoints(keepaside)
@test size(lpts,2) == N

# DevelopPartialPairwise must only modify the second dimension of proposal distribution on X2
@test norm(X2lpts[1,:] - lpts[1,:]) < 1e-10
@test norm(X2lpts[2,:] - lpts[2,:]) > 1e-10 # 10*N*0.5 # why so big?
memcheck = getVal(v2)
@test norm(X2lpts - memcheck) < 1e-10

# end

@testset "test belief prediction with partials..." begin

global v2, fg

# partial prior
X2pts = getVal(v2)
val, = predictbelief(fg, v2, [f4], N=N)
@test norm(X2pts[2,:] - val[2,:]) < 1e-10
@test 0.0 < norm(X2pts[1,:] - val[1,:])
@test norm(Statistics.mean(val[1,:]) .+ 20.0) < 0.75

# partial pairwise
X2pts = getVal(v2)
val, = predictbelief(fg, v2, [f3], N=N)
@test norm(X2pts[1,:] - val[1,:]) < 1e-10
@test 0.0 < norm(X2pts[2,:] - val[2,:])
@test abs(Statistics.mean(val[2,:] - getVal(v1)[2,:]) .- 10.0) < 0.75

# combination of partials
val, = predictbelief(fg, v2, [f3;f4], N=N)
# plotKDE(kde!(val),levels=3)
@test norm(Statistics.mean(val,dims=2)[1] .- [-20.0]) < 2.0
@test norm(Statistics.mean(val,dims=2)[2] .- [10.0]) < 2.0
@test (Statistics.std(val,dims=2)[1] .- 1.0) < 3.0
@test (Statistics.std(val,dims=2)[2] .- 1.0) < 3.0

tree = wipeBuildNewTree!(fg )#, drawpdf=true)
# run(`evince bt.pdf`)

inferOverTree!(fg,tree, N=N)

pts = getVal(fg, :x1)
@test norm(Statistics.mean(pts,dims=2)[1] .- [0.0]) < 0.5
@test norm(Statistics.mean(pts,dims=2)[2] .- [0.0]) < 0.5

pts = getVal(fg, :x2)
@test norm(Statistics.mean(pts,dims=2)[1] .- [-20.0]) < 2.0
@test norm(Statistics.mean(pts,dims=2)[2] .- [10.0]) < 2.0
@test (Statistics.std(pts,dims=2)[1]-1.0) < 3.0
@test (Statistics.std(pts,dims=2)[2]-1.0) < 3.0

end

# plotKDE(getVertKDE(fg, :x2),levels=3)

# spyCliqMat(tree, :x2)

#


text:
struct RocketInfo1{T,R,RZ,Rθ}
	ip_dyn::InteriorPoint{T,R,RZ,Rθ}
	ip_proj::InteriorPoint{T,R,RZ,Rθ}
	idx_z::IndicesZ
	idx_θ::Indicesθ
	idx3::Vector{Int}
	h::T
	u_max::T
	du_dyn_cache::Matrix{T}
	du_proj_cache::Matrix{T}
end

function RocketInfo(rocket, u_max, h, r_func, rz_func, rθ_func, r_func_proj, rz_func_proj, rθ_func_proj)
	nx = rocket.nq
	nu = rocket.nu
	nw = rocket.nw

	nz = nx
	nθ = nx + nu + 1

	solver_opts=InteriorPointOptions(
				r_tol=1e-8,
				κ_tol=1.0,
				max_ls=25,
				ϵ_min=0.25,
				diff_sol=true,
				verbose=false)

	z0 = zeros(nx)
	θ0 = zeros(nθ)

	idx_opt = RoboDojo.indices_optimization(rocket)

	ip_dyn = interior_point(
		z0,
		θ0,
		idx=idx_opt,
		r! = r_func,
		rz! = rz_func,
		rθ! = rθ_func,
		rz=zeros(nz, nz),
		rθ=zeros(nz, nθ),
		opts=solver_opts)

	idx_z = RoboDojo.indices_z(rocket)
	idx_θ = RoboDojo.indices_θ(rocket)

	# projection
	nz_proj = 10
	nθ_proj = 4

	idx_opt_proj = IndicesOptimization(
						10,
						10,
						[[5, 3], [6, 4]],
						[[5, 3], [6, 4]],
						[[collect([3, 1, 2]), collect([10, 8, 9])]],
						[[collect([3, 1, 2]), collect([10, 8, 9])]],
						collect(1:5),
						collect(6:7),
						collect(8:10),
						[collect(8:10)],
						collect(6:10))

		z0_proj = zeros(nz_proj)
		θ0_proj = zeros(nθ_proj)

	ip_proj = interior_point(
		z0_proj,
		θ0_proj,
		idx=idx_opt_proj,
		r! = r_func_proj,
		rz! = rz_func_proj,
		rθ! = rθ_func_proj,
		rz=zeros(nz_proj, nz_proj),
		rθ=zeros(nz_proj, nθ_proj),
		opts=InteriorPointOptions(
			r_tol=1e-8,
			κ_tol=1.0e-4,
			max_ls=25,
			ϵ_min=0.0,
			undercut=Inf,
			γ_reg=0.0,
			κ_reg=0.0,
			diff_sol=false,
			verbose=false))

	RocketInfo1(
		ip_dyn,
		ip_proj,
		idx_z,
		idx_θ,
		collect(1:3),
		h,
		u_max,
		zeros(nx, nu),
		zeros(nu, nu)
	)
end

function f_rocket(d, info::RocketInfo1, x, u, w)
    # initialize
    info.ip_dyn.z[info.idx_z.q] .= x
    info.ip_dyn.θ[info.idx_θ.q1] .= x
    info.ip_dyn.θ[info.idx_θ.u] .= u
    info.ip_dyn.θ[info.idx_θ.h] .= info.h
    # solve
    info.ip_dyn.opts.diff_sol = false
    interior_point_solve!(info.ip_dyn)
    # solution
	sol = @views info.ip_dyn.z[info.idx_z.q]
    d .= sol
    return d
end

# d = zeros(nx)
# dx = zeros(nx, nx)
# du = zeros(nx, nu)

# x = rand(nx)
# u = rand(nu)
# w = rand(nw)

# using BenchmarkTools
# f_rocket(d, info, x, u, w)
# @benchmark f_rocket($d, $info, $x, $u, $w)

# fx_rocket(dx, info, x, u, w)
# @benchmark fx_rocket($dx, $info, $x, $u, $w)

# fu_rocket(du, info, x, u, w)
# @benchmark fu_rocket($du, $info, $x, $u, $w)

function fx_rocket(dx, info::RocketInfo1, x, u, w)
    # initialize
    info.ip_dyn.z[info.idx_z.q] .= x
    info.ip_dyn.θ[info.idx_θ.q1] .= x
    info.ip_dyn.θ[info.idx_θ.u] .= u
    info.ip_dyn.θ[info.idx_θ.h] .= info.h
    # solve
    info.ip_dyn.opts.diff_sol = true
    interior_point_solve!(info.ip_dyn)
    # solution
	sol = @views info.ip_dyn.δz[info.idx_z.q, info.idx_θ.q1]
    dx .= sol
    return dx
end

function fu_rocket(du, info::RocketInfo1, x, u, w)
    # initialize
    info.ip_dyn.z[info.idx_z.q] .= x
    info.ip_dyn.θ[info.idx_θ.q1] .= x
    info.ip_dyn.θ[info.idx_θ.u] .= u
    info.ip_dyn.θ[info.idx_θ.h] .= info.h
    # solve
    info.ip_dyn.opts.diff_sol = true
    interior_point_solve!(info.ip_dyn)

    # solution
	sol = @views info.ip_dyn.δz[info.idx_z.q, info.idx_θ.u]
    du .= sol
    return du
end

# info = RocketInfo(rocket, eval(r_func), eval(rz_func), eval(rθ_func), eval(r_func_proj), eval(rz_func_proj), eval(rθ_func_proj))

function soc_projection(x, info::RocketInfo1)
	info.ip_proj.z .= 0.1
    info.ip_proj.z[3] += 1.0
    info.ip_proj.z[10] += 1.0
    info.ip_proj.z[7] = 0.0

	info.ip_proj.θ[info.idx3] .= x
	info.ip_proj.θ[4] = info.u_max

	info.ip_proj.opts.diff_sol = false
	status = interior_point_solve!(info.ip_proj)

    # !status && (@warn "projection failure (res norm: $(norm(info.ip_proj.r, Inf))) \n
	# 	               z = $(info.ip_proj.z), \n
	# 				   θ = $(info.ip_proj.θ)")

	sol = @views info.ip_proj.z[info.idx3]
	return sol
end
# soc_projection(u, info)
# @benchmark soc_projection($u, $info)

function soc_projection_gradient(x, info::RocketInfo1)
	info.ip_proj.z .= 0.1
    info.ip_proj.z[3] += 1.0
    info.ip_proj.z[10] += 1.0
    info.ip_proj.z[7] = 0.0

	# info.ip_proj.θ .= [x; uu]
	info.ip_proj.θ[info.idx3] .= x
	info.ip_proj.θ[4] = info.u_max

	info.ip_proj.opts.diff_sol = true
	status = interior_point_solve!(info.ip_proj)

    # !status && (@warn "projection failure (res norm: $(norm(info.ip_proj.r, Inf))) \n
	# 	               z = $(info.ip_proj.z), \n
	# 				   θ = $(info.ip_proj.θ)")

	sol = @views info.ip_proj.δz[info.idx3, info.idx3]

	return sol
end

# soc_projection_gradient(u, info)
# @benchmark soc_projection_gradient($u, $info)

function f_rocket_proj(d, info::RocketInfo1, x, u, w)
    # initialize
    info.ip_dyn.z[info.idx_z.q] .= x
    info.ip_dyn.θ[info.idx_θ.q1] .= x
    info.ip_dyn.θ[info.idx_θ.u] .= soc_projection(u, info)
    info.ip_dyn.θ[info.idx_θ.h] .= info.h
    # solve
    info.ip_dyn.opts.diff_sol = false
    interior_point_solve!(info.ip_dyn)
    # solution
    sol = @views info.ip_dyn.z[info.idx_z.q]
	d .= sol
    return d
end

# f_rocket_proj(d, info, x, u, w)
# @benchmark f_rocket_proj($d, $info, $x, $u, $w)

# fx_rocket_proj(dx, info, x, u, w)
# @benchmark fx_rocket_proj($dx, $info, $x, $u, $w)

# fu_rocket_proj(du, info, x, u, w)
# @benchmark fu_rocket_proj($du, $info, $x, $u, $w)

function fx_rocket_proj(dx, info::RocketInfo1, x, u, w)
    # initialize
    info.ip_dyn.z[info.idx_z.q] .= x
    info.ip_dyn.θ[info.idx_θ.q1] .= x
    info.ip_dyn.θ[info.idx_θ.u] .= soc_projection(u, info)
    info.ip_dyn.θ[info.idx_θ.h] .= info.h
    # solve
    info.ip_dyn.opts.diff_sol = true
    interior_point_solve!(info.ip_dyn)
    # solution
    sol = @views info.ip_dyn.δz[info.idx_z.q, info.idx_θ.q1]
    dx .= sol
    return dx
end

function fu_rocket_proj(du, info::RocketInfo1, x, u, w)
    # initialize
    info.ip_dyn.z[info.idx_z.q] .= x
    info.ip_dyn.θ[info.idx_θ.q1] .= x
    info.ip_dyn.θ[info.idx_θ.u] .= soc_projection(u, info)
    info.ip_dyn.θ[info.idx_θ.h] .= info.h
    # solve
    info.ip_dyn.opts.diff_sol = true
    interior_point_solve!(info.ip_dyn)
    # solution
    info.du_dyn_cache .= @views info.ip_dyn.δz[info.idx_z.q, info.idx_θ.u]
	info.du_proj_cache .= soc_projection_gradient(u, info)
    # du .= sol_dyn * sol_proj
	mul!(du, info.du_dyn_cache, info.du_proj_cache)
    return du
end


text:
# assert file to test polynomial implementation
using Test
using LinearAlgebra
using Polynomials
using SpecialFunctions
using RecipesBase: apply_recipe

import SparseArrays: sparse, nnz
using OffsetArrays

@testset "Standard basis" begin include("StandardBasis.jl") end
@testset "ChebyshevT" begin include("ChebyshevT.jl") end
@testset "Poly, Pade (compatability)" begin include("Poly.jl") end


text:
#console_app.jl: Simple command line interface for AudioStick package.

#==Constants
===============================================================================#
const APPNAME = "AudioStick"

#Entries in config file (currently not supported):
const CFGID_FSPLROOT = "FILESYSTEMPLAYLIST_ROOTFOLDER"

const DFLT_TGTROOT = "Z:\\"

#==Types
===============================================================================#
mutable struct ASOptions
	srcpath::String #Input playlist file
	tgtroot::String #Target path (excluding playlist subfolder)
end
ASOptions(srcpath) = ASOptions(srcpath, DFLT_TGTROOT)

#==Helper functions
===============================================================================#

#Parse & validate command-line arguments:
function parseargs(args)
	if length(args) < 1
		@error("Must set source playlist in `ARGS`")
		error("Insufficient arguments.") #throws
	end

	options = ASOptions(joinpath(pwd(), strip(args[1])))
	if !isfile(options.srcpath)
		error("Playlist not found: $(options.srcpath)")
	end

	if length(args) >= 2
		options.tgtroot = strip(args[2])
	end

	return options
end

#==Main "console" entry point
===============================================================================#

function run_console(args)
	options = parseargs(args)
	srcfilename = basename(options.srcpath)

	println_ul(APPNAME)
	println("\nSource playlist: $srcfilename\n\n")

	#Ask for target path:
	prompt = "Target path (excluding playlist subfolder)"
	options.tgtroot = input(String, prompt, options.tgtroot)
	println()

	destfoldername = splitext(srcfilename)
	if length(destfoldername) < 1
		error("Playlist does not have a name.")
	end
	destfoldername = destfoldername[1]

	#Create filesystem-based playlist:
	destpath = joinpath(options.tgtroot, destfoldername)
	playlist = read_m3u(options.srcpath)
	synchronize(playlist, destpath)
end

"""`run_console()`

Create/synchronize a filesystem-based media playlist from the playlist *file*
specified in `ARGS`.
"""
run_console() = run_console(ARGS)

#Last line

max_stars_repo_name: ma-laforge/AudioStick

max_issues_repo_name: ma-laforge/AudioStick

max_forks_repo_name: ma-laforge/AudioStick


text:
@render_program LineProgram
@vao LineVao

LineUploader() = Uploader{LineProgram}()
shaders(::Type{LineProgram}) = line_shaders()

function set_entity_uniforms_func(render_program::LineProgram, system::System)
    prog     = render_program.program
    comp(T)  = component(system, T)
    modelmat = comp(ModelMat)
    line     = comp(Line)
    return e -> begin
        gluniform(prog, :modelmat, modelmat[e].modelmat)
        gluniform(prog, :thickness, line[e].thickness)
        gluniform(prog, :MiterLimit, line[e].miter)
    end
end

struct LineRenderer <: AbstractRenderSystem end

function Overseer.requested_components(::LineRenderer)
    return (LineVao, LineProgram, ModelMat, Material, PointLight, Spatial, Camera3D,
            IOTarget, LineOptions)
end

function Overseer.update(::LineRenderer, m::AbstractLedger)
    fbo  = singleton(m, IOTarget)
    prog = singleton(m, LineProgram)
    idc  = m[Selectable]
    bind(fbo)
    draw(fbo)
    glDisable(GL_BLEND)
    glEnable(GL_DEPTH_TEST)
    glDepthFunc(GL_LEQUAL)

    bind(prog)
    light, ucolor, spat, modelmat, cam = m[PointLight], m[UniformColor], m[Spatial],
                                         m[ModelMat], m[Camera3D]
    for e in @entities_in(light && ucolor && spat)
        gluniform(prog, light[e], ucolor[e], spat[e])
    end
    for e in @entities_in(spat && cam)
        gluniform(prog, spat[e], cam[e])
    end
    gluniform(prog, :Viewport, Vec2f0(size(m[IOTarget][1])))
    vao, modelmat, line, vis = m[LineVao], m[ModelMat], m[LineOptions], m[Visible]
    for e in @entities_in(vao && modelmat && line)
        evao = vao[e]
        e_line = line[e]
        if vis[e].visible
            gluniform(prog, :modelmat, modelmat[e].modelmat)
            gluniform(prog, :thickness, e_line.thickness)
            gluniform(prog, :MiterLimit, e_line.miter)
            if e in idc
                gluniform(prog, :object_id_color, idc[e].color)
            else
                gluniform(prog, :object_id_color, RGBf0(1.0, 1.0, 1.0))
            end
            GLA.bind(evao)
            GLA.draw(evao)
        end
    end
end


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# Attributes #
# ========================== #
"""
    AbstractQuasiGradientAttribute

Abstract supertype for attribute objects specific to the L-shaped algorithm.
"""
abstract type AbstractQuasiGradientAttribute <: AbstractStructuredOptimizerAttribute end
"""
    SubProblems

An optimizer attribute for specifying if subproblems should be smoothed. Options are:

- [`Unaltered`](@ref):  Subproblems are solved in their original state (default)
- [`Smoothed`](@ref):  Subproblems are smoothed using Moreau envelopes ?SmoothSubProblem for parameter descriptions.
"""
struct SubProblems <: AbstractQuasiGradientAttribute end
"""
    Prox

An optimizer attribute for specifying a prox policy to be used in the quasi-gradient algorithm. Options are:

- [`NoProx`](@ref): Unconstrained.
- [`Polyhedron`](@ref): QP projection on polyhedral space ?PolyhedronProjection for parameter descriptions. (default)
- [`AndersonAcceleration`](@ref):  Anderson acceleration of inner prox step ?AndersonAcceleratedProximal for parameter descriptions.
- [`Nesterov`](@ref):  Nesterov acceleration of inner prox step ?NesterovProximal for parameter descriptions.
- [`DryFriction`](@ref):  Dry-friction acceleration of inner prox step ?DryFrictionProximal for parameter descriptions.
"""
struct Prox <: AbstractQuasiGradientAttribute end
"""
    StepSize

An optimizer attribute for specifying an aggregation procedure to be used in the quasi-gradient algorithm. Options are:

- [`Constant`](@ref):  Constant step size ?ConstantStep for parameter descriptions (default)
- [`Diminishing`](@ref):  Diminishing step size ?DiminishingStep for parameter descriptions.
- [`Polyak`](@ref):  Polyak step size ?PolyakStep for parameter descriptions.
- [`BB`](@ref):  Barzilai-Borwein step size ?BBStep for parameter descriptions.
"""
struct StepSize <: AbstractQuasiGradientAttribute end
"""
    Termination

An optimizer attribute for specifying a termination criterion to be used in the quasi-gradient algorithm. Options are:

- [`AfterMaximumIterations`](@ref):  Terminate after set number of iterations ?MaximumIterations for parameter descriptions (default)
- [`AtObjectiveThreshold`](@ref):  Terminate after reaching reference objective ?ObjectiveThreshold for parameter descriptions.
- [`AtGradientThreshold`](@ref):  Terminate after reaching zero gradient ?GradientThreshold for parameter descriptions.
"""
struct Termination <: AbstractQuasiGradientAttribute end
"""
    ProxParameter

Abstract supertype for prox-specific attributes.
"""
abstract type ProxParameter <: AbstractQuasiGradientAttribute end
"""
    StepParameter

Abstract supertype for step-specific attributes.
"""
abstract type StepParameter <: AbstractQuasiGradientAttribute end
"""
    TerminationParameter

Abstract supertype for termination-specific attributes.
"""
abstract type TerminationParameter <: AbstractQuasiGradientAttribute end


text:
__precompile__()

module GridOperations

# export interp_missing,
export  interp_missing!, create_mask, interp2d,
        colorize,
        hist_filter, pmf_filter, vosselman_filter,
        im_thresholding, im_segments_cleanup, im_segments_reduce

using Colors
using ColorSchemes
using Compat
using GeoStats
using ImageFiltering
using Images
using NearestNeighbors
using ProgressMeter

include("utils.jl")
include("interpolation.jl")
include("segments.jl")
include("filter.jl")
include("color.jl")

end


text:
"""
    char(numbers...)

Convert a list of numbers to a string with the corresponding characters

```jldoctest
julia> using QuerySQLite

julia> char(65, 90)
"AZ"
```
"""
char(numbers...) = string(map(Char, numbers)...)
export char

"""
    if_else(switch, yes, no)

`ifelse` that you can add methods to.

```jldoctest
julia> using QuerySQLite

julia> if_else(true, 1, 0)
1

julia> if_else(false, 1, 0)
0
```
"""
function if_else(switch, yes, no)
    ifelse(switch, yes, no)
end
export if_else

"""
    instr(haystack, needle)

Find the first index of `needle` in `haystack`.

```jldoctest
julia> using QuerySQLite

julia> instr("QuerySQLite", "SQL")
6
```
"""
function instr(haystack, needle)
    first(findfirst(needle, haystack))
end
export instr

"""
    type_of(it)

`typeof` that you can add methods to.

```jldoctest
julia> using QuerySQLite

julia> type_of('a')
Char
```
"""
function type_of(it)
    typeof(it)
end
export type_of

"""
    hex(it)

Uppercase hexadecimal representation

```jldoctest
julia> using QuerySQLite

julia> hex("hello")
"68656C6C6F"
```
"""
function hex(it::Number)
    uppercase(string(it, base=16))
end
function hex(it::AbstractString)
    join(hex(byte) for byte in codeunits(it))
end
export hex


text:
function squeezenet_layers()
  weight = Metalhead.weights("squeezenet.bson")
  weights = Dict{Any ,Any}()
  for ele in keys(weight)
    weights[string(ele)] = weight[ele]
  end
    c_1 = Conv(weights["conv10_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["conv10_b_0"], stride=(1, 1), pad=(0, 0), dilation = (1, 1))
    c_2 = Dropout(0.5f0)
    c_3 = Conv(weights["fire9/expand1x1_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire9/expand1x1_b_0"], stride=(1, 1), pad=(0, 0), dilation = (1, 1))
    c_4 = Conv(weights["fire9/squeeze1x1_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire9/squeeze1x1_b_0"], stride=(1, 1), pad=(0, 0), dilation = (1, 1))
    c_5 = Conv(weights["fire8/expand1x1_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire8/expand1x1_b_0"], stride=(1, 1), pad=(0, 0), dilation = (1, 1))
    c_6 = Conv(weights["fire8/squeeze1x1_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire8/squeeze1x1_b_0"], stride=(1, 1), pad=(0, 0), dilation = (1, 1))
    c_7 = Conv(weights["fire7/expand1x1_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire7/expand1x1_b_0"], stride=(1, 1), pad=(0, 0), dilation = (1, 1))
    c_8 = Conv(weights["fire7/squeeze1x1_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire7/squeeze1x1_b_0"], stride=(1, 1), pad=(0, 0), dilation = (1, 1))
    c_9 = Conv(weights["fire6/expand1x1_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire6/expand1x1_b_0"], stride=(1, 1), pad=(0, 0), dilation = (1, 1))
    c_10 = Conv(weights["fire6/squeeze1x1_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire6/squeeze1x1_b_0"], stride=(1, 1), pad=(0, 0), dilation = (1, 1))
    c_11 = Conv(weights["fire5/expand1x1_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire5/expand1x1_b_0"], stride=(1, 1), pad=(0, 0), dilation = (1, 1))
    c_12 = Conv(weights["fire5/squeeze1x1_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire5/squeeze1x1_b_0"], stride=(1, 1), pad=(0, 0), dilation = (1, 1))
    c_13 = Conv(weights["fire4/expand1x1_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire4/expand1x1_b_0"], stride=(1, 1), pad=(0, 0), dilation = (1, 1))
    c_14 = Conv(weights["fire4/squeeze1x1_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire4/squeeze1x1_b_0"], stride=(1, 1), pad=(0, 0), dilation = (1, 1))
    c_15 = Conv(weights["fire3/expand1x1_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire3/expand1x1_b_0"], stride=(1, 1), pad=(0, 0), dilation = (1, 1))
    c_16 = Conv(weights["fire3/squeeze1x1_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire3/squeeze1x1_b_0"], stride=(1, 1), pad=(0, 0), dilation = (1, 1))
    c_17 = Conv(weights["fire2/expand1x1_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire2/expand1x1_b_0"], stride=(1, 1), pad=(0, 0), dilation = (1, 1))
    c_18 = Conv(weights["fire2/squeeze1x1_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire2/squeeze1x1_b_0"], stride=(1, 1), pad=(0, 0), dilation = (1, 1))
    c_19 = Conv(weights["conv1_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["conv1_b_0"], stride=(2, 2), pad=(0, 0), dilation = (1, 1))
    c_20 = Conv(weights["fire2/expand3x3_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire2/expand3x3_b_0"], stride=(1, 1), pad=(1, 1), dilation = (1, 1))
    c_21 = Conv(weights["fire3/expand3x3_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire3/expand3x3_b_0"], stride=(1, 1), pad=(1, 1), dilation = (1, 1))
    c_22 = Conv(weights["fire4/expand3x3_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire4/expand3x3_b_0"], stride=(1, 1), pad=(1, 1), dilation = (1, 1))
    c_23 = Conv(weights["fire5/expand3x3_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire5/expand3x3_b_0"], stride=(1, 1), pad=(1, 1), dilation = (1, 1))
    c_24 = Conv(weights["fire6/expand3x3_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire6/expand3x3_b_0"], stride=(1, 1), pad=(1, 1), dilation = (1, 1))
    c_25 = Conv(weights["fire7/expand3x3_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire7/expand3x3_b_0"], stride=(1, 1), pad=(1, 1), dilation = (1, 1))
    c_26 = Conv(weights["fire8/expand3x3_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire8/expand3x3_b_0"], stride=(1, 1), pad=(1, 1), dilation = (1, 1))
    c_27 = Conv(weights["fire9/expand3x3_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire9/expand3x3_b_0"], stride=(1, 1), pad=(1, 1), dilation = (1, 1))

    ls = Chain(Conv(weights["conv1_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["conv1_b_0"], stride=(2, 2), pad=(0, 0), dilation = (1, 1)),
            x -> relu.(x), x->maxpool(x, (3,3), pad=(0,0), stride=(2,2)),
            Conv(weights["fire2/squeeze1x1_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire2/squeeze1x1_b_0"], stride=(1, 1), pad=(0, 0), dilation = (1, 1)),
            x -> relu.(x), x->cat(relu.(c_17(x)), relu.(c_20(x)), dims=3),
            Conv(weights["fire3/squeeze1x1_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire3/squeeze1x1_b_0"], stride=(1, 1), pad=(0, 0), dilation = (1, 1)),
            x -> relu.(x), x->cat(relu.(c_15(x)), relu.(c_21(x)), dims=3),
            x->maxpool(x, (3, 3), pad=(0, 0), stride=(2, 2)),
            Conv(weights["fire4/squeeze1x1_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire4/squeeze1x1_b_0"], stride=(1, 1), pad=(0, 0), dilation = (1, 1)),
            x -> relu.(x), x->cat(relu.(c_13(x)), relu.(c_22(x)), dims=3),
            Conv(weights["fire5/squeeze1x1_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire5/squeeze1x1_b_0"], stride=(1, 1), pad=(0, 0), dilation = (1, 1)),
            x -> relu.(x), x->cat(relu.(c_11(x)), relu.(c_23(x)), dims=3),
            x->maxpool(x, (3, 3), pad=(0, 0), stride=(2, 2)),
            Conv(weights["fire6/squeeze1x1_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire6/squeeze1x1_b_0"], stride=(1, 1), pad=(0, 0), dilation = (1, 1)),
            x -> relu.(x), x->cat(relu.(c_9(x)), relu.(c_24(x)), dims=3),
            Conv(weights["fire7/squeeze1x1_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire7/squeeze1x1_b_0"], stride=(1, 1), pad=(0, 0), dilation = (1, 1)),
            x -> relu.(x), x->cat(relu.(c_7(x)), relu.(c_25(x)), dims=3),
            Conv(weights["fire8/squeeze1x1_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire8/squeeze1x1_b_0"], stride=(1, 1), pad=(0, 0), dilation = (1, 1)),
            x -> relu.(x), x->cat(relu.(c_5(x)), relu.(c_26(x)), dims=3),
            Conv(weights["fire9/squeeze1x1_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["fire9/squeeze1x1_b_0"], stride=(1, 1), pad=(0, 0), dilation = (1, 1)),
            x -> relu.(x), x->cat(relu.(c_3(x)), relu.(c_27(x)), dims=3),
            Dropout(0.5f0),
            Conv(weights["conv10_w_0"][end:-1:1,:,:,:][:,end:-1:1,:,:], weights["conv10_b_0"], stride=(1, 1), pad=(0, 0), dilation = (1, 1)),
            x -> relu.(x), x->mean(x, dims=[1,2]),
            vec, softmax
            )
#end
  Flux.testmode!(ls)
  return ls
end

struct SqueezeNet <: ClassificationModel{ImageNet.ImageNet1k}
  layers::Chain
end

SqueezeNet() = SqueezeNet(squeezenet_layers())

Base.show(io::IO, ::SqueezeNet) = print(io, "SqueezeNet()")

@treelike SqueezeNet

(m::SqueezeNet)(x) = m.layers(x)


text:
include("../../../PlotUtils/PlotUtils.jl")
include("../../../PlotUtils/imports.jl")

using Distributed
rmprocs(filter(w -> w > 1, workers()))
addprocs(9)

if length(ARGS) == 0
  results_dir = "/scratchdata/alui2/cytof/results/repfam/test-sim-6-7-14/"
else
  results_dir = ARGS[1]  # path to results directory
end

@everywhere include("../../../PlotUtils/PlotUtils.jl")
@everywhere include("../../../PlotUtils/imports.jl")

@everywhere function plot_params(samples, simdat, imgdir)
  # Create a directory for images / txt if needed.
  mkpath("$(imgdir)/txt/")

  # Extract samples
  extract(s) = map(o -> o[s], samples)
  Zs = extract(:theta__Z)
  Ws = extract(:theta__W)
  lams = extract(:theta__lam)
  alphas = extract(:theta__alpha)
  vs = extract(:theta__v)
  deltas = extract(:theta__delta)
  sig2s = extract(:theta__sig2)
  etas = extract(:theta__eta)
  omegas = extract(:omega)
  ps = [MCMC.sigmoid.(o) for o in omegas]

  # Number of samples
  I = length(sig2s[1])

  # Print Zmean
  open("$(imgdir)/txt/Zmean.txt", "w") do io
    writedlm(io, mean(Zs))
  end

  # Print Wmean
  open("$(imgdir)/txt/Wmean.txt", "w") do io
    writedlm(io, mean(Ws))
  end

  # Print R
  computeR(w) = vec(sum(w .> 0, dims=2))
  open("$(imgdir)/txt/Rcounts.txt", "w") do io
    I, K = size(Ws[1])
    Rs = map(computeR, Ws)
    Rcounts = zeros(Int, I, K)
    for R in Rs
      for i in 1:I
        Rcounts[i, R[i]] += 1
      end
    end
    writedlm(io, Rcounts)
  end

  # Plot parameters
  PlotUtils.plot_W(Ws, imgdir=imgdir, W_true=simdat[:W])
  PlotUtils.plot_v(vs, imgdir)
  PlotUtils.plot_alpha(alphas, imgdir)
  PlotUtils.plot_mus(deltas, imgdir)
  PlotUtils.plot_sig2(sig2s, imgdir, sig2_true=simdat[:sig2])
  PlotUtils.plot_p(ps, imgdir)

  # Plot y / Z
  PlotUtils.make_yz(simdat[:y], Zs, Ws, lams, imgdir, vlim=(-4,4),
                    Z_true=simdat[:Z], w_thresh=0.0)
end

@everywhere function makeplots(path_to_output)
  output_dir = splitdir(path_to_output)[1]
  println(path_to_output)

  # Create img directory
  imgdir = joinpath(output_dir, "img")
  mkpath("$(imgdir)/txt")

  # Load simulated data
  path_to_simdat = joinpath(output_dir, "simdat.bson")
  simdat = BSON.load(path_to_simdat)[:simdat]

  # Load outputs
  output = BSON.load(path_to_output)
  samples = output[:samples][1]

  # Plot parameters
  plot_params(samples, simdat, imgdir)

  # Number of burn in iterations
  nburn = output[:nburn]

  # Plot loglikelihood
  println("Loglikes ...")
  if :ll in keys(output)
    ll = output[:ll]
    PlotUtils.plot_loglike(ll, imgdir, fname="loglike.pdf")
    PlotUtils.plot_loglike(ll[(nburn + 1):end],
                           imgdir, fname="loglike_postburn.pdf")
  end

  # Extract samples
  println("Extracting samples ...")
  extract(s) = map(o -> o[s], samples)
  Zs = extract(:theta__Z)
  Ws = extract(:theta__W)
  deltas = extract(:theta__delta)
  sig2s = extract(:theta__sig2)
  etas = extract(:theta__eta)

  # Number of samples
  I = length(sig2s[1])

  # Plot missing mechanism
  println("Plotting missmechs ...")
  for i in 1:I
    PlotUtils.plot_missmech(output[:c].constants.beta, i, xlim=[-5, 0])
    plt.savefig("$(imgdir)/missmech_$(i).pdf", bbox_inches="tight")
    plt.close()
  end

  # Print missing mechanism
  open(joinpath(imgdir, "txt/beta.txt"), "w") do io
    writedlm(io, output[:c].constants.beta, ',')
  end

  # Plot data density
  println("Plotting dden ...")
  PlotUtils.plot_dden(ddens=output[:dden],
                      etas=etas, Ws=Ws, Zs=Zs, sig2s=sig2s, deltas=deltas,
                      ygrid=output[:c].constants.y_grid,
                      imgdir=imgdir, simdat=simdat)
end

### MAIN ###

# Name of output file
OUTPUT_FILE = "output.bson"

# Path to all output files
output_paths = [joinpath(root, OUTPUT_FILE)
                for (root, _, files) in walkdir(results_dir)
                if OUTPUT_FILE in files]

# Make plots in parallel
error_msg = pmap(makeplots, output_paths, on_error=identity)
println(error_msg)

println("DONE!")


text:
module Objective

include("./../constraints/geometric/constants.jl")
using JuMP
using .Constants

export objective

function objective(m, mode::Int64)

    #Other weights
    print(mode)
    if (mode == 0)
        @expression(m, minner_og_maxer,
            0
        )
    elseif (mode == 2)
        @expression(m, minner_og_maxer,
            0.1 * sum(m[:ds_mb_sum_of_boxes_below_obj]) +
            0.1 * sum(m[:middle_objective])
        )
    elseif (mode == 3)
        @expression(m, minner_og_maxer,
            0.1 * sum(m[:ds_mb_sum_of_boxes_below_obj])
        )
    elseif (mode == 4)
        @expression(m, minner_og_maxer,
            0.1 * sum(m[:middle_objective])
        )
    else
        @expression(m, minner_og_maxer,
            0.001 * sum(m[:guillotine_objective])
        )
    end

    #Objective function
    @objective(m, Min, m[:maxHeight] + 0.001 * m[:minner_og_maxer])

    return m
end

end


text:
using SpikingNeuralNetworks
using Test
SNN.@load_units

if VERSION > v"1.1"
include("ctors.jl")
end
include("chain.jl")
include("hh_net.jl")
include("hh_neuron.jl")
include("if_net.jl")
include("if_neuron.jl")
include("iz_net.jl")
include("iz_neuron.jl")
include("oja.jl")
include("rate_net.jl")
include("stdp_demo.jl")
include("stdp_demo.jl")
include("adexp_neuron.jl")
include("adexp_net.jl")

max_stars_repo_name: russelljjarvis/backupSNN

max_issues_repo_name: russelljjarvis/backupSNN

max_forks_repo_name: russelljjarvis/backupSNN


text:
function fl07()
#
#  generalization of Farmer-Loizou example:
#
#    (x-1)^4k * (x-2)^3k * (x-3)^2K * (X-4)^k
#
    fl(7)
end


text:
"""
Scaling matrix.
"""
struct ScalMat{T<:Real} <: AbstractPDMat{T}
    dim::Int
    value::T
    inv_value::T
end

ScalMat(d::Int,v::Real) = ScalMat{typeof(inv(v))}(d, v, inv(v))

### Conversion
Base.convert(::Type{ScalMat{T}}, a::ScalMat) where {T<:Real} = ScalMat(a.dim, T(a.value), T(a.inv_value))
Base.convert(::Type{AbstractArray{T}}, a::ScalMat) where {T<:Real} = convert(ScalMat{T}, a)

### Basics

dim(a::ScalMat) = a.dim
Base.Matrix(a::ScalMat) = Matrix(Diagonal(fill(a.value, a.dim)))
LinearAlgebra.diag(a::ScalMat) = fill(a.value, a.dim)
LinearAlgebra.cholesky(a::ScalMat) = cholesky(Diagonal(fill(a.value, a.dim)))

### Arithmetics

function pdadd!(r::Matrix, a::Matrix, b::ScalMat, c)
    @check_argdims size(r) == size(a) == size(b)
    if r === a
        _adddiag!(r, b.value * c)
    else
        _adddiag!(copyto!(r, a), b.value * c)
    end
    return r
end

*(a::ScalMat, c::T) where {T<:Real} = ScalMat(a.dim, a.value * c)
/(a::ScalMat{T}, c::T) where {T<:Real} = ScalMat(a.dim, a.value / c)
*(a::ScalMat, x::StridedVecOrMat) = a.value * x
\(a::ScalMat, x::StridedVecOrMat) = a.inv_value * x
Base.kron(A::ScalMat, B::ScalMat) = ScalMat( dim(A) * dim(B), A.value * B.value )

### Algebra

Base.inv(a::ScalMat) = ScalMat(a.dim, a.inv_value, a.value)
LinearAlgebra.logdet(a::ScalMat) = a.dim * log(a.value)
LinearAlgebra.eigmax(a::ScalMat) = a.value
LinearAlgebra.eigmin(a::ScalMat) = a.value

### whiten and unwhiten

function whiten!(r::StridedVecOrMat, a::ScalMat, x::StridedVecOrMat)
    @check_argdims dim(a) == size(x, 1)
    mul!(r, x, sqrt(a.inv_value))
end

function unwhiten!(r::StridedVecOrMat, a::ScalMat, x::StridedVecOrMat)
    @check_argdims dim(a) == size(x, 1)
    mul!(r, x, sqrt(a.value))
end

### quadratic forms

quad(a::ScalMat, x::AbstractVector) = sum(abs2, x) * a.value
invquad(a::ScalMat, x::AbstractVector) = sum(abs2, x) * a.inv_value

quad!(r::AbstractArray, a::ScalMat, x::StridedMatrix) = colwise_sumsq!(r, x, a.value)
invquad!(r::AbstractArray, a::ScalMat, x::StridedMatrix) = colwise_sumsq!(r, x, a.inv_value)

### tri products

function X_A_Xt(a::ScalMat, x::StridedMatrix)
    @check_argdims dim(a) == size(x, 2)
    lmul!(a.value, x * transpose(x))
end

function Xt_A_X(a::ScalMat, x::StridedMatrix)
    @check_argdims dim(a) == size(x, 1)
    lmul!(a.value, transpose(x) * x)
end

function X_invA_Xt(a::ScalMat, x::StridedMatrix)
    @check_argdims dim(a) == size(x, 2)
    lmul!(a.inv_value, x * transpose(x))
end

function Xt_invA_X(a::ScalMat, x::StridedMatrix)
    @check_argdims dim(a) == size(x, 1)
    lmul!(a.inv_value, transpose(x) * x)
end


text:
struct Graph{T}
    adjacency_lists :: Dict{T, Set{T}}
    edge_weights :: Dict{Tuple{T, T}, Float64}
end

function Graph{T}() where T
    return Graph{T}(Dict(), Dict())
end

function Graph()
    return Graph{Int}()
end

"""
    add_edge!(G, u, v, w)

Add edge (u, v) with weight w to graph G.
"""
function add_edge!(G, u, v, w)
    G.edge_weights[(u, v)] = w
    if !(haskey(G.adjacency_lists, u))
        G.adjacency_lists[u] = Set()
    end
    push!(G.adjacency_lists[u], v)
end

"""
    get_edge_weight(G, u, v)

Return edge (u, v) weight from graph G.
"""
function get_edge_weight(G, u, v)
    return G.edge_weights[(u, v)]
end

"""
    get_adjacent_nodes(G, u)

Return a set of adjacent nodes to u in graph G.
"""
function get_adjacent_nodes(G::Graph{T}, u) where T
    return get(G.adjacency_lists, u, Set{T}())
end

"""
    generate_random_graph(nvertices, edge_density, min_weight, max_weight)

Generate random graph for test purposes.
"""
function generate_random_graph(nvertices, edge_density, min_weight, max_weight)
    G = Graph()
    for u in 1:nvertices
        nadj = round(Int, nvertices * edge_density)
        while length(get_adjacent_nodes(G, u)) < nadj
            v = rand(1:nvertices)
            u == v && continue
            w = rand(min_weight:max_weight)
            add_edge!(G, u, v, w)
        end
    end
    return G
end


text:
module RegularizedGLM

using Reexport
@reexport using PenaltyFunctions

# linreg: μ = identity,                 μ′ = x -> 1
# logreg: μ = x -> 1 / (1 + exp(-x)),   μ′ = x -> exp(-x) / (1 + exp(-x))^2
# poisreg: μ = exp,                     μ′ = exp

function glmfit(x, y, μ::Function, μ′::Function)
    n, p = size(x)
    β = zeros(p)
    xβ = zeros(n)
    r = zeros(n)
    xtr = zeros(p)
    x2 = vec(sum(abs2, x, 1))
    for iter in 1:100
        for i in 1:n
            r[i] = y[i] - μ(xβ[i])
        end
        At_mul_B!(xtr, x, r)
        for j in 1:p
            β[j] = β[j] + p * xtr[j] / x2[j]
        end
        A_mul_B!(xβ, x, β)
    end
    β
end

function logreg(x, y)
    n, p = size(x)
    β = zeros(p)
    buffer_n = zeros(n)
    buffer_p = zeros(p)
    buffer_p2 = zeros(p)
    xtx_inv = 4inv(Symmetric(x'x))
    for iter in 1:50
        for i in 1:n
            buffer_n[i] = y[i] - 1 / (1 + exp(-buffer_n[i]))
        end
        At_mul_B!(buffer_p, x, buffer_n)
        A_mul_B!(buffer_p2, xtx_inv, buffer_p)
        β[:] = β + buffer_p2
        A_mul_B!(buffer_n, x, β)
    end
    β
end

# (exp(x) / (1 + exp(x))) / (1 + exp(x))

end # module


text:
# parameters
mass_bin_name = "1540_1560" #ARGS[1]#
tslice = "t1"
path_to_working_folder = "/afs/cern.ch/user/m/mimikhas/work/private/pwa_data"#"data"
path_wavelist = "src"

######################################################
# @show ARGS
push!(LOAD_PATH,"src")
using amplitudes_compass
using PWAHelper
using SDMHelper

basisfunc_mc = joinpath(path_to_working_folder,"functions_$(mass_bin_name)_$(tslice)_mc.bin")
const PsiMC = read_precalc_basis(basisfunc_mc);
Nd, Nwaves = size(PsiMC)
#
const Nbstrap_attempts = 100
const allBmatMC = [fill(0.0im,Nwaves,Nwaves) for b in 1:Nbstrap_attempts]
#
for b in 1:Nbstrap_attempts
    @show "Calculating integrals"
    randMCentries = rand(1:Nd,Nd) # here is the bootstrap
    @time allBmatMC[b] .= [sum(PsiMC[e,i]'*PsiMC[e,j] for e in randMCentries)
        for i=1:Nwaves, j=1:Nwaves] ./ Nd;
end
#

for (b,_BmatMC) in enumerate(allBmatMC)
    write_cmatrix(_BmatMC,
        joinpath(path_to_working_folder,"integrmat_$(mass_bin_name)_$(tslice)_mc_b$(b).txt"))
end

# using Plots
# using LinearAlgebra
# plot(real.([diag(BmatMC) diag(allBmatMC[1]) diag(allBmatMC[2])]))
# plot(real.([BmatMC[50,:] allBmatMC[1][50,:] allBmatMC[2][50,:]]))


text:
import Base.convert

##############################################################################
##
## References
##
##############################################################################

# - https://en.bitcoin.it/wiki/Protocol_specification

##############################################################################
##
## Notes
##
##############################################################################

# MESSAGE STRUCTURE - A typical message has the following structure:
#
# Field size | Description | Data type | Comments
# ----------------------------------------------------------------------------
#          4 | magic       | uint32_t  | Value indicating origin network
#         12 | command     | char[12]  | ASCII string identifying packet content
#          4 | length      | uint32_t  | Length of payload in number of bytes
#          4 | checksum    | uint32_t  | First 4 bytes of sha256(sha256(payload))
#          ? | payload     | uchar[?]  | The actual data
#
# Almost all integers are encoded in little endian. Only IP or port number are
# encoded big endian.

# MESSAGE TYPES - The following are the current message types
#
# version:    When a node creates an outgoing connection, it will immediately
#             advertise its version. The remote node will respond with its version.
#             No further communication is possible until both peers have exchanged
#             their version.
# verack:     The verack message is sent in reply to version. This message consists
#             of only a message header with the command string "verack".
# addr:       Provide information on known nodes of the network. Non-advertised
#             nodes should be forgotten after typically 3 hours.
# inv:        Allows a node to advertise its knowledge of one or more objects.
#             It can be received unsolicited, or in reply to getblocks.
# getdata:    getdata is used in response to inv, to retrieve the content of a
#             specific object, and is usually sent after receiving an inv packet,
#             after filtering known elements. It can be used to retrieve
#             transactions, but only if they are in the memory pool or relay set -
#             arbitrary access to transactions in the chain is not allowed to
#             avoid having clients start to depend on nodes having full
#             transaction indexes (which modern nodes do not).
# notfound:   notfound is a response to a getdata, sent if any requested data
#             items could not be relayed, for example, because the requested
#             transaction was not in the memory pool or relay set.
# getblocks
# getheaders
# tx
# block
# headers
# getaddr
# mempool
# checkorder
# submitorder
# reply
# ping
# pong
# reject
# filterload
# filteradd
# filterclear
# merkleblock
# alert

# const commands = ["version", "verack", "addr", "inv", "getdata", "notfound",
#                   "getblocks", "getheaders", "tx", "block", "headers",
#                   "getaddr", "mempool", "checkorder", "submitorder", "reply",
#                   "ping", "pong", "reject", "filterload", "filteradd",
#                   "filterclear", "merkleblock", "alert"]

include("tx.jl")
include("version.jl")

# Define the known magic values
const magic_mainnet  = 0xd9b4bef9
const magic_testnet  = 0xdab5bffa
const magic_testnet3 = 0x0709110b
const magic_namecoin = 0xfeb4bef9

type Message
  magic::Uint32               # Network identifier
  command::Array{Uint8}       # Message command, right padded with \0 to 12 bytes
  payload::Array{Uint8}       # Payload data

  function Message(magic::Uint32, command::String, payload::Array{Uint8})
    # Right pad the command string with \0 to create 12-byte array
    append!(command.data, zeros(Char, 12 - length(command.data)))

    # Calculate checksum
    checksum = get_checksum(payload)

    new(magic, command.data, payload)
  end
end

function convert(::Type{Array{Uint8}}, x::Message)
  result = Array(Uint8, 0)

  append!(result, Crypto.int2oct(reverse_endian(x.magic)))
  append!(result, x.command)
  append!(result, Crypto.int2oct(reverse_endian(uint32(length(x.payload)))))
  append!(result, payload)

  return result
end


text:
addprocs(1)

using JLD
using BaseBenchmarks
using BenchmarkTools

# re-tune the entire suite

function rewrite_params_file(paramsgroup)
    jldopen(BaseBenchmarks.PARAMS_PATH, "w") do file
        for (id, suite) in paramsgroup
            write(file, id, suite)
        end
    end
end

function rewrite_params_file(paramsgroup, id)
    old = BenchmarkTools.load(BaseBenchmarks.PARAMS_PATH)
    jldopen(BaseBenchmarks.PARAMS_PATH, "w") do file
        for (oldid, oldsuite) in old
            if oldid == id
                write(file, oldid, paramsgroup)
            else
                write(file, oldid, oldsuite)
            end
        end
    end
end

# retune the whole suite #
#------------------------#
BaseBenchmarks.loadall!(tune = false)
warmup(BaseBenchmarks.SUITE)
tune!(BaseBenchmarks.SUITE; verbose = true)
rewrite_params_file(params(BaseBenchmarks.SUITE))

# retune an individual group (using "linalg" as an example) #
#-----------------------------------------------------------#
# id = "linalg"
# BaseBenchmarks.load!(id; tune = false)
# group = BaseBenchmarks.SUITE[id]
# warmup(group)
# tune!(group; verbose = true)
# rewrite_params_file(params(group), id)


text:
using Zygote

Zygote.@adjoint argmax(xs; dims = :) = argmax(xs; dims = dims), _ -> nothing


text:
mutable struct MeanSquaredError <: DiffableFunction
    grad_field :: GradField
end

function forward(::MeanSquaredError, y_true, y_pred)
    diff  = y_true .- y_pred
    loss = Base.sum((x -> x^2), diff) / length(diff)
    return loss
end

function backward(f::MeanSquaredError, gy)
    y_true, y_pred = f.grad_field.inputs
    diff = y_true .- y_pred
    gy = broadcast_to(gy, size(diff.values))
    gx1 = gy .* diff .* (2. / length(diff.values))
    gx2 = -gx1
    return (gx1, gx2)
end

"""
    mean_squared_error(y_true, y_pred)

# Examples
```julia-repl
julia> y_true = Variable([10, 20, 30])
name: nothing
values: [10, 20, 30]
creator: User-Defined(nothing)

julia> y_pred = Variable([12, 20, 25])
name: nothing
values: [12, 20, 25]
creator: User-Defined(nothing)

julia> mean_squared_error(y_true, y_pred)
name: nothing
values: 9.666666666666666
creator: JITrench.MeanSquaredError
```
"""
mean_squared_error(y_true, y_pred) = MeanSquaredError(GradField())(y_true, y_pred)


text:
"""
    AbstractZero <: AbstractDifferential

This is zero-like differential types.
If a AD system encounter a propagator taking as input only subtypes of `AbstractZero` then
it can stop performing any AD operations, as all propagator are linear functions, and thus
the final result will be zero.

All `AbstractZero` subtypes are singleton types.
There are two of them [`Zero()`](@ref) and [`DoesNotExist()`](@ref).
"""
abstract type AbstractZero <: AbstractDifferential end
Base.iszero(::AbstractZero) = true

Base.iterate(x::AbstractZero) = (x, nothing)
Base.iterate(::AbstractZero, ::Any) = nothing

Base.Broadcast.broadcastable(x::AbstractZero) = Ref(x)
Base.Broadcast.broadcasted(::Type{T}) where T<:AbstractZero = T()

# Linear operators
Base.adjoint(z::AbstractZero) = z
Base.transpose(z::AbstractZero) = z
Base.:/(z::AbstractZero, ::Any) = z

"""
    Zero() <: AbstractZero

The additive identity for differentials.
This is basically the same as `0`.
A derivative of `Zero()`. does not propagate through the primal function.
"""
struct Zero <: AbstractZero end

extern(x::Zero) = false  # false is a strong 0. E.g. `false * NaN = 0.0`

Base.eltype(::Type{Zero}) = Zero

Base.zero(::AbstractDifferential) = Zero()
Base.zero(::Type{<:AbstractDifferential}) = Zero()

"""
    DoesNotExist() <: AbstractZero

This differential indicates that the derivative does not exist.
It is the differential for a Primal type that is not differentiable.
Such an Integer, or Boolean (when not being used as a represention of a value that normally
would be a floating point.)
The only valid way to pertube such a values is to not change it at all.
As such, `DoesNotExist` is functionally identical to `Zero()`,
but provides additional semantic information.

If you are adding this differential to a primal then something is wrong.
A optimization package making use of this might like to check for such a case.

!!! note:
    This does not indicate that the derivative it is not implemented,
    but rather that mathematically it is not defined.

This mostly shows up as the deriviative with respect to dimension, index, or size
arguments.
```
    function rrule(fill, x, len::Int)
        y = fill(x, len)
        fill_pullback(ȳ) = (NO_FIELDS, @thunk(sum(Ȳ)), DoesNotExist())
        return y, fill_pullback
    end
```
"""
struct DoesNotExist <: AbstractZero end

function extern(x::DoesNotExist)
    throw(ArgumentError("Derivative does not exit. Cannot be converted to an external type."))
end


text:
## Example 2 from paper
# "Lasso Regularization Paths for NARMAX Models via Coordinate Descent"
# Antonio H. Ribeiro and Luis A. Aguirre
using NarmaxLasso
using Plots
plotlyjs()

## Signal generation
# Defining nonlinear system
function nonlinear_function(y, u, v)
    #return 0.5*y[1] + 0.2*u[1] + 0.1*v[1]
    return ((0.8 - 0.5*exp(-y[1]^2))*y[1] - (0.3 + 0.9*exp(-y[1]^2))*y[2] + u[1]
            + 0.2*u[2] + 0.1*u[1]*u[2] + 0.3*v[1] + 0.3*v[2])
end
yterms0 = [1, 2]
uterms0 = [1, 2]
vterms0 = [1, 2]

# Generate signals
u = generate_random_input(1000, 5)
y = simulate(u, nonlinear_function, yterms0, uterms0, vterms0;
             σv=0.5, seed=1)

## Model Estimation
# Generate all candidate terms
ny = 3
nu = 3
nv = 2
order = collect(1:2)
mdl = generate_all(NarmaxRegressors, Monomial, ny, nu, nv, order)

# Lasso regression
result = narmax_lasso(y, u, mdl; λ_min_ratio=1e-6)

# Plot result
p1 = plot(result, xscale=:log10, ylims=(-1.2, 1.2))

## Validation
# Generate signals
uv = generate_random_input(500, 5, seed=10)
yv = simulate(uv, nonlinear_function, yterms0, uterms0, vterms0;
                σv=0.5, seed=6)

# Compute Mean Absolute Erro
n = length(result.λ)
mae = Vector{Float64}(n)
for i = 1:n
    ysv = simulate(uv, result.mdl, result.β[:, i])
    mae[i] = mean(abs(ysv-yv))
end
vline!(p1, [result.λ[indmin(mae)]], color=:black, line=:dash)

## Plot number of nonzero parameters
p2 = plot(result.λ, result.df, line=:steppost,
          xscale=:log10, lab="non-zero terms", xflip=true, legend=false,
          ylabel="Non-zero terms",
          xlabel="lambda", grid=false, lw=0.7,
          size=(550, 300))
vline!(p2, [result.λ[indmin(mae)]], color=:black, line=:dash)

## Mean Absolut Error
p3 = plot(result.λ, mae, xscale=:log10, lab="test set", ylims=(1, 1.3),
         legend=false, ylabel="MAE", xlabel="lambda",
         grid=false, lw=0.7, size=(550, 300), xflip=true)
vline!(p3, [result.λ[indmin(mae)]], color=:black, line=:dash)

plot(p1, p2, p3, layout=(3, 1))