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
#############################################################################
# Defines functions used to parse the input optimization problem into
# a solvable form including routines used to classify input problems as a
# LP, SOCP, MILP, MISOCP, and convex problem types.
#############################################################################

"""
    add_objective!

Adds objective function (if any) to the parsed problem.
"""
add_objective!(m::ParsedProblem, f) = nothing
add_objective!(m::ParsedProblem, f::VI) = (m._objective = f; nothing)
add_objective!(m::ParsedProblem, f::SAF) = (m._objective = AffineFunctionIneq(f, LT_ZERO); nothing)
add_objective!(m::ParsedProblem, f::SQF) = (m._objective = BufferedQuadraticIneq(f, LT_ZERO); nothing)

"""
    add_nonlinear_evaluator!

    Adds a Evaluator structure if nonlinear terms are attached.
"""
add_nonlinear_evaluator!(m::GlobalOptimizer, evaluator::Nothing) = nothing
function add_nonlinear_evaluator!(m::GlobalOptimizer, d::JuMP.NLPEvaluator)
    m._working_problem._relaxed_evaluator = Evaluator()
    relax_eval = m._working_problem._relaxed_evaluator
    relax_eval.user_operators = OperatorRegistry(d.model.nlp_data.user_operators)
    relax_eval.subgrad_tol    = m._parameters.subgrad_tol
    m._nonlinear_evaluator_created = true
    return
end
add_nonlinear_evaluator!(m::GlobalOptimizer, nldata) = add_nonlinear_evaluator!(m, nldata.evaluator)
add_nonlinear_evaluator!(m::GlobalOptimizer) = add_nonlinear_evaluator!(m, m._input_problem._nlp_data)

function link_subexpression_dicts!(m::GlobalOptimizer)
    evaluator = m._working_problem._relaxed_evaluator
    n_subexpr = length(evaluator.subexpressions)

    dn = Dict{Int,Float64}()
    din = Dict{Int,Bool}()
    for i = 1:n_subexpr
        dn[i] = 0.0
        din[i] = false
    end
    for ex in evaluator.subexpressions
        mctyp = mc_type(ex)
        ds = Dict{Int,mctyp}()
        di = Dict{Int,mctyp}()
        for i = 1:n_subexpr
            ds[i] = zero(mctyp)
            di[i] = zero(mctyp)
        end
        copy_subexpr!(ex.relax_cache, ds, dn, din, di)
    end
    if m._working_problem._objective isa BufferedNonlinearFunction
        ex = m._working_problem._objective.ex
        mctyp = mc_type(ex)
        ds = Dict{Int,mctyp}()
        di = Dict{Int,mctyp}()
        for i = 1:n_subexpr
            ds[i] = zero(mctyp)
            di[i] = zero(mctyp)
        end
        copy_subexpr!(ex.relax_cache, ds, dn, din, di)
    end
    for f in m._working_problem._nonlinear_constr
        mctyp = mc_type(f.ex)
        ds = Dict{Int,mctyp}()
        di = Dict{Int,mctyp}()
        for i = 1:n_subexpr
            ds[i] = zero(mctyp)
            di[i] = zero(mctyp)
        end
        copy_subexpr!(f.ex.relax_cache, ds, dn, din, di)
    end
    return
end

add_nonlinear!(m::GlobalOptimizer, evaluator::Nothing) = nothing
function add_nonlinear!(m::GlobalOptimizer, evaluator::JuMP.NLPEvaluator)

    add_nonlinear_evaluator!(m, m._input_problem._nlp_data.evaluator)

    nlp_data = m._input_problem._nlp_data

    MOI.initialize(evaluator, Symbol[:Grad, :Jac, :ExprGraph])

    user_operator_registry = OperatorRegistry(evaluator.model.nlp_data.user_operators)

    # set nlp data structure
    m._working_problem._nlp_data = nlp_data
    mul_relax = m._parameters.mul_relax_style
    if mul_relax == 1
        rtype = Relax()
        renum = STD_RELAX
        ruse_apriori = true
    elseif mul_relax == 2
        rtype = RelaxAA()
        renum = MC_AFF_RELAX
        ruse_apriori = true
    elseif mul_relax == 3
        rtype = RelaxMulEnum()
        renum = MC_ENUM_RELAX
        ruse_apriori = true
    else
        rtype = Relax()
        renum = STD_RELAX
        ruse_apriori = false
    end

    # add subexpressions (assumes they are already ordered by JuMP)
    # creates a dictionary that lists the subexpression sparsity
    # by search each node for variables dict[2] = [2,3] indicates
    # that subexpression 2 depends on variables 2 and 3
    # this is referenced when subexpressions are called by other
    # subexpressions or functions to determine overall sparsity
    # the sparsity of a function is the collection of indices
    # in all participating subexpressions and the function itself
    # it is necessary to define this as such to enable reverse
    # McCormick constraint propagation
    relax_evaluator = m._working_problem._relaxed_evaluator
    relax_evaluator.relax_type = renum
    dict_sparsity = Dict{Int,Vector{Int}}()
    if length(evaluator.model.nlp_data.nlexpr) > 0      # should check for nonlinear objective, constraint
        for i = 1:length(evaluator.subexpressions)
            subexpr = evaluator.subexpressions[i]
            nlexpr = NonlinearExpression!(m._auxillary_variable_info, rtype, subexpr, MOI.NLPBoundsPair(-Inf, Inf),
                                          dict_sparsity, i, evaluator.subexpression_linearity,
                                          user_operator_registry, evaluator.model.nlp_data.nlparamvalues,
                                          m._parameters.relax_tag, ruse_apriori; is_sub = true)
            push!(relax_evaluator.subexpressions, nlexpr)
        end
    end

    # scrubs udf functions using Cassette to remove odd data structures...
    # alternatively convert udfs to JuMP scripts...
    m._parameters.presolve_scrubber_flag && Script.scrub!(m._working_problem._nlp_data)
    if m._parameters.presolve_to_JuMP_flag
        Script.udf_loader!(m)
    end

    # add nonlinear objective
    if evaluator.has_nlobj
        m._working_problem._objective = BufferedNonlinearFunction(m._auxillary_variable_info, rtype, evaluator.objective, MOI.NLPBoundsPair(-Inf, Inf),
                                                                 dict_sparsity, evaluator.subexpression_linearity,
                                                                 user_operator_registry,
                                                                 evaluator.model.nlp_data.nlparamvalues,
                                                                 m._parameters.relax_tag, ruse_apriori)
    end

    # add nonlinear constraints
    constraint_bounds = m._working_problem._nlp_data.constraint_bounds
    for i = 1:length(evaluator.constraints)
        constraint = evaluator.constraints[i]
        bnds = constraint_bounds[i]
        push!(m._working_problem._nonlinear_constr, BufferedNonlinearFunction(m._auxillary_variable_info, rtype, constraint, bnds, dict_sparsity,
                                                                              evaluator.subexpression_linearity,
                                                                              user_operator_registry,
                                                                              evaluator.model.nlp_data.nlparamvalues,
                                                                              m._parameters.relax_tag, ruse_apriori))
    end

    relax_evaluator.subexpressions_eval = fill(false, length(relax_evaluator.subexpressions))
    return
end
add_nonlinear!(m::GlobalOptimizer, nldata) = add_nonlinear!(m, nldata.evaluator)

"""
$(TYPEDSIGNATURES)

Adds a Evaluator, nonlinear functions, and populates each appropriately.
"""
add_nonlinear!(m::GlobalOptimizer) = add_nonlinear!(m, m._input_problem._nlp_data)

function reform_epigraph_min!(m::GlobalOptimizer, d::ParsedProblem, f::VI)
    flag = m._input_problem._optimization_sense == MOI.MAX_SENSE
    d._objective = AffineFunctionIneq(f, is_max = flag)
    d._objective_saf = SAF([SAT(flag ? -1.0 : 1.0, VI(f.value))], 0.0)
    return
end
function reform_epigraph_min!(m::GlobalOptimizer, d::ParsedProblem, f::AffineFunctionIneq)
    d._objective_saf = m._input_problem._objective
    if m._input_problem._optimization_sense == MOI.MAX_SENSE
        MOIU.operate!(-, Float64, d._objective_saf)
        d._objective = AffineFunctionIneq([(-i,j) for (i,j) in f.terms], -f.constant, f.len)
    else
        d._objective = f
    end
    return
end

function add_η!(m::ParsedProblem)
    m._variable_count += 1
    push!(m._variable_info, VariableInfo{Float64}())
    return m._variable_count
end

set_variable_values!(wp, v) = set_variable_values!(wp._relaxed_evaluator, v)
reform_epigraph_min!(m::GlobalOptimizer, d::ParsedProblem, f::Nothing) = nothing
function reform_epigraph_min!(m::GlobalOptimizer, d::ParsedProblem, f::BufferedQuadraticIneq)
    ip = m._input_problem
    wp = m._working_problem

    vi = wp._variable_info
    q = _variable_num(FullVar(), m)
    v = VariableValues{Float64}(x = mid.(vi), x0 = mid.(vi),
                                lower_variable_bounds = lower_bound.(vi),
                                upper_variable_bounds = upper_bound.(vi),
                                node_to_variable_map = [i for i in 1:q],
                                variable_to_node_map = [i for i in 1:q])
    set_variable_values!(wp, v)

    n = NodeBB(lower_bound.(vi), upper_bound.(vi), is_integer.(vi))
    m._current_node = n
    set_node!(wp._relaxed_evaluator, n)

    ηi = add_η!(d)
    if !isnothing(ip._nlp_data)
        @variable(ip._nlp_data.evaluator.model, η)
        @objective(ip._nlp_data.evaluator.model, Min, η)
    end

    sqf_obj = copy(m._input_problem._objective)
    if !_is_input_min(m)
        MOIU.operate!(-, Float64, sqf_obj)
    end
    d._objective_saf = SAF([SAT(1.0, VI(ηi))], 0.0)
    push!(sqf_obj.affine_terms, SAT(-1.0, VI(ηi)))
    _constraints(ip, SQF, LT)[CI{SQF,LT}(ip._constraint_count += 1)] = (sqf_obj, LT(0.0))
    push!(wp._sqf_leq, BufferedQuadraticIneq(sqf_obj, LT(0.0)))

    f.buffer[ηi] = 0.0
    f.len += 1
    if !_is_input_min(m)
        MOIU.operate!(-, Float64, f.func)
    end
    MOIU.operate!(-, Float64, f.func, VI(ηi))
    push!(f.saf.terms, SAT(-1.0, VI(ηi)))
    m._obj_var_slack_added = true
    return
end
function reform_epigraph_min!(m::GlobalOptimizer, d::ParsedProblem, f::BufferedNonlinearFunction)
    ip = m._input_problem
    wp = m._working_problem

    vi = wp._variable_info
    q = _variable_num(FullVar(), m)
    v = VariableValues{Float64}(x = mid.(vi), x0 = mid.(vi),
                                lower_variable_bounds = lower_bound.(vi),
                                upper_variable_bounds = upper_bound.(vi),
                                node_to_variable_map = [i for i in 1:q],
                                variable_to_node_map = [i for i in 1:q])

    ηi = add_η!(d)
    @variable(ip._nlp_data.evaluator.model, η)
    wp._objective_saf = SAF([SAT(1.0, VI(ηi))], 0.0)

    nd = ip._nlp_data.evaluator.model.nlp_data.nlobj.nd
    if !_is_input_min(m)
        pushfirst!(nd, NodeData(JuMP._Derivatives.CALLUNIVAR, 2, 1))
        pushfirst!(nd, NodeData(JuMP._Derivatives.CALL, 2, -1))
        nd[3] = NodeData(nd[3].nodetype, nd[3].index, 2)
        for i = 4:length(nd)
            @inbounds nd[i] = NodeData(nd[i].nodetype, nd[i].index, nd[i].parent + 2)
        end
    else
        pushfirst!(nd, NodeData(JuMP._Derivatives.CALL, 2, -1))
        nd[2] = NodeData(nd[2].nodetype, nd[2].index, 1)
        for i = 3:length(nd)
            @inbounds nd[i] = NodeData(nd[i].nodetype, nd[i].index, nd[i].parent + 1)
        end
    end
    push!(nd, NodeData(JuMP._Derivatives.VARIABLE, ηi, 1))
    nlobj = ip._nlp_data.evaluator.model.nlp_data.nlobj
    nlexpr = JuMP._NonlinearExprData(copy(nlobj.nd), copy(nlobj.const_values))
    nlcons = JuMP._NonlinearConstraint(nlexpr, -Inf, 0.0)

    ip._nlp_data.evaluator.model.nlp_data.nlobj = nothing
    ip._nlp_data.evaluator.has_nlobj = false
    push!(ip._nlp_data.evaluator.model.nlp_data.nlconstr, nlcons)
    constraint_bounds = ip._nlp_data.constraint_bounds
    push!(constraint_bounds, MOI.NLPBoundsPair(-Inf, 0.0))
    ip._nlp_data = MOI.NLPBlockData(constraint_bounds, ip._nlp_data.evaluator, false)

    empty!(d._relaxed_evaluator.subexpressions)
    empty!(d._nonlinear_constr)
    add_nonlinear!(m)
    m._obj_var_slack_added = true
    return
end

"""
$(TYPEDSIGNATURES)

Performs an epigraph reformulation assuming the working_problem is a minimization problem.
"""
function reform_epigraph_min!(m::GlobalOptimizer)
    ip = m._input_problem
    m._branch_variables = fill(false, m._working_problem._variable_count)
    m._obj_mult = (ip._optimization_sense == MOI.MAX_SENSE) ? -1.0 : 1.0
    reform_epigraph_min!(m, m._working_problem, m._working_problem._objective)
end

function check_set_is_fixed(v::VariableInfo)
    v.is_fixed && return true
    v.is_fixed = x.lower_bound === x.upper_bound
    return v.is_fixed
end

"""
$(TYPEDSIGNATURES)

Detects any variables set to a fixed value by equality or inequality constraints
and populates the `_fixed_variable` storage array.
"""
label_fixed_variables!(m::GlobalOptimizer) = map!(x -> check_set_is_fixed(x), m._fixed_variable, m._working_problem._variable_info)

"""
$(TYPEDSIGNATURES)

Detects any variables participating in nonconvex terms and populates the
`_branch_variables` storage array.
"""
function label_branch_variables!(m::GlobalOptimizer)

    wp = m._working_problem
    m._user_branch_variables = !isempty(m._parameters.branch_variable)
    if m._user_branch_variables
        m._branch_variables = m._parameters.branch_variable
    else
        m._branch_variables = fill(false, m._working_problem._variable_count)
        for f in wp._sqf_leq, t in f.func.quadratic_terms
            m._branch_variables[t.variable_1.value] = true
            m._branch_variables[t.variable_2.value] = true
        end
        for f in wp._sqf_eq, t in f.func.quadratic_terms
            m._branch_variables[t.variable_1.value] = true
            m._branch_variables[t.variable_2.value] = true
        end
        for f in wp._nonlinear_constr, i in sparsity(f)
             m._branch_variables[i] = true
        end
        if wp._objective isa BufferedQuadraticIneq
            for t in wp._objective.func.quadratic_terms
                m._branch_variables[t.variable_1.value] = true
                m._branch_variables[t.variable_2.value] = true
            end
        elseif wp._objective isa BufferedNonlinearFunction
            for i in sparsity(wp._objective)
                m._branch_variables[i] = true
            end
        end
    end

    # add a map of branch/node index to variables in the continuous solution
    for i = 1:wp._variable_count
        if is_fixed(wp._variable_info[i])
            m._branch_variables[i] = false
        elseif m._branch_variables[i] || wp._variable_info[i].is_integer
            push!(m._branch_to_sol_map, i)
        elseif i == wp._variable_count
            m._branch_variables[i] = false
        end
    end

    # creates reverse map
    m._sol_to_branch_map = zeros(wp._variable_count)
    for i = 1:length(m._branch_to_sol_map)
        j = m._branch_to_sol_map[i]
        m._sol_to_branch_map[j] = i
    end

    # adds branch solution to branch map to evaluator
    vnum = wp._variable_count
    initialize!(m._branch_cost, length(m._branch_to_sol_map))
    l = lower_bound.(m._working_problem._variable_info)
    u = upper_bound.(m._working_problem._variable_info)
    v = VariableValues{Float64}(x = zeros(vnum),
                                x0 = zeros(vnum),
                                lower_variable_bounds = l,
                                upper_variable_bounds = u,
                                node_to_variable_map = m._branch_to_sol_map,
                                variable_to_node_map = m._sol_to_branch_map)

    wp._relaxed_evaluator.variable_values = v
    (wp._objective isa BufferedNonlinearFunction) && set_variable_storage!(wp._objective, v)
    foreach(i -> set_variable_storage!(i, v), wp._nonlinear_constr)
    foreach(i -> set_variable_storage!(i, v), wp._relaxed_evaluator.subexpressions)
    return
end

function variable_load_parse!(m::Optimizer, ::Type{VI}, ::Type{T}) where T
    wp = m._global_optimizer._working_problem = m._working_problem
    for (i, v) in enumerate(values(_constraints(m, VI, T)))
        wp._variable_info[v[1].value] = VariableInfo(wp._variable_info[v[1].value], v[2])
    end
    return
end

"""
$(TYPEDSIGNATURES)

Translates input problem to working problem. Checks and optional manipulation is left to the presolve stage.
"""
function initial_parse!(m::Optimizer{R,S,T}) where {R,S,T}

    # reset initial time and solution statistics
    m._global_optimizer._time_left = m._parameters.time_limit

    ip = m._global_optimizer._input_problem = m._input_problem
    wp = m._global_optimizer._working_problem = m._working_problem
    m._global_optimizer._parameters = m._parameters

    # add variables to working model
    wp._variable_info = VariableInfo{Float64}[VariableInfo{Float64}() for i=1:ip._variable_count]
    variable_load_parse!(m, VI, LT)
    variable_load_parse!(m, VI, GT)
    variable_load_parse!(m, VI, ET)
    variable_load_parse!(m, VI, ZO)
    variable_load_parse!(m, VI, MOI.Integer)
    wp._variable_count = ip._variable_count

    for (f, s) in values(ip._linear_leq_constraints)
        push!(wp._saf_leq, AffineFunctionIneq(f, s))
    end
    for (f, s) in values(ip._linear_geq_constraints)
        push!(wp._saf_leq, AffineFunctionIneq(f, s))
    end
    for (f, s) in values(ip._linear_eq_constraints)
        push!(wp._saf_eq, AffineFunctionEq(f, s))
    end

    for (f, s) in values(ip._quadratic_leq_constraints)
        push!(wp._sqf_leq, BufferedQuadraticIneq(f, s))
    end
    for (f, s) in values(ip._quadratic_geq_constraints)
        push!(wp._sqf_leq, BufferedQuadraticIneq(f, s))
    end
    for (f, s) in values(ip._quadratic_eq_constraints)
        push!(wp._sqf_eq, BufferedQuadraticEq(f, s))
    end

    add_objective!(wp, ip._objective)    # set objective function
    add_nonlinear!(m._global_optimizer)  # add nonlinear constraints, evaluator, subexpressions

    # converts a maximum problem to a minimum problem (internally) if necessary
    # this is placed after adding nonlinear functions as this prior routine
    # copies the nlp_block from the input_problem to the working problem
    reform_epigraph_min!(m._global_optimizer)
    label_fixed_variables!(m._global_optimizer)
    label_branch_variables!(m._global_optimizer)

    link_subexpression_dicts!(m._global_optimizer)

    # updates run and parse times
    new_time = time() - m._global_optimizer._start_time
    m._global_optimizer._parse_time = new_time
    m._global_optimizer._run_time = new_time
    return
end

"""
Classifies the problem type
"""
function parse_classify_problem!(m::GlobalOptimizer)

    ip = m._input_problem

    nl_expr_num = 0
    if isnothing(ip._objective) && !isnothing(ip._nlp_data)
        if ip._nlp_data.has_objective
            nl_expr_num += 1
            has_objective = true
        else
            has_objective = false
        end
    elseif isnothing(ip._objective) && isnothing(ip._nlp_data)
        has_objective = false
    end
    nl_expr_num += length(m._working_problem._nonlinear_constr)
    cone_constr_num = length(ip._conic_second_order)
    quad_constr_num = length(ip._quadratic_leq_constraints) +
                      length(ip._quadratic_geq_constraints) +
                      length(ip._quadratic_eq_constraints)

    if !isnothing(ip._objective)
        has_objective = true
    end
    has_int_var = !iszero(length(ip._vi_zo_constraints) + length(ip._vi_int_constraints))

    lin_or_sv_obj = (ip._objective isa VI || ip._objective isa SAF || !has_objective)
    relaxed_supports_soc = false

    if (cone_constr_num == 0) && (quad_constr_num == 0) && (nl_expr_num == 0) && lin_or_sv_obj && !has_int_var
        m._working_problem._problem_type = LP()
    elseif (cone_constr_num == 0) && (quad_constr_num == 0) && (nl_expr_num == 0) && lin_or_sv_obj
        m._working_problem._problem_type = MILP()
    elseif (quad_constr_num == 0) && relaxed_supports_soc && (nl_expr_num == 0) && lin_or_sv_obj && !has_int_var
        m._working_problem._problem_type = SOCP()
    else
        m._working_problem._problem_type = MINCVX()
    end
    return
end

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
# Call Syntax
person = @shared_var Citizen(
            name::String = "Amin",
            number::Float64 = 20.0,
        )

person2 = Citizen(name = "Not-Amin", number = 1)

@test person.name == person2.name
@test person2.number == person2.number

ext: jl


text:
function ChiSquaredDistribution(arg0::jdouble)
    return ChiSquaredDistribution((jdouble,), arg0)
end

function ChiSquaredDistribution(arg0::jdouble, arg1::jdouble)
    return ChiSquaredDistribution((jdouble, jdouble), arg0, arg1)
end

function cumulative_probability(obj::ChiSquaredDistribution, arg0::jdouble)
    return jcall(obj, "cumulativeProbability", jdouble, (jdouble,), arg0)
end

function density(obj::ChiSquaredDistribution, arg0::jdouble)
    return jcall(obj, "density", jdouble, (jdouble,), arg0)
end

function get_degrees_of_freedom(obj::ChiSquaredDistribution)
    return jcall(obj, "getDegreesOfFreedom", jdouble, ())
end

function get_numerical_mean(obj::ChiSquaredDistribution)
    return jcall(obj, "getNumericalMean", jdouble, ())
end

function get_numerical_variance(obj::ChiSquaredDistribution)
    return jcall(obj, "getNumericalVariance", jdouble, ())
end

function get_support_lower_bound(obj::ChiSquaredDistribution)
    return jcall(obj, "getSupportLowerBound", jdouble, ())
end

function get_support_upper_bound(obj::ChiSquaredDistribution)
    return jcall(obj, "getSupportUpperBound", jdouble, ())
end

function is_support_connected(obj::ChiSquaredDistribution)
    return jcall(obj, "isSupportConnected", jboolean, ())
end

function log_density(obj::ChiSquaredDistribution, arg0::jdouble)
    return jcall(obj, "logDensity", jdouble, (jdouble,), arg0)
end


text:
#=
Julia (v1) function to get frequencies of words in a directory of TXT files.
Earl K. Brown, ekbrown byu edu (add appropriate characters to create email)
=#

using DataFrames

function get_frequency(pathway, sort_by = "freq")
    #=
    The function takes as input a string with the pathway to the directory with .txt files and returns as output a Julia DataFrame.

    pathway: A string with the (absolute or relative) pathway to the directory with the .txt files; other file types are ignored.
    sort_by: A string indicating how the frequencies should be ordered, from among "freq" (frequency) in descending order, or "alpha" (alphabetical order) in ascending order (default is "freq").
    return value: A Julia DataFrame with three columns: word, raw frequency, log frequency (base 10).
    =#

    # get filenames
    cd(pathway)
    filenames = [f for f in readdir() if occursin(r"\.txt$"i, f)]

    # loop over files, collecting frequencies along the way
    freqs = Dict{String, Int64}()
    for filename in filenames
        open(filename) do infile
            cur_file = uppercase(read(infile, String))  # convert to uppercase
            wds = split(cur_file, r"[^-'A-ZÁÉÍÓÚÜÑ]+"i)  # get words in current file
            for wd in wds
                freqs[wd] = get(freqs, wd, 0) + 1
            end  # next word
        end  # close file connection
    end  # next filename

    # push dictionary to data frame
    freqs_df = DataFrame(word = String[], freq = Int64[], log_freq = Float64[])
    for (k, v) in freqs
        push!(freqs_df, [k, v, log10(v)])
    end

    # sort as desired by user of function
    if sort_by == "freq"
        sort!(freqs_df, [order(:freq, rev = true), :word])
    elseif sort_by == "alpha"
        sort!(freqs_df, [order(:word)])
    end

    return freqs_df
end

### test the function
dir = "/Users/ekb5/Corpora/USA/California/Salinas/transcripts/"
sort_by = "freq"
@time freqs = get_frequency(dir, sort_by)
println(freqs[1:10,:])


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, BinaryBuilderBase, Pkg

name = "HMMER"
version_string = "3.1b2"
#version = v"3.1.0"
version = let ver = VersionNumber(version_string)
    VersionNumber(ver.major, ver.minor, ver.patch)
end
# Collection of sources required to complete build
sources = [
    ArchiveSource("http://eddylab.org/software/hmmer/hmmer-$(version_string).tar.gz",
        "dd16edf4385c1df072c9e2f58c16ee1872d855a018a2ee6894205277017b5536")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/hmmer-*/

export CPPFLAGS="-I${includedir}"
update_configure_scripts
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} \
    --with-gsl --enable-threads --enable-sse
make -j${nproc}
make install

cd easel && make install
mv LICENSE LICENSE.easel

install_license LICENSE.easel
install_license ../LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
#
# Notes:
# - configure demands either sse (i686/x86_64) or vmx (powerpc) vector
#   instructions
# - windows compile fails, can't find syslog.h
# - compile fails on powerpc64le:
#   "HMMER3 Altivec/VMX only supports bigendian platforms: e.g. ppc64 not ppc64le"
platforms = supported_platforms(; exclude = p -> Sys.iswindows(p) || proc_family(p) != "intel")

# The products that we will ensure are always built
products = [
    ExecutableProduct("alimask", :alimask),
    #ExecutableProduct("easel", :easel),
    ExecutableProduct("esl-afetch", :esl_afetch),
    ExecutableProduct("esl-alimanip", :esl_alimanip),
    ExecutableProduct("esl-alimap", :esl_alimap),
    ExecutableProduct("esl-alimask", :esl_alimask),
    ExecutableProduct("esl-alimerge", :esl_alimerge),
    ExecutableProduct("esl-alipid", :esl_alipid),
    #ExecutableProduct("esl-alirev", :esl_alirev),
    ExecutableProduct("esl-alistat", :esl_alistat),
    ExecutableProduct("esl-compalign", :esl_compalign),
    ExecutableProduct("esl-compstruct", :esl_compstruct),
    ExecutableProduct("esl-construct", :esl_construct),
    ExecutableProduct("esl-histplot", :esl_histplot),
    ExecutableProduct("esl-mask", :esl_mask),
    #ExecutableProduct("esl-mixdchlet", :esl_mixdchlet),
    ExecutableProduct("esl-reformat", :esl_reformat),
    ExecutableProduct("esl-selectn", :esl_selectn),
    ExecutableProduct("esl-seqrange", :esl_seqrange),
    ExecutableProduct("esl-seqstat", :esl_seqstat),
    ExecutableProduct("esl-sfetch", :esl_sfetch),
    ExecutableProduct("esl-shuffle", :esl_shuffle),
    ExecutableProduct("esl-ssdraw", :esl_ssdraw),
    #ExecutableProduct("esl-translate", :esl_translate),
    ExecutableProduct("esl-weight", :esl_weight),
    ExecutableProduct("hmmalign", :hmmalign),
    ExecutableProduct("hmmbuild", :hmmbuild),
    ExecutableProduct("hmmconvert", :hmmconvert),
    ExecutableProduct("hmmemit", :hmmemit),
    ExecutableProduct("hmmfetch", :hmmfetch),
    ExecutableProduct("hmmlogo", :hmmlogo),
    ExecutableProduct("hmmpgmd", :hmmpgmd),
    #ExecutableProduct("hmmpgmd_shard", :hmmpgmd_shard),
    ExecutableProduct("hmmpress", :hmmpress),
    ExecutableProduct("hmmscan", :hmmscan),
    ExecutableProduct("hmmsearch", :hmmsearch),
    ExecutableProduct("hmmsim", :hmmsim),
    ExecutableProduct("hmmstat", :hmmstat),
    ExecutableProduct("jackhmmer", :jackhmmer),
    ExecutableProduct("makehmmerdb", :makehmmerdb),
    ExecutableProduct("nhmmer", :nhmmer),
    ExecutableProduct("nhmmscan", :nhmmscan),
    ExecutableProduct("phmmer", :phmmer),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="GSL_jll", uuid="1b77fbbe-d8ee-58f0-85f9-836ddc23a7a4"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")


text:
# these constants come from and are inspired by the pyprocessing project
# (https://code.google.com/p/pyprocessing/source/browse/trunk/pyprocessing/constants.py)
# with some slight adjustments

# shapes
const POINTS = GL_POINTS
export POINTS
const LINES = GL_LINES
export LINES
const TRIANGLES = GL_TRIANGLES
export TRIANGLES
const TRIANGLE_FAN = GL_TRIANGLE_FAN
export TRIANGLE_FAN
const TRIANGLE_STRIP = GL_TRIANGLE_STRIP
export TRIANGLE_STRIP
const QUADS = GL_QUADS
export QUADS
const QUAD_STRIP = GL_QUAD_STRIP
export QUAD_STRIP
const CLOSE = 1
export CLOSE

#perlin noise table
const perlin = [
151,160,137,91,90,15,131,13,201,95,96,53,194,233,7,225,140,36,103,
30,69,142,8,99,37,240,21,10,23,190,6,148,247,120,234,75,0,26,197,
62,94,252,219,203,117,35,11,32,57,177,33,88,237,149,56,87,174,20,
125,136,171,168,68,175,74,165,71,134,139,48,27,166,77,146,158,231,
83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244,102,
143,54,65,25,63,161,1,216,80,73,209,76,132,187,208,89,18,169,200,
196,135,130,116,188,159,86,164,100,109,198,173,186,3,64,52,217,226,
250,124,123,5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,
58,17,182,189,28,42,223,183,170,213,119,248,152,2,44,154,163,70,
221,153,101,155,167,43,172,9,129,22,39,253,19,98,108,110,79,113,
224,232,178,185,112,104,218,246,97,228,251,34,242,193,238,210,144,
12,191,179,162,241,81,51,145,235,249,14,239,107,49,192,214,31,181,
199,106,157,184,84,204,176,115,121,50,45,127,4,150,254,138,236,
205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180,
151,160,137,91,90,15,131,13,201,95,96,53,194,233,7,225,140,36,103,
30,69,142,8,99,37,240,21,10,23,190,6,148,247,120,234,75,0,26,197,
62,94,252,219,203,117,35,11,32,57,177,33,88,237,149,56,87,174,20,
125,136,171,168,68,175,74,165,71,134,139,48,27,166,77,146,158,231,
83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244,102,
143,54,65,25,63,161,1,216,80,73,209,76,132,187,208,89,18,169,200,
196,135,130,116,188,159,86,164,100,109,198,173,186,3,64,52,217,226,
250,124,123,5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,
58,17,182,189,28,42,223,183,170,213,119,248,152,2,44,154,163,70,
221,153,101,155,167,43,172,9,129,22,39,253,19,98,108,110,79,113,
224,232,178,185,112,104,218,246,97,228,251,34,242,193,238,210,144,
12,191,179,162,241,81,51,145,235,249,14,239,107,49,192,214,31,181,
199,106,157,184,84,204,176,115,121,50,45,127,4,150,254,138,236,
205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180]
export perlin


text:
"""
    @test_forall argument_data, prop

Checks `prop` on each tuples of values specified in `argument_data`.
Returns `Pass` if `prop` holds for all values tested,
`Fail` if `prop` does not hold for some tested value,
and `Error` if a part of code could not be evaluated.

## Arguments

* `argument_data`: a sequence of at least one comma-separated expression(s),
  specifying the sets of values for individual variables.
  Expressions must have the form `x in iter`,
  where `x` is any symbol that will be treated as a dummy variable,
  and `iter` is any iterable expression,
  i.e. one valid in constructing a for loop: `"for x in iter"`.
  In particular, `iter` can be `a:b`,
  where `a` & `b` are integer-valued expressions; an array; or a set.

* `prop` is boolean-valued expression of the proposition to be checked for
  arguments in specified ranges,
  with reference to the dummy variables named in `argument_data`.

## Details

`@test_forall` constructs a nested sequence of for-loops,
and checks that `prop` evaluates to `true`
for all combinations of values of variables specified by the loops.
In set-theoretic terms, `@test_forall` tests
`prop` on the Cartesian product of the sets of variables
specified in `argument_data`.

Loops are generated from `argument_data` inductively,
starting with the outermost loop as the leftmost expression,
to the innermost loop as the rightmost.
Variable ranges may depend on the values of previously referenced
variables (see third example).

## Examples

```julia
julia> @test_forall x in -1:1, x*(x-1)*(x+1) == 0
Test Passed
  Expression: (x in -1:1,x * (x - 1) * (x + 1) == 0)

julia> @test_forall x in -1:1, x*(x-2)*(x+1) == 0
Test Failed
  Expression: (x in -1:1,x * (x - 2) * (x + 1) == 0)

julia> @test_forall x in [0,1,2], y in x:4, y+4 > 2*x
Test Passed
  Expression: (x in [0,1,2],y in x:4,(y + 4 > 2x) == true)

julia> @test_forall x in (1-1, 0+1, sqrt(4)), y in x:4, y+4 > 2*x
Test Passed
  Expression: (x in (1 - 1,0 + 1,sqrt(4)),y in x:4,(y + 4 > 2x) == true)

julia> @test_forall x in ["a","b"], y in ["z","w"], x*y in Set(["az","aw","bz","bw"])
Test Passed
  Expression: (x in ["a","b"],y in ["z","w"],x * y in Set(["az","aw","bz","bw"]))

# Bad design: throws UndefVarError: x not defined
# The outer loop for `y` refers to `x`, which comes later in `argument_data`.
julia> @test_forall y in x:4, x in 0:2, y+4 > 2*x
Error During Test
  Test threw an exception of type UndefVarError
  Expression: (y in x:4,x in 0:2,(y + 4 > 2x) == true)
  UndefVarError: x not defined
```
"""
macro test_forall(factex::Expr)
    ex = begin
        try
            if factex.head != :tuple
                error("Expression of unsupported format for @test_forall: $factex")
            end
            # Defining ex
            iterate = Expr(:block)
            for p in factex.args[1:(end-1)]
                #TODO: Rename these intelligibly
                key = p.args[2]
                val = p.args[3]
                push!(iterate.args,:($(esc(key))=$(esc(val))))
            end
            last_value  = nothing
            prop =
                quote
                    #Why on Earth does this work?
                    #Is there a clearer way to express this?
                    last_value = $(esc(factex.args[end].args[2].args[2]))
                    res = res && $(esc(factex.args[end]))
                    if !res
                        fail_data = $(esc(factex.args[end]))
                        more_data = $(esc(last_value))
                        break
                    end
                end
            Expr(:for,iterate,prop)
        catch err
            quote
                throw($err)
            end
        end
    end
    quote
        res = true
        fail_data = false
        more_data = false
        result =
            try
                $ex
                if res
                    @show 1
                    Pass(:test,$(Expr(:quote, factex)), nothing, nothing)
                else
                    @show 2
                    @show fail_data
                    @show more_data
                    Fail(:test,$(Expr(:quote, factex)), [fail_data, more_data],
                        nothing)
                end
            catch err
                Error(:test_error,$(Expr(:quote, factex)), err, catch_backtrace())
            end
        record(get_testset(), result)
    end
end


text:
using LocalBinaryPatterns: uniform_encoding_table

@testset "uniform_encoding" begin
    lookup = uniform_encoding_table(UInt8, 0:255, 2)
    new_lookup = uniform_encoding_table(UInt8, 0:255, 2)

    # repeat calls are cached in global dictionary
    @test lookup === new_lookup
    @test lookup === uniform_encoding_table(UInt8, 8, 2)

    X_UInt8 = uniform_encoding_table(UInt8, 0:255, 2)
    @test X_UInt8[0b0000_0001] == 0b0000_0001
    @test X_UInt8[0b0000_0101] == 9

    X_UInt8 = uniform_encoding_table(UInt8, 0:255, 4)
    @test X_UInt8[0b0000_0001] == 0b0000_0001
    @test X_UInt8[0b0000_0101] == 0b0000_0101
    @test X_UInt8[0b0000_1101] == 0b0000_1101
    @test X_UInt8[0b0100_1101] == 9

    X_UInt16 = uniform_encoding_table(UInt16, 0:255, 2)
    @test X_UInt16[0b0000_0001] == 0b0000_0001
    @test X_UInt16[0b0000_0101] == 17

    @test_throws ArgumentError uniform_encoding_table(UInt8, 0:65535, 2)
end


text:
using Base.Test
using ThreePhotons

#Also tested in test_structure.jl

function test_Dlms(l, m, mb)
    l1 = l
    l2 = l
    m1 = m
    m2 = m
    m1b = mb
    m2b = mb

    sum = 0
    da = 0.1
    for theta = 0:da:pi
      for phi = 0:da:2*pi
        for gamma = 0:da:2*pi
          sum += sin(theta)* Dlms(l1, m1, m1b, theta, phi, gamma)* conj(Dlms(l2, m2, m2b, theta, phi, gamma)) * da^3
        end
      end
    end

    @test real(sum) ≈ (8 * pi ^ 2) / (2l + 1) atol=0.02
end

test_Dlms(4,2,2)


text:
const typeinfo = [
      ( "i1",  "integer(kind=int8)"  , "i1rand" ),
      ( "i2",  "integer(kind=int16)" , "i2rand" ),
      ( "i4",  "integer(kind=int32)" , "i4rand" ),
      ( "i8",  "integer(kind=int64)" , "i8rand" ),
      ( "r4",  "real(kind=real32)"   , "r4rand" ),
      ( "r8",  "real(kind=real64)"   , "r8rand" ),
      ( "c8",  "complex(kind=real32)", "c8rand" ),
      ( "c16", "complex(kind=real64)", "c16rand" ),
      ( "str", "character(len=*)"    , "strrand" ) ]

function genrandom()
   for (name, ftype, func) in typeinfo
      cmd = `./ftl-expand name=$name type=$ftype func=$func`
      run(pipeline(cmd, stdin="random.ftl", stdout="random$(name).f90"))
   end
end

genrandom()


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

"""
    AbstractScalarShape{T} <: AbstractValueShape
"""
abstract type AbstractScalarShape{T} <: AbstractValueShape end

export AbstractScalarShape

@inline Base.size(::AbstractScalarShape) = ()
@inline Base.length(::AbstractScalarShape) = 1

@inline default_unshaped_eltype(shape::AbstractScalarShape{T}) where {T<:Real} =
    default_datatype(T)

@inline default_unshaped_eltype(shape::AbstractScalarShape{Complex}) = default_datatype(Real)

@inline default_unshaped_eltype(shape::AbstractScalarShape{<:Complex{T}}) where {T} =
    default_unshaped_eltype(_valshapeoftype(T))

@inline shaped_type(shape::AbstractScalarShape{<:Real}, ::Type{T}) where {T<:Real} = T

@inline shaped_type(shape::AbstractScalarShape{<:Complex}, ::Type{T}) where {T<:Real} = Complex{T}

"""
    ScalarShape{T} <: AbstractScalarShape{T}

An `ScalarShape` describes the shape of scalar values of a given type.

Constructor:

    ScalarShape{T::Type}()

T may be an abstract type of Union, or a specific type, e.g.

    ScalarShape{Real}()
    ScalarShape{Integer}()
    ScalarShape{Float32}()
    ScalarShape{Complex}()

Scalar shapes may have a total number of degrees of freedom
(see [`totalndof`](@ref)) greater than one, e.g. shapes of complex-valued
scalars:

    totalndof(ScalarShape{Real}()) == 1
    totalndof(ScalarShape{Complex}()) == 2

See also the documentation of [`AbstractValueShape`](@ref).
"""
struct ScalarShape{T} <: AbstractScalarShape{T} end

export ScalarShape

@inline _valshapeoftype(T::Type{<:Number}) = ScalarShape{T}()

@inline totalndof(::ScalarShape{T}) where {T <: Real} = 1

@inline function totalndof(::ScalarShape{T}) where {T <: Any}
    if @generated
        fieldtypes = ntuple(i -> fieldtype(T, i), Val(fieldcount(T)))
        field_flatlenghts = sum(U -> totalndof(_valshapeoftype(U)), fieldtypes)
        l = prod(field_flatlenghts)
        quote $l end
    else
        fieldtypes = ntuple(i -> fieldtype(T, i), Val(fieldcount(T)))
        field_flatlenghts = sum(U -> totalndof(_valshapeoftype(U)), fieldtypes)
        l = prod(field_flatlenghts)
        l
    end
end

(shape::ScalarShape{T})(::UndefInitializer) where {T<:Number} = zero(default_datatype(T))

@static if VERSION < v"1.3"
    # Workaround for Julia issue #14919
    @inline (shape::ScalarShape)(data::AbstractVector{<:Real}) =
        _apply_shape_to_data(shape, data)
end

const ScalarAccessor{T} = ValueAccessor{ScalarShape{T}} where {T}

@inline view_idxs(idxs::AbstractUnitRange{<:Integer}, va::ScalarAccessor{<:Real}) = first(idxs) + va.offset

# ToDo: view_idxs for scalars with dof greater than 1 (complex, etc.)

Base.@propagate_inbounds vs_getindex(data::AbstractVector{<:Real}, va::ScalarAccessor) =
    data[view_idxs(axes(data, 1), va)]

Base.@propagate_inbounds vs_unsafe_view(data::AbstractVector{<:Real}, va::ScalarAccessor) =
    Base.unsafe_view(data, view_idxs(axes(data, 1), va))

Base.@propagate_inbounds vs_setindex!(data::AbstractVector{<:Real}, v, va::ScalarAccessor) =
    setindex!(data, v, view_idxs(axes(data, 1), va))

Base.@propagate_inbounds function _bcasted_view(data::AbstractVectorOfSimilarVectors{<:Real}, va::ScalarAccessor)
    flat_data = flatview(data)
    idxs = view_idxs(axes(flat_data, 1), va)
    view(flat_data, idxs, :)
end

Base.copy(instance::VSBroadcasted2{typeof(getindex), AbstractVectorOfSimilarVectors{<:Real},Ref{<:ScalarAccessor}}) =
    copy(_bcasted_view(instance.args[1], instance.args[2][]))

Base.copy(instance::VSBroadcasted2{typeof(view), AbstractVectorOfSimilarVectors{<:Real},Ref{<:ScalarAccessor}}) =
    _bcasted_view(instance.args[1], instance.args[2][])


text:
export Operator
export bandwidths, bandrange, \, periodic
export neumann
export ldirichlet,rdirichlet,lneumann,rneumann
export ldiffbc,rdiffbc,diffbcs
export domainspace,rangespace

abstract type Operator{T} end #T is the entry type, Float64 or Complex{Float64}

eltype(::Operator{T}) where {T} = T
eltype(::Type{Operator{T}}) where {T} = T
eltype(::Type{OT}) where {OT<:Operator} = eltype(supertype(OT))

# default entry type
# we assume entries depend on both the domain and the basis
# realdomain case doesn't use

prectype(sp::Space) = promote_type(prectype(domaintype(sp)),eltype(rangetype(sp)))

 #Operators are struct
copy(A::Operator) = A

BroadcastStyle(::Type{<:Operator}) = DefaultArrayStyle{2}()
broadcastable(A::Operator) = A

## We assume operators are T->T
rangespace(A::Operator) = error("Override rangespace for $(typeof(A))")
domainspace(A::Operator) = error("Override domainspace for $(typeof(A))")
spaces(A::Operator) = (rangespace(A), domainspace(A)) # order is consistent with size(::Matrix)
domain(A::Operator) = domain(domainspace(A))

isconstspace(_) = false
## Functionals
isafunctional(A::Operator) = size(A,1)==1 && isconstspace(rangespace(A))

isonesvec(A) = A isa AbstractFill && getindex_value(A) == 1
# block lengths of a space are 1
hastrivialblocks(A::Space) = isonesvec(blocklengths(A))
hastrivialblocks(A::Operator) = hastrivialblocks(domainspace(A)) &&
                                hastrivialblocks(rangespace(A))

# blocklengths are constant lengths
hasconstblocks(A::Space) = isa(blocklengths(A),AbstractFill)
hasconstblocks(A::Operator) = hasconstblocks(domainspace(A)) && hasconstblocks(rangespace(A)) &&
                                getindex_value(blocklengths(domainspace(A))) == getindex_value(blocklengths(rangespace(A)))

macro functional(FF)
    quote
        Base.size(A::$FF,k::Integer) = k==1 ? 1 : dimension(domainspace(A))
        ApproxFunBase.rangespace(F::$FF) = ConstantSpace(eltype(F))
        ApproxFunBase.isafunctional(::$FF) = true
        ApproxFunBase.blockbandwidths(A::$FF) = 0,hastrivialblocks(domainspace(A)) ? bandwidth(A,2) : ∞
        function ApproxFunBase.defaultgetindex(f::$FF,k::Integer,j::Integer)
            @assert k==1
            f[j]::eltype(f)
        end
        function ApproxFunBase.defaultgetindex(f::$FF,k::Integer,j::AbstractRange)
            @assert k==1
            f[j]
        end
        function ApproxFunBase.defaultgetindex(f::$FF,k::Integer,j)
            @assert k==1
            f[j]
        end
        function ApproxFunBase.defaultgetindex(f::$FF,k::AbstractRange,j::Integer)
            @assert k==1:1
            f[j]
        end
        function ApproxFunBase.defaultgetindex(f::$FF,k::AbstractRange,j::AbstractRange)
            @assert k==1:1
            reshape(f[j],1,length(j))
        end
        function ApproxFunBase.defaultgetindex(f::$FF,k::AbstractRange,j)
            @assert k==1:1
            reshape(f[j],1,length(j))
        end
    end
end

blocksize(A::Operator,k) = k==1 ? length(blocklengths(rangespace(A))) : length(blocklengths(domainspace(A)))
blocksize(A::Operator) = (blocksize(A,1),blocksize(A,2))

Base.size(A::Operator) = (size(A,1),size(A,2))
Base.size(A::Operator,k::Integer) = k==1 ? dimension(rangespace(A)) : dimension(domainspace(A))
Base.length(A::Operator) = size(A,1) * size(A,2)

# used to compute "end" for last index
function lastindex(A::Operator, n::Integer)
    if n > 2
        1
    elseif n==2
        size(A,2)
    elseif isinf(size(A,2)) || isinf(size(A,1))
        ∞
    else
        size(A,1)
    end
end
lastindex(A::Operator) = size(A,1)*size(A,2)

Base.ndims(::Operator) = 2

## bandrange and indexrange
isbandedbelow(A::Operator) = isfinite(bandwidth(A,1))
isbandedabove(A::Operator) = isfinite(bandwidth(A,2))
isbanded(A::Operator) = isbandedbelow(A) && isbandedabove(A)

isbandedblockbandedbelow(_) = false
isbandedblockbandedabove(_) = false

isbandedblockbanded(A::Operator) = isbandedblockbandedabove(A) && isbandedblockbandedbelow(A)

# this should be determinable at compile time
#TODO: I think it can be generalized to the case when the domainspace
# blocklengths == rangespace blocklengths, in which case replace the definition
# of p with maximum(blocklength(domainspace(A)))
function blockbandwidths(A::Operator)
    hastrivialblocks(A) && return bandwidths(A)

    if hasconstblocks(A)
        a,b = bandwidths(A)
        p = getindex_value(blocklengths(domainspace(A)))
        return (-fld(-a,p),-fld(-b,p))
    end

    #TODO: Generalize to finite dimensional
    if size(A,2) == 1
        rs = rangespace(A)

        if hasconstblocks(rs)
            a = bandwidth(A,1)
            p = getindex_value(blocklengths(rs))
            return (-fld(-a,p),0)
        end
    end

    return (length(blocklengths(rangespace(A)))-1,length(blocklengths(domainspace(A)))-1)
end

# assume dense blocks
subblockbandwidths(K::Operator) = maximum(blocklengths(rangespace(K)))-1, maximum(blocklengths(domainspace(K)))-1

isblockbandedbelow(A) = isfinite(blockbandwidth(A,1))
isblockbandedabove(A) = isfinite(blockbandwidth(A,2))
isblockbanded(A::Operator) = isblockbandedbelow(A) && isblockbandedabove(A)

israggedbelow(A::Operator) = isbandedbelow(A) || isbandedblockbanded(A) || isblockbandedbelow(A)

blockbandwidth(K::Operator, k::Integer) = blockbandwidths(K)[k]
subblockbandwidth(K::Operator,k::Integer) = subblockbandwidths(K)[k]

bandwidth(A::Operator, k::Integer) = bandwidths(A)[k]
# we are always banded by the size
bandwidths(A::Operator) = (size(A,1)-1,size(A,2)-1)
bandwidths(A::Operator, k::Integer) = bandwidths(A)[k]

## Strides
# lets us know if operators decouple the entries
# to split into sub problems
# A diagonal operator has essentially infinite stride
# which we represent by a factorial, so that
# the gcd with any number < 10 is the number
stride(A::Operator) =
    isdiag(A) ? factorial(10) : 1

isdiag(A::Operator) = bandwidths(A)==(0,0)
istriu(A::Operator) = bandwidth(A, 1) == 0
istril(A::Operator) = bandwidth(A, 2) == 0

## Construct operators

include("SubOperator.jl")

#
# sparse(B::Operator,n::Integer)=sparse(BandedMatrix(B,n))
# sparse(B::Operator,n::AbstractRange,m::AbstractRange)=sparse(BandedMatrix(B,n,m))
# sparse(B::Operator,n::Colon,m::AbstractRange)=sparse(BandedMatrix(B,n,m))
# sparse(B::Operator,n::AbstractRange,m::Colon)=sparse(BandedMatrix(B,n,m))

## geteindex

getindex(B::Operator,k,j) = defaultgetindex(B,k,j)
getindex(B::Operator,k) = defaultgetindex(B,k)
getindex(B::Operator,k::Block{2}) = B[Block.(k.n)...]

## override getindex.

defaultgetindex(B::Operator,k::Integer) = error("Override [k] for $(typeof(B))")
defaultgetindex(B::Operator,k::Integer,j::Integer) = error("Override [k,j] for $(typeof(B))")

# Ranges

defaultgetindex(op::Operator,kr::AbstractRange) = eltype(op)[op[k] for k in kr]
defaultgetindex(B::Operator,k::Block,j::Block) = AbstractMatrix(view(B,k,j))
defaultgetindex(B::Operator,k::AbstractRange,j::Block) = AbstractMatrix(view(B,k,j))
defaultgetindex(B::Operator,k::Block,j::AbstractRange) = AbstractMatrix(view(B,k,j))
defaultgetindex(B::Operator,k::AbstractRange,j::AbstractRange) = AbstractMatrix(view(B,k,j))

defaultgetindex(op::Operator,k::Integer,jr::AbstractRange) = eltype(op)[op[k,j] for j in jr]
defaultgetindex(op::Operator,kr::AbstractRange,j::Integer) = eltype(op)[op[k,j] for k in kr]

defaultgetindex(B::Operator,k::Block,j::BlockRange) = AbstractMatrix(view(B,k,j))
defaultgetindex(B::Operator,k::BlockRange,j::BlockRange) = AbstractMatrix(view(B,k,j))

defaultgetindex(op::Operator,k::Integer,jr::BlockRange) = eltype(op)[op[k,j] for j in jr]
defaultgetindex(op::Operator,kr::BlockRange,j::Integer) = eltype(op)[op[k,j] for k in kr]

# Colon casdes
defaultgetindex(A::Operator,kj::CartesianIndex{2}) = A[kj[1],kj[2]]
defaultgetindex(A::Operator,kj::CartesianIndex{1}) = A[kj[1]]
defaultgetindex(A::Operator,k,j) = view(A,k,j)

# TODO: finite dimensional blocks
blockcolstart(A::Operator,J::Integer) = Block(max(1,J-blockbandwidth(A,2)))
blockrowstart(A::Operator,K::Integer) = Block(max(1,K-blockbandwidth(A,1)))
blockcolstop(A::Operator,J::Integer) = Block(min(J+blockbandwidth(A,1),blocksize(A,1)))
blockrowstop(A::Operator,K::Integer) = Block(min(K+blockbandwidth(A,2),blocksize(A,2)))

blockrows(A::Operator,K::Integer) = blockrange(rangespace(A),K)
blockcols(A::Operator,J::Integer) = blockrange(domainspace(A),J)

# default is to use bandwidth
# override for other shaped operators
#TODO: Why size(A,2) in colstart?
banded_colstart(A::Operator, i::Integer) = min(max(i-bandwidth(A,2), 1), size(A, 2))
banded_colstop(A::Operator, i::Integer) = max(0,min(i+bandwidth(A,1), size(A, 1)))
banded_rowstart(A::Operator, i::Integer) = min(max(i-bandwidth(A,1), 1), size(A, 1))
banded_rowstop(A::Operator, i::Integer) = max(0,min(i+bandwidth(A,2), size(A, 2)))

blockbanded_colstart(A::Operator, i::Integer) =
        blockstart(rangespace(A), block(domainspace(A),i)-blockbandwidth(A,2))
blockbanded_colstop(A::Operator, i::Integer) =
    min(blockstop(rangespace(A), block(domainspace(A),i)+blockbandwidth(A,1)),
        size(A, 1))
blockbanded_rowstart(A::Operator, i::Integer) =
        blockstart(domainspace(A), block(rangespace(A),i)-blockbandwidth(A,1))
blockbanded_rowstop(A::Operator, i::Integer) =
    min(blockstop(domainspace(A), block(rangespace(A),i)+blockbandwidth(A,2)),
        size(A, 2))

function bandedblockbanded_colstart(A::Operator, i::Integer)
    ds = domainspace(A)
    B = block(ds,i)
    ξ = i - blockstart(ds,B) + 1  # col in block
    bs = blockstart(rangespace(A), B-blockbandwidth(A,2))
    max(bs,bs + ξ - 1 - subblockbandwidth(A,2))
end

function bandedblockbanded_colstop(A::Operator, i::Integer)
    i ≤ 0 && return 0
    ds = domainspace(A)
    rs = rangespace(A)
    B = block(ds,i)
    ξ = i - blockstart(ds,B) + 1  # col in block
    Bend = B+blockbandwidth(A,1)
    bs = blockstart(rs, Bend)
    min(blockstop(rs,Bend),bs + ξ - 1 + subblockbandwidth(A,1))
end

function bandedblockbanded_rowstart(A::Operator, i::Integer)
    rs = rangespace(A)
    B = block(rs,i)
    ξ = i - blockstart(rs,B) + 1  # row in block
    bs = blockstart(domainspace(A), B-blockbandwidth(A,1))
    max(bs,bs + ξ - 1 - subblockbandwidth(A,1))
end

function bandedblockbanded_rowstop(A::Operator, i::Integer)
    ds = domainspace(A)
    rs = rangespace(A)
    B = block(rs,i)
    ξ = i - blockstart(rs,B) + 1  # row in block
    Bend = B+blockbandwidth(A,2)
    bs = blockstart(ds, Bend)
    min(blockstop(ds,Bend),bs + ξ - 1 + subblockbandwidth(A,2))
end

unstructured_colstart(A, i) = 1
unstructured_colstop(A, i) = size(A,1)
unstructured_rowstart(A, i) = 1
unstructured_rowstop(A, i) = size(A,2)

function default_colstart(A::Operator, i::Integer)
    if isbandedabove(A)
        banded_colstart(A,i)
    elseif isbandedblockbanded(A)
        bandedblockbanded_colstart(A, i)
    elseif isblockbanded(A)
        blockbanded_colstart(A, i)
    else
        unstructured_colstart(A, i)
    end
end

function default_colstop(A::Operator, i::Integer)
    if isbandedbelow(A)
        banded_colstop(A,i)
    elseif isbandedblockbanded(A)
        bandedblockbanded_colstop(A, i)
    elseif isblockbanded(A)
        blockbanded_colstop(A, i)
    else
        unstructured_colstop(A, i)
    end
end

function default_rowstart(A::Operator, i::Integer)
    if isbandedbelow(A)
        banded_rowstart(A,i)
    elseif isbandedblockbanded(A)
        bandedblockbanded_rowstart(A, i)
    elseif isblockbanded(A)
        blockbanded_rowstart(A, i)
    else
        unstructured_rowstart(A, i)
    end
end

function default_rowstop(A::Operator, i::Integer)
    if isbandedabove(A)
        banded_rowstop(A,i)
    elseif isbandedblockbanded(A)
        bandedblockbanded_rowstop(A, i)
    elseif isblockbanded(A)
        blockbanded_rowstop(A, i)
    else
        unstructured_rowstop(A, i)
    end
end

for OP in (:colstart,:colstop,:rowstart,:rowstop)
    defOP = Meta.parse("default_"*string(OP))
    @eval begin
        $OP(A::Operator,i::Integer) = $defOP(A,i)
        $OP(A::Operator,i::Infinity) = ∞
    end
end

function defaultgetindex(A::Operator,::Type{FiniteRange},::Type{FiniteRange})
    if isfinite(size(A,1)) && isfinite(size(A,2))
        A[1:size(A,1),1:size(A,2)]
    else
        error("Only exists for finite operators.")
    end
end

defaultgetindex(A::Operator,k::Type{FiniteRange},J::Block) = A[k,blockcols(A,J)]
function defaultgetindex(A::Operator,::Type{FiniteRange},jr::AbstractVector{Int})
    cs = (isbanded(A) || isblockbandedbelow(A)) ? colstop(A,maximum(jr)) : mapreduce(j->colstop(A,j),max,jr)
    A[1:cs,jr]
end

function defaultgetindex(A::Operator,::Type{FiniteRange},jr::BlockRange{1})
    cs = (isbanded(A) || isblockbandedbelow(A)) ? blockcolstop(A,maximum(jr)) : mapreduce(j->blockcolstop(A,j),max,jr)
    A[Block(1):cs,jr]
end

function view(A::Operator,::Type{FiniteRange},jr::AbstractVector{Int})
    cs = (isbanded(A) || isblockbandedbelow(A)) ? colstop(A,maximum(jr)) : mapreduce(j->colstop(A,j),max,jr)
    view(A,1:cs,jr)
end

function view(A::Operator,::Type{FiniteRange},jr::BlockRange{1})
    cs = (isbanded(A) || isblockbandedbelow(A)) ? blockcolstop(A,maximum(jr)) : mapreduce(j->blockcolstop(A,j),max,jr)
    view(A,Block(1):cs,jr)
end

defaultgetindex(A::Operator,K::Block,j::Type{FiniteRange}) = A[blockrows(A,K),j]
defaultgetindex(A::Operator,kr,::Type{FiniteRange}) =
    A[kr,1:rowstop(A,maximum(kr))]

## Composition with a Fun, LowRankFun, and ProductFun

getindex(B::Operator,f::Fun) = B*Multiplication(domainspace(B),f)
getindex(B::Operator,f::LowRankFun{S,M,SS,T}) where {S,M,SS,T} = mapreduce(i->f.A[i]*B[f.B[i]],+,1:rank(f))
getindex(B::Operator{BT},f::ProductFun{S,V,SS,T}) where {BT,S,V,SS,T} =
    mapreduce(i->f.coefficients[i]*B[Fun(f.space[2],[zeros(promote_type(BT,T),i-1);
                                            one(promote_type(BT,T))])],
                +,1:length(f.coefficients))

# Convenience for wrapper ops
unwrap_axpy!(α,P,A) = BLAS.axpy!(α,view(parent(P).op,P.indexes[1],P.indexes[2]),A)
iswrapper(_) = false
haswrapperstructure(_) = false

# use this for wrapper operators that have the same structure but
# not necessarily the same entries
#
#  Ex: c*op or real(op)
macro wrapperstructure(Wrap)
    ret = quote
        haswrapperstructure(::$Wrap) = true
    end

    for func in (:(ApproxFunBase.bandwidths),:(LinearAlgebra.stride),
                 :(ApproxFunBase.isbandedblockbanded),:(ApproxFunBase.isblockbanded),
                 :(ApproxFunBase.israggedbelow),:(Base.size),:(ApproxFunBase.isbanded),
                 :(ApproxFunBase.blockbandwidths),:(ApproxFunBase.subblockbandwidths),
                 :(LinearAlgebra.issymmetric))
        ret = quote
            $ret

            $func(D::$Wrap) = $func(D.op)
        end
    end

     for func in (:(ApproxFunBase.bandwidth),:(ApproxFunBase.colstart),:(ApproxFunBase.colstop),
                     :(ApproxFunBase.rowstart),:(ApproxFunBase.rowstop),:(ApproxFunBase.blockbandwidth),
                     :(Base.size),:(ApproxFunBase.subblockbandwidth))
         ret = quote
             $ret

             $func(D::$Wrap,k::Integer) = $func(D.op,k)
             $func(A::$Wrap,i::ApproxFunBase.Infinity) = $func(D.op,k)
         end
     end

    esc(ret)
end

# use this for wrapper operators that have the same entries but
# not necessarily the same spaces
#
macro wrappergetindex(Wrap)
    ret = quote
        Base.getindex(OP::$Wrap,k::Integer...) =
            OP.op[k...]::eltype(OP)

        Base.getindex(OP::$Wrap,k::Union{Number,AbstractArray,Colon}...) = OP.op[k...]
        Base.getindex(OP::$Wrap,k::ApproxFunBase.InfRanges, j::ApproxFunBase.InfRanges) = view(OP, k, j)
        Base.getindex(OP::$Wrap,k::ApproxFunBase.InfRanges, j::Colon) = view(OP, k, j)
        Base.getindex(OP::$Wrap,k::Colon, j::ApproxFunBase.InfRanges) = view(OP, k, j)
        Base.getindex(OP::$Wrap,k::Colon, j::Colon) = view(OP, k, j)

        BLAS.axpy!(α,P::ApproxFunBase.SubOperator{T,OP},A::AbstractMatrix) where {T,OP<:$Wrap} =
            ApproxFunBase.unwrap_axpy!(α,P,A)

        ApproxFunBase.mul_coefficients(A::$Wrap,b) = ApproxFunBase.mul_coefficients(A.op,b)
        ApproxFunBase.mul_coefficients(A::ApproxFunBase.SubOperator{T,OP,Tuple{UnitRange{Int},UnitRange{Int}}},b) where {T,OP<:$Wrap} =
            ApproxFunBase.mul_coefficients(view(parent(A).op,S.indexes[1],S.indexes[2]),b)
        ApproxFunBase.mul_coefficients(A::ApproxFunBase.SubOperator{T,OP},b) where {T,OP<:$Wrap} =
            ApproxFunBase.mul_coefficients(view(parent(A).op,S.indexes[1],S.indexes[2]),b)
    end

    for TYP in (:(ApproxFunBase.BandedMatrix),:(ApproxFunBase.RaggedMatrix),
                :Matrix,:Vector,:AbstractVector)
        ret = quote
            $ret

            $TYP(P::ApproxFunBase.SubOperator{T,OP}) where {T,OP<:$Wrap} =
                $TYP(view(parent(P).op,P.indexes[1],P.indexes[2]))
            $TYP(P::ApproxFunBase.SubOperator{T,OP,NTuple{2,UnitRange{Int}}}) where {T,OP<:$Wrap} =
                $TYP(view(parent(P).op,P.indexes[1],P.indexes[2]))
        end
    end

    ret = quote
        $ret

        # fast converts to banded matrices would be based on indices, not blocks
        function ApproxFunBase.BandedMatrix(S::ApproxFunBase.SubOperator{T,OP,NTuple{2,ApproxFunBase.BlockRange1}}) where {T,OP<:$Wrap}
            A = parent(S)
            ds = domainspace(A)
            rs = rangespace(A)
            KR,JR = parentindices(S)
            ApproxFunBase.BandedMatrix(view(A,
                              ApproxFunBase.blockstart(rs,first(KR)):ApproxFunBase.blockstop(rs,last(KR)),
                              ApproxFunBase.blockstart(ds,first(JR)):ApproxFunBase.blockstop(ds,last(JR))))
        end

        # if the spaces change, then we need to be smarter
        function ApproxFunBase.BlockBandedMatrix(S::ApproxFunBase.SubOperator{T,OP}) where {T,OP<:$Wrap}
            P = parent(S)
            if ApproxFunBase.blocklengths(domainspace(P)) === ApproxFunBase.blocklengths(domainspace(P.op)) &&
                    ApproxFunBase.blocklengths(rangespace(P)) === ApproxFunBase.blocklengths(rangespace(P.op))
                ApproxFunBase.BlockBandedMatrix(view(parent(S).op,S.indexes[1],S.indexes[2]))
            else
                ApproxFunBase.default_BlockBandedMatrix(S)
            end
        end

        function ApproxFunBase.PseudoBlockMatrix(S::ApproxFunBase.SubOperator{T,OP}) where {T,OP<:$Wrap}
            P = parent(S)
            if ApproxFunBase.blocklengths(domainspace(P)) === ApproxFunBase.blocklengths(domainspace(P.op)) &&
                    ApproxFunBase.blocklengths(rangespace(P)) === ApproxFunBase.blocklengths(rangespace(P.op))
                ApproxFunBase.PseudoBlockMatrix(view(parent(S).op,S.indexes[1],S.indexes[2]))
            else
                ApproxFunBase.default_blockmatrix(S)
            end
        end

        function ApproxFunBase.BandedBlockBandedMatrix(S::ApproxFunBase.SubOperator{T,OP}) where {T,OP<:$Wrap}
            P = parent(S)
            if ApproxFunBase.blocklengths(domainspace(P)) === ApproxFunBase.blocklengths(domainspace(P.op)) &&
                    ApproxFunBase.blocklengths(rangespace(P)) === ApproxFunBase.blocklengths(rangespace(P.op))
                ApproxFunBase.BandedBlockBandedMatrix(view(parent(S).op,S.indexes[1],S.indexes[2]))
            else
                ApproxFunBase.default_BandedBlockBandedMatrix(S)
            end
        end

        ApproxFunBase.@wrapperstructure($Wrap) # structure is automatically inherited
    end

    esc(ret)
end

# use this for wrapper operators that have the same spaces but
# not necessarily the same entries or structure
#
macro wrapperspaces(Wrap)
    ret = quote  end

    for func in (:(ApproxFunBase.rangespace),:(ApproxFunBase.domain),
                 :(ApproxFunBase.domainspace),:(ApproxFunBase.isconstop))
        ret = quote
            $ret

            $func(D::$Wrap) = $func(D.op)
        end
    end

    esc(ret)
end

# use this for wrapper operators that have the same entries and same spaces
#
macro wrapper(Wrap)
    ret = quote
        ApproxFunBase.@wrappergetindex($Wrap)
        ApproxFunBase.@wrapperspaces($Wrap)

        ApproxFunBase.iswrapper(::$Wrap) = true
    end

    esc(ret)
end

## Standard Operators and linear algebra

include("ldiv.jl")

include("spacepromotion.jl")
include("banded/banded.jl")
include("general/general.jl")

include("functionals/functionals.jl")
include("almostbanded/almostbanded.jl")

include("systems.jl")

include("qr.jl")
include("nullspace.jl")

## Conversion

zero(::Type{Operator{T}}) where {T<:Number} = ZeroOperator(T)
zero(::Type{O}) where {O<:Operator} = ZeroOperator(eltype(O))

Operator(L::UniformScaling) = ConstantOperator(L, UnsetSpace())
Operator(L::UniformScaling, s::Space) = ConstantOperator(L, s)
Operator(L::UniformScaling{Bool}, s::Space) = L.λ ? IdentityOperator(s) : ZeroOperator(s)
Operator(L::UniformScaling, d::Domain) = Operator(L, Space(d))

Operator{T}(f::Fun) where {T} =
    norm(f.coefficients)==0 ? zero(Operator{T}) : convert(Operator{T}, Multiplication(f))

Operator(f::Fun) = norm(f.coefficients)==0 ? ZeroOperator() : Multiplication(f)

convert(::Type{O}, f::Fun) where O<:Operator = O(f)
Operator{T}(A::Operator) where T = convert(Operator{T}, A)

## Promotion

promote_rule(::Type{N},::Type{Operator}) where {N<:Number} = Operator{N}
promote_rule(::Type{UniformScaling{N}},::Type{Operator}) where {N<:Number} =
    Operator{N}
promote_rule(::Type{Fun{S,N,VN}},::Type{Operator}) where {S,N<:Number,VN} = Operator{N}
promote_rule(::Type{N},::Type{O}) where {N<:Number,O<:Operator} =
    Operator{promote_type(N,eltype(O))}  # float because numbers are promoted to Fun
promote_rule(::Type{UniformScaling{N}},::Type{O}) where {N<:Number,O<:Operator} =
    Operator{promote_type(N,eltype(O))}
promote_rule(::Type{Fun{S,N,VN}},::Type{O}) where {S,N<:Number,O<:Operator,VN} =
    Operator{promote_type(N,eltype(O))}

promote_rule(::Type{BO1},::Type{BO2}) where {BO1<:Operator,BO2<:Operator} =
    Operator{promote_type(eltype(BO1),eltype(BO2))}

## Wrapper

#TODO: Should cases that modify be included?
const WrapperOperator = Union{SpaceOperator,MultiplicationWrapper,DerivativeWrapper,IntegralWrapper,
                                    ConversionWrapper,ConstantTimesOperator,TransposeOperator}

# The following support converting an Operator to a Matrix or BandedMatrix

## BLAS and matrix routines
# We assume that copy may be overriden

BLAS.axpy!(a, X::Operator, Y::AbstractMatrix) = (Y .= a .* AbstractMatrix(X) .+ Y)
copyto!(dest::AbstractMatrix, src::Operator) = copyto!(dest, AbstractMatrix(src))

# this is for operators that implement copy via axpy!

BandedMatrix(::Type{Zeros}, V::Operator) = BandedMatrix(Zeros{eltype(V)}(size(V)), bandwidths(V))
Matrix(::Type{Zeros}, V::Operator) = Matrix(Zeros{eltype(V)}(size(V)))
BandedBlockBandedMatrix(::Type{Zeros}, V::Operator) =
    BandedBlockBandedMatrix(Zeros{eltype(V)}(size(V)),
                            blocklengths(rangespace(V)), blocklengths(domainspace(V)),
                            blockbandwidths(V), subblockbandwidths(V))
BlockBandedMatrix(::Type{Zeros}, V::Operator) =
    BlockBandedMatrix(Zeros{eltype(V)}(size(V)),
                      AbstractVector{Int}(blocklengths(rangespace(V))),
                       AbstractVector{Int}(blocklengths(domainspace(V))),
                      blockbandwidths(V))
RaggedMatrix(::Type{Zeros}, V::Operator) =
    RaggedMatrix(Zeros{eltype(V)}(size(V)),
                 Int[max(0,colstop(V,j)) for j=1:size(V,2)])

convert_axpy!(::Type{MT}, S::Operator) where {MT <: AbstractMatrix} =
        BLAS.axpy!(one(eltype(S)), S, MT(Zeros, S))

BandedMatrix(S::Operator) = default_BandedMatrix(S)

function BlockBandedMatrix(S::Operator)
    if isbandedblockbanded(S)
        BlockBandedMatrix(BandedBlockBandedMatrix(S))
    else
        default_BlockBandedMatrix(S)
    end
end

function default_BlockMatrix(S::Operator)
    ret = PseudoBlockArray(zeros(size(S)),
                        AbstractVector{Int}(blocklengths(rangespace(S))),
                        AbstractVector{Int}(blocklengths(domainspace(S))))
    ret .= S
    ret
end

function PseudoBlockMatrix(S::Operator)
    if isbandedblockbanded(S)
        PseudoBlockMatrix(BandedBlockBandedMatrix(S))
    elseif isblockbanded(S)
        PseudoBlockMatrix(BlockBandedMatrix(S))
    else
        default_BlockMatrix(S)
    end
end

# TODO: Unify with SubOperator
for TYP in (:RaggedMatrix, :Matrix)
    def_TYP = Meta.parse("default_" * string(TYP))
    @eval function $TYP(S::Operator)
        if isinf(size(S,1)) || isinf(size(S,2))
            error("Cannot convert $S to a $TYP")
        end

        if isbanded(S)
            $TYP(BandedMatrix(S))
        else
            $def_TYP(S)
        end
    end
end

function Vector(S::Operator)
    if size(S,2) ≠ 1  || isinf(size(S,1))
        error("Cannot convert $S to a AbstractVector")
    end

    eltype(S)[S[k] for k=1:size(S,1)]
end

convert(::Type{AA}, B::Operator) where AA<:AbstractArray = AA(B)

# TODO: template out fully
arraytype(::Operator) = Matrix
function arraytype(V::SubOperator{T,B,Tuple{KR,JR}}) where {T, B, KR <: Union{BlockRange, Block}, JR <: Union{BlockRange, Block}}
    P = parent(V)
    isbandedblockbanded(P) && return BandedBlockBandedMatrix
    isblockbanded(P) && return BlockBandedMatrix
    return PseudoBlockMatrix
end

function arraytype(V::SubOperator{T,B,Tuple{KR,JR}}) where {T, B, KR <: Block, JR <: Block}
    P = parent(V)
    isbandedblockbanded(V) && return BandedMatrix
    return Matrix
end

function arraytype(V::SubOperator)
    P = parent(V)
    isbanded(P) && return BandedMatrix
    # isbandedblockbanded(P) && return BandedBlockBandedMatrix
    isinf(size(P,1)) && israggedbelow(P) && return RaggedMatrix
    return Matrix
end

AbstractMatrix(V::Operator) = arraytype(V)(V)
AbstractVector(S::Operator) = Vector(S)

# default copy is to loop through
# override this for most operators.
function default_BandedMatrix(S::Operator)
    Y=BandedMatrix{eltype(S)}(undef, size(S), bandwidths(S))

    for j=1:size(S,2),k=colrange(Y,j)
        @inbounds inbands_setindex!(Y,S[k,j],k,j)
    end

    Y
end

# default copy is to loop through
# override this for most operators.
function default_RaggedMatrix(S::Operator)
    data=Array{eltype(S)}(undef, 0)
    cols=Array{Int}(undef, size(S,2)+1)
    cols[1]=1
    for j=1:size(S,2)
        cs=colstop(S,j)
        K=cols[j]-1
        cols[j+1]=cs+cols[j]
        resize!(data,cols[j+1]-1)

        for k=1:cs
            data[K+k]=S[k,j]
        end
    end

    RaggedMatrix(data,cols,size(S,1))
end

function default_Matrix(S::Operator)
    n, m = size(S)
    if isinf(n) || isinf(m)
        error("Cannot convert $S to a Matrix")
    end

    eltype(S)[S[k,j] for k=1:n, j=1:m]
end

# The diagonal of the operator may not be the diagonal of the sub
# banded matrix, so the following calculates the row of the
# Banded matrix corresponding to the diagonal of the original operator

diagindshift(S,kr,jr) = first(kr)-first(jr)
diagindshift(S::SubOperator) = diagindshift(S,parentindices(S)[1],parentindices(S)[2])

#TODO: Remove
diagindrow(S,kr,jr) = bandwidth(S,2)+first(jr)-first(kr)+1
diagindrow(S::SubOperator) = diagindrow(S,parentindices(S)[1],parentindices(S)[2])


text:
include("Mustache_test.jl")
include("multiple_nested_sections.jl")
include("test_index.jl")
#include("data-frame-test.jl")


text:
function part1(input)
    instructions = parse_line.(eachline(input))
    grid = falses(1000, 1000)
    for (action, x1, y1, x2, y2) in instructions
        if action == :on
            grid[y1:y2, x1:x2] .= true
        elseif action == :off
            grid[y1:y2, x1:x2] .= false
        else
            grid[y1:y2, x1:x2] .= .!grid[y1:y2, x1:x2]
        end
    end
    return count(grid)
end

function part2(input)
    instructions = parse_line.(eachline(input))
    grid = fill(0, 1000, 1000)
    for (action, x1, y1, x2, y2) in instructions
        if action == :on
            grid[y1:y2, x1:x2] .+= 1
        elseif action == :off
            grid[y1:y2, x1:x2] .= max.(0, grid[y1:y2, x1:x2] .- 1)
        else
            grid[y1:y2, x1:x2] .+= 2
        end
    end
    return sum(grid)
end

function parse_line(line)
    parts = split(line)
    x1, y1 = parse.(Int, split(parts[end - 2], ",")) .+ 1
    x2, y2 = parse.(Int, split(parts[end], ",")) .+ 1
    return Symbol(parts[end - 3]), x1, y1, x2, y2
end


text:
using Libdl

# Define a few things that the Clang wrapper doesn't quite take care of for us
function rcutils_get_default_allocator end
function rcutils_reallocf end
const SIZE_MAX = ~Csize_t(0)
# const __WCHAR_MIN = Cwchar_t <: Unsigned ? Cwchar_t(0) : ~Cwchar_t(0)
# const __WCHAR_MAX = Cwchar_t <: Unsigned ? ~Cwchar_t(0) : Cwchar_t(~0) >>> 1
RCUTILS_S_TO_NS( seconds ) = ( seconds * ( 1000 * 1000 * 1000 ) ) # Convenience macro to convert milliseconds to nanoseconds.
RCUTILS_MS_TO_NS( milliseconds ) = ( milliseconds * ( 1000 * 1000 ) ) # Convenience macro to convert microseconds to nanoseconds.
RCUTILS_US_TO_NS( microseconds ) = ( microseconds * 1000 ) # Convenience macro to convert nanoseconds to seconds.
RCUTILS_NS_TO_S( nanoseconds ) = ( nanoseconds ÷ ( 1000 * 1000 * 1000 ) ) # Convenience macro to convert nanoseconds to milliseconds.
RCUTILS_NS_TO_MS( nanoseconds ) = ( nanoseconds ÷ ( 1000 * 1000 ) ) # Convenience macro to convert nanoseconds to microseconds.
RCUTILS_NS_TO_US( nanoseconds ) = ( nanoseconds ÷ 1000 ) # A single point in time, measured in nanoseconds since the Unix epoch.
const RCUTILS_WARN_UNUSED = nothing
# const __locale_t = Cstring  # HACK: the autogen creates stubs for string manipulation functions that we don't actually want/nede, and it needs this to defined

# dynamically load the ROS2 libraries
const librmw = find_library(["librmw.$dlext"])
const librcl = find_library(["librcl.$dlext"])
const librcutils = find_library(["librcutils.$dlext"])

# const librmw = dlopen(librmw_path)
# const librcutils = dlopen(librcutils_path)
# const librcl = dlopen(librcl_path)

# include the auto-gen-ed wrappers for the C functions/structs/constants
include("gen_librcl_h.jl")
include("gen_librcl.jl")

struct RclError <: Exception
    code :: Cint
    msg :: AbstractString
end

function checkcall(code :: rcl_ret_t; suppress_throw=false)
    code == RCL_RET_OK && return nothing

    msg = rcutils_get_error_string()
    @assert msg != Cstring(C_NULL)
    err = RclError(code, unsafe_string(msg))

    rcutils_reset_error()

    suppress_throw || throw(err)
    err
end


text:
# a basic create robot type node example
addprocs(3)
using IncrementalInference, CloudGraphs

# switch IncrementalInference to use CloudGraphs (Neo4j) data layer
# Caesar.useCloudGraphsDataLayer()
    # connect to the server, CloudGraph stuff
    configuration = CloudGraphs.CloudGraphConfiguration("localhost", 7474, "", "", "localhost", 27017, false, "", "");

    cloudGraph = connect(configuration);

    # register types of interest in CloudGraphs
    CloudGraphs.registerPackedType!(cloudGraph, VariableNodeData, PackedVariableNodeData, encodingConverter=VNDencoder, decodingConverter=VNDdecoder);
    CloudGraphs.registerPackedType!(cloudGraph, FunctionNodeData{Obsv2}, FunctionNodeData{PackedObsv2}, encodingConverter=FNDencode, decodingConverter=FNDdecode)
    CloudGraphs.registerPackedType!(cloudGraph, FunctionNodeData{Odo}, FunctionNodeData{PackedOdo}, encodingConverter=FNDencode, decodingConverter=FNDdecode)
    CloudGraphs.registerPackedType!(cloudGraph, FunctionNodeData{GenericMarginal}, FunctionNodeData{GenericMarginal}, encodingConverter=FNDencode, decodingConverter=FNDdecode)
    CloudGraphs.registerPackedType!(cloudGraph, FunctionNodeData{Ranged}, FunctionNodeData{Ranged}, encodingConverter=FNDencode, decodingConverter=FNDdecode)

IncrementalInference.setCloudDataLayerAPI!()

# this is being replaced by cloudGraph, added here for development period
fg = emptyFactorGraph()
fg.cg = cloudGraph

# Robot navigation and inference type stuff
N=200
doors = [-100.0;0.0;100.0;300.0]'
cov = [3.0]

# robot style, add first pose vertex
v1 = addNode!(fg,"x1",doors,N=N)

# add a prior for the initial position of the robot
f0  = addFactor!(fg,[v1], Obsv2(doors, cov', [1.0]))

# add second pose vertex
tem = 2.0*randn(1,N)+getVal(v1)+50.0
v2 = addNode!(fg, "x2", tem, N=N)

# now add the odometry factor between them
f1 = addFactor!(fg,[v1;v2],Odo([50.0]',[2.0]',[1.0]))

v3=addNode!(fg,"x3",4.0*randn(1,N)+getVal(v2)+50.0, N=N)
addFactor!(fg,[v2;v3],Odo([50.0]',[4.0]',[1.0]))
f2 = addFactor!(fg,[v3], Obsv2(doors, cov', [1.0]))

v4=addNode!(fg,"x4",2.0*randn(1,N)+getVal(v3)+50.0, N=N)
addFactor!(fg,[v3;v4],Odo([50.0]',[2.0]',[1.0]))

if true
    l1=addNode!(fg, "l1", 0.5*randn(1,N)+getVal(v3)+64.0, N=N)
    addFactor!(fg, [v3,l1], Ranged([64.0],[0.5],[1.0]))
    addFactor!(fg, [v4,l1], Ranged([16.0],[0.5],[1.0]))
end

v5=addNode!(fg,"x5",2.0*randn(1,N)+getVal(v4)+50.0, N=N)
addFactor!(fg,[v4;v5],Odo([50.0]',[2.0]',[1.0]))

v6=addNode!(fg,"x6",1.25*randn(1,N)+getVal(v5)+40.0, N=N)
addFactor!(fg,[v5;v6],Odo([40.0]',[1.25]',[1.0]))

v7=addNode!(fg,"x7",2.0*randn(1,N)+getVal(v6) +60.0, N=N)
addFactor!(fg,[v6;v7],Odo([60.0]',[2.0]',[1.0]))

f3 = addFactor!(fg,[v7], Obsv2(doors, cov', [1.0]))

# Now operate with the data in the DB

tree = prepBatchTree!(fg,drawpdf=true)

# recursive solving (single process, easy stack trace for debugging)
inferOverTree!(fg, tree)

# get vertex back from DB
# x1neoID = fg.cgIDs[fg.IDs["x1"]]
# cv1r = CloudGraphs.get_vertex(fg.cg, x1neoID, false)

# Get neighbors
# neighs = CloudGraphs.get_neighbors(fg.cg, cv1r)

gt = Dict{ASCIIString, Array{Float64,2}}()
gt["x1"]=([0.0;1.97304 ]')' # -0.0342366
gt["x2"]=([50.0; 2.83153 ]')' # 49.8797
gt["x3"]=([100.0; 1.65557 ]')' # 99.8351
gt["x4"]=([150.0; 1.64945 ]')' # 148.637
gt["x5"]=([200.0; 1.77992 ]')' # 198.62
gt["x6"]=([240.0; 2.20466 ]')' # 238.492
gt["x7"]=([300.0; 2.14353 ]')' # 298.467
gt["l1"]=([165.0; 1.17284 ]')' # 164.102

# draw all beliefs
using Gadfly

DOYTICKS = false
xx,ll = ls(fg)
msgPlots = drawHorBeliefsList(fg, xx, gt=gt,nhor=2);
evalstr = ""
for i in 1:length(msgPlots)
    evalstr = string(evalstr, ",msgPlots[$(i)]")
end
pl = eval(parse(string("vstack(",evalstr[2:end],")")));
Gadfly.draw(PDF("cg4doors.pdf",15cm,20cm),pl) # can also do PNG


text:
using ITensors
using ITensorInfiniteMPS
using Random

# Parameters
cutoff = 1e-8
maxdim = 100
tol = 1e-8
maxiter = 20
outer_iters = 3

N = 2
model = Model"ising"()

function space_shifted(::Model"ising", q̃sz)
  return [QN("SzParity", 1 - q̃sz, 2) => 1, QN("SzParity", 0 - q̃sz, 2) => 1]
end

space_ = fill(space_shifted(model, 1), N)
s = infsiteinds("S=1/2", N; space=space_)
initstate(n) = "↑"
ψ = InfMPS(s, initstate)

model_params = (J=-1.0, h=0.9)
vumps_kwargs = (multisite_update_alg="sequential", tol=tol, maxiter=maxiter)
subspace_expansion_kwargs = (cutoff=cutoff, maxdim=maxdim)

H = InfiniteITensorSum(model, s; model_params...)
# Alternate steps of running VUMPS and increasing the bond dimension
ψ1 = vumps(H, ψ; vumps_kwargs...)
for _ in 1:outer_iters
  ψ1 = subspace_expansion(ψ1, H; subspace_expansion_kwargs...)
  ψ1 = vumps(H, ψ1; vumps_kwargs...)
end

Hmpo = InfiniteMPOMatrix(model, s; model_params...)
# Alternate steps of running VUMPS and increasing the bond dimension
ψ2 = vumps(Hmpo, ψ; vumps_kwargs...)
for _ in 1:outer_iters
  ψ2 = subspace_expansion(ψ2, Hmpo; subspace_expansion_kwargs...)
  ψ2 = vumps(Hmpo, ψ2; vumps_kwargs...)
end

SzSz = prod(op("Sz", s[1]) * op("Sz", s[2]))


text:
using Random
export rand_supremacy2d
export pair_supremacy, print_square_bond, cz_entangler
"""
control-Z entangler.
"""
cz_entangler(n::Int, pairs) = chain(n, control(n, [ctrl], target=>Z) for (ctrl, target) in pairs)

"""
    rand_supremacy2d(nx::Int, ny::Int, depth::Int; seed=1) -> AbstactBlock

random supremacy circuit.

NOTE: the restriction to `T` gate is removed.
"""
function rand_supremacy2d(nx::Int, ny::Int, depth::Int; seed=1)
    Random.seed!(seed)
    nbits = nx*ny
    entanglers = map(pair->cz_entangler(nbits, pair), pair_supremacy(nx, ny))
    gateset = [ConstGate.T, SqrtX, SqrtY]
    c = chain(nbits, repeat(nbits, H, 1:nbits))
    pre_occ = (1:nbits...,)
    pre_gates = Vector{PrimitiveBlock}(fill(H,nbits))
    hastgate = zeros(Bool, nbits)
    for i=1:depth-2
        ent = entanglers[mod1(i,length(entanglers))]
        unit = chain(nbits, [ent])
        occ = occupied_locs(ent)
        for loc in setdiff(pre_occ, occ)
            g1 = hastgate[loc] ? rand(setdiff(gateset, (pre_gates[loc],))) : (hastgate[loc]=true; T)
            push!(unit, put(nbits, loc=>g1))
            pre_gates[loc] = g1
        end
        pre_occ = occ
        push!(c, unit)
    end
    depth!=1 && push!(c, repeat(nbits, H, 1:nbits))
    return c
end

"""obtain supremacy pairing patterns"""
function pair_supremacy(nx::Int, ny::Int; periodic=false)
    Kx = [0.25 0.5; 0.25 -0.5]
    Ky = [0.5 0.25; -0.5 0.25]
    out = Vector[]
    li = LinearIndices((nx, ny))

    for (dx, dy, isxdirection) in [(2,0,true), (0,0,true), (0,3,false), (0,1,false),
                                    (3,0,true), (1,0,true), (0,0,false), (0,2,false)]
        res = Pair{Int, Int}[]
        for i = 1:nx, j=1:ny
            if (periodic || (isxdirection ? i<nx : j<ny)) && all(mod.((isxdirection ? Kx : Ky)*[i-dx-1, j-dy-1], 1) .≈ 0)
                i_ = isxdirection ? mod1(i+1, nx) : i
                j_ = isxdirection ? j : mod1(j+1, ny)
                push!(res, li[i, j] => li[i_, j_])
            end
        end
        push!(out, res)
    end

    return out
end

print_square_bond(nx::Int, ny::Int, bonds) = print_square_bond(stdout, nx, ny, bonds)
function print_square_bond(io::IO, nx::Int, ny::Int, bonds)
    li = LinearIndices((nx, ny))
    for j=1:ny
        for i=1:nx
            print(io, "∘")
            i!=nx && print(io, ((li[i,j] => li[i+1,j]) in bonds || (li[i+1,j] => li[i,j]) in bonds) ? "---" : "   ")
        end
        println(io)
        j!=ny && for i=1:nx
                print(io, ((li[i,j] => li[i,j+1]) in bonds || (li[i,j] => li[i,j+1]) in bonds) ? "|" : " ")
                j!=ny && print(io, "   ")
        end
        println(io)
    end
end

max_stars_repo_name: QuantumBFS/YaoExtensions

max_issues_repo_name: QuantumBFS/YaoExtensions

max_forks_repo_name: QuantumBFS/YaoExtensions


text:
@testset "generate" begin
    root = "assets"

    sec = DemoSection(joinpath(root, "section", "default"))
    @test_reference joinpath("references", "generate_section.txt") generate(sec)

    page = DemoPage(joinpath(root, "page", "default"))
    @test_reference joinpath("references", "generate_page.txt") generate(page)

    abs_root = joinpath(pwd(), root, "page", "default")
    mktempdir() do root
        @test @suppress_err makedemos(abs_root, root=root) == joinpath("democards", "default", "index.md")
        @test cardtheme(root=root) == joinpath("democards", "cardtheme.css")
    end
end


text:
# specializations for multi-output GPs
# names: MultiOuputGaussian, MultiOutputProcess, VectorGaussian, VectorGaussianProcess
struct VectorGaussian{M, S} <: AbstractGaussian{M, S}
    μ::M
    Σ::S
	d::Int # output dimension
end
VectorGaussian(k, d::Int) = VectorGaussian(zero_mean(d), k, d)

(G::VectorGaussian)(x) = Gaussian(mean(G)(x), cov(G)(x, x))

Statistics.var(G::VectorGaussian) = x->max.(diag(G.Σ(x, x)), 0)
Statistics.std(G::VectorGaussian) = x->sqrt.(var(G)(x))
Statistics.mean(G::VectorGaussian, x::AbstractVector) = marginal(mean(G), x, G.d)
# TODO: call var(G).(x) for efficiency?
# Statistics.var(G::VectorGaussian, x::AbstractVector) = max.(diag(gramian(cov(G), x)), 0)
# # function var(k::MultiKernel, x::AbstractVector{<:Real}) end
# Statistics.std(G::VectorGaussian, x::AbstractVector) = sqrt.(var(G, x))

# getindex for multi-output GP selects specific ouput dimenion
function Base.getindex(G::VectorGaussian, i)
	μ(x) = mean(G)(x)[i]
	Σ(x, y) = Σ isa MultiKernel ? cov(G)[i, i](x, y) : cov(G)(x, y)[i, i]
	Gaussian(μ, Σ)
end

############################ vector conditional ################################
# iterative_trait(::AbstractMatOrFac) = Val(false)
# iterative_trait(::BlockFactorization) = Val(true)
function conditional(G::VectorGaussian, x::AbstractVector,
                     Gy::Union{AbstractVector{<:Number}, Multivariate}; tol::Real = tol)
    length(x) > 0 || return G # return if an empty vector is passed
    K = gramian(cov(G), x)
    K = length(x) == 1 ? Matrix(K) : factorize(K, tol = tol)
    y = Gy
    if Gy isa Gaussian
        y = mean(Gy)
        K = K isa Matrix ? K + cov(Gy) : LazyMatrixSum(K, cov(Gy)) # add noise variance to kernel matrix
    end
	d = length(Gy) ÷ length(x)
    μ = VectorConditionalMean(G, x, y, d, K)
    Σ = VectorConditionalKernel(G, x, K)
    return VectorGaussian(μ, Σ, d)
end

# conditional helpers
function conditional(G::VectorGaussian, x, y::Union{AbstractMatrix, VecOfVec}; tol::Real = tol)
	conditional(G, x, _vec(y); tol = tol)
end

# if Y is passed as a matrix, Σ is interpreted as a matricified version of the corresponding diagonal covariance
function conditional(G::VectorGaussian, x, y::MT, Σ::MT; tol::Real = tol) where {MT <: Union{AbstractMatrix, VecOfVec}}
	Gy = Gaussian(_vec(y), Diagonal(_vec(Σ)))
	conditional(G, x, Gy, tol = tol)
end

function conditional(G::VectorGaussian, x, y::Union{AbstractMatrix, AbstractVector}, σ²::Real; tol::Real = tol)
	Gy = Gaussian(_vec(y), σ²)
	conditional(G, x, Gy; tol = tol)
end

function conditional(G::VectorGaussian, x, y::AbstractVector{<:Number}, Σ; tol = tol)
	conditional(G, x, Gaussian(y, Σ), tol = tol)
end

# conditional with matrix input X
# unnecessary if we make first argument AbstractGaussian in main conditional file
function conditional(G::VectorGaussian, X::AbstractMatrix, y::Union{AbstractVector, Multivariate}; tol = tol)
	conditional(G, [x for x in eachcol(X)], y, tol = tol)
end

############################### Conditional Mean ###############################
# here, T would be the output type ..., or should it be field type?
struct VectorConditionalMean{M, K, U<:AbstractVector, V<:AbstractVector}
	μ::M # prior mean of predictions
	kxu::K # covariance between prediction and observation points (not necessarily prior cov for multi-GP)
    u::U
    α::V # α = K \ (y - μ(x)) IDEA: sparsity?
	d::Int
end
# converting Gaussian to separate mean and covariance function
VectorConditionalMean(G::VectorGaussian, x...) = VectorConditionalMean(G.μ, G.Σ, G.μ, x...)
# μx is mean prediction
# μu is mean observation
# kxu is covariance between predictions and observations
# y are observation at locations u
# Kuu is the factorized covariance matrix between observations
function VectorConditionalMean(μx, kxu, μu, u::AbstractVector, y::AbstractVector{<:Real}, d::Int, Kuu)
    μ = marginal(μu, u, d)
    α = Kuu \ difference(y, μ)
    VectorConditionalMean(μx, kxu, u, α, d)
end
(M::VectorConditionalMean)(x) = M.μ(x) + gramian(M.kxu, [x], M.u) * M.α

function Base.getindex(M::VectorConditionalMean, i::Int)
	kxu(x, y) = Matrix(M.kxu(x, y))[i, :] # covariance of ith predictive with all observation dimensions
	# kxu(x, y) = M.kxu(x, y)[i, :] # possible without Matrix?
	mean_component(x) = M.μ(x)[i, :] + gramian(kxu, [x], M.u) * M.α
end
# IDEA, but change-intensive to implement
# function Base.getindex(M::VectorConditionalMean, i::Int)
# 	k(x, y) = M.kxu(x, y)[i, :] # covariance of ith predictive with all observation dimensions
# 	ConditionalMean(M.μ, k, M.x, M.α)
# end
function marginal(M::VectorConditionalMean, x::AbstractVector)
    μ = marginal(M.μ, x, M.d)
    K = gramian(M.kxu, x, M.u)
    return mul!(μ, K, M.α, 1, 1)
end

########################## Conditional Kernel ##################################
# corresponds to the kernel of the conditional GP, with observations at x
# ASSUMPTION: x is vector of data, even if it is just a single element
struct VectorConditionalKernel{T, K, KX, U<:AbstractVector,
                    		   S<:AbstractMatOrFac{T}} <: MultiKernel{T}
    kxx::K # prior covariance function of predictions
	kxu::KX # covariance function between observations and predictions
	# kuy::KY # TODO: covariance function between observations and predictions
    u::U # input data on which we are conditioning
    Kuu::S # covariance between x and x
end
# default to having same covariance between observations and predictions
function VectorConditionalKernel(k, u::AbstractVector, Kuu::AbstractMatOrFac)
	VectorConditionalKernel(k, k, u, Kuu)
end

# Kxx might be constructed based on kernel type
function VectorConditionalKernel(k::MultiKernel, u::AbstractVector)
    Kuu = factorize(gramian(k, u))
    VectorConditionalKernel(k, u, Kuu)
end
function VectorConditionalKernel(G::VectorGaussian, u::AbstractVector, Kuu::AbstractMatOrFac)
    VectorConditionalKernel(cov(G), u, Kuu)
end

function (C::VectorConditionalKernel)(x::AbstractVector, y::AbstractVector)
    Kxu = gramian(C.kxu, [x], C.u) # why Matrix?
    Kyu = (x ≡ y) ? Kxu : gramian(C.kxu, [y], C.u)
	K = Matrix(C.kxx(x, y)) # need to convert this to Matrix to allocate, or do Woodbury
	K .-= Kxu * (C.Kuu \ Matrix(Kyu)') # (Kux', inverse(C.Kxx), Kuy)
	# Woodbury(C.kxx(x, y), Kxu, inverse(C.Kuu), Kyu') # inverse should already work
end

Base.getindex(C::VectorConditionalKernel, i, j) = (x, y)->C(x, y, i, j)
# only calculates output variance of single component
function (C::VectorConditionalKernel)(x::AbstractVector, y::AbstractVector, i, j)
	Kux = Matrix(gramian(C.kxu, [x], C.u))'[:, i] # can we do this without matrix?
	Kuy = (x ≡ y) && (i == j) ? Kux : Matrix(gramian(C.kxu, [y], C.u))'[:, j]
	K = C.kxx(x, y)[i, j]
	M = Kux' * (C.Kuu \ Kuy)
	K .- M # *(Kux', inverse(C.Kxx), Kuy)
end

# for efficiency, gramian should dispatch on the type of the factorization of c
# does it make sense to return a lazy matrix? given each element requires O(n^2) cost ...
import CovarianceFunctions: gramian
function gramian(C::VectorConditionalKernel, x::AbstractVector, y::AbstractVector)
   Kxu = gramian(C.kxu, x, C.u)
   Kyu = (x ≡ y) ? Kxu : gramian(C.kxu, y, C.u)
   Woodbury(gramian(C.kxx, x, y), Kux', inverse(C.Kxx), Kuy, -1)
end

# IDEA: take inspiration from finite-dimensional linear operator conditonal in conditiona.jl
# # # global optimization via deconvolution
# function conditional(G::Gaussian, H::VectorGaussian, X::AbstractMatrix, y::AbstractVector)
# 	k, h = cov(G), cov(H)
#     s = k + h # sum of kernels corresponds to convolution of GPs
#     G = gramian(s, X)
#     F = factorize(G)
#     α = F \ y
#     x = [c for c in eachcol(X)]
# 	VectorConditionalMean()
# 	VectorConditionalKernel()
# end

# doodles
################################################################################
# constructing multi-output GP with separable covariance kernel
# function Base.:*(A::Union{Real, AbstractMatOrFacOrUni}, G::Gaussian)
#     μ = x->A*mean(G, x)
#     Σ = Kernel.Separable(A, k)
#     Gaussian(μ, Σ)
# end
################################################################################

# multi-output process with nd input and nd ouput
# this doesn't work for conditional kernel
# function conditional(G::Gaussian{<:Any, <:MultiKernel}, x::AbstractVector{<:Real}, y::AbstractVector{<:Real},
#     Σ = nothing, ::Val{true} = ismulti(cov(G)); tol::Real = 1e-12)
#     conditional(G, [x], [y], Σ, tol = tol)
# end

# if we don't rely on Kernel.jl, this is still a bit of a headache
# function conditional(G::Gaussian, x::AbstractVector{<:Real}, y::AbstractVector{<:Real},
#     Σ = nothing, ::Val{true} = ismulti(cov(G)); tol::Real = 1e-12)
#     conditional(G, [x], [y], Σ, tol = tol)
# end
# function ismulti end
# ismulti(::AbstractKernel) = Val(false)
# ismulti(::MultiKernel) = Val(true)
# ismulti(k::Any) = error("to enable single-datapoint conditioning with a multi-ouput process, please define ismulti(::$(typeof(k)) = Val(true)")
# ismulti(::ConditionalKernel{K}) where K = ismulti(K)


text:
using DelimitedFiles
using JuMP
using Ipopt
using LinearAlgebra

debug = false

println("Definiendo el problema...")

model = Model(Ipopt.Optimizer)

@variable(model, alpha[1:20] >= 0)

@objective(model, Min, -sum(alpha))

println("Variable y FO definidas...")

m = readdlm("./features/txt/features.txt")
z = 1

for col in eachcol(m)
    global z

    constr1 = sum((alpha[i]*col[i])^2 for i = 1:20)

    if debug
        println(constr)
    end

    @NLconstraint(model, constr1 <= z^2)

    println("Restriccion ", z, " implementada...")
    z = z+1
end

println("Optimizando!...")
optimize!(model)

# Resultados
println("Resultados!...")

println("alpha1: ", JuMP.value(alpha[1]))
println("alpha2: ", JuMP.value(alpha[2]))
println("alpha3: ", JuMP.value(alpha[3]))
println("alpha4: ", JuMP.value(alpha[4]))
println("alpha5: ", JuMP.value(alpha[5]))
println("alpha6: ", JuMP.value(alpha[6]))
println("alpha7: ", JuMP.value(alpha[7]))
println("alpha8: ", JuMP.value(alpha[8]))
println("alpha9: ", JuMP.value(alpha[9]))
println("alpha10: ", JuMP.value(alpha[10]))
println("alpha11: ", JuMP.value(alpha[11]))
println("alpha12: ", JuMP.value(alpha[12]))
println("alpha13: ", JuMP.value(alpha[13]))
println("alpha14: ", JuMP.value(alpha[14]))
println("alpha15: ", JuMP.value(alpha[15]))
println("alpha16: ", JuMP.value(alpha[16]))
println("alpha17: ", JuMP.value(alpha[17]))

open("outputs/optimization/spotify_estimator_alpha_actual.txt", "w") do io
    write(io, string([JuMP.value(alpha[i]) for i = 1:20]))
end;

max_stars_repo_name: RodrigoNazar/SuccessPredictor

max_issues_repo_name: RodrigoNazar/SuccessPredictor

max_forks_repo_name: RodrigoNazar/SuccessPredictor


text:
function eta=EtaFromZ(x,y,z,Param)
etaOld=1.e-7;
F=-Param.Grav*z+PhiBaroWave(x,y,etaOld,Param);
FStart=F;
for i=1:100
  %F=-Param.Grav*z+PhiBaroWave(x,y,etaOld,Param);
  FPrim=-Param.Rd/etaOld*TBaroWave(x,y,etaOld,Param);
  eta=etaOld-F/FPrim;
  if abs(eta-etaOld)<=1.e-13
    break
  end
  etaOld=eta;
  F=-Param.Grav*z+PhiBaroWave(x,y,etaOld,Param);
%   if abs(F)<=abs(FStart)*1.e-12
%     break
%   end
end
end


text:
struct NullCache <:AbstractCacheStrategy end
function make_cache(T, ::NullCache, ax_list)
    axes = reduce(ax_list; init=()) do x, y
        (x..., y...)
    end
    NullCacheArray{T, length(axes), typeof(axes)}(axes)
end

# cache implementation details
struct NullCacheArray{T, N, AT} <: CachedView{T, N}
    axes::AT
end
Base.axes(C::NullCacheArray) = C.axes
Base.size(C::NullCacheArray) = map(length, axes(C))


text:
# Test that https://www.google.com successfully connects
curl = curl_easy_init()
curl == C_NULL && error("curl_easy_init() failed")

# Setup the callback function to recv data
function curl_write_cb(curlbuf::Ptr{Cvoid}, s::Csize_t, n::Csize_t, p_ctxt::Ptr{Cvoid})
    sz = s * n
    data = Array{UInt8}(undef, sz)

    ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, UInt64), data, curlbuf, sz)

    sz::Csize_t
end

@testset "SSL verify" begin
    # Set up the write function to consume the curl output so we don't see it in the
    # test output
    c_curl_write_cb = @cfunction(
        curl_write_cb,
        Csize_t,
        (Ptr{Cvoid}, Csize_t, Csize_t, Ptr{Cvoid})
    )
    curl_easy_setopt(curl, CURLOPT_WRITEFUNCTION, c_curl_write_cb)

    # Set up our SSL options
    curl_easy_setopt(curl, CURLOPT_URL, "https://www.google.com")
    curl_easy_setopt(curl, CURLOPT_USE_SSL, CURLUSESSL_ALL)
    curl_easy_setopt(curl, CURLOPT_SSL_VERIFYHOST, 2)
    curl_easy_setopt(curl, CURLOPT_SSL_VERIFYPEER, 1)

    @testset "SSL Success" begin
        res = curl_easy_perform(curl)
        # NOTE: This will fail on macOS due to the cross compilation not knowing where the
        # macOS CACerts would be
        # The same issue would exist with Windows and FreeBSD as well
        if Sys.isbsd() || Sys.iswindows()
            @test_broken res == CURLE_OK
        else
            @test res == CURLE_OK
        end
    end

end

curl_easy_cleanup(curl)


text:
############################################################################################
# 				     	CONTROL BOUND CONSTRAINTS 										   #
############################################################################################
"""
	ControlBoundConstraint{P,M,T}

Linear bound constraint on controls
# Constructors
```julia
ControlBoundConstraint(m; u_min, u_max)
```
Any of the bounds can be ±∞. The bound can also be specifed as a single scalar, which applies the bound to all controls.
"""
struct ControlBoundConstraint{P,M,T} <: TrajectoryOptimization.ControlConstraint
	m::Int
	u_max::SVector{M,T}
	u_min::SVector{M,T}
	i_max::Vector{Int}
	i_min::Vector{Int}
	inds::SVector{P,Int}
end

Base.copy(bnd::ControlBoundConstraint{P,m,T}) where {P,m,T} =
	ControlBoundConstraint(bnd.m, bnd.u_max, bnd.u_min,
		copy(bnd.i_max), copy(bnd.i_min), bnd.inds)

function ControlBoundConstraint(m; u_max=Inf*(@SVector ones(m)), u_min=-Inf*(@SVector ones(m)))

	# Check and convert bounds
	u_max, u_min = checkBounds(m, u_max, u_min)

	# Concatenate bounds
	b = [-u_max; u_min]
	inds = findall(isfinite, b)
	inds = SVector{length(inds)}(inds)

	# Get linear indices of 1s of Jacobian
	a_max = findall(isfinite, u_max)
	a_min = findall(isfinite, u_min)
	u = length(a_max)
	l = length(a_min)
	carts_u = [CartesianIndex(i,   j) for (i,j) in enumerate(a_max)]
	carts_l = [CartesianIndex(i+u, j) for (i,j) in enumerate(a_min)]
	∇c = zeros(u+l, m)
	linds_u = LinearIndices(zeros(u+l,m))[carts_u]
	linds_l = LinearIndices(zeros(u+l,m))[carts_l]

	ControlBoundConstraint(m, u_max, u_min, linds_u, linds_l, inds)
end

function TrajectoryOptimization.con_label(con::ControlBoundConstraint, ind::Int)
	i = con.inds[ind]
	m = control_dim(con)
	if 1 <= i <= m
		return "u max $i"
	elseif m < i <= 2m
		j = i - m
		return "u min $j"
	else
		throw(BoundsError())
	end
end

function checkBounds(n::Int, u::AbstractVector, l::AbstractVector)
	if all(u .>= l)
		return SVector{n}(u), SVector{n}(l)
	else
		throw(ArgumentError("Upper bounds must be greater than or equal to lower bounds"))
	end
end

checkBounds(n::Int, u::Real, l::Real) =
	checkBounds(n, (@SVector fill(u,n)), (@SVector fill(l,n)))
checkBounds(n::Int, u::AbstractVector, l::Real) = checkBounds(n, u, fill(l,n))
checkBounds(n::Int, u::Real, l::AbstractVector) = checkBounds(n, fill(u,n), l)

@inline TrajectoryOptimization.control_dim(con::ControlBoundConstraint) = con.m
@inline TrajectoryOptimization.is_bound(::ControlBoundConstraint) = true
@inline TrajectoryOptimization.lower_bound(bnd::ControlBoundConstraint) = bnd.u_min
@inline TrajectoryOptimization.upper_bound(bnd::ControlBoundConstraint) = bnd.u_max
@inline TrajectoryOptimization.sense(::ControlBoundConstraint) = Inequality()
@inline Base.length(con::ControlBoundConstraint) = length(con.i_max) + length(con.i_min)

function TrajectoryOptimization.primal_bounds!(zL, zU, bnd::ControlBoundConstraint)
	for i = 1:length(zL)
		zL[i] = max(bnd.u_min[i], zL[i])
		zU[i] = min(bnd.u_max[i], zU[i])
	end
	return true
end

function TrajectoryOptimization.evaluate(bnd::ControlBoundConstraint, z::AbstractKnotPoint)
	[(control(z) - bnd.u_max); (bnd.u_min - control(z))][bnd.inds]
end

function TrajectoryOptimization.jacobian!(∇c, bnd::ControlBoundConstraint{U,L}, z::AbstractKnotPoint) where {U,L}
	for i in bnd.i_max
		∇c[i]  = 1
	end
	for i in bnd.i_min
		∇c[i] = -1
	end
	return true
end

TrajectoryOptimization.∇jacobian!(G, con::ControlBoundConstraint, z::AbstractKnotPoint, λ::AbstractVector) = true # zeros

function TrajectoryOptimization.change_dimension(con::ControlBoundConstraint, m::Int, iu=1:m)
	m0 = con.m
	u_max = fill(Inf,m)
	u_min = fill(Inf,m)
	u_max[iu] = con.u_max[1:m0]
	u_min[iu] = con.u_min[1:m0]
	ControlBoundConstraint(m, u_max=u_max, u_min=u_min)
end


text:
#__precompile__()
module SparseUtils

import SparseArrays
import SparseArrays: SparseMatrixCSC, nnz
import LinearAlgebra
#import DataUtils
import StatsBase
import Printf

const Callable = Union{Function, Type}

export
    c_to_julia_index,
    c_to_julia_index!,
    density,
    nnzcounts,
    summarystats,
    nrows,
    ncols,
    prunecols!,
    prunecols,
    prunerows,
    renumbercols,
    renumberrows,
    renumberrowscols,
    hasemptycols,
    hasemptyrows

include("findrepeated.jl")
include("nzval.jl")
include("ijv.jl")
include("csc.jl")
include("coo.jl")
include("coorenumber.jl")
include("abstract.jl")
include("countmap.jl")

import .NZVal, .IJV

end # module SparseUtils

#  LocalWords:  nnzcounts countmap sp


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct Timestamp <: ProtoType
    seconds::Int64
    nanos::Int32
    Timestamp(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Timestamp

export Timestamp


text:
module Lib

export ackermann

function ackermann(m, n)
  if m == 0
    return n + 1
  elseif n == 0
    return ackermann(m-1, 1)
  else
    return ackermann(m-1, ackermann(m, n-1))
  end
end

end

max_stars_repo_name: ashwinvis/experiments-with-haskell

max_issues_repo_name: ashwinvis/experiments-with-haskell

max_forks_repo_name: ashwinvis/experiments-with-haskell